text:
"""
Variational Autoencoder (VAE) implementation
"""
# MIT License
#
# Copyright (c) 2019 David Sandberg
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import tensorflow as tf
import numpy as np
import math

class VAE(object):

    def __init__(self, input_dims, learning_rate, warmup_temp_ph, nrof_stochastic_units, nrof_mlp_units):

        assert len(nrof_stochastic_units)==len(nrof_mlp_units)
        self.nrof_stochastic_units = nrof_stochastic_units
        self.nrof_mlp_units = nrof_mlp_units

        self.input_dims = input_dims

        self.learning_rate = learning_rate
        self.warmup_temp_ph = warmup_temp_ph

    def softlimit(self, x, limit=0.1):
        return tf.math.log(tf.math.exp(x) + 1.0 + limit)

    def dense(self, x, nrof_units, activation=None, training=True, use_batch_norm=True):
        x = tf.compat.v1.layers.Flatten()(x)
        x = tf.compat.v1.layers.Dense(units=nrof_units)(x)
        if use_batch_norm:
            x = tf.compat.v1.layers.BatchNormalization()(x, training=training)
        x = x if activation is None else activation(x)
        return x

    def mlp(self, x, nrof_units, activation, nrof_layers=1, training=True):
        for _ in range(nrof_layers):
            x = self.dense(x, nrof_units=nrof_units, activation=activation, training=training)
        return x

    def sample(self, mu, sigma):
        epsilon = tf.random.normal(tf.shape(mu), mean=0.0, stddev=1.0)
        return mu + sigma*epsilon

    def decoder_sample(self, z, mu, sigma, draw_dec_sample=False):
        if draw_dec_sample:
            return self.sample(mu, sigma)
        return z

    def log_normal2(self, x, mean, log_var, eps=1e-5):
        c = - 0.5 * math.log(2*math.pi)
        return c - log_var/2 - (x - mean)**2 / (2 * tf.math.exp(log_var) + eps)

    def build_graph(self, x, is_training):

        o = dict()
        dbg = dict()

        reuse = None if is_training else True

        with tf.compat.v1.variable_scope('model', reuse=reuse):

            nrof_layers = len(self.nrof_mlp_units)

            p_mu, p_sigma = [None]*nrof_layers, [None]*nrof_layers
            q_mu, q_sigma = [None]*nrof_layers, [None]*nrof_layers
            z = [None]*nrof_layers

            # Create encoder
            h = x
            for l in range(nrof_layers):
                dbg['enc_prevh_%d' % l ] = h
                h = self.mlp(h, self.nrof_mlp_units[l], activation=tf.nn.leaky_relu, nrof_layers=nrof_layers, training=is_training)
                dbg['h_%d' % l] = h

                """ prepare for bidirectional inference """
                q_mu[l] = self.dense(h, self.nrof_stochastic_units[l], training=is_training)
                q_sigma[l] = self.dense(h, self.nrof_stochastic_units[l], activation=self.softlimit, training=is_training)
                z[l] = self.sample(q_mu[l], q_sigma[l])
                h = z[l]

                dbg['z_%d' % l] = z[l]
                dbg['q_mu_%d' % l] = q_mu[l]
                dbg['q_sigma_%d' % l] = q_sigma[l]

            # Create decoder
            for l in range(nrof_layers-1, -1, -1):

                if l == nrof_layers-1:
                    """ At the highest latent layer, mu & sigma are identical to those output from encoder.
                        And making actual z is not necessary for the highest layer."""
                    mu, sigma = q_mu[l], q_sigma[l]
                    p_mu[l], p_sigma[l] = tf.zeros_like(mu), tf.ones_like(sigma)
                else:
                    """ prior is developed from z of the above layer """
                    p_mu[l] = self.dense(h, self.nrof_stochastic_units[l], training=is_training)
                    p_sigma[l] = self.dense(h, self.nrof_stochastic_units[l], activation=self.softlimit, training=is_training)

                dbg['p_mu_%d' % l] = p_mu[l]
                dbg['p_sigma_%d' % l] = p_sigma[l]
                dbg['z_%d' % l]  = z[l]

                h = self.mlp(z[l], self.nrof_mlp_units[l], activation=tf.nn.leaky_relu, nrof_layers=nrof_layers, training=is_training)
                dbg['h_%d' % l]  = h

            nrof_features = np.prod(self.input_dims)
            x_reconst = self.dense(h, nrof_features, activation=tf.nn.sigmoid, training=is_training, use_batch_norm=False) # Reconstruction

            if len(self.input_dims)==3:
                x_reconst = tf.reshape(x_reconst, (-1, self.input_dims[0], self.input_dims[1], self.input_dims[2]))

            dbg['x_reconst']  = x_reconst
            dbg['x_orig']  = x

            # Calculate reconstruction loss
            eps = 1e-5
            x_reconst_clip = tf.clip_by_value(x_reconst, eps, 1-eps)
            log_pxz = -tf.reduce_sum(tf.keras.losses.binary_crossentropy(x, x_reconst_clip), axis=[1,2])
            o['log_px'] = tf.reduce_mean(log_pxz)

            # Calculate ELBO, i.e. log P(X|Z) + temp*(log(P(Z|X)) - log(q(Z)))
            log_qz_list = []
            log_pz_list = []
            for l in range(nrof_layers):
                log_qz_l = self.log_normal2(z[l], q_mu[l], tf.math.log(q_sigma[l])*2)
                log_pz_l = self.log_normal2(z[l], p_mu[l], tf.math.log(p_sigma[l])*2)
                dbg['log_pz_l_%d' % l] = log_pz_l
                dbg['log_qz_l_%d' % l] = log_qz_l
                log_pz_list += [ tf.reduce_sum(log_pz_l, axis=1) ]  # Sum over log probabilities
                log_qz_list += [ tf.reduce_sum(log_qz_l, axis=1) ]  # Sum over log probabilities
            log_pz = tf.add_n(log_pz_list)
            log_qz = tf.add_n(log_qz_list)
            elbo = log_pxz + self.warmup_temp_ph * (log_pz - log_qz)

            o['log_pz'] = [ tf.reduce_mean(w) for w in log_pz_list ]  # Average over the batch
            o['log_qz'] = [ tf.reduce_mean(w) for w in log_qz_list ]  # Average over the batch

            o['kl_tot'] = -tf.reduce_mean(log_pz - log_qz)  # Average over the batch
            o['kl'] = [ -tf.reduce_mean(lpz-lqz) for (lpz, lqz) in zip(log_pz_list, log_qz_list) ]

            """ set losses """
            o['elbo'] = tf.reduce_mean(elbo)  # Average over the batch

            loss = -o['elbo']

            train_op = None
            if is_training:
                optimizer = tf.compat.v1.train.AdamOptimizer(learning_rate=self.learning_rate, beta1=0.5)
                grads = optimizer.compute_gradients(loss)  # Minimize -elbo
                for i,(g,v) in enumerate(grads):
                    if g is not None:
                        grads[i] = (tf.clip_by_norm(g,5),v) # clip gradients
                update_ops = tf.compat.v1.get_collection(tf.compat.v1.GraphKeys.UPDATE_OPS)
                with tf.control_dependencies(update_ops):
                    train_op = optimizer.apply_gradients(grads)

            return train_op, o, dbg

max_stars_repo_name: davidsandberg/LadderVAE

max_issues_repo_name: davidsandberg/LadderVAE

max_forks_repo_name: davidsandberg/LadderVAE


text:
"""monte_carlo_simulations Script for running monte carlo simulations

Script used for running monte carlo simulations.
"""
import numpy as np
import scipy
from stonesoup.types.state import GaussianState
from matplotlib import pyplot as plt

from utils import open_object
from utils import calc_metrics

from utils.scenario_generator import generate_scenario_2
from trackers.kalman_filter_view_AIS_as_measurement import kalman_filter_ais_as_measurement
from trackers.kalman_filter_independent_fusion import kalman_filter_independent_fusion
from trackers.kalman_filter_dependent_fusion import kalman_filter_dependent_fusion

from utils.save_figures import save_figure
from utils.latex_utils import populate_latex_table

# seeds
start_seed = 0
end_seed = 100  # normally 500
num_mc_iterations = end_seed - start_seed

# params
save_fig = True
print_latex = True

# scenario parameters
sigma_process_list = [0.05, 0.05, 0.05, 0.5, 0.5, 0.5, 3, 3, 3]  # [0.05, 0.05]  #
sigma_meas_radar_list = [1, 5, 20, 1, 5, 20, 1, 5, 20]  # [5, 20]  #
sigma_meas_ais_list = sigma_meas_radar_list
num_steps = 200
num_steps_metrics = 150

# dicts to store final results for printing in a latex friendly way
aanees_overall = {}
rmse_overall = {}

for sigma_process, sigma_meas_radar, sigma_meas_ais in zip(sigma_process_list, sigma_meas_radar_list,
                                                           sigma_meas_ais_list):
    stats_overall = []
    for seed in range(start_seed, end_seed):
        # generate scenario
        generate_scenario_2(seed=seed, permanent_save=False, sigma_process=sigma_process,
                            sigma_meas_radar=sigma_meas_radar, sigma_meas_ais=sigma_meas_ais, timesteps=num_steps)

        folder = "temp"  # temp instead of seed, as it is not a permanent save

        # load ground truth and the measurements
        data_folder = "../scenarios/scenario2/" + folder + "/"
        ground_truth = open_object.open_object(data_folder + "ground_truth.pk1")
        measurements_radar = open_object.open_object(data_folder + "measurements_radar.pk1")
        measurements_ais = open_object.open_object(data_folder + "measurements_ais.pk1")

        # load start_time
        start_time = open_object.open_object(data_folder + "start_time.pk1")

        # prior
        initial_covar = np.diag([sigma_meas_radar * sigma_meas_ais, sigma_meas_radar*sigma_process,
                                 sigma_meas_radar * sigma_meas_ais, sigma_meas_radar*sigma_process]) ** 2
        prior = GaussianState([1, 1.1, -1, 0.9], initial_covar, timestamp=start_time)

        # trackers
        kf_ais_as_measurement = kalman_filter_ais_as_measurement(measurements_radar, measurements_ais, start_time,
                                                                 prior,
                                                                 sigma_process=sigma_process,
                                                                 sigma_meas_radar=sigma_meas_radar,
                                                                 sigma_meas_ais=sigma_meas_ais)
        kf_independent_fusion = kalman_filter_independent_fusion(start_time,
                                                                 prior,
                                                                 sigma_process_radar=sigma_process,
                                                                 sigma_process_ais=sigma_process,
                                                                 sigma_meas_radar=sigma_meas_radar,
                                                                 sigma_meas_ais=sigma_meas_ais)
        kf_dependent_fusion = kalman_filter_dependent_fusion(measurements_radar, measurements_ais, start_time, prior,
                                                             sigma_process_radar=sigma_process,
                                                             sigma_process_ais=sigma_process,
                                                             sigma_meas_radar=sigma_meas_radar,
                                                             sigma_meas_ais=sigma_meas_ais)

        tracks_fused_independent, _, _ = kf_independent_fusion.track(measurements_radar, measurements_ais)
        tracks_fused_dependent, _, _ = kf_dependent_fusion.track()
        tracks_fused_ais_as_measurement, _ = kf_ais_as_measurement.track()

        # fix the length of the fusions to dependent (as it is a bit shorter)
        num_tracks = len(tracks_fused_dependent)
        # get the num_tracks last elements
        tracks_fused_independent = tracks_fused_independent[-num_tracks:]
        tracks_fused_ais_as_measurement = tracks_fused_ais_as_measurement[-num_tracks:]
        # remove the first element of ground_truth (because we only fuse the n-1 last with dependent fusion)
        ground_truth = ground_truth[-num_tracks:]

        # Calculate some metrics
        stats_individual = {}
        # calculate NEES
        stats_individual["NEES"] = calc_metrics.calc_nees(tracks_fused_independent[-num_steps_metrics:],
                                                          ground_truth[-num_steps_metrics:])
        # calculate ANEES
        stats_individual["ANEES"] = calc_metrics.calc_anees(stats_individual["NEES"])
        # calculate RMSE
        stats_individual["RMSE"] = calc_metrics.calc_rmse(tracks_fused_independent[-num_steps_metrics:],
                                                          ground_truth[-num_steps_metrics:])
        stats_individual["seed"] = seed
        stats_individual["fusion_type"] = "independent"
        stats_overall.append(stats_individual)

        stats_individual = {}
        # calculate NEES
        stats_individual["NEES"] = calc_metrics.calc_nees(tracks_fused_dependent[-num_steps_metrics:],
                                                          ground_truth[-num_steps_metrics:])
        # calculate ANEES
        stats_individual["ANEES"] = calc_metrics.calc_anees(stats_individual["NEES"])
        # calculate RMSE
        stats_individual["RMSE"] = calc_metrics.calc_rmse(tracks_fused_dependent[-num_steps_metrics:],
                                                          ground_truth[-num_steps_metrics:])
        stats_individual["seed"] = seed
        stats_individual["fusion_type"] = "dependent"
        stats_overall.append(stats_individual)

        stats_individual = {}
        # calculate NEES
        stats_individual["NEES"] = calc_metrics.calc_nees(tracks_fused_ais_as_measurement[-num_steps_metrics:],
                                                          ground_truth[-num_steps_metrics:])
        # calculate ANEES
        stats_individual["ANEES"] = calc_metrics.calc_anees(stats_individual["NEES"])
        # calculate RMSE
        stats_individual["RMSE"] = calc_metrics.calc_rmse(tracks_fused_ais_as_measurement[-num_steps_metrics:],
                                                          ground_truth[-num_steps_metrics:])
        stats_individual["seed"] = seed
        stats_individual["fusion_type"] = "ais as measurement"
        stats_overall.append(stats_individual)

    # plot some results
    num_tracks = num_steps_metrics
    alpha = 0.95
    ci_nees = scipy.stats.chi2.interval(alpha, 4)
    ci_anees = np.array(scipy.stats.chi2.interval(alpha, 4 * num_tracks)) / num_tracks

    # plot ANEES and confidence interval
    fig_ci_anees = plt.figure(figsize=(10, 6))
    ax_ci_anees = fig_ci_anees.add_subplot(1, 1, 1)
    ax_ci_anees.set_xlabel("MC seed")
    ax_ci_anees.set_ylabel("ANEES")

    # plot upper and lower confidence intervals
    ax_ci_anees.plot([start_seed, end_seed], [ci_anees[0], ci_anees[0]], color='red')
    ax_ci_anees.plot([start_seed, end_seed], [ci_anees[1], ci_anees[1]], color='red')

    anees_independent = [stat['ANEES'] for stat in stats_overall if stat["fusion_type"] == "independent"]
    anees_dependent = [stat['ANEES'] for stat in stats_overall if stat["fusion_type"] == "dependent"]
    anees_ais_as_measurement = [stat['ANEES'] for stat in stats_overall if stat["fusion_type"] == "ais as measurement"]

    rmse_independent = np.mean([stat['RMSE'] for stat in stats_overall if stat["fusion_type"] == "independent"])
    rmse_dependent = np.mean([stat['RMSE'] for stat in stats_overall if stat["fusion_type"] == "dependent"])
    rmse_ais_as_measurement = np.mean([stat['RMSE'] for stat in stats_overall if stat["fusion_type"] == "ais as "
                                                                                "measurement"])

    # plot the ANEES values
    ax_ci_anees.plot(list(range(start_seed, end_seed)), anees_independent, marker='+', ls='None', color='blue',
                     label='Independent')
    ax_ci_anees.plot(list(range(start_seed, end_seed)), anees_dependent, marker='+', ls='None', color='red',
                     label='Dependent')
    ax_ci_anees.plot(list(range(start_seed, end_seed)), anees_ais_as_measurement, marker='+', ls='None', color='green',
                     label='AIS as measurement')

    ax_ci_anees.legend()
    # fig_ci_anees.show()

    # save_figure("../results/scenario2/1996", "monte_carlo_all_trackers_same_params.svg", fig_ci_anees)

    # plot to show how the average ANEES approaches some limit. Plot a changing confidence interval?
    fig_ci_average_anees = plt.figure(figsize=(10, 6))
    ax_ci_average_anees = fig_ci_average_anees.add_subplot(1, 1, 1)
    ax_ci_average_anees.set_xlabel("MC iterations")
    ax_ci_average_anees.set_ylabel("ANEES")
    ax_ci_average_anees.set_ylim(2.5, 6.5)

    ci_average_anees = [
        np.array(scipy.stats.chi2.interval(alpha, 4 * num_tracks * num_MC_it)) / (num_tracks * num_MC_it)
        for num_MC_it in range(1, end_seed - start_seed + 1)]
    ci_average_anees_lower = [ci[0] for ci in ci_average_anees]
    ci_average_anees_upper = [ci[1] for ci in ci_average_anees]

    # calculate average anees
    # not the quickest, but works
    average_anees_independent_list = [np.average(anees_independent[:idx]) for idx in range(1, num_mc_iterations + 1)]
    average_anees_dependent_list = [np.average(anees_dependent[:idx]) for idx in range(1, num_mc_iterations + 1)]
    average_anees_ais_as_measurement_list = [np.average(anees_ais_as_measurement[:idx])
                                             for idx in range(1, num_mc_iterations + 1)]

    # store anees and rmse in dicts
    aanees_overall[(sigma_process, sigma_meas_radar, "indep")] = average_anees_independent_list[-1]
    aanees_overall[(sigma_process, sigma_meas_radar, "dep")] = average_anees_dependent_list[-1]
    aanees_overall[(sigma_process, sigma_meas_radar, "kf")] = average_anees_ais_as_measurement_list[-1]
    rmse_overall[(sigma_process, sigma_meas_radar, "indep")] = rmse_independent
    rmse_overall[(sigma_process, sigma_meas_radar, "dep")] = rmse_dependent
    rmse_overall[(sigma_process, sigma_meas_radar, "kf")] = rmse_ais_as_measurement

    # plot upper and lower confidence intervals
    ax_ci_average_anees.plot(list(range(start_seed, end_seed)), ci_average_anees_lower, color='black',
                             label='Confidence Intervals')
    ax_ci_average_anees.plot(list(range(start_seed, end_seed)), ci_average_anees_upper, color='black')

    # plot the average ANEES
    ax_ci_average_anees.plot(list(range(start_seed, end_seed)), average_anees_independent_list, color='blue',
                             label='Independent')
    ax_ci_average_anees.plot(list(range(start_seed, end_seed)), average_anees_dependent_list, color='red',
                             label='Dependent')
    ax_ci_average_anees.plot(list(range(start_seed, end_seed)), average_anees_ais_as_measurement_list, color='green',
                             label='KF fusion')

    # set title
    title = "ANEES with parameters $\sigma_{process}=$" + str(sigma_process) + ",  $\sigma_{AIS}=$" + \
            str(sigma_meas_ais) + ", $\sigma_{Radar}=$" + str(sigma_meas_radar)
    ax_ci_average_anees.set_title(title, fontsize=20)

    ax_ci_average_anees.legend(prop={'size': 12})
    fig_ci_average_anees.show()

    if save_fig:
        fig_name = "process_" + str(sigma_process) + "_AIS_" + str(sigma_meas_ais) + "_Radar_" \
                   + str(sigma_meas_radar) + ".pdf"
        save_figure("../results/final_results/scenario1/mc_average_anees", fig_name, fig_ci_average_anees)

    # print some results
    print("Process: " + str(sigma_process) + " AIS: " + str(sigma_meas_ais) + " Radar: " + str(sigma_meas_radar))
    print("Average ANEES independent: " + str(average_anees_independent_list[-1]))
    print("Average ANEES dependent: " + str(average_anees_dependent_list[-1]))
    print("Average ANEES AIS as measurement: " + str(average_anees_ais_as_measurement_list[-1]))
    print("RMSE independent: " + str(rmse_independent))
    print("RMSE dependent: " + str(rmse_dependent))
    print("RMSE AIS as measurement: " + str(rmse_ais_as_measurement))
    print("")

if print_latex:
    fusion_types = ["indep", "dep", "kf"]
    for fusion_type in fusion_types:
        print("AANEES")
        print(fusion_type + "\n" + populate_latex_table(aanees_overall, fusion_type))
        print("")
        print("RMSE")
        print(fusion_type + "\n" + populate_latex_table(rmse_overall, fusion_type))
        print("")

print("CI intervals: " + str(ci_average_anees_lower[-1]) + ", " + str(ci_average_anees_upper[-1]))

ext: py

max_stars_repo_name: jonassagild/Track-to-Track-Fusion

max_issues_repo_name: jonassagild/Track-to-Track-Fusion

max_forks_repo_name: jonassagild/Track-to-Track-Fusion

reason: import numpy,import scipy


text:
#!/usr/bin/env python3
# Copyright (c) 2014-2021 Megvii Inc. All rights reserved.
import megengine as mge
import megengine.module as M
import numpy as np
import pytest

from basecls.layers import SE, DropPath, conv2d, gap2d, linear, norm2d, pool2d

@pytest.mark.parametrize("w_in", [4])
@pytest.mark.parametrize("w_out", [8])
@pytest.mark.parametrize("k", [3, 5])
@pytest.mark.parametrize("stride", [1, 2])
@pytest.mark.parametrize("dilation", [1, 2])
@pytest.mark.parametrize("groups", [1, 2])
@pytest.mark.parametrize("bias", [True, False])
def test_conv2d(w_in, w_out, k, stride, dilation, groups, bias):
    m = conv2d(w_in, w_out, k, stride=stride, dilation=dilation, groups=groups, bias=bias)
    assert isinstance(m, M.Conv2d)

    m(mge.random.normal(size=(2, 4, 8, 8)))

@pytest.mark.parametrize("drop_prob", [0.0, 0.5])
def test_drop_path(drop_prob):
    m = DropPath(drop_prob)
    assert isinstance(m, M.Module)

    m.training = True
    m(mge.random.normal(size=(2, 4, 8, 8)))

    m.training = False
    x = np.random.rand(2, 4, 8, 8).astype("float32")
    y = m(mge.Tensor(x)).numpy()
    np.testing.assert_allclose(y, x, rtol=1e-4, atol=1e-6)

@pytest.mark.parametrize("shape", [1, (7, 7)])
def test_gap2d(shape):
    m = gap2d(shape)
    assert isinstance(m, M.AdaptiveAvgPool2d)

    m(mge.random.normal(size=(2, 4, 8, 8)))

@pytest.mark.parametrize("w_in", [4])
@pytest.mark.parametrize("w_out", [8])
@pytest.mark.parametrize("bias", [True, False])
def test_linear(w_in, w_out, bias):
    m = linear(w_in, w_out, bias=bias)
    assert isinstance(m, M.Linear)

    m(mge.random.normal(size=(2, 8, 4)))

# TODO: "GN", "IN" and "LN" need different hyper-parameters
@pytest.mark.parametrize("name", [None, "BN", "SyncBN"])
@pytest.mark.parametrize("w_in", [4])
def test_norm2d(name, w_in):
    m = norm2d(name, w_in)
    assert isinstance(m, M.Module)

    m(mge.random.normal(size=(2, 4, 8, 8)))

@pytest.mark.parametrize("k", [3, 5])
@pytest.mark.parametrize("stride", [1, 2])
@pytest.mark.parametrize("name", ["avg", "max"])
def test_pool2d(k, stride, name):
    m = pool2d(k, stride=stride, name=name)
    assert isinstance(m, M.Module)

    m(mge.random.normal(size=(2, 4, 8, 8)))

@pytest.mark.parametrize("w_in", [8])
@pytest.mark.parametrize("w_se", [4])
@pytest.mark.parametrize("act_name", ["relu", "silu"])
def test_se(w_in, w_se, act_name):
    m = SE(w_in, w_se, act_name)
    assert isinstance(m, M.Module)

    m(mge.random.normal(size=(2, 8, 8, 8)))

max_stars_repo_name: megvii-research/basecls

max_issues_repo_name: megvii-research/basecls

max_forks_repo_name: megvii-research/basecls

reason: import numpy


text:
import oneflow as flow
import numpy as np
import oneflow.typing as tp
from typing import Tuple
import oneflow.nn as nn
import math
import time
import datetime

class KnowledgeConsistentAttention(object):
    def __init__(self, patch_size=3, propagate_size=3, stride=1,ti=""):
        super(KnowledgeConsistentAttention, self).__init__()
        self.patch_size = patch_size
        self.propagate_size = propagate_size
        self.stride = stride
        self.prop_kernels = None
        self.att_scores_prev = None
        self.masks_prev = None
        self.time =ti
        # self.ratio = nn.Parameter(torch.ones(1))

    def buildnet(self, foreground, masks):
        #time = datetime.datetime.now().strftime('%Y-%m-%d-%H_%M_%S_%f')
        bz, nc, w, h = foreground.shape
        if masks.shape[3] != foreground.shape[3]:
            # masks = F.interpolate(masks, foreground.size()[2:])
            if foreground.shape[3] >= masks.shape[3]:
                beishu = (foreground.shape[3]) / (masks.shape[3])
                masks = flow.layers.upsample_2d(masks, size=(beishu, beishu))
            else:
                kk = masks.shape[3] - foreground.shape[3]
                strid = int((masks.shape[3] - 4) / (masks.shape[3] - kk - 1))
                initializer = flow.truncated_normal(0.1)
                masks = flow.layers.conv2d(inputs=masks, filters=masks.shape[1], kernel_size=4,
                                           kernel_initializer=initializer, name="Conv2d" + self.time, padding='VALID',
                                           strides=strid)
        # background = foreground.clone()
        background = foreground
        background = background
        conv_kernels_all = flow.reshape(x=background, shape=(bz, nc, w * h, 1, 1))
        # conv_kernels_all = conv_kernels_all.permute(0, 2, 1, 3, 4)
        conv_kernels_all = flow.transpose(a=conv_kernels_all, perm=[0, 2, 1, 3, 4])
        output_tensor = []
        att_score = []
        for i in range(bz):
            time = datetime.datetime.now().strftime('%Y-%m-%d-%H_%M_%S_%f')
            # feature_map = foreground[i:i+1]
            feature_map = flow.slice(x=foreground, begin=[i, None, None, None], size=[1, None, None, None])
            # conv_kernels = conv_kernels_all[i] + 0.0000001
            conv_kernels = flow.slice(x=conv_kernels_all, begin=[i, None, None, None, None],
                                      size=[1, None, None, None, None])
            conv_kernels = flow.squeeze(input=conv_kernels, axis=[0])
            constant11 = flow.constant_like(like=conv_kernels, value=0.0000001, dtype=flow.float32)
            conv_kernels = flow.math.add(x=conv_kernels, y=constant11)
            # norm_factor = torch.sum(conv_kernels**2, [1, 2, 3], keepdim = True)**0.5
            constant1 = flow.constant_like(like=conv_kernels, value=2, dtype=flow.float32)
            conv_22 = flow.math.pow(x=conv_kernels, y=constant1)
            norm_factor = flow.math.reduce_sum(conv_22, (1, 2, 3), keepdims=True)
            constant = flow.constant_like(like=norm_factor, value=0.5, dtype=flow.float32)
            norm_factor = flow.math.pow(x=norm_factor, y=constant)
            conv_kernels = conv_kernels / norm_factor

            # conv_result = F.conv2d(feature_map, conv_kernels, padding = self.patch_size//2)
            conv_result = flow.nn.conv2d(feature_map, conv_kernels,
                                         padding=[0, 0, self.patch_size // 2, self.patch_size // 2], strides=1,name="nn"+self.time+str(i))
            if self.propagate_size != 1:
                if self.prop_kernels is None:
                    constant_blob = flow.constant(value=1.5, shape=(
                    conv_result.shape[1], 1, self.propagate_size, self.propagate_size), dtype=flow.float)
                    self.prop_kernels = flow.ones_like(like=constant_blob)
                    # self.prop_kernels = torch.ones([conv_result.size(1), 1, self.propagate_size, self.propagate_size])
                    # self.prop_kernels.requires_grad = False
                    # self.prop_kernels = self.prop_kernels.cuda()
                # conv_result = F.avg_pool2d(conv_result, 3, 1, padding = 1)*9
                conv_result = flow.nn.avg_pool2d(conv_result, ksize=3, strides=1, padding=[0, 0, 1, 1]) * 9  #
            attention_scores = flow.nn.softmax(conv_result, axis=1)
            if self.att_scores_prev is not None:
                self.ratio = flow.get_variable(shape=(1), name="ratio", initializer=init, trainable=True)
                # attention_scores = (self.att_scores_prev[i:i+1]*self.masks_prev[i:i+1] + attention_scores * (torch.abs(self.ratio)+1e-7))/(self.masks_prev[i:i+1]+(torch.abs(self.ratio)+1e-7))
                attention_scores = (self.att_scores_prev[i:i + 1] * self.masks_prev[i:i + 1] + attention_scores * (
                        flow.math.abs(self.ratio) + 1e-7)) / (
                                               self.masks_prev[i:i + 1] + (flow.math.abs(self.ratio) + 1e-7))

            att_score.append(attention_scores)
            # feature_map = F.conv_transpose2d(attention_scores, conv_kernels, stride = 1, padding = self.patch_size//2)
            # feature_map = F.conv_transpose2d(attention_scores, conv_kernels, stride=1, padding=self.patch_size // 2)

            weight = conv_kernels

            feature_map = flow.nn.torch_conv2d_transpose(value=attention_scores, filter=weight, padding_needed=(
            self.patch_size // 2 * 2, self.patch_size // 2 * 2), strides=1,
                                                         name="deconv" + self.time +str(i), output_padding=(0, 0))
            final_output = feature_map
            output_tensor.append(final_output)
        self.att_scores_prev = flow.concat(att_score, axis=0)
        self.att_scores_prev = flow.reshape(x=self.att_scores_prev, shape=[bz, h * w, h, w])
        # self.masks_prev = masks.view(bz, 1, h, w)
        self.masks_prev = flow.reshape(x=masks, shape=[bz, 1, h, w])
        # return torch.cat(output_tensor, dim = 0)
        return flow.concat(output_tensor, axis=0)

class AttentionModule(object):

    def __init__(self, inchannel, patch_size_list=[1], propagate_size_list=[3], stride_list=[1],tim=""):
        assert isinstance(patch_size_list,
                          list), "patch_size should be a list containing scales, or you should use Contextual Attention to initialize your module"
        assert len(patch_size_list) == len(propagate_size_list) and len(propagate_size_list) == len(
            stride_list), "the input_lists should have same lengths"
        super(AttentionModule, self).__init__()

        self.inchannel = inchannel
        self.patch_size_list = patch_size_list
        self.propagate_size_list = propagate_size_list
        self.stride_list = stride_list
        self.attt = KnowledgeConsistentAttention(self.patch_size_list[0], self.propagate_size_list[0],
                                                 self.stride_list[0], ti=tim)

        self.num_of_modules = len(patch_size_list)
        self.time =tim
    def buildnet(self, foreground, mask):
        #time = datetime.datetime.now().strftime('%Y-%m-%d-%H_%M_%S_%f')
        # self.attt = KnowledgeConsistentAttention(self.patch_size_list[0], self.propagate_size_list[0],self.stride_list[0])

        outputs = self.attt.buildnet(foreground, mask)
        outputs = flow.concat(inputs=[outputs, foreground], axis=1)
        # self.combiner = nn.Conv2d(inchannel * 2, inchannel, kernel_size=1)
        # outputs = self.combiner(outputs)
        initializer = flow.truncated_normal(0.1)
        outputs = flow.layers.conv2d(inputs=outputs, filters=self.inchannel, kernel_size=1, name="conv11" + self.time,
                                     kernel_initializer=initializer)
        return outputs

# @flow.global_function()
# def deconv2d_Job(x: tp.Numpy.Placeholder((2, 512, 6, 6)),y: tp.Numpy.Placeholder((2, 1, 6, 6))
# ) -> tp.Numpy:
# deconv = AttentionModule(512)
# output = deconv.buildnet(x,y)
# return output
# x = np.random.randn(2, 512, 6, 6).astype(np.float32)
# y = np.random.randn(2, 1, 6, 6).astype(np.float32)
# out = deconv2d_Job(x,y)
# print(out.shape)

# input =torch.randn(1, 512, 32, 32)
# filter=torch.randn(1, 1, 128, 128)
# deconv =AttentionModule(512)
# output =deconv(input,filter)
# print(output.size())


text:
# -*- coding: utf-8 -*-
"""Histogram-based Outlier Detection (HBOS)
"""
# Author: Yue Zhao <zhaoy@cmu.edu>
# License: BSD 2 clause

from __future__ import division
from __future__ import print_function

import numpy as np
from numba import njit
from sklearn.utils import check_array
from sklearn.utils.validation import check_is_fitted

from ..utils.utility import check_parameter
from ..utils.utility import invert_order

from .base import BaseDetector

class HBOS(BaseDetector):
    """Histogram- based outlier detection (HBOS) is an efficient unsupervised
    method. It assumes the feature independence and calculates the degree
    of outlyingness by building histograms. See :cite:`goldstein2012histogram`
    for details.

    Parameters
    ----------
    n_bins : int, optional (default=10)
        The number of bins.

    alpha : float in (0, 1), optional (default=0.1)
        The regularizer for preventing overflow.

    tol : float in (0, 1), optional (default=0.5)
        The parameter to decide the flexibility while dealing
        the samples falling outside the bins.

    contamination : float in (0., 0.5), optional (default=0.1)
        The amount of contamination of the data set,
        i.e. the proportion of outliers in the data set. Used when fitting to
        define the threshold on the decision function.

    Attributes
    ----------
    bin_edges_ : numpy array of shape (n_bins + 1, n_features )
        The edges of the bins.

    hist_ : numpy array of shape (n_bins, n_features)
        The density of each histogram.

    decision_scores_ : numpy array of shape (n_samples,)
        The outlier scores of the training data.
        The higher, the more abnormal. Outliers tend to have higher
        scores. This value is available once the detector is fitted.

    threshold_ : float
        The threshold is based on ``contamination``. It is the
        ``n_samples * contamination`` most abnormal samples in
        ``decision_scores_``. The threshold is calculated for generating
        binary outlier labels.

    labels_ : int, either 0 or 1
        The binary labels of the training data. 0 stands for inliers
        and 1 for outliers/anomalies. It is generated by applying
        ``threshold_`` on ``decision_scores_``.
    """

    def __init__(self, n_bins=10, alpha=0.1, tol=0.5, contamination=0.1):
        super(HBOS, self).__init__(contamination=contamination)
        self.n_bins = n_bins
        self.alpha = alpha
        self.tol = tol

        check_parameter(alpha, 0, 1, param_name='alpha')
        check_parameter(tol, 0, 1, param_name='tol')

    def fit(self, X, y=None):
        """Fit detector. y is ignored in unsupervised methods.

        Parameters
        ----------
        X : numpy array of shape (n_samples, n_features)
            The input samples.

        y : Ignored
            Not used, present for API consistency by convention.

        Returns
        -------
        self : object
            Fitted estimator.
        """
        # validate inputs X and y (optional)
        X = check_array(X)
        self._set_n_classes(y)

        n_samples, n_features = X.shape[0], X.shape[1]
        self.hist_ = np.zeros([self.n_bins, n_features])
        self.bin_edges_ = np.zeros([self.n_bins + 1, n_features])

        # build the histograms for all dimensions
        for i in range(n_features):
            self.hist_[:, i], self.bin_edges_[:, i] = \
                np.histogram(X[:, i], bins=self.n_bins, density=True)
            # the sum of (width * height) should equal to 1
            assert (np.isclose(1, np.sum(
                self.hist_[:, i] * np.diff(self.bin_edges_[:, i])), atol=0.1))

        # outlier_scores = self._calculate_outlier_scores(X)
        outlier_scores = _calculate_outlier_scores(X, self.bin_edges_,
                                                   self.hist_,
                                                   self.n_bins,
                                                   self.alpha, self.tol)

        # invert decision_scores_. Outliers comes with higher outlier scores
        self.decision_scores_ = invert_order(np.sum(outlier_scores, axis=1))
        self._process_decision_scores()
        return self

    def decision_function(self, X):
        """Predict raw anomaly score of X using the fitted detector.

        The anomaly score of an input sample is computed based on different
        detector algorithms. For consistency, outliers are assigned with
        larger anomaly scores.

        Parameters
        ----------
        X : numpy array of shape (n_samples, n_features)
            The training input samples. Sparse matrices are accepted only
            if they are supported by the base estimator.

        Returns
        -------
        anomaly_scores : numpy array of shape (n_samples,)
            The anomaly score of the input samples.
        """
        check_is_fitted(self, ['hist_', 'bin_edges_'])
        X = check_array(X)

        # outlier_scores = self._calculate_outlier_scores(X)
        outlier_scores = _calculate_outlier_scores(X, self.bin_edges_,
                                                   self.hist_,
                                                   self.n_bins,
                                                   self.alpha, self.tol)
        return invert_order(np.sum(outlier_scores, axis=1))

@njit
def _calculate_outlier_scores(X, bin_edges, hist, n_bins, alpha,
                              tol):  # pragma: no cover
    """The internal function to calculate the outlier scores based on
    the bins and histograms constructed with the training data. The program
    is optimized through numba. It is excluded from coverage test for
    eliminating the redundancy.

    Parameters
    ----------
    X : numpy array of shape (n_samples, n_features)
        The input samples.

    bin_edges : numpy array of shape (n_bins + 1, n_features )
        The edges of the bins.

    hist : numpy array of shape (n_bins, n_features)
        The density of each histogram.

    n_bins : int, optional (default=10)
        The number of bins.

    alpha : float in (0, 1), optional (default=0.1)
        The regularizer for preventing overflow.

    tol : float in (0, 1), optional (default=0.1)
        The parameter to decide the flexibility while dealing
        the samples falling outside the bins.

    Returns
    -------
    outlier_scores : numpy array of shape (n_samples, n_features)
        Outlier scores on all features (dimensions).
    """

    n_samples, n_features = X.shape[0], X.shape[1]
    outlier_scores = np.zeros(shape=(n_samples, n_features))

    for i in range(n_features):

        # Find the indices of the bins to which each value belongs.
        # See documentation for np.digitize since it is tricky
        # >>> x = np.array([0.2, 6.4, 3.0, 1.6, -1, 100, 10])
        # >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
        # >>> np.digitize(x, bins, right=True)
        # array([1, 4, 3, 2, 0, 5, 4], dtype=int64)

        bin_inds = np.digitize(X[:, i], bin_edges[:, i], right=True)

        # Calculate the outlying scores on dimension i
        # Add a regularizer for preventing overflow
        out_score_i = np.log2(hist[:, i] + alpha)

        for j in range(n_samples):

            # If the sample does not belong to any bins
            # bin_ind == 0 (fall outside since it is too small)
            if bin_inds[j] == 0:
                dist = bin_edges[0, i] - X[j, i]
                bin_width = bin_edges[1, i] - bin_edges[0, i]

                # If it is only slightly lower than the smallest bin edge
                # assign it to bin 1
                if dist <= bin_width * tol:
                    outlier_scores[j, i] = out_score_i[0]
                else:
                    outlier_scores[j, i] = np.min(out_score_i)

            # If the sample does not belong to any bins
            # bin_ind == n_bins+1 (fall outside since it is too large)
            elif bin_inds[j] == n_bins + 1:
                dist = X[j, i] - bin_edges[-1, i]
                bin_width = bin_edges[-1, i] - bin_edges[-2, i]

                # If it is only slightly larger than the largest bin edge
                # assign it to the last bin
                if dist <= bin_width * tol:
                    outlier_scores[j, i] = out_score_i[n_bins - 1]
                else:
                    outlier_scores[j, i] = np.min(out_score_i)
            else:
                outlier_scores[j, i] = out_score_i[bin_inds[j] - 1]

    return outlier_scores

max_issues_repo_name: durgeshsamariya/pyod

max_forks_repo_name: durgeshsamariya/pyod

reason: import numpy,from numba


text:
from mpl_toolkits.basemap import Basemap,shiftgrid
import matplotlib.pyplot as plt
import numpy as np
import glob, math
import pandas as pd
import warnings, matplotlib.cbook
warnings.filterwarnings("ignore", category=FutureWarning)
from rfsks_support.plotting_libs import shoot, equi, plot_topo, plot_events_loc, plot_merc
import os

# plot_h = 1
# plot_kappa = 1

def plot_h_kappa(h_k_file = "rfsks_support/h-kappa-values.txt",all_stationsfile = "results/InfoRF/all_stations_rf_retrieved.txt",plot_h = 1,plot_kappa = 1):
    fig_loc = all_stationsfile.split("/")[:-1]
    fig_loc_str = "/".join(fig_loc)

    ## read the stations inventory
    h_kappa_df = pd.read_csv(h_k_file,sep=",",header=None,names=['network','station', 'latitude', 'longitude','thickness','kappa'])
    if h_kappa_df.shape[0]:
        all_stations_df = pd.read_csv(all_stationsfile,sep="|")
        outfig_h = fig_loc_str+"/"+'all_stations_thickness_map.png'
        outfig_k = fig_loc_str+"/"+'all_stations_kappa_map.png'
        if plot_h and not os.path.exists(outfig_h):
            fig = plt.figure(figsize=(20,12))
            map = plot_merc(resolution='h',llcrnrlon=int(all_stations_df['Longitude'].min())-1, llcrnrlat=int(all_stations_df['Latitude'].min())-1,urcrnrlon=math.ceil(all_stations_df['Longitude'].max())+1, urcrnrlat=math.ceil(all_stations_df['Latitude'].max())+1,topo=False)

            stalon = h_kappa_df['longitude'].values
            stalat = h_kappa_df['latitude'].values
            x,y = map(stalon, stalat)
            sc = plt.scatter(x, y, c=h_kappa_df['thickness'].values, s=30,edgecolors='k',linewidth=0.1)
            cb = plt.colorbar(sc,fraction=0.025, pad=0.04)
            cb.ax.set_ylabel('Thickness in km',labelpad=15,fontsize=18)

            x1,y1 = map(all_stations_df['Longitude'].values, all_stations_df['Latitude'].values)
            plt.scatter(x1, y1, s=30, facecolors='none', edgecolors='k',linewidth=0.1)
            plt.tight_layout()

            # plt.legend(loc=4,fontsize=8)
            plt.savefig(outfig_h,dpi=300,bbox_inches='tight')
            plt.close('all')

        if plot_kappa and not os.path.exists(outfig_k):
            fig = plt.figure(figsize=(20,16))
            map = plot_merc(resolution='h',llcrnrlon=int(all_stations_df['Longitude'].min())-1, llcrnrlat=int(all_stations_df['Latitude'].min())-1,urcrnrlon=math.ceil(all_stations_df['Longitude'].max())+1, urcrnrlat=math.ceil(all_stations_df['Latitude'].max())+1,topo=False)

            stalon = h_kappa_df['longitude'].values
            stalat = h_kappa_df['latitude'].values
            x,y = map(stalon, stalat)
            sc = plt.scatter(x, y, c=h_kappa_df['kappa'].values, s=30,edgecolors='k',linewidth=0.1)
            cb = plt.colorbar(sc,fraction=0.025, pad=0.04)
            cb.set_label(r'$\frac{V_p}{V_s}$',labelpad=15,rotation=0,fontsize=18)

            x1,y1 = map(all_stations_df['Longitude'].values, all_stations_df['Latitude'].values)
            plt.scatter(x1, y1, s=30, facecolors='none', edgecolors='k',linewidth=0.1)
            plt.tight_layout()

            # plt.legend(loc=4,fontsize=8)
            plt.savefig(outfig_k,dpi=300,bbox_inches='tight')
            plt.close('all')
if __name__=="__main__":
    plot_h_kappa()

max_stars_repo_name: earthinversion/STADIUM-Py

max_issues_repo_name: cplegendre/STADIUM-Py

max_forks_repo_name: earthinversion/STADIUM-Py


text:
# ==========================================================================
#
#   Copyright NumFOCUS
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#          https://www.apache.org/licenses/LICENSE-2.0.txt
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.
#
# ==========================================================================*/

# also test the import callback feature

import sys
import os
import numpy as np
import pathlib

import itk

filename = sys.argv[1]

# xarray conversion
try:
    import xarray as xr

    print("Testing xarray conversion")

    image = itk.imread(filename)
    image.SetSpacing((0.1, 0.2))
    image.SetOrigin((30.0, 44.0))
    theta = np.radians(30)
    cosine = np.cos(theta)
    sine = np.sin(theta)
    rotation = np.array(((cosine, -sine), (sine, cosine)))
    image.SetDirection(rotation)
    image["MyMeta"] = 4.0

    data_array = itk.xarray_from_image(image)
    # Default name
    assert data_array.name == "image"
    image.SetObjectName("test_image")
    data_array = itk.xarray_from_image(image)
    assert data_array.name == "test_image"
    assert data_array.dims[0] == "y"
    assert data_array.dims[1] == "x"
    assert data_array.dims[2] == "c"
    assert np.array_equal(data_array.values, itk.array_from_image(image))
    assert len(data_array.coords["x"]) == 256
    assert len(data_array.coords["y"]) == 256
    assert len(data_array.coords["c"]) == 3
    assert data_array.coords["x"][0] == 30.0
    assert data_array.coords["x"][1] == 30.1
    assert data_array.coords["y"][0] == 44.0
    assert data_array.coords["y"][1] == 44.2
    assert data_array.coords["c"][0] == 0
    assert data_array.coords["c"][1] == 1
    assert data_array.attrs["direction"][0, 0] == cosine
    assert data_array.attrs["direction"][0, 1] == sine
    assert data_array.attrs["direction"][1, 0] == -sine
    assert data_array.attrs["direction"][1, 1] == cosine
    assert data_array.attrs["MyMeta"] == 4.0

    round_trip = itk.image_from_xarray(data_array)
    assert round_trip.GetObjectName() == "test_image"
    assert np.array_equal(itk.array_from_image(round_trip), itk.array_from_image(image))
    spacing = round_trip.GetSpacing()
    assert np.isclose(spacing[0], 0.1)
    assert np.isclose(spacing[1], 0.2)
    origin = round_trip.GetOrigin()
    assert np.isclose(origin[0], 30.0)
    assert np.isclose(origin[1], 44.0)
    direction = round_trip.GetDirection()
    assert np.isclose(direction(0, 0), cosine)
    assert np.isclose(direction(0, 1), -sine)
    assert np.isclose(direction(1, 0), sine)
    assert np.isclose(direction(1, 1), cosine)
    assert round_trip["MyMeta"] == 4.0

    wrong_order = data_array.swap_dims({"y": "z"})
    try:
        round_trip = itk.image_from_xarray(wrong_order)
        assert False
    except ValueError:
        pass

    # Check empty array
    empty_array = np.array([], dtype=np.uint8)
    empty_array.shape = (0, 0, 0)
    empty_image = itk.image_from_array(empty_array)
    empty_da = itk.xarray_from_image(empty_image)
    empty_image_round = itk.image_from_xarray(empty_da)

    # Check order
    arr = np.random.randint(0, 255, size=(4, 5, 6), dtype=np.uint8)
    data_array = xr.DataArray(arr, dims=["z", "y", "x"])
    image = itk.image_from_xarray(data_array)
    assert np.allclose(arr, itk.array_view_from_image(image))
    assert np.allclose(arr.shape, itk.array_view_from_image(image).shape)

    data_array = xr.DataArray(arr, dims=["x", "y", "z"])
    image = itk.image_from_xarray(data_array)
    assert np.allclose(arr.transpose(), itk.array_view_from_image(image))
    assert np.allclose(arr.shape[::-1], itk.array_view_from_image(image).shape)

    data_array = xr.DataArray(arr, dims=["y", "x", "c"])
    image = itk.image_from_xarray(data_array)
    assert np.allclose(arr, itk.array_view_from_image(image))
    assert np.allclose(arr.shape, itk.array_view_from_image(image).shape)

    data_array = xr.DataArray(arr, dims=["c", "x", "y"])
    image = itk.image_from_xarray(data_array)
    assert np.allclose(arr.transpose(), itk.array_view_from_image(image))
    assert np.allclose(arr.shape[::-1], itk.array_view_from_image(image).shape)

    # Test in-place "view" where xarray data is valid only as long as ITK image
    data_array = itk.xarray_from_image(image, view=True)
    assert data_array.name == image.GetObjectName()
    # verify we can run a computation on xarray data
    assert float(data_array.min().compute()) == 1.0

    data_array = xr.DataArray(arr, dims=["q", "x", "y"])
    try:
        image = itk.image_from_xarray(data_array)
        assert False
    except ValueError:
        pass

    if "(<itkCType unsigned char>, 4)" in itk.Image.GetTypesAsList():
        arr = np.random.randint(0, 255, size=(4, 5, 6, 3), dtype=np.uint8)
        data_array = xr.DataArray(arr, dims=["t", "z", "y", "x"])
        image = itk.image_from_xarray(data_array)
        assert np.allclose(arr, itk.array_view_from_image(image))
        assert np.allclose(arr.shape, itk.array_view_from_image(image).shape)

except ImportError:
    print("xarray not imported. Skipping xarray conversion tests")
    pass

max_stars_repo_name: khangthk/ITK

max_issues_repo_name: khangthk/ITK

max_forks_repo_name: khangthk/ITK


text:
# -*- coding: UTF-8 -*-
import numpy as np
import pandas as pd
import time
import json
import os
import requests

from apdt.config import _config

class DarkSkyAPI():
    '''DarkSkyAPI is a spider on darsky.net. Help you download weather history data for a set of sites using a set of accounts.

    darksky_secret_key_num and darksky_secret_key_0 should be decleared in config.py.

    Method
    ------
        __init__: Initialize.
            sites: a list of sites to collect data on. e.g. ['1001A', '1414A'].
            start_date: the start date (include) to collect data. e.g. '2019-01-01'.
            end_date: the end date (include) to collect data. e.g. '2019-12-31'.
            DailyAccountQuota: the maximum call quota for each accounts per day. e.g. (default) 990.
            save_path: the path to put download data. e.g. (default) 'data/weather/'.

        run: excute download.
            flag: set True to start download. default: False.
            shift: claim a start point for task, this can be useful for continue a job from previous start. default: 0.

    '''
    def __init__(self, sites, start_date, end_date, **kwarg):
        if type(sites) is str:
            sites = [sites]
        if 'DailyAccountQuota' not in kwarg.keys():
            kwarg['DailyAccountQuota'] = 990
        if type(kwarg['DailyAccountQuota']) is not list:
            kwarg['DailyAccountQuota'] = [kwarg['DailyAccountQuota'] for _ in range(_config['darksky_secret_key_num'])]
        if 'save_path' not in kwarg.keys():
            kwarg['save_path'] = 'data/weather/'

        self.site_list = pd.read_csv(_config['nms_data_path'] + 'site_location.csv')
        self.sites = sites
        self.start_date = start_date
        self.end_date = end_date
        self.quota = kwarg['DailyAccountQuota']
        self.save_path = kwarg['save_path']
        self.req_pool = self._requests_pool()
        self.account_num = _config['darksky_secret_key_num']
        self.secret_key_list = [_config['darksky_secret_key_' + str(i)] for i in range(self.account_num)]
        self.entry_num = len(self.req_pool)
        self.est_days = self.entry_num // sum(kwarg['DailyAccountQuota']) + 1
        ...

    def _requests_pool(self):
        req_pool = []
        for s in self.sites:
            start_time = time.mktime((time.strptime(self.start_date+' 12:00:00','%Y-%m-%d %H:%M:%S')))
            end_time = time.mktime((time.strptime(self.end_date+' 12:00:00','%Y-%m-%d %H:%M:%S')))
            lon = self.site_list.loc[self.site_list['监测点编码'] == s, '经度'].values[0]
            lat = self.site_list.loc[self.site_list['监测点编码'] == s, '纬度'].values[0]
            req_loc = '{0:f},{1:f},'.format(lat, lon)
            time_now = start_time
            while time_now <= end_time:
                req_0 = 'https://api.darksky.net/forecast/'
                req_1 = '/'+req_loc+str(round(time_now))
                file = self.save_path + s + '/' + time.strftime('%Y-%m-%d', time.gmtime(time_now)) + '.json'
                req_pool.append([req_0, req_1, file])
                time_now += 3600*24
        return req_pool

    def run(self, flag=False, shift=0, wait_hour=24):
        est_days = (self.entry_num - shift) // sum(self.quota) + 1
        print('This task include', self.entry_num - shift, ' entries and need ', est_days, ' days using ', self.account_num, ' accounts.')
        if flag:
            key = 0
            key_count = 0
            for s in self.sites:
                if not os.path.exists(self.save_path + s + '/'):
                    os.mkdir(self.save_path + s + '/')
            for i in range(shift, self.entry_num):
                if (i+1)%500 == 0:
                    print(i+1, '/', self.entry_num, '...')
                if key_count >= self.quota[key]:
                    # print('key ', key, ' run out of quota')
                    key_count = 0
                    key = key + 1
                    if key >= self.account_num:
                        key = 0
                        print('Waiting', wait_hour,' hours...')
                        print('Current shift:', i)
                        time.sleep(wait_hour*3600)
                req = self.req_pool[i]
                url = req[0] + self.secret_key_list[key] + req[1]
                file = req[2]
                try:
                    r = requests.get(url)
                    with open(file, 'w') as result_file:
                        result_file.write(r.text)
                except:
                    print('Error: ', file)
                    with open('ERROR_LOG_FILE', 'a') as ERROR_LOG_file:
                        ERROR_LOG_file.write(url+'\n')
                        ERROR_LOG_file.write(file+'\n')
                key_count = key_count + 1
        else:
            print('Use DarkSkyAPI.run(flag = True) to confirm and run task.')

    def retry_error(self, flag=False, shift=0, wait_hour=24):
        with open('ERROR_LOG_FILE', 'r') as ERROR_LOG_file:
            lines = [x.replace('\n','') for x in ERROR_LOG_file.readlines()]
            req_pool = [['https://api.darksky.net/forecast/', lines[2*i].rsplit('/',2)[-1], lines[2*i+1]] for i in range(len(lines)//2)]

        with open('ERROR_LOG_FILE_BCK', 'w') as ERROR_LOG_file:
            for line in lines:
                ERROR_LOG_file.write(line+'\n')
        with open('ERROR_LOG_FILE', 'w') as ERROR_LOG_file:
            ...

        entry_num = len(req_pool)
        est_days = (entry_num - shift) // sum(self.quota) + 1
        print('This task include', entry_num - shift, ' entries and need ', est_days, ' days using ', self.account_num, ' accounts.')

        if flag:
            key = 0
            key_count = 0
            for s in self.sites:
                if not os.path.exists(self.save_path + s + '/'):
                    os.mkdir(self.save_path + s + '/')
            for i in range(shift, entry_num):
                if (i+1)%500 == 0:
                    print(i+1, '/', entry_num, '...')
                if key_count >= self.quota[key]:
                    key_count = 0
                    key = key + 1
                    if key >= self.account_num:
                        key = 0
                        print('Waiting', wait_hour,' hours...')
                        print('Current shift:', i)
                        time.sleep(wait_hour*3600)
                req = req_pool[i]
                url = req[0] + self.secret_key_list[key] + '/' + req[1]
                file = req[2]
                try:
                    r = requests.get(url)
                    with open(file, 'w') as result_file:
                        result_file.write(r.text)
                except:
                    print('Error: ', file)
                    with open('ERROR_LOG_FILE', 'a') as ERROR_LOG_file:
                        ERROR_LOG_file.write(url+'\n')
                        ERROR_LOG_file.write(file+'\n')
                key_count = key_count + 1
        else:
            print('Use DarkSkyAPI.retry_error(flag = True) to confirm and run task.')

max_stars_repo_name: Zhiyuan-Wu/apdt

max_issues_repo_name: Zhiyuan-Wu/apdt

max_forks_repo_name: Zhiyuan-Wu/apdt


text:
#!/usr/bin/env python

"""
ADRMonitor - A program to query the temperature controller.

Purpose:
* To monitor an ADR's state via the Lakeshore 370 controller.
* To store that state in a directory of binary files ("dirfile").
* To offer the latest temperature information to local client programs.

As of Feb 6, 2014, this is still a regular script, but I will investigate
making it into a daemon (system service) that runs all the time (though
maybe not started by default at bootup??).

Author:  Joe Fowler, NIST Boulder Labs
Started: January 30, 2014
"""

import sys, time, os
import struct
import numpy as np
import pygetdata
import zmq

DEFAULT_CONFIGFILE = "/etc/adr_monitor_conf.xml"

class TemperatureDevice(object):
    """Interface to the Lakeshore 370 (or whatever your device might be
    in the future)."""

    def __init__(self, lakeshore_pad):
        import lakeshore370
        self.lakeshore = lakeshore370.Lakeshore370(pad=lakeshore_pad)

    def read_data(self):
        """Read and return (temperature,heater,time)."""
        currentTemperature = self.lakeshore.getTemperature(1)*1e3
        currentHeaterOut = self.lakeshore.getHeaterOut()
        currentTime = time.time()
        return currentTemperature, currentHeaterOut, currentTime

def current_yyyymmdd_hhmmss():
    t = time.localtime()
    return "%4.4d%2.2d%2.2d_%2.2d%2.2d%2.2d"%(
            t.tm_year, t.tm_mon, t.tm_mday, t.tm_hour, t.tm_min, t.tm_sec)

class TemperatureFaker(object):
    """Use this in place of TemperatureDevice class for testing with fake data."""
    PERIOD = 10 # seconds
    def read_data(self):
        t = time.time()
        frac_t = (t/self.PERIOD)%1
        temp = 70+2*np.sin(np.pi*2*frac_t) + np.random.standard_normal(1)[0]
        curr = 30+5*np.cos(np.pi*2*frac_t) + np.random.standard_normal(1)[0]
        return temp,curr,t

class DataWriter(object):
    """Interface to the output data writer. Currently, a dirfile writer.

    Implementation detail: for now I am using the pygetdata package to write
    the dirfile metadata (that is, the format file), but it seems easier to
    write the binary data directly through file objects with the struct.pack(...)
    function instead of using the pygetdata facilities for writing data."""

    def __init__(self, directory):
        directory = os.path.expandvars(directory)
        print "Trying directory ", directory
        basename = current_yyyymmdd_hhmmss()
        filename = os.path.join(directory, basename)
        self.dirfile = pygetdata.dirfile(filename, pygetdata.CREAT|pygetdata.RDWR)
        self.filename = filename

        # Not sure whether curfiles or soft links are now preferred in
        # kst2, so for now, set up one of each!

        # Set up a "curfile" that always points to the current dirfile.
        curfilename = os.path.join(directory, "temperatures.cur")
        prevfilename = os.path.join(directory, "temperatures.prev")
        if os.path.exists(prevfilename):
            os.unlink(prevfilename)
        if os.path.exists(curfilename):
            os.rename(curfilename, prevfilename)
        fp = open(curfilename, "w")
        fp.write("%s\n"%basename)
        fp.close()

        # Also set up a soft link to the current dirfile
        softlink = os.path.join(directory, "temperatures.lnk")
        if os.path.exists(softlink):
            os.unlink(softlink)
        os.symlink(basename, softlink)

        # Set up the dirfile fields
        fields={
            'temperature':{'type':pygetdata.FLOAT32,'spf':1},
            'current':{'type':pygetdata.FLOAT32,'spf':1},
            'ctime':{'type':pygetdata.FLOAT64,'spf':1}
                }
        for name,params in fields.iteritems():
            entry = pygetdata.entry(pygetdata.RAW_ENTRY, name, 0, params)
            self.dirfile.add(entry)
        self.dirfile.metaflush()

        # Open file handles for the raw data
        # Do this instead of self.dirfile.putdata(), because that's annoying
        self.dirfile.raw_close()
        self.raw_files = {}
        for f in fields.keys():
            fp = open(os.path.join(self.filename,f), "wb")
            self.raw_files[f] = fp

    def close(self):
        """Cleanly close the output files and flush their buffers."""
        for fp in self.raw_files.values():
            fp.close()
        self.dirfile.close()
        self.dirfile = None

    def __del__(self):
        """Make sure to call close on the output, if not already done."""
        if self.dirfile is not None:
            self.close()

    def store_data(self, temperature, current, ctime):
        results = locals()
        for f,fmt in zip(('temperature','current','ctime'),('f','f','d')):
            fp = self.raw_files[f]
            str = struct.pack(fmt, results[f])
            fp.write(str)
            fp.flush()

class ZMQServer(object):
    """A very simple data server, using ZeroMessageQueue (0MQ)
    in the interprocess communication mode, to answer requests
    received on a TCP-like socket from other local processes.

    The system avoids threading by cooperatively polling the
    socket for queries and sleeping. I hope this isn't a performance
    problem, because it sure is simple. A big drawback is that
    the self.serve_and_delay(X) blocks for X seconds. That means
    you can only use it if you know you have nothing better to
    do during that time. Anyway, it works for the ADR cryo temperature
    monitor, because the only other work is to poll the GPIB every
    X seconds (typically, X is between 0.3 and 3).

    Example usage:
    The following will bind 'ipc:///home/pcuser/wherever.../socket_name',
    which the client program or programs will also need to know.

    server = ZMQServer('/home/pcuser/wherever.../socket_name')
    while True:
        x,y,z = lakeshore_device_object.query()
        server.newdata(x,y,z)
        server.serve_and_delay(1.0)

    """

    def __init__(self, socket_name):
        '''Sets up an ipc socket at ipc://<socket_name> where
        <socket_name> must be a file-like pathname in the computer's
        filesystem, to which this process has write access.'''

        self.socket_name = "ipc://%s"%os.path.expandvars(socket_name)
        self.data = 0,0,0
        self.context = zmq.Context()
        self.socket = self.context.socket(zmq.REP)
        print "Binding to %s"%self.socket_name
        self.socket.bind(self.socket_name)

        self.poller = zmq.Poller()
        self.poller.register(self.socket, zmq.POLLIN)

    def __del__(self):
        self.socket.close()

    def newdata(self, temperature, current, ctime):
        self.data = (temperature, current, ctime)

    def serve_and_delay(self, delay_sec):
        '''For a time period of <delay_sec> seconds, poll the socket for
        queries

        '''
        t_end = time.time() + delay_sec
        while True:
            poll_interval_msec = (t_end - time.time())*1e3
            if poll_interval_msec <= 0:
                break

            poll_result = dict(self.poller.poll(poll_interval_msec))
            if self.socket in poll_result and poll_result[self.socket] == zmq.POLLIN:
                msg = self.socket.recv()
                print "\nReceived request: ", msg
                msg = "%.5f %.3f %.5f"%self.data
                self.socket.send(msg)

class ADRMonitor(object):
    """Run the ADR monitor by reading from a TemperatureDevice and
    writing to a DataWriter, plus handling a ZMQServer.

    Usage:
    m = ADRMonitor.ADRMonitor('/etc/adr_monitor_conf.xml')
    m.run()
    """

    def __init__(self, xml_config_file, simulate=False):
        from lxml import etree
        try:
            document = etree.parse(xml_config_file)
#             print etree.tostring(root, pretty_print=True)
        except:
            print "Could not parse %s"%xml_config_file

        self.doc = document
        self.root = document.getroot()
        assert self.root.get('version') == '1'
        for node in self.root:
            self.__dict__[node.tag] = node.text

        # This check helps make errors in the settings file be less cryptic to the user.
        REQUIRED_SETTINGS = ("dirfile_location","polling_period_msec",
                             "lakeshore_pad", "socket_name")
        for s in REQUIRED_SETTINGS:
            if s not in self.__dict__:
                raise IOError("xml config file %s does not contain a '%s' tag"%(
                        xml_config_file, s))

        if not simulate:
            self.device = TemperatureDevice(lakeshore_pad = int(self.lakeshore_pad))
        else:
            self.device = TemperatureFaker()
        self.writer = DataWriter(self.dirfile_location)
        self.server = ZMQServer(self.socket_name)

    def run(self):
        delay = 0.001*int(self.polling_period_msec)
        while True:
            a,b,c = self.device.read_data()
            print a,b,c
            self.writer.store_data(a,b,c)
            self.server.newdata(a,b,c)
            self.server.serve_and_delay(delay)

    def __del__(self):
        self.writer.close()

def create_settings_file(filename, **kwargs):
    """A quick way to generate a settings file. Normally not used, but
you can run this from a Python session to make a settings file.

If you are having trouble doing that, then perhaps you can just copy
the following lines into /etc/adr_monitor_conf.xml (the default place
is given by global variable DEFAULT_CONFIGFILE in this program).

<ADRMonitor version="1">
  <polling_period_msec>500</polling_period_msec>
  <dirfile_location>$HOME/data/cryo_logs</dirfile_location>
  <lakeshore_pad>13</lakeshore_pad>
  <socket_name>$HOME/data/cryo_logs/cryo_monitor</socket_name>
</ADRMonitor>
    """

    # Start with default settings, and update using any arguments here.
    settings = {'lakeshore_pad':13,
                'polling_period_msec':500,
                'dirfile_location':'$HOME/data/cryo_logs',
                'socket_name':'$HOME/data/cryo_logs/cryo_monitor'}
    settings.update(kwargs)

    # Generate an XML representation
    from lxml import etree
    root = etree.Element("ADRMonitor", version="1")
    for k,v in settings.iteritems():
        child = etree.SubElement(root, k)
        child.text = str(v)

    # Store XML to a file
    fp = open(filename, "w")
    fp.write(etree.tostring(root, pretty_print=True))
    fp.close()

def main():

    config_file = DEFAULT_CONFIGFILE
    if len(sys.argv)==2:
        config_file = sys.argv[1]
    elif len(sys.argv) > 2:
        print "usage: ADRMonitor.py [config_file]"
        print "    The default config file is %s"%DEFAULT_CONFIGFILE
        print "    If you need to construct one, see create_settings_file()"
        print "    Here is its docstring\n%s\n   "%(75*'-'),
        print create_settings_file.__doc__
        print 75*'-'
        return

    simulate = True  # Obviously, you remove this when working with a real LS370.

    if not os.path.isfile(config_file):
        if config_file == DEFAULT_CONFIGFILE:
            print "The default configuration "\
                "file %s is missing"%DEFAULT_CONFIGFILE
            print "Possibly helpful docstring from create_settings_file():\n"
            print create_settings_file.__doc__
        else:
            print "The specified configuration file %s doesn't exist."%config_file
            print "Perhaps you can try omitting the config file name and fall back"
            print "to the default?"
        return

    monitor = ADRMonitor(config_file, simulate=simulate)
    monitor.run()

if __name__=='__main__':
    main()


text:
from __future__ import division
from __future__ import print_function
from __future__ import absolute_import

import sys
if sys.version > '3':
    long = int

import numpy as np
import copy
import cPickle
from scipy.stats import chi2

#-----------------------------------------------------------------------
def airtovac(wave_air):
    """
	__author__ = 'Kyle B. Westfall'

    Wavelengths are corrected for the index of refraction of air under
    standard conditions.  Wavelength values below 2000 A will not be
    altered.  Uses formula from Ciddor 1996, Applied Optics 62, 958.

    Args:
        wave_air (int or float): Wavelength in Angstroms, scalar or
            vector. If this is the only parameter supplied, it will be
            updated on output to contain double precision vacuum
            wavelength(s).

    Returns:
        numpy.float64 : The wavelength of the line in vacuum.

    Example:
        If the air wavelength is  W = 6056.125 (a Krypton line), then
        :func:`airtovac` returns vacuum wavelength of W = 6057.8019.

    *Revision history*:
        | Written W. Landsman                November 1991
        | Use Ciddor (1996) formula for better accuracy in the infrared
        |   Added optional output vector, W Landsman Mar 2011
        | Iterate for better precision W.L./D. Schlegel  Mar 2011
        | Transcribed to python, K.B. Westfall Apr 2015

    .. note::
        Take care within 1 A of 2000 A.   Wavelengths below 2000 A *in
        air* are not altered.

    """

    # Copy the data
    wave_vac = wave_air.astype(np.float64) if hasattr(wave_air, "__len__") else float(wave_air)
    g = wave_vac > 2000.0                           # Only modify above 2000 A
    Ng = np.sum(g)

    if Ng > 0:
        # Handle both arrays and scalars
        if hasattr(wave_air, "__len__"):
            _wave_air = wave_air[g].astype(np.float64)
            _wave_vac = wave_vac[g]
        else:
            _wave_air = float(wave_air)
            _wave_vac = float(wave_vac)

        for i in range(0,2):
            sigma2 = np.square(1.0e4/_wave_vac)     #Convert to wavenumber squared
            fact = 1.0 + 5.792105e-2/(238.0185 - sigma2) + 1.67917e-3/(57.362 - sigma2)
            _wave_vac = _wave_air*fact

        if hasattr(wave_air, "__len__"):        # Save the result
            wave_vac[g] = _wave_vac
        else:
            wave_vac = _wave_vac

    return wave_vac

#-----------------------------------------------------------------------
def bisect_array(array):
	"""
	It takes an array as input and returns the bisected array :
	bisected array[i] = (array[i] + array[i+1] )/2. Its lenght is one less than the array.

	:param array: input array
	"""
	bisected_array = np.zeros(len(array) - 1)
	for ai in range(len(bisected_array)):
		bisected_array[ai] = (array[ai] + array[ai + 1])/2.0
	return bisected_array

#-----------------------------------------------------------------------
def max_pdf(probs,property,sampling):
	"""
	determines the maximum of a pdf of a property for a given sampling

	:param probs: probabilities
	:param  property: property
	:param  sampling: sampling of the property
	"""
	lower_limit 	= np.min(property)
	upper_limit 	= np.max(property)
	if upper_limit==lower_limit:
		return np.asarray(property),np.ones(len(probs))/np.size(probs)

	property_pdf_int= np.arange(lower_limit, upper_limit * 1.001, (upper_limit-lower_limit) /sampling ) + ( upper_limit - lower_limit) * 0.00000001
	prob_pdf 		= np.zeros(len(property_pdf_int))

	for p in range(len(property_pdf_int)-1):
		match_prop = np.where( (property <= property_pdf_int[p+1]) & (property > property_pdf_int[p]) )
		if np.size(match_prop) == 0:
			continue
		else:
			prob_pdf[p] = np.max( probs[match_prop] )

	property_pdf = bisect_array(property_pdf_int)
	return property_pdf,prob_pdf[:-1]/np.sum(prob_pdf)

#-----------------------------------------------------------------------
def convert_chis_to_probs(chis,dof):
	"""
	Converts chi squares to probabilities.

	:param chis: array containing the chi squares.
	:param dof: array of degrees of freedom.
	"""
	chis = chis / np.min(chis) * dof
	prob =  1.0 - chi2.cdf(chis,dof)
	prob = prob / np.sum(prob)
	return prob

#-----------------------------------------------------------------------
def light_weights_to_mass(light_weights,mass_factors):
	"""
	Uses the data/model mass-to-light ratio to convert
	SSP contribution (weights) by light into
	SSP contributions by mass.

	:param light_weights: light (luminosity) weights obtained when model fitting
	:param mass_factors: mass factors obtained when normalizing the spectrum
	"""
	mass_weights 	= np.zeros(np.shape(light_weights))
	unnorm_mass 	= np.zeros(np.shape(light_weights))
	for w in range(len(light_weights)):
		unnorm_mass[w]	= light_weights[w] * mass_factors
		mass_weights[w] = unnorm_mass[w] / np.sum(unnorm_mass[w])
	return unnorm_mass,mass_weights

#-----------------------------------------------------------------------
def find_closest(A, target):
	"""
	returns the id of the target in the array A.
	:param A: Array, must be sorted
	:param target: target value to be located in the array.
	"""
	idx = A.searchsorted(target)
	idx = np.clip(idx, 1, len(A)-1)
	left = A[idx-1]
	right = A[idx]
	idx -= target - left < right - target
	return idx

#-----------------------------------------------------------------------
def averages_and_errors(probs,prop,sampling):
	"""
	determines the average and error of a property for a given sampling

	returns : an array with the best fit value, +/- 1, 2, 3 sigma values.

	:param probs: probabilities
	:param  property: property
	:param  sampling: sampling of the property
	"""
	xdf,y = max_pdf(probs,prop,sampling)
	cdf = np.zeros(np.shape(y))
	cdf_probspace = np.zeros(np.shape(y))
	for m in range(len(y)):
		cdf[m] = np.sum(y[:m])

	cdf = cdf / np.max(cdf)
	area_probspace = y*(xdf[1]-xdf[0])
	area_probspace = area_probspace/np.sum(area_probspace)
	indx_probspace = np.argsort(area_probspace)[::-1]
	desc_probspace = np.sort(area_probspace)[::-1]

	cdf_probspace = np.zeros(np.shape(desc_probspace))
	for m in range(len(desc_probspace)):
		cdf_probspace[m] = np.sum(desc_probspace[:m])

	av_sigs = [0.6827,0.9545,0.9973] # Median, + / - 1 sig, + / - 2 sig, + / - 3 sig

	# Sorts results by likelihood and calculates confidence intervals on sorted space
	index_close = find_closest(cdf_probspace, av_sigs)

	best_fit 					= xdf[indx_probspace[0]]
	upper_onesig,lower_onesig 	= np.max(xdf[indx_probspace[:index_close[0]]]),np.min(xdf[indx_probspace[:index_close[0]]])
	upper_twosig,lower_twosig 	= np.max(xdf[indx_probspace[:index_close[1]]]),np.min(xdf[indx_probspace[:index_close[1]]])
	upper_thrsig,lower_thrsig 	= np.max(xdf[indx_probspace[:index_close[2]]]),np.min(xdf[indx_probspace[:index_close[2]]])

	if np.size(xdf) == 0:
		raise Exception('No solutions found??? FIREFLY error (see statistics.py)')

	return [best_fit,upper_onesig,lower_onesig,upper_twosig,lower_twosig,upper_thrsig,lower_thrsig]

#-----------------------------------------------------------------------
def calculate_averages_pdf(probs,light_weights,mass_weights,unnorm_mass,age,metal,sampling,dist_lum):

	"""
	Calculates light- and mass-averaged age and metallicities.
	Also outputs stellar mass and mass-to-light ratios.
	And errors on all of these properties.

	It works by taking the complete set of probs-properties and
	maximising over the parameter range (such that solutions with
	equivalent values but poorer probabilities are excluded). Then,
	we calculate the median and 1/2 sigma confidence intervals from
	the derived 'max-pdf'.

	NB: Solutions with identical SSP component contributions
	are re-scaled such that the sum of probabilities with that
	component = the maximum of the probabilities with that component.
	i.e. prob_age_ssp1 = max(all prob_age_ssp1) / sum(all prob_age_ssp1)
	This is so multiple similar solutions do not count multiple times.

	Outputs a dictionary of:
	- light_[property], light_[property]_[1/2/3]_sigerror
	- mass_[property], mass_[property]_[1/2/3]_sigerror
	- stellar_mass, stellar_mass_[1/2/3]_sigerror
	- mass_to_light, mass_to_light_[1/2/3]_sigerror
	- maxpdf_[property]
	- maxpdf_stellar_mass
	where [property] = [age] or [metal]

	:param probs: probabilities
	:param light_weights: light (luminosity) weights obtained when model fitting
	:param mass_weights: mass weights obtained when normalizing models to data
	:param unnorm_mass: mass weights obtained from the mass to light ratio
	:param age: age
	:param metal: metallicity
	:param sampling: sampling of the property
	:param dist_lum: luminosity distance in cm
	"""

	# Sampling number of max_pdf (100:recommended) from options
	log_age = np.log10(age)
	log_age[np.isnan(log_age)|np.isinf(log_age)] = 0.0

	av = {} # dictionnary where values are stored :
	av['light_age'],av['light_age_1_sig_plus'],av['light_age_1_sig_minus'], av['light_age_2_sig_plus'], av['light_age_2_sig_minus'], av['light_age_3_sig_plus'], av['light_age_3_sig_minus'] = averages_and_errors(probs,np.dot(light_weights,log_age),sampling)

	av['light_metal'], av['light_metal_1_sig_plus'], av['light_metal_1_sig_minus'], av['light_metal_2_sig_plus'], av['light_metal_2_sig_minus'], av['light_metal_3_sig_plus'], av['light_metal_3_sig_minus'] = averages_and_errors(probs, np.dot(light_weights, metal), sampling)

	av['mass_age'], av['mass_age_1_sig_plus'], av['mass_age_1_sig_minus'], av['mass_age_2_sig_plus'], av['mass_age_2_sig_minus'], av['mass_age_3_sig_plus'], av['mass_age_3_sig_minus'] = averages_and_errors(probs, np.dot(mass_weights, log_age), sampling)

	av['mass_metal'], av['mass_metal_1_sig_plus'], av['mass_metal_1_sig_minus'], av['mass_metal_2_sig_plus'], av['mass_metal_2_sig_minus'], av['mass_metal_3_sig_plus'], av['mass_metal_3_sig_minus'] = averages_and_errors(probs, np.dot(mass_weights, metal), sampling)

	conversion_factor 	= 10.0**(-17) * 4 * np.pi * dist_lum**2.0 # unit 1e-17 cm2
	tot_mass = np.log10(np.sum(unnorm_mass, 1) * conversion_factor)
	av['stellar_mass'], av['stellar_mass_1_sig_plus'], av['stellar_mass_1_sig_minus'], av['stellar_mass_2_sig_plus'], av['stellar_mass_2_sig_minus'], av['stellar_mass_3_sig_plus'], av['stellar_mass_3_sig_minus'] = averages_and_errors(probs,tot_mass,sampling)

	return av

#-----------------------------------------------------------------------
def normalise_spec(data_flux,model_flux):
	"""
	Normalises all models to the median value of the spectrum.
	Saves the factors for later use.

	Outputs : normed models and translation factors.

	:param data_flux: observed flux in the data
	:param model_flux: flux from the models
	"""
	data_norm 				= np.median(data_flux)
	num_mods 				= len(model_flux)
	model_norm,mass_factor 	= np.zeros(num_mods),np.zeros(num_mods)
	normed_model_flux 		= np.zeros((num_mods,len(model_flux[0])))

	for m in range(len(model_flux)):
		model_norm[m] 			= np.median(model_flux[m])
		mass_factor[m] 			= data_norm/model_norm[m]
		normed_model_flux[m] 	= model_flux[m] / model_norm[m] * data_norm

	return normed_model_flux,mass_factor

#-----------------------------------------------------------------------
def match_data_models( data_wave_int, data_flux_int, data_flags, error_flux_int, model_wave_int, model_flux_int, min_wave_in, max_wave_in, saveDowngradedModel = True, downgradedModelFile = "DGmodel.txt"):
	"""
	 * 0.Take data and models as inputs
	 * 1. interpolate data and model to the lowest sampled array.
		* 1.1. Defines the wavelength range on the model and on the data
		* 1.2. Downgrades the array, model or data, that has most sampling
	 	* 1.3. integrate between them to output a matched resolution array for data and model
	 * 2. Returns the matched wavelength array, the corresponding data, error and model arrays : matched_wave,matched_data,matched_error,matched_model

	:param data_wave_int: data wavelength array in the restframe
	:param data_flux_int: data flux array
	:param data_flags: data quality flag array : 1 for good data
	:param error_flux_int: data flux error array
	:param model_wave_int: model wavelength array (in the rest frame)
	:param model_flux_int: model flux array
	:param min_wave_in: minimum wavelength to be considered
	:param max_wave_in: maximum wavelength to be considered
	:param saveDowngradedModel: if True it will save the downgraded models
	:param downgradedModelFile: location where downgreaded models will be saved
	"""
	# 1. interpolate onto the bisection of lowest sampled array.
	num_models = len(model_flux_int)
	# 1.1. Defines the wavelength range on the model and on the data
	min_wave = np.max([np.min(data_wave_int[np.where(data_flags==1)]), np.min(model_wave_int),min_wave_in])
	max_wave = np.min([np.max(data_wave_int[np.where(data_flags==1)]), np.max(model_wave_int),max_wave_in])
	#print np.min(data_wave_int[np.where(data_flags==1)]), np.min(model_wave_int), min_wave_in
	#print np.max(data_wave_int[np.where(data_flags==1)]), np.max(model_wave_int), max_wave_in
	loc_model 	= np.array(( model_wave_int <= max_wave) & (model_wave_int >= min_wave))
	if np.sum(loc_model)==0:
		raise ValueError("The wavelength range input is below or above model wavelength coverage!")
	model_wave 	= model_wave_int[loc_model]
	num_mod  	= np.sum(loc_model)
	model_flux 	= np.zeros((num_models,num_mod))

	for m in range(num_models):
		model_flux[m] = model_flux_int[m][loc_model]

	loc_data 	= np.array(( data_wave_int <= max_wave) & (data_wave_int >= min_wave))
	if np.sum(loc_data)==0:
		raise ValueError("The wavelength range input is below or above data wavelength coverage!")
	num_dat  	= np.sum(loc_data)
	data_wave 	= data_wave_int[loc_data]
	data_flux 	= data_flux_int[loc_data]
	error_flux 	= error_flux_int[loc_data]
	# 1.2. Downgrades the array, model or data, that has most sampling
	if num_mod >= num_dat:
		#print "More model points than data points! Downgrading models to data sampling ..."
		bisect_data = bisect_array(data_wave) + np.min(data_wave)*0.0000000001
		matched_model = np.zeros((num_models,len(bisect_data) - 1))
		for m in range(num_models):
			model_flux_bounds 	= np.interp(bisect_data, model_wave, model_flux[m])
			combined_wave_int 	= np.concatenate((model_wave,bisect_data))
			combined_flux_int 	= np.concatenate((model_flux[m],model_flux_bounds))
			sort_indices 		= np.argsort(combined_wave_int)

			combined_wave 		= np.sort(combined_wave_int)
			boundary_indices 	= np.searchsorted(combined_wave,bisect_data)
			combined_flux 		= combined_flux_int[sort_indices]

			len_combo = len(combined_flux)
			# 1.3. produces a matched resolution array
			for l in range(len(boundary_indices) - 1):
				if boundary_indices[l + 1] >= len_combo:
					matched_model[m][l] = matched_model[m][l - 1]
				else:
					matched_model[m][l] = np.trapz(combined_flux[boundary_indices[l] : boundary_indices[l + 1] +  1], x=combined_wave[boundary_indices[l] :boundary_indices[l + 1] + 1]) / (combined_wave[boundary_indices[l + 1]] -  combined_wave[boundary_indices[l] ])

		matched_wave = data_wave[1:-1]
		matched_data = data_flux[1:-1]
		matched_error = error_flux[1:-1]
		# OPTION : saves the downgraded models.
		if saveDowngradedModel:
			#print "saving downgraded models to ",downgradedModelFile
			f.open(downgradedModelFile,'w')
			cPickle.dump([matched_wave, matched_data, matched_error],f)
			f.close()

	else:
		#print "More data points than model points! Downgrading data to model sampling ..."
		bisect_model = bisect_array(model_wave) + np.min(model_wave)*0.0000000001
		boundaries 	= np.searchsorted(data_wave,bisect_model)

		data_flux_bounds 	= np.interp(bisect_model, data_wave, data_flux)
		error_flux_bounds 	= np.interp(bisect_model, data_wave, error_flux)
		combined_wave_int 	= np.concatenate((data_wave,bisect_model))
		combined_flux_int 	= np.concatenate((data_flux,data_flux_bounds))
		combined_error_int 	= np.concatenate((data_flux,error_flux_bounds))
		sort_indices 		= np.argsort(combined_wave_int)

		combined_wave 		= np.sort(combined_wave_int)
		boundary_indices 	= np.searchsorted(combined_wave,bisect_model)
		combined_flux 		= combined_flux_int[sort_indices]
		combined_error 		= combined_error_int[sort_indices]

		# 1.3. produces a matched resolution array
		matched_data,matched_error= np.zeros(len(boundary_indices) - 1),np.zeros(len(boundary_indices) - 1)

		len_combo = len(combined_flux)
		for l in range(len(boundary_indices) - 1):
			if boundary_indices[l + 1] >= len_combo:
				matched_data[l] 	= matched_data[l - 1]
				matched_error[l] 	= matched_error[l - 1]
			else:
				matched_data[l] 	= np.trapz(combined_flux[boundary_indices[l]:boundary_indices[l + 1] + 1], x=combined_wave[boundary_indices[l]: boundary_indices[l + 1] + 1])/  (combined_wave[boundary_indices[l + 1]] - combined_wave[boundary_indices[l]])
				matched_error[l] 	= np.trapz(combined_error[boundary_indices[l]:boundary_indices[l + 1] + 1], x=combined_wave[boundary_indices[l]:boundary_indices[l + 1] + 1])/ (combined_wave[boundary_indices[l + 1]] - combined_wave[boundary_indices[l]])

		matched_wave 		= model_wave[1:-1]
		matched_model 		= np.zeros((num_models,len(matched_wave)))
		for m in range(num_models):
			matched_model[m][:] = model_flux[m][1:-1]

	return matched_wave,matched_data,matched_error,matched_model

max_stars_repo_name: AndresSixtos/pyeBOSS

max_issues_repo_name: AndresSixtos/pyeBOSS

max_forks_repo_name: AndresSixtos/pyeBOSS

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
from __future__ import division
from __future__ import print_function

import os
import sys
import unittest

import numpy as np

from pyhealth.evaluation.binaryclass import evaluator as binary_evalua
from pyhealth.evaluation.multiclass import evaluator as multiclass_evalua
from pyhealth.evaluation.multilabel import evaluator as multilabel_evalua
from pyhealth.evaluation.evaluator import func as evaluator
from pyhealth.data.simulation_data import generate_label

class TestEvaluation(unittest.TestCase):

    def _simulate_binary_data(self, N_case):
        y = generate_label(N_case, task = 'binaryclass').reshape(-1, 1)
        hat_y = np.random.rand(N_case).reshape(-1, 1)
        return hat_y, y

    def _simulate_multiclass_data(self, N_case, N_class):
        y = generate_label(N_case, task = 'multiclass', n_class = N_class)
        hat_y = np.random.rand(N_case, N_class)
        return hat_y, y

    def _simulate_multilabel_data(self, N_case, N_class):
        y = generate_label(N_case, task = 'multilabel', n_class = N_class)
        hat_y = np.random.rand(N_case, N_class)
        return hat_y, y

    def test_binary(self):
        hat_y, y = self._simulate_binary_data(100)
        assert np.shape(hat_y) == np.shape(y)
        assert int(np.sum(np.sum(y, 1))) > 0 and int(np.sum(np.sum(y, 1))) < 100
        results = binary_evalua(hat_y, y)

    def test_multiclass(self):
        hat_y, y = self._simulate_multiclass_data(100, 3)
        assert np.shape(hat_y) == np.shape(y)
        assert int(np.sum(np.sum(y, 1))/100) == 1
        results = multiclass_evalua(hat_y, y)

    def test_multilabel(self):
        hat_y, y = self._simulate_multilabel_data(100, 3)
        assert np.shape(hat_y) == np.shape(y)
        assert np.sum(np.sum(y, 1))/100. > 1.
        results = multilabel_evalua(hat_y, y)

    def test_evaluator(self):
        hat_y, y = self._simulate_binary_data(100)
        unified_r = evaluator(hat_y, y)
        be_r = binary_evalua(hat_y, y)
        for key in unified_r.keys():
            assert unified_r[key] == be_r[key]

        hat_y, y = self._simulate_multiclass_data(100, 3)
        unified_r = evaluator(hat_y, y)
        me_r = multiclass_evalua(hat_y, y)
        for key in unified_r.keys():
            assert unified_r[key] == me_r[key]

        hat_y, y = self._simulate_multilabel_data(100, 3)
        unified_r = evaluator(hat_y, y)
        me_r = multilabel_evalua(hat_y, y)
        for key in unified_r.keys():
            assert unified_r[key] == me_r[key]

max_stars_repo_name: rkalahasty/PyHealth

max_issues_repo_name: rkalahasty/PyHealth

max_forks_repo_name: rkalahasty/PyHealth


text:
import sys
import nltk
import string, re
import csv, pandas
import numpy

"""feature-extract.py: extract data with a features list"""
__author__ = "YuanSun"

# main
def main(feature_list_file, input_file, output_file):
    df_feature = pandas.read_csv(feature_list_file)
    feature_list = df_feature.values[:,0]
    feature_list = numpy.unique(feature_list).astype(numpy.str)
    print len(feature_list)

    colnames = ['label', 'id','message','source']
    data = pandas.read_csv(input_file, names=colnames)
    message = data.message.tolist()
    label = data.label.tolist()
    ids = data.id.tolist()

    # clean data, change urls to the string "URLS" and numbers to the string "NUMBERS"
    text = [re.sub("http\S+", 'URLS', t) for t in message]
    me_list = [re.sub('[%s]' % string.digits, ' NUMBERS', s) for s in text ]

    fex_ma = []
    for me in me_list:
        fextract_list = []
        for fe in feature_list:
            if fe in me:
                fextract_list.extend("1")
            else:
                fextract_list.extend("0")
        fex_ma.append(fextract_list)

    print len(label)
    print len(ids)
    print len(fex_ma)

    all_list = numpy.column_stack((label, ids, fex_ma))
    print all_list
    with open(output_file, 'w') as fp:
        a = csv.writer(fp, delimiter=',')
        a.writerows(all_list)

if __name__ == '__main__':
    main(sys.argv[1], sys.argv[2], sys.argv[3])


text:
from cjax.utils.math_trees import *
from jax.experimental.optimizers import l2_norm
import pickle
from jax.flatten_util import ravel_pytree
from jax.tree_util import tree_flatten
import matplotlib.pyplot as plt
import numpy as onp
from jax.experimental import stax
from jax.nn.initializers import zeros, orthogonal, delta_orthogonal, uniform
from jax.nn import sigmoid, hard_tanh, soft_sign, relu
from jax.experimental.stax import Dense, Sigmoid, Softplus, Tanh, LogSoftmax,Relu
from cjax.optimizer.optimizer import GDOptimizer, AdamOptimizer, OptimizerCreator
import numpy.random as npr
from jax import random
from jax import jit, vmap, grad
from cjax.utils.abstract_problem import AbstractProblem
from jax.tree_util import tree_map
from cjax.utils.math_trees import *
from cjax.utils.evolve_utils import *
from cjax.continuation.states.state_variables import StateWriter
import math
import mlflow
import json
from cjax.utils.custom_nn import constant_2d, HomotopyDense, v_2d
from cjax.utils.datasets import mnist, get_mnist_data, meta_mnist, mnist_preprocess_cont
import pickle
import jax.numpy as np
from examples.torch_data import get_data

def accuracy(params, bparam, batch):
    x, targets = batch
    x = np.reshape(x, (x.shape[0], -1))
    target_class = np.argmax(targets, axis=-1)
    predicted_class = np.argmax(predict_fun(params, x, bparam=bparam[0], activation_func=relu), axis=-1)
    return np.mean(predicted_class == target_class)

def objective(params, bparam, batch) -> float:
    x, targets = batch
    x = np.reshape(x, (x.shape[0], -1))
    logits = predict_fun(params, x, bparam=bparam[0], activation_func=relu)
    loss = -np.mean(np.sum(logits * targets, axis=1))
    loss += 5e-7 * (l2_norm(params)) #+ l2_norm(bparam))
    return loss

path = '/opt/ml/mlruns/5/5a518d1762cb4b1eafc73e5cfdb11760/artifacts/output/parc_6_i3.pkl'
path_std = "/opt/ml/mlruns/5/ddb66f98e43647c2a0c65489569ba9f7/artifacts/output/params.pkl"
with open(path_std, 'rb') as fp:
    my_params = pickle.load(fp)

# init_fun, predict_fun = stax.serial(
#         HomotopyDense(out_dim=18, W_init=orthogonal(), b_init=zeros),
#         Dense(out_dim=10, W_init=orthogonal(), b_init=zeros), LogSoftmax
#     )

init_fun, predict_fun = stax.serial(
        Dense(out_dim=18), Relu,
        Dense(out_dim=10), LogSoftmax
    )

bparam = [0.019664965569972992]
bparam = [1.1531]
bparam = [1.0]
train_images, train_labels, test_images, test_labels = mnist(
            permute_train=False, resize=True, filter=False
        )

print(len(test_labels))
val_loss = objective(my_params, bparam, (test_images, test_labels))
val_acc = accuracy(my_params, bparam, (test_images, test_labels))
print(val_loss, val_acc)

# _, _, test_images, test_labels = mnist(
#             permute_train=False, resize=True, filter=True
#         )
ae_shape, ae_params = init_fun(random.PRNGKey(0), (10000, 36))
print(tree_util.tree_structure(ae_params))

print(tree_util.tree_structure(my_params[0]))
new_params = []
print(len(my_params[0][0]))
print(len(my_params[1]))
for i in my_params[0][0]:
    new_params.append(np.asarray(i))
W = np.asarray(new_params).reshape(36,18)
b = np.asarray(my_params[0][1]).reshape(18)
one  = (W,b)
print(tree_util.tree_structure(one))
new_params = []
p1 = my_params[1]
for i in my_params[2][0]:
    new_params.append(np.asarray(i))
W = np.asarray(new_params).reshape(18,10)
b = np.asarray(my_params[2][1]).reshape(10)
two  = (W,b)
print(tree_util.tree_structure(two))
three = my_params[3]

final_params = [one, p1,two, three]
print(tree_util.tree_structure(final_params))

# _, def_tree = ravel_pytree(ae_params)
# leaves, _ = ravel_pytree(params)
#
# params = def_tree(leaves)
# print(tree_util.tree_structure(params))
#
val_loss = objective(final_params, bparam, (test_images, test_labels))
val_acc = accuracy(final_params, bparam, (test_images, test_labels))
print(val_loss, val_acc)

reason: import numpy,import jax,from jax