text:
import math

import cv2 as cv
import numpy as np
import torch
from PIL import Image
from scipy.stats import norm
from torchvision import transforms

from align_faces import get_reference_facial_points, warp_and_crop_face
from config import device
from config import image_h, image_w
from mtcnn.detector import detect_faces

data_transforms = {
    'train': transforms.Compose([
        transforms.RandomHorizontalFlip(),
        transforms.ToTensor(),
        transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])
        # transforms.Normalize([0.5, 0.5, 0.5], [0.5, 0.5, 0.5])
    ]),
    'val': transforms.Compose([
        transforms.ToTensor(),
        transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])
        # transforms.Normalize([0.5, 0.5, 0.5], [0.5, 0.5, 0.5])
    ]),
}
transformer = data_transforms['val']

checkpoint = 'BEST_checkpoint.tar'
print('loading model: {}...'.format(checkpoint))
checkpoint = torch.load(checkpoint)
model = checkpoint['model']
model = model.to(device)
model.eval()

# model params
threshold = 73.18799151798612
mu_0 = 89.6058
sigma_0 = 4.5451
mu_1 = 43.5357
sigma_1 = 8.83

class FaceNotFoundError(Exception):
    """Base class for other exceptions"""
    pass

def align_face(img_fn, facial5points):
    raw = cv.imread(img_fn, True)
    facial5points = np.reshape(facial5points, (2, 5))

    crop_size = (image_h, image_w)

    default_square = True
    inner_padding_factor = 0.25
    outer_padding = (0, 0)
    output_size = (image_h, image_w)

    # get the reference 5 landmarks position in the crop settings
    reference_5pts = get_reference_facial_points(
        output_size, inner_padding_factor, outer_padding, default_square)

    # dst_img = warp_and_crop_face(raw, facial5points)
    dst_img = warp_and_crop_face(raw, facial5points, reference_pts=reference_5pts, crop_size=crop_size)
    return dst_img

def get_central_face_attributes(full_path):
    try:
        img = Image.open(full_path).convert('RGB')
        bounding_boxes, landmarks = detect_faces(img)

        if len(landmarks) > 0:
            i = select_central_face(img.size, bounding_boxes)
            return True, [bounding_boxes[i]], [landmarks[i]]

    except KeyboardInterrupt:
        raise
    except:
        pass
    return False, None, None

def get_all_face_attributes(full_path):
    img = Image.open(full_path).convert('RGB')
    bounding_boxes, landmarks = detect_faces(img)
    return bounding_boxes, landmarks

def select_central_face(im_size, bounding_boxes):
    width, height = im_size
    nearest_index = -1
    nearest_distance = 100000
    for i, b in enumerate(bounding_boxes):
        x_box_center = (b[0] + b[2]) / 2
        y_box_center = (b[1] + b[3]) / 2
        x_img = width / 2
        y_img = height / 2
        distance = math.sqrt((x_box_center - x_img) ** 2 + (y_box_center - y_img) ** 2)
        if distance < nearest_distance:
            nearest_distance = distance
            nearest_index = i

    return nearest_index

def draw_bboxes(img, bounding_boxes, facial_landmarks=[]):
    for b in bounding_boxes:
        cv.rectangle(img, (int(b[0]), int(b[1])), (int(b[2]), int(b[3])), (255, 255, 255), 1)

    for p in facial_landmarks:
        for i in range(5):
            cv.circle(img, (int(p[i]), int(p[i + 5])), 1, (0, 255, 0), -1)

        break  # only first

    return img

def get_image(filename):
    has_face, bboxes, landmarks = get_central_face_attributes(filename)
    if not has_face:
        raise FaceNotFoundError(filename)

    img = align_face(filename, landmarks)
    img = transforms.ToPILImage()(img)
    img = transformer(img)
    img = img.to(device)

    print('drawing bboxes: {}'.format(filename))
    bboxes, landmarks = get_all_face_attributes(filename)
    pic = cv.imread(filename)
    pic = draw_bboxes(pic, bboxes, landmarks)
    cv.imwrite(filename, pic)

    return img

def resize(filename):
    img = cv.imread(filename)
    h, w = img.shape[:2]
    ratio_w = w / 1280
    ratio_h = h / 720
    if ratio_w > 1 or ratio_h > 1:
        ratio = max(ratio_w, ratio_h)
        new_w = int(w / ratio)
        new_h = int(h / ratio)
        img = cv.resize(img, (new_w, new_h))
        cv.imwrite(filename, img)

def compare(fn_0, fn_1):
    print('fn_0: ' + fn_0)
    print('fn_1: ' + fn_1)
    img0 = get_image(fn_0)
    img1 = get_image(fn_1)
    imgs = torch.zeros([2, 3, 112, 112], dtype=torch.float)
    imgs[0] = img0
    imgs[1] = img1

    with torch.no_grad():
        output = model(imgs)

        feature0 = output[0].cpu().numpy()
        feature1 = output[1].cpu().numpy()
        x0 = feature0 / np.linalg.norm(feature0)
        x1 = feature1 / np.linalg.norm(feature1)
        cosine = np.dot(x0, x1)
        theta = math.acos(cosine)
        theta = theta * 180 / math.pi

    print('theta: ' + str(theta))
    prob = get_prob(theta)
    print('prob: ' + str(prob))
    return prob, theta < threshold

def get_prob(theta):
    prob_0 = norm.pdf(theta, mu_0, sigma_0)
    prob_1 = norm.pdf(theta, mu_1, sigma_1)
    total = prob_0 + prob_1
    return prob_1 / total

def ensure_folder(folder):
    import os
    if not os.path.isdir(folder):
        os.mkdir(folder)

if __name__ == "__main__":
    compare('id_card.jpg', 'photo_1.jpg')
    compare('id_card.jpg', 'photo_2.jpg')
    compare('id_card.jpg', 'photo_3.jpg')
    compare('id_card.jpg', 'photo_4.jpg')

max_stars_repo_name: foamliu/Twins

max_issues_repo_name: foamliu/Twins

max_forks_repo_name: foamliu/Twins


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

# My imports
from __future__ import division
import numpy as np
import pandas as pd
import argparse
from utils import _update_par as updateBatch
from utils import _run_moog as runMoog
from utils import Readmoog
from interpolation import interpolator
import os

def solar_abundance(atom):
    '''Give atomic number and return solar abundance from Asplund et al. 2009

    Input
    -----
    atom : int
      The atomic number

    Output
    ------
    abundance : float
      The solar abundance of the atom
    '''
    if not isinstance(atom, int):
        raise ValueError('Atomic number need to be an integer')
    solar = [12.00, 10.93, 1.05, 1.38, 2.70, 8.43, 7.83, 8.69, 4.56, 7.93,
             6.24, 7.60, 6.45, 7.51, 5.41, 7.12, 5.50, 6.40, 5.03, 6.34,
             3.15, 4.95, 3.93, 5.64, 5.43, 7.47, 4.99, 6.22, 4.19, 4.56,
             3.04, 3.65, 2.30, 3.34, 2.54, 3.25, 2.52, 2.87, 2.21, 2.58,
             1.46, 1.88, -5.00, 1.75, 0.91, 1.57, 0.94, 1.71, 0.80, 2.04,
             1.01, 2.18, 1.55, 2.24, 1.08, 2.18, 1.10, 1.58, 0.72, 1.42,
             -5.00, 0.96, 0.52, 1.07, 0.30, 1.10, 0.48, 0.92, 0.10, 0.84,
             0.10, 0.85, -0.12, 0.85, 0.26, 1.40, 1.38, 1.62, 0.92, 1.17,
             0.90, 1.75, 0.65, -5.00, -5.00, -5.00, -5.00, -5.00, -5.00,
             0.02, -5.00, -0.54, -5.00, -5.00, -5.00]
    return solar[atom-1]

def recalSingleLine(line, params=None, version=2014, maxiter=40, driver='abfind'):
    '''Recalibrate a single line and return the new loggf

    Inputs
    ------
    line : list
      The line containing (wavelength, element, EP, loggf, EW) in that order
    params : list/tuple
      The parameters (Teff, logg, [Fe/H], vt)
    version : int
      The version of MOOG
    driver : str
      The MOOG driver to use (abfind or ewfind)

    Output
    ------
    loggf : float
      The new recalibrated loggf
    '''

    def moogAbund(loggf, ewdriver=False):
        line[3] = loggf
        np.savetxt('temporary.moog', line[:, np.newaxis].T, fmt=fmt, header=header)
        runMoog()
        if ewdriver:
            d = np.loadtxt('summary.out', skiprows=5, usecols=(6,))
            out = d-line[4]
        else:
            m = Readmoog(params=params, version=version)
            _, abund = m.elements()
            solar = solar_abundance(int(line[1]))
            out = round(abund[0] - solar, 3)
        return out

    ewdriver = True if driver == 'ewfind' else False

    fmt = ('%9.3f', '%10.1f', '%9.2f', '%9.3f', '%28.1f')
    header = 'Wavelength     ele    EP     loggf        EW'
    np.savetxt('temporary.moog', line[:, np.newaxis].T, fmt=fmt, header=header)
    loggf_old = line[3]
    a, b = loggf_old-5, loggf_old+5  # extreme values of loggf
    c = (a+b)/2
    for _ in range(maxiter):
        if c == 0:  # Don't evaluate at loggf = 0
            c += (abs(a) + abs(b)) / 10
        fa = moogAbund(a, ewdriver=ewdriver)
        fc = moogAbund(c, ewdriver=ewdriver)
        if fc == 0:
            return c
        elif fa*fc < 0:
            b = c
        else:
            a = c
        c = (a+b)/2
    return c

def _parser():
    parser = argparse.ArgumentParser(description='Recalibrate the loggf for a set of settings')
    parser.add_argument('input', help='Input linelist from rawLinelist folder')
    parser.add_argument('output', help='Name of output file (saved in rawLinelist)')
    parser.add_argument('-m', '--model', help='Model atmosphere', default='kurucz95', choices=['kurucz95', 'apogee_kurucz', 'marcs'])
    parser.add_argument('-v', '--moogversion', help='MOOG version', default=2014)
    parser.add_argument('-d', '--damping', help='Damping to be used in MOOG', default=1, choices=map(str, [1, 2]))
    parser.add_argument('-dr', '--driver', help='Which driver to use', default='abfind', choices=['abfind', 'ewfind'])
    parser.add_argument('-p', '--parameters', help='Atmospheric parameters, Teff, logg, [Fe/H], vt', nargs='+', default=None)
    args = parser.parse_args()
    return args

if __name__ == '__main__':  # pragma: no cover

    args = _parser()

    fname = args.input
    fout1 = 'rawLinelist/%s' % args.output
    fout2 = 'linelist/%s' % args.output.replace('.ares', '.moog')
    lines = pd.read_csv(fname, skiprows=2, delimiter=r'\s+',
                        names=['WL', 'num', 'EP', 'loggf', 'ele', 'EW'])

    if args.parameters is None:
        params = [5777, 4.44, 0.00, 1.00]
    else:
        params = map(float, args.parameters)
        params[0] = int(params[0])

    interpolator(params=params, atmtype=args.model, save=True)

    cols = ['WL', 'num', 'EP', 'loggf', 'EW']
    fmt2 = ('%9.3f', '%10.1f', '%9.2f', '%9.3f', '%28.1f')
    header1 = 'WL         num       E.P.     loggf         ele     EWsun\n'
    header1 += '-------    ----      ----     ------        ----    -----'
    header2 = 'Wavelength     ele       EP      loggf                          EW'
    x = lines[cols].values[0][:, np.newaxis].T
    np.savetxt('temporary.moog', x, fmt=fmt2, header=header2)

    options = {'driver': args.driver,
               'damping': args.damping}
    updateBatch(line_list='temporary.moog', **options)

    newloggf = np.zeros(lines.shape[0])
    for i, line in enumerate(lines[cols].values):
        print 'Wavelength: %.3f' % line[0]
        print 'Old loggf: %.3f' % line[3]
        zz = recalSingleLine(line, params=params, version=args.moogversion, driver=args.driver)
        zz = np.log10(zz) if zz > 0 else zz
        newloggf[i] = zz
        print 'New loggf: %.3f\n' % newloggf[i]

    lines['newloggf'] = pd.Series(newloggf)
    X = lines[['WL', 'num', 'EP', 'newloggf', 'ele', 'EW']]
    fmt1 = ('%7.2f', '%7.1f', '%9.2f', '%10.3f', '%10s', '%9.1f')
    print 'Saving results to: %s' % fout1
    np.savetxt(fout1, X, fmt=fmt1, header=header1, comments='')

    X = lines[['WL', 'num', 'EP', 'newloggf', 'EW']]
    print 'Saving results to: %s' % fout2
    np.savetxt(fout2, X, fmt=fmt2, header=header2)
    os.remove('temporary.moog')

ext: py

max_stars_repo_name: DanielAndreasen/pymoog

max_issues_repo_name: DanielAndreasen/pymoog

max_forks_repo_name: DanielAndreasen/pymoog

reason: import numpy


text:
import numpy as _np

def mi_to_m(d,out=None):
	"""Converts miles to meters.

	Examples
	--------
	>>> dist_meters = mi_to_m(dist_miles)

	Parameters
	----------
	d : array_like,scalar
		array / value in miles to convert to meters
	out : array_like, optional
		buffer for output values to go to.

	"""
	return _np.multiply(1609.344,d,out=out)

def ft_to_m(d,out=None):
	"""Converts feet to meters.

	Examples
	--------
	>>> dist_meters = mi_to_m(dist_feet)

	Parameters
	----------
	d : array_like,scalar
		array / value in feet to convert to meters
	out : array_like, optional
		buffer for output values to go to.

	"""
	return _np.multiply(0.3048,d,out=out)

def km_to_m(d,out=None):
	"""Converts kilometers to meters.

	Examples
	--------
	>>> dist_meters = mi_to_m(dist_kilometers)

	Parameters
	----------
	d : array_like,scalar
		array / value in milometers to convert to meters
	out : array_like, optional
		buffer for output values to go to.

	"""
	return _np.multiply(1000.0,d,out=out)

def m_to_mi(d,out=None):
	"""Converts meters to miles.

	Examples
	--------
	>>> dist_meiles = mi_to_m(dist_meters)

	Parameters
	----------
	d : array_like,scalar
		array / value in meters to convert to miles
	out : array_like, optional
		buffer for output values to go to.

	"""
	return _np.divide(d,1609.344,out=out)

def m_to_ft(d,out=None):
	"""Converts meters to feet.

	Examples
	--------
	>>> dist_feet = mi_to_m(dist_meters)

	Parameters
	----------
	d : array_like,scalar
		array / value in meters to convert to feet
	out : array_like, optional
		buffer for output values to go to.

	"""
	return _np.divide(d,0.3048,out=out)

def m_to_km(d,out=None):
	"""Converts meters to kilometers.

	Examples
	--------
	>>> dist_kilometers = mi_to_m(dist_meters)

	Parameters
	----------
	d : array_like,scalar
		array / value in meters to convert to kilometers
	out : array_like, optional
		buffer for output values to go to.

	"""
	return _np.divide(d,1000.0,out=out)


text:
import os
import numpy as np
import pandas as pd
import mne
from mne.preprocessing import ICA, create_ecg_epochs, create_eog_epochs
from autoreject import AutoReject
from scipy.io import loadmat, savemat
#from brainpipe import feature
from src.saflow_params import BIDS_PATH, LOGS_DIR
from mne.io import read_raw_fif, read_raw_ctf
#from hytools.meg_utils import get_ch_pos
from src.utils import get_SAflow_bids
from src.behav import find_logfile, get_VTC_from_file
import random
from matplotlib.pyplot import close
from mne.time_frequency import psd_multitaper, psd_welch
import pickle

def find_rawfile(subj, bloc, BIDS_PATH):
    filepath = '/sub-{}/ses-recording/meg/'.format(subj)
    files = os.listdir(BIDS_PATH + filepath)
    for file in files:
        if file[-8] == bloc:
            filename = file
    return filepath, filename

def saflow_preproc(filepath, savepath, reportpath, ica=True):

    report = mne.Report(verbose=True)
    raw_data = read_raw_ctf(filepath, preload=True)
    raw_data = raw_data.apply_gradient_compensation(grade=3) #required for source reconstruction
    picks = mne.pick_types(raw_data.info, meg=True, eog=True, exclude='bads')
    fig = raw_data.plot(show=False);
    report.add_figs_to_section(fig, captions='Time series', section='Raw data')
    close(fig)
    fig = raw_data.plot_psd(average=False, picks=picks, show=False);
    report.add_figs_to_section(fig, captions='PSD', section='Raw data')
    close(fig)

    ## Filtering
    high_cutoff = 200
    low_cutoff = 0.5
    raw_data.filter(low_cutoff, high_cutoff, fir_design="firwin")
    raw_data.notch_filter(np.arange(60, high_cutoff+1, 60), picks=picks, filter_length='auto',phase='zero', fir_design="firwin")
    fig = raw_data.plot_psd(average=False, picks=picks, fmax=120, show=False);
    report.add_figs_to_section(fig, captions='PSD', section='Filtered data')
    close(fig)
    if ica == False :
        report.save(reportpath, open_browser=False, overwrite=True);
        raw_data.save(savepath, overwrite=True)
        del report
        del raw_data
        del fig

    elif ica == True :
        ## ICA
        ica = ICA(n_components=20, random_state=0).fit(raw_data, decim=3)
        fig = ica.plot_sources(raw_data, show=False);
        report.add_figs_to_section(fig, captions='Independent Components', section='ICA')
        close(fig)

        ## FIND ECG COMPONENTS
        ecg_threshold = 0.50
        ecg_epochs = create_ecg_epochs(raw_data, ch_name='EEG059')
        ecg_inds, ecg_scores = ica.find_bads_ecg(ecg_epochs, ch_name='EEG059', method='ctps', threshold=ecg_threshold)
        fig = ica.plot_scores(ecg_scores, ecg_inds, show=False);
        report.add_figs_to_section(fig, captions='Correlation with ECG (EEG059)', section='ICA - ECG')
        close(fig)
        fig = list()
        try:
            fig = ica.plot_properties(ecg_epochs, picks=ecg_inds, image_args={'sigma': 1.}, show=False);
            for i, figure in enumerate(fig):
                report.add_figs_to_section(figure, captions='Detected component ' + str(i), section='ICA - ECG')
                close(figure)
        except:
            print('No component to remove')

        ## FIND EOG COMPONENTS
        eog_threshold = 4
        eog_epochs = create_eog_epochs(raw_data, ch_name='EEG057')
        eog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name='EEG057', threshold=eog_threshold)
        #TODO : if eog_inds == [] then eog_inds = [index(max(abs(eog_scores)))]
        fig = ica.plot_scores(eog_scores, eog_inds, show=False);
        report.add_figs_to_section(fig, captions='Correlation with EOG (EEG057)', section='ICA - EOG')
        close(fig)
        fig = list()
        try:
            fig = ica.plot_properties(eog_epochs, picks=eog_inds, image_args={'sigma': 1.}, show=False);
            for i, figure in enumerate(fig):
                report.add_figs_to_section(figure, captions='Detected component ' + str(i), section='ICA - EOG')
                close(figure)
        except:
            print('No component to remove')

        ## EXCLUDE COMPONENTS
        ica.exclude = ecg_inds
        ica.apply(raw_data)
        ica.exclude = eog_inds
        ica.apply(raw_data)
        fig = raw_data.plot(show=False); # Plot the clean signal.
        report.add_figs_to_section(fig, captions='After filtering + ICA', section='Raw data')
        close(fig)
        ## SAVE PREPROCESSED FILE
        report.save(reportpath, open_browser=False, overwrite=True);
        raw_data.save(savepath, overwrite=True)
        del ica
        del report
        del raw_data
        del fig

def segment_files(bids_filepath, tmin=0, tmax=0.8):
    raw = read_raw_fif(bids_filepath, preload=True)
    picks = mne.pick_types(raw.info, meg=True, ref_meg=True, eeg=False, eog=False, stim=False)
    ### Set some constants for epoching
    baseline = None #(None, -0.05)
    #reject = {'mag': 4e-12}
    try:
        events = mne.find_events(raw, min_duration=1/raw.info['sfreq'], verbose=False)
    except ValueError:
        events = mne.find_events(raw, min_duration=2/raw.info['sfreq'], verbose=False)
    event_id = {'Freq': 21, 'Rare': 31}
    epochs = mne.Epochs(raw, events=events, event_id=event_id, tmin=tmin,
                    tmax=tmax, baseline=baseline, reject=None, picks=picks, preload=True)
    ar = AutoReject(n_jobs=6)
    epochs_clean, autoreject_log = ar.fit_transform(epochs, return_log=True)
    return epochs_clean, autoreject_log

def remove_errors(logfile, events):
    '''
    Takes as input the raw events, including both stimuli and responses.
    Outputs an event vector containing only the correct trials, and indices of their position in the list
    of all stims
    '''

    # load RT vector
    data = loadmat(logfile)
    df_response = pd.DataFrame(data['response'])
    RT_array= np.asarray(df_response.loc[:,4])

    # get events vector with only stimuli events
    new_events = []
    for ev in events:
        if ev[2] != 99:
            new_events.append(ev)
    events = np.array(new_events)

    # check for each trial if it has a response or not
    events_noerr = []

    events_comerr = []
    events_omerr = []

    for idx, event in enumerate(events):
        if event[2] == 21:
            if RT_array[idx] != 0:
                events_noerr.append(event)
            else:
                events_omerr.append(event)
        if event[2] == 31:
            if RT_array[idx] == 0:
                events_noerr.append(event)
            else:
                events_comerr.append(event)
    events_noerr = np.array(events_noerr)
    events_comerr = np.array(events_comerr)
    events_omerr = np.array(events_omerr)

    return events_noerr, events_comerr, events_omerr

def trim_events(events_noerr, events_artrej):
    '''
    This function compares the events vectors of correct epochs (events_noerr)
    and of kept epochs after auto-reject (events_artrej).
    Returns a list of intersecting epochs, + their idx in the clean epochs vector
    '''
    events_trimmed = []
    idx_trimmed = []
    for idx, event in enumerate(events_artrej):
        if event[0] in events_noerr[:,0]:
            events_trimmed.append(event)
            idx_trimmed.append(idx)

    events_trimmed = np.array(events_trimmed)
    idx_trimmed = np.array(idx_trimmed)

    print('N events in clean epochs : {}'.format(len(events_artrej)))
    print('N events in correct epochs : {}'.format(len(events_noerr)))
    print('N events in intersection : {}'.format(len(idx_trimmed)))
    return events_trimmed, idx_trimmed

def trim_INOUT_idx(INidx, OUTidx, events_trimmed, events):
    '''
    With INidx_trimmed refering to indices in events_artrej
    '''
    # get events vector with only stimuli events
    new_events = []
    for ev in events:
        if ev[2] != 99:
            new_events.append(ev)
    events = np.array(new_events)

    INidx_trimmed = []
    OUTidx_trimmed = []
    # compare trimmed events with all_events, and store corresponding indices
    for idx, ev in enumerate(events):
        for idx_trim, ev_trim in enumerate(events_trimmed):
            if ev[0] == ev_trim[0]:
                if idx in INidx:
                    INidx_trimmed.append(idx_trim)
                if idx in OUTidx:
                    OUTidx_trimmed.append(idx_trim)
    INidx_trimmed = np.array(INidx_trimmed)
    OUTidx_trimmed = np.array(OUTidx_trimmed)

    return INidx_trimmed, OUTidx_trimmed

def get_odd_epochs(BIDS_PATH, LOGS_DIR, subj, bloc, stage='-epo'):
    '''
    Returns an array of indices of Freqs and Rares epochs. Retains only clean epochs.
    '''
    ### Get events after artifact rejection have been performed
    epo_path, epo_filename = get_SAflow_bids(BIDS_PATH, subj, bloc, stage=stage, cond=None)
    events_artrej = mne.read_events(epo_filename, verbose=False) # get events from the epochs file (so no resp event)

    ### Get original events from the raw file, to compare them to the events left in the epochs file
    events_fname, events_fpath = get_SAflow_bids(BIDS_PATH, subj, bloc, stage='preproc_raw', cond=None)
    raw = read_raw_fif(events_fpath, preload=False, verbose=False)#, min_duration=2/epochs.info['sfreq'])
    try:
        events = mne.find_events(raw, min_duration=1/raw.info['sfreq'], verbose=False)
    except ValueError:
        events = mne.find_events(raw, min_duration=2/raw.info['sfreq'], verbose=False)

    # Get the list of hits/miss events
    log_file = LOGS_DIR + find_logfile(subj,bloc,os.listdir(LOGS_DIR))
    events_noerr, events_comerr, events_omerr = remove_errors(log_file, events)

    # Keep only events that are clean, and split them by condition
    # Start with correct events
    events_noerr_trimmed, idx_noerr_trimmed = trim_events(events_noerr, events_artrej)
    freqs_hits_idx = np.array([idx_noerr_trimmed[i] for i, x in enumerate(events_noerr_trimmed) if x[2] == 21])
    rares_hits_idx = np.array([idx_noerr_trimmed[i] for i, x in enumerate(events_noerr_trimmed) if x[2] == 31])

    # Then commission errors
    if events_comerr.size > 0:
        events_comerr_trimmed, idx_comerr_trimmed = trim_events(events_comerr, events_artrej)
        rares_miss_idx = np.array(idx_comerr_trimmed)
    else:
        rares_miss_idx = np.array([])
    # And finally ommission errors
    if events_omerr.size > 0:
        events_omerr_trimmed, idx_omerr_trimmed = trim_events(events_omerr, events_artrej)
        freqs_miss_idx = np.array(idx_omerr_trimmed)
    else:
        freqs_miss_idx = np.array([])

    return freqs_hits_idx, freqs_miss_idx, rares_hits_idx, rares_miss_idx

def get_VTC_epochs(BIDS_PATH, LOGS_DIR, subj, bloc, stage='-epo', lobound=None, hibound=None, save_epochs=False, filt_order=3, filt_cutoff=0.1):
    '''
    This functions allows to use the logfile to split the epochs obtained in the epo.fif file.
    It works by comparing the timestamps of IN and OUT events to the timestamps in the epo file events
    It returns IN and OUT indices that are to be used in the split_PSD_data function

    '''
    ### Get events after artifact rejection have been performed
    epo_path, epo_filename = get_SAflow_bids(BIDS_PATH, subj, bloc, stage=stage, cond=None)
    events_artrej = mne.read_events(epo_filename, verbose=False) # get events from the epochs file (so no resp event)

    ### Find logfile to extract VTC
    log_file = LOGS_DIR + find_logfile(subj,bloc,os.listdir(LOGS_DIR))
    VTC, INbounds, OUTbounds, INidx, OUTidx, RT_array = get_VTC_from_file(log_file, lobound=lobound, hibound=hibound, filt=True, filt_order=filt_order, filt_cutoff=filt_cutoff)

    ### Get original events and split them using the VTC
    events_fname, events_fpath = get_SAflow_bids(BIDS_PATH, subj, bloc, stage='preproc_raw', cond=None)
    raw = read_raw_fif(events_fpath, preload=False, verbose=False)#, min_duration=2/epochs.info['sfreq'])
    try:
        events = mne.find_events(raw, min_duration=1/raw.info['sfreq'], verbose=False)
    except ValueError:
        events = mne.find_events(raw, min_duration=2/raw.info['sfreq'], verbose=False)

    events_noerr, events_comerr, events_omerr = remove_errors(log_file, events)
    # Keep only events that are correct and clean
    events_trimmed, idx_trimmed = trim_events(events_noerr, events_artrej)
    # Write INidx and OUTidx as indices of clean events
    INidx, OUTidx = trim_INOUT_idx(INidx, OUTidx, events_trimmed, events)

    VTC_epo = np.array([VTC[idx] for idx in idx_trimmed])

    return INidx, OUTidx, VTC_epo, idx_trimmed

def compute_PSD(epochs, freqlist=None, method='multitaper'):
    epochs_psds = []
    if freqlist == None:
        freqlist = [ [4, 8], [8, 12], [12, 20], [20, 30], [30, 60], [60, 90], [90, 120] ]
    # Compute PSD
    if method == 'multitaper':
        psds, freqs = psd_multitaper(epochs, fmin=min(min(freqlist)), fmax=max(max(freqlist)), n_jobs=1)
    if method == 'pwelch':
        psds, freqs = psd_welch(epochs, average='median', fmin=min(min(freqlist)), fmax=max(max(freqlist)), n_jobs=1)
    psds = 10. * np.log10(psds) # Convert power to dB scale.
    # Average in freq bands
    for low, high in freqlist:
        freq_idx = [i for i, x in enumerate(freqs) if x >=low and x <= high]
        psd = np.mean(psds[:,:,freq_idx], axis=2)
        epochs_psds.append(psd)
    epochs_psds = np.array(epochs_psds).swapaxes(2,0).swapaxes(1,0)
    return epochs_psds

def compute_PSD_hilbert(raw, ARlog, freqlist=None, tmin=0, tmax=0.8):
    epochs_envelopes = []
    if freqlist == None:
        freqlist = [ [4, 8], [8, 12], [12, 20], [20, 30], [30, 60], [60, 90], [90, 120] ]
    for low, high in freqlist:
        # Filter continuous data
        data = raw.copy().filter(low, high)
        # Segment them
        picks = mne.pick_types(raw.info, meg=True, ref_meg=True, eeg=False, eog=False, stim=False)
        try:
            events = mne.find_events(raw, min_duration=1/raw.info['sfreq'], verbose=False)
        except ValueError:
            events = mne.find_events(raw, min_duration=2/raw.info['sfreq'], verbose=False)
        event_id = {'Freq': 21, 'Rare': 31}
        epochs = mne.Epochs(raw, events=events, event_id=event_id, tmin=tmin,
                        tmax=tmax, baseline=None, reject=None, picks=picks, preload=True)
        # Drop bad epochs
        epochs.drop(ARlog.bad_epochs)
        # Apply Hilbert transform to obtain envelope
        hilbert = epochs.apply_hilbert(envelope=True)
        hilbert_pow = hilbert.copy()
        hilbert_pow._data = hilbert._data**2
        epochs_envelopes.append(hilbert_pow)
        del hilbert_pow
        del hilbert
        del data
        del epochs
    return epochs_envelopes

def compute_TFR(epochs, baseline=False):
    decim = 3
    fmin = 2
    fmax = 150
    n_bins=30
    freqs = np.logspace(*np.log10([fmin, fmax]), num=n_bins)
    #freqs = np.arange(2,15,1)  # define frequencies of interest
    n_cycles = freqs / freqs[0]
    zero_mean = False
    this_tfr = mne.time_frequency.tfr_morlet(epochs, freqs, n_cycles=n_cycles,
                      decim=decim, average=False, zero_mean=zero_mean,
                      return_itc=False)
    if baseline:
        this_tfr.apply_baseline(mode='ratio', baseline=(None, -0.05))
    #this_power = this_tfr.data[:, :, :, :]  # we only have one channel.
    return this_tfr

def load_PSD_data(BIDS_PATH, SUBJ_LIST, BLOCS_LIST, time_avg=True, stage='PSD'):
    '''
    Returns a list containing n_subj lists of n_blocs matrices of shape n_freqs X n_channels X n_trials
    '''
    PSD_alldata = []
    for subj in SUBJ_LIST:
        all_subj = [] ## all the data of one subject
        for run in BLOCS_LIST:
            SAflow_bidsname, SAflow_bidspath = get_SAflow_bids(BIDS_PATH, subj, run, stage=stage, cond=None)
            mat = loadmat(SAflow_bidspath)['PSD']
            if time_avg == True:
                mat = np.mean(mat, axis=2) # average PSDs in time across epochs
            all_subj.append(mat)
        PSD_alldata.append(all_subj)
    return PSD_alldata

def load_VTC_data(BIDS_PATH, LOGS_DIR, SUBJ_LIST, BLOCS_LIST):
    VTC_alldata = []
    for subj in SUBJ_LIST:
        all_subj = [] ## all the data of one subject
        for run in BLOCS_LIST:
            # get events from epochs file
            epo_path, epo_filename = get_SAflow_bids(BIDS_PATH, subj, run, 'epo', cond=None)
            events_epoched = mne.read_events(epo_filename, verbose=False) # get events from the epochs file (so no resp event)
            # get events from original file (only 599 events)
            events_fname, events_fpath = get_SAflow_bids(BIDS_PATH, subj, run, stage='preproc_raw', cond=None)
            raw = read_raw_fif(events_fpath, preload=False, verbose=False)#, min_duration=2/epochs.info['sfreq'])
            all_events = mne.find_events(raw, min_duration=2/raw.info['sfreq'], verbose=False)
            stim_idx = []
            for i in range(len(all_events)):
                 if all_events[i,2] in [21,31]:
                     stim_idx.append(i)
            all_events = all_events[stim_idx]
            # compute VTC for all trials
            log_file = find_logfile(subj,run,os.listdir(LOGS_DIR))
            VTC, INbounds, OUTbounds, INzone, OUTzone = get_VTC_from_file(LOGS_DIR + log_file, lobound=None, hibound=None)
            epochs_VTC = []
            for event_time in events_epoched[:,0]:
                idx = list(all_events[:,0]).index(event_time)
                epochs_VTC.append(VTC[idx])
            all_subj.append(np.array(epochs_VTC))
        VTC_alldata.append(all_subj)
    return VTC_alldata

def split_TFR(BIDS_PATH, subj, bloc, by='VTC', lobound=None, hibound=None, stage='1600TFR', filt_order=3, filt_cutoff=0.05):
    if by == 'VTC':
        event_id = {'IN': 1, 'OUT': 0}
        INidx, OUTidx, VTC_epochs, idx_trimmed = get_VTC_epochs(LOGS_DIR, subj, bloc, lobound=lobound, hibound=hibound, stage=stage[:-3]+'epo', save_epochs=False, filt_order=filt_order, filt_cutoff=filt_cutoff)
        epo_path, epo_filename = get_SAflow_bids(BIDS_PATH, subj, bloc, stage=stage[:-3]+'epo', cond=None)
        epo_events = mne.read_events(epo_filename, verbose=False) # get events from the epochs file (so no resp event)
        TFR_path, TFR_filename = get_SAflow_bids(BIDS_PATH, subj, bloc, stage=stage, cond=None)
        TFR = mne.time_frequency.read_tfrs(TFR_filename)
        for i, event in enumerate(epo_events):
            if i in INidx:
                TFR[0].events[i,2] = 1
            if i in OUTidx:
                TFR[0].events[i,2] = 0
        TFR[0].event_id = event_id
    return TFR

def split_PSD_data(BIDS_PATH, SUBJ_LIST, BLOCS_LIST, by='VTC', lobound=None, hibound=None, stage='PSD', filt_order=3, filt_cutoff=0.1):
    '''
    This func splits the PSD data into two conditions. It returns a list of 2 (cond1 and cond2), each containing a list of n_subject matrices of shape n_freqs X n_channels X n_trials
    '''
    PSD_alldata = load_PSD_data(BIDS_PATH, SUBJ_LIST, BLOCS_LIST, time_avg=True, stage=stage)
    PSD_cond1 = []
    PSD_cond2 = []
    for subj_idx, subj in enumerate(SUBJ_LIST):
        subj_cond1 = []
        subj_cond2 = []
        for bloc_idx, bloc in enumerate(BLOCS_LIST):
            print('Splitting sub-{}_run-{}'.format(subj, bloc))

            # Obtain indices of the two conditions
            if by == 'VTC':
                INidx, OUTidx, VTC_epochs, idx_trimmed = get_VTC_epochs(LOGS_DIR, subj, bloc, lobound=lobound, hibound=hibound, save_epochs=False, filt_order=filt_order, filt_cutoff=filt_cutoff)
                cond1_idx = INidx
                cond2_idx = OUTidx
            if by == 'odd':
                # Get indices of freq and rare events
                ev_fname, ev_fpath = get_SAflow_bids(BIDS_PATH, subj, bloc, stage='epo')
                events_artrej = mne.read_events(ev_fpath)
                log_file = LOGS_DIR + find_logfile(subj,bloc,os.listdir(LOGS_DIR))
                events_fname, events_fpath = get_SAflow_bids(BIDS_PATH, subj, bloc, stage='preproc_raw', cond=None)
                raw = read_raw_fif(events_fpath, preload=False, verbose=False)#, min_duration=2/epochs.info['sfreq'])
                try:
                    events = mne.find_events(raw, min_duration=1/raw.info['sfreq'], verbose=False)
                except ValueError:
                    events = mne.find_events(raw, min_duration=2/raw.info['sfreq'], verbose=False)

                events_noerr, events_comerr, events_omerr = remove_errors(log_file, events)
                events_trimmed, idx_trimmed = trim_events(events_noerr, events_artrej)
                cond1_idx = []
                cond2_idx = []
                for idx, ev in enumerate(events_trimmed):
                    if ev[2] == 21: # Frequent events
                        cond1_idx.append(idx)
                    if ev[2] == 31:
                        cond2_idx.append(idx)
                cond1_idx = np.array(cond1_idx)
                cond2_idx = np.array(cond2_idx)
                # Add this to keep the same number of trials in both conditions
                random.seed(0)
                cond1_idx = random.choices(cond1_idx, k=len(cond2_idx))
                print('N trials retained for each condition : {}'.format(len(cond2_idx)))
            # Pick the data of each condition
            if bloc_idx == 0: # if first bloc, init ndarray size using the first matrix
                subj_cond1 = PSD_alldata[subj_idx][bloc_idx][:,:,cond1_idx]
                subj_cond2 = PSD_alldata[subj_idx][bloc_idx][:,:,cond2_idx]
            else: # if not first bloc, just concatenate along the trials dimension
                subj_cond1 = np.concatenate((subj_cond1, PSD_alldata[subj_idx][bloc_idx][:,:,cond1_idx]), axis=2)
                subj_cond2 = np.concatenate((subj_cond2, PSD_alldata[subj_idx][bloc_idx][:,:,cond2_idx]), axis=2)
        PSD_cond1.append(subj_cond1)
        PSD_cond2.append(subj_cond2)
    splitted_PSD = [PSD_cond1, PSD_cond2]
    return splitted_PSD

def split_trials(BIDS_PATH, LOGS_DIR, subj, run, stage='PSD', by='VTC',
                    keep_errors=False, equalize=False,
                    lobound=None, hibound=None,
                    filt_order=3, filt_cutoff=0.1,
                    freq_names=None, oddball='hits'):
    '''
    If stage is 'env', freq_names must be specified.
    '''
    # Split epochs indices
    if by == 'VTC':
        INidx, OUTidx, VTC_epochs, idx_trimmed = get_VTC_epochs(BIDS_PATH, LOGS_DIR, subj, run, lobound=lobound, hibound=hibound, save_epochs=False, filt_order=filt_order, filt_cutoff=filt_cutoff)
        condA_idx = INidx
        condB_idx = OUTidx
        print('{} IN epochs'.format(len(INidx)))
        print('{} OUT epochs'.format(len(OUTidx)))
    elif by == 'odd':
        freqs_hits_idx, freqs_miss_idx, rares_hits_idx, rares_miss_idx = get_odd_epochs(BIDS_PATH, LOGS_DIR, subj, run, stage='-epo')
        if oddball == 'all':
            condA_idx = np.sort(np.concatenate((freqs_hits_idx, freqs_miss_idx)))
            condB_idx = np.sort(np.concatenate((freqs_hits_idx, rares_miss_idx)))
        elif oddball == 'hits':
            condA_idx = freqs_hits_idx
            condB_idx = rares_hits_idx
        elif oddball == 'miss':
            condA_idx = freqs_miss_idx
            condB_idx = rares_miss_idx
        print('{} Freqs epochs'.format(len(condA_idx)))
        print('{} Rare epochs'.format(len(condB_idx)))
    elif by == 'resp':
        freqs_hits_idx, freqs_miss_idx, rares_hits_idx, rares_miss_idx = get_odd_epochs(BIDS_PATH, LOGS_DIR, subj, run, stage='-epo')
        condA_idx = np.sort(np.concatenate((freqs_hits_idx, rares_miss_idx)))
        condB_idx = np.sort(np.concatenate((freqs_miss_idx, rares_hits_idx)))
        print('{} Resp epochs'.format(len(condA_idx)))
        print('{} NoResp epochs'.format(len(condB_idx)))
        condA_idx = [int(x) for x in condA_idx]
        condB_idx = [int(x) for x in condB_idx]
    # Load epochs
    if stage == 'PSD':
        fname, fpath = get_SAflow_bids(BIDS_PATH, subj, run, stage, cond=None)
        with open(fpath, 'rb') as f:
            data = pickle.load(f)
        condA = data[condA_idx]
        condB = data[condB_idx]
    elif 'env' in stage:
        condA = []
        condB = []
        for freq in freq_names:
            fname, fpath = get_SAflow_bids(BIDS_PATH, subj, run, stage, cond=freq)
            data = mne.read_epochs(fpath)
            condA.append(data[condA_idx])
            condB.append(data[condB_idx])
    return condA, condB

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
"""
.. Authors
    Novimir Pablant <npablant@pppl.gov>
    James Kring <jdk0026@tigermail.auburn.edu>
    Yevgeniy Yakusevich <eugenethree@gmail.com>
"""

import numpy as np
from scipy.stats import cauchy
import scipy.constants as const

from xicsrt.util import profiler
from xicsrt.tools import xicsrt_voigt
from xicsrt.tools import xicsrt_spread
from xicsrt.tools.xicsrt_doc import dochelper
from xicsrt.objects._RayArray import RayArray
from xicsrt.objects._GeometryObject import GeometryObject

@dochelper
class XicsrtSourceGeneric(GeometryObject):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.filter_objects = []

    def default_config(self):
        """
        xsize
          The size of this element along the xaxis direction.

        ysize
          The size of this element along the yaxis direction.

        zsize
          The size of this element along the zaxis direction.

        angular_dist : string ('isotropic')
          | The type of angular distribution to use for the emitted rays.
          | Available distributions (default is 'isotropic'):
          |
          | isotropic
          |   Isotropic emission (uniform spherical) emitted in a cone (circular
          |   cross-section) with a half-angle of 'spread'. The axis of the
          |   emission cone is aligned along the z-axis. 'spread' must be a
          |   single value (scalar or 1-element array).
          | isotropic_xy
          |   Isotropic emission (uniform spherical) emitted in a truncated-cone
          |   (rectangular cross-section) with different x and y half-angles.
          |   'spread' can contain either 1, 2 or 4 values:
          |   s or [s]
          |     A single value that will be used for both the x and y directions.
          |   [x, y]
          |     Two values values that will be used for the x and y directions.
          |   [xmin, xmax, ymin, ymax]
          |     For values that define the asymmetric exent in x and y directions.
          |     Example: [-0.1, 0.1, -0.5, 0.5]
          | flat
          |   Flat emission (uniform planar) emitted in a cone (circular cross-
          |   section) with a half-angle of 'spread'.
          | flat_xy
          |   Flat emission (uniform planar) emitted in a truncated-cone
          |   (rectangular cross-section) with different x and y half-angles.
          |   'spread' can contain either 1, 2 or 4 values, see above.
          | gaussian
          |   Emission with angles away from the z-axis having a Gaussian
          |   distribution (circular cross-section). The 'spread' defines the
          |   Half-Width-at-Half-Max (HWHM) of the distribution. 'spread' must
          |   be a single value (scalar or 1-element array).
          | gaussian_flat
          |   !! Not implemented !!
          |   Cross-section of emission (intersection with constant-z plane) will
          |   have a Gaussian distribution.

        spread : float or array (np.pi) [radians]
          The angular spread for the emission cone. The spread defines the
          half-angle of the emission cone. See 'angular_dist' for detailed
          documentation.

        intensity : int or float
          The number of rays for this source to emit. This should be an
          integer value unless `use_poisson = True`.

          Note: If filters are attached, this will be the number of rays
          emitted before filtering.

        use_poisson : bool (False)
          If `True` the `intenisty` will be treated as the expected value for
          a Poisson distribution and the number of rays will be randomly
          picked from a Poisson distribution. This is setting is typically
          only used internally for Plasma sources.

        wavelength_dist : str ('voigt')
          The type of wavelength distribution for this source.
          Possible values are: 'voigt', 'uniform', 'monochrome'.

          Note: A monochrome distribution can also be achieved by using a
          'voigt' distribution with zero linewidth and temperature.

        wavelength : float (1.0) [angstroms]
          Only used if `wavelength_dist = "monochrome" or "voigt"`
          Central wavelength of the distribution, in Angstroms.

        wavelength_range: tuple [angstroms]
          Only used if `wavelength_dist = "uniform"`
          The wavelength range of the distribution, in Angstroms.
          Must be a 2 element tuple, list or array: (min, max).

        linewidth : float (0.0) [1/s]
          Only used if `wavelength_dist = "voigt"`
          The natural width of the emission line.
          This will control the Lorentzian contribution to the the overall Voigt
          profile. If linewidth == 0, the resulting wavelength distribution will
          be gaussian.

          To convert from a fwhm in [eV]:
          linewidth = 2*pi*e/h*fwhm_ev

          To translate from linewidth to gamma in the Voigt equation:
          gamma = linewidth * wavelength**2 / (4*pi*c*1e10)

        mass_number : float (1.0) [au]
          Only used if `wavelength_dist = "voigt"`
          The mass of the emitting atom in atomic units (au). This mass in used
          to convert temperature into line width. See temperature option.

        temperature : float (0.0) [eV]
          Only used if `wavelength_dist = "voigt"`
          The temperature of the emission line.
          This will control the Gaussian contribution to the overall Voigt
          profile. If temperature == 0, the resulting wavelength distribution
          will be Lorentzian.

          To translate from temperature to sigma in the Voigt equation:
          sigma = np.sqrt(temperature/mass_number/amu_kg/c**2*ev_J)*wavelength

        velocity
          No documentation yet. Please help improve XICSRT!

        filters
          No documentation yet. Please help improve XICSRT!

        """
        config = super().default_config()

        config['xsize'] = 0.0
        config['ysize'] = 0.0
        config['zsize'] = 0.0

        config['intensity']        = 0.0
        config['use_poisson']      = False
        config['angular_dist']      = 'isotropic'
        config['spread']           = np.pi

        # Possible values: 'monochrome', 'voigt', 'uniform
        config['wavelength_dist']  = 'voigt'

        # Only used for wavelength_dist = 'voigt' or 'monochrome'
        config['wavelength']       = 1.0

        # Only used for wavelength_dist = 'voigt'
        config['mass_number']      = 1.0
        config['linewidth']        = 0.0
        config['temperature']      = 0.0
        config['velocity']         = np.array([0.0, 0.0, 0.0])

        # Only used for wavelength_dist = 'uniform'
        config['wavelength_range'] = np.array([0.0, 0.0])

        config['filters']          = []

        return config

    def initialize(self):
        super().initialize()

        if self.param['use_poisson']:
            self.param['intensity'] = np.random.poisson(self.param['intensity'])
        else:
            if self.param['intensity'] < 1:
                raise ValueError('intensity of less than one encountered. Turn on poisson statistics.')
        self.param['intensity'] = int(self.param['intensity'])

    def generate_rays(self):
        rays = RayArray()
        profiler.start('generate_rays')

        profiler.start('generate_origin')
        rays['origin'] = self.generate_origin()
        profiler.stop('generate_origin')

        profiler.start('generate_direction')
        rays['direction'] = self.generate_direction(rays['origin'])
        profiler.stop('generate_direction')

        profiler.start('generate_wavelength')
        rays['wavelength'] = self.generate_wavelength(rays['direction'])
        profiler.stop('generate_wavelength')

        profiler.start('generate_weight')
        rays['weight'] = self.generate_weight()
        profiler.stop('generate_weight')

        profiler.start('generate_mask')
        rays['mask'] = self.generate_mask()
        profiler.stop('generate_mask')

        profiler.start('filter_rays')
        rays = self.ray_filter(rays)
        profiler.stop('filter_rays')

        profiler.stop('generate_rays')
        return rays

    def generate_origin(self):
        # generic origin for isotropic rays
        x_offset = np.random.uniform(-1 * self.param['xsize']/2 ,  self.param['xsize']/2, self.param['intensity'])
        y_offset = np.random.uniform(-1 * self.param['ysize']/2, self.param['ysize']/2, self.param['intensity'])
        z_offset = np.random.uniform(-1 * self.param['zsize']/2 ,  self.param['zsize']/2, self.param['intensity'])

        origin = (self.origin
                  + np.einsum('i,j', x_offset, self.xaxis)
                  + np.einsum('i,j', y_offset, self.yaxis)
                  + np.einsum('i,j', z_offset, self.zaxis))
        return origin

    def generate_direction(self, origin):
        normal = self.make_normal()
        D = self.random_direction(normal)
        return D

    def make_normal(self):
        array = np.empty((self.param['intensity'], 3))
        array[:] = self.param['zaxis']
        normal = array / np.linalg.norm(array, axis=1)[:, np.newaxis]
        return normal

    def random_direction(self, normal):

        spread = self.param['spread']
        dir_local  = xicsrt_spread.vector_distribution(
            spread,
            self.param['intensity'],
            name=self.param['angular_dist'],
            )

        # Generate some basis vectors that are perpendicular
        # to the normal. The orientation does not matter here.
        o_1  = np.cross(normal, np.array([0,0,1])) + np.cross(normal, np.array([0,1,0]))
        o_1 /=  np.linalg.norm(o_1, axis=1)[:, np.newaxis]
        o_2  = np.cross(normal, o_1)
        o_2 /=  np.linalg.norm(o_2, axis=1)[:, np.newaxis]

        R        = np.empty((self.param['intensity'], 3, 3))
        R[:,0,:] = o_1
        R[:,1,:] = o_2
        R[:,2,:] = normal

        direction = np.einsum('ij,ijk->ik', dir_local, R)
        return direction

    def generate_wavelength(self, direction):
        wtype = str.lower(self.param['wavelength_dist'])
        if wtype == 'monochrome':
            wavelength  = np.ones(self.param['intensity'], dtype = np.float64)
            wavelength *= self.param['wavelength']
        elif wtype == 'uniform':
            wavelength = np.random.uniform(
                self.param['wavelength_range'][0]
                ,self.param['wavelength_range'][1]
                ,self.param['intensity']
                )
        elif wtype == 'voigt':
            #random_wavelength = self.random_wavelength_normal
            #random_wavelength = self.random_wavelength_cauchy
            random_wavelength = self.random_wavelength_voigt
            wavelength = random_wavelength(self.param['intensity'])

            #doppler shift
            c = const.physical_constants['speed of light in vacuum'][0]
            wavelength *= 1 - (np.einsum('j,ij->i', self.param['velocity'], direction) / c)
        else:
            raise Exception(f'Wavelength distribution {wtype} unknown')

        return wavelength

    def random_wavelength_voigt(self, size=None):
        #Units: wavelength (angstroms), natural_linewith (1/s), temperature (eV)

        # Check for the trivial case.
        if (self.param['linewidth']  == 0.0 and self.param['temperature'] == 0.0):
            return np.ones(size)*self.param['wavelength']
        # Check for the Lorentzian case.
        if (self.param['temperature'] == 0.0):
            # I need to update the cauchy routine first.
            #raise NotImplementedError('Random Lorentzian distribution not implemented.')

            # TEMPORARY:
            # The raytracer cannot currently accept a zero temperature, so just add 1eV for now.
            self.param['temperature'] += 1.0

        # Check for the Gaussian case.
        if (self.param['linewidth']  == 0.0):
            return self.random_wavelength_normal(size)

        c = const.physical_constants['speed of light in vacuum'][0]
        amu_kg = const.physical_constants['atomic mass unit-kilogram relationship'][0]
        ev_J = const.physical_constants['electron volt-joule relationship'][0]

        # Natural line width.
        gamma = (self.param['linewidth'] * self.param['wavelength']**2 / (4 * np.pi * c * 1e10))

        # Doppler broadened line width.
        sigma = (np.sqrt(self.param['temperature'] / self.param['mass_number'] / amu_kg / c**2 * ev_J)
                  * self.param['wavelength'] )

        rand_wave  = xicsrt_voigt.voigt_random(gamma, sigma, size)
        rand_wave += self.param['wavelength']
        return rand_wave

    def random_wavelength_normal(self, size=None):
        #Units: wavelength (angstroms), temperature (eV)
        c       = const.physical_constants['speed of light in vacuum'][0]
        amu_kg  = const.physical_constants['atomic mass unit-kilogram relationship'][0]
        ev_J    = const.physical_constants['electron volt-joule relationship'][0]

        # Doppler broadened line width.
        sigma = ( np.sqrt(self.param['temperature'] / self.param['mass_number'] / amu_kg / c**2 * ev_J)
                  * self.param['wavelength'] )

        rand_wave = np.random.normal(self.param['wavelength'], sigma, size)
        return rand_wave

    def random_wavelength_cauchy(self, size=None):
        # This function needs to be updated to use the same definitions
        # as random_wavelength_voigt.
        #
        # As currently writen natual_linewidth is not used in a way
        # consistent with physical units.
        #
        # It also may make sense to add some sort of cutoff here.
        # the extreme tails of the distribution are not really useful
        # for ray tracing.
        fwhm = self.param['linewidth']
        rand_wave = cauchy.rvs(loc=self.param['wavelength'], scale=fwhm, size=size)
        return rand_wave

    def generate_weight(self):
        # Weight is not currently used within XICSRT but might be useful
        # in the future.
        w = np.ones((self.param['intensity']), dtype=np.float64)
        return w

    def generate_mask(self):
        m = np.ones((self.param['intensity']), dtype=np.bool)
        return m

    def ray_filter(self, rays):
        for filter in self.filter_objects:
            rays = filter.filter(rays)
        return rays

max_stars_repo_name: PrincetonUniversity/xicsrt

max_issues_repo_name: PrincetonUniversity/xicsrt

max_forks_repo_name: PrincetonUniversity/xicsrt

reason: import numpy,import scipy,from scipy


text:
"""
An example of drawing an adjustable ROI on a spectrogram view

* Click and drag on the image to select a region
* Adjust region size
* Click away from region to clear selection
* TODO: show selected range on screen
* TODO: show currently hovered frequency and time at all times
"""

import numpy as np
import pyqtgraph as pg
import PyQt5.QtWidgets as widgets
from PyQt5 import QtGui

from soundsep.gui.main import run_app
from soundsep.gui.stft_view import ScrollableSpectrogram, ScrollableSpectrogramConfig
from soundsep.gui.components.selection_box import SelectionBox
from soundsep.gui.components.spectrogram_view_box import SpectrogramViewBox

example_file = "/home/kevin/Data/from_pepe/2018_02_11 16-31-00/ch0.wav"
example_channel = 0

class SelectionBoxExample(widgets.QWidget):
    """Selection box example. Uses the spectrogram example as a background to draw on

    Draws when click and dragged
    """
    def __init__(self):
        super().__init__()
        self.title = "Selection Box Example"

        self.current_selection_roi = None
        self._init_spectrogram_view()

    def _init_spectrogram_view(self):
        spectrogram_viewbox = SpectrogramViewBox()
        self.spectrogram_panel = ScrollableSpectrogram(
            filename=example_file,
            channel=example_channel,
            config=ScrollableSpectrogramConfig(
                window_size=500,
                window_step=22,
                spectrogram_size=22050 * 5,
                cmap="turbo",
            ),
            viewBox=spectrogram_viewbox
        )

        spectrogram_viewbox.dragComplete.connect(self.on_drag_complete)
        spectrogram_viewbox.dragInProgress.connect(self.on_drag_in_progress)
        spectrogram_viewbox.clicked.connect(self.on_click)

        layout = widgets.QVBoxLayout()
        layout.addWidget(self.spectrogram_panel)
        self.setLayout(layout)
        self.spectrogram_panel.scroll_to(0)

    def on_drag_complete(self, from_, to):
        pass

    def on_drag_in_progress(self, from_, to):
        if self.current_selection_roi is None:
            self.draw_selection_box(from_, to)
        else:
            line = to - from_
            self.current_selection_roi.setPos([
                min(to.x(), from_.x()),
                min(to.y(), from_.y())
            ])
            self.current_selection_roi.setSize([
                np.abs(line.x()),
                np.abs(line.y())
            ])

    def on_click(self, at):
        self.clear_selection_box()

    def clear_selection_box(self):
        if self.current_selection_roi is not None:
            self.spectrogram_panel.removeItem(self.current_selection_roi)
            self.current_selection_roi.deleteLater()
            self.current_selection_roi = None

    def draw_selection_box(self, from_, to):
        self.current_selection_roi = SelectionBox(
            pos=from_,
            size=to - from_,
            pen=(156, 156, 100),
            rotatable=False,
            removable=False,
            maxBounds=self.spectrogram_panel.image.boundingRect()
        )
        self.current_selection_roi.sigRegionChangeFinished.connect(
            self.on_selection_changed
        )
        self.spectrogram_panel.addItem(self.current_selection_roi)

    def on_selection_changed(self, roi):
        x_bounds, y_bounds = roi.get_scaled_bounds(
            self.spectrogram_panel.image,
            self.spectrogram_panel.xlim(),
            self.spectrogram_panel.ylim(),
        )

if __name__ == "__main__":
    run_app(SelectionBoxExample)


text:
import os
from typing import Any, Callable, List, Optional, Tuple

import numpy as np
import pandas as pd
from dask_ml.compose import make_column_transformer
from dask_ml.preprocessing import Categorizer, OrdinalEncoder
from git_root import git_root
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer

from ..util import exposure_and_offset_to_weights

# taken from https://github.com/lorentzenchr/Tutorial_freMTPL2/blob/master/glm_freMTPL2_example.ipynb  # noqa: B950
# Modified to generate data sets of different sizes

def create_insurance_raw_data(verbose=False) -> None:
    """Do some basic processing on the data that we will later transform into our \
    benchmark data sets."""
    # load the datasets
    # first row (=column names) uses "", all other rows use ''
    # use '' as quotechar as it is easier to change column names
    df = pd.read_csv(
        "https://www.openml.org/data/get_csv/20649148/freMTPL2freq.arff", quotechar="'"
    )

    # rename column names '"name"' => 'name'
    df.rename(lambda x: x.replace('"', ""), axis="columns", inplace=True)
    df["IDpol"] = df["IDpol"].astype(np.int64)
    df.set_index("IDpol", inplace=True)

    df_sev = pd.read_csv(
        "https://www.openml.org/data/get_csv/20649149/freMTPL2sev.arff", index_col=0
    )

    # join ClaimAmount from df_sev to df:
    #   1. cut ClaimAmount at 100_000
    #   2. aggregate ClaimAmount per IDpol
    #   3. join by IDpol
    df_sev["ClaimAmountCut"] = df_sev["ClaimAmount"].clip(upper=100_000)
    df = df.join(df_sev.groupby(level=0).sum(), how="left")
    df.fillna(value={"ClaimAmount": 0, "ClaimAmountCut": 0}, inplace=True)

    # Check if there are IDpol in df_sev that do not match any IDPol in df.
    df2 = pd.merge(
        df_sev,
        df.loc[:, ["ClaimNb"]],
        left_index=True,
        right_index=True,
        how="outer",
        indicator=True,
    )
    if verbose:
        print(
            "There are {} rows in freMTPL2sev that do not have a matching IDpol in "
            "freMTPL2freq. They have a ClaimAmountCut of {}.".format(
                df2[df2._merge == "left_only"].shape[0],
                df2.ClaimAmountCut[df2._merge == "left_only"].sum(),
            )
        )

    round(df_sev.ClaimAmountCut.sum() - df.ClaimAmountCut.sum(), 2)

    if verbose:
        print(
            "Number or rows with ClaimAmountCut > 0 and ClaimNb == 0: {}".format(
                df[(df.ClaimAmountCut > 0) & (df.ClaimNb == 0)].shape[0]
            )
        )

    # 9116 zero claims
    if verbose:
        print(
            "Number or rows with ClaimAmountCut = 0 and ClaimNb >= 1: {}".format(
                df[(df.ClaimAmountCut == 0) & (df.ClaimNb >= 1)].shape[0]
            )
        )

    # Note: Zero claims must be ignored in severity models, because the support is
    # (0, inf) not [0, inf). Therefore, we define the number of claims with positive
    # claim amount for later use.
    df["ClaimNb_pos"] = df["ClaimNb"]
    df.loc[(df.ClaimAmount <= 0) & (df.ClaimNb >= 1), "ClaimNb_pos"] = 0

    # correct for unreasonable observations (that might be data error)
    # see case study paper
    df["ClaimNb"] = df["ClaimNb"].clip(upper=4)
    df["ClaimNb_pos"] = df["ClaimNb_pos"].clip(upper=4)
    df["Exposure"] = df["Exposure"].clip(upper=1)

    out_path = git_root("data/insurance.parquet")
    os.makedirs(os.path.dirname(out_path), exist_ok=True)
    df.to_parquet(out_path)

def get_categorizer(col_name: str, name="cat") -> Tuple[str, Categorizer]:
    """
    Get a dask_ml Categorizer.

    Categorizer only operates on object columns unless you explictily pass the column
    name.
    """
    return name, Categorizer(columns=[col_name])

def func_returns_df(
    fn: Callable[[pd.DataFrame], np.ndarray]
) -> Callable[[pd.DataFrame], pd.DataFrame]:
    """
    Take a function that takes a dataframe and returns a Numpy array, and return a \
    function that takes a dataframe and returns a dataframe.

    fn: Function that takes a dataframe and returns a numpy array
    Returns: Function that takes a dataframe and returns a dataframe with the values
             determined by the original function, and the index and columns of the original
             dataframe.
    """
    return lambda x: x.assign(**{x.columns[0]: fn(x)})

def gen_col_trans() -> Tuple[Any, List[str]]:
    """Generate a ColumnTransformer and list of names.

    The transformer corresponds to the GLM of the case study paper.

    Encodes k categories with k binary features (redundant).
    """
    column_trans = make_column_transformer(
        # VehPower 4, 5, 6, 7, 8, 9, drop=4
        (
            Pipeline(
                [
                    (
                        "cut_9",
                        FunctionTransformer(lambda x: np.minimum(x, 9), validate=False),
                    ),
                    get_categorizer("VehPower"),
                ]
            ),
            ["VehPower"],
        ),
        # VehAge intervals [0,1), [1, 10], (10, inf), drop=[1,10]
        (
            Pipeline(
                [
                    (
                        "bin",
                        FunctionTransformer(
                            func_returns_df(
                                lambda x: np.digitize(
                                    np.where(x == 10, 9, x), bins=[1, 10]
                                )
                            ),
                            validate=False,
                        ),
                    ),
                    get_categorizer("VehAge"),
                ]
            ),
            ["VehAge"],
        ),
        # DrivAge intervals [18,21), [21,26), [26,31), [31,41), [41,51), [51,71),[71,∞),
        # drop=[41,51)
        (
            Pipeline(
                [
                    (
                        "bin",
                        FunctionTransformer(
                            func_returns_df(
                                lambda x: np.digitize(x, bins=[21, 26, 31, 41, 51, 71])
                            ),
                            validate=False,
                        ),
                    ),
                    get_categorizer("DrivAge"),
                ]
            ),
            ["DrivAge"],
        ),
        (
            Pipeline(
                [
                    (
                        "cutat150",
                        FunctionTransformer(
                            lambda x: np.minimum(x, 150), validate=False
                        ),
                    )
                ]
            ),
            ["BonusMalus"],
        ),
        (Pipeline([get_categorizer("VehBrand")]), ["VehBrand"]),
        (Pipeline([get_categorizer("VehGas")]), ["VehGas"]),
        (FunctionTransformer(np.log, validate=False), ["Density"]),
        (Pipeline([get_categorizer("Region")]), ["Region"]),
        (
            Pipeline(
                [
                    get_categorizer("Area"),
                    ("OE", OrdinalEncoder()),
                    (
                        "plus_1",
                        FunctionTransformer(lambda x: x + 1, validate=False),
                    ),
                ]
            ),
            ["Area"],
        ),
        remainder="drop",
    )
    column_trans_names = [
        "VehPower_4",
        "VehPower_5",
        "VehPower_6",
        "VehPower_7",
        "VehPower_8",
        "VehPower_9",
        "VehAge_[0,1)",
        "VehAge_[1, 10]",
        "VehAge_(10,inf)",
        "DrivAge_[18,21)",
        "DrivAge_[21,26)",
        "DrivAge_[26,31)",
        "DrivAge_[31,41)",
        "DrivAge_[41,51)",
        "DrivAge_[51,71)",
        "DrivAge_[71,inf)",
        "BonusMalus",
        "VehBrand_B10",
        "VehBrand_B11",
        "VehBrand_B12",
        "VehBrand_B13",
        "VehBrand_B14",
        "VehBrand_B1",
        "VehBrand_B2",
        "VehBrand_B3",
        "VehBrand_B4",
        "VehBrand_B5",
        "VehBrand_B6",
        "VehGas_Diesel",
        "VehGas_Regular",
        "Density_log",
        "Region_R11",
        "Region_R21",
        "Region_R22",
        "Region_R23",
        "Region_R24",
        "Region_R25",
        "Region_R26",
        "Region_R31",
        "Region_R41",
        "Region_R42",
        "Region_R43",
        "Region_R52",
        "Region_R53",
        "Region_R54",
        "Region_R72",
        "Region_R73",
        "Region_R74",
        "Region_R82",
        "Region_R83",
        "Region_R91",
        "Region_R93",
        "Region_R94",
        "Area_ord",
    ]
    column_trans_names = [
        i
        for i in column_trans_names
        if i
        not in [
            "VehPower_4",
            "VehAge_[1, 10]",
            "DrivAge_[41,51)",
            "VehBrand_B1",
            "VehGas_Diesel",
            "Region_R24",
        ]
    ]
    return column_trans, column_trans_names

def add_noise(df: pd.DataFrame, noise: float) -> pd.DataFrame:
    """Add noise by swapping out data points."""
    np.random.seed(43212)
    for col in df.columns:
        if col in ["ClaimNb", "Exposure", "ClaimAmountCut", "ClaimNb_pos"]:
            continue
        swap = np.random.uniform(size=len(df)) < noise
        shuffle = np.random.choice(df[col], size=len(df))
        df.loc[swap, col] = shuffle[swap]

    return df

def compute_y_exposure(df, distribution):
    """Compute y and exposure depending on distribution."""
    if distribution == "poisson":
        exposure = df["Exposure"].values
        y = df["ClaimNb"].values / exposure
    elif distribution in ["gamma", "gaussian"]:
        exposure = df["ClaimNb_pos"].values
        y = df["ClaimAmountCut"].values / exposure
    elif "tweedie" in distribution:
        exposure = df["Exposure"].values
        y = df["ClaimAmountCut"].values / exposure
    elif distribution == "binomial":
        exposure = df["Exposure"].values
        df["HasClaim"] = (df["ClaimNb"] > 0).astype(np.int32)
        y = df["HasClaim"].values
    else:
        raise ValueError(
            "distribution must be one of ['poisson', 'gamma', 'tweedie', 'gaussian', "
            f"'binomial'] not {distribution}."
        )

    return y, exposure

def _read_insurance_data(
    num_rows: Optional[int], noise: Optional[float], distribution: str
) -> pd.DataFrame:
    path = git_root("data/insurance.parquet")
    if not os.path.exists(path):
        create_insurance_raw_data()
    df = pd.read_parquet(path)

    if distribution in ["gamma", "gaussian"]:
        df = df.query("ClaimAmountCut > 0")

    if num_rows is not None:
        # if we're oversampling, set default value for noise to 0.05
        # can be turned off by setting noise to zero
        if noise is None and num_rows > len(df):
            noise = 0.05
        df = df.sample(n=num_rows, replace=True, random_state=12345)

    if noise is not None:
        df = add_noise(df, noise=noise)
    return df

def generate_narrow_insurance_dataset(
    num_rows=None, noise=None, distribution="poisson"
) -> Tuple[pd.DataFrame, np.ndarray, np.ndarray]:
    """Generate the tutorial data set from the sklearn fork and save it to disk."""
    df = _read_insurance_data(num_rows, noise, distribution)

    col_trans_GLM1, _ = gen_col_trans()
    y, exposure = compute_y_exposure(df, distribution)

    return col_trans_GLM1.fit_transform(df), y, exposure

def generate_real_insurance_dataset(
    num_rows=None, noise=None, distribution="poisson"
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Load real insurance data set."""
    df = pd.read_parquet(git_root("data", "outcomes.parquet"))
    X = pd.read_parquet(git_root("data", "X.parquet"))

    if distribution != "poisson":
        raise NotImplementedError("distribution must be poisson")

    # restrict X and df to train set
    train_set = df["sample"] == "train"
    X = X.loc[train_set].reset_index(drop=True)
    df = df.loc[train_set].reset_index(drop=True)

    # subsample
    if num_rows is not None:
        idx = df.sample(n=num_rows).index
        df = df.loc[idx].reset_index(drop=True)
        X = X.loc[idx].reset_index(drop=True)

    # account for exposure and offsets
    y, weights = exposure_and_offset_to_weights(
        power=1,
        y=df["sanzkh02"],
        exposure=df["je"],
        offset=df["offset_kh_sach_frequenz"],
    )

    return X, y, weights

def generate_wide_insurance_dataset(
    num_rows=None, noise=None, distribution="poisson"
) -> Tuple[pd.DataFrame, np.ndarray, np.ndarray]:
    """Generate a version of the tutorial data set with many features."""
    df = _read_insurance_data(num_rows, noise, distribution)
    cat_cols = [
        "Area",
        "VehPower",
        "VehAge",
        "DrivAge",
        "BonusMalus",
        "VehBrand",
        "VehGas",
        "Region",
    ]

    transformer = make_column_transformer(
        (
            FunctionTransformer(),
            lambda x: [
                elmt
                for elmt in x.select_dtypes(["number"]).columns
                if elmt not in cat_cols
            ],
        ),
        (
            Pipeline([get_categorizer(col, "cat_" + col) for col in cat_cols]),
            cat_cols,
        ),
        remainder="drop",
    )
    y, exposure = compute_y_exposure(df, distribution)

    return transformer.fit_transform(df), y, exposure

def generate_intermediate_insurance_dataset(
    num_rows=None, noise=None, distribution="poisson"
) -> Tuple[pd.DataFrame, np.ndarray, np.ndarray]:
    """Generate the tutorial data set from the sklearn fork and save it to disk."""
    df = _read_insurance_data(num_rows, noise, distribution)
    df["BonusMalusClipped"] = df["BonusMalus"].clip(50, 100)

    col_trans_GLM1, _ = gen_col_trans()
    col_trans_GLM1.transformers.append(
        (
            "BonusMalusClipped",
            Pipeline([get_categorizer("BonusMalusClipped")]),
            ["BonusMalusClipped"],
        )
    )
    y, exposure = compute_y_exposure(df, distribution)

    return col_trans_GLM1.fit_transform(df), y, exposure

max_stars_repo_name: readthedocs-assistant/glum

max_issues_repo_name: readthedocs-assistant/glum

max_forks_repo_name: readthedocs-assistant/glum


text:
import os, inspect, glob, random, shutil
import numpy as np
from tensorflow.examples.tutorials.mnist import input_data

PACK_PATH = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+"/.."

class DataSet(object):

    def __init__(self):

        self.dataset = input_data.read_data_sets('MNIST_data', one_hot=True)

        self.num_tr = self.dataset.train.num_examples
        self.num_te = self.dataset.test.num_examples
        print("Training: %d" %(self.num_tr))
        print("Test    : %d" %(self.num_te))
        self.height, self.width = 28, 28

    def next_train(self, batch_size):

        x, y = self.dataset.train.next_batch(batch_size)

        return x, y

    def next_test(self, batch_size):

        x, y = self.dataset.test.next_batch(batch_size)

        return x, y

max_stars_repo_name: YeongHyeon/DCRAW

max_issues_repo_name: YeongHyeon/DCRAW

max_forks_repo_name: YeongHyeon/DCRAW


text:
import logging
from functools import partial
from math import inf
from typing import List, Tuple

import numpy
from cached_property import cached_property

from cu_pass.dpa_calculator.aggregate_interference_calculator.aggregate_interference_calculator import \
    AggregateInterferenceCalculator
from cu_pass.dpa_calculator.aggregate_interference_calculator.aggregate_interference_calculator_winnforum.support.maximum_move_list_distance_calculator import \
    MaximumMoveListDistanceCalculator
from cu_pass.dpa_calculator.cbsd.cbsd import CbsdCategories
from cu_pass.dpa_calculator.utilities import get_distance_between_two_points, get_dpa_calculator_logger, get_dpa_center, \
    Point
from reference_models.common import mpool
from reference_models.common.data import CbsdGrantInfo, ProtectedEntityType, ProtectionConstraint
from reference_models.dpa.dpa_mgr import Dpa
from reference_models.dpa.move_list import DpaType, findGrantsInsideNeighborhood, formInterferenceMatrix, \
    HIGH_FREQ_COCH, \
    LOW_FREQ_COCH

COCHANNEL_BANDWIDTH_IN_MEGAHERTZ = 10
HERTZ_IN_MEGAHERTZ = 1e6

class AggregateInterferenceCalculatorWinnforum(AggregateInterferenceCalculator):
    def calculate(self, distance: float, cbsd_category: CbsdCategories) -> float:
        maximum_move_distance_calculator = self._get_maximum_move_distance_calculator(distance, cbsd_category)
        return maximum_move_distance_calculator.get_max_distance()[cbsd_category]

    def get_expected_interference(self, distance: float, cbsd_category: CbsdCategories) -> float:
        maximum_move_distance_calculator = self._get_maximum_move_distance_calculator(neighborhood_distance=distance,
                                                                                cbsd_category=cbsd_category)
        self._logger.info(f'\t\tThreshold used: {maximum_move_distance_calculator.interference_threshold} dBm')
        self._log_beamwidth()
        return maximum_move_distance_calculator.get_expected_interference()

    def _log_beamwidth(self) -> None:
        self._logger.info(f'\t\tBeamwidth: {self._dpa.beamwidth} degrees')

    @property
    def _logger(self) -> logging.Logger:
        return get_dpa_calculator_logger()

    def _get_maximum_move_distance_calculator(self,
                                              neighborhood_distance: float,
                                              cbsd_category: CbsdCategories) -> MaximumMoveListDistanceCalculator:
        neighbor_grants_info = self._get_neighbor_grants_info(neighborhood_distance=neighborhood_distance,
                                                              cbsd_category=cbsd_category)
        return MaximumMoveListDistanceCalculator(
            dpa=self._dpa,
            grant_distances=self._grant_distances,
            grants_with_inband_frequencies=self._grants_with_inband_frequencies,
            interference_matrix_info=self._interference_matrix_info,
            neighbor_grants_info=neighbor_grants_info)

    def _get_neighbor_grants_info(self, neighborhood_distance: float, cbsd_category: CbsdCategories) -> Tuple[List[CbsdGrantInfo], List[int]]:
        neighbor_distances = (inf, neighborhood_distance, 0, 0) \
            if cbsd_category == CbsdCategories.B \
            else (neighborhood_distance, inf, 0, 0)
        return findGrantsInsideNeighborhood(
            grants=self._grants_with_inband_frequencies,
            constraint=self._protection_constraint,
            dpa_type=self._dpa_type,
            neighbor_distances=neighbor_distances)

    @cached_property
    def _interference_matrix_info(self) -> Tuple[numpy.ndarray, List[int], numpy.ndarray]:
        grant_indexes = list(range(len(self._grants_with_inband_frequencies)))
        return formInterferenceMatrix(
            grants=self._grants_with_inband_frequencies,
            grants_ids=grant_indexes,
            constraint=self._protection_constraint,
            inc_ant_height=self._dpa.radar_height,
            num_iter=Dpa.num_iteration,
            dpa_type=self._dpa_type)

    @cached_property
    def _grant_distances(self) -> List[float]:
        pool = mpool.Pool()
        distanceFunction = partial(
            get_distance_between_two_points,
            point2=self._dpa_center)

        return pool.map(distanceFunction, [Point(latitude=grant.latitude, longitude=grant.longitude) for grant in self._grants_with_inband_frequencies])

    @cached_property
    def _grants_with_inband_frequencies(self) -> List[CbsdGrantInfo]:
        """
        Set inband frequency to 1 MHz to avoid scaling, since EIRPs are already given per 10 MHz
        """
        return [grant._replace(low_frequency=self._low_inband_frequency,
                               high_frequency=self._high_inband_frequency)
                for index, grant in enumerate(self._grants)]

    @property
    def _high_inband_frequency(self) -> float:
        return self._low_inband_frequency + self._frequency_in_hertz(COCHANNEL_BANDWIDTH_IN_MEGAHERTZ)

    @property
    def _low_inband_frequency(self) -> float:
        return self._frequency_in_hertz(self._first_inband_channel[0])

    @property
    def _first_inband_channel(self) -> Tuple[float, float]:
        return next(channel for channel in self._dpa._channels if self._channel_is_cochannel(channel))

    def _channel_is_cochannel(self, channel: Tuple[float, float]) -> bool:
        return self._frequency_in_hertz(channel[0]) >= LOW_FREQ_COCH and self._frequency_in_hertz(channel[1]) <= HIGH_FREQ_COCH

    @staticmethod
    def _frequency_in_hertz(frequency_in_megahertz: float) -> float:
        return frequency_in_megahertz * HERTZ_IN_MEGAHERTZ

    @property
    def _grants(self) -> List[CbsdGrantInfo]:
        return [cbsd.to_grant() for cbsd in self._cbsds]

    @property
    def _protection_constraint(self) -> ProtectionConstraint:
        return ProtectionConstraint(latitude=self._dpa_center.latitude,
                                    longitude=self._dpa_center.longitude,
                                    low_frequency=self._low_inband_frequency,
                                    high_frequency=self._high_inband_frequency,
                                    entity_type=ProtectedEntityType.DPA)

    @property
    def _dpa_center(self) -> Point:
        return get_dpa_center(dpa=self._dpa)

    @property
    def _dpa_type(self) -> DpaType:
        return DpaType.CO_CHANNEL

max_stars_repo_name: NSF-Swift/Spectrum-Access-System

max_issues_repo_name: NSF-Swift/Spectrum-Access-System

max_forks_repo_name: NSF-Swift/Spectrum-Access-System


text:
import numpy as np
import pickle
import torch
from torch import nn
from torch.utils.data import DataLoader
from torch.utils.data.dataloader import default_collate
from torch.nn import functional as F
from torchvision import transforms
import copy
import data
from encoder import Classifier
from vae_models import AutoEncoder

###################
## Loss function ##
###################

def loss_fn_kd(scores, target_scores, T=2.):
    """Compute knowledge-distillation (KD) loss given [scores] and [target_scores].

    Both [scores] and [target_scores] should be tensors, although [target_scores] should be repackaged.
    'Hyperparameter': temperature"""

    device = scores.device

    log_scores_norm = F.log_softmax(scores / T, dim=1)
    targets_norm = F.softmax(target_scores / T, dim=1)

    # if [scores] and [target_scores] do not have equal size, append 0's to [targets_norm]
    n = scores.size(1)
    if n>target_scores.size(1):
        n_batch = scores.size(0)
        zeros_to_add = torch.zeros(n_batch, n-target_scores.size(1))
        zeros_to_add = zeros_to_add.to(device)
        targets_norm = torch.cat([targets_norm.detach(), zeros_to_add], dim=1)

    # Calculate distillation loss (see e.g., Li and Hoiem, 2017)
    KD_loss_unnorm = -(targets_norm * log_scores_norm)
    KD_loss_unnorm = KD_loss_unnorm.sum(dim=1)                      #--> sum over classes
    KD_loss_unnorm = KD_loss_unnorm.mean()                          #--> average over batch

    # normalize
    KD_loss = KD_loss_unnorm * T**2

    return KD_loss

def loss_fn_kd_binary(scores, target_scores, T=2.):
    """Compute binary knowledge-distillation (KD) loss given [scores] and [target_scores].

    Both [scores] and [target_scores] should be tensors, although [target_scores] should be repackaged.
    'Hyperparameter': temperature"""

    device = scores.device

    scores_norm = torch.sigmoid(scores / T)
    targets_norm = torch.sigmoid(target_scores / T)

    # if [scores] and [target_scores] do not have equal size, append 0's to [targets_norm]
    n = scores.size(1)
    if n>target_scores.size(1):
        n_batch = scores.size(0)
        zeros_to_add = torch.zeros(n_batch, n-target_scores.size(1))
        zeros_to_add = zeros_to_add.to(device)
        targets_norm = torch.cat([targets_norm, zeros_to_add], dim=1)

    # Calculate distillation loss
    KD_loss_unnorm = -( targets_norm * torch.log(scores_norm) + (1-targets_norm) * torch.log(1-scores_norm) )
    KD_loss_unnorm = KD_loss_unnorm.sum(dim=1)      #--> sum over classes
    KD_loss_unnorm = KD_loss_unnorm.mean()          #--> average over batch

    # normalize
    KD_loss = KD_loss_unnorm * T**2

    return KD_loss

##-------------------------------------------------------------------------------------------------------------------##

#############################
## Data-handling functions ##
#############################

def get_data_loader(dataset, batch_size, cuda=False, collate_fn=None, drop_last=False, augment=False):
    '''Return <DataLoader>-object for the provided <DataSet>-object [dataset].'''

    # If requested, make copy of original dataset to add augmenting transform (without altering original dataset)
    if augment:
        dataset_ = copy.deepcopy(dataset)
        dataset_.transform = transforms.Compose([dataset.transform, *data.AVAILABLE_TRANSFORMS['augment']])
    else:
        dataset_ = dataset

    # Create and return the <DataLoader>-object
    return DataLoader(
        dataset_, batch_size=batch_size, shuffle=True,
        collate_fn=(collate_fn or default_collate), drop_last=drop_last,
        **({'num_workers': 0, 'pin_memory': True} if cuda else {})
    )

def label_squeezing_collate_fn(batch):
    x, y = default_collate(batch)
    return x, y.long().squeeze()

def to_one_hot(y, classes):
    '''Convert a nd-array with integers [y] to a 2D "one-hot" tensor.'''
    c = np.zeros(shape=[len(y), classes], dtype='float32')
    c[range(len(y)), y] = 1.
    c = torch.from_numpy(c)
    return c

##-------------------------------------------------------------------------------------------------------------------##

##########################################
## Object-saving and -loading functions ##
##########################################

def save_object(object, path):
    with open(path + '.pkl', 'wb') as f:
        pickle.dump(object, f, pickle.HIGHEST_PROTOCOL)

def load_object(path):
    with open(path + '.pkl', 'rb') as f:
        return pickle.load(f)

##-------------------------------------------------------------------------------------------------------------------##

################################
## Model-inspection functions ##
################################

def count_parameters(model, verbose=True):
    '''Count number of parameters, print to screen.'''
    total_params = learnable_params = fixed_params = 0
    for param in model.parameters():
        n_params = index_dims = 0
        for dim in param.size():
            n_params = dim if index_dims==0 else n_params*dim
            index_dims += 1
        total_params += n_params
        if param.requires_grad:
            learnable_params += n_params
        else:
            fixed_params += n_params
    if verbose:
        print("--> this network has {} parameters (~{} million)"
              .format(total_params, round(total_params / 1000000, 1)))
        print("      of which: - learnable: {} (~{} million)".format(learnable_params,
                                                                     round(learnable_params / 1000000, 1)))
        print("                - fixed: {} (~{} million)".format(fixed_params, round(fixed_params / 1000000, 1)))
    return total_params, learnable_params, fixed_params

def print_model_info(model, title="MODEL"):
    '''Print information on [model] onto the screen.'''
    print("Model-name: \"" + model.name + "\"")
    print(40*"-" + title + 40*"-")
    print(model)
    print(90*"-")
    _ = count_parameters(model)
    print(90*"-" + "\n\n")

##-------------------------------------------------------------------------------------------------------------------##

#################################
## Custom-written "nn-Modules" ##
#################################

class Identity(nn.Module):
    '''A nn-module to simply pass on the input data.'''
    def forward(self, x):
        return x

    def __repr__(self):
        tmpstr = self.__class__.__name__ + '()'
        return tmpstr

class Reshape(nn.Module):
    '''A nn-module to reshape a tensor to a 4-dim "image"-tensor with [image_channels] channels.'''
    def __init__(self, image_channels):
        super().__init__()
        self.image_channels = image_channels

    def forward(self, x):
        batch_size = x.size(0)   # first dimenstion should be batch-dimension.
        image_size = int(np.sqrt(x.nelement() / (batch_size*self.image_channels)))
        return x.view(batch_size, self.image_channels, image_size, image_size)

    def __repr__(self):
        tmpstr = self.__class__.__name__ + '(channels = {})'.format(self.image_channels)
        return tmpstr

class ToImage(nn.Module):
    '''Reshape input units to image with pixel-values between 0 and 1.

    Input:  [batch_size] x [in_units] tensor
    Output: [batch_size] x [image_channels] x [image_size] x [image_size] tensor'''

    def __init__(self, image_channels=1):
        super().__init__()
        # reshape to 4D-tensor
        self.reshape = Reshape(image_channels=image_channels)
        # put through sigmoid-nonlinearity
        self.sigmoid = nn.Sigmoid()

    def forward(self, x):
        x = self.reshape(x)
        x = self.sigmoid(x)
        return x

    def image_size(self, in_units):
        '''Given the number of units fed in, return the size of the target image.'''
        image_size = np.sqrt(in_units/self.image_channels)
        return image_size

class Flatten(nn.Module):
    '''A nn-module to flatten a multi-dimensional tensor to 2-dim tensor.'''
    def forward(self, x):
        batch_size = x.size(0)   # first dimenstion should be batch-dimension.
        return x.view(batch_size, -1)

    def __repr__(self):
        tmpstr = self.__class__.__name__ + '()'
        return tmpstr

##-------------------------------------------------------------------------------------------------------------------##

max_stars_repo_name: sanketvmehta/continual-learning

max_issues_repo_name: XiXiRuPan/continual-learning

max_forks_repo_name: XiXiRuPan/continual-learning


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

from collections import OrderedDict
from numpy import linspace, isnan
from numpy.testing import assert_equal
from pandas import DataFrame, Series
from json_tricks import dumps, loads
from tests.test_bare import nonpdata

COLUMNS = OrderedDict((
	('name', ('Alfa', 'Bravo', 'Charlie', 'Delta', 'Echo', 'Foxtrot', 'Golf',
		'Hotel', 'India', 'Juliett',)),
	('count', linspace(0, 10, 10, dtype=int)),
	('real', linspace(0, 7.5, 10, dtype=float)),
	('special', (float('NaN'), float('+inf'), float('-inf'), float('+0'),
		float('-0'), 1, 2, 3, 4, 5)),
	#todo: other types?
))

def test_pandas_dataframe():
	df = DataFrame(COLUMNS, columns=tuple(COLUMNS.keys()))
	txt = dumps(df, allow_nan=True)
	back = loads(txt)
	assert isnan(back.ix[0, -1])
	assert (df.equals(back))
	assert (df.dtypes == back.dtypes).all()
	df = DataFrame(COLUMNS, columns=tuple(COLUMNS.keys()))
	txt = dumps(df, primitives=True, allow_nan=True)
	back = loads(txt)
	assert isinstance(back, dict)
	assert isnan(back['special'][0])
	assert all(df.index.values == tuple(back.pop('index')))
	for name, col in back.items():
		assert name in COLUMNS
		assert_equal(list(COLUMNS[name]), col)

def test_pandas_series():
	for name, col in COLUMNS.items():
		ds = Series(data=col, name=name)
		txt = dumps(ds, allow_nan=True)
		back = loads(txt)
		assert (ds.equals(back))
		assert ds.dtype == back.dtype
	for name, col in COLUMNS.items():
		ds = Series(data=col, name=name)
		txt = dumps(ds, primitives=True, allow_nan=True)
		back = loads(txt)
		assert isinstance(back, dict)
		assert_equal(ds.index.values, back['index'])
		assert_equal(ds.values, back['data'])

def test_pandas_mixed_with_other_types():
	df = DataFrame(COLUMNS, columns=tuple(COLUMNS.keys()))
	mixed = dict(
		complex=1+42j,
		frames=[df, df],
		**nonpdata
	)
	txt = dumps(mixed, allow_nan=True)
	back = loads(txt)
	assert mixed['frames'][0].equals(back['frames'][0]) and mixed['frames'][1].equals(back['frames'][1])
	del mixed['frames'], back['frames']  # cannot compare dataframes with '=='
	assert mixed == back

reason: from numpy