text:
"""

Provides project-specific plotting functions.

The plotting library consists of one main plotting function `plot`, which makes use of several subfunctions
to combine/overlay plots of different data types:

- plot_fenics_function: for FENICS 2D data, with subfunctions

    - plot_fenics_function_scalar: for FENICS 2D scalar field data
    - plot_fenics_function_vector: for FENICS 2D vector field data

- plot_sitk_image: for 2D image data

Using this general plot function, multiple default layouts are defined for plotting the following data fields in 2D:

- plot_concentration
- plot_growth
- plot_proliferation
- plot_displacement

"""

import os
from builtins import classmethod
from datetime import time

import SimpleITK as sitk
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import tri
from glimslib.visualisation.helpers import mesh2triang, interpolate_over_grid, get_ranges_colormap, exclude_from_data, \
    add_colorbar

import glimslib.utils.file_utils as fu
from glimslib.visualisation import helpers as vh
from glimslib.fenics_local import Function, Point, cells, parameters

#==== PLOTTING SUBFUNCTIONS

#--- FENICS FUNCTIONS

def plot_fenics_function_vector(ax, f, mode='quiver', plot_nth=None,
                                interpolate=False, n_interpolate=100,
                                range_f=None, cmap='gist_earth', norm=None, norm_ref=None, n_cmap_levels=None,
                                exclude_below=None, exclude_above=None, exclude_min_max=False, exclude_around=None,
                                color=None, **kwargs):
    """
    Subfunction for plotting FENICS 2D vectorfields,
    returns plot axis.
    """
    mesh = f.function_space().mesh()
    gdim = mesh.geometry().dim()
    w0 = f.compute_vertex_values(mesh)
    nv = mesh.num_vertices()
    if len(w0) != gdim * nv:
        raise AttributeError('Vector length must match geometric dimension.')
    #-- rearrange coordinates
    X = mesh.coordinates()
    X_coords = [X[:, i] for i in range(gdim)]
    U = [w0[i * nv : (i + 1) * nv ] for i in range(gdim)]

    # -- quiver and streamplot assume that x and y axes are increasing.
    #    -> swap sign of Ux or Uy if one of the axes is ordered in different direction
    ax_x_lim = ax.get_xlim()
    ax_y_lim = ax.get_ylim()
    if ax_x_lim[0] > ax_x_lim[1]:
        print("-- fenics vector plot: x axis is decreasing; change sign of vector x value for plotting")
        U[0] = - U[0]
    if ax_y_lim[0] > ax_y_lim[1]:
        print("-- fenics vector plot: y axis is decreasing; change sign of vector y value for plotting")
        U[1] = - U[1]

    if mode=='quiver':
        # -- catch kwargs
        if 'density' in kwargs:
            kwargs.pop('density')

        if interpolate:
            X, Y, UX, UY = interpolate_over_grid(X_coords[0], X_coords[1], U[0], U[1],
                                                 n=n_interpolate,
                                                 return_coords='meshgrid', method='cubic')
            #plot = ax.quiver(xi, yi, ui, vi, [C], **kwargs)
        else:
            X, Y, UX, UY = [X_coords[0], X_coords[1], U[0], U[1]]

        if color: # ignore cmap
            if 'cbar_label' in kwargs:
                kwargs.pop('cbar_label')
            plot = ax.quiver(X, Y, UX, UY, color=color, **kwargs)
        else:
            # -- Compute magnitude
            C2 = UX ** 2 + UY ** 2
            C = np.sqrt(C2)

            # -- colormap and range settings
            min_f, max_f, colormap, norm = get_ranges_colormap(C,
                                                               range=range_f, cmap=cmap, norm=norm, norm_ref=norm_ref,
                                                               n_cmap_levels=n_cmap_levels)

            # -- exclude data from plotting
            mask = exclude_from_data(C, min_f, max_f,
                                     exclude_below=exclude_below, exclude_above=exclude_above,
                                     exclude_min_max=exclude_min_max, exclude_around=exclude_around,
                                     data_type='standard')
            C[mask] = np.nan
            plot = ax.quiver(X, Y, UX, UY, C, cmap=colormap, norm=norm, **kwargs)

    elif mode=='streamlines':
        xi, yi, ui, vi = interpolate_over_grid(X_coords[0], X_coords[1], U[0], U[1],
                                               n=n_interpolate,
                                               return_coords='linspace', method='cubic')
        plot = ax.streamplot(xi, yi, ui, vi, **kwargs)
    return plot

def plot_fenics_function_scalar(ax, f, showmesh=True,
                                range_f=None, cmap='gist_earth', norm=None, norm_ref=None, n_cmap_levels=None,
                                exclude_below=None, exclude_above=None, exclude_min_max=False, exclude_around=None,
                                plot_nth=None, shading='flat', alpha=1,
                                **kwargs):
    """
    Subfunction for plotting FENICS 2D scalar field,
    returns plot axis.
    """
    mesh    = f.function_space().mesh()
    nv      = mesh.num_vertices()
    gdim    = mesh.geometry().dim()
    mesh_coordinates = mesh.coordinates().reshape((nv, gdim))
    triangulation = mesh2triang(mesh)
    values = np.asarray([f(point) if len(mesh.bounding_box_tree().compute_collisions(Point(point)))>0 # check if point in domain
                                      else np.nan                                        # assign nan otherwise
                                for point in mesh_coordinates])
    #-- colormap and range settings
    min_f, max_f, colormap, norm = get_ranges_colormap(values,
                                                       range=range_f, cmap=cmap, norm=norm, norm_ref=norm_ref,
                                                       n_cmap_levels=n_cmap_levels)

    #-- exclude data from plotting
    data = values[triangulation.triangles]
    mask = exclude_from_data(data, min_f, max_f,
                             exclude_below=exclude_below, exclude_above=exclude_above,
                             exclude_min_max=exclude_min_max, exclude_around=exclude_around,
                             data_type='triangulation')
    triangulation.set_mask(mask)

    #-- plot
    if showmesh:
        plot = ax.tripcolor(triangulation, values, cmap=colormap, norm=norm, shading=shading, edgecolors='k', linewidth=0.05,
                          vmin=min_f, vmax=max_f, alpha=alpha, **kwargs)
    else:
        plot = ax.tripcolor(triangulation, values, cmap=colormap, norm=norm, shading=shading,
                           vmin=min_f, vmax=max_f, alpha=alpha, edgecolors="none", linewidth=0.0, **kwargs)
    return plot

def plot_fenics_function(ax, f, mode='quiver', **kwargs):
    """
    Subfunction for plotting FENICS 2D vector or scalarfield,
    returns plot axis.
    """
    mesh = f.function_space().mesh()
    gdim = mesh.geometry().dim()
    tdim = mesh.topology().dim()
    w0 = f.compute_vertex_values(mesh)
    #nv = mesh.num_vertices()
    value_dim = f.function_space().element().value_dimension(0)
    #if len(w0) != gdim * nv:
    #    raise AttributeError('Vector length must match geometric dimension.')

    if gdim == 2:
        # -- scalar
        if value_dim == 1:
            if 'color' in kwargs:
                kwargs.pop('color')
            plot = plot_fenics_function_scalar(ax, f, **kwargs)
            return plot
        elif value_dim == 2:
            if 'showmesh' in kwargs:
                kwargs.pop('showmesh')
            if 'shading' in kwargs:
                kwargs.pop('shading')
            plot = plot_fenics_function_vector(ax, f, mode=mode, **kwargs)
            return plot
        else:
            print("Plot function not defined for value dimension > 2. Here: %i" % tdim)

    else:
        print("Plot function not defined for geometric dimension %i" % gdim)

#--- SITK IMAGE
def plot_sitk_image(ax, image, segmentation=None, contour=False,
                    range=None, cmap='gist_earth', norm=None, norm_ref=0, n_cmap_levels=None,
                    label_alpha=1,
                    origin='lower',
                    **kwargs):
    """
    Subfunction for plotting image data (SimpleITK data format),
    returns plot axis.
    """
    #-- Convert Image Types
    img_type = sitk.sitkUInt8
    image_rec = sitk.Cast(sitk.RescaleIntensity(image), img_type)
    #-- Segmentation
    if segmentation:
        image_label_rec = sitk.Cast(segmentation, img_type)
        if contour:
            img_ol = sitk.LabelOverlay(image_rec, sitk.LabelContour(image_label_rec), opacity=label_alpha)
        else:
            img_ol = sitk.LabelOverlay(image_rec, image_label_rec, opacity=label_alpha)
    else:
        img_ol = image_rec
    #-- Prepare Figure with image
    nda = sitk.GetArrayFromImage(img_ol)
    min_f, max_f, colormap, norm = get_ranges_colormap(nda,
                                                       range=range, cmap=cmap, norm=norm, norm_ref=norm_ref,
                                                       n_cmap_levels=n_cmap_levels)

    img_origin  = img_ol.GetOrigin()
    img_spacing = img_ol.GetSpacing()
    img_size    = img_ol.GetSize()
    # extent -> (left, right, bottom, top)
    extent = (img_origin[0], img_origin[0] + img_size[0] * img_spacing[0],
              img_origin[1] + img_size[1] * img_spacing[1], img_origin[1])
    plot = ax.imshow(nda, interpolation=None, extent=extent, origin=origin,
                     cmap=colormap, norm=norm, vmin=min_f, vmax=max_f)
    return plot

#--- MAIN PLOT FUNCTION

def plot(plot_object_list, dpi=100, plot_range=None, margin=0.02, cbarwidth=0.05,
         save_path=None, show=True, xlabel='x position [mm]', ylabel='y position [mm]',
         show_axes=True, show_cbar=True, show_title=True, show_ticks=True,
         cbar_size='5%', cbar_pad=0.05, cbar_fontsize=None,
         **kwargs):
    """
    Each element in `plot_object_list` is a dictionary of the form::

        { 'object' : the object to be plotted,
          'param1' : one plot specific parameter,
          'param2' : another plot specific parameter,
          ...
        }

    See :py:meth:`show_img_seg_f()` for examples.
    Unless a 'zorder' argument is specified, elements are plotted in the order of occurrence in `plot_object_list`.

    """

    # -- create Figure
    fig, ax = plt.subplots(dpi=dpi)
    ax.set_aspect('equal')

    # -- if an image is provided, the axes will be oriented as in the image,
    #   otherwise use xlim, ylim for fixing display orientation
    # -- for discussion about imshow/extent https://matplotlib.org/tutorials/intermediate/imshow_extent.html
    if plot_range:
        ax.set_xlim(*plot_range[0:2])
        ax.set_ylim(*plot_range[2:4])

    if 'title' in kwargs:
        title = kwargs.pop('title')
        if (not title is None) and show_title:
            fig.suptitle(title)
    #-- Check if only a single plot_object has been provided, if so, transform to dict
    if not type(plot_object_list)==list:
        plot_object_dict = {}
        plot_object_dict['object'] = plot_object_list
        plot_object_dict.update(kwargs)
        plot_object_list = [plot_object_dict]
    # -- Iterate through plot objects
    cbar_ax_list = [ax]
    for plot_object_dict_orig in plot_object_list:
        plot_object_dict = plot_object_dict_orig.copy()
        plot_object = plot_object_dict.pop('object')
        # check for 'color'
        if 'color' in plot_object_dict:
            color = plot_object_dict['color']
        else:
            color=False
        cbar = False
        if ('cbar_label' in plot_object_dict) and not color:
            cbar_label = plot_object_dict.pop('cbar_label')
            cbar=True

        #-- use global kwargs as reference and overwrite with settings that are specific to this plot_object
        params = kwargs.copy()
        params.update(plot_object_dict)
        #-- plot
        if type(plot_object)==Function:
            plot = plot_fenics_function(ax, plot_object, **params)
        elif type(plot_object)==sitk.Image:
            plot = plot_sitk_image(ax, plot_object, **params)
        else:
            print("The plot_object is of type '%s' -- not supported."%(type(plot_object)))
            raise Exception

        if cbar and show_cbar:
            #cbax = add_colorbar(fig, cbar_ax_list[0], plot, cbar_label)
            cbax = add_colorbar(fig, cbar_ax_list[0], plot, cbar_label, size=cbar_size, pad=cbar_pad, fontsize=cbar_fontsize)
            cbar_ax_list.append(cbax)
            #fig.subplots_adjust(left=margin, bottom=margin, right=1 - margin - cbarwidth, top=1 - 2 * margin)

    if not show_axes:
        ax.spines['right'].set_visible(False)
        ax.spines['top'].set_visible(False)
        ax.spines['bottom'].set_visible(False)
        ax.spines['left'].set_visible(False)

    if show_ticks:
        if xlabel:
            ax.set_xlabel(xlabel)
        if ylabel:
            ax.set_ylabel(ylabel)
    else:
        ax.get_xaxis().set_ticks([])
        ax.get_yaxis().set_ticks([])
        ax.tick_params(bottom="off", left="off")

    if save_path:
        fu.ensure_dir_exists(os.path.dirname(save_path))
        fig.savefig(save_path, bbox_inches='tight',  dpi=dpi)
        print("  - Saved figure to '%s'"%save_path)

    if show:
        vh.show_plot()

def show_img_seg_f(image=None, segmentation=None, function=None, contour=False, alpha=1.0, dpi=100,
                   showmesh=True, shading='flat', colormap=None, n_cmap_levels=20, alpha_f=0.5, range_f=None,
                   exclude_below=None, exclude_above=None, exclude_min_max=False, exclude_around=None,
                   title=None, label=None, show=True, path=None, cmap_ref=None, norm=None, color=None,
                   plot_range=None,
                   show_axes=True, show_cbar=True, show_title=True, show_ticks=True,
                   cbar_size='5%', cbar_pad=0.05, cbar_fontsize=None,
                   add_plot_object_pre=None, add_plot_object_post=None):
    """
    Convenience function providing default settings for certain types of plots
    """
    plot_list = []
    if add_plot_object_pre is not None:
        plot_list.append(add_plot_object_pre)
    if image:
        plot_obj_img = {'object': image,
                        'segmentation': segmentation,
                        'label_alpha': alpha,
                        'contour': True,
                        'cmap': 'Greys',
                        'origin': 'upper'
                        }

        plot_list.append(plot_obj_img)

    if function:
        plot_obj_function = {'object': function,
                           'cbar_label': label,
                           'exclude_below': exclude_below,
                           'exclude_above': exclude_above,
                           'exclude_min_max': exclude_min_max,
                           'exclude_around': exclude_around,
                           'cmap': colormap,
                           'n_cmap_levels': n_cmap_levels,
                           'range_f': range_f,
                           'showmesh': showmesh,
                           'shading': shading,
                           'alpha': alpha_f,
                           'norm': norm,
                           'norm_ref': cmap_ref,
                           'color': color
                           }

        plot_list.append(plot_obj_function)
    if add_plot_object_post is not None:
        plot_list.append(add_plot_object_post)

    plot(plot_list, title=title, save_path=path, show=show, dpi=dpi, plot_range=plot_range,
         show_axes=show_axes, show_cbar=show_cbar, show_title=show_title, show_ticks=show_ticks,
         cbar_size=cbar_size, cbar_pad=cbar_pad, cbar_fontsize=cbar_fontsize)

def plot_concentration(image, label, fun, title, path=None, show=False, plot_range=None):
    """
    Convenience function for plotting concentration fields.
    """
    show_img_seg_f(image, label, fun, contour=True, showmesh=False, alpha_f=1,
                       range_f=[0.001, 1.01], exclude_min_max=True, colormap='viridis',
                       n_cmap_levels=20, title=title, label="concentration",
                       path=path, show=show, plot_range=plot_range)

def plot_growth(image, label, fun, title, path=None, show=False):
    """
    Convenience function for plotting growth fields.
    """
    show_img_seg_f(image, label, fun, contour=True, showmesh=False, alpha_f=1,
                       range_f=[0.0, 0.2], exclude_as_range=True, colormap='viridis',
                       n_cmap_levels=20,title=title, label="growth",
                       path=path, show=show)

def plot_proliferation(image, label, fun, title, path=None, show=False):
    """
    Convenience function for plotting proliferation fields.
    """
    show_img_seg_f(image, label, fun, contour=True, showmesh=False, alpha_f=1,
                       exclude_around=(0, 0.0001), range_f=[-0.02, 0.1],
                       title=title, label="proliferation term",
                       path=path, show=show, colormap='RdBu_r', n_cmap_levels=20, cmap_ref=0.0)

def plot_displacement(image, label, fun, title, path=None, show=False):
    """
    Convenience function for plotting displacement fields.
    """
    show_img_seg_f(image, label, fun, contour=True, showmesh=False, alpha_f=1,
                       exclude_around=None, range_f=[0.0, 20], exclude_min_max=True, exclude_below=0.5,
                       title=title, label="displacement", color=None,
                       path=path, show=show, colormap='viridis', n_cmap_levels=20)

def sitk_show(img, title=None, margin=0.05, dpi=40):
    nda = sitk.GetArrayFromImage(img)
    spacing = img.GetSpacing()
    figsize = (1 + margin) * nda.shape[0] / dpi, (1 + margin) * nda.shape[1] / dpi
    extent = (0, nda.shape[1] * spacing[1], nda.shape[0] * spacing[0], 0)
    fig = plt.figure(figsize=figsize, dpi=dpi)
    ax = fig.add_axes([margin, margin, 1 - 2 * margin, 1 - 2 * margin])

    plt.set_cmap("gray")
    ax.imshow(nda, extent=extent, interpolation=None)

    if title:
        plt.title(title)

    vh.show_plot()

def show_img_seg(image, segmentation, contour=False, alpha=1.0, dpi=100):
        # == Convert image types to uint8
        img_type = sitk.sitkUInt8
        image_label_rec = sitk.Cast(segmentation, img_type)
        image_rec = sitk.Cast(image, img_type)
        if contour:
            sitk_show(sitk.LabelOverlay(image_rec, sitk.LabelContour(image_label_rec), alpha), dpi=dpi)
        else:
            sitk_show(sitk.LabelOverlay(image_rec, image_label_rec, alpha), dpi=dpi)

def plotmat(function, title=''):
    p = plot(function)
    p.set_cmap("viridis")
    p.set_clim(function.vector().min(), function.vector().max())
    plt.colorbar(p)
    plt.title(title)
    vh.show_plot()

def plot_plt(field_var, showmesh=True,  shading='flat', contours=[],
             colormap=plt.cm.jet, norm=None, cmap_ref=None, range_f=None,alpha_f=1,
             title=None, path=None, show=True, dpi=300):
    mesh = field_var.function_space().mesh()
    n = mesh.num_vertices()
    d = mesh.geometry().dim()
    value_dim = field_var.function_space().num_sub_spaces()
    #-- only 2D
    if d == 2:
        # Create triangulation
        mesh_coordinates = mesh.coordinates().reshape((n, d))
        triangles = np.asarray([cell.entities(0) for cell in cells(mesh)])
        triangulation = tri.Triangulation(mesh_coordinates[:, 0],
                                                  mesh_coordinates[:, 1],
                                                  triangles)
        plt.figure()
        parameters['allow_extrapolation'] = True
        z = np.asarray([field_var(point) if mesh.bounding_box_tree().collides_entity(Point(point)) # check if point in domain
                                          else np.nan                                        # assign nan otherwise
                                        for point in mesh_coordinates])
        #z = np.asarray([field_var(point) for point in mesh_coordinates])
        if range_f is None:
            min_f = min(z)
            max_f = max(z)
        else:
            min_f = range_f[0]
            max_f = range_f[1]
        if type(colormap)==str:
            cmap = plt.cm.get_cmap(colormap)
        elif colormap==None:
            cmap = plt.cm.get_cmap('gist_earth')
        else:
            cmap = colormap

        if (norm is None) and (not cmap_ref is None):
            norm = vh.MidpointNormalize(midpoint=cmap_ref, vmin=min_f, vmax=max_f)

        #-- scalar function
        if value_dim == 1:

            if showmesh:
                 plt.tripcolor(triangulation, z, cmap=cmap, norm=norm, shading=shading, edgecolors='k',linewidth=0.1,
                                vmin = min_f, vmax = max_f, alpha = alpha_f)
            else:
                 plt.tripcolor(triangulation, z, cmap=cmap, norm=norm, shading=shading,
                                vmin = min_f, vmax = max_f, alpha = alpha_f)

        elif value_dim == 2:
            plot()

        plt.colorbar()
        if type(contours)==int:
            plt.tricontour(triangulation, z, contours, colors = 'k')
        elif len(contours)>0:
            plt.tricontour(triangulation, z, levels=contours, colors='k')
    elif d == 3:
        plot(field_var)

    if title is not None:
        plt.suptitle(title, y=1.08)
    if path is not None:
        fu.ensure_dir_exists(os.path.dirname(path))
        plt.savefig(path, dpi=dpi,bbox_inches='tight')
        print("-- saved figure to '%s'"%path)
    if show:
        vh.show_plot()

max_stars_repo_name: danielabler/glimslib

max_issues_repo_name: danielabler/glimslib

max_forks_repo_name: danielabler/glimslib


text:
from __future__ import division, print_function
import time
import numpy as np
import collections
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import BaggingClassifier, AdaBoostClassifier
from sklearn.metrics import confusion_matrix
import matplotlib.pyplot as plt

def visualize(digits):
    print(digits.keys())
    print(digits['data'].shape)
    print(digits['target'].shape)
    print(collections.Counter(digits['target']))
    print(digits['images'].shape)
    plt.figure(figsize=(3, 3))
    plt.imshow(digits['images'][-1], cmap=plt.cm.gray_r, interpolation='nearest')
    plt.savefig('figures/dataset_demo.png', dpi = 300)

def plot_curve(accuracy, fname):
    plt.figure(figsize = (8, 5))
    plt.plot(range(1, 21), accuracy)
    plt.xlabel('n_estimators')
    plt.ylabel('Accuracy / %')
    plt.tight_layout()
    plt.savefig('figures/%spng' % fname, dpi = 300)

def plot_confusion_matrix(conf_arr, fname):
    norm_conf = conf_arr / np.sum(conf_arr, axis = 0)

    fig = plt.figure()
    plt.clf()
    ax = fig.add_subplot(111)
    ax.set_aspect(1)
    res = ax.imshow(norm_conf, cmap=plt.cm.jet, interpolation='nearest')

    for x in range(10):
        for y in range(10):
            ax.annotate(str(conf_arr[x, y]), xy=(y, x), horizontalalignment='center', verticalalignment='center')

    plt.xticks(range(10), range(10))
    plt.yticks(range(10), range(10))
    plt.savefig('figures/%s.png' % fname, dpi = 300)

def singleLR(digits):
    LR = LogisticRegression()
    X_train, X_test, y_train, y_test = train_test_split(digits['data'], digits['target'], test_size = 0.2, random_state = 817)
    begin_time = time.time()
    LR.fit(X_train, y_train)
    print('Training takes %0.2f seconds.' % (time.time() - begin_time))
    begin_time = time.time()
    y_pred = LR.predict(X_test)
    print('Predicting takes %0.3f seconds.' % (time.time() - begin_time))
    conf_arr = confusion_matrix(y_test, y_pred, labels = range(10))
    print('Prediction accuracy: %0.2f' % (LR.score(X_test, y_test) * 100))
    plot_confusion_matrix(conf_arr, 'confusion_matrix')

def baggingLR(digits, n_estimators = 10):
    BLR = BaggingClassifier(base_estimator = LogisticRegression(), n_estimators = n_estimators, random_state = 817)
    X_train, X_test, y_train, y_test = train_test_split(digits['data'], digits['target'], test_size = 0.2, random_state = 817)
    begin_time = time.time()
    BLR.fit(X_train, y_train)
    print('Training takes %0.2f seconds.' % (time.time() - begin_time))
    begin_time = time.time()
    y_pred = BLR.predict(X_test)
    print('Predicting takes %0.3f seconds.' % (time.time() - begin_time))
    accuracy = BLR.score(X_test, y_test) * 100
    print('Prediction accuracy: %0.2f' %  accuracy)
    return accuracy
    conf_arr = confusion_matrix(y_test, y_pred, labels = range(10))
    plot_confusion_matrix(conf_arr, 'confusion_matrix_bagging')

def boostingLR(digits, n_estimators = 10):
    BLR = AdaBoostClassifier(base_estimator = LogisticRegression(), n_estimators = n_estimators, random_state = 817)
    X_train, X_test, y_train, y_test = train_test_split(digits['data'], digits['target'], test_size = 0.2, random_state = 817)
    begin_time = time.time()
    BLR.fit(X_train, y_train)
    print('Training takes %0.2f seconds.' % (time.time() - begin_time))
    begin_time = time.time()
    y_pred = BLR.predict(X_test)
    print('Predicting takes %0.3f seconds.' % (time.time() - begin_time))
    accuracy = BLR.score(X_test, y_test) * 100
    print('Prediction accuracy: %0.2f' %  accuracy)
    return accuracy
    conf_arr = confusion_matrix(y_test, y_pred, labels = range(10))
    plot_confusion_matrix(conf_arr, 'confusion_matrix_bagging')

np.random.seed(817)
digits = datasets.load_digits()
visualize(digits)
singleLR(digits)
accuracy = [baggingLR(digits, i) for i in range(1, 21)]
plot_curve(accuracy, 'bagging_accuracy')
accuracy = [boostingLR(digits, i) for i in range(1, 21)]
plot_curve(accuracy, 'boosting_accuracy')

ext: py

max_stars_repo_path: Ensembling Models on Digits Dataset/cross_validation.py

reason: import numpy


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

# ------------------------------------------------------------------------------
#
# Author: Gabriele Girelli
# Email: gigi.ga90@gmail.com
# Version: 1.0.0
# Description: assign bed rows to regions and collapse rows.
# ------------------------------------------------------------------------------

# DEPENDENCIES =================================================================

import argparse
import numpy as np
import os
import pandas as pd
import sys

# Loaded local bed-tools-gg python library
sys.path.append(os.path.dirname(os.path.realpath(__file__)) + '/../lib/')
import bed_lib as bd

# Change pandas default options
pd.options.mode.chained_assignment = None  # default='warn'

# PARAMETERS ===================================================================

# Add script description
parser = argparse.ArgumentParser(description = '''
Assign bed rows to Region of Interest (ROIs) and collapse them.
Every row in the bedfile is assigned to a ROI from the regfile.
Then, the file is collapsed to have a single row per ROI.
Rows can be collapsed in different ways: sum, max, min, median, mean, count.
The output is in bed format.
''')

# Add params
parser.add_argument('regfile', type = str, nargs = 1,
	help = 'Path to bedfile, containing regions to be assigned to.')
parser.add_argument('bedfile', type = str, nargs = 1,
	help = 'Path to bedfile, containing rows to be assigned.')

# Add flags
parser.add_argument('-c', '--collapse', type = str, nargs = 1,
	choices = ['min', 'mean', 'median', 'max', 'count', 'sum'],
	help = '''Collapse method. Default: sum''',
	default = ['sum'])
parser.add_argument('-u',
	action = 'store_const', const = True, default = False,
	help = 'Keep bedfile rows that do not match any region.')
parser.add_argument('-m',
	action = 'store_const', const = True, default = False,
	help = 'Assign to bedfile rows that partially match a region.')
parser.add_argument('-l',
	action = 'store_const', const = True, default = False,
	help = 'Assign to bedfile rows that include a region.')
parser.add_argument('-o', metavar = 'outfile', type = str, nargs = 1,
	default = [False],
	help = 'Output file (not a bed). Output to stdout if not specified.')
parser.add_argument('--float',
	action = 'store_const', dest = 'f',
	const = True, default = False,
	help = 'Value column as floats.')
parser.add_argument('--no-header',
	action = 'store_const', dest = 'header',
	const = True, default = False,
	help = 'Bed file has no header.')

# Parse arguments
args = parser.parse_args()

# Retrieve arguments
regfile = args.regfile[0]
bedfile = args.bedfile[0]
selected_collapse = args.collapse[0]
keep_unassigned_rows = args.u
keep_marginal_overlaps = args.m
keep_including = args.l
use_name = True
outfile = args.o[0]
floatValues = args.f
noHeader = args.header

# Default variables
bedcolnames = ['chr', 'start', 'end', 'name', 'score']

# RUN ==========================================================================

# Read regions file
rois = pd.read_csv(regfile, '\t', names = bedcolnames)

# Read bed file
if noHeader:
	bed = pd.read_csv(bedfile, '\t', names = bedcolnames)
else:
	bed = pd.read_csv(bedfile, '\t', names = bedcolnames, skiprows = [0])

# Assign rois to bed rows
rois = bd.assign_to_rois(rois, bed, keep_unassigned_rows,
	keep_marginal_overlaps, keep_including, use_name, selected_collapse,
	floatValues = floatValues)

# Output
if False == outfile:
	for i in range(rois.shape[0]):
		print('\t'.join(rois.iloc[i, :].astype('str').tolist()))
else:
	rois.to_csv(outfile, sep = '\t', header = False, index = False)

# END ==========================================================================

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

max_stars_repo_name: ggirelli/gg-bed-tools

max_issues_repo_name: ggirelli/gg-bed-tools

max_forks_repo_name: ggirelli/gg-bed-tools


text:
import numpy as np
import gym
import tensorflow as tf
import random
from collections import deque

class agent(object):
    def __init__(self,state_size=8,actions=4,gamma=0.99,replay=int(1e6)):
        '''
        double dueling deep Q learning network for lunar lander
        '''

        self.replay_memory = deque(maxlen=replay)
        self.burn_in = 10000
        self.decay_len = 5e4                            #number of steps to decay epsilon over
        self.epsilon = 1.0                              #start epsilon
        self.epsilon_min = 0.02                         #end epsilon
        self.epsilon_step = (1.0-self.epsilon_min)/self.decay_len
        self.total_steps = 0
        self.gamma = gamma

        #inputs
        self.state = tf.placeholder(dtype=tf.float32, shape=[None,state_size])
        self.action = tf.placeholder(dtype=tf.int32, shape=[None])
        self.target = tf.placeholder(dtype=tf.float32, shape=[None])
        batch_size = tf.reduce_sum(tf.ones_like(self.action))

        #fast network -- used to select next action during training
        with tf.variable_scope('fast'):

            #shared hidden layer
            self.hidden = tf.layers.dense(self.state,512,activation=tf.nn.elu,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())

            #calculate value of state
            self.value_dense = tf.layers.dense(self.hidden,256,activation=tf.nn.elu,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.value = tf.layers.dense(self.value_dense,1,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())

            #calculate relative advantage of each action
            self.adv_dense = tf.layers.dense(self.hidden,512,activation=tf.nn.elu,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.advantage = tf.layers.dense(self.adv_dense,n_actions,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())

            #predicted Q values of actions
            self.Q_values = tf.squeeze(self.value + tf.subtract(
                            self.advantage,tf.reduce_mean(self.advantage,axis=1,keepdims=True)))

        #slow network -- used to calculate q value of next action during training
        with tf.variable_scope('slow'):
            self.hidden_ = tf.layers.dense(self.state,512,activation=tf.nn.elu,trainable=False,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.value_dense_ = tf.layers.dense(self.hidden_,256,activation=tf.nn.elu,trainable=False,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.value_ = tf.layers.dense(self.value_dense_,1,trainable=False,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.adv_dense_ = tf.layers.dense(self.hidden_,512,activation=tf.nn.elu,trainable=False,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.advantage_ = tf.layers.dense(self.adv_dense_,n_actions,trainable=False,
                            kernel_initializer=tf.contrib.layers.xavier_initializer())
            self.Q_values_ = tf.squeeze(self.value_ + tf.subtract(
                            self.advantage_,tf.reduce_mean(self.advantage_,axis=1,keepdims=True)))

        #get Q values of actions taken
        self.Q_pred = tf.gather_nd(self.Q_values, tf.stack((tf.range(batch_size),self.action), axis=1))

        #loss and optimizer
        self.loss = tf.reduce_mean(tf.losses.huber_loss(self.target,self.Q_pred))
        self.fast_vars = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES,scope='fast')
        self.slow_vars = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES,scope='slow')
        self.optimizer = tf.train.AdamOptimizer(5e-5).minimize(self.loss,var_list=self.fast_vars)

        self.sess = tf.Session()
        self.sess.run(tf.global_variables_initializer())

    def add_memory(self,observation,action,reward,next_observation,done):
        '''
        add experience to replay memory
        '''

        self.replay_memory.append((observation, action, reward, next_observation,done))

    def update_slow(self):
        '''
        copy network weights from fast network to slow network
        '''

        update_slow_target_ops = []
        for slow_var,fast_var in zip(self.slow_vars,self.fast_vars):
            update_slow_target_op = slow_var.assign(fast_var)
            update_slow_target_ops.append(update_slow_target_op)
        update_slow_target_op = tf.group(*update_slow_target_ops)

        self.sess.run(update_slow_target_op)

    def train(self,batch_size=1024):
        '''
        train model using samples from replay memory
        '''

        #burn in
        if len(self.replay_memory) < self.burn_in:
            return 0,self.epsilon

        #update slow network every 100 steps
        if self.total_steps % 100 == 0:
            self.update_slow()

        #create training batch from replay memory
        batch = random.sample(self.replay_memory,batch_size)
        states = np.asarray([elem[0] for elem in batch])
        actions = np.asarray([elem[1] for elem in batch])
        rewards = np.asarray([elem[2] for elem in batch])
        next_states = np.asarray([elem[3] for elem in batch])
        terminal = np.asarray([elem[4] for elem in batch])

        #calculate target rewards
        d = {self.state:next_states}
        action_v,q_v = self.sess.run([self.Q_values,self.Q_values_],feed_dict=d)
        action_v = np.argmax(action_v,1)
        next_q = q_v[range(batch_size),action_v]
        targets = rewards + terminal * self.gamma * next_q

        #train
        d = {self.state:states,self.action:actions,self.target:targets}
        loss,_ = self.sess.run([self.loss,self.optimizer],feed_dict=d)

        #decay epsilon
        self.total_steps += 1
        if self.total_steps < self.decay_len:
            self.epsilon -= self.epsilon_step

        return loss,self.epsilon

    def predict(self,state):
        '''
        predict next action given an input state
        '''

        #explore
        if np.random.random() < self.epsilon:
            action = np.random.randint(n_actions)

        #exploit
        else:
            Q_pred = self.sess.run(self.Q_values,feed_dict={self.state:np.expand_dims(state,0)})
            action = np.argmax(np.squeeze(Q_pred))

        return action

#init environment
env = gym.make('LunarLander-v2')
state_dim = np.prod(np.array(env.observation_space.shape))
n_actions = env.action_space.n
dqn_agent = agent(state_size=state_dim,actions=n_actions)
viz = False

#run episodes
for ep in range(5000):
    total_reward = 0
    steps_in_ep = 0
    observation = env.reset()
    mean_loss = []

    #run each episode for 1000 steps
    for t in range(1000):

        #get next action
        action = dqn_agent.predict(observation)
        next_observation, reward, done, _info = env.step(action)
        total_reward += reward

        #visualize
        if viz:
            env.render('human')

        #normalize reward
        reward /= 10

        #update memory and train
        dqn_agent.add_memory(observation, action, reward, next_observation, 0.0 if done else 1.0)
        loss,epsilon = dqn_agent.train()
        mean_loss.append(loss)

        observation = next_observation
        steps_in_ep += 1

        #terminate episode if done
        if done:
            break

    #show episode results
    mean_loss = np.mean(mean_loss)
    print('Episode %2i, Reward: %6.2f, Steps: %4i, Loss: %5.4f, Next eps: %4.3f' %
         (ep,total_reward,steps_in_ep,mean_loss,epsilon))

    #enable vizualization once model gets over 200 score
    if total_reward >= 200 and ep > 500:
        viz = True

env.close()


text:
#!/usr/bin/python
import sys
import numpy as np
from nwb import nwb_file
from nwb import nwb_utils as ut

"""
Store extracellular ephys data

"""

OUTPUT_DIR = "../created_nwb_files/"
file_name = __file__[0:-3] + ".nwb"
########################################################################
# create a new NWB file
# several settings are specified when doing so. these can be supplied within
#   the NWB constructor or defined in a dict, as in in this example
settings = {}
settings["file_name"] = OUTPUT_DIR + file_name

# each file should have a descriptive globally unique identifier
#   that specifies the lab and this experiment session
# the function nwb.create_identifier() is recommended to use as it takes
#   the string and appends the present date and time
#- settings["identifier"] = nwb.create_identifier("extracellular spikes example")
settings["identifier"] = ut.create_identifier("extracellular spikes example")

# indicate that it's OK to overwrite exting file
settings["mode"] = "w"

# specify the start time of the experiment. all times in the NWB file
#   are relative to experiment start time
# if the start time is not specified the present time will be used
settings["start_time"] = "Sat Jul 04 2015 3:14:16"

# provide one or two sentences that describe the experiment and what
#   data is in the file
settings["description"] = "Test file demonstrating a simple extracellular ephys recording"

# create the NWB file object. this manages the file
print("Creating " + settings["file_name"])
f = nwb_file.open(**settings)

########################################################################
# create two electrical series, one with a single electrode and one with many
# then create a spike event series

# first describe the device(s).  Assume each probe is using a separate device
# names of devices are fictional
f.set_dataset("<device_X>", "Probe p0 device description", name="acme_model_23")
f.set_dataset("<device_X>", "Probe p1 device description", name="FooBar-X1")

# create the electrode map
# example simulated recording is made from two 2-electrode probes named
#   'p0' and 'p1'. we need to define the locations of the electrodes
#   relative to each probe, and the location of the probes
# electrode coordinates are in meters and their positions
#   are relative to each other. the location of the probe itself is
#   stored separately. using absolute coordinates here, if they are known,
#   is still OK
electrode_map = [[0, 0, 0], [0, 1.5e-6, 0], [0, 0, 0], [0, 3.0e-5, 0]]
electrode_group = [ "p0", "p0", "p1", "p1" ]

# make the group used to store extracellular ephys data.  This group
# will be created inside group /general
ex_ephys = f.make_group("extracellular_ephys")
ex_ephys.set_dataset("electrode_map", electrode_map)
ex_ephys.set_dataset("electrode_group", electrode_group)

# set electrode impedances and describe filtering
# impedances are stored as text in case there is a range
ex_ephys.set_dataset("impedance", [ "1e6", "1.1e6", "1.2e6", "1.3e6" ])
ex_ephys.set_dataset("filtering", "description of filtering")

# specify the description, location and device for each probe
p0 = ex_ephys.make_group("<electrode_group_X>", "p0")
p0.set_dataset("description", "Description of p0")
p0.set_dataset("location", "CA1, left hemisphere, stereotactic coordinates xx, yy")
p0.set_dataset("device", "acme_model_23")
p1 = ex_ephys.make_group("<electrode_group_X>", "p1")
p1.set_dataset("description", "Description of p1")
p1.set_dataset("location", "CA3, left hemisphere, stereotactic coordinates xx, yy")
p1.set_dataset("device", "FooBar-X1")

########################################################################
# the example is of two 2-electrode probes. the electrode data from these
#   probes can be stored individually, grouped as probes (eg, 2-electrode
#   pair) or all stored together. these approaches are all exampled here

# create time series with all electrode data stored together
quad = f.make_group("<ElectricalSeries>", "quad", path="/acquisition/timeseries")
quad.set_attr("comments", "Data corresponds to four electrodes (two probes)")
quad.set_dataset("data", np.zeros((10000, 4)), attrs={"resolution": 1.2345e-6})
quad.set_dataset("timestamps", np.arange(10000) * 0.0001)
# indicate that we're recording from the first electrode defined in the
#   above map (electrode numbers start at zero, so electrodes are
#   0, 1, 2 and 3
quad.set_dataset("electrode_idx", [0, 1, 2, 3])

########################################################################
# spikes can be reported by hardware or be detected by software
# in both cases, they are considered to be processed data and so belong
#   in a processing module

# create the module
spike_mod = f.make_group("<Module>", "my spikes")

# create an interface that stores the events. here they will be stored
#   with their waveforms, such as would be the input to a spike-sorting
#   algorithm
spike_iface = spike_mod.make_group("EventWaveform")
spike_iface.set_attr("source", "Data from device FooBar-X1 using dynamic multi-phasic threshold of 5xRMS")

# the event waveform interface publishes a SpikeEventSeries. make
#   that series
spike = spike_iface.make_group("<SpikeEventSeries>", "my waveforms")
spike.set_attr("comments", "Snapshots of spike events pulled from a recording")
spike.set_dataset("electrode_idx", [2, 3])    # probe 'p1'
# describe the source of the data (may be redundant w/ interface source)
spike.set_attr("source", "Data from device FooBar-X1 using dynamic multi-phasic threshold of 5xRMS")
# make some bogus simulated data
# this is 20 events all having the same shape and a pseudorandom time
evt = np.zeros((8,2))
evt[3][0] = 0.01
evt[4][0] = -0.005
evt[3][1] = 0.005
evt[4][1] = -0.0025
data = []
t = []
last = 1.0
for i in range(20):
    data.append(evt)
    last = last + (i * 17) % 29
    t.append(last)
# save the simulated data and time
spike.set_dataset("data", data, attrs={"resolution": 1.2345e-6})
spike.set_dataset("timestamps", t)

# if data were stored in another unit such as microvolts, it would be
#   necessary to specify a converstion between that unit and Volts.
#   that would be done by including the conversion when setting the "data"
#   dataset.  Call would be:
# spike.set_dataset("data", data, attrs={"resolution": 1.2345e-6, "conversion":1.0e-6})

# close file, otherwise it will fail to write properly
f.close()

max_stars_repo_name: bendichter/api-python

max_issues_repo_name: bendichter/api-python

max_forks_repo_name: bendichter/api-python


text:
from delicious_loader import load_dataset_hierarchical

import numpy as np
from sklearn.metrics import precision_recall_fscore_support, roc_auc_score

from keras.models import Model, model_from_json
from keras.layers import Dense, Input, Embedding, GlobalAveragePooling1D, TimeDistributed, LSTM, Dropout, Flatten
from keras.callbacks import EarlyStopping, ModelCheckpoint, Callback
from keras.layers.wrappers import Bidirectional
from keras import regularizers

import tensorflow as tf

def f1_score(y_true, y_pred):
	"""
	Compute the micro f(b) score with b=1.
	"""
	y_true = tf.cast(y_true, "float32")
	y_pred = tf.cast(tf.round(y_pred), "float32") # implicit 0.5 threshold via tf.round
	y_correct = y_true * y_pred

	sum_true = tf.reduce_sum(y_true, axis=1)
	sum_pred = tf.reduce_sum(y_pred, axis=1)
	sum_correct = tf.reduce_sum(y_correct, axis=1)

	precision = sum_correct / sum_pred
	recall = sum_correct / sum_true
	f_score = 2 * precision * recall / (precision + recall)
	f_score = tf.where(tf.is_nan(f_score), tf.zeros_like(f_score), f_score)

	return tf.reduce_mean(f_score)

def build_model(num_features,
	num_classes,
	embedding_dims,
	maxlen,
	max_sentence_len):
	"""
	"""

	input_layer = Input(shape=(maxlen,max_sentence_len,),
		dtype='int32')
	sentence_input = Input(shape=(max_sentence_len,),
		dtype='int32')

	embeddings = Embedding(num_features,
		embedding_dims,
		input_length=max_sentence_len,
		embeddings_regularizer=regularizers.l1(1e-6))(sentence_input)

	avg_layer = GlobalAveragePooling1D()(embeddings)
	sentEncoder = Model(inputs=sentence_input,
		outputs=avg_layer)
	sentEncoder.summary()
	textEncoder = TimeDistributed(sentEncoder)(input_layer)

	global_avg_layer = Flatten()(textEncoder)

	global_avg_layer = Dropout(0.5)(global_avg_layer)
	predictions = Dense(num_classes,
		activation='sigmoid',
		kernel_regularizer=regularizers.l1(1e-5))(global_avg_layer)

	model = Model(inputs=input_layer,
		outputs=predictions)
	model.compile(loss='binary_crossentropy',
              optimizer='adam',
              metrics=[f1_score])

	model.summary()

	return model

def load_model():
	"""
	"""

	json_file = open('hierarchical_model.json', 'r')
	loaded_model_json = json_file.read()
	json_file.close()
	model = model_from_json(loaded_model_json)

	model.load_weights('hierarchical_model.h5')
	print("Loaded model from disk")

	model.summary()

	model.compile(loss='binary_crossentropy',
		optimizer='adam',
		metrics=[f1_score])

	return model

if __name__ == '__main__':

	ngram_range = 1
	maxlen = 20
	max_sentence_len = 10
	batch_size = 32
	embedding_dims = 50
	epochs = 500
	num_classes = 20

	X_train,y_train,X_val,y_val,X_test,y_test,word_index = load_dataset_hierarchical(maxlen,max_sentence_len)

	num_features = len(word_index)
	print('Found %d words' % num_features)

	'''
	model = build_model(num_features,num_classes,embedding_dims,maxlen,max_sentence_len)

	model_json = model.to_json()
	with open("hierarchical_model.json", "w") as json_file:
		json_file.write(model_json)

	early_stopping =EarlyStopping(monitor='val_f1_score',
		patience=15,
		mode='max')
	bst_model_path = 'hierarchical_model.h5'
	model_checkpoint = ModelCheckpoint(bst_model_path,
		monitor='val_f1_score',
		verbose=1,
		save_best_only=True,
		mode='max',
		save_weights_only=True)

	model.fit(X_train, y_train,
		batch_size=batch_size,
		epochs=epochs,
		validation_data=(X_val, y_val),
		callbacks=[model_checkpoint,early_stopping])

	'''
	model = load_model()
	y_pred = model.predict(X_test)

	print 'AUC:',roc_auc_score(y_test, y_pred)
	y_pred[y_pred > 0.25] = 1
	y_pred[y_pred <= 0.25] = 0

	for i in range(10):
		pred,lab = y_pred[i],y_test[i]
		print np.where(pred == 1), np.where(lab == 1)

 	print precision_recall_fscore_support(y_test, y_pred, average='micro')
	print precision_recall_fscore_support(y_test, y_pred, average='macro')

max_stars_repo_name: AlexGidiotis/Advanced-ML-techniques

max_issues_repo_name: AlexGidiotis/Advanced-ML-techniques

max_forks_repo_name: AlexGidiotis/Advanced-ML-techniques


text:
import math
import pandas as pd
import numpy as np

class AllAnglePairs():
    def __init__(self, data_df, metadata_df):
        self.data_df = data_df
        self.metadata_df = metadata_df
        self.result_df = pd.DataFrame()

        # List of d-spacing pairs used to avoid duplicate plane comparisons
        self.pair_list = []
        self.angle_for_all()

    def angle_for_all(self):
        for i, rowi in self.data_df.iterrows():
            self.rowi_vs_all(rowi)

    def rowi_vs_all(self, rowi):
        for j, rowj in self.data_df.iterrows():
            if rowi['D-REF'] == rowj['D-REF']:
                pass
            else:
                self.rowi_vs_rowj(rowi, rowj)

    def rowi_vs_rowj(self, rowi, rowj):
        temp_df = pd.DataFrame()
        planei = [rowi['H'], rowi['K'], rowi['L']]
        planej = [rowj['H'], rowj['K'], rowj['L']]
        angle = angle_for_pair(planei, planej, self.metadata_df)
        data_dict = {'D-REF1': [rowi['D-REF']], 'D-REF2': [rowj['D-REF']],
                     'D-SPACING1': [rowi['D-SPACING']],
                     'D-SPACING2': [rowj['D-SPACING']], 'H1': [rowi['H']],
                     'K1': [rowi['K']], 'L1': [rowi['L']], 'H2': [rowj['H']],
                     'K2': [rowj['K']], 'L2': [rowj['L']], 'angle': [angle]}
        temp_df = pd.DataFrame(data=data_dict)
        pair_id = set([temp_df['D-SPACING1'].values[0],
                       temp_df['D-SPACING2'].values[0]])
        if self.result_df.shape[0] == 0:
            self.result_df = self.result_df.append(temp_df, ignore_index=True)
            self.pair_list.append(pair_id)
        elif pair_id not in self.pair_list:
            self.result_df = self.result_df.append(temp_df, ignore_index=True)
            self.pair_list.append(pair_id)

def angle_for_pair(plane1, plane2, metadata_df):
    """
    This is the function used to calculate the angle between
    two planes in miller indices,
    by calculating arccosine of dot product divide by length
    product plane1 in format of [h1, k1, l1],
    plane2 in format of [h2, k2, l2]. lattice_parameter should
    be written in format of [a, b, c] and alpha,
    beta, gamma are in degree
    """

    alpha = metadata_df['alpha']*math.pi/180
    beta = metadata_df['beta']*math.pi/180
    gamma = metadata_df['gamma']*math.pi/180
    V_square = ((
        metadata_df['a']*metadata_df['b']*metadata_df['c'])**2
    )*(1+2*math.cos(alpha)*math.cos(beta) *
       math.cos(gamma)-math.cos(alpha)**2-math.cos(beta)**2-math.cos(gamma)**2)
    temp = np.zeros((3, 3))
    temp[0][0] = (
        metadata_df['b']**2*metadata_df['c']**2*math.sin(alpha)**2
    )/V_square
    temp[0][1] = (
        metadata_df['a']*metadata_df['b']*metadata_df['c']**2*(
            math.cos(alpha)*math.cos(beta)-math.cos(gamma)
        ))/V_square
    temp[0][2] = (
        metadata_df['a']*metadata_df['b']**2*metadata_df['c']*(
            math.cos(alpha)*math.cos(gamma)-math.cos(beta)
        ))/V_square
    temp[1][0] = (
        metadata_df['a']*metadata_df['b']*metadata_df['c']**2*(
            math.cos(alpha)*math.cos(beta)-math.cos(gamma)
        ))/V_square
    temp[1][1] = (
        metadata_df['a']**2*metadata_df['c']**2*math.sin(beta)**2
    )/V_square
    temp[1][2] = (
        metadata_df['a']**2*metadata_df['b']*metadata_df['c']*(
            math.cos(beta)*math.cos(gamma)-math.cos(alpha)
        ))/V_square
    temp[2][0] = (
        metadata_df['a']*metadata_df['b']**2*metadata_df['c']*(
            math.cos(alpha)*math.cos(gamma)-math.cos(beta)
        ))/V_square
    temp[2][1] = (
        metadata_df['a']**2*metadata_df['b']*metadata_df['c']*(
            math.cos(beta)*math.cos(gamma)-math.cos(alpha)
        ))/V_square
    temp[2][2] = (
        metadata_df['a']**2*metadata_df['b']**2*math.sin(gamma)**2
    )/V_square
    d1_array = np.zeros(3)
    d2_array = np.zeros(3)
    for i in range(3):
        for j in range(3):
            d1_array[i] = d1_array[i]+plane1[j]*temp[j][i]
    for i in range(3):
        for j in range(3):
            d2_array[i] = d2_array[i]+plane2[j]*temp[j][i]
    d1_value = 0
    d2_value = 0
    for i in range(3):
        d1_value = d1_value+d1_array[i]*plane1[i]
    for i in range(3):
        d2_value = d2_value+d2_array[i]*plane2[i]

    numerator_array = np.zeros(3)
    numerator = 0
    for i in range(3):
        for j in range(3):
            numerator_array[i] = numerator_array[i]+plane1[j]*temp[j][i]
    for i in range(3):
        numerator = numerator+numerator_array[i]*plane2[i]
    denominator = math.sqrt(d1_value*d2_value)
    rad = math.acos(numerator/denominator)
    degree = 180*rad/math.pi
    return degree

max_stars_repo_name: crystalmaths/crystalmaths

max_issues_repo_name: crystalmaths/crystalmaths

max_forks_repo_name: crystalmaths/crystalmaths


text:
# python standard library
import datetime
from datetime import timedelta

# third party
import numpy

# this package
from theape import BaseClass, ApeError
from theape.infrastructure.strings import BOLD, RESET

DEBUG = 'debug'
INFO = 'info'
STAT_STRING = 'Time Stats -- Min: {min}, Q1: {q1}, Med: {med}, Q3: {q3}, Max: {max}, Mean: {mean}, StD: {std}'
ELAPSED_STRING = '{0}Elapsed Time:{1} {{0}}'.format(BOLD, RESET)
TOTAL_ELAPSED_STRING = '{0}Total Elapsed Time:{1} {{0}}'.format(BOLD, RESET)
ESTIMATED_REMAINING = '{0}Estimated Remaining Time:{1} {{0}}'.format(BOLD,
                                                                     RESET)
REPETITIONS_REMAINING = "{0}Repetitions Remaining:{1} {{0}}".format(BOLD, RESET)
TIME_LIMIT = "{0}Time-Limit Reached:{1} {{0}}".format(BOLD, RESET)
TIME_REMAINING = "{0}Time-Out Request in:{1} {{0}}".format(BOLD, RESET)
HARD_TIMEOUT = "{0}Absolute Quit at:{1} {{0}}".format(BOLD, RESET)
END_TIME = "{0}End-Time Reached:{1} {{0}}".format(BOLD, RESET)
MIN_PERCENTILE = 0
Q1_PERCENTILE = 25
MEDIAN_PERCENTILE = 50
Q3_PERCENTILE = 75
MAX_PERCENTILE = 100

CONTINUE = True
STOP = False
NOT_SET = None
UNSET = None
DECREMENT = -1
FINISHED = 0
ANNIHILATE = None

class TimeTracker(BaseClass):
    """
    A tracker of elapsed time
    """
    def __init__(self, log_level=DEBUG):
        """
        :param: `log_level`: level at which to report elapsed times (default='debug')
        """
        super(TimeTracker, self).__init__()
        self._logger = None
        self.log_level = log_level
        self.start = None
        self._times = None
        self._log = None
        return

    @property
    def times(self):
        """
        collection of elapsed times
        """
        if self._times is None:
            self._times = []
        return self._times

    @times.setter
    def times(self, times):
        """
        :param: ``times`` - collection
        :postcondition: self._times set to times
        """
        self._times = times
        return

    @property
    def log(self):
        """
        The logger method indicated by the log_level

        :return: logger.debug or logger.info
        """
        if self._log is None:
            if self.log_level == INFO:
                self._log = self.logger.info
            elif self.log_level == DEBUG:
                self._log = self.logger.debug
            else:
                raise ApeError("Unknown log level: {0}".format(self.log_level))
        return self._log

    def append(self, item):
        """
        Appends the item to the times array

        :param:

         - `item`: item to append to self.times (a numpy array)

        :postcondition: self.times contains item
        """
        self.times = numpy.append(self.times, [item])
        return

    def percentile(self, percentile):
        """
        calculates the percentile (e.g. 50 gets the median (the 50% item))

        :return: value for percintile of self.times as a timedelta
        """
        return timedelta(seconds=numpy.percentile(self.times, percentile))

    def log_update(self, elapsed):
        """
        Outputs to the log the most recent elapsed time information

        :param:

         - `elapsed`: timedelta
        """

        elapsed_string = ELAPSED_STRING.format(elapsed)
        self.log(elapsed_string)
        self.log(STAT_STRING.format(min=self.percentile(MIN_PERCENTILE),
                                    q1=self.percentile(Q1_PERCENTILE),
                                    med=self.percentile(MEDIAN_PERCENTILE),
                                    q3=self.percentile(Q3_PERCENTILE),
                                    max=self.percentile(MAX_PERCENTILE),
                                    mean=timedelta(seconds=numpy.mean(self.times)),
                                    std=timedelta(seconds=numpy.std(self.times))))
        return

    def __call__(self):
        """
        The main interface - starts and stops (toggles) the timer

        :return: True if starting, False if stopping
        :postcondition: elapsed time logged and added to self.times
        """
        if self.start is None:
            self.start = datetime.datetime.now()
            return True
        elapsed = datetime.datetime.now() - self.start
        # numpy can't handle timedeltas
        self.append(elapsed.total_seconds())
        self.start = None
        self.log_update(elapsed)
        return False

class CountdownTimer(TimeTracker):
    """
    A time-tracker that counts down
    """
    def __init__(self, repetitions=1, end_time=None, total_time=None,
                 *args, **kwargs):
        """
        :param:

         - ``repetitions``: number of calls to accept before stopping
         - ``end_time``: datetime to stop
         - ``total_time``: timedelta for amount of time to run
        """
        super(CountdownTimer, self).__init__(*args, **kwargs)
        self.repetitions = repetitions
        self.end_time = end_time
        self.total_time = total_time
        self.last_time = None
        return

    def time_remains(self):
        """
        Evaluates if there is still time (or repetitions) remaining

        :precondition: if total_time is set, self.start is set
        :return: True if reps or time remains, False otherwise
        """
        this_time = datetime.datetime.now()
        # check that the parameters have been set
        if not any((self.end_time, self.total_time, self.repetitions)):
            return STOP

        if self.end_time is not UNSET:
            if this_time >= self.end_time:
                # end-time takes first-precedence
                self.log(END_TIME.format(self.end_time))
                return STOP
            self.log(HARD_TIMEOUT.format(self.end_time))

        if self.total_time is not UNSET:
            # total_time is a relative time so it has to be on the LHS (need to fix this)
            if self.total_time <= (this_time - self.start):
                # total-time takes precedence over repetitions
                self.log(TIME_LIMIT.format(self.total_time))
                return STOP
            self.log(TIME_REMAINING.format(self.total_time - (this_time - self.start)))

        if self.repetitions is not UNSET:
            self.repetitions += DECREMENT
            if self.repetitions <= FINISHED:
                return STOP
            self.log(REPETITIONS_REMAINING.format(self.repetitions))

        # of no stopping-condtions are met, assume time-remains
        return CONTINUE

    def log_estimated_time_remaining(self):
        """
        Log an estitmated remaining time based on the median and repetitions  x xx
        """
        this_time = datetime.datetime.now()
        if not any((self.end_time, self.total_time, self.repetitions)):
            self.log(ESTIMATED_REMAINING.format(0))
            return

        estimated_end = estimated_total = estimated_reps = timedelta.max

        if self.end_time is not UNSET:
            estimated_end = self.end_time - this_time
        if self.total_time is not UNSET:
            estimated_total = self.total_time - (this_time - self.start)
        if self.repetitions is not UNSET:
            estimated_reps = (self.repetitions *
                              self.percentile(MEDIAN_PERCENTILE))

        self.log(ESTIMATED_REMAINING.format(min(estimated_end,
                                                estimated_total,
                                                estimated_reps)))
        return

    def __call__(self):
        """
        The main interface. Each call re-starts the timer and decrements the repetitions

        :return: True if repetitions > 0, False otherwise
        """
        call_time = datetime.datetime.now()

        if self.start is NOT_SET:
            self.start = self.last_time = call_time
            return CONTINUE

        # convert to seconds so numpy can calculate statistics
        elapsed, self.last_time = call_time - self.last_time, call_time

        self.append(elapsed.total_seconds())
        self.log_update(elapsed)

        if self.time_remains():
            self.log_estimated_time_remaining()
            return CONTINUE

        # out of time or repetitions, tear it down
        self.log(TOTAL_ELAPSED_STRING.format(call_time-self.start))
        self.close()
        return STOP

    def close(self):
        """
        Resets the attributes (makes __call__ always evaluate to false, not once like default)

        :postcondition:

         - ``self.end_time`` is None
         - ``self.total_time`` is None
         - ``self.repetitions`` is 0
         - ``self.start`` is None
         - ``self._times`` is None
        """
        self.start = UNSET
        self.end_time = UNSET
        self.total_time = UNSET
        self.repetitions = 0
        self._times = None
        return
# end class CountdownTimer

max_stars_repo_name: rsnakamura/theape

max_issues_repo_name: rsnakamura/theape

max_forks_repo_name: rsnakamura/theape


text:
from yo_fluq_ds import *
from unittest import TestCase
import os
import pandas as pd
import numpy as np

def path(*args):
    folder = FileIO.find_root_folder('yo.root')
    dir = os.path.join(folder, 'yo_fluq_ds__tests', *args)
    return dir


text:
import cmocean
import healpy as hp
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import AxesGrid

def StokesMatrix(n):
    if n not in [0, 1, 2, 3]:
        raise Exception("Input must be an integer in [0,1,2,3]")

    if n == 0:
        p = np.array([[1.0, 0], [0.0, 1.0]])
    elif n == 1:
        p = np.array([[1.0, 0], [0, -1.0]])
    elif n == 2:
        p = np.array([[0, 1.0], [1.0, 0]])
    elif n == 3:
        p = np.array([[0.0, -1j], [1j, 0]])

    return p

def MuellerMatrixElement(J, i, j):

    Pi = StokesMatrix(i)
    Pj = StokesMatrix(j)

    M_ij = np.einsum("...ab,...bc,...cd,...ad", Pi, J, Pj, J.conj()) / 2.0

    M_ij = np.real(M_ij)

    return M_ij

def PlotMueller(jones, rot_angle=90):
    npix = jones.shape[0]
    nside = hp.npix2nside(npix)
    xsize = 1600
    reso = 120 * 180 * np.sqrt(2.0) / np.pi / (xsize - 1)
    LambProj = hp.projector.AzimuthalProj(
        xsize=xsize, reso=reso, lamb=True, half_sky=True, rot=[0, rot_angle]
    )

    # Mind the negative sign on the z-coordinate. This came from trial and error for
    # making a particular plot. Its possible it should be removed.
    def p2v(x, y, z):
        return hp.vec2pix(nside, x, y, -z)

    logthresh = 4
    linscale = 2
    fig = plt.figure(figsize=(12, 12))
    grid = AxesGrid(
        fig,
        (1, 1, 1),
        nrows_ncols=(4, 4),
        axes_pad=(1.0, 0.5),
        label_mode="all",
        share_all=False,
        cbar_location="right",
        cbar_mode="each",
        cbar_size="5%",
        cbar_pad="1%",
    )
    for i in range(4):
        for j in range(4):
            M_ij = MuellerMatrixElement(jones, i, j)
            img_d = LambProj.projmap(M_ij, p2v)

            if i == j == 0:
                cmap = "viridis"
                vmin = 0
                vmax = 1

                tick_locs = list(np.linspace(0, 1, 7, endpoint=True))
                tick_labels = [
                    r"$ < 10^{-6}$",
                    r"$10^{-5}$",
                    r"$10^{-4}$",
                    r"$10^{-3}$",
                    r"$10^{-2}$",
                    r"$10^{-1}$",
                    r"$10^{0}$",
                ]

            elif i != j:
                cmap = "RdBu_r"
                vmin = -0.05
                vmax = 0.05

                d = np.log10(5) * np.diff(np.linspace(vmax * 1e-6, vmax, 7))[0]
                q = np.linspace(vmax * 1e-6, vmax, 7)[0::2] - d
                tick_locs = list(np.r_[-np.flipud(q)[:-1], [0], q[1:]])
                tick_labels = [
                    r"$-10^{-2}$",
                    r"$-10^{-4}$",
                    r"$-10^{-6}$",
                    r"$< 5 \times 10^{-8}$",
                    r"$10^{-6}$",
                    r"$10^{-4}$",
                    r"$10^{-2}$",
                ]

            else:
                # cmap='RdBu_r'
                cmap = cmocean.cm.delta
                vmin = -1.0
                vmax = 1

                q = np.linspace(vmax * 1e-6, vmax, 7)[0::2]
                tick_locs = list(np.r_[-np.flipud(q)[:-1], [0], q[1:]])
                tick_labels = [
                    r"$-10^{0}$",
                    r"$-10^{-2}$",
                    r"$-10^{-4}$",
                    r"$< 10^{-6}$",
                    r"$10^{-4}$",
                    r"$10^{-2}$",
                    r"$10^{0}$",
                ]
            n = 4 * i + j
            im = grid[n].imshow(
                img_d,
                interpolation="none",
                cmap=cmap,
                aspect="equal",
                vmin=vmin,
                vmax=vmax,
            )

            grid[n].set_xticks([])
            grid[n].set_yticks([])

            grid.cbar_axes[n].colorbar(im, ticks=tick_locs)
            grid.cbar_axes[n].set_yticklabels(tick_labels)

            im.set_norm(
                matplotlib.colors.SymLogNorm(
                    10 ** -logthresh, linscale, vmin=vmin, vmax=vmax
                )
            )
    plt.tight_layout(w_pad=0.5, h_pad=1.0)
    plt.show()


text:
"""Analysis for counting regions.
"""
import numpy as np
import svgwrite
from experiments.experiment import Experiment
from pysyrenn import LinesClassifier

class LinearityHypothesisExperiment(Experiment):
    """Empirically investigating Linear Explanation of Adversarial Examples.

    Corresponds to Section 5 of [1]. A number of different results are taken
    from this experiment:

    - Expository plot (Figure 4)
    - FGSM density/random density for {normal,diffai,pgd}x{mnist,cifar10}
      (Tables 3ab).
    - Normal/{DiffAI,PGD}x{mnist,cifar10} densities for {FGSM,random}
      directions (Tables 4ab).
    - Weighted relative gradient change (lines 240 - 250).
    """
    @staticmethod
    def perturb_l_inf(image, signs, l_inf):
        """Maximally perturbs (un-processed) @image within some l_inf bound.
        """
        distances = np.floor(l_inf * image)
        perturbed = image + (distances * signs)
        return perturbed.clip(0, 255)

    def fgsm(self, network, image, to_label, l_inf, process):
        """Returns an FGSM-perturbed version of @image.

        @process should be the pre-processing function for @network, while
        @to_label is the target label. @l_inf is the l_inf norm we respect.
        """
        signs = np.sign(network.compute_gradients(
            [process(image)], to_label))[0]
        return self.perturb_l_inf(image, signs, l_inf)

    def random_perturbation(self, image, l_inf, seed):
        """Returns a randomly-perturbed version of @image.

        @seed should be based on the image itself so that all networks get the
        same random images.
        """
        np.random.seed(seed)
        signs = np.sign(np.random.choice([-1, 1], image.shape))
        return self.perturb_l_inf(image, signs, l_inf)

    @staticmethod
    def gradient_error(network, line, transformed):
        """Returns the mean error in the gradient assumption made by FGSM.

        Effectively, this computes the relative error between the gradient
        computed at @line["start"] and the gradient within each linear region
        on @line. The Linearity Hypothesis works if these are all similar, i.e.
        if this function returns a small value.
        """
        start = line["start"]
        end = line["end"]
        delta = end - start
        label = line["target_fgsm_label"]

        start_gradient = network.compute_gradients([start], label)[0]

        region_lengths = []
        midpoint_preimages = []
        for start_ratio, end_ratio in zip(transformed[:-1], transformed[1:]):
            midpoint_ratio = (start_ratio + end_ratio) / 2.0
            midpoint_preimages.append(
                line["start"] + (midpoint_ratio * delta))
            region_lengths.append(end_ratio - start_ratio)

        along_gradients = network.compute_gradients(midpoint_preimages, label)
        gradient_errors = [np.mean(np.abs((start_gradient - along_gradient) /
                                          (along_gradient + 1e-12)))
                           for along_gradient in along_gradients]
        return np.average(gradient_errors, weights=region_lengths)

    def run_for_network(self, network_name, dataset):
        """Runs experiments for a particular network.

        @dataset should be the result of calling Experiment.load_input_data.
        """
        l_inf = 0.03
        network = self.load_network(network_name)
        process = dataset["process"]
        inputs = dataset["raw_inputs"]
        labels = dataset["labels"]

        # First, we generate line descriptions for each line we want to pass
        # through this network.
        lines = []
        for i, image in enumerate(inputs):
            label = labels[i]

            actual_label = np.argmax(network.compute([process(image)])[0])
            if actual_label != label:
                # We only consider points which are correctly classified by the
                # network.
                continue

            fgsm_label = (label + 1) % 10
            targets = {
                "fgsm": self.fgsm(network, image, fgsm_label, l_inf, process),
            }
            for j in range(4):
                targets["random%d" % j] = self.random_perturbation(
                    image, l_inf, ((20 * i) + j))

            image_dir = "%s/%03d" % (network_name, i)

            # Save the original image.
            self.record_artifact(self.rgbify_image(image),
                                 "%s/start" % image_dir, "rgb_image")

            # Then make lines to & save all of the other target images.
            for target_name, target_image in targets.items():
                target_image_key = "%s/%s" % (image_dir, target_name)
                self.record_artifact(self.rgbify_image(image),
                                     target_image_key, "rgb_image")

                target_label = np.argmax(
                    network.compute([process(target_image)])[0])
                length = np.linalg.norm(process(target_image) - process(image))

                lines.append({
                    "image": i,
                    "artifact_dir": image_dir,
                    "start": process(image),
                    "end": process(target_image),
                    "type": target_name,
                    "start_label": label,
                    "end_label": target_label,
                    "target_fgsm_label": fgsm_label,
                    "length": length,
                    "gradient_error": None,
                })

        transformed_lines = network.exactlines(
            [(line["start"], line["end"]) for line in lines],
            # We want post to get the 'break' points in the analysis.
            compute_preimages=False, include_post=True)

        for line, transformed in zip(lines, transformed_lines):
            # transformed is just a Numpy array of endpoint ratios because we
            # passed compute_preimages=False and include_post=False
            artifact_key = "%s/%s_transformed" % (
                line["artifact_dir"], line["type"])
            self.record_artifact(transformed, artifact_key, "pickle")
            # We use the gradients from the FGSM direction for the inline
            # mention of the "extended experiment."
            if line["type"] == "fgsm":
                line["gradient_error"] = self.gradient_error(network, line,
                                                             transformed[0])

        line_data_key = "%s/line_data" % network_name
        line_data_out = self.begin_csv(
            line_data_key,
            ["image", "artifact_dir", "type", "start_label", "end_label",
             "target_fgsm_label", "gradient_error", "length"],
            extrasaction="ignore")
        for line in lines:
            self.write_csv(line_data_out, line)
        self.record_artifact(line_data_key, line_data_key, "csv")

    def run(self):
        """Runs the linearity-hypothesis experiment.

        In the run phase, we select a number of lines (some shared by all
        models and some model-specific), then transform them and record the
        relevant statistics.
        """
        dataset_types = ["cifar10", "mnist"]
        network_types = ["relu_convsmall",
                         "relu_convsmall_diffai",
                         "relu_convsmall_pgd"]
        n_networks = input("Number of models (# or * for all): ")
        n_networks = -1 if n_networks == "*" else int(n_networks)
        networks_finished = 0
        network_names = []
        for dataset_type in dataset_types:
            dataset_name = "%s_test" % dataset_type
            for network_type in network_types:
                network_name = "%s_%s" % (dataset_type, network_type)
                network_names.append(network_name)
                print("Running Experiment for Network:", network_name)
                dataset = self.load_input_data(dataset_name,
                                               "conv" in network_name)
                self.run_for_network(network_name, dataset)
                networks_finished += 1
                if networks_finished == n_networks:
                    break
            if networks_finished == n_networks:
                break
        self.record_artifact(network_names, "networks", "pickle")

    @staticmethod
    def group_by_image(lines):
        """Groups lines by their starting ("natural") image.
        """
        groups = {}
        for line in lines:
            try:
                groups[int(line["image"])].append(line)
            except KeyError:
                groups[int(line["image"])] = [line]
        return groups

    def compute_density(self, line):
        """Computes the density for a line.
        """
        endpoints, postimages = self.read_artifact(
            "{artifact_dir}/{type}_transformed".format(**line))
        return len(endpoints) / float(line["length"])

    def classification_regions(self, line):
        """Computes the classification regions for a line.
        """
        transformed = self.read_artifact(
            "{artifact_dir}/{type}_transformed".format(**line))
        classifier = LinesClassifier.from_exactlines([transformed])
        return classifier.compute()[0]

    @staticmethod
    def corresponding_line(line, other_lines):
        """Finds a line in @other_lines that corresponds to @line.

        This is used for, eg., finding the same random-line in another model's
        results.
        """
        try:
            return next(other_line for other_line in other_lines
                        if (other_line["type"] == line["type"] and
                            other_line["image"] == line["image"]))
        except StopIteration:
            return None

    def analyze(self):
        """Analyzes the result of the experiment.
        """
        networks = self.read_artifact("networks")

        # First, we get the FGSM/Random data for all networks.
        print("Table 3 Data (FGSM/Random):")
        for network in networks:
            print("Network:", network)
            lines = self.read_artifact("%s/line_data" % network)
            groups = self.group_by_image(lines)
            ratios = []
            for image_lines in groups.values():
                # There should only be one FGSM line.
                fgsm_line = next(line for line in image_lines
                                 if line["type"] == "fgsm")
                fgsm_density = self.compute_density(fgsm_line)
                # There may be many random lines.
                random_lines = [line for line in image_lines
                                if line["type"].startswith("random")]
                random_densities = map(self.compute_density, random_lines)
                ratios.extend(fgsm_density / random_density
                              for random_density in random_densities)
            print("\t%s" % self.summarize(ratios))
        print("")

        print("Table 4 Data (Normal/{DiffAI,PGD})")
        normal_networks = [network for network in networks
                           if not ("pgd" in network or "diffai" in network)]
        for normal_network in normal_networks:
            print("Normal Network:", normal_network)
            normal_lines = self.read_artifact("%s/line_data" % normal_network)
            normal_fgsm_lines = [line for line in normal_lines
                                 if line["type"] == "fgsm"]
            normal_random_lines = [line for line in normal_lines
                                   if line["type"].startswith("random")]
            for suffix in ["diffai", "pgd"]:
                other_network = "%s_%s" % (normal_network, suffix)
                if other_network not in networks:
                    continue
                other_lines = self.read_artifact("%s/line_data" % other_network)
                other_fgsm_lines = [line for line in other_lines
                                    if line["type"] == "fgsm"]
                other_random_lines = [line for line in other_lines
                                      if line["type"].startswith("random")]
                ratios = []
                for normal_line in normal_fgsm_lines:
                    other_line = self.corresponding_line(normal_line,
                                                         other_fgsm_lines)
                    if other_line is None:
                        continue
                    ratios.append(self.compute_density(normal_line) /
                                  self.compute_density(other_line))
                print("\tFGSM Normal/%s:" % suffix, self.summarize(ratios))
                ratios = []
                for normal_line in normal_random_lines:
                    other_line = self.corresponding_line(normal_line,
                                                         other_random_lines)
                    if other_line is None:
                        continue
                    ratios.append(self.compute_density(normal_line) /
                                  self.compute_density(other_line))
                print("\tRandom Normal/%s:" % suffix, self.summarize(ratios))

            # Reproduce the gradient-error percentage.
            gradient_errors = [float(line["gradient_error"])
                               for line in normal_fgsm_lines]
            print("\tGradient error:", self.summarize(gradient_errors))

        print("Generating versions of Figure 4...")
        self.figure_4(networks)
        return True

    def figure_4(self, networks):
        """Builds expository figures for all valid lines.
        """
        if not ("cifar10_relu_convsmall" in networks and
                "cifar10_relu_convsmall_diffai" in networks):
            print("CIFAR10 convsmall {normal,diffai} networks not run, " +
                  "can't make Figure 4")
            return
        normal_net = "cifar10_relu_convsmall"
        diffai_net = "cifar10_relu_convsmall_diffai"
        normal_lines = self.read_artifact("%s/line_data" % normal_net)
        diffai_lines = self.read_artifact("%s/line_data" % diffai_net)
        normal_images = set(int(line["image"]) for line in normal_lines)
        diffai_images = set(int(line["image"]) for line in diffai_lines)
        both_images = normal_images & diffai_images
        for image in both_images:
            normal_fgsm = next(line for line in normal_lines
                               if (int(line["image"]) == image and
                                   line["type"] == "fgsm"))
            normal_random = next(line for line in normal_lines
                                 if (int(line["image"]) == image and
                                     line["type"] == "random3"))
            diffai_random = self.corresponding_line(normal_random, diffai_lines)
            fig4a, fig4b = self.figure_4_for_lines(normal_fgsm, normal_random,
                                                   diffai_random)
            self.record_artifact(fig4a, "figure_4a_%03d" % image, "svg")
            self.record_artifact(fig4b, "figure_4b_%03d" % image, "svg")

    def figure_4_for_lines(self, normal_fgsm, normal_random, diffai_random):
        """Builds the expository figure for a particular image.
        """
        fig4a = svgwrite.Drawing(profile="full")
        fig4b = svgwrite.Drawing(profile="full")

        def add_line(fig, line, left, top, width, height, color,
                     break_color="blue"):
            endpoints, _ = self.read_artifact(
                "{artifact_dir}/{type}_transformed".format(**line))
            mid = top + (height / 2)
            fig.add(fig.line((left, mid), (left + width, mid),
                             stroke=color, stroke_width=0.5))
            for endpoint in endpoints[1:-1]:
                x = left + (endpoint * width)
                fig.add(fig.line((x, top), (x, top + height),
                                 stroke=color, stroke_width=0.5))
            breaks = self.classification_regions(line)[0][1:]
            for break_distance in breaks:
                x = left + (break_distance[0] * width)
                fig.add(fig.line((x, top - 8), (x, top + height + 8),
                                 stroke=break_color, stroke_width=1.0))

        line_width = 125
        vpad = 3

        def add_start_image(fig, line, y):
            start_key = "{artifact_dir}/start".format(**line)
            start_image = self.read_artifact(start_key)
            start_datauri = self.image_to_datauri(start_image)
            fig.add(fig.image(start_datauri, (0, y), size=(32, 32)))

        def add_end_image(fig, line, y, labels, label_colors):
            end_key = "{artifact_dir}/{type}".format(**line)
            end_image = self.read_artifact(end_key)
            end_datauri = self.image_to_datauri(end_image)
            left = 32 + (2 * vpad) + line_width
            fig.add(fig.image(end_datauri, (left, y), size=(32, 32)))
            for i, label in enumerate(labels):
                fig.add(fig.text(label, insert=(left + 32 + 3, y + (i * 6) + 5),
                                 font_size="5", fill=label_colors[i]))

        labels = ["Airplane", "Automobile", "Bird", "Cat", "Deer", "Dog",
                  "Frog", "Horse", "Ship", "Truck"]
        # 4a top line is Normal-FGSM.
        add_start_image(fig4a, normal_fgsm, 16)
        fig4a.add(fig4a.text(labels[int(normal_fgsm["start_label"])],
                             insert=(12, 16 + 38), font_size="5"))
        add_end_image(fig4a, normal_fgsm, 0, [
            "Training Method: Normal",
            "Direction: FGSM +3%",
            "Density: %0.2f" % self.compute_density(normal_fgsm),
            "Label: %s" % labels[int(normal_fgsm["end_label"])],
        ], ["black", "black", "black", "black"])
        add_line(fig4a, normal_fgsm, 32 + vpad, 9, line_width, 14, "black")
        # 4a bottom line is Normal-Random
        add_end_image(fig4a, normal_random, 32 + vpad, [
            "Training Method: Normal",
            "Direction: Random +3%",
            "Density: %0.2f" % self.compute_density(normal_random),
            "Label: %s" % labels[int(normal_random["end_label"])],
        ], ["black", "black", "black", "black"])
        add_line(fig4a, normal_random, 32 + vpad, 35 + 9, line_width, 14,
                 "black")

        # 4b top line is Normal-Random.
        add_start_image(fig4b, normal_fgsm, 16)
        fig4b.add(fig4b.text(labels[int(normal_fgsm["start_label"])],
                             insert=(12, 16 + 38), font_size="5"))
        add_end_image(fig4b, normal_fgsm, 0, [
            "Training Method: Normal",
            "Direction: Random +3%",
            "Density: %0.2f" % self.compute_density(normal_random),
            "Label: %s" % labels[int(normal_random["end_label"])],
        ], ["black", "black", "black", "black"])
        add_line(fig4b, normal_random, 32 + vpad, 9, line_width, 14, "black")
        # 4b bottom line is DiffAI-Random
        add_end_image(fig4b, diffai_random, 32 + vpad, [
            "Training Method: DiffAI",
            "Direction: Random +3%",
            "Density: %0.2f" % self.compute_density(diffai_random),
            "Label: %s" % labels[int(diffai_random["end_label"])],
        ], ["green", "green", "green", "green"])
        add_line(fig4b, diffai_random, 32 + vpad, 35 + 9, line_width, 14,
                 "green")
        return fig4a, fig4b

if __name__ == "__main__":
    LinearityHypothesisExperiment("linearity_hypothesis").main()