text:
import transformers
from sentence_transformers import SentenceTransformer
from sentence_transformers.evaluation import SentenceEvaluator
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader
from torch.utils.tensorboard import SummaryWriter
from torch.optim import Optimizer
from torch import Tensor
import torch.nn.functional as F
from typing import Iterable, Dict, Tuple, Type, Callable
from sklearn.cluster import AgglomerativeClustering, SpectralClustering
import os
from tqdm.autonotebook import trange
import GPUtil
from clearml import Logger

class CustomSentenceTransformer(SentenceTransformer):

    def fit(self,
            train_objectives: Iterable[Tuple[DataLoader, nn.Module]],
            evaluator: SentenceEvaluator = None,
            test_evaluator: SentenceEvaluator = None,
            epochs: int = 1,
            steps_per_epoch = None,
            scheduler: str = 'WarmupLinear',
            warmup_steps: int = 10000,
            optimizer_class: Type[Optimizer] = transformers.AdamW,
            optimizer_params : Dict[str, object]= {'lr': 2e-5},
            weight_decay: float = 0.01,
            evaluation_steps: int = 0,
            output_path: str = None,
            save_best_model: bool = True,
            max_grad_norm: float = 1,
            use_amp: bool = False,
            callback: Callable[[float, int, int], None] = None,
            show_progress_bar: bool = True
            ):
        tensorboard_writer = SummaryWriter('./tensorboard_logs')
        if use_amp:
            from torch.cuda.amp import autocast
            scaler = torch.cuda.amp.GradScaler()

        self.to(self._target_device)
        GPUtil.showUtilization()

        if output_path is not None:
            os.makedirs(output_path, exist_ok=True)

        dataloaders = [dataloader for dataloader, _ in train_objectives]

        # Use smart batching
        for dataloader in dataloaders:
            dataloader.collate_fn = self.smart_batching_collate

        loss_models = [loss for _, loss in train_objectives]
        for loss_model in loss_models:
            loss_model.to(self._target_device)

        self.best_score = -9999999

        if steps_per_epoch is None or steps_per_epoch == 0:
            steps_per_epoch = min([len(dataloader) for dataloader in dataloaders])

        num_train_steps = int(steps_per_epoch * epochs)

        # Prepare optimizers
        optimizers = []
        schedulers = []
        for loss_model in loss_models:
            param_optimizer = list(loss_model.named_parameters())

            no_decay = ['bias', 'LayerNorm.bias', 'LayerNorm.weight']
            optimizer_grouped_parameters = [
                {'params': [p for n, p in param_optimizer if not any(nd in n for nd in no_decay)],
                 'weight_decay': weight_decay},
                {'params': [p for n, p in param_optimizer if any(nd in n for nd in no_decay)], 'weight_decay': 0.0}
            ]

            optimizer = optimizer_class(optimizer_grouped_parameters, **optimizer_params)
            scheduler_obj = self._get_scheduler(optimizer, scheduler=scheduler, warmup_steps=warmup_steps,
                                                t_total=num_train_steps)

            optimizers.append(optimizer)
            schedulers.append(scheduler_obj)

        global_step = 0
        data_iterators = [iter(dataloader) for dataloader in dataloaders]

        num_train_objectives = len(train_objectives)

        skip_scheduler = False

        config = {'epochs': epochs, 'steps_per_epoch': steps_per_epoch}
        for epoch in trange(config.get('epochs'), desc="Epoch", disable=not show_progress_bar):
            training_steps = 0
            running_loss_0 = 0.0

            for loss_model in loss_models:
                loss_model.zero_grad()
                loss_model.train()

            for _ in trange(config.get('steps_per_epoch'), desc="Iteration", smoothing=0.05, disable=not show_progress_bar):
                for train_idx in range(num_train_objectives):
                    loss_model = loss_models[train_idx]
                    optimizer = optimizers[train_idx]
                    scheduler = schedulers[train_idx]
                    data_iterator = data_iterators[train_idx]

                    try:
                        data = next(data_iterator)
                    except StopIteration:
                        data_iterator = iter(dataloaders[train_idx])
                        data_iterators[train_idx] = data_iterator
                        data = next(data_iterator)

                    features, labels = data

                    if use_amp:
                        with autocast():
                            loss_value = loss_model(features, labels)

                        scale_before_step = scaler.get_scale()
                        scaler.scale(loss_value).backward()
                        scaler.unscale_(optimizer)
                        torch.nn.utils.clip_grad_norm_(loss_model.parameters(), max_grad_norm)
                        scaler.step(optimizer)
                        scaler.update()
                        if train_idx == 0:
                            running_loss_0 += loss_value.item()

                        skip_scheduler = scaler.get_scale() != scale_before_step
                    else:
                        loss_value = loss_model(features, labels)
                        if train_idx == 0:
                            running_loss_0 += loss_value.item()
                        loss_value.backward()
                        torch.nn.utils.clip_grad_norm_(loss_model.parameters(), max_grad_norm)
                        optimizer.step()

                    optimizer.zero_grad()

                    if not skip_scheduler:
                        scheduler.step()

                training_steps += 1
                global_step += 1

                if evaluation_steps > 0 and training_steps % evaluation_steps == 0:
                    tensorboard_writer.add_scalar('training_loss', running_loss_0/evaluation_steps, global_step)
                    #logger.report_scalar('Loss', 'training_loss', iteration=global_step, value=running_loss_0/evaluation_steps)
                    running_loss_0 = 0.0
                    #self._eval_during_training(evaluator, output_path, save_best_model, epoch, training_steps, callback)
                    if evaluator is not None:
                        score = evaluator(self, output_path=output_path, epoch=epoch, steps=training_steps)
                        tensorboard_writer.add_scalar('val_ARI', score, global_step)
                        #logger.report_scalar('Training progress', 'val_ARI', iteration=global_step, value=score)
                        if callback is not None:
                            callback(score, epoch, training_steps)
                        if score > self.best_score:
                            self.best_score = score
                            if save_best_model:
                                print('Saving model at: '+output_path)
                                self.save(output_path)
                    for loss_model in loss_models:
                        loss_model.zero_grad()
                        loss_model.train()

            #self._eval_during_training(evaluator, output_path, save_best_model, epoch, -1, callback)
            #tensorboard_writer.add_scalar('training_loss', running_loss_0 / evaluation_steps, global_step)
            #logger.report_scalar('Loss', 'training_loss', iteration=global_step, value=running_loss_0 / evaluation_steps)
            if evaluator is not None:
                score = evaluator(self, output_path=output_path, epoch=epoch, steps=training_steps)
                tensorboard_writer.add_scalar('val_ARI', score, global_step)
                #logger.report_scalar('Training progress', 'val_ARI', iteration=global_step, value=score)
                if callback is not None:
                    callback(score, epoch, training_steps)
                if score > self.best_score:
                    self.best_score = score
                    if save_best_model:
                        self.save(output_path)
            if test_evaluator is not None:
                best_model = SentenceTransformer(output_path)
                device = self.device
                if torch.cuda.is_available():
                    self.to(torch.device('cpu'))
                    best_model.to(device)
                    test_ari = best_model.evaluate(test_evaluator)
                    best_model.to(torch.device('cpu'))
                    self.to(device)
                else:
                    test_ari = best_model.evaluate(test_evaluator)
                tensorboard_writer.add_scalar('test_ARI', test_ari, global_step)
                #logger.report_scalar('Training progress', 'test_ARI', iteration=global_step, value=test_ari)

        if evaluator is None and output_path is not None:  # No evaluator, but output path: save final model version
            self.save(output_path)

def euclid_dist(x):
    dist_mat = torch.norm(x[:, None] - x, dim=2, p=2)
    return dist_mat

def euclid_affinity(x, beta):
    dist_mat = torch.norm(x[:, None] - x, dim=2, p=2)
    affinity_mat = torch.exp(-beta*dist_mat/torch.std(dist_mat, unbiased=False))
    return affinity_mat

def clustering(batch_pairscore_matrix, num_clusters):
    batch_adjacency_matrix = np.zeros(batch_pairscore_matrix.shape)
    num_batch = batch_pairscore_matrix.shape[0]
    clustering_labels = []
    for i in range(num_batch):
        cl = AgglomerativeClustering(n_clusters=num_clusters[i], affinity='precomputed', linkage='average')
        cluster_label = cl.fit_predict(batch_pairscore_matrix[i])
        #clustering_labels.append(torch.from_numpy(cluster_label))
        clustering_labels.append(cluster_label)
        for m in range(cluster_label.shape[0]):
            for n in range(cluster_label.shape[0]):
                if cluster_label[m] == cluster_label[n]:
                    batch_adjacency_matrix[i][m][n] = 1.0
    return batch_adjacency_matrix, clustering_labels

class OptimCluster(torch.autograd.Function):

    @staticmethod
    def forward(ctx, batch_pairscore_matrix, lambda_val, num_clusters):
        ctx.lambda_val = lambda_val
        ctx.num_clusters = num_clusters
        ctx.batch_pairscore_matrix = batch_pairscore_matrix.detach().cpu().numpy()
        ctx.batch_adj_matrix, _ = clustering(ctx.batch_pairscore_matrix, ctx.num_clusters)
        return torch.from_numpy(ctx.batch_adj_matrix).float().to(batch_pairscore_matrix.device)

    @staticmethod
    def backward(ctx, grad_output):
        grad_output_numpy = grad_output.detach().cpu().numpy()
        batch_pairscore_matrix_prime = np.maximum(ctx.batch_pairscore_matrix + ctx.lambda_val * grad_output_numpy, 0.0)
        better_batch_adj_matrix, _ = clustering(batch_pairscore_matrix_prime, ctx.num_clusters)
        gradient = -(ctx.batch_adj_matrix - better_batch_adj_matrix) / ctx.lambda_val
        return torch.from_numpy(gradient.astype(np.float32)).to(grad_output.device), None, None

class BBClusterLossModel(nn.Module):

    def __init__(self, model: SentenceTransformer, device, lambda_val: float, reg_const: float):

        super(BBClusterLossModel, self).__init__()
        self.model = model
        self.lambda_val = lambda_val
        self.reg = reg_const
        self.optim = OptimCluster()
        self.device = device

    def true_adj_mat(self, label):
        n = label.numel()
        adj_mat = torch.zeros((n, n))
        for i in range(n):
            for j in range(n):
                if i == j or label[i] == label[j]:
                    adj_mat[i][j] = 1.0
        return adj_mat

    def forward(self, passage_features: Iterable[Dict[str, Tensor]], labels: Tensor):

        # passage_features will be a list of feature dicts -> [all 1st passages in batch, all 2nd passages in batch, ...]
        # labels shape: batch X maxpsg
        batch_size = labels.shape[0]
        n = labels.shape[1]
        ks = [torch.unique(labels[i]).numel() for i in range(batch_size)]
        true_adjacency_mats = torch.stack([self.true_adj_mat(labels[i]) for i in range(batch_size)]).to(self.device)

        # embeddings shape: batch X maxpsg X emb
        embeddings = torch.stack([self.model(passages)['sentence_embedding'] for passages in passage_features], dim=1)
        embeddings_dist_mats = torch.stack([euclid_dist(embeddings[i]) for i in range(batch_size)])
        #print('embedding dist mats device: '+str(embeddings_dist_mats.device))
        #print('true adj device: '+str(true_adjacency_mats.device))
        mean_similar_dist = (embeddings_dist_mats * true_adjacency_mats).sum() / true_adjacency_mats.sum()
        mean_dissimilar_dist = (embeddings_dist_mats * (1.0 - true_adjacency_mats)).sum() / (
                    1 - true_adjacency_mats).sum()
        adjacency_mats = self.optim.apply(embeddings_dist_mats, self.lambda_val, ks).to(self.device)

        #p = torch.sum(true_adjacency_mats, dim=(1,2)) - n
        #adjacency_wt_mats = torch.stack([(1.0 - true_adjacency_mats[i])*p[i]/(n*(n-1)) + true_adjacency_mats[i]*(1.0-p[i]/(n*(n-1))) for i in range(batch_size)])

        #weighted_err_mats = adjacency_wt_mats * (adjacency_mats * (1.0 - true_adjacency_mats) + (1.0 - adjacency_mats) * true_adjacency_mats)
        weighted_err_mats = adjacency_mats * (1.0 - true_adjacency_mats) + (1.0 - adjacency_mats) * true_adjacency_mats
        weighted_err_mean = weighted_err_mats.mean(dim=0).sum()

        #pprint('Weighted err mat mean: %.5f, mean similar dist: %.5f, mean dissimilar dist: %.5f, reg value: %.5f' %
        #       (weighted_err_mean, mean_similar_dist, mean_dissimilar_dist, 20*(mean_similar_dist/mean_dissimilar_dist)))

        loss = weighted_err_mean + self.reg*(mean_similar_dist - mean_dissimilar_dist)
        #loss = weighted_err_mean
        #pprint('Loss: %.5f' % loss.item())
        #print('Loss: '+str(loss.device))
        return loss

class BBClusterLossReluRegModel(nn.Module):

    def __init__(self, model: SentenceTransformer, device, lambda_val: float, reg_const: float):

        super(BBClusterLossReluRegModel, self).__init__()
        self.model = model
        self.lambda_val = lambda_val
        self.reg = reg_const
        self.margin = 5.0
        self.optim = OptimCluster()
        self.device = device

    def true_adj_mat(self, label):
        n = label.numel()
        adj_mat = torch.zeros((n, n))
        for i in range(n):
            for j in range(n):
                if i == j or label[i] == label[j]:
                    adj_mat[i][j] = 1.0
        return adj_mat

    def forward(self, passage_features: Iterable[Dict[str, Tensor]], labels: Tensor):

        batch_size = labels.shape[0]
        n = labels.shape[1]
        ks = [torch.unique(labels[i]).numel() for i in range(batch_size)]
        true_adjacency_mats = torch.stack([self.true_adj_mat(labels[i]) for i in range(batch_size)]).to(self.device)

        embeddings = torch.stack([self.model(passages)['sentence_embedding'] for passages in passage_features], dim=1)
        embeddings_dist_mats = torch.stack([euclid_dist(embeddings[i]) for i in range(batch_size)])

        mean_similar_dist = (embeddings_dist_mats * true_adjacency_mats).sum() / true_adjacency_mats.sum()
        mean_dissimilar_dist = (embeddings_dist_mats * (1.0 - true_adjacency_mats)).sum() / (
                    1 - true_adjacency_mats).sum()
        adjacency_mats = self.optim.apply(embeddings_dist_mats, self.lambda_val, ks).to(self.device)

        weighted_err_mats = adjacency_mats * (1.0 - true_adjacency_mats) + (1.0 - adjacency_mats) * true_adjacency_mats
        weighted_err_mean = weighted_err_mats.mean(dim=0).sum()

        loss = weighted_err_mean + self.reg*F.relu(mean_similar_dist - mean_dissimilar_dist + self.margin)

        return loss

class DBCLossModel(nn.Module):

    def __init__(self, model: SentenceTransformer, device):

        super(DBCLossModel, self).__init__()
        self.model = model
        self.device = device

    def true_adj_mat(self, label):
        n = label.numel()
        adj_mat = torch.zeros((n, n))
        for i in range(n):
            for j in range(n):
                if i == j or label[i] == label[j]:
                    adj_mat[i][j] = 1.0
        return adj_mat

    def forward(self, passage_features: Iterable[Dict[str, Tensor]], labels: Tensor):

        # passage_features will be a list of feature dicts -> [all 1st passages in batch, all 2nd passages in batch, ...]
        # labels shape: batch X maxpsg
        batch_size = labels.shape[0]
        n = labels.shape[1]
        ks = [torch.unique(labels[i]).numel() for i in range(batch_size)]
        true_adjacency_mats = torch.stack([self.true_adj_mat(labels[i]) for i in range(batch_size)]).to(self.device)

        # embeddings shape: batch X maxpsg X emb
        embeddings = torch.stack([self.model(passages)['sentence_embedding'] for passages in passage_features], dim=1)
        embeddings_dist_mats = torch.stack([euclid_dist(embeddings[i]) for i in range(batch_size)])
        mean_norm_similar_dist_squared = (embeddings_dist_mats * embeddings_dist_mats * true_adjacency_mats).sum() / true_adjacency_mats.sum()

        return mean_norm_similar_dist_squared

class BinaryLoss(nn.Module):

    def __init__(self, model: SentenceTransformer, margin: float = 5):
        super(BinaryLoss, self).__init__()
        self.model = model
        self.margin = margin
        self.distance_metric = lambda x, y: F.pairwise_distance(x, y, p=2)

    def forward(self, sentence_features: Iterable[Dict[str, Tensor]], labels: Tensor):
        vecs = [self.model(sentence_feature)['sentence_embedding'] for sentence_feature in sentence_features]

        vecs1, vecs2 = vecs
        pair_distances = self.distance_metric(vecs1, vecs2)
        distance_pos = pair_distances * labels
        distance_neg = pair_distances * (1-labels)

        losses = F.relu(distance_pos - distance_neg + self.margin)
        return losses.mean()

####################################################
# More experiments required for spectral clustering
# with proper similarity metrics and affinity matrix
####################################################

def spectral_clustering(batch_pairscore_matrix, num_clusters):
    batch_adjacency_matrix = np.zeros(batch_pairscore_matrix.shape)
    num_batch = batch_pairscore_matrix.shape[0]
    clustering_labels = []
    '''
    for i in range(num_batch):
        print('pairscore matrix shape: ' + str(batch_pairscore_matrix[i].shape))
    '''
    for i in range(num_batch):
        cl = SpectralClustering(n_clusters=num_clusters[i], affinity='precomputed')
        cluster_label = cl.fit_predict(batch_pairscore_matrix[i])
        #clustering_labels.append(torch.from_numpy(cluster_label))
        clustering_labels.append(cluster_label)
        for m in range(cluster_label.shape[0]):
            for n in range(cluster_label.shape[0]):
                if cluster_label[m] == cluster_label[n]:
                    batch_adjacency_matrix[i][m][n] = 1.0
    return batch_adjacency_matrix, clustering_labels

class OptimSpectralCluster(torch.autograd.Function):

    @staticmethod
    def forward(ctx, batch_pairscore_matrix, lambda_val, num_clusters):
        ctx.lambda_val = lambda_val
        ctx.num_clusters = num_clusters
        ctx.batch_pairscore_matrix = batch_pairscore_matrix.detach().cpu().numpy()
        ctx.batch_adj_matrix, _ = spectral_clustering(ctx.batch_pairscore_matrix, ctx.num_clusters)
        return torch.from_numpy(ctx.batch_adj_matrix).float().to(batch_pairscore_matrix.device)

    @staticmethod
    def backward(ctx, grad_output):
        grad_output_numpy = grad_output.detach().cpu().numpy()
        batch_pairscore_matrix_prime = np.maximum(ctx.batch_pairscore_matrix + ctx.lambda_val * grad_output_numpy, 0.0)
        better_batch_adj_matrix, _ = spectral_clustering(batch_pairscore_matrix_prime, ctx.num_clusters)
        gradient = -(ctx.batch_adj_matrix - better_batch_adj_matrix) / ctx.lambda_val
        return torch.from_numpy(gradient.astype(np.float32)).to(grad_output.device), None, None

class BBSpectralClusterLossModel(nn.Module):

    def __init__(self, model: SentenceTransformer, device, lambda_val: float, reg_const: float, beta: float):

        super(BBSpectralClusterLossModel, self).__init__()
        self.model = model
        self.lambda_val = lambda_val
        self.reg = reg_const
        self.beta = beta
        self.optim = OptimSpectralCluster()
        self.device = device

    def true_adj_mat(self, label):
        n = label.numel()
        adj_mat = torch.zeros((n, n))
        for i in range(n):
            for j in range(n):
                if i == j or label[i] == label[j]:
                    adj_mat[i][j] = 1.0
        return adj_mat

    def forward(self, passage_features: Iterable[Dict[str, Tensor]], labels: Tensor):

        # passage_features will be a list of feature dicts -> [all 1st passages in batch, all 2nd passages in batch, ...]
        # labels shape: batch X maxpsg
        batch_size = labels.shape[0]
        n = labels.shape[1]
        ks = [torch.unique(labels[i]).numel() for i in range(batch_size)]
        true_adjacency_mats = torch.stack([self.true_adj_mat(labels[i]) for i in range(batch_size)]).to(self.device)

        # embeddings shape: batch X maxpsg X emb
        embeddings = torch.stack([self.model(passages)['sentence_embedding'] for passages in passage_features], dim=1)
        embeddings_affinity_mats = torch.stack([euclid_affinity(embeddings[i], self.beta) for i in range(batch_size)])

        mean_similar_affinity = (embeddings_affinity_mats * true_adjacency_mats).sum() / true_adjacency_mats.sum()
        mean_dissimilar_affinity = (embeddings_affinity_mats * (1.0 - true_adjacency_mats)).sum() / (
                    1 - true_adjacency_mats).sum()
        adjacency_mats = self.optim.apply(embeddings_affinity_mats, self.lambda_val, ks).to(self.device)

        p = torch.sum(true_adjacency_mats, dim=(1,2)) - n
        adjacency_wt_mats = torch.stack([(1.0 - true_adjacency_mats[i])*p[i]/(n*(n-1)) +
                                         true_adjacency_mats[i]*(1.0-p[i]/(n*(n-1))) for i in range(batch_size)])

        weighted_err_mats = adjacency_wt_mats * (adjacency_mats * (1.0 - true_adjacency_mats) + (1.0 - adjacency_mats) * true_adjacency_mats)
        weighted_err_mean = weighted_err_mats.mean(dim=0).sum()

        #pprint('Weighted err mat mean: %.5f, mean similar dist: %.5f, mean dissimilar dist: %.5f, reg value: %.5f' %
        #       (weighted_err_mean, mean_similar_dist, mean_dissimilar_dist, 20*(mean_similar_dist/mean_dissimilar_dist)))

        loss = weighted_err_mean + self.reg*(mean_dissimilar_affinity - mean_similar_affinity)
        #loss = 20*(mean_similar_dist/mean_dissimilar_dist)
        #pprint('Loss: %.5f' % loss.item())
        #print('Loss: '+str(loss.device))
        return loss


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import os

import pytest
import numpy as np

from qtpy import QtCore
from glue.core import roi
from cubeviz.tools.moment_maps import MomentMapsGUI

from ...tests.helpers import (toggle_viewer, select_viewer, left_click,
                      left_button_press, right_button_press, enter_slice_text,
                      assert_all_viewer_indices, assert_slice_text)

DATA_LABELS = ['018.DATA', '018.NOISE']

@pytest.fixture(scope='module')
def moment_maps_gui(cubeviz_layout):
    cl = cubeviz_layout

    mm = MomentMapsGUI(cl._data, cl.session.data_collection, parent=cl)

    return mm

def assert_red_stylesheet(widget):
    assert widget.styleSheet() == "color: rgba(255, 0, 0, 128)"

def test_moment_maps_1(moment_maps_gui, cubeviz_layout):
    # Create GUI
    mm = moment_maps_gui
    mm.display()
    mm.order_combobox.setCurrentIndex(0)
    mm.data_combobox.setCurrentIndex(0)

    # Call calculate function and get result
    mm.calculate_callback()
    moment_component_id = [str(x) for x in cubeviz_layout._data.container_2d.component_ids() if str(x).startswith('018.DATA-moment-1')][0]
    np_result = cubeviz_layout._data.container_2d[moment_component_id]

    # Expected result
    np_data = cubeviz_layout._data[DATA_LABELS[0]]
    import spectral_cube
    cube = spectral_cube.SpectralCube(np_data, wcs=cubeviz_layout._data.coords.wcs)
    order = int(mm.order_combobox.currentText())
    cube_moment = np.asarray(cube.moment(order=order, axis=0))

    assert np.allclose(cube_moment, np_result, rtol=0.01, equal_nan=True)
    assert cubeviz_layout.cube_views[1]._widget.cubeviz_unit.unit == 'm'

def test_moment_maps_2(moment_maps_gui, cubeviz_layout):
    # Create GUI
    mm = moment_maps_gui
    mm.display()
    mm.order_combobox.setCurrentIndex(1)
    mm.data_combobox.setCurrentIndex(0)

    # Call calculate function and get result
    mm.calculate_callback()
    moment_component_id = [str(x) for x in cubeviz_layout._data.container_2d.component_ids() if str(x).startswith('018.DATA-moment-2')][0]
    np_result = cubeviz_layout._data.container_2d[moment_component_id]

    # Expected result
    np_data = cubeviz_layout._data[DATA_LABELS[0]]
    import spectral_cube
    cube = spectral_cube.SpectralCube(np_data, wcs=cubeviz_layout._data.coords.wcs)
    order = int(mm.order_combobox.currentText())
    cube_moment = np.asarray(cube.moment(order=order, axis=0))

    # cube_moment - np_result <= atol + rtol * absolute(np_result)
    assert np.allclose(cube_moment, np_result, rtol=0.01, equal_nan=True)
    assert cubeviz_layout.cube_views[1]._widget.cubeviz_unit.unit == 'm2'

ext: py

max_stars_repo_name: duytnguyendtn/cubeviz

max_issues_repo_name: duytnguyendtn/cubeviz

max_forks_repo_name: duytnguyendtn/cubeviz

reason: import numpy


text:
import numpy as np
import pandas as pd
import pytest
from sklearn.datasets import load_boston, load_breast_cancer

from vivid.env import Settings

@pytest.fixture(scope='function', autouse=True)
def stop_logging():
    before = Settings.LOG_LEVEL
    Settings.LOG_LEVEL = 'WARNING'
    yield
    Settings.LOG_LEVEL = before

@pytest.fixture
def regression_data() -> [pd.DataFrame, np.ndarray]:
    X, y = load_boston(True)
    return pd.DataFrame(X), y

@pytest.fixture
def binary_data() -> [pd.DataFrame, np.ndarray]:
    x, y = load_breast_cancer(True)
    return pd.DataFrame(x), y

max_stars_repo_name: upura/vivid

max_issues_repo_name: upura/vivid

max_forks_repo_name: upura/vivid


text:
from datetime import datetime, timezone

import numpy as np
import xarray as xr
from scipy import optimize
from scipy.ndimage import gaussian_filter1d

import carbonplan_trace.v1.glas_height_metrics as ht
import carbonplan_trace.v1.utils as utils

SPEED_OF_LIGHT = 299792458  # m/s
SECONDS_IN_NANOSECONDS = 10 ** -9

def calculate_derived_variables(data, tiles):
    """
    Calculate derived variables in a xarray dataset containing glas data
    """
    # convert receiving waveform digital bins from 0-543 to corresponding distance from satellite
    data["rec_wf_sample_dist"] = (
        (data.rec_wf_sample_loc + data.rec_wf_response_end_time - data.tx_wf_peak_time)
        * SECONDS_IN_NANOSECONDS
        * SPEED_OF_LIGHT
    ) / 2

    # waveform extent here is calculated based on the signal beginning and end identified in the GLAH14 data
    data = ht.get_height_metric_value(ds=data, metric="wf_extent", recalc=True)
    data = ht.get_dist_metric_value(ds=data, metric="gaussian_fit_dist")

    data['glas_elev'] = data.elevation + data.elevation_correction

    srtm = utils.open_srtm_data(tiles=tiles)
    if srtm is not None:
        srtm_raw = utils.find_matching_records(data=srtm, lats=data.lat, lons=data.lon)
        srtm_elev = srtm_raw.srtm
    else:
        srtm_elev = np.nan

    data['srtm_elev'] = srtm_elev + data.delta_ellipse + data.geoid
    del srtm

    return data

def process_coordinates(ds):
    """
    Process lat/lon to get xy from Landsat images, process time from "seconds since 2000/1/1" to unix/epoch timestamp
    All inputs are xr dataarrays
    """
    # original time format is seconds elapsed since Jan 1 2000 12:00:00 UTC
    d0 = datetime(2000, 1, 1, 12, 0, 0, tzinfo=timezone.utc).timestamp()
    ds['time'] = ds.time + d0

    # get datetime object
    ds['datetime'] = xr.apply_ufunc(
        datetime.fromtimestamp, ds.time.fillna(d0), vectorize=True, dask='parallelized'
    )
    ds['datetime'] = ds.datetime.astype(np.datetime64)

    return ds

def get_mask(ds):
    """
    True in mask = records to use
    """
    # all non nulls in the GLAH14 dataset
    mask = ds.lat.notnull()
    mask.name = 'mask'

    # Harris et al 2021 filtering conditions listed in Supplementary Information
    mask = (
        mask
        & (ds.num_gaussian_peaks >= 2)  # have at least two peaks
        # max amplitude of waveform greater than 2x baseline noise
        & (ds.rec_wf.max(dim='rec_bin') >= (ds.noise_mean * 2))
        # discrepancy bt SRTM and GLAS derived elevation less than 30m
        & (abs(ds.srtm_elev.fillna(ds.glas_elev) - ds.glas_elev) <= 30)
        # signal beginning is less than 70m (otherwise indicates potential inference of signal)
        & (abs(ds.sig_begin_offset) <= 70)
        # signal end is less than 20 and greater than 1m (otherwise indicates sig end is improperly captured)
        & (abs(ds.sig_end_offset) <= 20)
        & (abs(ds.sig_end_offset) >= 1)
    )

    return mask.load()

def filter_large_leading_and_trailing_edge_extent(ds):
    metrics = [
        "leading_edge_extent",
        "trailing_edge_extent",
    ]
    ds = ht.get_all_height_metrics(
        ds=ds,
        metrics=metrics,
    )

    # waveform extent here is calculated based on the signal beginning and end identified in the GLAH14 data
    # and leading/trailing edge extents were calculated using the sum of up to 6 gaussian peaks (modeled waveform)
    mask = (
        # leading edge <= 50% of wf extent, otherwise indicates large differences in canopy height
        (ds.leading_edge_extent <= (ds.wf_extent * 0.5))
        # trailing edge <= 35% of wf extent, otherwise indicates impacts from high slope
        & (ds.trailing_edge_extent <= (ds.wf_extent * 0.35))
    )
    total = ds.dims['unique_index']
    ds = ds.where(mask, drop=True)
    remained = ds.dims['unique_index']
    print(
        f'after edge extent filtering, {remained} valid shots remained out of {total} ({100 - round(remained/total*100, 2)}%) filtered'
    )
    return ds

def gaussian(x, amplitude, mean, stddev):
    return amplitude * np.exp(-1 * ((x - mean) ** 2) / (2 * (stddev ** 2)))

def find_gaussian_fit_sigma(wf, default=3):
    """
    wf is a 1D array, not a vectorized function
    """
    x = np.arange(len(wf))
    y = wf - wf.min()  # optimizer can't deal with negative number
    try:
        popt, _ = optimize.curve_fit(gaussian, x, y, p0=[0.5, 25, default])
        sigma = popt[2]
    except RuntimeError:
        # print("Using default sigma")
        sigma = default

    return sigma

def smooth_wf(rec_wf, tx_wf, verbose=False):
    """
    Find sigma from transmitted waveform, and apply gaussian filter smoothing on the recieved waveform with said sigma
    """
    if np.any(np.isnan(rec_wf)):
        if verbose:
            print('skipping record in smooothing due to nans')
        return rec_wf

    # note: in Farina et al 2018 this value is directly taken from GLAH05 data Data_40HZ/Transmit_Energy/d_sigmaTr
    # re fitting here since the data would be in `ns` unit and requires some translation to be used directly here
    sigma = find_gaussian_fit_sigma(tx_wf)

    return gaussian_filter1d(input=rec_wf, sigma=sigma)

def select_valid_area(bins, wf, signal_begin_dist, signal_end_dist):
    """
    vectorized
    """
    # get mask of valid area
    valid = (bins > signal_begin_dist) & (bins < signal_end_dist)

    # set all invalid area to 0s and clip lower at 0
    wf = wf.where(valid, other=0).clip(min=0)

    return wf

def get_modeled_waveform(ds):
    """
    Sum up the gaussian peaks in GLAH14
    """
    # we have gaussian_fit_dist in meters, gaussian_amp in volts
    # gaussian sigma is in nanoseconds, switch to meters and divide by 2 to get 1 side stdev
    gaussian_sigma_in_m = (ds.gaussian_sigma * SECONDS_IN_NANOSECONDS * SPEED_OF_LIGHT) / 2
    modeled_wf = gaussian(
        x=ds.rec_wf_sample_dist,
        amplitude=ds.gaussian_amp,
        mean=ds.gaussian_fit_dist,
        stddev=gaussian_sigma_in_m,
    ).sum(dim='n_gaussian_peaks')

    # denoise
    modeled_wf = modeled_wf - ds.noise_mean

    # set the energy outside of signal begin/end to 0
    modeled_wf = select_valid_area(
        bins=ds.rec_wf_sample_dist,
        wf=modeled_wf,
        signal_begin_dist=ds.sig_begin_dist,
        signal_end_dist=ds.sig_end_dist,
    )

    dims = ds.rec_wf.dims
    modeled_wf = modeled_wf.transpose(dims[0], dims[1])

    return modeled_wf

def get_smoothed_actual_waveform(ds):
    """
    Smooth and de-noise actual waveform from GLAH01
    """
    # apply gaussian filter to smooth
    processed_wf = xr.apply_ufunc(
        smooth_wf,
        ds.rec_wf,
        ds.tx_wf,
        input_core_dims=[["rec_bin"], ["tx_bin"]],
        output_core_dims=[["rec_bin"]],
        vectorize=True,
        dask="parallelized",
        dask_gufunc_kwargs={"allow_rechunk": 1},
        output_dtypes=['float'],
    )

    # denoise
    processed_wf = processed_wf - ds.noise_mean

    # set the energy outside of signal begin/end to 0
    processed_wf = select_valid_area(
        bins=ds.rec_wf_sample_dist,
        wf=processed_wf,
        signal_begin_dist=ds.sig_begin_dist,
        signal_end_dist=ds.sig_end_dist,
    )

    dims = ds.rec_wf.dims
    processed_wf = processed_wf.transpose(dims[0], dims[1])

    return processed_wf

def preprocess(ds, min_lat, max_lat, min_lon, max_lon):
    # find a list of 10x10 degree tile names covering the bounding box
    # the ancillary data used in preprocess are stored as these 10x10 degree tiles
    tiles = utils.find_tiles_for_bounding_box(min_lat, max_lat, min_lon, max_lon)
    # calculate variables used in the rest of the preprocess
    ds = calculate_derived_variables(data=ds, tiles=tiles)

    # stack the record index and shot number together so we have a ~1D tabular data structure
    ds = ds.stack(unique_index=("record_index", "shot_number"))

    # apply filtering
    mask = get_mask(ds)
    total = ds.dims['unique_index']
    ds = ds.where(mask, drop=True)
    remained = ds.dims['unique_index']
    print(
        f'after initial filtering, {remained} valid shots remained out of {total} ({100 - round(remained/total*100, 2)}%) filtered'
    )

    # smooth and denoise waveform
    # print('before smoothing: current time is ', time.strftime("%H:%M:%S", time.localtime()))
    ds["processed_wf"] = get_smoothed_actual_waveform(ds)
    ds["modeled_wf"] = get_modeled_waveform(ds)
    # leading and trailing edge extent calculations require processed wf (which is time consuming to calculate)
    # thus, filtering is broken into two parts
    ds = filter_large_leading_and_trailing_edge_extent(ds)

    # preprocess the coordinate variables
    # print('before coordinates: current time is ', time.strftime("%H:%M:%S", time.localtime()))
    ds = process_coordinates(ds)

    return ds

max_stars_repo_name: carbonplan/trace

max_issues_repo_name: carbonplan/trace

max_forks_repo_name: carbonplan/trace

reason: import numpy,from scipy


text:
# Copyright 2019-2021 Cambridge Quantum Computing
#
# 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.

from pytket import Circuit, Qubit, Bit
from qermit.spam.full_spam_correction import (  # type: ignore
    gen_full_tomography_spam_circuits_task,
    gen_full_tomography_spam_characterisation_task,
    gen_full_tomography_spam_correction_task,
    gen_get_bit_maps_task,
)
from qermit.spam import (  # type: ignore
    CorrectionMethod,
)
from pytket.extensions.qiskit import AerBackend  # type: ignore
import numpy as np

def test_gen_full_tomography_spam_circuits_task():
    b = AerBackend()
    b._characterisation = dict()
    task = gen_full_tomography_spam_circuits_task(
        b, 5, [[Qubit(0), Qubit(1)], [Qubit(2), Qubit(3)]]
    )
    assert task.n_in_wires == 1
    assert task.n_out_wires == 3

    c0 = Circuit(3).CX(0, 1).X(2).measure_all()
    c1 = Circuit(2).X(0).X(1).measure_all()
    wire = [(c0, 10), (c1, 20)]
    res = task([wire])
    assert len(res) == 3
    assert res[0] == wire
    assert len(res[1]) == len(res[2])
    assert len(res[1][0][0].get_commands()) == 6
    assert len(res[1][1][0].get_commands()) == 8
    assert len(res[1][2][0].get_commands()) == 8
    assert len(res[1][3][0].get_commands()) == 10

def test_full_tomography_spam_characterisation_task_gen():
    b = AerBackend()
    b._characterisation = dict()
    qb_subsets = [[Qubit(0), Qubit(1)], [Qubit(2), Qubit(3)]]

    c0 = Circuit(3).CX(0, 1).X(2).measure_all()
    c1 = Circuit(2).X(0).X(1).measure_all()
    wire = [(c0, 10), (c1, 20)]
    spam_info = gen_full_tomography_spam_circuits_task(b, 5, qb_subsets)([wire])

    spam_circs = [c[0] for c in spam_info[1]]

    handles = b.process_circuits(spam_circs, 5)
    results = b.get_results(handles)

    task = gen_full_tomography_spam_characterisation_task(b, qb_subsets)
    assert task.n_in_wires == 2
    assert task.n_out_wires == 1

    task_res = task([results, spam_info[2]])
    assert task_res == (True,)
    char = b.backend_info.misc["FullCorrelatedSpamCorrection"]
    assert char[0] == qb_subsets
    assert char[1] == {
        Qubit(0): (0, 0),
        Qubit(1): (0, 1),
        Qubit(2): (1, 0),
        Qubit(3): (1, 1),
    }
    assert char[2][0].all() == np.identity(2).all()
    assert char[2][1].all() == np.identity(2).all()

def test_gen_get_bit_maps_task():
    task = gen_get_bit_maps_task()
    assert task.n_in_wires == 1
    assert task.n_out_wires == 2

    c0 = Circuit(3).CX(0, 1).X(2).measure_all()
    c1 = Circuit(2, 2).X(0).Measure(0, 0).X(1).SWAP(0, 1).Measure(0, 1)
    wire = [(c0, 10), (c1, 50)]
    res = task([wire])
    assert len(res) == 2
    assert res[0] == wire
    comp0 = (c0.qubit_to_bit_map, {})
    comp1 = (c1.qubit_to_bit_map, {Bit(0): Qubit(0)})
    assert res[1][0] == comp0
    assert res[1][1] == comp1

def test_gen_full_tomography_spam_correction_task():
    # characterise noiseless matrix
    # use prior experiment
    b = AerBackend()
    b._characterisation = dict()
    qb_subsets = [[Qubit(0), Qubit(1)], [Qubit(2), Qubit(3)]]
    c0 = Circuit(3).CX(0, 1).X(2).measure_all()
    c1 = Circuit(2, 2).X(0).Measure(0, 0).X(1).SWAP(0, 1).Measure(0, 1)
    wire = [(c0, 10), (c1, 20)]
    spam_info = gen_full_tomography_spam_circuits_task(b, 5, qb_subsets)([wire])
    spam_circs = [c[0] for c in spam_info[1]]
    handles = b.process_circuits(spam_circs, 5)
    results = b.get_results(handles)
    # just returns bool
    gen_full_tomography_spam_characterisation_task(b, qb_subsets)(
        [results, spam_info[2]]
    )

    handles1 = b.process_circuits([c0, c1], 20)
    results1 = b.get_results(handles1)
    q_b_maps = [(c0.qubit_to_bit_map, {}), (c1.qubit_to_bit_map, {Bit(0): Qubit(0)})]
    task = gen_full_tomography_spam_correction_task(b, CorrectionMethod.Invert)
    assert task.n_in_wires == 3
    assert task.n_out_wires == 1

    wire = [results1, q_b_maps, True]
    corrected_results = task(wire)[0]
    assert len(corrected_results) == 2
    assert corrected_results[0].get_counts()[(0, 0, 1)] == 20
    assert corrected_results[1].get_counts()[(1, 1)] == 20

if __name__ == "__main__":
    test_gen_full_tomography_spam_circuits_task()
    test_full_tomography_spam_characterisation_task_gen()
    test_gen_full_tomography_spam_correction_task()
    test_gen_get_bit_maps_task()

max_stars_repo_name: CQCL/qermit

max_issues_repo_name: CQCL/qermit

max_forks_repo_name: CQCL/qermit


text:
import numpy as np
from copy import deepcopy

from mushroom_rl.utils.callbacks.callback import Callback
from mushroom_rl.utils.table import EnsembleTable

class CollectQ(Callback):
    """
    This callback can be used to collect the action values in all states at the
    current time step.

    """
    def __init__(self, approximator):
        """
        Constructor.

        Args:
            approximator ([Table, EnsembleTable]): the approximator to use to
                predict the action values.

        """
        self._approximator = approximator

        super().__init__()

    def __call__(self, dataset):
        if isinstance(self._approximator, EnsembleTable):
            qs = list()
            for m in self._approximator.model:
                qs.append(m.table)
            self._data_list.append(deepcopy(np.mean(qs, 0)))
        else:
            self._data_list.append(deepcopy(self._approximator.table))

max_stars_repo_name: PuzeLiu/mushroom-rl

max_issues_repo_name: PuzeLiu/mushroom-rl

max_forks_repo_name: PuzeLiu/mushroom-rl


text:
"""The Brownian motion kernel embeddings."""

# Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved.
# SPDX-License-Identifier: Apache-2.0

from typing import Union

import numpy as np

from ...quadrature.interfaces.standard_kernels import IBrownian, IProductBrownian
from ..kernels import LebesgueEmbedding, QuadratureKernel, QuadratureProductKernel
from ..measures import IntegrationMeasure, LebesgueMeasure

class QuadratureBrownian(QuadratureKernel):
    r"""Base class for a Brownian motion kernel augmented with integrability.

    .. math::
        k(x, x') = \sigma^2 \operatorname{min}(x, x')\quad\text{with}\quad x, x' \geq 0,

    where :math:`\sigma^2` is the ``variance`` property.

    .. note::
        This class is compatible with the standard kernel :class:`IBrownian`.
        Each child of this class implements an embedding w.r.t. a specific integration measure.

    .. seealso::
       * :class:`emukit.quadrature.interfaces.IBrownian`
       * :class:`emukit.quadrature.kernels.QuadratureKernel`

    :param brownian_kernel: The standard EmuKit Brownian motion kernel.
    :param measure: The integration measure.
    :param variable_names: The (variable) name(s) of the integral.

    :raises ValueError: If ``measure`` has wrong dimensionality.
    :raises ValueError: If the reasonable box of the measure allows for negative values.

    """

    def __init__(
        self,
        brownian_kernel: IBrownian,
        measure: IntegrationMeasure,
        variable_names: str,
    ) -> None:

        if measure.input_dim != 1:
            raise ValueError(
                "Integration measure for Brownian motion kernel must be 1-dimensional. Current dimesnion is ({}).".format(
                    measure.input_dim
                )
            )

        super().__init__(kern=brownian_kernel, measure=measure, variable_names=variable_names)

        if any(self.reasonable_box.lower_bounds < 0):
            raise ValueError(
                "The domain defined by the reasonable box seems to allow negative values. "
                "Brownian motion is only defined for positive input values."
            )

    @property
    def variance(self) -> float:
        r"""The scale :math:`\sigma^2` of the kernel."""
        return self.kern.variance

class QuadratureBrownianLebesgueMeasure(QuadratureBrownian, LebesgueEmbedding):
    """A Brownian motion kernel augmented with integrability w.r.t. the standard Lebesgue measure.

    .. seealso::
       * :class:`emukit.quadrature.interfaces.IBrownian`
       * :class:`emukit.quadrature.kernels.QuadratureBrownian`
       * :class:`emukit.quadrature.measures.LebesgueMeasure`

    :param brownian_kernel: The standard EmuKit Brownian motion kernel.
    :param measure: The Lebesgue measure.
    :param variable_names: The (variable) name(s) of the integral.

    """

    def __init__(self, brownian_kernel: IBrownian, measure: LebesgueMeasure, variable_names: str = "") -> None:
        super().__init__(brownian_kernel=brownian_kernel, measure=measure, variable_names=variable_names)

    def qK(self, x2: np.ndarray) -> np.ndarray:
        lb = self.measure.domain.lower_bounds[None, :]
        ub = self.measure.domain.upper_bounds[None, :]
        kernel_mean = ub * x2 - 0.5 * x2**2 - 0.5 * lb**2
        return (self.variance * self.measure.density) * kernel_mean.T

    def qKq(self) -> float:
        lb = self.measure.domain.lower_bounds[None, :]
        ub = self.measure.domain.upper_bounds[None, :]
        qKq = 0.5 * ub * (ub**2 - lb**2) - (ub**3 - lb**3) / 6 - 0.5 * lb**2 * (ub - lb)
        return (self.measure.density**2 * self.variance) * float(qKq)

    def dqK_dx(self, x2: np.ndarray) -> np.ndarray:
        ub = self.measure.domain.upper_bounds[None, :]
        return (self.measure.density * self.variance) * (ub - x2).T

class QuadratureProductBrownian(QuadratureProductKernel):
    r"""Base class for a product Brownian kernel augmented with integrability.

    The kernel is of the form :math:`k(x, x') = \sigma^2 \prod_{i=1}^d k_i(x, x')` where

    .. math::
        k_i(x, x') = \operatorname{min}(x_i-c, x_i'-c)\quad\text{with}\quad x_i, x_i' \geq c,

    :math:`d` is the input dimensionality,
    :math:`\sigma^2` is the ``variance`` property
    and :math:`c` is the ``offset`` property.

    .. note::
        This class is compatible with the standard kernel :class:`IProductBrownian`.
        Each subclass of this class implements an embedding w.r.t. a specific integration measure.

    .. seealso::
       * :class:`emukit.quadrature.interfaces.IProductBrownian`
       * :class:`emukit.quadrature.kernels.QuadratureProductKernel`

    :param brownian_kernel: The standard EmuKit product Brownian kernel.
    :param measure: The integration measure.
    :param variable_names: The (variable) name(s) of the integral.

    :raises ValueError: If the reasonable box of the measure allows for values smaller than the offset :math:`c`.

    """

    def __init__(
        self,
        brownian_kernel: IProductBrownian,
        measure: IntegrationMeasure,
        variable_names: str,
    ) -> None:
        super().__init__(kern=brownian_kernel, measure=measure, variable_names=variable_names)

        lower_bounds_x = self.reasonable_box.lower_bounds
        if any(lower_bounds_x < self.offset):
            raise ValueError(
                f"The domain defined by the reasonable box seems allow to for values smaller than the offset "
                f"({self.offset}). Brownian motion is only defined for input values larger than the offset."
            )

    @property
    def variance(self) -> float:
        r"""The scale :math:`\sigma^2` of the kernel."""
        return self.kern.variance

    @property
    def offset(self) -> float:
        r"""The offset :math:`c` of the kernel."""
        return self.kern.offset

class QuadratureProductBrownianLebesgueMeasure(QuadratureProductBrownian, LebesgueEmbedding):
    """A product Brownian kernel augmented with integrability w.r.t. the standard Lebesgue measure.

    .. seealso::
       * :class:`emukit.quadrature.interfaces.IProductBrownian`
       * :class:`emukit.quadrature.kernels.QuadratureProductBrownian`
       * :class:`emukit.quadrature.measures.LebesgueMeasure`

    :param brownian_kernel: The standard EmuKit product Brownian kernel.
    :param measure: The Lebesgue measure.
    :param variable_names: The (variable) name(s) of the integral.

    """

    def __init__(self, brownian_kernel: IProductBrownian, measure: LebesgueMeasure, variable_names: str = "") -> None:
        super().__init__(brownian_kernel=brownian_kernel, measure=measure, variable_names=variable_names)

    def _scale(self, z: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
        return self.variance * z

    def _get_univariate_parameters(self, dim: int) -> dict:
        return {
            "domain": self.measure.domain.bounds[dim],
            "offset": self.offset,
            "normalize": self.measure.is_normalized,
        }

    def _qK_1d(self, x: np.ndarray, **parameters) -> np.ndarray:
        a, b = parameters["domain"]
        offset = parameters["offset"]
        normalization = 1 / (b - a) if parameters["normalize"] else 1.0
        kernel_mean = b * x - 0.5 * x**2 - 0.5 * a**2
        return (kernel_mean.T - offset * (b - a)) * normalization

    def _qKq_1d(self, **parameters) -> float:
        a, b = parameters["domain"]
        offset = parameters["offset"]
        normalization = 1 / (b - a) if parameters["normalize"] else 1.0
        qKq = 0.5 * b * (b**2 - a**2) - (b**3 - a**3) / 6 - 0.5 * a**2 * (b - a)
        return (float(qKq) - offset * (b - a) ** 2) * normalization**2

    def _dqK_dx_1d(self, x: np.ndarray, **parameters) -> np.ndarray:
        a, b = parameters["domain"]
        normalization = 1 / (b - a) if parameters["normalize"] else 1.0
        return (b - x).T * normalization

max_stars_repo_name: amzn/emukit

max_issues_repo_name: amzn/emukit

max_forks_repo_name: amzn/emukit


text:
from copy import Error
from typing import Dict, List, Set
from common.game import Game
from common.state import State
from common.network_wrapper import NNWrapper
import numpy as np
import math

class NNMCTS():
    def __init__(self, game: Game, nnet: NNWrapper, n_sims: int, cpuct=1) -> None:
        self.game: Game = game
        self.nnet = nnet
        self.n_sims = n_sims
        self.cpuct = cpuct
        self.visited: Set[tuple] = set()
        # self.qsa: Dict[tuple, float] = {}
        self.nsa: Dict[tuple, int] = {}
        self.ns: Dict[tuple, int] = {}
        self.ps: Dict[tuple, np.ndarray] = {}
        self.wsa: Dict[tuple, int] = {}  # Wins
        self.dsa: Dict[tuple, int] = {}  # draws
        self.lsa: Dict[tuple, int] = {}  # losses

    def search(self, state: State) -> np.ndarray:
        return self._search(state,0)

    def _search(self, state: State,depth:int) -> np.ndarray:
        if state.is_game_over():
            wdl = state.game_result()
            return wdl[::-1]

        s = state.to_short()
        if s not in self.visited:
            self.visited.add(s)
            obs = state.to_obs()
            self.ps[s], wdl = self.nnet.predict(obs)
            return wdl[::-1]
        max_u: float
        best_a: int
        max_u,  best_a = -float("inf"), -1
        legal_actions = state.get_legal_actions()
        for a, is_legal in enumerate(legal_actions):
            if not is_legal:
                continue
            if (*s, a) not in self.nsa:
                self.nsa[(*s, a)] = 0
            if s not in self.ns:
                self.ns[s] = 0
            if (*s, a) in self.wsa:
                wsa = self.wsa[(*s, a)]
                lsa = self.lsa[(*s, a)]
                nsa = self.nsa[(*s, a)]
                qsa = (wsa - lsa) / nsa
                u = qsa + self.cpuct * self.ps[s][a] * \
                    math.sqrt(self.ns[s]) / (1 + nsa)
            else:
                u = self.cpuct * \
                    self.ps[s][a] * \
                    math.sqrt(self.ns[s]) / (1 + self.nsa[(*s, a)])

            if u > max_u:
                max_u = u
                best_a = a

        if best_a == -1:
            raise Error(f'Monte carlo tree search gave a non existing best action of {best_a}.')
        a = best_a
        new_state = state.move(a)
        wdl = self._search(new_state,depth=depth+1)
        if (*s, a) in self.wsa:
            self.wsa[(*s, a)] += wdl[0]
            self.dsa[(*s, a)] += wdl[1]
            self.lsa[(*s, a)] += wdl[2]
            self.nsa[(*s, a)] += 1
        else:
            self.wsa[(*s, a)] = wdl[0]
            self.dsa[(*s, a)] = wdl[1]
            self.lsa[(*s, a)] = wdl[2]
            self.nsa[(*s, a)] = 1
        self.ns[s] += 1
        return wdl[::-1]

    def get_probs(self, state: State, temperature: float = 1):
        assert not state.is_game_over()
        for _ in range(self.n_sims):
            self._search(state,0)

        s = state.to_short()

        counts: List[float] = [self.nsa[(*s, a)] if(*s, a)
                               in self.nsa else 0 for a in range(self.game.n_actions)]

        if temperature == 0:
            best_actions = np.array(np.argwhere(
                counts == np.max(counts))).flatten()
            best_action = np.random.choice(best_actions)
            probs: List[float] = [0] * len(counts)
            probs[best_action] = 1
            return probs

        counts = [x**(1/temperature) for x in counts]
        counts_sum = float(sum(counts))
        probs = [x/counts_sum for x in counts]
        return probs

max_stars_repo_name: rusenburn/PBT

max_issues_repo_name: rusenburn/PBT

max_forks_repo_name: rusenburn/PBT


text:
#!/usr/bin/env python3
__author__ = "Shivchander Sudalairaj"
__license__ = "MIT"

'''
model.py: Definition of the izhhikevich neuron

This is a python implementation of the Izhikevich Neuron Model
Credits: Eugene M. Izhikevich (neuron_RS2.m)
'''

import numpy as np

class INeuron:
    def __init__(self):
        self.a = 0.02
        self.b = 0.25
        self.c = -65
        self.d = 6

        self.steps = 1000
        self.tau = 0.25  # step size
        self.timesteps = np.arange(0, self.steps + self.tau, self.tau)  # time steps

        self.v = np.zeros(len(self.timesteps))  # voltage history
        self.v[0] = -64
        self.vt = 30.0  # 30mV threshold
        self.u = np.zeros(len(self.timesteps))  # Recovery history
        self.u[0] = self.b * self.v[0]

        self.T1 = 50  # T1 is the time at which the step input rises
        self.spike_ts = np.zeros(len(self.timesteps))  # spike history
        self.I = np.zeros(len(self.timesteps))         # current history

    def simulate(self, I=None):
        """
        :param I: Current mA, default I = 1mA
        :return: None

        Simulates a single neuron and drives it with I
        """
        if I is None:
            self.I = [1.0] * len(self.timesteps)        # default I = 1.0
        elif isinstance(I, float):
            self.I = [I] * len(self.timesteps)          # user defined I
        else:
            self.I = I                                  # I vector

        for t in range(1, len(self.timesteps)):

            if self.v[t-1] < self.vt:
                # calc membrane potential
                dV = (0.04 * self.v[t - 1] + 5) * self.v[t - 1] + 140 - self.u[t - 1]
                self.v[t] = self.v[t - 1] + (dV + self.I[t - 1]) * self.tau
                # calc recovery variable
                du = self.a * (self.b * self.v[t - 1] - self.u[t - 1])
                self.u[t] = self.u[t - 1] + self.tau * du

            else:  # spike
                self.v[t-1] = self.vt  # reset to spike value
                self.spike_ts[t-1] = 1
                self.v[t] = self.c  # reset membrane voltage
                self.u[t] = self.u[t - 1] + self.d  # reset recovery

    def get_v(self):
        """
        :return: vector which accumulates time series V
        """
        return self.v

    def get_timesteps(self):
        """
        :return: returns the time steps
        """
        return self.timesteps

    def mean_spike_rate(self):
        """
        :return: Mean spike rate

        Calculates the mean spike rate for the corresponding I value
        """
        return sum(self.spike_ts[800:])/800


text:
#!/usr/bin/env python
from netCDF4 import Dataset
from tempfile import TemporaryFile
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib as mpl
from matplotlib.colors import LogNorm
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic
import numpy as np

deg = u'\xb0'

CH_list = ['C30H30', 'C30H50', 'C50H30', 'C50H50']

# Open CNN (1981-2017)
f = open('/home/jhkim/cnn/output/ver4f_18month_SODA/combination.gdat','r')
cnn = np.fromfile(f, dtype=np.float32)

# Open New SINTEX (1984-2017)
f = open('/home/jhkim/data/SINTEX_ensemble2/prediction_em/nino3.4/nino3.4.18month_12.gdat','r')
sin = np.fromfile(f, dtype=np.float32)

# Open observation (GODAS, 1981-2017)
f = Dataset('/home/jhkim/data/GODAS/godas.nino.1981_2017.nc','r')
obs = f.variables['pr'][:,0,0,0]

cnn = cnn/np.std(cnn)
sin = sin/np.std(sin)
obs = obs/np.std(obs)

# open SST & HC map
f = Dataset('/home/jhkim/data/GODAS/godas.sst_t300.1980_2017.144x73.nc','r')
osst = np.mean(f.variables['sst'][:,4:7,:,:],axis=1)
ot300 = np.mean(f.variables['t300'][:,4:7,:,:],axis=1)

chimera = np.zeros(osst.shape, dtype=np.float32)
for i in range(osst.shape[0]):
  chimera[i,:,:] = osst[i,:,:]
  chimera[i,28:45,48:113] = ot300[i,28:45,48:113]

# Compute correlation coefficient (1984-2017)
cor_cnn = np.round(np.corrcoef(obs[3:],cnn[3:])[0,1],2)
cor_sin = np.round(np.corrcoef(obs[3:],sin)[0,1],2)

# Open Heatmap of each case (1981-2016)
heat_each = np.zeros((4,36,18,6), dtype=np.float32)
for i in range(4):
  f = open('/home/jhkim/analysis/heatmap-ver4f_18month_SODA/each_case/'+CH_list[i]+'/heatmap.gdat','r')
  heat_each[i,:,:,:] = np.fromfile(f, dtype=np.float32).reshape(37,18,6)[1:37,:,:]

heat_each = np.swapaxes(heat_each,2,3)
heat_each_mean = np.mean(heat_each, axis=0)

# 0-360 -> 0-380
# chimera: 38x73x144 -> 38x73x152
ext_chimera = np.append(chimera,chimera[:,:,0:8],axis=2)

# heatmap: 36x6x18 -> 36x6x19
ext_heatmap = np.append(heat_each_mean,heat_each_mean[:,:,0:1],axis=2)

# standard deviation (36x6x19 -> 6x19)
std_heatmap = np.std(ext_heatmap,axis=0)

# mean heatmap (36x6x19 -> 6x19)
mean_heatmap = np.mean(ext_heatmap,axis=0)

# significant test
mask_heatmap = np.zeros((36,6,19),dtype=np.float32)
for i in range(36):
  for j in range(6):
    for k in range(19):
      level = abs(ext_heatmap[i,j,k]-mean_heatmap[j,k])/(std_heatmap[j,k]/np.sqrt(37))
      if level > 2.56:
        mask_heatmap[i,j,k] = ext_heatmap[i,j,k]

# Draw Figure
# Figure 3-(a)
plt.subplot(2,1,1)
x = np.arange(1,38)
y = np.arange(4,38)
lines = plt.plot(x, obs, 'black', x, cnn, 'orangered', y, sin, 'dodgerblue')
my_plot = plt.gca()
line0 = my_plot.lines[0]
line1 = my_plot.lines[1]
line2 = my_plot.lines[2]
plt.setp(line0,linewidth=2)
plt.setp(line1,linewidth=0.5, marker='o', markersize=2)
plt.setp(line2,linewidth=0.5, marker='v', markersize=2)
plt.legend(('Observation','CNN(Cor='+str(cor_cnn)+')', 'SINTEX-F(Cor='+str(cor_sin)+')'),loc='upper right', prop={'size':7}, ncol=3)
plt.xlim([0,37])
plt.ylim([-3,3.5])
plt.xticks(np.arange(2,39,2),np.arange(1982, 2019, 2), fontsize=6.5)
plt.tick_params(labelsize=6.,direction='in',length=2,width=0.3,color='black')

plt.yticks(np.arange(-3,3.51,1), fontsize=6.5)
plt.grid(linewidth=0.2, alpha=0.7)
plt.axhline(0,color='black',linewidth=0.5)
plt.title('(a) 18-month lead prediction for DJF Nino3.4', fontsize=8, x=0.268, y=0.967)
plt.xlabel('Year', fontsize=7)
plt.ylabel('DJF Nino3.4', fontsize=7)

# Figure 3-(b)
plt.subplot(2,1,2)

x, y  = np.meshgrid(np.arange(0,380,2.5), np.arange(-91.25,91.25,2.5))

# Contour
CS = plt.contour(x,y,ext_chimera[1996-1980,:,:],np.arange(0.5,1.501,0.5),colors='black',linewidths=0.6)
plt.clabel(CS,fontsize=5, fmt='%1.1f',inline=False)
[txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=0, alpha=0.7)) for txt in CS.labelTexts]
zc = CS.collections[0]
plt.setp(zc, linewidth=0.6)
for line in CS.collections:
  if line.get_linestyle() == [(None, None)]:
    print("Solid Line")
  else:
    line.set_linestyle([(0, (5, 3))])

CS = plt.contour(x,y,ext_chimera[1996-1980,:,:],np.arange(-1.5,0,0.5),colors='black',linewidths=0.6)
plt.clabel(CS,fontsize=5, fmt='%1.1f',inline=False)
[txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=0, alpha=0.7)) for txt in CS.labelTexts]
zc = CS.collections[0]
plt.setp(zc, linewidth=0.6)
for line in CS.collections:
  if line.get_linestyle() == [(None, None)]:
    print("Solid Line")
  else:
    line.set_linestyle([(0, (5, 3))])

# shade
cax = plt.imshow(mask_heatmap[1996-1981,:,:], cmap='RdBu_r',clim=[-20,20], extent=[0,380,60,-55],zorder=1)
plt.gca().invert_yaxis()
map = Basemap(projection='cyl', llcrnrlat=-55,urcrnrlat=59, resolution='c',
              llcrnrlon=20, urcrnrlon=380)
map.drawcoastlines(linewidth=0.2)
map.drawparallels(np.arange( -90., 90.,30.),labels=[1,0,0,0],fontsize=6.5,
                  color='grey', linewidth=0.2)
map.drawmeridians(np.arange(0.,380.,60.),labels=[0,0,0,1],fontsize=6.5,
                  color='grey', linewidth=0.2)
map.fillcontinents(color='silver', zorder=3)
space = '                                                               '
plt.title('(b) MJJ 1996 Heatmap'+space+'[97/98 El Niño Case]',fontsize=8, y=0.962,x=0.5)

x = [  120,   280,   280,   120,   120]
y = [  -16,   -16,    23,    23,   -16]
plt.plot(x,y,'black',zorder=4,linewidth=0.9)

cax = plt.axes([0.08, 0.073, 0.72, 0.013])
cbar = plt.colorbar(cax=cax, orientation='horizontal')
cbar.ax.tick_params(labelsize=6.5,direction='out',length=2,width=0.4,color='black')

plt.tight_layout(h_pad=0,w_pad=-0.6)
plt.subplots_adjust(bottom=0.10, top=0.9, left=0.08, right=0.8)
plt.savefig('/home/jhkim/analysis/fig/figure3.png', dpi=1000)
plt.close()


text:
''' Translate input text with trained model. '''

import torch
import torch.utils.data
import argparse
from tqdm import tqdm
import os,json
from models import get_model
from misc.run import get_loader, run_eval, get_forword_results
from sklearn.decomposition import PCA
from sklearn import manifold
import numpy as np
import os

import torch.nn.functional as F
import matplotlib.pyplot as plt
import json
def plot(mode, jsonfile, path, filename, time, specific):
    info = json.load(open(jsonfile))['id_to_category']
    #info = json.load(open("/home/yangbang/VideoCaptioning/MSRVTT/hidden_topic_mining/tc10s10/20_15/annotation.json"))['id_to_category']
    if mode == 'validate':
        start = 6513
    else:
        start = 7010

    data = np.load(os.path.join(path, filename))
    if time == -1:
        print("======= mean =======")
        x = data[:, 0]
        y = data[:, 1]
    else:
        x = data[:, time, 0]
        y = data[:, time, 1]
    color = []
    idx = []
    for i in range(x.shape[0]):
        category = info[str(i + start)]
        #category = int(info['video' + str(i + start)])
        value = category * 0.05
        if specific != -1 and category != specific:
            continue
        idx.append(i)
        color.append(value)

    plt.scatter(x[idx], y[idx], c=color)
    plt.show()

def plot_several_category(mode, jsonfile, path, filename, time,
    specific=[
        3,
        4,
        9,
        #11,
        #14,
        16,
        17
        ], n_frames=8):

    def cal_dis(x, y, num_category=20, all=False, n_frames=8):
        ave_dis = []
        for s in range(num_category):
            idx = []
            for i in range(x.shape[0]):
                category = info[str((i // n_frames if all else i) + start)]
                #category = int(info['video' + str(i + start)])
                value = category * 0.05
                if s != -1 and category != s:
                    continue
                idx.append(i)
            cluster = torch.stack([torch.from_numpy(x[idx]), torch.from_numpy(y[idx])], dim=1)
            center = cluster.mean(0).unsqueeze(0)
            ave_dis.append(F.pairwise_distance(center, cluster, p=2).view(-1))
        ave_dis = torch.cat(ave_dis, 0).mean(0)
        print(ave_dis)

    info = json.load(open(jsonfile))['id_to_category']
    #info = json.load(open("/home/yangbang/VideoCaptioning/MSRVTT/hidden_topic_mining/tc10s10/20_15/annotation.json"))['id_to_category']
    if mode == 'validate':
        start = 6513
    else:
        start = 7010

    data = np.load(os.path.join(path, filename))
    if time == -1:
        print("======= mean =======")
        #x = data[:, 0]
        #y = data[:, 1]
        x = [data[:2990, 0], data[2990:, 0]]
        y = [data[:2990, 1], data[2990:, 1]]
    elif time == -2:
        print("======= all =======")
        x = data[:, 0]
        y = data[:, 1]
    else:
        x = data[:, time, 0]
        y = data[:, time, 1]
    color = ['r', 'g', 'b', 'purple', 'orange', 'mediumturquoise', 'orange']
    label = ['sports', 'news', 'vehicles', 'food', 'cooking']
    #label = ['sports', 'news', 'vehicles', 'tarvel', 'kids', 'food', 'cooking']

    if isinstance(x, list):
        index = 1
        for item1, item2 in zip(x, y):
            print(item1.shape, item2.shape)
            print(item1[0], item2[0])
            ax = plt.subplot(1, 2, index)
            index += 1

            for j, s in enumerate(specific):
                idx = []
                for i in range(item1.shape[0]):
                    tmp = i // n_frames if time == -2 else i
                    category = info[str(tmp + start)]
                    #category = int(info['video' + str(i + start)])
                    value = category * 0.05
                    if s != -1 and category != s:
                        continue
                    idx.append(i)
                ax.scatter(item1[idx], item2[idx], c=color[j], label=label[j])
            ax.legend()
    else:
        for j, s in enumerate(specific):
            idx = []
            for i in range(x.shape[0]):
                tmp = i // n_frames if time == -2 else i
                category = info[str(tmp + start)]
                #category = int(info['video' + str(i + start)])
                value = category * 0.05
                if s != -1 and category != s:
                    continue
                idx.append(i)
            plt.scatter(x[idx], y[idx], c=color[j], label=label[j])

        cal_dis(x, y, all=(time==-2), n_frames=n_frames)

    plt.legend()
    plt.show()

def visualize(opt, option):
    '''
    plot(opt.em, option['info_json'], opt.pca_path, '%s_%d_%d.npy' % (opt.pca_name, opt.length_id, opt.gate_id),
        -1 if opt.mean else opt.t,
        opt.specific
        )
    '''
    plot_several_category(opt.em, option['info_json'], opt.pca_path, '%s_%d_%d.npy' % (opt.pca_name, opt.length_id, opt.gate_id),
        -1 if opt.mean else (-2 if opt.all else opt.t)
        )

def cal_centers(mode, outputfile, data, jsonfile, num_category=20):

    info = json.load(open(jsonfile))['id_to_category']
    #info = json.load(open("/home/yangbang/VideoCaptioning/MSRVTT/hidden_topic_mining/tc10s10/20_15/annotation.json"))['id_to_category']
    ind = [[] for _ in range(num_category)]

    start = 6513 if mode == 'validate' else 7010

    for i in range(data.shape[0]):
        category = info[str(i + start)]
        #category = int(info['video' + str(i + start)])
        ind[category].append(i)

    file = open(outputfile, 'w')
    n_frames = data.shape[1]
    assert len(data.shape) == 3
    header = ['time']
    for i in range(num_category):
        header.append('%d' % i)
    file.write(','.join(header) + '\n')

    inner_dis = []
    inter_dis = []
    for nf in range(n_frames):
        d = data[:, nf, :]
        inner = []
        centers = []
        for i in range(num_category):
            cluster = d[ind[i]]
            center = cluster.mean(0).unsqueeze(0)
            centers.append(center)
            ave_dis = F.pairwise_distance(center, cluster, p=2).mean(0)
            inner.append(ave_dis.view(1))
        inner_dis.append(torch.cat(inner, dim=0))

        inter = []
        for i in range(num_category):
            this_center = centers[i]
            other_centers = []
            for j in range(num_category):
                if j == i: continue
                other_centers.append(centers[j])
            other_centers = torch.cat(other_centers, dim=0)
            ave_dis = F.pairwise_distance(this_center, other_centers, p=2).mean(0)
            inter.append(ave_dis.view(1))
        inter_dis.append(torch.cat(inter, dim=0))

    inner_dis = torch.stack(inner_dis, dim=0)
    inter_dis = torch.stack(inter_dis, dim=0)

    for t in range(n_frames):
        write_data = ['%d' % t]
        for i in range(num_category):
            write_data.append('%.3f' % inner_dis[t, i])
        file.write(','.join(write_data) + '\n')

    for t in range(n_frames):
        write_data = ['%d' % t]
        for i in range(num_category):
            write_data.append('%.3f' % inter_dis[t, i])
        file.write(','.join(write_data) + '\n')

    file.close()

def main(opt):
    device = torch.device('cuda' if not opt.no_cuda else 'cpu')
    opt_pth = os.path.join(opt.model_path, 'opt_info.json')
    option = json.load(open(opt_pth, 'r'))
    option.update(vars(opt))
    #print(option)

    model = get_model(option)
    checkpoint = torch.load(os.path.join(opt.model_path, 'best', opt.model_name))
    model.load_state_dict(checkpoint['state_dict'])
    model.to(device)

    loader = get_loader(option, mode=opt.em, print_info=False, specific=opt.specific)
    vocab = loader.dataset.get_vocab()

    # rec length predicted results
    rec = {}
    length = len(option['modality']) - sum(option['skip_info'])
    num_gates = 3
    gate_data = [[[] for _ in range(num_gates)] for __ in range(length)]

    opt.pca_name = opt.pca_name + '_%s' % opt.em + ('_mean' if opt.mean else '') + ('_all' if opt.all else '')
    opt.pca_path = os.path.join(opt.pca_path, opt.model_path.split('/')[-2])
    print(opt.pca_path)
    print(opt.model_path.split('/'))
    if not os.path.exists(opt.pca_path):
        os.makedirs(opt.pca_path)

    if opt.plot:
        visualize(opt, option)
    elif opt.oe:
        metric = run_eval(option, model, None, loader, vocab, device,
                        json_path=opt.json_path, json_name=opt.json_name, print_sent=opt.print_sent, no_score=opt.ns, save_videodatainfo=opt.sv)
        print(metric)
    else:
        for data in tqdm(loader, ncols=150, leave=False):
            with torch.no_grad():
                results, _, _ = get_forword_results(option, model, data, device=device, only_data=True)
                gate = results['gate']
                assert len(gate) == length
                assert len(gate[0]) == num_gates - 1
                assert isinstance(results['enc_output'], list)
                assert len(results['enc_output']) == length
                for i in range(length):
                    gate[i].append(results['enc_output'][i])
                    for j in range(num_gates):
                        gate_data[i][j].append(gate[i][j])

        for i in range(length):
            for j in range(num_gates):
                gate_data[i][j] = torch.cat(gate_data[i][j], dim=0) #[len_dataset, n_frames, dim_hidden]
                print(i, j, gate_data[i][j][0, 0, :10].tolist())

                if i == 0:
                    data = torch.cat([gate_data[0][j], torch.cat(gate_data[1][j], dim=0)], dim=0).cpu().numpy()
                else:
                    data = gate_data[i][j].cpu().numpy()
                name = '%s_%d_%d.npy' % (opt.pca_name, i, j)
                if opt.mean:
                    data = data.mean(1)
                    pca = manifold.TSNE(n_components=opt.pca)
                    #pca = PCA(n_components=opt.pca)     #加载PCA算法，设置降维后主成分数目为2
                    #collect = pca.fit_transform(data) #对样本进行降维
                    #print(pca.explained_variance_ratio_)
                    collect = pca.fit_transform(data) #对样本进行降维
                elif opt.all:
                    bsz, seq_len, dim = data.shape
                    data = data.reshape(bsz * seq_len, dim)
                    pca = manifold.TSNE(n_components=opt.pca)
                    collect = pca.fit_transform(data) #对样本进行降维
                else:
                    assert len(data.shape) == 3
                    seq_len = data.shape[1]
                    collect = []
                    for nf in range(seq_len):
                        x = data[:, nf, :]
                        pca = manifold.TSNE(n_components=opt.pca)
                        #pca = PCA(n_components=opt.pca)     #加载PCA算法，设置降维后主成分数目为2
                        reduced_x = pca.fit_transform(x) #对样本进行降维
                        collect.append(reduced_x)
                    collect = np.stack(collect, 1)
                print(name, collect.shape)
                np.save(os.path.join(opt.pca_path, name), collect)

            #print('--------')
            #print(i, resetgate[i].max(2)[0].max(0)[0].tolist())
            #print(i, inputgate[i].max(2)[0].max(0)[0].tolist())
            #print('--------')
            #print(i, resetgate[i].min(2)[0].min(0)[0].tolist())
            #print(i, inputgate[i].min(2)[0].min(0)[0].tolist())

        '''
        metric = run_eval(option, model, None, loader, vocab, device,
                        json_path=opt.json_path, json_name=opt.json_name, print_sent=opt.print_sent, no_score=opt.ns, save_videodatainfo=opt.sv)
        print(metric)
        '''

def test(opt):
    device = torch.device('cuda' if not opt.no_cuda else 'cpu')
    opt_pth = os.path.join(opt.model_path, 'opt_info.json')
    option = json.load(open(opt_pth, 'r'))
    option.update(vars(opt))

    model = get_model(option)
    checkpoint = torch.load(os.path.join(opt.model_path, 'best', opt.model_name))
    model.load_state_dict(checkpoint['state_dict'])
    model.to(device)

    loader = get_loader(option, mode=opt.em, print_info=False, specific=opt.specific)
    vocab = loader.dataset.get_vocab()

    # rec length predicted results
    rec = {}
    length = len(option['modality']) - sum(option['skip_info'])
    num_gates = 3
    gate_data = [[[] for _ in range(num_gates)] for __ in range(length)]

    opt.pca_name = opt.pca_name + '_%s' % opt.em
    if not os.path.exists(opt.pca_path):
        os.makedirs(opt.pca_path)

    for data in tqdm(loader, ncols=150, leave=False):
        with torch.no_grad():
            results, _, _ = get_forword_results(option, model, data, device=device, only_data=True)
            gate = results['gate']
            assert len(gate) == length
            assert len(gate[0]) == num_gates
            for i in range(length):
                for j in range(num_gates):
                    gate_data[i][j].append(gate[i][j])

    for i in range(length):
        for j in range(num_gates):
            gate_data[i][j] = torch.cat(gate_data[i][j], dim=0) #[len_dataset, n_frames, dim_hidden]
            data = gate_data[i][j]
            name = '%s_%d_%d.csv' % (opt.pca_name, i, j)
            outputfile = os.path.join(opt.pca_path, name)
            cal_centers(opt.em, outputfile, data, option['info_json'])

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='translate.py')
    parser.add_argument('-model_path', type=str, default="/home/yangbang/VideoCaptioning/0105save/MSRVTT/GRU_LSTM/ADD0_WA0_EBN1_SS1_WC0_Mi-Im_type2/")

    parser.add_argument('-model_name', default='0028_175313_179067_182774_182189_180426.pth.tar', type=str)
    #"/home/yangbang/VideoCaptioning/0105save/Youtube2Text/GRU_LSTM/ADD0_WA0_EBN1_SS1_WC0_Mi-Im_type2/best/0047_248567_257792_258102_257767_255302.pth.tar"
    #"/home/yangbang/VideoCaptioning/0105save/MSRVTT/GRU_LSTM/ADD0_WA0_EBN1_SS1_WC0_Mi-Im_type2/best/0028_175313_179067_182774_182189_180426.pth.tar"
    parser.add_argument('-beam_size', type=int, default=5,
                        help='Beam size')
    parser.add_argument('-beam_alpha', type=float, default=1.0)
    parser.add_argument('-batch_size', type=int, default=128,
                        help='Batch size')
    parser.add_argument('-topk', type=int, default=1,
                        help="""If verbose is set, will output the n_best
                        decoded sentences""")
    parser.add_argument('-max_len', type=int, default=30)
    parser.add_argument('-no_cuda', action='store_true')
    parser.add_argument('-em', type=str, default='test')
    parser.add_argument('-print_sent', action='store_true')
    parser.add_argument('-func', type=int, default=1)
    parser.add_argument('-json_path', type=str, default='')
    parser.add_argument('-json_name', type=str, default='')
    parser.add_argument('-ns', default=False, action='store_true')
    parser.add_argument('-sv', default=False, action='store_true')
    parser.add_argument('-rgi', '--return_gate_info', default=False, action='store_true')
    parser.add_argument('-specific', default=-1, type=int)
    parser.add_argument('-pca', default=2, type=int)
    parser.add_argument('--pca_path', default='./pca_results')
    parser.add_argument('-pn', '--pca_name', default='')
    parser.add_argument('-mean', default=False, action='store_true')
    parser.add_argument('-all', default=False, action='store_true')

    parser.add_argument('-plot', default=False, action='store_true')
    parser.add_argument('-t', default=0, type=int)
    parser.add_argument('-length_id', default=0, type=int)
    parser.add_argument('-gate_id', default=0, type=int)

    parser.add_argument('-oe', default=False,action='store_true')
    parser.add_argument('-dfi', '--dummy_feats_i', default=False, action='store_true')
    parser.add_argument('-dfm', '--dummy_feats_m', default=False, action='store_true')
    parser.add_argument('-dfa', '--dummy_feats_a', default=False, action='store_true')
    opt = parser.parse_args()

    main(opt)
    #test(opt)

max_stars_repo_name: ybCliff/VideoCaptioning

max_issues_repo_name: ybCliff/VideoCaptioning

max_forks_repo_name: ybCliff/VideoCaptioning


text:
import numpy
import ctypes

DTW=ctypes.CDLL("/usr/local/lib/dtw_python.so")

# void generates_distances(int len_x, int len_y, int dim, float *distances, float *x, float *y)
generates_distances = DTW.generates_distances
generates_distances.restype = None
generates_distances.argtypes = [
    ctypes.c_int32,
    ctypes.c_int32,
    ctypes.c_int32,
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS'),
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS'),
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS')]

def distances_matrix(x,y):

    assert len(x.shape)==len(y.shape), "Inputs must have the same number of columns"

    if len(x.shape)==1:
        dim = ctypes.c_int32(1)
    else:
        dim = ctypes.c_int32(x.shape[1])

    len_x = ctypes.c_int32(x.shape[0])
    len_y = ctypes.c_int32(y.shape[0])

    distances = numpy.empty((x.shape[0],y.shape[0]),dtype=numpy.float32)

    generates_distances(len_x, len_y, dim, distances, x.astype(dtype=numpy.float32), y.astype(dtype=numpy.float32))

    return distances#.reshape((x.shape[0],y.shape[0]))

# void generates_accumulated_cost(int len_x, int len_y, float *accumulated_cost, float *distances)
generates_accumulated_cost = DTW.generates_accumulated_cost
generates_accumulated_cost.restype = None
generates_accumulated_cost.argtypes = [
    ctypes.c_int32,
    ctypes.c_int32,
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS'),
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS')]

def accumulated_cost_matrix(x,y,distances):
    len_x = ctypes.c_int32(x.shape[0])
    len_y = ctypes.c_int32(y.shape[0])
    accumulated_cost = numpy.empty((x.shape[0],y.shape[0]),dtype=numpy.float32)

    generates_accumulated_cost(len_x, len_y, accumulated_cost, distances)

    return accumulated_cost

# float path_cost(int len_x, int len_y, float *accumulated_cost, float *distances)
path_cost = DTW.path_cost
path_cost.restype = ctypes.c_float
path_cost.argtypes = [
    ctypes.c_int32,
    ctypes.c_int32,
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS'),
    numpy.ctypeslib.ndpointer(dtype=numpy.float32, flags='ALIGNED,C_CONTIGUOUS')]

# int returns_path_length()
returns_path_length = DTW.returns_path_length
returns_path_length.restype = ctypes.c_int32
returns_path_length.argtypes = []

# void returns_path(int *output, int size)
returns_path = DTW.returns_path
returns_path.restype = None
returns_path.argtypes = [
    numpy.ctypeslib.ndpointer(dtype=numpy.int32, flags='ALIGNED,C_CONTIGUOUS'),
    ctypes.c_int32]

def path_cost_calculation(x,y,distances,accumulated_cost):
    len_x = ctypes.c_int32(x.shape[0])
    len_y = ctypes.c_int32(y.shape[0])

    cost = path_cost(len_x, len_y, accumulated_cost, distances)

    path_length = returns_path_length()

    path = numpy.empty((path_length,2),dtype=numpy.int32)

    returns_path(path, path_length)

    return path,cost

max_stars_repo_name: meenamishra/DTW

max_issues_repo_name: meenamishra/DTW

max_forks_repo_name: meenamishra/DTW


text:
# Launch with
#
# gunicorn -D --threads 4 -b 0.0.0.0:5000 --access-logfile server.log --timeout 60 server:app glove.6B.300d.txt bbc

from flask import Flask, render_template
import sys,os
import pandas as pd
import numpy as np

app = Flask(__name__)

@app.route("/")
def Users():
    """Show a list of article titles"""

    return render_template('Home.html', fav_songs_genre = user_list)

@app.route("/users/<user_id>")
def article(user_id):
    """
    Show an article with relative path filename. Assumes the BBC structure of
    topic/filename.txt so our URLs follow that.
    """
    return render_template('Songs.html', customized = customized_songs[user_id], lyrics = list(lyrics_songs[user_id]))

page1_data = pd.read_csv('page1_users.csv').sort_values('avg_play_count', ascending=False)
unique_users = list(page1_data.user_id.unique())
user_dct = {user_id:[] for user_id in unique_users}
user_genre = {user_id:'' for user_id in unique_users}

for i in range(page1_data.shape[0]):
    user_dct[page1_data.loc[i,'user_id']].append(page1_data.loc[i,'title'])
    user_genre[page1_data.loc[i,'user_id']] = page1_data.loc[i,'genre']

user_list = []
for i in user_dct:
    user_dct[i] = user_dct[i][:5]
    #user_dct[i].append(user_genre[i])
    user_list.append([i,user_dct[i],user_genre[i]])

page2_data = pd.read_csv('page2.csv')
customized_songs = {i:[] for i in user_dct}

for i in range(page2_data.shape[0]):
    customized_songs[page2_data.loc[i,'user_id']].append([page2_data.loc[i,'title'],
                                                          page2_data.loc[i,'artist_name'],
                                                          page2_data.loc[i,'prediction']])

page2_lyrics = pd.read_csv('page2_lyrics.csv')

lyrics_songs = {i:set() for i in user_dct}

for i in range(page2_lyrics.shape[0]):
    lyrics_songs[page2_lyrics.loc[i,'user_id']].add((page2_lyrics.loc[i,'title'],
                                                          page2_lyrics.loc[i,'artist_name'],
                                                          page2_lyrics.loc[i,'popularity_prediction']))

for i in lyrics_songs:
    if len(lyrics_songs[i])<10:
        n_random = 10 - len(lyrics_songs[i])
        np.random.seed(10)
        list_ind = np.random.randint(low = 1 , high = 3160, size = n_random)
        for j in list_ind:
            lyrics_songs[i].add((page2_lyrics.loc[j,'title'],
                                                          page2_lyrics.loc[j,'artist_name'],
                                                          page2_lyrics.loc[j,'popularity_prediction']))

app.run(debug = True)


text:
import os
import sys
import json
import numpy as np
import threading
import time
import copy
import random
import glob
import shutil

from termcolor import colored
from sacred import Ingredient, Experiment

from alfred.env.thor_env import ThorEnv
from alfred.gen import constants
from alfred.gen.utils import augment_util, video_util
from alfred.utils import helper_util, model_util

args_ingredient = Ingredient('args')
ex = Experiment('render_trajs', ingredients=[args_ingredient])

@args_ingredient.config
def cfg_args():
    # dataset folder to dump frames to
    data_output = 'generated_2.1.0'
    # dataset folder to load jsons from
    data_input = 'json_2.1.0'
    # smooth naviagation (like the original data)
    smooth_nav = True
    # time delays (like the original data)
    time_delays = True
    # whether to shuffle the order of augmenting
    shuffle = False
    # number of threads to start in parallel
    num_threads = 8
    # frame size to render
    render_size = 300
    # X server number
    x_display = '0'
    # render and save RGB images
    render_frames = True
    # render and save depth images
    render_depth = False
    # render and save class segmentation masks
    render_class_masks = False
    # render and save instance segmentation masks
    render_instance_masks = False
    # save object bounding boxes
    save_detections = False
    # partitions to render data for
    partitions = ('train', 'valid_seen', 'valid_unseen')
    # whether to overwrite data folder if it already exists
    overwrite = False

def setup_task(env, traj_data, args):
    # scene setup
    scene_num = traj_data['scene']['scene_num']
    object_poses = traj_data['scene']['object_poses']
    object_toggles = traj_data['scene']['object_toggles']
    dirty_and_empty = traj_data['scene']['dirty_and_empty']
    # reset
    scene_name = 'FloorPlan%d' % scene_num
    env.reset(scene_name, silent=True)
    env.restore_scene(object_poses, object_toggles, dirty_and_empty)
    env.step(dict(traj_data['scene']['init_action']))
    print("Task: %s (%s)" % (traj_data['task_type'], traj_data['task_id']))
    # setup task
    env.set_task(traj_data, reward_type='dense')
    augment_util.check_image(env.last_event.frame)

def augment_traj(env, json_file, args, video_saver, render_settings):
    # load json data
    with open(json_file) as f:
        traj_data = json.load(f)
    # remember images corresponding to low-level actions and create a fresh list
    action_images_orig = [None] * len(traj_data['plan']['low_actions'])
    for image_dict in traj_data['images']:
        if action_images_orig[image_dict['low_idx']] is None:
            action_images_orig[image_dict['low_idx']] = image_dict
    traj_data['images'] = list()

    root_dir_to, rendered_images_dir, save_settings = augment_util.prepare_for_traj(
        json_file, args)
    setup_task(env, traj_data, args)
    rewards, img_count = [], 0

    for ll_idx, ll_action in enumerate(traj_data['plan']['low_actions']):
        # check the allignment of the old json and the rendered images
        alligned_image = action_images_orig[ll_idx]
        if alligned_image['high_idx'] != ll_action['high_idx']:
            print(colored('high_idxs are not alligned', 'red'))
            return False

        # next cmd under the current hl_action
        cmd = ll_action['api_action']
        hl_action = traj_data['plan']['high_pddl'][ll_action['high_idx']]
        # remove unnecessary keys
        cmd = {k: cmd[k] for k in [
            'action', 'objectId', 'receptacleObjectId',
            'placeStationary', 'forceAction'] if k in cmd}

        if ('MoveAhead' in cmd['action'] or
            'Rotate' in cmd['action'] or
            'Look' in cmd['action']):
            event, img_count = augment_util.env_navigate(
                cmd, env, save_settings, root_dir_to,
                render_settings, args.smooth_nav, img_count)
            if event is None:
                return False
        # handle the exception for CoolObject tasks where the actual
        # 'CoolObject' action is actually 'CloseObject'
        elif "CloseObject" in cmd['action'] and \
             "CoolObject" in hl_action['planner_action']['action'] and \
             "OpenObject" in traj_data['plan']['low_actions'][ll_idx + 1][
                 'api_action']['action']:
            cool_action = hl_action['planner_action']
            event, img_count = augment_util.env_interact(
                cmd, env, save_settings, root_dir_to,
                args.time_delays, img_count, action_dummy=cool_action)
        else:
            event, img_count = augment_util.env_interact(
                cmd, env, save_settings, root_dir_to, args.time_delays, img_count)

        # update image list
        img_count_before = len(traj_data['images'])
        for j in range(img_count - img_count_before):
            traj_data['images'].append({
                'low_idx': ll_idx,
                'high_idx': ll_action['high_idx'],
                'image_name': '%09d.png' % int(img_count_before + j)
            })
        if not event.metadata['lastActionSuccess']:
            print(colored("Replay Failed: %s" % (
                env.last_event.metadata['errorMessage']), 'red'))
            return False
        reward, _ = env.get_transition_reward()
        rewards.append(reward)

    # save 1 frame in the end and increase the counter by 10
    # (to be alligned with the train data)
    augment_util.save_image(env.last_event, root_dir_to, save_settings, img_count)
    img_count += 10
    # store color to object type dictionary
    color_to_obj_id_type = {}
    all_objects = env.last_event.metadata['objects']
    for color, object_id in env.last_event.color_to_object_id.items():
        for obj in all_objects:
            if object_id == obj['objectId']:
                color_to_obj_id_type[str(color)] = {
                    'objectID': obj['objectId'],
                    'objectType': obj['objectType']
                }
    augmented_traj_data = copy.deepcopy(traj_data)
    augmented_traj_data['scene']['color_to_object_type'] = color_to_obj_id_type
    augmented_traj_data['task'] = {'rewards': rewards,
                                   'reward_upper_bound': sum(rewards)}
    #FIXME: POSITION CHECK (save distance between final )
    # simulator metadata
    # write an updated traj_data.json (updated images, colors and rewards)
    with open(os.path.join(root_dir_to, 'traj_data.json'), 'w') as aj:
        json.dump(augmented_traj_data, aj, sort_keys=True, indent=4)

    # save video
    images_path = os.path.join(rendered_images_dir, '*.png')
    video_save_path = os.path.join(root_dir_to, 'video.mp4')
    video_saver.save(images_path, video_save_path)
    # write compressed frames to the disk
    augment_util.write_compressed_images(args, root_dir_to)
    return True

def start_worker(worker_id, traj_list, args, lock, processed_files_path):
    '''
    worker loop
    '''
    if isinstance(args.x_display, (list, tuple)):
        x_display = args.x_display[worker_id % len(args.x_display)]
    else:
        x_display = args.x_display
    env = ThorEnv(x_display=x_display,
                  player_screen_width=args.render_size,
                  player_screen_height=args.render_size)
    video_saver = video_util.VideoSaver()
    render_settings = {
        'renderImage': True, # otherwise other images won't be rendered as well
        'renderDepthImage': args.render_depth,
        'renderObjectImage': args.render_instance_masks,
        'renderClassImage': args.render_class_masks}

    while len(traj_list) > 0:
        lock.acquire(timeout=120)
        json_file = traj_list.pop()
        json_path = os.path.join(args.data_input, json_file)
        jsons_left = len(traj_list)
        lock.release()

        print ('Rendering {} ({} left)'.format(json_path, jsons_left))
        augment_success = augment_traj(
            env, json_path, args, video_saver, render_settings)

        # update processed_files on the disk
        lock.acquire(timeout=120)
        with open(processed_files_path, 'a') as f:
            f.write('{};{}'.format(json_file, int(augment_success)) + '\n')
        model_util.update_log(
            args.data_output, stage='augment', update='increase', progress=1)
        lock.release()

    env.stop()
    print("Finished.")

@ex.automain
def main(args):
    args = helper_util.AttrDict(**args)
    if args.data_output is None:
        raise RuntimeError('Please, specify the name of output dataset')
    if (not args.render_frames and not args.render_depth
        and not args.render_instance_masks and not args.render_class_masks):
        raise RuntimeError('At least one type of images should be rendered')

    # set up the paths
    args.data_input = os.path.join(constants.ET_DATA, args.data_input)
    print('Creating a dataset {} using data from {}'.format(
        args.data_output, args.data_input))
    if not os.path.isdir(args.data_input):
        raise RuntimeError('The input dataset {} does not exist'.format(
            args.data_input))
    args.data_output = os.path.join(constants.ET_DATA, args.data_output)
    processed_files_path = os.path.join(args.data_output, 'processed.txt')
    if os.path.exists(args.data_output) and args.overwrite:
            print('Erasing the old directory')
            shutil.rmtree(args.data_output)
    os.makedirs(args.data_output, exist_ok=True)

    # make a list of all the traj_data json files
    traj_list = []
    print('Indexing images in {}'.format(args.partitions))
    for partition in args.partitions:
        for dir_name in sorted(
                glob.glob(os.path.join(args.data_input, partition, '*/*'))):
            if 'trial_' in os.path.basename(dir_name):
                json_path = os.path.join(dir_name, 'traj_data.json')
                if not os.path.isfile(json_path):
                    continue
                traj_list.append('/'.join(json_path.split('/')[-4:]))
    num_files, num_processed_files = len(traj_list), 0

    # remove jsons that were already processed
    if os.path.exists(processed_files_path):
        with open(processed_files_path) as f:
            processed_files = set(
                [line.strip().split(';')[0] for line in f.readlines()])
            # check whether which files are in the desired partitions
            processed_files = set(
                [f for f in processed_files if f.split('/')[0] in args.partitions])
        traj_list = [traj for traj in traj_list if traj not in processed_files]
        num_processed_files += len(processed_files)
    print('{} jsons were already processed'.format(num_processed_files))
    print(colored('The total number of triajectories to process is {}'.format(
        len(traj_list)), 'yellow'))
    model_util.save_log(args.data_output, progress=num_processed_files,
             total=num_files, stage='augment')

    # random shuffle
    if args.shuffle:
        random.shuffle(traj_list)

    lock = threading.Lock()
    if args.num_threads > 0:
        # start threads
        threads = []
        for worker_id in range(min(args.num_threads, len(traj_list))):
            thread = threading.Thread(
                target=start_worker,
                args=(worker_id, traj_list, args, lock, processed_files_path))
            threads.append(thread)
            thread.start()
            time.sleep(1)
        for thread in threads:
            thread.join()
    else:
        # run in the main thread
        start_worker(0, traj_list, args, lock, processed_files_path)
        return

    with open(processed_files_path) as f:
        num_processed_files = len(f.readlines())
    if num_files != num_processed_files:
        print(colored('{} trajectories were skipped'.format(
            num_files - num_processed_files), 'red'))
    else:
        print(colored('All trajectories were successfully recorded', 'green'))

    print('Copying tests folders')
    if not os.path.exists(os.path.join(args.data_output, 'tests_seen')):
        shutil.copytree(os.path.join(args.data_input, 'tests_seen'),
                        os.path.join(args.data_output, 'tests_seen'))
    if not os.path.exists(os.path.join(args.data_output, 'tests_unseen')):
        shutil.copytree(os.path.join(args.data_input, 'tests_unseen'),
                        os.path.join(args.data_output, 'tests_unseen'))

    print('The generated dataset is saved to {}'.format(args.data_output))


text:
"""
1) Generemos un array de 1 dimension con 1000 elementos con distribución normal de media 5 y desvío 2, inicialicemos la semilla en el valor 4703.

2) Usando algunas de las funciones de Numpy listadas en Métodos matemáticos y estadísticos, calculemos la media y el desvío de los elementos del array que construimos en el punto 1.

3) Generemos otro array de dimensiones 100 filas, 20 columnas con distribución normal de media 5 y desvío 2.

4) Usando las mismas funciones que en 2) calculemos la media y el desvío de cada fila del resultado de 3.

5) Usando las mismas funciones que en 2) calculemos la media y el desvío de cada columna del resultado de 3.
"""
import numpy as np
from random_distribuciones import random_binomial

def print_stats(a, axis):
    print(f"Media: {a.mean(axis=axis)}")
    print(f"Desvio: {a.std(axis=axis)}")

array_normal_1 = random_binomial(seed=4703, size=1000, n=5, p=2)
print(array_normal_1[:5])

print_stats(array_normal_1, 0)

array_normal_2 = random_binomial(seed=4703, size=(100, 20), n=5, p=2)
print(array_normal_2[:5, :2])

print_stats(array_normal_2, 0)

print_stats(array_normal_2, 1)