text:
import time
import torch
import torch.nn as nn
import numpy as np
import json
from torch.utils.data import DataLoader
from torch.utils.tensorboard import SummaryWriter
from pathlib import Path
from tqdm import tqdm
from tabulate import tabulate
from statistics import mean

import utils.metrics as metrics
from utils.utils import *

class Evaluator:
    def __init__(
        self, cnn: nn.Module, data_loader: DataLoader, device: torch.device, name: str,
    ):
        self.cnn = cnn
        self.data_loader = data_loader
        self.device = device
        self.name = name
        self.predictions = {}

    def predict(self):

        # Turn on evaluation for network
        self.cnn.model.eval()
        results = {"labels": [], "logits": [], "predictions": []}

        with torch.no_grad():
            for i, (spatial_data, temporal_data, label, metadata) in enumerate(
                tqdm(self.data_loader, desc="Prediction", leave=False, unit="sample")
            ):

                spatial_data = spatial_data.to(self.device)
                temporal_data = temporal_data.to(self.device)
                label = label.to(self.device)

                ape_id = metadata["ape_id"][0]
                start_frame = metadata["start_frame"][0]
                video = metadata["video"][0]

                logits = self.cnn.model(spatial_data, temporal_data)

                # Accumulate predictions against ground truth labels
                prediction = logits.argmax().item()

                # Logits are required for Top-k accuracy
                logits = (logits.detach().cpu()).numpy()

                # Collect reulsts for metrics
                results["labels"].append(label.item())
                results["logits"].append(logits[0].tolist())
                results["predictions"].append(prediction)

                # Insert results to dictionary
                if video not in self.predictions.keys():
                    self.predictions[video] = []

                self.predictions[video].append(
                    {
                        "ape_id": ape_id.item(),
                        "label": label.item(),
                        "prediction": prediction,
                        "start_frame": start_frame.item(),
                    }
                )

        # Get accuracy by checking for correct predictions across all predictions
        top1, top3 = metrics.compute_topk_accuracy(
            torch.LongTensor(results["logits"]), torch.LongTensor(results["labels"]), topk=(1, 3)
        )

        # Get per class accuracies and sort by label value (0...9)
        class_accuracy = metrics.compute_class_accuracy(results["labels"], results["predictions"])
        class_accuracy_average = mean(class_accuracy.values())

        print("==> Per Class Results")
        per_class_results = [
            [
                "Class",
                "camera_interaction",
                "climbing_down",
                "climbing_up",
                "hanging",
                "running",
                "sitting",
                "sitting_on_back",
                "standing",
                "walking",
            ],
            [
                "Accuracy",
                f"{class_accuracy[0]:2.2f}",
                f"{class_accuracy[1]:2.2f}",
                f"{class_accuracy[2]:2.2f}",
                f"{class_accuracy[3]:2.2f}",
                f"{class_accuracy[4]:2.2f}",
                f"{class_accuracy[5]:2.2f}",
                f"{class_accuracy[6]:2.2f}",
                f"{class_accuracy[7]:2.2f}",
                f"{class_accuracy[8]:2.2f}",
            ],
        ]

        print(tabulate(per_class_results, tablefmt="fancy_grid"))

        print("==> Overall Results")
        test_results = [
            ["Average Class Accuracy:", f"{class_accuracy_average:2f}"],
            ["Top1 Accuracy:", f"{top1.item():.2f}"],
            ["Top3 Accuracy:", f"{top3.item():.2f}"],
        ]

        print(tabulate(test_results, tablefmt="fancy_grid"))

        return self.predictions

max_stars_repo_name: fznsakib/great-ape-behaviour-detector

max_stars_repo_licenses: FTL CECILL-B

max_issues_repo_name: fznsakib/great-ape-behaviour-detector

max_issues_repo_licenses: FTL CECILL-B

max_forks_repo_name: fznsakib/great-ape-behaviour-detector

max_forks_repo_licenses: FTL CECILL-B


text:
'''
Apply FiRE to the pbmc  68k dataset and visualization

@author:chenxofhit@gmail.com
@since: 2019-01-10
'''
import warnings
warnings.filterwarnings("ignore", message="numpy.dtype size changed")
warnings.filterwarnings("ignore", message="numpy.ufunc size changed")
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
plt.rcParams['svg.fonttype'] = 'none'

# set font
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = 'Roboto'

# set the style of the axes and the text color
plt.rcParams['axes.edgecolor']='#333F4B'
plt.rcParams['axes.linewidth']=0.8
plt.rcParams['xtick.color']='#333F4B'
plt.rcParams['ytick.color']='#333F4B'
plt.rcParams['text.color']='#333F4B'
plt.rcParams['legend.fontsize']='small'

import time
import numpy as np
from scipy import io

import sys
sys.path.append('utils')
import preprocess as pp
import dao as dao
import pickle

method = "fire"

filename = "./data/dat.cell_expression.mm"
X = io.mmread(filename)

X_array = X.toarray()

genes = np.arange(1, X_array.shape[1]+1)

datasave = "output_68k"
dao.create_folder(datasave)
preprocessedData, selGenes = pp.ranger_preprocess(X_array, genes, optionToSave=False, dataSave=datasave)
pickle.dump(preprocessedData, open("./output_68k/preprocessedData.p", "wb"))
pickle.dump(selGenes, open("./output_68k/selGenes.p", "wb"))

import fire
predictions, scores, duration = fire.fire(preprocessedData)
pickle.dump(predictions, open("./output_68k/predictions_{}.p".format(method), "wb"))
pickle.dump(scores, open("./output_68k/scores_{}.p".format(method), "wb"))

from MulticoreTSNE import MulticoreTSNE as TSNE
t0 = time.time()
tsne = TSNE(n_jobs=16, random_state=1)
X_embedded = tsne.fit_transform(preprocessedData)
t1 = time.time()
durationTsne = round(t1 - t0, ndigits=4)
print("Total running tsne time is :" + str(durationTsne) + " s")

pickle.dump(X_embedded, open("./output_68k/X_embedded.p", "wb"))

import pickle
X_embedded = pickle.load(open("./output_68k/X_embedded.p", "rb"))
predictions = pickle.load(open("./output_68k/predictions_{}.p".format(method), "rb"))

import pandas as pd

annotation_name = "./data/dat.cell_annotation.csv"
annotation = pd.read_csv(annotation_name, sep=",", header=0)
pickle.dump(annotation, open("./output_68k/annotation.p", "wb"))

fig, ax = plt.subplots()
plt.scatter(X_embedded[:, 0], X_embedded[:, 1], s=5, c='lightgray', edgecolors=None, marker='.', linewidths=0)
rareind = np.where(predictions == 1)[0]
plt.scatter(X_embedded[rareind, 0], X_embedded[rareind, 1], s=5, c='blue', edgecolors=None, marker='.', linewidths=0, label='rare')

# change the style of the axis spines
#ax.spines['top'].set_color('none')
#ax.spines['right'].set_color('none')
# ax.spines['left'].set_smart_bounds(True)
# ax.spines['bottom'].set_smart_bounds(True)
#
# # set the spines position
# ax.spines['bottom'].set_position(('axes', -0.04))
# ax.spines['left'].set_position(('axes', 0.015))

#remove axis
#plt.axis('off')
plt.xlabel("tSNE-1")
plt.ylabel("tSNE-2")
#remove white margin of the plots
plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0, hspace=0.2)

plt.savefig("./output_68k/plot_tsne_{}.svg".format(method), format='svg', dpi=100)
plt.savefig("./output_68k/plot_tsne_{}.png".format(method), dpi=100)

plt.show()

ext: py

max_stars_repo_name: CSUBioGroup/DoRC

max_issues_repo_name: CSUBioGroup/DoRC

max_forks_repo_name: CSUBioGroup/DoRC

reason: import numpy,from scipy


text:
"""Score proset classifier trained on the 'continuous XOR' problem with 4 features.

Copyright by Nikolaus Ruf
Released under the MIT license - see LICENSE file for details
"""

import gzip
import os
import pickle

import numpy as np
from sklearn.metrics import classification_report, log_loss, roc_auc_score

import proset.utility as utility

print("* Apply user settings")
input_path = "scripts/results"
output_path = "scripts/reports"
input_file = "xor_4f_2d_95_model.gz"
export_file = input_file.replace(".gz", "_export.xlsx")
model_name = input_file.replace(".gz", "")

print("* Load model fit results")
with gzip.open(os.path.join(input_path, input_file), mode="rb") as file:
    result = pickle.load(file)

print("* Generate prototype report")
scale = np.sqrt(result["model"]["transform"].var_)
offset = result["model"]["transform"].mean_
export = result["model"]["model"].export(
    feature_names=result["data"]["feature_names"],
    scale=scale,
    offset=offset
)
utility.write_report(file_path=os.path.join(output_path, export_file), report=export)

print("* Show results")
train_features = result["model"]["transform"].transform(result["data"]["X_train"])
train_labels = result["data"]["y_train"]
test_features = result["model"]["transform"].transform(result["data"]["X_test"])
test_labels = result["data"]["y_test"]
prediction = result["model"]["model"].predict(test_features)
probabilities = result["model"]["model"].predict_proba(test_features)
active_features = result["model"]["model"].set_manager_.get_active_features()
misclassified = prediction != test_labels
print("- Hyperparameter selection")
utility.print_hyperparameter_report(result)
print("-  Final model")
print("log-loss = {:.2f}".format(log_loss(y_true=test_labels, y_pred=probabilities)))
print("roc-auc  = {:.2f}".format(roc_auc_score(y_true=test_labels, y_score=probabilities[:, 1])))
print("number of active features = {}".format(active_features.shape[0]))
print("number of prototypes = {}\n".format(result["model"]["model"].set_manager_.get_num_prototypes()))
print("- Selected features and weights")
utility.print_feature_report(model=result["model"]["model"], feature_names=result["data"]["feature_names"])
print("- Classification report")
print(classification_report(y_true=test_labels, y_pred=prediction))
utility.plot_select_results(result=result, model_name=model_name)
plotter = utility.ClassifierPlots(
    model=result["model"]["model"],
    model_name=model_name,
    feature_names=result["data"]["feature_names"],
    scale=scale,
    offset=offset
)
ix = np.prod(result["data"]["X_train"][:, 2:], axis=1) >= 0
# select test samples which have identical class based on the first two features
x_range, y_range = plotter.plot_surface(
    features=train_features[ix, :],
    target=train_labels[ix],
    baseline=np.ones((1, 4)) * 0.5,  # fix remaining features to positive sign
    plot_index=np.array([0, 1]),
    comment="training samples",
    use_proba=True
)
ix = np.prod(result["data"]["X_test"][:, 2:], axis=1) >= 0
# noinspection PyUnresolvedReferences
plotter.plot_surface(
    features=test_features[ix, :],
    target=test_labels[ix],
    baseline=np.ones((1, 4)) * 0.5,  # fix remaining features to positive sign
    plot_index=np.array([0, 1]),
    comment="test samples",
    highlight=misclassified[ix],
    highlight_name="misclassified",
    x_range=x_range,
    y_range=y_range
)
x_range, y_range = plotter.plot_batch_map(
    batch=1,
    features=train_features,
    target=train_labels,
    comment="training samples",
    show_index=False
)
plotter.plot_batch_map(
    batch=1,
    features=test_features,
    target=test_labels,
    comment="test samples",
    highlight=misclassified,
    highlight_name="misclassified",
    x_range=x_range,
    y_range=y_range,
    show_index=False
)
plotter.plot_batch_map(batch=1, x_range=x_range, y_range=y_range, show_index=False)
# feature plot is pointless as no pair of features carries enough information to determine the class

print("* Done")

reason: import numpy


text:
import os
import hashlib
import tarfile
import zipfile
import requests

import torch
import numpy as np
import pandas as pd
from torch import nn
from d2l import torch as d2l

DATA_HUB = dict()
DATA_URL = "http://d2l-data.s3-accelerate.amazonaws.com/"

DATA_HUB = {
    "kaggle_house_train": (
        DATA_URL + "kaggle_house_pred_train.csv",
        "585e9cc93e70b39160e7921475f9bcd7d31219ce",
    ),
    "kaggle_house_test": (
        DATA_URL + "kaggle_house_pred_test.csv",
        "fa19780a7b011d9b009e8bff8e99922a8ee2eb90",
    ),
}

def download(name, cache_dir=os.path.join("..", "data")):
    """Download a file inserted into DATA_HUB, return the local filename."""
    assert name in DATA_HUB, f"{name} does not exist in {DATA_HUB}."
    url, sha1_hash = DATA_HUB[name]
    os.makedirs(cache_dir, exist_ok=True)
    file_name = os.path.join(cache_dir, url.split("/")[-1])
    if os.path.exists(file_name):
        sha1 = hashlib.sha1()
        with open(file_name, "rb", encoding="UTF-8") as input_file:
            while True:
                data = input_file.read(1048576)
                if not data:
                    break
                sha1.update(data)
        if sha1.hexdigest() == sha1_hash:
            return file_name  # Hit cache
    print(f"Downloading {file_name} from {url}...")
    request_get = requests.get(url, stream=True, verify=True)
    with open(file_name, "wb", encoding="UTF-8") as output_file:
        output_file.write(request_get.content)
    return file_name

def download_and_extract_zip_file(name, folder=None):
    """Download and extract a zip/tar file."""
    file_name = download(name)
    base_dir = os.path.dirname(file_name)
    data_dir, ext = os.path.splitext(file_name)
    if ext == ".zip":
        fp = zipfile.ZipFile(file_name, "r")
    elif ext in (".tar", ".gz"):
        fp = tarfile.open(file_name, "r")
    else:
        assert False, "Only zip/tar files can be extracted."
    fp.extractall(base_dir)
    return os.path.join(base_dir, folder) if folder else data_dir

def download_all(data_hub=DATA_HUB):
    """Download all files in the DATA_HUB."""
    for name in data_hub:
        download(name)

def get_net():
    net = nn.Sequential(nn.Linear(in_features, 1))
    return net

def log_rmse(net, features, labels):
    """
    This leads to the following root-mean-squared-error between
    the logarithm of the predicted price and the logarithm of the label price
    """
    # To further stabilize the value when the logarithm is taken, set the
    # value less than 1 as 1
    clipped_preds = torch.clamp(net(features), 1, float("inf"))
    rmse = torch.sqrt(nn.MSELoss(torch.log(clipped_preds), torch.log(labels)))
    return rmse.item()

def train(
    net,
    train_features,
    train_labels,
    test_features,
    test_labels,
    num_epochs,
    learning_rate,
    weight_decay,
    batch_size,
):
    """
    Defining the Training Loop
    """
    train_log_rmse, test_log_rmse = [], []
    train_iter = d2l.load_array((train_features, train_labels), batch_size)
    # The Adam optimization algorithm is used here
    optimizer = torch.optim.Adam(
        net.parameters(), lr=learning_rate, weight_decay=weight_decay
    )
    for _ in range(num_epochs):
        for X, y in train_iter:
            optimizer.zero_grad()
            train_loss = nn.MSELoss(net(X), y)
            train_loss.backward()
            optimizer.step()
        train_log_rmse.append(log_rmse(net, train_features, train_labels))
        if test_labels is not None:
            test_log_rmse.append(log_rmse(net, test_features, test_labels))
    return train_log_rmse, test_log_rmse

def get_data_in_k_fold_cross_validation_procedure(k, i, X, y):
    """
    returns the  i-th  fold of the data in a  K-fold cross-validation procedure.
    It proceeds by slicing out the  i-th  segment as validation data and returning the rest as training data.
    """
    assert k > 1
    fold_size = X.shape[0] // k
    X_train, y_train = None, None
    for j in range(k):
        idx = slice(j * fold_size, (j + 1) * fold_size)
        X_part, y_part = X[idx, :], y[idx]
        if j == i:
            X_valid, y_valid = X_part, y_part
        elif X_train is None:
            X_train, y_train = X_part, y_part
        else:
            X_train = torch.cat([X_train, X_part], 0)
            y_train = torch.cat([y_train, y_part], 0)
    return X_train, y_train, X_valid, y_valid

def k_fold(
    k, X_train, y_train, num_epochs, learning_rate, weight_decay, batch_size
):
    """
    The training and verification error averages are returned when we train  K  times in the  K -fold cross-validation.
    """
    train_l_sum, valid_l_sum = 0, 0
    for i in range(k):
        data = get_data_in_k_fold_cross_validation_procedure(
            k, i, X_train, y_train
        )
        net = get_net()
        train_ls, valid_ls = train(
            net, *data, num_epochs, learning_rate, weight_decay, batch_size
        )
        train_l_sum += train_ls[-1]
        valid_l_sum += valid_ls[-1]
        if i == 0:
            d2l.plot(
                list(range(1, num_epochs + 1)),
                [train_ls, valid_ls],
                xlabel="epoch",
                ylabel="rmse",
                xlim=[1, num_epochs],
                legend=["train", "valid"],
                yscale="log",
            )
        print(
            f"fold {i + 1}, train log rmse {float(train_ls[-1]):f}, "
            f"valid log rmse {float(valid_ls[-1]):f}"
        )
    return train_l_sum / k, valid_l_sum / k

def train_and_pred(
    train_features,
    test_feature,
    train_labels,
    test_data,
    num_epochs,
    lr,
    weight_decay,
    batch_size,
):
    """
    train and predict on Kaggle
    """
    net = get_net()
    train_ls, _ = train(
        net,
        train_features,
        train_labels,
        None,
        None,
        num_epochs,
        lr,
        weight_decay,
        batch_size,
    )
    d2l.plot(
        np.arange(1, num_epochs + 1),
        [train_ls],
        xlabel="epoch",
        ylabel="log rmse",
        xlim=[1, num_epochs],
        yscale="log",
    )
    print(f"train log rmse {float(train_ls[-1]):f}")
    # Apply the network to the test set
    preds = net(test_features).detach().numpy()
    # Reformat it to export to Kaggle
    test_data["SalePrice"] = pd.Series(preds.reshape(1, -1)[0])
    submission = pd.concat([test_data["Id"], test_data["SalePrice"]], axis=1)
    submission.to_csv("submission.csv", index=False)

if __name__ == "__main__":
    train_data = pd.read_csv(download("kaggle_house_train"))
    test_data = pd.read_csv(download("kaggle_house_test"))

    print(train_data.shape)
    print(test_data.shape)
    print(train_data.iloc[0:4, [0, 1, 2, 3, -3, -2, -1]])
    all_features = pd.concat((train_data.iloc[:, 1:-1], test_data.iloc[:, 1:]))
    # If test data were inaccessible, mean and standard deviation could be
    # calculated from training data
    numeric_features = all_features.dtypes[
        all_features.dtypes != "object"
    ].index
    all_features[numeric_features] = all_features[numeric_features].apply(
        lambda x: (x - x.mean()) / (x.std())
    )
    # After standardizing the data all means vanish, hence we can set missing
    # values to 0
    all_features[numeric_features] = all_features[numeric_features].fillna(0)
    # `Dummy_na=True` considers "na" (missing value) as a valid feature value, and
    # creates an indicator feature for it
    all_features = pd.get_dummies(all_features, dummy_na=True)
    all_features.shape
    n_train = train_data.shape[0]
    train_features = torch.tensor(
        all_features[:n_train].values, dtype=torch.float32
    )
    test_features = torch.tensor(
        all_features[n_train:].values, dtype=torch.float32
    )
    train_labels = torch.tensor(
        train_data.SalePrice.values.reshape(-1, 1), dtype=torch.float32
    )
    in_features = train_features.shape[1]

    k, num_epochs, lr, weight_decay, batch_size = 5, 100, 5, 0, 64
    train_l, valid_l = k_fold(
        k,
        train_features,
        train_labels,
        num_epochs,
        lr,
        weight_decay,
        batch_size,
    )
    print(
        f"{k}-fold validation: avg train log rmse: {float(train_l):f}, "
        f"avg valid log rmse: {float(valid_l):f}"
    )

    train_and_pred(
        train_features,
        test_features,
        train_labels,
        test_data,
        num_epochs,
        lr,
        weight_decay,
        batch_size,
    )


text:
import sys
import os
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + "/../")
try:
    from utils.priority_queue import PriorityQueue
except:
    raise

import numpy as np

def test_priority_queue(arr):
    output = sorted(arr)
    pq = PriorityQueue()

    for i in arr:
        pq.put(i, i)

    check = []
    for j in pq.len():
        check.append(pq.get())

    print(check, output)
    assert(check == output)

arr = np.random.shuffle(range(50))

max_stars_repo_name: Liam-Stow/novaRoverDemos

max_issues_repo_name: Liam-Stow/novaRoverDemos

max_forks_repo_name: Liam-Stow/novaRoverDemos


text:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize

from eugene import simulate_outbreak

from grid_parallel import params

np.random.seed(2019)

fig, ax = plt.subplots(figsize=(4, 3))

R0_grid = np.linspace(0.7, 2.5, 10)

cmap = lambda x: plt.cm.viridis((x - R0_grid.min())/R0_grid.ptp())

for j in range(R0_grid.shape[0]):
    for i in range(20):
        parameters = dict(
            R0=R0_grid[j],
            k=1,
            n=1,
            D=10,
            gamma_shape=2,
            max_time=90,
            days_elapsed_max=52+7,
            max_cases=1e4
        )

        times, cumulative_incidence = simulate_outbreak(**parameters)

        days_elapsed_min = params['days_elapsed_min']
        days_elapsed = parameters['days_elapsed_max']
        min_number_cases = params['min_number_cases']
        max_number_cases = params['max_number_cases']

        delta_t = (np.array(days_elapsed_min) -
                   max(days_elapsed_min))
        cases_at_measurement_times = np.interp(days_elapsed +
                                               delta_t,
                                               times, cumulative_incidence)

        accept = ((np.asarray(min_number_cases) <
                   cases_at_measurement_times) &
                  (cases_at_measurement_times <
                   np.asarray(max_number_cases))).all()
        cax = ax.semilogy(times, cumulative_incidence, '.-',
                          color=cmap(R0_grid[j]),
                          alpha=1.0 if accept else 0.25,
                          zorder=1 if not accept else 5)

plot_kwargs = dict(fmt='s', color='k', zorder=10, ecolor='k')
ax.errorbar(48, 2890, xerr=14, yerr=[[2700], [2700]], **plot_kwargs)
ax.errorbar(52, 5000, xerr=14, yerr=[[4000], [4700]], **plot_kwargs)

norm = Normalize(vmin=R0_grid.min(), vmax=R0_grid.max())
cbar = plt.colorbar(mappable=ScalarMappable(norm=norm, cmap=plt.cm.viridis),
                    label='$\mathcal{R}_0$')

ax.set_xticks(np.arange(0, 60, 1), minor=True)

for sp in ['right', 'top']:
    ax.spines[sp].set_visible(False)

ax.set_xlabel('Time [days]')
ax.set_ylabel('Cumulative Incidence')
fig.tight_layout()
fig.savefig('plots/trajectories.pdf', bbox_inches='tight')
plt.show()


text:
"""Tests for modularity code.
"""

from __future__ import print_function

#-----------------------------------------------------------------------------
# Library imports
#-----------------------------------------------------------------------------

# Stdlib
import os
import copy
import math
import time

# Third party
import networkx as nx
import numpy as np
import nose.tools as nt
import numpy.testing as npt

# Our own
from brainx import modularity as mod
from brainx import util

# While debugging the library, reload everything
#map(reload,[mod,util])

#-----------------------------------------------------------------------------
# Local utility functions
#-----------------------------------------------------------------------------

def betweenness_to_modularity(g,ppart):
    """Function to convert between betweenness fractions and modularity
    Parameters:
    ----------
    g = graph object
    ppart = perfect partition

    Returns:
    --------
    mod = best modularity associated with this graph object
    """

    graph_partition = mod.GraphPartition(g,ppart)
    return graph_partition.modularity()

#-----------------------------------------------------------------------------
# Tests
#-----------------------------------------------------------------------------
def test_graphpartition():
    """ test GraphPartition correctly handles graph whose
    nodes are strings"""
    graph = nx.Graph()
    graph.add_edge('a','b')
    graph.add_edge('c','d')
    index = {0:set([0,1]), 1:set([2,3])}
    gpart = mod.GraphPartition(graph, index)
    assert gpart._node_set == set([0,1,2,3])
    # test raise error if matrix unweighted
    jnk = np.random.random((10,10))
    jnk = np.triu(jnk,1)
    graph = nx.from_numpy_matrix(jnk, nx.Graph(weighted=False))
    npt.assert_raises(ValueError, mod.GraphPartition, graph, index)

def test_find_unconnected_nodes():
    jnk = np.zeros((10,10))
    jnk[:6,:6] = 1
    jnk = np.triu(jnk,1)
    graph = nx.from_numpy_matrix(jnk>0, nx.Graph(weighted=False))
    index = {0:set([0,1,2,3,4,5,6,7,8,9])}
    graph_partition = mod.GraphPartition(graph, index)
    solitary_nodes = graph_partition.find_unconnected_nodes()
    npt.assert_equal(solitary_nodes, [6,7,8,9])

def test_index_as_node_names():
    graph = nx.Graph()
    graph.add_edge('a','b')
    graph.add_edge('c','d')

    ## NOTE network x does not store names as added
    nodes = graph.nodes()

    index = {0:set([0,2]), 1:set([1,3])}
    gpart = mod.GraphPartition(graph, index)
    named_index = gpart.index_as_node_names()
    assert [nodes[0], nodes[2]] in named_index

def test_random_modular_graph_between_fraction():
    """Test for graphs with non-zero between_fraction"""
    # We need to measure the degree within/between modules
    nnods = 120, 240, 360
    nmods = 2, 3, 4
    av_degrees = 8, 10, 16
    btwn_fractions = 0, 0.1, 0.3, 0.5

    for nnod in nnods:
        for nmod in nmods:
            for av_degree in av_degrees:
                for btwn_fraction in btwn_fractions:
                    g = mod.random_modular_graph(nnod, nmod, av_degree,
                                                 btwn_fraction)

                    # First, check the average degree.
                    av_degree_actual = np.mean(list(g.degree().values()))
                    # Since we are generating random graphs, the actual average
                    # degree we get may be off from the reuqested one by a bit.
                    # We allow it to be off by up to 1.
                    #print 'av deg:',av_degree, av_degree_actual  # dbg
                    nt.assert_true (abs(av_degree-av_degree_actual)<1.5,
                          """av deg: %.2f  av deg actual: %.2f -
                          This is a stochastic test - repeat to confirm.""" %
                                          (av_degree, av_degree_actual))

                    # Now, check the between fraction
                    mat = nx.adj_matrix(g)

                    #compute the total number of edges in the real graph
                    nedg = nx.number_of_edges(g)

                     # sanity checks:
                    if nnod%nmod:
                        raise ValueError("nmod must divide nnod evenly")

                    #Compute the of nodes per module
                    nnod_mod = nnod/nmod

                    #compute what the values are in the real graph
                    blocks = [np.ones((nnod_mod, nnod_mod))] * nmod
                    mask = util.diag_stack(blocks)
                    mask[mask==0] = 2
                    mask = np.triu(mask,1)
                    btwn_real = np.sum(mat[mask == 2].flatten())
                    btwn_real_frac = btwn_real / nedg

                    #compare to what the actual values are
                    nt.assert_almost_equal(btwn_fraction,
                                           btwn_real_frac, 1,
                    "This is a stochastic test, repeat to confirm failure")

                    nt.assert_true(abs(btwn_fraction - btwn_real_frac) < 0.1,
                          """btwn fraction: %.2f  real btwn frac: %.2f -
                          This is a stochastic test - repeat to confirm.""" %
                                          (btwn_fraction, btwn_real_frac))

def test_modularity():
    """Test the values that go into the modularity calculation after randomly
    creating a graph"""
    # Given a partition with the correct labels of the input graph, verify that
    # the modularity returns 1

    # We need to measure the degree within/between modules
    nnods = 120, 240, 360
    nmods = 2, 3, 4
    av_degrees = 8, 10, 16

    for nnod in nnods:
        for nmod in nmods:
            for av_degree in av_degrees:
                g = mod.random_modular_graph(nnod, nmod, av_degree)
                #Compute the of nodes per module
                nnod_mod = nnod/nmod
                #Make a "correct" partition for the graph
                part = mod.perfect_partition(nmod,nnod_mod)
                #Make a graphpartition object
                graph_partition = mod.GraphPartition(g,part)
                #call modularity
                mod_meas = graph_partition.modularity()
                mod_true = 1.0 - 1.0 / nmod
                npt.assert_almost_equal(mod_meas, mod_true, 1)

def test_apply_module_merge():
    """Test the GraphPartition operation that merges modules so that it returns
    a change in modularity that reflects the difference between the modularity
    of the new and old parititions"""

    # nnod_mod, av_degrees, nmods
    networks = [ [3, [2], [3, 4]],
                 [4, [2, 3], [2, 4, 6]],
                 [8, [4, 6], [4, 6, 8]] ]

    for nnod_mod, av_degrees, nmods in networks:
        for nmod in nmods:
            nnod = nnod_mod*nmod
            for av_degree in av_degrees:

                g = mod.random_modular_graph(nnod, nmod, av_degree)

                #Make a "correct" partition for the graph
                part = mod.perfect_partition(nmod,nnod/nmod)

                #Make a random partition for the graph
                part_rand = dict()
                while len(part_rand) <= 1: #check if there is only one module
                    part_rand = mod.rand_partition(g)

                #List of modules in the partition
                r_mod=list(range(len(part)))

                #Loop through pairs of modules
                for i in range(1): # DB: why is this necessary?
                    #select two modules to merge
                    mod_per = np.random.permutation(r_mod)
                    m1 = mod_per[0]; m2 = mod_per[1]

                    #make a graph partition object
                    graph_partition = mod.GraphPartition(g,part)

                    #index of nodes within the original module (before merge)
                    n1_init = list(graph_partition.index[m1])
                    n2_init = list(graph_partition.index[m2])
                    n_all_init = n1_init+n2_init

                    #calculate modularity before merging
                    mod_init = graph_partition.modularity()

                    #merge modules
                    merge_module,e1,a1,delta_energy_meas,type,m1,m2,m2 = \
                                  graph_partition.compute_module_merge(m1,m2)

                    graph_part2 = copy.deepcopy(graph_partition)
                    graph_part2.apply_module_merge(m1,m2,merge_module,e1,a1)
                    #index of nodes within the modules after merging
                    n_all = list(graph_part2.index[min(m1,m2)])

                    # recalculate modularity after splitting
                    mod_new = graph_part2.modularity()

                    # difference between new and old modularity
                    delta_energy_true = -(mod_new - mod_init)

                    # Test the measured difference in energy against the
                    # function that calculates the difference in energy
                    npt.assert_almost_equal(delta_energy_meas,
                                                  delta_energy_true)
                    # Check that the list of nodes in the two original modules
                    # is equal to the list of nodes in the merged module
                    n_all_init.sort()
                    n_all.sort()
                    npt.assert_equal(n_all_init, n_all)

                    # Test that the keys are equivalent after merging modules
                    npt.assert_equal(r_mod[:-1],
                                     sorted(graph_part2.index.keys()))

                    # Test that the values in the mod_e and mod_a matrices for
                    # the merged module are correct.
                    npt.assert_equal(graph_part2.mod_e[min(m1,m2)],e1)
                    npt.assert_equal(graph_part2.mod_a[min(m1,m2)],a1)

def test_rename_keys():
    a = {0:0,1:1,2:2,4:4,5:5}
    mod.rename_keys(a, 3)
    npt.assert_equal(a, {0:0,1:1,2:2,3:4,4:5})

    a = {0:0,1:1,3:3,}
    mod.rename_keys(a, 2)
    npt.assert_equal(a, {0:0,1:1,2:3})

    # If called with the last key in dict, it should leave the input alone
    a = {0:0,1:1,2:2,3:3}
    mod.rename_keys(a, 3)
    npt.assert_equal(a, a)

def danon_benchmark():
    """This test comes from Danon et al 2005. It will create the line plot of
    Mututal Information vs. betweenness fraction to assess the performance of
    the simulated annealing algorithm."""
    networks = [[32, [16], [6]]]
    btwn_fracs = [float(i)/100 for i in range(0,80,3)]

    temperature = 0.1
    temp_scaling = 0.9995
    tmin=1e-4

    num_reps = range(1)
    mi_arr=np.empty((len(btwn_fracs),len(num_reps)))

    #keep time
    for rep in num_reps:
        t1 = time.clock()
        for nnod_mod, av_degrees, nmods in networks:
            for nmod in nmods:
                nnod = nnod_mod*nmod
                for av_degree in av_degrees:
                    x_mod = []
                    for ix,btwn_frac in enumerate(btwn_fracs):
                        print('btwn_frac: ',btwn_frac)
                        g = mod.random_modular_graph(nnod, nmod, av_degree,btwn_frac)
                        #Compute the # of nodes per module
                        nnod_mod = nnod/nmod
                        #Make a "correct" partition for the graph
                        ppart = mod.perfect_partition(nmod,nnod_mod)

                        graph_out, graph_dict =mod.simulated_annealing(g,
                        temperature = temperature,temp_scaling = temp_scaling,
                        tmin=tmin, extra_info = True)

                        #print "SA partition",graph_out.index
                        mi = mod.mutual_information(ppart,graph_out.index)
                        t2 = time.clock()
                        print('Elapsed time: ', (float(t2-t1)/60), ' minutes')
                        print('partition similarity: ', mi)
                        mi_arr[ix,rep] = mi
                        ## plot_partition(g,graph_out.index,'mi: '+ str(mi),'danon_test_6mod'+str(btwn_frac)+'_graph.png')
                        x_mod.append(betweenness_to_modularity(g,ppart))

                    ## mi_arr_avg = np.mean(mi_arr,1)
                    ## plt.figure()
                    ## plt.plot(btwn_fracs,mi_arr_avg)
                    ## plt.xlabel('Betweenness fraction')
                    ## plt.ylabel('Mutual information')
                    ## plt.savefig('danon_test_6mod/danontest_btwn.png')

                    ## plt.figure()
                    ## plt.plot(x_mod,mi_arr_avg)
                    ## plt.xlabel('Modularity')
                    ## plt.ylabel('Mutual information')
                    ## plt.savefig('danon_test_6mod/danontest_mod.png')

    #plt.figure()
    #plt.plot(graph_dict['energy'], label = 'energy')
    #plt.plot(graph_dict['temperature'], label = 'temperature')
    #plt.xlabel('Iteration')

    return mi_arr

def SA():
    """ Test the simulated annealing script"""
    #nnod_mod, av_degrees, nmods
    #networks = [ [4, [2, 3], [2, 4, 6]]]#,
    #networks =  [ [8, [4, 6], [4, 6, 8]]]
    #networks = [[40, [20], [2]]]
    networks = [[32, [16], [4]]]
    #networks = [[64, [12], [6]]]
    btwn_fracs = [0]
    temperature = 10
    temp_scaling = 0.9995
    tmin=1e-4
    nochange_ratio_min=0.01

    #keep time

    for nnod_mod, av_degrees, nmods in networks:
        for nmod in nmods:
            nnod = nnod_mod*nmod
            for av_degree in av_degrees:
                for btwn_frac in btwn_fracs:
                    t1=time.clock()
                    g = mod.random_modular_graph(nnod, nmod, av_degree,btwn_frac)
                    #Compute the # of nodes per module
                    nnod_mod = nnod/nmod
                    #Make a "correct" partition for the graph
                    ppart = mod.perfect_partition(nmod,nnod_mod)

                    graph_out, energy_array, rej_array, temp_array =mod.simulated_annealing(g,
                    temperature = temperature,temp_scaling = temp_scaling,
                    tmin=tmin, nochange_ratio_min = nochange_ratio_min)

                    print("perfect partition", ppart)
                    print("SA partition",graph_out.index)

                    t2 = time.clock()
                    print('Elapsed time: ', float(t2-t1)/60, 'minutes')
                    print('partition similarity: ',
                          mod.mutual_information(ppart,graph_out.index))
                    return graph_out, g, energy_array, rej_array, ppart, temp_array

def test_mutual_information_simple():
    """MI computations with hand-validated values.
    """
    from math import log
    # Define two simple partitions that are off by one assignment
    a = {0:[0, 1], 1:[2, 3], 2:[4, 5]}
    b = {0:[0, 1], 1:[2, 3, 4], 2:[5]}
    N_true = np.array([ [2,0,0], [0,2,0], [0,1,1] ], dtype=float)
    N = mod.confusion_matrix(a, b)
    # test confusion matrix
    npt.assert_equal(N, N_true)
    # Now compute mi by hand
    num = -6*log(3)-4*log(2)
    den = -(3*log(2)+8*log(3)+log(6))
    mi_true = num/den
    mi = mod.mutual_information(a, b)
    npt.assert_almost_equal(mi, mi_true)
    # Let's now flip the labels and confirm that the computation is impervious
    # to module labels
    b2 = {2:[0, 1], 0:[2, 3, 4], 1:[5]}
    npt.assert_almost_equal(mod.mutual_information(b, b2), 1)
    npt.assert_almost_equal(mod.mutual_information(a, b2), mi)

def test_mutual_information_empty():
    """Validate that empty modules don't affect MI.
    """
    # Define two simple partitions that are off by one assignment
    a = {0:[0, 1], 1:[2, 3], 2:[4, 5]}
    b = {0:[0, 1], 1:[2, 3], 2:[4, 5], 3:[]}

    try:
        mod.mutual_information(a, b)
    except ValueError as e:
        nt.assert_equals(e.args[0], "Empty module in second partition.")

    try:
        mod.mutual_information(b, a)
    except ValueError as e:
        nt.assert_equals(e.args[0], "Empty module in first partition.")

def test_mutual_information():
    """ Test the function which returns the mutual information in two
    partitions

    XXX - This test is currently incomplete - it only checks the most basic
    case of MI(x, x)==1, but doesn't do any non-trivial checks.
    """

    # nnod_mod, av_degrees, nmods
    networks = [ [4, [2, 3], [2, 4, 6]],
                 [8, [4, 6], [4, 6, 8]],
                 [40, [20], [2]] ]

    for nnod_mod, av_degrees, nmods in networks:
        for nmod in nmods:
            nnod = nnod_mod*nmod
            for av_degree in av_degrees:
                #make a graph object
                g = mod.random_modular_graph(nnod, nmod, av_degree)

                #Compute the of nodes per module
                nnod_mod = nnod//nmod
                #Make a "correct" partition for the graph
                ppart = mod.perfect_partition(nmod,nnod_mod)

                #graph_out, mod_array =mod.simulated_annealing(g, temperature =
                #temperature,temp_scaling = temp_scaling, tmin=tmin)

                #test the perfect case for now: two of the same partition
                #returns 1
                mi_orig  = mod.mutual_information(ppart,ppart)
                npt.assert_equal(mi_orig,1)

                #move one node and test that mutual_information comes out
                #correctly
                graph_partition = mod.GraphPartition(g,ppart)
                graph_partition.node_update(0,0,1)

                mi = mod.mutual_information(ppart,graph_partition.index)
                npt.assert_array_less(mi, mi_orig)
                ## NOTE: CORRECTNESS NOT TESTED YET

                #merge modules and check that mutual information comes out
                #correctly/lower
                graph_partition2 = mod.GraphPartition(g,ppart)
                merged_module, e_new, a_new, d,t,m1,m2,x = graph_partition2.compute_module_merge(0,1)
                graph_partition2.apply_module_merge(m1,m2,merged_module,e_new,a_new)
                mi2 = mod.mutual_information(ppart,graph_partition2.index)
                npt.assert_array_less(mi2,mi_orig)
                ## NOTE: CORRECTNESS NOT TESTED YET

                #split modules and check that mutual information comes out
                #correclty/lower
                graph_partition3 = mod.GraphPartition(g,ppart)
                n1 = set(list(graph_partition3.index[0])[::2])
                n2 = set(list(graph_partition3.index[0])[1::2])

                (split_modules, e_new,
                 a_new, d, t, m,
                 n1,n2) = graph_partition3.compute_module_split(0,n1,n2)
                graph_partition3.apply_module_split(m, n1, n2,
                                                    split_modules,
                                                    e_new, a_new)
                mi3 = mod.mutual_information(ppart, graph_partition3.index)
                npt.assert_array_less(mi3,mi_orig)
                ## NOTE: CORRECTNESS NOT TESTED YET

def test_random_mod():
    """ Test the GraphPartition operation that selects random modules
    to merge and split
    XXX not working yet"""

    #nnod_mod, av_degrees, nmods
    networks = [ [4, [2, 3], [2, 4, 6]],
                 [8, [4, 6], [4, 6, 8]] ]

    n_iter = 100
    for nnod_mod, av_degrees, nmods in networks:
        for nmod in nmods:
            nnod = nnod_mod*nmod
            for av_degree in av_degrees:

                g = mod.random_modular_graph(nnod, nmod, av_degree)
                part = dict()
                while (len(part) <= 1) or (len(part) == nnod):
                    part = mod.rand_partition(g)

                graph_partition = mod.GraphPartition(g,part)

                for i in range(n_iter):
                    graph_partition.random_mod()

                    #check that the partition has > 1 modules
                    true = len(graph_partition)>1
                    npt.assert_equal(true,1)

                    #check that the partition has < nnod modules
                    true = len(graph_partition)<nnod
                    npt.assert_equal(true,1)

def test_random_nod():
    """ Test the GraphPartition operation that selects random nodes to move
    between modules """

def test_decide_if_keeping():
    """ Test the function which decides whether or not to keep the new
    partition"""

    dEs = [-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5]
    temperatures = [.01,1,10,100]
    iter = 1000
    tolerance = 1

    for temp in temperatures:
        for dE in dEs:
            keep_list = np.empty(iter)
            for i in range(iter):
                keep_list[i] = float(mod.decide_if_keeping(dE,temp))

            if dE <= 0:
                keep_correct = np.ones(iter)
                npt.assert_equal(keep_list,keep_correct)
            else:
                mean_keep = np.mean(keep_list)
                mean_correct = math.exp(-dE/temp)
                npt.assert_almost_equal(mean_keep,mean_correct, tolerance)

def test_sim_anneal_simple():
    """Very simple simulated_annealing test with a small network"""

    #
    nnod, nmod, av_degree, btwn_frac = 24, 3, 4, 0
    g = mod.random_modular_graph(nnod, nmod, av_degree, btwn_frac)

    #Compute the # of nodes per module
    ## nnod_mod = nnod/nmod
    #Make a "correct" partition for the graph
    ## ppart = mod.perfect_partition(nmod,nnod_mod)

    temperature = 10
    temp_scaling = 0.95
    tmin=1

    graph_out, graph_dict = mod.simulated_annealing(g,
               temperature = temperature, temp_scaling = temp_scaling,
               tmin=tmin, extra_info = True, debug=True)

    # Ensure that there are no empty modules
    util.assert_no_empty_modules(graph_out.index)

    ## mi = mod.mutual_information(ppart, graph_out.index)
    #nt.assert_equal(mi, 1)

def test_apply_module_split():
    """Test the GraphPartition operation that splits modules so that it returns
    a change in modularity that reflects the difference between the modularity
    of the new and old parititions.
    Also test that the module that was split now contains the correct nodes,
    the correct modularity update, the correct energy,and that no empty modules
    result from it."""

    # nnod_mod, av_degrees, nmods
    networks = [ [3, [2], [2, 3, 4]],
                 [4, [2, 3], [2, 4, 6]],
                 [8, [4, 6], [4, 6, 8]] ]

    for nnod_mod, av_degrees, nmods in networks:
        for nmod in nmods:
            nnod = nnod_mod*nmod
            for av_degree in av_degrees:

                g = mod.random_modular_graph(nnod, nmod, av_degree)

                # Make a "correct" partition for the graph
                ## part = mod.perfect_partition(nmod,nnod/nmod)

                # Make a random partition for the graph
                part_rand = mod.rand_partition(g, nnod//2)

                #List of modules in the partition that have two or more nodes
                r_mod = []
                for m, nodes in part_rand.items():
                    if len(nodes)>2:
                        r_mod.append(m)

                # Note: The above can be written in this more compact, if
                # slightly less intuitively clear, list comprehension:
                # r_mod = [ m for m, nodes in part_rand.items() if
                #          len(nodes)>2 ]

                #Module that we are splitting
                for m in r_mod:

                    graph_partition = mod.GraphPartition(g,part_rand)

                    #index of nodes within the original module (before split)
                    n_init = list(graph_partition.index[m])

                    #calculate modularity before splitting
                    mod_init = graph_partition.modularity()

                    # assign nodes to two groups
                    n1_orig,n2_orig = graph_partition.determine_node_split(m)

                    # make sure neither of these is empty
                    nt.assert_true(len(n1_orig)>= 1)
                    nt.assert_true(len(n2_orig)>= 1)

                    #make sure that there are no common nodes between the two
                    node_intersection = set.intersection(n1_orig,n2_orig)
                    nt.assert_equal(node_intersection,set([]))

                    #make sure that sum of the two node sets equals the
                    #original set
                    node_union = set.union(n1_orig,n2_orig)
                    npt.assert_equal(np.sort(list(node_union)),np.sort(n_init))

                    # split modules
                    split_modules,e1,a1,delta_energy_meas,type,m,n1,n2 = \
                               graph_partition.compute_module_split(m,n1_orig,n2_orig)

                    #note: n1 and n2 are output from this function (as well as
                    #inputs) because the function is called from within another
                    #(rand_mod) def but then output to the simulated script, so
                    #the node split needs to be passed along.

                    #as a simple confirmation, can make sure they match
                    npt.assert_equal(n1_orig,n1)
                    npt.assert_equal(n2_orig,n2)

                    #split_moduels should be a dictionary with two modules
                    #(0,1) that contain the node sets n1 and n2 respectively.
                    #test this.
                    npt.assert_equal(split_modules[0],n1)
                    npt.assert_equal(split_modules[1],n2)

                    #make a new graph partition equal to the old one and apply
                    #the module split to it (graph_part2)
                    graph_part2 = copy.deepcopy(graph_partition)
                    graph_part2.apply_module_split(m,n1,n2,split_modules,e1,a1)

                    #make a third graph partition using only the original graph
                    #and the partition from graph_part2
                    graph_part3 = mod.GraphPartition(g,graph_part2.index)

                    #index of nodes within the modules after splitting
                    n1_new = list(graph_part2.index[m])
                    n2_new = list(graph_part2.index[len(graph_part2)-1])
                    n_all = n1_new + n2_new

                    # recalculate modularity after splitting
                    mod_new = graph_part2.modularity()
                    mod_new_3 = graph_part3.modularity()

                    # difference between new and old modularity
                    delta_energy_true = -(mod_new - mod_init)

                    # Test that the measured change in energy by splitting a
                    # module is equal to the function output from module_split
                    npt.assert_almost_equal(delta_energy_meas,
                                                  delta_energy_true)

                    # Test that the nodes in the split modules are equal to the
                    # original nodes of the module
                    nt.assert_equal(sorted(list(n1)), sorted(n1_new))
                    nt.assert_equal(sorted(list(n2)), sorted(n2_new))

                    n_init.sort()
                    n_all.sort()
                    # Test that the initial list of nodes in the module are
                    # equal to the nodes in m1 and m2 (split modules)
                    npt.assert_equal(n_init,n_all)

                    # Test that the computed modularity found when
                    # apply_module_split is used is equal to the modularity you
                    # would find if using that partition and that graph
                    npt.assert_almost_equal(mod_new,mod_new_3)

                    # Check that there are no empty modules in the final
                    # partition
                    for m in graph_part2.index:
                        nt.assert_true(len(graph_part2.index[m]) > 0)

def test_apply_node_move():
    """Test the GraphPartition operation that moves a single node so that it
    returns a change in modularity that reflects the difference between the
    modularity of the new and old parititions"""

    # nnod_mod, av_degrees, nmods
    #networks = [ [3, [2], [2, 3, 4]],
    #             [4, [2, 3], [2, 4, 6]],
    #             [8, [4, 6], [4, 6, 8]] ]
    networks = [ [4, [2, 3], [2, 4, 6]],
                 [8, [4, 6], [4, 6, 8]] ]

    for nnod_mod, av_degrees, nmods in networks:
        for nmod in nmods:
            nnod = nnod_mod*nmod
            for av_degree in av_degrees:
                print(nnod_mod,nmod,av_degree)
                g = mod.random_modular_graph(nnod, nmod, av_degree)

                #Make a "correct" partition for the graph
                #part = mod.perfect_partition(nmod,nnod/nmod)

                #Make a random partition for the graph
                part_rand = dict()
                while len(part_rand) <= 1: #check if there is only one module
                    part_rand = mod.rand_partition(g)

                #List of modules in the partition
                r_mod=range(len(part_rand))

                #Make a graph_partition object
                graph_partition = mod.GraphPartition(g,part_rand)

                #select two modules to change node assignments
                mod_per = np.random.permutation(r_mod)
                m1 = mod_per[0]; m2 = mod_per[1]
                while len(graph_partition.index[m1]) <= 1:
                    mod_per = np.random.permutation(r_mod)
                    m1 = mod_per[0]
                    m2 = mod_per[1]

                #pick a random node to move between modules m1 and m2
                node_list=list(graph_partition.index[m1])
                nod_per = np.random.permutation(node_list)
                n = nod_per[0]

                #list of nodes within the original modules (before node move)
                ## n1_init = list(nod_per) #list(graph_partition.index[m1])
                ## n2_init = list(graph_partition.index[m2])
                ## n1_new = copy.deepcopy(n1_init)
                ## n2_new = copy.deepcopy(n2_init)

                # calculate modularity before node move
                mod_init = graph_partition.modularity()

                # move node from m1 to m2
                node_moved_mods,e1,a1,delta_energy_meas,n,m1,m2 = \
                              graph_partition.compute_node_update(n,m1,m2)

                graph_part2 = copy.deepcopy(graph_partition)
                m2_new = graph_part2.apply_node_update(n,m1,m2,node_moved_mods,e1,a1)

                #if the keys get renamed, the m1,m2 numbers are no longer the same

                #test that m2 now contains n
                nt.assert_true(n in graph_part2.index[m2_new])
                #if n not in graph_part2.index[m2_new]:
                #    1/0

                # recalculate modularity after splitting
                mod_new = graph_part2.modularity()

                # difference between new and old modularity
                delta_energy_true = -(mod_new - mod_init)
                #print delta_energy_meas,delta_energy_true

                # Test that the measured change in energy is equal to the true
                # change in energy calculated in the node_update function
                npt.assert_almost_equal(delta_energy_meas,delta_energy_true)

def test_adjust_partition():
    e = np.loadtxt(os.path.join(os.path.dirname(__file__), 'jazz.net'),
                   skiprows=3, dtype=int)[:, :2] - 1
    g = nx.Graph()
    g.add_edges_from(e)

    p0 = mod.newman_partition(g)
    #p1 = mod.adjust_partition(g, p0, max_iter=6)

    ## This doesnt test what we want to test FIXME
    #npt.assert_(p0 > 0.38)
    #npt.assert_(p1 > 0.42)

def test_empty_graphpartition():
    g = nx.Graph()
    g.add_node(1)
    npt.assert_raises(ValueError, mod.GraphPartition, g, {1: set(g.nodes())})

def test_badindex_graphpartition():
    """ when making a GraphPArtition, check index is valid"""
    ## index should be dict of sets
    e = np.loadtxt(os.path.join(os.path.dirname(__file__), 'jazz.net'),
                                skiprows=3, dtype=int)[:, :2] - 1
    g = nx.Graph()
    g.add_edges_from(e)
    index = {0: set(g.nodes()[:100]), 1: set(g.nodes()[100:])}
    gp = mod.GraphPartition(g, index)
    nt.assert_true(gp.index == index)
    npt.assert_raises(TypeError, mod.GraphPartition, g, {0: g.nodes()})
    npt.assert_raises(ValueError, mod.GraphPartition, g,
                      {0:set(g.nodes()[:-1])})
    npt.assert_raises(TypeError, mod.GraphPartition, g, g.nodes())

def test_newman_partition():
    """ Test Newman Partition function """
    tmpmat = np.random.random((10,10))
    tmpmat[tmpmat < .5] = 0
    graph = nx.from_numpy_matrix(tmpmat, nx.Graph(weighted = False))
    npt.assert_raises(ValueError, mod.newman_partition, graph)
    tmpmat[:] = 0
    # test that no edges raises error (from GraphPartition)
    graph = nx.from_numpy_matrix(tmpmat, nx.Graph(weighted = False))
    npt.assert_raises(ValueError, mod.newman_partition, graph)
    tmpmat[:] = 1
    util.fill_diagonal(tmpmat, 0)
    graph = nx.from_numpy_matrix(tmpmat, nx.Graph(weighted=False))
    part = mod.newman_partition(graph)
    ## if all edges are connected expect only one partition
    expected_part = {0: set([0,1,2,3,4,5,6,7,8,9])}
    nt.assert_equal(part.index, expected_part)

if __name__ == "__main__":
    npt.run_module_suite()

max_stars_repo_name: daniel-ge/brainx

max_issues_repo_name: daniel-ge/brainx

max_forks_repo_name: daniel-ge/brainx

reason: import numpy,import networkx


text:
from PIL import Image
import torch
import torch.nn as nn
import segmentation_models_pytorch as smp
from collections import namedtuple
import numpy as np
from scipy.ndimage.morphology import binary_fill_holes, binary_closing, binary_erosion
import cv2
from skimage.exposure import match_histograms
from typing import List

MODELPATH = [
    "../model/skull_stripper_timm-efficientnet-b4_3D.pth",
    "../model/skull_stripper_resnet18_3D.pth",
]

MODELS = [
    smp.Unet(
        encoder_name="timm-efficientnet-b4",
        encoder_weights="imagenet",
        in_channels=1,
        classes=1,
    ),
    smp.Unet(
        encoder_name="resnet18",
        encoder_weights="imagenet",
        in_channels=1,
        classes=1,
    ),
]

Segmented = namedtuple("Segmented", "brain mask")

# the reference images use to make histograms matching
REF = np.asarray(Image.open("../example/y_slice_9.png")) / 255

class Brook:
    def __init__(self, model_path: List[str] = MODELPATH, models: List[str] = MODELS):
        """[Core class of pybrook package. That class strip a 2D SLICE OF NUMPY TENSOR]
        Args:
            model_path (List[str], optional): [Path of pretrained model weight, SHOULD BE
            CHANGE by the user]. Defaults to MODELPATH.
            models (List[str], optional): [Models architecture (should be in the following
            order [timm-efficientnet-b4,resnet18])]. Defaults to MODELS.
        """
        self.models = models
        for model_idx in range(len(models)):
            self.models[model_idx].load_state_dict(torch.load(model_path[model_idx]))
            self.models[model_idx].eval()

    def strip(
        self, x: np.ndarray, pre_process: bool = False, post_process: bool = False
    ) -> Segmented:
        """[Strip a slice of array,using stacking of 2 pretrained models.]

        Args:
            x (np.ndarray): [2D tensor (slice of the MRI)]
            pre_process (bool, optional): [True to make histogram matching beetween new MRI and
            training MRI]. Defaults to False.
            post_process (bool, optional): [True to apply morphological operator to fill holes
            and remove artefacts]. Defaults to False.

        Returns:
            Segmented: [Namedtuple that contain, the segmented brain images and the segmentation mask]
        """

        H, W = x.shape
        cp = x.copy()
        if pre_process:
            ref = cv2.resize(REF, dsize=(H, W), interpolation=cv2.INTER_CUBIC)
            x = match_histograms(x, ref, multichannel=False)
        if x.max() > 1:
            x = x / x.max()
        x = torch.tensor(x, dtype=torch.float32).unsqueeze(0).unsqueeze(0)
        mask = torch.ones_like(x)
        for model in self.models:
            logits = model(x)
            mask *= nn.Sigmoid()(logits)

        mask = mask ** (1 / len(self.models))  # geometric mean
        mask = torch.where(mask < 0.5, 0, 1)
        mask = mask.squeeze(0).squeeze(0).numpy()
        # remove little artefact on edge
        if post_process:
            mask = binary_closing(
                mask, structure=np.ones((H // 3, W // 3))
            )  # fill big holes in mask
            mask = binary_fill_holes(mask).astype(int)  # fill small holes in mask
            mask = binary_erosion(mask, structure=np.ones((H // 25, W // 25)))

        return Segmented(cp * mask, mask)


text:
from sklearn.externals import joblib
import numpy as np

def getData(data):
    labeled_features = data.item(0)["labeled_features"]
    frequencies = data.item(0)["frequencies"]
    labels = []
    features = []
    for labeled_feature in labeled_features:
        labels.append(labeled_feature['label'])
        features.append(labeled_feature['features'].flatten())
    return labels, features

model_file = 'QDA.pkl' #You can replace this by any model produced by multipleClassifiers file
#preprocess_recordings('19_06_05_07_2017_freq_19.mat')
data_test2 = np.load('19_06_05_07_2017_freq_19.mat.npy')
labels2, features2 = getData(data_test2) #Replace this with online data
labels2 = np.array(labels2)

classifier = joblib.load(model_file)
out_test = classifier.predict(features2)
accuracy = 100 - 100 * np.sqrt(np.sum(np.power(labels2 - out_test, 2))) / len(labels2)
print("Classifier accuracy: " + str(accuracy))

max_stars_repo_name: gumpy-hybridBCI/gumpy-Realtime

max_issues_repo_name: swxie/SSVEP-Brain-Computer-Interface

max_forks_repo_name: swxie/SSVEP-Brain-Computer-Interface


text:
import numpy as np
import cv2
import time
import Motor
import sys

NAO_ENCONTROU = 1
FINALIZOU = 2
FATOR_MOTOR = 80

lower_red1 = np.array([0,100,80])#np.array([30,150,50])
upper_red1 = np.array([10,256,256])#np.array([255,255,180])
lower_red2 = np.array([170,100,80])#np.array([30,150,50])
upper_red2 = np.array([180,256,256])#np.array([255,255,180])

#Defica cam como objeto da camera
cam = cv2.VideoCapture(0)
##cam = cv2.VideoCapture('/home/pi/Projeto RAMLC/RMA Oficial/V_20170812_101753.mp4')
#resolução da camera em 320x240
heigth=320
width=240

cam.set(3,heigth)
cam.set(4,width)

#motor = Motor.motores(24,23,17,22,26,19)

#corte da visao da camera
r = 0
h = 240
c = 80
w = 200

def deploy():
    cam.release()
    cv2.destroyAllWindows()

    #motor.clean()
    #del motor

def varrerAteFaixa():
    Executa = True
    Tentativa = 0
    left = 0
    right = 0

    ret, frame = cam.read()
    #Corta a imagem capiturada
    roi = frame[r:r+h, c:c+w]

    #Da imagem trasforma as configurações de cores de RGB para HSV
    hsv=cv2.cvtColor(roi,cv2.COLOR_BGR2HSV)
    #Efeito de blur para retirar imperfeição
    #blur = cv2.blur(hsv,(5,5))
    #blur = cv2.GaussianBlur(hsv,(3,3), 0)

    #Cria uma mascara para colher imagem no intervalo de HSV conforme os dois array
    image_mask1=cv2.inRange(hsv,lower_red1,upper_red1)
    image_mask2=cv2.inRange(hsv,lower_red2,upper_red2)
    image_mask = image_mask1 + image_mask2
    output=cv2.bitwise_and(roi,roi,mask=image_mask)

    #Efeito de blur para retirar imperfeição
    smoothed = cv2.blur(output,(1,1))
    #Realiaza verificação de linhas
    edges = cv2.Canny(smoothed,50,150,apertureSize=3)

    lines = cv2.HoughLines(edges,1,np.pi/180,80)
    if lines is not None:
        Tentativa = 0
        for rho,theta in lines[0]:
            a = np.cos(theta)
            b = np.sin(theta)
            x0 = a*rho
            y0 = b*rho
            x1 = int(x0 + 1000*(-b))#direita
            y1 = int(y0 + 1000*(a))
            x2 = int(x0 - 1000*(-b))#Equerda
            y2 = int(y0 - 1000*(a))
            cv2.line(edges,(x1,y1),(x2,y2),(0,0,255),2)

            left = y2
            right= y1
            #if(y2<=0):
            #    left = y2*-1
            #if(y1<= 0):
            #    right= y1*-1

    #cv2.imshow('Original',frame)
    #cv2.imshow('Canny',edges)

    #cv2.imshow('filter', roi)

    #cv2.imshow('Output',edges)

    return right, left

#########################
#         MAIN          #
#########################

#motor.start(100,255)

#faz algumas leituras de frames antes de consierar a analise
#motivo: algumas camera podem demorar mais para se "acosumar a luminosidade" quando ligam, capturando frames consecutivos com muita variacao de luminosidade. Para nao levar este efeito ao processamento de imagem, capturas sucessivas sao feitas fora do processamento da imagem, dando tempo para a camera "se acostumar" a luminosidade do ambiente
#for i in range(0,20):
#    ret, Frame = cam.read()

while ( True ):
    try:
        right, left = varrerAteFaixa()
        if(True):
            #Esquerda y2    #Direita y1
            print("ESQUERDA: " + str(left), "DIREITA: " + str(right))

            if( ((right - left) >= 0) and (FATOR_MOTOR >= (right - left) and (right - left) <= FATOR_MOTOR)):
                print("Motor frente 1")

                #motor.frente(0.2)
                #motor.frente(10)

            elif( ((left - right) >= 0) and (FATOR_MOTOR >= (left - right) and (left - right) <= FATOR_MOTOR)):
                print("Motor frente 2")

                #motor.frente(0.2)
                #motor.frente(10)
            elif(right < left):
                print("Motor frente esquerda")
                #motor.frenteEsq(0.2)

            elif(right > left):
                print("Motor frente direita")
                #motor.frenteDir(0.2)
            elif( right <= 0  or left <= 0 ):
                print("Parar Motor 1")
                #motor.stop()

            ##elif( y1 < 0 or y2 < 0 ):
             ##   print("Localizado ponto vertical")

            #if(status == NAO_ENCONTROU)
            print("EXECUTANDO NOVO CICLO")
            if cv2.waitKey(33) == 27:
               break

    except Exception as e:
        print("Houve um problema inesperado, operação reexecutará em 10 segundos, segue detalhe do erro:")
        print (str(e))
        time.sleep(10)
        sys.exc_info()
        pass
    finally:
        deploy()
        exit(1)

max_stars_repo_name: viniciuspgng/Herbie

max_issues_repo_name: viniciuspgng/Herbie

max_forks_repo_name: viniciuspgng/Herbie


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# import the necessary packages

import numpy as np
import dlib
import cv2
import imutils
from imutils import face_utils
from scipy.spatial import distance as dist

# 论文 Real-Time Eye Blink Detection using Facial Landmarks
def eye_aspect_ratio(eye):
    # 眼睛垂直距离
    A = dist.euclidean(eye[1], eye[5])
    B = dist.euclidean(eye[2], eye[4])
    # 眼睛水平距离
    C = dist.euclidean(eye[0], eye[3])
    # 使用两者的比值来度量眨眼
    ear = (A + B) / (2.0 * C)
    return ear

# 眼睛长宽比阈值
EYE_AR_THRESH = 0.2
# 噪声阈值 -- 可自行调节
# 自测发现 2 帧时检测眨眼效果比 3 帧好，但噪声较大
EYE_AR_CONSEC_FRAMES = 2
# 疲劳阈值
EYE_AR_CONSEC_FRAMES_SLEEP = 10

# 统计噪声
COUNTER = 0
# 眨眼次数
TOTAL = 0

detector = dlib.get_frontal_face_detector()
predictor = dlib.shape_predictor('shape_predictor_68_face_landmarks.dat')

# 左右眼索引
(lStart, lEnd) = face_utils.FACIAL_LANDMARKS_68_IDXS["left_eye"]
(rStart, rEnd) = face_utils.FACIAL_LANDMARKS_68_IDXS["right_eye"]

cap = cv2.VideoCapture(0)

while True:
    _, frame = cap.read()

    # 转换成灰度图
    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    rects = detector(gray, 0)

    for rect in rects:
        # 找到 bindingbox
        shape = predictor(gray, rect)
        shape = face_utils.shape_to_np(shape)

        leftEye = shape[lStart:lEnd]
        rightEye = shape[rStart:rEnd]
        leftEAR = eye_aspect_ratio(leftEye)
        rightEAR = eye_aspect_ratio(rightEye)
        # 双眼平均
        ear = (leftEAR + rightEAR) / 2.0

        leftEyeHull = cv2.convexHull(leftEye)
        rightEyeHull = cv2.convexHull(rightEye)
        cv2.drawContours(frame, [leftEyeHull], -1, (0, 255, 0), 1)
        cv2.drawContours(frame, [rightEyeHull], -1, (0, 255, 0), 1)

        # 疲劳检测
        if COUNTER >= EYE_AR_CONSEC_FRAMES_SLEEP:
            cv2.putText(frame, "DANGER!!!", (70,320), 0, 3, (0,0,255), 5)

        if ear < EYE_AR_THRESH:
            COUNTER += 1
        else:
            # 消除抖动噪声
            if COUNTER >= EYE_AR_CONSEC_FRAMES:
                TOTAL += 1

            COUNTER = 0

    cv2.putText(frame, "COUNTER: {}".format(COUNTER), (150, 30),
    cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 0, 255), 2)
    cv2.putText(frame, "Blinks: {}".format(TOTAL), (10, 30),
    cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 0, 255), 2)
    cv2.putText(frame, "EAR: {:.2f}".format(ear), (300, 30),
    cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 0, 255), 2)

    cv2.imshow("Frame", frame)

    if cv2.waitKey(33) & 0xFF == ord('q'):
        break

cv2.destroyAllWindows()

max_stars_repo_name: chinasanshi/faceLandmarkRect

max_issues_repo_name: chinasanshi/faceLandmarkRect

max_forks_repo_name: chinasanshi/faceLandmarkRect


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 20 07:58:00 2019

@author: davidrounce
"""

import os
import cartopy
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import numpy as np
import pandas as pd

#pd.read_csv('')

#%% X-Y PLOT
#ds = pd.read_csv(...)

# X,Y values
x_values = [0,1,2,3]
y_values = [2,-1,1,5]
y2_values = [3, 2, 3, 0]

# Set up your plot (and/or subplots)
fig, ax = plt.subplots(1, 1, squeeze=False, sharex=False, sharey=False, gridspec_kw = {'wspace':0.4, 'hspace':0.15})

# Plot
#  zorder controls the order of the plots (higher zorder plots on top)
#  label used to automatically generate legends (legends can be done manually for more control)
ax[0,0].plot(x_values, y_values, color='k', linewidth=1, zorder=2, label='plot1')
ax[0,0].plot(x_values, y2_values, color='b', linewidth=1, zorder=2, label='plot2')

# Fill between
#  fill between is useful for putting colors between plots (e.g., error bounds)
#ax[0,0].fill_between(x_values, y_values, y2_values, facecolor='k', alpha=0.2, zorder=1)

# Text
#  text can be used to manually add labels or to comment on plot
#  transform=ax.transAxes means the x and y are between 0-1
ax[0,0].text(0.5, 0.99, '[insert text]', size=10, horizontalalignment='center', verticalalignment='top',
             transform=ax[0,0].transAxes)

# X-label
ax[0,0].set_xlabel('X LABEL', size=12)
#ax[0,0].set_xlim(time_values_annual[t1_idx:t2_idx].min(), time_values_annual[t1_idx:t2_idx].max())
#ax[0,0].xaxis.set_tick_params(labelsize=12)
#ax[0,0].xaxis.set_major_locator(plt.MultipleLocator(50))
#ax[0,0].xaxis.set_minor_locator(plt.MultipleLocator(10))
#ax[0,0].set_xticklabels(['2015','2050','2100'])

# Y-label
ax[0,0].set_ylabel('Y LABEL', size=12)
#ax[0,0].set_ylim(0,1.1)
#ax[0,0].yaxis.set_major_locator(plt.MultipleLocator(0.2))
#ax[0,0].yaxis.set_minor_locator(plt.MultipleLocator(0.05))

# Tick parameters
#  controls the plotting of the ticks
#ax[0,0].yaxis.set_ticks_position('both')
#ax[0,0].tick_params(axis='both', which='major', labelsize=12, direction='inout')
#ax[0,0].tick_params(axis='both', which='minor', labelsize=12, direction='inout')

# Example Legend
# Option 1: automatic based on labels
ax[0,0].legend(loc=(0.05, 0.05), fontsize=10, labelspacing=0.25, handlelength=1, handletextpad=0.25, borderpad=0,
               frameon=False)
# Option 2: manually define legend
#leg_lines = []
#labels = ['plot1', 'plot2']
#label_colors = ['k', 'b']
#for nlabel, label in enumerate(labels):
#    line = Line2D([0,1],[0,1], color=label_colors[nlabel], linewidth=1)
#    leg_lines.append(line)
#ax[0,0].legend(leg_lines, labels, loc=(0.05,0.05), fontsize=10, labelspacing=0.25, handlelength=1,
#               handletextpad=0.25, borderpad=0, frameon=False)

# Save figure
#  figures can be saved in any format (.jpg, .png, .pdf, etc.)
fig.set_size_inches(4, 4)
figure_fp = os.getcwd() + '/../Output/'
if os.path.exists(figure_fp) == False:
    os.makedirs(figure_fp)
figure_fn = 'example_plot_xy.png'
fig.savefig(figure_fp + figure_fn, bbox_inches='tight', dpi=300)

#%% MAP PLOT
xtick = 5
ytick = 5
xlabel = 'Longitude [$^\circ$]'
ylabel = 'Latitude [$^\circ$]'
labelsize = 12
west = -180
east = 180
south = -90
north = 90

rgiO1_shp_fn = os.getcwd() + '/../RGI/rgi60/00_rgi60_regions/00_rgi60_O1Regions.shp'

# Time, Latitude, Longitude
lons = [7.5,86.8,43.2,85.6,6.9,-69.9,10.6,-70.0,170]
lats = [46.0,28.0,42.8,28.2,45.8,-33.6,46.5,-30.1,-43]
values = [2665,5471,3050,4076,2030,3459,2623,4570,900]
#sizes = [100, 50, 20]

#var_change = temp_change

# Create the projection
fig, ax = plt.subplots(1, 1, figsize=(10,5), subplot_kw={'projection':cartopy.crs.PlateCarree()})
# Add country borders for reference
ax.add_feature(cartopy.feature.BORDERS, alpha=0.3, zorder=1)
ax.add_feature(cartopy.feature.COASTLINE)

# Set the extent
ax.set_extent([east, west, south, north], cartopy.crs.PlateCarree())
# Label title, x, and y axes
ax.set_xticks(np.arange(east,west+1,xtick), cartopy.crs.PlateCarree())
ax.set_yticks(np.arange(south,north+1,ytick), cartopy.crs.PlateCarree())
ax.set_xlabel(xlabel, size=labelsize)
ax.set_ylabel(ylabel, size=labelsize)

# Add regions
#  facecolor='none' just plots the lines
#group_shp = cartopy.io.shapereader.Reader(rgiO1_shp_fn)
#group_feature = cartopy.feature.ShapelyFeature(group_shp.geometries(), cartopy.crs.PlateCarree(),
#                                               edgecolor='black', facecolor='grey', alpha=0.2, linewidth=1)
#ax.add_feature(group_feature,zorder=2)

# Add colorbar
cmap = 'RdYlBu_r'
norm = plt.Normalize(int(np.min(values)), np.ceil(np.max(values)))
var_label = 'Elevation [m a.s.l.]'
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm._A = []
plt.colorbar(sm, ax=ax, fraction=0.024, pad=0.01)
fig.text(1, 0.5, var_label, va='center', ha='center', rotation='vertical', size=labelsize)

ax.scatter(lons, lats, c=values, s=50, cmap=cmap, norm=norm, edgecolors='k', linewidth=0.5, alpha=0.9, zorder=3)

#ax.pcolormesh(lons, lats, var_change, cmap=cmap, norm=norm, zorder=2, alpha=0.8)

# Title
#ax.set_title('[INSERT TITLE]')

# Save figure
fig.set_size_inches(6,4)
figure_fp = os.getcwd() + '/../Output/'
if os.path.exists(figure_fp) == False:
    os.makedirs(figure_fp)
fig_fn = 'example_plot_map.png'
fig.savefig(figure_fp + fig_fn, bbox_inches='tight', dpi=300)

#%%
import xarray as xr

ds1 =


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

from absl import app
from absl import flags
from jaxline import base_config
from jaxline import platform
from examples.autoencoder_mnist import experiment
from ml_collections import config_dict

Experiment = experiment.AutoencoderMnistExperiment

def get_config() -> config_dict.ConfigDict:
  """Creates the config for the experiment."""
  config = base_config.get_base_config()
  config.random_seed = 123109801
  config.training_steps = None
  config.interval_type = None
  config.logging_interval_type = "steps"
  config.log_train_data_interval = 10
  config.log_tensors_interval = 1
  config.checkpoint_interval_type = "steps"
  config.save_checkpoint_interval = 100
  config.checkpoint_dir = "/tmp/kfac_jax_jaxline/"
  config.train_checkpoint_all_hosts = False

  # Experiment config.
  config.experiment_kwargs = config_dict.ConfigDict(
      dict(
          config=dict(
              l2_reg=1e-5,
              training=dict(
                  steps=5_000,
                  epochs=None,
              ),
              batch_size=dict(
                  train=dict(
                      total=60_000,
                      per_device=-1,
                  ),
                  eval=dict(
                      total=10_000,
                      per_device=-1,
                  ),
              ),
              optimizer=dict(
                  name="kfac",
                  kfac=dict(
                      inverse_update_period=5,
                      damping_adaptation_interval=5,
                      num_burnin_steps=5,
                      curvature_ema=0.95,
                      use_adaptive_damping=True,
                      use_adaptive_learning_rate=True,
                      use_adaptive_momentum=True,
                      damping_adaptation_decay=0.95,
                      initial_damping=150.0,
                      min_damping=1e-5,
                      max_damping=1000.0,
                  ),
              )
          )
      )
  )

  config.lock()
  return config

if __name__ == "__main__":
  flags.mark_flag_as_required("config")
  app.run(functools.partial(platform.main, Experiment))

reason: from jax


text:
from pulp import *
import sys
import os
import copy
import numpy as np
sys.path.insert(1, os.path.join(sys.path[0], '..'))
from vaccine_alloc_instance import *
from allocation_solution import *

class LPOnlineMaxFlowMaxUtility:
	def __init__(self, vaccine_obj):
		self.vaccine_obj = copy.deepcopy(vaccine_obj)

	def solve(self):

		lmda=0.5*(1/np.sum(self.vaccine_obj.U_nxd))

		#Storing the quotas locally to update after every iteration
		local_q = self.vaccine_obj.q
		local_Q_d = self.vaccine_obj.Q_d.copy()
		local_Q_c = self.vaccine_obj.Q_c.copy()
		local_Q_cxd = self.vaccine_obj.Q_cxd.copy()
		vaccination_status = [False for i in range(self.vaccine_obj.n)]

		#Create a solution instance to keep track of agents who gets vaccinated
		sol = AllocationSolution(self.vaccine_obj.n,self.vaccine_obj.d)

		for d_i in range(self.vaccine_obj.d):
			catagories = ["c_"+str(i) for i in range(self.vaccine_obj.c)]
			agents = []
			for i in range(self.vaccine_obj.n):
				if (self.vaccine_obj.availability[i][d_i]==1 and vaccination_status[i]==False):
					agents.append("a_"+str(i))
			catagory_nodes = ["c_"+str(d_i)+"_"+str(j) for j in range(self.vaccine_obj.c)]
			nodes = ["s","v","d_"+str(d_i)]+catagory_nodes+agents+["t"]

				# supply or demand of nodes
			            #NodeID : [Supply,Demand]
			nodeData = {x:[0,0] for x in nodes}
			nodeData["s"] = [local_q,0]
			nodeData["t"] = [0,local_q]

			# Arcs List
			arcs = [("s","v"),("v","d_"+str(d_i))]

			for j in range(self.vaccine_obj.c):
				arcs.append(("d_"+str(d_i),"c_"+str(d_i)+"_"+str(j)))

			for i in range(self.vaccine_obj.n):
				if (self.vaccine_obj.availability[i][d_i]==1 and vaccination_status[i]==False):
					for k in range(self.vaccine_obj.c):
						if (self.vaccine_obj.belongsToCatagory[i][k]):
							arcs.append(("c_"+str(d_i)+"_"+str(k),"a_"+str(i)))

			for i in range(self.vaccine_obj.n):
				if(vaccination_status[i]==False):
					arcs.append(("a_"+str(i),"t"))

			#Arcs Data
			# arcs cost, lower bound and capacity
			   		#Arc : [Cost,MinFlow,MaxFlow]
			arcData = {("s","v"):[0,0,local_q]}
			arcData[("v","d_"+str(d_i))]=[0,0,local_Q_d[d_i]]

			for j in range(self.vaccine_obj.c):
				arcData[("d_"+str(d_i),"c_"+str(d_i)+"_"+str(j))]=[0,0,local_Q_cxd[j][d_i]]

			for i in range(self.vaccine_obj.n):
				if (self.vaccine_obj.availability[i][d_i]==1 and vaccination_status[i]==False):
					for k in range(self.vaccine_obj.c):
						if (self.vaccine_obj.belongsToCatagory[i][k]):
							arcData[("c_"+str(d_i)+"_"+str(k),"a_"+str(i))]=[self.vaccine_obj.U_nxd[i][d_i],0,1]

			for i in range(self.vaccine_obj.n):
				if(vaccination_status[i]==False):
					arcData[("a_"+str(i),"t")]=[0,0,1]

			# Splits the dictionaries to be more understandable
			(supply, demand) = splitDict(nodeData)
			(costs, mins, maxs) = splitDict(arcData)

			## BEGIN: Setup LP
			# Creates the boundless Variables as Continues_valued_variables
			vars = LpVariable.dicts("Route",arcs,None,None,cat='Continuous')#LpInteger)

			# Creates the upper and lower bounds on the variables
			for a in arcs:
			    vars[a].bounds(mins[a], maxs[a])

			##LP
			# Creates the 'prob' variable to contain the problem data
			prob = LpProblem(name="Online_maxFlow_max_Utility_Qd_Qdxc",sense=LpMaximize)

			# Creates the objective function
			obj = [vars[arcs[0]]]+[vars[a]* costs[a]*lmda for a in arcs[1:]]
			prob += lpSum(obj), "Total Cost of Transport"

			# Creates all problem constraints - this ensures the amount going into each node is equal to the amount leaving
			for k in nodes[1:-1]: #Except for last node node t and first node s
			    prob += (supply[k]+ lpSum([vars[(i,j)] for (i,j) in arcs if j == k]) ==
			             demand[k]+ lpSum([vars[(i,j)] for (i,j) in arcs if i == k])), \
			            "Flow Conservation in Node %s"%k

			# The problem data is written to an .lp file
			# prob.writeLP("simple_MCFP.lp")

			## END: Setup LP
			# The problem is solved using PuLP's choice of Solver
			prob.solve()

			sol.status = prob.status
			sol.total_flow += value(vars[("s","v")])
			local_q = local_q - value(vars[("s","v")])

			for i in range(self.vaccine_obj.n):
				if (self.vaccine_obj.availability[i][d_i]==1 and vaccination_status[i] == False):
					for k in range(self.vaccine_obj.c):
						if (self.vaccine_obj.belongsToCatagory[i][k]==1):
							if ((value(vars[("c_"+str(d_i)+"_"+str(k),"a_"+str(i))]))==1.0):
								sol.Allocation[i][d_i]=k
								sol.total_utility += self.vaccine_obj.U_nxd[i][d_i]
								local_Q_d[d_i] -= 1
								local_Q_c[k] -= 1
								local_Q_cxd[k][d_i] -= 1
								vaccination_status[i] = True

			# print("Today's vaccine",value(vars[("s","v")]),value(vars[("v","d_"+str(d_i))]))
			# print("Agents",todays_agents)
			# print("MF: Day: ",d_i,"- utility:",todays_util)

		return sol

max_stars_repo_name: severus-tux/vaccine-alloc

max_issues_repo_name: severus-tux/vaccine-alloc

max_forks_repo_name: severus-tux/vaccine-alloc


text:
import pandas as pd
import numpy as np

from tqdm import tqdm_notebook

def cyclical_encoding(series, period):
    features = pd.concat([np.sin((2 * np.pi * series / period)),
                          np.cos((2 * np.pi * series / period))], axis=1)
    features.columns = [f'sin_{series.name}', f'cos_{series.name}']
    return features

def create_time_features(
    fld, keep_frac_only=False, include_additional=False, cyc_encode=False, timezone=None,
):
    """
    Create features out of a series of datetimes.
    Time zones are converted to local time if specified.

    :param fld: series of datetimes
    :param keep_frac_only: maintain only fractional times
    :param include_additional: whether to include additional attributes
    :param cyc_encode: whether to cyclically encode time and date features
    :param timezone: string for the time zone. if passed, times are converted to local

    :return df: dataframe with added date and time columns
    """
    # Convert to a series (in case of index)
    fld = pd.to_datetime(pd.Series(fld))

    # Create a dataframe with index as original times
    df = fld.to_frame().drop(columns=[fld.name])

    # Used for naming the columns
    prefix = fld.name
    prefix += "_"

    # Convert to local time and then remove time zone information
    if timezone:
        fld = fld.dt.tz_convert(timezone).dt.tz_localize(None)
        df["local"] = fld

    # Basic attributes
    attr = ["second", "minute", "hour", "year", "month",
            "week", "day", "dayofweek", "dayofyear"]

    if include_additional:
        # Additional attributes to extract
        attr = attr + [
            "is_month_end",
            "is_month_start",
            "is_quarter_end",
            "is_quarter_start",
            "is_year_end",
            "is_year_start",
            "days_in_month",
        ]

    # iterate through each attribute and add it to the dataframe
    for n in attr:
        df[prefix + n] = getattr(fld.dt, n)

    # Add fractional time of day converting to hours
    df[prefix + "fracday"] = (
        df[prefix + "hour"]
        + df[prefix + "minute"] / 60
        + df[prefix + "second"] / 60 / 60
    ) / 24

    # Add fractional time of week
    df[prefix + "fracweek"] = (
        df[prefix + "dayofweek"] + df[prefix + "fracday"]
    ) / 7

    # Add fractional time of month
    df[prefix + "fracmonth"] = (
        (df[prefix + "day"] - 1) + df[prefix + "fracday"]
    ) / (
        fld.dt.days_in_month
    )  # Use fld days_in_month in case this is not
    # one of the attributes specified

    # Calculate days in year (accounting for leap year rules)
    days_in_year = np.where(
        (df[prefix + "year"] % 4 == 0)
        & ( ( df[prefix + "year"] % 100 != 0) | (df[prefix + "year"] % 400 == 0)),
        366,
        365,
    )

    # Add fractional time of year
    df[prefix + "fracyear"] = (
        (df[prefix + "dayofyear"] - 1) + df[prefix + "fracday"]
    ) / days_in_year

    if cyc_encode:
        df = pd.concat([df, cyclical_encoding(
            df[prefix + 'hour'], 24)], axis=1)
        df = pd.concat([df, cyclical_encoding(
            df[prefix + 'dayofweek'], 6)], axis=1)
        df = pd.concat([df, cyclical_encoding(df[prefix + 'day'], 31)], axis=1)
        df = pd.concat([df, cyclical_encoding(
            df[prefix + 'month'], 12)], axis=1)
        df = pd.concat([df] + [cyclical_encoding(df[c], 1)
                               for c in df if 'frac' in c], axis=1)

    if keep_frac_only:
        df = df.drop(
            [
                prefix + c
                for c in [
                    "second",
                    "minute",
                    "hour",
                    "year",
                    "month",
                    "week",
                    "day",
                    "dayofweek",
                    "dayofyear",
                ]
            ],
            axis=1,
        )

    df = df.set_index(fld).sort_index()

    return df

def monthly_validation(data, model, track=False):
    train_stops = np.unique(data.index[data.index.is_month_end].date)

    X = data.copy()
    y = X.pop('energy')
    weighted_score = 0
    total_possible = 0
    train_points = []
    test_points = []
    scores = []

    for date in train_stops:
        y_train, y_test = y[:date], y[date:]
        X_train, X_test = X[:date], X[date:]

        model.fit(X_train, y_train)
        y_hat = model.predict(X_test)

        test_start, test_end = X_test.index.min().date(), X_test.index.max().date()
        n_days = (test_end - test_start).days
        score = 100 - mape(y_test, y_hat)

        if track:
            print(
                f'Accuracy: {score:.2f}% testing from {test_start} to {test_end} ({n_days} days).')
        weighted_score += score * len(X_test)
        total_possible += 100 * len(X_test)
        train_points.append(len(X_train))
        test_points.append(len(X_test))
        scores.append(score)

    model.fit(X, y)

    importance_df = None
    if hasattr(model, 'feature_importances_'):
        importance_df = pd.DataFrame(
            dict(features=X.columns, importance=model.feature_importances_))
    final_score = weighted_score / total_possible
    results_df = pd.DataFrame(
        dict(train_points=train_points, test_points=test_points, score=scores))
    return dict(results=results_df, importances=importance_df, score=final_score)

def mape(y_true, y_pred):
    return 100 * np.mean(np.abs((y_pred - y_true) / y_true))

def data_reading(filename):
    data = pd.read_csv(filename, parse_dates=['timestamp'])
    data = data.dropna(subset=['energy'])
    freq_counts = data['timestamp'].diff(1).value_counts()
    freq = round(freq_counts.idxmax().total_seconds() / 60)
    data = data.set_index('timestamp').sort_index()
    return data, freq, len(data)

def data_testing(filename, model):
    building_id = filename.split('_')[-1].split('.csv')[0]
    data, freq, dpoints = data_reading(filename)
    results = test_time_features(data, model)
    results['freq'] = freq
    results['dpoints'] = dpoints
    results['building_id'] = building_id
    return results

def test_time_features(data, model):

    data = pd.concat([data, create_time_features(data.index, cyc_encode=True)], axis=1)

    scores = []
    methods = []

    y = data.pop('energy')

    normal_features = ['timestamp_' + t for t in ['hour',
                                                  'dayofweek', 'month', 'dayofyear', 'year']]
    normal_cyc_features = ['sin_' + t for t in normal_features if t not in ['timestamp_dayofyear', 'timestamp_year']
                           ] + ['cos_' + t for t in normal_features if t not in ['timestamp_dayofyear', 'timestamp_year']]

    frac_features = ['timestamp_' +
                     t for t in ['fracday', 'fracweek', 'fracmonth', 'fracyear']]
    frac_cyc_features = ['sin_' + t for t in frac_features] + \
        ['cos_' + t for t in frac_features]

    data_normal = data[normal_features].copy()
    data_normal_cyc = data[normal_cyc_features].copy()
    data_frac = data[frac_features].copy()
    data_frac_cyc = data[frac_cyc_features].copy()

    results = {}
    dataset_names = ['normal', 'normal_cyc', 'frac', 'frac_cyc']

    for dataset, name in zip([data_normal,
                              data_normal_cyc,
                              data_frac,
                              data_frac_cyc],
                             dataset_names):

        to_drop = dataset.columns[(dataset.nunique() == 1)
                                  | (dataset.nunique() == len(dataset))]

        dataset = dataset.drop(columns=to_drop)
        dataset['energy'] = y.copy()
        try:
            data_results = monthly_validation(dataset, model)
            scores.append(data_results['score'])
            methods.append(name)
        except Exception as e:
            print(e, name)

    results = pd.DataFrame(dict(score=scores, method=methods))
    return results

max_stars_repo_name: thimotyb/Data-Analysis

max_issues_repo_name: AbhishekChikara/Data-Analysis

max_forks_repo_name: AbhishekChikara/Data-Analysis


text:
"""
============================================
Hierarchical Merging of Region Boundary RAGs
============================================

This example demonstrates how to perform hierarchical merging on region
boundary Region Adjacency Graphs (RAGs). Region boundary RAGs can be
constructed with the :py:func:`skimage.future.graph.rag_boundary` function.
The regions with the lowest edge weights are successively merged until there
is no edge with weight less than ``thresh``. The hierarchical merging is done
through the :py:func:`skimage.future.graph.merge_hierarchical` function.
For an example of how to construct region boundary based RAGs, see
:any:`plot_rag_boundary`.

"""
#%% Clean Memory
for i in list(globals().keys()):
    if(i[0] != '_'):
        exec('del {}'.format(i))

#%%
from well_plate_project.config import data_dir, reportings_dir

# Dataset paths
raw_data_dir = data_dir / 'raw'
interm_data_dir = raw_data = data_dir / 'intermediate'
describe_data_dir = reportings_dir / 'description'

import numpy as np
import cv2
import matplotlib.pyplot as plt

from skimage import data, segmentation, filters, color
from skimage.future import graph
#from matplotlib import pyplot as plt

def weight_boundary(graph, src, dst, n):
    """
    Handle merging of nodes of a region boundary region adjacency graph.

    This function computes the `"weight"` and the count `"count"`
    attributes of the edge between `n` and the node formed after
    merging `src` and `dst`.

    Parameters
    ----------
    graph : RAG
        The graph under consideration.
    src, dst : int
        The vertices in `graph` to be merged.
    n : int
        A neighbor of `src` or `dst` or both.

    Returns
    -------
    data : dict
        A dictionary with the "weight" and "count" attributes to be
        assigned for the merged node.

    """
    default = {'weight': 0.0, 'count': 0}

    count_src = graph[src].get(n, default)['count']
    count_dst = graph[dst].get(n, default)['count']

    weight_src = graph[src].get(n, default)['weight']
    weight_dst = graph[dst].get(n, default)['weight']

    count = count_src + count_dst
    return {
        'count': count,
        'weight': (count_src * weight_src + count_dst * weight_dst)/count
    }

def merge_boundary(graph, src, dst):
    """Call back called before merging 2 nodes.

    In this case we don't need to do any computation here.
    """
    pass

def segm_rag_merg(img, gray):
    edges = filters.sobel(gray)
    labels = segmentation.slic(img, compactness=30, n_segments=700, start_label=1)
    g = graph.rag_boundary(labels, edges)

    #graph.show_rag(labels, g, img)
    #plt.title('Initial RAG')

    labels2 = graph.merge_hierarchical(labels, g, thresh=0.01, rag_copy=False,
                                    in_place_merge=True,
                                    merge_func=merge_boundary,
                                    weight_func=weight_boundary)

    #graph.show_rag(labels, g, img)
    #plt.title('RAG after hierarchical merging')
    #plt.figure()

    out = color.label2rgb(labels2, img, kind='avg', bg_label=0)

    return out

# %%
def clear_all():
    """Clears all the variables from the workspace of the application."""
    gl = globals().copy()
    for var in gl:
        if var[0] == '_': continue
        if 'func' in str(globals()[var]): continue
        if 'module' in str(globals()[var]): continue

        del globals()[var]

def load_test_file():
    image_file_name = 'a2_a_cropped.jpg'
    from well_plate_project.config import data_dir
    path = data_dir / 'raw' / 'Cropped'
    image_file = path /  image_file_name
    assert image_file.is_file()

    img = cv2.imread(str(image_file))
    plt.imshow(img)
    plt.show()
    return img

if __name__ == "__main__":
    clear_all()
    original_image = load_test_file()

    print("Testing ... ")

    img = cv2.cvtColor(original_image,cv2.COLOR_BGR2RGB)
    gray = cv2.cvtColor(original_image,cv2.COLOR_BGR2GRAY)
    #gray=color.rgb2gray(img)
    out = segm_rag_merg(img, gray)

    print("Plotting... ")
    plt.imshow(out)
    plt.title('Final segmentation')
    plt.show()