text:
# -*- coding: utf-8 -*-
#
#  This file is part of Sequana software
#
#  Copyright (c) 2016-2020 - Sequana Development Team
#
#  File author(s):
#      Thomas Cokelaer <thomas.cokelaer@pasteur.fr>
#      Etienne Kornobis <etienne.kornobis@pasteur.fr>
#
#  Distributed under the terms of the 3-clause BSD license.
#  The full license is in the LICENSE file, distributed with this software.
#
#  website: https://github.com/sequana/sequana
#  documentation: http://sequana.readthedocs.io
#
##############################################################################
import sys
import os

from pathlib import Path
from jinja2 import Environment, PackageLoader
import subprocess
import seaborn as sns
from itertools import combinations
import os

from sequana.lazy import pandas as pd
from sequana.lazy import pylab
from sequana.lazy import numpy as np

from sequana.gff3 import GFF3
from sequana.viz import Volcano
from sequana.enrichment import PantherEnrichment, KeggPathwayEnrichment

import colorlog

logger = colorlog.getLogger(__name__)

import glob

__all__ = ["RNADiffAnalysis", "RNADiffResults", "RNADiffTable", "RNADesign"]

def strip(text):
    try:
        return text.strip()
    except AttributeError:
        return text

class RNADesign:
    """Simple RNA design handler"""

    def __init__(self, filename, sep="\s*,\s*", reference=None):
        self.filename = filename
        # \s to strip the white spaces
        self.df = pd.read_csv(filename, sep=sep, engine="python",
                comment="#",
                dtype={"label": str})

        self.reference = reference

    def _get_conditions(self):
        return sorted(self.df.condition.unique())

    conditions = property(_get_conditions)

    def _get_comparisons(self):
        conditions = self.conditions
        if self.reference is None:
            import itertools

            comps = list(itertools.combinations(conditions, 2))
        else:
            # only those versus reference
            comps = [(x, self.reference) for x in conditions if x != self.reference]
        return sorted(comps)

    comparisons = property(_get_comparisons)

class RNADiffAnalysis:
    """A tool to prepare and run a RNA-seq differential analysis with DESeq2

    :param counts_file: Path to tsv file out of FeatureCount with all samples together.
    :param design_file: Path to tsv file with the definition of the groups for each sample.
    :param condition: The name of the column from groups_tsv to use as condition. For more
        advanced design, a R function of the type 'condition*inter' (without the '~') could
        be specified (not tested yet). Each name in this function should refer to column
        names in groups_tsv.
    :param comparisons: A list of tuples indicating comparisons to be made e.g A vs B would be [("A", "B")]
    :param batch: None for no batch effect or name of a column in groups_tsv to add a batch effec.
    :param fit_type: Default "parametric".
    :param beta_prior: Default False.
    :param independent_filtering: To let DESeq2 perform the independentFiltering or not.
    :param cooks_cutoff: To let DESeq2 decide for the CooksCutoff or specifying a value.
    :param gff: Path to the corresponding gff3 to add annotations.
    :param fc_attribute: GFF attribute used in FeatureCounts.
    :param fc_feature: GFF feaure used in FeatureCounts.
    :param annot_cols: GFF attributes to use for results annotations
    :param threads: Number of threads to use
    :param outdir: Path to output directory.
    :param sep_counts: The separator used in the input count file.
    :param sep_design: The separator used in the input design file.

    This class reads a :class:`sequana.featurecounts.`

    ::

        r = rnadiff.RNADiffAnalysis("counts.csv", "design.csv",
                condition="condition", comparisons=[(("A", "B"), ('A', "C")],

    """

    _template_file = "rnadiff_light_template.R"
    _template_env = Environment(loader=PackageLoader("sequana", "resources/templates"))
    template = _template_env.get_template(_template_file)

    def __init__(
        self,
        counts_file,
        design_file,
        condition,
        comparisons,
        batch=None,
        fit_type="parametric",
        beta_prior=False,
        independent_filtering=True,
        cooks_cutoff=None,
        gff=None,
        fc_attribute=None,
        fc_feature=None,
        annot_cols=None,
        # annot_cols=["ID", "Name", "gene_biotype"],
        threads=4,
        outdir="rnadiff",
        sep_counts=",",
        sep_design="\s*,\s*",
    ):

        self.outdir = Path(outdir)

        self.usr_counts = counts_file
        self.usr_design = design_file

        self.counts_filename = self.outdir / "counts.csv"
        self.design_filename = self.outdir / "design.csv"

        self.counts, self.design = self.check_and_save_input_tables(
            sep_counts, sep_design
        )

        self.condition = condition
        self.check_condition()

        self.comparisons = comparisons
        self.check_comparisons()

        self.comparisons_str = (
            f"list({', '.join(['c' + str(x) for x in self.comparisons])})"
        )
        self.batch = batch
        self.model = f"~{batch + '+' + condition if batch else condition}"
        self.fit_type = fit_type
        self.beta_prior = "TRUE" if beta_prior else "FALSE"
        self.independent_filtering = "TRUE" if independent_filtering else "FALSE"
        self.cooks_cutoff = cooks_cutoff if cooks_cutoff else "TRUE"
        self.gff = gff
        self.fc_feature = fc_feature
        self.fc_attribute = fc_attribute
        self.annot_cols = annot_cols
        self.threads = threads

    def __repr__(self):
        info = f"RNADiffAnalysis object:\n\
- {self.counts.shape[1]} samples.\n\
- {len(self.comparisons)} comparisons.\n\n\
Counts overview:\n\
{self.counts.head()}\n\n\
Design overview:\n\
{self.design.head()}"

        return info

    def check_and_save_input_tables(self, sep_counts, sep_design):
        try:
            self.outdir.mkdir()
        except:
            pass

        counts = pd.read_csv(
            self.usr_counts, sep=sep_counts, index_col="Geneid", comment="#"
        ).sort_index(axis=1)

        design = RNADesign(self.usr_design, sep=sep_design)
        design = design.df.set_index("label").sort_index()

        if list(counts.columns) != list(design.index):
            logger.error(
                f"Counts columns and design rows does not match (after sorting)."
            )
            logger.error(counts.columns)
            logger.error(design.index)
            sys.exit(1)

        counts.to_csv(self.counts_filename)
        design.to_csv(self.design_filename)

        return counts, design

    def check_condition(self):
        # set condition of the statistical model

        columns = ",".join(self.design.columns)
        if self.condition not in columns:
            logger.error(
                f"""Your condition named '{condition}' is expected to
be in the header of your design file but was not found. Candidates are:
    {columns}"""
            )
            sys.exit(1)

    def check_comparisons(self):
        # let us check the consistenct of the design and comparisons
        valid_conditions = ",".join(set(self.design[self.condition].values))
        for item in [x for y in self.comparisons for x in y]:
            if item not in self.design[self.condition].values:
                logger.error(
                    f"""{item} not found in the design. Fix the design
or comparisons. possible values are {valid_conditions}"""
                )
                sys.exit(1)

    def run(self):
        """Create outdir and a DESeq2 script from template for analysis. Then execute
        this script.

        :return: a :class:`RNADiffResults` instance
        """
        logger.info("Running DESeq2 analysis. Please wait")

        rnadiff_script = self.outdir / "rnadiff_light.R"

        with open(rnadiff_script, "w") as f:
            f.write(RNADiffAnalysis.template.render(self.__dict__))

        logger.info("Starting differential analysis with DESeq2...")

        # capture_output is valid for py3.7 and above so we will use
        # stdout/stderr to be back compatible with py3.6
        # Capture rnadiff output
        p = subprocess.Popen(
            f"Rscript {rnadiff_script}",
            shell=True,
            universal_newlines=True,
            stderr=subprocess.PIPE,
            stdout=subprocess.PIPE,
        )
        p.wait()
        stdout, stderr = p.stdout.read(), p.stderr.read()

        # Capture rnadiff output, Unfortunately, R code mixes stdout/stderr
        # FIXME
        with open(self.outdir / "rnadiff.err", "w") as f:
            f.write(stderr)
        with open(self.outdir / "rnadiff.out", "w") as f:
            f.write(stdout)

        # if os.path.exists():
        #    logger.error(f"stderr of R call is not empty: {p.stderr}")
        #    sys.exit(1)

        results = RNADiffResults(
            self.outdir,
            self.design_filename,
            condition=self.condition,
            gff=self.gff,
            fc_feature=self.fc_feature,
            fc_attribute=self.fc_attribute,
            annot_cols=self.annot_cols,
        )
        return results

class RNADiffTable:
    def __init__(self, path, alpha=0.05, log2_fc=0, sep=",", condition="condition"):
        """ A representation of the results of a single rnadiff comparison """
        self.path = Path(path)
        self.name = self.path.stem.replace("_degs_DESeq2", "").replace("-", "_")

        self._alpha = alpha
        self._log2_fc = log2_fc

        self.df = pd.read_csv(self.path, index_col=0, sep=sep)
        self.condition = condition

        self.filt_df = self.filter()
        self.set_gene_lists()

    @property
    def alpha(self):
        return self._alpha

    @alpha.setter
    def alpha(self, value):
        self._alpha = value
        self.filt_df = self.filter()
        self.set_gene_lists()

    @property
    def log2_fc(self):
        return self._log2_fc

    @log2_fc.setter
    def log2_fc(self, value):
        self._log2_fc = value
        self.filt_df = self.filter()
        self.set_gene_lists()

    def filter(self):
        """filter a DESeq2 result with FDR and logFC thresholds"""

        fc_filt = self.df["log2FoldChange"].abs() < self._log2_fc
        fdr_filt = self.df["padj"] > self._alpha
        outliers = self.df["padj"].isna()

        filt_df = self.df.copy()
        filt_df[fc_filt.values | fdr_filt.values | outliers] = np.NaN
        return filt_df

    def set_gene_lists(self):

        only_drgs_df = self.filt_df.dropna(how="all")

        self.gene_lists = {
            "up": list(only_drgs_df.query("log2FoldChange > 0").index),
            "down": list(only_drgs_df.query("log2FoldChange < 0").index),
            "all": list(only_drgs_df.index),
        }

    def summary(self):

        return pd.DataFrame(
            {
                "log2_fc": self._log2_fc,
                "alpha": self._alpha,
                "up": len(self.gene_lists["up"]),
                "down": len(self.gene_lists["down"]),
                "all": len(self.gene_lists["all"]),
            },
            index=[self.name],
        )

    def plot_volcano(
        self,
        padj=0.05,
        add_broken_axes=False,
        markersize=4,
        limit_broken_line=[20, 40],
        plotly=False,
        annotations=None,
    ):
        """

        .. plot::
            :include-source:

            from sequana.rnadiff import RNADiffResults
            from sequana import sequana_data

            r = RNADiffResults(sequana_data("rnadiff/rnadiff_onecond_1"))
            r.plot_volcano()

        """

        if plotly:
            from plotly import express as px

            df = self.df.copy()

            if annotations is not None:
                try:
                    df = pd.concat([df, annotations.annotation], axis=1)
                except Exception as err:
                    logger.warning(
                        f"Could not merge rnadiff table with annotation. Full error is: {err}"
                    )
            df["log_adj_pvalue"] = -pylab.log10(df.padj)
            df["significance"] = [
                "<{}".format(padj) if x else ">={}".format(padj) for x in df.padj < padj
            ]

            if "Name" in df.columns:
                hover_name = "Name"
            elif "gene_id" in df.columns:
                hover_name = "gene_id"
            elif "locus_tag" in df.columns:
                hover_name = "locus_tag"
            elif "ID" in df.columns:
                hover_name = "ID"
            else:
                hover_name = None
            fig = px.scatter(
                df,
                x="log2FoldChange",
                y="log_adj_pvalue",
                hover_name=hover_name,
                hover_data=["baseMean"],
                log_y=False,
                opacity=0.5,
                color="significance",
                height=600,
                labels={"log_adj_pvalue": "log adjusted p-value"},
            )
            # axes[0].axhline(
            # -np.log10(0.05), lw=2, ls="--", color="r", label="pvalue threshold (0.05)"
            # i)
            # in future version of plotly, a add_hlines will be available. For
            # now, this is the only way to add axhline
            fig.update_layout(
                shapes=[
                    dict(
                        type="line",
                        xref="x",
                        x0=df.log2FoldChange.min(),
                        x1=df.log2FoldChange.max(),
                        yref="y",
                        y0=-pylab.log10(padj),
                        y1=-pylab.log10(padj),
                        line=dict(color="black", width=1, dash="dash"),
                    )
                ]
            )

            return fig

        from brokenaxes import brokenaxes

        M = max(-pylab.log10(self.df.padj.dropna()))

        br1, br2 = limit_broken_line
        if M > br1:
            if add_broken_axes:
                bax = brokenaxes(ylims=((0, br1), (M - 10, M)), xlims=None)
            else:
                bax = pylab
        else:
            bax = pylab

        d1 = self.df.query("padj>@padj")
        d2 = self.df.query("padj<=@padj")
        bax.plot(
            d1.log2FoldChange,
            -np.log10(d1.padj),
            marker="o",
            alpha=0.5,
            color="k",
            lw=0,
            markersize=markersize,
        )
        bax.plot(
            d2.log2FoldChange,
            -np.log10(d2.padj),
            marker="o",
            alpha=0.5,
            color="r",
            lw=0,
            markersize=markersize,
        )

        bax.grid(True)
        try:
            bax.set_xlabel("fold change")
            bax.set_ylabel("log10 adjusted p-value")
        except:
            bax.xlabel("fold change")
            bax.ylabel("log10 adjusted p-value")

        m1 = abs(min(self.df.log2FoldChange))
        m2 = max(self.df.log2FoldChange)
        limit = max(m1, m2)
        try:
            bax.set_xlim([-limit, limit])
        except:
            bax.xlim([-limit, limit])
        try:
            y1, _ = bax.get_ylim()
            ax1 = bax.axs[0].set_ylim([br2, y1[1] * 1.1])
        except:
            y1, y2 = bax.ylim()
            bax.ylim([0, y2])
        bax.axhline(
            -np.log10(0.05), lw=2, ls="--", color="r", label="pvalue threshold (0.05)"
        )
        return bax

        if colors is None:
            colors = {}
            for sample in self.sample_names:
                colors[sample] = self.colors[self.get_cond_from_sample(sample)]

        if plotly is True:
            assert n_components == 3
            variance = p.plot(
                n_components=n_components,
                colors=colors,
                show_plot=False,
                max_features=max_features,
            )
            from plotly import express as px

            df = pd.DataFrame(p.Xr)
            df.columns = ["PC1", "PC2", "PC3"]
            df["names"] = self.sample_names
            df["colors"] = [colors[x] for x in self.sample_names]
            df["size"] = [10] * len(df)
            df[self.condition] = [
                self.get_cond_from_sample(sample) for sample in self.sample_names
            ]
            fig = px.scatter_3d(
                df,
                x="PC1",
                y="PC2",
                z="PC3",
                color=self.condition,
                labels={
                    "PC1": "PC1 ({}%)".format(round(100 * variance[0], 2)),
                    "PC2": "PC2 ({}%)".format(round(100 * variance[1], 2)),
                    "PC3": "PC3 ({}%)".format(round(100 * variance[2], 2)),
                },
                height=800,
                text="names",
            )
            return fig
        else:
            variance = p.plot(
                n_components=n_components, colors=colors, max_features=max_features
            )

        return variance

    def plot_pvalue_hist(self, bins=60, fontsize=16, rotation=0):

        pylab.hist(self.df.pvalue.dropna(), bins=bins, ec="k")
        pylab.xlabel("raw p-value")
        pylab.ylabel("Occurences")

    def plot_pvalue_hist(self, bins=60, fontsize=16, rotation=0):
        pylab.hist(self.df.pvalue.dropna(), bins=bins, ec="k")
        pylab.grid(True)
        pylab.xlabel("raw p-value", fontsize=fontsize)
        pylab.ylabel("Occurences", fontsize=fontsize)
        try:
            pylab.tight_layout()
        except:
            pass

    def plot_padj_hist(self, bins=60, fontsize=16):
        pylab.hist(self.df.padj.dropna(), bins=bins, ec="k")
        pylab.grid(True)
        pylab.xlabel("Adjusted p-value", fontsize=fontsize)
        pylab.ylabel("Occurences", fontsize=fontsize)
        try:
            pylab.tight_layout()
        except:
            pass

class RNADiffResults:
    """The output of a RNADiff analysis"""

    def __init__(
        self,
        rnadiff_folder,
        design_file=None,
        gff=None,
        fc_attribute=None,
        fc_feature=None,
        pattern="*vs*_degs_DESeq2.csv",
        alpha=0.05,
        log2_fc=0,
        palette=sns.color_palette(desat=0.6),
        condition="condition",
        annot_cols=None,
        # annot_cols=["ID", "Name", "gene_biotype"],
    ):
        """

        :rnadiff_folder:

        """
        self.path = Path(rnadiff_folder)
        self.files = [x for x in self.path.glob(pattern)]

        self.counts_raw = pd.read_csv(
            self.path / "counts_raw.csv", index_col=0, sep=","
        )
        self.counts_raw.sort_index(axis=1, inplace=True)

        self.counts_norm = pd.read_csv(
            self.path / "counts_normed.csv", index_col=0, sep=","
        )
        self.counts_norm.sort_index(axis=1, inplace=True)

        self.counts_vst = pd.read_csv(
            self.path / "counts_vst_norm.csv", index_col=0, sep=","
        )
        self.counts_vst.sort_index(axis=1, inplace=True)

        self.dds_stats = pd.read_csv(
            self.path / "overall_dds.csv", index_col=0, sep=","
        )
        self.condition = condition
        self.design_df = self._get_design(
            design_file, condition=self.condition, palette=palette
        )

        # FIXME: This block was commented because:
        # Now RNADiffAnalysis create a sorted design file according to the count matrix
        # Not sure that creating a design from scratch at this stage is relevant.

        # read different results and sort by sample name all inputs
        # TODO make this a function to be reused in RNADiffAnalysis for example
        # if design_file == None:
        #     conditions = []
        #     labels = self.counts_raw.columns
        #     for label in labels:
        #         condition = input(
        #             f"Please give use a condition name for the {label} label: "
        #         )
        #         conditions.append(condition)
        #     df = pd.DataFrame({"label": labels, "condition": conditions})
        #     df.set_index("label", inplace=True)
        #     df.sort_index(inplace=True)
        #     col_map = dict(zip(df.loc[:, condition].unique(), palette))
        #     df["group_color"] = df.loc[:, condition].map(col_map)
        #     self.design_df = df
        # else:
        #     self.design_df = self._get_design(
        #         design_file, condition=condition, palette=palette
        #     )
        #     self.design_df.sort_index(inplace=True)
        #     self.design = RNADesign(design_file)

        # optional annotation
        self.fc_attribute = fc_attribute
        self.fc_feature = fc_feature
        self.annot_cols = annot_cols
        if gff:
            if fc_feature is None or fc_attribute is None:
                logger.error(
                    "Since you provided a GFF filem you must provide the feature and attribute to be used."
                )
            self.annotation = self.read_annot(gff)
        else:
            self.annotation = None

        # some filtering attributes
        self._alpha = alpha
        self._log2_fc = log2_fc

        self.comparisons = self.import_tables()

        self.df = self._get_total_df()
        self.filt_df = self._get_total_df(filtered=True)

    @property
    def alpha(self):
        return self._alpha

    @alpha.setter
    def alpha(self, value):
        self._alpha = value
        self.comparisons = self.import_tables()
        self.filt_df = self._get_total_df(filtered=True)

    @property
    def log2_fc(self):
        return self._log2_fc

    @log2_fc.setter
    def log2_fc(self, value):
        self._log2_fc = value
        self.comparisons = self.import_tables()
        self.filt_df = self._get_total_df(filtered=True)

    def to_csv(self, filename):
        self.df.to_csv(filename)

    def read_csv(self, filename):
        self.df = pd.read_csv(filename, index_col=0, header=[0, 1])

    def import_tables(self):
        from easydev import AttrDict

        data = {
            compa.stem.replace("_degs_DESeq2", "").replace("-", "_"): RNADiffTable(
                compa,
                alpha=self._alpha,
                log2_fc=self._log2_fc,
                condition=self.condition,
                # gff=self.annotation.annotation,
            )
            for compa in self.files
        }

        return AttrDict(**data)

    def read_annot(self, gff_filename):
        """Get a properly formatted dataframe from the gff."""

        gff = GFF3(gff_filename)
        df = gff.df

        if self.annot_cols is None:
            lol = [
                list(x.keys())
                for x in df.query("type==@self.fc_feature")["attributes"].values
            ]
            annot_cols = list(set([x for item in lol for x in item]))
        else:
            annot_cols = self.annot_cols

        df = df.query("type == @self.fc_feature").loc[:, annot_cols]
        df.drop_duplicates(inplace=True)

        df.set_index(self.fc_attribute, inplace=True)

        # It may happen that a GFF has duplicated IDs ! For instance ecoli
        # has 20 duplicated ID that are part 1 and 2 of the same gene
        df = df[~df.index.duplicated(keep="last")]

        df.columns = pd.MultiIndex.from_product([["annotation"], df.columns])

        return df

    def _get_total_df(self, filtered=False):
        """Concatenate all rnadiff results in a single dataframe.

        FIXME: Columns relative to significative comparisons are not using
        self.log2_fc and self.alpha
        """

        dfs = []

        for compa, res in self.comparisons.items():
            df = res.filt_df if filtered else res.df
            df = df.transpose().reset_index()
            df["file"] = res.name
            df = df.set_index(["file", "index"])
            dfs.append(df)

        df = pd.concat(dfs, sort=True).transpose()

        # Add number of comparisons which are significative for a given gene
        num_sign_compa = (df.loc[:, (slice(None), "padj")] < 0.05).sum(axis=1)
        df.loc[:, ("statistics", "num_of_significative_comparisons")] = num_sign_compa

        # Add list of comparisons which are significative for a given gene
        df_sign_padj = df.loc[:, (slice(None), "padj")] < 0.05
        sign_compa = df_sign_padj.loc[:, (slice(None), "padj")].apply(
            # Extract column names (comparison names) for significative comparisons
            lambda row: {col_name[0] for sign, col_name in zip(row, row.index) if sign},
            axis=1,
        )
        df.loc[:, ("statistics", "significative_comparisons")] = sign_compa

        if self.annotation is not None and self.fc_attribute and self.fc_feature:
            df = pd.concat([self.annotation, df], axis=1)
        else:
            logger.warning(
                "Missing any of gff, fc_attribute or fc_feature. No annotation will be added."
            )

        return df

    def summary(self):
        return pd.concat(res.summary() for compa, res in self.comparisons.items())

    def report(self):

        template_file = "rnadiff_report.html"
        template_env = Environment(
            loader=PackageLoader("sequana", "resources/templates")
        )
        template = template_env.get_template(template_file)

        with open("rnadiff_report.html", "w") as f:
            f.write(
                template.render(
                    {"table": self.summary().to_html(classes="table table-striped")}
                )
            )

    def run_enrichment_go(
        self, taxon, annot_col="Name", out_dir="enrichment"
    ):  # pragma: no cover

        out_dir = Path(out_dir) / "figures"
        out_dir.mkdir(exist_ok=True, parents=True)

        gene_lists_dict = self.get_gene_lists(
            annot_col=annot_col, Nmax=2000, dropna=True
        )
        enrichment = {}
        ontologies = {"GO:0003674": "BP", "GO:0008150": "MF", "GO:0005575": "CC"}
        failed_enrichments = []

        for compa in self.comparisons:
            gene_lists = gene_lists_dict[compa]
            pe = PantherEnrichment(gene_lists, taxon)
            pe.compute_enrichment(ontologies=ontologies.keys(), progress=False)

            for direction in ["up", "down", "all"]:
                if not pe.enrichment[direction]:
                    logger.warning(
                        f"No enrichment computed, so no plots computed for {compa} {direction} {ontology}"
                    )
                    failed_enrichments.append(
                        {
                            "comparison": compa,
                            "direction": direction,
                            "GO": "all",
                            "reason": "no enrichment computed",
                        }
                    )
                    continue

                for ontology in ontologies.keys():
                    pylab.figure()
                    enrichment_df = pe.plot_go_terms(
                        direction, ontology, compute_levels=False
                    )
                    if enrichment_df.empty:
                        failed_enrichments.append(
                            {
                                "comparison": compa,
                                "direction": direction,
                                "GO": ontology,
                                "reason": "no enrichment found",
                            }
                        )
                    else:
                        enrichment[(compa, direction, ontology)] = enrichment_df
                        pylab.tight_layout()
                        pylab.savefig(
                            out_dir
                            / f"go_{compa}_{direction}_{ontologies[ontology]}.pdf"
                        )
                        pe.save_chart(
                            enrichment_df,
                            out_dir
                            / f"chart_{compa}_{direction}_{ontologies[ontology]}.png",
                        )

            logger.info(f"Panther enrichment for {compa} DONE.")

        df = pd.concat(enrichment).sort_index()
        df.index.rename(
            ["comparison", "direction", "GO_category", "index"], inplace=True
        )

        self.enrichment_go = df
        self.failed_go_enrichments = pd.DataFrame(failed_enrichments)

        # Export results (should be moved to enrichment.py at some point I think)
        with pd.ExcelWriter(out_dir.parent / "enrichment_go.xlsx") as writer:
            df = self.enrichment_go.copy()
            df.reset_index(inplace=True)
            df.to_excel(writer, "go", index=False)
            ws = writer.sheets["go"]
            try:
                ws.autofilter(0, 0, df.shape[0], df.shape[1] - 1)
            except:
                logger.warning("XLS formatting issue.")

    def run_enrichment_kegg(
        self, organism, annot_col="Name", out_dir="enrichment"
    ):  # pragma: no cover

        out_dir = Path(out_dir) / "figures"
        out_dir.mkdir(exist_ok=True, parents=True)

        gene_lists_dict = self.get_gene_lists(annot_col=annot_col, dropna=True)
        enrichment = {}

        for compa in self.comparisons:
            gene_lists = gene_lists_dict[compa]
            ke = KeggPathwayEnrichment(gene_lists, organism, progress=False)
            ke.compute_enrichment()

            for direction in ["up", "down", "all"]:
                enrichment[(compa, direction)] = ke._get_final_df(
                    ke.enrichment[direction].results, nmax=10000
                )
                pylab.figure()
                ke.scatterplot(direction)
                pylab.tight_layout()
                pylab.savefig(out_dir / f"kegg_{compa}_{direction}.pdf")
                pylab.savefig(out_dir / f"kegg_{compa}_{direction}.png")

            logger.info(f"KEGG enrichment for {compa} DONE.")

        df = pd.concat(enrichment).sort_index()
        df.index.rename(["comparison", "direction", "index"], inplace=True)

        self.enrichment_kegg = df

        # Export results (should be moved to enrichment.py at some point I think)
        with pd.ExcelWriter(out_dir.parent / "enrichment_kegg.xlsx") as writer:
            df = self.enrichment_kegg.copy()
            df.reset_index(inplace=True)
            df.to_excel(writer, "kegg", index=False)
            ws = writer.sheets["kegg"]
            try:
                ws.autofilter(0, 0, df.shape[0], df.shape[1] - 1)
            except:
                logger.warning("Fixme")

    def get_gene_lists(
        self, annot_col="index", Nmax=None, dropna=False
    ):  # pragma: no cover

        gene_lists_dict = {}

        for compa in self.comparisons.keys():
            df = self.df.loc[:, ["annotation", compa]].copy()
            df = df.droplevel(0, axis=1)

            fc_filt = df["log2FoldChange"].abs() >= self._log2_fc
            fdr_filt = df["padj"] <= self._alpha

            df = df[fc_filt.values & fdr_filt.values]
            df.reset_index(inplace=True)

            if Nmax:
                df.sort_values("log2FoldChange", ascending=False, inplace=True)
                up_genes = list(df.query("log2FoldChange > 0")[annot_col])[:Nmax]

                df.sort_values("log2FoldChange", ascending=True, inplace=True)
                down_genes = list(df.query("log2FoldChange < 0")[annot_col])[:Nmax]

                all_genes = list(
                    list(
                        df.sort_values("log2FoldChange", key=abs, ascending=False)[
                            annot_col
                        ]
                    )[:Nmax]
                )

            else:
                up_genes = list(df.query("log2FoldChange > 0")[annot_col])
                down_genes = list(df.query("log2FoldChange < 0")[annot_col])
                all_genes = list(df.loc[:, annot_col])

            gene_lists_dict[compa] = {
                "up": up_genes,
                "down": down_genes,
                "all": all_genes,
            }

            if dropna:
                for direction in gene_lists_dict[compa]:
                    gl = gene_lists_dict[compa][direction]
                    if not gl:
                        continue
                    perc_unannotated = gl.count(None) / len(gl) * 100
                    logger.warning(
                        f"{compa} {direction}: Removing {perc_unannotated:.0f}% of the data for enrichment analysis due to missing identifiers in annotation."
                    )
                    gene_lists_dict[compa][direction] = [x for x in gl if x != None]

        return gene_lists_dict

    def _get_design(self, design_file, condition, palette):
        """Import design from a table file and add color groups following the
        groups defined in the column 'condition' of the table file.
        """

        df = RNADesign(design_file).df.set_index("label")

        col_map = dict(zip(df.loc[:, condition].unique(), palette))
        df["group_color"] = df.loc[:, condition].map(col_map)
        return df

    def _format_plot(self, title="", xlabel="", ylabel="", rotation=0):
        pylab.title(title)
        pylab.xticks(rotation=rotation, ha="right")
        pylab.xlabel(xlabel)
        pylab.ylabel(ylabel)

    def get_specific_commons(self, direction, compas=None, annot_col="index"):
        """From all the comparisons contained by the object, extract gene lists which
        are common (but specific, ie a gene appear only appears in the
        combination considered) comparing all combinations of comparisons.

        :param direction: The regulation direction (up, down or all) of the gene
        lists to consider

        :param compas: Specify a list of comparisons to consider (Comparisons
        names can be found with self.comparisons.keys()).

        """

        common_specific_dict = {}

        total_gene_lists = self.get_gene_lists(annot_col=annot_col)

        if not compas:
            compas = self.comparisons.keys()

        for size in range(1, len(compas)):
            for compa_group in combinations(compas, size):
                gene_lists = [
                    total_gene_lists[compa][direction] for compa in compa_group
                ]
                commons = set.intersection(
                    *[set(gene_list) for gene_list in gene_lists]
                )
                other_compas = [compa for compa in compas if compa not in compa_group]
                genes_in_other_compas = {
                    x
                    for other_compa in other_compas
                    for x in total_gene_lists[other_compa][direction]
                }

                commons = commons - genes_in_other_compas
                common_specific_dict[compa_group] = commons

        return common_specific_dict

    def plot_count_per_sample(self, fontsize=12, rotation=45):
        """Number of mapped and annotated reads (i.e. counts) per sample. Each color
        for each replicate

        .. plot::
            :include-source:

            from sequana.rnadiff import RNADiffResults
            from sequana import sequana_data

            r = RNADiffResults(sequana_data("rnadiff/rnadiff_onecond_1"))
            r.plot_count_per_sample()

        """
        pylab.clf()
        df = self.counts_raw.sum().rename("total_counts")
        df = pd.concat([self.design_df, df], axis=1)

        pylab.bar(
            df.index,
            df.total_counts / 1000000,
            color=df.group_color,
            lw=1,
            zorder=10,
            ec="k",
            width=0.9,
        )

        pylab.xlabel("Samples", fontsize=fontsize)
        pylab.ylabel("reads (M)", fontsize=fontsize)
        pylab.grid(True, zorder=0)
        pylab.title("Total read count per sample", fontsize=fontsize)
        pylab.xticks(rotation=rotation, ha="right")
        # pylab.xticks(range(N), self.sample_names)
        try:
            pylab.tight_layout()
        except:
            pass

    def plot_percentage_null_read_counts(self):
        """Bars represent the percentage of null counts in each samples.  The dashed
        horizontal line represents the percentage of feature counts being equal
        to zero across all samples

        .. plot::
            :include-source:

            from sequana.rnadiff import RNADiffResults
            from sequana import sequana_data

            r = RNADiffResults(sequana_data("rnadiff/rnadiff_onecond_1"))
            r.plot_percentage_null_read_counts()

        """
        pylab.clf()
        # how many null counts ?
        df = (self.counts_raw == 0).sum() / self.counts_raw.shape[0] * 100
        df = df.rename("percent_null")
        df = pd.concat([self.design_df, df], axis=1)

        pylab.bar(
            df.index, df.percent_null, color=df.group_color, ec="k", lw=1, zorder=10
        )

        all_null = (self.counts_raw == 0).all(axis=1).sum() / self.counts_raw.shape[0]

        pylab.axhline(all_null, ls="--", color="black", alpha=0.5)

        pylab.xticks(rotation=45, ha="right")
        pylab.ylabel("Proportion of null counts (%)")
        pylab.grid(True, zorder=0)
        pylab.tight_layout()

    def plot_pca(
        self,
        n_components=2,
        colors=None,
        plotly=False,
        max_features=500,
        genes_to_remove=[],
    ):

        """

        .. plot::
            :include-source:

            from sequana.rnadiff import RNADiffResults
            from sequana import sequana_data

            path = sequana_data("rnadiff/rnadiff_onecond_1")
            r = RNADiffResults(path)

            colors = {
                'surexp1': 'r',
                'surexp2':'r',
                'surexp3':'r',
                'surexp1': 'b',
                'surexp2':'b',
                'surexp3':'b'}
            r.plot_pca(colors=colors)
        """
        from sequana.viz import PCA

        # Get most variable genes (n=max_features)
        top_features = (
            self.counts_vst.var(axis=1)
            .sort_values(ascending=False)
            .index[:max_features]
        )

        if genes_to_remove:
            top_features = [x for x in top_features if x not in genes_to_remove]

        counts_top_features = self.counts_vst.loc[top_features, :]

        p = PCA(counts_top_features)

        if plotly is True:
            assert n_components == 3
            variance = p.plot(
                n_components=n_components,
                colors=colors,
                show_plot=False,
                max_features=max_features,
            )
            from plotly import express as px

            df = pd.DataFrame(p.Xr)
            df.index = p.df.columns
            df.columns = ["PC1", "PC2", "PC3"]
            df["size"] = [10] * len(df)  # same size for all points ?

            df = pd.concat([df, self.design_df], axis=1)
            df["label"] = df.index
            df["group_color"] = df[self.condition]

            fig = px.scatter_3d(
                df,
                x="PC1",
                y="PC2",
                z="PC3",
                color="group_color",
                labels={
                    "PC1": "PC1 ({}%)".format(round(100 * variance[0], 2)),
                    "PC2": "PC2 ({}%)".format(round(100 * variance[1], 2)),
                    "PC3": "PC3 ({}%)".format(round(100 * variance[2], 2)),
                },
                height=800,
                text="label",
            )
            return fig
        else:
            variance = p.plot(
                n_components=n_components,
                colors=self.design_df.group_color,
                max_features=max_features,
            )

        return variance

    def plot_mds(self, n_components=2, colors=None, clf=True):
        """IN DEV, not functional"""

        from sequana.viz.mds import MDS

        p = MDS(self.counts_vst)  # [self.sample_names])
        # if colors is None:
        #    colors = {}
        #    for sample in self.sample_names:
        #        colors[sample] = self.colors[self.get_cond_from_sample(sample)]
        p.plot(n_components=n_components, colors=self.design_df.group_color, clf=clf)

    def plot_isomap(self, n_components=2, colors=None):
        """IN DEV, not functional"""

        from sequana.viz.isomap import Isomap

        p = Isomap(self.counts_vst)
        # if colors is None:
        #    colors = {}
        #    for sample in self.sample_names:
        #        colors[sample] = self.colors[self.get_cond_from_sample(sample)]
        p.plot(n_components=n_components, colors=self.design_df.group_color)

    def plot_density(self):
        import seaborn

        seaborn.set()
        for sample in self.counts_raw.columns:
            seaborn.kdeplot(pylab.log10(self.counts_raw[sample].clip(lower=1)))

        self._format_plot(
            title="Count density distribution",
            xlabel="Raw counts (log10)",
            ylabel="Density",
        )

    def plot_feature_most_present(self):
        """"""

        df = []

        for x, y in self.counts_raw.idxmax().iteritems():

            most_exp_gene_count = self.counts_raw.stack().loc[y, x]
            total_sample_count = self.counts_raw.sum().loc[x]

            df.append(
                {
                    "label": x,
                    "gene_id": y,
                    "count": most_exp_gene_count,
                    "total_sample_count": total_sample_count,
                    "most_exp_percent": most_exp_gene_count / total_sample_count * 100,
                }
            )

        df = pd.DataFrame(df).set_index("label")
        df = pd.concat([self.design_df, df], axis=1)

        pylab.clf()
        p = pylab.barh(
            df.index,
            df.most_exp_percent,
            color=df.group_color,
            zorder=10,
            lw=1,
            ec="k",
            height=0.9,
        )

        for idx, rect in enumerate(p):
            pylab.text(
                2,  # * rect.get_height(),
                idx,  # rect.get_x() + rect.get_width() / 2.0,
                df.gene_id.iloc[idx],
                ha="center",
                va="center",
                rotation=0,
                zorder=20,
            )

        self._format_plot(
            # title="Counts monopolized by the most expressed gene",
            # xlabel="Sample",
            xlabel="Percent of total reads",
        )
        pylab.tight_layout()

    def plot_dendogram(
        self,
        max_features=5000,
        transform_method="log",
        method="ward",
        metric="euclidean",
    ):
        # for info about metric and methods: https://tinyurl.com/yyhk9cl8

        assert transform_method in ["log", "anscombe", None]
        # first we take the normalised data
        from sequana.viz import clusterisation
        from sequana.viz import dendogram

        cluster = clusterisation.Cluster(self.counts_norm)
        # cluster = clusterisation.Cluster(self.df[self.sample_names])
        if transform_method is not None:
            data = cluster.scale_data(
                transform_method=transform_method, max_features=max_features
            )
            df = pd.DataFrame(data[0])
            df.index = data[1]
            df.columns = self.counts_norm.columns
        else:
            df = pd.DataFrame(self.counts_norm)
            # df.index = data[1]
            df.columns = self.counts_norm.columns

        d = dendogram.Dendogram(
            df.T,
            metric=metric,
            method=method,
            side_colors=list(self.design_df.group_color.unique()),
        )

        # Convert groups into numbers for Dendrogram category
        group_conv = {
            group: i for i, group in enumerate(self.design_df[self.condition].unique())
        }
        d.category = self.design_df[self.condition].map(group_conv).to_dict()
        d.plot()

    def plot_boxplot_rawdata(self, fliersize=2, linewidth=2, rotation=0, **kwargs):
        import seaborn as sbn

        ax = sbn.boxplot(
            data=self.counts_raw.clip(1),
            linewidth=linewidth,
            fliersize=fliersize,
            palette=self.design_df.group_color,
            **kwargs,
        )
        pos, labs = pylab.xticks()
        pylab.xticks(pos, labs, rotation=rotation)
        ax.set_ylabel("Counts (raw) in log10 scale")
        ax.set_yscale("log")
        self._format_plot(ylabel="Raw count distribution")
        pylab.tight_layout()

    def plot_boxplot_normeddata(self, fliersize=2, linewidth=2, rotation=0, **kwargs):
        import seaborn as sbn

        ax = sbn.boxplot(
            data=self.counts_norm.clip(1),
            linewidth=linewidth,
            fliersize=fliersize,
            palette=self.design_df.group_color,
            **kwargs,
        )
        pos, labs = pylab.xticks()
        pylab.xticks(pos, labs, rotation=rotation)
        ax.set(yscale="log")
        self._format_plot(ylabel="Normalised count distribution")
        pylab.tight_layout()

    def plot_dispersion(self):

        pylab.plot(
            self.dds_stats.baseMean,
            self.dds_stats.dispGeneEst,
            "ok",
            label="Estimate",
            ms=1,
        )
        pylab.plot(
            self.dds_stats.baseMean,
            self.dds_stats.dispersion,
            "ob",
            label="final",
            ms=1,
        )
        pylab.plot(
            self.dds_stats.baseMean, self.dds_stats.dispFit, "or", label="Fit", ms=1
        )
        pylab.legend()
        ax = pylab.gca()
        ax.set(yscale="log")
        ax.set(xscale="log")

        self._format_plot(
            title="Dispersion estimation",
            xlabel="Mean of normalized counts",
            ylabel="Dispersion",
        )

    def heatmap(self, comp, log2_fc=1, padj=0.05):
        assert comp in self.comparisons.keys()
        from sequana.viz import heatmap

        h = heatmap.Clustermap(
            self.counts_norm.loc[
                self.comparisons[comp]
                .df.query(
                    "(log2FoldChange<-@log2_fc or log2FoldChange>@log2_fc) and padj<@padj"
                )
                .index
            ]
        ).plot()

max_stars_repo_name: ddesvillechabrol/sequana

max_issues_repo_name: ddesvillechabrol/sequana

max_forks_repo_name: ddesvillechabrol/sequana


text:
import numpy as np
# from scipy.linalg import logm, expm
from math import cos, sin, pi
import matplotlib.pyplot as plt
import utils

def generate_point_cloud():
    mean = np.array([0, 0, 0.])

    # R = utils.euler_angle_to_rotation(0.4, -1., 6.6666)
    R = utils.euler_angle_to_rotation(0.3, 0.5, 0.0)
    # Transformation for covariance
    cov = R @ np.diag([0.001, 1, 9.]) @ R.transpose()
    num_point = 100
    points = np.random.multivariate_normal(mean, cov, size=num_point)

    points_mean = np.mean(points, axis=0)
    print('points_mean:', points_mean)
    points = points - points_mean

    np.testing.assert_almost_equal(R.transpose() @ R, np.identity(3))
    print('R_gt: ', R)
    print('cov_gt:', cov)
    print('numpy.linalg.eig(cov):', np.linalg.eig(cov))

    return points.transpose()

def take_n_cols(R, n):
    return R[:, :n]

class PCA_SO3_projection_minimization:
    """
    https://stats.stackexchange.com/questions/10251/what-is-the-objective-function-of-pca
    """
    def __init__(self, points):
        self.num_data = points.shape[1]
        self.num_residuals = points.size
        points_mean = np.mean(points, axis=1)
        self.points = (points.transpose() - points_mean).transpose()
        self.points_covariance = points @ points.transpose() / self.num_data
        # print(self.points_covariance)
        self.var_SO3 = np.identity(3)
        self.rank = 1
        self.var_projection = np.zeros([self.rank, self.num_data])

    def residaul(self, var_SO3, var_projection):
        """
          r_i = p_i - first_k_cols(R) * w
        """
        # self.point.shape == (3, n)
        r = self.points - take_n_cols(var_SO3, self.rank) @ var_projection
        # make r col major
        return r.transpose().flatten()

    def compute_projection(self):
        self.var_projection = take_n_cols(self.var_SO3, self.rank).transpose() @ self.points

    def numerical_jacobi_wrt_so3(self):
        """
            r_i = p_i - first_k_cols(R) * w
        """
        DELTA = 1e-8

        jacobian = np.zeros([self.num_residuals, 3])

        w_so3_local = np.array([0, 0, 0.])
        curret_params = w_so3_local.copy()
        for p_idx in range(3):
            params_plus = curret_params.copy()
            params_plus[p_idx] += DELTA
            residual_plus = self.residaul(self.var_SO3 @ utils.so3_exp(params_plus), self.var_projection)

            params_minus = curret_params.copy()
            params_minus[p_idx] -= DELTA
            residual_minus = self.residaul(self.var_SO3 @ utils.so3_exp(params_minus), self.var_projection)

            dr_dpidx = (residual_plus - residual_minus) / (2. * DELTA)
            jacobian[:, p_idx] = dr_dpidx

        return jacobian

    def jacobi_wrt_so3(self):
        """
        The derivation is here,
            r_i(w) = p_i - first_k_cols(R) * proj
                   = p_i - first_k_cols(R * exp(w)) * proj
                     (3,1)            (3, k)        (k, 1)
            let f3 = - first_k_cols(R * exp(w)) * proj
                f2 = first_k_cols(R * exp(w))
            f3 = - f2 * proj
            f3 =  - (proj.T * f2.T).T

            df3/ dw = - (proj.T * f2dw.T).T,  (3, 3)
                         (1,k)   ((k,3), 3) = (( 3*1), 3)

            let f2 = first_k_cols(R * exp(w))
                f1 = R * exp(w)

            df2/dw = df2/df1 * df1/dw, ((3,k), 3)

            df2/df1 = [1 when taking the element, 0 otherwise],   ((3,k), (3,3))

            df1/dw = [R * G1 , R * G2, R * G3], ((3,3), 3)
        """

        G1 = np.array([
            [0, 0, 0],
            [0, 0, -1],
            [0, 1, 0.],
        ])
        G2 = np.array([
            [0, 0, 1],
            [0, 0, 0],
            [-1, 0, 0.],
        ])
        G3 = np.array([
            [0, -1, 0],
            [1, 0, 0],
            [0, 0, 0.],
        ])
        f1 = np.zeros([3,3,3])
        f1[:, :, 0] = self.var_SO3 @ G1
        f1[:, :, 1] = self.var_SO3 @ G2
        f1[:, :, 2] = self.var_SO3 @ G3

        # ((k, 3), 3)
        f2_t = np.zeros([self.rank, 3, 3])
        for i in range(3):
            f2_t[:, :, i] = f1[:, :self.rank, i].transpose()

        jacobi_all_points = np.zeros([self.num_data, 3, 3])
        for i in range(3):
            #          ((n,3), 3)      =                (n, k)                  @    ((k, 3), 3)
            jacobi_all_points[:, :, i] = - (self.var_projection.transpose() @ f2_t[:, :, i])

        jacobi = jacobi_all_points.reshape(self.num_data * 3, 3)
        return jacobi

    def solve_normal_equation_and_update_wrt_so3(self):
        """
        """
        jacobi = self.jacobi_wrt_so3()
        # jacobi = self.numerical_jacobi_wrt_so3()
        # print('jacobian', jacobi)
        r = self.residaul(self.var_SO3, self.var_projection)

        # rhs is invertable when rank == 1
        regulization = 1e-6
        rhs = jacobi.transpose() @ jacobi + regulization * np.identity(3)
        lhs = - jacobi.transpose() @ r
        # print('rhs:', rhs)
        # print('lhs:', lhs)
        delta_so3 = np.linalg.solve(rhs, lhs)
        print('delta_so3:', delta_so3)
        self.var_SO3 = self.var_SO3 @ utils.so3_exp(delta_so3)

    def cost(self):
        r = self.residaul(self.var_SO3, self.var_projection)
        return r.transpose() @ r

    def print_variable(self):
        print('cost:', self.cost())
        np.testing.assert_almost_equal(self.var_SO3.transpose() @ self.var_SO3, np.identity(3))
        print('R_k:', take_n_cols(self.var_SO3, self.rank).transpose())
        for i in range(self.rank):
            p = self.var_projection[i, :]
            print(p.transpose() @ p / self.num_data)

    def minimize(self):
        for iter in range(10):
            if self.cost() < 1e-8:
                break
            self.print_variable()
            self.compute_projection()
            self.solve_normal_equation_and_update_wrt_so3()

            np.testing.assert_almost_equal(self.var_SO3.transpose() @ self.var_SO3, np.identity(3))

def main():
    points = generate_point_cloud()

    # pca = PCA_SO3_Covariance(points)
    # print("pca.cost():", pca.cost())
    # print("pca jocibian:", pca.numerical_jacobi())
    # pca.minimize()

    pca_2 = PCA_SO3_projection_minimization(points)

    pca_2.compute_projection()
    # pca_2.var_projection = np.ones([1, 1])
    # j1 = pca_2.jacobi_wrt_so3()
    # j2 =pca_2.numerical_jacobi_wrt_so3()
    # print('j1:', j1)
    # print('j2:', j2)
    # 1 / 0

    print("pca_2.cost()", pca_2.cost())
    pca_2.minimize()
    print('finial so3:', pca_2.var_SO3)

    cov_stats = points @ points.transpose() / points.shape[1]
    # print('cov_stats:', cov_stats)
    e, v = np.linalg.eig(cov_stats)
    idx = np.argsort(e)[::-1]
    e = e[idx]
    v = v[:,idx]
    print('e(cov_stats):', e)
    print('v(cov_stats):', v)

    pca_using_eig = PCA_SO3_projection_minimization(points)
    pca_using_eig.var_SO3 = v
    pca_using_eig.compute_projection()
    pca_using_eig.print_variable()
    print("pca_using_eig.cost() by svd:", pca_using_eig.cost())

if __name__ == "__main__":
    main()

ext: py

max_stars_repo_name: yimuw/yimu-blog

max_issues_repo_name: yimuw/yimu-blog

max_forks_repo_name: yimuw/yimu-blog

reason: import numpy,from scipy


text:
import time
from collections import Counter
import warnings; warnings.filterwarnings('ignore')

import numpy as np
import random
import pandas as pd
import matplotlib.pyplot as plt

from algorithms import ShapeletTransformer
from extractors.extractor import MultiGeneticExtractor
from data.load_all_datasets import load_data_train_test

from sklearn.metrics import accuracy_score, log_loss
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV

from tslearn.shapelets import ShapeletModel

def grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_shapelets, l, r):
    base_size = int(l * ts_sz)
    d = {}
    for sz_idx in range(r):
        shp_sz = base_size * (sz_idx + 1)
        d[shp_sz] = n_shapelets
    return d

np.random.seed(1337)
random.seed(1337)

TRAIN_PATH = 'data/partitioned/MoteStrain/MoteStrain_train.csv'
TEST_PATH = 'data/partitioned/MoteStrain/MoteStrain_test.csv'

# Load the training and testing dataset (features + label vector)
train_df = pd.read_csv(TRAIN_PATH)
test_df = pd.read_csv(TEST_PATH)
X_train = train_df.drop('target', axis=1).values
y_train = train_df['target']
X_test = test_df.drop('target', axis=1).values
y_test = test_df['target']

map_dict = {}
for j, c in enumerate(np.unique(y_train)):
    map_dict[c] = j
y_train = y_train.map(map_dict)
y_test = y_test.map(map_dict)

print(set(y_train), set(y_test))

y_train = y_train.values
y_test = y_test.values

nr_shap, l, r, reg, max_it = 0.05, 0.075, 2, 0.01, 5000

measurements = []

for _nr_shap in np.arange(0.05, 0.55, 0.05):

    shapelet_dict = grabocka_params_to_shapelet_size_dict(
            X_train.shape[0], X_train.shape[1], int(_nr_shap*X_train.shape[1]), l, r
    )
    clf = ShapeletModel(n_shapelets_per_size=shapelet_dict,
                        max_iter=max_it, verbose_level=0, batch_size=8,
                        optimizer='sgd', weight_regularizer=reg)

    start = time.time()
    clf.fit(
        np.reshape(
            X_train,
            (X_train.shape[0], X_train.shape[1], 1)
        ),
        y_train
    )
    learning_time = time.time() - start

    X_distances_train = clf.transform(X_train)
    X_distances_test = clf.transform(X_test)

    lr = GridSearchCV(LogisticRegression(), {'penalty': ['l1', 'l2'], 'C': [0.001, 0.01, 0.1, 1.0, 10.0]})
    lr.fit(X_distances_train, y_train)

    print('[K]', _nr_shap, accuracy_score(y_test, lr.predict(X_distances_test)))

    measurements.append([_nr_shap, accuracy_score(y_test, lr.predict(X_distances_test))])

plt.figure()
plt.plot([x[0] for x in measurements], [x[1] for x in measurements])
plt.title('Sensitive of LTS on hyper-parameter K')
plt.show()

measurements_df = pd.DataFrame(measurements, columns=['nr_shap', 'accuracy'])
measurements_df.to_csv('lts_param_K.csv')

measurements = []

for _l in np.arange(0.075, 0.4, 0.075):

    shapelet_dict = grabocka_params_to_shapelet_size_dict(
            X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), _l, r
    )
    clf = ShapeletModel(n_shapelets_per_size=shapelet_dict,
                        max_iter=max_it, verbose_level=0, batch_size=8,
                        optimizer='sgd', weight_regularizer=reg)

    start = time.time()
    clf.fit(
        np.reshape(
            X_train,
            (X_train.shape[0], X_train.shape[1], 1)
        ),
        y_train
    )
    learning_time = time.time() - start

    X_distances_train = clf.transform(X_train)
    X_distances_test = clf.transform(X_test)

    lr = GridSearchCV(LogisticRegression(), {'penalty': ['l1', 'l2'], 'C': [0.001, 0.01, 0.1, 1.0, 10.0]})
    lr.fit(X_distances_train, y_train)

    print('[L', _l, accuracy_score(y_test, lr.predict(X_distances_test)))

    measurements.append([_l, accuracy_score(y_test, lr.predict(X_distances_test))])

plt.figure()
plt.plot([x[0] for x in measurements], [x[1] for x in measurements])
plt.title('Sensitive of LTS on hyper-parameter L (ItalyPowerDemand)')
plt.show()

measurements_df = pd.DataFrame(measurements, columns=['l', 'accuracy'])
measurements_df.to_csv('lts_param_L.csv')

measurements = []

for _r in range(1, 6):

    shapelet_dict = grabocka_params_to_shapelet_size_dict(
            X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, _r
    )
    clf = ShapeletModel(n_shapelets_per_size=shapelet_dict,
                        max_iter=max_it, verbose_level=0, batch_size=8,
                        optimizer='sgd', weight_regularizer=reg)

    start = time.time()
    clf.fit(
        np.reshape(
            X_train,
            (X_train.shape[0], X_train.shape[1], 1)
        ),
        y_train
    )
    learning_time = time.time() - start

    X_distances_train = clf.transform(X_train)
    X_distances_test = clf.transform(X_test)

    lr = GridSearchCV(LogisticRegression(), {'penalty': ['l1', 'l2'], 'C': [0.001, 0.01, 0.1, 1.0, 10.0]})
    lr.fit(X_distances_train, y_train)

    print('[R]', _r, accuracy_score(y_test, lr.predict(X_distances_test)))

    measurements.append([_r, accuracy_score(y_test, lr.predict(X_distances_test))])

plt.figure()
plt.plot([x[0] for x in measurements], [x[1] for x in measurements])
plt.title('Sensitive of LTS on hyper-parameter R (ItalyPowerDemand)')
plt.show()

measurements_df = pd.DataFrame(measurements, columns=['r', 'accuracy'])
measurements_df.to_csv('lts_param_R.csv')

measurements = []

for _reg in np.arange(0., 0.55, 0.05):

    shapelet_dict = grabocka_params_to_shapelet_size_dict(
            X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
    )
    clf = ShapeletModel(n_shapelets_per_size=shapelet_dict,
                        max_iter=max_it, verbose_level=0, batch_size=8,
                        optimizer='sgd', weight_regularizer=_reg)

    start = time.time()
    clf.fit(
        np.reshape(
            X_train,
            (X_train.shape[0], X_train.shape[1], 1)
        ),
        y_train
    )
    learning_time = time.time() - start

    X_distances_train = clf.transform(X_train)
    X_distances_test = clf.transform(X_test)

    lr = GridSearchCV(LogisticRegression(), {'penalty': ['l1', 'l2'], 'C': [0.001, 0.01, 0.1, 1.0, 10.0]})
    lr.fit(X_distances_train, y_train)

    print('[REG]', _reg, accuracy_score(y_test, lr.predict(X_distances_test)))

    measurements.append([_reg, accuracy_score(y_test, lr.predict(X_distances_test))])

plt.figure()
plt.plot([x[0] for x in measurements], [x[1] for x in measurements])
plt.title('Sensitive of LTS on hyper-parameter lambda (ItalyPowerDemand)')
plt.show()

measurements_df = pd.DataFrame(measurements, columns=['reg', 'accuracy'])
measurements_df.to_csv('lts_param_lambda.csv')

measurements = []

for _max_it in range(1000, 11000, 1000):

    shapelet_dict = grabocka_params_to_shapelet_size_dict(
            X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
    )
    clf = ShapeletModel(n_shapelets_per_size=shapelet_dict,
                        max_iter=_max_it, verbose_level=0, batch_size=8,
                        optimizer='sgd', weight_regularizer=reg)

    start = time.time()
    clf.fit(
        np.reshape(
            X_train,
            (X_train.shape[0], X_train.shape[1], 1)
        ),
        y_train
    )
    learning_time = time.time() - start

    X_distances_train = clf.transform(X_train)
    X_distances_test = clf.transform(X_test)

    lr = GridSearchCV(LogisticRegression(), {'penalty': ['l1', 'l2'], 'C': [0.001, 0.01, 0.1, 1.0, 10.0]})
    lr.fit(X_distances_train, y_train)

    print('[Max Iter]', _max_it, accuracy_score(y_test, lr.predict(X_distances_test)))

    measurements.append([_max_it, accuracy_score(y_test, lr.predict(X_distances_test))])

plt.figure()
plt.plot([x[0] for x in measurements], [x[1] for x in measurements])
plt.title('Sensitive of LTS on hyper-parameter max. iter. (ItalyPowerDemand)')
plt.show()

measurements_df = pd.DataFrame(measurements, columns=['maxiter', 'accuracy'])
measurements_df.to_csv('lts_param_maxiter.csv')

max_stars_repo_name: GillesVandewiele/pyShapelets

max_issues_repo_name: GillesVandewiele/pyShapelets

max_forks_repo_name: GillesVandewiele/pyShapelets

reason: import numpy


text:
# -*- coding: utf-8 -*-
'''
#-------------------------------------------------------------------------------
# NATIONAL UNIVERSITY OF SINGAPORE - NUS
# SINGAPORE INSTITUTE FOR NEUROTECHNOLOGY - SINAPSE
# Singapore
# URL: http://www.sinapseinstitute.org
#-------------------------------------------------------------------------------
# Neuromorphic Engineering Group
# Author: Andrei Nakagawa-Silva, MSc
# Contact: nakagawa.andrei@gmail.com
#-------------------------------------------------------------------------------
# Description: Class for handling socket communication with the HD Neuromorphic
# Array. The original software was written in C#. Therefore, we are using sockets
# to read the array from C# and then sending the data via sockets to a Python
# server where shape recognition is implemented.
#-------------------------------------------------------------------------------
'''
#-------------------------------------------------------------------------------
import os, sys, glob
sys.path.append('../general')
import numpy as np
import pyqtgraph as pg
from collections import deque #necessary for acquisition
from copy import copy #useful for copying data structures
from threading import Lock #control access in threads
from threadhandler import ThreadHandler #manage threads
from sockethandler import * #import socket communication handlers
#-------------------------------------------------------------------------------
class HDNerArray():
    def __init__(self,_ip='127.0.0.1',_port=8888):
        self.ip = _ip #ip of host
        self.port = _port #port number
        self.server = None
        self.NROWS = 64
        self.NCOLS = 64
    #start the server
    def start(self):
        self.server = SocketServerHandler(self.ip,self.port)
    #retrieve the data queue, process the incoming data
    #and convert it to a queue containing spiking events
    #with the following structure
    #[SPIKE_TIME][POS_X][POS_Y][SPIKE_POLARITY]
    def getData(self):
        if self.server is not None:
            q = self.server.getDataQueue()
            self.dataQueue = deque()
            numData = len(q)
            if numData > 0:
                for k in range(numData):
                    #retrieve the spike events contained in the data queue
                    #the join function concatenates to a list the
                    #byte values converted to chr, forming a string
                    #that is being split by new lines (\n)
                    spikeEvents = "".join(map(chr,q.popleft())).split('\n')
                    #for every event, do a post-processing in the data
                    for z in range(len(spikeEvents)-1):
                        correctedData = []
                        #convert the data to int
                        [correctedData.append(int(x)) for x in spikeEvents[z].split(' ')]
                        #add to the class queue
                        self.dataQueue.append(correctedData)
                        #debugging
                        #print(correctedData) #print the read values
                return copy(self.dataQueue) #return a copy of the queued spike event
            else:
                return False
#-------------------------------------------------------------------------------
#Run the appo
if __name__ == '__main__':
    #create an object
    hdarray = HDNerArray()
    def read():
        global hdarray
        z = hdarray.getData()
        n = len(z)
        for k in range(n):
            print(z.popleft())
    #thread to read
    th = ThreadHandler(read)
    #waits for input to start
    input('press ENTER to start...')
    hdarray.start()
    th.start()
    #wait for input to finish
    input('press ENTER to finish...')
    th.kill()


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import tensorflow as tf

from vocab import Vocab
from lib.models import NN

class BaseSimpleSrlers(NN):

    def __call__(self, dataset, moving_paras=None):
        raise NotImplementedError

    def parse_argmax(self):
        raise NotImplementedError

    # =============================================================
    @property
    def input_idxs(self):
        return (0, 1, 2, 3, 4)

    @property
    def target_idxs(self):
        return (5)

max_stars_repo_name: Impavidity/neural-srl

max_issues_repo_name: Impavidity/neural-srl

max_forks_repo_name: Impavidity/neural-srl


text:
# -*- coding: utf-8 -*-

"""
fairsearchcore.mtable_generator
~~~~~~~~~~~~~~~
Contains the mechanics for creating an mtable
"""

import pandas as pd
import scipy.stats as stats

from fairsearchcore import fail_prob

class MTableGenerator:

    def __init__(self, k, p, alpha, adjust_alpha):
        # assign parameters
        self.k = k
        self.p = p
        self.alpha = alpha
        self.adjust_alpha = adjust_alpha

        if self.adjust_alpha:
            fail_prob_pair = fail_prob.RecursiveNumericFailProbabilityCalculator(k, p, alpha).adjust_alpha()
            self.adjusted_alpha = fail_prob_pair.alpha
            self._mtable = fail_prob_pair.mtable
        else:
            self.adjusted_alpha = alpha
            self._mtable = self._compute_mtable()

    def mtable_as_list(self):
        return [int(i) for i in self._mtable.m.tolist()]

    def mtable_as_dataframe(self):
        return self._mtable

    def m(self, k):
        if k < 1:
            raise ValueError("Parameter k must be at least 1")
        elif k > self.k:
            raise ValueError("Parameter k must be at most {0}".format(self.k))

        result = stats.binom.ppf(self.adjusted_alpha if self.adjust_alpha else self.alpha, k, self.p)
        return 0 if result < 0 else result

    def _compute_mtable(self):
        """ Computes a table containing the minimum number of protected elements
            required at each position
        """
        mtable = pd.DataFrame(columns=["m"])
        for i in range(1, self.k + 1):
            if i % 2000 == 0:
                print("Computing m: {:.0f} of {:.0f}".format(i, self.k))
            mtable.loc[i] = [self.m(i)]
        return mtable

def compute_aux_mtable(mtable):
    """
    Stores the inverse of an mTable entry and the size of the block with respect to the inverse
    """
    if not (isinstance(mtable, pd.DataFrame)):
        raise TypeError("Internal mtable must be a DataFrame")

    aux_mtable = pd.DataFrame(columns=["inv", "block"])
    last_m_seen = 0
    last_position = 0
    for position in range(1, len(mtable)):
        if position % 2000 == 0:
            print("Computing m inverse: {:.0f} of {:.0f}".format(position, len(mtable)))
        if mtable.at[position, "m"] == last_m_seen + 1:
            last_m_seen += 1
            aux_mtable.loc[position] = [position, position - last_position]
            last_position = position
        elif mtable.at[position, "m"] != last_m_seen:
            raise RuntimeError("Inconsistent mtable")

    return aux_mtable

max_stars_repo_name: utkarshmttl/fairsearch-fair-python

max_issues_repo_name: utkarshmttl/fairsearch-fair-python

max_forks_repo_name: utkarshmttl/fairsearch-fair-python

reason: import scipy


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

class SolveLLS():

    def __init__(self, y_tr_norm, X_tr_norm, y_te_norm, X_te_norm, sy, my):
        self.y_tr_norm = y_tr_norm
        self.X_tr_norm = X_tr_norm
        self.y_te_norm = y_te_norm
        self.X_te_norm = X_te_norm
        self.sy = sy
        self.my = my

    def LLS_With_Diagrams(self):
        w_hat = np.linalg.inv(self.X_tr_norm.T@self.X_tr_norm)@(self.X_tr_norm.T@self.y_tr_norm)
        y_hat_te_norm = self.X_te_norm@w_hat
        MSE_norm = np.mean((y_hat_te_norm - self.y_te_norm)**2)
        MSE=self.sy**2*MSE_norm

        # plot the regression line
        y_hat_te=(self.X_te_norm@w_hat)*self.sy+self.my
        y_te=self.y_te_norm*self.sy+self.my
        plt.figure(figsize=(6,4))
        plt.plot(y_te,y_hat_te,'.')
        v=plt.axis()
        plt.plot([v[0],v[1]],[v[0],v[1]],'r',linewidth=2)
        plt.xlabel(r'$y$')
        plt.ylabel(r'$\^y$')
        plt.grid()
        plt.title('Linear Least Squares: test')
        plt.tight_layout()
        plt.savefig('./LLS-yhat_vs_y.png')
        plt.show()

        # plot the error histogram
        E_tr = (self.y_tr_norm-self.X_tr_norm@w_hat)*self.sy# training
        E_te = (self.y_te_norm-self.X_te_norm@w_hat)*self.sy# test
        e = [E_tr,E_te]
        plt.figure(figsize=(6,4))
        plt.hist(e,bins=50,density=True, histtype='bar',
                 label=['training','test'])
        plt.xlabel(r'$e=y-\^y$')
        plt.ylabel('P(error in bin)')
        plt.legend()
        plt.grid()
        plt.title('LLS-Error histograms')
        plt.tight_layout()
        plt.savefig('./LLS-hist.png')
        plt.show()

        E_tr_mu = E_tr.mean()
        E_tr_sig = E_tr.std()
        E_tr_MSE = np.mean(E_tr**2)
        y_tr = self.y_tr_norm * self.sy + self.my
        R2_tr = 1 - E_tr_sig**2 / np.mean(y_tr**2)
        E_te_mu = E_te.mean()
        E_te_sig = E_te.std()
        E_te_MSE = np.mean(E_te**2)
        y_te = self.y_te_norm * self.sy + self.my
        R2_te = 1 - E_te_sig**2 / np.mean(y_te**2)
        rows = ['Training', 'test']
        cols = ['mean', 'std', 'MSE', 'R^2']
        p = np.array([[E_tr_mu, E_tr_sig, E_tr_MSE, R2_tr],
                      [E_te_mu, E_te_sig, E_te_MSE, R2_te]])
        results = pd.DataFrame(p, columns =cols, index=rows)
        print(results)

max_stars_repo_path: Gaussian Process Regression on Parkinson disease/bobonazarov_s267331_report_2_code/sub/algorithms.py

max_stars_repo_name: Russeal/ICT-for-Health

max_issues_repo_name: Russeal/ICT-for-Health

max_forks_repo_name: Russeal/ICT-for-Health


text:
# Standard Libraries
import numpy.random
import tensorflow as tf
import random
import os

# 3rd Party Libraries
import cv2

def makedirs(path):
    """ Try to create the directory, pass if the directory exists already, fails otherwise.
    :param path:            <string>            directory path, that should be created

    """
    try:
        os.makedirs(path, exist_ok=True)
    except OSError:
        if not os.path.isdir(path):
            raise

def layer_to_index(target_layer, model):
    """
    Turns the identifier of an layer into a
    layer index, regardless if it is a str or already index

    :param target_layer: <str or int> identifier for the target layer. None will return None
    :param model: <keras.Model> The model containing the target layer

    :returns: <int> index of target layer
    """
    if target_layer is None:
        return None
    if isinstance(target_layer, str):
        for idx, layer in enumerate(model.layers):
            if layer.name == target_layer:
                return idx
    elif isinstance(target_layer, int):
        return target_layer
    else:
        raise TypeError("layer has to be int or str")

    raise Exception("Layer %s could not be found"%(str(target_layer)))

def tb_write_images(callback, names, imgs):
    """
    :param callback: <tensorflow.python.keras.callbacks.TensorBoard> Tensorboard callback
    :param names: <list of str> headings for the iamges
    :param imgs: <list of numpy.array> Images as Bitmaps
    """
    for name, img in zip(names, imgs):
        tf_img_enc = cv2.imencode('.jpg', img)[1].tostring()
        tf_img_tensor = tf.image.decode_jpeg(tf_img_enc, channels=3)
        tf_img_tensor = tf.expand_dims(tf_img_tensor, axis=0)

        imgsumarry=tf.summary.image(name, tf_img_tensor, max_outputs=3, collections=None)
        callback.writer.add_summary(imgsumarry.eval(session=tf.Session()))
        callback.writer.flush()

def tb_write_texts(callback, names, texts):
    """
    :param callback: <tensorflow.python.keras.callbacks.TensorBoard> Tensorboard callback
    :param names: <list of str> headings for the texts
    :param imgs: <list of str> Contents
    """

    for name, text in zip(names, texts):
        txtsummary = tf.summary.text(name, tf.convert_to_tensor(text))
        callback.writer.add_summary(txtsummary.eval(session=tf.keras.backend.get_session()))
        callback.writer.flush()

def initialize_seed(seed=0):
    """
    This makes experiments more comparable by
    forcing the random number generator to produce
    the same numbers in each run
    """
    random.seed(a=seed)
    numpy.random.seed(seed)

    if hasattr(tf, 'set_random_seed'):
        tf.set_random_seed(seed)
    elif hasattr(tf.random, 'set_random_seed'):
        tf.random.set_random_seed(seed)
    elif hasattr(tf.random, 'set_seed'):
        tf.random.set_seed(seed)
    else:
        raise AttributeError("Could not set seed for TensorFlow")

def get_session(gpu_usage=None):
    """ Construct a modified tf session.

    :param gpu_usage:       <float>             GPU memory usage from 0 to 1, None for dynamic growth

    :return tf.Session:     <tf.Session>        Tensorflow Session object
    """

    if gpu_usage:
        gpu_options = tf.compat.v1.GPUOptions(per_process_gpu_memory_fraction=gpu_usage, allow_growth=True)
    else:
        gpu_options = tf.compat.v1.GPUOptions(allow_growth=True)

    config = tf.compat.v1.ConfigProto(gpu_options=gpu_options)

    return tf.compat.v1.Session(config=config)

def output_index_by_name(model, output_name):
    """
    :param model: the keras model
    :param output_name: the string name for a specific output

    :returns: <int> specifying the index of the requested output
    """
    name_to_index = {name:i for i,name in enumerate(model.output_names)}
    return name_to_index[output_name]

def input_index_by_name(model, input_name):
    """
    :param model: the keras model
    :param output_name: the string name for a specific output

    :returns: <int> specifying the index of the requested output
    """
    name_to_index = {name:i for i,name in enumerate(model.input_names)}
    return name_to_index[input_name]


text:
# coding: utf-8

# ## DIGIT RECOGNITION CHALLENGE

# In[1]:

#import necessary libraries
import numpy as np
import pandas as pd
import matplotlib as mpl
import sklearn as skl

import matplotlib.pyplot as plt

# In[2]:

#import datasets. It is already separated in test and train. I am going to use only the "train" set, and then divide it,
#as I want to be able to know how good is the performance of my model ("test" does not include labels)
dataset = pd.read_csv('C:/Users/Marta/Desktop/Digit Recognition mnist Challenge/train.csv')

# In[3]:

#I want to see the size
print(dataset.shape)

# In[4]:

#look at the datasets
print(type(dataset))
print(dataset.head(5))
print(dataset.shape)
print(dataset.dtypes)
print(dataset.describe())
#we realize it is a set of pixels with a number that indicates if it is coloured and how intensely
#The numbers that indicates coloured is an integer

# In[5]:

#Many pixels has never been coloured (is expected, noone draw a number in the vey corner of the space e.g.)
print(dataset['pixel5'].values.sum())
print(dataset['pixel500'].values.sum())
dataset.sum()

# In[6]:

#Lets check if there are missing values, null, na or nan. There are not!
dataset.isnull().sum().sum()

# In[7]:

#Lets divide now our dataset into a train set and a test set, 80-20%
np.random.seed(7) #I want to obtain the same results when I rerun this, so I put a seed
indicator = np.random.rand(len(dataset),) < 0.8

# In[8]:

#I divide the dataset
train = dataset[indicator]
test = dataset[~indicator]

# In[9]:

#I want to see if the distribution of actual numbers is even
labelsarray = dataset.iloc[:,0]
labelsarray.hist()
plt.show()

# In[10]:

#A fun way to see the numbers!!

a = dataset.loc[40].values[1:]  #select an instance from the dataset, without label
type(a) #it is an array
a.shape #because of the size, 784=2*2*2*2*7*7, I imagine that they are painting in a square 28*28
b = a.reshape(28,28)  #we put them in a matrix (but when printing this, it cuts and looks weird...)
#print(b)
for i in range(28):  #little loop to keep the shape! you could also print *  instead of numbers
    for j in range(28):
        if j==0:
            print("\n",end="")
        if b[i,j] == 0:
            print("   ",end="")
        else:
            print(b[i,j],end=" ")

            #note that the numbers in the border tends to be smaller than the ones in the center, because of the pen preasure

# In[11]:

#I will prepare my dataframe for the classifier algorithm; I separate between features and labels
print(train.shape)
features_train = train.iloc[:,1:]
print(features_train.shape)
labels_train = train.iloc[:,0]
print(labels_train.shape)

features_test = test.iloc[:,1:]
labels_test = test.iloc[:,0]

# ## First approach; Decision Tree

# In[12]:

#I import my method
from sklearn.tree import DecisionTreeClassifier

# In[13]:

#lets define a decision tree algorithm for this case   (THE PARAMETERS I AM USING, I CHOOSED THEM BY A RULE OF THUMB, TRYING)
mytree = DecisionTreeClassifier(criterion = 'entropy',random_state = 7, max_depth=12, min_samples_leaf=9)

# In[14]:

#Train the model!
mytree.fit(features_train,labels_train)

# In[15]:

#Lets try to predict with it
predictions = mytree.predict(features_test)

# In[16]:

#just a first look...
labels_test.values

# In[17]:

#Accuracy within the training set AS THE LABELS ARE EVENLY DISTRIBUTED, THIS IS A GOOD MEASURE OF HOW GOOD IS MY MODEL
actrain = sum(mytree.predict(features_train)==labels_train.values)/len(mytree.predict(features_train))
#Accuracy
actest = sum(predictions==labels_test.values)/len(predictions)
print(actrain)
print(actest)

# ##### I am going to change the dataset; I do not want to see a scale of intensity anymore, I will change it to 0 and 1

# In[18]:

#preparing dataset again!
newdataset = dataset[:]
newdataset[newdataset > 200] = 1 #This avoid little noise around image and also avoid changes in labels

# In[19]:

#Dividing dataset again
np.random.seed(7) #I want to obtain the same results when I rerun this, so I put a seed
newindicator = np.random.rand(len(dataset),) < 0.8

# In[20]:

#I divide the dataset
newtrain = newdataset[newindicator]
newtest = newdataset[~newindicator]

# In[21]:

#preparing dataset again!
features_newtrain = newtrain.iloc[:,1:]
labels_newtrain = newtrain.iloc[:,0]

features_newtest = newtest.iloc[:,1:]
labels_newtest = newtest.iloc[:,0]

# In[22]:

#Train the new model!
mytree.fit(features_newtrain,labels_newtrain)

# In[23]:

#Lets try to predict with it
newpredictions = mytree.predict(features_newtest)

# In[24]:

#Accuracy within the new training set AS THE LABELS ARE EVENLY DISTRIBUTED, THIS IS A GOOD MEASURE OF HOW GOOD IS MY MODEL
newactrain = sum(mytree.predict(features_newtrain)==labels_newtrain.values)/len(mytree.predict(features_newtrain))
#Accuracy
newactest = sum(newpredictions==labels_newtest.values)/len(newpredictions)
print(newactrain)
print(newactest)

# ##### Same accuracy...

# ## Just for fun, I would like to see if doing a k-mean method it cluster the digits

# In[25]:

#I am going to use just the new dataset (without the labels), as I think it will perform better that way
features_newdataset = newdataset.iloc[:,1:]

# In[26]:

#I import my method
from sklearn.cluster import KMeans

# In[27]:

#I adjust the method
kmeans = KMeans(n_clusters=10)

# In[28]:

#training!
kmeans.fit(features_newdataset)

# In[29]:

kmeans.labels_

# In[30]:

labels_newdataset = newdataset.iloc[:,0]

# In[31]:

for i in range(10):
    labels_newdataset[kmeans.labels_ == i].hist()
    plt.show()

#  not at all what we expected. This is not a good method for this problem.

## Lets try to do a neural network

#I need to modify my set of labels; right now I have a vector with numbers 0 to 9,
#What I would need is a matrix of wide 10 so 0= (1,0,0,...,0), 1=(0,1,0,...,0)
labels_array = np.asarray(labels_newtrain)
labels_matrix = np.zeros((len(labels_newtrain),10))
#len(labels_newtrain)
#labels_matrix[0,0]

for i in range(0,len(labels_newtrain)):
    for j in labels_array[[i]]:
        labels_matrix[i,j] = 1.

labels_matrix

#Importing libraries
from keras.models import Sequential
from keras.layers import Dense

# create model
from keras.layers import Dropout
model = Sequential()
model.add(Dense(70, input_dim=784, activation='relu')) #784 is the number of features, pixels in my case
model.add(Dropout(0.2))
model.add(Dense(30, activation='relu'))  #2 hidden layers,50 and 20 newron respectively Relu is the activation function
model.add(Dropout(0.2)) ####Lets add a drop out, that I do not know how to use!! :D
model.add(Dense(10, activation='sigmoid')) #10 is the output layer, the dense mens that they are fully connect
#sigmoid is the last activation function

# Compile model
model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])

# Fit the model
model.fit(features_newtrain, labels_matrix, epochs=50, batch_size=64)

##Cool accuracy! In then training set
scores = model.evaluate(features_newtrain, labels_matrix)
print("\n%s: %.2f%%" % (model.metrics_names[1], scores[1]*100))

#Now, lets try to predict
predictionsNN = model.predict(features_newtest)

#Preparing the test labels in a comparable way, as we did before:
#I need to modify my set of labels; right now I have a vector with numbers 0 to 9,
#What I would need is a matrix of wide 10 so 0= (1,0,0,...,0), 1=(0,1,0,...,0)
labels_array_test = np.asarray(labels_newtest)
labels_matrix_test = np.zeros((len(labels_newtest),10))

for i in range(0,len(labels_newtest)):
    for j in labels_array_test[[i]]:
        labels_matrix_test[i,j] = 1.

##Cool accuracy in the test set!
scores = model.evaluate(features_newtest, labels_matrix_test)
print("\n%s: %.2f%%" % (model.metrics_names[1], scores[1]*100))

max_stars_repo_name: martabraojos/digit-recognition-MNIST

max_issues_repo_name: martabraojos/digit-recognition-MNIST

max_forks_repo_name: martabraojos/digit-recognition-MNIST


text:
import torch
import torch.nn as nn
#import numpy as np

####################################################################################################################################################

# Context-aware loss function (see SoccerNetLoss)

####################################################################################################################################################

class SegmentationLoss(torch.nn.Module):

    def __init__(self, K, hit_radius = 0.1, miss_radius = 0.9):

        super(SegmentationLoss,self).__init__()

        self.K = K
        self.hit_radius = float(hit_radius)
        self.miss_radius = float(miss_radius)
        self.weight = torch.tensor([0.00059148, 0.0011937, 0.0257837, 0.02792131, 0.29968935, 0.02532632, 0.09931299, 0.00722039, 0.04672078, 0.01889945, 0.24348754, 0.03107723, 0], dtype=torch.float).cuda()

    def forward(self, gt_label, pred_score):
        # pred_score shape: torch.Size([32, 48, 13]) # Batch_size, Chunk_size, num_classes

        pred_score = torch.clamp(pred_score, 10e-3, 1 - 10e-3)
        # loss = nn.CrossEntropyLoss(weight = self.weight)

#        sum_loss=0;
#        for i in range(pred_score.shape[1]):
#            if (i>(pred_score.shape[1]/3)) and (i<2*(pred_score.shape[1]/3)):
#                temp_loss=loss(pred_score[:,i,:],torch.max(gt_label[:,i,],1)[1])
#            else:
#                temp_loss=0
#            sum_loss=sum_loss+temp_loss

        loss = self.cross_entropy_no_softmax(gt_label, pred_score)
        return loss

    # def forward(self, gt_label, pred_score):
    #     print('gt',gt_label.shape)
    #     loss =nn.CrossEntropyLoss()
    #     print(pred_score.shape)

    #     return loss(gt_label,pred_score)

    def cross_entropy_no_softmax(self, gt_label, pred_score):
        loss_elementwise = -gt_label * torch.log(pred_score) - (1 - gt_label) * torch.log(1 - pred_score)
        loss_summed = torch.sum(loss_elementwise, dim=1)
        loss_summed = torch.sum(loss_summed, dim=0) * self.weight
        loss_summed = torch.sum(loss_summed)
        return loss_summed

####################################################################################################################################################

# Spotting loss

####################################################################################################################################################

class SpottingLoss(torch.nn.Module):

    def __init__(self, lambda_coord, lambda_noobj):
        super(SpottingLoss,self).__init__()

        self.lambda_coord = lambda_coord
        self.lambda_noobj = lambda_noobj

    def forward(self,y_true, y_pred):
        y_pred = self.permute_ypred_for_matching(y_true,y_pred)
        loss = torch.sum(y_true[:,:,0]*self.lambda_coord*torch.square(y_true[:,:,1]-y_pred[:,:,1])  +  y_true[:,:,0]*torch.square(y_true[:,:,0]-y_pred[:,:,0]) +  (1-y_true[:,:,0])*self.lambda_noobj*torch.square(y_true[:,:,0]-y_pred[:,:,0]) +  y_true[:,:,0]*torch.sum(torch.square(y_true[:,:,2:]-y_pred[:,:,2:]),axis=-1)) #-y_true[:,:,0]*torch.sum(y_true[:,:,2:]*torch.log(y_pred[:,:,2:]),axis=-1)
        return loss

    def permute_ypred_for_matching(self, y_true, y_pred):

        alpha = y_true[:,:,0]
        x = y_true[:,:,1]
        p = y_pred[:,:,1]
        nb_pred = x.shape[-1]

        D = torch.abs(x.unsqueeze(-1).repeat(1,1,nb_pred) - p.unsqueeze(-2).repeat(1,nb_pred,1))
        D1 = 1-D
        Permut = 0*D

        alpha_filter = alpha.unsqueeze(-1).repeat(1,1,nb_pred)

        v_filter = alpha_filter
        h_filter = 0*v_filter + 1
        D2 = v_filter * D1

        for i in range(nb_pred):
            D2 = v_filter * D2
            D2 = h_filter * D2
            A = torch.nn.functional.one_hot(torch.argmax(D2,axis=-1),nb_pred)
            B = v_filter * A * D2
            C = torch.nn.functional.one_hot(torch.argmax(B,axis=-2),nb_pred).permute(0, 2, 1)
            E = v_filter * A * C
            Permut = Permut + E
            v_filter = (1-torch.sum(Permut,axis=-1))*alpha
            v_filter = v_filter.unsqueeze(-1).repeat(1,1,nb_pred)
            h_filter = 1-torch.sum(Permut, axis=-2)
            h_filter = h_filter.unsqueeze(-2).repeat(1,nb_pred,1)

        v_filter = 1-alpha_filter
        D2 = v_filter * D1
        D2 = h_filter * D2

        for i in range(nb_pred):
            D2 = v_filter * D2
            D2 = h_filter * D2
            A = torch.nn.functional.one_hot(torch.argmax(D2,axis=-1),nb_pred)
            B = v_filter * A * D2
            C = torch.nn.functional.one_hot(torch.argmax(B,axis=-2),nb_pred).permute(0, 2, 1)
            E = v_filter * A * C
            Permut = Permut + E
            v_filter = (1-torch.sum(Permut,axis=-1))*(1-alpha) #here comes the change
            v_filter = v_filter.unsqueeze(-1).repeat(1,1,nb_pred)
            h_filter = 1-torch.sum(Permut, axis=-2)
            h_filter = h_filter.unsqueeze(-2).repeat(1,nb_pred,1)

        permutation = torch.argmax(Permut,axis=-1)
        permuted = torch.gather(y_pred, 1, permutation.unsqueeze(-1).repeat(1,1,y_true.shape[-1]))

        return permuted

if __name__ == "__main__":

    #K_V1 = torch.FloatTensor([[-20,-20,-40],[-10,-10,-20],[60,10,10],[90,20,20]])
    labels = torch.load("labels.pt")
    output_segmentation = torch.load("output_segmentation.pt")

    #criterion_segmentation = ContextAwareLoss(K=K_V1)

    print(labels)
    print(output_segmentation)

    print(output_segmentation.max())
    print(output_segmentation.min())
    loss = criterion_segmentation(labels,output_segmentation)
    print(loss)

    torch.save(labels, "labels.pt")
    torch.save(output_segmentation, "output_segmentation.pt")

"""
def SymmetricLoss(pred_score, gt_label, K, hit_radius = 0.1, miss_radius = 0.9, coef_0 = 1., coef_1 = 1.):

    zeros = torch.zeros(pred_score.size()).to(pred_score.device).type(torch.float)
    pred_score = 1.-pred_score

    case1 = DownStep(gt_label, self.K[0]) * torch.max(zeros, - torch.log(pred_score) + torch.log(zeros + miss_radius))
    case2 = Interval(gt_label, self.K[0], self.K[1]) * torch.max(zeros, - torch.log(pred_score + (1.-pred_score)*(PartialIdentity(gt_label,self.K[0],self.K[1])-self.K[0])/(self.K[1]-self.K[0])) + torch.log(zeros + miss_radius))
    #note cases 1 and 2 could be regrouped with DownStep(gt_label, self.K[1]) * formula loss of case2
    case3 = Interval(gt_label, self.K[1], 0.) * torch.max(zeros, - torch.log(1.-pred_score + pred_score*(PartialIdentity(gt_label,self.K[1],0.)-0.)/(self.K[1]-0.)) + torch.log(zeros + 1.-hit_radius))
    case4 = Interval(gt_label, 0., self.K[2]) * torch.max(zeros, - torch.log(1.-pred_score + pred_score*(PartialIdentity(gt_label,0.,self.K[2])-0.)/(self.K[2]-0.)) + torch.log(zeros + 1.-hit_radius))
    case5 = Interval(gt_label, self.K[2], self.K[3]) * torch.max(zeros, - torch.log(pred_score + (1.-pred_score)*(PartialIdentity(gt_label,self.K[2],self.K[3])-self.K[3])/(self.K[2]-self.K[3])) + torch.log(zeros + miss_radius))
    case6 = UpStep(gt_label, self.K[3]) * torch.max(zeros, - torch.log(pred_score) + torch.log(zeros + miss_radius))
    #note cases 5 and 6 could be regrouped

    L = coef_0*case1 + coef_0*case2 + coef_1*case3 + coef_1*case4 + coef_0*case5 + coef_0*case6

    return torch.sum(L)

def MirrorLoss(pred_score, gt_label, K, hit_radius = 0.1, miss_radius = 0.9, coef_0 = 1., coef_1 = 1.):

    zeros = torch.zeros(pred_score.size()).to(pred_score.device).type(torch.float)
    pred_score = 1.-pred_score

    case1 = DownStep(gt_label, self.K[0]) * torch.max(zeros, - torch.log(pred_score) + torch.log(zeros + miss_radius))
    case2 = Interval(gt_label, self.K[0], self.K[1]) * torch.max(zeros, - torch.log(pred_score + (1.-pred_score)*(PartialIdentity(gt_label,self.K[0],self.K[1])-self.K[0])/(self.K[1]-self.K[0])) + torch.log(zeros + miss_radius))
    #note cases 1 and 2 could be regrouped with DownStep(gt_label, self.K[1]) * formula loss of case2
    case3 = Interval(gt_label, self.K[1], 0.) * torch.max(zeros, - torch.log(1.-pred_score + pred_score*(PartialIdentity(gt_label,self.K[1],0.)-0.)/(self.K[1]-0.)) + torch.log(zeros + 1.-hit_radius))
    case4 = Interval(gt_label, 0., self.K[2]) * zeros
    case5 = Interval(gt_label, self.K[2], self.K[3]) * torch.max(zeros, - torch.log(pred_score + (1.-pred_score)*(PartialIdentity(gt_label,self.K[2],self.K[3])-self.K[3])/(self.K[2]-self.K[3])) + torch.log(zeros + miss_radius))
    case6 = UpStep(gt_label, self.K[3]) * torch.max(zeros, - torch.log(pred_score) + torch.log(zeros + miss_radius))
    #note cases 5 and 6 could be regrouped

    L = coef_0*case1 + coef_0*case2 + coef_1*case3 + coef_1*case4 + coef_0*case5 + coef_0*case6

    return torch.sum(L)
"""

#Code for testing

#import numpy as np

#np.random.seed(25)
#
#y_true = np.random.randint(-20.,20.,size=100)
#y_pred = np.random.uniform(0.,1.,size=100)
#
#params = np.array([-10.,-5.,5.,10.])
#hit_radius = 0.1
#miss_radius = 0.9
#
#L1 = SoccerNetLoss(torch.Tensor(y_pred), torch.Tensor(y_true), params, hit_radius, miss_radius, 1., 1.)
#L2 = SymmetricLoss(torch.Tensor(y_pred), torch.Tensor(y_true), params, hit_radius, miss_radius, 1., 1.)
#L3 = MirrorLoss(torch.Tensor(y_pred), torch.Tensor(y_true), params, hit_radius, miss_radius, 1., 1.)
#
#print(L1)
#print(L2)
#print(L3)

#For testing

#import numpy as np
#
#nb_predict = 7
#batch_size = 2
#num_classes = 3
#
#x = np.random.uniform(0,1,(batch_size, nb_predict, 2+num_classes))
#x = np.around(x,2)
#x[...,0] = np.round(x[...,0])
#x[...,2:] = 0
#X = torch.tensor(x)
#
#
#p = np.random.uniform(0,1,(batch_size, nb_predict, 2+num_classes))
#p = np.around(p,2)
#p[...,0] = 0
#p[...,2:] = 0
#P = torch.tensor(p)
#
#
#print(x)
#print(p)
#permutation = torch_permute_ypred_for_matching(X,P)
#
#print(permutation)


text:
"""Tests for clone_using."""
import numpy as np
import pytest
import xarray as xr

import janitor  # noqa: F401

@pytest.mark.xarray
def test_successful_cloning_coords(da):
    """Test that clone_using coordinates works correctly."""

    # with copying coords
    new_da: xr.DataArray = da.clone_using(np.random.randn(*da.data.shape))

    with pytest.raises(AssertionError):
        np.testing.assert_equal(new_da.data, da.data)

    assert all(
        (
            new_coord == old_coord
            for new_coord, old_coord in zip(new_da.coords, da.coords)
        )
    )
    assert new_da.dims == da.dims

@pytest.mark.xarray
def test_successful_cloning_no_coords(da):
    """Test that cloning works without coordinates."""

    new_da: xr.DataArray = da.clone_using(
        np.random.randn(*da.data.shape), use_coords=False
    )

    with pytest.raises(AssertionError):
        np.testing.assert_equal(new_da.data, da.data)

    assert new_da.dims == da.dims

@pytest.mark.xarray
def test_metadata_cloning(da):
    """Test that metadata gets cloned over."""
    new_da: xr.DataArray = da.clone_using(
        np.random.randn(*da.data.shape), use_attrs=True, new_name="new_name"
    )

    assert new_da.name != da.name
    assert new_da.attrs == da.attrs

@pytest.mark.xarray
def test_no_coords_errors(da: xr.DataArray):
    """Test that errors are raised when dims do not match."""
    # number of dims should match
    with pytest.raises(ValueError):
        da.clone_using(np.random.randn(10, 10, 10), use_coords=False)

    # shape of each axis does not need to match
    da.clone_using(np.random.randn(10, 10), use_coords=False)

@pytest.mark.xarray
def test_coords_errors(da: xr.DataArray):
    # number of dims should match
    with pytest.raises(ValueError):
        da.clone_using(np.random.randn(10, 10, 10), use_coords=False)

    # shape of each axis must match when using coords
    with pytest.raises(ValueError):
        da.clone_using(np.random.randn(10, 10), use_coords=True)

max_stars_repo_name: richardqiu/pyjanitor

max_issues_repo_name: richardqiu/pyjanitor

max_forks_repo_name: richardqiu/pyjanitor