text:
import numpy as np
import arviz as az
import seaborn as sns
import pandas as pd
import pickle as pkl
import importlib
import anndata as ad
import ast

import matplotlib.pyplot as plt
from scdcdm.util import result_classes as res
from scdcdm.util import multi_parameter_sampling as mult
from scdcdm.util import multi_parameter_analysis_functions as ana
from scdcdm.util import data_generation as gen

pd.set_option('display.max_columns', 500)
pd.set_option('display.max_rows', 500)
sns.set_style("ticks")

#%%
importlib.reload(ana)
# Get data

path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\overall_benchmark"

#%%
# Find relations between numerical values and composition/increase vectors

b = []
for y1_0 in [20, 30, 50, 75, 115, 180, 280, 430, 667, 1000]:
    b_i = np.round(gen.counts_from_first(y1_0, 5000, 5), 3)
    b.append(np.round(np.log(b_i / 5000), 2))
print(b)

b_counts = dict(zip([b_i[0] for b_i in b], [20, 30, 50, 75, 115, 180, 280, 430, 667, 1000]))

b2 = []
for y1_0 in [20, 30, 50, 75, 115, 180, 280, 430, 667, 1000]:
    b_i = np.round(gen.counts_from_first(y1_0, 5000, 5), 3)
    b2.append(b_i)

b_w_dict = {}
i = 0
w_all = []
for b_i in b2:
    b_t = np.round(np.log(b_i / 5000), 3)
    print(b_t)
    w_d = {}
    for change in [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 400, 600, 800, 1000]:
        _, w = gen.b_w_from_abs_change(b_i, change, 5000)
        w_0 = np.round(w[0], 3)
        w_d[w_0] = change
    b_w_dict[b_t[0]] = w_d
    i += 1
print(b_w_dict)

#%%
# all_study_params: One line per data point
# all_study_params_agg: Aggregate identical simulation parameters
p = []
params = []
for b_i in b:
    print(b_i[0])
    _, all_study_params, all_study_params_agg = ana.multi_run_study_analysis_prepare(path,
                                                                      file_identifier=str(b_i[0]))

    all_study_params_agg = ana.get_scores(all_study_params_agg)
    all_study_params = ana.get_scores(all_study_params)
    params.append(all_study_params_agg)
    p.append(all_study_params)

#%%
# Get some new metrics that make plotting more convenient

all_study_params_agg_2 = pd.concat(params)
all_study_params_agg_2["n_controls"] = [ast.literal_eval(x)[0] for x in all_study_params_agg_2["n_samples"].tolist()]
all_study_params_agg_2["n_cases"] = [ast.literal_eval(x)[1] for x in all_study_params_agg_2["n_samples"].tolist()]
all_study_params_agg_2["n_total"] = all_study_params_agg_2["n_total"].astype("float")
all_study_params_agg_2["w"] = [ast.literal_eval(x)[0][0] for x in all_study_params_agg_2["w_true"]]
all_study_params_agg_2["b_0"] = [ast.literal_eval(x)[0] for x in all_study_params_agg_2["b_true"]]
all_study_params_agg_2["b_count"] = [b_counts[np.round(x, 2)] for x in all_study_params_agg_2["b_0"]]

bs = all_study_params_agg_2["b_0"].tolist()
ws = all_study_params_agg_2["w"].tolist()
increases = []
for i in range(len(bs)):
    increases.append(b_w_dict[bs[i]][ws[i]])
all_study_params_agg_2["num_increase"] = increases
all_study_params_agg_2["log_fold_increase"] = np.log2((all_study_params_agg_2["num_increase"] +
                                                       all_study_params_agg_2["b_count"]) /
                                                      all_study_params_agg_2["b_count"])

print(all_study_params_agg_2)

#%%
# Get some new metrics that make plotting more convenient

all_study_params_2 = pd.concat(p)
all_study_params_2["n_controls"] = [ast.literal_eval(x)[0] for x in all_study_params_2["n_samples"].tolist()]
all_study_params_2["n_cases"] = [ast.literal_eval(x)[1] for x in all_study_params_2["n_samples"].tolist()]
all_study_params_2["n_total"] = all_study_params_2["n_total"].astype("float")
all_study_params_2["w"] = [ast.literal_eval(x)[0][0] for x in all_study_params_2["w_true"]]
all_study_params_2["b_0"] = [ast.literal_eval(x)[0] for x in all_study_params_2["b_true"]]
all_study_params_2["b_count"] = [b_counts[np.round(x, 2)] for x in all_study_params_2["b_0"]]

bs = all_study_params_2["b_0"].tolist()
ws = all_study_params_2["w"].tolist()
increases = []
for i in range(len(bs)):
    increases.append(b_w_dict[bs[i]][ws[i]])
all_study_params_2["num_increase"] = increases
all_study_params_2["log_fold_increase"] = np.log2((all_study_params_2["num_increase"] +
                                                       all_study_params_2["b_count"]) /
                                                      all_study_params_2["b_count"])

print(all_study_params_2)

#%%
# Save that stuff
result_path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\benchmark_results"

with open(result_path + "\\normal_results_aggregated.pkl", "wb") as f:
    pkl.dump(all_study_params_agg_2, file=f, protocol=4)
with open(result_path + "\\normal_results.pkl", "wb") as f:
    pkl.dump(all_study_params_2, file=f, protocol=4)

#%%
# Exploration...
sns.heatmap(data=all_study_params_agg_2[["b_count", "num_increase", "mcc"]].pivot_table("mcc", "b_count", "num_increase"), vmin=-1, vmax=1)
plt.show()

#%%

sns.heatmap(data=all_study_params_agg_2[["n_cases", "n_controls", "mcc"]].pivot_table("mcc", "n_controls", "n_cases"), vmin=-1, vmax=1)
plt.show()

#%%

print(all_study_params_agg_2[all_study_params_agg_2["n_samples"]=="[10, 10]"])

#%%
# One heatmap for every combination of b, w
def draw_heatmap(*args, **kwargs):
    data = kwargs.pop('data')
    d = data.pivot(index=args[1], columns=args[0], values=args[2])
    sns.heatmap(d, **kwargs, vmin=-1, vmax=1)

fg = sns.FacetGrid(all_study_params_agg_2, col='num_increase', row="b_count")
fg.map_dataframe(draw_heatmap, 'n_controls', 'n_cases', 'mcc', cbar=False)
plt.show()

#%%

sns.lineplot(data=all_study_params_agg_2, x="b_count", y="mcc")
plt.show()

#%%
sns.lineplot(data=all_study_params_agg_2, x="num_increase", y="mcc", hue="b_count")
plt.show()

#%%
sns.lineplot(data=all_study_params_agg_2, x="log_fold_increase", y="mcc", hue="b_count")
plt.show()

#%%
# Inclusion probability - No results here!

def recalculate_inclusion_probability(results, threshold):

    for r in results:
        for res in r.mcmc_results.values():
            res[1]["final_parameter"] = np.where(res[1]["inclusion_prob"] > threshold,
                                                 res[1]["mean_nonzero"],
                                                 0)
        r.get_discovery_rates()

    # Generate all_study_params
    all_study_params = pd.concat([r.parameters for r in results])
    simulation_parameters = ["cases", "K", "n_total", "n_samples", "b_true", "w_true", "num_results"]
    all_study_params[simulation_parameters] = all_study_params[simulation_parameters].astype(str)

    # Aggregate over identical parameter sets
    all_study_params_agg = all_study_params.groupby(simulation_parameters).sum()
    all_study_params_agg = ana.get_scores(all_study_params_agg.reset_index())

    return all_study_params_agg

def multiple_threshold_compare(results, thresholds, results_df):

    for t in thresholds:
        print("threshold: " + str(t))
        t_results = recalculate_inclusion_probability(results, t)

        t_results["threshold"] = t

        results_df = results_df.append(t_results)

    return results_df

import os

def multiple_thresholds_with_load(path, thresholds, file_identifier="result_"):
    files = os.listdir(path)

    results = []

    print("Calculating discovery rates...")
    i = 0

    # For all files:
    for f in files:
        i += 1

        print("Preparing: ", i / len(files))
        if file_identifier in f:
            # Load file
            r = ana.renamed_load(open(path + "/" + f, "rb"))

            results.append(r)

    results_df = pd.DataFrame()

    results_df = multiple_threshold_compare(results, thresholds, results_df)

    return results_df

id = "b_[[-2.01 -1.53 -1.53 -1.53 -1.53]]_w_[[[0.034, 0.0, 0.0, 0.0, 0.0]]]"
results_df = multiple_thresholds_with_load(path, [0.5, 0.7], file_identifier=id)

#%%
# Calculate metrics for different threshold values

thresholds = np.round(np.arange(0.05, 1, 0.05), 2)

total_df = pd.DataFrame()
for b_i in b:
    print(b_i[0])
    results_df = multiple_thresholds_with_load(path, thresholds, file_identifier=str(b_i[0]))

    total_df = total_df.append(results_df)

#%%

sns.lineplot(data=total_df, x="threshold", y="mcc", hue="")
plt.show()

#%%
# extra metrics

total_df_2 = total_df
total_df_2["n_controls"] = [ast.literal_eval(x)[0] for x in total_df_2["n_samples"].tolist()]
total_df_2["n_cases"] = [ast.literal_eval(x)[1] for x in total_df_2["n_samples"].tolist()]
total_df_2["n_total"] = total_df_2["n_total"].astype("float")
total_df_2["w"] = [ast.literal_eval(x)[0][0] for x in total_df_2["w_true"]]
total_df_2["b_0"] = [ast.literal_eval(x)[0] for x in total_df_2["b_true"]]
total_df_2["b_count"] = [b_counts[np.round(x, 2)] for x in total_df_2["b_0"]]

bs = total_df_2["b_0"].tolist()
ws = total_df_2["w"].tolist()
increases = []
for i in range(len(bs)):
    increases.append(b_w_dict[bs[i]][ws[i]])
total_df_2["num_increase"] = increases
total_df_2["log_fold_increase"] = np.log2((total_df_2["num_increase"] +
                                           total_df_2["b_count"]) /
                                          total_df_2["b_count"])

print(total_df_2)

#%%
# Saving
result_path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\benchmark_results"

with open(result_path + "\\thresholds.pkl", "wb") as f:
    pkl.dump(total_df_2, file=f, protocol=4)

#%%
# more plots
fig, ax = plt.subplots(1, 2, figsize=(13, 5))
sns.lineplot(data=total_df_2, x="threshold", y="tpr", hue="num_increase", ax=ax[0])
plt.axvline(x=0.56)
sns.lineplot(data=total_df_2, x="threshold", y="tnr", hue="num_increase", ax=ax[1])
plt.axvline(x=0.56)

plt.show()

sns.lineplot(data=total_df_2.loc[total_df_2["num_increase"]>50], x="threshold", y="mcc")
plt.axvline(x=0.56)
plt.show()

#%%

# Negative effect testing
# Same setup as before
importlib.reload(ana)

path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\negative_benchmark"

_, all_study_params_neg, all_study_params_neg_agg = ana.multi_run_study_analysis_prepare(path)

all_study_params_neg_agg = ana.get_scores(all_study_params_neg_agg)
all_study_params_neg = ana.get_scores(all_study_params_neg)

b = []
for y1_0 in [115, 280, 1000]:
    b_i = np.round(gen.counts_from_first(y1_0, 5000, 5), 3)
    b.append(np.round(np.log(b_i / 5000), 2))
print(b)

b_counts = dict(zip([b_i[0] for b_i in b], [115, 280, 1000]))

b2 = []
for y1_0 in [115, 280, 1000]:
    b_i = np.round(gen.counts_from_first(y1_0, 5000, 5), 3)
    b2.append(b_i)

b_w_dict = {}
i = 0
w_all = []
for b_i in b2:
    b_t = np.round(np.log(b_i / 5000), 3)
    print(b_t)
    w_d = {}
    for change in [-10, -50, -100]:
        _, w = gen.b_w_from_abs_change(b_i, change, 5000)
        w_0 = np.round(w[0], 3)
        w_d[w_0] = change
    b_w_dict[b_t[0]] = w_d
    i += 1
print(b_w_dict)

all_study_params_neg = all_study_params_neg
all_study_params_neg["n_controls"] = [ast.literal_eval(x)[0] for x in all_study_params_neg["n_samples"].tolist()]
all_study_params_neg["n_cases"] = [ast.literal_eval(x)[1] for x in all_study_params_neg["n_samples"].tolist()]
all_study_params_neg["n_total"] = all_study_params_neg["n_total"].astype("float")
all_study_params_neg["w"] = [ast.literal_eval(x)[0][0] for x in all_study_params_neg["w_true"]]
all_study_params_neg["b_0"] = [ast.literal_eval(x)[0] for x in all_study_params_neg["b_true"]]
all_study_params_neg["b_count"] = [b_counts[np.round(x, 2)] for x in all_study_params_neg["b_0"]]

bs = all_study_params_neg["b_0"].tolist()
ws = all_study_params_neg["w"].tolist()
increases = []
for i in range(len(bs)):
    increases.append(b_w_dict[bs[i]][ws[i]])
all_study_params_neg["num_increase"] = increases
all_study_params_neg["log_fold_increase"] = np.log2((all_study_params_neg["num_increase"] +
                                                       all_study_params_neg["b_count"]) /
                                                      all_study_params_neg["b_count"])

print(all_study_params_neg)

all_study_params_neg_agg["n_controls"] = [ast.literal_eval(x)[0] for x in all_study_params_neg_agg["n_samples"].tolist()]
all_study_params_neg_agg["n_cases"] = [ast.literal_eval(x)[1] for x in all_study_params_neg_agg["n_samples"].tolist()]
all_study_params_neg_agg["n_total"] = all_study_params_neg_agg["n_total"].astype("float")
all_study_params_neg_agg["w"] = [ast.literal_eval(x)[0][0] for x in all_study_params_neg_agg["w_true"]]
all_study_params_neg_agg["b_0"] = [ast.literal_eval(x)[0] for x in all_study_params_neg_agg["b_true"]]
all_study_params_neg_agg["b_count"] = [b_counts[np.round(x, 2)] for x in all_study_params_neg_agg["b_0"]]

bs = all_study_params_neg_agg["b_0"].tolist()
ws = all_study_params_neg_agg["w"].tolist()
increases = []
for i in range(len(bs)):
    increases.append(b_w_dict[bs[i]][ws[i]])
all_study_params_neg_agg["num_increase"] = increases
all_study_params_neg_agg["log_fold_increase"] = np.log2((all_study_params_neg_agg["num_increase"] +
                                                       all_study_params_neg_agg["b_count"]) /
                                                      all_study_params_neg_agg["b_count"])

print(all_study_params_neg_agg)

#%%
# Save that stuff
result_path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\benchmark_results"

with open(result_path + "\\results_negative.pkl", "wb") as f:
    pkl.dump(all_study_params_neg, file=f, protocol=4)
with open(result_path + "\\results_negative_aggregated.pkl", "wb") as f:
    pkl.dump(all_study_params_neg_agg, file=f, protocol=4)
#%%
# Plots
result_path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\benchmark_results"

def draw_heatmap(*args, **kwargs):
    data = kwargs.pop('data')
    d = data.pivot(index=args[1], columns=args[0], values=args[2])
    sns.heatmap(d, **kwargs, vmin=-1, vmax=1)

fg = sns.FacetGrid(all_study_params_neg_agg, col='num_increase', row="b_count")
fg.map_dataframe(draw_heatmap, 'n_controls', 'n_cases', 'mcc', cbar=False)
plt.savefig(result_path + "\\negative_heatmaps.png")
plt.show()

#%%

# For comparison: Positive heatmaps
result_path = "C:\\Users\\Johannes\\Documents\\Uni\\Master's_Thesis\\compositionalDiff-johannes_tests_2\\data\\benchmark_results"

with open(result_path + "\\results_aggregated.pkl", "rb") as f:
    all_study_params_agg_pos = pkl.load(file=f)

#%%

fg = sns.FacetGrid(all_study_params_agg_pos.loc[(all_study_params_agg_pos["b_count"].isin([115, 280, 1000])) &
                                              (all_study_params_agg_pos["num_increase"].isin([10, 50, 100]))],
                   col='num_increase', row="b_count")
fg.map_dataframe(draw_heatmap, 'n_controls', 'n_cases', 'mcc', cbar=False)
plt.savefig(result_path + "\\pos_heatmaps_for_negative.png")
plt.show()

# -> For same absolute cell count change, negative effects are detected even a little better!

#%%

def draw_heatmap(*args, **kwargs):
    data = kwargs.pop('data')
    d = data.pivot(index=args[1], columns=args[0], values=args[2])
    sns.heatmap(d, **kwargs, vmin=-1, vmax=1)

fg = sns.FacetGrid(all_study_params_agg_pos.loc[(all_study_params_agg_pos["b_count"].isin([50, 180, 430])) &
                                              (all_study_params_agg_pos["num_increase"].isin([40, 80, 200]))],
                   col='num_increase', row="b_count")
fg.map_dataframe(draw_heatmap, 'n_controls', 'n_cases', 'mcc', cbar=False)
plt.savefig(result_path + "\\negative_heatmaps.png")
plt.show()

max_stars_repo_name: neurogenomics/scCODA

max_issues_repo_name: neurogenomics/scCODA

max_forks_repo_name: neurogenomics/scCODA


text:
"""
Synchronization module
======================

This module provides functions to process the synchronization data
acquired with Thor Sync during imaging.
"""
import warnings
import json

import numpy as np
import h5py
import scipy.signal

import utils2p.main as main

class SynchronizationError(Exception):
    """The input data is not consistent with synchronization assumption."""

def get_lines_from_h5_file(file_path, line_names):
    warnings.warn(
        "get_lines_from_h5_file is deprecated use get_lines_from_sync_file instead",
        DeprecationWarning)
    return get_lines_from_sync_file(file_path, line_names)

def get_lines_from_sync_file(file_path, line_names):
    """
    This function returns the values of the requested lines save in
    an h5 generated by ThorSync.

    Parameters
    ----------
    file_path : string
        Path to h5 file.
    line_names : list of strings
        List of the ThorSync line names to be returned.

    Returns
    -------
    lines : tuple
        Line arrays in the same order as given in line_names.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_z_stack")
    >>> line_names = ["Frame Counter", "Capture On"]
    >>> frame_counter, capture_on = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> type(frame_counter)
    <class 'numpy.ndarray'>
    >>> frame_counter.shape
    (54000,)
    >>> type(capture_on)
    <class 'numpy.ndarray'>
    >>> capture_on.shape
    (54000,)
    """
    lines = []

    with h5py.File(file_path, "r") as f:
        for name in line_names:
            lines_with_this_name = []
            for line_type in ("DI", "CI", "AI"):
                try:
                    lines_with_this_name.append(
                        f[line_type][name][:].squeeze())
                except KeyError:
                    pass
            if len(lines_with_this_name) == 1:
                lines.append(lines_with_this_name[0])
            elif len(lines_with_this_name) == 0:
                DI_keys = list(f["DI"].keys())
                CI_keys = list(f["CI"].keys())
                AI_keys = list(f["AI"].keys())
                raise KeyError(
                    f"No line named '{name}' exists. The digital lines are " +
                    f"{DI_keys}, the continuous lines are {CI_keys}, " +
                    f"and  the analogue inputs are {AI_keys}.")
            else:
                DI_keys = list(f["DI"].keys())
                CI_keys = list(f["CI"].keys())
                AI_keys = list(f["AI"].keys())
                raise KeyError(
                    f"Multiple lines named '{name}' exist. " +
                    f"The digital lines are {DI_keys}, the continuous lines " +
                    + f"are {CI_keys}, and  the analogue inputs are {AI_keys}."
                )
    return tuple(lines)

def get_times(length, freq):
    """
    This function returns the time point of each tick
    for a given sequence length and tick frequency.

    Parameters
    ----------
    length : int
        Length of sequence.
    freq : float
        Frequency in Hz.

    Returns
    -------
    times : array
        Times in seconds.

    Examples
    --------
    >>> import utils2p.synchronization
    >>> utils2p.synchronization.get_times(5, 20)
    array([0.  , 0.05, 0.1 , 0.15, 0.2 ])
    """
    times = np.arange(0, length / freq, 1 / freq)
    return times

def edges(line, size=0, correct_possible_split_edges=True):
    """
    Returns the indices of edges in a line. An
    edge is change in value of the line. A size
    argument can be specified to filter for changes
    of specific magnitude. By default only rising
    edges (increases in value) are returned.

    Parameters
    ----------
    line : numpy array
        Line signal from h5 file.
    size : float or tuple
        Size of the rising edge. If float it is used as minimum.
        Tuples specify a range. To get falling edges use negative values.
        Only one boundary can be applied using np.inf as one of the values.
        All boundaries are excluding the specified value.
    correct_possible_split_edges : boolean
        The rise or fall of an edge can in some cases be spread over
        several ticks. If `True` these "blurry" edges are sharpened
        with :func:`utils2p.synchronization.correct_split_edges`.
        Default is True.

    Returns
    -------
    indices : list
        Indices of the rising edges.

    Examples
    --------
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> binary_line = np.array([0, 1, 1, 0, 1, 1])
    >>> utils2p.synchronization.edges(binary_line)
    (array([1, 4]),)
    >>> utils2p.synchronization.edges(binary_line, size=2)
    (array([], dtype=int64),)
    >>> utils2p.synchronization.edges(binary_line, size=(-np.inf, np.inf))
    (array([1, 3, 4]),)
    >>> continuous_line = np.array([0, 0, 3, 3, 3, 5, 5, 8, 8, 10, 10, 10])
    >>> utils2p.synchronization.edges(continuous_line)
    (array([2, 5, 7, 9]),)
    >>> utils2p.synchronization.edges(continuous_line, size=2)
    (array([2, 7]),)
    >>> utils2p.synchronization.edges(continuous_line, size=(-np.inf, 3))
    (array([5, 9]),)
    """
    if correct_possible_split_edges:
        line = correct_split_edges(line)
    diff = np.diff(line.astype(np.float64))
    if isinstance(size, tuple):
        zero_elements = np.isclose(diff, np.zeros_like(diff))
        edges_in_range = np.logical_and(diff > size[0], diff < size[1])
        valid_edges = np.logical_and(edges_in_range,
                                     np.logical_not(zero_elements))
        indices = np.where(valid_edges)
    else:
        indices = np.where(diff > size)
    indices = tuple(i + 1 for i in indices)
    return indices

def correct_split_edges(line):
    """
    This function corrects edges that are spread over multiple ticks.

    Parameters
    ----------
    line : numpy array
        The line for which the edges should be corrected.

    Returns
    -------
    line : numpy array
        Line with corrected edges.

    Examples
    --------
    >>> import numpy as np
    >>> import utils2p.synchronization
    >>> line = np.array([0, 0, 0, 1, 2, 3, 3, 3, 2, 1, 0, 0, 0])
    >>> utils2p.synchronization.correct_split_edges(line)
    array([0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 0, 0, 0])
    """
    rising_edges = np.where(np.diff(line) > 0)[0] + 1
    falling_edges = np.where(np.diff(line) < 0)[0]

    split_rising_edges = np.where(np.diff(rising_edges) == 1)[0]
    split_falling_edges = np.where(np.diff(falling_edges) == 1)[0]

    if len(split_rising_edges) == 0 and len(split_falling_edges) == 0:
        return line

    first_halfs_rising = rising_edges[split_rising_edges]
    second_halfs_rising = rising_edges[split_rising_edges + 1]
    line[first_halfs_rising] = line[second_halfs_rising]

    first_halfs_falling = falling_edges[split_falling_edges]
    second_halfs_falling = falling_edges[split_falling_edges + 1]
    line[second_halfs_falling] = line[first_halfs_falling]

    # Recursive to get edges spread over more than two ticks
    return correct_split_edges(line)

def get_start_times(line, times, zero_based_counter=False):
    """
    Get the start times of a digital signal,
    i.e. the times of the rising edges.
    If the line is a zero based counter, such as the processed
    `frame_counter` or the processed `cam_line`, there is a
    possibility that the first element in line is already zero.
    This corresponds to the case where the acquisition of the
    first frame was triggered before ThorSync started.
    If `zero_based_counter` is `False` this frame will be
    dropped, i.e. no time for the frame is returned, since
    there is no rising edge corresponding to the frame.

    Parameters
    ----------
    line : numpy array
        Line signal from h5 file.
    times : numpy array
        Times returned by :func:`utils2p.synchronization.get_times`
    zero_based_counter : boolean
        Indicates whether the line is a zero based counter.

    Returns
    -------
    time_points : list
        List of the start times.

    Examples
    --------
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> binary_line = np.array([0, 1, 1, 0, 1, 1])
    >>> times = utils2p.synchronization.get_times(len(binary_line), freq=20)
    >>> times
    array([0.  , 0.05, 0.1 , 0.15, 0.2 , 0.25])
    >>> utils2p.synchronization.get_start_times(binary_line, times)
    array([0.05, 0.2 ])
    """
    indices = edges(line, size=(0, np.inf))
    if zero_based_counter and line[0] >= 0:
        if line[0] > 0:
            warnings.warn(f"The counter start with value {line[0]}")
        indices_with_first_frame = np.zeros(len(indices[0]) + 1, dtype=int)
        indices_with_first_frame[1:] = indices[0]
        indices = (indices_with_first_frame, )
    time_points = times[indices]
    return time_points

def _capture_metadata(n_frames, dropped_frames=None):
    """
    Returns a dictionary as it is usually saved by the seven
    camera setup in the "capture_metadata.json" file.
    It assumes that no frames where dropped.

    Parameters
    ----------
    n_frames : list of integers
        Number of frames for each camera.
    dropped_frames : list of list of integers
        Frames that were dropped for each camera.
        Default is None which means no frames where
        dropped.

    Returns
    -------
    capture_info : dict
        Default metadata dictionary for the seven camera
        system.
    """
    if dropped_frames is None:
        dropped_frames = [[] for i in range(len(n_frames))]
    capture_info = {"Frame Counts": {}}
    for cam_idx, n in enumerate(n_frames):
        frames_dict = {}
        current_frame = 0
        for i in range(n):
            while current_frame in dropped_frames[cam_idx]:
                current_frame += 1
            frames_dict[str(i)] = current_frame
            current_frame += 1
        capture_info["Frame Counts"][str(cam_idx)] = frames_dict
    return capture_info

def process_cam_line(line, seven_camera_metadata):
    """
    Removes superfluous signals and uses frame numbers in array.
    The cam line signal form the h5 file is a binary sequence.
    Rising edges mark the acquisition of a new frame.
    The setup keeps producing rising edges after the acquisition of the
    last frame. These rising edges are ignored.
    This function converts the binary line to frame numbers using the
    information stored in the metadata file of the seven camera setup.
    In the metadata file the keys are the indices of the file names
    and the values are the grabbed frame numbers. Suppose the 3
    frame was dropped. Then the entries in the dictionary will
    be as follows:
    "2": 2
    "3": 4
    "4": 5

    Parameters
    ----------
    line : numpy array
        Line signal from h5 file.
    seven_camera_metadata : string
        Path to the json file saved by our camera software.
        This file is usually located in the same folder as the frames
        and is called 'capture_metadata.json'. If None, it is assumed
        that no frames were dropped.

    Returns
    -------
    processed_line : numpy array
        Array with frame number for each time point.
        If no frame is available for a given time,
        the value is -9223372036854775808.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> seven_camera_metadata = utils2p.find_seven_camera_metadata_file("data/mouse_kidney_raw")
    >>> line_names = ["Basler"]
    >>> (cam_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> set(np.diff(cam_line))
    {0, 8, 4294967288}
    >>> processed_cam_line = utils2p.synchronization.process_cam_line(cam_line, seven_camera_metadata)
    >>> set(np.diff(processed_cam_line))
    {0, 1, -9223372036854775808, 9223372036854775749}
    >>> cam_line = np.array([0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0])
    >>> utils2p.synchronization.process_cam_line(cam_line, seven_camera_metadata=None)
    array([-9223372036854775808,                    0,                    0,
                              0,                    0,                    0,
                              1,                    1,                    1,
                              1,                    1])
    """
    # Check that sequence is binary
    if len(set(line)) > 2:
        raise ValueError("Invalid line argument. Sequence is not binary.")

    # Find indices of the start of each frame acquisition
    rising_edges = edges(line, (0, np.inf))[0]

    # Load capture metadata or generate default
    if seven_camera_metadata is not None:
        with open(seven_camera_metadata, "r") as f:
            capture_info = json.load(f)
    else:
        capture_info = _capture_metadata([
            len(rising_edges),
        ])

    # Find the number of frames for each camera
    n_frames = []
    for cam_idx in capture_info["Frame Counts"].keys():
        max_in_json = max(capture_info["Frame Counts"][cam_idx].values())
        n_frames.append(max_in_json + 1)

    # Ensure all cameras acquired the same number of frames
    if len(np.unique(n_frames)) > 1:
        raise SynchronizationError(
            "The frames across cameras are not synchronized.")

    # Last rising edge that corresponds to a frame
    last_tick = max(n_frames)

    # check that there is a rising edge for every frame
    if len(rising_edges) < last_tick:
        raise ValueError(
            "The provided cam line and metadata are inconsistent. " +
            "cam line has less frame acquisitions than metadata.")

    # Ensure correct handling if no rising edges are present after last frame
    if len(rising_edges) == int(last_tick):
        average_frame_length = int(np.mean(np.diff(rising_edges)))
        last_rising_edge = rising_edges[-1]
        additional_edge = last_rising_edge + average_frame_length
        if additional_edge > len(line):
            additional_edge = len(line)
        rising_edges = list(rising_edges)
        rising_edges.append(additional_edge)
        rising_edges = np.array(rising_edges)

    processed_line = np.ones_like(line) * np.nan

    current_frame = 0
    first_camera_used = sorted(list(capture_info["Frame Counts"].keys()))[0]
    for i, (start, stop) in enumerate(
            zip(rising_edges[:last_tick], rising_edges[1:last_tick + 1])):
        if capture_info["Frame Counts"][first_camera_used][str(current_frame +
                                                               1)] <= i:
            current_frame += 1
        processed_line[start:stop] = current_frame
    return processed_line.astype(np.int)

def process_frame_counter(line, metadata=None, steps_per_frame=None):
    """
    Converts the frame counter line to an array with frame numbers for each
    time point.

    Parameters
    ----------
    line : numpy array
        Line signal from h5 file.
    metadata : :class:`utils2p.Metadata`
        :class:`utils2p.Metadata` object that holding the 2p imaging
        metadata for the experiment. Optional. If metadata is not
        given steps_per_frame has to be set.
    steps_per_frame : int
        Number of steps the frame counter takes per frame.
        This includes fly back frames and averaging, i.e. if you
        acquire one frame and flyback frames is set to 3 this number
        should be 4.

    Returns
    -------
    processed_frame_counter : numpy array
        Array with frame number for each time point.
        If no frame was recorded at a time point,
        the value is -9223372036854775808.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_z_stack")
    >>> line_names = ["Frame Counter",]
    >>> (frame_counter,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> set(frame_counter)
    {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30}
    >>> metadata_file = utils2p.find_metadata_file("data/mouse_kidney_z_stack")
    >>> metadata = utils2p.Metadata(metadata_file)
    >>> processed_frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, metadata)
    >>> set(processed_frame_counter)
    {0, -9223372036854775808}
    >>> steps_per_frame = metadata.get_n_z() * metadata.get_n_averaging()
    >>> steps_per_frame
    30
    >>> processed_frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, steps_per_frame=steps_per_frame)
    >>> set(processed_frame_counter)
    {0, -9223372036854775808}

    By default the function treat volumes as frames.
    If you want to treat every slice of the volume as a separate frame,
    you can do so by `steps_per_frame`. The example has three steps in z.

    >>> steps_per_frame = metadata.get_n_averaging()
    >>> steps_per_frame
    10
    >>> processed_frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, steps_per_frame=steps_per_frame)
    >>> set(processed_frame_counter)
    {0, 1, 2, -9223372036854775808}
    """
    if metadata is not None and steps_per_frame is not None:
        warnings.warn("metadata argument will be ignored " +
                      "because steps_per_frame argument was set.")
    if metadata is not None and not isinstance(metadata, main.Metadata):
        raise TypeError(
            "metadata argument must be of type utils2p.Metadata or None.")
    if steps_per_frame is not None and not isinstance(steps_per_frame, int):
        raise TypeError(f"steps_per_frame has to be of type int not {type(steps_per_frame)}")

    if metadata is not None and steps_per_frame is None:
        if metadata.get_value("Streaming", "zFastEnable") == "0":
            steps_per_frame = 1
        else:
            steps_per_frame = metadata.get_n_z()
            if metadata.get_value("Streaming", "enable") == "1":
                steps_per_frame += metadata.get_n_flyback_frames()
        if metadata.get_value(
                "LSM",
                "averageMode") == "1" and metadata.get_area_mode() not in [
                    "line", "kymograph"
                ]:
            steps_per_frame = steps_per_frame * metadata.get_n_averaging()
    elif steps_per_frame is None:
        raise ValueError("If no metadata object is given, " +
                         "the steps_per_frame argument has to be set.")

    processed_frame_counter = np.ones_like(line) * np.nan
    rising_edges = edges(line, (0, np.inf))[0]

    # Case of one frame/volume only
    if len(rising_edges) <= steps_per_frame:
        processed_frame_counter[rising_edges[0]:] = 0
        return processed_frame_counter.astype(np.int)

    for i, index in enumerate(
            range(0,
                  len(rising_edges) - steps_per_frame, steps_per_frame)):
        processed_frame_counter[
            rising_edges[index]:rising_edges[index + steps_per_frame]] = i
    processed_frame_counter[rising_edges[-1 * steps_per_frame]:] = (
        processed_frame_counter[rising_edges[-1 * steps_per_frame] - 1] + 1)
    return processed_frame_counter.astype(np.int)

def process_stimulus_line(line):
    """
    This function converts the stimulus line to an array with
    0s and 1s for stimulus off and on respectively. The raw
    stimulus line can contain values larger than 1.

    Parameters
    ----------
    line : numpy array
        Line signal from h5 file.

    Returns
    -------
    processed_frame_counter : numpy array
        Array with binary stimulus state for each time point.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> line_names = ["CO2_Stim"]
    >>> (stimulus_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> set(stimulus_line)
    {0, 4}
    >>> processed_stimulus_line = utils2p.synchronization.process_stimulus_line(stimulus_line)
    >>> set(processed_stimulus_line)
    {0, 1}
    """
    processed_stimulus_line = np.zeros_like(line)
    indices = np.where(line > 0)
    processed_stimulus_line[indices] = 1
    return processed_stimulus_line.astype(np.int)

def process_optical_flow_line(line):
    """
    This function converts the optical flow line
    into a step function. The value corresponds
    to the index of optical flow value at this
    time point. If the value is -9223372036854775808, no optical flow
    value was recorded for this time point.

    Note: Due to the time it takes to transfer the data
    from the Arduino to the computer it is possible that
    the last optical flow data point is missing, i.e.
    the processed optical flow line indicates one more
    data point than the text file contains. This can be
    solved by cropping all lines before the acquisition
    of the last optical flow data point. Lines can be
    cropped with :func:`crop_lines`.

    Parameters
    ----------
    line : numpy array
        Line signal for h5 file.

    Returns
    -------
    processed_optical_flow_line : numpy array
        Array with monotonically increasing step
        function.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> line_names = ["OpFlow"]
    >>> (optical_flow_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> set(optical_flow_line)
    {0, 16}
    >>> processed_optical_flow_line = utils2p.synchronization.process_optical_flow_line(optical_flow_line)
    >>> len(set(processed_optical_flow_line))
    1409
    """
    processed_optical_flow_line = np.ones_like(line) * np.nan
    rising_edges = edges(line, (0, np.inf))[0]
    for i in range(0, len(rising_edges) - 1):
        processed_optical_flow_line[rising_edges[i]:rising_edges[i + 1]] = i
    processed_optical_flow_line[rising_edges[-1]:] = (
        processed_optical_flow_line[rising_edges[-1] - 1] + 1)
    return processed_optical_flow_line.astype(np.int)

def crop_lines(mask, lines):
    """
    This function crops all lines based on a binary signal/mask.
    The 'Capture On' line of the h5 file can be used as a mask.

    Parameters
    ----------
    mask : numpy array
        Mask that is used for cropping.
    lines : list of numpy arrays
        List of the lines that should be cropped.

    Returns
    -------
    cropped_lines : tuple of numpy arrays
        Tuple of cropped lines in same order as in input list.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> line_names = ["Frame Counter", "Capture On", "CO2_Stim", "OpFlow"]
    >>> (frame_counter, capture_on, stimulus_line, optical_flow_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, steps_per_frame=4)
    >>> len(frame_counter), len(capture_on), len(stimulus_line), len(optical_flow_line)
    (117000, 117000, 117000, 117000)
    >>> mask = np.logical_and(frame_counter >= 0, capture_on)
    >>> np.sum(mask)
    105869
    >>> (frame_counter, capture_on, stimulus_line, optical_flow_line,) = utils2p.synchronization.crop_lines(mask, (frame_counter, capture_on, stimulus_line, optical_flow_line,))
    >>> len(frame_counter), len(capture_on), len(stimulus_line), len(optical_flow_line)
    (105869, 105869, 105869, 105869)
    >>> line = np.arange(10)
    >>> mask = np.ones(10, dtype=np.bool)
    >>> mask[0] = False
    >>> mask[-1] = False
    >>> mask[4] = False
    >>> utils2p.synchronization.crop_lines(mask, (line,))
    (array([1, 2, 3, 4, 5, 6, 7, 8]),)
    """
    indices = np.where(mask)[0]
    first_idx = indices[0]
    last_idx = indices[-1]
    cropped_lines = []
    for line in lines:
        cropped_lines.append(line[first_idx:last_idx + 1])
    return tuple(cropped_lines)

def beh_idx_to_2p_idx(beh_indices, cam_line, frame_counter):
    """
    This functions converts behaviour frame numbers into the corresponding
    2p frame numbers.

    Parameters
    ----------
    beh_indices : numpy array
        Indices of the behaviour frames to be converted.
    cam_line : numpy array
        Processed cam line.
    frame_counter : numpy array
        Processed frame counter.

    Returns
    -------
    indices_2p : numpy array
        Corresponding 2p frame indices.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> line_names = ["Frame Counter", "Basler"]
    >>> (frame_counter, cam_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, steps_per_frame=4)
    >>> seven_camera_metadata = utils2p.find_seven_camera_metadata_file("data/mouse_kidney_raw")
    >>> cam_line = utils2p.synchronization.process_cam_line(cam_line, seven_camera_metadata)
    >>> utils2p.synchronization.beh_idx_to_2p_idx(np.array([0,]), cam_line, frame_counter)
    array([-9223372036854775808])
    >>> utils2p.synchronization.beh_idx_to_2p_idx(np.array([10,]), cam_line, frame_counter)
    array([0])
    >>> utils2p.synchronization.beh_idx_to_2p_idx(np.arange(30), cam_line, frame_counter)
    array([-9223372036854775808,                    0,                    0,
                              0,                    0,                    0,
                              0,                    0,                    0,
                              0,                    0,                    0,
                              0,                    0,                    0,
                              0,                    0,                    0,
                              0,                    0,                    0,
                              0,                    1,                    1,
                              1,                    1,                    1,
                              1,                    1,                    1])
    """
    thor_sync_indices = edges(cam_line)[0]
    if not cam_line[0] < 0:
        thor_sync_indices = np.append(np.array([0]), thor_sync_indices)

    indices_2p = np.ones(len(beh_indices), dtype=np.int) * np.nan

    first_frame_of_cam_line = np.min(cam_line[np.where(cam_line >= 0)])

    for i, frame_num in enumerate(beh_indices):

        # This is necessary for cropped lines that don't start at 0
        frame_num = frame_num - first_frame_of_cam_line
        if frame_num < 0:
            raise ValueError(f"{frame_num + first_frame_of_cam_line} is smaller than first frame in cam_line ({first_frame_of_cam_line})")

        thor_sync_index = thor_sync_indices[frame_num]
        indices_2p[i] = frame_counter[thor_sync_index]

    return indices_2p.astype(np.int)

def reduce_during_2p_frame(frame_counter, values, function):
    """
    Reduces all values occurring during the acquisition of a
    2-photon frame to a single value using the `function`
    given by the user.

    Parameters
    ----------
    frame_counter : numpy array
        Processed frame counter.
    values : numpy array
        Values upsampled to the frequency of ThorSync,
        i.e. 1D numpy array of the same length as
        `frame_counter`.
    function : function
        Function used to reduce the value,
        e.g. np.mean.

    Returns
    -------
    reduced : numpy array
        Numpy array with value for each 2p frame.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> line_names = ["Frame Counter", "CO2_Stim"]
    >>> (frame_counter, stimulus_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, steps_per_frame=1)
    >>> stimulus_line = utils2p.synchronization.process_stimulus_line(stimulus_line)
    >>> np.max(frame_counter)
    4
    >>> stimulus_during_2p_frames = utils2p.synchronization.reduce_during_2p_frame(frame_counter, stimulus_line, np.mean)
    >>> len(stimulus_during_2p_frames)
    5
    >>> np.max(stimulus_during_2p_frames)
    0.7136134613556422
    >>> stimulus_during_2p_frames = utils2p.synchronization.reduce_during_2p_frame(frame_counter, stimulus_line, np.max)
    >>> len(stimulus_during_2p_frames)
    5
    >>> set(stimulus_during_2p_frames)
    {0.0, 1.0}
    """
    warnings.warn(
        "reduce_during_2p_frame is deprecated use reduce_during_frame instead",
        DeprecationWarning)
    return reduce_during_frame(frame_counter, values, function)

def reduce_during_frame(line, values, function):
    """
    Reduces all values occurring during the acquisition of a
    frame to a single value using the `function` given by the user.
    The line function should be of the resolution of
    the ThorSync ticks and have the frame index as values.
    Possible choices are the processed frame_counter line or the
    processed cam_line.

    Parameters
    ----------
    line : numpy array
        Line holding frame indices.
    values : numpy array
        Values upsampled to the frequency of ThorSync,
        i.e. 1D numpy array of the same length as
        `frame_counter`.
    function : function
        Function used to reduce the value,
        e.g. np.mean.

    Returns
    -------
    reduced : numpy array
        Numpy array with value for each 2p frame.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> import numpy as np
    >>> h5_file = utils2p.find_sync_file("data/mouse_kidney_raw")
    >>> line_names = ["Frame Counter", "CO2_Stim"]
    >>> (frame_counter, stimulus_line,) = utils2p.synchronization.get_lines_from_h5_file(h5_file, line_names)
    >>> frame_counter = utils2p.synchronization.process_frame_counter(frame_counter, steps_per_frame=1)
    >>> stimulus_line = utils2p.synchronization.process_stimulus_line(stimulus_line)
    >>> np.max(frame_counter)
    4
    >>> stimulus_during_2p_frames = utils2p.synchronization.reduce_during_frame(frame_counter, stimulus_line, np.mean)
    >>> len(stimulus_during_2p_frames)
    5
    >>> np.max(stimulus_during_2p_frames)
    0.7136134613556422
    >>> stimulus_during_2p_frames = utils2p.synchronization.reduce_during_frame(frame_counter, stimulus_line, np.max)
    >>> len(stimulus_during_2p_frames)
    5
    >>> set(stimulus_during_2p_frames)
    {0.0, 1.0}
    """
    if len(line) != len(values):
        raise ValueError("line and values need to have the same length.")

    thor_sync_indices = tuple(edges(line, (0, np.inf))[0])

    starts = thor_sync_indices
    stops = thor_sync_indices[1:] + (len(line), )

    if not line[0] == -9223372036854775808:
        starts = (0, ) + starts
        stops = (thor_sync_indices[0], ) + stops

    dtype = values.dtype
    if np.issubdtype(dtype, np.number):
        dtype = np.float
    else:
        dtype = np.object
    reduced = np.empty(len(starts), dtype=dtype)

    for i, (start, stop) in enumerate(zip(starts, stops)):
        reduced[i] = function(values[start:stop])

    return reduced

class SyncMetadata(main._XMLFile):
    """
    Class for managing ThorSync metadata.
    Loads metadata file 'ThorRealTimeDataSettings.xml'
    and returns the root of an ElementTree.

    Parameters
    ----------
    path : string
        Path to xml file.

    Returns
    -------
    Instance of class Metadata
        Based on given xml file.

    Examples
    --------
    >>> import utils2p.synchronization
    >>> metadata = utils2p.synchronization.SyncMetadata("data/mouse_kidney_raw/2p/Sync-025/ThorRealTimeDataSettings.xml")
    >>> type(metadata)
    <class 'utils2p.synchronization.SyncMetadata'>
    """
    def get_active_devices(self):
        active_devices = []
        for device in self.get_value("DaqDevices", "AcquireBoard"):
            if device.attrib["active"] == "1":
                active_devices.append(device)
        return active_devices

    def get_freq(self):
        """
        Returns the frequency of the ThorSync
        value acquisition, i.e. the sample rate.

        Returns
        -------
        freq : integer
            Sample frequency in Hz.

        Examples
        --------
        >>> import utils2p.synchronization
        >>> metadata = utils2p.synchronization.SyncMetadata("data/mouse_kidney_raw/2p/Sync-025/ThorRealTimeDataSettings.xml")
        >>> metadata.get_freq()
        30000
        """
        sample_rate = -1
        for device in self.get_active_devices():
            set_for_device = False
            for element in device.findall("SampleRate"):
                if element.attrib["enable"] == "1":
                    if set_for_device:
                        raise ValueError(
                            "Invalid metadata file. Multiple sample rates " +
                            f"are enabled for device {device.type}")
                    if sample_rate != -1:
                        raise ValueError("Multiple devices are enabled.")
                    sample_rate = int(element.attrib["rate"])
                    set_for_device = True
        return sample_rate

def get_processed_lines(sync_file,
                        sync_metadata_file,
                        metadata_2p_file,
                        seven_camera_metadata_file=None):
    """
    This function extracts all the standard lines and processes them.
    It works for both microscopes.

    Parameters
    ----------
    sync_file : str
        Path to the synchronization file.
    sync_metadata_file : str
        Path to the synchronization metadata file.
    metadata_2p_file : str
        Path to the ThorImage metadata file.
    seven_camera_metadata_file : str
        Path to the metadata file of the 7 camera system.

    Returns
    -------
    processed_lines : dictionary
        Dictionary with all processed lines.

    Examples
    --------
    >>> import utils2p
    >>> import utils2p.synchronization
    >>> experiment_dir = "data/mouse_kidney_raw/"
    >>> sync_file = utils2p.find_sync_file(experiment_dir)
    >>> metadata_file = utils2p.find_metadata_file(experiment_dir)
    >>> sync_metadata_file = utils2p.find_sync_metadata_file(experiment_dir)
    >>> seven_camera_metadata_file = utils2p.find_seven_camera_metadata_file(experiment_dir)
    >>> processed_lines = utils2p.synchronization.get_processed_lines(sync_file, sync_metadata_file, metadata_file, seven_camera_metadata_file)
    """
    processed_lines = {}
    processed_lines["Capture On"], processed_lines[
        "Frame Counter"] = get_lines_from_sync_file(
            sync_file, ["Capture On", "Frame Counter"])

    try:
        # For microscope 1
        processed_lines["CO2"], processed_lines["Cameras"], processed_lines[
            "Optical flow"] = get_lines_from_sync_file(sync_file, [
                "CO2_Stim",
                "Basler",
                "OpFlow",
            ])
    except KeyError:
        # For microscope 2
        processed_lines["CO2"], processed_lines[
            "Cameras"] = get_lines_from_h5_file(sync_file, [
                "CO2",
                "Cameras",
            ])

    processed_lines["Cameras"] = process_cam_line(processed_lines["Cameras"],
                                                  seven_camera_metadata_file)

    metadata_2p = main.Metadata(metadata_2p_file)
    processed_lines["Frame Counter"] = process_frame_counter(
        processed_lines["Frame Counter"], metadata_2p)

    processed_lines["CO2"] = process_stimulus_line(processed_lines["CO2"])

    if "Optical flow" in processed_lines.keys():
        processed_lines["Optical flow"] = process_optical_flow_line(
            processed_lines["Optical flow"])

    mask = np.logical_and(processed_lines["Capture On"],
                          processed_lines["Frame Counter"] >= 0)

    # Make sure the clipping start just before the
    # acquisition of the first frame
    indices = np.where(mask)[0]
    mask[max(0, indices[0] - 1)] = True

    for line_name, _ in processed_lines.items():
        processed_lines[line_name] = crop_lines(mask, [
            processed_lines[line_name],
        ])[0]

    # Get times of ThorSync ticks
    metadata = SyncMetadata(sync_metadata_file)
    freq = metadata.get_freq()
    times = get_times(len(processed_lines["Frame Counter"]), freq)
    processed_lines["Times"] = times

    return processed_lines

def epoch_length_filter(line, cut_off):
    """
    This function filters a binary based on the length
    of each event.

    Parameters
    ----------
    line : numpy array of type bool
        Binary trace that is filtered.
    cut_off : int
        The minimal event length. All event shorter
        than `cut_off` are set to `False`.

    Returns
    -------
    filtered : numpy array of type bool
        The filtered binary trace.
    """
    diff = np.diff(np.pad(line.astype(int), 1, "constant", constant_values=0))
    rising_edges = np.where(diff > 0)[0]
    falling_edges = np.where(diff < 0)[0]
    epoch_length = falling_edges - rising_edges

    discarded_epochs = (epoch_length < cut_off)

    discarded_rising_edges = rising_edges[discarded_epochs]
    discarded_falling_edges = falling_edges[discarded_epochs]

    filtered = line.copy()
    for start, stop in zip(discarded_rising_edges, discarded_falling_edges):
        filtered[start:stop] = 0

    return filtered.astype(bool)

def process_odor_line(line,
                      freq=30000,
                      arduino_commands=(
                          "None",
                          "Odor1",
                          "Odor2",
                          "Odor3",
                          "Odor4",
                          "Odor5",
                          "Odor6",
                          "Odor1R",
                          "Odor2R",
                          "Odor1L",
                          "Odor2L",
                          "Odor1B",
                          "Odor2B",
                          "WaterB",
                          "bubbleMFC_R0",
                          "MFC1_R2",
                          "MFC2_L1",
                      ),
                      step_size=0.2703,
                      filter_only=False):
    """
    The odor line is based on a PWM signal for the Arduino controlling the
    valves. This function applies a Butterworth filter and converts the
    resulting voltages to level indices. The corresponding the setting of the
    valves are given by the `arduino_commands` argument.

    Parameters
    ----------
    line : numpy array
        Unprocessed odor line from h5 file.
    freq : int
        Frequency of ThorSync. Necessary for the Butterworth filter.
    arduino_commands : list of strings
        Description of the valve settings for commands sent to arduino.
        Note: The order matters since the serial communications between
        computer and Arduino is based on the index in the list.
        This index is converted to a PWM signal that is recorded by ThorSync.
    step_size : float
        The voltage step size between different levels of the PWM. This is used
        to convert the voltage to indices.
    filter_only : bool
        If `True`, only the filtered line is returned instead of the odors
        based on the `arduino_commands`. This is useful for determining
        the `step_size`.

    Returns
    -------
    numpy array of strings
    """
    b, a = scipy.signal.butter(3, 10, fs=freq)
    filtered_line = scipy.signal.filtfilt(b, a, line)
    if filter_only:
        return filtered_line
    indices = np.rint(filtered_line / step_size).astype(int)
    for index in np.unique(indices):
        mask = (indices == index)
        filtered_mask = epoch_length_filter(mask, freq)
        indices[mask & ~filtered_mask] = 0
    return np.array(arduino_commands)[indices]

def event_based_frame_indices(event_indicator):
    """
    Calculates frame indices based on events.
    Frames before an event have negative numbers.
    The event onset has frame number 0 and the frames
    count up for the duration of the event.
    To be able to distinguish multiple events in the
    `event_indicator` an array with event numbers is
    returned.

    Parameters
    ----------
    event_indicator : numpy array of type bool
        True indicates some event happening.

    Returns
    -------
    event_based_indices : numpy array of type int
        Event based indices as described above.
    event_number : numpy array of type int
        Array of the same length as `event_based_indices`
        counting the number of events in event indicator.
    """
    mask = event_indicator.astype(bool)
    inv_mask = ~mask
    inv_mask = inv_mask[::-1]
    mask = mask.astype(np.int8)
    inv_mask = inv_mask.astype(np.int8)
    mask = np.concatenate(([
        0,
    ], mask))
    inv_mask = np.concatenate((inv_mask, [
        0,
    ]))

    event_numbers = np.cumsum(np.clip(np.diff(mask), 0, None))
    inv_event_numbers = np.cumsum(np.clip(np.diff(inv_mask), 0, None))

    mask = mask[1:]
    inv_mask = inv_mask[:-1]

    # Count up from zero during the event
    event_frame_indices = np.cumsum(mask)
    inv_event_frame_indices = np.cumsum(inv_mask)
    n_events = max(event_numbers)
    for event in np.arange(1, n_events + 1):
        i = np.where(event_numbers == event)
        event_frame_indices[i] = event_frame_indices[i] - \
            event_frame_indices[i[0][0]]
    event_frame_indices[~mask.astype(np.bool)] = 0

    # Count down from zero before each event
    n_inv_event = max(inv_event_numbers)
    for inv_event in np.arange(1, n_inv_event + 1):
        i = np.where(inv_event_numbers == inv_event)
        inv_event_frame_indices[i] = inv_event_frame_indices[i] - \
            inv_event_frame_indices[i[0][0]]
    inv_event_frame_indices[~inv_mask.astype(np.bool)] = 0
    inv_event_frame_indices = -inv_event_frame_indices[::-1]

    event_frame_indices[~mask.astype(bool)] = inv_event_frame_indices[
        ~mask.astype(bool)]

    event_numbers = np.cumsum(
        -1 * np.clip(np.diff(np.concatenate(
            ([2], event_frame_indices))), -1, 0))
    # Make sure the last frames are not counted as the pre-event
    # frames of a new event
    n_events = max(event_numbers)
    last_event = np.where(event_numbers == n_events)
    if np.all(event_frame_indices[last_event] < 0):
        event_numbers[last_event] = -1

    return event_frame_indices, event_numbers

ext: py

reason: import numpy,import scipy


text:
# coding=utf-8
# name=hu_yang_jie
#coding=utf-8
import cv2
import numpy as np

img = cv2.imread("bili.jpg")    #载入图像
h, w = img.shape[:2]      #获取图像的高和宽
cv2.imshow("Origin", img)     #显示原始图像

blured = cv2.blur(img,(5,5))    #进行滤波去掉噪声
cv2.imshow("Blur", blured)     #显示低通滤波后的图像

mask = np.zeros((h+2, w+2), np.uint8)  #掩码长和宽都比输入图像多两个像素点，满水填充不会超出掩码的非零边缘
#进行泛洪填充
cv2.floodFill(blured, mask, (w-1,h-1), (255,255,255), (2,2,2),(3,3,3),8)
cv2.imshow("floodfill", blured)

#得到灰度图
gray = cv2.cvtColor(blured,cv2.COLOR_BGR2GRAY)
cv2.imshow("gray", gray)

#定义结构元素
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(50, 50))
#开闭运算，先开运算去除背景噪声，再继续闭运算填充目标内的孔洞
opened = cv2.morphologyEx(gray, cv2.MORPH_OPEN, kernel)
closed = cv2.morphologyEx(opened, cv2.MORPH_CLOSE, kernel)
cv2.imshow("closed", closed)

#求二值图
ret, binary = cv2.threshold(closed,250,255,cv2.THRESH_BINARY)
cv2.imshow("binary", binary)

#找到轮廓
_,contours, hierarchy = cv2.findContours(binary,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
#绘制轮廓

cv2.drawContours(img,contours,-1,(0,0,255),3)
#绘制结果
cv2.imshow("result", img)

cv2.waitKey(0)
cv2.destroyAllWindows()

reason: import numpy


text:
import cv2
import os
import math
import torch
import numpy as np
from PIL import Image
from torchvision import transforms
from torch.utils.data import Dataset

RANDOM_CROP = 0.5
MIN_AREA = 256
SIGMA = 0.02

class ExtremePointsDataset(Dataset):

    def __init__(self, coco, image_folder, is_training=False, training_size=None):
        """
        During the training I use square images.
        But during the evaluation I only pad images with zeros.
        And I want all my image sizes be divisible by 32.

        Arguments:
            coco: an instance of COCO.
            image_folder: a string.
            is_training: a boolean.
            training_size: an integer or None.
        """
        self.labels = coco.getCatIds(catNms=['person'])
        self.ids = coco.getImgIds(catIds=self.labels)
        self.coco = coco
        self.image_folder = image_folder
        self.is_training = is_training

        if is_training:
            assert training_size % 32 == 0
            self.training_size = training_size

        color_transforms = transforms.RandomChoice([
            transforms.ColorJitter(brightness=0.4, contrast=0.4, saturation=0.2, hue=0.1),
            transforms.RandomGrayscale(p=0.01)
        ])
        self.color_transforms = color_transforms

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

    def __getitem__(self, index):
        """
        Returns:
            image: a float tensor with shape [3, h, w].
                It represents a RGB image with
                pixel values in the range [0, 1].
            heatmaps: a float tensor with shape [5, h/4, w/4].
            offsets: a float tensor with shape [10, h/4, w/4].
            masks: a float tensor with shape [2, h/4, w/4].
            num_boxes: a long tensor with shape [].
        """

        # LOAD AN IMAGE AND ANNOTATIONS

        image_id = self.ids[index]
        filename = self.coco.loadImgs(image_id)[0]['file_name']
        image = cv2.imread(os.path.join(self.image_folder, filename))

        # sometimes images are gray
        if len(image.shape) == 2:
            image = cv2.cvtColor(image, cv2.COLOR_GRAY2BGR)

        image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
        height, width, _ = image.shape

        # load all person annotations
        annotations_ids = self.coco.getAnnIds(imgIds=image_id, catIds=self.labels, iscrowd=None)
        annotations = self.coco.loadAnns(annotations_ids)

        # OPTIONALLY DO A RANDOM CROP

        if self.is_training:

            # choose a random window
            min_dimension = min(height, width)
            min_size = int(min_dimension * RANDOM_CROP)
            size = np.random.randint(min_size, min_dimension)
            x = np.random.randint(0, width - size)
            y = np.random.randint(0, height - size)
            height, width = size, size  # new size

            ymin, xmin, ymax, xmax = y, x, y + size, x + size
            image = image[ymin:ymax, xmin:xmax]
        else:
            ymin, xmin, ymax, xmax = 0, 0, height, width
            # it is whole image window

        # DETECT EXTREME POINTS AND COLLECT ANNOTATIONS INTO ARRAYS

        segmentation_mask = np.zeros([height, width], dtype='bool')
        not_ignore_mask = np.ones([height, width], dtype='bool')  # mask for loss

        extreme_points = []
        for a in annotations:

            mask = self.coco.annToMask(a)
            mask = mask[ymin:ymax, xmin:xmax]

            unannotated = a['iscrowd'] == 1
            too_small = mask.sum() < MIN_AREA  # whether area is too small

            if too_small or unannotated:
                not_ignore_mask = np.logical_and(mask == 0, not_ignore_mask)
                continue

            extreme_points.append(get_extreme_points(mask))
            segmentation_mask = np.logical_or(segmentation_mask, mask == 1)

        # number of persons on the image
        num_boxes = len(extreme_points)

        if num_boxes == 0:
            extreme_points = np.empty([0, 4, 2], dtype='int32')
        else:
            extreme_points = np.stack(extreme_points, axis=0)
            # it has shape [num_boxes, 4, 2]

        masks = np.stack([segmentation_mask, not_ignore_mask], axis=2)  # shape [height, width, 2]
        image_and_masks = np.concatenate([image, masks], axis=2)  # shape [height, width, 5]

        # RESIZE THE IMAGE AND ANNOTATIONS (AND OPTIONALLY AUGMENT DATA)

        if self.is_training:

            image_and_masks, extreme_points = random_flip(image_and_masks, extreme_points)
            image_and_masks, extreme_points = random_pad(image_and_masks, extreme_points)

            height, width, _ = image_and_masks.shape  # it could change
            w, h = self.training_size, self.training_size

            image_and_masks = cv2.resize(image_and_masks, (w, h), cv2.INTER_NEAREST)
            scaler = np.array([w/width, h/height])
            extreme_points = (extreme_points * scaler).astype('int32')
        else:
            # pad the image and masks with zeros
            w, h = 32 * math.ceil(width/32), 32 * math.ceil(height/32)
            empty = np.zeros((h, w, 5), dtype='uint8')
            empty[:height, :width] = image_and_masks
            image_and_masks = empty

        # DOWNSAMPLE ANNOTATIONS AND CREATE HEATMAPS FOR REGRESSION

        # find centers
        t, b, l, r = np.split(extreme_points, 4, axis=1)
        cx, cy = (l[:, :, 0] + r[:, :, 0]) // 2, (t[:, :, 1] + b[:, :, 1]) // 2
        # they have shape [num_boxes, 1]

        center = np.stack([cx, cy], axis=2)  # shape [num_boxes, 1, 2]
        extreme_points = np.concatenate([extreme_points, center], axis=1)

        image, masks = np.split(image_and_masks, indices_or_sections=[3], axis=2)
        w, h = w // 4, h // 4  # downsampled size
        masks = cv2.resize(masks, (w, h), cv2.INTER_NEAREST)
        offsets = (extreme_points / 4.0) - (extreme_points // 4)  # shape [num_boxes, 5, 2]
        extreme_points = extreme_points // 4  # floor

        heights = extreme_points[:, 1, 1] - extreme_points[:, 0, 1]
        widths = extreme_points[:, 3, 0] - extreme_points[:, 2, 0]
        sigmas = SIGMA * np.maximum(heights, widths)
        sigmas = np.clip(sigmas, 1.0, None)
        # they have shape [num_boxes]

        heatmaps, vectormaps = [], []
        for i in range(5):

            # get a particular point type
            points = extreme_points[:, i]

            heatmaps.append(generate_heatmap(points, h, w, sigmas))
            # don't forget to convert to (y, x) format:
            values = torch.FloatTensor(offsets[:, i, [1, 0]])
            index = torch.LongTensor(points[:, [1, 0]])
            vectormap = torch.sparse.FloatTensor(index.t(), values, [h, w, 2]).to_dense()
            vectormaps.append(vectormap)

        # CONVERT TO PYTORCH TENSORS

        if self.is_training:
            image = Image.fromarray(image)
            image = self.color_transforms(image)
            image = np.array(image)

        image = torch.FloatTensor(image).permute(2, 0, 1).div(255.0)
        heatmaps = torch.FloatTensor(np.stack(heatmaps, axis=0))
        offsets = torch.cat(vectormaps, dim=2).permute(2, 0, 1)
        masks = torch.FloatTensor(masks).permute(2, 0, 1)
        num_boxes = torch.tensor(num_boxes)

        labels = {
            'heatmaps': heatmaps, 'offsets': offsets,
            'masks': masks, 'num_boxes': num_boxes
        }
        return image, labels

def get_extreme_points(mask):
    """
    Also see the original implementation:
    https://github.com/xingyizhou/ExtremeNet/blob/master/tools/gen_coco_extreme_points.py

    Arguments:
        mask: a numpy uint8 array with shape [h, w].
    Returns:
        a numpy int array with shape [4, 2].
        It represents points [top, bottom, left, right] in format [x, y].
        Coordinates are in ranges [0, w - 1] and [0, h - 1].
    """
    threshold = 0.02

    points = np.where(mask > 0)
    points = np.asarray(points, dtype='int32').transpose()  # shape [num_points, 2]
    y, x = points[:, 0], points[:, 1]

    ymin, ymax = y.min(), y.max()
    xmin, xmax = x.min(), x.max()
    h, w = ymax - ymin, xmax - xmin

    # top
    near_x = x[np.abs(y - ymin) <= threshold * h]
    tx = (near_x.max() + near_x.min()) // 2
    ty = ymin

    # bottom
    near_x = x[np.abs(y - ymax) <= threshold * h]
    bx = (near_x.max() + near_x.min()) // 2
    by = ymax

    # left
    near_y = y[np.abs(x - xmin) <= threshold * w]
    lx = xmin
    ly = (near_y.max() + near_y.min()) // 2

    # right
    near_y = y[np.abs(x - xmax) <= threshold * w]
    rx = xmax
    ry = (near_y.max() + near_y.min()) // 2

    return np.array([[tx, ty], [bx, by], [lx, ly], [rx, ry]], dtype='int32')

def generate_heatmap(points, h, w, sigmas):
    """
    Arguments:
        points: a numpy int array with shape [num_points, 2].
            It is in the form [x, y].
            And it is assumed that `0 <= x < w` and `0 <= y < h`.
        h, w: integers.
        sigmas: a numpy float array with shape [num_points].
    Returns:
        a numpy float array with shape [h, w].
    """

    k = 11  # sigmas mustn't be too large!
    a = np.arange(-k, k + 1, dtype='float32')
    X, Y = np.meshgrid(a, a)  # they have shape [2*k + 1, 2*k + 1]
    D = 0.5 * (X**2 + Y**2)
    heatmap = np.zeros([h + 2*k, w + 2*k], dtype='float32')

    for (x, y), sigma in zip(points, sigmas):
        g = np.exp(-D/(sigma**2))
        ymin, ymax = y, y + 2*k + 1
        xmin, xmax = x, x + 2*k + 1
        heatmap[ymin:ymax, xmin:xmax] = np.maximum(g, heatmap[ymin:ymax, xmin:xmax])

    return heatmap[k:-k, k:-k]

def random_flip(image_and_masks, extreme_points):

    if np.random.rand() > 0.5:

        image_and_masks = cv2.flip(image_and_masks, 1)
        height, width, _ = image_and_masks.shape
        x, y = np.split(extreme_points, 2, axis=2)
        extreme_points = np.concatenate([width - x - 1, y], axis=2)

        # switch left and right
        correct_order = [0, 1, 3, 2]
        extreme_points = extreme_points[:, correct_order]

    return image_and_masks, extreme_points

def random_pad(image_and_masks, extreme_points):

    height, width, _ = image_and_masks.shape
    ymin = extreme_points[:, 0, 1]
    ymax = extreme_points[:, 1, 1]
    sizes = ymax - ymin  # shape [num_boxes]

    if len(sizes) == 0 or 0.25 * height > sizes.mean():
        # if mean person height is
        # smaller than quarter image height,
        # in other words if image
        # contains only small people
        return image_and_masks, extreme_points

    # scale the image by this number
    scaler = np.random.uniform(1.75, 2.75)

    # calculate new image size
    h = int(scaler * height)
    w = int(scaler * width)
    empty = np.zeros((h, w, 5), dtype='uint8')

    y = np.random.randint(0, h - height)
    x = np.random.randint(0, w - width)
    empty[y:(y + height), x:(x + width)] = image_and_masks

    image_and_masks = empty
    extreme_points += np.array([x, y])

    return image_and_masks, extreme_points

max_stars_repo_name: TropComplique/ExtremeNet

max_issues_repo_name: TropComplique/ExtremeNet

max_forks_repo_name: TropComplique/ExtremeNet


text:
# https://udemy.com/recommender-systems
# https://deeplearningcourses.com/recommender-systems
from __future__ import print_function, division
from builtins import range, input
# Note: you may need to update your version of future
# sudo pip install -U future

import pickle
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.utils import shuffle

from keras.models import Model
from keras.layers import Input, Embedding, Dot, Add, Flatten, Dense, Concatenate
from keras.layers import Dropout, BatchNormalization, Activation
from keras.regularizers import l2
from keras.optimizers import SGD, Adam

# load in the data
df = pd.read_csv('../large_files/movielens-20m-dataset/edited_rating.csv')

N = df.userId.max() + 1 # number of users
M = df.movie_idx.max() + 1 # number of movies

# split into train and test
df = shuffle(df)
cutoff = int(0.8*len(df))
df_train = df.iloc[:cutoff]
df_test = df.iloc[cutoff:]

# initialize variables
K = 10 # latent dimensionality
mu = df_train.rating.mean()
epochs = 15
reg = 0. # regularization penalty

# keras model
u = Input(shape=(1,))
m = Input(shape=(1,))
u_embedding = Embedding(N, K)(u) # (N, 1, K)
m_embedding = Embedding(M, K)(m) # (N, 1, K)

##### main branch
u_bias = Embedding(N, 1)(u) # (N, 1, 1)
m_bias = Embedding(M, 1)(m) # (N, 1, 1)
x = Dot(axes=2)([u_embedding, m_embedding]) # (N, 1, 1)
x = Add()([x, u_bias, m_bias])
x = Flatten()(x) # (N, 1)

##### side branch
u_embedding = Flatten()(u_embedding) # (N, K)
m_embedding = Flatten()(m_embedding) # (N, K)
y = Concatenate()([u_embedding, m_embedding]) # (N, 2K)
y = Dense(400)(y)
y = Activation('elu')(y)
# y = Dropout(0.5)(y)
y = Dense(1)(y)

##### merge
x = Add()([x, y])

model = Model(inputs=[u, m], outputs=x)
model.compile(
  loss='mse',
  # optimizer='adam',
  # optimizer=Adam(lr=0.01),
  optimizer=SGD(lr=0.08, momentum=0.9),
  metrics=['mse'],
)

r = model.fit(
  x=[df_train.userId.values, df_train.movie_idx.values],
  y=df_train.rating.values - mu,
  epochs=epochs,
  batch_size=128,
  validation_data=(
    [df_test.userId.values, df_test.movie_idx.values],
    df_test.rating.values - mu
  )
)

# plot losses
plt.plot(r.history['loss'], label="train loss")
plt.plot(r.history['val_loss'], label="test loss")
plt.legend()
plt.show()

# plot mse
plt.plot(r.history['mean_squared_error'], label="train mse")
plt.plot(r.history['val_mean_squared_error'], label="test mse")
plt.legend()
plt.show()

max_stars_repo_name: JouniVatanen/NLP-and-Deep-Learning


text:
import numpy as np
import tensorflow as tf

DEFAULT_FILTERS_COUNT = 1
DEFAULT_KERNEL_SIZE = 5
DEFAULT_POOL_SIZE = 2
DEFAULT_STRIDES = 2
DEFAULT_FC_SIZE = 128

def net(in_data, layers, fc_dropout, conv_dropout):
    """
    Build neural net of architecture described in layers list. Layers types
    include convolutional, pool and feed-forward (deep) layers, eg.:

    layers = [
        {"type": "conv", "filters_count": 32, "kernel_size": 5},
        {"type": "pool", "strides": 4},
        {"type": "conv", "filters_count": 128, "kernel_size": 5},
        {"type": "deep", "out_size": 100},
        {"type": "deep", "out_size": 7}
    ]

    :param x:       input placeholder
    :param layers:  list of dicts describing layers
    :return:        output variable of a net
    """
    x = in_data

    # build convolutional and pool layers
    for i, layer in enumerate(layers):
        layer_type = layer.get("type")

        if layer_type is "conv":
            filters_count = layer.get("filters_count") or DEFAULT_FILTERS_COUNT
            kernel_size = layer.get("kernel_size") or DEFAULT_KERNEL_SIZE

            x = tf.layers.conv2d(
                inputs=x,
                filters=filters_count,
                kernel_size=[kernel_size, kernel_size],
                padding="same",
                activation=tf.nn.relu
            )
        elif layer_type is "pool":
            poll_size = layer.get("pool_size") or DEFAULT_POOL_SIZE
            strides = layer.get("strides") or DEFAULT_STRIDES

            x = tf.layers.max_pooling2d(
                inputs=x,
                pool_size=poll_size,
                strides=strides
            )
            print(x.shape)
            x = tf.nn.dropout(
                x, 1 - conv_dropout,
                noise_shape=[tf.shape(x)[0], 1, 1, tf.shape(x)[-1]]
            )
            print(x.shape)
        else:
            break

    x_shape = x.get_shape().as_list()

    # unroll input to fit deep layer
    if 3 < len(x_shape):
        # multiply all dimensions except for the first that represents the
        # number of input vectors
        in_size = np.product(x_shape[1:])
        x = tf.reshape(x, [-1, in_size])

    for j, layer in enumerate(layers[i:]):
        x = build_fc_layer("fc_{}".format(j), x, in_size, layer["out_size"],
                           fc_dropout, j == len(layers[i:]) - 1)
        in_size = layer.get("out_size", DEFAULT_FC_SIZE)

    return x

def build_fc_layer(name, in_data, in_size, out_size, dropout, is_out):
    """ Build a fc layer on the top of previous output. """
    with tf.variable_scope(name):
        W = tf.get_variable("W", shape=[in_size, out_size],
                            initializer=tf.contrib.layers.xavier_initializer())
        b = tf.get_variable("b", shape=[out_size, ],
                            initializer=tf.constant_initializer())
    if not is_out:
        return tf.nn.dropout(tf.nn.relu(tf.matmul(in_data, W) + b), 1 - dropout)
    else:
        return tf.matmul(in_data, W) + b

max_stars_repo_name: AffectiveComputing/temp

max_issues_repo_name: AffectiveComputing/temp

max_forks_repo_name: AffectiveComputing/emu


text:
import random
from collections import deque

import numpy as np

class BaseBuffer:
    """
    Base class for replay buffers.
    """

    def __init__(self, size, initial_size=None, batch_size=32):
        """
        Initialize replay buffer.
        Args:
            size: Buffer maximum size.
            initial_size: Buffer initial size to be filled before training starts.
                To be used by caller.
            batch_size: Size of the batch that should be used in get_sample() implementation.
        """
        assert (
            initial_size is None or initial_size > 0
        ), f'Buffer initial size should be > 0, got {initial_size}'
        assert size > 0, f'Buffer size should be > 0,  got {size}'
        assert batch_size > 0, f'Buffer batch size should be > 0, got {batch_size}'
        assert (
            batch_size <= size
        ), f'Buffer batch size `{batch_size}` should be <= size `{size}`'
        if initial_size:
            assert size >= initial_size, 'Buffer initial size exceeds max size'
        self.size = size
        self.initial_size = initial_size or size
        self.batch_size = batch_size
        self.current_size = 0

    def append(self, *args):
        """
        Add experience to buffer.
        Args:
            *args: Items to store, types are implementation specific.

        """
        raise NotImplementedError(
            f'append() should be implemented by {self.__class__.__name__} subclasses'
        )

    def get_sample(self):
        """
        Sample from stored experience.

        Returns:
            Sample as numpy array.
        """
        raise NotImplementedError(
            f'get_sample() should be implemented by {self.__class__.__name__} subclasses'
        )

class ReplayBuffer1(BaseBuffer):
    """
    deque-based replay buffer that holds state transitions
    """

    def __init__(self, size, **kwargs):
        """
        Initialize replay buffer.
        Args:
            size: Buffer maximum size.
            **kwargs: kwargs passed to BaseBuffer.
        """
        super(ReplayBuffer1, self).__init__(size, **kwargs)
        self.main_buffer = deque(maxlen=size)
        self.temp_buffer = []

    def append(self, *args):
        """
        Append experience and auto-allocate to temp buffer / main buffer(self)
        Args:
            *args: Items to store

        Returns:
            None
        """
        self.main_buffer.append(args)
        self.current_size = len(self.main_buffer)

    def get_sample(self):
        """
        Sample from stored experience.

        Returns:
            Same number of args passed to append, having self.batch_size as
            first shape.
        """
        memories = random.sample(self.main_buffer, self.batch_size)
        if self.batch_size > 1:
            return [np.array(item) for item in zip(*memories)]
        return memories[0]

class ReplayBuffer2(BaseBuffer):
    """
    numpy-based replay buffer added for compatibility with tensorflow shortcomings.
    """

    def __init__(self, size, slots, **kwargs):
        """
        Initialize replay buffer.

        Args:
            size: Buffer maximum size.
            slots: Number of args that will be passed to self.append()
            **kwargs: kwargs passed to BaseBuffer.
        """
        super(ReplayBuffer2, self).__init__(size, **kwargs)
        self.slots = [np.array([])] * slots
        self.current_size = 0

    def append(self, *args):
        """
        Add experience to buffer.
        Args:
            *args: Items to store.

        Returns:
            None
        """
        for i, arg in enumerate(args):
            if not isinstance(arg, np.ndarray):
                arg = np.array([arg])
            if not self.slots[i].shape[0]:
                self.slots[i] = np.zeros((self.size, *arg.shape), arg.dtype)
            self.slots[i][self.current_size % self.size] = arg.copy()
        if self.current_size < self.size:
            self.current_size += 1

    def get_sample(self):
        """
        Sample from stored experience.

        Returns:
            Same number of args passed to append, having self.batch_size as
            first shape.
        """
        indices = np.random.randint(
            0, min(self.current_size, self.size), self.batch_size
        )
        return [slot[indices] for slot in self.slots]


text:
def agent2Dposition(x_info,agent_position):
    '''
    This function is to prepare the dataset that gets used to plot the animated visualisation of the results.
    This function puts all agents positions for each iteration and indicates the family that agent belongs to.
    Parameters:
    ---------------
    x_info: This is the dataframe which contains each agents's characteristics including details about the family they belong to.
    agent_position: This is the dataframe that contains all the updated positions for each agent at different iterations.

    Returns:
    ---------------
    result : This is the dataframe with all individuals position details at different iterations including details about the family they belong to.
    '''
    import numpy as np
    import pandas as pd
    fam_id = x_info[['Family']]
    fam_id = fam_id.reset_index()

    df = agent_position.reset_index()
    iterations = np.unique(df['iter'])

    ind_pos = pd.DataFrame()
    for i in iterations:
        df2 = df[df['iter']==i]
        df2 = df2.drop(['iter'],axis=1)
        df2 = df2.set_index('position')
        df2 = np.transpose(df2)
        a=i.split('_')
        df2['iter'] = int(a[1])
        ind_pos = ind_pos.append(df2)
        df2.head()
    ind_pos = ind_pos.reset_index()
    ind_pos = ind_pos.sort_values(by=['iter'], ascending=True)
    result = pd.merge(ind_pos,fam_id, on='index')
    result = result.sort_values(by=['iter'], ascending=True)
    return result

def position2Dplot(result,lb,ub,render):
    '''
    This function uses the ploty package to plot the results.
    Parameters:
    -------------
    result:  This is the dataframe with all individuals position details at different iterations including details about the
                family they belong to. Prepared using agent2Dposition((x_info,agent_position).

    Returns:
    ------------
    fig1 : This first figure shows the position of the agent within the whole community
    fig2 : This second figure shows the position of the agent within their respective families.
    '''
    import plotly.express as px
    import plotly.io as pio
    import plotly.graph_objs as go
    from plotly.subplots import make_subplots
    import numpy as np
    import pandas as pd
    pio.renderers.default = render
    fig = px.scatter(result, x="IndPosition0",y = "IndPosition1",
           animation_frame="iter",
           animation_group="index",color="Family",hover_name="index",title='
Individual Position
',log_x=False,
           range_x=[(lb-5),(ub+5)],range_y=[(lb-5),(ub+5)])

    fig2 = px.scatter(result, x="IndPosition0",y = "IndPosition1",
           animation_frame="iter",
           animation_group="index",color="Family",hover_name="index", facet_col='Family',title='
Family Position Updates
',log_x=False,
           range_x=[(lb-5),(ub+5)],range_y=[(lb-5),(ub+5)])

    return fig.show(),fig2.show()

def animated2Dvisualisation(x_info,x_pos,lb,ub,render='notebook_connected'):
    '''
    This function plots the animated plots showing how each agent moved around in the community and within their family in search of the optimal solution.
    Parameters:
    ---------------
    x_info: This is the dataframe which contains each agents's characteristics including details about the family they belong to.
    agent_position: This is the dataframe that contains all the updated positions for each agent at different iterations.
    lb : Lower Bounds
    ub : Upper Bounds
    render      : This is the render options as defined by help(plotly.io.renderers)

    Returns:
    ------------
    fig1 : This first figure shows the position of the agent within the whole community
    fig2 : This second figure shows the position of the agent within their respective families.
    '''
    import plotly.express as px
    import plotly.io as pio
    import plotly.graph_objs as go
    from plotly.subplots import make_subplots
    import numpy as np
    import pandas as pd
    pio.renderers.default = render
    result = agent2Dposition(x_info,x_pos)
    fig1,fig2 = position2Dplot(result,lb,ub,render)
    return fig1,fig2

def convergenceRatePlot(ind_payoffs,objFunction,render='notebook_connected'):
    '''
    This function plots the convergence curve of the obejective function being investigated.
    Parameters:
    -------------
    ind_payoffs : These are the objective values of all agents in the community.
    objFunction : This is the objective function being solved for.
    render      : This is the render options as defined by help(plotly.io.renderers)

    Returns:
    ------------
    fig : This is the convergence rate plot.
    '''
    import plotly.express as px
    import plotly.io as pio
    import plotly.graph_objs as go
    from plotly.subplots import make_subplots
    import numpy as np
    import pandas as pd
    pio.renderers.default = render
    data = ind_payoffs.reset_index(drop=True)
    data['Objective Value'] = np.mean(data,axis=1)
    data['Iteration'] = data.index
    fig = px.line(data, x='Iteration',y = "Objective Value",
              title=objFunction.__name__ + '
Convergence Rate
',
        line_shape="spline",render_mode="svg")

    return fig.show()


text:
import torch
import numpy as np

def calc_floor_ceil_delta(x):
    '''
    x_fl = torch.floor(x + 1e-8)
    x_ce = torch.ceil(x - 1e-8)
    x_ce_fake = torch.floor(x) + 1

    dx_ce = x - x_fl
    dx_fl = x_ce_fake - x
    return [x_fl.long(), dx_fl], [x_ce.long(), dx_ce]
    '''

    x_fl = torch.floor(x)
    x_ce = x_fl + 1

    dx_ce = x - x_fl
    dx_fl = x - x_ce

    return [x_fl.long(), dx_fl], [x_ce.long(), dx_ce]

def create_update(x, y, t, dt, p, vol_size, device="cpu"):
    '''
    assert (x>=0).byte().all()
    assert (x<vol_size[2]).byte().all()
    assert (y>=0).byte().all()
    assert (y<vol_size[1]).byte().all()
    assert (t>=0).byte().all()
    '''
    #assert (t<vol_size[0] // 2).byte().all()

    ''''
    if not (t < vol_size[0] // 2).byte().all():
        print(t[t >= vol_size[0] // 2])
        print(vol_size)
        raise AssertionError()
    '''

    vol_mul = torch.where(p < 0,
                          torch.ones(p.shape, dtype=torch.long).to(device) * vol_size[0] // 2,
                          torch.zeros(p.shape, dtype=torch.long).to(device))

    inds = (vol_size[1]*vol_size[2]) * (t + vol_mul)\
         + (vol_size[2])*y\
         + x

    vals = dt

    return inds, vals

def gen_discretized_event_volume(x, y, t, p, vol_size, device="cpu"):
    # volume is [t, x, y]
    # events are Nx4
    npts = x.shape[0]
    volume = torch.zeros(vol_size)

    t_min = t.min()
    t_max = t.max()
    t_scaled = torch.clip((t-t_min) / (t_max - t_min) - 1e-5, 0) * (vol_size[0] // 2-1)

    ts_fl, ts_ce = calc_floor_ceil_delta(t_scaled.squeeze())

    inds_fl, vals_fl = create_update(x, y,
                                     ts_fl[0], ts_fl[1],
                                     p,
                                     vol_size,
                                     device=device)
    inds_ce, vals_ce = create_update(x, y,
                                     ts_ce[0], ts_ce[1],
                                     p,
                                     vol_size,
                                     device=device)

    volume.view(-1).put_(torch.stack((inds_fl, inds_ce)),
                        torch.stack((vals_fl, vals_ce)), accumulate=True)
    volume.view(-1).put_(inds_ce, vals_ce, accumulate=True)
    return volume

def normalize_event_volume(event_volume):
    event_volume_flat = event_volume.view(-1)
    nonzero = torch.nonzero(event_volume_flat)
    nonzero_values = event_volume_flat[nonzero]
    if nonzero_values.shape[0]:
        lower = torch.kthvalue(nonzero_values,
                               max(int(0.02 * nonzero_values.shape[0]), 1),
                               dim=0)[0][0]
        upper = torch.kthvalue(nonzero_values,
                               max(int(0.98 * nonzero_values.shape[0]), 1),
                               dim=0)[0][0]
        max_val = max(abs(lower), upper)
        event_volume = torch.clamp(event_volume, -max_val, max_val)
        event_volume /= max_val
    return event_volume

max_stars_repo_name: ZiyunClaudeWang/AtlasNet

max_issues_repo_name: ZiyunClaudeWang/AtlasNet

max_forks_repo_name: ZiyunClaudeWang/AtlasNet