text:
"""
Module that contains fonction for adcp data quality control.

Based on the adcp2nc package by jeanlucshaw: https://github.com/jeanlucshaw/adcp2nc/

Notes
-----
   Tests return `True` where cells fail a test.

   - Cells that haven't fail any test are given a flag value of `1` (probably_good_value).
   - Velocities in any direction are set to NaN if greater than 15 meters per second gets a 4.
   - Failing amplitude, correlation, percentgood, roll, pitch, side_lobe, horizontal
   velocity or vertical velocity test returns a flag_value of `3` (probably_bad_value)
   for the corresponding velocity cells.
   - Temperatures outside [-2, 32] Celsius have a flag_value of `4` (bad_value).
   - Pressures outside [0, 180] dbar value have a flag_value of `4` (bad_value).
   and the corresponding velocity cells have a flag_value of `3` (probably_bad_value)
   - The amplitude, correlation and percentgood thresholds are also applied to
   sentinelV fifth beam data.

   SeaDataNet Quality Control Flags Value
   * 0: no_quality_control
   * 1: good_value
   * 2: probably_good_value
   * 3: probably_bad_value
       - Unusual data value, inconsistent with real phenomena.
   * 4: bad_value
       - Obviously erroneous data value.
   * 5: changed_value
   * 6: value_below_detection
   * 7: value_in_excess
   * 8: interpolated_value
   * 9: missing_value

NOTE
IML flags meaning : (Basically the same)
   * 0: no_quality_control
   * 1: value_seems_correct
   * 2: value_appears_inconsistent_with_other_values
   * 3: values_seems_doubtful
   * 4: value_seems_erroneous
   * 5: value_was_modified
   * 9: value_missing

"""
import typing as tp

import numpy as np
import xarray as xr
from magtogoek.tools import circular_distance
from magtogoek.utils import Logger
from pandas import Timestamp
from scipy.stats import circmean

# Brand dependent quality control defaults
#    rti_qc_defaults = dict(amp_th=20)
#    rdi_qc_defaults = dict(amp_th=0)

IMPLAUSIBLE_VEL_TRESHOLD = 15  # meter per second
MIN_TEMPERATURE = -2  # Celcius
MAX_TEMPERATURE = 32  # Celcius
MIN_PRESSURE = 0  # dbar
MAX_PRESSURE = 10000  # dbar (mariana trench pressure)

l = Logger(level=0)
FLAG_REFERENCE = "BODC SeaDataNet"
FLAG_VALUES = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
FLAG_MEANINGS = (
    "no_quality_control",
    "good_value",
    "probably_good_value",
    "probably_bad_value",
    "bad_value",
    "changed_value",
    "value_below_detection",
    "value_in_excess",
    "interpolated_value",
    "missing_value",
)

def no_adcp_quality_control(dataset: xr.Dataset):
    """Adds var_QC ancillary variables to dataset with value 0.

    AAncillary variables:  temperature, pres, u, v, w.

    SeaDataNet Quality Control Flags Value
    * 0: no_quality_control

    Parameters
    ----------
    dataset :
        ADCP dataset formatted as done by adcp_init.
    """
    l.reset()
    l.section("No Quality Controlled", t=True)

    l.log("No quality control carried out")
    variables = ["temperature", "pres", "u", "v", "w"]
    for var in variables:
        if var in dataset:
            dataset[var + "_QC"] = dataset[var].copy().astype("int8") * 0

    dataset.attrs["flags_reference"] = FLAG_REFERENCE
    dataset.attrs["flags_values"] = FLAG_VALUES
    dataset.attrs["flags_meanings"] = FLAG_MEANINGS
    dataset.attrs["quality_comments"] = "No quality control."

def adcp_quality_control(
    dataset: xr.Dataset,
    amp_th: int = 30,
    corr_th: int = 64,
    pg_th: int = 90,
    roll_th: float = 20,
    pitch_th: float = 20,
    horizontal_vel_th: float = 5,
    vertical_vel_th: float = 5,
    error_vel_th: float = 5,
    motion_correction_mode: str = "",
    sidelobes_correction: bool = False,
    bottom_depth: float = None,
):
    """
    Perform ADCP quality control.

    This was adapted from jeanlucshaw adcp2nc package.

    Parameters
    ----------
    dataset :
        ADCP dataset formatted as done by adcp_init.
    amp_th :
        Require more than this amplitude values.
    pg_th :
        Require more than this percentage of good 4-beam transformations.
    corr_th :
        Require more than this beam correlation value.
    roll_th :
        Require roll values be smaller than this value (degrees).
    pitch_th :
        Require pitch values be smaller than this value (degrees).
    horizontal_vel_th:
        Require u, v  values be smaller than this value (meter per seconds).
    vertical_vel_th :
        Require w values be smaller than this value (meter per seconds).
    error_vel_th :
        Require e values be smaller than this value (meter per seconds).
    motion_correction_mode :
        If 'nav' or 'bt' will corrected velocities from the platform motion,
        will either correct u,v,w with navigation ('nav') or bottom track ('bt')
        data. No motion correction si carried out if motion_correction == 'off'.
    sidelobes_correction :
        Use fixed depth or bottom track range to remove side lobe
        contamination. Set to either "dep" or "bt" or None.
    bottom_depth :
        If not `None`, this depth used for removing side lobe contamination.
    Notes
    -----
       Tests return `True` where cells fail a test.

       - Cells that haven't fail any test are given a flag value of `2` (probably_good_value).
       - Failing amplitude, correlation, percentgood, roll, pitch, side_lobe, horizontal
       velocity or vertical velocity test returns a flag_value of `3` (probably_bad_value)
       for the corresponding velocity cells.
       - Temperatures outside [-2, 32] Celsius have a flag_value of `4` (bad_value).
       - Pressures outside [0, 180] dbar value have a flag_value of `4` (bad_value).
       and the corresponding velocity cells have a flag_value of `3` (probably_bad_value)
       - The amplitude, correlation and percentgood thresholds are also applied to
       sentinelV fifth beam data.

       SeaDataNet Quality Control Flags Value
       * 0: no_quality_control
       * 1: good_value
       * 2: probably_good_value
       * 3: probably_bad_value
           - Unusual data value, inconsistent with real phenomena.
       * 4: bad_value
           - Obviously erroneous data value.
       * 5: changed_value
       * 6: value_below_detection
       * 7: value_in_excess
       * 8: interpolated_value
       * 9: missing_value
    """
    l.reset()
    l.section("Quality Control")

    if motion_correction_mode in ["bt", "nav"]:
        motion_correction(dataset, motion_correction_mode)

    # vel_flags = set_implausible_vel_to_nan(dataset, thres=IMPLAUSIBLE_VEL_TRESHOLD)

    vel_flags = np.ones(dataset.depth.shape + dataset.time.shape).astype(int)

    vel_qc_test = []

    if amp_th:
        l.log(f"amplitude threshold {amp_th}")
        amp_flag = amplitude_test(dataset, amp_th)
        vel_flags[amp_flag] = 3
        vel_qc_test.append(f"amplitude_threshold:{amp_th}")

    if corr_th:
        l.log(f"correlation threshold {corr_th}")
        corr_flag = correlation_test(dataset, corr_th)
        vel_flags[corr_flag] = 3
        vel_qc_test.append(f"correlation_threshold:{corr_th}")

    if pg_th:
        l.log(f"percentgood threshold {pg_th}")
        pg_flag = percentgood_test(dataset, pg_th)
        vel_flags[pg_flag] = 3
        vel_qc_test.append(f"percentgood_threshold:{pg_th}")

    if horizontal_vel_th:
        l.log(f"horizontal velocity threshold {horizontal_vel_th} m/s")
        horizontal_vel_flag = horizontal_vel_test(dataset, horizontal_vel_th)
        vel_flags[horizontal_vel_flag] = 3
        vel_qc_test.append(f"horizontal_velocity_threshold:{horizontal_vel_th} m/s")

    if vertical_vel_th:
        l.log(f"vertical velocity threshold {vertical_vel_th} m/s")
        vertical_vel_flag = vertical_vel_test(dataset, vertical_vel_th)
        vel_flags[vertical_vel_flag] = 3
        vel_qc_test.append(f"vertical_velocity_threshold:{vertical_vel_th} m/s")

    if error_vel_th:
        l.log(f"error velocity threshold {error_vel_th} m/s")
        error_vel_flag = error_vel_test(dataset, error_vel_th)
        vel_flags[error_vel_flag] = 3
        vel_qc_test.append(f"velocity_error_threshold:{error_vel_th} m/s")

    if roll_th:
        l.log(f"roll threshold {roll_th} degree")
        roll_flag = np.tile(roll_test(dataset, roll_th), dataset.depth.shape + (1,))
        vel_flags[roll_flag] = 3
        vel_qc_test.append(f"roll_threshold:{roll_th} degree")

    if pitch_th:
        l.log(f"pitch threshold {pitch_th} degree")
        pitch_flag = np.tile(pitch_test(dataset, pitch_th), dataset.depth.shape + (1,))
        vel_flags[pitch_flag] = 3
        vel_qc_test.append(f"pitch_threshold:{pitch_th} degree")

    if sidelobes_correction:
        sidelobe_flag = sidelobe_test(dataset, bottom_depth)
        if isinstance(sidelobe_flag, np.ndarray):
            l.log("Sidelobe correction carried out.")
            vel_flags[sidelobe_flag] = 3
            vel_qc_test.append("sidelobes")

    if "pres" in dataset:
        l.log(f"Good pressure range {MIN_PRESSURE} to {MAX_PRESSURE} dbar")
        pressure_QC = np.ones(dataset.pres.shape)
        pressure_flags = pressure_test(dataset)
        pressure_QC[pressure_flags] = 4
        dataset["pres_QC"] = (["time"], pressure_QC)
        vel_flags[np.tile(pressure_flags, dataset.depth.shape + (1,))] = 3
        dataset["pres_QC"].attrs[
            "quality_test"
        ] = f"pressure_threshold: less than {MIN_PRESSURE} dbar and greater than {MAX_PRESSURE} dbar"
        vel_qc_test.append(dataset["pres_QC"].attrs["quality_test"])

    if "temperature" in dataset:
        l.log(f"Good temperature range {MIN_TEMPERATURE} to {MAX_TEMPERATURE} celsius")
        temperature_QC = np.ones(dataset.temperature.shape)
        temperature_QC[temperature_test(dataset)] = 4
        dataset["temperature_QC"] = (["time"], temperature_QC)
        dataset["temperature_QC"].attrs[
            "quality_test"
        ] = f"temperature_threshold: less than {MIN_TEMPERATURE} Celsius and greater than {MAX_TEMPERATURE} celsius"
    if "vb_vel" in dataset:
        l.log(
            "Fifth beam quality control carried out with"
            + f"amplitude_threshold: {amp_th}" * ("vb_amp" in dataset)
            + f", correlation_threshold: {corr_th}" * ("vb_corr" in dataset)
            + f", percentgood_threshold: {pg_th}" * ("vb_pg" in dataset)
            + "."
        )
        vb_flag = vertical_beam_test(dataset, amp_th, corr_th, pg_th)
        dataset["vb_vel_QC"] = (["depth", "time"], vb_flag * 3)
        dataset["vb_vel_QC"].attrs["quality_test"] = (
            f"amplitude_threshold: {amp_th}\n" * ("vb_amp" in dataset)
            + f"correlation_threshold: {corr_th}\n" * ("vb_corr" in dataset)
            + f"percentgood_threshold: {pg_th}\n" * ("vb_pg" in dataset)
        )

    vel_flags[flag_implausible_vel(dataset, threshold=IMPLAUSIBLE_VEL_TRESHOLD)] = 4

    missing_vel = np.bitwise_or(
        *(~np.isfinite(dataset[v].values) for v in ("u", "v", "w"))
    )
    vel_flags[missing_vel] = 9

    for v in ("u", "v", "w"):
        dataset[v + "_QC"] = (["depth", "time"], vel_flags)
        dataset[v + "_QC"].attrs["quality_test"] = "\n".join(vel_qc_test)

    for var in list(dataset.variables):
        if "_QC" in var:
            dataset[var].attrs["quality_date"] = Timestamp.now().strftime("%Y-%m-%d")
            dataset[var].attrs["flag_meanings"] = FLAG_MEANINGS
            dataset[var].attrs["flag_values"] = FLAG_VALUES
            dataset[var].attrs["flag_reference"] = FLAG_REFERENCE

    dataset.attrs["quality_comments"] = l.logbook[len('[Quality Control]'):]
    l.log(f"Quality Control was carried out with {l.w_count} warnings")
    dataset.attrs["logbook"] += l.logbook

    dataset.attrs["flags_reference"] = FLAG_REFERENCE
    dataset.attrs["flags_values"] = FLAG_VALUES
    dataset.attrs["flags_meanings"] = FLAG_MEANINGS

def motion_correction(dataset: xr.Dataset, mode: str):
    """Carry motion correction on velocities.

    If mode is 'bt' the motion correction is along x, y, z.
    If mode is 'nav' the motion correction is along x, y.
    """
    if mode == "bt":
        if all(f"bt_{v}" in dataset for v in ["u", "v", "w"]):
            for field in ["u", "v", "w"]:
                dataset[field].values -= dataset[f"bt_{field}"].values
            l.log("Motion correction carried out with bottom track")
        else:
            l.warning(
                "Motion correction aborted. Bottom velocity (bt_u, bt_v, bt_w) missing"
            )
    elif mode == "nav":
        if all(f"{v}_ship" in dataset for v in ["u", "v"]):
            for field in ["u", "v"]:
                dataset[field] += np.tile(
                    dataset[field + "_ship"].where(np.isfinite(dataset.lon.values), 0),
                    (dataset.depth.size, 1),
                )
            l.log("Motion correction carried out with navigation")
        else:
            l.warning(
                "Motion correction aborted. Navigation velocity (u_ship, v_ship) missing"
            )
    else:
        l.warning(
            "Motion correction aborted. Motion correction mode invalid. ('bt' or 'nav')"
        )

def flag_implausible_vel(
    dataset: xr.Dataset, threshold: float = 15
) -> tp.Type[np.array]:
    """Values greater than `threshold` return True"""
    return (
        (dataset.v > threshold) & (dataset.u > threshold) & (dataset.w > threshold)
    ).data

def set_implausible_vel_to_nan(dataset: xr.Dataset, threshold: float = 15):
    """Set bin with improbable values to Nan."""
    for v in ["u", "v", "w"]:
        plausible = (dataset[v] > -threshold) & (dataset[v] < threshold)
        dataset["u"] = dataset["u"].where(plausible)
        dataset["v"] = dataset["v"].where(plausible)
        dataset["w"] = dataset["w"].where(plausible)

def correlation_test(dataset: xr.Dataset, threshold: int):
    """FIXME
    Value must be greater than the threshold to be good. (True fails)
    NOTE JeanLucShaw used absolute but is it needed ?"""

    if all(f"corr{i}" in dataset for i in range(1, 5)):
        return (
            (dataset.corr1 < threshold)
            & (dataset.corr2 < threshold)
            & (dataset.corr3 < threshold)
            & (dataset.corr4 < threshold)
        ).data
    else:
        l.warning("Correlation test aborted. Missing one or more corr data")
        return np.full(dataset.depth.shape + dataset.time.shape, False)

def amplitude_test(dataset: xr.Dataset, threshold: int):
    """FIXME
    Value must be greater than the threshold to be good. (True fails)
    NOTE JeanLucShaw used absolute but is it needed ?"""
    if all(f"amp{i}" in dataset for i in range(1, 5)):
        return (
            (dataset.amp1 < threshold)
            & (dataset.amp2 < threshold)
            & (dataset.amp3 < threshold)
            & (dataset.amp4 < threshold)
        ).data
    else:
        l.warning("Amplitude test aborted. Missing one or more corr data")
        return np.full(dataset.depth.shape + dataset.time.shape, False)

def percentgood_test(dataset: xr.Dataset, threshold: int):
    """FIXME
    Value must be greater than the threshold to be good. (True fails)
    NOTE JeanLucShaw used absolute but is it needed ?"""
    if "pg" in dataset:
        return (dataset.pg < threshold).data
    else:
        l.warning("Percent Good test aborted. Missing one or more corr data")
        return np.full(dataset.depth.shape + dataset.time.shape, False)

def roll_test(dataset: xr.Dataset, threshold: float) -> tp.Type[np.array]:
    """FIXME
    Roll conditions (True fails)
    Distance from mean"""
    if "roll_" in dataset:
        roll_mean = circmean(dataset.roll_.values, low=-180, high=180)
        roll_from_mean = circular_distance(dataset.roll_.values, roll_mean, units="deg")
        return roll_from_mean > threshold
    else:
        l.warning("Roll test aborted. Missing one or more corr data")
        return np.full(dataset.depth.shape + dataset.time.shape, False)

def pitch_test(dataset: xr.Dataset, threshold: float) -> tp.Type[np.array]:
    """FIXME
    Pitch conditions (True fails)
    Distance from Mean
    """
    if "pitch" in dataset:
        pitch_mean = circmean(dataset.pitch.values, low=-180, high=180)
        pitch_from_mean = circular_distance(
            dataset.pitch.values, pitch_mean, units="deg"
        )
        return pitch_from_mean > threshold

    else:
        l.warning("Pitch test aborted. Missing one or more corr data")
        return np.full(dataset.depth.shape + dataset.time.shape, False)

def horizontal_vel_test(
    dataset: xr.Dataset, threshold: float
) -> tp.Type[np.array]:
    """FIXME
    None finite value value will also fail"""

    horizontal_velocity = np.sqrt(dataset.u ** 2 + dataset.v ** 2)

    return np.greater(
        horizontal_velocity.values, threshold, where=np.isfinite(horizontal_velocity),
    )

def vertical_vel_test(dataset: xr.Dataset, threshold: float) -> tp.Type[np.array]:
    """FIXME
    None finite value value will also fail"""
    return np.greater(
        abs(dataset.w.values), threshold, where=np.isfinite(dataset.w.values),
    )

def error_vel_test(dataset: xr.Dataset, threshold: float) -> tp.Type[np.array]:
    """FIXME
    None finite value value will also fail"""
    return np.greater(
        abs(dataset.e.values), threshold, where=np.isfinite(dataset.w.values),
    )

def vertical_beam_test(
    dataset: xr.Dataset, amp_threshold: float, corr_threshold: float, pg_threshold: float
) -> tp.Type[np.array]:
    """FIXME"""
    vb_test = np.full(dataset.depth.shape + dataset.time.shape, False)
    if "vb_amp" in dataset.variables and amp_threshold:
        vb_test[dataset.vb_amp < amp_threshold] = True
    if "vb_corr" in dataset.variables and corr_threshold:
        vb_test[dataset.vb_corr < corr_threshold] = True
    if "vb_pg" in dataset.variables and pg_threshold:
        vb_test[dataset.vb_pg < pg_threshold] = True

    return vb_test

def sidelobe_test(dataset: xr.Dataset, bottom_depth: float = None):
    """FIXME
    Test for sidelobe contamination (True fails).

    Returns a boolean array or a False statement if the test cannot be carried out.

    Downward equation:
        bin_depth + 0.5*bin_size > XducerDepth + (bottom_depth - XducerDepth)*cos(beam_angle)
    Upward equation:
        bin_depth + 0.5*bin_size < XducerDepth*( 1 - cos(beam_angle))

    Parameters
    ----------
    dataset :
    bottom_depth : optional
        Fixed bottom depth to use for sidelobe correction
    """
    if dataset.attrs["beam_angle"] and dataset.attrs["orientation"]:
        angle_cos = np.cos(np.radians(dataset.attrs["beam_angle"]))
        depth_array = (
            np.tile(dataset.depth.data, dataset.time.shape + (1,))
            + dataset.attrs["bin_size_m"] / 2
        )

        if "xducer_depth" in dataset.attrs:
            xducer_depth = np.tile(dataset.attrs["xducer_depth"], dataset.time.shape)
        elif "xducer_depth" in dataset:
            xducer_depth = dataset["xducer_depth"].data
        else:
            l.warning(
                "Sidelobes correction aborted. Adcp depth `xducer_depth` not provided."
            )
            return False

        if dataset.attrs["orientation"] == "down":
            if "bt_depth" in dataset:
                bottom_depth = dataset.bt_depth.data
            elif not bottom_depth:
                l.warning(
                    "Sidelobes correction aborted. Bottom depth not found or provided."
                )
                return False

            max_depth = xducer_depth + (bottom_depth - xducer_depth) * angle_cos
            return depth_array.T > max_depth

        elif dataset.attrs["orientation"] == "up":
            min_depth = xducer_depth * (1 - angle_cos)
            return depth_array.T < min_depth.T

        else:
            l.warning(
                "Can not correct for sidelobes, `adcp_orientation` parameter not set. Must be `up` or `down`."
            )
            return False
    else:
        l.warning("Can not correct for sidelobes, beam_angle not found in raw file.")
        return False

def temperature_test(dataset: xr.Dataset):
    """FIXME"""
    return np.bitwise_or(
        dataset.temperature > MAX_TEMPERATURE, dataset.temperature < MIN_TEMPERATURE
    ).data

def pressure_test(dataset: xr.Dataset):
    """FIXME"""
    return np.bitwise_or(dataset.pres > MAX_PRESSURE, dataset.pres < MIN_PRESSURE).data

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    from magtogoek.adcp.loader import load_adcp_binary

    sillex_path = "/media/jeromejguay/5df6ae8c-2af4-4e5b-a1e0-a560a316bde3/home/jeromejguay/WorkSpace_2019/Data/Raw/ADCP/"
    sillex_fns = ["COR1805-ADCP-150kHz009_000001", "COR1805-ADCP-150kHz009_000002"]

    v50_files = (
        "/media/jeromejguay/Bruno/TREX2020/V50/TREX2020_V50_20200911T121242_003_*.ENX"
    )

    pd0_sw_path = "/home/jeromejguay/ImlSpace/Projects/magtogoek/test/files/sw_300_4beam_20deg_piston.pd0"

    ens_sw_path = (
        "/home/jeromejguay/ImlSpace/Projects/magtogoek/test/files/rowetech_seawatch.ens"
    )

    test = "SW_PD0"

    if test == "SV":
        _dataset = load_adcp_binary(v50_files, sonar="sv", yearbase=2020, orientation="down",)

    if test == "ENX":
        _dataset = load_adcp_binary(
            #            [sillex_path + fn + ".ENX" for fn in sillex_fns],
            sillex_path + "COR1805-ADCP-150kHz009.ENX",
            sonar="os",
            yearbase=2018,
            orientation="down",
        )
    if test == "SW_PD0":
        _dataset = load_adcp_binary(
            pd0_sw_path, sonar="sw_pd0", yearbase=2020, orientation="down"
        )
    if test == "ENS":
        _dataset = load_adcp_binary(
            ens_sw_path,
            sonar="sw",
            yearbase=2020,
            orientation="down",
            leading_index=None,
            trailing_index=None,
        )

    adcp_quality_control(_dataset, roll_th=20, pitch_th=20, horizontal_vel_th=2, vertical_vel_th=0.1,
                         motion_correction_mode="bt", sidelobes_correction=True, bottom_depth=None)

    _dataset.u.where(_dataset.u_QC == 1).plot()

    plt.show()


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.
import numpy as np

from nn_meter.utils.graph_tool import ModelGraph
from .frozenpb_parser import FrozenPbParser
from .shape_inference import ShapeInference
from .shape_fetcher import ShapeFetcher

class FrozenPbConverter:
    def __init__(self, file_name):
        self.model_graph = ModelGraph()

        # Parse pb to graph
        parser = FrozenPbParser(file_name)
        parser.parse_graph(self.model_graph)
        dynamic_fetcher = ShapeFetcher(parser.graph)

        # Change split to more firendly scheme
        parser.fix_split_naming(self.model_graph)

        # Get the static shape
        ShapeInference(self.model_graph, dynamic_fetcher)

        # Strip constant and indentity nodes
        parser.strip_useless_nodes(self.model_graph)

    def get_flatten_graph(self):
        def np_encoder(d):
            for k, v in d.items():
                if isinstance(v, dict):
                    np_encoder(v)
                else:
                    if isinstance(v, np.ndarray):
                        d[k] = v.tolist()
                    if isinstance(v, (bytes, bytearray)):
                        d[k] = v.decode("utf-8")

        np_encoder(self.model_graph.get_graph())
        return self.model_graph.get_graph()

ext: py

max_stars_repo_name: kaleid-liner/nn-Meter

max_issues_repo_name: kaleid-liner/nn-Meter

max_forks_repo_name: kaleid-liner/nn-Meter

reason: import numpy


text:
import numpy as np
import itertools
import functools

# region Octahedral rotation group
"""
Implementation of Octahedral Rotation group with matrix representation from [2]

[2]:
@article{Lenz2009,
author = {Lenz, Reiner},
file = {:home/anonymous/data/library/papers/Lenz - 2009 - Octahedral Filters for 3D Image Processing.pdf:pdf},
journal = {Proceedings SSBA 2009},
mendeley-groups = {group theory},
pages = {109--112},
title = {{Octahedral Filters for 3D Image Processing}},
year = {2009}
}
"""

class OctahedralMatrixRotationGroup2Generators:
    """
    Implementation of Octahedral rotation group. Isomorphic to S4 freely generated with two elements:

        G = < A, D | DDD = AAAA = ADDA(1/D)= ADA(1/D)(1/A)(1/A)(1/D) = e

    """
    # The two generator elements for matrix
    _R1 = np.array(((0, 0, 1), (1, 0, 0), (0, 1, 0)))
    _R2 = np.array(((0, 1, 0), (-1, 0, 0), (0, 0, 1)))

    # Identifiers of generators in string word
    _R1_id = 'D'
    _R2_id = 'A'

    # Complete representation as shortest words
    _elements_as_words = ['', 'A', 'D', 'AA', 'AD', 'DA', 'DD', 'AAA', 'AAD', 'ADA', 'ADD', 'DAA', 'DAD', 'DDA',
                          'AADA', 'AADD', 'ADAA', 'ADAD', 'DADA', 'DADD', 'DDAA', 'DDAD', 'DADAA', 'DADAD']

    @property
    def elements(self):
        return self._elements_as_words

    @property
    def generators(self):
        return [self._R1_id, self._R2_id]

    @classmethod
    def word_to_matrix(cls, word):
        matrix_word = []
        for c in word:
            if c == 'D':
                matrix_word.append(cls._R1)
            if c == 'A':
                matrix_word.append(cls._R2)

        return functools.reduce(np.matmul, matrix_word)

    def __iter__(self):
        return iter(self._elements_as_words)

# endregion

# region Lebedev orbits
"""
Class names were given with respect to [1].
In this a Lebedev point is a tuple of length 2 where
    point[0] -> string identifier of the class to which the lebedev point contains.
    point[1] -> its identifier within the class.

[1]: https://en.wikipedia.org/wiki/Lebedev_quadrature#Construction

This region implements:
    1. Lebedev point sets a1, a2, a3

    2. A Lebedev Grid with 26 points.

"""

class _LebedevOrbitBase:
    """
    Base class of a LebedevOrbit. A Lebedev Orbit is a point set on S^2 which is invariant
    under the Octahedral Symetry Group.
    """
    # Possible point numbers for points of child class. Should be [str]
    _point_numbers = None

    # id connecting a LebedevPoint with its class. Should be str
    _string_id = None
    _points = None

    def __init__(self):
        self._spherical_dict = self._initial_state_spherical_dict()

    @staticmethod
    def _initial_state_spherical_dict():
        raise NotImplementedError("Abstract method.")

    @classmethod
    def _get_points(cls):
        for num in cls._point_numbers:
            yield (cls._string_id, num)

    def _check_point(self, point):
        point = tuple(point)

        if len(point) != 2:
            raise ValueError("{} is no point of a Lebedev point set".format(point))

        if point[0] != self._string_id:
            raise ValueError("{} does not contain to {}.".format(point, type(self).__name__))

        if point[1] not in self._point_numbers:
            raise ValueError("{} has no point with number {}.".format(type(self).__name__, point.number))

        return point

    @staticmethod
    def __snap_to_zero_one(value):
        if np.isclose([value], [1]):
            return 1
        elif np.isclose([value], [0]):
            return 0
        else:
            return value

    def to_cartesian(self, point: tuple)->tuple:
        """
        Returns cartesian coordinates of point.

        Parameters
        ----------
        point : Lebedev Point.

        Returns
        -------
            tuple. (x,y,z) cartesian coordinates of point.
        """
        point = self._check_point(point)

        polar = self.to_spherical(point)
        theta = polar[0] * np.pi / 180
        phi = polar[1] * np.pi / 180

        x = np.cos(theta) * np.sin(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(phi)

        return tuple(map(self.__snap_to_zero_one, [x, y, z]))

    def to_spherical(self, point: tuple):
        """
        Returns spherical coordinates of point.

        Parameters
        ----------
        point :
            tuple. Lebedev Point.

        Returns
        -------
            tuple. (theta, phi) spherical coordinates in degree of point:

                (theta, phi) \in [-180, 180] x [0, 180]
        """
        point = self._check_point(point)
        point_number = point[1]
        return self._spherical_dict[point_number]

    def point_permutation_by_generator(self, octahedral_rotation_group_matrix_implementation):
        """
        Calculates the permutation on the orbit induced by the generator elements of
        octahedral_rotation_group_implementation.

        Parameters
        ----------
        octahedral_rotation_group_matrix_implementation :
            type. A class which implements a generative representation of S4 in matrix form.
                Expected call implementation:
                    self.generators
                    self.word_to_matrix(...)

        Returns
        -------
            dict. return_value['xy'][point_id] = (Matrix Representation of 'xy')(point_id)

        """
        O = octahedral_rotation_group_matrix_implementation()
        return_value = {}
        generator_words = O.generators
        for generator in generator_words:
            mapping = {}
            for p in self:
                generator_matrix = O.word_to_matrix(generator)

                rotated_point = np.matmul(generator_matrix, self.to_cartesian(p))

                for pp in self:
                    if np.isclose(self.to_cartesian(pp), rotated_point).all():
                        mapping[p] = pp
                        break
            return_value[generator] = mapping
        return return_value

    @property
    def points(self)->list:
        """
        Gives the points of the orbit.

        Returns
        -------
            list.
        """
        return list(self._get_points())

    def __iter__(self):
        return self._get_points()

    def __len__(self):
        return len(self._point_numbers)

class LebedevOrbit_a1(_LebedevOrbitBase):
    """
    Orbit of (1, 0, 0).
    """
    _point_numbers = [str(i + 1) for i in range(6)]
    _string_id = 'a1'

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

    @staticmethod
    def _initial_state_spherical_dict():
        return {'1': (0, 90), '2': (180, 90), '3': (90, 90),
                '4': (-90, 90), '5': (90, 0), '6': (90, 180)}

    def to_cartesian(self, point):
        cart = super().to_cartesian(point)
        return tuple([int(c) for c in cart])

class LebedevOrbit_a2(_LebedevOrbitBase):
    """
    Orbit of 1/sqrt(2) * (1, 1, 0).
    """
    _point_numbers = [str(i + 1) for i in range(12)]
    _string_id = 'a2'

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

    @staticmethod
    def _initial_state_spherical_dict():
        return {'1': (90, 45),
                '2': (90, 135),
                '3': (-90, 45),
                '4': (-90, 135),
                '5': (0, 45),
                '6': (0, 135),
                '7': (180, 45),
                '8': (180, 135),
                '9': (45, 90),
                '10': (-45, 90),
                '11': (135, 90),
                '12': (-135, 90)}

class LebedevOrbit_a3(_LebedevOrbitBase):
    """
    Orbit of 1/sqrt(3) * (1, 1, 1)
    """
    _point_numbers = [str(i + 1) for i in range(8)]
    _string_id = 'a3'

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

    @staticmethod
    def _initial_state_spherical_dict():
        return {'1': (45, 54.735610317245346),
                '2': (45, 125.264389682754654),
                '3': (-45, 54.735610317245346),
                '4': (-45, 125.264389682754654),
                '5': (135, 54.735610317245346),
                '6': (135, 125.264389682754654),
                '7': (-135, 54.735610317245346),
                '8': (-135, 125.264389682754654)}

def LebedevOrbit_b_k_Meta(typical_point: tuple):
    """
    class LebedevPoints_b_k(_LebedevOrbitBase):
        pass

    return LebedevPoints_b_k
    """
    raise NotImplementedError()

def LebedevOrbit_c_k_Meta(typical_point: tuple):
    raise NotImplementedError()

def LebedevOrbit_d_k_Meta(typical_point: tuple):
    raise NotImplementedError()

# endregion

# region LebedevGrids
"""
A LebedevGrid is a union of distinct Lebedev point sets.
"""
def LebedevGridMeta(lebedev_point_sets: [_LebedevOrbitBase])->type:
    """

    Parameters
    ----------
    lebedev_point_sets :
        LebedevPointBase. The LebedevPointBase derivations from which the grid will be built.

    Returns
    -------
        type. A LebedevGrid class built from the chosen LebedevPointSets.
    """
    class LebedevGrid:
        _lebedev_point_sets_types = lebedev_point_sets

        def __init__(self):
            self._point_set_instances = {}
            for point_set_type in self._lebedev_point_sets_types:
                self._point_set_instances[point_set_type._string_id] = point_set_type()

        def point_permutation_by_generator(self, octahedral_rotation_group_matrix_implementation):
            """
            Calculates the permutation on the grid induced by the generator elements of
            octahedral_rotation_group_implementation.

            Parameters
            ----------
            octahedral_rotation_group_matrix_implementation :
                type. A class which implements a generative representation of S4 in matrix form.
                    Expected call implementation:
                        self.generators
                        self.word_to_matrix(...)

            Returns
            -------
                dict. return_value['xy'][point_id] = (Matrix Representation of 'xy')(point_id)

            """
            point_set_instances = self._point_set_instances.values()
            word_to_point_mappings = []
            for psi in point_set_instances:
                wpm = psi.point_permutation_by_generator(octahedral_rotation_group_matrix_implementation)
                word_to_point_mappings.append(wpm)

            G = OctahedralMatrixRotationGroup2Generators()
            generator_words = G.generators

            return_value = {}
            for w in generator_words:
                return_value[w] = {}

            for wpm in word_to_point_mappings:
                for w in generator_words:
                        try:
                            len_before_update = len(return_value[w])
                            return_value[w].update(wpm[w])
                            if len_before_update + len(wpm[w]) != len(return_value[w]):
                                raise AssertionError(
                                    """
                                    It seems that two orbits are not disjoint.
                                    """
                                )

                        except KeyError:
                            # If we get here than one of the point set instances does not generate a permutation
                            # dict for w
                            raise AssertionError(
                                """
                                {} seems not to act on one of the orbits.
                                """.format(w)
                            )

            return return_value

        @property
        def points(self):
            return list(self.__iter__())

        def _check_point(self, point):
            point = tuple(point)

            if point[0] not in self._point_set_instances.keys():
                raise ValueError('{} does not contain to this grid.'.format(point))

            return point

        def to_spherical(self, point):
            """
            Returns spherical coordinates of point.
            Parameters
            ----------
            point :
                tuple. Lebedev Point.

            Returns
            -------
                tuple. (theta, phi) spherical coordinates in degree of point:

                    (theta, phi) \in [-180, 180] x [0, 180]
            """
            point = self._check_point(point)
            point_class_id = point[0]
            return self._point_set_instances[point_class_id].to_spherical(point)

        def to_cartesian(self, point):
            """
            Returns cartesian coordinates of point.

            Parameters
            ----------
            point : Lebedev Point.

            Returns
            -------
                tuple. (x,y,z) cartesian coordinates of point.
            """
            point = self._check_point(point)
            point_class_id = point[0]
            return self._point_set_instances[point_class_id].to_cartesian(point)

        def __iter__(self):
            return itertools.chain(*self._point_set_instances.values())

        def __len__(self):
            return sum([len(point_set_instance) for point_set_instance in self._point_set_instances.values()])

    return LebedevGrid

LebedevGrid26 = LebedevGridMeta([LebedevOrbit_a1, LebedevOrbit_a2, LebedevOrbit_a3])

# endregion

# region Lebedev Integration
"""
implemented after https://people.sc.fsu.edu/~jburkardt/datasets/sphere_lebedev_rule/sphere_lebedev_rule.html
"""

class _Lebedev26Integrator:
    """
    Implementation for 26 point rule. Using weights from
    https://people.sc.fsu.edu/~jburkardt/datasets/sphere_lebedev_rule/lebedev_007.txt
    """
    _grid = LebedevGrid26()
    _w = {'a1': 0.047619047619048,
          'a2': 0.038095238095238,
          'a3': 0.032142857142857}

    @classmethod
    def _check_function(cls, function):
        if set(cls._grid) != function.keys():
            raise ValueError('{} is not a valid function from the LebedevGrid26 to R. ')

    @classmethod
    def weight(cls, lebedev_point):
        return cls._w[lebedev_point[0]]

    @classmethod
    def integrate(cls, function: dict):
        cls._check_function(function)

        return 4 * np.pi * sum([f_value * cls.weight(leb_point) for leb_point, f_value in function.items()])

def lebedev_26_integration(function: dict)->float:
    """
    Integration of function residing on 26 point Lebedev grid.

    Parameters
    ----------
    function :
        dict. Expects Lebedev points, i.e., ('a1', '1') ... ('a3', '8'), as keys.

    Returns
    -------
        float.

    """
    return _Lebedev26Integrator.integrate(function)

# endregion

# region Octahedral rotation group action on the set of lebedev grid functions
"""
A Lebedev Grid Function, f, is a dict with keys in a LebedevGrid.
The octahedral rotation group, O, induces an action, sigma, M = {f: f is a Lebedev Grid Function}
by setting

    sigma(f, rho) := f(rho(.)).

This region implements:
    1. A representation of O on Lebedev point classes.

    2. sigma
"""

class ActionOctahedralRotationGroupOnLebedevGridFunctions:
    def __init__(self, lebedev_grid: type, octahedral_rotation_group_matrix_implementation: type):
        self._lebedev_orbit_instance = lebedev_grid()
        self._group_instance = octahedral_rotation_group_matrix_implementation()
        self._recursion_tails = self._generate_recursion_tails(octahedral_rotation_group_matrix_implementation)

    def _generate_recursion_tails(self,
                                  octahedral_rotation_group_matrix_implementation):
        grid = self._lebedev_orbit_instance
        permutation_by_generator = grid.point_permutation_by_generator(octahedral_rotation_group_matrix_implementation)

        recursion_tails = {}

        def get_recursion_tail_for_generator_word(permut):
            def recursion_tail_for_generator_word(f):
                f_new = {}
                for k, v in permut.items():
                    f_new[v] = f[k]

                return f_new
            return recursion_tail_for_generator_word

        for word, permutation in permutation_by_generator.items():
            recursion_tails[word] = get_recursion_tail_for_generator_word(permutation)

        return recursion_tails

    def _check_f(self, f):
        if set(f.keys()) != set(self._lebedev_orbit_instance.points):
            raise ValueError(
             """
             {} is not a function from a Lebedev Grid.
             """.format(f)
            )

    def _check_w(self, word):
        generators = self._group_instance.generators

        for c in word:
            if c not in generators:
                raise ValueError(
                    """
                    {} does not contain to {}. It is no combination of its generators.
                    """.format(word, type(self._group_instance))
                )

    def _execute(self, f, w):
        """
        This is a right action on the set of functions from the lebedev grid/orbit.
        """
        if w == '':
            return f

        if w in self._recursion_tails:
            return self._recursion_tails[w](f)
        else:
            return self._execute(self._execute(f, w[0]), w[1:])

    def __call__(self, f: dict, word: str):
        self._check_f(f)
        self._check_w(word)
        return self._execute(f, word)

# endregion


text:
import config
from db import text_search, items
from utils import resultproxy_to_dict
import os
from flask_cors import CORS, cross_origin
import insightface
import urllib
import urllib.request
import cv2
import json
import numpy as np
import uuid
import os
from face_searcher import face_searcher
from flask import Flask, request, render_template, send_file

if not os.path.exists(config.images_dir):
    os.makedirs(config.images_dir)

print("Preparing model")
model = insightface.app.FaceAnalysis(
    rec_name="arcface_r100_v1", det_name="retinaface_mnet025_v2"
)
ctx_id = -1
model.prepare(ctx_id=ctx_id, nms=0.4)
searcher = face_searcher()
app = Flask(__name__)

@app.route("/images/<path:path>")
@cross_origin()
def images(path):
    print(path)
    return send_file(os.path.join(config.images_dir, path), mimetype="image/JPG")

def url_to_image(url):
    resp = urllib.request.urlopen(url)
    image = np.asarray(bytearray(resp.read()), dtype="uint8")
    image = cv2.imdecode(image, cv2.IMREAD_COLOR)
    return image

@app.route("/people/search")
@cross_origin()
def search_people():
    name = request.args.get("name")
    people = resultproxy_to_dict(text_search.get_people_by_name(name))
    result = {"people": people}
    return json.dumps(result)

@app.route("/people/<id>")
@cross_origin()
def get_person(id):
    person = resultproxy_to_dict(items.get_person_by_id(id))
    films = resultproxy_to_dict(items.get_related_films(id))
    result = {"person": person, "films": films}
    return json.dumps(result)

@app.route("/films/search")
@cross_origin()
def search_films():
    title = request.args.get("title")
    films = resultproxy_to_dict(text_search.get_films_by_title(title))
    result = {"films": films}
    return json.dumps(result)

@app.route("/films/<id>")
@cross_origin()
def get_film_info(id):
    film = resultproxy_to_dict(items.get_film_by_id(id))
    people = resultproxy_to_dict(items.get_related_people(id))
    result = {"film": film, "people": people}
    return json.dumps(result)

@app.route("/actors_recognition", methods=["POST"])
@cross_origin()
def recognize_actors():
    image_undecoded = request.get_data()
    img = cv2.imdecode(np.frombuffer(image_undecoded, np.uint8), -1)
    img_out = img.copy()
    faces = model.get(img)
    faces_arr = []
    img_name = str(uuid.uuid4())
    for idx, face in enumerate(faces):
        face_d = {}
        face_d["index"] = idx
        face_d["age"] = face.age
        gender = "Male"
        if face.gender == 0:
            gender = "Female"
        face_d["gender"] = gender
        face_search_res = searcher.search_face(face.embedding, 1)
        face_d["search_results"] = face_search_res
        bboxes = face.bbox.astype(np.int)
        face_name = img_name + "_" + str(idx) + ".jpg"
        face_d["face_image_url"] = "http://127.0.0.1:5000/images/" + face_name
        cv2.imwrite(
            "images/" + face_name,
            img[int(bboxes[1]) : int(bboxes[3]), int(bboxes[0]) : int(bboxes[2])],
        )
        face_d["bbox"] = str(bboxes.flatten())
        img_out = cv2.rectangle(
            img_out, (bboxes[0], bboxes[1]), (bboxes[2], bboxes[3]), (0, 255, 0), 2
        )
        img_out = cv2.putText(
            img_out,
            str(face.age),
            (int(bboxes[0] - 10), int(bboxes[1] - 10)),
            cv2.FONT_HERSHEY_SIMPLEX,
            1.1,
            (0, 255, 0),
            2,
        )
        img_out = cv2.putText(
            img_out,
            str(gender),
            (int(bboxes[0] + 40), int(bboxes[1] - 10)),
            cv2.FONT_HERSHEY_SIMPLEX,
            1.1,
            (0, 0, 255),
            2,
        )
        img_out = cv2.putText(
            img_out,
            str(idx),
            (int(bboxes[0] - 40), int(bboxes[1] - 10)),
            cv2.FONT_HERSHEY_SIMPLEX,
            1.1,
            (0, 0, 255),
            2,
        )
        faces_arr.append(face_d)
    cv2.imwrite(os.path.join(config.images_dir, img_name + ".jpg"), img_out)
    resp = {
        "faces": faces_arr,
        "result_image": "http://127.0.0.1:5000/images/" + img_name + ".jpg",
    }
    return json.dumps(resp, ensure_ascii=False), 200

if __name__ == "__main__":
    print("Starting server on port 5000")
    app.run(host="0.0.0.0", port=5000)

max_stars_repo_name: svetinsideu/monorepo

max_issues_repo_name: svetinsideu/monorepo

max_forks_repo_name: svetinsideu/monorepo


text:
""" Amazon Access Challenge Code for ensemble

Marios Michaildis script for Amazon .

xgboost on input data

based on Paul Duan's Script.

"""
from __future_
_ import division
import numpy as np
from sklearn import  preprocessing
from sklearn.metrics import roc_auc_score
import XGBoostClassifier as xg
from sklearn.cross_validation import StratifiedKFold

SEED = 42  # always use a seed for randomized procedures

def load_data(filename, use_labels=True):
    """
    Load data from CSV files and return them as numpy arrays
    The use_labels parameter indicates whether one should
    read the first column (containing class labels). If false,
    return all 0s.
    """

    # load column 1 to 8 (ignore last one)
    data = np.loadtxt(open( filename), delimiter=',',
                      usecols=range(1, 9), skiprows=1)
    if use_labels:
        labels = np.loadtxt(open( filename), delimiter=',',
                            usecols=[0], skiprows=1)
    else:
        labels = np.zeros(data.shape[0])
    return labels, data

def save_results(predictions, filename):
    """Given a vector of predictions, save results in CSV format."""
    with open(filename, 'w') as f:
        f.write("id,ACTION\n")
        for i, pred in enumerate(predictions):
            f.write("%d,%f\n" % (i + 1, pred))

def bagged_set(X_t,y_c,model, seed, estimators, xt, update_seed=True):

   # create array object to hold predictions
   baggedpred=[ 0.0  for d in range(0, (xt.shape[0]))]
   #loop for as many times as we want bags
   for n in range (0, estimators):
        #shuff;e first, aids in increasing variance and forces different results
        #X_t,y_c=shuffle(Xs,ys, random_state=seed+n)

        if update_seed: # update seed if requested, to give a slightly different model
            model.set_params(random_state=seed + n)
        model.fit(X_t,y_c) # fit model0.0917411475506
        preds=model.predict_proba(xt)[:,1] # predict probabilities
        # update bag's array
        for j in range (0, (xt.shape[0])):
                baggedpred[j]+=preds[j]
   # divide with number of bags to create an average estimate
   for j in range (0, len(baggedpred)):
                baggedpred[j]/=float(estimators)
   # return probabilities
   return np.array(baggedpred)

# using numpy to print results
def printfilcsve(X, filename):

    np.savetxt(filename,X)

def main():
    """
    Fit models and make predictions.
    We'll use one-hot encoding to transform our categorical features
    into binary features.
    y and X will be numpy array objects.
    """

    filename="main_xgboost" # nam prefix
    #model = linear_model.LogisticRegression(C=3)  # the classifier we'll use

    model=xg.XGBoostClassifier(num_round=1000 ,nthread=25,  eta=0.12, gamma=0.01,max_depth=12, min_child_weight=0.01, subsample=0.6,
                                   colsample_bytree=0.7,objective='binary:logistic',seed=1)
    # === load data in memory === #
    print "loading data"
    y, X = load_data('train.csv')
    y_test, X_test = load_data('test.csv', use_labels=False)

    # === one-hot encoding === #
    # we want to encode the category IDs encountered both in
    # the training and the test set, so we fit the encoder on both
    encoder = preprocessing.OneHotEncoder()
    encoder.fit(np.vstack((X, X_test)))
    X = encoder.transform(X)  # Returns a sparse matrix (see numpy.sparse)
    X_test = encoder.transform(X_test)

    # if you want to create new features, you'll need to compute them
    # before the encoding, and append them to your dataset after

    #create arrays to hold cv an dtest predictions
    train_stacker=[ 0.0  for k in range (0,(X.shape[0])) ]

    # === training & metrics === #
    mean_auc = 0.0
    bagging=20 # number of models trained with different seeds
    n = 5  # number of folds in strattified cv
    kfolder=StratifiedKFold(y, n_folds= n,shuffle=True, random_state=SEED)
    i=0
    for train_index, test_index in kfolder: # for each train and test pair of indices in the kfolder object
        # creaning and validation sets
        X_train, X_cv = X[train_index], X[test_index]
        y_train, y_cv = np.array(y)[train_index], np.array(y)[test_index]
        #print (" train size: %d. test size: %d, cols: %d " % ((X_train.shape[0]) ,(X_cv.shape[0]) ,(X_train.shape[1]) ))

        # if you want to perform feature selection / hyperparameter
        # optimization, this is where you want to do it

        # train model and make predictions
        preds=bagged_set(X_train,y_train,model, SEED , bagging, X_cv, update_seed=True)

        # compute AUC metric for this CV fold
        roc_auc = roc_auc_score(y_cv, preds)
        print "AUC (fold %d/%d): %f" % (i + 1, n, roc_auc)
        mean_auc += roc_auc

        no=0
        for real_index in test_index:
                 train_stacker[real_index]=(preds[no])
                 no+=1
        i+=1

    mean_auc/=n
    print (" Average AUC: %f" % (mean_auc) )
    print (" printing train datasets ")
    printfilcsve(np.array(train_stacker), filename + ".train.csv")

    # === Predictions === #
    # When making predictions, retrain the model on the whole training set
    preds=bagged_set(X, y,model, SEED, bagging, X_test, update_seed=True)

    #create submission file
    printfilcsve(np.array(preds), filename+ ".test.csv")
    #save_results(preds, filename+"_submission_" +str(mean_auc) + ".csv")

if __name__ == '__main__':
    main()


text:
# Z-test for model compare

# input:
#          | model1 | model2 | ...
# dataset1 |  acc   |  acc   |
# dataset2 | ...
# ...

# output:
# model compare crosstable

import numpy as np
import math

def Z_test(acc, datan, threshold=1.96):
    n=acc.shape[0]
    result=np.zeros([n,n,3]) # win-loss-draw
    for i in range(n):
        for j in range(n):
            for k in range(acc.shape[1]):
                pa=acc[i,k]
                pb=acc[j,k]
                p=(pa+pb)/2
                Z=(pa-pb)/math.sqrt(2*p*(1-p)/datan[k])
                if Z>threshold:
                    result[i,j,0]+=1
                elif Z<-threshold:
                    result[i,j,1]+=1
                else:
                    result[i,j,2]+=1
            print("%d-%d-%d "%tuple(result[i,j,:]),end='')
        print()
    return result

if __name__=="__main__":
    acc=[
    [92.09,79.62,87.19],
    [85.51,76.81,84.78],
    [81.95,58.05,70.73],
    [95.14,95.99,96.42],
    [76.24,83.50,84.49],
    [85.80,77.54,85.07],
    [72.40,75.91,76.82],
    [70.90,74.70,74.40],
    [67.29,48.59,59.81],
    [80.00,84.07,83.70],
    [81.94,83.23,87.10],
    [85.33,78.80,72.61],
    [85.33,78.80,82.61],
    [94.67,95.33,96.00],
    [78.95,94.74,92.98],
    [73.34,73.16,73.56],
    [77.03,83.11,86.49],
    [74.35,76.04,76.95],
    [78.85,69.71,76.92],
    [83.72,70.04,98.33],
    [71.04,45.04,74.94],
    [94.38,96.63,98.88],
    [93.07,93.07,96.04],
    ]
    datan=[898,690,205,699,303,690,768,1000,214,270,155,368,351,150,57,3200,148,768,208,958,846,178,101]
    ret=Z_test(np.array(acc).T/100, np.array(datan))


text:
# Create fake income/age clusters for N people in k clusters
# generate random integer values
from random import seed
from random import randint
import numpy as np
import decimal

f = open("LaptopPrice.csv", "w")
f.write("Price,Speed,Ram,HDD,Class\n")
Class = ['Personal', 'Business', 'Gaming']

def createClusteredData(N, k):
    pointsPerCluster = float(N)/k
    X = []
    y = []
    for i in range(k):
        Speed = i*5/k + 1 + np.random.uniform(1, .1)
        # np.random.uniform(1000, 32000)
        Ram = i*4000/k + 2000 + np.random.uniform(0, 500)
        HDD = np.random.uniform(1000, 2000)

        for j in range(int(pointsPerCluster)):
            SpeedI = round(np.random.normal(Speed, .7), 2)
            RamI = round(np.random.normal(Ram, 600))
            HDDI = round(np.random.normal(HDD, 300))
            PriceS = SpeedI*75
            PriceR = RamI/25
            PriceH = HDDI/15
            print('PriceS ' + str(PriceS) + ' PriceR ' +
                  str(PriceR) + ' PriceH ' + str(PriceH))
            Price = PriceH + PriceR + PriceS
            a = str(round(np.random.uniform(Price, 100), 2)) + ","\
                + str(SpeedI) + ","\
                + str(RamI) + ","\
                + str(HDDI) + "," \
                + str(Class[i]) + "\n"
            # print(a)
            f.write(a)
    return
createClusteredData(100, 3)

max_stars_repo_name: lavishthomas/MLexamples

max_issues_repo_name: lavishthomas/MLexamples

max_forks_repo_name: lavishthomas/MLexamples


text:
import numpy as np
import torch

class CTRMF(torch.nn.Module):
    def __init__(self,
                 ratings_data,
                 recipe_feature,
                 n_hidden=80,
                 reg_term=0.01,
                 device=None,
                 verbose=False):
        """
                ratings_data : ((n_users, n_recipes) 2D array)
                        Collection of recipe ratings for each user.
                recipe_feature : ((n_recipes, n_features) 2D array)
                        Recipe-feature matrix, where each element is 1 if the recipe
                        contains the corresponding feature, 0 otherwise.
                n_hidden : (Integer)
                        Number of latent dimensions
                reg_term : (Double)
                        Regularization term
                verbose : (Boolean)
                        Prints helpful training progress messages if True
                """
        super().__init__()

        self.ratings_data = ratings_data
        self.recipe_feature = torch.FloatTensor(recipe_feature).to(device)
        self.n_hidden = n_hidden
        self.n_users, self.n_recipes = ratings_data.shape
        self.n_features = recipe_feature.shape[1]
        self.reg_term = reg_term
        self.verbose = verbose

        # Set learned matrices as Embeddings
        self.user_hidden = torch.nn.Embedding(
            self.n_users,
            self.n_hidden
        ).to(device)
        self.hidden_feature = torch.nn.Embedding(
            self.n_hidden,
            self.n_features
        ).to(device)
        self.user_biases = torch.nn.Embedding(self.n_users, 1).to(device)
        self.recipe_biases = torch.nn.Embedding(
            self.n_recipes, 1).to(device)

        # Initialize learned matrices
        torch.nn.init.xavier_uniform_(self.user_hidden.weight)
        torch.nn.init.xavier_uniform_(self.hidden_feature.weight)
        self.user_biases.weight.data.fill_(0.)
        self.recipe_biases.weight.data.fill_(0.)

        # Calculate the mean of the ratings data
        self.global_bias = torch.FloatTensor([np.mean(
            self.ratings_data[np.where(self.ratings_data != 0)]
        )]).to(device)

        if self.verbose:
            print("Initializing ctrmf")
            print("==================")
            print(f"n_users: {self.n_users}")
            print(f"n_recipes: {self.n_recipes}")
            print(f"recipe_feature.shape: {recipe_feature.shape}")
            print(f"n_features: {self.n_features}")
            print(f"n_hidden: {self.n_hidden}")
            print(
                f"user_hidden.shape: ({self.user_hidden.num_embeddings},{self.user_hidden.embedding_dim})")
            print(
                f"hidden_feature.shape: {self.hidden_feature.weight.shape}")
            print('\n')

    def forward(self, user, recipe):
        pred = torch.matmul(self.user_hidden(user), self.hidden_feature.weight)
        pred = torch.matmul(pred, self.recipe_feature[recipe].T)
        pred += self.user_biases(user) + \
            self.recipe_biases(recipe) + self.global_bias
        return pred

max_stars_repo_name: NathanRiviere/Recipe-Recommender

max_issues_repo_name: NathanRiviere/Recipe-Recommender

max_forks_repo_name: NathanRiviere/Recipe-Recommender


text:
import torch

import numpy as np
import sklearn.neighbors
import time

def get_knn_idx_dist(pos:torch.FloatTensor, query:torch.FloatTensor, k, offset=0):
    """
    :param  pos:     (B, N, F)
    :param  query:   (B, M, F)
    :return knn_idx: (B, M, k)
    """
    B, N, F = tuple(pos.size())
    M = query.size(1)

    pos = pos.unsqueeze(1).expand(B, M, N, F)
    query  = query.unsqueeze(2).expand(B, M, N, F)   # B * M * N * F
    dist = torch.sum((pos - query) ** 2, dim=3, keepdim=False)   # B * M * N
    knn_idx = torch.argsort(dist, dim=2)[:, :, offset:k+offset]   # B * M * k
    knn_dist = torch.gather(dist, dim=2, index=knn_idx)           # B * M * k

    return knn_idx, knn_dist

def get_knn_idx(pos, query, k, offset=0):
    """
    :param  pos:     (B, N, F)
    :param  query:   (B, M, F)
    :return knn_idx: (B, M, k)
    """
    knn_idx, _ = get_knn_idx_dist(pos=pos, query=query, k=k, offset=offset)

    return knn_idx

def group(x:torch.FloatTensor, idx:torch.LongTensor):
    """
    :param  x:      (B, N, F)
    :param  idx:    (B, M, k)
    :return (B, M, k, F)
    """
    B, N, F = tuple(x.size())
    _, M, k = tuple(idx.size())

    x = x.unsqueeze(1).expand(B, M, N, F)
    idx = idx.unsqueeze(3).expand(B, M, k, F)

    return torch.gather(x, dim=2, index=idx)

def gather(x:torch.FloatTensor, idx:torch.LongTensor):
    """
    :param  x:      (B, N, F)
    :param  idx:    (B, M)
    :return (B, M, F)
    """
    # x       : B * N * F
    # idx     : B * M
    # returns : B * M * F
    B, N, F = tuple(x.size())
    _, M    = tuple(idx.size())

    idx = idx.unsqueeze(2).expand(B, M, F)

    return torch.gather(x, dim=1, index=idx)

def feature_interp(k, feat, pos, pos_new, avg_dist, feat_new=None, avg_feat_diff=None):
    """
    :param  feat:     (B, N, F)
    :param  pos:      (B, N, 3)
    :param  pos_new:  (B, M, 3)
    :param  feat_new: (B, M, F)
    :return (B, M, F)
    """
    knn_idx = get_knn_idx(pos, pos_new, k=k, offset=0)
    pos_grouped = group(pos, idx=knn_idx)    # (B, M, k, 3)
    feat_grouped = group(feat, idx=knn_idx)  # (B, M, k, F)

    d_pos = ((pos_grouped - pos_new.unsqueeze(-2).expand_as(pos_grouped)) ** 2).sum(dim=-1)     # (B, M, k)
    weight = - d_pos / (avg_dist ** 2)

    if feat_new is not None:
        d_feat = ((feat_grouped - feat_new.unsqueeze(-2).expand_as(feat_grouped)) ** 2).sum(dim=-1) # (B, M, k)
        weight = weight - d_feat / (avg_feat_diff ** 2)

    weight = weight.softmax(dim=-1)   # (B, M, k)
    return (feat_grouped * weight.unsqueeze(-1).expand_as(feat_grouped)).sum(dim=-2)

def get_1d_mesh(steps, start=-0.2, end=0.2):
    return torch.linspace(start=start, end=end, steps=steps).unsqueeze(-1)

def get_2d_mesh(steps, start=-0.2, end=0.2):
    mesh_1d = get_1d_mesh(steps=steps, start=start, end=end).flatten()
    return torch.cartesian_prod(mesh_1d, mesh_1d)

def get_mesh(dim, steps, start=-0.2, end=0.2):
    assert dim in (1, 2)
    if dim == 1:
        return get_1d_mesh(steps, start=start, end=end)
    elif dim == 2:
        return get_2d_mesh(steps, start=start, end=end)

def get_sample_points(dim, samples, num_points=1, num_batch=None, start=-0.3, end=0.3):
    length = end - start
    if num_batch is None:
        size = [samples * num_points, dim]
    else:
        size = [num_batch, samples * num_points, dim]
    return (torch.rand(size) * length) - (length / 2)

if __name__ == '__main__':
    nbatch = 1
    npts = 5
    mesh_dim = 2
    bth = torch.arange(0, nbatch * npts * 3).float().reshape([nbatch, npts, 3]).repeat(1, 4, 1)
    print(bth)
    mesh = get_2d_mesh(2)
    print(mesh)

    mesh = mesh.unsqueeze(-1).repeat(1, 1, npts).transpose(1, 2).reshape(1, -1, mesh_dim).repeat(nbatch, 1, 1)
    print(mesh.size())
    # expanded = torch.cat([bth, mesh], dim=-1)
    # print(torch.cat([bth, mesh], dim=-1))

class Aggregator(torch.nn.Module):

    def __init__(self, oper):
        super().__init__()
        assert oper in ('mean', 'sum', 'max')
        self.oper = oper

    def forward(self, x, dim=2):
        if self.oper == 'mean':
            return x.mean(dim=dim, keepdim=False)
        elif self.oper == 'sum':
            return x.sum(dim=dim, keepdim=False)
        elif self.oper == 'max':
            ret, _ = x.max(dim=dim, keepdim=False)
            return ret

class FullyConnected(torch.nn.Module):

    def __init__(self, in_features, out_features, bias=True, activation=None):
        super().__init__()

        self.linear = torch.nn.Linear(in_features, out_features, bias=bias)

        if activation is None:
            self.activation = torch.nn.Identity()
        elif activation == 'relu':
            self.activation = torch.nn.ReLU()
        elif activation == 'elu':
            self.activation = torch.nn.ELU(alpha=1.0)
        elif activation == 'lrelu':
            self.activation = torch.nn.LeakyReLU(0.1)
        else:
            raise ValueError()

    def forward(self, x):
        return self.activation(self.linear(x))

def normalize_point_cloud(pc, center=None, scale=None):
    """
    :param  pc: (B, N, 3)
    :return (B, N, 3)
    """
    if center is None:
        center = torch.mean(pc, dim=-2, keepdim=True).expand_as(pc)
    if scale is None:
        scale, _ = torch.max(pc.reshape(pc.size(0), -1).abs(), dim=1, keepdim=True)
        scale = scale.unsqueeze(-1).expand_as(pc)
    norm = (pc - center) / scale
    return norm, center, scale

def denormalize_point_cloud(pc, center, scale):
    """
    :param  pc: (B, N, 3)
    :return (B, N, 3)
    """
    return pc * scale + center

max_stars_repo_name: RRemixx/DMRDenoise

max_issues_repo_name: LifeBeyondExpectations/DMRDenoise

max_forks_repo_name: LifeBeyondExpectations/DMRDenoise


text:
import theano
import theano.tensor as T
from theano.sandbox.rng_mrg import MRG_RandomStreams
from theano.tensor.nnet.conv import conv2d
from theano.tensor.signal.downsample import max_pool_2d
from theano.tensor.shared_randomstreams import RandomStreams

import numpy as np
import scipy.io
import time
import sys
import logging
import copy

from toolbox import *
from modelbase import *

class Draw_sgru1(ModelULBase):
    def __init__(self, data, hp):
        super(Draw_sgru1, self).__init__(self.__class__.__name__, data, hp)

        n_x = self.data['n_x']
        n_h = 256
        n_t = 12
        n_zpt = 32
        gates = 4
        self.n_t = n_t
        self.n_z = n_t * n_zpt
        self.sample_steps = False

        self.params = Parameters()

        if hp.load_model and os.path.isfile(self.filename):
            self.params.load(self.filename)
        else:
            with self.params:
                W1 = shared_normal((gates*1, n_x + n_h, n_h), scale=hp.init_scale)
                W11 = shared_normal((gates*1, n_h, n_h), scale=hp.init_scale)
                W4 = shared_normal((gates*1, n_zpt, n_h), scale=hp.init_scale)
                W44 = shared_normal((gates*1, n_h, n_h), scale=hp.init_scale)
                b1 = shared_zeros((gates*1, n_h,))
                b4 = shared_zeros((gates*1, n_h,))
                b10_h = shared_zeros((n_h,))
                b40_h = shared_zeros((n_h,))

                W2 = shared_normal((n_h, n_zpt), scale=hp.init_scale)
                W3 = shared_normal((n_h, n_zpt), scale=hp.init_scale)
                W5 = shared_normal((n_h, n_x), scale=hp.init_scale)
                b2 = shared_zeros((n_zpt,))
                b3 = shared_zeros((n_zpt,))
                b5 = shared_zeros((n_x,))
                ex0 = shared_zeros((n_x,))

        def rnn(X, h, W, U, b, t):
            return T.tanh(T.dot(X,W[0]) + T.dot(h,U[0]) + b[0])

        #def gru(X, h, W, U, b):
        #    z_t = T.nnet.sigmoid(T.dot(X,W[:,:n_h]) + T.dot(h,U[:,:n_h]) + b[n_h])
        #    r_t = T.nnet.sigmoid(T.dot(X,W[:,n_h:2*n_h]) + T.dot(h,U[:,n_h:2*n_h]) + b[n_h:2*n_h])
        #    h_t = T.tanh(T.dot(X,W[:,2*n_h:3*n_h]) + r_t * T.dot(h,U[:,2*n_h:3*n_h]) + b[2*n_h:3*n_h])
        #    return (1 - z_t) * h + z_t * h_t

        def sgru(X, h, W, U, b, t):
            t = 0
            z_t = T.tanh(T.dot(X,W[t*2+0]) + T.dot(h,U[t*2+0]) + b[t*2+0])
            r_t = T.tanh(T.dot(X,W[t*2+1]) + T.dot(h,U[t*2+1]) + b[t*2+1])
            return z_t * r_t

        def sgru2(X, h, W, U, b, t):
            t = 0
            z_t = T.tanh(T.dot(X,W[t*2+0]) + b[t*2+0])
            r_t = (T.dot(h,U[t*2+0]) + b[t*2+1])
            return T.tanh(T.dot(z_t*r_t,T.transpose(U[t*2+1])))

        def sgru3(X, h, W, U, b, t):
            t = 0
            z_t = T.tanh(T.dot(X,W[t*2+0]) + b[t*2+0])
            r_t = (T.dot(h,U[t*2+0]) + b[t*2+1])
            z_t2 = (T.dot(X,W[t*2+1]) + b[t*2+2])
            r_t2 = T.tanh(T.dot(h,U[t*2+1]) + b[t*2+3])
            return T.tanh(T.dot(z_t*r_t,T.transpose(U[t*2+2])) + T.dot(z_t2*r_t2,T.transpose(U[t*2+3])))

        frnn = sgru

        # Encoder
        p = self.params

        x = self.X  # binomial(self.X)
        input_size = x.shape[0]
        ex = p.ex0
        h_encoder_h = p.b10_h
        h_decoder_h = T.zeros((input_size, p.b40_h.shape[0])) + p.b40_h
        log_qpz = 0

        for t in xrange(0, n_t):
            x_e = x - T.nnet.sigmoid(ex)
            h_x = concatenate([x_e, h_decoder_h], axis=1)
            h_encoder_h = frnn(h_x, h_encoder_h, p.W1, p.W11, p.b1, t)
            mu_encoder_t = T.dot(h_encoder_h, p.W2) + p.b2
            log_sigma_encoder_t = 0.5*(T.dot(h_encoder_h, p.W3) + p.b3)
            log_qpz += -0.5* T.sum(1 + 2*log_sigma_encoder_t - mu_encoder_t**2 - T.exp(2*log_sigma_encoder_t))

            eps = srnd.normal(mu_encoder_t.shape, dtype=theano.config.floatX)
            z = mu_encoder_t + eps*T.exp(log_sigma_encoder_t)
            h_decoder_h = frnn(z, h_decoder_h, p.W4, p.W44, p.b4, t)
            ex += T.dot(h_decoder_h, p.W5) + p.b5

        pxz = T.nnet.sigmoid(ex)
        log_pxz = T.nnet.binary_crossentropy(pxz, x).sum()
        cost = log_pxz + log_qpz

        # Generate
        z = self.Z.reshape((-1, n_t, n_zpt), ndim=3)
        input_size = z.shape[0]
        s_ex = p.ex0
        s_h_decoder_h = p.b40_h

        if self.sample_steps:
            a_pxz = T.zeros((n_t + 1, input_size, n_x))
        else:
            a_pxz = T.zeros((1, input_size, n_x))

        for t in xrange(0, n_t):
            if self.sample_steps:
                a_pxz = T.set_subtensor(a_pxz[t,:,:], T.nnet.sigmoid(s_ex))
            s_h_decoder_h = frnn(z[:,t,:], s_h_decoder_h, p.W4, p.W44, p.b4, t)
            s_ex +=  T.dot(s_h_decoder_h, p.W5) + p.b5

        a_pxz = T.set_subtensor(a_pxz[-1,:,:], T.nnet.sigmoid(s_ex))

        self.compile(log_pxz, log_qpz, cost, a_pxz)

max_stars_repo_name: MilesQLi/Theano-Lights

max_issues_repo_name: MilesQLi/Theano-Lights

max_forks_repo_name: Ivaylo-Popov/Theano-Lights

reason: import numpy,import scipy,import theano,from theano


text:
## Automatically adapted for numpy.oldnumeric Apr 14, 2008 by -c

import sys
import os
import Gnuplot
import numpy as numpy
"""
Script to creat an animation of Osher solution dat file
"""
print "Making animations"

datFile = sys.argv[1]
totalNumberOfIndex = int(sys.argv[2])
T_start = float(sys.argv[3])
T_stop  = float(sys.argv[4])

FudgeFactor = 1e-5
i = 0

while (i <= totalNumberOfIndex):

    time  = T_start + i*( (T_stop-T_start)/totalNumberOfIndex )

    g = Gnuplot.Gnuplot(debug=1)
    g('set style data lines')
    g('set key bottom left')
    g('set yrange [-0.0:1.0]')
    g.xlabel('x')
    g.ylabel('S_e')
    PlotTitle = "S_e vs x at time: " + str(time)
    g.title(PlotTitle)
    g('set term postscript eps enhanced color solid')
    Fraction = float(time)/float(T_stop) + FudgeFactor
    outFile = 'set output \"' + str('Se') + 'plotFraction' + str(Fraction) + 'MakeAnimations.eps\"'
    g(outFile)
    Plot = 'plot \'' + str(datFile) + '\' index ' + str(i) + ' title \"\"'
    g(Plot)

    i += 1

i = 0
while(i < 1000000):
    # Counting delay
    i = i + 1

Convert = 'convert -loop 1 *MakeAnimations.eps OsherAnimation.gif'
print Convert
os.system(Convert)

max_stars_repo_name: robertsawko/proteus

max_issues_repo_name: robertsawko/proteus

max_forks_repo_name: robertsawko/proteus


text:
import numpy as numpy

a = numpy.array([[0.0, 0.0, 0.0],[10.0, 10.0, 10.0],[20.0 ,20.0, 20.0],[30.0, 30.0, 30.0]])
b = numpy.array([1.0,2.0,3.0])

print ('First array:')
print (a)
print ('\n')

print ('Second array:')
print (b)
print ('\n')

print ('First Array + Second Array')
print (a + b)


text:
"""Training Neural A*
Author: Ryo Yonetani
Affiliation: OSX
"""

import argparse
from datetime import datetime
import subprocess

import numpy as np
from tqdm import tqdm
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.tensorboard import SummaryWriter

from neural_astar.planner import NeuralAstar, VanillaAstar
from neural_astar.utils.data import create_dataloader
from neural_astar.utils.training import (set_global_seeds, run_planner,
                                         calc_metrics, visualize_results,
                                         Metrics)

def main():

    # argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("--data_npz",
                        "-d",
                        type=str,
                        default="data/mazes_032_moore_c8.npz")
    parser.add_argument("--logdir", "-l", type=str, default="log")
    parser.add_argument("--encoder-arch",
                        "-e",
                        type=str,
                        default="CNN",
                        choices=["CNN", "Unet"])
    parser.add_argument("--batch-size", "-b", type=int, default=100)
    parser.add_argument("--num-epochs", "-n", type=int, default=50)
    parser.add_argument("--seed", "-s", type=int, default=1234)
    args = parser.parse_args()

    # dataloaders
    set_global_seeds(args.seed)
    train_loader = create_dataloader(args.data_npz,
                                     "train",
                                     args.batch_size,
                                     shuffle=True)
    val_loader = create_dataloader(args.data_npz,
                                   "valid",
                                   args.batch_size,
                                   shuffle=False)
    test_loader = create_dataloader(args.data_npz,
                                    "test",
                                    args.batch_size,
                                    shuffle=False)

    # planners
    device = "cuda" if torch.cuda.is_available() else "cpu"
    neural_astar = NeuralAstar(encoder_arch=args.encoder_arch)
    neural_astar.to(device)
    vanilla_astar = VanillaAstar()
    vanilla_astar.to(device)

    # training setup
    opt = optim.RMSprop(neural_astar.parameters(), lr=0.001)
    criterion = nn.L1Loss()

    # logger setup
    now = str(datetime.now())
    logdir = f"{args.logdir}/{now}"
    writer = SummaryWriter(f"{logdir}/tb")
    h_mean_best = -1.0

    for e in range(args.num_epochs):
        train_loss, val_loss, p_opt, p_exp, h_mean = 0., 0., 0., 0., 0.

        # training
        for batch in tqdm(train_loader, desc="training", ncols=60):
            neural_astar.train()
            loss, na_outputs = run_planner(batch, neural_astar, criterion)
            train_loss += loss.item()
            opt.zero_grad()
            loss.backward()
            opt.step()
        train_loss /= len(train_loader)

        # validation
        with torch.no_grad():
            for batch in tqdm(val_loader, desc="validation", ncols=60):
                neural_astar.eval()
                loss, na_outputs = run_planner(batch, neural_astar, criterion)
                _, va_outputs = run_planner(batch, vanilla_astar, criterion)
                metrics = calc_metrics(na_outputs, va_outputs)
                val_loss += loss
                p_opt += metrics.p_opt
                p_exp += metrics.p_exp
                h_mean += metrics.h_mean
        val_loss /= len(val_loader)
        p_opt /= len(val_loader)
        p_exp /= len(val_loader)
        h_mean /= len(val_loader)

        # logging
        print(
            f"[epoch:{e:03d}] train_loss:{train_loss:.2e}, val_loss:{val_loss:.2e}, ",
            Metrics(p_opt, p_exp, h_mean))

        writer.add_scalar("metrics/train_loss", train_loss, e)
        writer.add_scalar("metrics/val_loss", val_loss, e)
        writer.add_scalar("metrics/p_opt", p_opt, e)
        writer.add_scalar("metrics/p_exp", p_exp, e)
        writer.add_scalar("metrics/h_mean", h_mean, e)

        va_results = visualize_results(batch[0], va_outputs)
        na_results = visualize_results(batch[0], na_outputs)
        writer.add_image("vis/astar", va_results, e, dataformats="HWC")
        writer.add_image("vis/neural-astar", na_results, e, dataformats="HWC")

        # checkpointing
        if h_mean > h_mean_best:
            print(f"best score updated: {h_mean_best:0.3f} -> {h_mean:0.3f}")
            h_mean_best = h_mean
            subprocess.run(["rm", "-rf", f"{logdir}/best.pt"])
            torch.save(neural_astar.state_dict(), f"{logdir}/best.pt")
    writer.close()

    # testing
    neural_astar.load_state_dict(torch.load(f"{logdir}/best.pt"))
    p_opt, p_exp, h_mean = 0., 0., 0.
    with torch.no_grad():
        for batch in tqdm(test_loader, desc="test", ncols=60):
            neural_astar.eval()
            loss, na_outputs = run_planner(batch, neural_astar, criterion)
            _, va_outputs = run_planner(batch, vanilla_astar, criterion)
            metrics = calc_metrics(na_outputs, va_outputs)
            p_opt += metrics.p_opt
            p_exp += metrics.p_exp
            h_mean += metrics.h_mean
        p_opt /= len(test_loader)
        p_exp /= len(test_loader)
        h_mean /= len(test_loader)

    print(f"[final performance]", Metrics(p_opt, p_exp, h_mean))
    np.savetxt(f"{logdir}/score.txt", [p_opt, p_exp, h_mean])

if __name__ == "__main__":
    main()


text:
#%% variables are initialized in const.py file
"""
    AGENT = how many samples to collect in one step= 1000
    EPSILON_START = 0.50
    CLASSES = 2 = as we are loading miniboones dataset(mb) by default, it has 2 classes(0 and 1)
    FEATURE_DIM = 50 = as we are loading miniboones dataset(mb) by default, it has 50 features
    ACTION_DIM = FEATURE_DIM + CLASSES
    STATE_DIM = FEATURE_DIM * 2( multiply by 2 because of the mask)
              = First two coloum of the SATE_DIM specify classes(assuming miniboone dataset is loaded by deffault)
    MAX_MASK_CONST = int(1e6) = for not considering the actions that are already performed
"""
#%%
import numpy as np
from env import Environment
from IPython.core.debugger import set_trace
from consts import *
from ipdb import set_trace

class Agent():
    def __init__(self, env, pool, brain):
        self.env  = env
        self.pool = pool
        self.brain = brain

        self.epsilon = EPSILON_START
        self.s = self.env.reset()  # initialize state to zeros.

    def store(self, x):
        """
        Params:
            x = sample expreiences = (s, a, s', r)
        """
        self.pool.put(x)

    def act(self, s):
        """
        Params:
            s = batch of input states(dimension = FEATURE_DIM *2)
        Outputs:
            a = batch of actions
        """
        m = np.zeros((AGENTS, ACTION_DIM))    # create max_mask
        m[:, CLASSES:] = s[:, FEATURE_DIM:]

        if self.epsilon < 1.0:
            p = self.brain.predict_np(s) - MAX_MASK_CONST * m     # select an action not considering those  already performed
            a = np.argmax(p, axis=1)
        else:
            a = np.zeros(AGENTS, dtype=np.int32)

        # override with random action
        rand_agents = np.where( np.random.rand(AGENTS) < self.epsilon )[0] #returns a touple , [0] = first index of the touple
        rand_number = np.random.rand(len(rand_agents))

        for i in range(len(rand_agents)):
            agent = rand_agents[i]

            possible_actions = np.where( m[agent] == 0. )[0]     # select a random action, don't repeat an action
            w = int(rand_number[i] * len(possible_actions))
            a[agent] = possible_actions[w]

        return a

    def step(self):
        """
        store (states, actions, rewards, next states)
        """
        a = self.act(self.s)
        s_, r = self.env.step(a)

        self.store( (self.s, a, r, s_) )

        self.s = s_

    def update_epsilon(self, epoch):
        if epoch >= EPSILON_EPOCHS:
            self.epsilon = EPSILON_END
        else:
            self.epsilon = EPSILON_START + epoch * (EPSILON_END - EPSILON_START) / EPSILON_EPOCHS


text:
import sys
import os
sys.path.insert(1, os.path.join(sys.path[0], '..'))

import numpy as np

from base.least_square_estimator import LeastSquareEstimator

class MagnetometerCalibrator:
    """
    Online calibration of the hard-iron and the soft-iron effects.
    FIXME: using a sphere model now, which produces the hard-iron calibration only.
           we may need to upgrade to an ellipsoid model to include the soft-iron
           calibration.
    """
    MIN_SAMPLE_COUNT = 24
    COARSE_FACTOR = 4
    INITIAL_COVAR = 1e3
    ERROR_MARGIN_RATIO = 0.1

    def __init__(self, center_guess):
        """
        :param center_guess: 3-by-1 numpy array; initial guess for the static bias
        :param radius_guess: initial guess for the magnetic field intensity
        """
        self._samples = {}
        self._coarse_samples = {}
        initial_guess = np.array([center_guess[0],
                                  center_guess[1],
                                  center_guess[2],
                                  0.0], dtype=np.float)
        initial_covar = np.eye(4) * MagnetometerCalibrator.INITIAL_COVAR
        self._estimator = LeastSquareEstimator(initial_guess, initial_covar)

    def update(self, mag):
        """
        Update the mag bias estimate with the new raw measurement.
        :param mag: 3-by-1 numpy array; raw measurement from the magnetometer.
        """
        self._add_sample(mag)
        (phi, y) = self._convert_to_estimator_sample(mag)
        self._estimator.update(phi, y)

        #estimate = self._estimator.get_estimate_mean()
        #print("[MagCal] center={}, sample size={}, coarse sample size={}".format(\
        #        estimate[0:3], len(self._samples), len(self._coarse_samples)))

    @property
    def bias(self):
        """
        :return static bias of the magnetometer; or None if not yet available.
        """
        if (not self._ready()):
            return None
        else:
            estimate = self._estimator.get_estimate_mean()
            center = np.array(estimate[0:3], dtype=np.float)
            return center

    @property
    def intensity(self):
        """
        :return intensity of the external magnetic field; or None if not yet available.
        """
        if (not self._ready()):
            return None
        else:
            estimate = self._estimator.get_estimate_mean()
            center = np.array(estimate[0:3], dtype=np.float)
            radius = np.sqrt(estimate[3] + np.dot(center, center))
            return radius

    def calibrate_measurement(self, mag, error_margin_ratio=None):
        """
        :return the calibrated measurement, or None if not yet calibrated, or the raw
                measurement does not fit into the model within the specified error margin.
        """
        if (not self._ready()):
            return None

        if (error_margin_ratio is None):
            error_margin_ratio = MagnetometerCalibrator.ERROR_MARGIN_RATIO

        center = self.bias
        radius = self.intensity
        calibrated = mag - center

        # measurement is good if it approximately lies on the estimated sphere
        error = np.abs(radius - np.linalg.norm(calibrated))
        if (error < error_margin_ratio * radius):
            return calibrated
        else:
            # possibly a corrupted measurement
            return None

    def _convert_to_estimator_sample(self, mag):
        """
        convert raw magnetometer measurements to estimator samples
        :return (phi, y) the sample input and the sample output.
        """
        # a sphere centered at [x_c, y_c, z_c] with radius R is described by
        # (x - x_c)**2 + (y - y_c)**2 + (z - z_c)**2 == R**2
        # rearranging, we have
        # [2*x, 2*y, 2*z, 1] * [x_c, y_c, z_c, d].T == x**2 + y**2 + z**2
        # where d = R**2 - x_c**2 - y_c**2 - z_c**2

        phi = np.array([mag[0] * 2.0,
                        mag[1] * 2.0,
                        mag[2] * 2.0,
                        1.0], dtype=np.float)
        y = np.dot(mag, mag)
        return (phi, y)

    def _ready(self):
        k = MagnetometerCalibrator.MIN_SAMPLE_COUNT
        return ((len(self._samples) > k) and \
                (len(self._coarse_samples) > k))

    def _add_sample(self, mag):
        """
        Add the provided sample to sample history
        """
        hash_value = self._get_sample_hash(mag)
        self._samples[hash_value] = mag

        hash_value = self._get_coarse_sample_hash(mag)
        self._coarse_samples[hash_value] = mag

    def _get_sample_hash(self, mag):
        """
        Calculate the hash value for the given mag value.
        :param mag: 3-by-1 numpy array.
        """
        x = np.int(np.round(mag[0]))
        y = np.int(np.round(mag[1]))
        z = np.int(np.round(mag[2]))
        return x ^ y ^ z

    def _get_coarse_sample_hash(self, mag):
        """
        similar to _get_sample_hash, but uses a larger bin
        """
        k = 1.0 / MagnetometerCalibrator.COARSE_FACTOR
        coarse_mag = [mag[0] * k, mag[1] * k, mag[2] * k]
        return self._get_sample_hash(coarse_mag)

if (__name__ == "__main__"):
    pass

max_stars_repo_name: lonelycorn/AHRS

max_issues_repo_name: lonelycorn/AHRS

max_forks_repo_name: lonelycorn/AHRS


text:
import pytest
import numpy as np

from kdtrees import _utils as utils
from .test_fixtures import KDSubType, BadType

@pytest.mark.parametrize("args,dim_exp", [
	([[0],[1]], 1),
	([[0]], 1),
	([[0],[1],[2]], 1),
	([[0]], 1),
	([[0,1]], 2),
	([[[0],[1]]], 1),
	([[[0,1]]], 2),
	([[[0,1]]], 2),
	([[[0,1],[2,3]]], 2),
	([np.asarray([0,1])], 2),
])

def test_check_dimensionality(args, dim_exp):
	if len(args) > 1:
		assert utils.check_dimensionality(*args) == dim_exp
	else:
		assert utils.check_dimensionality(args[0]) == dim_exp

def test_check_dimensionality_accept():
	assert utils.check_dimensionality(KDSubType(1,1), accept=KDSubType) == 1

def test_check_dimensionality_no_args():
	with pytest.raises(ValueError):
		utils.check_dimensionality()

def test_check_dimensionality_mismatch():
	with pytest.raises(ValueError):
		utils.check_dimensionality([0], [0,1])

def test_check_dimensionality_no_dim():
	with pytest.raises(AttributeError):
		utils.check_dimensionality(BadType(), accept=BadType)

def test_check_dimensionality_no_shape():
	with pytest.raises(AttributeError):
		utils.check_dimensionality([0], accept=BadType)

def test_check_dimensionality_mismatch_accept():
	with pytest.raises(ValueError):
		utils.check_dimensionality(KDSubType(1,1), KDSubType(2,1), accept=KDSubType)

max_stars_repo_name: paradoxysm/kdtrees

max_issues_repo_name: paradoxysm/kdtrees

max_forks_repo_name: paradoxysm/kdtrees


text:
#!/usr/bin/env python2
import sys
import rospy
import numpy as np

from geometry_msgs.msg import Twist, PoseStamped
from nav_msgs.msg import Odometry
from gazebo_msgs.srv import GetModelState
from tf.transformations import euler_from_quaternion, quaternion_from_euler

from RVO import RVO_update, reach, compute_V_des, reach, distance

from potential_field import PotentialField

def twist_to_xyyaw(data):
        position = data.pose.position
        x = position.x
        y = position.y
        orientation = data.pose.orientation
        quaternion = [orientation.x,
                      orientation.y,
                      orientation.z,
                      orientation.w]
        euler = euler_from_quaternion(quaternion)
        yaw = euler[2]

        return x, y, yaw

def get_robot_pose(robot_name):
    ''' Ground truth
    '''
    rospy.wait_for_service('/gazebo/get_model_state')
    try:
        robot_state_service = rospy.ServiceProxy('/gazebo/get_model_state', GetModelState)
        robot_state_res = robot_state_service(robot_name, "")
    except rospy.ServiceException as e:
        print("Service call failed: %s"%e)

    return twist_to_xyyaw(robot_state_res)

def map_xy_to_xz(v_robot, t_diff, k):
    if k * t_diff == 0:
        return 0, 0
    # Ori.z
    desired_theta = np.arctan2(v_robot[1], v_robot[0])
    omega = desired_theta / (k * (t_diff))
    # Vel.x
    vel = np.sqrt((v_robot[0])**2 + (v_robot[1])**2)

    return vel, omega

# Another way to encounter while loop:
# def float_range(start, stop, step):
#     while start < stop:
#         yield float(start)
#         start += decimal.Decimal(step)

def rvo_test():
    # ROS node
    rospy.init_node('RVO_Bot')
    rate = rospy.Rate(10)
    velocity_publisher = rospy.Publisher('gopher_1/base_controller/cmd_vel', Twist, queue_size=1)
    visual_publisher = rospy.Publisher('gopher_1/rvo', PoseStamped, queue_size=1)
    vel_msg = Twist()

    # RVO setting
    # define workspace model
    ws_model = dict()
    # robot radius
    ws_model['robot_radius'] = 0.2
    # circular obstacles, format [x, y, rad]
    obstacles = [[-3, 5, 1]]
    ws_model['circular_obstacles'] = obstacles

    # Robot settings
    # initialization
    X = [[0, 0]]
    V = [[0, 0]]
    # robot pos & goal
    goal = [[-1.5, 7.5]]
    # robot initial vel & max velocity
    V_max = [1, 1]

    # Simulation starts
    last_time = rospy.Time.now().to_sec()

    pf = PotentialField()
    pf.set_goal(goal[0])
    pf.add_obstacles(obstacles)

    while not rospy.is_shutdown():
        rospy.loginfo("Start simulation")

        # Another way to encounter while loop:
        # while round(distance(X[0],goal[0]),2) not in list(float_range(0, 0.2, '0.01')):

        while not (distance(X[0], goal[0]) < 0.5):
            # Get the latest position of the robot
            x, y, yaw = get_robot_pose("gopher_1")
            X = [[x, y]]

            # Run RVO
            # compute desired vel to goal
            V_des = compute_V_des(X, goal, V_max)

            # compute force
            force = pf.compute_force(X[0])
            # print("force: %.2f"%(force))

            # next step
            descent_rate = 1.0

            V_potential = descent_rate * force
            # delta_x = next_move[0]
            # delta_y = next_move[1]

            # compute the optimal vel to avoid collision
            # V = RVO_update(X, V_des, V, ws_model)

            V = RVO_update(X, V_potential, V, ws_model)

            # Map vx, vy to v_pos.x, v_ori.z
            # delta_t
            current_time = rospy.Time.now().to_sec()
            time_diff = current_time - last_time
            last_time = current_time
            # w.r.t. robot

            # We are rotating the local coordinates of robot w.r.t global coordinates (Rotation matrix * [x,y])
            v_robot = [ np.cos(yaw)*V[0][0] + np.sin(yaw)*V[0][1],
                       -np.sin(yaw)*V[0][0] + np.cos(yaw)*V[0][1]]
            # to v.x and o.z
            vel, omega = map_xy_to_xz(v_robot, time_diff, 3)

            # Publish to robot
            vel_msg.linear.x = float(vel)
            vel_msg.angular.z = float(omega)
            velocity_publisher.publish(vel_msg)

            # Visualize RVO result
            vis = PoseStamped()
            vis.header.stamp = rospy.Time.now()
            vis.header.frame_id = "map"
            vis.pose.position.x = x
            vis.pose.position.y = y
            [vis.pose.orientation.x, vis.pose.orientation.y, vis.pose.orientation.z,
            vis.pose.orientation.w] = quaternion_from_euler(0, 0, np.arctan2(V[0][1], V[0][0]))
            visual_publisher.publish(vis)

            print("Pos_X: %.2f, Pos_Y: %.2f" %(X[0][0], X[0][1]))
            print("Vx: %.02f, Vy: %.2f"%(v_robot[0], v_robot[1]))
            print("time_diff %.4f"%(time_diff))
            print("Vel: %.2f, Omega: %.2f"%(vel, omega))
            # print

            rate.sleep()

        # Stop robot
        velocity_publisher.publish(Twist())
        break

if __name__ == '__main__':
    try:
        rvo_test()
    except rospy.ROSInterruptException:
        rospy.loginfo('Shutting down')


text:
"""Print a data summary text file."""
import logging
from pathlib import Path
from gammapy.data import DataStore
import astropy.units as u
from .conf import config

log = logging.getLogger(__name__)

def make_summary_data():
    lines = make_summary_data_lines()

    filename = "results/summary/data.md"
    log.info(f"Writing {filename}")
    Path(filename).write_text("\n".join(lines))

def make_summary_data_lines():
    yield "# Data summary\n"
    yield ""

    for name in config.all_datasets:
        for line in make_summary_for_dataset(name):
            yield line

def make_summary_for_dataset(name):
    dataset = config.datasets[name]
    spec = dataset.get_SpectrumObservationList().stack()
    stats = spec.total_stats

    yield f"## {dataset.name}"
    yield "\nInfos from dataset:\n"

    energy_min = dataset.energy_range[0].to("TeV").value
    yield f"- energy_min: {energy_min:.5f} TeV"

    energy_max = dataset.energy_range[1].to("TeV").value
    yield f"- energy_max: {energy_max:.5f} TeV"

    r_on = dataset.on_radius.to("deg").value
    yield f"- r_on: {r_on:.5f} deg"

    yield "\nInfos from spec:\n"
    livetime = spec.livetime.to("hour").value
    yield f"- livetime: {livetime:.2f} hours"

    # n_on = spec.on_vector.data.data.value.sum()
    yield f"- n_on: {stats.n_on}"
    yield f"- excess: {stats.excess:.1f}"
    yield f"- background: {stats.background:.1f}"

    if dataset.name != "fermi":
        data_store = DataStore.from_dir(f"data/{dataset.name}")
        yield "\nInfos from datastore:\n"

        n_obs = len(dataset.obs_ids)
        yield f"- n_obs: {n_obs}"

        obstime = data_store.obs_table["ONTIME"].sum() / 3600
        yield f"- obstime: {obstime:.2f} hours"

    yield ""

def make_summary_latex():
    """Make summary LaTeX table for the paper."""
    path = "results/summary/data.tex"
    log.info(f"Writing {path}")

    table = open(path, "w")
    table.write(r"\begin{tabular}{lrrrrrl}")
    table.write("\n")
    table.write(r"\hline")
    table.write("\n")
    table.write(
        r"Dataset &  $T_{\rm obs}$ & $E_{\rm min}$ & $E_{\rm max}$ & $N_{\rm on}$ & $N_{\rm bkg}$ & $R_{\rm on}$  \\"
    )
    table.write("\n")
    table.write(
        r"        &       & TeV           & TeV           &              &         & deg        \\ \hline"
    )
    table.write("\n")

    for name in config.all_datasets:
        dataset = config.datasets[name]
        e_min_list = []
        for _spec in dataset.get_SpectrumObservationList():
            e_min_list.append(_spec.lo_threshold)
        spec = dataset.get_SpectrumObservationList().stack()
        stats = spec.total_stats

        row_name = rf"\{name} & "
        if name == "fermi":
            T_obs = r"$\sim$7 yr & "
            e_min = dataset.energy_range[0].to("TeV").value
            E_min = f"{e_min:.2f} & "
        else:
            data_store = DataStore.from_dir(f"data/{dataset.name}")
            ontime = sum(data_store.obs_table["ONTIME"]) * u.s
            ontime = ontime.to("h").value
            T_obs = f"{ontime:.2f} h & "
            # in case of the IACT e_min is taken from the staked obs
            e_min = min(e_min_list)
            E_min = f"{e_min.to('TeV').value:.2f} & "

        e_max = dataset.energy_range[1].to("TeV").value
        E_max = f"{e_max:.0f} & "
        N_on = f"{stats.n_on} & "
        N_bkg = f"{stats.background:.1f} & "
        r_on = dataset.on_radius.to("deg").value
        R_on = rf"{r_on:.2f} \\"

        table.write(row_name + T_obs + E_min + E_max + N_on + N_bkg + R_on)
        table.write("\n")

    table.write("\hline")
    table.write("\n")
    table.write("\end{tabular}")
    table.close()

max_stars_repo_name: open-gamma-ray-astro/joint-crab

max_issues_repo_name: open-gamma-ray-astro/joint-crab

max_forks_repo_name: open-gamma-ray-astro/joint-crab

reason: import astropy


text:
#!/usr/bin/env python

import pickle
import argparse

import plotting

import numpy as np

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("--color-bar", action='store_true')
    parser.add_argument("activations")
    args = parser.parse_args()

    with open(args.activations, 'rb') as f:
        activations = pickle.load(f)

    fig_titles = []
    for a in activations:
        nb_zeros = np.sum(a==0.0)
        fig_titles.append(", ".join([str(x) for x in [nb_zeros, np.mean(a), np.max(a)]]))

    color_setup = {
        'cmap': 'gnuplot',
        'vmin': 0,
        'vmax': 30,
        'colorbar': args.color_bar,
    }

    plotting.grid_plot(activations, lambda x:x, "activations", fig_titles, coloring=color_setup)

max_stars_repo_name: entn-at/BrnoLM

max_issues_repo_name: entn-at/BrnoLM

max_forks_repo_name: entn-at/BrnoLM


text:
#import use_ngmca as use
import pyfits as pf
import numpy as np
import matplotlib.pyplot as plt
from MuSCADeT import MCA
from MuSCADeT import wave_transform as mw
import scipy.stats as sc
from MuSCADeT import colour_subtraction as cs

## Initialise your data with a cube with size (nb,n1,n2) where nb is the number of bands in your observations
cube = ##To be filed

hdus = pf.PrimaryHDU(cube)
lists = pf.HDUList([hdus])
lists.writeto('Simu_For_User/Cube.fits', clobber=True)

num,n,n = np.shape(cube)

## If A is unknown, set it to zero.
Aprior = 0 #pf.open('Simu_simple/Simu_A.fits')[0].data

## Input parameters
pca = 'PCA'     #Estimation of the mixing coefficients from PCA. If different from PCA it will use the array provided in Aprior
n = 100        #Number of iterations: increase if the separation is not good enough!
nsig = 5        #Threshold in units of noise standard deviation: Can be lowered down to 3 but 5 should be fine.
ns = 2          #Number of sources:
angle = 10      #Resolution angle for the PCA colour estimation (start with 15 then adjust empirically)
alpha = [0,0]   #If automated estimation of PCA coefficients fails, chose adequate alphas. See readme for more details
plot = False     #option to plot the PCA coefficients of the SEDs in the image. This option is usefull if one wants to make sure that SEDs have been correctly estimated. In automated mode, keep this option at False. In case the SEDs have to be refined, set plot to True, identify the features (alignements) on the plot that stand for different SEDs and use this to give values for alpha. (see readme.)

## Running MuSCADeT
S,A = MCA.mMCA(cube, Aprior, nsig,n, PCA=[ns,angle], mode=pca, alpha = [0,0])

## MuSCADeT estimates ns source, which means, variable S contains ns images
## Saves the sources in a fits file
hdus = pf.PrimaryHDU(S)
lists = pf.HDUList([hdus])
lists.writeto('Simu_For_User/Sources_'+str(n)+'.fits', clobber=True)

## Saves the estimated mixing matrix, A, in a fits file
hdus = pf.PrimaryHDU(A)
lists = pf.HDUList([hdus])
lists.writeto('Simu_For_User/Estimated_A.fits', clobber=True)

## This command shows the result of the separation in various formats:
## The command needs : A fits file with the extracted source and mixing coefficients (MuSCADeT's outputs)
##                     A fits file with the original data
##                     An optional prefix for file names
##                     Cuts: optional, these are the lower and upper cuts for Red, Green, Blue bands and for the sources to be displayed
## This command will show in DS9 the sources as extrated by MuSCADeT, the original images in RGB, the subtraction of blue (or red) component from the original data, the subtraction of red (or bleu) component from the original data and finally the residuals. If the separation is successful, the residuals should present only nonise.
cs.make_colour_sub('Simu_For_User/Sources_'+str(n)+'.fits',
                   'Simu_For_User/Estimated_A.fits',
                   './Simu_For_User/Cube.fits',
                   prefix = './Simu_For_User/',
                   cuts = ['-0.1','0.6','-0.05','0.3','-0.02','0.1', '0','0.5'])

max_stars_repo_name: aymgal/MuSCADeT

max_issues_repo_name: aymgal/MuSCADeT

max_forks_repo_name: aymgal/MuSCADeT

reason: import numpy,import scipy