text:
import pickle
import numpy as np

pad = '<pad>'
eos = '<eos>'
mask = [f'm_{num}' for num in range(1)]

special_tokens = [pad, eos]

time_signature_token = ['4/4', '3/4', '2/4', '6/8']

program_num = [f'i_{num}' for num in range(128)]

tempo_token = [f't_{i}' for i in range(7)]

track_num = [f'track_{num}' for num in range(3)]

structure_token = ['bar'] + track_num

song_token = time_signature_token + tempo_token + program_num

rests = ['rest_e', 'rest_s']
step_token = [f'e_{num}' for num in range(16)]
duration_multi = ['whole', 'half', 'quarter', 'eighth', 'sixteenth']
duration_single = [f'n_{num}' for num in range(1,33)]

pitch_tokens = [f'p_{num}' for num in range(21, 109)]

duration_to_time = {'whole':4,'half':2,'quarter':1,'eighth':0.5,'sixteenth':0.25}

all_key_names = ['C major', 'G major', 'D major', 'A major',
                 'E major', 'B major', 'F major', 'B- major',
                 'E- major', 'A- major', 'D- major', 'G- major',
                 'A minor', 'E minor', 'B minor', 'F# minor',
                 'C# minor', 'G# minor', 'D minor', 'G minor',
                 'C minor', 'F minor', 'B- minor', 'E- minor',
                 ]

all_major_names = np.array(['C major', 'D- major', 'D major', 'E- major',
                   'E major', 'F major', 'G- major', 'G major',
                   'A- major', 'A major', 'B- major', 'B major'])

all_minor_names = np.array(['A minor', 'B- minor', 'B minor', 'C minor',
                   'C# minor', 'D minor', 'E- minor', 'E minor',
                   'F minor', 'F# minor', 'G minor', 'G# minor'])

key_token = [f'k_{num}' for num in range(len(all_key_names))]
key_to_token = {name: f'k_{i}' for i, name in enumerate(all_key_names)}
token_to_key = {v: k for k, v in key_to_token.items()}

track_note_density_token = [f'd_{num}' for num in range(10)]
track_occupation_rate_token = [f'o_{num}' for num in range(10)]
track_polyphony_rate_token = [f'y_{num}' for num in range(10)]
track_pitch_register_token = [f'r_{num}' for num in range(8)]

tensile_strain_token = [f's_{num}' for num in range(12)]
diameter_token = [f'a_{num}' for num in range(12)]
tension_token = [f'l_{num}' for num in range(12)]

control_bins = np.arange(0, 1, 0.1)
tensile_bins = np.arange(0, 2.1, 0.2).tolist() + [4]
diameter_bins = np.arange(0, 4.1, 0.4).tolist() + [5]

tempo_bins = np.array([0] + list(range(60, 190, 30)) + [200])
tension_bin = np.arange(0,6.5,0.5)
tension_bin[-1] = 6.5

class WordVocab(object):
    def __init__(self, mode,control_list):
        super(WordVocab, self).__init__()
        #control list is []
        if mode == 0:
            # SMER
            duration_only_token = duration_multi
            duration_tokens = duration_only_token + rests + ['continue']

        else:
            # REMI
            duration_only_token = duration_single
            duration_tokens = step_token + duration_only_token

        note_tokens = pitch_tokens + duration_tokens

        basic_tokens = special_tokens + mask + structure_token + \
                     song_token + note_tokens

        all_tokens = basic_tokens + \
                     track_note_density_token + \
                     track_polyphony_rate_token + \
                     track_occupation_rate_token + \
                     key_token + tensile_strain_token + \
                     diameter_token

        self.pad_index = 0
        self.eos_index = 1
        self.char_lst = all_tokens
        self.basic_tokens = basic_tokens

        self._char2idx = {
            '<pad>': self.pad_index,
            '<eos>': self.eos_index,
        }

        for char in self.char_lst:
            if char not in self._char2idx:
                self._char2idx[char] = len(self._char2idx)
        self._idx2char = dict((idx, char) for char, idx in self._char2idx.items())
        print(f'vocab size: {self.vocab_size}')

        self.token_class_ranges = {}
        self.name_to_tokens = {}
        self.structure_indices = [self._char2idx[name] for name in structure_token]
        self.pitch_indices = [self._char2idx[name] for name in pitch_tokens]
        self.mask_indices = [self._char2idx[name] for name in mask]
        self.duration_indices = [self._char2idx[name] for name in duration_tokens]
        self.duration_only_indices = [self._char2idx[name] for name in duration_only_token]
        self.program_indices = [self._char2idx[name] for name in program_num]
        self.tempo_indices = [self._char2idx[name] for name in tempo_token]
        self.time_signature_indices = [self._char2idx[name] for name in time_signature_token]
        self.rest_indices = []
        self.control_indices = {}
        self.control_tokens = []

        if mode == 0:
            self.rest_indices = [self._char2idx[name] for name in rests]
            self.continue_index = self._char2idx['continue']
        if mode == 1:
            self.step_indices = [self._char2idx[name] for name in step_token]

        for index in self.program_indices:
            self.token_class_ranges[index] = 'program'
            if 'program' in self.name_to_tokens:
                self.name_to_tokens['program'].append(self._idx2char[index])
            else:
                self.name_to_tokens['program'] = [self._idx2char[index]]

        for index in self.rest_indices:
            self.token_class_ranges[index] = 'rests'
            if 'rests' in self.name_to_tokens:
                self.name_to_tokens['rests'].append(self._idx2char[index])
            else:
                self.name_to_tokens['rests'] = [self._idx2char[index]]
        for index in self.tempo_indices:
            self.token_class_ranges[index] = 'tempo'
            if 'tempo' in self.name_to_tokens:
                self.name_to_tokens['tempo'].append(self._idx2char[index])
            else:
                self.name_to_tokens['tempo'] = [self._idx2char[index]]
        for index in self.time_signature_indices:
            self.token_class_ranges[index] = 'time_signature'
            if 'time_signature' in self.name_to_tokens:
                self.name_to_tokens['time_signature'].append(self._idx2char[index])
            else:
                self.name_to_tokens['time_signature'] = [self._idx2char[index]]

        for index in self.structure_indices:
            self.token_class_ranges[index] = 'structure'
            if 'structure' in self.name_to_tokens:
                self.name_to_tokens['structure'].append(self._idx2char[index])
            else:
                self.name_to_tokens['structure'] = [self._idx2char[index]]
        for index in self.pitch_indices:
            self.token_class_ranges[index] = 'pitch'
            if 'pitch' in self.name_to_tokens:
                self.name_to_tokens['pitch'].append(self._idx2char[index])
            else:
                self.name_to_tokens['pitch'] = [self._idx2char[index]]
        for index in self.duration_indices:
            self.token_class_ranges[index] = 'duration'
            if 'duration' in self.name_to_tokens:
                self.name_to_tokens['duration'].append(self._idx2char[index])
            else:
                self.name_to_tokens['duration'] = [self._idx2char[index]]

        self.token_class_ranges[self.eos_index] = 'eos'

        self.name_to_tokens['eos'] = self._idx2char[self.eos_index]

        self.class_names = set(self.token_class_ranges.values())

    def char2index(self, token):
        # print(token)
        if token not in self._char2idx:
            print('invalid')

        return self._char2idx.get(token)

    def index2char(self, idxs):

        return self._idx2char.get(idxs)

    def get_token_classes(self, idx):

        return self.token_class_ranges[idx]

    @property
    def vocab_size(self):
        return len(self._char2idx)

    def save_vocab(self, vocab_path):
        with open(vocab_path, "wb") as f:
            pickle.dump(self, f)

    @staticmethod
    def load_vocab(vocab_path: str) -> 'WordVocab':
        with open(vocab_path, "rb") as f:
            return pickle.load(f)

max_stars_repo_name: ruiguo-bio/smer

max_issues_repo_name: ruiguo-bio/smer

max_forks_repo_name: ruiguo-bio/smer


text:
"""
Analog data acquisition for QuSpin's OPMs via National Instruments' cDAQ unit
The following assumes:
"""

# Imports
import matplotlib.pyplot as plt
import numpy as np

import nidaqmx
from nidaqmx.stream_readers import AnalogMultiChannelReader
from nidaqmx import constants
# from nidaqmx import stream_readers  # not needed in this script
# from nidaqmx import stream_writers  # not needed in this script

import threading
import pickle
from datetime import datetime
import scipy.io

# Parameters
sampling_freq_in = 1000  # in Hz
buffer_in_size = 100
bufsize_callback = buffer_in_size
buffer_in_size_cfg = round(buffer_in_size * 1)  # clock configuration
chans_in = 1  # set to number of active OPMs (x2 if By and Bz are used, but that is not recommended)
refresh_rate_plot = 10  # in Hz
crop = 10  # number of seconds to drop at acquisition start before saving
my_filename = 'test_3_opms'  # with full path if target folder different from current folder (do not leave trailing /)

# Initialize data placeholders
buffer_in = np.zeros((chans_in, buffer_in_size))
data = np.zeros((chans_in, 1))  # will contain a first column with zeros but that's fine

# Definitions of basic functions
def ask_user():
    global running
    input("Press ENTER/RETURN to stop acquisition and coil drivers.")
    running = False

def cfg_read_task(acquisition):  # uses above parameters
    acquisition.ai_channels.add_ai_voltage_chan("Dev1/ai1:3")  # has to match with chans_in
    acquisition.timing.cfg_samp_clk_timing(rate=sampling_freq_in, sample_mode=constants.AcquisitionType.CONTINUOUS,
                                           samps_per_chan=buffer_in_size_cfg)

def reading_task_callback(task_idx, event_type, num_samples, callback_data):  # bufsize_callback is passed to num_samples
    global data
    global buffer_in

    if running:
        # It may be wiser to read slightly more than num_samples here, to make sure one does not miss any sample,
        # see: https://documentation.help/NI-DAQmx-Key-Concepts/contCAcqGen.html
        buffer_in = np.zeros((chans_in, num_samples))  # double definition ???
        stream_in.read_many_sample(buffer_in, num_samples, timeout=constants.WAIT_INFINITELY)

        data = np.append(data, buffer_in, axis=1)  # appends buffered data to total variable data

    return 0  # Absolutely needed for this callback to be well defined (see nidaqmx doc).

# Configure and setup the tasks
task_in = nidaqmx.Task()
cfg_read_task(task_in)
stream_in = AnalogMultiChannelReader(task_in.in_stream)
task_in.register_every_n_samples_acquired_into_buffer_event(bufsize_callback, reading_task_callback)

# Start threading to prompt user to stop
thread_user = threading.Thread(target=ask_user)
thread_user.start()

# Main loop
running = True
time_start = datetime.now()
task_in.start()

# Plot a visual feedback for the user's mental health
f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex='all', sharey='none')
while running:  # make this adapt to number of channels automatically
    ax1.clear()
    ax2.clear()
    ax3.clear()
    ax1.plot(data[0, -sampling_freq_in * 5:].T)  # 5 seconds rolling window
    ax2.plot(data[1, -sampling_freq_in * 5:].T)
    ax3.plot(data[2, -sampling_freq_in * 5:].T)
    # Label and axis formatting
    ax3.set_xlabel('time [s]')
    ax1.set_ylabel('voltage [V]')
    ax2.set_ylabel('voltage [V]')
    ax3.set_ylabel('voltage [V]')
    xticks = np.arange(0, data[0, -sampling_freq_in * 5:].size, sampling_freq_in)
    xticklabels = np.arange(0, xticks.size, 1)
    ax3.set_xticks(xticks)
    ax3.set_xticklabels(xticklabels)

    plt.pause(1/refresh_rate_plot)  # required for dynamic plot to work (if too low, nulling performance bad)

# Close task to clear connection once done
task_in.close()
duration = datetime.now() - time_start

# Final save data and metadata ... first in python reloadable format:
filename = my_filename
with open(filename, 'wb') as f:
    pickle.dump(data, f)
'''
Load this variable back with:
with open(name, 'rb') as f:
    data_reloaded = pickle.load(f)
'''
# Human-readable text file:
extension = '.txt'
np.set_printoptions(threshold=np.inf, linewidth=np.inf)  # turn off summarization, line-wrapping
with open(filename + extension, 'w') as f:
    f.write(np.array2string(data.T, separator=', '))  # improve precision here!
# Now in matlab:
extension = '.mat'
scipy.io.savemat(filename + extension, {'data':data})

# Some messages at the end
num_samples_acquired = data[0,:].size
print("\n")
print("OPM acquisition ended.\n")
print("Acquisition duration: {}.".format(duration))
print("Acquired samples: {}.".format(num_samples_acquired - 1))

# Final plot of whole time course the acquisition
plt.close('all')
f_tot, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex='all', sharey='none')
ax1.plot(data[0, 10:].T)  # note the exclusion of the first 10 iterations (automatically zoomed in plot)
ax2.plot(data[1, 10:].T)
ax3.plot(data[2, 10:].T)
# Label formatting ...
ax3.set_xlabel('time [s]')
ax1.set_ylabel('voltage [V]')
ax2.set_ylabel('voltage [V]')
ax3.set_ylabel('voltage [V]')
xticks = np.arange(0, data[0, :].size, sampling_freq_in)
xticklabels = np.arange(0, xticks.size, 1)
ax3.set_xticks(xticks)
ax3.set_xticklabels(xticklabels)
plt.show()

ext: py

max_stars_repo_name: VasoTracker/VasoTracker

max_issues_repo_name: VasoTracker/VasoTracker

max_forks_repo_name: VasoTracker/VasoTracker

reason: import numpy,import scipy


text:
def data():
    import pickle

    with open('embedding_labels.pkl', 'rb') as lf:
        labels = pickle.load(lf)

    with open('embedding_posts.pkl', 'rb') as pf:
        posts = pickle.load(pf)

    len(labels)

    # Tokenizing data
    from keras.preprocessing.text import Tokenizer
    from keras.preprocessing import sequence

    max_features = 10000
    max_len = 50
    batch_size = 32
    train_samples = 329885 # Roughly about 70% of Data
    max_words = 10000

    tokenizer = Tokenizer(num_words=max_words)
    tokenizer.fit_on_texts(posts)
    sequences = tokenizer.texts_to_sequences(posts)
    word_index = tokenizer.word_index
    print('Found %s unique tokens' % len(word_index))

    data = sequence.pad_sequences(sequences, maxlen=max_len)
    label_dictionary = {
        'I': 0,
        'E': 1
    }

    labels = [label_dictionary[label[:1]] for label in labels]

    import numpy as np
    labels = np.asarray(labels)

    indices = np.arange(data.shape[0])
    np.random.shuffle(indices)
    data = data[indices]
    labels = labels[indices]

    x_test = data[train_samples:]
    y_test = labels[train_samples:]

    x_train = data[:train_samples]
    y_train = labels[:train_samples]

    return x_train, y_train, x_test, y_test

class_weights = {
    0:1.0,
    1:3.27
}

from keras import models
from keras import layers
from keras.regularizers import l1_l2
from keras.layers import Embedding, Dropout, LSTM, Dense

from hyperas import optim
from hyperas.distributions import choice, uniform

def create_model(x_train, y_train, x_test, y_test):
    model = models.Sequential()
    model.add(Embedding(max_features, max_len))
    model.add(LSTM({{choice([8,16,32])}}, return_sequences=True))
    model.add(LSTM({{choice([8,16,32])}}))
    if {{choice(['three','two'])}} == 'three':
        model.add(Dense({{choice([8,16,32])}}))
    else:
        pass

    model.add(Dense(1, activation='sigmoid'))

    # model.layers[0].set_weights([embedding_matrix])
    # model.layers[0].trainable = False

    model.compile(optimizer={{choice(['rmsprop','adam','sgd'])}},
                  loss='binary_crossentropy',
                  metrics=['accuracy'])

    model.fit(x_train, y_train, batch_size={{choice([64, 128])}},
             epochs={{uniform(10, 20)}}, validation_data=(x_test, y_test),
             class_weight=class_weights)

    score, acc = model.evaluate(x_test_, y_test)
    print('Test accuracy:', acc)
    return {'loss':-acc, 'status': STATUS_OK, 'model': model}

from hyperopt import Trials, STATUS_OK, tpe

if __name__ == '__main__':
    best_run, best_model = optim.minimize(model=create_model,
                                          data=data,
                                          algo=tpe.suggest,
                                          max_evals=30,
                                          trials=Trials())
    X_train, Y_train, X_test, Y_test = data()
    print("Evaluation of best performing model:")
    print(best_model.evaluate(X_test, Y_test))
    print("Best performing model chosen hyper-parameters:")
    print(best_run)

max_stars_repo_name: PaulKMandal/MBTI-Prediction

max_issues_repo_name: PaulKMandal/MBTI-Prediction

max_forks_repo_name: PaulKMandal/MBTI-Prediction

reason: import numpy


text:
"""

.. module:: statistics
	:platform: Unix
	:synopsis: This module implements a set of statistical operations on ensembles of weak lensing maps (shear/convergence)

.. moduleauthor:: Andrea Petri <apetri@phys.columbia.edu>

"""

from __future__ import division

from operator import add
from functools import reduce

from .index import Indexer
from ..utils import pcaHandler as PCA

import numpy as np
import scipy.io as sio
from scipy import sparse

from emcee.ensemble import _function_wrapper

try:
	import pandas as pd
	pd = pd
except ImportError:
	pd = None

##########################################################################
#########Useful for bootstrap estimates of the covariance matrix###########
##########################################################################

def _bsp_covariance(data):
	sub = data - data.mean(0)[None]
	return np.dot(sub.T,sub) / (data.shape[0] - 1.0)

##########################################################################
#########Useful for bootstrap estimates of the covariance matrix###########
##########################################################################

def _bsp_covariance(data):
	sub = data - data.mean(0)[None]
	return np.dot(sub.T,sub) / (data.shape[0] - 1.0)

##########################################
########Ensemble class####################
##########################################

class Ensemble(object):

	"""
	A class that handles statistical operations on weak lensing maps; an ensemble is a collection of different statistical realization of the same random variable. This class has an attribute 'data' that is a numpy array which first axis corresponds to the realization number.

	>>> from lenstools.statistics import Ensemble

	"""

	def __init__(self,file_list=list(),data=None,num_realizations=0,metric="chi2"):

		self.file_list = file_list
		self.data = data
		self.num_realizations = num_realizations
		self.metric = metric

	####################################
	#############I/O####################
	####################################

	@classmethod
	def fromfilelist(cls,file_list):

		"""
		Builds the ensemble from a file list: each file corresponds to a different realization

		:param file_list: List of files on which to define the ensemble
		:type file_list: list of str.

		"""

		#See how many realizations are there in the ensemble
		num_realizations = len(file_list)
		assert num_realizations>0,"There are no realizations in your ensemble!!"

		#Build the ensemble instance and return it
		return cls(file_list=file_list,num_realizations=num_realizations)

	@classmethod
	def fromdata(cls,npy_data):

		"""
		Builds the ensemble from data in numpy array format

		:param npy_data: numpy array with the data, the first dimension must be the number of realizations
		:type npy_data: array

		"""

		return cls(num_realizations=npy_data.shape[0],data=npy_data)

	@classmethod
	def read(cls,filename,callback_loader=None,**kwargs):

		"""
		Reads a numpy file into an Ensemble

		:param filename: name of the file to read
		:type filename: str.

		:param callback_loader: This function gets executed on each of the files in the list and populates the ensemble. If None provided, it performs a numpy.load on the specified file. Must return a numpy array with the loaded data
		:type callback_loader: function

		:param kwargs: Any additional keyword arguments to be passed to callback_loader
		:type kwargs: dict.

		:returns: Ensemble instance read from the file

		"""

		new_ensemble = cls.fromfilelist([filename])
		new_ensemble.load(from_old=True,callback_loader=callback_loader,**kwargs)
		return new_ensemble

	@classmethod
	def readall(cls,filelist,callback_loader=None,**kwargs):

		"""
		Reads a list of files into an Ensemble

		:param filelist: list of files to read
		:type filelist: list.

		:param callback_loader: This function gets executed on each of the files in the list and populates the ensemble. If None provided, it performs a numpy.load on the specified file. Must return a numpy array with the loaded data
		:type callback_loader: function

		:param kwargs: Any additional keyword arguments to be passed to callback_loader
		:type kwargs: dict.

		:returns: Ensemble instance read from the file

		"""

		return reduce(add,[ cls.read(f,callback_loader,**kwargs) for f in filelist ])

	def load(self,callback_loader=None,pool=None,from_old=False,**kwargs):
		"""
		Loads the ensemble into memory, can spread the calculations on multiple processors using a MPI pool

		:param callback_loader: This function gets executed on each of the files in the list and populates the ensemble. If None provided, it performs a numpy.load on the specified files. Must return a numpy array with the loaded data
		:type callback_loader: function

		:param pool: MPI pool for multiprocessing (imported from emcee https://github.com/dfm/emcee)
		:type pool: MPI pool object

		:param from_old: If True, the loaded data are interpreted as an old, already existing ensemble, which means that only one file (in which the old ensemble is saved) is loaded, the first dimension of the data is 1 and hence it is discarded
		:type from_old: bool.

		:param kwargs: Any additional keyword arguments to be passed to callback_loader
		:type kwargs: dict.

		>>> from lenstools import Ensemble
		>>> from lenstools.statistics import default_callback_loader

		>>> map_list = ["conv1.fit","conv2.fit","conv3.fit"]
		>>> l_edges = np.arange(200.0,50000.0,200.0)

		>>> conv_ensemble = Ensemble.fromfilelist(map_list)
		>>> conv_ensemble.load(callback_loader=default_callback_loader,pool=pool,l_edges=l_edges)

		"""

		if callback_loader is None:
			callback_loader = lambda f: np.load(f)

		self.pool = pool

		#Build a function wrapper of the callback loader, so it becomes pickleable
		_callback_wrapper = _function_wrapper(callback_loader,args=tuple(),kwargs=kwargs)

		#Execute the callback on each file in the list (spread calculations with MPI pool if it is not none)
		if pool is not None:
			M = pool.map
		else:
			M = map

		full_data = np.array(M(_callback_wrapper,self.file_list))

		assert type(full_data) == np.ndarray
		assert full_data.shape[0] == self.num_realizations

		if from_old:
			full_data = full_data[0]

		self.num_realizations = full_data.shape[0]
		self.data = full_data

	def save(self,filename,format=None,**kwargs):

		"""
		Save ensemble data in an external file (in arbitrary format)

		:param filename: file name of the external file
		:type filename: str.

		:format: format in which to save the ensemble; if None the format is auto detected from the filename
		:type format: str.or callable

		:param kwargs: the keyword arguments are passed to the saver (or to format if callable)
		:type kwargs: dict.

		"""

		#Auto detect format
		if format is None:
			if filename.endswith(".npy"):
				format = "numpy"
			elif filename.endswith(".mat"):
				format = "matlab"
			else:
				raise ValueError("Format not recognized!")

		#Proceed to the saving procedure
		if format=="numpy":
			np.save(filename,self.data)
		elif format=="matlab":
			sio.savemat(filename,{"data": self.data},**kwargs)
		else:
			format(self,filename,**kwargs)

	#Convert into pandas DataFrame
	def toPandas(self,column_label=None):

		"""
		Convert the Ensemble into a pandas DataFrame, interpreting the first dimension as observation and the second dimension as field

		:param column_label: columns of the DataFrame
		:type column_label: list. or Index

		:returns: DataFrame

		"""

		#Check if pandas is installed
		if pd is None:
			raise ImportError("pandas needs to be installed to use this feature!")

		#Convert into DataFrame
		row_index = pd.Index(np.arange(self.num_realizations),name="realization")
		return pd.DataFrame(self.data,index=row_index,columns=column_label)

	####################################
	#############Operations#############
	####################################

	def mean(self):

		"""
		Computes the ensemble average over realizations

		:returns: ndarray with the averages, has the same shape as self.data[0]

		"""

		self._mean = self.data.mean(axis=0)

		return self._mean

	def scale(self,weights):

		"""
		Set manually the units of each row of the ensemble by multiplying it by a row of weights

		:param weights: row of weights used to scale the ensemble, must have the same length as the number of rows
		:type weights: array

		"""

		if isinstance(weights,np.ndarray):

			assert len(weights)==self.data.shape[0]
			self.data *= weights.reshape(self.data.shape)

		else:

			self.data *= weights

	def group(self,group_size,kind="sparse",inplace=True):

		"""
		Sometimes it happens that different realizations in the ensemble need to be grouped together, for example when they belong to different subfields of the same observation field. With this function you can group different realizations together by taking the mean, and reduce the total number of realizations in the ensemble

		:param group_size: how many realizations to put in a group, must divide exactly the total number of realizations in the ensemble
		:type group_size: int

		:param kind: specifies how to do the grouping; if set to "sparse" the groups are formed by taking one realizations every num_realizations/group_size (for example ([1,1001,...,9001],[2,1002,...,9002]) if num_realizations=10000 and group_size=10). If set to "contiguous" then the realizations are grouped as ([1,2,...,10],[11,12,...,20]). Otherwise you can set kind to your own sparse matrix scheme
		:type kind: str. or scipy.sparse

		"""

		assert not(self.num_realizations%group_size),"The group size must divide exactly the number of realizations!"
		num_groups = self.num_realizations//group_size

		#Build the appropriate grouping scheme
		if isinstance(kind,sparse.csr.csr_matrix):

			assert kind.dtype == np.bool,"The scheme type must be bool, only 0 and 1!"
			scheme = kind

		elif kind=="sparse":

			scheme = reduce(add,[ sparse.eye(num_groups,self.num_realizations,k=num_groups*i,dtype=np.int8) for i in range(group_size)])

		elif kind=="contiguous":

			row = np.array(reduce(add,[ (i,) * group_size  for i in range(num_groups) ]))
			col = np.arange(self.num_realizations,dtype=np.int)
			dat = np.ones(self.num_realizations,dtype=np.int8)

			scheme = sparse.csr_matrix((dat,(row,col)),shape=(num_groups,self.num_realizations),dtype=np.int8)

		else:
			raise TypeError("The scheme kind you inputed is not valid!")

		self._scheme = scheme

		if inplace:

			#Dot the data with the scheme to produce the groups, and take the mean in every group
			self.num_realizations = num_groups
			self.data = scheme.dot(self.data) / group_size

		else:

			#Dot the data with the scheme to produce the groups, and take the mean in every group
			return self.__class__.fromdata(scheme.dot(self.data) / group_size)

	def differentiate(self,step=None,order=1):

		"""
		Compute a new ensemble, in which the feature is differentiated with respect to its label (e.g. the differentiation of the first minkowski functional is the PDF)

		:param step: measure unit on the feature label axis
		:type step: float or array

		:returns: new Ensemble instance with the differentiated feature

		"""

		#differentiate
		diff_data = np.diff(self.data,n=order,axis=1)

		#Optionally scale the measure units
		if step is not None:
			diff_data /= step

		#return the new ensemble
		return self.__class__(data=diff_data,num_realizations=diff_data.shape[0])

	def cut(self,min=None,max=None,feature_label=None,inplace=True):

		"""
		Allows to manually cut the ensemble along the second axis, if you want to select a feature subset; you better know what you are doing if you are using this function

		:param min: left extreme of the cut, included; if a list of indices is passed, the cut is performed on those indices, on the second axis and the remaining parameters are ignored
		:type min: int or float

		:param max: right extreme of the cut, included
		:type max: int or float

		:param feature_label: if not None, serves as a reference for min and max, in which the ensemble is cut according to the position of the min and max elements in the feature_label array
		:type feature_label: array

		:param inplace: if True, the Ensemble cut is made in place, otherwise a new ensemble is returned
		:type inplace: bool.

		:returns: the min and max cut indices if feature_label is None, otherwise it returns the array of the cut feature; if min was a list of indices, these are returned back. Returns the new Ensemble if inplace is False

		"""

		#Sanity checks
		assert self.data.ndim==2,"Only one dimensional feature cuts implemented so far!"
		assert (min is not None) or (max is not None),"No cutting extremes selected!"

		if type(min)==list:

			new_data = self.data[:,min]

			if inplace:

				self.data = new_data

				#Return
				return min

			else:
				return self.__class__.fromdata(new_data)

		else:

			if feature_label is not None:

				#Look for the corresponding indices in the feature label
				min_idx = np.abs(feature_label-min).argmin()
				max_idx = np.abs(feature_label-max).argmin()

			else:

				#Interpret max and min as indices between which to perform the cut
				min_idx = min
				max_idx = max

			new_data = self.data[:,min_idx:max_idx+1]

			if inplace:

				self.data = new_data

				#Return
				if feature_label is not None:
					return feature_label[min_idx:max_idx+1]
				else:
					return min_idx,max_idx

			else:
				return self.__class__.fromdata(new_data)

	def subset(self,realizations):

		"""
		Returns a new ensemble that contains only the selected realizations

		:param realizations: realizations to keep in the subset Ensemble, must be compatible with numpy array indexing syntax
		:type realizations: int. or array of int.

		:returns: new Ensemble instance that contains only the selected realizations

		"""

		data = self.data[realizations]

		return self.__class__.fromdata(data)

	def transform(self,transformation,inplace=False,**kwargs):

		"""
		Allows a general transformation on the Ensemble by calling a function on its data

		:param transformation: callback function called on the ensemble data
		:type transformation: callable

		:param inplace: if True the transformation is performed in place, otherwise a new Ensemble is created
		:type inplace: bool.

		:param kwargs: the keyword arguments are passed to the transformation callable
		:type kwargs: dict.

		:returns: new Ensemble instance if the transformation is nor performed in place

		"""

		#Apply the transformation
		transformed_data = transformation(self.data,**kwargs)

		#Return the new Ensemble
		if inplace:

			self.data = transformed_data
			self.num_realizations = transformed_data.shape[0]

		else:
			return self.__class__.fromdata(transformed_data)

	def covariance(self,bootstrap=False,**kwargs):

		"""
		Computes the ensemble covariance matrix

		:param bootstrap: if True the covariance matrix is computed with a bootstrap estimate
		:type bootstrap: bool.

		:param kwargs: the keyword arguments are passed to the bootstrap method
		:type kwargs: dict.

		:returns: ndarray with the covariance matrix, has shape (self.data[1],self.data[1])

		"""

		assert self.num_realizations>1, "I can't compute a covariance matrix with one realization only!!"
		assert self.data.dtype == np.float, "This operation is unsafe with non float numbers!!"

		if bootstrap:
			return self.bootstrap(_bsp_covariance,**kwargs)

		else:
			if not hasattr(self,"_mean"):
				self.mean()

			subtracted = self.data - self._mean[np.newaxis,:]
			return np.dot(subtracted.transpose(),subtracted) / (self.num_realizations - 1.0)

	def correlation(self):

		"""
		Computes the ensemble correlation matrix

		:returns: ndarray with the correlation matrix, has shape (num_realizations,num_realizations)

		"""

		assert self.data.dtype == np.float, "This operation is unsafe with non float numbers!!"
		if self.num_realizations==1:
			return np.ones((1,1))
		else:

			if not hasattr(self,"_mean"):
				self.mean()

			subtracted = self.data - self._mean[np.newaxis,:]
			std = np.sqrt((subtracted**2).sum(-1))

			return np.dot(subtracted,subtracted.T) / np.outer(std,std)

	def bootstrap(self,callback,bootstrap_size=10,resample=10,seed=None):

		"""
		Computes a custom statistic on the Ensemble using the bootstrap method

		:param callback: statistic to compute on the ensemble; takes the resampled Ensemble data as an input
		:type callback: callable

		:param bootstrap_size: size of the resampled ensembles used in the bootstraping; must be less than or equal to the number of realizations in the Ensemble
		:type bootstrap_size: int.

		:param resample: number of times the Ensemble is resampled
		:type resample: int.

		:param seed: if not None, this is the random seed of the random resamples
		:type seed: int.

		:returns: the bootstraped Ensemble statistic

		"""

		#Safety check
		assert bootstrap_size<=self.num_realizations,"The size of the resampling cannot exceed the original number of realizations"

		#Set the random seed
		if seed is not None:
			np.random.seed(seed)

		#TODO: Parallelize
		M = map

		#Construct the randomization matrix
		randomizer = np.random.randint(self.num_realizations,size=(resample,bootstrap_size))

		#Compute the statistic with the callback
		statistic = np.array(M(callback,self[randomizer]))

		#Return the bootstraped statistic expectation value
		return statistic.mean(0)

	def principalComponents(self):

		"""
		Computes the principal components of the Ensemble

		:returns: pcaHandler instance

		"""

		pca = PCA()
		pca.fit(self.data)
		return pca

	def compare(self,rhs,**kwargs):

		"""
		Computes a comparison score between two Ensembles: computes a chi2-style difference between two different ensembles to assert how different they are

		"""

		assert isinstance(rhs,self.__class__)
		assert self.metric == rhs.metric,"The two ensemble instances must have the same metric!!"

		if self.metric=="chi2":

			mean1 = self.mean()

			if "covariance" in kwargs.keys():
				covariance = kwargs["covariance"]
			else:
				covariance = self.covariance()

			mean2 = rhs.mean()

			return np.dot(mean1 - mean2,np.dot(np.linalg.inv(covariance),mean1 - mean2))

		else:

			raise ValueError("Only chi2 metric implemented so far!!")

	def selfChi2(self):

		"""
		Computes the Ensemble distribution of chi squared values, defined with respect to the same Ensemble mean and covariance

		:returns: array with the self chi squared for each realization

		"""

		#Compute mean and covariance
		mean = self.mean()
		covariance = self.covariance()
		difference = self.data - mean[None]

		#Compute the chi2 for each realization
		return (difference * np.linalg.solve(covariance,difference.T).T).sum(-1)

	def shuffle(self,seed=None):

		"""
		Changes the order of the realizations in the Ensemble

		:param seed: random seed for the random shuffling
		:type seed: int.

		"""

		if seed is not None:
			np.random.seed(seed)

		#Shuffle
		np.random.shuffle(self.data)

	def split(self,index):

		"""
		Inverse of the * operator: this method uses an Indexer instance to break down a multiple descriptor ensemble in many, smaller, single descriptor ensembles

		:param index: index of descriptors with which to perform the split
		:type index: Indexer instance

		:returns: list of Ensemble instances, one for each element in index

		:raises: AssertionError if shape of the ensemble data is not suitable

		"""

		assert isinstance(index,Indexer)

		splitted = list()

		for n in range(index.num_descriptors):

			splitted.append(self.__class__(file_list=self.file_list,num_realizations=self.num_realizations,data=self.data[:,index[n].first:index[n].last],metric=self.metric))

		return splitted

	#####################################################################################################################

	def __add__(self,rhs):

		"""
		Overload of the sum operator: a sum between ensembles is another ensemble whose data are vstacked. This operation is useful if you have two ensembles with different, independent realizations and you want to create an ensemble which is the union of the two. The data must be vstackable

		"""

		#Safety checks
		assert isinstance(rhs,self.__class__)
		assert self.metric == rhs.metric,"The two ensemble instances must have the same metric!!"

		if self.data is not None:
			new_data = np.vstack((self.data,rhs.data))
		else:
			new_data = rhs.data

		return self.__class__(file_list=self.file_list+rhs.file_list,data=new_data,num_realizations=self.num_realizations+rhs.num_realizations,metric=self.metric)

	def __mul__(self,rhs):

		"""
		Overload of the multiplication operator: a multiplication of two ensembles is another ensemble whose data are hstacked. This operation is useful if you have two ensembles with the same realizations but different statistical descriptor, and you want to create an ensemble in which both descriptors are included. The data must be hstackable

		"""

		#Safety checks
		assert isinstance(rhs,self.__class__)
		assert self.metric == rhs.metric
		assert self.num_realizations == rhs.num_realizations

		new_data = np.hstack((self.data,rhs.data))

		return self.__class__(file_list=list(set(self.file_list + rhs.file_list)),data=new_data,num_realizations=self.num_realizations,metric=self.metric)

	def __getitem__(self,n):

		"""
		Retrieves the n-th realization of the ensemble (starting from 0)

		:raises: IndexError if out of bounds

		"""

		return self.data[n]

	#Protect Ensemble against changes in the data
	def __setattr__(self,a,x):

		"""
		This overload protects the Ensemble against inplace changes in its data

		"""

		#Call parent method
		super(Ensemble,self).__setattr__(a,x)

		#Recompute means if ensemble data is modified
		if a=="data" and hasattr(self,"_mean"):
			self.mean()

max_stars_repo_name: asabyr/LensTools

max_issues_repo_name: asabyr/LensTools

max_forks_repo_name: asabyr/LensTools

reason: import numpy,import scipy,from scipy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import scipy.stats
import scipy.cluster.hierarchy

def cluster_gnn_map(X, by="population"):
    """Cluster GNN values.

    Args:
        X (:class:`~pd.DataFrame`) :
            A DataFrame consisting of gnn proportions. Each row
            corresponds to a chromosome from an individual in a
            population

    """
    # Zscore normalise
    for col in list(X):
        X[col] = scipy.stats.zscore(X[col])

    row_linkage = scipy.cluster.hierarchy.linkage(
        X, method="average", optimal_ordering=True
    )
    order = scipy.cluster.hierarchy.leaves_list(row_linkage)
    x_pop = X.index.values[order]
    X = X.reindex(x_pop)

    return X, row_linkage

max_stars_repo_name: percyfal/wg-genealogy-smk

max_issues_repo_name: percyfal/wg-genealogy-smk

max_forks_repo_name: percyfal/wg-genealogy-smk

reason: import scipy


text:
# -*- coding: utf-8 -*-
# pc_type           lenovo
# create_time:      2019/11/10 9:13
# file_name:        4_1.py

import numpy as np
import cv2
import matplotlib.pyplot as plt
import seaborn as sns
# 设置中文字体和负号正常显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
sns.set_context("paper") # 背景
sns.set_style('whitegrid') # 主题
sns.set(font='SimHei')  # 解决Seaborn中文显示问题，这一句必须放在前两后面

img1 = cv2.imread("img/test7.tif",flags=0)
img2 = cv2.imread("img/test8.tif",flags=0)

def my_grad(img):
    height,width = img.shape[0],img.shape[1]
    Grad = np.zeros(img.shape,dtype=np.int8)
    # 不对右、下边界作处理、影响不大
    for h in range(height-1):
        for w in range(width-1):
            gx = int(img[h,w+1])- int(img[h,w])
            gy = int(img[h+1,w]) -int(img[h,w])
            Grad[h][w] = abs(gx) + abs(gy)
    return Grad

def my_roberts(img):
    # Roberts算子
    kernelx = np.array([[-1, 0], [0, 1]], dtype=int)
    kernely = np.array([[0, -1], [1, 0]], dtype=int)
    x = cv2.filter2D(img, cv2.CV_16S, kernelx)
    y = cv2.filter2D(img, cv2.CV_16S, kernely)
    # 转uint8
    absX = cv2.convertScaleAbs(x)
    absY = cv2.convertScaleAbs(y)
    Roberts = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)
    return Roberts

def my_prewitt(img):
    # Prewitt算子
    kernelx = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]], dtype=int)
    kernely = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], dtype=int)
    x = cv2.filter2D(img, cv2.CV_16S, kernelx)
    y = cv2.filter2D(img, cv2.CV_16S, kernely)
    # 转uint8
    absX = cv2.convertScaleAbs(x)
    absY = cv2.convertScaleAbs(y)
    Prewitt = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)
    return Prewitt

def my_sobel(img):
    # Sobel算子
    x = cv2.Sobel(img, cv2.CV_16S, 1, 0)  # 对x求一阶导
    y = cv2.Sobel(img, cv2.CV_16S, 0, 1)  # 对y求一阶导
    absX = cv2.convertScaleAbs(x)
    absY = cv2.convertScaleAbs(y)
    Sobel = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)
    return Sobel

plt.figure()
plt.subplot(2,5,1)
plt.imshow(img1,cmap='gray')
plt.title("原图")
plt.subplot(2,5,2)
mg1 = my_grad(img1)
plt.imshow(mg1,cmap='gray')
plt.title("梯度算子")
plt.subplot(2,5,3)
mr1 = my_roberts(img1)
plt.imshow(mr1,cmap='gray')
plt.title("Roberts算子")
plt.subplot(2,5,4)
mp1 = my_prewitt(img1)
plt.imshow(mp1,cmap='gray')
plt.title("Prewitt算子")
plt.subplot(2,5,5)
ms1 = my_sobel(img1)
plt.imshow(ms1,cmap='gray')
plt.title("Sobel算子")

plt.subplot(2,5,6)
plt.imshow(img2,cmap='gray')
plt.title("原图")
plt.subplot(2,5,7)
mg2 = my_grad(img2)
plt.imshow(mg2,cmap='gray')
plt.title("梯度算子")
plt.subplot(2,5,8)
mr2 = my_roberts(img2)
plt.imshow(mr2,cmap='gray')
plt.title("Roberts算子")
plt.subplot(2,5,9)
mp2 = my_prewitt(img2)
plt.imshow(mp2,cmap='gray')
plt.title("Prewitt算子")
plt.subplot(2,5,10)
ms2 = my_sobel(img2)
plt.imshow(ms2,cmap='gray')
plt.title("Sobel算子")

def my_laplacian(img):
    dst = cv2.Laplacian(img, cv2.CV_16S, ksize=3)
    Laplacian = cv2.convertScaleAbs(dst)
    return Laplacian

plt.figure()
plt.subplot(2,2,1)
plt.imshow(img1,cmap='gray')
plt.title("原图")
plt.subplot(2,2,2)
ml1 = my_laplacian(img1)
plt.imshow(ml1,cmap='gray')
plt.title("Laplacian算子")
plt.subplot(2,2,3)
plt.imshow(img2,cmap='gray')
plt.title("原图")
plt.subplot(2,2,4)
ml2 = my_laplacian(img2)
plt.imshow(ml2,cmap='gray')
plt.title("Laplacian算子")

plt.figure()

def my_adaptivethreshold(img):
    result_img = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY,11,2)
    return result_img

plt.subplot(2,5,1)
plt.imshow(mg1,cmap='gray')
plt.title("梯度算子")
plt.subplot(2,5,2)
plt.imshow(mr1,cmap='gray')
plt.title("Roberts算子")
plt.subplot(2,5,3)
plt.imshow(mp1,cmap='gray')
plt.title("Prewitt算子")
plt.subplot(2,5,4)
plt.imshow(ms1,cmap='gray')
plt.title("Sobel算子")
plt.subplot(2,5,5)
plt.imshow(ml1,cmap='gray')
plt.title("Laplacian算子")
plt.subplot(2,5,6)
mg1 = np.uint8(mg1)
plt.imshow(my_adaptivethreshold(mg1),cmap='gray')
plt.title("梯度算子(自适应阈值处理)")
plt.subplot(2,5,7)
plt.imshow(my_adaptivethreshold(mr1),cmap='gray')
plt.title("Roberts算子(自适应阈值处理)")
plt.subplot(2,5,8)
plt.imshow(my_adaptivethreshold(mp1),cmap='gray')
plt.title("Prewitt算子(自适应阈值处理)")
plt.subplot(2,5,9)
plt.imshow(my_adaptivethreshold(ms1),cmap='gray')
plt.title("Sobel算子(自适应阈值处理)")
plt.subplot(2,5,10)
plt.imshow(my_adaptivethreshold(ml1),cmap='gray')
plt.title("Laplacian算子(自适应阈值处理)")
plt.show()

max_stars_repo_name: inspurer/ImageProcess

max_issues_repo_name: inspurer/ImageProcess

max_forks_repo_name: inspurer/ImageProcess


text:
from abc import ABC
from functools import cached_property, lru_cache
from typing import Union, Sequence

import numpy as np

from pyqumo.cqumo.randoms import Variable, RandomsFactory
from pyqumo.chains import ContinuousTimeMarkovChain, DiscreteTimeMarkovChain
from pyqumo.matrix import cbdiag, order_of, \
    check_markovian_arrival, fix_markovian_arrival, str_array, \
    is_subinfinitesimal, is_pmf, fix_infinitesimal, fix_stochastic
from pyqumo.random import Distribution, Exponential
from pyqumo.errors import MatrixShapeError

class RandomProcess(ABC, Distribution):

    @lru_cache
    def lag(self, n: int) -> float:
        """
        Get auto-correlation coefficient with lag n of the random process.

        Notes
        -----

        :math:`r_k = (E[X_{t+n} - m_1][X_{t} - m_1]) / s^2`, where `m_1` -
        mean value and :math:`s^2` - variance (dispersion).

        Parameters
        ----------
        n : int
            Time lag (number of steps).

        Returns
        -------
        value : float

        Raises
        ------
        ValueError
            raised if n is not an integer or is non-positive
        """
        if n < 0 or (n - np.floor(n)) > 0:
            raise ValueError(f'positive integer expected, but {n} found')
        if n == 0:
            return 1
        return self._lag(n)

    def _lag(self, n: int) -> float:
        """
        Get lag-n autocorrelation. In this method it can be assumed that
        `n` is a non-zero positive integer.
        """
        raise NotImplementedError

    def copy(self) -> 'RandomProcess':
        raise NotImplementedError

class GIProcess(RandomProcess):
    """
    GI-arrival model. Samples of this kind of process are built from a
    known distribution, and they don't change over lifetime.

    Poisson process is an example of GI-process with exponential arrivals.
    """
    def __init__(self, dist: Distribution, factory: RandomsFactory = None):
        """
        Constructor.

        Parameters
        ----------
        dist : random distribution
        """
        super().__init__(factory)
        if dist is None:
            raise ValueError('distribution required')
        self._dist = dist

    @property
    def dist(self) -> Distribution:
        """
        Get random distribution.
        """
        return self._dist

    @cached_property
    def mean(self) -> float:
        return self._dist.mean

    @cached_property
    def var(self) -> float:
        return self._dist.var

    @cached_property
    def std(self) -> float:
        return self._dist.std

    @cached_property
    def cv(self) -> float:
        return self._dist.std / self._dist.mean

    @property
    def rnd(self) -> Variable:
        return self._dist.rnd

    def _moment(self, n: int) -> float:
        return self._dist.moment(n)

    def _lag(self, k: int) -> float:
        return 0.0  # always 0.0

    def _eval(self, size: int) -> np.ndarray:
        return self._dist(size)

    def copy(self) -> 'GIProcess':
        return GIProcess(self._dist.copy())

    def __repr__(self):
        return f'(GI: f={self.dist})'

class Poisson(GIProcess):
    """
    Custom case of GI-process with exponential arrivals.
    """
    def __init__(self, rate: float, factory: RandomsFactory = None):
        """
        Constructor.

        Parameters
        ----------
        rate : float
            Rate (1 / mean) of exponential distribution.
        """
        if rate <= 0.0:
            raise ValueError(f"positive rate expected, {rate} found")
        super().__init__(Exponential(rate), factory)

    def __repr__(self):
        return f'(Poisson: r={self.rate:.3g})'

class MarkovArrival(RandomProcess):
    """
    Markovian arrival process (MAP).
    """
    def __init__(
            self,
            d0: Union[np.ndarray, Sequence[Sequence[float]]],
            d1: Union[np.ndarray, Sequence[Sequence[float]]],
            safe: bool = False,
            tol: float = 1e-3,
            factory: RandomsFactory = None):
        """
        Create MAP process with the given D0 and D1 matrices.

        If `safe = False` is set (this is default), we validate matrices
        D0 and D1. These matrices must comply three requirements:

        1) Sum `D0 + D1` must be an infinitesimal matrix
        2) All elements of D1 must be non-negative
        3) All elements of D0 except the main diagonal must be non-negative

        To avoid problems with "1e-12 is not zero", constructor accepts
        parameters `tol` (tolerance). If matrices D0 and D1 fail the check,
        it tries to fix them using `fix_markovian_arrival()` call with this
        tolerance.

        Sometimes it is useful to avoid matrices validation (e.g., when MAP
        matrices are obtained as a result of another computation). To disable
        validation, set `safe = True`.

        Parameters
        ----------
        d0 : array_like
            Matrix D0, it must be subinfinitesimal
        d1 : array_like
            Matrix D1, all its elements must be non-negative
        safe : bool, optional
            Flag indicating whether not to validate or try to fix matrices D0
            and D1 if they don't satisfy requirements. By default, `False`.
        tol : float, optional
            If `safe = False` and matrices D0 and D1 violate requirements,
            try to fix errors those are less then `tol`. If `safe = True`
            this parameter is ignored. Default: 1e-3.
        """
        super().__init__(factory)
        need_copy = [False, False]
        matrices = [d0, d1]

        for i in range(2):
            if not isinstance(matrices[i], np.ndarray):
                matrices[i] = np.asarray(matrices[i])
            else:
                need_copy[i] = True

        if not safe and not check_markovian_arrival(matrices):
            matrices, _ = fix_markovian_arrival(matrices, tol=tol)
            # Since matrices are re-built, no need to copy them:
            for i in range(2):
                need_copy[i] = False

        # Copy matrices if needed:
        for i in range(2):
            if need_copy[i]:
                matrices[i] = matrices[i].copy()

        # Extract matrices:
        self._matrices = tuple(matrices)
        self._generator = sum(self._matrices)
        self._order = order_of(self._matrices[0])
        self._inv_d0 = np.linalg.inv(self._matrices[0])

        # Since we will use Mmap for data generation, we need to generate
        # special matrices for transitions probabilities and rates.

        # Store rates:
        self._rates = -self._matrices[0].diagonal()

        # Build embedded DTMC and CTMC:
        self._ctmc = ContinuousTimeMarkovChain(
            self._generator,
            safe=True
        )
        self._dtmc = DiscreteTimeMarkovChain(
            -self._inv_d0.dot(self.d1),
            safe=True
        )
        self._states = [
            Exponential(rate, factory=factory)
            for rate in self._rates
        ]

    @staticmethod
    def erlang(shape: int, rate: float) -> 'MarkovArrival':
        """
        MAP representation of Erlang process with the given shape and rate.

        Parameters
        ----------
        shape : int
            Number of phases
        :param : float
            Rate at each phase
        """
        d0 = cbdiag(shape, [
            (0, np.asarray([[-rate]])),
            (1, np.asarray([[rate]]))
        ])
        d1 = np.zeros((shape, shape))
        d1[shape-1, 0] = rate
        return MarkovArrival(d0, d1)

    @staticmethod
    def poisson(rate: float) -> 'MarkovArrival':
        """
        MAP representation of a Poisson process with the given rate.

        Parameters
        ----------
        rate : float
            Exponential distribution rate
        """
        return MarkovArrival([[-rate]], [[rate]])

    @staticmethod
    def phase_type(
            s: Sequence[Sequence[float]],
            p: Sequence[float]) -> 'MarkovArrival':
        """
        MAP representation of a PH distribution.

        Parameters
        ----------
        s : square matrix
            PH distribution generator (order NxN)
        p : vector
            PH distribution initial probabilities (order N)
        """
        if not isinstance(s, np.ndarray):
            s = np.asarray(s)
        if not isinstance(p, np.ndarray):
            p = np.asarray(p)

        # Validate:
        order = order_of(s)
        if order != order_of(p):
            raise MatrixShapeError(f'({order},)', p.shape, 'PMF')
        if not is_subinfinitesimal(s):
            s = fix_infinitesimal(s, sub=True)
        if not is_pmf(p):
            p = fix_stochastic(p)

        # Build:
        d0 = s
        d1 = (-s.dot(np.ones((order, 1)))).dot(p.reshape((1, order)))
        return MarkovArrival(d0, d1)

    def copy(self) -> 'MarkovArrival':
        """
        Build a new MAP with the same matrices D0 and D1 without validation.
        """
        return MarkovArrival(self.d0, self.d1, safe=True)

    @property
    def d0(self) -> np.ndarray:
        """
        Get matrix D0.
        """
        return self._matrices[0]

    @property
    def d1(self) -> np.ndarray:
        """
        Get matrix D1.
        """
        return self._matrices[1]

    def d(self, n: int) -> np.ndarray:
        """
        Get matrix Dn for n = 0 or n = 1.

        Parameters
        ----------
        n : int (0 or 1)
        """
        return self._matrices[n]

    @cached_property
    def generator(self) -> np.ndarray:
        return self._generator

    @cached_property
    def order(self) -> int:
        return self._order

    def scale(self, k) -> 'MarkovArrival':
        """
        Get a new MAP with mean = (this MAP) x K.
        """
        return MarkovArrival(self.d0 / k, self.d1 / k)

    @lru_cache
    def d0n(self, k: int) -> np.ndarray:
        """
        Returns :math:`(-D0)^{k}`.
        """
        if k == -1:
            return -self._inv_d0
        if k == 0:
            return np.eye(self.order)
        if k > 0:
            return self.d0n(k - 1).dot(-self.d0)
        # If we are here, degree <= -2
        return self.d0n(k + 1).dot(self.d0n(-1))

    @lru_cache
    def _moment(self, n: int) -> float:
        pi = self.dtmc.steady_pmf
        x = np.math.factorial(n)*pi.dot(self.d0n(-n)).dot(np.ones(self.order))
        return x.item()

    @lru_cache
    def _lag(self, n: int) -> float:
        #
        # Computing lag-k as:
        #
        #   r^2 * pi * (-D0)^(-1) * P^k * (-D0)^(-1) * 1s - 1
        #   -------------------------------------------------- ,
        #   2 * r^2 * pi * (-D0)^(-2) * 1s - 1
        #
        # where r - rate (\lambda), pi - stationary distribution of the
        # embedded DTMC, 1s - vector of ones of MAP order
        #
        dtmc_matrix_k = self._pow_dtmc_matrix(n)
        pi = self.dtmc.steady_pmf
        rate2 = pow(self.rate, 2.0)
        e = np.ones(self.order)
        d0ni = self.d0n(-1)
        d0ni2 = self.d0n(-2)

        num = rate2 * pi.dot(d0ni).dot(dtmc_matrix_k).dot(d0ni).dot(e) - 1
        den = (2 * rate2 * pi.dot(d0ni2).dot(e) - 1)
        return num / den

    @property
    def ctmc(self) -> ContinuousTimeMarkovChain:
        """
        Get background continuous-time Markov chain (CTMC).
        """
        return self._ctmc

    @property
    def dtmc(self) -> DiscreteTimeMarkovChain:
        """
        Get embedded (at arrivals) discrete-time Markov chain (DTMC).
        """
        return self._dtmc

    @cached_property
    def rnd(self) -> Variable:
        # We need to store cumulative transition probabilities P.
        #    We store them in a stochastic matrix of shape N x (K*N):
        #      P[I, J] is a probability, that:
        #      - new state will be J (mod N), and
        #      - if J < N, then no packet is generated, or
        #      - if J >= N, then packet of type J // N is generated.
        vars = [state.rnd for state in self._states]
        trans_pmf = np.hstack((
            self._matrices[0] + np.diag(self._rates),
            self._matrices[1]
        )) / self._rates[:, None]
        return self.factory.createSemiMarkovArrivalVariable(
            vars, self.dtmc.steady_pmf, trans_pmf)

    def compose(self, other: 'MarkovArrival') -> 'MarkovArrival':
        # TODO:  write unit tests
        if not isinstance(other, MarkovArrival):
            raise TypeError(f'expected MarkovArrivalProcess, '
                            f'{type(other)} found')
        self_eye = np.eye(self.order)
        other_eye = np.eye(other.order)
        d0_out = np.kron(self.d0, other_eye) + np.kron(other.d0, self_eye)
        d1_out = np.kron(self.d1, other_eye) + np.kron(other.d1, self_eye)
        return MarkovArrival(d0_out, d1_out)

    @lru_cache
    def _pow_dtmc_matrix(self, k):
        if k == 0:
            return np.eye(self.order)
        elif k > 0:
            return self._pow_dtmc_matrix(k - 1).dot(self.dtmc.matrix)
        raise ValueError(f"k={k} must be non-negative")

    def __repr__(self):
        return f'(MAP: d0={str_array(self.d0)}, d1={str_array(self.d1)})'

max_stars_repo_name: larioandr/thesis-queues

max_issues_repo_name: larioandr/thesis-queues

max_forks_repo_name: larioandr/thesis-queues


text:
#-*- coding: utf-8 -*-

import os
import torch
import argparse
import numpy as np
from tqdm import tqdm
from transformers import  BertTokenizer
from dataset import DualSample, TokenizedSample, OriginalDataset

def tokenize_data(data, mode='train'):
    max_forward_asp_query_length = 0
    max_forward_opi_query_length = 0
    max_sentiment_query_length = 0

    max_aspect_num = 0
    tokenized_sample_list = []

    header_fmt = 'Tokenize data {:>5s}'
    for sample in tqdm(data, desc=f"{header_fmt.format(mode.upper())}"):
        forward_queries = []
        forward_answers = []
        sentiment_queries = []
        sentiment_answers = []

        forward_queries_seg = []
        sentiment_queries_seg = []

        if int(len(sample.sentiment_queries)) > max_aspect_num:
            max_aspect_num = int(len(sample.sentiment_queries))

        for idx in range(len(sample.forward_queries)):
            temp_query = sample.forward_queries[idx]
            temp_text = sample.text
            temp_answer = sample.forward_answers[idx]
            temp_query_to = ['[CLS]'] + temp_query + ['[SEP]'] + temp_text
            temp_query_seg = [0] * (len(temp_query) + 2) + [1] * len(temp_text)
            temp_answer[0] = [-1] * (len(temp_query) + 2) + temp_answer[0]
            temp_answer[1] = [-1] * (len(temp_query) + 2) + temp_answer[1]

            assert len(temp_answer[0]) == len(temp_answer[1]) == len(temp_query_to) == len(temp_query_seg)

            if len(temp_query_to) > max_forward_asp_query_length:
                max_forward_asp_query_length = len(temp_query_to)

            forward_queries.append(temp_query_to)
            forward_answers.append(temp_answer)
            forward_queries_seg.append(temp_query_seg)

        for idx in range(len(sample.sentiment_queries)):
            temp_query = sample.sentiment_queries[idx]
            temp_text = sample.text
            temp_answer = sample.sentiment_answers[idx]
            temp_query_to = ['[CLS]'] + temp_query + ['[SEP]'] + temp_text
            temp_query_seg = [0] * (len(temp_query) + 2) + [1] * len(temp_text)

            assert len(temp_query_to) == len(temp_query_seg)

            if len(temp_query_to) > max_sentiment_query_length:
                max_sentiment_query_length = len(temp_query_to)

            sentiment_queries.append(temp_query_to)
            sentiment_answers.append(temp_answer)
            sentiment_queries_seg.append(temp_query_seg)

        # import numpy as np
        # print(f"forward_queries: {np.shape(forward_queries)} {type(forward_queries)} | {forward_queries}")
        # print(f"forward_answers: {np.shape(forward_answers)}")
        # print(f"sentiment_queries: {np.shape(sentiment_queries)} | {sentiment_queries}")
        # print(f"sentiment_answers: {np.shape(sentiment_answers)} | {sentiment_answers}")
        # print(f"forward_queries_seg: {np.shape(forward_queries_seg)} | {forward_queries_seg}")
        # print(f"sentiment_queries_seg: {np.shape(sentiment_queries_seg)}")

        temp_sample = TokenizedSample(
            sample.original_sample, forward_queries,
            forward_answers, sentiment_queries,
            sentiment_answers, forward_queries_seg,
            sentiment_queries_seg
        )
        # print(temp_sample)
        tokenized_sample_list.append(temp_sample)

    max_attributes = {
        'mfor_asp_len': max_forward_asp_query_length,
        'max_sent_len': max_sentiment_query_length,
        'max_aspect_num': max_aspect_num
    }
    return tokenized_sample_list, max_attributes

def preprocessing(sample_list, max_len, mode='train'):
    _tokenizer = BertTokenizer.from_pretrained('bert-base-uncased')
    _forward_asp_query = []
    _forward_asp_answer_start = []
    _forward_asp_answer_end = []
    _forward_asp_query_mask = []
    _forward_asp_query_seg = []

    _sentiment_query = []
    _sentiment_answer = []
    _sentiment_query_mask = []
    _sentiment_query_seg = []

    _aspect_num = []

    header_fmt = 'Preprocessing {:>5s}'
    for instance in tqdm(sample_list, desc=f"{header_fmt.format(mode.upper())}"):
        f_query_list = instance.forward_queries
        f_answer_list = instance.forward_answers
        f_query_seg_list = instance.forward_seg

        s_query_list = instance.sentiment_queries
        s_answer_list = instance.sentiment_answers
        s_query_seg_list = instance.sentiment_seg

        # _aspect_num: 1/2/3/...
        _aspect_num.append(int(len(s_query_list)))

        # Forward
        # Aspect
        # query
        assert len(f_query_list[0]) == len(f_answer_list[0][0]) == len(f_answer_list[0][1])
        f_asp_pad_num = max_len['mfor_asp_len'] - len(f_query_list[0])

        _forward_asp_query.append(_tokenizer.convert_tokens_to_ids(
            [word.lower() if word not in ['[CLS]', '[SEP]'] else word for word in f_query_list[0]]))
        _forward_asp_query[-1].extend([0] * f_asp_pad_num)

        # query_mask
        _forward_asp_query_mask.append([1 for i in range(len(f_query_list[0]))])
        _forward_asp_query_mask[-1].extend([0] * f_asp_pad_num)

        # answer
        _forward_asp_answer_start.append(f_answer_list[0][0])
        _forward_asp_answer_start[-1].extend([-1] * f_asp_pad_num)
        _forward_asp_answer_end.append(f_answer_list[0][1])
        _forward_asp_answer_end[-1].extend([-1] * f_asp_pad_num)

        # seg
        _forward_asp_query_seg.append(f_query_seg_list[0])
        _forward_asp_query_seg[-1].extend([1] * f_asp_pad_num)

        # Sentiment
        single_sentiment_query = []
        single_sentiment_query_mask = []
        single_sentiment_query_seg = []
        single_sentiment_answer = []
        for j in range(len(s_query_list)):
            sent_pad_num = max_len['max_sent_len'] - len(s_query_list[j])
            single_sentiment_query.append(_tokenizer.convert_tokens_to_ids(
                [word.lower() if word not in ['[CLS]', '[SEP]'] else word for word in s_query_list[j]]))
            single_sentiment_query[-1].extend([0] * sent_pad_num)

            single_sentiment_query_mask.append([1 for i in range(len(s_query_list[j]))])
            single_sentiment_query_mask[-1].extend([0] * sent_pad_num)

            # query_seg
            single_sentiment_query_seg.append(s_query_seg_list[j])
            single_sentiment_query_seg[-1].extend([1] * sent_pad_num)

            single_sentiment_answer.append(s_answer_list[j])

        _sentiment_query.append(single_sentiment_query)
        _sentiment_query[-1].extend(
            [[0 for i in range(max_len['max_sent_len'])]] * (max_len['max_aspect_num'] - _aspect_num[-1]))

        _sentiment_query_mask.append(single_sentiment_query_mask)
        _sentiment_query_mask[-1].extend(
            [[0 for i in range(max_len['max_sent_len'])]] * (max_len['max_aspect_num'] - _aspect_num[-1]))

        _sentiment_query_seg.append(single_sentiment_query_seg)
        _sentiment_query_seg[-1].extend(
            [[0 for i in range(max_len['max_sent_len'])]] * (max_len['max_aspect_num'] - _aspect_num[-1]))

        _sentiment_answer.append(single_sentiment_answer)
        _sentiment_answer[-1].extend([-1] * (max_len['max_aspect_num'] - _aspect_num[-1]))

    result = {
        "_forward_asp_query": _forward_asp_query,
        "_forward_asp_answer_start": _forward_asp_answer_start,
        "_forward_asp_answer_end": _forward_asp_answer_end,
        "_forward_asp_query_mask": _forward_asp_query_mask,
        "_forward_asp_query_seg": _forward_asp_query_seg,
        "_sentiment_query": _sentiment_query,
        "_sentiment_answer": _sentiment_answer,
        "_sentiment_query_mask": _sentiment_query_mask,
        "_sentiment_query_seg": _sentiment_query_seg,
        "_aspect_num": _aspect_num,
    }

    return OriginalDataset(result)

if __name__ == '__main__':
    parser = argparse.ArgumentParser()

    parser.add_argument('--data_path', type=str, default='./data/14rest/preprocess',
                        help='Path to the processed data from `data_process.py`')
    parser.add_argument('--output_path', type=str, default='./data/14rest/preprocess',
                        help='Path to the saved data.')

    args = parser.parse_args()

    train_data_path = f"{args.data_path}/train_DUAL.pt"
    dev_data_path = f"{args.data_path}/dev_DUAL.pt"
    test_data_path = f"{args.data_path}/test_DUAL.pt"

    train_data = torch.load(train_data_path)
    dev_data = torch.load(dev_data_path)
    test_data = torch.load(test_data_path)

    train_tokenized, train_max_len = tokenize_data(train_data, mode='train')
    dev_tokenized, dev_max_len = tokenize_data(dev_data, mode='dev')
    test_tokenized, test_max_len = tokenize_data(test_data, mode='test')

    print(f"\nMax attributes")
    print(f"train_max_len : {train_max_len}")
    print(f"dev_max_len : {dev_max_len}")
    print(f"test_max_len : {test_max_len}\n")

    train_preprocess = preprocessing(train_tokenized, train_max_len, mode='train')
    dev_preprocess = preprocessing(dev_tokenized, dev_max_len, mode='dev')
    test_preprocess = preprocessing(test_tokenized, test_max_len, mode='test')

    if not os.path.exists(args.output_path):
        os.makedirs(args.output_path)

    output_path = f"{args.output_path}/data.pt"
    print(f"Saved data : `{output_path}`.")
    torch.save({
        'train': train_preprocess,
        'dev': dev_preprocess,
        'test': test_preprocess
    }, output_path)

max_stars_repo_name: phanxuanphucnd/ABSA-Pair-Extraction

max_issues_repo_name: phanxuanphucnd/ABSA-Pair-Extraction

max_forks_repo_name: phanxuanphucnd/ABSA-Pair-Extraction


text:
import unittest
from unittest import TestCase
import giant.ray_tracer.rays as g_rays
import giant.rotations as at
import numpy as np
import copy

class TestRay(TestCase):

    def setUp(self):

        self.ray_start = np.array([1, 2, 3])
        self.ray_direction = np.array([4, 5, 6])
        self.ray = g_rays.Rays(self.ray_start, self.ray_direction)

        self.rays_start = np.array([[1, 3], [4, 5], [6, 7]])
        self.rays_direction = np.array([[4, 5], [5, 6], [7, 8]])
        self.rays = g_rays.Rays(self.rays_start, self.rays_direction)
        self.num_rays = 2

        self.rotation = at.Rotation([np.pi, np.pi / 2, np.pi / 4])

    def test_creation_one_ray(self):

        for seq_type in [list, np.array, tuple]:

            with self.subTest(seq_type=seq_type):

                ray = g_rays.Rays(seq_type(self.ray_start), seq_type(self.ray_direction))

                np.testing.assert_array_equal(ray._start, self.ray_start)
                np.testing.assert_array_equal(ray._direction, self.ray_direction)

                self.assertEqual(ray.num_rays, 1)

    def test_creation_multi_ray(self):

        for seq_type in [list, np.array, tuple]:

            with self.subTest(seq_type=seq_type):

                rays = g_rays.Rays(seq_type(self.rays_start), seq_type(self.rays_direction))

                np.testing.assert_array_equal(rays._start, self.rays_start)
                np.testing.assert_array_equal(rays._direction, self.rays_direction)

                self.assertEqual(rays.num_rays, self.num_rays)

    def test_single_iter(self):

        for ray in self.ray:
            self.assertIsInstance(ray, g_rays.Rays)

            np.testing.assert_array_equal(ray._start, self.ray_start)
            np.testing.assert_array_equal(ray._direction, self.ray_direction)

    def test_multi_iter(self):

        for ind, ray in enumerate(self.rays):
            self.assertIsInstance(ray, g_rays.Rays)

            np.testing.assert_array_equal(ray._start, self.rays_start[:, ind])
            np.testing.assert_array_equal(ray._direction, self.rays_direction[:, ind])

    def test_single_getitem(self):

        with self.assertRaises(ValueError):
            _ = self.ray[0]

    def test_multi_getitem(self):

        second = self.rays[1]

        np.testing.assert_array_equal(second._start, self.rays_start[:, 1])
        np.testing.assert_array_equal(second._direction, self.rays_direction[:, 1])

    def test_single_len(self):

        self.assertEqual(len(self.ray), 1)

    def test_multi_len(self):

        self.assertEqual(len(self.rays), 2)

    # TODO: Multiplication tests -- maybe

    def test_single_rotate(self):

        ray_copy = copy.deepcopy(self.ray)

        ray_copy.rotate(self.rotation)

        np.testing.assert_array_almost_equal(ray_copy._start, np.matmul(self.rotation.matrix, self.ray._start),
                                             err_msg="Rotation rotation ray start")
        np.testing.assert_array_almost_equal(ray_copy._direction,
                                             np.matmul(self.rotation.matrix, self.ray._direction),
                                             err_msg="Rotation rotation ray direction")

        ray_copy = copy.deepcopy(self.ray)

        ray_copy.rotate(self.rotation.matrix)

        np.testing.assert_array_almost_equal(ray_copy._start, np.matmul(self.rotation.matrix, self.ray._start),
                                             err_msg="Matrix rotation ray start")
        np.testing.assert_array_almost_equal(ray_copy._direction,
                                             np.matmul(self.rotation.matrix, self.ray._direction),
                                             err_msg="Matrix rotation ray direction")

    def test_multi_rotate(self):

        ray_copy = copy.deepcopy(self.rays)

        ray_copy.rotate(self.rotation)

        np.testing.assert_array_almost_equal(ray_copy._start, np.matmul(self.rotation.matrix, self.rays._start),
                                             err_msg="Rotation rotation ray start")
        np.testing.assert_array_almost_equal(ray_copy._direction,
                                             np.matmul(self.rotation.matrix, self.rays._direction),
                                             err_msg="Rotation rotation ray direction")

        ray_copy = copy.deepcopy(self.rays)

        ray_copy.rotate(self.rotation.matrix)

        np.testing.assert_array_almost_equal(ray_copy._start, np.matmul(self.rotation.matrix, self.rays._start),
                                             err_msg="Matrix rotation ray start")
        np.testing.assert_array_almost_equal(ray_copy._direction,
                                             np.matmul(self.rotation.matrix, self.rays._direction),
                                             err_msg="Matrix rotation ray direction")

    def test_single_translate(self):

        ray_copy = copy.deepcopy(self.ray)

        ray_copy.translate(5)

        np.testing.assert_array_equal(ray_copy._start, self.ray_start + 5,
                                      err_msg="Scalar Translation")

        ray_copy = copy.deepcopy(self.ray)

        ray_copy.translate(self.ray_start)

        np.testing.assert_array_equal(ray_copy._start, self.ray_start + self.ray_start,
                                      err_msg="Vector Translation")

    def test_multi_translate(self):

        # ray_copy = copy.deepcopy(self.rays)
        #
        # ray_copy.translate(5)
        #
        # np.testing.assert_array_equal(ray_copy._start, self.rays_start + 5,
        #                               err_msg="Scalar Translation")

        ray_copy = copy.deepcopy(self.rays)

        ray_copy.translate(self.ray_start)

        np.testing.assert_array_equal(ray_copy._start, self.rays_start + self.ray_start[..., np.newaxis],
                                      err_msg="Vector Translation")

        ray_copy = copy.deepcopy(self.rays)

        ray_copy.translate(self.rays_start)

        np.testing.assert_array_equal(ray_copy._start, self.rays_start + self.rays_start,
                                      err_msg="Matrix Translation")

max_stars_repo_name: nasa/giant

max_issues_repo_name: nasa/giant

max_forks_repo_name: nasa/giant


text:
import numpy as np
import cv2 as cv

class Camera(object):

    def __init__(self):

        # Implement this constructor that opens a webcam and stores it in self._camera
        self._camera = cv.VideoCapture(0)
        if not self._camera.isOpened():
            print("Not able to open camera")
            exit()

    def capture(self):

        # Implement this function that grabs an image from the webcam and returns a numpy array
        ret, img = self._camera.read()
        if not ret:
            print("Error reading from camera")

        return img

    def __del__(self):

        # Implement this destructor. Remember to free the camera.
        self._camera.release()
        cv.destroyAllWindows()

if __name__ == "__main__":
    cam = Camera()
    while True:
        img = cam.capture()
        print(len(img.tostring()), img.dtype, img.shape)
        cv.imshow("image", img)
        cv.waitKey(1)


text:
# Copyright (c) 2020 Branislav Holländer. All rights reserved.
# See the file LICENSE for copying permission.

import pytest

import jax
import jax.tree_util
import jax.numpy as jnp

import piper.functional as func
import piper.distributions as dist
from piper import test_util as tu

def test_kl_normal_normal_one_dimensional():
    n1 = dist.normal(jnp.array(0.), jnp.array(1.))
    n2 = dist.normal(jnp.array(0.), jnp.array(1.))
    n3 = dist.normal(jnp.array(1.), jnp.array(1.))
    n4 = dist.normal(jnp.array(1.), jnp.array(2.))

    assert func.compute_kl_div(n1, n2) < 1e-6
    assert func.compute_kl_div(n1, n3) - jnp.array(0.5) \
        < jnp.array(1e-6)
    assert func.compute_kl_div(n3, n4) \
        - jnp.array(0.09657359028) < jnp.array(1e-6)

def test_kl_normal_normal_multi_dimensional():
    n1 = dist.normal(jnp.array([0., 0.]), jnp.array([1., 1.]))
    n2 = dist.normal(jnp.array([0., 0.]), jnp.array([1., 1.]))
    n3 = dist.normal(jnp.array([1., 1.]), jnp.array([1., 1.]))
    n4 = dist.normal(jnp.array([1., 1.]), jnp.array([2., 2.]))

    assert jnp.all(func.compute_kl_div(n1, n2) < 1e-6)
    assert jnp.all(func.compute_kl_div(n1, n3) - 0.5 < 1e-6)
    assert jnp.all(func.compute_kl_div(n3, n4) - 0.09657359 < 1e-6)

def test_sample_normal():
    def model1(key):
        n = func.sample('n', dist.normal(jnp.array(0.), jnp.array(1.)), key)
        return n

    keys = jax.random.split(jax.random.PRNGKey(123), 100)
    samples = jax.vmap(lambda k: model1(k))(keys)

    assert abs(jnp.mean(samples)) < 0.2

    def model2(key):
        n = func.sample('n', dist.normal(jnp.array(10.), jnp.array(1.)), key)
        return n

    keys = jax.random.split(jax.random.PRNGKey(123), 100)
    samples = jax.vmap(lambda k: model2(k))(keys)

    assert abs(jnp.mean(samples)) - 10. < 0.2

def test_sample_joint_normal():
    def model(key):
        keys = jax.random.split(key)
        weight = func.sample('weight',
                             dist.normal(jnp.array(0.), jnp.array(1.)),
                             keys[0])
        measurement = func.sample('measurement',
                                  dist.normal(weight, jnp.array(1.)),
                                  keys[1])
        return measurement

    keys = jax.random.split(jax.random.PRNGKey(123), 100)
    samples = jax.vmap(lambda k: model(k))(keys)
    assert abs(jnp.mean(samples)) < 0.3

def test_incompatible_dimensions():
    def model(key):
        keys = jax.random.split(key)
        weight = func.sample('weight',
                             dist.normal(
                                 jnp.array([0., 0.]),
                                 jnp.array([1., 1.])),
                             keys[0])
        measurement = func.sample('measurement',
                                  dist.normal(weight, jnp.array([1.])), keys[1])
        return measurement

    with pytest.raises(ValueError):
        key = jax.random.PRNGKey(123)
        model(key)

def test_sample_conditioned():
    def model(key):
        keys = jax.random.split(key)
        weight = func.sample('weight',
                             dist.normal(jnp.array(0.), jnp.array(1.)),
                             keys[0])
        measurement = func.sample('measurement',
                                  dist.normal(weight, jnp.array(1.)),
                                  keys[1])
        return measurement

    conditioned_model = func.condition(model, {'weight': jnp.array(0.)})
    keys = jax.random.split(jax.random.PRNGKey(123), 100)
    samples = jax.vmap(lambda k: conditioned_model(k))(keys)

    assert abs(jnp.mean(samples)) < 0.2

def test_log_prob_normal():
    n = dist.normal(jnp.array(0.), jnp.array(1.))
    log_prob_mean = n.log_prob(jnp.array(0.))
    tu.check_close(log_prob_mean, -0.9189385)
    log_prob_std = n.log_prob(jnp.array(1.))
    tu.check_close(log_prob_std, -1.4189385)

reason: import jax


text:
import math as math

import numba
import numpy as np

SQRT_2PI = np.sqrt(2.0 * np.pi)

@numba.vectorize(['float64(float64, float64, float64)'], nopython=True)
def gaussian_pdf(x, mean, sigma):
    '''Compute the value of a Gaussian probability density function at x with
    given mean and sigma.'''
    return math.exp(-0.5 * ((x - mean) / sigma)**2) / (sigma * SQRT_2PI)

@numba.njit(parallel=True, fastmath=True, nogil=True)
def product_kde2(eval_points, samples, bandwidths):
    n_samples, n_bandwidths = samples.shape
    n_eval = eval_points.shape[0]

    assert(eval_points.shape[1] == samples.shape[1])
    assert(eval_points.shape[1] == bandwidths.shape[0])

    results = np.zeros((n_eval,))

    for eval_ind in numba.prange(n_eval):
        eval_point = eval_points[eval_ind]
        for sample_ind in range(n_samples):
            product_kde = 1.0
            sample = samples[sample_ind]
            for bandwidth_ind in range(n_bandwidths):
                bandwidth = bandwidths[bandwidth_ind]
                product_kde *= (
                    np.exp(-0.5 *
                           ((eval_point[bandwidth_ind] - sample[bandwidth_ind])
                            / bandwidth)**2) / (bandwidth * SQRT_2PI)
                    / bandwidth)
            results[eval_ind] += product_kde
        results[eval_ind] /= n_samples

    return results

@numba.njit(parallel=True, fastmath=True)
def product_kde(eval_points, samples, bandwidths):
    result = 0
    n_samples = len(samples)

    for sample_ind in numba.prange(n_samples):
        result += np.prod(
            gaussian_pdf(eval_points, samples[sample_ind], bandwidths)
        )

    return result / (n_samples * np.prod(bandwidths))

@numba.njit(parallel=True, nogil=True)
def gauss(x, mean, sigma):
    '''Compute the value of a Gaussian probability density function at x with
    given mean and sigma.'''
    return math.exp(-0.5 * ((x - mean) / sigma)**2) / (sigma * SQRT_2PI)

@numba.njit(parallel=True, nogil=True)
def kde(test_mark, test_position, train_marks, train_position,
        place_bandwidth, mark_bandwidth, n_marks, n_position_dims,
        n_samples):
    joint_density = 0.0

    for sample_ind in range(n_samples):
        product_kde = 1.0
        for x, mean in zip(test_mark, train_marks[sample_ind]):
            product_kde *= gauss(x, mean, mark_bandwidth) / mark_bandwidth

        for x, mean in zip(test_position, train_position[sample_ind]):
            product_kde *= gauss(x, mean, mark_bandwidth) / mark_bandwidth
        joint_density += product_kde

    return joint_density / n_samples

@numba.njit(parallel=True, nogil=True)
def get_log_likelihood(marks, place_bin_centers, occupancy,
                       ground_intensity, train_marks, train_positions,
                       mean_rates, place_bandwidth=2.0, mark_bandwidth=20.0):
    '''

    Parameters
    ----------
    marks : ndarray, shape (n_time, n_marks, n_tetrodes, )
    place_bin_centers : ndarray, shape (n_bins, n_position_dims)
    bandwidths : ndarray, shape (n_marks + n_position_dims, )
    occupancy : ndarray, shape (n_bins,)
    ground_intensity : ndarray, shape (n_tetrodes, n_bins)
    train_data : list of ndarray, shape (n_samples, n_marks +
                                            n_position_dims)
    mean_rates : ndarray, shape (n_tetrodes,)

    Returns
    -------
    log_likelihood : ndarray, shape (n_time, n_bins)

    '''
    n_time, n_marks, n_tetrodes = marks.shape
    n_bins, n_position_dims = place_bin_centers.shape
    log_likelihood = np.zeros((n_time, n_bins))

    for time_ind in numba.prange(n_time):
        for tetrode_ind in numba.prange(n_tetrodes):
            is_spike = np.any(~np.isnan(marks[time_ind, :, tetrode_ind]))
            if is_spike:
                test_mark = marks[time_ind, :, tetrode_ind]
                train_mark = train_marks[tetrode_ind]
                mean_rate = np.log(mean_rates[tetrode_ind])
                n_train = train_marks[tetrode_ind].shape[0]

                for bin_ind in numba.prange(n_bins):
                    position_bin = place_bin_centers[bin_ind, :]
                    log_likelihood[time_ind, bin_ind] += np.log(kde(
                        test_mark, position_bin, train_mark,
                        train_positions[tetrode_ind], place_bandwidth,
                        mark_bandwidth, n_marks, n_position_dims, n_train))

                    log_likelihood[time_ind, bin_ind] += mean_rate
                    log_likelihood[time_ind, bin_ind] -= np.log(
                        occupancy[bin_ind] + np.spacing(1))
                    log_likelihood[time_ind, bin_ind] -= (
                        ground_intensity[tetrode_ind][bin_ind])

            else:
                for bin_ind in range(n_bins):
                    log_likelihood[time_ind, bin_ind] -= (
                        ground_intensity[tetrode_ind][bin_ind])

    return log_likelihood

reason: import numpy,import numba


text:
# _*_ coding: utf-8 _*_
__author__ = 'Di Meng'
__date__ = '2/5/2018 9:37 PM'

import tensorflow as tf
import numpy as np

x_data = np.random.rand(100).astype(np.float32)
y_data = x_data * 0.1 + 0.3

weights = tf.Variable(tf.random_uniform([1], -1.0, 1.0))
biases = tf.Variable(tf.zeros([1]))

y = weights * x_data + biases

loss = tf.reduce_mean(tf.square(y - y_data))

# choose number less than 1
optimizer = tf.train.GradientDescentOptimizer(0.5)
train = optimizer.minimize(loss)

# need to init before run
init = tf.initialize_all_variables()

session = tf.Session()
session.run(init)

for step in range(201):
    session.run(train)
    if step % 20 == 0:
        print(step, session.run(weights), session.run(biases))

max_stars_repo_name: leonsariel/python

max_issues_repo_name: leonsariel/python

max_forks_repo_name: leonsariel/python