text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_storage(path::AbstractString, inputs::Dict,setup::Dict, EP::Model)

Function for writing the capacities of different storage technologies, including hydro reservoir, flexible storage tech etc.
"""
function write_storage(path::AbstractString, inputs::Dict,setup::Dict, EP::Model)
	dfGen = inputs["dfGen"]
	T = inputs["T"]     # Number of time steps (hours)
	G = inputs["G"]

	# Storage level (state of charge) of each resource in each time step
	dfStorage = DataFrame(Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone])
	s = zeros(G,T)
	storagevcapvalue = zeros(G,T)
	for i in 1:G
		if i in inputs["STOR_ALL"]
			s[i,:] = value.(EP[:vS])[i,:]
		elseif i in inputs["HYDRO_RES"]
			s[i,:] = value.(EP[:vS_HYDRO])[i,:]
		elseif i in inputs["FLEX"]
			s[i,:] = value.(EP[:vS_FLEX])[i,:]
		end
	end

	# Incorporating effect of Parameter scaling (ParameterScale=1) on output values
	for y in 1:G
		if setup["ParameterScale"]==1
			storagevcapvalue[y,:] = s[y,:].*ModelScalingFactor
		else
			storagevcapvalue[y,:] = s[y,:]
		end
	end

	dfStorage = hcat(dfStorage, DataFrame(storagevcapvalue, :auto))
	auxNew_Names=[Symbol("Resource");Symbol("Zone");[Symbol("t$t") for t in 1:T]]
	rename!(dfStorage,auxNew_Names)
	CSV.write(joinpath(path, "storage.csv"), dftranspose(dfStorage, false), writeheader=false)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using LinearAlgebra
using MKL
import Fermi.Orbitals: AbstractRestrictedOrbitals

function orbital_eval(orb::AbstractRestrictedOrbitals, bset::BasisSet, x, y, z, n)

    atoms = orb.molecule.atoms

    # Initialize array with basis function value
    χ = zeros(bset.nbas)

    i = 1
    for a in eachindex(atoms)
        basis_on_atom_a = bset[a]
        x0, y0, z0 = atoms[a].xyz
        for b in eachindex(basis_on_atom_a)
            bf = basis_on_atom_a[b]
            if bf.l == 1
                ml_vals = [1, -1, 0]
            else
                ml_vals = collect(-bf.l:1:bf.l)
            end
            for ml in ml_vals
                χ[i] = orbital_eval(bf, x,y,z, ml, x0=x0, y0=y0, z0=z0)
                i += 1
            end
        end
    end

    Cμi = orb.C[:,n]

    return sum(Cμi .* χ)
end

function isoplot(orb::AbstractRestrictedOrbitals, min, max, N)

    fig = Figure(resolution=(900,800))
    lscene = LScene(fig[1,1], tellheight=false, height=700, width=600)

    # Plot atoms
    atoms = orb.molecule.atoms
    plot!(lscene, atoms)

    gui_grid = fig[1,2] = GridLayout(width = 250, tellheigth=false)

    # Orbital select menu
    Label(gui_grid[1,1][1,1], "Orbital")
    menu = Menu(gui_grid[1,1][1,2], options = collect(1:size(orb.C,1)))

    # Isovalue slider
    ls = labelslider!(fig, "Isovalue", 0:0.01:0.5)#, tellheight = false)
    gui_grid[2,1] = ls.layout

    # Toggles: Show axis
    togglegrid = gui_grid[3,1] = GridLayout(tellwidth=false)
    ax_t = Toggle(togglegrid[1,1])
    Label(togglegrid[1,2], "Show axis")

    ax_t.active[] = true
    on(ax_t.active) do visible
        lscene.scene.plots[1].plots[1].visible[] = visible
        lscene.scene.plots[1].plots[2].visible[] = visible
    end

    # Save button
    sbut = Button(gui_grid[4,1], label="Save", tellwidth=false)
    on(sbut.clicks) do n
        lscene.scene.center = false
        save("atomicorbital.png", lscene.scene)
    end

    set_close_to!(ls.slider, 0.48)
    menu.selection[] = 1

    # Positive and negative isovalues conected as nodes
    iv = ls.slider.value
    negiv = lift(iv) do x
        -x
    end

    # Get basis set info
    bset = BasisSet(orb.basis, atoms)
    nbas = bset.nbas

    # Allocate array for basis set grids
    χ = zeros(Float32, N, N, N, size(orb.C,1))
    Cμi = Float32.(orb.C)

    bas_per_atom = Int[]
    for a in eachindex(atoms)
        basis_on_atom_a = bset[a]
        nbf = 0
        for b in basis_on_atom_a
            nbf += 2*b.l + 1
        end
        push!(bas_per_atom, nbf)
    end

    println(bas_per_atom)
    println(nbas)
    @assert sum(bas_per_atom) == nbas

    offset = [sum(bas_per_atom[1:(a-1)]) for a in eachindex(atoms)]

    # Compute χ(x,y,z,μ) rank-4 tensor: x,y,z are the cartesian coordinates
    # μ is the basis function index. Thus, this tensor stores the value of all
    # basis functions over the whole grid
    r = LinRange(min, max, N)
    Threads.@threads for ix in 1:N
        x = r[ix]
        @inbounds for iy in 1:N, iz in 1:N
        y, z = r[iy], r[iz]
        for a in eachindex(atoms)
            basis_on_atom_a = bset[a]
            x0, y0, z0 = atoms[a].xyz
            μ = offset[a] + 1
            for b in eachindex(basis_on_atom_a)
                bf = basis_on_atom_a[b]
                if bf.l == 1
                    ml_vals = [1, -1, 0]
                    #ml_vals = [1, 0, -1]
                else
                    ml_vals = collect(-bf.l:1:bf.l)
                end
                for ml in ml_vals
                    χ[ix,iy,iz,μ] = orbital_eval(bf, x,y,z, ml, x0=x0, y0=y0, z0=z0)
                    μ += 1
                end
            end
        end
        end
    end

    # Perform the AO -> MO transformation as a matrix multiplication
    # C(x,y,z;i) = χ(x,y,z;μ) * C(μ;i)
    χ = reshape(χ, (N^3, nbas))
    Cmo = reshape(χ * Cμi, (N, N, N, nbas))

    # The menu for orbital selection simply slice
    # the full tensor
    wf = lift(menu.selection) do n
        Cmo[:,:,:,n]
    end

    # Plot isosurfaces
    volume!(r, r, r, wf, algorithm = :iso, isorange = 0.01, tellheight=false,
            isovalue =  iv, show_axis=true, colormap = colormap("reds"), colorrange = (0, 0.5))
    volume!(r, r, r, wf, algorithm = :iso, isorange = 0.01, tellheight=false,
            isovalue = negiv, show_axis=true, colormap = colormap("blues"), colorrange = (0, -0.5))

    return fig
end

ext: jl

max_stars_repo_name: gustavojra/Alchemy

max_issues_repo_name: gustavojra/Alchemy

max_forks_repo_name: gustavojra/Alchemy


text:
using DashBootstrapComponents, DashHtmlComponents

text_input = html_div([
    dbc_input(id="input", placeholder="Type something...", type="text"),
    html_br(),
    html_p(id="output"),
]);

callback!(app, Output("output", "children"), Input("input", "value")) do value
    return value
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using AccelerationBenchmark, DataFrames, CSV, Optim, LineSearches, CUTEst, JLD
const run_tests_more = true # Run the Moré et al. tests from the O-ACCEL paper
const run_tests_cutest = true # Run CUTEst tests

const savejld = true
const savecsv = true

function lstests()
    lsm = MoreThuente()
    lsh = HagerZhang()
    lsbt2 = BackTracking(order=2)
    lsbt3 = BackTracking(order=3)

    ag = InitialStatic(alpha=1.0, scaled=false)

    solvers = [
        LBFGS(alphaguess = ag, linesearch = lsm, scaleinvH0=true),
        LBFGS(alphaguess = ag, linesearch = lsm, scaleinvH0=false),
        LBFGS(alphaguess = ag, linesearch = lsh, scaleinvH0=true),
        LBFGS(alphaguess = ag, linesearch = lsh, scaleinvH0=false),
        LBFGS(alphaguess = ag, linesearch = lsbt2, scaleinvH0=true),
        LBFGS(alphaguess = ag, linesearch = lsbt2, scaleinvH0=false),
        LBFGS(alphaguess = ag, linesearch = lsbt3, scaleinvH0=true),
        LBFGS(alphaguess = ag, linesearch = lsbt3, scaleinvH0=false),
    ]
    solvernames = [
        "(MT,true)", "(MT,false)",
        "(HZ,true)", "(HZ,false)",
        "(BT2,true)", "(BT2,false)",
        "(BT3,true)", "(BT3,false)",
    ]
    #stoptype = :FunctionTolRelative
    #stoptol  = 1e-10
    stoptype = :GradientTolRelative
    stoptol  = 1e-10
    timelog  = 3 # timelog = 1 can cause very inaccurate timings due to compilation and garbage collection
    maxiter  = 5000
    timelimit = NaN
    TestSetup(solvers,solvernames,stoptype,stoptol,timelog,maxiter,timelimit)
end

testset = lstests()
savebase = "data/"

if run_tests_more
    tests = [(:A,100), (:A,200),
             (:B,100), (:B,200),
             (:C,100), (:C,200),
             (:D,500), (:D,1000),
             (:E,100), (:E,200),
             (:F,200), (:F,500),
             (:G,100), (:G,200),
             (:D, 50000), (:D, 100000),
             (:E, 50000), (:E, 100000)]
    seeds = 0:99

    OACCEL2017.runmany(tests,seeds, testset;
                       savejld = savejld, savecsv = savecsv,
                       savebase = savebase)
end

if run_tests_cutest
    min_var = 0
    max_var = 1000
    cutestnames = sort(CUTEst.select(contype=:unc, max_var=max_var,
                                     min_var=min_var))
    oruns = createmeasures(cutestnames, testset)
    mdf = createmeasuredataframe(oruns)

    if savejld
        save(savebase*"/cutest_$(min_var)_$(max_var).jld", "mdf", mdf)
    end
    if savecsv
        CSV.write(savebase*"/cutest_$(min_var)_$(max_var).csv", mdf)
    end
end


text:
using Pkg
Pkg.activate(@__DIR__)

using Documenter, Agents
using Science_ABM
using Literate

# %% Literate convertion
# indir = joinpath(@__DIR__, "..", "examples")
outdir = joinpath(@__DIR__, "src")
# for file in ("example1.jl", "example2.jl")
# 	Literate.markdown(joinpath(indir, file), outdir)
# end

# %%
makedocs(modules = [Science_ABM,],
sitename= "Science_ABM.jl",
authors = "Ali R. Vahdati and Carlos Melian.",
doctest = false,
format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    ),
pages = [
	"Introduction" => "index.md",
	"Model description" => "model_description.md",
	# "Tutorial" => "tutorial.md",
	# "Examples" => [
	#   "Simple Wright-Fisher" => "example1.md",
	#   "Predator prey" => "example2.md"
	# 	],
    ],
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(repo = "github.com/kavir1698/Science_ABM.jl.git",
               target = "build")
end


text:
"""
mutable struct CellInfo
    regime::String
    ispdf::Bool
end
"""

function rc!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{ControlVolume2D2F,2}}

    if KS.set.interpOrder == 1
        return
    end

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    #--- conservative variables ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        #swL = extract_last(ctr[1, j].sw, 1, mode = :view)
        #reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j]))

        swR = KitBase.extract_last(ctr[nx, j].sw, 1, mode = :view)
        reconstruct2!(
            swR,
            ctr[nx-1, j].w,
            ctr[nx, j].w,
            0.5 * (dx[nx-1, j] + dx[nx, j]),
        )
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sw = KitBase.extract_last(ctr[i, j].sw, 1, mode = :view)
            reconstruct3!(
                sw,
                ctr[i-1, j].w,
                ctr[i, j].w,
                ctr[i+1, j].w,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 1:ny
        for i = 1:nx
            sw = KitBase.extract_last(ctr[i, j].sw, 2, mode = :view)
            reconstruct3!(
                sw,
                ctr[i, j-1].w,
                ctr[i, j].w,
                ctr[i, j+1].w,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

    #--- particle distribution function ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        #shL = extract_last(ctr[1, j].sh, 1, mode = :view)
        #reconstruct2!(shL, ctr[1, j].h, ctr[2, j].h, 0.5 * (dx[1, j] + dx[2, j]))
        #sbL = extract_last(ctr[1, j].sb, 1, mode = :view)
        #reconstruct2!(sbL, ctr[1, j].b, ctr[2, j].b, 0.5 * (dx[1, j] + dx[2, j]))

        shR = KitBase.extract_last(ctr[nx, j].sh, 1, mode = :view)
        reconstruct2!(
            shR,
            ctr[nx-1, j].h,
            ctr[nx, j].h,
            0.5 * (dx[nx-1, j] + dx[nx, j]),
        )
        sbR = KitBase.extract_last(ctr[nx, j].sb, 1, mode = :view)
        reconstruct2!(
            sbR,
            ctr[nx-1, j].b,
            ctr[nx, j].b,
            0.5 * (dx[nx-1, j] + dx[nx, j]),
        )
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sh = KitBase.extract_last(ctr[i, j].sh, 1, mode = :view)
            reconstruct3!(
                sh,
                ctr[i-1, j].h,
                ctr[i, j].h,
                ctr[i+1, j].h,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )

            sb = KitBase.extract_last(ctr[i, j].sb, 1, mode = :view)
            reconstruct3!(
                sb,
                ctr[i-1, j].b,
                ctr[i, j].b,
                ctr[i+1, j].b,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 1:ny
        for i = 1:nx
            sh = KitBase.extract_last(ctr[i, j].sh, 2, mode = :view)
            reconstruct3!(
                sh,
                ctr[i, j-1].h,
                ctr[i, j].h,
                ctr[i, j+1].h,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )

            sb = KitBase.extract_last(ctr[i, j].sb, 2, mode = :view)
            reconstruct3!(
                sb,
                ctr[i, j-1].b,
                ctr[i, j].b,
                ctr[i, j+1].b,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
# Autogenerated wrapper script for OpenGLMathematics_jll for any
export glm_hpp

JLLWrappers.@generate_wrapper_header("OpenGLMathematics")
JLLWrappers.@declare_file_product(glm_hpp)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        glm_hpp,
        "include/glm/glm.hpp",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module PrintOpTrees

using FillArrays
using Gridap.Arrays
using Gridap.Geometry
using Gridap.CellData
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Arrays: CompressedArray
using Gridap.Visualization

io = IOBuffer()

a = zeros(3)
print_op_tree(io,a)

b = Fill(4,10)
print_op_tree(io,b)

c = CompressedArray([1.0,2.0],[1,2,1,1,1,2,2,1])
print_op_tree(io,c,showid=true)

d = apply(-,collect(c))
d = apply(+,d,d)
print_op_tree(io,d)

r = reindex(d,[4,3,3,1])
print_op_tree(io,r,showid=true)

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

trian = Triangulation(model)
quad = CellQuadrature(trian,2)

u(x) = x[1]
cf = CellField(u, trian)

xe = get_cell_coordinates(trian)
q = get_coordinates(quad)

uq = evaluate(cf,q)

V = FESpace(model=model,valuetype=Float64,reffe=:Lagrangian,order=1,conformity=:H1)

vh = FEFunction(V,rand(num_free_dofs(V)))

vh_q = evaluate(vh,q)
∇vh_q = evaluate(∇(vh),q)

print_op_tree(io,vh_q)
print_op_tree(io,∇vh_q)

print_op_tree(io,xe)
print_op_tree(io,q)

print_op_tree(io,uq)

data = read(seekstart(io), String)

end # module


text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

### Data structures that are additionally available for media <: PartialPureSubstance -----------------------

const X_dummy = fill(0.0,0)

"""
    state = setState_pT(medium, p,T)

Generate a state object for medium `medium::PureSubstance` for
pressure `p` [Pa] and temperature `T` [K].
"""
setState_pT( m::PureSubstance,p,T) = setState_pTX(m,p,T,X_dummy)

"""
    setState_pT!(state, p,T)

Update the `state::ThermodynamicState` of a `PureSubstance` medium with
pressure `p` [Pa] and temperature `T` [K].
"""
setState_pT!(state::ThermodynamicState,p,T) = setState_pTX!(state,p,T,X_dummy)

"""
    state = setState_ph(medium, p,h)

Generate a state object for medium `medium::PureSubstance` for
pressure `p` [Pa] and specific enthalpy `h` [J/kg].
"""
setState_ph(m::PureSubstance,p,h) = setState_phX(m,p,h,X_dummy)

"""
    setState_ph!(state, p,h)

Update the `state::ThermodynamicState` of a `PureSubstance` medium with
pressure `p` [Pa] and specific enthalpy `h` [J/kg]].
"""
setState_ph!(state::ThermodynamicState,p,h) = setState_phX!(state,p,h,X_dummy)

"""
    state = setState_ps(medium, p,s)

Generate a state object for medium `medium::PureSubstance` for
pressure `p` [Pa] and specific entropy `s` [J/(kg*K)].
"""
setState_ps( m::PureSubstance,p,s) = setState_psX(m,p,s,X_dummy)

"""
    setState_ps!(state, p,s)

Update the `state::ThermodynamicState` of a `PureSubstance` medium with
pressure `p` [Pa] and specific entropy `s` [J/(kg*K)].
"""
setState_ps!(state::ThermodynamicState,p,s) = setState_psX!(state,p,s,X_dummy)

"""
    state = setState_dT(medium, d,T)

Generate a state object for medium `medium::PureSubstance` for
density `d` [kg/m^3] and temperature `T` [K].
"""
setState_dT( m::PureSubstance,d,T) = setState_dTX(m,d,T,X_dummy)

"""
    setState_dT!(state, d,T)

Update the `state::ThermodynamicState` of a `PureSubstance` medium with
density `d` [kg/m^3] and temperature `T` [K].
"""
setState_dT!(state::ThermodynamicState,d,T) = setState_dTX!(state,d,T,X_dummy)

### Medium functions that are additionally available for media <: PartialPureSubstance -----------------------

"""
    density_ph(medium,p,h)

Return density in [kg/m^3] for `medium::PureSubstance` from
pressure `p` [Pa] and specific enthalpy `h` [J/kg].
"""
density_ph(m::PureSubstance, p,h) = density(setState_ph(m,p,h))

"""
    temperature_ph(medium,p,h)

Return temperature in [K] for `medium::PureSubstance` from
pressure `p` [Pa] and specific enthalpy `h` [J/kg].
"""
temperature_ph(m::PureSubstance, p,h) = temperature(setState_ph(m,p,h))

"""
    pressure_dT(medium,d,T)

Return pressure in [Pa] for `medium::PureSubstance` from density `d` in [kg/m^3]
and temperature `T` [K].
"""
pressure_dT(m::PureSubstance, d,T) = pressure(setState_dT(m,d,T))

"""
    specificEnthalpy_dT(medium,d,T)

Return specific enthalpy in [J/kg] for `medium::PureSubstance` from
density `d` [kg/m^3] and and temperature `T` [K].
"""
specificEnthalpy_dT(m::PureSubstance, d,T) = specificEnthalpy(setState_dT(m,d,T))

### Derivatives needed for mass and energy balance ------------------------------------------------------

density_ph_der_1(m::PureSubstance, p, h) = 0.0
density_ph_der_2(m::PureSubstance, p, h) = undefinedFunction("density_ph_der_2", m)
density_ph_der_3(m::PureSubstance, p, h) = undefinedFunction("density_ph_der_3", m)


text:
#

# ### **MAE 103 - Elementary Fluid Mechanics**
#


text:
# This file is a part of EasyPkg, licensed under the MIT License (MIT).

export include_all_sources

"""
`include_all_sources([basedir])`

Recursively includes all sources in the specified directory and it's
subdirectories, by searching for and including all ".jl" files. Directories
are processed before files.

`basedir` defaults to the directory containing the current source file
(or `pwd()`, if run from the REPL), and is interpreted relative to it
if `basedir` is not an absolute path.

`include_all_sources()` reduces the average Julia project's "Project.jl" to

```julia
import EasyPkg
EasyPkg.include_all_sources()

function __init__()
    ...
end
```

For this to work, the structure of the sources files must not depend on
a specific inclusion order - so the source files should not contain any
initialization steps that depend on another source file. For precompiled
packages, this shouldn't be an issue, as all initialization steps will
reside in `__init()__` anyhow.
"""

function include_all_sources(basedir::AbstractString = include_base_path())
    const abs_basedir = isabspath(basedir) ? basedir : joinpath(include_base_path(), basedir)
    const sources = find_sources(abs_basedir)

    for fname in sources
        const absfname = joinpath(abs_basedir, fname)
        include(absfname)
    end
end


text:
using Plots

# Some plots utils, multiple dispatch is used.
function plotSets(title, set1, labels)
    x1 = Vector{Float64}()
    y1 = Vector{Float64}()
    for el in set1
        push!(x1, el[1])
        push!(y1, el[2])
    end
return scatter([x1], [y1], title=title, label = labels[1], legend=:outertopright, reuse = false)
end

function plotSets(title, set1, set2, labels)
    x1 = Vector{Float64}()
    y1 = Vector{Float64}()
    for el in set1
        push!(x1, el[1])
        push!(y1, el[2])
    end

    x2 = Vector{Float64}()
    y2 = Vector{Float64}()
    for el in set2
        push!(x2, el[1])
        push!(y2, el[2])
    end

    return scatter([x1,x2], [y1,y2], title=title, marker=[:o :square], c=[:skyblue2 :orange], label = [labels[1] labels[2]],legend=:outertopright, reuse = false)
end

function plotSets(title, set1, set2, set3, labels)
    x1 = Vector{Float64}()
    y1 = Vector{Float64}()
    for el in set1
        push!(x1, el[1])
        push!(y1, el[2])
    end

    x2 = Vector{Float64}()
    y2 = Vector{Float64}()
    for el in set2
        push!(x2, el[1])
        push!(y2, el[2])
    end

    x3 = Vector{Float64}()
    y3 = Vector{Float64}()
    for el in set3
        push!(x3, el[1])
        push!(y3, el[2])
    end

    return scatter([x1,x2,x3], [y1,y2,y3], title = title, marker=[:o :o :square], c=[:lightgray :orange :green], label = [labels[1] labels[2] labels[3]], legend=:outertopright, reuse = false)
end


text:
# This file is a part of HistTools.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    HistTools

Bayesian Blocks algorithm described in [Scargle 2021](https://doi.org/10.1088/0004-637X/764/2/167).
"""
module HistTools

using StatsBase

import StatsBase: reshape

include("core.jl")

end # module


text:
using SafeTestsets

##############################################
# - Run unit tests
##############################################

# Each function has its own .jl file containing a SafeTestset
# - Using the 'SafeTestsets' approach requires each test set to run as a
#   stand alone file (e.g. contains its own Using ... packages code for each function).
# - SafeTestsets runs each test set in its own clean environment - no leakage between
#   environments and test sets.

# (1) Test 'plot_qq' function
@safetestset "test plot_qq" begin include("test_plot_qq.jl") end

# (2) Test 'plot_fit_res' function
@safetestset "test plot_fit_res" begin include("test_plot_fit_res.jl") end

# (3) Test 'plot_scale_loc' function
@safetestset "test plot_scale_loc" begin include("test_plot_scale_loc.jl") end

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
#include("CVRPInstanceModule.jl")
#include("MyGraph.jl")
#include("DagFlowInstance.jl")

module CVRPGraph

using ..CVRPInstanceModule
using ..MyGraphModule
using ..DagFlowInstanceModule

function constructCVRPGraph(instance::CVRPInstance)
    #println(instance.c)
    n = instance.n
    nn = instance.nn
    Q = instance.Q
    demand = instance.demand
    println("n=$n, Q=$Q")
    # nodeDemandMap[i,d] should contain the id of the node in the extended graph corresponding to visiting node i in the original
    # graph on a route with total demand d
    nodeDemandMap = zeros(Int32, n, Q)
    for i=1:n
        if demand[i] <= 0
            println("Node $i has zero or negative demand. We cannot construct extended CVRP graph")
            throw("Cannot construct extended CVRP graph. Demand too low")
        end
        if demand[i] > Q
            println("Node $i has a demand that is higher than the vehicle capacity. We cannot construct extended CVRP graph")
            throw("Cannot construct extended CVRP graph. Demand too high")
        end
    end
    graph=MyGraphOut()
    startNode = addNode!(graph)
    endNode = addNode!(graph)
    for i=1:n
        node = recursiveConstruct(convert(Int32, i), demand[i], graph, nodeDemandMap, instance, endNode)
        addArcNoDup!(graph,startNode,node, instance.c[nn,i])
    end
    nNodes, nArcs = MyGraphModule.size(graph)
    println("Extended graph contains $nNodes nodes and $nArcs arcs.")

    nodeInfoDict = Dict{Int32,String}()
    nodePartition = []
    node2Partion = zeros(Float64,nNodes)
    for i=1:n
        nodeSet = Vector{Int32}()
        partitionId = length(nodePartition)+1
        for q=1:Q
            nodeId = nodeDemandMap[i,q]
            if nodeId != 0
                push!(nodeSet, nodeId)
                push!(nodeInfoDict, nodeId => "node $i acc demand $q")
                node2Partion[nodeId] = partitionId
            end
        end
        push!(nodePartition, nodeSet)
    end
    push!(nodeInfoDict, startNode => "Start")
    push!(nodeInfoDict, endNode => "End")

    # Create DagFlowInstance
    gIO = MyGraphModule.createMyGraph(graph)
    # right now the CVRP instance does not contain a max number of vehicles so we just set it to n.
    maxVeh = n

    nodeOrder=Vector{Int32}()
    push!(nodeOrder,startNode)
    for q=1:Q
        for i=1:n
            if nodeDemandMap[i,q] != 0
                push!(nodeOrder, nodeDemandMap[i,q])
            end
        end
    end

    dfi = DagFlowInstance(gIO, nodePartition, node2Partion, startNode, endNode, nodeOrder, [OutFlowLimit(startNode, maxVeh)], Set(startNode), nodeInfoDict)
    return dfi
end

function recursiveConstruct(node::Int32, totalDemand::Int32, graph::MyGraphOut, nodeDemandMap, instance::CVRPInstance, endNode::Int32)
    n = instance.n
    nn = instance.nn
    Q = instance.Q
    demand = instance.demand
    newNode = addNode!(graph)
    nodeDemandMap[node,totalDemand] = newNode
    # connect new node to the end node:
    addArcNoDup!(graph, newNode, endNode, instance.c[node,nn])
    for i=1:n
        if i != node
            newDemand = totalDemand + demand[i]
            if newDemand <= Q
                otherNode = nodeDemandMap[i,newDemand]
                # have we already constructed this extended node?
                if otherNode == 0
                    # no, we need to do that
                    otherNode = recursiveConstruct(convert(Int32, i),newDemand, graph, nodeDemandMap, instance, endNode)
                end
                # connect extended nodes
                addArcNoDup!(graph, newNode, otherNode, instance.c[node,i])
            end
        end
    end
    return newNode
end

#instance = CVRPInstanceModule.readInstance("Data/CVRP/Data/A-n32-k5.vrp", true)
#instance = CVRPInstanceModule.readInstance("Data/CVRP/Data/M-n200-k16.vrp", true)
#graph = constructCVRPGraph(instance)

end


text:
function digital_root(a::Integer)::Int
    while a >= 10
        a = sum(digits(a))
    end
    a
end

function digital_root(a::Vector{<:Integer})::Int
    a = cumprod(big.(a))
    digital_root(sum(a))
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
module TreesHeaps
using DataStructures, AbstractTrees, Reexport
import Base:pairs, insert!, delete!, fieldnames, getproperty, ==, findmin, findmax
import AbstractTrees: printnode, print_tree, PostOrderDFS, PreOrderDFS, TreeIterator

export
    # Abstract types
    AbstractTree, AbstractNode, AbstractBinaryNode,

    # Nodes
    NullNode, SimpleBinaryNode, SBN, HeightBinaryNode, HBN, RedBlackBinaryNode, RBN,

    # Trees
    BinarySearchTree, BST, AVLTree, AVL, SplayTree, Splay, RedBlackTree, RBT,

    # preoperty
    height, isnull,
    # not imported: eltype, size and length

    # operations
    search, splay!, topdownsplay!, insert!, topdowninsert!,
    delete!, topdowndelete!, findmin, findmax, findrightmin, findleftmax,

    # internal operations
    split!, topdownsplit!, join!,

    # initialize plotting
    plot_init

@reexport using AbstractTrees

include("Nodetree.jl")
include("Core.jl")
include("Balance.jl")
include("Interface.jl")

# Avoid importing Plots as default
plot_init() = include(joinpath(@__DIR__(), "Plot.jl"))

__init__() = precompile(plot_init,())

end # module


text:
#=
NOTE: This module is still incomplete, do not use in simualtions.

References:

Mark L. Henle1 and Alex J. Levine. Hydrodynamics in curved membranes: The effect of geometry on particulate mobility.
=#

include("../geometry/sphere.jl")
include("./fluid.jl")
include("../utils/data.jl")

using Symbolics
using DataFrames
using CSV
using Interpolations

#=
(θ_0,ϕ_0) is the location of the force
(θ,ϕ) is the location of the observation point
f is the force magnitude
α is the orientation angle of force from ∂ϕ_0
β is the angle from ∂ϕ along which we measure the component of velocity
=#
@variables θ_0 ϕ_0 θ ϕ f α β dl

∂θ = Differential(θ)
∂ϕ = Differential(ϕ)
∂θ_0 = Differential(θ_0)
∂ϕ_0 = Differential(ϕ_0)

grad_θϕ( f ) = [ ∂θ(f) csc(θ)*∂ϕ(f) ]'
grad_θ_0ϕ_0( f ) = [ ∂θ_0(f) csc(θ_0)*∂ϕ_0(f) ]'

γ = acos( cos(θ)*cos(θ_0) + sin(θ)*sin(θ_0)*cos(ϕ - ϕ_0) )
∂γ = Differential(γ)

twist = [ 0   1; -1  0 ]
tgrad_θϕ( f ) = twist * grad_θϕ(f)
tgrad_θ_0ϕ_0( f ) = twist * grad_θ_0ϕ_0(f)

#=
α is the angle that the force at (θ_0,ϕ_0) makes with ∂/∂ϕ_0 vector at the point
γ is the geodesic length between (θ,ϕ) and (θ_0,ϕ_0)
β is the angle with ∂/∂ϕ vector along which we measure the component of velocity

v_t = \gradient_{θ,ϕ} ( force_t ⋅ \gradient_{θ_0,ϕ_0} P(γ) )
    = d/dγ( P(γ) ) \gradient_{θ,ϕ} ( force_t ⋅ \gradient_{θ_0,ϕ_0} γ )
v is the velocity vector at point (θ,ϕ) due to force at point (θ',ϕ')
force_t is force rotated by π/2
i.e. force_t' = force' ⋅ [ ( cos(π/2), sin(π/2) ), ( - sin(π/2), cos(π/2) ) ]
              = force' ⋅ [ ( 0, 1 ), ( - 1, 0 ) ]
v_t is v rotated by π/2
i.e. v = [ ( 0, 1 )', ( - 1, 0 )' ] ⋅ v_t

v_t = d/dγ( P(γ) ) \gradient_{θ,ϕ} ( force_t ⋅ (∂/∂θ_0(γ), cosec(θ_0) ∂/∂ϕ_0(γ) )
    = d/dγ( P(γ) ) \gradient_{θ,ϕ} ( force' ⋅ [ ( 0, 1 ), ( - 1, 0 ) ] ⋅ (∂/∂θ_0(γ), cosec(θ_0) ∂/∂ϕ_0(γ) )
    = d/dγ( P(γ) ) \gradient_{θ,ϕ} ( force' ⋅ ( cosec(θ_0) ∂/∂ϕ_0(γ), -∂/∂θ_0(γ) ) )
    = d/dγ( P(γ) ) \gradient_{θ,ϕ} ( force_θ_0 cosec(θ_0) ∂/∂ϕ_0(γ) - force_ϕ_0 ∂/∂θ_0(γ) )
    = d/dγ( P(γ) ) ( force_θ_0 cosec(θ_0) ∂/∂θ( ∂/∂ϕ_0(γ) ) - force_ϕ_0  ∂/∂θ( ∂/∂θ_0(γ) ),
                       force_θ_0 cosec(θ) cosec(θ_0) ∂/∂ϕ( ∂/∂ϕ_0(γ) ) - force_ϕ_0 cosec(θ) ∂/∂ϕ( ∂/∂θ_0(γ) )
                     )

v = d/dγ( P(γ) ) ( force_θ_0 cosec(θ) cosec(θ_0) ∂/∂ϕ( ∂/∂ϕ_0(γ) ) - force_ϕ_0 cosec(θ) ∂/∂ϕ( ∂/∂θ_0(γ),
                   - force_θ_0 cosec(θ_0) ∂/∂θ( ∂/∂ϕ_0(γ) ) + force_ϕ_0 ∂/∂θ( ∂/∂θ_0(γ) ) )
                 )
  = d/dγ( P(γ) ) ( f sin(α) cosec(θ) cosec(θ_0) ∂/∂ϕ( ∂/∂ϕ_0(γ) ) - f cos(α)  cosec(θ) ∂/∂ϕ( ∂/∂θ_0(γ),
                   - f sin(α) cosec(θ_0) ∂/∂θ( ∂/∂ϕ_0(γ) ) + f cos(α) ∂/∂θ( ∂/∂θ_0(γ) ) )
                 )

 v_th = ∂γ( P ) * (   f * sin(α) * csc(θ) * csc(θ_0) * ∂ϕ( ∂ϕ_0(γ) )  -  f * cos(α) * csc(θ) * ∂ϕ( ∂θ_0(γ) ) )
 v_ph = ∂γ( P ) * ( - f * sin(α)          * csc(θ_0) * ∂θ( ∂ϕ_0(γ) )  +  f * cos(α) *          ∂θ( ∂θ_0(γ) ) )
=#

# P as a function of γ is stored in file
pData = DataFrame(CSV.File("./data/pointforce_import_lc.tsv"))
# we compute the first few derivatives and store them in global variables to avoid
# recomputing the derivatives each time are encountered in final expression
pDataInterpolated = LinearInterpolation( pData.x, pData.y; extrapolation_bc = Line())
pDataInterpolated1 = dataDerivative( ( pDataInterpolated, pData.x ) )[1]
pDataInterpolated2 = dataDerivative( ( pDataInterpolated1, pData.x ) )[1]
pDataInterpolated3 = dataDerivative( ( pDataInterpolated2, pData.x ) )[1]
pDataInterpolated4 = dataDerivative( ( pDataInterpolated3, pData.x ) )[1]
pDataInterpolated5 = dataDerivative( ( pDataInterpolated4, pData.x ) )[1]

function Pf(n::Int64, x::Float64)
  if n==0
    pDataInterpolated(x)
  elseif n==1
    pDataInterpolated1(x)
  elseif n==2
    pDataInterpolated2(x)
  elseif n==3
    pDataInterpolated3(x)
  elseif n==4
    pDataInterpolated4(x)
  elseif n==5
    pDataInterpolated5(x)
  else
    throw( sprint(n) + "-derivative of P not precomputed" )
  end
end

# register the pre-computed derivates in the symbolics framework
Symbolics.derivative( ::typeof(Pf), args::NTuple{2,Any}, ::Val{2} ) = SymbolicUtils.Term( Pf, [ args[1]+1, args[2] ] )
Symbolics.@register Pf( data::Int64, x )

P = Pf( 0, γ )

v =  tgrad_θϕ( ( [ f*sin(α) f*cos(α) ] * tgrad_θ_0ϕ_0( P ) )[ 1 ] )
# v = ∂γ( dP ) * tgrad_θϕ( ( [ f*sin(α) f*cos(α) ] * tgrad_θ_0ϕ_0( γ ) )[ 1 ] ) +
#     ∂γ(∂γ( P )) * tgrad_θϕ( γ ) * ( [ f*sin(α) f*cos(α) ] * tgrad_θ_0ϕ_0( γ ) )[ 1 ]

# gradient matrix of v at (θ,ϕ)
#vg = [ ∂θ(v_θ) csc(θ)∂ϕ(v_θ); ∂θ(v_ϕ) csc(θ)∂ϕ(v_ϕ)]
vg = [ ∂θ(v[1]) csc(θ)*∂ϕ(v[1]); ∂θ(v[2]) csc(θ)*∂ϕ(v[2]) ]

v_simplified = simplify.( expand_derivatives.( v ) )
vg_simplified = simplify.( expand_derivatives.( vg ) )

pointForceV = build_function( v_simplified, [θ_0, ϕ_0, θ, ϕ, f, α ], expression=Val{false} )
pointForceVG = build_function( vg_simplified, [θ_0, ϕ_0, θ, ϕ, f, α ], expression=Val{false} )

#========================================
# Interface functions
========================================#

struct PointForce <: Object
    location :: SphericalCoord  # (θ_0,ϕ_0) location of
    orientation :: Float64  # angle that the force makes with ∂/∂ϕ_0
    magnitude :: Float64  # magnitude of the force
end

# velocity field generated by point force
function velocityField( destination::SphericalCoord, source::PointForce, fluid::FluidParams )::SphericalCoord
    θs = source.location.theta
    ϕs = source.location.phi
    θd = destination.theta
    ϕd = destination.phi
    if θs == θd && ϕs == ϕd
        return( SphericalCoord( 0, 0 ) )
    else
        velocity = pointForceV[1]([ θs, ϕs, θd, ϕd, source.magnitude, source.orientation ])
        return( SphericalCoord( velocity[ 1 ], velocity[ 2 ] ) )
    end
end

# gradient of velocity field generated by dipole force
# NOTE: not yet implemented
function velocityCurlField( destination::SphericalCoord, source::PointForce, fluid::FluidParams )::Float64
    return( 0 )
end

# gradient of velocity field generated by point force
function velocityGradientField( destination::SphericalCoord, source::PointForce, fluid::FluidParams )::Array{Float64,2}
    θs = source.location.theta
    ϕs = source.location.phi
    θd = destination.theta
    ϕd = destination.phi
    if θs == θd && ϕs == ϕd
        return( [0 0; 0 0] )
    else
        velocityGrad = pointForceVG[1]([ θs, ϕs, θd, ϕd, source.magnitude, source.orientation ])
        return( velocityGrad )
    end
end

# how point force moves in the presence of external velocity field
# NOTE: not yet implemented
function stepForward( object::PointForce, v::SphericalCoord, vc::Float64, dt::Float64 )::PointForce
    # right now we don't change the location of pointforce based on external field
    return( object )
end

#plot the point forces on 3d chart
function plot3d!( threeDPlot::Union{Scene,LScene}, pointForce::Observable{PointForce} )::Nothing
    pointForcePosState = @lift( $pointForce.location )
    pointForceVectorMap( f, α ) = SphericalCoord( f*sin(α), f*cos(α) )
    pointForceForceState = @lift( pointForceVectorMap(  $pointForce.magnitude,  $pointForce.orientation ) )
    pfps = @lift( sphericalToCartesian( $pointForcePosState ) )
    pfpx = @lift( [$pfps.x] )
    pfpy = @lift( [$pfps.y] )
    pfpz = @lift( [$pfps.z] )
    pffs = @lift( sphericalToCartesianVelocity( $pointForcePosState, $pointForceForceState ) )
    pffx = @lift( [$pffs.x] )
    pffy = @lift( [$pffs.y] )
    pffz = @lift( [$pffs.z] )
    arrows!(threeDPlot, pfpx, pfpy, pfpz, pffx, pffy, pffz, arrowsize = 0.05, lengthscale = 0.05 )
    nothing
end

# show point forces on theta phi chart
function plot2d!( twoDPlot::Union{Scene,LScene}, pointForce::Observable{PointForce} )::Nothing
    pointForcePosState = @lift( $pointForce.location )
    pointForceVectorMap( f, α ) = SphericalCoord( f*sin(α), f*cos(α) )
    pointForceForceState = @lift( pointForceVectorMap( $pointForce.magnitude,  $pointForce.orientation ) )
    pfptx = @lift( [$pointForcePosState.phi] )
    pfpty = @lift( [$pointForcePosState.theta] )
    pfftx = @lift( [$pointForceForceState.phi * csc($pointForcePosState.theta)] )
    pffty = @lift( [$pointForceForceState.theta] )
    arrows!( twoDPlot, pfptx, pfpty, pfftx, pffty, arrowsize = 0.05, lengthscale = 0.05 )
    nothing
end

#========================================
# Tests
========================================#

#=
Rotational symmetry check

If we have a point force at ( θ0, ϕ0 ) making an angle α with ϕ axis,
then the look at component of velocity induced at point ( θ, ϕ )  making an angle β with ϕ axis

This should be same as component of velocity induced at point ( θ0, ϕ0 ) making an angle β - rot with ϕ axis
by a point force at ( θ, ϕ ) making an angle α + rot

where rot is π + the rotation a vector suffers when parallel transported from ( θ0, ϕ0 ) to ( θ, ϕ )

the operations involved in this symmetry are exchange force and velocity locations,
parallely transporting the force and velocity vectors along the geodesic connecting them,
and then rotating the frame 180 degrees about the diameter of sphere passing through midpoint of geodesic
=#

v_th = v[ 1 ]
v_ph = v[ 2 ]
v_beta = expand_derivatives( v_ph * cos(β) + v_th * sin(β) )

tangent_angle = expand_derivatives( atan( ∂θ(γ), csc(θ)∂ϕ(γ) ) )
tangent_0_angle = expand_derivatives( atan( ∂θ_0(γ), csc(θ_0)∂ϕ_0(γ) ) )
rot = tangent_angle - tangent_0_angle

args = Dict(
      θ => π/2.5
    , ϕ => 0.4
    , θ_0 => 2π/3
    , ϕ_0 => -π/4
    , f => 1
    , α => 0.2
    , β => 0.4
)

exchange = Dict( θ_0 => θ, ϕ_0 => ϕ, θ => θ_0, ϕ => ϕ_0, α => α + rot, β => β - rot )
v_ex = substitute( v_beta, exchange )
sym_check = simplify( v_ex - v_beta )  # this should be zero, atleast in domain of interest
sym_check_eval = substitute( sym_check, args ) # this should be zero for all args as sym_check should be zero
@assert abs(sym_check_eval) < 1e-10

tangent_ex_angle = substitute( tangent_angle, exchange )
tangent_check = simplify( tangent_ex_angle - tangent_0_angle )
tangent_check_eval = substitute( tangent_check, args ) # this should be zero
@assert abs(tangent_check_eval) < 1e-10


text:
immutable Skellam <: DiscreteUnivariateDistribution
    μ1::Float64
    μ2::Float64

    function Skellam(μ1::Float64, μ2::Float64)
        μ1 > 0.0 && μ2 > 0.0 || error("μ1 and μ2 must be positive.")
        new(μ1, μ2)
    end

    @compat Skellam(μ1::Real, μ2::Real) = Skellam(Float64(μ1), Float64(μ2))

    Skellam(μ::Real) = Skellam(μ, μ)

    Skellam() = new(1.0, 1.0)
end

@distr_support Skellam -Inf Inf

### Parameters

params(d::Skellam) = (d.μ1, d.μ2)

### Statistics

mean(d::Skellam) = d.μ1 - d.μ2

var(d::Skellam) = d.μ1 + d.μ2

skewness(d::Skellam) = mean(d) / (var(d)^1.5)

kurtosis(d::Skellam) = 1.0 / var(d)

### Evaluation

function logpdf(d::Skellam, x::Int)
    μ1, μ2 = params(d)
    - (μ1 + μ2) + (x / 2.0) * log(μ1 / μ2) + log(besseli(x, 2.0 * sqrt(μ1) * sqrt(μ2)))
end

pdf(d::Skellam, x::Int) = exp(logpdf(d, x))

function mgf(d::Skellam, t::Real)
    μ1, μ2 = params(d)
    exp(-(μ1 + μ2) + μ1 * exp(t) + μ2 * exp(-t))
end

function cf(d::Skellam, t::Real)
    μ1, μ2 = params(d)
    exp(-(μ1 + μ2) + μ1 * cis(t) + μ2 * cis(-t))
end

### Sampling

rand(d::Skellam) = rand(Poisson(d.μ1)) - rand(Poisson(d.μ2))


text:
@inline Base.signbit(num::OraNumber) = isnegative(num)

@inline function Base.exponent(num::OraNumber) :: Int
    return Int(decode_exponent_byte(num))
end

function Base.show(io::IO, num::OraNumber)
    print(io, "OraNumber(\"")
    print_number_string(io, num)
    print(io, "\")")
end

function Base.string(num::OraNumber)
    io = IOBuffer()
    print_number_string(io, num)
    return String(take!(io))
end

@inline Base.zero(::Type{OraNumber}) = OraNumber(0x01, 0x80, (0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
@inline Base.one(::Type{OraNumber}) = OraNumber(0x02, 0xc1, (0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))

@inline OraNumber(::Nothing) = OraNumber(0xff, 0x00, (0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))

function Base.:(==)(num1::OraNumber, num2::OraNumber)
    @assert isnormalized(num1) && isnormalized(num2) "Can't compare non-normalized OraNumbers."

    if num1.len != num2.len || num1.ex != num2.ex
        return false
    end

    # compare mantissa
    for i in 1:sizeof_mantissa(num1)
        if num1.mantissa[i] != num2.mantissa[i]
            return false
        end
    end

    return true
end


text:
"""
    Controller{}

Controller is an abstract type that represents a controller acting on Task state.
"""
abstract type Controller{} end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
module EAGOBranchBound

using IntervalArithmetic: Interval, diam
using EAGOIntervalArithmetic: MCInterval, diam

import IntervalRootFinding.bisect

export BnBSolver, BnBModel, set_Branch_Scheme!, set_Bisect_Func!,
       set_Verbosity!, set_to_default!, solveBnB!, getsolution, getobjval,
       getfeasibility, getobjbound, LBDtime, UBDtime, default_pre, default_post

include("src/types/BnBModel.jl")
include("src/types/BnBSolver.jl")
include("src/schemes/Node_Select.jl")
include("src/schemes/Branch_Method.jl")
include("src/schemes/Bisect_Method.jl")
include("src/utils/Access.jl")
include("src/utils/Display.jl")
include("src/utils/Fathom.jl")
include("src/utils/Processing.jl")
include("src/utils/Checks.jl")
include("src/utils/Set_Options.jl")
include("src/solve.jl")

function __init__()
end

end # module


text:
using Test
import VTUFileHandler

include("test_VTU.jl")

max_stars_repo_name: baxmittens/VTUFileHandler

max_issues_repo_name: baxmittens/VTUFileHandler

max_forks_repo_name: baxmittens/VTUFileHandler


text:
function align_fst!(fst::FST, opts::Options)
    is_leaf(fst) && return
    const_idxs = Int[]
    assignment_idxs = Int[]
    pair_arrow_idxs = Int[]

    for (i, n) in enumerate(fst.nodes)
        if is_leaf(n)
            continue
        elseif opts.align_struct_field && (n.typ === Struct || n.typ === Mutable)
            align_struct!(n)
        elseif opts.align_conditional && n.typ === Conditional
            align_conditional!(n)
        else
            align_fst!(n, opts)
        end

        if is_assignment(n) || n.typ === Kw
            # Gather all assignments within the current code block
            # they will be aligned at the end
            push!(assignment_idxs, i)
        elseif n.typ === Binary && op_kind(n) === Tokens.PAIR_ARROW
            push!(pair_arrow_idxs, i)
        end
    end

    if opts.align_assignment
        align_binaryopcalls!(fst, assignment_idxs)
    end

    if opts.align_pair_arrow
        align_binaryopcalls!(fst, pair_arrow_idxs)
    end
end

"""
Group of FST node indices and required metadata to potentially align them.

- `node_idxs`. Indices of FST nodes affected by alignment.
- `line_offsets`. Line offset of the character nodes may be aligned to
in the source file.
- `lens`. Length of the FST node prior to the alignment character. Used
to calculate extra whitespace padding.
- `whitespaces`. Number of whitespaces between the alignment character and
the prior FST node. If this is > 1 it signifies additional whitespace was
manually added by the user since the formatter would only use 0 or 1 whitespaces.

"""
struct AlignGroup
    node_idxs::Vector{Int}
    line_offsets::Vector{Int}
    lens::Vector{Int}
    whitespaces::Vector{Int}
end
AlignGroup() = AlignGroup(Int[], Int[], Int[], Int[])

function Base.push!(g::AlignGroup, idx::Int, line_offset::Int, len::Int, ws::Int)
    push!(g.node_idxs, idx)
    push!(g.line_offsets, line_offset)
    push!(g.lens, len)
    push!(g.whitespaces, ws)
    return
end

function align_to(g::AlignGroup)::Union{Nothing,Int}
    # determine whether alignment might be warranted
    max_len, max_idx = findmax(g.lens)
    max_idxs = findall(==(g.line_offsets[max_idx]), g.line_offsets)
    length(max_idxs) > 1 || return nothing

    # Is there custom whitespace?
    # Formatter would only add 0 or 1 whitespaces.
    # > 2 implies a manual edit in the source file.
    for i in max_idxs
        g.whitespaces[i] > 1 && return max_len
    end

    return nothing
end

function align_binaryopcall!(fst::FST, diff::Int)
    # insert whitespace before and after operator
    fidx = findfirst(x -> x.typ === WHITESPACE, fst.nodes)
    lidx = findlast(x -> x.typ === WHITESPACE, fst.nodes)

    if fidx === nothing
        insert!(fst, 2, Whitespace(diff))
    else
        fst[fidx] = Whitespace(diff)
    end

    if lidx === nothing
        insert!(fst, 4, Whitespace(1))
    end
end

"""
    align_struct!(fst::FST)

Aligns struct fields.
"""
function align_struct!(fst::FST)
    idx = findfirst(n -> n.typ === Block, fst.nodes)
    idx === nothing && return
    length(fst[idx]) == 0 && return

    block_fst = fst[idx]
    prev_endline = block_fst[1].endline
    groups = AlignGroup[]
    g = AlignGroup()

    for (i, n) in enumerate(block_fst.nodes)
        if n.typ === Binary
            if n.startline - prev_endline > 1
                push!(groups, g)
                g = AlignGroup()
            end

            nlen = length(n[1])
            idx = findfirst(x -> x.typ === OPERATOR, n.nodes)
            ws = n[idx].line_offset - (n.line_offset + nlen)

            push!(g, i, n[idx].line_offset, nlen, ws)
            prev_endline = n.endline
        end
    end
    push!(groups, g)

    for g in groups
        align_len = align_to(g)
        align_len === nothing && continue
        for (i, idx) in enumerate(g.node_idxs)
            diff = align_len - g.lens[i] + 1
            align_binaryopcall!(block_fst[idx], diff)
            block_fst[idx].nest_behavior = NeverNest
        end
    end
end

"""
    align_binaryopcalls!(fst::FST, op_idxs::Vector{Int})

Aligns binary operator expressions.

Additionally handles the case where a keyword such as `const` is used
prior to the binary op call.
"""
function align_binaryopcalls!(fst::FST, op_idxs::Vector{Int})
    length(op_idxs) > 1 || return
    prev_endline = fst[op_idxs[1]].endline
    groups = AlignGroup[]
    g = AlignGroup()

    for i in op_idxs
        n = fst[i]
        if n.startline - prev_endline > 1
            push!(groups, g)
            g = AlignGroup()
        end

        binop, nlen, ws = if n.typ === Binary || n.typ === Kw
            nlen = length(n[1])
            n, nlen, (n[3].line_offset - n.line_offset) - nlen
        else
            binop = n[3]
            nlen = length(binop[1]) + length(fst[i][1]) + length(fst[i][2])
            binop, nlen, (binop[3].line_offset - n.line_offset) - nlen
        end

        # @info "" binop.typ nlen ws
        push!(g, i, binop[3].line_offset, nlen, ws)

        prev_endline = n.endline
    end
    push!(groups, g)

    for g in groups
        align_len = align_to(g)
        align_len === nothing && continue

        for (i, idx) in enumerate(g.node_idxs)
            diff = align_len - g.lens[i] + 1

            typ = fst[idx].typ
            if typ === Binary || typ === Kw
                align_binaryopcall!(fst[idx], diff)
            else
                align_binaryopcall!(fst[idx][3], diff)
            end
            fst[idx].nest_behavior = NeverNest
        end
    end

    return
end

"""
    align_conditional!(fst::FST)

Aligns a conditional expression.
"""
function align_conditional!(fst::FST)
    nodes = flatten_conditionalopcall(fst)

    cond_group = AlignGroup()
    cond_prev_endline = 0

    colon_group = AlignGroup()
    colon_prev_endline = 0

    for (i, n) in enumerate(nodes)
        if n.typ === OPERATOR && n.val == "?"
            if cond_prev_endline != n.endline
                nlen = length(nodes[i-2])
                ws = n.line_offset - (nodes[i-2].line_offset + nlen)
                push!(cond_group, i, n.line_offset, nlen, ws)
            end
            cond_prev_endline = n.endline
        elseif n.typ === OPERATOR && n.val == ":"
            if colon_prev_endline != n.endline
                nlen = length(nodes[i-2])
                ws = n.line_offset - (nodes[i-2].line_offset + nlen)
                push!(colon_group, i, n.line_offset, nlen, ws)
            end
            colon_prev_endline = n.endline
        end
    end
    length(cond_group.lens) > 1 || return

    cond_len = align_to(cond_group)
    colon_len = align_to(colon_group)

    cond_len === nothing && colon_len === nothing && return

    if cond_len !== nothing
        for (i, idx) in enumerate(cond_group.node_idxs)
            diff = cond_len - cond_group.lens[i] + 1
            nodes[idx-1] = Whitespace(diff)
        end
    end

    for (i, idx) in enumerate(colon_group.node_idxs)
        # the placeholder would be i+1 if not for a possible inline comment
        nidx = findnext(n -> n.typ === PLACEHOLDER, nodes, idx + 1)
        if nodes[nidx+1].startline != nodes[nidx].startline
            nodes[nidx] = Newline(nest_behavior = AlwaysNest)
        end

        if colon_len !== nothing
            diff = colon_len - colon_group.lens[i] + 1
            nodes[idx-1] = Whitespace(diff)
        end
    end

    fst.nodes = nodes
    fst.nest_behavior = NeverNest
    fst.indent = fst.line_offset - 1

    return
end


text:
macro date(_x) :(println("$(now()) "*$(string(_x)));flush(STDOUT);@time $(esc(_x))) end

# typealias Net Vector{Any}       # vector of model weights

# using Knet

# cpucopy(a::Array)=copy(a)
# cpucopy(a::KnetArray)=Array(a)
# cpucopy(net::Net)=map(cpucopy,net)
# gpucopy(a::Array)=KnetArray(a)
# gpucopy(a::KnetArray)=copy(a)
# gpucopy(net::Net)=map(gpucopy,net)
# testnet(net::Net)=cpucopy(net)

# findprev does not exist in v0.3
# if VERSION < v"0.4-"
# function findprev(A, v, start)
#     for i = start:-1:1
#         A[i] == v && return i
#     end
#     0
# end
# end

# This is "distribute" with an extra procs argument:
# will be fixed in Julia 0.4
# VERSION >= v"0.4-" && eval(Expr(:using,:DistributedArrays))

function distproc(a::AbstractArray, procs)
    owner = myid()
    rr = RemoteRef()
    put!(rr, a)
    d = DArray(size(a), procs) do I
        remotecall_fetch(owner, ()->fetch(rr)[I...])
    end
    # Ensure that all workers have fetched their localparts.
    # Else a gc in between can recover the RemoteRef rr
    for chunk in d.chunks
        wait(chunk)
    end
    d
end

# This is "sortperm!" with AbstractVector argument.
# this has been fixed in Julia 0.4:
# VERSION >= v"0.4.0-" && eval(Expr(:import,:Base,:sortperm!))
using Base: Algorithm, Ordering, DEFAULT_UNSTABLE, Forward, Perm, ord
function sortpermx{I<:Integer}(x::AbstractVector{I}, v::AbstractVector; alg::Algorithm=DEFAULT_UNSTABLE,
                               lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward,
                               initialized::Bool=false)
    length(x) != length(v) && throw(ArgumentError("Index vector must be the same length as the source vector."))
    !initialized && @inbounds for i = 1:length(v); x[i] = i; end
    sort!(x, alg, Perm(ord(lt,by,rev,order),v))
end

# Compute relevant scores
function evalparse(parsers, sentences; ispunct=zn11punct, postag=[])
    ph = map(p->p.head, parsers)
    sh = map(s->s.head, sentences)
    pd = map(p->p.deprel, parsers)
    sd = map(s->s.deprel, sentences)
    wf = map(s->s.word, sentences)
    sp = map(s->s.postag, sentences)
    phcat = vcat2(ph)
    shcat = vcat2(sh)
    pdcat = vcat2(pd)
    sdcat = vcat2(sd)
    wfcat = vcat2(wf)
    spcat = vcat2(sp)
    uas = mean(phcat .== shcat) # unlabeled attachment score
    las = mean((phcat .== shcat) & (pdcat .== sdcat)) # labeled attachment score
    uem = mean(ph .== sh)                             # unlabeled exact match
    lem = mean((ph .== sh) & (pd .== sd))             # labeled exact match
    if ispunct == kcc08punct
        @assert !isempty(postag)
        isword = map(p->!ispunct(p), postag[spcat])
    else
        isword = map(w->!ispunct(w), wfcat)
    end
    uas2 = mean(phcat[isword] .== shcat[isword])
    las2 = mean((phcat[isword] .== shcat[isword]) & (pdcat[isword] .== sdcat[isword]))
    return (uas, uas2, las, las2, uem, lem)
end

function vcat2{T}(v::Vector{Vector{T}})
    w = similar(v[1], 0)
    for vi in v; append!(w, vi); end
    return w
end

# Everybody means something else by "excluding punctuation":
# The buggy conll07 eval.pl script looks for non-alpha characters, but excludes ` and $, also won't match -LRB-, -RRB- etc.
c07punct(w)=(ismatch(r"^\W+$",w) && !ismatch(r"^[`$]+$",w))
# ZN11 eval.py uses his own regexp for punctuation (I think both for chinese and english)
zn11punct(w)=ismatch(r"^[,?!:;]$|^-LRB-$|^-RRB-$|^[.]+$|^[`]+$|^[']+$|^（$|^）$|^、$|^。$|^！$|^？$|^…$|^，$|^；$|^／$|^：$|^“$|^”$|^「$|^」$|^『$|^』$|^《$|^》$|^一一$",w)
# KCC08 uses postags rather than wordforms, but for some reason skips #, $, -LRB-, -RRB- etc.
kcc08punct(p)=in(p,["``", "''", ":", ",", "."])

# Copying net between cpu and gpu

# import Base: copy
# KUnet.GPU && eval(Expr(:using,:CUDArt))

# function copy(net::Net, to::Symbol)
#     net = KUnet.strip!(net)
#     ((to == :gpu) ? gpucopy(net) :
#      (to == :cpu) ? cpucopy(net) :
#      error("Don't know how to copy net to $to"))
# end

# function testnet(net)
#     net = cpucopy(net)
#     KUnet.strip!(net)
#     for l in net
#         for n in fieldnames(l)
#             if isdefined(l,n) && isa(l.(n), KUparam)
#                 # Get rid of all the training fields
#                 l.(n) = KUparam(l.(n).arr)
#             end
#         end
#     end
#     return net
# end


text:
module DualTest

using Test
using Printf
using Random
using ForwardDiff
using ForwardDiff: Partials, Dual, value, partials

using NaNMath, SpecialFunctions, LogExpFunctions
using DiffRules

import Calculus

struct TestTag end

samerng() = MersenneTwister(1)

# By lower-bounding the Int range at 2, we avoid cases where differentiating an
# exponentiation of an Int value would cause a DomainError due to reducing the
# exponent by one
intrand(V) = V == Int ? rand(2:10) : rand(V)

dual_isapprox(a, b) = isapprox(a, b)
dual_isapprox(a::Dual{T,T1,T2}, b::Dual{T,T3,T4}) where {T,T1,T2,T3,T4} = isapprox(value(a), value(b)) && isapprox(partials(a), partials(b))
dual_isapprox(a::Dual{T,T1,T2}, b::Dual{T3,T4,T5}) where {T,T1,T2,T3,T4,T5} = error("Tags don't match")

ForwardDiff.:≺(::Type{TestTag()}, ::Int) = true
ForwardDiff.:≺(::Int, ::Type{TestTag()}) = false

for N in (0,3), M in (0,4), V in (Int, Float32)
    println("  ...testing Dual{TestTag(),$V,$N} and Dual{TestTag(),Dual{TestTag(),$V,$M},$N}")

    PARTIALS = Partials{N,V}(ntuple(n -> intrand(V), N))
    PRIMAL = intrand(V)
    FDNUM = Dual{TestTag()}(PRIMAL, PARTIALS)

    PARTIALS2 = Partials{N,V}(ntuple(n -> intrand(V), N))
    PRIMAL2 = intrand(V)
    FDNUM2 = Dual{TestTag()}(PRIMAL2, PARTIALS2)

    PARTIALS3 = Partials{N,V}(ntuple(n -> intrand(V), N))
    PRIMAL3 = intrand(V)
    FDNUM3 = Dual{TestTag()}(PRIMAL3, PARTIALS3)

    M_PARTIALS = Partials{M,V}(ntuple(m -> intrand(V), M))
    NESTED_PARTIALS = convert(Partials{N,Dual{TestTag(),V,M}}, PARTIALS)
    NESTED_FDNUM = Dual{TestTag()}(Dual{TestTag()}(PRIMAL, M_PARTIALS), NESTED_PARTIALS)

    M_PARTIALS2 = Partials{M,V}(ntuple(m -> intrand(V), M))
    NESTED_PARTIALS2 = convert(Partials{N,Dual{TestTag(),V,M}}, PARTIALS2)
    NESTED_FDNUM2 = Dual{TestTag()}(Dual{TestTag()}(PRIMAL2, M_PARTIALS2), NESTED_PARTIALS2)

    ################
    # Constructors #
    ################

    @test Dual{TestTag()}(PRIMAL, PARTIALS...) === FDNUM
    @test Dual(PRIMAL, PARTIALS...) === Dual{Nothing}(PRIMAL, PARTIALS...)
    @test Dual(PRIMAL) === Dual{Nothing}(PRIMAL)

    @test typeof(Dual{TestTag()}(widen(V)(PRIMAL), PARTIALS)) === Dual{TestTag(),widen(V),N}
    @test typeof(Dual{TestTag()}(widen(V)(PRIMAL), PARTIALS.values)) === Dual{TestTag(),widen(V),N}
    @test typeof(Dual{TestTag()}(widen(V)(PRIMAL), PARTIALS...)) === Dual{TestTag(),widen(V),N}
    @test typeof(NESTED_FDNUM) == Dual{TestTag(),Dual{TestTag(),V,M},N}

    #############
    # Accessors #
    #############

    @test value(PRIMAL) == PRIMAL
    @test value(FDNUM) == PRIMAL
    @test value(NESTED_FDNUM) === Dual{TestTag()}(PRIMAL, M_PARTIALS)

    @test partials(PRIMAL) == Partials{0,V}(tuple())
    @test partials(FDNUM) == PARTIALS
    @test partials(NESTED_FDNUM) === NESTED_PARTIALS

    for i in 1:N
        @test partials(FDNUM, i) == PARTIALS[i]
        for j in 1:M
            @test partials(NESTED_FDNUM, i, j) == partials(NESTED_PARTIALS[i], j)
        end
    end

    @test ForwardDiff.npartials(FDNUM) == N
    @test ForwardDiff.npartials(typeof(FDNUM)) == N
    @test ForwardDiff.npartials(NESTED_FDNUM) == N
    @test ForwardDiff.npartials(typeof(NESTED_FDNUM)) == N

    @test ForwardDiff.valtype(FDNUM) == V
    @test ForwardDiff.valtype(typeof(FDNUM)) == V
    @test ForwardDiff.valtype(NESTED_FDNUM) == Dual{TestTag(),V,M}
    @test ForwardDiff.valtype(typeof(NESTED_FDNUM)) == Dual{TestTag(),V,M}

    #####################
    # Generic Functions #
    #####################

    @test FDNUM === copy(FDNUM)
    @test NESTED_FDNUM === copy(NESTED_FDNUM)

    if V != Int
        @test eps(FDNUM) === eps(PRIMAL)
        @test eps(typeof(FDNUM)) === eps(V)
        @test eps(NESTED_FDNUM) === eps(PRIMAL)
        @test eps(typeof(NESTED_FDNUM)) === eps(V)

        @test floor(Int, FDNUM) === floor(Int, PRIMAL)
        @test floor(Int, FDNUM2) === floor(Int, PRIMAL2)
        @test floor(Int, NESTED_FDNUM) === floor(Int, PRIMAL)

        @test floor(FDNUM) === floor(PRIMAL)
        @test floor(FDNUM2) === floor(PRIMAL2)
        @test floor(NESTED_FDNUM) === floor(PRIMAL)

        @test ceil(Int, FDNUM) === ceil(Int, PRIMAL)
        @test ceil(Int, FDNUM2) === ceil(Int, PRIMAL2)
        @test ceil(Int, NESTED_FDNUM) === ceil(Int, PRIMAL)

        @test ceil(FDNUM) === ceil(PRIMAL)
        @test ceil(FDNUM2) === ceil(PRIMAL2)
        @test ceil(NESTED_FDNUM) === ceil(PRIMAL)

        @test trunc(Int, FDNUM) === trunc(Int, PRIMAL)
        @test trunc(Int, FDNUM2) === trunc(Int, PRIMAL2)
        @test trunc(Int, NESTED_FDNUM) === trunc(Int, PRIMAL)

        @test trunc(FDNUM) === trunc(PRIMAL)
        @test trunc(FDNUM2) === trunc(PRIMAL2)
        @test trunc(NESTED_FDNUM) === trunc(PRIMAL)

        @test round(Int, FDNUM) === round(Int, PRIMAL)
        @test round(Int, FDNUM2) === round(Int, PRIMAL2)
        @test round(Int, NESTED_FDNUM) === round(Int, PRIMAL)

        @test round(FDNUM) === round(PRIMAL)
        @test round(FDNUM2) === round(PRIMAL2)
        @test round(NESTED_FDNUM) === round(PRIMAL)

        @test div(FDNUM, FDNUM2) === div(PRIMAL, PRIMAL2)
        @test div(FDNUM, PRIMAL2) === div(PRIMAL, PRIMAL2)
        @test div(PRIMAL, FDNUM2) === div(PRIMAL, PRIMAL2)

        @test div(NESTED_FDNUM, NESTED_FDNUM2) === div(PRIMAL, PRIMAL2)
        @test div(NESTED_FDNUM, PRIMAL2) === div(PRIMAL, PRIMAL2)
        @test div(PRIMAL, NESTED_FDNUM2) === div(PRIMAL, PRIMAL2)

        if VERSION ≥ v"1.4"
            @test div(FDNUM, FDNUM2, RoundUp) === div(PRIMAL, PRIMAL2, RoundUp)
            @test div(NESTED_FDNUM, NESTED_FDNUM2, RoundUp) === div(PRIMAL, PRIMAL2, RoundUp)
        end

        @test Base.rtoldefault(typeof(FDNUM)) ≡ Base.rtoldefault(typeof(PRIMAL))
        @test Dual{TestTag()}(PRIMAL-eps(V), PARTIALS) ≈ FDNUM
        @test Base.rtoldefault(typeof(NESTED_FDNUM)) ≡ Base.rtoldefault(typeof(PRIMAL))
    end

    @test hash(FDNUM) === hash(PRIMAL)
    @test hash(FDNUM, hash(PRIMAL)) === hash(PRIMAL, hash(PRIMAL))
    @test hash(NESTED_FDNUM) === hash(PRIMAL)
    @test hash(NESTED_FDNUM, hash(PRIMAL)) === hash(PRIMAL, hash(PRIMAL))

    TMPIO = IOBuffer()
    write(TMPIO, FDNUM)
    seekstart(TMPIO)
    @test read(TMPIO, typeof(FDNUM)) === FDNUM
    seekstart(TMPIO)
    write(TMPIO, FDNUM2)
    seekstart(TMPIO)
    @test read(TMPIO, typeof(FDNUM2)) === FDNUM2
    seekstart(TMPIO)
    write(TMPIO, NESTED_FDNUM)
    seekstart(TMPIO)
    @test read(TMPIO, typeof(NESTED_FDNUM)) === NESTED_FDNUM
    close(TMPIO)

    @test zero(FDNUM) === Dual{TestTag()}(zero(PRIMAL), zero(PARTIALS))
    @test zero(typeof(FDNUM)) === Dual{TestTag()}(zero(V), zero(Partials{N,V}))
    @test zero(NESTED_FDNUM) === Dual{TestTag()}(Dual{TestTag()}(zero(PRIMAL), zero(M_PARTIALS)), zero(NESTED_PARTIALS))
    @test zero(typeof(NESTED_FDNUM)) === Dual{TestTag()}(Dual{TestTag()}(zero(V), zero(Partials{M,V})), zero(Partials{N,Dual{TestTag(),V,M}}))

    @test one(FDNUM) === Dual{TestTag()}(one(PRIMAL), zero(PARTIALS))
    @test one(typeof(FDNUM)) === Dual{TestTag()}(one(V), zero(Partials{N,V}))
    @test one(NESTED_FDNUM) === Dual{TestTag()}(Dual{TestTag()}(one(PRIMAL), zero(M_PARTIALS)), zero(NESTED_PARTIALS))
    @test one(typeof(NESTED_FDNUM)) === Dual{TestTag()}(Dual{TestTag()}(one(V), zero(Partials{M,V})), zero(Partials{N,Dual{TestTag(),V,M}}))

    if V <: Integer
        @test rand(samerng(), FDNUM) == rand(samerng(), value(FDNUM))
        @test rand(samerng(), NESTED_FDNUM) == rand(samerng(), value(NESTED_FDNUM))
    elseif V <: AbstractFloat
        @test rand(samerng(), typeof(FDNUM)) === Dual{TestTag()}(rand(samerng(), V), zero(Partials{N,V}))
        @test rand(samerng(), typeof(NESTED_FDNUM)) === Dual{TestTag()}(Dual{TestTag()}(rand(samerng(), V), zero(Partials{M,V})), zero(Partials{N,Dual{TestTag(),V,M}}))
        @test randn(samerng(), typeof(FDNUM)) === Dual{TestTag()}(randn(samerng(), V), zero(Partials{N,V}))
        @test randn(samerng(), typeof(NESTED_FDNUM)) === Dual{TestTag()}(Dual{TestTag()}(randn(samerng(), V), zero(Partials{M,V})),
        zero(Partials{N,Dual{TestTag(),V,M}}))
        @test randexp(samerng(), typeof(FDNUM)) === Dual{TestTag()}(randexp(samerng(), V), zero(Partials{N,V}))
        @test randexp(samerng(), typeof(NESTED_FDNUM)) === Dual{TestTag()}(Dual{TestTag()}(randexp(samerng(), V), zero(Partials{M,V})),
        zero(Partials{N,Dual{TestTag(),V,M}}))
    end

    # Predicates #
    #------------#

    @test ForwardDiff.isconstant(zero(FDNUM))
    @test ForwardDiff.isconstant(one(FDNUM))
    @test ForwardDiff.isconstant(FDNUM) == (N == 0)

    @test ForwardDiff.isconstant(zero(NESTED_FDNUM))
    @test ForwardDiff.isconstant(one(NESTED_FDNUM))
    @test ForwardDiff.isconstant(NESTED_FDNUM) == (N == 0)

    @test isequal(FDNUM, Dual{TestTag()}(PRIMAL, PARTIALS2))
    @test isequal(PRIMAL, PRIMAL2) == isequal(FDNUM, FDNUM2)

    @test isequal(NESTED_FDNUM, Dual{TestTag()}(Dual{TestTag()}(PRIMAL, M_PARTIALS2), NESTED_PARTIALS2))
    @test isequal(PRIMAL, PRIMAL2) == isequal(NESTED_FDNUM, NESTED_FDNUM2)

    @test FDNUM == Dual{TestTag()}(PRIMAL, PARTIALS2)
    @test (PRIMAL == PRIMAL2) == (FDNUM == FDNUM2)
    @test (PRIMAL == PRIMAL2) == (NESTED_FDNUM == NESTED_FDNUM2)

    @test isless(Dual{TestTag()}(1, PARTIALS), Dual{TestTag()}(2, PARTIALS2))
    @test !(isless(Dual{TestTag()}(1, PARTIALS), Dual{TestTag()}(1, PARTIALS2)))
    @test !(isless(Dual{TestTag()}(2, PARTIALS), Dual{TestTag()}(1, PARTIALS2)))

    @test isless(Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS), Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS2), NESTED_PARTIALS2))
    @test !(isless(Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS), Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2)))
    @test !(isless(Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS), NESTED_PARTIALS), Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2)))

    @test Dual{TestTag()}(1, PARTIALS) < Dual{TestTag()}(2, PARTIALS2)
    @test !(Dual{TestTag()}(1, PARTIALS) < Dual{TestTag()}(1, PARTIALS2))
    @test !(Dual{TestTag()}(2, PARTIALS) < Dual{TestTag()}(1, PARTIALS2))

    @test Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) < Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS2), NESTED_PARTIALS2)
    @test !(Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) < Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2))
    @test !(Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS), NESTED_PARTIALS) < Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2))

    @test Dual{TestTag()}(1, PARTIALS) <= Dual{TestTag()}(2, PARTIALS2)
    @test Dual{TestTag()}(1, PARTIALS) <= Dual{TestTag()}(1, PARTIALS2)
    @test !(Dual{TestTag()}(2, PARTIALS) <= Dual{TestTag()}(1, PARTIALS2))

    @test Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) <= Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS2), NESTED_PARTIALS2)
    @test Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) <= Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2)
    @test !(Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS), NESTED_PARTIALS) <= Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2))

    @test Dual{TestTag()}(2, PARTIALS) > Dual{TestTag()}(1, PARTIALS2)
    @test !(Dual{TestTag()}(1, PARTIALS) > Dual{TestTag()}(1, PARTIALS2))
    @test !(Dual{TestTag()}(1, PARTIALS) > Dual{TestTag()}(2, PARTIALS2))

    @test Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS), NESTED_PARTIALS) > Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2)
    @test !(Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) > Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2))
    @test !(Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) > Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS2), NESTED_PARTIALS2))

    @test Dual{TestTag()}(2, PARTIALS) >= Dual{TestTag()}(1, PARTIALS2)
    @test Dual{TestTag()}(1, PARTIALS) >= Dual{TestTag()}(1, PARTIALS2)
    @test !(Dual{TestTag()}(1, PARTIALS) >= Dual{TestTag()}(2, PARTIALS2))

    @test Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS), NESTED_PARTIALS) >= Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2)
    @test Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) >= Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS2), NESTED_PARTIALS2)
    @test !(Dual{TestTag()}(Dual{TestTag()}(1, M_PARTIALS), NESTED_PARTIALS) >= Dual{TestTag()}(Dual{TestTag()}(2, M_PARTIALS2), NESTED_PARTIALS2))

    @test isnan(Dual{TestTag()}(NaN, PARTIALS))
    @test !(isnan(FDNUM))

    @test isnan(Dual{TestTag()}(Dual{TestTag()}(NaN, M_PARTIALS), NESTED_PARTIALS))
    @test !(isnan(NESTED_FDNUM))

    @test isfinite(FDNUM)
    @test !(isfinite(Dual{TestTag()}(Inf, PARTIALS)))

    @test isfinite(NESTED_FDNUM)
    @test !(isfinite(Dual{TestTag()}(Dual{TestTag()}(NaN, M_PARTIALS), NESTED_PARTIALS)))

    @test isinf(Dual{TestTag()}(Inf, PARTIALS))
    @test !(isinf(FDNUM))

    @test isinf(Dual{TestTag()}(Dual{TestTag()}(Inf, M_PARTIALS), NESTED_PARTIALS))
    @test !(isinf(NESTED_FDNUM))

    @test isreal(FDNUM)
    @test isreal(NESTED_FDNUM)

    @test isinteger(Dual{TestTag()}(1.0, PARTIALS))
    @test isinteger(FDNUM) == (V == Int)

    @test isinteger(Dual{TestTag()}(Dual{TestTag()}(1.0, M_PARTIALS), NESTED_PARTIALS))
    @test isinteger(NESTED_FDNUM) == (V == Int)

    @test iseven(Dual{TestTag()}(2))
    @test !(iseven(Dual{TestTag()}(1)))

    @test iseven(Dual{TestTag()}(Dual{TestTag()}(2)))
    @test !(iseven(Dual{TestTag()}(Dual{TestTag()}(1))))

    @test isodd(Dual{TestTag()}(1))
    @test !(isodd(Dual{TestTag()}(2)))

    @test isodd(Dual{TestTag()}(Dual{TestTag()}(1)))
    @test !(isodd(Dual{TestTag()}(Dual{TestTag()}(2))))

    ########################
    # Promotion/Conversion #
    ########################

    WIDE_T = widen(V)

    @test promote_type(Dual{TestTag(),V,N}, V) == Dual{TestTag(),V,N}
    @test promote_type(Dual{TestTag(),V,N}, WIDE_T) == Dual{TestTag(),WIDE_T,N}
    @test promote_type(Dual{TestTag(),WIDE_T,N}, V) == Dual{TestTag(),WIDE_T,N}
    @test promote_type(Dual{TestTag(),V,N}, Dual{TestTag(),V,N}) == Dual{TestTag(),V,N}
    @test promote_type(Dual{TestTag(),V,N}, Dual{TestTag(),WIDE_T,N}) == Dual{TestTag(),WIDE_T,N}
    @test promote_type(Dual{TestTag(),WIDE_T,N}, Dual{TestTag(),Dual{TestTag(),V,M},N}) == Dual{TestTag(),Dual{TestTag(),WIDE_T,M},N}

    # issue #322
    @test promote_type(Bool, Dual{TestTag(),V,N}) == Dual{TestTag(),promote_type(Bool, V),N}
    @test promote_type(BigFloat, Dual{TestTag(),V,N}) == Dual{TestTag(),promote_type(BigFloat, V),N}

    WIDE_FDNUM = convert(Dual{TestTag(),WIDE_T,N}, FDNUM)
    WIDE_NESTED_FDNUM = convert(Dual{TestTag(),Dual{TestTag(),WIDE_T,M},N}, NESTED_FDNUM)

    @test typeof(WIDE_FDNUM) === Dual{TestTag(),WIDE_T,N}
    @test typeof(WIDE_NESTED_FDNUM) === Dual{TestTag(),Dual{TestTag(),WIDE_T,M},N}

    @test value(WIDE_FDNUM) == PRIMAL
    @test value(WIDE_NESTED_FDNUM) == PRIMAL

    @test convert(Dual, FDNUM) === FDNUM
    @test convert(Dual, NESTED_FDNUM) === NESTED_FDNUM
    @test convert(Dual{TestTag(),V,N}, FDNUM) === FDNUM
    @test convert(Dual{TestTag(),Dual{TestTag(),V,M},N}, NESTED_FDNUM) === NESTED_FDNUM
    @test convert(Dual{TestTag(),WIDE_T,N}, PRIMAL) === Dual{TestTag()}(WIDE_T(PRIMAL), zero(Partials{N,WIDE_T}))
    @test convert(Dual{TestTag(),Dual{TestTag(),WIDE_T,M},N}, PRIMAL) === Dual{TestTag()}(Dual{TestTag()}(WIDE_T(PRIMAL), zero(Partials{M,WIDE_T})), zero(Partials{N,Dual{TestTag(),V,M}}))
    @test convert(Dual{TestTag(),Dual{TestTag(),V,M},N}, FDNUM) === Dual{TestTag()}(convert(Dual{TestTag(),V,M}, PRIMAL), convert(Partials{N,Dual{TestTag(),V,M}}, PARTIALS))
    @test convert(Dual{TestTag(),Dual{TestTag(),WIDE_T,M},N}, FDNUM) === Dual{TestTag()}(convert(Dual{TestTag(),WIDE_T,M}, PRIMAL), convert(Partials{N,Dual{TestTag(),WIDE_T,M}}, PARTIALS))

    ##############
    # Arithmetic #
    ##############

    # Addition/Subtraction #
    #----------------------#

    @test FDNUM + FDNUM2 === Dual{TestTag()}(value(FDNUM) + value(FDNUM2), partials(FDNUM) + partials(FDNUM2))
    @test FDNUM + PRIMAL === Dual{TestTag()}(value(FDNUM) + PRIMAL, partials(FDNUM))
    @test PRIMAL + FDNUM === Dual{TestTag()}(value(FDNUM) + PRIMAL, partials(FDNUM))

    @test NESTED_FDNUM + NESTED_FDNUM2 === Dual{TestTag()}(value(NESTED_FDNUM) + value(NESTED_FDNUM2), partials(NESTED_FDNUM) + partials(NESTED_FDNUM2))
    @test NESTED_FDNUM + PRIMAL === Dual{TestTag()}(value(NESTED_FDNUM) + PRIMAL, partials(NESTED_FDNUM))
    @test PRIMAL + NESTED_FDNUM === Dual{TestTag()}(value(NESTED_FDNUM) + PRIMAL, partials(NESTED_FDNUM))

    @test FDNUM - FDNUM2 === Dual{TestTag()}(value(FDNUM) - value(FDNUM2), partials(FDNUM) - partials(FDNUM2))
    @test FDNUM - PRIMAL === Dual{TestTag()}(value(FDNUM) - PRIMAL, partials(FDNUM))
    @test PRIMAL - FDNUM === Dual{TestTag()}(PRIMAL - value(FDNUM), -(partials(FDNUM)))
    @test -(FDNUM) === Dual{TestTag()}(-(value(FDNUM)), -(partials(FDNUM)))

    @test NESTED_FDNUM - NESTED_FDNUM2 === Dual{TestTag()}(value(NESTED_FDNUM) - value(NESTED_FDNUM2), partials(NESTED_FDNUM) - partials(NESTED_FDNUM2))
    @test NESTED_FDNUM - PRIMAL === Dual{TestTag()}(value(NESTED_FDNUM) - PRIMAL, partials(NESTED_FDNUM))
    @test PRIMAL - NESTED_FDNUM === Dual{TestTag()}(PRIMAL - value(NESTED_FDNUM), -(partials(NESTED_FDNUM)))
    @test -(NESTED_FDNUM) === Dual{TestTag()}(-(value(NESTED_FDNUM)), -(partials(NESTED_FDNUM)))

    # Multiplication #
    #----------------#

    @test FDNUM * FDNUM2 === Dual{TestTag()}(value(FDNUM) * value(FDNUM2), ForwardDiff._mul_partials(partials(FDNUM), partials(FDNUM2), value(FDNUM2), value(FDNUM)))
    @test FDNUM * PRIMAL === Dual{TestTag()}(value(FDNUM) * PRIMAL, partials(FDNUM) * PRIMAL)
    @test PRIMAL * FDNUM === Dual{TestTag()}(value(FDNUM) * PRIMAL, partials(FDNUM) * PRIMAL)

    @test NESTED_FDNUM * NESTED_FDNUM2 === Dual{TestTag()}(value(NESTED_FDNUM) * value(NESTED_FDNUM2), ForwardDiff._mul_partials(partials(NESTED_FDNUM), partials(NESTED_FDNUM2), value(NESTED_FDNUM2), value(NESTED_FDNUM)))
    @test NESTED_FDNUM * PRIMAL === Dual{TestTag()}(value(NESTED_FDNUM) * PRIMAL, partials(NESTED_FDNUM) * PRIMAL)
    @test PRIMAL * NESTED_FDNUM === Dual{TestTag()}(value(NESTED_FDNUM) * PRIMAL, partials(NESTED_FDNUM) * PRIMAL)

    # Division #
    #----------#

    if M > 0 && N > 0
        @test Dual{1}(FDNUM) / Dual{1}(PRIMAL) === Dual{1}(FDNUM / PRIMAL)
        @test Dual{1}(PRIMAL) / Dual{1}(FDNUM) === Dual{1}(PRIMAL / FDNUM)
        @test_broken Dual{1}(FDNUM) / FDNUM2 === Dual{1}(FDNUM / FDNUM2)
        @test_broken FDNUM / Dual{1}(FDNUM2) === Dual{1}(FDNUM / FDNUM2)
        # following may not be exact, see #264
        @test Dual{1}(FDNUM / PRIMAL, FDNUM2 / PRIMAL) ≈ Dual{1}(FDNUM, FDNUM2) / PRIMAL
    end

    @test dual_isapprox(FDNUM / FDNUM2, Dual{TestTag()}(value(FDNUM) / value(FDNUM2), ForwardDiff._div_partials(partials(FDNUM), partials(FDNUM2), value(FDNUM), value(FDNUM2))))
    @test dual_isapprox(FDNUM / PRIMAL, Dual{TestTag()}(value(FDNUM) / PRIMAL, partials(FDNUM) / PRIMAL))
    @test dual_isapprox(PRIMAL / FDNUM, Dual{TestTag()}(PRIMAL / value(FDNUM), (-(PRIMAL) / value(FDNUM)^2) * partials(FDNUM)))

    @test dual_isapprox(NESTED_FDNUM / NESTED_FDNUM2, Dual{TestTag()}(value(NESTED_FDNUM) / value(NESTED_FDNUM2), ForwardDiff._div_partials(partials(NESTED_FDNUM), partials(NESTED_FDNUM2), value(NESTED_FDNUM), value(NESTED_FDNUM2))))
    @test dual_isapprox(NESTED_FDNUM / PRIMAL, Dual{TestTag()}(value(NESTED_FDNUM) / PRIMAL, partials(NESTED_FDNUM) / PRIMAL))
    @test dual_isapprox(PRIMAL / NESTED_FDNUM, Dual{TestTag()}(PRIMAL / value(NESTED_FDNUM), (-(PRIMAL) / value(NESTED_FDNUM)^2) * partials(NESTED_FDNUM)))

    # Exponentiation #
    #----------------#
    # If V == Int, the LHS terms are Int's. Large inputs cause integer overflow
    # within the generic fallback of `isapprox`, resulting in a DomainError.
    # Promote to Float64 to avoid issues.
    @test dual_isapprox(1.0 * FDNUM^FDNUM2, exp(FDNUM2 * log(FDNUM)))
    @test dual_isapprox(1.0 * FDNUM^PRIMAL, exp(PRIMAL * log(FDNUM)))
    @test dual_isapprox(1.0 * PRIMAL^FDNUM, exp(FDNUM * log(PRIMAL)))

    @test dual_isapprox(1.0 * NESTED_FDNUM^NESTED_FDNUM2, exp(NESTED_FDNUM2 * log(NESTED_FDNUM)))
    @test dual_isapprox(1.0 * NESTED_FDNUM^PRIMAL, exp(PRIMAL * log(NESTED_FDNUM)))
    @test dual_isapprox(1.0 * PRIMAL^NESTED_FDNUM, exp(NESTED_FDNUM * log(PRIMAL)))

    @test partials(NaNMath.pow(Dual{TestTag()}(-2.0, 1.0), Dual{TestTag()}(2.0, 0.0)), 1) == -4.0

    ###################################
    # General Mathematical Operations #
    ###################################

    @test conj(FDNUM) === FDNUM
    @test conj(NESTED_FDNUM) === NESTED_FDNUM

    @test transpose(FDNUM) === FDNUM
    @test transpose(NESTED_FDNUM) === NESTED_FDNUM

    @test abs(-FDNUM) === FDNUM
    @test abs(FDNUM) === FDNUM
    @test abs(-NESTED_FDNUM) === NESTED_FDNUM
    @test abs(NESTED_FDNUM) === NESTED_FDNUM

    if V != Int
        for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing)
            if f in (:hankelh1, :hankelh1x, :hankelh2, :hankelh2x, :/, :rem2pi)
                continue  # Skip these rules
            elseif !(isdefined(@__MODULE__, M) && isdefined(getfield(@__MODULE__, M), f))
                continue  # Skip rules for methods not defined in the current scope
            end
            println("       ...auto-testing $(M).$(f) with $arity arguments")
            if arity == 1
                deriv = DiffRules.diffrule(M, f, :x)
                modifier = if in(f, (:asec, :acsc, :asecd, :acscd, :acosh, :acoth))
                    one(V)
                elseif in(f, (:log1mexp, :log2mexp))
                    -one(V)
                else
                    zero(V)
                end
                @eval begin
                    x = rand() + $modifier
                    dx = $M.$f(Dual{TestTag()}(x, one(x)))
                    @test value(dx) == $M.$f(x)
                    @test partials(dx, 1) == $deriv
                end
            elseif arity == 2
                derivs = DiffRules.diffrule(M, f, :x, :y)
                @eval begin
                    x, y = rand(1:10), rand()
                    dx = $M.$f(Dual{TestTag()}(x, one(x)), y)
                    dy = $M.$f(x, Dual{TestTag()}(y, one(y)))
                    actualdx = $(derivs[1])
                    actualdy = $(derivs[2])
                    @test value(dx) == $M.$f(x, y)
                    @test value(dy) == value(dx)
                    if isnan(actualdx)
                        @test isnan(partials(dx, 1))
                    else
                        @test partials(dx, 1) ≈ actualdx
                    end
                    if isnan(actualdy)
                        @test isnan(partials(dy, 1))
                    else
                        @test partials(dy, 1) ≈ actualdy
                    end
                end
            end
        end
    end

    # Special Cases #
    #---------------#

    @test_broken dual_isapprox(hypot(FDNUM, FDNUM2, FDNUM), sqrt(2*(FDNUM^2) + FDNUM2^2))
    @test_broken dual_isapprox(hypot(FDNUM, FDNUM2, FDNUM3), sqrt(FDNUM^2 + FDNUM2^2 + FDNUM3^2))

    @test all(map(dual_isapprox, ForwardDiff.sincos(FDNUM), (sin(FDNUM), cos(FDNUM))))

    if VERSION >= v"1.6.0-DEV.292"
        @test all(map(dual_isapprox, sincospi(FDNUM), (sinpi(FDNUM), cospi(FDNUM))))
    end

    if V === Float32
        @test typeof(sqrt(FDNUM)) === typeof(FDNUM)
        @test typeof(sqrt(NESTED_FDNUM)) === typeof(NESTED_FDNUM)
    end

    for f in (fma, muladd)
        @test dual_isapprox(f(FDNUM, FDNUM2, FDNUM3),   Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2 + PRIMAL2*PARTIALS + PARTIALS3))
        @test dual_isapprox(f(FDNUM, FDNUM2, PRIMAL3),  Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2 + PRIMAL2*PARTIALS))
        @test dual_isapprox(f(PRIMAL, FDNUM2, FDNUM3),  Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2 + PARTIALS3))
        @test dual_isapprox(f(PRIMAL, FDNUM2, PRIMAL3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2))
        @test dual_isapprox(f(FDNUM, PRIMAL2, FDNUM3),  Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL2*PARTIALS + PARTIALS3))
        @test dual_isapprox(f(FDNUM, PRIMAL2, PRIMAL3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL2*PARTIALS))
        @test dual_isapprox(f(PRIMAL, PRIMAL2, FDNUM3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PARTIALS3))
    end
end

@testset "Exponentiation of zero" begin
    x0 = 0.0
    x1 = Dual{:t1}(x0, 1.0)
    x2 = Dual{:t2}(x1, 1.0)
    x3 = Dual{:t3}(x2, 1.0)
    pow = ^  # to call non-literal power
    @test pow(x3, 2) === x3^2 === x3 * x3
    @test pow(x2, 1) === x2^1 === x2
    @test pow(x1, 0) === x1^0 === Dual{:t1}(1.0, 0.0)
end

@testset "Type min/max" begin
    d1 = Dual(1.0)
    dinf = typemax(typeof(d1))
    dminf = typemin(typeof(d1))
    @test dminf < d1 < dinf
    @test typeof(dminf) === typeof(d1)
    @test typeof(dinf) === typeof(d1)
    @test !isfinite(dminf)
    @test !isfinite(dinf)

    dfmin = floatmin(typeof(d1))
    dfmax = floatmax(typeof(d1))
    @test dfmin < d1 < dfmax
    @test typeof(dfmin) === typeof(d1)
    @test typeof(dfmax) === typeof(d1)
    @test isfinite(dfmin)
    @test isfinite(dfmax)
end

@testset "Integer" begin
    x = Dual(1.0,0)
    @test Int(x) ≡ Integer(x) ≡ convert(Int,x) ≡ convert(Integer,x) ≡ 1
    x = Dual(1.0,0,0)
    @test Int(x) ≡ Integer(x) ≡ convert(Int,x) ≡ convert(Integer,x) ≡ 1
    @test_throws InexactError Int(Dual(1.5,1.2))
    @test_throws InexactError Integer(Dual(1.5,1.2))
    @test_throws InexactError Int(Dual(1,1))
    @test_throws InexactError Integer(Dual(1,1,2))
    @test length(UnitRange(Dual(1.5), Dual(3.5))) == 3
    @test length(UnitRange(Dual(1.5,1), Dual(3.5,3))) == 3
end

if VERSION >= v"1.6.0-rc1"
    @testset "@printf" begin
        for T in (Float16, Float32, Float64, BigFloat)
            d1 = Dual(one(T))
            @test_nowarn @printf("Testing @printf: %.2e\n", d1)
            @test @sprintf("Testing @sprintf: %.2e\n", d1) == "Testing @sprintf: 1.00e+00\n"
        end
    end
end

@testset "float" begin # issue #492
    @test float(Dual{Nothing, Int, 2}) === Dual{Nothing, Float64, 2}
    @test float(Dual(1)) isa Dual{Nothing, Float64, 0}
    @test value.(float.(Dual.(1:4, 2:5, 3:6))) isa Vector{Float64}
    @test ForwardDiff.derivative(float, 1)::Float64 === 1.0
end

end # module