text:
using KernelDensity
using PyPlot
using PyCall
using Printf
import Random.randperm
@pyimport matplotlib.patches as patches
@pyimport matplotlib.legend_handler as lh
@pyimport matplotlib.lines as lines

function scale_outer_intervals(intervals, scale=1.0)
    int_min = min([int[1] for int in intervals]...)
    int_max = max([int[2] for int in intervals]...)
    int_scaled_half = (int_max - int_min) * scale / 2.0
    int_mid = (int_min + int_max) / 2.0
    int_mid - int_scaled_half, int_mid + int_scaled_half
end

function my_linspace(bounds::Tuple{Float64, Float64}, length::Int64)
    range(bounds[1], stop=bounds[2], length=length)
end

# This script relies on data (emu_out, sim_out) that is generated by running smc-abc-script.jl
function plot_emulation_vs_simulation(emu_out, sim_out, last_pop=-1)
            if last_pop == -1
                last_pop = length(sim_out.population)
            end
            grid_size = emu_out.n_params
            emu_handle = nothing
            sim_handle = nothing
            kernel_bandwidth_scale = 0.09
            bounds_scale = 1.2
            # population_colors = ["#B1E9DE", "#007731"]
            population_colors = ["#DDF4F7", "#B1E9DE", "#63D3BB", "#00BD8B", "#007731"]

            contour_colors = ["white", "#FFE9EC", "#FFBBC5", "#FF8B9C", "#FF5D75", "#FF2F4E", "#D0001F", "#A20018", "#990017", "#800013"]
            simulation_color = "#08519c"
            emulation_color = "#ff6600"
            sim_count = 30
            h_sim_joint = nothing
            h_emu_pops = nothing
            h_sim_marg = nothing
            h_emu_marg = nothing
            h_true = nothing
            for i in 1:grid_size
                for j in 1:grid_size
                    subplot2grid((grid_size, grid_size), (i - 1, j - 1))
                    if j < i
                        x_data_emu = emu_out.population[last_pop][:,j]
                        y_data_emu = emu_out.population[last_pop][:,i]
                        x_data_sim = sim_out.population[last_pop][:,j]
                        y_data_sim = sim_out.population[last_pop][:,i]
                        sim_size = size(sim_out.population[last_pop], 1)
                        if sim_size > sim_count
                            idx = randperm(sim_size)[1:sim_count]
                            x_data_sim = x_data_sim[idx]
                            y_data_sim = y_data_sim[idx]
                        end
                        x_extr_emu = extrema(x_data_emu)
                        y_extr_emu = extrema(y_data_emu)
                        x_extr_sim = extrema(x_data_sim)
                        y_extr_sim = extrema(y_data_sim)

                        xlim(scale_outer_intervals([x_extr_emu, x_extr_sim], bounds_scale))
                        ylim(scale_outer_intervals([y_extr_emu, y_extr_sim], bounds_scale))

                        bandwidth = (-kernel_bandwidth_scale * -(x_extr_emu...), -kernel_bandwidth_scale * -(y_extr_emu...))
                        kde_joint = kde((x_data_emu, y_data_emu), bandwidth=bandwidth)
                        contour_x = my_linspace(scale_outer_intervals([x_extr_emu], bounds_scale), 100)
                        contour_y = my_linspace(scale_outer_intervals([y_extr_emu], bounds_scale), 100)
                        contour_z = pdf(kde_joint, contour_x, contour_y)
                        contourf(contour_x, contour_y, contour_z, 8, colors=contour_colors, zorder=1)

                        h_sim_joint = scatter(x_data_sim, y_data_sim, marker="x", color=simulation_color, zorder=2)
                        true_param_x = true_params[param_indices[j]]
                        true_param_y = true_params[param_indices[i]]
                        plot([true_param_x, true_param_x], [0, true_param_y], color=:black, linestyle=:dashed, linewidth=0.5)
                        plot([0, true_param_x], [true_param_y, true_param_y], color=:black, linestyle=:dashed, linewidth=0.5)
                    elseif j > i
                        h_emu_pops = []
                        for iter_idx in 1:last_pop
                            h_emu_pop = scatter(emu_out.population[iter_idx][:,j], emu_out.population[iter_idx][:,i],
                                color=population_colors[iter_idx], zorder=iter_idx, marker=".")
                            push!(h_emu_pops, h_emu_pop)
                        end
                        xlims = scale_outer_intervals([extrema(pop[:, j]) for pop in emu_out.population])
                        xlim(xlims)
                        ylims = scale_outer_intervals([extrema(pop[:, i]) for pop in emu_out.population])
                        ylim(ylims)
                        true_param_x = true_params[param_indices[j]]
                        true_param_y = true_params[param_indices[i]]
                        plot([true_param_x, true_param_x], [ylims[1], true_param_y], color=:black, linestyle=:dashed, linewidth=0.5, zorder=length(emu_out.population)+1)
                        plot([xlims[1], true_param_x], [true_param_y, true_param_y], color=:black, linestyle=:dashed, linewidth=0.5, zorder=length(emu_out.population)+1)
                    else#if i == j
                        emu_data = emu_out.population[last_pop][:,i]
                        sim_data = sim_out.population[last_pop][:,i]
                        extr_emu = extrema(emu_data)
                        extr_sim = extrema(sim_data)
                        kde_emu = kde(emu_data, bandwidth=-kernel_bandwidth_scale * -(extr_emu...))
                        kde_sim = kde(sim_data, bandwidth=-kernel_bandwidth_scale * -(extr_sim...))
                        x_bounds = scale_outer_intervals([extr_emu, extr_sim], bounds_scale)
                        xlim(x_bounds)
                        x_plot = my_linspace(x_bounds, 100)
                        y_emu_plot = pdf(kde_emu,x_plot)
                        y_sim_plot = pdf(kde_sim,x_plot)
                        h_emu_marg = plot(x_plot, y_emu_plot, color=emulation_color)
                        h_sim_marg = plot(x_plot, y_sim_plot, color=simulation_color)
                        yticks([])
                        true_param = true_params[param_indices[i]]
                        max_pdf = max(pdf(kde_emu, true_param), pdf(kde_sim, true_param))
                        h_true = plot([true_param, true_param], [0, max_pdf], color=:black, linestyle=:dashed, linewidth=0.5)
                    end
                end
            end
            # https://matplotlib.org/gallery/text_labels_and_annotations/legend_demo.html
            h_emu_pops = Tuple(h_emu_pops[max(1, length(h_emu_pops)-3):end])
            h_emu = (h_emu_marg[1],
                lines.Line2D([], [], linestyle="None", marker="s", color=contour_colors[7],
                markerfacecoloralt=contour_colors[5], fillstyle="right", markeredgewidth=0.0))
            h_sim = (h_sim_marg[1], h_sim_joint)
            figlegend([h_sim; h_emu_pops; h_emu; h_true],
                ["Simulation result", "Emulated populations", "Emulation result", "True value"],
                ncol=2, loc="lower center",
                handler_map=Dict(h_sim=>lh.HandlerTuple(ndivide=nothing),
                        h_emu_pops=>lh.HandlerTuple(ndivide=nothing),
                        h_emu=>lh.HandlerTuple(ndivide=nothing))
            )
end

ion()
fig = figure()
ioff()
plot_emulation_vs_simulation(emu_out, sim_out)
subplots_adjust(
left    =  0.08,
bottom  =  0.2,
right   =  0.96,
top     =  0.97,
wspace  =  0.35,
hspace  =  0.26)
show(fig)

# savefig("/bmm/home/et517/gaussian_processes/GaussianProcesses/gp_results/retraining_1_2_3_7/figure_retraining_prev_pop_thresh_1.0_var_100.pdf")

function plot_ode_ref_emu_sim(emu_out, sim_out, idxs)
                        kernel_bandwidth_scale = 0.09
                        bounds_scale=1.2
                        param_sim = zeros(length(idxs))
                        param_emu = zeros(length(idxs))
                        for i in 1:length(idxs)
                            idx = idxs[i]
                            emu_data = emu_out.population[end][:,idx]
                            sim_data = sim_out.population[end][:,idx]
                            extr_emu = extrema(emu_data)
                            extr_sim = extrema(sim_data)
                            kde_emu = kde(emu_data, bandwidth=-kernel_bandwidth_scale * -(extr_emu...))
                            kde_sim = kde(sim_data, bandwidth=-kernel_bandwidth_scale * -(extr_sim...))
                            x_bounds = scale_outer_intervals([extr_emu, extr_sim], bounds_scale)
                            x_plot = my_linspace(x_bounds, 100)
                            param_sim[i] = x_plot[findmax(pdf(kde_sim,x_plot))[2]]
                            param_emu[i] = x_plot[findmax(pdf(kde_emu,x_plot))[2]]
                        end
                        dst_sim = simulate_distance(param_indices[idxs], param_sim)
                        dst_emu = simulate_distance(param_indices[idxs], param_emu)
                        plt_ref = Plots.plot(reference_data,
                            # title=LaTeXString("Reference; \$\\theta_$(param_indices[idx]) = $(true_params[param_indices[idx]])\$"));
                            title="Reference");
                        plt_emu = Plots.plot(simulator_function(param_indices[idxs], param_emu),
                            # title=LaTeXString("Emulation; \$\\theta_$(param_indices[idx]) = $(@sprintf("%.2f", param_emu)); d=$(@sprintf("%.2f", dst_emu))\$"));
                            title="Emulation; d=$(@sprintf("%.2f", dst_emu))");
                        plt_sim = Plots.plot(simulator_function(param_indices[idxs], param_sim),
                            # title=LaTeXString("Simulation; \$\\theta_$(param_indices[idx]) = $(@sprintf("%.2f", param_sim)); d=$(@sprintf("%.2f", dst_sim))\$"));
                            title="Simulation; d=$(@sprintf("%.2f", dst_sim))");
                        Plots.plot(plt_ref, plt_sim, plt_emu, layout=(3,1))
end
# import Plots
# Plots.pyplot()
# using LaTeXStrings
# plot_ode_ref_emu_sim(emu_out, sim_out, 1:4)
# Plots.savefig("/bmm/home/et517/gaussian_processes/GaussianProcesses/gp_results/retraining_1_2_3_7/diff_param_7_inc_var_1_eps.pdf")


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    SlackBridge{T, F, G}

The `SlackBridge` converts an objective function of type `G` into a
[`MOI.VariableIndex`](@ref) objective by creating a slack variable and a
`F`-in-[`MOI.LessThan`](@ref) constraint for minimization or
`F`-in-[`MOI.LessThan`](@ref) constraint for maximization where `F` is
`MOI.Utilities.promote_operation(-, T, G, MOI.VariableIndex}`.
Note that when using this bridge, changing the optimization sense
is not supported. Set the sense to `MOI.FEASIBILITY_SENSE` first
to delete the bridge in order to change the sense, then re-add the objective.
"""
struct SlackBridge{
    T,
    F<:MOI.AbstractScalarFunction,
    G<:MOI.AbstractScalarFunction,
} <: AbstractBridge
    slack::MOI.VariableIndex
    constraint::Union{
        MOI.ConstraintIndex{F,MOI.LessThan{T}},
        MOI.ConstraintIndex{F,MOI.GreaterThan{T}},
    }
end

const Slack{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{SlackBridge{T},OT}

function bridge_objective(
    ::Type{SlackBridge{T,F,G}},
    model::MOI.ModelLike,
    func::G,
) where {T,F,G<:MOI.AbstractScalarFunction}
    slack = MOI.add_variable(model)
    f = MOIU.operate(-, T, func, slack)
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
        set = MOI.LessThan(zero(T))
    elseif MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        set = MOI.GreaterThan(zero(T))
    else
        error(
            "Set `MOI.ObjectiveSense` before `MOI.ObjectiveFunction` when",
            " using `MOI.Bridges.Objective.SlackBridge`.",
        )
    end
    constraint = MOIU.normalize_and_add_constraint(model, f, set)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), slack)
    return SlackBridge{T,F,G}(slack, constraint)
end

function supports_objective_function(
    ::Type{<:SlackBridge},
    ::Type{MOI.VariableIndex},
)
    return false
end

function supports_objective_function(
    ::Type{<:SlackBridge},
    ::Type{<:MOI.AbstractScalarFunction},
)
    return true
end

MOIB.added_constrained_variable_types(::Type{<:SlackBridge}) = Tuple{Type}[]

function MOIB.added_constraint_types(::Type{<:SlackBridge{T,F}}) where {T,F}
    return Tuple{Type,Type}[(F, MOI.GreaterThan{T}), (F, MOI.LessThan{T})]
end

function MOIB.set_objective_function_type(::Type{<:SlackBridge})
    return MOI.VariableIndex
end

function concrete_bridge_type(
    ::Type{<:SlackBridge{T}},
    G::Type{<:MOI.AbstractScalarFunction},
) where {T}
    F = MOIU.promote_operation(-, T, G, MOI.VariableIndex)
    return SlackBridge{T,F,G}
end

# Attributes, Bridge acting as a model
MOI.get(::SlackBridge, ::MOI.NumberOfVariables)::Int64 = 1

function MOI.get(bridge::SlackBridge, ::MOI.ListOfVariableIndices)
    return [bridge.slack]
end

function MOI.get(
    bridge::SlackBridge{T,F},
    ::MOI.NumberOfConstraints{F,S},
)::Int64 where {T,F,S<:Union{MOI.GreaterThan{T},MOI.LessThan{T}}}
    return bridge.constraint isa MOI.ConstraintIndex{F,S} ? 1 : 0
end

function MOI.get(
    bridge::SlackBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,S},
) where {T,F,S<:Union{MOI.GreaterThan{T},MOI.LessThan{T}}}
    if bridge.constraint isa MOI.ConstraintIndex{F,S}
        return [bridge.constraint]
    else
        return MOI.ConstraintIndex{F,S}[]
    end
end

function MOI.delete(model::MOI.ModelLike, bridge::SlackBridge)
    MOI.delete(model, bridge.constraint)
    MOI.delete(model, bridge.slack)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOIB.ObjectiveFunctionValue{G},
    bridge::SlackBridge{T,F,G},
) where {T,F,G}
    slack = MOI.get(
        model,
        MOIB.ObjectiveFunctionValue{MOI.VariableIndex}(attr.result_index),
    )
    # There may be a gap between the value of the original objective `g` and the
    # value of `bridge.slack`. Since `bridge.constraint` is `g - slack`, we can
    # get the value of the original objective by summing the value of `slack`
    # with the `ConstraintPrimal` of the constraint.
    obj_slack_constant =
        MOI.get(model, MOI.ConstraintPrimal(), bridge.constraint)
    # The constant was moved to the set as it is a scalar constraint.
    constant =
        MOI.constant(MOI.get(model, MOI.ConstraintSet(), bridge.constraint))
    return obj_slack_constant + slack - constant
end

_constant_term(set::MOI.LessThan) = set.upper
_constant_term(set::MOI.GreaterThan) = set.lower

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ObjectiveFunction{G},
    bridge::SlackBridge{T,F,G},
) where {T,F,G<:MOI.AbstractScalarFunction}
    func = MOI.get(model, MOI.ConstraintFunction(), bridge.constraint)
    set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
    f = MOIU.operate(-, T, func, _constant_term(set))
    return MOIU.convert_approx(G, MOIU.remove_variable(f, bridge.slack))
end

ext: jl

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using DelimitedFiles

# Unspecified scientific constancs
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\constants.jl")

# Table of masses of nuclides
masses = readdlm("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\masses.txt")
atomic_weight = Dict()
for i = 1:size(masses)[1]
    atomic_weight[convert(Int64, masses[i, 1])] = masses[i, 2]
end
masses = 0

# Table of natural abundances of nuclides
abundances = readdlm("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\abundances.txt")

# Table of atomic numbers and symbols
symbols = readdlm("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\atomic_symbols.txt")

function to_ZA(nuclide_str::String)
    #=
    Converts a string representing a nuclide into an Int64 in ZA notation

    :param nuclide_str: string representing a nuclide. "92235" or "U-235"
    :return ZA: Int64 representing the nuclide in ZA notation. 1000*Z + A
    =#
    if isnumeric(nuclide_str[1])
        ZA =  parse(Int64, nuclide_str)
    else
        hyphen = findfirst("-", nuclide_str)[]
        symb = nuclide_str[1 : hyphen-1]
        Z = findfirst(x -> x == symb, symbols[:,2])
        A_str = nuclide_str[hyphen+1 : end]
        if A_str == "nat"
            A = 0
        else
            A = parse(Int64, A_str)
        end
        ZA = 1000*Z + A
    end
    return ZA
end

function natural_nuclide(ZA::Int)
    #=
    Gather information needed to implement natural isotopes

    :param ZA: Int representing the nuclide in ZA notation. 1000*Z + A
    :return natural_weight: Float, average mass of atom given natural abundances
    :return nuclides: list of Int for ZA of nuclides present
    :fracs: List of Float for fractions of each nuclide naturally occuring
    =#
    Z = ZA ÷ 1000
    indices = findall(x -> x == Z, abundances[:, 1])
    As = abundances[indices, 2]
    nuclides = convert.(Int64, As.+Z*1000)
    fracs = abundances[indices, 4]
    weights = abundances[indices, 3]
    natural_weight = sum(fracs.*weights)
    mass_fracs = fracs.*weights / natural_weight
    return natural_weight, nuclides, fracs
end

# Add natural isotope atomic weights to table
for Z in 1:92
    atomic_weight[Z*1000] = natural_nuclide(Z*1000)[1]
end

mutable struct Material
    #=
    Structure to contain information about a material

    nuclides: List of Int to represent ZA of each nuclide in material
    atomic_dens: List of Float for the number denisity of each nuclide (1/cm^3)
    temperature: Float of temperature of material in K
    =#
    nuclides::Vector{Int64}
    atomic_dens::Vector{Float64}
    temperature::Float64
    color::Vector{Int64}
end

function define_Material(;
    density::Float64 = 0.0,
    number_density::Float64 = 0.0,
    nuclides::Vector{String} = String[],
    mass_fracs::Vector{Float64} = Float64[],
    mass_dens::Vector{Float64} = Float64[],
    atomic_fracs::Vector{Float64} = Float64[],
    atomic_dens::Vector{Float64} = Float64[],
    temperature::Float64 = 293.6,
    normalize_fracs::Bool = true,
    color::Vector{Int64} = [-1, -1, -1])
    #=
    This function allows for a variety of methods to define a material. As long
    as the provided information is sufficient to deduce the namber densities of
    each nuclide, the function will work unless there is conflicting
    information.

    :param density: density of material in g/cm^3
    :param number_density: number density of material in atoms/cm^3
    :param nuclides: List of nuclides present. Format may be "92235" or "U-235".
    Natural isotopes are indicated as "92000" or "U-nat".
    :param mass_fracs: fraction of mass from each nuclide
    :param mass_dens: density of each nuclide in g/cm^3
    :param atomic_fracs: fraction of total atoms from each nuclide
    :param atomic_dens: List of the number denisity of each nuclide (1/cm^3)
    :param temperature: temperature of material in K
    :param normalize_fracs: if false and fractions don't add to 1, the total
    number density or density will be inconsistent with the input. Setting to
    true will maintain nuclide ratios and the input total density or number
    density.
    :return material: Material struct based on the inputs
    =#
     if sum(length.([mass_fracs, mass_dens, atomic_fracs, atomic_dens]).>0) > 1
         throw(ArgumentError("Confliciting information in define_Material. Only
         one of mass_fracs, mass dens, atomic_fracs, or atomic_dens can be
         provided."))
     end

     nuclides = to_ZA.(nuclides)

     if length(mass_fracs) > 0

         if length(mass_fracs) != length(nuclides)
             throw(ArgumentError("Mismatch in define_Material. Number of nuclides
             and number of mass fractions provided must match."))
         end

         if !(density > 0)
             throw(ArgumentError("Incomplete information in define_Material. If
             using mass fractions, total density must be provided."))
         end

         if normalize_fracs
             mass_dens = mass_fracs * density / sum(mass_fracs)
         else
             mass_dens = mass_fracs * density
         end
     end

     if length(mass_dens) > 0

         if length(mass_dens) != length(nuclides)
             throw(ArgumentError("Mismatch in define_Material. Number of nuclides
             and number of mass densities provided must match."))
         end

         atomic_weights = [atomic_weight[nuclide] for nuclide in nuclides]
         atomic_dens = N_A * mass_dens ./ atomic_weights
     end

     if length(atomic_fracs) > 0

         if length(atomic_fracs) != length(nuclides)
             throw(ArgumentError("Mismatch in define_Material. Number of nuclides
             and number of atomic_dens fractions provided must match."))
         end

         if !(number_density > 0) & !(density > 0)
             throw(ArgumentError("Incomplete information in define_Material. If
             using atomic fractions, total number density or density must be provided."))
         end

         if (number_density > 0) & (density > 0)
             throw(ArgumentError("Conflicting information in define_Material. Only
             one of density and number density can be provided."))
         end

         if density > 0
             atomic_weights = [atomic_weight[nuclide] for nuclide in nuclides]
             number_density = sum(atomic_fracs) * N_A * density / sum(atomic_weights.*atomic_fracs)
         end

         if normalize_fracs
             atomic_dens = number_density * atomic_fracs / sum(atomic_fracs)
         else
             atomic_dens = number_density * atomic_fracs
         end
     end

     # Replace natural nuclides with their constituents
     full_nuclides = []
     full_atomic_dens = []
     for i = 1:length(nuclides)
         if nuclides[i] % 1000 == 0
             AW, ZA, fracs = natural_nuclide(nuclides[i])
             for item in ZA
                 push!(full_nuclides, item)
             end
             for item in fracs
                 push!(full_atomic_dens, atomic_dens[i]*item)
             end
         else
             push!(full_nuclides, nuclides[i])
             push!(full_atomic_dens, atomic_dens[i])
         end
     end

    order = sortperm(full_nuclides)

    return Material(full_nuclides[order],
                    full_atomic_dens[order],
                    temperature,
                    color)
end

function density(material::Material)
    #=
    Returns the density of a material in g/cm^3

    :param material: Material struct
    =#
    dens = 0
    for i = 1:length(material.nuclides)
        dens += material.atomic_dens[i] * atomic_weight[material.nuclides[i]] / N_A
    end
    return dens
end

function mass_frac(material::Material, nuclide)
    #=
    Returns the fraction of a material's mass from a given nuclide

    :param material: Material struct
    :param nuclide: Nuclide identifier String or ZA Int
    =#
    if typeof(nuclide) == String
        nuclide = to_ZA(nuclide)
    end
    nuclide_dens = 0
    dens = 0
    for i = 1:length(material.nuclides)
        partial_dens = material.atomic_dens[i] * atomic_weight[material.nuclides[i]] / N_A
        dens += partial_dens
        if material.nuclides[i] == nuclide
            nuclide_dens += partial_dens
        end
    end
    return nuclide_dens / dens
end

function mass_frac_elem(material::Material, element)
    #=
    Returns the fraction of a material's mass from a given element

    :param material: Material struct
    :param element: String or Int to identify element. "U" or 92.
    =#
    if typeof(element) == String
        element = findfirst(x -> x == element, symbols[:,2])
    end
    elem_dens = 0
    dens = 0
    for i = 1:length(material.nuclides)
        partial_dens = material.atomic_dens[i] * atomic_weight[material.nuclides[i]] / N_A
        dens += partial_dens
        if material.nuclides[i] ÷ 1000 == element
            elem_dens += partial_dens
        end
    end
    return elem_dens / dens
end

function mass_dens(material::Material, nuclide)
    #=
    Returns the density of a nuclide in a material in g/cm^3

    :param material: Material struct
    :param nuclide: Nuclide identifier String or ZA Int
    =#
    if typeof(nuclide) == String
        nuclide = to_ZA(nuclide)
    end
    i = findfirst(x -> x == nuclide, material.nuclides)
    if i == nothing
        return 0
    else
        return atomic_weight[nuclide] * material.atomic_dens[i] / N_A
    end
end

function mass_dens_elem(material::Material, element)
    #=
    Returns the density of an element in a material in g/cm^3

    :param material: Material struct
    :param element: String or Int to identify element. "U" or 92.
    =#
    if typeof(element) == String
        element = findfirst(x -> x == element, symbols[:,2])
    end
    elem_dens = 0
    for i = 1:length(material.nuclides)
        if material.nuclides[i] ÷ 1000 == element
            elem_dens += material.atomic_dens[i] * atomic_weight[material.nuclides[i]] / N_A
        end
    end
    return elem_dens
end

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
using JuMP
using ConstraintSolver
using JSON

CS = ConstraintSolver
include("../../test/sudoku_fcts.jl")

function parseJSON(json_sums)
    sums = []
    for s in json_sums
        indices = Tuple[]
        for ind in s["indices"]
            push!(indices, tuple(ind...))
        end

        push!(sums, (result=s["result"], indices=indices))
    end
    return sums
end

function solve_all(filenames; benchmark=false, single_times=true)
    ct = time()
    for (i,filename) in enumerate(filenames)
        sums = parseJSON(JSON.parsefile("data/$(filename)"))

        m = Model(with_optimizer(CS.Optimizer))
        @variable(m, 1 <= x[1:9,1:9] <= 9, Int)

        for s in sums
            @constraint(m, sum([x[ind...] for ind in s.indices]) == s.result)
            # @constraint(m, [x[ind...] for ind in s.indices] in CS.AllDifferentSet(length(s.indices)))
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        if single_times
            GC.enable(false)
            t = time()
            optimize!(m)
            status = JuMP.termination_status(m)
            t = time()-t
            GC.enable(true)
            println(i-1,", ", t)
        else
            GC.enable(false)
            optimize!(m)
            status = JuMP.termination_status(m)
            GC.enable(true)
        end
        if !benchmark
            @show JuMP.backend(m).optimizer.model.inner.info
            println("Status: ", status)
            @assert jump_fulfills_sudoku_constr(JuMP.value.(x))
        end
        com = nothing
        GC.gc()
    end
    println("")
    tt = time()-ct
    println("total time: ", tt)
    println("avg: ", tt/length(filenames))
end

function main(; benchmark=false, single_times=true)
    solve_all(["niallsudoku_5500", "niallsudoku_5501", "niallsudoku_5502", "niallsudoku_5503"]; benchmark=benchmark, single_times=single_times)
    # solve_all(from_file("hardest.txt"), "hardest")
end


text:
using SumProductNetworks
using Test

@testset "utility functions" begin
    q = rand(2)
    @test sum(projectToPositiveSimplex!(q)) == 1.0

    # check corner cases
    @test sum(projectToPositiveSimplex!(zeros(2))) == 1.0

    # check corner cases
    @test projectToPositiveSimplex!([0.5, 0.5]) == [0.5, 0.5]
end


text:
struct SMaxHandler{N} end

gethandler(::Bool,
           ::Type{typeof(max)},
           x::Type{<:SBitstreamLike},
           y::Type{<:SBitstreamLike},
           zs::Type{<:SBitstreamLike}...) =
    SMaxHandler{length(zs) + 2}()
init_state(::SMaxHandler) = (id = 0,)

function (handler::SMaxHandler{N})(buffer, netlist, state, inputs, outputs) where N
    @assert length(inputs) == N "Cannot apply SMaxHandler{$N} to length(inputs) == $(length(inputs))"

    # update netlist with inputs
    foreach(inputs) do input
        setsigned!(netlist, input, true)
    end

    # compute broadcast naming
    names = handle_broadcast_name(name.(inputs), netsize.(inputs), netsize(outputs[1]))

    # update netlist with output
    setsigned!(netlist, outputs[1], true)

    if handler.broadcasted
        num_elements = join(netsize(outputs[1]), "*")
        push!(netlist, Net(name = "max_$(N)_$(state.id)_inputs", size = (N, netsize(outputs[1])...), signed = true))
        write(buffer, """
            $stdcomment
            // BEGIN max_$(N)_$(state.id)
            """)
        write_bcast_instantiation(buffer, "max_$(N)_$(state.id)", netsize(outputs[1]), """
            assign max_$(N)_$(state.id)_inputs_p[((max_$(N)_$(state.id)_i + 1)*$N - 1) -: $N] = {$(join(map(fname -> fname("_p"), names), "[max_$(N)_$(state.id)_i], "))[max_$(N)_$(state.id)_i]};
            assign max_$(N)_$(state.id)_inputs_m[((max_$(N)_$(state.id)_i + 1)*$N - 1) -: $N] = {$(join(map(fname -> fname("_m"), names), "[max_$(N)_$(state.id)_i], "))[max_$(N)_$(state.id)_i]};

            stoch_signed_nmax #(
                    .NUM_INPUTS($N)
                ) max_$(N)_$(state.id) (
                    .CLK(CLK),
                    .nRST(nRST),
                    .as_p(max_$(N)_$(state.id)_inputs_p[((max_$(N)_$(state.id)_i + 1)*$N - 1) -: $N]),
                    .as_m(max_$(N)_$(state.id)_inputs_m[((max_$(N)_$(state.id)_i + 1)*$N - 1) -: $N]),
                    .y_p($(name(outputs[1]))_p[max_$(N)_$(state.id)_i]),
                    .y_m($(name(outputs[1]))_m[max_$(N)_$(state.id)_i])
                );""")
        write(buffer, """
            // END max_$(N)_$(state.id)
            \n""")
    else
        push!(netlist, Net(name = "max_$(N)_$(state.id)_inputs", size = (N, 1), signed = true))
        write(buffer, """
            $stdcomment
            // BEGIN max_$(N)_$(state.id)
            """)
        write(buffer, "assign max_$(N)_$(state.id)_inputs_p = {")
        write(buffer, join(map(fname -> fname("_p"), names), ", "))
        write(buffer, "};\n")
        write(buffer, "assign max_$(N)_$(state.id)_inputs_m = {")
        write(buffer, join(map(fname -> fname("_m"), names), ", "))
        write(buffer, "};\n")
        write(buffer, """
            stoch_signed_nmax #(
                    .NUM_INPUTS($N)
                ) max_$(N)_$(state.id) (
                    .CLK(CLK),
                    .nRST(nRST),
                    .as_p(max_$(N)_$(state.id)_inputs_p),
                    .as_m(max_$(N)_$(state.id)_inputs_m),
                    .y_p($(name(outputs[1]))_p),
                    .y_m($(name(outputs[1]))_m)
                );
            // END max_$(N)_$(state.id)
            \n""")
    end

    return buffer, (id = state.id + 1,)
end


text:
import Escher: Elem, Tile, render
import Base: convert

type VegaPlot <: Tile
    json::AbstractString
end

convert(::Type{Tile}, v::VegaVisualization) = VegaPlot(JSON.json(tojs(v)))

render(plot::VegaPlot, state) =
    Elem(:"vega-plot", attributes = Dict(:json=>plot.json))

max_stars_repo_name: jpthor/vega


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 6e487246-5d73-11eb-2f43-1f159557c68c
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
	using LinearAlgebra:dot
end

# ╔═╡ 41796f22-5d73-11eb-0365-f35b80ba9bbb
md"""
# Chapter 13 Short Corridor
"""

# ╔═╡ 8fa18536-5d73-11eb-398d-b93cbd63aedc
begin
	Base.@kwdef mutable struct ShortCorridorEnv <: AbstractEnv
		position::Int = 1
	end

	RLBase.state_space(env::ShortCorridorEnv) = Base.OneTo(4)
	RLBase.action_space(env::ShortCorridorEnv) = Base.OneTo(2)

	function (env::ShortCorridorEnv)(a)
		if env.position == 1 && a == 2
			env.position += 1
		elseif env.position == 2
			env.position += a == 1 ? 1 : -1
		elseif env.position == 3
			env.position += a == 1 ? -1 : 1
		end
		nothing
	end

	function RLBase.reset!(env::ShortCorridorEnv)
		env.position = 1
		nothing
	end

	RLBase.state(env::ShortCorridorEnv) = env.position
	RLBase.is_terminated(env::ShortCorridorEnv) = env.position == 4
	RLBase.reward(env::ShortCorridorEnv) =  env.position == 4 ? 0.0 : -1.0
end

# ╔═╡ 9a269140-5d73-11eb-01a6-a78c2b583887
world = ShortCorridorEnv()

# ╔═╡ aab4885a-5d73-11eb-33c7-b97e45de4016
ns, na = length(state_space(world)), length(action_space(world))

# ╔═╡ ae7c57c4-5d73-11eb-1734-f531d5b4cc6e
function run_once(A)
    avg_rewards = []
    for ϵ in A
        p = TabularRandomPolicy(;table=Dict(s => [1-ϵ, ϵ] for s in 1:ns))
        env = ShortCorridorEnv()
        hook=TotalRewardPerEpisode()
        run(p, env, StopAfterEpisode(1000),hook)
        push!(avg_rewards, mean(hook.rewards[end-100:end]))
    end
    avg_rewards
end

# ╔═╡ 2159d8de-5d74-11eb-0f9b-9d4fbb340ad6
X = 0.05:0.05:0.95

# ╔═╡ c2a474fe-5d73-11eb-0073-2f4afc5a1ada
plot(X, mean([run_once(X) for _ in 1:10]), legend=nothing)

# ╔═╡ d39bd212-5d73-11eb-0baf-b54fd44c2337
md"""
## REINFORCE Policy

Based on descriptions in Chapter 13.1, we need to define a new customized approximator.
"""

# ╔═╡ 33c22aaa-5d74-11eb-2271-e1c21336c28a
begin
	Base.@kwdef struct LinearPreferenceApproximator{F,O} <: AbstractApproximator
		weight::Vector{Float64}
		feature_func::F
		actions::Int
		opt::O
	end

	function (A::LinearPreferenceApproximator)(s)
		h = [dot(A.weight, A.feature_func(s, a)) for a in 1:A.actions]
		softmax(h)
	end

	function RLBase.update!(A::LinearPreferenceApproximator, correction::Pair)
		(s, a), Δ = correction
		w, x = A.weight, A.feature_func
		w̄ = -Δ .* (x(s,a) .- sum(A(s) .* [x(s, b) for b in 1:A.actions]))
		Flux.Optimise.update!(A.opt, w, w̄)
	end
end

# ╔═╡ 4d00b372-5d74-11eb-0760-11feb1439567
begin
	Base.@kwdef struct ReinforcePolicy{A<:AbstractApproximator} <: AbstractPolicy
		approximator::A
		γ::Float64
	end

	(p::ReinforcePolicy)(env::AbstractEnv) = prob(p, state(env)) |> WeightedExplorer(;is_normalized=true)

	RLBase.prob(p::ReinforcePolicy, s) = p.approximator(s)

	function RLBase.update!(
		p::ReinforcePolicy,
		t::AbstractTrajectory,
		::AbstractEnv,
		::PostEpisodeStage
	)
		S, A, R = t[:state], t[:action], t[:reward]
		Q, γ = p.approximator, p.γ
		G = 0.

		for i in 1:length(R)
			s,a,r = S[end-i], A[end-i], R[end-i+1]
			G = r + γ*G

			update!(Q, (s, a) => G)
		end
	end

	function RLBase.update!(
		t::AbstractTrajectory,
		::ReinforcePolicy,
		::AbstractEnv,
		::PreEpisodeStage
	)
		empty!(t)
	end
end

# ╔═╡ 5f0926f8-5d74-11eb-114c-77cfe8c57a11
function run_once_RL(α)
    agent = Agent(
        policy=ReinforcePolicy(
            approximator=LinearPreferenceApproximator(
                weight=[-1.47, 1.47],  # init_weight
                feature_func=(s,a) -> a == 1 ? [0, 1] : [1, 0],
                actions=na,
                opt=Descent(α)
            ),
            γ=1.0
        ),
        trajectory=VectorSARTTrajectory()
    )

    env = ShortCorridorEnv()
    hook = TotalRewardPerEpisode()
    run(agent,env,StopAfterEpisode(1000;is_show_progress=false),hook)
    hook.rewards
end

# ╔═╡ 6310c5c4-5d74-11eb-331a-8339e03cd161
begin
	fig_13_1 = plot(legend=:bottomright)
	for x in [-13, -14]  # for -12, it seems not that easy to converge in short time
		plot!(fig_13_1, mean(run_once_RL(2. ^ x) for _ in 1:50), label="alpha = 2^$x")
	end
	fig_13_1
end

# ╔═╡ 9bf6ff86-5d74-11eb-2235-1bada5c181d4
md"""
Interested in how to reproduce figure 13.2? A PR is warmly welcomed! See you there!
"""

# ╔═╡ Cell order:
# ╟─41796f22-5d73-11eb-0365-f35b80ba9bbb
# ╠═6e487246-5d73-11eb-2f43-1f159557c68c
# ╠═8fa18536-5d73-11eb-398d-b93cbd63aedc
# ╠═9a269140-5d73-11eb-01a6-a78c2b583887
# ╠═aab4885a-5d73-11eb-33c7-b97e45de4016
# ╠═ae7c57c4-5d73-11eb-1734-f531d5b4cc6e
# ╠═2159d8de-5d74-11eb-0f9b-9d4fbb340ad6
# ╠═c2a474fe-5d73-11eb-0073-2f4afc5a1ada
# ╟─d39bd212-5d73-11eb-0baf-b54fd44c2337
# ╠═33c22aaa-5d74-11eb-2271-e1c21336c28a
# ╠═4d00b372-5d74-11eb-0760-11feb1439567
# ╠═5f0926f8-5d74-11eb-114c-77cfe8c57a11
# ╠═6310c5c4-5d74-11eb-331a-8339e03cd161
# ╟─9bf6ff86-5d74-11eb-2235-1bada5c181d4


text:
# using Test, Distances, TransferEntropy

# d1 = rand(50)
# d2 = rand(50)
# r = rand(50)
# tol = 1e-8

# # Directly from time series
# @test standard_te(d1, r, estimator = :tetogrid) >= 0 - tol
# @test standard_te(d2, r, estimator = :tetogrid) >= 0 - tol
# @test standard_te(d1, r, estimator = :tefreq) >= 0 - tol
# @test standard_te(d2, r, estimator = :tefreq) >= 0 - tol
# @test standard_te(d1, r, estimator = :tekNN) >= 0 - tol
# @test standard_te(d2, r, estimator = :tekraskov) >= 0 - tol

# #########################
# # Surrogates
# #########################
# @test standard_te(d1, r, which_is_surr = :driver) >= 0 - tol
# @test standard_te(d1, r, which_is_surr = :response) >= 0 - tol
# @test standard_te(d1, r, which_is_surr = :none) >= 0 - tol

# @test standard_te(d1, r, which_is_surr = :both, surr_func = aaft) >= 0 - tol
# @test standard_te(d1, r, which_is_surr = :both, surr_func = iaaft) >= 0 - tol
# @test standard_te(d1, r, which_is_surr = :both, surr_func = randomphases) >= 0 - tol
# @test standard_te(d1, r, which_is_surr = :both, surr_func = randomamplitudes) >= 0 - tol
# @test standard_te(d1, r, which_is_surr = :both, surr_func = randomshuffle) >= 0 - tol

# #########################
# # Forward prediction lags
# #########################
# @test standard_te(d1, r, ν = 1) >= 0 - tol
# @test standard_te(d1, r, ν = 5) >= 0 - tol

# #########################
# # Embedding lags and dimensions
# #########################
# d1 = rand(100)
# d2 = rand(100)
# r = rand(100)
# @test standard_te(d1, r, ν = 1, τ = 2, dim = 4) >= 0 - tol
# @test standard_te(d1, r, ν = 3, τ = 3, dim = 5) >= 0 - tol

# #########################
# # Tuning the bins
# #########################
# @test standard_te(d1, r, n_ϵ = 10) >= 0 - tol
# @test standard_te(d1, r, n_ϵ = 10, max_numbins = 5, min_numbins = 2) >= 0 - tol

##########################
# Regular TE
##########################

using Test, TransferEntropy, CausalityToolsBase
k, l, m = 1, 1, 1
τ = 1
η = 1
n_subdivs = 3
te_res1 = te_reg(rand(1000), rand(1000), k, l, m, η = η, τ = τ, n_subdivs = n_subdivs)
te_res2 = te_reg(rand(1000), rand(1000), k + 1, l, m + 1, η = η, τ = τ + 1, n_subdivs = n_subdivs)

@test te_res1 isa Vector{<:Real}
@test te_res2 isa Vector{<:Real}
@test length(te_res1) == 4

@test te_reg(rand(100), rand(100), 1, 1, 1, RectangularBinning(4)) |> typeof <: Real
@test te_reg(rand(100), rand(100), 1, 1, 1, [RectangularBinning(4), RectangularBinning(0.2)]) isa Vector{<:Real}

##########################
# Conditional TE
##########################

# Default discretization schemes
using Test
k, l, m, n = 1, 1, 1, 1
τ = 1
η = 1
n_subdivs = 3
te_res1 = te_cond(rand(1000), rand(1000), rand(1000), k, l, m, n, η = η, τ = τ, n_subdivs = n_subdivs)
te_res2 = te_cond(rand(1000), rand(1000), rand(1000), k + 1, l, m + 1, n, η = η, τ = τ + 1, n_subdivs = n_subdivs)

@test te_res1 isa Vector{<:Real}
@test te_res2 isa Vector{<:Real}
@test length(te_res1) == 4

# User-provided binning schemes
using Test
x, y, z = rand(100), rand(100), rand(100)
@test te_cond(x, y, z, 1, 1, 1, 1, RectangularBinning(4)) |> typeof <: Real
@test te_cond(x, y, z, 1, 1, 1, 1, [RectangularBinning(4), RectangularBinning(0.2)]) isa Vector{<:Real}
@test te_cond(x, y, z, 1, 1, 1, 1, RectangularBinning(3), estimator = TransferOperatorGrid()) |> typeof <: Real


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite
# These files are JuMP implementations of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# This file will run all the HS test problems, allowing for easy testing
#############################################################################

if Pkg.installed("Ipopt") == nothing
    println("Cannot run NLP tests because Ipopt is not installed.")
    exit(0)
end

HS_path = dirname(@__FILE__)

files = readdir(HS_path)

println("Running HS Tests")
for file in files
    contains(file, "runhs") && continue
    try
        include(file)
        println("  $file  PASSED")
    catch
        println("  $file  FAILED")
    end
end


text:
#############################################################
#
# Saved Searches Functions
#
#############################################################

function get_saved_searches_show_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/saved_searches/show/$(id).json", options)

    return r.status == 200 ? Requests.json(r) : r

end

function post_saved_searches_create(query::AbstractString; options=Dict{AbstractString, AbstractString}())

    options["query"] = query

    r = post_oauth("https://api.twitter.com/1.1/saved_searches/create.json", options)

    return r.status == 200 ? Requests.json(r) : r

end

function post_saved_searches_destroy_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = post_oauth("https://api.twitter.com/1.1/saved_searches/destroy/$(id).json", options)

    return r.status == 200 ? Requests.json(r) : r

end


text:
# Pi reduction - vectorized
#
# Numeric integration to estimate pi
# Asks the user to select a device at runtime
# Vector size must be present as a CLI argument
#
# History: C version written by Tim Mattson, May 2010
#          Ported to the C++ Wrapper API by Benedict R. Gaster, September 2011
#          C++ version Updated by Tom Deakin and Simon McIntosh-Smith, October 2012
#          Ported to Python by Tom Deakin, July 2013
#          Ported to Julia by Jake Bolewski, Nov 2013

using OpenCL

# get the directory of this file
# (used for test runner)
src_dir = dirname(Base.source_path())

# Some constant values
INSTEPS = 512 * 512 * 512
ITERS = -1
WGS = -1
NAME = ""

if length(ARGS) < 1
    info("Usage: julia pi_vocl.jl [num] (where num = 1, 4, or 8)")
    exit(1)
end
vector_size = parse(Int, ARGS[1])

if vector_size == 1
        ITERS = 262144
        WGS = 8
elseif vector_size == 4
        ITERS = 65536 # (262144/4)
        WGS = 32
elseif vector_size == 8
        ITERS = 32768 # (262144/8)
        WGS = 64
else
    warn("Invalid vector size")
    exit(1)
end

# Set some default values:
# Default number of steps (updated later to device prefereable)
in_nsteps = INSTEPS

# Default number of iterations
niters = ITERS
work_group_size = WGS

# Create context, queue and build program

#---------------------------------------------
# Uncomment to switch between gpu/cpu devices
#---------------------------------------------
#device = first(cl.devices(:gpu))
#device = first(cl.devices(:cpu))
device, ctx, queue = cl.create_compute_context()

kernelsource = open(readall, joinpath(src_dir, "pi_vocl.cl"))
program = cl.Program(ctx, source=kernelsource) |> cl.build!

if vector_size == 1
    pi_kernel = cl.Kernel(program, "pi")
elseif vector_size == 4
    pi_kernel = cl.Kernel(program, "pi_vec4")
elseif vector_size == 8
    pi_kernel = cl.Kernel(program, "pi_vec8")
end

# Now that we know the size of the work_groups, we can set the number of work
# groups, the actual number of steps, and the step size
nwork_groups = in_nsteps ÷ (work_group_size * niters)

# get the max work group size for the kernel on our device
if vector_size == 1
    max_size = cl.work_group_info(pi_kernel, :size, device)
elseif vector_size == 4
    max_size = cl.work_group_info(pi_kernel, :size, device)
elseif vector_size == 8
    max_size = cl.work_group_info(pi_kernel, :size, device)
end

if max_size > work_group_size
    work_group_size = max_size
    nwork_groups = in_nsteps ÷ (work_group_size * niters)
end

if nwork_groups < 1
    nwork_groups = device[:max_compute_units]
    work_group_size = in_nsteps ÷ (nwork_groups * niters)
end

nsteps = work_group_size * niters * nwork_groups
step_size = 1.0 / nsteps

# vector to hold partial sum
h_psum = Array(Float32, nwork_groups)

println("$nwork_groups work groups of size $work_group_size.")
println("$nsteps integration steps")

d_partial_sums = cl.Buffer(Float32, ctx, :w, length(h_psum))

# start timer
rtime = time()

# Execute the kernel over the entire range of our 1d input data et
# using the maximum number of work group items for this device
# Set the global and local size as tuples
global_size = (nwork_groups * work_group_size,)
local_size  = (work_group_size,)
localmem    = cl.LocalMem(Float32, work_group_size)

cl.call(queue, pi_kernel, global_size, local_size,
        Int32(niters), Float32(step_size), localmem, d_partial_sums)

cl.copy!(queue, h_psum, d_partial_sums)

# complete the sum and compute final integral value
pi_res = sum(h_psum) * step_size

# stop the timer
rtime = time() - rtime

println("The calculation ran in $rtime secs")
println("pi=$pi_res for $nsteps steps")


text:
"""
    cl3(x::Float64)::Float64

Returns the value of the Clausen function ``\\operatorname{Cl}_3(x)``
for a real angle ``x`` of type `Float64`.  This function is defined as

```math
\\operatorname{Cl}_3(x) = \\Re[\\operatorname{Li}_3(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\cos(kx)}{k^3}
```

Author: Alexander Voigt

License: MIT

# Example
```jldoctest; setup = :(using ClausenFunctions), output = false
julia> cl3(1.0)
0.44857300728001737
```
"""
function cl3(x::Float64)::Float64
    zeta3 = 1.2020569031595943
    pi28 = pi*pi/8.0

    x = range_reduce_odd(x)

    if x == zero(x)
        return zeta3
    end

    if x < 0.5*pi
        P = (-7.5000000000000001e-01,  1.5707637881835541e-02,
             -3.5426736843494423e-05, -2.4408931585123682e-07)
        Q = (1.0000000000000000e+00, -2.5573146805410089e-02,
             1.5019774853075050e-04, -1.0648552418111624e-07)
        y = x*x
        y2 = y*y
        p = P[1] + y * P[2] + y2 * (P[3] + y * P[4])
        q = Q[1] + y * Q[2] + y2 * (Q[3] + y * Q[4])
        zeta3 + y*(p/q + 0.5*log(x))
    else
        P = (-4.9017024647634973e-01, 4.1559155224660940e-01,
             -7.9425531417806701e-02, 5.9420152260602943e-03,
             -1.8302227163540190e-04, 1.8027408929418533e-06)
        Q = (1.0000000000000000e+00, -1.9495887541644712e-01,
             1.2059410236484074e-02, -2.5235889467301620e-04,
             1.0199322763377861e-06,  1.9612106499469264e-09)
        y = pi - x
        z = y*y - pi28
        z2 = z*z
        z4 = z2*z2
        p = P[1] + z * P[2] + z2 * (P[3] + z * P[4]) + z4 * (P[5] + z * P[6])
        q = Q[1] + z * Q[2] + z2 * (Q[3] + z * Q[4]) + z4 * (Q[5] + z * Q[6])
        p/q
    end
end


text:
using Plots, Resample, Dates

times_long = DateTime(2019):Minute(5):DateTime(2019, 1, 1, 23, 55)
times_short = DateTime(2019):Hour(1):DateTime(2019, 1, 1, 23)

power_long = rand(length(times_long))
energy_long = power_long / 12

power_short = resample(power_long, times_long, times_short)
energy_short = resample(energy_long, times_long, times_short, Sum())

plot(times_long, power_long, label = "5 min", st = :step)
plot!(times_short, power_short, label = "1 h", st = :step)
png(joinpath(@__DIR__, "power"))

plot(times_long, cumsum(energy_long), label = "5 min", st = :step)
plot!(times_short, cumsum(energy_short), label = "1 h", st = :step)
png(joinpath(@__DIR__, "energy"))


text:
using OpenQuantumSystems
import QuantumOpticsBase, LinearAlgebra, OrdinaryDiffEq, QuadGK, DelayDiffEq, Interpolations

function interpolate_with_tspan(itp, tspan, t)
    t0 = tspan[1]
    t1 = tspan[end]
    N = length(tspan)
    if t < t0
        t_i = 1
    elseif t > t1
        t_i = N
    else
        t_i = (t-t0)/t1 * (N-1) + 1
    end
    return itp[t_i]
end

function normalize_bath(W_bath, aggCore, aggTools, aggOperators)
    elLen = aggCore.molCount+1
    indicesMap = aggTools.indicesMap
    W_bath_tr = trace_bath(W_bath, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
    for b=1:elLen
        for a=1:elLen
            a1 = indicesMap[a][1]
            a2 = indicesMap[a][end]
            b1 = indicesMap[b][1]
            b2 = indicesMap[b][end]
            W_bath[a1:a2, b1:b2] /= W_bath_tr[a, b]
        end
    end
    return W_bath
end

function W_aabb_1_bath_core(t, s, p, tmp1, tmp2; bath_evolution=:none, K_rtol=1e-12, K_atol=1e-12)
    aggCore, aggTools, aggOperators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, _ = p

    # W_1_bath = deepcopy(W0_bath.data)
    elLen = aggCore.molCount+1
    indicesMap = aggTools.indicesMap
    Ham_0 = aggOperators.Ham_0
    Ham = aggOperators.Ham

    W_1_bath = zeros(ComplexF64, aggTools.bSize, aggTools.bSize)
    K_aabb_s = K_aabb_W_bath_intp(s, p, tmp1, tmp2; rtol=K_rtol, atol=K_atol)
    rho_s = interpolate_with_tspan(rho_0_int_t_itp, tspan, s)
    W_bath_s = interpolate_with_tspan(W_0_bath_t_itp, tspan, s)
    W_1_bath[:, :] = W_bath_s[:, :]
    for b=2:elLen
        if bath_evolution == :interaction_picture
            U_bb_ = evolution_el_part(Ham.data, s, b, b, indicesMap)
            U_0_bb = evolution_el_part(Ham_0.data, s, b, b, indicesMap)
            U_bb = U_0_bb * U_bb_
        elseif bath_evolution == :shroedinger_picture
            U_bb = evolution_el_part(Ham.data, s, b, b, indicesMap)
        end
        for a=2:elLen
            a1 = indicesMap[a][1]
            a2 = indicesMap[a][end]
            b1 = indicesMap[b][1]
            b2 = indicesMap[b][end]
            if bath_evolution == :interaction_picture
                U_aa_ = evolution_el_part(Ham.data, t-s, a, a, indicesMap)
                U_0_aa = evolution_el_part(Ham_0.data, t-s, a, a, indicesMap)
                U_aa = U_0_aa * U_aa_
            elseif bath_evolution == :shroedinger_picture
                U_aa = evolution_el_part(Ham.data, t-s, a, a, indicesMap)
            end
            # 1 / rho_s[a, a] is already in K rate constant
            if bath_evolution == :none
                W_1_bath[a1:a2, a1:a2] = W_1_bath[a1:a2, a1:a2] +
                    K_aabb_s[a, b] * rho_s[b, b] * W_bath_s[b1:b2, b1:b2]
            else
                U_ = U_aa * U_bb
                W_1_bath[a1:a2, a1:a2] = W_1_bath[a1:a2, a1:a2] +
                    K_aabb_s[a, b] * rho_s[b, b]  * U_ * W_bath_s[b1:b2, b1:b2] * adjoint(U_)
            end
        end
    end
    return W_1_bath
end

function W_abcd_1_bath_core(t, s, p, tmp1, tmp2; bath_evolution=:none, K_rtol=1e-12, K_atol=1e-12)
    aggCore, aggTools, aggOperators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, _ = p

    # W_1_bath = deepcopy(W0_bath.data)
    elLen = aggCore.molCount+1
    indicesMap = aggTools.indicesMap
    Ham_0 = aggOperators.Ham_0
    Ham = aggOperators.Ham

    W_1_bath = zeros(ComplexF64, aggTools.bSize, aggTools.bSize)
    K_abcd_s = K_abcd_W_bath_intp(s, p, tmp1, tmp2; rtol=K_rtol, atol=K_atol)
    rho_s = interpolate_with_tspan(rho_0_int_t_itp, tspan, s)
    W_bath_s = interpolate_with_tspan(W_0_bath_t_itp, tspan, s)
    W_1_bath[:, :] = W_bath_s[:, :]
    for c=2:elLen, d=2:elLen
        c1 = indicesMap[c][1]
        c2 = indicesMap[c][end]
        d1 = indicesMap[d][1]
        d2 = indicesMap[d][end]
        if bath_evolution == :interaction_picture
            U_cd_ = evolution_el_part(Ham.data, s, c, d, indicesMap)
            U_0_cd = evolution_el_part(Ham_0.data, s, c, d, indicesMap)
            U_cd = U_0_cd * U_cd_
        elseif bath_evolution == :shroedinger_picture
            U_cd = evolution_el_part(Ham.data, s, c, d, indicesMap)
        end
        for a=2:elLen, b=1:elLen
            a1 = indicesMap[a][1]
            a2 = indicesMap[a][end]
            b1 = indicesMap[b][1]
            b2 = indicesMap[b][end]
            if bath_evolution == :interaction_picture
                U_ab_ = evolution_el_part(Ham.data, t-s, a, b, indicesMap)
                U_0_ab = evolution_el_part(Ham_0.data, t-s, a, b, indicesMap)
                U_ab = U_0_ab * U_ab_
            elseif bath_evolution == :shroedinger_picture
                U_aa = evolution_el_part(Ham.data, t-s, a, b, indicesMap)
            end
            # 1 / rho_s[a, b] is already in K rate constant
            if bath_evolution == :none
                W_1_bath[a1:a2, b1:b2] = W_1_bath[a1:a2, b1:b2] +
                    K_abcd_s[a, b, c, d] * rho_s[c, d] * W_bath_s[c1:c2, d1:d2]
            else
                U_ = U_aa * U_bb
                W_1_bath[a1:a2, a1:a2] = W_1_bath[a1:a2, a1:a2] +
                K_abcd_s[a, b, c, d] * rho_s[c, d] * U_ * W_bath_s[c1:c2, d1:d2]  * adjoint(U_)
            end
        end
    end
    return W_1_bath
end

function W_1_bath(
        t, p, tmp1, tmp2;
        bath_evolution=:none, bath_ansatz=:population, normalize=false,
        W_1_rtol=1e-12, W_1_atol=1e-12, K_rtol=1e-12, K_atol=1e-12
    )
    aggCore, aggTools, aggOperators, W0, W0_bath, _ = p

    if bath_ansatz == :population
        W_1_diff, err = QuadGK.quadgk(
            s -> W_aabb_1_bath_core(t, s, p, tmp1, tmp2; bath_evolution=bath_evolution, K_rtol=K_rtol, K_atol=K_atol),
            0,
            t,
            rtol = W_1_rtol,
            atol = W_1_atol,
        )
    elseif bath_ansatz == :population_coherences
        W_1_diff, err = QuadGK.quadgk(
            s -> W_abcd_1_bath_core(t, s, p, tmp1, tmp2; bath_evolution=bath_evolution, K_rtol=K_rtol, K_atol=K_atol),
            0,
            t,
            rtol = W_1_rtol,
            atol = W_1_atol,
        )
    end

    W_1_bath = deepcopy(W0_bath.data) + W_1_diff
    if normalize
        W_1_bath = normalize_bath(W_1_bath, aggCore, aggTools, aggOperators)
    end
    return W_1_bath
end

function W_1_bath(t, W0, W0_bath, agg::Aggregate; W_1_rtol=1e-12, W_1_atol=1e-12, K_rtol=1e-12, K_atol=1e-12)
    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))
    return W_1_bath(t, p, tmp1, tmp2; W_1_rtol=W_1_rtol, W_1_atol=W_1_atol, K_rtol=K_rtol, K_atol=K_atol)
end

function QME_sI_iterative(
    W0::T,
    rho_0_int_t,
    W_0_bath_t,
    tspan::Array,
    agg::Aggregate;
    bath_evolution=:none,
    bath_ansatz=:population,
    normalize=false,
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    W_1_rtol::AbstractFloat = 1e-12,
    W_1_atol::AbstractFloat = 1e-12,
    K_rtol::AbstractFloat = 1e-12,
    K_atol::AbstractFloat = 1e-12,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)

    # Calculate and interpolate rho_0_int_t, W_0_bath_t, W_1_bath_t
    if ndims(rho_0_int_t) == 1 && typeof(rho_0_int_t[1]) <: Operator
        rho_0_int_t = operator_recast(rho_0_int_t)
    end
    if ndims(rho_0_int_t) == 3
        rho_0_int_t = [rho_0_int_t[i, :, :] for i=1:size(rho_0_int_t, 1)]
    end
    rho_0_int_t_itp = Interpolations.interpolate(
        rho_0_int_t,
        Interpolations.BSpline(Interpolations.Linear())
    )

    if ndims(rho_0_int_t) == 1 && typeof(rho_0_int_t[1]) <: Operator
        W_0_bath_t = operator_recast(W_0_bath_t)
    end
    if ndims(W_0_bath_t) == 3
        W_0_bath_t = [W_0_bath_t[i, :, :] for i=1:size(W_0_bath_t, 1)]
    end
    W_0_bath_t_itp = Interpolations.interpolate(
        W_0_bath_t,
        Interpolations.BSpline(Interpolations.Linear())
    )

    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, eltype(W0))
    elLen = agg.core.molCount+1
    W_1_bath_t = []
    for t_i=1:length(tspan)
        t = tspan[t_i]
        W_1_bath_ = W_1_bath(
            t, p, tmp1, tmp2;
            bath_evolution=bath_evolution, bath_ansatz=bath_ansatz, normalize=normalize,
            W_1_rtol=W_1_rtol, W_1_atol=W_1_atol, K_rtol=K_rtol, K_atol=K_atol
        )
        push!(W_1_bath_t, W_1_bath_)
    end
    W_1_bath_itp = Interpolations.interpolate(W_1_bath_t, Interpolations.BSpline(Interpolations.Linear()))

    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, W_1_bath_itp, tspan, eltype(W0))

    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_iterative(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    tspan, rho_int_1_t = OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
    return tspan, rho_int_1_t, W_1_bath_t
end

function dQME_sI_iterative(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, _, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_iterative(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_iterative(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, W0, W0_bath, W_1_bath_itp, tspan, _ = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    W0_int_s = interpolate_with_tspan(W_1_bath_itp, tspan, s)
    tmp1[:, :] = ad(rho, W0_int_s, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
mutable struct Preamble
  ins::Int
  preceding::Array{Int,1}
  Preamble(i,p) = 1 <= i <= length(p) ? new(i,p) : error("Bad insert index")
end

function Base.push!(preamble::Preamble, value::Int)
  preamble.preceding[preamble.ins] = value
  preamble.ins = preamble.ins % length(preamble.preceding) + 1
  preamble
end

function follows(next::Int, preamble::Preamble)
  for i in 1:length(preamble.preceding), j in i+1:length(preamble.preceding)
      if preamble.preceding[i] + preamble.preceding[j] == next
        return true
      end
  end
  false
end

input = parse.(Int, readlines("a9.txt"))

function findnonfollowing()
  buffer = Preamble(1, input[1:25])
  for value in input[26:end]
    if !follows(value, buffer)
      return value
    end
    push!(buffer, value)
  end
end

part1 = findnonfollowing()
println("$part1")

function findcontinuoussum(target)
  first = 1
  last = 2
  cs = input[first] + input[last]
  while cs != target && (first < last || last < length(input))
    while cs < target
      last += 1
      cs += input[last]
    end
    while cs > target
      cs -= input[first]
      first += 1
    end
  end
  cs == target && sum(extrema(input[first:last]))
end

part2 = findcontinuoussum(part1)
println("$part2")


text:
using GitHub
using JLD
using Base.Test

event_request = JLD.load(joinpath(dirname(@__FILE__), "commit_comment.jld"), "request")

event_json = Requests.json(event_request)

event = GitHub.event_from_payload!("commit_comment", event_json)

################
# WebhookEvent #
################

@test get(event.repository.name) == "BenchmarkTrackers.jl"
@test get(event.sender.login) == "jrevels"

#################
# EventListener #
#################

@test !(GitHub.has_valid_secret(event_request, "wrong"))
@test GitHub.has_valid_secret(event_request, "secret")
@test !(GitHub.is_valid_event(event_request, ["wrong"]))
@test GitHub.is_valid_event(event_request, ["commit_comment"])
@test !(GitHub.from_valid_repo(event, ["JuliaWeb/GitHub.jl"]))
@test GitHub.from_valid_repo(event, ["JuliaCI/BenchmarkTrackers.jl"])
@test GitHub.handle_event_request(event_request, x -> true,
                                  secret = "secret",
                                  events = ["commit_comment"],
                                  repos = ["JuliaCI/BenchmarkTrackers.jl"])

@test begin
    listener = EventListener(x -> true;
                             secret = "secret",
                             repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                             events = ["commit_comment"],
                             forwards = ["http://bob.com", HttpCommon.URI("http://jim.org")])
    r = listener.server.http.handle(HttpCommon.Request(),HttpCommon.Response())
    r.status == 400
end

###################
# CommentListener #
###################

result = GitHub.handle_comment((e, m) -> m, event, GitHub.AnonymousAuth(), r"`RunBenchmarks\(.*?\)`", false)

@test result.match == "`RunBenchmarks(\"binary\", \"unary\")`"

@test begin
    listener = CommentListener((x, y) -> true, r"trigger";
                               secret = "secret",
                               repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                               forwards = ["http://bob.com", HttpCommon.URI("http://jim.org")],
                               check_collab = false)
    r = listener.listener.server.http.handle(HttpCommon.Request(), HttpCommon.Response())
    r.status == 400
end


text:
include("simple.jl")
include("contrast.jl")
include("scalesigned.jl")
include("tile.jl")
include("annotations.jl")
include("statusbar.jl")
include("test4d.jl")
ImageView.imshow(img)

include("newtests.jl")

ImageView.closeall()


text:
lines = readlines("input.txt")

edges = [Tuple(split(line, '-')) for line = lines]
nodes = unique(vcat(collect.(edges)...))
N = length(nodes)
adjacency = falses((N,N))
node2idx = n -> findfirst(nodes .== n)
idx2node = i -> nodes[i]

function isupper(s)
    return isuppercase.(collect(s)) |> all
end

function DFS(from, to, cond, state)
    state = [state..., from]
    if from == to
        names = idx2node.(state)
        return [state]
    else
        C_n = findall(adjacency[from, :] .== true)
        paths = [DFS(n, to, cond, state) for n ∈ C_n if cond(from, n, state)]
        return length(paths) > 0 ? reduce(vcat, paths) : []
    end
end

for (from, to) = edges
    ix, iy = node2idx(from), node2idx(to)
    adjacency[ix,iy] = adjacency[iy,ix] = true
end

# big oof.
function twice(c)
    return (from, to, state) -> from != node2idx("end") && (
        isupper(idx2node(to))
        || !(to ∈ state)
        || (to == c && count(==(c), state) < 2)
    )
end

paths = []
for c = nodes
    if c in ["start", "end"]
        continue
    end
    append!(paths, DFS(node2idx("start"), node2idx("end"), twice(node2idx(c)), []))
end
println(length(Set(vcat(paths))))


text:
"""
TestProcessingCore.jl defines the TestProcessingCore types and methods

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

using Boid: Boid
using Boid: AbstractProcessingCore

# --- Type definitions

mutable struct TestProcessingCore <: AbstractProcessingCore
    # Data
    field_1::Float64
    field_2::Float64

    # Default constructor
    TestProcessingCore() = new(0.1, 1)
end

# --- Method definitions

function Boid.process_data!(processing_core::TestProcessingCore, data::Vector)
    processing_core.field_1 = data[1]
    processing_core.field_2 = data[2]
    return processing_core.field_1 + processing_core.field_2
end


text:
module TestDataFrame

using GarishPrint
using DataFrames
using Test

struct ABC{T1, T2, T3}
    hee::T1
    haa::T2
    hoo::T3
end

struct Example{T1, T2}
    field_a::T1
    field_b::T2
    abc::ABC
end

df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"])

x = Example(
    Dict(
        "a"=>Example(
            [1, 2, 3],
            2.0,
            ABC(1, 2.0im, 3.12f0),
        ),
        "str" => Set([1, 2, 3]),
        "boolean"=> false,
        "missing" => missing,
        "empty set" => Set(),
        "set" => Set([1, 2, 3]),
        "set{any}" => Set(Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),
        "type" => Any,
        "matrix" => rand(10, 10),
        "dataframe" => df,
        "named_tuple" => (name="ABC", age=25),
        "nested" => Example(
            Dict(
                "a"=>Example(
                    [1, 2, 3],
                    2.0,
                    ABC(1, 2.0im, 3.12f0),
                ),
            ),
            undef,
            ABC(nothing, 1.2+2.1im, π),
        )
    ),
    undef,
    ABC(nothing, 1.2+2.1im, π),
)

pprint(x)
pprint(x; color=false)
pprint(x; compact=true)
pprint(x; show_indent=false)

pprint(IOContext(stdout, :color=>false), x)

end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo
using Random
using LinearAlgebra
using JLD2
include(joinpath(@__DIR__, "algorithms/ars.jl")) # augmented random search

# ## Ant
env = get_environment(:ant,
    representation=:minimal,
    gravity=-9.81,
    timestep=0.05,
    damper=50.0,
    spring=25.0,
    friction_coefficient=0.5,
    contact_feet=true,
    contact_body=true);

obs = reset(env)
initialize!(env.mechanism, :ant,
    body_position=[0.0, 0.0, 1.0],
    body_orientation=[0.0, 0.0, 0.0])
env.state .= get_minimal_state(env.mechanism)
render(env)

# ## Open visualizer
initialize!(env.mechanism, :ant)
open(env.vis)

# ## Set up policy
hp = HyperParameters(
        main_loop_size=100,
        horizon=150,
        n_directions=6,
        b=6,
        step_size=0.02)
input_size = length(obs)
output_size = length(env.input_previous)
normalizer = Normalizer(input_size)

# ## Training
train_times = Float64[]
rewards = Float64[]
policies = Matrix{Float64}[]
N = 5
for i = 1:N
    ## Reset environment
    env = get_environment(:ant,
        representation=:minimal,
        gravity=-9.81,
        timestep=0.05,
        damper=50.0,
        spring=25.0,
        friction_coefficient=0.5,
        contact_feet=true,
        contact_body=true)
    obs = reset(env)

    ## Random policy
    Random.seed!(i)
    hp = HyperParameters(
        main_loop_size=100,
        horizon=150,
        n_directions=6,
        b=6,
        step_size=0.02)
    input_size = length(obs)
    output_size = length(env.input_previous)
    normalizer = Normalizer(input_size)
    policy = Policy(input_size, output_size, hp)

    ## Train policy
    train_time = @elapsed train(env, policy, normalizer, hp)

    ## Evaluate policy
    reward = rollout_policy(policy.θ, env, normalizer, hp)

    ## Cache
    push!(train_times, train_time)
    push!(rewards, reward)
    push!(policies, policy.θ)
end

## @save joinpath(@__DIR__, "results/ant_rl.jld2") train_times rewards policies
@load joinpath(@__DIR__, "results/ant_rl.jld2") train_times rewards policies

# ## Training statistics
N_best = 3
max_idx = sortperm(rewards,
    lt=Base.isgreater)
train_time_best = (train_times[max_idx])[1:N_best]
rewards_best = (rewards[max_idx])[1:N_best]
policies_best = (policies[max_idx])[1:N_best]

@show rewards
@show mean(train_time_best)
@show std(train_time_best)
@show mean(rewards)
@show std(rewards)

# ## Save/Load policy
## θ = policy.θ
## @save joinpath(@__DIR__, "ant_policy.jld2") θ
## @load joinpath(@__DIR__, "ant_policy.jld2") θ

# ## Recover policy
hp = HyperParameters(
        main_loop_size=30,
        horizon=80,
        n_directions=6,
        b=6,
        step_size=0.02)
input_size = length(obs)
output_size = length(env.input_previous)
normalizer = Normalizer(input_size)
θ = policies_best[1]

# ## Visualize policy
## traj = display_random_policy(env, hp)
traj = display_policy(env, Policy(hp, θ), normalizer, hp)
visualize(env, traj)
open(env.vis)


text:
#' ---
#' title: Rocket Control
#' ---

#' **Originally Contributed by**: Iain Dunning

#' This tutorial shows how to solve a nonlinear rocketry control problem.
#' The problem was drawn from the [COPS3](http://www.mcs.anl.gov/~more/cops/cops3.pdf) benchmark.

#' Our goal is to maximize the final altitude of a vertically launched rocket.
#' We can control the thrust of the rocket, and must take account of
#' the rocket mass, fuel consumption rate, gravity, and aerodynamic drag.

#' Let us consider the basic description of the model (for the full description,
#' including parameters for the rocket, see the COPS3 PDF)

#' ### Overview
#' We will use a discretized model of time, with a fixed number of time steps, $n$.
#' We will make the time step size $\Delta t$, and thus the final time $t_f = n \cdot \Delta t$, a variable in the problem.
#' To approximate the derivatives in the problem we will use the [trapezoidal rule](http://en.wikipedia.org/wiki/Trapezoidal_rule).

#' ### State and Control
#' We will have three state variables:
#'
#' * Velocity, $v$
#' * Altitude, $h$
#' * Mass of rocket and remaining fuel, $m$
#'
#' and a single control variable, thrust $T$.
#' Our goal is thus to maximize $h(t_f)$.
#' Each of these corresponds to a JuMP variable indexed by the time step.

#' ### Dynamics
#' We have three equations that control the dynamics of the rocket:
#'
#' Rate of ascent: $$h^\prime = v$$
#' Acceleration: $$v^\prime = \frac{T - D(h,v)}{m} - g(h)$$
#' Rate of mass loss: $$m^\prime = -\frac{T}{c}$$
#'
#' where drag $D(h,v)$ is a function of altitude and velocity, and gravity $g(h)$ is a function of altitude.

#' These forces are defined as
#'
#' $$D(h,v) = D_c v^2 exp\left( -h_c \left( \frac{h-h(0)}{h(0)} \right) \right)$$
#' and
#' $$g(h) = g_0 \left( \frac{h(0)}{h} \right)^2$$
#'
#' The three rate equations correspond to JuMP constraints,
#' and for convenience we will represent the forces with nonlinear expressions.

using JuMP, Ipopt

# Create JuMP model, using Ipopt as the solver
rocket = Model(with_optimizer(Ipopt.Optimizer, print_level=0))

# Constants
# Note that all parameters in the model have been normalized
# to be dimensionless. See the COPS3 paper for more info.
h_0 = 1    # Initial height
v_0 = 0    # Initial velocity
m_0 = 1    # Initial mass
g_0 = 1    # Gravity at the surface

T_c = 3.5  # Used for thrust
h_c = 500  # Used for drag
v_c = 620  # Used for drag
m_c = 0.6  # Fraction of initial mass left at end

c     = 0.5 * sqrt(g_0 * h_0)  # Thrust-to-fuel mass
m_f   = m_c * m_0            # Final mass
D_c   = 0.5 * v_c * m_0 / g_0    # Drag scaling
T_max = T_c * g_0 * m_0        # Maximum thrust

n = 800   # Time steps

@variables(rocket, begin
    Δt ≥ 0, (start = 1/n) # Time step
    # State variables
    v[1:n] ≥ 0            # Velocity
    h[1:n] ≥ h_0          # Height
    m_f ≤ m[1:n] ≤ m_0    # Mass
    # Control
    0 ≤ T[1:n] ≤ T_max    # Thrust
end)

# Objective: maximize altitude at end of time of flight
@objective(rocket, Max, h[n])

# Initial conditions
@constraints(rocket, begin
    v[1] == v_0
    h[1] == h_0
    m[1] == m_0
    m[n] == m_f
end)

# Forces
# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )
@NLexpression(rocket, drag[j = 1:n], D_c * (v[j]^2) * exp(-h_c * (h[j] - h_0) / h_0))
# Grav(h)   = go * (h0 / h)^2
@NLexpression(rocket, grav[j = 1:n], g_0 * (h_0 / h[j])^2)
# Time of flight
@NLexpression(rocket, t_f, Δt * n)

# Dynamics
for j in 2:n
    # h' = v

    # Rectangular integration
    # @NLconstraint(rocket, h[j] == h[j - 1] + Δt * v[j - 1])

    # Trapezoidal integration
    @NLconstraint(rocket,
        h[j] == h[j - 1] + 0.5 * Δt * (v[j] + v[j - 1]))

    # v' = (T-D(h,v))/m - g(h)

    # Rectangular integration
    # @NLconstraint(rocket, v[j] == v[j - 1] + Δt *(
    #                 (T[j - 1] - drag[j - 1]) / m[j - 1] - grav[j - 1]))

    # Trapezoidal integration
    @NLconstraint(rocket,
        v[j] == v[j-1] + 0.5 * Δt * (
            (T[j] - drag[j] - m[j] * grav[j]) / m[j] +
            (T[j - 1] - drag[j - 1] - m[j - 1] * grav[j - 1]) / m[j - 1]))

    # m' = -T/c

    # Rectangular integration
    # @NLconstraint(rocket, m[j] == m[j - 1] - Δt * T[j - 1] / c)

    # Trapezoidal integration
    @NLconstraint(rocket,
        m[j] == m[j - 1] - 0.5 * Δt * (T[j] + T[j-1]) / c)
end

# Solve for the control and state
println("Solving...")
status = optimize!(rocket)

# Display results
# println("Solver status: ", status)
println("Max height: ", objective_value(rocket))

#+

# Can visualize the state and control variables
using Gadfly

#+

h_plot = plot(x = (1:n) * value.(Δt), y = value.(h)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Altitude"))
m_plot = plot(x = (1:n) * value.(Δt), y = value.(m)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Mass"))
v_plot = plot(x = (1:n) * value.(Δt), y = value.(v)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Velocity"))
T_plot = plot(x = (1:n) * value.(Δt), y = value.(T)[:], Geom.line,
                Guide.xlabel("Time (s)"), Guide.ylabel("Thrust"))
draw(SVG(6inch, 6inch), vstack(hstack(h_plot, m_plot), hstack(v_plot, T_plot)))


text:
"""
    from_dict(::Type{T}, d::AbstractDict{String}; kw...) where T

Convert dictionary `d` to an option type `T`, the value of valid fields of `T`
in this dictionary `d` can be override by keyword arguments.

# Example

```julia-repl
julia> @option struct OptionA
           name::String = "Sam"
           age::Int = 25
       end

julia> d = Dict{String, Any}(
           "name" => "Roger",
           "age" => 10,
       );

julia> from_dict(OptionA, d; age=25)
OptionA(;
    name = "Roger",
    age = 25,
)
```
"""
function from_dict(::Type{OptionType}, d::AbstractDict{String}; kw...) where {OptionType}
    if !isempty(kw)
        d = from_underscore_kwargs!(deepcopy(d), OptionType; kw...)
    end
    return from_dict_dynamic(OptionType, d)
end

struct OptionField{name} end
OptionField(name::Symbol) = OptionField{name}()

function from_dict(::Type{OptionType}, ::OptionField, ::Type{T}, x) where {OptionType, T}
    return from_dict(OptionType, T, x)
end

function from_dict(::Type{OptionType}, ::OptionField{field}, ::Type{T}, x) where {OptionType, field, T <: AbstractVector}
    if eltype(T) isa Union
        return map(x) do each
            from_dict_union_type(OptionType, field, eltype(T), each)
        end
    elseif is_option(eltype(T))
        return map(x) do each
            from_dict_option_type(OptionType, field, eltype(T), each)
        end
    else
        return map(x) do each
            from_dict_other_type(OptionType, field, eltype(T), each)
        end
    end
end

function from_dict(::Type{OptionType}, ::Type{T}, x) where {OptionType, T}
    is_option(T) && return from_dict(T, x)
    # TODO: deprecate convert_to_option
    # then just use the following
    # return convert(T, x)
    return deprecated_conversion(OptionType, T, x)
    # return convert(T, x)
end

from_dict(::Type, ::Type{VersionNumber}, x) = VersionNumber(x)

function deprecated_conversion(::Type{OptionType}, ::Type{T}, x) where {OptionType, T}
    ret = convert_to_option(OptionType, T, x)
    if ret === ConvertNotFound()
        return convert(T, x)
    else
        return ret
    end
end

function from_dict_dynamic(::Type{OptionType}, d::AbstractDict{String}) where {OptionType}
    isconcretetype(OptionType) || throw(ArgumentError("expect concrete type, got $OptionType"))

    nf = fieldcount(OptionType)
    args = ntuple(nf) do f_idx
        f_name = fieldname(OptionType, f_idx)
        f_type = fieldtype(OptionType, f_idx)
        f_default = field_default(OptionType, f_name)
        key = string(f_name)
        if !haskey(d, key)
            f_default === no_default && error("expect key: $key")
            return f_default
        end

        value = d[key]
        # empty dict is treated as nothing
        if value isa AbstractDict && isempty(value) && f_default === nothing
            return nothing
        end

        return if is_option(f_type)
            from_dict_option_type(OptionType, f_name, f_type, value)
        elseif f_type isa Union
            from_dict_union_type(OptionType, f_name, f_type, value)
        elseif f_type <: Reflect
            type_alias(OptionType) == value || parse_jltype(value) <: OptionType ||
                throw(ArgumentError("type mismatch, expect $OptionType got $value"))
            Reflect()
        else
            from_dict_other_type(OptionType, f_name, f_type, value)
        end
    end
    return OptionType(args...)
end

function from_dict_option_type(::Type{OptionType}, f_name::Symbol, ::Type{T}, value) where {OptionType, T}
    value isa AbstractDict || error("expect an AbstractDict, got $(typeof(value))")
    return from_dict_dynamic(T, value)
end

function from_dict_other_type(::Type{OptionType}, f_name::Symbol, ::Type{T}, value) where {OptionType, T}
    return from_dict(OptionType, OptionField(f_name), T, value)
end

function from_dict_union_type(::Type{OptionType}, f_name::Symbol, ::Type{FieldType}, value) where {OptionType, FieldType}
    assert_duplicated_alias_union(FieldType)

    for T in Base.uniontypes(FieldType)
        if is_option(T)
            alias = type_alias(T)
            reflect_idx = find_relfect_field(T)
            # type can be determined by alias
            if haskey(value, alias)
                return from_dict_dynamic(T, value[alias])
            elseif reflect_idx !== nothing # type can be determined by reflect field
                reflect_f_name = fieldname(T, reflect_idx)
                reflect_key = string(reflect_f_name)
                if haskey(value, reflect_key)
                    if alias == value[reflect_key] # find alias in reflect
                        return from_dict_dynamic(T, value)
                    else
                        type = tryparse_jltype(value[reflect_key])
                        type === nothing && continue
                        # NOTE: type is always more specialized
                        type <: T && return from_dict_dynamic(type, value)
                    end
                end
                continue
            end

            if alias === nothing && reflect_idx === nothing
                # not much information, try parse
                try
                    return from_dict_dynamic(T, value)
                catch
                    continue
                end
            end
        else
            try
                return from_dict(OptionType, OptionField(f_name), T, value)
            catch
                continue
            end
        end
    end
    error("cannot parse field $f_name, expect $FieldType, got $(typeof(value))")
end

function find_relfect_field(::Type{OptionType}) where {OptionType}
    nf = fieldcount(OptionType)
    for f_idx in 1:nf
        f_type = fieldtype(OptionType, f_idx)
        if f_type <: Reflect
            return f_idx
        end
    end
    return
end

function assert_duplicated_alias_union(::Type{UnionType}) where {UnionType}
    set = String[]
    for T in Base.uniontypes(UnionType)
        if is_option(T)
            alias = type_alias(T)
            if alias !== nothing
                alias in set && throw(DuplicatedAliasError(alias))
                push!(set, alias)
            end
        end
    end
    return
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

###############
# Prepare FMU #
###############

cd(dirname(@__FILE__))
pathToFMU = joinpath(pwd(), "../model/IO.fmu")

myFMU = fmiLoad(pathToFMU)

#########################
# Testing Setter/Getter #
#########################

# test reference conversion
@test fmi2String2ValueReference(myFMU, "p_real") == 16777216
@test fmi2String2ValueReference(myFMU, "p_integer") == 16777217
@test fmi2ValueReference2String(myFMU, Cint(16777217)) == ["p_integer"]

# create arrays for in place getter
vR = zeros(Real, 3)
vI = zeros(Integer, 3)
vB = zeros(Bool, 3)
vS = ["test","test"]

c1 = fmiInstantiate!(myFMU; loggingOn=true)

fmiEnterInitializationMode(myFMU)
@test fmiSetReal(myFMU, "p_real", 5.0) == 0
@test fmiSetReal(myFMU, ["u_real", "p_real"], [7.0, 8.0]) == 0

@test fmiSetInteger(myFMU, "p_integer", 6) == 0
@test fmiSetInteger(myFMU, ["u_integer"], [4]) == 0

@test fmiSetBoolean(myFMU, "p_boolean", false) == 0
@test fmiSetBoolean(myFMU, ["u_boolean"], [true]) == 0

@test fmiSetString(myFMU, ["p_string"], ["OldString"]) == 0
@test fmiSetString(myFMU, "p_string", "New String") == 0

fmiGetReal!(myFMU, ["p_real", "y_real", "u_real"], vR)
@test vR == [8.0, 7.0, 7.0]
@test fmiGetReal(myFMU, "p_real") == 8.0

fmiGetInteger!(myFMU, ["u_integer", "y_integer", "p_integer"], vI)
@test vI == [4, 4, 6]
@test fmiGetInteger(myFMU, "y_integer") == 4

fmiGetBoolean!(myFMU, ["u_boolean", "y_boolean", "p_boolean"], vB)
@test vB == [true, true, false]
@test fmiGetBoolean(myFMU, "p_boolean") == false

fmiGetString!(myFMU, ["p_string", "p_string"], vS)
@test vS == ["New String", "New String"]
@test fmiGetString(myFMU, "p_string") == "New String"
fmiExitInitializationMode(myFMU)

############
# Clean up #
############

fmiReset(myFMU)
fmiTerminate(myFMU)
fmiUnload(myFMU)


text:
using Dates, ReachabilityBenchmarks

# print current time
now()

# Van der Pol benchmark
include("VanDerPol/vanderpol_benchmark.jl")

# Laub-Loomis benchmark
include("LaubLoomis/laubloomis_benchmark.jl")

# Quadrotor benchmark
include("Quadrotor/quadrotor_benchmark.jl")

nothing

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
module Simulations

export Simulation, run!

import Base: show

using OrderedCollections: OrderedDict
using Oceananigans: AbstractDiagnostic, AbstractOutputWriter

using Oceananigans.Models
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.TimeSteppers
using Oceananigans.Utils

mutable struct Simulation{M, Δ, C, I, T, W, R, D, O, P, F, Π}
                 model :: M
                    Δt :: Δ
         stop_criteria :: C
        stop_iteration :: I
             stop_time :: T
       wall_time_limit :: W
              run_time :: R
           diagnostics :: D
        output_writers :: O
              progress :: P
    progress_frequency :: F
            parameters :: Π
end

"""
    Simulation(model; Δt,
         stop_criteria = Function[iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded],
        stop_iteration = Inf,
             stop_time = Inf,
       wall_time_limit = Inf,
           diagnostics = OrderedDict{Symbol, AbstractDiagnostic}(),
        output_writers = OrderedDict{Symbol, AbstractOutputWriter}(),
              progress = nothing,
    progress_frequency = 1,
            parameters = nothing)

Construct an Oceananigans.jl `Simulation` for a `model` with time step `Δt`.

Keyword arguments
=================
- `Δt`: Required keyword argument specifying the simulation time step. Can be a `Number`
  for constant time steps or a `TimeStepWizard` for adaptive time-stepping.
- `stop_criteria`: A list of functions or callable objects (each taking a single argument,
  the `simulation`). If any of the functions return `true` when the stop criteria is
  evaluated the simulation will stop.
- `stop_iteration`: Stop the simulation after this many iterations.
- `stop_time`: Stop the simulation once this much model clock time has passed.
- `wall_time_limit`: Stop the simulation if it's been running for longer than this many
   seconds of wall clock time.
- `progress`: A function with a single argument, the `simulation`. Will be called every
  `progress_frequency` iterations. Useful for logging simulation health.
- `progress_frequency`: How often to update the time step, check stop criteria, and call
  `progress` function (in number of iterations).
- `parameters`: Parameters that can be accessed in the `progress` function.
"""
function Simulation(model; Δt,
        stop_criteria = Any[iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded],
       stop_iteration = Inf,
            stop_time = Inf,
      wall_time_limit = Inf,
          diagnostics = OrderedDict{Symbol, AbstractDiagnostic}(),
       output_writers = OrderedDict{Symbol, AbstractOutputWriter}(),
             progress = nothing,
   progress_frequency = 1,
           parameters = nothing)

   if stop_iteration == Inf && stop_time == Inf && wall_time_limit == Inf
       @warn "This simulation will run forever as stop iteration = stop time " *
             "= wall time limit = Inf."
   end

   if Δt isa TimeStepWizard && progress_frequency == 1
       @warn "You have used the default progress_frequency=1. This simulation will " *
             "recalculate the time step every iteration which can be slow."
   end

   run_time = 0.0

   return Simulation(model, Δt, stop_criteria, stop_iteration, stop_time, wall_time_limit,
                     run_time, diagnostics, output_writers, progress, progress_frequency,
                     parameters)
end

function stop(sim)
    time_before = time()

    for sc in sim.stop_criteria
        if sc(sim)
            time_after = time()
            sim.run_time += time_after - time_before
            return true
        end
    end

    time_after = time()
    sim.run_time += time_after - time_before

    return false
end

function iteration_limit_exceeded(sim)
    if sim.model.clock.iteration >= sim.stop_iteration
          @info "Simulation is stopping. Model iteration $(sim.model.clock.iteration) " *
                "has hit or exceeded simulation stop iteration $(sim.stop_iteration)."
          return true
    end
    return false
end

function stop_time_exceeded(sim)
    if sim.model.clock.time >= sim.stop_time
          @info "Simulation is stopping. Model time $(prettytime(sim.model.clock.time)) " *
                "has hit or exceeded simulation stop time $(prettytime(sim.stop_time))."
          return true
    end
    return false
end

function wall_time_limit_exceeded(sim)
    if sim.run_time >= sim.wall_time_limit
          @info "Simulation is stopping. Simulation run time $(prettytime(sim.run_time)) " *
                "has hit or exceeded simulation wall time limit $(prettytime(sim.wall_time_limit))."
          return true
    end
    return false
end

get_Δt(Δt) = Δt
get_Δt(wizard::TimeStepWizard) = wizard.Δt
get_Δt(simulation::Simulation) = get_Δt(simulation.Δt)

"""
    run!(simulation)

Run a `simulation` until one of the stop criteria evaluates to true. The simulation
will then stop.
"""
function run!(sim)
    model = sim.model
    clock = model.clock

    while !stop(sim)
        time_before = time()

        if clock.iteration == 0
            [run_diagnostic(sim.model, diag) for diag in values(sim.diagnostics)]
            [write_output(sim.model, out)    for out  in values(sim.output_writers)]
        end

        for n in 1:sim.progress_frequency
            euler = clock.iteration == 0 || (sim.Δt isa TimeStepWizard && n == 1)
            time_step!(model, get_Δt(sim.Δt), euler=euler)

            [time_to_run(clock, diag) && run_diagnostic(sim.model, diag) for diag in values(sim.diagnostics)]
            [time_to_run(clock, out)  && write_output(sim.model, out)    for out  in values(sim.output_writers)]
        end

        sim.progress isa Function && sim.progress(sim)
        sim.Δt isa TimeStepWizard && update_Δt!(sim.Δt, model)

        time_after = time()
        sim.run_time += time_after - time_before
    end

    return nothing
end

Base.show(io::IO, s::Simulation) =
    print(io, "Simulation{$(typeof(s.model).name){$(typeof(s.model.architecture)), $(eltype(s.model.grid))}}\n",
            "├── Model clock: time = $(prettytime(s.model.clock.time)), iteration = $(s.model.clock.iteration) \n",
            "├── Next time step ($(typeof(s.Δt))): $(prettytime(get_Δt(s.Δt))) \n",
            "├── Progress frequency: $(s.progress_frequency)\n",
            "├── Stop criteria: $(s.stop_criteria)\n",
            "├── Run time: $(prettytime(s.run_time)), wall time limit: $(s.wall_time_limit)\n",
            "├── Stop time: $(prettytime(s.stop_time)), stop iteration: $(s.stop_iteration)\n",
            "├── Diagnostics: $(ordered_dict_show(s.diagnostics, "│"))\n",
            "└── Output writers: $(ordered_dict_show(s.output_writers, "│"))")

end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2018 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
# Function: generate_problem_dictionary
# Description: Holds simulation and model parameters as key => value pairs in a Julia Dict()
# Generated on: 2018-03-15T00:00:56.939
#
# Output arguments:
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model and simulation parameters as key => value pairs
# ----------------------------------------------------------------------------------- #
function generate_problem_dictionary()

	# Load the stoichiometric network from disk -
	path_to_network_file = joinpath(_PATH_TO_CONFIG,"Network.dat")
	stoichiometric_matrix = readdlm(path_to_network_file);

	# What is the system dimension? -
	(number_of_species,number_of_reactions) = size(stoichiometric_matrix)

	# Setup the flux bounds array -
	path_to_flux_file = joinpath(_PATH_TO_CONFIG,"flux_bounds_array.dat")
	flux_bounds_array = readdlm(path_to_flux_file);
	# TODO: update the flux_bounds_array for each reaction in your network
	# TODO: col 1 => lower bound
	# TODO: col 2 => upper bound
	# TODO: each row is a reaction

	# Setup default species bounds array -
	species_bounds_array = zeros(number_of_species,2)
	# TODO: NOTE - we by default assume Sv = 0, so species_bounds_array should be a M x 2 array of zeros
	# TODO: however, if you formulate the problem differently you *may* need to change this

	# Min/Max flag - default is minimum -
	is_minimum_flag = true

	# Setup the objective coefficient array -
	objective_coefficient_array = zeros(number_of_reactions)
	objective_coefficient_array = [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; -1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
	# TODO: update me to maximize Urea production (Urea leaving the virtual box)
	# TODO: if is_minimum_flag = true => put a -1 in the index for Urea export

	# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{String,Any}()
	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
	data_dictionary["objective_coefficient_array"] = objective_coefficient_array
	data_dictionary["flux_bounds_array"] = flux_bounds_array;
	data_dictionary["species_bounds_array"] = species_bounds_array
	data_dictionary["is_minimum_flag"] = is_minimum_flag
	data_dictionary["number_of_species"] = number_of_species
	data_dictionary["number_of_reactions"] = number_of_reactions
	# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
	return data_dictionary
end


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

d2 = OrderedDict(
  :AFF => [:ALN, :APA, :CDR],
  :AIS => [:AFF, :EGC, :SUS],
  :ALN => [:APA, :DET, :FTW, :PER, :SUS],
  :CDR => [:DET],
  :EGC => [:HOS],
  :FTW => [:DET, :EGC],
  :PER => [:DET],
  :SAN => [:AFF, :AIS, :ALN, :APA, :CDR],
  :SUS => [:EGC, :FTW, :HOS]
)

dag = DAG("unix", d2)

fname = joinpath(mktempdir(), "unix.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)


text:
function scale(reference::Array{Float64,1}, X::Array{Float64,1})
    return scale(reshape(reference, size(reference,1),1), reshape(X, size(X,1),1))
end

function scale(reference::Array{Float64,2}, X::Array{Float64,2})
    # We solve the scaling problem by min || aX - reference - b||^2,
    # where (a,b) are the scale and offset parameters

    @assert size(reference,2) == 1
    @assert size(reference) == size(X)

    n::Int64 = size(X,1)
    a::Float64 = dot(X,reference) / (dot(X,X) - sum(X)/(n * sum(reference)))
    b::Float64 = (a * sum(X) - sum(reference)) / n

    return a*X - b
end

function CCC(X::Array{Float64,2}, Y::Array{Float64,2})::Float64
    if size(X,2) == 1 & size(Y,2) == 1
        X = reshape(X, size(X,1), )
        Y = reshape(Y, size(Y,1), )
        return CCC(X,Y)
    else
        warn("Check dimensions")
    end
end

function CCC(X::Array{Float64,1}, Y::Array{Float64,1})::Float64
    ρ = 2cov(X, Y, corrected=true) / (var(X) + var(Y) + (mean(X) - mean(Y)) ^ 2)
end

function CCCweights(annotations::DataFrame)
    # Precompute the number of annotators
    nannotators = size(annotations, 2)

    CCCs = zeros(Float64, nannotators, nannotators)

    for i in 1:nannotators, j in 1:nannotators
        CCCs[i,j] = CCC(annotations[:,i], annotations[:,j])
    end

    return sum(CCCs, dims=2) .- 1
end

function agreements(triplets1::T, triplets2::T) where T <: Triplets

    n1, n2 = maximum(triplets1), maximum(triplets2)

    n1 == n2 || throw(ArgumentError("Number of items in each triplet set does not match"))

    intersection = intersect(Set(triplets1), Set(triplets2))

    return length(intersection) / ( n1 * binomial(n1 - 1, 2))
end

function dense(data::DataFrame, raters::Symbol, items::Symbol, ratings::Symbol)
    groups = groupby(select(data, [items, raters, ratings]), raters)
    df = select(convert(DataFrame, groups[1]), [items, ratings])

    for i = 2:length(groups)
        aux = select(convert(DataFrame, groups[i]), [items, ratings])
        df = outerjoin(df, aux, on=items, makeunique=true)
    end

    annotators = unique(data[!, raters])
    rename!(df, [items; Symbol.(annotators)])

    return df
end

function columncounts(df::AbstractDataFrame, column::Symbol)
    return Dict([(i, count(x -> x == i, df[!,column])) for i in unique(df[!,column])])
end

"""
    function anonymize(df::DataFrame, raters::Symbol)

Anonymize the raters' names in columns `raters`. Returns
the same dataframe `df` with hashed values in the
`raters` column.
"""
function anonymize(df::DataFrame, raters::Symbol)
    annotators = unique(df[!,raters])
    hashed = Dict(annotators .=> hashsha256.(annotators))
    df[!,raters] = [hashed[annotator] for annotator in df[!,raters]]
    return df
end

function anonymize!(df::DataFrame, raters::Symbol)
    annotators = unique(df[!,raters])
    hashed = Dict(annotators .=> hashsha256.(annotators))
    df[!,raters] = [hashed[annotator] for annotator in df[!,raters]]
end

function hashsha256(s::AbstractString; nchars::Int = 8)
    return bytes2hex(sha256(s))[1:nchars]
end

function anonymize(df::DataFrame)
    annotators = names(df)
    hashed = Dict(annotators .=> hashsha256.(annotators))
    rename!(df, hashed)
    return df
end

"""
    function name(method::FusionMethod)

Return the name of an annotation fusion method from its signature.

This method should be used to create column names in a DataFrame.
"""
function name(method::FusionMethod)
    if :scaling in fieldnames(typeof(method))
        return lowercase(string(typeof(method), "_", method.scaling))
    else
        return lowercase(string(typeof(method)))
    end
end

"""
    function fillmissing!(row, value::T) where T <: Real

Fill missing values for a vector, using `value` to replace the missing values..
"""
function fillmissing!(row, value::T) where T <: Real
    for i in eachindex(row)
        row[i] = ismissing(row[i]) ? value : row[i]
    end
end

"""
    function fillmissing!(df::DataFrame, f::Function)

Fill missing values per row, using the function f applied row-wise.
"""
function fillmissing(df::DataFrame, index::Symbol, f::Function)
    imputed = copy(df)

    for name in names(imputed[!, Not(index)])
        imputed[!, name] = convert(Vector{Union{Float64,Missing}}, imputed[!, name])
    end

    for row in eachrow(imputed)
        fillmissing!(row[Not(index)], f(row[Not(index)]))
    end

    return imputed
end

function fillmissing(df::AbstractMatrix, f::Function)
    imputed = convert(Matrix{Union{Missing, Float64}}, df)

    for row in eachrow(imputed)
        fillmissing!(row, f(row))
    end

    return imputed
end