text:
"""Pure python example on quartic potential, for illustration and to identify
implementation errors / pitfalls that occur independent of the complexity of
OpenMM (numerical precision, constraints, CustomIntegrator compilation, etc.).

Current problems:
* Shadow work in first step appears to be negative for some integrator schemes
  (OVRVO and RVOVR)
* Weird jump in accumulated shadow work after first step, even when using
  null midpoint operator (for all schemes). For the null midpoint operator,
  it should be in steady state / the first step shouldn't be statistically
  any different from any of the following steps...
"""

import numpy as np
import matplotlib
from numba import jit

matplotlib.use('agg')
import matplotlib.pyplot as plt
from functools import partial
from scipy.stats import entropy
import statsmodels.api as sm

# import pymbar
np.random.seed(0)

figure_directory = "figures/"  # relative to script
figure_format = ".jpg"

# Define system
beta = 1.0  # inverse temperature
dim = 1  # system dimension

def potential(x):
    return x**4

def reduced_potential(x): return potential(x) * beta

@jit
def log_q(x): return - reduced_potential(x)

@jit
def q(x): return np.exp(log_q(x))

def force(x):
    return - 4.0 * x ** 3

# normalized density
x = np.linspace(-3, 3, 1000)
Z = np.trapz(map(q, x), x)

def p(x): return q(x) / Z

def plot(x, y, color='blue', padding=0):
    plt.plot(x, y, alpha=0.5, linewidth=2, color=color)
    plt.fill_between(x, y, alpha=0.25, color=color)
    span = max(y) - min(y)
    plt.ylim(min(y) - (padding * span), max(y) + (padding * span))

def savefig(name):
    plt.savefig("{}{}{}".format(figure_directory, name, figure_format), dpi=300)

def plot_system():
    print("Plotting system...")

    plt.figure()
    plt.subplot(1, 3, 1)
    plot(x, map(reduced_potential, x))
    plt.title('Reduced potential\n$u(x)=\\beta U(x)$')

    plt.subplot(1, 3, 2)
    plot(x, map(p, x))
    plt.title('Equilibrium\ndistribution\n$p(x)\propto e^{-u(x)}$')

    plt.subplot(1, 3, 3)
    plot(x, map(force, x))
    plt.title('Force\n$F(x)=-\\nabla U(x)$')

    plt.tight_layout()
    savefig("quartic_system")
    plt.close()

# example initial conditions
x_0, v_0 = np.zeros(dim), np.ones(dim)

m = 10.0  # mass
velocity_scale = np.sqrt(1.0 / (beta * m))
timestep = 1.0
gamma = 100.0

def draw_velocities():
    return velocity_scale * np.random.randn(dim)

@jit
def kinetic_energy(v):
    return np.sum(0.5 * m * v ** 2)

@jit
def total_energy(x, v):
    return kinetic_energy(v) + potential(x)

@jit
def R_map(x, v, h):  # linear "drift"
    x_ = x + (h * v)
    v_ = v
    return x_, v_

@jit
def V_map(x, v, h):  # linear "kick"
    x_ = x
    v_ = v + (h * force(x) / m)
    return x_, v_

@jit
def O_map(x, v, h, gamma):  # Ornstein-Uhlenbeck
    x_ = x
    a = np.exp(-gamma * h)
    b = np.sqrt(1 - np.exp(-2 * gamma * h))
    v_ = (a * v) + b * draw_velocities()
    return x_, v_

splitting_map = {'R': R_map, 'V': V_map, 'O': O_map}

@jit
def langevin_map(x, v, h, splitting, gamma):
    n_R = sum([l == "R" for l in splitting])
    n_V = sum([l == "V" for l in splitting])
    n_O = sum([l == "O" for l in splitting])

    for step in splitting:
        if step == "R":
            x, v = R_map(x, v, h / n_R)
        elif step == "V":
            x, v = V_map(x, v, h / n_V)
        elif step == "O":
            x, v = O_map(x, v, h / n_O, gamma)
    return x, v

# now, we just wrap each of the R, V, O updates with a couple lines to measure
# the energy difference during the substep
@jit
def bookkeeping_R_map(x, v, h):
    ''' deterministic position update '''
    x_, v_ = R_map(x, v, h)
    e_diff = potential(x_) - potential(x)
    return x_, v_, e_diff

@jit
def bookkeeping_V_map(x, v, h):
    ''' deterministic velocity update '''
    x_, v_ = V_map(x, v, h)
    e_diff = kinetic_energy(v_) - kinetic_energy(v)
    return x_, v_, e_diff

@jit
def bookkeeping_O_map(x, v, h, gamma):
    ''' stochastic velocity update'''
    x_, v_ = O_map(x, v, h, gamma)
    e_diff = kinetic_energy(v_) - kinetic_energy(v)
    return x_, v_, e_diff

bookkeeping_map = {'R': bookkeeping_R_map,
                   'V': bookkeeping_V_map,
                   'O': bookkeeping_O_map}

def bookkeeping_langevin_map(x, v, h, splitting, gamma):
    W_shad, Q = 0, 0
    n_R = sum([l == "R" for l in splitting])
    n_V = sum([l == "V" for l in splitting])
    n_O = sum([l == "O" for l in splitting])

    for i in range(len(splitting)):
        substep = splitting[i]

        if substep == "R":
            h_substep = h / n_R
            x, v, e_diff = bookkeeping_R_map(x, v, h_substep)

        elif substep == "V":
            h_substep = h / n_V
            x, v, e_diff = bookkeeping_V_map(x, v, h_substep)

        elif substep == "O":
            h_substep = h / n_O
            x, v, e_diff = bookkeeping_O_map(x, v, h_substep, gamma)

        if substep in "RV":
            W_shad += e_diff
        else:
            Q += e_diff

    return x, v, W_shad, Q

def bookkeeping_langevin_factory(n_steps, splitting='VRORV'):
    def multistep_langevin(x, v, h, gamma):
        xs, vs = np.zeros((n_steps, len(x))), np.zeros((n_steps, len(x)))
        xs[0], xs[0] = x, v
        W_shad, Q = np.zeros(n_steps), np.zeros(n_steps)
        for i in range(1, n_steps):
            xs[i], vs[i], W_shad[i], Q[i] = \
                bookkeeping_langevin_map(xs[i - 1], vs[i - 1], h, splitting, gamma)
        return xs, vs, W_shad, Q

    return multistep_langevin

def test_scheme(x_0, v_0, n_steps, scheme, h, gamma):
    langevin = bookkeeping_langevin_factory(n_steps, scheme)
    xs, vs, W_shad, Q = langevin(x_0, v_0, h, gamma)

    # plot the x marginal
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.title('Sampled x ({})'.format(scheme))
    plt.hist(xs.flatten(), bins=50, normed=True, histtype='stepfilled', alpha=0.5);
    x = np.linspace(-3, 3, 1000)
    plt.plot(x, map(p, x))
    plt.xlabel(r'$x$')
    plt.ylabel('Sampled P(x)')

    # plot the trajectory
    plt.subplot(1, 2, 2)
    plt.title('Trajectory ({})'.format(scheme))
    plt.plot(xs[:, 0], vs[:, 0] * m, linewidth=0.5)
    plt.xlabel(r'$x$')
    plt.ylabel(r'$mv$')

    plt.tight_layout()

    savefig("quartic_samples_{}".format(scheme))
    plt.close()

    # plot the accumulated shadow work
    plt.figure()
    plt.plot(np.cumsum(W_shad))
    plt.xlabel('Step')
    plt.ylabel('Accumulated W_shad')
    plt.title('Accumulated shadow work ({})'.format(scheme))
    savefig("quartic_shadow_work_{}".format(scheme))
    plt.close()

    return xs, vs, W_shad, Q

#schemes = "VRORV RVOVR OVRVO".split()
#results = dict()
#print("Simulating long trajectories using {} integrators".format(schemes))
#for scheme in schemes:
#    results[scheme] = test_scheme(np.random.randn(dim), draw_velocities(),
#                                  n_steps=int(1e4), scheme=scheme, h=timestep,
#                                  gamma=gamma
#                                  )

# let's confirm that energy change accumulated by the bookkeeper is the same as the
# actual energy change
#print("Confirming that bookkeeping is self-consistent...")

def check_bookkeeping(results):
    for scheme in results:
        xs, vs, W_shad, Q = results[scheme]
        es = np.array([total_energy(xs[i], vs[i]) for i in range(len(xs))])

        # accumulated
        delta_e_bk = np.cumsum(W_shad + Q)
        accumulated = delta_e_bk - delta_e_bk[0]

        actual = es - es[0]
        discrepancy = accumulated - actual

        print('\t{}: Mean squared discrepancy between '
              'accumulated and actual energy change: {:.5f}'.format(scheme, np.mean(np.linalg.norm(discrepancy))))

#check_bookkeeping(results)

### Estimate nonequilibrium free energy difference
def collect_equilibrium_samples():
    print("Collecting equilibrium samples...")
    # collect equilibrium samples
    import emcee

    n_walkers = 2 * dim
    sampler = emcee.EnsembleSampler(nwalkers=n_walkers, dim=dim, lnpostfn=log_q)
    sampler.run_mcmc(pos0=np.random.randn(n_walkers, dim), N=50000)
    equilibrium_samples = sampler.flatchain[1000:]
    draw_configuration = lambda: equilibrium_samples[np.random.randint(len(equilibrium_samples))]
    plt.figure()
    plt.title('Equilibrium samples')
    plt.hist(equilibrium_samples.flatten(), bins=50, normed=True, histtype='stepfilled', alpha=0.5);
    x = np.linspace(-3, 3, 1000)
    plt.plot(x, map(p, x))
    plt.xlabel(r'$x$')
    plt.ylabel('Sampled P(x)')
    savefig("quartic_equilibrium_samples")
    plt.close()

def run_protocol(integrator, midpoint_operator):
    # simulate starting from equilibrium sample
    x_traj, v_traj, W_shad, Q = integrator(draw_configuration(), draw_velocities())
    W_shad_F = np.cumsum(W_shad)

    # apply midpoint operator
    if midpoint_operator == "null":
        x_1, v_1 = x_traj[-1], v_traj[-1]
    elif midpoint_operator == "randomize-velocity":
        x_1, v_1 = x_traj[-1], draw_velocities()
    elif midpoint_operator == "randomize-positions":
        x_1, v_1 = draw_configuration(), v_traj[-1]
    else:
        raise (ValueError("Invalid midpoint operator"))

    # simulate starting from out-of-equilibrium ensemble
    _, _, W_shad, _ = integrator(x_1, v_1)
    W_shad_R = np.cumsum(W_shad)

    return W_shad_F, W_shad_R, x_1

def perform_benchmark(integrator, midpoint_operator, n_samples):
    samples = [run_protocol(integrator, midpoint_operator) for _ in range(n_samples)]
    W_shads_F, W_shads_R, x = [s[0] for s in samples], [s[1] for s in samples], [s[2] for s in samples]
    return W_shads_F, W_shads_R, x

def estimate_nonequilibrium_free_energy(W_shads_F, W_shads_R, N_eff=None, verbose=True):
    if N_eff == None: N_eff = len(W_shads_F)
    DeltaF_neq = 0.5 * (np.mean(W_shads_F) - np.mean(W_shads_R))
    # squared uncertainty = [var(W_0->M) + var(W_M->2M) - 2 cov(W_0->M, W_M->2M)] / (4 N_eff)
    sq_uncertainty = (np.var(W_shads_F) + np.var(W_shads_R) - 2 * np.cov(W_shads_F, W_shads_R)[0, 1]) / (4 * N_eff)
    if verbose:
        print(
        "\t<W_F>: {:.5f} +/- {:.5f}".format(np.mean(W_shads_F), 1.96 * np.std(W_shads_F) / np.sqrt(len(W_shads_F))))
        print(
        "\t<W_R>: {:.5f} +/- {:.5f}".format(np.mean(W_shads_R), 1.96 * np.std(W_shads_R) / np.sqrt(len(W_shads_R))))
        print("\tcov(W_F, W_R): {:.5f}".format(np.cov(W_shads_F, W_shads_R)[0, 1]))
    return DeltaF_neq, sq_uncertainty

## KDE-based, via <U> - S
def estimate_entropy_kde_1d(samples):
    dens = sm.nonparametric.KDEUnivariate(samples)
    dens.fit()
    return dens.entropy

def estimate_F(X, U):
    '''

    Define free energy, given an energy function U and samples
     from P(x) = 1/Z e^{-U(x)} as the
    average energy - the entropy:

    F_P = \langle U \rangle_P - S_P

    Parameters
    ----------
    X : samples
    U : reduced potential function

    Returns
    -------
    F_hat : <U> - S
    '''

    # compute average energy
    U_X = np.array([U(x) for x in X])
    mean_U = np.mean(U_X)

    # compute entropy
    S = estimate_entropy_kde_1d(X)

    return mean_U - S

def estimate_dFNeq_kde(X_neq, X_eq, U):
    return estimate_F(X_neq, U) - estimate_F(X_eq, U)

## work-based

# def BAR_estimator(W_shad, W_ss):
#    return pymbar.BAR(W_shad,W_ss)[0]

def plot_mean_and_95_confidence_bands(array, x=None, color="blue", label=None):
    """Given a list of trajectories, plot the mean trajectory +/- 1.96 * standard error"""
    mean = np.mean(array, 0)
    band = 1.96 * np.std(array, 0) / np.sqrt(len(array))

    plt.plot(mean, color=color, label=label)
    if x == None: x = range(len(array[0]))
    plt.fill_between(x, mean - band, mean + band, color=color, alpha=0.3)

def unpack_W_shads(W_shads_F, W_shads_R):
    return [F[-1] for F in W_shads_F], [R[-1] for R in W_shads_R]

# okay, now we should also do the comparison -- apply each estimator

if __name__ == "__main__":
    n_steps = 20
    results = dict()
    midpoint_operators = ["null", "randomize-velocity", "randomize-positions"]
    colors = {"null": "blue",
              "randomize-velocity": "green",
              "randomize-positions": "orange"
              }

    for scheme in schemes:
        print("\n\nScheme: {}".format(scheme))

        integrator = partial(bookkeeping_langevin_factory(n_steps, scheme), h=timestep, gamma=gamma)

        results[scheme] = dict()
        plt.figure()
        for midpoint_operator in midpoint_operators:
            print("\nMidpoint operator: {}".format(midpoint_operator))
            W_shads_F, W_shads_R, x = perform_benchmark(integrator, midpoint_operator, n_samples=50000)
            results[scheme][midpoint_operator] = W_shads_F, W_shads_R, x

            ax1 = plt.subplot(1, 2, 1)
            plot_mean_and_95_confidence_bands(W_shads_F, color=colors[midpoint_operator], label=midpoint_operator)
            plt.title("{} for {} steps,\nstarting in equilibrium".format(scheme, len(W_shads_F[0])))
            plt.xlabel('steps')
            plt.ylabel('W_shad')
            ax2 = plt.subplot(1, 2, 2, sharey=ax1)
            plot_mean_and_95_confidence_bands(W_shads_R[1:], color=colors[midpoint_operator], label=midpoint_operator)

            title = "{} for another {} steps".format(scheme, len(W_shads_R[0]))
            if midpoint_operator != "null":
                title += ",\n after applying midpoint operator"
            plt.title(title)
            plt.xlabel('steps')

            DeltaF_neq, sq_uncertainty = estimate_nonequilibrium_free_energy(*unpack_W_shads(W_shads_F, W_shads_R))

            plt.legend(fancybox=True, loc='best')

            print("\tDeltaF_neq: {:.5f} +/- {:.5f}".format(DeltaF_neq, np.sqrt(sq_uncertainty)))

        savefig("benchmark_{}".format(scheme))
        plt.close()

        # plot midprotocol x samples
        x_samples = np.vstack([results[scheme][midpoint_operator][-1] for midpoint_operator in midpoint_operators])
        plt.figure()
        plt.title('Sampled x ({})'.format(scheme))
        plt.hist(x_samples.flatten(), bins=50, normed=True, histtype='stepfilled', alpha=0.5);
        x = np.linspace(-3, 3, 1000)
        plt.plot(x, map(p, x))
        plt.xlabel(r'$x$')
        plt.ylabel('Sampled P(x)')
        savefig("quartic_midprotocol_samples_{}".format(scheme))
        plt.close()

        # let's also do the comparison in section 5 of https://arxiv.org/pdf/1203.5428.pdf

max_stars_repo_name: choderalab/integrator-benchmark

max_issues_repo_name: choderalab/integrator-benchmark

max_forks_repo_name: choderalab/integrator-benchmark

reason: import numpy,from scipy,from numba,import statsmodels


text:
import argparse
import pickle
import os
import pandas as pd

from DEwNF.utils import get_split_idx_on_day
from time import time

import numpy as np
from statsmodels.nonparametric.kernel_density import KDEMultivariateConditional
from sklearn.preprocessing import StandardScaler

def main(args):
    # Notebook experiment settings
    experiment_name = args.experiment_name
    experiment_results_folder = args.results_folder
    results_path = os.path.join("../", experiment_results_folder)
    data_folder = args.data_folder
    data_file = args.data_file
    extra_data_file = args.extra_data_file

    file_name_test = f"{experiment_name}_test.pickle"
    file_path_test = os.path.join(results_path, file_name_test)

    print(f"Saving: {file_name_test}")
    with open(file_path_test, 'wb') as f:
        pickle.dump("test", f)

    # Data settings
    obs_cols = args.obs_cols

    # Load data
    csv_path = os.path.join(data_folder, data_file)
    donkey_df = pd.read_csv(csv_path, parse_dates=[4, 11])

    csv_path = os.path.join(data_folder, extra_data_file)
    extra_df = pd.read_csv(csv_path, parse_dates=[4, 12])

    # Data prep
    train_idx, test_idx = get_split_idx_on_day(donkey_df)

    # Create full data
    test_data = donkey_df.loc[test_idx, obs_cols]
    train_data = pd.concat((donkey_df.loc[train_idx, obs_cols], extra_df.loc[:, obs_cols]))

    # Normalize data
    obs_scaler = StandardScaler().fit(train_data)
    scaled_train_data = obs_scaler.transform(train_data)

    # Create conditional variable
    hours = pd.concat((donkey_df.loc[train_idx, :], extra_df.loc[:, :])).merge_date.dt.hour.values
    hours = np.expand_dims(hours, 1)

    scaled_test_data = obs_scaler.transform(test_data)

    statsmods = KDEMultivariateConditional(endog=scaled_train_data, exog=hours, indep_type='o', dep_type='cc',
                                           bw='cv_ml')

    results_dict = {'model': statsmods, }

    file_name = f"{experiment_name}.pickle"
    file_path = os.path.join(results_path, file_name)

    print(f"Saving: {file_name}")
    with open(file_path, 'wb') as f:
        pickle.dump(results_dict, f)

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    # Folder args
    parser.add_argument("--experiment_name", help="Name of the experiment, used in naming outputs")
    parser.add_argument("--results_folder", help="folder for results")
    parser.add_argument("--data_folder", help="Folder with the data")
    parser.add_argument("--data_file", help="filename for the data")
    parser.add_argument("--extra_data_file", help="file name for the extra data for semisupervised learning")

    # Data args
    parser.add_argument("--obs_cols", nargs="+", help="The column names for the observation data")

    args = parser.parse_args()

    start = time()
    main(args)
    end = time()
    print(f"finished elapsed time: {end-start}")

ext: py

max_stars_repo_name: MathiasNT/DEwNF

max_issues_repo_name: MathiasNT/DEwNF

max_forks_repo_name: MathiasNT/DEwNF

reason: import numpy,from statsmodels


text:
def avi_to_fits(data_list_entry):

	import os
	from astropy.io import fits
	import numpy as np
	import cv2
	from PIL import Image

	l = 1

	if data_list_entry.data_mode == "single":
		n = 1
		if data_list_entry.data_filedata.endswith(".avi"):
			filepath = os.path.dirname(os.path.realpath(data_list_entry.data_filedata))
			actual_file = os.path.basename(data_list_entry.data_filedata)
			video_capture = cv2.VideoCapture(data_list_entry.data_filedata)

			filename = actual_file.replace(".avi", "")
			folder = filepath+"/"+filename

			if not os.path.exists(folder):
				os.mkdir(folder)

			if data_list_entry.data_type == "raw":
				imtype = "frame"
			else:
				imtype = data_list_entry.data_type

			while True:
				filename_initial = folder+"/"+str(imtype)+"_"+str(l)+"_"+str(n)+".fits"
				if os.path.isfile(filename_initial):
					l = l + 1
				else:
					break

			while True:
				ret, frame = video_capture.read()
				if (ret != True):
					break

				filename_tif = folder+"/"+str(imtype)+"_"+str(l)+"_"+str(n)+".tif"
				filename_fits = filename_tif.replace(".tif", ".fits")

				cv2.imwrite(filename_tif, frame)

				hdu = fits.PrimaryHDU()
				im = Image.open(filename_tif)
				hdu.data = np.array(im)
				hdu.writeto(filename_fits, overwrite=True)

				if data_list_entry.state == False:
					os.remove(filename_tif)

				n += 1
		else:
			wrn_dialog = Gtk.Dialog(self, 0, Gtk.MessageType.WARNING, Gtk.ButtonsType.OK, "File not AVI")
			wrn_dialog.format_secondary_text("The selected file is not an avi and cannot be split into FITS")
			wrn_dialog.run()
			wrn_dialog.destroy()
			return(1)

	elif data_list_entry.data_mode == "group":
		for file in os.listdir(data_list_entry.data_filedata):
			n = 1

			if file.endswith(".avi"):
				filepath = data_list_entry.data_filedata
				actual_file = data_list_entry.data_filedata+"/"+file
				video_capture = cv2.VideoCapture(actual_file)

				filename = file.replace(".avi", "")
				folder = filepath+"/"+filename

				if not os.path.exists(folder):
					os.mkdir(folder)

				if data_list_entry.data_type == "raw":
					imtype = "frame"
				else:
					imtype = data_list_entry.data_type

				while True:
					filename_initial = folder+"/"+str(imtype)+"_"+str(l)+"_"+str(n)+".fits"
					if os.path.isfile(filename_initial):
						l = l + 1
					else:
						break

				while True:
					ret, frame = video_capture.read()
					if (ret != True):
						break

					filename_tif = folder+"/"+str(imtype)+"_"+str(l)+"_"+str(n)+".tif"
					filename_fits = filename_tif.replace(".tif", ".fits")

					cv2.imwrite(filename_tif, frame)

					hdu = fits.PrimaryHDU()
					im = Image.open(filename_tif)
					hdu.data = np.array(im)
					hdu.writeto(filename_fits, overwrite=True)

					if data_list_entry.state == False:
						os.remove(filename_tif)

					n += 1
				l += 1

			else:
				pass

	else:
		print("Data mode error")
		return(2)

	return(0)

max_stars_repo_name: cjduffy/HighRes

max_issues_repo_name: cjduffy/HighRes

max_forks_repo_name: cjduffy/HighRes

reason: import numpy,from astropy


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

# Copyright (c) 2016-2017, Cabral, Juan; Luczywo, Nadia
# All rights reserved.

# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:

# * Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.

# * Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.

# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.

# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

# =============================================================================
# DOCS
# =============================================================================

"""Methods based on an aggregating function representing
“closeness to the ideal”.

"""

__all__ = ['TOPSIS']

# =============================================================================
# IMPORTS
# =============================================================================

import numpy as np

from .. import rank
from ..validate import MIN, MAX
from ..utils.doc_inherit import doc_inherit

from ._dmaker import DecisionMaker

# =============================================================================
# Function
# =============================================================================

def topsis(nmtx, ncriteria, nweights):

    # apply weights
    wmtx = np.multiply(nmtx, nweights)

    # extract mins and maxes
    mins = np.min(wmtx, axis=0)
    maxs = np.max(wmtx, axis=0)

    # create the ideal and the anti ideal arrays
    ideal = np.where(ncriteria == MAX, maxs, mins)
    anti_ideal = np.where(ncriteria == MIN, maxs, mins)

    # calculate distances
    d_better = np.sqrt(np.sum(np.power(wmtx - ideal, 2), axis=1))
    d_worst = np.sqrt(np.sum(np.power(wmtx - anti_ideal, 2), axis=1))

    # relative closeness
    closeness = d_worst / (d_better + d_worst)

    # compute the rank and return the result
    return rank.rankdata(closeness, reverse=True), ideal, anti_ideal, closeness

# =============================================================================
# OO
# =============================================================================

class TOPSIS(DecisionMaker):
    """TOPSIS is based on the concept that the chosen alternative should have
    the shortest geometric distance from the ideal solution
    and the longest euclidean distance from the worst solution.

    An assumption of TOPSIS is that the criteria are monotonically increasing
    or decreasing, and also allow trade-offs between criteria, where a poor
    result in one criterion can be negated by a good result in another
    criterion.

    Parameters
    ----------

    mnorm : string, callable, optional (default="vector")
        Normalization method for the alternative matrix.

    wnorm : string, callable, optional (default="sum")
        Normalization method for the weights array.

    Returns
    -------

    Decision : :py:class:`skcriteria.madm.Decision`
        With values:

        - **kernel_**: None
        - **rank_**: A ranking (start at 1) where the i-nth element represent
          the position of the i-nth alternative.
        - **best_alternative_**: The index of the best alternative.
        - **alpha_solution_**: True
        - **beta_solution_**: False
        - **gamma_solution_**: True
        - **e_**: Particular data created by this method.

          - **e_.closeness**: Array where the i-nth element represent the
            closenees of the i-nth alternative to ideal and worst solution.

    References
    ----------

    .. [1] Yoon, K., & Hwang, C. L. (1981). Multiple attribute decision
       making: methods and applications. SPRINGER-VERLAG BERLIN AN.
    .. [2] TOPSIS. In Wikipedia, The Free Encyclopedia. Retrieved
       from https://en.wikipedia.org/wiki/TOPSIS
    .. [3] Tzeng, G. H., & Huang, J. J. (2011). Multiple attribute decision
       making: methods and applications. CRC press.

    """

    def __init__(self, mnorm="vector", wnorm="sum"):
        super(TOPSIS, self).__init__(mnorm=mnorm, wnorm=wnorm)

    @doc_inherit
    def solve(self, ndata):
        nmtx, ncriteria, nweights = ndata.mtx, ndata.criteria, ndata.weights
        rank, ideal, anti_ideal, closeness = topsis(nmtx, ncriteria, nweights)
        extra = {
            "ideal": ideal,
            "anti_ideal": anti_ideal,
            "closeness": closeness}
        return None, rank, extra

max_stars_repo_name: dand-oss/scikit-criteria

max_issues_repo_name: dand-oss/scikit-criteria

max_forks_repo_name: dand-oss/scikit-criteria

reason: import numpy


text:
import string

from flask import Blueprint
from root.modules.consolidations.dao.consolidation_dao_impl import getFileDetails, getFileList, getProjectDetails, saveConsolidation
import os
from flask import Flask, render_template, url_for, json
import re
import pandas as pd
from flask import request, jsonify
import numpy as np

mod = Blueprint ('consolidations',__name__)

@mod.route('/getConsolidationFile', methods=['GET', 'POST'])
def consolidationList():
  try:
    if request.method == "POST":
        srcGen_ac_token = request.get_json()['srcGen_ac_token']
        user_id = request.get_json()['user_id']

        if user_id == '' and srcGen_ac_token == '':
          return jsonify({'message': 'Mandatory field validation failed!'})
        else:
          getData = getFileDetails(srcGen_ac_token, user_id)
          #fileName =  getData['FILE_PATH'] + getData['FILE_NAME']

          json_url = os.path.join(os.path.abspath(getData['FILE_PATH']), "", getData['FILE_NAME'])
          data = json.load(open(json_url))
          return jsonify({'FILE_CONTENT':data,'FILE_NAME':getData['FILE_NAME'],'FILE_ID':getData['FILE_ID']})

  except Exception as e:
    return jsonify({'error': str(e)})

@mod.route('/fileList', methods=['GET', 'POST'])
def fileList():
  try:
    if request.method == "POST":
        project_ac_token = request.get_json()['project_ac_token']
        user_id = request.get_json()['user_id']
        project_details = getProjectDetails(project_ac_token, user_id)
        project_id = project_details.PKEY_PROJECTS
        if project_id !='' and user_id !='' :
          projectFileList = getFileList(project_id, user_id)
          return jsonify(projectFileList)

  except Exception as e:
    return jsonify({'error': str(e)})

@mod.route('/addMoreFiles', methods=['GET', 'POST'])
def addMoreFiles():
    try:
        if request.method == "POST":
            project_ac_token = request.get_json()['project_ac_token']
            user_id = request.get_json()['user_id']
            files = request.get_json()['files']
            #print(files)

            con_attribute = list()
            rowarray_list = []
            resultarray_list = []
            for value in files:
              getData = getFileDetails(value, user_id)
              json_url = os.path.join(os.path.abspath(getData['FILE_PATH']), "", getData['FILE_NAME'])
              foo = open(json_url)
              data = json.load(foo)
              foo.close()
              #dataKeys = list()

              for i in data[0].keys():
                #dataKeys.append(i)
                if i in con_attribute:
                    y = "test"
                else:
                    con_attribute.append(i)
              #print(dataKeys)
              #t = {'FILE_AC_TOKEN':value,'KEYS':dataKeys,'NO_KEYS':len(dataKeys),'FILE_NAME':getData['FILE_NAME']}
              #rowarray_list.append(t)

            #t= {'CON_ATTR':con_attribute,'NO_CON_ATTR':len(con_attribute), 'DATA':rowarray_list}
            #resultarray_list.append(t)

            for value in files:
              getData = getFileDetails(value, user_id)
              json_url = os.path.join(os.path.abspath(getData['FILE_PATH']), "", getData['FILE_NAME'])
              foo = open(json_url)
              data = json.load(foo)
              foo.close()
              dataKeys = list()
              data_list = list()

              for i in data[0].keys():
                dataKeys.append(i)

              for conval in con_attribute:
                  if conval in dataKeys:
                      data_list.append(conval)
                  else:
                      data_list.append("")
              #print(dataKeys)
              t = {'FILE_AC_TOKEN':value,'KEYS':data_list,'NO_KEYS':len(data_list),'FILE_NAME':getData['FILE_NAME']}
              rowarray_list.append(t)

            t= {'CON_ATTR':con_attribute,'NO_CON_ATTR':len(con_attribute), 'DATA':rowarray_list}
            resultarray_list.append(t)

            return jsonify(resultarray_list)

    except Exception as e:
        return jsonify({'error': str(e)})

@mod.route('/viewPreview', methods=['GET', 'POST'])
def viewPreview():
    try:
        if request.method == "POST":
            project_ac_token = request.get_json()['project_ac_token']
            user_id = request.get_json()['user_id']
            actual_file = request.get_json()['actual_file']
            files = request.get_json()['files']
            #con_attributes = request.get_json()['con_attributes']

            getData = getFileDetails(actual_file, user_id)
            json_url = os.path.join(os.path.abspath(getData['FILE_PATH']), "", getData['FILE_NAME'])
            foo = open(json_url)
            file1 = json.load(foo)
            foo.close()

            column1 = list()
            result_data = list();

            ###   get columns from 1st Array
            for i in file1[0].keys():
                column1.append(i)

            for value in files:
                if actual_file!=value:
                    column2 = list()
                    common_col = list()

                    ###   Read 2nd File
                    getData2 = getFileDetails(value, user_id)
                    json_url2 = os.path.join(os.path.abspath(getData2['FILE_PATH']), "", getData2['FILE_NAME'])
                    foo2 = open(json_url2)
                    file2 = json.load(foo2)

                    foo2.close()

                    ###   get columns from 2nd Array
                    for i in file2[0].keys():
                        column2.append(i)
                    new_array = np.intersect1d(column1, column2)

                    for i in new_array:
                        common_col.append(i)
    #                tuple_common_col = tuple(common_col)

                    ###   Merge Both array
                    A = pd.DataFrame(file1)
                    B = pd.DataFrame(file2)

    #                frames = [A, B]
    #                merge_data = pd.concat(frames,str(common_col))
                    print(common_col)
                    merge_data = pd.DataFrame.merge(A, B, how='inner', on=common_col)
                    column1 = common_col
                    print(merge_data)
                    result_data = pd.DataFrame.to_json(merge_data, orient='records',lines=False)
                    file1 = result_data

            return result_data

    except Exception as e:
        return jsonify({'error': str(e)})

@mod.route('/savePreview', methods=['GET', 'POST'])
def savePreview():
    try:
        if request.method == "POST":
            project_ac_token = request.get_json()['project_ac_token']
            user_id = request.get_json()['user_id']
            actual_file = request.get_json()['actual_file']
            files = request.get_json()['files']
            con_file_name = request.get_json()['con_file_name']
            result_data = saveConsolidation(project_ac_token,user_id,actual_file,files,con_file_name)

            return result_data

    except Exception as e:
        return jsonify({'error': str(e)})

max_stars_repo_name: itcodertech/freedom

max_issues_repo_name: itcodertech/freedom

max_forks_repo_name: itcodertech/freedom


text:
import pickle
import numpy as np
import pandas as pd
pred_file = "/home/mhxia/whou/workspace/my_models/adversarial_attack/0/checkpoint-best/predict_results.txt"
preds=pickle.load(open(pred_file,'rb'))
preds = np.argmax(preds, axis=1)
result_file="./result.csv"
outf = open(result_file,'w')
df = pd.DataFrame({'pred':preds})
df.to_csv(result_file, header= False, index=True, sep='\t')


text:
from sympy.holonomic.recurrence import RecurrenceOperators, RecurrenceOperator
from sympy import symbols, QQ

def test_RecurrenceOperator():
    n = symbols("n", integer=True)
    R, Sn = RecurrenceOperators(QQ.old_poly_ring(n), "Sn")
    assert Sn * n == (n + 1) * Sn
    assert Sn * n ** 2 == (n ** 2 + 1 + 2 * n) * Sn
    assert Sn ** 2 * n ** 2 == (n ** 2 + 4 * n + 4) * Sn ** 2
    p = (Sn ** 3 * n ** 2 + Sn * n) ** 2
    q = (
        (n ** 2 + 3 * n + 2) * Sn ** 2
        + (2 * n ** 3 + 19 * n ** 2 + 57 * n + 52) * Sn ** 4
        + (n ** 4 + 18 * n ** 3 + 117 * n ** 2 + 324 * n + 324) * Sn ** 6
    )
    assert p == q

def test_RecurrenceOperatorEqPoly():
    n = symbols("n", integer=True)
    R, Sn = RecurrenceOperators(QQ.old_poly_ring(n), "Sn")
    rr = RecurrenceOperator([n ** 2, 0, 0], R)
    rr2 = RecurrenceOperator([n ** 2, 1, n], R)
    assert not rr == rr2

    # polynomial comparison issue, see https://github.com/sympy/sympy/pull/15799
    # should work once that is solved
    # d = rr.listofpoly[0]
    # assert rr == d

    d2 = rr2.listofpoly[0]
    assert not rr2 == d2

max_stars_repo_name: msgoff/sympy

max_issues_repo_name: msgoff/sympy

max_forks_repo_name: msgoff/sympy

reason: from sympy


text:
""" wrapper for gaussian mixture models of sklearn"""

# author: Eugene J. Ragasa <eragasa@ufl.edu>
# License: MIT

import os

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

from sklearn.mixture import GaussianMixture

import pypospack.utils
from pypospack.pyposmat.data import PyposmatDataFile
from pypospack.pyposmat.data import PyposmatConfigurationFile
from pypospack.pyposmat.analysis import AbstractAnalysis

class GmmClusterAnalysis(AbstractAnalysis):

    def __init__(self,
                 pyposmat_configuration,
                 pyposmat_data,
                 names,
                 n_components):
        AbstractAnalysisClass.__init__(
                self,
                pyposmat_configuration=pyposmat_configuration,
                pyposmat_data=pyposmat_data,
                names=names)

        assert isinstance(n_components,int)
        assert n_components > 0
        self.n_components = n_components
        self.initialize_model(n_components=n_components)

    def initialize_model(self,n_components):
        self.model = GaussianMixture(n_components=n_components,
                                     covariance_type='full',
                                     random_state=0)
        self.model.fit(self.data.df[self.names])
        self.data.df['cluster_id'] = self.model.predict(self.data.df[self.names])

def plot_ellipse(position,covariance,ax=None,**kwargs):
    if ax is None:
        fig, ax = plt.subplots(1,1)

    if covariance.shape == (2,2):
        U, s, Vt = np.linalg.svd(covariance)
        angle = np.degrees(np.arctan2(U[1,0],U[0,0]))
        width, height = 2 * np.sqrt(s)
    else:
        angle = 0
        width, height = 2 * np.sqrt(covariance)

    # draw ellipse

    for nsig in range(1,4):
        ax.add_patch(Ellipse(position,nsig*width,nsig*height,angle,**kwargs))

def plot_gmm(gmm,X,label=True,ax=None,dpi=1200,filename=None):
    if ax is None:
        fig, ax = plt.subplots(1,1)

    cluster_id = gmm.fit(X).predict(X)

    if isinstance(X,np.ndarray):
        x = X[:,0]
        y = X[:,1]
    elif isinstance(X,pd.DataFrame):
        x = X[X.columns[0]]
        y = X[X.columns[1]]

    if label:
        ax.scatter(x,y,c=cluster_id,s=1,cmap='viridis',zorder=2)
    else:
        ax.scatter(x,y,s=40,zorder=2)

    w_factor = 0.2 / gmm.weights_.max()
    for pos, covar, w in zip(gmm.means_,gmm.covariances_,gmm.weights_):
        plot_ellipse(pos,covar,alpha=w*w_factor,ax=ax)

    ax.set_xlabel(X.columns[0])
    ax.set_ylabel(X.columns[1])

    if filename is None:
        plt.show()
    else:
        fig.savefig(filename,dpi=dpi)

class GmmInformationCriteriaAnalysis(AbstractAnalysis):

    def __init__(self,
                 pyposmat_configuration,
                 pyposmat_data,
                 names,
                 max_components=20):

        AbstractAnalysisClass.__init__(
                pyposmat_configuration=pyposmat_configuration,
                pyposmat_data=pyposmat-data,
                names=names)

        assert isinstance(max_components,int)
        assert max_components > 0
        self.max_components = max_components

        self.information_criteria = None

    def do_analysis(self):
        n_components = np.arange(1,max_components)
        self.models = [
                GaussianMixture(n_components=n,covariance_type='full',
                                random_state=0).fit(
                                    self.data.df[self.names]) \
                for n in n_components
                ]

        aic_criteria = [m.aic(self.data.df[self.names]) for m in self.models]
        bic_criteria = [m.bic(self.data.df[self.names]) for m in self.models]
        self.information_criteria = pd.DataFrame(data=OrderedDict([
                ('n_components',n_components),
                ('aic',aic_criteria)
                ('bic',bic_criteria)
                ]))

    def write(self,filename='aic_bic.csv'):
        with open(filename,'w') as f:
            f.write(self.information.criteria.to_csv)

    def plot(self,filename,dpi=1200):
        plt.close("all")

class GmmDataFile():
    pass

def gmm_aic_bic_analysis(
        config_fn,
        data_fn,
        names,
        output_directory='gmm_analysis',
        max_components=20):
    assert isinstance(config_fn,str)
    assert isinstance(data_fn,str)
    assert os.path.isfile(config_fn)
    assert os.path.isfile(data_fn)

    o_config = PyposmatConfigurationFile()
    o_config.read(filename=config_fn)

    o_data = PyposmatDataFile()
    o_data.read(filename=data_fn)
    o_data.create_normalized_errors(
            normalize_type='by_qoi_target',
            qoi_targets=o_config.qoi_targets)
    o_data.df['score'] = o_data.df[o_config.normalized_error_names].abs().sum(axis=1)

    data = o_data.df[names]

    n_components = np.arange(1,max_components)
    models = [GaussianMixture(n_components=n,covariance_type='full',random_state=0).fit(data) for n in n_components]

    # AIC analysis
    aic, aic_idx = min((val,idx) for (idx,val) in enumerate([m.aic(data) for m in models]))
    aic_n_components = n_components[aic_idx]
    aic_criteria = [m.aic(data) for m in models]
    # BIC analysis
    bic, bic_idx = min((val,idx) for (idx,val) in enumerate([m.bic(data) for m in models]))
    bic_n_components = n_components[bic_idx]
    bic_criteria = [m.bic(data) for m in models]

    #plot the criteria
    print('bic_n_components:{}'.format(bic_n_components))
    print('aic_n_components:{}'.format(aic_n_components))
    if not os.path.exists(output_directory):
        os.mkdir(output_directory)
    plot_fn=os.path.join(output_directory,'aic_bic_plot.eps')
    plot_gmm_aic_bic(
            filename=plot_fn,
            n_components=n_components,
            aic_criteria=aic_criteria,
            bic_criteria=bic_criteria,
            aic_n_components=aic_n_components,
            bic_n_components=bic_n_components)
    return aic_n_components, bic_n_components

max_stars_repo_name: eragasa/pypospack

max_issues_repo_name: eragasa/pypospack

max_forks_repo_name: eragasa/pypospack


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

class ParticleApproximation:
    def __init__(self, num_particles, prior, init=True):
        self.num_particles = num_particles
        self.prior = prior

        # create an initial Monte Carlo approximation of the prior
        if init:
            self.particles = prior.rvs(size=num_particles)
            self.weights = np.full(num_particles, 1.0/num_particles)

    @staticmethod
    def load(filename, prior):
        """
        load(filename, prior)

        Loads a previously saved particle approximation (see ParticleApproximation.save).

        Parameters
        ----------
        filename : string
            Name of the file in which the particle approximation was previously saved.
        prior
            Prior distribution used for this particle approximation.
        Returns
        -------
        load(filename, prior) : ParticleApproximation
            The particle approximation that was previously saved into the file.
        """
        data = np.load(filename)
        approximation = ParticleApproximation(data['particles'].size, prior, init=False)
        approximation.particles = data['particles']
        approximation.weights = data['weights']
        return approximation

    def save(self, filename):
        """
        save(filename)

        Saves the particle and weights of the particle approximation in a file for later use.

        Parameters
        ----------
        filename : string
            Name to use for the file in which to save to particle approximation.
        """
        np.savez(filename, particles=self.particles, weights=self.weights)

    def hist(self, **kwargs):
        """
        hist(**kwargs)

        Plots a histrogram of the approximated posterior distribution.

        Parameters
        ----------
        kwargs
            Additional keyword arguments to pass to Matplotlib.hist when plotting the approximation.
        """
        plt.hist(self.particles, weights=self.weights, **kwargs)

    def integrate(self, f):
        """
        integrate(f)

        Estimates the integral of the function with respect to the particle approximation.

        Parameters
        ----------
        f : (float) -> float
            Function to integrate.

        Returns
        -------
        integrate(f) : float
            Integral of this function with respect to the particle approximation.
        """
        fv = np.vectorize(f)
        return np.dot(fv(self.particles), self.weights)

    def sample(self, size):
        """
        sample(int) : np.ndarray

        Generate a pseudo random sample from the approximated distribution.

        Parameters
        ----------
        size : int
            Size of the generated sample.

        Returns
        -------
        sample(size) : np.ndarray
            A sample generated from the approximated distribution of size `size`.
        """
        return np.random.choice(self.particles, size=size, p=self.weights, replace=True)

    def resample(self):
        """
        resample()

        Resamples the set of particles from the current approximation of the target distribution. This
        will also reset the weighting of the particles to an uniform weighting.
        """
        self.particles = np.random.choice(self.particles, size=self.num_particles, p=self.weights, replace=True)
        self.weights = np.full(self.num_particles, 1/self.num_particles)

    def effective_sample_size(self):
        """
        effective_sample_size() : float

        Computes the effective sample size to approximate the effective number of samples of this
        approximation based of the deviation of the approximate weights variance.

        Returns
        -------
        effective_sample_size : float
            The effective sample size of the current particle approximation
        """
        return 1 / np.dot(self.weights, self.weights)

    def reweight(self, importance_potential, target_potential):
        """
        reweight(importance_potential, target_potential)

        Given the potential functions of a importance and target distribution, and assuming that the
        current ParticleApproximation approximates the importance distribution, this function updates
        the weights of the particles to approximate the target distribution.

        Parameters
        ----------
        importance_potential : (np.ndarray) -> np.ndarray
            Potential function of the importance distribution.
        target_potential : (np.ndarray) -> np.ndarray
            Potential function of the target distribution.
        """
        # Update is done in log-scale
        self.weights = np.log(self.weights)

        # Compute log-importance-weights and update current weights
        importance_weights = target_potential(self.particles) - importance_potential(self.particles)
        self.weights += importance_weights

        # Return to linear-scale to normalize weights
        self.weights = np.exp(self.weights)
        self.weights /= self.weights.sum()

    def mh_correction(self, target_potential, proposal_kernel, n_steps):
        """
        mh_correction(target_potential, proposal_kernel, n_steps) : float

        Improve the particle approximation by sampling several steps of the Metropolis-Hastings
        Markov Chain constructed using a given proposal kernel and the with the given target
        distribution as limiting distribution.

        Parameters
        ----------
        target_potential: (np.ndarray) -> np.ndarray
            Potential function of the target distribution.
        proposal_kernel: (np.ndarray) -> np.ndarray
            Function to sample proposals for the Metropolis-Hastings algorithm, conditioned
            on the current particles, passed as a parameter.
        n_steps : int, optional
            Number of steps of sampling from the Markov Chain.

        Returns
        -------
        mh_correction : float
            The average acceptance rate over the `n_steps`.
        """
        total_accepted = 0

        # Sample from the proposal kernel, conditioned on currect particles
        proposals = proposal_kernel(self.particles)

        proposal_potentials = target_potential(proposals)
        current_potentials = target_potential(self.particles)

        for i in range(n_steps):
            # Compute the log acceptance ratio
            potential_ratio = proposal_potentials - current_potentials
            prior_ratio = self.prior.logpdf(proposals) - self.prior.logpdf(self.particles)
            acceptance_ratio = np.exp(potential_ratio + prior_ratio)

            # Randomly accept the transitions based on the log acceptance ratio
            accepted = np.random.uniform(size=self.num_particles) < acceptance_ratio
            self.particles[accepted] = proposals[accepted]

            total_accepted += np.sum(accepted)

            # Recompute necessary potentials for next step
            if i < n_steps - 1:
                # Update current potentials
                current_potentials[accepted] = proposal_potentials[accepted]

                # Sample new proposals
                proposals = proposal_kernel(self.particles)
                proposal_potentials = target_potential(proposals)

        return total_accepted / (n_steps * self.num_particles)

    def smc_update(self, importance_potential, target_potential, proposal_kernel, correction_steps, ess_ratio):
        """
        smc_update(importance_potential, target_potential, proposal_kernel, correction_steps, ess_ratio)

        Updates the current particle approximation using the SMC algorithm with Metropolis-Hastings MCMC correction
        steps.

        Parameters
        ----------
        importance_potential : (np.ndarray) -> np.ndarray
            Potential function of the importance distribution.
        target_potential: (np.ndarray) -> np.ndarray
            Potential function of the target distribution.
        proposal_kernel: (np.ndarray) -> np.ndarray
            Function to sample proposals for the Metropolis-Hastings algorithm, conditioned
            on the current particles, passed as a parameter.
        correction_steps : int, optional
            Number of steps of sampling from the Markov Chain.
        ess_ratio: int, optional
            Specifies the lower bound on the effective sample size. A resampling step
            is then performed if `ESS < self.num_particles/ess_ratio`.

        Returns
        -------
        smc_update : (float, float)
            (Acceptance ratio of the correction steps, ESS after reweighting)
        """
        self.reweight(importance_potential, target_potential)
        acceptance_ratio = self.mh_correction(target_potential, proposal_kernel, n_steps=correction_steps)

        ess = self.effective_sample_size()
        if ess < self.num_particles/ess_ratio:
            self.resample()

        return (acceptance_ratio, ess)

max_stars_repo_name: BayesianLearning/PenduSMC

max_issues_repo_name: BayesianLearning/PenduSMC

max_forks_repo_name: BayesianLearning/PenduSMC


text:
import math

import numpy as np
import zarr
from netCDF4 import Dataset
from src.concurrency import threaded
from concurrent.futures import ThreadPoolExecutor
from concurrent.futures import ProcessPoolExecutor

# Convert NetCDF files to Zarr store
def netcdf_to_zarr(datasets, store, append_axis):
    root = zarr.group(store=store, overwrite = True)
    for i, ds in enumerate(datasets):
        if i == 0:
            __set_meta(ds,root)
            __set_dims(ds,root)
            __set_vars(ds,root)
        else:
            __append_vars(ds, root, append_axis)

# Convert non-json-encodable types to built-in types
def __json_encode(val):
    if isinstance(val, np.integer):
        return int(val)
    elif isinstance(val, np.floating):
        return float(val)
    elif isinstance(val, np.ndarray):
        return val.tolist()
    else:
        return val

# Return attributes as dict
def __dsattrs(dataset):
        # JSON encode attributes so they can be serialized
        return {key: __json_encode(getattr(dataset, key)) for key in dataset.ncattrs() }

def __get_dataset(ds):
    return Dataset(ds)

# Set file metadata
def __set_meta(ds, group):
    print("Set meta")
    group.attrs.put(__dsattrs(Dataset(ds)));

def __set_dim(ds, group, name):
    print("Set dim")
    dataset = Dataset(ds)
    dim = dataset.dimensions[name]
    group.create_dataset(name, \
        data=np.arange(dim.size), \
        shape=(dim.size,), \
        chunks=(1<<16,) if dim.isunlimited() else (dim.size,), \
        dtype=np.int32 \
    )
    # Set dimension attrs
    group[name].attrs['_ARRAY_DIMENSIONS'] = [name]

# Set dimensions
def __set_dims(ds, group):
    dataset = Dataset(ds)
    with ProcessPoolExecutor(max_workers=8) as executor:
        for name, dim in dataset.dimensions.items():
            executor.submit(__set_dim, ds, group, name)

# Calculate chunk size for variable
def __get_var_chunks(var, max_size):
    chunks = []
    # TODO: Improve chunk size calculation
    for i, dim in enumerate(var.shape):
        dim_chunk_length = min(math.floor(max_size ** (1/(len(var.shape)-i))), dim)
        max_size //= dim_chunk_length
        chunks.append(dim_chunk_length)

    return tuple(chunks)

# Set variable data, including dimensions and metadata

def __set_var(ds, group, name):
    print("Setting " + name)
    dataset = Dataset(ds)
    var = dataset.variables[name]
    group.create_dataset(name, \
        data=var, \
        shape=var.shape, \
        chunks=(__get_var_chunks(var, 2<<24)), \
        dtype=var.dtype \
    )
    attrs = __dsattrs(var)
    attrs['_ARRAY_DIMENSIONS'] = list(var.dimensions)
    group[name].attrs.put(attrs);

def __set_vars(ds, group):
    dataset = Dataset(ds)
    print("Set vars")
    with ProcessPoolExecutor(max_workers=8) as executor:
        for name, var in dataset.variables.items():
            executor.submit(__set_var, ds, group, name)

# Append data to existing variable

def __append_var(ds, group, name, dim_name):
    print("Appending " + name)
    dataset = Dataset(ds)
    var = dataset.variables[name]
    if dim_name in var.dimensions:
        axis = group[name].attrs['_ARRAY_DIMENSIONS'].index(dim_name)
        group[name].append(var, axis)

def __append_vars(ds, group, dim):
    print("Append vars")
    dataset = Dataset(ds)
    group[dim].append(np.arange(group[dim].size, group[dim].size + dataset.dimensions[dim].size))

    with ProcessPoolExecutor(max_workers=8) as executor:
        for name, var in dataset.variables.items():
            executor.submit(__append_var, ds, group, name, dim)


text:
from sklearn.cluster import KMeans
from sklearn.metrics.cluster import normalized_mutual_info_score
import numpy as np

def evaluation(X, Y, Kset):
    print('----test------')
    print(X)
    print('-------')
    print(Y)
    print('---------------')
    num = X.shape[0]
    classN = np.max(Y)+1
    kmax = np.max(Kset)
    recallK = np.zeros(len(Kset))
    #compute NMI
    kmeans = KMeans(n_clusters=classN).fit(X)
    nmi = normalized_mutual_info_score(Y, kmeans.labels_, average_method='arithmetic')
    #compute Recall@K
    sim = X.dot(X.T)
    minval = np.min(sim) - 1.
    sim -= np.diag(np.diag(sim))
    sim += np.diag(np.ones(num) * minval)
    indices = np.argsort(-sim, axis=1)[:, : kmax]
    YNN = Y[indices]
    for i in range(0, len(Kset)):
        pos = 0.
        for j in range(0, num):
            if Y[j] in YNN[j, :Kset[i]]:
                pos += 1.
        recallK[i] = pos/num
    return nmi, recallK

max_stars_repo_name: cameronchoi/SoftTriple

max_issues_repo_name: cameronchoi/SoftTriple

max_forks_repo_name: cameronchoi/SoftTriple


text:
import argparse

import numpy as np
import torch
from model import NetworkEvaluator
from model import TCN
from sklearn.model_selection import RandomizedSearchCV
from utils import data_generator
from utils import make_experiment_dir
from utils import plot_and_save
from utils import plot_and_save_prediction
from utils import save_args

parser = argparse.ArgumentParser(description='Sequence Modeling - Simply Cozy')
parser.add_argument('--cuda', action='store_false',
                    help='use CUDA (default: True)')
parser.add_argument('--dropout', type=float, default=0.25,
                    help='dropout applied to layers (default: 0.25)')
parser.add_argument('--clip', type=float, default=0.2,
                    help='gradient clip, -1 means no clip (default: 0.2)')
parser.add_argument('--epochs', type=int, default=100,
                    help='upper epoch limit (default: 100)')
parser.add_argument('--ksize', type=int, default=5,
                    help='kernel size (default: 5)')
parser.add_argument('--levels', type=int, default=4,
                    help='# of levels (default: 4)')
parser.add_argument('--log-interval', type=int, default=1, metavar='N',
                    help='report interval (default: 1')
parser.add_argument('--lr', type=float, default=1e-3,
                    help='initial learning rate (default: 1e-3)')
parser.add_argument('--lr_adapt', type=str, default='valid',
                    help='learn rate adaptation (default: valid)')
parser.add_argument('--optim', type=str, default='Adam',
                    help='optimizer to use (default: Adam)')
parser.add_argument('--nhid', type=int, default=150,
                    help='number of hidden units per layer (default: 150)')
parser.add_argument('--series_x', type=str, default='10000000000000000000000000',
                    help='the series of the dataset to use for training (default: 10000000000000000000000000)')
parser.add_argument('--series_y', type=str, default='10000000000000000000000000',
                    help='the series of the dataset to predict (default: 10000000000000000000000000)')
parser.add_argument('--seed', type=int, default=-1,
                    help='random seed (default: -1 = random seed)')
parser.add_argument('--validseqlen', type=int, default=1,
                    help='valid sequence length (# of Samples for Loss computation) (default: 1)')
parser.add_argument('--seqstepwidth', type=int, default=1,
                    help='stepwidth for batch drawing (default: 1)')
parser.add_argument('--seq_len', type=int, default=100,
                    help='total sequence length (default: 200)')
parser.add_argument('--batch_size', type=int, default=12, metavar='N',
                    help='batch size (default: 12)')
parser.add_argument('--documentation', action='store_false',
                    help='enable documentation of experiments (default: True)')
parser.add_argument('--pca', action='store_true',
                    help='use PCA on exegenous series (default: False)')
parser.add_argument('--dataset_path', type=str, default=r'./data/Synthetic/synthetic_data.p',
                    help='path to the dataset (default: ./data/Synthetic/synthetic_data.p)')

args = parser.parse_args()
print(args)

dataset_config = {
    'test_size': 0.2,
    'sampling_size': 1,
    'data_file_name': args.dataset_path,
    'sample_frequency': 1,
    'max_data': 9999999,
    'target_variables': ['Soll'],
    'seed': None if args.seed == -1 else args.seed,
    'train_valid_split': [0.8, 0.2],
    'series_x': args.series_x,
    'series_y': args.series_y,
    'pca': args.pca,
    'batch_size': args.batch_size,
    'cuda': args.cuda,
}

test_hyperparams = {
    "dropout": [-1] + list(np.linspace(0.05, 0.6, 12)),
    "clip": [-1] + list(np.linspace(0.05, 0.6, 12)),
    "ksize": [2, 4, 6, 8, 10],
    "levels": [2, 4, 6, 8],
    "lr": [1e-2, 1e-3, 1e-4, 1e-5, 1e-6],
    "nhid": [50, 150, 200, 250],
    "batch_size": [16, 32, 64, 96],
    "seq_len": [100, 300, 500, 700],
    "validseqlen": [1]
}

def run(optimize=False):
    if args.seed != -1:
        torch.manual_seed(args.seed)
        torch.cuda.manual_seed(args.seed)
    else:
        torch.cuda.seed_all()
    if torch.cuda.is_available():
        if not args.cuda:
            print("WARNING: You have a CUDA device, so you should probably run with --cuda")

    data, pca_scaler = data_generator(dataset_config)
    data.x_valid = data.x_valid[0]
    data.y_valid = data.y_valid[0]
    if args.cuda:
        data.x_train = data.x_train.cuda()
        data.x_test = data.x_test.cuda()
        data.x_valid = data.x_valid.cuda()
        data.y_train = data.y_train.cuda()
        data.y_test = data.y_test.cuda()
        data.y_valid = data.y_valid.cuda()
    input_size = data.x_train.size(1)
    output_size = data.y_train.size(1)

    if optimize:
        optimizations_path = "./optimizations/"
        optimization_experiment_folder = make_experiment_dir(optimizations_path)
        save_args(optimization_experiment_folder, args)
        print("Saving best Random Search Field...")
        with open(optimization_experiment_folder + "random_search_field.txt", "w") as file:
            for key, value in test_hyperparams.items():
                line = str(key) + "     " + str(value) + "\n"
                file.write(line)
        file.close()
        print('Saved as %s' % optimization_experiment_folder + "random_search_field.txt")
        rs = RandomizedSearchCV(NetworkEvaluator(input_size, output_size, args.cuda), test_hyperparams, n_jobs=1)
        rs.fit(data.x_train.cpu().numpy(), y=data.y_train.cpu().numpy())
        best_estimator = rs.best_estimator_
        print("Best Params: \n")
        print(rs.best_params_)
        print("Saving best Parameters...")
        with open(optimization_experiment_folder + "best_params.txt", "w") as file:
            for key, value in rs.best_params_.items():
                line = str(key) + "     " + str(value) + "\n"
                file.write(line)
        file.close()
        print('Saved as %s' % optimization_experiment_folder + "best_params.txt")
    else:
        experiment_path = "./experiments/"
        experiment_folder = make_experiment_dir(experiment_path)

        model = TCN(input_size, output_size, cuda=args.cuda, ksize=args.ksize, dropout=args.dropout, clip=args.clip,
                    epochs=args.epochs, levels=args.levels, log_interval=args.log_interval, lr=args.lr,
                    optim=args.optim, nhid=args.nhid, validseqlen=args.validseqlen, seqstepwidth=args.seqstepwidth,
                    seq_len=args.seq_len, batch_size=args.batch_size)
        all_epoch_train_losses, all_epoch_valid_losses, all_epoch_test_losses = model.fit_and_test(x_train=data.x_train,
                                                                                y_train=data.y_train,
                                                                                x_test=data.x_test,
                                                                                y_test=data.y_test,
                                                                                x_valid=data.x_valid,
                                                                                y_valid=data.y_valid,
                                                                                lr_adapt=args.lr_adapt,
                                                                                experiment_folder=experiment_folder)
        save_args(experiment_folder, args)
        plot_and_save_prediction(model, data.x_test, data.y_test, args, path=experiment_folder)
        plot_and_save([all_epoch_train_losses, all_epoch_valid_losses, all_epoch_test_losses],
                      ['Train Loss over Epochs', 'Validation Loss over Epochs', 'Test Loss over Epochs'],
                      experiment_folder + 'losses.svg',
                      x_names=['Epoch', 'Epoch', 'Epoch'])

if __name__ == "__main__":
    run(False)

max_stars_repo_name: rmiddelanis/master-thesis

max_issues_repo_name: rmiddelanis/master-thesis

max_forks_repo_name: rmiddelanis/master-thesis


text:
import numpy as np
from scipy import sparse as sp

def check_clusterings(labels_true, labels_pred):
    """Check that the two clusterings matching 1D integer arrays."""
    labels_true = np.asarray(labels_true)
    labels_pred = np.asarray(labels_pred)

    # input checks
    if labels_true.ndim != 1:
        raise ValueError(
            "labels_true must be 1D: shape is %r" % (labels_true.shape,))
    if labels_pred.ndim != 1:
        raise ValueError(
            "labels_pred must be 1D: shape is %r" % (labels_pred.shape,))
    if labels_true.shape != labels_pred.shape:
        raise ValueError(
            "labels_true and labels_pred must have same size, got %d and %d"
            % (labels_true.shape[0], labels_pred.shape[0]))
    return labels_true, labels_pred

def contingency_matrix(labels_true, labels_pred, eps=None, sparse=False):
    if eps is not None and sparse:
        raise ValueError("Cannot set 'eps' when sparse=True")

    classes, class_idx = np.unique(labels_true, return_inverse=True)
    clusters, cluster_idx = np.unique(labels_pred, return_inverse=True)
    n_classes = classes.shape[0]
    n_clusters = clusters.shape[0]
    # Using coo_matrix to accelerate simple histogram calculation,
    # i.e. bins are consecutive integers
    # Currently, coo_matrix is faster than histogram2d for simple cases
    contingency = sp.coo_matrix((np.ones(class_idx.shape[0]),
                                 (class_idx, cluster_idx)),
                                shape=(n_classes, n_clusters),
                                dtype=np.int)
    if sparse:
        contingency = contingency.tocsr()
        contingency.sum_duplicates()
    else:
        contingency = contingency.toarray()
        if eps is not None:
            # don't use += as contingency is integer
            contingency = contingency + eps
    return contingency

def fowlkes_mallows_score(labels_true, labels_pred, sparse=False):
    labels_true, labels_pred = check_clusterings(labels_true, labels_pred)
    n_samples, = labels_true.shape

    c = contingency_matrix(labels_true, labels_pred, sparse=True)
    tk = np.dot(c.data, c.data) - n_samples
    pk = np.sum(np.asarray(c.sum(axis=0)).ravel() ** 2) - n_samples
    qk = np.sum(np.asarray(c.sum(axis=1)).ravel() ** 2) - n_samples
    return tk / pk, tk / qk, tk / np.sqrt(pk * qk) if tk != 0. else 0.

reason: import numpy,from scipy


text:
from tensorflow.python.platform import flags
from tensorflow.python.framework import ops
import tensorflow as tf
import tensorflow.contrib.slim as slim
from tensorflow.contrib import layers
from tensorflow.contrib.layers.python.layers import regularizers

from tensorflow.python.ops import array_ops
from tensorflow.python.ops import math_ops
from tensorflow.python.ops import nn_ops
from tensorflow.python.ops import variable_scope
from tensorflow.python.util import nest
from tensorflow.python.ops import variable_scope as vs
from tensorflow.python.ops import init_ops
from tensorflow.python.platform import flags

import numpy as np
_WEIGHTS_VARIABLE_NAME = "kernel"
_BIAS_VARIABLE_NAME = "bias"

Flags = flags.FLAGS

def orthogonal_initializer(shape, dtype=tf.float32, *args, **kwargs):
  del args
  del kwargs
  flat_shape = (shape[0], np.prod(shape[1:]))
  w = np.random.randn(*flat_shape)
  u, _, v = np.linalg.svd(w, full_matrices=False)
  w = u if u.shape == flat_shape else v
  return tf.constant(w.reshape(shape), dtype=dtype)

def linear(args,
           output_size,
           bias,
           bias_initializer=None,
           kernel_initializer=None):
  """Linear map: sum_i(args[i] * W[i]), where W[i] is a variable.

  Args:
    args: a 2D Tensor or a list of 2D, batch x n, Tensors.
    output_size: int, second dimension of W[i].
    bias: boolean, whether to add a bias term or not.
    bias_initializer: starting value to initialize the bias
      (default is all zeros).
    kernel_initializer: starting value to initialize the weight.

  Returns:
    A 2D Tensor with shape [batch x output_size] equal to
    sum_i(args[i] * W[i]), where W[i]s are newly created matrices.

  Raises:
    ValueError: if some of the arguments has unspecified or wrong shape.
  """
  if args is None or (nest.is_sequence(args) and not args):
    raise ValueError("`args` must be specified")
  if not nest.is_sequence(args):
    args = [args]

  # Calculate the total size of arguments on dimension 1.
  total_arg_size = 0
  shapes = [a.get_shape() for a in args]
  for shape in shapes:
    if shape.ndims != 2:
      raise ValueError("linear is expecting 2D arguments: %s" % shapes)
    if shape[1].value is None:
      raise ValueError("linear expects shape[1] to be provided for shape %s, "
                       "but saw %s" % (shape, shape[1]))
    else:
      total_arg_size += shape[1].value

  dtype = [a.dtype for a in args][0]

  # Now the computation.
  scope = vs.get_variable_scope()
  with vs.variable_scope(scope) as outer_scope:
    # linear_regularizer = regularizers.l2_regularizer(Flags.weight_decay)
    weights = vs.get_variable(
        _WEIGHTS_VARIABLE_NAME, [total_arg_size, output_size],
        dtype=dtype,
        initializer=kernel_initializer,)
    if len(args) == 1:
      res = math_ops.matmul(args[0], weights)
    else:
      res = math_ops.matmul(array_ops.concat(args, 1), weights)
    if not bias:
      return res
    with vs.variable_scope(outer_scope) as inner_scope:
      inner_scope.set_partitioner(None)
      if bias_initializer is None:
        bias_initializer = init_ops.constant_initializer(0.0, dtype=dtype)
      biases = vs.get_variable(
          _BIAS_VARIABLE_NAME, [output_size],
          dtype=dtype,
          initializer=bias_initializer,)
    return nn_ops.bias_add(res, biases)


text:
import numpy as np
import os
from keras.models import Model
from keras.layers import Dense, Embedding, Activation, Permute
from keras.layers import Input, Flatten, Dropout
from keras.layers.recurrent import LSTM
from keras.layers.wrappers import TimeDistributed, Bidirectional
from .custom_recurrents import AttentionDecoder

def simpleNMT(pad_length=100,
              n_chars=105,
              n_labels=6,
              embedding_learnable=False,
              encoder_units=256,
              decoder_units=256,
              trainable=True,
              return_probabilities=False):
    """
    Builds a Neural Machine Translator that has alignment attention
    :param pad_length: the size of the input sequence
    :param n_chars: the number of characters in the vocabulary
    :param n_labels: the number of possible labelings for each character
    :param embedding_learnable: decides if the one hot embedding should be refinable.
    :return: keras.models.Model that can be compiled and fit'ed

    *** REFERENCES ***
    Lee, Jason, Kyunghyun Cho, and Thomas Hofmann.
    "Neural Machine Translation By Jointly Learning To Align and Translate"
    """
    input_ = Input(shape=(pad_length,), dtype='float32')
    input_embed = Embedding(n_chars, n_chars,
                            input_length=pad_length,
                            trainable=embedding_learnable,
                            weights=[np.eye(n_chars)],
                            name='OneHot')(input_)

    rnn_encoded = Bidirectional(LSTM(encoder_units, return_sequences=True),
                                name='bidirectional_1',
                                merge_mode='concat',
                                trainable=trainable)(input_embed)

    y_hat = AttentionDecoder(decoder_units,
                             name='attention_decoder_1',
                             output_dim=n_labels,
                             return_probabilities=return_probabilities,
                             trainable=trainable)(rnn_encoded)

    model = Model(inputs=input_, outputs=y_hat)

    return model

if __name__ == '__main__':
    model = simpleNMT()
    model.summary()

max_stars_repo_name: FrancescoZ/Chemception

max_issues_repo_name: FrancescoZ/Chemception

max_forks_repo_name: FrancescoZ/Chemception


text:
"""Classes of random variables."""

import numpy as np
from numpy import sqrt
from struct_tools import NicePrint

from dapper.tools.matrices import CovMat

class RV(NicePrint):
    "Class to represent random variables."

    printopts = NicePrint.printopts.copy()
    printopts["ordering"] = "linenumber"
    printopts["reverse"] = True

    def __init__(self, M, **kwargs):
        """
         - M    <int>     : ndim
         - is0  <bool>    : if True, the random variable is identically 0
         - func <func(N)> : use this sampling function. Example:
                            RV(M=4,func=lambda N: rand(N,4)
         - file <str>     : draw from file. Example:
                            RV(M=4,file=dpr.rc.dirs.data/'tmp.npz')
        The following kwords (versions) are available,
        but should not be used for anything serious
        (use instead subclasses, like GaussRV).
         - icdf <func(x)> : marginal/independent  "inverse transform" sampling.
                            Example: RV(M=4,icdf = scipy.stats.norm.ppf)
         - cdf <func(x)>  : as icdf, but with approximate icdf, from interpolation.
                            Example: RV(M=4,cdf = scipy.stats.norm.cdf)
         - pdf  <func(x)> : "acceptance-rejection" sampling. Not implemented.
        """
        self.M = M
        for key, value in kwargs.items():
            setattr(self, key, value)

    def sample(self, N):
        if getattr(self, 'is0', False):
            # Identically 0
            E = np.zeros((N, self.M))
        elif hasattr(self, 'func'):
            # Provided by function
            E = self.func(N)
        elif hasattr(self, 'file'):
            # Provided by numpy file with sample
            data   = np.load(self.file)
            sample = data['sample']
            N0     = len(sample)
            if 'w' in data:
                w = data['w']
            else:
                w = np.ones(N0)/N0
            idx = np.random.choice(N0, N, replace=True, p=w)
            E   = sample[idx]
        elif hasattr(self, 'icdf'):
            # Independent "inverse transform" sampling
            icdf = np.vectorize(self.icdf)
            uu   = np.random.rand(N, self.M)
            E    = icdf(uu)
        elif hasattr(self, 'cdf'):
            # Like above, but with inv-cdf approximate, from interpolation
            if not hasattr(self, 'icdf_interp'):
                # Define inverse-cdf
                from scipy.interpolate import interp1d
                from scipy.optimize import fsolve
                cdf    = self.cdf
                Left,  = fsolve(lambda x: cdf(x) - 1e-9, 0.1)
                Right, = fsolve(lambda x: cdf(x) - (1-1e-9), 0.1)
                xx     = np.linspace(Left, Right, 1001)
                uu     = np.vectorize(cdf)(xx)
                icdf   = interp1d(uu, xx)
                self.icdf_interp = np.vectorize(icdf)
            uu = np.random.rand(N, self.M)
            E  = self.icdf_interp(uu)
        elif hasattr(self, 'pdf'):
            # "acceptance-rejection" sampling
            raise NotImplementedError
        else:
            raise KeyError
        assert self.M == E.shape[1]
        return E

# TODO 4: improve constructor (treatment of arg cases is too fragile).
class RV_with_mean_and_cov(RV):
    """Generic multivariate random variable characterized by mean and cov.

    This class must be subclassed to provide sample(),
    i.e. its main purpose is provide a common convenience constructor.
    """

    def __init__(self, mu=0, C=0, M=None):
        """Init allowing for shortcut notation."""

        if isinstance(mu, CovMat):
            raise TypeError("Got a covariance paramter as mu. "
                            + "Use kword syntax (C=...) ?")

        # Set mu
        mu = np.atleast_1d(mu)
        assert mu.ndim == 1
        if len(mu) > 1:
            if M is None:
                M = len(mu)
            else:
                assert len(mu) == M
        else:
            if M is not None:
                mu = np.ones(M)*mu

        # Set C
        if isinstance(C, CovMat):
            if M is None:
                M = C.M
        else:
            if np.isscalar(C) and C == 0:
                pass  # Assign as pure 0!
            else:
                if np.isscalar(C):
                    M = len(mu)
                    C = CovMat(C*np.ones(M), 'diag')
                else:
                    C = CovMat(C)
                    if M is None:
                        M = C.M

        # Validation
        if len(mu) not in (1, M):
            raise TypeError("Inconsistent shapes of (M,mu,C)")
        if M is None:
            raise TypeError("Could not deduce the value of M")
        try:
            if M != C.M:
                raise TypeError("Inconsistent shapes of (M,mu,C)")
        except AttributeError:
            pass

        # Assign
        self.M  = M
        self.mu = mu
        self.C  = C

    def sample(self, N):
        """Sample N realizations. Returns N-by-M (ndim) sample matrix.

        Example::

        >>> plt.scatter(*(UniRV(C=randcov(2)).sample(10**4).T))
        """
        if self.C == 0:
            D = np.zeros((N, self.M))
        else:
            D = self._sample(N)
        return self.mu + D

    def _sample(self, N):
        raise NotImplementedError("Must be implemented in subclass")

class GaussRV(RV_with_mean_and_cov):
    """Gaussian (Normal) multivariate random variable."""

    def _sample(self, N):
        R = self.C.Right
        D = np.random.randn(N, len(R)) @ R
        return D

class LaplaceRV(RV_with_mean_and_cov):
    """
    Laplace (double exponential) multivariate random variable.
    This is an elliptical generalization. Ref:
    Eltoft (2006) "On the Multivariate Laplace Distribution".
    """

    def _sample(self, N):
        R = self.C.Right
        z = np.random.exponential(1, N)
        D = np.random.randn(N, len(R))
        D = z[:, None]*D
        return D @ R / sqrt(2)

class LaplaceParallelRV(RV_with_mean_and_cov):
    """
    A NON-elliptical multivariate generalization of
    the Laplace (double exponential) random variable.
    """

    def _sample(self, N):
        # R = self.C.Right   # contour: sheared rectangle
        R = self.C.sym_sqrt  # contour: rotated rectangle
        D = np.random.laplace(0, 1, (N, len(R)))
        return D @ R / sqrt(2)

class StudRV(RV_with_mean_and_cov):
    """
    Student-t multivariate random variable.
    Assumes the covariance exists,
    which requires degreee-of-freedom (dof) > 1+ndim.
    Also requires that dof be integer,
    since chi2 is sampled via Gaussians.
    """

    def __init__(self, dof, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.dof = dof

    def _sample(self, N):
        R = self.C.Right
        nu = self.dof
        r = nu/np.sum(np.random.randn(N, nu)**2, axis=1)  # InvChi2
        D = sqrt(r)[:, None]*np.random.randn(N, len(R))
        return D @ R * sqrt((nu-2)/nu)

class UniRV(RV_with_mean_and_cov):
    """
    Uniform multivariate random variable.
    with an elliptic-shape support.
    Ref: Voelker et al. (2017) "Efficiently sampling
    vectors and coordinates from the n-sphere and n-ball"
    """

    def _sample(self, N):
        R = self.C.Right
        D = np.random.randn(N, len(R))
        r = np.random.rand(N)**(1/len(R)) / np.sqrt(np.sum(D**2, axis=1))
        D = r[:, None]*D
        return D @ R * 2

class UniParallelRV(RV_with_mean_and_cov):
    """
    Uniform multivariate random variable,
    with a parallelogram-shaped support, as determined by the cholesky factor
    applied to the (corners of) the hypercube.
    """

    def _sample(self, N):
        R = self.C.Right
        D = np.random.rand(N, len(R))-0.5
        return D @ R * sqrt(12)

max_stars_repo_name: myying/DAPPER

max_issues_repo_name: myying/DAPPER

max_forks_repo_name: myying/DAPPER

reason: import numpy,from numpy,from scipy


text:
import simcompyl as sim

import holoviews as hv
import numpy as np
import pytest

pytestmark = [pytest.mark.integration,
              pytest.mark.filterwarnings(
                  "ignore:numpy.dtype size changed.*:RuntimeWarning")]

def test_execute(engine):
    from . import walk

    model = walk.ComplexWalk()
    alloc = walk.Simulation() + walk.BasicDistance()

    exec = engine(model, alloc)

    out = exec.run(initial_energy=10)

    limit = alloc.n_samples.value / 4

    assert set(out.columns) == {'x', 'y', 'energy'}
    assert len(out.query('x**2 + y**2 < 10**2')) > limit
    assert len(out.query('energy <= 0')) > limit

    out = exec.run()
    assert len(out.query('energy <= 0')) < limit

def test_nontracing():
    from . import walk

    model = walk.ComplexWalk()
    alloc = walk.Simulation() + walk.BasicDistance()

    exec = sim.engine.NonTracingNumbaEngine(model, alloc)

    out = exec.run(initial_energy=10)

    limit = alloc.n_samples.value / 4

    assert set(out.columns) == {'x', 'y', 'energy'}
    assert len(out.query('x**2 + y**2 < 10**2')) > limit
    assert len(out.query('energy <= 0')) > limit

    out = exec.run()
    assert len(out.query('energy <= 0')) < limit

def test_tracing(engine):
    from . import walk

    model = walk.ComplexWalk()
    alloc = walk.Simulation() + walk.BasicDistance()

    exec = engine(model, alloc)

    tr = model['x', 'y'].take(12)
    with exec.trace(tr, skip=6) as pos:
        exec.run()

    assert len(pos) == (alloc.n_steps.value // 6 + 1) * 12
    assert list(pos.columns) == ['x', 'y']

    qr = (model(dist=lambda x, y: np.sqrt(x ** 2 + y ** 2))
          .quantile([.25, .5, .75]))

    with exec.trace(tr, qr, target=sim.trace.Holotrace) as (pos, dist):
        exec.run()

    assert isinstance(dist, hv.streams.Buffer)
    assert len(dist.data) == (alloc.n_steps.value + 1) * 3
    assert list(dist.data.columns) == ['dist']

max_stars_repo_name: gameduell/simcompyl

max_issues_repo_name: gameduell/simcompyl

max_forks_repo_name: gameduell/simcompyl