text:
import os
import torch
import pickle
import gensim
import numpy as np
import pandas as pd
from config import model_config as config
from gensim.scripts.glove2word2vec import glove2word2vec

from sklearn.metrics import confusion_matrix, accuracy_score, f1_score, precision_score, recall_score

import itertools
import matplotlib.pyplot as plt

def generate_word_embeddings(vocab):
    if not os.path.exists('{}gensim.glove.6B.{}d.txt'.format(
            config['embeddings_dir'], config['embedding_dim'])):
        glove2word2vec(glove_input_file='{}glove.6B.{}d.txt'.format(
            config['embeddings_dir'], config['embedding_dim']),
            word2vec_output_file='{}gensim.glove.6B.{}d.txt'.format(
            config['embeddings_dir'], config['embedding_dim']))

    embeddings_all = gensim.models.KeyedVectors.load_word2vec_format(
        '{}gensim.glove.6B.{}d.txt'.format(config['embeddings_dir'],
                                           config['embedding_dim']))
    print('Loaded original embeddings')

    # initialize word embeddings matrix
    combined_word_embeddings = np.zeros((vocab.size,
                                         config['embedding_dim']))
    for index, word in vocab.index2word.items():
        try:
            if index < 4:  # deal with special tokens
                combined_word_embeddings[index] = np.random.normal(
                    size=(config['embedding_dim'], ))
                continue
            combined_word_embeddings[index] = embeddings_all[word]
        except KeyError as e:
            print('KeyError triggered for {}'.format(word))
            combined_word_embeddings[index] = np.random.normal(
                size=(config['embedding_dim'], ))
    print('Created combined + filtered embeddings.')
    with open('{}saved_{}d_word_embeddings.pkl'.format(
            config['embeddings_dir'], config['embedding_dim']), 'wb') as f:
        pickle.dump(combined_word_embeddings, f)
    combined_word_embeddings = torch.from_numpy(combined_word_embeddings).float()
    return combined_word_embeddings

def load_word_embeddings():
    with open('{}saved_{}d_word_embeddings.pkl'.format(
            config['embeddings_dir'], config['embedding_dim']), 'rb') as f:
        combined_word_embeddings = pickle.load(f)
        return torch.from_numpy(combined_word_embeddings).float()

def zero_padding(l, fillvalue=config['<PAD>']):
    return list(itertools.zip_longest(*l, fillvalue=fillvalue))

def binary_matrix(l, value=config['<PAD>']):
    m = []
    for i, seq in enumerate(l):
        m.append([])
        for token in seq:
            if token == 0:
                m[i].append(0)
            else:
                m[i].append(1)
    return m

# Returns padded input sequence tensor and lengths
def input_var(l, vocab):
    indexes_batch = [indexes_from_sentence(vocab, sentence) for sentence in l]
    for idx, indexes in enumerate(indexes_batch):
        indexes_batch[idx] = indexes_batch[idx] + [config['<EOS>']]
    lengths = torch.tensor([len(indexes) for indexes in indexes_batch])
    pad_list = zero_padding(indexes_batch)
    pad_var = torch.LongTensor(pad_list)
    return pad_var, lengths

def indexes_from_sentence(vocab, sentence):
    indexes = []
    for word in sentence.strip().split():
        try:
            indexes.append(vocab.word2index[word])
        except KeyError as e:
            indexes.append(config['<UNK>'])
    return indexes[:config['max_sequence_length']]

def load_data(batched=True, test=False, file_dir='../../data/t2e/'):
    # Load vocab
    with open(config['vocab_path'], 'rb') as f:
        vocab = pickle.load(f)

    bs = config['batch_size']
    ftype = 'test' if test else 'train'

    df = pd.read_csv('{}text_{}.csv'.format(file_dir, ftype))
    data = (np.array(list(df['transcription'])), np.array(df['label']))

    data = list(zip(data[0], data[1]))
    data.sort(key=lambda x: len(x[0].split()), reverse=True)

    n_iters = len(data) // bs

    if test:
        input_batch = []
        output_batch = []
        for e in data:
            input_batch.append(e[0])
            output_batch.append(e[1])
        inp, lengths = input_var(input_batch, vocab)
        return [inp, lengths, torch.LongTensor(output_batch)]

    batches = []
    for i in range(1, n_iters + 1):
        input_batch = []
        output_batch = []
        for e in data[bs * (i-1):bs * i]:
            input_batch.append(e[0])
            output_batch.append(e[1])
        inp, lengths = input_var(input_batch, vocab)
        batches.append([inp, lengths,
                        torch.LongTensor(output_batch)])
    return batches

def evaluate(targets, predictions):
    performance = {
        'acc': accuracy_score(targets, predictions),
        'f1': f1_score(targets, predictions, average='macro'),
        'precision': precision_score(targets, predictions, average='macro'),
        'recall': recall_score(targets, predictions, average='macro')}
    return performance

def plot_confusion_matrix(targets, predictions, classes,
                          normalize=False,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    # plt.figure(figsize=(8,8))
    cm = confusion_matrix(targets, predictions)
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45)
    plt.yticks(tick_marks, classes)

    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    print(cm)

    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, cm[i, j],
                 horizontalalignment="center",
                 color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')

max_forks_repo_name: Demfier/speech-emotion-recognition


text:
import unittest
import numpy as np
from numpy.core.fromnumeric import shape
import numpy.testing as npt

from LogisticRegression import sigmoid, J, y_estim, log_regression, y_estim_all
from plotTools import plot_features
from NeuralNetwork import NeuralNetwork

class TestNode(unittest.TestCase):

    def test_a(self):
        self.assertEqual(1, 1, "should be 1")

    def test_sigmoid(self):
        npt.assert_almost_equal(
            np.array([0.9525741268224331, 4.539786870243442e-05, 0.5]),
            sigmoid(np.array([3, -10, 0]))
        )

    def test_J(self):
        npt.assert_equal(
            17,
            J(np.array([1, 3]), np.array([4, -2]))
        )

    def test_y_estim(self):
        npt.assert_almost_equal(
            0.549833997312478,
            y_estim(np.array([1, 2, 3]), np.array([4, 1, -3]), 3.2)
        )

if __name__ == '__main__':

    # w = np.array([2, 1])
    # b = -15
    # x_data = 10*np.random.rand(100, 2)
    # y_data = np.zeros(shape=(100, 1))
    # for i in range(x_data.shape[0]):
    #     i_data = x_data[i]
    #     if np.dot(i_data, w) + b > 0:
    #         y_data[i][0] = 1

    # w_hat, b_hat = log_regression(x_data, y_data)
    # print('w_hat: ', w_hat)
    # print('b_hat: ', b_hat)

    # y_hat = np.zeros(shape=(100, 1))

    # for i in range(x_data.shape[0]):
    #     i_data = x_data[i]
    #     if sigmoid(np.dot(i_data, w_hat) + b_hat) > 0.5:
    #         y_hat[i][0] = 1

    # print('error: ', np.sum((y_hat - y_data)**2))

    nn = NeuralNetwork(alpha=1, nh=(3, 5, 2), tmax=10)
    print(nn)

    print(nn.forward_propagation(np.array([[[1], [2], [3]], [[1], [2], [3]], [[1], [2], [3]], [[1], [2], [3]]])))

    print(nn.backward_propagation(np.array([[[1], [2], [3]], [[1], [2], [3]], [[1], [2], [3]], [[1], [2], [3]]])))

    # unittest.main()

ext: py

max_stars_repo_name: troncosoae/myNN

max_issues_repo_name: troncosoae/myNN

max_forks_repo_name: troncosoae/myNN

reason: import numpy,from numpy


text:
import warnings
warnings.filterwarnings("ignore")

import copy
import math
from collections import defaultdict
from contextlib import suppress

import toolz
import torch
import numpy as np
import pandas as pd

import pinot
from pinot.generative import SemiSupervisedNet
from pinot.metrics import _independent
from pinot.active.acquisition import thompson_sampling, _pi, _ei_analytical, _ucb
from pinot.active.biophysical_acquisition import (biophysical_thompson_sampling,
                                                  _sample_and_marginalize_delta_G)

######################
# Utilities

def _get_thompson_values(net, data, q=1):
    """ Gets the value associated with the Thompson Samples.
    """
    # unpack data
    gs, _ = data

    # get predictive posterior
    distribution = _independent(net.condition(gs))

    # obtain samples from posterior
    thetas = distribution.sample((q,)).detach()

    # find the max of each sample
    thompson_values = torch.max(thetas, axis=1).values
    return thompson_values

def _get_biophysical_thompson_values(
    net,
    data,
    acq_func,
    q=1,
    concentration=20,
    dG_samples=10):
    """ Generates m Thompson samples and maximizes them.
    """
    # fill batch
    thompson_values = []
    for _ in range(q):

        # for each thompson sample,
        # get max, marginalizing across all distributions
        thompson_value = _sample_and_marginalize_delta_G(
            net,
            data,
            concentration=concentration,
            dG_samples=dG_samples,
            n_samples=1
        ).max().item()
        thompson_values.append(thompson_value)

    # convert to tensor
    thompson_values = torch.Tensor(thompson_values)
    return thompson_values

def thompson_sample(net, data, unseen, num_samples=1, **kwargs):
    """ Perform retrospective and prospective Thompson Sampling
        to check model beliefs about y_max.

    Returns
    -------
    thompson_samples : dict
        Dictionary with keys 'prospective' and 'retrospective'
        Values are tensors.
    """
    def _ts(net, data, num_samples=1):
        """ Get Thompson samples.
        """
        if isinstance(net.output_regressor, pinot.regressors.BiophysicalRegressor):
            ts_values = _get_biophysical_thompson_values(net, data, q=num_samples)
        else:
            ts_values = _get_thompson_values(net, data, q=num_samples)
        return ts_values

    # save log sigma
    sigma = torch.exp(net.output_regressor.log_sigma)

    # construct measurement error distribution
    epsilon = torch.distributions.Normal(0., sigma)

    # make predictive posterior noiseless by setting sigma to 0
    zero = torch.tensor(-10.)
    if torch.cuda.is_available():
        zero = zero.cuda()
    net.output_regressor.log_sigma = torch.nn.Parameter(zero)

    # set net to eval
    net.eval()

    # initialize thompson sampling dict
    thompson_samples = {}

    # get prospective thompson samples
    if unseen:

        # prospective evaluates only unseen data
        unseen_data = pinot.active.experiment._slice_fn_tuple(data, unseen)
        thompson_samples['prospective'] = _ts(
            net,
            unseen_data,
            num_samples=num_samples
        ) + epsilon.sample((num_samples,))

    # get retrospective thompson samples on all data
    thompson_samples['retrospective'] = _ts(
        net,
        data,
        num_samples=num_samples
    ) + epsilon.sample((num_samples, ))

    return thompson_samples

def _max_utility(net, data, unseen, utility_func):
    """ Finds max of the beliefs of a network according to some utility function
    """
    # unpack data
    gs, ys = data

    # get y_best
    seen = [s for s in range(len(ys)) if s not in unseen]
    y_best_round = ys[seen].max().item()
    y_best_global = ys.max().item()

    # set net to eval
    net.eval()

    # get predictive posterior
    distribution = _independent(net.condition(gs))

    # initialize beliefs dict
    beliefs = {}

    # get prospective beliefs
    if unseen:

        # prospective evaluates only unseen data
        unseen_data = pinot.active.experiment._slice_fn_tuple(data, unseen)
        beliefs['prospective'] = torch.max(
            utility_func(
                distribution, y_best=y_best_round,
            )
        ).detach().unsqueeze(0)

    # get retrospective thompson samples on all data
    beliefs['retrospective'] = torch.max(
        utility_func(
            distribution, y_best=y_best_global,
        )
    ).detach().unsqueeze(0)

    return beliefs

def max_probability_of_improvement(net, data, unseen, **kwargs):
    """ Computes the belief about the max using PI utility function.
    """
    beliefs = _max_utility(net, data, unseen, _pi)
    return beliefs

def max_upper_confidence_bound(net, data, unseen, **kwargs):
    """ Computes the belief about the max using UCB utility function.
    """
    beliefs = _max_utility(net, data, unseen, _ucb)
    return beliefs

def max_expected_improvement(net, data, unseen, **kwargs):
    """ Computes the belief about the max using EI utility function.
    """
    beliefs = _max_utility(net, data, unseen, _ei_analytical)
    return beliefs

######################
# Function definitions

class BeliefActivePlot():

    def __init__(self, net, config,
                 lr, optimizer_type,
                 data, acquisition, num_samples, num_thompson_samples, q,
                 beliefs,
                 device, num_trials, num_rounds, num_epochs):

        # net config
        self.net = net
        self.config = config

        # optimizer config
        self.lr = lr
        self.optimizer_type = optimizer_type

        # experiment config
        self.data = data
        self.acquisition = acquisition
        self.num_samples = num_samples
        self.num_thompson_samples = num_thompson_samples
        self.q = q
        self.train = pinot.app.experiment.Train

        # beliefs
        self.belief_functions = beliefs
        self.belief_functions_dict = {
            'ThompsonSampling': thompson_sample,
            'MaxProbabilityOfImprovement': max_probability_of_improvement,
            'MaxUCB': max_upper_confidence_bound,
            'MaxExpectedImprovement': max_expected_improvement,
        }

        # housekeeping
        self.device = torch.device(device)
        self.num_trials = num_trials
        self.num_rounds = num_rounds
        self.num_epochs = num_epochs

        # recording data
        self.results = []
        self.prospective_beliefs = []
        self.retrospective_beliefs = []

    def generate(self):
        """
        Performs experiment loops.
        """
        ds = self.generate_data()

        # get results for each trial
        final_results, prospective_beliefs, retrospective_beliefs = self.run_trials(ds)

        # create pandas dataframe to play nice with seaborn
        best_df = pd.DataFrame(final_results)
        pro_df = pd.DataFrame(prospective_beliefs)
        retro_df = pd.DataFrame(retrospective_beliefs)

        return best_df, pro_df, retro_df

    def run_trials(self, ds):
        """
        Plot the results of an active training loop

        Parameters
        ----------
        results : defaultdict of dict
            Empty dict in which to store results.
        ds : list of tuple
            Output of `pinot.data` and batched. Contains DGLGraph gs and Tensor ys.
        num_trials : int
            number of times to run each acquisition function
        limit : int
            Number of runs of bayesian optimization.
        Returns
        -------
        results
        """
        # get the real solution
        ds = self.generate_data()[0]
        acq_fn = self.get_acquisition()

        for i in range(self.num_trials):

            print(i)

            # make fresh net and optimizer
            net = self.get_net().to(self.device)

            optimizer = pinot.app.utils.optimizer_translation(
                opt_string=self.optimizer_type,
                lr=self.lr,
                weight_decay=0.01,
                kl_loss_scaling=1.0/float(len(ds[1]))
            )

            # instantiate experiment
            self.bo = pinot.active.experiment.BayesOptExperiment(
                net=net,
                data=ds,
                optimizer=optimizer,
                acquisition=acq_fn,
                num_epochs=self.num_epochs,
                early_stopping=False,
                q=self.q,
                slice_fn=pinot.active.experiment._slice_fn_tuple,
                collate_fn=pinot.active.experiment._collate_fn_graph,
                train_class=self.train
            )

            # run experiment
            acquisitions = self.bo.run(num_rounds=self.num_rounds)

            # process acquisitions into pandas-friendly format
            self.results = self.process_results(acquisitions, ds, i)

            # get beliefs
            self.beliefs = self.get_beliefs(net, ds, acquisitions, methods=self.belief_functions)

            # generate long-form records for pandas
            self.prospective_beliefs.extend(
                self.process_beliefs(self.beliefs['prospective'], i)
            )

            # generate long-form records for pandas
            self.retrospective_beliefs.extend(
                self.process_beliefs(self.beliefs['retrospective'], i)
            )

        return self.results, self.prospective_beliefs, self.retrospective_beliefs

    def get_beliefs(self, net, ds, acquisitions, methods=['ThompsonSampling']):
        """ Gets Thompson samples at each round

        Parameters
        ----------
        net : pinot.Net
            The network used for Bayesian optimization

        ds : pinot.data.Dataset
            The dataset used

        acquisitions : list of dict
            Choice of acquired and unacquired points at each round

        methods : list of str
            The method used to estimate model beliefs

        Returns
        -------
        beliefs : dict
            Values are 'prospective' and 'retrospective'
            Keys are

        """
        # loop through rounds
        round_beliefs = []
        for idx, state in self.bo.states.items():

            # unpack acquisitions
            seen, unseen = acquisitions[idx]

            # load network states
            net.load_state_dict(state)

            # set _x_tr and _y_tr if exact GP
            if isinstance(net.output_regressor, pinot.regressors.ExactGaussianProcessRegressor):
                seen_gs, seen_ys = pinot.active.experiment._slice_fn_tuple(ds, seen)
                net.g_last, net.y_last = seen_gs, seen_ys

            # gather pro and retro beliefs for each belief function
            for belief_name, belief_function in self.belief_functions_dict.items():
                if belief_name in methods:
                    round_belief = belief_function(
                        net, ds, unseen,
                        num_samples=self.num_thompson_samples
                    )
                    round_beliefs.append({belief_name: round_belief})

        # convert from list of dicts to dict of lists
        beliefs = {'prospective': defaultdict(list), 'retrospective': defaultdict(list)}
        for round_record in round_beliefs:
            for m in methods:
                for direction in beliefs.keys():
                    with suppress(KeyError):
                        belief = round_record[m][direction]
                        beliefs[direction][m].append(belief)
        return beliefs

    def process_beliefs(self, beliefs_dict, i):
        """ Processes the output of BayesOptExperiment.

        Parameters
        ----------
        beliefs_dict : tuple
            Dictionary of beliefs
        i : int
            Index of loop
        methods : list of str
            The method used to estimate model beliefs

        Returns
        -------
        ts_list : dict
            Long-style belief records
        """
        # instantiate belief list
        belief_list = []

        # record results
        for belief_name, step_beliefs in beliefs_dict.items():

            for step, beliefs in enumerate(step_beliefs):

                for belief_index, belief in enumerate(beliefs):

                    belief_list.append(
                        {'Acquisition Function': self.acquisition,
                         'Belief Function': belief_name,
                         'Trial': i,
                         'Round': step,
                         'Index': belief_index,
                         'Value': belief.item()}
                    )

        return belief_list

    def process_results(self, acquisitions, ds, i):
        """ Processes the output of BayesOptExperiment.

        Parameters
        ----------
        seen : list of int
            Items chosen by BayesOptExperiment object
        ds : tuple
            Dataset object
        i : int
            Index of loop

        Returns
        -------
        self.results : defaultdict
            Processed data
        """
        # get the final acquired points
        seen, _ = acquisitions[-1]

        # get max f(x)
        gs, ys = ds
        actual_sol = torch.max(ys).item()

        # pad if experiment stopped early
        # candidates_acquired = limit + 1 because we begin with a blind pick
        results_size = len(ys) + 1

        if self.net == 'multitask':
            results_data = actual_sol*np.ones((results_size, ys.size(1)))
            output = ys[seen]
            output[torch.isnan(output)] = -np.inf
        else:
            results_data = actual_sol*np.ones(results_size)
            output = ys[seen]

        results_data[:len(seen)] = np.maximum.accumulate(output.cpu().squeeze())

        # record results
        for step, result in enumerate(results_data):
            self.results.append({'Acquisition Function': self.acquisition,
                                 'Trial': i, 'Datapoint Acquired': step,
                                 'Round': step // self.q, 'Value': result})
        return self.results

    def generate_data(self):
        """ Generate data, put on GPU if possible.
        """
        # Load data
        ds = getattr(pinot.data, self.data)()

        # Limit to first two fields of tuple
        if len(ds[0]) > 2:
            ds = [(d[0], d[1]) for d in ds]

        # Batch data
        ds = pinot.data.utils.batch(ds, len(ds), seed=None)

        # Move data to GPU
        ds = [tuple([i.to(self.device) for i in ds[0]])]

        return ds

    def get_acquisition(self):
        """ Retrieve acquisition function and prepare for BO Experiment
        """
        acquisitions = {
            'ExpectedImprovement': pinot.active.acquisition.expected_improvement_analytical,
            'ProbabilityOfImprovement': pinot.active.acquisition.probability_of_improvement,
            'UpperConfidenceBound': pinot.active.acquisition.upper_confidence_bound,
            'Uncertainty': pinot.active.acquisition.uncertainty,
            'Human': pinot.active.acquisition.temporal,
            'Random': pinot.active.acquisition.random,
            'ThompsonSampling': pinot.active.acquisition.thompson_sampling,
        }

        return acquisitions[self.acquisition]

    def get_net(self):
        """
        Retrive GP using representation provided in args.
        """
        representation = pinot.representation.SequentialMix(config=self.config)

        if hasattr(pinot.regressors, self.net):
            output_regressor = getattr(pinot.regressors, self.net)
            net = pinot.Net(
                representation=representation,
                output_regressor_class=output_regressor,
            )

        if self.net == 'semi':
            output_regressor = pinot.regressors.NeuralNetworkRegressor
            net = SemiSupervisedNet(
                representation=representation,
                output_regressor=output_regressor
            )

        elif self.net == 'semi_gp':
            output_regressor = pinot.regressors.ExactGaussianProcessRegressor

            net = SemiSupervisedNet(
                representation=representation,
                output_regressor=output_regressor,
            )

        elif self.net == 'multitask':
            output_regressor = pinot.regressors.ExactGaussianProcessRegressor
            net = pinot.multitask.MultitaskNet(
                representation=representation,
                output_regressor=output_regressor,
            )

        return net

if __name__ == '__main__':

    # Running functions
    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument('--net', type=str, default='ExactGaussianProcessRegressor')
    parser.add_argument('--config', nargs='+', type=str,
        default=['GraphConv', '32', 'activation', 'tanh',
                 'GraphConv', '32', 'activation', 'tanh',
                 'GraphConv', '32', 'activation', 'tanh'])

    parser.add_argument('--lr', type=float, default=1e-4)
    parser.add_argument('--optimizer', type=str, default='Adam')

    parser.add_argument('--data', type=str, default='esol')
    parser.add_argument('--acquisition', type=str, default='ExpectedImprovement')
    parser.add_argument('--num_samples', type=int, default=1000)
    parser.add_argument('--q', type=int, default=1)

    parser.add_argument('--beliefs', nargs='+', type=str,
        default=['ThompsonSampling', 'MaxProbabilityOfImprovement', 'MaxUCB', 'MaxExpectedImprovement']
    )

    parser.add_argument('--device', type=str, default='cuda:0')
    parser.add_argument('--num_trials', type=int, default=1)
    parser.add_argument('--num_thompson_samples', type=int, default=500)
    parser.add_argument('--num_rounds', type=int, default=100)
    parser.add_argument('--num_epochs', type=int, default=10)

    parser.add_argument('--index_provided', type=bool, default=True)
    parser.add_argument('--index', type=int, default=0)

    args = parser.parse_args()

    plot = BeliefActivePlot(
        # net config
        net=args.net,
        config=args.config,

        # optimizer config
        optimizer_type=args.optimizer,
        lr=args.lr,

        # experiment config
        data=args.data,
        acquisition=args.acquisition,
        num_samples=args.num_samples,
        num_thompson_samples=args.num_thompson_samples,
        q=args.q,

        # beliefs
        beliefs=args.beliefs,

        # housekeeping
        device=args.device,
        num_trials=args.num_trials,
        num_rounds=args.num_rounds,
        num_epochs=args.num_epochs,
    )

    # run experiment
    best_df, pro_ts_df, retro_ts_df = plot.generate()

    # write to disk
    beliefs_string = '_'.join(args.beliefs)
    filename = f'{args.net}_{args.optimizer}_{args.data}_num_epochs{args.num_epochs}_{args.acquisition}_q{args.q}_beliefs{beliefs_string}_{args.index}.csv'

    best_df.to_csv(f'best_{filename}')
    pro_ts_df.to_csv(f'pro_{filename}')
    retro_ts_df.to_csv(f'retro_{filename}')

max_stars_repo_name: choderalab/pinot

max_issues_repo_name: choderalab/pinot

max_forks_repo_name: choderalab/pinot

reason: import numpy


text:
# Make frames of exoplanet orbits.

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

plotpar = {'axes.labelsize': 18,
           'font.size': 10,
           'legend.fontsize': 18,
           'xtick.labelsize': 18,
           'ytick.labelsize': 18,
           'text.usetex': True}
plt.rcParams.update(plotpar)
plt.rcParams['axes.facecolor'] = 'black'

def get_orbit(a, radius, period, frames, xcenter=0, ycenter=0,
              slow_factor=30):
    """
    nframes = snapshots per orbit.
    if period long snapshots per orbit goes up.
    if period short, snapshots per orbit goes down.
    snapshots per orbit (nframes) is proportional to the period and
    1/angular frequency.
    """
    random_phase = np.random.uniform(0, 2*np.pi)
    x = a * np.cos(np.pi/(period*slow_factor) * np.arange(nframes)
                   + random_phase)
    y = a * np.sin(np.pi/(period*slow_factor) * np.arange(nframes)
                   + random_phase)
    return x, y

if __name__ == "__main__":

    # Read in the data
    df = pd.read_csv("planets.csv", skiprows=69)

    # Select only transiting planets.
    m = []
    for i, method in enumerate(list(df.pl_discmethod.values)):
        if method == "Transit":
            m.append(i)
    df = df.iloc[np.array(m)]

    R_j = 6.991e7
    radius_ratio = R_j/6.371e6  # Convert to Earth radii
    nframes = 800
    assert nframes % 2 == 0, "nframes must be an even number"
    nstar = 800

    a_s = df.pl_orbsmax.values #* 1.496e11 / R_j
    radii = df.pl_radj.values * radius_ratio
    periods = df.pl_orbper.values
    xpositions = np.zeros((nstar, nframes))
    ypositions = np.zeros((nstar, nframes))
    for star in range(nstar):
        x, y = get_orbit(a_s[star], radii[star], periods[star], nframes)
        xpositions[star, :] = x
        ypositions[star, :] = y

    lim_min = -np.logspace(-2, np.log10(.5), nframes)
    lim_max = np.logspace(-2, np.log10(.5), nframes)
    for pos in range(nframes):
        plt.clf()
        plt.figure(figsize=(16, 9))
        for star in range(nstar):
            plt.scatter(xpositions[star, pos], ypositions[star, pos],
                        s=radii[star]*1./(lim_max[pos]-lim_min[pos]),
                        c="w", alpha=.5)
        plt.xlim(lim_min[pos], lim_max[pos])
        plt.ylim(lim_min[pos], lim_max[pos])
        plt.gca().set_aspect('equal', adjustable='box')
        plt.savefig("planet_cloud_movie/frame_{}".format(str(pos).zfill(4)))
        plt.close()

    # Make the movie file
    import os
    # framerate = 30
    framerate = 20
    quality = 25
    os.system("/Applications/ffmpeg -r {0} -f image2 -s 1920x1080 -i "\
            "planet_cloud_movie/frame_%04d.png -vcodec libx264 -crf {1} "\
             "-pix_fmt yuv420p planet_cloud_movie.mp4".format(framerate,
                                                              quality))

max_stars_repo_name: RuthAngus/astro-viz

max_issues_repo_name: RuthAngus/astro-viz

max_forks_repo_name: RuthAngus/astro-viz


text:
import levitate.arrays
import levitate.hardware
import numpy as np
import pytest

# Tests created with these air properties
from levitate.materials import air
air.c = 343
air.rho = 1.2

# Tests were mostly written before the layout was transposed,
# so a lot of the "expected" positions etc are transposed in the test hardcoding.

def test_rectangular_grid():
    # positions, normals = levitate.arrays.rectangular_grid(shape=(5, 3), spread=10e-3)
    array = levitate.arrays.RectangularArray(shape=(5, 3), transducer_size=10e-3)
    expected_positions = np.array([
        [-0.02, -0.01, 0.],
        [-0.01, -0.01, 0.],
        [0.   , -0.01, 0.],
        [ 0.01, -0.01, 0.],
        [ 0.02, -0.01, 0.],
        [-0.02,  0.  , 0.],
        [-0.01,  0.  , 0.],
        [ 0.  ,  0.  , 0.],
        [ 0.01,  0.  , 0.],
        [ 0.02,  0.  , 0.],
        [-0.02,  0.01, 0.],
        [-0.01,  0.01, 0.],
        [ 0.  ,  0.01, 0.],
        [ 0.01,  0.01, 0.],
        [ 0.02,  0.01, 0.]]).T
    expected_normals = np.array([
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.]]).T
    np.testing.assert_allclose(array.positions, expected_positions)
    np.testing.assert_allclose(array.normals, expected_normals)

@pytest.mark.parametrize('rectangular_class, shape, spread', [
    (levitate.hardware.DragonflyArray, (16, 16), 10.47e-3),
    (levitate.hardware.AcoustophoreticBoard, (16, 16), 10.47e-3),
])
def test_rectangular_variants(rectangular_class, shape, spread):
    # Check that we do place all the indices.
    all_indices = np.sort(rectangular_class.grid_indices.flatten())
    np.testing.assert_equal(all_indices, list(range(np.prod(shape))))
    array = rectangular_class()
    # Test that every element is where it should be
    for idx, (x, y, _) in enumerate(zip(*array.positions)):
        y_idx, x_idx = np.argwhere(array.grid_indices == idx)[0]
        np.testing.assert_allclose(x / array.spread + 7.5, x_idx)
        np.testing.assert_allclose(7.5 - y / array.spread, y_idx)
    # Test that this is indeed a 16x16 rectangular array with spread 10.47 mm
    rectangular = levitate.arrays.RectangularArray(shape=shape, spread=spread)
    differences = rectangular.positions[:, :, None] - array.positions[:, None, :]
    distances = np.sum(differences**2, axis=0)
    assert all([(np.min(distances, axis=axis) == 0).all() for axis in range(distances.ndim)])

def test_array_offset():
    array = levitate.arrays.RectangularArray(shape=(4, 2), offset=(0.1, -0.2, 1.4))
    expected_positions = np.array([
        [0.085, -0.205, 1.4],
        [0.095, -0.205, 1.4],
        [0.105, -0.205, 1.4],
        [0.115, -0.205, 1.4],
        [0.085, -0.195, 1.4],
        [0.095, -0.195, 1.4],
        [0.105, -0.195, 1.4],
        [0.115, -0.195, 1.4]]).T
    expected_normals = np.array([
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.],
        [0., 0., 1.]]).T
    np.testing.assert_allclose(array.positions, expected_positions)
    np.testing.assert_allclose(array.normals, expected_normals)

@pytest.mark.parametrize('offset', [(0., 0., 0.), (1.85194e-3, 28.2839e-3, 35.2830e-3)])
@pytest.mark.parametrize('radius', [0.4, 0.05])
@pytest.mark.parametrize('rings', [4, 15])
@pytest.mark.parametrize('normal', [(0., 0., 1.), (-2., 5., 4.37)])
@pytest.mark.parametrize('packing', ['distance', 'count'])
@pytest.mark.parametrize('spread', [10e-3, 5e-3])
def test_spherical_cap(radius, rings, normal, offset, packing, spread):
    normal, offset = np.asarray(normal), np.asarray(offset)
    array = levitate.arrays.SphericalCapArray(radius=radius, rings=rings, spread=spread, packing=packing, normal=normal, offset=offset)
    focus = normal / np.sum(normal**2)**0.5 * radius + offset
    diff = focus[:, None] - array.positions
    dist = np.sum(diff**2, 0)**0.5
    np.testing.assert_allclose(dist, radius)
    np.testing.assert_allclose(np.sum(diff * array.normals, 0) / dist, 1.)

    array = levitate.arrays.DoublesidedArray(
        levitate.arrays.SphericalCapArray(
            radius=radius, rings=rings, spread=spread, packing=packing, normal=normal
        ),
        separation=radius * 2, normal=normal, offset=offset
    )
    diff = offset[:, None] - array.positions
    dist = np.sum(diff**2, 0)**0.5
    np.testing.assert_allclose(dist, radius)
    np.testing.assert_allclose(np.sum(diff * array.normals, 0) / dist, 1.)

def test_array_normal():
    array = levitate.arrays.RectangularArray(shape=(4, 2), normal=(2, 3, 4))
    expected_positions = np.array([
        [-1.321925551875e-02, -2.328883278124e-03, +8.356290217967e-03],
        [-4.010697510416e-03, -3.516046265624e-03, +4.642383454426e-03],
        [+5.197860497917e-03, -4.703209253125e-03, +9.284766908853e-04],
        [+1.440641850625e-02, -5.890372240625e-03, -2.785430072656e-03],
        [-1.440641850625e-02, +5.890372240625e-03, +2.785430072656e-03],
        [-5.197860497917e-03, +4.703209253125e-03, -9.284766908853e-04],
        [+4.010697510416e-03, +3.516046265624e-03, -4.642383454426e-03],
        [+1.321925551875e-02, +2.328883278124e-03, -8.356290217967e-03]]).T
    expected_normals = np.array([
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01],
        [+3.713906763541e-01, +5.570860145312e-01, +7.427813527082e-01]]).T
    np.testing.assert_allclose(array.positions, expected_positions)
    np.testing.assert_allclose(array.normals, expected_normals)

def test_array_rotation():
    array = levitate.arrays.RectangularArray(shape=(4, 2), rotation=1, normal=(-1, 4, -2))
    expected_positions = np.array([
        [-8.747019947780e-03, +4.075784833732e-03, +1.252507964135e-02],
        [-9.564871486238e-04, +2.940455719910e-03, +6.359155014132e-03],
        [+6.834045650533e-03, +1.805126606089e-03, +1.932303869110e-04],
        [+1.462457844969e-02, +6.697974922671e-04, -5.972694240310e-03],
        [-1.462457844969e-02, -6.697974922671e-04, +5.972694240310e-03],
        [-6.834045650533e-03, -1.805126606089e-03, -1.932303869110e-04],
        [+9.564871486238e-04, -2.940455719910e-03, -6.359155014132e-03],
        [+8.747019947780e-03, -4.075784833732e-03, -1.252507964135e-02]]).T
    expected_normals = np.array([
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01],
        [-2.182178902360e-01, +8.728715609440e-01, -4.364357804720e-01]]).T
    np.testing.assert_allclose(array.positions, expected_positions)
    np.testing.assert_allclose(array.normals, expected_normals)

def test_double_sided_grid():
    array = levitate.arrays.DoublesidedArray(levitate.arrays.RectangularArray, separation=0.5, shape=(2, 2), spread=5e-3, offset=(0.1, -0.2, 1.4), normal=(0.2, -1.4, 2), rotation=1)

    expected_positions = np.array([
        [+8.02489979e-02, -5.99269663e-02, +1.19384001e+00],
        [+8.30486836e-02, -5.64068310e-02, +1.19602413e+00],
        [+7.61264873e-02, -5.78193656e-02, +1.19572758e+00],
        [+7.89261731e-02, -5.42992303e-02, +1.19791170e+00],
        [+1.21073827e-01, -3.45700770e-01, +1.60208830e+00],
        [+1.23873513e-01, -3.42180634e-01, +1.60427242e+00],
        [+1.16951316e-01, -3.43593169e-01, +1.60397587e+00],
        [+1.19751002e-01, -3.40073034e-01, +1.60615999e+00]]).T
    expected_normals = np.array([
        [+8.16496581e-02, -5.71547607e-01, +8.16496581e-01],
        [+8.16496581e-02, -5.71547607e-01, +8.16496581e-01],
        [+8.16496581e-02, -5.71547607e-01, +8.16496581e-01],
        [+8.16496581e-02, -5.71547607e-01, +8.16496581e-01],
        [-8.16496581e-02, +5.71547607e-01, -8.16496581e-01],
        [-8.16496581e-02, +5.71547607e-01, -8.16496581e-01],
        [-8.16496581e-02, +5.71547607e-01, -8.16496581e-01],
        [-8.16496581e-02, +5.71547607e-01, -8.16496581e-01]]).T
    np.testing.assert_allclose(array.positions, expected_positions)
    np.testing.assert_allclose(array.normals, expected_normals)
    expected_signature = np.array([0, 0, 0, 0, np.pi, np.pi, np.pi, np.pi])
    np.testing.assert_allclose(array.signature(stype='doublesided'), expected_signature)

def test_Array_basics():
    pos, norm = np.random.normal(size=(2, 3, 8))
    array = levitate.arrays.TransducerArray(pos, norm)
    array.omega = 200000
    np.testing.assert_allclose(2 * np.pi * array.freq, array.omega)
    array.k = 730
    np.testing.assert_allclose(2 * np.pi / array.k, array.wavelength)
    array.wavelength = 8.5e-3
    np.testing.assert_allclose(air.c, array.wavelength * array.freq)
    array.freq = 41e3
    np.testing.assert_allclose(array.transducer.freq, 41e3)

    from levitate.transducers import PlaneWaveTransducer
    array = levitate.arrays.RectangularArray(shape=(4, 4), transducer=PlaneWaveTransducer)
    pos = np.array([7e-3, -3e-3, 70e-3])
    np.testing.assert_allclose(array.focus_phases(pos), np.array([+2.069588782645e+00, -2.511641902621e+00, -1.794667262208e+00, -2.103144086214e+00, +2.763870843565e+00, -1.794667262208e+00, -1.067679552580e+00, -1.380498791671e+00, +2.465230031099e+00, -2.103144086214e+00, -1.380498791671e+00, -1.691434660688e+00, +1.189733359830e+00, +2.863786956082e+00, -2.714709695304e+00, -3.017860322253e+00]))
    np.testing.assert_allclose(array.signature(pos, stype='twin'), np.array([-1.570796326795e+00, -1.570796326795e+00, +1.570796326795e+00, +1.570796326795e+00, -1.570796326795e+00, -1.570796326795e+00, -1.570796326795e+00, +1.570796326795e+00, -1.570796326795e+00, -1.570796326795e+00, -1.570796326795e+00, +1.570796326795e+00, -1.570796326795e+00, -1.570796326795e+00, -1.570796326795e+00, +1.570796326795e+00]))
    np.testing.assert_allclose(array.signature(pos, stype='vortex'), np.array([-2.642245931910e+00, -2.356194490192e+00, -1.735945004210e+00, -9.827937232473e-01, -3.050932766389e+00, -2.976443976175e+00, -2.356194490192e+00, -2.449786631269e-01, +2.792821650006e+00, +2.553590050042e+00, +1.815774989922e+00, +7.853981633974e-01, +2.455863142684e+00, +2.158798930342e+00, +1.681453547969e+00, +1.152571997216e+00]))
    np.testing.assert_allclose(array.signature(pos, stype='bottle'), np.array([+3.141592653590e+00, +3.141592653590e+00, +0.000000000000e+00, +0.000000000000e+00, +3.141592653590e+00, +0.000000000000e+00, +0.000000000000e+00, +0.000000000000e+00, +3.141592653590e+00, +0.000000000000e+00, +0.000000000000e+00, +0.000000000000e+00, +3.141592653590e+00, +3.141592653590e+00, +3.141592653590e+00, +3.141592653590e+00]))
    pos = np.array([0.1, -0.2, 0.3])
    np.testing.assert_allclose(array.focus_phases(pos), np.array([-1.4782875, 0.70451472, 2.70433793, -1.76613547, 1.07535199, -3.05482743, -1.08272084, 0.70451472, -2.79668059, -0.67375749, 1.27038374, 3.03192953, -0.52217161, 1.57048339, -2.79668059, -1.0609514]))
    np.testing.assert_allclose(array.signature(pos, stype='twin'), np.array([-1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633]))
    np.testing.assert_allclose(array.signature(stype='vortex'), np.array([-2.35619449, -1.89254688, -1.24904577, -0.78539816, -2.8198421, -2.35619449, -0.78539816, -0.32175055, 2.8198421, 2.35619449, 0.78539816, 0.32175055, 2.35619449, 1.89254688, 1.24904577, 0.78539816]))
    np.testing.assert_allclose(array.signature(stype='bottle'), np.array([3.14159265, 0., 0., 3.14159265, 0., 0., 0., 0., 0., 0., 0., 0., 3.14159265, 0., 0., 3.14159265]))

def test_Array_visualizer():
    array = levitate.arrays.RectangularArray(shape=2)
    array.visualize()
    pos = np.array([0, 0, 0.05])
    signature = array.signature(angle=np.pi, stype='twin')
    phase = array.focus_phases(pos) + signature
    amps = levitate.utils.complex(phase)
    from levitate.utils import find_trap
    trap_pos = find_trap(array, pos, amps)
    np.testing.assert_allclose(pos, trap_pos, atol=0.1e-3)
    np.testing.assert_allclose(signature, array.signature(pos, phase))

def test_Array_calculations():
    array = levitate.arrays.RectangularArray(shape=2)
    pos = np.array([0.1, -0.2, 0.3])
    expected_result = np.array([
        [-1.60298320e+01+1.39433272e+00j,  7.60005836e+00+1.43149212e+01j, 1.23972894e+01+9.89785204e+00j,  4.72749101e+00-1.52603890e+01j],
        [-2.75579656e+02-3.30839037e+03j, -2.69691942e+03+1.41912500e+03j, -2.02250900e+03+2.51457673e+03j,  2.82526527e+03+8.86498882e+02j],
        [ 5.11790789e+02+6.14415354e+03j,  5.53578197e+03-2.91294080e+03j, 3.94870805e+03-4.90941170e+03j, -6.09662505e+03-1.91297127e+03j],
        [-7.87370445e+02-9.45254391e+03j, -8.51658764e+03+4.48144738e+03j, -5.77859714e+03+7.18450493e+03j,  8.92189032e+03+2.79947015e+03j],
        [ 6.80382495e+05-8.33698641e+04j, -2.91485806e+05-4.94127262e+05j, -5.27798435e+05-3.91154134e+05j, -1.38372830e+05+5.31781607e+05j],
        [ 2.35305285e+06-2.10377336e+05j, -1.13689503e+06-2.12990534e+06j, -1.95769683e+06-1.55833572e+06j, -7.53077687e+05+2.44212524e+06j],
        [ 5.57294329e+06-4.54866720e+05j, -2.65207242e+06-5.06161450e+06j, -4.17057990e+06-3.36464336e+06j, -1.64672816e+06+5.21935967e+06j],
        [-1.26844169e+06+9.63139996e+04j,  5.40041580e+05+1.04492376e+06j, 9.92856867e+05+8.10438183e+05j,  3.62769003e+05-1.12739219e+06j],
        [ 1.95144875e+06-1.48175384e+05j, -8.30833200e+05-1.60757501e+06j, -1.45296127e+06-1.18600710e+06j, -5.30881467e+05+1.64984223e+06j],
        [-3.62411911e+06+2.75182856e+05j,  1.70539446e+06+3.29975923e+06j, 2.83673391e+06+2.31553766e+06j,  1.14558632e+06-3.56018587e+06j],
        [ 2.86670148e+07+1.39550752e+08j,  8.79550534e+07-6.44855683e+07j, 7.10265677e+07-1.14387829e+08j, -1.01799225e+08-1.57149418e+07j],
        [-9.48653965e+07-9.00596323e+08j, -8.15465047e+08+4.51662413e+08j, -6.08347649e+08+7.86174017e+08j,  9.81016375e+08+2.88086328e+08j],
        [ 2.69225703e+08+3.28515432e+09j,  3.00523878e+09-1.57526192e+09j, 1.95396625e+09-2.42520390e+09j, -3.05538963e+09-9.62345944e+08j],
        [-2.90779477e+07-2.61000235e+08j, -1.91908617e+08+1.10366719e+08j, -1.57582477e+08+2.07891701e+08j,  2.12034769e+08+5.76457627e+07j],
        [ 4.47353042e+07+4.01538824e+08j,  2.95244026e+08-1.69794952e+08j, 2.30608503e+08-3.04231758e+08j, -3.10294784e+08-8.43596528e+07j],
        [ 3.80717090e+07+4.85924318e+08j,  4.02816731e+08-2.09323496e+08j, 3.21279107e+08-3.94767783e+08j, -4.51078135e+08-1.44502368e+08j],
        [ 1.08776311e+08+1.38835519e+09j,  1.27205284e+09-6.61021567e+08j, 9.17940305e+08-1.12790795e+09j, -1.42445727e+09-4.56323268e+08j],
        [ 8.12193376e+07+1.15079185e+09j,  9.57197835e+08-4.88115774e+08j, 6.93574595e+08-8.40914649e+08j, -9.64000841e+08-3.15741830e+08j],
        [-1.50835913e+08-2.13718486e+09j, -1.96477450e+09+1.00192185e+09j, -1.35412183e+09+1.64178574e+09j,  2.08021234e+09+6.81337634e+08j],
        [-4.85644305e+07-7.48335747e+08j, -6.23978731e+08+3.13792174e+08j, -4.77252163e+08+5.71923550e+08j,  6.57524678e+08+2.19514891e+08j]])
    np.testing.assert_allclose(array.pressure_derivs(pos), expected_result)

    expected_result = np.array([
        [-5.682427477075e+01 + 4.942780801222e+00j, +2.694150540810e+01 + 5.074507426259e+01j, +4.394724655282e+01 + 3.508697193457e+01j, +1.675851928357e+01 - 5.409667034377e+01j],
        [-3.803836147518e+01 - 1.657048080374e+01j, +1.279975203258e+00 + 4.120938631452e+01j, +1.710554711045e+01 + 3.829641433944e+01j, +2.807511601090e+01 - 3.087116808379e+01j],
        [-6.597796686300e+00 - 7.920790438223e+01j, -7.136502785396e+01 + 3.755243658313e+01j, -4.842194588566e+01 + 6.020279674447e+01j, +7.476127507154e+01 + 2.345825276365e+01j],
        [-3.477261873074e+01 + 2.263543161484e+01j, +3.323968212216e+01 + 2.439211736111e+01j, +4.107318751625e+01 + 8.497550262354e+00j, -5.405546885483e+00 - 4.137657816630e+01j],
        [-1.686407000850e+01 - 2.178756862198e+01j, -1.125090492261e+01 + 2.455396689049e+01j, -1.726622397846e+00 + 2.850517035519e+01j, +2.594291002246e+01 - 1.071176900518e+01j],
        [+3.031007376630e+01 - 6.821025513107e+01j, -7.468896294954e+01 + 1.768761675083e+00j, -6.814162421776e+01 + 2.984827855687e+01j, +5.477558195724e+01 + 5.054666208556e+01j],
        [+5.988534535282e+01 - 4.547151978544e+00j, -2.861052213724e+01 - 5.535835638677e+01j, -4.322596009999e+01 - 3.528400690252e+01j, -1.773137889725e+01 + 5.510453754651e+01j],
        [-4.026173202706e+01 - 6.285166991374e+01j, -4.464287578055e+01 + 5.990478567180e+01j, -1.559702090250e+01 + 7.273880410627e+01j, +7.397427410032e+01 - 9.117905355749e+00j],
        [-1.338098961723e+01 + 2.408410514370e+01j, +2.653751067525e+01 + 5.024010230645e+00j, +2.758190339425e+01 - 7.400308576910e+00j, -1.482787772910e+01 - 2.383087530295e+01j],
        [-2.707354746999e+00 - 1.746661744378e+01j, -1.332639233207e+01 + 1.070479817804e+01j, -9.733075343865e+00 + 1.606625396807e+01j, +1.821285983077e+01 + 9.712939338454e-01j],
        [+4.676829722253e+01 - 3.538544878575e+01j, -5.237549857426e+01 - 2.470351811875e+01j, -5.977893107220e+01 - 4.270891236968e+00j, +2.203974234614e+01 + 5.507262673066e+01j],
        [+7.892432023664e+01 + 3.611381986030e+01j, -1.111183761150e+00 - 8.806141599513e+01j, -3.294243372319e+01 - 7.747677728413e+01j, -5.866451648153e+01 + 6.219821126638e+01j],
        [+9.287877976849e-01 + 1.431181429618e+01j, +1.559873321909e+01 - 7.844434688683e+00j, +5.510582917790e+00 - 6.603704269743e+00j, -1.080601629861e+01 - 3.607593105422e+00j],
        [+7.359464973840e+01 - 4.601177935703e+01j, -7.002975115181e+01 - 5.340301180676e+01j, -8.212235075402e+01 - 1.854115596653e+01j, +9.533619685308e+00 + 8.496618783049e+01j],
        [-5.094954128769e+01 - 2.904389528692e+01j, -1.140045085313e+01 + 5.677575520142e+01j, +1.490296577659e+01 + 5.804879605259e+01j, +5.070404281442e+01 - 3.078708333006e+01j],
        [-4.271062938279e-01 + 1.767003324465e+01j, +1.653992786475e+01 - 4.314652112229e+00j, +1.406441459177e+01 - 1.245188798563e+01j, -1.514361128043e+01 - 1.016477803538e+01j]])
    np.testing.assert_allclose(array.spherical_harmonics(pos, orders=3), expected_result)

max_stars_repo_name: AppliedAcousticsChalmers/levitate

max_issues_repo_name: AppliedAcousticsChalmers/levitate

max_forks_repo_name: AppliedAcousticsChalmers/levitate


text:
#!/usr/bin/env python
"""

"""

from __future__ import division, absolute_import
import logging
import numpy as np
import scipy as sp
import scipy.interpolate as spinterp
import time
import datetime as dt

import matplotlib.pyplot as plt
try:
    from mayavi import mlab
except ImportError:
    mlab = None
except (ValueError,RuntimeError) as e:
    mlab = None
    print('Mayavi not imported due to {}'.format(e))

try:
    plt.get_cmap('viridis')
    defmap3d = 'viridis'
except ValueError:
    defmap3d = 'jet'
    #%%
def plot3Dslice(geodata, surfs, vbounds, titlestr='', time=0, gkey=None, cmap=defmap3d,
                ax=None, fig=None, method='linear', fill_value=np.nan, view=None, units='',
                colorbar=False, outimage=False):
    """
    This function create 3-D slice image given either a surface or list of
    coordinates to slice through.
    Inputs:
    geodata - A geodata object that will be plotted in 3D
    surfs - This is a three element list. Each element can either be
        altlist - A list of the altitudes that RISR parameter slices will be taken at
        xyvecs- A list of x and y numpy arrays that have the x and y coordinates that the
                data will be interpolated over.
                ie, xyvecs=[np.linspace(-100.0,500.0), np.linspace(0.0,600.0)]
    vbounds = a list of bounds for the geodata objec's parameters. ie, vbounds=[500,2000]
    title - A string that holds for the overall image
    ax - A handle for an axis that this will be plotted on.

    Returns a mayavi image with a surface
    """

    if mlab is None:
        print('mayavi was not successfully imported')
        return

    assert geodata.coordnames.lower() == 'cartesian'

    datalocs = geodata.dataloc

    xvec = sp.unique(datalocs[:, 0])
    yvec = sp.unique(datalocs[:, 1])
    zvec = sp.unique(datalocs[:, 2])

    assert len(xvec)*len(yvec)*len(zvec) == datalocs.shape[0]

    #determine if the ordering is fortran or c style ordering
    diffcoord = sp.diff(datalocs, axis=0)

    if diffcoord[0, 1] != 0.0:
        ar_ord = 'f'
    elif diffcoord[0, 2] != 0.0:
        ar_ord = 'c'
    elif diffcoord[0, 0] != 0.0:
        if len(np.where(diffcoord[:, 1])[0]) == 0:
            ar_ord = 'f'
        elif len(np.where(diffcoord[:, 2])[0]) == 0:
            ar_ord = 'c'

    matshape = (len(yvec), len(xvec), len(zvec))
    # reshape the arrays into a matricies for plotting
    x, y, z = [sp.reshape(datalocs[:, idim], matshape, order=ar_ord) for idim in range(3)]

    if gkey is None:
        gkey = geodata.datanames()[0]
    porig = geodata.data[gkey][:, time]

    mlab.figure(fig)
    #determine if list of slices or surfaces are given

    islists = isinstance(surfs[0], list)
    if isinstance(surfs[0], np.ndarray):
        onedim = surfs[0].ndim == 1
    #get slices for each dimension out
    surflist = []
    if islists or onedim:
        p = np.reshape(porig, matshape, order=ar_ord)
        xslices = surfs[0]
        for isur in xslices:
            indx = sp.argmin(sp.absolute(isur-xvec))
            xtmp = x[:, indx]
            ytmp = y[:, indx]
            ztmp = z[:, indx]
            ptmp = p[:, indx]
            pmask = sp.zeros_like(ptmp).astype(bool)
            pmask[sp.isnan(ptmp)] = True
            surflist.append(mlab.mesh(xtmp, ytmp, ztmp, scalars=ptmp, vmin=vbounds[0],
                                      vmax=vbounds[1], colormap=cmap, mask=pmask))
            surflist[-1].module_manager.scalar_lut_manager.lut.nan_color = 0, 0, 0, 0

        yslices = surfs[1]
        for isur in yslices:
            indx = sp.argmin(sp.absolute(isur-yvec))
            xtmp = x[indx]
            ytmp = y[indx]
            ztmp = z[indx]
            ptmp = p[indx]
            pmask = sp.zeros_like(ptmp).astype(bool)
            pmask[sp.isnan(ptmp)] = True
            surflist.append(mlab.mesh(xtmp, ytmp, ztmp, scalars=ptmp, vmin=vbounds[0],
                                      vmax=vbounds[1], colormap=cmap, mask=pmask))
            surflist[-1].module_manager.scalar_lut_manager.lut.nan_color = 0, 0, 0, 0
        zslices = surfs[2]
        for isur in zslices:
            indx = sp.argmin(sp.absolute(isur-zvec))
            xtmp = x[:, :, indx]
            ytmp = y[:, :, indx]
            ztmp = z[:, :, indx]
            ptmp = p[:, :, indx]
            pmask = sp.zeros_like(ptmp).astype(bool)
            pmask[sp.isnan(ptmp)] = True
            surflist.append(mlab.mesh(xtmp, ytmp, ztmp, scalars=ptmp, vmin=vbounds[0],
                                      vmax=vbounds[1], colormap=cmap, mask=pmask))
            surflist[-1].module_manager.scalar_lut_manager.lut.nan_color = 0, 0, 0, 0
    else:
        # For a general surface.
        xtmp, ytmp, ztmp = surfs[:]
        gooddata = ~np.isnan(porig)
        curparam = porig[gooddata]
        curlocs = datalocs[gooddata]
        new_coords = np.column_stack((xtmp.flatten(), ytmp.flatten(), ztmp.flatten()))
        ptmp = spinterp.griddata(curlocs, curparam, new_coords, method, fill_value)
        pmask = sp.zeros_like(ptmp).astype(bool)
        pmask[sp.isnan(ptmp)] = True
        surflist.append(mlab.mesh(xtmp, ytmp, ztmp, scalars=ptmp, vmin=vbounds[0],
                                  vmax=vbounds[1], colormap=cmap, mask=pmask))
        surflist[-1].module_manager.scalar_lut_manager.lut.nan_color = 0, 0, 0, 0
    mlab.title(titlestr, color=(0, 0, 0))
    #mlab.outline(color=(0,0,0))
    mlab.axes(color=(0, 0, 0), x_axis_visibility=True, xlabel='x in km', y_axis_visibility=True,
              ylabel='y in km', z_axis_visibility=True, zlabel='z in km')

    mlab.orientation_axes(xlabel='x in km', ylabel='y in km', zlabel='z in km')

    if view is not None:
        # order of elevation is changed between matplotlib and mayavi
        mlab.view(view[0], view[1])
    if colorbar:
        if units == '':
            titlestr = gkey
        else:
            titlstr = gkey +' in ' +units
        mlab.colorbar(surflist[-1], title=titlstr, orientation='vertical')

    if outimage:
        arr = mlab.screenshot(fig, antialiased=True)
        mlab.close(fig)
        return arr
    else:
        return surflist

max_stars_repo_name: jswoboda/GeoDataPython

max_issues_repo_name: jswoboda/GeoDataPython

max_forks_repo_name: jswoboda/GeoDataPython

reason: import numpy,import scipy


text:
# Fully convolutional neural network, adapted from the implementation from Fawaz et. al
# https://github.com/hfawaz/dl-4-tsc/blob/master/classifiers/fcn.py
#
# Network originally proposed by:
#
# @inproceedings{wang2017time,
#   title={Time series classification from scratch with deep neural networks: A strong baseline},
#   author={Wang, Zhiguang and Yan, Weizhong and Oates, Tim},
#   booktitle={2017 International joint conference on neural networks (IJCNN)},
#   pages={1578--1585},
#   year={2017},
#   organization={IEEE}
# }

__author__ = "James Large"

import keras
import numpy as np
import pandas as pd

from sktime_dl.classifiers.deeplearning._base import BaseDeepClassifier

class FCNClassifier(BaseDeepClassifier):

    def __init__(self,
                 random_seed=0,
                 verbose=False):
        self.verbose = verbose

        # calced in fit
        self.classes_ = None
        self.nb_classes = -1
        self.input_shape = None
        self.model = None
        self.history = None

        # predefined
        self.nb_epochs = 2000
        self.batch_size = 16
        self.callbacks = None

        self.random_seed = random_seed
        self.random_state = np.random.RandomState(self.random_seed)

    def build_model(self, input_shape, nb_classes, **kwargs):
        input_layer = keras.layers.Input(input_shape)

        conv1 = keras.layers.Conv1D(filters=128, kernel_size=8, padding='same')(input_layer)
        conv1 = keras.layers.normalization.BatchNormalization()(conv1)
        conv1 = keras.layers.Activation(activation='relu')(conv1)

        conv2 = keras.layers.Conv1D(filters=256, kernel_size=5, padding='same')(conv1)
        conv2 = keras.layers.normalization.BatchNormalization()(conv2)
        conv2 = keras.layers.Activation('relu')(conv2)

        conv3 = keras.layers.Conv1D(128, kernel_size=3, padding='same')(conv2)
        conv3 = keras.layers.normalization.BatchNormalization()(conv3)
        conv3 = keras.layers.Activation('relu')(conv3)

        gap_layer = keras.layers.pooling.GlobalAveragePooling1D()(conv3)

        output_layer = keras.layers.Dense(nb_classes, activation='softmax')(gap_layer)

        model = keras.models.Model(inputs=input_layer, outputs=output_layer)

        model.compile(loss='categorical_crossentropy', optimizer=keras.optimizers.Adam(),
                      metrics=['accuracy'])

        reduce_lr = keras.callbacks.ReduceLROnPlateau(monitor='loss', factor=0.5, patience=50,
                                                      min_lr=0.0001)

        # file_path = self.output_directory + 'best_model.hdf5'
        # model_checkpoint = keras.callbacks.ModelCheckpoint(filepath=file_path, monitor='loss',
        #                                                   save_best_only=True)
        # self.callbacks = [reduce_lr, model_checkpoint]
        self.callbacks = [reduce_lr]

        return model

    def fit(self, X, y, **kwargs):

        if isinstance(X, pd.DataFrame):
            if X.shape[1] > 1 or not isinstance(X.iloc[0, 0], pd.Series):
                raise TypeError(
                    "Input should either be a 2d numpy array, or a pandas dataframe with a single column of Series objects (CNN cannot yet handle multivariate problems")
            else:
                X = np.asarray([a.values for a in X.iloc[:, 0]])

        if len(X.shape) == 2:
            # add a dimension to make it multivariate with one dimension
            X = X.reshape((X.shape[0], X.shape[1], 1))

        y_onehot = self.convert_y(y)
        self.input_shape = X.shape[1:]

        self.batch_size = int(min(X.shape[0] / 10, self.batch_size))

        self.model = self.build_model(self.input_shape, self.nb_classes)

        if self.verbose:
            self.model.summary()

        self.history = self.model.fit(X, y_onehot, batch_size=self.batch_size, epochs=self.nb_epochs,
                                      verbose=self.verbose, callbacks=self.callbacks)

max_issues_repo_name: abdkg/sktime-dl

max_forks_repo_name: abdkg/sktime-dl


text:
import numpy as np

def cal_chi(fm, x):
    """
    Calculates ChI with given fuzzy measure and input

    :param fm: Fuzzy measure
    :param x: Input
    :return: Single value Chi output
    """
    pi_i = np.argsort(-x) + 1 # Arg sort of input, with the smallest index being 1
    ch = x[pi_i[0] - 1] * (fm[tuple(pi_i[:1])])
    for i in range(1, len(x)):
        latt_pti = np.sort(pi_i[:i+1])
        latt_ptimin1 = np.sort(pi_i[:i])
        ch = ch + x[pi_i[i] - 1] * (fm[tuple(latt_pti)] - fm[tuple(latt_ptimin1)])
    return ch

def get_cal_chi(fm):
    return lambda x: cal_chi(fm, x)

def get_cal_chi_multi(fm):
    return lambda arr, axis: np.apply_along_axis(get_cal_chi(fm), axis, arr)


text:
'''
Lucas-Kanade sparse optical flow tracker. Uses goodFeaturesToTrack
for track initialization and back-tracking for match verification
between frames.
'''
from math import sqrt
import numpy as np
import cv2

class Tracker:
    """Lucas-Kanade sparse optical flow tracker"""

    def __init__(self):
        self.track_len = 5
        self.tracks = []
        self.lk_params = dict(winSize=(15, 15),
                              maxLevel=2,
                              criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))
        self.feature_params = dict(maxCorners=500,
                                   qualityLevel=0.3,
                                   minDistance=7,
                                   blockSize=7)

    def update_tracks(self, img_old, img_new):
        """Update tracks."""
        # Get old points, using the latest one.
        points_old = np.float32([track[-1]
                                 for track in self.tracks]).reshape(-1, 1, 2)

        # Get new points from old points.
        points_new, _st, _err = cv2.calcOpticalFlowPyrLK(
            img_old, img_new, points_old, None, **self.lk_params)

        # Get inferred old points from new points.
        points_old_inferred, _st, _err = cv2.calcOpticalFlowPyrLK(
            img_new, img_old, points_new, None, **self.lk_params)

        # Compare between old points and inferred old points
        error_term = abs(
            points_old - points_old_inferred).reshape(-1, 2).max(-1)
        point_valid = error_term < 1

        new_tracks = []
        for track, (x, y), good_flag in zip(self.tracks, points_new.reshape(-1, 2), point_valid):
            # Track is good?
            if not good_flag:
                continue

            # New point is good, add to track.
            track.append((x, y))

            # Need to drop first old point?
            if len(track) > self.track_len:
                del track[0]

            # Track updated, add to track groups.
            new_tracks.append(track)

        # New track groups got, do update.
        self.tracks = new_tracks

    def get_new_tracks(self, frame, roi):
        """Get new tracks every detect_interval frames."""
        # Using mask to determine where to look for feature points.
        mask = np.zeros_like(frame)
        mask[roi[0]:roi[1], roi[2]:roi[3]] = 255

        # Get good feature points.
        feature_points = cv2.goodFeaturesToTrack(
            frame, mask=mask, **self.feature_params)

        if feature_points is not None:
            for x, y in np.float32(feature_points).reshape(-1, 2):
                self.tracks.append([(x, y)])

    def get_average_track_length(self):
        """Get the average track length"""
        length = 0
        tracks = np.array(self.tracks)
        def distance(track):
            """Get distance between the first and last point."""
            delta_x = abs(track[-1][0] - track[0][0])
            delta_y = abs(track[-1][1] - track[0][1])
            return sqrt(delta_x*delta_x + delta_y*delta_y)
        for track in tracks:
            length += distance(track)
        return length / len(tracks)

    def draw_track(self, image):
        """Draw track lines on image."""
        cv2.polylines(image, [np.int32(track)
                              for track in self.tracks], False, (0, 255, 0))

def main():
    """Test code"""
    import sys
    try:
        video_src = sys.argv[1]
    except:
        video_src = 0

    tracker = Tracker()

    cam = cv2.VideoCapture(video_src)
    detect_interval = 5
    frame_idx = 0

    prev_gray = cam.read()
    while True:
        _ret, frame = cam.read()
        frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

        # Update tracks.
        if len(tracker.tracks) > 0:
            tracker.update_tracks(prev_gray, frame_gray)

        # Get new tracks every detect_interval frames.
        target_box = [100, 400, 100, 400]
        if frame_idx % detect_interval == 0:
            tracker.get_new_tracks(frame_gray, target_box)

        # Draw tracks
        tracker.draw_track(frame)

        frame_idx += 1
        prev_gray = frame_gray
        cv2.imshow('lk_track', frame)
        ch = cv2.waitKey(1)
        if ch == 27 & ord("q"):
            break

if __name__ == '__main__':
    main()


text:
"""Run the Sample ACE problem from [Breiman85]_."""

import numpy.random
import scipy.special

from ace import ace

def build_sample_ace_problem_breiman85(N=200):
    """Sample problem from Breiman 1985."""
    x_cubed = numpy.random.standard_normal(N)
    x = scipy.special.cbrt(x_cubed)
    noise = numpy.random.standard_normal(N)
    y = numpy.exp((x ** 3.0) + noise)
    return [x], y

def build_sample_ace_problem_breiman2(N=500):
    """Build sample problem y(x) = exp(sin(x))."""
    x = numpy.linspace(0, 1, N)
    # x = numpy.random.uniform(0, 1, size=N)
    noise = numpy.random.standard_normal(N)
    y = numpy.exp(numpy.sin(2 * numpy.pi * x)) + 0.0 * noise
    return [x], y

def run_breiman85():
    """Run Breiman 85 sample."""
    x, y = build_sample_ace_problem_breiman85(200)
    ace_solver = ace.ACESolver()
    ace_solver.specify_data_set(x, y)
    ace_solver.solve()
    try:
        ace.plot_transforms(ace_solver, 'sample_ace_breiman85.png')
    except ImportError:
        pass
    return ace_solver

def run_breiman2():
    """Run Breiman's other sample problem."""
    x, y = build_sample_ace_problem_breiman2(500)
    ace_solver = ace.ACESolver()
    ace_solver.specify_data_set(x, y)
    ace_solver.solve()
    try:
        plt = ace.plot_transforms(ace_solver, None)
    except ImportError:
        pass

    plt.subplot(1, 2, 1)
    phi = numpy.sin(2.0 * numpy.pi * x[0])
    plt.plot(x[0], phi, label='analytic')
    plt.legend()
    plt.subplot(1, 2, 2)
    y = numpy.exp(phi)
    plt.plot(y, phi, label='analytic')
    plt.legend(loc='lower right')
    # plt.show()
    plt.savefig('no_noise_linear_x.png')

    return ace_solver

if __name__ == '__main__':
    run_breiman2()

max_stars_repo_name: partofthething/ace

max_issues_repo_name: Jimmy-INL/ace

max_forks_repo_name: Jimmy-INL/ace


text:
import datetime

import pandas
import statsmodels.api
import yaml
from cassandra.cqlengine import connection

from toolbox.cassandra_object_mapper_models import PlmnProcessedCord
from .utils import parse_check_date

def calculate_cluster_decomposition(start_date, end_date, kpi_basename, cord_id, acronym, frequency):
    start_date = parse_check_date(start_date)
    end_date = parse_check_date(end_date)

    # frequency = options.get('frequency')
    if not frequency:
        frequency = 31
    else:
        frequency = int(frequency)

    if not start_date and not end_date:
        return {"error": "Incorrect dates."}, 400

    with open("config.yml", 'r') as yml_file:
        config = yaml.load(yml_file)['database_options']

    connection.setup([config['address']], config['keyspace'])
    step = datetime.timedelta(days=1)
    kpi_basename = kpi_basename.upper()
    frame_setup = {'values': [], 'dates': []}
    while start_date < end_date:
        result = PlmnProcessedCord.objects.filter(cord_id=cord_id).filter(date=start_date) \
            .filter(kpi_basename=kpi_basename).filter(acronym=acronym)
        start_date += step
        for row in result:
            frame_setup['values'].append(row.value)
            frame_setup['dates'].append(row.date)

    if not len(frame_setup['values']):
        return {"error": "No data found for given parameters."}, 400

    if frequency >= len(frame_setup['values']):
        return {"error": "Frequency can't be higher than amount of found data."}, 400

    data_frame = pandas.DataFrame(frame_setup)
    data_frame = data_frame.set_index(['dates'])

    decomp = statsmodels.api.tsa.seasonal_decompose(data_frame, freq=frequency, model='additive')

    cut_nan = int(frequency / 2)

    data = {"cord_id": cord_id,
            "acronym": acronym,
            "kpi_basename": kpi_basename,
            "start_date": start_date,
            "end_date": end_date,
            "seasonal_values": decomp.seasonal['values'].tolist(),
            "seasonal_dates": decomp.seasonal.index.tolist(),
            "trend_values": decomp.trend['values'].tolist()[cut_nan:-cut_nan],
            "trend_dates": decomp.trend.index.tolist()[cut_nan:-cut_nan],
            "observed_values": decomp.observed['values'].tolist(),
            "observed_dates": decomp.observed.index.tolist()
            }
    return data, 200

reason: import statsmodels


text:
import math
import numpy as np
import itertools
import matplotlib.pyplot as plt

from .train_results import FitResult

def tensors_as_images(
    tensors, nrows=1, figsize=(8, 8), titles=[], wspace=0.1, hspace=0.2, cmap=None
):
    """
    Plots a sequence of pytorch tensors as images.

    :param tensors: A sequence of pytorch tensors, should have shape CxWxH
    """
    assert nrows > 0

    num_tensors = len(tensors)

    ncols = math.ceil(num_tensors / nrows)

    fig, axes = plt.subplots(
        nrows=nrows,
        ncols=ncols,
        figsize=figsize,
        gridspec_kw=dict(wspace=wspace, hspace=hspace),
        subplot_kw=dict(yticks=[], xticks=[]),
    )
    axes_flat = axes.reshape(-1)

    # Plot each tensor
    for i in range(num_tensors):
        ax = axes_flat[i]

        image_tensor = tensors[i]
        assert image_tensor.dim() == 3  # Make sure shape is CxWxH

        image = image_tensor.numpy()
        image = image.transpose(1, 2, 0)
        image = image.squeeze()  # remove singleton dimensions if any exist

        # Scale to range 0..1
        min, max = np.min(image), np.max(image)
        image = (image - min) / (max - min)

        ax.imshow(image, cmap=cmap)

        if len(titles) > i and titles[i] is not None:
            ax.set_title(titles[i])

    # If there are more axes than tensors, remove their frames
    for j in range(num_tensors, len(axes_flat)):
        axes_flat[j].axis("off")

    return fig, axes

def dataset_first_n(
    dataset, n, show_classes=False, class_labels=None, random_start=True, **kw
):
    """
    Plots first n images of a dataset containing tensor images.
    """

    if random_start:
        start = np.random.randint(0, len(dataset) - n)
        stop = start + n
    else:
        start = 0
        stop = n

    # [(img0, cls0), ..., # (imgN, clsN)]
    first_n = list(itertools.islice(dataset, start, stop))

    # Split (image, class) tuples
    first_n_images, first_n_classes = zip(*first_n)

    if show_classes:
        titles = first_n_classes
        if class_labels:
            titles = [class_labels[cls] for cls in first_n_classes]
    else:
        titles = []

    return tensors_as_images(first_n_images, titles=titles, **kw)

def plot_fit(fit_res: FitResult, fig=None, log_loss=False, legend=None):
    """
    Plots a FitResult object.
    Creates four plots: train loss, test loss, train acc, test acc.
    :param fit_res: The fit result to plot.
    :param fig: A figure previously returned from this function. If not None,
        plots will the added to this figure.
    :param log_loss: Whether to plot the losses in log scale.
    :param legend: What to call this FitResult in the legend.
    :return: The figure.
    """
    if fig is None:
        fig, axes = plt.subplots(
            nrows=2, ncols=2, figsize=(16, 10), sharex="col", sharey=False
        )
        axes = axes.reshape(-1)
    else:
        axes = fig.axes

    for ax in axes:
        for line in ax.lines:
            if line.get_label() == legend:
                line.remove()

    p = itertools.product(["train", "test"], ["loss", "acc"])
    for idx, (traintest, lossacc) in enumerate(p):
        ax = axes[idx]
        attr = f"{traintest}_{lossacc}"
        data = getattr(fit_res, attr)
        h = ax.plot(np.arange(1, len(data) + 1), data, label=legend)
        ax.set_title(attr)
        if lossacc == "loss":
            ax.set_xlabel("Iteration #")
            ax.set_ylabel("Loss")
            if log_loss:
                ax.set_yscale("log")
                ax.set_ylabel("Loss (log)")
        else:
            ax.set_xlabel("Epoch #")
            ax.set_ylabel("Accuracy (%)")
        if legend:
            ax.legend()
        ax.grid(True)

    return fig, axes


text:
import numpy as np
import biorbd

from biorbd_optim import (
    Instant,
    OptimalControlProgram,
    ProblemType,
    Objective,
    Constraint,
    Bounds,
    QAndQDotBounds,
    InitialConditions,
    ShowResult,
    InterpolationType,
)

def prepare_ocp(biorbd_model_path, number_shooting_points, final_time, initial_guess=InterpolationType.CONSTANT):
    # --- Options --- #
    # Model path
    biorbd_model = biorbd.Model(biorbd_model_path)
    nq = biorbd_model.nbQ()
    nqdot = biorbd_model.nbQdot()
    ntau = biorbd_model.nbGeneralizedTorque()
    torque_min, torque_max, torque_init = -100, 100, 0

    # Add objective functions
    objective_functions = {"type": Objective.Lagrange.MINIMIZE_TORQUE, "weight": 100}

    # Dynamics
    problem_type = ProblemType.torque_driven

    # Constraints
    constraints = (
        {"type": Constraint.ALIGN_MARKERS, "instant": Instant.START, "first_marker_idx": 0, "second_marker_idx": 1,},
        {"type": Constraint.ALIGN_MARKERS, "instant": Instant.END, "first_marker_idx": 0, "second_marker_idx": 2,},
    )

    # Path constraint
    X_bounds = QAndQDotBounds(biorbd_model)
    X_bounds.min[1:6, [0, -1]] = 0
    X_bounds.max[1:6, [0, -1]] = 0
    X_bounds.min[2, -1] = 1.57
    X_bounds.max[2, -1] = 1.57

    # Define control path constraint
    U_bounds = Bounds([torque_min] * ntau, [torque_max] * ntau,)

    # Initial guesses
    if initial_guess == InterpolationType.CONSTANT:
        x = [0] * (nq + nqdot)
        u = [torque_init] * ntau
    elif initial_guess == InterpolationType.CONSTANT_WITH_FIRST_AND_LAST_DIFFERENT:
        x = np.array([[1.0, 0.0, 0.0, 0, 0, 0], [1.5, 0.0, 0.785, 0, 0, 0], [2.0, 0.0, 1.57, 0, 0, 0]]).T
        u = np.array([[1.45, 9.81, 2.28], [0, 9.81, 0], [-1.45, 9.81, -2.28]]).T
    elif initial_guess == InterpolationType.LINEAR:
        x = np.array([[1.0, 0.0, 0.0, 0, 0, 0], [2.0, 0.0, 1.57, 0, 0, 0]]).T
        u = np.array([[1.45, 9.81, 2.28], [-1.45, 9.81, -2.28]]).T
    elif initial_guess == InterpolationType.EACH_FRAME:
        x = np.random.random((nq + nqdot, number_shooting_points + 1))
        u = np.random.random((ntau, number_shooting_points))
    else:
        raise RuntimeError("Initial guess not implemented yet")
    X_init = InitialConditions(x, interpolation_type=initial_guess)
    U_init = InitialConditions(u, interpolation_type=initial_guess)
    # ------------- #

    return OptimalControlProgram(
        biorbd_model,
        problem_type,
        number_shooting_points,
        final_time,
        X_init,
        U_init,
        X_bounds,
        U_bounds,
        objective_functions,
        constraints,
    )

if __name__ == "__main__":
    for initial_guess in InterpolationType:
        print(f"Solving problem using {initial_guess} initial guess")
        ocp = prepare_ocp("cube.bioMod", number_shooting_points=30, final_time=2, initial_guess=initial_guess)
        sol = ocp.solve()
        print("\n")

    # Save the last ocp
    ocp.save(sol, "cube_ocp_sol")
    ocp_loaded, sol_loaded = OptimalControlProgram.load("cube_ocp_sol.bo")

    # Rerun the loaded OCP
    sol_loaded_and_reran = ocp_loaded.solve()

    # Print the last solution
    result_plot = ShowResult(ocp, sol)
    result_plot.graphs()

max_stars_repo_name: Dangzilla/BiorbdOptim

max_issues_repo_name: Dangzilla/BiorbdOptim

max_forks_repo_name: Dangzilla/BiorbdOptim


text:
"""
utility functions for narps analysis
"""

import os
import glob
import nilearn.input_data
import numpy
import pandas
import nibabel
from scipy.stats import norm, t
import scipy.stats
from datetime import datetime
from sklearn.metrics import cohen_kappa_score

def stringify_dict(d):
    """create a pretty version of arguments for printing"""
    if 'self' in d:
        del d['self']
    s = 'Arguments:\n'
    for k in d:
        if not isinstance(d[k], str):
            d[k] = str(d[k])
        s = s + '%s: %s\n' % (k, d[k])
    return(s)

def log_to_file(fname, s, flush=False,
                add_timestamp=True,
                also_print=True,
                headspace=0):
    """ save string to log file"""
    if flush and os.path.exists(fname):
        os.remove(fname)
    if not isinstance(s, str):
        s = str(s)
    # add spacing before line
    if headspace > 0:
        s = os.linesep*headspace + s
    with open(fname, 'a+') as f:
        if also_print:
            print(s)
        f.write(s + os.linesep)
        if flush and add_timestamp:
            f.write(datetime.isoformat(
                datetime.now()) + 2 * os.linesep)

def get_masked_data(hyp, mask_img, output_dir,
                    imgtype='unthresh', dataset='zstat'):
    """
    load data from within mask
    """
    if imgtype == 'unthresh':
        hmaps = glob.glob(os.path.join(
            output_dir,
            '%s/*/hypo%d_unthresh.nii.gz' % (dataset, hyp)))
    elif imgtype == 'thresh':
        hmaps = glob.glob(os.path.join(
            output_dir,
            '%s/*/hypo%d_thresh.nii.gz' % (dataset, hyp)))
    else:
        raise Exception('bad imgtype argument')
    hmaps.sort()
    masker = nilearn.input_data.NiftiMasker(mask_img=mask_img)
    maskdata = masker.fit_transform(hmaps)  # combined_data
    maskdata = numpy.nan_to_num(maskdata)
    if imgtype == 'thresh':
        maskdata = (maskdata > 1e-4).astype('float')
    labels = [os.path.basename(os.path.dirname(i)).split('_')[1]
              for i in hmaps]
    return(maskdata, labels)

# load concatenated data - this is meant to replace
# get_masked_data()
def get_concat_data(hyp, mask_img, output_dir,
                    imgtype='unthresh', dataset='zstat',
                    vox_mask_thresh=None,
                    logfile=None):
    """
    load data from within mask
    if vox_mask_thresh is specified, then the relevant
    file is loaded and only voxels with at least
    this proportion of teams present will be included
    """

    concat_file = os.path.join(
        output_dir,
        '%s_concat_%s' % (imgtype, dataset),
        'hypo%d.nii.gz' % hyp)
    assert os.path.exists(concat_file)

    labelfile = concat_file.replace('.nii.gz', '.labels')
    assert os.path.exists(labelfile)
    labels = []
    with open(labelfile, 'r') as f:
        for l in f.readlines():
            l_s = l.strip().split()
            labels.append(l_s[0])
    masker = nilearn.input_data.NiftiMasker(mask_img=mask_img)
    maskdata = masker.fit_transform(concat_file)  # combined_data
    maskdata = numpy.nan_to_num(maskdata)
    if imgtype == 'thresh':
        maskdata = (maskdata > 1e-4).astype('float')
    if vox_mask_thresh is not None:
        assert vox_mask_thresh >= 0 and vox_mask_thresh <= 1
        mask_file = concat_file.replace(
            '.nii.gz', '_voxelmap.nii.gz'
        )
        voxmaskdata = masker.fit_transform(mask_file)
        maskdata = maskdata[:, voxmaskdata[0, :] >= vox_mask_thresh]
        if logfile is not None:
            log_to_file(
                logfile,
                'hyp %d: number of nonzero voxels (%s %s): %d' %
                (hyp, imgtype, dataset,
                 numpy.sum(voxmaskdata[0, :] >= vox_mask_thresh)))
    return(maskdata, labels)

def get_metadata(metadata_file,
                 index_var='teamID'):
    """ get team metadata"""
    metadata = pandas.read_excel(metadata_file, header=1)
    metadata.teamID = [i.strip() for i in metadata.teamID]
    metadata.index = metadata[index_var].values

    # fix issues with metadata
    metadata['used_fmriprep_data'] = [
        i.strip().split(',')[0]
        for i in metadata['used_fmriprep_data']]
    # manual fixes to textual responses
    metadata['used_fmriprep_data'] = metadata[
        'used_fmriprep_data'].replace({'Yas': 'Yes'})
    metadata.loc['E6R3', 'n_participants'] = 108
    metadata.loc['J7F9', 'n_participants'] = 107
    metadata['n_participants'] = [
        int(i.split('\n')[0]) if isinstance(i, str)
        else i for i in metadata['n_participants']]
    metadata['NV_collection_string'] = [
        os.path.basename(i.strip('/')) for i in
        metadata['NV_collection_link']]
    return(metadata)

def get_map_metadata(map_metadata_file):
    """
    get selected metadata from map metadata file
    """
    map_info = pandas.read_csv(
        map_metadata_file,
        names=['timestamp', 'teamID', 'software', 'unthresh_type',
               'thresh_type', 'MNItemplate', 'hyp5_direction',
               'hyp6_direction', 'hyp9_direction', 'comments'],
        skiprows=1)
    # manual fixes
    map_info.teamID = [i.upper() for i in map_info.teamID]
    del map_info['timestamp']
    map_info.index = map_info.teamID
    map_info = map_info.drop_duplicates(
        subset='teamID', keep='last')
    map_info.loc[:, 'unthresh_type'] = [
        i.split('values')[0].strip()
        for i in map_info.unthresh_type]
    # manual fixes
    map_info.loc['E3B6', 'unthresh_type'] = 't'
    map_info.loc['5G9K', 'unthresh_type'] = 't'
    map_info.loc['DC61', 'unthresh_type'] = 't'
    # for those that don't fit, set to NA
    map_info.loc[:, 'unthresh_type'] = [
        i if i in ['t', 'z'] else 'NA' for i in map_info.unthresh_type]

    return(map_info)

def get_decisions(decisions_file, tidy=False):
    colnames = ['teamID']
    for hyp in range(1, 10):
        colnames += ['Decision%d' % hyp,
                     'Confidence%d' % hyp,
                     'Similar%d' % hyp]
    colnames += ['collection']
    decisions = pandas.read_excel(decisions_file, skiprows=1,
                                  encoding='utf-8')
    decisions.columns = colnames

    # make a tidy version
    if tidy:
        decisions_long = pandas.melt(
            decisions,
            id_vars=['teamID'],
            value_vars=decisions.columns.values[1:28])
        decisions_long['vartype'] = [
            i[:-1] for i in decisions_long['variable']]
        decisions_long['varnum'] = [
            i[-1] for i in decisions_long['variable']]
        del decisions_long['variable']

        Decision_df = decisions_long.query('vartype =="Decision"')
        Similar_df = decisions_long.query('vartype =="Similar"')
        Confidence_df = decisions_long.query('vartype =="Confidence"')

        decision_df = Decision_df.merge(
            Similar_df,
            'left',
            on=['teamID', 'varnum'],
            suffixes=['_decision', '_similar']).merge(
                Confidence_df, 'left', on=['teamID', 'varnum'])

        del decision_df['vartype_decision']
        del decision_df['vartype_similar']
        del decision_df['vartype']

        decision_df.columns = ['teamID', 'Decision',
                               'varnum', 'Similar', 'Confidence']
        decision_df['Decision'] = (
            decision_df['Decision'] == 'Yes').astype('int')
        decision_df['Similar'] = decision_df['Similar'].astype('int')
        decision_df['Confidence'] = decision_df['Confidence'].astype('int')
        decision_df.head()
        return(decision_df)
    else:
        return(decisions)

def get_merged_metadata_decisions(metadata_file, decisions_file,):
    """ get all metadata in tidy format"""
    metadata = get_metadata(metadata_file)
    decision_df = get_decisions(decisions_file, tidy=True)
    alldata_df = decision_df.merge(metadata, on='teamID', how='left')
    return(alldata_df)

def get_teamID_to_collectionID_dict(metadata):
    """create dictionary mapping from teamID to collection ID"""
    teamid_to_collectionID_dict = {}
    for i in metadata.index:
        idx = i.split()[0]
        teamid_to_collectionID_dict[idx] = '_'.join(
            [metadata.loc[i, 'NV_collection_string'].strip(), idx])
    return(teamid_to_collectionID_dict)

def TtoZ(tmapfile, outfile, df):
    """
    takes a nibabel file object and converts from z to t
    using Hughett's transform
    adapted from:
    https://github.com/vsoch/TtoZ/blob/master/TtoZ/scripts.py
    """

    mr = nibabel.load(tmapfile)
    data = mr.get_data()

    # Select just the nonzero voxels
    nonzero = data[data != 0]

    # We will store our results here
    Z = numpy.zeros(len(nonzero))

    # Select values less than or == 0, and greater than zero
    c = numpy.zeros(len(nonzero))
    k1 = (nonzero <= c)
    k2 = (nonzero > c)

    # Subset the data into two sets
    t1 = nonzero[k1]
    t2 = nonzero[k2]

    # Calculate p values for <=0
    p_values_t1 = t.cdf(t1, df=df)
    z_values_t1 = norm.ppf(p_values_t1)

    # Calculate p values for > 0
    p_values_t2 = t.cdf(-t2, df=df)
    z_values_t2 = -norm.ppf(p_values_t2)
    Z[k1] = z_values_t1
    Z[k2] = z_values_t2

    # Write new image to file
    empty_nii = numpy.zeros(mr.shape)
    empty_nii[mr.get_data() != 0] = Z
    Z_nii_fixed = nibabel.nifti1.Nifti1Image(
        empty_nii,
        affine=mr.get_affine(),
        header=mr.get_header())
    nibabel.save(Z_nii_fixed, outfile)

def t_corr(y, res_mean=None, res_var=None, Q=None):
    """
    perform a one-sample t-test on correlated data
    y = data (n observations X n vars)
    res_mean = Common mean over voxels and results
    res_var  = Common variance over voxels and results
    Q = "known" correlation across observations
    - (use empirical correlation based on maps)
    """

    npts = y.shape[0]
    X = numpy.ones((npts, 1))

    if res_mean is None:
        res_mean = 0

    if res_var is None:
        res_var = 1

    if Q is None:
        Q = numpy.eye(npts)

    VarMean = res_var * X.T.dot(Q).dot(X) / npts**2

    # T  =  mean(y,0)/s-hat-2
    # use diag to get s_hat2 for each variable
    T = (numpy.mean(y, 0)-res_mean
         )/numpy.sqrt(VarMean)*numpy.sqrt(res_var) + res_mean

    # Assuming variance is estimated on whole image
    # and assuming infinite df
    p = 1 - scipy.stats.norm.cdf(T)

    return(T, p)

def randn_from_shape(shape):
    """
    take in a tuple defining a 4d matrix shape,
    and return a random matrix of that shape
    """
    assert len(shape) == 4
    return(
        numpy.random.randn(
            shape[0],
            shape[1],
            shape[2],
            shape[3]))

def matrix_kappa_score(d):
    """
    compute cohen's kappa for each combination
    of rows
    """
    score = numpy.eye(d.shape[0])
    for i in range(d.shape[0]):
        for j in range(i + 1, d.shape[0]):
            score[i, j] = cohen_kappa_score(d[i, :], d[j, :])
            score[j, i] = score[i, j]
    return(score)

def matrix_pct_agreement(d):
    """
    compute cohen's kappa for each combination
    of rows
    """
    score = numpy.eye(d.shape[0])
    for i in range(d.shape[0]):
        for j in range(i + 1, d.shape[0]):
            score[i, j] = numpy.mean(d[i, :] == d[j, :])
            score[j, i] = score[i, j]
    return(score)

max_stars_repo_name: nicholst/narps

max_issues_repo_name: nicholst/narps

max_forks_repo_name: nicholst/narps

reason: import numpy,import scipy,from scipy


text:
# Polynomial Regression

# Importing the libraries
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Importing the dataset
dataset = pd.read_csv('Position_Salaries.csv')
X = dataset.iloc[:, 1:2].values  #here if we put [:,1]-then its treated as a vector and not as a matrix so we put [:,1:2]- so 10 rows and 1 column
y = dataset.iloc[:, 2].values  #here salary,X is the level
#There is a non-linear relation between the level and salary so we use polynomial regression for this data
#Here there is only one dependent column and one independent column-level is independent and salary is depenedent,based on the level we try to predict the salary

# Splitting the dataset into the Training set and Test set
"""from sklearn.cross_validation import train_test_split #we dont have much observations to split our data into training set and split test so thats why we dont split the dataset
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 0)""" #it wont make much sense if we split the dataset

# Feature Scaling
"""from sklearn.preprocessing import StandardScaler
sc_X = StandardScaler()
X_train = sc_X.fit_transform(X_train)
X_test = sc_X.transform(X_test)"""  #we use the same library as the simple,multiple linear regression so no need to do any feature scaling here

# Fitting Linear Regression to the dataset
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)  # fit the regressor to the dataset-so its X,y and not X_train and y_train

# Fitting Polynomial Regression to the dataset
# Its like we are fitting a regression line for the 4th degree polynomial-this is nothing but polynomial regression
from sklearn.preprocessing import PolynomialFeatures  # everytime its the same import some class(PolynomialFeatures) from some library and then create an object of that class
poly_reg = PolynomialFeatures(degree = 4)  # First fit the 4th degree model to X-degree 4 so we will have sa x,x^2,x^3 till x^4(poly terms of the independent variable)
X_poly = poly_reg.fit_transform(X)  # fit and transform-so add the additional polynomial terms
poly_reg.fit(X_poly, y)  # poly_reg object automatically includes a column of 1 into the matrix(50 ones to the start)
lin_reg_2 = LinearRegression()
lin_reg_2.fit(X_poly, y)  # we build 2 models here-linear regression and polynomial regression just to compare both the models.

# Visualising the Linear Regression results
plt.scatter(X, y, color = 'red')
plt.plot(X, lin_reg.predict(X), color = 'blue')  #x coordinates so X,y coordinates so the predicted salaries by the LinReg model
plt.title('Truth or Bluff (Linear Regression)')
plt.xlabel('Position level')
plt.ylabel('Salary')
plt.show()  #So in red we have the obs points and the predicted points in blue(the straight line)

# Visualising the Polynomial Regression results
plt.scatter(X, y, color = 'red')
plt.plot(X, lin_reg_2.predict(poly_reg.fit_transform(X)), color = 'blue') #we just cant put lin_reg_2 here because its of the class LinReg
plt.title('Truth or Bluff (Polynomial Regression)') #here again the predicted values of the polynomial regression so poly_reg.fit_transform(X)
plt.xlabel('Position level')
plt.ylabel('Salary')
plt.show()  #we get to know that the polynomial model fits much better than the linear regression model for the dataset

# Visualising the Polynomial Regression results (for higher resolution and smoother curve)
X_grid = np.arange(min(X), max(X), 0.1)
X_grid = X_grid.reshape((len(X_grid), 1))  #reshape the vector into the matrix,so 1 column and 90 levels,all the levels incremented by 0.1
plt.scatter(X, y, color = 'red')
plt.plot(X_grid, lin_reg_2.predict(poly_reg.fit_transform(X_grid)), color = 'blue') # This is nothing but the linear regression model but fitting with for a higher power function(which is nothing but Polynomial Regression)
plt.title('Truth or Bluff (Polynomial Regression)')
plt.xlabel('Position level')
plt.ylabel('Salary')
plt.show()

# Predicting a new result with Linear Regression
lin_reg.predict(6.5)  #gives the corresponding salary for the level 6.5 in the linear model

# Predicting a new result with Polynomial Regression
lin_reg_2.predict(poly_reg.fit_transform(6.5)) #gives the corresponding salary for the level 6.5 in the polynomial model
# This gives more accurate results than the linear regression model-so this model is preferred here
# After seeing the visualization we can see that our curve has actually overfitted the model but it seems fine because it has given a decent prediction


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# $Id: test_maximal_independent_set.py 577 2011-03-01 06:07:53Z lleeoo $
#    Copyright (C) 2004-2018 by
#    Leo Lopes <leo.lopes@monash.edu>
#    Aric Hagberg <hagberg@lanl.gov>
#    Dan Schult <dschult@colgate.edu>
#    Pieter Swart <swart@lanl.gov>
#    All rights reserved.
#    BSD license.
#
# Author: Leo Lopes <leo.lopes@monash.edu>
"""
Tests for maximal (not maximum) independent sets.

"""

from nose.tools import *
import networkx as nx
import random

class TestMaximalIndependantSet(object):
    def setup(self):
        self.florentine = nx.Graph()
        self.florentine.add_edge('Acciaiuoli', 'Medici')
        self.florentine.add_edge('Castellani', 'Peruzzi')
        self.florentine.add_edge('Castellani', 'Strozzi')
        self.florentine.add_edge('Castellani', 'Barbadori')
        self.florentine.add_edge('Medici', 'Barbadori')
        self.florentine.add_edge('Medici', 'Ridolfi')
        self.florentine.add_edge('Medici', 'Tornabuoni')
        self.florentine.add_edge('Medici', 'Albizzi')
        self.florentine.add_edge('Medici', 'Salviati')
        self.florentine.add_edge('Salviati', 'Pazzi')
        self.florentine.add_edge('Peruzzi', 'Strozzi')
        self.florentine.add_edge('Peruzzi', 'Bischeri')
        self.florentine.add_edge('Strozzi', 'Ridolfi')
        self.florentine.add_edge('Strozzi', 'Bischeri')
        self.florentine.add_edge('Ridolfi', 'Tornabuoni')
        self.florentine.add_edge('Tornabuoni', 'Guadagni')
        self.florentine.add_edge('Albizzi', 'Ginori')
        self.florentine.add_edge('Albizzi', 'Guadagni')
        self.florentine.add_edge('Bischeri', 'Guadagni')
        self.florentine.add_edge('Guadagni', 'Lamberteschi')

    def test_random_seed(self):
        G = nx.complete_graph(5)
        for node in G:
            assert_equal(nx.maximal_independent_set(G, [node], seed=1), [node])

    def test_K5(self):
        """Maximal independent set: K5"""
        G = nx.complete_graph(5)
        for node in G:
            assert_equal(nx.maximal_independent_set(G, [node]), [node])

    def test_K55(self):
        """Maximal independent set: K55"""
        G = nx.complete_graph(55)
        for node in G:
            assert_equal(nx.maximal_independent_set(G, [node]), [node])

    def test_exception(self):
        """Bad input should raise exception."""
        G = self.florentine
        assert_raises(nx.NetworkXUnfeasible,
                      nx.maximal_independent_set, G, ["Smith"])
        assert_raises(nx.NetworkXUnfeasible,
                      nx.maximal_independent_set, G, ["Salviati", "Pazzi"])

    def test_digraph_exception(self):
        G = nx.DiGraph([(1, 2), (3, 4)])
        assert_raises(nx.NetworkXNotImplemented, nx.maximal_independent_set, G)

    def test_florentine_family(self):
        G = self.florentine
        indep = nx.maximal_independent_set(G, ["Medici", "Bischeri"])
        assert_equal(sorted(indep),
                     sorted(["Medici", "Bischeri", "Castellani", "Pazzi",
                             "Ginori", "Lamberteschi"]))

    def test_bipartite(self):
        G = nx.complete_bipartite_graph(12, 34)
        indep = nx.maximal_independent_set(G, [4, 5, 9, 10])
        assert_equal(sorted(indep), list(range(12)))

    def test_random_graphs(self):
        """Generate 50 random graphs of different types and sizes and
        make sure that all sets are independent and maximal."""
        for i in range(0, 50, 10):
            G = nx.random_graphs.erdos_renyi_graph(i * 10 + 1, random.random())
            IS = nx.maximal_independent_set(G)
            assert_false(list(G.subgraph(IS).edges()))
            neighbors_of_MIS = set.union(*(set(G.neighbors(v)) for v in IS))
            for v in set(G.nodes()).difference(IS):
                assert_true(v in neighbors_of_MIS)

max_stars_repo_name: rakschahsa/networkx

max_issues_repo_name: rakschahsa/networkx

max_forks_repo_name: rakschahsa/networkx

reason: import networkx


text:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
import cloudpickle

kinematic = cloudpickle.load(open('./swing_kinematic.dll', 'rb'))

def swing_animation(sol, param0):
    fig, ax = plt.subplots(figsize=(6, 6))
    plt.axis('equal')
    plt.axis([-3, 3, -3, 3])

    line1, = ax.plot([], [], lw=1, color='k', linestyle='-', marker='o', ms=5)
    line2, = ax.plot([], [], lw=2, color='b', linestyle='-', marker='o', ms=3)

    def init():
        line1.set_data([], [])
        line2.set_data([], [])
        return (line1, line2)

    def animate(i):
        y = sol[i, :]
        p1, p2, p3, p21, p31 = kinematic(sol[i, :], param0)
        line1.set_data([0, p1[0]], [0, p1[1]])
        line2.set_data([p21[0], p2[0], p3[0], p31[0]], [p21[1], p2[1], p3[1], p31[1]])
        return (line1,line2)
    # call the animator. blit=True means only re-draw the parts that have changed.
    anim = animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=1000, interval=10, blit=True)
    return anim

max_stars_repo_name: HuiminHe/PyDy

max_issues_repo_name: HuiminHe/PyDy

max_forks_repo_name: HuiminHe/PyDy