text:
#!/usr/bin/env python
"""
Example script that can compare diagnostics from two different netCDF
collections.  Similar to compute_diagnostics.ipynb, but can be used
without having to open a Jupyter notebook.
"""

# Imports
import os
import sys
import warnings
import numpy as np
import xarray as xr
from yaml import load as yaml_load_file
import gcpy.benchmark as bmk
import gcpy.constants as constants
import gcpy.util as util

# Tell matplotlib not to look for an X-window
os.environ["QT_QPA_PLATFORM"] = "offscreen"

# Suppress harmless run-time warnings (mostly about underflow in division)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)

# ======================================================================
# Methods
# ======================================================================

def create_dirs(config):
    """
    Create directories for plots and weights if they do not exist.

    Arguments:
        config : dict
            Configuration information read from a YAML file
    """

    # Extract fields from the config object
    do_level_plot = config["options"]["level_plot"]["create_plot"]
    do_zonal_mean = config["options"]["zonal_mean"]["create_plot"]
    plotsdir = config["paths"]["plots_dir"]
    weightsdir = config["paths"]["weights_dir"]

    # Create plotsdir and weightsdir if necessary
    if do_level_plot or do_zonal_mean:
        if not os.path.isdir(plotsdir):
            os.mkdir(plotsdir)
            if not os.path.isdir(weightsdir):
                os.mkdir(weightsdir)

def read_data(config):
    """
    Read data from the Ref and Dev datasets

    Arguments:
        config : dict
            Configuration information read from a YAML file

    Returns:
        data : dict
            Contains Ref and Dev data as xarray Dataset fields.
    """

    # If we are using the gcpy_test data, use the gcpy_test_dir
    if config["options"]["gcpy_test"]:
        rootdir = config["paths"]["test_data_dir"]
    else:
        rootdir = config["paths"]["main_dir"]

    # Define paths to Ref & Dev files
    ref_file = os.path.join(
        rootdir,
        config["data"]["ref"]["dir"],
        config["data"]["ref"]["subdir"],
        config["data"]["ref"]["file"]
    )
    dev_file = os.path.join(
        rootdir,
        config["data"]["dev"]["dir"],
        config["data"]["dev"]["subdir"],
        config["data"]["dev"]["file"]
    )

    # Read Ref data
    try:
        refdata = xr.open_dataset(
            ref_file,
            drop_variables=constants.skip_these_vars
        )
    except Exception:
        msg = "Error reading " + ref_file
        raise Exception(msg)

    # Read Dev data
    try:
        devdata = xr.open_dataset(
            dev_file,
            drop_variables=constants.skip_these_vars
        )
    except Exception:
        msg = "Error reading " + dev_file
        raise Exception(msg)

    # Define dictionary for return
    data = {
        "ref": refdata,
        "dev": devdata
    }

    return data

def print_totals_and_diffs(config, refdata, devdata, varlist):
    """
    Prints totals and differences of data.  This is a quick way
    to see if there are nonzero differences between Ref and Dev
    datasets.

    Arguments:
        config : dict
            Configuration information read from a YAML file
        refdata : xarray Dataset
            Contains data from the Ref model run
        devdata : xarray Dataset
            Contains data from the Dev model run
        varlist : list of str
            Contains a list of data variables.
    """

    # Get quantities from the config object
    filename = config["options"]["totals_and_diffs"]["filename"]
    do_screen = config["options"]["totals_and_diffs"]["print_to_screen"]

    # Determine if we will print to a file
    do_file = len(filename) > 0
    if not do_file and not do_screen:
        print("...Printing to screen and file are both disabled!")
        return

    # Open a text file for output if a filename is specified
    if do_file:
        pathname = os.path.join(
            config["paths"]["plots_dir"],
            filename
        )
        f = open(pathname, 'w')

    # If we are printing only variables with zero diffs,
    # then alert the user (print to screen & file)
    if config["options"]["totals_and_diffs"]["skip_zero_diffs"]:
        line = "... Only showing variables with nonzero differences"
        if do_file:
            print(line, file=f)
        else:
            print(line)

    # Print top header
    line = "{} Ref={} Dev={} {}".format(
        "Variable".ljust(22),
        config["data"]["ref"]["label"].ljust(20),
        config["data"]["dev"]["label"].ljust(20),
        "Dev-Ref"
    )
    if do_file:
        print(line, file=f)
    if do_screen:
        print(line)

    # Always print nonzero differences, but only print zero differences
    # if the configuration option "skip_zero_diffs" is False.
    for v in varlist:
        refsum = np.sum(refdata[v].values)
        devsum = np.sum(devdata[v].values)
        diff = devsum - refsum

        if np.abs(diff) > 0.0:
            line = "{} : {} | {} | {} ".format(
                v.ljust(20),
                str(refsum).ljust(22),
                str(devsum).ljust(22),
                diff
            )
            if do_file:
                print(line, file=f)
            if do_screen:
                print(line)
        else:
            if not config["options"]["totals_and_diffs"]["skip_zero_diffs"]:
                line = "{} : {} | {} | {} ".format(
                    v.ljust(20),
                    str(refsum).ljust(22),
                    str(devsum).ljust(22),
                    diff
                )
                if do_file:
                    print(line, file=f)
                if do_screen:
                    print(line)

    # Close file
    if do_file:
        f.close()

def compare_data(config, data):
    """
    Compares data from two different xarray datasets.

    Args:
        data : dict
            Contains Ref and Dev data as xarray Dataset fields.
    """

    # Get xarray datasets from the data object
    refdata = data["ref"]
    devdata = data["dev"]

    # For each variable in refdata, but not non devdata, add an
    # array of NaN values to refdata. Ditto for devdata.  This will
    # allow us to show that the variable is missing in either
    # refdata or devdata.
    [refdata, devdata] = util.add_missing_variables(refdata, devdata)

    # Get the list of common variable names
    verbose = config["options"]["verbose"]
    quiet = not verbose
    vardict = util.compare_varnames(refdata, devdata, quiet=quiet)
    varlist_level = vardict["commonvars2D"] + vardict["commonvars3D"]
    varlist_zonal = vardict["commonvars3D"]

    # Restrict variables to those containing a given substring
    restrict_vars = config["options"]["restrict_vars"]
    if len(restrict_vars) > 0:
        varlist_level = [v for v in varlist_level if v in restrict_vars]
        varlist_zonal = [v for v in varlist_zonal if v in restrict_vars]

    # ==================================================================
    # Generate the single level comparison plot
    # ==================================================================
    if config["options"]["level_plot"]["create_plot"]:
        print('... Creating single-level plots')
        pdfname = os.path.join(
            config["paths"]["plots_dir"],
            config["options"]["level_plot"]["pdfname"]
        )
        bmk.compare_single_level(
            refdata,
            config["data"]["ref"]["label"],
            devdata,
            config["data"]["dev"]["label"],
            ilev=config["options"]["level_plot"]["level_to_plot"],
            varlist=varlist_level,
            pdfname=pdfname,
            weightsdir=config["paths"]["weights_dir"],
            verbose=verbose
        )

    # ==================================================================
    # Generate the zonal mean comparison plot
    # ==================================================================
    if config["options"]["zonal_mean"]["create_plot"]:
        print('... Creating zonal mean plots')
        pdfname = os.path.join(
            config["paths"]["plots_dir"],
            config["options"]["zonal_mean"]["pdfname"]
        )
        bmk.compare_zonal_mean(
            refdata,
            config["data"]["ref"]["label"],
            devdata,
            config["data"]["dev"]["label"],
            varlist=varlist_zonal,
            pdfname=pdfname,
            weightsdir=config["paths"]["weights_dir"],
            verbose=verbose
        )

    # ==================================================================
    # Print totals for each quantity
    # ==================================================================
    if config["options"]["totals_and_diffs"]["create_table"] or \
       config["options"]["totals_and_diffs"]["print_to_screen"]:
        print('... Printing totals and differences')
        print_totals_and_diffs(
            config,
            refdata,
            devdata,
            varlist_level
        )

def main():
    """
    Main program, reads data and calls compare_data to make plots.
    """

    # Take the config file as the 2nd argument (or use a default)
    # NOTE: sys.argv[0] is always the program name!
    config_filename = sys.argv[1] if len(sys.argv) == 2 else "compare_diags.yml"

    # Get paths and options from the configuration file
    config = util.read_config_file(config_filename)

    # Create dirs for plots & weights (if necessary)
    create_dirs(config)

    # Read and compare data
    compare_data(config, read_data(config))

if __name__ == "__main__":
    main()

max_stars_repo_name: NanaWei-Nevaeh/gcpy

max_stars_repo_licenses: NCSA Apache-2.0 MIT

max_issues_repo_name: NanaWei-Nevaeh/gcpy

max_issues_repo_licenses: NCSA Apache-2.0 MIT

max_forks_repo_name: NanaWei-Nevaeh/gcpy

max_forks_repo_licenses: NCSA Apache-2.0 MIT


text:
import time
import sys
import pandas as pd
import os
import shutil
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from IPython.display import clear_output
from IPython.display import display
import category_encoders as ce
from itertools import combinations
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.preprocessing import PolynomialFeatures

recycle = pd.DataFrame()
checkpoint_ = {}

class DataObjects:
    """Initializing the Dataobjects
    arg:
    train: train DataFrame
    test:  test DataFrame
    submission: Sample Submission Data Frame
    data: Merged train and test
    """

    def __init__(self):
        self.train = None
        self.test = None
        self.submission = None
        self.data = None
        self.encoded_data = None
        self.categories = None
        self.bundle = None
        self.checkpoint_ = {}

    def merge(self, info=True):
        """
        merge two datasets for evaluation
        """
        self.train['source'] = 'train'
        self.test['source'] = 'test'
        self.data = pd.concat([self.train, self.test], ignore_index=True)
        if info:
            print(f'Test Shape: {self.test.shape},\nTrain Shape: {self.train.shape}, \nData Shape: {self.data.shape}')

    def split(self):
        """
        split datasets into train and test
        """
        self.train = self.data[self.data.source == "train"]
        self.test = self.data[self.data.source == "test"]
        print(self.train.shape, self.test.shape)
        del self.train['source']
        del self.test['source']
        return self.train, self.test

    def checkpoint(self, save=False, index=0):
        """
        creates a checkpoint with the option of saving to the disk
        """
        try:
            os.mkdir('checkpoints')
        except:
            print('No New Directory was made..')
        self.checkpoint_.update({(f'checkpoint{len(self.checkpoint_)}'): self.data})
        print('New Checkpoint Created')
        if save:
            self.data.to_csv(f"checkpoints/Checkpoint{len(os.listdir('checkpoints'))}.csv", index=None)

    def checkpoints(self, restore=None, load=False):
        """ restores a previous saved checkpoint!"""
        if load:
            file_ = []
            for csv in os.listdir('checkpoints'):
                file_.append(pd.read_csv('checkpoints/' + csv))
            for csv in file_:
                self.checkpoint_.update({('checkpoint{}f'.format(len(self.checkpoint_), )): csv})
        if restore is None:
            print(self.checkpoint_.keys())
        if restore is not None:
            self.data = self.checkpoint_[restore]

    def encodings(self, one_hot=False, catboost=True):
        self.encoded_data = self.data
        cat_feats = [column for column in self.encoded_data.select_dtypes(include='object').columns]
        if (one_hot == False) & (catboost == True):
            extra = [column for column in self.encoded_data.select_dtypes(
                include='int64').columns if len(self.encoded_data[column].unique()) < 10]
            cat_feats.extend(extra)
            self.categories = cat_feats
            self.categories.remove('source')
        if (one_hot == False) & (catboost == False):
            for column in cat_feats:
                self.catcode(column, encoded=True)
        if (one_hot == True) & (catboost == False):
            for column in cat_feats:
                self.encoded_data = pd.concat(
                    [self.encoded_data, pd.get_dummies(self.encoded_data[column], prefix=column)], 1)
                self.encoded_data.drop(column, axis=1, inplace=True)

    def bundles(self, verbose=False):
        """
        this creates a bundle of the train and test set, almost similar to a catboost pool method
        """
        if verbose:
            print('Bundling..')
        self.bundle = [self.train, self.test]

    def categoricals(self, categories, bundle, target):
        """
        creates a bunch of categorical encodings using the method learnt from the kaggle feature
        engineering class- Count encodings are also things to consider adding in the near future.
        """
        target_enc = ce.CatBoostEncoder(cols=categories)
        target_enc.fit(self.bundle[0][categories], self.bundle[0][target])
        train_ = self.bundle[0].join(self.bundle[0].transform(self.bundle[0][categories]).add_suffix('_cb'))
        test_ = self.bundle[1].join(target_enc.transform(self.bundle[1][categories]).add_suffix('_cb'))
        return train_, test_

    def add_interactions(self, numerical=True, categories=False):
        """
        generate mathematical feature interactions on numerical columns.
        multiplying, dividing and applying the crabby formula.

        also generate categorical features by combining existing categorical features
        using concatenation
        """
        if numerical:
            numbs = [column for column in self.data.select_dtypes('number') if len(self.data[column].unique()) > 10]
            cardinal = 0
            for column in numbs:
                for index in range(cardinal, len(numbs)):
                    if index != numbs.index(column):
                        self.data['Multiple{}-{}-{}'.format(index, column, numbs[index])] = self.data[column] * \
                                                                                self.data[numbs[index]]
                        self.data['Ratio{}-{}-{}'.format(index, column, numbs[index])] = self.data[column] / self.data[
                            numbs[index]]
                        self.data['Crabby{}-{}-{}'.format(index, column, numbs[index])] = (
                                (self.data[column] + self.data[numbs[index]]) / (self.data[numbs[index]] ** 2))
                    else:
                        pass
                cardinal = cardinal + 1
            cardinal = 0

        if categories:
            numbs = [column for column in self.data.select_dtypes('object').columns]
            for column in numbs:
                for index in range(cardinal, len(numbs)):
                    if index != numbs.index(column):
                        self.data['{}-{}-{}'.format(index, column, numbs[index])] = self.data[column] + '_' + self.data[
                            numbs[index]]
                    else:
                        pass
                cardinal = cardinal + 1

    def viewmissing(self, heat=False):
        if heat == False:
            mis_val = self.data.isnull().sum()
            mis_val_percent = 100 * self.data.isnull().sum() / len(self.data)
            mis_val_table = pd.concat([mis_val, mis_val_percent], axis=1)
            mis_val_table_ren_columns = mis_val_table.rename(
                columns={0: 'Missing Values', 1: '% of Total Values'})
            mis_val_table_ren_columns = mis_val_table_ren_columns[
                mis_val_table_ren_columns.iloc[:, 1] != 0].sort_values(
                '% of Total Values', ascending=False).round(1)
            if len(mis_val_table_ren_columns) > 1:
                print("Your selected dataframe has " + str(self.data.shape[1]) + " columns.\n"
                                                                                "There are " + str(
                    mis_val_table_ren_columns.shape[0]) +
                      " columns that have missing values.")
            else:
                print("No missing Values!")
            return mis_val_table_ren_columns
        elif heat:
            sns.heatmap(self.data.isnull(), cbar=False)
        else:
            print('Error 404! Lol, you set the heat wrongly though!')

    def hour_split(self, hour):
        self.data['hours {}'.format(hour)] = 0
        self.data.loc[self.data[hour] < 12, 'hours {}'.format(hour)] = 'morning'
        self.data.loc[(self.data[hour] >= 12) & (self.data[hour] < 16), 'hours {}'.format(hour)] = 'afternoon'
        self.data.loc[(self.data[hour] >= 16) & (self.data[hour] < 20), 'hours {}'.format(hour)] = 'evening'
        self.data.loc[(self.data[hour] >= 19) & (self.data[hour] <= 24), 'hours {}'.format(hour)] = 'night'

    def date_split(self, date):
        a = pd.to_datetime(self.data[date])
        self.data['weekday'] = a.dt.dayofweek
        self.data['hour'] = a.dt.hour
        self.data['year'] = a.dt.dayofyear
        self.data['day'] = a.dt.day
        self.data['quarter'] = a.dt.quarter
        self.data['is_weekend'] = 0
        self.data['is_monthend'] = 0
        self.data['is_monthstart'] = 0
        self.data.loc[(self.data['day'] >= 26), 'is_monthend'] = 1
        self.data.loc[(self.data['day'] <= 5), 'is_monthstart'] = 1
        self.data.loc[(self.data['weekday'] >= 4), 'is_weekend'] = 1

    def catcode(self, column, encoded=False):
        if encoded == False:
            self.data[column] = self.data[column].astype('category')
            self.data[column] = self.data[column].cat.codes
        else:
            self.encoded_data[column] = self.encoded_data[column].astype('category')
            self.encoded_data[column] = self.encoded_data[column].cat.codes
        return 'done'

    def quality_report(self):
        dtypes = self.data.dtypes
        nuniq = self.data.T.apply(lambda x: x.nunique(), axis=1)
        total = self.data.isnull().sum().sort_values(ascending=False)
        percent = (self.data.isnull().sum() / self.data.isnull().count() * 100).sort_values(ascending=False)
        quality_df = pd.concat([total, percent, nuniq, dtypes], axis=1, keys=['Total', 'Percent', 'Nunique', 'Dtype'])
        display(quality_df)

    def eda(self, target):
        for col in self.data.drop(target, axis=1).columns:
            if self.data[col].dtype == 'object' or self.data[col].nunique():
                xx = self.data.groupby(col)[target].value_counts().unstack(1)
                per_not_promoted = xx.iloc[:, 0] * 100 / xx.apply(lambda x: x.sum(), axis=1)
                per_promoted = xx.iloc[:, 1] * 100 / xx.apply(lambda x: x.sum(), axis=1)
                xx['%_0'] = per_not_promoted
                xx['%_1'] = per_promoted
                display(xx)

    def drop(self, column):
        recycle[column] = self.data[column].copy()
        self.data.drop(columns=column, inplace=True)

    def pick(self, column):
        self.data[column] = recycle[column]

    def evaluate(self, target, model, test_size=.25, stratify=True, metric=accuracy_score):
        if stratify:
            x_train, x_test, y_train, y_test = train_test_split(self.data.drop(target, axis=1), self.data[target],
                                                                stratify=self.data[target], test_size=test_size)
            model.fit(x_train, y_train)
            pred = model.predict(x_test)
            metric(pred, y_test)
        else:
            x_train, x_test, y_train, y_test = train_test_split(self.data.drop(target, axis=1), self.data[target],
                                                                test_size=test_size)
            model.fit(x_train, y_train)
            pred = model.predict(x_test)
            metric(pred, y_test)
        return pred

ext: py

max_stars_repo_name: onuohasilver/quickml

max_issues_repo_name: onuohasilver/quickml

max_forks_repo_name: onuohasilver/quickml

reason: import numpy


text:
import numpy as np

def accuracy(predict_tagging, true_tagging):
	if len(predict_tagging) != len(true_tagging):
		return 0, 0, 0
	cnt = 0
	for i in range(len(predict_tagging)):
		if predict_tagging[i] == true_tagging[i]:
			cnt += 1
	total_correct = cnt
	total_words = len(predict_tagging)
	if total_words == 0:
		return 0, 0, 0
	return total_correct, total_words, total_correct*1.0/total_words


text:
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.patches import Circle
import ephem
import sqlite3 as lite
from progressbar import ProgressBar

# Altitude and Azimuth of a single field at t (JD) in rad
def Fields_local_coordinate(Field_ra, Field_dec, t, Site):

    # date and time
    Site.date = t
    curr_obj = ephem.FixedBody()
    curr_obj._ra = Field_ra * np.pi / 180
    curr_obj._dec = Field_dec * np.pi / 180
    curr_obj.compute(Site)
    altitude = curr_obj.alt
    azimuth = curr_obj.az
    return altitude, azimuth

def update_moon(t, Site) :
    Moon = ephem.Moon()
    Site.date = t
    Moon.compute(Site)
    X, Y = AltAz2XY(Moon.alt, Moon.az)
    r = Moon.size / 3600 * np.pi / 180  *2
    return X, Y, r, Moon.alt

def AltAz2XY(Alt, Az) :
    X = np.cos(Alt) * np.cos(Az) * -1
    Y = np.cos(Alt) * np.sin(Az)
    #Y = Alt * 2/ np.pi
    #X = Az / (2*np.pi)

    return Y, -1*X

def Visualize(Date, PlotID = 1,FPS = 15,Steps = 20,MP4_quality = 300, Name = "LSST Scheduler Simulator.mp4", showClouds = True):

    # Import data
    All_Fields = np.loadtxt("NightDataInLIS/Constants/fieldID.lis", unpack = True)
    N_Fields   = len(All_Fields[0])

    if showClouds:
        Time_slots = np.loadtxt("NightDataInLIS/TimeSlots{}.lis".format(int(ephem.julian_date(Date))), unpack = True)
        All_Cloud_cover = np.loadtxt("NightDataInLIS/Clouds{}.lis".format(int(ephem.julian_date(Date))), unpack = True)

    Site            = ephem.Observer()
    Site.lon        = -1.2320792
    Site.lat        = -0.517781017
    Site.elevation  = 2650
    Site.pressure   = 0.
    Site.horizon    = 0.

    #Initialize date and time
    lastN_start = float(Date) -1;   lastN_end = float(Date)
    toN_start = float(Date);        toN_end = float(Date) + 1

    #Connect to the History data base
    con = lite.connect('FBDE.db')
    cur = con.cursor()

    # Prepare to save in MP4 format
    FFMpegWriter = animation.writers['ffmpeg']
    metadata = dict(title='LSST Simulation', artist='Elahe', comment='Test')
    writer = FFMpegWriter(fps=FPS, metadata=metadata)

    #Progress bar initialization

    pbar = ProgressBar()

    # Initialize plot
    Fig = plt.figure()
    if PlotID == 1:
        ax = plt.subplot(111, axisbg = 'black')
    if PlotID == 2:
        ax = plt.subplot(211, axisbg = 'black')

    unobserved, Observed_lastN, Obseved_toN,\
    ToN_History_line, last_10_History_line,\
    Horizon, airmass_horizon, S_Pole,\
    LSST,\
    Clouds\
        = ax.plot([], [], '*',[], [], '*',[], [], '*',
                  [], [], '-',[], [], '*',
                  [], [], '-',[], [], '-',[], [], 'D',
                  [], [], 'o',
                  [], [], 'o')

    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    ax.set_aspect('equal', adjustable = 'box')
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)

    # Coloring
    Horizon.set_color('white'); airmass_horizon.set_color('red')
    S_Pole.set_markersize(3);   S_Pole.set_markerfacecolor('red')
    star_size = 4

    unobserved.set_color('dimgray');        unobserved.set_markersize(star_size)
    Observed_lastN.set_color('blue');       Observed_lastN.set_markersize(star_size)
    Obseved_toN.set_color('chartreuse');    Obseved_toN.set_markersize(star_size +2)
    Clouds.set_color('white');              Clouds.set_markersize(10)

    ToN_History_line.set_color('orange');   ToN_History_line.set_lw(.5)
    last_10_History_line.set_color('red');  last_10_History_line.set_lw(.5)

    LSST.set_color('red'); LSST.set_markersize(8)

    if PlotID == 2:
        freqAX = plt.subplot(212)
        cur.execute('SELECT N_visit, Last_visit, Second_last_visit, Third_last_visit, Fourth_last_visit From FieldsStatistics')
        row = cur.fetchall()
        N_visit     = [x[0] for x in row]
        Last_visit   = [x[1] for x in row]
        Second_last_visit = [x[2] for x in row]
        Third_last_visit  = [x[3] for x in row]
        Fourth_last_visit = [x[4] for x in row]

        initHistoricalcoverage = N_visit
        for index, id in enumerate(All_Fields):
            if Last_visit[index] > toN_start:
                initHistoricalcoverage[index] -= 1
                if Second_last_visit[index] > toN_start:
                    initHistoricalcoverage[index] -= 1
                    if Third_last_visit > toN_start:
                        initHistoricalcoverage[index] -= 1

        covering,current_cover = freqAX.plot(All_Fields[0],initHistoricalcoverage,'-',[],[],'o')

        freqAX.set_xlim(0,N_Fields)
        freqAX.set_ylim(0,np.max(initHistoricalcoverage)+5)
        covering.set_color('chartreuse');   covering.set_markersize(2)
        current_cover.set_color('red');     current_cover.set_markersize(6)

    cur.execute('SELECT Night_count, T_start, T_end FROM NightSummary WHERE T_start BETWEEN (?) AND (?)',(toN_start, toN_end))
    row = cur.fetchone()
    vID = row[0]
    t_start = row[1]
    t_end   = row[2]
    t = t_start

    # Figure labels and fixed elements
    Phi = np.arange(0, 2* np.pi, 0.05)
    Horizon.set_data(1.01*np.cos(Phi), 1.01*np.sin(Phi))
    ax.text(-1.3, 0, 'West', color = 'white', fontsize = 7)
    ax.text(1.15, 0 ,'East', color = 'white', fontsize = 7)
    ax.text( 0, 1.1, 'North', color = 'white', fontsize = 7)
    airmass_horizon.set_data(np.cos(np.pi/4) * np.cos(Phi), np.cos(np.pi/4) *  np.sin(Phi))
    ax.text(-.3, 0.6, 'Acceptable airmass horizon', color = 'white', fontsize = 5, fontweight = 'bold')
    Alt, Az = Fields_local_coordinate(180, -90, t, Site)
    x, y = AltAz2XY(Alt,Az)
    S_Pole.set_data(x, y)
    ax.text(x+ .05, y, 'S-Pole', color = 'white', fontsize = 7)

    # Observed last night fields
    cur.execute('SELECT Field_id FROM Schedule WHERE ephemDate BETWEEN (?) AND (?)',(lastN_start, lastN_end))
    row = cur.fetchall()
    if row is not None:
        F1 = [x[0] for x in row]
    else:
        F1 = []

    # Tonight observation path
    cur.execute('SELECT Field_id, ephemDate FROM Schedule WHERE ephemDate BETWEEN (?) AND (?)',(toN_start, toN_end))
    row = cur.fetchall()
    if row[0][0] is not None:
        F2 = [x[0] for x in row]
        F2_timing = [x[1] for x in row]
    else:
        F2 = []; F2_timing = []

    # Sky elements
    Moon = Circle((0, 0), 0, color = 'silver', zorder = 3)
    ax.add_patch(Moon)
    Moon_text = ax.text([], [], 'Moon', color = 'white', fontsize = 7)

    with writer.saving(Fig, Name, MP4_quality) :
        for t in pbar(np.linspace(t_start, t_end, num = Steps)):

            # Find the index of the current time
            time_index = 0
            while t > F2_timing[time_index]:
                time_index += 1
            if showClouds:
                Slot_n = 0
                while t > Time_slots[Slot_n]:
                    Slot_n += 1

            visit_index = 0
            visited_field = 0

            # Object fields: F1)Observed last night F2)Observed tonight F3)Unobserved F4)Covered by clouds
            F1_X = []; F1_Y = []; F2_X = []; F2_Y = []; F3_X = []; F3_Y = []; F4_X = []; F4_Y = []

            # F1  coordinate:
            for i in F1:
                Alt, Az = Fields_local_coordinate(All_Fields[1,i-1], All_Fields[2,i-1], t, Site)
                if Alt > 0:
                    X, Y    = AltAz2XY(Alt,Az)
                    F1_X.append(X); F1_Y.append(Y)

            # F2  coordinate:
            for i,tau in zip(F2,F2_timing):
                Alt, Az = Fields_local_coordinate(All_Fields[1,i-1], All_Fields[2,i-1], t, Site)
                if Alt > 0:
                    X, Y    = AltAz2XY(Alt,Az)
                    F2_X.append(X); F2_Y.append(Y)
                    if t >= tau:
                        visit_index = len(F2_X) -1
                        visited_field = i

            # F3  coordinate:
            for i in range(0,N_Fields):
                if True:
                    Alt, Az = Fields_local_coordinate(All_Fields[1,i], All_Fields[2,i], t, Site)
                    if Alt > 0:
                        X, Y    = AltAz2XY(Alt,Az)
                        F3_X.append(X); F3_Y.append(Y)

            # F4 coordinates
            if showClouds:
                for i in range(0,N_Fields):
                    if All_Cloud_cover[i, Slot_n] == 1:
                        Alt, Az = Fields_local_coordinate(All_Fields[1,i], All_Fields[2,i], t, Site)
                    if Alt > 0:
                        X, Y    = AltAz2XY(Alt,Az)
                        F4_X.append(X); F4_Y.append(Y)

            # Update plot
            unobserved.set_data([F3_X,F3_Y])
            Observed_lastN.set_data([F1_X,F1_Y])
            Obseved_toN.set_data([F2_X[0:visit_index],F2_Y[0:visit_index]])

            ToN_History_line.set_data([F2_X[0:visit_index], F2_Y[0:visit_index]])
            last_10_History_line.set_data([F2_X[visit_index - 10: visit_index], F2_Y[visit_index - 10: visit_index]])
            LSST.set_data([F2_X[visit_index],F2_Y[visit_index]])
            Clouds.set_data([F4_X,F4_Y])

            # Update Moon
            X, Y, r, alt = update_moon(t, Site)
            Moon.center = X, Y
            Moon.radius = r
            if alt > 0:
                #Moon.set_visible(True)
                Moon_text.set_visible(True)
                Moon_text.set_x(X+.002); Moon_text.set_y(Y+.002)
            else :
                Moon.set_visible(False)
                Moon_text.set_visible(False)

            #Update coverage
            if PlotID == 2:
                Historicalcoverage = np.zeros(N_Fields)
                for i,tau in zip(F2, F2_timing):
                    if tau <= t:
                        Historicalcoverage[i -1] += 1
                    else:
                        break
                tot = Historicalcoverage + initHistoricalcoverage
                current_cover.set_data(visited_field -1,tot[visited_field -1])
                covering.set_data(All_Fields[0], tot)

            #Observation statistics
            leg = plt.legend([Observed_lastN, Obseved_toN],
                       ['Visited last night', time_index])
            for l in leg.get_texts():
                l.set_fontsize(6)
            date = ephem.date(t)
            Fig.suptitle('Top view of the LSST site on {}, GMT'.format(date))

            '''
            # progress
            perc= int(100*(t - t_start)/(t_end - t_start))
            if perc <= 100:
                print('{} %'.format(perc))
            else:
                print('100 %')
            '''
            #Save current frame
            writer.grab_frame()

'''
Date = ephem.Date('2016/09/3 12:00:00.00') # times are in UT
# Animation specifications
FPS = 10            # Frame per second
Steps = 100          # Simulation steps
MP4_quality = 300   # MP4 size and quality

PlotID = 2        # 1 for one Plot, 2 for including covering pattern
Visualize(Date, PlotID ,FPS, Steps, MP4_quality, 'Visualizations/LSST1plot{}.mp4'.format(int(ephem.julian_date(Date))), showClouds= False)
'''


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

from mesostat.utils.signals.resample import bin_data
from mesostat.utils.signals.filter import approx_decay_conv, zscore

x = np.random.normal(0,1,(2000, 3))

x = approx_decay_conv(x, 100, 1)
y = bin_data(x, 4, axis=1)
xz = zscore(x, axis=0)
yz = zscore(y, axis=0)

plt.figure()
plt.plot(xz[:, 1])
plt.plot(yz[:, 1])
plt.show()

max_stars_repo_name: HelmchenLabSoftware/mesostat-dev

max_issues_repo_name: HelmchenLabSoftware/mesostat-dev

max_forks_repo_name: HelmchenLabSoftware/mesostat-dev


text:
import json
import pickle
from pathlib import Path
from textwrap import fill

import networkx as nx
from gridengineapp import entry, GridParser

from cascade.core import getLoggers
from cascade.core.db import use_local_odbc_ini
from cascade.executor.execution_context import make_execution_context
from cascade.executor.job_definitions import job_graph_from_settings
from cascade.input_data.db.configuration import json_settings_to_frozen_settings
from cascade.input_data.db.configuration import load_settings
from cascade.input_data.db.locations import location_hierarchy
from cascade.runner.application_config import application_config
from cascade.runner.cascade_logging import logging_config

CODELOG, MATHLOG = getLoggers(__name__)

class DismodAT:
    """
    Responsible for management of settings and creation of job graphs.

    All arguments have default None, which is the typical way to
    instantiate this, unless it is under test. Using arguments here
    makes it unnecessary to use ``create_settings`` or ``load_settings``.

    Args:
        locations (nx.DiGraph): Graph of locations.
        settings (SimpleNamespace): Settings for the whole run.
        execution_context (ExecutionContext): defines the environment.
    """
    def __init__(
            self,
            locations=None,
            settings=None,
            execution_context=None,
            args=None
    ):
        self.locations = locations
        self.settings = settings
        self.execution_context = execution_context
        self.args = args
        self._job_graph = None

    @staticmethod
    def add_arguments(parser=None):
        """Add arguments to an argument parser. These arguments are relevant
        to the application but not to how it is run.

        Args:
            parser (ArgumentParser): If not supplied, a parser is created.

        Returns:
            ArgumentParser: The one that is created, or the one passed in.
        """
        if parser is None:
            parser = GridParser()
        parser.add_argument(
            "--meid", type=int,
            help="Modelable entity ID. This identifies the cause of disease.",
        )
        parser.add_argument(
            "--mvid", type=int,
            help=("Model version ID. There are multiple model versions for "
                  "each model entity"),
        )
        context_parser = parser.add_argument_group(
            "context",
            "Settings to change the process directories and environment."
        )
        context_parser.add_argument(
            "--infrastructure", action="store_true",
            help=fill("Whether we are running as infrastructure component. "
                      "Controls whether this tries to store files in the "
                      "base directory or whether it works in the"
                      "current directory."),
        )
        context_parser.add_argument(
            "--base-directory", type=Path, default=".",
            help="Directory in which to find and store files.",
        )

        pdb_parser = parser.add_argument_group(
            "debug",
            "These affect how this executes, for debugging."
        )
        pdb_parser.add_argument(
            "--no-upload", action="store_true",
            help=fill("This turns off all writing to databases, so that a run "
                      "will not affect outputs"),
        )
        pdb_parser.add_argument(
            "--db-only", action="store_true",
            help=fill("Run until it creates a Dismod-AT db file, and then quit. "
                      "If one recipe runs with --db-only then the next "
                      "recipe will fail."),
        )

        data_parser = parser.add_argument_group(
            "data",
            "Parsers to change what data is used for input."
        )
        data_parser.add_argument(
            "--settings-file", type=Path,
            help="Read settings from this file.",
        )
        data_parser.add_argument(
            "-b", "--bundle-file", type=Path,
            help=fill("The bundle normally comes from the databases but this "
                      "lets you specify a particular file as input. If this "
                      "isn't given, it will look for the bundle in the db or "
                      "in a known input directory."),
        )
        data_parser.add_argument(
            "-s", "--bundle-study-covariates-file", type=Path,
            help=fill("Read the study covariates from a separate file. "
                      "If this is not specified, then they are read from "
                      "either the databases or from the known file location."),
        )

        graph_parser = parser.add_argument_group(
            "graph",
            "Changes to the graph of work that is done."
        )
        graph_parser.add_argument(
            "--skip-cache", action="store_true",
            help=fill("Don't save bundle data to tier 3. Instead, read "
                      "it directly from tier 2."),
        )
        graph_parser.add_argument("--num-samples", type=int, help="Override number of samples.")

        sub_graph = parser.add_argument_group(
            "sub_graph",
            "These arguments select a subset of nodes to run."
        )
        sub_graph.add_argument("--location-id", type=int, help="location ID for this work")
        sub_graph.add_argument("--sex", type=str, help="sex as male, female, both")
        sub_graph.add_argument("--recipe", type=str, help="name of the recipe")
        sub_graph.add_argument("--name", type=str, help="job within the recipe")

        config = application_config()["DataLayout"]
        root_dir = Path(config["root-directory"]).resolve()
        code_log = config["code-log-directory"]
        epiviz_log = config["epiviz-log-directory"]
        log_parse = parser.add_argument_group(
            "logs",
            "Options that affect logging",
        )
        log_parse.add_argument(
            "-v", "--verbose", action="count", default=0,
            help="Increase verbosity of logging")
        log_parse.add_argument(
            "-q", "--quiet", action="count", default=0,
            help="Decrease verbosity of logging")
        log_parse.add_argument(
            "--logmod", action="append", default=[],
            help="Set logging to debug for submodule")
        log_parse.add_argument(
            "--modlevel", type=str, default="debug",
            help="Log level for specified modules")
        log_parse.add_argument(
            "--epiviz-log", type=Path, default=epiviz_log,
            help="Directory for EpiViz log")
        log_parse.add_argument(
            "--code-log", type=Path, default=code_log,
            help="Directory for code log")
        log_parse.add_argument(
            "--root-dir", type=Path, default=root_dir,
            help="Directory to use as root for logs.")
        return parser

    @staticmethod
    def job_id_to_arguments(job_id):
        # Our job_id is a custom class that knows what its arguments are.
        return job_id.arguments

    def initialize(self, args):
        self.args = args
        use_local_odbc_ini()
        logging_config(args)
        self.create_settings(args)

    def create_settings(self, args):
        # We need a sort-of-correct execution context when we first load
        # and then it gets refined after settings are loaded.
        execution_context = execution_context_without_settings(args)
        # If the application was created with settings and locations,
        # then keep them.
        if self.settings is None:
            self.settings = load_settings(
                execution_context, args.meid, args.mvid, args.settings_file
            )
        if self.locations is None:
            self.locations = location_hierarchy(
                location_set_version_id=self.settings.location_set_version_id,
                gbd_round_id=self.settings.gbd_round_id
            )
        configure_execution_context_from_settings(
            execution_context, self.settings
        )
        self.execution_context = execution_context

    def load_settings(self, args):
        # The execution context isn't part of the settings, so it is
        # rebuilt here when settings are loaded.
        self.args = args
        self.execution_context = execution_context_without_settings(args)
        base = self.execution_context.model_base_directory(0)
        setting_file = base / "settings.json"
        settings_dict = json.load(setting_file.open("r"))
        self.settings = json_settings_to_frozen_settings(settings_dict)
        location_file = base / "locations.pickle"
        self.locations = pickle.load(location_file.open("rb"))
        CODELOG.info(f"Loading settings from {setting_file} and "
                     f"locations from {location_file}")
        configure_execution_context_from_settings(
            self.execution_context, self.settings
        )

    def save_settings(self):
        base = self.execution_context.model_base_directory(0)
        base.mkdir(exist_ok=True, parents=True)
        setting_file = base / "settings.json"
        json.dump(self.settings.to_dict(), setting_file.open("w"), indent=4)
        location_file = base / "locations.pickle"
        pickle.dump(self.locations, location_file.open("wb"))
        CODELOG.info(f"Saving settings to {setting_file} "
                     f"and locations to {location_file}")

    def job_graph(self):
        if self._job_graph is None:
            self._job_graph = job_graph_from_settings(
                self.locations, self.settings, self.args, self.execution_context)
            CODELOG.info(f"Job graph has {len(self._job_graph)} nodes.")
        return self._job_graph

    def job_identifiers(self, args):
        job_graph = self.job_graph()
        nodes = job_graph.nodes

        used_queries = list()
        for search in ["location_id", "recipe", "sex", "name"]:
            if search in args and getattr(args, search) is not None:
                nodes = [n for n in nodes if getattr(n, search) == getattr(args, search)]
                used_queries.append(f"{search}={getattr(args, search)}")

        if len(nodes) == 0:
            message = f"Job chose no nodes using search {used_queries} from {list(job_graph.nodes)}"
            CODELOG.error(message)
            raise RuntimeError(message)

        sub_graph = nx.subgraph(job_graph, nodes)
        sub_graph.graph["execution_context"] = self.execution_context
        CODELOG.info(f"Execution graph has {len(sub_graph)} nodes.")
        return sub_graph

    def job(self, identifier):
        return self.job_graph().nodes[identifier]["job"]

def execution_context_without_settings(args):
    execution_context = make_execution_context(
        gbd_round_id=6
    )
    if args.infrastructure:
        execution_context.parameters.organizational_mode = "infrastructure"
    else:
        execution_context.parameters.organizational_mode = "local"

    execution_context.parameters.base_directory = args.base_directory

    if args.meid:
        execution_context.parameters.modelable_entity_id = args.meid
    if args.mvid:
        execution_context.parameters.model_version_id = args.mvid

    return execution_context

def configure_execution_context_from_settings(execution_context, settings):
    """
    This later configuration exists because the application is
    started by telling it the model version ID but *not telling it the
    modelable entity ID.* That modelable entity ID is used to decide
    where files are on disk, so we need it early.
    """
    for param in ["modelable_entity_id", "model_version_id"]:
        if hasattr(settings, "model") and hasattr(settings.model, param):
            setattr(execution_context.parameters, param,
                    getattr(settings.model, param))

def cascade_entry():
    app = DismodAT()
    entry(app)

if __name__ == "__main__":
    cascade_entry()

max_stars_repo_name: adolgert/cascade

max_issues_repo_name: adolgert/cascade

max_forks_repo_name: adolgert/cascade

reason: import networkx


text:
"""Return True if two arrays are element-wise equal within a tolerance."""
import numpy
import numpoly

from .common import implements

@implements(numpy.allclose)
def allclose(a, b, rtol=1e-5, atol=1e-8, equal_nan=False):
    """
    Return True if two arrays are element-wise equal within a tolerance.

    The tolerance values are positive, typically very small numbers.  The
    relative difference (`rtol` * abs(`b`)) and the absolute difference
    `atol` are added together to compare against the absolute difference
    between `a` and `b`.

    If either array contains one or more NaNs, False is returned.
    Infs are treated as equal if they are in the same place and of the same
    sign in both arrays.

    Args:
        a, b (numpoly.ndpoly):
            Input arrays to compare.
        rtol (float):
            The relative tolerance parameter (see Notes).
        atol : float
            The absolute tolerance parameter (see Notes).
        equal_nan : bool
            Whether to compare NaN's as equal.  If True, NaN's in `a` will be
            considered equal to NaN's in `b` in the output array.

    Returns:
        (bool):
            Returns True if the two arrays are equal within the given
            tolerance; False otherwise.

    Notes:
        If the following equation is element-wise True, then allclose returns
        True.

        absolute(`a` - `b`) <= (`atol` + `rtol` * absolute(`b`))

        The above equation is not symmetric in `a` and `b`, so that
        ``allclose(a, b)`` might be different from ``allclose(b, a)`` in some
        rare cases.

        The comparison of `a` and `b` uses standard broadcasting, which means
        that `a` and `b` need not have the same shape in order for
        ``allclose(a, b)`` to evaluate to True.  The same is true for `equal`
        but not `array_equal`.

    Examples:
        >>> x, y = numpoly.symbols("x y")
        >>> numpoly.allclose([1e10*x, 1e-7], [1.00001e10*x, 1e-8])
        False
        >>> numpoly.allclose([1e10*x, 1e-8], [1.00001e10*x, 1e-9])
        True
        >>> numpoly.allclose([1e10*x, 1e-8], [1.00001e10*y, 1e-9])
        False
        >>> numpoly.allclose([x, numpy.nan],
        ...                  [x, numpy.nan], equal_nan=True)
        True

    """
    a, b = numpoly.align_polynomials(a, b)
    for coeff1, coeff2 in zip(a.coefficients, b.coefficients):
        if not numpy.allclose(
                coeff1, coeff2, atol=atol, rtol=rtol, equal_nan=equal_nan):
            return False
    return True

max_stars_repo_name: FredrikMeyer/numpoly

max_issues_repo_name: FredrikMeyer/numpoly

max_forks_repo_name: FredrikMeyer/numpoly


text:
import copy
import itertools
import math
import random

import numpy as np
import open3d as o3d
from matplotlib import pyplot as plt
from sklearn import linear_model

import modeling.collision_model as cm
import modeling.geometric_model as gm
import basis.trimesh as trimesh
# import trimesh.sample as ts
import math_utils as mu
import basis.robot_math as rm
import o3dhelper as o3d_helper
import envloader as el

def get_objpcd(objcm, objmat4=np.eye(4), sample_num=100000, toggledebug=False):
    objpcd = np.asarray(objcm.sample_surface(nsample=sample_num,toggle_option=None))
    objpcd = trans_pcd(objpcd, objmat4)

    if toggledebug:
        objpcd = o3d_helper.nparray2o3dpcd(copy.deepcopy(objpcd))
        objpcd.paint_uniform_color([1, 0.706, 0])
        o3d.visualization.draw_geometries([objpcd])

    return objpcd

def remove_pcd_zeros(pcd):
    return pcd[np.all(pcd != 0, axis=1)]

def trans_pcd(pcd, transmat):
    homopcd = np.ones((4, len(pcd)))
    homopcd[:3, :] = pcd.T
    realpcd = np.dot(transmat, homopcd).T
    return realpcd[:, :3]

def trans_p(p, transmat=None):
    if transmat is None:
        return p
    return trans_pcd(np.asarray([p]), transmat)[0]

def show_pcd(pcd, rgba=(1, 1, 1, 1)):
    pcd_gm = gm.GeometricModel(initor=pcd)
    pcd_gm.set_rgba(rgba)
    pcd_gm.attach_to(base)

def show_pcd_withrbt(pcd, rgba=(1, 1, 1, 1), rbtx=None, toggleendcoord=False):
    rbt, rbtmg, rbtball = el.loadUr3e()
    env = el.Env_wrs(boundingradius=7.0)
    env.reparentTo(base.render)

    if rbtx is not None:
        for armname in ["lft", "rgt"]:
            tmprealjnts = rbtx.getjnts(armname)
            print(armname, tmprealjnts)
            rbt.movearmfk(tmprealjnts, armname)

    rbtmg.genmnp(rbt, toggleendcoord=toggleendcoord).reparentTo(base.render)

    pcddnp = base.pg.genpointcloudnp(pcd, colors=rgba)
    pcddnp.reparentTo(base.render)

def get_pcd_center(pcd):
    return np.array((np.mean(pcd[:, 0]), np.mean(pcd[:, 1]), np.mean(pcd[:, 2])))

def get_pcd_tip(pcd, axis=0):
    """
    get the smallest point along an axis

    :param pcd:
    :param axis: 0-x,1-y,2-z
    :return: 3D point
    """

    return pcd[list(pcd[:, axis]).index(min(list(pcd[:, axis])))]

def get_pcd_w_h(objpcd_std):
    def __sort_w_h(a, b):
        if a > b:
            return a, b
        else:
            return b, a

    return __sort_w_h(max(objpcd_std[:, 0]) - min(objpcd_std[:, 0]), max(objpcd_std[:, 1]) - min(objpcd_std[:, 1]))

def get_org_convexhull(pcd, color=(1, 1, 1), transparency=1, toggledebug=False):
    """
    create CollisionModel by pcd

    :param pcd:
    :param color:
    :param transparency:
    :return: CollisionModel
    """

    convexhull = trimesh.Trimesh(vertices=pcd)
    convexhull = convexhull.convex_hull
    obj = cm.CollisionModel(initor=convexhull, type="ball")
    if toggledebug:
        obj.set_rgba(color[0], color[1], color[2], transparency)
        obj.reparentTo(base.render)
        obj.showlocalframe()

    return obj

def get_std_convexhull(pcd, origin="center", color=(1, 1, 1), transparency=1, toggledebug=False, toggleransac=True):
    """
    create CollisionModel by pcd, standardized rotation

    :param pcd:
    :param origin: "center" or "tip"
    :param color:
    :param transparency:
    :param toggledebug:
    :param toggleransac:
    :return: CollisionModel, position of origin
    """

    rot_angle = get_rot_frompcd(pcd, toggledebug=toggledebug, toggleransac=toggleransac)
    center = get_pcd_center(pcd)
    origin_pos = np.array(center)

    pcd = pcd - np.array([center]).repeat(len(pcd), axis=0)
    pcd = trans_pcd(pcd, rm.homobuild((0, 0, 0), rm.rodrigues((0, 0, 1), -rot_angle)))

    convexhull = trimesh.Trimesh(vertices=pcd)
    convexhull = convexhull.convex_hull
    obj = cm.CollisionModel(initor=convexhull)
    obj_w, obj_h = get_pcd_w_h(pcd)

    if origin == "tip":
        tip = get_pcd_tip(pcd, axis=0)
        origin_pos = center + \
                     trans_pcd(np.array([tip]), rm.homobuild((0, 0, 0), rm.rodrigues((0, 0, 1), rot_angle)))[0]
        pcd = pcd - np.array([tip]).repeat(len(pcd), axis=0)

        convexhull = trimesh.Trimesh(vertices=pcd)
        convexhull = convexhull.convex_hull
        obj = cm.CollisionModel(initor=convexhull)

    if toggledebug:
        obj.set_rgba(color[0], color[1], color[2], transparency)
        obj.reparentTo(base.render)
        obj.showlocalframe()
        print("Rotation angle:", rot_angle)
        print("Origin:", center)
        base.pggen.plotSphere(base.render, center, radius=2, rgba=(1, 0, 0, 1))

    return obj, obj_w, obj_h, origin_pos, rot_angle

def get_rot_frompcd(pcd, toggledebug=False, toggleransac=True):
    """

    :param pcd:
    :param toggledebug:
    :param toggleransac: use ransac linear regression or not
    :return: grasp center and rotation angle
    """

    if max(pcd[:, 0]) - min(pcd[:, 0]) > max(pcd[:, 1]) - min(pcd[:, 1]):
        X = pcd[:, 0].reshape((len(pcd), 1))
        y = pcd[:, 1]
    else:
        X = pcd[:, 1].reshape((len(pcd), 1))
        y = pcd[:, 0]
    ransac = linear_model.RANSACRegressor()
    ransac.fit(X, y)
    ransac_coef = ransac.estimator_.coef_

    lr = linear_model.LinearRegression()
    lr.fit(X, y)
    lr_coef = lr.coef_

    if toggledebug:
        inlier_mask = ransac.inlier_mask_
        outlier_mask = np.logical_not(inlier_mask)

        # Predict data of estimated models
        line_X = np.arange(X.min(), X.max())[:, np.newaxis]
        line_y = lr.predict(line_X)
        line_y_ransac = ransac.predict(line_X)

        # Compare estimated coefficients
        print("Estimated coefficients (linear regression, RANSAC):", lr.coef_, ransac.estimator_.coef_)

        plt.scatter(X[inlier_mask], y[inlier_mask], color='yellowgreen', marker='.',
                    label='Inliers')
        plt.scatter(X[outlier_mask], y[outlier_mask], color='gold', marker='.',
                    label='Outliers')
        plt.plot(line_X, line_y, color='navy', linewidth=2, label='Linear regressor')
        plt.plot(line_X, line_y_ransac, color='cornflowerblue', linewidth=2,
                 label='RANSAC regressor')
        plt.legend(loc='lower right')
        plt.xlabel("Input")
        plt.ylabel("Response")
        plt.show()

    if toggleransac:
        coef = ransac_coef[0]
    else:
        coef = lr_coef[0]

    if max(pcd[:, 0]) - min(pcd[:, 0]) > max(pcd[:, 1]) - min(pcd[:, 1]):
        return math.degrees(math.atan(coef))
    else:
        return math.degrees(math.atan(1 / coef))

def reconstruct_surface(pcd, radii=[.005], toggledebug=False):
    print("---------------reconstruct surface bp---------------")
    pcd = np.asarray(pcd)
    tmmesh = o3d_helper.reconstructsurfaces_bp(pcd, radii=radii, doseparation=False)
    obj = cm.CollisionModel(initor=tmmesh)
    if toggledebug:
        obj.set_rgba(1, 1, 1, 1)
        obj.reparentTo(base.render)
    return obj

def reconstruct_surface_list(pcd, radii=[5], color=(1, 1, 1), transparency=1, toggledebug=False):
    pcd = np.asarray(pcd)
    tmmeshlist = o3d_helper.reconstructsurfaces_bp(pcd, radii=radii, doseparation=True)
    obj_list = []
    for tmmesh in tmmeshlist:
        obj = cm.CollisionModel(initor=tmmesh)
        obj_list.append(obj)
        if toggledebug:
            obj.set_rgba(color[0], color[1], color[2], transparency)
            obj.reparentTo(base.render)
    return obj_list

def get_pcdidx_by_pos(pcd, realpos, diff=10, dim=3):
    idx = 0
    distance = 100
    result_point = None
    for i in range(len(pcd)):
        point = pcd[i]
        if realpos[0] - diff < point[0] < realpos[0] + diff and realpos[1] - diff < point[1] < realpos[1] + diff:
            temp_distance = np.linalg.norm(realpos[:dim] - point[:dim])
            # print(i, point, temp_distance, distance)
            if temp_distance < distance:
                distance = temp_distance
                idx = i
                result_point = point
    return idx, result_point

def get_objpcd_withnrmls(objcm, objmat4=np.eye(4), sample_num=100000, toggledebug=False, sample_edge=False):
    objpcd_nrmls = []
    faces = objcm.objtrm.faces
    vertices = objcm.objtrm.vertices
    nrmls = objcm.objtrm.face_normals

    if sample_num is not None:
        objpcd, faceid = objcm.sample_surface(nsample=sample_num)
        objpcd = list(objpcd)
        for i in faceid:
            objpcd_nrmls.append(np.array(nrmls[i]))
    else:
        objpcd = []

    v_temp = []
    for i, face in enumerate(faces):
        for j, v in enumerate(face):
            if v not in v_temp:
                v_temp.append(v)
                objpcd.append(vertices[v])
                objpcd_nrmls.append(nrmls[i])

    if sample_edge:
        for i, face in enumerate(faces):
            for v_pair in itertools.combinations(face, 2):
                edge_plist = mu.linear_interp_3d(vertices[v_pair[0]], vertices[v_pair[1]], step=.5)
                objpcd.extend(edge_plist)
                objpcd_nrmls.extend(np.repeat([nrmls[i]], len(edge_plist), axis=0))

    objpcd = np.asarray(objpcd)
    objpcd = trans_pcd(objpcd, objmat4)
    # objpcd_nrmls = np.asarray([-n if np.dot(n, np.asarray([0, 0, 1])) < 0 else n for n in objpcd_nrmls])

    if toggledebug:
        # objpcd = o3d_helper.nparray2o3dpcd(copy.deepcopy(objpcd))
        # objpcd.paint_uniform_color([1, 0.706, 0])
        # o3d.visualization.draw_geometries([objpcd])

        # objcm.sethomomat(objmat4)
        # objcm.set_rgba(1, 1, 1, 0.7)
        # objcm.reparentTo(base.render)
        # show_pcd(objpcd, rgba=(1, 0, 0, 1))
        for i, n in enumerate(objpcd_nrmls):
            import random
            v = random.choice(range(0, 10000))
            if v == 1:
                base.pggen.plotArrow(base.render, spos=objpcd[i], epos=objpcd[i] + 10 * n)
                base.pggen.plotSphere(base.render, pos=objpcd[i], rgba=(1, 0, 0, 1))
        base.run()

    return objpcd, np.asarray(objpcd_nrmls)

def get_objpcd_partial(objcm, objmat4=np.eye(4), sample_num=100000, toggledebug=False):
    objpcd = np.asarray(ts.sample_surface(objcm.trimesh, count=sample_num))
    objpcd = trans_pcd(objpcd, objmat4)

    grid = {}
    for p in objpcd:
        x = round(p[0], 0)
        y = round(p[1], 0)
        if str((x, y)) in grid.keys():
            grid[str((x, y))].append(p)
        else:
            grid[str((x, y))] = [p]
    objpcd_new = []
    for k, v in grid.items():
        z_max = max(np.array(v)[:, 2])
        for p in v:
            objpcd_new.append([p[0], p[1], z_max])
    objpcd_new = np.array(objpcd_new)

    print("Length of org pcd", len(objpcd))
    print("Length of partial pcd", len(objpcd_new))

    if toggledebug:
        objpcd = o3d_helper.nparray2o3dpcd(copy.deepcopy(objpcd))
        objpcd.paint_uniform_color([1, 0.706, 0])
        o3d.visualization.draw_geometries([objpcd])

        objpcd_partial = o3d_helper.nparray2o3dpcd(copy.deepcopy(objpcd_new))
        objpcd_partial.paint_uniform_color([0, 0.706, 1])
        o3d.visualization.draw_geometries([objpcd_partial])

        # pcddnp = base.pg.genpointcloudnp(objpcd)
        # pcddnp.reparentTo(base.render)

    return objpcd_new

def get_objpcd_partial_bycampos(objcm, objmat4=np.eye(4), sample_num=100000, cam_pos=np.array([860, 80, 1780]),
                                toggledebug=False):
    def __get_angle(x, y):
        lx = np.sqrt(x.dot(x))
        ly = np.sqrt(y.dot(y))
        cos_angle = x.dot(y) / (lx * ly)
        angle = np.arccos(cos_angle)
        return angle * 360 / 2 / np.pi
        # return cos_angle

    def sigmoid(angle):
        return 1 / (1 + np.exp((angle - 90) / 90)) - 0.5

    objpcd = np.asarray(ts.sample_surface(objcm.trimesh, count=sample_num))
    objpcd_center = get_pcd_center(objpcd)
    face_num = len(objcm.trimesh.face_normals)

    objpcd_new = []
    area_list = objcm.trimesh.area_faces
    area_sum = sum(area_list)

    for i, n in enumerate(objcm.trimesh.face_normals):
        n = np.dot(n, objmat4[:3, :3])
        angle = __get_angle(n, np.array(cam_pos - objpcd_center))

        if angle > 90:
            continue
        else:
            objcm_temp = copy.deepcopy(objcm)
            # print(i, angle, sigmoid(angle))
            mask_temp = [False] * face_num
            mask_temp[i] = True
            objcm_temp.trimesh.update_faces(mask_temp)
            objpcd_new.extend(np.asarray(objcm_temp.sample_surface(
                nsample=int(sample_num / area_sum * area_list[i] * sigmoid(angle) * 100),
                toggle_option=None)))
    if len(objpcd_new) > sample_num:
        objpcd_new = random.sample(objpcd_new, sample_num)
    objpcd_new = np.array(objpcd_new)
    objpcd_new = trans_pcd(objpcd_new, objmat4)

    # print("Length of org pcd", len(objpcd))
    # print("Length of source pcd", len(objpcd_new))

    if toggledebug:
        objpcd = o3d_helper.nparray2o3dpcd(copy.deepcopy(objpcd))
        objpcd.paint_uniform_color([1, 0.706, 0])
        o3d.visualization.draw_geometries([objpcd])

        objpcd_partial = o3d_helper.nparray2o3dpcd(copy.deepcopy(objpcd_new))
        objpcd_partial.paint_uniform_color([0, 0.706, 1])
        o3d.visualization.draw_geometries([objpcd_partial])

        # pcddnp = base.pg.genpointcloudnp(objpcd)
        # pcddnp.reparentTo(base.render)

    return objpcd_new

def get_nrmls(pcd, camera_location=(800, -200, 1800), toggledebug=False):
    pcd_o3d = o3d_helper.nparray2o3dpcd(pcd)
    pcd_o3d.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=20, max_nn=1000))
    # for n in np.asarray(pcd.normals)[:10]:
    #     print(n)
    # print("----------------")
    o3d.geometry.PointCloud.orient_normals_towards_camera_location(pcd_o3d, camera_location=camera_location)
    # for n in np.asarray(pcd.normals)[:10]:
    #     print(n)
    pcd_nrmls = np.asarray(pcd_o3d.normals)
    pcd_nrmls = np.asarray([-n if np.dot(n, np.asarray([0, 0, 1])) < 0 else n for n in pcd_nrmls])

    if toggledebug:
        for i, n in enumerate(pcd_nrmls):
            import random
            v = random.choice(range(0, 100))
            if v == 1:
                base.pggen.plotArrow(base.render, spos=pcd[i], epos=pcd[i] + 10 * n)
        base.run()
    return pcd_nrmls

def get_plane(pcd, dist_threshold=0.2, toggledebug=False):
    pcd_o3d = o3d_helper.nparray2o3dpcd(pcd)
    plane, inliers = pcd_o3d.segment_plane(distance_threshold=dist_threshold,
                                           ransac_n=3, num_iterations=1000)
    plane_pcd = pcd[inliers]
    center = get_pcd_center(plane_pcd)
    if toggledebug:
        show_pcd(pcd[inliers], rgba=(1, 0, 0, 1))
        base.pggen.plotArrow(base.render, spos=center, epos=plane[:3] * plane[3], rgba=(1, 0, 0, 1))

    return plane[:3], plane[3]

if __name__ == '__main__':
    '''
    set up env and param
    '''
    base, env = el.loadEnv_wrs()
    objcm = el.loadObj("pentip.stl")

    # source_pcd = np.asarray(ts.sample_surface(objcm.trimesh, count=10000))
    source = o3d_helper.nparray2o3dpcd(source_pcd[source_pcd[:, 2] > 5])
    # source.paint_uniform_color([0, 0.706, 1])
    # o3d.visualization.draw_geometries([source])

    # get_objpcd_partial(objcm, sample_num=10000, toggledebug=True)

    # inithomomat = pickle.load(
    #     open(el.root + "/graspplanner/graspmap/pentip_cover_objmat4_list.pkl", "rb"))[1070]
    #
    # get_normals(get_objpcd(objcm, sample_num=10000))
    # pcd, pcd_normals = get_objpcd_withnormals(objcm, sample_num=100000)
    # for i, p in enumerate(pcd):
    #     base.pggen.plotArrow(base.render, spos=p, epos=p + 10 * pcd_normals[i])
    # base.run()
    get_objpcd_partial_bycampos(objcm, sample_num=10000, toggledebug=True)

    # pcd = pickle.load(open(el.root + "/dataset/pcd/a_lft_0.pkl", "rb"))
    # amat = pickle.load(open(el.root + "/camcalib/data/phoxi_calibmat_0117.pkl", "rb"))
    # pcd = transform_pcd(remove_pcd_zeros(pcd), amat)
    # print(len(pcd))
    # obj = get_org_surface(pcd)
    # base.run()

max_stars_repo_name: takuya-ki/wrs

max_issues_repo_name: takuya-ki/wrs

max_forks_repo_name: takuya-ki/wrs


text:
import numpy as np
import pandas as pd
import sklearn
from sklearn import datasets,svm,metrics
from sklearn.model_selection import KFold
from sklearn.neighbors import KNeighborsClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import *

def reduce_features(solution, features):
    selected_elements_indices = np.where(solution ==1)[0]
    reduced_features = features[:, selected_elements_indices]
    return reduced_features

def classification_accuracy(labels, predictions):
    correct = np.where(labels == predictions)[0]
    accuracy = correct.shape[0]/labels.shape[0]
    return accuracy

def cal_pop_fitness(pop, features, labels, train_indices,val_indices,classifier):
    test_accuracies = np.zeros(pop.shape[0])
    val_accuracies = np.zeros(pop.shape[0])
    idx = 0

    val_pop_pred = np.zeros(shape=(pop.shape[0],val_indices.shape[0]))
    for i,curr_solution in enumerate(pop):

        reduced_features = reduce_features(curr_solution, features)
        train_data = reduced_features[train_indices, :]
        val_data=reduced_features[val_indices,:]

        train_labels = labels[train_indices]
        val_labels=labels[val_indices]
        if classifier=='SVM':
          SV_classifier = sklearn.svm.SVC(kernel='rbf',gamma='scale',C=5000)
          SV_classifier.fit(X=train_data, y=train_labels)
          val_predictions = SV_classifier.predict(val_data)
          val_accuracies[idx] = classification_accuracy(val_labels, val_predictions)
          idx = idx + 1
        elif classifier == 'KNN':
          knn=KNeighborsClassifier(n_neighbors=8)
          knn.fit(train_data,train_labels)
          val_predictions=knn.predict(val_data)
          val_accuracies[idx]=classification_accuracy(val_labels,predictions)
          idx = idx + 1
        else :
          mlp = MLPClassifier()
          mlp.fit(train_data,train_labels)
          val_predictions=mlp.predict(val_data)
          val_accuracies[idx]=classification_accuracy(val_labels,predictions)
          idx = idx + 1
        val_pop_pred[i] = val_predictions

    return val_accuracies,val_pop_pred

def get_vector(labels,pop_preds):
    vector = np.zeros(shape=(pop_preds.shape[0],4))
    for i in range(pop_preds.shape[0]):
        preds = pop_preds[i]
        acc = classification_accuracy(labels,preds)
        pre = precision_score(labels,preds,average="macro")
        rec = recall_score(labels,preds,average="macro")
        f1 = f1_score(labels,preds,average="macro")
        vector[i] = np.array([acc,pre,rec,f1])
    return vector

def OBHSA(data_inputs,data_outputs,input_chromosome,
          HMCR=0.9,
          PAR=0.35,
          classifier="SVM",
          num_generations = 10 #Number of generations in each fold
          ):
    popSize = input_chromosome.shape[0]
    print("\nOPPOSITION-BASED HARMONY SEARCH:\n")
    population_output = np.zeros(shape = (popSize,data_inputs.shape[1],5))
    vector_output = np.zeros(shape=(popSize,4,5)) #4 because acc,pre,rec,f1; and 5 because folds=5

    num_samples = data_inputs.shape[0]
    num_feature_elements = data_inputs.shape[1]
    HM_shape=(popSize,num_feature_elements)
    #harmony_memory=np.random.randint(low=0,high=2,size=HM_shape)
    NCHV = np.ones((1, num_feature_elements))
    best_outputs = []
    best_opp_outputs = []

    kf=KFold(5,True,random_state=1)
    fold=0
    for train_indices,test_val_indices in kf.split(data_inputs):
      print("Fold : ",fold+1)
      val_indices,test_indices=train_test_split(test_val_indices,test_size=0.5,shuffle=True,random_state=8)
      best_test_outputs=[]

      #harmony_memory=np.random.randint(low=0,high=2,size=HM_shape)
      harmony_memory = input_chromosome[:,:,fold]
      opposite_memory=1-harmony_memory
      total_memory=np.concatenate((harmony_memory,opposite_memory),axis=0)
      total_fitness,_ = cal_pop_fitness(total_memory,data_inputs,data_outputs,train_indices,val_indices,classifier)
      fit_ind = np.argpartition(total_fitness, -popSize)[-popSize:]
      harmony_memory=total_memory[fit_ind,:]

      gen_fit = np.array([-1])
      for currentIteration in range(num_generations):
          NCHV = np.ones((1, num_feature_elements))
          print("Generation : ", currentIteration+1)

          fitness,val_pop_preds=cal_pop_fitness(harmony_memory,data_inputs,data_outputs,train_indices,val_indices,classifier)
          best_outputs.append(np.max(fitness))
          print("Best validation result : ", max(best_outputs))

          if max(fitness)>max(gen_fit):
              gen_fit = fitness
              gen_labels = data_outputs[val_indices]
              gen_preds = val_pop_preds

          for i in range(num_feature_elements):
              ran = np.random.rand()
              if ran < HMCR:
                  index = np.random.randint(0, popSize)
                  NCHV[0, i] = harmony_memory[index, i]
                  pvbran = np.random.rand()
                  if pvbran < PAR:
                      pvbran1 = np.random.rand()
                      result = NCHV[0, i]
                      if pvbran1 < 0.5:
                          result =1-result

              else:
                  NCHV[0, i] = np.random.randint(low=0,high=2,size=1)

          new_fitness,_ = cal_pop_fitness(NCHV, data_inputs, data_outputs, train_indices, val_indices,classifier)
          if new_fitness > min(fitness):
              min_fit_idx = np.where(fitness == min(fitness))
              harmony_memory[min_fit_idx, :] = NCHV
              fitness[min_fit_idx] = new_fitness

          opp_NCHV=1-NCHV
          new_opp_fitness,_=cal_pop_fitness(opp_NCHV,data_inputs, data_outputs, train_indices, val_indices,classifier)
          if new_opp_fitness > min(fitness):
              min_fit_idx = np.where(fitness == min(fitness))
              harmony_memory[min_fit_idx, :] = opp_NCHV
              fitness[min_fit_idx] = new_opp_fitness

      fitness,_ = cal_pop_fitness(harmony_memory, data_inputs, data_outputs, train_indices,val_indices,classifier)

      best_match_idx = np.where(fitness == np.max(fitness))[0]
      best_match_idx = best_match_idx[0]

      best_solution = harmony_memory[best_match_idx, :]
      best_solution_indices = np.where(best_solution == 1)[0]
      best_solution_num_elements = best_solution_indices.shape[0]
      best_solution_fitness = np.max(fitness)

      #print("best_match_idx : ", best_match_idx)
      #print("best_solution : ", best_solution)
      #print("Selected indices : ", best_solution_indices)
      print("Number of selected elements : ", best_solution_num_elements)

      vector_output[:,:,fold] = get_vector(gen_labels,gen_preds)
      population_output[:,:,fold] = harmony_memory

      fold=fold+1


text:
from __future__ import print_function
import boto3
import traceback
import datetime
import numpy as np

batch_client = boto3.client('batch')
s3_client = boto3.client('s3')
sqs_client = boto3.client('sqs')

def gradient_handler(event, context):

    try:
        iteration = event['iterator']['count']
        batchsize = event['batchsize']
        num_shots = event['num_shots']
        num_nodes = event['num_nodes']
        num_cores = event['omp_num_threads']
        memory = event['memory']
        bucket = event['bucket_name']
        partial_gradient_path = event['partial_gradient_path']
        full_gradient_path = event['full_gradient_path']
        gradient_name = event['gradient_name']
        model_path = event['model_path']
        velocity_name = event['velocity_name']
        water_name = event['water_name']
        variable_path = event['variable_path']
        variable_name = event['variable_name']
        script_path = event['script_path']
        script_name = event['script_name']
        data_path = event['data_path']
        data_name = event['data_name']
        sqs_queue = event['sqs_queue']
        num_chunks = event['num_chunks']
        step_length = event['step_length']
        step_scaling = event['step_scaling']
        container = event['container']
        batch_queue = event['batch_queue']
        num_retry = event['num_retry']
        thread_pinning = event['thread_pinning']
        instance_type = event['instance_type']
        user_id = event['user_id']
        #min_index = event['min_index']
        #max_index = event['max_index']

        # Purge queues for new iteration
        queue_list = sqs_client.list_queues()
        print(queue_list['QueueUrls'])
        for j in range(1, num_chunks+1):
            queue_name = sqs_queue + str(j)
            print(queue_name)
            for url in queue_list['QueueUrls']:
                if url[-len(queue_name):] == queue_name:
                    queue_url = url
            try:
                sqs_client.purge_queue(QueueUrl=queue_url)
            except:
                raise Exception('Error purging queue')

        # shot_range = np.arange(min_index, max_index+1,step=2)
        #shot_range = np.arange(min_index, max_index+1)
        #shot_indices = np.random.permutation(shot_range)[0:batchsize]

        for i in range(batchsize):

            # Register job definition for current iteration
            job_def = batch_client.register_job_definition(
                jobDefinitionName='isotropic-multi-node-job-' + str(i),
                type='multinode',
                nodeProperties={
                    'numNodes': num_nodes,
                    'mainNode': 0,
                    'nodeRangeProperties': [
                        {
                            'targetNodes': '0:' + str(num_nodes-1),
                            'container': {
                                'image': container,
                                'vcpus': num_cores,
                                'memory': memory,
                                'environment': [
                                    {
                                        'name': 'JUDI_ITERATION',
                                        'value': str(iteration)
                                    },
                                    {
                                        'name': 'NUM_SHOTS',
                                        'value': str(num_shots)
                                    },
                                    {
                                        'name': 'BATCHSIZE',
                                        'value': str(batchsize)
                                    },
                                    {
                                        'name': 'NUM_CHUNKS',
                                        'value': str(num_chunks)
                                    },
                                    {
                                        'name': 'NUM_NODES',
                                        'value': str(num_nodes)
                                    },
                                    {
                                        'name': 'NUM_CORES',
                                        'value': str(num_cores)
                                    },
                                    {
                                        'name': 'S3_BUCKET',
                                        'value': bucket
                                    },
                                    {
                                        'name': 'GRAD_PATH_PARTIAL',
                                        'value': partial_gradient_path
                                    },
                                    {
                                        'name': 'GRAD_PATH_FULL',
                                        'value': full_gradient_path
                                    },
                                    {
                                        'name': 'GRAD_NAME',
                                        'value': gradient_name
                                    },
                                    {
                                        'name': 'MODEL_PATH',
                                        'value': model_path
                                    },
                                    {
                                        'name': 'VELOCITY_NAME',
                                        'value': velocity_name
                                    },
                                    {
                                        'name': 'WATER_NAME',
                                        'value': water_name
                                    },
                                    {
                                        'name': 'VARIABLE_PATH',
                                        'value': variable_path
                                    },
                                    {
                                        'name': 'VARIABLE_NAME',
                                        'value': variable_name
                                    },
                                    {
                                        'name': 'SCRIPT_PATH',
                                        'value': script_path
                                    },
                                    {
                                        'name': 'SCRIPT_NAME',
                                        'value': script_name
                                    },
                                    {
                                        'name': 'DATA_PATH',
                                        'value': data_path
                                    },
                                    {
                                        'name': 'DATA_NAME',
                                        'value': data_name
                                    },
                                    {
                                        'name': 'SQS_QUEUE',
                                        'value': sqs_queue
                                    },
                                    {
                                        'name': 'STEP_LENGTH',
                                        'value': str(step_length)
                                    },
                                    {
                                        'name': 'STEP_SCALING',
                                        'value': str(step_scaling)
                                    },
                                    {
                                        'name': 'THREAD_PINNING',
                                        'value': thread_pinning
                                    },
                                    {
                                        'name': 'USER_ID',
                                        'value': user_id
                                    },
                                ],
                                'volumes': [
                                    {
                                        'host': {
                                            'sourcePath': '/efs'
                                        },
                                        'name': 'efs'
                                    },
                                ],
                                'mountPoints': [
                                    {
                                        'containerPath': '/efs',
                                        'sourceVolume': 'efs'
                                    },
                                ],
                                'instanceType': instance_type
                            }
                        },
                    ]
                },
                retryStrategy={
                    'attempts': num_retry
                }
            )

            # Submit job
            revision = job_def['revision']
            response = batch_client.submit_job(
                jobName='multi-node-batch-' + str(i) + '-iter-' + str(iteration),
                jobQueue=batch_queue,
                jobDefinition='isotropic-multi-node-job-' + str(i) + ':' + str(revision),
            )

            # Save job id
            if i == 0:
                job_ids = response['jobId']
            else:
                job_ids += ('&' + response['jobId'])

        # Write job ids to s3 (too big to save in event)
        timestamp = datetime.datetime.now()
        timestamp = str(timestamp).replace(' ','')  # convert to string and remove spaces
        key = 'logs/multinode_job_ids_' + timestamp
        s3_client.put_object(Body=job_ids.encode(), Bucket=bucket, Key=key)
        event['jobIdFile'] = key

        return event

    except Exception as e:
        traceback.print_exc()
        raise e

max_stars_repo_name: slimgroup/ServerlessImagingAWS

max_issues_repo_name: slimgroup/ServerlessImagingAWS

max_forks_repo_name: slimgroup/ServerlessImagingAWS


text:
"""MIT - CSAIL - Gifford Lab - seqgra

Helper class for functions operating on DNA

@author: Konstantin Krismer
"""
import re
import logging
from typing import List

import numpy as np

class DNAHelper:
    @staticmethod
    def convert_dense_to_one_hot_encoding(seq: str):
        nt_to_num = dict({"A": 0, "C": 1, "G": 2, "T": 3})
        seq = list(seq.upper())
        one_hot_encoded_seq = np.zeros((len(seq), len(nt_to_num)))
        for i, char in enumerate(seq):
            if char in nt_to_num:
                one_hot_encoded_seq[i, nt_to_num[char]] = 1
            else:
                one_hot_encoded_seq[i, :] = 0.25
        return one_hot_encoded_seq

    @staticmethod
    def convert_one_hot_to_dense_encoding(seq: str):
        densely_encoded_seq = ["N"] * seq.shape[0]
        for i in range(seq.shape[0]):
            if all(seq[i, :] == [1, 0, 0, 0]):
                densely_encoded_seq[i] = "A"
            elif all(seq[i, :] == [0, 1, 0, 0]):
                densely_encoded_seq[i] = "C"
            elif all(seq[i, :] == [0, 0, 1, 0]):
                densely_encoded_seq[i] = "G"
            elif all(seq[i, :] == [0, 0, 0, 1]):
                densely_encoded_seq[i] = "T"
        return "".join(densely_encoded_seq)

    @staticmethod
    def check_sequence(seqs: List[str]) -> bool:
        logger = logging.getLogger(__name__)
        is_valid: bool = True
        for seq in seqs:
            if not re.match("^[actgnNACGT]*$", seq):
                logger.warning("example with invalid DNA sequence "
                               "(only 'A', 'C', 'G', 'T' allowed): %s", seq)
                is_valid = False
        return is_valid

max_stars_repo_name: gifford-lab/seqgra

max_issues_repo_name: gifford-lab/seqgra

max_forks_repo_name: gifford-lab/seqgra


text:
from __future__ import absolute_import, division

from functools import cmp_to_key
import logging
from math import ceil, sqrt

from fastcluster import linkage
from matplotlib import pyplot as plt
import numpy as np
from scipy.cluster.hierarchy import dendrogram
import sklearn

import pyproteome as pyp

from . import clusterer

LOGGER = logging.getLogger('pyp.cluster.plot')
COLOR_MAP = plt.cm.rainbow
CORR_COLOR_MAP = plt.cm.coolwarm

def hierarchical_clusters(data, y_pred):
    ds = data['ds']
    z = data['z']

    ss = sorted(
        set(y_pred),
        key=cmp_to_key
        (
            lambda x, y:
            np.corrcoef(
                ds[y_pred == x].mean(axis=0),
                ds[y_pred == y].mean(axis=0),
            )[0, 1]
        )
    )

    mapping = {
        val: key
        for key, val in enumerate(ss)
    }

    link = linkage(
        np.array([
            z[y_pred.values == x].mean(axis=0)
            for x in ss
        ])
    )

    leaves = dendrogram(
        link,
        labels=[str(x) for x in ss],
        no_plot=True,
    )['leaves']

    mapping = {
        val: key
        for key, val in enumerate(leaves)
    }

    new_ind = np.argsort(np.vectorize(lambda x: mapping[x])(y_pred))
    cn = np.corrcoef(data['data'].values[new_ind])

    cmapping = {
        cluster_n: index
        for index, cluster_n in enumerate(
            sorted(mapping, key=lambda x: mapping[x])
        )
    }

    return mapping, new_ind, cmapping, cn

def cluster_corrmap(
    data, y_pred,
    colorbar=True,
    f=None,
    ax=None,
    div_scale=None,
    show_names=None,
):
    if ax is None:
        f, ax = plt.subplots(figsize=(13, 12))
    elif f is None:
        f = ax.get_figure()

    if div_scale is None:
        div_scale = max([data['data'].shape[0] // 1000, 1])

    if len(set(y_pred)) > 1:
        mapping, new_ind, _, cn = hierarchical_clusters(data, y_pred)
        sorted_clusters = sorted(mapping, key=lambda x: mapping[x])
        y_pred = y_pred.values[new_ind]
    else:
        sorted_clusters = sorted(set(y_pred))
        cn = data['c']

    if div_scale > 1:
        cn = cn[::div_scale, ::div_scale]

    mesh = ax.pcolormesh(
        cn,
        cmap=CORR_COLOR_MAP,
        vmin=-1, vmax=1,
    )

    if show_names:
        sorted_clusters = [
            i
            for i in sorted_clusters
            if i in show_names
        ]

    for cluster_n in sorted_clusters:
        ind = np.arange(
            0,
            data['data'].shape[0],
        )[y_pred == cluster_n]
        xy = np.median(ind) / div_scale

        if show_names and isinstance(show_names, dict):
            cluster_n = show_names.get(cluster_n)

        ax.text(
            s=str(cluster_n),
            x=xy,
            y=xy,
            # fontsize=sqrt(ind.shape[0]),
            color='k',
            ha='center',
            va='center',
        )

    if colorbar:
        f.colorbar(mesh)

    return f

def plot_cluster(
    data, y_pred, cluster_n,
    f=None,
    ax=None,
    div_scale=None,
    ylabel=True,
    title=None,
    color=None,
):
    if ax is None:
        f, ax = plt.subplots(
            dpi=pyp.DEFAULT_DPI,
            figsize=(8, 4),
        )
    elif f is None:
        f = ax.get_figure()

    z = data['z']
    names = data['names']

    dad = z[y_pred.values == cluster_n]
    n = dad.shape[0]

    if div_scale is None:
        div_scale = max([dad.shape[0] // 100, 1])

    if div_scale > 1:
        dad = dad[::div_scale]

    prev_ind = 0
    means = np.mean(dad, axis=0)

    for ind, v in enumerate(data['classes']):
        n_v, p_v = None, None

        if ind > 0:
            p_v = data['classes'][ind - 1]

        if ind + 1 < len(data['classes']):
            n_v = data['classes'][ind + 1]

        if v != p_v and p_v is not None:
            ax.axvline(
                x=ind - .5,
                color='k',
                linestyle=':',
            )

        if v == n_v:
            continue

        ax.plot(
            np.tile(
                np.arange(prev_ind, ind + 1),
                (dad.shape[0], 1),
            ).T,
            dad[:, prev_ind:ind + 1].T,
            alpha=min([5 / dad.shape[0], 1]),
            color=color or COLOR_MAP(cluster_n / len(set(y_pred))),
        )

        ax.plot(
            range(prev_ind, ind + 1),
            means[prev_ind:ind + 1],
            color='k',
        )

        prev_ind = ind + 1

    if ylabel:
        ax.set_ylabel('Z-scored Change')

    ax.set_facecolor((.9,) * 3)
    ax.set_xticks(range(dad.shape[1]))
    ax.set_xticklabels(names, rotation=45, horizontalalignment='right')
    ax.set_title(
        'Cluster #{}, N={}'.format(cluster_n, n)
        if title is None else
        title
    )

    return ax

def plot_all_clusters(
    data, y_pred,
    cols=4,
):
    ss = sorted(set(y_pred))
    rows = int(np.ceil(len(ss) / cols))
    f, axes = plt.subplots(
        rows, cols,
        figsize=(cols * 3, rows * 3),
        sharey=True,
        sharex=True,
    )

    ax_iter = iter(axes.ravel())

    for ind, (cluster_n, ax) in enumerate(zip(ss, ax_iter)):
        plot_cluster(
            data, y_pred, cluster_n,
            ax=ax,
            ylabel=ind % cols == 0,
        )

    for ax in ax_iter:
        ax.axis('off')

    f.tight_layout(h_pad=1)

    return f

def show_cluster(
    data, y_pred,
    seq=None,
    protein=None,
    ylabel=True,
    f=None,
    ax=None,
    color=None,
    div_scale=None,
):
    if ax is None:
        f, ax = plt.subplots(figsize=(6, 6))
    elif f is None:
        f = ax.get_figure()

    ds = data['ds']
    z = data['z']
    dp = ds.copy()
    mod = ''
    mask = None

    if seq is not None:
        mask = dp['Sequence'] == seq
        mod = dp[mask].iloc[0]['Modifications']
        protein = ' / '.join(dp.filter(series=mask).genes)
        cluster = y_pred[mask].iloc[0]
    elif protein is not None:
        mask = dp['Proteins'] == protein
        mod = dp[mask].iloc[0]['Modifications']
        protein = ' / '.join(dp.filter(series=mask).genes)
        cluster = y_pred[mask].iloc[0]

    mod_str = mod.get_mods([(None, 'Phospho')]).__str__(prot_index=0)

    title = '{}{}'.format(
        (protein[:15] + ' ...') if len(protein) > 15 else protein,
        ' {}'.format(mod_str) if mod_str else '',
    )
    plot_cluster(
        data, y_pred, cluster,
        ax=ax,
        color=color,
        div_scale=div_scale,
        ylabel=ylabel,
        title=title,
    )

    if mask is not None:
        z_m = z[mask]
        prev_ind = 0

        for ind, v in enumerate(data['classes']):
            n_v = None

            if ind + 1 < len(data['classes']):
                n_v = data['classes'][ind + 1]

            if v == n_v:
                continue

            ax.plot(
                np.tile(
                    np.arange(prev_ind, ind + 1),
                    (z_m.shape[0], 1),
                ).T,
                z_m[:, prev_ind:ind + 1].T,
                color='k',
                linestyle=':',
                linewidth=5,
            )

            prev_ind = ind + 1

    dp.psms = dp[y_pred == cluster]

    return dp, f, ax

def show_peptide_clusters(
    data, y_pred,
    filters,
    new_colors=False,
    div_scale=None,
    cols=4,
):
    rows = int(ceil(len(filters) / cols))
    f, axes = plt.subplots(
        rows, cols,
        figsize=(cols * 2, rows * 2),
        sharex=True,
        sharey=True,
    )
    ax_iter = iter(axes.ravel())

    clusters = [
        sorted(
            set(
                y_pred[
                    data['ds']['Sequence'] == fil['seq']
                    if 'seq' in fil else
                    data['ds']['Proteins'] == fil['protein']
                ]
            )
        )
        for fil in filters
    ]

    clusters = [
        i[0]
        for i in clusters
        if i
    ]

    ss = sorted(set(clusters))

    colors = [
        COLOR_MAP(ss.index(cluster_n) / len(ss))
        if new_colors else
        COLOR_MAP(cluster_n / len(set(y_pred)))
        for cluster_n in clusters
    ]

    for ind, (fil, color, ax) in enumerate(zip(filters, colors, ax_iter)):
        show_cluster(
            data, y_pred,
            ax=ax,
            ylabel=ind % cols == 0,
            color=color,
            div_scale=div_scale,
            save=False,
            **fil
        )

    for ax in ax_iter:
        ax.axis('off')

    return f, axes

def pca(data):
    classes = data['classes']
    z = data['z']
    names = data['names']

    f, ax = plt.subplots(1, 1, figsize=(4, 4))
    x = sklearn.decomposition.PCA().fit_transform(z.T)

    for ind, label in enumerate(data['labels']):
        ax.scatter(
            x[classes == ind, 0],
            x[classes == ind, 1],
            label=label,
        )

    offset = 0
    for ind, column in enumerate(names):
        ax.text(
            x=x[ind, 0] + offset,
            y=x[ind, 1] + offset,
            s=column,
        )

    ax.legend()
    ax.set_title('PCA {}'.format(data['ds'].name))
    ax.set_xlabel('Component 1')
    ax.set_ylabel('Component 2')

    return f, ax

def cluster_range(
    data,
    min_clusters=2, max_clusters=20,
    cols=3,
):
    clusters = range(min_clusters, max_clusters + 1)

    rows = int(np.ceil(len(clusters) / cols))

    f, axes = plt.subplots(
        rows,
        cols,
        figsize=(4 * cols, 4 * rows),
    )

    for n, ax in zip(clusters, axes.ravel()):
        LOGGER.info('Generating plots for cluster {}'.format(n))

        _, y_pred = clusterer.cluster(
            data,
            n_clusters=n,
        )

        pyp.cluster.plot.cluster_corrmap(
            data, y_pred,
            ax=ax,
            colorbar=False,
        )

        ax.set_title('{} Clusters'.format(n))

    for ax in axes.ravel()[len(clusters):]:
        ax.axis('off')

    return f, axes

max_stars_repo_name: white-lab/pyproteome

max_issues_repo_name: white-lab/pyproteome

max_forks_repo_name: white-lab/pyproteome

reason: import numpy,from scipy