text:
"""
modified from https://github.com/lllyasviel/sketchKeras
"""
import cv2
import numpy as np
import os
from keras.models import load_model
from scipy import ndimage
from .xdog_blend import *

def get_light_map_single(img):
    gray = img
    gray = gray[None]
    gray = gray.transpose((1,2,0))
    blur = cv2.GaussianBlur(gray, (0, 0), 3)
    gray = gray.reshape((gray.shape[0],gray.shape[1]))
    highPass = gray.astype(int) - blur.astype(int)
    highPass = highPass.astype(np.float)
    highPass = highPass / 128.0
    return highPass

def normalize_pic(img):
    img = img / np.max(img)
    return img

def resize_img_512_3d(img):
    zeros = np.zeros((1,3,512,512), dtype=np.float)
    zeros[0 , 0 : img.shape[0] , 0 : img.shape[1] , 0 : img.shape[2]] = img
    return zeros.transpose((1,2,3,0))

def to_keras_enhanced(img):
    mat = img.astype(np.float)
    mat[mat<0.1] = 0
    mat = - mat + 1
    mat = mat * 255.0
    mat[mat < 0] = 0
    mat[mat > 255] = 255
    mat=mat.astype(np.uint8)
    mat = ndimage.median_filter(mat, 1)

    return mat

def get_light_map(img):
    from_mat = img
    width = float(from_mat.shape[1])
    height = float(from_mat.shape[0])
    new_width = 0
    new_height = 0

    if (width > height):
        if width != 512:
            from_mat = cv2.resize(from_mat, (512, int(512 / width * height)), interpolation=cv2.INTER_AREA)
        new_width = 512
        new_height = int(512 / width * height)
    else:
        if height != 512:
            from_mat = cv2.resize(from_mat, (int(512 / height * width), 512), interpolation=cv2.INTER_AREA)
        new_width = int(512 / height * width)
        new_height = 512

    from_mat = from_mat.transpose((2, 0, 1))
    light_map = np.zeros(from_mat.shape, dtype=np.float)
    for channel in range(3):
        light_map[channel] = get_light_map_single(from_mat[channel])

    light_map = normalize_pic(light_map)
    light_map = resize_img_512_3d(light_map)

    return new_height, new_width, light_map

mod = None
def get_mod():
    global mod
    if mod is None:
        mod = load_model(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'sketchKeras.h5'))
    return mod

def get_keras_enhanced(img):
    mod = get_mod()
    new_height, new_width, light_map = get_light_map(img)

    # TODO: batch this!
    line_mat = mod.predict(light_map, batch_size=1)

    line_mat = line_mat.transpose((3, 1, 2, 0))[0]
    line_mat = line_mat[0:int(new_height), 0:int(new_width), :]
    line_mat = np.amax(line_mat, 2)

    keras_enhanced = to_keras_enhanced(line_mat)
    return keras_enhanced

def get_keras_high_intensity(img, intensity=1.7):
    keras_img = get_keras_enhanced(img)
    return add_intensity(keras_img, intensity)

def batch_keras_enhanced(img_list):
    mod = get_mod()
    light_maps = list(map(get_light_map, img_list))

    hw_list = [(h, w) for h, w, _ in light_maps]
    light_map = [l for _, _, l in light_maps]

    batch_light_map = np.concatenate(light_map, axis=0)
    batch_line_mat = mod.predict(batch_light_map, batch_size=len(img_list))

    mat_list = np.array_split(batch_line_mat, len(img_list))
    batch_line_mat = map(lambda map: map.transpose((3, 1, 2, 0))[0], mat_list)

    result_list = []
    for line_mat, (new_height, new_width) in zip(batch_line_mat, hw_list):
        mat = line_mat[0:int(new_height), 0:int(new_width), :]
        mat = np.amax(mat, 2)
        keras_enhanced = to_keras_enhanced(mat)
        result_list.append(keras_enhanced)

    return result_list


text:
# coding: utf-8
import numpy as np
from numpy import pi, sqrt, sin, cos, deg2rad, rad2deg, arctan2

##### Constants ###############################################
G0 = 9.80665

### WGS84 ###
A = 6378137.0           # Equatorial radius [m]
ONE_F = 298.257223563   # Inverse flattening [-]
GM = 3.986004418e14     # Gravitational constant [m3/s2]
OMEGA = 7.292115e-5     # Angular velocity of the Earth [rad/s]

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

# lat: geodetic latitude
# latc: geocentric latitude

def posECEFfromLLH(posLLH_):
    """
    Args:
        posLLH_ (np.array 3x1) : position in LLH coordinate [deg, deg, m]
    Return:
        (np.array 3x1) : position in ECEF coordinate [m, m, m]
    """
    lat = deg2rad(posLLH_[0])
    lon = deg2rad(posLLH_[1])
    alt = posLLH_[2]
    a = A                           # Equatorial radius [m]
    f = 1.0 / ONE_F                 # Flattening [-]
    e2 = 2.0 * f - f * f            # Square of eccentricity e
    W = sqrt(1.0 - e2 * sin(lat) * sin(lat))
    N = a / W                       # Radius of prime vertical circle [m]

    posECEF_ = np.empty(3)
    posECEF_[0] = (N + alt) * cos(lat) * cos(lon)
    posECEF_[1] = (N + alt) * cos(lat) * sin(lon)
    posECEF_[2] = (N * (1 - e2) + alt) * sin(lat)

    return posECEF_

def posLLHfromECEF(posECEF_):
    """
    This function is an approximation.
    Please refer https://www.enri.go.jp/~fks442/K_MUSEN/1st/1st060428rev2.pdf.

    Args:
        posECEF_ (np.array 3x1) : [deg, deg, m]
    Return:
        (np.array 3x1) : position in LLH coordinate [deg, deg, m]
    """
    a = A                           # Equatorial radius [m]
    f = 1.0 / ONE_F                 # Flattening [-]
    e2 = 2.0 * f - f * f            # Square of eccentricity e
    b = a * (1.0 - f)               # Polar radius [m]
    ed2 = e2 / (1.0 - e2)           # Square of second eccentricity e'
    p = sqrt(posECEF_[0]**2 + posECEF_[1]**2)    # Distance from the rotation axis of the Earth [m]
    theta = arctan2(posECEF_[2] * a, p * b)      # [rad]

    posLLH = np.empty(3)
    posLLH[0] = rad2deg(arctan2((posECEF_[2] + ed2 * b * sin(theta)**3), p - e2 * a * cos(theta)**3))
    posLLH[1] = rad2deg(arctan2(posECEF_[1], posECEF_[0]))
    posLLH[2] = p / cos(deg2rad(posLLH[0])) - a / sqrt(1.0 - e2 * sin(deg2rad(posLLH[0]))**2)

    return posLLH  # deg, deg, m

def dcmECEF2NEDfromLLH(posLLH_):  # deg, deg, m
    """
    Args:
        posLLH_ (np.array 3x1) : [deg, deg, m]
    Return:
        dcm (np.array 3x3) : DCM from ECEF to NED
    """
    lat = deg2rad(posLLH_[0])
    lon = deg2rad(posLLH_[1])

    dcm = np.array([[-sin(lat) * cos(lon), -sin(lat)*sin(lon),  cos(lat)],
                    [-sin(lon),             cos(lon),           0.      ],
                    [-cos(lat) * cos(lon), -cos(lat)*sin(lon), -sin(lat)]])
    return dcm

def dcmNED2ECEFfromLLH(posLLH_):
    """
    Args:
        posLLH_ (np.array 3x1) : [deg, deg, m]
    Return:
        dcm (np.array 3x3) : DCM from NED to ECEF
    """
    return dcmECEF2NEDfromLLH(posLLH_).T

def dcmECI2ECEF(second):
    """
    Args:
        second (double) : time from reference time[s]
    Return:
        dcm (np.array 3x3) : DCM from ECI to ECEF
    """
    omega = OMEGA
    theta = omega * second

    dcm = np.array([[ cos(theta), sin(theta), 0.],
                    [-sin(theta), cos(theta), 0.],
                    [0.,          0.,         1.]])
    return dcm

def dcmECEF2ECI(second):
    """
    Args:
        second (double) : time from reference time[s]
    Return:
        dcm (np.array 3x3) : DCM from ECEF to ECI
    """
    return dcmECI2ECEF(second).T

def dcmNED2BODY(azimuth_deg, elevation_deg, roll_deg):
    """
    Args:
        azimuth_deg (double) : azimuth of attitude [deg]
        elevation_deg (double) : elevation of attitude [deg]
        roll_deg (double) : roll of attitude [deg]
    Return:
        dcm (np.array 3x3) : DCM from NED to Body
    """
    az = deg2rad(azimuth_deg)
    el = deg2rad(elevation_deg)
    ro = deg2rad(roll_deg)

    dcm = np.array([
        [cos(el)*cos(az), cos(el)*sin(az), -sin(el)],
        [-cos(ro)*sin(az)+sin(ro)*sin(el)*cos(az), cos(ro)*cos(az)+sin(ro)*sin(el)*sin(az), sin(ro)*cos(el)],
        [sin(ro)*sin(az)+cos(ro)*sin(el)*cos(az),  -sin(ro)*cos(az)+cos(ro)*sin(el)*sin(az),  cos(ro)*cos(el)]
        ])
    return dcm

def dcmBODY2NED(azimuth_deg, elevation_deg, roll_deg):
    """
    Args:
        azimuth_deg (double) : azimuth of attitude [deg]
        elevation_deg (double) : elevation of attitude [deg]
        roll_deg (double) : roll of attitude [deg]
    Return:
        dcm (np.array 3x3) : DCM from Body to NED
    """
    dcm_inv = dcmNED2BODY(azimuth_deg, elevation_deg, roll_deg)
    return dcm_inv.T

def posECIfromECEF(posECEF_, second):
    """
    Args:
        posECEF_ (np.array 3x1) : position in ECEF coordinate [m, m, m]
        second (double) : time from reference time [s]
    Return:
        (np.array 3x1) : position in ECI coordinate [m, m, m]
    """
    return np.matmul(dcmECEF2ECI(second), posECEF_)

def posECEFfromECI(posECI_, second):
    return np.matmul(dcmECI2ECEF(second), posECI_)

def velECEFfromECI(posECI_, velECI_, second):
    omega = OMEGA
    omegaECI2ECEF_ = np.array([[0.,   -omega, 0.],
                               [omega, 0.,    0.],
                               [0.,    0.,    0.]])
    dcmECI2ECEF_ = dcmECI2ECEF(second)

    return np.matmul(dcmECI2ECEF_, velECI_ - np.matmul(omegaECI2ECEF_, posECI_))

def velECIfromECEF(posECEF_, velECEF_, second):
    omega = OMEGA
    omegaECI2ECEF_ = np.array([[0.,   -omega, 0.],
                               [omega, 0.,    0.],
                               [0.,    0.,    0.]])
    dcmECEF2ECI_ = dcmECEF2ECI(second)

    return np.matmul(dcmECEF2ECI_, velECEF_ + np.matmul(omegaECI2ECEF_, posECEF_))

def velECEFfromNED(posLLH_, velNED_):
    dcmNED2ECEF_ = dcmNED2ECEFfromLLH(posLLH_)
    return np.matmul(dcmNED2ECEF_, velNED_)

def velNEDfromECEF(posECEF_, velECEF_):
    posLLH_ = posLLHfromECEF(posECEF_)
    dcmECEF2NED_ = dcmECEF2NEDfromLLH(posLLH_)
    return np.matmul(dcmECEF2NED_, velECEF_)

def posLLHfromECI(posECI_, second):
    posECEF_ = posECEFfromECI(posECI_, second)
    posLLH_ = posLLHfromECEF(posECEF_)
    return posLLH_

def posECIfromLLH(posLLH_, second):
    posECEF_ = posECEFfromLLH(posLLH_)
    posECI_ = posECIfromECEF(posECEF_, second)
    return posECI_

# def velNEDfromECI(posECI_, velECI_, dcmECI2NED_):
#     omega = OMEGA
#     omegaECI2ECEF_ = np.array([[0.,   -omega, 0.],
#                                [omega, 0.,    0.],
#                                [0.,    0.,    0.]])
#
#     return np.matmul(dcmECI2NED_, (velECI_ - np.matmul(omegaECI2ECEF_, posECI_)))

def velNEDfromECI(posECI_, velECI_, second):
    posECEF_ = posECEFfromECI(posECI_, second)
    velECEF_ = velECEFfromECI(posECI_, velECI_, second)
    velNED_ = velNEDfromECEF(posECEF_, velECEF_)
    return velNED_

def velECIfromNED(posLLH_, velNED_, second):
    posECEF_ = posECEFfromLLH(posLLH_)
    velECEF_ = velECEFfromNED(posLLH_, velNED_)
    velECI_ = velECIfromECEF(posECEF_, velECEF_, second)
    return velECI_

def posLLH_IIP(posLLH_, velNED_):
    g0 = G0
    vel_north_ = velNED_[0]
    vel_east_ = velNED_[1]
    vel_up_ = - velNED_[2]
    h = posLLH_[2]
    if h > 0.0:
        tau = 1.0 / g0 * (vel_up_ + sqrt(vel_up_**2 + 2 * h * g0))
    else:
        tau = 0.0
    dist_IIP_from_now_NED = np.array([vel_north_ * tau, vel_east_ * tau, h])
    posECEF_ = posECEFfromLLH(posLLH_)
    dcmECEF2NED_ = dcmECEF2NEDfromLLH(posLLH_)
    dcmNED2ECEF_ = dcmECEF2NED_.T
    posECEF_IIP_ = posECEF_ + np.matmul(dcmNED2ECEF_, dist_IIP_from_now_NED)
    return posLLHfromECEF(posECEF_IIP_)

def posLLH_IIP_fromECEF(posECEF_, velECEF_):
    posLLH_ = posLLHfromECEF(posECEF_)
    velNED_ = velNEDfromECEF(posECEF_, velECEF_)
    return posLLH_IIP(posLLH_, velNED_)

### utility functions ###

def distance_surface(posLLH1, posLLH2):
    earth_radius = A            # Earth radius ~ Equatorial radius
    if len(posLLH1) == 2:
        posLLH1 = np.append(posLLH1, 0.)
    if len(posLLH2) == 2:
        posLLH2 = np.append(posLLH2, 0.)
    posECEF1 = posECEFfromLLH(posLLH1)
    posECEF2 = posECEFfromLLH(posLLH2)

    costheta = np.dot(posECEF1, posECEF2) / (np.linalg.norm(posECEF1) * np.linalg.norm(posECEF2))
    costheta = np.clip(costheta, -1., 1.)
    theta = np.arccos(costheta)     # Angle at the center of the Earth
    return earth_radius * theta

def posfromDistanceAndAzimuth(latlon, disaz):  # ([deg], [deg]), ([m], [deg])
    a = A                           # Equatorial radius [m]
    f = 1.0 / ONE_F                 # Flattening [-]
    e2 = 2.0 * f - f * f            # Square of eccentricity e

    lat, lon = latlon
    d, az = disaz

    theta_lat = deg2rad(lat)
    N = a / sqrt(1.0 - e2 * sin(theta_lat)**2)  # Prime vertical radius of curvature
    Re = N * sqrt(cos(theta_lat)**2 + ((1.0 - e2) * sin(theta_lat))**2)

    theta = d / Re
    phi = pi - deg2rad(az)

    x = sin(theta) * cos(phi)
    y = sin(theta) * sin(phi)
    z = cos(theta)

    alpha = pi * 0.5 - deg2rad(lat)

    X = x * cos(alpha) + z * sin(alpha)
    Y = y
    Z = -x * sin(alpha) + z * cos(alpha)

    theta_new = np.arccos(Z)
    phi_new = np.sign(Y) * np.arccos(X / sqrt(X**2 + Y**2))

    lat_ans = rad2deg(pi * 0.5 - theta_new)
    lon_ans = lon + rad2deg(phi_new)

    lat_ans = (lat_ans + 90) % 180 - 90
    lon_ans = (lon_ans + 180) % 360 - 180

    return (lat_ans, lon_ans)

def lonlat2xyratio(lat):  # [lon[deg], lat[deg]] to [x[m], y[m]] ratio
    a = A                           # Equatorial radius [m]
    f = 1.0 / ONE_F                 # Flattening [-]
    e2 = 2.0 * f - f * f            # Square of eccentricity e

    theta_lat = deg2rad(lat)
    N = a / sqrt(1.0 - e2 * sin(theta_lat)**2)  # Prime vertical radius of curvature

    dxdlon = N * cos(theta_lat)
    dydlat = N * (1.0 - e2) / (1.0 - e2 * sin(theta_lat)**2)

    return np.array([dxdlon, dydlat]).T * deg2rad(1.)

# points in hull must not be overlapping (hull[0] != hull[-1]), and ordered in anti-clockwise
def extend_lonlat_hull(lonlat_hull, distance, dlondlat=[0, 0], max_deg=30):  # [deg, deg], [m], [deg, deg]
    lonlat_ave = np.mean(lonlat_hull, axis=0)
    lola2xy = lonlat2xyratio(lonlat_ave[1])
    dxdy = dlondlat * lola2xy
    xy_hull = lonlat_hull * lola2xy
    xy_hull_new = extend_xy_hull(xy_hull, distance, dxdy, max_deg)
    lonlat_hull_new = xy_hull_new / lola2xy
    return lonlat_hull_new

def extend_xy_hull(xy_hull, distance, dxdy=[0, 0], max_deg=30):
    xy0 = xy_hull
    xy1 = np.roll(xy_hull, -1, axis=0)
    dxy = xy1 - xy0
    udxy = dxy / np.array([np.linalg.norm(dxy, axis=1)]).T
    ndxy = np.matmul(udxy, [[0, -1], [1, 0]])
    ddxy = distance * ndxy
    xy0n = xy0 + ddxy + dxdy
    xy1n = xy1 + ddxy + dxdy

    costheta = np.sum(ndxy * np.roll(ndxy, -1, axis=0), axis=1).clip(-1, 1)
    theta = np.arccos(costheta)
    xy_corners = []
    for theta_, ddxy_, xy1_ in zip(theta, ddxy, xy1):
        N_ = (theta_ // deg2rad(max_deg))
        n_ = np.arange(0, N_, 1)
        theta_n = theta_ * (n_ + 1) / (N_ + 1)
        An_T = np.array([[cos(theta_n), -sin(theta_n)], [sin(theta_n), cos(theta_n)]]).transpose(2, 1, 0)
        ddxy_n = np.matmul(ddxy_, An_T)
        xy1_nn = xy1_ + ddxy_n + dxdy
        xy_corners += [xy1_nn]

    xy_ans_list = []
    for xy0n_, xy1n_, xycs_ in zip(xy0n, xy1n, xy_corners):
        xy_ans_list += [xy0n_, xy1n_]
        xy_ans_list += list(xycs_)
    return np.array(xy_ans_list)

ext: py

max_stars_repo_name: KazufumiIkemoto/OpenTsiolkovsky

max_issues_repo_name: KazufumiIkemoto/OpenTsiolkovsky

max_forks_repo_name: KazufumiIkemoto/OpenTsiolkovsky

reason: import numpy,from numpy


text:
# Wilcoxon signed-rank test
from numpy.random import seed
from numpy.random import randn
from scipy.stats import wilcoxon
# seed the random number generator
seed(1)
# generate two independent samples
data1 = 5 * randn(100) + 50
data2 = 5 * randn(100) + 51
# compare samples
stat, p = wilcoxon(data1, data2)
print('Statistics=%.3f, p=%.3f' % (stat, p))
# interpret
alpha = 0.05
if p > alpha:
	print('Same distribution (fail to reject H0)')
else:
	print('Different distribution (reject H0)')

reason: from numpy,from scipy


text:
import sys

import numpy
import PIL.Image
import PIL.ImageDraw
import PIL.ImageFont

__all__ = [
    'char_to_image',
    'char_to_image_dict',
    'chars_to_images',
    'filename_to_font'
]

FOREGROUND_COLOR = 255
BACKGROUND_COLOR = 0

def char_to_image(char, font, *, size):
    empty_image = PIL.Image.new("L", size, color=BACKGROUND_COLOR)

    try:
        image = empty_image.copy()
        PIL.ImageDraw.Draw(image).text((0, 0),
                                       char,
                                       font=font,
                                       fill=FOREGROUND_COLOR)
        return numpy.array(image, dtype=numpy.uint8)
    except UnicodeEncodeError:
        print(("Could not render the unicode character \\u{:04X}"
               .format(ord(char))),
              file=sys.stderr)
        return None

def char_to_image_dict(chars, font):
    size = tuple(map(max, zip(*[font.getsize(char) for char in chars])))
    pairs = {char: char_to_image(char, font, size=size) for char in chars}
    return {char: image.tolist() for char, image in pairs.items()
            if image is not None and (image != BACKGROUND_COLOR).any()}

def chars_to_images(chars, font, unknown_char='\uFFFD'):
    char_to_image = char_to_image_dict(set(chars), font)

    return [(char_to_image[char]
             if char in char_to_image else
             char_to_image[unknown_char])
            for char in chars]

def filename_to_font(filename, size):
    return PIL.ImageFont.truetype(filename, size=size)

reason: import numpy


text:
import numpy as np
import sys
import torch
from collections import Counter
from skimage.metrics import structural_similarity as ssim
from sklearn.metrics import precision_score, recall_score, f1_score

sys.path.insert(0, './')

from utils.metric import Metric, MetricWithSimpleState

def _intersection_over_union_boxes(boxes_preds, boxes_labels, box_format="midpoint"):
    """
    Calculates intersection over union
    Parameters:
        boxes_preds (tensor): Predictions of Bounding Boxes (BATCH_SIZE, 4)
        boxes_labels (tensor): Correct Labels of Boxes (BATCH_SIZE, 4)
        box_format (str): midpoint/corners, if boxes (x,y,w,h) or (x1,y1,x2,y2)
    Returns:
        tensor: Intersection over union for all examples
    """

    # Slicing idx:idx+1 in order to keep tensor dimensionality
    # Doing ... in indexing if there would be additional dimensions
    # Like for Yolo algorithm which would have (N, S, S, 4) in shape
    if box_format == "midpoint":
        box1_x1 = boxes_preds[..., 0:1] - boxes_preds[..., 2:3] / 2
        box1_y1 = boxes_preds[..., 1:2] - boxes_preds[..., 3:4] / 2
        box1_x2 = boxes_preds[..., 0:1] + boxes_preds[..., 2:3] / 2
        box1_y2 = boxes_preds[..., 1:2] + boxes_preds[..., 3:4] / 2
        box2_x1 = boxes_labels[..., 0:1] - boxes_labels[..., 2:3] / 2
        box2_y1 = boxes_labels[..., 1:2] - boxes_labels[..., 3:4] / 2
        box2_x2 = boxes_labels[..., 0:1] + boxes_labels[..., 2:3] / 2
        box2_y2 = boxes_labels[..., 1:2] + boxes_labels[..., 3:4] / 2

    elif box_format == "corners":
        box1_x1 = boxes_preds[..., 0:1]
        box1_y1 = boxes_preds[..., 1:2]
        box1_x2 = boxes_preds[..., 2:3]
        box1_y2 = boxes_preds[..., 3:4]
        box2_x1 = boxes_labels[..., 0:1]
        box2_y1 = boxes_labels[..., 1:2]
        box2_x2 = boxes_labels[..., 2:3]
        box2_y2 = boxes_labels[..., 3:4]

    x1 = torch.max(box1_x1, box2_x1)
    y1 = torch.max(box1_y1, box2_y1)
    x2 = torch.min(box1_x2, box2_x2)
    y2 = torch.min(box1_y2, box2_y2)

    # Need clamp(0) in case they do not intersect, then we want intersection to be 0
    intersection = (x2 - x1).clamp(0) * (y2 - y1).clamp(0)
    box1_area = abs((box1_x2 - box1_x1) * (box1_y2 - box1_y1))
    box2_area = abs((box2_x2 - box2_x1) * (box2_y2 - box2_y1))

    return intersection / (box1_area + box2_area - intersection + 1e-6)

def _mAP_coco_evaluator(pred_boxes, targets, dataset):
    raise DeprecationWarning
    outputs = [{k: v.cpu() for k, v in t.items()} for t in pred_boxes]
    res = {target["image_id"].item(): output for target, output in zip(targets, outputs)}

    coco = get_coco_api_from_dataset(dataset)
    iou_types = ["bbox"]
    coco_evaluator = CocoEvaluator(coco, iou_types)

    coco_evaluator.update(res)

    coco_evaluator.synchronize_between_processes()
    coco_evaluator.accumulate()
    coco_evaluator.summarize()

    for v in coco_evaluator.coco_eval.values():
        return v.stats[0]

def _mAP(
        pred_boxes_batch, true_boxes_batch, iou_threshold=0.5, box_format="midpoint", num_classes=91,
        method="continuous"
):
    """
    Calculates mean average precision
    Parameters:
        pred_boxes (list): list of lists containing all bboxes with each bboxes
        specified as [train_idx, class_prediction, prob_score, x1, y1, x2, y2]
        true_boxes (list): Similar as pred_boxes except all the correct ones
        iou_threshold (float): threshold where predicted bboxes is correct
        box_format (str): "midpoint" or "corners" used to specify bboxes
        num_classes (int): number of classes
        method (str): 'interp' or 'continuous' for area under curve calculation
    Returns:
        float: mAP value across all classes given a specific IoU threshold
    """

    # pred_boxes = rcnn_2_boxes(pred_boxes)
    # true_boxes = rcnn_2_boxes(true_boxes)

    # used for numerical stability later on
    epsilon = 1e-6

    # list storing all mAP for respective batch idx
    map_batch = []
    for image_idx, (pred_boxes, true_boxes) in enumerate(zip(pred_boxes_batch, true_boxes_batch)):
        # list storing all AP for respective classes
        average_precisions = []
        for c in range(num_classes):
            detections = []
            ground_truths = []

            # Go through all predictions and targets,
            # and only add the ones that belong to the
            # current class c and current batch_idx
            for detection in pred_boxes:
                if detection[1] == c:
                    detections.append(detection)

            for true_box in true_boxes:
                if true_box[1] == c:
                    ground_truths.append(true_box)

            # find the amount of ground truth bboxes for the
            # training example and build empty tensor of length
            amount_bboxes = torch.zeros(len(ground_truths))

            # sort by box probabilities which is index 2
            detections.sort(key=lambda x: x[2], reverse=True)
            TP = torch.zeros((len(detections)))
            FP = torch.zeros((len(detections)))
            total_true_bboxes = len(ground_truths)

            # If none exists for this class then we can safely skip
            if total_true_bboxes == 0:
                continue

            for detection_idx, detection in enumerate(detections):
                best_iou = 0

                for idx, gt in enumerate(ground_truths):
                    iou = _intersection_over_union_boxes(
                        torch.tensor(detection[3:]),
                        torch.tensor(gt[3:]),
                        box_format=box_format,
                    )

                    if iou > best_iou:
                        best_iou = iou
                        best_gt_idx = idx

                if best_iou > iou_threshold:
                    # only detect ground truth detection once
                    if amount_bboxes[best_gt_idx] == 0:
                        # true positive and add this bounding box to seen
                        TP[detection_idx] = 1
                        amount_bboxes[best_gt_idx] = 1
                    else:
                        FP[detection_idx] = 1

                # if IOU is lower then the detection is a false positive
                else:
                    FP[detection_idx] = 1

            TP_cumsum = torch.cumsum(TP, dim=0)
            FP_cumsum = torch.cumsum(FP, dim=0)
            recalls = TP_cumsum / (total_true_bboxes + epsilon)
            precisions = TP_cumsum / (TP_cumsum + FP_cumsum + epsilon)
            precisions = torch.cat((torch.tensor([0]), precisions, torch.tensor([0])))
            recalls = torch.cat((torch.tensor([0]), recalls, torch.tensor([1])))
            # AP from area under curve
            if method == 'interp':
                ap = torch.trapz(precisions, recalls)
            else:
                i = torch.where(recalls[1:] - recalls[:-1])[0]  # points where x axis (recall) changes
                ap = torch.sum((recalls[i + 1] - recalls[i]) * precisions[i + 1])  # area under curve
            average_precisions.append(ap)
        mean_ap_image = sum(average_precisions) / (len(average_precisions) + 1e-10)
        map_batch.append(mean_ap_image.item() if not isinstance(mean_ap_image, float) else mean_ap_image)
    return map_batch

class mAP(MetricWithSimpleState):
    target_format_metric = 'boxes'

    @Metric.try_me
    @Metric.adapt_labels
    def __call__(self, *args, **kwargs):
        true_boxes, pred_boxes = args
        iou_threshold = kwargs.get('iou_threshold', self.params.get('iou_threshold', 0.5))
        box_format = kwargs.get('box_format', self.params.get('box_format', 'midpoint'))
        num_classes = kwargs.get('num_classes', self.params.get('num_classes', 91))
        method = kwargs.get('method', self.params.get('method', 'continuous'))
        batch = _mAP(pred_boxes, true_boxes, iou_threshold=iou_threshold, box_format=box_format,
                     num_classes=num_classes,
                     method=method)

        self._set_state(batch)
        return batch

    def __str__(self):
        return f'mAP[IoU={self.params.get("iou_threshold", 0.5):0.2f}]'

max_stars_repo_name: rickdr/Data-analysis-DNN-testing

max_issues_repo_name: rickdr/Data-analysis-DNN-testing

max_forks_repo_name: rickdr/Data-analysis-DNN-testing


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import lightgbm as lgb

from sklearn.model_selection import KFold
from sklearn.metrics import roc_auc_score
from sklearn.preprocessing import LabelEncoder

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

import gc

def kde_target(var_name, df):
    """
    Plot the correlation between the target and given variable.
    :param var_name: String, the name of the given variable.
    :param df: Dataframe that contains both the given variable and the 'TARGET' column.
    :return:
    """
    # Calculate the correlation coeffienct between the new variable and the target
    corr = df['TARGET'].corr(df[var_name])

    # Calculate medians for repaid vs not repaid
    avg_repaid = df.loc[df['TARGET'] == 0, var_name].median()
    avg_not_repaid = df.loc[df['TARGET'] == 1, var_name].median()

    plt.figure(figsize=(12, 6))

    # Plot the distribution for target == 0 and target == 1
    sns.kdeplot(df.loc[df['TARGET'] == 0, var_name].dropna(), label='TARGET == 0')
    sns.kdeplot(df.loc[df['TARGET'] == 1, var_name].dropna(), label='TARGET == 1')

    # Lable the plot
    plt.xlabel(var_name)
    plt.ylabel('Density')
    plt.title('%s Distribution' % var_name)
    plt.legend()

    # Print out the correlation
    print('The correlation between %s and the TARGET is %0.4f' % (var_name, corr))
    # Print out average values
    print('Median value for loan that was not repaid =\t %0.4f' % avg_not_repaid)
    print('Median value for loan that was repaid =\t %0.4f' % avg_repaid)

def agg_numeric(df, group_var, df_name):
    """Aggregates the numeric values in a dataframe. This can
    be used to create features for each instance of the grouping variable.

    Parameters
    --------
        df (dataframe):
            the dataframe to calculate the statistics on
        group_var (string):
            the variable by which to group df
        df_name (string):
            the variable used to rename the columns

    Return
    --------
        agg (dataframe):
            a dataframe with the statistics aggregated for
            all numeric columns. Each instance of the grouping variable will have
            the statistics (mean, min, max, sum; currently supported) calculated.
            The columns are also renamed to keep track of features created.

    """
    # Remove id variables other than grouping variable
    for col in df:
        if col != group_var and 'SK_ID' in col:
            df = df.drop(columns=col)

    group_ids = df[group_var]
    numeric_df = df.select_dtypes('number')
    numeric_df[group_var] = group_ids

    # Group by the specified variable and calculate the statistics
    agg = numeric_df.groupby(group_var).agg(['count', 'mean', 'max', 'min', 'sum']).reset_index()

    # Need to create new column names
    columns = [group_var]

    # Iterate through the variables names
    for var in agg.columns.levels[0]:
        # Skip the grouping variable
        if var != group_var:
            # Iterate through the stat names
            for stat in agg.columns.levels[1][:-1]:
                # Make a new column name for the variable and stat
                columns.append('%s_%s_%s' % (df_name, var, stat))

    agg.columns = columns
    return agg

def count_categorical(df, group_var, df_name):
    """Computes counts and normalized counts for each observation
    of `group_var` of each unique category in every categorical variable

    Parameters
    --------
    df : dataframe
        The dataframe to calculate the value counts for.

    group_var : string
        The variable by which to group the dataframe. For each unique
        value of this variable, the final dataframe will have one row

    df_name : string
        Variable added to the front of column names to keep track of columns

    Return
    --------
    categorical : dataframe
        A dataframe with counts and normalized counts of each unique category in every categorical variable
        with one row for every unique value of the `group_var`.

    """

    # Select the categorical columns
    categorical = pd.get_dummies(df.select_dtypes('object'))

    # Make sure to put the identifying id on the column
    categorical[group_var] = df[group_var]

    # Groupby the group var and calculate the sum and mean
    categorical = categorical.groupby(group_var).agg(['sum', 'mean'])

    column_names = []

    # Iterate through the columns in level 0
    for var in categorical.columns.levels[0]:
        # Iterate through the stats in level 1
        for stat in ['count', 'count_norm']:
            # Make a new column name
            column_names.append('%s_%s_%s' % (df_name, var, stat))

    categorical.columns = column_names

    return categorical

def missing_values_table(df):
    """
    Function to calculate missing values by column
    :param df:
    :return:
    """
    # Total missing values
    mis_val = df.isnull().sum()

    # Percentage of missing values
    mis_val_percent = 100 * df.isnull().sum() / len(df)

    # Make a table with the results
    mis_val_table = pd.concat([mis_val, mis_val_percent], axis=1)

    # Rename the columns
    mis_val_table_ren_columns = mis_val_table.rename(
        columns={0: 'Missing Values', 1: '% of Total Values'})

    # Sort the table by percentage of missing descending
    mis_val_table_ren_columns = mis_val_table_ren_columns[
        mis_val_table_ren_columns.iloc[:, 1] != 0].sort_values(
        '% of Total Values', ascending=False).round(1)

    # Print some summary information
    print("Your selected dataframe has " + str(df.shape[1]) + " columns.\n"
                                                              "There are " + str(mis_val_table_ren_columns.shape[0]) +
          " columns that have missing values.")

    # Return the dataframe with missing information
    return mis_val_table_ren_columns

def plot_feature_importances(df):
    """
        Plot importances returned by a model. This can work with any measure of
        feature importance provided that higher importance is better.

        Args:
            df (dataframe): feature importances. Must have the features in a column
            called `features` and the importances in a column called `importance

        Returns:
            shows a plot of the 15 most importance features

            df (dataframe): feature importances sorted by importance (highest to lowest)
            with a column for normalized importance
    """

    # Sort features according to importance
    df = df.sort_values('importance', ascending= False).reset_index()

    # Normalize the feature importances to add up to one
    df['importance_normalized'] = df['importance'] / df['importance'].sum()

    # Make a horizontal bar chart of feature importances
    plt.figure(figsize=(10, 6))
    ax = plt.subplot()

    # Need to reverse the index to plot most important on top
    ax.barh(list(reversed(list(df.index[:15]))),
            df['importance_normalized'].head(15),
            align='center', edgecolor='k')

    # Set the yticks and labels
    ax.set_yticks(list(reversed(list(df.index[:15]))))
    ax.set_yticklabels(df['feature'].head(15))

    # Plot labeling
    plt.xlabel('Normalized Importance')
    plt.title('Feature Importances')
    plt.show()

    return df


text:
# -*- coding: utf-8 -*-
"""
model_440_basicChSa.py: erosion model using depth-dependent cubic diffusion
with a soil layer, basic stream power, and discharge proportional to drainage
area.

Model 440 BasicChSa

Landlab components used: FlowRouter, DepressionFinderAndRouter,
                         FastscapeStreamPower, DepthDependentCubicDiffuser,
                         ExponentialWeatherer

@author: gtucker
@author: Katherine Barnhart
"""

from terrainbento.base_class import ErosionModel
from landlab.components import (FlowAccumulator, DepressionFinderAndRouter,
                                FastscapeEroder, DepthDependentTaylorDiffuser,
                                ExponentialWeatherer)
import numpy as np

class BasicChSa(ErosionModel):
    """
    A BasicChSa model computes erosion using depth-dependent cubic diffusion
    with a soil layer, basic stream power, and Q~A.
    """

    def __init__(self, input_file=None, params=None,
                 BaselevelHandlerClass=None):
        """Initialize the BasicChSa model."""

        # Call ErosionModel's init
        super(BasicChSa, self).__init__(input_file=input_file,
                                        params=params,
                                        BaselevelHandlerClass=BaselevelHandlerClass)

        self.K_sp = self.get_parameter_from_exponent('K_sp')
        linear_diffusivity = (self._length_factor**2.)*self.get_parameter_from_exponent('linear_diffusivity') # has units length^2/time
        try:
            initial_soil_thickness = (self._length_factor)*self.params['initial_soil_thickness'] # has units length
        except KeyError:
            initial_soil_thickness = 1.0  # default value
        soil_transport_decay_depth = (self._length_factor)*self.params['soil_transport_decay_depth']  # has units length
        max_soil_production_rate = (self._length_factor)*self.params['max_soil_production_rate'] # has units length per time
        soil_production_decay_depth = (self._length_factor)*self.params['soil_production_decay_depth']   # has units length

        # Create soil thickness (a.k.a. depth) field
        if 'soil__depth' in self.grid.at_node:
            soil_thickness = self.grid.at_node['soil__depth']
        else:
            soil_thickness = self.grid.add_zeros('node', 'soil__depth')

        # Create bedrock elevation field
        if 'bedrock__elevation' in self.grid.at_node:
            bedrock_elev = self.grid.at_node['bedrock__elevation']
        else:
            bedrock_elev = self.grid.add_zeros('node', 'bedrock__elevation')

        soil_thickness[:] = initial_soil_thickness
        bedrock_elev[:] = self.z - initial_soil_thickness

        # Instantiate a FlowAccumulator with DepressionFinderAndRouter using D8 method
        self.flow_router = FlowAccumulator(self.grid,
                                           flow_director='D8',
                                           depression_finder = DepressionFinderAndRouter)

        # Instantiate a FastscapeEroder component
        self.eroder = FastscapeEroder(self.grid,
                                      K_sp=self.K_sp,
                                      m_sp=self.params['m_sp'],
                                      n_sp=self.params['n_sp'])

        # Instantiate a weathering component
        self.weatherer = ExponentialWeatherer(self.grid,
                                              max_soil_production_rate=max_soil_production_rate,
                                              soil_production_decay_depth=soil_production_decay_depth)

        # Instantiate a soil-transport component
        self.diffuser = DepthDependentTaylorDiffuser(self.grid,
                                                    linear_diffusivity=linear_diffusivity,
                                                    slope_crit=self.params['slope_crit'],
                                                    soil_transport_decay_depth=soil_transport_decay_depth,
                                                    nterms=11)

    def run_one_step(self, dt):
        """
        Advance model for one time-step of duration dt.
        """

        # Route flow
        self.flow_router.run_one_step()

        # Get IDs of flooded nodes, if any
        flooded = np.where(self.flow_router.depression_finder.flood_status==3)[0]

        # Do some erosion (but not on the flooded nodes)
        # (if we're varying K through time, update that first)
        if self.opt_var_precip:
            self.eroder.K = (self.K_sp
                             * self.pc.get_erodibility_adjustment_factor(self.model_time))

        self.eroder.run_one_step(dt, flooded_nodes=flooded)

        # We must also now erode the bedrock where relevant. If water erosion
        # into bedrock has occurred, the bedrock elevation will be higher than
        # the actual elevation, so we simply re-set bedrock elevation to the
        # lower of itself or the current elevation.
        b = self.grid.at_node['bedrock__elevation']
        b[:] = np.minimum(b, self.grid.at_node['topographic__elevation'])

        # Calculate regolith-production rate
        self.weatherer.calc_soil_prod_rate()

        # Do some soil creep
        self.diffuser.run_one_step(dt,
                                   dynamic_dt=True,
                                   if_unstable='raise',
                                   courant_factor=0.1)

        # calculate model time
        self.model_time += dt

        # Lower outlet
        self.update_outlet(dt)

        # Check walltime
        self.check_walltime()

def main():
    """Executes model."""
    import sys

    try:
        infile = sys.argv[1]
    except IndexError:
        print('Must include input file name on command line')
        sys.exit(1)

    cdsp = BasicChSa(input_file=infile)
    cdsp.run()

if __name__ == '__main__':
    main()

max_stars_repo_name: mcflugen/terrainbento

max_issues_repo_name: mcflugen/terrainbento

max_forks_repo_name: mcflugen/terrainbento


text:
import os
import glob
import shutil
import math
import logging
import numpy as np
import pandas as pd
import xml.etree.ElementTree as etree

import sphere.IRDIS as IRDIS
import sphere.IFS as IFS
import sphere.SPARTA as SPARTA

from pathlib import Path
from astropy.io import fits
from astropy.time import Time

_log = logging.getLogger(__name__)

def process_mainFiles(mainFiles, files, logger=_log):
    '''
    Process top-level file association XML from the ESO archive

    Parameters
    ----------
    mainFiles : etree element
        Main file element

    files : list
        List where the files will be appended

    logger : logHandler object
        Log handler for the reduction. Default is root logger

    '''
    # append file names to the list
    for file in mainFiles:
        fname = file.attrib['name']
        files.append(fname)

        logger.debug(' ==> {0}'.format(fname))

def process_association(tree, files, logger=_log):
    '''
    Process file association XML from the ESO archive

    Parameters
    ----------
    tree : etree element
        Associated file element

    files : list
        List where the files will be appended

    logger : logHandler object
        Log handler for the reduction. Default is root logger

    '''
    catg = tree.attrib['category']

    logger.debug(catg)

    # skip unused calibrations
    if (catg == 'IFS_STD_ASTROM') or (catg == 'IFS_STD_PHOT') or \
       (catg == 'IFS_DIST') or (catg == 'IRD_CLI_PHOT') or \
       (catg == 'IRD_DIST'):
        logger.debug(' ==> skipping')
        return

    # process differently mainFiles from associatedFiles
    for elt in tree:
        if elt.tag == 'mainFiles':
            logger.debug('mainFiles')
            process_mainFiles(elt, files)
        elif elt.tag == 'associatedFiles':
            logger.debug('associatedFiles')
            for nelt in elt:
                process_association(nelt, files, logger=logger)

def sort_files_from_xml(path, logger=_log):
    '''Sort files downloaded from the ESO archive with associated raw
       calibrations

    When choosing this option from the archive, the data package
    includes xml files that provide the association between the
    science files and the calibration files. The method uses these
    xml files to copy the science and associated raw files into
    dedicated directories.

    For a given executed OB, the method will extract the target
    name from the OBS.NAME keyword and the OB ID from the OBS.ID
    keyword. It will create a directory {OBS.NAME}_id={OBS.ID} and
    copy the FITS files. Note that any space in the name will be
    replaced by an underscore to avoid subsequent issues with
    esorex in the data reduction.

    The files are *copied* because science data acquired on the
    same night with the same setup can share some calibrations. It
    is therefore safer to copy the needed fiels instead of moving
    them. At the end, all the original files are more to the
    all_files/ subdirectory.

    Parameters
    ----------
    path : str
        Path where to look for XML files

    logger : logHandler object
        Log handler for the reduction. Default is root logger

    '''

    xml_files = list(path.glob('*.xml'))

    logger.info('Sort data based on XML files (ESO automated calibration selection)')
    logger.info(' * {0} XML files\n'.format(len(xml_files)))

    # sort files
    for file in xml_files:
        tree = etree.parse(file)
        root = tree.getroot()

        logger.info(' * {}'.format(file.name))

        # process only IFS and IRDIS science data
        catg = root.attrib['category']
        if (catg.find('ACQUISITION') != -1):
            continue

        # get target name from first mainFile element
        scifiles = root.find('mainFiles')
        filename = scifiles[0].attrib['name']

        # Mac OS X replaces : by _ in file names...
        if not (path / '{}.fits'.format(filename)).exists():
            filename = filename.replace(':', '_')

        if not (path / '{}.fits'.format(filename)).exists():
            logger.info('   ==> file {} does not exist. Skipping'.format(filename))
            continue

        fpath = path / '{}.fits'.format(filename)
        hdr = fits.getheader(fpath)

        # target and arm
        target = hdr['HIERARCH ESO OBS NAME']
        obs_id = hdr['HIERARCH ESO OBS ID']
        mjd    = Time(math.floor(float(hdr['MJD-OBS']+0.5))-0.5, format='mjd')
        night  = mjd.isot[:10]
        if catg == 'SCIENCE_OBJECT_AO':
            instrument = 'SPARTA'
        else:
            try:
                arm = hdr['HIERARCH ESO SEQ ARM']
            except KeyError:
                logger.error('No \'HIERARCH ESO SEQ ARM\' keyword in {}'.format(fpath))
                continue

            if arm == 'IRDIS':
                instrument = 'IRDIS'
            elif arm == 'IFS':
                instrument = 'IFS'
            else:
                logger.error('Unknown arm {0}'.format(arm))
                continue

        # get files
        files = []
        process_association(root, files, logger=logger)

        # target path
        directory = '{0}_id={1}'.format(target, obs_id)
        directory = '_'.join(directory.split())
        target_path = path / directory / night / instrument / 'raw'
        target_path.mkdir(parents=True, exist_ok=True)

        # copy files
        for filename in files:
            fpath = path / '{}.fits'.format(filename)

            # Mac OS X replaces : by _ in file names...
            if not fpath.exists():
                filename = filename.replace(':', '_')
                fpath  = path / '{}.fits'.format(filename)

            # check if file actually exists
            if not fpath.exists():
                logger.info(' ==> file {} does not exist. Skipping.'.format(fpath))
                continue

            # copy if needed
            nfpath = target_path / '{}.fits'.format(filename)
            if not nfpath.exists():
                shutil.copy(fpath, nfpath)

        # print status
        logger.debug('{0} - id={1}'.format(target, obs_id))
        logger.debug(' ==> found {0} files'.format(len(files)))
        logger.debug(' ==> copied to {0}'.format(target_path))

    # move all files
    path_new = path / 'all_files'
    path_new.mkdir(parents=True, exist_ok=True)

    files = []
    files.extend(list(path.glob('*.fits')))
    files.extend(list(path.glob('*.xml')))
    files.extend(list(path.glob('*.txt')))

    if len(files) != 0:
        for file in files:
            file.rename(path_new / file.name)

def sort_files_from_fits(path, logger=_log):
    '''Sort FITS files based only based on their headers

    Contrary to sort_files_from_xml(), this method is dumb in the
    sense that it does not try to keep any file associated. It
    just sorts the FITS files as a function of the
    sub-system. Still very convenient when manually downloading
    data and calibrations.

    And contrary to sort_files_from_xml(), the files are directly
    moved to their sub-system directory. This makes this method
    much faster than sort_files_from_xml(). At the end, all
    unsorted files are moved to an unsorted_files/ subdirectory.

    Parameters
    ----------
    path : str
        Path where to look for FITS files

    logger : logHandler object
        Log handler for the reduction. Default is root logger

    '''

    fits_files = list(path.glob('*.fits'))

    logger.info('Sort data based on FITS files')
    logger.info(' * {0} FITS files\n'.format(len(fits_files)))

    # sort files
    for file in fits_files:
        logger.info(' * {}'.format(file.name))

        # target and arm
        hdr = fits.getheader(file)

        try:
            target = hdr['HIERARCH ESO OBS NAME']
            obs_id = hdr['HIERARCH ESO OBS ID']
            dpr_type = hdr['HIERARCH ESO DPR TYPE']
        except KeyError:
            logger.error('Missing ESO HIERARCH keywords in {}'.format(file))
            continue

        if dpr_type == 'OBJECT,AO':
            instrument = 'SPARTA'
        else:
            try:
                arm = hdr['HIERARCH ESO SEQ ARM']
            except KeyError:
                logger.error('No \'HIERARCH ESO SEQ ARM\' keyword in {}'.format(file))
                continue

            if arm == 'IRDIS':
                instrument = 'IRDIS'
            elif arm == 'IFS':
                instrument = 'IFS'
            else:
                logger.error('Unknown arm {0}'.format(arm))
                continue

        # target path
        target_path = path / instrument / 'raw'
        target_path.mkdir(parents=True, exist_ok=True)

        # move file
        file.rename(target_path / file.name)

        # print status
        logger.debug('{0} - id={1}'.format(target, obs_id))
        logger.debug(' ==> copied to {0}'.format(target_path))

    # move all files
    path_new = path / 'unsorted_files'
    path_new.mkdir(parents=True, exist_ok=True)

    files = []
    files.extend(list(path.glob('*.fits')))
    files.extend(list(path.glob('*.txt')))

    if len(files) != 0:
        for file in files:
            file.rename(path_new / file.name)

def classify_irdis_dataset(path, logger=_log):
    '''Classify an IRDIS dataset based on the science files

    Parameters
    ----------
    path : str
        Path to the directory containing the dataset

    logger : logHandler object
        Log handler for the reduction. Default is root logger

    Returns
    -------
    mode : str
        Generic string representing the name of the mode. None in case
        of failure.

    '''

    # zeroth-order reduction validation
    raw = path / 'raw'
    if not raw.exists():
        logger.error('No raw/ subdirectory. {0} is not a valid reduction path!'.format(path))
        return None

    # list all fits files
    files = list(raw.glob('*.fits'))
    if len(files) == 0:
        return None

    # search for science files
    modes = []
    for file in files:
        hdr = fits.getheader(file)

        dpr_catg = hdr.get('HIERARCH ESO DPR CATG')
        mode     = hdr.get('HIERARCH ESO INS1 MODE')

        if dpr_catg == 'SCIENCE':
            modes.append(mode)

    modes = np.array(modes)

    n_imaging = (modes == 'NIROBS').sum() + (modes == 'EXT').sum() + \
        (modes == 'DBI').sum() + (modes == 'CI').sum()
    n_pola    = (modes == 'DPI').sum()
    n_spectro = (modes == 'LSS').sum()

    if (n_imaging >= n_pola) and (n_imaging >= n_spectro):
        return 'imaging'
    elif (n_pola >= n_imaging) and (n_pola >= n_spectro):
        return 'pola'
    else:
        return 'spectro'

class Dataset:
    '''
    A SPHERE dataset for a given object
    '''

    ######################
    # Constructor
    ######################
    def __init__(self, path, log_level='info'):
        '''
        Initialization code for a SPHERE dataset

        Parameters
        ----------
        path : str
            Path to the SPHERE data
        '''

        if not isinstance(path, str):
            raise ValueError('path must be a string')

        # path
        path = Path(path).expanduser().resolve()
        self._path = path

        # configure logging
        logger = logging.getLogger(str(path))
        logger.setLevel(log_level.upper())
        if logger.hasHandlers():
            for hdlr in logger.handlers:
                logger.removeHandler(hdlr)

        handler = logging.FileHandler(self.path / 'dataset.log', mode='w', encoding='utf-8')
        formatter = logging.Formatter('%(asctime)s\t%(levelname)8s\t%(message)s')
        formatter.default_msec_format = '%s.%03d'
        handler.setFormatter(formatter)
        logger.addHandler(handler)

        self._log_level = log_level
        self._handler   = handler
        self._logger    = logger

        self._logger.info('Looking for SPHERE data sets at path {}'.format(path))

        # list of reductions
        self._IFS_reductions    = []
        self._IRDIS_reductions  = []
        self._SPARTA_reductions = []

        # search for data with calibrations downloaded from ESO archive
        xml_files = list(path.glob('*.xml'))
        if len(xml_files) != 0:
            sort_files_from_xml(path, logger=self._logger)

        # directly search for data
        fits_files = list(path.glob('*.fits'))
        if len(fits_files) != 0:
            sort_files_from_fits(path, logger=self._logger)

        # recursively look for valid reduction
        self._create_reductions()

    ##################################################
    # Representation
    ##################################################

    def __repr__(self):
        return '<SPHERE datasets: {0} IFS, {1} IRDIS, {2} SPARTA>'.format(len(self.IFS_reductions), len(self.IRDIS_reductions), len(self.SPARTA_reductions))

    ##################################################
    # Properties
    ##################################################

    @property
    def reductions(self):
        return self._reductions

    @property
    def IRDIS_reductions(self):
        return self._IRDIS_reductions

    @property
    def IFS_reductions(self):
        return self._IFS_reductions

    @property
    def SPARTA_reductions(self):
        return self._SPARTA_reductions

    @property
    def path(self):
        return self._path

    ##################################################
    # Generic class methods
    ##################################################

    def init_reduction(self):
        '''
        Sort files and frames, perform sanity check
        '''

        for r in self.reductions:
            self._logger.info('Init: {}'.format(str(r)))

            r.init_reduction()

    def create_static_calibrations(self):
        '''
        Create static calibrations with esorex
        '''

        for r in self.reductions:
            self._logger.info('Static calibrations: {}'.format(str(r)))

            r.create_static_calibrations()

    def preprocess_science(self):
        '''
        Clean and collapse images
        '''

        for r in self.reductions:
            self._logger.info('Science pre-processing: {}'.format(str(r)))

            r.preprocess_science()

    def process_science(self):
        '''
        Perform star center, combine cubes into final (x,y,time,lambda)
        cubes, correct anamorphism and scale the images
        '''

        for r in self.reductions:
            self._logger.info('Science processing: {}'.format(str(r)))

            r.process_science()

    def clean(self):
        '''
        Clean the reduction directory, leaving only the raw and products
        sub-directory
        '''

        for r in self.reductions:
            print(r)
            self._logger.info('Clean-up: {}'.format(str(r)))

            r.clean()

    def full_reduction(self):
        '''
        Performs a full reduction of a data set, from the static
        calibrations to the final (x,y,time,lambda) cubes
        '''

        for r in self.reductions:
            self._logger.info('###########################################################################')
            self._logger.info('# Full reduction: {}'.format(str(r)))
            self._logger.info('###########################################################################')

            r.full_reduction()

    ##################################################
    # Class methods
    ##################################################

    def _create_reductions(self):
        '''
        Detect and create valid reductions in path
        '''

        self._logger.info('Create reductions from sorted data')

        wpath = os.walk(self.path)
        for w in wpath:
            subs = w[1]
            if 'raw' in subs:
                # if directory has a raw/ sub-directory, make sure it
                # has FITS files and that they are from a valid
                # sub-system
                reduction_path = Path(w[0])
                fits_files = list((reduction_path / 'raw').glob('*.fits'))
                if len(fits_files) != 0:
                    hdr = fits.getheader(fits_files[0])

                    if hdr['HIERARCH ESO DPR TYPE'] == 'OBJECT,AO':
                        arm = 'SPARTA'
                    else:
                        try:
                            arm = hdr['HIERARCH ESO SEQ ARM']
                        except KeyError:
                            self._logger.error('No \'HIERARCH ESO SEQ ARM\' keyword in {}'.format(fits_files[0]))

                    if arm == 'IRDIS':
                        mode = classify_irdis_dataset(reduction_path, logger=self._logger)

                        # an error occured in dataset classification
                        if mode is None:
                            continue

                        if mode == 'imaging':
                            self._logger.info(' * IRDIS imaging reduction at path {}'.format(reduction_path))
                            reduction  = IRDIS.ImagingReduction(reduction_path, log_level=self._log_level,
                                                                sphere_handler=self._handler)
                        elif mode == 'polar':
                            self._logger.warning('IRDIS DPI not supported yet')
                        elif mode == 'spectro':
                            self._logger.info(' * IRDIS spectro <reduction at path {}'.format(reduction_path))
                            reduction  = IRDIS.SpectroReduction(reduction_path, log_level=self._log_level,
                                                                sphere_handler=self._handler)

                        # save if reduction was successfully created
                        if reduction is not None:
                            self.IRDIS_reductions.append(reduction)
                    elif arm == 'IFS':
                        self._logger.info(' * IFS reduction at path {}'.format(reduction_path))
                        reduction  = IFS.Reduction(reduction_path, log_level=self._log_level,
                                                   sphere_handler=self._handler)

                        # save if reduction was successfully created
                        if reduction is not None:
                            self.IFS_reductions.append(reduction)
                    elif arm == 'SPARTA':
                        self._logger.info(' * SPARTA reduction at path {}'.format(reduction_path))
                        reduction  = SPARTA.Reduction(reduction_path, log_level=self._log_level,
                                                      sphere_handler=self._handler)

                        # save if reduction was successfully created
                        if reduction is not None:
                            self.SPARTA_reductions.append(reduction)
                    else:
                        self._logger.error('Unknown arm {0}'.format(arm))
                        continue

                    # self._logger.info(reduction_path)
                    self._logger.info('   ==> {} files'.format(len(fits_files)))

        # merge all reductions into a single list
        self._reductions = self.IFS_reductions + self.IRDIS_reductions

max_stars_repo_name: tomasstolker/SPHERE

max_issues_repo_name: tomasstolker/SPHERE

max_forks_repo_name: tomasstolker/SPHERE

reason: import numpy,from astropy


text:
from os import PathLike
from typing import Callable, Union

import numpy as np

PathOrStr = Union[PathLike, str]
ScoreFunction = Callable[[np.ndarray, np.ndarray], float]

max_stars_repo_name: agkphysics/emotion

max_issues_repo_name: agkphysics/emotion

max_forks_repo_name: agkphysics/emotion


text:
import glob
import json
import os
import sys
import numpy as np
import re
from collections import OrderedDict as od
import operator
from keras.datasets import imdb
import random
import pickle
from bs4 import BeautifulSoup
from keras.preprocessing.text import Tokenizer, text_to_word_sequence
from keras.utils.np_utils import to_categorical
import shutil
from nltk.stem.porter import PorterStemmer
import nltk

pattern = re.compile(r'\b(' + r'|'.join(nltk.corpus.stopwords.words('english')) + r')\b\s*')
keep_pattern = re.compile('([\W])')

def clean_str(string):
    """
    Tokenization/string cleaning for dataset
    Every dataset is lower cased except
    """
    string = re.sub(r"\\", "", string)
    string = re.sub(r"\'", "", string)
    string = re.sub(r"\"", "", string)
    return string.strip().lower()

def split_data(X_train,y_train,X_test,y_test, word_index,
               MAX_NB_WORDS, MAX_SENTS, MAX_SENT_LENGTH,
                singlefn, BATCH_TRAIN_SIZE=20, BATCH_TEST_SIZE=20):

    dir, fname = os.path.split(singlefn)

    newdir=dir+'/'+fname.split('.')[0]
    print(newdir)

    if not os.path.isdir(newdir):
        os.mkdir(newdir)

    train_dir=newdir+'/train'
    print(train_dir)
    if not os.path.isdir(train_dir):
        os.mkdir(train_dir)
    else:
        shutil.rmtree(train_dir)
        os.mkdir(train_dir)

    test_dir = newdir + '/test'
    print(test_dir)
    if not os.path.isdir(test_dir):
        os.mkdir(test_dir)
    else:
        shutil.rmtree(test_dir)
        os.mkdir(test_dir)
    num_sample_train=X_train.shape[0]
    num=0

    print('dump train data...')

    while num<num_sample_train:
        bfxname=train_dir+'/Xy_train.{}.pkl'.format(num)
        pickle.dump((X_train[num:num+BATCH_TRAIN_SIZE],y_train[num:num+BATCH_TRAIN_SIZE]), open(bfxname,'wb'))
        num+=BATCH_TRAIN_SIZE

    num=0
    num_sample_test = X_test.shape[0]

    print('dump test data')

    while num<num_sample_test:
        bfxname=test_dir+'/Xy_test.{}.pkl'.format(num)
        pickle.dump((X_test[num:num+BATCH_TEST_SIZE],y_test[num:num+BATCH_TEST_SIZE]), open(bfxname,'wb'))
        num += BATCH_TEST_SIZE

    pickle.dump((word_index, MAX_NB_WORDS, MAX_SENTS, MAX_SENT_LENGTH,
                 BATCH_TRAIN_SIZE, BATCH_TEST_SIZE,
                 num_sample_train, num_sample_test),
                open(newdir+'/meta.pkl','wb'))

    print('done!')

def preprocess_rawbig_imdb_ha(MAX_NB_WORDS, MAX_SENTS, MAX_SENT_LENGTH,
                              inputfn='./data/data.json',
                              output_dir='./data/big_imdb_prep_ha50200.pkl'):
    objs=json.load(open(inputfn))
    X=[]
    y=[]
    for docobj in objs:
        X.append(docobj['review'])
        y.append(docobj['rating']-1)
    print(len(X))
    print(X[:10])
    print(y[:10])

    num_sam = len(X)
    num_train = int(num_sam * 0.9)
    index_shuf = list(range(num_sam))
    print(index_shuf[:10])
    random.shuffle(index_shuf)
    print(index_shuf[:10])
    print(len(index_shuf))
    X_train=[]
    y_train=[]
    for i in index_shuf:
        X_train.append(X[i])
        y_train.append(y[i])

    X_test = X_train[num_train:]
    X_train = X_train[:num_train]
    y_test = y_train[num_train:]
    y_train = y_train[:num_train]

    from nltk import tokenize

    reviews = []
    labels = []
    texts = []
    p_stemmer = PorterStemmer()
    print('start get train text...')
    for idx in range(len(X_train)):
        text = BeautifulSoup(X_train[idx], "html.parser")
        text = clean_str(text.get_text())
        # text = pattern.sub('', text)
        # sentences = tokenize.sent_tokenize(text)
        tokens = re.split(keep_pattern, text)
        # stemmed_tokens = [p_stemmer.stem(i)  for i in tokens if i is not ' ' and i is not '']
        text = ' '.join(tokens).strip()
        texts.append(text)
        sentences = tokenize.sent_tokenize(text)
        reviews.append(sentences)
        if idx % 1000 == 0:
            print('done {}'.format(idx))

        labels.append(y_train[idx])

    print(reviews[0][:10])
    print(reviews[1][:10])

    print('number of docs: {}'.format(len(reviews)))
    print('number of labels: {}'.format(len(labels)))

    tokenizer = Tokenizer(num_words=MAX_NB_WORDS, filters='')
    tokenizer.fit_on_texts(texts)

    data = np.zeros((len(texts), MAX_SENTS, MAX_SENT_LENGTH), dtype='int32')
    test_add_word_count = 0
    word2ind = {}
    list_num_sen = []
    list_num_word = []
    for i, sentences in enumerate(reviews):
        avg_sen = 0
        for j, sent in enumerate(sentences):
            avg_sen += 1
            if j < MAX_SENTS:
                wordTokens = text_to_word_sequence(sent)
                k = 1
                list_num_word.append(len(wordTokens))
                for _, word in enumerate(reversed(wordTokens)):
                    if k < MAX_SENT_LENGTH and tokenizer.word_index.get(word) \
                            and tokenizer.word_index[word] < MAX_NB_WORDS:
                        data[i, j, -k] = tokenizer.word_index[word]
                        word2ind[word] = tokenizer.word_index[word]
                        test_add_word_count += 1
                        k = k + 1
        list_num_sen.append(avg_sen)

    word_index = tokenizer.word_index
    print('Total %s unique tokens.' % len(word_index))
    print('num word add count {}'.format(test_add_word_count))
    print('num word in limit text {}'.format(len(word2ind)))

    print('=====some words======')
    print((sorted(word2ind.items(), key=operator.itemgetter(1)))[:10])

    print('sen statistic')
    print(min(list_num_sen))
    print(max(list_num_sen))
    print(sum(list_num_sen) / (1.0 * len(list_num_sen)))

    print('words statistic')
    print(min(list_num_word))
    print(max(list_num_word))
    print(sum(list_num_word) / (1.0 * len(list_num_word)))

    labels = to_categorical(np.asarray(labels))
    print('Shape of data tensor:', data.shape)
    print('Shape of label tensor:', labels.shape)

    data_train = data
    labels_train = labels

    reviews = []
    labels = []
    texts = []

    print('start get test text...')
    for idx in range(len(X_test)):
        text = BeautifulSoup(X_test[idx], "html.parser")
        text = clean_str(text.get_text())
        # text = pattern.sub('', text)
        # texts.append(text)
        # sentences = tokenize.sent_tokenize(text)
        tokens = re.split(keep_pattern, text)
        # stemmed_tokens = [p_stemmer.stem(i) for i in tokens if i is not ' ' and i is not '']
        text = ' '.join(tokens).strip()
        texts.append(text)
        sentences = tokenize.sent_tokenize(text)
        reviews.append(sentences)

        if idx % 1000 == 0:
            print('done {}'.format(idx))

        labels.append(y_test[idx])

    print(reviews[0][:10])
    print(reviews[1][:10])

    print('number of docs: {}'.format(len(reviews)))
    print('number of labels: {}'.format(len(labels)))

    data = np.zeros((len(texts), MAX_SENTS, MAX_SENT_LENGTH), dtype='int32')
    test_add_word_count = 0
    strange_words = set()
    list_num_sen = []
    list_num_word = []
    for i, sentences in enumerate(reviews):
        avg_sen = 0
        for j, sent in enumerate(sentences):
            avg_sen += 1
            if j < MAX_SENTS:
                wordTokens = text_to_word_sequence(sent)
                k = 1
                list_num_word.append(len(wordTokens))
                for _, word in enumerate(reversed(wordTokens)):
                    if k < MAX_SENT_LENGTH and word2ind.get(word):
                        data[i, j, -k] = tokenizer.word_index[word]
                        test_add_word_count += 1
                        k = k + 1
                    if not word2ind.get(word):
                        strange_words.add(word)

        list_num_sen.append(avg_sen)

    print('num word add count {}'.format(test_add_word_count))
    print('num word strange {}'.format(len(strange_words)))

    print('sen statistic')
    print(min(list_num_sen))
    print(max(list_num_sen))
    print(sum(list_num_sen) / (1.0 * len(list_num_sen)))

    print('words statistic')
    print(min(list_num_word))
    print(max(list_num_word))
    print(sum(list_num_word) / (1.0 * len(list_num_word)))

    labels = to_categorical(np.asarray(labels))
    print('Shape of data tensor:', data.shape)
    print('Shape of label tensor:', labels.shape)

    data_test = data
    labels_test = labels

    print('=======================')
    print(data_train.shape)
    print(labels_train.shape)
    print(data_test.shape)
    print(labels_test.shape)

    print('start dump....')
    split_data(data_train, labels_train,data_test, labels_test,
                  word2ind,MAX_NB_WORDS,MAX_SENTS,MAX_SENT_LENGTH,
                  output_dir,20,20)
    print('done!!!')

def gen_imdb_tf_from_obj(X_train, y_train, X_test, y_test,
                         train_filename, test_filename):

    train_fn_word=train_filename+'.tok'
    test_fn_word = test_filename + '.tok'

    print('done load')
    print('start flatten..')

    print(y_train[:10])
    print(y_test[:10])

    if not isinstance(y_train,list):
        y_train = np.argmax(y_train, 1)
    if not isinstance(y_test, list):
        y_test = np.argmax(y_test, 1)

    print(y_train[:10])
    print(y_test[:10])

    print(X_train[0])
    print(X_train[1])
    print(X_test[0])
    print(X_test[1])

    print('Done load imdb data!')
    n_train = len(y_train)
    n_test = len(y_test)
    print('num train samples: {}'.format(n_train))
    print('num test samples: {}'.format(n_test))
    global_count = {}
    train_docs = []
    test_docs = []

    print('start counting training...')
    for i in range(n_train):
        local_count = {}
        words = X_train[i]
        if len(words.shape)>1:
            words=words.flatten()
        for w in words:
            if w != 0:
                if not local_count.get(w):
                    local_count[w] = 0
                local_count[w] += 1
                if not global_count.get(w):
                    global_count[w] = 0
                global_count[w] += 1
        train_docs.append(local_count)

    print('done counting trainining...')
    print('start write train data')
    with open(train_filename, 'w') as f:
        c = 0
        for doc in train_docs:
            doc = od(sorted(doc.items()))
            y = y_train[c]
            f.write(str(y) + ' ')
            for k, v in doc.items():
                f.write(str(k) + ':' + str(v) + ' ')
            f.write('\n')
            c += 1
    print('done write train data...')
    print('start counting testing...')

    for i in range(n_test):
        local_count = {}
        words = X_test[i]
        if len(words.shape)>1:
            words=words.flatten()
        for w in words:
            # print(w)
            if w != 0:
                if not local_count.get(w):
                    local_count[w] = 0
                local_count[w] += 1
                if not global_count.get(w):
                    global_count[w] = 0
                global_count[w] += 1
        test_docs.append(local_count)

    print('done counting test...')
    print('start write test data...')

    with open(test_filename, 'w') as f:
        c = 0
        for doc in test_docs:
            doc = od(sorted(doc.items()))
            y = y_test[c]
            f.write(str(y) + ' ')
            for k, v in doc.items():
                f.write(str(k) + ':' + str(v) + ' ')
            f.write('\n')
            c += 1
    print('done write test data...')

def gen_imdb_tf_form_batch(train_filename, test_filename, input_dir):
    all_files = glob.glob(input_dir + '/train/*.pkl')
    X_train=None
    y_train=None
    print('start get obj train')
    for i,fn in enumerate(all_files):
        (Xs, ys) = pickle.load(open(fn, 'rb'))
        if X_train is None:
            X_train=Xs
            y_train=ys
        else:
            X_train=np.concatenate((X_train,Xs), axis=0)
            y_train = np.concatenate((y_train, ys), axis=0)
        if i%100==10:
            print('done {}'.format(i))
            # break

    all_files = glob.glob(input_dir + '/test/*.pkl')
    X_test = None
    y_test = None
    print('start get obj test')
    for i, fn in enumerate(all_files):
        (Xs, ys) = pickle.load(open(fn, 'rb'))
        if X_test is None:
            X_test = Xs
            y_test = ys
        else:
            X_test = np.concatenate((X_test,Xs), axis=0)
            y_test = np.concatenate((y_test,ys), axis=0)
        if i%100==10:
            print('done {}'.format(i))
            # break

    gen_imdb_tf_from_obj(X_train, y_train, X_test, y_test, train_filename, test_filename)

if __name__ == '__main__':
    mode='PREPARE_HA'

    if len(sys.argv) >= 2:
        mode=sys.argv[1]

    if mode=='PREPARE_HA':
        maxnbw = 30000
        maxsenl = 200
        maxsen = 50
        inputfn = './data/data.json'
        output_dir = './data/big_imdb_prep_ha50200.pkl'
        if len(sys.argv) >= 3:
            maxnbw = sys.argv[2]
        if len(sys.argv) >= 4:
            maxsenl = sys.argv[3]
        if len(sys.argv) >= 5:
            maxsen = sys.argv[4]
        if len(sys.argv) >= 6:
            inputfn = sys.argv[5]
        if len(sys.argv) >= 7:
            output_dir = sys.argv[6]

        preprocess_rawbig_imdb_ha(MAX_NB_WORDS=int(maxnbw),
                          MAX_SENT_LENGTH=int(maxsenl),
                          MAX_SENTS=int(maxsen),
                          inputfn=inputfn,
                          output_dir=output_dir)
    elif mode=='PREPARE_SDR':
        output_train_fn='./data/big_imdb_count_data/big_imdb_raw_train.data'
        output_test_fn='./data/big_imdb_count_data/big_imdb_raw_test.data'
        input_dir='./data/big_imdb_prep_ha50200/'
        if len(sys.argv) >= 3:
            output_train_fn = sys.argv[4]
        if len(sys.argv) >= 4:
            output_test_fn = sys.argv[5]
        if len(sys.argv) >= 5:
            input_dir = sys.argv[6]

        gen_imdb_tf_form_batch(output_train_fn, output_test_fn, input_dir)

max_stars_repo_name: zizuzi/DeepEmbedding

max_issues_repo_name: zizuzi/DeepEmbedding

max_forks_repo_name: thaihungle/DeepTopicEmbedding


text:
"""Module containing various simulation diagnostic tools."""
import pathlib
from copy import copy
import numpy as np
from matplotlib import pyplot, colors, cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import unyt as u
import colorcet as cc

def phase_space_hist(iteration, tseries, *, nbuz=200, nbz=200, uzmin=40.0, uzmax=None):
    """
    Compute the phase-space uz vs z 2D histogram,
    with given number of bins.
    """
    if uzmax is None:
        hist_range = None
    else:
        hist_range = [None, [uzmin, uzmax]]

    z, uz, w = tseries.get_particle(
        var_list=["z", "uz", "w"],
        species="electrons",
        iteration=iteration,
        select={"uz": [uzmin, uzmax]},
    )
    # convert from m to um
    z_um = z * 1e6

    H, zedges, uzedges = np.histogram2d(
        z_um,
        uz,
        bins=(nbz, nbuz),
        weights=w,  # convert to count of real electrons
        range=hist_range,
    )
    Z, UZ = np.meshgrid(zedges, uzedges)

    H = H.T  # Let each row list bins with common x range.

    return H, Z, UZ

def phase_space_plot(
    iteration,
    tseries,
    uzmax=None,
    vmax=None,
    save_path=pathlib.Path.cwd(),
):
    """
    Plot the longitudinal (uz v z) electron phase space.
    """

    H, Z, UZ = phase_space_hist(iteration, tseries, uzmax=uzmax)

    fig, ax = pyplot.subplots(figsize=(7, 5))

    current_time = (tseries.current_t * u.second).to(u.picosecond)
    ax.set_title(f"t = {current_time:.2f} (iteration {iteration:,g})")

    img = ax.pcolormesh(
        Z,
        UZ,
        H,
        cmap=cm.get_cmap("plasma"),
        norm=colors.SymLogNorm(linthresh=1, vmin=0.0, vmax=vmax, base=10),
    )

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="4%", pad=0.02)
    cbar = ax.figure.colorbar(
        img,
        cax=cax,
    )
    cbar.set_label(r"number of electrons")

    ax.set_xlabel(r"$z$ ($\mathrm{\mu m}$)")
    ax.set_ylabel(r"$u_z$ ($m_e c$)")

    filename = pathlib.Path(save_path) / f"phasespace{iteration:06d}.png"
    fig.savefig(filename)
    pyplot.close(fig)

def laser_density_plot(
    iteration,
    tseries,
    rho_field_name="rho_electrons",
    laser_polarization="x",
    save_path=pathlib.Path.cwd(),
    n_c=1.7419595910637713e27,  # 1/m^3
    E0=4013376052599.5396,  # V/m
    ylim=None,
):
    """
    Plot on the same figure the laser pulse envelope and the electron density.
    """

    laser_cmap = copy(cc.m_kr)
    laser_cmap.set_under("black", alpha=0)

    rho, rho_info = tseries.get_field(
        field=rho_field_name,
        iteration=iteration,
    )
    envelope, env_info = tseries.get_laser_envelope(
        iteration=iteration, pol=laser_polarization
    )
    # get longitudinal field
    e_z_of_z, e_z_of_z_info = tseries.get_field(
        field="E",
        coord="z",
        iteration=iteration,
        slice_across="r",
    )
    # the field "rho" has (SI) units of charge/volume (Q/V), C/(m^3)
    # the initial density n_e has units of N/V, N = electron number
    # multiply by electron charge q_e to get (N e) / V
    # so we get Q / N e, which is C/C, i.e. dimensionless
    # Note: one can also normalize by the critical density n_c

    fig, ax = pyplot.subplots(figsize=(10, 6))

    R = rho / (u.electron_charge.to_value("C") * n_c)
    E = envelope / E0

    im_rho = ax.imshow(
        R,
        extent=rho_info.imshow_extent * 1e6,  # conversion to microns
        origin="lower",
        norm=colors.SymLogNorm(
            vmin=-0.3, vmax=6.0, linthresh=1.0e-3, linscale=0.15, base=10
        ),
        cmap=cm.get_cmap("cividis"),
    )
    im_envelope = ax.imshow(
        E,
        extent=env_info.imshow_extent * 1e6,
        origin="lower",
        norm=colors.Normalize(vmin=0.0, vmax=5.0),
        cmap=laser_cmap,
    )
    im_envelope.set_clim(vmin=1.0)

    # plot longitudinal field
    ax.plot(e_z_of_z_info.z * 1e6, e_z_of_z / E0 * 15 - 15, color="0.75")
    ax.axhline(-15, color="0.65", ls="-.")

    cbaxes_rho = inset_axes(
        ax,
        width="3%",  # width = 10% of parent_bbox width
        height="46%",  # height : 50%
        loc=2,
        bbox_to_anchor=(1.01, 0.0, 1, 1),
        bbox_transform=ax.transAxes,
        borderpad=0,
    )
    cbaxes_env = inset_axes(
        ax,
        width="3%",  # width = 5% of parent_bbox width
        height="46%",  # height : 50%
        loc=3,
        bbox_to_anchor=(1.01, 0.0, 1, 1),
        bbox_transform=ax.transAxes,
        borderpad=0,
    )
    cbar_env = fig.colorbar(
        mappable=im_envelope,
        orientation="vertical",
        ticklocation="right",
        cax=cbaxes_env,
    )
    cbar_rho = fig.colorbar(
        mappable=im_rho, orientation="vertical", ticklocation="right", cax=cbaxes_rho
    )
    cbar_env.set_label(r"$eE_{x} / m c \omega_\mathrm{L}$")
    cbar_rho.set_label(r"$n_{e} / n_\mathrm{cr}$")

    ax.set_ylabel(r"${} \;(\mu m)$".format(rho_info.axes[0]))
    ax.set_xlabel(r"${} \;(\mu m)$".format(rho_info.axes[1]))

    if ylim is not None:
        ax.set_ylim(*ylim)

    current_time = (tseries.current_t * u.second).to(u.picosecond)
    ax.set_title(f"t = {current_time:.2f} (iteration {iteration:,g})")

    filename = save_path / f"rho{iteration:06d}.png"

    fig.subplots_adjust(right=0.85)
    fig.savefig(filename)
    pyplot.close(fig)

def particle_energy_histogram(
    tseries,
    iteration: int,
    energy_min=1,
    energy_max=500,
    delta_energy=1,
    cutoff=35,  # CHANGEME
    species=None,
):
    """
    Compute the weighted particle energy histogram from ``tseries`` at step ``iteration``.

    :param tseries: whole simulation time series
    :param iteration: time step in the simulation
    :param energy_min: lower energy threshold (MeV)
    :param energy_max: upper energy threshold (MeV)
    :param delta_energy: size of each energy bin (MeV)
    :param cutoff: upper threshold for the histogram, in pC / MeV
    :return: histogram values and bin edges
    """
    nbins = (energy_max - energy_min) // delta_energy
    energy_bins = np.linspace(start=energy_min, stop=energy_max, num=nbins + 1)

    ux, uy, uz, w = tseries.get_particle(
        ["ux", "uy", "uz", "w"], iteration=iteration, species=species
    )

    if (np.ndim(w) == 1) and (w.size == 1):
        fill_value = w.item(0)
        w = np.full_like(ux, fill_value)

    energy = (u.electron_mass * u.clight ** 2).to_value("MeV") * np.sqrt(
        1 + ux ** 2 + uy ** 2 + uz ** 2
    )

    # Explanation of weights:
    #     1. convert electron charge from C to pC (factor 1e12)
    #     2. multiply by weight w to get real number of electrons
    #     3. divide by energy bin size delta_energy to get charge / MeV
    hist, _ = np.histogram(
        energy,
        bins=energy_bins,
        weights=u.elementary_charge.to_value("pC") / delta_energy * w,
    )

    # cut off histogram
    np.clip(hist, a_min=None, a_max=cutoff, out=hist)

    return hist, energy_bins, nbins

def main():
    import random
    from openpmd_viewer.addons import LpaDiagnostics
    import signac

    random.seed(42)

    proj = signac.get_project(search=False)
    job = next(iter(proj))

    h5_path = pathlib.Path(job.ws) / "diags" / "hdf5"
    time_series = LpaDiagnostics(h5_path)

    it = random.choice(time_series.iterations.tolist())
    print(f"job {job.id}, iteration {it}")

    # compute 1D histogram
    energy_hist, bin_edges, _ = particle_energy_histogram(
        tseries=time_series,
        iteration=it,
        species="electrons",
        cutoff=np.inf,  # no cutoff
    )

    laser_density_plot(iteration=it, tseries=time_series)
    phase_space_plot(iteration=it, tseries=time_series)

if __name__ == "__main__":
    main()

max_stars_repo_name: berceanu/signac-driven-fbpic

max_issues_repo_name: berceanu/signac-driven-fbpic

max_forks_repo_name: berceanu/signac-driven-fbpic


text:
import numpy as np

def bounding_box_fit_in_image_boundary_(bbox_: np.ndarray, image_size):
    min_x = 0
    min_y = 0
    max_x = image_size[0] - 1
    max_y = image_size[1] - 1
    bbox_[0] = np.minimum(np.maximum(bbox_[0], min_x), max_x)
    bbox_[1] = np.minimum(np.maximum(bbox_[1], min_y), max_y)
    bbox_[2] = np.minimum(np.maximum(bbox_[2], min_x), max_x)
    bbox_[3] = np.minimum(np.maximum(bbox_[3], min_y), max_y)


text:
'''
Created on Feb 23, 2016

@author: jason
'''
import numpy as np
from sklearntools.sklearntools import StagedEstimator, MaskedEstimator,\
    ColumnSubsetTransformer, NonNullSubsetFitter, safe_assign_column
from sklearn.linear_model.base import LinearRegression
from sklearn.linear_model.logistic import LogisticRegression
from sklearntools.calibration import CalibratedEstimatorCV, ResponseTransformingEstimator,\
    LogTransformer, PredictorTransformer, HazardToRiskEstimator,\
    MovingAverageSmoothingEstimator, ThresholdClassifier, ProbaPredictingEstimator

from sklearntools.feature_selection import SingleEliminationFeatureImportanceEstimatorCV,\
    BackwardEliminationEstimator, UnivariateFeatureImportanceEstimatorCV,\
    BestKFeatureSelector
from numpy.testing.utils import assert_raises
from sklearntools.glm import GLM
import statsmodels.api as sm
import warnings
import pandas
from sklearntools.model_selection import ModelSelector
from sklearntools.scoring import log_loss_metric
from sklearn.ensemble.forest import RandomForestRegressor
from numpy.ma.testutils import assert_array_almost_equal
from sklearntools.earth import Earth
from sklearntools.kfold import CrossValidatingEstimator
from sklearn.metrics.regression import r2_score
from sklearn.model_selection import KFold
from nose.tools import assert_list_equal
warnings.simplefilter("error")

def test_safe_assign_column():
    data = pandas.DataFrame({'A': [1,2,3], 'B': [4,5,6]})
    safe_assign_column(data, 'A', [7,8,9])
    assert_list_equal(list(sorted(data.columns)), ['A', 'B'])

def test_single_elimination_feature_importance_estimator_cv():
    np.random.seed(0)
    m = 100000
    n = 6
    factor = .9

    X = np.random.normal(size=(m,n))
    beta = 100 * np.ones(shape=n)
    for i in range(1, n):
        beta[i] = factor * beta[i-1]
    beta = np.random.permutation(beta)[:,None]

    y = np.dot(X, beta) + 0.01 * np.random.normal(size=(m, 1))

    target_sequence = np.ravel(np.argsort(beta ** 2, axis=0))
    model1 = SingleEliminationFeatureImportanceEstimatorCV(LinearRegression())
    model1.fit(X, y)
    fitted_sequence = np.ravel(np.argsort(model1.feature_importances_, axis=0))

    np.testing.assert_array_equal(fitted_sequence, target_sequence)

def test_univariate_feature_importance_estimator_cv():
    np.random.seed(0)
    m = 100000
    n = 6
    factor = .9

    X = np.random.normal(size=(m,n))
    beta = 100 * np.ones(shape=n)
    for i in range(1, n):
        beta[i] = factor * beta[i-1]
    beta = np.random.permutation(beta)[:,None]

    y = np.dot(X, beta) + 0.01 * np.random.normal(size=(m, 1))

    target_sequence = np.ravel(np.argsort(beta ** 2, axis=0))
    model1 = UnivariateFeatureImportanceEstimatorCV(LinearRegression())
    model1.fit(X, y)
    fitted_sequence = np.ravel(np.argsort(model1.feature_importances_, axis=0))

    np.testing.assert_array_equal(fitted_sequence, target_sequence)

def test_k_best_feature_selector():
    np.random.seed(0)
    m = 100000
    n = 6
    factor = .9

    X = np.random.normal(size=(m,n))
    beta = 100 * np.ones(shape=n)
    for i in range(1, n):
        beta[i] = factor * beta[i-1]
    beta = np.random.permutation(beta)[:,None]
#     beta = np.random.normal(size=(n,1))

    y = np.dot(X, beta) + 0.01 * np.random.normal(size=(m, 1))

    target_vars = np.ravel(np.argsort(beta ** 2, axis=0))[::-1][:3]
    target_support = np.zeros(shape=n, dtype=bool)
    target_support[target_vars] = True

    model1 = BestKFeatureSelector(UnivariateFeatureImportanceEstimatorCV(LinearRegression()), k=3)
    model1.fit(X, y)

    np.testing.assert_array_equal(model1.support_, target_support)

def test_backward_elimination_estimation():
    np.random.seed(0)
    m = 100000
    n = 6
    factor = .9

    X = np.random.normal(size=(m,n))
    beta = 100 * np.ones(shape=n)
    for i in range(1, n):
        beta[i] = factor * beta[i-1]
    beta = np.random.permutation(beta)[:,None]
#     beta = np.random.normal(size=(n,1))

    y = np.dot(X, beta) + 0.01 * np.random.normal(size=(m, 1))

    target_sequence = np.ravel(np.argsort(beta ** 2, axis=0))
    model1 = BackwardEliminationEstimator(SingleEliminationFeatureImportanceEstimatorCV(LinearRegression()))
    model1.fit(X, y)

#     model2 = BRFE(FeatureImportanceEstimatorCV(LinearRegression()))
#     model2.fit(X, y)

    np.testing.assert_array_equal(model1.elimination_sequence_, target_sequence)

def test_multiple_response_regressor():
    np.random.seed(1)
    m = 100000
    n = 10

    X = np.random.normal(size=(m,n))
    beta1 = np.random.normal(size=(n,1))
    beta2 = np.random.normal(size=(n,1))

    y1 = np.dot(X, beta1)
    p2 = 1. / (1. + np.exp( - np.dot(X, beta2)))
    y2 = np.random.binomial(n=1, p=p2)
    y = np.concatenate([y1, y2], axis=1)

    model = MaskedEstimator(LinearRegression(), [True, False]) & MaskedEstimator(ProbaPredictingEstimator(LogisticRegression()), [False, True])
#     MultipleResponseEstimator([('linear', np.array([True, False], dtype=bool), LinearRegression()),
#                                        ('logistic', np.array([False, True], dtype=bool), ProbaPredictingEstimator(LogisticRegression()))])
    model.fit(X, y)

    assert np.mean(beta1 - model.estimators_[0].estimator_.coef_) < .01
    assert np.mean(beta2 - model.estimators_[1].estimator_.estimator_.coef_) < .01
    model.get_params()
    model.predict(X)

def test_calibration():
    np.random.seed(1)
    m = 10000
    n = 10

    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    y_lin = np.dot(X, beta)
    y_clas = np.random.binomial( 1, 1. / (1. + np.exp(-y_lin)) )
    y = np.concatenate([y_lin, y_clas], axis=1)
    estimator = MaskedEstimator(LinearRegression(), np.array([True, False], dtype=bool))
    calibrator = MaskedEstimator(LogisticRegression(), [False, True])
#     estimator = linear_regressor & calibrator
#     MultipleResponseEstimator([('estimator', np.array([True, False], dtype=bool), LinearRegression())])
#     calibrator = MultipleResponseEstimator([('calibrator', np.array([False, True], dtype=bool), LogisticRegression())])
    model = CalibratedEstimatorCV(estimator, calibrator, cv=KFold(n_splits=4, shuffle=True), n_jobs=1)
    model.fit(X, y)
    assert np.max(beta[:, 0] - model.estimator_.estimator_.coef_) < .000001
    assert np.max(model.calibrator_.estimator_.coef_ - 1.) < .1

def test_predictor_transformer_calibration():
    np.random.seed(1)
    m = 10000
    n = 10

    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    y_lin = np.dot(X, beta)
    y_clas = np.random.binomial( 1, 1. / (1. + np.exp(-y_lin)) )
    y = np.concatenate([y_lin, y_clas], axis=1)
    estimator = MaskedEstimator(LinearRegression(), np.array([True, False], dtype=bool))
    calibrator = MaskedEstimator(LogisticRegression(), [False, True])
#     estimator = linear_regressor & calibrator
#     MultipleResponseEstimator([('estimator', np.array([True, False], dtype=bool), LinearRegression())])
#     calibrator = MultipleResponseEstimator([('calibrator', np.array([False, True], dtype=bool), LogisticRegression())])
    model = PredictorTransformer(estimator) >> calibrator
    model.fit(X, y)
    assert np.max(beta[:, 0] - model.intermediate_stages_[0].estimator_.estimator_.coef_) < .000001
    assert np.max(model.final_stage_.estimator_.coef_ - 1.) < .1

def test_pipeline():
    np.random.seed(1)
    m = 10000
    n = 10

    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    beta[np.random.binomial(p=2.0/float(n), n=1, size=n).astype(bool)] = 0
    y = np.dot(X, beta) + 0.5 * np.random.normal(size=(m, 1))
    beta_reduced = beta[beta != 0]

    model = BackwardEliminationEstimator(SingleEliminationFeatureImportanceEstimatorCV(LinearRegression()))
    model >>= LinearRegression()

    model.fit(X, y)
    assert np.max(np.abs(model.final_stage_.coef_ - beta_reduced)) < .1

def test_response_transforming_estimator():
    np.random.seed(1)
    m = 10000
    n = 10

    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    sigma = .1
    y_pre = np.dot(X, beta) + sigma * np.random.normal(size=(m,1))
    y_post = np.exp(y_pre)

    model = ResponseTransformingEstimator(LinearRegression(), LogTransformer(offset=0.))
    model.fit(X, y_post)

    assert np.abs(np.mean(model.predict(X) - y_pre)) < .01

    # Because LinearRegression has no transform method
    assert_raises(AttributeError, lambda: model.transform(X))

def test_hazard_to_risk():
    np.random.seed(1)

    m = 10000
    n = 10

    # Simulate an event under constant hazard, with hazard = X * beta and
    # iid exponentially distributed exposure times.
    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    hazard = np.exp(np.dot(X, beta))
    exposure = np.random.exponential(size=(m,1))
    rate = np.random.poisson(hazard * exposure) / exposure

    model = CalibratedEstimatorCV(GLM(sm.families.Gaussian(sm.families.links.log), add_constant=False),
                                  ProbaPredictingEstimator(ThresholdClassifier(HazardToRiskEstimator(LogisticRegression()))))

    model.fit(X, rate, exposure=exposure)

    y_pred = model.predict(X, exposure)
    assert np.abs((np.sum(y_pred) - np.sum(rate > 0)) / np.sum(rate > 0))  < .1
    assert np.max(np.abs(model.estimator_.coef_ - beta[:,0])) < .1

def test_hazard_to_risk_staged():
    np.random.seed(1)

    m = 10000
    n = 10

    # Simulate an event under constant hazard, with hazard = X * beta and
    # iid exponentially distributed exposure times.
    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    hazard = np.exp(np.dot(X, beta))
    exposure = np.random.exponential(size=(m,1))
    rate = np.random.poisson(hazard * exposure) / exposure

    model = CalibratedEstimatorCV(GLM(sm.families.Gaussian(sm.families.links.log), add_constant=False),
                                  ProbaPredictingEstimator(ThresholdClassifier(HazardToRiskEstimator(LogisticRegression()))))

    model.fit(X, rate, exposure=exposure)

    y_pred = model.predict(X, exposure)
    assert np.abs((np.sum(y_pred) - np.sum(rate > 0)) / np.sum(rate > 0))  < .1
    assert np.max(np.abs(model.estimator_.coef_ - beta[:,0])) < .1

def test_moving_average_smoothing_estimator():
    np.random.seed(1)

    m = 10000
    n = 10

    # Simulate an event under constant hazard, with hazard = X * beta and
    # iid exponentially distributed exposure times.
    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    hazard = np.exp(np.dot(X, beta))
    exposure = np.random.exponential(size=(m,1))
    rate = np.random.poisson(hazard * exposure) / exposure

    model = CalibratedEstimatorCV(GLM(sm.families.Gaussian(sm.families.links.log), add_constant=False),
                                  ThresholdClassifier(HazardToRiskEstimator(MovingAverageSmoothingEstimator(RandomForestRegressor()))))

    model.fit(X, rate, exposure=exposure)

    y_pred = model.predict(X, exposure)
    assert np.abs((np.sum(y_pred) - np.sum(rate > 0)) / np.sum(rate > 0))  < .1
    assert np.max(np.abs(model.estimator_.coef_ - beta[:,0])) < .1

def test_staged_estimator():
    np.random.seed(1)
    m = 10000
    n = 10

    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    beta[np.random.binomial(p=2.0/float(n), n=1, size=n).astype(bool)] = 0
    y = np.dot(X, beta) + 0.5 * np.random.normal(size=(m, 1))
    beta_reduced = beta[beta != 0]

    stage0 = BackwardEliminationEstimator(SingleEliminationFeatureImportanceEstimatorCV(LinearRegression()))
    stage1 = LinearRegression()
    model = StagedEstimator([stage0, stage1])

    model.fit(X, y)
    assert np.max(np.abs(model.final_stage_.coef_ - beta_reduced)) < .1
#
#     y_lin = np.dot(X, beta)
#     y_clas = np.random.binomial( 1, 1. / (1. + np.exp(-y_lin)) )
#     y = np.concatenate([y_lin, y_clas], axis=1)
#     estimator = mask_estimator(LinearRegression(), np.array([True, False], dtype=bool))
#     calibrator = mask_estimator(LogisticRegression(), [False, True])
# #     estimator = linear_regressor & calibrator
# #     MultipleResponseEstimator([('estimator', np.array([True, False], dtype=bool), LinearRegression())])
# #     calibrator = MultipleResponseEstimator([('calibrator', np.array([False, True], dtype=bool), LogisticRegression())])
#     model = CalibratedEstimatorCV(estimator, calibrator)
#     model.fit(X, y)
#     assert np.max(beta[:, 0] - model.estimator_.estimators_[0][2].coef_) < .000001
#     assert np.max(model.calibrator_.estimators_[0][2].coef_ - 1.) < .1

def test_column_subset_transformer():
    m = 1000
    n = 10
    X = np.random.normal(size=(m,n))
    x_cols = [0,3,4,5]
    y_cols = 9
    sample_weight_cols = 8
    exposure_cols = 7

    subsetter1 = ColumnSubsetTransformer(x_cols=x_cols, y_cols=y_cols,
                                         sample_weight_cols=sample_weight_cols,
                                         exposure_cols=exposure_cols)
    np.testing.assert_array_equal(subsetter1.transform(X), X[:, x_cols])
    args = {'X': X}
    subsetter1.update(args)
    np.testing.assert_array_equal(args['X'], X[:, x_cols])
    np.testing.assert_array_equal(args['y'], X[:, y_cols])
    np.testing.assert_array_equal(args['sample_weight'], X[:, sample_weight_cols])
    np.testing.assert_array_equal(args['exposure'], X[:, exposure_cols])

    X_ = pandas.DataFrame(X, columns=['x%d' % n for n in range(10)])
    x_cols_ = ['x%d' % n for n in x_cols]
    y_cols_ = 'x%d' % y_cols
    sample_weight_cols_ = 'x%d' % sample_weight_cols
    exposure_cols_ = 'x%d' % exposure_cols
    subsetter2 = ColumnSubsetTransformer(x_cols=x_cols_, y_cols=y_cols_,
                                         sample_weight_cols=sample_weight_cols_,
                                         exposure_cols=exposure_cols_)
    np.testing.assert_array_equal(subsetter2.transform(X_), X[:, x_cols])
    args_ = {'X': X_}
    subsetter2.update(args_)
    np.testing.assert_array_equal(args_['X'], X[:, x_cols])
    np.testing.assert_array_equal(args_['y'], X[:, y_cols])
    np.testing.assert_array_equal(args_['sample_weight'], X[:, sample_weight_cols])
    np.testing.assert_array_equal(args_['exposure'], X[:, exposure_cols])

    lin = ColumnSubsetTransformer(x_cols=x_cols_, y_cols=y_cols_) >> LinearRegression()
    lin.fit(X_)
    lin.predict(X_.loc[:, x_cols_])
    lin.score(X_)

def test_model_selector():
    np.random.seed(1)

    m = 10000
    n = 10

    # Simulate an event under constant hazard, with hazard = X * beta and
    # iid exponentially distributed exposure times.
    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    hazard = np.exp(np.dot(X, beta))
    exposure = np.random.exponential(size=(m,1))
    rate = np.random.poisson(hazard * exposure) / exposure
    best_subset = np.ravel(np.argsort(np.abs(beta))[::-1][:3])
    worst_subset = np.ravel(np.argsort(np.abs(beta))[:3])

    basic_model = CalibratedEstimatorCV(GLM(sm.families.Gaussian(sm.families.links.log), add_constant=False),
                                  ProbaPredictingEstimator(ThresholdClassifier(HazardToRiskEstimator(LogisticRegression()))))

    model1 = CrossValidatingEstimator(ColumnSubsetTransformer(x_cols=best_subset) >> basic_model, metric=log_loss_metric)
    model2 = CrossValidatingEstimator(ColumnSubsetTransformer(x_cols=worst_subset) >> basic_model, metric=log_loss_metric)

    model = ModelSelector([model1, model2])

    model.fit(X, rate, exposure=exposure)
    np.testing.assert_array_equal(model.best_estimator_.estimator_.intermediate_stages_[0].x_cols, best_subset)

def test_cross_validating_estimator():
    np.random.seed(1)

    m = 1000
    n = 5

    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    y = np.ravel(np.dot(X, beta)) + np.random.normal(.5, size=m)

    model = CrossValidatingEstimator(LinearRegression(), cv=KFold(n_splits=4, shuffle=True), n_jobs=2)

    y_pred_cv = model.fit_predict(X, y)

    y_pred = model.predict(X)

    assert r2_score(np.ravel(y_pred_cv), np.ravel(y_pred)) > .98

def test_non_null_row_subset_fitter():
    np.random.seed(1)

    m = 10000
    n = 10

    # Simulate an event under constant hazard, with hazard = X * beta and
    # iid exponentially distributed exposure times.
    X = np.random.normal(size=(m,n))
    beta = np.random.normal(size=(n,1))
    y = np.ravel(np.dot(X, beta))

    missing = np.random.binomial(p=.001, n=1, size=(m,n)) == 1
    X[missing] = None

    model = NonNullSubsetFitter(LinearRegression())
    model.fit(X, y)
    assert np.max(np.abs(np.ravel(beta) - model.estimator_.coef_)) < .001

def test_linear_transformation():
    np.random.seed(1)

    m = 10000
    n = 10

    X = np.random.normal(size=(m,n))
    beta1 = np.random.normal(size=(n,1))
    y1 = np.ravel(np.dot(X, beta1))

    beta2 = np.random.normal(size=(n,1))
    y2 = np.ravel(np.dot(X, beta2))

    model1 = (Earth() >> LinearRegression()).fit(X, y1)
    model2 = Earth().fit(X, y2)

    combination = 2*model1 - model2

    assert_array_almost_equal(combination.predict(X), 2 * np.ravel(model1.predict(X)) - np.ravel(model2.predict(X)))

if __name__ == '__main__':
    import sys
    import nose
    # This code will run the test in this file.'
    module_name = sys.modules[__name__].__file__

    result = nose.run(argv=[sys.argv[0],
                            module_name,
                            '-s', '-v'])

max_stars_repo_name: jcrudy/sklearntools

max_issues_repo_name: jcrudy/sklearntools

max_forks_repo_name: jcrudy/sklearntools

reason: import numpy,from numpy,import statsmodels