text:
import os
import struct
import numpy as np

def load_mnist(path,kind='train',kind1='t10k'):
	labels_path = os.path.join(path,'%s-labels-idx1-ubyte'%kind)
	images_path = os.path.join(path,'%s-images-idx3-ubyte'%kind)
	labels_path1 = os.path.join(path,'%s-labels-idx1-ubyte'%kind1)
	images_path1 = os.path.join(path,'%s-images-idx3-ubyte'%kind1)

	with open(labels_path,'rb') as lbpath:
		magic,n = struct.unpack('>II',
			lbpath.read(8))
		labels = np.fromfile(lbpath,dtype=np.uint8)

	with open(images_path,'rb') as imgpath:
		magic,num,rows,cols =struct.unpack(">IIII",
			imgpath.read(16))
		images = np.fromfile(imgpath,
			dtype=np.uint8).reshape(len(labels),784)
		images = ((images/255.)-.5)*2

	with open(labels_path1,'rb') as lbpath1:
		magic,n = struct.unpack('>II',
			lbpath1.read(8))
		labels1 = np.fromfile(lbpath1,dtype=np.uint8)
	with open(images_path1,'rb') as imgpath1:
		magic,num,rows,cols =struct.unpack(">IIII",
			imgpath1.read(16))
		images1 = np.fromfile(imgpath1,
			dtype=np.uint8).reshape(len(labels1),784)
		images1 = ((images1/255.)-.5)*2

	return images,labels,images1,labels1

X_train, y_train,X_test,y_test = load_mnist('', kind='train',kind1='t10k')
print('Rows: %d, columns: %d'
	% (X_train.shape[0], X_train.shape[1]))
print('Rows: %d, columns: %d'
	% (X_test.shape[0], X_test.shape[1]))

np.savez_compressed('mnist_scaled.npz',
	X_train=X_train,
	y_train=y_train,
	X_test=X_test,
	y_test=y_test)

mnist = np.load('mnist_scaled.npz')
X_test = mnist['X_test']
y_test = mnist['y_test']
print(y_test.shape[0])
print(X_test.shape)


text:
import numpy as np
import scipy.constants as c
import math, cmath
from sympy import *
from math import e
import numba
from numba import jit
import sympy as sp
r, x, a, theta = symbols('r x a theta')
init_printing(use_unicode=True)

@jit(nopython=True, cache=True, parallel=True)
def spherical_to_cartesian(r, theta, phi):
    x = r * np.sin(theta) * np.cos(phi)
    y = r * np.sin(theta) * np.sin(phi)
    z = r * np.cos(theta)
    return x, y, z

@jit(nopython=True, cache=True, parallel=True)
def cartesian_to_spherical(x, y, z):
    r = math.sqrt(x ** 2 + y ** 2 + z ** 2)
    theta = np.arccos(z / math.sqrt(x ** 2 + y ** 2 + z ** 2)) if r != 0 else 0
    if x == 0:
        phi = 0 if (y == 0) else 1.5708
    else:
        phi = np.arctan(y / x)
    return r, theta, phi

@jit(nopython=True, cache=True, parallel=True)
def absolute(number):
    return (number.real**2 + number.imag**2)**0.5

def P_l(l):
    ans = 1 / (2 ** l * factorial(l)) * sp.diff((x ** 2 - 1) ** l, x, int(l))
    return ans

def Pm_l(m, l):
    ans = ((1 - x ** 2) ** (abs(m) / 2)) * sp.diff(P_l(l), x, abs(int(m)))
    return ans

def angular_wave_func(m, l, theta_value, phi):
    if m > 0:
        E = (-1)**m
    else:
        E = 1
    A_factor_1 = (2*l+1)/(4*math.pi)
    A_factor_2 = math.factorial(l-abs(m)) / math.factorial(l+abs(m))
    A = math.sqrt(A_factor_1*A_factor_2)
    B = cmath.exp(m*phi*1j)
    C = (Pm_l(m,l).subs(x,cos(theta)).subs(theta, theta_value)).evalf()
    long_ans = complex(E * A * B * C)
    ans = round(long_ans.real, 5) + round(long_ans.imag, 5)*1j

    return ans

def Lq(q):
    ans = e ** x * sp.diff(e ** -x * x ** q, x, int(q))
    return ans

def assoc_Lq(p, q):
    ans = (-1) ** p * sp.diff(Lq(q), x, int(p))
    return ans

bohr=c.physical_constants['Bohr radius'][0]

def radial_wave_func(n, l, radius):
    if radius == 0:
        return np.nan
    else:
        A_factor_1 = (2/(n*a))**3
        A_factor_2 = (math.factorial(n-l-1)) / (2*n*(math.factorial(n+l))**3)
        A = (A_factor_1 * A_factor_2)**0.5
        B = e**(-radius/(n*a))
        C = ((2*radius)/(n*a))**l
        p = 2*l + 1
        q = n-l-1+p
        D = assoc_Lq(p,q)
        expression = A * B * C * D / (a**(-3/2))
        subbed = expression.subs(x, 2*r/(n*a)).subs(r, radius).subs(a, bohr)
        return round(subbed.evalf(), 5)

#@jit
def linspace(start, stop, num=50):
    increment = (stop-start)/(num-1)
    current = start
    output = []
    for i in range(0,int((stop-start)/increment)+1):
        output.append(round(float(current),5))
        current+=increment
    return output

def meshgrid(x,y,z):
    output = [[],[],[]]
    output_0, output_1, output_2 = [],[],[]
    x_list, y_list, z_list = [],[],[]
    z_inner = []
    for k,z_i in enumerate(z):
        z_inner.append(z_i)
    for i,x_i in enumerate(x):
        x_list.append([x_i,x_i])
        z_list.append(z_inner)
    for j,y_i in enumerate(y):
        y_inner = [[y_i,y_i] for i in range(len(x))]
        y_list.append(y_inner)
        output_0.append(x_list)
        output_2.append(z_list)

    return output_0, y_list, output_2

cartesian_to_spherical_vector = (np.vectorize(cartesian_to_spherical))
angular_wave_vector = (np.vectorize(angular_wave_func))
radial_wave_vector = (np.vectorize(radial_wave_func))
absolute_vector = (np.vectorize(absolute))

def hydrogen_wave_func(n, l, m, roa, Nx, Ny, Nz):
    x_space = np.linspace(-roa, roa, Nx)
    y_space = np.linspace(-roa, roa, Ny)
    z_space = np.linspace(-roa, roa, Nz)

    xx, yy, zz = np.meshgrid(y_space, x_space, z_space)

    r, theta, phi = (cartesian_to_spherical_vector(yy, xx, zz))

    if m == 0:
        angular = angular_wave_vector(m, l, theta, phi)
    elif m < 0:
        angular = (1j / math.sqrt(2)) * (angular_wave_vector(m, l, theta, phi) - (-1) ** m * angular_wave_vector(-m, l, theta, phi))
    elif m > 0:
        angular = (1 / math.sqrt(2)) * (angular_wave_vector(-m, l, theta, phi) + (-1) ** m * angular_wave_vector(m, l, theta, phi))

    radial = radial_wave_vector(n, l, r * a)

    mag = absolute_vector(radial * angular) ** 2

    return np.array(yy), np.array(xx), np.array(zz), np.round(mag, 5), radial * angular

#don't use?
def hydrogen_wave_func_cross(n, l, m, roa, Nx, Ny, Nz):
    x_space = np.linspace(-roa, roa, int(Nx/2))
    y_space = np.linspace(-roa, roa, Ny)
    z_space = np.linspace(-roa, roa, Nz)

    xx, yy, zz = np.meshgrid(y_space, x_space, z_space)

    r, theta, phi = (cartesian_to_spherical_vector(yy, xx, zz))

    if m == 0:
        angular = angular_wave_vector(m, l, theta, phi)
    elif m < 0:
        angular = (1j / math.sqrt(2)) * (angular_wave_vector(m, l, theta, phi) - (-1) ** m * angular_wave_vector(-m, l, theta, phi))
    elif m > 0:
        angular = (1 / math.sqrt(2)) * (angular_wave_vector(-m, l, theta, phi) + (-1) ** m * angular_wave_vector(m, l, theta, phi))

    radial = radial_wave_vector(n, l, r * a)
    mag_true = radial * angular
    mag = absolute_vector(mag_true) ** 2

    return np.array(yy), np.array(xx), np.array(zz), np.round(mag, 5), mag_true

ext: py

max_stars_repo_name: tlkh/atomic-orbitals

max_issues_repo_name: tlkh/atomic-orbitals

max_forks_repo_name: tlkh/atomic-orbitals

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


text:
import chainer
import numpy as np

from onnx_chainer.functions.opset_version import support
from onnx_chainer import onnx_helper

@support((9,))
def convert_SoftmaxCrossEntropy(
        func, opset_version, input_names,
        output_names, context, parameters):
    # obtain input variable
    if not isinstance(func, chainer.FunctionNode):
        raise NotImplementedError(
            'SoftmaxCrossEntropy is currently supported for Chainer>=6.0.0a1.')

    x_var, t_var = func.get_retained_inputs()
    if len(x_var.shape) != 2:
        raise NotImplementedError(
            'ONNX-Chainer currently handles SoftmaxCrossEntropy only when '
            'the dimension of input variable x is exactly two.')
    if np.any(t_var.array == func.ignore_label):
        raise NotImplementedError(
            'ONNX-Chainer currently handles SoftmaxCrossEntropy only when '
            'ignore_label is not used in input variable t.')
    if (not func.normalize) or (func.class_weight is not None) or\
       (func.ignore_label != -1) or (func.reduce != 'mean'):
        raise NotImplementedError(
            'ONNX-Chainer currently handles SoftmaxCrossEntropy only when '
            'argument parameters are default setting.')

    # create intermediate values
    gb = onnx_helper.GraphBuilder()
    x, t = input_names
    y_log = gb.op('LogSoftmax', [x])
    depth = context.add_const(np.array([x_var.shape[1]], dtype=np.int32),
                              'depth')
    zeroone = context.add_const(np.array([0, 1], dtype=x_var.dtype), 'zeroone')
    th = gb.op('OneHot', [t, depth, zeroone])
    s0 = gb.op('Mul', [y_log, th])
    sn = gb.op('Neg', [s0])
    sr = gb.op('ReduceSum', [sn], axes=[1], keepdims=0)
    gb.op_output_named('ReduceMean', [sr], output_names, axes=[0], keepdims=0)

    return gb.nodes()

max_stars_repo_name: disktnk/onnx-chainer

max_issues_repo_name: disktnk/onnx-chainer

max_forks_repo_name: disktnk/onnx-chainer

reason: import numpy


text:
#imports
import scipy.io as sio
import csv
import pandas as pd
import numpy as np

num_iter = 10 #number of times the strategies were run

if __name__ == "__main__":

    dataprep = []

    for iter in range(0, num_iter, 1):

        human_cost_for_failures_atta = 0 #sum of the cost incurred for the tasks the human failed on
        robot_cost_for_failures_atta = 0 #sum of the cost incurred for the tasks the robot failed on
        human_cost_for_failures_tsarouchi = 0 #sum of the cost incurred for the tasks the human failed on and were discarded to the human
        robot_cost_for_failures_tsarouchi = 0 #sum of the cost incurred for the tasks the robot failed on and were discarded to the robot
        human_cost_for_failures_random = 0 #sum of the cost incurred for the tasks the human failed on
        robot_cost_for_failures_random = 0 #sum of the cost incurred for the tasks the robot failed on

        load_file_name = "./results/atta/atta_caseI_eta50_" + str(iter) + ".mat"
        output_atta_mat = sio.loadmat(load_file_name)

        human_successes_from_atta = output_atta_mat["human_successes"] #tasks assigned to the human that were successes
        human_failures_from_atta = output_atta_mat["human_failures"] #tasks assigned to the human that were failures
        robot_successes_from_atta = output_atta_mat["robot_successes"] #tasks assigned to the robot that were successes
        robot_failures_from_atta = output_atta_mat["robot_failures"] #tasks assigned to the robot that were failures
        human_num_tasks_from_atta = output_atta_mat["human_num_tasks"][0][0]
        robot_num_tasks_from_atta = output_atta_mat["robot_num_tasks"][0][0]

        num_successes_from_atta = human_successes_from_atta.shape[1] + robot_successes_from_atta.shape[1]
        num_failures_from_atta = human_failures_from_atta.shape[1] + robot_failures_from_atta.shape[1]
        #human_num_tasks_from_atta
        #robot_num_tasks_from_atta
        #human_successes_from_atta.shape[1]
        #human_failures_from_atta.shape[1]
        #robot_successes_from_atta.shape[1]
        #robot_failures_from_atta.shape[1]

        human_num_failures_from_atta = len(human_failures_from_atta[0]) #the number of tasks the human failed on
        for i in range(human_num_failures_from_atta):
            human_cost_for_failures_atta += (human_failures_from_atta[0][i] + human_failures_from_atta[1][i])/(3.0) #compute the human cost for every task the human failed on. This is the same human cost equation used during allocation.

        robot_num_failures_from_atta = len(robot_failures_from_atta[0]) #the number of tasks the robot failed on
        for i in range(robot_num_failures_from_atta):
            robot_cost_for_failures_atta += (robot_failures_from_atta[0][i] + robot_failures_from_atta[1][i])/(8.0) #compute the robot cost for every task the robot failed on. This is the same robot cost equation used during allocation.

        load_file_name = "./results/tsarouchi/tsarouchi_caseI_eta50_" + str(iter) + ".mat"
        output_tsarouchi_mat = sio.loadmat(load_file_name)

        human_successes_from_tsarouchi = output_tsarouchi_mat["human_successes"] #tasks assigned to the human that were successes
        human_failures_from_tsarouchi = output_tsarouchi_mat["human_failures"] #tasks assigned to the human that were failures
        robot_successes_from_tsarouchi = output_tsarouchi_mat["robot_successes"] #tasks assigned to the robot that were successes
        robot_failures_from_tsarouchi = output_tsarouchi_mat["robot_failures"] #tasks assigned to the robot that were failures
        human_num_tasks_from_tsarouchi = output_tsarouchi_mat["human_num_tasks"][0][0]
        robot_num_tasks_from_tsarouchi = output_tsarouchi_mat["robot_num_tasks"][0][0]
        discarded_num_tasks_from_tsarouchi = output_tsarouchi_mat["discarded_num_tasks"][0][0]
        discarded_num_robotCost_from_tsarouchi = output_tsarouchi_mat["discarded_num_robotCost"][0][0]
        discarded_num_humanCost_from_tsarouchi = output_tsarouchi_mat["discarded_num_humanCost"][0][0]

        num_successes_from_tsarouchi = human_successes_from_tsarouchi.shape[1] + robot_successes_from_tsarouchi.shape[1]
        num_failures_from_tsarouchi = human_failures_from_tsarouchi.shape[1] + robot_failures_from_tsarouchi.shape[1] + discarded_num_tasks_from_tsarouchi
        #human_num_tasks_from_tsarouchi
        #robot_num_tasks_from_tsarouchi
        #discarded_num_tasks_from_tsarouchi
        #discarded_num_humanCost_from_tsarouchi
        #discarded_num_robotCost_from_tsarouchi
        #human_successes_from_tsarouchi.shape[1]
        #human_failures_from_tsarouchi.shape[1]
        #robot_successes_from_tsarouchi.shape[1]
        #robot_failures_from_tsarouchi.shape[1]

        human_num_failures_from_tsarouchi = len(human_failures_from_tsarouchi[0]) #the number of tasks the human failed on
        for i in range(human_num_failures_from_tsarouchi):
            human_cost_for_failures_tsarouchi += (human_failures_from_tsarouchi[0][i] + human_failures_from_tsarouchi[1][i])/(3.0) #compute the human cost for every task the human failed on. This is the same human cost equation used during allocation.

        for i in range(discarded_num_humanCost_from_tsarouchi):
            human_cost_for_failures_tsarouchi += (0.65)/(3.0) #add the human cost for every task discarded to the human.
            #Note: the task requirements for tasks discarded to the human and counted as failures were not recorded during the simulation. However, one \lambda_i for every task must be > 0.65.

        robot_num_failures_from_tsarouchi = len(robot_failures_from_tsarouchi[0]) #the number of tasks the robot failed on
        for i in range(robot_num_failures_from_tsarouchi):
            robot_cost_for_failures_tsarouchi += (robot_failures_from_tsarouchi[0][i] + robot_failures_from_tsarouchi[1][i])/(8.0) #compute the robot cost for every task the robot failed on. This is the same robot cost equation used during allocation.

        for i in range(discarded_num_robotCost_from_tsarouchi):
            robot_cost_for_failures_tsarouchi += (0.65)/(8.0) #add the robot cost for every task discarded to the robot.
            #Note: the task requirements for tasks discarded to the robot and counted as failures were not recorded during the simulation. However, one \lambda_i for every task must be > 0.65.

        load_file_name = "./results/random/random_caseI_eta50_" + str(iter) + ".mat"
        output_random_mat = sio.loadmat(load_file_name)

        human_successes_from_random = output_random_mat["human_successes"] #tasks assigned to the human that were successes
        human_failures_from_random = output_random_mat["human_failures"] #tasks assigned to the human that were failures
        robot_successes_from_random = output_random_mat["robot_successes"] #tasks assigned to the robot that were successes
        robot_failures_from_random = output_random_mat["robot_failures"] #tasks assigned to the robot that were failures
        human_num_tasks_from_random = output_random_mat["human_num_tasks"][0][0]
        robot_num_tasks_from_random = output_random_mat["robot_num_tasks"][0][0]

        num_successes_from_random = human_successes_from_random.shape[1] + robot_successes_from_random.shape[1]
        num_failures_from_random = human_failures_from_random.shape[1] + robot_failures_from_random.shape[1]
        #human_num_tasks_from_random
        #robot_num_tasks_from_random
        #human_successes_from_random.shape[1]
        #human_failures_from_random.shape[1]
        #robot_successes_from_random.shape[1]
        #obot_failures_from_random.shape[1]

        human_num_failures_from_random = len(human_failures_from_random[0]) #the number of tasks the human failed on
        for i in range(human_num_failures_from_random):
            human_cost_for_failures_random += (human_failures_from_random[0][i] + human_failures_from_random[1][i])/(3.0) #compute the human cost for every task the human failed on. This is the same human cost equation used during allocation.

        robot_num_failures_from_random = len(robot_failures_from_random[0]) #the number of tasks the robot failed on
        for i in range(robot_num_failures_from_random):
            robot_cost_for_failures_random += (robot_failures_from_random[0][i] + robot_failures_from_random[1][i])/(8.0) #compute the robot cost for every task the robot failed on. This is the same robot cost equation used during allocation.

        row = [num_successes_from_atta, num_failures_from_atta, num_successes_from_random, num_failures_from_random,
        num_successes_from_tsarouchi, num_failures_from_tsarouchi,
        human_num_tasks_from_atta, robot_num_tasks_from_atta, human_num_tasks_from_random, robot_num_tasks_from_random,
        human_num_tasks_from_tsarouchi, robot_num_tasks_from_tsarouchi, discarded_num_tasks_from_tsarouchi, discarded_num_humanCost_from_tsarouchi, discarded_num_robotCost_from_tsarouchi,
        human_successes_from_atta.shape[1], human_failures_from_atta.shape[1], robot_successes_from_atta.shape[1], robot_failures_from_atta.shape[1],
        human_successes_from_random.shape[1], human_failures_from_random.shape[1], robot_successes_from_random.shape[1], robot_failures_from_random.shape[1],
        human_successes_from_tsarouchi.shape[1], human_failures_from_tsarouchi.shape[1], human_failures_from_tsarouchi.shape[1] + discarded_num_humanCost_from_tsarouchi,
        robot_successes_from_tsarouchi.shape[1], robot_failures_from_tsarouchi.shape[1], robot_failures_from_tsarouchi.shape[1] + discarded_num_robotCost_from_tsarouchi,
        human_cost_for_failures_atta, robot_cost_for_failures_atta, human_cost_for_failures_random, robot_cost_for_failures_random, human_cost_for_failures_tsarouchi, robot_cost_for_failures_tsarouchi]

        dataprep.append(row)

    df = pd.DataFrame(dataprep, columns=['ATTA Num Succ', 'ATTA Num Fail', 'Random Num Succ', 'Random Num Fail', 'Tsarouchi Num Succ', 'Tsarouchi Num Fail',
    'ATTA H Num Tasks', 'ATTA R Num Tasks', 'Random H Num Tasks', 'Random R Num Tasks',
    'Tsarouchi H Num Tasks', 'Tsarouchi R Num Tasks', 'Tsarouchi Discarded Num Tasks', 'Tsarouchi Discarded Num to H', 'Tsarouchi Discarded Num to R',
    'ATTA H Num Succ', 'ATTA H Num Fail', 'ATTA R Num Succ', 'ATTA R Num Fail', 'Random H Num Succ', 'Random H Num Fail', 'Random R Num Succ', 'Random R Num Fail',
    'Tsarouchi H Num Succ', 'Tsarouchi H Num Fail', 'Tsarouchi H Num Fail + Discarded', 'Tsarouchi R Num Succ', 'Tsarouchi R Num Fail', 'Tsarouchi R Num Fail + Discarded',
    'ATTA H Cost for Failures', 'ATTA R Cost for Failures', 'Random H Cost for Failures', 'Random R Cost for Failures', 'Tsarouchi MIN H Cost for Failures', 'Tsarouchi MIN R Cost for Failures'])
    df.to_csv('./results/spss/other_data_caseI.csv', index=False)

reason: import numpy,import scipy


text:
from __future__ import annotations

import inspect
import warnings
from typing import Any, Callable, Iterable, Mapping, Optional, Set, Union

import numpy as np
from numpy.typing import ArrayLike, NDArray

from .generation_surface import GenerationSurface
from .utils import get_column, get_table

class Weighter:
    """
    Abstract base class from which all weighers derive.

    Weighters will take a file object as input and calculate the weights of the events in the file
    for a given flux. As well as helper functions for all columns in the file. Weighters will keep
    track of the generation surface for the Monte Carlo in question. Weighters can be
    added together to form samples with different simulation parameters
    """

    def __init__(self, data: Iterable[tuple[Any, Mapping]], surface: GenerationSurface):
        colnames: Set[str] = set()
        for _, event_map in data:
            keys = set(event_map.keys())
            if colnames:
                colnames = colnames.intersection(keys)
            else:
                colnames = keys

        self.data = list(data)
        self.surface = surface
        self.colnames = sorted(colnames)
        self.__cache: dict = {}

    def get_column(self, table: str, column: str) -> NDArray[np.floating]:
        """
        Helper function to get a specific column from the file
        """
        retval: NDArray = np.array([])
        for datafile, _ in self.data:
            retval = np.append(retval, get_column(get_table(datafile, table), column))
        print("RETVAL", retval)
        return retval

    def get_weight_column(self, name: str) -> NDArray[np.floating]:
        """
        Helper function to get a column needed in the weight calculation
        """
        if name in self.__cache:
            return self.__cache[name]

        if name == "cos_zen":
            retval = np.cos(self.get_weight_column("zenith"))

        else:
            retval = np.array([])
            for datafile, event_map in self.data:
                if event_map[name] is None:
                    tablename, columnname = event_map["energy"]
                    fileval = np.full(len(get_column(get_table(datafile, tablename), columnname)), 1)
                else:
                    tablename, columnname = event_map[name]
                    fileval = get_column(get_table(datafile, tablename), columnname)
                retval = np.append(retval, fileval)
            if name == "pdgid":
                retval = retval.astype(np.int32)

        self.__cache[name] = retval
        return retval

    def get_weights(self, flux: Any) -> NDArray[np.floating]:
        """
        Calculate the weights for the sample for the given flux.

        Args:
          flux: can be any one of several types:

            * An instance of :py:class:`nuflux.FluxFunction` from
              `nuflux <https://github.com/icecube/nuflux>`_
            * A callable where the names of the arguments are match the weight objects weighting columns.
            * An iterable the same length os the datasample. If you have other means of calculating
              the flux for each event than the above options this can be useful.
            * A scalar number. This calculates the unrealistic scenario where all events have the same
              flux, this can be useful for testing or calculating effective areas. If the value is 1 then
              the return value will be the well known quantity OneWeight.
        """

        event_col = {k: self.get_weight_column(k) for k in ["energy", "pdgid", "cos_zen"]}

        # do nothing if everything is empty
        if event_col["pdgid"].shape == (0,):
            return np.array([])

        epdf = self.surface.get_epdf(**event_col)
        event_weight = self.get_weight_column("event_weight")

        # calculate the flux based on which type of flux it is
        if hasattr(flux, "getFlux"):
            # this is a nuflux model
            assert callable(flux.getFlux)
            flux_val = flux.getFlux(
                self.get_weight_column("pdgid"),
                self.get_weight_column("energy"),
                self.get_weight_column("cos_zen"),
            )
        elif callable(flux):
            # this is a cosmic ray flux model or just a function
            keys = inspect.signature(flux).parameters.keys()
            try:
                arguments = {k: self.get_weight_column(k) for k in keys}
            except KeyError as missing_params:
                raise ValueError(
                    f"get_weights() was passed callable {repr(flux)} which has parameters {list(keys)}. "
                    "The weight columns which are available are {repr(self.colnames)}"
                ) from missing_params
            flux_val = flux(**arguments)
        elif hasattr(flux, "__len__"):
            # this is an array with a length equal to the number of events
            flux_val = np.array(flux, dtype=float)
        elif np.isscalar(flux):
            # this is a scalar
            flux_val = np.full(epdf.shape, flux)
        else:
            raise ValueError(f"I do not understand what to do with flux {flux}")
        assert flux_val.shape == epdf.shape

        # Getting events with epdf=0 indicates some sort of mismatch between the
        # the surface and the dataset that can't be solved here so print a
        # warning and ignore the events
        mask = epdf > 0

        if not np.all(mask):
            warnings.warn(
                f"simweights :: {np.logical_not(mask).sum()} events out of {mask.size} were found to be "
                "outside the generation surface. This could indicate a problem with this dataset."
            )
        weights = np.zeros_like(epdf)
        weights[mask] = (event_weight * flux_val)[mask] / epdf[mask]
        return weights

    def effective_area(
        self, energy_bins: ArrayLike, cos_zenith_bins: ArrayLike, mask: Optional[ArrayLike] = None
    ) -> NDArray[np.floating]:
        r"""
        Calculate The effective area for the given energy and zenith bins.

        This is accomplished by histogramming the generation surface the simulation sample
        in energy and zenith bins and dividing by the size of the energy and solid angle of each bin.
        If mask is passed as a parameter, only events which are included in the mask are used.
        Effective areas are given units of :math:`\mathrm{m}^2`

        .. Note ::

            If the sample contains more than one type of primary particle, then the result will be
            averaged over the number of particles. This is usually what you want. However, it can
            cause strange behavior if there is a small number of one type. In this case, the mask
            should be used to select the particle types individually.

        Args:
            energy_bins(array_like): A length N+1 array of energy bin edges
            coz_zenith_bins(array_like): A length M+1 array of cos(zenith) bin edges
            mask(array_like): boolean array where 1 indicates to use the event in the calculation.
                Must have the same length as the simulation sample.

        Returns:
            array_like: An NxM array of effective areas. Where N is the number of energy bins and
                M is the number of cos(zenith) bins.

        """

        energy_bins = np.array(energy_bins)
        cos_zenith_bins = np.array(cos_zenith_bins)

        assert energy_bins.ndim == 1
        assert cos_zenith_bins.ndim == 1
        assert len(energy_bins) >= 2
        assert len(cos_zenith_bins) >= 2

        energy = self.get_weight_column("energy")
        cos_zen = np.cos(self.get_weight_column("zenith"))

        weights = self.get_weights(1e-4)
        if mask is None:
            maska = np.full(weights.size, 1, dtype=bool)
        else:
            maska = np.asarray(mask, dtype=bool)

        assert maska.shape == weights.shape

        hist_val, czbin, enbin = np.histogram2d(
            cos_zen[maska],
            energy[maska],
            weights=weights[maska],
            bins=[cos_zenith_bins, energy_bins],
        )

        nspecies = len(np.unique(self.get_weight_column("pdgid")[maska]))

        assert np.array_equal(enbin, energy_bins)
        assert np.array_equal(czbin, cos_zenith_bins)
        e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin))
        return hist_val / (e_width * 2 * np.pi * z_width * nspecies)

    def __add__(self, other: Weighter) -> Weighter:
        return Weighter(self.data + other.data, self.surface + other.surface)

    def tostring(self, flux: Union[None, object, Callable, ArrayLike] = None) -> str:
        """
        Creates a string with important information about this weighting object:
        generation surface, event map, number of events, and effective area.
        if optional flux is provided the event rate and livetime are added as well.

        """
        output = str(self.surface) + "\n"
        output += f"Weight Columns   : {', '.join(self.colnames)}\n"
        output += f"Number of Events : {len(self.get_weights(1)):8d}\n"
        eff_area = self.effective_area(
            self.surface.get_energy_range(None), self.surface.get_cos_zenith_range(None)
        )
        output += f"Effective Area   : {eff_area[0][0]:8.6g} m²\n"
        if flux:
            weights = self.get_weights(flux)
            output += f"Using flux model : {flux.__class__.__name__}\n"
            output += f"Event Rate       : {weights.sum():8.6g} Hz\n"
            output += f"Livetime         : {weights.sum() / (weights ** 2).sum():8.6g} s\n"
        return output

    def __str__(self) -> str:
        return self.tostring()

max_stars_repo_name: icecube/simweights

max_issues_repo_name: icecube/simweights

max_forks_repo_name: icecube/simweights

reason: import numpy,from numpy


text:
#!/usr/bin/env python

from __future__ import division

import argparse
import io
import math
import os
import random

import numpy as np
import tensorflow as tf

SCRIPT_PATH = os.path.dirname(os.path.abspath(__file__))

# Default data paths.
DEFAULT_LABEL_CSV = os.path.join(SCRIPT_PATH, '../image-data/labels-map.csv')
DEFAULT_LABEL_FILE = os.path.join(SCRIPT_PATH,
                                  '../labels/2350-common-hangul.txt')
DEFAULT_OUTPUT_DIR = os.path.join(SCRIPT_PATH, '../tfrecords-output')
DEFAULT_NUM_SHARDS_TRAIN = 3
DEFAULT_NUM_SHARDS_TEST = 1

def _int64_feature(value):
    return tf.train.Feature(int64_list=tf.train.Int64List(value=[value]))

def _bytes_feature(value):
    return tf.train.Feature(bytes_list=tf.train.BytesList(value=[value]))

class TFRecordsConverter(object):
    """Class that handles converting images to TFRecords."""

    def __init__(self, labels_csv, label_file, output_dir,
                 num_shards_train, num_shards_test):

        self.output_dir = output_dir
        self.num_shards_train = num_shards_train
        self.num_shards_test = num_shards_test

        if not os.path.exists(self.output_dir):
            os.makedirs(self.output_dir)

        # Get lists of images and labels.
        self.filenames, self.labels = \
            self.process_image_labels(labels_csv, label_file)

        # Counter for total number of images processed.
        self.counter = 0

    def process_image_labels(self, labels_csv, label_file):
        """This will constuct two shuffled lists for images and labels.

        The index of each image in the images list will have the corresponding
        label at the same index in the labels list.
        """
        labels_csv = io.open(labels_csv, 'r', encoding='utf-8')
        labels_file = io.open(label_file, 'r',
                              encoding='utf-8').read().splitlines()

        # Map characters to indices.
        label_dict = {}
        count = 0
        for label in labels_file:
            label_dict[label] = count
            count += 1

        # Build the lists.
        images = []
        labels = []
        for row in labels_csv:
            file, label = row.strip().split(',')
            images.append(file)
            labels.append(label_dict[label])

        # Randomize the order of all the images/labels.
        shuffled_indices = list(range(len(images)))
        random.seed(12121)
        random.shuffle(shuffled_indices)
        filenames = [images[i] for i in shuffled_indices]
        labels = [labels[i] for i in shuffled_indices]

        return filenames, labels

    def write_tfrecords_file(self, output_path, indices):
        """Writes out TFRecords file."""
        writer = tf.python_io.TFRecordWriter(output_path)
        for i in indices:
            filename = self.filenames[i]
            label = self.labels[i]
            with tf.gfile.GFile(filename, 'rb') as f:
                im_data = f.read()

            # Example is a data format that contains a key-value store, where
            # each key maps to a Feature message. In this case, each Example
            # contains two features. One will be a ByteList for the raw image
            # data and the other will be an Int64List containing the index of
            # the corresponding label in the labels list from the file.
            example = tf.train.Example(features=tf.train.Features(feature={
                'image/class/label': _int64_feature(label),
                'image/encoded': _bytes_feature(tf.compat.as_bytes(im_data))}))
            writer.write(example.SerializeToString())
            self.counter += 1
            if not self.counter % 1000:
                print('Processed {} images...'.format(self.counter))
        writer.close()

    def convert(self):
        """This function will drive the conversion to TFRecords.

        Here, we partition the data into a training and testing set, then
        divide each data set into the specified number of TFRecords shards.
        """

        num_files_total = len(self.filenames)

        # Allocate about 15 percent of images to testing
        num_files_test = int(num_files_total * .15)

        # About 85 percent will be for training.
        num_files_train = num_files_total - num_files_test

        print('Processing training set TFRecords...')

        files_per_shard = int(math.ceil(num_files_train /
                                        self.num_shards_train))
        start = 0
        for i in range(1, self.num_shards_train):
            shard_path = os.path.join(self.output_dir,
                                      'train-{}.tfrecords'.format(str(i)))
            # Get a subset of indices to get only a subset of images/labels for
            # the current shard file.
            file_indices = np.arange(start, start+files_per_shard, dtype=int)
            start = start + files_per_shard
            self.write_tfrecords_file(shard_path, file_indices)

        # The remaining images will go in the final shard.
        file_indices = np.arange(start, num_files_train, dtype=int)
        final_shard_path = os.path.join(self.output_dir,
                                        'train-{}.tfrecords'.format(
                                            str(self.num_shards_train)))
        self.write_tfrecords_file(final_shard_path, file_indices)

        print('Processing testing set TFRecords...')

        files_per_shard = math.ceil(num_files_test / self.num_shards_test)
        start = num_files_train
        for i in range(1, self.num_shards_test):
            shard_path = os.path.join(self.output_dir,
                                      'test-{}.tfrecords'.format(str(i)))
            file_indices = np.arange(start, start+files_per_shard, dtype=int)
            start = start + files_per_shard
            self.write_tfrecords_file(shard_path, file_indices)

        # The remaining images will go in the final shard.
        file_indices = np.arange(start, num_files_total, dtype=int)
        final_shard_path = os.path.join(self.output_dir,
                                        'test-{}.tfrecords'.format(
                                            str(self.num_shards_test)))
        self.write_tfrecords_file(final_shard_path, file_indices)

        print('\nProcessed {} total images...'.format(self.counter))
        print('Number of training examples: {}'.format(num_files_train))
        print('Number of testing examples: {}'.format(num_files_test))
        print('TFRecords files saved to {}'.format(self.output_dir))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--image-label-csv', type=str, dest='labels_csv',
                        default=DEFAULT_LABEL_CSV,
                        help='File containing image paths and corresponding '
                             'labels.')
    parser.add_argument('--label-file', type=str, dest='label_file',
                        default=DEFAULT_LABEL_FILE,
                        help='File containing newline delimited labels.')
    parser.add_argument('--output-dir', type=str, dest='output_dir',
                        default=DEFAULT_OUTPUT_DIR,
                        help='Output directory to store TFRecords files.')
    parser.add_argument('--num-shards-train', type=int,
                        dest='num_shards_train',
                        default=DEFAULT_NUM_SHARDS_TRAIN,
                        help='Number of shards to divide training set '
                             'TFRecords into.')
    parser.add_argument('--num-shards-test', type=int,
                        dest='num_shards_test',
                        default=DEFAULT_NUM_SHARDS_TEST,
                        help='Number of shards to divide testing set '
                             'TFRecords into.')
    args = parser.parse_args()
    converter = TFRecordsConverter(args.labels_csv,
                                   args.label_file,
                                   args.output_dir,
                                   args.num_shards_train,
                                   args.num_shards_test)
    converter.convert()


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 24 15:50:24 2020

@author: petrapoklukar

Computes recall metric.
"""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from disentanglement_lib.evaluation.metrics import utils
from disentanglement_lib.evaluation.metrics import iprd_score as iprd
import numpy as np
from six.moves import range
import gin.tf
from disentanglement_lib.data.ground_truth import named_data
import tensorflow as tf
from sklearn.decomposition import PCA
import pickle

@gin.configurable(
    "recall_on_holdout",
    blacklist=["holdout_ground_truth_data", "encoder_fn", "repr_transform_fn",
               "decoder_fn", "random_state", "artifact_dir"])
def compute_recall_on_holdout(holdout_ground_truth_data,
                              encoder_fn,
                              repr_transform_fn,
                              decoder_fn,
                              random_state,
                              artifact_dir=None,
                              nhood_sizes=gin.REQUIRED,
                              num_interventions_per_latent_dim=gin.REQUIRED,
                              pca_components=gin.REQUIRED
                              ):
  """TBA

  Args:
    ground_truth_data: GroundTruthData to be sampled from.

    random_state: Numpy random state used for randomness.
    artifact_dir: Optional path to directory where artifacts can be saved.
  """
  del artifact_dir
  num_recall_samples = holdout_ground_truth_data.data_size
  dummy_input = holdout_ground_truth_data.sample_observations(1, random_state)
  dummy_mean, dummy_var = encoder_fn(dummy_input)

  # Samples from the normal prior
  latent_dim = repr_transform_fn(*encoder_fn(dummy_input)).shape[-1]
  latent_shape = [num_recall_samples, latent_dim]
  latent_prior_samples_np = np.random.normal(size=latent_shape)

  # Ground truth samples
  gt_samples = holdout_ground_truth_data.images
  gt_repr = repr_transform_fn(*encoder_fn(gt_samples))
  gt_repr_mean = np.mean(gt_repr, axis=0) # latent_shape
  gt_repr_std = np.std(gt_repr, axis=0)
  gt_repr_min = np.min(gt_repr, axis=0)
  gt_repr_max = np.max(gt_repr, axis=0)

  # The predetermined set of interventions from the estimated training prior
  fixed_trained_prior_samples_np = np.random.normal(loc=gt_repr_mean,
                                                    scale=gt_repr_std,
                                                    size=latent_shape)
  print(fixed_trained_prior_samples_np.shape, fixed_trained_prior_samples_np)

  pca_results_per_comp = []
  result_d_gen = {'nhoods': nhood_sizes,
              'gt_repr_mean': list(gt_repr_mean),
              'gt_repr_std': list(gt_repr_std),
              'gt_repr_min': list(gt_repr_min),
              'gt_repr_max': list(gt_repr_max)
              }

  subset_interventions = np.random.choice(
      np.arange(num_recall_samples), size=num_interventions_per_latent_dim,
      replace=False)
  result_d_gen['subset_interventions'] = list(subset_interventions)

  for num_comp in pca_components:
    # Load the pretrained gt pca
    pca_path = "backbone/pca/pca_3dshapes_model_all_{0}comp.pkl".format(num_comp)
    print(pca_path)
    with open(pca_path, 'rb') as f:
      gt_pca = pickle.load(f)

    result_d = result_d_gen.copy()
    result_d['pca_components'] = num_comp
    sess = tf.Session()
    with sess.as_default():
      # Choose a subset of interventions

      print('\n\n\n Computing the total recall...')
      # Sample ground truth data and vae process it
      decoded_gt_samples = decoder_fn(gt_repr)
      decoded_gt_samples = decoded_gt_samples.reshape(num_recall_samples, -1)
      reduced_decoded_gt_samples = gt_pca.transform(decoded_gt_samples)

      # Generated samples from the normal prior, processed with gt PCA
      generated_prior_samples = decoder_fn(latent_prior_samples_np)
      generated_prior_samples = generated_prior_samples.reshape(num_recall_samples, -1)
      reduced_generated_prior_samples = gt_pca.transform(generated_prior_samples)
      assert(reduced_generated_prior_samples.shape == reduced_decoded_gt_samples.shape)

      # Generated samples from the estimated training prior, processed with gt PCA
      generated_trained_prior_samples = decoder_fn(fixed_trained_prior_samples_np)
      generated_trained_prior_samples = generated_trained_prior_samples.reshape(num_recall_samples, -1)
      reduced_generated_trained_prior_samples = gt_pca.transform(generated_trained_prior_samples)
      assert(reduced_generated_prior_samples.shape == reduced_generated_trained_prior_samples.shape)

      # --- Original sharp images - discarded to now
  #    gt_train_samples = gt_samples.reshape(num_recall_samples, -1)
  #    gt_pca = PCA(n_components=n_comp)
  #    reduced_gt_samples = gt_pca.fit_transform(gt_samples)
  #    assert(reduced_gt_samples.shape == reduced_decoded_gt_samples.shape)

  #    # compute model recall: gt vs generated
  #    gt_generated_result = iprd.knn_precision_recall_features(
  #        reduced_gt_samples,
  #        reduced_generated_prior_samples,
  #        nhood_sizes=nhood_sizes,
  #        row_batch_size=500, col_batch_size=100, num_gpus=1)
  #    update_result_dict(result_d, ['gt_generated_', gt_generated_result])
      # ----

      # compute model recall: model(gt) vs normal prior generated
      decoded_gt_prior_generated_result = iprd.knn_precision_recall_features(
          reduced_decoded_gt_samples,
          reduced_generated_prior_samples,
          nhood_sizes=nhood_sizes,
          row_batch_size=500, col_batch_size=100, num_gpus=1)
      update_result_dict(result_d, ['decoded_gt_prior_generated_',
                                    decoded_gt_prior_generated_result])

      # compute model recall: model(gt) vs estimated training prior generated
      decoded_gt_trained_prior_generated_result = iprd.knn_precision_recall_features(
          reduced_decoded_gt_samples,
          reduced_generated_trained_prior_samples,
          nhood_sizes=nhood_sizes,
          row_batch_size=500, col_batch_size=100, num_gpus=1)
      update_result_dict(result_d, ['decoded_gt_trained_prior_generated_',
                                    decoded_gt_trained_prior_generated_result])

      # compute model recall:normal prior generated vs estimated training prior generated
      prior_generated_trained_prior_generated_result = iprd.knn_precision_recall_features(
          reduced_generated_prior_samples,
          reduced_generated_trained_prior_samples,
          nhood_sizes=nhood_sizes,
          row_batch_size=500, col_batch_size=100, num_gpus=1)
      update_result_dict(result_d, ['prior_generated_trained_prior_generated_',
                                    prior_generated_trained_prior_generated_result])

      # Pick a latent dimension
      for dim in range(latent_dim):
        print('\n\n\n Computing the recall for latent dim ', dim)

        agg_fix_one_vs_prior_generated_result = {'precision': [], 'recall': []}
        agg_fix_one_vs_trained_prior_generated_result = {'precision': [], 'recall': []}
        agg_fix_one_vs_decoded_gt_result = {'precision': [], 'recall': []}

        agg_vary_one_vs_prior_generated_result = {'precision': [], 'recall': []}
        agg_vary_one_vs_trained_prior_generated_result = {'precision': [], 'recall': []}
        agg_vary_one_vs_decoded_gt_result =  {'precision': [], 'recall': []}

        # intervene several times
        for intervention in range(num_interventions_per_latent_dim):
          inter_id = subset_interventions[intervention]
          print('n\n\n Intervention num', intervention)
          print(' Intervention row', inter_id)

          # --- fix one, vary the rest
          latent_intervention = float(fixed_trained_prior_samples_np[inter_id, dim])
          fix_one_latent_from_trained_prior_samples = np.copy(fixed_trained_prior_samples_np)
          fix_one_latent_from_trained_prior_samples[:, dim] = latent_intervention
          # decode the samples and trasform them with the PCA
          gen_fix_one_latent_from_trained_prior_samples = decoder_fn(
              fix_one_latent_from_trained_prior_samples).reshape(num_recall_samples, -1)
          reduced_gen_fix_one_latent_from_trained_prior_samples = gt_pca.transform(
              gen_fix_one_latent_from_trained_prior_samples)
          assert(reduced_gen_fix_one_latent_from_trained_prior_samples.shape == \
                 reduced_generated_prior_samples.shape)

          # - Calculate the relative recall
          # compare to the generated normal prior samples
          fix_one_vs_prior_generated_result = iprd.knn_precision_recall_features(
              reduced_generated_prior_samples,
              reduced_gen_fix_one_latent_from_trained_prior_samples,
              nhood_sizes=nhood_sizes,
              row_batch_size=500, col_batch_size=100, num_gpus=1)
          print(agg_fix_one_vs_prior_generated_result)
          print(fix_one_vs_prior_generated_result)
          agg_fix_one_vs_prior_generated_result = agg_recall_dict(
              agg_fix_one_vs_prior_generated_result,
              fix_one_vs_prior_generated_result,
              intervention)
          print(agg_fix_one_vs_prior_generated_result)

          # compare to the generated trained prior samples
          fix_one_vs_trained_prior_generated_result = iprd.knn_precision_recall_features(
              reduced_generated_trained_prior_samples,
              reduced_gen_fix_one_latent_from_trained_prior_samples,
              nhood_sizes=nhood_sizes,
              row_batch_size=500, col_batch_size=100, num_gpus=1)
          agg_fix_one_vs_trained_prior_generated_result = agg_recall_dict(
              agg_fix_one_vs_trained_prior_generated_result,
              fix_one_vs_trained_prior_generated_result,
              intervention)

          # compare to the decoded gt
          fix_one_vs_decoded_gt_result = iprd.knn_precision_recall_features(
              reduced_decoded_gt_samples,
              reduced_gen_fix_one_latent_from_trained_prior_samples,
              nhood_sizes=nhood_sizes,
              row_batch_size=500, col_batch_size=100, num_gpus=1)
          agg_fix_one_vs_decoded_gt_result = agg_recall_dict(
              agg_fix_one_vs_decoded_gt_result,
              fix_one_vs_decoded_gt_result,
              intervention)

          # --- vary one, fix the rest
          latent_variation = np.copy(fixed_trained_prior_samples_np[:, dim])
          print(latent_variation.shape, latent_variation)
          vary_one_latent_from_trained_prior_samples = np.copy(
              fixed_trained_prior_samples_np[inter_id]).reshape(1, latent_dim)
          vary_one_latent_from_trained_prior_samples = np.full(
              latent_shape, vary_one_latent_from_trained_prior_samples)
          print(vary_one_latent_from_trained_prior_samples.shape)
          print(vary_one_latent_from_trained_prior_samples)

          vary_one_latent_from_trained_prior_samples[:, dim] = latent_variation
          print(vary_one_latent_from_trained_prior_samples.shape)
          print(vary_one_latent_from_trained_prior_samples)

          # decode the samples and trasform them with the PCA
          gen_vary_one_latent_from_trained_prior_samples = decoder_fn(
              vary_one_latent_from_trained_prior_samples).reshape(num_recall_samples, -1)
          reduced_gen_vary_one_latent_from_trained_prior_samples = gt_pca.transform(
              gen_vary_one_latent_from_trained_prior_samples)
          assert(reduced_gen_vary_one_latent_from_trained_prior_samples.shape == \
                 reduced_generated_prior_samples.shape)

          # - Calculate the recall
          # compare to the generated normal prior samples
          vary_one_vs_prior_generated_result = iprd.knn_precision_recall_features(
              reduced_generated_prior_samples,
              reduced_gen_vary_one_latent_from_trained_prior_samples,
              nhood_sizes=nhood_sizes,
              row_batch_size=500, col_batch_size=100, num_gpus=1)
          print(agg_vary_one_vs_prior_generated_result)
          print(vary_one_vs_prior_generated_result)
          agg_vary_one_vs_prior_generated_result = agg_recall_dict(
              agg_vary_one_vs_prior_generated_result,
              vary_one_vs_prior_generated_result,
              intervention)
          print(agg_vary_one_vs_prior_generated_result)

          # compare to the generated trained prior samples
          vary_one_vs_trained_prior_generated_result = iprd.knn_precision_recall_features(
              reduced_generated_trained_prior_samples,
              reduced_gen_vary_one_latent_from_trained_prior_samples,
              nhood_sizes=nhood_sizes,
              row_batch_size=500, col_batch_size=100, num_gpus=1)
          agg_vary_one_vs_trained_prior_generated_result = agg_recall_dict(
              agg_vary_one_vs_trained_prior_generated_result,
              vary_one_vs_trained_prior_generated_result,
              intervention)

          # compare to the decoded gt
          vary_one_vs_decoded_gt_result = iprd.knn_precision_recall_features(
              reduced_decoded_gt_samples,
              reduced_gen_vary_one_latent_from_trained_prior_samples,
              nhood_sizes=nhood_sizes,
              row_batch_size=500, col_batch_size=100, num_gpus=1)
          agg_vary_one_vs_decoded_gt_result = agg_recall_dict(
              agg_vary_one_vs_decoded_gt_result,
              vary_one_vs_decoded_gt_result,
              intervention)

        update_result_dict_with_agg(
            result_d,
            [str(dim) + '_fix_one_vs_prior_generated_', agg_fix_one_vs_prior_generated_result],
            [str(dim) + '_fix_one_vs_trained_prior_generated_', agg_fix_one_vs_trained_prior_generated_result],
            [str(dim) + '_fix_one_vs_decoded_gt_', agg_fix_one_vs_decoded_gt_result],
            [str(dim) + '_vary_one_vs_prior_generated_', agg_vary_one_vs_prior_generated_result],
            [str(dim) + '_vary_one_vs_trained_prior_generated_', agg_vary_one_vs_trained_prior_generated_result],
            [str(dim) + '_vary_one_vs_decoded_gt_', agg_vary_one_vs_decoded_gt_result])
    pca_results_per_comp.append(result_d)
  return pca_results_per_comp

@gin.configurable(
    "recall",
    blacklist=["ground_truth_data", "encoder_fn", "repr_transform_fn",
               "decoder_fn", "random_state", "artifact_dir"])
def compute_recall(ground_truth_data,
                   encoder_fn,
                   repr_transform_fn,
                   decoder_fn,
                   random_state,
                   artifact_dir=None,
#                   num_recall_samples=gin.REQUIRED,
                   nhood_sizes=gin.REQUIRED,
                   num_interventions_per_latent_dim=gin.REQUIRED,
                   num_pca_components=gin.REQUIRED
                   ):
  """TBA

  Args:

    random_state: Numpy random state used for randomness.
    artifact_dir: Optional path to directory where artifacts can be saved.
  """
  del artifact_dir
  train_ground_truth_data, test_ground_truth_data = ground_truth_data
  ground_truth_data = train_ground_truth_data
  num_recall_samples = train_ground_truth_data.data_size
  dummy_input = ground_truth_data.sample_observations(1, random_state)
  dummy_mean, dummy_var = encoder_fn(dummy_input)

  # Samples from the normal prior
  latent_dim = repr_transform_fn(*encoder_fn(dummy_input)).shape[-1]
  latent_shape = [num_recall_samples, latent_dim]
  latent_prior_samples_np = np.random.normal(size=latent_shape)

  # Ground truth samples
  gt_train_samples = ground_truth_data.sample_observations(num_recall_samples, random_state)
  gt_train_repr = repr_transform_fn(*encoder_fn(gt_train_samples))
  gt_train_repr_mean = np.mean(gt_train_repr, axis=0) # latent_shape
  gt_train_repr_std = np.std(gt_train_repr, axis=0)
  gt_train_repr_min = np.min(gt_train_repr, axis=0)
  gt_train_repr_max = np.max(gt_train_repr, axis=0)

  # The predetermined set of interventions from the estimated training prior
  fixed_trained_prior_samples_np = np.random.normal(loc=gt_train_repr_mean,
                                              scale=gt_train_repr_std,
                                              size=latent_shape)
  print(fixed_trained_prior_samples_np.shape, fixed_trained_prior_samples_np)

  result_d = {'nhoods': nhood_sizes,
              'gt_train_repr_mean': list(gt_train_repr_mean),
              'gt_train_repr_std': list(gt_train_repr_std),
              'gt_train_repr_min': list(gt_train_repr_min),
              'gt_train_repr_max': list(gt_train_repr_max)}
  sess = tf.Session()
  with sess.as_default():
    n_comp = min(num_recall_samples, num_pca_components)

    print('\n\n\n Computing the total recall...')
    # Sample ground truth data and vae process it
    decoded_gt_samples = decoder_fn(repr_transform_fn(*encoder_fn(gt_train_samples)))
    decoded_gt_samples = decoded_gt_samples.reshape(num_recall_samples, -1)
    decoded_gt_pca = PCA(n_components=n_comp)
    reduced_decoded_gt_samples = decoded_gt_pca.fit_transform(decoded_gt_samples)

    # Generated samples from the normal prior, processed with gt PCA
    generated_prior_samples = decoder_fn(latent_prior_samples_np)
    generated_prior_samples = generated_prior_samples.reshape(num_recall_samples, -1)
    reduced_generated_prior_samples = decoded_gt_pca.transform(generated_prior_samples)
    assert(reduced_generated_prior_samples.shape == reduced_decoded_gt_samples.shape)

    # Generated samples from the estimated training prior, processed with gt PCA
    generated_trained_prior_samples = decoder_fn(fixed_trained_prior_samples_np)
    generated_trained_prior_samples = generated_trained_prior_samples.reshape(num_recall_samples, -1)
    reduced_generated_trained_prior_samples = decoded_gt_pca.transform(generated_trained_prior_samples)
    assert(reduced_generated_prior_samples.shape == reduced_generated_trained_prior_samples.shape)

    # --- Original sharp images - discarded to now
#    gt_train_samples = gt_samples.reshape(num_recall_samples, -1)
#    gt_pca = PCA(n_components=n_comp)
#    reduced_gt_samples = gt_pca.fit_transform(gt_samples)
#    assert(reduced_gt_samples.shape == reduced_decoded_gt_samples.shape)

#    # compute model recall: gt vs generated
#    gt_generated_result = iprd.knn_precision_recall_features(
#        reduced_gt_samples,
#        reduced_generated_prior_samples,
#        nhood_sizes=nhood_sizes,
#        row_batch_size=500, col_batch_size=100, num_gpus=1)
#    update_result_dict(result_d, ['gt_generated_', gt_generated_result])
    # ----

    # compute model recall: model(gt) vs normal prior generated
    decoded_gt_prior_generated_result = iprd.knn_precision_recall_features(
        reduced_decoded_gt_samples,
        reduced_generated_prior_samples,
        nhood_sizes=nhood_sizes,
        row_batch_size=500, col_batch_size=100, num_gpus=1)
    update_result_dict(result_d, ['decoded_gt_prior_generated_',
                                  decoded_gt_prior_generated_result])

    # compute model recall: model(gt) vs estimated training prior generated
    decoded_gt_trained_prior_generated_result = iprd.knn_precision_recall_features(
        reduced_decoded_gt_samples,
        reduced_generated_trained_prior_samples,
        nhood_sizes=nhood_sizes,
        row_batch_size=500, col_batch_size=100, num_gpus=1)
    update_result_dict(result_d, ['decoded_gt_trained_prior_generated_',
                                  decoded_gt_trained_prior_generated_result])

    # compute model recall:normal prior generated vs estimated training prior generated
    prior_generated_trained_prior_generated_result = iprd.knn_precision_recall_features(
        reduced_generated_prior_samples,
        reduced_generated_trained_prior_samples,
        nhood_sizes=nhood_sizes,
        row_batch_size=500, col_batch_size=100, num_gpus=1)
    update_result_dict(result_d, ['prior_generated_trained_prior_generated_',
                                  prior_generated_trained_prior_generated_result])

    # Choose a subset of interventions
    subset_interventions = np.random.choice(
          np.arange(num_recall_samples), size=num_interventions_per_latent_dim,
          replace=False)
    result_d['subset_interventions'] = list(subset_interventions)
    result_d['num_pca_comp'] = n_comp

    # Pick a latent dimension
    for dim in range(latent_dim):
      print('\n\n\n Computing the recall for latent dim ', dim)

      agg_fix_one_vs_prior_generated_result = {'precision': [], 'recall': []}
      agg_fix_one_vs_trained_prior_generated_result = {'precision': [], 'recall': []}
      agg_fix_one_vs_decoded_gt_result = {'precision': [], 'recall': []}

      agg_vary_one_vs_prior_generated_result = {'precision': [], 'recall': []}
      agg_vary_one_vs_trained_prior_generated_result = {'precision': [], 'recall': []}
      agg_vary_one_vs_decoded_gt_result =  {'precision': [], 'recall': []}

      # intervene several times
      for intervention in range(num_interventions_per_latent_dim):
        inter_id = subset_interventions[intervention]
        print('n\n\n Intervention num', intervention)
        print(' Intervention row', inter_id)

        # --- fix one, vary the rest
        latent_intervention = float(fixed_trained_prior_samples_np[inter_id, dim])
        fix_one_latent_from_trained_prior_samples = np.copy(fixed_trained_prior_samples_np)
        fix_one_latent_from_trained_prior_samples[:, dim] = latent_intervention
        # decode the samples and trasform them with the PCA
        gen_fix_one_latent_from_trained_prior_samples = decoder_fn(
            fix_one_latent_from_trained_prior_samples).reshape(num_recall_samples, -1)
        reduced_gen_fix_one_latent_from_trained_prior_samples = decoded_gt_pca.transform(
            gen_fix_one_latent_from_trained_prior_samples)
        assert(reduced_gen_fix_one_latent_from_trained_prior_samples.shape == \
               reduced_generated_prior_samples.shape)

        # - Calculate the relative recall
        # compare to the generated normal prior samples
        fix_one_vs_prior_generated_result = iprd.knn_precision_recall_features(
            reduced_generated_prior_samples,
            reduced_gen_fix_one_latent_from_trained_prior_samples,
            nhood_sizes=nhood_sizes,
            row_batch_size=500, col_batch_size=100, num_gpus=1)
        print(agg_fix_one_vs_prior_generated_result)
        print(fix_one_vs_prior_generated_result)
        agg_fix_one_vs_prior_generated_result = agg_recall_dict(
            agg_fix_one_vs_prior_generated_result,
            fix_one_vs_prior_generated_result,
            intervention)
        print(agg_fix_one_vs_prior_generated_result)

        # compare to the generated trained prior samples
        fix_one_vs_trained_prior_generated_result = iprd.knn_precision_recall_features(
            reduced_generated_trained_prior_samples,
            reduced_gen_fix_one_latent_from_trained_prior_samples,
            nhood_sizes=nhood_sizes,
            row_batch_size=500, col_batch_size=100, num_gpus=1)
        agg_fix_one_vs_trained_prior_generated_result = agg_recall_dict(
            agg_fix_one_vs_trained_prior_generated_result,
            fix_one_vs_trained_prior_generated_result,
            intervention)

        # compare to the decoded gt
        fix_one_vs_decoded_gt_result = iprd.knn_precision_recall_features(
            reduced_decoded_gt_samples,
            reduced_gen_fix_one_latent_from_trained_prior_samples,
            nhood_sizes=nhood_sizes,
            row_batch_size=500, col_batch_size=100, num_gpus=1)
        agg_fix_one_vs_decoded_gt_result = agg_recall_dict(
            agg_fix_one_vs_decoded_gt_result,
            fix_one_vs_decoded_gt_result,
            intervention)

        # --- vary one, fix the rest
        latent_variation = np.copy(fixed_trained_prior_samples_np[:, dim])
        print(latent_variation.shape, latent_variation)
        vary_one_latent_from_trained_prior_samples = np.copy(
            fixed_trained_prior_samples_np[inter_id]).reshape(1, latent_dim)
        vary_one_latent_from_trained_prior_samples = np.full(
            latent_shape, vary_one_latent_from_trained_prior_samples)
        print(vary_one_latent_from_trained_prior_samples.shape)
        print(vary_one_latent_from_trained_prior_samples)

        vary_one_latent_from_trained_prior_samples[:, dim] = latent_variation
        print(vary_one_latent_from_trained_prior_samples.shape)
        print(vary_one_latent_from_trained_prior_samples)

        # decode the samples and trasform them with the PCA
        gen_vary_one_latent_from_trained_prior_samples = decoder_fn(
            vary_one_latent_from_trained_prior_samples).reshape(num_recall_samples, -1)
        reduced_gen_vary_one_latent_from_trained_prior_samples = decoded_gt_pca.transform(
            gen_vary_one_latent_from_trained_prior_samples)
        assert(reduced_gen_vary_one_latent_from_trained_prior_samples.shape == \
               reduced_generated_prior_samples.shape)

        # - Calculate the recall
        # compare to the generated normal prior samples
        vary_one_vs_prior_generated_result = iprd.knn_precision_recall_features(
            reduced_generated_prior_samples,
            reduced_gen_vary_one_latent_from_trained_prior_samples,
            nhood_sizes=nhood_sizes,
            row_batch_size=500, col_batch_size=100, num_gpus=1)
        print(agg_vary_one_vs_prior_generated_result)
        print(vary_one_vs_prior_generated_result)
        agg_vary_one_vs_prior_generated_result = agg_recall_dict(
            agg_vary_one_vs_prior_generated_result,
            vary_one_vs_prior_generated_result,
            intervention)
        print(agg_vary_one_vs_prior_generated_result)

        # compare to the generated trained prior samples
        vary_one_vs_trained_prior_generated_result = iprd.knn_precision_recall_features(
            reduced_generated_trained_prior_samples,
            reduced_gen_vary_one_latent_from_trained_prior_samples,
            nhood_sizes=nhood_sizes,
            row_batch_size=500, col_batch_size=100, num_gpus=1)
        agg_vary_one_vs_trained_prior_generated_result = agg_recall_dict(
            agg_vary_one_vs_trained_prior_generated_result,
            vary_one_vs_trained_prior_generated_result,
            intervention)

        # compare to the decoded gt
        vary_one_vs_decoded_gt_result = iprd.knn_precision_recall_features(
            reduced_decoded_gt_samples,
            reduced_gen_vary_one_latent_from_trained_prior_samples,
            nhood_sizes=nhood_sizes,
            row_batch_size=500, col_batch_size=100, num_gpus=1)
        agg_vary_one_vs_decoded_gt_result = agg_recall_dict(
            agg_vary_one_vs_decoded_gt_result,
            vary_one_vs_decoded_gt_result,
            intervention)

      update_result_dict_with_agg(
          result_d,
          [str(dim) + '_fix_one_vs_prior_generated_', agg_fix_one_vs_prior_generated_result],
          [str(dim) + '_fix_one_vs_trained_prior_generated_', agg_fix_one_vs_trained_prior_generated_result],
          [str(dim) + '_fix_one_vs_decoded_gt_', agg_fix_one_vs_decoded_gt_result],
          [str(dim) + '_vary_one_vs_prior_generated_', agg_vary_one_vs_prior_generated_result],
          [str(dim) + '_vary_one_vs_trained_prior_generated_', agg_vary_one_vs_trained_prior_generated_result],
          [str(dim) + '_vary_one_vs_decoded_gt_', agg_vary_one_vs_decoded_gt_result])
  print(result_d)
  return [result_d]

def update_result_dict(result_d, *args):
  for arg in args:
    update_key = arg[0]
    update_d = {update_key + key: list(value) for key, value in arg[1].items()}
    result_d.update(update_d)
  return result_d

def update_result_dict_with_agg(result_d, *args):
  for arg in args:
    update_key = arg[0]
    update_recall = arg[1]['recall']
    print('\n\n\n\n\n', arg)
#    update_d = {update_key + 'recall_mean': list(np.mean(value, axis=0)) for key, value in arg[1].items()}
    update_d = {update_key + 'recall_mean': list(np.mean(update_recall, axis=0)),
                update_key + 'recall_std': list(np.std(update_recall, axis=0)),
                update_key + 'recall_sum': list(np.sum(update_recall, axis=0))
                }
    result_d.update(update_d)
  return result_d

def agg_recall_dict(agg_d, new_d, inter_id):
  if inter_id == 0:
    return new_d
  else:
    for k, v in new_d.items():
      agg_d[k] = np.vstack([agg_d[k], new_d[k]])
    return agg_d


text:
from __init__ import *

import sys
import subprocess
import numpy as np
from fractions import Fraction

sys.path.insert(0, ROOT)

from compiler import *
from constructs import *

def test_fusion_pipe(pipe_data):

    R = Parameter(Int, "R")
    C = Parameter(Int, "C")
    x = Variable(Int, "x")
    y = Variable(Int, "y")

    pipe_data['R'] = R
    pipe_data['C'] = C

    row = Interval(Int, 0, R+1)
    col = Interval(Int, 0, C+1)

    cond = Condition(x, '>=', 1) & Condition(x, '<=', R) & \
           Condition(y, '<=', C) & Condition(y, '>=', 1)

    img = Image(Float, "img", [R+2, C+2])

    Iy = Function(([x, y], [row, col]), Float, "Iy")
    Iy.defn = [ Case(cond, img(x-1, y-1)*(-1.0/12) + img(x-1, y+1)*(1.0/12) + \
                           img(x, y-1)*(-2.0/12) + img(x, y+1)*(2.0/12) + \
                           img(x+1, y-1)*(-1.0/12) + img(x+1, y+1)*(1.0/12)) ]

    Ix = Function(([x, y], [row, col]), Float, "Ix")
    Ix.defn = [ Case(cond, img(x-1, y-1)*(-1.0/12) + img(x+1, y-1)*(1.0/12) + \
                           img(x-1, y)*(-2.0/12) + img(x+1, y)*(2.0/12) + \
                           img(x-1, y+1)*(-1.0/12) + img(x+1, y+1)*(1.0/12)) ]

    return [Ix, Iy]


text:
import birdvoxdetect as bvd
import os
import soundfile as sf
import numpy as np
import noisereduce as nr
import pandas as pd

def get_mono(file_path):
    """
    Returns the data in a given file after converting it to mono channel.
    """
    with sf.SoundFile(file_path, "r") as f:
        a = f.read(always_2d=True)
        sr = f.samplerate
    a = np.mean(a, -1)
    return a, sr

def get_temp_name(outdir):
    i = 0
    b = True
    temp_csv_name = "temp{}.csv".format(i)
    temp_wav_name = "temp{}.wav".format(i)
    while b:
        b = False
        while os.path.exists(os.path.join(outdir, temp_csv_name)):
            i += 1
            temp_csv_name = "temp{}.csv".format(i)
        temp_wav_name = "temp{}.wav".format(i)
        if os.path.exists(os.path.join(outdir, temp_wav_name)):
            b = True
            i += 1
    return os.path.join(outdir, temp_wav_name), os.path.join(outdir, temp_csv_name)

def run_birdvoxdetect_segmentation(
    file_names=None,
    dataset_dir=None,
    noise_reduction=False,
    noise_file=None,
    threshold=50.0,
    outdir=None,
):
    """
    Runs birdvoxdetect for getting timestamps of bird calls.

    Parameters
    ----------
    file_names : list of str, optionals
        The list of file names for input.
        If not given then all files in dataset directory are taken.

    dataset_dir : str, optional
        The dataset directory. If not given then it is assumed that file_names contains
        the complete path to the files.

    noise_reduction : boolean, optional
        If True, then noise reduction is done.
        By default set to False.

    noise_file : string, optional
        The path to the file containing the sample expected noise.
        Required if noise_reduction is True.

    threshold : float, optional
        The threshold for birdvoxdetect. A between from 10 to 90, both included.
        Higher value means greater precision at the cost of recall.
        By default set to 50.0.

    Returns
    -------
    pairs : list of tuple (str, list of str, float)
        A list, where each element is a tuple of file_path, timestamps and sample rate,
        where timestamps is a list of timestamps where a bird call is detected,
        as a str in format hh:mm:ss.ss
    """
    found_birds = []
    # giving a random name so as to prevent such a file from existing, and
    # saving it in the same folder as this file. This ensures the checklist
    # generated is also saved in this folder.
    save_path, out_path = get_temp_name(outdir)
    threshold = min(max(threshold, 10.0), 90.0)
    if file_names is None:
        assert (
            dataset_dir is not None
        ), "Need at least one of file_names and dataset_dir"
        file_names = os.listdir(dataset_dir)
    for file_name in file_names:
        file_path = file_name
        if dataset_dir is not None:
            file_path = os.path.join(dataset_dir, file_path)
        wave, samplerate = get_mono(file_path)
        with sf.SoundFile(save_path, "w", samplerate=samplerate, channels=1) as f:
            f.write(wave)
        df = bvd.process_file(save_path, threshold=threshold)
        timestamps = np.array(df["Time (hh:mm:ss)"])
        if noise_reduction:
            assert noise_file is not None, "Need noise file to apply noise reduction"
            with sf.SoundFile(noise_file, "r") as f:
                noise = f.read()
            wave = nr.reduce_noise(wave, noise)
            with sf.SoundFile(save_path, "w", samplerate=samplerate, channels=1) as f:
                f.write(wave)
            df = bvd.process_file(save_path, threshold=threshold)
            timestamps = np.concatenate((timestamps, np.array(df["Time (hh:mm:ss)"])))
            timestamps = np.sort(timestamps, kind="mergesort")
        if timestamps.shape[0] > 0:
            found_birds.append((file_name, timestamps, samplerate))
    if os.path.isfile(save_path):
        os.remove(save_path)
    if os.path.isfile(out_path):
        os.remove(out_path)
    return found_birds

max_stars_repo_name: jainris/FreeByrd

max_issues_repo_name: jainris/FreeByrd

max_forks_repo_name: jainris/FreeByrd


text:
from scipy.io import loadmat
import argparse
import pickle
from numpy.linalg import inv
import numpy as np
import pandas as pd
from scipy.spatial import distance
import math

def scipy_distance(v, u):
    return distance.euclidean(v, u)

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Reads Raw data in .mat format and converts to dictionary of class names and their attributes, saves as pickle")
    parser.add_argument("-d", "--data_set", required=True,help=("Provide the dataset name"))

    args = vars(parser.parse_args())

    #Setting location variables
    data_dir = "/home/hd71992/thesis/new_blob/data"
    data_set = "/"+ args['data_set']
    data_loc = data_dir+data_set
    out_dir = "/home/hd71992/thesis/new_blob/outputs"

    #Loading data attribute features and class names
    mat_attr = loadmat(data_loc+'/att_splits.mat')
    #mat_attr.keys()
    m_classnames = mat_attr['allclasses_names']
    m_classatt = pd.DataFrame(mat_attr['att']).T

    #Creating Classname and Attribure dictionary
    att_dict = {}
    for i in np.arange(0,len(m_classnames)):
        class_name = m_classnames[i][0][0]
        class_att = m_classatt.iloc[i,:].values
        att_dict[class_name] = class_att

    #Saving Classname and Attribure dictionary for further use
    f = open(data_loc+"/att_dict.pkl","wb")
    pickle.dump(att_dict,f)
    f.close()

    #Generating close word dictionary for prediction phase
#     dict_keys = list(att_dict.keys())
#     closeWords_Count = 6
#     closeWord_dict = {}
#
#     for word in dict_keys:
#         distance_dict = {}
#         for fast_word in dict_keys:
#             dist = scipy_distance(att_dict[word],att_dict[fast_word])
#             distance_dict[fast_word] = dist
#         closeWords_dict = {k: v for k, v in sorted(distance_dict.items(), key=lambda item: item[1])[:closeWords_Count]}
#         closeWord_dict[word] = list(closeWords_dict.keys())
#     pickle.dump(closeWord_dict, open(data_loc+'/closeWord_dict.pkl', 'wb'))

    print('Program Completed')

reason: import numpy,from numpy,from scipy


text:
from collections import OrderedDict
import numpy as np

def check_arrays(*arrays):
    """
    Left for consistency, version of `sklearn.validation.check_arrays`
    :param list[iterable] arrays: arrays with same length of first dimension.
    """
    assert len(arrays) > 0, 'The number of array must be greater than zero'
    checked_arrays = []
    shapes = []
    for arr in arrays:
        if arr is not None:
            checked_arrays.append(np.array(arr))
            shapes.append(checked_arrays[-1].shape[0])
        else:
            checked_arrays.append(None)
    assert np.sum(np.array(shapes) == shapes[0]) == len(shapes), 'Different shapes of the arrays {}'.format(
        shapes)
    return checked_arrays

def get_efficiencies(prediction, spectator, sample_weight=None, bins_number=20,
                     thresholds=None, errors=False, ignored_sideband=0.0):
    """
    Construct efficiency function dependent on spectator for each threshold
    Different score functions available: Efficiency, Precision, Recall, F1Score,
    and other things from sklearn.metrics
    :param prediction: list of probabilities
    :param spectator: list of spectator's values
    :param bins_number: int, count of bins for plot
    :param thresholds: list of prediction's threshold
        (default=prediction's cuts for which efficiency will be [0.2, 0.4, 0.5, 0.6, 0.8])
    :return:
        if errors=False
        OrderedDict threshold -> (x_values, y_values)
        if errors=True
        OrderedDict threshold -> (x_values, y_values, y_err, x_err)
        All the parts: x_values, y_values, y_err, x_err are np.arrays of the same length.
    """
    prediction, spectator, sample_weight = \
        check_arrays(prediction, spectator, sample_weight)

    spectator_min, spectator_max = weighted_quantile(spectator, [ignored_sideband, (1. - ignored_sideband)])
    mask = (spectator >= spectator_min) & (spectator <= spectator_max)
    spectator = spectator[mask]
    prediction = prediction[mask]
    bins_number = min(bins_number, len(prediction))
    sample_weight = sample_weight if sample_weight is None else np.array(sample_weight)[mask]

    if thresholds is None:
        thresholds = [weighted_quantile(prediction, quantiles=1 - eff, sample_weight=sample_weight)
                      for eff in [0.2, 0.4, 0.5, 0.6, 0.8]]

    binner = Binner(spectator, bins_number=bins_number)
    if sample_weight is None:
        sample_weight = np.ones(len(prediction))
    bins_data = binner.split_into_bins(spectator, prediction, sample_weight)

    bin_edges = np.array([spectator_min] + list(binner.limits) + [spectator_max])
    xerr = np.diff(bin_edges) / 2.
    result = OrderedDict()
    for threshold in thresholds:
        x_values = []
        y_values = []
        N_in_bin = []
        for num, (masses, probabilities, weights) in enumerate(bins_data):
            y_values.append(np.average(probabilities > threshold, weights=weights))
            N_in_bin.append(np.sum(weights))
            if errors:
                x_values.append((bin_edges[num + 1] + bin_edges[num]) / 2.)
            else:
                x_values.append(np.mean(masses))

        x_values, y_values, N_in_bin = check_arrays(x_values, y_values, N_in_bin)
        if errors:
            result[threshold] = (x_values, y_values, np.sqrt(y_values * (1 - y_values) / N_in_bin), xerr)
        else:
            result[threshold] = (x_values, y_values)
    return result

def weighted_quantile(array, quantiles, sample_weight=None, array_sorted=False, old_style=False):
    """Computing quantiles of array. Unlike the np.percentile, this function supports weights,
    but it is inefficient and performs complete sorting.
    :param array: distribution, array of shape [n_samples]
    :param quantiles: floats from range [0, 1] with quantiles of shape [n_quantiles]
    :param sample_weight: optional weights of samples, array of shape [n_samples]
    :param array_sorted: if True, the sorting step will be skipped
    :param old_style: if True, will correct output to be consistent with np.percentile.
    :return: array of shape [n_quantiles]
    Example:
    >>> weighted_quantile([1, 2, 3, 4, 5], [0.5])
    Out: array([ 3.])
    >>> weighted_quantile([1, 2, 3, 4, 5], [0.5], sample_weight=[3, 1, 1, 1, 1])
    Out: array([ 2.])
    """
    array = np.array(array)
    quantiles = np.array(quantiles)
    sample_weight = check_sample_weight(array, sample_weight)
    assert np.all(quantiles >= 0) and np.all(quantiles <= 1), 'Percentiles should be in [0, 1]'

    if not array_sorted:
        array, sample_weight = reorder_by_first(array, sample_weight)

    weighted_quantiles = np.cumsum(sample_weight) - 0.5 * sample_weight
    if old_style:
        # To be convenient with np.percentile
        weighted_quantiles -= weighted_quantiles[0]
        weighted_quantiles /= weighted_quantiles[-1]
    else:
        weighted_quantiles /= np.sum(sample_weight)
    return np.interp(quantiles, weighted_quantiles, array)

def check_sample_weight(y_true, sample_weight):
    """Checks the weights, if None, returns array.
    :param y_true: labels (or any array of length [n_samples])
    :param sample_weight: None or array of length [n_samples]
    :return: np.array of shape [n_samples]
    """
    if sample_weight is None:
        return np.ones(len(y_true), dtype=np.float)
    else:
        sample_weight = np.array(sample_weight, dtype=np.float)
        assert len(y_true) == len(sample_weight), \
            "The length of weights is different: not {0}, but {1}".format(len(y_true), len(sample_weight))
        return sample_weight

def reorder_by_first(*arrays):
    """
    Applies the same permutation to all passed arrays,
    permutation sorts the first passed array
    """
    arrays = check_arrays(*arrays)
    order = np.argsort(arrays[0])
    return [arr[order] for arr in arrays]

class Binner(object):
    def __init__(self, values, bins_number):
        """
        Binner is a class that helps to split the values into several bins.
        Initially an array of values is given, which is then splitted into 'bins_number' equal parts,
        and thus we are computing limits (boundaries of bins).
        """
        percentiles = [i * 100.0 / bins_number for i in range(1, bins_number)]
        self.limits = np.percentile(values, percentiles)

    def get_bins(self, values):
        """Given the values of feature, compute the index of bin
        :param values: array of shape [n_samples]
        :return: array of shape [n_samples]
        """
        return np.searchsorted(self.limits, values)

    def set_limits(self, limits):
        """Change the thresholds inside bins."""
        self.limits = limits

    @property
    def bins_number(self):
        """:return: number of bins"""
        return len(self.limits) + 1

    def split_into_bins(self, *arrays):
        """
        :param arrays: data to be splitted, the first array corresponds
        :return: sequence of length [n_bins] with values corresponding to each bin.
        """
        values = arrays[0]
        for array in arrays:
            assert len(array) == len(values), "passed arrays have different length"
        bins = self.get_bins(values)
        result = []
        for bin in range(len(self.limits) + 1):
            indices = bins == bin
            result.append([np.array(array)[indices] for array in arrays])
        return result


text:
# -------------------------------------------------------------------
# Full moment tensor
#
# refs
# Tape, W., & Tape, C. (2015). A uniform parametrization of moment tensors.
# Geophysical Journal International, 202(3), 2074–2081.
#
# Author: Liang Ding
# Email: myliang.ding@mail.utoronto.ca
# -------------------------------------------------------------------

from seisgen.MTTools.DLune import DCreate_Triple
from seisgen.MTTools.DOrientations import DRotationMatrix, DRotz
import numpy as np

gRMax = DRotz(np.deg2rad(-45))
ginvRMax = np.linalg.inv(gRMax)

def DMomentTensor(strike, dip, rake, colatitude, longitude):
    '''
    Calculate the FULL moment tensor for the given strike, dip, rake and the colatitude, longitude in Lune.

    :param strike:  The strike in radians, k in Carl's paper.  ranging [0, 2pi]
    :param dip:     The dip in radians, thetha in Carl's paper. ranging [0, pi/2]
    :param rake:    The slip (rake) in radians, sigma in Carl's paper. ranging [-pi/2, pi/2]
    :param colatitude: The colatitude in the Lune in radians (Beta), [0, PI]
    :param longitude:  The longitude in the Lune in radians (Gamma), [-PI/6, PI/6]
    :return: The full moment tensor.
    '''

    triple = DCreate_Triple(colatitude, longitude)
    triple_matrix = np.zeros((3, 3))
    triple_matrix[0, 0] = triple[0]
    triple_matrix[1, 1] = triple[1]
    triple_matrix[2, 2] = triple[2]

    U = DRotationMatrix(strike, dip, rake)
    return np.dot(np.dot(U, triple_matrix), np.linalg.inv(U))

def DMT_enz(strike, dip, rake, colatitude, longitude):
    ''' return the moment tensor in e-n-z coord.'''

    mt_nwz = DMomentTensor(strike, dip, rake, colatitude, longitude)
    mt_enz = np.ones(6)
    mt_enz[0] = mt_nwz[1][1]
    mt_enz[1] = mt_nwz[0][0]
    mt_enz[2] = mt_nwz[2][2]

    mt_enz[3] = -1.0 * mt_nwz[0][1]
    mt_enz[4] = -1.0 * mt_nwz[1][2]
    mt_enz[5] = mt_nwz[0][2]

    return mt_enz

max_stars_repo_name: Liang-Ding/seisgen

max_issues_repo_name: Liang-Ding/seisgen

max_forks_repo_name: Liang-Ding/seisgen