text:
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: all,-execution,-papermill,-trusted
#     formats: ipynb,py//py:percent
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.7.1
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Description

# %% [markdown]
# This notebook generates a heatmap using a subset of the $\hat{\mathbf{M}}$ matrix (projectiong of S-MultiXcan results into latent gene expression representation). Traits in columns (organized by clusters of traits) and LVs in rows. The only difference with `05-global_clustermap.ipynb` is that here LVs shown are those detected to be driving clusters _and_ *NOT* well-aligned (FDR < 0.05) with any pathways (from the MultiPLIER models), thus representing potentially novel mechanisms.

# %% [markdown]
# # Environment variables

# %%
from IPython.display import display

import conf

N_JOBS = min((conf.GENERAL["N_JOBS"], 3))
display(N_JOBS)

# %%
# %env MKL_NUM_THREADS=$N_JOBS
# %env OPEN_BLAS_NUM_THREADS=$N_JOBS
# %env NUMEXPR_NUM_THREADS=$N_JOBS
# %env OMP_NUM_THREADS=$N_JOBS

# %% [markdown]
# # Modules loading

# %%
# %load_ext autoreload
# %autoreload 2

# %%
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

import conf
from entity import Trait
from data.cache import read_data

# %% [markdown]
# # Settings

# %%
OUTPUT_FIGURES_DIR = Path(conf.MANUSCRIPT["FIGURES_DIR"], "clustering").resolve()
display(OUTPUT_FIGURES_DIR)
OUTPUT_FIGURES_DIR.mkdir(parents=True, exist_ok=True)

# %% [markdown]
# # Data loading

# %% [markdown]
# ## MultiPLIER summary

# %%
multiplier_model_summary = read_data(conf.MULTIPLIER["MODEL_SUMMARY_FILE"])

# %%
multiplier_model_summary.shape

# %%
multiplier_model_summary.head()

# %%
well_aligned_lvs = multiplier_model_summary[
    (multiplier_model_summary["FDR"] < 0.05) | (multiplier_model_summary["AUC"] >= 0.75)
]

display(well_aligned_lvs.shape)
display(well_aligned_lvs.head())

# %% [markdown]
# ## S-MultiXcan results

# %%
smultixcan_results_filename = conf.PHENOMEXCAN[
    "SMULTIXCAN_EFO_PARTIAL_MASHR_ZSCORES_FILE"
]
display(smultixcan_results_filename)

# %%
smultixcan_results = pd.read_pickle(smultixcan_results_filename)

# %%
smultixcan_results.shape

# %%
smultixcan_results.head()

# %% [markdown]
# ## S-MultiXcan projection (`z_score_std`)

# %% tags=[]
INPUT_SUBSET = "z_score_std"

# %% tags=[]
INPUT_STEM = "projection-smultixcan-efo_partial-mashr-zscores"

# %% tags=[]
input_filepath = Path(
    conf.RESULTS["DATA_TRANSFORMATIONS_DIR"],
    INPUT_SUBSET,
    f"{INPUT_SUBSET}-{INPUT_STEM}.pkl",
).resolve()
display(input_filepath)

assert input_filepath.exists(), "Input file does not exist"

input_filepath_stem = input_filepath.stem
display(input_filepath_stem)

# %% tags=[]
data = read_data(input_filepath)

# %% tags=[]
data.shape

# %% tags=[]
data.head()

# %% [markdown]
# ## Clustering results

# %%
CONSENSUS_CLUSTERING_DIR = Path(
    conf.RESULTS["CLUSTERING_DIR"], "consensus_clustering"
).resolve()

display(CONSENSUS_CLUSTERING_DIR)

# %%
input_file = Path(CONSENSUS_CLUSTERING_DIR, "best_partitions_by_k.pkl").resolve()
display(input_file)

# %%
best_partitions = pd.read_pickle(input_file)

# %%
# keep selected partitions only
best_partitions = best_partitions[best_partitions["selected"]]

# %%
best_partitions.shape

# %%
best_partitions.head()

# %% [markdown]
# # Get all clusters information

# %%
partition_k = 29

# %%
part = best_partitions.loc[partition_k, "partition"]

# %%
part.shape

# %%
part_cluster_size = pd.Series(part).value_counts()
display(part_cluster_size)

# %%
clusters_labels = {
    12: "Asthma",
    27: "Well-being",
    22: "Nutrients",
    15: "Heel bone",
    4: "Reticulocytes",
    2: "Red blood volumes",
    5: "Erythrocytes",
    23: "Read blood DW",
    20: "BMI",
    18: "Spirometry",
    19: "Height",
    10: "Keratometry",
    1: "Platelets",
    6: "Skin/hair color",
    13: "Autoimmune #13",
    26: "Autoimmune #26",
    8: "Autoimmune #8",
    17: "HTN/high cholesterol",
    25: "Other traits #25",
    21: "Eczema/IBD/SCZ",
    28: "Env. factors #28",
    11: "CAD/breast cancer",
    16: "Lipids/AD/Chronotype",
    14: "Other heart diseases",
    7: "Monocytes",
    24: "Lymphocytes",
    9: "Neutrophils",
    3: "Eosinophils",
}

# %%
all_clusters = []

for c_idx in part_cluster_size[part_cluster_size < 50].index:
    _df = pd.DataFrame(
        {
            "trait": data.index[part == c_idx].tolist(),
            "cluster_label": clusters_labels[c_idx],
            "cluster_idx": c_idx,
        }
    )

    all_clusters.append(_df)

# %%
# combine all
df = pd.concat(all_clusters, ignore_index=True)

# %%
df.shape

# %%
df.head()

# %% [markdown]
# # Get top LVs for each cluster

# %%
CLUSTER_LV_DIR = (
    conf.RESULTS["CLUSTERING_INTERPRETATION"]["BASE_DIR"]
    / "cluster_lvs"
    / f"part{partition_k}"
)
display(CLUSTER_LV_DIR)

# %%
def _get_lvs_data(cluster_idx):
    cluster_lvs = pd.read_pickle(
        CLUSTER_LV_DIR / f"cluster_interpreter-part{partition_k}_k{cluster_idx}.pkl"
    )

    # keep "novel" lvs only
    cluster_lvs = cluster_lvs[
        ~cluster_lvs.index.astype(str).isin(well_aligned_lvs["LV index"])
    ]

    cluster_lvs["cluster"] = cluster_idx

    return pd.Series([set(cluster_lvs["name"]), cluster_lvs])

# %%
clusters_df = df.groupby("cluster_idx")[["cluster_label", "cluster_idx"]].head(1)

# add top lvs for each cluster
clusters_df[["lvs_names", "lvs_data"]] = clusters_df["cluster_idx"].apply(_get_lvs_data)

clusters_df = clusters_df.set_index("cluster_label")
assert clusters_df.index.is_unique

display(clusters_df)

# %%
lvs_data = (
    pd.concat(clusters_df["lvs_data"].tolist())
    .reset_index()
    .set_index(["cluster", "idx"])
)

# %%
lvs_data.shape

# %%
lvs_data.head(10)

# %%
_traits = df["trait"]
_lvs = set(lvs_data["name"])

# %%
_lvs

# %%
data_subset = data.loc[_traits, _lvs].rename_axis("trait")
data_subset = data_subset.assign(
    cluster_name=df.set_index("trait")["cluster_label"].astype("category")
)

# %%
data_subset.shape

# %%
data_subset

# %%
data_subset["trait_description"] = data_subset.apply(
    lambda x: Trait.get_trait(full_code=x.name).description
    if not Trait.is_efo_label(x.name)
    else x.name,
    axis=1,
)

# %%
data_subset

# %%
data_subset = data_subset.replace(
    {
        "trait_description": {
            "platelet count": "Platelet count",
            "Red blood cell (erythrocyte) count": "Erythrocyte count",
            "erythrocyte count": "Erythrocyte count",
            "Red blood cell (erythrocyte) distribution width": "Erythrocyte distribution width",
            "reticulocyte count": "Reticulocyte count",
        }
    }
)

# %%
_cluster_column = data_subset.set_index("trait_description")["cluster_name"].to_dict()

# %%
_cluster_column

# %%
data_subset

# %%
# calculate mean of traits with the same description
_tmp = data_subset.groupby("trait_description").mean()
_tmp = _tmp.assign(cluster_name=_tmp.apply(lambda x: _cluster_column[x.name], axis=1))

# %%
_tmp.shape

# %%
data_subset = _tmp.sort_values("cluster_name")

# %%
data_subset.shape

# %%
data_subset.head()

# %%
data_subset.describe()

# %%
data_subset["cluster_name"].unique()

# %%
# this is the order from the clustering tree
cluster_order = [
    "Asthma",
    #     "Well-being",
    #     "Nutrients",
    "Heel bone",
    "Reticulocytes",
    "Red blood volumes",
    "Erythrocytes",
    "Read blood DW",
    "BMI",
    "Spirometry",
    "Height",
    "Keratometry",
    "Platelets",
    "Skin/hair color",
    "Autoimmune #13",
    "Autoimmune #26",
    "Autoimmune #8",
    "HTN/high cholesterol",
    #     "Other traits #25",
    "Eczema/IBD/SCZ",
    #     "Other traits #28",
    "CAD/breast cancer",
    "Lipids/AD/Chronotype",
    "Other heart diseases",
    "Monocytes",
    "Lymphocytes",
    "Neutrophils",
    "Eosinophils",
]

# %%
# _cluster = data_subset["cluster_name"].unique()
lut = dict(zip(cluster_order, sns.color_palette("Paired", len(cluster_order))))

# %%
row_colors = data_subset[["cluster_name"]].assign(
    color=data_subset["cluster_name"].map(lut)
)

# %%
row_colors.head(20)

# %% [markdown]
# # Plot data

# %%
_tmp = data_subset.reset_index().set_index(["cluster_name", "trait_description"])

# %%
plot_data = _tmp.loc[cluster_order].droplevel("cluster_name")

# %%
plot_data = plot_data.assign(
    color=row_colors["color"], cluster_name=row_colors["cluster_name"]
)

# %%
plot_data.shape

# %%
plot_data.head()

# %%
_tmp = plot_data.groupby("cluster_name").count().iloc[:, 0]

# %%
_tmp = _tmp.loc[cluster_order]

# %%
_tmp

# %%
_tmp.cumsum()

# %%
cluster_order_start_idx = np.array([0.5] + _tmp.cumsum()[:-1].tolist())
# cluster_order = _tmp.index.tolist()

# %%
cluster_order_start_idx

# %%
# slightly change the "Blood DW" and "Reticulocytes" index since they overlap a bit in the figure
cluster_label_order_start_idx = cluster_order_start_idx.copy()
# cluster_label_order_start_idx[10] = cluster_label_order_start_idx[10] - 0.2
cluster_label_order_start_idx[4] = cluster_label_order_start_idx[4] - 1.0
cluster_label_order_start_idx[6] = cluster_label_order_start_idx[6] + 2.0

cluster_label_order_start_idx[13] = cluster_label_order_start_idx[13] - 1.0
cluster_label_order_start_idx[15] = cluster_label_order_start_idx[15] + 2.0

# %%
cluster_label_order_start_idx

# %%
cluster_order_start_idx.shape

# %%
len(cluster_order)

# %%
pd.Series(plot_data.drop(columns=["cluster_name", "color"]).values.flatten()).quantile(
    [0.90, 0.95, 0.97, 0.98, 0.99]
)

# %%
pd.Series(plot_data.drop(columns=["cluster_name", "color"]).values.flatten()).describe()

# %% [markdown]
# # Plotting functions

# %%
def _get_lv_pathway(lv_name):
    p = (
        well_aligned_lvs_df.loc[lv_name]
        .sort_values("FDR")
        .index.get_level_values("pathway")
    )
    return ", ".join(p)

def _get_lv_label(lv_name, n_top=2):
    lv_info = all_results[lv_name]

    cell_types_label = ""
    tissues_label = ""

    if "attr" in lv_info:
        cell_types_label = ", ".join(lv_info["attr"][:n_top])

    if "tissue" in lv_info:
        tissues_label = ", ".join(lv_info["tissue"][:n_top])

    final_label = ""

    if cell_types_label is not None:
        final_label += cell_types_label

    if tissues_label is not None:
        final_label += f" | {tissues_label}"

    return final_label

# %%
def plot_clustermap(rc, lv_labeling_function, output_filepath):
    with sns.plotting_context("paper", font_scale=1.25, rc=rc):
        g = sns.clustermap(
            data=plot_data.drop(columns=["cluster_name", "color"]).T,
            col_colors=plot_data["color"].rename("Cluster"),
            vmin=0.0,
            vmax=7.0,
            row_cluster=True,
            col_cluster=False,
            cmap="YlOrBr",
            dendrogram_ratio=0.05,
            cbar_pos=None,  # (-0.06, 0.8, 0.05, 0.18),
            xticklabels=False,
            yticklabels=True,
            figsize=(18, 38),
            method="complete",
            linewidths=0.10,
            rasterized=True,
            colors_ratio=0.01,
        )

        g.ax_heatmap.set_xlabel("")
        g.ax_heatmap.set_ylabel("")

        g.ax_col_colors.set_xticks(cluster_label_order_start_idx)
        g.ax_col_colors.set_xticklabels(
            cluster_order,
            rotation=45,
            horizontalalignment="left",
        )
        g.ax_col_colors.xaxis.set_tick_params(size=0)
        g.ax_col_colors.xaxis.tick_top()

        g.ax_heatmap.get_xaxis().set_ticklabels(
            g.ax_heatmap.get_xaxis().get_ticklabels(),
            rotation=45,
            horizontalalignment="right",
        )

        if lv_labeling_function is not None:
            new_y_labels = []
            for t in g.ax_heatmap.get_yaxis().get_ticklabels():
                t.set_text(f"{t.get_text()}: {lv_labeling_function(t.get_text())}")
                t.set_fontsize(0.75 * t.get_fontsize())
                new_y_labels.append(t)

            g.ax_heatmap.get_yaxis().set_ticklabels(
                new_y_labels,
            )

        g.ax_heatmap.vlines(
            cluster_order_start_idx[1:], *g.ax_heatmap.get_ylim(), colors="black"
        )

        display(output_filepath)
        plt.savefig(
            output_filepath,
            dpi=600,
            bbox_inches="tight",
            facecolor="white",
        )

# %% [markdown]
# # Simple plot

# %%
rc = {
    "axes.labelsize": 14,
    "xtick.labelsize": 14,
    "ytick.labelsize": 8,
}

plot_clustermap(rc, None, OUTPUT_FIGURES_DIR / "global_clustermap-novel-plain.svg")

# %%

max_stars_repo_name: greenelab/phenoplier

max_issues_repo_name: greenelab/phenoplier

max_forks_repo_name: greenelab/phenoplier


text:
import pandas as pd
import numpy as np
import tensorflow as tf
from sklearn.model_selection import train_test_split
import tensorflow_datasets as tfds
import os
from PIL import Image
from tqdm import tqdm
from pdb import set_trace as st

def load_CH_MNIST():
    """
    Loads CH_MNIST dataset and maps it to Target Model and Shadow Model.
    :param model_mode: one of "TargetModel" and "ShadowModel".
    :return: Tuple of numpy arrays:'(x_train, y_train), (x_test, y_test), member'.
    :raise: ValueError: in case of invalid `model_mode`.
    """

    # Initialize Data
    images, labels = tfds.load('colorectal_histology', split='train', batch_size=-1, as_supervised=True)

    x_train, x_test, y_train, y_test = train_test_split(images.numpy(), labels.numpy(), train_size=0.8,
                                                        random_state=1,
                                                        stratify=labels.numpy())

    x_train = tf.image.resize(x_train, (64, 64))
    y_train = tf.keras.utils.to_categorical(y_train-1, num_classes=8)
    m_train = np.ones(y_train.shape[0])

    x_test = tf.image.resize(x_test, (64, 64))
    y_test = tf.keras.utils.to_categorical(y_test-1, num_classes=8)
    m_test = np.zeros(y_test.shape[0])
    st()
    return (x_train, y_train), (x_test, y_test)

load_CH_MNIST()

ext: py

max_stars_repo_name: ziqi-zhang/FedML

max_issues_repo_name: ziqi-zhang/FedML

max_forks_repo_name: ziqi-zhang/FedML

reason: import numpy


text:
"""This module contains simple helper functions """
from __future__ import print_function
import torch
import numpy as np
from PIL import Image
import os
from typing import Any, List, Tuple, Union
import sys
import random

def set_seed(seed=None):
    if seed is not None:
        torch.manual_seed(seed)
        np.random.seed(seed)
        torch.backends.cudnn.deterministic = True
        random.seed(seed)
        torch.cuda.manual_seed_all(seed)
        torch.backends.cudnn.benchmark = True

class Logger(object):
    """Redirect stderr to stdout, optionally print stdout to a file, and optionally force flushing on both stdout and the file."""

    def __init__(self, file_name: str = None, file_mode: str = "a", should_flush: bool = True, append=False):
        self.file = None

        if append:
            file_mode = 'a'
        else:
            file_mode = 'w'
        self.file = open(file_name, file_mode)

        self.should_flush = should_flush
        self.stdout = sys.stdout
        self.stderr = sys.stderr

        sys.stdout = self
        sys.stderr = self

    def __enter__(self) -> "Logger":
        return self

    def __exit__(self, exc_type: Any, exc_value: Any, traceback: Any) -> None:
        self.close()

    def write(self, text: str) -> None:
        """Write text to stdout (and a file) and optionally flush."""
        if len(text) == 0: # workaround for a bug in VSCode debugger: sys.stdout.write(''); sys.stdout.flush() => crash
            return

        if self.file is not None:
            self.file.write(text)

        self.stdout.write(text)

        if self.should_flush:
            self.flush()

    def flush(self) -> None:
        """Flush written text to both stdout and a file, if open."""
        if self.file is not None:
            self.file.flush()

        self.stdout.flush()

    def close(self) -> None:
        """Flush, close possible files, and remove stdout/stderr mirroring."""
        self.flush()

        # if using multiple loggers, prevent closing in wrong order
        if sys.stdout is self:
            sys.stdout = self.stdout
        if sys.stderr is self:
            sys.stderr = self.stderr

        if self.file is not None:
            self.file.close()

def write_loss(iterations, trainer, train_writer, prefix):
    members = [attr for attr in dir(trainer) \
               if not callable(getattr(trainer, attr)) and not attr.startswith("__") and (
                       'loss' in attr or 'grad' in attr or 'nwd' in attr ) and 'name' not in attr and 'pool' not in attr]
    for m in members:
        train_writer.add_scalar(prefix+'/'+m, getattr(trainer, m), iterations + 1)

def format_time(seconds: Union[int, float]) -> str:
    """Convert the seconds to human readable string with days, hours, minutes and seconds."""
    s = int(np.rint(seconds))

    if s < 60:
        return "{0}s".format(s)
    elif s < 60 * 60:
        return "{0}m {1:02}s".format(s // 60, s % 60)
    elif s < 24 * 60 * 60:
        return "{0}h {1:02}m {2:02}s".format(s // (60 * 60), (s // 60) % 60, s % 60)
    else:
        return "{0}d {1:02}h {2:02}m".format(s // (24 * 60 * 60), (s // (60 * 60)) % 24, (s // 60) % 60)

def tensor2im(input_image, imtype=np.uint8):
    """"Converts a Tensor array into a numpy image array.

    Parameters:
        input_image (tensor) --  the input image tensor array
        imtype (type)        --  the desired type of the converted numpy array
    """
    if not isinstance(input_image, np.ndarray):
        if isinstance(input_image, torch.Tensor):  # get the data from a variable
            image_tensor = input_image.data
        else:
            return input_image
        image_numpy = image_tensor[0].cpu().float().numpy()  # convert it into a numpy array
        if image_numpy.shape[0] == 1:  # grayscale to RGB
            image_numpy = np.tile(image_numpy, (3, 1, 1))
        image_numpy = (np.transpose(image_numpy, (1, 2, 0)) + 1) / 2.0 * 255.0  # post-processing: tranpose and scaling
    else:  # if it is a numpy array, do nothing
        image_numpy = input_image
    return image_numpy.astype(imtype)

def str2bool(x):
    return x.lower() in ('true')

def diagnose_network(net, name='network'):
    """Calculate and print the mean of average absolute(gradients)

    Parameters:
        net (torch network) -- Torch network
        name (str) -- the name of the network
    """
    mean = 0.0
    count = 0
    for param in net.parameters():
        if param.grad is not None:
            mean += torch.mean(torch.abs(param.grad.data))
            count += 1
    if count > 0:
        mean = mean / count
    print(name)
    print(mean)

def save_image(image_numpy, image_path, aspect_ratio=1.0):
    """Save a numpy image to the disk

    Parameters:
        image_numpy (numpy array) -- input numpy array
        image_path (str)          -- the path of the image
    """

    image_pil = Image.fromarray(image_numpy)
    h, w, _ = image_numpy.shape

    if aspect_ratio > 1.0:
        image_pil = image_pil.resize((h, int(w * aspect_ratio)), Image.BICUBIC)
    if aspect_ratio < 1.0:
        image_pil = image_pil.resize((int(h / aspect_ratio), w), Image.BICUBIC)
    image_pil.save(image_path)

def print_numpy(x, val=True, shp=False):
    """Print the mean, min, max, median, std, and size of a numpy array

    Parameters:
        val (bool) -- if print the values of the numpy array
        shp (bool) -- if print the shape of the numpy array
    """
    x = x.astype(np.float64)
    if shp:
        print('shape,', x.shape)
    if val:
        x = x.flatten()
        print('mean = %3.3f, min = %3.3f, max = %3.3f, median = %3.3f, std=%3.3f' % (
            np.mean(x), np.min(x), np.max(x), np.median(x), np.std(x)))

def get_model_list(dirname, key, exclude='latest'):
    if os.path.exists(dirname) is False:
        return None
    gen_models = [os.path.join(dirname, f) for f in os.listdir(dirname) if
                  os.path.isfile(os.path.join(dirname, f)) and key in f and ".pt" in f and exclude not in f]
    if gen_models is None:
        return None
    gen_models.sort()
    last_model_name = gen_models[-1]
    return last_model_name

def mkdirs(paths):
    """create empty directories if they don't exist

    Parameters:
        paths (str list) -- a list of directory paths
    """
    if isinstance(paths, list) and not isinstance(paths, str):
        for path in paths:
            mkdir(path)
    else:
        mkdir(paths)

def mkdir(path):
    """create a single empty directory if it didn't exist

    Parameters:
        path (str) -- a single directory path
    """
    if not os.path.exists(path):
        os.makedirs(path)

def adjust_dynamic_range(data, drange_in, drange_out):
    if drange_in != drange_out:
        scale = (np.float32(drange_out[1]) - np.float32(drange_out[0])) / (np.float32(drange_in[1]) - np.float32(drange_in[0]))
        bias = (np.float32(drange_out[0]) - np.float32(drange_in[0]) * scale)
        data = data * scale + bias
    return data

def create_image_grid(images, grid_size=None):
    assert images.ndim == 3 or images.ndim == 4
    num, img_w, img_h = images.shape[0], images.shape[-1], images.shape[-2]

    if grid_size is not None:
        grid_w, grid_h = tuple(grid_size)
    else:
        grid_w = max(int(np.ceil(np.sqrt(num))), 1)
        grid_h = max((num - 1) // grid_w + 1, 1)

    grid = np.zeros(list(images.shape[1:-2]) + [grid_h * img_h, grid_w * img_w], dtype=images.dtype)
    for idx in range(num):
        x = (idx % grid_w) * img_w
        y = (idx // grid_w) * img_h
        grid[..., y : y + img_h, x : x + img_w] = images[idx]
    return grid

def convert_to_pil_image(image, drange=[-1,1]):
    assert image.ndim == 2 or image.ndim == 3
    if image.ndim == 3:
        if image.shape[0] == 1:
            image = image[0] # grayscale CHW => HW
        elif image.shape[1]>image.shape[0]:
            image = image.transpose(1, 2, 0) # CHW -> HWC

    image = adjust_dynamic_range(image, drange, [0,255])
    image = np.rint(image).clip(0, 255).astype(np.uint8)
    fmt = 'RGB' if image.ndim == 3 else 'L'
    return Image.fromarray(image, fmt)

def save_images(image, filename, drange=[-1,1], quality=95):
    img = convert_to_pil_image(image, drange)
    if '.jpg' in filename:
        img.save(filename, "JPEG", quality=quality, optimize=True)
    else:
        img.save(filename)
def to_var( x):
    """Converts numpy to variable."""
    if torch.cuda.is_available():
        x = x.cuda()
    return torch.autograd.Variable(x)

def to_data(x):
    """Converts variable to numpy."""
    if torch.cuda.is_available():
        x = x.cpu()
    return x.data.numpy()
def save_image_grid(images, filename, drange=[-1,1], grid_size=None):
    convert_to_pil_image(create_image_grid(images, grid_size), drange).save(filename)

max_stars_repo_name: Mid-Push/IrwGAN

max_issues_repo_name: Mid-Push/IrwGAN

max_forks_repo_name: Mid-Push/IrwGAN


text:
#!/usr/bin/env python3
'''
Primjer pokretanja:
-sam gradi poravnanja od parova gena
python3 matrix_generator -g gene_pairs.txt 10

-poravnanja predana u datoteci
python3 matrix_generator -a algnmt_pairs.txt
'''
import sys
import numpy as np
from Bio import Align
__author__ = "Stjepan Jakovac"
__credits__ = ["Jurij Kos", "Frane Polić"]
__license__ = "MIT"

#emiss_i = {'-':0, 'A':1, 'B':2, 'C':3}  #mapira baze na indexe za matrice
emiss_i = {'-':0, 'A':1, 'C':2, 'G':3, 'T':4}
emiss_matrix = np.zeros((len(emiss_i),len(emiss_i)))
trans_matrix = np.zeros((5,5)) # 5 jer start, gap_x, gap_y, miss/match, end

def generate_empty_matrices():
	return np.zeros((len(emiss_i),len(emiss_i))), np.zeros((5,5))  #emiss, trans

def cleaned_fasta(file_name, chars_to_remove = ('x','N', 'Y')):
	'''
	Iz datoteke ucita genom, izbaci iz nj. s pojave znakov iz tupla chars_to_remove.
	Vrati string reprezentaciju ociscenog genoma.
	'''
	genome = None
	with open("../data/"+file_name, 'r') as source:
		genome = source.readline()
	replace_table = str.maketrans({c:None for c in chars_to_remove})

	genome = genome.translate(replace_table)
	return genome.strip()

def populate_from_pair(g_1, g_2, trans_m, emiss_m, N = 10):
	'''
	Popunjava trans i emis matricu iz najvise N optimalnih poravnanja koja se dobivaju
	od genoma g1 i g2. Tablicu popunjavaju pojavama(count) te se ne radi pretovrba u vjerojatnosu
	matricu
	'''
	aligner = Align.PairwiseAligner()
	alignments = aligner.align(g_1, g_2)

	for i,alignment in enumerate(alignments):
		#print(alignment)
		if N is not None and i>=N:
			break
		populate_from_aligment(alignment, trans_m, emiss_m)

def populate_from_aligment(alignment, trans_m, emiss_m, emiss_i = emiss_i):
	'''
	Iz jednog poravnanja popunjuje trans i emiss matricu
	'''
	s_algmnt = alignment.format()
	algmnt_len = len(s_algmnt)//3-1
	prev_state = 0 #start | gap_a/emit_y | gap_b/emit_x | m (match/miss) | end
	state = None
	for i in range(algmnt_len):
		c1 = s_algmnt[i]
		c2 = s_algmnt[i+(algmnt_len+1)*2]

		if c1 == '-':
			state = 1
		elif c2 == '-':
			state = 2
		elif c2 != '\n':
			state = 3
		trans_m[prev_state][state] += 1
		prev_state = state

		i1 = emiss_i[c1]; i2 = emiss_i[c2]
		emiss_m[i1][i2] += 1
	trans_m[prev_state][4] += 1

def populate_from_raw_aligment(raw_1, raw_2, trans_m, emiss_m, emiss_i = emiss_i):
	'''
	Iz dva pojedinacna string poravnanja popunjava matrice transmisje i emisije
	'''
	prev_state = 0 #start | gap_a | gap_b | m (match/miss) | end
	state = None
	for c1, c2 in zip(raw_1, raw_2):
		if c1 == '-' and c2 == '-':
			continue

		if c1 == '-':
			state = 1
		elif c2 == '-':
			state = 2
		else:
			state = 3
		trans_m[prev_state][state] += 1
		prev_state = state

		i1 = emiss_i[c1]; i2 = emiss_i[c2]
		emiss_m[i1][i2] += 1
	trans_m[prev_state][4] += 1

def proba_from_counts(trans_m, emiss_m):
	'''
	Matrice u kojima je biljezen broj pojava transformia u vjerojatnosne matrice.
	(in-place zasada)
	'''
	row_sums = trans_m.sum(axis=1)
	row_sums[4] = 1 #da zadnji ne bude ZeroDiv (zadnji su P(#nesto|end_state))
	trans_m /= row_sums[:,None]

	gap_a_sum = emiss_matrix[0,:].sum()
	gap_b_sum = emiss_matrix[:,0].sum()
	mm_sum = emiss_matrix.sum() - gap_a_sum - gap_b_sum
	emiss_matrix[0,:] /= gap_a_sum
	emiss_matrix[:,0] /= gap_b_sum
	emiss_matrix[1:,1:] /= mm_sum

mode = 0 if sys.argv[1] == '-g' else 1
pairs = []
config_file_name = sys.argv[2]
with open(config_file_name, 'r') as config_file:
	for line in config_file.readlines():
		pairs.append(line.split())

#iz datoteke zadane preko komandne linije ucitava parove iz kojih ce se graditi matrice
if mode == 0:
	N = int(sys.argv[3])
	#iterira po parovima i popunjava tablice
	for fn_1, fn_2 in pairs:
		g_1 = cleaned_fasta(fn_1)
		g_2 = cleaned_fasta(fn_2)
		populate_from_pair(g_1, g_2, trans_matrix, emiss_matrix, N = N)

#popunjava tablice iz gotovog poravnanja
else:
	for fn_1, fn_2 in pairs:
		g_1 = cleaned_fasta(fn_1)
		g_2 = cleaned_fasta(fn_2)
		populate_from_raw_aligment(g_1, g_2, trans_matrix, emiss_matrix)

'''
print(trans_matrix)
print()
print(emiss_matrix)
'''
proba_from_counts(trans_matrix, emiss_matrix)

formating = {
				'prefix' : '+', \
				'separator' : ' ', \
				'formatter' : {'float' : (lambda x : ' {:.6f}'.format(x))}
			}

tm_file = open("../data/tran_mat.txt", 'w')
tm_file.write(np.array2string(trans_matrix, **formating).replace('[', '').replace(']', ''))
tm_file.close()

em_file = open("../data/emis_mat.txt", 'w')
em_file.write(np.array2string(emiss_matrix, **formating).replace('[', '').replace(']', ''))
em_file.close()

max_stars_repo_name: jurijkos/bioinformatika

max_issues_repo_name: jurijkos/bioinformatika

max_forks_repo_name: jurijkos/bioinformatika


text:
import os
import math
import numpy
import pyopencl as cl
import pyopencl.array
import time
from PIL import Image

TASKS = 1048576
CL_TASKS = int(TASKS / 4)

if __name__ == '__main__':
    print('load program from cl source file')
    f = open('blur.cl', 'r', encoding='utf-8')
    kernels = ''.join(f.readlines())
    f.close()

    print('prepare data ... ')
    filename = '7680x4320'
    ext = '.jpg'
    mask_size = 15
    img = Image.open(os.path.join(os.path.dirname(__file__), filename + ext))
    if img.mode != 'RGBA':
        img = img.convert('RGBA')
    img_width = img.size[0]
    img_height = img.size[1]
    img_size = img_width * img_height

    start_time = time.time()
    # prepare host memory for OpenCL
    img_bytes = img.tobytes()

    time_hostdata_loaded = time.time()

    print('create context ...')
    ctx = cl.create_some_context()
    print('create command queue ...')
    queue = cl.CommandQueue(ctx, properties=cl.command_queue_properties.PROFILING_ENABLE)
    time_ctx_queue_creation = time.time()

    dev_image_format = cl.ImageFormat(cl.channel_order.RGBA,
                                      cl.channel_type.UNSIGNED_INT8)
    # prepare device memory for OpenCL
    print('prepare device memory for input / output')
    input_image = cl.Image(ctx,
                           cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR,
                           dev_image_format,
                           img.size,
                           None,
                           img_bytes)
    output_image = cl.Image(ctx,
                            cl.mem_flags.WRITE_ONLY,
                            dev_image_format,
                            img.size)
    time_devicedata_loaded = time.time()

    print('compile kernel code')
    prg = cl.Program(ctx, kernels).build()
    time_kernel_compilation = time.time()

    np_masksize = numpy.int32(mask_size)
    print('execute kernel programs')
    evt = prg.to_blur(queue, (img_size, ), (1, ),
                      np_masksize,
                      input_image, output_image)
    print('wait for kernel executions')
    evt.wait()
    elapsed = 1e-9 * (evt.profile.end - evt.profile.start)

    time_before_readback = time.time()
    buffer = numpy.zeros(img_width * img_height * 4, numpy.uint8)
    origin = (0, 0, 0)
    region = (img_width, img_height, 1)
    cl.enqueue_read_image(queue, output_image,
                          origin, region, buffer).wait()
    time_after_readback = time.time()

    print('Prepare host data took       : {}'.format(time_hostdata_loaded - start_time))
    print('Create CTX/QUEUE took        : {}'.format(time_ctx_queue_creation-time_hostdata_loaded))
    print('Upload data to device took   : {}'.format(time_devicedata_loaded - time_ctx_queue_creation))
    print('Compile kernel took          : {}'.format(time_kernel_compilation-time_devicedata_loaded))
    print('OpenCL elapsed time          : {}'.format(elapsed))
    print('Offload data from device took: {}'.format(time_after_readback - time_before_readback))

    out_filename = os.path.join('out_' + filename + ext)
    out_im = Image.frombytes("RGBA", img.size, buffer.tobytes())
    if out_im.mode == 'RGBA':
        out_im = out_im.convert('RGB')
    out_im.save(out_filename)

    print('Results is OK')

max_stars_repo_name: sixAndjune/pyopencl-examples

max_issues_repo_name: sixAndjune/pyopencl-examples

max_forks_repo_name: sixAndjune/pyopencl-examples


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 10 10:46:48 2020

@author: RogelioTESI
"""
import numpy as np
import matplotlib.pyplot as plt
Em=10
Ec=20
fm=5e3
fc=100e3
Tm=1/fm
t=np.linspace(0,2*Tm,num=2000)
Vm=Em*np.sin(2*np.pi*fm*t)
Vc=Ec*np.sin(2*np.pi*fc*t)
Vmod=(Ec+Em*np.sin(2*np.pi*fm*t))*np.sin(2*np.pi*fc*t)
fig,axes=plt.subplots()
axes.plot(t*1000,Vm)
axes.plot(t*1000,Vc)
axes.set_xlabel('Time(ms)')
axes.set_ylabel('Amplitud(V)')
plt.grid(True)
fig1,axes1=plt.subplots()
axes1.plot(t*1000,Vmod)
axes1.set_xlabel('Time(ms)')
axes1.set_ylabel('Amplitud(V)')
plt.grid(True)


text:
"""Based on SIAF transforms PATHS
NIRCAM
------
science --> ideal --> V2,V3
V2,V3 --> ideal --> science
"""
import sys

import numpy as np
from astropy.modeling.models import Polynomial2D, Shift

#t = ascii.read("NIRCam_SIAF_2016-09-29.csv",header_start=1)

def get_siaf_transform(row, aperture, from_system, to_system, degree):
    """
    This reads in the file with transformations that the TEL team
    is using to construct the SIAF file. These transformations are
    defined as two polynomials (one in x, and one in y) transforming coordinates
    from one coordinate system to another.
    Parameters
    ----------
    aperture: str
        Name of aperture on NIRCam, composed of the detector name followed
        by an underscore and the subarray name. (e.g. "NRCA1_FULL",
        "NRCB5_SUB160")
    from_system : str
        Starting system (e.g. "science", "ideal")
    to_system : str
        Ending coordinate system (e.g. "science" , "ideal")
    degree : int
        Degree of polynomial
    Returns
    -------
    x_model : astropy.modeling.Model
        Correction in x
    y_model : astropy.modeling.Model
        Correction in y
    from_units : str
        Units in the starting system
    to_units : str
        Units in the ending system
    Examples
    --------
    >>> get_siaf_transform('NRCA1_FULL', "science", "ideal", 5)
    """
    #read in csv file of coefficients
    #t = ascii.read("NIRCam_SIAF_2016-09-29.csv",header_start=1)
    #t = ascii.read(coefffile,header_start=1)

    #from_system and to_system are very limited. Can only be "ideal" for the
    #distortion-free coords, and "science" for distorted coords
    from_system = from_system.lower()
    if from_system not in ['ideal','science']:
        print("Requested from_system of {} not recognized.".format(from_system))
        sys.exit()

    to_system = to_system.lower()
    if to_system not in ['ideal','science']:
        print("Requested to_system of {} not recognized.".format(to_system))
        sys.exit()

    #Generate the string corresponding to the requested coefficient labels
    if from_system == 'ideal' and to_system == 'science':
        label = 'Idl2Sci'
        from_units = 'arcsec'
        to_units = 'distorted pixels'
    elif from_system == 'science' and to_system == 'ideal':
        label = 'Sci2Idl'
        from_units = 'distorted pixels'
        to_units = 'arcsec'

    #Find the row that matches the requested aperture
    #match = t['AperName'] == aperture
    #if np.any(match) == False:
    #    print("Aperture name {} not found in input CSV file.".format(aperture))
    #    sys.exit()
    #row = t[match]

    #Get the coefficients for "science" to "ideal" transformation (and back)
    #"science" is distorted pixels. "ideal" is undistorted arcsec from the
    #the reference pixel location.

    #we need the parity value, which describes the relationship between the
    #v2,v3 coordinate system and the science pixel coordinate system
    parity = row['VIdlParity'].data[0]

    #Then create the model for the transformation
    X_cols = [c for c in row.colnames if label+'X' in c]
    #for kw in X_cols:
    #    ele = row[kw]
    #    ele *= parity
    #    row[kw] = ele
    x_coeffs = row[X_cols]
    X_model = to_model(x_coeffs,degree)

    Y_cols = [c for c in row.colnames if label+'Y' in c]
    y_coeffs = row[Y_cols]
    Y_model = to_model(y_coeffs,degree)

    return X_model, Y_model, from_units, to_units

def to_model(coeffs, degree=5):
    """
    Creates an astropy.modeling.Model object
    Parameters
    ----------
    coeffs : array like
        Coefficients from the ISIM transformations file.
    degree : int
        Degree of polynomial.
        Default is 5 as in the ISIM file but many of the polynomials are of
        a smaller degree.
    Returns
    -------
    poly : astropy.modeling.Polynomial2D
        Polynomial model transforming one coordinate (x or y) between two systems.
    """

    #map Colin's coefficients into the order expected by Polynomial2D
    c = {}
    for cname in coeffs.colnames:
        siaf_i = int(cname[-2])
        siaf_j = int(cname[-1])
        name = 'c{0}_{1}'.format(siaf_i-siaf_j,siaf_j)
        c[name] = coeffs[cname].data[0]

    #0,0 coefficient should not be used, according to Colin's TR
    #JWST-STScI-001550
    c['c0_0'] = 0

    return Polynomial2D(degree, **c)

def get_siaf_v2v3_transform(row,aperture,from_system='v2v3',to_system='v2v3'):
    """
    Generate transformation model to go to/from V2/V3 from
    undistorted angular distnaces from the reference pixel ("ideal")
    """
    #read in csv file of coefficients
    #t = ascii.read("NIRCam_SIAF_2016-09-29.csv",header_start=1)
    #t = ascii.read(coefffile,header_start=1)

    from_system = from_system.lower()
    to_system = to_system.lower()

    if from_system != 'v2v3' and to_system != 'v2v3':
        print("WARNING, either from_system or to_system must be 'v2v3'")
        sys.exit()

    #Find the row that matches the requested aperture
    #match = t['AperName'] == aperture
    #if np.any(match) == False:
    #    print("Aperture name {} not found in input CSV file.".format(aperture))
    #    sys.exit()
    #row = t[match]

    #Then create the model for the transformation
    parity = row['VIdlParity'].data[0]
    v3_ideal_y_angle = row['V3IdlYAngle'].data[0] * np.pi / 180.

    #print("parity and angle are {}, {}".format(parity,v3_ideal_y_angle))

    X_model, Y_model = v2v3_model(from_system,to_system,parity,v3_ideal_y_angle)

    return X_model, Y_model

def v2v3_model(from_sys, to_sys, par, angle):
    """
    Creates an astropy.modeling.Model object
    for the undistorted ("ideal") to V2V3 coordinate translation
    """
    if from_sys != 'v2v3' and to_sys != 'v2v3':
        print("This function is designed to generate the transformation either to or from V2V3.")
        sys.exit()

    #cast the transform functions as 1st order polynomials
    xc = {}
    yc = {}
    if to_sys == 'v2v3':
        xc['c1_0'] = par * np.cos(angle)
        xc['c0_1'] = np.sin(angle)
        yc['c1_0'] = (0.-par) * np.sin(angle)
        yc['c0_1'] = np.cos(angle)

    if from_sys == 'v2v3':
        xc['c1_0'] = par * np.cos(angle)
        xc['c0_1'] = par * (0. - np.sin(angle))
        yc['c1_0'] = np.sin(angle)
        yc['c0_1'] = np.cos(angle)

    #0,0 coeff should never be used.
    xc['c0_0'] = 0
    yc['c0_0'] = 0

    #print("coeffs for v2v3 transform:")
    #for key in xc:
    #    print("{} {}".format(key,xc[key]))
    #sys.exit()

    xmodel = Polynomial2D(1, **xc)
    ymodel = Polynomial2D(1, **yc)

    return xmodel, ymodel

def get_refpix(row):
    #Return the reference location within the given aperture
    xref = row['XSciRef'].data[0]
    yref = row['YSciRef'].data[0]
    return Shift(-xref) & Shift(-yref)

def get_v2v3ref(row):
    #Return v2 and v3 at the reference location
    #These are arcsec in the SIAF file. Convert to degrees
    v2ref = row['V2Ref'].data[0]
    v3ref = row['V3Ref'].data[0]
    return Shift(v2ref) & Shift(v3ref)

max_stars_repo_name: mtakahiro/mirage

max_issues_repo_name: mtakahiro/mirage

max_forks_repo_name: mtakahiro/mirage

reason: import numpy,from astropy


text:
import numpy as np

from .utils import _checkTime
from .vectors import calcNae
from .vectors import calcDelta
from .vectors import calcXae
from .vectors import calcXa
from .vectors import calcNhat
from .vectors import calcR1
from .vectors import calcR2
from .projections import cartesianToGnomonic
from .coordinates import transformCoordinates

__all__ = ["TestOrbit"]

class TestOrbit:
    """
    TestOrbit: Class that calculates and stores the rotation matrices
    for a guess of heliocentric distance and velocity. To be used in
    tandem with the Cell class.

    Parameters
    ----------
    elements : `~numpy.ndarray` (6)
        Cartesian ecliptic orbital elements with postions in units of AU
        and velocities in units of AU per day.
    t0 : `~astropy.time.core.Time` (1)
        Epoch at which orbital elements are defined.
    """
    def __init__(self, elements, epoch):
        _checkTime(epoch, "epoch")

        self.elements = elements
        self.epoch = epoch

    def prepare(self, verbose=True):
        """
        Calculate rotation matrices.

        Populates the following class properties:
            n_hat : vector normal to the plane of orbit
            R1 : rotation matrix to rotate towards x-y plane
            R2 : rotation matrix to rotate towards x-axis
            M : final rotation matrix

        Parameters
        ----------
        verbose : bool, optional
            Print progress statements.
            [Default = True]

        Returns
        -------
        None
        """
        if verbose is True:
            print("Calculating vector normal to plane of orbit...")
        self.n_hat = calcNhat(self.elements[:3])

        if verbose is True:
            print("Calculating R1 rotation matrix...")
        self.R1 = calcR1(self.elements[:3], self.n_hat)
        self.x_a_xy = np.array(self.R1 @ self.elements[:3])[0]

        if verbose is True:
            print("Calculating R2 rotation matrix...")
        self.R2 = calcR2(self.x_a_xy)

        if verbose is True:
            print("Calculating final rotation matrix...")
        self.M = self.R2 @ self.R1

        if verbose is True:
            print("Done.")
            print("")
        return

    def applyToObservations(self, observations, verbose=True):
        """
        Apply the prepared rotations to the given observations. Adds the gnomonic
        plane coordinates to observations (columns: theta_x_deg, theta_y_deg)

        Parameters
        ----------
        observations : `~pandas.DataFrame`
            DataFrame of observations defined at the same epoch as this test orbit,
            to project into the test orbit's frame.
        verbose : bool, optional
            Print progress statements?
            [Default = True]

        Returns
        -------
        None
        """

        if verbose is True:
            print("Applying rotation matrices to observations...")
            print("Converting to ecliptic coordinates...")

        #velocities_present = False
        #if "vRAcosDec" in observations.columns and "vDec" in observations.columns:
        #    coords_eq_r = observations[["RA_deg", "Dec_deg"]].values
        #    coords_eq_v = observations[["vRAcosDec", "vDec"]].values
        #    coords_eq_v[:, 0] /= np.cos(np.radians(coords_eq_r[:, 1]))
        #    coords_eq = np.hstack([
        #        np.ones((len(coords_eq_r), 1)),
        #        coords_eq_r,
        #        np.zeros((len(coords_eq_r), 1)),
        #        coords_eq_v
        #    ])
        #    velocities_present = True

        #else:
        coords_eq = observations[["RA_deg", "Dec_deg"]].values
        coords_eq = np.hstack([np.ones((len(coords_eq), 1)), coords_eq])
        coords_ec = transformCoordinates(coords_eq,
            "equatorial",
            "ecliptic",
            representation_in="spherical",
            representation_out="spherical"
        )

        if verbose is True:
            print("Calculating object to observer unit vector...")
        n_ae = calcNae(coords_ec[:, 1:3])
        x_e = observations[["obs_x", "obs_y", "obs_z"]].values

        if verbose is True:
            print("Calculating object to observer distance assuming r = {} AU...".format(np.linalg.norm(self.elements[:3])))
        delta = np.zeros(len(n_ae))
        for i in range(len(delta)):
            delta[i] = calcDelta(np.linalg.norm(self.elements[:3]), x_e[i, :], n_ae[i, :])

        if verbose is True:
            print("Calculating object to observer position vector...")
        x_ae = np.zeros([len(delta), 3])
        for i, (delta_i, n_ae_i) in enumerate(zip(delta, n_ae)):
            x_ae[i] = calcXae(delta_i, n_ae_i)

        if verbose is True:
            print("Calculating heliocentric object position vector...")
        x_a = np.zeros([len(x_ae), 3])
        for i, (x_ae_i, x_e_i) in enumerate(zip(x_ae, x_e)):
            x_a[i] = calcXa(x_ae_i, x_e_i)

        if verbose is True:
            print("Applying rotation matrix M to heliocentric object position vector...")
        coords_cart_rotated = np.array(self.M @ x_a.T).T

        if verbose is True:
            print("Performing gnomonic projection...")
        gnomonic_coords = cartesianToGnomonic(coords_cart_rotated)

        observations["obj_x"] = x_a[:, 0]
        observations["obj_y"] = x_a[:, 1]
        observations["obj_z"] = x_a[:, 2]
        observations["theta_x_deg"] = np.degrees(gnomonic_coords[:, 0])
        observations["theta_y_deg"] = np.degrees(gnomonic_coords[:, 1])
        observations["test_obj_x"] = self.elements[0]
        observations["test_obj_y"] = self.elements[1]
        observations["test_obj_z"] = self.elements[2]
        observations["test_obj_vx"] = self.elements[3]
        observations["test_obj_vy"] = self.elements[4]
        observations["test_obj_vz"] = self.elements[5]

        if verbose is True:
            print("Done.")
            print("")
        return

    def applyToEphemeris(self, ephemeris, verbose=True):
        """
        Apply the prepared rotations to the given ephemerides. Adds the gnomonic
        plane coordinates to observations (columns: theta_x_deg, theta_y_deg, vtheta_x, and vtheta_y)

        Parameters
        ----------
        ephemeris : `~pandas.DataFrame`
            DataFrame of ephemeris generated by a THOR backend defined at the same epoch as this test orbit,
            to project into the test orbit's frame.
        verbose : bool, optional
            Print progress statements?
            [Default = True]

        Returns
        -------
        None
        """
        coords_cart = ephemeris[["obj_x", "obj_y", "obj_z", "obj_vx", "obj_vy", "obj_vz"]].values
        coords_cart_rotated = np.zeros_like(coords_cart)

        if verbose is True:
            print("Applying rotation matrix M to heliocentric object position vector...")
        coords_cart_rotated[:, :3] = np.array(self.M @ coords_cart[:, :3].T).T

        if verbose is True:
            print("Applying rotation matrix M to heliocentric object velocity vector...")
        # Calculate relative velocity, then rotate to projected frame
        coords_cart[:, 3:] = coords_cart[:, 3:] - self.elements[3:].reshape(1, -1)
        coords_cart_rotated[:, 3:] = np.array(self.M @ coords_cart[:, 3:].T).T

        if verbose is True:
            print("Performing gnomonic projection...")
        gnomonic_coords = cartesianToGnomonic(coords_cart_rotated)

        ephemeris["theta_x_deg"] = np.degrees(gnomonic_coords[:, 0])
        ephemeris["theta_y_deg"] = np.degrees(gnomonic_coords[:, 1])
        ephemeris["vtheta_x_deg"] = np.degrees(gnomonic_coords[:, 2])
        ephemeris["vtheta_y_deg"] = np.degrees(gnomonic_coords[:, 3])
        ephemeris["test_obj_x"] = self.elements[0]
        ephemeris["test_obj_y"] = self.elements[1]
        ephemeris["test_obj_z"] = self.elements[2]
        ephemeris["test_obj_vx"] = self.elements[3]
        ephemeris["test_obj_vy"] = self.elements[4]
        ephemeris["test_obj_vz"] = self.elements[5]

        if verbose is True:
            print("Done.")
            print("")
        return


text:
#Optimizer Function Library

#General Purpose Imports
import os, sys
import numpy as np
import scipy as sp

try:
    import pandas as pd
except:
    import pip
    pip.main(['install', 'pandas'])
    import pandas as pd

#Import Parser
from parsers import parse

from lumerical_lumapi import lumapi

#Lumopt Imports for Optimization
from lumopt.utilities.wavelengths import Wavelengths
from lumopt.geometries.polygon import FunctionDefinedPolygon
from lumopt.utilities.materials import Material
from lumopt.figures_of_merit.modematch import ModeMatch
from lumopt.optimizers.generic_optimizers import ScipyOptimizers
from lumopt.optimization import Optimization

#Insert Full Path for Datafile
datafile = os.path.join(os.path.dirname(os.path.realpath(__file__)), "y_branch_specifications.yaml")

# Default Simulation Variables
mesh_x=20e-9;
mesh_y=20e-9;
mesh_z=20e-9;
c = 3.0e8;

# CLASS FOR Y BRANCH OPTIMIZATION
class y_branch_optimization:

    #Initialize
    def __init__(self,fdtd):
        """
        Initializes FDTD session, exracts component physical parameters
        and optimization settings from YAML design intent file

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product

        Returns
        -------
        None.

        """
        #Design intent parameters from datafile
        self._fdtd = fdtd
        self._parseparse = parse()
        self._df = self._parseparse.extract_YAML(datafile)
        #Component Parsing
        self._sub_df = self._df.iloc[0]['component']
        self._sub_df = pd.DataFrame(self._sub_df)
        #Simulation parameters from datafile
        self._sim_df = self._df.iloc[0]['optimization_variables']
        self._sim_df = pd.DataFrame(self._sim_df)

    def extract_x_points(self):
        """
        Extract available space in between the input and output waveguides
        for Y Branch Creation (X axis)

        Returns
        -------
        coordinates : int array
            Range of coordinates for available space in the x axis

        """
        coordinates =  []
        for i in range(len(self._sub_df)):
            if self._sub_df.loc[i,'name'] == "input_wg" or self._sub_df.loc[i,'name'] == "output_wg_top":
                coordinates.append(self._sub_df.loc[i,'x_span'])
                coordinates.append(self._sub_df.loc[i,'x'])
        return coordinates

    #Extract available space for Y Branch Creation (Y axis)
    def extract_y_points(self):
        """
       Extract available space in between the input and output waveguides
        for Y Branch Creation (Y axis)

        Returns
        -------
        coordinates : int array
            Range of coordinates for available space in the y axis

        """
        coordinates =  []
        for i in range(len(self._sub_df)):
            if self._sub_df.loc[i,'name'] == "input_wg" or self._sub_df.loc[i,'name'] == "output_wg_top":
                coordinates.append(self._sub_df.loc[i,'y_span'])
                coordinates.append(self._sub_df.loc[i,'y'])
        return coordinates

    def calculate_x_points(self, x_coordinates):
        """
        Uses the range of coordinates for available space in the x axis and
        returns evenly spaced sequence in a specified interval

        Parameters
        ----------
        x_coordinates : int array
            Range of coordinates for available space in the x axis

        Returns
        -------
        initial_points_x: int array
            Array of evenly spaced range of values (x-axis)

        """
        self.starting_point_left = x_coordinates[1] + x_coordinates[0]/2
        self.starting_point_right = x_coordinates[3] - x_coordinates[2]/2
        self.initial_points_x = np.linspace(self.starting_point_left, self.starting_point_right, 10)
        return self.initial_points_x

    def calculate_y_points(self, y_coordinates):
        """
        Uses the range of coordinates for available space in the y axis and
        returns evenly spaced sequence in a specified interval

        Parameters
        ----------
        y_coordinates : int array
            Range of coordinates for available space in the x axis

        Returns
        -------
        initial_points_y: int array
            Array of evenly spaced range of values (y-axis)

        """
        self.starting_point_up = y_coordinates[1] + y_coordinates[0]/2
        self.starting_point_down = y_coordinates[3] + y_coordinates[2]/2
        self.initial_points_y = np.linspace(self.starting_point_up, self.starting_point_down, 10)
        return self.initial_points_y

    def opt_parameter_setup_polygon(self,initial_points_y, initial_points_x):
        """
        Defines the polygon of the Y-Splitter, also takes in the draw_y_splitter

        Parameters
        ----------
        initial_points_y : int array
            Array of evenly spaced range of values (y-axis).
        initial_points_x : int array
            Array of evenly spaced range of values (x-axis)

        Returns
        -------

        polygon: int array
            Coordinates that outline the shape of a Y branch
        """
        initial_points_x = initial_points_x
        initial_points_y = initial_points_y

        def draw_y_splitter(params):
            points_x = np.concatenate(([initial_points_x.min() - 0.01e-6], initial_points_x, [initial_points_x.max() + 0.01e-6]))
            points_y = np.concatenate(([initial_points_y.min()], params, [initial_points_y.max()]))
            n_interpolation_points = 100
            polygon_points_x = np.linspace(min(points_x), max(points_x), n_interpolation_points)
            interpolator = sp.interpolate.interp1d(points_x, points_y, kind = 'cubic')
            polygon_points_y = interpolator(polygon_points_x)
            polygon_points_up = [(x, y) for x, y in zip(polygon_points_x, polygon_points_y)]
            polygon_points_down = [(x, -y) for x, y in zip(polygon_points_x, polygon_points_y)]
            polygon_points = np.array(polygon_points_up[::-1] + polygon_points_down)
            return polygon_points
        try:
            prev_results = np.loadtxt('2D_parameters.txt')
        except:
            print("Couldn't find the file containing 2D optimization parameters. Starting with default parameters")
            prev_results = initial_points_y

        bounds = [(0, 1.2e-6)] * initial_points_y.size
        #Changing mesh orders to simulate Si
        eps_in = Material(name = 'Si: non-dispersive', mesh_order = 2)
        eps_out = Material(name = 'SiO2: non-dispersive', mesh_order = 3)
        #Defining the polygon
        polygon = FunctionDefinedPolygon(func = draw_y_splitter,
                                 initial_params = prev_results,
                                 bounds = bounds,
                                 z = 0.0,
                                 depth = 220.0e-9,
                                 eps_out = eps_out, eps_in = eps_in,
                                 edge_precision = 5,
                                 dx = 1.0e-9)
        return polygon

    def def_param_and_run(self, polygon):
        """
        Defines optimizer parameters and runs optimizer, saving optimized results in file

        Parameters
        ----------
        polygon: int array
            Coordinates that outline the shape of the intial Y branch shape

        Returns
        -------
        results int array
            Array that outlines the shape the optimized Y branch

        """

        wavelengths = Wavelengths(start = self._sim_df.loc[0,'start'],
                                  stop = self._sim_df.loc[0,'stop'],
                                  points = self._sim_df.loc[0,'points'])

        fom = ModeMatch(monitor_name = 'fom',
                mode_number = self._sim_df.loc[0,'mode_selection'],
                direction = self._sim_df.loc[0,'direction'],
                target_T_fwd = lambda wl: np.ones(wl.size),
                norm_p = 1)

        #scaling_factor = self._sim_df.loc[0,'scaling_factor']
        optimizer = ScipyOptimizers(max_iter = self._sim_df.loc[0,'max_iter'],
                            method = self._sim_df.loc[0,'method'],
                            scaling_factor = 1.0e6,
                            pgtol = self._sim_df.loc[0,'pgtol'],
                            ftol = self._sim_df.loc[0,'ftol'],
                            scale_initial_gradient_to = 0.0)

        opt = Optimization(base_script = __call__,
                   wavelengths = wavelengths,
                   fom = fom,
                   geometry = polygon,
                   optimizer = optimizer,
                   use_var_fdtd = False,
                   hide_fdtd_cad = False,
                   use_deps = True,
                   plot_history = True,
                   store_all_simulations = False)

        results = opt.run()

        np.savetxt('../3D_parameters.txt', results[1])

        return results[1]

    def main(self):
        """
        Main function

        Returns
        -------
        everything : int list
            Returns list of arrays with the following information:
                - Array that outlines the shape the optimized Y branch
                - Array of evenly spaced range of values in which the Y branch is placed in (y-axis)
                - Array of evenly spaced range of values in which the Y branch is placed in (x-axis)

        """

        self.x_coordinates = self.extract_x_points()
        self.y_coordinates = self.extract_y_points()
        self.initial_points_x = self.calculate_x_points(self.x_coordinates)
        self.initial_points_y = self.calculate_y_points(self.y_coordinates)
        #self.polygon = opt_parameter_setup_polygon(self.initial_points_y)
        self.polygon = self.opt_parameter_setup_polygon(self.initial_points_y,
                                                        self.initial_points_x)

        self.results = self.def_param_and_run(self.polygon)
        #append initial_points_x, initial_points_y to results
        everything = np.vstack((self.results,self.initial_points_x, self.initial_points_y))

        #return self.initial_points_y
        return everything

#CLASS FOR ESTABLISHING SIMULATION RECIPE IN FDTD
class recipe:

    @staticmethod
    def y_branch_init_(fdtd):
        """
        Setup FDTD session, clears any existing data

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product

        Returns
        -------
        None.

        """
        fdtd.selectall()
        fdtd.delete()

    #Extract component variables
    @staticmethod
    def create_df_component():
        """
        Extract physical component variables (waveguide info)
        from design intent YAML file

        Returns
        -------
        sub_df : panda dataframe
            Dataframe with component physical parameters
        """
            parseparse = parse()
            df = parseparse.extract_YAML(datafile)
            sub_df = df.iloc[0]['component']
            sub_df = pd.DataFrame(sub_df)
            return sub_df

    @staticmethod
    def define_material(fdtd):
        """
        Simulates/builds material that Y Branch is drawn with/built on

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product

        Returns
        -------
        Nothing signifcant to the code (sub_df is just the layering information)

        """

        parseparse = parse()
        df = parseparse.extract_YAML(datafile)
        sub_df = df.iloc[0]['layers_used']
        sub_df = pd.DataFrame(sub_df)
        _sim_df = df.iloc[0]['optimization_variables']
        _sim_df = pd.DataFrame(_sim_df)

        for i in range(len(sub_df)):

            opt_material= fdtd.addmaterial(sub_df.loc[i,'material_type'])
            opt_material = str(opt_material)
            fdtd.setmaterial(opt_material,'name',sub_df.loc[i,'name'])
            n_opt = fdtd.getindex(sub_df.loc[i,'material_name'],c/_sim_df.loc[0,'center_wavelength'])
            fdtd.setmaterial(sub_df.loc[i,'name'],'Refractive Index',n_opt)
        return sub_df

    @staticmethod
    def draw_waveguide(df,fdtd):
        """
        Draws waveguides specified in dataframe with component physical parameters

        Parameters
        ----------
        df : panda dataframe
            Dataframe with component physical parameters
        fdtd : method
            calls the contructor for the Lumerical FDTD product

        Returns
        -------
        None.

        """

        for i in range(len(df)):

            fdtd.addrect();
            fdtd.set('name',df.loc[i,'name'])
            fdtd.set('x span',df.loc[i,'x_span'])
            fdtd.set('y span',df.loc[i,'y_span'])
            fdtd.set('z span',df.loc[i,'z_span'])
            fdtd.set('y',df.loc[i,'y'])
            fdtd.set('x',df.loc[i,'x'])
            fdtd.set('z',df.loc[i,'z'])
            fdtd.set('material',df.loc[i,'material'])

            if df.loc[i,'material'] == "SiO2: non-dispersive":
                fdtd.set('override mesh order from material database',1)
                fdtd.set('mesh order',df.loc[i,'mesh_order'])
                fdtd.set('alpha',df.loc[i,'alpha'])

    #Extract simulation variables
    @staticmethod
    def extract_sim_param():
        """
        Extract simulation variables from Design Intent YAML

        Returns
        -------
        _sim_df : panda dataframe
            Dataframe with simulation parameters from Design Intent YAML

        """
        _parseparse = parse()
        _df = _parseparse.extract_YAML(datafile)
        _sim_df = _df.iloc[0]['optimization_variables']
        _sim_df = pd.DataFrame(_sim_df)
        return _sim_df

    @staticmethod
    def configure_fdtd(fdtd,component_df,sim_df):
        """
        Configures FDTD Eigensolver Dimensions and settings

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product
        component_df : panda dataframe
            Dataframe with component physical parameters
        sim_df : panda dataframe
            Dataframe with simulation parameters from Design Intent YAML

        Returns
        -------
        None.

        """
        fdtd.addfdtd();
        fdtd.set('mesh accuracy',2)
        fdtd.set('dimension',"3D")
        fdtd.set('x', 0)
        fdtd.set('x span',(component_df.loc[1,'x']+component_df.loc[1,'x_span']/2))
        fdtd.set('y',0)
        fdtd.set('y span',component_df.loc[1,'y']*10)
        fdtd.set('z',0)
        fdtd.set('z span',component_df.loc[1,'z_span']*7)
        fdtd.set('force symmetric y mesh',1)

    # Add mode source
    @staticmethod
    def add_mode_source(fdtd,component_df,sim_df):
        """
       Configures FDTD Mode Source Dimensions and Settings

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product
        component_df : panda dataframe
            Dataframe with component physical parameters
        sim_df : panda dataframe
            Dataframe with simulation parameters from Design Intent YAML

        Returns
        -------
        None.

        """
        fdtd.addmode();
        fdtd.set('direction',sim_df.loc[0,'direction'])
        fdtd.set('injection axis',sim_df.loc[0,'injection_axis'])
        fdtd.set('y',0)
        fdtd.set("y span",component_df.loc[1,'y_span']*10);
        fdtd.set('x',component_df.loc[0,'x']+(component_df.loc[0,'x_span']/4))
        fdtd.set('z span',component_df.loc[1,'z_span']*7)
        fdtd.set('center wavelength',sim_df.loc[0,'center_wavelength'])
        fdtd.set('wavelength span',0)
        fdtd.set('mode selection',sim_df.loc[0,'mode_selection'])
        #normally you would just select the fundamental TE mode, but for some reason that
        #doesn't work

    @staticmethod
    def add_port(fdtd,df):
        """
        Adds ports for S-parameter sweeps

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product
        df :panda dataframe
            Dataframe with component physical parameters.

        Returns
        -------
        None.

        """
        for i in range(len(df)):
            if df.loc[i,'name'] != "Si02":
                fdtd.addport()
                fdtd.set('name',df.loc[i,'name'])
                fdtd.set('x span', 0)
                if df.loc[i,'x']<0:
                    fdtd.set('x', (df.loc[i,'x']+(df.loc[i,'x_span']/2)))
                    #fdtd.set('x', -1e-06 )
                    fdtd.set('y', df.loc[i,'y'])
                    fdtd.set('y span', df.loc[i,'y_span'])
                else:
                    fdtd.set('x', (df.loc[i,'x']-(df.loc[i,'x_span']/2)))
                    #fdtd.set('x', 1e-06 )
                    fdtd.set('y', df.loc[i,'y'])
                    fdtd.set('y span', df.loc[i,'y_span']+0.2e-6)
                fdtd.set('z', df.loc[i,'z'])
                fdtd.set('z span',df.loc[i,'z_span'])
                if df.loc[i,'x'] <0:
                    fdtd.set('direction','Forward')
                else:
                    fdtd.set('direction','Backward')
                fdtd.set('injection axis','x-axis')
                fdtd.set('mode selection','fundamental TE mode');
                #dtd.select('Calculate Modes')
                #fdtd.select('Select Mode(s)')

    @staticmethod
    def add_mesh(fdtd,df):
        """
        Adds Mesh Region for Better Simulation Accuracy

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product
        df :panda dataframe
            Dataframe with component physical parameters.

        Returns
        -------
        None.

        """
        fdtd.addmesh();
        fdtd.set('x',0);
        fdtd.set('x span',(df.loc[1,'x']-(df.loc[1,'x_span']/2))*(9/5)*2)
        fdtd.set('y',0);
        fdtd.set('y span',(5/6)* 3e-6)
        #fdtd.set('y span',(df.loc[1,'y']+(df.loc[1,'y_span']/2))*(6/5)*2)
        fdtd.set('z',0);
        fdtd.set('z span',(df.loc[1,'z_span']/2)*30)
        fdtd.set('dx',mesh_x);
        fdtd.set('dy',mesh_y);
        fdtd.set('dz',mesh_z);

    @staticmethod
    def add_monitors(fdtd,df):
        """
        Add monitors to provide time-domain information for field components

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product
        df :panda dataframe
            Dataframe with component physical parameters.

        Returns
        -------
        None.

        """
        fdtd.addpower();
        fdtd.set('name','opt_fields');
        fdtd.set('monitor type','3D');
        fdtd.set('x',0);
        fdtd.set('x span',(df.loc[1,'x']-(df.loc[1,'x_span']/2))*(9/5)*2);
        fdtd.set('y',0);
        fdtd.set('y span',(5/6)* 3e-6);
        fdtd.set('z',0);
        fdtd.set('z span',0.4e-6);

    @staticmethod
    def add_fom(fdtd,df):
        """
        Add monitor for figure of merit

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product
        df :panda dataframe
            Dataframe with component physical parameters.

        Returns
        -------
        None.

        """
        fdtd.addpower();
        fdtd.set('name','fom');
        fdtd.set('monitor type','2D X-Normal');
        fdtd.set('x',(df.loc[1,'x']-(df.loc[1,'x_span']/2))*(9/5));
        fdtd.set('y',0);
        fdtd.set('y span',3e-6);
        fdtd.set('z',0);
        fdtd.set('z span',1.2e-6)

    # Main function
    @staticmethod
    def main(fdtd):
        """
        Main function

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product

        Returns
        -------
        Nothing signifcant to the code (layer_df is just the layering information)

        """
        #Initialize FDTD Session
        recipe.y_branch_init_(fdtd)

        component_df = recipe.create_df_component()

        layer_df = recipe.define_material(fdtd)

        sim_df = recipe.extract_sim_param()

        draw_waveguides = recipe.draw_waveguide(component_df,fdtd)

        configure_fdtd = recipe.configure_fdtd(fdtd,component_df,sim_df)

        add_mode = recipe.add_mode_source(fdtd,component_df,sim_df)

        add_mesh = recipe.add_mesh(fdtd,component_df)

        add_monitors = recipe.add_monitors(fdtd,component_df)

        add_fom = recipe.add_fom(fdtd,component_df)

        fdtd.save("y_branch_3D")

        return layer_df

    def insert_ports(fdtd):
        """
        Calls the function that adds ports for FDTD

        Parameters
        ----------
        fdtd : method
            calls the contructor for the Lumerical FDTD product

        Returns
        -------
        None.

        """

        component_df = recipe.create_df_component()
        add_ports = recipe.add_port(fdtd,component_df)

# This function is called by the optimizer and constructs the base script
def __call__(fdtd):
    """
    This function is called by the optimizer and constructs the base script
    Parameters

    ----------
    fdtd : method
            calls the contructor for the Lumerical FDTD product

    Returns
    -------
    None.

    """

    FDTD = recipe.main(fdtd)

#This class retrieves the optimization requirements
class retrieve_reqs:

    @staticmethod
    def parse_reqs():
        """
        Retrieves optimization requirements from Design Intent YAML

        Returns
        -------
        req_arr : int array
            Array that lists optimization requirements

        """
        _parseparse = parse()
        _df = _parseparse.extract_YAML(datafile)
        _opt_req_df = _df.iloc[0]['optimization_reqs']
        _opt_req_df = pd.DataFrame(_opt_req_df)

        req_arr = []
        for i in range(len(_opt_req_df)):
            req_dict.append(_opt_req_df.loc[i,'gain_through_single_port'])

        return req_arr

reason: import numpy,import scipy


text:
import numpy as np
import cv2

# 3840 2160
cap = cv2.VideoCapture("data/02.mov")
feature_params = dict( maxCorners = 100,qualityLevel = 0.3,minDistance = 7,blockSize = 7 )
lk_params = dict(winSize  = (15,15),maxLevel = 2,criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))
color = np.random.randint(0,255,(100,3))
ret, old_frame = cap.read()
old_frame = cv2.resize(old_frame, (3840//4 , 2160//4))
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
p0 = cv2.goodFeaturesToTrack(old_gray, mask = None, **feature_params)
mask = np.zeros_like(old_frame)

while(1):
    ret,frame = cap.read()
    if frame is None: break
    frame = cv2.resize(frame, (3840//4 , 2160//4))
    frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
    if(p1 is not None):
        good_new = p1[st==1]
        good_old = p0[st==1]

        for i,(new,old) in enumerate(zip(good_new,good_old)):
            a,b = new.ravel()
            c,d = old.ravel()
            mask = cv2.line(mask, (a,b),(c,d), color[i].tolist(), 2)
            frame = cv2.circle(frame,(a,b),5,color[i].tolist(),-1)
        img = cv2.add(frame,mask)
    else:
        img = frame

    cv2.imshow('frame',img)
    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1,1,2)

cv2.destroyAllWindows()
cap.release()

'''
frame_1 = cv2.imread("data/01.jpg")
frame_2 = cv2.imread("data/02.jpg")
frame_gray_1 = cv2.cvtColor(frame_1, cv2.COLOR_BGR2GRAY)
frame_gray_2 = cv2.cvtColor(frame_2, cv2.COLOR_BGR2GRAY)
feature_params = dict( maxCorners = 100,qualityLevel = 0.3,minDistance = 7,blockSize = 7 )
lk_params = dict(winSize  = (15,15),maxLevel = 2,criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))
color = np.random.randint(0,255,(100,3))

p0 = cv2.goodFeaturesToTrack(frame_gray_1, mask = None, **feature_params)#选取好的特征点，返回特征点列表
mask = np.zeros_like(frame_1)

p1, st, err = cv2.calcOpticalFlowPyrLK(frame_gray_1, frame_gray_2, p0, None, **lk_params)#计算新的一副图像中相应的特征点额位置
good_new = p1[st==1]
good_old = p0[st==1]

for i,(new,old) in enumerate(zip(good_new,good_old)):
    a,b = new.ravel() #ravel()函数用于降维并且不产生副本
    c,d = old.ravel()
    mask = cv2.line(mask, (a,b),(c,d), color[i].tolist(), 2)
    frame = cv2.circle(frame_2,(a,b),5,color[i].tolist(),-1)
img = cv2.add(frame,mask)

cv2.imshow('frame',img)
cv2.waitKey(0)
'''

max_stars_repo_name: Onicc/OpticalFlow

max_issues_repo_name: Onicc/OpticalFlow

max_forks_repo_name: Onicc/OpticalFlow


text:
import numpy as np
import matplotlib.pyplot as plt
import pyflamestk.potential as potential

# ---- test for one symbol
symbols = ['Ni']
expfunc = potential.ExponentialDensityFunction(symbols)
print('--- test exponential density function for one symbol')
print(expfunc.potential_type)
print(expfunc.parameter_names)

# --- test for two symbols
symbols = ['Ni','Al']
expfunc = potential.ExponentialDensityFunction(symbols)
print('--- test exponential density function for two symbols')
print(expfunc.potential_type)
print(expfunc.parameter_names)

# --- parameters for test
symbols = ['Ni']
param_dict = {}
param_dict['Ni_rho0'] = 1
param_dict['Ni_beta'] = 1
param_dict['Ni_r0'] = 1

r_low = 0.5
r_high = 10
N_r = 1000
rcut = 4
h = 1
r = np.linspace(r_low, r_high, N_r)

# ---- test with no cutoff function
print('--- calculate without cutoff')
symbol = 'Ni'
expfunc_nocut = potential.ExponentialDensityFunction(symbols)
vals_nocut = expfunc_nocut.evaluate(r=r,symbol='Ni',params=param_dict)

print('--- calculate with cutoff')
symbol = 'Ni'
expfunc_cut = potential.ExponentialDensityFunction(symbols)
vals_cut = expfunc_cut.evaluate(r=r,symbol='Ni',params=param_dict,
        rcut=rcut,h=h)

fig = plt.figure()
ax = plt.subplot(111)
ax.plot(r,vals_nocut,color='b',label='no_cut')
ax.plot(r,vals_cut,color='r',label='cut')
fig.savefig('exp_dens_funct.png',dpi=800)
plt.close(fig)

max_stars_repo_name: eragasa/pypospack

max_issues_repo_name: eragasa/pypospack

max_forks_repo_name: eragasa/pypospack


text:
"""
Given an image in which  two pages of a book are present, it returns two images where
each page contains one page of the book
"""
import math
import os
import matplotlib.pyplot as plt
import cv2
import numpy as np
import pandas as pd

from utils import rotate_image
from constants import get_constants

class SplitImage:
    def __init__(self, image_files):
        self._image_files = image_files
        self._origin = (0, 0)
        self._y_crop_fraction = 0.01
        self._x_crop_fraction = 0.01
        self._buffer = 40

        self._canny_apertureSize = 3
        self._vertical_text_pointcount_min = 400

        # This will be the directory relative to the location where files are present.
        self._data_directory = 'split_images'
        self._constants = get_constants()['split_image']
        self._max_angle_deviation = 5

    def _get_edges(self, img):
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        gray_cropped = self._rectangle_crop(gray)
        invert = 255 - gray_cropped
        edges = cv2.Canny(invert, 50, 150, apertureSize=self._canny_apertureSize)
        return edges

    def _rectangle_crop(self, img):
        skip_x = int(img.shape[0] / 2 * self._x_crop_fraction)
        skip_y = int(img.shape[1] / 2 * self._y_crop_fraction)
        img = img[skip_x:-skip_x, skip_y:-skip_y, ]
        return img

    def _get_bin_count(self, rotated_edges):
        non_zero = np.count_nonzero(rotated_edges, axis=0)
        df = pd.Series(non_zero, index=list(range(0, rotated_edges.shape[1])))
        df = df.groupby(np.arange(len(df)) // self._constants['bin_size'], axis=0).sum()
        return df

    def _get_y_from_bin(self, bin):
        return (bin + 0.5) * self._constants['bin_size']

    def _get_y_wrt_center(self, y_wrt_top_left_corner, img):
        return y_wrt_top_left_corner - img.shape[1] // 2

    def _get_y_wrt_left_top_corner(self, y_wrt_center, img):
        return y_wrt_center + img.shape[1] // 2

    def _get_start_end(self, df, index, vertical_text_pointcount_min=None):
        if vertical_text_pointcount_min is None:
            vertical_text_pointcount_min = self._vertical_text_pointcount_min

        start = index
        end = index
        while start > 0:
            if df.iloc[start] <= vertical_text_pointcount_min:
                start = start - 1
            else:
                break

        while end < df.shape[0] - 1:
            if df.iloc[end] <= vertical_text_pointcount_min:
                end = end + 1
            else:
                break
        return (start, end)

    def _near_center(self, index, size):
        mid = size // 2
        deviation = int(0.1 * size)
        return index in range(mid - deviation, mid + deviation)

    def _get_y_min_bin_index(self, bin_count):
        """
        Given the histogram (bin_count) of how many points lie in a y-axis range, it returns a middle bin
        """
        test_bin_count = bin_count.copy()
        while test_bin_count.shape[0] > 0:
            idx = test_bin_count.idxmin()
            if self._near_center(idx, bin_count.shape[0]):
                break
            test_bin_count = test_bin_count.drop([idx], axis=0)

        start, end = self._get_start_end(bin_count, idx, vertical_text_pointcount_min=test_bin_count.loc[idx])
        return (start + end) // 2

    def _get_score(self, rotated_edges):
        bin_count = self._get_bin_count(rotated_edges)
        y_min_bin_index = self._get_y_min_bin_index(bin_count)
        start, end = self._get_start_end(bin_count, y_min_bin_index)
        deriv = bin_count.diff().abs()
        der_left = deriv.iloc[start - 3:start + 3].max()
        der_right = deriv.iloc[end - 3:end + 3].max()
        return max(der_left, der_right)

    def _find_angle(self, edges):
        plt.figure(figsize=(20, 15))

        score_df = pd.DataFrame([], columns=['score'])
        # plot_index = 1
        for angle in range(-self._max_angle_deviation, self._max_angle_deviation):

            rotated_edges = rotate_image(edges, angle)
            score_df.loc[angle] = self._get_score(rotated_edges)

            # plt.subplot(self._max_angle_deviation, 2, plot_index)
            # plt.imshow(rotated_edges)

            # bin_count = self._get_bin_count(rotated_edges)
            # y_val = self._get_y_from_bin(self._get_y_min_bin_index(bin_count))

            # x_min = 0
            # x_max = rotated_edges.shape[0]

            # plt.plot([y_val, y_val], [x_min, x_max])
            # plt.title('Angle {} Score:{}'.format(angle * 180 / math.pi, score_df.loc[angle].values[0]))
            # plot_index += 1

        # plt.show()
        score_df = score_df.sort_values('score')
        return score_df.index[-1]

    def _splitted_filenames(self, image_file):
        tokens = image_file.split('.')
        name = '.'.join(tokens[:-1])
        fname1 = name + '_1.{extension}'.format(extension=tokens[-1])
        fname2 = name + '_2.{extension}'.format(extension=tokens[-1])

        # Create a diretory and move these filepaths inside that directory.
        new_directory = os.path.dirname(image_file) + '/{}/'.format(self._data_directory)
        if not os.path.exists(new_directory):
            os.mkdir(new_directory)

        fname1 = new_directory + os.path.basename(fname1)
        fname2 = new_directory + os.path.basename(fname2)

        return (fname1, fname2)

    def _split_file(self, image_file):
        print('Splitting file {}'.format(image_file))
        img = cv2.imread(image_file)
        img1, img2 = self._split_image(img)
        fname1, fname2 = self._splitted_filenames(image_file)
        cv2.imwrite(fname1, img1)
        cv2.imwrite(fname2, img2)
        return (fname1, fname2)

    def _split_image(self, img):
        # First work is to find the optimal angle to rotate
        edges = self._get_edges(img)
        angle = self._find_angle(edges)

        # After angle is found
        rotated_edges = rotate_image(edges, angle)

        bin_count = self._get_bin_count(rotated_edges)
        y_min_bin_index = self._get_y_min_bin_index(bin_count)
        y_val_in_rectangle = int(self._get_y_from_bin(y_min_bin_index))

        y = self._get_y_wrt_center(y_val_in_rectangle, rotated_edges)

        dst = rotate_image(img, angle)
        y_val = self._get_y_wrt_left_top_corner(y, dst)

        first_image = dst[:, :(y_val + self._buffer), :]
        second_image = dst[:, (y_val - self._buffer):, :]

        return (first_image, second_image)

    def split(self):
        output = {}
        for i, image_file in enumerate(self._image_files):
            f1, f2 = self._split_file(image_file)
            output[image_file] = (f1, f2)
            if i % 10 == 1:
                print('[ SplitImage ]: {}% complete'.format((i * 100) // len(self._image_files)))

        print('[ SplitImage ]: 100% complete')
        return output


text:
# Copyright 2019 Graphcore Ltd.
import numpy as np
import torch
from torch import nn
import popart
import onnx

from bert_model import BertConfig, Bert
from tests.torch_bert import BertConfig as TorchBertConfig, BertEmbeddings, BertLMPredictionHead, BertLayerNorm
from tests.utils import run_py, copy_weights_to_torch, run_fwd_model, check_tensors, check_model, TestFailureError

'''
Tests the embedding with a projection. This is the case for PRETRAINING.
'''
TORCH_TO_ONNX = {
    "embeddings.word_embeddings.weight": "Embedding_Dict",
    "embeddings.LayerNorm.weight": "Gamma",
    "embeddings.LayerNorm.bias": "Beta",
    "cls.transform.dense.weight": "CLS/LMPredictionW",
    "cls.transform.dense.bias": "CLS/LMPredictionB",
    "cls.transform.LayerNorm.weight": "CLS/Gamma",
    "cls.transform.LayerNorm.bias": "CLS/Beta"
}
TRANSPOSE_WEIGHTS = {
    "embeddings.word_embeddings.weight": np.transpose,
    "cls.transform.dense.weight": np.transpose
}

class BertEmbeddingsNoPosition(nn.Module):
    def __init__(self, config):
        super(BertEmbeddingsNoPosition, self).__init__()
        self.word_embeddings = nn.Embedding(
            config.vocab_size, config.hidden_size, padding_idx=0)
        self.LayerNorm = BertLayerNorm(
            config.hidden_size, eps=config.layer_norm_eps)

    def forward(self, input_ids):
        seq_length = input_ids.size(1)
        words_embeddings = self.word_embeddings(input_ids)
        words_embeddings = words_embeddings.detach()
        embeddings = self.LayerNorm(words_embeddings)
        return embeddings

class EmbeddingProjectionModel(nn.Module):
    def __init__(self, config):
        super().__init__()
        self.embeddings = BertEmbeddingsNoPosition(config)
        self.cls = BertLMPredictionHead(config)
        self.tie_weights()

    def tie_weights(self):
        self.cls.decoder.weight = self.embeddings.word_embeddings.weight
        self.embeddings.word_embeddings.weight.data = \
            self.cls.decoder.weight.data.detach()

    def forward(self, indices):
        sequence_output = self.embeddings(indices)
        projected = self.cls(sequence_output)
        return projected

def test_embedding_projection_fwd(custom_ops):
    #  ------------------- PopART --------------------
    builder = popart.Builder(opsets={
        "ai.onnx": 9,
        "ai.onnx.ml": 1,
        "ai.graphcore": 1
    })
    config = BertConfig(vocab_length=9728,
                        projection_serialization_steps=4,
                        batch_size=1,
                        hidden_size=768,
                        sequence_length=128,
                        activation_type='relu',
                        popart_dtype="FLOAT",
                        no_dropout=True,
                        inference=True)
    popart_model = Bert(config, builder=builder)

    sequence_info = popart.TensorInfo(
        "UINT32", [config.batch_size * config.sequence_length])
    indices = builder.addInputTensor(sequence_info)
    data = {
        indices:
        np.random.randint(0, config.vocab_length,
                          (config.batch_size * config.sequence_length)).astype(
                              np.uint32)
    }

    x = popart_model.gather(
        indices, config.vocab_length, "Embedding_Dict")
    x = popart_model.norm(x)
    x = popart_model.dropout(x)
    with popart_model.builder.nameScope("CLS"):
        x = popart_model.lm_prediction_head(x)
    output = popart_model.projection(x)

    proto = builder.getModelProto()

    outputs, post_proto = run_py(proto, data, output)

    # ----------------- PopART -> PyTorch ----------------
    proto = onnx.load_model_from_string(proto)

    inputs = [data[indices].reshape(
        config.batch_size, config.sequence_length).astype(np.int32)]

    #  ------------------- PyTorch -------------------------
    torch_model = EmbeddingProjectionModel(
        TorchBertConfig(config.vocab_length,
                        config.hidden_size,
                        max_position_embeddings=config.max_positional_length,
                        layer_norm_eps=config.layer_norm_eps))

    torch_model.eval()

    copy_weights_to_torch(torch_model, proto, TORCH_TO_ONNX,
                          TRANSPOSE_WEIGHTS)
    torch_model.tie_weights()

    torch_outputs = run_fwd_model(inputs, torch_model)

    check_tensors(torch_outputs, outputs)

def test_embedding_projection_bwd(custom_ops):
    l1_lambda = 0.1

    #  ------------------- PopART --------------------
    builder = popart.Builder(opsets={
        "ai.onnx": 9,
        "ai.onnx.ml": 1,
        "ai.graphcore": 1
    })
    config = BertConfig(vocab_length=9728,
                        projection_serialization_steps=4,
                        batch_size=1,
                        hidden_size=768,
                        sequence_length=128,
                        activation_type='relu',
                        popart_dtype="FLOAT",
                        no_dropout=True,
                        # Currently updating embedding dict with projection is only
                        # available with momentum. And PopART != Pytorch momentum
                        # due to a bootstrapping step on iter 0.
                        update_embedding_dict=False)
    popart_model = Bert(config, builder=builder)

    sequence_info = popart.TensorInfo(
        "UINT32", [config.batch_size * config.sequence_length])
    indices = builder.addInputTensor(sequence_info)
    data = {
        indices:
        np.random.randint(0, config.vocab_length,
                          (config.batch_size * config.sequence_length)).astype(
                              np.uint32)
    }

    x = popart_model.gather(
        indices, config.vocab_length, "Embedding_Dict")
    x = popart_model.norm(x)
    x = popart_model.dropout(x)
    with popart_model.device_scope(nameScope="CLS"):
        x = popart_model.lm_prediction_head(x)
    output = popart_model.projection(x)
    l1 = builder.aiGraphcore.l1loss(
        [output], l1_lambda, debugPrefix="l1LossVal", reduction=popart.ReductionType.Sum)

    proto = builder.getModelProto()

    optimizer = popart.ConstSGD(0.01)

    outputs, post_proto = run_py(proto,
                                 data, output,
                                 loss=l1,
                                 optimizer=optimizer)

    # ----------------- PopART -> PyTorch ----------------
    proto = onnx.load_model_from_string(proto)

    inputs = [data[indices].reshape(
        config.batch_size, config.sequence_length).astype(np.int32)]

    #  ------------------- PyTorch -------------------------

    torch_model = EmbeddingProjectionModel(
        TorchBertConfig(config.vocab_length,
                        config.hidden_size,
                        max_position_embeddings=config.max_positional_length,
                        layer_norm_eps=config.layer_norm_eps))
    # Turn off dropout
    torch_model.eval()

    copy_weights_to_torch(torch_model,
                          proto,
                          TORCH_TO_ONNX,
                          transform=TRANSPOSE_WEIGHTS)

    optim = torch.optim.SGD(torch_model.parameters(),
                            0.01,
                            weight_decay=0.0,
                            momentum=0.0)

    torch_output = torch_model(*[torch.from_numpy(t).long() for t in inputs])
    torch_loss = l1_lambda * torch.norm(torch_output, 1)
    torch_loss.backward()
    optim.step()

    check_tensors([torch_output.detach().numpy()], outputs, margin=1e-5)

    check_model(torch_model,
                post_proto,
                TORCH_TO_ONNX,
                transform=TRANSPOSE_WEIGHTS)


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 16 20:23:01 2020

@author: Salomé
"""
# Standard library imports
import numpy as np

# Mosqito functions import
from mosqito.sq_metrics.tonality.tone_to_noise_ecma._critical_band import _critical_band

def _find_highest_tone(freqs, spec_db, index, nb_tones, ind):
    """

    Method to find the two highest tones in a given spectrum from a given index
    according to their critical band

    Parameters
    ----------
    freqs : numpy.array
        frequency axis
    spec_db : numpy.array
        signal spectrum in dB
    index : numpy.array
        list of candidate tones index
    index : numpy.array
        list of candidate tones index
    nb_tones : integer
        number of candidate tones non examinated

    Returns
    -------
    ind_p : integer
        index of the highest tone in the critical band
    ind_p : integer
        index of the second highest tone in the critical band
    index : numpy.array
        list of candidate tones index updated
    nb_tones : integer
        number of candidate tones non examinated updated
    """

    f = freqs[ind]
    # critical band centered on f
    f1, f2 = _critical_band(f)
    low_limit_idx = np.argmin(np.abs(freqs - f1))
    high_limit_idx = np.argmin(np.abs(freqs - f2))

    # Other tones in the critical band centered on f tones
    multiple_idx = index[index > low_limit_idx]
    multiple_idx = multiple_idx[multiple_idx < high_limit_idx]

    if len(multiple_idx) > 1:
        sort_spec = np.argsort(-1 * spec_db[multiple_idx])

        # highest tones in the critical band
        ind_p = multiple_idx[sort_spec[0]]
        ind_s = multiple_idx[sort_spec[1]]

        # suppression of the lower values
        for s in sort_spec[2:]:
            sup = np.where(index == multiple_idx[s])[0]
            index = np.delete(index, sup)
            nb_tones -= 1

        if ind_p != ind:
            # screening to find the highest value in the critical band centered on fp
            ind_p, ind_s, index, nb_tones = _find_highest_tone(
                freqs, spec_db, index, nb_tones, ind_p
            )

    else:
        ind_p = ind
        ind_s = None

    return ind_p, ind_s, index, nb_tones

max_stars_repo_name: MitchellAcoustics/MoSQITo

max_issues_repo_name: MitchellAcoustics/MoSQITo

max_forks_repo_name: MitchellAcoustics/MoSQITo


text:
#! -*- coding: utf-8 -*-
# 代码合集

import six
import logging
import numpy as np
import re

if not six.PY2:
    basestring = str

def is_string(s):
    """判断是否是字符串
    """
    return isinstance(s, basestring)

def strQ2B(ustring):
    """全角符号转对应的半角符号
    """
    rstring = ''
    for uchar in ustring:
        inside_code = ord(uchar)
        # 全角空格直接转换
        if inside_code == 12288:
            inside_code = 32
        # 全角字符（除空格）根据关系转化
        elif (inside_code >= 65281 and inside_code <= 65374):
            inside_code -= 65248
        rstring += unichr(inside_code)
    return rstring

def string_matching(s, keywords):
    """判断s是否至少包含keywords中的至少一个字符串
    """
    for k in keywords:
        if re.search(k, s):
            return True
    return False

class Progress:
    """显示进度，自己简单封装，比tqdm更可控一些
    iterable: 可迭代的对象；
    period: 显示进度的周期；
    steps: iterable可迭代的总步数，相当于len(iterable)
    """
    def __init__(self, iterable, period=1, steps=None, desc=None):
        self.iterable = iterable
        self.period = period
        if hasattr(iterable, '__len__'):
            self.steps = len(iterable)
        else:
            self.steps = steps
        self.desc = desc
        if self.steps:
            self._format_ = u'%s/%s passed' % ('%s', self.steps)
        else:
            self._format_ = u'%s passed'
        if self.desc:
            self._format_ = self.desc + ' - ' + self._format_
        self.logger = logging.getLogger()

    def __iter__(self):
        for i, j in enumerate(self.iterable):
            if (i + 1) % self.period == 0:
                self.logger.info(self._format_ % (i + 1))
            yield j

def parallel_apply(func,
                   iterable,
                   workers,
                   max_queue_size,
                   callback=None,
                   dummy=False):
    """多进程或多线程地将func应用到iterable的每个元素中。
    注意这个apply是异步且无序的，也就是说依次输入a,b,c，但是
    输出可能是func(c), func(a), func(b)。
    参数：
        dummy: False是多进程/线性，True则是多线程/线性；
        callback: 处理单个输出的回调函数；
    """
    if dummy:
        from multiprocessing.dummy import Pool, Queue
    else:
        from multiprocessing import Pool, Queue

    in_queue, out_queue = Queue(max_queue_size), Queue()

    def worker_step(in_queue, out_queue):
        # 单步函数包装成循环执行
        while True:
            d = in_queue.get()
            r = func(d)
            out_queue.put(r)

    # 启动多进程/线程
    pool = Pool(workers, worker_step, (in_queue, out_queue))

    if callback is None:
        results = []

    # 后处理函数
    def process_out_queue():
        out_count = 0
        for _ in range(out_queue.qsize()):
            d = out_queue.get()
            out_count += 1
            if callback is None:
                results.append(d)
            else:
                callback(d)
        return out_count

    # 存入数据，取出结果
    in_count, out_count = 0, 0
    for d in iterable:
        in_count += 1
        while True:
            try:
                in_queue.put(d, block=False)
                break
            except six.moves.queue.Full:
                out_count += process_out_queue()
        if in_count % max_queue_size == 0:
            out_count += process_out_queue()

    while out_count != in_count:
        out_count += process_out_queue()

    pool.terminate()

    if callback is None:
        return results

def get_all_attributes(something):
    """获取类下的所有属性和方法
    """
    return {
        name: getattr(something, name)
        for name in dir(something) if name[0] != '_'
    }

def sequence_padding(inputs, length=None, padding=0):
    """Numpy函数，将序列padding到同一长度
    """
    if length is None:
        length = max([len(x) for x in inputs])

    outputs = np.array([
        np.concatenate([x, [padding] * (length - len(x))])
        if len(x) < length else x[:length] for x in inputs
    ])
    return outputs

def is_one_of(x, ys):
    """判断x是否在ys之中
    等价于x in ys，但有些情况下x in ys会报错
    """
    for y in ys:
        if x is y:
            return True
    return False