text:
using Documenter, CovidData

makedocs(;
    modules=[CovidData],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/schrimpf/CovidData.jl/blob/{commit}{path}#L{line}",
    sitename="CovidData.jl",
    authors="Paul Schrimpf <paul.schrimpf@gmail.com>"
         )

deploydocs(repo="github.com/schrimpf/CovidData.jl.git")


text:
using Lazy
using StatsBase
using Random

import Base.get, Base.get!

ext: jl


text:
# julia --project -O3 --check-bounds=no ./scripts-part1/diffusion3D_work_precision.jl

# Initialize Globals
const USE_GPU = true
const BENCHMARK = false
const VISUALIZE = false
const STEADY = true

# Include main script
include("diffusion3D_xpu.jl")

# Grid sizes for the work precision diagram
grid_sizes = 16 * 2 .^ (1:5)
iters = zeros(5)

# local point for observation
local_point = (5,5,5)
lp_vals = zeros(5)

# loop over grid sizes
for it = 1:5
    # run simulation with current grid size
    iters[it], xc, H = diffusion_3D(grid_sizes[it])

    # get the index of the nearest local point of interest
    lp_index = findmin(abs.(xc.-local_point[1]))[2], findmin(abs.(xc.-local_point[2]))[2], findmin(abs.(xc.-local_point[3]))[2]
    # save the value
    lp_vals[it] = H[lp_index[1],lp_index[2],lp_index[3]]
end

# get the well converged solution
iter, xc, H_wc = diffusion_3D(128; tol=1e-24)

# tolerances for the work precision diagram
tols_exp = (6:23)
tols = (1.0)./((10.0).^tols_exp)

# allocate arrays for convergence behavior
rel_conv = zeros(length(tols))
abs_conv = zeros(length(tols))

# loop over tolerances
for it = 1:length(tols)
    # get the solution with current tolerance
    iter, xc, H_c = diffusion_3D(128; tol=tols[it])

    # Calculate absolute difference
    abs_diff = H_wc - H_c
    abs_conv[it] = norm(abs_diff)

    # Calculate relative difference
    rel_diff = norm(abs_diff)/norm(H_wc)
    rel_conv[it] = rel_diff
end

# Visualization of Results
gr()
ENV["GKSwstype"]="nul"

# Plot: Number of grid points VS. Iterations
plot(grid_sizes.^3, iters, xlabel="Number of grid points", xaxis=:log, ylabel="Iterations", label=false)
png("$(@__DIR__)/../docs/img/diffusion3D_workprecision_1")

plot(grid_sizes, iters, xlabel="cbrt(Number of grid points)", ylabel="Iterations", label=false)
png("$(@__DIR__)/../docs/img/diffusion3D_workprecision_2")

# Plot: Number of grid points VS. Value at local_point
plot(grid_sizes.^3, lp_vals, xlabel="Number of grid points", xaxis=:log, ylabel="Value at $(local_point)", label=false)
png("$(@__DIR__)/../docs/img/diffusion3D_workprecision_3")

plot(grid_sizes, lp_vals, xlabel="cbrt(Number of grid points)", ylabel="Value at $(local_point)", label=false)
png("$(@__DIR__)/../docs/img/diffusion3D_workprecision_4")

# Plot: Tolerance VS. Relative/Absolute difference to well converged solution
plot(tols_exp, rel_conv, xlabel="tolerance [1e-x]", ylabel="rel diff to well converged sol", yaxis=:log, label=false)
png("$(@__DIR__)/../docs/img/diffusion3D_workprecision_5")

plot(tols_exp, abs_conv, xlabel="tolerance [1e-x]", ylabel="abs diff to well converged sol", yaxis=:log, label=false)
png("$(@__DIR__)/../docs/img/diffusion3D_workprecision_6")


text:
function outer_vals(v::Vector{<:Real}, bounds)
    length(bounds) == 2 || error("bounds must contain 2 values")
    (bounds[1] < v[1] || bounds[2] > v[end]) && error("$bounds is outside of vector range")
    lower = v[findfirst(x -> x > bounds[1], v)-1]
    upper = v[findfirst(x -> x >= bounds[2], v)]
    lower, upper
end

function inner_vals(v::Vector{<:Real}, bounds)
    length(bounds) == 2 || error("bounds must contain 2 values")
    (bounds[1] < v[1] || bounds[2] > v[end]) && error("$bounds is outside of vector range")
    lower = v[findfirst(x -> x >= bounds[1], v)]
    upper = v[findfirst(x -> x > bounds[2], v)-1]
    lower, upper
end

function inner_area(area::Vector{<:Real}, gridres::Real)
    alons = -180.0:gridres:180.0 |> collect
    alats = -90.0:gridres:90.0 |> collect

    inlons = inner_vals(alons, (area[2], area[4]))
    inlats = inner_vals(alats, (area[3], area[1]))
    round.([inlats[2], inlons[1], inlats[1], inlons[2]], digits=6)
end

function outer_area(area::Vector{<:Real}, gridres::Real)
    alons = -180.0:gridres:180.0 |> collect
    alats = -90.0:gridres:90.0 |> collect

    inlons = outer_vals(alons, (area[2], area[4]))
    inlats = outer_vals(alats, (area[3], area[1]))
    round.([inlats[2], inlons[1], inlats[1], inlons[2]], digits=6)
end

round_area(area::Vector{<:Real}, mult=1) = return [ceil(area[1]*mult)/mult, floor(area[2]*mult)/mult, floor(area[3]*mult)/mult, ceil(area[4]*mult)/mult]


text:
#=##############################################################################
# DESCRIPTION
    Plotting and outputting-processing functions for monitoring vehicle and
    component performance. The functions in this module generate
    monitor-functions that can be passed to the simulation engine through the
    optional arguments `extra_runtime_function`.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

"""
    `generate_monitor_rotors( rotors::Array{vlm.Rotor, 1}, J_ref::Real,
rho_ref::Real, RPM_ref::Real, nsteps_sim::Int)`

Generate monitor of load distribution and performance of a rotor system.
`J_ref` and `rho_ref` are the reference advance ratio and air density used for
calculating propulsive efficiency and coefficients. `RPM_ref` is the reference
RPM used to estimate the age of the wake. `nsteps_sim` is the number of time
steps by the end of the simulation (used for generating the color gradient).
"""
function generate_monitor_rotors( rotors::Array{vlm.Rotor, 1},
                                    J_ref::Real, rho_ref::Real, RPM_ref::Real,
                                    nsteps_sim::Int;
                                    t_scale=1.0,                    # Time scaling factor
                                    t_lbl="Simulation time (s)",    # Time-axis label
                                    # OUTPUT OPTIONS
                                    save_path=nothing,
                                    run_name="rotor",
                                    figname="monitor_rotor",
                                    disp_conv=true,
                                    conv_suff="_convergence.csv",
                                    nsteps_savefig=10,
                                    colors="rgbcmy"^100,
                                    stls="o^*.px"^100, )

    fcalls = 0                  # Number of function calls

    # Name of convergence file
    if save_path!=nothing
        fname = joinpath(save_path, run_name*conv_suff)
    end

    # Call figure
    if disp_conv
        fig = figure(figname, figsize=(7*3,5*2))
        axs = fig.subplots(2, 3)
        axs = [axs[1], axs[3], axs[5], axs[2], axs[4], axs[6]]
    end

    # Function for run_vpm! to call on each iteration
    function extra_runtime_function(sim::Simulation{V, M, R},
                                    PFIELD::vpm.ParticleField,
                                    T::Real, DT::Real
                                   ) where{V<:AbstractVLMVehicle, M, R}

        # rotors = vcat(sim.vehicle.rotor_systems...)
        angle = T*360*RPM_ref/60
        t_scaled = T*t_scale

        if fcalls==0
            # Format subplots
            if disp_conv
                ax = axs[1]
                ax.title.set_text("Circulation Distribution")
                ax.set_xlabel("Element index")
                ax.set_ylabel(L"Circulation $\Gamma$ (m$^2$/s)")
                ax.grid(true, color="0.8", linestyle="--")
                ax = axs[2]
                ax.title.set_text("Plane-of-rotation Normal Force")
                ax.set_xlabel("Element index")
                ax.set_ylabel(L"Normal Force $N_p$ (N/m)")
                ax.grid(true, color="0.8", linestyle="--")
                ax = axs[3]
                ax.title.set_text("Plane-of-rotation Tangential Force")
                ax.set_xlabel("Element index")
                ax.set_ylabel(L"Tangential Force $T_p$ (N/m)")
                ax.grid(true, color="0.8", linestyle="--")
                ax = axs[4]
                ax.set_xlabel(t_lbl)
                ax.set_ylabel(L"Thrust $C_T$")
                ax.grid(true, color="0.8", linestyle="--")
                ax = axs[5]
                ax.set_xlabel(t_lbl)
                ax.set_ylabel(L"Torque $C_Q$")
                ax.grid(true, color="0.8", linestyle="--")
                ax = axs[6]
                ax.set_xlabel(t_lbl)
                ax.set_ylabel(L"Propulsive efficiency $\eta$")
                ax.grid(true, color="0.8", linestyle="--")
            end

            # Convergence file header
            if save_path!=nothing
                f = open(fname, "w")
                print(f, "ref age (deg),T,DT")
                for (i, rotor) in enumerate(rotors)
                    print(f, ",RPM_$i,CT_$i,CQ_$i,eta_$i")
                end
                print(f, "\n")
                close(f)
            end
        end

        # Write rotor position and time on convergence file
        if save_path!=nothing
            f = open(fname, "a")
            print(f, angle, ",", T, ",", DT)
        end

        # Plot circulation and loads distributions
        if disp_conv

            cratio = PFIELD.nt/nsteps_sim
            cratio = cratio > 1 ? 1 : cratio
            clr = fcalls==0 && false ? (0,0,0) : (1-cratio, 0, cratio)
            stl = fcalls==0 && false ? "o" : "-"
            alpha = fcalls==0 && false ? 1 : 0.5

            # Circulation distribution
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, [vlm.get_blade(rotor, j).sol["Gamma"] for j in 1:rotor.B]...)
            end
            axs[1].plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)

            # Np distribution
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, rotor.sol["Np"]["field_data"]...)
            end
            axs[2].plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)

            # Tp distribution
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, rotor.sol["Tp"]["field_data"]...)
            end
            axs[3].plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)
        end

        # Plot performance parameters
        for (i,rotor) in enumerate(rotors)
            CT, CQ = vlm.calc_thrust_torque_coeffs(rotor, rho_ref)
            eta = J_ref*CT/(2*pi*CQ)

            if disp_conv
                axs[4].plot([t_scaled], [CT], "$(stls[i])", alpha=alpha, color=clr)
                axs[5].plot([t_scaled], [CQ], "$(stls[i])", alpha=alpha, color=clr)
                axs[6].plot([t_scaled], [eta], "$(stls[i])", alpha=alpha, color=clr)
            end

            if save_path!=nothing
                print(f, ",", rotor.RPM, ",", CT, ",", CQ, ",", eta)
            end
        end

        if disp_conv
            # Save figure
            if fcalls%nsteps_savefig==0 && fcalls!=0 && save_path!=nothing
                fig.savefig(joinpath(save_path, run_name*"_convergence.png"),
                                                    transparent=false, dpi=300)
            end
        end

        # Close convergence file
        if save_path!=nothing
            print(f, "\n")
            close(f)
        end

        fcalls += 1

        return false
    end

    return extra_runtime_function
end

"""
    `generate_monitor_wing(wing, Vinf::Function, b_ref::Real, ar_ref::Real,
rho_ref::Real, qinf_ref::Real, nsteps_sim::Int)`

Generate monitor of load distribution and performance of a wing.
"""
function generate_monitor_wing(wing, Vinf::Function, b_ref::Real, ar_ref::Real,
                                rho_ref::Real, qinf_ref::Real, nsteps_sim::Int;
                                lencrit_f=0.5,      # Factor for critical length to ignore horseshoe forces
                                L_dir=[0,0,1],      # Direction of lift component
                                D_dir=[-1,0,0],     # Direction of drag component
                                # OUTPUT OPTIONS
                                out_Lwing=nothing,
                                out_Dwing=nothing,
                                out_CLwing=nothing,
                                out_CDwing=nothing,
                                save_path=nothing,
                                run_name="wing",
                                figname="monitor_wing",
                                disp_plot=true,
                                title_lbl="",
                                CL_lbl=L"Lift Coefficient $C_L$",
                                CD_lbl=L"Drag Coefficient $C_D$",
                                y2b_i=2,
                                y2b_ref=nothing,
                                ClCL_ref=nothing, CdCD_ref=nothing,
                                CL_ref=nothing, CD_ref=nothing,
                                ref_lbl="Reference", ref_stl="ok",
                                conv_suff="_convergence.csv",
                                figsize_factor=5/6,
                                nsteps_plot=10,
                                nsteps_savefig=10)

    fcalls = 0                  # Number of function calls

    # y2b = 2*wing._ym/b_ref          # Initial span position
    y2b = 2*[vlm.getControlPoint(wing, i)[y2b_i] for i in 1:vlm.get_m(wing)]/b_ref
    prev_wing = nothing

    # Critical length to ignore horseshoe forces
    meanchord = b_ref/ar_ref
    lencrit = lencrit_f*meanchord/vlm.get_m(wing)

    # Name of convergence file
    if save_path!=nothing
        fname = joinpath(save_path, run_name*conv_suff)
    end

    if disp_plot
        fig1 = figure(figname, figsize=[7*2, 5*2]*figsize_factor)
        fig1.suptitle(title_lbl)
        axs1 = fig1.subplots(2, 2)
        axs1 = [axs1[1], axs1[3], axs1[2], axs1[4]]
        ax = axs1[1]
        # xlim([-1,1])
        ax.set_xlabel(L"$\frac{2y}{b}$")
        ax.set_ylabel(L"$\frac{Cl}{CL}$")
        ax.title.set_text("Spanwise lift distribution")

        ax = axs1[2]
        # xlim([-1,1])
        ax.set_xlabel(L"$\frac{2y}{b}$")
        ax.set_ylabel(L"$\frac{Cd}{CD}$")
        ax.title.set_text("Spanwise drag distribution")

        ax = axs1[3]
        ax.set_xlabel("Simulation time (s)")
        ax.set_ylabel(CL_lbl)

        ax = axs1[4]
        ax.set_xlabel("Simulation time (s)")
        ax.set_ylabel(CD_lbl)

        fig2 = figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
        axs2 = fig2.subplots(1, 2)
        ax = axs2[1]
        ax.set_xlabel(L"$\frac{2y}{b}$")
        ax.set_ylabel(L"Circulation $\Gamma$")
        ax = axs2[2]
        ax.set_xlabel(L"$\frac{2y}{b}$")
        ax.set_ylabel(L"Effective velocity $V_\infty$")
    end

    function extra_runtime_function(sim, PFIELD, T, DT)

        aux = PFIELD.nt/nsteps_sim
        clr = (1-aux, 0, aux)

        if fcalls==0
            # Convergence file header
            if save_path!=nothing
                f = open(fname, "w")
                print(f, "T,CL,CD\n")
                close(f)
            end
        end

        if PFIELD.nt>2

            # Force at each VLM element
            Ftot = calc_aerodynamicforce(wing, prev_wing, PFIELD, Vinf, DT,
                                                            rho_ref; t=PFIELD.t,
                                                            lencrit=lencrit)
            L, D, S = decompose(Ftot, L_dir, D_dir)
            vlm._addsolution(wing, "L", L)
            vlm._addsolution(wing, "D", D)
            vlm._addsolution(wing, "S", S)

            # Force per unit span at each VLM element
            ftot = calc_aerodynamicforce(wing, prev_wing, PFIELD, Vinf, DT,
                                        rho_ref; t=PFIELD.t, per_unit_span=true,
                                        lencrit=lencrit)
            l, d, s = decompose(ftot, L_dir, D_dir)
            vlm._addsolution(wing, "l", l)
            vlm._addsolution(wing, "d", d)
            vlm._addsolution(wing, "s", s)

            # Lift of the wing
            Lwing = norm(sum(L))
            CLwing = Lwing/(qinf_ref*b_ref^2/ar_ref)
            ClCL = norm.(l) / (Lwing/b_ref)

            # Drag of the wing
            Dwing = norm(sum(D))
            CDwing = Dwing/(qinf_ref*b_ref^2/ar_ref)
            CdCD = [sign(dot(this_d, [1,0,0])) for this_d in d].*norm.(d) / (Dwing/b_ref) # Preserves the sign of drag

            vlm._addsolution(wing, "Cl/CL", ClCL)
            vlm._addsolution(wing, "Cd/CD", CdCD)

            if out_Lwing!=nothing; push!(out_Lwing, Lwing); end;
            if out_Dwing!=nothing; push!(out_Dwing, Dwing); end;
            if out_CLwing!=nothing; push!(out_CLwing, CLwing); end;
            if out_CDwing!=nothing; push!(out_CDwing, CDwing); end;

            if PFIELD.nt%nsteps_plot==0 && disp_plot

                ax = axs1[1]
                if y2b_ref!=nothing && ClCL_ref!=nothing
                    ax.plot(y2b_ref, ClCL_ref, ref_stl, label=ref_lbl)
                end
                ax.plot(y2b, ClCL, "-", alpha=0.5, color=clr)

                ax = axs1[2]
                if y2b_ref!=nothing && CdCD_ref!=nothing
                    ax.plot(y2b_ref, CdCD_ref, ref_stl, label=ref_lbl)
                end
                ax.plot(y2b, CdCD, "-", alpha=0.5, color=clr)

                ax = axs1[3]
                if CL_ref!=nothing
                    ax.plot([0, T], CL_ref*ones(2), ":k", label=ref_lbl)
                end
                ax.plot([T], [CLwing], "o", alpha=0.5, color=clr)

                ax = axs1[4]
                if CD_ref!=nothing
                    ax.plot([0, T], CD_ref*ones(2), ":k", label=ref_lbl)
                end
                ax.plot([T], [CDwing], "o", alpha=0.5, color=clr)

                ax = axs2[1]
                ax.plot(y2b, wing.sol["Gamma"], "-", alpha=0.5, color=clr)
                ax = axs2[2]
                if "Vkin" in keys(wing.sol)
                    ax.plot(y2b, norm.(wing.sol["Vkin"]), "-", alpha=0.5, color=[clr[1], 1, clr[3]])
                end
                if "Vvpm" in keys(wing.sol)
                    ax.plot(y2b, norm.(wing.sol["Vvpm"]), "-", alpha=0.5, color=clr)
                elseif "Vind" in keys(wing.sol)
                    ax.plot(y2b, norm.(wing.sol["Vind"]), "-", alpha=0.5, color=clr)
                end
                ax.plot(y2b, [norm(Vinf(vlm.getControlPoint(wing, i), T)) for i in 1:vlm.get_m(wing)],
                                                            "-k", alpha=0.5)

                if save_path!=nothing
                    # Save figures
                    if PFIELD.nt%nsteps_savefig==0
                        fig1.savefig(joinpath(save_path, run_name*"_convergence.png"),
                                                    transparent=false, dpi=300)

                        fig2.savefig(joinpath(save_path, run_name*"_convergence2.png"),
                                                    transparent=false, dpi=300)
                    end
                end
            end

            if save_path != nothing
                # Write rotor position and time on convergence file
                f = open(fname, "a")
                print(f, T, ",", CLwing, ",", CDwing, "\n")
                close(f)
            end
        end

        prev_wing = deepcopy(wing)
        fcalls += 1
        return false
    end
end

function generate_monitor_statevariables(; figname="monitor_statevariables",
                                           save_path=nothing,
                                           run_name="",
                                           nsteps_savefig=10)

    fig = figure(figname, figsize=[7*2, 5*1])
    axs = fig.subplots(1, 3)
    ax = axs[1]
    ax.set_xlabel("Simulation time")
    ax.set_ylabel("Velocity")
    Vlbls = [L"V_x", L"V_y", L"V_z"]
    ax = axs[2]
    ax.set_xlabel("Simulation time")
    ax.set_ylabel(L"Angular velocity ($^\circ/t$)")
    Wlbls = [L"\Omega_x", L"\Omega_y", L"\Omega_z"]
    ax = axs[3]
    ax.set_xlabel("Simulation time")
    ax.set_ylabel(L"$O$ position")
    Olbls = [L"O_x", L"O_y", L"O_z"]

    function extra_runtime_function(sim, PFIELD, T, DT)
        for j in 1:3
            axs[1].plot(sim.t, sim.vehicle.V[j], ".", label=Vlbls[j], alpha=0.8,
                                                            color=clrs[j])
            axs[2].plot(sim.t, sim.vehicle.W[j], ".", label=Wlbls[j], alpha=0.8,
                                                            color=clrs[j])
            axs[3].plot(sim.t, sim.vehicle.system.O[j], ".", label=Olbls[j], alpha=0.8,
                                                            color=clrs[j])
        end
        if sim.nt==0
            for j in 1:2
                axs[j].legend(loc="best", frameon=false)
            end
        end

        if save_path!=nothing
            # Save figures
            if PFIELD.nt%nsteps_savefig==0
                fig.savefig(joinpath(save_path, run_name*"statevariables.png"),
                                                    transparent=false, dpi=300)
            end
        end

        return false
    end

    return extra_runtime_function
end

function generate_monitor_enstrophy(; save_path=nothing,
                                     figname="monitor_enstrophy", run_name="",
                                     nsteps_savefig=10,
                                     nsteps_plot=10)

    fig = figure(figname, figsize=[7*1, 5*1])
    ax = fig.gca()
    ax.set_xlabel("Simulation time")
    ax.set_ylabel(L"Enstrophy ($\mathrm{m}^6/\mathrm{s}^2$)")

    enstrophy = []
    ts = []

    function extra_runtime_function(sim, PFIELD, T, DT)
        vpm.monitor_enstrophy(PFIELD, T, DT; save_path=save_path,
                                               run_name=run_name, out=enstrophy)

        push!(ts, T)

        if PFIELD.nt%nsteps_plot == 0
            ax.plot(ts, enstrophy[end-length(ts)+1:end], ".k")
            ts = []
        end

        if save_path!=nothing
            if PFIELD.nt%nsteps_savefig==0
                fig.savefig(joinpath(save_path, run_name*"enstrophy.png"),
                                                    transparent=false, dpi=300)
            end
        end

        return false
    end

    return extra_runtime_function
end

max_stars_repo_name: byuflowlab/FlightVehicleSim

max_issues_repo_name: password-now/FLOWUnsteady

max_forks_repo_name: password-now/FLOWUnsteady


text:
using DataFrames
using CSV, DelimitedFiles
using JuMP, CPLEX

#Preparing optimization model
model = Model(CPLEX.Optimizer)
set_optimizer_attribute(model, "CPX_PARAM_EPINT", 1e-8)

const N_GEN = 3
const N_HRS = 4

dfGenerator = CSV.read("data_generators.csv", DataFrame)
dfDemand = CSV.read("demand_reserves.csv", DataFrame)

#println(dfGenerator);
#println(dfDemand);
for g=1:N_GEN
  println("\t Generator at Time $g: ", dfGenerator.Uini[g])
end

@variable(model, genOnOff[1:N_GEN, 1:N_HRS], Bin)
@variable(model, genOut[1:N_GEN, 1:N_HRS] >= 0)
@variable(model, genStartUp[1:N_GEN, 1:N_HRS], Bin)
@variable(model, genShutDown[1:N_GEN, 1:N_HRS], Bin)

#Setting the objective function
@objective(model, Min, sum(sum(dfGenerator.FixedCost[g]*genOnOff[g,t] + dfGenerator.VariableCost[g]*genOut[g,t] + dfGenerator.StartUpCost[g]*genStartUp[g,t] + dfGenerator.ShutdownCost[g]*genShutDown[g, t] for g in 1:N_GEN) for t in 1:N_HRS) )

#Initial constraints
#Generators should have an initial and defined value t=0
@constraint(model, constInit[g=1:N_GEN], genOnOff[g,1] == dfGenerator.Uini[g])

#Initial values of Generator shutdown and startup values should be defined
@constraint(model, constShutDn[g=1:N_GEN], genShutDown[g,1] == 0)
@constraint(model, constShutUp[g=1:N_GEN], genStartUp[g,1] == 0)

#Demand constraint
#for t in 1:N_HRS
@constraint(model, constDemand[t=1:N_HRS], sum(genOut[g, t] for g=1:N_GEN) == dfDemand.DemandNode01[t])
#end

#Power bounds
@constraint(model, constMaxPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] <= dfGenerator.MaxPowerOut[g]*genOnOff[g,t] )
@constraint(model, constMinPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] >= dfGenerator.MinPowerOut[g]*genOnOff[g,t] )

#Logical Conditions
@constraint(model, constLog1[t=2:N_HRS, g=1:N_GEN], (genStartUp[g,t] - genShutDown[g,t]) == (genOnOff[g,t] - genOnOff[g,t-1]))
@constraint(model, constLog2[t=1:N_HRS, g=1:N_GEN], genStartUp[g,t] + genShutDown[g, t] <= 1)

#Ramp limits
@constraint(model, constRampUp[t=2:N_HRS, g=1:N_GEN], genOut[g,t] - genOut[g,t-1] <= dfGenerator.RampUpLimit[g]*genOnOff[g, t-1] + dfGenerator.RampStartUpLimit[g]*genStartUp[g,t])
@constraint(model, constRampDown[t=2:N_HRS, g=1:N_GEN], (genOut[g,t-1] - genOut[g,t]) <= (dfGenerator.RampDownLimit[g]*genOnOff[g,t] + dfGenerator.RampShutDownLimit[g]*genShutDown[g,t]))

#Reserves
@constraint(model, constReserves[t=1:N_HRS], sum(dfGenerator.MaxPowerOut[g]*genOnOff[g,t] for g=1:N_GEN) >= (dfDemand.DemandNode01[t] + dfDemand.Reserves[t]))

optimize!(model)

println("------------------------------------")
println("------- OPTIMAL SOLUTION -------")
for t=1:N_HRS
  print("Hour $t","\t")
end
println("")
for g=1:N_GEN
  for t=1:N_HRS
    print(round(JuMP.value(genOut[g,t]); digits=3 ),"\t")
  #  println(round(pi; digits = 3))

  end
  println("\t Generator: $g")
end

println("------------------------------------")
println("------- GENERATORS COMMITED --------")
for t=1:N_HRS
  print("Hour $t","\t")
end
println("")

for g=1:N_GEN
  for t=1:N_HRS
    #print(JuMP.value(genOnOff[g,t]),"\t")
    print(convert(Int32, JuMP.value(genOnOff[g,t])),"\t")
  end
  println("\t Generator: $g")
end

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
module TestShow

include("preamble.jl")

xforms = [
    Cat(),
    Count(),
    Count(2),
    Count(3, 4),
    GetIndex([0]),
    # GetIndex{true}([0]),
    Iterated(sqrt, 0.1),
    # Iterated(sqrt, 0.1, Number),
    Map(sin),
    Filter(isfinite),
    Scan(*),
    TeeZip(Filter(isfinite)),
    TeeZip(Filter(isfinite) |> Map(sin)),
    TeeZip(Filter(isfinite) |> Map(sin)) |> Map(sum),
    let xf = TeeZip(Filter(isfinite) |> Map(sin))
        xf |> Map(sum) |> xf
    end,
    let xf = Map(first) |> Map(last)
        xf = TeeZip(TeeZip(xf) |> Map(identity)) |> xf
        xf |> TeeZip(xf) |> xf
    end,
    # Zip(Map(sin), Map(cos), Map(tan)) |> Map(prod),
]

@testset "smoke test summary(xf)" begin
    for xf in xforms
        @test sprint(summary, xf) isa String
    end
end

# Note: Use `ENV["JULIA_DEBUG"] = "Main"` to enable debugging
@testset "eval(show(xf))" begin
    @testset "$(summary(xf))" for xf in xforms
        code = sprint(show, xf)
        @debug """
        show(xf) =
        $code
        """
        xf2 = include_string(@__MODULE__, code)
        @test xf == xf2
    end
end
@testset "eval(show(text/plain, xf))" begin
    @testset "$(summary(xf))" for xf in xforms
        code = sprint(show, "text/plain", xf)
        @debug """
        show("text/plain", xf) =
        $code
        """
        xf2 = include_string(@__MODULE__, code)
        @test xf == xf2
    end
end

end  # module


text:
@testset "Properties" begin
    @test iszero(GF256(0))
    @test !iszero(GF256(1))
    @test !iszero(GF256(255))
    @test isone(GF256(1))
    @test !isone(GF256(0))
    @test !isone(GF256(255))
    @test zero(GF256) == GF256(0)
    @test zero(GF256) == zero(GF256(1))
    @test one(GF256) == GF256(1)
    @test one(GF256) == one(GF256(0))
end

@testset "Logarithms and exponentiation" begin
    a = GF256(0)
    @test log(exp(a)) == a
    for x in 1:254
        a = GF256(x)
        @test log(exp(a)) == a
        @test exp(log(a)) == a
    end
    a = GF256(255)
    @test exp(log(a)) == a
end

@testset "Division and multiplication" begin
    β = GF256(43) # chosen arbitrarily
    @test β * zero(GF256) == zero(GF256)
    @test β * one(GF256) == β
    for x in 1:255
        @test (β * GF256(x)) / GF256(x) == β
        @test (β / GF256(x)) * GF256(x) == β
    end
end

@testset "Addition and subtraction" begin
    β = GF256(43) # chosen arbitrarily
    @test β + zero(GF256) == β
    @test β - zero(GF256) == β
    for x in 0:255
        a = GF256(x)
        @test iszero(a+a)
        @test iszero(a-a)
    end
end

# TODO: remove
# @test GF256(2) + GF256(10) - GF256(10) == GF256(2)
# @test GF256(2) * GF256(10) / GF256(10) == GF256(2)
# @test GF256(1) + GF256(2) * GF256(10) - GF256(10) - GF256(10) == GF256(21)
# @test GF256(13) * GF256(1) == GF256(13)
# @test GF256(13) / GF256(1) == GF256(13)
# @test GF256(0) * GF256(1) == GF256(0)
# @test GF256(0) / GF256(1) == GF256(0)
# @test [GF256(i) for i in 0:255].*zero(GF256) == zeros(GF256, 256)
# @test [GF256(i) for i in 0:255].*one(GF256) == [GF256(i) for i in 0:255]
# @test [GF256(i) for i in 0:255].*GF256(2) == [GF256(2)*GF256(i) for i in 0:255]
# # @test GF256(11) + false == GF256(11)
# # @test GF256(11) + true == GF256(11) + GF256(1)
# # @test GF256(11) * false == GF256(0)
# # @test GF256(11) * true == GF256(11)
# # @test GF256(11) / true == GF256(11)
# @test GF256(123)^0 == GF256(1)
# @test GF256(123)^1 == GF256(123)
# @test GF256(123)^2 == GF256(123) * GF256(123)
# @test GF256(46)^3 == GF256(46) * GF256(46) * GF256(46)
# @test GF256(10)^20 == reduce(*, (GF256(10) for _ in 1:20))
# @test GF256(56)^30 == reduce(*, (GF256(56) for _ in 1:30))

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
using Base.Test
using NeuralNets
using PyPlot

include("../apma2821v/plotting.jl")

# very strong all-to-all connections
function strong(pre_N, post_N, stype::SynapseType, stdp_enabled)
    syn = Synapse[]
    for i=1:(pre_N*post_N)
        w = 1
        delay = 5
        (s,t) = ind2sub((pre_N,post_N),i)
        s = Synapse(stype[i],s,t,float(w),float(delay))
        push!(syn,s)
    end
    SpikingConnection(0,0,syn,stdp_enabled,20.0,10.0,5.0)
end

# only stim the E1 neuron
function E1_stim(l::Int64, N::Int64, i::Int64)
    Iinj = zeros(N)
    if l==1
        dt = 0.01
        Iinj = 30 * (int(i*dt/20)%2)* ones(N)
    end
    Iinj
end

function genplots(out::RunResults, time_trace)
    figure()
    subplot(1,2,1)
    lnames=["E1","E2"]
    ax = cell(2)
    for (i,lname) in enumerate(lnames)
        v = out.data[lname]["v"]'
        didSpike = out.data[lname]["didSpike"]'
        ax[i] = plot(time_trace,v)
        plot(time_trace[find(didSpike)],v[find(didSpike)],"r.")
    end
    legend([ax[1],ax[2]],lnames)
    ylim([-100,10])
    title("v")

    subplot(1,2,2)
    for (i,lname) in enumerate(lnames)
        g = out.data[lname]["gsyn"]'
        plot(time_trace,g)
        #ylim([-10,10])
        title("Presynaptic Conductance")
    end
    show()
end

T = 300
dt = 0.01
time_trace = 0:dt:T
measurements = [E=>["v","didSpike","gsyn"] for E in ["E1","E2"]]
settings = RunSettings(T,dt,E1_stim,measurements)
plot_stim(settings,["E1","E2"])

# no STDP
net = SpikingNet()
add_layer!(net,HHLayer("E1",1))
add_layer!(net,HHLayer("E2",1))
connect!(net, "E1", "E2", strong, EXC, false)
out = runSpikingNet(net,settings)
genplots(out,time_trace)

@test out.data["E1"]["v"] != out.data["E2"]["v"]

# now, enable STDP
# Now we enable STDP in the single synapse from E1->E2. STDP is triggered
# on a synapse when either it's source or target neuron fires. When green (above) spikes,
# the duration Δt = green - blue is short, so LTP is strong. When blue spikes,
# Δt = green - blue < 0, but since the duration is long, the LTP effect is weak.
# There is no reciprocal connection from E2->E1, so potentiation only occurs on one synapse.

net2 = SpikingNet()
add_layer!(net2,HHLayer("E1",1))
add_layer!(net2,HHLayer("E2",1))
connect!(net2, "E1", "E2", strong, EXC, true)

println("starting weight: ", net2.pre[2][1].synapses[1].w)
out = runSpikingNet(net2,settings)
println("ending weight: ", net2.pre[2][1].synapses[1].w)

genplots(out,time_trace)

println("Press Enter to Continue...")
readline()


text:
module MarketData

using CSV
using Dates
using Dates: datetime2unix
using HTTP
using Random
using Reexport
using TimeSeries
using JSON3

@reexport using TimeSeries

export AAPL, BA, CAT, DELL, EBAY, F, GE, TX,
       cl, op, ohlc, ohlcv, datetime1, datetime2, mdata,
       o, h, l, c, v
# downloads.jl
export AbstractQueryOpt,
       YahooOpt,
       yahoo, fred, ons

###### include ##################

include("utils.jl")
include("const.jl")
include("helpers.jl")
include("downloads.jl")

end


text:
using Test
using Random
using LogicCircuits

x1, x2, x3 = bdd_var(1), bdd_var(2), bdd_var(3)
X = Bdd[x1, x2, x3]

@testset "Valuations" begin
  Random.seed!(101)
  for i ∈ 1:50
    Sc = Random.randsubseq(collect(1:10), 0.5)
    n = length(Sc)
    local p::Int64 = 0
    for (j, x) ∈ enumerate(valuations(Sc))
      @test Set{Int}(keys(x)) == Set{Int}(Sc)
      z = count_zeros(p) - 64 + n
      @test count(y -> y.second, x) == n - z
      @test count(y -> !y.second, x) == z
      p += 1
    end
  end
end

@testset "Conjunctions" begin
  y1, y2, y3 = bdd_var(1), bdd_var(2), bdd_var(3)
  Sc = Int[1, 2, 3]
  E = Bdd[¬y1 ∧ ¬y2 ∧ ¬y3, y1 ∧ ¬y2 ∧ ¬y3, ¬y1 ∧ y2 ∧ ¬y3, y1 ∧ y2 ∧ ¬y3, ¬y1 ∧ ¬y2 ∧ y3,
              y1 ∧ ¬y2 ∧ y3, ¬y1 ∧ y2 ∧ y3, y1 ∧ y2 ∧ y3]
  for (i, α) ∈ enumerate(conjunctions(Sc))
    @test α == E[i]
    for x ∈ valuations(Sc) @test α|x == E[i]|x end
  end
end

@testset "Convals" begin
  y1, y2, y3 = bdd_var(1), bdd_var(2), bdd_var(3)
  Sc = Int[1, 2, 3]
  E = Bdd[¬y1 ∧ ¬y2 ∧ ¬y3, y1 ∧ ¬y2 ∧ ¬y3, ¬y1 ∧ y2 ∧ ¬y3, y1 ∧ y2 ∧ ¬y3, ¬y1 ∧ ¬y2 ∧ y3,
              y1 ∧ ¬y2 ∧ y3, ¬y1 ∧ y2 ∧ y3, y1 ∧ y2 ∧ y3]
  for (i, P) ∈ enumerate(convals(Sc))
    α, x = P
    @test α == E[i]
    for y ∈ valuations(Sc) @test α|y == E[i]|y end
    n = length(Sc)
    z = count_zeros(i-1) - 64 + n
    @test Set{Int}(keys(x)) == Set{Int}(Sc)
    @test count(y -> y.second, x) == n - z
    @test count(y -> !y.second, x) == z
  end
end

@testset "Terminal ⊤" begin
  @test is_⊤(⊤)
  @test !is_⊥(⊤)
  @test is_term(⊤)
  @test !is_var(⊤)
  @test !isdefined(⊤, :high)
  @test !isdefined(⊤, :low)
  @test ⊤.value
  @test ⊤.index == -1
  @test terminal(true) == ⊤
  @test terminal(false) == ¬⊤
  @test !is_lit(⊤)
  @test sign(⊤) == 0
  @test to_int(⊤) == 0
  @test is_atom(⊤)
end

@testset "Terminal ⊥" begin
  @test is_⊥(⊥)
  @test !is_⊤(⊥)
  @test is_term(⊥)
  @test !is_var(⊥)
  @test !isdefined(⊥, :high)
  @test !isdefined(⊥, :low)
  @test !⊥.value
  @test ⊥.index == -1
  @test terminal(false) == ⊥
  @test terminal(true) == ¬⊥
  @test !is_lit(⊥)
  @test sign(⊥) == 0
  @test to_int(⊥) == 0
  @test is_atom(⊥)
end

@testset "Variable" begin
  for (i, v) ∈ enumerate(X)
    @test is_var(v)
    @test is_lit(v)
    @test !is_term(v)
    @test !is_⊤(v)
    @test !is_⊥(v)
    @test v.index == i
    @test v.high.value
    @test !v.low.value
    @test v == i
    @test i == v
    @test v != i+1
    @test i+1 != v
    @test ¬v == ¬i
    @test v == ¬¬i
    @test ¬v == ¬¬¬i
    @test sign(v) == 1
    @test sign(¬v) == -1
    @test to_int(v) == i
    @test to_int(¬v) == -i
    @test is_atom(v)
    @test is_atom(¬v)
  end
end

@testset "Forced construction" begin
  for i ∈ 3:100
    @test Bdd(1, i, ⊥, ⊤).id == i
  end
end

@testset "Reduce" begin
  @test reduce!(⊤) == ⊤
  @test reduce!(⊥) == ⊥

  for v ∈ X @test reduce!(v) == v end

  a = Bdd(3, ⊥, ⊤)
  b = Bdd(3, ⊥, ⊤)
  c = Bdd(2, b, a)
  d = Bdd(2, ⊥, b)
  e = Bdd(1, d, c)
  R = reduce!(e)
  E = Union{Int, Bool}[1, 2, false, 3, true]
  i = 1
  foreach(function(α::Bdd)
            v = E[i]
            if is_term(α) @test v == α.value
            else @test v == α.index end
            i += 1
          end, R)
end

@testset "Restrict" begin
  Y = Dict{Int, Bool}(1 => true, 2 => false, 3 => true)

  @test ⊤|Y == ⊤
  @test ⊥|Y == ⊥

  for (i, x) ∈ enumerate(X)
    if Y[i] @test x|Y == ⊤
    else @test x|Y == ⊥ end
  end

  mkeg(P::Vararg{Pair{Int, Bool}})::Vector{Dict{Int, Bool}} = collect(Dict{Int, Bool}.(P))
  mkeg!(P::Vararg{Pair{Int, Bool}})::Dict{Int, Bool} = Dict{Int, Bool}(P)

  f1 = x1 ∧ ¬x2
  Y = mkeg(1 => false, 2 => false)
  @test f1|Y[1] == ⊥
  @test f1|Y[2] == x1
  @test f1|merge(Y...) == ⊥

  Y = mkeg(1 => false, 2 => true)
  @test f1|Y[1] == ⊥
  @test f1|Y[2] == ⊥
  @test f1|merge(Y...) == ⊥

  Y = mkeg(1 => true, 2 => false)
  @test f1|Y[1] == ¬x2
  @test f1|Y[2] == x1
  @test f1|merge(Y...) == ⊤

  Y = mkeg(1 => true, 2 => true)
  @test f1|Y[1] == ¬x2
  @test f1|Y[2] == ⊥
  @test f1|merge(Y...) == ⊥

  f2 = x2 ∧ ¬x3
  Y = mkeg!(1 => false, 4 => false)
  @test f2|Y == f2
  Y = mkeg!(1 => false, 4 => true)
  @test f2|Y == f2
  Y = mkeg!(1 => true, 4 => false)
  @test f2|Y == f2
  Y = mkeg!(1 => true, 4 => true)
  @test f2|Y == f2

  f3 = ¬x2 ∨ x3
  E = Bdd[⊤, ⊤, ⊥, ⊥, ⊤, ⊤, ⊤, ⊤, ⊤, ⊤, ⊥, ⊥, ⊤, ⊤, ⊤, ⊤]
  for (i, y) ∈ enumerate(valuations(1:4)) f3|y == E[i] end

  f4 = x1 ∨ ¬x2
  @test f4|1 == ⊤
  @test f4|-1 == ¬x2
  @test f4|2 == x1
  @test f4|-2 == ⊤
  @test f4|[1, 2] == ⊤
  @test f4|[-1, 2] == ⊥
  @test f4|[1, -2] == ⊤
  @test f4|[-1, -2] == ⊤

  f5 = ¬x1 ∧ x2
  @test f5|1 == ⊥
  @test f5|-1 == x2
  @test f5|2 == ¬x1
  @test f5|-2 == ⊥
  @test f5|[1, 2] == ⊥
  @test f5|[-1, 2] == ⊤
  @test f5|[1, -2] == ⊥
  @test f5|[-1, -2] == ⊥

  Φ = [(1 ∧ ¬2) ∨ (2 ∧ 3), (1 ∨ ¬2) ∧ (2 ∨ 3), 1 ∧ ¬2 ∨ ((3 ∧ 2) ∨ (5 ∧ ¬4)),
       (1 ∧ 3 ∧ 2) ∨ (1 ∧ ¬3 ∧ 2) ∨ (1 ∧ ¬ 2 ∧ 3)]
  for ϕ ∈ Φ
    S = scope(ϕ)
    for i ∈ 1:2^length(S)
      shuffle!(S)
      α, β, γ, δ = copy(ϕ), copy(ϕ), copy(ϕ), copy(ϕ)
      for x ∈ S
        α, β = α|x, β|-x
        @test α == γ|x
        @test β == δ|-x
        γ, δ = copy(α), copy(β)
      end
    end
  end

  for i ∈ 1:1000
    for ϕ ∈ Φ
      X = map(x -> rand(Bool) ? x : -x, randsubseq(scope(ϕ), 0.5))
      α, β = copy(ϕ), copy(ϕ)
      for x ∈ X
        α = α|x
        @test α == β|x
        β = copy(α)
      end
    end
  end
end

@testset "Evaluation" begin
  for (i, x) ∈ enumerate(X)
    @test x(i)
    @test !x(-i)
    @test !(¬x)(i)
    @test (¬x)(-i)
  end

  n = 3
  Y = collect(valuations(1:n))
  dict2vec(C::Dict{Int, Bool})::Vector{Int} = map(p -> p.second ? p.first : -p.first, collect(C))
  dict2bitvec(C::Dict{Int, Bool})::BitVector = (V = BitVector(undef, n); foreach(p -> V[p.first] = p.second, collect(C)); V)
  dict2sub1int(C::Dict{Int, Bool})::SubArray{Int} = view(dict2vec(C), :)
  dict2sub1bool(C::Dict{Int, Bool})::SubArray{Bool} = view(dict2bitvec(C), :)
  dict2sub2int(C::Dict{Int, Bool})::SubArray{Int} = view(dict2vec(C), :, :)
  dict2sub2bool(C::Dict{Int, Bool})::SubArray{Bool} = view(dict2bitvec(C), :, :)

  ⋀ = and(collect(1:n))
  for y ∈ Y
    if all(values(y))
      @test ⋀(y)
      @test ⋀(dict2vec(y))
      @test ⋀(dict2bitvec(y))
      @test ⋀(dict2sub1int(y))
      @test ⋀(dict2sub1bool(y))
      @test ⋀(dict2sub2int(y))
      @test ⋀(dict2sub2bool(y))
    else
      @test !⋀(y)
      @test !⋀(dict2vec(y))
      @test !⋀(dict2bitvec(y))
      @test !⋀(dict2sub1int(y))
      @test !⋀(dict2sub1bool(y))
      @test !⋀(dict2sub2int(y))
      @test !⋀(dict2sub2bool(y))
    end
  end

  ⋁ = or(collect(1:n))
  for y ∈ Y
    if all(p -> !p, values(y))
      @test !⋁(y)
      @test !⋁(dict2vec(y))
      @test !⋁(dict2bitvec(y))
      @test !⋁(dict2sub1int(y))
      @test !⋁(dict2sub1bool(y))
      @test !⋁(dict2sub2int(y))
      @test !⋁(dict2sub2bool(y))
    else
      @test ⋁(y)
      @test ⋁(dict2vec(y))
      @test ⋁(dict2bitvec(y))
      @test ⋁(dict2sub1int(y))
      @test ⋁(dict2sub1bool(y))
      @test ⋁(dict2sub2int(y))
      @test ⋁(dict2sub2bool(y))
    end
  end
end

@testset "Negate" begin
  @test ¬false == !false
  @test ¬true == !true

  @test ⊤ == ¬⊥
  @test ⊥ == ¬⊤

  X′ = Bdd[¬x for x ∈ X]
  for (i, x) ∈ enumerate(X′)
    @test !is_term(x)
    @test x.index == i
    @test !x.high.value
    @test x.low.value
  end

  a = Bdd(3, ⊥, ⊤)
  b = Bdd(2, ⊥, a)
  c = Bdd(1, b, a)
  d = ¬c
  E = Union{Int, Bool}[1, 2, true, 3, false, 3]
  i = 1
  foreach(function(α::Bdd)
            v = E[i]
            if is_term(α) @test v == α.value
            else @test v == α.index end
            i += 1
          end, d)
end

@testset "Conjunction" begin
  @test false ∧ false == false & false
  @test false ∧ true == false & true
  @test true ∧ false == true & false
  @test true ∧ true == true & true

  # Idempotency.
  @test ⊤ ∧ ⊤ == ⊤
  @test ⊥ ∧ ⊥ == ⊥
  @test x1 ∧ x1 == x1
  @test x2 ∧ x2 == x2
  @test x3 ∧ x3 == x3
  a = Bdd(3, ⊥, ⊤)
  b = Bdd(2, ⊥, a)
  c = Bdd(1, b, a)
  @test c ∧ c == c

  # Commutative.
  @test ⊥ ∧ ⊤ == ⊤ ∧ ⊥
  @test x1 ∧ x2 == x2 ∧ x1
  @test (x1 ∧ x2) ∧ x3 == x3 ∧ (x1 ∧ x2)
  @test c ∧ (x2 ∧ ¬x3) == (x2 ∧ ¬x3) ∧ c

  # Associative.
  @test x1 ∧ (x2 ∧ x3) == (x1 ∧ x2) ∧ x3
  @test ¬x1 ∧ (x2 ∧ x3) == (¬x1 ∧ x2) ∧ x3
  @test x1 ∧ (¬x2 ∧ x3) == (x1 ∧ ¬x2) ∧ x3
  @test x1 ∧ (x2 ∧ ¬x3) == (x1 ∧ x2) ∧ ¬x3

  # Neutral element.
  @test x1 ∧ ⊤ == x1
  @test x2 ∧ ⊤ == x2
  @test x3 ∧ ⊤ == x3
  @test ¬x1 ∧ ⊤ == ¬x1
  @test ¬x2 ∧ ⊤ == ¬x2
  @test ¬x3 ∧ ⊤ == ¬x3
  @test (x1 ∧ x2 ∧ x3) ∧ ⊤ == x1 ∧ x2 ∧ x3
  @test (¬x1 ∧ ¬x2 ∧ ¬x3) ∧ ⊤ == ¬x1 ∧ ¬x2 ∧ ¬x3
  @test c ∧ ⊤ == c

  # Opposite element.
  @test x1 ∧ ¬x1 == ⊥
  @test ¬x2 ∧ x2 == ⊥
  @test ¬x3 ∧ x3 == ⊥
  @test c ∧ ¬c == ⊥
  @test ¬(x1 ∧ ¬x3 ∧ x2) ∧ (x1 ∧ ¬x3 ∧ x2) == ⊥
  @test ¬x1 ∧ ¬¬x1 == ⊥
  @test ¬(¬x2) ∧ ¬x2 == ⊥
  @test ¬(¬x3) ∧ ¬x3 == ⊥
  @test ¬c ∧ ¬¬c == ⊥
  @test ¬¬(x1 ∧ ¬x3 ∧ x2) ∧ ¬(x1 ∧ ¬x3 ∧ x2) == ⊥

  # Distributive conjunction over disjunction.
  @test x1 ∧ (x2 ∨ x3) == (x1 ∧ x2) ∨ (x1 ∧ x3)
  @test c ∧ (x2 ∨ x3) == (c ∧ x2) ∨ (c ∧ x3)
  @test x1 ∧ (c ∨ x3) == (x1 ∧ c) ∨ (x1 ∧ x3)
  @test x1 ∧ (x2 ∨ c) == (x1 ∧ x2) ∨ (x1 ∧ c)

  @test x1 ∧ (x2 ∧ x3) == and(x1, and(x2, x3))
  @test (x2 ∧ x3) ∧ x1 == and(and(x2, x3), x1)
  @test x3 ∧ x2 ∧ x1 == and(and(x3, x2), x1)

  @test 1 ∧ ¬2 == x1 ∧ ¬x2
  @test ¬(3 ∧ 2) == ¬(x3 ∧ x2)
  @test 1 ∧ ¬2 ∧ 3 == x1 ∧ ¬x2 ∧ x3
  @test ¬(1 ∧ (2 ∧ ¬3)) == ¬(x1 ∧ (x2 ∧ ¬x3))
  @test 1 ∧ (2 ∨ 3) ∧ (¬3 ∨ ¬1) == x1 ∧ (x2 ∨ x3) ∧ (¬x3 ∨ ¬x1)

  @test 1 ∧ (2 ∧ 3) == and(1, and(2, 3))
  @test (2 ∧ 3) ∧ 1 == and(and(2, 3), 1)
  @test 3 ∧ 2 ∧ 1 == and(and(3, 2), 1)

  @test !is_atom(x1 ∧ x3)
  @test !is_atom(x1 ∧ x2 ∨ x3)

  @test and(1, 2, -3) == 1 ∧ 2 ∧ ¬3
  @test and([1, 2, -3]) == 1 ∧ 2 ∧ ¬3
  @test and(x1, x2, ¬x3) == 1 ∧ 2 ∧ ¬3
  @test and([x1, x2, ¬x3]) == 1 ∧ 2 ∧ ¬3
  @test and(-1, 2, -3) == ¬1 ∧ 2 ∧ ¬3
  @test and([-1, 2, -3]) == ¬1 ∧ 2 ∧ ¬3
  @test and(¬x1, x2, ¬x3) == ¬1 ∧ 2 ∧ ¬3
  @test and([¬x1, x2, ¬x3]) == ¬1 ∧ 2 ∧ ¬3
  @test and(1, -2, -3) == 1 ∧ ¬2 ∧ ¬3
  @test and([1, -2, -3]) == 1 ∧ ¬2 ∧ ¬3
  @test and(x1, ¬x2, ¬x3) == 1 ∧ ¬2 ∧ ¬3
  @test and([x1, ¬x2, ¬x3]) == 1 ∧ ¬2 ∧ ¬3
  @test and(1 ∨ 2, 3 ∨ ¬2, ¬1 ∨ 3) == (1 ∨ 2) ∧ (3 ∨ ¬2) ∧ (¬1 ∨ 3)
  @test and(1 ∨ 2, 3, ¬1 ∨ 3, 1) == (1 ∨ 2) ∧ 3 ∧ (¬1 ∨ 3) ∧ 1
end

@testset "Disjunction" begin
  @test false ∨ false == false | false
  @test false ∨ true == false | true
  @test true ∨ false == true | false
  @test true ∨ true == true | true

  # Idempotency.
  @test ⊤ ∨ ⊤ == ⊤
  @test ⊥ ∨ ⊥ == ⊥
  @test x1 ∨ x1 == x1
  @test x2 ∨ x2 == x2
  @test x3 ∨ x3 == x3
  a = Bdd(3, ⊥, ⊤)
  b = Bdd(2, ⊥, a)
  c = Bdd(1, b, a)
  @test c ∨ c == c

  # Commutative.
  @test ⊥ ∨ ⊤ == ⊤ ∨ ⊥
  @test x1 ∨ x2 == x2 ∨ x1
  @test (x1 ∨ x2) ∨ x3 == x3 ∨ (x1 ∨ x2)
  @test c ∨ (x2 ∨ ¬x3) == (x2 ∨ ¬x3) ∨ c

  # Associative.
  @test x1 ∨ (x2 ∨ x3) == (x1 ∨ x2) ∨ x3
  @test ¬x1 ∨ (x2 ∨ x3) == (¬x1 ∨ x2) ∨ x3
  @test x1 ∨ (¬x2 ∨ x3) == (x1 ∨ ¬x2) ∨ x3
  @test x1 ∨ (x2 ∨ ¬x3) == (x1 ∨ x2) ∨ ¬x3

  # Neutral element.
  @test x1 ∨ ⊥ == x1
  @test x2 ∨ ⊥ == x2
  @test x3 ∨ ⊥ == x3
  @test ¬x1 ∨ ⊥ == ¬x1
  @test ¬x2 ∨ ⊥ == ¬x2
  @test ¬x3 ∨ ⊥ == ¬x3
  @test (x1 ∨ x2 ∨ x3) ∨ ⊥ == x1 ∨ x2 ∨ x3
  @test (¬x1 ∨ ¬x2 ∨ ¬x3) ∨ ⊥ == ¬x1 ∨ ¬x2 ∨ ¬x3
  @test c ∨ ⊥ == c

  # Opposite element.
  @test x1 ∨ ¬x1 == ⊤
  @test ¬x2 ∨ x2 == ⊤
  @test ¬x3 ∨ x3 == ⊤
  @test c ∨ ¬c == ⊤
  @test ¬(x1 ∨ ¬x3 ∨ x2) ∨ (x1 ∨ ¬x3 ∨ x2) == ⊤
  @test ¬x1 ∨ ¬¬x1 == ⊤
  @test ¬(¬x2) ∨ ¬x2 == ⊤
  @test ¬(¬x3) ∨ ¬x3 == ⊤
  @test ¬c ∨ ¬¬c == ⊤
  @test ¬¬(x1 ∨ ¬x3 ∨ x2) ∨ ¬(x1 ∨ ¬x3 ∨ x2) == ⊤

  # Distributive disjunction over conjunction.
  @test x1 ∨ (x2 ∧ x3) == (x1 ∨ x2) ∧ (x1 ∨ x3)
  @test c ∨ (x2 ∧ x3) == (c ∨ x2) ∧ (c ∨ x3)
  @test x1 ∨ (c ∧ x3) == (x1 ∨ c) ∧ (x1 ∨ x3)
  @test x1 ∨ (x2 ∧ c) == (x1 ∨ x2) ∧ (x1 ∨ c)

  @test x1 ∨ (x2 ∨ x3) == or(x1, or(x2, x3))
  @test (x2 ∨ x3) ∨ x1 == or(or(x2, x3), x1)
  @test x3 ∨ x2 ∨ x1 == or(or(x3, x2), x1)

  @test 1 ∨ ¬2 == x1 ∨ ¬x2
  @test ¬(3 ∨ 2) == ¬(x3 ∨ x2)
  @test 1 ∨ ¬2 ∨ 3 == x1 ∨ ¬x2 ∨ x3
  @test ¬(1 ∨ (2 ∨ ¬3)) == ¬(x1 ∨ (x2 ∨ ¬x3))
  @test 1 ∨ (2 ∧ 3) ∨ (¬3 ∧ ¬1) == x1 ∨ (x2 ∧ x3) ∨ (¬x3 ∧ ¬x1)

  @test 1 ∨ (2 ∨ 3) == or(1, or(2, 3))
  @test (2 ∨ 3) ∨ 1 == or(or(2, 3), 1)
  @test 3 ∨ 2 ∨ 1 == or(or(3, 2), 1)

  @test !is_atom(x1 ∨ ¬x3)
  @test !is_atom(x2 ∨ x3 ∧ x1)

  @test or(1, 2, -3) == 1 ∨ 2 ∨ ¬3
  @test or([1, 2, -3]) == 1 ∨ 2 ∨ ¬3
  @test or(x1, x2, ¬x3) == 1 ∨ 2 ∨ ¬3
  @test or([x1, x2, ¬x3]) == 1 ∨ 2 ∨ ¬3
  @test or(-1, 2, -3) == ¬1 ∨ 2 ∨ ¬3
  @test or([-1, 2, -3]) == ¬1 ∨ 2 ∨ ¬3
  @test or(¬x1, x2, ¬x3) == ¬1 ∨ 2 ∨ ¬3
  @test or([¬x1, x2, ¬x3]) == ¬1 ∨ 2 ∨ ¬3
  @test or(1, -2, -3) == 1 ∨ ¬2 ∨ ¬3
  @test or([1, -2, -3]) == 1 ∨ ¬2 ∨ ¬3
  @test or(x1, ¬x2, ¬x3) == 1 ∨ ¬2 ∨ ¬3
  @test or([x1, ¬x2, ¬x3]) == 1 ∨ ¬2 ∨ ¬3
  @test or(1 ∧ 2, 3 ∧ ¬2, ¬1 ∧ 3) == (1 ∧ 2) ∨ (3 ∧ ¬2) ∨ (¬1 ∧ 3)
  @test or(1 ∧ 2, 3, ¬1 ∧ 3, 1) == (1 ∧ 2) ∨ 3 ∨ (¬1 ∧ 3) ∨ 1
end

@testset "XOR" begin
  # Idempotency.
  @test ⊤ ⊻ ⊤ == ⊥
  @test ⊥ ⊻ ⊥ == ⊥
  @test x1 ⊻ x1 == ⊥
  @test x2 ⊻ x2 == ⊥
  @test x3 ⊻ x3 == ⊥
  a = Bdd(3, ⊥, ⊤)
  b = Bdd(2, ⊥, a)
  c = Bdd(1, b, a)
  @test c ⊻ c == ⊥

  # Commutative.
  @test ⊥ ⊻ ⊤ == ⊤ ⊻ ⊥
  @test x1 ⊻ x2 == x2 ⊻ x1
  @test (x1 ⊻ x2) ⊻ x3 == x3 ⊻ (x1 ⊻ x2)
  @test c ⊻ (x2 ⊻ ¬x3) == (x2 ⊻ ¬x3) ⊻ c

  # Associative.
  @test x1 ⊻ (x2 ⊻ x3) == (x1 ⊻ x2) ⊻ x3
  @test ¬x1 ⊻ (x2 ⊻ x3) == (¬x1 ⊻ x2) ⊻ x3
  @test x1 ⊻ (¬x2 ⊻ x3) == (x1 ⊻ ¬x2) ⊻ x3
  @test x1 ⊻ (x2 ⊻ ¬x3) == (x1 ⊻ x2) ⊻ ¬x3

  # Neutral element.
  @test x1 ⊻ ⊥ == x1
  @test x2 ⊻ ⊥ == x2
  @test x3 ⊻ ⊥ == x3
  @test ¬x1 ⊻ ⊥ == ¬x1
  @test ¬x2 ⊻ ⊥ == ¬x2
  @test ¬x3 ⊻ ⊥ == ¬x3
  @test (x1 ⊻ x2 ⊻ x3) ⊻ ⊥ == x1 ⊻ x2 ⊻ x3
  @test (¬x1 ⊻ ¬x2 ⊻ ¬x3) ⊻ ⊥ == ¬x1 ⊻ ¬x2 ⊻ ¬x3
  @test c ⊻ ⊥ == c

  # Opposite element.
  @test x1 ⊻ ¬x1 == ⊤
  @test ¬x2 ⊻ x2 == ⊤
  @test ¬x3 ⊻ x3 == ⊤
  @test c ⊻ ¬c == ⊤
  @test ¬(x1 ⊻ ¬x3 ⊻ x2) ⊻ (x1 ⊻ ¬x3 ⊻ x2) == ⊤
  @test ¬x1 ⊻ ¬¬x1 == ⊤
  @test ¬(¬x2) ⊻ ¬x2 == ⊤
  @test ¬(¬x3) ⊻ ¬x3 == ⊤
  @test ¬c ⊻ ¬¬c == ⊤
  @test ¬¬(x1 ⊻ ¬x3 ⊻ x2) ⊻ ¬(x1 ⊻ ¬x3 ⊻ x2) == ⊤

  @test bdd_var(1) ⊻ bdd_var(2) == x1 ⊻ x2
  @test bdd_var(1) ⊻ ¬3 == x1 ⊻ ¬x3
  @test ¬(¬2 ⊻ 1) == ¬(¬x2 ⊻ x1)
  @test ¬(1 ⊻ ¬3 ⊻ 2) == ¬(x1 ⊻ ¬x3 ⊻ x2)
  @test ¬(1 ⊻ (¬2 ∨ 3) ⊻ ¬(3 ∧ 2)) ∧ (3 ∧ ¬2) == ¬(x1 ⊻ (¬x2 ∨ x3) ⊻ ¬(x3 ∧ x2)) ∧ (x3 ∧ ¬x2)
end

@testset "Implication" begin
  @test (false → false) == true
  @test (false → true) == true
  @test (true → false) == false
  @test (true → true) == true

  for i ∈ 1:100
    @test (i → (i+1)) == ((¬i) ∨ (i+1))
    @test (bdd_var(i) → bdd_var(i+1)) == ((¬i) ∨ (i+1))
  end

  α = (1 ∨ ¬2) ∧ (¬3 ∨ 4) ∧ (¬2 ∨ ¬4)
  β = (5 ∧ 1) ∨ (¬3 ∧ 2) ∨ (¬1 ∧ 4)
  @test (α → β) == ((¬α) ∨ β)
  @test (α → 2) == ((¬α) ∨ 2)
  @test (-3 → α) == (3 ∨ α)
end

@testset "Equality and inequality" begin
  @test isequal(⊥, ⊥)
  @test isequal(⊤, ⊤)
  @test !isequal(⊥, ⊤)
  @test !isequal(⊤, ⊥)

  @test !isequal(⊥, ⊥) == (⊥ != ⊥)
  @test !isequal(⊤, ⊤) == (⊤ != ⊤)
  @test !isequal(⊥, ⊤) == (⊥ != ⊤)
  @test !isequal(⊤, ⊥) == (⊤ != ⊥)

  ϕ = ¬(x1 ∨ (x2 ∧ x3)) ∨ ((x3 ∧ x1) ∨ x2)
  @test isequal(ϕ, ϕ) == (ϕ == ϕ)
  @test isequal(ϕ, ⊥) == (ϕ == ⊥)
  @test !isequal(ϕ, ϕ) == (ϕ != ϕ)
  @test !isequal(ϕ, ⊥) == (ϕ != ⊥)

  @test !isequal(ϕ, 2) == (ϕ != 2)
  @test !isequal(2, ϕ) == (2 != ϕ)
  @test !isequal(ϕ, ¬2) == (ϕ != ¬2)
  @test isequal(ϕ, ¬(1 ∨ (2 ∧ 3)) ∨ ((3 ∧ 1) ∨ 2))
end

@testset "Iterators" begin
  Φ = Bdd[¬(x1 ∨ (x2 ∧ x3)) ∨ ((x3 ∧ x1) ∨ x2), x1 ∧ x2 ∨ x3, (x1 ∧ x2) ∨ (¬x2 ∧ ¬x3),
              (x1 ∧ x2) ∨ (x2 ∧ ¬x3 ∧ ¬x1)]
  for ϕ ∈ Φ
    U, V, W = Vector{Bdd}(), Vector{Bdd}(), Vector{Bdd}()
    for v ∈ ϕ push!(V, v) end
    foreach(x -> push!(U, x), ϕ)
    W = collect(ϕ)
    @test U == V
    @test V == W
    # Transitive, therefore expect U == W if both conditions satisfy.
    @test U == W
  end
end

@testset "Hash function" begin
  Φ = Bdd[¬(x1 ∨ (x2 ∧ x3)) ∨ ((x3 ∧ x1) ∨ x2), x1 ∧ x2 ∨ x3, (x1 ∧ x2) ∨ (¬x2 ∧ ¬x3),
       (x1 ∧ x2) ∨ (x2 ∧ ¬x3 ∧ ¬x1)]
  test_hash(ϕ::Bdd) = foreach((x -> @test shallowhash(x) == hash((x.id, x.value, x.index))), ϕ)
  test_hash.(Φ)

  @test shallowhash(⊤) != shallowhash(x1)
  @test shallowhash(⊥) != shallowhash(x1)
  @test shallowhash(⊤) != shallowhash(⊥)

  H = hash.(Φ)
  @test allunique(H)

  @test hash(bdd_var(1)) == hash(x1)
  @test hash(bdd_var(2)) == hash(x2)
  @test hash(bdd_var(3)) == hash(x3)
  @test hash(⊤) != hash(x1)
  @test hash(⊥) != hash(x1)
  @test hash(⊤) != hash(⊥)
  @test hash(⊥) == hash(1 ∧ ¬1)
  @test hash(⊤) == hash(1 ∨ ¬1)
  @test hash((1 ∧ 2) ∧ 3) == hash(1 ∧ (2 ∧ 3))
  @test hash((1 ∨ 2) ∨ 3) == hash(1 ∨ (2 ∨ 3))
end

@testset "Shannon decomposition" begin
  function test_formula(ϕ::Bdd, E::Vector{Tuple{Bdd, Bdd, Bdd, Bdd}})
    for (i, e) ∈ enumerate(E)
      u, α, v, β = shannon(ϕ, i)
      @test u == e[1]
      @test α == e[2]
      @test v == e[3]
      @test β == e[4]
    end
  end

  test_formula(x1 ∧ x2 ∨ x3, [(x1, x2 ∨ x3, ¬x1, x3), (x2, x1 ∨ x3, ¬x2, x3), (x3, ⊤, ¬x3, x1 ∧ x2)])
  test_formula((x1 ∧ x2) ∨ (¬x2 ∧ ¬x3), [(x1, x2 ∨ ¬x3, ¬x1, ¬x2 ∧ ¬x3), (x2, x1, ¬x2, ¬x3), (x3, x1 ∧ x2, ¬x3, x1 ∨ ¬x2)])

  ϕ = (x1 ∧ x2) ∨ (x2 ∧ ¬x3 ∧ ¬x1)
  E = [⊥, ⊥, ¬x3, ⊤]
  D = shannon(ϕ, [1, 2])
  for (i, x) ∈ enumerate(conjunctions([1, 2]))
    @test D[i][1] == x
    @test D[i][2] == E[i]
  end

  D = shannon!(ϕ, [1, 2])
  @test length(D) == 2
  @test D[1][1] == ¬x1 ∧ x2
  @test D[1][2] == ¬x3
  @test D[2][1] == x1 ∧ x2
  @test D[2][2] == ⊤
end

@testset "Copy and deep copy" begin
  @test x1 == copy(x1)
  @test x2 == copy(x2)
  @test x3 == copy(x3)
  @test ⊤ == copy(⊤)
  @test ⊥ == copy(⊥)
  @test ¬x1 == ¬copy(x1)
  @test ¬x2 == ¬copy(x2)
  @test ¬x3 == ¬copy(x3)

  Φ = Bdd[¬(x1 ∨ (x2 ∧ x3)) ∨ ((x3 ∧ x1) ∨ x2), x1 ∧ x2 ∨ x3, (x1 ∧ x2) ∨ (¬x2 ∧ ¬x3),
              (x1 ∧ x2) ∨ (x2 ∧ ¬x3 ∧ ¬x1)]
  for ϕ ∈ Φ
    @test ϕ == deepcopy(ϕ)
  end
end

@testset "String representation" begin
  @test string(⊤) == "@ (value=true, id=$(⊤.id))\n"
  @test string(⊥) == "@ (value=false, id=$(⊥.id))\n"
  for x ∈ X
    @test string(x) == "@ (index=$(x.index), id=$(x.id))\n|  - (value=false, id=$(x.low.id))\n|  + (value=true, id=$(x.high.id))\n"
  end

  ϕ = x1 ∧ x2 ∨ x3
  @test string(x1 ∧ x2 ∨ x3) === """@ (index=$(x1.index), id=$(ϕ.id))
|  - (index=$(x3.index), id=$(ϕ.low.id))
|  |  - (value=$(x3.low.value), id=$(ϕ.low.low.id))
|  |  + (value=$(x3.high.value), id=$(ϕ.high.high.id))
|  + (index=$(x2.index), id=$(ϕ.high.id))
|  |  - (index=$(x3.index), id=$(ϕ.low.id))
|  |  |  - (value=$(x3.low.value), id=$(ϕ.high.low.low.id))
|  |  |  + (value=$(x3.high.value), id=$(ϕ.high.low.high.id))
|  |  + (value=$(x2.high.value), id=$(ϕ.high.high.id))
"""
end

@testset "Elimination" begin
  elimtest(α::Bdd, Sc::UnitRange{Int}) = for v ∈ Sc @test eliminate(α, v) == ((α|v) ∨ (α|-v)) end
  elimtest(⊥, 1:2)
  elimtest(⊤, 1:2)
  elimtest(bdd_var(1), 1:2)
  elimtest(¬1, 1:2)
  elimtest(1 ∨ 2 ∧ 3, 1:3)
  elimtest(1 ∧ ¬2 ∨ 3, 1:4)
  elimtest((1 ∧ ¬2) ∨ (¬3 ∧ ¬4), 1:6)
  elimtest((1 ∨ ¬2) ∧ (¬3 ∨ ¬4), 1:6)
  elimtest((1 ∨ 2 ∧ ¬3) ∧ 2 ∨ (¬4 ∧ ¬5 ∨ 6) ∧ (1 ∨ ¬2), 1:8)
  elimtest(2 ∧ (2 ∨ ¬3) ∧ (¬3 ∨ ¬4), 1:4)

  @test eliminate(⊤, 1) == ⊤
  @test eliminate(⊥, 1) == ⊥

  Φ = [(1 ∨ 2) ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4), and(collect(1:5)), or(collect(1:5)),
       (1 ∨ ¬2) ∧ (3 ∨ ¬4) ∧ (3 ∨ 5) ∧ (4 ∨ ¬5), atmost(5, collect(1:10)),
       atleast(5, collect(1:10)), exactly(5, collect(1:10)), (1 ∧ 2) ∨ (3 ∧ 4),
       (1 ∧ ¬2) ∨ (¬3 ∧ 4), (1 ∧ ¬2) ∨ (3 ∧ ¬4) ∨ (3 ∧ 5) ∨ (4 ∧ ¬5)]
  for (i, ϕ) ∈ enumerate(Φ)
    sc = scope(ϕ)
    for x ∈ sc
      f, g = ϕ|x, ϕ|-x
      @test eliminate(ϕ, x) == (f ∨ g)
    end
    α = ϕ
    rem = Int[]
    rem_set = BitSet()
    while !isempty(sc)
      x = pop!(sc)
      f, g = α|x, α|-x
      β = eliminate(α, x)
      @test β == (f ∨ g)
      α = β
      push!(rem, x)
      push!(rem_set, x)
      @test eliminate(ϕ, rem) == α
      @test eliminate(ϕ, rem_set) == α
    end
  end
end

@testset "Marginalization" begin
  margtest(α::Bdd, ⊕, Sc::UnitRange{Int}) = for v ∈ Sc
    @test marginalize(α, v, ⊕) == (α|v) ⊕ (α|-v)
  end
  for o ∈ Function[⊻, ∨, ∧]
    margtest(⊥, o, 1:2)
    margtest(⊤, o, 1:2)
    margtest(bdd_var(1), o, 1:2)
    margtest(¬1, o, 1:2)
    margtest(1 ∨ 2 ∧ 3, o, 1:3)
    margtest(1 ∧ ¬2 ∨ 3, o, 1:4)
    margtest((1 ∧ ¬2) ∨ (¬3 ∧ ¬4), o, 1:6)
    margtest((1 ∨ ¬2) ∧ (¬3 ∨ ¬4), o, 1:6)
    margtest((1 ∨ 2 ∧ ¬3) ∧ 2 ∨ (¬4 ∧ ¬5 ∨ 6) ∧ (1 ∨ ¬2), o, 1:8)
  end
end

@testset "Scope" begin
  F = Tuple{Bdd, Set{Int}}[
    (⊥, Set{Int}()),
    (⊤, Set{Int}()),
    (bdd_var(1), Set{Int}([1])),
    (bdd_var(2), Set{Int}([2])),
    (¬3, Set{Int}([3])),
    (¬4, Set{Int}([4])),
    (1 ∧ 2, Set{Int}([1, 2])),
    (2 ∧ 4, Set{Int}([2, 4])),
    (2 ∨ 3, Set{Int}([2, 3])),
    (3 ∨ 5, Set{Int}([3, 5])),
    (1 ∧ 2 ∨ ¬3, Set{Int}(collect(1:3))),
    (2 ∨ ¬1 ∧ 3, Set{Int}(collect(1:3))),
    ((1 ∧ ¬2) ∨ (¬3 ∧ ¬4), Set{Int}(collect(1:4))),
    ((1 ∨ ¬2) ∧ (¬3 ∨ ¬4), Set{Int}(collect(1:4))),
    ((1 ∧ 2 ∨ ¬3) ∧ (¬2 ∧ 1 ∨ ¬4) ∧ (¬3 ∨ 5 ∧ 6), Set{Int}(collect(1:6)))
  ]
  for (α, Sc) ∈ F
    @test scopeset(α) == Sc
    @test Set{Int}(scope(α)) == Sc
  end
end

@testset "Mentions" begin
  Sc = collect(1:10)
  T = [collect(1:4), [1, 3, 5, 7, 9], [2, 4, 6, 8], [1, 2, 3, 6, 7, 8], [1, 5, 8], collect(1:10)]
  Φ = [(1 ∧ 2) ∨ (¬2 ∧ 3) ∨ (¬3 ∧ 4), (1 ∧ 3 ∧ 5) ∨ (7 ∧ 9), 2 ∨ 4 ∨ 6 ∨ 8, 1 ∧ 2 ∨ 3 ∧ 6 ∨ 7 ∧ 8,
       (1 ∧ 5) ∨ (¬1 ∧ ¬5) ∨ (5 ∧ 8) ∨ (¬8 ∧ 1), and(collect(1:10))]
  C = [Int32, Int64]
  C_v = [Vector{Int32}, Vector{Int64}]
  for (t, ϕ) ∈ zip(T, Φ)
    for x ∈ Sc
      for cast ∈ C
        y = cast(x)
        @test mentions(ϕ, y) == (y ∈ t)
        @test !mentions(ϕ, y) == (y ∉ t)
        @test (y ∈ ϕ) == (y ∈ t)
        @test (y ∉ ϕ) == (y ∉ t)
      end
    end
    for cast ∈ C_v
      u = convert(cast, t)
      @test mentions(ϕ, u)
      @test !mentions(ϕ, setdiff(Sc, u))
      @test mentions(ϕ, u[begin:3])
      @test mentions(ϕ, u[2:end])
      @test u ∈ ϕ
      @test setdiff(Sc, u) ∉ ϕ
      @test u[begin:3] ∈ ϕ
      @test u[2:end] ∈ ϕ
    end
  end
end

@testset "Literal vector" begin
  Φ = [1∧2∧3∧4∧5, 1∧¬2∧3∧4∧¬5, 4∧¬3∧7∧¬1, 5∧4∧¬2∧¬1, ¬1∧¬5∧¬4∧¬2]
  E = [trues(5), [1,0,1,1,0], [0,0,1,1], [0,0,1,1], falses(4)]
  for (ϕ, e) ∈ zip(Φ, E)
    X, S = lit_vec(ϕ)
    @test X == e
    @test S == sort(scope(ϕ))
  end
end

@testset "Literal conversion" begin
  m = 50
  E_p, E_n = 1:m, -m:-1
  p = bdd_var.(E_p)
  n = bdd_var.(E_n)
  for i ∈ 1:m
    @test to_lit(p[i]) == E_p[i]
    @test to_lit(n[i]) == E_n[i]
  end
end

@testset "At least" begin
  B = [3, 5, 10, 13, 25, 37, 48, 59, 72, 81, 97]
  L = [collect(1:100), collect(1:200), -collect(1:100), -collect(1:200)]
  h_p, h_n = x -> count(values(x)), x -> length(x)-count(values(x))
  how = [h_p, h_p, h_n, h_n]
  for (l, h) ∈ zip(L, how)
    for b ∈ B
      α = atleast(b, l)
      for x ∈ valuations(l)
        if h(x) >= b
          @test α(x) == true
        else
          @test α(x) == false
        end
      end
    end
  end
end

@testset "At most" begin
  B = [3, 5, 10, 13, 25, 37, 48, 59, 72, 81, 97]
  L = [collect(1:100), collect(1:200), -collect(1:100), -collect(1:200)]
  h_p, h_n = x -> count(values(x)), x -> length(x)-count(values(x))
  how = [h_p, h_p, h_n, h_n]
  for (l, h) ∈ zip(L, how)
    for b ∈ B
      α = atmost(b, l)
      for x ∈ valuations(l)
        if h(x) <= b
          @test α(x) == true
        else
          @test α(x) == false
        end
      end
    end
  end
end

@testset "Exactly" begin
  B = [3, 5, 10, 13, 25, 37, 48, 59, 72, 81, 97]
  L = [collect(1:100), collect(1:200), -collect(1:100), -collect(1:200)]
  h_p, h_n = x -> count(values(x)), x -> length(x)-count(values(x))
  how = [h_p, h_p, h_n, h_n]
  for (l, h) ∈ zip(L, how)
    for b ∈ B
      α = atmost(b, l) ∧ atleast(b, l)
      β = exactly(b, l)
      for x ∈ valuations(l)
        if h(x) == b
          @test α(x) == true
          @test β(x) == true
        else
          @test α(x) == false
          @test β(x) == false
        end
      end
    end
  end
end

@testset "Size" begin
  Φ = [⊤, ⊥, bdd_var(1), bdd_var(-1), 1 ∧ 2, 1 ∨ ¬2, 1 ∧ 2 ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4),
       (1 ∧ 2) ∨ (¬2 ∧ 3) ∨ (¬3 ∧ 4), (1 ∧ 3 ∧ 5) ∨ (7 ∧ 9), 2 ∨ 4 ∨ 6 ∨ 8, 1 ∧ 2 ∨ 3 ∧ 6 ∨ 7 ∧ 8,
       (1 ∧ 5) ∨ (¬1 ∧ ¬5) ∨ (5 ∧ 8) ∨ (¬8 ∧ 1), and(collect(1:10))]
  for ϕ ∈ Φ
    @test size(ϕ) == length(collect(ϕ))
  end
end

@testset "Valuations as a Matrix" begin
  Random.seed!(101)
  for i ∈ 1:50
    Sc = Random.randsubseq(collect(1:10), 0.5)
    n = length(Sc)
    M = all_valuations(Sc)
    for (j, X) ∈ enumerate(valuations(Sc))
      Y = BitVector(undef, n)
      K = sort(collect(keys(X)))
      for (l, k) ∈ enumerate(K) Y[l] = X[k] end
      @test M[j,:] == Y
    end
  end
end

@testset "Value of literal in conjunction" begin
  Φ = [1∧2∧3∧4∧5, 1∧¬2∧3∧4∧¬5, 4∧¬3∧7∧¬1, 5∧4∧¬2∧¬1, ¬1∧¬5∧¬4∧¬2]
  C = [UInt16, UInt32, UInt64, Int32, Int64]
  O = [11, 20, 30, 40, 50, 60]
  for ϕ ∈ Φ
    X, S = lit_vec(ϕ)
    for (i, x) ∈ enumerate(S)
      for c ∈ C
        y = c(x)
        @test lit_val(ϕ, y) == !is_⊥(ϕ|x)
        for o ∈ O @test lit_val(ϕ, o) == false end
        if c <: Signed
          @test lit_val(ϕ, -y) == is_⊥(ϕ|x)
          for o ∈ O @test lit_val(ϕ, -o) == false end
        end
      end
    end
  end
end

@testset "Pretty print conjunctions" begin
  Φ = [1∧2∧3∧4∧5, 1∧¬2∧3∧4∧¬5, 4∧¬3∧7∧¬1, 5∧4∧¬2∧¬1, ¬1∧¬5∧¬4∧¬2]
  E = ["1 ∧ 2 ∧ 3 ∧ 4 ∧ 5", "1 ∧ ¬2 ∧ 3 ∧ 4 ∧ ¬5", "¬1 ∧ ¬3 ∧ 4 ∧ 7", "¬1 ∧ ¬2 ∧ 4 ∧ 5",
       "¬1 ∧ ¬2 ∧ ¬4 ∧ ¬5"]
  for (i, ϕ) ∈ enumerate(Φ)
    @test E[i] == print_conjunction(ϕ; out = false)
  end
end

@testset "Parallelism" begin
  n = 8
  m = 2^n
  Φ = Vector{Bdd}(undef, m)
  M = all_valuations(1:n)
  E = and.(map.(x -> x[2] ? x[1] : -x[1], enumerate.(eachrow(M))))
  Threads.@threads for i ∈ 1:m
    Φ[i] = and(map(x -> x[2] ? x[1] : -x[1], enumerate(M[i,:])))
  end
  for i ∈ 1:m @test Φ[i] == E[i] end
  Ψ = Vector{Bdd}(undef, m)
  F = Vector{Bdd}(undef, m)
  Threads.@threads for i ∈ 1:m
    α = (2 ∧ ¬10 ∧ ¬11 ∧ ¬3)
    Ψ[i] = Φ[i] ∨ α
    F[i] = E[i] ∨ α
  end
  for i ∈ 1:m @test Ψ[i] == F[i] end
  Φ = [(1 ∨ 2) ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4), and(collect(1:5)), or(collect(1:5)),
       (1 ∨ ¬2) ∧ (3 ∨ ¬4) ∧ (3 ∨ 5) ∧ (4 ∨ ¬5), atmost(5, collect(1:10)),
       atleast(5, collect(1:10)), exactly(5, collect(1:10)), (1 ∧ 2) ∨ (3 ∧ 4),
       (1 ∧ ¬2) ∨ (¬3 ∧ 4), (1 ∧ ¬2) ∨ (3 ∧ ¬4) ∨ (3 ∧ 5) ∨ (4 ∧ ¬5)]
  for i ∈ 1:length(Φ)
    ϕ = Φ[i]
    S = scope(ϕ)
    ψ = [forget(ϕ, S[1:k]) for k ∈ 1:length(S)]
    α = Vector{Bdd}(undef, length(S))
    Threads.@threads for j ∈ 1:length(S)
      α[j] = forget(ϕ, S[1:j])
    end
    for j ∈ 1:length(S) @test α[j] == ψ[j] end
  end
end

@testset "Post-order" begin
  Φ = [1∧2∧3∧4∧5, 1∧¬2∧3∧4∧¬5, 4∧¬3∧7∧¬1, 5∧4∧¬2∧¬1, ¬1∧¬5∧¬4∧¬2, ⊤, ⊥, bdd_var(1), bdd_var(-1),
       1 ∧ 2, 1 ∨ ¬2, 1 ∧ 2 ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4),
       (1 ∧ 2) ∨ (¬2 ∧ 3) ∨ (¬3 ∧ 4), (1 ∧ 3 ∧ 5) ∨ (7 ∧ 9), 2 ∨ 4 ∨ 6 ∨ 8, 1 ∧ 2 ∨ 3 ∧ 6 ∨ 7 ∧ 8,
       (1 ∧ 5) ∨ (¬1 ∧ ¬5) ∨ (5 ∧ 8) ∨ (¬8 ∧ 1), and(collect(1:10))]
  for ϕ ∈ Φ
    T = postorder(ϕ)
    if is_term(ϕ)
      @test length(T) == 1
      @test first(T) == ϕ
      continue
    end
    V = Set{UInt64}()
    @test is_⊥(T[1]) || is_⊤(T[1])
    @test is_⊥(T[2]) || is_⊤(T[2])
    @test T[1] != T[2]
    push!(V, objectid(T[1]))
    push!(V, objectid(T[2]))
    foreach(function(α::Bdd)
              o = objectid(α)
              @test o ∉ V
              @test objectid(α.low) ∈ V
              @test objectid(α.high) ∈ V
              push!(V, o)
            end, T[3:end])
  end
end

@testset "Load and save CNF" begin
  Φ = [(1 ∨ 2) ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4), and(collect(1:5)), or(collect(1:5)),
       (1 ∨ ¬2) ∧ (3 ∨ ¬4) ∧ (3 ∨ 5) ∧ (4 ∨ ¬5), atmost(5, collect(1:10)),
       atleast(5, collect(1:10)), exactly(5, collect(1:10))]
  for (i, ϕ) ∈ enumerate(Φ)
    f = tempname() * ".cnf"
    save_bdd(ϕ, f)
    α, c, n_c = ⊤, 0, -1
    open(f, "r") do input
      for (j, line) ∈ enumerate(eachline(input))
        l = lstrip(line)
        if l[1] == 'c' continue end
        if l[1] == 'p'
          n, m = map(x -> parse(Int, x), split(l[6:end]))
          @test n == length(scopeset(ϕ))
          n_c = m
          continue
        end
        c += 1
        α = α ∧ reduce(∨, BitSet(map(x -> parse(Int, x), split(l[begin:end-1]))))
      end
    end
    @test α == ϕ
    @test n_c > 0
    @test c == n_c
    β = load(Bdd, f)
    @test β == ϕ
  end
end

@testset "Load and save DNF" begin
  Φ = [(1 ∧ 2) ∨ (3 ∧ 4), (1 ∧ ¬2) ∨ (¬3 ∧ 4), and(collect(1:5)), or(collect(1:5)),
       (1 ∧ ¬2) ∨ (3 ∧ ¬4) ∨ (3 ∧ 5) ∨ (4 ∧ ¬5), atmost(5, collect(1:10)),
       atleast(5, collect(1:10)), exactly(5, collect(1:10))]
  for (i, ϕ) ∈ enumerate(Φ)
    f = tempname() * ".dnf"
    save_bdd(ϕ, f)
    α, c, n_c = ⊥, 0, -1
    open(f, "r") do input
      for (j, line) ∈ enumerate(eachline(input))
        l = lstrip(line)
        if l[1] == 'c' continue end
        if l[1] == 'p'
          n, m = map(x -> parse(Int, x), split(l[6:end]))
          @test n == length(scopeset(ϕ))
          n_c = m
          continue
        end
        c += 1
        α = α ∨ reduce(∧, BitSet(map(x -> parse(Int, x), split(l[begin:end-1]))))
      end
    end
    @test α == ϕ
    @test n_c > 0
    @test c == n_c
    β = load(Bdd, f)
    @test β == ϕ
  end
end

@testset "Save and load BDD" begin
  Φ = [(1 ∨ 2) ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4), and(collect(1:5)), or(collect(1:5)),
       (1 ∨ ¬2) ∧ (3 ∨ ¬4) ∧ (3 ∨ 5) ∧ (4 ∨ ¬5), atmost(5, collect(1:10)),
       atleast(5, collect(1:10)), exactly(5, collect(1:10)), (1 ∧ 2) ∨ (3 ∧ 4),
       (1 ∧ ¬2) ∨ (¬3 ∧ 4), (1 ∧ ¬2) ∨ (3 ∧ ¬4) ∨ (3 ∧ 5) ∨ (4 ∧ ¬5)]
  for ϕ ∈ Φ
    f = tempname() * ".bdd"
    save_bdd(ϕ, f)
    α = load(Bdd, f)
    @test ϕ == α
  end
end

@testset "Forget" begin
  @test forget(⊤, 1) == ⊤
  @test forget(⊥, 1) == ⊥

  Φ = [(1 ∨ 2) ∧ (3 ∨ 4), (1 ∨ ¬2) ∧ (¬3 ∨ 4), and(collect(1:5)), or(collect(1:5)),
       (1 ∨ ¬2) ∧ (3 ∨ ¬4) ∧ (3 ∨ 5) ∧ (4 ∨ ¬5), atmost(5, collect(1:10)),
       atleast(5, collect(1:10)), exactly(5, collect(1:10)), (1 ∧ 2) ∨ (3 ∧ 4),
       (1 ∧ ¬2) ∨ (¬3 ∧ 4), (1 ∧ ¬2) ∨ (3 ∧ ¬4) ∨ (3 ∧ 5) ∨ (4 ∧ ¬5)]
  for (i, ϕ) ∈ enumerate(Φ)
    sc = scope(ϕ)
    for x ∈ sc
      f, g = ϕ|x, ϕ|-x
      @test forget(ϕ, x) == (f ∨ g)
    end
    α = ϕ
    while !isempty(sc)
      x = pop!(sc)
      f, g = α|x, α|-x
      β = forget(α, x)
      @test β == (f ∨ g)
      α = β
    end
  end
end

@testset "Variable frequencies" begin
  Φ = [and(collect(1:10)), (1 ∧ 2) ∨ (¬2 ∧ 3) ∨ (¬3 ∧ ¬4), or(collect(1:10)),
       (1 ∨ 2) ∧ (¬2 ∨ 3) ∧ (¬3 ∨ ¬4), bdd_var(1), ¬1, 1 ∨ 2, ¬1 ∨ 2]
  F = Dict{Int, Int}.([[i => 1 for i ∈ 1:10], [1 => 1, 2 => 2, 3 => 3, 4 => 2],
                       [i => 1 for i ∈ 1:10], [1 => 1, 2 => 2, 3 => 3, 4 => 2], [1 => 1], [1 => 1],
                       [1 => 1, 2 => 1], [1 => 1, 2 => 1]])
  for i ∈ 1:length(Φ)
    G = culledfreqs(Φ[i])
    @test G == F[i]
  end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
module DimensionalData

# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end DimensionalData

# Standard lib
using Dates,
      LinearAlgebra,
      Random,
      Statistics,
      SparseArrays

using Base.Broadcast: Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle,
      Unknown

using Base: tail, OneTo, @propagate_inbounds

# Ecosystem
import Adapt,
       ArrayInterface,
       ConstructionBase,
       RecipesBase,
       Tables

using RecipesBase: @recipe

include("Dimensions/Dimensions.jl")

using .Dimensions
using .Dimensions.LookupArrays
using .Dimensions: StandardIndices, DimOrDimType, DimTuple, DimType, AllDims
import .LookupArrays: metadata, set, _set, rebuild, basetypeof,
    order, span, sampling, locus, val, index, bounds, hasselection, units
import .Dimensions: dims, refdims, name, lookup, dimstride, kwdims, hasdim, label

export LookupArrays, Dimensions

# Dimension
export X, Y, Z, Ti, Dim, Coord

# Selector
export At, Between, Contains, Near, Where

export AbstractDimArray, DimArray

export AbstractDimStack, DimStack

export AbstractDimTable, DimTable

export DimIndices, DimKeys, DimPoints

# getter methods
export dims, refdims, metadata, name, bounds

# Dimension/Lookup primitives
export dimnum, hasdim, hasselection, otherdims

# utils
export set, rebuild, reorder, modify, dimwise, dimwise!

const DD = DimensionalData

# Common
include("interface.jl")
include("name.jl")

# Arrays
include("array/array.jl")
include("array/indexing.jl")
include("array/methods.jl")
include("array/matmul.jl")
include("array/broadcast.jl")
include("array/show.jl")
# Stacks
include("stack/stack.jl")
include("stack/indexing.jl")
include("stack/methods.jl")
include("stack/show.jl")
# Other
include("dimindices.jl")
include("tables.jl")
# Combined (easier to work on these in one file)
include("plotrecipes.jl")
include("utils.jl")
include("set.jl")

function _precompile()
    precompile(DimArray, (Array{Int,1}, typeof(X())))
    precompile(DimArray, (Array{Int,2}, Tuple{typeof(X()), typeof(Y())}))
    precompile(DimArray, (Array{Float64,1}, typeof(X())))
    precompile(DimArray, (Array{Float64,2}, Tuple{typeof(X()), typeof(Y())}))
    precompile(DimArray, (Array{Float64,3}, Tuple{typeof(X()), typeof(Y()), typeof(Z())}))
end

_precompile()

end


text:
using JSON

module Bokehjs

const RealVect  = Union{AbstractVector{Int}, AbstractVector{Float64}}
const RealMat   = Union{AbstractMatrix{Int}, AbstractMatrix{Float64}}
const RealArray = Union{RealMat, RealVect}
# would be nice to parameterize, but more important to constrain dims to 1 or 2
# RealArray{N} = Union{AbstractArray{Int, N}, AbstractArray{Float64, N}}

# like nothing except omitted from json rather than being null
mutable struct Omit
    v::AbstractString
    Omit() = new("__omitted from json__")
end
const omit = Omit()

# in case we want to restrict value types in future:
const BkAny      = Any # Union{Dict, Array, AbstractString, Number, Bool, Void, UUID}
const NullDict   = Union{Void, Dict{AbstractString, BkAny}}
const OmitDict   = Union{Omit, Dict{Symbol, BkAny}}
const NullString = Union{Void, AbstractString}
const OmitString = Union{Omit, AbstractString}
const NullSymbol = Union{Void, Symbol}
const OmitSymbol = Union{Omit, Symbol}
const NullFloat  = Union{Float64, Void}
const NullInt    = Union{Int, Void}

const uuid4 = Base.Random.uuid4
const UUID  = Base.Random.UUID

abstract type PlotObject end

abstract type BkRange <: PlotObject end

const NullBkRange = Union{Void, BkRange}

abstract type Renderer <: PlotObject end

abstract type Axis <: PlotObject end

mutable struct TypeID
    plotob::Union{PlotObject, Void}
end

function TypeID()
    TypeID(nothing)
end

mutable struct Plot <: PlotObject
    uuid::UUID
    title::AbstractString
    tools::Vector{BkAny}
    plot_height::Int
    plot_width::Int
    x_range::TypeID
    y_range::TypeID
    # could be Vector{TypeID}?:
    renderers::Vector{BkAny}
    above::Vector{TypeID}
    below::Vector{TypeID}
    left::Vector{TypeID}
    right::Vector{TypeID}
    data_sources::Vector{BkAny}
end

mutable struct ColumnDataSource <: PlotObject
    uuid::UUID
    column_names::Vector{Symbol}
    selected::Vector{Any}
    discrete_ranges::Dict{Symbol, BkAny}
    cont_ranges::Dict{Symbol, BkAny}
    data::Dict{Symbol, Vector}
end

function ColumnDataSource(data::Dict{Symbol, Vector})
    ColumnDataSource(uuid4(),
                     collect(keys(data)),
                     BkAny[],
                     Dict{Symbol, BkAny}(),
                     Dict{Symbol, BkAny}(),
                     data)
end

mutable struct DataRange1d <: BkRange
    uuid::UUID
    sources::Vector{BkAny}
end

function DataRange1d(cdss::Vector{ColumnDataSource}, columns::Vector{AbstractString})
    source(cds) = Dict{AbstractString, BkAny}("columns" => columns,
                                      "source" => TypeID(cds))
    sources = map(source, cdss)
    DataRange1d(uuid4(), sources)
end

mutable struct TickFormatter <: PlotObject
    uuid::UUID
    _type_name::Symbol
    format::OmitDict
    function TickFormatter(name::Symbol)
        @assert name in (:BasicTickFormatter, :DatetimeTickFormatter,
                         :LogTickFormatter)
        # format only seems to occur for DatetimeTickFormatter and even then is empty
        format = name == :DatetimeTickFormatter ? Dict{Symbol, BkAny}() : omit
        new(uuid4(), name, format)
    end
end

mutable struct Ticker <: PlotObject
    uuid::UUID
    _type_name::Symbol
    num_minor_ticks::Int64
    Ticker(name::Symbol) = new(uuid4(), name, 5)
end

mutable struct LinearAxis <: Axis
    uuid::UUID
    dimension::Int
    bounds::AbstractString
    location::AbstractString
    formatter::TypeID
    ticker::TypeID
    plot::TypeID
end

function LinearAxis(dimension::Int, tf::TickFormatter, t::Ticker, plot::Plot)
    LinearAxis(uuid4(),
               dimension,
               "auto",
               "min",
               TypeID(tf),
               TypeID(t),
               TypeID(plot))
end

mutable struct Grid <: Renderer
    uuid::UUID
    dimension::Int
    plot::TypeID
    ticker::TypeID
end

function Grid(dimension::Int, plot::Plot, ticker::Ticker)
    Grid(uuid4(), dimension, TypeID(plot), TypeID(ticker))
end

mutable struct Legend <: Renderer
    uuid::UUID
    plot::TypeID
    legends::Vector{Tuple}
    orientation::OmitSymbol
end

function Legend(plot::Plot, legends::Vector{Tuple})
    Legend(plot, legends, nothing)
end

function Legend(plot::Plot, legends::Vector{Tuple}, orientation::NullSymbol)
    orientation = orientation == nothing ? omit : orientation
    @assert orientation in (omit, :top_left, :top_center, :top_right,
                            :right_center, :bottom_right, :bottom_center,
                            :bottom_left, :left_center, :center)
    Legend(uuid4(),
           TypeID(plot),
           [(l, [TypeID(r)]) for (l, r) in legends],
           orientation)
end

mutable struct Glyph <: PlotObject
    uuid::UUID
    _type_name::Symbol
    line_color::OmitDict
    line_width::OmitDict
    line_alpha::OmitDict
    fill_color::OmitDict
    fill_alpha::OmitDict
    size::OmitDict
    line_dash::Union{Omit, Vector{Int64}}
    x::Dict{Symbol, Symbol}
    y::Dict{Symbol, Symbol}
end

mutable struct GlyphRenderer <: Renderer
    uuid::UUID
    data_source::TypeID
    nonselection_glyph::TypeID
    selection_glyph::TypeID
    glyph::TypeID
    name::Union{Void, AbstractString}
    server_data_source::NullDict
end

const NullGlyph = Union{Void, Glyph}

function GlyphRenderer(coldata::ColumnDataSource, nonsel_g::NullGlyph,
                       sel_g::NullGlyph, glyph::Glyph)
    GlyphRenderer(uuid4(),
                  TypeID(coldata),
                  TypeID(nonsel_g),
                  TypeID(sel_g),
                  TypeID(glyph),
                  nothing,
                  nothing
                 )
end

mutable struct Metatool <: PlotObject
    uuid::UUID
    _type_name::AbstractString
    plot::TypeID
    dimensions::Union{Vector{AbstractString}, Void, Omit}
end

function Metatool(typename::AbstractString, plot::Plot, dimensions)
    plot = TypeID(plot)
    Metatool(uuid4(), typename, plot, dimensions)
end

function Metatool(typename::AbstractString, plot::Plot)
    Metatool(typename, plot, omit)
end

# TODO: these are duplicated as HEIGHT and WIDTH in Bokeh.jl
_DEFAULT_HEIGHT = 600
_DEFAULT_WIDTH = 600

function Plot()
    Plot(uuid4(),
         "",
         Void[],
         _DEFAULT_HEIGHT,
         _DEFAULT_WIDTH,
         TypeID(),
         TypeID(),
         Void[],

         Void[],
         Void[],
         Void[],
         Void[],

         Void[]
        )
end

function Plot(plot::Plot,
              xrange::BkRange,
              yrange::BkRange,
              renderers::Array{PlotObject,1},
              axes,#::Dict{Symbol, Array{PlotObject,1}},
              tools::Array{PlotObject,1},
              title::AbstractString="Bokeh Plot",
              height::Int=_DEFAULT_HEIGHT,
              width::Int=_DEFAULT_WIDTH)
    data_sources = BkAny[]# [TypeID(coldata)]

    Plot(plot.uuid,
         title,
         map(TypeID, tools),
         height,
         width,
         TypeID(xrange),
         TypeID(yrange),
         map(TypeID, renderers),
         map(TypeID, axes[:above]),
         map(TypeID, axes[:below]),
         map(TypeID, axes[:left]),
         map(TypeID, axes[:right]),
         data_sources
        )
end

mutable struct PlotContext <: PlotObject
    uuid::UUID
    children::Vector{TypeID}
end

function PlotContext(plot::Plot)
    PlotContext(uuid4(),[TypeID(plot)])
end

end  # module

# extract useful types from Bokehjs module
const RealVect   = Bokehjs.RealVect
const RealMat    = Bokehjs.RealMat
const RealArray  = Bokehjs.RealArray
const BkAny      = Bokehjs.BkAny
const omit       = Bokehjs.omit
const Glyph      = Bokehjs.Glyph
const NullString = Bokehjs.NullString
const NullSymbol = Bokehjs.NullSymbol
const NullFloat  = Bokehjs.NullFloat
const NullInt    = Bokehjs.NullInt

JSON.lower(uuid::Bokehjs.UUID) = string(uuid)

function JSON.lower(tid::Bokehjs.TypeID)
    tid.plotob == nothing && return nothing
    attrs = fieldnames(tid.plotob)
    obtype = in(:_type_name, attrs) ? tid.plotob._type_name : typeof(tid.plotob)
    Dict{AbstractString, BkAny}("type" => obtype, "id" => tid.plotob.uuid)
end

JSON.lower(::Type{T}) where {T<:Bokehjs.PlotObject} = string(T.name.name)


text:
module ComplexRegions

DEFAULT = Dict(:tol=>1e-12)
"""
	ComplexRegions.default()
Return a dictionary of global default settings for the ComplexRegions package.

	ComplexRegions.default(key=value)
Change a global default setting in the running instance of the ComplexRegions package.
"""
function default(;kw...)
	if isempty(kw)
		return DEFAULT
	end
	for (sym,val) in kw
		if !haskey(DEFAULT,sym)
			@error "Unrecognized default setting `$sym`"
		else
			DEFAULT[sym] = val
			@info "Default value of `$sym` set to $val."
		end
	end
	return nothing
end

using Reexport
@reexport using ComplexValues
using LinearAlgebra,StaticArrays

AnyComplex = Union{Complex{S},Polar{S},Spherical{S}} where {S<:AbstractFloat}

export Polar,Spherical
import Base: +,-,*,/,!,∘,sign,inv,angle,real,imag,conj,show,iterate,eltype,length,getindex,isapprox,isfinite,intersect,union,truncate,reverse,in
export +,-,*,/,!,∘,sign,inv,angle,real,imag,conj,show,iterate,eltype,length,getindex,isapprox,isfinite,intersect,union,truncate,reverse,in

include("utilities.jl")

export point,arclength,slope,dist,closest,isleft,plotdata,isright,reflect,tangent,unittangent,normal,arg,isinside,isoutside
export Curve,ClosedCurve,Circle,Line,Arc,Segment,Ray
include("curves.jl")

export Path,ClosedPath,Polygon,CircularPolygon
export curve,curves,vertex,vertices,side,sides,isfinite,sign,angle,angles,winding,ispositive,rectangle,n_gon
include("paths.jl")

export SimplyConnectedRegion,ConnectedRegion,ExteriorRegion,disk,unitdisk,Annulus,PolygonalRegion
export halfplane,upperhalfplane,lowerhalfplane,lefthalfplane,righthalfplane
export RegionIntersection,RegionUnion
export region,interior,exterior,between,boundary,innerboundary,outerboundary
include("regions.jl")
include("simplyconnected.jl")

export Möbius,Mobius
include("mobius.jl")

include("docs.jl")
include("plotrecipes.jl")

end # module


text:
module MosekBlocks
import Mosek

export Model,variable!

mutable struct ManagedIndexes
    next       :: Vector{Int}
    prev       :: Vector{Int}
    first_used :: Int
    last_used  :: Int
    first_free :: Int
    last_free  :: Int

    nfree      :: Int
end

ManagedIndexes() = ManagedIndexes( Int[], Int[], 0, 0, 0, 0, 0 )

Base.length(mi :: ManagedIndexes) = length(mi.next)

function ensure!(self :: ManagedIndexes, num :: Int)
    if self.nfree < num
        addnum = num-self.nfree
        oldnum = length(self.next)
        newnum = oldnum + addnum
        resize!(self.next,newnum)
        resize!(self.prev,newnum)

        for i in oldnum+1:newnum
            self.next[i] = i+1
            self.prev[i] = i-1
        end
        first = oldnum+1
        last  = newnum

        self.next[last] = 0
        self.prev[first] = self.last_free
        if self.last_free == 0
            self.first_free = first
            self.last_free = last
        else
            self.next[self.last_free] = first
            self.last_free = last
        end
        self.nfree += addnum
    end
end

function alloc!(self :: ManagedIndexes, num :: Int)
    ensure!(self,num)
    first = self.first_free
    last = first
    for i in 1:num-1
        last = self.next[last]
    end

    if self.next[last] > 0
        self.prev[self.next[last]] = 0
    end

    self.first_free = self.next[last]
    if self.first_free == 0
        self.last_free = 0
    end

    self.next[last] = 0
    self.prev[first] = self.last_used
    if self.last_used == 0
        self.first_used = first
    else
        self.next[self.last_used] = first
    end
    self.last_used = last
    self.nfree -= num

    first
end

function free(mi :: ManagedIndexes, index :: Int, num :: Int)
    first = index
    last = index
    for i in 1:num-1
        last = mi.next[last]
    end

    if mi.prev[first] == 0
        if mi.next[last] == 0
            mi.first_used = 0
            mi.last_used = 0
        else
            mi.first_used = mi.next[last]
            mi.prev[mi.first_used] = 0
        end
    else
        if mi.next[last] == 0
            mi.last_used = mi.prev[first]
            mi.next[mi.last_used] = 0
        else
            mi.prev[mi.next[last]] = mi.prev[first]
            mi.next[mi.prev[first]] = mi.next[last]
        end
    end

    if mi.first_free == 0
        mi.first_free = first
        mi.last_free = last
    else
        mi.next[last] = mi.first_free
        mi.prev[mi.first_free] = last
        mi.first_free = first
    end

    mi.nfree += num
end

function get_indexes(self::ManagedIndexes, head::Int, num::Int, res::Vector{Int})
    ii = head
    res[1] = ii
    for i in 2:num
        ii = self.next[ii]
        res[i] = ii
    end
end
#----------------------------
mutable struct ManagedBlocks
    ii :: ManagedIndexes
    blocklist :: ManagedIndexes

    blocksize :: Vector{Int}
    blockhead :: Vector{Int}
end

ManagedBlocks() = ManagedBlocks(ManagedIndexes(),ManagedIndexes(),Int[],Int[])

function alloc!(self :: ManagedBlocks, num :: Int)
    p       = alloc!(self.ii,num)
    blockid = alloc!(self.blocklist,1)

    numblocks = length(self.blocklist)
    if length(self.blocksize) < numblocks
        resize!(self.blocksize,numblocks)
        resize!(self.blockhead,numblocks)
    end
    self.blocksize[blockid] = num
    self.blockhead[blockid] = p

    blockid
end

function free!(self :: ManagedBlocks, blockid :: Int)
    free!(self.ii, self.blockhead[blockid], self.blocksize[blockid])
end

Base.length(self::ManagedBlocks) = length(self.ii)

num_blocks(self::ManagedBlocks) = length(self.blocklist)

function get_indexes(self::ManagedBlocks, blockid :: Int)
    num = self.blocksize[blockid]
    head = self.blockhead[blockid]
    res = zeros(Int,num)
    get_indexes(self.ii, head, num, res)
    res
end

type Model
    task :: Mosek.MSKtask

    varblocks :: ManagedBlocks
    conblocks :: ManagedBlocks

    conbfix   :: Vector{Float64}
    conslack  :: Vector{Int}
end

Model() = Model(Mosek.maketask(),
                ManagedBlocks(),
                ManagedBlocks(),
                Float64[],
                Float64[] )

type Variable{N}
    shape   :: NTuple{N,Int}
    indexes :: Vector{Int}
end

type Constraint{N}
    shape :: NTuple{N,Int}
    indexes :: Vector{Int}
end

type Expr{N}
    shape :: NTuple{N,Int}
    ptr   :: Vector{Int}
    nidx  :: Vector{Int}
    cof   :: Vector{Float64}
    bfix  :: Vector{Float64}
end

function alloc_var!(m::Model, num :: Int)
    numvar = Mosek.getnumvar(m.task)
    blockid = alloc!(m.varblocks,num)
    if length(m.varblocks) < numvar
        nadd = numvar-length(m.varblocks)
        Mosek.appendvars(m.task,nadd)
        Mosek.putvarboundslice(m.task,numvar+1,numvar+nadd+1,
                               fill(Mosek.MSK_BK_FR,nadd),
                               zeros(Float64,nadd),
                               zeros(Float64,nadd))
    end
    blockid
end

function alloc_con!(m::Model, num :: Int)
    numcon = Mosek.getnumcon(m.task)
    blockid = alloc!(m.varblocks,num)
    if length(m.conblocks) < numvar
        nadd = numcon-length(m.conblocks)
        newnumcon = numcon+nadd
        Mosek.appendcons(m.task,nadd)
        Mosek.putconboundslice(m.task,numcon+1,numcon+nadd+1,
                               fill(Mosek.MSK_BK_FR,nadd),
                               zeros(Float64,nadd),
                               zeros(Float64,nadd))
        resize!(m.conbfix,newnumcon)
        resize!(m.conslack,newnumcon)

        m.conbfix[numcon+1:]  = 0.0
        m.conslack[numcon+1:] = 0
    end
    blockid
end

function variable!{N}(m::Model, shape :: NTuple{N,Int})
    num = prod(shape)
    blockid = alloc_var!(m,num)
    idxs = get_indexes(m.varblocks, blockid)
    Variable( shape, idxs )
end

variable!(m::Model, shape :: Int) = variable!(m,(shape,))

function constraint!(m:Model, shape :: NTuple{N,Int}, expr :: Expr, dom ::D)
    where N, D <: AbstractDomain

end

###############
# TEST

m = MosekBlocks.Model()

v1 = MosekBlocks.variable!(m, (2,2))
v2 = MosekBlocks.variable!(m,2)


text:
using NonconvexCore, LinearAlgebra, Test, Zygote

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

m = Model(f)
addvar!(m, [0.0, 0.0], [10.0, 10.0])
add_ineq_constraint!(m, x -> g(x, 2, 0))
add_ineq_constraint!(m, x -> g(x, -1, 1))

@test getmin(m) == fill(0, 2)
@test getmax(m) == fill(10, 2)
@test getmin(m, 1) == 0.0
@test getmax(m, 1) == 10.0
@test getmin(m, 2) == 0.0
@test getmax(m, 2) == 10.0

p = [1.234, 2.345]

val0, grad0 = NonconvexCore.value_gradient(NonconvexCore.getobjective(m), p)
@test val0 == f(p)
@test grad0 == Zygote.gradient(f, p)[1]

val1, grad1 = NonconvexCore.value_gradient(NonconvexCore.getineqconstraint(m, 1), p)
@test val1 == g(p, 2, 0)
@test grad1 == Zygote.gradient(p -> g(p, 2, 0), p)[1]

val2, grad2 = NonconvexCore.value_gradient(NonconvexCore.getineqconstraint(m, 2), p)
@test val2 == g(p, -1, 1)
@test grad2 == Zygote.gradient(p -> g(p, -1, 1), p)[1]

vals, jac = NonconvexCore.value_jacobian(NonconvexCore.getineqconstraints(m), p)
@test [val1, val2] == vals
@test [grad1 grad2]' == jac

vals, jac = NonconvexCore.value_jacobian(NonconvexCore.getobjectiveconstraints(m), p)
@test [val0, val1, val2] == vals
@test [grad0 grad1 grad2]' == jac


text:
## UNDER CONSTRUCTION!
"""Matrix A such as A ⋅ b = n × b"""
cross(n) = begin
   nx,ny,nz = n
    [0 -nz ny;
     nz 0 -nx;
    -ny nx 0]
end
"""Rodrigues' formula: Rotation matrix that when applied rotates with respect to "n" in an angle θ anti clock-wise"""
Un(θ,n) = I * cos(θ) + sin(θ) * cross(n) + (1-cos(θ)) * (n * n')
"""
Slab oriented along the x axis.

Bar, L., & Sochen, N. (2015). A spectral framework for NMR signal with restricted diffusion. Concepts in Magnetic Resonance Part A, 44(1), 16–53. doi:10.1002/cmr.a.21326
Grebenkov, D.S. (2008), Laplacian eigenfunctions in NMR. I. A numerical tool. Concepts Magn. Reson., 32A: 277-301. doi:10.1002/cmr.a.20117
"""
function Planes(L,D=2e-9,M=30)
    #Normalization
    ϵ(m) = m==0 ? 1 : sqrt(2)
    #Eigen values
    λ(m) = π^2*m^2
    Λ =  D/L^2 * diagm([λ(m) for m=0:M])
    #Matrix A
    Ax = L * [m!=n ? ((-1)^(m+n)-1)*ϵ(m)*ϵ(n)*(λ(m)+λ(n))/(λ(m)-λ(n))^2 : 1/2 for m=0:M, n=0:M]
    A = cat(Ax,0*Ax,0*Ax; dims=3)
    n = [1 0 0;
         0 1 0;
         0 0 1]
    μ = (Λ, A, n)
end
"""
Infinite cylinder oriented along the z axis.

Bar, L., & Sochen, N. (2015). A spectral framework for NMR signal with restricted diffusion. Concepts in Magnetic Resonance Part A, 44(1), 16–53. doi:10.1002/cmr.a.21326
Grebenkov, D.S. (2008), Laplacian eigenfunctions in NMR. I. A numerical tool. Concepts Magn. Reson., 32A: 277-301. doi:10.1002/cmr.a.20117
"""
function Cylinder(R,D=2e-9,M=20)
    #J'n(αnk) = 0
    α = [1.841184, 3.054237, 3.831706, 4.201189, 5.317553,
         5.331443, 6.415616, 6.706133, 7.015587, 7.501266,
         8.015237, 8.536316, 8.577836, 9.282396, 9.647422,
         9.969468, 10.17347, 10.51986, 10.71143, 11.34592]
    n = [1,2,0,3,4, 1,5,2,0,6, 3,1,7,4,8, 2,0,5,9,3]
    #Eigen values
    λ(nk) = nk==0 ? 0 : α[nk]^2
    Λ =  D/R^2 * diagm([λ(nk) for nk=0:M])
    #βnk
    β(nk) = nk==0 ? 1 : sqrt(λ(nk)/(λ(nk)-n[nk]^2))
    #Some definitions
    cond1(mk,nk) = abs(n[mk] - n[nk]) == 1
    ϵ1(mk,nk) = sqrt( 1 + (n[mk] == 0) + (n[nk] == 0) )
    cond2(mk,nk) = (n[mk] == n[nk]-1) - (n[mk] == n[nk]+1)
    #Matrix A
    Ax = R * [cond1(i,j)*ϵ1(i,j)*β(i)*β(j)*(λ(i)+λ(j)-2*n[i]*n[j])/(λ(i)-λ(j))^2 for i=0:M, j=0:M]
    Ay = 1im*R * [cond2(i,j)*β(i)*β(j)*(λ(i)+λ(j)-2*n[i]*n[j])/(λ(i)-λ(j))^2 for i=0:M, j=0:M]
    A = (Ax,Ay,0*Ax)
    μ = (Λ, A)
end

"""
Sphere of radius R.

Bar, L., & Sochen, N. (2015). A spectral framework for NMR signal with restricted diffusion. Concepts in Magnetic Resonance Part A, 44(1), 16–53. doi:10.1002/cmr.a.21326
Grebenkov, D.S. (2008), Laplacian eigenfunctions in NMR. I. A numerical tool. Concepts Magn. Reson., 32A: 277-301. doi:10.1002/cmr.a.20117
"""
function Sphere(R,D=2e-9,M=20)
    #j'n(αlk) = 0
    α = [2.081576, 3.342094, 4.493409, 4.514100, 5.646704,
         5.940370, 6.756456, 7.289932, 7.725252, 7.851078,
         8.583755, 8.934839, 9.205840, 9.840446, 10.01037,
         10.61386, 10.90412, 11.07021, 11.07942, 11.97273]
    l = [1,2,0,3,4, 1,5,2,0,6, 3,7,1,4,8, 2,0,5,9,3]
    m = [0,0,1,0,0, 1,0,1,2,0, 1,0,2,1,0, 2,3,1,0,2]
    #Eigen values
    λ(lm) = lm==0 ? 0 : α[lm]^2
    Λ =  D/R^2 * diagm([λ(lm) for lm=0:M])
    #βnk
    β(lm) = lm==0 ? sqrt(3/2) : sqrt((2*l[lm]+1)*λ(lm)/(λ(lm)-l[lm]*(l[lm]+1)))
    # TODO: Finish matrix Ax
    δlδm(i,j) = (abs(l[i]-l[j])==1)*(abs(m[i]-m[j])==1)
    ϵ1(mk,nk) = (1+n[mk]+n[nk])/((2*n[mk]+1)*(2*n[nk]+1))
    ϵ2(mk,nk) = n[mk]*(n[nk]+1) + n[nk]*(n[mk]+1) + 1
    #Matrix A
    Ax = R * [cond1(mk,nk)*ϵ1(mk,nk)*β(mk)*β(nk)*(λ(mk)+λ(nk)-ϵ2(mk,nk))/(λ(mk)-λ(nk))^2 for mk=0:M, nk=0:M]
    A = (Ax,Ay,Az)
    μ = (Λ, A)
end

function SignalE(μ, seq)
    𝒊 = 1im;
    M, N = size(seq.GR)
    G = getproperty.(seq.GR,:A)
    δ = getproperty.(seq.GR[1,:],:T)
    # E = [ Π exp( -(Λ + iγ Gn⋅A) ⋅ δn ) ]_{0,0}
    A = [exp(-( μ[1] .+ 𝒊*2π*γ* .+([μ[2][:,:,m]*G[m,n] for m=1:M]...) )*δ[n]) for n=1:N]
    E = *(A...)[1,1]
end


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf
using ImPlot
import CImGui.LibCImGui: ImGuiCond_Always, ImGuiCond_Once

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

window = GLFW.CreateWindow(1280, 720, "ImPlot Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync; set to 0 to benchmark

ctx = CImGui.CreateContext()
ctxp = ImPlot.CreateContext()
ImPlot.SetImGuiContext(ctx)

CImGui.StyleColorsDark()

fonts_dir = joinpath(@__DIR__, "..", "fonts")
fonts = CImGui.GetIO().Fonts
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)

ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    show_another_window = false
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    # make up some data
    xs1 = Float64.(collect(1:100))
    ys1 = rand(6000)
    noise = rand(length(xs1))
    bar_vals = Vector{Vector{Float64}}(undef, 5)
    bar_maxes = zeros(5)
    bar_mins = zeros(5)
    for i in 1:length(bar_vals)
        bar_vals[i] = Float64.(collect(range(rand(1:10), step=rand([-2,-1,1,2]), length=120)))
        bar_maxes[i] = maximum(bar_vals[i])
        bar_mins[i] = minimum(bar_vals[i])
    end
    bar_max = maximum(bar_maxes)
    bar_min = minimum(bar_mins)
    bar_counter = 1
    # main loop
    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()

        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        CImGui.Begin("Example Plots")

        @c CImGui.Checkbox("Show Examples", &show_another_window)
        CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)",
                             1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))

        CImGui.End()

        if show_another_window
            @c CImGui.Begin("Examples Window", &show_another_window)
            if CImGui.CollapsingHeader("Line plots")

                ys1 .= rand(6000)
                ImPlot.SetNextPlotLimits(0.0, 6000, 0.0, 1.0, ImGuiCond_Always)
                # Using '##' in the label name hides the plot label, but lets
                # us keep the label ID unique for modifying styling etc.
                if ImPlot.BeginPlot("##line", "x1", "y1", CImGui.ImVec2(-1,300))
                    ImPlot.PlotLine(ys1)
                    ImPlot.EndPlot()
                end
            end
            if CImGui.CollapsingHeader("Scatter plot")
                noise .= xs1 .+ rand(-5.0:0.1:5.0, length(xs1))
                ImPlot.SetNextPlotLimits(0,100,-5,105, ImGuiCond_Always)
                if ImPlot.BeginPlot("##scatter", "x2", "y2", CImGui.ImVec2(-1,300))
                    ImPlot.PlotScatter(xs1, noise)
                    ImPlot.EndPlot()
                end
            end
            if CImGui.CollapsingHeader("Bar plot")

                bar_val_step = [bar_vals[j][bar_counter] for j in 1:length(bar_vals)]

                ImPlot.SetNextPlotLimits(-0.5,4.5,bar_min, bar_max, ImGuiCond_Always)
                if ImPlot.BeginPlot("##bars", "", "", CImGui.ImVec2(-1,300))
                    ImPlot.PlotBars(bar_val_step)
                    ImPlot.EndPlot()
                end

                if bar_counter == 120
                    bar_counter = 1
                else
                    bar_counter += 1
                end
            end
            if CImGui.CollapsingHeader("Shaded plot")
                x = 1:1000
                y1 = [sin(x) for x in range(0,2π, length = length(x))]
                y_ref = -2.0
                ImPlot.SetNextPlotLimits(0,1000,-2,1, ImGuiCond_Always)
                if ImPlot.BeginPlot("##shaded", "", "", CImGui.ImVec2(-1,300))
                    ImPlot.PlotShaded(x, y1, y_ref)
                    ImPlot.EndPlot()
                end
            end
            CImGui.End()
        end

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end
catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    ImPlot.DestroyContext(ctxp)
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end


text:
using QuantLib

function main()
  todaysDate = Date(2007, 10, 16)
  set_eval_date!(settings, todaysDate)

  bbCurveRate = 0.055
  bbDayCount = QuantLib.Time.ActualActualBond()
  bbIR = InterestRate(bbCurveRate, bbDayCount, CompoundedCompounding(), QuantLib.Time.Semiannual())
  ts = FlatForwardTermStructure(todaysDate, bbIR.rate, bbIR.dc, bbIR.comp, bbIR.freq)

  # set up call schedule
  callSchedule = CallabilitySchedule(24)
  callPrice = 100.0
  callDate = Date(2006, 9, 15)

  for i in eachindex(callSchedule)
    nullCal = QuantLib.Time.NullCalendar()
    myPrice = Price(callPrice, CleanCall())

    callSchedule[i] = Callability(myPrice, Call(), callDate)
    callDate = QuantLib.Time.advance(Dates.Month(3), nullCal, callDate)
  end

  # set up callable bond
  dated = Date(2004, 9, 16)
  issue = dated
  mat = Date(2012, 9, 15)

  settlementDays = 3
  bondCal = QuantLib.Time.USGovernmentBondCalendar()
  coupon = 0.0465
  freq = QuantLib.Time.Quarterly()
  redemption = 100.0
  faceAmount = 100.0

  bondDayCount = QuantLib.Time.ActualActualBond()
  accrualConvention = QuantLib.Time.Unadjusted()
  paymentConvention = QuantLib.Time.Unadjusted()

  schedule = QuantLib.Time.Schedule(dated, mat, QuantLib.TenorPeriod(freq), accrualConvention, accrualConvention, QuantLib.Time.DateGenerationBackwards(), false, bondCal)

  maxIter = 1000
  accuracy = 1e-8
  gridIntervals = 40
  reversionParameter = 0.03

  sigma = eps()
  hw0 = HullWhite(ts, reversionParameter, sigma)
  engine0 = TreeCallableFixedRateEngine(hw0, gridIntervals)

  callableBond = CallableFixedRateBond(settlementDays, faceAmount, schedule, coupon, bondDayCount, paymentConvention, redemption, issue, callSchedule, engine0)

  println(clean_price(callableBond))
  println(QuantLib.yield(callableBond, bondDayCount, CompoundedCompounding(), freq, accuracy, maxIter))
end