text:
"""
Utils for doing data analysis

"""
import pandas as pd
import numpy as np

def impute_empty_years(
    yearly_stats: pd.DataFrame, min_year: int = None, max_year: int = None
) -> pd.DataFrame:
    """
    Imputes zero values for years without data

    Args:
        yearly_stats: A dataframe with a 'year' column and other columns with data
        min_year: Lower bound for years to keep; can be smaller than yearly_stats.year.min()
        max_year: Higher bound for years; can be larger than yearly_stats.year.min()

    Returns:
        A data frame with imputed 0s for years with no data
    """
    min_year, max_year = set_def_min_max_years(yearly_stats, min_year, max_year)
    return (
        pd.DataFrame(data={"year": range(min_year, max_year + 1)})
        .merge(yearly_stats, how="left")
        .fillna(0)
        .astype(yearly_stats.dtypes)
    )

def set_def_min_max_years(df: pd.DataFrame, min_year: int, max_year: int) -> (int, int):
    """Set the default values for min and max years"""
    if min_year is None:
        min_year = df.year.min()
    if max_year is None:
        max_year = df.year.max()
    return min_year, max_year

### GtR specific utils

def gtr_deduplicate_projects(gtr_docs: pd.DataFrame) -> pd.DataFrame:
    """
    Deduplicates projects that have the same title and description. This can be used
    to report the overall funding amount of the whole project when it has
    received funding in separate installments across different years.

    Args:
        gtr_docs: A dataframe with columns with GtR project 'title', 'description',
            'amount' (research funding) and other data

    Returns:
        A dataframe where projects with the exact same title and description have
        been merged and their funding has been summed up
    """
    gtr_docs_summed_amounts = (
        gtr_docs.groupby(["title", "description"])
        .agg(amount=("amount", "sum"))
        .reset_index()
    )
    # Add the summed up amounts to the project and keep the earliest instance
    # of the duplicates
    return (
        gtr_docs.drop("amount", axis=1)
        .merge(gtr_docs_summed_amounts, on=["title", "description"], how="left")
        .sort_values("start")
        .drop_duplicates(["title", "description"], keep="first")
        .reset_index(drop=True)
        # Restore previous column order
    )[gtr_docs.columns]

def gtr_funding_per_year(
    gtr_docs: pd.DataFrame, min_year: int = None, max_year: int = None
) -> pd.DataFrame:
    """
    Given a table with projects and their funding, return an aggregation by year

    Args:
        gtr_docs: A dataframe with columns for 'start', 'project_id' and 'amount'
            (research funding) among other project data
        min_year: Lower bound for years to keep
        max_year: Higher bound for years to keep

    Returns:
        A dataframe with the following columns:
            'year',
            'no_of_projects' - number of new projects in a given year,
            'amount_total' - total amount of research funding in a given year,
            'amount_median' - median project funding in a given year

    """
    # Convert project start dates to years
    gtr_docs = gtr_docs.copy()
    gtr_docs["year"] = pd.to_datetime(gtr_docs.start).dt.year
    # Set min and max years for aggregation
    min_year, max_year = set_def_min_max_years(gtr_docs, min_year, max_year)
    # Group by year
    return (
        gtr_docs.groupby("year")
        .agg(
            # Number of new projects in a given year
            no_of_projects=("project_id", "count"),
            # Total amount of research funding in a given year
            amount_total=("amount", "sum"),
            # Median project funding in a given year
            amount_median=("amount", np.median),
        )
        .reset_index()
        # Limit results between min and max years
        .query(f"year>={min_year}")
        .query(f"year<={max_year}")
        # Convert to thousands
        .assign(
            amount_total=lambda x: x.amount_total / 1000,
            amount_median=lambda x: x.amount_median / 1000,
        )
        # Add zero values for years without data
        .pipe(impute_empty_years, min_year=min_year, max_year=max_year)
    )

def gtr_get_all_timeseries(
    gtr_docs: pd.DataFrame, min_year: int = None, max_year: int = None
) -> pd.DataFrame:
    """
    Calculates all typical time series from a list of GtR projects and return
    as one combined table

    Args:
        gtr_docs: A dataframe with columns for 'start', 'project_id' and 'amount'
            (research funding) among other project data
        min_year: Lower bound for years to keep
        max_year: Higher bound for years to keep

    Returns:
        Dataframe with columns for 'year', 'no_of_projects', 'amount_total' and
        'amount_median'
    """
    # Deduplicate projects. This is used to report the number of new projects
    # started each year, accounting for cases where the same project has received
    # additional funding in later years
    gtr_docs_dedup = gtr_deduplicate_projects(gtr_docs)
    # Number of new projects per year
    time_series_projects = gtr_funding_per_year(gtr_docs_dedup, min_year, max_year)[
        ["year", "no_of_projects"]
    ]
    # Amount of research funding per year (note: here we use the non-duplicated table,
    # to account for additional funding for projects that might have started in earlier years
    time_series_funding = gtr_funding_per_year(gtr_docs, min_year, max_year)
    # Join up both tables
    time_series_funding["no_of_projects"] = time_series_projects["no_of_projects"]
    return time_series_funding

### Crunchbase specific utils

def cb_orgs_founded_per_year(
    cb_orgs: pd.DataFrame, min_year: int = None, max_year: int = None
) -> pd.DataFrame:
    """
    Calculates the number of Crunchbase organisations founded in a given year

    Args:
        cb_orgs: A dataframe with columns for 'id' and 'founded_on' among other data
        min_year: Lower bound for years to keep
        max_year: Higher bound for years to keep

    Returns:
        A dataframe with the following columns:
            'year',
            'no_of_orgs_founded' - number of new organisations founded in a given year
    """
    # Remove orgs that don't have year when they were founded
    cb_orgs = cb_orgs[-cb_orgs.founded_on.isnull()].copy()
    # Convert dates to years
    cb_orgs["year"] = pd.to_datetime(cb_orgs.founded_on).dt.year
    # Set min and max years for aggregation
    min_year, max_year = set_def_min_max_years(cb_orgs, min_year, max_year)
    # Group by year
    return (
        cb_orgs.groupby("year")
        .agg(no_of_orgs_founded=("id", "count"))
        .reset_index()
        .pipe(impute_empty_years, min_year=min_year, max_year=max_year)
    )

def cb_investments_per_year(
    cb_funding_rounds: pd.DataFrame, min_year: int = None, max_year: int = None
) -> pd.DataFrame:
    """
    Aggregates the raised investment amount and number of deals across all orgs

    Args:
        cb_funding_rounds: A dataframe with columns for 'funding_round_id', 'raised_amount_usd'
            'raised_amount_gbp' and 'announced_on' among other data
        min_year: Lower bound for years to keep
        max_year: Higher bound for years to keep

    Returns:
        A dataframe with the following columns:
            'year',
            'no_of_rounds' - number of funding rounds (deals) in a given year
            'raised_amount_usd_total' - total raised investment (USD) in a given year
            'raised_amount_gbp_total' - total raised investment (GBP) in a given year
    """
    # Convert dates to years
    cb_funding_rounds["year"] = pd.to_datetime(cb_funding_rounds.announced_on).dt.year
    # Set min and max years for aggregation
    min_year, max_year = set_def_min_max_years(cb_funding_rounds, min_year, max_year)
    # Group by year
    return (
        cb_funding_rounds.groupby("year")
        .agg(
            no_of_rounds=("funding_round_id", "count"),
            raised_amount_usd_total=("raised_amount_usd", "sum"),
            raised_amount_gbp_total=("raised_amount_gbp", "sum"),
        )
        .reset_index()
        .query(f"year>={min_year}")
        .query(f"year<={max_year}")
        .pipe(impute_empty_years, min_year=min_year, max_year=max_year)
    )

def cb_get_all_timeseries(
    cb_orgs: pd.DataFrame,
    cb_funding_rounds: pd.DataFrame,
    min_year: int = None,
    max_year: int = None,
) -> pd.DataFrame:
    """
    Calculates all typical time series from a list of GtR projects and return
    as one combined table

    Args:
        cb_orgs: A dataframe with columns for 'id' and 'founded_on' among other data
        cb_funding_rounds: A dataframe with columns for 'funding_round_id', 'raised_amount_usd'
            'raised_amount_gbp' and 'announced_on' among other data
        min_year: Lower bound for years to keep
        max_year: Higher bound for years to keep

    Returns:
        A dataframe with the following columns:
            'year',
            'no_of_rounds' - number of funding rounds (deals) in a given year
            'raised_amount_usd_total' - total raised investment (USD) in a given year
            'raised_amount_gbp_total' - total raised investment (GBP) in a given year
            'no_of_orgs_founded' - number of new organisations founded in a given year
    """
    # Number of new companies per year
    time_series_orgs_founded = cb_orgs_founded_per_year(cb_orgs, min_year, max_year)
    # Amount of raised investment per year
    time_series_investment = cb_investments_per_year(
        cb_funding_rounds, min_year, max_year
    )
    # Join up both tables
    time_series_investment["no_of_orgs_founded"] = time_series_orgs_founded[
        "no_of_orgs_founded"
    ]
    return time_series_investment

### Time series trends

def moving_average(
    timeseries_df: pd.DataFrame, window: int = 3, replace_columns: bool = False
) -> pd.DataFrame:
    """
    Calculates rolling mean of yearly timeseries (not centered)

    Args:
        timeseries_df: Should have a 'year' column and at least one other data column
        window: Window of the rolling mean
        rename_cols: If True, will create new set of columns for the moving average
            values with the name pattern `{column_name}_sma{window}` where sma
            stands for 'simple moving average'; otherwise this will replace the original columns

    Returns:
        Dataframe with moving average values
    """
    # Rolling mean
    df_ma = timeseries_df.rolling(window, min_periods=1).mean().drop("year", axis=1)
    # Create new renamed columns
    if not replace_columns:
        column_names = timeseries_df.drop("year", axis=1).columns
        new_column_names = ["{}_sma{}".format(s, window) for s in column_names]
        df_ma = df_ma.rename(columns=dict(zip(column_names, new_column_names)))
        return pd.concat([timeseries_df, df_ma], axis=1)
    else:
        return pd.concat([timeseries_df[["year"]], df_ma], axis=1)

def magnitude(time_series: pd.DataFrame, year_start: int, year_end: int) -> pd.Series:
    """
    Calculates signals' magnitude (i.e. mean across year_start and year_end)

    Args:
        time_series: A dataframe with a columns for 'year' and other data
        year_start: First year of the trend window
        year_end: Last year of the trend window

    Returns:
        Series with magnitude estimates for all data columns
    """
    magnitude = time_series.set_index("year").loc[year_start:year_end, :].mean()
    return magnitude

def percentage_change(initial_value, new_value):
    """Calculates percentage change from first_value to second_value"""
    return (new_value - initial_value) / initial_value * 100

def smoothed_growth(
    time_series: pd.DataFrame, year_start: int, year_end: int, window: int = 3
) -> pd.Series:
    """Calculates a growth estimate by using smoothed (rolling mean) time series

    Args:
        time_series: A dataframe with a columns for 'year' and other data
        year_start: First year of the trend window
        year_end: Last year of the trend window
        window: Moving average windows size (in years) for the smoothed growth estimate

    Returns:
        Series with smoothed growth estimates for all data columns
    """
    # Smooth timeseries
    ma_df = moving_average(time_series, window, replace_columns=True).set_index("year")
    # Percentage change
    return percentage_change(
        initial_value=ma_df.loc[year_start, :], new_value=ma_df.loc[year_end, :]
    )

def estimate_magnitude_growth(
    time_series: pd.DataFrame, year_start: int, year_end: int, window: int = 3
) -> pd.DataFrame:
    """
    Calculates signals' magnitude, estimates their growth and returns a combined dataframe

    Args:
        time_series: A dataframe with a columns for 'year' and other data
        year_start: First year of the trend window
        year_end: Last year of the trend window
        window: Moving average windows size (in years) for the smoothed growth estimate

    Returns:
        Dataframe with magnitude and growth trend estimates; magnitude is in
        absolute units (e.g. GBP 1000s if analysing research funding) whereas
        growth is expresed as a percentage
    """
    magnitude_df = magnitude(time_series, year_start, year_end)
    growth_df = smoothed_growth(time_series, year_start, year_end, window)
    combined_df = (
        pd.DataFrame([magnitude_df, growth_df], index=["magnitude", "growth"])
        .reset_index()
        .rename(columns={"index": "trend"})
    )
    return combined_df


text:
from collections import OrderedDict, Sized
from scipy.stats import t as scipy_stats_t

import inspect
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import scipy.signal as signal

from .fit import fit

class fit2d(fit):
    def __init__(self, f, yx_data, z_data, scatter_kwargs={}, plot_kwargs={}, imshow_kwargs={},
                 contour_kwargs={}, plotting_function=None, fitting_function=None,
                 x_rescale=1, y_rescale=1, **kwargs):

        self.f = f
        self.yx_data = yx_data
        self.z_data = z_data
        self.scatter_kwargs = {**{'zorder':1, 'color':'green', 'alpha':0.5}, **scatter_kwargs}
        self.plot_kwargs = {**{'zorder':3}, **plot_kwargs}
        self.imshow_kwargs = {**{'aspect':y_rescale/x_rescale}, **imshow_kwargs}
        self.contour_kwargs = {**{'colors':'w', 'alpha':0.4}, **contour_kwargs}
        self.x_rescale = x_rescale
        self.y_rescale = y_rescale

        if not plotting_function:
            plotting_function = self.plotting_function_2d

        if not fitting_function:
            fitting_function = self.curve_fit_wrap

        super(fit2d, self).__init__(f, yx_data, np.ravel(z_data), fitting_function=fitting_function,
                                    plotting_function=plotting_function, **kwargs)

    def curve_fit_wrap(self, f, *args, **kwargs):
        def fit_function_wrap(*args, **kwargs):
            return np.ravel(f(*args, **kwargs))
        return scipy.optimize.curve_fit(fit_function_wrap, **kwargs)

    def plotting_function_2d(self, f, val, yx, z):

        z = z.reshape(yx[0].shape)

        left, width = 0.1, 0.65
        bottom, height = 0.1, 0.65
        bottom_h = left_h = left + width + 0.02

        image_ratio = z.shape[0]/z.shape[1]

        rect_imshow = [left, bottom, width, height]
        rect_projx = [left, bottom_h, width, 0.2]
        rect_projy = [left_h, bottom, 0.2, height]
        rect_colorbar = [left_h, bottom_h, 0.1, 0.2]

        if image_ratio>1:
            self.fig = plt.figure(1, figsize=(8*self.x_rescale/image_ratio, 8*self.y_rescale))
        else:
            self.fig = plt.figure(1, figsize=(12*self.x_rescale, 12*image_ratio*self.y_rescale))

        self.ax_imshow = plt.axes(rect_imshow)
        self.ax_projx = plt.axes(rect_projx)
        self.ax_projy = plt.axes(rect_projy)
        self.ax_colorbar = plt.axes(rect_colorbar)

        self.ax_projx.xaxis.set_major_formatter(plt.NullFormatter())
        self.ax_projy.yaxis.set_major_formatter(plt.NullFormatter())

        z_data_imshow = self.ax_imshow.imshow(z, **self.imshow_kwargs)
        self.fit_result_data = f(yx, *val).reshape(z.shape)
        self.ax_imshow.contour(yx[1], yx[0], self.fit_result_data, 10, **self.contour_kwargs)

        self.ax_projy.scatter(np.sum(z, 1), np.arange(z.shape[0]), **self.scatter_kwargs)
        self.ax_projy.plot(np.sum(self.fit_result_data, 1), np.arange(self.fit_result_data.shape[0]), **self.plot_kwargs)

        self.ax_projx.scatter(np.arange(z.shape[1]), np.sum(z, 0),**self.scatter_kwargs)
        self.ax_projx.plot(np.arange(self.fit_result_data.shape[1]), np.sum(self.fit_result_data, 0), **self.plot_kwargs)

        plt.colorbar(z_data_imshow, cax=self.ax_colorbar)

        self.ax_projx.set_xlim(-0.5, z.shape[1]-0.5)
        self.ax_projx.grid(alpha=0.3)
        self.ax_projy.set_ylim(z.shape[0]-0.5, -0.5)
        self.ax_projy.grid(alpha=0.3)
        self.ax_imshow.grid(alpha=0.3)

ext: py

max_stars_repo_name: AndreaAmico/fitwrap

max_issues_repo_name: AndreaAmico/fitwrap

max_forks_repo_name: AndreaAmico/fitwrap

reason: import numpy,import scipy,from scipy


text:
import os.path as osp
import glob
import numpy as np
import json

import neuron.ops as ops
from neuron.config import registry
from .dataset import SeqDataset

__all__ = ['UAV123']

@registry.register_module
class UAV123(SeqDataset):
    """`UAV123 <https://ivul.kaust.edu.sa/Pages/pub-benchmark-simulator-uav.aspx>`_ Dataset.

    Publication:
        ``A Benchmark and Simulator for UAV Tracking``,
        M. Mueller, N. Smith and B. Ghanem, ECCV 2016.

    Args:
        root_dir (string): Root directory of dataset where sequence
            folders exist.
        version (integer or string): Specify the benchmark version, specify as one of
            ``UAV123`` and ``UAV20L``.
    """
    def __init__(self, root_dir=None, version='UAV123'):
        assert version.upper() in ['UAV20L', 'UAV123']
        if root_dir is None:
            root_dir = osp.expanduser('~/data/UAV123')
        self.root_dir = root_dir
        self.version = version.upper()

        # initialize the dataset
        super(UAV123, self).__init__(
            name=self.version,
            root_dir=self.root_dir,
            version=self.version)

    def _construct_seq_dict(self, root_dir, version):
        # sequence meta information
        meta_file = osp.join(
            osp.dirname(__file__), 'uav123.json')
        with open(meta_file) as f:
            seq_metas = json.load(f)

        # image and annotation paths
        anno_files = sorted(glob.glob(
            osp.join(root_dir, 'anno/%s/*.txt' % version)))
        seq_names = [
            osp.basename(f)[:-4] for f in anno_files]
        seq_dirs = [osp.join(
            root_dir, 'data_seq/UAV123/%s' % \
                seq_metas[version][n]['folder_name'])
            for n in seq_names]

        # construct seq_dict
        seq_dict = {}
        for s, seq_name in enumerate(seq_names):
            # valid frame range
            start_frame = seq_metas[version][
                seq_names[s]]['start_frame']
            end_frame = seq_metas[version][
                seq_names[s]]['end_frame']
            img_files = [osp.join(
                seq_dirs[s], '%06d.jpg' % f)
                for f in range(start_frame, end_frame + 1)]

            # load annotations
            anno = np.loadtxt(anno_files[s], delimiter=',')
            anno[:, 2:] = anno[:, :2] + anno[:, 2:] - 1

            # meta information
            seq_len = len(img_files)
            img0 = ops.read_image(img_files[0])
            meta = {
                'width': img0.shape[1],
                'height': img0.shape[0],
                'frame_num': seq_len,
                'target_num': 1,
                'total_instances': seq_len}

            # update seq_dict
            seq_dict[seq_name] = {
                'img_files': img_files,
                'target': {
                    'anno': anno,
                    'meta': meta}}

        return seq_dict

max_stars_repo_name: tsingqguo/ABA

max_issues_repo_name: tsingqguo/ABA

max_forks_repo_name: tsingqguo/ABA

reason: import numpy


text:
# noinspection PyInterpreter
import logging
import pickle
import time

import cloudpickle
import numpy as np
from mpi4py import MPI

import abcpy.backends.mpimanager
from abcpy.backends import BDS, PDS, Backend, NestedParallelizationController

class NestedParallelizationControllerMPI(NestedParallelizationController):
    def __init__(self, mpi_comm):
        self.logger = logging.getLogger(__name__)
        self.logger.info("#### Initialize NPC ####")
        self.loop_workers = True
        self.mpi_comm = mpi_comm
        self.nested_func = "NoFunction"
        self.func_args = ()
        self.func_kwargs = {}
        self.result = None
        if self.mpi_comm.Get_rank() != 0:
            self.nested_execution()

    def communicator(self):
        return self.mpi_comm

    def nested_execution(self):
        rank = self.mpi_comm.Get_rank()
        self.logger.debug("Starting nested loop on rank {}".format(rank))
        while self.loop_workers:
            self.mpi_comm.barrier()
            self.loop_workers = self.mpi_comm.bcast(self.loop_workers, root=0)
            if self.loop_workers == False:
                return
            func_p = None
            func_args_p = None
            func_kwargs_p = None
            if self.mpi_comm.Get_rank() == 0:
                self.logger.debug("Start pickling func on rank {}".format(rank))
                func_p = cloudpickle.dumps(self.nested_func, pickle.HIGHEST_PROTOCOL)
                func_args_p = cloudpickle.dumps(self.func_args, pickle.HIGHEST_PROTOCOL)
                func_kwargs_p = cloudpickle.dumps(self.func_kwargs, pickle.HIGHEST_PROTOCOL)

            self.logger.debug("Broadcasting function {} on rank {}".format(self.nested_func, rank))
            func_p = self.mpi_comm.bcast(func_p, root=0)
            func_args_p = self.mpi_comm.bcast(func_args_p, root=0)
            func_kwargs_p = self.mpi_comm.bcast(func_kwargs_p, root=0)
            self.nested_func = cloudpickle.loads(func_p)
            self.func_args = cloudpickle.loads(func_args_p)
            self.func_kwargs = cloudpickle.loads(func_kwargs_p)

            func = self.nested_func
            self.logger.debug("Starting map function {} on rank {}".format(func.__name__, self.mpi_comm.Get_rank()))
            self.func_kwargs['mpi_comm'] = self.mpi_comm
            self.mpi_comm.barrier()
            self.result = func(*self.func_args, **self.func_kwargs)
            self.logger.debug("Ending map function on rank {}".format(self.mpi_comm.Get_rank()))
            self.mpi_comm.barrier()
            if self.mpi_comm.Get_rank() == 0:
                return
        self.loop_workers = True
        self.logger.debug("Ending nested loop on rank {}".format(self.mpi_comm.Get_rank()))

    def run_nested(self, func, *args, **kwargs):
        self.logger.debug("Executing nested function {}.".format(func.__name__))
        self.nested_func = func
        self.func_args = args
        self.func_kwargs = kwargs
        self.nested_execution()
        self.logger.debug("Return from nested execution of master rank")
        self.nested_func = None
        self.func_args = ()
        self.func_kwargs = {}
        self.logger.info(self.result)
        return self.result

    def __del__(self):
        rank = self.mpi_comm.Get_rank()
        self.logger.debug("Stopping npc on rank {}".format(rank))
        self.loop_workers = False
        if rank == 0:
            self.mpi_comm.barrier()
            self.loop_workers = self.mpi_comm.bcast(self.loop_workers, root=0)
        self.logger.debug(">>>>>>>> NPC stopped on rank {}".format(rank))

class BackendMPIScheduler(Backend):
    """Defines the behavior of the scheduler process

    This class defines the behavior of the scheduler process (The one
    with rank==0) in MPI.

    """

    # Define some operation codes to make it more readable
    OP_PARALLELIZE, OP_MAP, OP_COLLECT, OP_BROADCAST, OP_DELETEPDS, OP_DELETEBDS, OP_FINISH = [1, 2, 3, 4, 5, 6, 7]
    finalized = False

    def __init__(self, chunk_size=1):
        """
        Parameters
        ----------
        chunk_size: Integer
            size of one block of data to be sent to free
            execution teams
       """

        # Initialize the current_pds_id and bds_id
        self.__current_pds_id = 0
        self.__current_bds_id = 0

        # Initialize a BDS store for both scheduler & team.
        self.bds_store = {}
        self.pds_store = {}

        # Initialize a store for the pds data that
        # .. hasn't been sent to the teams yet
        self.pds_pending_store = {}

        self.chunk_size = chunk_size

    def __command_teams(self, command, data):
        """Tell teams to enter relevant execution block
        This method handles the sending of the command to the teams
        telling them what operation to perform next.

        Parameters
        ----------
        command: operation code of OP_xxx
            One of the operation codes defined in the class definition as OP_xxx
            which tell the teams what operation they're performing.
        data:  tuple
            Any of the data required for the operation which needs to be bundled
            in the data packet sent.
        """

        if command == self.OP_PARALLELIZE:
            # In parallelize we receive data as (pds_id)
            data_packet = (command, data[0])

        elif command == self.OP_MAP:
            # In map we receive data as (pds_id,pds_id_new,func)
            # Use cloudpickle to dump the function into a string.
            function_packed = cloudpickle.dumps(data[2], pickle.HIGHEST_PROTOCOL)
            data_packet = (command, data[0], data[1], function_packed)

        elif command == self.OP_BROADCAST:
            data_packet = (command, data[0])

        elif command == self.OP_COLLECT:
            # In collect we receive data as (pds_id)
            data_packet = (command, data[0])

        elif command == self.OP_DELETEPDS or command == self.OP_DELETEBDS:
            # In deletepds we receive data as (pds_id) or bds_id
            data_packet = (command, data[0])

        elif command == self.OP_FINISH:
            data_packet = (command,)

        _ = self.mpimanager.get_scheduler_communicator().bcast(data_packet, root=0)

    def __generate_new_pds_id(self):
        """
        This method generates a new pds_id to associate a PDS with it's remote counterpart
        that teams use to store & index data based on the pds_id they receive

        Returns
        -------
        Returns a unique integer id.

        """

        self.__current_pds_id += 1
        return self.__current_pds_id

    def __generate_new_bds_id(self):
        """
        This method generates a new bds_id to associate a BDS with it's remote counterpart
        that teams use to store & index data based on the bds_id they receive

        Returns
        -------
        Returns a unique integer id.

        """

        self.__current_bds_id += 1
        return self.__current_bds_id

    def parallelize(self, python_list):
        """
        This method distributes the list on the available teams and returns a
        reference object.

        The list is split into number of teams many parts as a numpy array.
        Each part is sent to a separate team node using the MPI scatter.

        scheduler: python_list is the real data that is to be split up

        Parameters
        ----------
        list: Python list
            the list that should get distributed on the leader nodes of the teams

        Returns
        -------
        PDSMPI class (parallel data set)
            A reference object that represents the parallelized list
        """

        # Tell the teams to enter parallelize()
        pds_id = self.__generate_new_pds_id()
        self.__command_teams(self.OP_PARALLELIZE, (pds_id,))

        # Don't send any data. Just keep it as a queue we're going to pop.
        self.pds_store[pds_id] = list(python_list)

        pds = PDSMPI([], pds_id, self)

        return pds

    def orchestrate_map(self, pds_id):
        """Orchestrates the teams to perform a map function

        This works by keeping track of the teams who haven't finished executing,
        waiting for them to request the next chunk of data when they are free,
        responding to them with the data and then sending them a Sentinel
        signalling that they can exit.
        """
        is_map_done = [True if i in self.mpimanager.get_scheduler_node_ranks() else False for i in
                       range(self.mpimanager.get_scheduler_size())]
        status = MPI.Status()

        # Copy it to the pending. This is so when scheduler accesses
        # the PDS data it's not empty.
        self.pds_pending_store[pds_id] = list(self.pds_store[pds_id])

        # While we have some ranks that haven't finished
        while sum(is_map_done) < self.mpimanager.get_scheduler_size():
            # Wait for a request from anyone
            data_request = self.mpimanager.get_scheduler_communicator().recv(
                source=MPI.ANY_SOURCE,
                tag=MPI.ANY_TAG,
                status=status,
            )
            request_from_rank = status.source

            if data_request != pds_id:
                print("Ignoring stale PDS data request from",
                      request_from_rank, ":", data_request, "/", pds_id)
                continue

            # Pointer so we don't have to keep doing dict lookups
            current_pds_items = self.pds_pending_store[pds_id]
            num_current_pds_items = len(current_pds_items)

            # Everyone's already exhausted all the data.
            # Send a sentinel and mark the node as finished
            if num_current_pds_items == 0:
                self.mpimanager.get_scheduler_communicator().send(None, dest=request_from_rank, tag=pds_id)
                is_map_done[request_from_rank] = True
            else:
                # Create the chunk of data to send. Pop off items and tag them with an id.
                # so we can sort them later
                chunk_to_send = []
                for i in range(self.chunk_size):
                    chunk_to_send += [(num_current_pds_items - i, current_pds_items.pop())]
                    self.mpimanager.get_scheduler_communicator().send(chunk_to_send, dest=request_from_rank, tag=pds_id)

    def map(self, func, pds):
        """
        A distributed implementation of map that works on parallel data sets (PDS).

        On every element of pds the function func is called.

        Parameters
        ----------
        func: Python func
            A function that can be applied to every element of the pds
        pds: PDS class
            A parallel data set to which func should be applied

        Returns
        -------
        PDSMPI class
            a new parallel data set that contains the result of the map
        """

        # Tell the teams to enter the map() with the current pds_id & func.
        # Get pds_id of dataset we want to operate on
        pds_id = pds.pds_id

        # Generate a new pds_id to be used by the teams for the resultant PDS
        pds_id_new = self.__generate_new_pds_id()

        data = (pds_id, pds_id_new, func)
        self.__command_teams(self.OP_MAP, data)
        self.orchestrate_map(pds_id)

        pds_res = PDSMPI([], pds_id_new, self)

        return pds_res

    def collect(self, pds):
        """
        Gather the pds from all the teams,
            send it to the scheduler and return it as a standard Python list.

        Parameters
        ----------
        pds: PDS class
            a parallel data set

        Returns
        -------
        Python list
            all elements of pds as a list
        """

        # Tell the teams to enter collect with the pds's pds_id
        self.__command_teams(self.OP_COLLECT, (pds.pds_id,))

        # all_data = self.world_communicator.gather(pds.python_list, root=0)
        all_data = self.mpimanager.get_scheduler_communicator().gather(pds.python_list, root=0)

        # Initialize lists to accumulate results
        all_data_indices, all_data_items = [], []

        for node_data in all_data:
            for index, item in node_data:
                if isinstance(item, Exception):
                    raise item
                all_data_indices.append(index)
                all_data_items.append(item)

        # Sort the accumulated data according to the indices we tagged
        # them with when distributing
        rdd_sorted = [all_data_items[i] for i in np.argsort(all_data_indices)]
        return rdd_sorted

    def broadcast(self, value):
        """
        Sends a data to all leaders and workers
        First instruction is sent to leaders which then send it to their workers
        Then every process enters a broadcast to receive data from scheduler
        """
        # Tell the teams to enter broadcast()
        bds_id = self.__generate_new_bds_id()
        self.__command_teams(self.OP_BROADCAST, (bds_id,))

        _ = self.mpimanager.get_world_communicator().bcast(value, root=0)

        bds = BDSMPI(value, bds_id, self)
        return bds

    def delete_remote_pds(self, pds_id):
        """
        A public function for the PDS objects on the scheduler to call when they go out of
        scope or are deleted in order to ensure the same happens on the teams.

        Parameters
        ----------
        pds_id: int
            A pds_id identifying the remote PDS on the teams to delete.
        """

        if not self.finalized:
            self.__command_teams(self.OP_DELETEPDS, (pds_id,))

    def delete_remote_bds(self, bds_id):
        """
        Public function for the BDS objects on the scheduler to call when they go
        out of score or are deleted in order to ensure they are deleted
        ont he teams as well.

        Parameters
        ----------
        bds_id: int
            A bds_id identifying the remote BDS on the teams to delete.
        """

        if not self.finalized:
            # The scheduler deallocates it's BDS data. Explicit because
            # .. bds_store and BDSMPI object are disconnected.
            del backend.bds_store[bds_id]
            self.__command_teams(self.OP_DELETEBDS, (bds_id,))

    def __del__(self):
        """
        Overriding the delete function to explicitly call MPI.finalize().
        This is also required so we can tell the teams to get out of the
        while loop they are in and exit gracefully and they themselves call
        finalize when they die.
        """

        # Tell the teams they can exit gracefully.
        self.__command_teams(self.OP_FINISH, None)

        # Finalize the connection because the teams should have finished.
        MPI.Finalize()
        self.finalized = True

class BackendMPIWorker(Backend):
    """
    Workers are processes that are used to execute (maybe MPI) models
    There is one communicator per model to execute, compounded of one leader and workers
    Leaders receives instructions from the scheduler which then transmit them to workers
    Leaders are themselves workers
    """

    OP_PARALLELIZE, OP_MAP, OP_COLLECT, OP_BROADCAST, OP_DELETEPDS, OP_DELETEBDS, OP_FINISH = [1, 2, 3, 4, 5, 6, 7]

    def __init__(self):
        """ No parameter, just call worker_run """
        self.logger = logging.getLogger(__name__)
        self.__worker_run()

    def run_function(self, function_packed, data_item):
        """
        Receives a serialized function unpack it and run it
        Passes the model communicator if ther is more than one process per model
        """
        func = cloudpickle.loads(function_packed)
        res = None
        try:
            if self.mpimanager.get_model_size() > 1:
                npc = NestedParallelizationControllerMPI(self.mpimanager.get_model_communicator())
                if self.mpimanager.get_model_communicator().Get_rank() == 0:
                    self.logger.debug("Executing map function on master rank 0.")
                    res = func(data_item, npc=npc)
                del npc
            else:
                res = func(data_item)
        except Exception as e:
            msg = "Exception occured while calling the map function {}: ".format(func.__name__)
            res = type(e)(msg + str(e))
        return res

    def __worker_run(self):
        """
        Workers enter an infinite loop and waits for instructions from their leader
        """
        while True:
            data = self.mpimanager.get_model_communicator().bcast(None, root=0)
            op = data[0]
            if op == self.OP_MAP:
                # Receive data from scheduler of the model
                function_packed = self.mpimanager.get_model_communicator().bcast(None, root=0)[0]
                data_item = self.mpimanager.get_model_communicator().bcast(None, root=0)[0]
                self.run_function(function_packed, data_item)
            elif op == self.OP_BROADCAST:
                self._bds_id = data[1]
                self.broadcast(None)
            elif op == self.OP_FINISH:
                quit()
            else:
                raise Exception("worker model received unknown command code")

    def collect(self):
        pass

    def map(self):
        pass

    def parallelize():
        pass

    def broadcast(self, value):
        """
        Receives data from scheduler
        """
        value = self.mpimanager.get_world_communicator().bcast(None, root=0)
        self.bds_store[self._bds_id] = value

class BackendMPILeader(BackendMPIWorker):
    """Defines the behavior of the leader processes

    This class defines how the leaders should behave during operation.
    leaders are those processes(not nodes like Spark) that have rank==0 in the model communicator
    """

    OP_PARALLELIZE, OP_MAP, OP_COLLECT, OP_BROADCAST, OP_DELETEPDS, OP_DELETEBDS, OP_FINISH = [1, 2, 3, 4, 5, 6, 7]

    def __init__(self):
        """ No parameter, just call leader_run """
        self.logger = logging.getLogger(__name__)
        self.__leader_run()

    def __leader_run(self):
        """
        This method is the infinite loop a leader enters directly from init.
        It makes the leader wait for a command to perform from the scheduler and
        then calls the appropriate function.

        This method also takes care of the synchronization of data between the
        scheduler and the leaders by matching PDSs based on the pds_ids sent by the scheduler
        with the command.

        Commands received from the scheduler are of the form of a tuple.
        The first component of the tuple is always the operation to be performed
        and the rest are conditional on the operation.

        (op,pds_id) where op == OP_PARALLELIZE for parallelize
        (op,pds_id, pds_id_result,func) where op == OP_MAP for map.
        (op,pds_id) where op == OP_COLLECT for a collect operation
        (op,pds_id) where op == OP_DELETEPDS for a delete of the remote PDS on leaders
        (op,) where op==OP_FINISH for the leader to break out of the loop and terminate
        """

        # Initialize PDS data store here because only teams need to do it.
        self.pds_store = {}

        while True:
            data = self.mpimanager.get_scheduler_communicator().bcast(None, root=0)

            op = data[0]
            if op == self.OP_PARALLELIZE:
                pds_id = data[1]
                self._rec_pds_id = pds_id
                pds_id, pds_id_new = self.__get_received_pds_id()
                self.pds_store[pds_id] = None

            elif op == self.OP_MAP:
                pds_id, pds_id_result, function_packed = data[1:]
                self._rec_pds_id, self._rec_pds_id_result = pds_id, pds_id_result

                # Enter the map so we can grab data and perform the func.
                # Func sent before and not during for performance reasons
                pds_res = self.map(function_packed)

                # Store the result in a newly generated PDS pds_id
                self.pds_store[pds_res.pds_id] = pds_res

            elif op == self.OP_BROADCAST:
                self._bds_id = data[1]
                # relay command and data into model communicator
                self.mpimanager.get_model_communicator().bcast(data, root=0)
                self.broadcast(None)

            elif op == self.OP_COLLECT:
                pds_id = data[1]

                # Access an existing PDS from data store
                pds = self.pds_store[pds_id]

                self.collect(pds)

            elif op == self.OP_DELETEPDS:
                pds_id = data[1]
                del self.pds_store[pds_id]

            elif op == self.OP_DELETEBDS:
                bds_id = data[1]
                del self.bds_store[bds_id]

            elif op == self.OP_FINISH:
                # tells other processes of the worker to finish
                self.mpimanager.get_model_communicator().bcast([self.OP_FINISH], root=0)
                quit()
            else:
                raise Exception("team received unknown command code")

    def __get_received_pds_id(self):
        """
        Function to retrieve the pds_id(s) we received from the scheduler to associate
        our team's created PDS with the scheduler's.
        """

        return self._rec_pds_id, self._rec_pds_id_result

    def __leader_run_function(self, function_packed, data_item):
        """
        This function sends data and serialized function to workers and executes it
        """
        self.mpimanager.get_model_communicator().bcast([self.OP_MAP], root=0)
        self.mpimanager.get_model_communicator().bcast([function_packed], root=0)
        self.mpimanager.get_model_communicator().bcast([data_item], root=0)
        return self.run_function(function_packed, data_item)

    def parallelize(self):
        pass

    def map(self, function_packed):
        """
        A distributed implementation of map that works on parallel data sets (PDS).
        On every element of pds the function func is called.
        We consider that process 0 of each MPI model should return the final result.

        Parameters
        ----------
        func: Python function_packed
            A serialized function that can be applied to every element of the pds

        Returns
        -------
        PDSMPI class
            a new parallel data set that contains the result of the map
        """

        map_start = time.time()

        # Get the PDS id we operate on and the new one to store the result in
        pds_id, pds_id_new = self.__get_received_pds_id()

        rdd = []
        while True:
            # Ask for a chunk of data since it's free
            data_chunks = self.mpimanager.get_scheduler_communicator().sendrecv(pds_id, 0, pds_id)

            # If it receives a sentinel, it's done and it can exit
            if data_chunks is None:
                break

            # Accumulate the indices and *processed* chunks
            for chunk in data_chunks:
                data_index, data_item = chunk
                res = self.__leader_run_function(function_packed, data_item)
                rdd += [(data_index, res)]

        pds_res = PDSMPI(rdd, pds_id_new, self)

        return pds_res

    def collect(self, pds):
        """
        Gather the pds from all the leaders,
        send it to the scheduler and return it as a standard Python list.

        Parameters
        ----------
        pds: PDS class
            a parallel data set

        Returns
        -------
        Python list
            all elements of pds as a list
        """

        # Send the data we have back to the scheduler
        _ = self.mpimanager.get_scheduler_communicator().gather(pds.python_list, root=0)

class BackendMPITeam(BackendMPILeader if abcpy.backends.mpimanager.get_mpi_manager().is_leader() else BackendMPIWorker):
    """
    A team is compounded by workers and a leader. One process per team is a leader, others are workers
    """

    OP_PARALLELIZE, OP_MAP, OP_COLLECT, OP_BROADCAST, OP_DELETEPDS, OP_DELETEBDS, OP_FINISH = [1, 2, 3, 4, 5, 6, 7]

    def __init__(self):
        # Define the vars that will hold the pds ids received from scheduler to operate on
        self._rec_pds_id = None
        self._rec_pds_id_result = None

        # Initialize a BDS store for both scheduler & team.
        self.bds_store = {}

        # print("In BackendMPITeam, rank : ", self.rank, ", model_rank_global : ", globals()['model_rank_global'])

        self.logger = logging.getLogger(__name__)
        super().__init__()

class BackendMPI(BackendMPIScheduler if abcpy.backends.mpimanager.get_mpi_manager().is_scheduler() else BackendMPITeam):
    """A backend parallelized by using MPI

    The backend conditionally inherits either the BackendMPIScheduler class
    or the BackendMPIteam class depending on it's rank. This lets
    BackendMPI have a uniform interface for the user but allows for a
    logical split between functions performed by the scheduler
    and the teams.
    """

    def __init__(self, scheduler_node_ranks=[0], process_per_model=1):
        """
        Parameters
        ----------
        scheduler_node_ranks: Python list
            list of scheduler nodes

        process_per_model: Integer
            number of MPI processes to allocate to each model
        """
        # get mpimanager instance from the mpimanager module (which has to be setup before calling the constructor)
        self.logger = logging.getLogger(__name__)
        self.mpimanager = abcpy.backends.mpimanager.get_mpi_manager()

        if self.mpimanager.get_world_size() < 2:
            raise ValueError('A minimum of 2 ranks are required for the MPI backend')

        # Set the global backend
        globals()['backend'] = self

        # Call the appropriate constructors and pass the required data
        super().__init__()

    def size(self):
        """ Returns world size """
        return self.mpimanager.get_world_size()

    def scheduler_node_ranks(self):
        """ Returns scheduler node ranks """
        return self.mpimanager.get_scheduler_node_ranks()

    @staticmethod
    def disable_nested(mpi_comm):
        if mpi_comm.Get_rank() != 0:
            mpi_comm.Barrier()

    @staticmethod
    def enable_nested(mpi_comm):
        if mpi_comm.Get_rank() == 0:
            mpi_comm.Barrier()

class PDSMPI(PDS):
    """
    This is an MPI wrapper for a Python parallel data set.
    """

    def __init__(self, python_list, pds_id, backend_obj):
        self.python_list = python_list
        self.pds_id = pds_id
        self.backend_obj = backend_obj

    def __del__(self):
        """
        Destructor to be called when a PDS falls out of scope and/or is being deleted.
        Uses the backend to send a message to destroy the teams' copy of the pds.
        """
        try:
            self.backend_obj.delete_remote_pds(self.pds_id)
        except AttributeError:
            # Catch "delete_remote_pds not defined" for teams and ignore.
            pass

class BDSMPI(BDS):
    """
    This is a wrapper for MPI's BDS class.
    """

    def __init__(self, object, bds_id, backend_obj):
        # The BDS data is no longer saved in the BDS object.
        # It will access & store the data only from the current backend
        self.bds_id = bds_id
        backend.bds_store[self.bds_id] = object

    def value(self):
        """
        This method returns the actual object that the broadcast data set represents.
        """
        return backend.bds_store[self.bds_id]

    def __del__(self):
        """
        Destructor to be called when a BDS falls out of scope and/or is being deleted.
        Uses the backend to send a message to destroy the teams' copy of the bds.
        """

        try:
            backend.delete_remote_bds(self.bds_id)
        except AttributeError:
            # Catch "delete_remote_pds not defined" for teams and ignore.
            pass

class BackendMPITestHelper:
    """
    Helper function for some of the test cases to be able to access and verify class members.
    """

    def check_pds(self, k):
        """Checks if a PDS exists in the pds data store. Used to verify deletion and creation
        """
        return k in backend.pds_store.keys()

    def check_bds(self, k):
        """Checks if a BDS exists in the BDS data store. Used to verify deletion and creation
        """
        return k in backend.bds_store.keys()

max_stars_repo_name: vishalbelsare/abcpy

max_issues_repo_name: vishalbelsare/abcpy

max_forks_repo_name: vishalbelsare/abcpy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 16 13:46:54 2015

@author: egendron
"""

import pyfits
import numpy as np
import gral_createModalBasis

import matplotlib.pyplot as plt

def checkSMAandCo(path, filename):
    """
    Input args:
    <path> directory where files are
    <filename> filename of the slopes+voltage data set

    """

    # reading data
    a = pyfits.getdata(path+'/'+filename)
    vHO = a.field(5).T
    vTT = a.field(6).T
    slopes = a.field(4).T

    # reading other data
    AWFbasis = pyfits.getdata(path+'/'+'HOCtr.AWF_IM_KERNEL.fits')
    SMAbasis = pyfits.getdata(path+'/'+'HOCtr.SMA_BASIS.fits')
    HO2TT = pyfits.getdata(path+'/'+'HOCtr.HO_TO_TT.fits')
    TT2HO = pyfits.getdata(path+'/'+'HOCtr.TT_TO_HO.fits')
    pistonMode = pyfits.getdata(path+'/'+'HOCtr.PRA_PISTON_MODE.fits')
    pistonProj = pyfits.getdata(path+'/'+'HOCtr.PRA_PISTON_PROJECTION.fits')
    S2M = pyfits.getdata(path+'/'+'S2M.fits')
    M2V = pyfits.getdata(path+'/'+'M2V.fits')

    # projecting voltages onto piston mode
    pistonProj /= np.dot(pistonProj, pistonMode)  # normalizing things (just as sparta does internally)
    piston = np.dot(pistonProj, vHO) * np.sqrt(np.sum(pistonMode**2))  # projecting onto piston mode
    print 'This is PISTON component : %f ± rms %f ' % (np.average(piston) , np.std(piston))

    # transforming TT signals into HO .. not sure this is useful
    # but on the present data set, one cannot check
    vHO += np.dot(TT2HO, vTT)

    # computing voltages increments between 2 frames
    dvHO = vHO[:,1:] - vHO[:,:-1]
    dvTT = vTT[:,1:] - vTT[:,:-1]  # useless
    slopes = slopes[:,1:]          # drop first frame

    # cutting data set to avoid the beginning, which is openloop
    Nstart = 0
    dvHO = dvHO[:,Nstart:]
    dvTT = dvTT[:,Nstart:]
    slopes = slopes[:,Nstart:]

    # Here i'm trying to re-generate the command increments that
    # occured during the closed loop, and see how well they match the
    # command increments computed using the SPARTA voltages.
    #
    MC = gral_createModalBasis.computeSystemControlMatrix( S2M, M2V, TT2HO*0 )
    sm = np.dot(MC, slopes)

    # now i'm doing a linear regression, on all actuators, to
    # find what was the close-loop gain that was applied
    na     = 60
    gain   = np.zeros(na)
    offset = np.zeros(na)
    for i in range(na):
        gain[i], offset[i] = np.polyfit(sm[i,:], -dvHO[i,:] ,1)
    avgain = np.average(gain)
    print 'Average gain found is %f ± rms %f' % (avgain, np.std(gain))
    print 'Offset variation is rms %f' % np.std(offset)

    # Subtraction between SPARTA incr-voltage conmmands, and
    # our calculation of what should be the incr-volt commands
    errcomm = sm[0:na,:] * avgain + dvHO
    print 'Error between sparta volts and computed volts : %f' % np.average(np.std(errcomm,axis=1))

    # Using the M2V, we're gonna re-generate the projection matrix onto
    # the non-commanded modes
    na = 60
    nuseful = M2V.shape[1]
    M2VHO = np.zeros((na, nuseful))
    M2VHO[:,0:2] = TT2HO
    M2VHO[:,2: ] = M2V[:na,2: ]
    projAWF_parall = np.dot(M2VHO, gral_createModalBasis.pseudoinv( M2VHO, 0 ))
    projAWF_kernel = np.identity(na) - projAWF_parall
    # using sparta 'V2 = AWF_IM_KERNEL' to check we're retrieving the
    # same projection matrix. We use the <AWFbasis> passed to sparta, and use the
    # internal sparta algorithm (the famous V2.transpose(V2)) for generating it.
    # If the <AWFbasis> is correct, both computation should give the same
    # result.
    projAWF = np.dot(AWFbasis, AWFbasis.T)
    # computing difference between sparta proj matrix and our proj matrix
    errAWF = np.sum((projAWF-projAWF_kernel)**2) / np.sum((projAWF)**2)
    print 'Error on AWF matrix is %f' % errAWF

    # now computing the fraction of the command which
    # belongs to the commanded space
    dpHO = np.dot(projAWF_parall, dvHO)
    # and the command lying out of the commanded space (belong to kernel)
    dkHO = np.dot(projAWF_kernel, dvHO)
    # see where kernel-commands were sent
    kpow = np.sqrt(np.sum(dkHO**2, axis=0)) > 5e-7
    plt.plot( kpow )
    print 'Nb of occurences SMA sent kernel-commands: %d' % np.sum(kpow)

max_stars_repo_name: soylentdeen/CIAO-commissioning-tools

max_issues_repo_name: soylentdeen/CIAO-commissioning-tools

max_forks_repo_name: soylentdeen/CIAO-commissioning-tools


text:
import cupy
import cupyx

from cupy import core
from cupy.core import _routines_math as _math
from cupy.core import fusion
from cupy.lib import stride_tricks

_dot_kernel = core.ReductionKernel(
    'T x1, T x2',
    'T y',
    'x1 * x2',
    'a + b',
    'y = a',
    '0',
    'dot_product'
)

def convolve(a, v, mode='full'):
    """Returns the discrete, linear convolution of two one-dimensional sequences.

    Args:
        a (cupy.ndarray): first 1-dimensional input.
        v (cupy.ndarray): second 1-dimensional input.
        mode (str, optional): `valid`, `same`, `full`

    Returns:
        cupy.ndarray: Discrete, linear convolution of a and v.

    .. seealso:: :func:`numpy.convolve`

    """
    if a.size == 0:
        raise ValueError('a cannot be empty')
    if v.size == 0:
        raise ValueError('v cannot be empty')
    if v.ndim > 1:
        raise ValueError('v cannot be multidimensional array')
    if v.size > a.size:
        a, v = v, a
    a = a.ravel()
    v = v.ravel()

    method = cupyx.scipy.signal.choose_conv_method(a, v, mode)
    if method == 'direct':
        out = _dot_convolve(a, v, mode)
    elif method == 'fft':
        out = _fft_convolve(a, v, mode)
    else:
        raise ValueError('Unsupported method')
    return out

def _fft_convolve(a1, a2, mode):

    offset = 0
    if a1.size < a2.size:
        a1, a2 = a2, a1
        offset = 1 - a2.size % 2

    if a1.dtype.kind == 'c' or a2.dtype.kind == 'c':
        fft, ifft = cupy.fft.fft, cupy.fft.ifft
    else:
        fft, ifft = cupy.fft.rfft, cupy.fft.irfft

    dtype = cupy.result_type(a1, a2)
    n1, n2 = a1.size, a2.size
    out_size = n1 + n2 - 1
    fa1 = fft(a1, out_size)
    fa2 = fft(a2, out_size)
    out = ifft(fa1 * fa2, out_size)

    if mode == 'full':
        start, end = None, None
    elif mode == 'same':
        start = (n2 - 1) // 2 + offset
        end = start + n1
    elif mode == 'valid':
        start, end = n2 - 1, n1
    else:
        raise ValueError(
            'acceptable mode flags are `valid`, `same`, or `full`.')

    out = out[start: end]

    if dtype.kind in 'iu':
        out = cupy.around(out)

    return out.astype(dtype, copy=False)

def _dot_convolve(a1, a2, mode):

    offset = 0
    if a1.size < a2.size:
        a1, a2 = a2, a1
        offset = 1 - a2.size % 2

    dtype = cupy.result_type(a1, a2)
    n1, n2 = a1.size, a2.size
    a1 = a1.astype(dtype, copy=False)
    a2 = a2.astype(dtype, copy=False)

    if mode == 'full':
        out_size = n1 + n2 - 1
        a1 = cupy.pad(a1, n2 - 1)
    elif mode == 'same':
        out_size = n1
        pad_size = (n2 - 1) // 2 + offset
        a1 = cupy.pad(a1, (n2 - 1 - pad_size, pad_size))
    elif mode == 'valid':
        out_size = n1 - n2 + 1

    stride = a1.strides[0]
    a1 = stride_tricks.as_strided(a1, (out_size, n2), (stride, stride))
    output = _dot_kernel(a1, a2[::-1], axis=1)
    return output

def clip(a, a_min=None, a_max=None, out=None):
    """Clips the values of an array to a given interval.

    This is equivalent to ``maximum(minimum(a, a_max), a_min)``, while this
    function is more efficient.

    Args:
        a (cupy.ndarray): The source array.
        a_min (scalar, cupy.ndarray or None): The left side of the interval.
            When it is ``None``, it is ignored.
        a_max (scalar, cupy.ndarray or None): The right side of the interval.
            When it is ``None``, it is ignored.
        out (cupy.ndarray): Output array.

    Returns:
        cupy.ndarray: Clipped array.

    .. seealso:: :func:`numpy.clip`

    """
    if fusion._is_fusing():
        return fusion._call_ufunc(_math.clip,
                                  a, a_min, a_max, out=out)

    # TODO(okuta): check type
    return a.clip(a_min, a_max, out=out)

# sqrt_fixed is deprecated.
# numpy.sqrt is fixed in numpy 1.11.2.
sqrt = sqrt_fixed = core.sqrt

cbrt = core.create_ufunc(
    'cupy_cbrt',
    ('e->e', 'f->f', 'd->d'),
    'out0 = cbrt(in0)',
    doc='''Elementwise cube root function.

    .. seealso:: :data:`numpy.cbrt`

    ''')

square = core.create_ufunc(
    'cupy_square',
    ('b->b', 'B->B', 'h->h', 'H->H', 'i->i', 'I->I', 'l->l', 'L->L', 'q->q',
     'Q->Q', 'e->e', 'f->f', 'd->d', 'F->F', 'D->D'),
    'out0 = in0 * in0',
    doc='''Elementwise square function.

    .. seealso:: :data:`numpy.square`

    ''')

absolute = core.absolute

# TODO(beam2d): Implement it
# fabs

_unsigned_sign = 'out0 = in0 > 0'
_complex_sign = '''
if (in0.real() == 0) {
  out0 = (in0.imag() > 0) - (in0.imag() < 0);
} else {
  out0 = (in0.real() > 0) - (in0.real() < 0);
}
'''
sign = core.create_ufunc(
    'cupy_sign',
    ('b->b', ('B->B', _unsigned_sign), 'h->h', ('H->H', _unsigned_sign),
     'i->i', ('I->I', _unsigned_sign), 'l->l', ('L->L', _unsigned_sign),
     'q->q', ('Q->Q', _unsigned_sign), 'e->e', 'f->f', 'd->d',
     ('F->F', _complex_sign), ('D->D', _complex_sign)),
    'out0 = (in0 > 0) - (in0 < 0)',
    doc='''Elementwise sign function.

    It returns -1, 0, or 1 depending on the sign of the input.

    .. seealso:: :data:`numpy.sign`

    ''')

_float_preamble = '''
#ifndef NAN
#define NAN __int_as_float(0x7fffffff)
#endif
'''
_float_maximum = ('out0 = (isnan(in0) | isnan(in1)) ? out0_type(NAN) : '
                  'out0_type(max(in0, in1))')
maximum = core.create_ufunc(
    'cupy_maximum',
    ('??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
     'LL->L', 'qq->q', 'QQ->Q',
     ('ee->e', _float_maximum),
     ('ff->f', _float_maximum),
     ('dd->d', _float_maximum),
     ('FF->F', _float_maximum),
     ('DD->D', _float_maximum)),
    'out0 = max(in0, in1)',
    preamble=_float_preamble,
    doc='''Takes the maximum of two arrays elementwise.

    If NaN appears, it returns the NaN.

    .. seealso:: :data:`numpy.maximum`

    ''')

_float_minimum = ('out0 = (isnan(in0) | isnan(in1)) ? out0_type(NAN) : '
                  'out0_type(min(in0, in1))')
minimum = core.create_ufunc(
    'cupy_minimum',
    ('??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
     'LL->L', 'qq->q', 'QQ->Q',
     ('ee->e', _float_minimum),
     ('ff->f', _float_minimum),
     ('dd->d', _float_minimum),
     ('FF->F', _float_minimum),
     ('DD->D', _float_minimum)),
    'out0 = min(in0, in1)',
    preamble=_float_preamble,
    doc='''Takes the minimum of two arrays elementwise.

    If NaN appears, it returns the NaN.

    .. seealso:: :data:`numpy.minimum`

    ''')

fmax = core.create_ufunc(
    'cupy_fmax',
    ('??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
     'LL->L', 'qq->q', 'QQ->Q',
     ('ee->e', 'out0 = fmax(in0, in1)'),
     ('ff->f', 'out0 = fmax(in0, in1)'),
     ('dd->d', 'out0 = fmax(in0, in1)'),
     'FF->F', 'DD->D'),
    'out0 = max(in0, in1)',
    doc='''Takes the maximum of two arrays elementwise.

    If NaN appears, it returns the other operand.

    .. seealso:: :data:`numpy.fmax`

    ''')

fmin = core.create_ufunc(
    'cupy_fmin',
    ('??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
     'LL->L', 'qq->q', 'QQ->Q',
     ('ee->e', 'out0 = fmin(in0, in1)'),
     ('ff->f', 'out0 = fmin(in0, in1)'),
     ('dd->d', 'out0 = fmin(in0, in1)'),
     'FF->F', 'DD->D'),
    'out0 = min(in0, in1)',
    doc='''Takes the minimum of two arrays elementwise.

    If NaN appears, it returns the other operand.

    .. seealso:: :data:`numpy.fmin`

    ''')

_nan_to_num_preamble = '''
template <class T>
__device__ T nan_to_num(T x, T large) {
    if (isnan(x))
        return 0;
    if (isinf(x))
        return copysign(large, x);
    return x;
}

template <class T>
__device__ complex<T> nan_to_num(complex<T> x, T large) {
    T re = nan_to_num(x.real(), large);
    T im = nan_to_num(x.imag(), large);
    return complex<T>(re, im);
}
'''

nan_to_num = core.create_ufunc(
    'cupy_nan_to_num',
    ('?->?', 'b->b', 'B->B', 'h->h', 'H->H',
     'i->i', 'I->I', 'l->l', 'L->L', 'q->q', 'Q->Q',
     ('e->e',
      'out0 = nan_to_num(in0, float16(32 * 0x7FF))'),
     ('f->f',
      'out0 = nan_to_num(in0, __int_as_float(0x7F800000 - 1))'),
     ('d->d',
      'out0 = nan_to_num(in0, __longlong_as_double(0x7FF0000000000000 - 1))'),
     ('F->F',
      'out0 = nan_to_num(in0, __int_as_float(0x7F800000 - 1))'),
     ('D->D',
      'out0 = nan_to_num(in0, __longlong_as_double(0x7FF0000000000000 - 1))')),
    'out0 = in0',
    preamble=_nan_to_num_preamble,
    doc='''Elementwise nan_to_num function.

    .. seealso:: :data:`numpy.nan_to_num`

    ''')

# TODO(okuta): Implement real_if_close

# TODO(okuta): Implement interp

reason: import cupy,from cupy


text:
import sqlalchemy as sa
from enum import Enum
import inspect
import numpy as np
from astropy import units as u, wcs

from baselayer.app.env import load_env

from . import facility_apis

from sncosmo.bandpasses import (
    Bandpass, _BANDPASSES, _BANDPASS_INTERPOLATORS, read_bandpass)
from sncosmo.magsystems import _MAGSYSTEMS

_, cfg = load_env()

def tophat_bandpass_um(ctr, width, name=None):
    """Create a tophat Bandpass centered at `ctr` with width `width` (both
    in microns)."""

    wave = np.array([ctr - width / 2.0, ctr + width / 2.0])
    trans = np.array([1.0, 1.0])
    return Bandpass(wave, trans, wave_unit=u.micron, name=name)

    trans = np.array([1.0, 1.0])
    return Bandpass(wave, trans, wave_unit=u.nano, name=name)

# GRANDMA FILTERS

# find values for C and w
grandmafilters_meta = {
    'filterset': 'grandma filters',
    'retrieved': '15 Nov 2021',
    'description': 'grandma filters from Alldata.txt'}
for name, ctr, width in [('grandma::B',  0.445, 0.01),
                    ('grandma::C',  0.01, 0.01),
                    ('grandma::G', 0.464, 0.01),
                    ('grandma::I', 0.806, 0.01),
                    ('grandma::L', 3.450, 0.01),
                    ('grandma::N', 10.5, 0.01),
                    ('grandma::R', 0.658, 0.01),
                    ('grandma::V', 0.551, 0.01),
                    ('grandma::w',  0.01, 0.01)]:
    _BANDPASSES.register_loader(name, tophat_bandpass_um,
                                args=(ctr, width), meta=grandmafilters_meta)

def force_render_enum_markdown(values):
    return ', '.join(list(map(lambda v: f'`{v}`', values)))

ALLOWED_SPECTRUM_TYPES = tuple(
    cfg.get('spectrum_types.types', ['source', 'host', 'host_center'])
)
ALLOWED_MAGSYSTEMS = tuple(val['name'] for val in _MAGSYSTEMS.get_loaders_metadata())
ALLOWED_BANDPASSES = tuple(val['name'] for val in _BANDPASSES.get_loaders_metadata())
THUMBNAIL_TYPES = (
    'new',
    'ref',
    'sub',
    'sdss',
    'dr8',
    'ps1',
    "new_gz",
    'ref_gz',
    'sub_gz',
)
INSTRUMENT_TYPES = ('imager', 'spectrograph', 'imaging spectrograph')
FOLLOWUP_PRIORITIES = ('1', '2', '3', '4', '5')
FOLLOWUP_HTTP_REQUEST_ORIGINS = ('remote', 'skyportal')
LISTENER_CLASSNAMES = [
    k
    for k, v in facility_apis.__dict__.items()
    if inspect.isclass(v)
    and issubclass(v, facility_apis.Listener)
    and v is not facility_apis.Listener
]

LISTENER_CLASSES = [getattr(facility_apis, c) for c in LISTENER_CLASSNAMES]

allowed_spectrum_types = sa.Enum(
    *ALLOWED_SPECTRUM_TYPES, name='spectrumtypes', validate_strings=True
)
default_spectrum_type = cfg.get('spectrum_types.default', "source")

allowed_magsystems = sa.Enum(
    *ALLOWED_MAGSYSTEMS, name="magsystems", validate_strings=True
)
allowed_bandpasses = sa.Enum(
    *ALLOWED_BANDPASSES, name="bandpasses", validate_strings=True
)
thumbnail_types = sa.Enum(
    *THUMBNAIL_TYPES, name='thumbnail_types', validate_strings=True
)
instrument_types = sa.Enum(
    *INSTRUMENT_TYPES, name='instrument_types', validate_strings=True
)

followup_priorities = sa.Enum(
    *FOLLOWUP_PRIORITIES, name='followup_priorities', validate_strings=True
)

api_classnames = sa.Enum(
    *[
        k
        for k, v in facility_apis.__dict__.items()
        if inspect.isclass(v)
        and issubclass(v, facility_apis.FollowUpAPI)
        and v is not facility_apis.FollowUpAPI
    ],
    name='followup_apis',
    validate_strings=True,
)

listener_classnames = sa.Enum(
    *LISTENER_CLASSNAMES,
    name='followup_listeners',
    validate_strings=True,
)

py_allowed_spectrum_types = Enum('spectrumtypes', ALLOWED_SPECTRUM_TYPES)
py_allowed_magsystems = Enum('magsystems', ALLOWED_MAGSYSTEMS)
py_allowed_bandpasses = Enum('bandpasses', ALLOWED_BANDPASSES)
py_thumbnail_types = Enum('thumbnail_types', THUMBNAIL_TYPES)
py_followup_priorities = Enum('priority', FOLLOWUP_PRIORITIES)

sqla_enum_types = [
    allowed_spectrum_types,
    allowed_bandpasses,
    thumbnail_types,
    instrument_types,
    followup_priorities,
    api_classnames,
    listener_classnames,
]

max_stars_repo_name: julietteBrt/skyportal

max_issues_repo_name: julietteBrt/skyportal

max_forks_repo_name: julietteBrt/skyportal

reason: import numpy,from astropy


text:
import numpy as np

def apply_corrections(data, operator, station_rigidity, data_keys):
    """
    main function for the corrections module. Takes data and an operator string, finds what corrections to apply and
    applies them
    :param station_rigidity: cutoff rigidity if GV of station taking measurements
    :param data: dataframe containing data to be corrected and correction data
    :param operator: string specifying the network operator
    :return: dataframe containing corrected data
    """
    # get the keys of data to correct
    correction_key_dict = set_corr_keys(operator)
    # get the correction factors
    correction_factors = get_corr_factors(data, correction_key_dict, station_rigidity)
    # correct the data
    for key in data_keys:
        data[key] = data[key] * correction_factors['p_corr'] * correction_factors['h_corr']
    return data

def set_corr_keys(operator):
    """
    function to accept an operator string and return the dataframe keys in a dictionary
    :param operator: string specifying the operator, must be 'COSMOS-UK' or 'COSMOS-US'
    :return: dictionary containing the pressure and humidity correction keys

    """
    if "COSMOS-US" == operator:
        key_dict = {'p_corr': 'PRESS', 'h_corr': None}
    elif "COSMOS-UK" == operator:
        key_dict = {'p_corr': 'PA', 'h_corr': 'Q'}
    else:
        raise KeyError('%s is not a valid operator key' % operator)

    return key_dict

def get_corr_factors(data, key_dict, station_rigidity):
    """
    takes some data and an operator and returns correction factors
    :param station_rigidity: cutoff rigidity of station taking measurements (GV)
    :param data: dataframe containing data and correction data
    :param key_dict: dictionary containing correction keys
    :return: dictionary containing the correction factors
    """
    if key_dict['p_corr'] is not None:
        p_corr = pressure_correction(data[key_dict['p_corr']], station_rigidity)
    else:
        p_corr = np.full(len(data), 1)

    if key_dict['h_corr'] is not None:
        h_corr = humidity_correction(data[key_dict['h_corr']])
    else:
        h_corr = np.full(len(data), 1)

    ret_dict = {'p_corr': p_corr, 'h_corr': h_corr}

    return ret_dict

def pressure_correction(pressure, rigidity):
    """
    function to get pressure correction factors, given a pressure time series and rigidity value for the station
    :param pressure: time series of pressure values over the time of the data observations
    :param rigidity: cut-off rigidity of the station making the observations
    :return: series of correction factors
    """
    p_0 = np.nanmean(pressure)

    pressure_diff = pressure - p_0
    # g cm^-2. See Desilets & Zreda 2003
    mass_attenuation_length = attenuation_length(p_0, rigidity)

    exponent = pressure_diff * mass_attenuation_length

    pressure_corr = np.exp(exponent)

    return pressure_corr

def humidity_correction(humidity):
    """

    :param humidity:
    :return:
    """

    mean_hum = np.nanmean(humidity)
    hum_change = humidity - mean_hum

    hum_corr = 1 + (0.0054 * hum_change)

    return hum_corr

def attenuation_length(pressure, rigidity):
    """

    :param pressure:
    :param rigidity:
    :return:
    """
    # constants from Desilets & Zreda 2003, Earth & Planetary Science Letters 206
    c_0 = 5.4196 * (10 ** (-3))
    c_1 = 2.2082 * (10 ** (-4))
    c_2 = -5.1952 * (10 ** (-7))
    c_3 = 7.2062 * (10 ** (-6))
    c_4 = -1.9702 * (10 ** (-6))
    c_5 = -9.8334 * (10 ** (-9))
    c_6 = 3.4201 * (10 ** (-9))
    c_7 = 4.9898 * (10 ** (-12))
    c_8 = -1.7192 * (10 ** (-12))
    depth = pressure * 0.981
    beta = c_0 \
           + c_1 * rigidity \
           + c_2 * rigidity ** 2 \
           + (c_3 + c_4 * rigidity) * depth \
           + (c_5 + (c_6 * rigidity)) * depth ** 2 \
           + (c_7 + c_8 * rigidity) * (depth ** 3)

    return beta


text:
# coding: utf-8

import chainer

# Network definition

class A(chainer.Chain):

    def __init__(self):
        super(A, self).__init__()

    def forward(self, x):
        return x

# ======================================

import ch2o

if __name__ == '__main__':
    import numpy as np
    np.random.seed(314)

    model = A()

    x = np.random.rand(5, 7).astype(np.float32)
    x = [x]
    ch2o.generate_testcase(model, x)

max_stars_repo_name: disktnk/chainer-compiler

max_issues_repo_name: disktnk/chainer-compiler

max_forks_repo_name: disktnk/chainer-compiler


text:
import numpy as np
import torch
import torch.nn.functional as F
from torch import nn

class shifted_softplus(torch.nn.Module):

    def __init__(self):
        super(shifted_softplus, self).__init__()

    def forward(self, input):
        return F.softplus(input) - np.log(2.0)

class Swish(torch.nn.Module):
    def __init__(self):
        super().__init__()

    def forward(self, x):
        return F.silu(x)

class LearnableSwish(torch.nn.Module):
    def __init__(self,
                 alpha=1.0,
                 beta=1.702):
        super().__init__()
        self.alpha_inv = nn.Parameter(
            torch.log(
                torch.exp(torch.Tensor([alpha])) - 1
            )
        )
        self.beta_inv = nn.Parameter(
            torch.log(
                torch.exp(torch.Tensor([beta])) - 1
            )
        )

    @property
    def alpha(self):
        return F.softplus(self.alpha_inv)

    @property
    def beta(self):
        return F.softplus(self.beta_inv)

    def forward(self, x):
        device = x.device
        alpha = self.alpha.to(device)
        beta = self.beta.to(device)
        output = alpha / beta * F.silu(beta * x)

        return output

max_stars_repo_name: jkaraguesian/NeuralForceField

max_issues_repo_name: jkaraguesian/NeuralForceField

max_forks_repo_name: jkaraguesian/NeuralForceField


text:
# Authors: Douglas Buisson, Ryan-Rhys Griffiths
"""
This script generates simulated light curves. Timmer and Konig's algorithm generates Gaussian flux distributions.
"""

import pickle

import numpy as np
from matplotlib import pyplot as p

from fourier_methods import psd
from ts_gen import ts_gen

RAW_DATA_PATH= '../raw_data/mkn335_xrt_w2_lc.dat'
PROCESSED_DATA_PATH = '../processed_data/'

if __name__ == '__main__':

    # Save the gapped and full lightcurves in these files.

    gapped_output_xray_dat = 'sim_curves/xray_lightcurves_new.dat'
    full_output_xray_dat = 'sim_curves/xray_lightcurves_no_gaps_new.dat'

    gapped_output_uv_dat = f'sim_curves/w2_lightcurves.dat'
    full_output_uv_dat = f'sim_curves/w2_lightcurves_no_gaps.dat'

    # Read data
    with open('../processed_data/xray/x_ray_times.pickle', 'rb') as handle:
        xtg = pickle.load(handle)  # timings for x-ray simulatons
    with open('../processed_data/uv/uv_fl_times.pickle', 'rb') as handle:
        wtg = pickle.load(handle)  # timings for UVW2 simulations
    with open('../processed_data/xray/x_ray_band_count_rates.pickle', 'rb') as handle:
        xrg = pickle.load(handle)  # x-ray band count rates
    with open('../processed_data/uv/uv_band_flux.pickle', 'rb') as handle:
        wrg = pickle.load(handle)  # UVW2 flux values
        wrg = wrg/5.976e-16  # convert flux to count rate (erg/cm^2/s/Angstrom to cts/s)

    # save the filtered timings. This step is probably now unnecessary as the filtering is done in process_data.py
    with open(PROCESSED_DATA_PATH + 'xray_simulations/x_ray_sim_times.pickle', 'wb') as handle:
        pickle.dump(xtg, handle, protocol=pickle.HIGHEST_PROTOCOL)

    # X-ray PSD

    splits = np.where(xtg[1:] - xtg[:-1] > 16)[0] + 1
    tregs = np.array([])
    xregs = np.array([])
    prev_split = 0
    for split in splits:
        if split - prev_split > 10:
            tsec = xtg[prev_split:split]
            xsec = xrg[prev_split:split]
            treg = np.arange(tsec[0], tsec[-1], 4)
            xreg = np.interp(treg, tsec, xsec)
            tregs = np.concatenate([tregs, treg])
            xregs = np.concatenate([xregs, xreg])
            #p.psd(xreg, Fs=0.25, color='k')
            #p.plot(treg, xreg, 'k')
            prev_split = split

    xf, xp, xe = psd(tregs, np.log(xregs))  # calculate psd in log space

    # W2 PSD

    # save the filtered timings - this step is probably unnecessary now
    with open(PROCESSED_DATA_PATH + '/uv_simulations/uv_sim_times.pickle', 'wb') as handle:
        pickle.dump(wtg, handle, protocol=pickle.HIGHEST_PROTOCOL)

    splits = np.where(wtg[1:] - wtg[:-1] > 16)[0]+1
    tregs = np.array([])
    wregs = np.array([])
    prev_split = 0
    for split in splits:
        if split - prev_split > 10:
            tsec = wtg[prev_split:split]
            wsec = wrg[prev_split:split]
            treg = np.arange(tsec[0], tsec[-1], 4)
            wreg = np.interp(treg, tsec, wsec)
            tregs = np.concatenate([tregs, treg])
            wregs = np.concatenate([wregs, wreg])
            # p.psd(wreg, Fs=0.25, color='b')
            # p.plot(treg, wreg, 'b')
            prev_split = split

    wf, wp, we = psd(tregs, wregs)

    #p.xscale('log')
    #p.show()

    # Fit PSDs

    from scipy import optimize

    def linear_res(mc, freq, pow):
        return (mc[0]*freq + mc[1]) - pow

    xresult = optimize.least_squares(linear_res, [2, -3], args=(np.log(xf), np.log(xp)))
    wresult = optimize.least_squares(linear_res, [2, -3], args=(np.log(wf), np.log(wp)))

    # Plot PSDs

    p.errorbar(wf/86400, wf*wp, wf*we, fmt='b')
    p.errorbar(xf/86400, xf*xp, xf*xe, fmt='k')
    p.xscale('log')
    p.yscale('log')
    p.xlabel('Frequency (Hz)')
    p.ylabel(r'Power ($fP(f)$)')
    #p.show()

    # Simulate lightcurves

    # X-ray

    nsims = 1000

    # Shape of power spectrum
    freq = np.array([1e-15, 1e5])
    pow = freq**xresult.x[0]*np.exp(xresult.x[1])

    # Size/scale of real data
    mean_rate = np.mean(np.log(xrg))
    var_rate = np.std(np.log(xrg))

    # Time of each observation (nearest day/10)
    time = np.round(xtg*10)/10
    time -= min(time)

    # Exposure of data
    exposure = 1000.

    # Work with timesteps of tenths of days

    xlightcurves = np.empty([nsims, len(xrg)])  # preallocate the array of count rates
    xlightcurves_no_gaps = np.empty([nsims, 4390])  # obtained by checking size of rates

    for i in range(nsims):

        rate_inc = (ts_gen(2**19, freq=freq, pow=pow, loginterp=1))[500:(50000+500)]
        rates_all = np.exp((rate_inc - np.mean(rate_inc)) / np.std(rate_inc)*var_rate + mean_rate)

        gap_rates = rates_all[(time*10).astype(int)]
        np.random.seed(i)
        gap_rateso = np.float_(np.random.poisson(gap_rates * exposure)) / exposure
        xlightcurves[i, :] = gap_rates

        full_rates = rates_all[:(time[-1].astype(int) * 10 + 1)]
        np.random.seed(i)
        full_rates = np.float_(np.random.poisson(full_rates * exposure)) / exposure
        xlightcurves_no_gaps[i, :] = full_rates[:43891: 10]  # output in units of band count rate

        print(i)

    gap_file = open(gapped_output_xray_dat, 'wb')
    pickle.dump(xlightcurves, gap_file)
    gap_file.close()
    #ascii.write(xlightcurves, output_xray_dat)
    no_gap_file = open(full_output_xray_dat, 'wb')
    pickle.dump(xlightcurves_no_gaps, no_gap_file)
    no_gap_file.close()
    #ascii.write(xlightcurves_no_gaps, output_xray_dat)

    # W2 band

    # Shape of power spectrum
    freq = np.array([1e-15, 1e5])
    pow = freq**wresult.x[0]*np.exp(wresult.x[1])

    # Size/scale of real data
    mean_rate = np.mean(np.log(wrg))
    var_rate = np.std(np.log(wrg))

    # Time of each observation (nearest day/10)
    time = np.round(wtg*10)/10
    time -= min(time)

    # Exposure of data
    exposure = 1000.

    # Work with timesteps of tenths of days

    wlightcurves = np.empty([nsims, len(wrg)])
    wlightcurves_no_gaps = np.empty([nsims, 4390])

    for i in range(nsims):
        rate_inc = (ts_gen(2**19, freq=freq, pow=pow, loginterp=1))[500:(50000+500)]
        rates_all = np.exp((rate_inc - np.mean(rate_inc)) / np.std(rate_inc)*var_rate + mean_rate)

        rates = rates_all[(time*10).astype(int)]
        np.random.seed(i)
        rates = np.float_(np.random.poisson(rates * exposure)) / exposure
        wlightcurves[i, :] = rates

        rates = rates_all[:(time[-1].astype(int) * 10 + 1)]
        np.random.seed(i)
        rates = np.float_(np.random.poisson(rates * exposure)) / exposure
        wlightcurves_no_gaps[i, :] = rates[:43891: 10]  # every day, units of count rate, 4390 values

        print(i)

    # we use pickle because ascii can't deal with large files

    gap_file = open(gapped_output_uv_dat, 'wb')
    wlightcurves *= 5.976e-16  # Convert back to units of flux
    pickle.dump(wlightcurves, gap_file)
    gap_file.close()
    no_gap_file = open(full_output_uv_dat, 'wb')
    wlightcurves_no_gaps *= 5.976e-16
    pickle.dump(wlightcurves_no_gaps, no_gap_file)
    no_gap_file.close()

reason: import numpy,from scipy


text:
## utils for reading the emails from various sources
#https://github.com/tuanavu/udacity-course/blob/master/intro_to_machine_learning/lesson/lesson_10_text_learning/vectorize_text.py
#https://gtraskas.github.io/post/spamit/
#https://towardsdatascience.com/k-means-clustering-8e1e64c1561c
import os
import pickle
from itertools import chain
import random
import nltk
import numpy as np
from nltk import tokenize

ENRON_PATH = 'data/Enron/'

def load_enron_vectorized(enron_ids, n, shared=False):

    enrons = [load_text(enron_id) for enron_id in enron_ids]
    flattened_enron = chain.from_iterable(chain.from_iterable(enrons))
    mfw = most_freq_words(flattened_enron, n)

    count = sum((len(h)+len(s) for h,s in enrons)) #number of emails

    X = np.ndarray([count, n])
    Y = np.ndarray([count, 2])

    i = 0
    for nom,anom in enrons:
        X[i:i+len(nom),:] = vectorize_text(nom, mfw)
        Y[i:i+len(nom),:] = vectorize_class(nom, is_anom=False)
        i += len(nom)

        X[i:i+len(anom),:] = vectorize_text(anom, mfw)
        Y[i:i+len(anom),:] = vectorize_class(anom, is_anom=True)
        i += len(anom)

    return X,Y

def vectorize_text(all_text, most_freq_words):

    n = len(all_text)
    m = len(most_freq_words)

    V = np.zeros([n, m], dtype='float32')
    for i,text in enumerate(all_text):
        words = tokenize.word_tokenize(text)
        l_words = {w.lower() for w in words}

        for j,mfw in enumerate(most_freq_words):
            if mfw in l_words:
                V[i][j] = 1.0

    return V

def vectorize_class(all_text, is_anom):

    n = len(all_text)
    V = np.zeros([n, 2])

    if is_anom:
        V[:,1]=1.0
    else:
        V[:,0]=1.0
    return V

def load_text(enron_id):

    if not (isinstance(enron_id, int) and 1<=enron_id<=6):
        raise ValueError('enron_id must be an int from 1 to 6')

    nom = []
    anom = []

    nom_path = ENRON_PATH+'enron{}/ham/'.format(enron_id)
    try:
        for filepath in os.listdir(nom_path):
            try:
                with open(nom_path+filepath, 'r') as fd:
                    txt = fd.read()
                    nom.append(txt)
            except UnicodeDecodeError:
                print('utf-8 error. Skipping file {}'.format(nom_path+filepath))
    except FileNotFoundError:
        print('No enron nom dataset found for id {}'.format(enron_id))
        raise

    anom_path = ENRON_PATH+'enron{}/spam/'.format(enron_id)
    try:
        for filepath in os.listdir(anom_path):
            try:
                with open(anom_path+filepath, 'r') as fd:
                    txt = fd.read()
                    anom.append(txt)
            except UnicodeDecodeError:
                print('utf-8 error. Skipping file {}'.format(anom_path+filepath))
    except FileNotFoundError:
        print('No enron anom dataset found for id {}'.format(enron_id))
        raise

    return (nom, anom)

def most_freq_words(all_text, n):

    freq_dist = nltk.FreqDist()

    for i,text in enumerate(all_text):
        words = tokenize.word_tokenize(text)

        l_words = [w.lower() for w in words
                   if w.isalnum() and not w.isnumeric()]

        freq_dist.update(l_words)

    top_n = [w[0] for w in freq_dist.most_common(n)]
    top_n.sort()
    return top_n


text:
from collections.abc import Iterator, Sequence, Sized

import numpy as np

class BatchIterator(Iterator, Sized):

    def __init__(self, iterator, size):
        if not isinstance(size, int):
            raise TypeError("iterator size must be an int value")
        self._iterator = iterator
        self._size = size

    def __iter__(self):
        return self

    def __next__(self):
        return next(self._iterator)

    def __len__(self):
        return self._size

def _get_row_iterator(samples, batches):
    for indices in batches:
        yield np.take(samples, indices, axis=0)

def _get_col_iterator(columns, batches):
    for indices in batches:
        yield tuple(_take(column, indices) for column in columns)

def _take(column, indices):
    if isinstance(column, (list, tuple)) \
            and isinstance(column[0], np.ndarray):
        return type(column)(column[idx] for idx in indices)
    else:
        return np.take(column, indices, axis=0)

class Dataset(Sequence):
    """Immutable Class"""

    def __init__(self, *samples):
        self._n_cols = len(samples)
        if self._n_cols == 0:
            self._columns = ([],)
            self._n_cols = 1
        elif self._n_cols == 1:
            if isinstance(samples[0], (list, tuple)):
                columns = tuple(
                    np.array(column)
                    if isinstance(column[0], np.ndarray)
                    and column[0].ndim <= 1 else column
                    for column in map(tuple, zip(*samples[0])))
                self._columns = columns
                self._n_cols = len(columns)
            else:
                self._columns = (samples[0],)
        elif self._n_cols > 1:
            self._columns = samples
        self._samples = list(map(tuple, zip(*self._columns)))
        self._len = len(self._samples)

    def __iter__(self):
        return (sample for sample in self._samples)

    def batch(self, size, shuffle=False, colwise=True):
        indices = np.arange(self._len)
        if shuffle:
            np.random.shuffle(indices)
        batches = [indices[i: i + size]
                   for i in range(0, self._len, size)]
        if colwise:
            iterator = _get_col_iterator(self._columns, batches)
        else:
            iterator = _get_row_iterator(self._samples, batches)
        return BatchIterator(iterator, size=len(batches))

    def __len__(self):
        return self._len

    def __getitem__(self, key):
        return self._samples.__getitem__(key)

    def take(self, indices):
        return np.take(self._samples, indices, axis=0)

    def __reversed__(self):
        return self._samples.__reversed__()

    def __contains(self, sample):
        return self._samples.__contains__(sample)

    def index(self, sample):
        return self._samples.index(sample)

    def count(self, sample):
        return self._samples.count(sample)

    def __repr__(self):
        return repr(self._samples)

    def __str__(self):
        return str(self._samples)

    def cols(self):
        return self._columns

    @property
    def size(self):
        return self.__len__()

    @property
    def n_cols(self):
        return self._n_cols

class BucketDataset(Dataset):

    def __init__(self, *samples, key=0, equalize_by_key=False):
        super().__init__(*samples)
        self._key = key
        self._equalize = equalize_by_key

    def batch(self, size, shuffle=False, colwise=True):
        key = self._key
        lengths = ((float(len(sample[key])), index)
                   for index, sample in enumerate(self))
        if shuffle:
            # Add noise to the lengths so that the sorting isn't deterministic.
            lengths = ((length + np.random.random(), index)
                       for length, index in lengths)
            # Randomly choose the order from (ASC, DESC).
            lengths = sorted(lengths, key=lambda x: x[0],
                             reverse=np.random.random() > .5)
        # Bucketing
        if self._equalize:
            buckets = []
            bucket = []
            accum_length = 0
            for length, index in lengths:
                if accum_length + int(length) > size:
                    buckets.append(np.array(bucket))
                    bucket = []
                    accum_length = 0
                bucket.append(index)
                accum_length += int(length)
            if bucket:
                buckets.append(np.array(bucket))
        else:
            lengths = list(lengths)
            buckets = [np.array([index for _, index in lengths[i: i + size]])
                       for i in range(0, self._len, size)]
        if shuffle:
            np.random.shuffle(buckets)
        if colwise:
            iterator = _get_col_iterator(self._columns, buckets)
        else:
            iterator = _get_row_iterator(self._samples, buckets)
        return BatchIterator(iterator, size=len(buckets))


text:
import pickle

import numpy as np
import torch

from detection.src.yolov3.utils.datasets import ListDataset

class DetectionSetDataManager():
    """
    Data Manager used for YOLOMAML
    """
    def __init__(self, n_way, n_support, n_query, n_episode, image_size):
        """

        Args:
            n_way (int): number of different classes in a detection class
            n_support (int): number of images in the support set with an instance of one class,
            for each of the n_way classes
            n_query (int): number of images in the query set with an instance of one class,
            for each of the n_way classes
            n_episode (int): number of episodes per epoch
            image_size (int): size of images (square)
        """
        self.n_way = n_way
        self.n_support = n_support
        self.n_query = n_query
        self.n_episode = n_episode
        self.image_size = image_size

    def get_data_loader(self, path_to_data_file, path_to_images_per_label=None):
        """

        Args:
            path_to_data_file (str): path to file containing paths to images
            path_to_images_per_label (str): path to pkl file containing images_per_label dictionary (optional)

        Returns:
            DataLoader: samples data in the shape of a detection task
        """
        dataset = ListDataset(path_to_data_file, img_size=self.image_size)
        sampler = DetectionTaskSampler(
            dataset,
            self.n_way,
            self.n_support,
            self.n_query,
            self.n_episode,
            path_to_images_per_label,
        )
        data_loader = torch.utils.data.DataLoader(dataset,
                                                  batch_sampler=sampler,
                                                  num_workers=12,
                                                  collate_fn=dataset.collate_fn_episodic,
                                                  )
        return data_loader

def create_dict_images_per_label(data_source):
    """
            Compute and returns dictionary of images per label
            Args:
                data_source (ListDataset) : The data set containing the images
            Returns:
                dict: each key maps to a list of the images which contain at least one target which label is the key
            """
    images_per_label={}

    for index in range(len(data_source)):
        try:
            targets = data_source[index][2]
            if targets is not None:
                for target in targets:
                    label = int(target[1])
                    if label not in images_per_label:
                        images_per_label[label] = []
                    if len(images_per_label[label]) == 0 or images_per_label[label][-1] != index:
                        images_per_label[label].append(index)
                if index % 100 == 0:
                    print('{index}/{length_data_source} images considered'.format(
                        index=index,
                        length_data_source=len(data_source))
                    )
        except OSError:
            print('Corrupted image : {image_index}'.format(image_index=index))
    return images_per_label

class DetectionTaskSampler(torch.utils.data.Sampler):
    """
    Samples elements in detection episodes of defined shape.
    """
    def __init__(self, data_source, n_way, n_support, n_query, n_episodes, path_to_images_per_label=None):
        """

        Args:
            data_source (ListDataset): source dataset
            n_way (int): number of different classes in a detection class
            n_support (int): number of images in the support set with an instance of one class,
            for each of the n_way classes
            n_query (int): number of images in the query set with an instance of one class,
            for each of the n_way classes
            n_episodes (int): number of episodes per epoch
            path_to_images_per_label (str): path to a pickle file containing a dictionary of images per label
        """
        self.data_source = data_source
        self.n_way = n_way
        self.n_support = n_support
        self.n_query = n_query
        self.n_episodes = n_episodes

        self.images_per_label = self.get_images_per_label(path_to_images_per_label)
        self.label_list = self.get_label_list()

    def get_images_per_label(self, path):
        """
        Returns dictionary of images per label from a file if specified or compute it from scratch
        Args:
            path (str) : path to a pickle file containing a dictionary of images per label
        Returns:
            dict: each key maps to a list of the images which contain at least one target which label is the key
        """
        if path:
            with open(path, 'rb') as dictionary_file:
                images_per_label = pickle.load(dictionary_file)
        else:
            images_per_label = create_dict_images_per_label(self.data_source)

        return images_per_label

    def get_label_list(self):
        """

        Returns:
            list: list of appropriate labels, i.e. labels that are present in at least n_support+n_query images
        """
        label_list = []
        for label in self.images_per_label:
            if len(self.images_per_label[label]) >= self.n_support + self.n_query:
                label_list.append(label)
        return label_list

    def sample_labels(self):
        """

        Returns:
            numpy.ndarray: n_way labels sampled at random from all available labels
        """
        labels = np.random.choice(self.label_list, self.n_way, replace=False)
        return labels

    def sample_images_from_labels(self, labels):
        """
        For each label in labels, samples n_support+n_query images containing at least one box associated with label
        The first n_way elements of the returned tensor will be used to determine the sampled labels
        Args:
            labels (numpy.ndarray): labels from which images will be sampled

        Returns:
            torch.Tensor: length = n_way*(1+n_support+n_query) information about the labels,
            and indices of images constituting an episode
        """
        #TODO: images can appear twice
        images_indices = list(-labels-1)
        for label in labels:
            images_from_label = np.random.choice(
                self.images_per_label[label],
                self.n_support+self.n_query,
                replace=False
            )
            images_indices.extend(images_from_label)
        return torch.tensor(images_indices, dtype=torch.int32)

    def __len__(self):
        return self.n_episodes

    def __iter__(self):
        for i in range(self.n_episodes):
            labels = self.sample_labels()
            yield self.sample_images_from_labels(labels)

max_stars_repo_name: artificially-ai/FewShotVision

max_issues_repo_name: artificially-ai/FewShotVision

max_forks_repo_name: ebennequin/FewShotVision