text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Python version: 3.8

import os
import copy
import json
import datetime, time
import numpy as np
import operator
from tensorboardX import SummaryWriter
from statsmodels.stats.weightstats import DescrStatsW
from tqdm import tqdm
from sklearn.datasets import make_blobs

import config
from options import args_parser
from update import LocalUpdate, test_inference
from models import MLP, CNNMnist, CNNFashion_Mnist, CNNCifar, LR
from utils import get_dataset_cluster, average_weights, mixing_weights, exp_details_cluster, get_device, init_seed,\
    complete, star, ring, independent, lr_decay, \
    plot_label_dis, plot_user_dis, optimal_mixing_cycles, update_topo

from simulator.wsn.network import Network
from simulator.wsn.utils import *
from simulator.wsn.fcm import *
from statsmodels.stats.weightstats import DescrStatsW

if __name__ == '__main__':
    start_time = time.time()

    args = args_parser()
    print(args)
    init_seed(args.seed)
    exp_details_cluster(args)
    device = get_device(args)
    fake = args.fake
    budget = config.energy_budget

    # Define paths
    tag = 'Dataset[{}]_Cluster[{}]_Sim[{:.2f}]_Topo[{}]_Mix[{}]_lr[{:.3f}]_frac[{:.1f}].pkl'. \
        format(args.dataset, args.num_clusters, args.cluster_similarity, args.topo, args.mix_ep, args.lr, args.frac)
    cur_datetime = datetime.datetime.now().strftime('%Y-%m-%d_%H:%M:%S')
    file_name = '%s_[ENERGY(%d)]_%s' % (cur_datetime, budget, tag) if not fake else \
        '%s_[ENERGY(%d)]_%s.fake' % (cur_datetime, budget, tag)
    logger_dir = '../logs/%s/' % file_name
    logger = SummaryWriter(logger_dir)

    if fake:
        file_list = os.listdir('../logs')
        file_list = sorted(file_list, reverse=True)  # recent first, file name already includes time information
        is_find = False
        for name in file_list:
            if tag in name and 'fake' not in name:  # fake result can not be used for fake training
                fake_path = os.path.join('../logs', name)
                is_find = True
                break
        if is_find:
            print("fake path:", fake_path)
            with open(os.path.join(fake_path, 'convergence.json'), 'r') as f:
                fake_rs = json.load(f)

    # Build topo
    num_clusters = args.num_clusters
    topo_matrix = None
    if args.topo == 'complete':
        topo_matrix = complete(num_clusters)
    elif args.topo == 'star':
        topo_matrix = star(num_clusters)
    elif args.topo == 'ring':
        topo_matrix = ring(num_clusters, centroids=config.centroids)
    elif args.topo == 'independent':
        topo_matrix = independent(num_clusters)
    else:
        exit('Error: unrecognized topology')
    spectral_norm = np.linalg.norm(topo_matrix - 1 / num_clusters, ord=2) ** 2
    p = 1 - spectral_norm
    print("\nTopo (p=%.5f):" % p)
    for row in topo_matrix:
        print(','.join(map(lambda x: "%.3f" % x, row)))

    # Setup wsn networks and routing protocol
    traces = {}
    init_nodes, y = make_blobs(n_samples=args.num_users, centers=config.centroids, n_features=2,
                               random_state=args.seed, cluster_std=15)
    network = Network(init_nodes=init_nodes)
    network.set_topo(topo_matrix)
    routing_protocol_class = eval(config.routing_topology)
    network.init_routing_protocol(routing_protocol_class())
    users_groups = network.get_cluster_members()
    print("Clusters: ", users_groups)
    plot_clusters(network, logger_dir)

    # Load dataset and split users into clusters.
    # If config.users_groups is None, the number of users in each cluster is equal.
    # Otherwise, the size is the same with config.users_groups.
    train_dataset, test_dataset, split = get_dataset_cluster(args, users_groups=users_groups)
    user_groups_list = split['train'][0]
    test_idxs_list = split['test'][1]

    num_classes = args.num_classes
    print("\nLabel distribution in clusters:")
    for c, user_groups in enumerate(user_groups_list):
        targets = np.array(train_dataset.targets)[np.array(list(user_groups.values())).flatten()].tolist()
        dis = [(i, targets.count(i)) for i in range(num_classes)]
        my_users = [user for user in user_groups]
        print("-Cluster %d: " % c, dis)
        # print("%d users:" % len(my_users), my_users)
        fig1 = plot_label_dis([i for i, _ in dis], [c for _, c in dis])
        logger.add_figure('Cluster/%d/Label/distribution' % c, figure=fig1)  # label distribution in a cluster
    fig2 = plot_user_dis(train_dataset, user_groups_list, num_classes, args.num_users, save_path=logger_dir)  # all user distribution
    logger.add_figure('All/User/distribution', figure=fig2)

    # Build model
    if args.model == 'cnn':
        # Convolutional neural netork
        if args.dataset == 'mnist':
            global_model = CNNMnist(args=args)
        elif args.dataset == 'fmnist':
            global_model = CNNFashion_Mnist(args=args)
        elif args.dataset == 'cifar':
            global_model = CNNCifar(args=args)
    elif args.model == 'mlp':
        # Multi-layer preceptron
        img_size = train_dataset[0][0].shape
        len_in = 1
        for x in img_size:
            len_in *= x
        print(len_in)
        global_model = MLP(dim_in=len_in, dim_hidden=64,
                           dim_out=args.num_classes)
    elif args.model == 'lr':
        # Logistic regression
        img_size = train_dataset[0][0].shape
        len_in = 1
        for x in img_size:
            len_in *= x
        global_model = LR(dim_in=len_in, dim_out=args.num_classes)
    else:
        exit('Error: unrecognized model')

    # Set the model to train and send it to device.
    global_model.to(device)
    global_model.train()
    print(global_model)
    total_params = sum([w.numel() for w in global_model.parameters()])
    print("num of params: ", total_params)

    # Copy weights to each cluster.
    global_weights = global_model.state_dict()
    cluster_wegihts_list = [copy.deepcopy(global_weights) for i in range(num_clusters)]

    # Observations in each cluster.
    cluster_params_list = [None for _ in range(num_clusters)]
    cluster_gradients_list = [None for _ in range(num_clusters)]
    cluster_data_ratio_list = [None for _ in range(num_clusters)]
    cluster_agg_delta_list = [0 for i in range(num_clusters)]
    cluster_mix_delta_list = [0 for i in range(num_clusters)]

    # Training and testing results
    train_loss_list, train_accuracy_list = [[] for _ in range(num_clusters)], [[] for _ in range(num_clusters)]
    test_loss_list, test_accuracy_list = [[] for _ in range(num_clusters)], [[] for _ in range(num_clusters)]

    global_params_list = []
    opt_tau, mixing_timer, num_mixing = 1, 1, 0
    intra_energy_list, inter_energy_list = [], []
    opt_tau_list = []
    total_energy = network.get_remaining_energy()

    print_every = 2

    for epoch in tqdm(range(args.epochs)):
        print(f'\n | Global Training Round : {epoch+1} |\n')

        print("Remain energy:", network.get_remaining_energy())
        logger.add_scalar('Observation/energy_consumption', total_energy - network.get_remaining_energy(), global_step=epoch)
        if network.get_remaining_energy() <= total_energy - budget:
            print("\n Terminated for using up energy budget (%d J)!" % budget)
            break

        if fake and args.mix_ep > 0:
            assert operator.eq(users_groups, config.users_groups)
            train_loss = fake_rs['train_loss'][epoch]
            test_loss = fake_rs['test_loss'][epoch]
            test_accuracy = fake_rs['test_acc'][epoch]
            logger.add_scalar('All/Test/acc', test_accuracy, global_step=epoch)
            logger.add_scalar('All/Test/loss', test_loss, global_step=epoch)
            logger.add_scalar('All/Train/loss', train_loss, global_step=epoch)
            final_epoch = epoch

            if args.mix_ep > 0 and (epoch + 1) % args.mix_ep == 0:
                network.activate_mix()
            else:
                network.deactivate_mix()
            traces[config.routing_topology] = network.simulate_one_round()

            # record energy
            intra_energy_list.append(network.energy_dis['intra-comm'] + network.energy_dis['local-update'])
            inter_energy_list.append(network.energy_dis['inter-comm'])
            inter_intra_ratio = np.mean(inter_energy_list) / np.mean(intra_energy_list)
            print("******", np.mean(inter_energy_list), np.mean(intra_energy_list))
            logger.add_scalar('Observation/intra', intra_energy_list[-1], global_step=epoch)
            logger.add_scalar('Observation/inter', inter_energy_list[-1], global_step=epoch)
            logger.add_scalar('Observation/inter_intra_ratio', inter_intra_ratio, global_step=epoch)

            continue

        lr = lr_decay(args.lr, epoch)
        logger.add_scalar('Observation/lr', lr, global_step=epoch)
        K = len(train_dataset) / args.num_users / args.local_bs  # 100 is the number of users

        # Used in dynamic mixing when decreasing to 0
        mixing_timer -= 1

        # Aggregation in each cluster
        total_num_participation = 0
        for c, user_groups in enumerate(user_groups_list):
            local_weights, local_losses = [], []
            global_model.load_state_dict(cluster_wegihts_list[c])

            # test
            test_acc, test_loss = test_inference(args, global_model, test_dataset, test_idxs_list[c])
            test_accuracy_list[c].append(test_acc)
            test_loss_list[c].append(test_loss)

            # train
            global_model.train()

            # get alive users
            # note that the clusters sequence number in network is not constant
            # we modified FCM to stabilize the clustering results.
            my_users_total = list(user_groups.keys())
            my_users = [node.id for node in network.get_nodes_by_membership(c)]
            total, alive = len(my_users_total), len(my_users)
            print("-" * 40)
            if alive > 0:
                print("Cluster %d is running..." % c)
                print("alive: %d, total: %d, alive-ratio: %.2f" % (alive, total, alive/total))
            else:
                print("Cluster %d is dead !" % c)

            # select participating users, which should be alive
            m = max(int(args.frac * total), 1)
            m = min(m, alive)
            idxs_users = np.random.choice(my_users, m, replace=False)

            # train models and estimate params
            params_groups = []
            init_grad_groups = []
            local_weight_delta = 0

            for idx in idxs_users:
                # initialize local model in each participating user
                local_model = LocalUpdate(args=args, dataset=train_dataset, idxs=user_groups[idx], logger=logger)

                # local updates and estimate params
                rs_params = {}
                w, loss = local_model.update_weights(model=copy.deepcopy(global_model), global_round=epoch,
                                                     rt_params=rs_params, num_params=total_params)
                local_weights.append(copy.deepcopy(w))
                local_losses.append(copy.deepcopy(loss))

                params_groups.append(rs_params['estmt_params'])

                # # approximate
                # init_grad_vector = []
                # for name in cluster_wegihts_list[c]:
                #     approxmimate_grad = (cluster_wegihts_list[c][name] - w[name]) / lr / K
                #     init_grad_vector.extend(approxmimate_grad.cpu().numpy().reshape(-1).tolist())
                # init_grad_groups.append(init_grad_vector)

                init_grad_groups.append(rs_params['init_grad_vector'])

                local_weight_delta += rs_params['weight_delta']
            local_weight_delta /= len(idxs_users)

            # aggregate models and params in a cluster
            agg_weights = average_weights(local_weights)
            cluster_gradients_list[c] = np.mean(init_grad_groups, axis=0)
            cluster_data_ratio_list[c] = len(idxs_users)
            partial_params = np.mean(params_groups, axis=0)

            gamma = np.var(init_grad_groups, axis=0).sum() / total_params
            G = np.linalg.norm(cluster_gradients_list[c]) ** 2 / total_params

            agg_weight_delta, local_grad_delta, index = 0, 0, 0
            for name in agg_weights:
                agg_weight_delta += (agg_weights[name] - cluster_wegihts_list[c][name]).norm() ** 2
                local_grad_delta += np.linalg.norm(
                    (cluster_wegihts_list[c][name] - agg_weights[name]).cpu().numpy().reshape(-1) / lr
                    - K * cluster_gradients_list[c][index:index + agg_weights[name].numel()]) ** 2
                index += agg_weights[name].numel()
            L = np.sqrt(local_grad_delta / local_weight_delta / K)

            cluster_params_list[c] = np.concatenate([partial_params, [gamma, G, L]])
            cluster_agg_delta_list[c] = agg_weight_delta.item()
            cluster_wegihts_list[c] = agg_weights
            loss_avg = sum(local_losses) / len(local_losses)
            train_loss_list[c].append(loss_avg)
            logger.add_scalar('Cluster/%d/Train/loss' % c, loss_avg, global_step=epoch)
            total_num_participation += len(idxs_users)

        logger.add_scalar('Observation/agg', np.mean(cluster_agg_delta_list), global_step=epoch)

        # Collect params from all clusters and obtain global params
        agg_params = DescrStatsW(cluster_params_list, weights=cluster_data_ratio_list)
        agg_gradients = DescrStatsW(cluster_gradients_list, weights=cluster_data_ratio_list)

        psi = agg_gradients.var.sum() / total_params  # TODO
        epsilon = np.linalg.norm(agg_gradients.mean) ** 2

        global_params = np.concatenate([agg_params.mean, [psi, epsilon]])
        global_params_list.append(global_params)

        print("Global estimated params:")
        for i, name in enumerate(['loss', 'sigma', 'gamma', 'G', 'L', 'psi', 'epsilon']):
            logger.add_scalar('Estimation/%s' % name, global_params[i], global_step=epoch)

        # Mix between clusters
        cluster_wegihts_list_mixing = []
        if args.mix_ep > 0 and (epoch + 1) % args.mix_ep == 0 or \
                args.mix_ep == 0 and mixing_timer <= 0:

            network.activate_mix()
            num_mixing += 1

            for i in range(num_clusters):
                cluster_wegihts_list_mixing.append(mixing_weights(cluster_wegihts_list, topo_matrix[:, i]))
                delta = 0
                for name in cluster_wegihts_list_mixing[i]:
                    delta += (cluster_wegihts_list_mixing[i][name] - cluster_wegihts_list[i][name]).norm() ** 2
                cluster_mix_delta_list[i] = delta.item()
            cluster_wegihts_list = cluster_wegihts_list_mixing
        else:
            network.deactivate_mix()
            cluster_mix_delta_list = [0 for i in range(num_clusters)]
        logger.add_scalar('Observation/mix', np.mean(cluster_mix_delta_list), global_step=epoch)

        # simulate one round in wsn
        traces[config.routing_topology] = network.simulate_one_round()
        if network.is_mix:
            logger.add_scalar('Observation/mix', np.mean(cluster_mix_delta_list), global_step=epoch)

        # record energy
        intra_energy_list.append(network.energy_dis['intra-comm'] + network.energy_dis['local-update'])
        inter_energy_list.append(network.energy_dis['inter-comm'])
        inter_intra_ratio = np.mean(inter_energy_list) / np.mean(intra_energy_list)
        logger.add_scalar('Observation/intra', intra_energy_list[-1], global_step=epoch)
        logger.add_scalar('Observation/inter', inter_energy_list[-1], global_step=epoch)
        logger.add_scalar('Observation/inter_intra_ratio', inter_intra_ratio, global_step=epoch)

        # mixing_timer = min(mixing_timer, opt_tau)
        if args.mix_ep == 0 and mixing_timer <= 0:

            intra_energy_per_round = intra_energy_list[-1] / len(intra_energy_list)
            inter_energy_per_round = inter_energy_list[-1] / num_mixing

            # calculate optimal mixing cycles
            loss, sigma, gamma, G, L, psi, epsilon = np.mean(global_params_list, axis=0)
            init_loss = global_params_list[0][0]
            opt_tau = optimal_mixing_cycles(L, G, sigma, gamma, psi, lr, K, p, init_loss, inter_energy_per_round, budget)
            mixing_timer = opt_tau
            opt_tau_list.append(opt_tau)
            logger.add_scalar('Observation/tau', opt_tau, global_step=epoch)

        # update mixing matrix
        alive_clusters = network.get_alive_clusters()
        topo_matrix = update_topo(topo_matrix, alive_clusters)
        spectral_norm = np.linalg.norm(topo_matrix - 1 / num_clusters, ord=2) ** 2
        p = 1 - spectral_norm
        logger.add_scalar('Estimation/p', p, global_step=epoch)

        # Print global training loss after every 'i' rounds
        train_loss = DescrStatsW(np.array(train_loss_list)[:, -1], weights=cluster_data_ratio_list).mean
        test_loss = DescrStatsW(np.array(test_loss_list)[:, -1], weights=cluster_data_ratio_list).mean
        test_accuracy = DescrStatsW(np.array(test_accuracy_list)[:, -1], weights=cluster_data_ratio_list).mean

        if (epoch+1) % print_every == 0:
            print(f' \nAvg Training Stats after {epoch+1} global rounds:')
            print(f'Training Loss : {train_loss}')
            print(f' \nAvg Testing Stats after {epoch+1} global rounds:')
            print(f'Testing Accuracy: {100*test_accuracy}%')
            print(f'Testing Loss: {test_loss}')

        logger.add_scalar('All/Test/acc', test_accuracy, global_step=epoch)
        logger.add_scalar('All/Test/loss', test_loss, global_step=epoch)
        logger.add_scalar('All/Train/loss', train_loss, global_step=epoch)

        # Print global testing result after every 'i' rounds
        # test_accuracy, test_loss = test_inference(args, global_model, test_dataset)
        # test_accuracy_list.append(test_accuracy)
        # test_loss_list.append(test_loss)

        # if (epoch+1) % print_every == 0:
        #     print(f' \nAvg Testing Stats after {epoch+1} global rounds:')
        #     print('Testing Accuracy: {:.2f}% \n'.format(100*test_accuracy))
        #     logger.add_scalar('All/Test/acc', test_accuracy, global_step=epoch)

    if not fake:
        with open(os.path.join(logger_dir, 'convergence.json'), 'w') as f:
            rs = {
                'train_loss': DescrStatsW(train_loss_list, cluster_data_ratio_list).mean.tolist(),
                'test_loss': DescrStatsW(test_loss_list, cluster_data_ratio_list).mean.tolist(),
                'test_acc': DescrStatsW(test_accuracy_list, cluster_data_ratio_list).mean.tolist()
            }
            f.write(json.dumps(rs, indent=4))
            # json.dump(dict, f, indent=4)
    else:
        with open(os.path.join(logger_dir, 'convergence.json'), 'w') as f:
            rs = {
                'train_loss': fake_rs['train_loss'][:final_epoch],
                'test_loss': fake_rs['test_loss'][:final_epoch],
                'test_acc': fake_rs['test_acc'][:final_epoch]
            }
            f.write(json.dumps(rs, indent=4))

    with open(os.path.join(logger_dir, 'user_groups_list.json'), 'w') as f:
        f.write(json.dumps(network.get_cluster_members()))

    if args.mix_ep == 0:
        with open(os.path.join(logger_dir, 'opt_tau_list.json'), 'w') as f:
            f.write(json.dumps(opt_tau_list))

    # plot cluster information
    # plot_clusters(network, logger_dir)
    # plot_time_of_death(network, logger_dir)
    # plot_traces(traces, logger_dir)

    print('\n Total Run Time: {0:0.4f}'.format(time.time()-start_time))

reason: import numpy,from statsmodels


text:
"""Description of classes Section, SectionL."""
from typing import NoReturn
import numpy

from cryspy.A_functions_base.function_2_crystallography_base import \
    calc_atoms_in_unit_cell

from cryspy.B_parent_classes.cl_1_item import ItemN
from cryspy.B_parent_classes.cl_2_loop import LoopN

from cryspy.C_item_loop_classes.cl_1_cell import Cell
from cryspy.C_item_loop_classes.cl_1_atom_site import AtomSiteL
from cryspy.C_item_loop_classes.cl_1_space_group_symop import \
    SpaceGroupSymop, SpaceGroupSymopL

class Section(ItemN):
    """Section description.

    Describe information concerning the density point.
    """

    ATTR_MANDATORY_NAMES = ("id", "size_x", "size_y", "atom_center",
                            "atom_along_axis_x", "atom_x_operation_xyz",
                            "atom_along_axis_y", "atom_y_operation_xyz",
                            "points_x", "points_y", "url_out")
    ATTR_MANDATORY_TYPES = (str, float, float, str, str, str, str, str, int,
                            int, str)
    ATTR_MANDATORY_CIF = ("id", "size_x", "size_y", "atom_center",
                          "atom_along_axis_x", "atom_x_operation_xyz",
                          "atom_along_axis_y", "atom_y_operation_xyz",
                          "points_x", "points_y", "url_out")

    ATTR_OPTIONAL_NAMES = ("field_x", "field_y", "field_z")
    ATTR_OPTIONAL_TYPES = (float, float, float)
    ATTR_OPTIONAL_CIF = ("field_x", "field_y", "field_z")

    ATTR_NAMES = ATTR_MANDATORY_NAMES + ATTR_OPTIONAL_NAMES
    ATTR_TYPES = ATTR_MANDATORY_TYPES + ATTR_OPTIONAL_TYPES
    ATTR_CIF = ATTR_MANDATORY_CIF + ATTR_OPTIONAL_CIF

    ATTR_INT_NAMES = ()
    ATTR_INT_PROTECTED_NAMES = ()

    # parameters considered are refined parameters
    ATTR_REF = ()
    ATTR_SIGMA = tuple([f"{_h:}_sigma" for _h in ATTR_REF])
    ATTR_CONSTR_FLAG = tuple([f"{_h:}_constraint" for _h in ATTR_REF])
    ATTR_REF_FLAG = tuple([f"{_h:}_refinement" for _h in ATTR_REF])

    # formats if cif format
    D_FORMATS = {"size_x": "{:.2f}", "size_y": "{:.2f}"}

    # constraints on the parameters
    D_CONSTRAINTS = {}

    # default values for the parameters
    D_DEFAULT = {}
    for key in ATTR_SIGMA:
        D_DEFAULT[key] = 0.
    for key in (ATTR_CONSTR_FLAG + ATTR_REF_FLAG):
        D_DEFAULT[key] = False

    PREFIX = "section"

    def __init__(self, **kwargs) -> NoReturn:
        super(Section, self).__init__()

        # defined for any integer and float parameters
        D_MIN = {}

        # defined for ani integer and float parameters
        D_MAX = {}

        self.__dict__["D_MIN"] = D_MIN
        self.__dict__["D_MAX"] = D_MAX
        for key, attr in self.D_DEFAULT.items():
            setattr(self, key, attr)
        for key, attr in kwargs.items():
            setattr(self, key, attr)

    def calc_axes_x_y_z(self, cell: Cell, atom_site: AtomSiteL):
        """
        Calculate three vectors of axes: pos_x, pos_y, pos_z.

        Arguments
        ---------
            - cell
            - atom_site (loop)
        """
        atom_o_label = self.atom_center
        _item_a_s_o = atom_site[atom_o_label]
        fract_a_o_xyz = numpy.array([_item_a_s_o.fract_x, _item_a_s_o.fract_y,
                                     _item_a_s_o.fract_z], dtype=float)
        pos_a_o_xyz = cell.calc_position_by_coordinate(
            fract_a_o_xyz[0], fract_a_o_xyz[1], fract_a_o_xyz[2])
        pos_a_o_xyz = numpy.array(pos_a_o_xyz, dtype=float)

        atom_x_label = self.atom_along_axis_x
        atom_x_operation_xyz = self.atom_x_operation_xyz
        s_g_s_a_x = SpaceGroupSymop(operation_xyz=atom_x_operation_xyz)
        _item_a_s_x = atom_site[atom_x_label]
        fract_a_x_xyz = numpy.array([_item_a_s_x.fract_x, _item_a_s_x.fract_y,
                                     _item_a_s_x.fract_z], dtype=float)
        fract_a_x_xyz = numpy.dot(s_g_s_a_x.r, fract_a_x_xyz) + s_g_s_a_x.b
        pos_a_x_xyz = cell.calc_position_by_coordinate(
            fract_a_x_xyz[0], fract_a_x_xyz[1], fract_a_x_xyz[2])
        pos_a_x_xyz = numpy.array(pos_a_x_xyz, dtype=float)

        atom_y_label = self.atom_along_axis_y
        atom_y_operation_xyz = self.atom_y_operation_xyz
        s_g_s_a_y = SpaceGroupSymop(operation_xyz=atom_y_operation_xyz)
        _item_a_s_y = atom_site[atom_y_label]
        fract_a_y_xyz = numpy.array([_item_a_s_y.fract_x, _item_a_s_y.fract_y,
                                     _item_a_s_y.fract_z], dtype=float)
        fract_a_y_xyz = numpy.dot(s_g_s_a_y.r, fract_a_y_xyz) + s_g_s_a_y.b
        pos_a_y_xyz = cell.calc_position_by_coordinate(
            fract_a_y_xyz[0], fract_a_y_xyz[1], fract_a_y_xyz[2])
        pos_a_y_xyz = numpy.array(pos_a_y_xyz, dtype=float)

        v_pos_x = (pos_a_x_xyz-pos_a_o_xyz)/((
            numpy.square(pos_a_x_xyz-pos_a_o_xyz)).sum())**0.5
        v_pos_a_y = (pos_a_y_xyz-pos_a_o_xyz)/((
            numpy.square(pos_a_y_xyz-pos_a_o_xyz)).sum())**0.5
        v_pos_y_not_norm = v_pos_a_y - (v_pos_a_y * v_pos_x).sum() * v_pos_x
        v_pos_y = v_pos_y_not_norm / (
            numpy.square(v_pos_y_not_norm).sum())**0.5

        v_pos_z = numpy.cross(v_pos_x, v_pos_y)

        return v_pos_x, v_pos_y, v_pos_z

    def calc_fractions(self, cell: Cell, atom_site: AtomSiteL) -> \
            numpy.ndarray:
        """Give a numpy.nd_array of fractions: fractions_xyz.

        Arguments
        ---------
            - cell
            - atom_site
        """
        size_x, size_y = self.size_x, self.size_y
        points_x, points_y = self.points_x, self.points_y

        atom_o_label = self.atom_center
        _item_a_s_o = atom_site[atom_o_label]
        fract_a_o_xyz = numpy.array([_item_a_s_o.fract_x, _item_a_s_o.fract_y,
                                     _item_a_s_o.fract_z], dtype=float)
        pos_a_o_xyz = cell.calc_position_by_coordinate(
            fract_a_o_xyz[0], fract_a_o_xyz[1], fract_a_o_xyz[2])
        pos_a_o_xyz = numpy.array(pos_a_o_xyz, dtype=float)

        v_pos_x, v_pos_y, v_pos_z = self.calc_axes_x_y_z(cell, atom_site)

        v_delta_pos_x = v_pos_x * size_x / float(points_x)
        v_delta_pos_y = v_pos_y * size_y / float(points_y)
        np_x_2d, np_y_2d = numpy.meshgrid(range(-points_x//2, points_x//2),
                                          range(-points_y//2, points_y//2),
                                          indexing="ij")
        np_x_1d, np_y_1d = np_x_2d.flatten(), np_y_2d.flatten()

        pos_xyz = np_x_1d[numpy.newaxis, :]*v_delta_pos_x[:, numpy.newaxis] + \
            np_y_1d[numpy.newaxis, :] * v_delta_pos_y[:, numpy.newaxis] + \
            pos_a_o_xyz[:, numpy.newaxis]
        fract_x, fract_y, fract_z = cell.calc_coordinate_by_position(
            pos_xyz[0, :], pos_xyz[1, :], pos_xyz[2, :])
        return fract_x, fract_y, fract_z

    def calc_atoms(self, cell: Cell, atom_site: AtomSiteL,
                   space_group_symop: SpaceGroupSymopL, distance_min=0.3):
        """Calculate position of atoms in the section plane.

        Argument
        --------
            - cell
            - atom_site
            - space_group_symop
            - distance_min is minimal distance to search atom in angstrems
              (default value is 0.3)

        Output
        ------
            - atom_x
            - atom_y
            - atom_label
        """

        r_11 = numpy.array(space_group_symop.r_11, dtype=int)
        r_12 = numpy.array(space_group_symop.r_12, dtype=int)
        r_13 = numpy.array(space_group_symop.r_13, dtype=int)
        r_21 = numpy.array(space_group_symop.r_21, dtype=int)
        r_22 = numpy.array(space_group_symop.r_22, dtype=int)
        r_23 = numpy.array(space_group_symop.r_23, dtype=int)
        r_31 = numpy.array(space_group_symop.r_31, dtype=int)
        r_32 = numpy.array(space_group_symop.r_32, dtype=int)
        r_33 = numpy.array(space_group_symop.r_33, dtype=int)

        b_1 = numpy.array(space_group_symop.b_1, dtype=float)
        b_2 = numpy.array(space_group_symop.b_2, dtype=float)
        b_3 = numpy.array(space_group_symop.b_3, dtype=float)

        r_ij = (r_11, r_12, r_13, r_21, r_22, r_23, r_31, r_32, r_33)
        b_i = (b_1, b_2, b_3)

        fract_atom_auc_x = numpy.array(atom_site.fract_x, dtype=float)
        fract_atom_auc_y = numpy.array(atom_site.fract_y, dtype=float)
        fract_atom_auc_z = numpy.array(atom_site.fract_z, dtype=float)
        fract_atom_auc_xyz = (fract_atom_auc_x, fract_atom_auc_y,
                              fract_atom_auc_z)
        label_atom_auc = numpy.array(atom_site.label, dtype=str)

        fract_atom_uc_x, fract_atom_uc_y, fract_atom_uc_z, label_atom_uc = \
            calc_atoms_in_unit_cell(r_ij, b_i, fract_atom_auc_xyz,
                                    label_atom_auc)

        size_x = self.size_x
        size_y = self.size_y

        atom_center = self.atom_center
        atom_site_center = atom_site[atom_center]
        center_fract_x = atom_site_center.fract_x
        center_fract_y = atom_site_center.fract_y
        center_fract_z = atom_site_center.fract_z

        center_pos_x, center_pos_y, center_pos_z = \
            cell.calc_position_by_coordinate(center_fract_x, center_fract_y,
                                             center_fract_z)

        v_pos_x, v_pos_y, v_pos_z = self.calc_axes_x_y_z(cell, atom_site)

        pos_atom_uc_x, pos_atom_uc_y, pos_atom_uc_z = \
            cell.calc_position_by_coordinate(fract_atom_uc_x, fract_atom_uc_y,
                                             fract_atom_uc_z)

        pos_atom_loc_x = v_pos_x[0]*(pos_atom_uc_x - center_pos_x) + \
            v_pos_x[1]*(pos_atom_uc_y - center_pos_y) + \
            v_pos_x[2]*(pos_atom_uc_z - center_pos_z)
        pos_atom_loc_y = v_pos_y[0]*(pos_atom_uc_x - center_pos_x) + \
            v_pos_y[1]*(pos_atom_uc_y - center_pos_y) + \
            v_pos_y[2]*(pos_atom_uc_z - center_pos_z)
        pos_atom_loc_z = v_pos_z[0]*(pos_atom_uc_x - center_pos_x) + \
            v_pos_z[1]*(pos_atom_uc_y - center_pos_y) + \
            v_pos_z[2]*(pos_atom_uc_z - center_pos_z)

        flag_x = numpy.abs(pos_atom_loc_x) < 0.5*size_x
        flag_y = numpy.abs(pos_atom_loc_y) < 0.5*size_y
        flag_z = numpy.abs(pos_atom_loc_z) < distance_min
        flag_xyz = numpy.logical_and(flag_x, numpy.logical_and(flag_y, flag_z))

        atom_x = pos_atom_loc_x[flag_xyz]
        atom_y = pos_atom_loc_y[flag_xyz]
        atom_label = label_atom_uc[flag_xyz]
        return atom_x, atom_y, atom_label

class SectionL(LoopN):
    """Description of sections.

    Describe information concerning the density point.
    """

    ITEM_CLASS = Section
    ATTR_INDEX = "id"

    def __init__(self, loop_name=None) -> NoReturn:
        super(SectionL, self).__init__()
        self.__dict__["items"] = []
        self.__dict__["loop_name"] = loop_name

# s_cont = """
# loop_
# _section_id
# _section_size_x
# _section_size_y
# _section_atom_center
# _section_atom_along_axis_x
# _section_atom_along_axis_x_symop
# _section_atom_along_axis_y
# _section_atom_along_axis_y_symop
# _section_points_x
# _section_points_y
# _section_url_out
#   core 5 5 Ho1 O1 x,y,z O2 x,y,z 100 100 s_ho1.dat
#   """

# obj = SectionL.from_cif(s_cont)
# print(obj, end="\n\n")
# print(obj["core"], end="\n\n")

ext: py

max_stars_repo_name: eandklahn/cryspy

max_issues_repo_name: eandklahn/cryspy

max_forks_repo_name: eandklahn/cryspy

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
ES 1

"""
import numpy as np

"""

Punto A: la matrice A non ammette la fattorizzazione LU senza pivoting non avendo tutti i determinanti
            delle sottomatrici di testa != 0
            Metre la matrice B rispetta questa condizione dunque ammette tale fattorizzazione

"""

"""
Punto B
"""

def LUnoPivot(M):
    m, n = M.shape
    if m != n:
        print("Matrice non quadrata")
        return [], [], [], 1
    U = M.copy()
    for i in range(n):
        if U[i, i] == 0:
            print("Elemento diagonale nullo ")
            return [], [], [], 1
    for k in range(n - 1):
        U[k + 1 : n, k] = U[k + 1 : n, k] / U[k,k]
        U[k + 1 : n, k + 1] = U[k + 1 : n, k + 1] - (U[k + 1 : n, k] * U[k, k + 1 : n])
    L = np.tril(U, -1) + np.eye(n)
    U = np.triu(U)
    return L, U, np.eye(n), 0

B = np.array([[5, -2, 2, 0],
              [-2, 5, 0, 1],
              [2, 0, 5, 1],
              [0, 1, 1, 5]])
L, U, P, flag = LUnoPivot(B)
p1 = P.copy()
m, n = U.shape
p1[[0, n-1]] = p1[[n-1, 0]]
p1[[1, 2]] = p1[[2, 1]]
print(p1)
det = 1
det2 = 1
U2 = U * p1
for k in range(n):
    det *= U[k,k]


text:
#!/usr/bin/env python

# Copyright (c) 2018, Simon Brodeur
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors
#    may be used to endorse or promote products derived from this software without
#    specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
# NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.

import os
import logging
import numpy as np
import cv2
import scipy.signal

import matplotlib.pyplot as plt

from h5utils import Hdf5Dataset
from optparse import OptionParser

logger = logging.getLogger(__name__)

def is_cv2():
    import cv2 as lib
    return lib.__version__.startswith("2.")

def is_cv3():
    import cv2 as lib
    return lib.__version__.startswith("3.")

clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))

class FlowPlotter(object):

    def __init__(self, x, y, scale):

        h, w = x.shape
        self.x = x
        self.y = y

        # NOTE: the y-axis needs to be inverted to be in image-coordinates
        fig = plt.figure(figsize=(20,10), facecolor='white')
        ax = fig.add_subplot(121)
        q = ax.quiver(x, y, np.zeros((h,w)), np.zeros((h,w)), edgecolor='k', scale=1, angles='xy', scale_units='xy')
        ax.invert_yaxis()
        plt.axis('off')

        self.fig = fig
        self.ax = ax
        self.q = q

        ax2 = fig.add_subplot(122)
        m = ax2.imshow(np.zeros((int(h/scale), int(w/scale))), vmin = 0, vmax = 255, cmap = plt.get_cmap('gray'))
        plt.axis('off')
        self.ax2 = ax2
        self.m = m

        plt.ion()

    def update(self, flow, img):

        # NOTE: the y-axis needs to be negated to be in image-coordinates
        self.q.set_UVC(flow[:,:,0], -flow[:,:,1])
        self.m.set_data(img)
        self.fig.canvas.draw()

def computeSparseFlow(img, imgPrev, gridShape, border):

    # Define the fixed grid where optical flow is calculated
    h, w = img.shape[:2]
    y, x = np.meshgrid(np.linspace(border * h, (1.0 - border) * h, gridShape[0], dtype=np.int),
                       np.linspace(border * w, (1.0 - border) * w, gridShape[1], dtype=np.int),
                       indexing='ij')
    p0 = np.stack((y,x), axis=-1).astype(np.float32)

    # Calculate optical flow
    if is_cv2():
        p1, _, _ = cv2.calcOpticalFlowPyrLK(imgPrev, img, p0.reshape((-1,2)), None, winSize=(63,63), maxLevel=5)
    else:
        p1, _, _ = cv2.calcOpticalFlowPyrLK(imgPrev, img, p0.reshape((-1,2)), None, winSize=(63,63), maxLevel=5)

    flow = np.zeros((h, w, 2), dtype=np.float32)
    flow[y, x] = (p1.reshape(p0.shape) - p0)
    # TODO: test if all valid with p1[st==1] ?
    f = flow[y, x]

    return f, x, y

def computeCompactFlow(img, imgPrev, gridShape, border, filtering=False, filterSize=16):

    # Calculate optical flow
    if is_cv2():
        flow = cv2.calcOpticalFlowFarneback(imgPrev, img, pyr_scale=0.5, levels=5, winsize=63, iterations=3, poly_n=7, poly_sigma=1.5, flags=cv2.OPTFLOW_FARNEBACK_GAUSSIAN)
    else:
        flow = np.zeros((img.shape[0], img.shape[1], 2), dtype=np.float32)
        flow = cv2.calcOpticalFlowFarneback(imgPrev, img, flow, pyr_scale=0.5, levels=5, winsize=63, iterations=3, poly_n=7, poly_sigma=1.5, flags=cv2.OPTFLOW_FARNEBACK_GAUSSIAN)

    flow = np.array(flow, np.float32)

    h, w = flow.shape[:2]
    y, x = np.meshgrid(np.linspace(border * h, (1.0 - border) * h, gridShape[0], dtype=np.int),
                       np.linspace(border * w, (1.0 - border) * w, gridShape[1], dtype=np.int),
                       indexing='ij')

    if filtering:
        kernel = np.ones((filterSize, filterSize)) / float(filterSize * filterSize)
        flow = np.stack([scipy.signal.convolve2d(flow[:,:,k], kernel, mode='same', boundary='symm')
                         for k in range(flow.shape[-1])], axis=-1)
    f = flow[y, x]
    return f, x, y

def preprocessImage(img, scale=0.5):
    imggray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    imggray = clahe.apply(imggray)

    imggray = np.array(imggray * 255, dtype=np.uint8)
    imggray = cv2.resize(imggray, None,fx=scale, fy=scale, interpolation=cv2.INTER_LINEAR)
    return imggray

def processImageSequence(raw, clock, shape, scale=0.25, gridShape=(12,16), border=0.1, filtering=True, filterSize=5, visualize=False, stepSize=1):

    if visualize:
        plotter = None

    flows = []
    for i in range(stepSize, raw.shape[0]):

        logger.info('Processing image %d (%d total)' % (i+1, raw.shape[0]))

        img = cv2.imdecode(raw[i, :shape[i,0]], flags=-1) #CV_LOAD_IMAGE_UNCHANGED
        img = preprocessImage(img, scale)

        imgPrev = cv2.imdecode(raw[i-stepSize, :shape[0,0]], flags=-1) # CV_LOAD_IMAGE_UNCHANGED
        imgPrev = preprocessImage(imgPrev, scale)

        flow, x, y = computeCompactFlow(img, imgPrev, gridShape, border, filtering, filterSize)
        #flow, x, y = computeSparseFlow(img, imgPrev, gridShape, border)

        # Divide by the time interval and scale to get flow in unit of pixels/sec relative to the full-size image
        dt = clock[i] - clock[i-stepSize]
        assert dt > 0
        flow /=  (dt * scale)

        if visualize:
            if plotter is None:
                plotter = FlowPlotter(x, y, scale)
            plotter.update(flow, img)
            plt.show()

        flows.append(flow)

    flows = np.stack(flows, axis=0)
    clock = clock[stepSize:]

#     # Temporal filtering
#     if filtering:
#         kernel = np.ones((filterSize,)) / float(filterSize)
#         for i in range(flow.shape[0]):
#             for j in range(flow.shape[1]):
#                 flow[i,j,:] = scipy.signal.convolve(flow[i,j,:], kernel, mode='same')

    return flows, clock

def main(args=None):

    parser = OptionParser()
    parser.add_option("-i", "--input", dest="input", type='string', default=None,
                      help='specify the path of the input HDF5 file')
    parser.add_option("-o", "--output", dest="output", type='string', default=None,
                      help='specify the path of the output HDF5 file')
    parser.add_option("-s", "--scale", dest="scale", type='float', default=0.5,
                      help='specify the rescaling of images')
    parser.add_option("-b", "--border", dest="border", type='float', default=0.1,
                      help='specify the border to keep')
    parser.add_option("-f", "--filter-size", dest="filterSize", type='int', default=16,
                      help='specify the size of the smoothing filter')
    parser.add_option("-m", "--step-size", dest="stepSize", type='int', default=1,
                      help='specify the number of frame between optical flow calculation')
    parser.add_option("-y", "--grid-shape-y", dest="gridShape_y", type='int', default=12,
                  help='specify the shape of the optical flow on the y-axis')
    parser.add_option("-x", "--grid-shape-x", dest="gridShape_x", type='int', default=16,
                  help='specify the shape of the optical flow on the x-axis')
    parser.add_option("-t", "--filtering", action="store_true", dest="filtering", default=False,
                      help='specify to visualize computed optical flow')
    parser.add_option("-v", "--visualize", action="store_true", dest="visualize", default=False,
                      help='specify to visualize computed optical flow')
    (options,args) = parser.parse_args(args=args)

    inputDatasetPath = os.path.abspath(options.input)
    logger.info('Using input HDF5 dataset file: %s' % (inputDatasetPath))

    outputDatasetPath = os.path.abspath(options.output)
    logger.info('Using output HDF5 dataset file: %s' % (outputDatasetPath))

    with Hdf5Dataset(outputDatasetPath, mode='w') as outHdf5:
        with Hdf5Dataset(inputDatasetPath, mode='r') as inHdf5:

            # Process states
            for state in inHdf5.getAllStates():
                name, group, raw, clock, shape = state

                # Write original data to output file
                if group is not None:
                    ngroup = '/' + group
                else:
                    ngroup = ''
                fs = 1.0/np.mean(clock[1:] - clock[:-1])
                logger.info('Writing to output HDF5 dataset file: %s, shape=%s, fs=%f Hz' % (ngroup + '/' + name, str(raw.shape), fs))

                outHdf5.addStates(name, raw, clock, group, shape)

                if group == 'video':
                    rawFlow, clockFlow = processImageSequence(raw, clock, shape,
                                                              scale=options.scale, gridShape=(options.gridShape_y, options.gridShape_x), border=options.border,
                                                              filtering=options.filtering, filterSize=options.filterSize,
                                                              visualize=options.visualize, stepSize=options.stepSize)
                    groupFlow = 'optical_flow'

                    # Write original data to output file
                    fsFlow = 1.0/np.mean(clockFlow[1:] - clockFlow[:-1])
                    logger.info('Writing to output HDF5 dataset file: %s, shape=%s, fs=%f Hz' % ('/' + groupFlow + '/' + name, str(rawFlow.shape), fsFlow))

                    outHdf5.addStates(name, rawFlow, clockFlow, groupFlow, shape=None)

    logger.info('All done.')

if __name__ == '__main__':
    logging.basicConfig(level=logging.INFO)
    main()

max_stars_repo_name: sbrodeur/CREATE-dataset

max_issues_repo_name: sbrodeur/CREATE-dataset

max_forks_repo_name: sbrodeur/CREATE-dataset

reason: import numpy,import scipy


text:
import pytest
import numpy as np
from numpy.testing import assert_allclose

from pyvdms import util

def test_verify_tuple_range():
    out = util.verify_tuple_range(None, allow_none=True, todict=True)
    assert out == dict(first=None, last=None, step=None, unit=None)

    out = util.verify_tuple_range((0., 10., 1., 'km'), todict=True)
    assert out == dict(first=0., last=10., step=1., unit='km')

    out = util.verify_tuple_range((0., 10., 'km'), todict=True)
    assert out == dict(first=0., last=10., step=None, unit='km')

    out = util.verify_tuple_range((0., 10., 1.), todict=True)
    assert out == dict(first=0., last=10., step=1., unit=None)

    out = util.verify_tuple_range((0., 10., 1.), unit=False, todict=True)
    assert out == dict(first=0., last=10., step=1., unit=None)

    out = util.verify_tuple_range((0., 10., 'km'), step=False, todict=True)
    assert out == dict(first=0., last=10., step=None, unit='km')

def test_strlist_contains():
    strlist = ['foo', 'bar', 'Hello world']

    out = util.strlist_contains(None, '')
    assert out == False

    out = util.strlist_contains([], '')
    assert out == False

    out = util.strlist_contains(strlist, 'foo')
    assert out == True

    out = util.strlist_contains(strlist, 'hello')
    assert out == False

    out = util.strlist_contains(strlist, 'hello', False)
    assert out == True

def test_strlist_extract():
    out = util.strlist_extract(None, '')
    assert out == []

    out = util.strlist_extract([], '')
    assert out == []

    strlist = ['foo', 'bar', 'Foo bar']

    out = util.strlist_extract(strlist, 'foo', False)
    assert out == strlist[:1]

    out = util.strlist_extract(strlist, 'bar', False)
    assert out == strlist[1:]

    out = util.strlist_extract(strlist, 'Bar', False)
    assert out == []

    out = util.strlist_extract(strlist, 'Bar', False, case_sensitive=False)
    assert out == strlist[1:]

    out = util.strlist_extract(strlist, 'Foo', split=True)
    assert out == strlist[1:2]

    out = util.strlist_extract(strlist, 'foo', offset=1)
    assert out == strlist[1:2]

    out = util.strlist_extract(strlist, 'foo', offset=2, split=True)
    assert out == strlist[-1:]

    out = util.strlist_extract(strlist, 'foo', offset=-2, case_sensitive=False)
    assert out == strlist[:1]

max_stars_repo_name: psmsmets/pyVDMS

max_issues_repo_name: psmsmets/PyVDMS

max_forks_repo_name: psmsmets/PyVDMS

reason: import numpy,from numpy


text:
import time
import gc

class gridrad:
    """
    Name: GRIDRAD Python Class
    Purpose: This module contains three functions for dealing with Gridded NEXRAD WSR-88D Radar (GridRad).
    Authors:
    Cameron R. Homeyer
    Randy Chase

    """

    def __init__(self,filename=None,filter=False,toxr=False):
        self.filename=filename
        self.filter_flag=0
        self.ds = None

        if self.filename is not None:
            #auto read the file
            stime = time.time()
            self.read_file()
            print('Time to read file: {}'.format(time.time()-stime))

        if filter:
            stime = time.time()
            self.filter()
            print('Time to filter: {}'.format(time.time()-stime))
            stime = time.time()
            self.remove_clutter(skip_weak_ll_echo=1)
            print('Time to filter2: {}'.format(time.time()-stime))

        if toxr:
            self.to_xarray()

    # GridRad read routine
    def read_file(self):

        """
        A method that uses netCDF4 to read the GridRad file.
        This method also reshapes the GridRad data into 3d coordinates.

        Inputs:

        infile: str, path to nc4 file you wish to open.

        Returns: dict, raw data from the file.

        """

        # Import python libraries
        import sys
        import os
        import numpy as np
        import netCDF4

        # Check to see if file exists
        if not os.path.isfile(self.filename):
            print('File "' + self.filename + '" does not exist.  Returning -2.')
            return -2

        # Check to see if file has size of zero
        if os.stat(self.filename).st_size == 0:
            print('File "' + infile + '" contains no valid data.  Returning -1.')
            return -1

        from netCDF4 import Dataset
        from netCDF4 import Variable
        # Open GridRad netCDF file
        id = Dataset(self.filename, "r", format="NETCDF4")

        # Read global attributes
        Analysis_time           = str(id.getncattr('Analysis_time'          ))
        Analysis_time_window    = str(id.getncattr('Analysis_time_window'   ))
        File_creation_date      = str(id.getncattr('File_creation_date'     ))
        Grid_scheme             = str(id.getncattr('Grid_scheme'            ))
        Algorithm_version       = str(id.getncattr('Algorithm_version'      ))
        Algorithm_description   = str(id.getncattr('Algorithm_description'  ))
        Authors                 = str(id.getncattr('Authors'                ))
        Project_sponsor         = str(id.getncattr('Project_sponsor'        ))
        Project_name            = str(id.getncattr('Project_name'           ))

        # Read list of merged radar sweeps
        sweeps_list   = (id.variables['sweeps_merged'])[:]
        sweeps_merged = ['']*(id.dimensions['Sweep'].size)
        for i in range(0,id.dimensions['Sweep'].size):
            for j in range(0,id.dimensions['SweepRef'].size):
                sweeps_merged[i] += str(sweeps_list[i,j])

        # Read longitude dimension
        x = id.variables['Longitude']
        x = {'values'    : x[:],             \
              'long_name' : str(x.long_name), \
              'units'     : str(x.units),     \
              'delta'     : str(x.delta),     \
              'n'         : len(x[:])}

        # Read latitude dimension
        y = id.variables['Latitude']
        y = {'values'    : y[:],             \
              'long_name' : str(y.long_name), \
              'units'     : str(y.units),     \
              'delta'     : str(y.delta),     \
              'n'         : len(y[:])}

        # Read altitude dimension
        z = id.variables['Altitude']
        z = {'values'    : z[:],             \
              'long_name' : str(z.long_name), \
              'units'     : str(z.units),     \
              'delta'     : str(z.delta),     \
              'n'         : len(z[:])}

        # Read observation and echo counts
        nobs  = (id.variables['Nradobs' ])[:]
        necho = (id.variables['Nradecho'])[:]
        index = (id.variables['index'   ])[:]

        # Read reflectivity at horizontal polarization
        Z_H  = id.variables['Reflectivity' ]
        wZ_H = id.variables['wReflectivity']

        # Create arrays to store binned values for reflectivity at horizontal polarization
        values    = np.zeros(x['n']*y['n']*z['n'])
        wvalues   = np.zeros(x['n']*y['n']*z['n'])
        values[:] = float('nan')

        # Add values to arrays
        values[index[:]]  =  (Z_H)[:]
        wvalues[index[:]] = (wZ_H)[:]

        # Reshape arrays to 3-D GridRad domain
        values  =  values.reshape((z['n'], y['n'] ,x['n']))
        wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

        Z_H = {'values'     : values,              \
                 'long_name'  : str(Z_H.long_name),  \
                 'units'      : str(Z_H.units),      \
                 'missing'    : float('nan'),        \
                 'wvalues'    : wvalues,             \
                 'wlong_name' : str(wZ_H.long_name), \
                 'wunits'     : str(wZ_H.units),     \
                 'wmissing'   : wZ_H.missing_value,  \
                 'n'          : values.size}

        # Read velocity spectrum width
        SW  = id.variables['SpectrumWidth' ]
        wSW = id.variables['wSpectrumWidth']

        # Create arrays to store binned values for velocity spectrum width
        values    = np.zeros(x['n']*y['n']*z['n'])
        wvalues   = np.zeros(x['n']*y['n']*z['n'])
        values[:] = float('nan')

        # Add values to arrays
        values[index[:]]  =  (SW)[:]
        wvalues[index[:]] = (wSW)[:]

        # Reshape arrays to 3-D GridRad domain
        values  =  values.reshape((z['n'], y['n'] ,x['n']))
        wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

        SW  = {'values'     : values,             \
                 'long_name'  : str(SW.long_name),  \
                 'units'      : str(SW.units),      \
                 'missing'    : float('nan'),       \
                 'wvalues'    : wvalues,            \
                 'wlong_name' : str(wSW.long_name), \
                 'wunits'     : str(wSW.units),     \
                 'wmissing'   : wSW.missing_value,  \
                 'n'          : values.size}

        if ('AzShear' in id.variables):
            # Read azimuthal shear
            AzShr  = id.variables['AzShear' ]
            wAzShr = id.variables['wAzShear']

            # Create arrays to store binned values for azimuthal shear
            values    = np.zeros(x['n']*y['n']*z['n'])
            wvalues   = np.zeros(x['n']*y['n']*z['n'])
            values[:] = float('nan')

            # Add values to arrays
            values[index[:]]  =  (AzShr)[:]
            wvalues[index[:]] = (wAzShr)[:]

            # Reshape arrays to 3-D GridRad domain
            values  =  values.reshape((z['n'], y['n'] ,x['n']))
            wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

            AzShr = {'values'     : values,                \
                       'long_name'  : str(AzShr.long_name),  \
                       'units'      : str(AzShr.units),      \
                       'missing'    : float('nan'),          \
                       'wvalues'    : wvalues,               \
                       'wlong_name' : str(wAzShr.long_name), \
                       'wunits'     : str(wAzShr.units),     \
                       'wmissing'   : wAzShr.missing_value,  \
                       'n'          : values.size}

            # Read radial divergence
            Div  = id.variables['Divergence' ]
            wDiv = id.variables['wDivergence']

            # Create arrays to store binned values for radial divergence
            values    = np.zeros(x['n']*y['n']*z['n'])
            wvalues   = np.zeros(x['n']*y['n']*z['n'])
            values[:] = float('nan')

            # Add values to arrays
            values[index[:]]  =  (Div)[:]
            wvalues[index[:]] = (wDiv)[:]

            # Reshape arrays to 3-D GridRad domain
            values  =  values.reshape((z['n'], y['n'] ,x['n']))
            wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

            Div = {'values'     : values,              \
                     'long_name'  : str(Div.long_name),  \
                     'units'      : str(Div.units),      \
                     'missing'    : float('nan'),        \
                     'wvalues'    : wvalues,             \
                     'wlong_name' : str(wDiv.long_name), \
                     'wunits'     : str(wDiv.units),     \
                     'wmissing'   : wDiv.missing_value,  \
                     'n'          : values.size}

        else:
            AzShr = -1
            Div   = -1

        if ('DifferentialReflectivity' in id.variables):
            # Read radial differential reflectivity
            Z_DR  = id.variables['DifferentialReflectivity' ]
            wZ_DR = id.variables['wDifferentialReflectivity']

            # Create arrays to store binned values for differential reflectivity
            values    = np.zeros(x['n']*y['n']*z['n'])
            wvalues   = np.zeros(x['n']*y['n']*z['n'])
            values[:] = float('nan')

            # Add values to arrays
            values[index[:]]  =  (Z_DR)[:]
            wvalues[index[:]] = (wZ_DR)[:]

            # Reshape arrays to 3-D GridRad domain
            values  =  values.reshape((z['n'], y['n'] ,x['n']))
            wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

            Z_DR = {'values'     : values,               \
                      'long_name'  : str(Z_DR.long_name),  \
                      'units'      : str(Z_DR.units),      \
                      'missing'    : float('nan'),         \
                      'wvalues'    : wvalues,              \
                      'wlong_name' : str(wZ_DR.long_name), \
                      'wunits'     : str(wZ_DR.units),     \
                      'wmissing'   : wZ_DR.missing_value,  \
                      'n'          : values.size}

            # Read specific differential phase
            K_DP  = id.variables['DifferentialPhase' ]
            wK_DP = id.variables['wDifferentialPhase']

            # Create arrays to store binned values for specific differential phase
            values    = np.zeros(x['n']*y['n']*z['n'])
            wvalues   = np.zeros(x['n']*y['n']*z['n'])
            values[:] = float('nan')

            # Add values to arrays
            values[index[:]]  =  (K_DP)[:]
            wvalues[index[:]] = (wK_DP)[:]

            # Reshape arrays to 3-D GridRad domain
            values  =  values.reshape((z['n'], y['n'] ,x['n']))
            wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

            K_DP = {'values'     : values,               \
                      'long_name'  : str(K_DP.long_name),  \
                      'units'      : str(K_DP.units),      \
                      'missing'    : float('nan'),         \
                      'wvalues'    : wvalues,              \
                      'wlong_name' : str(wK_DP.long_name), \
                      'wunits'     : str(wK_DP.units),     \
                      'wmissing'   : wK_DP.missing_value,  \
                      'n'          : values.size}

            # Read correlation coefficient
            r_HV  = id.variables['CorrelationCoefficient' ]
            wr_HV = id.variables['wCorrelationCoefficient']

            # Create arrays to store binned values for correlation coefficient
            values    = np.zeros(x['n']*y['n']*z['n'])
            wvalues   = np.zeros(x['n']*y['n']*z['n'])
            values[:] = float('nan')

            # Add values to arrays
            values[index[:]]  =  (r_HV)[:]
            wvalues[index[:]] = (wr_HV)[:]

            # Reshape arrays to 3-D GridRad domain
            values  =  values.reshape((z['n'], y['n'] ,x['n']))
            wvalues = wvalues.reshape((z['n'], y['n'] ,x['n']))

            r_HV = {'values'     : values,               \
                      'long_name'  : str(r_HV.long_name),  \
                      'units'      : str(r_HV.units),      \
                      'missing'    : float('nan'),         \
                      'wvalues'    : wvalues,              \
                      'wlong_name' : str(wr_HV.long_name), \
                      'wunits'     : str(wr_HV.units),     \
                      'wmissing'   : wr_HV.missing_value,  \
                      'n'          : values.size}

        else:
            Z_DR = -1
            K_DP = -1
            r_HV = -1

        # Close netCDF4 file
        id.close()
        del id
        gc.collect()

        #store data dict. in the class obj
        self.data =  {'name'                    : 'GridRad analysis for ' + Analysis_time, \
                  'x'                       : x, \
                  'y'                       : y, \
                  'z'                       : z, \
                  'Z_H'                     : Z_H, \
                  'SW'                      : SW, \
                  'AzShr'                   : AzShr, \
                  'Div'                     : Div, \
                  'Z_DR'                    : Z_DR, \
                  'K_DP'                    : K_DP, \
                  'r_HV'                    : r_HV, \
                  'nobs'                    : nobs, \
                  'necho'                   : necho, \
                  'file'                    : self.filename, \
                  'sweeps_merged'           : sweeps_merged, \
                  'Analysis_time'           : Analysis_time, \
                  'Analysis_time_window'    : Analysis_time_window, \
                  'File_creation_date'      : File_creation_date, \
                  'Grid_scheme'             : Grid_scheme, \
                  'Algorithm_version'       : Algorithm_version, \
                  'Algorithm_description'   : Algorithm_description, \
                  'Authors'                 : Authors, \
                  'Project_sponsor'         : Project_sponsor, \
                  'Project_name'            : Project_name}

    def filter(self):
        import copy
        if self.filter_flag == 0:
            data0 = copy.deepcopy(self.data)
            #del copy of data to clear out RAM
            del self.data
            gc.collect()
        elif self.filter_flag > 0:
            data0 = copy.deepcopy(self.data_filt)
            #del copy of data to clear out RAM
            del self.data_filt
            gc.collect()

        # Import python libraries
        import sys
        import os
        import numpy as np

        #Extract year from GridRad analysis time string
        year = int((data0['Analysis_time'])[0:4])

        wthresh     = 1.5												# Set default bin weight threshold for filtering by year (dimensionless)
        freq_thresh = 0.6												# Set echo frequency threshold (dimensionless)
        Z_H_thresh  = 15.0											# Reflectivity threshold (dBZ)
        nobs_thresh = 2												# Number of observations threshold

        # Extract dimension sizes
        nx = (data0['x'])['n']
        ny = (data0['y'])['n']
        nz = (data0['z'])['n']

        echo_frequency = np.zeros((nz,ny,nx))					# Create array to compute frequency of radar obs in grid volume with echo

        ipos = np.where(data0['nobs'] > 0)						# Find bins with obs
        npos = len(ipos[0])											# Count number of bins with obs

        if (npos > 0):
            echo_frequency[ipos] = (data0['necho'])[ipos]/(data0['nobs'])[ipos]		# Compute echo frequency (number of scans with echo out of total number of scans)

        inan = np.where(np.isnan((data0['Z_H'])['values']))				# Find bins with NaNs
        nnan = len(inan[0])														# Count number of bins with NaNs

        if (nnan > 0): ((data0['Z_H'])['values'])[inan] = 0.0

        # Find observations with low weight
        ifilter = np.where((((data0['Z_H'])['wvalues'] < wthresh    ) & ((data0['Z_H'])['values'] < Z_H_thresh)) |
                                  ((echo_frequency           < freq_thresh) &  (data0['nobs'] > nobs_thresh)))

        nfilter = len(ifilter[0])									# Count number of bins that need to be removed

        # Remove low confidence observations
        if (nfilter > 0):
            ((data0['Z_H'])['values'])[ifilter] = float('nan')
            ((data0['SW' ])['values'])[ifilter] = float('nan')

            if (type(data0['AzShr']) is dict):
                ((data0['AzShr'])['values'])[ifilter] = float('nan')
                ((data0['Div'  ])['values'])[ifilter] = float('nan')

            if (type(data0['Z_DR']) is dict):
                ((data0['Z_DR'])['values'])[ifilter] = float('nan')
                ((data0['K_DP'])['values'])[ifilter] = float('nan')
                ((data0['r_HV'])['values'])[ifilter] = float('nan')

        # Replace NaNs that were previously removed
        if (nnan > 0): ((data0['Z_H'])['values'])[inan] = float('nan')

        #free up some RAM
        del ipos,inan,ifilter,echo_frequency
        gc.collect()

        #this step dies with super big domains
        self.data_filt = copy.deepcopy(data0)

        #clear out some RAM
        del data0
        gc.collect()

        self.filter_flag = 1

    def remove_clutter(self,skip_weak_ll_echo=0):

        import copy
        if self.filter_flag == 0:
            data0 = copy.deepcopy(self.data)
            #del copy of data to clear out RAM
            del self.data
            gc.collect()
        elif self.filter_flag > 0:
            data0 = copy.deepcopy(self.data_filt)
            #del copy of data to clear out RAM
            del self.data_filt
            gc.collect()

        # Import python libraries
        import sys
        import os
        import numpy as np

        # Set fractional areal coverage threshold for speckle identification
        areal_coverage_thresh = 0.32

        # Extract dimension sizes
        nx = (data0['x'])['n']
        ny = (data0['y'])['n']
        nz = (data0['z'])['n']

        # Copy altitude array to 3 dimensions
        zzz = ((((data0['z'])['values']).reshape(nz,1,1)).repeat(ny, axis = 1)).repeat(nx, axis = 2)

        # Light pass at a correlation coefficient decluttering approach first
        if (type(data0['Z_DR']) is dict):
            ibad = np.where((((data0['Z_H'])['values'] < 40.0) & ((data0['r_HV'])['values'] < 0.8)) | \
                          (((data0['Z_H'])['values'] < 25.0) & ((data0['r_HV'])['values'] < 0.9) & (zzz >= 10.0)))
            nbad = len(ibad[0])

            if (nbad > 0):
                ((data0['Z_H' ])['values'])[ibad] = float('nan')
                ((data0['SW'  ])['values'])[ibad] = float('nan')
                ((data0['Z_DR'])['values'])[ibad] = float('nan')
                ((data0['K_DP'])['values'])[ibad] = float('nan')
                ((data0['r_HV'])['values'])[ibad] = float('nan')

                if (type(data0['AzShr']) is dict):
                    ((data0['AzShr'])['values'])[ibad] = float('nan')
                    ((data0['Div'  ])['values'])[ibad] = float('nan')

        # First pass at removing speckles
        fin = np.isfinite((data0['Z_H'])['values'])

        # Compute fraction of neighboring points with echo
        cover = np.zeros((nz,ny,nx))
        for i in range(-2,3):
            for j in range(-2,3):
                cover += np.roll(np.roll(fin, i, axis=2), j, axis=1)
        cover = cover/25.0

        # Find bins with low nearby areal echo coverage (i.e., speckles) and remove (set to NaN).
        ibad = np.where(cover <= areal_coverage_thresh)
        nbad = len(ibad[0])
        if (nbad > 0):
            ((data0['Z_H'])['values'])[ibad] = float('nan')
            ((data0['SW' ])['values'])[ibad] = float('nan')

            if (type(data0['AzShr']) is dict):
                ((data0['AzShr'])['values'])[ibad] = float('nan')
                ((data0['Div'  ])['values'])[ibad] = float('nan')

            if (type(data0['Z_DR']) is dict):
                ((data0['Z_DR'])['values'])[ibad] = float('nan')
                ((data0['K_DP'])['values'])[ibad] = float('nan')
                ((data0['r_HV'])['values'])[ibad] = float('nan')

        # Attempts to mitigate ground clutter and biological scatterers
        if (skip_weak_ll_echo == 0):
            # First check for weak, low-level echo
            inan = np.where(np.isnan((data0['Z_H'])['values']))				# Find bins with NaNs
            nnan = len(inan[0])															# Count number of bins with NaNs

            if (nnan > 0): ((data0['Z_H'])['values'])[inan] = 0.0

            # Find weak low-level echo and remove (set to NaN)
            ibad = np.where(((data0['Z_H'])['values'] < 10.0) & (zzz <= 4.0))
            nbad = len(ibad[0])
            if (nbad > 0):
                ((data0['Z_H'])['values'])[ibad] = float('nan')
                ((data0['SW' ])['values'])[ibad] = float('nan')

                if (type(data0['AzShr']) is dict):
                    ((data0['AzShr'])['values'])[ibad] = float('nan')
                    ((data0['Div'  ])['values'])[ibad] = float('nan')

                if (type(data0['Z_DR']) is dict):
                    ((data0['Z_DR'])['values'])[ibad] = float('nan')
                    ((data0['K_DP'])['values'])[ibad] = float('nan')
                    ((data0['r_HV'])['values'])[ibad] = float('nan')

            # Replace NaNs that were removed
            if (nnan > 0): ((data0['Z_H'])['values'])[inan] = float('nan')

            # Second check for weak, low-level echo
            inan = np.where(np.isnan((data0['Z_H'])['values']))				# Find bins with NaNs
            nnan = len(inan[0])															# Count number of bins with NaNs

            if (nnan > 0): ((data0['Z_H'])['values'])[inan] = 0.0

            refl_max   = np.nanmax( (data0['Z_H'])['values'],             axis=0)
            echo0_max  = np.nanmax(((data0['Z_H'])['values'] >  0.0)*zzz, axis=0)
            echo0_min  = np.nanmin(((data0['Z_H'])['values'] >  0.0)*zzz, axis=0)
            echo5_max  = np.nanmax(((data0['Z_H'])['values'] >  5.0)*zzz, axis=0)
            echo15_max = np.nanmax(((data0['Z_H'])['values'] > 15.0)*zzz, axis=0)

            # Replace NaNs that were removed
            if (nnan > 0): ((data0['Z_H'])['values'])[inan] = float('nan')

            # Find weak and/or shallow echo
            ibad = np.where(((refl_max   <  20.0) & (echo0_max  <= 4.0) & (echo0_min  <= 3.0)) | \
                                 ((refl_max   <  10.0) & (echo0_max  <= 5.0) & (echo0_min  <= 3.0)) | \
                                 ((echo5_max  <=  5.0) & (echo5_max  >  0.0) & (echo15_max <= 3.0)) | \
                                 ((echo15_max <   2.0) & (echo15_max >  0.0)))
            nbad = len(ibad[0])
            if (nbad > 0):
                kbad = (np.zeros((nbad))).astype(int)
                for k in range(0,nz):
                    ((data0['Z_H'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                    ((data0['SW' ])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')

                    if (type(data0['AzShr']) is dict):
                        ((data0['AzShr'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                        ((data0['Div'  ])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')

                    if (type(data0['Z_DR']) is dict):
                        ((data0['Z_DR'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                        ((data0['K_DP'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                        ((data0['r_HV'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')

        # Find clutter below convective anvils
        k4km = ((np.where((data0['z'])['values'] >= 4.0))[0])[0]
        fin  = np.isfinite((data0['Z_H'])['values'])
        ibad = np.where((          fin[k4km         ,:,:]          == 0) & \
                                 (np.sum(fin[k4km:(nz  -1),:,:], axis=0) >  0) & \
                                 (np.sum(fin[   0:(k4km-1),:,:], axis=0) >  0))
        nbad = len(ibad[0])
        if (nbad > 0):
            kbad = (np.zeros((nbad))).astype(int)
            for k in range(0,k4km+1):
                ((data0['Z_H'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                ((data0['SW' ])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')

                if (type(data0['AzShr']) is dict):
                    ((data0['AzShr'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                    ((data0['Div'  ])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')

                if (type(data0['Z_DR']) is dict):
                    ((data0['Z_DR'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                    ((data0['K_DP'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')
                    ((data0['r_HV'])['values'])[(k+kbad),ibad[0],ibad[1]] = float('nan')

        # Second pass at removing speckles
        fin = np.isfinite((data0['Z_H'])['values'])

        # Compute fraction of neighboring points with echo
        cover = np.zeros((nz,ny,nx))
        for i in range(-2,3):
            for j in range(-2,3):
                cover += np.roll(np.roll(fin, i, axis=2), j, axis=1)
        cover = cover/25.0

        # Find bins with low nearby areal echo coverage (i.e., speckles) and remove (set to NaN).
        ibad = np.where(cover <= areal_coverage_thresh)
        nbad = len(ibad[0])
        if (nbad > 0):
            ((data0['Z_H'])['values'])[ibad] = float('nan')
            ((data0['SW' ])['values'])[ibad] = float('nan')

            if (type(data0['AzShr']) is dict):
                ((data0['AzShr'])['values'])[ibad] = float('nan')
                ((data0['Div'  ])['values'])[ibad] = float('nan')

            if (type(data0['Z_DR']) is dict):
                ((data0['Z_DR'])['values'])[ibad] = float('nan')
                ((data0['K_DP'])['values'])[ibad] = float('nan')
                ((data0['r_HV'])['values'])[ibad] = float('nan')

        #free up some RAM
        del zzz,ibad,cover,fin,k4km
        gc.collect()

        self.data_filt = copy.deepcopy(data0)

        #del data0 to free up RAM
        del data0
        gc.collect()

        self.filter_flag = 2

    def to_xarray(self, key_list_in = ['Z_H','Z_DR','K_DP','r_HV','AzShr','Div','SW'],
                key_list_out = ['ZH','ZDR','KDP','RHV','VOR','DIV','SW']):
        for i,ii in enumerate(key_list_in):
            self._append_data_array(variable_key_in=ii,variable_key_out=key_list_out[i])

    def _append_data_array(self,variable_key_in='Z_H',variable_key_out='ZH'):

        if self.filter_flag == 0:
            data = self.data
        else:
            data = self.data_filt

        #need to convert time
        from pandas import to_datetime
        from netCDF4 import date2num
        time = date2num(to_datetime(data['Analysis_time']),
                          'seconds since 2001-01-01 00:00:00')

        import xarray as xr
        import numpy as np
        if (variable_key_in =='Z_DR') or (variable_key_in =='K_DP') or (variable_key_in =='r_HV'):
                if type(data[variable_key_in]) == int:
                        da = xr.DataArray(np.ones(data['Z_H']['values'][np.newaxis,:,:,:].shape)*-9999,
                          dims=['time','Altitude','Latitude','Longitude'],
                          coords={'time': [time],
                                  'Longitude': data['x']['values'],
                                  'Latitude': data['y']['values'],
                                  'Altitude': data['z']['values']})
                else:
                        da = xr.DataArray(data[variable_key_in]['values'][np.newaxis,:,:,:],
                          dims=['time','Altitude','Latitude','Longitude'],
                          coords={'time': [time],
                                  'Longitude': data['x']['values'],
                                  'Latitude': data['y']['values'],
                                  'Altitude': data['z']['values']})
        else:
            da = xr.DataArray(data[variable_key_in]['values'][np.newaxis,:,:,:],
                              dims=['time','Altitude','Latitude','Longitude'],
                              coords={'time': [time],
                                      'Longitude': data['x']['values'],
                                      'Latitude': data['y']['values'],
                                      'Altitude': data['z']['values']})
        da.fillna(value=-9999)

        if self.ds is None:
            self.ds = da.to_dataset(name = variable_key_out)
        else:
            self.ds[variable_key_out] = da

max_stars_repo_name: dopplerchase/GewitterGefahr

max_issues_repo_name: dopplerchase/GewitterGefahr

max_forks_repo_name: dopplerchase/GewitterGefahr


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np

'''create DataFrame
DataFrame 数据桢,数据表,;类似于excel
特点:
1. 他是Series的集合
2. 与Series的区别:
    2.1 series吧通过自定义index,当做标记,实现行,一维列表
    2.2 DataFrame通过在'吧自定义index当做标记实现行'上与Series是一致,
    2.3 DataFrame 除了行以外,还提供了columns(列), 每一列都是一个Series,所以DataFrame是Series的集合,通过colums将Series集合在一起
3. DataFrame 支持花式索引, 以及提供API方便我们处理数据
4. 每一个DataFrame(表),内部的数据最好是具有相同columns的数据, 行:往外代表数据的多少,列:往往代表数据的特点和结构
5. DataFrame: 每一行代表着一组完整的数据集,其中每一个数据有自己的属性.每一个columns 应该是一种数据类型,代表这该数据结构的每一个属性或数据类型
'''

# index:datatime_list, columns:list, data:ndarray
dates = pd.date_range('20130101', periods=6)
df = pd.DataFrame(np.random.randn(6, 4), index=dates, columns=list('ABCD'))
print(df)

# 通过dict 创建对象, 且个一个columns数据类型不同, 自动填充数据:
# data:dict_value, index:auto_create_index, columns:dict_key

data1 = {
    'A': 1.,
    'B': pd.Timestamp('20130102'),
    'C': pd.Series(1, index=list(range(5)), dtype='float32'),
    'D': np.array([3] * 5, dtype='int32'),
    'E': pd.Categorical(['test', 'train', 'test', 'train', 'train']),
    'F': 'foo'
}

df2 = pd.DataFrame(data1)
print(df2)
print(df2.dtypes)


text:
import numpy as np

from .metric import Metric
from nn_analysis import metrics as me
from nn_analysis import acts as ac

class Factorization(Metric):
    def __init__(self, variable_configs):
        self.variable_configs = variable_configs

    @property
    def required_acts(self):
        acts = []
        for _, variable_config in self.variable_configs.items():
            acts.append(variable_config['acts'])
        return acts

    def _evaluate_single(self, model_name, epoch, layer_name, variable, acts, n_pcs, metric_types, kwargs):
        result = {}

        X = ac.utils.load_data(model_name, epoch, acts['name'], acts['version'], layer_name=layer_name)[...,:n_pcs] # (...,n_pcs)
        scores = me.core.compute_factorization_metrics(X, **kwargs)

        for metric_type in metric_types:
            result[f'{metric_type}-{variable}'] = scores[metric_type]

        return result

    def evaluate(self, model_name, epoch, layer_name):
        result = {}
        for variable, variable_config in self.variable_configs.items():
            print(f"Evaluating variable {variable}")
            single_result = self._evaluate_single(model_name, epoch, layer_name, variable, **variable_config)
            result.update(single_result)

        return result


text:
import matplotlib.pyplot as plt
import numpy as np
import sys, os
sys.path.insert(0, "..")
import torch
from torch.utils.data import Dataset, DataLoader
from torchvision import transforms, utils
from datasets import PlanningDataset
from diff_gpmp2.env.env_2d import Env2D

plt.ion()
dataset = PlanningDataset("../datasets/dataset_files/dataset_2d_8/", mode='train', label_subdir='opt_trajs_gpmp2')
print len(dataset)

##Iterate through the dataset and visualize
for i in xrange(len(dataset)):
  print dataset[i]
  im = dataset[i]['im'][0,:,:]
  sdf = dataset[i]['sdf'][:,:]
  start = dataset[i]['start']
  goal = dataset[i]['goal']
  th_opt = dataset[i]['th_opt']
  env_params = dataset.meta_data['env_params']
  env = Env2D(env_params)
  env.initialize_from_image(im, sdf)
  path_f = []
  for j in xrange(th_opt.shape[0]):
    path_f.append(th_opt[j,0:2])
  print start[0][0:2], goal[0][0:2]
  env.initialize_plot(start[0][0:2], goal[0][0:2])
  env.plot_edge(path_f)
  # env.plot_signed_distance_transform()
  plt.show()
  raw_input('Curr datapoint = {}. Press enter to view next data point'.format(i))
  env.close_plot()

# def visualize_batch(sample_batch):
#   im_batch, sdf_batch, start_batch, goal_batch, th_opt_batch = \
#     sample_batch['im'], sample_batch['sdf'], sample_batch['start'], sample_batch['goal'], sample_batch['th_opt']
#   batch_size = len(im_batch)
#   im_size = im_batch.size(2)
#   sdf_batch = sdf_batch.unsqueeze(1)
#   grid_im = utils.make_grid(im_batch)
#   grid_sdf = utils.make_grid(sdf_batch)

#   plt.figure()
#   plt.imshow(grid_im.numpy().transpose((1,2,0)))
#   plt.figure()
#   plt.imshow(grid_sdf.numpy().transpose((1,2,0)))

# #Let's test batching
# dataloader = DataLoader(dataset, batch_size=4,
#                         shuffle=True, num_workers=4)

# for i_batch, sample_batched in enumerate(dataloader):
#   print(i_batch, sample_batched['im'].size(),
#         sample_batched['sdf'].size(), sample_batched['start'].size(),
#         sample_batched['goal'].size(), sample_batched['th_opt'].size())

#   if i_batch == 3:
#     visualize_batch(sample_batched)
#     plt.axis('off')
#     plt.ioff()
#     plt.show()
#     break


text:
import os
import random
import time
import traceback
from csv import reader
from os import path
from queue import Queue
from threading import Thread
from typing import Tuple, List, Any, Optional

import numpy as np
import tensorflow as tf
from tensorflow import Session
from tensorflow.summary import FileWriter
from tensorflow.train import Saver

from deepword.agents.utils import bert_nlu_input
from deepword.students.student_learner import NLUClassificationLearner
from deepword.students.utils import align_batch_str
from deepword.tokenizers import Tokenizer

tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.FATAL)
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'  # FATAL

def process_one_line(
        tokenizer: Tokenizer,
        start_str: str,
        ending_str: List[str],
        max_start_len: int,
        max_end_len: int) -> Tuple[List[int], int, np.ndarray, np.ndarray]:
    tj_ids = tokenizer.convert_tokens_to_ids(tokenizer.tokenize(start_str))
    tj_len = len(tj_ids)
    aligned_tj = align_batch_str(
        [tj_ids], str_len_allowance=max_start_len, padding_val_id=0,
        valid_len=[tj_len])
    tj_ids = aligned_tj[0][0]
    tj_len = aligned_tj[1][0]

    actions = [
        tokenizer.convert_tokens_to_ids(tokenizer.tokenize(a))
        for a in ending_str]
    action_matrix, action_len = align_batch_str(
        actions, str_len_allowance=max_end_len, padding_val_id=0,
        valid_len=[len(a) for a in actions])

    return tj_ids, tj_len, action_matrix, action_len

def load_swag_data(
        fn_swag: str) -> Tuple[List[str], List[np.ndarray], List[int]]:
    """read swag data, make sure every start sentence has four ends"""
    with open(fn_swag, "r") as f:
        lines = list(reader(f.readlines()))[1:]  # remove the first line

    assert len(lines[0]) == 12, "use train.csv or val.csv which contain labels"

    lines = np.asarray(lines, dtype=np.object)
    labels = list(lines[:, -1].astype(np.int32))
    start_str = list(lines[:, 3])
    ending_str = list(lines[:, 7:11])

    return start_str, ending_str, labels

def get_bert_input(
        start_str: List[str], ending_str: List[np.ndarray],
        tokenizer: Tokenizer, sep_val_id, cls_val_id,
        num_tokens, n_tokens_per_action
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    TODO: here we use n_tokens_per_actions for end str size
    """
    inp = []
    seg_tj_action = []
    inp_size = []
    num_tokens_tj = num_tokens - n_tokens_per_action - 3

    for tj, actions in zip(start_str, ending_str):
        tj_ids, tj_len, action_matrix, action_len = process_one_line(
            tokenizer, tj, list(actions), num_tokens_tj, n_tokens_per_action)
        _inp, _seg_tj_action, _inp_size = bert_nlu_input(
            action_matrix, action_len, tj_ids, tj_len,
            sep_val_id, cls_val_id, num_tokens)
        inp.append(_inp)
        seg_tj_action.append(_seg_tj_action)
        inp_size.append(_inp_size)

    inp = np.concatenate(inp, axis=0)
    seg_tj_action = np.concatenate(seg_tj_action, axis=0)
    inp_size = np.concatenate(inp_size, axis=0)
    return inp, seg_tj_action, inp_size

class SwagLearner(NLUClassificationLearner):
    def _test_impl(self, data: Tuple) -> np.ndarray:
        """
        Dummy method, won't be used.
        """
        raise NotImplementedError()

    def _prepare_training(
            self
    ) -> Tuple[Session, Any, Saver, FileWriter, int, Queue]:
        sess, model, saver, train_steps = self._prepare_model(
            "/device:GPU:0", training=True)

        # save the very first model to verify weight has been loaded
        if train_steps == 0:
            saver.save(
                sess, self.ckpt_prefix,
                global_step=tf.train.get_or_create_global_step(
                    graph=model.graph))
        else:
            pass

        sw_path = path.join(self.model_dir, "summaries", "train")
        sw = tf.summary.FileWriter(sw_path, sess.graph)

        queue = Queue(maxsize=100)
        t = Thread(
            target=self._add_batch,
            args=("{}/data/train.csv".format(self.train_data_dir), queue))
        t.setDaemon(True)
        t.start()

        return sess, model, saver, sw, train_steps, queue

    def _prepare_test(
            self, device_placement: str = "/device:GPU:0",
            restore_from: Optional[str] = None
    ) -> Tuple[Session, Any, Saver, int, Queue]:
        sess, model, saver, train_steps = self._prepare_model(
            device_placement, training=False, restore_from=restore_from)
        queue = Queue(maxsize=100)
        t = Thread(
            target=self._add_batch,
            args=(self.eval_data_path, queue, False))
        t.setDaemon(True)
        t.start()
        return sess, model, saver, train_steps, queue

    def test(
            self, device_placement: str = "/device:GPU:0",
            restore_from: Optional[str] = None) -> Tuple[float, int]:
        if self.sess is None:
            (self.sess, self.model, self.saver, self.train_steps, self.queue
             ) = self._prepare_test(device_placement, restore_from)

        wait_times = 10
        while wait_times > 0 and self.queue.empty():
            self.info("waiting data ... (retry times: {})".format(wait_times))
            time.sleep(10)
            wait_times -= 1

        acc = 0
        total = 0
        self.info("start test")
        i = 0
        for data in iter(self.queue.get, None):
            inp, seg_tj_action, inp_len, _, swag_labels = data

            if i % 100 == 0:
                self.debug("process a batch of {} .. {}".format(len(inp), i))
                self.debug("partial acc.: {}".format(
                    acc * 1. / total if total else "Nan"))

            q_actions_t = self.sess.run(
                self.model.q_actions,
                feed_dict={
                    self.model.src_: inp,
                    self.model.src_len_: inp_len,
                    self.model.seg_tj_action_: seg_tj_action
                })
            q_actions_t = q_actions_t.reshape((len(swag_labels), -1))
            predicted_swag_labels = np.argmax(q_actions_t, axis=-1)
            acc += np.count_nonzero(predicted_swag_labels == swag_labels)
            total += len(swag_labels)
            i += 1
        return acc * 1. / total, total

    def _add_batch(
            self, swag_path: str,
            queue: Queue, training: bool = True,
            append_new_data: bool = True) -> None:
        start_str, ending_str, labels = load_swag_data(swag_path)
        data = list(zip(start_str, ending_str, labels))
        while True:
            if training:
                random.shuffle(data)
            i = 0
            while i < len(data) // self.hp.batch_size:
                ss = i * self.hp.batch_size
                ee = min((i + 1) * self.hp.batch_size, len(data))
                try:
                    batch_data = data[ss:ee]
                    batch_start_str = [x[0] for x in batch_data]
                    batch_ending_str = [x[1] for x in batch_data]
                    batch_labels = [x[2] for x in batch_data]

                    inp, seg_tj_action, inp_size = get_bert_input(
                        batch_start_str, batch_ending_str, self.tokenizer,
                        self.hp.sep_val_id, self.hp.cls_val_id,
                        self.hp.num_tokens, self.hp.n_tokens_per_action)
                    queue.put(
                        (inp, seg_tj_action, inp_size, None, batch_labels))
                except Exception as e:
                    self.error("add_batch error: {}".format(e))
                    traceback.print_tb(e.__traceback__)
                    raise RuntimeError()
                i += 1

            if not training:
                self.queue.put(None)
                break


text:
#!/usr/bin/env python

"""
Reference:
https://ps1images.stsci.edu/ps1image.html
"""
import requests
import numpy as np
import matplotlib.pyplot as pl
from astropy import units as u
from astropy.table import Table
from PIL import Image
from io import BytesIO
from astropy.io import fits

__all__ = ['Panstarrs']

PIXEL_SCALE = 0.25*u.arcsec/u.pix

class Panstarrs:
    def __init__(self, ra, dec, fov=30*u.arcsec, filters="grizy", color=False,
                 format="fits", output_fov=None):
        """
        Attributes
        ----------
        ra, dec = position in degrees
        fov = field of view
        output_size = output (display) image size in pixels (default = size).
                      output_size has no effect for fits format images.
        filters = string with filters to include
        format = data format (options are "jpg", "png" or "fits")
        """
        self.ra = ra
        self.dec = dec
        self.fov = fov
        self.size = int((self.fov/PIXEL_SCALE).value)
        self.color = color
        self.filters= filters
        self.format = format
        self.output_fov = output_fov
        if self.output_fov is not None:
            self.output_size = self.output_fov/PIXEL_SCALE
        else:
            self.output_size = self.size
        msg = 'The PanSTARRS 3PI survey images cover the sky north of dec-30 deg.'
        if self.dec<-30:
            raise ValueError(msg)

    def get_images(self):
        """Query ps1filenames.py service to get a list of images
        Returns a table with the results
        """

        service = "https://ps1images.stsci.edu/cgi-bin/ps1filenames.py"
        url = (f"{service}?ra={self.ra}&dec={self.dec}&size={self.size}"
               f"&format={self.format}&filters={self.filters}")
        table = Table.read(url, format='ascii')
        return table

    def get_url(self):
        """Get URL for images in the table
        Returns a string with the URL
        """
        if self.color and self.format == "fits":
            raise ValueError("color images are available only for jpg or png formats")
        if self.format not in ("jpg","png","fits"):
            raise ValueError("format must be one of jpg, png, fits")
        table = self.get_images()
        url = ("https://ps1images.stsci.edu/cgi-bin/fitscut.cgi?"
               f"ra={self.ra}&dec={self.dec}&size={self.size}&format={self.format}")
        if self.output_size:
            url = url + f"&output_size={self.output_size}"
        # sort filters from red to blue
        flist = ["yzirg".find(x) for x in table['filter']]
        table = table[np.argsort(flist)]
        if self.color:
            if len(table) > 3:
                # pick 3 filters
                table = table[[0,len(table)//2,len(table)-1]]
            for i, param in enumerate(["red","green","blue"]):
                url = url + f"&{param}={table['filename'][i]}"
        else:
            urlbase = url + "&red="
            url = []
            for filename in table['filename']:
                url.append(urlbase+filename)
        assert len(url)>0, 'empty url'
        return url

    def get_color_img(self):
        """Get color image at a sky position
        Returns the image
        """
        if self.format not in ("jpg","png"):
            raise ValueError("format must be jpg or png")
        url = self.get_url()
        if not self.color and isinstance(url,list):
            r = requests.get(url[0])
        else:
            r = requests.get(url)
        im = Image.open(BytesIO(r.content))
        return im

    def get_gray_img(self, filter="g"):
        """Get grayscale image at a sky position
        Returns the image
        """
        self.color = False
        if self.format not in ("jpg","png"):
            raise ValueError("format must be jpg or png")
        if filter not in list("grizy"):
            raise ValueError("filter must be one of grizy")
        url = self.get_url()
        r = requests.get(url[0])
        im = Image.open(BytesIO(r.content))
        return im

    def get_fits(self, filter=None, verbose=True):
        """Get fits image at a sky position
        Returns the image and header
        """
        self.format = "fits"
        self.filters = filter if filter is not None else "g"
        fitsurl = self.get_url()
        hl = fits.open(fitsurl[0])
        if verbose:
            print(hl.info())
        img, hdr = hl[0].data, hl[0].header
        return img, hdr

if __name__=='__main__':
    from astropy.coordinates import SkyCoord
    from astropy.visualization import ZScaleInterval
    from astropy.wcs import WCS

    tcoord = SkyCoord.from_name('TRAPPIST-1')
    ps = Panstarrs(ra=tcoord.ra.deg, dec=tcoord.dec.deg,
                   filters="g", format="fits", color=False)
    interval = ZScaleInterval(contrast=0.5)
    img, hdr = ps.get_fits(verbose=False)
    zmin, zmax = interval.get_limits(img)

    fig = pl.figure(figsize=(5,5))
    ax = fig.add_subplot(111, projection=WCS(hdr))
    ax.imshow(img, vmin=zmin, vmax=zmax, cmap='viridis', origin='lower')
    ax.grid(color='white', ls='solid')
    pl.show()

reason: import numpy,from astropy


text:
import torch
import gym
import numpy as np
import argparse
import matplotlib.pyplot as plt
from agent import Agent, Policy
from cp_cont import CartPoleEnv
from parallel_env import ParallelEnvs
import pandas as pd

# Policy training function
def train(env_name="ContinuousCartPole-v0"):
    update_steps = 10

    # Create a Gym environment
    # This creates 64 parallel envs running in 8 processes (8 threads each)
    env = ParallelEnvs(env_name, processes=5, envs_per_process=5)

    # Get dimensionalities of actions and observations
    action_space_dim = env.action_space.shape[-1]
    observation_space_dim = env.observation_space.shape[-1]

    # Instantiate agent and its policy
    policy = Policy(observation_space_dim, action_space_dim)
    agent = Agent(policy)
    # Arrays to keep track of rewards
    reward_history, timestep_history = [], []
    average_reward_history = []
    total_steps = 0

    # Run actual training
    reward_sum, timesteps = 0, 0
    done = False
    # Reset the environment and observe the initial state
    observation = env.reset()

    # Loop forever
    while True:
        # Get action from the agent
        action, action_probabilities = agent.get_action(observation)
        previous_observation = observation

        # Perform the action on the environment, get new state and reward
        observation, reward, done, info = env.step(action.detach().numpy())

        # Store action's outcome (so that the agent can improve its policy)
        for i in range(len(info["infos"])):
            env_done = False
            # Check if the environment is finished; if so, store cumulative reward
            for envid, envreward in info["finished"]:
                if envid == i:
                    reward_history.append(envreward)
                    average_reward_history.append(np.mean(reward_history[-500:]))
                    env_done = True
                    break
            agent.store_outcome(previous_observation[i], observation[i],
                                action_probabilities[i], reward[i], env_done)

        # Draw the frame, if desired
        #if args.render_training:
        #    env.render()

        # Store total episode reward
        timesteps += 1
        total_steps += 1

        if total_steps % update_steps == update_steps-1:
            agent.episode_finished(0)

        plot_freq = 5000
        if total_steps % plot_freq == plot_freq-1:
            # Training is finished - plot rewards
            plt.plot(reward_history)
            plt.plot(average_reward_history)
            plt.legend(["Reward", "500-episode average"])
            plt.title("Non-episodic AC reward history")
            plt.savefig("rewards_%s.png" % env_name)
            plt.clf()
            torch.save(agent.policy.state_dict(), "model.mdl")
            print("%d: Plot and model saved." % total_steps)
    return 1

# Function to test a trained policy
def test(episodes, agent, env):
    test_reward, test_len = 0, 0
    for ep in range(episodes):
        done = False
        observation = env.reset()
        while not done:
            # Similar to the training loop above -
            # get the action, act on the environment, save total reward
            # (evaluation=True makes the agent always return what it thinks to be
            # the best action - there is no exploration at this point)
            action, _ = agent.get_action(observation, evaluation=True)
            action = action.detach().numpy()
            observation, reward, done, info = env.step(action)

            if args.render_test:
                env.render()
            test_reward += reward
            test_len += 1
    print("Average test reward:", test_reward/episodes, "episode length:", test_len/episodes)

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("--test", "-t", type=str, default=None, help="Model to be tested")
    parser.add_argument("--env", type=str, default="ContinuousCartPole-v0", help="Environment to use")
    parser.add_argument("--render_test", action='store_true', help="Render test")
    args = parser.parse_args()

    # Print some stuff
    print("Environment:", args.env)
    #print("Training device:", agent.train_device)
    #print("Observation space dimensions:", observation_space_dim)
    #print("Action space dimensions:", action_space_dim)

    # If no model was passed, train a policy from scratch.
    # Otherwise load the policy from the file and go directly to testing.
    if args.test is None:
        try:
            train(env_name=args.env)
        # Handle Ctrl+C - save model and go to tests
        except KeyboardInterrupt:
            print("Interrupted!")
        model_file = "%s_params.mdl" % args.env
        torch.save(policy.state_dict(), model_file)
        print("Model saved to", model_file)
    else:
        # Instantiate agent and its policy
        env = gym.make(args.env)

        # Get dimensionalities of actions and observations
        action_space_dim = env.action_space.shape[-1]
        observation_space_dim = env.observation_space.shape[-1]
        policy = Policy(observation_space_dim, action_space_dim)
        agent = Agent(policy)
        state_dict = torch.load(args.test)
        policy.load_state_dict(state_dict)
        print("Testing...")
        test(100, agent, env)

max_stars_repo_name: Whaith/reinforcement-learning-project

max_issues_repo_name: Whaith/reinforcement-learning-project

max_forks_repo_name: Whaith/reinforcement-learning-project