text:
import utils
import processing
import librosa
import config
import numpy as np
import argparse

def detect_emotions(model, song_path):
    """
    Score the given song using one of the trained models
    """
    raw_song = processing.get_raw_data(song_path)
    spectrogram = librosa.feature.melspectrogram(
        y=np.array(raw_song, dtype=np.float), sr=config.SAMPLING_RATE, hop_length=config.FFT_HOP)

    spectrogram = utils.normalise_by_max(spectrogram)
    predictions = utils.average_predictions(spectrogram, model)
    return predictions.numpy()

def display_emotions(scores):
    """
    display score for each emotion
    """
    emotions = config.EMOTIFY_EMOTIONS_ORDERED_LIST
    for i in range(len(emotions)):
        template = "{} : {:1.2f}"
        print(template.format(emotions[i], scores[0][i]))

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("model_path", help="path to model used to predict data")
    parser.add_argument("song_path", help="path to songs to predict scores from")
    args = parser.parse_args()
    model = utils.load_model(args.model_path)
    scores = detect_emotions(model, args.song_path)
    display_emotions(scores)

ext: py

reason: import numpy


text:
import sys, os
import time
import json
import pickle
import argparse
from tqdm import tqdm
import numpy as np
import torch
import torch.nn as nn
from torch.optim import Adam
import torch.utils.data as utils
# from models.SELDNet import Seldnet_vanilla, Seldnet_augmented
from SEDNet import Sednet_vanilla, Sednet_augmented, QSednet_augmented, PHMSednet_augmented, Full_PHMSednet_augmented
from utility_functions import load_model, save_model

import matplotlib.pyplot as plt

'''
Train our baseline model for the Task2 of the L3DAS21 challenge.
This script saves the best model checkpoint, as well as a dict containing
the results (loss and history). To evaluate the performance of the trained model
according to the challenge metrics, please use evaluate_baseline_task2.py.
Command line arguments define the model parameters, the dataset to use and
where to save the obtained results.
'''

def evaluate(model, device, criterion_sed, criterion_doa, dataloader):
    #compute loss without backprop
    model.eval()
    test_loss = 0.
    with tqdm(total=len(dataloader) // args.batch_size) as pbar, torch.no_grad():
        for example_num, (x, target) in enumerate(dataloader):
            target = target.to(device)
            x = x.to(device)
            t = time.time()
            # Compute loss for each instrument/model
            sed = model(x)
            loss = seld_loss(x, target, model, criterion_sed, criterion_doa)
            test_loss += (1. / float(example_num + 1)) * (loss - test_loss)
            pbar.set_description("Current loss: {:.4f}".format(test_loss))
            pbar.update(1)
    return test_loss

def seld_loss(x, target, model, criterion_sed, criterion_doa):
    '''
    compute seld loss as weighted sum of sed (BCE) and doa (MSE) losses
    '''
    target_sed = target[:,:,:args.output_classes*args.class_overlaps]

    #compute loss
    sed = model(x)
    sed = torch.flatten(sed, start_dim=1)
    target_sed = torch.flatten(target_sed, start_dim=1)

    loss_sed = criterion_sed(sed, target_sed) * args.sed_loss_weight

    return loss_sed

def main(args):
    if args.use_cuda:
        device = 'cuda:' + str(args.gpu_id)
    else:
        device = 'cpu'

    if args.fixed_seed:
        seed = 1
        np.random.seed(seed)
        torch.backends.cudnn.deterministic = True
        torch.backends.cudnn.benchmark = False
        torch.manual_seed(seed)
        if torch.cuda.is_available():
            torch.cuda.manual_seed_all(seed)

    #LOAD DATASET
    print ('\nLoading dataset')

    with open(args.training_predictors_path, 'rb') as f:
        training_predictors = pickle.load(f)
    with open(args.training_target_path, 'rb') as f:
        training_target = pickle.load(f)
    with open(args.validation_predictors_path, 'rb') as f:
        validation_predictors = pickle.load(f)
    with open(args.validation_target_path, 'rb') as f:
        validation_target = pickle.load(f)
    with open(args.test_predictors_path, 'rb') as f:
        test_predictors = pickle.load(f)
    with open(args.test_target_path, 'rb') as f:
        test_target = pickle.load(f)

    training_predictors = np.array(training_predictors)
    training_target = np.array(training_target)
    validation_predictors = np.array(validation_predictors)
    validation_target = np.array(validation_target)
    test_predictors = np.array(test_predictors)
    test_target = np.array(test_target)

    print ('\nShapes:')
    print ('Training predictors: ', training_predictors.shape)
    print ('Validation predictors: ', validation_predictors.shape)
    print ('Test predictors: ', test_predictors.shape)
    print ('Training target: ', training_target.shape)
    print ('Validation target: ', validation_target.shape)
    print ('Test target: ', test_target.shape)

    # Normalize training predictors with mean 0 and std 1
    train_mag_min = np.mean(training_predictors[:,:8,:,:])
    train_mag_std = np.std(training_predictors[:,:8,:,:])
    train_phase_min = np.mean(training_predictors[:,8:,:,:])
    train_phase_std = np.std(training_predictors[:,8:,:,:])

    # Input normalization
    training_predictors[:,:8,:,:] -= train_mag_min
    training_predictors[:,:8,:,:] /= train_mag_std
    training_predictors[:,8:,:,:] -= train_phase_min
    training_predictors[:,8:,:,:] /= train_phase_std

    # Normalize validation predictors with mean 0 and std 1
    val_mag_min = np.mean(validation_predictors[:,:8,:,:])
    val_mag_std = np.std(validation_predictors[:,:8,:,:])
    val_phase_min = np.mean(validation_predictors[:,8:,:,:])
    val_phase_std = np.std(validation_predictors[:,8:,:,:])
    validation_predictors[:,:8,:,:] -= val_mag_min
    validation_predictors[:,:8,:,:] /= val_mag_std
    validation_predictors[:,8:,:,:] -= val_phase_min
    validation_predictors[:,8:,:,:] /= val_phase_std

    # Normalize test predictors with mean 0 and std 1
    test_mag_min = np.mean(test_predictors[:,:8,:,:])
    test_mag_std = np.std(test_predictors[:,:8,:,:])
    test_phase_min = np.mean(test_predictors[:,8:,:,:])
    test_phase_std = np.std(test_predictors[:,8:,:,:])
    test_predictors[:,:8,:,:] -= test_mag_min
    test_predictors[:,:8,:,:] /= test_mag_std
    test_predictors[:,8:,:,:] -= test_phase_min
    test_predictors[:,8:,:,:] /= test_phase_std

    features_dim = int(test_target.shape[-2] * test_target.shape[-1])

    training_predictors = torch.tensor(training_predictors).float()
    validation_predictors = torch.tensor(validation_predictors).float()
    test_predictors = torch.tensor(test_predictors).float()
    training_target = torch.tensor(training_target).float()
    validation_target = torch.tensor(validation_target).float()
    test_target = torch.tensor(test_target).float()
    #build dataset from tensors
    tr_dataset = utils.TensorDataset(training_predictors, training_target)
    val_dataset = utils.TensorDataset(validation_predictors, validation_target)
    test_dataset = utils.TensorDataset(test_predictors, test_target)

    #build data loader from dataset
    tr_data = utils.DataLoader(tr_dataset, args.batch_size, shuffle=True, pin_memory=True)
    val_data = utils.DataLoader(val_dataset, args.batch_size, shuffle=False, pin_memory=True)
    test_data = utils.DataLoader(test_dataset, args.batch_size, shuffle=False, pin_memory=True)

    #LOAD MODEL
    if args.architecture == 'sednet_vanilla':
        n_time_frames = test_predictors.shape[-1]
        model = Sednet_vanilla(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    n_cnn_filters=args.n_cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose)
    if args.architecture == 'sednet_augmented':
        n_time_frames = test_predictors.shape[-1]
        model = Sednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose)
    if args.architecture == 'qsednet':
        n_time_frames = test_predictors.shape[-1]
        model = QSednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose)
    if args.architecture == 'phmsednet_n2':
        n_time_frames = test_predictors.shape[-1]
        model = PHMSednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose, n=2)
    if args.architecture == 'phmsednet_n4':
        n_time_frames = test_predictors.shape[-1]
        model = PHMSednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose, n=4)
    if args.architecture == 'phmsednet_n8':
        n_time_frames = test_predictors.shape[-1]
        model = PHMSednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose, n=8)
    if args.architecture == 'phmsednet_n16':
        n_time_frames = test_predictors.shape[-1]
        model = PHMSednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose, n=16)
    if args.architecture == 'full_phmsednet_n4':
        n_time_frames = test_predictors.shape[-1]
        model = Full_PHMSednet_augmented(time_dim=n_time_frames, freq_dim=args.freq_dim, input_channels=args.input_channels,
                    output_classes=args.output_classes, pool_size=args.pool_size,
                    pool_time=args.pool_time, rnn_size=args.rnn_size, n_rnn=args.n_rnn,
                    fc_size=args.fc_size, dropout_perc=args.dropout_perc,
                    cnn_filters=args.cnn_filters, class_overlaps=args.class_overlaps,
                    verbose=args.verbose, n=4)

    if args.use_cuda:
        print("Moving model to gpu")
    model = model.to(device)

    #compute number of parameters
    model_params = sum([np.prod(p.size()) for p in model.parameters()])
    print ('Total paramters: ' + str(model_params))

    #set up the loss functions
    criterion_sed = nn.BCELoss()
    criterion_doa = nn.MSELoss()

    #set up optimizer
    optimizer = Adam(params=model.parameters(), lr=args.lr)

    #set up training state dict that will also be saved into checkpoints
    state = {"step" : 0,
             "worse_epochs" : 0,
             "epochs" : 0,
             "best_loss" : np.Inf}

    #load model checkpoint if desired
    if args.load_model is not None:
        print("Continuing training full model from checkpoint " + str(args.load_model))
        state = load_model(model, optimizer, args.load_model, args.use_cuda)

    #TRAIN MODEL
    print('TRAINING START')
    train_loss_hist = []
    val_loss_hist = []
    epoch = 1
#     while state["worse_epochs"] < args.patience:
    for curr_epoch in range(args.epochs):
#         print("Training epoch " + str(epoch))
        print("Training epoch " + str(curr_epoch))
        avg_time = 0.
        model.train()
        train_loss = 0.
        with tqdm(total=len(tr_dataset) // args.batch_size) as pbar:
            for example_num, (x, target) in enumerate(tr_data):

                target = target.to(device)
                x = x.to(device)
                t = time.time()
                # Compute loss for each instrument/model
                optimizer.zero_grad()
                # sed, doa = model(x)
                sed = model(x)

                # loss = seld_loss(x, target, model, criterion_sed, criterion_doa)

                target_sed = target[:,:,:args.output_classes*args.class_overlaps]
                sed_flat = torch.flatten(sed, start_dim=1)
                target_sed_flat = torch.flatten(target_sed, start_dim=1)
                loss = criterion_sed(sed_flat, target_sed_flat)

                loss.backward()

                train_loss += (1. / float(example_num + 1)) * (loss - train_loss)
                optimizer.step()
                state["step"] += 1
                t = time.time() - t
                avg_time += (1. / float(example_num + 1)) * (t - avg_time)

                pbar.update(1)

            #PASS VALIDATION DATA
            val_loss = evaluate(model, device, criterion_sed, criterion_doa, val_data)
            print("VALIDATION FINISHED: LOSS: " + str(val_loss))

            # EARLY STOPPING CHECK
            checkpoint_path = args.checkpoint_dir + "/checkpoint_" + args.architecture

            if val_loss >= state["best_loss"]:
                state["worse_epochs"] += 1
            else:
                print("MODEL IMPROVED ON VALIDATION SET!")
                state["worse_epochs"] = 0
                state["best_loss"] = val_loss
                state["best_checkpoint"] = checkpoint_path

                # CHECKPOINT
                print("Saving model...")
                save_model(model, optimizer, state, checkpoint_path)

            state["epochs"] += 1
            #state["worse_epochs"] = 200
            train_loss_hist.append(train_loss.cpu().detach().numpy())
            val_loss_hist.append(val_loss.cpu().detach().numpy())
            epoch += 1

    #LOAD BEST MODEL AND COMPUTE LOSS FOR ALL SETS
    print("TESTING")
    # Load best model based on validation loss
    state = load_model(model, None, state["best_checkpoint"], args.use_cuda)
    #compute loss on all set_output_size
    train_loss = evaluate(model, device, criterion_sed, criterion_doa, tr_data)
    val_loss = evaluate(model, device, criterion_sed, criterion_doa, val_data)
    test_loss = evaluate(model, device, criterion_sed, criterion_doa, test_data)

    #PRINT AND SAVE RESULTS
    results = {'train_loss': train_loss.cpu().detach().numpy(),
               'val_loss': val_loss.cpu().detach().numpy(),
               'test_loss': test_loss.cpu().detach().numpy(),
               'train_loss_hist': train_loss_hist,
               'val_loss_hist': val_loss_hist}

    print ('RESULTS')
    for i in results:
        if 'hist' not in i:
            print (i, results[i])
    out_path = os.path.join(args.results_path, 'results_dict.json')
    np.save(out_path, results)

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    #saving/loading parameters
    parser.add_argument('--results_path', type=str, default='RESULTS/Task2',
                        help='Folder to write results dicts into')
    parser.add_argument('--checkpoint_dir', type=str, default='RESULTS/Task2',
                        help='Folder to write checkpoints into')
    parser.add_argument('--load_model', type=str, default=None,
                        help='Reload a previously trained model (whole task model)')
    #dataset parameters
    parser.add_argument('--training_predictors_path', type=str, default='DATASETS/processed/task2_predictors_train.pkl')
    parser.add_argument('--training_target_path', type=str, default='DATASETS/processed/task2_target_train.pkl')
    parser.add_argument('--validation_predictors_path', type=str, default='DATASETS/processed/task2_predictors_validation.pkl')
    parser.add_argument('--validation_target_path', type=str, default='DATASETS/processed/task2_target_validation.pkl')
    parser.add_argument('--test_predictors_path', type=str, default='DATASETS/processed/task2_predictors_test.pkl')
    parser.add_argument('--test_target_path', type=str, default='DATASETS/processed/task2_target_test.pkl')
    #training parameters
    parser.add_argument('--gpu_id', type=int, default=0)
    parser.add_argument('--use_cuda', type=str, default='True')
    parser.add_argument('--early_stopping', type=str, default='True')
    parser.add_argument('--fixed_seed', type=str, default='False')

    parser.add_argument('--epochs', type=int, default=200)
    parser.add_argument('--lr', type=float, default=0.00001)
    parser.add_argument('--batch_size', type=int, default=3,
                        help="Batch size")
    parser.add_argument('--sr', type=int, default=32000,
                        help="Sampling rate")
    parser.add_argument('--patience', type=int, default=100,
                        help="Patience for early stopping on validation set")
    #model parameters
    #the following parameters produce a prediction for each 100-msecs frame
    parser.add_argument('--architecture', type=str, default='seldnet_augmented',
                        help="model's architecture, [seldnet_vanilla, seldnet_augmented, qseldnet, phmseldnet_n2, phmseldnet_n4, phmseldnet_n8]")
    parser.add_argument('--input_channels', type=int, default=8,
                        help="4/8 for 1/2 mics, multiply x2 if using also phase information")
    parser.add_argument('--class_overlaps', type=int, default=3,
                        help= 'max number of simultaneous sounds of the same class')

    parser.add_argument('--time_dim', type=int, default=4800)
    parser.add_argument('--freq_dim', type=int, default=256)
    parser.add_argument('--output_classes', type=int, default=14)
    parser.add_argument('--pool_size', type=str, default='[[8,2],[8,2],[2,2],[1,1]]')
    parser.add_argument('--cnn_filters', type=str, default='[64,128,256,512]',
                        help= 'only for seldnet augmented')
    parser.add_argument('--pool_time', type=str, default='True')
    parser.add_argument('--rnn_size', type=int, default=256)
    parser.add_argument('--n_rnn', type=int, default=3)
    parser.add_argument('--fc_size', type=int, default=1024)
    parser.add_argument('--dropout_perc', type=float, default=0.3)
    parser.add_argument('--n_cnn_filters', type=float, default=64,
                        help= 'only for seldnet vanilla')
    parser.add_argument('--verbose', type=str, default='False')
    parser.add_argument('--sed_loss_weight', type=float, default=1.)
    parser.add_argument('--doa_loss_weight', type=float, default=5.)

    args = parser.parse_args()

    #eval string bools and lists
    args.use_cuda = eval(args.use_cuda)
    args.early_stopping = eval(args.early_stopping)
    args.fixed_seed = eval(args.fixed_seed)
    args.pool_size= eval(args.pool_size)
    args.pool_time = eval(args.pool_time)
    args.cnn_filters = eval(args.cnn_filters)
    args.verbose = eval(args.verbose)

    main(args)


text:
import pandas as pd
import numpy as np
import ast
import scipy.sparse as sp
from scipy.spatial.distance import pdist, squareform
import pdb
import pygsp as pg
from scipy import sparse
import scipy.sparse.linalg
from matplotlib import pyplot as plt
from pyunlocbox import functions, solvers
import networkx as nx
from error import error_func

class Our_Graph():
    def __init__(self, adjacency):
        self.adjacency = adjacency
        self.n_nodes = len(adjacency)
        self.n_edges = np.count_nonzero(adjacency)//2
        self.D = np.diag(adjacency.sum(axis=1))
        self.D_norm = np.power(np.linalg.inv(self.D),0.5)
        self.laplacian_combinatorial = self.D - self.adjacency # combinatorial laplacian
        self.laplacian_normalized = self.D_norm @ self.laplacian_combinatorial @ self.D_norm
        #self.compute_gradient(normalized=False)
        self.compute_fourier_basis(normalized=False)

    def get_laplacian(self, normalized=False):
        self.laplacian = self.laplacian_combinatorial  # combinatorial laplacian
        if normalized:
            self.laplacian = self.laplacian_normalized
        return self.laplacian

    def compute_gradient(self, normalized=False):
        # Find incidence matrix. Since our graph is undirected, we chose to only consider the upper right
        # triangle of our adjacency matrix when finding the gradient.
        self.S = np.zeros((self.n_nodes, self.n_edges))
        self.S_normalized = np.zeros((self.n_nodes, self.n_edges))
        edge_idx = 0
        for i in range(self.n_nodes):
            for k in range(i):
                if self.adjacency[i, k] != 0:
                    self.S[i, edge_idx] = np.sqrt(self.adjacency[i,k])
                    self.S[k, edge_idx] = -np.sqrt(self.adjacency[i,k])
                    self.S_normalized[i, edge_idx] = np.sqrt(self.adjacency[i,k]) * self.D_norm[i,i]
                    self.S_normalized[k, edge_idx] = -np.sqrt(self.adjacency[i,k]) * self.D_norm[k,k]
                    edge_idx += 1

        assert np.allclose(self.S @ self.S.T, self.laplacian_combinatorial), "Wrong incidence matrix"
        assert np.allclose(self.S_normalized @ self.S_normalized.T, self.laplacian_normalized), "Wrong normalized incidence matrix"

        if normalized:
            return self.S
        else:
            return self.S_normalized

    def compute_fourier_basis(self, normalized=False):
        # e: Ordered Laplacian eigenvalues
        # U: Ordered graph Fourier basis
        self.e, self.U = np.linalg.eigh(self.laplacian_combinatorial)
        self.e_normalized, self.U_normalized = np.linalg.eigh(self.laplacian_normalized)
        if normalized:
            return self.e_normalized, self.U_normalized
        else:
            return self.e, self.U

    def get_laplacian_eigenmaps(self, new_dim=10, use_normalized=True):
        if use_normalized:
            laplacian = self.laplacian_normalized
        else:
            laplacian = self.laplacian_combinatorial
        eigenvalues, eigenvectors = sparse.linalg.eigsh(laplacian, which='SM', k=new_dim + 1)
        sort_ind = np.argsort(eigenvalues)
        eigenvalues = eigenvalues[sort_ind]
        eigenvectors = eigenvectors[:,sort_ind]
        embeddings = eigenvectors[:, 1:new_dim + 1]
        if use_normalized:
            return (embeddings.T * np.diag(self.D_norm)).T
        return embeddings

    def spectral_clustering(self, test_pos, labels_bin):
        fiedler_vector_c = self.get_laplacian_eigenmaps(new_dim=1, use_normalized=False).squeeze()
        fiedler_vector_n = self.get_laplacian_eigenmaps(new_dim=1, use_normalized=True).squeeze()

        # We do clustering according to the fiedler vector
        clusters_c = np.zeros([self.n_nodes, ])
        clusters_c[fiedler_vector_c > 0] = 1
        clusters_c[fiedler_vector_c <= 0] = -1

        clusters_n = np.zeros([self.n_nodes, ])
        clusters_n[fiedler_vector_n > 0] = 1
        clusters_n[fiedler_vector_n <= 0] = -1

        clusters_c = clusters_c[test_pos]
        clusters_n = clusters_n[test_pos]

        print('Percentage Error with fiedler vector combinatorial: {:.2f}'.format(error_func(labels_bin, clusters_c)))
        print('Percentage Error with fiedler vector normalized: {:.2f}'.format(error_func(labels_bin, clusters_n)))

    def graph_pnorm_interpolation(self, gradient, P, w, labels_bin, x0=None, p=1., **kwargs):
        r"""
        Solve an interpolation problem via gradient p-norm minimization.

        A signal :math:`x` is estimated from its measurements :math:`y = A(x)` by solving
        :math:`\text{arg}\underset{z \in \mathbb{R}^n}{\min}
        \| \nabla_G z \|_p^p \text{ subject to } Az = y`
        via a primal-dual, forward-backward-forward algorithm.

        Parameters
        ----------
        gradient : array_like
            A matrix representing the graph gradient operator
        P : callable
            Orthogonal projection operator mapping points in :math:`z \in \mathbb{R}^n`
            onto the set satisfying :math:`A P(z) = A z`.
        x0 : array_like, optional
            Initial point of the iteration. Must be of dimension n.
            (Default is `numpy.random.randn(n)`)
        p : {1., 2.}
        labels_bin : array_like
            A vector that holds the binary labels.
        kwargs :
            Additional solver parameters, such as maximum number of iterations
            (maxit), relative tolerance on the objective (rtol), and verbosity
            level (verbosity). See :func:`pyunlocbox.solvers.solve` for the full
            list of options.

        Returns
        -------
        x : array_like
            The solution to the optimization problem.

        """

        grad = lambda z: gradient.dot(z)
        div = lambda z: gradient.transpose().dot(z)

        # Indicator function of the set satisfying :math:`y = A(z)`
        f = functions.func()
        f._eval = lambda z: 0
        f._prox = lambda z, gamma: P(z, w, labels_bin)

        # :math:`\ell_1` norm of the dual variable :math:`d = \nabla_G z`
        g = functions.func()
        g._eval = lambda z: np.sum(np.abs(grad(z)))
        g._prox = lambda d, gamma: functions._soft_threshold(d, gamma)

        # :math:`\ell_2` norm of the gradient (for the smooth case)
        h = functions.norm_l2(A=grad, At=div)

        stepsize = (0.9 / (1. + scipy.sparse.linalg.norm(gradient, ord='fro'))) ** p

        solver = solvers.mlfbf(L=grad, Lt=div, step=stepsize)

        if p == 1.:
            problem = solvers.solve([f, g, functions.dummy()], x0=x0, solver=solver, **kwargs)
            return problem['sol']
        if p == 2.:
            problem = solvers.solve([f, functions.dummy(), h], x0=x0, solver=solver, **kwargs)
            return problem['sol']
        else:
            return x0

    def P(self, a, w, labels_bin):
        mask_pos = np.where(w == 1)
        b = a.copy()
        b[mask_pos] = labels_bin[mask_pos]
        return b

    def transductive_learning(self, w, thresholds, n_trials, labels_bin, test_ind, p, **kwargs):
        for _, threshold in enumerate(thresholds):
            # Simulate n_trials times
            for _ in range(n_trials):
                # Subsample randomly
                subsampled = labels_bin * w
                # Solve p-norm interpolation
                sol = self.graph_pnorm_interpolation(sparse.csr_matrix(self.S).T, self.P, w, labels_bin,
                                                x0=subsampled.copy(), p=p, **kwargs)
                # Threshold to -1 and 1
                sol_bin = (sol > threshold).astype(int)
                sol_bin[sol_bin == 0] = -1
                pdb.set_trace()
                print('Percentage Error with threshold: {:.2f}'.format(error_func(labels_bin[test_ind], sol_bin[test_ind])))

reason: import numpy,import scipy,from scipy,import networkx


text:
"""Run a demo

Usage:
1. Single image: python demo.py image  --model PATH --input IMAGE_PATH
2. Video or webcam: python demo.py video --model PATH --input PATH_OR_CAM_ID
"""
import os
import colorsys
import random
import time
import cv2
import argparse
import numpy as np
import torch
import torchvision
from torchvision.transforms import functional as F
from skimage.measure import find_contours
from PIL import Image

from core.model_factory import get_instance_segmentation_model
from core.tracking.sort import Sort

CLASS_MAP = {
    1: 'tyre',
    2: 'car',
    3: 'man',
    4: 'box',
    5: 'big_box',
    6: 'small_box',
    7: 'suitcase',
    8: 'bag',
    9: 'metal_suitcase',
    10: 'cylinder'
}

def random_colors(N, bright=True):
    """
    Generate random colors.
    To get visually distinct colors, generate them in HSV space then
    convert to RGB.
    """
    brightness = 1.0 if bright else 0.7
    hsv = [(i / N, 1, brightness) for i in range(N)]
    colors = list(map(lambda c: colorsys.hsv_to_rgb(*c), hsv))
    random.shuffle(colors)
    return colors

def apply_mask(image, mask, color, alpha=0.5):
    """Apply the given mask to the image.
    """
    for c in range(3):
        image[:, :, c] = np.where(mask == 1,
                                  image[:, :, c] *
                                  (1 - alpha) + alpha * color[c] * 255,
                                  image[:, :, c])
    return image

def draw_instances(image, boxes, masks, classes, scores, colors):
    """
    boxes: [num_instance, (y1, x1, y2, x2, class_id)] in image coordinates.
    masks: [num_instances, height, width]
    class_ids: [num_instances]
    labels: list of class names of the dataset
    scores: (optional) confidence scores for each box
    colors: TODO
    """
    # Number of instances
    N = boxes.shape[0]
    if not N:
        print("\n*** No instances to display *** \n")
    else:
        assert boxes.shape[0] == masks.shape[0]

    # Show area outside image boundaries.
    height, width = image.shape[:2]

    masked_image = image.copy()
    for i in range(N):
        score = scores[i] if scores is not None else None
        color = colors[i]
        if not np.any(boxes[i]):
            # Skip this instance. Has no bbox. Likely lost in image cropping.
            continue

        x1, y1, x2, y2 = boxes[i]
        cv2.rectangle(masked_image, (x1, y1), (x2, y2),
                      [int(x*255) for x in (color)], 2)

        # Label
        class_id = classes[i]
        label = CLASS_MAP[class_id]
        text = "{}: {}%".format(label, int(score*100))

        yyy = y1 - 16
        if yyy < 0:
            yyy = 0

        cv2.putText(masked_image, text, (x1, yyy), cv2.FONT_HERSHEY_SIMPLEX,
                    1., [int(x*255) for x in (color)], 1, cv2.LINE_AA)

        # Mask
        mask = masks[i, :, :]
        masked_image = apply_mask(masked_image, mask, color)

        # Mask Polygon
        # Pad to ensure proper polygons for masks that touch image edges.
        padded_mask = np.zeros(
            (mask.shape[0] + 2, mask.shape[1] + 2), dtype=np.uint8)
        padded_mask[1:-1, 1:-1] = mask
        contours = find_contours(padded_mask, 0.5)
        for verts in contours:
            # Subtract the padding and flip (y, x) to (x, y)
            verts = np.fliplr(verts) - 1
            pts = np.array(verts.tolist(), np.int32)
            pts = pts.reshape((-1, 1, 2))
            cv2.polylines(masked_image, [pts], True,
                          [int(x*255) for x in (color)], 2)
    return masked_image.astype(np.uint8)

def detect(model, image, device, conf_threshold=0.5, mask_threshold=0.5):
    """Detect on single image"""
    image = F.to_tensor(image)
    model.eval()
    with torch.no_grad():
        prediction = model([image.to(device)])

    scores = prediction[0]['scores'].cpu().numpy()
    conf_mask = scores > conf_threshold

    masks = (prediction[0]['masks'] > mask_threshold).byte().cpu().numpy()
    masks = np.array([m[0] for m in masks])
    boxes = prediction[0]['boxes'].to(torch.int64).cpu().numpy()
    classes = prediction[0]['labels'].byte().cpu().numpy()

    boxes = boxes[conf_mask]
    masks = masks[conf_mask]
    classes = classes[conf_mask]
    scores = scores[conf_mask]
    return boxes, masks, classes, scores

def detect_image(model, image, device,
                 conf_threshold=0.5, mask_threshold=0.5, out_file=None):
    """Detect and display result on a single image."""
    start = time.time()
    boxes, masks, classes, scores = detect(model, image, device,
                                           conf_threshold, mask_threshold)
    print('Detection time: {:.3f}s'.format(time.time() - start))
    image = np.array(image)
    cmap = random_colors(len(CLASS_MAP) + 1)
    colors = [cmap[class_id] for class_id in classes]
    masked_image = draw_instances(image, boxes, masks, classes, scores, colors)

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

    if out_file:
        cv2.imwrite(out_file, masked_image)
    print(f'Saved the masked image as {out_file}')

def detect_video(model, video, device, conf_threshold=0.5,
                 mask_threshold=0.5, out_file=None):
    """Detect and display result on a video or webcam"""
    cmap = random_colors(len(CLASS_MAP) + 1)
    cap = cv2.VideoCapture(video)
    writer = None
    while True:
        ret, image_raw = cap.read()
        if not ret:
            break

        image = Image.fromarray(image_raw).convert('RGB')
        start = time.time()
        boxes, masks, classes, scores = detect(model, image, device,
                                               conf_threshold,
                                               mask_threshold)
        print('FPS: {:3d}'.format(int(1. / (time.time() - start))))
        colors = [cmap[class_id] for class_id in classes]
        masked_image = draw_instances(image_raw, boxes, masks,
                                      classes, scores, colors)
        cv2.imshow('img', masked_image)
        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

        if out_file and writer is None:
            h, w = image_raw.shape[:2]
            writer = cv2.VideoWriter(out_file,
                                     cv2.VideoWriter_fourcc('M', 'J', 'P', 'G'),
                                     24, (w, h))
        if writer:
            writer.write(masked_image)
    cap.release()
    if writer:
        writer.release()
        print(f'Saved masked video as {out_file}')
    cv2.destroyAllWindows()

def detect_video_with_tracking(model, video, device, conf_threshold=0.5,
                               mask_threshold=0.5, out_file=None):
    """Detect, track direction and display result on a video or webcam."""
    cmap = random_colors(len(CLASS_MAP) + 1)
    cap = cv2.VideoCapture(video)
    writer = None

    tracker = Sort()
    while True:
        ret, image_raw = cap.read()
        if not ret:
            break

        image = Image.fromarray(image_raw).convert('RGB')
        start = time.time()
        boxes, masks, classes, scores = detect(model, image,
                                               device, mask_threshold)
        print('FPS: {:3d}'.format(int(1. / (time.time() - start))))

        # Update multi-tracker
        colors = [cmap[class_id] for class_id in classes]

        if boxes.shape[0] > 0:
            track_bb_ids, traces = tracker.update(boxes)
            for i, box_and_id in enumerate(track_bb_ids):
                box = list(map(int, box_and_id[:4]))
                track_id = int(box_and_id[-1])
                trace = traces[i]
                cv2.putText(image_raw, str(track_id), (box[2] + 10, box[1]),
                            cv2.FONT_HERSHEY_SIMPLEX, 1.,
                            colors[i], 1, cv2.LINE_AA)

                # calculate changing direction and draw history
                if len(trace) >= 5:
                    top_trace = [b[1] for b in list(trace)[-5:]]
                    direction_hist = (np.diff(top_trace) > 0).sum()
                    direction = 'forward' if direction_hist < 3 else 'backward'
                    cv2.putText(image_raw, direction, (box[2] + 10, box[3]),
                                cv2.FONT_HERSHEY_SIMPLEX, 1, (0, 0, 255),
                                1, cv2.LINE_AA)
                    for b in trace:
                        b = list(map(int, b))
                        center = ((b[0] + b[2]) // 2, (b[1] + b[3]) // 2)
                        cv2.circle(image_raw, center, 2, (0, 0, 255), -1)

        masked_image = draw_instances(image_raw, boxes, masks,
                                      classes, scores, colors)
        cv2.imshow('img', masked_image)
        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

        if out_file and writer is None:
            h, w = image_raw.shape[:2]
            writer = cv2.VideoWriter(out_file,
                                     cv2.VideoWriter_fourcc('M', 'J', 'P', 'G'),
                                     24, (w, h))
        if writer:
            writer.write(masked_image)
    cap.release()
    if writer:
        writer.release()
        print(f'Saved masked video as {out_file}')
    cv2.destroyAllWindows()

def main():
    parent_parser = argparse.ArgumentParser(add_help=False)
    parent_parser.add_argument('--model', default=None, type=str,
                               help='Path to model')
    parent_parser.add_argument('--vis_threshold', default=0.5, type=float,
                               help='Visualization threshold')
    parent_parser.add_argument('--mask_threshold', default=0.5, type=float,
                               help='Mask threshold')
    parent_parser.add_argument('--conf_threshold', default=0.5, type=float,
                               help='Confidence threshold')
    parent_parser.add_argument('--out_dir',
                               default='saved_models/final_output', type=str,
                               help='Result output dir if specified')

    parser = argparse.ArgumentParser(add_help=False)
    subparser = parser.add_subparsers()
    video_parser = subparser.add_parser('video', parents=[parent_parser])
    video_parser.add_argument('--input', default='0', type=str,
                              help='Path to video input or webcam id')
    video_parser.add_argument('--with_tracking',
                              action='store_true',
                              help='Whether to tracking')
    video_parser.set_defaults(input_type='video')

    image_parser = subparser.add_parser('image', parents=[parent_parser])
    image_parser.add_argument('--input', default='', type=str,
                              help='Path to image input')
    image_parser.set_defaults(input_type='image')
    args = parser.parse_args()

    if not os.path.exists(args.model):
        raise FileNotFoundError(f'Not found model file {args.model}')

    # load model
    num_classes = 11
    device = (torch.device('cuda')
              if torch.cuda.is_available() else torch.device('cpu'))
    model = get_instance_segmentation_model(num_classes)
    model.to(device)

    print(f'Loading model from {args.model}')
    ckpt = torch.load(args.model, map_location=device)
    model.load_state_dict(ckpt['model'])

    input_type = args.input_type
    if input_type == 'video':
        video_out_file = None
        if args.out_dir:
            video_out_dir = os.path.join(args.out_dir, 'video_outputs')
            os.makedirs(video_out_dir, exist_ok=True)
            video_out_file = os.path.join(video_out_dir, 'output.avi')
        video = args.input
        if video.isnumeric():
            video = int(video)
        if isinstance(video, str) and not os.path.exists(video):
            raise FileNotFoundError(f'Not found video file {video}')

        if args.with_tracking:
            detect_video_with_tracking(model, video, device,
                                       args.conf_threshold,
                                       args.mask_threshold,
                                       video_out_file)
        else:
            detect_video(model, video, device,
                         args.conf_threshold,
                         args.mask_threshold, video_out_file)
    elif input_type == 'image':
        if not os.path.exists(args.input):
            raise FileNotFoundError(f'Not found image file {args.image}')
        image = Image.open(args.input)
        image_out_file = None
        if args.out_dir:
            image_out_dir = os.path.join(args.out_dir, 'image_outputs')
            os.makedirs(image_out_dir, exist_ok=True)
            image_out_file = os.path.join(image_out_dir, 'output.jpg')
        detect_image(model, image, device,
                     args.conf_threshold, args.mask_threshold,
                     image_out_file)
    else:
        raise ValueError('Unsupported input type')

if __name__ == '__main__':
    main()

max_stars_repo_name: kkthej/Pytorch-Mask-R-CNN-on-UMCD

max_issues_repo_name: kkthej/Pytorch-Mask-R-CNN-on-UMCD

max_forks_repo_name: kkthej/Pytorch-Mask-R-CNN-on-UMCD


text:
# PyVot Python Variational Optimal Transportation
# Author: Liang Mi <icemiliang@gmail.com>
# Date: April 28th 2020
# Licence: MIT

import numpy as np
from scipy.spatial.distance import cdist
from scipy.optimize import minimize
import imageio
import warnings
import utils

class Vot:
    """ variational optimal transportation """

    def __init__(self, data_p, data_e, label_p=None, label_e=None,
                 weight_p=None, weight_e=None, thres=1e-3, verbose=True):
        """ set up parameters

        p are centroids or source samples
        e are empirical or target samples
        In some literature, definitions of source and target are swapped.

        Throughout PyVot, the term "weight" is referred to the pre-defined value
        for each sample; the term "mass" of a p sample is referred to the weighted summation of
        all the e samples that are indexed to that p

        Args:
            data_p (numpy ndarray): coordinates of p
            data_e (numpy ndarray): coordinates of e
            label_p (numpy ndarray): labels of p
            label_e (numpy ndarray): labels of e
            weight_p (numpy ndarray): weights of p
            weight_e (numpy ndarray): weights of e
            thres (float): threshold to break loops
            verbose (bool): console output verbose flag

        Atts:
            y (numpy ndarray): coordinates of p
            x (numpy ndarray): coordinates of e
            label_y (numpy ndarray): labels of p
            label_e (numpy ndarray): labels of e
            weight_p (numpy ndarray): weight of p
            weight_e (numpy ndarray): weight of e
            mass_p (numpy ndarray): mass of p
            thres    (float): Threshold to break loops
            verbose   (bool): console output verbose flag
        """

        if not isinstance(data_p, np.ndarray):
            raise Exception('y is not a numpy ndarray')
        if not isinstance(data_e, np.ndarray):
            raise Exception('x is not a numpy ndarray')

        if label_p is not None and not isinstance(label_p, np.ndarray):
            raise Exception('label_y is not a numpy ndarray')
        if label_e is not None and not isinstance(label_e, np.ndarray):
            raise Exception('label_e is not a numpy ndarray')

        if weight_p is not None and not isinstance(weight_p, np.ndarray):
            raise Exception('weight_p is not a numpy ndarray')
        if weight_e is not None and not isinstance(weight_e, np.ndarray):
            raise Exception('weight_e is not a numpy ndarray')

        # deep copy all the data?
        self.data_p = data_p
        self.data_e = data_e
        self.data_p_original = self.data_p.copy()

        num_p = data_p.shape[0]
        num_e = data_e.shape[0]

        self.label_p = label_p
        self.label_e = label_e

        self.thres = thres
        self.verbose = verbose

        self.weight_p = weight_p if weight_p is not None else np.ones(num_p) / num_p
        self.weight_e = weight_e if weight_e is not None else np.ones(num_e) / num_e

        utils.assert_boundary(self.data_p)
        utils.assert_boundary(self.data_e)

    def cluster(self, lr=0.5, max_iter_p=10, max_iter_h=5000, lr_decay=500, early_stop=-1):
        """ compute Wasserstein clustering

        Args:
            reg_type   (int): specify regulazation term, 0 means no regularization
            reg        (int): regularization weight
            max_iter_p (int): max num of iteration of clustering
            max_iter_h (int): max num of updating h
            lr       (float): GD learning rate
            lr_decay (float): learning rate decay

        Returns:
            idx (numpy ndarray): assignment of e to p
            pred_label_e (numpy ndarray): labels of e that come from nearest p

        See Also
        --------
        update_p : update p
        update_map: compute optimal transportation
        """
        e_idx, pred_label_e = None, None
        for iter_p in range(max_iter_p):
            dist = cdist(self.data_p, self.data_e) ** 2
            e_idx, pred_label_e = self.update_map(dist, max_iter_h, lr=lr, lr_decay=lr_decay, early_stop=early_stop)
            if self.update_p(e_idx, iter_p):
                break
        return e_idx, pred_label_e

    def update_map(self, dist, max_iter=3000, lr=0.5, beta=0.9, lr_decay=200, early_stop=200):
        """ update assignment of each e as the map to p

        Args:
            dist (numpy ndarray): dist matrix across p and e
            max_iter   (int): max num of iterations
            lr       (float): gradient descent learning rate
            beta     (float): GD momentum
            lr_decay (int): learning rate decay frequency
            early_stop (int): early_stop check frequency

        Returns:
            idx (numpy ndarray): assignment of e to p
            pred_label_e (numpy ndarray): labels of e that come from nearest p
        """

        num_p = self.data_p.shape[0]
        dh = 0
        e_idx = None
        running_median, previous_median = [], 0

        for i in range(max_iter):
            # find nearest p for each e and add mass to p
            e_idx = np.argmin(dist, axis=0)
            mass_p = np.bincount(e_idx, weights=self.weight_e, minlength=num_p)
            # gradient descent with momentum and decay
            dh = beta * dh + (1-beta) * (mass_p - self.weight_p)
            if i != 0 and i % lr_decay == 0:
                lr *= 0.5
            # update dist matrix
            dist += lr * dh[:, None]

            # check if converge
            max_change = np.max((mass_p - self.weight_p)/self.weight_p)
            if max_change.size > 1:
                max_change = max_change[0]
            max_change *= 100

            if self.verbose and ((i < 100 and i % 10 == 0) or i % 100 == 0):
                print("{0:d}: mass diff {1:.2f}%".format(i, max_change))

            if max_change < 1:
                if self.verbose:
                    print("{0:d}: mass diff {1:.2f}%".format(i, max_change))
                break

            # early stop if loss does not decrease TODO better way to early stop?
            if early_stop >= 1:
                running_median.append(max_change)
                if len(running_median) >= early_stop:
                    if previous_median != 0 and\
                            np.abs(np.median(np.asarray(running_median))-previous_median) / previous_median < 0.02:
                        if self.verbose:
                            print("loss saturated, early stopped")
                        break
                    else:
                        previous_median = np.median(np.asarray(running_median))
                        running_median = []

        # labels come from p
        pred_label_e = self.label_p[e_idx] if self.label_p is not None else None

        return e_idx, pred_label_e

    @staticmethod
    def update_p_base(e_idx, data_p, data_e):
        """ base function to update each p to the centroids of its cluster

        Args:
            e_idx (numpy ndarray): assignment of e to p
            data_p (numpy ndarray): cluster centroids, p
            data_e (numpy ndarray): empirical samples, e
            p0 (numpy ndarray): iteration index

        Returns:
            p0 (numpy ndarray): new p
            max_change_pct (float): max_change
        """
        p0 = np.zeros_like(data_p)
        num_p = data_p.shape[0]

        max_change_pct = 0.0
        # update p to the centroid of its clustered e samples
        bincount = np.bincount(e_idx, minlength=num_p)
        if 0 in bincount:
            print('Empty cluster found, optimal transport probably did not converge\n'
                  'Try larger lr or max_iter after checking the measures.')
            # return False
        eps = 1e-8
        for i in range(data_p.shape[1]):
            # update p to the centroid of their correspondences one dimension at a time
            p_target = np.bincount(e_idx, weights=data_e[:, i], minlength=num_p) / (bincount+eps)
            change_pct = np.max(np.abs((data_p[:, i] - p_target) / (data_p[:, i])+eps))
            max_change_pct = max(max_change_pct, change_pct)
            p0[:, i] = p_target

        # replace nan by original data TODO replace nan by nn barycenter?
        mask = np.isnan(p0).any(axis=1)
        p0[mask] = data_p[mask].copy()

        return p0, max_change_pct

    def update_p(self, e_idx, iter_p=0):
        """ update each p to the centroids of its cluster

        Args:
            e_idx (numpy ndarray): assignment of e to p
            iter_p (int): iteration index

        Returns:
            (bool): convergence or not, determined by max p change
        """

        p0, max_change_pct = self.update_p_base(e_idx, self.data_p, self.data_e)

        self.data_p = p0
        if self.verbose:
            print("it {0:d}: max centroid change {1:.2f}".format(iter_p, max_change_pct))

        # return convergence or not
        return True if max_change_pct < self.thres else False

class VotReg(Vot):
    """ variational optimal transportation with regularization on sample supports"""

    def __init__(self, data_p, data_e, label_p=None, label_e=None,
                 weight_p=None, weight_e=None, thres=1e-3, verbose=True):
        super(VotReg, self).__init__(data_p, data_e, label_p=label_p, label_e=label_e,
                                     weight_p=weight_p, weight_e=weight_e, thres=thres, verbose=verbose)

    def cluster(self, reg_type=0, reg=0.01, lr=0.5, max_iter_p=10, max_iter_h=3000, lr_decay=200, early_stop=-1):
        """ compute Wasserstein clustering

        Args:
            reg_type   (int): specify regulazation term, 0 means no regularization
            reg      (float): regularization weight
            lr       (float): GD learning rate
            max_iter_p (int): max num of iteration of clustering
            max_iter_h (int): max num of updating h
            lr_decay   (int): learning rate decay interval

        See Also
        --------
        update_p : update p
        update_map: compute optimal transportation
        """
        e_idx, pred_label_e = None, None
        for iter_p in range(max_iter_p):
            dist = cdist(self.data_p, self.data_e) ** 2
            e_idx, pred_label_e = self.update_map(dist, max_iter_h, lr=lr, lr_decay=lr_decay, early_stop=early_stop)
            if self.update_p(e_idx, iter_p, reg_type, reg):
                break
        return e_idx, pred_label_e

    def update_p(self, e_idx, iter_p=0, reg_type=0, reg=0.01):
        """ update p

        Args:
            e_idx (numpy ndarray): assignment of e to p
            iter_p (int): iteration index
            reg_type (int or string): regularization type
            reg (float): regularizer weight

        Returns:
            bool: convergence or not
        """

        if reg_type == 1 or reg_type == 'potential':
            return self.update_p_reg_potential(e_idx, iter_p, reg)
        elif reg_type == 2 or reg_type == 'transform':
            return self.update_p_reg_transform(e_idx, iter_p, reg)
        else:
            return self.update_p_noreg(e_idx, iter_p)

    def update_p_noreg(self, e_idx, iter_p=0):
        """ update each p to the centroids of its cluster

        Args:
            e_idx (numpy ndarray): assignment of e to p
            iter_p (int): iteration index

        Returns:
            bool: convergence or not, determined by max p change
        """

        p0, max_change_pct = self.update_p_base(e_idx, self.data_p, self.data_e)
        self.data_p = p0
        # replace nan by original data
        mask = np.isnan(self.data_p).any(axis=1)
        self.data_p[mask] = self.data_p_original[mask].copy()
        if self.verbose:
            print("it {0:d}: max centroid change {1:.2f}%".format(iter_p, 100 * max_change_pct))
        # return max p coor change
        return True if max_change_pct < self.thres else False

    def update_p_reg_potential(self, e_idx, iter_p=0, reg=0.01):
        """ update each p to the centroids of its cluster,
            regularized by intra-class distances

        Args:
            e_idx (numpy ndarray): assignment of e to p
            iter_p (int): index of the iteration of updating p
            reg (float): regularizer weight

        Returns:
            bool: convergence or not, determined by max p change
        """

        def f(p, p0, label=None, reg=0.01):
            """ objective function incorporating labels

            Args:
                p  (numpy ndarray): p
                p0 (numpy ndarray): centroids of e
                label (numpy ndarray): labels of p
                reg (float): regularizer weight

            Returns:
                float: f = sum(|p-p0|^2) + reg * sum(1(li == lj)*|pi-pj|^2)
            """

            p = p.reshape(p0.shape)
            reg_term = 0.0
            for l in np.unique(label):
                p_sub = p[label == l, :]
                # pairwise distance with smaller memory burden
                # |pi - pj|^2 = pi^2 + pj^2 - 2*pi*pj
                reg_term += np.sum((p_sub ** 2).sum(axis=1, keepdims=True) +
                                   (p_sub ** 2).sum(axis=1) -
                                   2 * p_sub.dot(p_sub.T))

            return np.sum((p - p0) ** 2.0) + reg * reg_term

        if np.unique(self.label_p).size == 1:
            warnings.warn("All known samples belong to the same class")

        p0, max_change_pct = self.update_p_base(e_idx, self.data_p, self.data_e)

        if self.verbose:
            print("it {0:d}: max centroid change {1:.2f}".format(iter_p, max_change_pct))

        # regularize
        res = minimize(f, self.data_p, method='BFGS', args=(p0, self.label_p, reg))
        self.data_p = res.x.reshape(p0.shape)
        # return convergence or not
        return True if max_change_pct < self.thres else False

    def update_p_reg_transform(self, e_idx, iter_p=0, reg=0.01):
        """ update each p to the centroids of its cluster,
            regularized by an affine transformation
            which is estimated from the OT map.

        Args:
            e_idx (numpy ndarray): assignment of e to p
            iter_p (int): index of the iteration of updating p
            reg (float): regularizer weight

        Returns:
            bool: convergence or not, determined by max p change
        """

        assert self.data_p.shape[1] == 2, "dim has to be 2 for geometric transformation"

        p0, max_change_pct = self.update_p_base(e_idx, self.data_p, self.data_e)

        if self.verbose:
            print("it {0:d}: max centroid change {1:.2f}".format(iter_p, max_change_pct))

        pt = self.data_p.copy()
        pt = utils.estimate_transform_target(pt, p0)

        # regularize within each label
        # pt = np.zeros(p0.shape)
        # for label in np.unique(self.label_y):
        #     idx_p_label = self.label_y == label
        #     p_sub = self.y[idx_p_label, :]
        #     p0_sub = p0[idx_p_label, :]
        #     T = tf.EuclideanTransform()
        #     # T = tf.AffineTransform()
        #     # T = tf.ProjectiveTransform()
        #     T.estimate(p_sub, p0_sub)
        #     pt[idx_p_label, :] = T(p_sub)
        #
        # pt = self.y.copy()
        # T = tf.EuclideanTransform()
        # T.estimate(pt, p0)
        # pt = T(pt)

        self.data_p = 1 / (1 + reg) * p0 + reg / (1 + reg) * pt
        # return convergence
        return True if max_change_pct < self.thres else False

class VOTAP:
    """
        y are the centroids
        x are the area samples
        This is a minimum class for area-preserving maps
    """

    def __init__(self, data, sampling='square', label=None, nu=None, thres=1e-5, ratio=100, verbose=False):
        """ set up parameters
        """

        if not isinstance(data, np.ndarray):
            raise Exception('input is not a numpy ndarray')

        if label is not None and not isinstance(label, np.ndarray):
            raise Exception('label is neither a numpy array not a numpy ndarray')

        if nu is not None and not isinstance(nu, np.ndarray):
            raise Exception('label is neither a numpy array not a numpy ndarray')

        self.y = data
        self.data_p_original = self.y.copy()
        self.K = self.y.shape[0]

        self.label_y = label
        self.weight_p = nu if nu is not None else np.ones(self.K) / self.K

        self.thres = thres
        self.verbose = verbose

        utils.assert_boundary(self.y)

        self.N0 = int(ratio * self.K)
        ndim = self.y.shape[1]
        self.x, _ = utils.random_sample(self.N0, ndim, sampling=sampling)

        self.dist = cdist(self.y, self.x, 'sqeuclidean')

    def map(self, plot_filename=None, beta=0.9, max_iter=1000, lr=0.5, lr_decay=200, early_stop=100):
        """ map y into the area
        """

        imgs = []
        dh = 0

        idx = None
        running_median, previous_median = [], 0

        for i in range(max_iter):
            # find nearest p for each e
            idx = np.argmin(self.dist, axis=0)

            # calculate total mass of each cell
            mass_p = np.bincount(idx, minlength=self.K) / self.N0
            # gradient descent with momentum and decay
            dh = beta * dh + (1-beta) * (mass_p - self.weight_p)
            if i != 0 and i % lr_decay == 0:
                lr *= 0.9
            self.dist += lr * dh[:, None]

            # plot to gif, TODO this is time consuming, got a better way?
            if plot_filename and i % 10 == 0:
                fig = utils.plot_map(self.x, idx / (self.K - 1))
                img = utils.fig2data(fig)
                imgs.append(img)

            # check if converge
            max_change = np.max((mass_p - self.weight_p) / self.weight_p)
            if max_change.size > 1:
                max_change = max_change[0]
            max_change *= 100

            if self.verbose and ((i < 100 and i % 10 == 0) or i % 100 == 0):
                print("{0:d}: mass diff {1:.2f}%".format(i, max_change))

            if max_change < 1:
                if self.verbose:
                    print("{0:d}: mass diff {1:.2f}%".format(i, max_change))
                break

            if early_stop > 0:
                # early stop if loss does not decrease TODO better way to early stop?
                running_median.append(max_change)
                if len(running_median) >= early_stop:
                    if previous_median != 0 and \
                            np.abs(np.median(np.asarray(running_median)) - previous_median) / previous_median < 0.02:
                        if self.verbose:
                            print("loss saturated, early stopped")
                        break
                    else:
                        previous_median = np.median(np.asarray(running_median))
                        running_median = []

            if max_change <= 1:
                break
        if plot_filename and imgs:
            imageio.mimsave(plot_filename, imgs, fps=4)
        # labels come from y
        pred_label_x = self.label_y[idx] if self.label_y is not None else None

        # update coordinates of y
        bincount = np.bincount(idx, minlength=self.K)
        if 0 in bincount:
            print('Empty cluster found, optimal transport probably did not converge\nTry larger lr or max_iter')
            # return
        for i in range(self.y.shape[1]):
            # update p to the centroid of their correspondences
            self.y[:, i] = np.bincount(idx, weights=self.x[:, i], minlength=self.K) / bincount

        return idx, pred_label_x

class VOT:
    def __init__(self, y, x, nu=None, mu=None, lam=None, label_y=None, label_x=None, tol=1e-4, verbose=True):

        # marginals (x, mu)
        # centroids (y, nu)

        if type(x) is np.ndarray:
            if x.ndim == 2:
                self.x = [x]
            elif x.ndim == 3:
                self.x = [x[i] for i in range(x.shape[0])]
        else:
            self.x = x

        self.y = y.copy()
        self.y_original = y

        self.K = y.shape[0]  # number of centroids
        self.n = y.shape[1]  # number of dimensions
        self.N = len(self.x)  # number of marginals

        self.tol = tol
        self.verbose = verbose

        self.lam = lam if lam is not None else np.ones(self.N) / self.N

        self.idx = []
        self.mu = []
        self.sum_mu = []
        if mu is not None:
            # copy mu
            if type(mu) is np.ndarray:
                self.mu = [mu]
            else:
                self.mu = mu
            for m in self.mu:
                self.idx.append(np.ones_like(m, dtype=np.int64))
                self.sum_mu.append(np.sum(m))
        else:
            # create uniform mu
            self.mu = []
            self.idx = []
            for i in range(self.N):
                N_i = self.x[i].shape[0]
                self.mu.append(1. / N_i)
                self.idx.append(np.zeros(N_i, dtype=np.int64))
                self.sum_mu.append(1.)

        if nu is not None:
            self.nu = nu
            self.sum_nu = np.sum(self.nu)
            if abs(self.sum_nu - 1) > 1e-3:
                self.nu /= self.sum_nu
                self.sum_nu = 1
                self.mu = [m / self.sum_nu for m in self.mu]
        else:
            self.nu = 1. / self.K
            self.sum_nu = 1.

        self.label_y = label_y
        self.label_x = []

        # all data should be in (-1, 1) in each dimension
        utils.assert_boundary(self.y)
        for i in range(self.N):
            utils.assert_boundary(self.x[i])

    def cluster(self, lr=0.5, max_iter_y=10, max_iter_h=3000, lr_decay=200, stop=-1, beta=0, reg=0., keep_idx=False, space='euclidean', icp=False):
        """ compute Wasserstein clustering
        """

        lrs = [lr / m for m in self.sum_mu]
        idxs = []
        for it in range(max_iter_y):
            for i in range(self.N):
                print("solving marginal #" + str(i))
                if space == 'spherical':
                    dist = np.matmul(self.y, self.x[i].T)
                else:
                    dist = cdist(self.y, self.x[i], 'sqeuclidean')
                output = self.update_map(i, dist, max_iter_h, lr=lrs[i], lr_decay=lr_decay, beta=beta, stop=stop, reg=reg, keep_idx=keep_idx, space=space)
                self.idx[i] = output['idx']
                if keep_idx:
                    idxs.append(output['idxs'])

            if icp:
                if self.update_x(it):
                    break
            elif self.update_y(it, space=space):
                break
        output = dict()
        output['idxs'] = idxs

        # pass label from y to x
        if self.label_y is not None:
            for i in range(self.N):
                self.label_x.append(self.label_y[self.idx[i]])

        # compute W_2^2
        twd = 0
        wds = []
        for i in range(self.N):
            tmp = (self.y[self.idx[i], :] - self.x[i]) ** 2
            wd = np.sum(np.sum(tmp, axis=1) * self.mu[i])
            twd += wd
            wds.append(wd)

        output['wd'] = twd
        output['wds'] = wds
        return output

    def update_map(self, i, dist, max_iter=3000, lr=0.5, beta=0, lr_decay=200, stop=200, reg=0., keep_idx=False, space='euclidean'):
        """ update assignment of each e as the ot_map to y
        """

        dh = 0
        idx = None
        idxs = []
        running_median, previous_median = [], 0

        h = np.ones(self.K) if space == 'spherical' else None

        dist_original = 0 if reg == 0 else dist.copy()

        for it in range(max_iter):
            # find nearest y for each x and add mass to y
            if space == 'spherical':
                idx = np.argmin(dist / np.cos(h)[:, None], axis=0)
            else:
                idx = np.argmin(dist, axis=0)
            if keep_idx:
                idxs.append(idx)
            if isinstance(self.mu[i], float):
                mass = np.bincount(idx, minlength=self.K) * self.mu[i]
            else:
                mass = np.bincount(idx, weights=self.mu[i], minlength=self.K)

            # gradient descent with momentum and decay
            dh = beta * dh + (1 - beta) * (mass - self.nu)
            if it != 0 and it % lr_decay == 0:
                lr *= 0.5
            # update dist matrix
            dh *= lr
            if space == 'spherical':
                h += dh
            else:
                dist += dh[:, None]

            # check if converge
            if self.verbose and it % 1000 == 0:
                print(dh)
            max_change = np.max((mass - self.nu) / self.nu)
            if max_change.size > 1:
                max_change = max_change[0]
            max_change *= 100

            if self.verbose and ((i < 20 and i % 1 == 0) or i % 200 == 0):
                print("{0:d}: mass diff {1:.2f}%".format(it, max_change))

            if max_change < 1:
                if self.verbose:
                    print("{0:d}: mass diff {1:.2f}%".format(it, max_change))
                break

            # early stop if loss does not decrease TODO better way to early stop?
            if stop >= 1:
                running_median.append(max_change)
                if len(running_median) >= stop:
                    if previous_median != 0 and\
                            np.abs(np.median(np.array(running_median))-previous_median) / previous_median < 0.02:
                        if self.verbose:
                            print("loss saturated, early stopped")
                        break
                    else:
                        previous_median = np.median(np.array(running_median))
                        running_median = []

        if reg != 0.:
            idx = np.argmin(reg / (1 + reg) * dist + 1 / (1 + reg) * dist_original, axis=0)

        output = dict()
        output['idx'] = idx
        output['idxs'] = idxs
        return output

    @staticmethod
    def update_y_base(idx, y, x):
        """ base function to update each y to the centroids of its cluster
        """

        new_y = np.zeros_like(y)
        max_change_pct = 0.0
        K, ndim = y.shape

        bincount = np.bincount(idx, minlength=K)
        if 0 in bincount:
            print('Empty cluster found, OT probably did not converge\n'
                  'Try a different lr or max_iter assuming the input is correct.')
            # return False
        eps = 1e-8

        # update y to the centroid of their correspondences one dimension at a time
        # for spherical domains, use Euclidean barycenter to approximate and project it to the surface
        for n in range(ndim):
            mass_center = np.bincount(idx, weights=x[:, n], minlength=K) / (bincount + eps)
            change_pct = np.max(np.abs((y[:, n] - mass_center) / (y[:, n]) + eps))
            max_change_pct = max(max_change_pct, change_pct)
            new_y[:, n] = mass_center

        # replace nan by original data TODO replace nan by nn barycenter?
        mask = np.isnan(new_y).any(axis=1)
        new_y[mask] = y[mask].copy()

        return new_y, max_change_pct

    def update_y(self, it=0, idx=None, space='euclidean', icp=False):
        """ update each y to the centroids of its cluster
        """
        if idx is None:
            idx = self.idx
        max_change_pct = 1e9

        y = np.zeros((self.N, self.K, self.n))
        if icp:
            yR = np.zeros((self.N, self.K, self.n))
        for i in range(self.N):
            y[i], change = self.update_y_base(idx[i], self.y, self.x[i])
            max_change_pct = max(max_change_pct, change)
            if icp:
                yR[i] = utils.estimate_transform_target(self.y, y[i])

        if icp:
            y = yR

        self.y = np.sum(y * self.lam[:, None, None], axis=0)

        if space == 'spherical':
            self.y /= np.linalg.norm(self.y, axis=1, keepdims=True)

        if self.verbose:
            print("iter {0:d}: max centroid change {1:.2f}%".format(it, 100 * max_change_pct))

        return True if max_change_pct < self.tol else False

    def update_x(self, it=0, idx=None):
        """ update each x
        """
        if idx is None:
            idx = self.idx
        max_change_pct = 1e9

        y = np.zeros((self.N, self.K, self.n))
        for i in range(self.N):
            y[i], change = self.update_y_base(idx[i], self.y, self.x[i])
            max_change_pct = max(max_change_pct, change)

        self.y = np.sum(y * self.lam[:, None, None], axis=0)

        for i in range(self.N):
            r, t = utils.estimate_inverse_transform(self.y, y[i])
            self.x[i] = (np.matmul(r, self.x[i].T) + t).T

        if self.verbose:
            print("iter {0:d}: max centroid change {1:.2f}%".format(it, 100 * max_change_pct))

        return True if max_change_pct < self.tol else False

class VOTREG(VOT):
    """ variational optimal transportation with regularization on sample supports"""

    def map(self, reg_type, reg, lr=0.5, max_iter_y=10, max_iter_h=3000, lr_decay=200, stop=-1, keep_idx=False):
        """ compute Wasserstein clustering
        """
        lrs = [lr / m for m in self.sum_mu]
        idxs = []
        for iter_y in range(max_iter_y):
            dist = cdist(self.y, self.x[0], 'sqeuclidean')
            output = self.update_map(0, dist, max_iter_h, lr=lrs[0], lr_decay=lr_decay, stop=stop)
            self.idx[0] = output['idx']
            if keep_idx:
                idxs.append(output['idxs'])
            if reg_type == 1 or reg_type == 'potential':
                if self.update_y_potential(iter_y, reg):
                    break
            elif reg_type == 2 or reg_type == 'transform':
                if self.update_p_transform(iter_y, reg):
                    break
            else:
                raise Exception('regularization type not defined')

        # pass label from y to x
        if self.label_y is not None:
            for i in range(self.N):
                self.label_x.append(self.label_y[self.idx[i]])

        output = dict()
        output['idxs'] = idxs
        return output

    def update_y_potential(self, iter_y=0, reg=0.01):
        """ update each p to the centroids of its cluster,
        """

        def f(p, p0, label=None, reg=0.01):
            """ objective function incorporating labels

            Args:
                p  (numpy ndarray): p
                p0 (numpy ndarray): centroids of e
                label (numpy ndarray): labels of p
                reg (float): regularizer weight

            Returns:
                float: f = sum(|p-p0|^2) + reg * sum(1(li == lj)*|pi-pj|^2)
            """

            p = p.reshape(p0.shape)
            reg_term = 0.0
            for l in np.unique(label):
                p_sub = p[label == l, :]
                # pairwise distance with smaller memory burden
                # |pi - pj|^2 = pi^2 + pj^2 - 2*pi*pj
                reg_term += np.sum((p_sub ** 2).sum(axis=1, keepdims=True) +
                                   (p_sub ** 2).sum(axis=1) -
                                   2 * p_sub.dot(p_sub.T))

            return np.sum((p - p0) ** 2.0) + reg * reg_term

        if np.unique(self.label_y).size == 1:
            warnings.warn("All known samples belong to the same class")

        y0, max_change_pct = self.update_y_base(self.idx[0], self.y, self.x[0])

        if self.verbose:
            print("it {0:d}: max centroid change {1:.2f}".format(iter_y, max_change_pct))

        # regularize
        res = minimize(f, self.y, method='BFGS', args=(y0, self.label_y, reg))
        self.y = res.x.reshape(y0.shape)

        return True if max_change_pct < self.tol else False

    def update_p_transform(self, iter_p=0, reg=0.01):
        """ update each p to the centroids of its cluster,
        """

        assert self.y.shape[1] == 2, "dim has to be 2 for geometric transformation"

        p0, max_change_pct = self.update_y_base(self.idx[0], self.y, self.x[0])

        if self.verbose:
            print("it {0:d}: max centroid change {1:.2f}".format(iter_p, max_change_pct))

        pt = np.zeros_like(self.y)
        pt = utils.estimate_transform_target(pt, p0)

        # regularize within each label
        # pt = np.zeros(p0.shape)
        # for label in np.unique(self.label_y):
        #     idx_p_label = self.label_y == label
        #     p_sub = self.y[idx_p_label, :]
        #     p0_sub = p0[idx_p_label, :]
        #     T = tf.EuclideanTransform()
        #     # T = tf.AffineTransform()
        #     # T = tf.ProjectiveTransform()
        #     T.estimate(p_sub, p0_sub)
        #     pt[idx_p_label, :] = T(p_sub)
        #
        # pt = self.y.copy()
        # T = tf.EuclideanTransform()
        # T.estimate(pt, p0)
        # pt = T(pt)

        self.y = pt
        # self.y = 1 / (1 + reg) * p0 + reg / (1 + reg) * pt
        # return convergence
        return True if max_change_pct < self.tol else False

max_stars_repo_name: icemiliang/pyvot

max_issues_repo_name: icemiliang/pyvot

max_forks_repo_name: icemiliang/pyvot

reason: import numpy,from scipy


text:
import numpy as np
import tensorflow as tf

class MAML:
	def __init__(self, d, c, nway, meta_lr=1e-3, train_lr=1e-2):
		"""

		:param d:
		:param c:
		:param nway:
		:param meta_lr:
		:param train_lr:
		"""
		self.d = d
		self.c = c
		self.nway = nway
		self.meta_lr = meta_lr
		self.train_lr = train_lr

		print('img shape:', self.d, self.d, self.c, 'meta-lr:', meta_lr, 'train-lr:', train_lr)

	def build(self, support_xb, support_yb, query_xb, query_yb, K, meta_batchsz, mode='train'):
		"""

		:param support_xb:   [b, setsz, 84*84*3]
		:param support_yb:   [b, setsz, n-way]
		:param query_xb:     [b, querysz, 84*84*3]
		:param query_yb:     [b, querysz, n-way]
		:param K:           train update steps
		:param meta_batchsz:tasks number
		:param mode:        train/eval/test, for training, we build train&eval network meanwhile.
		:return:
		"""
		# create or reuse network variable, not including batch_norm variable, therefore we need extra reuse mechnism
		# to reuse batch_norm variables.
		self.weights = self.conv_weights()
		# TODO: meta-test is sort of test stage.
		training = True if mode is 'train' else False

		def meta_task(input):
			"""
			map_fn only support one parameters, so we need to unpack from tuple.
			:param support_x:   [setsz, 84*84*3], e.g. 1-shot, setsz=1
			:param support_y:   [setsz, n-way], e.g. 1-shot, setsz=1
			:param query_x:     [querysz, 84*84*3], e.g. querysz=15 (default case)
			:param query_y:     [querysz, n-way]
			:param training:    training or not, for batch_norm
			:return:
			"""
			support_x, support_y, query_x, query_y = input
			# to record the op in t update step.
			query_preds, query_losses, query_accs = [], [], []

			# ==================================
			# REUSE       True        False
			# Not exist   Error       Create one
			# Existed     reuse       Error
			# ==================================
			# That's, to create variable, you must turn off reuse
			support_pred = self.forward(support_x, self.weights, training)
			support_loss = tf.nn.softmax_cross_entropy_with_logits(logits=support_pred, labels=support_y)
			support_acc = tf.contrib.metrics.accuracy(tf.argmax(tf.nn.softmax(support_pred, dim=1), axis=1),
			                                             tf.argmax(support_y, axis=1))
			# compute gradients
			grads = tf.gradients(support_loss, list(self.weights.values()))
			# grad and variable dict
			gvs = dict(zip(self.weights.keys(), grads))

			# theta_pi = theta - alpha * grads
			fast_weights = dict(zip(self.weights.keys(),
			                        [self.weights[key] - self.train_lr * gvs[key] for key in self.weights.keys()]))
			# use theta_pi to forward meta-test
			query_pred = self.forward(query_x, fast_weights, training)
			# meta-test loss
			query_loss = tf.nn.softmax_cross_entropy_with_logits(logits=query_pred, labels=query_y)
			# record T0 pred and loss for meta-test
			query_preds.append(query_pred)
			query_losses.append(query_loss)

			# continue to build T1-TK steps graph
			for _ in range(1, K):
				# T_k loss on meta-train
				# we need meta-train loss to fine-tune the task and meta-test loss to update theta
				loss = tf.nn.softmax_cross_entropy_with_logits(logits=self.forward(support_x, fast_weights, training),
				                                               labels=support_y)
				# compute gradients
				grads = tf.gradients(loss, list(fast_weights.values()))
				# compose grad and variable dict
				gvs = dict(zip(fast_weights.keys(), grads))
				# update theta_pi according to varibles
				fast_weights = dict(zip(fast_weights.keys(), [fast_weights[key] - self.train_lr * gvs[key]
				                         for key in fast_weights.keys()]))
				# forward on theta_pi
				query_pred = self.forward(query_x, fast_weights, training)
				# we need accumulate all meta-test losses to update theta
				query_loss = tf.nn.softmax_cross_entropy_with_logits(logits=query_pred, labels=query_y)
				query_preds.append(query_pred)
				query_losses.append(query_loss)

			# compute every steps' accuracy on query set
			for i in range(K):
				query_accs.append(tf.contrib.metrics.accuracy(tf.argmax(tf.nn.softmax(query_preds[i], dim=1), axis=1),
					                                            tf.argmax(query_y, axis=1)))
			# we just use the first step support op: support_pred & support_loss, but igonre these support op
			# at step 1:K-1.
			# however, we return all pred&loss&acc op at each time steps.
			result = [support_pred, support_loss, support_acc, query_preds, query_losses, query_accs]

			return result

		# return: [support_pred, support_loss, support_acc, query_preds, query_losses, query_accs]
		out_dtype = [tf.float32, tf.float32, tf.float32, [tf.float32] * K, [tf.float32] * K, [tf.float32] * K]
		result = tf.map_fn(meta_task, elems=(support_xb, support_yb, query_xb, query_yb),
		                   dtype=out_dtype, parallel_iterations=meta_batchsz, name='map_fn')
		support_pred_tasks, support_loss_tasks, support_acc_tasks, \
			query_preds_tasks, query_losses_tasks, query_accs_tasks = result

		if mode is 'train':
			# average loss
			self.support_loss = support_loss = tf.reduce_sum(support_loss_tasks) / meta_batchsz
			# [avgloss_t1, avgloss_t2, ..., avgloss_K]
			self.query_losses = query_losses = [tf.reduce_sum(query_losses_tasks[j]) / meta_batchsz
			                                        for j in range(K)]
			# average accuracy
			self.support_acc = support_acc = tf.reduce_sum(support_acc_tasks) / meta_batchsz
			# average accuracies
			self.query_accs = query_accs = [tf.reduce_sum(query_accs_tasks[j]) / meta_batchsz
			                                        for j in range(K)]

			# # add batch_norm ops before meta_op
			# update_ops = tf.get_collection(tf.GraphKeys.UPDATE_OPS)
			# with tf.control_dependencies(update_ops):
			# 	# TODO: the update_ops must be put before tf.train.AdamOptimizer,
			# 	# otherwise it throws Not in same Frame Error.
			# 	meta_loss = tf.identity(self.query_losses[-1])

			# meta-train optim
			optimizer = tf.train.AdamOptimizer(self.meta_lr, name='meta_optim')
			# meta-train gradients, query_losses[-1] is the accumulated loss across over tasks.
			gvs = optimizer.compute_gradients(self.query_losses[-1])
			# meta-train grads clipping
			gvs = [(tf.clip_by_norm(grad, 10), var) for grad, var in gvs]
			# update theta
			self.meta_op = optimizer.apply_gradients(gvs)

		else: # test & eval

			# average loss
			self.test_support_loss = support_loss = tf.reduce_sum(support_loss_tasks) / meta_batchsz
			# [avgloss_t1, avgloss_t2, ..., avgloss_K]
			self.test_query_losses = query_losses = [tf.reduce_sum(query_losses_tasks[j]) / meta_batchsz
			                                        for j in range(K)]
			# average accuracy
			self.test_support_acc = support_acc = tf.reduce_sum(support_acc_tasks) / meta_batchsz
			# average accuracies
			self.test_query_accs = query_accs = [tf.reduce_sum(query_accs_tasks[j]) / meta_batchsz
			                                        for j in range(K)]

		# NOTICE: every time build model, support_loss will be added to the summary, but it's different.
		tf.summary.scalar(mode + '：support loss', support_loss)
		tf.summary.scalar(mode + '：support acc', support_acc)
		for j in range(K):
			tf.summary.scalar(mode + '：query loss, step ' + str(j + 1), query_losses[j])
			tf.summary.scalar(mode + '：query acc, step ' + str(j + 1), query_accs[j])

	def conv_weights(self):
		weights = {}

		conv_initializer = tf.contrib.layers.xavier_initializer_conv2d()
		fc_initializer = tf.contrib.layers.xavier_initializer()
		k = 3

		with tf.variable_scope('MAML', reuse= tf.AUTO_REUSE):

			weights['conv0']      = tf.get_variable('conv0w', [7, 7, 3, 64],  initializer=conv_initializer)
			weights['b0']         = tf.get_variable('conv0b', initializer=tf.zeros([64]))

			weights['conv1']      = tf.get_variable('conv1w', [k, k, 64, 64],  initializer=conv_initializer)
			weights['b1']         = tf.get_variable('conv1b', initializer=tf.zeros([64]))
			weights['conv1_2']    = tf.get_variable('conv1w_2', [k, k, 64, 64],  initializer=conv_initializer)
			weights['b1_2']       = tf.get_variable('conv1b_2', initializer=tf.zeros([64]))
			weights['conv1_3']    = tf.get_variable('conv1w_3', [1, 1, 64, 64],  initializer=conv_initializer)
			weights['b1_3']       = tf.get_variable('conv1b_3', initializer=tf.zeros([64]))

			weights['conv2']      = tf.get_variable('conv2w', [k, k, 64, 128], initializer=conv_initializer)
			weights['b2']         = tf.get_variable('conv2b', initializer=tf.zeros([128]))
			weights['conv2_2']    = tf.get_variable('conv2w_2', [k, k, 128, 128], initializer=conv_initializer)
			weights['b2_2']       = tf.get_variable('conv2b_2', initializer=tf.zeros([128]))
			weights['conv2_3']    = tf.get_variable('conv2w_3', [1, 1, 64, 128], initializer=conv_initializer)
			weights['b2_3']       = tf.get_variable('conv2b_3', initializer=tf.zeros([128]))

			weights['conv3']      = tf.get_variable('conv3w', [k, k, 128, 256], initializer=conv_initializer)
			weights['b3']         = tf.get_variable('conv3b', initializer=tf.zeros([256]))
			weights['conv3_2']    = tf.get_variable('conv3w_2', [k, k, 256, 256], initializer=conv_initializer)
			weights['b3_2']       = tf.get_variable('conv3b_2', initializer=tf.zeros([256]))
			weights['conv3_3']    = tf.get_variable('conv3w_3', [1, 1, 128, 256], initializer=conv_initializer)
			weights['b3_3']       = tf.get_variable('conv3b_3', initializer=tf.zeros([256]))

			weights['conv4']      = tf.get_variable('conv4w', [k, k, 256, 512], initializer=conv_initializer)
			weights['b4']         = tf.get_variable('conv4b', initializer=tf.zeros([512]))
			weights['conv4_2']    = tf.get_variable('conv4w_2', [k, k, 512, 512], initializer=conv_initializer)
			weights['b4_2']       = tf.get_variable('conv4b_2', initializer=tf.zeros([512]))
			weights['conv4_3']    = tf.get_variable('conv4w_3', [1, 1, 256, 512], initializer=conv_initializer)
			weights['b4_3']       = tf.get_variable('conv4b_3', initializer=tf.zeros([512]))
			# assumes max pooling
			final_imgsize = 1#int(self.d/16)
			# print(final_imgsize)
			weights['w5']       = tf.get_variable('fc1w', [512 * final_imgsize * final_imgsize, self.nway], initializer=fc_initializer)
			weights['b5']       = tf.get_variable('fc1b', initializer=tf.zeros([self.nway]))

			return weights

	def conv_block(self, input_tensor, weight, bias, strides, scope, training):
		"""
		build a block with conv2d->batch_norm->pooling
		:param x:
		:param weight:
		:param bias:
		:param scope:
		:param training:
		:return:
		"""
		# conv
		x = tf.nn.conv2d(input_tensor, weight[0], [1, strides, strides, 1], 'SAME', name=scope + '_conv2d') + bias[0]
		# batch norm, activation_fn=tf.nn.relu,
		# NOTICE: must have tf.layers.batch_normalization
		# x = tf.contrib.layers.batch_norm(x, activation_fn=tf.nn.relu)
		with tf.variable_scope('MAML'):
			# train is set to True ALWAYS, please refer to https://github.com/cbfinn/maml/issues/9
			# when FLAGS.train=True, we still need to build evaluation network
			x = tf.layers.batch_normalization(x, training=True, name=scope + '_bn', reuse=tf.AUTO_REUSE)
		# relu
		x = tf.nn.relu(x, name=scope + '_relu')
		x = tf.nn.conv2d(x, weight[1], [1, 1, 1, 1], 'SAME', name=scope + '_conv2d_2') + bias[1]
		with tf.variable_scope('MAML'):
			x = tf.layers.batch_normalization(x, training=True, name=scope + '_bn_2', reuse=tf.AUTO_REUSE)
		x = tf.nn.relu(x, name=scope + '_relu_2')
		# # pooling
		# x = tf.nn.max_pool(x, [1, 2, 2, 1], [1, 2, 2, 1], 'VALID', name=scope + '_pool')
		shortcut = tf.nn.conv2d(input_tensor, weight[2], [1, strides, strides, 1], 'VALID', name=scope + '_conv2d_3') + bias[2]
		with tf.variable_scope('MAML'):
			shortcut = tf.layers.batch_normalization(shortcut, training=True, name=scope + '_bn_3', reuse=tf.AUTO_REUSE)

		x = tf.add(shortcut, x, name=scope + '_add_shortcut')
		x = tf.nn.relu(x, name=scope + '_relu_3')
		return x

	def forward(self, x, weights, training):
		"""

		:param x:
		:param weights:
		:param training:
		:return:
		"""
		# [b, 84, 84, 3]
		x = tf.reshape(x, [-1, self.d, self.d, self.c], name='reshape1')
		padding = np.array([[0,0], [3,3], [3,3],[0,0]])
		x = tf.pad(x, padding)
		hidden0 = tf.nn.conv2d(x, weights['conv0'], [1, 2, 2, 1], 'VALID', name='conv00' + '_conv2d') + weights['b0']
		with tf.variable_scope('MAML'):
			# train is set to True ALWAYS, please refer to https://github.com/cbfinn/maml/issues/9
			# when FLAGS.train=True, we still need to build evaluation network
			hidden0 = tf.layers.batch_normalization(hidden0, training=True, name='conv00' + '_bn', reuse=tf.AUTO_REUSE)
		hidden0 = tf.nn.relu(hidden0, name='conv00' + '_relu')
		hidden0 = tf.nn.max_pool(hidden0, [1, 2, 2, 1], [1, 2, 2, 1], 'VALID', name='conv00' + '_pool')
		# print("="*50)
		# print(hidden0.shape)
		# print("="*50)
		hidden1 = self.conv_block(hidden0,  [weights['conv1'], weights['conv1_2'], weights['conv1_3']], [weights['b1'], weights['b1_2'], weights['b1_3']], 1, 'conv0', training)
		hidden2 = self.conv_block(hidden1,  [weights['conv2'], weights['conv2_2'], weights['conv2_3']], [weights['b2'], weights['b2_2'], weights['b2_3']], 2, 'conv1', training)
		hidden3 = self.conv_block(hidden2,  [weights['conv3'], weights['conv3_2'], weights['conv3_3']], [weights['b3'], weights['b3_2'], weights['b3_3']], 2, 'conv2', training)
		hidden4 = self.conv_block(hidden3,  [weights['conv4'], weights['conv4_2'], weights['conv4_3']], [weights['b4'], weights['b4_2'], weights['b4_3']], 2, 'conv3', training)
		# print(hidden4.shape)
		# hidden4 = tf.nn.max_pool(hidden4, [1, 2, 2, 1], [1, 2, 2, 1], 'SAME', name='conv3' + '_pool')

		hidden4 = tf.reduce_mean(hidden4, axis=[1,2])
		# get_shape is static shape, (5, 5, 5, 32)
		# print('flatten:', hidden4.get_shape())
		# flatten layer
		hidden4 = tf.reshape(hidden4, [-1, np.prod([int(dim) for dim in hidden4.get_shape()[1:]])], name='reshape2')

		output = tf.add(tf.matmul(hidden4, weights['w5']), weights['b5'], name='fc1')

		return output

max_stars_repo_name: ncble/MAML

max_issues_repo_name: ncble/MAML

max_forks_repo_name: ncble/MAML


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Time    : 18-5-21 上午11:38
# @Author  : Luo Yao
# @Site    : http://icode.baidu.com/repos/baidu/personal-code/Luoyao
# @File    : lanenet_hnet_model.py
# @IDE: PyCharm Community Edition
"""
LaneNet中的HNet模型
"""
import tensorflow as tf

from encoder_decoder_model import cnn_basenet
from lanenet_model import lanenet_hnet_loss

class LaneNetHNet(cnn_basenet.CNNBaseModel):
    """
    实现lanenet中的hnet模型
    """
    def __init__(self, phase):
        """

        :param phase:
        """
        super(LaneNetHNet, self).__init__()
        self._train_phase = tf.constant('train', tf.string)
        self._phase = phase
        self._is_training = self._init_phase()

        return

    def _init_phase(self):
        """

        :return:
        """
        return tf.equal(self._phase, self._train_phase)

    def _conv_stage(self, inputdata, out_channel, name):
        """

        :param inputdata:
        :param out_channel:
        :param name:
        :return:
        """
        with tf.variable_scope(name):
            conv = self.conv2d(inputdata=inputdata, out_channel=out_channel, kernel_size=3, use_bias=False, name='conv')
            bn = self.layerbn(inputdata=conv, is_training=self._is_training, name='bn')
            relu = self.relu(inputdata=bn, name='relu')

        return relu

    def _build_model(self, input_tensor, name):
        """

        :param input_tensor:
        :param name:
        :return:
        """
        with tf.variable_scope(name):
            conv_stage_1 = self._conv_stage(inputdata=input_tensor, out_channel=16, name='conv_stage_1')
            conv_stage_2 = self._conv_stage(inputdata=conv_stage_1, out_channel=16, name='conv_stage_2')
            maxpool_1 = self.maxpooling(inputdata=conv_stage_2, kernel_size=2, stride=2, name='maxpool_1')
            conv_stage_3 = self._conv_stage(inputdata=maxpool_1, out_channel=32, name='conv_stage_3')
            conv_stage_4 = self._conv_stage(inputdata=conv_stage_3, out_channel=32, name='conv_stage_4')
            maxpool_2 = self.maxpooling(inputdata=conv_stage_4, kernel_size=2, stride=2, name='maxpool_2')
            conv_stage_5 = self._conv_stage(inputdata=maxpool_2, out_channel=64, name='conv_stage_5')
            conv_stage_6 = self._conv_stage(inputdata=conv_stage_5, out_channel=64, name='conv_stage_6')
            maxpool_3 = self.maxpooling(inputdata=conv_stage_6, kernel_size=2, stride=2, name='maxpool_3')
            fc = self.fullyconnect(inputdata=maxpool_3, out_dim=1024, use_bias=False, name='fc')
            fc_relu = self.relu(inputdata=fc, name='fc_relu')
            output = self.fullyconnect(inputdata=fc_relu, out_dim=6, use_bias=False, name='fc_output')
            output = self.squeeze(inputdata=output, axis=0)

        return output

    def compute_loss(self, input_tensor, gt_label_pts, name):
        """
        计算hnet损失函数
        :param input_tensor: 原始图像[n, h, w, c]
        :param gt_label_pts: 原始图像对应的标签点集[x, y, 1]
        :param name:
        :return:
        """
        with tf.variable_scope(name):
            transformation_coefficient = self._build_model(input_tensor, name='transfomation_coefficient')
            loss = lanenet_hnet_loss.hnet_loss(gt_pts=gt_label_pts,
                                               transformation_coeffcient=transformation_coefficient,
                                               name='hnet_loss')

            return loss

    def inference(self, input_tensor, name):
        """

        :param input_tensor:
        :param name:
        :return:
        """
        with tf.variable_scope(name):
            return self._build_model(input_tensor, name='transfomation_coefficient')

if __name__ == '__main__':
    tensor_in = tf.placeholder(dtype=tf.float32, shape=[1, 64, 128, 3])
    gt_label_pts = tf.placeholder(dtype=tf.float32, shape=[None, 3])

    net = LaneNetHNet(phase=tf.constant('train', tf.string))
    c_loss = net.compute_loss(tensor_in, gt_label_pts=gt_label_pts, name='hnet')

    from data_provider import lanenet_hnet_data_processor
    import numpy as np
    import cv2
    try:
        from cv2 import cv2
    except ImportError:
        pass
    train_dataset = lanenet_hnet_data_processor.DataSet(
        ['/home/baidu/DataBase/Semantic_Segmentation/TUSimple_Lane_Detection/training/label_data_0531.json'])

    with tf.Session() as sess:
        sess.run(tf.global_variables_initializer())

        image, label_pts = train_dataset.next_batch(1)
        image = [cv2.resize(tmp, (128, 64), interpolation=cv2.INTER_LINEAR) for tmp in image]
        label_pts = np.concatenate((label_pts[0], np.ones(shape=[len(label_pts[0]), 1], dtype=np.float32)), axis=1)
        loss_val = sess.run(c_loss, feed_dict={tensor_in: image, gt_label_pts: label_pts})
        print(loss_val)


text:
import numpy as np

class GaussianNB():
    '''
    A class with simple implementation of Gaussian Naive Bayes.

    Args:
        smoothing   Smoothning factor for variance to avoid divide by zero.
                    Default: 1e-9
    '''
    def __init__(self, smoothing= 1e-9):
        self.smoothing = smoothing

    def fit(self, X, y):
        self.classes_ = np.unique(y)
        num_classes = len(self.classes_)
        num_samples, num_features = X.shape

        self.mean_ = np.zeros((num_classes, num_features))
        self.var_  = np.zeros((num_classes, num_features))
        self.counts_ = np.zeros(num_classes)

        for i, c in enumerate(self.classes_):
            X_c = X[y == c]
            self.mean_[i, :] = np.mean(X_c, axis=0)
            self.var_[i, :] = np.var(X_c, axis=0)
            self.counts_[i] = X_c.shape[0]

        self.var_[:,:] += self.smoothing
        self.priors_ = self.counts_ / self.counts_.sum()

    def _calc_log_likelihood(self, X_test):
        likelihood = []

        for i in range(len(self.classes_)):
            log_prior = np.log(self.priors_[i])
            gaus = -0.5 * np.sum(np.log(2 * np.pi * self.var_[i, :]))
            gaus -= 0.5 * np.sum(((X_test - self.mean_[i, :] )** 2) / (self.var_[i, :]), axis= 1)
            likelihood.append(log_prior+gaus)

        return np.array(likelihood).T

    def predict(self, X_test):
        y_preds = self._calc_log_likelihood(X_test)
        return self.classes_[np.argmax(y_preds, axis=1)]