text:
import neptune
from config import EnvConfig
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

comet_cfg = EnvConfig()
session = neptune.Session(api_token=comet_cfg.neptune_token)
project = session.get_project(project_qualified_name=comet_cfg.neptune_project_name)

experiments = project.get_experiments()[171:]
print(experiments)

params = []
lyap_inter = []
lyap_intra = []
neg_intra = []

for exp in experiments:

    properties = exp.get_properties()

    model_update = properties['target_model_update']
    mem_len = properties['memory_limit']
    alpha = properties['learning_rate']

    try:
        # print(exp.id)
        lyap_intra = exp.get_numeric_channels_values('lyap_exp_intra_ins_0','lyap_exp_intra_ins_1','lyap_exp_intra_ins_2').to_numpy()[:,1:]
        neg_intra.append(sum(sum(lyap_intra<0))/600)
    # print(lyap)
        params.append([model_update, mem_len, alpha])
    except:
        print('except '+exp.id)

params = np.array(params).astype(np.float)
lyap_intra = np.array(lyap_intra).astype(np.float)

plt.figure()
plt.subplot(1,3,1)
plt.scatter(params[:,0], neg_intra)
plt.xlabel('target_model_update')
plt.ylabel('Percentage of episodes w/ negative lyap exp')
plt.subplot(1,3,2)
plt.scatter(params[:,1], neg_intra)
plt.xlabel('memory_limit')
plt.subplot(1,3,3)
plt.scatter(params[:,2], neg_intra)
plt.xlabel('alpha')

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
img = ax.scatter(params[:,0],params[:,1], params[:,2], c=neg_intra, cmap=plt.hot())
ax.set_xlabel('Model_update'); ax.set_ylabel('Memory_limit'); ax.set_zlabel('Learning_rate');
fig.colorbar(img)
plt.title('3d plot of percentage of negative exponent episodes')

plt.figure()
plt.scatter(range(len(neg_intra)), neg_intra)
plt.show()

print('done')


text:
"""Unit tests for loader module."""

# -----------------------------------------------------------------------------
# Imports
# -----------------------------------------------------------------------------
import os
from collections import Counter

import numpy as np
import numpy.random as rnd
import pandas as pd
import shutil
from nose.tools import with_setup

from mock_data import setup as setup_klusters
from mock_data import (teardown, TEST_FOLDER, nspikes, nclusters, nsamples,
    nchannels, fetdim)
from kwiklib.dataio import (KwikLoader, Experiment, klusters_to_kwik,
    check_dtype, check_shape, get_array, select, get_indices)

# -----------------------------------------------------------------------------
# Fixtures
# -----------------------------------------------------------------------------
def setup():
    setup_klusters()
    klusters_to_kwik(filename='test', dir=TEST_FOLDER)

# -----------------------------------------------------------------------------
# Tests
# -----------------------------------------------------------------------------
def test_kwik_loader_1():

    # Open the mock data.
    dir = TEST_FOLDER
    xmlfile = os.path.join(dir, 'test.xml')
    l = KwikLoader(filename=xmlfile)

    # Get full data sets.
    features = l.get_features()
    # features_some = l.get_some_features()
    masks = l.get_masks()
    waveforms = l.get_waveforms()
    clusters = l.get_clusters()
    spiketimes = l.get_spiketimes()
    nclusters = len(Counter(clusters))

    # probe = l.get_probe()
    cluster_colors = l.get_cluster_colors()
    cluster_groups = l.get_cluster_groups()
    group_colors = l.get_group_colors()
    group_names = l.get_group_names()
    cluster_sizes = l.get_cluster_sizes()

    # Check the shape of the data sets.
    # ---------------------------------
    assert check_shape(features, (nspikes, nchannels * fetdim + 1))
    # assert features_some.shape[1] == nchannels * fetdim + 1
    assert check_shape(masks, (nspikes, nchannels * fetdim + 1))
    assert check_shape(waveforms, (nspikes, nsamples, nchannels))
    assert check_shape(clusters, (nspikes,))
    assert check_shape(spiketimes, (nspikes,))

    # assert check_shape(probe, (nchannels, 2))
    assert check_shape(cluster_colors, (nclusters,))
    assert check_shape(cluster_groups, (nclusters,))
    assert check_shape(group_colors, (4,))
    assert check_shape(group_names, (4,))
    assert check_shape(cluster_sizes, (nclusters,))

    # Check the data type of the data sets.
    # -------------------------------------
    assert check_dtype(features, np.float32)
    assert check_dtype(masks, np.float32)
    # HACK: Panel has no dtype(s) attribute
    # assert check_dtype(waveforms, np.float32)
    assert check_dtype(clusters, np.int32)
    assert check_dtype(spiketimes, np.float64)

    # assert check_dtype(probe, np.float32)
    assert check_dtype(cluster_colors, np.int32)
    assert check_dtype(cluster_groups, np.int32)
    assert check_dtype(group_colors, np.int32)
    assert check_dtype(group_names, object)
    assert check_dtype(cluster_sizes, np.int32)

    l.close()

def test_kwik_loader_control():
    # Open the mock data.
    dir = TEST_FOLDER
    xmlfile = os.path.join(dir, 'test.xml')
    l = KwikLoader(filename=xmlfile)

    # Take all spikes in cluster 3.
    spikes = get_indices(l.get_clusters(clusters=3))

    # Put them in cluster 4.
    l.set_cluster(spikes, 4)
    spikes_new = get_indices(l.get_clusters(clusters=4))

    # Ensure all spikes in old cluster 3 are now in cluster 4.
    assert np.all(np.in1d(spikes, spikes_new))

    # Change cluster groups.
    clusters = [2, 3, 4]
    group = 0
    l.set_cluster_groups(clusters, group)
    groups = l.get_cluster_groups(clusters)
    assert np.all(groups == group)

    # Change cluster colors.
    clusters = [2, 3, 4]
    color = 12
    l.set_cluster_colors(clusters, color)
    colors = l.get_cluster_colors(clusters)
    assert np.all(colors == color)

    # Change group name.
    group = 0
    name = l.get_group_names(group)
    name_new = 'Noise new'
    assert name == 'Noise'
    l.set_group_names(group, name_new)
    assert l.get_group_names(group) == name_new

    # Change group color.
    groups = [1, 2]
    colors = l.get_group_colors(groups)
    color_new = 10
    l.set_group_colors(groups, color_new)
    assert np.all(l.get_group_colors(groups) == color_new)

    # Add cluster and group.
    spikes = get_indices(l.get_clusters(clusters=3))[:10]
    # Create new group 100.
    l.add_group(100, 'New group', 10)
    # Create new cluster 10000 and put it in group 100.
    l.add_cluster(10000, 100, 10)
    # Put some spikes in the new cluster.
    l.set_cluster(spikes, 10000)
    clusters = l.get_clusters(spikes=spikes)
    assert np.all(clusters == 10000)
    groups = l.get_cluster_groups(10000)
    assert groups == 100
    l.set_cluster(spikes, 2)

    # Remove the new cluster and group.
    l.remove_cluster(10000)
    l.remove_group(100)
    assert np.all(~np.in1d(10000, l.get_clusters()))
    assert np.all(~np.in1d(100, l.get_cluster_groups()))

    l.close()

@with_setup(setup)
def test_kwik_save():
    """WARNING: this test should occur at the end of the module since it
    changes the mock data sets."""
    # Open the mock data.
    dir = TEST_FOLDER
    xmlfile = os.path.join(dir, 'test.xml')
    l = KwikLoader(filename=xmlfile)

    clusters = l.get_clusters()
    cluster_colors = l.get_cluster_colors()
    cluster_groups = l.get_cluster_groups()
    group_colors = l.get_group_colors()
    group_names = l.get_group_names()

    # Set clusters.
    indices = get_indices(clusters)
    l.set_cluster(indices[::2], 2)
    l.set_cluster(indices[1::2], 3)

    # Set cluster info.
    cluster_indices = l.get_clusters_unique()
    l.set_cluster_colors(cluster_indices[::2], 10)
    l.set_cluster_colors(cluster_indices[1::2], 20)
    l.set_cluster_groups(cluster_indices[::2], 1)
    l.set_cluster_groups(cluster_indices[1::2], 0)

    # Save.
    l.remove_empty_clusters()
    l.save()

    clusters = l.get_clusters()
    cluster_colors = l.get_cluster_colors()
    cluster_groups = l.get_cluster_groups()
    group_colors = l.get_group_colors()
    group_names = l.get_group_names()

    assert np.all(clusters[::2] == 2)
    assert np.all(clusters[1::2] == 3)

    assert np.all(cluster_colors[::2] == 10)
    assert np.all(cluster_colors[1::2] == 20)

    print cluster_groups

    assert np.all(cluster_groups[::2] == 1)
    assert np.all(cluster_groups[1::2] == 0)

    l.close()

ext: py

max_stars_repo_name: fiath/test

max_issues_repo_name: fiath/test

max_forks_repo_name: fiath/test

reason: import numpy


text:
""" Basic benchmarking of some Non-Cartesian MRI reconstruction cases. """
from __future__ import print_function, division

import numpy as np

from matplotlib import pyplot as plt

from mrrt.mri.operators.tests.test_mri_noncartesian import _test_mri_multi

def _split_timings(bench_tuple, print_timings=False):
    t_adj = [(k, v["MRI: adjoint"]) for k, v in bench_tuple.items()]
    adjoint_times = sorted(t_adj, key=lambda x: x[1], reverse=True)

    t_for = [(k, v["MRI: forward"]) for k, v in bench_tuple.items()]
    forward_times = sorted(t_for, key=lambda x: x[1], reverse=True)

    t_create = [(k, v["MRI: object creation"]) for k, v in bench_tuple.items()]
    creation_times = sorted(t_create, key=lambda x: x[1], reverse=True)

    t_norm = [(k, v["MRI: norm"]) for k, v in bench_tuple.items()]
    norm_times = sorted(t_norm, key=lambda x: x[1], reverse=True)

    if print_timings:
        print("creation times")
        for k, v in creation_times[::-1]:
            print("\t{}: {}".format(k, v))
        print("forward times")
        for k, v in forward_times[::-1]:
            print("\t{}: {}".format(k, v))
        print("adjoint times")
        for k, v in adjoint_times[::-1]:
            print("\t{}: {}".format(k, v))
        print("norm times")
        for k, v in norm_times[::-1]:
            print("\t{}: {}".format(k, v))
    return creation_times, forward_times, adjoint_times, norm_times

def _autolabel(ax, rects, fmt="%0.3g"):
    """
    Attach a text label to the right of each bar displaying its width
    """
    # use mean_width to determine text spacing so that it will be the same
    # for all bars
    mean_width = np.mean([rect.get_width() for rect in rects])
    max_width = 0
    for rect in rects:
        width = rect.get_width()
        t = ax.text(
            width + 0.02 * mean_width,
            rect.get_y() + rect.get_height() / 2.0,
            fmt % width,
            va="center",
            ha="left",
        )
        if width > max_width:
            t_max = t
            max_width = width

    # determine size of the text as in:
    #   http://stackoverflow.com/questions/5320205/matplotlib-text-dimensions
    bb = t_max.get_window_extent(renderer=ax.get_figure().canvas.get_renderer())
    bb_ax = ax.get_window_extent()
    stretch_factor = (bb_ax.width + bb.width) / bb_ax.width
    # stretch xlim a bit to accomodate the labels
    ax.set_xlim((ax.get_xlim()[0], stretch_factor * ax.get_xlim()[1]))

def _plot_timings(
    labels_times_tuple, title=None, xlabel=None, xlim=None, ax=None, fmt="%0.3g"
):
    if ax is None:
        fig, ax = plt.subplots(1, 1, figsize=[12, 6])
    nlabel = len(labels_times_tuple)
    labels = [a[0] for a in labels_times_tuple]
    labels = [l.replace("single", "s") for l in labels]
    labels = [l.replace("double", "d") for l in labels]
    labels = [l.replace("complex", "C") for l in labels]
    labels = [l.replace("real", "R") for l in labels]
    rects = ax.bar(
        np.zeros(nlabel),
        height=0.8,
        bottom=np.arange(nlabel),
        width=[a[1] for a in labels_times_tuple],
        tick_label=labels,
        orientation="horizontal",
        align="center",
    )
    _autolabel(ax, rects, fmt=fmt)
    if title is not None:
        ax.set_title(title)
    if xlabel is not None:
        ax.set_xlabel(xlabel)
    # ax.axis('tight')
    if xlim is not None:
        ax.set_xlim(xlim)
    #        mean_width = np.mean([r.get_width() for r in rects])
    #        for rect in rects:
    #            width = rect.get_width()
    #            if width < xlim[-1]:
    #                t = ax.text(xlim[-1] + 0.02*mean_width,
    #                            rect.get_y() + rect.get_height()/2.,
    #                            '%0.3g' % width,
    #                            va='center', ha='left')

    # center the text of any multi-line labels
    ax.yaxis.set_ticklabels(ax.yaxis.get_ticklabels(), multialignment="center")
    return ax

if False:
    # BART CPU without Toeplitz 2873 ms
    # BART CPU with Toeplitz 2154 ms
    # BART GPU with Toeplitz 2027 ms
    # BART GPU without Toeplitz 2027 msG
    # ('BART (CPU)', 2154),

    from matplotlib import rcParams

    rcParams["font.size"] = 14

    labels_times_tuple = (
        ("Michigan\nIRT", 5861),
        # fftn/ifftn with FFTW_MEAURE planner effort
        ("PythonCPU\n(16 cores)", 2670),
        ("BART (GPU)", 2027),
        ("gpuNUFFT", 353),
        ("PythonGPU\n(No Precomp)", 58.5),
        ("PythonGPU\n(Full_Precomp)", 56.2),
    )
    fig1, ax1 = plt.subplots(1, 1)
    _plot_timings(
        labels_times_tuple[::-1],
        title="96x96x96 with 16 coils",
        xlabel="$A^{H}A$ execution time (ms)",
        ax=ax1,
        fmt="%d",
    )
    plt.tight_layout()

    labels_times_tuple = (
        ("Single Coil", 122),
        ("Single Coil\n(no Precomp)", 138),
        ("16 Coils", 588),
        ("16 Coils\n(no Precomp)", 727),
        ("16 Coils +\nFieldmap\n(6 time segments)", 3090),
        ("16 Coils +\nFieldmap\n(6 time segments)\n(No Precomp)", 4080),
    )
    fig2, ax2 = plt.subplots(1, 1)
    _plot_timings(
        labels_times_tuple[::-1],
        title="PythonGPU\n(256x256x256 voxels)",
        xlabel="$A^{H}A$ execution time (ms)",
        ax=ax2,
        fmt="%d",
    )
    plt.tight_layout()

    labels_times_tuple = (
        ("SimpleITK\nCPU\n(16 cores))", 706),
        ("PythonGPU", 57.7),
        ("PythonGPU\n(warp field\npreloaded)", 37.8),
    )
    fig3, ax3 = plt.subplots(1, 1)
    _plot_timings(
        labels_times_tuple[::-1],
        title=("Warp with cubic interpolation" "\n(200x148x96 complex array)"),
        xlabel="execution time (ms)",
        ax=ax3,
        fmt="%d",
    )
    plt.tight_layout()

    out_dir = "/home/lee8rx/Dropbox/Grants/SSE/LaTex_NIH/Figures/"
    from os.path import join as pjoin

    fig1.savefig(pjoin(out_dir, "comparison2.png"))
    fig1.savefig(pjoin(out_dir, "comparison2.eps"))
    fig1.savefig(pjoin(out_dir, "comparison2.pdf"))
    fig2.savefig(pjoin(out_dir, "gpu_cases.png"))
    fig2.savefig(pjoin(out_dir, "gpu_cases.eps"))
    fig2.savefig(pjoin(out_dir, "gpu_cases.pdf"))
    fig3.savefig(pjoin(out_dir, "warp.png"))
    fig3.savefig(pjoin(out_dir, "warp.eps"))
    fig3.savefig(pjoin(out_dir, "warp.pdf"))
    # 96x96x96x16coils with 400 shots instead of 150
    # PythonGPU: 90.1 ms
    # BART CPU with Toeplitz: 2205 ms
    # BART GPU with Toeplitz: 2223 ms

    # for data of size 128x128x128
    ("ShearLab3D (dec)", 14.626),
    # ('ShearLab3D (thresh)', 5.264),
    ("ShearLab3D (rec)", 16.303),
    ("UDCT (dec)", 3.37),
    ("UDCT (rec)", 1.71),
    ("DAST (Matlab) (dec)", 1.298),  # (double precision)
    # ('DAST (Matlab) (thresh)', 8.145),
    ("DAST (Matlab) (rec)", 1.577),
    ("DAST (Python) (dec)", 0.845),  # (single precision)
    # ('DAST (Python) (dec)', 0.786),  # w/ batchFFT
    # ('DAST (Python) (thresh)', 1.668),
    ("DAST (Python) (rec)", 0.865),
    ("Surfacelet (dec, R=6.4, mex)", 1.538)
    ("Surfacelet (dec, R=6.4, pure Matlab)", 2.137)
    ("Surfacelet (rec, R=6.4, mex)", 1.334)
    ("Surfacelet (rec, R=6.4, pure Matlab)", 2.282)
    # ('DAST (Python) (rec)', 0.770),   # w/ batchFFT

def _plot_all_timings(
    creation_times,
    forward_times,
    adjoint_times,
    norm_times=None,
    title_substring="",
    uniform_xscale=False,
    plot_creation_times=True,
):
    nsub = 2
    if plot_creation_times:
        nsub += 1
    if norm_times is not None:
        nsub += 1
    if nsub == 4:
        fig, axes = plt.subplots(2, 2, figsize=[16, 16])
        axes = axes.ravel()
    else:
        fig, axes = plt.subplots(nsub, 1, figsize=[16, 12])
    max_time_create = np.max([t[1] for t in creation_times])
    max_time_forward = np.max([t[1] for t in forward_times])
    max_time_adjoint = np.max([t[1] for t in adjoint_times])
    max_time = np.maximum(max_time_forward, max_time_adjoint)
    if norm_times is not None:
        max_time_norm = np.max([t[1] for t in norm_times])
        max_time = np.maximum(max_time, max_time_norm)
    if uniform_xscale:
        max_time = np.maximum(max_time, max_time_create)
        xlim1 = [0, max_time]
        xlim = xlim1
    else:
        xlim1 = None
        xlim = None  # [0, max_time]

    print([0, max_time])
    n = 0
    if plot_creation_times:
        _plot_timings(
            creation_times,
            title="creation: " + title_substring,
            xlabel="time (s)",
            xlim=xlim1,
            ax=axes[n],
        )
        n += 1
    _plot_timings(
        forward_times,
        title="forward: " + title_substring,
        xlabel="time (s)",
        xlim=xlim,
        ax=axes[n],
    )
    n += 1
    _plot_timings(
        adjoint_times,
        title="adjoint: " + title_substring,
        xlabel="time (s)",
        xlim=xlim,
        ax=axes[n],
    )
    n += 1
    if norm_times is not None:
        _plot_timings(
            norm_times,
            title="norm: " + title_substring,
            xlabel="time (s)",
            xlim=xlim,
            ax=axes[n],
        )
        n += 1
    if uniform_xscale and plot_creation_times:
        axes[0].set_xticks([])
        axes[1].set_xticks([])
        axes[0].set_xlabel("")
        axes[1].set_xlabel("")
    if not plot_creation_times:
        axes[0].set_xticks([])
        axes[0].set_xlabel("")

    return axes

def _save_bench(save_dir, fname_base, timing_tuples, axes=None):
    import pickle
    from os.path import join as pjoin

    if axes is not None:
        axes[0].get_figure().savefig(pjoin(save_dir, fname_base + ".png"))
    with open(pjoin(save_dir, fname_base + ".p"), "wb") as f:
        pickle.dump(timing_tuples, f)

def _old_bar_graph(bench_tuple):
    t = bench_tuple
    fields = ["MRI: adjoint", "MRI: forward", "MRI: object creation"]
    for field in fields:
        sorted_keys = sorted(t.keys())
        plt.figure()
        y = [t[key][field] for key in sorted_keys]
        x = np.arange(len(y))
        width = 0.8
        plt.bar(
            x,
            y,
            tick_label=[s.replace(",", "\n") for s in sorted_keys],
            width=width,
        )
        plt.xticks(plt.xticks()[0] + width / 2)
        plt.title(field, fontdict=dict(fontsize=24))

# navg_time=10; plot_timings=True; include_GPU=True; include_CPU=True; print_timings=True; phasings=['real', 'complex']

def bench_mri_2d_nocoils_nofieldmap(
    navg_time=(4, 16),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=2,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=1,
        fieldmap_segments=None,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=navg_time,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="2D (no coils, nofieldmap)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_2d_001coil_001fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_2d_16coils_nofieldmap(
    navg_time=(4, 16),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=2,
        N0=128,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=16,
        fieldmap_segments=None,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=navg_time,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="2D (16 coils, nofieldmap)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_2d_016coil_001fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_2d_16coils_fieldmap(
    navg_time=(2, 8),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=2,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=16,
        fieldmap_segments=6,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=4,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="2D (16 coils, 6 fieldmap seg)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_2d_016coil_006fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_2d_16coils_fieldmap_multispectral(
    navg_time=(1, 4),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
    nspectra=2,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]
    # spectral_offsets = [0, 434.316]
    spectral_offsets = np.asarray(
        [0, 434.316, 332.124, 485.412, -76.644, 49.8186, 247.8156]
    )
    spectral_offsets = spectral_offsets[:nspectra]

    t = _test_mri_multi(
        ndim=2,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=16,
        fieldmap_segments=6,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=4,
        spectral_offsets=spectral_offsets,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        title_string = "2D (16 coils, 6 fieldmap seg, %d spectral)" % nspectra
        axes = _plot_all_timings(
            t_create, t_forward, t_adj, t_norm, title_substring=title_string
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_2d_016coil_006fmapseg_%dspectra" % nspectra
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_3d_nocoils_nofieldmap(
    navg_time=(2, 8),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=3,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=1,
        fieldmap_segments=None,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=navg_time,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="3D (no coils, nofieldmap)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_3d_001coil_001fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_3d_nocoils_fieldmap(
    navg_time=(2, 8),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=3,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=1,
        fieldmap_segments=6,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=navg_time,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="3D (no coils, 6 fieldmap seg)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_3d_001coil_006fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_3d_16coils_nofieldmap(
    navg_time=(1, 4),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=3,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=16,
        fieldmap_segments=None,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=navg_time,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="3D (no coils, nofieldmap)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_3d_016coil_001fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

def bench_mri_3d_16coils_fieldmap(
    navg_time=(1, 4),
    plot_timings=False,
    include_GPU=False,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir=None,
):

    recon_cases = []
    if include_GPU:
        recon_cases += ["GPU,Sp", "GPU,Tab"]
        precisions = ["single"]
    else:
        precisions = ["single", "double"]
    if include_CPU:
        recon_cases += ["CPU,Tab", "CPU,Sp"]

    t = _test_mri_multi(
        ndim=3,
        N0=64,
        grid_os_factor=1.5,
        J0=4,
        Ld=4096,
        n_coils=16,
        fieldmap_segments=6,
        precisions=precisions,
        phasings=phasings,
        recon_cases=recon_cases,
        compare_to_exact=False,
        navg_time=navg_time,
    )

    t_create, t_forward, t_adj, t_norm = _split_timings(t, print_timings)
    if plot_timings:
        axes = _plot_all_timings(
            t_create,
            t_forward,
            t_adj,
            t_norm,
            title_substring="3D (16 coils, 6 fieldmap seg)",
        )
    else:
        axes = None

    if save_dir is not None:
        fname_base = "bench_3d_016coil_006fmapseg"
        _save_bench(
            save_dir, fname_base, [t_create, t_forward, t_adj, t_norm], axes
        )

# plot_timings=True; include_GPU=True; include_CPU=True; print_timings=True; phasings=['real', 'complex']; save_dir='/tmp'

def run_all_bench(
    plot_timings=True,
    include_GPU=True,
    include_CPU=True,
    print_timings=True,
    phasings=["real", "complex"],
    save_dir="/tmp",
):
    common_kwargs = dict(
        plot_timings=plot_timings,
        include_GPU=include_GPU,
        include_CPU=include_CPU,
        print_timings=print_timings,
        phasings=phasings,
        save_dir=save_dir,
    )
    if include_GPU:
        import cupy
        from cupy.cuda import device

        device.get_cusparse_handle()
        cupy.fft.cache.enable()

    # 2d cases
    bench_mri_2d_nocoils_nofieldmap(navg_time=(4, 20), **common_kwargs)
    bench_mri_2d_16coils_nofieldmap(navg_time=(2, 8), **common_kwargs)
    bench_mri_2d_16coils_fieldmap(navg_time=(2, 8), **common_kwargs)
    bench_mri_2d_16coils_fieldmap_multispectral(
        navg_time=4, nspectra=2, **common_kwargs
    )

    # 3d cases
    bench_mri_3d_nocoils_nofieldmap(navg_time=(2, 4), **common_kwargs)
    bench_mri_3d_nocoils_fieldmap(navg_time=(2, 4), **common_kwargs)
    bench_mri_3d_16coils_nofieldmap(navg_time=(1, 4), **common_kwargs)
    bench_mri_3d_16coils_fieldmap(navg_time=(1, 4), **common_kwargs)

reason: import numpy,import cupy,from cupy


text:
import cv2
import copy
import math
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import generic_filter, uniform_filter

# 窗口大小
kernel_size = 3
popu_size = 30
max_iter = 20
dim = 4
F = 0.5
CR = 0.5

def generate_population(popu_size):
    x = []
    for i in range(popu_size):
        x.append(
            [
                random.uniform(0, 1.5),
                random.uniform(0, 0.5),
                random.random(),
                random.uniform(0.5, 1.5),
            ]
        )
    return x

# 窗口标准差
def window_stdev(X, window_size):
    (
        r,
        c,
    ) = X.shape
    X += np.random.rand(r, c) * 1e-6
    c1 = uniform_filter(X, window_size, mode="reflect")
    c2 = uniform_filter(X * X, window_size, mode="reflect")
    return np.sqrt(c2 - c1 * c1)

# 可以实现灰度变换
def ok_transform(img, a, b, c, k, k_size):
    new_img = copy.deepcopy(img).astype(np.float)
    D = np.mean(new_img)
    M = cv2.blur(new_img, (k_size, k_size))
    S = window_stdev(new_img, k_size)
    # S = generic_filter(new_img, np.std, size=k_size)
    new_img = k * D / (S + b) * (new_img - c * M) + np.power(M, a)
    return np.uint8(np.clip(new_img, 0, 255))  # 转换成0-255之间的数字

# 计算图像的熵
def calc_entropy(img):
    rows = img.shape[0]
    cols = img.shape[1]
    tmp = [0 for i in range(256)]
    val = 0
    entropy = 0
    # 统计不同灰度值像素的个数
    for i in range(rows):
        for j in range(cols):
            val = img[i][j]
            tmp[val] += 1
    # 计算每个灰度值像素出现的概率
    P = [0 for i in range(256)]
    P = [float(i / (rows * cols)) for i in tmp]
    for i in range(256):
        if P[i] == 0:
            continue
        else:
            entropy = float(entropy - P[i] * math.log2(P[i]) / 1)
    return entropy

# 计算适应度
def calc_fitness(img):
    rows, cols = img.shape
    E = np.sum(img)
    sobelX = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3)
    sobelY = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=3)
    absx = cv2.convertScaleAbs(sobelX)
    absy = cv2.convertScaleAbs(sobelY)
    sobel = cv2.addWeighted(absx, 0.5, absy, 0.5, 0)
    ret, threshold = cv2.threshold(
        sobel, 100, 255, cv2.THRESH_BINARY_INV
    )  # 转化为二值化图像，图像阈值为100
    n_edgels = len(threshold[threshold == 255])  # 二值化图像大于阈值的像素个数
    entropy = calc_entropy(img)
    return math.log(math.log(E)) * n_edgels / (rows * cols) * entropy

def DEIE(popu_size, max_iter, dim, img, calc_fitness, k_size, F, CR):
    X = np.array(generate_population(popu_size))  # 初始化粒子群位置
    V = [[0 for i in range(dim)] for j in range(popu_size)]
    U = [[0 for i in range(dim)] for j in range(popu_size)]

    print("====================初始化适应度====================")
    fitness = [
        calc_fitness(ok_transform(img, x[0], x[1], x[2], x[3], k_size)) for x in X
    ]
    print("====================适应度初始化完成====================")
    print("fitness", fitness)
    best = X[np.argmax(fitness)]  # 初始化全局最优位置
    best_history = [best.tolist()]
    figure_count = 0
    print("====================开始迭代====================")
    for iter in range(max_iter):
        print("====================第{}次迭代====================".format(iter))
        for i in range(popu_size):
            img_out = ok_transform(img, X[i][0], X[i][1], X[i][2], X[i][3], k_size)
            figure_count = figure_count + 1
            cv2.imwrite("result/rice/figure{}.png".format(figure_count), img_out)
            new_fitness = calc_fitness(img_out)
            if new_fitness > fitness[i]:
                best = X[i].copy()

            # mutation
            # r1 r2 r3 random and different
            r1 = random.randrange(popu_size)
            while r1 == i:
                r1 = random.randrange(popu_size)
            r2 = random.randrange(popu_size)
            while r2 == i or r2 == r1:
                r2 = random.randrange(popu_size)
            r3 = random.randrange(popu_size)
            while r3 == i or r3 == r1 or r3 == r2:
                r3 = random.randrange(popu_size)

            for j in range(dim):
                V[i][j] = X[r1][j] + F * (X[r2][j] - X[r3][j])
                # limit of bound
                # V[i][j] = (
                #     x_low_bound if V[i][j] < x_low_bound else V[i][j]
                # )
                # V[i][j] = (
                #     x_up_bound if V[i][j] > x_up_bound else V[i][j]
                # )

            # crossover
            for j in range(dim):
                U[i][j] = V[i][j] if random.random() <= CR or j == random.randrange(dim) else X[i][j]
                # limit of bound
                # U[i][j] = (x_low_bound if U[i][j] < x_low_bound else U[i][j]
                # )
                # U[i][j] = (
                #     x_up_bound if U[i][j] > x_up_bound else U[i][j]
                # )

            # select
            # x_fitness = calc_fitness(X[i])
            u_fitness = calc_fitness(ok_transform(img, U[i][0], U[i][1], U[i][2], U[i][3], k_size))
            if u_fitness > fitness[i]:
                X[i] = U[i].copy()
                # for j in range(dim):
                #     X[i][j] = U[i][j]

            # x_fitness = calc_fitness(X[i])
            # if x_fitness > global_fitness:
            #     global_fitness = x_fitness
            #     global_solution = X[i]
        print("best: ", best)
        best_history.append(best.tolist())
        print("更新后粒子的位置：")
        print(X)
        # X[:, 0] = np.clip(X[:, 0], 0, 1.5)
        # X[:, 1] = np.clip(X[:, 1], 0, 0.5)
        # X[:, 2] = np.clip(X[:, 2], 0, 1)
        # X[:, 3] = np.clip(X[:, 3], 0.5, 1.5)
        # print("修改范围后粒子的位置：")
        # print(X)
    print("history of best:")
    for best in best_history:
        print(best)
    return best

# 读入原始图像
img = cv2.imread("../images/rice.png", 0)
DEIE(popu_size, max_iter, dim, img, calc_fitness, kernel_size, F, CR)

# 通过窗口展示图片 第一个参数为窗口名 第二个为读取的图片变量
# cv2.imshow("img", img)

# cv2.waitKey(0)
# cv2.destroyAllWindows()

max_stars_repo_name: zxhyJack/image-enhancement

max_issues_repo_name: zxhyJack/image-enhancement

max_forks_repo_name: zxhyJack/image-enhancement

reason: import numpy,from scipy


text:
# Databricks notebook source
# MAGIC %md # W261 Final Project - Airline Delays and Weather Modeling Notebook
# MAGIC __`MIDS w261: Machine Learning at Scale | UC Berkeley School of Information | Summer 2020`__
# MAGIC ### Team 1:
# MAGIC * Sayan Das
# MAGIC * Kevin Hartman
# MAGIC * Hersh Solanki
# MAGIC * Nick Sylva

# COMMAND ----------

# MAGIC %md # Table of Contents
# MAGIC ### 6. Pre-Model Data Prep
# MAGIC #### a. Configure access to post-processed data from the silver data pipeline
# MAGIC #### b. Split into Train/Validation/Test
# MAGIC #### c. Apply specific cleanup required before modelling (e.g. imputation)
# MAGIC #### d. Save the model-ready datasets in silver
# MAGIC ### 7. Model Exploration (Pipeline)
# MAGIC #### a. Feature Selection
# MAGIC #### b. Feature Engineering
# MAGIC #### c. Transformations (Encoding & Scaling)
# MAGIC #### d. Evaluation
# MAGIC ### 8. Model Selection and Tuning

# COMMAND ----------

# MAGIC %md **Note that this notebook consists of our work modeling three main algorithms: Logistic Regression, Random Forest, and Gradient-Boosted Trees.**

# COMMAND ----------

# MAGIC %md ## Notebook Imports

# COMMAND ----------

!pip install networkx
from pyspark.sql import functions as f
from pyspark.sql.types import StructType, StructField, StringType, DoubleType, IntegerType, NullType, ShortType, DateType, BooleanType, BinaryType
from pyspark.sql import SQLContext
from math import radians, cos, sin, asin, sqrt
import matplotlib.pyplot as plt
import time
import numpy as np
import pandas as pd
import seaborn as sns
from pytz import timezone
from datetime import  datetime, timedelta
import os
import networkx as nx
from delta.tables import DeltaTable

# Model Imports
from pyspark.ml import Pipeline
from pyspark.ml.feature import StringIndexer, VectorAssembler, StandardScaler, OneHotEncoder, PCA, FeatureHasher
from pyspark.ml.classification import LogisticRegression
from pyspark.ml.classification import RandomForestClassifier
from pyspark.ml.classification import GBTClassifier

%matplotlib inline
sqlContext = SQLContext(sc)

# COMMAND ----------

# MAGIC %md ### 6a. Configure access to data pipeline staging areas

# COMMAND ----------

username = "kevin"
dbutils.widgets.text("username", username)
spark.sql(f"CREATE DATABASE IF NOT EXISTS airline_delays_{username}")
spark.sql(f"USE airline_delays_{username}")

flights_and_weather_pipeline_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_pipeline/"
flights_and_weather_train_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_train/"
flights_and_weather_validation_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_validation/"
flights_and_weather_test_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_test/"

spark.conf.set("spark.sql.shuffle.partitions", 80)

# COMMAND ----------

# MAGIC %md
# MAGIC #### If running this notebook for the first time, continue with section 6. Otherwise, jump to section 7.

# COMMAND ----------

# MAGIC %md
# MAGIC # 6. Pre-Model Data Prep

# COMMAND ----------

# MAGIC %md ## 6b. Split data into train, validation, test

# COMMAND ----------

train_data = spark.sql("SELECT * FROM flights_and_weather_pipeline_processed WHERE YEAR IN (2015, 2016, 2017)").drop('YEAR')
validation_data = spark.sql("SELECT * FROM flights_and_weather_pipeline_processed WHERE YEAR = 2018").drop('YEAR')
test_data = spark.sql("SELECT * FROM flights_and_weather_pipeline_processed WHERE YEAR = 2019").drop('YEAR')

# COMMAND ----------

# MAGIC %md ## 6c. Apply specific cleanup required before modeling (e.g. imputation)

# COMMAND ----------

# MAGIC %md ### Prepare an imputation dictionary to use to replace missing values

# COMMAND ----------

def make_imputation_dict(df):
  impute_dict = {}

  cat_cols = [item[0] for item in df.dtypes if item[1].startswith('string')]  # string
  num_cols = [item[0] for item in df.dtypes if item[1].startswith('int') | item[1].startswith('double')] # will select name of column with integer or double data type

  for x in cat_cols:
    mode = df.groupBy(x).count().sort(f.col("count").desc()).collect()
    impute_dict[x] = mode[0][0]

  # Fill the missing numerical values with the average of each #column
  for i in num_cols:
    mean_value = df.select(f.mean(i).cast(DoubleType())).collect()
    impute_dict[i] = mean_value[0][0]

  return impute_dict

# COMMAND ----------

# MAGIC %md Dictionary should be populated from ONLY the train set to prevent leakage into validation and test sets.

# COMMAND ----------

impute_dict = make_imputation_dict(train_data)

# COMMAND ----------

# MAGIC %md Impute our missing values using the mean and modes found in our impute dictionary

# COMMAND ----------

def impute_missing_values(df, impute_dict):
  missing_count_list = []
  for c in df.columns:
      if df.where(f.col(c).isNull()).count() > 0:
          tup = (c,int(df.where(f.col(c).isNull()).count()))
          missing_count_list.append(tup)

  missing_column_list = [x[0] for x in missing_count_list]
  missing_df = df.select(missing_column_list)

  missing_cat_columns = [item[0] for item in missing_df.dtypes if item[1].startswith('string')]  # string
  print("\nCategorical Columns with missing data:", missing_cat_columns)

  missing_num_columns = [item[0] for item in missing_df.dtypes if item[1].startswith('int') | item[1].startswith('double')] # will select name of column with integer or double data type
  print("\nNumerical Columns with missing data:", missing_num_columns)

  # Fill the missing categorical values with the most frequent category
  for x in missing_cat_columns:
    mode = impute_dict[x]
    if mode:
      df = df.withColumn(x, f.when(df[x].isNull(), f.lit(mode)).otherwise(df[x]))
    else:
      df = df.withColumn(x, f.when(df[x].isNull(), 'None').otherwise(df[x]))

  # Fill the missing numerical values with the average of each #column
  for i in missing_num_columns:
    mean_value = impute_dict[x]
    if mean_value:
        df = df.withColumn(i, f.when(df[i].isNull(), mean_value).otherwise(df[i]))
    else:
        df = df.withColumn(i, f.when(df[i].isNull(), 0).otherwise(df[i]))

  return df

# COMMAND ----------

print(impute_dict)

# COMMAND ----------

# MAGIC %md Impute our training data

# COMMAND ----------

imputed_train_data = impute_missing_values(train_data, impute_dict)

# COMMAND ----------

# MAGIC %md Let's save our train data that has been imputed to the Delta Lake.

# COMMAND ----------

flights_and_weather_train_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_train/"
dbutils.fs.rm(flights_and_weather_train_loc + 'processed', recurse=True)
imputed_train_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_train_loc + 'processed')

# COMMAND ----------

# MAGIC %sql
# MAGIC
# MAGIC DROP TABLE IF EXISTS flights_and_weather_train_processed;
# MAGIC
# MAGIC CREATE TABLE flights_and_weather_train_processed
# MAGIC USING DELTA
# MAGIC LOCATION "/airline_delays/$username/DLRS/flights_and_weather_train/processed"

# COMMAND ----------

# MAGIC %md Next let's impute our validation data based on values stored in our dictionary (this is unseen data)

# COMMAND ----------

imputed_validation_data = impute_missing_values(validation_data, impute_dict)

# COMMAND ----------

# MAGIC %md Let's save our validation data that has been imputed to the Delta Lake.

# COMMAND ----------

flights_and_weather_validation_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_validation/"
dbutils.fs.rm(flights_and_weather_validation_loc + 'processed', recurse=True)
imputed_validation_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_validation_loc + 'processed')

# COMMAND ----------

# MAGIC %sql
# MAGIC
# MAGIC DROP TABLE IF EXISTS flights_and_weather_validation_processed;
# MAGIC
# MAGIC CREATE TABLE flights_and_weather_validation_processed
# MAGIC USING DELTA
# MAGIC LOCATION "/airline_delays/$username/DLRS/flights_and_weather_validation/processed"

# COMMAND ----------

# MAGIC %md Next let's impute our test data based on values stored in our dictionary (this is unseen data)

# COMMAND ----------

imputed_test_data = impute_missing_values(test_data, impute_dict)

# COMMAND ----------

# MAGIC %md Let's save our test data that has been imputed to the Delta Lake.

# COMMAND ----------

flights_and_weather_test_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_test/"
dbutils.fs.rm(flights_and_weather_test_loc + 'processed', recurse=True)
imputed_test_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_test_loc + 'processed')

# COMMAND ----------

# MAGIC %sql
# MAGIC
# MAGIC DROP TABLE IF EXISTS flights_and_weather_test_processed;
# MAGIC
# MAGIC CREATE TABLE flights_and_weather_test_processed
# MAGIC USING DELTA
# MAGIC LOCATION "/airline_delays/$username/DLRS/flights_and_weather_test/processed"

# COMMAND ----------

imputed_train_data = spark.sql("select * from flights_and_weather_train_processed VERSION as of 0")
imputed_validation_data = spark.sql("select * from flights_and_weather_validation_processed VERSION as of 0")
imputed_test_data = spark.sql("select * from flights_and_weather_test_processed VERSION as of 0")

# COMMAND ----------

# MAGIC %md ### Replace all empty strings

# COMMAND ----------

def replace_empty_strings(df):
  cat_columns = [item[0] for item in df.dtypes if item[1].startswith('string')]  # string

  # Fill the empty string values with 'Empty'
  for x in cat_columns:
     df = df.withColumn(x, f.when(df[x] == '', f.lit("Empty")).otherwise(df[x]))

  return df

# COMMAND ----------

cleaned_train_data = replace_empty_strings(imputed_train_data)
cleaned_validation_data = replace_empty_strings(imputed_validation_data)
cleaned_test_data = replace_empty_strings(imputed_test_data)

# COMMAND ----------

cleaned_train_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_train_loc + 'processed')
cleaned_validation_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_validation_loc + 'processed')
cleaned_test_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_test_loc + 'processed')

# COMMAND ----------

# MAGIC %md ### Drop Numeric Columns that are all Null

# COMMAND ----------

def get_cols_to_drop(impute_dict, numeric_cols):
  colsToDrop = []
  numeric_cols = set(numeric_cols)
  for key, value in impute_dict.items():
    if value == None and key in numeric_cols:
      colsToDrop.append(key)
  return colsToDrop

numeric_cols = [item[0] for item in train_data.dtypes if item[1] in ['int', 'double']]
cols_to_drop = get_cols_to_drop(impute_dict, numeric_cols)

#cols_to_drop = ['ORIGIN_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'ORIGIN_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'ORIGIN_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'ORIGIN_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'DEST_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'DEST_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'DEST_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION']

# COMMAND ----------

train_data = train_data.drop(*cols_to_drop)
validation_data = validation_data.drop(*cols_to_drop)
test_data = test_data.drop(*cols_to_drop)

# COMMAND ----------

# MAGIC %md ## 6d. Save datasets into silver

# COMMAND ----------

train_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_train_loc + 'processed')
validation_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_validation_loc + 'processed')
test_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_test_loc + 'processed')

# COMMAND ----------

# MAGIC %md ### Drop Delay Columns (revisit if we want to use a linear regression model to predict the actual time delayed)

# COMMAND ----------

cols_to_drop = ['DEP_DELAY', 'DEP_DELAY_NEW', 'DEP_DELAY_GROUP', 'TAIL_NUM', 'ORIGIN_CITY_NAME', 'DEST_CITY_NAME', 'FL_DATE', 'ORIGIN_WEATHER_KEY', 'DEST_WEATHER_KEY']

# COMMAND ----------

train_data = train_data.drop(*cols_to_drop)
validation_data = validation_data.drop(*cols_to_drop)
test_data = test_data.drop(*cols_to_drop)

# COMMAND ----------

train_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_train_loc + 'processed')
validation_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_validation_loc + 'processed')
test_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_test_loc + 'processed')

# COMMAND ----------

# MAGIC %sql
# MAGIC DESCRIBE HISTORY flights_and_weather_train_processed

# COMMAND ----------

# MAGIC %sql
# MAGIC DESCRIBE HISTORY flights_and_weather_validation_processed

# COMMAND ----------

# MAGIC %sql
# MAGIC DESCRIBE HISTORY flights_and_weather_test_processed

# COMMAND ----------

# COMMAND ----------

# MAGIC %sql
# MAGIC select * from flights_and_weather_train_processed VERSION as of 3

# COMMAND ----------

# MAGIC %md # 7. Model Pipeline for Exploration

# COMMAND ----------

# MAGIC %md Let's read in our data from silver for modeling.

# COMMAND ----------

train_data = spark.sql("select * from flights_and_weather_train_processed")
validation_data = spark.sql("select * from flights_and_weather_validation_processed")
test_data = spark.sql("select * from flights_and_weather_test_processed")

cols_to_drop = ['ORIGIN', 'DEST', 'ORIGIN_IATA', 'DEST_IATA', 'DEP_DELAY', 'DEP_DELAY', 'DEP_DELAY_NEW', 'DEP_DELAY_GROUP', 'TAIL_NUM', 'ORIGIN_CITY_NAME', 'DEST_CITY_NAME', 'FL_DATE', 'ORIGIN_WEATHER_KEY', 'DEST_WEATHER_KEY', 'ORIGIN_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'ORIGIN_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'ORIGIN_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'ORIGIN_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'DEST_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'DEST_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'DEST_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION']

train_data = train_data.drop(*cols_to_drop)
validation_data = validation_data.drop(*cols_to_drop)
test_data = test_data.drop(*cols_to_drop)

# COMMAND ----------

# MAGIC %md ### Routines to encode and scale our features
# MAGIC
# MAGIC We encode and scale only on train data and set up the transformation pipeline with those stages so it may be applied on validate and test data later

# COMMAND ----------

# MAGIC %md ### Encoding Function

# COMMAND ----------

# Use the OneHotEncoderEstimator to convert categorical features into one-hot vectors
# Use VectorAssembler to combine vector of one-hots and the numerical features
# Append the process into the stages array to reproduce

def create_encoding_stages(data, label_name):

  cat_cols = [item[0] for item in data.dtypes if item[1].startswith('string')]
  numeric_cols = [item[0] for item in data.dtypes if item[1].startswith('int') | item[1].startswith('double')]
  cat_cols.remove(label_name)
  hash_cols = ['ORIGIN_AW1_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER', 'ORIGIN_AW2_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER', \
               'ORIGIN_AW3_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER', 'ORIGIN_AW4_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER',  \
               'DEST_AW1_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER', 'DEST_AW2_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER', \
               'DEST_AW3_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER', 'DEST_AW4_PRESENT_WEATHER_OBSERVATION_AUTOMATED_OCCURRENCE_IDENTIFIER']
  for c in hash_cols:
     cat_cols.remove(c)

  # Encode certain categorical columns with FeatureHasher
  hash_pipeline = Pipeline(stages=[
    FeatureHasher(inputCols=hash_cols, outputCol="hash_features")
  ])

  # One-Hot Encode Categorical Columns in the vector
  cat_pipeline = Pipeline(stages=[
    StringIndexer(inputCols=cat_cols, outputCols=[c + '_idx' for c in cat_cols] , handleInvalid = 'keep') ,
    OneHotEncoder(inputCols= [c + '_idx' for c in cat_cols], outputCols = [c + '_enc' for c in cat_cols]),
    VectorAssembler(inputCols=[c + '_enc' for c in cat_cols], outputCol="cat_features")
   ])

  # Deal with Numeric Features
  num_pipeline = Pipeline(stages=[
    VectorAssembler(inputCols = numeric_cols, outputCol = 'unscaled_features', handleInvalid = 'keep'),
    StandardScaler(inputCol = 'unscaled_features', outputCol = 'scaled_features', withStd = True, withMean = True)
  ])

  feature_assembler = Pipeline(stages=[
    hash_pipeline,
    cat_pipeline,
    num_pipeline,
    VectorAssembler(inputCols = ["hash_features", "cat_features", "scaled_features"], outputCol = 'features', handleInvalid = 'keep')
  ])

  label_string_indexer = StringIndexer(inputCol = label_name, outputCol = 'label')

  stages = [feature_assembler] + [label_string_indexer]

  return stages

# COMMAND ----------

# Use the OneHotEncoderEstimator to convert categorical features into one-hot vectors
# Use VectorAssembler to combine vector of one-hots and the numerical features
# Append the process into the stages array to reproduce

def create_encoding_stages(data, label_name):

  cat_cols = [item[0] for item in data.dtypes if item[1].startswith('string')]
  numeric_cols = [item[0] for item in data.dtypes if item[1].startswith('int') | item[1].startswith('double')]
  cat_cols.remove(label_name)

  # Encode certain categorical columns with FeatureHasher
  #cat_pipeline = Pipeline(stages=[
  #  FeatureHasher(inputCols=cat_cols, outputCol="cat_features")
  # ])

  # One-Hot Encode Categorical Columns in the vector
  cat_pipeline = Pipeline(stages=[
    StringIndexer(inputCols=cat_cols, outputCols=[c + '_idx' for c in cat_cols] , handleInvalid = 'keep') ,
    OneHotEncoder(inputCols= [c + '_idx' for c in cat_cols], outputCols = [c + '_enc' for c in cat_cols]),
    VectorAssembler(inputCols=[c + '_enc' for c in cat_cols], outputCol="cat_features")
   ])

  # Deal with Numeric Features
  num_pipeline = Pipeline(stages=[
    VectorAssembler(inputCols = numeric_cols, outputCol = 'unscaled_features', handleInvalid = 'keep'),
    StandardScaler(inputCol = 'unscaled_features', outputCol = 'scaled_features', withStd = True, withMean = True)
  ])

  feature_assembler = Pipeline(stages=[
    cat_pipeline,
    num_pipeline,
    VectorAssembler(inputCols = ["cat_features", "scaled_features"], outputCol = 'features', handleInvalid = 'keep')
  ])

  label_string_indexer = StringIndexer(inputCol = label_name, outputCol = 'label')

  stages = [feature_assembler] + [label_string_indexer]

  return stages

# COMMAND ----------

# MAGIC %md ### Routines to evaluate our model performance
# MAGIC
# MAGIC We train our models and evaluate performance throughout the pipeline

# COMMAND ----------

## This function will calculate evaluation metrics based on predictions
def evaluation_metrics(predictions, model_name):
  predictions = predictions.createOrReplaceTempView('predictions')
  display(predictions)
  true_positives = spark.sql("SELECT COUNT(*) FROM predictions WHERE predictions.label = 1 AND predictions.prediction = 1").collect()[0][0]
  true_negatives = spark.sql("SELECT COUNT(*) FROM predictions WHERE predictions.label = 0 AND predictions.prediction = 0").collect()[0][0]
  false_positives = spark.sql("SELECT COUNT(*) FROM predictions WHERE predictions.label = 0 AND predictions.prediction = 1").collect()[0][0]
  false_negatives = spark.sql("SELECT COUNT(*) FROM predictions WHERE predictions.label = 1 AND predictions.prediction = 0").collect()[0][0]

  # Now we should compute our main statistics
  accuracy = (true_positives + true_negatives)/(true_positives + true_negatives + false_positives + false_negatives)
  recall = (true_positives) / (true_positives + false_negatives)
  precision = (true_positives) / (true_positives + false_positives) if (true_positives + false_positives) != 0 else 0
  f1_score = (2 * recall * precision) / (recall + precision) if (recall + precision) != 0 else 0

  print("The accuracy is: %s" % np.round(accuracy, 4))
  print("The recall is: %s" % np.round(recall, 4))
  print("The precision is: %s" % np.round(precision, 4))
  print("The f1_score is: %s" % np.round(f1_score, 4))
  df_cols = ['Model','accuracy', 'precision', 'recall', 'f1_score', 'true_positives', 'true_negatives', 'false_positives', 'false_negatives']
  metrics_dataframe = pd.DataFrame(columns = df_cols)
  result_dict = {'Model': model_name,
                 'accuracy': accuracy,
                 'precision': precision,
                 'recall': recall,
                 'f1_score': f1_score,
                 'true_positives': true_positives,
                 'true_negatives': true_negatives,
                 'false_positives': false_positives,
                 'false_negatives': false_negatives}
  metrics_dataframe = metrics_dataframe.append(result_dict, ignore_index=True)

  return metrics_dataframe

def confusion_matrix(model_results,index):
  confusion_matrix = pd.DataFrame(columns = ['_','Predicted Delay', 'Predicted No Delay'])
  conf_matrix_delay = {'_': 'Actual Delay', 'Predicted Delay': model_results['true_positives'][index], 'Predicted No Delay': model_results['false_negatives'][index] }
  conf_matrix_no_delay = {'_': 'Actual  No Delay', 'Predicted Delay': model_results['false_positives'][index], 'Predicted No Delay': model_results['true_negatives'][index] }
  confusion_matrix = confusion_matrix.append(conf_matrix_delay,ignore_index=True)
  confusion_matrix = confusion_matrix.append(conf_matrix_no_delay,ignore_index=True)
  return confusion_matrix

def AUC(model, predictions):
  from pyspark.ml.evaluation import BinaryClassificationEvaluator
  import matplotlib.pyplot as plt
  evaluator = BinaryClassificationEvaluator()
  evaluation = evaluator.evaluate(predictions)
  print("evaluation (area under ROC): %f" % evaluation)

def ROC(model, predictions):
  plt.figure(figsize=(10,10))
  plt.plot([0, 1], [0, 1], 'r--')
  plt.plot(model.summary.roc.select('FPR').collect(),
         model.summary.roc.select('TPR').collect())
  plt.xlabel('FPR')
  plt.ylabel('TPR')
  plt.show()

def showPR(model_results):
  pr = model_results.toPandas()
  ax.plot(pr['recall'], pr['precision'])
  ax.set_xlabel('Precision')
  ax.set_ylabel('Recall')
  ax.set_title('Precision Recall')

# COMMAND ----------

# MAGIC %md # Model 1: Baseline Evaluation (Logistic Regression)
# MAGIC
# MAGIC Set up and run the pipeline for the baseline model (using select highest correlated features)

# COMMAND ----------

correlated_cols = ['ORIGIN_WND_SPEED_RATE', 'ORIGIN_AA2_LIQUID_PRECIPITATION_PERIOD_QUANTITY', 'ORIGIN_SNOW_DEPTH_DIMENSION',
 'ORIGIN_AL1_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'ORIGIN_RUNWAY_VISUAL_RANGE_OBSERVATION_DIRECTION_ANGLE',
 'ORIGIN_RUNWAY_VISUAL_RANGE_OBSERVATION_VISIBILITY_DIMENSION', 'ORIGIN_GA2_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION',
 'ORIGIN_SKY_CONDITION_OBSERVATION_LOWEST_CLOUD_BASE_HEIGHT_DIMENSION',
 'DEST_WND_SPEED_RATE', 'DEST_AA2_LIQUID_PRECIPITATION_PERIOD_QUANTITY', 'DEST_AA3_LIQUID_PRECIPITATION_DEPTH_DIMENSION',
 'DEST_SNOW_DEPTH_DIMENSION', 'DEST_AL1_SNOW_ACCUMULATION_DEPTH_DIMENSION',
 'DEST_GA2_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION',  'DEST_GD4_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION',
 'DEST_SKY_CONDITION_OBSERVATION_LOWEST_CLOUD_BASE_HEIGHT_DIMENSION', 'DEP_DEL15']

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data[correlated_cols],'DEP_DEL15'))
encoding_pipeline = encoding_pipeline.fit(train_data[correlated_cols])
# apply the transformations to our train data
transformed_train_data = encoding_pipeline.transform(train_data[correlated_cols])['features', 'label']

# COMMAND ----------

# train a model on our transformed train data
startTime = time.time()
lr = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.001, elasticNetParam = 0.25, standardization = False)
model = lr.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Logistic Regression on training data")
display(train_metrics)

# COMMAND ----------

AUC(model, train_preds)

# COMMAND ----------

ROC(model, train_preds)

# COMMAND ----------

# MAGIC %md #### Evaluate against validation

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
transformed_validation_data = encoding_pipeline.transform(validation_data[correlated_cols])['features', 'label']

# run the fitted model on the transformed validation data
validation_preds = model.transform(transformed_validation_data)

# display our evaluation metrics
validation_metrics = evaluation_metrics(validation_preds, "Logistic Regression on validation data")
display(validation_metrics)

# COMMAND ----------

# MAGIC %md # Parameter Search

# COMMAND ----------

# transformed_train_data_sample = transformed_train_data.sample(False, 0.0001)
from pyspark.ml import Pipeline
from pyspark.ml.classification import LogisticRegression
from pyspark.ml.evaluation import BinaryClassificationEvaluator
from pyspark.ml.tuning import CrossValidator, ParamGridBuilder

def parameter_search(model, featuresCol, labelCol, dictionary, data):
  mod = model(featuresCol = featuresCol, labelCol = labelCol)
  pipeline = Pipeline(stages=[mod])

  paramGrid = ParamGridBuilder()
  for i in range(len(dictionary)):
    key = list(dictionary.keys())[i]
    value = list(dictionary.values())[i]
    paramGrid = paramGrid.addGrid(getattr(mod,key) , value)

  paramGrid = paramGrid.build()

  crossval = CrossValidator(estimator=pipeline,
                            estimatorParamMaps=paramGrid,
                            evaluator=BinaryClassificationEvaluator(),
                            numFolds=2)

  cvModel = crossval.fit(data)
  return(cvModel.getEstimatorParamMaps()[np.argmax(cvModel.avgMetrics)])

# COMMAND ----------

parameters = {'maxIter': [1, 30], 'regParam' : [0.01, 0.001], 'elasticNetParam' : [0.25, 0.5, 0.1]}
parameter_search(LogisticRegression, 'features', 'label', parameters, transformed_train_data.sample(False, 0.005))

# COMMAND ----------

# MAGIC %md
# MAGIC     # Results
# MAGIC     {Param(parent='LogisticRegression_6a0dd178a09d', name='maxIter', doc='max number of iterations (>= 0).'): 30,
# MAGIC      Param(parent='LogisticRegression_6a0dd178a09d', name='regParam', doc='regularization parameter (>= 0).'): 0.01,
# MAGIC      Param(parent='LogisticRegression_6a0dd178a09d', name='elasticNetParam',
# MAGIC      doc='the ElasticNet mixing parameter, in range [0, 1]. For alpha = 0, the penalty is an L2 penalty. For alpha = 1, it is an L1 penalty.'): 0.5}

# COMMAND ----------

# MAGIC %md
# MAGIC ### LR Model using results from parameter search

# COMMAND ----------

# Using results from our parameter search
startTime = time.time()
lr = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.01, elasticNetParam = 0.5, standardization = False)
model = lr.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")
train_metrics = evaluation_metrics(train_preds, "Logistic Regression on training data")
display(train_metrics)

# COMMAND ----------

AUC(model, train_preds)
ROC(model, train_preds)

# COMMAND ----------

# run the fitted model on the transformed validation data
validation_preds = model.transform(transformed_validation_data)
# display our evaluation metrics
validation_metrics = evaluation_metrics(validation_preds, "Logistic Regression on validation data")
display(validation_metrics)

# COMMAND ----------

# MAGIC %md Our precision is much worse with these parameters. The 5% sample is not enough to conduct a proper CV. We will review a full CV when our final model is selected later.

# COMMAND ----------

# MAGIC %md # Feature Selection and Engineering for next models

# COMMAND ----------

# MAGIC %md ## PageRank Features (PR is based on train data only)

# COMMAND ----------

#PageRank
flights_and_weather_train_processed = spark.sql("SELECT * from flights_and_weather_train_processed")
airlineGraph = {'nodes': flights_and_weather_train_processed.select('ORIGIN', 'DEST').rdd.flatMap(list).distinct().collect(),
                'edges': flights_and_weather_train_processed.select('ORIGIN', 'DEST').rdd.map(tuple).collect()}

directedGraph = nx.DiGraph()
directedGraph.add_nodes_from(airlineGraph['nodes'])
directedGraph.add_edges_from(airlineGraph['edges'])

pageRank = nx.pagerank(directedGraph, alpha = 0.85)
pandasPageRank = pd.DataFrame(pageRank.items(), columns = ['Station', 'PageRank'])
pandasPageRank = spark.createDataFrame(pandasPageRank)

# COMMAND ----------

pagerank_origin = pandasPageRank.withColumnRenamed('Station','ORIGIN_IATA').withColumnRenamed('PageRank','ORIGIN_PAGERANK')
pagerank_dest= pandasPageRank.withColumnRenamed('Station','DEST_IATA').withColumnRenamed('PageRank','DEST_PAGERANK')

# COMMAND ----------

pagerank_loc = f"/airline_delays/{username}/DLRS/pagerank/"

dbutils.fs.rm(pagerank_loc + 'origin', recurse=True)
dbutils.fs.rm(pagerank_loc + 'dest', recurse=True)

# COMMAND ----------

pagerank_origin.write.option('mergeSchema', True).mode('overwrite').format('delta').save(pagerank_loc + 'origin')
pagerank_dest.write.option('mergeSchema', True).mode('overwrite').format('delta').save(pagerank_loc + 'dest')

# COMMAND ----------

# MAGIC %sql
# MAGIC
# MAGIC DROP TABLE IF EXISTS pagerank_origin;
# MAGIC
# MAGIC CREATE TABLE pagerank_origin
# MAGIC USING DELTA
# MAGIC LOCATION "/airline_delays/$username/DLRS/pagerank/origin"

# COMMAND ----------

# MAGIC %sql
# MAGIC
# MAGIC DROP TABLE IF EXISTS pagerank_dest;
# MAGIC
# MAGIC CREATE TABLE pagerank_dest
# MAGIC USING DELTA
# MAGIC LOCATION "/airline_delays/$username/DLRS/pagerank/dest"

# COMMAND ----------

# MAGIC %md
# MAGIC ## Rolling Average Features

# COMMAND ----------

train_data = spark.sql('SELECT * FROM flights_and_weather_train_ra_processed')
validation_data = spark.sql('SELECT * FROM flights_and_weather_validation_ra_processed')
test_data = spark.sql('SELECT * FROM flights_and_weather_test_ra_processed')

# COMMAND ----------

train_data = replace_empty_strings(train_data)
validation_data = replace_empty_strings(validation_data)
test_data = replace_empty_strings(test_data)

# COMMAND ----------

flights_and_weather_train_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_ra_train/"
flights_and_weather_validation_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_ra_validation/"
flights_and_weather_test_loc = f"/airline_delays/{username}/DLRS/flights_and_weather_ra_test/"

# COMMAND ----------

train_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_train_loc + 'processed')
validation_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_validation_loc + 'processed')
test_data.write.option('mergeSchema', True).mode('overwrite').format('delta').save(flights_and_weather_test_loc + 'processed')

# COMMAND ----------

cols_to_drop = ['ORIGIN', 'DEST', 'ORIGIN_IATA', 'DEST_IATA', 'DEP_DELAY', 'DEP_DELAY_NEW', 'DEP_DELAY_GROUP', 'TAIL_NUM', 'ORIGIN_CITY_NAME', 'DEST_CITY_NAME', 'FL_DATE', 'ORIGIN_WEATHER_KEY', 'DEST_WEATHER_KEY', 'ORIGIN_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'ORIGIN_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'ORIGIN_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'ORIGIN_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'DEST_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'DEST_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'DEST_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION']

train_data = train_data.drop(*cols_to_drop)
validation_data = validation_data.drop(*cols_to_drop)
test_data = test_data.drop(*cols_to_drop)

# COMMAND ----------

# COMMAND ----------

# MAGIC %md # Model Attempts

# COMMAND ----------

# MAGIC %md
# MAGIC ## Gradient Boosted Trees

# COMMAND ----------

# MAGIC %md #### GBT on limited correlated columns

# COMMAND ----------

startTime = time.time()
gbt = GBTClassifier(labelCol = "label", featuresCol = "features", maxIter = 10)
gbt_model = gbt.fit(transformed_train_data)
gbt_train_preds = gbt_model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Gradient-Boosted Tree model is: {(endTime - startTime) / (60)} minutes")
gbt_train_metrics = evaluation_metrics(gbt_train_preds, "Gradient-Boosted Trees on training data")
display(gbt_train_metrics)

# COMMAND ----------

AUC(model, train_preds)
ROC(model, train_preds)

# COMMAND ----------

gbt_validation_preds = gbt_model.transform(transformed_validation_data)
gbt_validation_metrics = evaluation_metrics(gbt_validation_preds, "Gradient-Boosted Trees on validation data")
display(gbt_validation_metrics)

# COMMAND ----------

AUC(model, train_preds)
ROC(model, train_preds)

# COMMAND ----------

# MAGIC %md #### GBT on full set of columns

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data,'DEP_DEL15'))
encoding_pipeline = encoding_pipeline.fit(train_data)
# apply the transformations to our train data
transformed_train_data = encoding_pipeline.transform(train_data)['features', 'label']

# COMMAND ----------

startTime = time.time()
gbt = GBTClassifier(labelCol = "label", featuresCol = "features", maxIter = 10)
gbt_model = gbt.fit(transformed_train_data)
gbt_train_preds = gbt_model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Gradient-Boosted Tree model is: {(endTime - startTime) / (60)} minutes")
gbt_train_metrics = evaluation_metrics(gbt_train_preds, "Gradient-Boosted Trees on training data")
display(gbt_train_metrics)

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
transformed_validation_data = encoding_pipeline.transform(validation_data)['features', 'label']

gbt_validation_preds = gbt_model.transform(transformed_validation_data)
gbt_validation_metrics = evaluation_metrics(gbt_validation_preds, "Gradient-Boosted Trees on validation data")
display(gbt_validation_metrics)

# COMMAND ----------

# MAGIC %md #### Gradient-Boosted Trees with PR

# COMMAND ----------

# MAGIC %md Running this example on the limited set of highly correlated columns

# COMMAND ----------

join_statement = " JOIN pagerank_origin op ON op.ORIGIN_IATA = f.ORIGIN JOIN pagerank_dest dp ON dp.DEST_IATA = f.DEST"

train_data_pr = spark.sql("SELECT * FROM flights_and_weather_train_processed f" + join_statement).drop(*cols_to_drop)
validation_data_pr = spark.sql("SELECT * FROM flights_and_weather_validation_processed f" + join_statement).drop(*cols_to_drop)
test_data_pr = spark.sql("SELECT * FROM flights_and_weather_test_processed f" + join_statement).drop(*cols_to_drop)

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data_pr[correlated_cols],'DEP_DEL15')).fit(train_data_pr[correlated_cols])

# apply the transformations to our train data
transformed_train_data = encoding_pipeline.transform(train_data_pr[correlated_cols])['features', 'label']

# COMMAND ----------

startTime = time.time()
gbt = GBTClassifier(labelCol = "label", featuresCol = "features", maxIter = 10)
gbt_model = gbt.fit(transformed_train_data)
train_preds = gbt_model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Gradient-Boosted Tree model is: {(endTime - startTime) / (60)} minutes")
train_metrics = evaluation_metrics(train_preds, "Gradient-Boosted Trees on training data")
display(train_metrics)

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
transformed_validation_data = encoding_pipeline.transform(validation_data_pr[correlated_cols])['features', 'label']

#transformed_validation_data = encoding_pipeline.transform(validation_data)['features', 'label']
validation_preds = gbt_model.transform(transformed_validation_data)
validation_metrics = evaluation_metrics(validation_preds, "Gradient-Boosted Trees on validation data")
display(validation_metrics)

# COMMAND ----------

# COMMAND ----------

# MAGIC %md ## Random Forest

# COMMAND ----------

# MAGIC %md Note that we ran this model in another notebook and the output is not available here. Our results are summarized in the overview notebook.

# COMMAND ----------

# create an encoding pipeline based on information from our training data
# encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data,'DEP_DEL15')).fit(train_data)

# apply the transformations to our train data
# transformed_train_data = encoding_pipeline.transform(train_data)['features', 'label']

# train a model on our transformed train data
startTime = time.time()
rf = RandomForestClassifier(labelCol = "label", featuresCol = "features", numTrees = 45,
                           featureSubsetStrategy = "auto", impurity = "gini", subsamplingRate = 0.6, maxDepth = 25, maxBins = 16)
model = rf.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Random Forest model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Random Forest on training data")
display(train_metrics)

# COMMAND ----------

#transformed_validation_data = encoding_pipeline.transform(validation_data)['features', 'label']
validation_preds = model.transform(transformed_validation_data)
validation_metrics = evaluation_metrics(validation_preds, "Random Forest on validation data")
display(validation_metrics)

# COMMAND ----------

# MAGIC %md ## PageRank Random Forest

# COMMAND ----------

#encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data,'DEP_DEL15')).fit(train_data)

# apply the transformations to our train data
#transformed_train_data = encoding_pipeline.transform(train_data)['features', 'label']

# COMMAND ----------

# train a model on our transformed train data
startTime = time.time()
rf = RandomForestClassifier(labelCol = "label", featuresCol = "features", numTrees = 45,
                           featureSubsetStrategy = "auto", impurity = "gini", subsamplingRate = 0.6, maxDepth = 25, maxBins = 16)
model = rf.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Random Forest model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Random Forest on training data")
display(train_metrics)

# COMMAND ----------

# MAGIC %md ### Validation Evaluation

# COMMAND ----------

#transformed_validation_data = encoding_pipeline.transform(validation_data)['features', 'label']
validation_preds = model.transform(transformed_validation_data)
validation_metrics = evaluation_metrics(validation_preds, "Random Forest on validation data")
display(validation_metrics)

# COMMAND ----------

# MAGIC %md # 8. Model Selection and Tuning

# COMMAND ----------

# MAGIC %md **We selected the Logistic Regression model for further fine-tuning for reasons presented in the summary notebook.**

# COMMAND ----------

# MAGIC %md ## PageRank Logistic Regression Model

# COMMAND ----------

join_statement = " JOIN pagerank_origin op ON op.ORIGIN_IATA = f.ORIGIN JOIN pagerank_dest dp ON dp.DEST_IATA = f.DEST"

train_data_pr = spark.sql("SELECT * FROM flights_and_weather_train_processed f" + join_statement).drop(*cols_to_drop)
validation_data_pr = spark.sql("SELECT * FROM flights_and_weather_validation_processed f" + join_statement).drop(*cols_to_drop)
test_data_pr = spark.sql("SELECT * FROM flights_and_weather_test_processed f" + join_statement).drop(*cols_to_drop)

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data_pr,'DEP_DEL15')).fit(train_data_pr)

# apply the transformations to our train data
transformed_train_data = encoding_pipeline.transform(train_data_pr)['features', 'label']

# train a model on our transformed train data
startTime = time.time()
lr = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.001)
model = lr.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

# MAGIC %md ### Evaluate our training results

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Logistic Regression for PageRank Train")
display(train_metrics)

# COMMAND ----------

# MAGIC %md ### Evaluate our valuation results

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
transformed_validation_data = encoding_pipeline.transform(validation_data_pr)['features', 'label']

# run the fitted model on the transformed validation data
validation_preds = model.transform(transformed_validation_data)

# display our evaluation metrics
validation_metrics = evaluation_metrics(validation_preds, "Logistic Regression for PageRank Validation")
display(validation_metrics)

# COMMAND ----------

# MAGIC %md
# MAGIC ## Logistic Regression RA's Full 586 Columns

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data,'DEP_DEL15'))
encoding_pipeline = encoding_pipeline.fit(train_data)
# apply the transformations to our train data
transformed_train_data = encoding_pipeline.transform(train_data)['features', 'label']

# COMMAND ----------

# train a model on our transformed train data
startTime = time.time()
lr = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.001, elasticNetParam = 0.25, standardization = False)
model = lr.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

# MAGIC %md ### Let's evaluate our training results

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Logistic Regression on training data")
display(train_metrics)

# COMMAND ----------

AUC(model, train_preds)

# COMMAND ----------

ROC(model, train_preds)

# COMMAND ----------

# MAGIC %md ### Train and Evaluate on Validation

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
transformed_validation_data = encoding_pipeline.transform(validation_data)['features', 'label']

# run the fitted model on the transformed validation data
validation_preds = model.transform(transformed_validation_data)

# display our evaluation metrics
validation_metrics = evaluation_metrics(validation_preds, "Logistic Regression on validation data")
display(validation_metrics)

# COMMAND ----------

AUC(model, validation_preds)

# COMMAND ----------

ROC(model, validation_preds)

# COMMAND ----------

# MAGIC %md
# MAGIC ## Logistic Regression on Rolling Average Engineered features, PCA

# COMMAND ----------

pca = PCA(k=100, inputCol="features", outputCol="pca_features")
pca_model = pca.fit(transformed_train_data)
pca_transformed_train_data = pca_model.transform(transformed_train_data)

# COMMAND ----------

# train a model on our transformed train data
startTime = time.time()
lr_pca = LogisticRegression(featuresCol = 'pca_features', labelCol = 'label', maxIter = 30, regParam = 0.001, elasticNetParam = 0.25, standardization = False)
lr_pca_model = lr_pca.fit(pca_transformed_train_data)
lr_pca_train_preds = model.transform(pca_transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

# MAGIC %md ### Let's evaluate our training data

# COMMAND ----------

lr_pca_train_metrics = evaluation_metrics(lr_pca_train_preds, "Logistic Regression on training data")
display(lr_pca_train_metrics)

# COMMAND ----------

#transform validation data
pca_transformed_validation_data = pca_model.transform(transformed_validation_data)

# COMMAND ----------

# run the fitted model on the transformed validation data
lr_pca_validation_preds = lr_pca_model.transform(pca_transformed_validation_data)

# COMMAND ----------

# MAGIC %md ### Let's evaluate our validation data

# COMMAND ----------

# display our evaluation metrics
lr_pca_validation_metrics = evaluation_metrics(lr_pca_validation_preds, "Logistic Regression on validation data")
display(lr_pca_validation_metrics)

# COMMAND ----------

# MAGIC %md
# MAGIC ## Logistic Regression on Rolling Average Engineered features, with highly correlated columns

# COMMAND ----------

cor_picked_cols = ['ORIGIN_WND_SPEED_RATE', 'ORIGIN_WIND_GUST_OBSERVATION_SPEED_RATE', 'DEST_WND_SPEED_RATE', 'ORIGIN_WND_DIRECTION_ANGLE_3_RA', 'ORIGIN_WND_SPEED_RATE_3_RA',
 'ORIGIN_CIG_CEILING_HEIGHT_DIMENSION_3_RA', 'ORIGIN_RUNWAY_VISUAL_RANGE_OBSERVATION_VISIBILITY_DIMENSION_3_RA', 'ORIGIN_GD4_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION_3_RA',
 'ORIGIN_WIND_GUST_OBSERVATION_SPEED_RATE_3_RA', 'DEST_WND_SPEED_RATE_3_RA', 'DEST_AA1_LIQUID_PRECIPITATION_PERIOD_QUANTITY_3_RA', 'DEST_RUNWAY_VISUAL_RANGE_OBSERVATION_VISIBILITY_DIMENSION_3_RA',
 'DEST_GA3_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION_3_RA', 'DEST_GD1_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION_3_RA', 'DEST_WIND_GUST_OBSERVATION_SPEED_RATE_3_RA', 'DEP_DEL15']

# COMMAND ----------

ra_train_data_2 = spark.sql('SELECT * FROM flights_and_weather_train_ra_processed')[cor_picked_cols]
ra_validation_data_2 = spark.sql('SELECT * FROM flights_and_weather_validation_ra_processed')[cor_picked_cols]
ra_test_data_2 = spark.sql('SELECT * FROM flights_and_weather_test_ra_processed')[cor_picked_cols]

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(ra_train_data_2,'DEP_DEL15'))
encoding_pipeline = encoding_pipeline.fit(ra_train_data_2)
# apply the transformations to our train data
transformed_train_data_2 = encoding_pipeline.transform(ra_train_data_2)['features', 'label']

# COMMAND ----------

# train a model on our transformed train data
startTime = time.time()
lr_2 = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.001, elasticNetParam = 0.25, standardization = False)
lr_model2 = lr_2.fit(transformed_train_data_2)
lr_2_train_preds = lr_model2.transform(transformed_train_data_2)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

# MAGIC %md ### Let's evaluate our training data

# COMMAND ----------

lr_2_train_metrics = evaluation_metrics(lr_2_train_preds, "Logistic Regression on training data")
display(lr_2_train_metrics)

# COMMAND ----------

# MAGIC %md ### Let's evaluate our validation data

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
transformed_validation_data_2 = encoding_pipeline.transform(ra_validation_data_2)['features', 'label']

# COMMAND ----------

# run the fitted model on the transformed validation data
lr_2_validation_preds = lr_model2.transform(transformed_validation_data_2)

# COMMAND ----------

# display our evaluation metrics
lr_2_validation_metrics = evaluation_metrics(lr_2_validation_preds, "Logistic Regression on validation data")
display(lr_2_validation_metrics)

# COMMAND ----------

# COMMAND ----------

# MAGIC %md
# MAGIC ## Logistic Regression with arrival delay noted 2 hours before departure

# COMMAND ----------

cols_to_drop = ['ORIGIN', 'DEST', 'ORIGIN_IATA', 'DEST_IATA', 'DEP_DELAY', 'DEP_DELAY', 'DEP_DELAY_NEW', 'DEP_DELAY_GROUP', 'TAIL_NUM', 'ORIGIN_CITY_NAME', 'DEST_CITY_NAME', 'FL_DATE', 'ORIGIN_WEATHER_KEY', 'DEST_WEATHER_KEY', 'ORIGIN_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'ORIGIN_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'ORIGIN_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'ORIGIN_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'DEST_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'DEST_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'DEST_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION']

# COMMAND ----------

join_statement = " JOIN linked_arr_delays lad ON lad.ORIGIN_WEATHER_KEY = f.ORIGIN_WEATHER_KEY AND lad.TAIL_NUM = f.TAIL_NUM"

train_data_arr = spark.sql("SELECT * FROM flights_and_weather_train_processed version as of 3 f" + join_statement).drop(*cols_to_drop).drop(*['ORIGIN_WEATHER_KEY', 'TAIL_NUM'])
validation_data_arr = spark.sql("SELECT * FROM flights_and_weather_validation_processed version as of 2 f" + join_statement).drop(*cols_to_drop).drop(*['ORIGIN_WEATHER_KEY', 'TAIL_NUM'])
test_data_arr = spark.sql("SELECT * FROM flights_and_weather_test_processed version as of 2 f" + join_statement).drop(*cols_to_drop).drop(*['ORIGIN_WEATHER_KEY', 'TAIL_NUM'])

# COMMAND ----------

join_statement = " JOIN linked_arr_delays lad ON lad.ORIGIN_WEATHER_KEY = f.ORIGIN_WEATHER_KEY AND lad.TAIL_NUM = f.TAIL_NUM"

train_data_arr = spark.sql("SELECT * FROM flights_and_weather_train_ra_processed f" + join_statement).drop(*cols_to_drop).drop(*['ORIGIN_WEATHER_KEY', 'TAIL_NUM'])
validation_data_arr = spark.sql("SELECT * FROM flights_and_weather_validation_ra_processed f" + join_statement).drop(*cols_to_drop).drop(*['ORIGIN_WEATHER_KEY', 'TAIL_NUM'])
test_data_arr = spark.sql("SELECT * FROM flights_and_weather_test_ra_processed f" + join_statement).drop(*cols_to_drop).drop(*['ORIGIN_WEATHER_KEY', 'TAIL_NUM'])

# COMMAND ----------

dataset_size=float(train_data_arr.select("DEP_DEL15").count())
numPositives=train_data_arr.select("DEP_DEL15").where('DEP_DEL15 == "1"').count()
per_ones=(float(numPositives)/float(dataset_size))*100
numNegatives=float(dataset_size-numPositives)
print('The number of ones are {}'.format(numPositives))
print('Percentage of ones are {}'.format(per_ones))

BalancingRatio= numNegatives/dataset_size
print('BalancingRatio = {}'.format(BalancingRatio))

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data_arr,'DEP_DEL15')).fit(train_data_arr)

# apply the transformations to our train data
transformed_train_data = encoding_pipeline.transform(train_data_arr)['features', 'label']

# COMMAND ----------

BalancingRatio = 0.82
transformed_train_data_with_weights = transformed_train_data.withColumn("classWeights", f.when(transformed_train_data.label == "1",BalancingRatio).otherwise(1-BalancingRatio))

# COMMAND ----------

# Using results from our parameter search
startTime = time.time()
lr = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.01, elasticNetParam = 0.75, standardization = False, weightCol="classWeights")
model = lr.fit(transformed_train_data_with_weights)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Logistic Regression on training data")
display(train_metrics)

# COMMAND ----------

AUC(model, train_preds)
ROC(model, train_preds)

# COMMAND ----------

# apply the encoding transformations from our pipeline to the validation data
#transformed_validation_data = encoding_pipeline.transform(validation_data_arr)['features', 'label']

# run the fitted model on the transformed validation data
validation_preds = model.transform(transformed_validation_data)
# display our evaluation metrics
validation_metrics = evaluation_metrics(validation_preds, "Logistic Regression on validation data")
display(validation_metrics)

# COMMAND ----------

# MAGIC %md ## Running best model on Test

# COMMAND ----------

# MAGIC %md ## Logistic Regression RA's Full 586 Columns

# COMMAND ----------

train_data = spark.sql('SELECT * FROM flights_and_weather_train_ra_processed')
validation_data = spark.sql('SELECT * FROM flights_and_weather_validation_ra_processed')
test_data = spark.sql('SELECT * FROM flights_and_weather_test_ra_processed')

cols_to_drop = ['DEP_DELAY', 'DEP_DELAY_NEW', 'DEP_DELAY_GROUP', 'TAIL_NUM', 'ORIGIN_CITY_NAME', 'DEST_CITY_NAME', 'FL_DATE', 'ORIGIN_WEATHER_KEY', 'DEST_WEATHER_KEY', 'ORIGIN_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'ORIGIN_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'ORIGIN_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'ORIGIN_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'ORIGIN_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'ORIGIN_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'ORIGIN_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_AL2_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_PERIOD_QUANTITY', 'DEST_AL3_SNOW_ACCUMULATION_DEPTH_DIMENSION', 'DEST_GA6_SKY_COVER_LAYER_BASE_HEIGHT_DIMENSION', 'DEST_GD5_SKY_COVER_SUMMATION_STATE_HEIGHT_DIMENSION', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_UPPER_RANGE_ATTRIBUTE', 'DEST_SKY_CONDITION_OBSERVATION_BASE_HEIGHT_LOWER_RANGE_ATTRIBUTE', 'DEST_GG1_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG2_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG3_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION', 'DEST_GG4_BELOW_STATION_CLOUD_LAYER_TOP_HEIGHT_DIMENSION']

train_data = train_data.drop(*cols_to_drop)
validation_data = validation_data.drop(*cols_to_drop)
test_data = test_data.drop(*cols_to_drop)

# COMMAND ----------

# create an encoding pipeline based on information from our training data
encoding_pipeline = Pipeline(stages = create_encoding_stages(train_data,'DEP_DEL15'))
encoding_pipeline = encoding_pipeline.fit(train_data)
transformed_train_data = encoding_pipeline.transform(train_data)

# COMMAND ----------

# train a model on our transformed train data
startTime = time.time()
lr = LogisticRegression(featuresCol = 'features', labelCol = 'label', maxIter = 30, regParam = 0.001, elasticNetParam = 0.25, standardization = False)
model = lr.fit(transformed_train_data)
train_preds = model.transform(transformed_train_data)
endTime = time.time()
print(f"The training time of the Logistic Regression model is: {(endTime - startTime) / (60)} minutes")

# COMMAND ----------

train_metrics = evaluation_metrics(train_preds, "Logistic Regression on training data")
display(train_metrics)

# COMMAND ----------

# apply the encoding transformations from our pipeline to the test data
transformed_test_data = encoding_pipeline.transform(test_data)['features', 'label']

# run the fitted model on the transformed test data
test_preds = model.transform(transformed_test_data)

# display our evaluation metrics
test_metrics = evaluation_metrics(test_preds, "Logistic Regression on test data")
display(test_metrics)

# COMMAND ----------

AUC(model, test_preds)
ROC(model, test_preds)

# COMMAND ----------

# COMMAND ----------

reason: import numpy,import networkx


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Convert a tdms file to an fcs file using nptdms and fcswrite

(requires nptdms>=0.23.0)
"""
import nptdms
import numpy as np
from os.path import dirname, abspath, join

import fcswrite

def read_tdms(tdms_file):
    """Read tdms file and return channel names and data"""
    tdms_file = nptdms.TdmsFile(tdms_file)
    ch_names = []
    ch_data = []

    for ch in tdms_file.groups()[0].channels():
        if ch.data is not None and len(ch.data):
            chn = ch.name
            if "unit_string" in ch.properties:
                unit = ch.properties["unit_string"]
                ch_names.append("{} [{}]".format(chn, unit))
            else:
                ch_names.append(chn)

            ch_data.append(ch.data)

    return ch_names, ch_data

def add_deformation(chn_names, data):
    """From circularity, compute the deformation

    This method is useful for RT-DC data sets that contain
    the circularity but not the deformation.
    """
    if "deformation" not in chn_names:
        for ii, ch in enumerate(chn_names):
            if ch == "circularity":
                chn_names.append("deformation")
                data.append(1-data[ii])

    return chn_names, data

def tdms2fcs(tdms_file):
    """Creates an fcs file for a given tdms file"""
    fcs_file = tdms_file[:-4]+"fcs"
    chn_names, data = read_tdms(tdms_file)
    chn_names, data = add_deformation(chn_names, data)
    fcswrite.write_fcs(filename=fcs_file,
                       chn_names=chn_names,
                       data=np.array(data).transpose())

if __name__ == "__main__":
    # Create an fcs file where this tdms file is located:
    tdms_file = join(dirname(abspath(__file__)), "data/example_rt-dc.tdms")
    tdms2fcs(tdms_file)

max_stars_repo_name: ZELLMECHANIK-DRESDEN/fcswrite

max_issues_repo_name: ZELLMECHANIK-DRESDEN/fcswrite

max_forks_repo_name: ZELLMECHANIK-DRESDEN/fcswrite