text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Sun May 12 19:11:09 2019

@author: Soriba
"""
import os
import numpy as np
from keras.models import load_model
from PIL import ImageTk, Image, ImageDraw, ImageChops
from process import get_px, get_burst_start

#if you use python2 or python3
try:
    import tkinter as tk
    import tkinter.filedialog as tkfd
except ImportError:
    import Tkinter as tk
    import tkFileDialog as tkfd

if __name__ == '__main__':
    #We load the architecture and the weights of the pretrained model
    model=load_model('model.h5')
    """
    Function to draw a rectangle from the coordinates
    """
    def drawrect(drawcontext, xy, outline=None, width=0):
        (x1, y1), (x2, y2) = xy
        points = (x1, y1), (x2, y1), (x2, y2), (x1, y2), (x1, y1)
        drawcontext.line(points, fill=outline, width=width)

    def forget():
            CB1.destroy()
            CB2.destroy()

    """
    Function to detect the white border of the image
    """

    def trim(im):
        bg = Image.new(im.mode, im.size, im.getpixel((0,0)))
        diff = ImageChops.difference(im, bg)
        diff = ImageChops.add(diff, diff, 2.0, -100)
        bbox = diff.getbbox()
        if bbox:
            # in the form (x1, y1, x2, y2)
            return bbox

    """
    Function to load an image from the file system and make a prediction
    """

    def select_image():
        global panelA, Label
        #the user specify the path to the spectrogram he wants to predict
        path=tkfd.askopenfilename(initialdir = "/mnt/disque2/IMT/base",
                                  title = "Select a spectrogram",
                        filetypes = (("png files","*.png"),
                                     ("all files","*.*")))

        filename=os.path.basename(path)
        #file containing the starts of the bursts
        filename_txt=filename.split('.')[0]+'.wav.txt'
        classes = np.load('classes.npy')
        #class of the image (name of the folder in which it is)
        groundtruth = os.path.basename(os.path.dirname(path))

        if path:
            #we adapt the size of the image for the neural network
            img=Image.open(path).resize((709,532))
            #we convert the image into a matrix
            pix=np.array(img)
            #we remove the white border
            box = trim(img)
            max_x=box[2]-box[0]+1
            max_y=box[3]-box[1]+1
            print('max_x : {}'.format(max_x))
            bandwidths_px, burst_durations_px = get_px(max_x,max_y)
            #starts of the bursts
            burst_starts = get_burst_start(groundtruth, filename_txt, max_y,
                                           box[1])
            burst_len = burst_durations_px[np.where(classes==groundtruth)[0][0]]
            #f0 = 0, it corresponds to the center of the image (on the x axis)
            f0_px=box[0] + np.ceil(max_x/2)
            bw=bandwidths_px[np.where(classes==groundtruth)[0][0]]
            #normalization
            processed = pix/255
            #prediction of the model
            prediction = model.predict_classes(processed[np.newaxis,:])[0]

        if var1.get()==1:
            for k in range(len(burst_starts)):
                draw = ImageDraw.Draw(img)
                drawrect(
                        draw,
                         [(f0_px-bw/2, burst_starts[k]),
                          (f0_px+bw/2, burst_starts[k]+burst_len)],
                          outline="blue",
                          width=2
                          )

        display = ImageTk.PhotoImage(img)
        #initialization of panelA
        if panelA is None:
            #panelA will display the image
            panelA = tk.Label(image=display)
            panelA.image = display
            panelA.pack(side='bottom', padx=10, pady=10)

        else:
            #if panelA is already initialized
            panelA.configure(image=display)
            panelA.image=display
            forget()

        #prediction on a green background if it is correct, red otherwise
        bg='green' if classes[prediction] == groundtruth  else 'red'

        text = "Predicted class : {} \nTrue class : {}".format(classes[prediction],
                                  groundtruth)

        if Label is None:
            Label= tk.Label(top_frame, text=text,
                       fg='white', bg=bg)
            Label.pack(fill='x')

        else:
            Label.configure(text=text, bg=bg)

    #main window
    root=tk.Tk()
    root.title('application')
    top_frame = tk.Frame(root).pack()
    bottom_frame = tk.Frame(root).pack(side = "bottom")

    panelA=None
    Label=None
    var1=tk.IntVar()
    var2=tk.IntVar()

    CB1=tk.Checkbutton(bottom_frame, text = "groundtruth detection",
                       variable=var1, command=forget)
    CB1.pack(side = 'right')
    CB2=tk.Checkbutton(bottom_frame, text = "predicted detection",
                       variable=var2, command=forget)
    CB2.pack(side= 'right')
    btn = tk.Button(bottom_frame, text = 'Choose a spectrogram',
                    command=select_image)
    btn.pack(side='bottom', fill='both', expand='yes', padx='10', pady='10')

    root.mainloop()

max_issues_repo_name: soribadiaby/CNN-signaldetection

max_forks_repo_name: soribadiaby/CNN-signaldetection


text:
import numpy as np

class Conv1DT:

    def __call__(self, x, weight, bias, stride, padding):
        x = x.T
        in_len, in_chan = x.shape
        in_chan, out_chan, kernel_size = weight.shape
        w = weight.reshape(in_chan, -1)

        temp = np.dot(x, w).reshape(in_len, out_chan,
                                    kernel_size).transpose(0, 2, 1)
        temp = temp.reshape(in_len*kernel_size, out_chan)
        temp = np.pad(temp, ((kernel_size//2, kernel_size//2), (0, 0)))
        temp = temp.reshape(in_len+1, 2, kernel_size//2, out_chan)
        output = temp.sum(1).reshape(-1, out_chan) + bias[None]

        return output[padding:-padding].T

ext: py

max_stars_repo_name: neblar/numpynn

max_issues_repo_name: neblar/numpynn

max_forks_repo_name: neblar/numpynn

reason: import numpy


text:
"""
This module includes helper functions for plotting.
"""
import matplotlib.pyplot as plt
import numpy as np

def plot_image(image, figsize=None, title=None, xlim=None, ylim=None, extent=None, colorbar=False, cmap=None):
    if cmap is None:
        cmap = 'nipy_spectral'
    if figsize is None:
        figsize = (9, 6)
    fig = plt.figure(figsize=figsize)
    if extent is None:
        plt.imshow(image, origin='lower', cmap=cmap)  # 'hot' 'nipy_spectral'
    else:
        ax = fig.add_subplot(1, 1, 1)
        ax.imshow(image, origin='lower', cmap=cmap, extent=extent)
    if title:
        plt.title(title)
    if xlim is not None:
        plt.xlim(xlim)
    if ylim is not None:
        plt.ylim(ylim)
    if colorbar:
        plt.colorbar(orientation="horizontal")
    plt.show()

def plot_images(images, figsize=None, title_list=None, xlim=None, ylim=None, extend=None, colorbar=False, cmap=None):
    if cmap is None:
        cmap = 'nipy_spectral'
    n_images = len(images)
    if figsize is None:
        figwidth = np.min((n_images * 6, 24))
        figsize = (figwidth, 6)
    fig, axes = plt.subplots(1, n_images, figsize=figsize)
    for count, i_image in enumerate(images):
        if extend is None:
            im = axes[count].imshow(i_image, origin='lower', cmap=cmap)
        else:
            im = axes[count].imshow(i_image, origin='lower', cmap=cmap, extend=extend)
        if title_list is not None:
            axes[count].set_title(title_list[count])
        if xlim is not None:
            axes[count].set_xlim(xlim)
        if ylim is not None:
            axes[count].set_ylim(ylim)
        if colorbar:
            fig.colorbar(im, ax=axes[count])
    plt.show()

def plot_graph(y, figsize=None, title=None, xlim=None, ylim=None):
    if figsize is None:
        figsize = (9, 6)
    plt.figure(figsize=figsize)
    plt.plot(y)
    if title:
        plt.title(title)
    if xlim is not None:
        plt.xlim(xlim)
    else:
        plt.xlim([0, len(y)])
    if ylim is not None:
        plt.ylim(ylim)
    plt.show()

def plot_graphs(graph_list, title_list=[], xlim=None, ylim=None):
    n_graphs = len(graph_list)
    fig_depth = 2 * n_graphs
    # fig_depth = np.max((4, fig_depth))
    fig_depth = np.min((24, fig_depth))
    plt.figure(figsize=(12, fig_depth))
    for i in range(n_graphs):
        plt.subplot(n_graphs, 1, i + 1)
        plt.plot(graph_list[i])
        if title_list:
            plt.title(title_list[i])
        if xlim is not None:
            plt.xlim(xlim)
        if ylim is not None:
            plt.ylim(ylim)
    plt.tight_layout()
    plt.show()


text:
"""
Simple Horten Wing as used by Richards. Baseline and simplified models

"""
import numpy as np
from cases.hangar.horten_wing import HortenWing
import sharpy.utils.algebra as algebra

class Baseline(HortenWing):

    def set_properties(self):

        # Wing geometry
        self.span = 20.0  # [m]
        self.sweep_LE = 20 * np.pi / 180  # [rad] Leading Edge Sweep
        self.c_root = 1.0  # [m] Root chord - Richards
        self.taper_ratio = 0.25  # Richards
        self.thrust_nodes = [self.n_node_fuselage - 1,
                             self.n_node_fuselage + self.n_node_wing + 1]

        self.loc_cg = 0.45 # CG position wrt to LE (from sectional analysis)
        # EA is the reference in NATASHA - defined with respect to the midchord. SHARPy is wrt to LE and as a pct of
        # local chord
        self.main_ea_root = 0.33
        self.main_ea_tip = 0.33
        self.n_mass = 2 * self.n_elem_wing

        # FUSELAGE GEOMETRY
        self.fuselage_width = 1.65/2
        self.c_fuselage = self.c_root

        # WASH OUT
        self.washout_root = 0*np.pi/180
        self.washout_tip = -2 * np.pi / 180

        # Horseshoe wake
        self.horseshoe = False
        self.wake_type = 2
        self.dt_factor = 1
        self.dt = 1 / self.M / self.u_inf * self.dt_factor

        # Dynamics
        self.n_tstep = int(self.physical_time/self.dt)
        self.gust_intensity = 0.1

        # Numerics
        self.tolerance = 1e-12
        self.fsi_tolerance = 1e-10
        self.relaxation_factor = 0.2

    def update_mass_stiffness(self, sigma=1., sigma_mass=1.):
        """
        Set's the mass and stiffness properties of the default wing

        Returns:

        """

        n_elem_fuselage = self.n_elem_fuselage
        n_elem_wing = self.n_elem_wing
        n_node_wing = self.n_node_wing
        n_node_fuselage = self.n_node_fuselage
        c_root = self.c_root
        taper_ratio = self.taper_ratio

        # Local chord to root chord initialisation
        c_bar_temp = np.linspace(c_root, taper_ratio * c_root, n_elem_wing)

        # Structural properties at the wing root section from Richards 2016
        ea = 1e6
        ga = 1e6
        gj = 4.24e5
        eiy = 3.84e5
        eiz = 2.46e7

        root_i_beam = IBeam()
        root_i_beam.build(c_root)
        root_i_beam.rotation_axes = np.array([0, self.main_ea_root-0.25, 0])

        root_airfoil = Airfoil()
        root_airfoil.build(c_root)
        root_i_beam.rotation_axes = np.array([0, self.main_ea_root-0.25, 0])

        mu_0 = root_i_beam.mass + root_airfoil.mass
        j_xx = root_i_beam.ixx + root_airfoil.ixx
        j_yy = root_i_beam.iyy + root_airfoil.iyy
        j_zz = root_i_beam.izz + root_airfoil.izz

        # Number of stiffnesses used
        n_stiffness = self.n_stiffness

        # Initialise the stiffness database
        base_stiffness = self.base_stiffness

        stiffness_root = sigma * np.diag([ea, ga, ga, gj, eiy, eiz])
        stiffness_tip = taper_ratio ** 2 * stiffness_root

        # Assume a linear variation in the stiffness. Richards et al. use VABS on the linearly tapered wing to find the
        # spanwise properties
        alpha = np.linspace(0, 1, self.n_elem_wing)
        for i_elem in range(0, self.n_elem_wing):
            base_stiffness[i_elem + 1, :, :] = stiffness_root*(1-alpha[i_elem]**2) + stiffness_tip*alpha[i_elem]**2
        base_stiffness[0] = base_stiffness[1]

        # Mass variation along the span
        # Right wing centre of mass - wrt to 0.25c
        cm = (root_airfoil.centre_mass * root_airfoil.mass + root_i_beam.centre_mass * root_i_beam.mass) \
             / np.sum(root_airfoil.mass + root_i_beam.mass)
        cg = np.array([0, -(cm[0] + 0.25 * self.c_root - self.main_ea_root), 0]) * 1
        n_mass = self.n_mass
        # sigma_mass = 1.25

        # Initialise database
        base_mass = self.base_mass

        mass_root_right = np.diag([mu_0, mu_0, mu_0, j_xx, j_yy, j_zz]) * sigma_mass
        mass_root_right[:3, -3:] = -algebra.skew(cg) * mu_0
        mass_root_right[-3:, :3] = algebra.skew(cg) * mu_0

        mass_root_left = np.diag([mu_0, mu_0, mu_0, j_xx, j_yy, j_zz]) * sigma_mass
        mass_root_left[:3, -3:] = -algebra.skew(-cg) * mu_0
        mass_root_left[-3:, :3] = algebra.skew(-cg) * mu_0

        mass_tip_right = taper_ratio * mass_root_right
        mass_tip_left = taper_ratio * mass_root_left

        ixx_dummy = []
        iyy_dummy = []
        izz_dummy = []

        for i_elem in range(self.n_elem_wing):
            # Create full cross section
            c_bar = self.c_root * ((1-alpha[i_elem]) + self.taper_ratio * alpha[i_elem])
            x_section = WingCrossSection(c_bar)
            print(i_elem)
            print('Section Mass: %.2f ' %x_section.mass)
            print('Linear Mass: %.2f' % (mu_0 * (1-alpha[i_elem]) + mu_0 * self.taper_ratio * alpha[i_elem]))
            print('Section Ixx: %.4f' % x_section.ixx)
            print('Section Iyy: %.4f' % x_section.iyy)
            print('Section Izz: %.4f' % x_section.izz)
            print('Linear Ixx: %.2f' % (j_xx * (1-alpha[i_elem]) + j_xx * self.taper_ratio * alpha[i_elem]))
            # base_mass[i_elem, :, :] = mass_root_right*(1-alpha[i_elem]) + mass_tip_right*alpha[i_elem]
            # base_mass[i_elem + self.n_elem_wing + self.n_elem_fuselage - 1] = mass_root_left*(1-alpha[i_elem]) + mass_tip_left*alpha[i_elem]

            base_mass[i_elem, :, :] = np.diag([x_section.mass, x_section.mass, x_section.mass,
                                               x_section.ixx, x_section.iyy, x_section.izz])
            cg = np.array([0, -(x_section.centre_mass[0] + (0.25 - self.main_ea_root) * c_bar / self.c_root), 0]) * 1
            base_mass[i_elem, :3, -3:] = -algebra.skew(cg) * x_section.mass
            base_mass[i_elem, -3:, :3] = algebra.skew(cg) * x_section.mass

            base_mass[i_elem + self.n_elem_wing + self.n_elem_fuselage - 1, :, :] = np.diag([x_section.mass, x_section.mass, x_section.mass,
                                               x_section.ixx, x_section.iyy, x_section.izz])
            cg = np.array([0, -(x_section.centre_mass[0] + (0.25 - self.main_ea_root) * c_bar / self.c_root), 0]) * 1
            base_mass[i_elem + self.n_elem_wing + self.n_elem_fuselage - 1, :3, -3:] = -algebra.skew(-cg) * x_section.mass
            base_mass[i_elem + self.n_elem_wing + self.n_elem_fuselage - 1, -3:, :3] = algebra.skew(-cg) * x_section.mass

            ixx_dummy.append(x_section.ixx)
            iyy_dummy.append(x_section.iyy)
            izz_dummy.append(x_section.izz)

            # for item in x_section.items:
            #     plt.plot(item.y, item.z)
            # plt.scatter(x_section.centre_mass[0], x_section.centre_mass[1])
            # plt.show()
            # print(x_section.centre_mass)
            # print(cg)

        # plt.plot(range(self.n_elem_wing), ixx_dummy)
        # plt.plot(range(self.n_elem_wing), iyy_dummy)
        # plt.plot(range(self.n_elem_wing), izz_dummy)
        # plt.show()
        # Lumped mass initialisation
        lumped_mass_nodes = self.lumped_mass_nodes
        lumped_mass = self.lumped_mass
        lumped_mass_inertia = self.lumped_mass_inertia
        lumped_mass_position = self.lumped_mass_position

        # Lumped masses nodal position
        # 0 - Right engine
        # 1 - Left engine
        # 2 - Fuselage
        lumped_mass_nodes[0] = 2
        lumped_mass_nodes[1] = n_node_fuselage + n_node_wing + 1
        lumped_mass_nodes[2] = 0

        # Lumped mass value from Richards 2013
        lumped_mass[0:2] = 51.445 / 9.81
        lumped_mass[2] = 150 / 9.81
        # lumped_mass_position[2] = [0, 0, -10.]

        # Lumped mass inertia
        lumped_mass_inertia[0, :, :] = np.diag([0.29547, 0.29322, 0.29547])
        lumped_mass_inertia[1, :, :] = np.diag([0.29547, 0.29322, 0.29547])
        lumped_mass_inertia[2, :, :] = np.diag([0.5, 1, 1]) * lumped_mass[2]

        # Define class attributes
        self.lumped_mass = lumped_mass * 1
        self.lumped_mass_nodes = lumped_mass_nodes * 1
        self.lumped_mass_inertia = lumped_mass_inertia * 1
        self.lumped_mass_position = lumped_mass_position * 1

        self.base_stiffness = base_stiffness
        self.base_mass = base_mass

class CrossSection(object):
    def __init__(self):
        self.rho = 2770

        self.rotation_axes = np.array([0, 0.33-0.25, 0])

        self.y = np.ndarray((2,))
        self.z = np.ndarray((2,))
        self.t = np.ndarray((2,))

    @property
    def mass(self):
        """
        Mass of the I beam per unit length
        """
        return np.sum(self.t * self.elem_length) * self.rho

    @property
    def ixx(self):
        ixx_ = np.sum(self.elem_length * self.t * self.rho * (self.elem_cm_y ** 2 + self.elem_cm_z ** 2))
        return ixx_ + self.mass * (self.centre_mass[0] - self.rotation_axes[1]) ** 2

    @property
    def elem_length(self):
        elem_length = np.sqrt(np.diff(self.y) ** 2 + np.diff(self.z) ** 2)
        return elem_length

    @property
    def elem_cm_y(self):
        elem_cm_y_ = np.ndarray((self.n_elem, ))
        elem_cm_y_[:] = 0.5 * (self.y[:-1] + self.y[1:])
        return elem_cm_y_

    @property
    def elem_cm_z(self):
        elem_cm_z_ = np.ndarray((self.n_elem, ))
        elem_cm_z_[:] = 0.5 * (self.z[:-1] + self.z[1:])
        return elem_cm_z_

    @property
    def centre_mass(self):
        y_cm = np.sum(self.elem_cm_y * self.elem_length) / np.sum(self.elem_length)
        z_cm = np.sum(self.elem_cm_z * self.elem_length) / np.sum(self.elem_length)

        return np.array([y_cm, z_cm])

    @property
    def iyy(self):
        x_dom = np.linspace(-0.5, 0.5, 100)
        x_cg = 0.5 * (x_dom[:-1].copy() + x_dom[1:].copy())
        dx = np.diff(x_dom)[0]
        iyy_ = 0
        for elem in range(len(self.elem_length)):
            z_cg = np.ones_like(x_cg) * self.elem_cm_z[elem]
            iyy_ += np.sum(self.elem_length[elem] * self.t[elem] * dx * self.rho * (x_cg ** 2 + z_cg ** 2))
        return iyy_ #np.sum(self.elem_length * self.t * self.rho * 1 * self.elem_cm_z ** 2)

    @property
    def izz(self):
        x_dom = np.linspace(-0.5, 0.5, 100)
        x_cg = 0.5 * (x_dom[:-1].copy() + x_dom[1:].copy())
        dx = np.diff(x_dom)[0]
        iyy_ = 0
        izz_ = 0
        for elem in range(len(self.elem_length)):
            y_cg = np.ones_like(x_cg) * self.elem_cm_y[elem]
            izz_ += np.sum(self.elem_length[elem] * self.t[elem] * dx * self.rho * (x_cg ** 2 + y_cg ** 2))
        return izz_ #np.sum(self.elem_length * self.t * self.rho * 1 * self.elem_cm_y ** 2)

    @property
    def n_node(self):
        return self.y.shape[0]

    @property
    def n_elem(self):
        return self.n_node - 1

    def build(self, y, z, t):
        self.y = y
        self.z = z
        self.t = t

class IBeam(CrossSection):

    def build(self, c_root):
        t_skin = 0.127e-2
        t_c = 0.12
        w_I = 10e-2 * c_root  # Width of the Ibeam

        self.rho = 2770

        self.y = np.ndarray((self.n_node, ))
        self.z = np.ndarray((self.n_node, ))
        self.t = np.ndarray((self.n_node, ))

        z_max = t_c * c_root
        y = np.array([-w_I/2, w_I/2, 0, 0, -w_I/2, w_I/2])
        z = np.array([z_max/2, z_max/2, z_max/2, -z_max/2, -z_max/2, -z_max/2])
        t = np.array([t_skin, 0, t_skin, 0, t_skin])

        self.y = y
        self.z = z
        self.t = t

class Airfoil(CrossSection):

    def build(self, c_root):
        t_c = 0.12
        t_skin = 0.127e-2 * 1.5

        y_dom = np.linspace(0, c_root, 100)
        y = np.concatenate((y_dom, y_dom[:-1][::-1]))
        z_dom = 5 * t_c * (0.2969 * np.sqrt(y_dom/c_root) -
                                              0.1260 * y_dom/c_root -
                                              0.3516 * (y_dom/c_root) ** 2 +
                                              0.2843 * (y_dom/c_root) ** 3 -
                                              0.1015 * (y_dom/c_root) ** 4) * c_root
        z = np.concatenate((z_dom, -z_dom[:-1][::-1]))

        self.y = y - 0.25 * c_root
        self.z = z
        self.t = t_skin * np.ones(self.n_elem)

class WingCrossSection:

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

        self.items = list()

        self.items.append(Airfoil())
        self.items.append(IBeam())

        for item in self.items:
            item.build(chord)

    @property
    def mass(self):
        return np.sum([item.mass for item in self.items])

    @property
    def centre_mass(self):
        y = np.sum([item.mass * item.centre_mass[0] for item in self.items]) / self.mass
        z = np.sum([item.mass * item.centre_mass[1] for item in self.items]) / self.mass
        return np.array([y, z])

    @property
    def ixx(self):
        return np.sum([item.ixx for item in self.items])

    @property
    def iyy(self):
        return np.sum([item.iyy for item in self.items])

    @property
    def izz(self):
        return np.sum([item.izz for item in self.items])

if __name__ == '__main__':

    ws = Baseline(M=4,
                  N=11,
                  Mstarfactor=5,
                  u_inf=28,
                  rho=1.225,
                  alpha_deg=4)

    # ws.clean_test_files()
    ws.update_mass_stiffness()
    ws.update_fem_prop()
    # ws.generate_fem_file()
    ws.update_aero_properties()
    # ws.generate_aero_file()
    # ws.set_default_config_dict()


text:
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
from collections import deque
from heap import *

class GraphNode:
    def __init__(self, idx=-1):
        self.idx = idx
        self.edges = [] # List of nodes that are adjacent
        self.dists = [] # Parallel list of distances of the edges to these nodes
        self.data = {}
        self.visited = False
        self.touched = False

    def __str__(self):
        return "{}".format(self.idx)

def draw_2d_graph(nodes, edges, draw_nodes=True, draw_labels=True, linewidth=1):
    ax = plt.gca()
    ax.set_facecolor((0.9, 0.9, 0.9))
    for (i, j, d) in edges:
        x1, y1 = nodes[i].data['x'], nodes[i].data['y']
        x2, y2 = nodes[j].data['x'], nodes[j].data['y']
        plt.plot([x1, x2], [y1, y2], linewidth=linewidth, c='k')
    for i, n in enumerate(nodes):
        if draw_nodes:
            sz = 50
            if n.visited:
                sz = 100
            c = 'g'
            if n.touched:
                if n.visited:
                    c = 'k'
                else:
                    c = 'r'
            plt.scatter(n.data['x'], n.data['y'], sz, c=c)
            if draw_labels:
                plt.text(n.data['x']+0.002, n.data['y']+0.002, "{}".format(i), zorder=10, c=c, fontsize='large')

def draw_grid_graph(nodes, edges, draw_nodes=True, draw_labels=True, draw_diags=True):
    N = int(np.sqrt(len(nodes)))
    ax = plt.gca()
    ax.set_facecolor((0.9, 0.9, 0.9))
    # Draw horizontal edges
    for i in range(N):
        plt.plot([i, i], [0, N-1], 'k')
        plt.plot([0, N-1], [i, i], 'k')
    # Draw diagonal edges
    if draw_diags:
        for i in range(N):
            plt.plot([0, i], [i, 0], 'k')
            plt.plot([N-1, i], [i, N-1], 'k')
            plt.plot([0, N-i-1], [i, N-1], 'k')
            plt.plot([N-1, N-i-1], [i, 0], 'k')
    X = []
    colors = []
    sizes = []
    if not draw_nodes:
        return

    for i, n in enumerate(nodes):
        sz = 50
        if n.visited:
            sz = 100
        c = [0, 0.5, 0]
        if n.touched:
            if n.visited:
                c = [0, 0, 0]
            else:
                c = [1, 0, 0]
        colors.append(c)
        sizes.append(sz)
        X.append([n.data['x'], n.data['y']])
    X = np.array(X)
    colors = np.array(colors)
    sizes = np.array(sizes)
    plt.scatter(X[:, 0], X[:, 1], s=sizes, c=colors, zorder=10)
    plt.axis("off")

def make_delaunay_graph(N):
    x = np.random.rand(N)
    y = np.random.rand(N)
    nodes = []
    for i in range(N):
        n = GraphNode(i)
        n.data = {'x':x[i], 'y':y[i]}
        nodes.append(n)
    tri = Delaunay(np.array([x, y]).T).simplices
    edges = set()
    for i in range(tri.shape[0]):
        for k in range(3):
            i1, i2 = tri[i, k], tri[i, (k+1)%3]
            nodes[i1].edges.append(nodes[i2])
            nodes[i2].edges.append(nodes[i1])
            d = np.sqrt(np.sum((x[i1]-x[i2])**2 + (y[i1]-y[i2])**2))
            edges.add((i1, i2, d))
    return nodes, list(edges)

def make_grid_graph(N, seed = 0, include_diags=True):
    """
    Parameters
    ----------
    N: int
        Resolution of grid
    """
    np.random.seed(seed)
    nodes = []
    for i in range(N):
        for j in range(N):
            n = GraphNode(i*N+j)
            n.data = {'x':j, 'y':i}
            nodes.append(n)
    edges = []
    neighbs = []
    for di in range(-1, 2):
        for dj in range(-1, 2):
            if di != 0 or dj != 0:
                if include_diags or abs(di)+abs(dj) < 2:
                    neighbs.append([di, dj])
    for i in range(N):
        for j in range(N):
            idx1 = i*N + j
            for [di, dj] in neighbs:
                ii = i + di
                jj = j + dj
                if ii >= 0 and jj >= 0 and ii < N and jj < N:
                    idx2 = ii*N + jj
                    x1, y1 = nodes[idx1].data['x'], nodes[idx1].data['y']
                    x2, y2 = nodes[idx2].data['x'], nodes[idx2].data['y']
                    dx = x1-x2
                    dy = y1-y2
                    d = (dx**2 + dy**2)**0.5
                    edges.append((idx1, idx2, d))
                    nodes[idx1].edges.append(nodes[idx2])
                    nodes[idx1].dists.append(d)
                    nodes[idx2].edges.append(nodes[idx1])
                    nodes[idx2].dists.append(d)
    return nodes, edges

def get_collection_str(coll):
    s = "["
    for i, n in enumerate(list(coll)):
        if (i+1) % 20 == 0:
            s += "\n"
        s += "%i"%n.idx
        if i < len(coll)-1:
            s += ","
    s += "]"
    return s

def dist_of_edge(e):
    return e[2]

def get_mst_kruskal(nodes, edges):
    from unionfind import UFFast
    edges = sorted(edges, key = dist_of_edge)
    djset = UFFast(len(nodes))
    new_edges = []
    for e in edges:
        (i, j, d) = e
        if not djset.find(i, j):
            djset.union(i, j)
            new_edges.append(e)
    # Update node adjacency
    for n in nodes:
        n.edges = []
        n.dists = []
    for (i, j, d) in new_edges:
        nodes[i].edges.append(nodes[j])
        nodes[i].dists.append(d)
        nodes[j].edges.append(nodes[i])
        nodes[i].dists.append(d)
    return new_edges

def do_bfs(nodes, edges, draw_labels=True, plot_fn=draw_2d_graph):
    queue = deque()
    nodes[0].touched = True
    queue.append(nodes[0])
    plt.figure(figsize=(8, 6))
    idx = 0
    while len(queue) > 0:
        plt.clf()
        plot_fn(nodes, edges, draw_labels=draw_labels)
        n = queue.popleft()
        queue.appendleft(n)
        if draw_labels:
            plt.title("Processing {}\n{}".format(n.idx, get_collection_str(queue)))
        plt.savefig("{}.png".format(idx), bbox_inches='tight')
        n = queue.popleft()
        n.visited = True
        for e in n.edges:
            if not e.visited and not e.touched:
                e.touched = True
                queue.append(e)
        idx += 1
    plt.clf()
    plot_fn(nodes, edges, draw_labels=draw_labels)
    plt.savefig("{}.png".format(idx), bbox_inches='tight')

def do_dijkstra(nodes, edges, draw_labels=True, plot_fn=draw_2d_graph):
    for n in nodes:
        n.touched = False
        n.visited = False
    queue = HeapTree()
    nodes[0].touched = True
    queue.push((0, nodes[0]))
    plt.figure(figsize=(8, 6))
    idx = 0
    while len(queue) > 0:
        (d, n) = queue.peek()
        if n.visited:
            queue.pop()
            continue
        plt.clf()
        plot_fn(nodes, edges, draw_labels=draw_labels)
        if draw_labels:
            s = get_collection_str([n[1].idx for n in queue._arr])
            plt.title("Processing {} at distance {}\n{}".format(n.idx, d, s))
        plt.savefig("{}.png".format(idx), bbox_inches='tight')
        (d, n) = queue.pop()
        n.visited = True
        for i in range(len(n.edges)):
            ni = n.edges[i]
            if not ni.visited:
                ni.touched = True
                di = n.dists[i]
                queue.push((d+di, ni))
        idx += 1
    plt.clf()
    plot_fn(nodes, edges, draw_labels=draw_labels)
    plt.savefig("{}.png".format(idx), bbox_inches='tight')

def do_dfs(nodes, edges, draw_labels=True, plot_fn=draw_2d_graph):
    stack = [nodes[0]]
    nodes[0].touched = True
    plt.figure(figsize=(8, 6))
    idx = 0
    while len(stack) > 0:
        plt.clf()
        plot_fn(nodes, edges, draw_labels=draw_labels)
        n = stack[-1]
        if draw_labels:
            plt.title("Processing {}\n{}".format(n.idx, get_collection_str(stack)))
        plt.savefig("{}.png".format(idx), bbox_inches='tight')
        n = stack.pop()
        n.visited = True
        for e in n.edges:
            if not e.visited and not e.touched:
                e.touched = True
                stack.append(e)
        idx += 1
    plt.clf()
    plot_fn(nodes, edges, draw_labels=draw_labels)
    plt.savefig("{}.png".format(idx), bbox_inches='tight')

def do_2d_tsp(nodes, edges, start_node = 0, do_plot=True, draw_labels=True):
    """
    Do a depth-first search on the MST
    """
    for n in nodes:
        n.touched = False
        n.visited = False
    edges = get_mst_kruskal(nodes, edges)
    stack = [nodes[start_node]]
    nodes[start_node].touched = True
    if do_plot:
        plt.figure(figsize=(12, 6))
    idx = 0
    path = []
    tsp_edges = []
    total_len = 0
    plot_fn=draw_2d_graph
    while len(stack) > 0:
        n = stack[-1]
        path.append(n)
        if len(path) > 1:
            [n1, n2] = path[-2::]
            dx = n1.data['x'] - n2.data['x']
            dy = n1.data['y'] - n2.data['y']
            d = (dx**2 + dy**2)**0.5
            tsp_edges.append((path[-2].idx, path[-1].idx, d))
            total_len += d
        if do_plot:
            plt.clf()
            plt.subplot(121)
            plot_fn(nodes, edges, draw_labels=draw_labels)
            if draw_labels:
                plt.title("Processing {}\n{}".format(n.idx, get_collection_str(stack)))
            plt.subplot(122)
            plot_fn(nodes, tsp_edges, draw_labels=False)
            plt.title("Total length = {:.3f}".format(total_len))
            plt.savefig("{}.png".format(idx), bbox_inches='tight')
        n = stack.pop()
        n.visited = True

        for e in n.edges:
            if not e.visited and not e.touched:
                e.touched = True
                stack.append(e)
        idx += 1

    # Loop path back around
    path.append(path[0])
    [n1, n2] = path[-2::]
    dx = n1.data['x'] - n2.data['x']
    dy = n1.data['y'] - n2.data['y']
    d = (dx**2 + dy**2)**0.5
    tsp_edges.append((path[-2].idx, path[-1].idx, d))
    total_len += d

    if do_plot:
        plt.clf()
        plt.subplot(121)
        plot_fn(nodes, edges, draw_labels=draw_labels)
        plt.subplot(122)
        plot_fn(nodes, tsp_edges, draw_labels=False)
        plt.title("Total length = {:.3f}".format(total_len))
        plt.savefig("{}.png".format(idx), bbox_inches='tight')

    return path, total_len

if __name__ == '__main__':
    """
    np.random.seed(1)
    nodes, edges = make_delaunay_graph(20)
    do_2d_tsp(nodes, edges, 11)
    """

    #"""
    np.random.seed(0)
    N = 20
    nodes, edges = make_grid_graph(N)
    idx = N*int(N/2) + int(N/2)
    nodes[0], nodes[idx] = nodes[idx], nodes[0]
    do_dijkstra(nodes, edges, draw_labels=False, plot_fn=draw_grid_graph)
    #"""

reason: import numpy,from scipy


text:
import numbers
import random

import numpy as np

class MeshCell(object):

    def __init__(self, coo, bounds = None):
        self.coo      = coo
        self.bounds   = bounds
        self.elements = []
        self.metadata = []

    def add(self, e, metadata = None):
        self.elements.append((e, metadata))

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

    def draw(self, replace = True):
        idx = random.randint(0, len(self.elements) - 1)
        if replace:
            return self.elements[idx]
        else:
            return self.elements.pop(idx)

class MeshGrid(object):

    def __init__(self, bounds, res):
        all(isinstance(e, numbers.Number) for e in bounds)
        self.bounds = bounds
        if isinstance(res, int):
            self.res = len(bounds)*[res]
        else:
            assert len(res) == len(bounds) and all(isinstance(e, int) for e in res)
            self.res = res

        self.dim = len(bounds)
        self._cells = {}
        self._cells[None] = MeshCell(None) # a cell for everything not inside the bounds
        self._size = 0
        self._nonempty_cells = []

    def _coo(self, p):
        assert len(p) == self.dim
        coo = []
        for pi, (si_min, si_max), res_i in zip(p, self.bounds, self.res):
            coo.append(int((pi - si_min)/(si_max - si_min)*res_i))
            if pi == si_max:
                coo[-1] == res_i - 1
            if si_min > pi or si_max < pi:
                return None
        return tuple(coo)

    def __len__(self):
        return self._size

    def add(self, p, metadata = None):
        assert len(p) == self.dim
        self._size += 1
        coo = self._coo(p)
        if not coo in self._cells:
            self._cells[coo] = MeshCell(coo)
        cell = self._cells[coo]
        cell.add(p, metadata)
        if len(cell) == 1:
            self._nonempty_cells.append(cell)

    def draw(self, replace = True):
        """Draw uniformly between existing (non-empty) cells"""
        try:
            idx = random.randint(0, len(self._nonempty_cells) - 1)
        except ValueError:
            raise ValueError("can't draw from an empty meshgrid")
        e, md = self._nonempty_cells[idx].draw(replace = replace)
        if not replace and len(self._nonempty_cells[idx]) == 0:
            self._nonempty_cells.pop(idx)
        return e, md

max_stars_repo_name: afcarl/toolbox

max_issues_repo_name: afcarl/toolbox

max_forks_repo_name: afcarl/toolbox


text:
"""
Copyright 2021 Patrick Schwab, Arash Mehrjou, GlaxoSmithKline plc; Andrew Jesson, University of Oxford; Ashkan Soleymani, MIT

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""
import numpy as np
from typing import AnyStr, List
from slingpy import AbstractDataSource
from slingpy.models.abstract_base_model import AbstractBaseModel
from genedisco.active_learning_methods.acquisition_functions.base_acquisition_function import (
    BaseBatchAcquisitionFunction,
)

class MarginSamplingAcquisition(BaseBatchAcquisitionFunction):
    def __call__(
        self,
        dataset_x: AbstractDataSource,
        select_size: int,
        available_indices: List[AnyStr],
        last_selected_indices: List[AnyStr] = None,
        model: AbstractBaseModel = None,
    ) -> List:
        avail_dataset_x = dataset_x.subset(available_indices)
        model_pedictions = model.predict(avail_dataset_x, return_std_and_margin=True)

        if len(model_pedictions) != 3:
            raise TypeError("The provided model does not output margins.")

        # print("Hi! This is Margin Sampling!")
        pred_mean, pred_uncertainties, pred_margins = model_pedictions
        # print(pred_margins)

        if len(pred_mean) < select_size:
            raise ValueError(
                "The number of query samples exceeds" "the size of the available data."
            )

        numerical_selected_indices = np.flip(np.argsort(pred_margins))[:select_size]
        selected_indices = [available_indices[i] for i in numerical_selected_indices]

        return selected_indices

max_stars_repo_name: PrateekMunjal/genedisco

max_issues_repo_name: PrateekMunjal/genedisco

max_forks_repo_name: PrateekMunjal/genedisco


text:
import numpy, yaml, itertools, h5py, math
import pysweep.utils.validate as validate
import pysweep.tests as tests
import pysweep
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import ticker
from matplotlib import cm
from matplotlib import tri
from mpl_toolkits import mplot3d
from matplotlib.colors import Normalize
import mpi4py
mpi4py.rc.recv_mprobe = False
import mpi4py.MPI as MPI
import os
from scipy.stats import linregress

#Global data
standardSizes = numpy.asarray([160, 320, 480, 640, 800, 960, 1120])
blockSizes = numpy.asarray([8,12,16,20,24,32])
best = []
worst = []
shares = numpy.arange(0,1.1,0.1)
combos = []
cutoffs = [0,]

def switchColorScheme(printParams=False):
    matplotlib.rcParams.update({'text.color' : "white",
                     'axes.labelcolor' : "white","axes.edgecolor":"white","xtick.color":"white","ytick.color":"white","grid.color":"white","savefig.facecolor":"#333333","savefig.edgecolor":"#333333","axes.facecolor":"#333333"})
    if printParams:
        for k in matplotlib.rcParams:
            if "color" in k:
                print(k)

def printArray(a):
    for row in a:
        for item in row:
            print("{:8.3f}".format(col), end=" ")
        print("")

def calcNumberOfPts():
    npts = 0
    global combos, cutoffs
    for ss in standardSizes:
        for bs in blockSizes:
            nblks = numpy.ceil(ss/bs)
            shareset = set()
            for sh in shares:
                numerator = float('%.5f'%(sh*nblks))
                tempShare = numpy.ceil(numerator)/nblks
                # print(tempShare,nblks,sh*nblks)
                # input()
                shareset.add(tempShare)
            shareset = list(shareset)
            shareset.sort()
            for sh in shareset:
                combos.append([ss,bs,sh])
            npts+=len(shareset)
        cutoffs.append(npts)
    return npts

def validateClusterEuler():
    #Plotting
    elev=45
    azim=25
    fig =  plt.figure()
    nrows = 2
    ncols = 3
    axes = [fig.add_subplot(nrows, ncols, i+1) for i in range(nrows*ncols)]
    fig.subplots_adjust(wspace=0.75,hspace=0.8,right=0.8)
    #Physical Colorbar
    cbounds=numpy.linspace(0.4,1,6)
    # cbar_ax = fig.add_axes([0.89, 0.39, 0.05, 0.51]) #left bottom width height
    cbar_ax = fig.add_axes([0.89, 0.11, 0.05, 0.76]) #left bottom width height
    ibounds = numpy.linspace(cbounds[0],cbounds[-1],100)
    cbar = fig.colorbar(cm.ScalarMappable(cmap=cm.inferno),cax=cbar_ax,boundaries=ibounds)
    cbar.ax.set_yticklabels([["{:0.1f}".format(i) for i in cbounds]])
    tick_locator = ticker.MaxNLocator(nbins=len(cbounds))
    cbar.locator = tick_locator
    cbar.update_ticks()
    cbar_ax.set_title('$\\rho(x,y,t)$\n$[kg/m^3]$',y=1.01)
    with h5py.File("./data/eulerOutput.hdf5","r") as f:
        data = f['data']
        times = [0,] #Times to plot
        times.append(numpy.shape(data)[0]//2)
        times.append(numpy.shape(data)[0]-1)
        #Euler stuff
        d = 0.1
        gamma = 1.4
        dx = 10/(data.shape[-1]-1)
        dt = d*dx
        #Numerical
        for i,ax in enumerate(axes[3:6]):
            ax.set_xlabel("X")
            ax.set_ylabel("Y")
            ax.yaxis.labelpad=-1
            ax.set_title('N:$t={:0.2f}$ [s]'.format(times[i]*dt))
            validate.eulerContourAx(ax,data[times[i],0],1,1)
            ax.grid('on',color="k")
        #Analytical
        npx = numpy.shape(data)[2]
        for i,ax in enumerate(axes[:3]):
            ax.set_xlabel("X")
            ax.set_ylabel("Y")
            ax.yaxis.labelpad=-1
            ax.set_title('A:$t={:0.2f}$ [s]'.format(times[i]*dt))
            adata = numpy.zeros((1,4,npx,npx))
            adata[0,:,:,:] = pysweep.equations.euler.getAnalyticalArray(npx,npx,t=(times[i])*dt)
            validate.eulerContourAx(ax,adata[0,0],1,1)
            ax.grid('on',color="k")
    plt.savefig("./plots/eulerValidate.pdf",bbox_inches='tight')
    plt.close('all')

def generateArraySizes():
    """Use this function to generate array sizes based on block sizes."""
    blocksizes = [8, 12, 16, 20, 24, 32] #blocksizes with most options
    arraystart = 0
    divisible =[False,False]
    arraysizes = []
    for i in range(7):
        arraystart += 32*20
        while not numpy.all(divisible):
            arraystart+=blocksizes[-1]
            divisible =[arraystart%bs==0 for bs in blocksizes]
        arraysizes.append(arraystart)
    return arraysizes

def getYamlData(file,equation):
    document = open(file,'r')
    yamlFile = yaml.load(document,Loader=yaml.FullLoader)
    data = []
    for key in yamlFile:
        if equation in yamlFile[key]['cpu']:
            swept = 1 if yamlFile[key]['swept'] else 0
            data.append([float(swept),float(yamlFile[key]['array_shape'][2]),float(yamlFile[key]['blocksize']),float(yamlFile[key]['share']),float(yamlFile[key]['runtime']),float(yamlFile[key]['time_per_step'])])
    i = 0
    while i < len(data):
        j = i+1
        while j < len(data):
            if data[i][:4] == data[j][:4]:
                data.pop(j)
            j+=1
        i+=1

    data = numpy.array(data,)
    indexes = numpy.lexsort((data[:,3],data[:,2],data[:,1],data[:,0]))
    sortedData = numpy.zeros(data.shape)
    for i,idx in enumerate(indexes):
        sortedData[i,:] = data[idx,:]
    return sortedData,standardSizes

def checkForAllPts(sortedData):
    ct = 0
    for j in range(len(sortedData)):

        if list(sortedData[j,1:4])!=combos[(j+ct)%452]:
            print("{:n}, {:n}, {:n}, {:.4f}".format(*sortedData[j,:4]))
            print("{:n}, {:n}, {:.4f}".format(*combos[(j+ct)%452]))
            ct+=1
            input()

def getContourData(data,arraysize,uBlocks,uShares):
    triang = tri.Triangulation(uBlocks, uShares)
    interpolator = tri.LinearTriInterpolator(triang, data)
    S,B = numpy.meshgrid(shares,blockSizes)
    Z = interpolator(B, S)
    return B,S,Z

def makeArrayContours(data,rBlocks,rShares,cbounds,cmap,cbt,fname,switch=False,printer=False,record=False):
    #Make speed up figure
    ai = lambda x: slice(int(data.shape[0]//len(standardSizes)*x),int(data.shape[0]//len(standardSizes)*(x+1)),1)
    fig, axes = plt.subplots(ncols=3,nrows=2)
    fig.subplots_adjust(wspace=0.55,hspace=0.4,right=0.75)
    axes = numpy.reshape(axes,(6,))
    cbar_ax = fig.add_axes([0.85, 0.11, 0.05, 0.77])
    ibounds = numpy.linspace(cbounds[0],cbounds[-1],100)
    cbar = fig.colorbar(cm.ScalarMappable(cmap=cmap),cax=cbar_ax,boundaries=ibounds)
    cbar.ax.set_yticklabels([["{:0.1f}".format(i) for i in cbounds]])
    tick_locator = ticker.MaxNLocator(nbins=len(cbounds))
    cbar.locator = tick_locator
    cbar.update_ticks()
    cbar_ax.set_title(cbt,y=1.01)

    if switch:
        mbc = ('k','w')
    else:
        mbc = ('w','k')
    for i in range(0,len(cutoffs)-2,1):
        lb = cutoffs[i+1]
        ub = cutoffs[i+2]
        xw= numpy.where(numpy.amax(data[lb:ub])==data[lb:ub])
        xb= numpy.where(numpy.amin(data[lb:ub])==data[lb:ub])
        #Contours
        axes[i].tricontour(rBlocks[lb:ub],rShares[lb:ub]*100,data[lb:ub],
                 colors=('k',),linestyles=('-',),linewidths=(0.25,),vmin=cbounds[0],vmax=cbounds[-1])
        axes[i].tricontourf(rBlocks[lb:ub],rShares[lb:ub]*100,data[lb:ub],cmap=cmap,vmin=cbounds[0],vmax=cbounds[-1])
        #Marking best and worst
        axes[i].plot(rBlocks[lb:ub][xb],rShares[lb:ub][xb]*100,linestyle=None,marker='o',markerfacecolor=mbc[0],markeredgecolor=mbc[1],markersize=6)
        axes[i].plot(rBlocks[lb:ub][xw],rShares[lb:ub][xw]*100,linestyle=None,marker='o',markerfacecolor=mbc[1],markeredgecolor=mbc[0],markersize=6)
        if record:
            worst.append((rBlocks[lb:ub][xb],rShares[lb:ub][xb]))
            best.append((rBlocks[lb:ub][xw],rShares[lb:ub][xw]))
        #Labels
        axes[i].set_title('Array Size: ${}$'.format(standardSizes[i+1]))
        axes[i].set_ylabel('Share [%]')
        axes[i].set_xlabel('Block Size')
        axes[i].grid(color='k', linewidth=1)
        axes[i].set_xticks([8,16,24,32])
        axes[i].set_yticks([0,25,50,75,100])
        axes[i].yaxis.labelpad = 0.5
        axes[i].xaxis.labelpad = 0.5
    plt.savefig(fname,bbox_inches='tight')
    plt.close('all')

def getContourYamlData(file,equation):
    data,standardSizes = getYamlData(file,equation)
    standarddata = data[:data.shape[0]//2,:]
    sweptdata = data[data.shape[0]//2:,:]
    speedup = standarddata[:,4]/sweptdata[:,4]
    return data,standarddata,sweptdata,standardSizes,speedup

def getDataLimits(data,npts=10):
    upperLimit = numpy.amax(data)
    lowerLimit = numpy.amin(data)
    upperLimit = math.ceil(upperLimit*10)/10
    lowerLimit = math.floor(lowerLimit*10)/10
    return numpy.linspace(lowerLimit,upperLimit,npts)

def getStudyContour(equation):
    #Get data
    newdata,newstandarddata,newsweptdata,newstandardSizes,newspeedup = getContourYamlData("./newHardware/log.yaml",equation)
    olddata,oldstandarddata,oldsweptdata,oldstandardSizes,oldspeedup = getContourYamlData("./oldHardware/log.yaml",equation)
    #Make contour
    #Speedup contours
    speedLimits = getDataLimits([newspeedup,oldspeedup])
    makeArrayContours(oldspeedup,oldsweptdata[:,2],oldsweptdata[:,3],speedLimits,cm.inferno,'Speedup',"./plots/speedUp{}.pdf".format(equation+"Old"),switch=True,record=True)
    makeArrayContours(newspeedup,newsweptdata[:,2],newsweptdata[:,3],speedLimits,cm.inferno,'Speedup',"./plots/speedUp{}.pdf".format(equation+"New"),switch=True,record=True)
    #Time Contours
    timeLimits = getDataLimits([oldsweptdata[:,4],newsweptdata[:,4],oldstandarddata[:,4],newstandarddata[:,4]])
    makeArrayContours(oldsweptdata[:,4],oldsweptdata[:,2],oldsweptdata[:,3],timeLimits,cm.inferno_r,'Clocktime [s]',"./plots/clockTimeSwept{}.pdf".format(equation+"Old"))
    makeArrayContours(newsweptdata[:,4],newsweptdata[:,2],newsweptdata[:,3],timeLimits,cm.inferno_r,'Clocktime [s]',"./plots/clockTimeSwept{}.pdf".format(equation+"New"))
    makeArrayContours(oldstandarddata[:,4],oldstandarddata[:,2],oldstandarddata[:,3],timeLimits,cm.inferno_r,'Clocktime [s]',"./plots/clockTimeStandard{}.pdf".format(equation+"Old"))
    makeArrayContours(newstandarddata[:,4],newstandarddata[:,2],newstandarddata[:,3],timeLimits,cm.inferno_r,'Clocktime [s]',"./plots/clockTimeStandard{}.pdf".format(equation+"New"))

    #Time per timestep contours
    timePerStepLimits = getDataLimits([oldsweptdata[:,5],newsweptdata[:,5],oldstandarddata[:,5],newstandarddata[:,5]])
    makeArrayContours(oldsweptdata[:,5],oldsweptdata[:,2],oldsweptdata[:,3],timePerStepLimits,cm.inferno_r,'Time/Step [s]',"./plots/timePerStepSwept{}.pdf".format(equation+"Old"))
    makeArrayContours(newsweptdata[:,5],newsweptdata[:,2],newsweptdata[:,3],timePerStepLimits,cm.inferno_r,'Time/Step [s]',"./plots/timePerStepSwept{}.pdf".format(equation+"New"))
    makeArrayContours(oldstandarddata[:,5],oldstandarddata[:,2],oldstandarddata[:,3],timePerStepLimits,cm.inferno_r,'Time/Step [s]',"./plots/timePerStepStandard{}.pdf".format(equation+"Old"))
    makeArrayContours(newstandarddata[:,5],newstandarddata[:,2],newstandarddata[:,3],timePerStepLimits,cm.inferno_r,'Time/Step [s]',"./plots/timePerStepStandard{}.pdf".format(equation+"New"))

    #Speedup based on time per timestep
    new_tpt_speedup = newstandarddata[:,5]/newsweptdata[:,5]
    old_tpt_speedup = oldstandarddata[:,5]/oldsweptdata[:,5]
    tptSpeedLimits = getDataLimits([new_tpt_speedup,old_tpt_speedup])
    makeArrayContours(new_tpt_speedup,newsweptdata[:,2],newsweptdata[:,3],tptSpeedLimits,cm.inferno,'Speedup',"./plots/SpeedupTPT{}.pdf".format(equation+"New"))
    makeArrayContours(old_tpt_speedup,oldsweptdata[:,2],oldsweptdata[:,3],tptSpeedLimits,cm.inferno,'Speedup',"./plots/SpeedupTPT{}.pdf".format(equation+"Old"))

    # makeArrayContours(oldstandarddata[:,5],oldstandarddata[:,2],oldstandarddata[:,3],timePerStepLimits,cm.inferno_r,'Time/Step [s]',"./plots/timePerStepStandard{}.pdf".format(equation+"Old"))
    # makeArrayContours(newstandarddata[:,5],newstandarddata[:,2],newstandarddata[:,3],timePerStepLimits,cm.inferno_r,'Time/Step [s]',"./plots/timePerStepStandard{}.pdf".format(equation+"New"))

    #Hardware Speedup
    hardwareSpeedUp = oldsweptdata[:,5]/newsweptdata[:,5]
    hardwareSpeedUpStand = oldstandarddata[:,5]/newstandarddata[:,5]
    # print("Hardware Speed Up Average {}:{:0.2f}".format(equation,numpy.mean(hardwareSpeedUp)))
    # print("Hardware Speed Up Min {}:{:0.2f}".format(equation,numpy.amin(hardwareSpeedUp)))
    # print("Hardware Speed Up Max {}:{:0.2f}\n".format(equation,numpy.amax(hardwareSpeedUp)))
    speedLimits = getDataLimits([hardwareSpeedUp,])
    # if equation == "heat":
    #     speedLimits = numpy.linspace(1,2,11)
    makeArrayContours(hardwareSpeedUp,oldsweptdata[:,2],oldsweptdata[:,3],speedLimits,cm.inferno,'Speedup',"./plots/hardwareSpeedUp{}.pdf".format(equation),switch=True)
    makeArrayContours(hardwareSpeedUpStand,oldstandarddata[:,2],oldstandarddata[:,3],speedLimits,cm.inferno,'Speedup',"./plots/hardwareStandardSpeedUp{}.pdf".format(equation),switch=True)

    # print("New Speed Up Average {}:{:0.2f}".format(equation,numpy.mean(newspeedup)))
    # print("Old Speed Up Average {}:{:0.2f}".format(equation,numpy.mean(oldspeedup)))
    print("New Speed Up Min,Max {}:{:0.2f},{:0.2f}".format(equation,numpy.amin(new_tpt_speedup),numpy.amax(new_tpt_speedup)))
    print("Old Speed Up Min,Max {}:{:0.2f},{:0.2f}".format(equation,numpy.amin(old_tpt_speedup),numpy.amax(old_tpt_speedup)))

    # 3D plot
    ThreeDimPlot(oldsweptdata[:,2],100*oldsweptdata[:,3],oldspeedup,speedLimits)
    return oldsweptdata,newsweptdata

def ThreeDimPlot(x,y,z,lims):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_trisurf(x,y,z,cmap=cm.inferno)
    ax.set_xticks([8,16,24,32])
    ax.set_yticks([0,25,50,75,100])
    ax.set_zlim([lims[0],lims[-1]])
    ax.yaxis.labelpad = 0.5
    ax.xaxis.labelpad = 0.5
    # plt.show()

def lineregressSlope(d1,d2,name):
    slope, intercept, r_value, p_value, std_err = linregress(d1, d2)
    print(name,slope)
    return slope

def ScalabilityPlots():
    plt.close()
    data,standardSizes = getYamlData("./scalability/scalability.yaml","euler")
    dl = len(data)
    eulersweptdata = data[dl//2:]
    eulerstandarddata = data[:dl//2]
    data,standardSizes = getYamlData("./scalability/scalability.yaml","heat")
    dl = len(data)
    heatsweptdata = data[dl//2:]
    heatstandarddata = data[:dl//2]
    lsty = "--"
    #Euler
    fig, ax = plt.subplots(ncols=1,nrows=1)
    ax.set_ylabel("Time per Timestep [s]")
    ax.set_xlabel("Spatial Points")
    # ax.set_title("Compressible Euler Equations")
    spacedata = eulersweptdata[:,1]*eulersweptdata[:,1] #2580*numpy.arange(1,5,1)#
    # ax.axis([10**6,10**7, 10, 30])
    l1 = ax.loglog(spacedata,eulersweptdata[:,5],marker="o",color='#d95f02')
    l2 = ax.loglog(spacedata,eulerstandarddata[:,5],marker="s",color="#7570b3",linestyle=lsty)
    lrgsw = lineregressSlope(spacedata,eulersweptdata[:,5],"Swept Euler ")
    lrgst = lineregressSlope(spacedata,eulerstandarddata[:,5],"Standard Euler ")
    # ax.annotate('m={:.2f}'.format(lrgsw), xy=(eulersweptdata[-2,1],eulersweptdata[-2,5]), xytext=(eulersweptdata[-2,1]-1000, eulersweptdata[-2,5]+500),arrowprops=dict(facecolor='black', arrowstyle="simple"))
    # ax.annotate('m={:.2f}'.format(lrgst), xy=(eulerstandarddata[-2,1],eulerstandarddata[-2,4]), xytext=(eulerstandarddata[-2,1]-500, eulerstandarddata[-2,4]+1000),arrowprops=dict(facecolor='black', arrowstyle="simple"))
    leg = ax.legend(["Swept","Standard"])
    plt.savefig("./plots/weakScalabilityEuler.pdf",bbox_inches='tight')
    plt.close('all')

    #Heat
    fig, ax = plt.subplots(ncols=1,nrows=1)
    # ax.axis([10**6,10**7, 3* 10**-1, 10**0])
    ax.yaxis.set_label_coords(-0.135,0.5)
    ax.set_ylabel("Time per Timestep [s]")
    ax.set_xlabel("Spatial Points")
    # ax.set_title("Heat Diffusion Equation")
    l1 = ax.loglog(spacedata,heatsweptdata[:,5],marker="o",color='#d95f02')
    l2 = ax.loglog(spacedata,heatstandarddata[:,5],marker="s",color="#7570b3",linestyle=lsty)
    lrgsw = lineregressSlope(spacedata,heatsweptdata[:,5],"Swept Heat ")
    lrgst = lineregressSlope(spacedata,heatstandarddata[:,5],"Standard Heat ")
    # ax.annotate('m={:.2f}'.format(lrgsw), xy=(heatsweptdata[-2,1],heatsweptdata[-2,5]), xytext=(heatsweptdata[-2,1]-1000, heatsweptdata[-2,5]),arrowprops=dict(facecolor='black', arrowstyle="simple"))
    # ax.annotate('m={:.2f}'.format(lrgst), xy=(heatstandarddata[-2,1],heatstandarddata[-2,4]), xytext=(heatstandarddata[-2,1], heatstandarddata[-2,4]-30),arrowprops=dict(facecolor='black', arrowstyle="simple"))
    leg = ax.legend(["Swept","Standard"])
    plt.savefig("./plots/weakScalabilityHeat.pdf",bbox_inches='tight')
    plt.close('all')
    # plt.show()

def getOccurences(vals,targets):
    L = len(vals)
    return [vals.count(i)/L*100 for i in targets]

def modePlots():
    global best, worst
    best = [(i[0],round(j[0]*10,0)/10) for i, j in best]
    worst = [(i[0],round(j[0]*10,0)/10) for i, j in worst]

    bBlocks,bShares = zip(*best)
    wBlocks,wShares = zip(*worst)
    bestBlockOccurences = getOccurences(bBlocks,blockSizes) #[:len(bBlocks//2)]
    worstBlockOccurences = getOccurences(wBlocks,blockSizes)
    bestShareOccurences = getOccurences(bShares,shares)
    worstShareOccurences = getOccurences(wShares,shares)
    fig =  plt.figure()
    nrows = 1
    ncols = 2
    axes = [fig.add_subplot(nrows, ncols, i+1) for i in range(nrows*ncols)]
    fig.subplots_adjust(wspace=0.3)

    ax1,ax2 = axes

    ax1.set_ylabel("Percentage Of Occurences[%]")
    ax1.set_xlabel("Block Size")
    ax1.plot(blockSizes,bestBlockOccurences,marker="o",color='#d95f02')
    ax1.plot(blockSizes,worstBlockOccurences,marker="o",color="#7570b3")
    leg = ax1.legend(["Best","Worst"])

    ax2.set_ylabel("Percentage Of Occurences[%]")
    ax2.set_xlabel("Share")
    ax2.plot(shares,bestShareOccurences,marker="o",color='#d95f02')
    ax2.plot(shares,worstShareOccurences,marker="o",color="#7570b3")
    leg = ax2.legend(["Best","Worst"])
    plt.savefig("./plots/caseModes.pdf",bbox_inches='tight')

def validateHeat(adata,ndata,times):
    #Plotting
    fig =  plt.figure()
    nrows = 2
    ncols = 3
    axes = [fig.add_subplot(nrows, ncols, i+1) for i in range(nrows*ncols)]
    fig.subplots_adjust(wspace=0.75,hspace=0.8,right=0.8)
    #Physical Colorbar
    cbounds=numpy.linspace(0.4,1,6)
    # cbar_ax = fig.add_axes([0.89, 0.39, 0.05, 0.51]) #left bottom width height
    cbar_ax = fig.add_axes([0.89, 0.11, 0.05, 0.76]) #left bottom width height
    ibounds = numpy.linspace(cbounds[0],cbounds[-1],100)
    cbar = fig.colorbar(cm.ScalarMappable(cmap=cm.inferno),cax=cbar_ax,boundaries=ibounds)
    cbar.ax.set_yticklabels([["{:0.1f}".format(i) for i in cbounds]])
    tick_locator = ticker.MaxNLocator(nbins=len(cbounds))
    cbar.locator = tick_locator
    cbar.update_ticks()
    cbar_ax.set_title('$T(x,y,t) [K]$',y=1.01)
    #Numerical stuff
    for i,ax in enumerate(axes[3:6]):
        ax.set_xlabel("X")
        ax.set_ylabel("Y")
        ax.yaxis.labelpad=-1
        ax.set_title('N:$t={:0.2f}$ [s]'.format(times[i]))
        validate.heatContourAx(ax,ndata[i],1,1)

    #Analytical
    npx = numpy.shape(adata)[-1]
    for i,ax in enumerate(axes[:3]):
        ax.set_xlabel("X")
        ax.set_ylabel("Y")
        ax.yaxis.labelpad=-1
        ax.set_title('N:$t={:0.2f}$ [s]'.format(times[i]))
        validate.heatContourAx(ax,adata[i],1,1)
    plt.savefig("./plots/heatValidate.pdf",bbox_inches='tight')
    plt.close('all')

def PresentationHDE():
    #MPI stuff
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()  #current rank

    if rank==0:
        validate.switchColorScheme() #Change color scheme

    npx = 120
    tsteps = 2000
    alpha = 0.01
    tf,dt = pysweep.equations.heat.getFinalTime(npx,tsteps,alpha=alpha)
    gifmod=100
    adata = numpy.zeros((3,npx,npx))
    ndata = numpy.zeros((3,npx,npx))
    ndata = numpy.zeros((3,npx,npx))
    times = numpy.zeros((3,))
    #create conditions
    pysweep.equations.heat.createInitialConditions(npx,npx)

    #Numerical
    tests.testHeatForwardEuler(npx,timesteps=tsteps,runTest=False)

    #Numerical plots
    with h5py.File("testingHeatFE.hdf5","r",driver="mpio",comm=comm) as f:
        data = f['data']
        finalsteps = data.shape[0]
        if rank == 0:
            validate.createSurface(data[:,0,:,:],0,1,1,1,gif=True,gmod=gifmod,filename="/home/anthony-walker/nrg-swept-project/pysweep-git/study/plots/numericalHeatSurface.pdf")
            validate.createContourf(data[:,0,:,:],0,1,1,1,gif=True,gmod=gifmod,filename="/home/anthony-walker/nrg-swept-project/pysweep-git/study/plots/numericalHeatContour.pdf")
            for i,ts in enumerate([finalsteps//10,finalsteps//2,finalsteps-1]):
                ndata[i,:,:] = data[ts,0,:,:]
        comm.Barrier()

    #Analytical solution
    if rank == 0:
        points = [(int(i),k) for i,k in enumerate(numpy.linspace(0,tf,tsteps))]
        split_points = numpy.array_split(points,comm.Get_size())
    else:
        split_points = None
    split_points = comm.scatter(split_points)

    #Analytical plots
    with h5py.File('temp.hdf5', 'w', driver="mpio",comm=comm) as f:
        data = f.create_dataset("data", (tsteps,1,npx,npx))
        for i,t in split_points:
            T,x,y = pysweep.equations.heat.analytical(npx,npx,t,alpha=0.01)
            data[int(i),0,:,:] = T[0,:,:]
        if rank==0:
            validate.createSurface(data[:,0,:,:],0,1,1,1,gif=True,gmod=gifmod,filename="/home/anthony-walker/nrg-swept-project/pysweep-git/study/plots/analyticalHeatSurface.pdf")
            validate.createContourf(data[:,0,:,:],0,1,1,1,gif=True,gmod=gifmod,filename="/home/anthony-walker/nrg-swept-project/pysweep-git/study/plots/analyticalHeatContour.pdf")
            for i,ts in enumerate([finalsteps//10,finalsteps//2,finalsteps-1]):
                adata[i,:,:] = data[ts,0,:,:]
                times[i] =  dt*ts
    comm.Barrier()

    #Clean up and other plots
    if rank == 0:
        validateHeat(adata,ndata,times)
        os.system("rm temp.hdf5 heatConditions.hdf5 testingHeatFE.hdf5")

def generateOtherEulerGIFs():
    npx = 48
    tsteps = 100
    gmod = 5
    mul=3
    data = numpy.zeros((int(mul*tsteps),4,npx,npx))
    sdata = numpy.zeros((tsteps,4,npx,npx))
    for i,t in enumerate(numpy.linspace(0,mul,mul*tsteps)):
        data[i,:,:,:] = pysweep.equations.euler.getAnalyticalArray(npx,npx,t)
    for i,t in enumerate(numpy.linspace(0,1,tsteps)):
        sdata[i,:,:,:] = pysweep.equations.euler.getPeriodicShock(npx,t)
    validate.createContourf(data[:,0,:,:],0,5,5,1,gif=True,gmod=mul*gmod,filename="analyticalVortex.pdf",LZn=0.4)
    validate.createSurface(sdata[:,0,:,:],0,1,1,1,gif=True,gmod=gmod,filename="analyticalShock.pdf")

def generateEuler(comm,npx,tsteps,gifmod):
    #Numerical
    rank = comm.Get_rank()  #current rank
    pysweep.equations.euler.createInitialConditions(npx,npx) #create initial conditions
    tests.testEulerVortex(npx,timesteps=tsteps,runTest=False) #Run test

    with h5py.File('testingVortex.hdf5', 'r', driver="mpio",comm=comm) as f:
        data = f['data']
        finalsteps = data.shape[0]
        if rank == 0:
            validate.createContourf(data[:,0,:,:],0,5,5,1,gif=True,gmod=gifmod,filename="/home/anthony-walker/nrg-swept-project/pysweep-git/study/plots/numericalVortex{:d}.pdf".format(npx),LZn=0.4,gridBool=True)
    comm.Barrier()

def PresentationEuler():
    """Use this function to produce presentation data"""
    #MPI stuff
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()  #current rank

    npx1 = 120
    tsteps = 1000
    gamma=1.4
    d=0.1
    tf,dt = pysweep.equations.euler.getFinalTime(npx1,tsteps,gamma=gamma)

    if rank==0:
        validate.switchColorScheme() #Change color scheme
    sizes = [npx1,360,720]
    for npx in sizes:
        dx = 10/npx
        dt = d*dx
        tsteps = tf//dt
        gifmod = int((tsteps//20))
        generateEuler(comm,npx,tsteps,gifmod)

    #Clean up and other plots
    if rank == 0:
        os.system("eulerConditions.hdf5 testingVortex.hdf5")

if __name__ == "__main__":
    #Find number of points there should be
    calcNumberOfPts()

    #Presentation bool
    pres = False
    #Switch colors to match presentation
    if pres:
        switchColorScheme()

    # Produce contour plots
    sweptEuler = getStudyContour('heat')
    sweptEuler = getStudyContour('euler')

    # Scalability
    ScalabilityPlots()

    # # validateClusterEuler()

    # # Produce presentation figures
    # if pres:
    #     PresentationEuler()
    #     PresentationHDE()
    #     generateOtherEulerGIFs()
    # # mode plots of best and worst cases
    # modePlots()

max_stars_repo_name: anthony-walker/pysweep

max_issues_repo_name: anthony-walker/pysweep

max_forks_repo_name: walkanth/pysweep


text:
"""Helper classes and functions for cl_performance_and_correlation.py.

Available functions:
- get_auc: returns an np.ndarray of auc values evenly spaced between
    specified AUC limits.
- extra_upper_diagonal: return the extra upper diagonal elements of the
    square symmetric matrix as a ndarray
- compute_mean_cond_corr_matrix: compute the average class conditional
    correlation matrix from rank data and class labels

Available Classes:
- ReplicateData: helper class for storing replicate auc values, and computing
    the mean and sem
- DataForPlots: organize data for plotting
- StatsTable: helper class for organize mean and sem values given parameter
    indexes
"""

import numpy as np

def get_auc(m, auc_limits):
    """Gustavo's method for generating uniformly distributed auc values.

    Args:
        m: (int) > 2, representing the number of base classifiers
        auc_limits: (len(2) list, tuple, or ndarray, i.e. Sequence)
            representing the (minimum AUC, maximum AUC) among the
            set of base classifiers. Each AUC value must be on
            interval (0, 1).

    Returns:
        ((m,) ndarray) of auc values.
    """
    if m < 2:
        raise ValueError(("The number of base classifiers, m, "
                        "must be greater than 2."))
    if len(auc_limits) != 2:
        raise ValueError("Must supply len(2) sequence of AUC values.")

    for auc_lim in auc_limits:
        if auc_lim < 0 or auc_lim > 1:
            raise ValueError("AUC values must be on interval [0,1].")

    auc_range = auc_limits[1] - auc_limits[0]

    auc = np.zeros(m)

    for i in range(m):
        auc[i] = auc_limits[0] + auc_range * i / (m-1)

    return auc

def extra_upper_diagonal(A):
    """Retrieve upper extra diagonal elements from a square symmetric matrix.

    Args:
        A: ((M, M) ndarray) representing a symmetric matrix

    Returns:
        out: ((M(M-1)/2,) ndarray) of values
    """
    if A.ndim != 2:
        raise ValueError("Input must be a square symmetric matrix.")
    elif A.shape[0] != A.shape[1]:
        raise ValueError("Input must be a square symmetric matrix.")
    if (A != A.T).any():
        raise ValueError("Input must be a square symmetric matrix.")

    M = A.shape[0]

    out = np.zeros(int(M*(M-1)/2))
    l = 0
    for i in range(M):
        for j in range(i+1, M):
            out[l] = A[i, j]
            l += 1

    return out

def compute_mean_cond_corr_matrix(R, y):
    """Compute and return the mean conditional correlation matrix.

    Args:
        R: ((M, N) ndarray) independent rows of sample ranks, no ties in row
        y: ((N,) ndarray) of sample class labels [0, 1].

    Returns:
        C: ((M, M) ndarray) average conditional correlation matrix
    """
    C = np.corrcoef(R[:, y == 0])
    C += np.corrcoef(R[:, y == 1])
    C = C / 2

    M = C.shape[0]

    for i in range(M):
        for j in range(i+1, M):
            C[j, i] = C[i, j]

    return C

class ReplicateData:
    """Store and compute statistics of experimental replicates.

    Args:
        nreps: (int) number of replicates
    """
    def __init__(self, nreps):
        if nreps < 2:
            raise ValueError("Number of replicates should be >= 2")
        self._size = nreps
        self._next = 0
        self._data = np.zeros(nreps)

    def __len__(self):
        return self._next

    @property
    def mean(self):
        """Mean over replicate values.

        Returns:
            (float) of the mean

        Raises:
            ZeroDivisionError if there are no replicate data entries
        """
        if self._next == 0:
            raise ZeroDivisionError("No data entries to compute mean.")
        elif self._next < self._size:
            print("Mean computed for less than total replicate number.")
            return np.mean(self._data[:self._next])

        return np.mean(self._data)

    @property
    def sem(self):
        """Standard error of mean over replicate values.

        The standard error of mean is computed as:

        sem = (unbiased standard deviation) / sqrt(N replicates)

        Returns:
            (float) of the standard error of mean

        Raises:
            ZeroDivisionError if there are either 0 or 1 replicate data
                entry.
        """
        if self._next < 2:
            raise ZeroDivisionError("Require >= 2 data points to compute S.E.M.")
        elif self._next < self._size:
            print("SEM computed for less than total replicate number.")
            return np.std(self._data[:self._next], ddof=1) / np.sqrt(len(self))

        return np.std(self._data[:self._next], ddof=1) / np.sqrt(len(self))

    def reset(self):
        self._data.fill(0)
        self._next = 0

    def append(self, val):
        """Add data value to replicates.

        Args:
            val: (int or float) to be added to replicate data set
        """
        if not isinstance(val, int) and not isinstance(val, float):
            raise ValueError("Input data must be an int or float.")

        self._data[self._next] = val

        self._next += 1

class DataForPlots:
    """Organize data for plotting.

    Args:
        label: (str) label to be used in figure legend
        color: (any matplotlib compatible color argument),
            designates the line and marker color associated
            with these data.
    """
    def __init__(self, label, color):
        self.label = label
        self.color = color

        self._data_attributes = ("x", "xerr", "y", "yerr")

    def reset_data(self):
        """Restore all attributes specified in _data_attributes to None."""
        for attribute in self._data_attributes:
            setattr(self, attribute, None)

    def update_data(self, xerr=None, yerr=None, *,x,y):
        """Update attribute values specified in _data_attributes.

        Args:
            x: (N length sequence) x-axis values for plots,
            y: (N length sequence) y-axis values for plots,
            xerr: (N length sequence) optional, containing error along x-axis
            yerr: (N length sequence) optional, containing error along y-axis
        """
        self.reset_data()

        self.x = x
        self.y = y

        self.xerr=xerr
        self.yerr=yerr

class StatsTable:
    """Store and retrieve mean and sem of data.

    Properties:
    - mean: returns ndarray of means
    - sem: returns ndarray of standard error of means (s.e.m.)

    Available Methods:
    - update: add mean and sem over replicates to the position
           specified by given index
    """
    def __init__(self, row_size, col_size):
        self._means = np.zeros(shape=(row_size, col_size))
        self._sem = np.zeros(shape=(row_size, col_size))

    @property
    def mean(self):
        """Return means as ndarray"""
        return self._means

    @property
    def sem(self):
        """Return standard error of means as ndarray"""
        return self._sem

    def update(self, rep_obj, idx):
        """Add mean and sem to corresponding ndarray at idx.

        Args:
            rep_obj: Instance of ReplicateData
            idx: tuple of indices to store mean and sem values.
        """
        self._means[idx] = rep_obj.mean
        self._sem[idx] = rep_obj.sem


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# This code includes modifications of the Vispy example at:
# https://github.com/vispy/vispy/tree/5098146746c098c7af9e4e09ea44ff6ef5ea2187/examples/demo/gloo/realtime_signals.py
# See the file "LICENSE_vispy.txt" for the copyright license associated with that original example.
#
# Otherwise, modifications made herein are considered copyright under the terms described in the file "LICENSE_actichampy.txt"
import sys
sys.path.append('./pylsl') # help python find pylsl relative to this example program
from pylsl import StreamInfo, StreamInlet, resolve_stream

from vispy import gloo
from vispy import app
import numpy as np
import scipy.fftpack
import scipy.special
import math
import OpenGL.GL as gl
import time

VERT_SHADER = """
#version 120
// y coordinate of the position.
attribute float a_position ;
// row , col , and time index.
attribute vec3 a_index ;
varying vec3 v_index ;
// Size of the table.
uniform float u_nrows ;
// Number of samples per signal.
uniform float u_n ;
// Varying variables used for clipping in the fragment shader.
varying vec2 v_position ;
varying vec4 v_ab ;
void main() {
	// Compute the x coordinate from the time index.
	float x = -1 + 2 * a_index.z / ( u_n-1 ) ;
	vec2 position = vec2( x , a_position ) ;
	// Find the affine transformation for the subplots.
	vec2 a = vec2( 1. , 1. / u_nrows ) * .99 ;
	vec2 b = vec2( -1 + 2 * ( a_index.x + .5 ) , -1 + 2 * ( a_index.y + .5 ) / u_nrows ) ;
	// Apply the static subplot transformation + scaling.
	gl_Position = vec4( a * position + b , 0.0 , 1.0 ) ;
	v_index = a_index ;
	// For clipping test in the fragment shader.
	v_position = gl_Position.xy ;
	v_ab = vec4( a , b ) ;
}
"""

FRAG_SHADER = """
#version 120
// Number of samples per signal.
uniform float u_alpha ;
varying vec3 v_index ;
varying vec2 v_position ;
varying vec4 v_ab ;
void main() {
	gl_FragColor = vec4( 1. , 1. , 1. , u_alpha ) ;
	// Discard the fragments between the signals ( emulate glMultiDrawArrays ).
	if ( ( fract( v_index.x ) > 0. ) || ( fract( v_index.y ) > 0. ) )
		discard ;
	// Clipping test.
	vec2 test = abs( ( v_position.xy-v_ab.zw ) / v_ab.xy ) ;
	if ( ( test.x > 1 ) || ( test.y > 1 ) )
		discard ;
}
"""

class Canvas( app.Canvas ):
	def update_index( self ):
		self.index = np.c_[ np.repeat( np.repeat( 0 , self.nrows ) , 2 ** self.samples_per_screen_pow ) ,
			  np.repeat( np.tile( np.arange( self.nrows ) , 1 ) , 2 ** self.samples_per_screen_pow ) ,
			  np.tile( np.arange( 2 ** self.samples_per_screen_pow ) , self.nrows ) ].astype( np.float32 )
	def __init__( self , data ):
		self.data = data
		self.alpha_shift = 12
		self.samples_per_screen_pow = np.floor( np.log2( self.data.shape[ 1 ] ) )
		self.show_latest = True
		self.last_sample_to_show = self.data.shape[ 1 ]#None

		streams = resolve_stream('name', 'RandomSpehricalData')
		streamInfo = streams[0]
		# create a new inlet to read from the stream
		self.inlet = StreamInlet(streamInfo)
		# Number of cols and rows in the table.
		self.nrows = streamInfo.channel_count()

		#self.nrows = self.data.shape[ 0 ]
		self.update_index( ) #creates self.index
		self.equalize_amplitudes = False
		self.bandpass = False
		app.Canvas.__init__( self , title = 'Bla' , keys = 'interactive' )
		self.program = gloo.Program( VERT_SHADER , FRAG_SHADER )
		self.program[ 'a_position' ] = self.data.reshape( -1 , 1 )
		self.program[ 'a_index' ] = self.index
		self.program[ 'u_nrows' ] = self.nrows
		self.program[ 'u_n' ] = 2 ** self.samples_per_screen_pow
		self.program[ 'u_alpha' ] = 1 - scipy.special.expit( ( self.samples_per_screen_pow - self.alpha_shift ) )#1.0 / np.log2( self.samples_per_screen_pow ** 2 )
		gloo.set_state( clear_color = 'black' , blend = True , blend_func = ( 'src_alpha' , 'one_minus_src_alpha' ) )
		gl.glEnable( gl.GL_BLEND )
		gl.glBlendFunc( gl.GL_SRC_ALPHA , gl.GL_ONE_MINUS_SRC_ALPHA )
		gl.glEnable( gl.GL_LINE_SMOOTH )
		gl.glHint( gl.GL_LINE_SMOOTH_HINT , gl.GL_NICEST )

	def on_resize( self , event ):
		self.width , self.height = event.size
		gloo.set_viewport( 0 , 0 , self.width , self.height )
	def on_key_press( self , event ):
		if event.key.name in [' ','Space']: #toggle show_latest on / off
			self.show_latest = not self.show_latest
		elif event.key.name in ['a','A']:
			self.equalize_amplitudes = not self.equalize_amplitudes
		elif event.key.name in ['b','B']:
			self.bandpass = not self.bandpass
		elif event.key.name == 'Left':
			self.samples_per_screen_pow += 1
		elif event.key.name == 'Right':
			self.samples_per_screen_pow -= 1
		elif event.key.name == 'Up':
			self.alpha_shift += 1
		elif event.key.name == 'Down':
			self.alpha_shift -= 1
		else:
			print event.key.name
	def on_mouse_wheel( self , event ):
		if not self.show_latest:
			self.last_sample_to_show = self.last_sample_to_show + int( ( 2 ** self.samples_per_screen_pow ) * ( event.delta[ 0 ] / 100 ) )
			if self.last_sample_to_show < ( 2 ** self.samples_per_screen_pow ):
				self.last_sample_to_show = ( 2 ** self.samples_per_screen_pow )
			elif self.last_sample_to_show > self.data.shape[ 1 ]:
				self.last_sample_to_show = self.data.shape[ 1 ]

	def on_draw( self , event ):
		start = time.time( )
		"""Add some data at the end of each signal ( real-time signals )."""
		k = int( 1000 / 60 ) #emulating 1kHz new data rate ( assuming this function gets called every refresh at 60Hz refresh rate )

		sampleSet = amplitudes * np.random.randn( self.nrows , k )

		sample, timestamp = self.inlet.pull_sample(0.0)
		if sample != None:
			while len(sample) > 0:
				k = k + 1
				sample, timestamp = self.inlet.pull_sample(0.0)
				sample = np.array(sample).reshape((1,self.nrows))
				np.c_[ sampleSet, sample ]

		# original
		#self.data = np.c_[ self.data , amplitudes * np.random.randn( self.nrows , k ) ]
		self.data = np.c_[ self.data ,  np.zeros ]

		if self.samples_per_screen_pow < 2:
			self.samples_per_screen_pow = 2
		while ( 2 ** self.samples_per_screen_pow ) > self.data.shape[ 1 ]:
			self.samples_per_screen_pow -= 1
		if self.show_latest:
			self.last_sample_to_show = self.data.shape[ 1 ]
			this_data = self.data[ : , -( 2 ** self.samples_per_screen_pow ): ].copy( )
		else:
			while ( 2 ** self.samples_per_screen_pow ) > self.last_sample_to_show:
				self.samples_per_screen_pow -= 1
			this_data = self.data[ : , self.last_sample_to_show-( 2 ** self.samples_per_screen_pow ):self.last_sample_to_show ].copy( )
		from_fft = scipy.fftpack.rfft( this_data )
		# power = 2.0 / from_fft.shape[ 1 ] * np.abs( from_fft[ 0:from_fft.shape[ 1 ] / 2 ] )
		if self.bandpass:
			#apply bandpass
			W = scipy.fftpack.fftfreq( this_data.shape[ 1 ] , d = 0.001 )
			from_fft[ : , ( W < 0.1 ) ] = 0
			from_fft[ : , ( W > 30 ) ] = 0
			this_data = scipy.fftpack.irfft( from_fft )
		if self.equalize_amplitudes:
			this_data = this_data.transpose( )
			this_data -= np.amin( this_data , axis = 0 )
			this_data /= np.amax( this_data , axis = 0 )
			this_data *= 2
			this_data -= 1
			this_data = this_data.transpose( )
		else:
			this_data -= np.amin( this_data )
			this_data /= np.amax( this_data )
			this_data *= 2
			this_data -= 1
		# print time.time( )-start
		self.update_index()
		self.program[ 'a_index' ] = self.index
		self.program[ 'u_n' ] = 2 ** self.samples_per_screen_pow
		self.program[ 'u_alpha' ] = 1 - scipy.special.expit( ( self.samples_per_screen_pow - self.alpha_shift ) )#1.0 / np.log2( self.samples_per_screen_pow ** 2 )
		self.program[ 'a_position' ] = this_data.reshape( -1 , 1 ).astype( np.float32 )
		self.update()
		gloo.clear()
		self.program.draw( 'line_strip' )

if __name__ == '__main__':
	nrows = 3
	n = 10000
	amplitudes = .1 + .2 * np.random.rand( nrows , 1 ).astype( np.float32 )
	y = amplitudes * np.random.randn( nrows , n ).astype( np.float32 )
	c = Canvas( data = y )
	c.show( )
	app.run( )

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon May 25 22:33:04 2020

@author: kkrao
"""

import pandas as pd
from init import dir_data, lc_dict, color_dict
import seaborn as sns
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import mannwhitneyu, pearsonr
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
from sklearn.metrics import plot_roc_curve, roc_auc_score
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm

color_dict['All forests'] = 'k'

sns.set(style='ticks',font_scale = 0.9)

def assemble_df(trait =  "p50"):
    df = pd.read_csv(os.path.join(dir_data, "fire_collection_median_with_climate_500m_variogram.csv"))
    # df = pd.read_csv(os.path.join(dir_data, "fire_collection_median_with_climate_500m_variogram_6_jul_2020.csv"))

    dfr = pd.read_csv(os.path.join(dir_data, "fire_collection_median_extra_lfmc_vars_500m_variogram.csv"))
    dfr = dfr[['lfmc_t_1_seasonal_mean_inside','lfmc_t_1_seasonal_mean_outside', 'lfmc_t_2_inside', 'lfmc_t_2_outside']]
    df = df.join(dfr)

    dfr = pd.read_csv(os.path.join(dir_data, "fire_collection_median_fwi_500m_variogram.csv"))
    dfr = dfr[['fwi_t_4_inside','fwi_t_4_outside']]
    df = df.join(dfr)

    dfr = pd.read_csv(os.path.join(dir_data, "fires_2016_2019_with_p50_KL_isohydricity_rootdepth_500m.csv"))
    dfr = dfr[[trait]]
    df = df.join(dfr)

    df = df.loc[df.landcover.isin(lc_dict.keys())]
    df['landcover'] = df.landcover.map(lc_dict)
    return df

def plot_hist(df, trait, xlab, bins = 10):
    fig, ax = plt.subplots(figsize = (3,3))
    df[trait].hist(ax = ax, bins = bins)
    ax.set_xlabel(xlab)
    ax.set_ylabel('Frequency')

    return ax

def find_nearest(array, value, trait):
    array = np.asarray(array)
    if trait=="rootdepth":
        idx = (np.abs(np.log10(array) - np.log10(value))).argmin()
    else:
        idx = (np.abs(array - value)).argmin()
    return array[idx]

def custom_round(x, base=2):
    return int(base * round(float(x)/base))

def myround(x, base=2):
    return np.ceil(x / float(base)) * base -1
#%% just lfmc first
def auc_by_trait(clf, ndf, kind = 'occurence', trait = "p50", sequence = np.arange(-12,1, 2)):
    auc = pd.DataFrame(index = [0],columns = sequence)
    for p in auc.columns:
        sub = ndf.loc[ndf[trait].apply(lambda x: find_nearest(sequence, x, trait))==p]
        if kind == "occurence":
            X = sub.drop(['fire',trait], axis = 1)
            y = sub['fire']
            try:
                auc.loc[0,p] = roc_auc_score(y, clf.predict(X))
            except:
                auc.loc[0,p] = np.nan
        else:
            X = sub.drop(['size',trait], axis = 1)
            y = sub['size']
            try:
                auc.loc[0,p] = roc_auc_score(y, clf.predict(X))
            except:
                auc.loc[0,p] = np.nan
    return auc

def calc_auc_size(dfsub, clf, trait = "p50",sequence = np.arange(-12,1, 2)):
    ndf = dfsub.copy()
    ndf = ndf.sample(frac=1).reset_index(drop=True)
    ndf.dropna(inplace = True)
    # print(ndf.columns)
    X = ndf.drop(['size',trait], axis = 1)
    y = ndf['size']
    # print(y.mean())
    # try:
    clf.fit(X, y)
        # rfc_disp = plot_roc_curve(clf, X, y, ax=ax,label = lc,color = color_dict[lc])
    auc = auc_by_trait(clf, ndf, kind = "size", trait = trait,sequence = np.arange(-12,1, 2))
        # print(roc_auc_score(y, clf.predict(X)))
    # except:
        # print("Could not fit RF")
    # print(auc)
    return auc

def calc_auc_occurence(dfsub,  clf, trait = "p50",sequence =np.arange(-12,1, 2)):
    df = dfsub.copy()

    ndf = pd.DataFrame()
    for var in ['outside','inside']:
        cols = [col for col in df.columns if var in col]+[trait]
        # cols.remove('lfmc_t_1_%s'%var)
        data = df[cols].copy()
        new_cols = [col.split('_')[0] for col in data.columns]
        data.columns = (new_cols)
        data['fire'] = int(var=='inside')
        ndf = pd.concat([ndf, data], axis = 0).reset_index(drop=True)

    ndf = ndf.sample(frac=1).reset_index(drop=True)
    ndf.dropna(inplace = True)

    X = ndf.drop(['fire',trait], axis = 1)
    y = ndf['fire']

    # try:
    clf.fit(X, y)
        # rfc_disp = plot_roc_curve(clf, X, y, ax=ax,label = lc,color = color_dict[lc])
    auc = auc_by_trait(clf, ndf, kind = "occurence", trait = trait,sequence = sequence)
    # except:
        # print("Could not fit RF")

    return auc

def ensemble_auc(dfsub, clf, iters = 100, kind = "occurence", trait = "p50",sequence =np.arange(-12,1, 2)):
    clf.random_state = 0
    # dummy = calc_auc_occurence(dfsub, category_dict, clf)
    # aucs = np.expand_dims(dummy.values, axis = 2)
    for itr in range(1, iters):
        clf.random_state = itr
        if itr ==1:
            if kind == 'occurence':
                auc = calc_auc_occurence(dfsub, clf, trait = trait,sequence = sequence)
            else:
                auc = calc_auc_size(dfsub, clf, trait = trait,sequence = sequence)
        else:
            if kind == 'occurence':
                auc = auc.append(calc_auc_occurence(dfsub, clf, trait = trait,sequence = sequence)).reset_index(drop=True)
            else:
                auc = auc.append(calc_auc_size(dfsub, clf, trait = trait,sequence = sequence)).reset_index(drop=True)
        # aucs = np.append(aucs,auc, axis = 2)
    # print("aucs ready")
    # dummy.loc[:,:] = np.nanmean(aucs.astype(float), axis = 2)
    # mean = dummy.copy()
    # dummy.loc[:,:] = np.nanstd(aucs.astype(float), axis = 2)
    # sd = dummy.copy()

    return auc

def calc_auc_diff(dfs, replace_by_random = False, kind = "occurence",trait = "p50",sequence =np.arange(-12,1, 2), iters = 100):
    df = dfs.copy()
    # allVars = pd.DataFrame(index = [0],columns = category_dict.keys())
    # onlyClimate = allVars.copy()
    cols = [col for col in df.columns if 'lfmc' in col]+[trait,'area']
    # cols = ['landcover']
    cols+=[col for col in df.columns if 'erc' in col]
    cols+=[col for col in df.columns if 'ppt' in col]
    cols+=[col for col in df.columns if 'vpd' in col]
    cols+=[col for col in df.columns if 'fwi' in col]

    df = df[cols]
    df['lfmc_t_1_inside_anomaly'] = df['lfmc_t_1_inside'] - df['lfmc_t_1_seasonal_mean_inside']
    df['lfmc_t_1_outside_anomaly'] = df['lfmc_t_1_outside'] - df['lfmc_t_1_seasonal_mean_outside']

    df['size'] = 0
    df.loc[df.area>4,'size'] = 1

    df.drop('area',axis = 1, inplace = True)

    df.drop(['lfmc_t_1_seasonal_mean_inside','lfmc_t_1_seasonal_mean_outside'],axis = 1, inplace = True)

    if kind == 'size':
        remove_outside = [col for col in df.columns if "outside" in col]
        df.drop(remove_outside, axis = 1,inplace = True)
    ###testing with random numbers instead of LFMC
    # df.loc[:,remove_lfmc] = np.zeros(shape = df.loc[:,remove_lfmc].shape)
    # clf = RandomForestClassifier(max_depth=15, min_samples_leaf = 5, random_state=0, oob_score = True,n_estimators = 50)

    if kind=='occurence':
        clf = RandomForestClassifier(max_depth=6, random_state=0, oob_score = True,n_estimators = 40)
    else:
        clf = RandomForestClassifier(max_depth=10, min_samples_leaf = 1, random_state=0, oob_score = True,n_estimators = 20)
    allVars = ensemble_auc(df, clf, kind = kind, trait = trait, sequence = sequence, iters = iters)

    # allVars = calc_auc(df, size_dict, clf)

    remove_lfmc = [col for col in df.columns if 'lfmc' in col]
    if replace_by_random:
        ###testing with random numbers instead of LFMC
        df.loc[:,remove_lfmc] = np.ones(shape = df.loc[:,remove_lfmc].shape)
        onlyClimate, s2 = ensemble_auc(df, clf, kind = kind, trait = trait, sequence = sequence)
    else:
        onlyClimate = ensemble_auc(df.drop(remove_lfmc, axis = 1), clf, kind = kind, trait = trait, sequence = sequence)

    diff = (allVars - onlyClimate).copy().astype(float).round(3)
    onlyClimate.index.name = "only climate"
    diff.index.name = "difference, mean"
    allVars.index.name = "all variables"

    # sd = (s1.pow(2)+s2.pow(2)).pow(0.5).astype(float).round(3)
    # sd.index.name = "difference, sd"
    # print(onlyClimate.astype(float).round(2))
    # print(allVars.astype(float).round(2))
    # print(diff.astype(float).round(2))
    # print(sd.astype(float).round(2))

    return allVars, onlyClimate

def plot_importance(allVars, onlyClimate, xlab = 'P50 (MPA)', xlim = [1,-15], xticks = [0,-5,-10,-15], lc = None, n = None):
    y = (allVars - onlyClimate)*100
    mean = y.mean()
    sd = y.std()

    fig, ax= plt.subplots(figsize = (3,3))
    ax.plot(mean,'o-',color = color_dict[lc], markeredgecolor = "grey")
    ax.errorbar(x = mean.index, y = mean, yerr = sd, fmt = 'o', color = "grey", capsize = 2, zorder = -1)
    ax.scatter(x = mean.index, y = mean,marker = 'o', edgecolor = "grey",color = color_dict[lc])

    ax.set_ylabel("LFMC Value (%)")
    ax.set_xlabel(xlab)
    ax.set_xlim(*xlim)
    ax.set_xticks(xticks)
    # ax1.set_xlim(0.5,1)
    ax.set_title("%s, N = %d"%(lc, n))

    return ax

#%%
# ITERS = 100
# trait = "isohydricity"
# df = assemble_df(trait)
# # df = df.loc[df.isohydricity>=0]
# plot_hist(df, trait = trait,xlab = "$\sigma$")
# corrs = pd.DataFrame(index = df.landcover.unique(), columns = ["r", "p"])
# for lc in df.landcover.unique():
#     dfsub = df.loc[df.landcover==lc].copy()
#     allVars, onlyClimate = calc_auc_diff(dfsub, replace_by_random = False,
#                                           kind = 'occurence',
#                                           trait = trait,
#                                           iters = ITERS,
#                                           sequence = np.linspace(-0.,1.0,5))
#     plot_importance(allVars, onlyClimate,
#                     xlab = "$\sigma$",
#                     xlim = [-0.25,1.25],
#                     xticks = np.linspace(-0.,1.0,5),
#                     lc = lc, n = dfsub.dropna().shape[0])
#     corrs.loc[lc] = pearsonr(np.repeat(allVars.columns.values, allVars.shape[0]),
#                              (allVars - onlyClimate).values.flatten(order = "F").astype(float))
# print(corrs)

#%%
# ITERS = 100
# trait = "p50"
# df = assemble_df(trait)
# # df = df.loc[df.isohydricity>=0]
# plot_hist(df, trait = trait,xlab = "P50 (Mpa)")
# corrs = pd.DataFrame(index = df.landcover.unique(), columns = ["r", "p"])
# for lc in ['Closed broadleaf\ndeciduous', 'Closed needleleaf\nevergreen','Mixed forest']:
#     dfsub = df.loc[df.landcover==lc].copy()
#     allVars, onlyClimate = calc_auc_diff(dfsub, replace_by_random = False,
#                                           kind = 'occurence',
#                                           trait = trait,
#                                           iters = ITERS,
#                                           sequence = [-13,-9,-7,-5,-3,-1])
#     plot_importance(allVars, onlyClimate,
#                     xlab = "P50 (Mpa)",
#                     xlim = [-14,0],
#                     xticks = np.linspace(-13,-1,7),
#                     lc = lc, n = dfsub.dropna().shape[0])
#     corrs.loc[lc] = pearsonr(np.repeat(allVars.columns.values, allVars.shape[0]),
#                               (allVars - onlyClimate).values.flatten(order = "F").astype(float))
# print(corrs)

# lc = "All forests"
# dfsub = df.loc[df.landcover.isin(['Closed broadleaf\ndeciduous', 'Closed needleleaf\nevergreen','Mixed forest'])].copy()
# allVars, onlyClimate = calc_auc_diff(dfsub, replace_by_random = False,
#                                       kind = 'occurence',
#                                       trait = trait,
#                                       iters = ITERS,
#                                       sequence = [-13,-9,-7,-5,-3,-1])
# plot_importance(allVars, onlyClimate,
#                 xlab = "P50 (Mpa)",
#                 xlim = [-14,0],
#                 xticks = np.linspace(-13,-1,7),
#                 lc = lc)

#%%
ITERS = 100
trait = "root_depth" #need to remove underscore because _inside _outside vars are filtered based on _
df = assemble_df(trait)
df = df.rename(columns = {'root_depth':'rootdepth'})
trait = "rootdepth"
ax = plot_hist(df, trait = trait, xlab = "Rooting depth (m)", bins = 100)
ax.set_xscale('log')
corrs = pd.DataFrame(index = df.landcover.unique(), columns = ["r", "p"])
for lc in df.landcover.unique():
    dfsub = df.loc[df.landcover==lc].copy()
    allVars, onlyClimate = calc_auc_diff(dfsub, replace_by_random = False,
                                          kind = 'occurence',
                                          trait = trait,
                                          iters = ITERS,
                                          sequence =np.array([0.01,0.1,1,10]))
    ax = plot_importance(allVars, onlyClimate,
                    xlab = "Rooting depth(m)",
                      xlim = [0.0005,20],
                    xticks = [0.001,0.01,0.1,1,10],
                    lc = lc, n = dfsub.dropna().shape[0])
    ax.set_xscale('log')

    corrs.loc[lc] = pearsonr(np.log10(np.repeat(allVars.columns.values, allVars.shape[0])),
                              (allVars - onlyClimate).values.flatten(order = "F").astype(float))
print(corrs)

reason: import numpy,from scipy,import statsmodels


text:
from skimage.util import pad
import numpy as np

from dataset.loader.DatasetLoader import DatasetLoader

class Pad(DatasetLoader):
    def __init__(self, baseDatasetLoader, padding, imageDimensions, labelDimensions):
        self.__baseDatasetLoader = baseDatasetLoader
        self.__padding = padding
        self.__imageDimensions = imageDimensions
        self.__labelDimensions = labelDimensions

    def loadData(self) -> ((np.ndarray, np.ndarray), (np.ndarray, np.ndarray), (np.ndarray, np.ndarray)):
        (xTrain, yTrain), (xVal, yVal), (xTest, yTest) = self.__baseDatasetLoader.loadData()
        xTrainPadded = self.__pad(xTrain)
        xValPadded = self.__pad(xVal)
        xTestPadded = self.__pad(xTrain)
        return (xTrainPadded, yTrain), (xValPadded, yVal), (xTestPadded, yTest)

    def __pad(self, images):
        return np.array([pad(image, self.__padding, 'constant') for image in images])

    def dataPointShape(self):
        return self.__imageDimensions, self.__labelDimensions

max_stars_repo_name: dhruvtapasvi/implementation

max_issues_repo_name: dhruvtapasvi/implementation

max_forks_repo_name: dhruvtapasvi/implementation