text:
import matplotlib as mpl
import numpy as np
from matplotlib.backends.backend_pgf import (
    common_texification as mpl_common_texification,
)

from . import _color

def _common_texification(string):
    # Work around <https://github.com/matplotlib/matplotlib/issues/15493>
    return mpl_common_texification(string).replace("&", "\\&")

class Axes:
    def __init__(self, data, obj):  # noqa: C901
        """Returns the PGFPlots code for an axis environment."""
        self.content = []

        # Are we dealing with an axis that hosts a colorbar? Skip then, those are
        # treated implicitily by the associated axis.
        self.is_colorbar = _is_colorbar_heuristic(obj)
        if self.is_colorbar:
            return

        # instantiation
        self.nsubplots = 1
        self.subplot_index = 0
        self.is_subplot = False

        if isinstance(obj, mpl.axes.Subplot):
            self._subplot(obj, data)

        self.axis_options = []

        # check if axes need to be displayed at all
        if not obj.axison:
            self.axis_options.append("hide x axis")
            self.axis_options.append("hide y axis")

        # get plot title
        title = obj.get_title()
        data["current axis title"] = title
        if title:
            title = _common_texification(title)
            self.axis_options.append(f"title={{{title}}}")

        # get axes titles
        xlabel = obj.get_xlabel()
        if xlabel:
            xlabel = _common_texification(xlabel)
            self.axis_options.append(f"xlabel={{{xlabel}}}")
            xrotation = obj.xaxis.get_label().get_rotation()
            if xrotation != 0:
                self.axis_options.append(f"xlabel style={{rotate={xrotation - 90}}}")
        ylabel = obj.get_ylabel()
        if ylabel:
            ylabel = _common_texification(ylabel)
            self.axis_options.append(f"ylabel={{{ylabel}}}")
            yrotation = obj.yaxis.get_label().get_rotation()
            if yrotation != 90:
                self.axis_options.append(f"ylabel style={{rotate={yrotation - 90}}}")

        # Axes limits.
        ff = data["float format"]
        xlim = list(obj.get_xlim())
        xlim0, xlim1 = sorted(xlim)
        ylim = list(obj.get_ylim())
        ylim0, ylim1 = sorted(ylim)
        # Sort the limits so make sure that the smaller of the two is actually *min.
        self.axis_options.append(f"xmin={xlim0:{ff}}, xmax={xlim1:{ff}}")
        self.axis_options.append(f"ymin={ylim0:{ff}}, ymax={ylim1:{ff}}")
        # When the axis is inverted add additional option
        if xlim != sorted(xlim):
            self.axis_options.append("x dir=reverse")
        if ylim != sorted(ylim):
            self.axis_options.append("y dir=reverse")

        # axes scaling
        if obj.get_xscale() == "log":
            self.axis_options.append("xmode=log")
            self.axis_options.append(
                f"log basis x={{{_try_f2i(obj.xaxis._scale.base)}}}"
            )
        if obj.get_yscale() == "log":
            self.axis_options.append("ymode=log")
            self.axis_options.append(
                f"log basis y={{{_try_f2i(obj.yaxis._scale.base)}}}"
            )

        # Possible values for get_axisbelow():
        #   True (zorder = 0.5):   Ticks and gridlines are below all Artists.
        #   'line' (zorder = 1.5): Ticks and gridlines are above patches (e.g.
        #                          rectangles) but still below lines / markers.
        #   False (zorder = 2.5):  Ticks and gridlines are above patches and lines /
        #                          markers.
        if not obj.get_axisbelow():
            self.axis_options.append("axis on top")

        # aspect ratio, plot width/height
        aspect = obj.get_aspect()
        if aspect in ["auto", "normal"]:
            aspect_num = None  # just take the given width/height values
        elif aspect == "equal":
            aspect_num = 1.0
        else:
            aspect_num = float(aspect)

        self._set_axis_dimensions(data, aspect_num, xlim, ylim)

        # axis positions
        xaxis_pos = obj.get_xaxis().label_position
        if xaxis_pos == "top":
            # default: "bottom"
            self.axis_options.append("axis x line=top")

        yaxis_pos = obj.get_yaxis().label_position
        if yaxis_pos == "right":
            # default: "left"
            self.axis_options.append("axis y line=right")

        self._ticks(data, obj)

        self._grid(obj, data)

        # axis line styles
        # Assume that the bottom edge color is the color of the entire box.
        axcol = obj.spines["bottom"].get_edgecolor()
        data, col, _ = _color.mpl_color2xcolor(data, axcol)
        if col != "black":
            self.axis_options.append(f"axis line style={{{col}}}")

        # background color
        bgcolor = obj.get_facecolor()

        data, col, _ = _color.mpl_color2xcolor(data, bgcolor)
        if col != "white":
            self.axis_options.append(f"axis background/.style={{fill={col}}}")

        # find color bar
        colorbar = _find_associated_colorbar(obj)
        if colorbar:
            self._colorbar(colorbar, data)

        if self.is_subplot:
            self.content.append("\n\\nextgroupplot")
        else:
            self.content.append(data["flavor"].start("axis"))

    def get_begin_code(self):
        content = self.content
        if self.axis_options:
            # Put axis_options in a deterministic order to avoid diff churn.
            self.axis_options.sort()
            content.append("[\n" + ",\n".join(self.axis_options) + "\n]\n")
        return content

    def get_end_code(self, data):
        if not self.is_subplot:
            return data["flavor"].end("axis") + "\n\n"
        elif self.is_subplot and self.nsubplots == self.subplot_index:
            data["is_in_groupplot_env"] = False
            return data["flavor"].end("groupplot") + "\n\n"

        return ""

    def _set_axis_dimensions(self, data, aspect_num, xlim, ylim):
        if data["axis width"] and data["axis height"]:
            # width and height overwrite aspect ratio
            self.axis_options.append("width=" + data["axis width"])
            self.axis_options.append("height=" + data["axis height"])
        elif data["axis width"]:
            # only data["axis width"] given. calculate height by the aspect ratio
            self.axis_options.append("width=" + data["axis width"])
            if aspect_num:
                alpha = aspect_num * (ylim[1] - ylim[0]) / (xlim[1] - xlim[0])
                if alpha == 1.0:
                    data["axis height"] = data["axis width"]
                else:
                    # Concatenate the literals, as data["axis width"] could as well
                    # be a LaTeX length variable such as \figurewidth.
                    data["axis height"] = str(alpha) + "*" + data["axis width"]
                self.axis_options.append("height=" + data["axis height"])
        elif data["axis height"]:
            # only data["axis height"] given. calculate width by the aspect ratio
            self.axis_options.append("height=" + data["axis height"])
            if aspect_num:
                alpha = aspect_num * (ylim[1] - ylim[0]) / (xlim[1] - xlim[0])
                if alpha == 1.0:
                    data["axis width"] = data["axis height"]
                else:
                    # Concatenate the literals, as data["axis height"] could as
                    # well be a LaTeX length variable such as \figureheight.
                    data["axis width"] = str(1.0 / alpha) + "*" + data["axis height"]
                self.axis_options.append("width=" + data["axis width"])
        else:
            # TODO keep an eye on https://tex.stackexchange.com/q/480058/13262
            pass

    def _ticks(self, data, obj):
        # get ticks
        self.axis_options.extend(
            _get_ticks(data, "x", obj.get_xticks(), obj.get_xticklabels())
        )
        self.axis_options.extend(
            _get_ticks(data, "y", obj.get_yticks(), obj.get_yticklabels())
        )
        self.axis_options.extend(
            _get_ticks(
                data,
                "minor x",
                obj.get_xticks(minor=True),
                obj.get_xticklabels(minor=True),
            )
        )
        self.axis_options.extend(
            _get_ticks(
                data,
                "minor y",
                obj.get_yticks(minor=True),
                obj.get_yticklabels(minor=True),
            )
        )

        try:
            l0 = obj.get_xticklines()[0]
        except IndexError:
            pass
        else:
            c0 = l0.get_color()
            data, xtickcolor, _ = _color.mpl_color2xcolor(data, c0)
            self.axis_options.append(f"xtick style={{color={xtickcolor}}}")

        try:
            l0 = obj.get_yticklines()[0]
        except IndexError:
            pass
        else:
            c0 = l0.get_color()
            data, ytickcolor, _ = _color.mpl_color2xcolor(data, c0)
            self.axis_options.append(f"ytick style={{color={ytickcolor}}}")

        # Find tick direction
        # For new matplotlib versions, we could replace the direction getter by
        # `get_ticks_direction()`, see
        # <https://github.com/matplotlib/matplotlib/pull/5290>.  Unfortunately, _tickdir
        # doesn't seem to be quite accurate. See
        # <https://github.com/matplotlib/matplotlib/issues/5311>.  For now, just take
        # the first tick direction of each of the axes.
        x_tick_dirs = [tick._tickdir for tick in obj.xaxis.get_major_ticks()]
        y_tick_dirs = [tick._tickdir for tick in obj.yaxis.get_major_ticks()]
        if x_tick_dirs or y_tick_dirs:
            if x_tick_dirs and y_tick_dirs:
                direction = x_tick_dirs[0] if x_tick_dirs[0] == y_tick_dirs[0] else None
            elif x_tick_dirs:
                direction = x_tick_dirs[0]
            else:
                # y_tick_dirs must be present
                direction = y_tick_dirs[0]

            if direction:
                if direction == "in":
                    # 'tick align=inside' is the PGFPlots default
                    pass
                elif direction == "out":
                    self.axis_options.append("tick align=outside")
                else:
                    assert direction == "inout"
                    self.axis_options.append("tick align=center")

        # Set each rotation for every label
        x_tick_rotation_and_horizontal_alignment = (
            self._get_label_rotation_and_horizontal_alignment(obj, data, "x")
        )
        if x_tick_rotation_and_horizontal_alignment:
            self.axis_options.append(x_tick_rotation_and_horizontal_alignment)

        y_tick_rotation_and_horizontal_alignment = (
            self._get_label_rotation_and_horizontal_alignment(obj, data, "y")
        )
        if y_tick_rotation_and_horizontal_alignment:
            self.axis_options.append(y_tick_rotation_and_horizontal_alignment)

        # Set tick position
        x_tick_position_string, x_tick_position = _get_tick_position(obj, "x")
        y_tick_position_string, y_tick_position = _get_tick_position(obj, "y")

        if x_tick_position == y_tick_position and x_tick_position is not None:
            self.axis_options.append(f"tick pos={x_tick_position}")
        else:
            self.axis_options.append(x_tick_position_string)
            self.axis_options.append(y_tick_position_string)

    def _grid(self, obj, data):
        # Don't use get_{x,y}gridlines for gridlines; see discussion on
        # <http://sourceforge.net/p/matplotlib/mailman/message/25169234/> Coordinate of
        # the lines are entirely meaningless, but styles (colors,...) are respected.

        try:
            # mpl 3.3.3+
            # <https://github.com/matplotlib/matplotlib/pull/18769>
            has_major_xgrid = obj.xaxis._major_tick_kw["gridOn"]
            has_minor_xgrid = obj.xaxis._minor_tick_kw["gridOn"]
            has_major_ygrid = obj.yaxis._major_tick_kw["gridOn"]
            has_minor_ygrid = obj.yaxis._minor_tick_kw["gridOn"]
        except KeyError:
            has_major_xgrid = obj.xaxis._gridOnMajor
            has_minor_xgrid = obj.xaxis._gridOnMinor
            has_major_ygrid = obj.yaxis._gridOnMajor
            has_minor_ygrid = obj.yaxis._gridOnMinor

        if has_major_xgrid:
            self.axis_options.append("xmajorgrids")
        if has_minor_xgrid:
            self.axis_options.append("xminorgrids")

        xlines = obj.get_xgridlines()
        if xlines:
            xgridcolor = xlines[0].get_color()
            data, col, _ = _color.mpl_color2xcolor(data, xgridcolor)
            if col != "black":
                self.axis_options.append(f"x grid style={{{col}}}")

        if has_major_ygrid:
            self.axis_options.append("ymajorgrids")
        if has_minor_ygrid:
            self.axis_options.append("yminorgrids")

        ylines = obj.get_ygridlines()
        if ylines:
            ygridcolor = ylines[0].get_color()
            data, col, _ = _color.mpl_color2xcolor(data, ygridcolor)
            if col != "black":
                self.axis_options.append(f"y grid style={{{col}}}")

    def _colorbar(self, colorbar, data):
        colorbar_styles = []

        orientation = colorbar.orientation
        limits = colorbar.mappable.get_clim()
        if orientation == "horizontal":
            self.axis_options.append("colorbar horizontal")

            colorbar_ticks = colorbar.ax.get_xticks()
            colorbar_ticks_minor = colorbar.ax.get_xticks(minor=True)
            axis_limits = colorbar.ax.get_xlim()

            # In matplotlib, the colorbar color limits are determined by get_clim(), and
            # the tick positions are as usual with respect to {x,y}lim. In PGFPlots,
            # however, they are mixed together.  Hence, scale the tick positions just
            # like {x,y}lim are scaled to clim.
            colorbar_ticks = (colorbar_ticks - axis_limits[0]) / (
                axis_limits[1] - axis_limits[0]
            ) * (limits[1] - limits[0]) + limits[0]
            colorbar_ticks_minor = (colorbar_ticks_minor - axis_limits[0]) / (
                axis_limits[1] - axis_limits[0]
            ) * (limits[1] - limits[0]) + limits[0]
            # Getting the labels via get_* might not actually be suitable:
            # they might not reflect the current state.
            colorbar_ticklabels = colorbar.ax.get_xticklabels()
            colorbar_ticklabels_minor = colorbar.ax.get_xticklabels(minor=True)

            colorbar_styles.extend(
                _get_ticks(data, "x", colorbar_ticks, colorbar_ticklabels)
            )
            colorbar_styles.extend(
                _get_ticks(
                    data, "minor x", colorbar_ticks_minor, colorbar_ticklabels_minor
                )
            )

        else:
            assert orientation == "vertical"

            self.axis_options.append("colorbar")
            colorbar_ticks = colorbar.ax.get_yticks()
            colorbar_ticks_minor = colorbar.ax.get_yticks(minor=True)
            axis_limits = colorbar.ax.get_ylim()

            # In matplotlib, the colorbar color limits are determined by get_clim(), and
            # the tick positions are as usual with respect to {x,y}lim. In PGFPlots,
            # however, they are mixed together.  Hence, scale the tick positions just
            # like {x,y}lim are scaled to clim.
            colorbar_ticks = (colorbar_ticks - axis_limits[0]) / (
                axis_limits[1] - axis_limits[0]
            ) * (limits[1] - limits[0]) + limits[0]
            colorbar_ticks_minor = (colorbar_ticks_minor - axis_limits[0]) / (
                axis_limits[1] - axis_limits[0]
            ) * (limits[1] - limits[0]) + limits[0]

            # Getting the labels via get_* might not actually be suitable:
            # they might not reflect the current state.
            colorbar_ticklabels = colorbar.ax.get_yticklabels()
            colorbar_ylabel = colorbar.ax.get_ylabel()
            colorbar_ticklabels_minor = colorbar.ax.get_yticklabels(minor=True)
            colorbar_styles.extend(
                _get_ticks(data, "y", colorbar_ticks, colorbar_ticklabels)
            )
            colorbar_styles.extend(
                _get_ticks(
                    data, "minor y", colorbar_ticks_minor, colorbar_ticklabels_minor
                )
            )
            colorbar_styles.append("ylabel={" + colorbar_ylabel + "}")

        mycolormap, is_custom_cmap = _mpl_cmap2pgf_cmap(
            colorbar.mappable.get_cmap(), data
        )
        if is_custom_cmap:
            self.axis_options.append("colormap=" + mycolormap)
        else:
            self.axis_options.append("colormap/" + mycolormap)

        ff = data["float format"]
        self.axis_options.append(f"point meta min={limits[0]:{ff}}")
        self.axis_options.append(f"point meta max={limits[1]:{ff}}")

        if colorbar_styles:
            self.axis_options.append(
                "colorbar style={{{}}}".format(",".join(colorbar_styles))
            )

    def _subplot(self, obj, data):
        # https://github.com/matplotlib/matplotlib/issues/7225#issuecomment-252173667
        geom = obj.get_subplotspec().get_topmost_subplotspec().get_geometry()

        self.nsubplots = geom[0] * geom[1]
        if self.nsubplots > 1:
            # Is this an axis-colorbar pair? No need for groupplot then.
            is_groupplot = self.nsubplots != 2 or not _find_associated_colorbar(obj)

            if is_groupplot:
                self.is_subplot = True
                # subplotspec geometry positioning is 0-based
                self.subplot_index = geom[2] + 1
                if "is_in_groupplot_env" not in data or not data["is_in_groupplot_env"]:
                    group_style = [f"group size={geom[1]} by {geom[0]}"]
                    group_style.extend(data["extra groupstyle options [base]"])
                    options = ["group style={{{}}}".format(", ".join(group_style))]
                    self.content.append(
                        data["flavor"].start("groupplot") + f"[{', '.join(options)}]"
                    )
                    data["is_in_groupplot_env"] = True
                    data["pgfplots libs"].add("groupplots")

    def _get_label_rotation_and_horizontal_alignment(self, obj, data, x_or_y):
        label_style = ""

        major_tick_labels = (
            obj.xaxis.get_majorticklabels()
            if x_or_y == "x"
            else obj.yaxis.get_majorticklabels()
        )

        if not major_tick_labels:
            return None

        tick_label_text_width_identifier = f"{x_or_y} tick label text width"
        if tick_label_text_width_identifier in self.axis_options:
            self.axis_options.remove(tick_label_text_width_identifier)

        values = []

        tick_labels_rotation = [label.get_rotation() for label in major_tick_labels]
        if len(set(tick_labels_rotation)) == 1:
            if tick_labels_rotation[0] != 0:
                values.append(f"rotate={tick_labels_rotation[0]}")
        else:
            values.append(
                "rotate={{{},0}}[\\ticknum]".format(
                    ",".join(str(x) for x in tick_labels_rotation)
                )
            )

        tick_labels_horizontal_alignment = [
            label.get_horizontalalignment() for label in major_tick_labels
        ]
        if len(set(tick_labels_horizontal_alignment)) == 1:
            anchor = {"right": "east", "left": "west", "center": "center"}[
                tick_labels_horizontal_alignment[0]
            ]
            if not (x_or_y == "x" and anchor == "center") and not (
                x_or_y == "y" and anchor == "east"
            ):
                values.append(f"anchor={anchor}")

        if values:
            label_style = "{}ticklabel style={{{}}}".format(x_or_y, ",".join(values))

        # Ignore horizontal alignment if no '{x,y} tick label text width' has been
        # passed in the 'extra' parameter
        # if tick_label_text_width:
        #     if is_tick_label_alignment_identical:
        #         values.append(f"align={tick_labels_horizontal_alignment[0]}")
        #         values.append(f"text width={tick_label_text_width}")
        #     else:
        #         for idx, x in enumerate(tick_labels_horizontal_alignment):
        #             label_style += f"{x_or_y}_tick_label_ha_{idx}/.initial = {x}"

        #         values.append(
        #             f"align=\\pgfkeysvalueof{{/pgfplots/{x_or_y}_tick_label_ha_\\ticknum}}"
        #         )
        #         values.append(f"text width={tick_label_text_width}")

        # label_style = (
        #     "every {} tick label/.style = {{\n"
        #     "{}\n"
        #     "}}".format(x_or_y, ",\n".join(values))
        # )

        return label_style

def _get_tick_position(obj, axes_obj):
    major_ticks = obj.xaxis.majorTicks if axes_obj == "x" else obj.yaxis.majorTicks

    major_ticks_bottom = [tick.tick1line.get_visible() for tick in major_ticks]
    major_ticks_top = [tick.tick2line.get_visible() for tick in major_ticks]

    major_ticks_bottom_show_all = False
    if len(set(major_ticks_bottom)) == 1 and major_ticks_bottom[0] is True:
        major_ticks_bottom_show_all = True

    major_ticks_top_show_all = False
    if len(set(major_ticks_top)) == 1 and major_ticks_top[0] is True:
        major_ticks_top_show_all = True

    major_ticks_position = None
    if not major_ticks_bottom_show_all and not major_ticks_top_show_all:
        position_string = f"{axes_obj}majorticks=false"
    elif major_ticks_bottom_show_all and major_ticks_top_show_all:
        major_ticks_position = "both"
    elif major_ticks_bottom_show_all:
        major_ticks_position = "left"
    elif major_ticks_top_show_all:
        major_ticks_position = "right"

    if major_ticks_position:
        position_string = f"{axes_obj}tick pos={major_ticks_position}"

    return position_string, major_ticks_position

def _get_ticks(data, xy, ticks, ticklabels):
    """Gets a {'x','y'}, a number of ticks and ticks labels, and returns the
    necessary axis options for the given configuration.
    """
    axis_options = []
    pgfplots_ticks = []
    pgfplots_ticklabels = []
    is_label_required = False
    for tick, ticklabel in zip(ticks, ticklabels):
        # store the label anyway
        label = ticklabel.get_text()
        if "," in label:
            label = "{" + label + "}"
        if ticklabel.get_visible():
            label = _common_texification(label)
            pgfplots_ticklabels.append(label)
        else:
            is_label_required = True
        # Check if the label is necessary. If one of the labels is, then all of them
        # must appear in the TikZ plot.
        if label:
            try:
                label_float = float(label.replace("\N{MINUS SIGN}", "-"))
                is_label_required = is_label_required or (label and label_float != tick)
            except ValueError:
                is_label_required = True
    # note: ticks may be present even if labels are not, keep them for grid lines
    for tick in ticks:
        pgfplots_ticks.append(tick)

    # if the labels are all missing, then we need to output an empty set of labels
    if len(ticklabels) == 0 and len(ticks) != 0:
        axis_options.append(f"{xy}ticklabels={{}}")
        # remove the multiplier too
        axis_options.append(f"scaled {xy} ticks=" + r"manual:{}{\pgfmathparse{#1}}")

    # Leave the ticks to PGFPlots if not in STRICT mode and if there are no explicit
    # labels.
    if data["strict"] or is_label_required:
        if pgfplots_ticks:
            ff = data["float format"]
            axis_options.append(
                "{}tick={{{}}}".format(
                    xy, ",".join([f"{el:{ff}}" for el in pgfplots_ticks])
                )
            )
        else:
            val = "{}" if "minor" in xy else "\\empty"
            axis_options.append(f"{xy}tick={val}")

        if is_label_required:
            length = sum(len(label) for label in pgfplots_ticklabels)
            sep = ("", ",", "") if length < 75 else ("\n  ", ",\n  ", "\n")
            string = sep[1].join(pgfplots_ticklabels)
            axis_options.append(f"{xy}ticklabels={{{sep[0]}{string}{sep[2]}}}")
    return axis_options

def _is_colorbar_heuristic(obj):
    """Find out if the object is in fact a color bar."""
    # TODO come up with something more accurate here
    # Might help:
    # TODO Are the colorbars exactly the l.collections.PolyCollection's?
    try:
        aspect = float(obj.get_aspect())
    except ValueError:
        # e.g., aspect == 'equal'
        return False

    # Assume that something is a colorbar if and only if the ratio is above 5.0
    # and there are no ticks on the corresponding axis. This isn't always true,
    # though: The ratio of a color can be freely adjusted by the aspect
    # keyword, e.g.,
    #
    #    plt.colorbar(im, aspect=5)
    #
    limit_ratio = 5.0

    return (aspect >= limit_ratio and len(obj.get_xticks()) == 0) or (
        aspect <= 1.0 / limit_ratio and len(obj.get_yticks()) == 0
    )

def _mpl_cmap2pgf_cmap(cmap, data):
    """Converts a color map as given in matplotlib to a color map as
    represented in PGFPlots.
    """
    if isinstance(cmap, mpl.colors.LinearSegmentedColormap):
        return _handle_linear_segmented_color_map(cmap, data)

    assert isinstance(
        cmap, mpl.colors.ListedColormap
    ), "Only LinearSegmentedColormap and ListedColormap are supported"
    return _handle_listed_color_map(cmap, data)

def _handle_linear_segmented_color_map(cmap, data):
    assert isinstance(cmap, mpl.colors.LinearSegmentedColormap)

    if cmap.is_gray():
        is_custom_colormap = False
        return ("blackwhite", is_custom_colormap)

    # For an explanation of what _segmentdata contains, see
    # http://matplotlib.org/mpl_examples/pylab_examples/custom_cmap.py
    # A key sentence:
    # If there are discontinuities, then it is a little more complicated.  Label the 3
    # elements in each row in the cdict entry for a given color as (x, y0, y1).  Then
    # for values of x between x[i] and x[i+1] the color value is interpolated between
    # y1[i] and y0[i+1].
    segdata = cmap._segmentdata
    red = segdata["red"]
    green = segdata["green"]
    blue = segdata["blue"]

    # Loop over the data, stop at each spot where the linear interpolations is
    # interrupted, and set a color mark there.
    #
    # Set initial color.
    k_red = 0
    k_green = 0
    k_blue = 0
    colors = []
    X = []
    while True:
        # find next x
        x = min(red[k_red][0], green[k_green][0], blue[k_blue][0])

        if red[k_red][0] == x:
            red_comp = red[k_red][1]
            k_red += 1
        else:
            red_comp = _linear_interpolation(
                x,
                (red[k_red - 1][0], red[k_red][0]),
                (red[k_red - 1][2], red[k_red][1]),
            )

        if green[k_green][0] == x:
            green_comp = green[k_green][1]
            k_green += 1
        else:
            green_comp = _linear_interpolation(
                x,
                (green[k_green - 1][0], green[k_green][0]),
                (green[k_green - 1][2], green[k_green][1]),
            )

        if blue[k_blue][0] == x:
            blue_comp = blue[k_blue][1]
            k_blue += 1
        else:
            blue_comp = _linear_interpolation(
                x,
                (blue[k_blue - 1][0], blue[k_blue][0]),
                (blue[k_blue - 1][2], blue[k_blue][1]),
            )

        X.append(x)
        colors.append((red_comp, green_comp, blue_comp))

        if x >= 1.0:
            break

    # The PGFPlots color map has an actual physical scale, like (0cm,10cm), and the
    # points where the colors change is also given in those units. As of now
    # (2010-05-06) it is crucial for PGFPlots that the difference between two successive
    # points is an integer multiple of a given unity (parameter to the colormap; e.g.,
    # 1cm).  At the same time, TeX suffers from significant round-off errors, so make
    # sure that this unit is not too small such that the round- off errors don't play
    # much of a role. A unit of 1pt, e.g., does most often not work.
    unit = "pt"

    # Scale to integer (too high integers will firstly be slow and secondly may produce
    # dimension errors or memory errors in latex)
    # 0-1000 is the internal granularity of PGFplots.
    # 16300 was the maximum value for pgfplots<=1.13
    X = _scale_to_int(np.array(X), 1000)

    color_changes = []
    ff = data["float format"]
    for k, x in enumerate(X):
        color_changes.append(
            f"rgb({x}{unit})="
            f"({colors[k][0]:{ff}},{colors[k][1]:{ff}},{colors[k][2]:{ff}})"
        )

    colormap_string = "{{mymap}}{{[1{}]\n  {}\n}}".format(
        unit, ";\n  ".join(color_changes)
    )
    is_custom_colormap = True
    return (colormap_string, is_custom_colormap)

def _handle_listed_color_map(cmap, data):
    assert isinstance(cmap, mpl.colors.ListedColormap)

    # check for predefined colormaps in both matplotlib and pgfplots
    from matplotlib import pyplot as plt

    cm_translate = {
        # All the rest are LinearSegmentedColorMaps. :/
        # 'autumn': 'autumn',
        # 'cool': 'cool',
        # 'copper': 'copper',
        # 'gray': 'blackwhite',
        # 'hot': 'hot2',
        # 'hsv': 'hsv',
        # 'jet': 'jet',
        # 'spring': 'spring',
        # 'summer': 'summer',
        "viridis": "viridis",
        # 'winter': 'winter',
    }
    for mpl_cm, pgf_cm in cm_translate.items():
        if cmap.colors == plt.get_cmap(mpl_cm).colors:
            is_custom_colormap = False
            return (pgf_cm, is_custom_colormap)

    unit = "pt"
    ff = data["float format"]
    if cmap.N is None or cmap.N == len(cmap.colors):
        colors = [
            f"rgb({k}{unit})=({rgb[0]:{ff}},{rgb[1]:{ff}},{rgb[2]:{ff}})"
            for k, rgb in enumerate(cmap.colors)
        ]
    else:
        reps = int(float(cmap.N) / len(cmap.colors) - 0.5) + 1
        repeated_cols = reps * cmap.colors
        colors = [
            f"rgb({k}{unit})=({rgb[0]:{ff}},{rgb[1]:{ff}},{rgb[2]:{ff}})"
            for k, rgb in enumerate(repeated_cols[: cmap.N])
        ]
    colormap_string = "{{mymap}}{{[1{}]\n {}\n}}".format(unit, ";\n  ".join(colors))
    is_custom_colormap = True
    return (colormap_string, is_custom_colormap)

def _scale_to_int(X, max_val):
    """Scales the array X such that it contains only integers."""
    # if max_val is None:
    #     X = X / _gcd_array(X)
    X = X / max(1 / max_val, _gcd_array(X))
    return [int(entry) for entry in X]

def _gcd_array(X):
    """
    Return the largest real value h such that all elements in x are integer
    multiples of h.
    """
    greatest_common_divisor = 0.0
    for x in X:
        greatest_common_divisor = _gcd(greatest_common_divisor, x)

    return greatest_common_divisor

def _gcd(a, b):
    """Euclidean algorithm for calculating the GCD of two numbers a, b.
    This algoritm also works for real numbers:
    Find the greatest number h such that a and b are integer multiples of h.
    """
    # Keep the tolerance somewhat significantly above machine precision as otherwise
    # round-off errors will be accounted for, returning 1.0e-10 instead of 1.0 for the
    # values
    #   [1.0, 2.0000000001, 3.0, 4.0].
    while a > 1.0e-5:
        a, b = b % a, a
    return b

def _linear_interpolation(x, X, Y):
    """Given two data points [X,Y], linearly interpolate those at x."""
    return (Y[1] * (x - X[0]) + Y[0] * (X[1] - x)) / (X[1] - X[0])

def _find_associated_colorbar(obj):
    """A rather poor way of telling whether an axis has a colorbar associated: Check the
    next axis environment, and see if it is de facto a color bar; if yes, return the
    color bar object.
    """
    for child in obj.get_children():
        try:
            cbar = child.colorbar
        except AttributeError:
            continue
        if cbar is not None:  # really necessary?
            # if fetch was successful, cbar contains
            # (reference to colorbar, reference to axis containing colorbar)
            return cbar
    return None

def _try_f2i(x):
    """If possible, convert float to int without rounding.
    Used for log base: if not used, base for log scale can be "10.0" (and then
    printed as such  by pgfplots).
    """
    return int(x) if int(x) == x else x

max_stars_repo_name: vthuongt/tikzplotlib

max_issues_repo_name: vthuongt/tikzplotlib

max_forks_repo_name: vthuongt/tikzplotlib


text:
# -*- coding: utf-8 -*-
"""
Created on Sun Apr  7 17:28:57 2019

@author: abhishekt
"""

"""
   plot_confusion_matrix: Plot the confusion matrix

   Input:  pred_rating - list of predictions
           actual_rating - list of actual target values
   Returns: None
"""
from sklearn.metrics import confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sn
import numpy as np

def plot_confusion_matrix(title, pred_rating, actual_rating):
    #class_names
    classes = ['0','1','2','3','4','5','6','7','8','9']

    #cnf_matrix = confusion_matrix(y_test, y_pred)
    cnf_matrix = confusion_matrix(np.asarray(actual_rating), np.asarray(pred_rating))
    #print(cnf_matrix)
    plt.figure(figsize = (10,6))
    sn.heatmap(cnf_matrix, annot=True)

    plt.title(title)
    plt.xlabel('Predicted label')
    plt.ylabel('Actual label')

    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45)
    plt.yticks(tick_marks, classes, rotation=90)

    plt.show()

"""
   plot_Accuracy: Function to plot Validation & Test Accuracy
"""
def plot_Accuracy(num_hidden_units, train_accu, test_accu, train_label, test_label):

    #create figure object for plotting
    plt.figure(figsize=(7, 5))

    #plot accuracy on Y-axis & log(C) on x-axis
    plt.plot(num_hidden_units, train_accu, color = 'blue', label = train_label)
    plt.plot(num_hidden_units, test_accu, color = 'green', label = test_label)
    plt.title('Train & Test Accuracy Plot')
    plt.xlabel("Number of hidden units (constant for each hidden layer) -->")
    plt.ylabel("Prediction Accuracies")
    plt.legend(loc='lower right', shadow=True, fontsize='large', numpoints=1)
    plt.show()

def relu(Z):
    """
    Implement the RELU function.
    Arguments:
    Z -- Output of the linear layer, of any shape
    Returns:
    A -- Post-activation parameter, of the same shape as Z
    cache -- a python dictionary containing "A" ; stored for computing the backward pass efficiently
    """

    A = np.maximum(0,Z)

    assert(A.shape == Z.shape)

    cache = Z
    return A, cache

def relu_backward(dA, cache):
    """
    Implement the backward propagation for a single RELU unit.
    Arguments:
    dA -- post-activation gradient, of any shape
    cache -- 'Z' where we store for computing backward propagation efficiently
    Returns:
    dZ -- Gradient of the cost with respect to Z
    """

    Z = cache
    dZ = np.array(dA, copy=True) # just converting dz to a correct object.

    # When z <= 0, you should set dz to 0 as well.
    dZ[Z <= 0] = 0

    assert (dZ.shape == Z.shape)

    return dZ

ext: py

reason: import numpy


text:
import math
import cv2
import numpy as np
import sys
from scipy import signal

def pascalSmooth(n):
    pascalSmooth = np.zeros([1,n],np.float32)
    for i in range(n):
        pascalSmooth[0][i] = math.factorial(n-1)/(math.factorial(i)*math.factorial(n-1-i))
    return pascalSmooth

def pascalDiff(n):
    pascalDiff = np.zeros([1,n], np.float32)
    pascalSmooth_previous = pascalSmooth(n-1)
    for i in range(n):
        if i==0:
            pascalDiff[0][i] = pascalSmooth_previous[0][i]
        elif i == n-1:
            pascalDiff[0][i] = -pascalSmooth_previous[0][i-1]
        else:
            pascalDiff[0][i] = pascalSmooth_previous[0][i] - pascalSmooth_previous[0][i-1]

    return pascalDiff

def getSobeKernel(n):
    pascalSmoothKernel = pascalSmooth(n)
    pascalDiffKernel = pascalDiff(n)

    sobelKernel_x = singal.convolve2d(pascalSmoothKernel.transpose(), pascalDiffKernel, mode='full')
    sobelKernel_y = singal.convolve2d(pascalSmoothKernel, pascalDiffKernel.transpose(), mode='full')
    return (sobelKernel_x, sobelKernel_y)

def sobel(image, n):
    wors,cols =  image.shape
    pascalSmoothKernel = pascalSmooth(n)
    pascalDiffKernel = pascalDiff(n)

    image_sobel_x = signal.convolve2d(image, pascalSmoothKernel.transpose(), mode='same')
    image_sobel_x = signal.convolve2d(image_sobel_x,pascalDiffKernel,mode='same')

    image_sobel_y = signal.convolve2d(image, pascalSmoothKernel, mode='same')
    image_sobel_y = signal.convolve2d(image_sobel_y, pascalDiffKernel.transpose(), mode='same')

    return (image_sobel_x, image_sobel_y)

if __name__ == "__main__":
    if len(sys.argv) > 1:
        image = cv2.imread(sys.argv[1], cv2.IMREAD_GRAYSCALE)
    else:
        print("usage: python sobel.py imageFile")
        sys.exit(0)

    image_sobel_x,image_sobel_y = sobel(image,7)

    edge = np.sqrt(np.power(image_sobel_x,2.0)+np.power(image_sobel_y,2.0))

    edge = edge/np.max(edge)
    edge = np.power(edge,1)
    edge*=255
    edge = edge.astype(np.uint8)
    cv2.imshow('sobel edge', edge)
    cv2.imwrite('sobel.jpg',edge)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

max_stars_repo_name: yugengde/img-process

max_issues_repo_name: yugengde/img-process

max_forks_repo_name: yugengde/img-process

reason: import numpy,from scipy


text:
from abc import ABC
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.statespace.sarimax import SARIMAX

# Construct the model
class StateSpaceModel(sm.tsa.statespace.MLEModel, ABC):
    def __init__(self, endog, exog, factors_x, factors_y):
        # Initialize the state space model
        endog_aug = pd.concat([exog, endog.to_frame()], axis=1)
        k_states = k_posdef = factors_x + factors_y
        super(StateSpaceModel, self).__init__(endog_aug, k_states=k_states, k_posdef=k_posdef,
                                              initialization='approximate_diffuse')
        self._covariates = endog_aug.columns
        self._factors_x = factors_x
        self._factors_y = factors_y
        self._param_cov_state_f_idx = self._params_cov_state_z_idx = self._params_cov_obs_idx = None
        self._params_phi_1_idx = self._params_phi_23_idx = self._params_phi_23_idx = None
        self._params_a_1_idx = self._params_a_2_idx = None

        # Setup the fixed components of the state space representation
        transition_f = np.hstack((np.eye(factors_x), np.zeros((factors_x, factors_y))))
        transition_z = np.hstack((np.zeros((factors_y, factors_x)), np.eye(factors_y)))
        transition = np.vstack((transition_f, transition_z))
        dims_x = endog_aug.shape[1] - 1
        dims_y = 1
        self._dims_x = dims_x
        self._dims_y = dims_y
        # Assume [x, y]'
        design_x = np.hstack((np.ones((dims_x, factors_x)), np.zeros((dims_x, factors_y))))
        design_y = np.hstack((np.zeros((dims_y, factors_x)), np.ones((dims_y, factors_y))))
        design = np.vstack((design_x, design_y))

        self.ssm['design'] = design.reshape((dims_x + 1, k_states, 1))
        self.ssm['transition'] = transition.reshape((k_states, k_states, 1))
        self.ssm['selection'] = np.eye(k_states)
        self.ssm['obs_intercept'] = np.zeros(dims_x + dims_y).reshape(-1, 1)

        # Cache some indices
        self._state_cov_idx = np.diag_indices(k_posdef)
        self._obs_cov_idx = np.diag_indices(dims_x + dims_y)

        # grid_transition_f = (np.repeat(np.arange(factors_x), factors_x),
        #                      np.tile(np.arange(factors_x), factors_x))
        grid_transition_f = np.diag_indices(factors_x)
        grid_transition_z = (np.repeat(np.arange(factors_x, k_states), k_states),
                             np.tile(np.arange(k_states), factors_y))
        self._transition_f_idx = grid_transition_f
        self._transition_z_idx = grid_transition_z

        grid_design_x = (np.repeat(np.arange(dims_x), factors_x),
                         np.tile(np.arange(factors_x), dims_x))
        grid_design_y = (np.repeat(np.arange(dims_x, dims_x + dims_y), factors_y),
                         np.tile(np.arange(factors_x, k_states), dims_y))
        self._design_x_idx = grid_design_x
        self._design_y_idx = grid_design_y

        self.init_param_indx()

    @staticmethod
    def get_position(idx, i, row_offset=0, col_offset=0):
        return idx[0][i]-row_offset, idx[1][i]-col_offset

    def init_param_indx(self):
        c = 0
        params_cov_obs = ['sigma2.%s' % i for i in self._covariates]
        self._params_cov_obs_idx = (c, c + len(params_cov_obs))
        c += len(params_cov_obs)
        params_cov_state_f = ['sigma2.f.%i' % i for i in range(self._factors_x)]
        self._param_cov_state_f_idx = (c, c + len(params_cov_state_f))
        c += len(params_cov_state_f)
        params_cov_state_z = ['sigma2.z.%i' % i for i in range(self._factors_y)]
        self._params_cov_state_z_idx = (c, c + len(params_cov_state_z))
        c += len(params_cov_state_z)
        params_cov = params_cov_state_f + params_cov_state_z + params_cov_obs

        params_phi_1 = ['phi.1.%i%i' % self.get_position(self._transition_f_idx, i) for i in range(len(self._transition_f_idx[0]))]
        self._params_phi_1_idx = (c, c+len(params_phi_1))
        c += len(params_phi_1)

        params_phi_23 = ['phi.23.%i%i' % self.get_position(self._transition_z_idx, i,
                                                           row_offset=self._factors_x) for i in range(len(self._transition_z_idx[0]))]
        self._params_phi_23_idx = (c, c + len(params_phi_23))
        c += len(params_phi_23)
        params_phi = params_phi_1 + params_phi_23

        params_a_1 = ['a.1.%i%i' % self.get_position(self._design_x_idx, i) for i in range(len(self._design_x_idx[0]))]
        self._params_a_1_idx = (c, c+len(params_a_1))
        c += len(params_a_1)
        params_a_2 = ['a.2.%i%i' % self.get_position(self._design_y_idx, i,
                                                     row_offset=self._dims_x,
                                                     col_offset=self._factors_x) for i in range(len(self._design_y_idx[0]))]
        self._params_a_2_idx = (c, c+len(params_a_2))
        c += len(params_a_2)
        params_a = params_a_1 + params_a_2

        return params_cov + params_phi + params_a

    @property
    def param_names(self):
        return self.init_param_indx()

    # Describe how parameters enter the model
    def update(self, params, *args, **kwargs):
        params = super(StateSpaceModel, self).update(params, *args, **kwargs)
        # Observation covariance
        self.ssm[('obs_cov',) + self._obs_cov_idx] = params[self._params_cov_obs_idx[0]:self._params_cov_obs_idx[1]]
        # State covariance
        self.ssm[('state_cov',) + self._state_cov_idx] = params[self._param_cov_state_f_idx[0]:self._params_cov_state_z_idx[1]]
        # Transition matrix
        self.ssm[('transition',) + self._transition_f_idx] = params[self._params_phi_1_idx[0]:self._params_phi_1_idx[1]]
        self.ssm[('transition',) + self._transition_z_idx] = params[self._params_phi_23_idx[0]:self._params_phi_23_idx[1]]
        # Design matrix
        self.ssm[('design',) + self._design_x_idx] = params[self._params_a_1_idx[0]:self._params_a_1_idx[1]]
        self.ssm[('design',) + self._design_y_idx] = params[self._params_a_2_idx[0]:self._params_a_2_idx[1]]

    # Specify start parameters and parameter names
    @property
    def start_params(self):
        design, obs_cov, state_cov, transition = self.generate_start_matrices()

        params_state_cov = state_cov[self._state_cov_idx]
        params_obs_cov = obs_cov[self._obs_cov_idx]

        params_phi = np.concatenate((transition[self._transition_f_idx],
                                     transition[self._transition_z_idx]), axis=0)
        params_a = np.concatenate((design[self._design_x_idx],
                                   design[self._design_y_idx]), axis=0)

        return np.concatenate((params_obs_cov, params_state_cov, params_phi, params_a))

    def generate_start_matrices(self):
        _exog = pd.DataFrame(self.endog[:, :-1], columns=self._covariates[:-1]).interpolate().fillna(0)
        _endog = pd.Series(self.endog[:, -1], name=self._covariates[-1]).interpolate().fillna(0)
        cov = _exog.cov()
        w, v = np.linalg.eig(cov)
        factors = pd.DataFrame(np.dot(_exog, v[:, :self._factors_x]), index=_exog.index)
        _model = SARIMAX(endog=_endog, exog=factors, order=(self._factors_y, 0, 0))
        res = _model.fit(disp=False, maxiter=100)
        params_arx = res.params

        phi1 = np.eye(self._factors_x)
        factors_coeff = params_arx.values[:self._factors_x].reshape(1, -1)
        ar_coeff = params_arx.values[self._factors_x:-1].reshape(1, -1)
        phi2 = np.vstack([factors_coeff, np.zeros((self._factors_y - 1, self._factors_x))])
        phi3 = np.vstack([ar_coeff, np.eye(self._factors_y)[:-1, :]])
        transition = np.vstack([np.hstack([phi1, np.zeros((self._factors_x, self._factors_y))]),
                                np.hstack([phi2, phi3])])

        a1 = v.T[:, :self._factors_x]
        a2 = np.eye(self._dims_y, self._factors_y)
        design_x = np.hstack([a1, np.zeros((self._dims_x, self._factors_y))])
        design_y = np.hstack([np.zeros((self._dims_y, self._factors_x)), a2])
        design = np.vstack([design_x, design_y])

        state_cov = np.eye(self.k_states)

        obs_cov = np.eye(len(self._covariates))
        obs_cov[-1, -1] = params_arx.values[-1]

        return design, obs_cov, state_cov, transition

    def transform_params(self, unconstrained):
        constrained = unconstrained.copy()
        for i1, i2 in [self._param_cov_state_f_idx, self._params_cov_state_z_idx, self._params_cov_obs_idx]:
            constrained[i1:i2] = unconstrained[i1:i2] ** 2
        return constrained

    def untransform_params(self, constrained):
        unconstrained = constrained.copy()
        for i1, i2 in [self._param_cov_state_f_idx, self._params_cov_state_z_idx, self._params_cov_obs_idx]:
            unconstrained[i1:i2] = constrained[i1:i2] ** 0.5
        return unconstrained

if __name__ == "__main__":
    import os
    from src.d01_data.dengue_data_api import DengueDataApi
    from src.d04_modeling.abstract_sm import AbstractSM
    os.chdir('../')
    dda = DengueDataApi(interpolate=False)
    x1, x2, y1, y2 = dda.split_data(random=False)

    factors_x = 3
    factors_y = 3
    abstract_model = AbstractSM(x_train=x1, y_train=y1, bias=False)
    city = 'sj'

    endog, exog = dda.format_data(x1.loc[city].copy(), y1.loc[city].copy(), interpolate=False)

    endog_mean = endog.mean(axis=0)
    endog_std = endog.std(axis=0)
    exog_mean = exog.mean(axis=0)
    exog_std = exog.std(axis=0)

    endog = (endog - endog_mean) / endog_std
    exog = (exog - exog_mean) / exog_std

    print("-------------- DFM Model --------------")
    model_dfmq = sm.tsa.DynamicFactorMQ(exog,
                                        factors=factors_x,
                                        factor_orders=1,
                                        idiosyncratic_ar1=False)
    results_dfmq = model_dfmq.fit(method='em')

    state_names = pd.Index(['f%i' % i for i in range(factors_x)])
    transition_df = pd.DataFrame(model_dfmq.ssm['transition'], index=state_names, columns=state_names)
    print(transition_df.round(4).to_string())
    design_df = pd.DataFrame(model_dfmq.ssm['design'],
                             index=exog.columns,
                             columns=state_names)
    print(design_df.round(4).to_string())

    print("-------------- SSM Model --------------")

    model = StateSpaceModel(endog=endog, exog=exog, factors_x=factors_x, factors_y=factors_y)
    model.update(model.start_params)
    state_names = pd.Index(['f%i' % i for i in range(factors_x)]).append(pd.Index(['z%i' % i for i in range(factors_y)]))
    transition_df = pd.DataFrame(model.ssm['transition'], index=state_names, columns=state_names)
    print(transition_df.round(4).to_string())
    design_df = pd.DataFrame(model.ssm['design'],
                             index=exog.columns.append(pd.Index([endog.name])),
                             columns=state_names)
    print(design_df.round(4).to_string())
    results = model.fit(maxiter=200)
    transition_df = pd.DataFrame(model.ssm['transition'], index=state_names, columns=state_names)
    print(transition_df.round(4).to_string())
    design_df = pd.DataFrame(model.ssm['design'],
                             index=exog.columns.append(pd.Index([endog.name])),
                             columns=state_names)
    print(design_df.round(4).to_string())

reason: import numpy,import statsmodels,from statsmodels


text:
"""main client script
controls microphone stream and organise all audio responses"""

import socket
import pickle
import pyaudio
import numpy as np
import concurrent.futures
# from sound import SoundBot
from random import random
import trio
from time import sleep

class Client:
    def __init__(self):
        self.running = True
        self.connected = False
        self.logging = False

        self.HOST = '127.0.0.1'  # Client IP (this)
        self.PORT = 65432
        # Port to listen on (non-privileged ports are > 1023)

        self.CHUNK = 2 ** 11
        self.RATE = 44100
        self.p = pyaudio.PyAudio()
        self.stream = self.p.open(format=pyaudio.paInt16,
                                  channels=1,
                                  rate=self.RATE,
                                  input=True,
                                  frames_per_buffer=self.CHUNK)

        # build send data dict
        self.send_data_dict = {'mic_level': 0,
                               'speed': 1,
                               'tempo': 0.1
                               }

        # init got dict
        self.got_dict = {}

    def snd_listen(self):
        print("mic listener: started!")
        while self.running:
            data = np.frombuffer(self.stream.read(self.CHUNK),
                                 dtype=np.int16)
            peak = np.average(np.abs(data)) * 2
            if peak > 2000:
                bars = "#" * int(50 * peak / 2 ** 16)
                print("%05d %s" % (peak, bars))
            self.send_data_dict['mic_level'] = peak / 30000

    def terminate(self):
        self.stream.stop_stream()
        self.stream.close()
        self.p.terminate()

    # def mincer(self, got_AI_data, rhythm_rate):
    #     print(f'in the mincer ===== {got_AI_data, rhythm_rate}')
    #     # tic = time()
    #     rnd_dur = random()
    #     duration = rnd_dur + rhythm_rate
    #
    #     # make a sound at calc duration
    #     self.snd.play_sound(got_AI_data, duration)
    #     # toc = time()

    # async def left(self):
    #     print('left bot: started')
    #     while self.connected:
    #         # get latest data
    #         left_master_data = self.got_dict['master_output']
    #         left_rhythm_rate = self.got_dict['rhythm_rate']
    #
    #         # calc temp timing
    #         rnd_dur = random()
    #         left_duration = rnd_dur + left_rhythm_rate
    #
    #         # make a sound at calc duration
    #         self.bot_left.play_sound(left_master_data, left_duration)
    #
    # async def right(self):
    #     print('right bot: started')
    #     while self.connected:
    #         # get latest data
    #         right_master_data = self.got_dict['master_output']
    #         right_rhythm_rate = self.got_dict['rhythm_rate']
    #
    #         # calc temp timing
    #         rnd_dur = random()
    #         right_duration = rnd_dur + right_rhythm_rate
    #
    #         # make a sound at calc duration
    #         self.bot_right.play_sound(right_master_data, right_duration)

    def client(self):
        print("client: started!")
        while self.running:
            print(f"client: connecting to {self.HOST}:{self.PORT} ..... start engine_server in iTERM NOW!!!!!!")
            with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
                s.bind((self.HOST, self.PORT))
                s.listen()
                client_stream, addr = s.accept()
                with client_stream:
                    print('Connected by', addr)
                    self.connected = True
                    while self.connected:
                        # get data from stream
                        data = client_stream.recv(1024)
                        self.got_dict = pickle.loads(data)
                        if self.logging:
                            print(f"receiver: got data {self.got_dict}")

                        # send it to the mincer for soundBot control
                        # NB play_with_simpleaudio does not hold thread
                        # master_data = self.got_dict['master_output']
                        # rhythm_rate = self.got_dict['rhythm_rate']
                        # self.mincer(master_data, rhythm_rate)

                        # send out-going data to server
                        send_data = pickle.dumps(self.send_data_dict, -1)
                        client_stream.sendall(send_data)

    # async def parent(self):
    #     print("parent: started!")
    #     while self.connected:
    #         async with trio.open_nursery() as nursery:
    #             # spawning left independent soundbot
    #             print("parent: spawning left bot ...")
    #             nursery.start_soon(self.left)
    #
    #             # spawning right independent soundbot
    #             print("parent: spawning right bot ...")
    #             nursery.start_soon(self.right)

    def parent_go(self):
        while self.running:
            if not self.connected:
                sleep(1)
            else:
                trio.run(self.parent)

    def main(self):
        # snd_listen and client need dependent threads.
        # All other IO is ok as a single Trio thread inside self.client
        tasks = [self.snd_listen, self.client]

        with concurrent.futures.ThreadPoolExecutor(max_workers=6) as executor:
            futures = {executor.submit(task): task for task in tasks}

if __name__ == '__main__':
    cl = Client()
    cl.main()


text:
# coding: utf-8
# created by deng on 7/30/2018

import numpy as np
import scipy as sp

import term_weighting_model.calc_weightings as cw
from utils.data_util import load_to_df
from utils.path_util import from_project_root

from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer
from sklearn.externals import joblib

MAX_FEATURES = 4000000
MIN_DF = 3
MAX_DF = 0.8
MAX_N = 3
TW_TYPE = 'dc'
COLUMN = 'word_seg'
TRAIN_URL = from_project_root('data/train_set.csv')

def transform(train_url=TRAIN_URL, test_url=None, column='word_seg', tw_type=TW_TYPE):
    """

    Args:
        column: column to use
        train_url: str, url to train data (with header)
        test_url: url to test data
        tw_type: str, term wighting type {idf, dc, bdc}

    Returns:
        X, y, X_test: vertorized data

    """
    data_url = from_project_root("processed_data\phrase_level_data.csv")
    if column == 'article':
        data_url = data_url.replace('phrase', 'word')

    if tw_type == 'idf':
        return tfidf_transform(train_url, test_url, column=column)
    elif tw_type == 'dc':
        dc_dict = cw.calc_dc(data_url, ngram=MAX_N)
        return dict_transform(dc_dict, train_url, test_url, column=column)
    elif tw_type == 'bdc':
        bdc_dict = cw.calc_bdc(data_url, ngram=MAX_N)
        return dict_transform(bdc_dict, train_url, test_url, column=column)

def tfidf_transform(train_url, test_url, column='word_seg', sublinear_tf=True, max_n=MAX_N,
                    min_df=MIN_DF, max_df=MAX_DF, max_features=MAX_FEATURES):
    """ vectorize use TfidfVectorizer

    Args:
        train_url: url to train data
        test_url: url to test data
        column: column to use
        sublinear_tf:
        max_n:
        min_df:
        max_df:
        max_features:

    Returns:
        X, X_test, y: vectorized data

    """
    # set token_pattern(default: (?u)\b\w\w+\b' to keep single char tokens
    vectorizer = TfidfVectorizer(min_df=min_df, max_df=max_df, max_features=max_features, ngram_range=(1, max_n),
                                 sublinear_tf=sublinear_tf, token_pattern='(?u)\w+')
    train_df = load_to_df(train_url)
    X = vectorizer.fit_transform(train_df[column])
    y = np.asarray(train_df['class'])
    X_test = None
    if test_url:
        X_test = vectorizer.transform(load_to_df(test_url)[column])
    return X, y, X_test

def dict_transform(tw_dict, train_url=TRAIN_URL, test_url=None, column='word_seg',
                   max_n=MAX_N, min_df=MIN_DF, max_df=MAX_DF, max_features=MAX_FEATURES,
                   normalize=True, sublinear_tf=True, re_weight=0):
    """ use offline dict to transform data into vector

    Args:
        train_url: url to train data (with header)
        test_url: url to test data (with header)
        sentences: list of sentence to be vectorized
        tw_dict: term weighting dict to use
        max_n: max_n for CountVectorizer
        min_df: min_df for CountVectorizer
        max_df: max_df for CountVectorizer
        normalize: normalize the vector or not
        sublinear_tf: use 1 + log(tf) instead of tf
        max_features: max_features for CountVectorizer
        re_weight: if re_weight > 0, use (1-re_weight) + (re_weight) * weights instead of weight
        column: column to use in dataframe

    Returns:
        X, y, X_test: vectorized data

    """
    print("transforming...")
    train_df = load_to_df(train_url)
    vectorizer = CountVectorizer(min_df=min_df, max_df=max_df, ngram_range=(1, max_n),
                                 token_pattern='(?u)\w+', max_features=max_features)
    X_train = vectorizer.fit_transform(train_df[column])  # use train data to get vocab
    y_train = np.asarray(train_df['class'])
    X_test = vectorizer.transform(load_to_df(test_url)[column]) if test_url else None

    # get words of all columns represent to
    words = vectorizer.get_feature_names()
    # get weights of words
    weights = np.array([tw_dict[word] for word in words])
    if re_weight > 0:
        weights = 1 + re_weight * weights

    for X in (X_train, X_test):
        if X is None:
            continue

        # sublinear_tf like tf-idf
        if sublinear_tf:
            X.data = np.log(X.data) + 1

        X = X.multiply(weights)  # can not use * to multiply
        if normalize:
            norm = sp.sparse.linalg.norm(X, axis=1)
            for i, row in enumerate(X.row):
                X.data[i] /= norm[row]

    return X_train, y_train, X_test

def main():
    print("data generating...")
    xy_url = from_project_root("processed_data/vector/{}_tf{}_{}gram_{}_XyN.pk"
                               .format(COLUMN, TW_TYPE, MAX_N, MAX_FEATURES))
    # test_url = None
    test_url = from_project_root('data/test_set.csv')
    if test_url:
        xy_url.replace('XyN', 'XyX_test')
    print("generated (X, y, X_test) will be saved at", xy_url)
    X, y, X_test = transform(TRAIN_URL, test_url, column=COLUMN, tw_type=TW_TYPE)
    joblib.dump((X, y, X_test), xy_url)
    pass

if __name__ == '__main__':
    main()

reason: import numpy,import scipy


text:
from copy import deepcopy
from scripts.fileReadWriteOperations import *
from scripts.removeRedundantTranscripts import *
import multiprocessing

import numpy as np

def extractSJ( eachinput ):
    """
    """
    bam_filename, output_filename, options = eachinput
    cmd = "regtools "
    cmd += " junctions "
    cmd += " extract "
    cmd += " -a 4 "
    cmd += " -m 20 "
    cmd += " -M 10000000000 "
    cmd += " -s 0 "
    cmd += " -o " + output_filename + ".temp "
    cmd += " > " + output_filename + ".output "
    cmd += " 2> " + output_filename + ".error "
    cmd += bam_filename
    if os.path.exists( bam_filename + ".bai" ) == False:
        cmd_samtools = "samtools index " + bam_filename
        os.system( cmd_samtools )
        cmd_samtools = "samtools index -c " + bam_filename
        os.system( cmd_samtools )
    if os.path.exists( output_filename + ".temp" ) == False:
        os.system( cmd )

    fhr = open( output_filename + ".temp", "r" )
    fhw = open( output_filename, "w" )
    for line in fhr:
        chromosome, chrom_start, chrom_end, name, read_support, strand, thick_start, thick_end, item_rgb, block_count, block_sizes, block_starts = line.strip().split( "\t" )
        block_sizes_int = list( map( int, block_sizes.split( "," ) ) )
        chrom_start, chrom_end = str( int( chrom_start ) + block_sizes_int[0] + 1 ), str( int( chrom_end ) - block_sizes_int[1] )
        fhw.write( "\t".join( [chromosome, chrom_start, chrom_end, name, read_support, strand, thick_start, thick_end, item_rgb, block_count, block_sizes, block_starts] ) + "\n" )
    fhw.close()
    fhr.close()

    cmd = "rm " + output_filename + ".temp"
    os.system( cmd )

def splitTranscriptsWithQuestionableSpliceJunctions( options ):
    """
    """
    pool = multiprocessing.Pool( processes = int( options.cpu ) )
    inputfilename = options.output_assemblies_psiclass_terminal_exon_length_modified + "/combined/combined_merged_transcripts.gtf"
    outputfilename = options.output_assemblies_psiclass_terminal_exon_length_modified + "/combined/combined_split_transcripts_with_bad_SJ.gtf"
    if os.path.exists( outputfilename ) == True:return
    transcript_info, useless1, useless2 = readAllTranscriptsFromGTFFileInParallel( [inputfilename, "dummy", "dummy"] )
    gene_to_transcripts = {}
    for chromosome in transcript_info:
        gene_to_transcripts[chromosome] = {}
        for transcript_id in transcript_info[chromosome]:
            gene_id = ".".join( transcript_id.split( "." )[:2] )
            if gene_id not in gene_to_transcripts[chromosome]:
                gene_to_transcripts[chromosome][gene_id] = []
            gene_to_transcripts[chromosome][gene_id].append( transcript_id )

    # Generate Splice junction information using regtools
    allinputs = []
    for condition in options.mrna_md:
        for Run in options.mrna_md[condition]:
            bam_filename = options.output_star + "/" + Run + "_for_psiclass.bam "
            output_filename = options.output_star + "/" + Run + "_SJ_regtools.bed"
            if os.path.exists( output_filename ) == False:
                allinputs.append( [bam_filename, output_filename, options] )
    pool.map( extractSJ, allinputs )

    all_SJ_info = {}
    allinputs = []
    for condition in options.mrna_md:
        all_SJ_info[condition] = {}
        for Run in options.mrna_md[condition]:
            all_SJ_info[condition][Run] = {}
            SJ_filename_regtools = options.output_star + "/" + Run + "_SJ_regtools.bed"
            allinputs.append( [SJ_filename_regtools, Run, condition, options] )
    results = pool.map( readFromRegtoolsOutput, allinputs )

    for result in results:
        SJ_info, Run, condition = result
        all_SJ_info[condition][Run] = SJ_info

    intron_info = []
    for chromosome in transcript_info:
        for transcript_id in transcript_info[chromosome]:
            # if transcript_id!="1.3.4":continue
            gene_id = ".".join( transcript_id.split( "." )[:2] )
            if len( gene_to_transcripts[chromosome][gene_id] ) == 1:continue
            if len( transcript_info[chromosome][transcript_id]["introns"] ) < 2:continue
            direction = transcript_info[chromosome][transcript_id]["direction"]
            present_in_num_conditions_direction_match = []
            present_in_num_conditions_direction_mismatch = []
            intron_coverages = []
            for condition in options.mrna_md:
                for Run in options.mrna_md[condition]:
                    intron_coverages_per_run = []
                    intron_sizes = []
                    for intron in transcript_info[chromosome][transcript_id]["introns"]:
                        intron_sizes.append( intron[1] - intron[0] + 1 )
                        if chromosome not in all_SJ_info[condition][Run]:continue
                        if str( intron[0] ) + "-" + str( intron[1] ) not in all_SJ_info[condition][Run][chromosome]:
                            intron_coverages_per_run.append( 0 )
                        else:
                            intron_coverages_per_run.append( all_SJ_info[condition][Run][chromosome][str( intron[0] ) + "-" + str( intron[1] )] )
                    if len( intron_coverages_per_run ) != 0:
                        intron_coverages.append( np.array( intron_coverages_per_run ) )
            """pprint.pprint(intron_coverages)
            sys.stdout.flush()"""
            if len( intron_coverages ) < 1:continue
            mean = list( np.mean( np.array( intron_coverages ), axis = 0 ) )
            std = []
            for ele in list( np.std( np.array( intron_coverages ), axis = 0 ) ):
                std.append( round( ele, 2 ) )
            if mean.index( min( mean ) ) != 0 and mean.index( min( mean ) ) != len( mean ) - 1:
                intron_info.append( [transcript_id, np.std( mean ), min( mean ), mean, intron_sizes] )
    all_stds = [row[1] for row in intron_info]
    all_min_mean = [row[2] for row in intron_info]
    ratios = [all_stds[i] / all_min_mean[i] for i in range( len( all_stds ) )]
    ratios = sorted( ratios )
    if len( ratios ) == 0:
        cmd = "cp " + inputfilename + " " + outputfilename
        os.system( cmd )
        return
    threshold = ratios[int( 0.75 * len( ratios ) )]
    mikado_results = {}
    """fhr=open("/work/LAS/rpwise-lab/sagnik/arabidopsis_gene_annotation/data/4_tissues_SE_PE/assemblies_psiclass_modified/combined/combined_merged_transcripts_mikado.tmap","r")
    for line in fhr:
        if line.strip().split("\t")[3] not in mikado_results:
            mikado_results[line.strip().split("\t")[3]]=[]
        mikado_results[line.strip().split("\t")[3]].append(line.strip().split("\t"))
    fhr.close()"""

    new_transcript_info = deepcopy( transcript_info )
    for row in intron_info:
        # print(row[1])
        transcript_id, std_of_intron_cov, min_mean, mean, intron_sizes = row
        # if std_of_intron_cov>=threshold and min_mean<10:
        ratio = std_of_intron_cov / min_mean
        if ratio > threshold:
            intron_indices_to_be_eliminated = [i for i in range( len( mean ) ) if mean[i] == min_mean]
            # intron_indices_to_be_eliminated.insert(0,0)
            # if len(intron_indices_to_be_eliminated)<=1:continue
            # intron_indices_to_be_eliminated.append(len(new_transcript_info[chromosome][transcript_id]["exons"])-1)

            """perfect_match="NO"
            for entry in mikado_results[transcript_id]:
                if "=" in entry[2]:
                    perfect_match="YES"
            print(" ".join([transcript_id,str(std_of_intron_cov),str(min_mean),perfect_match]),mean,intron_sizes,("large" if max(intron_sizes)>1000 else "small"))
            #print("Introns",intron_indices_to_be_eliminated,len(intron_indices_to_be_eliminated))
            """
            chromosome = transcript_id.split( "." )[0]
            previous_intron_list = new_transcript_info[chromosome][transcript_id]["introns"]
            new_exon_lists = []
            current_exon_chain = []

            for exon_num, exon in enumerate( new_transcript_info[chromosome][transcript_id]["exons"] ):
                if exon_num - 1 in intron_indices_to_be_eliminated:
                    new_exon_lists.append( current_exon_chain )
                    current_exon_chain = []
                    current_exon_chain.append( new_transcript_info[chromosome][transcript_id]["exons"][exon_num] )
                else:
                    current_exon_chain.append( new_transcript_info[chromosome][transcript_id]["exons"][exon_num] )
            new_exon_lists.append( current_exon_chain )
            """pprint.pprint(new_transcript_info[chromosome][transcript_id]["exons"])
            pprint.pprint(new_exon_lists)
            print("="*150)"""
            for exon_num, new_exon_list in enumerate( new_exon_lists ):
                new_transcript_id = transcript_id + "_isplit_" + str( exon_num )
                new_transcript_info[chromosome][new_transcript_id] = {"exons":new_exon_list,
                                                      "introns":[],
                                                      "cov":new_transcript_info[chromosome][transcript_id]["cov"],
                                                      "TPM":new_transcript_info[chromosome][transcript_id]["TPM"],
                                                      "FPKM":new_transcript_info[chromosome][transcript_id]["FPKM"],
                                                      "direction":new_transcript_info[chromosome][transcript_id]["direction"]
                                                      }
            del new_transcript_info[chromosome][transcript_id]

    writeTranscriptsToFile( [new_transcript_info, outputfilename] )
    removeRedundantTranscripts( outputfilename, outputfilename + "temp", options )
    cmd = "mv " + outputfilename + "temp " + outputfilename
    os.system( cmd )

    """cmd="mikado compare -r /work/LAS/rpwise-lab/sagnik/data/arath/transcriptome/Arabidopsis_thaliana.TAIR10.43.modified.gtf "
    cmd+=" -p "+options.output_assemblies_psiclass_terminal_exon_length_modified+"/combined/combined_split_transcripts_with_bad_SJ.gtf"
    cmd+=" -o "+options.output_assemblies_psiclass_terminal_exon_length_modified+"/combined/combined_split_transcripts_with_bad_SJ_mikado "
    os.system(cmd)"""


text:
# -*- coding: utf-8 -*-
import numpy as np

def get_axis(self):
    """Get the value of variables stored in Solution.

    Parameters
    ----------
    self : SolutionMat
        an SolutionMat object

    Returns
    -------
    axis_dict: dict
        a dict of axis names containing axis sizes

    """

    return self.axis


text:
# Author: Julien Mairal <julien.mairal@inria.fr>
#
# License: BSD 3 clause

from abc import abstractmethod
import numpy as np
import scipy.sparse
import cyanure_wrap

def preprocess(X,centering=False,normalize=True,columns=False):
    """ Perform in-place centering or normalization, either of columns or rows
    of the input matrix X

    Parameters
    ----------
    X : numpy array, or scipy sparse CSR matrix
        input matrix

    centering : boolean, default=False
        perform a centering operation

    normalize : boolean, default=True
        l2-normalization

    columns : boolean, default=False
        operates on rows (False) or columns (True).
    """

    if scipy.sparse.issparse(X):
        Xf=X.T
    else:
        Xf=np.asfortranarray(X.T)
    return cyanure_wrap.preprocess_(Xf,centering,normalize,not columns)

class ERM:
    """ The generic class for empirical risk minimization problems.
    For univariates problems, minimizes

        min_{w,b} (1/n) sum_{i=1}^n L( y_i, <w, x_i> + b)   + psi(w)

    """

    def __init__(self,loss='square',penalty='l2',fit_intercept=False):
        """ Initialization function of the ERM class.

        Parameters
        ----------
        loss : string, default='square'
            Loss function to be used. Possible choices are
            - 'square' =>  L(y,z) = 0.5 ( y-z)^2
            - 'logistic' => L(y,z) = log(1 + e^{-y z} )
            - 'sqhinge' or 'squared_hinge' => L(y,z) = 0.5 max( 0, 1- y z)^2
            - 'safe-logistic' => L(y,z) = e^{ yz - 1 } - y z  if yz <= 1 and 0 otherwise
            - 'multiclass-logistic' => multinomial logistic (see Latex documentation).
            Note that for binary classification, we assume the labels to be of the form {-1,+1}

        penalty: string, default='none'
            Regularization function psi. Possible choices are

            For univariate problems
            - 'none' => psi(w) = 0
            - 'l2' =>  psi{w) = (lambd/2) ||w||_2^2
            - 'l1' =>  psi{w) = lambd ||w||_1
            - 'elastic-net' =>  psi{w) = lambd ||w||_1 + (lambd2/2)||w||_2^2
            - 'fused-lasso' => psi(w) = lambd sum_{i=2}^p |w[i]-w[i-1]| + lambd2||w||_1 + (lambd3/2)||w||_2^2
            - 'l1-ball'     => encodes the constraint ||w||_1 <= lambd
            - 'l2-ball'     => encodes the constraint ||w||_2 <= lambd
            For multivariate problems, the previous penalties operate on each individual (e.g., class) predictor.
            In addition, multitask-group Lasso penalties are provided for multivariate problems (w is then a matrix)
            - 'l1l2' or 'l1linf', see Latex documentation

        fit_intercept: boolean, default='False'
            learns an unregularized intercept b  (or several intercepts for multivariate problems)
        """

        self.loss=loss
        if (loss=='squared_hinge'):
            self.loss='sqhinge'
        self.penalty=penalty
        self.fit_intercept=fit_intercept
        self.w=0
        self.b=0
        self.dual_variable=None;

    def fit(self,X,y,lambd=0,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_epochs=500,l_qning=20,f_restart=50,verbose=True,restart=False,univariate=True,nthreads=-1,seed=0):
        """
        The fitting function (the one that does the job)

        Parameters
        ----------

        X : numpy array, or scipy sparse CSR matrix
            input n x p numpy matrix; the samples are on the rows

        y : labels, numpy array.
            - vector of size n with real values for regression
            - vector of size n with {-1,+1} labels for binary classification, which will be automatically converted if labels in {0,1} are provided
            - matrix of size n x k for multivariate regression
            - vector of size n with entries in {0,1,k-1} for classification with k classes

        lambd: float, default=0
            first regularization parameter

        lambd2: float, default=0
            second regularization parameter, if needed

        lambd3: float, default=0
            third regularization parameter, if needed

        solver: string, default='auto'
            Optimization solver. Possible choices are
            - 'ista'
            - 'fista'
            - 'catalyst-ista'
            - 'qning-ista'  (proximal quasi-Newton method)
            - 'svrg'
            - 'catalyst-svrg' (accelerated SVRG with Catalyst)
            - 'qning-svrg'  (quasi-Newton SVRG)
            - 'acc-svrg'    (SVRG with direct acceleration)
            - 'miso'
            - 'catalyst-miso' (accelerated MISO with Catalyst)
            - 'qning-miso'  (quasi-Newton MISO)
            - 'auto'
            see the Latex documentation for more details. If you are unsure, use 'auto'

         tol: float, default='1e-3'
            Tolerance parameter. For almost all combinations of loss and
            penalty functions, this parameter is based on a duality gap.
            Assuming the (non-negative) objective function is "f" and its
            optimal value is "f^*", the algorithm stops with the guarantee

            f(x_t) - f^*  <=  tol f(x_t)

         max_epochs: int, default=500
            Maximum number of iteration of the algorithm in terms of passes over the data

         it0: int, default=10
            Frequency of duality-gap computation

         verbose: boolean, default=True
            Display information or not

         nthreads: int, default=-1
            maximum number of cores the method may use (-1 = all cores). Note that more cores
            is not always better.

         seed: int, default=0
            random seed

         restart: boolean, default=False
            use a restart strategy (useful for computing regularization path)

         univariate: boolean, default=True
            univariate or multivariate problems

         l_qning: int, default=20
            memory paramter for the qning method

         f_restart: int, default=50
            restart strategy for fista

        Returns
        -------

        test returns a numpy array carrying information about the optimization process
        (number of iterations, objective function values, duality gap)
        will be documented in the future if people ask me,
        """

        if y.ndim==0:
            print("Please provide label vector")
            return;
        Xf = X.T if scipy.sparse.issparse(X) else np.asfortranarray(X.T)
        p = X.shape[1]+1 if self.fit_intercept else X.shape[1]
        y=np.squeeze(y)
        if y.shape[0] != X.shape[0]:
            print("X and y should have the same number of observations")
            return;
        if univariate:
            w0 = np.zeros(p,dtype=Xf.dtype)
            yf=np.squeeze(y.astype(Xf.dtype))
        else:
            if y.squeeze().ndim > 1:
                nclasses=y.squeeze().shape[1]
                yf=np.asfortranarray(y.T)
            else:
                nclasses=int(np.max(y)+1)
                yf=np.squeeze(np.int32(y))
            print(nclasses)
            print(p)
            w0 = np.zeros([p,nclasses],dtype=Xf.dtype,order='F')

        if restart and np.any(self.w != 0):
            print("Restart")
            if self.fit_intercept:
                w0[-1,]=self.b
                w0[0:-1,]=self.w
            else:
                w0=self.w

        if restart and (solver=='auto' or solver=='miso' or solver=='catalyst-miso' or solver=='qning-miso'):
            n=X.shape[0];
            reset_dual = np.any(self.dual_variable==None)
            if not reset_dual and univariate:
                reset_dual = self.dual_variable.shape[0] != n
            if not reset_dual and not univariate:
                reset_dual = np.any(self.dual_variable.shape != [n,nclasses])
            if reset_dual and univariate:
                self.dual_variable=np.zeros(n,dtype=Xf.dtype)
            if reset_dual and not univariate:
                self.dual_variable=np.zeros([n,nclasses],dtype=Xf.dtype)

        w=np.copy(w0)
        optim_info=cyanure_wrap.erm_(Xf,yf,w0,w,dual_variable=self.dual_variable,loss=self.loss,penalty=self.penalty,solver=solver,lambd=float(lambd),lambd2=float(lambd2),lambd3=float(lambd3),intercept=bool(self.fit_intercept),tol=float(tol),it0=int(it0),nepochs=int(max_epochs),l_qning=int(l_qning),f_restart=int(f_restart),verbose=bool(verbose),univariate=bool(univariate),nthreads=int(nthreads),seed=int(seed))
        if self.fit_intercept:
            self.b=w[-1,]
            self.w=w[0:-1,]
        else:
            self.w=w
        return optim_info

    @abstractmethod
    def predict(self,X):
        """
        predict the labels given an input matrix X (same format as function fit)
        """
        pass

    def get_weights(self):
        """
        get the model parameters (either w or the tuple (w,b))
        """
        return (self.w,self.b) if self.fit_intercept else self.w

    def eval(self,X,y,lambd=0,lambd2=0,lambd3=0):
        """
        get the value of the objective function and computes a relative duality gap, see function fit for the format of parameters.
        """
        return self.fit(X,y,lambd=lambd,lambd2=lambd2,lambd3=lambd3,max_epochs=0,verbose=False,restart=True)

class BinaryClassifier(ERM):
    """ The binary classification class, which derives from ERM. The goal is to minimize
        the following objective

        .. math::
            \min_{w,b} \\frac{1}{n} \sum_{i=1}^n L\\left( y_i, w^\\top x_i + b\\right)   + \psi(w),

        where :math:`L` is a classification loss, :math:`\psi` is a regularization function (or constraint),
        :math:`w` is a p-dimensional vector representing model parameters,
        and b is an optional unregularized intercept. We expect binary labels in {-1,+1}.

        Parameters
        ----------
        loss: string, default='square'
            Loss function to be used. Possible choices are

            - 'square' =>  :math:`L(y,z) = \\frac{1}{2} ( y-z)^2`
            - 'logistic' => :math:`L(y,z) = \log(1 + e^{-y z} )`
            - 'sqhinge' or 'squared_hinge' => :math:`L(y,z) = \\frac{1}{2} \max( 0, 1- y z)^2`
            - 'safe-logistic' => :math:`L(y,z) = e^{ yz - 1 } - y z ~\\text{if}~ yz \leq 1~~\\text{and}~~0` otherwise

        penalty: string, default='l2'
            Regularization function psi. Possible choices are

            - 'none' => :math:`\psi(w) = 0`
            - 'l2' =>  :math:`\psi(w) = \\frac{\\lambda}{2} \|w\|_2^2`
            - 'l1' =>  :math:`\psi(w) = \\lambda \|w\|_1`
            - 'elastic-net' =>  :math:`\psi(w) = \\lambda \|w\|_1 + \\frac{\\lambda_2}{2}\|w\|_2^2`
            - 'fused-lasso' => :math:`\psi(w) = \\lambda \sum_{i=2}^p |w[i]-w[i-1]| + \\lambda_2\|w\|_1 + \\frac{\lambda_3}{2}\|w\|_2^2`
            - 'l1-ball'     => encodes the constraint :math:`\|w\|_1 \leq \lambda`
            - 'l2-ball'     => encodes the constraint :math:`\|w\|_2 \leq \lambda`

        fit_intercept: boolean, default='False'
            learns an unregularized intercept b

    """

    def fit(self,X,y,lambd=0,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_epochs=500,l_qning=20,f_restart=50,verbose=True,restart=False,nthreads=-1,seed=0):
        """
        The fitting function (the one that does the job)

        Parameters
        ----------

        X: numpy array, or scipy sparse CSR matrix
            input n x p numpy matrix; the samples are on the rows

        y: labels, numpy array.
            - vector of size n with {-1,+1} labels for binary classification, which will be automatically converted if labels in {0,1} are provided

        lambd: float, default=0
            first regularization parameter :math:`\lambda`

        lambd2: float, default=0
            second regularization parameter :math:`\lambda_2`, if needed

        lambd3: float, default=0
            third regularization parameter :math:`\lambda_3`, if needed

        solver: string, default='auto'
            Optimization solver. Possible choices are

            - 'ista'
            - 'fista'
            - 'catalyst-ista'
            - 'qning-ista'  (proximal quasi-Newton method)
            - 'svrg'
            - 'catalyst-svrg' (accelerated SVRG with Catalyst)
            - 'qning-svrg'  (quasi-Newton SVRG)
            - 'acc-svrg'    (SVRG with direct acceleration)
            - 'miso'
            - 'catalyst-miso' (accelerated MISO with Catalyst)
            - 'qning-miso'  (quasi-Newton MISO)
            - 'auto'
            see the Latex documentation for more details. If you are unsure, use 'auto'

         tol: float, default='1e-3'
            Tolerance parameter. For almost all combinations of loss and
            penalty functions, this parameter is based on a duality gap.
            Assuming the (non-negative) objective function is :math:`f` and its
            optimal value is :math:`f^*`, the algorithm stops with the guarantee

            .. math::
                f(x_t) - f^*  \leq  tol f(x_t)

         max_epochs: int, default=500
            Maximum number of iteration of the algorithm in terms of passes over the data

         it0: int, default=10
            Frequency of duality-gap computation

         verbose: boolean, default=True
            Display information or not

         nthreads: int, default=-1
            maximum number of cores the method may use (-1 = all cores). Note that more cores
            is not always better.

         seed: int, default=0
            random seed

         restart: boolean, default=False
            use a restart strategy (useful for computing regularization path)

         univariate: boolean, default=True
            univariate or multivariate problems

         l_qning: int, default=20
            memory parameter for the qning method

         f_restart: int, default=50
            restart strategy for fista

        Returns
        -------
        numpy array
            information about the optimization process (number of iterations,
            objective function values, duality gap) will be documented in the
            future if people ask me.
        """

        y=np.squeeze(y)
        uniq=np.unique(y)
        if uniq.shape[0] != 2:
            print(uniq.shape[0]+" classes detected; use MulticlassClassifier instead")
            return
        if not np.all(uniq == [-1, 1]):
            print("You have called BinaryClassifier, labels should be either -1 or 1, but they are")
            print(np.unique(y))
            print("Automatic conversion to [-1,1]")
            neg= y==uniq[0]
            y[neg]=-1
            y[np.logical_not(neg)]=1
        return super().fit(X,y,lambd=lambd,lambd2=lambd2,lambd3=lambd3,solver=solver,tol=tol,it0=it0,max_epochs=max_epochs,l_qning=l_qning,f_restart=f_restart,verbose=verbose,restart=restart,univariate=True,nthreads=nthreads,seed=seed)

    def predict(self,X):
        pred=X.dot(self.w)
        if self.fit_intercept: pred += self.b
        return np.sign(pred)

    def score(self,X,y):
        """ compute classification accuracy of the model for new test data (X,y)
        """
        y=np.squeeze(y)
        uniq=np.unique(y)
        if not np.all(uniq == [-1, 1]):
            print("You have called BinaryClassifier, labels should be either -1 or 1, but they are")
            print(np.unique(y))
            print("Automatic conversion to [-1,1]")
            neg= y==uniq[0]
            y[neg]=-1
            y[np.logical_not(neg)]=1
        pred=np.squeeze(self.predict(X))
        return np.sum(np.squeeze(y)==pred)/pred.shape[0]

class Regression(ERM):
    """ The regression class. The objective is the same as for the BinaryClassifier class,
        but we use a regression loss only (see below), and the targets will be real values.

        Parameters
        ----------
        loss: string, default='square'
            Only the square loss is implemented at this point

            - 'square' =>  :math:`L(y,z) = \\frac{1}{2} ( y-z)^2`

        penalty: string, default='l2'
            same as for the class BinaryClassifier

        fit_intercept: boolean, default='False'
            learns an unregularized intercept b

    """
    def __init__(self,loss='square',penalty='l2',fit_intercept=False):
        if self.loss!='square':
            print("square loss should be used")
            return
        super().__init__(self,loss=loss,penalty=penalty,fit_intercept=fit_intercept)

    def fit(self,X,y,lambd=0,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_epochs=500,l_qning=20,f_restart=50,verbose=True,restart=False,nthreads=-1,seed=0):
        """
        The fitting function is the same as for the class BinaryClassifier, except that we do not
        necessarily expect binary labels in y.
        """

        return super().fit(X,y,lambd=lambd,lambd2=lambd2,lambd3=lambd3,solver=solver,tol=tol,it0=it0,max_epochs=max_epochs,l_qning=l_qning,f_restart=f_restart,verbose=verbose,restart=restart,univariate=True,nthreads=nthreads,seed=seed)

    def predict(self,X):
        pred=X.dot(self.w)
        if self.fit_intercept: pred += self.b
        return pred

class MultiClassifier(ERM):
    """ The multi-class classification class. The goal is to minimize
        the following objective

        .. math::
            \min_{W,b} \\frac{1}{n} \sum_{i=1}^n L\\left( y_i, W^\\top x_i + b\\right)   + \psi(W),

        where :math:`L` is a classification loss, :math:`\psi` is a regularization function (or constraint),
        :math:`W=[w_1,\ldots,w_k]` is a (p x k) matrix that carries the k predictors, where k is the number
        of classes, and :math:`y_i` is a label in :math:`\{1,\ldots,k\}`.
        b is a k-dimensional vector representing an unregularized intercept (which is optional).

        Parameters
        ----------
        loss: string, default='square'
            Loss function to be used. Possible choices are

            - any loss function compatible with the class BinaryClassifier such as
              ('square', 'logistic', 'sqhinge', 'safe-logistic'). In such a
              case, the loss function encodes a one vs. all strategy based on
              the chosen binary-classification loss.
            - 'multiclass-logistic', which is also called multinomial or softmax logistic:

            .. math::
                L(y, W^\\top x + b) = \sum_{j=1}^k \log\left(e^{w_j^\\top + b_j} - e^{w_y^\\top + b_y} \\right)

        penalty: string, default='l2'
            Regularization function psi. Possible choices are

            - any penalty function compatible with the class BinaryClassifier such as
                ('none', 'l2', 'l1', 'elastic-net', 'fused-lasso', 'l1-ball', 'l2-ball'). In such a case,. the penalty is
                applied on each predictor :math:`w_j` individually:

            .. math::
                \psi(W) = \sum_{j=1}^k \psi(w_j).
            - 'l1l2', which is the multi-task group Lasso regularization

            .. math::
                \psi(W) = \sum_{j=1}^p \|W^j\|_2~~~~\\text{where}~W^j~\\text{is the j-th row of}~W.
            - 'l1linf'

            .. math::
                \psi(W) = \sum_{j=1}^p \|W^j\|_\infty.

        fit_intercept: boolean, default='False'
            learns an unregularized intercept b, which is a k-dimensional vector
    """
    def fit(self,X,y,lambd=0,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_epochs=500,l_qning=20,f_restart=50,verbose=True,restart=False,nthreads=-1,seed=0):
        """ same as BinaryClassifier, but y should be a vector a n-dimensional vector of integers
        """
        y=np.squeeze(y)
        if y.squeeze().ndim != 1 or np.any(y != y.astype(int)):
            print("y should be a n-dimensional vector of integers")
            return
        nclasses=np.max(y)+1
        uniqu=np.unique(y)
        if nclasses==2:
            print("Two classes detected, use BinaryClassifier instead")
            return
        if nclasses!= uniqu.shape[0] or not all(np.unique(y) == np.arange(nclasses)):
            print("Class labels should be of the form")
            print(np.arange(nclasses))
            print("but they are")
            print(uniqu)
            return
        return super().fit(X,y,lambd=lambd,lambd2=lambd2,lambd3=lambd3,solver=solver,tol=tol,it0=it0,max_epochs=max_epochs,l_qning=l_qning,f_restart=f_restart,verbose=verbose,restart=restart,univariate=False,nthreads=nthreads,seed=seed)

    def predict(self,X):
        """ predicts the class label
        """
        pred=X.dot(self.w)
        if self.fit_intercept: pred += self.b[newaxix,:]
        return np.argmax(pred,1)

    def score(self,X,y):
        """ gives a classification score on new test data
        """
        pred=np.squeeze(self.predict(X))
        return np.sum(np.squeeze(y)==pred)/pred.shape[0]

class MultiVariateRegression(ERM):
    """ The multivariate regression class. The objective is the same as for the MultiClassifier class,
        but we use a regression loss only (see below), and the targets :math:`y_i` are k-dimensional vectors.

        Parameters
        ----------
        loss: string, default='square'
            Only the square loss is implemented at this point. Given two k-dimensional vectors y,z:

            - 'square' =>  :math:`L(y,z) = \\frac{1}{2} \|y-z\|^2`

        penalty: string, default='l2'
            same as for the class MultiClassifier

        fit_intercept: boolean, default='False'
            learns an unregularized intercept b
    """
    def __init__(self,loss='square',penalty='l2',fit_intercept=False):
        if self.loss!='square':
            print("square loss should be used")
            return
        super().__init__(self,loss=loss,penalty=penalty,fit_intercept=fit_intercept)

    def fit(self,X,y,lambd=0,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_epochs=500,l_qning=20,f_restart=50,verbose=True,restart=False,nthreads=-1,seed=0):
        """ same as ERM.fit, but y should be n x k, where k is size of the target for each data point
        """
        if y.squeeze().ndim <= 1:
            print("y should be n x k, where k is size of the target for each data point")
            return
        return super().fit(X,y,lambd=lambd,lambd2=lambd2,lambd3=lambd3,solver=solver,tol=tol,it0=it0,max_epochs=max_epochs,l_qning=l_qning,f_restart=f_restart,verbose=verbose,restart=restart,univariate=False,nthreads=nthreads,seed=seed)

    def predict(self,X):
        """ predicts the targets
        """
        pred=X.dot(self.w)
        if self.fit_intercept: pred += self.b[newaxix,:]
        return pred

class LinearSVC(BinaryClassifier):
    """ A compatibility class for scikit-learn user, but only for square hinge loss
        It is perfectly equivalent to the BinaryClassifier class, but the regularization
        parameter (here "C") is provided during the class initialization.
        Note that :math:`C= \\frac{1}{2n \lambda}`

        Parameters
        ----------

        loss: should be 'sqhinge' or 'squared_hinge'

        penalty: same as BinaryClassifier

        fit_intercept: same as BinaryClassifier

        C: regularization parameter

        max_iter: maximum number of iterations for the optimization solver
    """
    def __init__(self,loss='sqhinge',penalty='l2',fit_intercept=False,C=1,max_iter=500):
        if loss != 'sqhinge' and loss != 'squared_hinge':
            print("LinearSVC is only compatible with squared hinge loss at the moment")
        super().__init__(self,loss='sqhinge',penalty=penalty,fit_intercept=fit_intercept,verbose=verbose,max_epochs=max_iter)
        self.C=C
        self.verbose=verbose
        self.max_iter=max_iter

    def fit(self,X,y,C=None,verbose=None,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_iter=None,l_qning=20,f_restart=50,restart=False,nthreads=-1,seed=0):
        """ same as BinaryClassification, but the parameter C replaces lambd, and max_iter replaces max_epochs """
        n=X.shape[0]
        if C != None:
            self.C=C
        if verbose != None:
            self.verbose=verbose
        if max_iter != None:
            self.max_iter=max_iter
        super().fit(self,X=X,y=y,lambd=1/(2*n*self.C),lambd2=lambd2,lambd3=lambd3,solver=solver,tol=tol,it0=it0,max_epochs=max_epochs,l_qning=l_qning,f_restart=f_rstart,verbose=verbose,restart=restart,nthreads=nthreads,seed=seed)

class LogisticRegression(BinaryClassifier):
    """ A compatibility class for scikit-learn user, but only for square hinge loss
        It is perfectly equivalent to the BinaryClassifier class, but the regularization
        parameter (here "C") is provided during the class initialization.
        Note that :math:`C= \\frac{1}{n \lambda}`

        Parameters
        ----------

        loss: should be 'sqhinge' or 'squared_hinge'

        penalty: same as BinaryClassifier

        fit_intercept: same as BinaryClassifier

        C: regularization parameter

        max_iter: maximum number of iterations for the optimization solver
    """

    def __init__(self,penalty='l2',fit_intercept=False,C=1,max_iter=500):
        super().__init__(self,loss='logistic',penalty=penalty,fit_intercept=fit_intercept,verbose=verbose,max_epochs=max_iter)
        self.C=C
        self.verbose=verbose
        self.max_iter=max_iter

    def fit(self,X,y,C=None,lambd2=0,lambd3=0,solver='auto',tol=1e-3,it0=10,max_iter=None,l_qning=20,f_restart=50,verbose=None,restart=False,nthreads=-1,seed=0):
        """ same as BinaryClassification, but the parameter C replaces lambd, and max_iter replaces max_epochs """
        n=X.shape[0]
        if C != None:
            self.C=C
        if verbose != None:
            self.verbose=verbose
        if max_iter != None:
            self.max_iter=max_iter
        super().fit(self,X=X,y=y,lambd=1/(n*self.C),lambd2=lambd2,lambd3=lambd3,solver=solver,tol=tol,it0=it0,max_epochs=max_epochs,l_qning=l_qning,f_restart=f_rstart,verbose=verbose,restart=restart,nthreads=nthreads,seed=seed)

max_stars_repo_name: tomMoral/cyanure

max_issues_repo_name: tomMoral/cyanure

max_forks_repo_name: tomMoral/cyanure