text:
#!/usr/bin/env python
"""
Configure folder for testing Cramer-Rao bounds calculations. Ideally
this should be a unit test, but calculating splines takes too long and
they are a bit big to include in a project, so such is life.

Hazen 10/17
"""
import argparse
import numpy
import os

import storm_analysis.sa_library.parameters as parameters
import storm_analysis.sa_library.sa_h5py as saH5Py

import storm_analysis.simulator.background as background
import storm_analysis.simulator.camera as camera
import storm_analysis.simulator.drift as drift
import storm_analysis.simulator.emitters_on_grid as emittersOnGrid
import storm_analysis.simulator.emitters_uniform_random as emittersUniformRandom
import storm_analysis.simulator.photophysics as photophysics
import storm_analysis.simulator.psf as psf
import storm_analysis.simulator.simulate as simulate

import storm_analysis.pupilfn.make_pupil_fn as makePupilFn

import storm_analysis.psf_fft.make_psf_from_pf as makePSFFromPF

import storm_analysis.spliner.measure_psf_beads as measurePSFBeads
import storm_analysis.spliner.psf_to_spline as psfToSpline

import storm_analysis.diagnostics.cramer_rao.settings as settings

#parser = argparse.ArgumentParser(description = 'Cramer-Rao diagnostics configuration.')
#args = parser.parse_args()

def configure():

    # Create PF for pupil function.
    #
    print("Creating pupil function.")
    pf_size = 2*(settings.spline_size - 1)
    makePupilFn.makePupilFunction("pupilfn.pfn",
                                  pf_size,
                                  settings.pixel_size * 1.0e-3,
                                  settings.zmn,
                                  z_offset = settings.z_offset)

    # Create PSF using pupil functions directly.
    #
    if False:
        print("Creating (theoritical) psf.")
        makePSFFromPF.makePSF("psf_fft.psf",
                              settings.spline_size,
                              settings.pixel_size * 1.0e-3,
                              settings.zmn,
                              settings.psf_fft_z_range,
                              settings.psf_fft_z_step)

        exit()

    # Localizations on a sparse parse grid for PSF
    # measurement for Spliner and PSF FFT.
    #
    print("Creating data for PSF measurement.")
    emittersOnGrid.emittersOnGrid("sparse_list.hdf5",
                                  6,
                                  3,
                                  1.5,
                                  40,
                                  0.0,
                                  settings.z_offset)

    # Create beads.txt file for spline measurement.
    #
    with saH5Py.SAH5Py("sparse_list.hdf5") as h5:
        locs = h5.getLocalizations()
        numpy.savetxt("beads.txt", numpy.transpose(numpy.vstack((locs['x'], locs['y']))))

    # Create drift file, this is used to displace the localizations in the
    # PSF measurement movie.
    #
    dz = numpy.arange(-settings.spline_z_range, settings.spline_z_range + 0.001, 0.01)
    drift_data = numpy.zeros((dz.size, 3))
    drift_data[:,2] = dz
    numpy.savetxt("drift.txt", drift_data)

    # Also create the z-offset file.
    #
    z_offset = numpy.ones((dz.size, 2))
    z_offset[:,1] = dz
    numpy.savetxt("z_offset.txt", z_offset)

    # Create simulated data for PSF measurement.
    #
    bg_f = lambda s, x, y, i3 : background.UniformBackground(s, x, y, i3, photons = 10)
    cam_f = lambda s, x, y, i3 : camera.Ideal(s, x, y, i3, 100.0)
    drift_f = lambda s, x, y, i3 : drift.DriftFromFile(s, x, y, i3, "drift.txt")
    pp_f = lambda s, x, y, i3 : photophysics.AlwaysOn(s, x, y, i3, 20000.0)
    psf_f = lambda s, x, y, i3 : psf.PupilFunction(s, x, y, i3, settings.pixel_size, settings.zmn, pf_size = pf_size)

    sim = simulate.Simulate(background_factory = bg_f,
                            camera_factory = cam_f,
                            drift_factory = drift_f,
                            photophysics_factory = pp_f,
                            psf_factory = psf_f,
                            x_size = settings.x_size,
                            y_size = settings.y_size)

    sim.simulate("psf.dax", "sparse_list.hdf5", dz.size)

    # Create spline for Spliner
    #

    # Measure the PSF for Spliner
    #
    print("Measuring PSF.")
    psf_name = "psf_spliner.psf"
    measurePSFBeads.measurePSFBeads("psf.dax",
                                    "z_offset.txt",
                                    "beads.txt",
                                    psf_name,
                                    aoi_size = int(settings.spline_size + 1),
                                    pixel_size = settings.pixel_size * 1.0e-3)

    # Measure the Spline.
    #

    # This is slow, sometimes you don't want to do it.
    if True:
        print("Measuring Spline.")
        psfToSpline.psfToSpline(psf_name, "psf.spline", settings.spline_size)

    # Create measured PSF for PSF FFT.
    #

    # Measure the PSF using spliner/measure_psf_beads.py
    #
    print("Measuring PSF.")
    measurePSFBeads.measurePSFBeads("psf.dax",
                                    "z_offset.txt",
                                    "beads.txt",
                                    "psf_fft.psf",
                                    aoi_size = int(settings.spline_size - 1),
                                    pixel_size = settings.pixel_size * 1.0e-3,
                                    z_range = settings.psf_fft_z_range,
                                    z_step = settings.psf_fft_z_step)

if (__name__ == "__main__"):
    configure()

max_stars_repo_name: bintulab/storm-analysis

max_issues_repo_name: bintulab/storm-analysis

max_forks_repo_name: bintulab/storm-analysis


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Created: 11-2020 - Carmelo Mordini <carmelo> <carmelo.mordini@unitn.it>

"""
Module docstring

"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.interpolate import RectBivariateSpline
from scipy.optimize import brentq

from becpy.data import hf_scan
from becpy.physics import trapped_hartree_fock as thf

mu0, T, N0, Nt, freq_x, freq_y = [hf_scan[k] for k in ['mu0', 'T', 'N0', 'Nt', 'freq_x', 'freq_y']]

omega_x = 2 * np.pi * freq_x
omega_rho = 2 * np.pi * freq_y
omega_ho = (omega_x * omega_rho**2)**(1 / 3)

N = N0 + Nt
bec_fraction = N0 / N

mu0_1 = mu0[:, 0]
T_1 = T[0]

mu0_lims = mu0_1.min(), mu0_1.max()

print('interpolating...')
interp_kwargs = dict(kx=1, ky=1, s=0)
_N = RectBivariateSpline(mu0_1, T_1, N, **interp_kwargs)
_N0 = RectBivariateSpline(mu0_1, T_1, N0, **interp_kwargs)

def __mu0(N, T):
    def fun(mu0, N, T):
        return _N(mu0, T*1e9, grid=False) - N
    return brentq(fun, *mu0_lims, args=(N, T))

print('done')
_mu0 = np.vectorize(__mu0)

def _lda_bec_fraction(N, T):
    return _N0(__mu0(N, T), T*1e9, grid=False) / N

lda_bec_fraction = np.vectorize(_lda_bec_fraction)

fig, (ax, ax2) = plt.subplots(1, 2, figsize=(12, 4))
c = ax.contourf(T, mu0, bec_fraction, 200)
ax.contour(T, mu0, bec_fraction, [0.05], colors='r')
plt.colorbar(c, ax=ax)
ax.set(xlabel='T [nK]', ylabel='mu0 [nK]')

c = ax2.contourf(T, mu0, N, levels=np.linspace(1e5, N.max(), 200))
plt.colorbar(c, ax=ax2)

fig, (ax, ax1) = plt.subplots(1, 2, figsize=(12, 4))
t = np.linspace(0.5, 1.3, 200)

f2, a2 = plt.subplots()
a2.hist(N.ravel(), 200)

for N in 1e6, 5e6, 1e7, 3e7:
    Tc = thf.critical_temperature(N, omega_ho)
    T = t * Tc

    eta = thf.eta(N, omega_ho)
    u0 = eta * Tc * 1e-9
    f0 = thf.bec_fraction(N, T, omega_ho, clip=False)

    u1 = _mu0(N, T)
    f1 = lda_bec_fraction(N, T)

    l, = ax.plot(t, f1, label=f"{N*1e-6:.1f} M")
    print(l)
    ax.plot(t, f0, '--', color=l.get_color())

    ax1.plot(T, u1)

ax.legend()
plt.show()

ext: py

max_stars_repo_name: BEC-Trento/becpy

max_issues_repo_name: BEC-Trento/becpy

max_forks_repo_name: BEC-Trento/becpy

reason: import numpy,from scipy


text:
import torch.utils.data as data
import PIL.Image as Image
import os
import json
import numpy as np
import torch

def make_dataset(root, seed, name):
    imgs = []
    img_labels = {}
    img_names = os.listdir(root +'/'+ name+'-image/')
    img_names = sorted(img_names, key=lambda i: int(i.split(".")[0]))

    for i in seed:
        img_name = img_names[i]
        img = os.path.join(root +'/'+ name+ '-image/', img_name)
        mask = os.path.join(root +'/'+ name+ '-mask/', img_name)
        imgs.append((img, mask))
    return imgs

def make_testset(root):
    imgs = []
    img_labels = {}
    img_names = os.listdir(root +'/test-image/')
    img_names = sorted(img_names, key=lambda i: int(i.split(".")[0]))

    for img_name in img_names:
        img = os.path.join(root +'/test-image/', img_name)
        mask = os.path.join(root +'/test-mask/', img_name)
        imgs.append((img, mask))
    return imgs

class TN3K(data.Dataset):
    def __init__(self, mode, transform=None, return_size=False, fold=0):
        self.mode = mode
        root = './data/tn3k/'
        trainval = json.load(open(root + 'tn3k-trainval-fold'+str(fold)+'.json', 'r'))
        if mode == 'train':
            imgs = make_dataset(root, trainval['train'], 'trainval')
        elif mode == 'val':
            imgs = make_dataset(root, trainval['val'], 'trainval')
        elif mode == 'test':
            imgs = make_testset(root)

        self.imgs = imgs
        self.transform = transform
        self.return_size = return_size

    def __getitem__(self, item):
        if True:
            image_path, label_path = self.imgs[item]
            assert os.path.exists(image_path), ('{} does not exist'.format(image_path))
            assert os.path.exists(label_path), ('{} does not exist'.format(label_path))

            image = Image.open(image_path).convert('RGB')
            label = np.array(Image.open(label_path).convert('L'))
            label = label / label.max()
            label = Image.fromarray(label.astype(np.uint8))

            w, h = image.size
            size = (h, w)

            sample = {'image': image, 'label': label}

            if self.transform:
                sample = self.transform(sample)
            if self.return_size:
                sample['size'] = torch.tensor(size)

            label_name = os.path.basename(label_path)
            sample['label_name'] = label_name

            return sample
        # else:
        #     image_path = self.imgs[item]
        #     image = Image.open(image_path).convert('RGB')
        #     label = Image.open(image_path)
        #     sample = {'image': image, 'label': label}
        #     w, h = image.size
        #     size = (h, w)
        #
        #     if self.transform:
        #         sample = self.transform(sample)
        #     if self.return_size:
        #         sample['size'] = torch.tensor(size)
        #
        #     label_name = os.path.basename(image_path)
        #     sample['label_name'] = label_name
        #
        #     return sample

    def __len__(self):
        return len(self.imgs)

reason: import numpy


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

import numpy as np
import h5py
import torch
from torch.utils.data import Dataset
import torchvision
import random

from h5transform import H5RandomCrop, H5RandomHorizontalFlip, H5RandomVerticalFlip, H5RandomRotate

MEAN_FILE = "npy/mean_new.npy"
STD_FILE = "npy/std_new.npy"

class H5Dataset(Dataset):
    def __init__(self, h5file, istrain=False, data_mode="s2", label_mode="default"):
        super(H5Dataset, self).__init__()
        fid = h5py.File(h5file, 'r')

        self.s2 = fid['sen2']
        self.labels = fid['label']
        self.len = self.labels.shape[0]

        self.data_mode  = data_mode
        if data_mode == "s2":
            self.mean = np.load(MEAN_FILE).tolist()[-10:]
            self.std  = np.load(STD_FILE).tolist()[-10:]
        elif data_mode == "s1+s2":
            self.s1 = fid['sen1']
            idx = np.concatenate([np.array([4,5]), np.arange(8,18)]).astype(np.int)
            self.mean = np.load(MEAN_FILE)[idx].tolist()
            self.std  = np.load(STD_FILE)[idx].tolist()
        else:
            raise ValueError("data mode not found!")
        self.Normalize = torchvision.transforms.Normalize(self.mean, self.std)

        assert len(self.std) == len(self.mean), "shape mismatch error"

        self.istrain = istrain
        if istrain:
            self.randomCrop  = H5RandomCrop(32, 4)
            self.randomHflip = H5RandomHorizontalFlip()
            self.randomVflip = H5RandomVerticalFlip()
            self.randomRotate = H5RandomRotate()

        self.label_mode = label_mode
        if label_mode == "cluster":
            self._get_label = np.array([0,1,2,0,0,3,4,5,4,4, 6,7,7,8,7,7,9],dtype=np.int)

    def __getitem__(self, index):
        '''
            data: [32x32x18] numpy.ndarray

            return:
            data_torch: [18x32x32] torch.FloatTensor
        '''
        if self.data_mode == "s2":
            data = self.s2[index]
        elif self.data_mode == "s1+s2":
            s1 = self.s1[index][:,:,4:6]
            s2 = self.s2[index]
            data = np.concatenate([s1,s2],axis=-1)

        label = self.labels[index].argmax()

        data_torch = self.transform(data)

        if self.label_mode == "only_urban_rural":
            label = 0 if label < 10 else 1
        elif self.label_mode == "cluster":
            label = self.getLabel(label)

        return data_torch, label

    def __len__(self):
        return self.len

    def transform(self, data):
        """
        Input:
            data: ndarray image [32x32x18]
        Output:
            data_torch: Tensor image [18x32x32]
        """
        assert isinstance(data, np.ndarray), "data should be np.ndarray type"

        data_torch = torch.from_numpy(data).float().permute(2,0,1)
        if self.istrain:
            data_torch = self.randomCrop(data_torch)
            data_torch = self.randomHflip(data_torch)
            data_torch = self.randomVflip(data_torch)
            data_torch = self.randomRotate(data_torch)

        # maybe normalization after randomcrop or randomflip will be better
        data_torch = self.Normalize(data_torch)

        return data_torch

    def getLabel(self, label):
        return self._get_label[label]

class RoundDataset(Dataset):
    def __init__(self, h5file, data_mode="s2"):
        super(RoundDataset, self).__init__()
        fid = h5py.File(h5file, 'r')

        self.s1 = fid['sen1']
        self.s2 = fid['sen2']
        self.len = self.s2.shape[0]

        self.data_mode  = data_mode
        if data_mode == "s2":
            self.mean = np.load(MEAN_FILE).tolist()[-10:]
            self.std  = np.load(STD_FILE).tolist()[-10:]
        elif data_mode == "s1+s2":
            idx = np.concatenate([np.array([4,5]), np.arange(8,18)]).astype(np.int)
            self.mean = np.load(MEAN_FILE)[idx].tolist()
            self.std  = np.load(STD_FILE)[idx].tolist()
        else:
            raise ValueError("data mode not found!")
        self.Normalize = torchvision.transforms.Normalize(self.mean, self.std)

    def __getitem__(self, index):
        '''
            data: [32x32x18] numpy.ndarray

            return:
            data_torch: [18x32x32] torch.FloatTensor
        '''
        if self.data_mode == "s2":
            data = self.s2[index]
        elif self.data_mode == "s1+s2":
            s1 = self.s1[index][:,:,4:6]
            s2 = self.s2[index]
            data = np.concatenate([s1,s2],axis=-1)
        data_torch = self.transform(data)

        return data_torch

    def __len__(self):
        return self.len

    def transform(self, data):
        """
        Input:
            data: ndarray image [32x32x18]
        Output:
            data_torch: Tensor image [18x32x32]
        """
        assert isinstance(data, np.ndarray), "data should be np.ndarray type"

        data_torch = torch.from_numpy(data).float().permute(2,0,1)
        # maybe normalization after randomcrop or randomflip will be better
        data_torch = self.Normalize(data_torch)

        return data_torch

if __name__ == "__main__":
    trainset = H5Dataset("data/training.h5", data_mode="s1+s2")
    trainLoader = torch.utils.data.DataLoader(trainset, batch_size=10)
    for idx, (input, label) in enumerate(trainLoader):
        print(type(input))
        print(input.shape)
        print(input.dtype)
        print(label)
        import ipdb; ipdb.set_trace()


text:
from keras import Sequential
from keras.layers import Conv2D, Activation, MaxPooling2D, GlobalAveragePooling2D, Flatten, Dropout, Dense, BatchNormalization
import numpy as np
import keras, os
from keras.callbacks import ModelCheckpoint
import matplotlib.pyplot as plt
import breedsHelpFunctions as hf
from keras.models import load_model
from sklearn.metrics import confusion_matrix, classification_report, accuracy_score

#Name for the Model that will be saved
name = 'CNNnoAug120'

#TXT file where the list of pictures for each set are saved
traintxt = os.getcwd() + "/BreedTrain.txt"
testtxt = os.getcwd() + "/BreedTest.txt"
valtxt = os.getcwd() + "/BreedValidation.txt"

#Pictures in color = 3, in GreyScale = 1
color = 3

#Size of the picture
picsize = [100,100]

#Batch Size, epochs to run and optimizer
batch_size = 100
epochs = 30
optimizer = 'rmsprop'
classes = 6

############################################################
#----------------------- CODE ------------------------------
############################################################

print('#####################################################')
print('--------------- LOADING DATASETS --------------------')
print('#####################################################')

''' ---------- JUST TO LOAD ------------ '''

'''
#Dataset load, using the help function created and choosing if we want to include more pictures than the basics
# extratrain: [['TRANS100','TRANS100']]
# validsets: [['GB100','GB100']]
X_train, X_test, X_val, y_train, y_test, y_val = \
    hf.datasetExpansion(traintxt, testtxt, valtxt, extratrain = '', validsets = '')
print('Training sets loaded')

#Creating the arrays that we will use as input for our CNN
X_trainPic = hf.MatrixToDataset(X_train, picsize[0], picsize[1], color)
X_testPic = hf.MatrixToDataset(X_test, picsize[0], picsize[1], color)
X_valPic = hf.MatrixToDataset(X_val, picsize[0], picsize[1], color)
print('Training sets preprocessed')

#Creating output classes and dummies arrays to use as output
classes, y_tr = hf.create_dummies(y_train)
_, y_te = hf.create_dummies(y_test)
_, y_va = hf.create_dummies(y_val)
'''

print('#################################################################')
print('------------------------ MODEL STRUCTURE ------------------------')
print('#################################################################')

#Creating the model structure
keras.backend.clear_session()

model = Sequential()

model.add(Conv2D(16, (4, 4), padding='same', use_bias=True, activation='relu', input_shape=(picsize[0], picsize[1], color)))
model.add(MaxPooling2D(pool_size=(2, 2), strides=(2, 2), padding='same'))
model.add(Dropout(0.5))

model.add(Conv2D(32, (4, 4), padding='same', use_bias=True, activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2), strides=(2, 2), padding='same'))
model.add(Dropout(0.5))

model.add(Conv2D(64, (4, 4), padding='same', use_bias=True, activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2), strides=(2, 2), padding='same'))
model.add(Dropout(0.5))
#model.add(Conv2D(128, (4, 4), padding='same', use_bias=True, activation='relu'))
#model.add(Dropout(0.5))

model.add(Flatten())

model.add(Dense(500, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(len(classes), activation='softmax'))
model.summary()

print('Model Structure Created')

print('#####################################################')
print('--------------- MODEL CALLBACKS ---------------------')
print('#####################################################')

#Creting paths to save the model with best accuracy and loss
bestModelPath = os.getcwd() + '/bestModel/'
bestModelPathAcc = bestModelPath + 'model_acc_' + name + '.hdf5'
bestModelPathLoss = bestModelPath + 'model_loss_' + name + '.hdf5'

#Creating the checkpoint saver for our model
bestModelAcc = ModelCheckpoint(bestModelPathAcc, monitor="val_acc",
                      save_best_only=True, save_weights_only=False)
bestModelLoss = ModelCheckpoint(bestModelPathLoss, monitor="val_loss",
                      save_best_only=True, save_weights_only=False)

history=[]
history = hf.AccuracyHistory()

print('#####################################################')
print('----------------- MODEL COMPILE ---------------------')
print('#####################################################')

#optimizer depending on our previous selection to include in the model
if optimizer == 'adam':
    optCompile = keras.optimizers.Adam(lr=0.001)
else: optCompile = 'rmsprop'

#Compiling the model before running it
model.compile(loss=keras.losses.categorical_crossentropy,
              optimizer= optCompile,
              metrics=['accuracy'])
print('Model Compiled')

print('#####################################################')
print('------------------- MODEL RUN -----------------------')
print('#####################################################')

# Model Run
modelHistory = model.fit(X_trainPic, y_tr,
          batch_size=batch_size,
          epochs=epochs,
          verbose=1,
          validation_data=(X_testPic, y_te),
          callbacks=[history, bestModelAcc, bestModelLoss])

print('#####################################################')
print('----------------- MODEL EVALUATION ---------------------')
print('#####################################################')

#PLOT ACCURACY AND LOSS EVOLUTION FOR TRAINING AND TESTING SETS
plt.plot(range(len(history.acc)), history.acc)
plt.plot(range(len(history.acc)), history.val_acc)
plt.legend(['Train', 'Test'], loc='upper left')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.show()

plt.plot(range(len(history.loss)), history.loss)
plt.plot(range(len(history.loss)), history.val_loss)
plt.legend(['Train', 'Test'], loc='upper left')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.show()

#Load Accuracy and Loss models
best_model_acc = load_model(bestModelPathAcc)
best_model_loss = load_model(bestModelPathLoss)
print('Saved models loaded')

#Print the model evaluation for the Test and Validation Sets
#Include the confusion matrix and accuracy of the model
#We use the modelEvaluation help function
print('------------------------------------------------------------------------')
print('----------------------------- TEST -------------------------------------')
print('------------------------------------------------------------------------')
print('TEST: BEST ACCURACY MODEL')
_,_,_ = hf.modelEvaluation(best_model_acc, X_testPic, y_test, classes)
print('------------------------------------------------------------------------')
print('TEST: BEST LOSS MODEL')
_,_,_ = hf.modelEvaluation(best_model_loss, X_testPic, y_test, classes)
#print('------------------------------------------------------------------------')
#print('LAST MODEL')
#_,_,_ = hf.modelEvaluation(model, X_testPic, y_test, classes)

print('------------------------------------------------------------------------')
print('-------------------------- VALIDATION ----------------------------------')
print('------------------------------------------------------------------------')
print('VALIDATION: BEST ACCURACY MODEL')
_,_,_ = hf.modelEvaluation(best_model_acc, X_valPic, y_val, classes)
print('------------------------------------------------------------------------')
print('VALIDATION: BEST LOSS MODEL')
_,_,_ = hf.modelEvaluation(best_model_loss, X_valPic, y_val, classes)
#print('------------------------------------------------------------------------')
#print('LAST MODEL')
#_,_,_ = hf.modelEvaluation(model, X_valPic, y_val, classes)

print('------------------------------------------------------------------------')
print('-------------------- CNN ENSEMBLE VALIDATION ---------------------------')
print('------------------------------------------------------------------------')
#Ensemble model creation and evaluation
val_prob_ens = np.zeros([len(y_val), len(classes)])

val_prob_ens += best_model_acc.predict(X_valPic, verbose = 0)
val_prob_ens += best_model_loss.predict(X_valPic, verbose = 0)
val_prob_ens /= 2
pred = []

for p in val_prob_ens:
    max_value = max(p)
    max_index = np.where(p==max_value)
    pred += [classes[max_index[0][0]]]

print('CNN ENSEMBLE MODEL')
accuracy = round(accuracy_score(y_val, pred)*100,2)
print('-----------------------------\n')
print('Accuracy: ', accuracy, '%\n')
print(confusion_matrix(y_val, pred))
print(classification_report(y_val, pred))
print('-----------------------------')


text:
import os
import numpy as np
import pandas as pd
from datetime import timedelta, datetime
from matplotlib import pyplot as plt
from matplotlib import colors
from PIL import Image

def mkdirs(path):
    r"""
    create a new folder named `path` if not exists.
    """
    path = f'{path}'
    if not os.path.exists(path):
        os.makedirs(path)

def get_next_day(current_day):
    r"""
    Get the next day given the current day.
    """
    return str(pd.to_datetime(current_day) + timedelta(days=1))[:10]

def reshape_mask(mask, shape):
    r"""reshape the mask to align the data

    Parameters:
        mask (array)  -- a detected mask
        shape (tuple) -- the shape of the data

    Returns:
        the aligned mask
    """
    mask = Image.fromarray(mask).resize(shape, Image.ANTIALIAS)
    mask = np.fliplr(np.array(mask).T)
    return mask.astype(int)

def convert_matlab_time(matlab_time):
    r"""
    convert the Matlab time to Python time
    """
    try:
        return (datetime.fromordinal(int(matlab_time)) + timedelta(days=matlab_time % 1) - timedelta(days=366)).date()
    except:
        return np.nan

def time2num(tm):
    r"""Convert time to number.
    fit on numbers and then convert back the numbers to time
    """
    time = tm.time()
    minutes = time.hour * 60 + time.minute
    return minutes

def num2time(date, num):
    r"""Convert number to time.
    """
    year, month, day = date.split('-')
    hour = int(num / 60)
    minute = int(num - hour*60)
    return datetime(int(year), int(month), int(day), hour, minute)

max_stars_repo_name: cvvsu/maskNPF

max_issues_repo_name: cvvsu/maskNPF

max_forks_repo_name: cvvsu/maskNPF


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

import numpy as np

def generate_time_stamps(number_of_time_stamps, time_step) :
  """ Генерує значення табульованих часових відліків """
  if number_of_time_stamps <= 0 or time_step <= 0 :
    return np.array([])
  time_stamps = np.arange(0, number_of_time_stamps * time_step, time_step)
  return time_stamps

def generate_distances(
    number_of_vehicles, number_of_lanes, vehicle_length,
    road_interval, time_stamps,
    mean_speed, speed_deviation,
    average_number_of_vehicles) :
  """ Генерує відстані, які пройшли авто за кожен часовий інтервал """
  number_of_time_stamps = time_stamps.size
  time_step = time_stamps[1] - time_stamps[0]
  mean_speed *= time_step
  speed_deviation *= time_step
  shift = 2.0

  distances = np.zeros((number_of_vehicles, number_of_time_stamps))
  starting_positions = __generate_starting_positions(number_of_vehicles, average_number_of_vehicles,
      vehicle_length, road_interval, number_of_lanes)
  positions = np.copy(starting_positions)

  first_unexistent_vehicle_index = __find_first_unexistent_vehicle(positions)

  for i in range(0, number_of_time_stamps) :
    if first_unexistent_vehicle_index < number_of_vehicles :
      first_unexistent_vehicle_index += int(np.random.uniform(0, number_of_lanes + 1))
      if first_unexistent_vehicle_index > number_of_vehicles :
        first_unexistent_vehicle_index = number_of_vehicles
    for j in range(0, first_unexistent_vehicle_index) :
      if positions[j] == road_interval :
        continue
      distances[j][i] = np.random.normal(mean_speed, speed_deviation)
      if distances[j][i] < 0 :
        distances[j][i] = 0
      if (positions[j] + distances[j][i]) <= road_interval :
        positions[j] += distances[j][i]
      else :
        distances[j][i] = road_interval - positions[j]
        positions[j] = road_interval

  return [distances, positions, starting_positions]

def __generate_starting_positions(number_of_vehicles, average_number_of_vehicles,
    vehicle_length, road_interval, number_of_lanes) :
  """ Генерує початкові позиції автомобілів """
  starting_positions = np.zeros(number_of_vehicles)
  mean_spacing = ((road_interval*number_of_lanes) - (vehicle_length*average_number_of_vehicles)) / average_number_of_vehicles

  for i in range(0, number_of_vehicles, number_of_lanes) :
    for j in range(0, number_of_lanes) :
      if (i + j) >= number_of_vehicles :
        break
      position_shift = np.random.normal(mean_spacing, mean_spacing*0.2)
      position_shift = position_shift if (position_shift >= 0) else 0
      if (i - number_of_lanes) >= 0 :
        position = starting_positions[i - number_of_lanes] + position_shift
      else :
        position = position_shift
      starting_positions[i + j] = position if (position <= road_interval) else 0

  first_unexistent_vehicle_index = __find_first_unexistent_vehicle(starting_positions)
  reverse_index = first_unexistent_vehicle_index - 1

  for i in range(0, first_unexistent_vehicle_index) :
    if i >= reverse_index :
      break
    temp = starting_positions[i]
    starting_positions[i] = starting_positions[reverse_index]
    starting_positions[reverse_index] = temp
    reverse_index -= 1

  return starting_positions

def __find_first_unexistent_vehicle(starting_positions) :
  """ Знаходить останній автомобіль, який знаходиться в даному інтервалі дороги """
  vehicle_index = 0

  while vehicle_index < starting_positions.size :
    if starting_positions[vehicle_index] != 0 :
      vehicle_index += 1
    else :
      break

  return vehicle_index

def generate_tabled_density_values(max_density) :
  """ Генерує табульовані значення густина для фундаментальної діаграми """
  step = 0.001
  return np.arange(0, max_density + step, step)

max_stars_repo_name: LunaLovegoood/traffic-flow

max_issues_repo_name: LunaLovegoood/traffic-flow

max_forks_repo_name: LunaLovegoood/traffic-flow


text:
from r2base.index.index_bases import FieldOpBase
from r2base.mappings import TermScoreMapping
from r2base.processors.pipeline import Pipeline
from typing import Dict, List, Optional, Union
import logging
import numpy as np

class InvertedField(FieldOpBase):
    logger = logging.getLogger(__name__)

    ALPHA = 10.0
    MAX_NUM_COPY = 100

    @classmethod
    def to_mapping(self, mapping: TermScoreMapping):
        if mapping.mode == 'float':
            return {"type": "rank_features"}

        elif mapping.mode == 'int':
            return {"type": "text", "index_options": "freqs",
                    "analyzer": "cutter_analyzer", 'similarity': "tf_alone"}
        else:
            raise Exception("Unknown term score mode={}".format(mapping.mode))

    @classmethod
    def to_add_body(cls, mapping: TermScoreMapping, ts: Dict):
        if mapping.mode == 'float':
            ts = {term: np.exp(s) - 1.0 for term, s in ts.items()}
            return ts

        elif mapping.mode == 'int':
            tokens = []
            for term, s in ts.items():
                tokens.extend([term] * min(cls.MAX_NUM_COPY, int(s * cls.ALPHA)))
            return '/'.join(tokens)
        else:
            raise Exception("Unknown term score mode={}".format(mapping.mode))

    @classmethod
    def to_query_body(cls, key: str, mapping: TermScoreMapping,
                      tokens: Union[List[str], str],
                      top_k: int,
                      json_filter: Optional[Dict],
                      from_:int):

        if type(tokens) is str:
            if mapping.q_processor is None:
                raise Exception("q_processor is required for string query")

            kwargs = {'model_id': mapping.model_id}
            pipe = Pipeline(mapping.q_processor)
            tokens = pipe.run(tokens, **kwargs)

        if mapping.mode == 'float':
            main_query = [{'rank_feature': {'field': '{}.{}'.format(key, t),
                                            "log": {"scaling_factor": 1.0}
                                            }} for t in tokens]
            if json_filter is None:
                es_query = {
                    "_source": False,
                    "query": {"bool": {"should": main_query}},
                    "size": top_k,
                    "from": from_
                }
            else:
                es_query = {
                    "_source": False,
                    "query": {"bool": {"should": main_query, "filter": json_filter}},
                    "size": top_k,
                    "from": from_
                }
            return es_query
        elif mapping.mode == 'int':
            if json_filter is None:
                es_query = {
                    "_source": False,
                    "query": {"match": {key: '/'.join(tokens)}},
                    "size": top_k,
                    "from": from_
                }
            else:
                value = '/'.join(tokens)
                es_query = {
                    "_source": False,
                    "query": {"bool": {"must": {"match": {key: value}}, "filter": json_filter}},
                    "size": top_k,
                    "from": from_
                }
            return es_query
        else:
            raise Exception("Unknown term score mode={}".format(mapping.mode))

    @classmethod
    def process_score(cls, mapping: TermScoreMapping, score):
        if mapping.mode == 'int':
            score /= InvertedField.ALPHA
        return score


text:
import os
import math
import glob
import random
import importlib
from pathlib import Path
from collections import defaultdict

import torch
import torch.nn as nn
import numpy as np
from tqdm import tqdm
from tensorboardX import SummaryWriter

from optimizers import get_optimizer
from schedulers import get_scheduler
from utility.helper import count_parameters

SAMPLE_RATE = 16000

class Runner():
    """
    Used to handle high-level concepts of a ML experiment
    eg. training loop, evaluation loop, upstream propagation, optimization, tensorboard logging, checkpoint saving
    """
    def __init__(self, args, config):
        self.args = args
        self.config = config
        self.logger = SummaryWriter(args.expdir)

        self.init_ckpt = torch.load(self.args.past_exp, map_location='cpu') if self.args.past_exp else {}
        self.upstream = self._get_upstream()
        self.downstream = self._get_downstream()

        # set up the downstream name used by Tensorboard
        self.downstream_name = self.args.downstream
        if hasattr(self.downstream, 'get_downstream_name'):
            self.downstream_name = self.downstream.get_downstream_name()

    def _get_upstream(self):
        Upstream = getattr(importlib.import_module('hubconf'), self.args.upstream)
        upstream = Upstream(
            feature_selection = self.args.upstream_feature_selection,
            model_config = self.args.upstream_model_config,
            refresh = self.args.upstream_refresh,
            ckpt = self.args.upstream_ckpt,
        ).to(self.args.device)

        assert hasattr(upstream, 'forward')
        assert hasattr(upstream, 'get_output_dim')
        assert hasattr(upstream, 'get_downsample_rate')

        print(f'[Runner] - Upstream model architecture: {upstream}')
        print(f'[Runner] - Upstream output dimension: {upstream.get_output_dim()}')
        downsample = upstream.get_downsample_rate()
        print(f'[Runner] - Upstream downsample rate: {downsample} ({downsample / SAMPLE_RATE * 1000} ms/frame)')

        init_upstream = self.init_ckpt.get('Upstream')
        if init_upstream:
            print('[Runner] - Loading upstream weights from the previous experiment')
            upstream.load_state_dict(init_upstream)
        return upstream

    def _get_downstream(self):
        module_path = f'downstream.{self.args.downstream}.expert'
        Downstream = getattr(importlib.import_module(module_path), 'DownstreamExpert')
        downstream = Downstream(
            upstream_dim = self.upstream.get_output_dim(),
            **self.config,
            **vars(self.args)
        ).to(self.args.device)

        print(f'[Runner] - Downstream model architecture: {downstream}')
        print(f'[Runner] - Downstream has {count_parameters(downstream)} parameters')

        assert hasattr(downstream, 'get_train_dataloader')
        assert hasattr(downstream, 'get_dev_dataloader')
        assert hasattr(downstream, 'get_test_dataloader')
        assert hasattr(downstream, 'forward')
        assert hasattr(downstream, 'log_records')

        init_downstream = self.init_ckpt.get('Downstream')
        if init_downstream:
            print('[Runner] - Loading downstream weights from the previous experiment')
            downstream.load_state_dict(init_downstream)
        return downstream

    def _get_optimizer(self, model_params):
        optimizer = get_optimizer(
            model_params,
            self.config['runner']['total_steps'],
            self.config['optimizer']
        )

        init_optimizer = self.init_ckpt.get('Optimizer')
        if init_optimizer:
            print('[Runner] - Loading optimizer weights from the previous experiment')
            optimizer.load_state_dict(init_optimizer)
        return optimizer

    def _get_scheduler(self, optimizer):
        scheduler = get_scheduler(
            optimizer,
            self.config['runner']['total_steps'],
            self.config['scheduler']
        )

        init_scheduler = self.init_ckpt.get('Scheduler')
        if init_scheduler:
            print('[Runner] - Loading scheduler weights from the previous experiment')
            scheduler.load_state_dict(init_scheduler)
        return scheduler

    def train(self):
        # set model train/eval modes
        self.downstream.train()
        self.upstream.eval()
        if self.args.upstream_trainable:
            self.upstream.train()

        # set optimizer
        model_params = [self.downstream]
        if self.args.upstream_trainable:
            model_params.append(self.upstream)
        optimizer = self._get_optimizer(model_params)

        # set scheduler
        scheduler = None
        if self.config.get('scheduler'):
            scheduler = self._get_scheduler(optimizer)

        # set progress bar
        pbar = tqdm(total=self.config['runner']['total_steps'], dynamic_ncols=True, desc='overall')
        init_step = self.init_ckpt.get('Step')
        if init_step:
            pbar.n = init_step

        # prepare data
        dataloader = self.downstream.get_train_dataloader()

        all_loss = []
        backward_steps = 0
        records = defaultdict(list)
        prefix = f'{self.downstream_name}/train-'

        while pbar.n < pbar.total:
            for batch_id, (wavs, *others) in enumerate(tqdm(dataloader, dynamic_ncols=True, desc='train')):
                # try/except block for forward/backward
                try:
                    if pbar.n >= pbar.total:
                        break
                    global_step = pbar.n + 1

                    wavs = [wav.to(self.args.device) for wav in wavs]
                    if self.args.upstream_trainable:
                        features = self.upstream(wavs)
                    else:
                        with torch.no_grad():
                            features = self.upstream(wavs)

                    loss = self.downstream(
                        features, *others,
                        records = records,
                        logger = self.logger,
                        prefix = prefix,
                        global_step = global_step,
                        log_step = self.config['runner']['log_step'],
                        batch_id = batch_id,
                        batch_num = len(dataloader),
                    )
                    gradient_accumulate_steps = self.config['runner'].get('gradient_accumulate_steps')
                    (loss / gradient_accumulate_steps).backward()

                except RuntimeError as e:
                    if 'CUDA out of memory' in str(e):
                        print(f'[Runner] - CUDA out of memory at step {global_step}')
                        with torch.cuda.device(self.args.device):
                            torch.cuda.empty_cache()
                        optimizer.zero_grad()
                        continue
                    else:
                        raise

                # record loss
                all_loss.append(loss.item())
                del loss

                # whether to accumulate gradient
                backward_steps += 1
                if backward_steps % gradient_accumulate_steps > 0:
                    continue

                # gradient clipping
                paras = list(self.downstream.parameters())
                if self.args.upstream_trainable:
                    paras += list(self.upstream.parameters())
                grad_norm = torch.nn.utils.clip_grad_norm_(paras, self.config['runner']['gradient_clipping'])

                # optimize
                if math.isnan(grad_norm):
                    print(f'[Runner] - grad norm is NaN at step {global_step}')
                else:
                    optimizer.step()
                optimizer.zero_grad()

                # adjust learning rate
                if scheduler:
                    scheduler.step()

                # logging
                if global_step % self.config['runner']['log_step'] == 0:
                    # log loss
                    average_loss = torch.FloatTensor(all_loss).mean().item()
                    self.logger.add_scalar(f'{prefix}loss', average_loss, global_step=global_step)
                    all_loss = []

                    # log customized contents
                    self.downstream.log_records(
                        records = records,
                        logger = self.logger,
                        prefix = prefix,
                        global_step = global_step,
                        log_step = self.config['runner']['log_step'],
                    )
                    records = defaultdict(list)

                # evaluation and save checkpoint
                save_names = []

                if global_step % self.config['runner']['eval_step'] == 0:
                    for split in self.config['runner']['eval_dataloaders']:
                        save_names += self.evaluate(split, global_step)

                if global_step % self.config['runner']['save_step'] == 0:
                    def check_ckpt_num(directory):
                        max_keep = self.config['runner']['max_keep']
                        ckpt_pths = glob.glob(f'{directory}/states-*.ckpt')
                        if len(ckpt_pths) >= max_keep:
                            ckpt_pths = sorted(ckpt_pths, key=lambda pth: int(pth.split('-')[-1].split('.')[0]))
                            for ckpt_pth in ckpt_pths[:len(ckpt_pths) - max_keep + 1]:
                                os.remove(ckpt_pth)
                    check_ckpt_num(self.args.expdir)
                    save_names.append(f'states-{global_step}.ckpt')

                if len(save_names) > 0:
                    all_states = {
                        'Downstream': self.downstream.state_dict(),
                        'Optimizer': optimizer.state_dict(),
                        'Step': global_step,
                        'Args': self.args,
                        'Config': self.config,
                    }

                    if scheduler:
                        all_states['Scheduler'] = scheduler.state_dict()

                    if self.args.upstream_trainable:
                        all_states['Upstream'] = self.upstream.state_dict()

                    save_paths = [os.path.join(self.args.expdir, name) for name in save_names]
                    tqdm.write(f'[Runner] - Save the checkpoint to:')
                    for i, path in enumerate(save_paths):
                        tqdm.write(f'{i + 1}. {path}')
                        torch.save(all_states, path)

                pbar.update(1)

        Path(f'{self.args.expdir}/train_finished').touch(exist_ok=True)
        pbar.close()

    def evaluate(self, split='test', global_step=0):
        # fix seed to guarantee the same evaluation protocol across steps
        random.seed(self.args.seed)
        np.random.seed(self.args.seed)
        torch.manual_seed(self.args.seed)
        if torch.cuda.is_available():
            torch.cuda.manual_seed_all(self.args.seed)
        with torch.cuda.device(self.args.device):
            torch.cuda.empty_cache()

        # record original train/eval states and set all models to eval
        downstream_training = self.downstream.training
        upstream_training = self.upstream.training
        self.downstream.eval()
        self.upstream.eval()

        # prepare data
        dataloader = eval(f'self.downstream.get_{split}_dataloader')()

        # main evaluation block
        all_loss = []
        records = defaultdict(list)
        prefix = f'{self.downstream_name}/{split}-'

        for batch_id, (wavs, *others) in enumerate(tqdm(dataloader, dynamic_ncols=True, desc=split)):

            wavs = [wav.to(self.args.device) for wav in wavs]
            with torch.no_grad():
                features = self.upstream(wavs)

                loss = self.downstream(
                    features, *others,
                    records = records,
                    logger = self.logger,
                    prefix = prefix,
                    global_step = global_step,
                    log_step = self.config['runner']['log_step'],
                    batch_id = batch_id,
                    batch_num = len(dataloader),
                )
                all_loss.append(loss.item())

        # log loss
        average_loss = torch.FloatTensor(all_loss).mean().item()
        self.logger.add_scalar(f'{prefix}loss', average_loss, global_step=global_step)
        all_loss = []

        # log customized contents
        save_names = self.downstream.log_records(
            records = records,
            logger = self.logger,
            prefix = prefix,
            global_step = global_step,
            log_step = self.config['runner']['log_step'],
        )
        records = defaultdict(list)

        # prepare back to training
        with torch.cuda.device(self.args.device):
            torch.cuda.empty_cache()
        if downstream_training:
            self.downstream.train()
        if upstream_training:
            self.upstream.train()

        return [] if type(save_names) is not list else save_names


text:
import numpy as np

from dnnweaver2.simulator.stats import Stats

class LoopStack(object):
    def __init__(self, size=1024):
        self.loop_stack = np.empty(size, dtype=Instruction)
        self.head_ptr = 0
        self.tail_ptr = 0
        self.loop_count = 0
        self.mem_read_count = 0
        self.mem_write_count = 0

        self.mem_ops = []
        self.compute_ops = []

    def insert_compute(self, op, *args):
        self.compute_ops.append(ComputeInstruction(op, *args))

    def insert_mem_read(self, namespace, addr, size, stride, level=0, name=None):

        if name is None:
            name = 'mem_rd_{0}'.format(self.mem_read_count)

        mem_rd = MemoryReadInstruction(name=name,
                                       namespace=namespace,
                                       addr=addr,
                                       size=size,
                                       stride=stride,
                                       level=level)
        self.mem_ops.append(mem_rd)
        self.insert_instruction(mem_rd)

    def insert_mem_write(self, namespace, addr, size, stride, level=0, name=None):
        if name is None:
            name = 'mem_wr_{0}'.format(self.mem_write_count)

        mem_wr = MemoryWriteInstruction(name=name,
                                        namespace=namespace,
                                        addr=addr,
                                        size=size,
                                        stride=stride,
                                        level=level)
        self.mem_ops.append(mem_wr)
        self.insert_instruction(mem_wr)

    def insert_loop(self, loop_count, stride, level=0, name=None):
        if name is None:
            name = 'loop_{}'.format(self.loop_count)
        self.loop_count += 1
        loop = LoopInstruction(name=name,
                               loop_count=loop_count,
                               stride=stride,
                               level=level)

        self.insert_instruction(loop)

    def insert_instruction(self, inst):
        level = inst.level
        if self.loop_stack[self.head_ptr] is None:
            self.loop_stack[self.tail_ptr] = inst
            self.tail_ptr = 1
        else:
            l = self.loop_stack[self.head_ptr]
            assert l.level <= level
            if l.level < level:
                l.set_loop_with_level(level, inst)
                self.loop_stack[self.tail_ptr] = inst
                self.tail_ptr += 1
            else:
                self.head_ptr = self.tail_ptr
                self.loop_stack[self.tail_ptr] = inst
                self.tail_ptr += 1

    def promote_mem_ops(self, sram):
        for op in self.mem_ops:
            old = op.outer_loop
            if old is None:
                continue
            else:
                old.inner_loop.remove(op)
            curr = old
            namespace = op.namespace
            while curr is not None and (curr.stride[namespace] == 0 or sram[namespace] > op.size * curr.loop_count):
                if curr.stride[namespace] != 0:
                    op.size *= curr.loop_count
                curr = curr.outer_loop
            if curr is not None:
                if isinstance(op, MemoryWriteInstruction):
                    curr.inner_loop.append(op)
                else:
                    curr.inner_loop.insert(0, op)
                op.outer_loop = curr
                op.level = curr.level + 1
            else:
                op.level = 0
                op.outer_loop = None

    def __str__(self):
        ret = ''
        for l in self.loop_stack:
            if l is not None and l.level == 0:
                ret += '*' * 50 + '\n'
                ret += l.__str__() + '\n'
        ret += '*' * 50 + '\n'
        return ret

    def get_stats(self, acc_obj, verbose=False):
        stats = {}
        stats['total'] = Stats()
        for l in self.loop_stack:
            if l is not None and l.level == 0 and isinstance(l, LoopInstruction):
                stats[l.name] = l.get_stats(acc_obj, self.mem_ops, self.compute_ops, top=True)
                stats['total'] += stats[l.name]
        return stats

class Instruction(object):
    def __init__(self, type=None):
        self.inner_loop = []
        self.outer_loop = None

    def __str__(self, l=0):
        raise NotImplementedError('STR not implemented')

    def get_stats(self, acc_obj, mem_ops, compute_ops, top=False):
        raise NotImplementedError('get_stats not implemented')

class ComputeInstruction(Instruction):
    def __init__(self, op, *args):
        Instruction.__init__(self)
        self.op = op
        self.args = args

    def __str__(self, l=0):
        ret = 'Compute Op: Args = {}'.format(self.args)
        return ret

    def get_stats(self, acc_obj):
        stats = self.op(*self.args)
        return stats

class MemoryReadInstruction(Instruction):
    def __init__(self, name, namespace, addr, size, stride, level=0):
        Instruction.__init__(self)
        assert namespace in ['act', 'wgt', 'out']
        self.name = name
        self.namespace = namespace
        self.base_addr = addr
        self.size = size
        self.stride = stride
        self.level = level

    def __str__(self, l=0):
        ret = ' | ' * l + ('{0}: addr {1}, size {2}, stride {3}').format(self.name,
                                                                         self.base_addr,
                                                                         self.size,
                                                                         self.stride)
        ret += ', level: {}'.format(self.level)
        return ret

    def get_stats(self, acc_obj, mem_ops, compute_ops, top=False):
        stats = Stats()
        read_cycles = acc_obj.get_mem_read_cycles(self.namespace, self.size)
        if read_cycles is not None:
            stats.total_cycles = read_cycles
            stats.mem_stall_cycles = read_cycles
            stats.writes[self.namespace] = self.size
            stats.reads['dram'] = self.size
        return stats

class MemoryWriteInstruction(Instruction):
    def __init__(self, name, namespace, addr, size, stride, level=0):
        Instruction.__init__(self)
        assert namespace in ['act', 'wgt', 'out']
        self.name = name
        self.namespace = namespace
        self.base_addr = addr
        self.size = size
        self.stride = stride
        self.level = level

    def __str__(self, l=0):
        ret = ' | ' * l + ('{0}: addr {1}, size {2}, stride {3}').format(self.name,
                                                                         self.base_addr,
                                                                         self.size,
                                                                         self.stride)
        ret += ', level: {}'.format(self.level)
        return ret

    def get_stats(self, acc_obj, mem_ops, compute_ops, top=False):
        stats = Stats()
        write_cycles = acc_obj.get_mem_write_cycles(self.namespace, self.size)
        if write_cycles is not None:
            stats.total_cycles = write_cycles
            stats.mem_stall_cycles = write_cycles
            stats.reads[self.namespace] = self.size
            stats.writes['dram'] = self.size
        return stats

class LoopInstruction(Instruction):
    def __init__(self, name, loop_count, stride, level=0):
        Instruction.__init__(self)
        self.name = name
        self.loop_count = loop_count
        self.stride = stride
        self.level = level

    def __str__(self, l=0):
        ret = ' | ' * l + '{0}: Range {1}, stride {2}'.format(self.name,
                                                              self.loop_count,
                                                              self.stride)
        for il in self.inner_loop:
            ret += '\n{}'.format(il.__str__(l + 1))
        return ret

    def set_loop_with_level(self, level, loop):
        assert level > self.level
        if level == self.level + 1:
            self.inner_loop.append(loop)
            loop.outer_loop = self
        else:
            if len(self.inner_loop) == 0:
                raise IndexError(
                    'Inserting loop {0} with level {1} into loop {2} with level {3}'.format(loop.name, loop.level,
                                                                                self.name, self.level))
            self.inner_loop[-1].set_loop_with_level(level, loop)

    def get_pipe_stats(self, acc_obj, mem_ops, compute_ops, top=False):

        rd_stats = Stats()
        pipe_stats = []
        compute_stats = Stats()
        wr_stats = Stats()

        for il in self.inner_loop:
            if isinstance(il, MemoryReadInstruction):
                rd_stats += il.get_stats(acc_obj, mem_ops, compute_ops) * self.loop_count
            elif isinstance(il, MemoryWriteInstruction):
                wr_stats += il.get_stats(acc_obj, mem_ops, compute_ops) * self.loop_count

        il_count = 0
        for il in self.inner_loop:
            if isinstance(il, LoopInstruction):
                il_count += 1

        if il_count == 0:
            for cop in compute_ops:
                compute_stats += cop.get_stats(acc_obj)
            pipeline = Pipeline(compute_stats.total_cycles, rd_stats.total_cycles//self.loop_count, wr_stats.total_cycles//self.loop_count, self.loop_count)
            total_stats = rd_stats + compute_stats * self.loop_count + wr_stats
            total_stats.total_cycles = pipeline.get_cycles()
            total_stats.mem_stall_cycles = total_stats.total_cycles - compute_stats.total_cycles * self.loop_count
        else:
            for il in self.inner_loop:
                count = 0
                if isinstance(il, LoopInstruction):
                    count += 1
                    assert count == 1
                    pipe, il_stats = il.get_pipe_stats(acc_obj, mem_ops, compute_ops)
                    pipeline = Pipeline(pipe, rd_stats.total_cycles//self.loop_count, wr_stats.total_cycles//self.loop_count, self.loop_count)
                    total_stats = rd_stats + il_stats * self.loop_count + wr_stats
                    total_stats.total_cycles = pipeline.get_cycles()
                    total_stats.mem_stall_cycles = total_stats.total_cycles - (il_stats.total_cycles - il_stats.mem_stall_cycles) * self.loop_count

        # print(self.name)
        # print(pipeline)
        # print(total_stats)

        return pipeline, total_stats

    def get_stats(self, acc_obj, mem_ops, compute_ops, top=True):
        pipe, stats = self.get_pipe_stats(acc_obj, mem_ops, compute_ops, top=True)
        if top:
            for mop in mem_ops:
                if mop.level == 0:
                    stats += mop.get_stats(acc_obj, mem_ops, compute_ops)
        return stats


text:
# -*- coding: utf-8 -*-
# /usr/bin/env python2.7 + brainvisa compliant env
#
#  This software and supporting documentation are distributed by
#      Institut Federatif de Recherche 49
#      CEA/NeuroSpin, Batiment 145,
#      91191 Gif-sur-Yvette cedex
#      France
#
# This software is governed by the CeCILL license version 2 under
# French law and abiding by the rules of distribution of free software.
# You can  use, modify and/or redistribute the software under the
# terms of the CeCILL license version 2 as circulated by CEA, CNRS
# and INRIA at the following URL "http://www.cecill.info".
#
# As a counterpart to the access to the source code and  rights to copy,
# modify and redistribute granted by the license, users are provided only
# with a limited warranty  and the software's author,  the holder of the
# economic rights,  and the successive licensors  have only  limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading,  using,  modifying and/or developing or reproducing the
# software by the user in light of its specific status of free software,
# that may mean  that it is complicated to manipulate,  and  that  also
# therefore means  that it is reserved for developers  and  experienced
# professionals having in-depth computer knowledge. Users are therefore
# encouraged to load and test the software's suitability as regards their
# requirements in conditions enabling the security of their systems and/or
# data to be ensured and,  more generally, to use and operate it in the
# same conditions as regards security.
#
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license version 2 and that you accept its terms.

"""Creating pickle file from T1 MRI datas

The aim of this script is to create dataset of cropped skeletons from MRIs
saved in a .pickle file.
Several steps are required: normalization, crop and .pickle generation

  Typical usage
  -------------
  You can use this program by first entering in the brainvisa environment
  (here brainvisa 5.0.0 installed with singurity) and launching the script
  from the terminal:
  >>> bv bash
  >>> python dataset_gen_pipe.py

  Alternatively, you can launch the script in the interactive terminal ipython:
  >>> %run dataset_gen_pipe.py

"""

from __future__ import division
from __future__ import print_function

import argparse
import sys
import os
from os import listdir
from os.path import join

import numpy as np

import six

from deep_folding.anatomist_tools.utils.logs import LogJson
from deep_folding.anatomist_tools.utils.load_bbox import compute_max_box
from deep_folding.anatomist_tools.utils.resample import resample
from deep_folding.anatomist_tools.utils.sulcus_side import complete_sulci_name
from deep_folding.anatomist_tools.load_data import fetch_data

_ALL_SUBJECTS = -1

_SIDE_DEFAULT = 'L'  # hemisphere 'L' or 'R'

_INTERP_DEFAULT = 'nearest'  # default interpolation for ApplyAimsTransform

_RESAMPLING_DEFAULT = None # if None, resampling method is AimsApplyTransform

_OUT_VOXEL_SIZE = (1, 1, 1) # default output voxel size for Bastien's resampling

# sulcus to encompass:
# its name depends on the hemisphere side
_SULCUS_DEFAULT = 'S.T.s.ter.asc.ant.'

# Input directories
# -----------------

# Input directory contaning the morphologist analysis of the HCP database
_SRC_DIR_DEFAULT = '/neurospin/hcp'

# Directory where subjects to be processed are stored.
# Default is for HCP dataset
_MORPHOLOGIST_DIR_DEFAULT = 'ANALYSIS/3T_morphologist'

# Directory that contains the transformation file
# from native to MNI through SPM
# These files have been created with spm_skeleton
_TRANSFORM_DIR_DEFAULT = '/neurospin/dico/data/deep_folding/data/transform'

# Directory containing bounding box json files
# default corresponds to bounding boxes computed for voxels of 1mm
_BBOX_DIR_DEFAULT = '/neurospin/dico/data/deep_folding/data/bbox'

# Output (target) directory
# -------------------------
_TGT_DIR_DEFAULT = '/neurospin/dico/data/deep_folding/test'

class DatasetCroppedSkeleton:
    """Generates cropped skeleton files and corresponding pickle file
    """

    def __init__(self, src_dir=_SRC_DIR_DEFAULT,
                 tgt_dir=_TGT_DIR_DEFAULT,
                 transform_dir=_TRANSFORM_DIR_DEFAULT,
                 bbox_dir=_BBOX_DIR_DEFAULT,
                 morphologist_dir=_MORPHOLOGIST_DIR_DEFAULT,
                 list_sulci=_SULCUS_DEFAULT,
                 side=_SIDE_DEFAULT,
                 interp=_INTERP_DEFAULT,
                 resampling=_RESAMPLING_DEFAULT,
                 out_voxel_size=_OUT_VOXEL_SIZE):
        """Inits with list of directories and list of sulci

        Args:
            src_dir: list of strings naming ful path source directories,
                    containing MRI images
            tgt_dir: name of target (output) directory with full path
            transform_dir: directory containing transformation files
                    (generated using transform.py)
            bbox_dir: directory containing bbox json files
                    (generated using bounding_box.py)
            list_sulci: list of sulcus names
            side: hemisphere side (either L for left, or R for right hemisphere)
            interp: string giving interpolation for AimsApplyTransform
        """

        self.src_dir = src_dir
        self.side = side
        # Transforms sulcus in a list of sulci
        self.list_sulci = ([list_sulci] if isinstance(list_sulci, str)
                           else list_sulci)
        self.list_sulci = complete_sulci_name(self.list_sulci, self.side)

        self.tgt_dir = tgt_dir
        self.transform_dir = transform_dir
        self.bbox_dir = bbox_dir
        self.morphologist_dir = morphologist_dir
        self.interp = interp
        self.resampling = resampling
        self.out_voxel_size = out_voxel_size

        # Morphologist directory
        self.morphologist_dir = join(self.src_dir, self.morphologist_dir)
        # default acquisition subdirectory
        self.acquisition_dir = "%(subject)s/t1mri/default_acquisition"
        # (input) name of normalized SPM file
        self.normalized_spm_file = "normalized_SPM_%(subject)s.nii"

        # Directory where to store cropped files
        self.cropped_dir = join(self.tgt_dir, self.side + 'crops')

        # Names of files in function of dictionary: keys -> 'subject' and 'side'
        # Files from morphologist pipeline
        self.normalized_spm_file = 'normalized_SPM_%(subject)s.nii'
        self.skeleton_file = 'default_analysis/segmentation/' \
                             '%(side)sskeleton_%(subject)s.nii.gz'

        # Names of files in function of dictionary: keys -> 'subject' and 'side'
        self.transform_file = 'natif_to_template_spm_%(subject)s.trm'
        self.cropped_file = '%(subject)s_normalized.nii.gz'

        # Initialization of bounding box coordinates
        self.bbmin = np.zeros(3)
        self.bbmax = np.zeros(3)

        # Creates json log class
        json_file = join(self.tgt_dir, self.side + 'dataset.json')
        self.json = LogJson(json_file)

    def crop_one_file(self, subject_id):
        """Crops one file

        Args:
            subject_id: string giving the subject ID
        """

        # Identifies 'subject' in a mapping (for file and directory namings)
        subject = {'subject': subject_id, 'side': self.side}
        print(subject_id)

        # Names directory where subject analysis files are stored
        subject_dir = \
            join(self.morphologist_dir, self.acquisition_dir % subject)

        # Transformation file name
        file_transform = join(self.transform_dir, self.transform_file % subject)

        # Normalized SPM file name
        file_SPM = join(subject_dir, self.normalized_spm_file % subject)

        # Skeleton file name
        file_skeleton = join(subject_dir, self.skeleton_file % subject)
        if os.path.exists(file_skeleton) and os.path.exists(file_transform):
            # Creates output (cropped) file name
            file_cropped = join(self.cropped_dir, self.cropped_file % subject)

            # Normalization and resampling of skeleton images
            if self.resampling:
                resample(file_skeleton,
                        file_cropped,
                        output_vs=self.out_voxel_size,
                        transformation=file_transform)

            else :
                cmd_normalize = 'AimsApplyTransform' + \
                                ' -i ' + file_skeleton + \
                                ' -o ' + file_cropped + \
                                ' -m ' + file_transform + \
                                ' -r ' + file_SPM + \
                                ' -t ' + self.interp
                os.system(cmd_normalize)

            # Take the coordinates of the bounding box
            bbmin = self.bbmin
            bbmax = self.bbmax
            xmin, ymin, zmin = str(bbmin[0]), str(bbmin[1]), str(bbmin[2])
            xmax, ymax, zmax = str(bbmax[0]), str(bbmax[1]), str(bbmax[2])

            # Crop of the images based on bounding box
            cmd_bounding_box = ' -x ' + xmin + ' -y ' + ymin + ' -z ' + zmin + \
                           ' -X ' + xmax + ' -Y ' + ymax + ' -Z ' + zmax
            cmd_crop = 'AimsSubVolume' + \
                   ' -i ' + file_cropped + \
                   ' -o ' + file_cropped + cmd_bounding_box
            os.system(cmd_crop)

    def crop_files(self, number_subjects=_ALL_SUBJECTS):
        """Crop nii files

        The programm loops over all subjects from the input (source) directory.

        Args:
            number_subjects: integer giving the number of subjects to analyze,
                by default it is set to _ALL_SUBJECTS (-1).
        """

        if number_subjects:

            # subjects are detected as the directory names under src_dir
            list_all_subjects = [dI for dI in os.listdir(self.morphologist_dir)\
             if os.path.isdir(os.path.join(self.morphologist_dir,dI))]

            # Gives the possibility to list only the first number_subjects
            list_subjects = (
                list_all_subjects
                if number_subjects == _ALL_SUBJECTS
                else list_all_subjects[:number_subjects])

            # Creates target and cropped directory
            if not os.path.exists(self.tgt_dir):
                os.makedirs(self.tgt_dir)
            if not os.path.exists(self.cropped_dir):
                os.makedirs(self.cropped_dir)

            # Writes number of subjects and directory names to json file
            dict_to_add = {'nb_subjects': len(list_subjects),
                           'src_dir': self.src_dir,
                           'transform_dir': self.transform_dir,
                           'bbox_dir': self.bbox_dir,
                           'side': self.side,
                           'interp': self.interp,
                           'list_sulci': self.list_sulci,
                           'bbmin': self.bbmin.tolist(),
                           'bbmax': self.bbmax.tolist(),
                           'tgt_dir': self.tgt_dir,
                           'cropped_dir': self.cropped_dir,
                           'resampling_type': 'AimsApplyTransform' if self.resampling is None else 'Bastien',
                           'out_voxel_size': self.out_voxel_size
                           }
            self.json.update(dict_to_add=dict_to_add)

            for subject in list_subjects:
                self.crop_one_file(subject)

    def dataset_gen_pipe(self, number_subjects=_ALL_SUBJECTS):
        """Main API to create pickle files

        The programm loops over all subjects from the input (source) directory.

        Args:
            number_subjects: integer giving the number of subjects to analyze,
                by default it is set to _ALL_SUBJECTS (-1).
        """

        self.json.write_general_info()

        # Generate cropped files
        if number_subjects:
            self.bbmin, self.bbmax = compute_max_box(sulci_list=self.list_sulci,
                                                     side=self.side,
                                                     talairach_box=False,
                                                     src_dir=self.bbox_dir)
        # Generate cropped files
        self.crop_files(number_subjects=number_subjects)
        # Creation of .pickle file for all subjects
        if number_subjects:
            fetch_data(cropped_dir=self.cropped_dir,
                       tgt_dir=self.tgt_dir,
                       side=self.side)

def parse_args(argv):
    """Function parsing command-line arguments

    Args:
        argv: a list containing command line arguments

    Returns:
        params: dictionary with keys: src_dir, tgt_dir, nb_subjects, list_sulci
    """

    # Parse command line arguments
    parser = argparse.ArgumentParser(
        prog='dataset_gen_pipe.py',
        description='Generates cropped and pickle files')
    parser.add_argument(
        "-s", "--src_dir", type=str, default=_SRC_DIR_DEFAULT,
        help='Source directory where the MRI data lies. '
             'Default is : ' + _SRC_DIR_DEFAULT)
    parser.add_argument(
        "-t", "--tgt_dir", type=str, default=_TGT_DIR_DEFAULT,
        help='Target directory where to store the cropped and pickle files. '
             'Default is : ' + _TGT_DIR_DEFAULT)
    parser.add_argument(
        "-r", "--transform_dir", type=str, default=_TRANSFORM_DIR_DEFAULT,
        help='Transform directory where transformation files from native '
             'to Talairach files have been stored. '
             'Default is : ' + _TRANSFORM_DIR_DEFAULT)
    parser.add_argument(
        "-b", "--bbox_dir", type=str, default=_BBOX_DIR_DEFAULT,
        help='Bounding box directory where json files containing '
             'bounding box coordinates have been stored. '
             'Default is : ' + _BBOX_DIR_DEFAULT)
    parser.add_argument(
        "-m", "--morphologist_dir", type=str, default=_MORPHOLOGIST_DIR_DEFAULT,
        help='Directory where subjects to be processed are stored')
    parser.add_argument(
        "-u", "--sulcus", type=str, default=_SULCUS_DEFAULT, nargs='+',
        help='Sulcus name around which we determine the bounding box. '
             'If there are several sulci, add all sulci '
             'one after the other. Example: -u sulcus_1 sulcus_2 '
             'Default is : ' + _SULCUS_DEFAULT)
    parser.add_argument(
        "-i", "--side", type=str, default=_SIDE_DEFAULT,
        help='Hemisphere side (either L or R). Default is : ' + _SIDE_DEFAULT)
    parser.add_argument(
        "-n", "--nb_subjects", type=str, default="all",
        help='Number of subjects to take into account, or \'all\'. '
             '0 subject is allowed, for debug purpose.'
             'Default is : all')
    parser.add_argument(
        "-e", "--interp", type=str, default=_INTERP_DEFAULT,
        help="Same interpolation type as for AimsApplyTransform. "
             "Type of interpolation used for Volumes: "
             "n[earest], l[inear], q[uadratic], c[cubic], quartic, "
             "quintic, six[thorder], seven[thorder]. "
             "Modes may also be specified as order number: "
             "0=nearest, 1=linear...")
    parser.add_argument(
        "-p", "--resampling", type=str, default=None,
        help='Method of resampling to perform. '
             'Type of resampling: '
             's[ulcus] for Bastien method'
             'If None, AimsApplyTransform is used.'
             'Default is : None')
    parser.add_argument(
        "-v", "--out_voxel_size", type=int, nargs='+', default=_OUT_VOXEL_SIZE,
        help='Voxel size of output images'
             'Default is : 1 1 1')

    params = {}

    args = parser.parse_args(argv)
    params['src_dir'] = args.src_dir
    params['tgt_dir'] = args.tgt_dir
    params['bbox_dir'] = args.bbox_dir
    params['transform_dir'] = args.transform_dir
    params['list_sulci'] = args.sulcus  # a list of sulci
    params['side'] = args.side
    params['interp'] = args.interp
    params['resampling'] = args.resampling
    params['out_voxel_size'] = tuple(args.out_voxel_size)
    params['morphologist_dir'] = args.morphologist_dir

    number_subjects = args.nb_subjects

    # Check if nb_subjects is either the string "all" or a positive integer
    try:
        if number_subjects == "all":
            number_subjects = _ALL_SUBJECTS
        else:
            number_subjects = int(number_subjects)
            if number_subjects < 0:
                raise ValueError
    except ValueError:
        raise ValueError(
            "number_subjects must be either the string \"all\" or an integer")
    params['nb_subjects'] = number_subjects

    return params

def dataset_gen_pipe(src_dir=_SRC_DIR_DEFAULT, tgt_dir=_TGT_DIR_DEFAULT,
                     transform_dir=_TRANSFORM_DIR_DEFAULT,
                     bbox_dir=_BBOX_DIR_DEFAULT,
                     morphologist_dir=_MORPHOLOGIST_DIR_DEFAULT,
                     side=_SIDE_DEFAULT, list_sulci=_SULCUS_DEFAULT,
                     number_subjects=_ALL_SUBJECTS, interp=_INTERP_DEFAULT,
                     resampling=_RESAMPLING_DEFAULT,
                     out_voxel_size=_OUT_VOXEL_SIZE):
    """Main program generating cropped files and corresponding pickle file
    """

    dataset = DatasetCroppedSkeleton(src_dir=src_dir, tgt_dir=tgt_dir,
                                     transform_dir=transform_dir,
                                     bbox_dir=bbox_dir,
                                     morphologist_dir=morphologist_dir,
                                     side=side, list_sulci=list_sulci,
                                     interp=interp, resampling=resampling,
                                     out_voxel_size=out_voxel_size)
    dataset.dataset_gen_pipe(number_subjects=number_subjects)

def main(argv):
    """Reads argument line and creates cropped files and pickle file

    Args:
        argv: a list containing command line arguments
    """

    # This code permits to catch SystemExit with exit code 0
    # such as the one raised when "--help" is given as argument
    try:
        # Parsing arguments
        params = parse_args(argv)
        # Actual API
        dataset_gen_pipe(src_dir=params['src_dir'],
                         tgt_dir=params['tgt_dir'],
                         transform_dir=params['transform_dir'],
                         bbox_dir=params['bbox_dir'],
                         morphologist_dir=params['morphologist_dir'],
                         side=params['side'],
                         list_sulci=params['list_sulci'],
                         interp=params['interp'],
                         number_subjects=params['nb_subjects'],
                         resampling=params['resampling'],
                         out_voxel_size=params['out_voxel_size'])
    except SystemExit as exc:
        if exc.code != 0:
            six.reraise(*sys.exc_info())

######################################################################
# Main program
######################################################################

if __name__ == '__main__':
    # This permits to call main also from another python program
    # without having to make system calls
    main(argv=sys.argv[1:])


text:
import json
import gzip, glob
import pandas as pd
import dask.dataframe as dd
import os, csv
import numpy as np
import multiprocessing
from collections import defaultdict
from functools import partial
import codecs
from dask.diagnostics import ProgressBar

def step_1():
    import sys
    print(sys.stdout.encoding)
    header = '/mnt/e/MAG/mag-2021-01-05/'
    papers = 'mag/Papers.txt'

    for chunk in pd.read_csv(header+papers, header=None, sep='\t', encoding='utf-8',
                             quoting=csv.QUOTE_NONE,
                            chunksize=100000):# , #error_bad_lines=False):
        valid_chunk = chunk[[0, 2, 7, 19]]

        valid_chunk.to_csv('data_temp/papers_filtered_1.txt', sep='\t', mode='a',
                           header=False, index=False)

def step_2():
    # sort --parallel=20 -Tdata_temp -t$'\t' -nk1 data_temp/papers_filtered_1.txt > data/papers_filtered_1_byPaperID.txt
    pass

def step_3():
    header = '/mnt/e/MAG/mag-2021-01-05/'
    paper_refs = 'mag/PaperReferences.txt'
    paper_refs = dd.read_csv(header+paper_refs, header=None, sep='\t', names = ['paper_id', 'paper_reference_id'])
#                          dtype={'paper_id': 'Int64', 'paper_reference_id': 'Int64'})

    papers = dd.read_csv('data/papers_filtered_1_byPaperID.txt', header=None, sep='\t',
                         names=['paper_id', 'doi', 'year', 'citation_count'],
                         dtype={'citation_count': 'object',
                           'year': 'object'})

    paper_refs = paper_refs.set_index('paper_id', sorted=True)
    papers = papers.set_index('paper_id', sorted=True)
    join_result = paper_refs.merge(papers, left_index=True, right_index=True)
    join_result.to_csv('data_temp/paper_refs_year.txt', sep='\t', single_file=True, header=False)

def step_4():
    # sort --parallel=20 -Tdata_temp -t$'\t' -nk2 data_temp/paper_refs_year.txt > data_temp/paper_refs_year_byReferenceID.txt
    pass

def step_5():
    # split -l 1000000 -a 5 -d /mnt/e/MAG/mag-2021-01-05/mag/PaperAuthorAffiliations.txt data_temp/AuthorAffili_split/author_affil_
    # split -l 1000000 -a 5 -d data_temp/paper_refs_year_byReferenceID.txt data_temp/PaperRefsYear_split/paper_ref_year_
    pass
# IF author_affil is .gz:

#     header = '/mnt/e/MAG/mag-2021-01-05/'
#     authors_affil = 'mag/PaperAuthorAffiliations.txt.gz'
#     i = 0
#     j = 0
#     with gzip.open(header+authors_affil, 'rb') as infile:
#         to_write = ''
#         for line in infile:
#             to_write += line.decode('utf-8')
#             i += 1

#             if i%1000000 == 0:
#                 f = open('data_temp/AuthorAffili_split/author_affili_%d.txt' % j, 'w')
#                 f.write(to_write)
#                 f.close()
#                 to_write = ''
#                 j += 1

#         if len(to_write) > 0:
#             f = open('data_temp/AuthorAffili_split/author_affili_%d.txt' % j, 'w')
#             f.write(to_write)
#             f.close()

def json_years(years):
    temp = defaultdict(lambda:0)
    for year in years:
        temp[str(year)[:4]] += 1
    temp = json.dumps(temp)
    return temp

def concat_json(json_1, values_2):
    json_dict = json.loads(json_1)
    for value in values_2:
        key = str(value)[:4]
        if key in json_dict:
            json_dict[key] += 1
        else:
            json_dict[key] = 1
    return json.dumps(json_dict)

def papers_citations(inputfile):
    i = int(inputfile.split('_')[-1])
    chunk = pd.read_csv('data_temp/PaperRefsYear_split/'+inputfile, header=None,
                sep='\t', names=['paper_id', 'reference_id', 'doi', 'year', 'citation'])

    output = open("data_temp/PaperCits_split/paper_cit_%05d" % i, 'w')

    current_reference = -1
    years = []
    t_citation = 0

    output_text = ''
    if current_reference == -1:
        current_reference = chunk.iloc[0,1]
    for idx,row in chunk.iterrows():
        if current_reference == row['reference_id']:
            years.append(str(row['year']))
            t_citation += 1
        else:
            output.write("%d\t%d\t%s\n" % (current_reference, t_citation, json_years(years)))
            current_reference = row['reference_id']
            years = [str(row['year'])]
            t_citation = 1

    output.write("%d\t%d\t%s" % (current_reference, t_citation, json_years(years)))
    output.close()

def step_6():
    pool = multiprocessing.Pool(16)
    inputs = os.listdir('data_temp/PaperRefsYear_split/')

    pool.map(papers_citations, inputs)
    pool.close()

def step_7():
    inputs = sorted(glob.glob('data_temp/PaperCits_split/*'))

    panda_concat_final = pd.read_csv(inputs[0], header=None,
                            sep='\t', names=['paper_id', 'cits', 'years'])

    for inputs_idx in inputs[1:]:
        panda_processed_i1 = pd.read_csv(inputs_idx, header=None,
                            sep='\t', names=['paper_id', 'cits', 'years'])

        last_id = panda_concat_final.iloc[-1, 0]
        first_id = panda_processed_i1.iloc[0, 0]
        if last_id == first_id:
            # união
            panda_processed_i1.iloc[0, 1] += panda_concat_final.iloc[-1, 1]
            panda_processed_i1.iloc[0, 2] = concat_json(panda_processed_i1.iloc[0, 2], panda_concat_final.iloc[-1, 2])
            panda_concat_final = panda_concat_final[:-1]

        panda_concat_final = panda_concat_final.append(panda_processed_i1, ignore_index=True)

    panda_concat_final.to_csv('data/paper_cits.csv', header=None, index=None, sep='\t')

def join_authors_of_paper(input_name):
    print(input_name)
    chunk = pd.read_csv(input_name, header=None, index_col=False, error_bad_lines=False,
                    encoding='utf-8', quoting=csv.QUOTE_NONE,
                    sep='\t', names=['paper_id', 'author_id', '11', '22', '33'])

    authors = []

    output_text = ''
    current_reference = chunk.iloc[0,0]
    i = int(input_name.split('_')[-1])
    output = open("data_temp/Colabs_split/paper_authors_%05d" % i, 'w')
    for idx,row in chunk.iterrows():
        if current_reference == row['paper_id']:
            authors.append(str(row['author_id']))
        else:
            output.write("%d\t%s\n" % (current_reference, ','.join(authors)))
            current_reference = row['paper_id']
            authors = [str(row['author_id'])]

    if len(authors) > 0:
        output.write("%d\t%s" % (current_reference, ','.join(authors)))

    output.close()

def step_8():
    pool = multiprocessing.Pool(16)

    inputs = glob.glob('data_temp/AuthorAffili_split/*')

#     input_erros = ['data_temp/AuthorAffili_split/author_affil_00066',
#     'data_temp/AuthorAffili_split/author_affil_00110',
#     'data_temp/AuthorAffili_split/author_affil_00090',
#     'data_temp/AuthorAffili_split/author_affil_00147',
#     'data_temp/AuthorAffili_split/author_affil_00137',
#     'data_temp/AuthorAffili_split/author_affil_00075',
#     'data_temp/AuthorAffili_split/author_affil_00264',
#     'data_temp/AuthorAffili_split/author_affil_00302',
#     'data_temp/AuthorAffili_split/author_affil_00293',
#     'data_temp/AuthorAffili_split/author_affil_00251',
#     'data_temp/AuthorAffili_split/author_affil_00378',
#     'data_temp/AuthorAffili_split/author_affil_00444']

    pool.map(join_authors_of_paper,inputs)
    pool.close()

def step_9():
    # TODO VERIFICAR SE ESTÁ SORTED SEM USAR SORTED
    inputs = glob.glob('data_temp/Colabs_split/*')
    inputs = sorted(inputs)

    panda_concat_final = pd.read_csv(inputs[0], header=None, index_col=False,
                            sep='\t', names=['paper_id', 'authors'])

    for inputs_idx in inputs[1:]:

        last_id = panda_concat_final.iloc[-1, 0]
        last_refs = panda_concat_final[panda_concat_final['paper_id'] == last_id]['authors']
        n_authors = len(last_refs)

        panda_processed_i1 = pd.read_csv(inputs_idx, header=None, index_col=False,
                            sep='\t', names=['paper_id', 'authors'])
        first_id = panda_processed_i1.iloc[0, 0]

        if last_id == first_id:
            # união
            new_authors = panda_processed_i1.iloc[0, 1] + ',' + ','.join(last_refs)
            panda_processed_i1.iloc[0, 1] = new_authors
            panda_concat_final = panda_concat_final[:-1]
            # contatena

        panda_concat_final = panda_concat_final.append(panda_processed_i1, ignore_index=True)

    panda_concat_final.to_csv('data/paper_authors.csv', header=None, index=None, sep='\t')

def step_10():
    paper_cits = dd.read_csv('data/paper_cits.csv', header=None, sep='\t',
                             names = ['paper_id', 'cits', 'years'],
                             dtype={'years': 'string'})
    paper_authors = dd.read_csv('data/paper_authors.csv', header=None, sep='\t',
                         names=['paper_id', 'authors'])
    paper_infos = dd.read_csv('data/papers_filtered_1_byPaperID.txt', header=None, sep='\t',
                        names=['paper_id', 'doi', 'year', 'citation_count'],
                        dtype={'citation_count': 'object', 'year': 'object'})[['paper_id', 'doi', 'year']]

    paper_cits = paper_cits.set_index('paper_id', sorted=True)
    paper_authors = paper_authors.set_index('paper_id', sorted=True)
    paper_infos = paper_infos.set_index('paper_id', sorted=True)
    join_result1 = paper_infos.merge(paper_authors, how='left', left_index=True, right_index=True)
    join_result2 = join_result1.merge(paper_cits, how='left', left_index=True, right_index=True)

    join_result2.to_csv('data/paper_complete_infos.csv', sep='\t', single_file=True, header=False)

if '__main__' == __name__:
#     step_1() # filter erros in papers.txt file to future open on dask
#     step_2()
#     step_3()
#     step_4()
    # step_5()
#     print('step 6')
#     step_6() # citações dos artigos
#     print('step 7')
#     step_7() # união dos arquivos de citações dos artigos
#     print('step 8')
#     step_8() # autores dos artigos
#     print('step 9')
#     step_9() # união das colaborações dos autores
#     print('step 10')
    step_10()

# 249 990 849 authors
# 249 265 454 papers


text:
from unittest import TestCase
from numpy import uint
from model.grid import Grid
from model.grid_pos import GridPos
from model.agent_state import AgentState
from controller.scheduler import Scheduler

class TestGrid(TestCase):
    """
    Author: Beil Benedikt
    """

    def test_is_fully_occupied_True(self):
        sut = self.__get_sut()
        sut.place_agents_on_the_field(uint(100), uint(100), 2523, 1, 0)
        self.assertTrue(sut.is_fully_occupied())

    def test_is_fully_occupied_False(self):
        sut = self.__get_sut()
        sut.place_agents_on_the_field(uint(10), uint(10), 325, 0.99, 0)
        self.assertFalse(sut.is_fully_occupied())

    def test_is_occupied(self):
        sut = self.__get_sut()
        sut.init_empty_grid(uint(100), uint(100))
        test_pos = GridPos(uint(0), uint(0))
        sut.spawn_agent(test_pos, AgentState.SUSCEPTIBLE)
        self.assertTrue(sut.is_occupied(test_pos))
        self.assertFalse(sut.is_occupied(GridPos(uint(0), uint(1))))
        self.assertFalse(sut.is_occupied(GridPos(uint(1), uint(0))))
        self.assertFalse(sut.is_occupied(GridPos(uint(1), uint(1))))
        self.assertFalse(sut.is_occupied(GridPos(uint(99), uint(99))))
        self.assertFalse(sut.is_occupied(GridPos(uint(99), uint(0))))
        self.assertFalse(sut.is_occupied(GridPos(uint(0), uint(99))))

    def __get_sut(self) -> Grid:
        return Grid(Scheduler())

reason: from numpy


text:
"""generate_bottleneck_features.py"""

import argparse
import json
import numpy as np
import os
import sys
from keras import backend
from keras.applications.vgg16 import VGG16
from keras.applications.vgg16 import preprocess_input
from keras.preprocessing import image
from keras.preprocessing.image import ImageDataGenerator

parser = argparse.ArgumentParser(description='Generate bottleneck features from'
                                ' APOD images.')

parser.add_argument('conf',
                    type=argparse.FileType('r'),
                    default=sys.stdin,
                    help='path to configuration file')

def main():
    """Saves bottleneck features for all images using the VGG16 pre-trained CNN model."""
    # Read configuration file
    conf = None
    args = parser.parse_args()
    with open(args.conf.name, 'r') as conf_file:
        conf = json.load(conf_file)

    if not conf:
        raise SystemExit('Invalid Configuration File')

    # Paths initialization
    src_dir = os.path.dirname(os.path.abspath(__file__))
    app_dir = os.path.dirname(src_dir)

    # Images directory
    im_dir = os.path.join(app_dir, conf['download_dir'])

    # Configurations for generating bottleneck features
    gen_features_conf = conf['bottleneck_features']

    # Image target size
    img_width = gen_features_conf['img_width']
    img_height = gen_features_conf['img_height']
    crop_size = (img_width, img_height)

    # Save features to
    features_dir = gen_features_conf['features_dir']
    features_file = gen_features_conf['features_file']
    features_path = os.path.join(app_dir, features_dir, features_file)

    # Arguments for ImageDataGenerator
    data_gen_args = dict(rescale=1./255)

    # Initialize image data generator
    image_datagen = ImageDataGenerator(**data_gen_args)

    # Get an image generator
    image_generator = image_datagen.flow_from_directory(
        im_dir,
        target_size=crop_size,
        class_mode=None,
        shuffle=False)

    # Load the pre-trained model
    model = VGG16(weights='imagenet', include_top=False)

    # Extract bottleneck features
    features = model.predict_generator(
        image_generator,
        verbose=gen_features_conf['predict_generator']['verbose'],
        use_multiprocessing=gen_features_conf['predict_generator']['use_multiprocessing'])

    # Save bottleneck features
    np.save(features_path, features)

if __name__ == "__main__":
    main()

max_stars_repo_name: dtcrout/apod-caption-generator

max_issues_repo_name: dtcrout/apod-caption-generator

max_forks_repo_name: dtcrout/apod-caption-generator