text:
"""Functions for working with Protein Structure Graphs."""
# %%
# Graphein
# Author: Arian Jamasb <arian@jamasb.io>, Eric Ma, Charlie Harris
# License: MIT
# Project Website: https://github.com/a-r-j/graphein
# Code Repository: https://github.com/a-r-j/graphein
from __future__ import annotations

import logging
import traceback
from functools import partial
from typing import Any, Callable, Dict, List, Optional, Tuple, Union

import networkx as nx
import numpy as np
import pandas as pd
from Bio.PDB.Polypeptide import three_to_one
from biopandas.pdb import PandasPdb
from tqdm.contrib.concurrent import process_map

from graphein.protein.config import (
    DSSPConfig,
    GetContactsConfig,
    ProteinGraphConfig,
)
from graphein.protein.edges.distance import compute_distmat
from graphein.protein.resi_atoms import BACKBONE_ATOMS, RESI_THREE_TO_1
from graphein.protein.subgraphs import extract_subgraph_from_chains
from graphein.protein.utils import (
    ProteinGraphConfigurationError,
    compute_rgroup_dataframe,
    filter_dataframe,
    get_protein_name_from_filename,
    three_to_one_with_mods,
)
from graphein.utils.utils import (
    annotate_edge_metadata,
    annotate_graph_metadata,
    annotate_node_metadata,
    compute_edges,
)

logging.basicConfig(level="DEBUG")
log = logging.getLogger(__name__)

def read_pdb_to_dataframe(
    pdb_path: Optional[str] = None,
    pdb_code: Optional[str] = None,
    verbose: bool = False,
    granularity: str = "CA",
) -> pd.DataFrame:
    """
    Reads PDB file to ``PandasPDB`` object.

    Returns ``atomic_df``, which is a dataframe enumerating all atoms and their cartesian coordinates in 3D space. Also
    contains associated metadata from the PDB file.

    :param pdb_path: path to PDB file. Defaults to None.
    :type pdb_path: str, optional
    :param pdb_code: 4-character PDB accession. Defaults to None.
    :type pdb_code: str, optional
    :param verbose: print dataframe?
    :type verbose: bool
    :param granularity: Specifies granularity of dataframe. See :class:`~graphein.protein.config.ProteinGraphConfig` for further
        details.
    :type granularity: str
    :returns: ``pd.DataFrame`` containing protein structure
    :rtype: pd.DataFrame
    """

    assert (pdb_code and not pdb_path) or (not pdb_code and pdb_path), (
        "Either a PDB ID or a path to a local PDB file"
        " must be specified to read a PDB"
    )

    atomic_df = (
        PandasPdb().read_pdb(pdb_path)
        if pdb_path is not None
        else PandasPdb().fetch_pdb(pdb_code)
    )

    # Assign Node IDs to dataframes
    atomic_df.df["ATOM"]["node_id"] = (
        atomic_df.df["ATOM"]["chain_id"].apply(str)
        + ":"
        + atomic_df.df["ATOM"]["residue_name"]
        + ":"
        + atomic_df.df["ATOM"]["residue_number"].apply(str)
    )
    if granularity == "atom":
        atomic_df.df["ATOM"]["node_id"] = (
            atomic_df.df["ATOM"]["node_id"]
            + ":"
            + atomic_df.df["ATOM"]["atom_name"]
        )
    if verbose:
        print(atomic_df)
    return atomic_df

def deprotonate_structure(df: pd.DataFrame) -> pd.DataFrame:
    """Remove protons from PDB dataframe.

    :param df: Atomic dataframe.
    :type df: pd.DataFrame
    :returns: Atomic dataframe with all ``atom_name == "H"`` removed.
    :rtype: pd.DataFrame
    """
    log.debug(
        "Deprotonating protein. This removes H atoms from the pdb_df dataframe"
    )
    return filter_dataframe(
        df, by_column="atom_name", list_of_values=["H"], boolean=False
    )

def convert_structure_to_centroids(df: pd.DataFrame) -> pd.DataFrame:
    """Overwrite existing ``(x, y, z)`` coordinates with centroids of the amino acids.

    :param df: Pandas Dataframe protein structure to convert into a dataframe of centroid positions.
    :type df: pd.DataFrame
    :return: pd.DataFrame with atoms/residues positions converted into centroid positions.
    :rtype: pd.DataFrame
    """
    log.debug(
        "Converting dataframe to centroids. This averages XYZ coords of the atoms in a residue"
    )

    centroids = calculate_centroid_positions(df)
    df = df.loc[df["atom_name"] == "CA"].reset_index(drop=True)
    df["x_coord"] = centroids["x_coord"]
    df["y_coord"] = centroids["y_coord"]
    df["z_coord"] = centroids["z_coord"]

    return df

def subset_structure_to_atom_type(
    df: pd.DataFrame, granularity: str
) -> pd.DataFrame:
    """
    Return a subset of atomic dataframe that contains only certain atom names.

    :param df: Protein Structure dataframe to subset.
    :type df: pd.DataFrame
    :returns: Subsetted protein structure dataframe.
    :rtype: pd.DataFrame
    """
    return filter_dataframe(
        df, by_column="atom_name", list_of_values=[granularity], boolean=True
    )

def remove_insertions(df: pd.DataFrame, keep: str = "first") -> pd.DataFrame:
    """
    This function removes insertions from PDB dataframes.

    :param df: Protein Structure dataframe to remove insertions from.
    :type df: pd.DataFrame
    :param keep: Specifies which insertion to keep. Options are ``"first"`` or ``"last"``.
        Default is ``"first"``
    :type keep: str
    :return: Protein structure dataframe with insertions removed
    :rtype: pd.DataFrame
    """
    # Catches unnamed insertions
    duplicates = df.duplicated(
        subset=["chain_id", "residue_number", "atom_name"], keep=keep
    )
    df = df[~duplicates]

    # Catches explicit insertions
    df = filter_dataframe(
        df, by_column="insertion", list_of_values=[""], boolean=True
    )

    # Remove alt_locs
    df = filter_dataframe(
        df, by_column="alt_loc", list_of_values=["", "A"], boolean=True
    )

    return df

def filter_hetatms(
    df: pd.DataFrame, keep_hets: List[str]
) -> List[pd.DataFrame]:
    """Return hetatms of interest.

    :param df: Protein Structure dataframe to filter hetatoms from.
    :type df: pd.DataFrame
    :param keep_hets: List of hetero atom names to keep.
    :returns: Protein structure dataframe with heteroatoms removed
    :rtype pd.DataFrame
    """
    return [df.loc[df["residue_name"] == hetatm] for hetatm in keep_hets]

def process_dataframe(
    protein_df: pd.DataFrame,
    atom_df_processing_funcs: Optional[List[Callable]] = None,
    hetatom_df_processing_funcs: Optional[List[Callable]] = None,
    granularity: str = "centroids",
    chain_selection: str = "all",
    insertions: bool = False,
    deprotonate: bool = True,
    keep_hets: List[str] = [],
    verbose: bool = False,
) -> pd.DataFrame:
    """
    Process ATOM and HETATM dataframes to produce singular dataframe used for graph construction.

    :param protein_df: Dataframe to process.
        Should be the object returned from :func:`~graphein.protein.graphs.read_pdb_to_dataframe`.
    :type protein_df: pd.DataFrame
    :param atom_df_processing_funcs: List of functions to process dataframe. These must take in a dataframe and return a
        dataframe. Defaults to None.
    :type atom_df_processing_funcs: List[Callable], optional
    :param hetatom_df_processing_funcs: List of functions to process the hetatom dataframe. These must take in a dataframe and return a dataframe
    :type hetatom_df_processing_funcs: List[Callable], optional
    :param granularity: The level of granularity for the graph. This determines the node definition.
        Acceptable values include: ``"centroids"``, ``"atoms"``,
        any of the atom_names in the PDB file (e.g. ``"CA"``, ``"CB"``, ``"OG"``, etc.).
        See: :const:`~graphein.protein.config.GRAPH_ATOMS` and :const:`~graphein.protein.config.GRANULARITY_OPTS`.
    :type granularity: str
    :param insertions: Whether or not to keep insertions.
    :param insertions: bool
    :param deprotonate: Whether or not to remove hydrogen atoms (i.e. deprotonation).
    :type deprotonate: bool
    :param keep_hets: Hetatoms to keep. Defaults to an empty list.
        To keep a hetatom, pass it inside a list of hetatom names to keep.
    :type keep_hets: List[str]
    :param verbose: Verbosity level.
    :type verbose: bool
    :param chain_selection: Which protein chain to select. Defaults to ``"all"``. Eg can use ``"ACF"``
        to select 3 chains (``A``, ``C`` & ``F``)
    :type chain_selection: str
    :return: A protein dataframe that can be consumed by
        other graph construction functions.
    :rtype: pd.DataFrame
    """
    # TODO: Need to properly define what "granularity" is supposed to do.
    atoms = protein_df.df["ATOM"]
    hetatms = protein_df.df["HETATM"]

    # This block enables processing via a list of supplied functions operating on the atom and hetatom dataframes
    # If these are provided, the dataframe returned will be computed only from these and the default workflow
    # below this block will not execute.
    if atom_df_processing_funcs is not None:
        for func in atom_df_processing_funcs:
            atoms = func(atoms)
        if hetatom_df_processing_funcs is None:
            return atoms

    if hetatom_df_processing_funcs is not None:
        for func in hetatom_df_processing_funcs:
            hetatms = func(hetatms)
        return pd.concat([atoms, hetatms])

    # Deprotonate structure by removing H atoms
    if deprotonate:
        atoms = deprotonate_structure(atoms)

    # Restrict DF to desired granularity
    if granularity == "atom":
        pass
    elif granularity == "centroids":
        atoms = convert_structure_to_centroids(atoms)
    else:
        atoms = subset_structure_to_atom_type(atoms, granularity)

    protein_df = atoms

    if keep_hets:
        hetatms_to_keep = filter_hetatms(atoms, keep_hets)
        protein_df = pd.concat([atoms, hetatms_to_keep])

    # Remove alt_loc residues
    if not insertions:
        protein_df = remove_insertions(protein_df)

    # perform chain selection
    protein_df = select_chains(
        protein_df, chain_selection=chain_selection, verbose=verbose
    )

    log.debug(f"Detected {len(protein_df)} total nodes")

    return protein_df

def assign_node_id_to_dataframe(
    protein_df: pd.DataFrame, granularity: str
) -> pd.DataFrame:
    """
    Assigns the node ID back to the ``pdb_df`` dataframe

    :param protein_df: Structure Dataframe
    :type protein_df: pd.DataFrame
    :param granularity: Granularity of graph. Atom-level,
        residue (e.g. ``CA``) or ``centroids``.
        See: :const:`~graphein.protein.config.GRAPH_ATOMS`
        and :const:`~graphein.protein.config.GRANULARITY_OPTS`.
    :type granularity: str
    :return: Returns dataframe with added ``node_ids``
    :rtype: pd.DataFrame
    """
    protein_df["node_id"] = (
        protein_df["chain_id"].apply(str)
        + ":"
        + protein_df["residue_name"]
        + ":"
        + protein_df["residue_number"].apply(str)
    )
    if granularity == "atom":
        protein_df[
            "node_id"
        ] = f'{protein_df["node_id"]}:{protein_df["atom_name"]}'

def select_chains(
    protein_df: pd.DataFrame, chain_selection: str, verbose: bool = False
) -> pd.DataFrame:
    """
    Extracts relevant chains from ``protein_df``.

    :param protein_df: pandas dataframe of PDB subsetted to relevant atoms
        (``CA``, ``CB``).
    :type protein_df: pd.DataFrame
    :param chain_selection: Specifies chains that should be extracted from
        the larger complexed structure.
    :type chain_selection: str
    :param verbose: Print dataframe?
    :type verbose: bool
    :return: Protein structure dataframe containing only entries in the
        chain selection.
    :rtype: pd.DataFrame
    """
    if chain_selection != "all":
        protein_df = filter_dataframe(
            protein_df,
            by_column="chain_id",
            list_of_values=list(chain_selection),
            boolean=True,
        )

    return protein_df

def initialise_graph_with_metadata(
    protein_df: pd.DataFrame,
    raw_pdb_df: pd.DataFrame,
    granularity: str,
    name: Optional[str] = None,
    pdb_code: Optional[str] = None,
    pdb_path: Optional[str] = None,
) -> nx.Graph:
    """
    Initializes the nx Graph object with initial metadata.

    :param protein_df: Processed Dataframe of protein structure.
    :type protein_df: pd.DataFrame
    :param raw_pdb_df: Unprocessed dataframe of protein structure for comparison and traceability downstream.
    :type raw_pdb_df: pd.DataFrame
    :param granularity: Granularity of the graph (eg ``"atom"``, ``"CA"``, ``"CB"`` etc or ``"centroid"``).
        See: :const:`~graphein.protein.config.GRAPH_ATOMS` and :const:`~graphein.protein.config.GRANULARITY_OPTS`.
    :type granularity: str
    :param name: specified given name for the graph. If None, the PDB code or the file name will be used to name the graph.
    :type name: Optional[str], defaults to ``None``
    :param pdb_code: PDB ID / Accession code, if the PDB is available on the PDB database.
    :type pdb_code: Optional[str], defaults to ``None``
    :param pdb_path: path to local PDB file, if constructing a graph from a local file.
    :type pdb_path: Optional[str], defaults to ``None``
    :return: Returns initial protein structure graph with metadata.
    :rtype: nx.Graph
    """

    assert (pdb_code and not pdb_path) or (not pdb_code and pdb_path), (
        "Either a PDB ID or a path to a local PDB file"
        " must be specified to read a PDB"
    )

    # Get name for graph if no name was provided
    if not name:
        if pdb_path:
            name = get_protein_name_from_filename(pdb_path)
        else:
            name = pdb_code

    G = nx.Graph(
        name=name,
        pdb_code=pdb_code,
        pdb_path=pdb_path,
        chain_ids=list(protein_df["chain_id"].unique()),
        pdb_df=protein_df,
        raw_pdb_df=raw_pdb_df,
        rgroup_df=compute_rgroup_dataframe(remove_insertions(raw_pdb_df)),
        coords=np.asarray(protein_df[["x_coord", "y_coord", "z_coord"]]),
    )

    # Create graph and assign intrinsic graph-level metadata
    G.graph["node_type"] = granularity

    # Add Sequences to graph metadata
    for c in G.graph["chain_ids"]:
        G.graph[f"sequence_{c}"] = (
            protein_df.loc[protein_df["chain_id"] == c]["residue_name"]
            .apply(three_to_one_with_mods)
            .str.cat()
        )
    return G

def add_nodes_to_graph(
    G: nx.Graph,
    protein_df: Optional[pd.DataFrame] = None,
    verbose: bool = False,
) -> nx.Graph:
    """Add nodes into protein graph.

    :param G: ``nx.Graph`` with metadata to populate with nodes.
    :type G: nx.Graph
    :protein_df: DataFrame of protein structure containing nodes & initial node metadata to add to the graph.
    :type protein_df: pd.DataFrame, optional
    :param verbose: Controls verbosity of this step.
    :type verbose: bool
    :returns: nx.Graph with nodes added.
    :rtype: nx.Graph
    """

    # If no protein dataframe is supplied, use the one stored in the Graph object
    if protein_df is None:
        protein_df = G.graph["pdb_df"]
    # Assign intrinsic node attributes
    chain_id = protein_df["chain_id"].apply(str)
    residue_name = protein_df["residue_name"]
    residue_number = protein_df["residue_number"]  # .apply(str)
    coords = np.asarray(protein_df[["x_coord", "y_coord", "z_coord"]])
    b_factor = protein_df["b_factor"]
    atom_type = protein_df["atom_name"]
    nodes = protein_df["node_id"]
    element_symbol = protein_df["element_symbol"]
    G.add_nodes_from(nodes)

    # Set intrinsic node attributes
    nx.set_node_attributes(G, dict(zip(nodes, chain_id)), "chain_id")
    nx.set_node_attributes(G, dict(zip(nodes, residue_name)), "residue_name")
    nx.set_node_attributes(
        G, dict(zip(nodes, residue_number)), "residue_number"
    )
    nx.set_node_attributes(G, dict(zip(nodes, atom_type)), "atom_type")
    nx.set_node_attributes(
        G, dict(zip(nodes, element_symbol)), "element_symbol"
    )
    nx.set_node_attributes(G, dict(zip(nodes, coords)), "coords")
    nx.set_node_attributes(G, dict(zip(nodes, b_factor)), "b_factor")

    # TODO: include charge, line_idx for traceability?
    if verbose:
        print(nx.info(G))
        print(G.nodes())

    return G

def calculate_centroid_positions(
    atoms: pd.DataFrame, verbose: bool = False
) -> pd.DataFrame:
    """
    Calculates position of sidechain centroids.

    :param atoms: ATOM df of protein structure.
    :type atoms: pd.DataFrame
    :param verbose: bool controlling verbosity.
    :type verbose: bool
    :return: centroids (df).
    :rtype: pd.DataFrame
    """
    centroids = (
        atoms.groupby("residue_number")
        .mean()[["x_coord", "y_coord", "z_coord"]]
        .reset_index()
    )
    if verbose:
        print(f"Calculated {len(centroids)} centroid nodes")
    log.debug(f"Calculated {len(centroids)} centroid nodes")
    return centroids

def compute_edges(
    G: nx.Graph,
    funcs: List[Callable],
    get_contacts_config: Optional[GetContactsConfig] = None,
) -> nx.Graph:
    """
    Computes edges for the protein structure graph. Will compute a pairwise
    distance matrix between nodes which is
    added to the graph metadata to facilitate some edge computations.

    :param G: nx.Graph with nodes to add edges to.
    :type G: nx.Graph
    :param funcs: List of edge construction functions.
    :type funcs: List[Callable]
    :param get_contacts_config: Config object for ``GetContacts`` if
        intramolecular edges are being used.
    :type get_contacts_config: graphein.protein.config.GetContactsConfig
    :return: Graph with added edges.
    :rtype: nx.Graph
    """
    # This control flow prevents unnecessary computation of the distance matrices
    if "config" in G.graph:
        if G.graph["config"].granularity == "atom":
            G.graph["atomic_dist_mat"] = compute_distmat(G.graph["raw_pdb_df"])
        else:
            G.graph["dist_mat"] = compute_distmat(G.graph["pdb_df"])

    for func in funcs:
        func(G)

    return G

def construct_graph(
    config: Optional[ProteinGraphConfig] = None,
    name: Optional[str] = None,
    pdb_path: Optional[str] = None,
    pdb_code: Optional[str] = None,
    chain_selection: str = "all",
    df_processing_funcs: Optional[List[Callable]] = None,
    edge_construction_funcs: Optional[List[Callable]] = None,
    edge_annotation_funcs: Optional[List[Callable]] = None,
    node_annotation_funcs: Optional[List[Callable]] = None,
    graph_annotation_funcs: Optional[List[Callable]] = None,
) -> nx.Graph:
    """
    Constructs protein structure graph from a ``pdb_code`` or ``pdb_path``.

    Users can provide a :class:`~graphein.protein.config.ProteinGraphConfig`
    object to specify construction parameters.

    However, config parameters can be overridden by passing arguments directly to the function.

    :param config: :class:`~graphein.protein.config.ProteinGraphConfig` object. If None, defaults to config in ``graphein.protein.config``.
    :type config: graphein.protein.config.ProteinGraphConfig, optional
    :param name: an optional given name for the graph. the PDB ID or PDB file name will be used if not specified.
    :type name: str, optional
    :param pdb_path: Path to ``pdb_file`` when constructing a graph from a local pdb file. Default is ``None``.
    :type pdb_path: Optional[str], defaults to ``None``
    :param pdb_code: A 4-character PDB ID / accession to be used to construct the graph, if available. Default is ``None``.
    :type pdb_code: Optional[str], defaults to ``None``
    :param chain_selection: String of polypeptide chains to include in graph. E.g ``"ABDF"`` or ``"all"``. Default is ``"all"``.
    :type chain_selection: str
    :param df_processing_funcs: List of dataframe processing functions. Default is ``None``.
    :type df_processing_funcs: List[Callable], optional
    :param edge_construction_funcs: List of edge construction functions. Default is ``None``.
    :type edge_construction_funcs: List[Callable], optional
    :param edge_annotation_funcs: List of edge annotation functions. Default is ``None``.
    :type edge_annotation_funcs: List[Callable], optional
    :param node_annotation_funcs: List of node annotation functions. Default is ``None``.
    :type node_annotation_funcs: List[Callable], optional
    :param graph_annotation_funcs: List of graph annotation function. Default is ``None``.
    :type graph_annotation_funcs: List[Callable]
    :return: Protein Structure Graph
    :type: nx.Graph
    """

    assert (pdb_code and not pdb_path) or (not pdb_code and pdb_path), (
        "Either a PDB ID or a path to a local PDB file"
        " must be specified to construct a graph"
    )

    # If no config is provided, use default
    if config is None:
        config = ProteinGraphConfig()

    # If config params are provided, overwrite them
    config.protein_df_processing_functions = (
        df_processing_funcs
        if config.protein_df_processing_functions is None
        else config.protein_df_processing_functions
    )
    config.edge_construction_functions = (
        edge_construction_funcs
        if config.edge_construction_functions is None
        else config.edge_construction_functions
    )
    config.node_metadata_functions = (
        node_annotation_funcs
        if config.node_metadata_functions is None
        else config.node_metadata_functions
    )
    config.graph_metadata_functions = (
        graph_annotation_funcs
        if config.graph_metadata_functions is None
        else config.graph_metadata_functions
    )
    config.edge_metadata_functions = (
        edge_annotation_funcs
        if config.edge_metadata_functions is None
        else config.edge_metadata_functions
    )

    raw_df = read_pdb_to_dataframe(
        pdb_path,
        pdb_code,
        verbose=config.verbose,
        granularity=config.granularity,
    )
    protein_df = process_dataframe(
        raw_df,
        chain_selection=chain_selection,
        granularity=config.granularity,
        insertions=config.insertions,
    )

    # Initialise graph with metadata
    g = initialise_graph_with_metadata(
        protein_df=protein_df,
        raw_pdb_df=raw_df.df["ATOM"],
        name=name,
        pdb_code=pdb_code,
        pdb_path=pdb_path,
        granularity=config.granularity,
    )
    # Add nodes to graph
    g = add_nodes_to_graph(g)

    # Add config to graph
    g.graph["config"] = config

    # Annotate additional node metadata
    if config.node_metadata_functions is not None:
        g = annotate_node_metadata(g, config.node_metadata_functions)

    # Compute graph edges
    g = compute_edges(
        g,
        funcs=config.edge_construction_functions,
        get_contacts_config=None,
    )

    # Annotate additional graph metadata
    if config.graph_metadata_functions is not None:
        g = annotate_graph_metadata(g, config.graph_metadata_functions)

    # Annotate additional edge metadata
    if config.edge_metadata_functions is not None:
        g = annotate_edge_metadata(g, config.edge_metadata_functions)

    return g

def _mp_graph_constructor(
    args: Tuple[str, str], use_pdb_code: bool, config: ProteinGraphConfig
) -> nx.Graph:
    """
    Protein graph constructor for use in multiprocessing several protein structure graphs.

    :param args: Tuple of pdb code/path and the chain selection for that PDB
    :type args: Tuple[str, str]
    :param use_pdb_code: Whether or not we are using pdb codes or paths
    :type use_pdb_code: bool
    :param config: Protein structure graph construction config
    :type config: ProteinGraphConfig
    :return: Protein structure graph
    :rtype: nx.Graph
    """
    log.info(f"Constructing graph for: {args[0]}. Chain selection: {args[1]}")
    func = partial(construct_graph, config=config)
    try:
        return (
            func(pdb_code=args[0], chain_selection=args[1])
            if use_pdb_code
            else func(pdb_path=args[0], chain_selection=args[1])
        )

    except Exception as ex:
        log.info(
            f"Graph construction error (PDB={args[0]})! {traceback.format_exc()}"
        )
        log.info(ex)
        return None

def construct_graphs_mp(
    pdb_code_it: Optional[List[str]] = None,
    pdb_path_it: Optional[List[str]] = None,
    chain_selections: Optional[list[str]] = None,
    config: ProteinGraphConfig = ProteinGraphConfig(),
    num_cores: int = 16,
    return_dict: bool = True,
    out_path: Optional[str] = None,
) -> Union[List[nx.Graph], Dict[str, nx.Graph]]:
    """
    Constructs protein graphs for a list of pdb codes or pdb paths using multiprocessing.

    :param pdb_code_it: List of pdb codes to use for protein graph construction
    :type pdb_code_it: Optional[List[str]], defaults to ``None``
    :param pdb_path_it: List of paths to PDB files to use for protein graph construction
    :type pdb_path_it: Optional[List[str]], defaults to ``None``
    :param chain_selections: List of chains to select from the protein structures (e.g. ``["ABC", "A", "L", "CD"...]``)
    :type chain_selections: Optional[List[str]], defaults to ``None``
    :param config: ProteinGraphConfig to use.
    :type config: graphein.protein.config.ProteinGraphConfig, defaults to default config params
    :param num_cores: Number of cores to use for multiprocessing. The more the merrier
    :type num_cores: int, defaults to ``16``
    :param return_dict: Whether or not to return a dictionary (indexed by pdb codes/paths) or a list of graphs.
    :type return_dict: bool, default to ``True``
    :param out_path: Path to save the graphs to. If None, graphs are not saved.
    :type out_path: Optional[str], defaults to ``None``
    :return: Iterable of protein graphs. None values indicate there was a problem in constructing the graph for this particular pdb
    :rtype: Union[List[nx.Graph], Dict[str, nx.Graph]]
    """
    assert (
        pdb_code_it is not None or pdb_path_it is not None
    ), "Iterable of pdb codes OR pdb paths required."

    if pdb_code_it is not None:
        pdbs = pdb_code_it
        use_pdb_code = True

    if pdb_path_it is not None:
        pdbs = pdb_path_it
        use_pdb_code = False

    if chain_selections is None:
        chain_selections = ["all"] * len(pdbs)

    constructor = partial(
        _mp_graph_constructor, use_pdb_code=use_pdb_code, config=config
    )

    graphs = list(
        process_map(
            constructor,
            [(pdb, chain_selections[i]) for i, pdb in enumerate(pdbs)],
            max_workers=num_cores,
        )
    )
    if out_path is not None:
        [
            nx.write_gpickle(
                g, str(f"{out_path}/" + f"{g.graph['name']}.pickle")
            )
            for g in graphs
        ]

    if return_dict:
        graphs = {pdb: graphs[i] for i, pdb in enumerate(pdbs)}

    return graphs

def compute_chain_graph(
    g: nx.Graph,
    chain_list: Optional[List[str]] = None,
    remove_self_loops: bool = False,
    return_weighted_graph: bool = False,
) -> Union[nx.Graph, nx.MultiGraph]:
    """Computes a chain-level graph from a protein structure graph.

    This graph features nodes as individual chains in a complex and edges as
    the interactions between constituent nodes in each chain. You have the
    option of returning an unweighted graph (multigraph,
    ``return_weighted_graph=False``) or a weighted graph
    (``return_weighted_graph=True``). The difference between these is the
    unweighted graph features and edge for each interaction between chains
    (ie the number of edges will be equal to the number of edges in the input
    protein structure graph), while the weighted graph sums these interactions
    to a single edge between chains with the counts stored as features.

    :param g: A protein structure graph to compute the chain graph of.
    :type g: nx.Graph
    :param chain_list: A list of chains to extract from the input graph.
        If ``None``, all chains will be used. This is provided as input to
        ``extract_subgraph_from_chains``. Default is ``None``.
    :type chain_list: Optional[List[str]]
    :param remove_self_loops: Whether to remove self-loops from the graph.
        Default is False.
    :type remove_self_loops: bool
    :return: A chain-level graph.
    :rtype: Union[nx.Graph, nx.MultiGraph]
    """
    # If we are extracting specific chains, do it here.
    if chain_list is not None:
        g = extract_subgraph_from_chains(g, chain_list)

    # Initialise new graph with Metadata
    h = nx.MultiGraph()
    h.graph = g.graph
    h.graph["node_type"] = "chain"

    # Set nodes
    nodes_per_chain = {chain: 0 for chain in g.graph["chain_ids"]}
    sequences = {chain: "" for chain in g.graph["chain_ids"]}
    for n, d in g.nodes(data=True):
        nodes_per_chain[d["chain_id"]] += 1
        sequences[d["chain_id"]] += RESI_THREE_TO_1[d["residue_name"]]

    h.add_nodes_from(g.graph["chain_ids"])

    for n, d in h.nodes(data=True):
        d["num_residues"] = nodes_per_chain[n]
        d["sequence"] = sequences[n]

    # Add edges
    for u, v, d in g.edges(data=True):
        h.add_edge(
            g.nodes[u]["chain_id"], g.nodes[v]["chain_id"], kind=d["kind"]
        )
    # Remove self-loops if necessary. Checks for equality between nodes in a given edge.
    if remove_self_loops:
        edges_to_remove: List[Tuple[str]] = [
            (u, v) for u, v in h.edges() if u == v
        ]
        h.remove_edges_from(edges_to_remove)

    # Compute a weighted graph if required.
    if return_weighted_graph:
        return compute_weighted_graph_from_multigraph(h)
    return h

def compute_weighted_graph_from_multigraph(g: nx.MultiGraph) -> nx.Graph:
    """Computes a weighted graph from a multigraph.

    This function is used to convert a multigraph to a weighted graph. The
    weights of the edges are the number of interactions between the nodes.

    :param g: A multigraph.
    :type g: nx.MultiGraph
    :return: A weighted graph.
    :rtype: nx.Graph
    """
    H = nx.Graph()
    H.graph = g.graph
    H.add_nodes_from(g.nodes(data=True))
    for u, v, d in g.edges(data=True):
        if H.has_edge(u, v):
            H[u][v]["weight"] += len(d["kind"])
            H[u][v]["kind"].update(d["kind"])
            for kind in list(d["kind"]):
                try:
                    H[u][v][kind] += 1
                except KeyError:
                    H[u][v][kind] = 1
        else:
            H.add_edge(u, v, weight=len(d["kind"]), kind=d["kind"])
            for kind in list(d["kind"]):
                H[u][v][kind] = 1
    return H

def number_groups_of_runs(list_of_values: List[Any]) -> List[str]:
    """Numbers groups of runs in a list of values.

    E.g. ``["A", "A", "B", "A", "A", "A", "B", "B"] ->
    ["A1", "A1", "B1", "A2", "A2", "A2", "B2", "B2"]``

    :param list_of_values: List of values to number.
    :type list_of_values: List[Any]
    :return: List of numbered values.
    :rtype: List[str]
    """
    df = pd.DataFrame({"val": list_of_values})
    df["idx"] = df["val"].shift() != df["val"]
    df["sum"] = df.groupby("val")["idx"].cumsum()
    return list(df["val"].astype(str) + df["sum"].astype(str))

def compute_secondary_structure_graph(
    g: nx.Graph,
    allowable_ss_elements: Optional[List[str]] = None,
    remove_non_ss: bool = True,
    remove_self_loops: bool = False,
    return_weighted_graph: bool = False,
) -> Union[nx.Graph, nx.MultiGraph]:
    """Computes a secondary structure graph from a protein structure graph.

    :param g: A protein structure graph to compute the secondary structure
        graph of.
    :type g: nx.Graph
    :param remove_non_ss: Whether to remove non-secondary structure nodes from
        the graph. These are denoted as ``"-"`` by DSSP. Default is True.
    :type remove_non_ss: bool
    :param remove_self_loops: Whether to remove self-loops from the graph.
        Default is ``False``.
    :type remove_self_loops: bool
    :param return_weighted_graph: Whether to return a weighted graph.
        Default is False.
    :type return_weighted_graph: bool
    :raises ProteinGraphConfigurationError: If the protein structure graph is
        not configured correctly with secondary structure assignments on all
        nodes.
    :return: A secondary structure graph.
    :rtype: Union[nx.Graph, nx.MultiGraph]
    """
    # Initialise list of secondary structure elements we use to build the graph
    ss_list: List[str] = []

    # Check nodes have secondary structure assignment & store them in list
    for _, d in g.nodes(data=True):
        if "ss" not in d.keys():
            raise ProteinGraphConfigurationError(
                "Secondary structure not defined for all nodes."
            )
        ss_list.append(d["ss"])

    # Number SS elements
    ss_list = pd.Series(number_groups_of_runs(ss_list))
    ss_list.index = list(g.nodes())

    # Remove unstructured elements if necessary
    if remove_non_ss:
        ss_list = ss_list[~ss_list.str.contains("-")]
    # Subset to only allowable SS elements if necessary
    if allowable_ss_elements:
        ss_list = ss_list[
            ss_list.str.contains("|".join(allowable_ss_elements))
        ]

    constituent_residues: Dict[str, List[str]] = ss_list.index.groupby(
        ss_list.values
    )
    constituent_residues = {
        k: list(v) for k, v in constituent_residues.items()
    }
    residue_counts: Dict[str, int] = ss_list.groupby(ss_list).count().to_dict()

    # Add Nodes from secondary structure list
    h = nx.MultiGraph()
    h.add_nodes_from(ss_list)
    nx.set_node_attributes(h, residue_counts, "residue_counts")
    nx.set_node_attributes(h, constituent_residues, "constituent_residues")
    # Assign ss
    for n, d in h.nodes(data=True):
        d["ss"] = n[0]

    # Add graph-level metadata
    h.graph = g.graph
    h.graph["node_type"] = "secondary_structure"

    # Iterate over edges in source graph and add SS-SS edges to new graph.
    for u, v, d in g.edges(data=True):
        try:
            h.add_edge(
                ss_list[u], ss_list[v], kind=d["kind"], source=f"{u}_{v}"
            )
        except KeyError as e:
            log.debug(
                f"Edge {u}-{v} not added to secondary structure graph. \
                Reason: {e} not in graph"
            )

    # Remove self-loops if necessary.
    # Checks for equality between nodes in a given edge.
    if remove_self_loops:
        edges_to_remove: List[Tuple[str]] = [
            (u, v) for u, v in h.edges() if u == v
        ]
        h.remove_edges_from(edges_to_remove)

    # Create weighted graph from h
    if return_weighted_graph:
        return compute_weighted_graph_from_multigraph(h)
    return h

max_stars_repo_name: avivko/graphein

max_issues_repo_name: avivko/graphein

max_forks_repo_name: avivko/graphein

reason: import numpy,import networkx


text:
import objax
from jax import vmap
import jax.numpy as np
from jax.scipy.linalg import cho_factor, cho_solve, block_diag
from .utils import scaled_squared_euclid_dist, softplus, softplus_inv, rotation_matrix
from warnings import warn

__all__ = ["Matern12","Matern32"]

class Kernel(objax.Module):
    """
    """

    def __call__(self, X, X2):
        return self.K(X, X2)

    def K(self, X, X2):
        raise NotImplementedError('kernel function not implemented')

    def measurement_model(self):
        raise NotImplementedError

    def inducing_precision(self):
        return None, None

    def kernel_to_state_space(self, R=None):
        raise NotImplementedError

    def spatial_conditional(self, R=None, predict=False):
        """
        """
        return None, None

class StationaryKernel(Kernel):
    """
    """

    def __init__(self,
                 variance=1.0,
                 lengthscale=1.0,
                 fix_variance=False,
                 fix_lengthscale=False):
        # check whether the parameters are to be optimised
        if fix_lengthscale:
            self.transformed_lengthscale = objax.StateVar(softplus_inv(np.array(lengthscale)))
        else:
            self.transformed_lengthscale = objax.TrainVar(softplus_inv(np.array(lengthscale)))
        if fix_variance:
            self.transformed_variance = objax.StateVar(softplus_inv(np.array(variance)))
        else:
            self.transformed_variance = objax.TrainVar(softplus_inv(np.array(variance)))

    @property
    def variance(self):
        return softplus(self.transformed_variance.value)

    @property
    def lengthscale(self):
        return softplus(self.transformed_lengthscale.value)

    def K(self, X, X2):
        r2 = scaled_squared_euclid_dist(X, X2, self.lengthscale)
        return self.K_r2(r2)

    def K_r2(self, r2):
        # Clipping around the (single) float precision which is ~1e-45.
        r = np.sqrt(np.maximum(r2, 1e-36))
        return self.K_r(r)

    @staticmethod
    def K_r(r):
        raise NotImplementedError('kernel not implemented')

    def kernel_to_state_space(self, R=None):
        raise NotImplementedError

    def measurement_model(self):
        raise NotImplementedError

    def state_transition(self, dt):
        raise NotImplementedError

    def stationary_covariance(self):
        raise NotImplementedError

class Matern12(StationaryKernel):
    """
    The Matern 1/2 kernel. Functions drawn from a GP with this kernel are not
    differentiable anywhere. The kernel equation is

    k(r) = σ² exp{-r}

    where:
    r  is the Euclidean distance between the input points, scaled by the lengthscales parameter ℓ.
    σ² is the variance parameter
    """

    @property
    def state_dim(self):
        return 1

    def K_r(self, r):
        return self.variance * np.exp(-r)

    def kernel_to_state_space(self, R=None):
        F = np.array([[-1.0 / self.lengthscale]])
        L = np.array([[1.0]])
        Qc = np.array([[2.0 * self.variance / self.lengthscale]])
        H = np.array([[1.0]])
        Pinf = np.array([[self.variance]])
        return F, L, Qc, H, Pinf

    def stationary_covariance(self):
        Pinf = np.array([[self.variance]])
        return Pinf

    def measurement_model(self):
        H = np.array([[1.0]])
        return H

    def state_transition(self, dt):
        """
        Calculation of the discrete-time state transition matrix A = expm(FΔt) for the exponential prior.
        :param dt: step size(s), Δtₙ = tₙ - tₙ₋₁ [scalar]
        :return: state transition matrix A [1, 1]
        """
        A = np.broadcast_to(np.exp(-dt / self.lengthscale), [1, 1])
        return A

class Matern32(StationaryKernel):
    """
    The Matern 3/2 kernel. Functions drawn from a GP with this kernel are once
    differentiable. The kernel equation is

    k(r) = σ² (1 + √3r) exp{-√3 r}

    where:
    r  is the Euclidean distance between the input points, scaled by the lengthscales parameter ℓ,
    σ² is the variance parameter.
    """

    @property
    def state_dim(self):
        return 2

    def K_r(self, r):
        sqrt3 = np.sqrt(3.0)
        return self.variance * (1.0 + sqrt3 * r) * np.exp(-sqrt3 * r)

    def kernel_to_state_space(self, R=None):
        lam = 3.0 ** 0.5 / self.lengthscale
        F = np.array([[0.0,       1.0],
                      [-lam ** 2, -2 * lam]])
        L = np.array([[0],
                      [1]])
        Qc = np.array([[12.0 * 3.0 ** 0.5 / self.lengthscale ** 3.0 * self.variance]])
        H = np.array([[1.0, 0.0]])
        Pinf = np.array([[self.variance, 0.0],
                         [0.0, 3.0 * self.variance / self.lengthscale ** 2.0]])
        return F, L, Qc, H, Pinf

    def stationary_covariance(self):
        Pinf = np.array([[self.variance, 0.0],
                         [0.0, 3.0 * self.variance / self.lengthscale ** 2.0]])
        return Pinf

    def measurement_model(self):
        H = np.array([[1.0, 0.0]])
        return H

    def state_transition(self, dt):
        """
        Calculation of the discrete-time state transition matrix A = expm(FΔt) for the Matern-3/2 prior.
        :param dt: step size(s), Δtₙ = tₙ - tₙ₋₁ [scalar]
        :return: state transition matrix A [2, 2]
        """
        lam = np.sqrt(3.0) / self.lengthscale
        A = np.exp(-dt * lam) * (dt * np.array([[lam, 1.0], [-lam**2.0, -lam]]) + np.eye(2))
        return A

ext: py

max_stars_repo_name: AntixK/parallel-SSVI

max_issues_repo_name: AntixK/parallel-SSVI

max_forks_repo_name: AntixK/parallel-SSVI

reason: import jax,from jax


text:
import numpy as np
from typing import List
from abc import ABCMeta, abstractmethod

from mot.utils import Registry
from mot.structures import Tracklet, Prediction

__all__ = ['Predictor', 'PREDICTOR_REGISTRY', 'build_predictor']

PREDICTOR_REGISTRY = Registry('predictors')

class Predictor(metaclass=ABCMeta):
    @abstractmethod
    def __init__(self):
        pass

    def __call__(self, tracklets: List[Tracklet], img: np.ndarray):
        return self.predict(tracklets, img)

    @abstractmethod
    def initiate(self, tracklets: List[Tracklet]) -> None:
        """
        Initiate tracklets' states that are used by the predictor.
        """
        pass

    @abstractmethod
    def update(self, tracklets: List[Tracklet]) -> None:
        """
        Update tracklets' states that are used by the predictor.
        """
        pass

    @abstractmethod
    def predict(self, tracklets: List[Tracklet], img: np.ndarray) -> List[Prediction]:
        """
        Predict state in the following time step.

        Args:
            tracklets: A list of tracklet objects. The tracklets that require motion prediction.

        Returns:
            A list of Prediction objects corresponding to the tracklets.
        """
        pass

def build_predictor(cfg):
    if cfg is not None:
        return PREDICTOR_REGISTRY.get(cfg.type)(**(cfg.to_dict()))
    else:
        return None

reason: import numpy


text:
#!/usr/bin/env python

#  Author: Laur Ferris
#  Email address: lnferris@alum.mit.edu
#  Website: https://github.com/lnferris/ocean_data_tools
#  Jun 2020; Last revision: 16-Jun-2020
#  Distributed under the terms of the MIT License

import datetime
import glob
import netCDF4
import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import gsw

# getStations() finds files in local directory that match full_path.
# Writes stations into a pandas dataframe and returns dataframe.
def getStations(full_path):
    filename_list = glob.glob(full_path) # Search for files in full_path.
    print('\nMatching files:\n'+'\n'.join(filename_list)) # Print a list of matching files.

    Hydro_Dataframe = pd.DataFrame(columns=["CTDPRS","CTDSAL","CTDTMP","CTDOXY","LAT","LON","UTC"]) # Make an empty dataframe to hold profile data.

    for filename in filename_list: # For each file in full_path...
        nc = netCDF4.Dataset(filename) # Call the Dataset constructor.

        LAT = nc.variables['latitude'][:]
        LON = nc.variables['longitude'][:]
        CTDTMP = nc.variables['temperature'][:]
        CTDSAL = nc.variables['salinity'][:]
        CTDPRS = nc.variables['pressure'][:]
        CTDOXY = nc.variables['oxygen'][:]
        UTC = nc.variables['woce_date'][:]

        Station_Dataframe = pd.DataFrame({"CTDPRS":[CTDPRS],"CTDSAL":[CTDSAL],"CTDTMP":[CTDTMP],"CTDOXY":[CTDOXY],"LAT":LAT,"LON":LON,"UTC":UTC},index=[0]) # Write station into a temporary dataframe.
        Hydro_Dataframe= pd.concat([Hydro_Dataframe, Station_Dataframe],axis=0,sort=False) # Combine temporary with general dataframe.

        nc.close() # Close the file.

    # Convert WOCE date (YYYMMDD integer format) into a datetime variable.
    Hydro_Dataframe['UTC'] =  pd.to_datetime(Hydro_Dataframe['UTC'], format='%Y%m%d')

    print(Hydro_Dataframe) # Print the dataframe.
    print(Hydro_Dataframe.shape)
    return(Hydro_Dataframe) # Return the list and the dataframe.

# map_dataframe() plots latitude vs. longitude of profiles over basemap.
def map_dataframe(Hydro_Dataframe,BasemapLimits):
    south, north, west, east = BasemapLimits[:]
    map = Basemap(projection='cyl',llcrnrlon=west,llcrnrlat=south,urcrnrlon=east,urcrnrlat=north,resolution='l') # c'l'ihf (Choose basemap resolution.)
    map.drawcoastlines()
    map.fillcontinents(color='black',lake_color='white')
    map.drawparallels(np.arange(south,north,5.0))
    map.drawmeridians(np.arange(west,east,5.0))
    plt.plot(Hydro_Dataframe.LON,Hydro_Dataframe.LAT,marker='.',linestyle='None')
    plt.show()

# Example of how plot oxygen vs. depth for each station in dataframe.
def vertical_profile(Hydro_Dataframe):
    for index, row in Hydro_Dataframe.iterrows():
        tracer = row["CTDOXY"]
        pres = row["CTDPRS"]
        plt.plot(tracer,-pres,marker='.',linestyle='None',markersize=1)
    plt.show()

# Calculate and plot density profiles using Thermodynamic Equation of Seawater 2010 (TEOS-10).
def density_profile(Hydro_Dataframe):
    for index, row in Hydro_Dataframe.iterrows():
        t = row["CTDTMP"] #ITS-90
        SP = row["CTDSAL"] # PSS-78
        p = row["CTDPRS"] # decibar
        lon = [row["LON"]]*t.size
        lat = [row["LAT"]]*t.size
        SA = gsw.SA_from_SP(SP, p, lon, lat) #gsw.SA_from_SP(SP, p, lon, lat)
        CT = gsw.CT_from_t(SA,t,p) # gsw.CT_from_t(SA, t, p)
        D = gsw.density.sigma0(SA, CT)# gsw.density.sigma0(SA, CT)
        plt.plot(D+1000,-p,marker='.',linestyle='None',markersize=1)
    plt.show()

# Example of how to subset dataframe by date of data collection.
def subset_by_date(Hydro_Dataframe):
    date1 = datetime.datetime(2015, 6, 15)
    date2 = datetime.datetime(2016, 6, 25)
    return Hydro_Dataframe[ (Hydro_Dataframe.UTC >= date1) & (Hydro_Dataframe.UTC <= date2 ) ]

# Example of how to subset dataframe by location of station.
def subset_by_region(Hydro_Dataframe):
    lat_s = 48.0 # [-90.0 90.0]
    lat_n = 60.0
    lon_w = -160.0 # [-180.0 180.0]
    lon_e = -145.0
    return Hydro_Dataframe[ (Hydro_Dataframe.LAT >= lat_s) & (Hydro_Dataframe.LAT <= lat_n ) & \
                            (Hydro_Dataframe.LON >= lon_w) & (Hydro_Dataframe.LON <= lon_e)]

#___________Example_Script______________________________________________________

full_path = '/Users/lnferris/Documents/data/goship/P16N/*.nc'
BasemapLimits = [40.0, 65.0, -155.0, -125.0] # Dimensions of lat/lon map.

# Write station data from files in full_path to dataframe.
Hydro_Dataframe = getStations(full_path)

# Map profiles in dataframe.
map_dataframe(Hydro_Dataframe,BasemapLimits)

# Access and plot vertical profile data.
vertical_profile(Hydro_Dataframe)

# Calculate and plot density profiles using Thermodynamic Equation of Seawater 2010 (TEOS-10).
density_profile(Hydro_Dataframe)

# Subset dataframe by date.
Small_Date_Frame = subset_by_date(Hydro_Dataframe)

# Subset dataframe by location of station.
Small_Region_Frame = subset_by_region(Hydro_Dataframe)


text:
import numpy as np
import pandas as pd
from agent.util import load_valid_data, get_sim_dir
from analyze.util import save_dict
from featnames import X_OFFER, NORM, SLR, LOOKUP, START_PRICE
from utils import safe_reindex
from agent.const import DELTA_BYR
from constants import IDX

def get_total_con(data=None):
    norm = data[X_OFFER][NORM].unstack()
    norm = norm[norm[1] < 1]
    norm.loc[:, IDX[SLR]] = 1 - norm.loc[:, IDX[SLR]]
    start_price = safe_reindex(data[LOOKUP][START_PRICE], idx=norm.index)
    for t in norm.columns:
        norm[t] *= start_price
    con = pd.DataFrame(index=norm.index)
    con[1] = norm[1]
    con[2] = start_price - norm[2]
    for t in range(3, 8):
        con[t] = np.abs(norm[t - 2] - norm[t])
    for t in con.columns:
        con[t] /= start_price
    return con.fillna(0).mean()

def main():
    df = pd.DataFrame()

    # human sellers
    data_obs = load_valid_data(byr=True, minimal=True)
    df['Humans'] = get_total_con(data=data_obs)

    # agent buyers
    for delta in DELTA_BYR:
        run_dir = get_sim_dir(byr=True, delta=delta)
        data_rl = load_valid_data(sim_dir=run_dir, minimal=True)
        df['$\\lambda = {}$'.format(delta)] = get_total_con(data=data_rl)

    d = {'stacked_con': df}
    save_dict(d, 'byrcon')

if __name__ == '__main__':
    main()


text:
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftn, ifftn, fftshift, ifftshift
from scipy.ndimage import uniform_filter

from ..compute.reconstruct_phase_util import *
from ..datastructures import IntensityData, StokesData, PhysicalData

class phase_reconstructor:

    '''

    phase_reconstructor contains methods to compute weak object transfer function
    and conduct 2D/2.5D/3D phase reconstruction with a through-focus intensity stack

    Parameters
    ----------
        img_dim         : tuple
                          shape of the computed 3D space with size of (N, M, N_defocus)

        lambda_illu     : float
                          wavelength of the incident light

        ps              : float
                          xy pixel size of the image space

        psz             : float
                          z step size of the image space

        NA_obj          : float
                          numerical aperture of the detection objective

        NA_illu         : float
                          numerical aperture of the illumination condenser

        focus_idx       : int
                          z index of the focus layer in the stack

        n_objective_media         : float
                          refractive index of the immersing media

        phase_deconv    : list
                          a list of string contains the deconvolution dimension on the data
                          '2D'      for 2D phase deconvolution
                          'semi-3D' for 2.5D phase deconvolution
                          '3D'      for 3D phase deconvolution

        ph_deconv_layer : int
                          number of layers included for each layer of semi-3D phase reconstruction

        pad_z           : int
                          number of z-layers to pad (padded with mean of the stack) for 3D phase reconstruction

        use_gpu         : bool
                          option to use gpu or not

        gpu_id          : int
                          number refering to which gpu will be used

    '''

    def __init__(self, img_dim, lambda_illu, ps, psz, NA_obj, NA_illu, focus_idx=None, \
                 n_objective_media=1, phase_deconv=['2D'], ph_deconv_layer=5, pad_z=0, use_gpu=False, gpu_id=0):

        # GPU/CPU

        self.use_gpu = use_gpu
        self.gpu_id = gpu_id

        if self.use_gpu:
            try:
                globals()['cp'] = __import__("cupy")
                cp.cuda.Device(self.gpu_id).use()
                print("cupy found, using GPU")
            except ModuleNotFoundError:
                print("cupy not installed, using CPU instead")
                self.use_gpu = False

        # Basic parameter
        self.N, self.M, self.N_defocus   = img_dim
        self.n_objective_media                     = n_objective_media
        self.lambda_illu                 = lambda_illu/n_objective_media
        self.ps                          = ps
        self.psz                         = np.abs(psz)
        self.pad_z                       = pad_z
        self.phase_deconv                = phase_deconv
        self.focus_idx                   = focus_idx

        if psz >0 and focus_idx is not None:
            self.z_defocus               = -(np.r_[:self.N_defocus]-focus_idx)*self.psz
        elif psz <0 and focus_idx is not None:
            self.z_defocus               = (np.r_[:self.N_defocus]-focus_idx)*self.psz
        elif psz >0 and focus_idx is None:
            self.z_defocus               = -(np.r_[:self.N_defocus]-self.N_defocus//2)*self.psz
        elif psz <0 and focus_idx is None:
            self.z_defocus               = (np.r_[:self.N_defocus]-self.N_defocus//2)*self.psz

        self.NA_obj      = NA_obj/n_objective_media
        self.NA_illu     = NA_illu/n_objective_media

        # setup coordinate systems
        self.xx, self.yy, self.fxx, self.fyy = gen_coordinate((self.N, self.M), ps)

        # detection setup
        self.Pupil_obj     = gen_Pupil(self.fxx, self.fyy, self.NA_obj, self.lambda_illu)
        self.Pupil_support = self.Pupil_obj.copy()

        # illumination setup
        self.Source = gen_Pupil(self.fxx, self.fyy, self.NA_illu, self.lambda_illu)

        # select either 2D or 3D model for deconvolution
        self.phase_deconv_setup(self.phase_deconv, ph_deconv_layer)

    def phase_deconv_setup(self, phase_deconv, ph_deconv_layer):

        '''

        setup transfer functions for phase deconvolution with the corresponding dimensions

        Parameters
        ----------
            phase_deconv    : list
                              a list of string contains the deconvolution dimension on the data
                              '2D'      for 2D phase deconvolution
                              'semi-3D' for 2.5D phase deconvolution
                              '3D'      for 3D phase deconvolution

            ph_deconv_layer : int
                              number of layers included for each layer of semi-3D phase reconstruction

        '''

        if '2D' in phase_deconv:

            Hz_det = gen_Hz_stack(self.fxx, self.fyy, self.Pupil_support, self.lambda_illu, self.z_defocus, type='Prop')
            self.gen_2D_WOTF(Hz_det)

        if 'semi-3D' in phase_deconv:

            self.ph_deconv_layer = ph_deconv_layer
            if self.z_defocus[0] - self.z_defocus[1] >0:
                z_deconv = -(np.r_[:self.ph_deconv_layer]-self.ph_deconv_layer//2)*self.psz
            else:
                z_deconv = (np.r_[:self.ph_deconv_layer]-self.ph_deconv_layer//2)*self.psz

            Hz_det  = gen_Hz_stack(self.fxx, self.fyy, self.Pupil_support, self.lambda_illu, z_deconv, type='Prop')
            G_fun_z = gen_Hz_stack(self.fxx, self.fyy, self.Pupil_support, self.lambda_illu, z_deconv, type='Green')
            self.gen_semi_3D_WOTF(Hz_det, G_fun_z)

        if '3D' in phase_deconv:

            self.N_defocus_3D = self.N_defocus + 2*self.pad_z

            if self.z_defocus[0] - self.z_defocus[1] >0:
                z = -ifftshift((np.r_[0:self.N_defocus_3D]-self.N_defocus_3D//2)*self.psz)
            else:
                z = ifftshift((np.r_[0:self.N_defocus_3D]-self.N_defocus_3D//2)*self.psz)

            Hz_det = gen_Hz_stack(self.fxx, self.fyy, self.Pupil_support, self.lambda_illu, z, type='Prop')
            G_fun_z = gen_Hz_stack(self.fxx, self.fyy, self.Pupil_support, self.lambda_illu, z, type='Green')
            self.gen_3D_WOTF(Hz_det, G_fun_z)

    def gen_2D_WOTF(self, Hz_det):

        '''

        setup transfer functions for 2D phase deconvolution

        Parameters
        ----------
            Hz_det : numpy.ndarray
                     propagation kernel with size of (N, M, N_defocus)

        '''

        self.Hu = np.zeros((self.N, self.M, self.N_defocus),complex)
        self.Hp = np.zeros((self.N, self.M, self.N_defocus),complex)

        for i in range(self.N_defocus):
            self.Hu[:,:,i], self.Hp[:,:,i] = WOTF_2D_compute(self.Source, self.Pupil_obj * Hz_det[:,:,i], \
                                                             use_gpu=self.use_gpu, gpu_id=self.gpu_id)

    def gen_semi_3D_WOTF(self, Hz_det, G_fun_z):

        '''

        setup transfer functions for semi-3D phase deconvolution

        Parameters
        ----------
            Hz_det  : numpy.ndarray
                      propagation kernel with size of (N, M, ph_deconv_layer)

            G_fun_z : numpy.ndarray
                      Weyl's representation of Green's function with size of (N, M, ph_deconv_layer)

        '''

        self.Hu_semi3D = np.zeros((self.N, self.M, self.ph_deconv_layer),complex)
        self.Hp_semi3D = np.zeros((self.N, self.M, self.ph_deconv_layer),complex)

        for i in range(self.ph_deconv_layer):
            self.Hu_semi3D[:,:,i], self.Hp_semi3D[:,:,i] = WOTF_semi_3D_compute(self.Source, self.Pupil_obj, Hz_det[:,:,i], \
                                                                                G_fun_z[:,:,i]*4*np.pi*1j/self.lambda_illu, \
                                                                                use_gpu=self.use_gpu, gpu_id=self.gpu_id)

    def gen_3D_WOTF(self, Hz_det, G_fun_z):

        '''

        setup transfer functions for 3D phase deconvolution

        Parameters
        ----------
            Hz_det  : numpy.ndarray
                      propagation kernel with size of (N, M, N_defocus_3D)

            G_fun_z : numpy.ndarray
                      Weyl's representation of Green's function with size of (N, M, N_defocus_3D)

        '''

        self.H_re = np.zeros((self.N, self.M, self.N_defocus_3D),dtype='complex64')
        self.H_im = np.zeros((self.N, self.M, self.N_defocus_3D),dtype='complex64')

        self.H_re, self.H_im = WOTF_3D_compute(self.Source.astype('float32'), self.Pupil_obj.astype('complex64'), Hz_det.astype('complex64'), \
                                               G_fun_z.astype('complex64'), self.psz, use_gpu=self.use_gpu, gpu_id=self.gpu_id)

    def Phase_solver_para_setter(self, denoiser_2D='Tikhonov',
                                       Tik_reg_abs_2D = 1e-6, Tik_reg_ph_2D = 1e-6, \
                                       TV_reg_abs_2D  = 1e-3, TV_reg_ph_2D  = 1e-3,
                                       rho_2D         = 1e-5, itr_2D        = 20, \
                                       denoiser_3D    ='Tikhonov',
                                       Tik_reg_ph_3D  = 1e-4, TV_reg_ph_3D  = 1e-3,
                                       rho_3D         = 1e-5, itr_3D        = 20, \
                                       verbose        = True, bg_filter     = True):

        '''

        setup parameters for phase deconvolution with the corresponding dimensions

        Parameters
        ----------
            denoiser_2D    : str
                             denoiser for 2D and semi-3D phase reconstruction
                             'Tikhonov' for Tikhonov denoiser
                             'TV'       for TV denoiser

            Tik_reg_abs_2D : float
                             Tikhonov regularization parameter for 2D and semi-3D absorption

            Tik_reg_ph_2D  : float
                             Tikhonov regularization parameter for 2D and semi-3D phase

            TV_reg_abs_2D  : float
                             TV regularization parameter for 2D and semi-3D absorption

            TV_reg_ph_2D   : float
                             TV regularization parameter for 2D and semi-3D absorption

            rho_2D         : float
                             augmented Lagrange multiplier for 2D amd semi-3D ADMM algorithm

            itr_2D         : int
                             number of iterations for 2D and semi-3D ADMM algorithm

            denoiser_3D    : str
                             denoiser for 3D phase reconstruction
                             'Tikhonov' for Tikhonov denoiser
                             'TV'       for TV denoiser

            Tik_reg_ph_3D  : float
                             Tikhonov regularization parameter for 3D phase

            TV_reg_ph_3D   : float
                             TV regularization parameter for 3D phase

            rho_3D         : float
                             augmented Lagrange multiplier for 3D ADMM algorithm

            itr_3D         : int
                             number of iterations for 3D ADMM algorithm

            verbose        : bool
                             option to display detailed progress of computations or not

            bg_filter      : bool
                             option for slow-varying 2D background normalization with uniform filter
        '''

        valid_denoiser_list = ['Tikhonov','TV']

        if denoiser_2D in valid_denoiser_list:
            self.denoiser_2D = denoiser_2D
        else:
            self.denoiser_2D = 'Tikhonov'
            print("denoiser_2D must be 'Tikhonov' or 'TV, set it to 'Tikhonov' in default")

        if denoiser_3D in valid_denoiser_list:
            self.denoiser_3D = denoiser_3D
        else:
            self.denoiser_3D = 'Tikhonov'
            print("denoiser_3D must be 'Tikohonov' or 'TV, set it to 'Tikhonov' in default")

        self.Tik_reg_abs_2D = Tik_reg_abs_2D
        self.Tik_reg_ph_2D  = Tik_reg_ph_2D

        self.Tik_reg_ph_3D  = Tik_reg_ph_3D

        self.TV_reg_abs_2D  = TV_reg_abs_2D
        self.TV_reg_ph_2D   = TV_reg_ph_2D
        self.rho_2D         = rho_2D
        self.itr_2D         = itr_2D

        self.TV_reg_ph_3D   = TV_reg_ph_3D
        self.rho_3D         = rho_3D
        self.itr_3D         = itr_3D

        self.verbose        = verbose
        self.bg_filter      = bg_filter

    def Phase_recon_2D(self, stoke_stack: StokesData):

        '''

        conduct 2D phase reconstruction

        Parameters
        ----------
            stoke_stack : StokesData
                          StokesData of the sample with each channel having (N, M, N_defocus)-sized numpy.ndarray

        Returns
        -------
            mu_sample  : numpy.ndarray
                         2D absorption reconstruction with the size of (N, M)

            phi_sample : numpy.ndarray
                         2D phase reconstruction with the size of (N, M)

        '''

        S0_stack = stoke_stack.s0

        if self.use_gpu:

            S0_stack = inten_normalization(cp.array(S0_stack), type='2D', bg_filter=self.bg_filter, use_gpu=True, gpu_id=self.gpu_id)

            Hu = cp.array(self.Hu, copy=True)
            Hp = cp.array(self.Hp, copy=True)

            S0_stack_f = cp.fft.fft2(S0_stack, axes=(0,1))

            AHA   = [cp.sum(cp.abs(Hu)**2, axis=2) + self.Tik_reg_abs_2D, cp.sum(cp.conj(Hu)*Hp, axis=2),\
                     cp.sum(cp.conj(Hp)*Hu, axis=2),                      cp.sum(cp.abs(Hp)**2, axis=2) + self.Tik_reg_ph_2D]

            b_vec = [cp.sum(cp.conj(Hu)*S0_stack_f, axis=2), \
                     cp.sum(cp.conj(Hp)*S0_stack_f, axis=2)]

        else:
            S0_stack   = inten_normalization(S0_stack, type='2D', bg_filter=self.bg_filter)
            S0_stack_f = fft2(S0_stack,axes=(0,1))

            AHA   = [np.sum(np.abs(self.Hu)**2, axis=2) + self.Tik_reg_abs_2D, np.sum(np.conj(self.Hu)*self.Hp, axis=2),\
                     np.sum(np.conj(self.Hp)*self.Hu, axis=2),                 np.sum(np.abs(self.Hp)**2, axis=2) + self.Tik_reg_ph_2D]

            b_vec = [np.sum(np.conj(self.Hu)*S0_stack_f, axis=2), \
                     np.sum(np.conj(self.Hp)*S0_stack_f, axis=2)]

        if self.denoiser_2D == 'Tikhonov':

            # Deconvolution with Tikhonov regularization

            mu_sample, phi_sample = WOTF_Tikhonov_deconv_2D(AHA, b_vec, use_gpu=self.use_gpu, gpu_id=self.gpu_id)

        elif self.denoiser_2D == 'TV':

            # ADMM deconvolution with anisotropic TV regularization

            mu_sample, phi_sample = WOTF_ADMM_TV_deconv_2D(AHA, b_vec, self.rho_2D, self.TV_reg_abs_2D, self.TV_reg_ph_2D, \
                                                           self.itr_2D, self.verbose, use_gpu=self.use_gpu, gpu_id=self.gpu_id)
        else:
            raise TypeError("denoiser type must be 'Tikhonov' or 'TV'")

        phi_sample -= phi_sample.mean()

        return mu_sample, phi_sample

    def Phase_recon_semi_3D(self, stoke_stack: StokesData):

        '''

        conduct semi-3D phase reconstruction

        Parameters
        ----------
            stoke_stack : StokesData
                          StokesData of the sample with each channel having (N, M, N_defocus)-sized numpy.ndarray

        Returns
        -------
            mu_sample  : numpy.ndarray
                         semi-3D absorption reconstruction with the size of (N, M, N_defocus)

            phi_sample : numpy.ndarray
                         semi-3D phase reconstruction with the size of (N, M, N_defocus)

        '''

        S0_stack   = stoke_stack.s0

        mu_sample  = np.zeros((self.N, self.M, self.N_defocus))
        phi_sample = np.zeros((self.N, self.M, self.N_defocus))

        for i in range(self.N_defocus):

            # determines which indices to slice for each semi-3D reconstruction
            if i <= self.ph_deconv_layer//2:
                tf_start_idx = self.ph_deconv_layer//2 - i
            else:
                tf_start_idx = 0

            obj_start_idx = np.maximum(0,i-self.ph_deconv_layer//2)

            if self.N_defocus -i -1 < self.ph_deconv_layer//2:
                tf_end_idx = self.ph_deconv_layer//2 + (self.N_defocus - i)
            else:
                tf_end_idx = self.ph_deconv_layer

            obj_end_idx = np.minimum(self.N_defocus,i+self.ph_deconv_layer-self.ph_deconv_layer//2)

            if self.verbose:
                print('TF_index = (%d,%d), obj_z_index=(%d,%d), consistency: %s'\
                      %(tf_start_idx,tf_end_idx, obj_start_idx, obj_end_idx, (obj_end_idx-obj_start_idx)==(tf_end_idx-tf_start_idx)))

            if self.use_gpu:
                S0_stack_sub = inten_normalization(cp.array(S0_stack[:,:,obj_start_idx:obj_end_idx]), type='2D', bg_filter=self.bg_filter, use_gpu=True, gpu_id=self.gpu_id)
                S0_stack_f   = cp.fft.fft2(S0_stack_sub, axes=(0,1))

                Hu = cp.array(self.Hu_semi3D[:,:,tf_start_idx:tf_end_idx], copy=True)
                Hp = cp.array(self.Hp_semi3D[:,:,tf_start_idx:tf_end_idx], copy=True)

                AHA   = [cp.sum(cp.abs(Hu)**2, axis=2) + self.Tik_reg_abs_2D, cp.sum(cp.conj(Hu)*Hp, axis=2),\
                         cp.sum(cp.conj(Hp)*Hu, axis=2),                      cp.sum(cp.abs(Hp)**2, axis=2) + self.Tik_reg_ph_2D]

                b_vec = [cp.sum(cp.conj(Hu)*S0_stack_f, axis=2), \
                         cp.sum(cp.conj(Hp)*S0_stack_f, axis=2)]

            else:
                S0_stack_sub = inten_normalization(S0_stack[:,:,obj_start_idx:obj_end_idx], type='2D', bg_filter=self.bg_filter)
                S0_stack_f   = fft2(S0_stack_sub,axes=(0,1))

                Hu = self.Hu_semi3D[:,:,tf_start_idx:tf_end_idx]
                Hp = self.Hp_semi3D[:,:,tf_start_idx:tf_end_idx]

                AHA   = [np.sum(np.abs(Hu)**2, axis=2) + self.Tik_reg_abs_2D, np.sum(np.conj(Hu)*Hp, axis=2),\
                         np.sum(np.conj(Hp)*Hu, axis=2),                      np.sum(np.abs(Hp)**2, axis=2) + self.Tik_reg_ph_2D]

                b_vec = [np.sum(np.conj(Hu)*S0_stack_f, axis=2), \
                         np.sum(np.conj(Hp)*S0_stack_f, axis=2)]

            if self.denoiser_2D == 'Tikhonov':

                # Deconvolution with Tikhonov regularization

                mu_sample_temp, phi_sample_temp = WOTF_Tikhonov_deconv_2D(AHA, b_vec, use_gpu=self.use_gpu, gpu_id=self.gpu_id)

            elif self.denoiser_2D == 'TV':

                # ADMM deconvolution with anisotropic TV regularization

                mu_sample_temp, phi_sample_temp = WOTF_ADMM_TV_deconv_2D(AHA, b_vec, self.rho_2D, self.TV_reg_abs_2D, self.TV_reg_ph_2D, \
                                                                         self.itr_2D, self.verbose, use_gpu=self.use_gpu, gpu_id=self.gpu_id)
            else:
                raise TypeError("denoiser type must be 'Tikhonov' or 'TV'")

            mu_sample[:,:,i]  = mu_sample_temp.copy()
            phi_sample[:,:,i] = phi_sample_temp - phi_sample_temp.mean()

        return mu_sample, phi_sample

    def Phase_recon_3D(self, stoke_stack: StokesData, absorption_ratio=0.0):

        '''

        conduct semi-3D phase reconstruction

        Parameters
        ----------
            stoke_stack      : StokesData
                               StokesData of the sample with each channel having (N, M, N_defocus)-sized numpy.ndarray

            absorption_ratio : float
                               absorption = absorption_ratio * phase (absorption_ratio = 0 refers to pure phase approximation)

        Returns
        -------
            f_real           : numpy.ndarray
                               3D phase reconstruction with the size of (N, M, N_defocus)

        '''

        S0_stack = stoke_stack.s0

        if self.pad_z == 0:
            S0_stack = inten_normalization(S0_stack, type='3D')
        else:
            S0_pad = np.pad(S0_stack,((0,0),(0,0),(self.pad_z,self.pad_z)), mode='constant',constant_values=S0_stack.mean())
            if self.pad_z < self.N_defocus:
                S0_pad[:,:,:self.pad_z] = (S0_stack[:,:,:self.pad_z])[:,:,::-1]
                S0_pad[:,:,-self.pad_z:] = (S0_stack[:,:,-self.pad_z:])[:,:,::-1]
            else:
                print('pad_z is larger than number of z-slices, use mean value padding (not effective) instead of reflection padding')

            S0_stack = inten_normalization(S0_pad, type='3D')

        H_eff = self.H_re + absorption_ratio*self.H_im

        if self.denoiser_3D == 'Tikhonov':

            f_real = WOTF_Tikhonov_deconv_3D(S0_stack, H_eff, self.Tik_reg_ph_3D, use_gpu=self.use_gpu, gpu_id=self.gpu_id)

        elif self.denoiser_3D == 'TV':

            f_real = WOTF_ADMM_TV_deconv_3D(S0_stack, H_eff, self.rho_3D, self.Tik_reg_ph_3D, self.TV_reg_ph_3D, self.itr_3D,\
                                            self.verbose, use_gpu=self.use_gpu, gpu_id=self.gpu_id)
        else:
            raise TypeError("denoiser type must be 'Tikhonov' or 'TV'")

        if self.pad_z != 0:
            f_real = f_real[:,:,self.pad_z:-(self.pad_z)]

        return -f_real*self.psz/(4*np.pi/self.lambda_illu)

max_stars_repo_name: czbiohub/reconstruct-order

max_issues_repo_name: czbiohub/ReconstructOrder

max_forks_repo_name: mehta-lab/reconstruct-order

reason: import numpy,from numpy,from scipy


text:
import numpy as np
import xarray as xr

import python2latex as p2l

import os
path = os.path.dirname(__file__)

def plot_mprime_best(ms, risks, d, delta):
    ms, risks = np.array(ms), np.array(risks)

    plot = p2l.Plot(plot_name=f'mprime_best{d=}',
                    plot_path=path+'/figures',
                    as_float_env=False,
                    width='7.45cm',
                    height='5.45cm',
                    marks='1.5pt',
                    lines=False,
                    palette='holi5')

    best_mprimes = xr.open_dataset(path+'/data/optimal_bound.nc')['mprime'].drop_sel(d=50)

    for m in ms:
        plot_mprimes = best_mprimes.sel(m=m, d=d, delta=delta)
        plot.add_plot(risks, plot_mprimes/m, legend=f'\\scriptsize $m={m}$')

    plot.x_min = -0.05
    plot.x_max = 0.55
    plot.x_ticks = 0, .1, .2, .3, .4, .5
    plot.x_label = r"$k/m$"
    plot.y_label = r"$\displaystyle\frac{m'_{\textrm{\scriptsize best}}}{m}$"
    plot.legend_position = 'north west'
    plot.axis.kwoptions['ylabel style'] = r'{rotate=-90, xshift=.3cm}'
    plot.axis.kwoptions['legend cell align'] = '{left}'

    return plot

if __name__ == "__main__":
    d = 10
    # d = 35
    ms = [100, 200, 300, 500, 1000]
    risks =  np.linspace(0, .5, 11)
    delta = 0.05

    filename = f'mprime_best_{d=}'
    doc = p2l.Document(filename, filepath=path, doc_type='standalone')

    doc += plot_mprime_best(ms, risks, d, delta)

    print('Building...')
    doc.build(delete_files='all', show_pdf=False)


text:
import codecs
import os

import numpy as np

NUM_TOP_WORDS = 30
VOCABS_PATH = os.path.join('Datasets', 'vocabularies.txt')
TOP_WORDS_PATH_FMT = os.path.join('Datasets', 'top-' + str(NUM_TOP_WORDS) + '_k{}.txt')

if __name__ == '__main__':
    beta = []
    with open('beta.txt', 'r') as f:
        for l in f:
            l = l.strip()
            if not len(l):
                continue

            new_line = l.replace('  ', ' ')
            while len(new_line) != len(l):
                l = new_line
                new_line = l.replace('  ', ' ')

            beta.append([ float(b) for b in l.split()[2:] ])
    beta = np.array(beta)
    idx = np.argsort(beta)

    vocabularies = [ l for l in codecs.open(VOCABS_PATH, 'r', encoding='utf-8') ]
    for k in range(idx.shape[0]):
        NUM_TOP_WORDS_IDX = NUM_TOP_WORDS + 1
        print(idx[k][-1:-NUM_TOP_WORDS_IDX:-1])
        with codecs.open(TOP_WORDS_PATH_FMT.format(k), 'w', encoding='utf-8') as f:
            top_words = u''.join([ vocabularies[v] for v in idx[k][-1:-NUM_TOP_WORDS_IDX:-1] ])
            f.write(top_words)

max_stars_repo_name: VahidHeidari/TopicModeling

max_issues_repo_name: VahidHeidari/TopicModeling

max_forks_repo_name: VahidHeidari/TopicModeling


text:
"""
This file defines some simple experiments to illustrate how Pytorch-Routing functions.
"""
import numpy as np
import tqdm
import torch
from PytorchRouting.DecisionLayers import REINFORCE, QLearning, SARSA, ActorCritic, GumbelSoftmax, PerTaskAssignment, \
    WPL, AAC, AdvantageLearning, RELAX, EGreedyREINFORCE, EGreedyAAC
from PytorchRouting.Examples.Models import PerTask_all_fc, RoutedAllFC, PerTask_1_fc, PerDecisionSingleAgent, \
    Dispatched
from PytorchRouting.Examples.Datasets import CIFAR100MTL

def compute_batch(model, batch):
    samples, labels, tasks = batch
    out, meta = model(samples, tasks=tasks)
    correct_predictions = (out.max(dim=1)[1].squeeze() == labels.squeeze()).cpu().numpy()
    accuracy = correct_predictions.sum()
    oh_labels = one_hot(labels, out.size()[-1])
    module_loss, decision_loss = model.loss(out, meta, oh_labels)
    return module_loss, decision_loss, accuracy

def one_hot(indices, width):
    indices = indices.squeeze().unsqueeze(1)
    oh = torch.zeros(indices.size()[0], width).to(indices.device)
    oh.scatter_(1, indices, 1)
    return oh

def run_experiment(model, dataset, learning_rates, routing_module_learning_rate_ratio):
    print('Loaded dataset and constructed model. Starting Training ...')
    for epoch in range(50):
        optimizers = []
        parameters = []
        if epoch in learning_rates:
            try:
                optimizers.append(torch.optim.SGD(model.routing_parameters(),
                                                  lr=routing_module_learning_rate_ratio*learning_rates[epoch]))
                optimizers.append(torch.optim.SGD(model.module_parameters(),
                                                  lr=learning_rates[epoch]))
                parameters = model.module_parameters() + model.module_parameters()
            except AttributeError:
                optimizers.append(torch.optim.SGD(model.parameters(), lr=learning_rates[epoch]))
                parameters = model.parameters()
        train_log, test_log = np.zeros((3,)), np.zeros((3,))
        train_samples_seen, test_samples_seen = 0, 0
        dataset.enter_train_mode()
        model.train()
        # while True:
        pbar = tqdm.tqdm(unit=' samples')
        while True:
            try:
                batch = dataset.get_batch()
            except StopIteration:
                break
            train_samples_seen += len(batch[0])
            pbar.update(len(batch[0]))
            module_loss, decision_loss, accuracy = compute_batch(model, batch)
            (module_loss + decision_loss).backward()
            torch.nn.utils.clip_grad_norm_(parameters, 40., norm_type=2)
            for opt in optimizers:
                opt.step()
            model.zero_grad()
            train_log += np.array([module_loss.tolist(), decision_loss.tolist(), accuracy])
        pbar.close()
        dataset.enter_test_mode()
        model.eval()
        model.start_logging_selections()
        while True:
            try:
                batch = dataset.get_batch()
            except StopIteration:
                break
            test_samples_seen += len(batch[0])
            module_loss, decision_loss, accuracy = compute_batch(model, batch)
            test_log += np.array([module_loss.tolist(), decision_loss.tolist(), accuracy])
        print('Epoch {} finished after {} train and {} test samples..\n'
              '    Training averages: Model loss: {}, Routing loss: {}, Accuracy: {}\n'
              '    Testing averages:  Model loss: {}, Routing loss: {}, Accuracy: {}'.format(
            epoch + 1, train_samples_seen, test_samples_seen,
            *(train_log/train_samples_seen).round(3), *(test_log/test_samples_seen).round(3)))
        model.stop_logging_selections_and_report()

if __name__ == '__main__':
    # MNIST
    # dataset = MNIST_MTL(64, data_files=['./Datasets/mnist.pkl.gz'])
    # model = PerTask_all_fc(1, 288, 2, dataset.num_tasks, dataset.num_tasks)
    # model = WPL_routed_all_fc(1, 288, 2, dataset.num_tasks, dataset.num_tasks)
    cuda = False
    # cuda = True

    # CIFAR
    dataset = CIFAR100MTL(10, data_files=['./Datasets/cifar-100-py/train', './Datasets/cifar-100-py/test'], cuda=cuda)
    model = RoutedAllFC(WPL, 3, 128, 5, dataset.num_tasks, dataset.num_tasks)
    # model = RoutedAllFC(RELAX, 3, 128, 5, dataset.num_tasks, dataset.num_tasks)
    # model = RoutedAllFC(EGreedyREINFORCE, 3, 128, 5, dataset.num_tasks, dataset.num_tasks)
    # model = RoutedAllFC(AdvantageLearning, 3, 128, 5, dataset.num_tasks, dataset.num_tasks)
    # model = PerDecisionSingleAgent(AdvantageLearning, 3, 128, 5, dataset.num_tasks, dataset.num_tasks)
    # model = Dispatched(AdvantageLearning, 3, 128, 5, dataset.num_tasks, dataset.num_tasks)

    learning_rates = {0: 3e-3, 5: 1e-3, 10: 3e-4}
    routing_module_learning_rate_ratio = 0.3
    if cuda:
        model.cuda()
    run_experiment(model, dataset, learning_rates, routing_module_learning_rate_ratio)

'''
WPL_routed_all_fc(3, 512, 5, dataset.num_tasks, dataset.num_tasks)
    Training averages: Model loss: 0.427, Routing loss: 8.864, Accuracy: 0.711
    Testing averages:  Model loss: 0.459, Routing loss: 9.446, Accuracy: 0.674
'''

max_stars_repo_name: oleksost/RoutingNetworks

max_issues_repo_name: oleksost/RoutingNetworks

max_forks_repo_name: oleksost/RoutingNetworks


text:
from pandas import Series
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as pl
import numpy as np

import os, sys
inserted_path = '/../github_helpers/'
sys.path.insert(0, os.path.realpath(os.path.dirname(__file__) + inserted_path))
import github_helpers

def Main():
  client = github_helpers.authenticate()
  keywords = raw_input("Please, enter keywords to search repositories: ")
  if keywords is '':
    keywords = 'javascript'
    print 'No keywords provided. It will use the keyword: ' + keywords
  search = client.search_repositories(keywords)
  first_page = search.get_page(0)

  languages = Series(r.language for r in first_page)
  languages = languages.dropna()
  languages.sort()

  percentages = (100.0 * languages.value_counts() / len(languages)).map('{:,.2f} %'.format)

  print 'Languages percentage:'
  print percentages

  # Create plot
  x = [int(r.stargazers_count) for r in first_page]
  y = [int(r.forks) for r in first_page]

  # Add one to every value for logarithmic scale
  x = [val + 1 for val in x]
  y = [val + 1 for val in y]

  area = [100 for r in first_page]
  names = [r.name for r in first_page]
  colors = np.random.rand(len(first_page))
  pl.scatter(x, y, s=area, c=colors, alpha=0.5)
  for i in range(0, len(x)):
    pl.annotate(names[i], (x[i], y[i]), fontsize=2)
  pl.title("All values are with addition of 1 (for the logarithmic scale)")
  pl.xlabel("Stars")
  pl.xscale("log")
  pl.yscale("log")
  pl.ylabel("Forks")
  pl.tight_layout()
  filepath = 'reports/APIs/github'
  if not os.path.isdir(filepath): os.makedirs(filepath)
  filepath += '/search_repositories.png'
  pl.savefig(filepath, figsize=(1020, 1020), dpi=300)
  pl.close()
  print('A chart with high resolution and small font size (to minimize overlaps) was created at ' +
    filepath)

if __name__ == "__main__":
  Main()

max_stars_repo_name: igncp/leaves

max_issues_repo_name: igncp/leaves

max_forks_repo_name: igncp/leaves


text:
#!/usr/bin/python3

import matplotlib.pyplot as plt
import numpy as np

folder = 'spitler'
datafiles = ['burst1','burst4','burst8','burst9','burst10','burst11']

for datafile in datafiles:
	data = np.genfromtxt(folder + '/' + datafile + '.csv', delimiter=' ')

	x = data[:,0] * 1e9 # hz
	restfreq = 1.4e9 #hz
	velocity = 3e8 * (x - restfreq) / restfreq # m/s
	velocity = velocity/1000

	y = data[:,1]

	plt.plot(velocity, y, 'k-')
	plt.xlabel('Velocity (km/s)')
	plt.ylabel('S/N')
	plt.title('Burst #' + datafile.split('burst')[1])
	plt.savefig(folder + '/' + datafile + '.pdf')
	plt.clf()
	# plt.show()

	# exit()


text:
"""
plot_obs_coverage.py

Copyright 2016 University of Melbourne.
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License.
You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and limitations under the License.
"""

import os
import numpy as np
import datetime
import matplotlib.pyplot as plt

import context
import fourdvar.user_driver as user
import fourdvar.util.date_handle as dt

obs = user.get_observed()

nrow = obs.grid_attr['NROWS']
ncol = obs.grid_attr['NCOLS']
nday = len( dt.get_datelist() )

tstep = obs.grid_attr['TSTEP']
tsec = (tstep//10000)*60*60 + ((tstep//100)%100)*60 + tstep%100
daysec = 24*60*60
nstep = daysec/tsec

coverage = np.zeros((nrow,ncol,),dtype=int)
obs_step = []

skipped_obs = 0
for weight in obs.weight_grid:
    #limit to obs the reach the surface
    if 0 not in [ k[2] for k in weight.keys() ]:
        skipped_obs += 1
        continue
    #assign obs to surface cell with largest weight
    _,obs_coord = max([ (v,k,) for k,v in weight.items() if k[2]==0 ])
    d,t,_,r,c,_ = obs_coord
    coverage[r,c] += 1
    date = datetime.datetime.strptime( str( d ), '%Y%m%d' ).date()
    obs_step.append( (date-dt.start_date).days*nstep + t )
if skipped_obs > 0:
    print '{:} non-surface obs were omitted.'

omax = coverage.max()

fig = plt.figure( figsize=(10,5) )
plt.subplot(121)
plt.title( '# of observations for each surface gridcell' )
cov = plt.pcolormesh( coverage, vmin=0, cmap='Blues',
                      edgecolor='k', linestyle=':', linewidth=.1 )
cbar = plt.colorbar()
cbar.set_label('# of observations')
cbar.set_ticks( range(0,omax,max(1,omax//5)) )

plt.subplot(122)
plt.title( '# of observations for each timestep' )
plt.hist( obs_step, bins=[ i+.5 for i in range(nday*nstep+1) ] )
plt.xlim(0,nday*nstep+1)
plt.xticks( [ nstep/2 + nstep*i for i in range(nday) ],
            [ dt.replace_date('<YYYY-MM-DD>',d) for d in dt.get_datelist() ] )
for i in range(1,nday):
    plt.axvline( i*nstep, ls=':', c='k' )
plt.show()


text:
import os
import glob
import pickle
import random
import json

import torch
from tqdm import tqdm
import torch.utils.data as data
import numpy as np
import pandas as pd

class SoundDataset(data.Dataset):

    def __init__(
        self, audio_files, labels=None,
        transform=None, is_noisy=None, clean_transform=None):

        self.transform = transform
        self.clean_transform = clean_transform
        self.audio_files = audio_files
        self.labels = labels
        self.is_noisy = is_noisy or np.zeros(len(self.audio_files))

    def __getitem__(self, index):

        sample = dict(
            filename=self.audio_files[index],
            is_noisy=self.is_noisy[index]
        )

        if self.labels is not None:
            sample["raw_labels"] = self.labels[index]

        if self.transform is not None:
            sample = self.transform(dataset=self, **sample)

        return sample

    def random_clean_sample(self):

        index = random.randint(0, len(self) - 1)

        sample = dict(
            filename=self.audio_files[index],
            is_noisy=self.is_noisy[index]
        )

        if self.labels is not None:
            sample["raw_labels"] = self.labels[index]

        if self.clean_transform is not None:
            sample = self.clean_transform(dataset=self, **sample)

        return sample

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


text:
import os
import cv2
import numpy as np
from math import (cos, sin, radians, floor, ceil, exp, log)

def translation(src, x, y):
    out = np.zeros(src.shape, src.dtype)
    X, Y = src.shape[:2]
    out[-X + x: X + x, -Y + y: Y + y] = src[-X - x: X - x, -Y - y: Y - y]
    return out

def resize(src, xfactor, yfactor):
    shape = src.shape
    rows, cols = src.shape[:2]
    if len(shape) == 2:
        shape = (floor(rows * yfactor), floor(cols * xfactor))
    else:
        shape = (floor(rows * yfactor), floor(cols * xfactor), shape[2])

    out = np.empty(shape, src.dtype)

    for i in range(shape[0]):
        for j in range(shape[1]):
            out[i, j] = src[floor(i // yfactor), floor(j // xfactor)]

    return out

def rotate(src, angle, pivot=None):
    """
    :params pivot: Center of rotation.
    """
    shape = src.shape
    out = np.zeros(shape, src.dtype)

    if pivot is None:
        pivot = (shape[0] // 2, shape[1] // 2)

    for i in range(shape[0]):  # rows
        for j in range(shape[1]):  # cols
            x, y = _rotate_pos((shape[1] - j - 1 , i), pivot, angle)
            x = shape[1] - x - 1

            x = round(x)
            y = round(y)

            if 0 <= x < shape[1] and 0 <= y < shape[1]:
                out[y,x] = src[i, j]

    return out

def _rotate_pos(point, center, angle):
    """:params angle: in degrees"""
    angle = radians(angle)
    x, y = point
    cx, cy = center
    _cos = cos(angle)
    _sin = sin(angle)

    x -= cx
    y -= cy

    nx = x * _cos - y * _sin
    ny = x * _sin + y * _cos

    return nx + cx, ny + cy

def exp_filter(x):
    return pow(x, 2.5) % 255

def exp_filter_B(x):
    return pow(x / 255, 2.5) * 255

@np.vectorize
def log_filter(x):
    if x == 0:
        return 0
    x /= 255
    return (log(x) * 255) % 255

@np.vectorize
def log_filter_B(x):
    if x == 0:
        return 0
    x /= 255
    return 0 if log(x) < log(0.5) else 255

@np.vectorize
def inverse_filter(x):
    return 255 - x

if __name__ == "__main__":
    img = cv2.imread('../imgs/cameraman.jpg', 0)

    os.makedirs('res', exist_ok=True)

    positive = translation(img, 50, 40)
    negative = translation(img, -50, -40)
    cv2.imwrite('res/translation_positive.jpg', positive)
    cv2.imwrite('res/translation_negative.jpg', negative)

    cv2.imwrite('res/rotate.jpg', rotate(img, 90))

    cv2.imwrite('res/upscale.jpg', resize(img, 3, 3))
    cv2.imwrite('res/downscale.jpg', resize(img, 0.4, 0.4))

    cv2.imwrite('res/exp.jpg', exp_filter(img))
    cv2.imwrite('res/expB.jpg', exp_filter_B(img))
    cv2.imwrite('res/log.jpg', log_filter(img))
    cv2.imwrite('res/logB.jpg', log_filter_B(img))
    cv2.imwrite('res/inverse.jpg', inverse_filter(img))

max_stars_repo_name: Fredy/UCSP-TopicsGraph

max_issues_repo_name: Fredy/UCSP-TopicsGraph

max_forks_repo_name: Fredy/UCSP-TopicsGraph


text:
from pioneer.common                                import linalg
from pioneer.das.api.datasources                   import AbstractDatasource
from pioneer.das.api.egomotion.egomotion_provider  import EgomotionProvider

import numpy as np
from scipy.interpolate import splprep, splev

class IMUEgomotionProvider(EgomotionProvider):
    """Ego Motion provider for datasets from the RAV-4 platform"""

    def __init__(self, referential_name:str, navposvel:AbstractDatasource, ekfeuler:AbstractDatasource, subsampling:int=100):
        super(IMUEgomotionProvider, self).__init__(referential_name, subsampling)
        self.navposvel = navposvel
        self.ekfeuler = ekfeuler
        self._tf_Global_from_EgoZero = self.get_Global_from_Ego_at(self.navposvel.timestamps[0])

    def get_Global_from_Ego_at(self, ts:int, dtype = np.float64) -> np.ndarray:
        """override"""
        nav = self.navposvel.get_at_timestamp(ts).raw  # latitude, longitude, altitude
        euler = self.ekfeuler.get_at_timestamp(ts).raw # roll, pitch, yaw
        return linalg.tf_from_poseENU(linalg.imu_to_utm(nav, euler), dtype = dtype)

    def get_timestamps_range(self) -> np.ndarray:
        """override"""
        ts = self.navposvel.timestamps
        return np.array([ts[0], ts[-1]])

    def trace_sampled_trajectory(self, sampling_distance:float=1.0, N:int=None) -> np.ndarray:
        xyz = np.empty((0,3))
        for ts in self.navposvel.timestamps:
            xyz = np.vstack([xyz, self.get_Global_from_Ego_at(ts)[:3,3]])

        distances = ((xyz[1:,0] - xyz[:-1,0])**2 + (xyz[1:,1] - xyz[:-1,1])**2)**0.5
        distances = np.insert(distances, 0, 0)
        total_distance = np.sum(distances)

        N = int(total_distance / sampling_distance) + 1 if N is None else N

        # To prevent a crash from scipy's splprep, we have to remove any identical consecutive vertices
        keep = np.where(~((xyz[:-1,0] == xyz[1:,0]) & (xyz[:-1,1] == xyz[1:,1])))[0] + 1
        keep = np.insert(keep, 0, 0)

        tck, u = splprep([xyz[keep,0], xyz[keep,1]])
        u = np.linspace(0,1,N)
        return np.array(splev(u, tck)).T

reason: import numpy,from scipy


text:
"""
现存的问题：可视化的结果是，补丁的宽高已知且固定。但实际上补丁的宽高收到缩放影响，因此在跟踪时才能确定。
版本改进20201222：因为我们训练了一个可以缩放的补丁，所以希望GT相对原始图像是尺寸固定的。这样可以避免补丁越来越大/小。
改进版本20201229：轨迹随机化
改进版本20210103：为其他数据集生成 FGT。
"""
import os
import random
import copy
import imageio
import numpy as np
from PIL import Image, ImageDraw

from videoanalyst.evaluation.got_benchmark.datasets import GOT10k, OTB, LaSOT
from videoanalyst.pipeline.utils import (xywh2cxywh, xywh2xyxy)

def get_crop(target_sz,
             z_size,
             context_amount=0.5):
    wc = target_sz[0] + context_amount * sum(target_sz)
    hc = target_sz[1] + context_amount * sum(target_sz)
    s_crop = np.sqrt(wc * hc)
    scale = z_size / s_crop
    return scale

def visualize(anno, patch_anno, resolution, video_index):
    img_list = []
    img = Image.new('RGB', resolution, (255, 255, 255))
    draw = ImageDraw.Draw(img)
    for anno_, patch_anno_ in zip(anno, patch_anno):
        anno_cxy = xywh2cxywh(anno_)[:2]
        patch_anno_cxy = xywh2cxywh(patch_anno_)[:2]
        print(patch_anno_cxy)
        w=2
        draw.ellipse([tuple(anno_cxy - w), tuple(anno_cxy + w)], fill='green')
        draw.ellipse([tuple(patch_anno_cxy - w), tuple(patch_anno_cxy + w)], fill='red')
        img_list.append(copy.deepcopy(img))
    file_path = "/tmp/anno/{}.gif".format(video_index)
    print(file_path)
    imageio.mimsave(file_path, img_list, fps=20)
    return

def check_boundary(box_xywh, resolution):
    img_w, img_h = resolution
    w, h = box_xywh[2:]
    x1, y1, x2, y2 = xywh2xyxy(box_xywh)

    if x1 < 0 or y1 < 0 or x1 >= img_w or y1 >= img_h:
        return [img_w-w, img_h-h, w, h]
    else:
        return box_xywh

def get_init_patch_box_in_original_img(gt_xywh_in_original_img):
    """
    :param gt_xywh_in_original_img: xywh
    """
    """START：读取原始边框"""
    gt_x1_ori, gt_y1_ori, gt_w_ori, gt_h_ori = gt_xywh_in_original_img
    """END：读取原始边框"""

    """START：计算补丁在原始图像的边框"""
    patch_x1_ori = gt_x1_ori - gt_w_ori
    patch_y1_ori = gt_y1_ori - gt_h_ori
    patch_gt_xywh_ori = [patch_x1_ori, patch_y1_ori, gt_w_ori, gt_h_ori]
    """END：计算补丁在原始图像的边框"""

    return patch_gt_xywh_ori

def generate_random_trajectory(start_y, length, bound_y):
    y = start_y
    result = []
    for _ in range(length):
        if y < 0:
            y += random.randint(0, 32)
        elif y > bound_y:
            y += random.randint(-32, 0)
        else:
            y += random.randint(-32, 32)
        result.append(y)
    return np.array(result)

def run_per_video(anno, resolution):
    """
    resolution: wh
    return: [xywh], 每帧的 FGT
    """
    """START：确定在第一帧中的起点"""
    patch_gt_xywh_start_in_original_img = get_init_patch_box_in_original_img(anno[0])
    x1_start, y1_start, patch_w, patch_h = patch_gt_xywh_start_in_original_img
    x2_start = x1_start + patch_w
    y2_start = y1_start + patch_h
    """END：确定在第一帧中的起点"""

    """START：确定补丁与物体的相对位置关系"""
    im_w, im_h = resolution
    direction = None
    if patch_w >= patch_h and y1_start >= im_h - y2_start:
        direction = 'up'
    elif patch_w >= patch_h and y1_start < im_h - y2_start:
        direction = 'down'
    elif patch_w < patch_h and x1_start >= im_w - x2_start:
        direction = 'left'
    elif patch_w < patch_h and x1_start < im_w - x2_start:
        direction = 'right'
    else:
        assert patch_gt_xywh_start_in_original_img
    delta = 16
    """END：确定补丁与物体的相对位置关系"""

    patch_anno_x1 = []
    patch_anno_y1 = []
    patch_anno_w = np.linspace(patch_w, 64, len(anno))
    patch_anno_h = np.linspace(patch_h, 64, len(anno))
    for i, gt in enumerate(anno):
        x1, y1, w, h = gt
        x2 = x1 + w
        y2 = y1 + h
        if direction == 'up':
            fgt_x1y1 = [x1, y1-delta-patch_anno_h[i]]
        elif direction == 'down':
            fgt_x1y1 = [x1, y2+delta]
        elif direction == 'left':
            fgt_x1y1 = [x1-delta-patch_anno_w[i], y1]
        elif direction == 'right':
            fgt_x1y1 = [x2+delta, y1]
        else:
            assert False, direction
        assert fgt_x1y1 is not None
        patch_anno_x1.append(fgt_x1y1[0])
        patch_anno_y1.append(fgt_x1y1[1])
    patch_anno_x1 = np.array(patch_anno_x1)
    patch_anno_y1 = np.array(patch_anno_y1)

    patch_anno = np.stack((patch_anno_x1, patch_anno_y1, patch_anno_w, patch_anno_h)).T

    return np.array(patch_anno)

def generate_patch_anno(save_root, visualize_flag):
    """"""
    """START：生成指定数据库"""
    if dataset_name == 'GOT-10k_Val':
        dataset = GOT10k('/home/zhbli/Dataset/data2/got10k', subset='val', return_meta=False)
    elif dataset_name == 'OTB_2015':
        dataset = OTB('/home/etvuz/projects/adversarial_attack/video_analyst/datasets/OTB/OTB2015',
                      version=2015, download=False)
    elif dataset_name == 'LaSOT':
        dataset = LaSOT('/home/zhbli/Dataset/data3/LASOT', subset='test', return_meta=False)
    else:
        assert False, dataset_name
    """END：生成指定数据库"""

    for s, (img_paths, anno) in enumerate(dataset):
        # anno: xywh
        resolution = Image.open(img_paths[0]).size  # wh
        patch_anno = run_per_video(anno, resolution)  # xywh

        """START：将补丁中心点保存为 txt"""
        save_dir = os.path.join(save_root, dataset_name)
        if not os.path.exists(save_dir):
            os.makedirs(save_dir)
        seq_name = dataset.seq_names[s]
        save_path = os.path.join(save_dir, seq_name + '.txt')
        print(save_path)
        np.savetxt(save_path, patch_anno, fmt='%d', delimiter=',')
        # 对 OTB2015 数据集 的 Jogging 视频进行特殊处理
        if seq_name == 'Jogging.1.txt':
            np.savetxt(os.path.join(save_dir, 'Jogging.txt'), patch_anno, fmt='%d', delimiter=',')
        if seq_name == 'Skating2.1.txt':
            np.savetxt(os.path.join(save_dir, 'Skating2.txt'), patch_anno, fmt='%d', delimiter=',')
        """END：将补丁中心点保存为 txt"""

        """START：可视化"""
        if visualize_flag:
            visualize(anno, patch_anno, resolution, s+1)
        """END：可视化"""

def test_get_crop():
    scale_x = get_crop(
        target_sz=(256, 256),  # 原始边框宽高
        z_size=127,
        context_amount=0.5,
    )
    return scale_x

if __name__ == '__main__':
    dataset_name = 'LaSOT'  # 'OTB_2015'
    generate_patch_anno(save_root='/home/etvuz/projects/adversarial_attack/patch_anno', visualize_flag=False)
    # test_get_crop()