text:
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import warnings
import seaborn as sns
from ipywidgets import widgets, interact
import rpy2.robjects as ro
from rpy2.robjects import numpy2ri
import copy
import cobra

from ipywidgets_New.interact import StaticInteract
from ipywidgets_New.widgets import RangeWidget, DropDownWidget

from Helper_methods import (
                             getSystemDistributionPerGraphLevel,
                             getSystemDistributionAcrossGraphLevels,
                             extractMetabolicSystems,
                             getListFrequencies,
                             getSystemFrequencies,
                             )
plt.style.use('seaborn')
ro.r['source']('Rfunctions.R')
FVA = ro.r['FVA']
sampleFluxCone = ro.r['sampleFluxCone']
getFluxOrders = ro.r['get_flux_orders']
numpy2ri.activate()

def plotDataOrder(data, color='#ffbf00', interactive=True, has_protein_costs=True, img_folder_name=None):
    if interactive and has_protein_costs:
        plot_func = plotDataOrderInteractiveWithCostsAndAlphas
    elif interactive and not has_protein_costs:
        plot_func = plotDataOrderInteractiveWithoutCostsAndAlphas
    else:
        plot_func = plotDataOrderStatic

    return plot_func(data, color=color, img_folder_name=img_folder_name)

def plotDataOrderStatic(plot_data, color='#25a9f4', img_folder_name=None):

    distributionSize = plot_data['distribution_size']
    cumSumDistributions = plot_data['cumulative']
    countsInOrderedPairs = plot_data['ordered']
    countsInSample = plot_data['sample']
    areaOrdered = plot_data['areaAboveZeroOfOrderedPairs']
    areaSample = plot_data['areaAboveZeroOfSample']
    pvalue = plot_data['pvalue']
    binEdges = plot_data['binEdges']
    dataName = plot_data['dataName']

    is_flux_data = 'flux data' in dataName.lower()
    if is_flux_data:
        delta_str = r'$\hat{{\delta}}$'
    else:
        delta_str = r'$\delta}$'

    fig = plt.figure(figsize=(18.5, 6.5))
    suptitle = r'{data}, area({delta} > 0) = {area:.2f}'.format(
        data=dataName,
        delta=delta_str,
        area=areaOrdered)
    plt.suptitle(suptitle, fontsize=16)

    plt.subplot(121)
    title = '$p$-value = {:.4f}'.format(pvalue)
    plotCumulativeDistributions(cumSumDistributions,
                                binEdges, color, title, xlabel=delta_str)
    plt.subplot(122)
    counts = {'Ordered': countsInOrderedPairs / countsInOrderedPairs.sum(),
              'Random': countsInSample / countsInSample.sum()}
    title = 'Size = {}'.format(distributionSize)
    colors = [color, '#afafaf']
    labels = [r'{}, area({} > 0) = {:.2f}'.format(type, delta_str, area)
              for type, area in zip(list(counts.keys()), [areaOrdered, areaSample])]

    plotHistograms(counts, binEdges, colors=colors, labels=labels,
                   title=title, xlabel=delta_str)
    plt.subplots_adjust(wspace=0.3, hspace=None)
    return fig

def plotDataOrderInteractiveWithoutCosts(plot_data, color='#ffbf00', img_folder_name=None):
    """
    Plot an interactive figure (of class StaticInteract) containig the distributions of
    data differences for the ordered pairs and the permutation sample. It takes as input
    plot_data: the dictionary returned by Orders.evaluateFluxOrders() and color: a str
    containing the color code to be used when plotting the ordered data distribution.
    """
    plt.rcParams.update({'figure.max_open_warning': 0})
    if img_folder_name is None:
        img_folder_name = plot_data[list(plot_data.keys())[0]]['dataName'].replace('(glucose)', '')

    def plot_func(carbonSource):
        warnings.filterwarnings('ignore')

        distributionSize = plot_data[carbonSource]['distribution_size']
        cumSumDistributions = plot_data[carbonSource]['cumulative']
        countsInOrderedPairs = plot_data[carbonSource]['ordered']
        countsInSample = plot_data[carbonSource]['sample']
        areaOrdered = plot_data[carbonSource]['areaAboveZeroOfOrderedPairs']
        areaSample = plot_data[carbonSource]['areaAboveZeroOfSample']
        pvalue = plot_data[carbonSource]['pvalue']
        binEdges = plot_data[carbonSource]['binEdges']
        dataName = plot_data[carbonSource]['dataName']

        fig = plt.figure(figsize=(20, 7))
        suptitle = r'{data}, area($\delta$ > 0) = {area:.2f}'.format(
            data=dataName,
            area=areaOrdered)
        plt.suptitle(suptitle, fontsize=16)

        plt.subplot(121)
        title = '$p$-value = {:.4f}'.format(pvalue)
        plotCumulativeDistributions(cumSumDistributions,
                                    binEdges, color, title, xlabel=r'$\delta$')
        plt.subplot(122)
        counts = {'Ordered': countsInOrderedPairs / countsInOrderedPairs.sum(),
                  'Random': countsInSample / countsInSample.sum()}
        title = 'Size = {}'.format(distributionSize)
        colors = [color, '#afafaf']
        labels = [r'{}, area($\delta$ > 0) = {:.2f}'.format(type, area)
                  for type, area in zip(list(counts.keys()), [areaOrdered, areaSample])]

        plotHistograms(counts, binEdges, colors=colors, labels=labels, title=title,
                       xlabel=r'$\delta$')
        plt.subplots_adjust(wspace=0.3, hspace=None)
        warnings.resetwarnings()
        return fig

    carbonSources = list(plot_data.keys())
    interactive_fig = StaticInteract(
        plot_func,
        carbonSource=DropDownWidget(carbonSources, description='Carbon source'),
        interact_name=img_folder_name
        )
    return interactive_fig

def plotDataOrderInteractiveWithoutCostsAndAlphas(plot_data, color='#ffbf00', img_folder_name=None):
    """
    Plot an interactive figure (of class StaticInteract) containig the distributions of
    data differences for the ordered pairs and the permutation sample. It takes as input
    plot_data: the dictionary returned by Orders.evaluateFluxOrders() and color: a str
    containing the color code to be used when plotting the ordered data distribution.
    """
    plt.rcParams.update({'figure.max_open_warning': 0})
    if img_folder_name is None:
        first_level_key_1 = list(plot_data.keys())[0]
        second_level_key_1 = list(plot_data[first_level_key_1].keys())[0]
        img_folder_name = plot_data[first_level_key_1][second_level_key_1]['dataName'].replace('(glucose)', '')

    def plot_func(a_carbonSource, b_alpha):
        warnings.filterwarnings('ignore')
        carbonSource = a_carbonSource
        alpha = b_alpha
        distributionSize = plot_data[carbonSource][alpha]['distribution_size']
        cumSumDistributions = plot_data[carbonSource][alpha]['cumulative']
        countsInOrderedPairs = plot_data[carbonSource][alpha]['ordered']
        countsInSample = plot_data[carbonSource][alpha]['sample']
        areaOrdered = plot_data[carbonSource][alpha]['areaAboveZeroOfOrderedPairs']
        areaSample = plot_data[carbonSource][alpha]['areaAboveZeroOfSample']
        pvalue = plot_data[carbonSource][alpha]['pvalue']
        binEdges = plot_data[carbonSource][alpha]['binEdges']
        dataName = plot_data[carbonSource][alpha]['dataName']

        fig = plt.figure(figsize=(20, 7))
        suptitle = r'{data}, area($\delta$ > 0) = {area:.2f}, alpha = {alpha}'.format(
            data=dataName,
            area=areaOrdered,
            alpha=alpha)
        plt.suptitle(suptitle, fontsize=16)

        plt.subplot(121)
        title = '$p$-value = {:.4f}'.format(pvalue)
        plotCumulativeDistributions(cumSumDistributions,
                                    binEdges, color, title, xlabel=r'$\delta$')
        plt.subplot(122)
        counts = {'Ordered': countsInOrderedPairs / countsInOrderedPairs.sum(),
                  'Random': countsInSample / countsInSample.sum()}
        title = 'Size = {}'.format(distributionSize)
        colors = [color, '#afafaf']
        labels = [r'{}, area($\delta$ > 0) = {:.2f}'.format(type, area)
                  for type, area in zip(list(counts.keys()), [areaOrdered, areaSample])]

        plotHistograms(counts, binEdges, colors=colors, labels=labels, title=title,
                       xlabel=r'$\delta$')
        plt.subplots_adjust(wspace=0.3, hspace=None)
        warnings.resetwarnings()
        return fig

    carbonSources = list(plot_data.keys())
    alphas = list(plot_data['glucose'].keys())
    interactive_fig = StaticInteract(
        plot_func,
        a_carbonSource=DropDownWidget(carbonSources, description='Carbon source', default='glucose'),
        b_alpha=DropDownWidget(alphas, description='alpha (biomass)', default=0.95),
        interact_name=img_folder_name
        )
    return interactive_fig

def plotDataOrderInteractiveWithCosts(plot_data, color='#ffbf00', img_folder_name=None):
    """
    Plot an interactive figure (of class StaticInteract) containig the distributions of
    data differences for the ordered pairs and the permutation sample. It takes as input
    plot_data: the dictionary returned by Orders.evaluateFluxOrders() and color: a str
    containing the color code to be used when plotting the ordered data distribution.
    """
    plt.rcParams.update({'figure.max_open_warning': 0})
    if img_folder_name is None:
        img_folder_name = plot_data['glucose'][0]['dataName'].replace('(glucose)', '')

    def plot_func(carbonSource, percentile):
        warnings.filterwarnings('ignore')
        distributionSize = plot_data[carbonSource][percentile]['distribution_size']
        cumSumDistributions = plot_data[carbonSource][percentile]['cumulative']
        countsInOrderedPairs = plot_data[carbonSource][percentile]['ordered']
        countsInSample = plot_data[carbonSource][percentile]['sample']
        areaOrdered = plot_data[carbonSource][percentile]['areaAboveZeroOfOrderedPairs']
        areaSample = plot_data[carbonSource][percentile]['areaAboveZeroOfSample']
        pvalue = plot_data[carbonSource][percentile]['pvalue']
        binEdges = plot_data[carbonSource][percentile]['binEdges']
        dataName = plot_data[carbonSource][percentile]['dataName']
        Plabel = '$P_{{}}$'.format(percentile)

        fig = plt.figure(figsize=(20, 7))
        suptitle = r'{data}, {plabel}, area($\hat{{\delta}}$ > 0) = {area:.2f}'.format(
            data=dataName,
            plabel=Plabel,
            area=areaOrdered)
        plt.suptitle(suptitle, fontsize=16)

        plt.subplot(121)
        title = '$p$-value = {:.4f}'.format(pvalue)
        plotCumulativeDistributions(cumSumDistributions,
                                    binEdges, color, title)
        plt.subplot(122)
        counts = {'Ordered': countsInOrderedPairs / countsInOrderedPairs.sum(),
                  'Random': countsInSample / countsInSample.sum()}
        title = 'Size = {}'.format(distributionSize)
        colors = [color, '#afafaf']
        labels = [r'{}, area($\hat{{\delta}}$ > 0) = {:.2f}'.format(type, area)
                  for type, area in zip(list(counts.keys()), [areaOrdered, areaSample])]

        plotHistograms(counts, binEdges, colors=colors, labels=labels, title=title)
        plt.subplots_adjust(wspace=0.3, hspace=None)
        warnings.resetwarnings()
        return fig

    carbonSources = list(plot_data.keys())
    percentiles = list(plot_data['glucose'].keys())
    PLabels = ['P' + str(p) for p in percentiles]
    interactive_fig = StaticInteract(
        plot_func,
        carbonSource=DropDownWidget(carbonSources, description='Carbon source', default='glucose'),
        percentile=DropDownWidget(percentiles, description='Percentile', labels=PLabels, default=0),
        interact_name=img_folder_name
        )
    return interactive_fig

def plotDataOrderInteractiveWithCostsAndAlphas(plot_data, color='#ffbf00', img_folder_name=None):
    """
    Plot an interactive figure (of class StaticInteract) containig the distributions of
    data differences for the ordered pairs and the permutation sample. It takes as input
    plot_data: the dictionary returned by Orders.evaluateFluxOrders() and color: a str
    containing the color code to be used when plotting the ordered data distribution.
    """
    plt.rcParams.update({'figure.max_open_warning': 0})
    if img_folder_name is None:
        img_folder_name = plot_data['glucose'][0][0]['dataName'].replace('(glucose)', '')

    def plot_func(a_carbonSource, b_percentile, c_alpha):
        warnings.filterwarnings('ignore')
        carbonSource = a_carbonSource
        percentile = b_percentile
        alpha = c_alpha

        if alpha in plot_data[carbonSource].keys():
            distributionSize = plot_data[carbonSource][alpha][percentile]['distribution_size']
            cumSumDistributions = plot_data[carbonSource][alpha][percentile]['cumulative']
            countsInOrderedPairs = plot_data[carbonSource][alpha][percentile]['ordered']
            countsInSample = plot_data[carbonSource][alpha][percentile]['sample']
            areaOrdered = plot_data[carbonSource][alpha][percentile]['areaAboveZeroOfOrderedPairs']
            areaSample = plot_data[carbonSource][alpha][percentile]['areaAboveZeroOfSample']
            pvalue = plot_data[carbonSource][alpha][percentile]['pvalue']
            binEdges = plot_data[carbonSource][alpha][percentile]['binEdges']
            dataName = plot_data[carbonSource][alpha][percentile]['dataName']
            Plabel = '$P_{{}}$'.format(percentile)

            fig = plt.figure(figsize=(20, 7))
            suptitle = r'{data}, {plabel}, area($\hat{{\delta}}$ > 0) = {area:.2f}, alpha = {alpha}'.format(
                data=dataName,
                plabel=Plabel,
                area=areaOrdered,
                alpha=alpha)
            plt.suptitle(suptitle, fontsize=16)

            plt.subplot(121)
            title = '$p$-value = {:.4f}'.format(pvalue)
            plotCumulativeDistributions(cumSumDistributions,
                                        binEdges, color, title)
            plt.subplot(122)
            counts = {'Ordered': countsInOrderedPairs / countsInOrderedPairs.sum(),
                      'Random': countsInSample / countsInSample.sum()}
            title = 'Size = {}'.format(distributionSize)
            colors = [color, '#afafaf']
            labels = [r'{}, area($\hat{{\delta}}$ > 0) = {:.2f}'.format(type, area)
                      for type, area in zip(list(counts.keys()), [areaOrdered, areaSample])]

            plotHistograms(counts, binEdges, colors=colors, labels=labels, title=title)
            plt.subplots_adjust(wspace=0.3, hspace=None)
        else:
            fig = plt.figure(figsize=(20, 7))
            plt.title(f'No data available for these flux-ordered pairs', fontsize=18)
            plt.plot(np.zeros(30), np.zeros(30), 'w-')

        warnings.resetwarnings()
        return fig

    carbonSources = list(plot_data.keys())
    percentiles = list(plot_data['glucose'][0].keys())
    PLabels = ['P' + str(p) for p in percentiles]
    alpha_values = list(plot_data['glucose'].keys())
    interactive_fig = StaticInteract(
        plot_func,
        a_carbonSource=DropDownWidget(carbonSources, description='Carbon source', default='glucose'),
        b_percentile=DropDownWidget(percentiles, description='Percentile', labels=PLabels, default=0),
        c_alpha=DropDownWidget(alpha_values, description='alpha (biomass)', default=0.95),
        interact_name=img_folder_name
        )
    return interactive_fig

def plotCumulativeDistributions(distributions=None, binEdges=None,
                                color='#ffbf00', title=None, xlabel=r'$\hat{\delta}$'):
    sampleDist = distributions['cumSumOfOrderedPairs']
    minDist = distributions['minCumSumsOfSamples']
    maxDist = distributions['maxCumSumsOfSamples']
    meanDist = distributions['meanCumSumsOfSamples']

    if title is not None:
        plt.title(title, fontsize=17)
    plt.xlabel(xlabel, fontsize=14)
    plt.ylabel('$area$', fontsize=14)
    plt.plot(binEdges[1:], sampleDist, color=color, linewidth=2.5)
    plt.plot(binEdges[1:], meanDist, color='black', linewidth=1.5)
    plt.fill_between(binEdges[1:], minDist, maxDist, color='#dddddd')

def plotHistograms(counts, bins, colors=None, labels=None,
                   title=None, xlabel=r'$\hat{\delta}$'):
    """
    Plot histograms from counts and bins, counts can be an array of
    arrays of counts, bins must be a single array of bin edges
    (same bins for all distributions)
    """
    if colors is None:
        colors = [None for _ in len(counts)]
    if labels is None:
        labels = [None for _ in len(counts)]
    if title is not None:
        plt.title(title, fontsize=17)
    plt.xlabel(xlabel, fontsize=14)
    plt.ylabel('$frequency$', fontsize=14)

    centroids = (bins[1:] + bins[:-1]) / 2
    names = list(counts.keys())
    for datatype, color, label in zip(names, colors, labels):

        sns.distplot(centroids, bins=bins, hist=True, kde=False,
                     kde_kws={'shade': True,
                              'linewidth': 3},
                     hist_kws={'weights': counts[datatype],
                               'range': (min(bins), max(bins)),
                               'edgecolor': color},
                     label=label, color=color)

    plt.legend(fontsize=12)

def plotValuesPerGraphLevel(values, labels=None, title=None):
    """
    Plots the distributions of number of regulatory events per hierarchy
    level
    """
    plt.figure(figsize=(12, 8))
    plt.boxplot(values, meanline=False, showmeans=True)
    plt.ylabel('values')
    if labels is not None:
        plt.xticks(list(range(1, len(values) + 1)), labels, fontsize=16)
    if title is not None:
        plt.title(title, fontsize=18, pad=20)
    locs, _ = plt.xticks()
    ax = plt.gca()
    sizes = [len(v) for v in values]
    for loc, size in zip(locs, sizes):
        plt.text(loc, 0.9, str(size),
                 transform=ax.get_xaxis_transform(),
                 fontsize=12)
    plt.show()

def plotLevelSubsystems(Model, graph_reactions, fig_dpi=100):
    """
    Computes the frequencies of the subsystems of the reactions appearing
    in the specified graph level which have the specified macrosystem
    """
    GEM = Model.GEM
    macrosystems = ['Amino acid metabolism', 'Carbohydrate metabolism', 'Cell wall biosynthesis',
                    'Cofactor and vitamin metabolism', 'Energy and maintenance', 'Lipid metabolism',
                    'Nucleotide metabolism', 'Transport']

    macrosystem_frequencies_per_level = getSystemDistributionPerGraphLevel(
        GEM, graph_reactions, 'macrosystem')

    macrosystem_frequencies_across_levels = getSystemDistributionAcrossGraphLevels(
        GEM, graph_reactions, 'macrosystem')

    def plot_func(macrosystem, macro_freq_type, level_number, save_fig):
        if macro_freq_type > 0:
            macro_frequencies = macrosystem_frequencies_per_level
            ylabel = 'frequency in level'
        else:
            macro_frequencies = macrosystem_frequencies_across_levels
            ylabel = 'frequency'
        total_subsystems = {}
        level_reactions = graph_reactions[level_number]
        fig, axs = plt.subplots(nrows=2, ncols=1, sharex=False, sharey=False, figsize=(14, 12))
        plt.subplots_adjust(wspace=None, hspace=0.3)

        # Plot macrosystems
        df1 = pd.DataFrame.from_dict(macro_frequencies[macrosystem])
        ax1 = df1.plot.bar(ax=axs[0], rot=0, fontsize=12)
        for p in ax1.patches:
            height = round(p.get_height(), 2)
            if height > 0:
                ax1.annotate(format(height, '.2f'),
                             (p.get_x() * 1.005, p.get_height() * 1.008))
        axs[0].set_title(macrosystem, fontsize=16)
        axs[0].set_ylabel(ylabel)
        axs[0].set_xlabel('graph level')

        # Plot subsystems
        for source in level_reactions.keys():
            subsystems = extractMetabolicSystems(GEM, level_reactions[source],
                                                 'subsystem', macrosystem)
            total_subsystems[source] = getListFrequencies(subsystems)

        df2 = pd.DataFrame(dict([(k, pd.Series(v)) for k, v in total_subsystems.items()]))
        try:
            df2.plot(ax=axs[1], kind='bar', rot=75, fontsize=12)
            axs[1].set_title('Subsystems in level ' + str(level_number), fontsize=16)
            axs[1].set_ylabel('frequency')
        except Exception:
            axs[1].set_title('No data available', fontsize=16)

        if save_fig:
            plt.savefig('figure.png', dpi=fig_dpi, bbox_inches="tight")

    interact(plot_func,
             macrosystem=macrosystems,
             macro_freq_type=widgets.Dropdown(
                 options=[('across graph levels', 0), ('per graph level', 1)],
                 value=0,
                 description='frequency'),
             level_number=widgets.IntSlider(
                 value=0, min=0,  max=20, step=1, description='graph level', readout=True),
             save_fig=widgets.ToggleButton(
                 value=False, description='Save figure', disabled=False,
                 layout=widgets.Layout(margin='3% 0 3% 8%'))
             )

def plotLevelSubsystemsStatic(Model, graph_reactions):
    """
    Computes the frequencies of the subsystems of the reactions appearing
    in the specified graph level which have the specified macrosystem
    """
    GEM = Model.GEM
    macrosystems = ['Amino acid metabolism', 'Carbohydrate metabolism', 'Cell wall biosynthesis',
                    'Cofactor and vitamin metabolism', 'Energy and maintenance', 'Lipid metabolism',
                    'Nucleotide metabolism', 'Transport']
    macrosystem_frequencies_per_level = getSystemDistributionPerGraphLevel(
        GEM, graph_reactions, 'macrosystem')
    plt.rcParams.update({'figure.max_open_warning': 0})

    def plot_func(a_macrosystem, b_level_number):

        total_subsystems = {}
        macro_frequencies = macrosystem_frequencies_per_level
        ylabel = 'frequency in level'
        level_reactions = graph_reactions[b_level_number]
        fig, axs = plt.subplots(nrows=2, ncols=1, sharex=False, sharey=False, figsize=(14, 12))
        plt.subplots_adjust(wspace=None, hspace=0.3)

        # Plot macrosystems
        df1 = pd.DataFrame.from_dict(macro_frequencies[a_macrosystem])
        ax1 = df1.plot.bar(ax=axs[0], rot=0, fontsize=12)
        for p in ax1.patches:
            height = round(p.get_height(), 2)
            if height > 0:
                ax1.annotate(format(height, '.2f'),
                             (p.get_x() * 1.005, p.get_height() * 1.008))
        axs[0].set_title(a_macrosystem, fontsize=16)
        axs[0].set_ylabel(ylabel)
        axs[0].set_xlabel('graph level')

        # Plot subsystems
        for source in level_reactions.keys():
            subsystems = extractMetabolicSystems(GEM, level_reactions[source],
                                                 'subsystem', a_macrosystem)
            total_subsystems[source] = getListFrequencies(subsystems)

        df2 = pd.DataFrame(dict([(k, pd.Series(v)) for k, v in total_subsystems.items()]))
        try:
            df2.plot(ax=axs[1], kind='bar', rot=75, fontsize=12)
            axs[1].set_title('Subsystems in level ' + str(b_level_number), fontsize=16)
            axs[1].set_ylabel('frequency in level')
        except Exception:
            axs[1].set_title('No data available', fontsize=16)

        return fig

    static_fig = StaticInteract(plot_func,
                                a_macrosystem=DropDownWidget(macrosystems,
                                                             description='macrosystem'),
                                b_level_number=RangeWidget(0, 20, step=1, default=0,
                                                           width=200, description='graph level'),
                                interact_name='DAG_systems'
                                )
    return static_fig

def plotTotalSystemsDistribution(Model, shared_reactions, exclusive_reactions,
                                 systemType='macrosystem',
                                 carbonSources=['glucose', 'acetate', 'glycerate']):
    """
    Generates a stacked barplot depicting the distribution of macro or subsystems found in
    the totality of the reactions present in the flux-order graph per carbon source.
    Returns the pandas dataframe employed to generate de plot.
    """
    # Compute frequencies
    GEM = Model.GEM
    total_shared_system_frequencies = getSystemFrequencies(
        GEM, shared_reactions, systemType)

    total_exclusive_system_frequencies = {}
    for carbonSource in carbonSources:
        total_exclusive_system_frequencies[carbonSource] = getSystemFrequencies(
            GEM, exclusive_reactions, systemType, carbonSource)

    # Plot stacked bar plots of total macrosystem frequencies
    df = pd.DataFrame.from_dict(total_exclusive_system_frequencies,
                                orient='index')
    df2 = pd.Series(total_shared_system_frequencies, name='Shared')
    df3 = df.append(df2, ignore_index=False)
    # Plot
    df3.plot.barh(stacked=True, rot=0, legend=False, figsize=(12, 8), fontsize=16)
    plt.xlabel('frequency', fontsize=16)
    plt.title('Macrosystems across levels', fontsize=18, pad=20)
    plt.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
    plt.subplots_adjust(right=0.6)
    plt.show()
    return df3

def plotNumberOfReactionsPerGraphLevel(Graph_Levels):
    """
    Plot barplot with the number of reactions per graph level for each carbon source
    """
    number_of_reactions_per_level, total_number_reactions = {}, {}
    level_with_P90_reactions = {}
    carbonSources = Graph_Levels.keys()
    max_number_of_levels = max([len(Graph_Levels[source]) for source in carbonSources])

    for carbonSource in carbonSources:
        number_of_reactions_per_level[carbonSource] = []
        total_number_reactions[carbonSource] = sum([len(Graph_Levels[carbonSource][L])
                                                    for L in Graph_Levels[carbonSource].keys()])
        for level_number in range(max_number_of_levels):
            try:
                n_reactions = len(
                    Graph_Levels[carbonSource][level_number])
            except Exception:
                n_reactions = 0
            number_of_reactions_per_level[carbonSource].append(n_reactions)
            cum_sum = sum(number_of_reactions_per_level[carbonSource])
            is_P90_level = ((cum_sum >= 0.9 * total_number_reactions[carbonSource])
                            & (carbonSource not in level_with_P90_reactions.keys()))
            if is_P90_level:
                level_with_P90_reactions[carbonSource] = level_number

    df = pd.DataFrame.from_dict(number_of_reactions_per_level)
    ax = df.plot.bar(rot=0, figsize=(12, 8))
    ax.set_xlabel('level', fontsize=16)
    ax.set_ylabel('counts', fontsize=16)
    ax.set_title('Number of reactions across DAG levels', fontsize=17, pad=20)
    # Add vertical line in level with >= 90% of reactions
    n_P90 = max(level_with_P90_reactions.values())
    x_P90_position = 1.15 * ax.patches[n_P90].get_x()
    ax.axvline(x_P90_position, color='k', linestyle='--')
    ax.text(1.1 * x_P90_position, 500, '$P_{{90}}$ (level: {})'.format(n_P90), fontsize=18)
    return df

def plotMacrosystemFrequencyPerGraphLevel(macrosystem_frequencies):
    """
    Plot barplot with the frequency of macrosystems per
    graph level for each carbon source
    """
    if isinstance(next(iter(macrosystem_frequencies.values())), dict):
        data_type = 'exclusive'
    else:
        data_type = 'shared'

    fig, axs = plt.subplots(nrows=4, ncols=2, sharex=False,
                            sharey=False, figsize=(14, 18))
    fig.text(0.5, 0.04, 'level', ha='center', fontsize=14)
    fig.text(0.04, 0.5, 'frequency per level', va='center', rotation='vertical', fontsize=14)
    fig.suptitle(('Metabolic macrosystems of '
                  + data_type
                  + ' ordered reaction pairs per hierarchy level'),
                 fontsize=14)
    plt.subplots_adjust(wspace=None, hspace=0.5)

    if data_type is 'shared':
        df = pd.DataFrame.from_dict(macrosystem_frequencies)
        if 'Unassigned' in df.columns:
            df.drop(['Unassigned'], axis=1)
        df.plot.bar(ax=axs, subplots=True, legend=False, rot=0, fontsize=8)
    else:
        i = 0
        for system in macrosystem_frequencies.keys():
            df = pd.DataFrame.from_dict(macrosystem_frequencies[system])
            if 'Unassigned' in df.columns:
                df.drop(['Unassigned'], axis=1)
            df.plot.bar(ax=axs[i // 2][i % 2], rot=0, title=system)
            i += 1

def plotOrderedCoupledPieCharts(data, ax, pie_name=None, colors=None):
    """
    Plots pie charts with the fraction of coupled reaction pairs
    """
    size = data['values']
    names = data['labels']
    _, texts = ax.pie(size, labels=names, colors=colors,
                      shadow=True, textprops={'fontsize': 14},
                      labeldistance=0.7, startangle=0,
                      explode=(0.05, 0.3, 0.1, 0.05))
    innerCircle = plt.Circle((0, 0), 0.5, color='white')
    ax.text(-0.27, 0, pie_name, {'fontsize': 16})
    for t in texts:
        t.set_fontsize(16)
    ax.add_artist(innerCircle)

def plotFractionOfCoupledOrdered(data):
    df = pd.DataFrame.from_dict(data)
    df.plot(kind='bar', rot=40, stacked=True, fontsize=18, figsize=(10, 9),
            color=['#dddddd', '#1d93f4'], legend=True)

    plt.ylabel('%', fontsize=16)
    plt.title('Ordered pairs per coupling type', fontsize=14)
    plt.xticks(np.arange(3), ('partial', 'directional'), fontsize=14)
    plt.show()

def plotOrderExample(min_flux=0, max_flux=1000, range_step=100,
                     default_value=600, description='$v_j^{ub}$'):

    def plot_example(fluxUpperBound):

        fig = plt.figure(figsize=(10, 8))
        plt.xlabel('$v_j$', fontsize=18)
        plt.ylabel('$v_i$', fontsize=18)
        x = np.linspace(0, 1000, 10)
        plt.vlines(fluxUpperBound, 1000, 0, colors='orange', linewidth=3)

        line1, = plt.plot(x, x, color='black', linestyle='--', label='$v_i=v_j$')
        line2, = plt.plot(x, 0.5*x + 200, color='blue', linewidth=3, label='$v_i=0.5v_j+200$')
        line3, = plt.plot(x, 2*x + 250, color='red', linewidth=3, label='$v_i=2v_j+250$')
        plt.ylim(0, 1000)
        plt.xlim(0, 1000)
        xfill = np.linspace(0, fluxUpperBound, 10)
        plt.fill_between(xfill, 0.5*xfill + 200, 2*xfill + 250, color='lightgrey')

        plt.legend(handles=[line1, line2, line3], loc=4, frameon=False)
        plt.text(300, 600, '$v_i>v_j$', fontsize=18)
        plt.text(600, 300, '$v_i<v_j$', fontsize=18)

        return fig

    interactive_fig = StaticInteract(plot_example,
                   fluxUpperBound=RangeWidget(min_flux, max_flux, step=range_step,
                                              default=default_value, width=200,
                                              description=description,
                                              show_range=False),
                   interact_name='plot_example'
                   )
    return interactive_fig

def plotSolutionSpace(model, v_i, v_j, v_i_max=None, n_points=100,
                      plot_diagonal=True, r_i_nickname=None, r_j_nickname=None,
                      plot_title=None):
    """
    Plots the feasible solution space of reactions v_i and v_j
    (v_i in the x axis) in a 2D plot.

    Arguments:
    ---------
    model: cobra model object, the GEM model of the system
    v_i, v_j: str or int, the reaction ids or indices in the GEM
    """
    warnings.filterwarnings('ignore')
    GEM = copy.deepcopy(model)
    rxn_i, rxn_j = GEM.reactions.get_by_any(v_i)[0], GEM.reactions.get_by_any(v_j)[0]
    if r_i_nickname is None:
        rxn_i_id = rxn_i.id
    else:
        rxn_i_id = r_i_nickname
    if r_j_nickname is None:
        rxn_j_id = rxn_j.id
    else:
        rxn_j_id = r_j_nickname
    if v_i_max is None:
        v_i_max = rxn_i.upper_bound

    GEM.objective = rxn_j
    v_i_seq = np.linspace(0, v_i_max, n_points)
    v_j_min_seq, v_j_max_seq = np.zeros(len(v_i_seq)), np.zeros(len(v_i_seq))
    for i, v in enumerate(v_i_seq):
        rxn_i.upper_bound = v
        GEM.objective.direction = 'min'
        sol_min = GEM.optimize()
        # print('min status ' + str(i) + ' ' + sol_min.status)
        if sol_min.status != 'infeasible':
            v_j_min_seq[i] = sol_min.objective_value

        GEM.objective.direction = 'max'
        sol_max = GEM.optimize()
        # print('max_status ' + str(i) + ' ' + sol_max.status)
        if sol_max.status != 'infeasible':
            v_j_max_seq[i] = sol_max.objective_value

    plt.figure(figsize=(12, 8))
    plt.xlabel('$v_{' + rxn_i_id + '}$')
    plt.ylabel('$v_{' + rxn_j_id + '}$')
    if plot_title is not None:
        plt.title(plot_title, fontsize=18)
    if plot_diagonal:
        plt.plot(v_i_seq, v_i_seq, color='green', linestyle='dashed',
                 label='$v_{' + rxn_i_id + '}= v_{' + rxn_j_id + '}$')
        plt.text(v_i_max / 3, 1.1 * v_i_max / 4,
                 '$v_{' + rxn_i_id + '}>' + 'v_{' + rxn_j_id + '}$', fontsize=18)
        plt.text(1.3 * v_i_max / 3, 3 * v_i_max / 4,
                 '$v_{' + rxn_i_id + '}<' + 'v_{' + rxn_j_id + '}$', fontsize=18)
    plt.plot(v_i_seq, v_j_min_seq, color='blue',
             label='$v^{min}_{' + rxn_j_id + '}$')
    plt.plot(v_i_seq, v_j_max_seq, color='orange',
             label='$v^{max}_{' + rxn_j_id + '}$')
    plt.fill_between(v_i_seq, v_j_min_seq, v_j_max_seq, color='lightgrey')
    plt.legend()
    plt.show()
    warnings.resetwarnings()

def plotFluxSampleOfOrderedChain(Model, ordered_chain):
    """
    Obtains a random flux sample and plots the flux values for the reactions
    in the supplied flux-ordered reaction chain. Model is a Model object from
    Model_module, ordered_chain a list with reaction ids of the cobra model.
    """
    def getLowerBound(GEM, rxn_id):
        return GEM.reactions.get_by_id(rxn_id).lower_bound

    fraction_of_optimum = 0.95
    target_rxn_id = ordered_chain[-1]
    old_rxn_lb = getLowerBound(Model.GEM, target_rxn_id)
    new_rxn_lb = fraction_of_optimum * cobra.flux_analysis.flux_variability_analysis(
        Model.GEM, reaction_list=target_rxn_id, fraction_of_optimum=0).maximum.item()
    try:
        Model.GEM.reactions.get_by_id(target_rxn_id).lower_bound = new_rxn_lb
        # fluxSample = iJO1366[example_source].getFluxSampleInRprogram(nsamples=500)
        # use this function (in R) if cobra sampling doesn't work
        fluxSample = Model.getFluxSample(nsamples=100)
    except Exception:
        print('Failed, try different lower bound')

    Model.GEM.reactions.get_by_id(target_rxn_id).lower_bound = old_rxn_lb

    # Plot flux sample for the reactions in the chain
    fluxSample[ordered_chain].transpose().plot(legend=False, figsize=(12, 8))
    plt.title('Flux sample of ordered reaction chain', fontsize=18, pad=20)
    plt.ylabel('v', fontsize=18)
    plt.xticks(range(len(ordered_chain)), ordered_chain, fontsize=16)
    plt.show()

def plotNumberFluxOrderedAtP0WithData(plot_data):
    """
    Bar plot with number of flux-ordered reaction pairs per carbon source
    across minimum biomass production cutoff values (alpha)
    """
    cost_P = 0
    alphas = list(plot_data['glucose'].keys())
    sources = plot_data.keys()
    data = dict.fromkeys(alphas, [])
    for alpha in alphas:
        counts = []
        for source in sources:
            try:
                counts.append(plot_data[source][alpha][cost_P]['distribution_size'])
            except Exception:
                counts.append(0)
        data[alpha] = counts
    df = pd.DataFrame.from_dict(data, orient='index', columns=sources)
    df.plot.bar(rot=0, figsize=(12,8))
    plt.xlabel('alpha')
    plt.ylabel('counts')
    plt.show()
    return df

def plotNumberOfFluxOrderedPairs(path_to_data=None, alphas=[0, 0.5, 0.75, 0.90, 0.925, 0.95, 1],
                                 sources = ['glucose', 'glycerate', 'acetate']):
    """
    Bar plot with the number of flux ordered pairs for the three carbon sources
    across parameter alpha (biomass lower bound) values
    """
    data = {}
    sources = ['glucose', 'glycerate', 'acetate']
    for alpha in alphas:
        counts = []
        for source in sources:
            directory = path_to_data + f'/{source}/iJO1366_{source[:3]}_A_alpha={alpha}.csv'
            A = np.genfromtxt(directory, delimiter=',')
            counts.append(sum(sum(A)))
        data[alpha] = counts

    df = pd.DataFrame.from_dict(data, orient='index', columns=sources)
    df.plot.bar(rot=0, figsize=(12,8))
    plt.xlabel('alpha')
    plt.ylabel('counts')
    plt.show()
    return df

max_stars_repo_name: Robaina/fluxOrders

max_issues_repo_name: Robaina/fluxOrders

max_forks_repo_name: Robaina/fluxOrders


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

from erebos import utils

GOES17_POS = (-28648.05623634, 0.0, -30937.08652789)
GOES17_LATLONH = (0.0, -137.2, 35786.0)
APPARENT_POS = (-5157.372532, 3567.386338, -1126.998925)
APPARENT_LATLONH = (34.04925, -102.326614, 0.0)
CORRECTED_POS = (-5167.923661, 3565.784003, -1140.388491)
CORRECTED_LATLONH = (33.971087, -102.443847, 10.0)

@pytest.fixture(
    params=(
        (GOES17_POS, GOES17_LATLONH),
        (APPARENT_POS, APPARENT_LATLONH),
        (CORRECTED_POS, CORRECTED_LATLONH),
    )
)
def known_positions(request):
    return request.param

def test_rotatedecrposition_from_geodetic(known_positions):
    pos = utils.RotatedECRPosition.from_geodetic(*known_positions[1])
    expected = utils.RotatedECRPosition(*known_positions[0])
    np.testing.assert_allclose(pos, expected, rtol=1e-5)

def test_rotatedecrposition_roundtrip_latlon(known_positions):
    exp = known_positions[1][:2]
    pos = utils.RotatedECRPosition.from_geodetic(*known_positions[1])
    back = pos.to_geodetic()
    np.testing.assert_allclose(exp, back)

def test_find_actual_cloud_position_no_height():
    app_cloud = utils.RotatedECRPosition(*APPARENT_POS)
    sat_pos = utils.RotatedECRPosition(*GOES17_POS)
    actual = utils.find_actual_cloud_position(sat_pos, app_cloud, 0)
    np.testing.assert_allclose(actual, app_cloud, rtol=1e-6)

def test_find_actual_cloud_position():
    app_cloud = utils.RotatedECRPosition(*APPARENT_POS)
    sat_pos = utils.RotatedECRPosition(*GOES17_POS)
    actual = utils.find_actual_cloud_position(sat_pos, app_cloud, CORRECTED_LATLONH[-1])
    exp = utils.RotatedECRPosition(*CORRECTED_POS)
    np.testing.assert_allclose(actual, exp)

def test_find_apparent_cloud_position():
    cloud = utils.RotatedECRPosition(*CORRECTED_POS)
    sat_pos = utils.RotatedECRPosition(*GOES17_POS)
    app = utils.find_apparent_cloud_position(sat_pos, cloud)
    exp = utils.RotatedECRPosition(*APPARENT_POS)
    np.testing.assert_allclose(app, exp, rtol=1e-6)

def test_find_apparent_cloud_position_terrain():
    cloud = utils.RotatedECRPosition(*CORRECTED_POS)
    sat_pos = utils.RotatedECRPosition(*GOES17_POS)
    app = utils.find_apparent_cloud_position(sat_pos, cloud, CORRECTED_LATLONH[-1])
    # no shift when terrain is as high as cloud
    np.testing.assert_allclose(app, cloud, rtol=1e-6)

def test_find_apparent_cloud_position_on_surface():
    cloud = utils.RotatedECRPosition(*APPARENT_POS)
    sat_pos = utils.RotatedECRPosition(*GOES17_POS)
    app = utils.find_apparent_cloud_position(sat_pos, cloud)
    np.testing.assert_allclose(app, cloud, rtol=1e-6)

def test_find_position_roundtrip():
    app_cloud = utils.RotatedECRPosition(*APPARENT_POS)
    sat_pos = utils.RotatedECRPosition(*GOES17_POS)
    actual = utils.find_actual_cloud_position(sat_pos, app_cloud, CORRECTED_LATLONH[-1])
    app = utils.find_apparent_cloud_position(sat_pos, actual)
    np.testing.assert_allclose(app, app_cloud, rtol=1e-6)

def test_get_solar_ecr_position():
    time = pd.Timestamp("20190726T125619-06:00")
    cloud = utils.RotatedECRPosition(*CORRECTED_POS)
    sol = utils.get_solar_ecr_position(
        time.value / 1e9, CORRECTED_LATLONH[0], CORRECTED_LATLONH[1]
    )
    pos = utils.find_apparent_cloud_position(sol, cloud)
    # sun should be at same latitude at this time
    np.testing.assert_allclose(pos.to_geodetic()[1], cloud.to_geodetic()[1])
    assert pos.to_geodetic()[0] > cloud.to_geodetic()[0]

ext: py

reason: import numpy


text:
import os
import argparse
import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torchvision
from torch.utils.data import DataLoader
from torchvision import models, datasets, transforms
import torch.optim as optim
import torch.functional as F
from torch.utils.tensorboard import SummaryWriter

parser = argparse.ArgumentParser()

parser.add_argument("-d", "--dataset", required=True, help="Location of the dataset folder")
parser.add_argument("-b", "--batch_size", required = True, help = "Batch Size of the dataset for training")
parser.add_argument("-ic","--image_channels", required=True, help="Channels in the input image")
parser.add_argument("-e", "--epochs", required = True, help = "Number of epochs to run the training loop")

args = vars(parser.parse_args())

class Discriminator(nn.Sequential):

    def __init__(self, channels_img, features_d):
        # Did not use BatchNorm in the last layer of the generator and the first layer of the
        # discriminator
        # Input: N x channels_img x 64 x 64
        modules = [nn.Conv2d(channels_img, features_d, kernel_size = 4, stride = 2, padding = 1), #32x32
                   nn.LeakyReLU(0.2, inplace=True),
                   self._block(features_d, features_d*2, 4, 2, 1),# 16x16
                   self._block(features_d*2, features_d*4, 4, 2, 1), #8x8
                   self._block(features_d*4, features_d*8, 4, 2, 1), #4x4
                   nn.Conv2d(features_d*8, 1, kernel_size = 4, stride = 2, padding = 0), #1x1(Prediction)
                   nn.Sigmoid()]

        super(Discriminator, self).__init__(*modules)

    def _block(self, in_channels, out_channels, kernel_size, stride, padding):
        return nn.Sequential(
            nn.Conv2d(in_channels, out_channels,kernel_size, stride, padding, bias=False),
            nn.BatchNorm2d(out_channels),
            nn.LeakyReLU(0.2, inplace=True)
        )

class Generator(nn.Sequential):

    # Here channels_img is nothing but the inputs channels
    # and features_g is nothing but the output channels
    def __init__(self, z_dim, channels_img, features_g):

        modules = [self._block(z_dim, features_g*16, 4, 1, 0),
                   self._block(features_g*16, features_g*8, 4, 2, 1),
                   self._block(features_g*8, features_g*4, 4, 2, 1),
                   self._block(features_g*4, features_g*2, 4, 2, 1),
                   nn.ConvTranspose2d(features_g*2, channels_img, 4, 2, 1),
                   nn.Tanh()]

        super(Generator, self).__init__(*modules)

    def _block(self, in_channels, out_channels, kernel_size, stride, padding):
        return nn.Sequential(
            nn.ConvTranspose2d(in_channels, out_channels, kernel_size, stride, padding, bias=False),
            nn.BatchNorm2d(out_channels),
            nn.ReLU(inplace=True)
        )

def weights_init(m):
    classname = m.__class__.__name__
    if classname.find('Conv') != -1:
        nn.init.normal_(m.weight.data, 0.0, 0.02)
    elif classname.find('BatchNorm') != -1:
        nn.init.normal_(m.weight.data, 1.0, 0.02)
        nn.init.constant_(m.bias.data, 0)

def weights_init(m):
    classname = m.__class__.__name__
    if classname.find('Conv') != -1:
        nn.init.normal_(m.weight.data, 0.0, 0.02)
    elif classname.find('BatchNorm') != -1:
        nn.init.normal_(m.weight.data, 1.0, 0.02)
        nn.init.constant_(m.bias.data, 0)

# def test():
#     N, in_channels, H, W = 8, 3, 64, 64
#     z_dim = 100
#     X = torch.randn((N, in_channels, H, W))
#     disc = Discriminator(in_channels,8)
#     assert disc(X).shape == (N, 1, 1, 1) # One Value per example
#     gen = Generator(z_dim, in_channels, 64)
#     z = torch.randn((N, z_dim, 1, 1))
#     assert gen(z).shape == (N, in_channels, H, W) # Ouput Generated image
#     print("Success")
# test()

## The training Setup
root = "C:\\Users\\shant\\img_align_celeba\\img_align_celeba\\"
LEARNING_RATE = 2e-4
BATCH_SIZE = args["batch_size"]
IMAGE_SIZE = 64
# Image Channels in the generator output and input to the discriminator
IMAGE_CHANNELS = args["image_channels"]

# Latent Space Dimensions
Z_DIM = 100
NUM_EPOCHS = args["epochs"]

# feature size in the discriminator
FEATURES_DISC = 64

# feature size in the generator
FEATURES_GEN = 64

# Setup Transforms
dataset = datasets.ImageFolder(os.path.join(args["dataset"]),transform = transforms.Compose([
    transforms.Resize(IMAGE_SIZE),
    transforms.CenterCrop(IMAGE_SIZE),
    transforms.ToTensor(),
    transforms.Normalize([0.5 for _ in range(IMAGE_CHANNELS)], [0.5 for _ in range(IMAGE_CHANNELS)])
]))

# DataLoader
dataloader = DataLoader(dataset, batch_size = BATCH_SIZE, shuffle=True, num_workers = 2)

device  = torch.device("cuda" if torch.cuda.is_available() else "cpu")
generator = Generator(Z_DIM, IMAGE_CHANNELS, FEATURES_GEN).to(device)
discriminator = Discriminator(IMAGE_CHANNELS,FEATURES_DISC).to(device)

generator.apply(weights_init)
discriminator.apply(weights_init)

generator_optimizer = optim.Adam(generator.parameters(), lr = LEARNING_RATE, betas = (0.5,0.999))
discriminator_optimizer = optim.Adam(discriminator.parameters(), lr = LEARNING_RATE, betas = (0.5,0.999))

# Loss
criterion = nn.BCELoss()

# Batch of Latent Vectors
fixed_noise = torch.randn((64, Z_DIM, 1, 1)).to(device)

writer_real = SummaryWriter(f"logs/real")
writer_fake = SummaryWriter(f"logs/fake")
loss_curves = SummaryWriter(f"logs/loss_curves")

step = 0 # Printing to tensorboard
# generator.train()
# discriminator.train()

for epoch in range(NUM_EPOCHS):

    # Unsupervised
    for batch_idx, (real, _) in enumerate(dataloader):
        discriminator.zero_grad()
        # Latent noise
        noise = torch.randn((BATCH_SIZE, Z_DIM, 1, 1)).to(device)
        # The real world images
        real = real.to(device)
        # Pass the latent vector through the generator
        fake = generator(noise)

        #####################################################
        # Train the discriminator max log(D(x)) + log(1-D(G(z)))
        #####################################################

        disc_real = discriminator(real).view(-1)
        # log(D(x)), y_n = ones hence only logx_n term is left. Refer to Pytorch Documentation
        loss_disc_real = criterion(disc_real, torch.ones_like(disc_real))

        loss_disc_real.backward()

        disc_fake = discriminator(fake.detach()).view(-1)
        # Subtracting
        loss_disc_fake = criterion(disc_fake, torch.zeros_like(disc_fake))

        # Addition of gradients from the all the real and fake samples
        loss_disc = (loss_disc_fake + loss_disc_real)

        loss_disc_fake.backward() # Reutilize fake, but pytorch removes intermeditate results
        discriminator_optimizer.step()

        ##################################
        # Train Generator max log(D(G(z)))
        ##################################
        generator.zero_grad()
        output = discriminator(fake).view(-1)
        loss_gen = criterion(output, torch.ones_like(output))
        loss_gen.backward()
        generator_optimizer.step()

        if batch_idx % 100 == 0:

            print(
            f"Epoch [{epoch}/{NUM_EPOCHS}] Batch {batch_idx}/{len(dataloader)} \
                  Loss D: {loss_disc:.4f}, loss G: {loss_gen:.4f}"
            )
            with torch.no_grad():
                fake = generator(fixed_noise)

                # The [:64] prints out the 4-D tensor BxCxHxW
                img_grid_real = torchvision.utils.make_grid(
                    real[:64], normalize = True)
                img_grid_fake = torchvision.utils.make_grid(
                    fake[:64], normalize = True)

                ##########################
                # TensorBoard Visualizations
                ##########################
                writer_real.add_image("Real", img_grid_real, global_step=step)
                writer_fake.add_image("Fake", img_grid_fake, global_step=step)
                loss_curves.add_scalar("generator", loss_gen, global_step=step)
                loss_curves.add_scalar("discriminator", loss_disc, global_step = step)
            step += 1 # See progression of images


text:
# randomly split a transactions file into test and train portions
# the transactions file is pickled and contains a Pandas DataFrame object
# the split files are also pickled and each contains a Pandas DataFrame object
# invocation options:
#  --test=FRACTION       : fraction to test, default is 0.10
#  --in PATHIN           : path to input transaction pickle file
#  --outtest PATHTEST    : path to test file (.pickle)
#  --outtrain PATHTRAIN  : path to training file (.pickle)
#
# The test data consists of a 10% random sample of all the data.
#
# Unlike the R version, the data are not stratified by sale month.

# import built-ins and libraries
import numpy as np
import datetime
import pandas as pd
import pdb
import sys
from sklearn import cross_validation
import warnings

# import my stuff
from Bunch import Bunch
from directory import directory
from Logger import Logger
import parse_command_line

def make_control(argv):
    # return a Bunch
    script_name = argv[0]

    b = Bunch(debugging=False,
              testing=False,
              now=datetime.datetime.now(),
              base_name=script_name.split('.')[0],
              me=script_name,
              random_seed=123,
              arg_fraction=parse_command_line.default(argv, '--fraction', 0.10),
              arg_in=parse_command_line.get_arg(argv, '--in'),
              arg_outtest=parse_command_line.get_arg(argv, '--outtest'),
              arg_outtrain=parse_command_line.get_arg(argv, '--outtrain'))
    return b

def randomly_split(df, fraction_to_testing, random_state):
    '''Randomly shuffly observations and split into train and test.'''

    ss = cross_validation.ShuffleSplit(n=df.shape[0],
                                       n_iter=1,
                                       test_size=fraction_to_testing,
                                       random_state=random_state)

    # extract the train and test indices
    # there should be exactly one set of such
    num_iterations = 0
    for train_index, test_index in ss:
        num_iterations = num_iterations + 1
        train_indices = train_index
        test_indices = test_index
    assert num_iterations == 1

    test = df.iloc[test_indices]
    train = df.iloc[train_indices]

    return test, train

def analyze(test, train):
    print 'test.shape', test.shape
    print 'train.shape', train.shape

    # check range of sale.year
    min_year = min(np.amin(test['sale.year']),
                   np.amin(train['sale.year']))
    max_year = max(np.amax(test['sale.year']),
                   np.amax(train['sale.year']))
    assert min_year == 2003
    print max_year
    assert max_year == 2009

    print 'year, month, #test, #train'
    for year in (2003, 2004, 2005, 2006, 2007, 2008, 2009):
        last_month_index = 12 if year != 2009 else 3
        for month_index in range(last_month_index):
            month = month_index + 1
            is_month_test = np.logical_and(test['sale.year'] == year,
                                           test['sale.month'] == month)
            in_month_test = test[is_month_test]
            is_month_train = np.logical_and(train['sale.year'] == year,
                                            train['sale.month'] == month)
            in_month_train = train[is_month_train]
            print year, month, in_month_test.shape[0], in_month_train.shape[0]

def main():
    warnings.filterwarnings('error')
    control = make_control(sys.argv)
    path = \
        directory('log') + \
        control.base_name + '.' + control.now.isoformat('T') + '.log'
    sys.stdout = Logger(logfile_path=path)  # print x now logs and prints x
    print control

    df = pd.read_pickle(control.arg_in)

    # make sure that sale.python_date is never NaN (null)
    if False:
        dates = df['sale.python_date']
        if dates.isnull().any():
            raise ValueError('at least one sale.python_date is null')

    if control.testing and False:
        pdb.set_trace()
        df = df[0:30]

    # print columns in df
    print 'df.shape', df.shape
    for column_name in df.columns:
        print 'df column name', column_name

    test, train = randomly_split(df=df,
                                 fraction_to_testing=control.arg_fraction,
                                 random_state=control.random_seed)
    analyze(test, train)

    test.to_pickle(control.arg_outtest)
    train.to_pickle(control.arg_outtrain)

    # log the control variables
    print control
    print 'done'

if __name__ == '__main__':
    main()

max_stars_repo_name: rlowrance/re-local-linear

max_issues_repo_name: rlowrance/re-local-linear

max_forks_repo_name: rlowrance/re-local-linear


text:
import random
import os

import numpy as np

from transformers import AutoTokenizer
import torch

def set_seed(seed=42):
    """
    Sets the seed of the entire notebook so results are the same every time we run.
    This is for REPRODUCIBILITY.
    """
    np.random.seed(seed)
    random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    # Torch backend controls the behavior of various backends that PyTorch supports
    # When running on the CuDNN backend, two further options must be set

    # only use deterministic convolution algorithms
    torch.backends.cudnn.deterministic = True
    # does not try different algorithms to choose the fastest
    torch.backends.cudnn.benchmark = False
    # Set a fixed value for the hash seed
    os.environ["PYTHONHASHSEED"] = str(seed)

CONFIG = dict(
    seed=42,
    model_name="../input/roberta-base",
    test_batch_size=64,
    max_length=128,
    num_classes=1,
    device=torch.device("cuda:0" if torch.cuda.is_available() else "cpu"),
)

set_seed(CONFIG["seed"])

CONFIG["tokenizer"] = AutoTokenizer.from_pretrained(CONFIG["model_name"])

max_stars_repo_name: TeamEpochGithub/Jigsaw

max_issues_repo_name: TeamEpochGithub/Jigsaw

max_forks_repo_name: TeamEpochGithub/Jigsaw


text:
# Copyright (c) 2016 The UUV Simulator Authors.
# All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import numpy as np
from scipy.misc import factorial

class BezierCurve(object):
    """
    Implementation of Bezier curves of orders 3, 4 and 5 based on [1].

    [1] Biagiotti, Luigi, and Claudio Melchiorri. Trajectory planning for
        automatic machines and robots. Springer Science & Business Media, 2008.
    """

    def __init__(self, pnts, order, tangents=None):
        assert order in [3, 4, 5], 'Invalid Bezier curve order'
        assert type(pnts) == list and len(pnts) >= 2, 'At least two points are needed to calculate the curve'

        self._pnts = list()
        for pnt in pnts:
            if type(pnt) == list:
                assert len(pnt) == 3, 'Point must have three elements'
                self._pnts.append(np.array(pnt))
            elif type(pnt) == np.ndarray:
                assert pnt.size == 3, 'Point must have three elements'
                self._pnts.append(pnt)
            else:
                raise TypeError('Point in list is neither a list or an array')

        if tangents is not None:
            assert type(tangents) == list and len(tangents) == 2, 'Tangent vectors must be provided'
            for t in tangents:
                if type(t) == list:
                    assert len(t) == 3, 'Tangent vector must have three elements'
                elif type(t) == np.ndarray:
                    assert t.size == 3, 'Tangent vector must have three elements'
                else:
                    raise TypeError('Tangent vector is neither a list or an array')

        self._control_pnts = [np.zeros(3) for _ in range(order + 1)]

        self._order = order
        if self._order == 3:
            assert len(self._pnts) == 2, 'Two points are needed for the curve to be computed'
            # Setting initial and last control points
            self._control_pnts[0] = self._pnts[0]
            self._control_pnts[3] = self._pnts[1]
            # Compute alpha
            a = 16 - np.linalg.norm(tangents[0] + tangents[1])**2
            b = 12 * np.dot(self._control_pnts[3] - self._control_pnts[0], tangents[0] + tangents[1])
            c = -36 * np.linalg.norm(self._control_pnts[3] - self._control_pnts[0])**2
            alpha = np.roots([a, b, c]).max()

            # Compute the rest of the control points
            self._control_pnts[1] = self._control_pnts[0] + (1.0 / 3) * alpha * tangents[0]
            self._control_pnts[2] = self._control_pnts[3] - (1.0 / 3) * alpha * tangents[1]
        elif self._order == 4:
            assert len(self._pnts) == 3, 'Three points are needed for the curve to be computed'
            # Setting initial control points
            self._control_pnts[0] = self._pnts[0]
            self._control_pnts[2] = self._pnts[1]
            self._control_pnts[4] = self._pnts[2]

            radius = np.linalg.norm(self._pnts[0] - self._pnts[1])
            tangents = list()
            tangents.append((self._pnts[1] - self._pnts[0]) / radius)
            tangents.append((self._pnts[2] - self._pnts[1]) / radius)

            # Compute alpha
            a = 4 - (1.0 / 4) * np.linalg.norm(tangents[0] + tangents[1])**2
            b = 3 * np.dot(self._control_pnts[4] - self._control_pnts[0], tangents[0] + tangents[1])
            c = -9 * np.linalg.norm(self._control_pnts[4] - self._control_pnts[0])**2
            alpha = np.roots([a, b, c]).max()

            # Compute the rest of the control points
            self._control_pnts[1] = self._control_pnts[0] + 0.25 * alpha * tangents[0]
            self._control_pnts[3] = self._control_pnts[4] - 0.25 * alpha * tangents[1]
        elif self._order == 5:
            assert len(self._pnts) == 3, 'Three points are needed for the curve to be computed'
            # Setting initial control points
            self._control_pnts[0] = self._pnts[0]
            self._control_pnts[5] = self._pnts[2]

            radius = np.linalg.norm(self._pnts[0] - self._pnts[1])
            tangents = list()
            tangents.append((self._pnts[1] - self._pnts[0]) / radius)
            tangents.append((self._pnts[2] - self._pnts[1]) / radius)

            # Compute alpha
            a = 256 - 49 * np.linalg.norm(tangents[0] + tangents[1])**2
            b = 420 * np.dot(self._control_pnts[5] - self._control_pnts[0], tangents[0] + tangents[1])
            c = -900 * np.linalg.norm(self._control_pnts[5] - self._control_pnts[0])**2
            alpha = np.roots([a, b, c]).max()

            # Compute the rest of the control points
            self._control_pnts[1] = self._control_pnts[0] + 0.2 * alpha * tangents[0]
            self._control_pnts[2] = 2 * self._control_pnts[1] - self._control_pnts[0]
            self._control_pnts[4] = self._control_pnts[5] - 0.2 * alpha * tangents[1]
            self._control_pnts[3] = 2 * self._control_pnts[4] - self._control_pnts[5]

    def control_pnts(self):
        return self._control_pnts

    def interpolate(self, u):
        u = max(u, 0)
        u = min(u, 1)

        b = np.zeros(3)
        n = len(self._control_pnts) - 1
        for i in range(len(self._control_pnts)):
            b = b + self._get_binomial(n, i) * (1 - u)**(n - i) * u**i * self._control_pnts[i]
        return b

    def get_length(self):
        return self._order * np.linalg.norm(self._control_pnts[1] - self._control_pnts[0])

    @staticmethod
    def _get_binomial(n, i):
        return factorial(n) / (factorial(i) * factorial(n - i))

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    print 'Test - Cubic Bezier curve'
    q_x = [0, 1, 2, 4, 5, 6]
    q_y = [0, 2, 3, 3, 2, 0]
    q_z = [0, 1, 0, 0, 2, 2]

    q = np.vstack((q_x, q_y, q_z)).T

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    ax.plot(q[:, 0], q[:, 1], q[:, 2], 'b--')
    ax.plot(q[:, 0], q[:, 1], q[:, 2], 'ro')

    # Compute the distance between points
    lengths = [np.linalg.norm(q[i + 1, :] - q[i, :]) for i in range(q.shape[0] - 1)]
    lengths = [0] + lengths
    u = np.cumsum(lengths) / np.sum(lengths)
    print q.T
    print 'u=', u

    tangents = [np.zeros(3) for _ in range(q.shape[0])]
    delta_u = lambda k: u[k] - u[k - 1]
    delta_q = lambda k: q[k, :] - q[k - 1, :]
    lamb_k = lambda k: delta_q(k) / delta_u(k)
    alpha_k = lambda k: delta_u(k) / (delta_u(k) + delta_u(k + 1))
    for i in range(1, len(u) - 1):
        tangents[i] = (1 - alpha_k(i)) * lamb_k(i) + alpha_k(i) * lamb_k(i + 1)
        if i == 1:
            tangents[0] = 2 * lamb_k(i) - tangents[1]

    tangents[-1] = 2 * lamb_k(len(u) - 1) - tangents[-2]

    print 'Tangents'
    for i in range(len(tangents)):
        tangents[i] = tangents[i] / np.linalg.norm(tangents[i])
        print '\t#%d - %.2f %.2f %.2f' % (i, tangents[i][0], tangents[i][1], tangents[i][2])
        ax.plot([q[i, 0], tangents[i][0] + q[i, 0]],
                [q[i, 1], tangents[i][1] + q[i, 1]],
                [q[i, 2], tangents[i][2] + q[i, 2]], 'c')

    segments = list()
    print 'Segments'
    for i in range(len(tangents) - 1):
        segments.append(BezierCurve([q[i, :], q[i + 1, :]], 3, tangents[i:i + 2]))
        print '\t', segments[-1]._control_pnts

    lengths = [seg.get_length() for seg in segments]
    lengths = [0] + lengths
    total_length = np.sum(lengths)

    u = np.cumsum(lengths) / total_length

    pnts = None
    for i in np.linspace(0, 1, 100):
        idx = (u - i >= 0).nonzero()[0][0]
        if idx == 0:
            u_k = 0
            pnts = segments[idx].interpolate(u_k)
        else:
            u_k = (i - u[idx - 1]) / (u[idx] - u[idx - 1])
            pnts = np.vstack((pnts, segments[idx - 1].interpolate(u_k)))

    ax.plot(pnts[:, 0], pnts[:, 1], pnts[:, 2], 'g')

    fig = plt.figure()

    u_hat = np.linspace(0, 1, 100)
    d_x = [(pnts[i, 0] - pnts[i - 1, 0]) / (u_hat[i] - u_hat[i - 1]) for i in range(1, pnts.shape[0])]
    d_y = [(pnts[i, 1] - pnts[i - 1, 1]) / (u_hat[i] - u_hat[i - 1]) for i in range(1, pnts.shape[0])]
    d_z = [(pnts[i, 2] - pnts[i - 1, 2]) / (u_hat[i] - u_hat[i - 1]) for i in range(1, pnts.shape[0])]

    ax = fig.add_subplot(311)
    ax.plot(u_hat[1::], d_x)
    ax.set_xlim([u_hat[1::].min(), u_hat[1::].max()])
    ax.grid(True)

    ax = fig.add_subplot(312)
    ax.plot(u_hat[1::], d_y)
    ax.set_xlim([u_hat[1::].min(), u_hat[1::].max()])
    ax.grid(True)

    ax = fig.add_subplot(313)
    ax.plot(u_hat[1::], d_z)
    ax.set_xlim([u_hat[1::].min(), u_hat[1::].max()])
    ax.grid(True)

    plt.show()

reason: import numpy,from scipy


text:
# This shows how to save ImageArray objects, such as those returned from
# volume renderings, to pngs with varying backgrounds.

# First we use the simple_volume_rendering.py recipe from above to generate
# a standard volume rendering.  The only difference is that we use
# grey_opacity=True with our TransferFunction, as the colored background
# functionality requires images with an opacity between 0 and 1.

# We have removed all the comments from the volume rendering recipe for
# brevity here, but consult the recipe for more details.

import yt
import numpy as np

ds = yt.load("Enzo_64/DD0043/data0043")
ad = ds.all_data()
mi, ma = ad.quantities.extrema("density")
tf = yt.ColorTransferFunction((np.log10(mi)+1, np.log10(ma)), grey_opacity=True)
tf.add_layers(5, w=0.02, colormap="spectral")
c = [0.5, 0.5, 0.5]
L = [0.5, 0.2, 0.7]
W = 1.0
Npixels = 512
cam = ds.camera(c, L, W, Npixels, tf)
im = cam.snapshot("original.png" % ds, clip_ratio=8.0)

# Our image array can now be transformed to include different background
# colors.  By default, the background color is black.  The following
# modifications can be used on any image array.

# write_png accepts a background keyword argument that defaults to 'black'.
# Other choices include:
# black (0.,0.,0.,1.)
# white (1.,1.,1.,1.)
# None  (0.,0.,0.,0.) <-- Transparent!
# any rgba list/array: [r,g,b,a], bounded by 0..1

# We include the clip_ratio=8 keyword here to bring out more contrast between
# the background and foreground, but it is entirely optional.

im.write_png('black_bg.png', background='black', clip_ratio=8.0)
im.write_png('white_bg.png', background='white', clip_ratio=8.0)
im.write_png('green_bg.png', background=[0.,1.,0.,1.], clip_ratio=8.0)
im.write_png('transparent_bg.png', background=None, clip_ratio=8.0)


text:
#!/usr/bin/env python3

import numpy as np
import copy as cp
import os
import pickle
import sys
import warnings

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

from lib import regression as reg
from lib import metrics
from lib import bootstrap as bs
from lib import cross_validation as cv
from lib import logistic_regression as logreg

import sklearn.model_selection as sk_modsel
import sklearn.preprocessing as sk_preproc
import sklearn.linear_model as sk_model
import sklearn.metrics as sk_metrics
import sklearn.utils as sk_utils

from tqdm import tqdm

from task_tools import read_t, load_pickle, save_pickle, print_parameters, \
    plot_accuracy_scores, retrieve_2d_ising_data

from matplotlib import rc, rcParams
rc("text", usetex=True)
rc("font", **{"family": "sans-serif", "serif": ["Computer Modern"]})
rcParams["font.family"] += ["serif"]

def task1c(sk=False, figure_folder="../fig"):
    """Task c) of project 2."""
    print("="*80)
    print("Task c: Logistic regression classification")

    training_size = 0.8
    data_size = 1000
    # options: float, inverse
    learning_rate = "inverse"
    max_iter = int(1e3)
    tolerance = 1e-6
    momentum = 0.0
    mini_batch_size = 20
    data_path = "../datafiles/MehtaIsingData"
    verbose=True
    try_get_pickle=True

    # Available solvers:
    # ["lr-gd", "gd", "cg", "sga", "sga-mb", "nr", "newton-cg"]
    solver = "lr-gd"
    # Available activation functions: sigmoid, softmax
    activation = "sigmoid"
    # Available penalites: l1, l2, elastic_net
    penalty = "elastic_net"

    # Define regularisation parameter
    lmbdas = np.logspace(-4, 3, 8)

    X, y = retrieve_2d_ising_data(data_path, data_size)

    # pick random data points from ordered and disordered states
    # to create the training and test sets
    X_train, X_test, y_train, y_test = sk_modsel.train_test_split(
        X, y, test_size=1-training_size)

    # full data set
    # X=np.concatenate((X_critical,X))
    # Y=np.concatenate((Y_critical,Y))

    print('X_train shape(samples, input):', X_train.shape)
    print('y_train shape(samples, input):', y_train.shape)

    run_sk_comparison(X_train, X_test, y_train, y_test,
                      lmbdas=lmbdas,
                      penalty="l2",
                      activation=activation,
                      solver=solver,
                      learning_rate=learning_rate,
                      momentum=momentum,
                      mini_batch_size=mini_batch_size,
                      max_iter=max_iter,
                      tolerance=tolerance,
                      try_get_pickle=try_get_pickle,
                      verbose=verbose)

    run_lambda_penalty(X_train, X_test, y_train, y_test,
                       lmbdas=lmbdas,
                       use_sk_learn=False,
                       penalties=["l1", "l2", "elastic_net"],
                       activation=activation,
                       solver=solver,
                       learning_rate=learning_rate,
                       momentum=momentum,
                       mini_batch_size=mini_batch_size,
                       max_iter=max_iter,
                       tolerance=tolerance,
                       try_get_pickle=try_get_pickle,
                       verbose=verbose)

    run_lambda_solver(X_train, X_test, y_train, y_test,
                      lmbdas=lmbdas,
                      penalty="l2",
                      activation=activation,
                      solvers=["lr-gd", "cg", "newton-cg"],
                      learning_rate="inverse",
                      momentum=momentum,
                      mini_batch_size=mini_batch_size,
                      max_iter=max_iter,
                      tolerance=tolerance,
                      try_get_pickle=try_get_pickle,
                      verbose=verbose)

    run_lambda_learning_rate_comparison(X_train, X_test, y_train, y_test,
                                        lmbdas=lmbdas,
                                        penalty="l2",
                                        activation=activation,
                                        solver="lr-gd",
                                        learning_rates=["inverse", 0.001,
                                                        0.005, 0.01, 0.05,
                                                        0.1, 0.5, 1.0],
                                        momentum=momentum,
                                        mini_batch_size=mini_batch_size,
                                        max_iter=max_iter,
                                        tolerance=tolerance,
                                        try_get_pickle=try_get_pickle,
                                        verbose=verbose)

    run_lambda_momentum(X_train, X_test, y_train, y_test,
                        lmbdas=lmbdas,
                        penalty="l2",
                        activation=activation,
                        solver="lr-gd",
                        learning_rate=learning_rate,
                        momentums=[0.0001, 0.0005, 0.001, 0.005, 0.01],
                        mini_batch_size=mini_batch_size,
                        max_iter=max_iter,
                        tolerance=tolerance,
                        try_get_pickle=try_get_pickle,
                        verbose=verbose)

def run_lambda_penalty(X_train, X_test, y_train, y_test,
                       lmbdas=None,
                       penalties=None,
                       try_get_pickle=True,
                       figure_folder="../fig",
                       **kwargs):

    param_dict = {"lmbdas": lmbdas,
                  "sklearn": True,
                  "figure_folder": figure_folder,
                  "penalties": penalties}
    param_dict.update(kwargs)
    print_parameters(**param_dict)

    test_accuracy_values = []
    train_accuracy_values = []

    for penalty in penalties:
        print("Regularisation:", penalty)

        pickle_fname = ("lambda_penalty_accuracy_penalty{}_actsigmoid"
                        "_solver{}_lrinverse_mom0.0_tol1e-06."
                        "pkl".format(penalty, kwargs["solver"]))

        if os.path.isfile(pickle_fname) and try_get_pickle:
            res_ = load_pickle(pickle_fname)
        else:
            res_ = logreg_core(X_train, X_test, y_train, y_test,
                               lmbdas=lmbdas,
                               penalty=penalty,
                               store_pickle=True,
                               pickle_fname=pickle_fname,
                               **kwargs)

        train_accuracy, test_accuracy, critical_accuracy = res_

        test_accuracy_values.append(test_accuracy)
        train_accuracy_values.append(train_accuracy)

    plot_accuracy_scores(lmbdas, train_accuracy_values, test_accuracy_values,
                         [r"$L^1$", r"$L^2$", r"Elastic net"],
                         "accuracy_regularisation_scores", r"$\lambda$",
                         r"Accuracy")

def run_lambda_solver(X_train, X_test, y_train, y_test,
                      lmbdas=None,
                      solvers=None,
                      try_get_pickle=True,
                      figure_folder="../fig",
                      **kwargs):

    param_dict = {"lmbdas": lmbdas,
                  "sklearn": True,
                  "figure_folder": figure_folder,
                  "solver": solvers}
    param_dict.update(kwargs)
    print_parameters(**param_dict)

    test_accuracy_values = []
    train_accuracy_values = []

    for solver in solvers:
        print("Solver:", solver)

        pickle_fname = ("lambda_solver_accuracy_penalty{}_actsigmoid"
                        "_solver{}_lr{}_mom0.0_tol1e-06."
                        "pkl".format(kwargs["penalty"], solver,
                                     kwargs["learning_rate"]))

        if os.path.isfile(pickle_fname) and try_get_pickle:
            res_ = load_pickle(pickle_fname)
        else:
            res_ = logreg_core(X_train, X_test, y_train, y_test,
                               use_sk_learn=False,
                               lmbdas=lmbdas,
                               solver=solver,
                               store_pickle=True,
                               pickle_fname=pickle_fname,
                               **kwargs)

        train_accuracy, test_accuracy, critical_accuracy = res_

        test_accuracy_values.append(test_accuracy)
        train_accuracy_values.append(train_accuracy)

    plot_accuracy_scores(lmbdas, train_accuracy_values, test_accuracy_values,
                         [r"Optimized Gradient Descent", r"Conjugate Gradient",
                          r"Newtons-CG"],
                         "accuracy_solver_scores", r"$\lambda$",
                         r"Accuracy")

def run_lambda_learning_rate_comparison(X_train, X_test, y_train, y_test,
                                        lmbdas=None,
                                        learning_rates=None,
                                        try_get_pickle=True,
                                        figure_folder="../fig",
                                        **kwargs):

    param_dict = {"lmbdas": lmbdas,
                  "sklearn": True,
                  "figure_folder": figure_folder,
                  "learning_rates": learning_rates}
    param_dict.update(kwargs)
    print_parameters(**param_dict)

    test_accuracy_values = []
    train_accuracy_values = []

    for lr in learning_rates:
        print ("Learning rate: ", lr)

        pickle_fname = ("lambda_lr_accuracy_penalty{}_actsigmoid"
                        "_solver{}_lr{}_mom0.0_tol1e-06."
                        "pkl".format(kwargs["solver"],
                                     kwargs["penalty"], str(lr)))

        if os.path.isfile(pickle_fname) and try_get_pickle:
            res_ = load_pickle(pickle_fname)
        else:
            res_ = logreg_core(X_train, X_test, y_train, y_test,
                               use_sk_learn=False,
                               lmbdas=lmbdas,
                               learning_rate=lr,
                               store_pickle=True,
                               pickle_fname=pickle_fname,
                               **kwargs)

        train_accuracy, test_accuracy, critical_accuracy = res_

        test_accuracy_values.append(test_accuracy)
        train_accuracy_values.append(train_accuracy)

    lr_labels = [r"Optimized"]
    lr_labels += [r"$\eta={0:.2f}$".format(lr) for lr in learning_rates[1:]]
    plot_accuracy_scores(lmbdas, train_accuracy_values, test_accuracy_values,
                         lr_labels, "accuracy_learning_rate_scores",
                         r"$\lambda$", r"Accuracy")

def run_lambda_momentum(X_train, X_test, y_train, y_test,
                        lmbdas=None,
                        momentums=None,
                        try_get_pickle=True,
                        figure_folder="../fig",
                        **kwargs):

    param_dict = {"lmbdas": lmbdas,
                  "sklearn": True,
                  "figure_folder": figure_folder,
                  "momentums": momentums}
    param_dict.update(kwargs)
    print_parameters(**param_dict)

    test_accuracy_values = []
    train_accuracy_values = []

    for momentum in momentums:
        print ("Momentum: ", momentum)
        pickle_fname = ("lambda_mom_accuracy_penalty{}_actsigmoid"
                        "_solverlr-gd_lrinverse_mom{}_tol1e-06."
                        "pkl".format(kwargs["penalty"], str(momentum)))

        if os.path.isfile(pickle_fname) and try_get_pickle:
            res_ = load_pickle(pickle_fname)
        else:
            res_ = logreg_core(X_train, X_test, y_train, y_test,
                               use_sk_learn=False,
                               lmbdas=lmbdas,
                               momentum=momentum,
                               store_pickle=True,
                               pickle_fname=pickle_fname,
                               **kwargs)

        train_accuracy, test_accuracy, critical_accuracy = res_

        test_accuracy_values.append(test_accuracy)
        train_accuracy_values.append(train_accuracy)

    plot_accuracy_scores(lmbdas, train_accuracy_values, test_accuracy_values,
                         [r"$\gamma={0:.1e}$".format(m) for m in momentums],
                         "accuracy_momentum_scores", r"$\lambda$",
                         r"Accuracy")

def run_sk_comparison(X_train, X_test, y_train, y_test,
                      lmbdas=None,
                      penalty=None,
                      activation=None,
                      solver=None,
                      learning_rate=None,
                      momentum=None,
                      mini_batch_size=None,
                      max_iter=None,
                      tolerance=None,
                      verbose=False,
                      try_get_pickle=True,
                      figure_folder="../fig"):
    """Runs a comparison between sk learn and our method."""

    param_dict = {"lmbdas": lmbdas,
                  "sklearn": True,
                  "penalty": penalty,
                  "activation": activation,
                  "solver": solver,
                  "learning_rate": learning_rate,
                  "momentum": momentum,
                  "mini_batch_size": mini_batch_size,
                  "max_iter": max_iter,
                  "tolerance": tolerance}
    print_parameters(**param_dict)

    pickle_fname = ("sk_comparison_accuracy_sklearn_penalty{}_actsigmoid"
                    "_solverlr-gd_lrinverse_mom0.0_tol1e-06"
                    ".pkl".format(penalty))

    if os.path.isfile(pickle_fname) and try_get_pickle:
        res_ = load_pickle(pickle_fname)
    else:
        res_ = logreg_core(X_train, X_test, y_train, y_test,
                           use_sk_learn=True,
                           lmbdas=lmbdas,
                           penalty=penalty,
                           activation=activation,
                           solver=solver,
                           learning_rate=learning_rate,
                           momentum=momentum,
                           mini_batch_size=mini_batch_size,
                           max_iter=max_iter,
                           tolerance=tolerance,
                           store_pickle=True,
                           pickle_fname=pickle_fname,
                           verbose=verbose)

    # Retrieves results
    train_accuracy, test_accuracy, critical_accuracy, \
        train_accuracy_SK, test_accuracy_SK, critical_accuracy_SK, \
        train_accuracy_SGD, test_accuracy_SGD, critical_accuracy_SGD = res_

    print('Mean accuracy: train, test')

    print("HomeMade: {0:0.4f} +/- {1:0.2f}, {2:0.4f} +/- {3:0.2f}".format(
        np.mean(train_accuracy), np.std(train_accuracy),
        np.mean(test_accuracy), np.std(test_accuracy)))

    print("SK: {0:0.4f} +/- {1:0.2f}, {2:0.4f} +/- {3:0.2f}".format(
        np.mean(train_accuracy_SK), np.std(train_accuracy_SK),
        np.mean(test_accuracy_SK), np.std(test_accuracy_SK)))

    print("SGD: {0:0.4f} +/- {1:0.2f}, {2:0.4f} +/- {3:0.2f}".format(
        np.mean(train_accuracy_SGD), np.std(train_accuracy_SGD),
        np.mean(test_accuracy_SGD), np.std(test_accuracy_SGD)))

    plot_accuracy_comparison(lmbdas, train_accuracy, test_accuracy,
                             train_accuracy_SK, test_accuracy_SK,
                             train_accuracy_SGD, test_accuracy_SGD,
                             "logistic_accuracy_sklearn_comparison",
                             figure_folder)

def run_logreg_fit(X_train, X_test, y_train, y_test, use_sk_learn=False,
                   lmbdas=[None], penalties=[None],
                   activations=[None], solvers=[None],
                   learning_rates=[None],
                   momentums=[None], mini_batch_sizes=[None],
                   max_iter=None, tolerances=[None], store_pickle=False,
                   verbose=False):
    """Method for retrieveing data for given lists of hyperparameters

    Args:
        X_train (ndarray)
        X_test (ndarray)
        y_train (ndarray)
        y_test (ndarray)
        use_sk_learn (bool): if we are to use SK-learn. Default is False.
        lmbdas (list(float)): list of lmbdas.
        penalties (list(str): list of penalty types. Choices: l1, l2,
            elastic_net.
        activations (list(str)): list of activation functions.
        solvers (list(str)): list of solver functions.
        learning_rates (list(str|float)): list of learning rates.
            Options: float, inverse.
        momentums (list(float)): momentum strengths, list of floats.
        mini_batch_sizes (list(float)): list of minibatch sizes.
        max_iter (int): maximum number of iterations
        tolerances (list(floats)): list of tolerances
        store_pickle (bool): store results as pickles
        verbose (bool): more verbose output. Default is False.

    Returns:
        Dictionary with logreg accuracy scores and times
        SK-learn dictionary with accuracy scores and times
        SGD-SK-learn dictionary with accuracy scores and times
    """

    for penalty in penalties:
        for activation in activations:
            for solver in solvers:
                for lr in learning_rates:
                    for momentum in momentums:
                        for mb in mini_batch_sizes:
                            for tol in tolerances:
                                logreg_core(X_train, X_test,
                                            y_train, y_test,
                                            use_sk_learn=use_sk_learn,
                                            lmbdas=lmbdas,
                                            penalty=penalty,
                                            activation=activation,
                                            learning_rate=lr,
                                            momentum=momentum,
                                            mini_batch_size=mb,
                                            max_iter=max_iter,
                                            tolerance=tol,
                                            store_pickle=store_pickle,
                                            verbose=False)

def logreg_core(X_train, X_test, y_train, y_test, use_sk_learn=False,
                lmbdas=[None], penalty=None,
                activation=None, solver=None, learning_rate=None,
                momentum=None, mini_batch_size=None, max_iter=None,
                tolerance=None, store_pickle=False, pickle_fname=None,
                verbose=False):
    """Method for retrieveing data for given lists of hyperparameters

    Args:
        X_train (ndarray)
        X_test (ndarray)
        y_train (ndarray)
        y_test (ndarray)
        use_sk_learn (bool): if we are to use SK-learn. Default is False.
        lmbdas (float): list of lmbdas
        penalty (str): penalty type. Choices: l1, l2, elastic_net
        activation (str): activation function.
        solver(str): solver function.
        learning_rate (str|float): learning rate. Options: float, inverse
        momentum (float): momentum strength
        mini_batch_size (float): minibatch size
        tolerance (float): tolerance, at what point we cut off the parameter
            search.
        store_pickle (bool): saves output as pickle
        pickle_fname (str): pickle filename
        verbose (bool): more verbose output. Default is False

    Returns:
        Dictionary with logreg accuracy scores and times
        SK-learn dictionary with accuracy scores and times
        SGD-SK-learn dictionary with accuracy scores and times
    """

    # Sets up data arrays
    train_accuracy = np.zeros(lmbdas.shape, np.float64)
    test_accuracy = np.zeros(lmbdas.shape, np.float64)
    critical_accuracy = np.zeros(lmbdas.shape, np.float64)

    if use_sk_learn:
        train_accuracy_SK = np.zeros(lmbdas.shape, np.float64)
        test_accuracy_SK = np.zeros(lmbdas.shape, np.float64)
        critical_accuracy_SK = np.zeros(lmbdas.shape, np.float64)

        train_accuracy_SGD = np.zeros(lmbdas.shape, np.float64)
        test_accuracy_SGD = np.zeros(lmbdas.shape, np.float64)
        critical_accuracy_SGD = np.zeros(lmbdas.shape, np.float64)

    # Loops over regularisation strength
    for i, lmbda in enumerate(lmbdas):

        if verbose:
            print("")
            print("="*80)
            print("Lambda = ", lmbda)

        if use_sk_learn:
            # Sets up lambda for SK learn logreg methods
            if lmbda == 0.0:
                sk_lmbda = 1.0 / 10000000
            else:
                sk_lmbda = 1.0 / lmbda

            # Set sk penalty
            if penalty == "elastic_net":
                sk_penalty = "l2"
            else:
                sk_penalty = penalty

            # Define SK-learn logistic regressor
            logreg_SK = sk_model.LogisticRegression(
                penalty=sk_penalty,
                fit_intercept=False, C=sk_lmbda,
                max_iter=max_iter, tol=tolerance)

            # SK learn fit
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                logreg_SK.fit(cp.deepcopy(X_train), cp.deepcopy(y_train))
            if verbose:
                print("SK-learn method done")

            # Retrieve accuracy scores for SK-learn
            train_accuracy_SK[i] = logreg_SK.score(X_train, y_train)
            test_accuracy_SK[i] = logreg_SK.score(X_test, y_test)

            # Sets up learning rate for SK-learn SGD
            if learning_rate == "inverse":
                sk_learning_rate = "optimal"
                eta0 = 0.0
            else:
                sk_learning_rate = "constant"
                eta0 = learning_rate

            # Sets up regularisation for SK-learn SGD
            if penalty == "elastic_net":
                sk_penalty = "elasticnet"

            # define SGD-based logistic regression
            logreg_SGD = sk_model.SGDClassifier(loss="log",
                                                penalty=sk_penalty,
                                                alpha=lmbda, max_iter=max_iter,
                                                shuffle=True, random_state=1,
                                                eta0=eta0,
                                                learning_rate=sk_learning_rate)

            # fit training data
            logreg_SGD.fit(X_train, y_train)
            if verbose:
                print("SGA SK-learn method done")

            # check accuracy
            train_accuracy_SGD[i] = logreg_SGD.score(X_train, y_train)
            test_accuracy_SGD[i] = logreg_SGD.score(X_test, y_test)
            # critical_accuracy_SGD[i]=logreg_SGD.score(X_critical,Y_critical)

        # Our implementation of logistic regression
        log_reg = logreg.LogisticRegression(penalty=penalty, solver=solver,
                                            activation=activation,
                                            tol=tolerance,
                                            alpha=lmbda, momentum=momentum,
                                            mini_batch_size=mini_batch_size,
                                            max_iter=max_iter)

        # Fit training data
        log_reg.fit(cp.deepcopy(X_train), cp.deepcopy(y_train),
                    eta=learning_rate)
        if verbose:
            print("HomeMade method done")

        # Accuracy score for our implementation
        train_accuracy[i] = log_reg.score(X_train, y_train)
        test_accuracy[i] = log_reg.score(X_test, y_test)
        # critical_accuracy[i]=log_reg.score(X_critical,Y_critical)

        # Prints result from single lambda run
        if verbose:
            print('Accuracy scores: train, test, critical')
            print('HomeMade: {0:0.4f}, {1:0.4f}, {2:0.4f}'.format(
                train_accuracy[i], test_accuracy[i], critical_accuracy[i]))

        if use_sk_learn and verbose:
            print('SK: {0:0.4f}, {1:0.4f}, {2:0.4f}'.format(
                train_accuracy_SK[i], test_accuracy_SK[i],
                critical_accuracy_SK[i]))

            print('SGD: %0.4f, %0.4f, %0.4f' % (
                train_accuracy_SGD[i],
                test_accuracy_SGD[i],
                critical_accuracy_SGD[i]))

        # Prints iteration values
        if verbose:
            print("Finished computing {}/11 iterations".format(i+1))

    if verbose:
        print("")

    results = [train_accuracy, test_accuracy, critical_accuracy]

    if use_sk_learn:
        results += [train_accuracy_SK, test_accuracy_SK,
                    critical_accuracy_SK, train_accuracy_SGD,
                    test_accuracy_SGD, critical_accuracy_SGD]

    if store_pickle:
        if use_sk_learn:
            sk_str = "sklearn_"
        else:
            sk_str = ""

        if isinstance(pickle_fname, type(None)):
            fname = ("accuracy_{}penalty{}_act{}_solver{}_lr{}_mom{}_tol{}."
                     "pkl".format(sk_str, penalty, activation, solver,
                                  str(learning_rate), str(momentum),
                                  str(tolerance)))
        else:
            fname = pickle_fname
        save_pickle(fname, results)

    return results

def plot_accuracy_comparison(lmbdas, train_accuracy, test_accuracy,
                             train_accuracy_SK, test_accuracy_SK,
                             train_accuracy_SGD, test_accuracy_SGD,
                             figname, figure_folder):
    """plot accuracy against regularisation strength."""
    fig = plt.figure()

    ax1 = fig.add_subplot(111)

    ax1.semilogx(lmbdas, train_accuracy,
                 marker="o", ls=(0, (3, 1, 1, 1)),  # Densely dashdotted
                 color="#7570b3",
                 label='Implemented train')
    ax1.semilogx(lmbdas, test_accuracy,
                 marker="x", ls=(0, (3, 1, 1, 1)),  # Densely dashdotted
                 color="#7570b3",
                 label='Implemented test')

    ax1.semilogx(lmbdas, train_accuracy_SK,
                 marker="o", ls=(0, (5, 1)),  # Densely dashed
                 color="#1b9e77",
                 label='SK-Learn train')
    ax1.semilogx(lmbdas, test_accuracy_SK,
                 marker="x", ls=(0, (5, 1)),  # Densely dashed
                 color="#1b9e77",
                 label='SK-Learn test')

    ax1.semilogx(lmbdas, train_accuracy_SGD,
                 marker="o", ls=(0, (3, 5, 1, 5)),  # Dashdotted
                 color="#d95f02",
                 label='SK-Learn SGD train')
    ax1.semilogx(lmbdas, test_accuracy_SGD,
                 marker="x", ls=(0, (3, 5, 1, 5)),  # Dashdotted
                 color="#d95f02",
                 label='SK-Learn SGD test')

    ax1.set_xlabel(r'$\lambda$')
    ax1.set_ylabel(r'$\mathrm{Accuracy}$')

    ax1.grid(True)
    ax1.legend(fontsize=8)

    figure_path = os.path.join(figure_folder, "{}.pdf".format(figname))
    fig.savefig(figure_path)
    print("Figure saved at {}".format(figure_path))
    plt.close(fig)

if __name__ == '__main__':
    task1c()


text:
import numpy as np
from flask import Flask, request, jsonify
app = Flask(__name__)

@app.route("/")
def help():
    return "Give the number of sides the die should have.\n"

@app.route("/<int:sides>")
def roll_die(sides):
    return str(np.random.randint(1,sides+1))+'\n'

@app.route("/json/<int:sides>")
def roll_die_json(sides):
    return jsonify({'sides': sides,'roll': np.random.randint(1,sides+1)})

max_stars_repo_name: perrindesign/data-science-class

max_issues_repo_name: perrindesign/data-science-class

max_forks_repo_name: perrindesign/data-science-class