text:
"""CLI submodule for creating a new dataset."""

from typing import List, Tuple
import argparse
import logging
import os

import numpy as np
import pandas as pd
import skimage.util

from ..io import EXTENSIONS
from ..io import basename
from ..io import grab_files
from ..io import load_image
from ..util import train_valid_split
from ._parseutil import CustomFormatter
from ._parseutil import FolderType
from ._parseutil import _add_utils

def _parse_args_create(
    subparsers: argparse._SubParsersAction, parent_parser: argparse.ArgumentParser,
):
    """Subparser for creation."""
    parser = subparsers.add_parser(
        "create",
        parents=[parent_parser],
        formatter_class=CustomFormatter,
        add_help=False,
        description=(
            "\U0001F4BE Creation submodule \U0001F4BE\n\n"
            "Create a custom dataset with raw files and corresponding labels. "
            "Relies on labeling output from FIJI that was saved with the provided macro "
            'described here "https://github.com/BBQuercus/deepBlink/wiki/Datasets".'
        ),
        help="\U0001F4BE Create a new dataset from raw files.",
    )
    group1 = parser.add_argument_group("Required")
    group1.add_argument(
        "-i",
        "--input",
        required=True,
        type=FolderType(),
        help=(
            "Path to the directory containing raw images. "
            "Note that only the specified filetypes will be processed. "
            f"[required] [filetypes: {', '.join(EXTENSIONS)}]"
        ),
    )
    group2 = parser.add_argument_group("Optional")
    group2.add_argument(
        "-l",
        "--labels",
        type=FolderType(),
        help=(
            "Path to the directory containing labels in csv format. "
            "The default path accounts for using the FIJI macro described on the wiki. "
            "[default: --INPUT/labels/]"
        ),
    )
    group2.add_argument(
        "-n",
        "--name",
        default="dataset",
        type=str,
        help=(
            "Custom dataset name. "
            'The file extension "npz" will be added automatically. '
            '[default: "dataset"]'
        ),
    )
    group2.add_argument(
        "-s",
        "--size",
        default=None,
        type=int,
        help=(
            "Image crop size. "
            "If given, crops all images into the specified size. "
            "Will crop non-overlapping and ignore areas that did not get covered."
            "deepBlink requires square images powers of 2, such as 256, 512... "
            "[default: None]"
        ),
    )
    group2.add_argument(
        "-vs",
        "--validsplit",
        default=0.2,
        type=float,
        help=(
            "Validation split. "
            "Split percentage (scaled between 0 - 1) of validation vs. train set. "
            "Note the validation split is done after splitting test and trainval. "
            "[default: 0.2]"
        ),
    )
    group2.add_argument(
        "-ts",
        "--testsplit",
        default=0.2,
        type=float,
        help=(
            "Testing split. "
            "Split percentage (scaled between 0 - 1) of test vs. trainval set. "
            "[default: 0.2]"
        ),
    )
    _add_utils(parser)

class HandleCreate:
    """Handle creation submodule for CLI.

    Args:
        arg_input: Path to folder with images.
        arg_labels: Path to folder with labels.
        arg_name: Name of dataset file to be saved.
        arg_size: Size of image to be cropped.
        arg_testsplit: Test vs. Trainval split percentage.
        arg_testsplit: Valid vs. Train split percentage.
        logger: Logger to log verbose output.
    """

    # pylint: disable=too-many-instance-attributes
    # We require many attributes for the splits

    def __init__(
        self,
        arg_input: str,
        arg_labels: str,
        arg_name: str,
        arg_size: int,
        arg_testsplit: int,
        arg_validsplit: int,
        logger: logging.Logger,
    ):
        self.raw_input = arg_input
        self.raw_labels = arg_labels
        self.raw_name = arg_name
        self.img_size = arg_size
        self.test_split = arg_testsplit
        self.valid_split = arg_validsplit
        self.logger = logger
        self.logger.info("\U0001F5BC starting creation submodule")

        self.abs_input = os.path.abspath(self.raw_input)
        self.extensions = EXTENSIONS

    def __call__(self):
        """Run dataset creation."""
        self.get_file_lists()
        self.train_val_test_split()
        self.save_npz()
        self.logger.info(f"\U0001F3C1 dataset created at {self.fname_out}")

    @property
    def abs_labels(self):
        """Return absolute path to directory with labels."""
        if self.raw_labels is not None:
            path = os.path.abspath(self.raw_labels)
            self.logger.debug(f"using provided label path at {path}")
        elif os.path.isdir(os.path.join(self.abs_input, "labels")):
            path = os.path.join(self.abs_input, "labels")
            self.logger.debug(f"using default label path at {path}")
        else:
            self.logger.debug(
                f"no labels found in default {self.abs_input} or input {self.raw_labels}"
            )
            raise ValueError(
                (
                    "\U0000274C No label path found.\t"
                    "Please use a directory called 'labels' in input or use the '--labels' flag."
                )
            )
        return path

    @property
    def fname_out(self):
        """Return the absolute path to the dataset."""
        raw_name = self.raw_name
        if raw_name is not None:
            if os.path.isfile(raw_name):
                path = raw_name
                self.logger.warning(
                    f"\U000026A0 input name {raw_name} is already a file"
                )
            elif os.path.isdir(raw_name):
                path = os.path.join(raw_name, "dataset.npz")
                self.logger.warning(
                    f"\U000026A0 input name {raw_name} is already a file"
                )
            else:
                fname = raw_name[:-4] if raw_name.endswith(".npz") else raw_name
                path = os.path.join(self.abs_input, f"{fname}.npz")
                self.logger.debug(f"using given name {raw_name}")
        else:
            path = os.path.join(self.raw_name, "dataset.npz")
            self.logger.debug(f"using default output at {path}")
        return path

    @property
    def image_label_lists(self):
        """Return lists with all images and labels."""
        fname_images = grab_files(self.abs_input, self.extensions)
        fname_labels = grab_files(self.abs_labels, extensions=("csv",))

        self.logger.debug(f"images - found {len(fname_images)} files: {fname_images}")
        self.logger.debug(f"labels - found {len(fname_labels)} files: {fname_labels}")

        size = self.img_size

        images = []
        labels = []
        for image, label in zip(fname_images, fname_labels):
            if basename(image) != basename(label):
                self.logger.warning(
                    f"\U0000274C file basenames do not match! {image} != {label}"
                )
            df = pd.read_csv(label, index_col=0)
            if len(df) <= 1:
                self.logger.warning(
                    f"\U000026A0 labels for {label} empty. will not be used"
                )
                continue

            img = load_image(image, is_rgb=False)

            if all(shape > size for shape in img.shape):
                images.append(load_image(image, is_rgb=False))
                labels.append(df)

        if len(images) != len(labels):
            raise ValueError(
                f"Number of images and labels must match. {len(images)} != {len(labels)}."
            )
        self.logger.debug(f"using {len(images)} non-empty files")

        return images, labels

    @staticmethod
    def convert_labels(image: np.ndarray, df: pd.DataFrame) -> pd.DataFrame:
        """Pre-processes labels to be used in deepBlink.

        Renames X/Y to c/r respectively for easier handling with rearrangement to r/c.
        Rounds coordinates on borders to prevent Fiji out-of bounds behavior.
        """
        if "X" in df.columns:
            df = df.rename(columns={"X": "c", "Y": "r"})[["r", "c"]]
        elif "POSITION_X" in df.columns:
            df = df.rename(columns={"POSITION_X": "c", "POSITION_Y": "r"})[["r", "c"]]
            df = df.reset_index(drop=True)
        else:
            raise ValueError(
                "Format of input labels not recognized. "
                "Requires X,Y or POSITION_X,POSITION_Y in columns. "
                f"Columns found are: {df.columns.to_list()}."
            )

        for name, var in zip(["r", "c"], image.shape):
            df[name] = df[name].where(df[name] < var, var)
            df[name] = df[name].where(df[name] > 0, 0)

        return df

    def crop_image(
        self, image: np.ndarray, df: pd.DataFrame
    ) -> Tuple[List[np.ndarray], List[pd.DataFrame]]:
        """Crop a image / label pair to a uniform size and scale labels accordingly."""
        size = self.img_size
        if size is None:
            self.logger.debug(f"using unchanged size in {image.shape}")
            return [image], [df]

        if any([size > s for s in image.shape]):
            self.logger.debug(f"skipping image due to small size {image.shape}")
            return [], []

        windows = skimage.util.view_as_windows(
            image, window_shape=(size, size), step=size
        )

        img_list = []
        df_list = []
        for r, cimages in enumerate(windows):
            for c, img in enumerate(cimages):
                img_list.append(img)

                r_min = r * size
                c_min = c * size
                df_slice = df.loc[
                    (
                        (df["r"] >= r_min)
                        & (df["r"] <= r_min + size)
                        & (df["c"] >= c_min)
                        & (df["c"] <= c_min + size)
                    )
                ]
                df_norm = df_slice.copy()
                df_norm.loc[:, "r"] = df_slice["r"] - r_min
                df_norm.loc[:, "c"] = df_slice["c"] - c_min
                df_list.append(df_norm)

        self.logger.debug(
            f"converted original size of {image.shape} into {len(df_list)} crops"
        )
        return img_list, df_list

    def get_file_lists(self):
        """Create a file list after cropping over-sized images."""
        self.image_list = []
        self.label_list = []

        for image, label in zip(*self.image_label_lists):
            label_normalized = self.convert_labels(image=image, df=label)
            image_crops, label_crops = self.crop_image(image=image, df=label_normalized)
            self.image_list.extend(image_crops)
            self.label_list.extend(label_crops)

        self.logger.debug(
            f"images converted: {len(self.image_list)} == {len(self.label_list)}"
        )

    def train_val_test_split(self):
        """Two-step split of the input data into train/valid/test."""
        x_trainval, self.x_test, y_trainval, y_test = train_valid_split(
            x_list=self.image_list, y_list=self.label_list, valid_split=self.test_split
        )
        self.x_train, self.x_valid, y_train, y_valid = train_valid_split(
            x_list=x_trainval, y_list=y_trainval, valid_split=self.valid_split
        )
        self.logger.info(
            f"\U0001F4A6 images split: {len(self.x_train)} train, {len(self.x_valid)} valid, {len(self.x_test)} test"
        )

        # Convert DataFrame to Numpy array
        self.y_train = [y.values for y in y_train]
        self.y_valid = [y.values for y in y_valid]
        self.y_test = [y.values for y in y_test]

    def save_npz(self):
        """Save dataset splits as single npz file."""
        np.savez_compressed(
            self.fname_out,
            x_train=self.x_train,
            y_train=self.y_train,
            x_valid=self.x_valid,
            y_valid=self.y_valid,
            x_test=self.x_test,
            y_test=self.y_test,
        )


text:
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from mbtr.losses import LinRegLoss, QuantileLoss, QuadraticQuantileLoss
# package required for plotting the tree; installation issues under windows
# from networkx.drawing.nx_agraph import graphviz_layout
from time import time
from numba import jit, float32
from mbtr import LOSS_MAP
from tqdm import tqdm
from mbtr.utils import set_figure
from pandas import read_feather

class Tree:
    """
    Tree class. Fits both univarite and multivariate targets. It implements histogram search for the decision of the
    splitting points.

    :param n_q: number of quantiles for the split search
    :param min_leaf: minimum number of observations in one leaf. This parameter greatly affect generalization abilities.
    :param loss_type: loss type for choosing the best splits. Currently the following losses are implemented:

                      mse: mean squared error loss, a.k.a. L2, ordinary least squares

                      time_smoother: mse with an additional penalization on the second order differences of the response
                      function. Requires to pass also lambda_smooth parameter.

                      latent_variable: generate response function of dimension n_t from an arbitrary linear combination
                      of n_r responses. This requires to pass also S and precision pars.

                      linear_regression: mse with linear response function. Using this loss function, when calling fit
                      and predict methods, one must also pass x_lr as additional argument, which is the matrix of
                      features used to train the linear response inside the leaf (which can be different from the
                      features used to grow the tree, x).

                      fourier: mse with linear response function, fitted on the first n_harmonics (where the
                      fundamental has wave-lenght equal to the target output). This requires to pass also the
                      n_harmonics parameter.

                      quantile: quantile loss function, a.k.a. pinball loss. This requires to pass also the alphas
                      parameter, a list of quantiles to be fitted.

                      quadratic_quantile: quadratic quantile loss function tailored for trees. It has a non-discontinuos
                      derivative. This requires to pass also the alphas parameter, a list of quantiles to be fitted.
    :param lambda_weights: coefficient for the quadratic regularization of the response's parameters
    :param lambda_leaves: coefficient for the quadratic regularization of the total number of leaves. This is only used
    when the Tree is used as a weak learner by MBT.
    :param loss_kwargs: possible additional arguments for the loss function

    """
    def __init__(self, n_q: int = 10, min_leaf: int = 100, loss_type: str = "mse",
                 lambda_weights: float = 0.1, lambda_leaves: float = 0.1, **loss_kwargs):
        self.n_q = n_q
        self.min_leaf = min_leaf
        self.loss_pars = {**{"loss_type": loss_type, "lambda_weights": lambda_weights, "lambda_leaves": lambda_leaves,
                             **loss_kwargs}}
        self.g = nx.DiGraph()
        self.g.add_node('origin', leaf=False, loss=np.inf)
        self.loss = LOSS_MAP[self.loss_pars['loss_type']](**self.loss_pars)

    def fit(self, x, y, hessian=None, learning_rate=1.0, x_lr=None):
        """
        Fits a tree using the features specified in the matrix :math:`x\\in\\mathbb{R}^{n_{obs} \\times n_{f}}`, in order
        to predict the targets in the matrix :math:`y\\in\\mathbb{R}^{n_{obs} \\times n_{t}}`, where :math:`n_{obs}` is
        the number of observations, :math:`n_{f}` the number of features and :math:`n_{t}` the dimension of the target.

        :param x: feature matrix, np.ndarray.
        :param y: target matrix, np.ndarray.
        :param hessian: diagonals of the hessians :math:`\\in\\mathbb{R}^{n_{obs} \\times n_{t}}`. If None, each entry
                        is set equal to one (this will result in the default behaviour under MSE loss). Default: None
        :param learning_rate: learning rate used by the MBT instance. Default: 1
        :param x_lr: features for fitting the linear response inside the leaves. This is only required if a LinearLoss
                      is being used.

        """
        # set loss dimensionality
        self.loss.set_dimension(y.shape[1])

        if hessian is None:
            hessian = np.ones_like(y)
        if type(self.loss) == LinRegLoss:
            self._fit_node_lr(x, y, x_lr, 'origin', learning_rate=learning_rate, hessian=hessian)
        else:
            self._fit_node(x, y, 'origin', learning_rate=learning_rate, hessian=hessian)

    def predict(self, x, x_lr=None):
        """
        Predicts the target based on the feature matrix x (and linear regression features x_lr).

        param: x: feature matrix, np.ndarray.
        param: x_lr: linear regression feature matrix, np.ndarray.

        :return: target's predictions
        """
        r = []
        if x_lr is None:
            for x_i in x:
                r.append(self._predict_node(x_i, 'origin'))
        else:
            for x_i, x_lr_i in zip(x,x_lr):
                r.append(self._predict_node(x_i, 'origin', x_lr_i))
        return np.vstack(r)

    def _predict_node(self, x:np.ndarray, node, x_lr:np.ndarray=None):
        """
        Recursively explores the tree and retrieve the response

        :param x: matrix of covariates, (n_obs, n_vars)
        :param node: current node of the tree
        :param x_lr: matrix of linear regression covariates. Only required if we are fitting linear regressions inside
         the nodes

        :return: response: matrix of responses
        """
        if self.g.nodes[node]['leaf']:
            if type(self.loss) == LinRegLoss:
                response = np.hstack([x_lr,1]) @ self.g.nodes[node]['response']
            else:
                response = self.g.nodes[node]['response']
        else:
            v = self.g.nodes[node]['variable']
            l = self.g.nodes[node]['level']
            if x[v] < l:
                node = [i for i in self.g.successors(node) if self.g.nodes[i]['side'] == 'l']
            else:
                node = [i for i in self.g.successors(node) if self.g.nodes[i]['side'] == 'r']
            response = self._predict_node(x, node[0], x_lr)
        return response

    def _fit_node(self, x:np.ndarray, y:np.ndarray, node, learning_rate=1.0, hessian=None):
        """
        Recursively fits the tree using second order approximation (gradient and hessian) of the objective function
        :param x: matrix of covariates, (n_obs, n_vars)
        :param y: matrix of target, (n_obs, n_t)
        :param node: current node of the tree
        :param: learning_rate: discount factor. Undershoot the weak learner response to dampen the learning process
        """
        if x.shape[0] < self.min_leaf:
            self.g.nodes[node]['leaf'] = True
            self.g.nodes[node]['response'] = self.loss.eval_optimal_response(
                np.sum(y, axis=0), np.sum(hessian, axis=0)) * learning_rate
            return
        else:
            is_leaf = True
            self.g.nodes[node]['loss'] = self.loss.eval_optimal_loss(np.sum(y,axis=0) ** 2, np.sum(hessian, axis=0))
            current_loss = np.copy(self.g.nodes[node]['loss'])

            # for all the variables, find best split point
            for i in range(x.shape[1]):
                # find quantiles
                qs_i = np.unique(np.quantile(x[:,i],np.linspace(1/self.n_q,1-1/self.n_q,self.n_q)))
                G2_left, G2_right = leaf_stats(y, np.hstack([-np.inf,qs_i,np.inf]), x[:, i], order=1)
                H_left, H_right = leaf_stats(hessian, np.hstack([-np.inf, qs_i, np.inf]), x[:, i], order=2)

                # for all the quantiles of the current variables, find best quantile split
                for j, q_i in enumerate(qs_i):
                    loss = self.compute_loss(G2_left, G2_right, H_left, H_right, j)
                    if loss < current_loss:
                        is_leaf = False
                        current_loss = np.copy(loss)
                        self.g.nodes[node]['variable'] = i
                        self.g.nodes[node]['level'] = q_i

            # if a split reducing the total loss has been found, split the node at best variable/level
            self.g.nodes[node]['leaf'] = is_leaf
            if is_leaf:
                self.g.nodes[node]['response'] = self.loss.eval_optimal_response(
                    np.sum(y, axis=0), np.sum(hessian, axis=0)) * learning_rate
            else:
                node_name = self._add_new_node(node)
                filt_1 = x[:, self.g.nodes[node]['variable']] < self.g.nodes[node]['level']
                filt_2 = x[:, self.g.nodes[node]['variable']] >= self.g.nodes[node]['level']
                if (np.sum(filt_1) == 0) or (np.sum(filt_2) == 0):
                    self.g.nodes[node]['response'] = self.loss.eval_optimal_response(
                        np.sum(y, axis=0), np.sum(hessian, axis=0)) * learning_rate
                    self.g.nodes[node]['leaf'] = True
                    return
                else:
                    self._fit_node(x[filt_1, :], y[filt_1, :], node_name, learning_rate, hessian[filt_1, :])
                    self._fit_node(x[filt_2, :], y[filt_2, :], node_name+1, learning_rate, hessian[filt_2, :])

    def _fit_node_lr(self, x, y, x_lr, node, learning_rate, hessian):
        if x.shape[0] < self.min_leaf:
            self.g.nodes[node]['leaf'] = True
            self.g.nodes[node]['response'] = self.loss.eval_optimal_response(
                y, x_lr) * learning_rate
            return
        else:
            is_leaf = True
            self.g.nodes[node]['loss'] = self.loss.eval_optimal_loss(y, x_lr)
            current_loss = np.copy(self.g.nodes[node]['loss'])

            # for all the variables, find best split point
            for i in range(x.shape[1]):
                # find quantiles
                qs_i = np.unique(np.quantile(x[:, i], np.linspace(1/self.n_q, 1-1/self.n_q, self.n_q)))
                # for all the quantiles of the current variables, find best quantile split
                for j, q_i in enumerate(qs_i):
                    filt = x[:, i] <= q_i
                    loss_l = self.loss.eval_optimal_loss(y[filt,:], x_lr[filt,:])
                    loss_r = self.loss.eval_optimal_loss(y[~filt,:], x_lr[~filt,:])
                    loss = loss_l + loss_r
                    if loss < current_loss:
                        is_leaf = False
                        current_loss = np.copy(loss)
                        self.g.nodes[node]['variable'] = i
                        self.g.nodes[node]['level'] = q_i

            # if a split reducing the total loss has been found, split the node at best variable/level
            self.g.nodes[node]['leaf'] = is_leaf
            if is_leaf:
                self.g.nodes[node]['response'] = self.loss.eval_optimal_response(
                    y, x_lr) * learning_rate
            else:
                node_name = self._add_new_node(node)
                filt_1 = x[:, self.g.nodes[node]['variable']] < self.g.nodes[node]['level']
                filt_2 = x[:, self.g.nodes[node]['variable']] >= self.g.nodes[node]['level']
                if (np.sum(filt_1) == 0) or (np.sum(filt_2) == 0):
                    self.g.nodes[node]['response'] = self.loss.eval_optimal_response(
                        y, x_lr) * learning_rate
                    self.g.nodes[node]['leaf'] = True
                    return
                else:
                    self._fit_node_lr(x[filt_1, :], y[filt_1, :], x_lr[filt_1,:],
                                      node_name, learning_rate,hessian[filt_1,:])
                    self._fit_node_lr(x[filt_2, :], y[filt_2, :], x_lr[filt_2,:],
                                      node_name + 1, learning_rate,hessian[filt_2,:])

    # this funciton requires graphviz, which has installation issues under windows
    '''
    def plot_tree(self):
        colors = [self.g.nodes[n]['loss'] for n in self.g.nodes()]
        pos = graphviz_layout(self.g,prog='dot')
        nc = nx.draw_networkx_nodes(self.g, pos, node_size=30, alpha=0.8, node_color=colors, cmap=plt.cm.plasma)
        nx.draw_networkx_edges(self.g, pos, alpha=0.2)
        nx.draw_networkx_labels(self.g,pos, font_size=10)
        plt.colorbar(nc,label='attribute')
    '''

    def _find_regions(self, x, y):
        leaves = nx.get_node_attributes(self.g, 'response').keys()
        filters = []
        for leaf in leaves:
            path = nx.shortest_path(self.g, 'origin', leaf)
            node_filter = np.ones(x.shape[0], dtype=bool)
            for i in range(len(path)-1):
                v_p = self.g.nodes[path[i]]['variable']
                level_p = self.g.nodes[path[i]]['level']
                side = self.g.nodes[path[i+1]]['side']
                if side == 'l':
                    node_filter_i = x[:,v_p] < level_p
                else:
                    node_filter_i = x[:, v_p] >= level_p
                node_filter &= node_filter_i
            self.g.nodes[leaf]['y'] = y[node_filter,:]
            self.g.nodes[leaf]['filter'] = node_filter
            filters.append(node_filter)
            del node_filter
        return filters

    def _refit(self, x, y, discount=1.0, x_lr=None):
        if "exact_response" not in dir(self.loss):
            return
        self._find_regions(x, y)
        leaves = nx.get_node_attributes(self.g, 'response').keys()
        for l in leaves:
            y_l = self.g.nodes[l]['y']
            filter_l = self.g.nodes[l]['filter']
            if x_lr is None:
                self.g.nodes[l]['response'] = self.loss.exact_response(y_l) * discount
            else:
                x_lr_l = x_lr[filter_l,:]
                self.g.nodes[l]['response'] = self.loss.exact_response(y_l, x_lr_l) * discount
            self.g.nodes[l]['y'] = None

    def _add_new_node(self, node):
        node_name = self.g.number_of_nodes() + 1
        self.g.add_node(node_name, loss=np.inf, side='l')
        self.g.add_node(node_name + 1, loss=np.inf, side='r')
        self.g.add_edges_from(((node, node_name), (node, node_name + 1)))
        return node_name

    @staticmethod
    def _compute_bin_sums(y, edges, x):
        accumarray_idx = np.digitize(x,edges)
        partial_sums = np.zeros((len(edges)+1, y.shape[1]))
        for i in range(y.shape[1]):
            partial_sums[:,i] = np.bincount(accumarray_idx, weights=y[:,i]).ravel()

        accumulated_bin_sums = np.cumsum(partial_sums, axis=0)
        accumulated_elements = np.cumsum(np.bincount(accumarray_idx), axis=0)
        yy_left = accumulated_bin_sums ** 2
        yy_right = (accumulated_bin_sums[-1,:]-accumulated_bin_sums) ** 2
        return yy_left, yy_right, accumulated_elements

    def compute_loss(self, G2_left, G2_right, H_left, H_right, j):
        loss_l = self.loss.eval_optimal_loss(G2_left[j], H_left[j])
        loss_r = self.loss.eval_optimal_loss(G2_right[j], H_right[j])
        loss = loss_l + loss_r
        return loss

@jit('float64[:,:](float64[:,:], float64[:], float64[:])',
     nopython=True, parallel=False)
def bin_sums(y, edges, x):
    edges = edges
    y_ij = np.zeros((len(edges)-1,y.shape[1]))
    for i in np.arange(len(edges)-1):
        index = ((edges[i] <= x) * (x < edges[i+1]))
        y_ij[i,:] = np.sum(y[index,:],axis=0) + y_ij[i-1,:]
    return y_ij

@jit('Tuple([float64[:,:],float64[:,:]])(float64[:,:], float64[:], float64[:], int16)',
     nopython=True, parallel=False)
def leaf_stats(y, edges, x, order):
    y_ij = bin_sums(y, edges, x)
    if order == 1:
        s_left = y_ij ** 2
        s_right = (y_ij[-1,:]-y_ij) ** 2
    elif order == 2:
        s_left = y_ij
        s_right = (y_ij[-1, :] - y_ij)
    else:
        s_left, s_right = None, None
    return s_left, s_right

class MBT:
    """
    Multivariate Boosted Tree class. Fits a multivariate tree using boosting.

    :param n_boosts: maximum number of boosting rounds. Default: 20
    :param early_stopping_rounds: if the total loss is non-decreasing after early_stopping_rounds, stop training.
                                  The final model is the one which achieved the lowest loss up to the final iteration.
                                  Default: 3.
    :param learning_rate: in [0, 1]. A learning rate < 1 helps reducing overfitting. Default: 0.1.
    :param val_ratio: in [0,1]. If provided, the early stop is triggered by the loss computed on a validation set,
                      randomly extracted from the training set. The length of the validation set is
                      val_ratio * len (training set). Default: 0.
    :param val_patj: path to x_val and y_val folder used as validation dataset. If so , the loss is calculated on
                     this validation dataset at each round. Allows to define k-folds externally. No default.
    :param n_q: number of quantiles for the split search. Default: 10.
    :param min_leaf: minimum number of observations in one leaf. This parameter greatly affect generalization abilities.
                     Default: 100.
    :param loss_type: loss type for choosing the best splits. Currently the following losses are implemented:

                      mse: mean squared error loss, a.k.a. L2, ordinary least squares

                      time_smoother: mse with an additional penalization on the second order differences of the response
                      function. Requires to pass also lambda_smooth parameter.

                      latent_variable: generate response function of dimension n_t from an arbitrary linear combination
                      of n_r responses. This requires to pass also S and precision pars.

                      linear_regression: mse with linear response function. Using this loss function, when calling fit
                      and predict methods, one must also pass x_lr as additional argument, which is the matrix of
                      features used to train the linear response inside the leaf (which can be different from the
                      features used to grow the tree, x).

                      fourier: mse with linear response function, fitted on the first n_harmonics (where the
                      fundamental has wave-lenght equal to the target output). This requires to pass also the
                      n_harmonics parameter.

                      quantile: quantile loss function, a.k.a. pinball loss. This requires to pass also the alphas
                      parameter, a list of quantiles to be fitted.

                      quadratic_quantile: quadratic quantile loss function tailored for trees. It has a non-discontinuos
                      derivative. This requires to pass also the alphas parameter, a list of quantiles to be fitted.
    :param lambda_weights: coefficient for the quadratic regularization of the response's parameters. Default: 0.1
    :param lambda_leaves: coefficient for the quadratic regularization of the total number of leaves. This is only used
                          when the Tree is used as a weak learner by MBT. Default: 0.1
    :param verbose: in {0,1}. If set to 1, the MBT return fitting information at each iteration.
    :param refit: if True, if the loss function has an "exact_response" method, use it to refit the tree
    :param loss_kwargs: possible additional arguments for the loss function
    """
    def __init__(self, n_boosts: int = 20, early_stopping_rounds: int = 3, learning_rate: float = 0.1,
                 val_ratio: int = 0, val_path: str = '~/workspace/bluecloud descriptor/data',
                 n_q: int = 10, min_leaf: int = 100, loss_type: str = "mse",
                 lambda_weights: float = 0.1, lambda_leaves: float = 0.1, verbose: int = 0, refit=True, **loss_kwargs):
        self.n_boosts = n_boosts
        self.early_stopping_rounds = early_stopping_rounds
        self.learning_rate = learning_rate
        self.val_ratio = val_ratio
        self.val_path = val_path
        self.refit = refit
        self.tree_pars = {**{'n_q': n_q,'min_leaf': min_leaf, 'loss_type': loss_type, 'lambda_weights': lambda_weights,
                             'lambda_leaves': lambda_leaves}, **loss_kwargs}

        self.trees = []
        self.y_0 = None
        self.verbose = verbose

    def fit(self, x, y, do_plot=False, x_lr=None):
        """
         Fits an MBT using the features specified in the matrix :math:`x\\in\\mathbb{R}^{n_{obs} \\times n_{f}}`, in
         order to predict the targets in the matrix :math:`y\\in\\mathbb{R}^{n_{obs} \\times n_{t}}`,
         where :math:`n_{obs}` is the number of observations, :math:`n_{f}` the number of features and :math:`n_{t}`
         the dimension of the target.

        :param x: feature matrix, np.ndarray.
        :param y: target matrix, np.ndarray.
        :param x_lr: features for fitting the linear response inside the leaves. This is only required if a LinearLoss
                      is being used.
        """
        # divide in training and validation sets (has effect only if pars['val_ratio'] or ['val_id'] was set
        x_tr, x_val, y_tr, y_val, x_lr_tr, x_lr_val = self._validation_split(x, y, x_lr)
        lowest_loss = np.inf
        best_iter = 0
        self.trees = []
        t_init = time()
        loss_track = []
        if do_plot:
            fig, ax = set_figure((5, 4))

        for iter in tqdm(range(self.n_boosts)):
            t0 = time()
            tree = Tree(**self.tree_pars)
            if iter == 0:
                y_hat_val = self._fit_initial_guess(tree, y)
                y_hat = self._fit_initial_guess(tree, y)
                neg_grad, hessian = self._get_neg_grad_and_hessian_diags(tree, y_tr, y_hat, 0, x)

            # fit the tree
            tree.fit(x_tr, neg_grad, learning_rate=self.learning_rate, hessian=hessian, x_lr=x_lr_tr)

            # if loss function has an "exact_response" method, use it to refit the tree
            if self.refit:
                tree._refit(x_tr, y_tr - y_hat, self.learning_rate, x_lr=x_lr_tr)

            self.trees.append(tree)
            y_hat = y_hat + tree.predict(x_tr,x_lr_tr)
            y_hat_val = y_hat_val + tree.predict(x_val, x_lr_val)
            se, regularization = tree.loss.eval(y_val, y_hat_val, self.trees)
            terminal_leaves = len(nx.get_node_attributes(tree.g, 'response').values())
            if self.verbose == 1:
                tqdm.write('   Iteration: {} fitted in {:0.1e} sec, total loss: {:0.3e}, squared err: {:0.2e},  '
                      'regularization: {:0.2e}, best iter: {}, terminal leaves: {}'.format(iter, time()-t0,
                                                                                se + regularization, se,
                                                                                regularization, best_iter,
                                                                                terminal_leaves))
            loss = se + regularization
            loss_track.append(loss)

            if loss < lowest_loss:
                best_iter = iter
                lowest_loss = np.copy(loss)
            if iter - best_iter > self.early_stopping_rounds:
                break

            neg_grad, hessian = self._get_neg_grad_and_hessian_diags(tree, y_tr, y_hat, iter+1, x)
            if do_plot:
                if type(tree.loss) in [QuantileLoss, QuadraticQuantileLoss]:
                    ax.cla()
                    n_q = y_hat.shape[1]
                    n_plot = 200
                    colors = plt.get_cmap('plasma', int(n_q))
                    for fl in np.arange(np.floor(n_q / 2), dtype=int):
                        q_low = np.squeeze(y_hat[:n_plot, fl])
                        q_up = np.squeeze(y_hat[:n_plot, n_q - fl - 1])
                        x_plot= np.arange(len(q_low))
                        ax.fill_between(x_plot, q_low, q_up, color=colors(fl), alpha=0.1 + 0.6 * fl / n_q,
                                        linewidth=0.0)

                    ax.plot(y[:n_plot, :], linewidth=2, label='target')
                    ax.legend(loc='upper right')
                    plt.title('Quantiles on first {} samples'.format(n_plot))
                else:
                    ax.cla()
                    ax.plot([np.min(y_tr[:, 0]), np.max(y_tr[:, 0])], [np.min(y_tr[:, 0]), np.max(y_tr[:, 0])], '--')
                    ax.scatter(y_tr[:, 0], y_hat[:, 0], marker='.', alpha=0.2)
                    ax.set_xlabel('observations')
                    ax.set_ylabel('predictions')
                    ax.set_title('Fit on first y dimension')

                plt.pause(0.1)

        print('#---------------- Model fitted in {:0.2e} min ----------------'.format((time()-t_init)/60))

        # keep trees up to best iteration
        self.trees = self.trees[:best_iter+1]
        if do_plot:
            plt.close(fig)

        return self, loss_track

    def predict(self, x, n=None, x_lr=None):
        """
        Predicts the target based on the feature matrix x (and linear regression features x_lr).

        :param x: feature matrix, np.ndarray.
        :param n: predict up to the nth fitted tree. If None, predict all the trees. Default: None
        :param x_lr: linear regression feature matrix, np.ndarray. Only required if LinearLoss has been used.

        :return: target's predictions
        """
        y_hat = self.y_0
        if n is None:
            for tree in self.trees:
                y_hat = y_hat + tree.predict(x, x_lr)
        else:
            for tree in self.trees[:n]:
                y_hat = y_hat + tree.predict(x, x_lr)
        return y_hat

    def _get_neg_grad_and_hessian_diags(self, tree, y: np.ndarray, y_hat: np.ndarray, iteration: int, x: np.ndarray):
        """
        Returns the negative gradient at current iteration for all the observations, and a matrix which ith row is
        the diagonal of the Hessian for the current observation. This matrix is later used by the loss functions to
        obtain the inverse of the Hessian in each leaf, with proper dimensions.

        :param tree: current fitted tree
        :param y: targets
        :param y_hat: predictions at current iteration
        :param iteration: current iteration
        :param x: feature treaning set (for finding tree regions, if needed by the loss,
        as for the quadratic_quantile loss)

        :return: negative gradient and Hessian's diags
        """
        if iteration == 0:
            leaves_idx = [np.ones(len(y), dtype=bool)]
        else:
            leaves_idx = tree._find_regions(x, y)

        grad, hessian_diags = tree.loss.get_grad_and_hessian_diags(y, y_hat, iteration, leaves_idx)

        return - grad, hessian_diags

    def _fit_initial_guess(self, tree, y) -> np.ndarray:
        y_0 = tree.loss.get_initial_guess(y)
        self.y_0 = y_0
        return y_0

    def _validation_split(self,x,y,x_lr):
        if self.val_path is not None:
            x_val = read_feather(self.val_path+'_X_val.feather').to_numpy()
            y_val = read_feather(self.val_path+'_Y_val.feather').to_numpy()
            x_lr_val = read_feather(self.val_path+'_X_lr_val.feather').to_numpy() if x_lr is not None else None
            x_tr, x_val = [x, x_val]
            y_tr, y_val = [y, y_val]
            x_lr_tr, x_lr_val = [x, x_lr_val] if x_lr is not None else [None, None]
        elif self.val_ratio == 0:
            x_tr = x_val = x
            y_tr = y_val = y
            x_lr_tr = x_lr_val = x_lr
        else:
            randidx = np.random.permutation(x.shape[0])
            n_val = int(self.val_ratio*x.shape[0])
            tr_idx, val_idx = [randidx[:-n_val],randidx[-n_val:]]
            x_tr, x_val = [x[tr_idx],x[val_idx]]
            y_tr, y_val = [y[tr_idx], y[val_idx]]
            x_lr_tr, x_lr_val = [x_lr[tr_idx], x_lr[val_idx]] if x_lr is not None else [None, None]

        return x_tr, x_val, y_tr, y_val, x_lr_tr, x_lr_val

if __name__ == '__main__':
    pass

ext: py

max_stars_repo_name: andrea-manzi/mbtr

max_issues_repo_name: andrea-manzi/mbtr

max_forks_repo_name: andrea-manzi/mbtr

reason: import numpy,from numba,import networkx,from networkx


text:
# -*- coding: utf-8 -*-
"""MMBT.ipynb

This py script trains and tests a Multi-Modal Bit Transformer from the Simple Transformers Library
"""

import pandas as pd
import numpy as np
from sklearn.metrics import matthews_corrcoef, confusion_matrix, precision_score, recall_score, f1_score, accuracy_score
from simpletransformers.classification.multi_modal_classification_model import MultiModalClassificationModel
from ekphrasis.classes.preprocessor import TextPreProcessor
from ekphrasis.classes.tokenizer import SocialTokenizer
from ekphrasis.dicts.emoticons import emoticons

#Create DataFrame with three columns text, image, labels

def combine_text(df):
    """
    Combines tweet and image text into one column

    df: Dataframe which holds the data
    """
    combined_text = []

    for row_num in range(len(df)):
        tweet_text = df.loc[row_num, "tweet_text"]
        image_text = df.loc[row_num, "img_text"]
        if type(image_text) == str:
            combined_text.append(tweet_text + image_text)
        else:
            combined_text.append(tweet_text)

    return combined_text

def clean_text(data, normalize_list, annotate_list):
    """
    This function preprocesses the text using the Ekphrasis library

    data: Pandas series object containing strings of text

    normalize_list: list of data features to clean

    annotate_list: list of data features to annotate
    """

    text_processor = TextPreProcessor(
        normalize= normalize_list,
        annotate= annotate_list,
        fix_html=True,
        segmenter="twitter",
        unpack_hashtags=True,
        unpack_contractions=True,
        spell_correct_elong=True,
        tokenizer=SocialTokenizer(lowercase=True).tokenize,
        dicts=[emoticons]
    )

    clean_data = data.map(lambda x: " ".join(text_processor.pre_process_doc(x)))

    return clean_data

def df_preparer(file_path):
    """
    file_path:
    """

    df = pd.read_csv(file_path, encoding='utf-8')
    df["text"] = combine_text(df)
    df["text"] = clean_text(df["text"], ['url', 'email', 'percent', 'money', 'phone', 'user', 'time', 'date', 'number'], ['hashtag', 'allcaps', 'elongated', 'repeated', 'emphasis', 'censored'])
    df["labels"] = df["Primary_numeric_gt"].map(lambda x: str(x)) #labels need to be strings
    df["images"] = df["image_number"].map(lambda x: str(x) + ".jpg")

    df = df[["labels", "text", "images"]]

    return df

def metrics(labels, preds, argmax_needed: bool = False):
    """
    Returns the Matthew's correlation coefficient, accuracy rate, true positive rate, true negative rate, false positive rate, false negative rate, precission, recall, and f1 score

    labels: list of correct labels

    pred: list of model predictions

    argmax_needed (boolean): converts logits to predictions. Defaulted to false.
    """
    labels = labels
    preds = preds

    if argmax_needed == True:
        preds = np.argmax(preds, axis=1).flatten()

    mcc = matthews_corrcoef(labels, preds)
    acc = accuracy_score(labels, preds)
    cm = confusion_matrix(labels, preds)

    f1 = f1_score(labels, preds, average= "weighted")
    precision = precision_score(labels, preds, average= "weighted")
    recall = f1_score(labels, preds, average= "weighted")

    results = {
        "mcc": mcc,
        "acc": acc,
        "confusion_matrix": cm,
        "precision": precision,
        "recall": recall,
        "f1": f1,
    }

    return results, labels, preds

# Prepare data
train = df_preparer("/content/drive/My Drive/Dog_Whistle_Code/Data/Train/dog_whistle_train.csv")
dev = df_preparer("/content/drive/My Drive/Dog_Whistle_Code/Data/Validation/dog_whistle_dev.csv")
test = df_preparer("/content/drive/My Drive/Dog_Whistle_Code/Data/Test/dog_whistle_test.csv")
test["int_labels"] = test["labels"].map(lambda x: int(x))

# Run Model
args = {"fp16": False, 'train_batch_size': 8, 'eval_batch_size': 8, 'learning_rate': 1e-5, 'num_train_epochs': 4, 'max_seq_length': 100, "weight_decay": 0.1, "warmup_ratio": 0.06,
        "use_early_stopping": True, "early_stopping_patience": 2, "early_stopping_delta": 0, "early_stopping_metric": "eval_loss", "overwrite_output_dir": True}

model = MultiModalClassificationModel("bert", "bert-large-uncased", label_list= ["0", "1", "2", "3"], args=args, use_cuda= True)
print("Loading finished")
model.train_model(train, eval_data=dev, image_path="/content/drive/My Drive/Dog_Whistle_Code/Data/Images/", output_dir= "/content/drive/My Drive/Dog_Whistle_Code/Fine_Tuned_Models/Multimodal/MMBT")
print("Training finished")
results, logits  = model.eval_model(test, image_path="/content/drive/My Drive/Dog_Whistle_Code/Data/Images/")
metric_vals, labels, preds = metrics(test["int_labels"], logits, argmax_needed=True)
print(metric_vals)

#Save
np.save("/content/drive/My Drive/Dog_Whistle_Code/Fine_Tuned_Models/Multimodal/MMBT/results.npy", metric_vals)

df = pd.DataFrame([labels, preds]).T
df = df.rename(columns ={"int_labels": "Labels", "Unnamed 0": "Predictions"})
df.to_csv("/content/drive/My Drive/Dog_Whistle_Code/Fine_Tuned_Models/Multimodal/MMBT/predictions.csv")

reason: import numpy


text:
import random
import collections

import numpy as np
import tensorflow as tf

class Agent:
    BUFFER_SIZE = 3000
    BATCH_SIZE = 200
    STATE_SIZE = 12
    DENSE_UNITS = 150
    N_HIDDEN_LAYERS = 3
    MOVE_OPTIONS = ["up", "down", "left", "right"]
    GAMMA = 0.9

    def __init__(self):
        self.epsilon = 0.5
        self.num_iterations = 0
        self.n_outputs = len(self.MOVE_OPTIONS)

        self.state_history = []

        self.state_buffer = collections.deque(maxlen=self.BUFFER_SIZE)
        self.reward_buffer = collections.deque(maxlen=self.BUFFER_SIZE)
        self.target_buffer = collections.deque(maxlen=self.BUFFER_SIZE)

        self.network = self.create_network()

    @staticmethod
    def get_state(game, output_format=None):
        state = dict()
        state["raspi_up"] = game.snake.y[0] > game.raspi.y
        state["raspi_down"] = game.snake.y[0] < game.raspi.y
        state["raspi_right"] = game.snake.x[0] < game.raspi.x
        state["raspi_left"] = game.snake.x[0] > game.raspi.x

        wall_up = (game.snake.y[0] - game.walls.y_min_inside) <= game.SNAKE_STEP
        wall_down = (game.walls.y_max_inside - game.snake.y[0]) <= game.SNAKE_STEP
        wall_left = (game.snake.x[0] - game.walls.x_min_inside) <= game.SNAKE_STEP
        wall_right = (game.walls.x_max_inside - game.snake.x[0]) <= game.SNAKE_STEP

        state["obstacle_up"] = wall_up or game.snake.is_body_up()
        state["obstacle_down"] = wall_down or game.snake.is_body_down()
        state["obstacle_left"] = wall_left or game.snake.is_body_left()
        state["obstacle_right"] = wall_right or game.snake.is_body_right()

        state["snake_moving_up"] = game.snake.direction == "up"
        state["snake_moving_down"] = game.snake.direction == "down"
        state["snake_moving_right"] = game.snake.direction == "right"
        state["snake_moving_left"] = game.snake.direction == "left"

        if output_format == "network":
            state = [float(val) for val in state.values()]
            state = np.array(state)[np.newaxis, :]

        return state

    @staticmethod
    def compute_reward(game):
        reward = (
            float(game.snake.has_lost(game.walls)) * (-20)
            + float(game.snake.has_eaten(game.raspi.x, game.raspi.y)) * 10
        )
        return reward

    def compute_target(self):
        if self.num_iterations > 0:
            q_function_last_state = self.network.predict(self.state_buffer[-2])

            q_function_current_state = self.network.predict(self.state_buffer[-1])

            last_reward = self.reward_buffer[-1]

            target_last_action = last_reward + self.GAMMA * np.max(
                q_function_current_state
            )
            global_target = q_function_last_state
            global_target[0, np.argmax(q_function_last_state)] = target_last_action
        else:
            global_target = np.zeros([1, self.n_outputs])

        return global_target

    def create_network(self):
        input_layer = tf.keras.layers.Input(shape=(self.STATE_SIZE,))
        layer = input_layer

        for i_layer in range(self.N_HIDDEN_LAYERS):
            layer = tf.keras.layers.Dense(self.DENSE_UNITS, activation="relu")(layer)

        output = tf.keras.layers.Dense(
            self.n_outputs, activation=tf.keras.activations.linear
        )(layer)

        network = tf.keras.Model(inputs=input_layer, outputs=output)

        network.compile(
            loss=tf.keras.losses.MeanSquaredError(),
            optimizer=tf.keras.optimizers.Adam(learning_rate=0.0005),
        )

        return network

    def train_network(self):
        n_epochs = 10
        if self.num_iterations > self.BATCH_SIZE:

            state_batch = np.squeeze(self.state_buffer)
            current_state_batch = state_batch[0:-1, :]

            reward_batch = np.array(self.reward_buffer)
            reward_batch = reward_batch[1:]

            next_state_batch = state_batch[1:, :]

            for i_epoch in range(n_epochs):
                target_batch = self.network.predict(current_state_batch)
                idx_action = np.argmax(target_batch, axis=1)
                target_batch[np.arange(target_batch.shape[0]), idx_action] = (
                    reward_batch
                    + self.GAMMA
                    * np.max(self.network.predict(next_state_batch), axis=1)
                )

                dataset = tf.data.Dataset.from_tensor_slices(
                    (current_state_batch, target_batch)
                )
                dataset = (
                    dataset.shuffle(200)
                    .prefetch(buffer_size=self.BATCH_SIZE)
                    .batch(self.BATCH_SIZE)
                )

                self.network.fit(dataset, verbose=False, epochs=1)

    def train_network_short_term(self):
        if self.num_iterations > 1:
            last_state = self.state_buffer[-2]
            target = self.target_buffer[-1]

            self.network.fit(x=last_state, y=target, verbose=False, epochs=1)

    def next_move(self, game):

        if random.random() > self.epsilon:
            state = self.get_state(game, output_format="network")
            q_function_output = self.network.predict(state)
            idx_max = np.argmax(q_function_output.ravel())
        else:
            idx_max = random.randint(0, self.n_outputs - 1)

        self.decrease_epsilon()
        next_move = self.MOVE_OPTIONS[idx_max]

        return next_move

    def store(self, game):
        state = Agent.get_state(game, output_format="network")
        self.state_history.append(state)
        self.state_buffer.append(state)

        reward = Agent.compute_reward(game)
        self.reward_buffer.append(reward)

        target = self.compute_target()
        self.target_buffer.append(target)

        self.num_iterations += 1

    def decrease_epsilon(self):
        self.epsilon *= 0.9995

    def summary(self, game):
        print(
            f"Game N°: {game.ind_game}, epsilon: {self.epsilon:1.2f}, "
            f"score: {game.score}, high score: {game.high_score}"
        )

max_stars_repo_name: DTicien/snake-game

max_issues_repo_name: DTicien/snake-game

max_forks_repo_name: DTicien/snake-game


text:
# -------------------------------------------------------------------------------
# This file is part of jHDF. A pure Java library for accessing HDF5 files.
#
# http://jhdf.io
#
# Copyright (c) 2022 James Mudd
#
# MIT License see 'LICENSE' file
# -------------------------------------------------------------------------------
import h5py

import numpy as np

# The idea of this test if to write variable length (ragged array) datasets
def write_vlen_datasets(f):
    # Unsigned int
    uint8_vlen_type = h5py.vlen_dtype(np.uint8)
    uint8_vlen_dataset = f.create_dataset("vlen_uint8_data", (3,), dtype=uint8_vlen_type)
    uint8_vlen_dataset[0] = [0]
    uint8_vlen_dataset[1] = [1, 2]
    uint8_vlen_dataset[2] = [3, 4, 5]

    uint16_vlen_type_chunked = h5py.vlen_dtype(np.uint16)
    uint16_vlen_dataset = f.create_dataset("vlen_uint16_data", (3,), dtype=uint16_vlen_type_chunked)
    uint16_vlen_dataset[0] = [0]
    uint16_vlen_dataset[1] = [1, 2]
    uint16_vlen_dataset[2] = [3, 4, 5]

    uint32_vlen_type = h5py.vlen_dtype(np.uint32)
    uint32_vlen_dataset = f.create_dataset("vlen_uint32_data", (3,), dtype=uint32_vlen_type)
    uint32_vlen_dataset[0] = [0]
    uint32_vlen_dataset[1] = [1, 2]
    uint32_vlen_dataset[2] = [3, 4, 5]

    uint64_vlen_type = h5py.vlen_dtype(np.uint64)
    uint64_vlen_dataset = f.create_dataset("vlen_uint64_data", (3,), dtype=uint64_vlen_type)
    uint64_vlen_dataset[0] = [0]
    uint64_vlen_dataset[1] = [1, 2]
    uint64_vlen_dataset[2] = [3, 4, 5]

    # Signed int
    int8_vlen_type = h5py.vlen_dtype(np.int8)
    int8_vlen_dataset = f.create_dataset("vlen_int8_data", (3,), dtype=int8_vlen_type)
    int8_vlen_dataset[0] = [0]
    int8_vlen_dataset[1] = [1, 2]
    int8_vlen_dataset[2] = [3, 4, 5]

    int16_vlen_type_chunked = h5py.vlen_dtype(np.int16)
    int16_vlen_dataset = f.create_dataset("vlen_int16_data", (3,), dtype=int16_vlen_type_chunked)
    int16_vlen_dataset[0] = [0]
    int16_vlen_dataset[1] = [1, 2]
    int16_vlen_dataset[2] = [3, 4, 5]

    int32_vlen_type = h5py.vlen_dtype(np.int32)
    int32_vlen_dataset = f.create_dataset("vlen_int32_data", (3,), dtype=int32_vlen_type)
    int32_vlen_dataset[0] = [0]
    int32_vlen_dataset[1] = [1, 2]
    int32_vlen_dataset[2] = [3, 4, 5]

    int64_vlen_type = h5py.vlen_dtype(np.int64)
    int64_vlen_dataset = f.create_dataset("vlen_int64_data", (3,), dtype=int64_vlen_type)
    int64_vlen_dataset[0] = [0]
    int64_vlen_dataset[1] = [1, 2]
    int64_vlen_dataset[2] = [3, 4, 5]

    # Floating point
    float32_vlen_type = h5py.vlen_dtype(np.float32)
    float32_vlen_dataset = f.create_dataset("vlen_float32_data", (3,), dtype=float32_vlen_type)
    float32_vlen_dataset[0] = [0]
    float32_vlen_dataset[1] = [1, 2]
    float32_vlen_dataset[2] = [3, 4, 5]

    float64_vlen_type = h5py.vlen_dtype(np.float64)
    float64_vlen_dataset = f.create_dataset("vlen_float64_data", (3,), dtype=float64_vlen_type)
    float64_vlen_dataset[0] = [0]
    float64_vlen_dataset[1] = [1, 2]
    float64_vlen_dataset[2] = [3, 4, 5]

    # https://github.com/jamesmudd/jhdf/issues/247
    int32_vlen_type = h5py.vlen_dtype(np.dtype(np.int32))
    int32_vlen_dataset = f.create_dataset('vlen_issue_247', (3,), dtype=int32_vlen_type)
    int32_vlen_dataset[0] = [1, 2, 3]
    int32_vlen_dataset[1] = []
    int32_vlen_dataset[2] = [1, 2, 3, 4, 5]

    # Chunked
    # Unsigned int
    uint8_vlen_type = h5py.vlen_dtype(np.uint8)
    uint8_vlen_dataset_chunked = f.create_dataset("vlen_uint8_data_chunked", (3,), dtype=uint8_vlen_type, chunks=(3,))
    uint8_vlen_dataset_chunked[0] = [0]
    uint8_vlen_dataset_chunked[1] = [1, 2]
    uint8_vlen_dataset_chunked[2] = [3, 4, 5]

    uint16_vlen_type_chunked = h5py.vlen_dtype(np.uint16)
    uint16_vlen_dataset = f.create_dataset("vlen_uint16_data_chunked", (3,), dtype=uint16_vlen_type_chunked,
                                           chunks=(3,))
    uint16_vlen_dataset[0] = [0]
    uint16_vlen_dataset[1] = [1, 2]
    uint16_vlen_dataset[2] = [3, 4, 5]

    uint32_vlen_type = h5py.vlen_dtype(np.uint32)
    uint32_vlen_dataset_chunked = f.create_dataset("vlen_uint32_data_chunked", (3,), dtype=uint32_vlen_type,
                                                   chunks=(3,))
    uint32_vlen_dataset_chunked[0] = [0]
    uint32_vlen_dataset_chunked[1] = [1, 2]
    uint32_vlen_dataset_chunked[2] = [3, 4, 5]

    uint64_vlen_type = h5py.vlen_dtype(np.uint64)
    uint64_vlen_dataset_chunked = f.create_dataset("vlen_uint64_data_chunked", (3,), dtype=uint64_vlen_type,
                                                   chunks=(3,))
    uint64_vlen_dataset_chunked[0] = [0]
    uint64_vlen_dataset_chunked[1] = [1, 2]
    uint64_vlen_dataset_chunked[2] = [3, 4, 5]

    # Signed int
    int8_vlen_type = h5py.vlen_dtype(np.int8)
    int8_vlen_dataset = f.create_dataset("vlen_int8_data_chunked", (3,), dtype=int8_vlen_type, chunks=(3,))
    int8_vlen_dataset[0] = [0]
    int8_vlen_dataset[1] = [1, 2]
    int8_vlen_dataset[2] = [3, 4, 5]

    int16_vlen_type_chunked = h5py.vlen_dtype(np.int16)
    int16_vlen_dataset = f.create_dataset("vlen_int16_data_chunked", (3,), dtype=int16_vlen_type_chunked, chunks=(3,))
    int16_vlen_dataset[0] = [0]
    int16_vlen_dataset[1] = [1, 2]
    int16_vlen_dataset[2] = [3, 4, 5]

    int32_vlen_type = h5py.vlen_dtype(np.int32)
    int32_vlen_dataset = f.create_dataset("vlen_int32_data_chunked", (3,), dtype=int32_vlen_type, chunks=(3,))
    int32_vlen_dataset[0] = [0]
    int32_vlen_dataset[1] = [1, 2]
    int32_vlen_dataset[2] = [3, 4, 5]

    int64_vlen_type = h5py.vlen_dtype(np.int64)
    int64_vlen_dataset = f.create_dataset("vlen_int64_data_chunked", (3,), dtype=int64_vlen_type, chunks=(3,))
    int64_vlen_dataset[0] = [0]
    int64_vlen_dataset[1] = [1, 2]
    int64_vlen_dataset[2] = [3, 4, 5]

    # Floating point
    float32_vlen_type = h5py.vlen_dtype(np.float32)
    float32_vlen_dataset_chunked = f.create_dataset("vlen_float32_data_chunked", (3,), dtype=float32_vlen_type,
                                                    chunks=(3,))
    float32_vlen_dataset_chunked[0] = [0]
    float32_vlen_dataset_chunked[1] = [1, 2]
    float32_vlen_dataset_chunked[2] = [3, 4, 5]

    float64_vlen_type = h5py.vlen_dtype(np.float64)
    float64_vlen_dataset_chunked = f.create_dataset("vlen_float64_data_chunked", (3,), dtype=float64_vlen_type,
                                                    chunks=(3,))
    float64_vlen_dataset_chunked[0] = [0]
    float64_vlen_dataset_chunked[1] = [1, 2]
    float64_vlen_dataset_chunked[2] = [3, 4, 5]

    # https://github.com/jamesmudd/jhdf/issues/247
    int32_vlen_type = h5py.vlen_dtype(np.dtype(np.int32))
    int32_vlen_dataset = f.create_dataset('vlen_issue_247_chunked', (3,), dtype=int32_vlen_type, chunks=(3,))
    int32_vlen_dataset[0] = [1, 2, 3]
    int32_vlen_dataset[1] = []
    int32_vlen_dataset[2] = [1, 2, 3, 4, 5]

    f.flush()
    f.close()

if __name__ == '__main__':
    print('Making variable length dataset test files...')

    f = h5py.File('test_vlen_datasets_earliest.hdf5', 'w', libver='earliest')
    write_vlen_datasets(f)
    print('created test_vlen_datasets_earliest.hdf5')

    f = h5py.File('test_vlen_datasets_latest.hdf5', 'w', libver='latest')
    write_vlen_datasets(f)
    print('created test_vlen_datasets_latest.hdf5')

max_stars_repo_name: raedma/jhdf

max_issues_repo_name: raedma/jhdf

max_forks_repo_name: raedma/jhdf


text:
from numpy.lib.function_base import median
import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread("crestamento_4.jpg")

image_RGB = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
image_gray = cv2.cvtColor(image_RGB, cv2.COLOR_RGB2GRAY)

median_blur = cv2.medianBlur(image_gray,9)

sobelx = cv2.Sobel(median_blur, cv2.CV_16S,1,0,ksize=-1)
sobelx_16S = np.absolute(sobelx)
sobelx_8U = np.uint8(sobelx_16S)

sobely = cv2.Sobel(median_blur, cv2.CV_16S,0,1,ksize=-1)
sobely_16S = np.absolute(sobely)
sobely_8U = np.uint8(sobely_16S)

result = np.sqrt(sobely**2 + sobelx**2)

# cálculo da quantidade de pixels da folha
_, imgBinaryLeaf = cv2.threshold(image_gray, 127, 255, cv2.THRESH_BINARY)

leaf_area = np.sum(imgBinaryLeaf == 0) #área da folha em pixels
# fim do cálculo da quantidade de pixels da folha

# cálculo da quantidade de pixels de borda
_, imgBinary = cv2.threshold(result, 127, 255, cv2.THRESH_BINARY)

border_area = np.sum(imgBinary == 255)  #área da borda em pixels

print('Border pixels: ', border_area)
print('Leaf area: ', leaf_area)
# fim do cálculo da quantidade de pixels de borda

#plt.imshow(imgBinary, cmap='gray')

cv2.imwrite('crestamento_sobel_bi.png', imgBinary)

plt.show()

max_stars_repo_name: joaocarlos-dev/filters-PDI

max_issues_repo_name: joaocarlos-dev/filters-PDI

max_forks_repo_name: joaocarlos-dev/filters-PDI

reason: import numpy,from numpy


text:
import settings
import model
from utils import VOC
import os
import time
import cv2
import matplotlib.pyplot as plt
import tensorflow as tf
import numpy as np

config = tf.ConfigProto()
config.gpu_options.visible_device_list = "0"

sess = tf.InteractiveSession(config = config)
model = model.Model(training = False)
sess.run(tf.global_variables_initializer())
saver = tf.train.Saver(tf.global_variables())
boundary1 = settings.cell_size * settings.cell_size * settings.num_class
boundary2 = boundary1 + settings.cell_size * settings.cell_size * settings.box_per_cell

tryOriginal = settings.try_original

try:
    if tryOriginal :
        saver.restore(sess, os.getcwd() + '/YOLO_small.ckpt')
        print('load from YOLO small pretrained')
    else :
        saver.restore(sess, os.getcwd() + '/model.ckpt')
        print('load from past checkpoint')
except:
    try:
        saver.restore(sess, os.getcwd() + '/YOLO_small.ckpt')
        print('load from YOLO small pretrained')
    except:
        print('you must train first, exiting..')
        exit(0)

def draw_result(img, result,color=(0, 255, 0)):
    for i in range(len(result)):
        x = int(result[i][1])
        y = int(result[i][2])
        w = int(result[i][3] / 2)
        h = int(result[i][4] / 2)
        cv2.rectangle(img, (x - w, y - h), (x + w, y + h), color, 2)
        cv2.rectangle(img, (x - w, y - h - 20), (x + w, y - h), (125, 125, 125), -1)
        cv2.putText(img, result[i][0] + ' : %.2f' % result[i][5], (x - w + 5, y - h - 7), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 0, 0), 1, cv2.LINE_AA)

def detect(img):
    img_h, img_w, _ = img.shape
    inputs = cv2.resize(img, (settings.image_size, settings.image_size))
    inputs = cv2.cvtColor(inputs, cv2.COLOR_BGR2RGB).astype(np.float32)
    inputs = (inputs / 255.0) * 2.0 - 1.0
    inputs = np.reshape(inputs, (1, settings.image_size, settings.image_size, 3))

    r1,r2 = detect_from_cvmat(inputs)
    result = r1[0]#detect_from_cvmat(inputs)[0]
    result2 = r2[0]
    #print(result)

    for i in range(len(result)):
        result[i][1] *= (1.0 * img_w / settings.image_size)
        result[i][2] *= (1.0 * img_h / settings.image_size)
        result[i][3] *= (1.0 * img_w / settings.image_size)
        result[i][4] *= (1.0 * img_h / settings.image_size)

        result2[i][0] *= (1.0 * img_w / settings.image_size)
        result2[i][1] *= (1.0 * img_h / settings.image_size)
        result2[i][2] *= (1.0 * img_w / settings.image_size)
        result2[i][3] *= (1.0 * img_h / settings.image_size)

    return result,result2

def detect_from_cvmat(inputs):
    net_output = sess.run(model.logits, feed_dict = {model.images: inputs})
    results = []
    results2 = []
    for i in range(net_output.shape[0]):
        r1,r2 = interpret_output(net_output[i])
        results.append(r1)
        results2.append(r2)
        #print(' -',i)
    '''print('dfc : ',results)
    print('dfc2 : ',results2)'''
    return results,results2

def iou(box1, box2):
    #tb=min(b1[x2
    tb = min(box1[0] + 0.5 * box1[2], box2[0] + 0.5 * box2[2]) - max(box1[0] - 0.5 * box1[2], box2[0] - 0.5 * box2[2])
    lr = min(box1[1] + 0.5 * box1[3], box2[1] + 0.5 * box2[3]) - max(box1[1] - 0.5 * box1[3], box2[1] - 0.5 * box2[3])
    if tb < 0 or lr < 0:
        intersection = 0
    else:
        intersection = tb * lr
    return intersection / (box1[2] * box1[3] + box2[2] * box2[3] - intersection)
    #return [intersection, (box1[2] * box1[3] + box2[2] * box2[3] - intersection)]

def iou_mine(box1, box2):
    tb = min(box1[0] + 0.5 * box1[2], box2[0] + 0.5 * box2[2]) - max(box1[0] - 0.5 * box1[2], box2[0] - 0.5 * box2[2])
    lr = min(box1[1] + 0.5 * box1[3], box2[1] + 0.5 * box2[3]) - max(box1[1] - 0.5 * box1[3], box2[1] - 0.5 * box2[3])
    if tb < 0 or lr < 0:
        intersection = 0
    else:
        intersection = tb * lr
    print(intersection,box1[2] * box1[3])
    return [intersection/(box1[2] * box1[3]) , (box1[2] * box1[3] + box2[2] * box2[3] - intersection)]

def interpret_output(output):
    probs = np.zeros((settings.cell_size, settings.cell_size, settings.box_per_cell, len(settings.classes_name)))
    class_probs = np.reshape(output[0 : boundary1], (settings.cell_size, settings.cell_size, settings.num_class))
    scales = np.reshape(output[boundary1 : boundary2], (settings.cell_size, settings.cell_size, settings.box_per_cell))
    boxes = np.reshape(output[boundary2 :], (settings.cell_size, settings.cell_size, settings.box_per_cell, 4))
    offset = np.transpose(np.reshape(np.array([np.arange(settings.cell_size)] * settings.cell_size * settings.box_per_cell), [settings.box_per_cell, settings.cell_size, settings.cell_size]), (1, 2, 0))

    boxes[:, :, :, 0] += offset
    boxes[:, :, :, 1] += np.transpose(offset, (1, 0, 2))
    boxes[:, :, :, :2] = 1.0 * boxes[:, :, :, 0:2] / settings.cell_size
    boxes[:, :, :, 2:] = np.square(boxes[:, :, :, 2:])

    boxes *= settings.image_size

    for i in range(settings.box_per_cell):
        for j in range(settings.num_class):
            probs[:, :, i, j] = np.multiply(class_probs[:, :, j], scales[:, :, i])

    filter_mat_probs = np.array(probs >= settings.threshold, dtype = 'bool')
    filter_mat_boxes = np.nonzero(filter_mat_probs)
    boxes_filtered = boxes[filter_mat_boxes[0], filter_mat_boxes[1], filter_mat_boxes[2]]
    probs_filtered = probs[filter_mat_probs]
    classes_num_filtered = np.argmax(filter_mat_probs, axis = 3)[filter_mat_boxes[0], filter_mat_boxes[1], filter_mat_boxes[2]]

    argsort = np.array(np.argsort(probs_filtered))[::-1]
    boxes_filtered = boxes_filtered[argsort]
    probs_filtered = probs_filtered[argsort]
    classes_num_filtered = classes_num_filtered[argsort]

    for i in range(len(boxes_filtered)):
        if probs_filtered[i] == 0:
            continue
        for j in range(i + 1, len(boxes_filtered)):
            if iou(boxes_filtered[i], boxes_filtered[j]) > settings.IOU_threshold:
                probs_filtered[j] = 0.0

    filter_iou = np.array(probs_filtered > 0.0, dtype = 'bool')
    boxes_filtered = boxes_filtered[filter_iou]
    probs_filtered = probs_filtered[filter_iou]
    classes_num_filtered = classes_num_filtered[filter_iou]

    result = []
    result2 = []
    for i in range(len(boxes_filtered)):
        result.append([settings.classes_name[classes_num_filtered[i]], boxes_filtered[i][0], boxes_filtered[i][1], boxes_filtered[i][2], boxes_filtered[i][3], probs_filtered[i]])
        result2.append([boxes_filtered[i][0], boxes_filtered[i][1], boxes_filtered[i][2], boxes_filtered[i][3],classes_num_filtered[i]])

    #print("The result : ",result)
    #print("The result2 : ",result2)
    return result,result2

def read_image(image, name,list_bb):
    result,result2 = detect(image)

    print(result,'\n',result2)
    #print('list_bb : ',list_bb[0])

    img_h, img_w, _ = image.shape

    list_bb_compare = []
    list_bb_compare2 = []

    for i in range(0,len(list_bb)) :
        #first convert the l_bb to the real image size

        list_bb[i][0] *= (1.0 * img_w / settings.image_size)
        list_bb[i][1] *= (1.0 * img_h / settings.image_size)
        list_bb[i][2] *= (1.0 * img_w / settings.image_size)
        list_bb[i][3] *= (1.0 * img_h / settings.image_size)

        #Now convert to the x,y,w,h format
        x = (list_bb[i][0] + list_bb[i][2])/2
        y = (list_bb[i][1] + list_bb[i][3])/2
        w = list_bb[i][2] - list_bb[i][0] ;
        h = list_bb[i][3] - list_bb[i][1] ;

        #Now add them
        list_bb_compare.append([x,y,w,h,list_bb[i][4]])
        list_bb_compare2.append(['test',x,y,w,h,0])

    print("lbb compare : ",list_bb_compare)

    draw_result(image, result,(0,255,0))
    draw_result(image, list_bb_compare2,(255,0,0))

    num_objects = len(list_bb_compare)

    num_cor = 0;
    num_false = 0;

    #classes_name =  ["van","truck", "bus", "car", "bicycle", "motorbike",'cyclist', "person"]
    #Now calculate the iou
    for i in range(len(list_bb_compare)) :
        if np.max(np.asarray([list_bb_compare[i][0],list_bb_compare[i][1],list_bb_compare[i][2],list_bb_compare[i][3]])) <= 0 :
            print("Skipping zero")
            continue
        else :
            print("Checking : ",list_bb_compare[i])

        t_threshold = .7
        if list_bb_compare[i][4] == 6 :
            t_threshold = .5
        elif  list_bb_compare[i][4]== 7 :
            t_threshold = .5

        bb_gt = [list_bb_compare[i][0],list_bb_compare[i][1],list_bb_compare[i][2],list_bb_compare[i][3]]
        det = False

        for j in range(len(result)) :
            bb_res = [result[j][1],result[j][2],result[j][3],result[j][4]]
            t_iou = iou_mine(bb_gt,bb_res)[0]

            print("\t \t IOU : ",t_iou)

            if t_iou > t_threshold :
                det = True;

        if det :
            num_cor += 1
        else :
            num_false+=1

    cv2.imshow('test',image)
    cv2.waitKey(0)

    #plt.imshow(image)
    #plt.savefig(os.getcwd() + '/' + name + 'output.png')
    return num_cor,num_false

if settings.output == 1:
    image = cv2.imread(settings.picture_name)
    read_image(image, settings.picture_name[-10:])

if settings.output == 2:
    labels = VOC('test').load_labels()
    #print(labels)

    totalData = 0
    t_tr = 0

    for i in range(len(labels)):
        #print(labels[i]['list_bb'])
        #print(labels[i]['imname'])
        image = cv2.imread(labels[i]['imname'])
        n_tr, n_fl = read_image(image, labels[i]['imname'][-10:],labels[i]['list_bb'])

        totalData+= (n_tr + n_fl)
        t_tr += n_tr

        print("True : ",n_tr,"False : ",n_fl,"/n***********************/n")

    print("total acc : ",(t_tr/totalData)*100,'%')

if settings.output == 3:
    cap = cv2.VideoCapture(settings.video_name)
    ret, _ = cap.read()
    while ret:
        ret, frame = cap.read()
        frame = cv2.resize(frame, (1280, 640))
        result = detect(frame)
        draw_result(frame, result)
        cv2.imshow('Camera', frame)
        cv2.waitKey(1)

max_stars_repo_name: deckyal/ImprovedYolo

max_issues_repo_name: deckyal/ImprovedYolo

max_forks_repo_name: deckyal/ImprovedYolo


text:
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon,Point
from .grids import GPS_to_grids,grids_centre
import math
import numpy as np

def sample_duration(data,col = ['Vehicleid','Time']):
    '''
    统计数据采样间隔

    输入
    -------
    data : DataFrame
        数据
    col : List
        列名，按[个体ID,时间]的顺序

    输出
    -------
    sample_duration : DataFrame
        一列的数据表，列名为duration，内容是数据的采样间隔，单位秒
    '''
    [Vehicleid,Time] = col
    data[Time] = pd.to_datetime(data[Time])
    data1 = data.copy()
    data1 = data1.sort_values(by = [Vehicleid,Time])
    data1[Vehicleid+'1'] = data1[Vehicleid].shift(-1)
    data1[Time+'1'] = data1[Time].shift(-1)
    data1['duration'] = (data1[Time+'1']-data1[Time]).dt.total_seconds()
    data1 = data1[data1[Vehicleid+'1']==data1[Vehicleid]]
    sample_duration = data1[['duration']]
    return sample_duration

def data_summary(data,col = ['Vehicleid','Time'],show_sample_duration = False):
    '''
    输入数据，打印数据概况

    输入
    -------
    data : DataFrame
        轨迹点数据
    col : List
        列名，按[个体ID，时间]的顺序
    sample_duration : bool
        是否输出个体采样间隔信息
    '''
    [Vehicleid,Time] = col
    print('数据量')
    print('-----------------')
    print('数据总量 :',len(data),'条')
    Vehicleid_count = data[Vehicleid].value_counts()
    print('个体总量 :',len(Vehicleid_count),'个')
    print('个体数据量均值 :',round(Vehicleid_count.mean(),2),'条')
    print('个体数据量上四分位 :',round(Vehicleid_count.quantile(0.75),2),'条')
    print('个体数据量中位数 :',round(Vehicleid_count.quantile(0.5),2),'条')
    print('个体数据量下四分位 :',round(Vehicleid_count.quantile(0.25),2),'条')
    print('')
    print('数据时间段')
    print('-----------------')
    print('开始时间 :',data[Time].min())
    print('结束时间 :',data[Time].max())
    print('')
    if show_sample_duration:
        sd = sample_duration(data, col=[Vehicleid, Time])
        print('个体采样间隔')
        print('-----------------')
        print('均值 :',round(sd['duration'].mean(),2),'秒')
        print('上四分位 :',round(sd['duration'].quantile(0.75),2),'秒')
        print('中位数 :',sd['duration'].quantile(0.5),'秒')
        print('下四分位 :',round(sd['duration'].quantile(0.25),2),'秒')

max_stars_repo_name: yan-yuchen/transbigdata

max_issues_repo_name: yan-yuchen/transbigdata

max_forks_repo_name: yan-yuchen/transbigdata


text:
#! /usr/bin/python3
# -*- coding: utf-8 -*-
"""
Attitude Filters Testing
========================

References
----------
.. [Crassidis] John L. Crassidis (2007) A Survey of Nonlinear Attitude
    Estimation Methods.
.. [Teage] Harris Teage (2016) Comparison of Attitude Estimation Techniques for
    Low-cost Unmanned Aerial Vehicles.
    https://arxiv.org/pdf/1602.07733.pdf
    http://ancs.eng.buffalo.edu/pdf/ancs_papers/2007/att_survey07.pdf
.. [Cirillo] A. Cirillo et al. (2016) A comparison of multisensor attitude
    estimation algorithms.
    https://www.researchgate.net/profile/Pasquale_Cirillo/publication/303738116_A_comparison_of_multisensor_attitude_estimation_algorithms/links/5750181208aeb753e7b4a0c0/A-comparison-of-multisensor-attitude-estimation-algorithms.pdf

"""

import numpy as np
import ahrs

RAD2DEG = ahrs.common.RAD2DEG
DEG2RAD = ahrs.common.DEG2RAD

class Test_Filter:
    def __init__(self, file_name, **kwargs):
        self.file = file_name
        self.data = ahrs.utils.io.load(self.file)

    def check_integrity(self, Q):
        sz = Q.shape
        qts_ok = not np.allclose(np.sum(Q, axis=0), sz[0]*np.array([1., 0., 0., 0.]))
        qnm_ok = np.allclose(np.linalg.norm(Q, axis=1).mean(), 1.0)
        return qts_ok and qnm_ok

    def allocate_arrays(self):
        self.Q = np.zeros((self.data.num_samples, 4))
        self.Q[:, 0] = 1.0

    def fourati(self):
        self.allocate_arrays()
        fourati = ahrs.filters.Fourati()
        for t in range(1, self.data.num_samples):
            self.Q[t] = fourati.update(DEG2RAD*self.data.gyr[t], self.data.acc[t], self.data.mag[t], self.Q[t-1])
        return self.check_integrity(self.Q)

    def ekf(self):
        self.allocate_arrays()
        ekf = ahrs.filters.EKF()
        for t in range(1, self.data.num_samples):
            self.Q[t] = ekf.update(DEG2RAD*self.data.gyr[t], self.data.acc[t], self.data.mag[t], self.Q[t-1])
        return self.check_integrity(self.Q)

    def mahony(self):
        self.allocate_arrays()
        mahony = ahrs.filters.Mahony()
        for t in range(1, self.data.num_samples):
            # self.Q[t] = mahony.updateIMU(DEG2RAD*self.data.gyr[t], self.data.acc[t], self.Q[t-1])
            self.Q[t] = mahony.updateMARG(DEG2RAD*self.data.gyr[t], self.data.acc[t], self.data.mag[t], self.Q[t-1])
        return self.check_integrity(self.Q)

    def madgwick(self):
        self.allocate_arrays()
        madgwick = ahrs.filters.Madgwick()
        for t in range(1, self.data.num_samples):
            # self.Q[t] = madgwick.updateIMU(DEG2RAD*self.data.gyr[t], self.data.acc[t], self.Q[t-1])
            self.Q[t] = madgwick.updateMARG(DEG2RAD*self.data.gyr[t], self.data.acc[t], self.data.mag[t], self.Q[t-1])
        return self.check_integrity(self.Q)


text:
from scipy.sparse import csr_matrix
import scipy.sparse as ssp
import networkx as nx
from .base import BaseGenerator
import numpy as np
from .. import register_feature

class Eigen(object):
    def __init__(self):
        pass

    def normalize_adj(self, adj):
        rowsum = np.array(adj.sum(1))
        d_inv_sqrt = np.power(rowsum, -0.5).flatten()
        d_inv_sqrt[np.isinf(d_inv_sqrt)] = 0.0
        d_inv_sqrt = ssp.diags(d_inv_sqrt)
        return adj.dot(d_inv_sqrt).transpose().dot(d_inv_sqrt)

    def forward(self, adj, d, use_eigenvalues=0, adj_norm=1):
        G = nx.from_scipy_sparse_matrix(adj)
        comp = list(nx.connected_components(G))
        results = np.zeros((adj.shape[0], d))
        for i in range(len(comp)):
            node_index = np.array(list(comp[i]))
            d_temp = min(len(node_index) - 2, d)
            if d_temp <= 0:
                continue
            temp_adj = adj[node_index, :][:, node_index].asfptype()
            if adj_norm == 1:
                temp_adj = self.normalize_adj(temp_adj)
            lamb, X = ssp.linalg.eigs(temp_adj, d_temp)
            lamb, X = lamb.real, X.real
            temp_order = np.argsort(lamb)
            lamb, X = lamb[temp_order], X[:, temp_order]
            for i in range(X.shape[1]):
                if np.sum(X[:, i]) < 0:
                    X[:, i] = -X[:, i]
            if use_eigenvalues == 1:
                X = X.dot(np.diag(np.sqrt(np.absolute(lamb))))
            elif use_eigenvalues == 2:
                X = X.dot(np.diag(lamb))
            results[node_index, :d_temp] = X
        return results

@register_feature("eigen")
class GeEigen(BaseGenerator):
    r"""concat Eigen features

    Notes
    -----
    An implementation of [1]_

    References
    ----------
    .. [1] Ziwei Zhang, Peng Cui, Jian Pei, Xin Wang, Wenwu Zhu:
        Eigen-GNN: A Graph Structure Preserving Plug-in for GNNs. CoRR abs/2006.04330 (2020)
        https://arxiv.org/abs/2006.04330

    Parameters
    ----------
    size : int
        EigenGNN hidden size
    """

    def __init__(self, size=32):
        super(GeEigen, self).__init__()
        self.size = size

    def _transform(self, data):

        adj = csr_matrix(
            (data.edge_weight, (data.edge_index[0, :], data.edge_index[1, :])),
            shape=(data.x.shape[0], data.x.shape[0]),
        )
        if np.max(adj - adj.T) > 1e-5:
            adj = adj + adj.T
        mf = Eigen()
        data.x = np.concatenate([data.x, mf.forward(adj, self.size)], axis=1)

        return data

    def _preprocess(self, data):
        if not hasattr(data, "edge_weight"):
            data.edge_weight = np.ones(data.edge_index.shape[1])

    def _postprocess(self, data):
        del data.edge_weight

max_issues_repo_name: lucylow/AutoGL

max_forks_repo_name: lucylow/AutoGL

reason: import numpy,import scipy,from scipy,import networkx


text:
#!/usr/bin/env python

"""APQBUNDLE.PY - This bundles the raw APOGEE frames and optionally compresses them.

The data will be output in this format:
 HDU0: header but NO data
 HDU1: header, read1 image  as UNSIGNED INTEGERS (BITPIX=16 or UINT)
 HDU2: header, read2 image  as UNSIGNED INTEGERS (BITPIX=16 or UINT)
 and so on for all the reads.

If /compress is set then APZIP will be used to compress the raw
APOGEE data.  Use APUNZIP to uncompress the data.
"""

from __future__ import print_function

__authors__ = 'David Nidever <dnidever@noao.edu>'
__version__ = '20180922'  # yyyymmdd

import os
import numpy as np
import warnings
from astropy.io import fits
from astropy.table import Table, Column
from glob import glob
from apogee.utils import yanny, apload
from sdss_access.path import path

# Ignore these warnings, it's a bug
warnings.filterwarnings("ignore", message="numpy.dtype size changed")
warnings.filterwarnings("ignore", message="numpy.ufunc size changed")

def apbundle(framenum,indir=None,outdir=None,compress=False,clobber=True,outfiles=None):
    """This bundles the raw APOGEE frames and optionally compresses them.

    The data will be output in this format:
     HDU0: header but NO data
     HDU1: header, read1 image  as UNSIGNED INTEGERS (BITPIX=16 or UINT)
     HDU2: header, read2 image  as UNSIGNED INTEGERS (BITPIX=16 or UINT)
     and so on for all the reads.

    If /compress is set then APZIP will be used to compress the raw
    APOGEE data.  Use APUNZIP to uncompress the data.

    Parameters
    ----------
    frameid : int or string
            The 8 character frame ID.
    indir : str, optional
          The input directory with the apRaw files.  If not given then
            the data model will be used to figure out the directory.
    outdir : str, optional
            The output directory for the apR-[abc]-frameid.fits files.
            If not given then the data model will be used.s
    compress : bool, optional
             Compress the output files.  False by default.
    clobber : bool, optional
            Overwrite files if they already exist. True by defalt.

    Returns
    -------
    outfiles : list
            The output files names.  Normally this is the bundles
            files, but if /compress is set then this is the
            compressed filenames.

    The raw bundles APOGEE files are names outdir/apR-[abc]-frameid.fits
    If /compress is set then compressed files will also be created and
    have names of outdir/apR-[abc]-frameid.apz.

    Example
    -------

    .. code-block:: python

        apbundle('00000117')

    """

    npix = 2048

    # ID8
    frameid = ("%08d" % np.int(framenum))
    mjd = apload.cmjd(np.int(framenum))

    # Get apRaw directory
    sdss_path = path.Path()
    rawdir = sdss_path.dir('apRaw',num=framenum,read=1,mjd=mjd)
    if observatory=='apo':
        prefix = 'ap'
    else:
        prefix = 'as'
    archive_dir = sdss_path.dir(prefix+'R',mjd=mjd,chip='a',num=1)
    if outdir is None: outdir = archive_dir

    # Find the files
    files = np.sort(glob(rawdir+'/'+prefix+'Raw-'+frameid+'-???.fits'))
    nfiles = len(files)
    if nfiles==0:
        raise Exception("No files for "+frameid)
    num = np.zeros(nfiles,dtype=int)
    for i in range(nfiles):
        base = os.path.basename(files[i])
        # apRaw-28190009-059.fits
        num[i] = np.int(base[15:18])
    maxind = np.argmax(num)
    print("Nfiles = "+str(nfiles)+" for "+frameid)

    # Read in the FIRST frame and initialize the output files
    #  this way the timestamp will be from the BEGINNING of
    #  the exposure
    minind = np.argmin(num)
    firstfile = files[minind]
    #FITS_READ,firstfile,im,head,exten=0,/no_abort,message=message
    im,head = fits.getdata(firstfile,0,header=True)

    chiptag = ['a','b','c']
    outfiles = []
    for ch in chiptag:
        o = os.path.join(outdir,os.path.basename(sdss_path.full(prefix+'R',mjd=mjd,chip=ch,num=frameid).replace('apz','fits')))
        outfiles.append(o)

    # Final Xsize
    shape = im.shape
    xsizemap = {6144: 2048, 8192: 2048+512}
    if shape[1] not in xsizemap:
        raise Exception('Non-standard Xsize = '+str(sz[1])+' not supported')
    xsize = xsizemap[shape[1]]

    # Check if the output files exist
    test = [os.path.exists(o) for o in outfiles]
    if (np.sum(test)>0) & (clobber is False):
      raise Exception('Output files exist already.  Set /clobber to overwrite')

    print('Writing bundled file to '+" ".join(outfiles))

    # Initialize the three bundled files
    ahead = head.copy()
    ahead["SIMPLE"] = "T"
    ahead["BITPIX"] = 16
    ahead["NAXIS"] = 0
    ahead["CHIP"] = "a"
    ahead["NREAD"] = nfiles
    remkeys = ["NAXIS1","NAXIS2","PCOUNT","GCOUNT","CHECKSUM","DATASUM","BZERO","BSCALE"]
    for k in remkeys:
        if ahead.get(k) is not None: ahead.remove(k)
    fits.writeto(outfiles[0],None,ahead,checksum=True,overwrite=True)

    bhead = ahead.copy()
    bhead["CHIP"] = "b"
    remkeys = ["CHEKSUM","DATASUM"]
    for k in remkeys:
        if bhead.get(k) is not None: bhead.remove(k)
    fits.writeto(outfiles[1],None,bhead,checksum=True,overwrite=True)

    chead = ahead
    chead["CHIP"] = "c"
    remkeys = ["CHEKSUM","DATASUM"]
    for k in remkeys:
        if chead.get(k) is not None: chead.remove(k)
    fits.writeto(outfiles[2],None,chead,checksum=True,overwrite=True)

    # Read in the frames
    for fpos in range(nfiles):
        # to use the fits_test_checksum, the data has to be read as a signed int first
        # this is because the pyfits checksum is done on signed values
        im,head = fits.getdata(files[fpos],0,header=True)
        #FITS_READ,files[fpos],im,head,exten=0,/no_abort,/no_unsigned,/noscale,/no_pdu,message=message
        #res = FITS_TEST_CHECKSUM(head, im, errmsg=errmsg)
        #if res eq -1 then begin
        #  # the checksum doesn't match -> send a warning and keep going
        #  message,/con,'BAD checksum for file '+files[fpos]
        #  return

        # Convert back to uint
        bzero = head.get("ZERO")
        if bzero is None: bzero=0
        if bzero!=0:
            im = np.uint(im)-np.uint(bzero[0])
        bscale = head.get("BSCALE")
        if bscale is None: bscale=1
        if bscale!=1:
            im /= bscale[0]

        # The image is: red, green, blue
        # if the 5th reference output is also there then
        # this is: green-window (not used), blue-ref, green-ref, red-ref

        # The RAW images:
        # Basic orientation is long wavelengths on the left, to shorter to the
        # right (higher column number).  readout direction is:
        #         Red                       Green                 Blue           WG RefB RefG RefR
        #|---->|<----|--->|<----||---->|<----|--->|<----||---->|<----|--->|<----||---->|<----|--->|<----|

        # Chip a, RED
        #--------------
        im1 = im[:,0:npix]
        # Add the 5th output, 4th column
        #  need to flip in X-direction so all three chips
        #  have the reference reading out left-to-right (--->)
        if xsize==2560:
            refim1 = im[:,6144+3*512:6144+4*512]
            refim1 = np.flip(refim1,axis=1)  # flip X-direction
            im1 = np.hstack((im1,refim1))
        ahead = head.copy()
        ahead.set("XTENSION","IMAGE",before="SIMPLE")
        bitpix = ahead.get("BITPIX")
        if bitpix is not None: ahead["BITPIX"] = bitpix
        ahead["NAXIS"] = 2
        ahead.set("NAXIS1",xsize,after="NAXIS")
        ahead.set("NAXIS2",len(im1[0,:]),after="NAXIS1")
        ahead.set("PCOUNT",0,after="NAXIS2")
        ahead.set("GCOUNT",1,after="PCOUNT")
        # delete SIMPLE if present, only allowed in PDU
        remkeys = ["SIMPLE","CHECKSUM","DATASUM"]
        for k in remkeys:
            if ahead.get(k) is not None: ahead.remove(k)
        ahead["CHIP"] = "a"
        ahead["NREAD"] = nfiles
        fits.append(outfiles[0],im1,ahead,checksum=True)

        # re-write the bzero and bscale to add the '/' that mwrfits automatically adds
        # so as to not mess up the checksum
        #if bzcount>0: sxaddpar,ahead,'BZERO',bzero
        #if bscount>0: sxaddpar,ahead,'BSCALE', bscale

        # Chip b, GREEN
        #---------------
        im2 = im[:,npix:2*npix-1]
        # Add the 5th output, 3rd column
        #  the reference readout direction is already left-to-right
        if xsize==2560:
            refim2 = im[:,6144+2*512:6144+3*512]
            im2 = np.hstack((im2,refim2))
        bhead = ahead.copy()
        bhead["NAXIS"] = 2
        bhead.set("NAXIS1",xsize,after="NAXIS")
        bhead.set("NAXIS2",len(im2[0,:]),after="NAXIS1")
        bhead.set("PCOUNT",0,after="NAXIS2")
        bhead.set("GCOUNT",1,after="PCOUNT")
        fits.append(outfiles[1],im2,bhead,checksum=True)

        # Chip c, BLUE
        #--------------
        im3 = im[:,2*npix:3*npix-1]
        # Add the 5th output, 2nd column
        #  need to flip in X-direction so all three chips
        #  have the reference reading out left-to-right (--->)
        if xsize==2560:
            refim3 = im[:,6144+512:6144+2*512]
            refim3 = np.flip(refim3,axis=1)  # flip X-direction
            im3 = np.hstack((im3,refim3))
        chead = ahead.copy()
        chead["NAXIS"] = 2
        chead.set("NAXIS1",xsize,after="NAXIS")
        chead.set("NAXIS2",len(im3[0,:]),after="NAXIS1")
        chead.set("PCOUNT",0,after="NAXIS2")
        chead.set("GCOUNT",1,after="PCOUNT")
        fits.append(outfiles[2],im3,chead,checksum=True)

        # Compress
        #----------
        if compress is True:
            print,'Compressing ',addslash(outdir)+file_basename(apogee_filename('R',num=frameid,chip=chiptag),'.apz')+'.fits'
            apzip(outfiles)

            # Return compressed filenames
            bundle_files = outfiles
            outfiles = file_dirname(bundle_files)+'/'+file_basename(bundle_files,'.fits')+'.apz'

    dt = systime(1)-t0
    print('dt = '+str(dt)+' sec')

    # 17 sec for a 60 read file.
    # 204 sec = 3.4 min for 60 read file WITH compression

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee

reason: import numpy,from astropy