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

# In[ ]:

# !pip install -U torch
# !pip install -U torchvision
# !pip install -U pillow==6.2.0
# !pip install -q monai
# !pip install -q git+https://github.com/ildoonet/pytorch-gradual-warmup-lr.git

# In[ ]:

import json

# Opening JSON file
with open('../settings.json') as json_file:
    settings = json.load(json_file)

from torch import nn
import numpy as np
import torch
import argparse
from torch import nn
import numpy as np
import torch
import sys
import json
my_parser = argparse.ArgumentParser(description='RSNA - Pulmonary Embolism Detection: trian 3dcnn')
my_parser.add_argument('--train',type=str)
my_parser.add_argument('--fast',type=int)
args = my_parser.parse_args()
do_train = args.train
DEBUG = True if args.fast==1 else False

CFG = {
    'image_target_cols': [
        'pe_present_on_image', # only image level
    ],
    'exam_target_cols': [
        'negative_exam_for_pe', # exam level
        'rv_lv_ratio_gte_1', # exam level
        'rv_lv_ratio_lt_1', # exam level
        'leftsided_pe', # exam level
        'chronic_pe', # exam level
        'rightsided_pe', # exam level
        'acute_and_chronic_pe', # exam level
        'central_pe', # exam level
        'indeterminate' # exam level
    ],
    'image_weight': 0.07361963,
    'exam_weights': [0.0736196319, 0.2346625767, 0.0782208589, 0.06257668712, 0.1042944785, 0.06257668712, 0.1042944785, 0.1877300613, 0.09202453988],
}
def rsna_metric(label, predicted ,bce_func = torch.nn.BCELoss(reduction='none'),CFG=CFG):

    y_pred_exam = predicted
    y_true_exam = label

    total_loss = torch.tensor(0, dtype=torch.float32).cuda()
    total_weights = torch.tensor(0, dtype=torch.float32).cuda()

    label_w = torch.tensor(CFG['exam_weights']).view(1, -1).cuda()

    exam_loss = bce_func(y_pred_exam, y_true_exam)
    exam_loss = torch.sum(exam_loss*label_w, 1)[0] # Kaggle us

    total_loss += exam_loss
    total_weights += label_w.sum()
    final_loss = total_loss.cuda()/total_weights.cuda()
    return final_loss

class RsnaLoss(nn.Module):
    def __init__(self):
        super(RsnaLoss, self).__init__()
        self.rsna_metric = rsna_metric
    def forward(self,predicted,label):
        rsna = self.rsna_metric(label, predicted)
        return rsna

# In[ ]:

import monai

# In[ ]:

import os
import sys
os.environ['CUDA_VISIBLE_DEVICES'] = '0,1,2' # specify GPUs locally

# libraries
import time
import numpy as np
import pandas as pd
import cv2
import PIL.Image
import matplotlib.pyplot as plt
from tqdm import tqdm as tqdm
from sklearn.metrics import cohen_kappa_score, confusion_matrix
from sklearn.model_selection import train_test_split, StratifiedKFold
import random

import torch
from torch.utils.data import TensorDataset, DataLoader,Dataset
import torch.nn as nn
import torch.nn.functional as F
import torchvision
import torchvision.transforms as transforms
import torch.optim as optim
from torch.optim import lr_scheduler
from torch.utils.data.sampler import SubsetRandomSampler, RandomSampler, SequentialSampler
from torch.optim.lr_scheduler import StepLR, ReduceLROnPlateau, CosineAnnealingLR
from warmup_scheduler import GradualWarmupScheduler
from sklearn.metrics import roc_auc_score
import albumentations

import monai
from monai.data import NiftiDataset
from monai.transforms import AddChannel, Compose, RandRotate90, Resize, ScaleIntensity, ToTensor

from apex import amp # I cannot install apex in Kagggle notebook

device = torch.device('cuda')

def set_seed(seed=0):
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
#     torch.backends.cudnn.deterministic = True

set_seed(7777)

# In[ ]:

if DEBUG:
    print("DEBUG mode")
# kernel_type = 'densenet121_rlc'
if not os.path.exists("../"+settings['MODEL_PATH3D']):
    os.makedirs("../"+settings['MODEL_PATH3D'])
print("../"+settings['MODEL_PATH3D'])
kernel_type = "../"+settings['MODEL_PATH3D']+'/densenet121_rlc'

image_size = 160
use_amp = False
data_dir = "../"+settings['jpeg_dir']#'../Datasets/RSNA/train256/'
num_workers = 32
init_lr = 1e-5
out_dim = 3
freeze_epo = 0
warmup_epo = 1
cosine_epo = 2 if DEBUG else 20
n_epochs = freeze_epo + warmup_epo + cosine_epo

# In[ ]:

target_cols = [
#         'negative_exam_for_pe', # exam level
#         'rv_lv_ratio_gte_1', # exam level
#         'rv_lv_ratio_lt_1', # exam level
        'leftsided_pe', # exam level
#         'chronic_pe', # exam level
        'rightsided_pe', # exam level
#         'acute_and_chronic_pe', # exam level
        'central_pe', # exam level
#         'indeterminate' # exam level
    ]

# In[ ]:

df = pd.read_csv("../"+settings['train_csv_path'])#'../Datasets/RSNA/train.csv')

# In[ ]:

# df = pd.read_csv('/media/medicalnfs/Datasets/RSNA/train.csv')
# df = df[(df.rv_lv_ratio_gte_1 == 1) | ( df.rv_lv_ratio_lt_1==1)].reset_index(drop=True)
# df.head()

# In[ ]:

from sklearn.model_selection import GroupKFold

np.random.seed(0)
group_kfold = GroupKFold(n_splits=5)
print(group_kfold)

df['fold'] = -1
for i, (_, val_index) in enumerate(group_kfold.split(df, groups=df.StudyInstanceUID)):
    df.loc[val_index, 'fold'] = i

df.fold.value_counts()

# In[ ]:

df_study = df.drop_duplicates('StudyInstanceUID')[['StudyInstanceUID','SeriesInstanceUID','fold']+target_cols]
if DEBUG:
    df_study = df_study.head(1000)

# In[ ]:

def preper(row):
    jpg_lst = sorted(glob(os.path.join(data_dir, row.StudyInstanceUID, row.SeriesInstanceUID, '*.jpg')))
    img_lst = [cv2.imread(jpg)[:,:,::-1] for jpg in jpg_lst]
    img = np.stack([image.astype(np.float32) for image in img_lst], axis=2).transpose(3,0,1,2)
    return row.StudyInstanceUID,img

# In[ ]:

# all_imgs ={}
# for index,row in tqdm(df_study.iterrows()):
#     jpg_lst = sorted(glob(os.path.join(data_dir, row.StudyInstanceUID, row.SeriesInstanceUID, '*.jpg')))
#     img_lst = [cv2.imread(jpg)[:,:,::-1] for jpg in jpg_lst]
#     img = np.stack([image.astype(np.float32) for image in img_lst], axis=2).transpose(3,0,1,2)
#     all_imgs[row.StudyInstanceUID] = img

# In[ ]:

from joblib import Parallel, delayed
from glob import glob
from monai.transforms import LoadNifti, Randomizable, apply_transform
from monai.transforms import AddChannel, Compose, RandRotate90, Resize, ScaleIntensity, ToTensor, RandAffine
from monai.utils import get_seed
from tqdm.auto import tqdm
# x = Parallel(n_jobs=32)(delayed(preper)(row) for index,row in tqdm(df_study.iterrows()))

# In[ ]:

# all_imgs ={}

# In[ ]:

# for row in x:
#     all_imgs[row[0]] = row[1]

# In[ ]:

class RSNADataset3D(torch.utils.data.Dataset, Randomizable):
    def __init__(self, csv, mode, transform=None):

        self.csv = csv.reset_index()
        self.mode = mode
        self.transform = transform

    def __len__(self):
        return self.csv.shape[0]

    def randomize(self) -> None:
        MAX_SEED = np.iinfo(np.uint32).max + 1
        self._seed = self.R.randint(MAX_SEED, dtype="uint32")

    def __getitem__(self, index):
        self.randomize()
        row = self.csv.iloc[index]
        jpg_lst = sorted(glob(os.path.join(data_dir, row.StudyInstanceUID, row.SeriesInstanceUID, '*.jpg')))
        img_lst = np.array([cv2.imread(jpg)[:,:,::-1] for jpg in jpg_lst])
#         print(img_lst.shape)
#         if np.random.randint(2)==0 and self.mode=='train':
#             img_lst = img_lst[::-1]
        if np.random.randint(2)==0 and self.mode=='train':
            z = np.random.randint(1,max(len(img_lst)//6,10))
            y = np.random.randint(3,30)
            x = np.random.randint(3,30)
            img_lst = img_lst[z:-z,x:-x,y:-y]
#             print("flip")
        img = np.stack([image.astype(np.float32) for image in img_lst], axis=2).transpose(3,0,1,2)

        if self.transform is not None:
            if isinstance(self.transform, Randomizable):
                self.transform.set_random_state(seed=self._seed)
            img = apply_transform(self.transform, img)

        if self.mode == 'test':
            return img
        else:
            return img, torch.tensor(row[target_cols]).float()

# In[ ]:

def default_collate(batch):
    data = torch.stack([item[0] for item in batch])
    target = torch.stack([item[1] for item in batch])  # image labels.
    return data, target

# In[ ]:

train_transforms = Compose([ScaleIntensity(),
                            Resize((image_size, image_size, image_size)),
                            RandAffine(
                                      prob=0.5,
#                                       rotate_range=(np.pi * 2, np.pi * 2, np.pi * 2),
                                      scale_range=(0.15, 0.15, 0.15),
                                      padding_mode='border'),
                            ToTensor()])
val_transforms = Compose([ScaleIntensity(),Resize((image_size, image_size, image_size)),ToTensor()])

# In[ ]:

dataset_show = RSNADataset3D(df_study.head(5), 'train', transform=val_transforms)
dataset_show_aug = RSNADataset3D(df_study.head(5), 'train', transform=train_transforms)
# from pylab import rcParams
# rcParams['figure.figsize'] = 20,5
# for i in range(5):
#     f, axarr = plt.subplots(1,6)
#     img, label = dataset_show[i]
#     for j in range(6):
#         if j<=2: axarr[j].imshow(img.numpy().transpose(1,2,3,0).mean(axis=j))
#         elif j==3: axarr[j].imshow(img.numpy().transpose(1,2,3,0)[image_size//2,:,:])
#         elif j==4: axarr[j].imshow(img.numpy().transpose(1,2,3,0)[:,image_size//2,:])
#         elif j==5: axarr[j].imshow(img.numpy().transpose(1,2,3,0)[:,:,image_size//2])
#         axarr[j].set_title(f"Orig {i}")
#     f, axarr = plt.subplots(1,6)
#     img, label = dataset_show_aug[i]
#     img = torch.nn.functional.interpolate(img.unsqueeze(0), size=(160,160,160))[0]
#     for j in range(6):
#         if j<=2: axarr[j].imshow(img.numpy().transpose(1,2,3,0).mean(axis=j))
#         elif j==3: axarr[j].imshow(img.numpy().transpose(1,2,3,0)[image_size//2,:,:])
#         elif j==4: axarr[j].imshow(img.numpy().transpose(1,2,3,0)[:,image_size//2,:])
#         elif j==5: axarr[j].imshow(img.numpy().transpose(1,2,3,0)[:,:,image_size//2])
#         axarr[j].set_title(f"Aug {i}")
#     plt.show()

# In[ ]:

# bce = nn.BCEWithLogitsLoss(weight= torch.tensor(CFG['exam_weights']).view(-1).cuda())
bce = nn.BCEWithLogitsLoss()
def criterion(logits, target):
    loss = bce(logits.cuda(), target.cuda())
    return loss

# In[ ]:

def train_epoch(model, loader, optimizer):

    model.train()
    train_loss = []
    bar = tqdm(loader)
    for (data, target) in bar:
        data, target = data.to(device), target.to(device)
#         data = torch.nn.functional.interpolate(data, size=(160,160,160))
        optimizer.zero_grad()
        logits = model(data)
        loss = criterion(logits, target)

        if not use_amp:
            loss.backward()
        else:
            with amp.scale_loss(loss, optimizer) as scaled_loss:
                scaled_loss.backward()

        optimizer.step()

        loss_np = loss.detach().cpu().numpy()
        train_loss.append(loss_np)
        smooth_loss = sum(train_loss[-100:]) / min(len(train_loss), 100)
        bar.set_description('loss: %.5f, smth: %.5f' % (loss_np, smooth_loss))
    return train_loss

def val_epoch(model, loader, is_ext=None, n_test=1, get_output=False):

    model.eval()
    val_loss = []
    LOGITS = []
    TARGETS = []

    with torch.no_grad():
        for (data, target) in tqdm(loader):
            data, target = data.to(device), target.to(device)
#             data = torch.nn.functional.interpolate(data, size=(160,160,160))
            logits = model(data)
            LOGITS.append(logits.detach().cpu())
            TARGETS.append(target.detach().cpu())

    val_loss = criterion(torch.cat(LOGITS), torch.cat(TARGETS)).cpu().numpy()
    PROBS = torch.sigmoid(torch.cat(LOGITS)).cpu().numpy().squeeze()
    LOGITS = torch.cat(LOGITS).cpu().numpy()
    TARGETS = torch.cat(TARGETS).cpu().numpy()

    if get_output:
        return LOGITS, PROBS, TARGETS
    else:
        acc = (PROBS.round() == TARGETS).mean() * 100.
        auc = roc_auc_score(TARGETS, LOGITS)
        return float(val_loss), acc, auc

# In[ ]:

class GradualWarmupSchedulerV2(GradualWarmupScheduler):
    def __init__(self, optimizer, multiplier, total_epoch, after_scheduler=None):
        super(GradualWarmupSchedulerV2, self).__init__(optimizer, multiplier, total_epoch, after_scheduler)
    def get_lr(self):
        if self.last_epoch > self.total_epoch:
            if self.after_scheduler:
                if not self.finished:
                    self.after_scheduler.base_lrs = [base_lr * self.multiplier for base_lr in self.base_lrs]
                    self.finished = True
                return self.after_scheduler.get_lr()
            return [base_lr * self.multiplier for base_lr in self.base_lrs]
        if self.multiplier == 1.0:
            return [base_lr * (float(self.last_epoch) / self.total_epoch) for base_lr in self.base_lrs]
        else:
            return [base_lr * ((self.multiplier - 1.) * self.last_epoch / self.total_epoch + 1.) for base_lr in self.base_lrs]

# In[ ]:

# monai.networks.nets.senet.se_resnext50_32x4d?

# In[ ]:

def run(fold):
    df_train = df_study[(df_study['fold'] != fold)]
    df_valid = df_study[(df_study['fold'] == fold)]

    dataset_train = RSNADataset3D(df_train, 'train', transform=train_transforms)
    dataset_valid = RSNADataset3D(df_valid, 'val', transform=val_transforms)
    train_loader = torch.utils.data.DataLoader(dataset_train, batch_size=9, sampler=RandomSampler(dataset_train), num_workers=num_workers)
    valid_loader = torch.utils.data.DataLoader(dataset_valid, batch_size=9, num_workers=num_workers)

#     model = monai.networks.nets.senet.se_resnext101_32x4d(spatial_dims=3, in_channels=3, num_classes=out_dim).to(device)
    model = monai.networks.nets.densenet.densenet121(spatial_dims=3, in_channels=3, out_channels=out_dim).to(device)

    val_loss_best = 1000
    model_file = f'{kernel_type}_best_fold{fold}.pth'

    optimizer = optim.Adam(model.parameters(), lr=init_lr)
    if use_amp:
        model, optimizer = amp.initialize(model, optimizer, opt_level="O1")
#     if len(os.environ['CUDA_VISIBLE_DEVICES'].split(',')) > 1:
    model = nn.DataParallel(model)
#     if os.path.exists(model_file):
#         model.load_state_dict(torch.load(model_file))
#         print("load")
#     if fold==1:
#         model.load_state_dict(torch.load('densenet121_best_fold1.pth'))
#         print("load")
    scheduler_cosine = torch.optim.lr_scheduler.CosineAnnealingWarmRestarts(optimizer, cosine_epo)
    scheduler_warmup = GradualWarmupSchedulerV2(optimizer, multiplier=10, total_epoch=warmup_epo, after_scheduler=scheduler_cosine)

    print(len(dataset_train), len(dataset_valid))

    for epoch in range(1, n_epochs+1):
        print(time.ctime(), 'Epoch:', epoch)
        scheduler_warmup.step(epoch-1)

        train_loss = train_epoch(model, train_loader, optimizer)
        val_loss, acc, auc = val_epoch(model, valid_loader)

        content = time.ctime() + ' ' + f'Fold {fold}, Epoch {epoch}, lr: {optimizer.param_groups[0]["lr"]:.7f}, train loss: {np.mean(train_loss):.5f}, valid loss: {(val_loss):.5f}, acc: {(acc):.4f}, auc: {(auc):.6f}'
        print(content)
#         with open(f'log_{kernel_type}.txt', 'a') as appender:
#             appender.write(content + '\n')

        if val_loss < val_loss_best:
            print('val_loss_best ({:.6f} --> {:.6f}).  Saving model ...'.format(val_loss_best, val_loss))
            torch.save(model.state_dict(), model_file)
            val_loss_best = val_loss

    torch.save(model.state_dict(), f'{kernel_type}_model_fold{fold}.pth')

if do_train == 'train':
    run(fold=1)
    run(fold=2)
    run(fold=0)
    run(fold=3)
    run(fold=4)

from sklearn.metrics import confusion_matrix

# In[ ]:

# if not os.path.exists(f"features_{kernel_type}"):
#     os.makedirs(f"features_{kernel_type}")
if not os.path.exists("../"+settings['features_densenet121_rlc']):
    os.makedirs("../"+settings['features_densenet121_rlc'])

fdir_f = "../"+settings['features_densenet121_rlc']
print(fdir_f)

# In[ ]:

torch.multiprocessing.set_sharing_strategy('file_system')

# In[ ]:

class RSNADataset3D(torch.utils.data.Dataset, Randomizable):
    def __init__(self, csv, mode, transform=None):

        self.csv = csv.reset_index()
        self.mode = mode
        self.transform = transform

    def __len__(self):
        return self.csv.shape[0]

    def randomize(self) -> None:
        MAX_SEED = np.iinfo(np.uint32).max + 1
        self._seed = self.R.randint(MAX_SEED, dtype="uint32")

    def __getitem__(self, index):
        self.randomize()
        row = self.csv.iloc[index]
        jpg_lst = sorted(glob(os.path.join(data_dir, row.StudyInstanceUID, row.SeriesInstanceUID, '*.jpg')))
        img_lst = [cv2.imread(jpg)[:,:,::-1] for jpg in jpg_lst]
        img = np.stack([image.astype(np.float32) for image in img_lst], axis=2).transpose(3,0,1,2)

        if self.transform is not None:
            if isinstance(self.transform, Randomizable):
                self.transform.set_random_state(seed=self._seed)
            img = apply_transform(self.transform, img)

        if self.mode == 'test':
            return img
        else:
            return img, torch.tensor(row[target_cols]).float(),row.StudyInstanceUID

# In[ ]:

for fold in range(0,5):
    df_valid = df_study[(df_study['fold'] == fold)]
    dataset_valid = RSNADataset3D(df_valid, 'val', transform=val_transforms)
    valid_loader = torch.utils.data.DataLoader(dataset_valid, batch_size=1, num_workers=num_workers)
    def load_model(model_file):
        model = monai.networks.nets.densenet.densenet121(spatial_dims=3, in_channels=3, out_channels=out_dim).to(device)

        try:  # single GPU model_file
            model.load_state_dict(torch.load(model_file), strict=True)
        except:  # multi GPU model_file
            state_dict = torch.load(model_file)
            state_dict = {k[7:] if k.startswith('module.') else k: state_dict[k] for k in state_dict.keys()}
            model.load_state_dict(state_dict, strict=True)

        model.eval()
        print()
        return model

    LOGITS = []
    PROBS = []
    try:
        model = load_model(f'{kernel_type}_best_fold{fold}.pth')
    except:
        model = load_model(f'{kernel_type}_model_fold{fold}.pth')
    model.eval()
    with torch.no_grad():
        for data,target,name in tqdm(valid_loader):
            data = data.to(device)
            target = target.numpy().reshape(-1)
            l1 = model(data)
            l = torch.sigmoid(l1)
            LOGITS.append(l.detach().cpu().numpy().reshape(-1))
            PROBS.append(target)
            l2 = model.features(data).cpu().numpy()[0]
            np.save(f"{fdir_f}/{name[0]}_3dcnn.npy",l2)
            np.save(f"{fdir_f}/{name[0]}_3dprob.npy",l1.cpu().numpy().reshape(-1))
    print(fold)
    print(confusion_matrix(np.array(PROBS)[:,0]>0.5,np.array(LOGITS)[:,0]>0.5))
    print(confusion_matrix(np.array(PROBS)[:,1]>0.5,np.array(LOGITS)[:,1]>0.5))
    print(confusion_matrix(np.array(PROBS)[:,2]>0.5,np.array(LOGITS)[:,2]>0.5))

# In[ ]:

reason: import numpy


text:
"""
Reads a hierarchical data format (h5) file containing a complexe frequency spectrogram and generates the audio file.
"""

__author__ = 'David Flury'
__email__ = "david@flury.email"

import os
import h5py
import librosa
import argparse
import numpy as np

def generate_audiofile(spectrograms, name, path, fft_window, sample_rate):
    file = os.path.join(path, name + ".wav")
    data = [librosa.istft((generate_stft(spectrogram))) for spectrogram in spectrograms]
    audio = np.array(data)
    print('Output audio file: %s' % file)
    librosa.output.write_wav(file, audio[0] if len(data) == 1 else audio, sample_rate, norm=False)

def generate_stft(spectrogram):
    real_part = spectrogram[:, :, 0]
    imag_part = spectrogram[:, :, 1]
    stft = real_part + imag_part * 1j
    return stft

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Reads h5 a file containing a complexe frequency spectrogram and generates the audio file.')
    parser.add_argument('--file', default='D:\\Data\\test_abs-angle-instrumental.h5',
        type=str, help='Hierarhical data format container file')
    parser.add_argument('--path', default='', type=str, help='Optional output folder')

    args = parser.parse_args()
    print('Arguments:', str(args))

    if not args.path:
        args.path = os.path.dirname(args.file)

    h5f = h5py.File(args.file,'r')
    file = h5f['file'].value
    fft_window = h5f['fft_window'].value
    sample_rate = h5f['sample_rate'].value
    channels = h5f['channels'].value

    spectrograms = []
    if channels > 1:
        spectrograms.append(h5f['spectrogram_left'][:, :, :])
        spectrograms.append(h5f['spectrogram_right'][:, :, :])
    else:
        spectrograms.append(h5f['spectrogram'][:, :, :])

    generate_audiofile(spectrograms, file, args.path, fft_window, sample_rate)

    h5f.close()
    print('Finished processing')

ext: py

max_stars_repo_name: unmix-io/tools

max_issues_repo_name: unmix-io/tools

max_forks_repo_name: unmix-io/tools


text:
# -*- coding: utf-8 -*-
"""
=== DFT_pyaudio_basics.py =================================================

 Demo für framebasierte DFT

 Eine Audio-Datei wird blockweise eingelesen, in numpy-Arrays umgewandelt
 dann wird die DFT von linkem und rechten Kanal berechnet, ein Teil der
 DFT-Bins wird zu Null gesetzt und davon die inverse DFT berechnet. Das
 Ergebnis wird wieder als Audiostream ausgegeben.

===========================================================================
"""
from __future__ import division, print_function, unicode_literals
import numpy as np
import os

from numpy import (pi, log10, exp, sqrt, sin, cos, tan, angle, arange,
                    linspace, array, zeros, ones)
from numpy.fft import fft, ifft, fftshift, ifftshift, fftfreq

import matplotlib.pyplot as plt
from matplotlib.pyplot import (figure, plot, stem, grid, xlabel, ylabel,
    subplot, title, clf, xlim, ylim)

import pyaudio
import wave

np_type = np.int16 # data type for audio samples
CHUNK = 256 # number of samples in one frame

# path = '/home/muenker/Daten/share/Musi/wav/'
path = '../_media/'
# filename = 'Ole_16bit.wav'
filename = 'SpaceRipple.wav'

wf = wave.open(os.path.join(path, filename))
n_chan = wf.getnchannels() # number of channels in wav-file
w_samp = wf.getsampwidth() # wordlength of samples
rate_in = wf.getframerate() # samplerate in wav-file

print("Channels:", n_chan, "\nSample width:",w_samp,"bytes\nSample rate:",rate_in)

wf = wave.open(os.path.join(path, filename))
p = pyaudio.PyAudio() # instantiate PyAudio + setup PortAudio system

# open a stream on the desired device with the desired audio parameters
# for reading or writing
stream = p.open(format=p.get_format_from_width(wf.getsampwidth()),
                channels=wf.getnchannels(),
                rate=wf.getframerate(),
                output=True)

# initialize arrays for samples
samples_in = samples_out = zeros(CHUNK*2, dtype=np_type) # stereo
samples_l  = samples_r = zeros(CHUNK, dtype=np_type) # mono

data_out = 'dummy'

while data_out:

# read CHUNK frames to string, convert to numpy array and split in R / L chan.:
# R / L samples are interleaved, each sample is 16 bit = 2 Bytes
    samples_in = np.fromstring(wf.readframes(CHUNK), dtype=np_type)

## Example for dtype = np.int8 (8 bits) = 1 ndarray element,
##    two consecutive bytes / ndarray elements = 1 sample
## Split signal into L and R channel:
#    samples_l[0::2] = samples_in[0::4]
#    samples_l[1::2] = samples_in[1::4]
#    samples_r[0::2] = samples_in[2::4]
#    samples_r[1::2] = samples_in[3::4]
#
## Do some numpy magic with samples_l and samples_r
#               ...
#
#  And re-combine L and R channel:
#    samples_out[0::4] = samples_l[0::2]
#    samples_out[1::4] = samples_l[1::2]
#    samples_out[2::4] = samples_r[0::2]
#    samples_out[3::4] = samples_r[1::2]
#---------------------------------------------------------------------------
## dtype = np.int16 (16 bits): 1 ndarray element = 1 sample :
    samples_l = samples_in[0::2]
    samples_r = samples_in[1::2]
    if len(samples_r) < 2:
        break # break out of the while loop when out of data
    # Check whether there was enough data for a full frame
    if len(samples_r) < CHUNK: # check whether frame has full length
        samples_out = samples_np = zeros(len(samples_in), dtype=np_type)
        samples_l = samples_l = zeros(len(samples_in)/2, dtype=np_type)

# ---- Numpy Magic happens here (suppress all FFT bins > 63) ---------------

    fft_l = fft(samples_in[0::2]) # convert to frequency domain
    fft_r = fft(samples_in[1::2])
    fft_l[64:-64] = 0 # set NFFT-bins between 63 ... NFFT-63 = 0
    fft_r[64:-64] = 0 # maintaining  symmetric spectrum -> real-valued time signal

    ifft_l = ifft(fft_l).astype(np_type).real # delete imaginary part, abs() doesn't
    ifft_r = ifft(fft_r).astype(np_type).real #   work with bipolar signals

    samples_out[0::2] = ifft_l # convert back to time domain
    samples_out[1::2] = ifft_r

#    data_out = np.chararray.tostring(samples_np.astype(np_type)) # convert back to string
    data_out = np.chararray.tostring(samples_out) # convert back to string
    stream.write(data_out) # play audio by writing audio data to the stream (blocking)

stream.stop_stream() # pause audio stream
stream.close() # close audio stream

p.terminate() # close PyAudio & terminate PortAudio system

# see: http://stackoverflow.com/questions/23370556/recording-24-bit-audio-with-pyaudio
# http://stackoverflow.com/questions/16767248/how-do-i-write-a-24-bit-wav-file-in-python?

reason: import numpy,from numpy


text:
import numpy as np
import tensorflow.keras.backend as K

# Use custom metrics for training and evalution

# Recall
def micro_recall(y_true, y_pred):
  y_true_and_pred = K.cast(y_true * y_pred, 'float')
  true_pos = K.sum(y_true_and_pred, axis=0)
  possible_pos = K.sum(y_true, axis=0)
  rec = K.sum(true_pos) / (K.sum(possible_pos) + K.epsilon())
  return rec

def macro_recall(y_true, y_pred):
  y_true_and_pred = K.cast(y_true * y_pred, 'float')
  true_pos = K.sum(y_true_and_pred, axis=0)
  possible_pos = K.sum(y_true, axis=0)
  rec = K.mean(true_pos / (possible_pos + K.epsilon()))
  return rec

def weighted_recall(y_true, y_pred):
  y_true_and_pred = K.cast(y_true * y_pred, 'float')
  true_pos = K.sum(y_true_and_pred, axis=0)
  possible_pos = K.sum(y_true, axis=0)
  wts = K.sum(y_true, axis=0) / K.sum(y_true)
  rec = K.mean((true_pos * wts) / (possible_pos + K.epsilon()))
  return rec

# Precision
def micro_precision(y_true, y_pred):
  y_true_and_pred = K.cast(y_true * y_pred, 'float')
  true_pos = K.sum(y_true_and_pred, axis=0)
  pred_pos = K.sum(y_pred, axis=0)
  prec = K.sum(true_pos) / (K.sum(pred_pos) + K.epsilon())
  return prec

def macro_precision(y_true, y_pred):
  y_true_and_pred = K.cast(y_true * y_pred, 'float')
  true_pos = K.sum(y_true_and_pred, axis=0)
  pred_pos = K.sum(y_pred, axis=0)
  prec = K.mean(true_pos / (pred_pos + K.epsilon()))
  return prec

def weighted_precision(y_true, y_pred):
  y_true_and_pred = K.cast(y_true * y_pred, 'float')
  true_pos = K.sum(y_true_and_pred, axis=0)
  pred_pos = K.sum(y_pred, axis=0)
  wts = K.sum(y_true, axis=0) / K.sum(y_true)
  rec = K.mean((true_pos * wts) / (pred_pos + K.epsilon()))
  return rec

#F1 score

# get one-hot rep of given tensor
def get_one_hot(y):
  return K.one_hot(K.argmax(y,axis=1), y.shape[1])

def fbeta(prec, rec, beta=1.0):
  f_score = ((1+beta) * prec * rec)/(beta*prec + rec + K.epsilon())
  return f_score

def micro_f1(y_true, y_pred):
  y_pred = get_one_hot(y_pred)
  prec = micro_precision(y_true, y_pred)
  rec = micro_recall(y_true, y_pred)
  f1 = fbeta(prec, rec)
  return f1

def macro_f1(y_true, y_pred):
  y_pred = get_one_hot(y_pred)
  prec = macro_precision(y_true, y_pred)
  rec = macro_recall(y_true, y_pred)
  f1 = fbeta(prec, rec)
  return f1

def weighted_f1(y_true, y_pred):
  y_pred = get_one_hot(y_pred)
  prec = weighted_precision(y_true, y_pred)
  rec = weighted_recall(y_true, y_pred)
  f1 = fbeta(prec, rec)
  return f1

max_stars_repo_name: wadpac/SleepStageClassification

max_issues_repo_name: wadpac/SleepStageClassification

max_forks_repo_name: wadpac/SleepStageClassification


text:
from __future__ import absolute_import
from collections import deque, namedtuple
import warnings
import random
import numpy as np
from rl.util import SumSegmentTree, MinSegmentTree

# This is to be understood as a transition: Given `state0`, performing `action`
# yields `reward` and results in `state1`, which might be `terminal`.
Experience = namedtuple('Experience', 'state0, action, reward, state1, terminal1')

def sample_batch_indexes(low, high, size):
    """Return a sample of (size) unique elements between low and high
        # Argument
            low (int): The minimum value for our samples
            high (int): The maximum value for our samples
            size (int): The number of samples to pick
        # Returns
            A list of samples of length size, with values between low and high
        """
    if high - low >= size:
        # We have enough data. Draw without replacement, that is each index is unique in the
        # batch. We cannot use `np.random.choice` here because it is horribly inefficient as
        # the memory grows. See https://github.com/numpy/numpy/issues/2764 for a discussion.
        # `random.sample` does the same thing (drawing without replacement) and is way faster.
        try:
            r = xrange(low, high)
        except NameError:
            r = range(low, high)
        batch_idxs = random.sample(r, size)
    else:
        # Not enough data. Help ourselves with sampling from the range, but the same index
        # can occur multiple times. This is not good and should be avoided by picking a
        # large enough warm-up phase.
        warnings.warn('Not enough entries to sample without replacement. Consider increasing your warm-up phase to avoid oversampling!')
        batch_idxs = np.random.random_integers(low, high - 1, size=size)
    assert len(batch_idxs) == size
    return batch_idxs

class RingBuffer(object):
    def __init__(self, maxlen):
        self.maxlen = maxlen
        self.start = 0
        self.length = 0
        self.data = [None for _ in range(maxlen)]

    def __len__(self):
        return self.length

    def __getitem__(self, idx):
        """Return element of buffer at specific index
        # Argument
            idx (int): Index wanted
        # Returns
            The element of buffer at given index
        """
        if idx < 0 or idx >= self.length:
            raise KeyError()
        return self.data[(self.start + idx) % self.maxlen]

    def append(self, v):
        """Append an element to the buffer
        # Argument
            v (object): Element to append
        """
        if self.length < self.maxlen:
            # We have space, simply increase the length.
            self.length += 1
        elif self.length == self.maxlen:
            # No space, "remove" the first item.
            self.start = (self.start + 1) % self.maxlen
        else:
            # This should never happen.
            raise RuntimeError()
        self.data[(self.start + self.length - 1) % self.maxlen] = v

def zeroed_observation(observation):
    """Return an array of zeros with same shape as given observation
    # Argument
        observation (list): List of observation
    # Return
        A np.ndarray of zeros with observation.shape
    """
    if hasattr(observation, 'shape'):
        return np.zeros(observation.shape)
    elif hasattr(observation, '__iter__'):
        out = []
        for x in observation:
            out.append(zeroed_observation(x))
        return out
    else:
        return 0.

class Memory(object):
    def __init__(self, window_length, ignore_episode_boundaries=False):
        self.window_length = window_length
        self.ignore_episode_boundaries = ignore_episode_boundaries

        self.recent_observations = deque(maxlen=window_length)
        self.recent_terminals = deque(maxlen=window_length)

    def sample(self, batch_size, batch_idxs=None):
        raise NotImplementedError()

    def append(self, observation, action, reward, terminal, training=True):
        self.recent_observations.append(observation)
        self.recent_terminals.append(terminal)

    def get_recent_state(self, current_observation):
        """Return list of last observations
        # Argument
            current_observation (object): Last observation
        # Returns
            A list of the last observations
        """
        # This code is slightly complicated by the fact that subsequent observations might be
        # from different episodes. We ensure that an experience never spans multiple episodes.
        # This is probably not that important in practice but it seems cleaner.
        state = [current_observation]
        idx = len(self.recent_observations) - 1
        for offset in range(0, self.window_length - 1):
            current_idx = idx - offset
            current_terminal = self.recent_terminals[current_idx - 1] if current_idx - 1 >= 0 else False
            if current_idx < 0 or (not self.ignore_episode_boundaries and current_terminal):
                # The previously handled observation was terminal, don't add the current one.
                # Otherwise we would leak into a different episode.
                break
            state.insert(0, self.recent_observations[current_idx])
        while len(state) < self.window_length:
            state.insert(0, zeroed_observation(state[0]))
        return state

    def get_config(self):
        """Return configuration (window_length, ignore_episode_boundaries) for Memory
        # Return
            A dict with keys window_length and ignore_episode_boundaries
        """
        config = {
            'window_length': self.window_length,
            'ignore_episode_boundaries': self.ignore_episode_boundaries,
        }
        return config

class SequentialMemory(Memory):
    def __init__(self, limit, **kwargs):
        super(SequentialMemory, self).__init__(**kwargs)

        self.limit = limit

        # Do not use deque to implement the memory. This data structure may seem convenient but
        # it is way too slow on random access. Instead, we use our own ring buffer implementation.
        self.actions = RingBuffer(limit)
        self.rewards = RingBuffer(limit)
        self.terminals = RingBuffer(limit)
        self.observations = RingBuffer(limit)

    def sample(self, batch_size, batch_idxs=None):
        """Return a randomized batch of experiences
        # Argument
            batch_size (int): Size of the all batch
            batch_idxs (int): Indexes to extract
        # Returns
            A list of experiences randomly selected
        """
        # It is not possible to tell whether the first state in the memory is terminal, because it
        # would require access to the "terminal" flag associated to the previous state. As a result
        # we will never return this first state (only using `self.terminals[0]` to know whether the
        # second state is terminal).
        # In addition we need enough entries to fill the desired window length.
        assert self.nb_entries >= self.window_length + 2, 'not enough entries in the memory'

        if batch_idxs is None:
            # Draw random indexes such that we have enough entries before each index to fill the
            # desired window length.
            batch_idxs = sample_batch_indexes(
                self.window_length, self.nb_entries - 1, size=batch_size)
        batch_idxs = np.array(batch_idxs) + 1
        assert np.min(batch_idxs) >= self.window_length + 1
        assert np.max(batch_idxs) < self.nb_entries
        assert len(batch_idxs) == batch_size

        # Create experiences
        experiences = []
        for idx in batch_idxs:
            terminal0 = self.terminals[idx - 2]
            while terminal0:
                # Skip this transition because the environment was reset here. Select a new, random
                # transition and use this instead. This may cause the batch to contain the same
                # transition twice.
                idx = sample_batch_indexes(self.window_length + 1, self.nb_entries, size=1)[0]
                terminal0 = self.terminals[idx - 2]
            assert self.window_length + 1 <= idx < self.nb_entries

            # This code is slightly complicated by the fact that subsequent observations might be
            # from different episodes. We ensure that an experience never spans multiple episodes.
            # This is probably not that important in practice but it seems cleaner.
            state0 = [self.observations[idx - 1]]
            for offset in range(0, self.window_length - 1):
                current_idx = idx - 2 - offset
                assert current_idx >= 1
                current_terminal = self.terminals[current_idx - 1]
                if current_terminal and not self.ignore_episode_boundaries:
                    # The previously handled observation was terminal, don't add the current one.
                    # Otherwise we would leak into a different episode.
                    break
                state0.insert(0, self.observations[current_idx])
            while len(state0) < self.window_length:
                state0.insert(0, zeroed_observation(state0[0]))
            action = self.actions[idx - 1]
            reward = self.rewards[idx - 1]
            terminal1 = self.terminals[idx - 1]

            # Okay, now we need to create the follow-up state. This is state0 shifted on timestep
            # to the right. Again, we need to be careful to not include an observation from the next
            # episode if the last state is terminal.
            state1 = [np.copy(x) for x in state0[1:]]
            state1.append(self.observations[idx])

            assert len(state0) == self.window_length
            assert len(state1) == len(state0)
            experiences.append(Experience(state0=state0, action=action, reward=reward,
                                          state1=state1, terminal1=terminal1))
        assert len(experiences) == batch_size
        return experiences

    def append(self, observation, action, reward, terminal, training=True):
        """Append an observation to the memory
        # Argument
            observation (dict): Observation returned by environment
            action (int): Action taken to obtain this observation
            reward (float): Reward obtained by taking this action
            terminal (boolean): Is the state terminal
        """
        super(SequentialMemory, self).append(observation, action, reward, terminal, training=training)

        # This needs to be understood as follows: in `observation`, take `action`, obtain `reward`
        # and weather the next state is `terminal` or not.
        if training:
            self.observations.append(observation)
            self.actions.append(action)
            self.rewards.append(reward)
            self.terminals.append(terminal)

    @property
    def nb_entries(self):
        """Return number of observations
        # Returns
            Number of observations
        """
        return len(self.observations)

    def get_config(self):
        """Return configurations of SequentialMemory
        # Returns
            Dict of config
        """
        config = super(SequentialMemory, self).get_config()
        config['limit'] = self.limit
        return config

class EpisodeParameterMemory(Memory):
    def __init__(self, limit, **kwargs):
        super(EpisodeParameterMemory, self).__init__(**kwargs)
        self.limit = limit

        self.params = RingBuffer(limit)
        self.intermediate_rewards = []
        self.total_rewards = RingBuffer(limit)

    def sample(self, batch_size, batch_idxs=None):
        """Return a randomized batch of params and rewards
        # Argument
            batch_size (int): Size of the all batch
            batch_idxs (int): Indexes to extract
        # Returns
            A list of params randomly selected and a list of associated rewards
        """
        if batch_idxs is None:
            batch_idxs = sample_batch_indexes(0, self.nb_entries, size=batch_size)
        assert len(batch_idxs) == batch_size

        batch_params = []
        batch_total_rewards = []
        for idx in batch_idxs:
            batch_params.append(self.params[idx])
            batch_total_rewards.append(self.total_rewards[idx])
        return batch_params, batch_total_rewards

    def append(self, observation, action, reward, terminal, training=True):
        """Append a reward to the memory
        # Argument
            observation (dict): Observation returned by environment
            action (int): Action taken to obtain this observation
            reward (float): Reward obtained by taking this action
            terminal (boolean): Is the state terminal
        """
        super(EpisodeParameterMemory, self).append(observation, action, reward, terminal, training=training)
        if training:
            self.intermediate_rewards.append(reward)

    def finalize_episode(self, params):
        """Append an observation to the memory
        # Argument
            observation (dict): Observation returned by environment
            action (int): Action taken to obtain this observation
            reward (float): Reward obtained by taking this action
            terminal (boolean): Is the state terminal
        """
        total_reward = sum(self.intermediate_rewards)
        self.total_rewards.append(total_reward)
        self.params.append(params)
        self.intermediate_rewards = []

    @property
    def nb_entries(self):
        """Return number of episode rewards
        # Returns
            Number of episode rewards
        """
        return len(self.total_rewards)

    def get_config(self):
        """Return configurations of SequentialMemory
        # Returns
            Dict of config
        """
        config = super(SequentialMemory, self).get_config()
        config['limit'] = self.limit
        return config

class PartitionedRingBuffer(object):
    """
    Buffer with a section that can be sampled from but never overwritten.
    Used for demonstration data (DQfD). Can be used without a partition,
    where it would function as a fixed-idxs variant of RingBuffer.
    """
    def __init__(self, maxlen):
        self.maxlen = maxlen
        self.length = 0
        self.data = [None for _ in range(maxlen)]
        self.permanent_idx = 0
        self.next_idx = 0

    def __len__(self):
        return self.length

    def __getitem__(self, idx):
        if idx < 0:
            raise KeyError()
        return self.data[idx % self.maxlen]

    def append(self, v):
        if self.length < self.maxlen:
            self.length += 1
        self.data[(self.permanent_idx + self.next_idx)] = v
        self.next_idx = (self.next_idx + 1) % (self.maxlen - self.permanent_idx)

    def load(self, load_data):
        assert len(load_data) < self.maxlen, "Must leave space to write new data."
        for idx, data in enumerate(load_data):
            self.length += 1
            self.data[idx] = data
            self.permanent_idx += 1

class PrioritizedMemory(Memory):
    def __init__(self, limit, alpha=.4, start_beta=1., end_beta=1., steps_annealed=1, **kwargs):
        super(PrioritizedMemory, self).__init__(**kwargs)

        #The capacity of the replay buffer
        self.limit = limit

        #Transitions are stored in individual RingBuffers, similar to the SequentialMemory.
        self.actions = PartitionedRingBuffer(limit)
        self.rewards = PartitionedRingBuffer(limit)
        self.terminals = PartitionedRingBuffer(limit)
        self.observations = PartitionedRingBuffer(limit)

        assert alpha >= 0
        #how aggressively to sample based on TD error
        self.alpha = alpha
        #how aggressively to compensate for that sampling. This value is typically annealed
        #to stabilize training as the model converges (beta of 1.0 fully compensates for TD-prioritized sampling).
        self.start_beta = start_beta
        self.end_beta = end_beta
        self.steps_annealed = steps_annealed

        #SegmentTrees need a leaf count that is a power of 2
        tree_capacity = 1
        while tree_capacity < self.limit:
            tree_capacity *= 2

        #Create SegmentTrees with this capacity
        self.sum_tree = SumSegmentTree(tree_capacity)
        self.min_tree = MinSegmentTree(tree_capacity)
        self.max_priority = 1.

        #wrapping index for interacting with the trees
        self.next_index = 0

    def append(self, observation, action, reward, terminal, training=True):\
        #super() call adds to the deques that hold the most recent info, which is fed to the agent
        #on agent.forward()
        super(PrioritizedMemory, self).append(observation, action, reward, terminal, training=training)
        if training:
            self.observations.append(observation)
            self.actions.append(action)
            self.rewards.append(reward)
            self.terminals.append(terminal)
            #The priority of each new transition is set to the maximum
            self.sum_tree[self.next_index] = self.max_priority ** self.alpha
            self.min_tree[self.next_index] = self.max_priority ** self.alpha

            #shift tree pointer index to keep it in sync with RingBuffers
            self.next_index = (self.next_index + 1) % self.limit

    def _sample_proportional(self, batch_size):
        #outputs a list of idxs to sample, based on their priorities.
        idxs = list()

        for _ in range(batch_size):
            mass = random.random() * self.sum_tree.sum(0, self.limit - 1)
            idx = self.sum_tree.find_prefixsum_idx(mass)
            idxs.append(idx)

        return idxs

    def sample(self, batch_size, beta=1.):
        idxs = self._sample_proportional(batch_size)

        #importance sampling weights are a stability measure
        importance_weights = list()

        #The lowest-priority experience defines the maximum importance sampling weight
        prob_min = self.min_tree.min() / self.sum_tree.sum()
        max_importance_weight = (prob_min * self.nb_entries)  ** (-beta)
        obs_t, act_t, rews, obs_t1, dones = [], [], [], [], []

        experiences = list()
        for idx in idxs:
            while idx < self.window_length + 1:
                idx += 1

            terminal0 = self.terminals[idx - 2]
            while terminal0:
                # Skip this transition because the environment was reset here. Select a new, random
                # transition and use this instead. This may cause the batch to contain the same
                # transition twice.
                idx = sample_batch_indexes(self.window_length + 1, self.nb_entries, size=1)[0]
                terminal0 = self.terminals[idx - 2]

            assert self.window_length + 1 <= idx < self.nb_entries

            #probability of sampling transition is the priority of the transition over the sum of all priorities
            prob_sample = self.sum_tree[idx] / self.sum_tree.sum()
            importance_weight = (prob_sample * self.nb_entries) ** (-beta)
            #normalize weights according to the maximum value
            importance_weights.append(importance_weight/max_importance_weight)

            # Code for assembling stacks of observations and dealing with episode boundaries is borrowed from
            # SequentialMemory
            state0 = [self.observations[idx - 1]]
            for offset in range(0, self.window_length - 1):
                current_idx = idx - 2 - offset
                assert current_idx >= 1
                current_terminal = self.terminals[current_idx - 1]
                if current_terminal and not self.ignore_episode_boundaries:
                    # The previously handled observation was terminal, don't add the current one.
                    # Otherwise we would leak into a different episode.
                    break
                state0.insert(0, self.observations[current_idx])
            while len(state0) < self.window_length:
                state0.insert(0, zeroed_observation(state0[0]))
            action = self.actions[idx - 1]
            reward = self.rewards[idx - 1]
            terminal1 = self.terminals[idx - 1]
            state1 = [np.copy(x) for x in state0[1:]]
            state1.append(self.observations[idx])

            assert len(state0) == self.window_length
            assert len(state1) == len(state0)
            experiences.append(Experience(state0=state0, action=action, reward=reward,
                                          state1=state1, terminal1=terminal1))
        assert len(experiences) == batch_size

        # Return a tuple whre the first batch_size items are the transititions
        # while -2 is the importance weights of those transitions and -1 is
        # the idxs of the buffer (so that we can update priorities later)
        return tuple(list(experiences)+ [importance_weights, idxs])

    def update_priorities(self, idxs, priorities):
        #adjust priorities based on new TD error
        for i, idx in enumerate(idxs):
            assert 0 <= idx < self.limit
            priority = priorities[i] ** self.alpha
            self.sum_tree[idx] = priority
            self.min_tree[idx] = priority
            self.max_priority = max(self.max_priority, priority)

    def calculate_beta(self, current_step):
        a = float(self.end_beta - self.start_beta) / float(self.steps_annealed)
        b = float(self.start_beta)
        current_beta = min(self.end_beta, a * float(current_step) + b)
        return current_beta

    def get_config(self):
        config = super(PrioritizedMemory, self).get_config()
        config['alpha'] = self.alpha
        config['start_beta'] = self.start_beta
        config['end_beta'] = self.end_beta
        config['beta_steps_annealed'] = self.steps_annealed

    @property
    def nb_entries(self):
        """Return number of observations
        # Returns
            Number of observations
        """
        return len(self.observations)

class PartitionedMemory(Memory):
    def __init__(self, limit, pre_load_data, alpha=.4, start_beta=1., end_beta=1., steps_annealed=1, **kwargs):
        super(PartitionedMemory, self).__init__(**kwargs)

        #The capacity of the replay buffer
        self.limit = limit

        #Transitions are stored in individual PartitionedRingBuffers.
        self.actions = PartitionedRingBuffer(limit)
        self.rewards = PartitionedRingBuffer(limit)
        self.terminals = PartitionedRingBuffer(limit)
        self.observations = PartitionedRingBuffer(limit)

        assert alpha >= 0
        #how aggressively to sample based on TD error
        self.alpha = alpha
        #how aggressively to compensate for that sampling.
        self.start_beta = start_beta
        self.end_beta = end_beta
        self.steps_annealed = steps_annealed

        #SegmentTrees need a leaf count that is a power of 2
        tree_capacity = 1
        while tree_capacity < self.limit:
            tree_capacity *= 2

        #Create SegmentTrees with this capacity
        self.sum_tree = SumSegmentTree(tree_capacity)
        self.min_tree = MinSegmentTree(tree_capacity)
        self.max_priority = 1.

        #unpack the expert transitions (assumes order recorded by the rl.utils.record_demo_data() method)
        demo_obs, demo_acts, demo_rews, demo_ts = [], [], [], []
        self.pre_load_data = pre_load_data
        for demo in self.pre_load_data:
            demo_obs.append(demo[0])
            demo_acts.append(demo[1])
            demo_rews.append(demo[2])
            demo_ts.append(demo[3])

        #pre-load the demonstration data
        self.observations.load(demo_obs)
        self.actions.load(demo_acts)
        self.rewards.load(demo_rews)
        self.terminals.load(demo_ts)

        self.permanent_idx = self.observations.permanent_idx
        assert self.permanent_idx == self.rewards.permanent_idx

        self.next_index = 0

        for idx in range(self.permanent_idx):
            self.sum_tree[idx] = (self.max_priority ** self.alpha)
            self.min_tree[idx] = (self.max_priority ** self.alpha)

    def append(self, observation, action, reward, terminal, training=True):
        #super() call adds to the deques that hold the most recent info, which is fed to the agent
        #on agent.forward()
        super(PartitionedMemory, self).append(observation, action, reward, terminal, training=training)
        if training:
            self.observations.append(observation)
            self.actions.append(action)
            self.rewards.append(reward)
            self.terminals.append(terminal)
            #The priority of each new transition is set to the maximum
            self.sum_tree[self.next_index + self.permanent_idx] = self.max_priority ** self.alpha
            self.min_tree[self.next_index + self.permanent_idx] = self.max_priority ** self.alpha
            #shift tree pointer index to keep it in sync with RingBuffers
            self.next_index = ((self.next_index + 1) % (self.limit - self.permanent_idx))

    def sample_proportional(self, batch_size):
        """
        Outputs a list of idxs to sample, based on their priorities.

        This function is public in this memory (vs. private in Sequential and
        Prioritized), because DQfD needs to be able to sample the same idxs
        twice (single step and n-step).
        """
        idxs = list()

        for _ in range(batch_size):
            mass = random.random() * self.sum_tree.sum(0, self.limit - 1)
            idx = self.sum_tree.find_prefixsum_idx(mass)
            idxs.append(idx)

        return idxs

    def sample_by_idxs(self, idxs, batch_size, beta=1., nstep=1, gamma=1):
        """
        Gathers transition data from the ring buffers. The PartitionedMemory
        separates generating the idxs and returning their transitions, allowing
        this method to be called multiple times with the same idxs.
        """
        #importance sampling weights are a stability measure
        importance_weights = list()

        #The lowest-priority experience defines the maximum importance sampling weight
        prob_min = self.min_tree.min() / self.sum_tree.sum()
        max_importance_weight = (prob_min * self.nb_entries)  ** (-beta)
        obs_t, act_t, rews, obs_t1, dones = [], [], [], [], []

        experiences = list()
        for idx in idxs:
            while idx < self.window_length + 1:
                idx += 1
            while idx + nstep > self.nb_entries and self.nb_entries < self.limit:
                # We are fine with nstep spilling back to the beginning of the buffer
                # once it has been filled.
                idx -= 1
            terminal0 = self.terminals[idx - 2]
            while terminal0:
                # Skip this transition because the environment was reset here. Select a new, random
                # transition and use this instead. This may cause the batch to contain the same
                # transition twice.
                idx = sample_batch_indexes(self.window_length + 1, self.nb_entries - nstep, size=1)[0]
                terminal0 = self.terminals[idx - 2]

            assert self.window_length + 1 <= idx < self.nb_entries

            #probability of sampling transition is the priority of the transition over the sum of all priorities
            prob_sample = self.sum_tree[idx] / self.sum_tree.sum()
            importance_weight = (prob_sample * self.nb_entries) ** (-beta)
            #normalize weights according to the maximum value
            importance_weights.append(importance_weight/max_importance_weight)

            #assemble the initial state from the ringbuffer.
            state0 = [self.observations[idx - 1]]
            for offset in range(0, self.window_length - 1):
                current_idx = idx - 2 - offset
                assert current_idx >= 1
                current_terminal = self.terminals[current_idx - 1]
                if current_terminal and not self.ignore_episode_boundaries:
                    # The previously handled observation was terminal, don't add the current one.
                    # Otherwise we would leak into a different episode.
                    break
                state0.insert(0, self.observations[current_idx])
            while len(state0) < self.window_length:
                state0.insert(0, zeroed_observation(state0[0]))

            action = self.actions[idx - 1]
            # N-step TD
            reward = 0
            nstep = nstep
            for i in range(nstep):
                reward += (gamma**i) * self.rewards[idx + i - 1]
                if self.terminals[idx + i - 1]:
                    #episode terminated before length of n-step rollout.
                    nstep = i
                    break

            terminal1 = self.terminals[idx + nstep - 1]

            # We assemble the second state in a similar way.
            state1 = [self.observations[idx + nstep - 1]]
            for offset in range(0, self.window_length - 1):
                current_idx = idx + nstep - 1 - offset
                assert current_idx >= 1
                current_terminal = self.terminals[current_idx - 1]
                if current_terminal and not self.ignore_episode_boundaries:
                    # The previously handled observation was terminal, don't add the current one.
                    # Otherwise we would leak into a different episode.
                    break
                state1.insert(0, self.observations[current_idx])
            while len(state1) < self.window_length:
                state1.insert(0, zeroed_observation(state0[0]))

            assert len(state0) == self.window_length
            assert len(state1) == len(state0)
            experiences.append(Experience(state0=state0, action=action, reward=reward,
                                          state1=state1, terminal1=terminal1))
        assert len(experiences) == batch_size

        return tuple(list(experiences)+ [importance_weights, idxs])

    def update_priorities(self, idxs, priorities):
        #adjust priorities based on new TD error
        for i, idx in enumerate(idxs):
            assert 0 <= idx < self.limit
            #expert transition priorities receive an extra boost
            if idx < self.permanent_idx:
                priority = (priorities[i] ** self.alpha) + .999
            else:
                priority = (priorities[i] ** self.alpha)
            self.sum_tree[idx] = priority
            self.min_tree[idx] = priority
            self.max_priority = max(self.max_priority, priority)

    def calculate_beta(self, current_step):
        a = float(self.end_beta - self.start_beta) / float(self.steps_annealed)
        b = float(self.start_beta)
        current_beta = min(self.end_beta, a * float(current_step) + b)
        return current_beta

    def get_config(self):
        config = super(PartitionedMemory, self).get_config()
        config['alpha'] = self.alpha
        config['start_beta'] = self.start_beta
        config['end_beta'] = self.end_beta
        config['beta_steps_annealed'] = self.steps_annealed
        config['pre_load_data'] = self.pre_load_data

    @property
    def nb_entries(self):
        """Return number of observations
        # Returns
            Number of observations
        """
        return len(self.observations)

max_stars_repo_name: rohanbadlani/keras-rl

max_issues_repo_name: rohanbadlani/keras-rl

max_forks_repo_name: rohanbadlani/keras-rl


text:
# coding=utf-8
import polyglot
from java.lang import System
import numpy as np
from random import random, randint, seed

from benchmark import Benchmark, time_phase, DEFAULT_BLOCK_SIZE_1D, DEFAULT_NUM_BLOCKS
from benchmark_result import BenchmarkResult

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

SQUARE_KERNEL = """
extern "C" __global__ void square(float* x, float* y, int n) {
    for(int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) {
        y[i] = x[i] * x[i];
    }
}
"""

DIFF_KERNEL = """
extern "C" __global__ void diff(const float* x, const float* y, float* z, int n) {
    for(int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) {
        z[i] = x[i] - y[i];
    }
}
"""

REDUCE_KERNEL = """

// From https://devblogs.nvidia.com/faster-parallel-reductions-kepler/

__inline__ __device__ float warp_reduce(float val) {
    int warp_size = 32;
    for (int offset = warp_size / 2; offset > 0; offset /= 2)
        val += __shfl_down_sync(0xFFFFFFFF, val, offset);
    return val;
}

__global__ void reduce(float *x, float *y, float* z, int N) {
    int warp_size = 32;
    float sum = float(0);
    for(int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) {
        sum += x[i] - y[i];
    }
    sum = warp_reduce(sum); // Obtain the sum of values in the current warp;
    if ((threadIdx.x & (warp_size - 1)) == 0) // Same as (threadIdx.x % warp_size) == 0 but faster
        atomicAdd(z, sum); // The first thread in the warp updates the output;
}
"""

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

class Benchmark1(Benchmark):
    """
    Compute the sum of difference of squares of 2 vectors, using multiple GrCUDA kernels.
    Structure of the computation:
       A: x^2 ──┐
                ├─> C: z=sum(x-y)
       B: x^2 ──┘
    """

    def __init__(self, benchmark: BenchmarkResult, nvprof_profile: bool = False):
        super().__init__("b1", benchmark, nvprof_profile)
        self.size = 0
        self.x = None
        self.y = None
        self.x1 = None
        self.y1 = None
        self.z = None
        self.res = None
        self.square_kernel = None
        self.diff_kernel = None
        self.reduce_kernel = None
        self.cpu_result = 0

        # self.num_blocks = DEFAULT_NUM_BLOCKS
        self.block_size = DEFAULT_BLOCK_SIZE_1D

    @time_phase("allocation")
    def alloc(self, size: int, block_size: dict = None) -> None:
        self.size = size
        self.block_size = block_size["block_size_1d"]

        # Allocate 2 vectors;
        self.x = polyglot.eval(language="grcuda", string=f"float[{size}]")
        self.y = polyglot.eval(language="grcuda", string=f"float[{size}]")
        self.x1 = polyglot.eval(language="grcuda", string=f"float[{size}]")
        self.y1 = polyglot.eval(language="grcuda", string=f"float[{size}]")

        # Allocate a support vector;
        self.res = polyglot.eval(language="grcuda", string=f"float[1]")

        # Build the kernels;
        build_kernel = polyglot.eval(language="grcuda", string="buildkernel")
        self.square_kernel = build_kernel(SQUARE_KERNEL, "square", "pointer, pointer, sint32")
        self.reduce_kernel = build_kernel(REDUCE_KERNEL, "reduce", "const pointer, const pointer, pointer, sint32")

    @time_phase("initialization")
    def init(self):
        self.random_seed = randint(0, 10000000)
        seed(self.random_seed)
        for i in range(self.size):
            if self.benchmark.random_init:
                self.x[i] = random()
                self.y[i] = 2 * random()
            else:
                self.x[i] = 1 / (i + 1)
                self.y[i] = 2 / (i + 1)

    @time_phase("reset_result")
    def reset_result(self) -> None:
        if self.benchmark.random_init:
            seed(self.random_seed)
            for i in range(self.size):
                self.x[i] = random()
                self.y[i] = 2 * random()
        else:
            for i in range(self.size):
                self.x[i] = 1 / (i + 1)
                self.y[i] = 2 / (i + 1)
        self.res[0] = 0.0

    def execute(self) -> object:
        self.block_size = self._block_size["block_size_1d"]
        start_comp = System.nanoTime()
        start = 0

        # A, B. Call the kernel. The 2 computations are independent, and can be done in parallel;
        self.execute_phase("square_1", self.square_kernel(self.num_blocks, self.block_size), self.x, self.x1, self.size)
        self.execute_phase("square_2", self.square_kernel(self.num_blocks, self.block_size), self.y, self.y1, self.size)

        # C. Compute the sum of the result;
        self.execute_phase("reduce", self.reduce_kernel(self.num_blocks, self.block_size), self.x1, self.y1, self.res, self.size)

        # Add a final sync step to measure the real computation time;
        if self.time_phases:
            start = System.nanoTime()
        result = self.res[0]
        end = System.nanoTime()
        if self.time_phases:
            self.benchmark.add_phase({"name": "sync", "time_sec": (end - start) / 1_000_000_000})
        self.benchmark.add_computation_time((end - start_comp) / 1_000_000_000)
        self.benchmark.add_to_benchmark("gpu_result", result)
        if self.benchmark.debug:
            BenchmarkResult.log_message(f"\tgpu result: {result:.4f}")

        return result

    def cpu_validation(self, gpu_result: object, reinit: bool) -> None:
        # Recompute the CPU result only if necessary;
        start = System.nanoTime()
        if self.current_iter == 0 or reinit:
            # Re-initialize the random number generator with the same seed as the GPU to generate the same values;
            seed(self.random_seed)
            if self.benchmark.random_init:
                x_g = np.zeros(self.size)
                y_g = np.zeros(self.size)
                for i in range(self.size):
                    x_g[i] = random()
                    y_g[i] = 2 * random()
            else:
                x_g = 1 / np.linspace(1, self.size, self.size)
                y_g = 2 / np.linspace(1, self.size, self.size)

            x_g = x_g ** 2
            y_g = y_g ** 2
            x_g -= y_g
            self.cpu_result = np.sum(x_g)
        cpu_time = System.nanoTime() - start
        difference = np.abs(self.cpu_result - gpu_result)
        self.benchmark.add_to_benchmark("cpu_time_sec", cpu_time)
        self.benchmark.add_to_benchmark("cpu_gpu_res_difference", difference)
        if self.benchmark.debug:
            BenchmarkResult.log_message(f"\tcpu result: {self.cpu_result:.4f}, " +
                                        f"difference: {difference:.4f}, time: {cpu_time:.4f} sec")


text:
import os
import gin
import logging
import random
import numpy as np
import torch
import torch.nn as nn
from pytorch_wrapper import DATA_DIR
import torch.nn.functional as F
from torch.utils.data import Dataset, DataLoader
from torch.utils.tensorboard import SummaryWriter
from tqdm import trange, tqdm

from pytorch_wrapper.utils import get_device, gin_wrap, set_seed
from pytorch_wrapper.data import DataModule
from pytorch_wrapper.models.nlp_models import NgramLM
from pytorch_wrapper.trainer import Trainer
from pytorch_wrapper.optimizer import OptimConfig

logger = logging.getLogger(__name__)

def load_tokens(path):
    with open(path, "r") as f:
        tokens_list = []
        for line in f:
            line = line.strip()
            if not line:
                continue
            tokens_list.append(["<bos>"] + line.split() + ["<eos>"])
    return tokens_list

def get_ngrams_arr(tokens_list, tok2idx, n, unknown_tok="<unk>"):
    ngram_rows = []
    for tokens in tokens_list:
        for i in range(len(tokens)-n+1):
            ngram = tokens[i:i+n]
            ngram_indices = [tok2idx.get(token, tok2idx[unknown_tok]) for token in ngram]
            ngram_rows.append(ngram_indices)
    return np.array(ngram_rows)

@gin.configurable
def get_datamodule(ngram_n, unknown_tok="<unk>"):
    train_tokens = load_tokens(os.path.join(DATA_DIR, "ptb.train.nounk.txt"))
    val_tokens = load_tokens(os.path.join(DATA_DIR, "ptb.valid.txt"))
    test_tokens = load_tokens(os.path.join(DATA_DIR, "ptb.valid.txt"))
    vocab = sorted(list({token for token_line in train_tokens for token in token_line})) + [unknown_tok]
    tok2idx = {vocab: i for i, vocab in enumerate(vocab)}
    idx2tok = dict(enumerate(vocab))
    train_arr = get_ngrams_arr(train_tokens, tok2idx, n=ngram_n)
    val_arr = get_ngrams_arr(val_tokens, tok2idx, n=ngram_n)
    test_arr = get_ngrams_arr(test_tokens, tok2idx, n=ngram_n)

    train_dataset = NgramDataset(train_arr)
    val_dataset = NgramDataset(val_arr)
    test_dataset = NgramDataset(test_arr)

    dm = DataModule(train_dataset=train_dataset, valid_dataset=val_dataset, test_dataset=test_dataset)
    return dm, vocab

class NgramDataset(Dataset):

    def __init__(self, tokenized_arr):
        self.tokenized_arr = tokenized_arr

    def __len__(self):
        return self.tokenized_arr.shape[0]

    def __getitem__(self, key):
        return self.tokenized_arr[key, :-1], self.tokenized_arr[key, -1]

@gin.configurable
def start_experiment(output_dir, max_epochs, seed=None):
    set_seed(seed=seed)
    writer = SummaryWriter(log_dir=os.path.join(output_dir, 'tb_logs'))
    device = get_device()
    dm, vocab = get_datamodule()

    model = NgramLM(vocab_size=len(vocab))
    optimizer = OptimConfig().create_optimizer(model)
    trainer = Trainer(optimizer, output_dir, max_epochs, tb_writer=writer, device=device)

    trainer.fit(model, dm)

if __name__ == "__main__":

    gin_wrap(start_experiment)


text:
from collections import defaultdict
import h5py
import numpy as np
import pandas as pd
from itertools import product

systems = ['Li2', 'Be2', 'B2', 'C2']
dets = [1, 3, 10, 30, 100]

def get_mean_err(energies):
    return energies.mean(), energies.mean(axis=0).std() / np.sqrt(energies.shape[1])

results = defaultdict(list)
with h5py.File(f'../data/raw/data_pub_diatomics.h5', 'r') as f:
    for system, d in product(systems, dets):
        results[system, d] = (lambda x: pd.Series({'energy': x[0], 'err': x[1],}))(
            get_mean_err(f[system][f'{d}det']['evaluate'][...][:, :, 0])
        )

results = (
    pd.concat(results, names=['system', 'ndet']).unstack().sort_index().reset_index()
)
results.to_csv('../data/final/diatomics.csv', index=False)

max_stars_repo_name: noegroup/dlqmc-project

max_issues_repo_name: noegroup/dlqmc-project

max_forks_repo_name: noegroup/dlqmc-project


text:
# Synthia (https://github.com/dmey/synthia).
# Copyright (c) 2020 D. Meyer and T. Nagler. Licensed under the MIT License.

import time
import numpy as np

import synthia as syn

def main():
    import pyvinecopulib as pv

    n_samples = 100
    n_features = 50
    n_synthetic_samples = 500
    input_data = np.random.normal(size=(n_samples, n_features))

    generator = syn.CopulaDataGenerator()

    t0 = time.time()
    generator.fit(input_data, syn.GaussianCopula())
    generator.generate(n_synthetic_samples)
    print(f'synthia: {time.time() - t0}s')

    t0 = time.time()
    generator.fit(input_data, syn.VineCopula(pv.FitControlsVinecop(
        family_set=[pv.gaussian], trunc_lvl=1, select_trunc_lvl=False)))
    generator.generate(n_synthetic_samples)
    print(f'pyvinecopulib: {time.time() - t0}s')

if __name__ == '__main__':
    main()

max_stars_repo_name: khinsen/synthia

max_issues_repo_name: khinsen/synthia

max_forks_repo_name: khinsen/synthia


text:
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from typing import Optional, Tuple
from tensorflow.contrib import rnn

def RNN(x, weights, biases, n_input, n_steps, n_hidden):
    # 准备数据形状以匹配RNN函数的要求
    # 当前数据输入形状: (batch_size, n_steps, n_input)
    # 要求形状: 'n_steps' 张量形状列表 (batch_size, n_input)
    # 置换 batch_size and n_steps
    x = tf.transpose(x, [1, 0, 2])
    # 重定义形状 (n_steps*batch_size, n_input)
    x = tf.reshape(x, [-1, n_input])
    # 分割得到列表 'n_steps' 形状张量 (batch_size, n_input)
    x = tf.split(x, n_steps, axis=0)
    # 用张量流定义一个LSTM单元
    lstm_cell = rnn.BasicLSTMCell(n_hidden, forget_bias=1.0)
    # 获取LSTM信元输出
    outputs, states = rnn.static_rnn(lstm_cell, x, dtype=tf.float32)
    # 线性激活函数，使用RNN内循环最后输出
    return tf.nn.bias_add(tf.matmul(outputs[-1], weights['out']), biases['out'])

def generate_sample(f: Optional[float] = 1.0, t0: Optional[float] = None, batch_size: int = 1,
                    predict: int = 50, samples: int = 100) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    生成数据样本。
    :参数 f: 用于所有时间序列或没有随机化的频率。
    :参数 t0: 用于所有时间序列或没有随机化的时间偏移。
    :参数 batch_size: 生成的时间序列的数量。
    :参数 predict: 要生成的未来样本的数量。
    :参数 samples: 要生成的过去（和当前）样本的数量。
    :返回: 元组包含过去的时间和值以及未来的时间和值。在所有输出中，每行表示批的一个时间序列。
    """
    Fs = 100

max_stars_repo_name: Zhousantu/adaptive-memory-curve

max_issues_repo_name: Zhousantu/adaptive-memory-curve

max_forks_repo_name: phil-zyx/adaptive-memory-curve


text:
import numpy as np
from src.features.build_features import Scaler
from scipy.spatial.distance import cdist
from scipy.stats import mode
import pandas as pd

class KNN():
    def __init__(self, n, metric='euclidean'):
        self.n = n
        self.xs_ = None
        self.y_ = None
        self.scaler_ = None
        self.metric = metric

    def fit(self, x, y):
        self.scaler_ = Scaler()
        self.scaler_.fit(x)
        self.xs_ = self.scaler_.transform(x)
        self.y_ = y
        return self

    def predict(self, x):
        # verify training has happened
        if self.xs_ is None:
            raise ValueError("You must train the model first")

        # scale x
        xs = self.scaler_.transform(x)

        # calculate the distance from every training vector to every x
        # result will be n_training by n_predict columns
        dists = cdist(self.xs_, xs, self.metric)

        # sort those to find out the closest to each x_predict
        inds = np.argsort(dists, axis=0)[:self.n, :]
        inds = inds.T

        # match those indices to column predictions
        l = mode(self.y_[inds], axis=1)

        labels = l.mode[:,0]

        return labels

    def correct(self, x, y):
        yp = self.predict(x).ravel()
        correct = np.argwhere(yp == y)
        return correct.shape[0]

    def confusion(self, x, y):
        yp = self.predict(x)

        tn = np.count_nonzero(np.where((y==2) & (yp==2)))
        fn = np.count_nonzero(np.where((y==4) & (yp==2)))
        fp = np.count_nonzero(np.where((y==2) & (yp==4)))
        tp = np.count_nonzero(np.where((y==4) & (yp==4)))

        r = np.array([
            [tn, fp],
            [fn, tp]
        ])

        return r

    def accuracy(self, x, y):
        r = self.confusion(x, y)
        return (r[0, 0] + r[1, 1]) / np.sum(r)

    def tpr(self, x, y):
        r = self.confusion(x, y)
        return r[1, 1] / (r[1, 1] + r[1, 0])

    def ppv(self, x, y):
        r = self.confusion(x, y)
        return r[1, 1] / (r[1, 1] + r[0, 1])

    def tnr(self, x, y):
        r = self.confusion(x, y)
        return r[0, 0] / (r[0, 0] + r[0, 1])

    def f1_score(self, x, y):
        num = 2 * self.ppv(x, y) * self.tpr(x, y)
        den = self.ppv(x, y) + self.tpr(x, y)
        return num / den

class Node():
    def __init__(self, impurity='entropy', min_impurity=0.1,
                 depth=0, max_depth=10):
        """Classifier with decision trees

        Implemented here with `pandas.DataFrame`s

        Paramters
        ---------
        impurity : {'entropy', 'gini', 'misclassification'}
        """
        impurity_functions = {
            'entropy': self.entropy,
            'gini': self.gini,
            'misclassification': self.misclassification
        }
        self.i_name = impurity
        self.impurity = impurity_functions[impurity]
        self.attribute_ = None
        self.value_ = None
        self.min_impurity = min_impurity
        self.depth = depth
        self.max_depth = max_depth
        self.left_ = None
        self.right_ = None
        self.is_leaf_ = False
        self.label_ = None
        self.confidence_ = None
        self.training_impurity_ = None

    def entropy(self, y):
        """Entropy of a class

        Satisfies equation 9.3 and 9.4 from the text:
        .. math::
            9.3: I = - sum_{i=1}^K p_m^i \log_2 p_m^i \\
            9.4: \phi(p, 1-p) = -p \log_2 p - (1-p) \log_2 (1-p)

        Parameters
        ----------
        y : pandas.DataFrame
            Training data encoded as one-hot vector

        Returns
        -------
        float
            Entropy of the training data
        """
        if y.shape[0] == 0:
            p = 0.0
        else:
            p = np.sum(y) / y.shape[0]

        if p == 0.0:
            fp = 0.0
        else:
            fp = -p * np.log(p) / np.log(2)

        if p == 1.0:
            sp = 0.0
        elif p > 1.0:
            raise Exception(f"p > 1.0: {p}")
        else:
            sp = (1 - p) * np.log(1 - p) / np.log(2)
        return fp - sp

    def gini(self, y):
        """Gini index for a class

        Satisfies equation 9.5 from the text:
        .. math::
            9.5: \phi(p, 1-p) = 2 p (1-p)

        Parameters
        ----------

        Returns
        -------

        """
        if y.shape[0] == 0:
            p = 0.0
        else:
            p = np.sum(y) / y.shape[0]
        return 2 * p * (1 - p)

    def misclassification(self, y):
        """Misclassification score of a class

        Implements equation 9.6 from the text:
        .. math::
            9.6: \phi(p, 1-p) = 1 - \mathrm{max}(p, 1-p)

        Parameters
        ----------

        Returns
        -------

        """
        if y.shape[0] == 0:
            p = 0.0
        else:
            p = np.sum(y) / y.shape[0]
        return 1 - np.max([p, 1-p])

    def fit(self, data):
        column_names = data.columns
        x = data[column_names[:-1]]
        y = data[column_names[-1]]
        # Calculate impurity
        impurity = self.impurity(y)
        self.training_impurity_ = impurity
        # Decide if I should split
        if impurity > self.min_impurity and self.depth < self.max_depth:
            # Calculate the split
            self.attribute_, self.value_ = self.split_attribute(x, y, impurity=self.impurity)
            self.left_ = Node(impurity=self.i_name, min_impurity=self.min_impurity,
                              depth=self.depth+1, max_depth=self.max_depth)
            self.right_ = Node(impurity=self.i_name, min_impurity=self.min_impurity,
                               depth=self.depth+1, max_depth=self.max_depth)
            left_data = data[data[self.attribute_] <= self.value_]
            right_data = data[data[self.attribute_] > self.value_]
            self.left_.fit(left_data)
            self.right_.fit(right_data)
        else:
            self.is_leaf_ = True
            self.label_ = mode(y)
            self.confidence_ = self.impurity(y)

    def predict(self, x):
        if self.is_leaf_:
            return self.label_.mode
        else:
            # Determine left or right
            if x[self.attribute_] <= self.value_:
                a = self.left_.predict(x)
                return a
            else:
                a = self.right_.predict(x)
                return a

    def split_attribute(self, x, y, impurity):
        min_entropy = np.finfo(np.float64).max # 1e308 on test system
        xt = x.copy()
        xt['Labels'] = y
        split_column = None
        split_value = None
        for c in x.columns:
            for v in x[c].unique():
                left_split = xt[xt[c] <= v]
                right_split = xt[xt[c] > v]
                left_impurity = impurity(left_split['Labels'])
                right_impurity = impurity(right_split['Labels'])
                e = left_impurity + right_impurity
                if e < min_entropy:
                    min_entropy = e
                    split_column = c
                    split_value = v + 0.5
        return split_column, split_value

    @np.vectorize
    def true_negative(self, pred, true):
        if pred == 0 and true == 0:
            return True
        else:
            return False

    @np.vectorize
    def false_negative(self, pred, true):
        if pred == 0 and true == 1:
            return True
        else:
            return False

    @np.vectorize
    def true_positive(self, pred, true):
        if pred == 1 and true == 1:
            return True
        else:
            return False

    @np.vectorize
    def false_positive(self, pred, true):
        if pred == 1 and true == 0:
            return True
        else:
            return False

    def confusion(self, pred, true):
        tp = np.sum(self.true_negative(pred, true))
        fp = np.sum(self.false_positive(pred, true))
        tn = np.sum(self.true_negative(pred, true))
        fn = np.sum(self.false_negative(pred, true))

        r = np.array([
            [tp, fp],
            [fn, tn]
        ])
        d = pd.DataFrame(r, index=[1, 0], columns=[1, 0])
        return d

    def accuracy(self, pred, true):
        tp = np.sum(self.true_negative(pred, true))
        fp = np.sum(self.false_positive(pred, true))
        tn = np.sum(self.true_negative(pred, true))
        fn = np.sum(self.false_negative(pred, true))
        return (tn + tp) / (tn + tp + fn + fp)

    def true_positive_rate(self, pred, true):
        """True positive rate

        Also known as the recall or sensitivity
        """
        tp = np.sum(self.true_negative(pred, true))
        fn = np.sum(self.false_negative(pred, true))

        return tp / (tp + fn)

    def positive_predictive_value(self, pred, true):
        """Positive predictive value

        Also known as precision
        """
        tp = np.sum(self.true_negative(pred, true))
        fp = np.sum(self.false_positive(pred, true))

        return tp / (tp + fp)

    def true_negative_rate(self, pred, true):
        """True negative rate

        Also known as specificity
        """
        fp = np.sum(self.false_positive(pred, true))
        tn = np.sum(self.true_negative(pred, true))

        return tn / (tn + fp)

    def f1(self, pred, true):
        """F1 score
        """
        ppv = self.positive_predictive_value(pred, true)
        tpr = self.true_positive_rate(pred, true)

        return 2 * ppv * tpr / (ppv + tpr)

    def html_print(self):
        s = "

"
        if self.is_leaf_:
            s += ""
            s += self.__repr__()
            s += "	"
        else:
            s += ""
            s += self.__repr__()
            s += "	"
            s += "

"
            s += ""
            s += self.left_.html_print()
            s += "	"
            s += self.right_.html_print()
            s += "	"
        s += "

"
        return s

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Nov  4 11:25:03 2019

@author: lealp
"""

import pandas as pd
pd.set_option('display.width', 50000)
pd.set_option('display.max_rows', 50000)
pd.set_option('display.max_columns', 5000)

import numpy as np
import xarray as xr

from dask.diagnostics import ProgressBar

import matplotlib.pyplot as plt

def binning_function(x, time,  distribution_type='Positive', b=False):

    y = np.where(np.abs(x)==np.inf, 0, x)

    y = np.where(np.isnan(y), 0, y)

    if np.all(y) == 0:

        return np.where(y==0, np.nan, np.nan)

    else:

        Classified = pd.qcut(y, np.linspace(0.0, 1, 11))

        codes = Classified.codes + 1

        return codes

def xarray_parse_extremes(ds, dim=['time'],
                          dask='allowed',
                          new_var_name='classes',
                          kwargs={'b': False,
                                  'distribution_type':'Positive'}):

    kwargs['time'] = np.unique(ds['time'].values)

    filtered = xr.apply_ufunc(binning_function,
                                  ds,
                                  dask=dask,
                                  vectorize=True,
                                  input_core_dims=[dim],
                                  #exclude_dims = [dim],
                                  output_core_dims=[dim],
                                  kwargs=kwargs,
                                  output_dtypes=[int],

                                  )

    filtered.name = new_var_name

    return filtered

if '__main__' == __name__:

    ds = xr.tutorial.open_dataset('rasm').load()

    def parse_datetime(time):
        return pd.to_datetime([str(x) for x in time])

    ds.coords['time'] = parse_datetime(ds.coords['time'].values)

    ds.isel({'time':0}).Tair.plot()

    #ds = ds.sel({'x':slice(154, 159), 'y':slice(32, 40) })

    with ProgressBar():
        da_binned = xarray_parse_extremes(ds['Tair'].fillna(0),
                                          ['time'],
                                         dask='parallelized')

    binned = da_binned.copy(deep=True)

    binned_ds = binned.to_dataset()

    ########### Plotting

    plt.figure()

    ds['Tair'].isel({'time':3}).plot(cmap='viridis', ax=plt.subplot(1,2,1))
    binned.isel({'time':3}).plot(cmap='viridis', ax=plt.subplot(1,2,2))

    import matplotlib.pyplot as plt
    import matplotlib as mpl

    cmap = mpl.cm.viridis

    fig3, axs = plt.subplots(3,3, sharex=True, sharey=True)
    axs = axs.ravel()

    for i, ax in enumerate(axs):
        print(i)
        temp = binned.where(binned==i+1, 1, 0).sum('time')

        temp.plot(ax=ax, add_colorbar=False)

        bounds = np.arange(0,temp.max())
        norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
        cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
                                    norm=norm,
                                    boundaries=bounds,
                                    extend='both',
                                    ticks=bounds,
                                    spacing='proportional',
                                    orientation='vertical')

    fig3.show()


text:
from os.path import exists, join, basename
import fasttext
import pickle
from nir.logger import log
import gc
import numpy as np
from gensim.models import KeyedVectors

class EmbeddingBase:
    def __init__(self, cache_folder, prefix_name, tokenizer, path):
        self.cache_folder = cache_folder
        self.prefix_name = prefix_name
        self.path = path
        self.tokenizer = tokenizer
        self.vocab_size = None
        self.embedding_size = None
        self.matrix = None
        self.name = Word2Vec.get_name(self.path, self.tokenizer.name)

    def _normalize_vector(self, x):
        return x/np.linalg.norm(x)

    def has_matrix(self):
        return self.matrix is not None

    def build_matrix(self, index_zero_reserved = True):
        raise NotImplementedError()

    def embedding_matrix(self):
        if not self.has_matrix():
            if not self.tokenizer.is_trained():
                raise "At this point the tokenizer should already have been trained"
            else:
                # build matrix and save
                self.build_matrix()

        return self.matrix

class Word2Vec(EmbeddingBase):
    def __init__(self, **kwargs):
        super(Word2Vec, self).__init__(**kwargs)
        self.name = Word2Vec.get_name(self.path, self.tokenizer.name)

    @staticmethod
    def get_name(file_path, tokenizer_name):
        return "WORD2VEC_embedding_{}_{}".format(basename(file_path).split(".")[0],
                                        tokenizer_name)

    @staticmethod
    def maybe_load(cache_folder, prefix_name, tokenizer, path):
        name = Word2Vec.get_name(path, tokenizer.name)
        cache_path = join(cache_folder, name)
        if exists(cache_path):
            print("[LOAD FROM CACHE] Load embedding matrix from", cache_path)
            with open(cache_path, "rb") as f:
                data = pickle.load(f)
                ft = Word2Vec(cache_folder=cache_folder, prefix_name=prefix_name, tokenizer=tokenizer, path=path)
                ft.matrix = data["matrix"]
                ft.vocab_size = data["vocab_size"]
                ft.embedding_size = data["embedding_size"]
                return ft
        else:
            print("Model not found, new Word2Vec instance was returned")
            return Word2Vec(cache_folder=cache_folder, prefix_name=prefix_name, tokenizer=tokenizer, path=path)

    def build_matrix(self, index_zero_reserved = True):
        print("[WORD2VEC] Creating embedding matrix")
        trained_ft_model = None
        try:
            log.info("[WORD2VEC] load model")
            trained_ft_model = KeyedVectors.load(self.path, mmap='r')
            log.info("[WORD2VEC] build embedding matrix")

            self.vocab_size = self.tokenizer.vocabulary_size() + (1 if index_zero_reserved else 0)
            self.embedding_size = trained_ft_model.vector_size

            self.matrix = np.random.normal(size=(self.vocab_size, self.embedding_size))

            if index_zero_reserved:
                self.matrix[0] = self._normalize_vector(self.matrix[0])

            for word in self.tokenizer.get_vocabulary():
                self.matrix[self.tokenizer.word_index[word]] = self._normalize_vector(trained_ft_model[word])

        except Exception as e:
            log.error(e)
            raise e
        finally:
            del trained_ft_model
            log.info("GC ({})".format(gc.collect()))

        # save after gc
        cache_path = join(self.cache_folder, self.name)
        with open(cache_path, "wb") as f:
            log.info("[FASTTEXT] save")
            data = {"vocab_size": self.vocab_size,
                    "embedding_size": self.embedding_size,
                    "matrix": self.matrix}
            pickle.dump(data, f, protocol=4)

class FastText(EmbeddingBase):
    def __init__(self, **kwargs):
        super(FastText, self).__init__(**kwargs)
        self.name = FastText.get_name(self.path, self.tokenizer.name)

    @staticmethod
    def get_name(file_path, tokenizer_name):
        return "embedding_{}_{}".format(basename(file_path).split(".")[0],
                                        tokenizer_name)

    @staticmethod
    def maybe_load(cache_folder, prefix_name, tokenizer, path):
        name = FastText.get_name(path, tokenizer.name)
        cache_path = join(cache_folder, name)
        if exists(cache_path):
            print("[LOAD FROM CACHE] Load embedding matrix from", cache_path)
            with open(cache_path, "rb") as f:
                data = pickle.load(f)
                ft = FastText(cache_folder=cache_folder, prefix_name=prefix_name, tokenizer=tokenizer, path=path)
                ft.matrix = data["matrix"]
                ft.vocab_size = data["vocab_size"]
                ft.embedding_size = data["embedding_size"]
                return ft
        else:
            print("Model not found, new FastText instance was returned")
            return FastText(cache_folder=cache_folder, prefix_name=prefix_name, tokenizer=tokenizer, path=path)

    def build_matrix(self, index_zero_reserved = True):
        print("[FASTTEXT] Creating embedding matrix")
        trained_ft_model = None
        try:
            log.info("[FASTTEXT] load model")
            trained_ft_model = fasttext.load_model(self.path)
            log.info("[FASTTEXT] build embedding matrix")

            self.vocab_size = self.tokenizer.vocabulary_size() + (1 if index_zero_reserved else 0)
            self.embedding_size = trained_ft_model.get_dimension()

            self.matrix = np.random.normal(size=(self.vocab_size, self.embedding_size))

            if index_zero_reserved:
                self.matrix[0] = self.matrix[0]/np.linalg.norm(self.matrix[0])

            for word in self.tokenizer.get_vocabulary():
                self.matrix[self.tokenizer.word_index[word]] = trained_ft_model[word]

            # normalize all entries, NOTE that previous iteration may miss some entries (out-of-vocabulary)
            for i in range(self.matrix.shape[0]):
                self.matrix[i] = self._normalize_vector(self.matrix[i])

        except Exception as e:
            log.error(e)
            raise e
        finally:
            del trained_ft_model
            log.info("GC ({})".format(gc.collect()))

        # save after gc
        cache_path = join(self.cache_folder, self.name)
        with open(cache_path, "wb") as f:
            log.info("[FASTTEXT] save")
            data = {"vocab_size": self.vocab_size,
                    "embedding_size": self.embedding_size,
                    "matrix": self.matrix}
            pickle.dump(data, f, protocol=4)

max_stars_repo_name: T-Almeida/DL-IR-tools

max_issues_repo_name: T-Almeida/DL-IR-tools

max_forks_repo_name: T-Almeida/DL-IR-tools


text:
import scipy, scipy.sparse,scipy.sparse.linalg, numpy as np

def computeAnatomicalDirectionsHeatEquation(mesh, apexPointId, valvesPointsId):
    """
    Computes the longitudinal direction, using a single orifice.
    """
    triangles = mesh.faces.reshape((-1, 4))[:, 1:]
    boundary = {v: 1. for v in valvesPointsId}
    boundary[apexPointId] = 0
    heat = solveLaplaceBeltrami(mesh.points, triangles, boundary)
    vLongitudinal = grad_3d(mesh, heat)
    vLongitudinal = vLongitudinal / np.linalg.norm(vLongitudinal, axis = 1).reshape((-1, 1))
    vCircumferential = np.cross(vLongitudinal, mesh.cell_normals)
    return vLongitudinal, vCircumferential

def solveLaplaceBeltrami(points, triangles, boundary):
    """
    Solve laplace-beltrami equation with boundary conditions
    """
    nPoints = len(points)
    nConstraints = len(boundary)
    L = ComputeCotangentLaplacian( points.T, triangles.T)
    lagrange_multipliers = np.zeros((len(boundary), nPoints))
    boundary_values = np.zeros(nConstraints)
    for i, (k, v) in enumerate(boundary.items()):
        lagrange_multipliers[i, k] = 1
        boundary_values[i] = v
    zeros = np.zeros((nConstraints, nConstraints))
    L_hat = scipy.sparse.bmat([[L, lagrange_multipliers.T], [lagrange_multipliers, zeros]])
    b_hat = np.concatenate([np.zeros(nPoints), boundary_values])
    x_hat = scipy.sparse.linalg.spsolve(L_hat, b_hat)
    return x_hat[:nPoints]

def ComputeCotangentLaplacian( vertex, faces ):
    """
    Calculates the laplacian of a mesh
    vertex 3xN numpy.array: vertices
    faces 3xM numpy.array: faces
    Return the Laplacian (sparse matrix probably).
    """
    n = vertex.shape[1]
    m = faces.shape[1]

    #compute mesh weight matrix
    W = scipy.sparse.coo_matrix((n,n))
    for i in np.arange(1,4,1):
        i1 = np.mod(i-1,3)
        i2 = np.mod(i  ,3)
        i3 = np.mod(i+1,3)
        pp = vertex[:,faces[i2,:]] - vertex[:,faces[i1,:]]
        qq = vertex[:,faces[i3,:]] - vertex[:,faces[i1,:]]
        #% normalize the vectors
        pp = pp / np.sqrt(np.sum(pp**2, axis=0))
        qq = qq / np.sqrt(np.sum(qq**2, axis=0))

        #% compute angles
        ang = np.arccos(np.sum(pp*qq, axis=0))
        W = W + scipy.sparse.coo_matrix( (1 / np.tan(ang),(faces[i2,:],faces[i3,:])), shape=(n, n) )
        W = W + scipy.sparse.coo_matrix( (1 / np.tan(ang),(faces[i3,:],faces[i2,:])), shape=(n, n) )

    #compute laplacian
    d = W.sum(axis=0)
    D = scipy.sparse.dia_matrix((d, 0), shape=(n,n) )
    L = D - W

    return L

def propagateDirections(meshRef, directions, meshTarget):
    pass

# Utils
# Implemented by Thomas D'Argent
def grad_3d(mesh, scalarField):
    cell_gradients = np.zeros((mesh.n_cells, 3))
    triangles = mesh.faces.reshape((-1, 4))[:, 1:]

    for i, (p0_i, p1_i, p2_i) in enumerate(triangles):
        p10 = mesh.points[p1_i] - mesh.points[p0_i]
        p20 = mesh.points[p2_i] - mesh.points[p0_i]

        gradient = np.zeros(3)
        # compute the difference in value between p{1,2} and p0
        t10 = scalarField[p1_i] - scalarField[p0_i]
        t20 = scalarField[p2_i] - scalarField[p0_i]

        # t10 = gradient@p10 (resp. t20 and p20)
        gradient = np.linalg.pinv([p10, p20, np.cross(p10, p20)])@np.array([t10, t20, 0])
        assert np.abs(gradient@p10 - t10) <= 1e-4
        assert np.abs(gradient@p20 - t20) <= 1e-4
        assert np.abs(gradient@np.cross(p10, p20)) <= 1e-4

        cell_gradients[i] = gradient
    return cell_gradients

max_stars_repo_name: gbernardino/rvmep

max_issues_repo_name: gbernardino/rvmep

max_forks_repo_name: gbernardino/rvmep

reason: import scipy


text:
import numpy as np
from ... import draw
from .Polygons import Polygons

class Arrows2D(draw.Arrows2D, Polygons):
    __doc__ = draw.Arrows2D.__doc__

    def render(self, *args, **kwargs):
        Polygons.render(self, *args, **kwargs)

max_stars_repo_name: hillarypan/plato

max_issues_repo_name: hillarypan/plato

max_forks_repo_name: hillarypan/plato