text:
#####
##### Model
#####
crc3, parameters, environment, alphas, stepsymbol, cl_stall = MT.build_crc3(; nwings)
crc3_model = MT.AircraftDynamics(crc3, parameters, environment, stepsymbol, cl_stall)

function X_2_Xdot(model, alpha, gamma; vinf = 1.0, rotors_on = false, U = zeros(length(crc3_model.aircraft.rotor_system.index)) * 2*pi/60.0)
    X = MT.steady_climb_state(alpha, gamma, vinf)
    Xdot = MT.dynamics(model, X, U; rotors_on)
    return X, Xdot
end

function test_xdot_cl_alpha(crc3_model;
    alphas = range(-40.0, stop=55.0, length=26) .* pi/180,
    gamma = 0 * pi/180,
    vinf = 1.0,
    rotors_on = false,
    U = zeros(length(crc3_model.aircraft.rotor_system.index)) * 2*pi/60.0,
    save_figures = true, fig_name = "test_dynamics.png", clear_figure = true, label_tag = ""
)
    # initialize vars
    Xs = zeros(7,length(alphas))
    Xdots = zeros(7,length(alphas))

    # build lift curves
    for (i,alpha) in enumerate(alphas)
        X, Xdot = X_2_Xdot(crc3_model, alpha, gamma; vinf, rotors_on, U)
        Xs[:,i] .= X
        Xdots[:,i] .= Xdot
    end

    # get forces to check (in global frame)
    _, CFs, CMs = test_cl_alpha(crc3_model; alphas, vinf, rotors_on, U, frame = MT.Global(), save_figures = false)

    # dimensionalize forces
    ρ = crc3_model.environment.ρ
    g = crc3_model.environment.g
    qinf = 0.5 * ρ * vinf^2
    m = crc3_model.aircraft.inertia_system.mass
    I = crc3_model.aircraft.inertia_system.inertia_y

    weight = MT.ST.@SVector [0.0, -g * m] # in global frame
    Fs = CFs .* (crc3_model.aircraft.wing_system.system.reference[1].S * qinf) .+ weight
    Ms = CMs .* (crc3_model.aircraft.wing_system.system.reference[1].S * crc3_model.aircraft.wing_system.system.reference[1].c * qinf)
    # predict dynamics: F = ma, T = Iα
    p_dots = deepcopy(Xs[3:4,:]) # in the body frame
    for (i,alpha) in enumerate(alphas)
        local theta = Xs[5,i]
        local R = MT.R_from_body(MT.Global(), alpha, theta)
        p_dots[:,i] .= R * p_dots[:,i] # rotate to global frame
    end
    v_dots = Fs ./ m # in the global frame
    for (i,alpha) in enumerate(alphas)
        local theta = Xs[5,i]
        local theta_dot = Xs[6,i]
        local R = MT.R_from_global(MT.Body(), alpha, theta)
        v_dots[:,i] .= R * v_dots[:,i] # rotate to body frame
            - theta_dot .* (MT.ST.@SArray [0. -1.; 1. 0.]) * Xs[3:4,i] # include coriolis force
    end
    θ_dots = Xs[6,:]
    θ_dot_dots = Ms ./ I
    Xdots_check = vcat(p_dots, v_dots, θ_dots', θ_dot_dots')

    fig_X_Xdot = PP.figure("test_Xdot")
    if clear_figure; fig_X_Xdot.clear(); end

    fig_X_Xdot.add_subplot(2, 3, 1, ylabel = MT.LS.L"P_x[m]")
    fig_X_Xdot.add_subplot(2, 3, 2, ylabel = MT.LS.L"V_x[m/s]")
    fig_X_Xdot.add_subplot(2, 3, 3, ylabel = MT.LS.L"\theta[^\circ]")
    fig_X_Xdot.add_subplot(2, 3, 4, ylabel = MT.LS.L"P_y[m]", xlabel = MT.LS.L"\alpha[^\circ]")
    fig_X_Xdot.add_subplot(2, 3, 5, ylabel = MT.LS.L"V_y[m/s]", xlabel = MT.LS.L"\alpha[^\circ]")
    fig_X_Xdot.add_subplot(2, 3, 6, ylabel = MT.LS.L"\dot{\theta}[^\circ/s]", xlabel = MT.LS.L"\alpha[^\circ]")

    axs_X_Xdot = fig_X_Xdot.get_axes()
    axs_X_Xdot[1].plot(alphas .* 180/pi, Xs[1,:], label = MT.LS.L"X" * label_tag)
    axs_X_Xdot[1].plot(alphas .* 180/pi, Xdots[1,:], label = MT.LS.L"\dot{X}" * label_tag)
    axs_X_Xdot[1].plot(alphas .* 180/pi, Xdots_check[1,:], label = "verify " * MT.LS.L"\dot{X}" * label_tag, "--")
    axs_X_Xdot[2].plot(alphas .* 180/pi, Xs[3,:], label = nothing)
    axs_X_Xdot[2].plot(alphas .* 180/pi, Xdots[3,:], label = nothing)
    axs_X_Xdot[2].plot(alphas .* 180/pi, Xdots_check[3,:], label = nothing, "--")
    axs_X_Xdot[3].plot(alphas .* 180/pi, Xs[5,:] .* 180/pi, label = nothing)
    axs_X_Xdot[3].plot(alphas .* 180/pi, Xdots[5,:] .* 180/pi, label = nothing)
    axs_X_Xdot[3].plot(alphas .* 180/pi, Xdots_check[5,:] .* 180/pi, label = nothing, "--")
    axs_X_Xdot[4].plot(alphas .* 180/pi, Xs[2,:], label = nothing)
    axs_X_Xdot[4].plot(alphas .* 180/pi, Xdots[2,:], label = nothing)
    axs_X_Xdot[4].plot(alphas .* 180/pi, Xdots_check[2,:], label = nothing, "--")
    axs_X_Xdot[5].plot(alphas .* 180/pi, Xs[4,:], label = nothing)
    axs_X_Xdot[5].plot(alphas .* 180/pi, Xdots[4,:], label = nothing)
    axs_X_Xdot[5].plot(alphas .* 180/pi, Xdots_check[4,:], label = nothing, "--")
    axs_X_Xdot[6].plot(alphas .* 180/pi, Xs[6,:] .* 180/pi, label = nothing)
    axs_X_Xdot[6].plot(alphas .* 180/pi, Xdots[6,:] .* 180/pi, label = nothing)
    axs_X_Xdot[6].plot(alphas .* 180/pi, Xdots_check[6,:] .* 180/pi, label = nothing, "--")

    fig_X_Xdot.legend(loc="upper left", bbox_to_anchor=(0.73,0.97))
    fig_X_Xdot.suptitle(MT.LS.L"\gamma = " * "$(gamma * 180/pi)" * MT.LS.L"^\circ")
    fig_X_Xdot.tight_layout()
    fig_X_Xdot.set_size_inches(12, 7, forward=true)

    if save_figures
        fig_X_Xdot.savefig(joinpath(MT.topdirectory, "data", "plots", MT.TODAY, fig_name))
    end

    return Xs, Xdots, Xdots_check
end

Xs, Xdots, Xdots_check = test_xdot_cl_alpha(crc3_model; gamma, fig_name = "test_dynamics_$(nwings)_wing_gamma_$(Int(round(gamma; digits = 0))).png")
_ = test_xdot_cl_alpha(crc3_model;
    alphas = range(-40.0, stop=55.0, length=26) .* pi/180,
    gamma = 0 * pi/180,
    vinf = 1.0,
    rotors_on = true,
    U = ones(length(crc3_model.aircraft.rotor_system.index)) * 3000 * 2*pi/60.0,
    save_figures = true, fig_name = "test_dynamics_$(nwings)_wing_gamma_$(Int(round(gamma; digits = 0)))_rotors_on.png", clear_figure = false, label_tag = " w/ rotors"
)

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
"""
    select(t::Table, which::Selection)

Select all or a subset of columns, or a single column from the table.

`Selection` is a type union of many types that can select from a table. It can be:

1. `Integer` -- returns the column at this position.
2. `Symbol` -- returns the column with this name.
3. `Pair{Selection => Function}` -- selects and maps a function over the selection, returns the result.
4. `AbstractArray` -- returns the array itself. This must be the same length as the table.
5. `Tuple` of `Selection` -- returns a table containing a column for every selector in the tuple.
6. `Regex` -- returns the columns with names that match the regular expression.
7. `Type` -- returns columns with elements of the given type.

# Examples:

    t = table(1:10, randn(10), rand(Bool, 10); names = [:x, :y, :z])

    # select the :x vector
    select(t, 1)
    select(t, :x)

    # map a function to the :y vector
    select(t, 2 => abs)
    select(t, :y => x -> x > 0 ? x : -x)

    # select the table of :x and :z
    select(t, (:x, :z))
    select(t, r"(x|z)")

    # map a function to the table of :x and :y
    select(t, (:x, :y) => row -> row[1] + row[2])
    select(t, (1, :y) => row -> row.x + row.y)
"""
function select(t::AbstractIndexedTable, which)
    ColDict(t)[which]
end

# optimization
@inline function select(t::IndexedTable, which::Union{Symbol, Int})
    getfield(columns(t), which)
end

"""
    selectkeys(x::NDSparse, sel)

Return an `NDSparse` with a subset of keys.
"""
function selectkeys(x::NDSparse, which; kwargs...)
    ndsparse(rows(keys(x), which), values(x); kwargs...)
end

"""
    selectvalues(x::NDSparse, sel)

Return an `NDSparse` with a subset of values
"""
function selectvalues(x::NDSparse, which; presorted=true, copy=false, kwargs...)
    ndsparse(keys(x), rows(values(x), which); presorted=presorted, copy=copy, kwargs...)
end

"""
    reindex(t::IndexedTable, by)
    reindex(t::IndexedTable, by, select)

Reindex table `t` with new primary key `by`, optionally keeping a subset of columns via
`select`.  For [`NDSparse`](@ref), use [`selectkeys`](@ref).

# Example

    t = table([2,1],[1,3],[4,5], names=[:x,:y,:z], pkey=(1,2))

    t2 = reindex(t, (:y, :z))

    pkeynames(t2)
"""
function reindex end

function reindex(T::Type, t, by, select; kwargs...)
    if isa(by, SpecialSelector)
        return reindex(T, t, lowerselection(t, by), select; kwargs...)
    end
    if !isa(by, Tuple)
        return reindex(T, t, (by,), select; kwargs...)
    end
    if T <: IndexedTable && !isa(select, Tuple) && !isa(select, SpecialSelector)
        return reindex(T, t, by, (select,); kwargs...)
    end
    perm = sortpermby(t, by)
    if isa(perm, Base.OneTo)
        convert(T, rows(t, by), rows(t, select); presorted=true, copy=false, kwargs...)
    else
        convert(T, rows(t, by)[perm], rows(t, select)[perm]; presorted=true, copy=false, kwargs...)
    end
end

function reindex(t::IndexedTable, by=pkeynames(t), select=excludecols(t, by); kwargs...)
    reindex(collectiontype(t), t, by, select; kwargs...)
end

function reindex(t::NDSparse, by=pkeynames(t), select=valuenames(t); kwargs...)
    reindex(collectiontype(t), t, by, select; kwargs...)
end

canonname(t, x::Symbol) = x
canonname(t, x::Int) = colnames(t)[colindex(t, x)]

"""
    map(f, t::IndexedTable; select)

Apply `f` to every item in `t` selected by `select` (see also the [`select`](@ref) function).
Returns a new table if `f` returns a tuple or named tuple.  If not, returns a vector.

# Examples

    t = table([1,2], [3,4], names=[:x, :y])

    polar = map(p -> (r = hypot(p.x, p.y), θ = atan(p.y, p.x)), t)

    back2t = map(p -> (x = p.r * cos(p.θ), y = p.r * sin(p.θ)), polar)
"""
function map(f, t::AbstractIndexedTable; select=nothing) end

function map(f, t::Dataset; select=nothing, copy=false, kwargs...)
    if isa(f, Tup) && select === nothing
        select = colnames(t)
    elseif select === nothing
        select = valuenames(t)
    end

    x = map_rows(f, rows(t, select))
    isa(x, Columns) ? table(x; copy=false, kwargs...) : x
end

missing_indxs(v::AbstractVector) = findall(!_ismissing, v)

function missing_indxs(t::StructArray)
    indxs = collect(1:length(t))
    for vec in getfield(t, :fieldarrays)
        filter!(i -> !_ismissing(vec[i]), indxs)
    end
    indxs
end

"""
    dropmissing(t        )
    dropmissing(t, select)

Drop rows of table `t` which contain missing values (either `Missing` or `DataValue`),
optionally only using the columns in `select`.  Column types will be converted to
non-missing types.  For example:

- `Vector{Union{Int, Missing}}` -> `Vector{Int}`
- `DataValueArray{Int}` -> Vector{Int}

# Example

    t = table([0.1,0.5,missing,0.7], [2,missing,4,5], [missing,6,missing,7], names=[:t,:x,:y])
    dropmissing(t)
    dropmissing(t, (:t, :x))
"""
function dropmissing(t::IndexedTable, sel = All())
    selection = lowerselection(t, sel)
    indxs = missing_indxs(rows(t, selection))
    t2 = subtable(t, indxs)
    d = ColDict(t2)
    for s in selection
        T = eltype(d[s])
        d[s] = convert(Vector{missingtype2type(T)}, d[s])
    end

    table(d[], copy=false, perms=t.perms, presorted=true)
end

dropmissing(t::NDSparse, sel=All()) = ndsparse(dropmissing(table(t), sel))

Base.@deprecate_binding dropna dropmissing

"""
    convertmissing(tbl, missingtype)

Convert the missing value representation in `tbl` to be of type `missingtype`.

# Example

    using IndexedTables, DataValues
    t = table([1,2,missing], [1,missing,3])
    IndexedTables.convertmissing(t, DataValue)
"""
function convertmissing(t::IndexedTable, ::Type{Missing})
    d = ColDict(t)
    for (k, v) in pairs(d)
        T = eltype(v)
        if T <: DataValue
            indxs = findall(!isna, v)
            y = Vector{Union{Missing, missingtype2type(T)}}(missing, length(v))
            y[indxs] = get.(v[indxs])
            d[k] = y
        end
    end
    subtable(d[], 1:length(t))
end
function convertmissing(t::IndexedTable, ::Type{DataValue})
    d = ColDict(t)
    for (k, v) in pairs(d)
        T = eltype(v)
        if Missing <: T
            indxs = findall(!ismissing, v)
            y = DataValueArray(Vector{Base.nonmissingtype(T)}(undef, length(v)), trues(length(v)))
            y[indxs] = v[indxs]
            d[k] = y
        end
    end
    subtable(d[], 1:length(t))
end
convertmissing(t::NDSparse, Typ) = ndsparse(convertmissing(table(t), Typ))

"""
    filter(f, t::Union{IndexedTable, NDSparse}; select)

Iterate over `t` and Return the rows for which `f(row)` returns true.  `select` determines
the rows that are given as arguments to `f` (see [`select`](@ref)).

`f` can also be a tuple of `column => function` pairs.  Returned rows will be those for
which all conditions are true.

# Example

    # filter iterates over ROWS of a IndexedTable
    t = table(rand(100), rand(100), rand(100), names = [:x, :y, :z])
    filter(r -> r.x + r.y + r.z < 1, t)

    # filter iterates over VALUES of an NDSparse
    x = ndsparse(1:100, randn(100))
    filter(val -> val > 0, x)
"""
function Base.filter(fn, t::Dataset; select=valuenames(t))
    x = rows(t, select)
    indxs = findall(fn, x)
    subtable(t, indxs, presorted=true)
end

function Base.filter(pred::Tuple, t::Dataset; select=nothing)
    indxs = [1:length(t);]
    x = select === nothing ? t : rows(t, select)
    for p in pred
        if isa(p, Pair)
            c, f = p
            filt_by_col!(f, rows(x, c), indxs)
        else
            filt_by_col!(p, x, indxs)
        end
    end
    subtable(t, indxs, presorted=true)
end

filt_by_col!(f, col, indxs) = filter!(i->f(col[i]), indxs)

function Base.filter(pred::Pair, t::Dataset; select=nothing)
    filter((pred,), t, select=select)
end

# We discard names of fields in a named tuple. keeps it consistent
# with map and reduce, we don't select using those
function Base.filter(pred::NamedTuple, t::Dataset; select=nothing)
    filter(astuple(pred), t, select=select)
end

ext: jl


text:
module LibMagick

using Color, FixedPointNumbers, ..ColorTypes

import Base: error, size
import Images: @Dict

export MagickWand,
    constituteimage,
    exportimagepixels!,
    getblob,
    getimagealphachannel,
    getimagecolorspace,
    getimagedepth,
    getnumberimages,
    importimagepixels,
    readimage,
    resetiterator,
    setimagecolorspace,
    setimagecompression,
    setimagecompressionquality,
    setimageformat,
    writeimage

# Find the library
depsfile = Pkg.dir("Images","deps","deps.jl")
versionfile = Pkg.dir("Images","deps","versioninfo.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("Images not properly installed. Please run Pkg.build(\"Images\") then restart Julia.")
end
if isfile(versionfile)
    include(versionfile)
end
const have_imagemagick = isdefined(:libwand)

# Initialize the library
function init()
    global libwand
    if have_imagemagick
        eval(:(ccall((:MagickWandGenesis, $libwand), Void, ())))
    else
        warn("ImageMagick utilities not found. Install for more file format support.")
    end
end

# Constants
# Storage types
const CHARPIXEL = 1
const DOUBLEPIXEL = 2
const FLOATPIXEL = 3
const INTEGERPIXEL = 4
const SHORTPIXEL = 7
IMStorageTypes = Union(Uint8, Uint16, Uint32, Float32, Float64)
storagetype(::Type{Uint8}) = CHARPIXEL
storagetype(::Type{Uint16}) = SHORTPIXEL
storagetype(::Type{Uint32}) = INTEGERPIXEL
storagetype(::Type{Float32}) = FLOATPIXEL
storagetype(::Type{Float64}) = DOUBLEPIXEL
storagetype{T<:Ufixed}(::Type{T}) = storagetype(FixedPointNumbers.rawtype(T))
storagetype{CV<:ColorValue}(::Type{CV}) = storagetype(eltype(CV))
storagetype{CV<:AbstractAlphaColorValue}(::Type{CV}) = storagetype(eltype(CV))

# Channel types
type ChannelType
    value::Uint32
end
const UndefinedChannel = ChannelType(0x00000000)
const RedChannel = ChannelType(0x00000001)
const GrayChannel = ChannelType(0x00000001)
const CyanChannel = ChannelType(0x00000001)
const GreenChannel = ChannelType(0x00000002)
const MagentaChannel = ChannelType(0x00000002)
const BlueChannel = ChannelType(0x00000004)
const YellowChannel = ChannelType(0x00000004)
const AlphaChannel = ChannelType(0x00000008)
const MatteChannel = ChannelType(0x00000008)
const OpacityChannel = ChannelType(0x00000008)
const BlackChannel = ChannelType(0x00000020)
const IndexChannel = ChannelType(0x00000020)
const CompositeChannels = ChannelType(0x0000002F)
const TrueAlphaChannel = ChannelType(0x00000040)
const RGBChannels = ChannelType(0x00000080)
const GrayChannels = ChannelType(0x00000080)
const SyncChannels = ChannelType(0x00000100)
const AllChannels = ChannelType(0x7fffffff)
const DefaultChannels = ChannelType( (AllChannels.value | SyncChannels.value) &~ OpacityChannel.value )

# Image type
const IMType = ["BilevelType", "GrayscaleType", "GrayscaleMatteType", "PaletteType", "PaletteMatteType", "TrueColorType", "TrueColorMatteType", "ColorSeparationType", "ColorSeparationMatteType", "OptimizeType", "PaletteBilevelMatteType"]
const IMTypedict = Dict([(IMType[i], i) for i = 1:length(IMType)])

const CStoIMTypedict = @Dict("Gray" => "GrayscaleType", "GrayAlpha" => "GrayscaleMatteType", "RGB" => "TrueColorType", "ARGB" => "TrueColorMatteType", "CMYK" => "ColorSeparationType")

# Colorspace
const IMColorspace = ["RGB", "Gray", "Transparent", "OHTA", "Lab", "XYZ", "YCbCr", "YCC", "YIQ", "YPbPr", "YUV", "CMYK", "sRGB"]
const IMColordict = Dict([(IMColorspace[i], i) for i = 1:length(IMColorspace)])

function nchannels(imtype::String, cs::String, havealpha = false)
    n = 3
    if beginswith(imtype, "Grayscale") || beginswith(imtype, "Bilevel")
        n = 1
        cs = havealpha ? "GrayAlpha" : "Gray"
    elseif cs == "CMYK"
        n = 4
    else
        cs = havealpha ? "ARGB" : "RGB" # only remaining variants supported by exportimagepixels
    end
    n + havealpha, cs
end

# channelorder = ["Gray" => "I", "GrayAlpha" => "IA", "RGB" => "RGB", "ARGB" => "ARGB", "RGBA" => "RGBA", "CMYK" => "CMYK"]

# Compression
const NoCompression = 1

type MagickWand
    ptr::Ptr{Void}
end

function MagickWand()
    wand = MagickWand(ccall((:NewMagickWand, libwand), Ptr{Void}, ()))
    finalizer(wand, destroymagickwand)
    wand
end

destroymagickwand(wand::MagickWand) = ccall((:DestroyMagickWand, libwand), Ptr{Void}, (Ptr{Void},), wand.ptr)

type PixelWand
    ptr::Ptr{Void}
end

function PixelWand()
    wand = PixelWand(ccall((:NewPixelWand, libwand), Ptr{Void}, ()))
    finalizer(wand, destroypixelwand)
    wand
end

destroypixelwand(wand::PixelWand) = ccall((:DestroyPixelWand, libwand), Ptr{Void}, (Ptr{Void},), wand.ptr)

const IMExceptionType = Array(Cint, 1)
function error(wand::MagickWand)
    pMsg = ccall((:MagickGetException, libwand), Ptr{Uint8}, (Ptr{Void}, Ptr{Cint}), wand.ptr, IMExceptionType)
    msg = bytestring(pMsg)
    relinquishmemory(pMsg)
    error(msg)
end
function error(wand::PixelWand)
    pMsg = ccall((:PixelGetException, libwand), Ptr{Uint8}, (Ptr{Void}, Ptr{Cint}), wand.ptr, IMExceptionType)
    msg = bytestring(pMsg)
    relinquishmemory(pMsg)
    error(msg)
end

function getsize(buffer, channelorder)
    if channelorder == "I"
        return size(buffer, 1), size(buffer, 2), size(buffer, 3)
    else
        return size(buffer, 2), size(buffer, 3), size(buffer, 4)
    end
end
getsize{C<:Union(ColorValue,AbstractAlphaColorValue)}(buffer::AbstractArray{C}, channelorder) = size(buffer, 1), size(buffer, 2), size(buffer, 3)

colorsize(buffer, channelorder) = channelorder == "I" ? 1 : size(buffer, 1)
colorsize{C<:Union(ColorValue,AbstractAlphaColorValue)}(buffer::AbstractArray{C}, channelorder) = 1

bitdepth{C<:ColorType}(buffer::AbstractArray{C}) = 8*eltype(C)
bitdepth{T}(buffer::AbstractArray{T}) = 8*sizeof(T)

# colorspace is included for consistency with constituteimage, but it is not used
function exportimagepixels!{T}(buffer::AbstractArray{T}, wand::MagickWand,  colorspace::ASCIIString, channelorder::ASCIIString; x = 0, y = 0)
    cols, rows, nimages = getsize(buffer, channelorder)
    ncolors = colorsize(buffer, channelorder)
    p = pointer(buffer)
    for i = 1:nimages
        status = ccall((:MagickExportImagePixels, libwand), Cint, (Ptr{Void}, Cssize_t, Cssize_t, Csize_t, Csize_t, Ptr{Uint8}, Cint, Ptr{Void}), wand.ptr, x, y, cols, rows, channelorder, storagetype(T), p)
        status == 0 && error(wand)
        nextimage(wand)
        p += sizeof(T)*cols*rows*ncolors
    end
    buffer
end

# function importimagepixels{T}(buffer::AbstractArray{T}, wand::MagickWand, colorspace::ASCIIString; x = 0, y = 0)
#     cols, rows = getsize(buffer, colorspace)
#     status = ccall((:MagickImportImagePixels, libwand), Cint, (Ptr{Void}, Cssize_t, Cssize_t, Csize_t, Csize_t, Ptr{Uint8}, Cint, Ptr{Void}), wand.ptr, x, y, cols, rows, channelorder[colorspace], storagetype(T), buffer)
#     status == 0 && error(wand)
#     nothing
# end

function constituteimage{T<:Unsigned}(buffer::AbstractArray{T}, wand::MagickWand, colorspace::ASCIIString, channelorder::ASCIIString; x = 0, y = 0)
    cols, rows, nimages = getsize(buffer, channelorder)
    ncolors = colorsize(buffer, channelorder)
    p = pointer(buffer)
    depth = bitdepth(buffer)
    for i = 1:nimages
        status = ccall((:MagickConstituteImage, libwand), Cint, (Ptr{Void}, Cssize_t, Cssize_t, Ptr{Uint8}, Cint, Ptr{Void}), wand.ptr, cols, rows, channelorder, storagetype(T), p)
        status == 0 && error(wand)
        setimagecolorspace(wand, colorspace)
        status = ccall((:MagickSetImageDepth, libwand), Cint, (Ptr{Void}, Csize_t), wand.ptr, depth)
        status == 0 && error(wand)
        p += sizeof(T)*cols*rows*ncolors
    end
    nothing
end

function getblob(wand::MagickWand, format::String)
    setimageformat(wand, format)
    len = Array(Csize_t, 1)
    ptr = ccall((:MagickGetImagesBlob, libwand), Ptr{Uint8}, (Ptr{Void}, Ptr{Csize_t}), wand.ptr, len)
    blob = pointer_to_array(ptr, int(len[1]))
    finalizer(blob, relinquishmemory)
    blob
end

function pingimage(wand::MagickWand, filename::String)
    status = ccall((:MagickPingImage, libwand), Cint, (Ptr{Void}, Ptr{Uint8}), wand.ptr, filename)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, filename::String)
    status = ccall((:MagickReadImage, libwand), Cint, (Ptr{Void}, Ptr{Uint8}), wand.ptr, filename)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, stream::IO)
    status = ccall((:MagickReadImageFile, libwand), Cint, (Ptr{Void}, Ptr{Void}), wand.ptr, CFILE(stream))
    status == 0 && error(wand)
    nothing
end

function writeimage(wand::MagickWand, filename::String)
    status = ccall((:MagickWriteImages, libwand), Cint, (Ptr{Void}, Ptr{Uint8}, Cint), wand.ptr, filename, true)
    status == 0 && error(wand)
    nothing
end

function size(wand::MagickWand)
    height = ccall((:MagickGetImageHeight, libwand), Csize_t, (Ptr{Void},), wand.ptr)
    width = ccall((:MagickGetImageWidth, libwand), Csize_t, (Ptr{Void},), wand.ptr)
    return int(width), int(height)
end

getnumberimages(wand::MagickWand) = int(ccall((:MagickGetNumberImages, libwand), Csize_t, (Ptr{Void},), wand.ptr))

nextimage(wand::MagickWand) = ccall((:MagickNextImage, libwand), Cint, (Ptr{Void},), wand.ptr) == 1

resetiterator(wand::MagickWand) = ccall((:MagickResetIterator, libwand), Void, (Ptr{Void},), wand.ptr)

newimage(wand::MagickWand, cols::Integer, rows::Integer, pw::PixelWand) = ccall((:MagickNewImage, libwand), Cint, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}), wand.ptr, cols, rows, pw.ptr) == 0 && error(wand)

# test whether image has an alpha channel
getimagealphachannel(wand::MagickWand) = ccall((:MagickGetImageAlphaChannel, libwand), Cint, (Ptr{Void},), wand.ptr) == 1

# # get number of colors in the image
# magickgetimagecolors(wand::MagickWand) = ccall((:MagickGetImageColors, libwand), Csize_t, (Ptr{Void},), wand.ptr)

# get the type
function getimagetype(wand::MagickWand)
    t = ccall((:MagickGetImageType, libwand), Cint, (Ptr{Void},), wand.ptr)
    # Apparently the following is necessary, because the type is "potential"
    ccall((:MagickSetImageType, libwand), Void, (Ptr{Void}, Cint), wand.ptr, t)
    1 <= t <= length(IMType) || error("Image type ", t, " not recognized")
    IMType[t]
end

# get the colorspace
function getimagecolorspace(wand::MagickWand)
    cs = ccall((:MagickGetImageColorspace, libwand), Cint, (Ptr{Void},), wand.ptr)
    1 <= cs <= length(IMColorspace) || error("Colorspace ", cs, " not recognized")
    IMColorspace[cs]
end

# set the colorspace
function setimagecolorspace(wand::MagickWand, cs::ASCIIString)
    status = ccall((:MagickSetImageColorspace, libwand), Cint, (Ptr{Void},Cint), wand.ptr, IMColordict[cs])
    status == 0 && error(wand)
    nothing
end

# set the compression
function setimagecompression(wand::MagickWand, compression::Integer)
    status = ccall((:MagickSetImageCompression, libwand), Cint, (Ptr{Void},Cint), wand.ptr, int32(compression))
    status == 0 && error(wand)
    nothing
end

function setimagecompressionquality(wand::MagickWand, quality::Integer)
    0 < quality <= 100 || error("quality setting must be in the (inclusive) range 1-100.\nSee http://www.imagemagick.org/script/command-line-options.php#quality for details")
    status = ccall((:MagickSetImageCompressionQuality, libwand), Cint, (Ptr{Void}, Cint), wand.ptr, quality)
    status == 0 && error(wand)
    nothing
end

# set the image format
function setimageformat(wand::MagickWand, format::ASCIIString)
    status = ccall((:MagickSetImageFormat, libwand), Cint, (Ptr{Void}, Ptr{Uint8}), wand.ptr, format)
    status == 0 && error(wand)
    nothing
end

# get the pixel depth
getimagedepth(wand::MagickWand) = int(ccall((:MagickGetImageDepth, libwand), Csize_t, (Ptr{Void},), wand.ptr))

# pixel depth for given channel type
getimagechanneldepth(wand::MagickWand, channelType::ChannelType) = int(ccall((:MagickGetImageChannelDepth, libwand), Csize_t, (Ptr{Void},Uint32), wand.ptr, channelType.value ))

pixelsetcolor(wand::PixelWand, colorstr::ByteString) = ccall((:PixelSetColor, libwand), Csize_t, (Ptr{Void},Ptr{Uint8}), wand.ptr, colorstr) == 0 && error(wand)

relinquishmemory(p) = ccall((:MagickRelinquishMemory, libwand), Ptr{Uint8}, (Ptr{Uint8},), p)

# get library information
# If you pass in "*", you get the full list of options
function queryoptions(pattern::String)
    nops = Cint[0]
    pops = ccall((:MagickQueryConfigureOptions, libwand), Ptr{Ptr{Uint8}}, (Ptr{Uint8}, Ptr{Cint}), pattern, nops)
    ret = Array(ASCIIString, nops[1])
    for i = 1:nops[1]
        ret[i] = bytestring(unsafe_load(pops, i))
    end
    ret
end

# queries the value of a particular option
function queryoption(option::String)
    p = ccall((:MagickQueryConfigureOption, libwand), Ptr{Uint8}, (Ptr{Uint8},), option)
    bytestring(p)
end

end

LibMagick.have_imagemagick


text:
# This scripts trains several feedforward and feedback network architectures on MNIST.
#
# These architectures include the ones trained in:
# Spoerer, McClure & Kriegeskorte (2017). Recurrent Convolutional Neural Networks:
# A Better Model of Biological Object Recognition, Frontiers in Psychology 8, 1551.
#
# It also includes additional architectures with more complex feedback operations.
# For a list of all architectures, see ModelFactory.jl
usegpu = false

using Base.Iterators: partition
using Flux, Statistics
using Flux: onehotbatch, onecold, crossentropy
using MLDatasets
using BSON: @save
using FeedbackNets
if usegpu
    using CuArrays
end

include("DataHelpers.jl")
using .DataHelpers

# load MNIST
@info("Preparing data ...")
images, labels = MNIST.traindata()
images = DataHelpers.standardize(Float32.(images))
trainimgs, trainlbls = DataHelpers.makebatches(images[:, :, 1:55000], labels[1:55000])
valimgs, vallbls = DataHelpers.makebatches(images[:, :, 55001:end], labels[55001:end])

if usegpu
    trainimgs = gpu.(trainimgs)
    trainlbls = gpu.(trainlbls)
    valimgs = gpu.(valimgs)
    vallbls = gpu.(vallbls)
end

# construct models
@info("Building models ...")
basenet = forwardconvnet(basenetparams)
featnet = forwardconvnet(featnetparams)
kernnet = forwardconvnet(kernnetparams)
slimnet = forwardconvnet(slimnetparams)
largenet = forwardconvnet(largenetparams)
models = basenet, featnet, kernnet, slimnet, largenet
if usegpu
    models = map(gpu, models)
end

function trainmodel(model, trainset, valset; epochs=100, opt=Nesterov(0.1), verbose=false)
    # helper functions
    loss(x, y) = crossentropy(model(x), y)
    accuracy(x, y) = mean([x == y for (x, y) in zip(onecold(model(x)), onecold(y))])
    # training loop
    verbose ? @info("initial accuracy: $(accuracy(valset...))") : nothing
    for e in 1:epochs
        Flux.train!(loss, params(model), trainset, opt)
        verbose ? @info("accuracy after episode $e: $(accuracy(valset...))") : nothing
    end
end

# do the actual training
for model in models
    trainmodel(model, zip(trainimgs, trainlbls), zip(valimgs, vallbls))
    @save "$(model)_MNIST.bson"
end


text:
using Klara

p = BasicContUnvParameter(:p, logtarget=p::Float64 -> -abs2(p))

model = likelihood_model(p, false)

sampler = SliceSampler(1., 5)

mcrange = BasicMCRange(nsteps=100000, burnin=10000)

v0 = Dict(:p=>3.11)

job = BasicMCJob(model, sampler, mcrange, v0)

run(job)

chain = output(job)

mean(chain)


text:
using OrderedCollections: LittleDict

macro enm(name, varargs...)
	modname = esc(name )
	names = map(string, varargs)
	basetype = Int32
	if isa(name,Expr) && name.head == :(::) && length(name.args) == 2 && isa(name.args[1],Symbol)
		modname = esc(name.args[1])
		basetype = name.args[2]
	end

	ex = quote
		module $modname
		@enum $(:Enm)::$(basetype) $(varargs...)
		__str_to_enum = LittleDict(zip(map(string,instances($(esc(:Enm)))),instances($(esc(:Enm)))))
		$(esc(:Enm))(x::String) = __str_to_enum[x]

		end
	end
	ex.head = :toplevel
	return ex
end


text:
"""
    optimize!(optimizer, nlp::AbstractNLPEvaluator, x0)

Use optimization routine implemented in `optimizer` to optimize
the optimal power flow problem specified in the evaluator `nlp`.
Initial point is specified by `x0`.

Return the solution as a named tuple, with fields
- `status::MOI.TerminationStatus`: Solver's termination status, as specified by MOI
- `minimum::Float64`: final objective
- `minimizer::AbstractVector`: final solution vector, with same ordering as the `Variables` specified in `nlp`.

    optimize!(optimizer, nlp::AbstractNLPEvaluator)

Wrap previous `optimize!` function and pass as initial guess `x0`
the initial value specified when calling `initial(nlp)`.

## Examples

```julia
nlp = ExaPF.ReducedSpaceEvaluator(datafile)
optimizer = Ipopt.Optimizer()
solution = ExaPF.optimize!(optimizer, nlp)

```

## Notes
By default, the optimization routine solves a minimization
problem.

"""
function optimize! end

# MOI-based solver
function build!(optimizer::MOI.AbstractOptimizer, nlp::AbstractNLPEvaluator)
    block_data = MOI.NLPBlockData(nlp)
    u♭, u♯ = ExaPF.bounds(nlp, ExaPF.Variables())
    n = ExaPF.n_variables(nlp)
    u = MOI.add_variables(optimizer, n)
    # Set bounds
    MOI.add_constraints(
        optimizer, u, MOI.LessThan.(u♯),
    )
    MOI.add_constraints(
        optimizer, u, MOI.GreaterThan.(u♭),
    )
    MOI.set(optimizer, MOI.NLPBlock(), block_data)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return u
end

function optimize!(optimizer::MOI.AbstractOptimizer, nlp::AbstractNLPEvaluator, x0)
    u = build!(optimizer, nlp)
    MOI.set(optimizer, MOI.VariablePrimalStart(), u, x0)
    MOI.optimize!(optimizer)
    x_opt = MOI.get(optimizer, MOI.VariablePrimal(), u)
    solution = (
        status=MOI.get(optimizer, MOI.TerminationStatus()),
        minimum=MOI.get(optimizer, MOI.ObjectiveValue()),
        minimizer=x_opt,
    )
    return solution
end
optimize!(optimizer, nlp::AbstractNLPEvaluator) = optimize!(optimizer, nlp, initial(nlp))


text:
# Script to calculate the availabily of wind power over the year
using DataFrames, FeatherFiles, Statistics

AREA = 0.25445
V_MIN = 6.0    # minimal wind speed at average height [m/s]

function calc_average_power(df)
    minpower = df[df.Wind_Speed_at_Hub_Height .> V_MIN, :Power_Generation]
    return AREA * mean(minpower) * calc_rel_availability(df, 0.0) / 100.0
end

function calc_peak_power(df)
    return AREA * maximum(df.Power_Generation)
end

function calc_capacity_factor(df)
    return 100 * calc_average_power(df) / calc_peak_power(df)
end

function calc_rel_availability(df, min_power)
    minpower = df[(df.Power_Generation * AREA .> min_power) .& (df.Wind_Speed_at_Hub_Height .> V_MIN), :Power_Generation]
    return 100 * size(minpower, 1) / size(df, 1)
end

if !isdefined(@__MODULE__, :df) dfw=DataFrame(load("data/windyear.feather")) end

av = calc_average_power(dfw)
println("Average wind power [W]:                      ", round(av, digits=1))
println("Capacity factor [%]:                          ", round(calc_capacity_factor(dfw), digits=1))
println("Percentage of wind power availability    [%]: ", round(calc_rel_availability(dfw, 0.0), digits=1))
println("Percentage of time for more than 50% av. [%]: ", round(calc_rel_availability(dfw, av/2), digits=1))
println("Percentage of time for more than average [%]: ", round(calc_rel_availability(dfw, av), digits=1))


text:
# ---
# title: 1432. Max Difference You Can Get From Changing an Integer
# id: problem1432
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/max-difference-you-can-get-from-changing-an-integer/description/>
# hidden: true
# ---
#
# You are given an integer `num`. You will apply the following steps exactly
# **two** times:
#
#   * Pick a digit `x (0 <= x <= 9)`.
#   * Pick another digit `y (0 <= y <= 9)`. The digit `y` can be equal to `x`.
#   * Replace all the occurrences of `x` in the decimal representation of `num` by `y`.
#   * The new integer **cannot** have any leading zeros, also the new integer **cannot** be 0.
#
# Let `a` and `b` be the results of applying the operations to `num` the first
# and second times, respectively.
#
# Return _the max difference_ between `a` and `b`.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 555
#     Output: 888
#     Explanation: The first time pick x = 5 and y = 9 and store the new integer in a.
#     The second time pick x = 5 and y = 1 and store the new integer in b.
#     We have now a = 999 and b = 111 and max difference = 888
#
#
# **Example 2:**
#
#
#
#     Input: num = 9
#     Output: 8
#     Explanation: The first time pick x = 9 and y = 9 and store the new integer in a.
#     The second time pick x = 9 and y = 1 and store the new integer in b.
#     We have now a = 9 and b = 1 and max difference = 8
#
#
# **Example 3:**
#
#
#
#     Input: num = 123456
#     Output: 820000
#
#
# **Example 4:**
#
#
#
#     Input: num = 10000
#     Output: 80000
#
#
# **Example 5:**
#
#
#
#     Input: num = 9288
#     Output: 8700
#
#
#
#
# **Constraints:**
#
#   * `1 <= num <= 10^8`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Autogenerated wrapper script for libbluray_jll for i686-linux-musl
export bd_info, bd_list_titles, bd_splice, libbluray

using XML2_jll
using FreeType2_jll
using Fontconfig_jll
using libudfread_jll
JLLWrappers.@generate_wrapper_header("libbluray")
JLLWrappers.@declare_executable_product(bd_info)
JLLWrappers.@declare_executable_product(bd_list_titles)
JLLWrappers.@declare_executable_product(bd_splice)
JLLWrappers.@declare_library_product(libbluray, "libbluray.so.2")
function __init__()
    JLLWrappers.@generate_init_header(XML2_jll, FreeType2_jll, Fontconfig_jll, libudfread_jll)
    JLLWrappers.@init_executable_product(
        bd_info,
        "bin/bd_info",
    )

    JLLWrappers.@init_executable_product(
        bd_list_titles,
        "bin/bd_list_titles",
    )

    JLLWrappers.@init_executable_product(
        bd_splice,
        "bin/bd_splice",
    )

    JLLWrappers.@init_library_product(
        libbluray,
        "lib/libbluray.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module AutoDiff

include( "dictionaries.jl" )

import Base.<, Base.<=, Base.==, Base.>=, Base.>
import Base.+, Base.-, Base.*, Base./, Base.^
import Base.cos, Base.sin, Base. tan, Base.acos, Base.asin, Base.atan
import Base.cosh, Base.sinh, Base.tanh, Base.acosh, Base.asinh, Base.atanh
import Base.sqrt, Base.exp, Base.log, Base.expm1
import Base.erf, Base.erfc
import Base.hypot, Base.atan2

import Base.getindex,Base.map
import Applicator.apply_unary, Applicator.apply_binary

using Identity, Scalar, Comparable, DataWrapper, Elementwise, Arithmetic, SpecFunable

export Dual, ArithmeticArray, AbstractComparable, AbstractDataWrapper, AbstractElementwise, data, unit
export Dual, re, du, <, <=,  ==, eq, neq
export +, -, *, /, ^
export apply_unary, apply_binary

abstract AbstractComparableScalar <: AbstractScalar
abstract AbstractDual <: AbstractComparableScalar

diffu( op ) = unary_diff[ op ]
diffb1( op ) = binary1_diff[ op ]
diffb2( op ) = binary2_diff[ op ]

type Dual{ Re, Du } <: AbstractScalar
    re::Re
    du::Du
end

Dual( x ) = Dual( x, additive_identity )

re( x::Dual ) = x.re
du( x::Dual ) = x.du

apply_unary( op, x::Dual ) = Dual( op( re( x )), du( x )*diffu( op )( re( x )))

apply_binary( op, x::Dual, y::Dual ) =
    Dual( op( re( x ), re( y )),
          du( x )*diffb1( op )( re( x ), re( y ))+
          du( y )*diffb2( op )( re( x ), re( y )))

apply_binary( op, x::Dual, y::Number ) =
    Dual( op( re( x ), y ),
          du( x )*diffb1( op )( re( x ), y ))

apply_binary( op, x::Number, y::Dual ) =
    Dual( op( x, re( y )),
          du( y )*diffb2( op )( x, re( y )))

<( x::Dual, y::Dual ) = re( x ) < re( y )
<( x::Dual, y::Number ) = re( x ) < y
<( x::Number , y::Dual ) = x < re( y )

@make_comparable Dual Number
@make_scalar Dual

for f in unary_function_list
    @eval @make_unary_specfun Dual $f
end

for f in binary_function_list
    @eval @make_binary_specfun Dual Number $f
end

end # module AutoDiff

max_stars_repo_name: sabjohnso/AutoDiff

max_issues_repo_name: sabjohnso/AutoDiff

max_forks_repo_name: sabjohnso/AutoDiff


text:
@testset "GR4J.jl" begin
    @testset "Constructors" begin
		model = GR4J()
		model2 = GR4J(model)
		@test model == model2
	end

    @testset "Check param or state values when creating new instance" begin
        @test_throws DomainError GR4J(-350.0, 0.0, 40.0, 0.5)
        @test_throws DomainError GR4J(350.0, 0.0, -40.0, 0.5)
        @test_throws DomainError GR4J(350.0, 0.0, 40.0, -0.5)
        @test_throws DomainError GR4J(350.0, 0.0, 40.0, 0.5, 351.0, 30.0)
        @test_throws DomainError GR4J(350.0, 0.0, 40.0, 0.5, 0.0, 41.0)
    end

    @testset "Check values when a parameter or state value is set" begin
        model = GR4J()
        @test_throws DomainError model.X1 = -350.0
        @test_throws DomainError model.X3 = -0.1
        @test_throws DomainError model.X4 = -0.1
        @test_throws DomainError model.Sp = 500.0
        @test_throws DomainError model.Sr = 50.0
    end

    @testset "Stash and recover a model for hot start" begin
        model = GR4J()
        d = stash(model)
        model2 = GR4J(d)
        @test model == model2
    end

    @testset "Print a model" begin
    	model = GR4J()
    	println(model)
	end

    @testset "Test one step running" begin
        model = GR4J()
        Q1, AET1 = simulate(model, 10.0, 5.0)
        Q2, AET2 = simulate(model, 0.0, 5.0)

        model = GR4J()
        Qs, AETs = simulate(model, [10.0, 0.0], [5.0, 5.0])

        @test Qs[2] == Q2
        @test AETs[2] == AET2
    end

    @testset "Check againt GR4J MS Excel modelling results" begin
        data = dataset("gr4j_sample")

        X1 = 320.0
        X2 = 2.5
        X3 = 69.5
        X4 = 1.5

        Sp0 = 0.60*X1
        Sr0 = 0.70*X3

        model = GR4J(X1, X2,X3, X4, Sp0, Sr0)

        Q, AET = simulate(model, data.P, data.PET)
        dQ = Q - data.Qsim
        @test maximum(abs.(dQ)) < 1E-3
  	end

end


text:
include("workspace.jl")

#make circuit
qpe = Circuit(4)
X!(qpe,4)
H!(qpe,1:3)
λ = π/4

reps = 1
for counting_qubit ∈ 1:3
    for i ∈ 1:reps
        CU!(qpe, λ, counting_qubit, 4)
    end
    reps *= 2
end

function qft_dagger!(C::Circuit, n::Int)
    #Don't forget the swaps!
    for qubit in 1:div(n,2)-1
        Swap!(C, qubit, n-qubit)
    end
    for j in 1:n-1
        for m in 1:j-1
            CU!(C,-π/(2^(j-m)), m, j)
        end
        H!(C,j)
    end
end

qft_dagger!(qpe,4)
data = shots(qpe,10_000)
graph!(data)

##qpe2
qpe2 = Circuit(4)

X!(qpe2,4)
H!(qpe2,1:3)

λ2 = 2*π/3
reps = 1
for counting_qubit ∈ 1:3
    for i ∈ 1:reps
        CU!(qpe2, λ2, counting_qubit, 4)
    end
    reps *= 2
end

qft_dagger!(qpe2,4)
data2 = shots(qpe2,10_000)
graph!(data2)

##qpe3
qpe3 = Circuit(6)

H!(qpe3,1:5)
X!(qpe3,6)

reps = 1
for counting_qubit ∈ 1:5
    for i ∈ 1:reps
        CU!(qpe3, λ2, counting_qubit, 6)
    end
    reps *= 2
end

qft_dagger!(qpe3, 6)
data3 = shots(qpe3,10_000)
graph!(data3)


text:
@testset "descendants" begin
    mtg = read_mtg("files/simple_plant.mtg");
    width_all = [nothing,nothing,1.0,6.0,nothing,7.0]
    @test descendants(mtg, :Width; type = Union{Nothing,Float64}) ==  width_all
    @test descendants(mtg, :Width) ==  width_all

    d = descendants(mtg, :Width, scale = 1)
    @test typeof(d) == Vector{Any}
    @test length(d) == 1
    @test d[1] == width_all[1]
    d_typed = descendants(mtg, :Width, type = Union{Nothing,Float64})
    @test typeof(d_typed) ==  Vector{Union{Nothing,Float64}}
    @test descendants(mtg, :Width, symbol = ("Leaf", "Internode")) == width_all[3:end]

    mtg2 = mtg[1][1][1][2]
    @test descendants(mtg2, :Width, symbol = "Leaf")[1] == width_all[end]

    @test descendants(mtg2, :Width, symbol = ("Leaf", "Internode"), self = true) ==
        width_all[end - 1:end]

    # Using the mutating version:
    @test descendants!(mtg, :Width) == descendants(mtg, :Width)
    @test descendants!(mtg2, :Width, symbol = ("Leaf", "Internode"), self = true) ==
        width_all[end - 1:end]

    clean_cache!(mtg)
    # Get the leaves values:
    @test descendants(mtg, :Width; filter_fun = isleaf) ==  width_all[[end - 2,end]]

    # descendants!(mtg, :Width, symbol = ("Leaf", "Internode"), self = true)
end

# using BenchmarkTools

# @benchmark descendants(mtg, :Width) # 19.000 μs
# @benchmark descendants!(mtg, :Width) # 19.000 μs
# @benchmark descendants(mtg, :Width; type = Union{Nothing,Float64})
# @benchmark descendants(mtg, :Width)


text:
using Parameters
include("modelFunctions.jl")
"""
An instruction is an object used to tell the experiment manager when and how
to run a step of the experiment. all instructions are subtypes of this
abstract type
"""
abstract type AbstractInstruction end

"""
An individual instruction tells the experiment manager to run a step of the
experiment that acts on one or more individuals.
--
An individual instruction can include the field "tag," which takes a Symbol.
When the individual instruction
"""
abstract type IndividualInstruction <: AbstractInstruction end
abstract type ModelInstruction <: AbstractInstruction end

@with_kw struct RandomGenome <: IndividualInstruction
    genomesize::Int
    interactiontype::String
    connectivity::Real
    # filter parameters
    tag::Union{Symbol,Nothing} = nothing
    value::Any = nothing
    exclude::Bool = false
    # associated function
    name::function = generate_genotype
end


text:
# Julia wrapper for header: /usr/include/libavcodec/vdpau.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_alloc_vdpaucontext,
    av_vdpau_hwaccel_get_render2,
    av_vdpau_hwaccel_set_render2,
    av_vdpau_bind_context,
    av_vdpau_get_surface_parameters,
    av_vdpau_alloc_context,
    av_vdpau_get_profile

function av_alloc_vdpaucontext()
    ccall((:av_alloc_vdpaucontext, libavcodec), Ptr{AVVDPAUContext}, ())
end

function av_vdpau_hwaccel_get_render2(arg1)
    ccall((:av_vdpau_hwaccel_get_render2, libavcodec), AVVDPAU_Render2, (Ptr{AVVDPAUContext},), arg1)
end

function av_vdpau_hwaccel_set_render2(arg1, arg2::AVVDPAU_Render2)
    ccall((:av_vdpau_hwaccel_set_render2, libavcodec), Cvoid, (Ptr{AVVDPAUContext}, AVVDPAU_Render2), arg1, arg2)
end

function av_vdpau_bind_context(avctx, device::VdpDevice, get_proc_address, flags::Integer)
    ccall((:av_vdpau_bind_context, libavcodec), Cint, (Ptr{AVCodecContext}, VdpDevice, Ptr{VdpGetProcAddress}, UInt32), avctx, device, get_proc_address, flags)
end

function av_vdpau_get_surface_parameters(avctx, _type, width, height)
    ccall((:av_vdpau_get_surface_parameters, libavcodec), Cint, (Ptr{AVCodecContext}, Ptr{VdpChromaType}, Ptr{UInt32}, Ptr{UInt32}), avctx, _type, width, height)
end

function av_vdpau_alloc_context()
    ccall((:av_vdpau_alloc_context, libavcodec), Ptr{AVVDPAUContext}, ())
end

function av_vdpau_get_profile(avctx, profile)
    ccall((:av_vdpau_get_profile, libavcodec), Cint, (Ptr{AVCodecContext}, Ptr{VdpDecoderProfile}), avctx, profile)
end


text:
function load_iorio_data(req::HTTP.Request)
    request_dictionary = Dict{String, String}()
    if length(req.body) > 0
        request_dictionary = JSON.parse(transcode(String, req.body))
    end
    host = "localhost"
    base_dir = joinpath(PROJECT_ROOT, "data", "iorio")
    force_reload = haskey(request_dictionary, "force_reload") ? request_dictionary["force_reload"] : false
    if Utils.experiment_registered("iorio")
        if force_reload
            proc_id = Utils.unregister_experiment("iorio")
            remotecall_wait(Data.delete_experiment, proc_id, "iorio")
        else
            req.response.status = 200
            req.response.body = create_response(JSON.json(Dict("status" => "success", "message" => "Iorio data set already exists, returning experiment ID. If you need to reload, call with parameter force_reload = true", "experiment_id" => "iorio")))
            return req.response
        end
    end

    # create the experiment
    new_proc_id = Utils.register_experiment("iorio")
    remotecall_wait(Data.create_experiment, new_proc_id, "iorio")

    #import the cell lines
    for f in readdir(joinpath(base_dir, "cell_lines"))
        if !startswith(basename(f), "json") continue end
        filename = joinpath(base_dir, "cell_lines", f)
        cl_cmd = `curl -X POST http://$host:$port/experiments/iorio/data/cell_line -d @$filename`
        run(cl_cmd)
    end

    #import the outcome, use AUC for the time being
    outcome_file = joinpath(PROJECT_ROOT, "data", "iorio", "json_AUCvalues.txt")
    outcome_add_cmd = `curl -X POST http://$host:$port/experiments/iorio/data/outcomes -d @$outcome_file`
    run(outcome_add_cmd)

    # add pathway info
    pw_file = joinpath(PROJECT_ROOT, "data", "json", "pathways_current.json")
    run(`curl -X POST http://$host:$port/experiments/iorio/data/pathways -d @$pw_file`)

    req.response.status = 200
    req.response.body = create_response(JSON.json(Dict("status" => "success", "message" => "Created Iorio data set", "experiment_id" => "iorio")))
    return req.response
end

function load_dream_challenge_data(req::HTTP.Request)
    request_dictionary = Dict{String, String}()
    if length(req.body) > 0
        request_dictionary = JSON.parse(transcode(String, req.body))
    end
    host = "localhost"
    base_dir = joinpath(PROJECT_ROOT, "data", "dream_challenge")
    force_reload = haskey(request_dictionary, "force_reload") ? request_dictionary["force_reload"] : false
    if Utils.experiment_registered("dream_challenge")
        if force_reload
            proc_id = Utils.unregister_experiment("dream_challenge")
            @spawnat proc_id Data.delete_experiment("dream_challenge")
        else
            req.response.status = 200
            req.response.body = create_response(JSON.json(Dict("status" => "success", "message" => "DREAM Challenge data set already exists, returning experiment ID. If you need to reload, call with parameter force_reload = true", "experiment_id" => "dream_challenge")))
            return req.response
        end
    end
    proc_id = Utils.register_experiment("dream_challenge")
    remotecall_wait(Data.load_dream_challenge_data, proc_id)

    req.response.status = 200
    req.response.body = create_response(JSON.json(Dict("status" => "success", "message" => "Created DREAM Challenge data set", "experiment_id" => "dream_challenge")))
    return req.response
end

# function run_dreamchallenge_data(pathways_file::AbstractString, dest_dir::AbstractString="results/")
#     @async start_server(8888)
#     srv_shutdown_cmd = `curl http://localhost:8888/stop`
#     experiment = load_dream_challenge_data()
#
#     # stop_cmd = `curl http://localhost:8888/stop`
#     # run(stop_cmd)
#
#     # do the different data settings
#     # 1) full exome sequencing data
#     #   a) full data, no pathways
#     #   b) variance filtered, no pathways
#     #   c) full data, full pathways
#     #   d) variance filtered, full pathways
#     #   e) full data, subsumed pathways
#     #   f) variance filtered, subsumed pathways
#     # 2) all the above but with exome seq data constant 1
#
#     # 1)
#     do_experiment_stack(joinpath(dest_dir, "full_exome"), experiment, pathways_file)
#
#     # normalize again with constant exome data
#     experiment.is_normalized = false
#     # overwrite function
#     get_measurement_value(d::ExomeSeq) = 1
#
#     # 2)
#     do_experiment_stack(joinpath(dest_dir, "const_exome"), experiment, pathways_file)
#
#     @async run(srv_shutdown_cmd)
#     nothing
# end
#
# function do_experiment_stack(dest_dir::String, experiment::Experiment, pathways_file::String)
#     add_pathways_cmd = `curl -X POST http://localhost:8888/experiments/dream_challenge/pathways -d @$pathways_file`
#     ic = InferenceConfiguration()
#     ic.do_gridsearch = true
#     ic.compute_wpc_index = true
#     # a
#     ic.do_variance_filtering = false
#     dep = create_drug_efficacy_predictor(experiment, ic)
#     ic.target_dir = joinpath(dest_dir, "full_data_no_pathways")
#     run_model(dep, inference_config = ic)
#
#     run(add_pathways_cmd)
#     # c
#     ic.subsume_pathways = false
#     dep = create_drug_efficacy_predictor(experiment, ic)
#     ic.target_dir = joinpath(dest_dir, "full_data_full_pathways")
#     run_model(dep, inference_config = ic)
#
#     # e
#     ic.subsume_pathways = true
#     dep = create_drug_efficacy_predictor(experiment, ic)
#     ic.target_dir = joinpath(dest_dir, "full_data_sub_pathways")
#     run_model(dep, inference_config = ic)
#
#     # remove pathways for experiment
#     empty!(experiment.pathway_information)
#
#     # b
#     ic.do_variance_filtering = true
#     dep = create_drug_efficacy_predictor(experiment, ic)
#     ic.target_dir = joinpath(dest_dir, "var_filter_data_no_pathways")
#     run_model(dep, inference_config = ic)
#
#     # add pathways again
#     run(add_pathways_cmd)
#     # d
#     ic.subsume_pathways = false
#     dep = create_drug_efficacy_predictor(experiment, ic)
#     ic.target_dir = joinpath(dest_dir, "var_filter_data_full_pathways")
#     run_model(dep, inference_config = ic)
#
#     # f
#     ic.subsume_pathways = true
#     dep = create_drug_efficacy_predictor(experiment, ic)
#     ic.target_dir = joinpath(dest_dir, "var_filter_data_sub_pathways")
#     run_model(dep, inference_config = ic)
#
# end


text:
# Some kernels
# Daniel Perry (dperry@cs.utah.edu)
# June 2014

# compute k using a gaussian kernel
# @param X1 - data set
# @param X2 - data set
# @param sigma - sigma used in gaussian kernel
function GaussianKernel!{T}(K::Array{T,2},X1::Array{T,2},X2::Array{T,2}, sigma)
		if !isa(sigma,T)
			sigma = convert(T,sigma)
		end
    #for row = 1:size(X1,1)
		row = 1
		while row <= size(X1,1)
        #for col = 1:size(X2,1)
				col = 1
				while col <= size(X2,1)
            @inbounds tmp = sum((X1[row,:]-X2[col,:]).^2)
            @inbounds K[row,col] = exp(-tmp/(2*sigma^2))
						col += 1
        end
				row += 1
    end
end
function GaussianKernel{T}(X1::Array{T,2},X2::Array{T,2}, sigma)
		if !isa(sigma,T)
			sigma = convert(T,sigma)
		end
    K = zeros(T,size(X1,1),size(X2,1))
		GaussianKernel!(K,X1,X2,sigma)
    return K
end

# compute k using a gaussian kernel
# @param X - data set
# @param sigma - sigma used in gaussian kernel
function GaussianKernel!(K, X, sigma)
    #for row = 1:size(X,1)
		row = 1
		while row <= size(X,1)
        #for col = 1:row
				col = 1
				while col <= row
            tmp = sum((X[row,:]-X[col,:]).^2)
            K[row,col] = exp(-tmp/(2*sigma^2))
            K[col,row] = K[row,col]
						col += 1
        end
				row += 1
    end
end
function GaussianKernel{T}(X::Array{T,2}, sigma)
	if !isa(sigma,T)
		sigma = convert(T,sigma)
	end
  K = zeros(T,size(X,1),size(X,1))
	GaussianKernel!(K,X,sigma)
	return K
end

# compute distance in "input" space from distance in feature space between X_i, X_j
# @param feature dot project of X_i and X_i
# @param feature dot project of X_j and X_j
# @param feature distance between X_i and X_j
# @param sigma kernel parameter
#
# @return "input" distance
function GaussianKernelInverseDistance(Kii,Kjj,feature_distance_squared, sigma)
	-sigma^2 * log( .5 * (Kii + Kjj - feature_distance_squared) )
end

# compute k using an ANOVA kernel
# @param X1 - data set
# @param X2 - data set
# @param sigma,d  - kernel params
function ANOVAKernel(X1,X2, sigma, d)
    K = zeros(size(X1,1),size(X2,1))
    for row = 1:size(X1,1)
        for col = 1:size(X2,1)
            K[row,col] = sum( exp( - sigma .* (X1[row,:]-X2[col,:]).^2 ).^d )
        end
    end
    return K
end
# compute k using an ANOVA kernel
# @param X - data set
# @param sigma,d - kernel params
function ANOVAKernel(X, sigma, d)
    K = zeros(size(X,1),size(X,1))
    for row = 1:size(X,1)
        for col = 1:row
            K[row,col] = sum( exp( - sigma .* (X[row,:]-X[col,:]).^2 ).^d )
            K[col,row] = K[row,col]
        end
    end
    return K
end

# compute k using a linear kernel
# @param X1,X2 - data set
# @param scale,offset,power - polynomial parameters
function PolynomialKernel(X1,X2, scale,offset,power)
  K = (scale .* (X1*X2') .+ offset).^power
end
# compute k using a linear kernel
# @param X - data set
# @param scale,offset,power - polynomial parameters
function PolynomialKernel(X, scale,offset,power)
    K = (scale .* X*X' .+ offset).^power
end

# compute k using a linear kernel
# @param X1,X2 - data set
function LinearKernel(X1,X2)
  K = X1*X2'
end
# compute k using a linear kernel
# @param X - data set
function LinearKernel(X)
    K = X*X'
end


text:
function esf_sum!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  @inbounds for col in 1:n
    for r in 1:col
        row = col - r + 1
        S[row+1] = x[col] * S[row] + S[row+1]
    end
  end
end

"""
    esf_sum(x)

Compute the elementary symmetric functions of order k = 1, ..., n
where n = length(x)

# Examples
```julia-repl
julia> esf_sum([3.5118, .6219, .2905, .8450, 1.8648])
6-element Array{Float64,1}:
  1.0
  7.134
 16.9493
 16.7781
  7.05289
 0.999736
```
"""
function esf_sum(x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  esf_sum!(S,x)
  return S
end

function esf_sum_reg!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  @inbounds for col in 1:n
    for r in 1:col
        row = col - r + 1
        S[row+1] = ((col-row)/col) * S[row+1] + (row/col) * x[col] * S[row]
    end
  end
end

"""
    esf_sum_reg(x)

Compute the elementary symmetric functions of order k = 1, ..., n
where n = length(x). Values are computed regularized by the binomial
coefficient binomial(n, k) to prevent over/under-flow.

# Examples
```julia-repl
julia> esf_sum_reg([3.5118, .6219, .2905, .8450, 1.8648])
6-element Array{Float64,1}:
  1.0
  1.4268
  1.69493
  1.67781
  1.41058
  0.999736
```
"""
function esf_sum_reg(x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  esf_sum_reg!(S,x)
  return S
end

function esf_sum_log!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,-Inf)
  S[1] = zero(T)
  @inbounds for col in 1:n
    for r in 1:col
        row = col - r + 1
        Sr = S[row] + log(x[col])
        Sr1 = S[row+1]
        if (Sr1 > Sr) && (Sr1 > zero(T))
          S[row+1] = Sr1 + log1p(exp(Sr - Sr1))
        elseif (Sr >= Sr1) && (Sr > zero(T))
          S[row+1] = Sr + log1p(exp(Sr1-Sr))
        else
          S[row+1] = log(exp(Sr1) + exp(Sr))
        end
    end
  end
end

function esf_sum_log(x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  esf_sum_log!(S,x)
  return S
end

#Regularized summation algorithm where one input is zeroed out (for computing derivatives)
function esf_sum_reg2!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  adj = 0
  @inbounds for col in 1:n
    if x[col] == 0.0
      adj += 1
      continue
    end
    for r in 1:(col-adj)
      row = (col-adj) - r + 1
      S[row+1] = (col-adj-row)/(col-adj+1) * S[row+1] + (row+1)/(col-adj+1) * x[col] * S[row]
    end
    S[1] *= (col-adj)/(col-adj+1)
  end
end

#Regularized summation algorithm where two inputs are zeroed out (for computing derivatives)
function esf_sum_reg3!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  adj = 0
  @inbounds for col in 1:n
    if x[col] == 0.0
      adj += 1
      continue
    end
    for r in 1:(col-adj)
      row = (col-adj) - r + 1
      S[row+1] = (col-adj-row)/(col-adj+2) * S[row+1] + (row+2)/(col-adj+2) * x[col] * S[row]
    end
    S[1] *= (col-adj) / (col-adj+2)
  end
end

function esf_sum_dervs_1(x::AbstractVector{T}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  P = Array{T,2}(undef,n,n+1)
  esf_sum!(S, x)
  esf_sum_dervs_1!(P, x)
  return S, P
end

function esf_sum_dervs_1!(P::AbstractArray{T,2}, x::AbstractVector{T}) where T <: Real
  n = length(x)
  xj=zero(T)
  @inbounds for j in 1:n
    xj = x[j]
    x[j] = zero(T)
    @views esf_sum!(P[j,:], x)
    x[j] = xj
  end
end
function esf_sum_dervs_1_reg(x::AbstractVector{T}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  P = Array{T,2}(undef,n,n+1)
  esf_sum_reg!(S, x)
  esf_sum_dervs_1_reg!(P, x)
  return S, P
end

function esf_sum_dervs_1_reg!(P::AbstractArray{T,2}, x::AbstractVector{T}) where T <: Real
  n = length(x)
  xj=zero(T)
  @inbounds for j in 1:n
    xj = x[j]
    x[j] = zero(T)
    @views esf_sum_reg2!(P[j,:], x)
    x[j] = xj
  end
end

function esf_sum_dervs_2(x::AbstractVector{T}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  H = Array{T,3}(undef,n,n,n+1)
  esf_sum!(S, x)
  esf_sum_dervs_2!(H, x)
  return S, H
end

function esf_sum_dervs_2!(H::AbstractArray{T,3}, x::AbstractVector{T}) where T <: Real
  n = length(x)
  xj=zero(T)
  @inbounds for j in 1:n
    xj = x[j]
    x[j] = zero(T)
    for k in j:n
      xk = x[k]
      x[k] = zero(T)
      @views esf_sum!(H[j,k,:], x)
      H[k,j,:] .= H[j,k,:]
      x[k] = xk
    end
    x[j] = xj
  end
end

function esf_sum_dervs_2_reg(x::AbstractVector{T}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  H = Array{T,3}(undef,n,n,n+1)
  esf_sum_reg!(S, x)
  esf_sum_dervs_2_reg!(H, x)
  return S, H
end

function esf_sum_dervs_2_reg!(H::AbstractArray{T,3}, x::AbstractVector{T}) where T <: Real
  n = length(x)
  xj=zero(T)
  @inbounds for j in 1:n
    xj = x[j]
    x[j] = zero(T)
    @views esf_sum_reg2!(H[j,j,:], x)
    for k in j+1:n
      xk = x[k]
      x[k] = zero(T)
      @views esf_sum_reg3!(H[j,k,:], x)
      H[k,j,:] .= H[j,k,:]
      x[k] = xk
    end
    x[j] = xj
  end
end

function esf_dc_fft!(S::AbstractArray{T,1}, tempS::AbstractArray{T,2},
                       x::AbstractArray{T,1}, si::AbstractArray{T,1},
                       group_sizes::AbstractArray{D,1},
                       group_start_idx::AbstractArray{D,1}) where {T <: Real, D <: Integer}
  n = length(x)
  M = size(tempS)[2]

  #convolve initial subsets
  @inbounds for g in 1:M
    @views esf_sum!(tempS[1:(group_sizes[g]+1),g],
                    x[group_start_idx[g]:(group_start_idx[g]+group_sizes[g]-1)])
    group_sizes[g] += 1
  end

  while M > 1
    next_avail_col = 1
    @inbounds for g in 1:2:M
      m = group_sizes[g] + group_sizes[g+1] - 1
      si .= zero(T)
      @views _filt!(S[1:m], tempS[1:group_sizes[g],g], tempS[1:m,g+1], si[1:(group_sizes[g]-1)])
      @views copyto!(tempS[1:m,next_avail_col], S[1:m])
      group_sizes[next_avail_col] = m
      next_avail_col += 1
    end
    M = div(M,2)
  end
end

function esf_dc_fft(x::AbstractArray{T,1}, k::D=2) where {T <: Real, D <: Integer}
  n = length(x)
  k = min(floor(D, log2(n)), k)
  M = 2^k
  L = n/M
  r = rem(n,M) / M

  group_sizes = [fill(floor(D, L), D(M*(1-r))); fill(ceil(D, L), D(M*r))]
  group_start_idx = cumsum(group_sizes) .- (group_sizes .- 1)

  S = Vector{T}(undef,n+1)
  tempS = zeros(T, n+1,M)
  si = zeros(T, n)

  esf_dc_fft!(S, tempS, x, si, group_sizes, group_start_idx)
  return S
end

function esf_dc_fft_reg!(S::AbstractArray{T,1}, tempS::AbstractArray{T,2},
                           x::AbstractArray{T,1}, si::AbstractArray{T,1},
                           group_sizes::AbstractArray{D,1},
                           group_start_idx::AbstractArray{D,1}) where {T <: Real, D <: Integer}
  n = length(x)
  M = size(tempS)[2]
  tempS .= zero(T)

  #convolve initial subsets
  @inbounds for g in 1:M
    @views esf_sum_reg!(tempS[1:(group_sizes[g]+1),g],
                        x[group_start_idx[g]:(group_start_idx[g]+group_sizes[g]-1)])
    group_sizes[g] += 1
  end

  while M > 1
    next_avail_col = 1
    @inbounds for g in 1:2:M
      m = group_sizes[g] + group_sizes[g+1] - 1
      si .= zero(T)
      @views _filt_reg!(S[1:m], tempS[1:group_sizes[g],g], tempS[1:m,g+1], si[1:(group_sizes[g]-1)], group_sizes[g]-1,group_sizes[g+1]-1)
      @views copyto!(tempS[1:m,next_avail_col], S[1:m])
      group_sizes[next_avail_col] = m
      next_avail_col += 1
    end
    M = div(M,2)
  end
end

function esf_dc_fft_reg(x::AbstractArray{T,1}, k::D=2) where {T <: Real, D <: Integer}
  n = length(x)
  k = min(floor(D, log2(n)), k)
  M = 2^k
  L = n/M
  r = rem(n,M) / M

  group_sizes = [fill(floor(D, L), D(M*(1-r))); fill(ceil(D, L), D(M*r))]
  group_start_idx = cumsum(group_sizes) .- (group_sizes .- 1)

  S = Vector{T}(undef,n+1)
  tempS = zeros(T, n+1,M)
  si = zeros(T, n)

  esf_dc_fft_reg!(S, tempS, x, si, group_sizes, group_start_idx)
  return S
end

function esf_diff!(S::AbstractArray{T,1}, P::AbstractArray{T,2}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S, zero(T))
  S[1] = one(T)
  S[2] = sum(x)
  P[:,1] .= one(T)
  @inbounds for i in 2:n
    for j in 1:n
      P[j,i] = S[i] - x[j] * P[j,i-1]
      S[i+1] = S[i+1] + x[j] * P[j,i]
    end
    S[i+1] /= i
  end
end
"""
    esf_diff(x)

Compute the elementary symmetric functions of order k = 1,...,n
where n = length(x) using the difference algorithm. Also computes
the matrix of first derivatives.

# Examples
```julia-repl
julia> esf_diff([3.5118, .6219, .2905, .8450, 1.8648])[1]
6-element Array{Float64,1}:
  1.0
  7.134
 16.9493
 16.7781
  7.05289
  0.999736

julia> esf_diff([3.5118, .6219, .2905, .8450, 1.8648])[2]
5×5 Array{Float64,2}:
  1.0  3.6222   4.22884   1.92728  0.284679
  1.0  6.5121  12.8994    8.75598  1.60755
  1.0  6.8435  14.9612   12.4319   3.44143
  1.0  6.289   11.6351    6.94648  1.18312
  1.0  5.2692   7.12328   3.49463  0.536109
```
"""
function esf_diff(x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  P = Matrix{T}(undef,n,n)
  esf_diff!(S,P,x)
  return S, P
end

function esf_diff_onlyP!(S::AbstractArray{T,1}, P::AbstractArray{T,2}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  P[:,1] .= one(T)
  @inbounds for i in 2:n
    for j in 1:n
      P[j,i] = S[i] - x[j] * S[j,i-1]
    end
  end
end

function esf_diff_reg!(S::AbstractArray{T,1}, P::AbstractArray{T,2}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S, zero(T))
  S[1] = one(T)
  S[2] = sum(x)/n
  P[:,1] .= one(T)/n
  @inbounds for i in 2:n
    for j in 1:n
      P[j,i] = (S[i] - x[j] * P[j,i-1]) * (i / (n+1-i))
      S[i+1] = S[i+1] + x[j] * P[j,i]
    end
    S[i+1] /= i
  end
end

function esf_diff_reg(x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  P = Matrix{T}(undef,n,n)
  esf_diff_reg!(S,P,x)
  return S, P
end

function esf_diff_updown(x::AbstractArray{T,1}) where T<:Real
  n = length(x)
  S = Vector{T}(undef,n+1)
  P = Matrix{T}(undef,n,n+1)
  err = zero(T)
  esf_diff_updown!(S,P,err,x)
  return S, P[:,1:end-1], err
end

function esf_diff_updown!(S::AbstractArray{T,1}, P::AbstractArray{T,2},
                          fw_bk_err::T,
                          x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S, zero(T))
  S[1] = one(T)
  S[2] = sum(x)
  S[n+1] = prod(x)
  P[:,1] .= one(T)
  P[:,end] .= zero(T)
  mid = div(n, 2) + 1
  @inbounds for i in 2:(mid-1)
    for j in 1:n
      P[j,i] = S[i] - x[j] * P[j,i-1]
      S[i+1] = S[i+1] + x[j] * P[j,i]
    end
    S[i+1] /= i
  end
  S_mid_est = S[mid]
  S[mid] = zero(T)
  @inbounds for ii in mid:n
    i = n - ii + mid
    for j in 1:n
      P[j,i] = (S[i+1] - P[j,i+1]) / x[j]
      S[i] += P[j,i]
    end
    S[i] /= (n - (i-1))
  end
  fw_bk_err = abs(one(T) - S_mid_est / S[mid])
end


text:
"""
Construct a Dict{Symbol, Any} that maps the arguments of the function
`ast` to the corresponding value in Vector `args`.
"""
function build_arg_name_map(args::Vector, ast::Expr)
    @assert is_function(ast)
    map = Dict()
    for (index, arg) in enumerate(ast.args[1].args[2:end])
        map[arg.args[1]] = args[index]
    end
    map
end

"""
Generates a loopnest with `body` as the body of the inner most loop using
`vars` as a list of loop variables and `ranges` as a list of ranges for each
loop nest.

Example:
    julia> gen_loop_nest(:(println((a, b))), [:a, :b], [1:2, 1:3])
    :(for b = 1:3
          for a = 1:2
              println((a, b))
          end
      end)
"""
function gen_loop_nest(body, vars, ranges)
    nest = body
    for (var, range) in zip(vars, ranges)
        nest = :(for $var = $range
            $nest
        end)
    end
    nest
end

"""
Synthesizes a loop nest around body based on the dimensionality of `buffer`.
We split each statement in `statements` into a separate synthesized loop nest
to facilitate pattern matching of statements.  If no pattern matching occurs,
the identical loop nests will be fused later.
"""
function append_neuron_loop!(body::Vector, statements::Vector, buffer::Array)
    N = ndims(buffer)
    for statement in statements
        vars   = [symbol(:_neuron_index_,i) for i in 1:N]
        ranges = [:(1:$(size(buffer,i)))    for i in 1:N]
        push!(body, gen_loop_nest(statement, vars, ranges))
    end
end


text:
#=
# 6. Variable free stream
In this notebook we will simulate the flow with a time-varying free stream past a
stationary body. To demonstrate this, we will solve for oscillatory flow past a
rectangular object, in which the $x$ component of the free stream is

$$U_\infty(t) = A \sin(\Omega t)$$
=#
using ViscousFlow
#-
using Plots

#=
### Problem specification
We will set the Reynolds number to 200
=#
my_params = Dict()
my_params["Re"] = 200

#=
In order to set a time-varying free stream, we have to define a function
that provides the instantaneous free stream components and pass that
function into the system definition. In this function, we will
use the `Sinusoid` function (available via the `RigidBodyTools` module)
to create the modulated free stream. To demonstrate its possibilities,
we will pass in the parameters for the sinusoid via the `my_params` dictionary.
The "freestream average" specifies a mean free stream, if desired.
=#
my_params["freestream average"] = 0.0
my_params["freestream frequency"]  = 2.0
my_params["freestream amplitude"] = 1.0
my_params["freestream phase"] = 0.0

#=
Now we define the function. We can call it anything we want,
but it has to have the argument signature as shown. The
`Sinusoid` function is used, with the shift operator `>>`
to apply any phase lag.
=#
function my_freestream(t,phys_params)
    U = phys_params["freestream average"]
    Ω = phys_params["freestream frequency"]
    Ax = phys_params["freestream amplitude"]
    ϕx = phys_params["freestream phase"]
    Vinfmag = Ax*(RigidBodyTools.Sinusoid(Ω) >> (ϕx/Ω))
    Vinf_angle = get(phys_params,"freestream angle",0.0)

    Uinf = (U + Vinfmag(t))*cos(Vinf_angle)
    Vinf = (U + Vinfmag(t))*sin(Vinf_angle)
    return Uinf, Vinf
end

#=
The freestream function is generically considered a forcing function,
so we pass it in via the "freestream" key in the forcing dictionary.
=#
forcing_dict = Dict("freestream" => my_freestream)

# Now let us carry on with the other usual steps:

# ### Discretize
xlim = (-2.0,2.0)
ylim = (-1.5,1.5)
my_params["grid Re"] = 4.0
g = setup_grid(xlim,ylim,my_params)

#=
### Set up bodies
Here, we will set up an ellipse in the center of the domain
=#
Δs = surface_point_spacing(g,my_params)
body = Ellipse(0.25,0.5,Δs)

T = RigidTransform((0.0,0.0),0.0)
T(body)
#-
plot(body,xlim=xlim,ylim=ylim)

#=
### Construct the system structure
This step is like the previous notebook, but now we also provide the body and the freestream:
=#
sys = viscousflow_system(g,body,phys_params=my_params,forcing=forcing_dict);
#=
### Initialize
Now, we initialize with zero vorticity
=#
u0 = init_sol(sys)
# and create the integrator
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

#=
### Solve
Now we are ready to solve the problem. Let's advance the solution to $t = 2.5$.
=#
@time step!(integrator,2.5)

#=
### Examine
Let's look at the flow field at the end of this interval
=#
sol = integrator.sol
plt = plot(layout = (2,2), size = (800, 600), legend=:false)
tsnap = 1.0:0.5:2.5
for (i, t) in enumerate(tsnap)
    plot!(plt[i],vorticity(sol,sys,t),sys,layers=false,title="t = $(round(t,digits=2))",clim=(-10,10),levels=range(-10,10,length=30),color = :RdBu)
    plot!(plt[i],surfaces(sol,sys,t))
end
plt

#=
#### Compute the force history
Just as we did for the stationary body in a constant free stream
=#
fx, fy = force(sol,sys,1);
# Plot them
plot(
plot(sol.t,2*fx,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_D\$",legend=:false),
plot(sol.t,2*fy,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=:false),
    size=(800,350)
)

# The mean drag and lift coefficients are
meanCD = GridUtilities.mean(2*fx[3:end])
#-
meanCL = GridUtilities.mean(2*fy[3:end])


text:
#!/usr/bin/env julia

using Luxor, Colors, Test, Random

Random.seed!(42)

function spirals()
    gsave()
    scale(0.3)
    r = 200
    setcolor("gray")
    for i = 0:pi/8:2pi
        gsave()
        translate(r * cos(i), r * sin(i))
        rotate(i)
        julialogo()
        grestore()
    end
    grestore()
    juliacircles(20)
end

function expandingspirals()
    gsave()
    scale(0.3)
    r = 200
    for i = pi:pi/12:6pi
        gsave()
        translate(i / 3 * r * cos(i), i / 3 * r * sin(i))
        scale(0.8)
        rotate(i)
        julialogo()
        grestore()
    end
    grestore()
end

function dropshadow()
    steps = 20
    # white-gray ramp
    gramp = range(
        colorant"white",
        stop = colorant"gray60",
        length = steps,
    )
    gsave()
    r = 200
    sethue("purple")
    rect(O, 5, 10, :stroke)
    setopacity(0.1)
    for i = 1:steps
        translate(-0.6, -0.5)
        julialogo(bodycolor = gramp[i], color = false)
    end
    julialogo()
    grestore()
end

function colorgrid()
    #cols = colormap("RdBu", 5; mid=0.5, logscale=false)
    #cols = sequential_palette(rand(10:360), 5, b=0.1)
    cols = distinguishable_colors(25)
    gsave()

    pagetiles = Tiler(500, 400, 5, 5)
    for (pos, n) in pagetiles
        gsave()
        setcolor(color(cols[n]))
        translate(pos)
        scale(0.3)
        julialogo(color = false)
        grestore()
    end
    grestore()
end

function boxes_and_rectangles(pt::Point)
    for i = 1:20
        randomhue()
        poly(
            box(pt + i, 10i, 10i, vertices = true),
            :stroke,
            close = true,
        )
    end
    for i = 40:10:150
        randomhue()
        translate(i, 0)
        poly(squircle(O, 30, 30, vertices = true), :stroke)
    end
end

function draw_julia_logos(fname)
    Drawing(1600, 1600, fname)
    origin()
    background("white")

    translate(-500, -200)
    spirals()

    translate(750, 0)
    expandingspirals()

    translate(-1000, 500)
    dropshadow()

    translate(800, 50)
    colorgrid()

    translate(-700, 300)
    boxes_and_rectangles(O)

    @test finish() == true
    println("...finished test: output in $(fname)")
end

# requires Luxor >= 2.0
function julialogodims(fname)
    Drawing(1600, 1600, fname)
    origin()
    background("white")
    sethue("black")
    julialogo(action = :path, centered = true)
    p = pathtopoly()
    lox, hix =
        extrema(first.(collect(Iterators.flatten(p))))
    loy, hiy = extrema(getfield.(
        collect(Iterators.flatten(p)),
        :y,
    ))
    bbx = BoundingBox(Point(lox, loy), Point(hix, hiy))
    sethue("grey58")
    box(bbx, :stroke)
    rule(boxmiddlecenter(bbx), boundingbox = bbx)
    rule(boxmiddlecenter(bbx), π / 2, boundingbox = bbx)
    dimension(
        boxtopleft(bbx),
        format = (d) -> string(round(d)),
        boxtopright(bbx),
        offset = -50,
        textgap = 0,
        textrotation = -π / 2,
        texthorizontaloffset = 20,
    )
    dimension(
        boxtopleft(bbx),
        format = (d) -> string(round(d)),
        boxbottomleft(bbx),
        textrotation = π,
        offset = 30,
        texthorizontaloffset = 20,
    )
    julialogo(centered = true)
    @test finish() == true
    println("...finished test: output in $(fname)")
end

function test1()
    randomhue()
    @layer begin
        translate(-100, 0)
        scale(0.5)
        julialogo(action = :path, centered = true)
        strokepath()
    end

    @layer begin
        translate(100, 0)
        scale(0.5)
        julialogo(action = :path, centered = true)
        strokepath()
    end
    text(
        "action=:path, centered=true)",
        halign = :center,
        Point(0, 100),
    )
end

function test2()
    randomhue()
    @layer begin
        translate(-100, 0)
        scale(0.5)
        julialogo(action = :clip, centered = true)
        paint()
        clipreset()
    end

    @layer begin
        translate(100, 0)
        scale(0.5)
        julialogo(action = :clip, centered = true)
        paint()
        clipreset()
    end

    text(
        "action=:clip, centered=true",
        halign = :center,
        Point(0, 100),
    )
end

function test3()
    sethue("blue")
    @layer begin
        translate(-100, 0)
        scale(0.5)
        julialogo(action = :stroke, centered = true)
    end
    @layer begin
        translate(100, 0)
        scale(0.5)
        julialogo(action = :stroke, centered = true)
    end
    text(
        "OG> action=:stroke, centered=true <NEW",
        halign = :center,
        Point(0, 100),
    )
end

function test4()
    randomhue()
    @layer begin
        translate(-100, 0)
        scale(0.5)
        julialogo(
            centered = true,
            color = true,
            bodycolor = "purple",
        )
    end

    @layer begin
        translate(100, 0)
        scale(0.5)
        julialogo(
            centered = true,
            color = true,
            bodycolor = "purple",
        )
    end

    text(
        "OG> centered=true, bodycolor=\"purple\" <NEW",
        halign = :center,
        Point(0, 100),
    )
end

function test5()
    randomhue()
    @layer begin
        translate(-100, 0)
        scale(0.5)
        julialogo(
            centered = true,
            color = true,
            bodycolor = "purple",
        )
    end

    @layer begin
        translate(100, 0)
        scale(0.5)
        julialogo(
            centered = true,
            color = true,
            bodycolor = "purple",
        )
    end

    text(
        "centered=true, color=true, bodycolor=\"purple\"",
        halign = :center,
        Point(0, 100),
    )
end

function test6()
    randomhue()
    @layer begin
        translate(-100, 0)
        scale(0.5)
        julialogo(
            centered = true,
            action = :path,
            color = true,
        )
    end
    fillpath()
    @layer begin
        translate(100, 0)
        scale(0.5)
        julialogo(
            centered = true,
            action = :path,
            color = true,
        )
    end
    fillpath()

    text(
        "centered=true, color=true, bodycolor=\"purple\"",
        halign = :center,
        Point(0, 100),
    )
end

function draw_julia_logos(fname)
    Drawing(1600, 1600, fname)
    origin()
    background("white")
    fontsize(12)
    tiles = Tiler(1200, 1200, 3, 2)
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            n == 1 && test1()
            n == 2 && test2()
            n == 3 && test3()
            n == 4 && test4()
            n == 5 && test5()
            n == 6 && test6()
        end
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
    finish()
end

draw_julia_logos("julia-logo-draw.png")
julialogodims("julia-logo-dims.png")
draw_julia_logos("julia-logo-draw-table.png")


text:
using LinearAlgebra, Statistics;
function Newtonfixedpointmap(f, f_prime, x_0, tolerance, maxiter)
    # setup the algorithm
    x_old = x_0
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - (f(x_old)/ f_prime(x_old)) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    return (x_old, normdiff, iter)
end

maxiter = 1000
tolerance = 1.0E-7
x_0 = 0.0 # initial condition

f(v) = (v-1)^3
f_prime(v)= 3*(v-1)^2

x_star, normdiff, iter = Newtonfixedpointmap(f, f_prime, x_0, tolerance, maxiter)
println("Fixed point = $x_star, and |f(x) - x| = $normdiff in $iter iterations")

f(v)= 27*v^3 - 3*v + 1
f_prime(v)= 81*v^2 - 3

x_star, normdiff, iter = Newtonfixedpointmap(f, f_prime, x_0, tolerance, maxiter)
println("Fixed point = $x_star, and |f(x) - x| = $normdiff in $iter iterations")

using ForwardDiff

# operator to get the derivative of this function using AD
D(f) = v-> ForwardDiff.derivative(f, v)

function Newtonfixedpointmap(f, x_0, tolerance, maxiter)
    # setup the algorithm
    x_old = x_0
    normdiff = Inf
    iter = 1
    f_prime = D(f)
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - (f(x_old)/ f_prime(x_old)) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    return (x_old, normdiff, iter)
end
maxiter = 1000
tolerance = 1.0E-7
x_0 = 0.0 # initial condition

f(v) = (v-1)^3
f_prime(v)= 3*(v-1)^2

x_star, normdiff, iter = Newtonfixedpointmap(f, f_prime, x_0, tolerance, maxiter)
println("Fixed point = $x_star, and |f(x) - x| = $normdiff in $iter iterations")

f(v)= 27*v^3 - 3*v + 1
f_prime(v)= 81*v^2 - 3

x_star, normdiff, iter = Newtonfixedpointmap(f, f_prime, x_0, tolerance, maxiter)
println("Fixed point = $x_star, and |f(x) - x| = $normdiff in $iter iterations")


text:
abstract type AbstractLogger end
abstract type AbstractTermination end

function terminate(::AbstractTermination,k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector) end

function termination_wrapper(k::Cint,fval::Cdouble,xbegin::Ptr{Cdouble},
                             xend::Ptr{Cdouble}, gbegin::Ptr{Cdouble},
                             termination_data::Ptr{Nothing})
    termination = unsafe_pointer_to_objref(termination_data)::AbstractTermination
    ptrdiff = Int(xend - xbegin)
    N = divrem(ptrdiff, sizeof(Cdouble))[1]
    x = unsafe_wrap(Array, xbegin, N)
    g = unsafe_wrap(Array, gbegin, N)
    return Cint(terminate(termination,k,fval,x,g))
end

function log(::AbstractLogger,k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector) end

function log_wrapper(k::Cint,fval::Cdouble,xbegin::Ptr{Cdouble},
                     xend::Ptr{Cdouble}, gbegin::Ptr{Cdouble},
                     logger_data::Ptr{Nothing})
    logger = unsafe_pointer_to_objref(logger_data)::AbstractLogger
    ptrdiff = Int(xend - xbegin)
    N = divrem(ptrdiff, sizeof(Cdouble))[1]
    x = unsafe_wrap(Array, xbegin, N)
    g = unsafe_wrap(Array, gbegin, N)
    log(logger,k,fval,x,g)
    return nothing
end

module Utility

using POLO
using POLO: AbstractTermination, AbstractLogger
using LinearAlgebra
using SparseArrays

import POLO: terminate

include("termination.jl")
include("logging.jl")

# TODO:
# 1. Maybe we can redesign `reader`s as above, i.e., `AbstractReader`.
# 2. Implement `svmreader` to support "dense" reading --- that will be more
#    efficient compared to the `push!`s.
# 3. One can still make improvements when `m` and `n` are known.
# 4. Add tests.
include("svmreader.jl")

end


text:
using JLD2
using GLMakie
using Statistics

tracer_case_j = 1
tracer_case_k = 1
# file output
# prefix = "fluxernathy_tracers_iteration9986400.jld2"
prefix = "fluxernathy_tracers_restarted_j1_k1_averages.jld2"
# prefix = "checkme.jld2"
# prefix = "fluxernathy_nh_ns_averages.jld2"
# prefix = "fluxernathy_averages.jld2"
# prefix = "fluxernathy_tracers_restarted_j1_k10_averages.jld2"
# prefix = "fluxernathy_tracers_restarted_j1_k20_averages.jld2"
# prefix = "check_hhh.jld2"
# prefix = "fluxernathy_tracers_restarted_j128_k1_averages.jld2"
prefix = "check_h200.jld2"
jl_file = jldopen(prefix, "r+")

function get_grid(field, jl_file; ghost = 3)

    yC = collect(jl_file["grid"]["yC"])[ghost+1:end-ghost]
    zC = collect(jl_file["grid"]["zC"])[ghost+1:end-ghost]

    yF = collect(jl_file["grid"]["yF"])[ghost+1:end-ghost]
    zF = collect(jl_file["grid"]["zF"])[ghost+1:end-ghost]

    y = size(field)[1] == size(yC)[1] ? yC : yF
    z = size(field)[2] == size(zC)[1] ? zC : zF

    return y, z
end

function get_grid2(field, jl_file; ghost = 3)

    yC = collect(jl_file["grid"]["yᵃᶜᵃ"])[ghost+1:end-ghost]
    zC = collect(jl_file["grid"]["zᵃᵃᶜ"])[ghost+1:end-ghost]

    yF = collect(jl_file["grid"]["yᵃᶠᵃ"])[ghost+1:end-ghost]
    zF = collect(jl_file["grid"]["zᵃᵃᶠ"])[ghost+1:end-ghost]

    y = size(field)[1] == size(yC)[1] ? yC : yF
    z = size(field)[2] == size(zC)[1] ? zC : zF

    return y, z
end

function get_field(field_string::String, time_index::Int, jl_file::JLD2.JLDFile; ghost = 3)
    t_keys = keys(jl_file["timeseries"][field_string])
    field = jl_file["timeseries"][field_string][t_keys[time_index]]
    return field[1,:,:]
end

function ∂y(field, y)
    δf = field[2:end, :] - field[1:end-1, :]
    δy = y[2:end] - y[1:end-1]
    δy = reshape(δy[:], (length(δy), 1))
    return δf ./ δy
end

function ∂z(field,z)
    δf = field[:, 2:end] - field[:, 1:end-1]
    δz = z[2:end] - z[1:end-1]
    δz = reshape(δz[:], (1, length(δz)))
    return δf ./ δz
end

function avg(field)
    return (field[2:end] + field[1:end-1]) .* 0.5
end
function avgy(field)
    return (field[2:end, :] + field[1:end-1, :]) .* 0.5
end
function avgz(field)
    return (field[:, 2:end] + field[:, 1:end-1]) .* 0.5
end

#=
field = get_field("c1", 9, jl_file)
y, z = get_grid(field, jl_file)

y̅ = avg(y)
z̅ = avg(z)
fieldz = ∂z(field,z)
# heatmap(y, z, field, colormap = :balance, interpolate = true)

heatmap(y, z̅, fieldz, colormap = :balance, interpolate = true)

field1 = get_field("c3", 9, jl_file)
field2 = get_field("c4", 9, jl_file)
field1_m = ∂z(field1, z)
field2_m = ∂z(field2, z)

fig, ax, hm = heatmap(y, z̅, field1_m, colormap = :balance, interpolate = true)
ax2 = fig[1,2] = Axis(fig)
heatmap!(ax2, y, z̅, field2_m, colormap = :balance, interpolate = true)

ax3 = fig[1,3] = Axis(fig)
cc = maximum(abs.(field1_m))
clims = (-cc, cc)
heatmap!(ax3, y, z̅, field1_m - field2_m, colormap = :balance, interpolate = true, colorrange = clims)
=#

# data pairs
#=
custom_clims = false
tracer_string = "c1"
time_index = 20
b = get_field("b", time_index, jl_file)
c = get_field(tracer_string, time_index, jl_file)
v = get_field("v", time_index, jl_file)
w = get_field("w", time_index, jl_file)
vc = get_field("v"*tracer_string, time_index, jl_file)
wc = get_field("w"*tracer_string, time_index, jl_file)
vcp = avgy(vc) - avgy(v) .* c
wcp = avgz(wc) - avgz(w) .* c

uq = 0.95 #upper quantile
lq = 0.05 #lower quantile

y,z = get_grid(c, jl_file; ghost = 3)

y̅ = avg(y)
z̅ = avg(z)

cz = ∂z(c, z)
cy = ∂y(c, y)

bz = ∂z(b, z)
by = ∂y(b, y)

if custom_clims
    clims = (0, 2e-5)
else
    clims = quantile.(Ref(cz[:]), [lq, uq])
end

fig, ax, hm = heatmap(y, z̅, cz, colormap = :balance, interpolate = true, colorrange = clims)
ax.title = "∂z(" * tracer_string * ")"
Colorbar(fig[1,2], hm)

if custom_clims
    clims2 = (-2e-9, 1e-8)
else
    clims2 = quantile.(Ref(cy[:]), [lq, uq])
end

ax2 = fig[1,3] = Axis(fig)
hm2 = heatmap!(ax2, y̅, z, cy, colormap = :balance, interpolate = true, colorrange = clims2)
ax2.title = "∂y(" * tracer_string * ")"
Colorbar(fig[1,4], hm2)

if custom_clims
    clims3 = (-1e-9, 2e-8)
else
    clims3 = quantile.(Ref(wcp[:]), [lq,uq])
end

ax3 = fig[2,1] = Axis(fig)
hm3 = heatmap!(ax3, y, z, wcp, colormap = :balance, interpolate = true, colorrange = clims3)
ax3.title = "w'"* tracer_string * "'"
Colorbar(fig[2,2], hm3)

if custom_clims
    clims4 = (-2e-5, 0)
else
    clims4 = quantile.(Ref(vcp[:]), [lq,uq])
end

ax4 = fig[2,3] = Axis(fig)
hm4 = heatmap!(ax4, y, z, vcp, colormap = :balance, interpolate = true, colorrange = clims4)
ax4.title = "v'"* tracer_string * "'"
Colorbar(fig[2,4], hm4)

##
vcp_avg = avgz(avgy(vcp))
wcp_avg = avgz(avgy(wcp))
bz_avg = avgy(bz)
by_avg = avgz(by)

dcp = vcp_avg .* by_avg + wcp_avg .* bz_avg
icp = vcp_avg .* bz_avg - wcp_avg .* by_avg

if custom_clims
    clims5 = quantile.(Ref(dcp[:]), [lq,uq])
else
    clims5 = quantile.(Ref(dcp[:]), [lq,uq])
end

if custom_clims
    clims6 = quantile.(Ref(icp[:]), [lq,uq])
    clims5 = clims6
else
    clims6 = quantile.(Ref(icp[:]), [lq,uq])
    clims5 = clims6
end

ax5 = fig[3,1] = Axis(fig)
hm5 = heatmap!(ax5, y̅, z̅, dcp , colormap = :balance, interpolate = true, colorrange = clims5)
ax5.title = "u⃗'"* tracer_string * "'" * "⋅∇b"
Colorbar(fig[3,2], hm5)

ax6 = fig[3,3] = Axis(fig)
hm6 = heatmap!(ax6, y̅, z̅, icp, colormap = :balance, interpolate = true, colorrange = clims6)
ax6.title = "u⃗'"* tracer_string * "'" * "⋅∇ᵖb"
Colorbar(fig[3,4], hm6)
=#

clim_options = [:custom, :symmetric]
eke_avail = true
diff_field = false
buoyancy_contour = true
check_grad = true
custom_clims =  clim_options[2]
tracer_string = "c4"
tracer_string2 = nothing # "c2" # "c3" #nothing # "c2" # nothing # "c2" # nothing
time_index = 14
b = get_field("b", time_index, jl_file)
if tracer_string2==nothing
    c = get_field(tracer_string, time_index, jl_file)
else
    c = get_field(tracer_string, time_index, jl_file) - get_field(tracer_string2, time_index, jl_file)
end
v = get_field("v", time_index, jl_file)
w = get_field("w", time_index, jl_file)
if tracer_string2 == nothing
    vc = get_field("v"*tracer_string, time_index, jl_file)
    wc = get_field("w"*tracer_string, time_index, jl_file)
else
    vc = get_field("v" * tracer_string, time_index, jl_file) - get_field("v"*tracer_string2, time_index, jl_file)
    wc = get_field("w" * tracer_string, time_index, jl_file) - get_field("w"*tracer_string2, time_index, jl_file)
end

vcp = avgy(vc) - avgy(v) .* c
wcp = avgz(wc) - avgz(w) .* c

uq = 0.95 #upper quantile
lq = 0.05 #lower quantile

y,z = get_grid(c, jl_file; ghost = 3)

y̅ = avg(y)
z̅ = avg(z)

cz = ∂z(c, z)
cy = ∂y(c, y)

bz = ∂z(b, z)
by = ∂y(b, y)

vcp_avg = avgz(avgy(vcp))
wcp_avg = avgz(avgy(wcp))
bz_avg = avgy(bz)
by_avg = avgz(by)

fig = Figure(resolution = (2800, 1400))

clims01 = quantile.(Ref(c[:]), [lq, uq])

ax01 = fig[1,1]= Axis(fig)
hm01 = heatmap!(ax01, y, z, c, colormap = :thermometer, interpolate = true, colorrange = clims01)
ax01.title = tracer_string
Colorbar(fig[1,2], hm01)

∇b = sqrt.(bz_avg .^2 + by_avg .^2)
clims03 = quantile.(Ref(∇b[:]), [lq, uq])
ax03 = fig[3,1]= Axis(fig)
hm03 = heatmap!(ax03, y, z, ∇b, colormap = :thermometer, interpolate = true, colorrange = clims03)
ax03.title = "|∇b|"
Colorbar(fig[3,2], hm03)

if eke_avail
    u = get_field("u", time_index, jl_file)
    v = get_field("v", time_index, jl_file)
    w = get_field("w", time_index, jl_file)
    uu = get_field("uu", time_index, jl_file)
    vv = get_field("vv", time_index, jl_file)
    ww = get_field("ww", time_index, jl_file)

    eke = uu .- (u .* u) + avgy(vv .- (v .* v) )  + avgz(ww .- (w .* w) )

    clims02 = quantile.(Ref(eke[:]), [lq, uq])
    ax02 = fig[2,1]= Axis(fig)
    hm02 = heatmap!(ax02, y, z, eke, colormap = :thermometer, interpolate = true, colorrange = clims02)
    ax02.title = "eke"
    Colorbar(fig[2,2], hm03)
else
    nothing
end

if custom_clims==:custom
    clims = (0, 2e-5)
elseif custom_clims==:symmetric
    clims = quantile.(Ref(cz[:]), [lq, uq])
    cu = maximum(abs.(clims))
    clims = (-cu, cu)
else
    clims = quantile.(Ref(cz[:]), [lq, uq])
end

ax = fig[1,1+2]= Axis(fig)
hm = heatmap!(ax, y, z̅, cz, colormap = :balance, interpolate = true, colorrange = clims)
ax.title = "∂z(" * tracer_string * ")"
Colorbar(fig[1,2+2], hm)

if custom_clims==:custom
    clims2 = (-2e-9, 1e-8)
elseif custom_clims==:symmetric
    clims2 = quantile.(Ref(cy[:]), [lq, uq])
    cu = maximum(abs.(clims2))
    clims2 = (-cu, cu)
else
    clims2 = quantile.(Ref(cy[:]), [lq, uq])
end

ax2 = fig[1,3+2] = Axis(fig)
hm2 = heatmap!(ax2, y̅, z, cy, colormap = :balance, interpolate = true, colorrange = clims2)
ax2.title = "∂y(" * tracer_string * ")"
Colorbar(fig[1,4+2], hm2)

if custom_clims==:custom
    clims3 = (-1e-9, 2e-8)
elseif custom_clims==:symmetric
    clims3 = quantile.(Ref(wcp[:]), [lq, uq])
    cu = maximum(abs.(clims3))
    clims3 = (-cu, cu)
else
    clims3 = quantile.(Ref(wcp[:]), [lq,uq])
end

ax3 = fig[2,1+2] = Axis(fig)
hm3 = heatmap!(ax3, y, z, wcp, colormap = :balance, interpolate = true, colorrange = clims3)
ax3.title = "w'"* tracer_string * "'"
Colorbar(fig[2,2+2], hm3)

if custom_clims==:custom
    clims4 = (-2e-5, 0)
elseif custom_clims==:symmetric
    clims4 = quantile.(Ref(vcp[:]), [lq, uq])
    cu = maximum(abs.(clims4))
    clims4 = (-cu, cu)
else
    clims4 = quantile.(Ref(vcp[:]), [lq,uq])
end

ax4 = fig[2,3+2] = Axis(fig)
hm4 = heatmap!(ax4, y, z, vcp, colormap = :balance, interpolate = true, colorrange = clims4)
ax4.title = "v'"* tracer_string * "'"
Colorbar(fig[2,4+2], hm4)

##
vcp_avg = avgz(avgy(vcp))
wcp_avg = avgz(avgy(wcp))
bz_avg = avgy(bz)
by_avg = avgz(by)

dcp = vcp_avg .* by_avg + wcp_avg .* bz_avg
icp = vcp_avg .* bz_avg - wcp_avg .* by_avg

if custom_clims==:custom
    clims5 = quantile.(Ref(dcp[:]), [lq,uq])
elseif custom_clims==:symmetric
    clims5 = quantile.(Ref(dcp[:]), [lq, uq])
    cu = maximum(abs.(clims5))
    clims5 = (-cu, cu)
else
    clims5 = quantile.(Ref(dcp[:]), [lq,uq])
end

if custom_clims==:custom
    clims6 = quantile.(Ref(icp[:]), [lq,uq])
    clims5 = clims6
elseif custom_clims==:symmetric
    clims6 = quantile.(Ref(icp[:]), [lq, uq])
    cu = maximum(abs.(clims6))
    clims6 = (-cu, cu)
    # clims5 = clims6
else
    clims6 = quantile.(Ref(icp[:]), [lq,uq])
    clims5 = clims6
end

ax5 = fig[3,1+2] = Axis(fig)
hm5 = heatmap!(ax5, y̅, z̅, dcp , colormap = :balance, interpolate = true, colorrange = clims5)
ax5.title = "u⃗'"* tracer_string * "'" * "⋅∇b"
Colorbar(fig[3,2+2], hm5)

ax6 = fig[3,3+2] = Axis(fig)
hm6 = heatmap!(ax6, y̅, z̅, icp, colormap = :balance, interpolate = true, colorrange = clims6)
ax6.title = "u⃗'"* tracer_string * "'" * "⋅∇ᵖb"
Colorbar(fig[3,4+2], hm6)

if buoyancy_contour
    blevels = b[end, 1:2:end]
    for axval in [ax01, ax02, ax03, ax, ax2, ax3, ax4, ax5, ax6]
        contour!(axval, y,z, b, levels = blevels, color = :black, linewidth = 1)
    end
end

if check_grad

    ∇c∇b  = avgz(cy) .* by_avg  + avgy(cz) .* bz_avg
    ∇c∇ᵖb = avgz(cy) .* bz_avg  - avgy(cz) .* by_avg

    if custom_clims==:custom
        clims7 = quantile.(Ref(∇c∇b[:]), [lq,uq])
    elseif custom_clims==:symmetric
        clims7 = quantile.(Ref(∇c∇b[:]), [lq, uq])
        cu = maximum(abs.(clims7))
        clims7 = (-cu, cu)
    else
        clims7 = quantile.(Ref(∇c∇b[:]), [lq,uq])
    end

    if custom_clims==:custom
        clims8 = quantile.(Ref(∇c∇ᵖb[:]), [lq,uq])
        clims7 = clims8
    elseif custom_clims==:symmetric
        clims8 = quantile.(Ref(∇c∇ᵖb[:]), [lq, uq])
        cu = maximum(abs.(clims8))
        clims8 = (-cu, cu)
        # clims5 = clims6
    else
        clims8 = quantile.(Ref(∇c∇ᵖb[:]), [lq,uq])
        clims7 = clims8
    end

    ax7 = fig[4,1+2+2] = Axis(fig)
    hm7 = heatmap!(ax7, y̅, z̅, -1 .* ∇c∇b , colormap = :balance, interpolate = true, colorrange = clims7)
    ax7.title = "-∇"* tracer_string * "⋅∇b = " * "-∇ᵖ"* tracer_string * "⋅∇ᵖb"
    Colorbar(fig[4,2+2+2], hm7)

    ax8 = fig[4,3+2-2] = Axis(fig)
    hm8 = heatmap!(ax8, y̅, z̅, ∇c∇ᵖb, colormap = :balance, interpolate = true, colorrange = clims8)
    ax8.title = "∇"* tracer_string * "⋅∇ᵖb = " * "-∇ᵖ"* tracer_string * "⋅∇b"
    Colorbar(fig[4,4+2-2], hm8)

    blevels = b[end, 1:2:end]

    # quick check

    ∇c = sqrt.(avgz(cy) .* avgz(cy)  + avgy(cz) .* avgy(cz))
    clims04 = quantile.(Ref( ∇c[:]), [lq, uq])
    ax04 = fig[4,1]= Axis(fig)
    hm04 = heatmap!(ax04, y̅, z̅, ∇c , colormap = :thermometer, interpolate = true, colorrange = clims04)
    ax04.title = "|∇"* tracer_string * "|"
    Colorbar(fig[4,2], hm04)

    for axval in [ax7, ax8, ax04]
        contour!(axval, y,z, b, levels = blevels, color = :black, linewidth = 1)
    end
end

for axval in [ax01, ax02, ax03, ax, ax2, ax3, ax4, ax5, ax6,      ax7, ax8, ax04]
    axval.titlesize = 40
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
"""
Google Cloud Identity Access Management (IAM) API
"""
module _iam

export iam

using ..api
using ...root

"""
Google Cloud IAM API root.
"""
iam = APIRoot(
    "https://iam.googleapis.com/v1",
    Dict(
        "cloud-platform" => "Full access to all resources and services in the specified Cloud Platform project.",
    );
    ServiceAccount=APIResource("projects/{project}/serviceAccounts";
        create=APIMethod(:POST, "", "Creates a ServiceAccount and returns it."),
        delete=APIMethod(:DELETE, "{serviceAccount}", "Deletes a ServiceAccount."),
        get=APIMethod(:GET, "{serviceAccount}", "Gets a ServiceAccount."),
        getIamPolicy=APIMethod(:POST, "{serviceAccount}:getIamPolicy", "Returns the IAM access control policy for a ServiceAccount."),
        list=APIMethod(:GET, "", "Lists ServiceAccounts for a project.";
            transform=(x, t) -> map(t, get(x, :accounts, []))
        ),
        setIamPolicy=APIMethod(:POST, "{serviceAccount}:setIamPolicy", "Sets the IAM access control policy for a ServiceAccount."),
        signBlob=APIMethod(:POST, "{serviceAccount}:signBlob", "Signs a blob using a service account's system-managed private key."),
        testIamPermissions=APIMethod(:POST, "{serviceAccount}:testIamPermissions", "Tests the specified permissions against the IAM access control policy for a ServiceAccount."),
        update=APIMethod(:PUT, "{serviceAccount}", "Updates a ServiceAccount."),
    ),
    Key=APIResource("projects/{project}/serviceAccounts/{serviceAccount}/keys";
        create=APIMethod(:POST, "", "Creates a ServiceAccountKey and returns it."),
        delete=APIMethod(:DELETE, "{key}", "Deletes a ServiceAccountKey."),
        get=APIMethod(:GET, "{key}", "Gets the ServiceAccountKey by key id."),
        list=APIMethod(:GET, "", "Lists ServiceAccountKeys.";
            transform=(x, t) -> map(t, get(x, :keys, []))
        ),
    ),
    Role=APIResource("roles";
        queryGrantableRoles=APIMethod(:POST, ":queryGrantableRoles", "Queries roles that can be granted on a particular resource."),
    ),
)

end


text:
@info "Visualization: [PoissonMatrix (5x5) (λ=5)]"
display(PoissonMatrix(5,5,rank=3,λ=5))
@info "Visualization: [GaussianMatrix (5x5) (μ=0,σ=1)]"
display(GaussianMatrix(5,5,rank=3,σ=1,μ=0))
@info "Visualization: [BinomialMatrix (5x5) (p=0.5)]"
display(BernoulliMatrix(5,5,rank=3,p=0.5))
@info "Visualization: [GammaMatrix (5,5) (α=5,θ=0.5)]"
display(GammaMatrix(5,5,rank=3,α=5,θ=0.5))