text:
function short_year(year)
    short_year_str = string(year)[end-1:end]
    (short_year_str, int(short_year_str))
end

# -------

function pacs_to_candidates(year)
    short_year_str, short_year_int = short_year(year)
    data_path = joinpath(get(_campaign_finance_data_sources, short_year_int, ""), "pacs$(short_year_str).txt")
    data_stream = OpenSecretsBuffer(Latin1Buffer(open(data_path, "r")))

    readtable(data_stream, separator = ',', header = false, quotemark = ['|'],
              names = [:Cycle, :FECRecNo, :PACID, :CID, :Amount, :Date, :RealCode, :Type, :DI, :FECCandID],
              eltypes = [Int, UTF8String, UTF8String, UTF8String, Float64, UTF8String, UTF8String, UTF8String,
              UTF8String, UTF8String])
end

# -------

function committees(year)
    short_year_str, short_year_int = short_year(year)
    data_path = joinpath(get(_campaign_finance_data_sources, short_year_int, ""), "cmtes$(short_year_str).txt")
    data_stream = OpenSecretsBuffer(Latin1Buffer(open(data_path, "r")))

    readtable(data_stream; separator = ',', header = false, quotemark = ['|'], truestrings = ["1"], falsestrings = ["0"],
              names = [:Cycle, :CmteID, :PACShort, :Affiliate, :Ultorg, :RecipID, :RecipCode,
              :FECCandID, :Party, :PrimCode, :Source, :Sensitive, :Foreign, :Active],
              eltypes = [Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
              UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, Int, Bool])
end

# -------

function candidates(year)
    short_year_str, short_year_int = short_year(year)
    data_path = joinpath(get(_campaign_finance_data_sources, short_year_int, ""), "cands$(short_year_str).txt")
    data_stream = OpenSecretsBuffer(Latin1Buffer(open(data_path, "r")))

    readtable(data_stream; separator = ',', header = false, quotemark = ['|'],
              names = [:Cycle, :FECCandID, :CID, :FirstLastP, :Party, :DistIDRunFor, :DistIDCurr, :CurrCand,
              :CycleCand, :CRPICO, :RecipCode, :NoPacs],
              eltypes = [Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
              UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String])
end

# -------

function pacs_to_other(year)
    short_year_str, short_year_int = short_year(year)
    data_path = joinpath(get(_campaign_finance_data_sources, short_year_int, ""), "pac_other$(short_year_str).txt")
    data_stream = OpenSecretsBuffer(Latin1Buffer(open(data_path, "r")))

    readtable(data_stream; separator = ',', header = false, quotemark = ['|'],
              names = [:Cycle, :FECRecNo, :Filerid, :DonorCmte, :ContribLendTrans, :City, :State, :Zip, :FECOccEmp,
              :Primcode, :Date, :Amount, :RecipID, :Party, :Otherid, :RecipCode, :RecipPrimcode, :Amend, :Report,
              :PG, :Microfilm, :Type, :RealCode, :Source],
              eltypes = [Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
              UTF8String, UTF8String, UTF8String, Float64, UTF8String, UTF8String, UTF8String, UTF8String,
              UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String])
end

# -------

function individual_contributions(year)
    short_year_str, short_year_int = short_year(year)
    data_path = joinpath(get(_campaign_finance_data_sources, short_year_int, ""), "indivs$(short_year_str).txt")
    data_stream = OpenSecretsBuffer(Latin1Buffer(open(data_path, "r")))

    # the OpenData naming scheme isn't really scalable, but for now, 90 and above for 1990 - 1999
    # and 00+ for 2000 and on.
    if 90 <= short_year_int || short_year_int <= 10
        field_names = [:Cycle, :FECTransID, :ContribID, :Contrib, :RecipID, :Orgname, :UltOrg, :RealCode, :Date,
                      :Amount, :Street, :City, :State, :Zip, :RecipCode, :Type, :CmteID, :OtherID, :Gender,
                      :FecOccEmp, :Microfilm, :Occ_EF, :Emp_EF, :Source]
        field_vals = [Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
                     UTF8String, Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
                     UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String]
    else
        field_names = [:Cycle, :FECTransID, :ContribID, :Contrib, :RecipID, :Orgname, :UltOrg, :RealCode, :Date,
                      :Amount, :Street, :City, :State, :Zip, :RecipCode, :Type, :CmteID, :OtherID, :Gender,
                      :Microfilm, :Occupation, :Employer, :Source]
        field_vals = [Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
                     UTF8String, Int, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String,
                     UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String, UTF8String]
    end

    readtable(data_stream; separator = ',', header = false, quotemark = ['|'],
              names = field_names, eltypes = field_vals)
end


text:
replacecoordinate!(x, i::Integer, val) = (x[i] = val; x)

replacecoordinate!(x::SVector{N,T}, i::Integer, val) where {N,T} =
    SVector{N,T}(_rpc(Tuple(x), i-1, T(val)))
@inline _rpc(t, i, val) = (ifelse(i == 0, val, t[1]), _rpc(tail(t), i-1, val)...)
_rps(::Tuple{}, i, val) = ()

ipcopy!(dest, src) = copy!(dest, src)
ipcopy!(dest::SVector, src) = src

"""
    ileft, iright, iclosest = splitindexes(box)

Among the evaluated split points in `box.parent`, which correspond to
the closest points to `box`. `ileft` is the index of the point to the "left,"
`iright` to the right, and `iclosest` is whichever of these actually lies on
one of the boundaries of `box`.
"""
function splitindexes(box)
    isroot(box) && error("root box doesn't have a parent")
    p = box.parent
    cindex = box.parent_cindex
    overhang_left = p.children[1].edge < p.xvalues[1]
    offset = (1-2*overhang_left)/4 + 1
    xn = (cindex - 1)/2 + offset
    floor(Int, xn), ceil(Int, xn), round(Int, xn)
end

"""
   xmin, fmin = qfit_min((xm, fm), (x0, f0), (xp, fp))

Given three points `xm < x0 < xp` and three corresponding values `fm`, `f0`, and `fp`,
fit a quadratic. If the quadratic is convex, return the position of the minimum and
the corresponding value of the quadratic.

If the quadratic is not convex, returns a point outside the interval by an amount equal to
the separation of the nearest pair (e.g., if f is increasing, returns a value `xm - (x0-xm)`).
"""
function qfit_min(xfm, xf0, xfp)
    xm, fm = xfm
    x0, f0 = xf0
    xp, fp = xfp
    @assert(xp > x0 && x0 > xm && isfinite(xm) && isfinite(xp))
    cm = fm/((xm-x0)*(xm-xp))  # coefficients of Lagrange polynomial
    c0 = f0/((x0-xm)*(x0-xp))
    cp = fp/((xp-xm)*(xp-x0))
    csum = cm+c0+cp
    qvalue(x) = cm*(x-x0)*(x-xp) + c0*(x-xm)*(x-xp) + cp*(x-xm)*(x-x0)
    if csum > 0
        # It's convex
        xmin = (cm*(x0+xp) + c0*(xm+xp) + cp*(xm+x0))/(2*csum)
        return xmin, qvalue(xmin)
    end
    if fm == f0 == fp
        return x0, f0
    end
    @assert(f0 >= fm || f0 >= fp)
    # This should only happen at a box edge. Just return something
    # to indicate the direction of the minimum
    if f0 > fm && fp > fm
        xmin = xm - (x0 - xm)
    else
        xmin = xp + (xp - x0)
    end
    return xmin, qvalue(xmin)
end

function Base.minimum(box::Box)
    isroot(box) && error("cannot compute minimum for the root box")
    ileft, iright, iclosest = splitindexes(box)
    p = box.parent
    (ileft == 0 || iright > length(p.xvalues)) && return p.fvalues[iclosest]
    # Compute value on "internal" edge (between xvalues) by linear interpolation
    frac = abs(box.edge - p.xvalues[iclosest])/(p.xvalues[iright] - p.xvalues[ileft])
    iother = iclosest == ileft ? iright : ileft
    fc, fo = p.fvalues[iclosest], p.fvalues[iother]
    return min(fc, fc + (fo-fc)*frac)
end

## Tree traversal
function get_root(box::Box)
    while !isroot(box)
        box = box.parent
    end
    box
end

struct DepthFirstIterator{T}
    root::Box{T}
end
struct DepthFirstState{T}
    cparent::Box{T}
    cindex::Int
end

function visit_leaves(root::Box)
    DepthFirstIterator(root)
end

function Base.start(iter::DepthFirstIterator)
    find_next_leaf(iter, DepthFirstState(iter.root, 0))
end
function Base.done(iter::DepthFirstIterator, state::DepthFirstState)
    !issplit(iter.root) && return true
    iter.root == state.cparent && state.cindex > length(iter.root.children)
end
function Base.next(iter::DepthFirstIterator, state::DepthFirstState)
    box = state.cparent.children[state.cindex]
    @assert(!issplit(box))
    return (box, find_next_leaf(iter, state))
end
function find_next_leaf(iter::DepthFirstIterator, state::DepthFirstState)
    !issplit(iter.root) && return state
    box, i = state.cparent, state.cindex+1
    if i > length(box.children)
        box, i = up(box, iter.root)
    end
    if i <= length(box.children) && issplit(box.children[i])
        return find_next_leaf(iter, DepthFirstState(box.children[i], 0))
    end
    DepthFirstState(box, i)
end
function up(box, root)
    local i
    while true
        box, i = box.parent, box.parent_cindex+1
        box == root && return (box, i)
        i <= length(box.children) && break
    end
    return (box, i)
end

ext: jl


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

export
    av_adts_header_parse

function av_adts_header_parse(buf, samples, frames)
    ccall((:av_adts_header_parse, libavcodec), Cint, (Ptr{UInt8}, Ptr{UInt32}, Ptr{UInt8}), buf, samples, frames)
end


text:
#!/usr/bin/env julia
"""
A small Julia module that defines a simple function, to be used in UCBjulia.py (with pyjulia).
"""

# Small Julia module to wrap the function that computes a UCB index
module UCBjulia
    function index(rewards, pulls, t, arm, alpha=4)
        if pulls[arm] < 1
            return Inf
        else
            return (rewards[arm] / pulls[arm]) + sqrt((alpha * log(t)) / (2 * pulls[arm]))
        end
    end
end

# Small Julia function that computes a UCB index
function index(rewards, pulls, t, arm, alpha=4)
    if pulls[arm] < 1
        return Inf
    else
        return (rewards[arm] / pulls[arm]) + sqrt((alpha * log(t)) / (2 * pulls[arm]))
    end
end


text:
# function nint(x::Float64)
#     return round(Int, x)
# end

# function geo_coordinate(x, y)
#     PI = 3.141592

#     deg = nint(x)
#     m = x - deg
#     latitude = PI * (deg + 5.0 * m / 3.0) / 180.0

#     deg = nint(y)
#     m = y - deg
#     longitude = PI * (deg + 5.0 * m / 3.0) / 180.0

#     return latitude, longitude
# end

# function distance2D(xi, yi, xj, yj; dist="EUC_2D")
#     if dist == "EUC_2D"
#         xd = xi - xj
#         yd = yi - yj
#         return nint(sqrt(xd*xd + yd*yd))
#     elseif dist == "MAN_2D"
#         xd = abs(xi - xj)
#         yd = abs(yi - yj)
#         return nint(xd + yd)
#     elseif dist == "MAX_2D"
#         xd = abs(xi - xj)
#         yd = abs(yi - yj)
#         return max(nint(xd), nint(yd))
#     elseif dist == "GEO"
#         lat_i, long_i = geo_coordinate(xi, yi)
#         lat_j, long_j = geo_coordinate(xj, yj)
#         RRR = 6378.388
#         q1 = cos(long_i - long_j)
#         q2 = cos(lat_i - lat_j)
#         q3 = cos(lat_i + lat_j)
#         dij =  RRR * acos( 0.5*((1.0+q1)*q2 - (1.0-q1)*q3) ) + 1.0
#         return floor(Int, dij)
#     else
#         error("Distance function $dist is not supported.")
#     end
# end

# function dist_matrix(x::Vector{Float64}, y::Vector{Float64}; dist="EUC_2D")
#     n_nodes = length(x)
#     @assert length(x) == length(y)

#     M = Matrix{Int}(undef, n_nodes, n_nodes)

#     for i in 1:n_nodes
#         for j in i:n_nodes
#             if i == j
#                 M[i, j] = 0
#             else
#                 M[i, j] = distance2D(x[i], y[i], x[j], y[j]; dist=dist)
#                 M[j, i] = M[i, j]
#             end
#         end
#     end
#     return M
# end


text:
__precompile__()

module FluxionArrays

include("types.jl")
export Unknown
export ArrayExpr

include("show.jl")

include("constructors.jl")
export lazy_array
export unknown_array
export unknown_random_state

include("operators.jl")
export random_uniform
export random_normal

include("function.jl")
export array_function

include("backend.jl")
export ArrayBackend
export BackendArray
export BackendArrayFunction
export compile_array_function
export random_state
export default_array_backend
export to_backend
export from_backend

include("backend_julia.jl")
export JuliaBackend

end


text:
######## Stan diagnose example  ###########

using StanDiagnose

bernoulli_model = "
data {
  int<lower=0> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
    y ~ bernoulli(theta);
}
"
bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

stanmodel = DiagnoseModel("bernoulli", bernoulli_model;
  method=StanDiagnose.Diagnose(StanDiagnose.Gradient(epsilon=1e-6)));

rc = stan_diagnose(stanmodel; data=bernoulli_data);

if success(rc)
  diags = read_diagnose(stanmodel)
  println()
  display(diags)
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/wd.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const libwd_ = "cd") : (const libwd_ = "libcd")  # Name of IUP shared lib.

# These used to be in wd_h.jl (meanwhile removed)
const wdVectorTextTransform = cdVectorTextTransform
const wdActivate = cdActivate
# Skipping MacroDefinition: wdClip(mode)cdClip(mode)
# Skipping MacroDefinition: wdBegin(mode)cdBegin(mode)
# Skipping MacroDefinition: wdEnd()cdEnd();
# Skipping MacroDefinition: wdMM2Pixel(mm_dx,mm_dy,dx,dy)cdMM2Pixel(mm_dx,mm_dy,dx,dy)
# Skipping MacroDefinition: wdPixel2MM(dx,dy,mm_dx,mm_dy)cdPixel2MM(dx,dy,mm_dx,mm_dy)

function wdCanvasWindow(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasWindow, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasGetWindow(canvas::Ptr{cdCanvas}, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetWindow, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasViewport(canvas::Ptr{cdCanvas}, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasViewport, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Cint, Cint), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasGetViewport(canvas::Ptr{cdCanvas}, xmin::Ptr{Cint}, xmax::Ptr{Cint}, ymin::Ptr{Cint}, ymax::Ptr{Cint})
  ccall( (:wdCanvasGetViewport, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasWorld2Canvas(canvas::Ptr{cdCanvas}, xw::Cdouble, yw::Cdouble, xv::Ptr{Cint}, yv::Ptr{Cint})
  ccall( (:wdCanvasWorld2Canvas, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), canvas, xw, yw, xv, yv)
end
function wdCanvasWorld2CanvasSize(canvas::Ptr{cdCanvas}, hw::Cdouble, vw::Cdouble, hv::Ptr{Cint}, vv::Ptr{Cint})
  ccall( (:wdCanvasWorld2CanvasSize, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), canvas, hw, vw, hv, vv)
end
function wdCanvasCanvas2World(canvas::Ptr{cdCanvas}, xv::Integer, yv::Integer, xw::Ptr{Cdouble}, yw::Ptr{Cdouble})
  ccall( (:wdCanvasCanvas2World, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xv, yv, xw, yw)
end
function wdCanvasSetTransform(canvas::Ptr{cdCanvas}, sx::Cdouble, sy::Cdouble, tx::Cdouble, ty::Cdouble)
  ccall( (:wdCanvasSetTransform, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, sx, sy, tx, ty)
end
function wdCanvasGetTransform(canvas::Ptr{cdCanvas}, sx::Ptr{Cdouble}, sy::Ptr{Cdouble}, tx::Ptr{Cdouble}, ty::Ptr{Cdouble})
  ccall( (:wdCanvasGetTransform, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, sx, sy, tx, ty)
end
function wdCanvasTranslate(canvas::Ptr{cdCanvas}, dtx::Cdouble, dty::Cdouble)
  ccall( (:wdCanvasTranslate, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, dtx, dty)
end
function wdCanvasScale(canvas::Ptr{cdCanvas}, dsx::Cdouble, dsy::Cdouble)
  ccall( (:wdCanvasScale, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, dsx, dsy)
end
function wdCanvasClipArea(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasClipArea, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasGetClipArea(canvas::Ptr{cdCanvas}, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetClipArea, libwd_), Cint, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasIsPointInRegion(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasIsPointInRegion, libwd_), Cint, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end
function wdCanvasOffsetRegion(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasOffsetRegion, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end
function wdCanvasGetRegionBox(canvas::Ptr{cdCanvas}, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetRegionBox, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasHardcopy(canvas::Ptr{cdCanvas}, ctx::Ptr{Cint}, data::Ptr{Void}, draw_func::Ptr{Void})
  ccall( (:wdCanvasHardcopy, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cint}, Ptr{Void}, Ptr{Void}), canvas, ctx, data, draw_func)
end
function wdCanvasPixel(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, color::Clong)
  ccall( (:wdCanvasPixel, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Clong), canvas, x, y, color)
end
function wdCanvasMark(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasMark, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end
function wdCanvasLine(canvas::Ptr{cdCanvas}, x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdCanvasLine, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, x1, y1, x2, y2)
end
function wdCanvasVertex(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasVertex, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end
function wdCanvasRect(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasRect, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasBox(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasBox, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end
function wdCanvasArc(canvas::Ptr{cdCanvas}, xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdCanvasArc, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xc, yc, w, h, angle1, angle2)
end
function wdCanvasSector(canvas::Ptr{cdCanvas}, xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdCanvasSector, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xc, yc, w, h, angle1, angle2)
end
function wdCanvasChord(canvas::Ptr{cdCanvas}, xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdCanvasChord, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xc, yc, w, h, angle1, angle2)
end
function wdCanvasText(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdCanvasText, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, x, y, s)
end
function wdCanvasPutImageRect(canvas::Ptr{cdCanvas}, image::Ptr{Cint}, x::Cdouble, y::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRect, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cint}, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, image, x, y, xmin, xmax, ymin, ymax)
end
function wdCanvasPutImageRectRGB(canvas::Ptr{cdCanvas}, iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRectRGB, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, iw, ih, r, g, b, x, y, w, h, xmin, xmax, ymin, ymax)
end
function wdCanvasPutImageRectRGBA(canvas::Ptr{cdCanvas}, iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, a::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRectRGBA, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, iw, ih, r, g, b, a, x, y, w, h, xmin, xmax, ymin, ymax)
end
function wdCanvasPutImageRectMap(canvas::Ptr{cdCanvas}, iw::Integer, ih::Integer, index::Ptr{Cuchar}, colors::Ptr{Clong}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRectMap, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Ptr{Clong}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, iw, ih, index, colors, x, y, w, h, xmin, xmax, ymin, ymax)
end
function wdCanvasPutBitmap(canvas::Ptr{cdCanvas}, bitmap::Ptr{Cint}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble)
  ccall( (:wdCanvasPutBitmap, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cint}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, bitmap, x, y, w, h)
end
function wdCanvasLineWidth(canvas::Ptr{cdCanvas}, width::Cdouble)
  ccall( (:wdCanvasLineWidth, libwd_), Cdouble, (Ptr{cdCanvas}, Cdouble), canvas, width)
end
function wdCanvasFont(canvas::Ptr{cdCanvas}, type_face::String, style::Integer, size::Cdouble)
  ccall( (:wdCanvasFont, libwd_), Cint, (Ptr{cdCanvas}, Ptr{UInt8}, Cint, Cdouble), canvas, type_face, style, size)
end
function wdCanvasGetFont(canvas::Ptr{cdCanvas}, type_face::String, style::Ptr{Cint}, size::Ptr{Cdouble})
  ccall( (:wdCanvasGetFont, libwd_), Void, (Ptr{cdCanvas}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cdouble}), canvas, type_face, style, size)
end
function wdCanvasMarkSize(canvas::Ptr{cdCanvas}, size::Cdouble)
  ccall( (:wdCanvasMarkSize, libwd_), Cdouble, (Ptr{cdCanvas}, Cdouble), canvas, size)
end
function wdCanvasGetFontDim(canvas::Ptr{cdCanvas}, max_width::Ptr{Cdouble}, height::Ptr{Cdouble}, ascent::Ptr{Cdouble}, descent::Ptr{Cdouble})
  ccall( (:wdCanvasGetFontDim, libwd_), Void, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, max_width, height, ascent, descent)
end
function wdCanvasGetTextSize(canvas::Ptr{cdCanvas}, s::String, width::Ptr{Cdouble}, height::Ptr{Cdouble})
  ccall( (:wdCanvasGetTextSize, libwd_), Void, (Ptr{cdCanvas}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, s, width, height)
end
function wdCanvasGetTextBox(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetTextBox, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, x, y, s, xmin, xmax, ymin, ymax)
end
function wdCanvasGetTextBounds(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String, rect::Ptr{Cdouble})
  ccall( (:wdCanvasGetTextBounds, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}), canvas, x, y, s, rect)
end
function wdCanvasStipple(canvas::Ptr{cdCanvas}, w::Integer, h::Integer, fgbg::Ptr{Cuchar}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdCanvasStipple, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Cdouble, Cdouble), canvas, w, h, fgbg, w_mm, h_mm)
end
function wdCanvasPattern(canvas::Ptr{cdCanvas}, w::Integer, h::Integer, color::Ptr{Clong}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdCanvasPattern, libwd_), Void, (Ptr{cdCanvas}, Cint, Cint, Ptr{Clong}, Cdouble, Cdouble), canvas, w, h, color, w_mm, h_mm)
end
function wdCanvasVectorTextDirection(canvas::Ptr{cdCanvas}, x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdCanvasVectorTextDirection, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, x1, y1, x2, y2)
end
function wdCanvasVectorTextSize(canvas::Ptr{cdCanvas}, size_x::Cdouble, size_y::Cdouble, s::String)
  ccall( (:wdCanvasVectorTextSize, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, size_x, size_y, s)
end
function wdCanvasGetVectorTextSize(canvas::Ptr{cdCanvas}, s::String, x::Ptr{Cdouble}, y::Ptr{Cdouble})
  ccall( (:wdCanvasGetVectorTextSize, libwd_), Void, (Ptr{cdCanvas}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, s, x, y)
end
function wdCanvasVectorCharSize(canvas::Ptr{cdCanvas}, size::Cdouble)
  ccall( (:wdCanvasVectorCharSize, libwd_), Cdouble, (Ptr{cdCanvas}, Cdouble), canvas, size)
end
function wdCanvasVectorText(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdCanvasVectorText, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, x, y, s)
end
function wdCanvasMultiLineVectorText(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdCanvasMultiLineVectorText, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, x, y, s)
end
function wdCanvasGetVectorTextBounds(canvas::Ptr{cdCanvas}, s::String, x::Cdouble, y::Cdouble, rect::Ptr{Cdouble})
  ccall( (:wdCanvasGetVectorTextBounds, libwd_), Void, (Ptr{cdCanvas}, Ptr{UInt8}, Cdouble, Cdouble, Ptr{Cdouble}), canvas, s, x, y, rect)
end
function wdCanvasGetVectorTextBox(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetVectorTextBox, libwd_), Void, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, x, y, s, xmin, xmax, ymin, ymax)
end
function wdWindow(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdWindow, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end
function wdGetWindow(xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdGetWindow, libwd_), Void, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), xmin, xmax, ymin, ymax)
end
function wdViewport(xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdViewport, libwd_), Void, (Cint, Cint, Cint, Cint), xmin, xmax, ymin, ymax)
end
function wdGetViewport(xmin::Ptr{Cint}, xmax::Ptr{Cint}, ymin::Ptr{Cint}, ymax::Ptr{Cint})
  ccall( (:wdGetViewport, libwd_), Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), xmin, xmax, ymin, ymax)
end
function wdWorld2Canvas(xw::Cdouble, yw::Cdouble, xv::Ptr{Cint}, yv::Ptr{Cint})
  ccall( (:wdWorld2Canvas, libwd_), Void, (Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), xw, yw, xv, yv)
end
function wdWorld2CanvasSize(hw::Cdouble, vw::Cdouble, hv::Ptr{Cint}, vv::Ptr{Cint})
  ccall( (:wdWorld2CanvasSize, libwd_), Void, (Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), hw, vw, hv, vv)
end
function wdCanvas2World(xv::Integer, yv::Integer, xw::Ptr{Cdouble}, yw::Ptr{Cdouble})
  ccall( (:wdCanvas2World, libwd_), Void, (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), xv, yv, xw, yw)
end
function wdClipArea(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdClipArea, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end
function wdGetClipArea(xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdGetClipArea, libwd_), Cint, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), xmin, xmax, ymin, ymax)
end
function wdGetClipPoly(n::Ptr{Cint})
  ccall( (:wdGetClipPoly, libwd_), Ptr{Cdouble}, (Ptr{Cint},), n)
end
function wdPointInRegion(x::Cdouble, y::Cdouble)
  ccall( (:wdPointInRegion, libwd_), Cint, (Cdouble, Cdouble), x, y)
end
function wdOffsetRegion(x::Cdouble, y::Cdouble)
  ccall( (:wdOffsetRegion, libwd_), Void, (Cdouble, Cdouble), x, y)
end
function wdRegionBox(xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdRegionBox, libwd_), Void, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), xmin, xmax, ymin, ymax)
end
function wdHardcopy(ctx::Ptr{Cint}, data::Ptr{Void}, cnv::Ptr{Cint}, draw_func::Ptr{Void})
  ccall( (:wdHardcopy, libwd_), Void, (Ptr{Cint}, Ptr{Void}, Ptr{Cint}, Ptr{Void}), ctx, data, cnv, draw_func)
end
function wdPixel(x::Cdouble, y::Cdouble, color::Clong)
  ccall( (:wdPixel, libwd_), Void, (Cdouble, Cdouble, Clong), x, y, color)
end
function wdMark(x::Cdouble, y::Cdouble)
  ccall( (:wdMark, libwd_), Void, (Cdouble, Cdouble), x, y)
end
function wdLine(x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdLine, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble), x1, y1, x2, y2)
end
function wdVertex(x::Cdouble, y::Cdouble)
  ccall( (:wdVertex, libwd_), Void, (Cdouble, Cdouble), x, y)
end
function wdRect(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdRect, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end
function wdBox(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdBox, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end
function wdArc(xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdArc, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), xc, yc, w, h, angle1, angle2)
end
function wdSector(xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdSector, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), xc, yc, w, h, angle1, angle2)
end
function wdChord(xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdChord, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), xc, yc, w, h, angle1, angle2)
end
function wdText(x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdText, libwd_), Void, (Cdouble, Cdouble, Ptr{UInt8}), x, y, s)
end
function wdPutImageRect(image::Ptr{Cint}, x::Cdouble, y::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRect, libwd_), Void, (Ptr{Cint}, Cdouble, Cdouble, Cint, Cint, Cint, Cint), image, x, y, xmin, xmax, ymin, ymax)
end
function wdPutImageRectRGB(iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRectRGB, libwd_), Void, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), iw, ih, r, g, b, x, y, w, h, xmin, xmax, ymin, ymax)
end
function wdPutImageRectRGBA(iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, a::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRectRGBA, libwd_), Void, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), iw, ih, r, g, b, a, x, y, w, h, xmin, xmax, ymin, ymax)
end
function wdPutImageRectMap(iw::Integer, ih::Integer, index::Ptr{Cuchar}, colors::Ptr{Clong}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRectMap, libwd_), Void, (Cint, Cint, Ptr{Cuchar}, Ptr{Clong}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), iw, ih, index, colors, x, y, w, h, xmin, xmax, ymin, ymax)
end
function wdPutBitmap(bitmap::Ptr{Cint}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble)
  ccall( (:wdPutBitmap, libwd_), Void, (Ptr{Cint}, Cdouble, Cdouble, Cdouble, Cdouble), bitmap, x, y, w, h)
end
function wdLineWidth(width::Cdouble)
  ccall( (:wdLineWidth, libwd_), Cdouble, (Cdouble,), width)
end
function wdFont(type_face::Integer, style::Integer, size::Cdouble)
  ccall( (:wdFont, libwd_), Void, (Cint, Cint, Cdouble), type_face, style, size)
end
function wdGetFont(type_face::Ptr{Cint}, style::Ptr{Cint}, size::Ptr{Cdouble})
  ccall( (:wdGetFont, libwd_), Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), type_face, style, size)
end
function wdMarkSize(size::Cdouble)
  ccall( (:wdMarkSize, libwd_), Cdouble, (Cdouble,), size)
end
function wdFontDim(max_width::Ptr{Cdouble}, height::Ptr{Cdouble}, ascent::Ptr{Cdouble}, descent::Ptr{Cdouble})
  ccall( (:wdFontDim, libwd_), Void, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), max_width, height, ascent, descent)
end
function wdTextSize(s::String, width::Ptr{Cdouble}, height::Ptr{Cdouble})
  ccall( (:wdTextSize, libwd_), Void, (Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), s, width, height)
end
function wdTextBox(x::Cdouble, y::Cdouble, s::String, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdTextBox, libwd_), Void, (Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), x, y, s, xmin, xmax, ymin, ymax)
end
function wdTextBounds(x::Cdouble, y::Cdouble, s::String, rect::Ptr{Cdouble})
  ccall( (:wdTextBounds, libwd_), Void, (Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}), x, y, s, rect)
end
function wdStipple(w::Integer, h::Integer, stipple::Ptr{Cuchar}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdStipple, libwd_), Void, (Cint, Cint, Ptr{Cuchar}, Cdouble, Cdouble), w, h, stipple, w_mm, h_mm)
end
function wdPattern(w::Integer, h::Integer, pattern::Ptr{Clong}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdPattern, libwd_), Void, (Cint, Cint, Ptr{Clong}, Cdouble, Cdouble), w, h, pattern, w_mm, h_mm)
end
function wdVectorTextDirection(x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdVectorTextDirection, libwd_), Void, (Cdouble, Cdouble, Cdouble, Cdouble), x1, y1, x2, y2)
end
function wdVectorTextSize(size_x::Cdouble, size_y::Cdouble, s::String)
  ccall( (:wdVectorTextSize, libwd_), Void, (Cdouble, Cdouble, Ptr{UInt8}), size_x, size_y, s)
end
function wdGetVectorTextSize(s::String, x::Ptr{Cdouble}, y::Ptr{Cdouble})
  ccall( (:wdGetVectorTextSize, libwd_), Void, (Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), s, x, y)
end
function wdVectorCharSize(size::Cdouble)
  ccall( (:wdVectorCharSize, libwd_), Cdouble, (Cdouble,), size)
end
function wdVectorText(x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdVectorText, libwd_), Void, (Cdouble, Cdouble, Ptr{UInt8}), x, y, s)
end
function wdMultiLineVectorText(x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdMultiLineVectorText, libwd_), Void, (Cdouble, Cdouble, Ptr{UInt8}), x, y, s)
end
function wdGetVectorTextBounds(s::String, x::Cdouble, y::Cdouble, rect::Ptr{Cdouble})
  ccall( (:wdGetVectorTextBounds, libwd_), Void, (Ptr{UInt8}, Cdouble, Cdouble, Ptr{Cdouble}), s, x, y, rect)
end


text:
#!/usr/local/bin/julia
# Assumes agent goes first always for now
using OpenAIGym
using go

const debug = false

model_dir = joinpath(splitdir(@__FILE__)[1], "../models")

function usage()
    println(STDOUT, "Usage: ./openai_gym.jl exp_name rounds {random|keras [model name]}")
    exit()
end

nargs = length(ARGS)
exp_name = ARGS[1]
nrounds = parse(Int, ARGS[2])

if nargs != 3
    if ARGS[3] == "random"
        policy = go.RandomPolicy()
    elseif ARGS[3] == "keras"
        if nargs != 4
            usage()
        end
        policy = go.KerasNetwork(model_dir, ARGS[4])
    else
        usage()
    end
else
    usage()
end

function getempty(board)
    reshape(board[3, :, :], (go.N, go.N)).'
end

function get_opponent_move(last_empty, cur_empty, last_move::Int=0)
    # @show last_empty
    # @show cur_empty
    move = 0
    for i = 1:(go.N * go.N)
        if last_empty[i] == 1 && cur_empty[i] == 0 && i != last_move
            move = i
        end
    end
    if move == 0
        @show cur_empty
        @show last_empty
        @show last_move
        return go.PASS_MOVE
    else
        return from_gym_coords(go.pointindex(move))
    end
end

"""
gym indexes in row order from upper left of the board
 while we index from lower left along columns.
Together with from_gym_coords, transform between them.
"""
function to_gym_coords(p::go.Point)
    go.Point(p[2], go.N-p[1] + 1)
end

function from_gym_coords(p::go.Point)
    go.Point(go.N - p[2] + 1, p[1])
end

function run_episode()
    state = reset(env)
    localboard = go.Board()

    # TODO: The coordinate transforms are hacky and uncommented, but they work.
    last_empty = getempty(state.observation)
    opponent_move = go.EMPTY_MOVE # Init here so it stays in scope for after loop
    while !state.done
        # Our move
        if opponent_move == go.PASS_MOVE
            move = go.PASS_MOVE
            state = step(env, 81)
            break
        else
            move = go.choose_move(localboard, policy)
            state = step(env, go.linearindex(to_gym_coords(move)) - 1)
        end
        go.play_move(localboard, move)

        # Opponent move
        board = state.observation
        empty = getempty(board)
        opponent_move = get_opponent_move(last_empty, empty, go.linearindex(to_gym_coords(move)))

        if debug
            @show move
            @show opponent_move
        end

        go.play_move(localboard, opponent_move)
        last_empty = empty

        if debug
            println(localboard)
            display(env)
        end
        flush(STDOUT)
        flush(STDERR)
    end
    println(localboard)
    display(env)
    @show opponent_move
    @show state.reward
    @show go.calculate_score(localboard)
end

env = Env("Go$(go.N)x$(go.N)-v0")
env.env[:monitor][:start](exp_name)

for episode in 1:nrounds
    println("Episode $(episode)/$(nrounds)")
    run_episode()
end

env.env[:monitor][:close]()


text:
# Autogenerated wrapper script for PCRE_jll for armv6l-linux-gnueabihf-cxx03
export libpcre

JLLWrappers.@generate_wrapper_header("PCRE")
JLLWrappers.@declare_library_product(libpcre, "libpcre.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpcre,
        "lib/libpcre.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function value_estimate(π::DiscreteNetwork, px, s) where {T <: DiscreteNetwork}
        pdfs = exp.(logpdf(px, s))
        sum(value(π, s) .* pdfs, dims=1)
end

#NOTE: Currently all of these assume that we get a reward (cost) ONCE at the end of an episode

function E_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        px = 𝒫[:px]
        if 𝒫[:use_likelihood_weights]
                return 𝒟[:likelihoodweight] .* (𝒟[:done] .* 𝒟[:r] .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp]))
        else
                return 𝒟[:done] .* 𝒟[:r] .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp])
        end
end

function CDF_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        rα = 𝒫[:rα][1]
        px = 𝒫[:px]

        if 𝒫[:use_likelihood_weights]
                return 𝒟[:likelihoodweight] .* (𝒟[:done] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp]))
        else
                return 𝒟[:done] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp])
        end
end

function CVaR_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        rα = 𝒫[:rα][1]
        px = 𝒫[:px]
        if 𝒫[:use_likelihood_weights]
                return 𝒟[:likelihoodweight] .* (𝒟[:done] .* 𝒟[:r] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp]))
        else
                return 𝒟[:done] .* 𝒟[:r] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp])
        end
end

function E_VaR_CVaR_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        [CDF_target(π.networks[1], 𝒫, 𝒟, γ; kwargs...), CVaR_target(π.networks[2], 𝒫, 𝒟, γ; kwargs...), E_target(π.networks[3], 𝒫, 𝒟, γ; kwargs...)]
end

function VaR_CVaR_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        [CDF_target(π.networks[1], 𝒫, 𝒟, γ; kwargs...), CVaR_target(π.networks[2], 𝒫, 𝒟, γ; kwargs...)]
end

function CERL_Discrete(;π::MixtureNetwork,
                        S,
                        N,
                        px,
                        prioritized=true,
                        use_likelihood_weights=true,
                        α,
                        𝒫=(;),
                        buffer_size=N,
                        ΔN=4,
                        pre_train_callback,
                        π_explore,
                        c_opt::NamedTuple=(;),
                        log::NamedTuple=(;),
                        c_loss,
                        kwargs...)

                    𝒫 = (;px, rα=[NaN], α, use_likelihood_weights, 𝒫...)
                    required_columns=[:logprob, :likelihoodweight]
                    agent = PolicyParams(π=π, π_explore=π_explore, π⁻=deepcopy(π), pa=px)
                    OffPolicySolver(;agent=agent,
                                     S=S,
                                     log=LoggerParams(;dir="log/cerl_dqn", period=100, fns=[log_episode_averages([:r], 100)], log...),
                                     𝒫=𝒫,
                                     N=N,
                                     ΔN=ΔN,
                                     pre_train_callback=pre_train_callback,
                                     buffer=ExperienceBuffer(S, agent.space, buffer_size, required_columns, prioritized=prioritized),
                                     c_opt = TrainingParams(;loss=c_loss, name="critic_", epochs=ΔN, c_opt...),
                                     target_fn=VaR_CVaR_target,
                                     kwargs...)
end


text:
using SuperTB, PyPlot, ProgressMeter, LinearAlgebra

# set up a model (amorphous system with AII-class hopping)
hopfun(ri,rj,per) = hopAgarwala(ri, rj, per,
                                1.0/6,  # R
                                -.5,    # M
                                0.25,   # t₂
                                0.5,    # λ
                                1.0/64) # a
tb = TBModel(hopfun)
tb.periodicity = trues(2) # periodic boundary conditions
tb.basisdim    = 2

# run over this parameter range
sitessqrt  = 24
sites      = sitessqrt^2
Nensembles = 5
getbottatindex = sites .+ (-(2*sitessqrt+1):(2*sitessqrt+1))

# preallocation
par  = zeros(sites*tb.basisdim, Nensembles)      # participation ratios
εeig = zeros(sites*tb.basisdim, Nensembles)      # energy eigenvalues
β    = zeros(size(getbottatindex,1), Nensembles) # bott indices

close("all")
# compute and visualize during...
@showprogress 1 "Calculating ensembles..." for i in 1:Nensembles
    tb.pos = rand(sites,2)
    # solve the system
    H = Hermitian(Matrix(setHopping!(tb)))
    F = eigen(H)

    # compute participation ratios & plot
    par[:,i]  = parratio(F.vectors)
	εeig[:,i] = F.values
	figure()
	plot(εeig[:,i], par[:,i], ".")

    # compute bott indices & plot
    β[:,i] = bott(F.vectors,  tb,getbottatindex)
	plot(εeig[getbottatindex,i], β[:,i]./10, ".") # scale to fit in same plot...
end

#PyPlot.svg(true)
#savefig("sdos_ensembleaveraged")

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
##  Eratosthenes' prime number sieve
function primesieve(n::Int64)
    if n <= 1
        return Integer[]
    end

    p = [i for i in 1:2:n]
    q = length(p)
    p[1] = 2
    if n >= 9
        for k in 3:2:isqrt(n)
            if p[(k+1)÷2] != 0
                p[(k*k+1).÷2:k:q] .= 0
            end
        end
    end

    return p[p .> 0]
end

##  Find all prime numbers in  a given interval
function primes2(n::Int64, m::Int64)
    if n > m
        error("Argument 'm' must be bigger than 'n'.")
    end
    if m <= 1000
        P = primesieve(m)
        return P[P .>= n]
    end

    myPrimes = primesieve(isqrt(m))
    N = [n:m]
    l = length(N)  # m-n+1
    A = zeros(Int8, l)
    if n == 1
        A[1] = -1
    end

    for p in myPrimes
        r = n % p
        if r == 0
            i = 1
        else
            i = p - r + 1
        end
        if i <= l && N[i] == p
            i = i + p
        end
        while i <= l
            A[i] = 1
            i = i + p
        end
    end
    return N[A .== 0]
end

##  Return unique prime factors of n sorted
function primefactors(n::Int64)
    if n <= 0
        error("Argument 'n' to be factored must be a positive integer.")
    end

    f = Primes.factor(n)
    return sort(collect(keys(f)))
end

##  Find prime number following n
function nextprime2(n::Int64)
    if n <= 1; return 2; end
    if n == 2; return 3; end
    if iseven(n)
        n += 1
    else
        n += 2
    end
    if isprime(n); return(n); end

    if mod(n, 3) == 1
        a = 4; b = 2
    elseif mod(n, 3) == 2
        a = 2; b = 4
    else
        n += 2
        a = 2; b = 4
    end

    p = n
    while !isprime(p)
        p += a
        if isprime(p); break; end
        p += b
    end

    return p
end

## Find prime number preceeding n
function prevprime2(n::Int64)
    if n <= 2
        return Array(typeof(n), 0)
    elseif n <= 3
        return 2
    end

    if iseven(n)
        n -= 1
    else
        n -= 2
    end
    if isprime(n); return n; end

    if mod(n, 3) == 1
        a = 2; b = 4
    elseif mod(n, 3) == 2
        a = 4; b = 2
    else
        n -= 2
        a = 2; b = 4
    end

    p = n
    while !isprime(p)
        p -= a
        if isprime(p); break; end
        p -= b
    end

    return p
end

##  Find all twin primes
function twinprimes(n::Int64, m::Int64)
    P = primes2(n, m)
    inds = findall(diff(P) .== 2)

    return hcat(P[inds], P[inds.+1])
end

##  Coprimality
function coprime(n::Int64, m::Int64)
    if n == 0 && m == 0
        return false
    end

    if gcd(n, m) > 1
        false
    else
        true
    end
end

##  Solve linear modulo equations
function linmod(a::Int64, b::Int64, m::Int64)
    m > 1 || error("Argument 'm' must be an integer greater 1")
    T = typeof(m)
    g, e, f = gcdx(a, m)

    x = T[]
    if mod(b, g) == 0
        x0 = mod(e * div(b, g), m)
        for i = 0:(g-1)
            x = [x, mod(x0 + i * div(m, g) , m)]
        end
    end

    return x
end

##  Order of the element n (in the ring) modulo m
function ordermod(n::Int64, m::Int64)
    if n <= 0 || m <= 0
        error("Arguments 'n' and 'm' must be positive integers.")
    end
    if m == 1 || gcd(n, m) > 1; return 0; end
    if n == 1; return 1; end

    r = mod(n, m)
    if r == 0; return 0; end

    k = 1
    while r != 1
        r = mod(n*r, m)
        k += 1
    end

    return k
end

##  Find a primitive root modulo m
function primroot(m::Int64)
    if !isprime(m)
        error("Argument 'm' must be a prime number")
    end
    if m == 2; return 1; end

    P = primefactors(m-1)
    for r = 2:(m-1)
        not_found = true
        for p in P
            if powermod(r, div(m-1, p), m) == 1
                not_found = false
            end
        end
        if not_found
            return r
        end
    end

    return 0
end

##  Euler's Phi (or: totient) function
function eulerphi(n::Int64)
    if n <= 0
        error("Argument 'n' must be a (positive) natural number.")
    end

    m = n
    for p in primefactors(n)    # must be unique
        m = m - div(m, p)       # m = m * (1 - 1/p)
    end

    return Integer(round(m))
end


text:
"""
    couple_models(aero::LiftingLine, stru::GEBT, dyn::RigidBody)

Create an aerostructural model using a lifting line aerodynamic model coupled
with a geometrically exact beam theory and a rigid-body dynamics model.  This
model introduces additional parameters corresponding to the freestream air density
``\\rho_\\infty``, followed by the external loads ``F_{x,i}, F_{y,i}, F_{z,i},
M_{x,i}, M_{y,i}, M_{z,i}`` or displacements ``u_{x,i}, u_{y,i}, u_{z,i},
\\theta_{x,i}, \\theta_{y,i}, \\theta_{z,i}`` for each node, followed by the
constant distributed loads ``f_{x,i}, f_{y,i}, f_{z,i}, m_{x,i}, m_{y,i},
m_{z,i}`` applied on each beam element (excluding aerodynamic loads).

**NOTE: When using this model, the local frame for each beam element should be
oriented with the x-axis along the beam's axis, the y-axis forward, and the
z-axis normal to the surface**
"""
couple_models(aero::LiftingLine, stru::GEBT, dyn::RigidBody) = (aero, stru, dyn)

# --- traits --- #

function number_of_additional_parameters(aero::LiftingLine, stru::GEBT, dyn::RigidBody)
    return 1 + 6*length(stru.icol_point) + 6*length(stru.icol_elem)
end

coupling_inplaceness(::Type{<:LiftingLine}, ::Type{<:GEBT}, ::Type{<:RigidBody}) = InPlace()

function coupling_rate_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}, ::Type{<:RigidBody}) where {NA,TA}
    model_types = TA.parameters
    if all(islinear.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

function coupling_state_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}, ::Type{<:RigidBody}) where {NA,TA}
    model_types = TA.parameters
    if all(islinear.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

coupling_parameter_jacobian_type(::Type{<:LiftingLine}, ::Type{<:GEBT}, ::Type{<:RigidBody}) = Nonlinear()

function coupling_time_gradient_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}, ::Type{<:RigidBody}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isconstant.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

# --- methods --- #

function get_coupling_inputs!(y, aero::LiftingLine{NA,TA}, stru::GEBT, dyn::RigidBody, dx, x, p, t) where {NA,TA}

    # extract number of state variables, inputs, and parameters
    nua = number_of_states(aero) # number of aerodynamic states
    nus = number_of_states(stru) # number of structural states
    nud = number_of_states(dyn) # number of rigid body states
    nya = number_of_inputs(aero) # number of aerodynamic inputs
    nys = number_of_inputs(stru) # number of structural inputs
    nyd = number_of_inputs(dyn) # number of rigid body inputs
    npa = number_of_parameters(aero) # number of aerodynamic parameters
    nps = number_of_parameters(stru) # number of structural parameters
    npd = number_of_parameters(dyn) # number of rigid body parameters
    npadd = number_of_additional_parameters(aero, stru, dyn) # number of additional parameters

    # get indices for state variables, inputs, and parameters
    iua = 1:nua # indices of aerodynamic states
    iya = 1:nya # indices of aerodynamic inputs
    ipa = 1:npa # indices of aerodynamic parameters
    ius = nua + 1 : nua + nus # indices of structural states
    iys = nya + 1 : nya + nys # indices of structural inputs
    ips = npa + 1 : npa + nps # indices of structural parameters
    iud = nua + nus + 1 : nua + nus + nud # indices of rigid body states
    iyd = nya + nys + 1 : nya + nys + nyd # indices of rigid body inputs
    ipd = npa + nps + 1 : npa + nps + npd # indices of rigid body parameters
    ipadd = npa + nps + npd + 1 : npa + nps + npd + npadd # indices of additional parameters
    iuas = state_indices(aero.models) # indices of aerodynamic states for each section
    iyas = input_indices(aero.models) # indices of aerodynamic inputs for each section
    ipas = parameter_indices(aero.models) # indices of aerodynamic inputs for each section

    # separate state variables, inputs, and parameters
    dua = view(dx, iua) # aerodynamic rates
    ua = view(x, iua) # aerodynamic states
    ya = view(y, iya) # aerodynamic inputs
    pa = view(p, ipa) # aerodynamic parameters
    dus = view(dx, ius) # structural rates
    us = view(x, ius) # structural states
    ys = view(y, iys) # structural inputs
    ps = view(p, ips) # structural parameters
    dud = view(dx, iud) # rigid body rates
    ud = view(x, iud) # rigid body states
    yd = view(y, iyd) # rigid body inputs
    pd = view(p, ipd) # rigid body parameters
    padd = view(p, ipadd) # additional parameters
    duas = view.(Ref(dua), iuas) # aerodynamic rates for each section
    uas = view.(Ref(ua), iuas) # aerodynamic states for each section
    yas = view.(Ref(ya), iyas) # aerodynamic inputs for each section
    pas = view.(Ref(pa), ipas) # aerodynamic parameters for each section

    # number of points and elements
    npoint = length(stru.icol_point)
    nelem = length(stru.icol_elem)

    # rigid body states
    xr, yr, zr, ϕr, θr, ψr, ur, vr, wr, pr, qr, rr = ud

    # rigid body state rates
    dxr, dyr, dzr, dϕr, dθr, dψr, dur, dvr, dwr, dpr, dqr, drr = dud

    # global parameters
    ρ = padd[1]

    # body linear and angular velocity
    V = SVector(ur, vr, wr)
    Ω = SVector(pr, qr, rr)

    # rigid body linear and angular acceleration
    dV = SVector(dur, dvr, dwr)
    dΩ = SVector(dpr, dqr, drr)

    # construct assembly from structural parameters
    assembly = gxbeam_assembly(ps, npoint, nelem, stru.start, stru.stop)

    # initialize total forces and moments
    Ftot = @SVector zeros(3)
    Mtot = @SVector zeros(3)

    # initialize rigid body properties
    mass = 0.0
    Ir = @SMatrix zeros(3,3)

    # save prescribed point loads/displacements
    for ip = 1:npoint
        yoff = 6*(ip-1)
        poff = 1 + 6*(ip-1)
        ys[yoff+1] = padd[poff+1]
        ys[yoff+2] = padd[poff+2]
        ys[yoff+3] = padd[poff+3]
        ys[yoff+4] = padd[poff+4]
        ys[yoff+5] = padd[poff+5]
        ys[yoff+6] = padd[poff+6]
    end

    # add prescribed loads to total forces and moments
    for ip = 1:npoint
        # point location
        point = assembly.points[ip]
        # is prescribed condition a force?
        prescribed_forces = SVector{6}(view(stru.displacement, :, ip)) .== false
        # point loads
        poff = 1 + 6*(ip-1)
        F1 = ifelse(prescribed_forces[1], padd[poff+1], zero(padd[poff+1]))
        F2 = ifelse(prescribed_forces[2], padd[poff+2], zero(padd[poff+2]))
        F3 = ifelse(prescribed_forces[3], padd[poff+3], zero(padd[poff+3]))
        M1 = ifelse(prescribed_forces[4], padd[poff+4], zero(padd[poff+4]))
        M2 = ifelse(prescribed_forces[5], padd[poff+5], zero(padd[poff+5]))
        M3 = ifelse(prescribed_forces[6], padd[poff+6], zero(padd[poff+6]))
        Fp = SVector(F1, F2, F3)
        Mp = SVector(M1, M2, M3)
        Ftot += Fp
        Mtot += cross(point, Fp) + Mp
    end

    # calculate aerodynamic inputs, distributed loads, and element properties
    for i = 1:NA
        # models for this section
        section_aero = aero.models[i]
        section_stru = LiftingLineSection()
        section_models = (section_aero, section_stru)

        # model dimensions for this section
        Nuai = number_of_states(section_aero)
        Nyai = number_of_inputs(section_aero)
        Npai = number_of_parameters(section_aero)
        Nusi = number_of_states(section_stru)
        Nysi = number_of_inputs(section_stru)
        Npsi = number_of_parameters(section_stru)

        # section properties
        icol = stru.icol_elem[i]
        element = assembly.elements[i]

        u_elem = SVector(us[icol], us[icol+1], us[icol+2]) # linear displacement
        θ_elem = SVector(us[icol+3], us[icol+4], us[icol+5]) # angular displacement
        F_elem = SVector(us[icol+6], us[icol+7], us[icol+8]) .* stru.force_scaling # internal force
        M_elem = SVector(us[icol+9], us[icol+10], us[icol+11]) .* stru.force_scaling # internal moment
        P_elem = SVector(us[icol+12], us[icol+13], us[icol+14]) .* stru.mass_scaling # linear momentum
        H_elem = SVector(us[icol+15], us[icol+16], us[icol+17]) .* stru.mass_scaling # angular momentum

        dF_elem = SVector(dus[icol+6], dus[icol+7], dus[icol+8]) .* stru.force_scaling
        dM_elem = SVector(dus[icol+9], dus[icol+10], dus[icol+11]) .* stru.force_scaling
        dP_elem = SVector(dus[icol+12], dus[icol+13], dus[icol+14]) .* stru.mass_scaling
        dH_elem = SVector(dus[icol+15], dus[icol+16], dus[icol+17]) .* stru.mass_scaling

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = GXBeam.rotation_parameter_scaling(θ_elem)
        θ_elem *= scaling # angular displacement (Wiener-Milenkovic parameters)

        # length and location of the element
        ΔL = element.L # element length
        pe = element.x + u_elem # location of element (in body frame)

        # local to body transformation
        CtCab = GXBeam.get_C(θ_elem)'*element.Cab # local to body transformation

        # local structural to aerodynamic transformation
        R = @SMatrix [0 -1 0; 1 0 0; 0 0 1]

        # local freestream linear and angular velocities
        vi = -R*(CtCab'*V + GXBeam.element_linear_velocity(element, P_elem, H_elem)) # local linear freestream velocity
        ωi = R*GXBeam.element_angular_velocity(element, P_elem, H_elem) # local angular freestream velocity

        # local freestream linear and angular accelerations
        dvi = -R*(CtCab'*dV + CtCab'*cross(dΩ, pe) + GXBeam.element_linear_velocity(element, dP_elem, dH_elem))
        dωi = R*CtCab'*dΩ + R*GXBeam.element_angular_velocity(element, dP_elem, dH_elem)

        poff = 3*npoint + 36*(i-1)
        μ = ps[poff + 31] # mass per unit length
        xm2 = ps[poff + 32] # center of mass location
        xm3 = ps[poff + 33] # center of mass location
        i22 = ps[poff + 34] # rotational inertia
        i33 = ps[poff + 35] # rotational inertia
        i23 = ps[poff + 36] # product of inertia

        me = ΔL*μ # mass
        Ie = ΔL*(@SMatrix [i22+i33 0 0; 0 i22 -i23; 0 -i23 i33]) # local frame inertia
        Ie = (CtCab*R')*Ie*(R*CtCab') # body frame inertia
        Ie = me*(pe'*pe*I - pe*pe') + Ie # body frame inertia about origin

        # NOTE: We assume that the local beam frame y-axis is oriented towards
        # the leading edge and the z-axis is oriented up
        # NOTE: We assume the local aerodynamic frame is oriented with the x-axis
        # in the chordwise direction and the z-axis in the (airfoil) normal
        # direction.

        # section rate variables
        duai = SVector{Nuai}(duas[i]) # aerodynamic rate variables
        dusi = vcat(dvi, dωi) # structural rate variables
        dui = vcat(duai, dusi)

        # section state variables
        uai = SVector{Nuai}(uas[i]) # aerodynamic state variables
        usi = vcat(vi, ωi) # structural state variables
        ui = vcat(uai, usi)

        # section parameters
        pai = SVector{Npai}(pas[i]) # aerodynamic parameters
        psi = SVector(ρ) # structural parameters
        pi = vcat(pai, psi)

        # section inputs
        yi = get_coupling_inputs(section_models, dui, ui, pi, t)

        # separate inputs
        yai = view(yi, 1:Nyai)
        ysi = view(yi, Nyai + 1 : Nyai + Nysi)

        # save local inputs
        y[iya[iyas[i]]] .= yai # aerodynamic inputs

        # section aerodynamic loads (in body frame)
        fi = CtCab*R'*SVector(ysi[1], ysi[2], ysi[3])
        mi = CtCab*R'*SVector(ysi[4], ysi[5], ysi[6])

        # add apparent forces due to body frame linear and angular acceleration
        fi += me*dV - cross(me*dΩ, pe)
        mi += me*dΩ

        # add constant distributed loads (in body frame)
        poff = 1 + 6*npoint + 6*(i-1)
        fi += SVector(padd[poff+1], padd[poff+2], padd[poff+3])
        mi += SVector(padd[poff+4], padd[poff+5], padd[poff+6])

        # save distributed loads for this element (in body frame)
        yoff = 6*npoint + 6*(i-1)
        y[iys[yoff+1 : yoff+3]] = fi
        y[iys[yoff+4 : yoff+6]] = mi

        # add distributed loads to total forces and moments
        Ftot += ΔL*fi
        Mtot += cross(pe, ΔL*fi) + ΔL*mi

        # add element mass and inertia to total mass and inertia
        mass += me
        Ir += Ie
    end

    # save body frame linear and angular velocities
    yoff = 6*npoint + 6*nelem
    y[iys[yoff+1]] = ur
    y[iys[yoff+2]] = vr
    y[iys[yoff+3]] = wr
    y[iys[yoff+4]] = pr
    y[iys[yoff+5]] = qr
    y[iys[yoff+6]] = rr

    # save rigid body model inputs
    y[iyd[1]] = mass
    y[iyd[2]] = Ixx = Ir[1,1]
    y[iyd[3]] = Iyy = Ir[2,2]
    y[iyd[4]] = Izz = Ir[3,3]
    y[iyd[5]] = Ixz = Ir[1,3]
    y[iyd[6]] = Ixy = Ir[1,2]
    y[iyd[7]] = Iyz = Ir[2,3]
    y[iyd[8:10]] = Ftot
    y[iyd[11:13]] = Mtot

    return y
end

# --- Performance Overloads --- #

# TODO:

# --- Convenience Methods --- #

function set_additional_parameters!(padd, model::LiftingLine, stru::GEBT,
    dyn::RigidBody; rho, point_conditions, element_loads)

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    padd[1] = rho

    for ip = 1:np
        padd[1+6*(ip-1)+1] = point_conditions[6*(ip-1)+1]
        padd[1+6*(ip-1)+2] = point_conditions[6*(ip-1)+2]
        padd[1+6*(ip-1)+3] = point_conditions[6*(ip-1)+3]
        padd[1+6*(ip-1)+4] = point_conditions[6*(ip-1)+4]
        padd[1+6*(ip-1)+5] = point_conditions[6*(ip-1)+5]
        padd[1+6*(ip-1)+6] = point_conditions[6*(ip-1)+6]
    end

    for ie = 1:ne
        padd[1+6*np+6*(ie-1)+1] = element_loads[6*(ie-1)+1]
        padd[1+6*np+6*(ie-1)+2] = element_loads[6*(ie-1)+2]
        padd[1+6*np+6*(ie-1)+3] = element_loads[6*(ie-1)+3]
        padd[1+6*np+6*(ie-1)+4] = element_loads[6*(ie-1)+4]
        padd[1+6*np+6*(ie-1)+5] = element_loads[6*(ie-1)+5]
        padd[1+6*np+6*(ie-1)+6] = element_loads[6*(ie-1)+6]
    end

    return padd
end

function separate_additional_parameters(model::LiftingLine, stru::GEBT,
    dyn::RigidBody, padd)

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    rho = padd[1]

    point_conditions = reshape(view(padd, 1 + 1 : 1 + 6*np), 6, np)

    element_loads = reshape(view(padd, 1 + 6*np + 1 : 1 + 6*np + 6*ne), 6, ne)

    return (rho = rho, point_conditions = point_conditions, element_loads = element_loads)
end


text:
####
#### Data assembling
####

function plotdata(x::AVM{<:CanMiss{<:Real}})
    isempty(x) && throw(ArgumentError("Cannot display an empty vector"))
    hasmissing = Missing <: eltype(x)
    hasmissing = hasmissing || any(isinf, x)
    hasmissing = hasmissing || any(isnan, x)
    return (data=zip(1:size(x, 1), eachcol(x)...),
            hasmissing=hasmissing,
            nobj=size(x, 2))
end
# NOTE: these typechecks within the function body is to avoid clashes/ambiguity with plotdata(x)
function plotdata(x, ys...)
    isempty(x) && throw(ArgumentError("Cannot display empty vectors."))
    x isa AV{<:CanMiss{<:Real}} || throw(ArgumentError("x has un-handled type $(typeof(x))"))
    nobj = 0
    hasmissing = Missing <: eltype(x)
    hasmissing = hasmissing || any(isinf, x)
    hasmissing = hasmissing || any(isnan, x)
    for y ∈ ys
        y isa AVM{<:CanMiss{<:Real}} || throw(ArgumentError("y has un-handled type $(typeof(y))"))
        size(y, 1) == length(x) || throw(DimensionMismatch("y data must match x"))
        nobj += size(y, 2)
        hasmissing = hasmissing || Missing <: eltype(y)
        hasmissing = hasmissing || any(isinf, y)
        hasmissing = hasmissing || any(isnan, y)
    end
    return (data=zip(x, (view(y, :, j) for y ∈ ys for j ∈ axes(y, 2))...),
            hasmissing=hasmissing, nobj=nobj)
end

function filldata(x::AVR, y1::Union{Real,AVR}, y2::Union{Real,AVR})
    isempty(x) && throw(ArgumentError("Cannot display empty vectors."))
    y1 isa AV || (y1 = fill(y1, length(x)))
    y2 isa AV || (y2 = fill(y2, length(x)))
    length(x) == length(y1) == length(y2) ||throw(DimensionMismatch("vectors must have " *
                                                                    "matching lengths"))
    return (data=zip(x, y1, y2))
end

function histdata(x::AV{<:CanMiss{<:Real}})
    isempty(x) && throw(ArgumentError("Cannot display an empty vector."))
    sx = skipmissing(x)
    return (data=zip(x),
            hasmissing=(Missing <: eltype(x)),
            nobs=sum(e->1, sx),
            range=fl(extrema(sx)))
end

#################################################################################
#################################################################################

####
#### plot, plot!
####

"""
    plot!(...)

Add a plot. Keyword arguments can be passed to specify the linestyle(s), label(s) and
markerstyle(s).

## Example

```julia
x = range(-2, 2, length=100)
y = @. exp(-abs(x)+sin(x))
plot(x, y, color="blue", lstyle="--", marker="o", lwidth=0.05, label="First plot")
```
"""
function plot!(x, ys...; axes=nothing, overwrite=false, o...)
    axes = check_axes(axes)
    if overwrite
        # are current axes empty? if so don't do anything as the user may have pre-specified
        # things like xlim etc and want the plot to appear with those
        # if it's not empty, reset the axes (will destroy xlim settings etc as it should)
        all(isempty, (axes.drawings, axes.objects)) || reset!(axes; parent=axes.parent)
    end
    pd = plotdata(x, ys...)
    scatter = Scatter2D(pd.data, pd.hasmissing, pd.nobj)
    set_properties!(scatter; defer_preview=true, o...)
    push!(axes.drawings, scatter)
    return DrawingHandle(scatter)
end

function plot!(f::Function, from, to; length=100, o...)
    x = range(from, stop=to, length=length)
    plot!(x, f.(x); o...)
end

"""
    plot(...)

Erase previous drawings and add a plot. See also [`plot!`](@ref).
"""
plot(a...; o...) = plot!(a...; overwrite=true, o...)

"""
    scatter!(...)

Add a scatter plot (no line joins the points). See also [`plot!`](@ref).
"""
scatter!(a...; o...) = plot!(a...; ls="none", marker=".", o...)

"""
    scatter(...)

Erase previous drawings and add a scatter plot (no line joins the points). See also [`plot`](@ref).
"""
scatter(a...; o...)  = scatter!(a...; overwrite=true, o...)

####
#### fill_between!, fill_between
####

"""
    fill_between!(...)

Add a fill plot between two lines. The arguments must not have missings but `y1` and/or `y2` can
be specified as single numbers (= horizontal line).
"""
function fill_between!(x, y1, y2; axes=nothing, overwrite=false, o...)
    axes = check_axes(axes)
    if overwrite
        all(isempty, (axes.drawings, axes.objects)) || reset!(axes; parent=axes.parent)
    end
    fill = Fill2D(data=filldata(x, y1, y2))
    set_properties!(fill; defer_preview=true, o...)
    push!(axes.drawings, fill)
    return DrawingHandle(fill)
end

"""
    fill_between(...)

Erase previous drawings and add a fill plot. See also [`fill_between!`](@ref).
"""
fill_between(a...; o...) = fill_between!(a...; overwrite=true, o...)

####
#### hist, hist!
####

"""
    hist!(...)

Add a histogram of `x` on the current axes.
"""
function hist!(x; axes=nothing, overwrite=false, o...)
    axes = check_axes(axes)
    if overwrite
        all(isempty, (axes.drawings, axes.objects)) || reset!(axes; parent=axes.parent)
    end
    hd = histdata(x)
    hist = Hist2D(data=hd.data, hasmissing=hd.hasmissing, nobs=hd.nobs, range=hd.range)
    set_properties!(hist; defer_preview=true, o...)
    push!(axes.drawings, hist)
    return DrawingHandle(hist)
end

"""
    hist(...)

Erase previous drawings and add a histogram. See also [`hist!`](@ref).
"""
hist(a...; o...)  = hist!(a...; overwrite=true, o...)

####
#### bar!, bar
####

"""
    bar!(...)

Add a bar plot.
"""
function bar!(x, ys...; axes=nothing, overwrite=false, o...)
    axes = check_axes(axes)
    if overwrite
        all(isempty, (axes.drawings, axes.objects)) || reset!(axes; parent=axes.parent)
    end
    bd = plotdata(x, ys...)
    bar = Bar2D(bd.data, bd.hasmissing, bd.nobj)
    set_properties!(bar; defer_preview=true, o...)
    push!(axes.drawings, bar)
    return DrawingHandle(bar)
end

"""
    bar(...)

Erase previous drawings and add a bar plot. See also [`bar!`](@ref).
"""
bar(a...; o...) =  bar!(a...; overwrite=true, o...)


text: divisiblebythree(str) = mapreduce( c -> c-'0', (a,b) -> (a+b)%3, str, init=0 ) == 0

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
@testset "SmoothingSpline" begin

    @testset "SmoothingSpline" begin

        # EURAAA_20191111
        @testset "EURAAA_20191111" begin
            euraaa_yields = [-0.602, -0.6059, -0.6096, -0.613, -0.6215, -0.6279, -0.6341, -0.6327, -0.6106, -0.5694, -0.5161, -0.456, -0.3932, -0.3305, -0.2698, -0.2123, -0.1091, 0.0159, 0.0818, 0.1601, 0.2524, 0.315]
            euraaa_maturities = [0.25, 0.333, 0.417, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 15, 17, 20, 25, 30]
            ss_param = SmoothingSpline(0.6440719852424944, -1.2135915867261, -1.8281745438894712, 0.1)
            ss = est_ss_params(euraaa_yields, euraaa_maturities)
            @test ss = ss_param
        end
    end

end


text:
"""
    tags(logger)

Returns a set of all the strings used as tags in messages serialized by
`logger`.

`logger` can be a `TBLogger` or the path of a valid TensorBoard logdir.
"""
function tags(tbl; kwargs...)
    tags = Set{String}()

    map_summaries(tbl; kwargs...) do tag, iter, val
        push!(tags, tag)
    end

    return tags
end

"""
    steps(logger)

Returns a set of all the steps used as tags in messages serialized by
`logger`.

`logger` can be a `TBLogger` or the path of a valid TensorBoard logdir.
"""
function steps(tbl; kwargs...)
    steps = BitSet()

    map_events(tbl; kwargs...) do ev
        push!(steps, ev.step)
    end

    return steps
end


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

using Markdown
using InteractiveUtils

# ╔═╡ bb6a9f88-2134-11ec-3337-2bdd2c6747a9
begin
	using PlutoUI
end

# ╔═╡ 7a46f79a-c4f4-4cea-8187-839337a9c962
PlutoUI.LocalResource("./figs/Rankine.png")

# ╔═╡ 88b9e9f5-f1ed-4146-a63d-f941be18c559
md"""

#### Example: Compute the steady-state rate of shaft work $\dot{W}_{sh}$ (and some other stuff)

Compute the steady-state rate of shaft work $\dot{W}_{sh}$ (units: kW), and the irreversible enthlapy change from an insulated turbine with efficiency of $\eta$ =86\%. The working fluid is water with a steady-state flow rate of $\dot{m}$ = 4.0 kg/s. Let the input temperature
T$_{2}$ = 500$^{\circ}$C and pressure P$_{2}$ = 11 MPa. The temperature of the exit stream from the turbine T$_{3}$ = 100$^{\circ}$C.

"""

# ╔═╡ 5cf4b776-e31a-48c3-be81-1d7517e9db8b
md"""
#### Approach

The steady-state open first law (energy balance) relates the heat, work and enthalpy for the system:

$$\dot{Q}+\dot{W}^{*}_{sh} + \sum_{s}^{\mathcal{S}}v_{s}\dot{m}_{s}H_{s} = 0$$

where $\dot{W}^{*}_{sh}$ denotes the rate of _actual_ shaft work recovered from the turbine. If the turbine is operating _reversibly_ then the rate of shaft work that we actually recover is equal to the theoretical maximum shaft work. However, in this problem, we know that the turbine is not operating reversibly, thus we use $\dot{W}^{*}_{sh}$.

Let's begin by throwing out terms in the open energy balance, and solving for the rate of shaft work $\dot{W}^{*}_{sh}$. The problem says the turbine is insulated, thus $\dot{Q}$ = 0. Second, we know that $\mathcal{S}$ = 2, and that the mass flow rate into the turbine equals mass flow rate out teh turbine $\dot{m}_{2}=\dot{m}_{3}\equiv{\dot{m}}$. Solving for the rate of shaft work gives the expression:

$$\dot{W}^{*}_{sh} = \dot{m}\left(H^{\prime}_{3} - H_{2}\right)_{I}$$

where the $\left(\star\right)_{I}$ term denotes the enthalpy difference for the _actual irreversible_ process, and $H^{\prime}_{3}$ denotes the _actual_ operating point 3. However, we do not know what the _actual irreversible_ enthalpy difference is i.e., we don't know $H^{\prime}_{3}$.

However, we do know:

* the turbine efficiency $\eta_{T}$ and
* that a reversible adiabatic process is _constant entropy_ ($S_{2}=S_{3}$).

##### How can we compute $\left(\star\right)_{I}$?

The turbine efficiency is given by:

$$\eta_{T} = \frac{\Delta{H}_{I}}{\Delta{H}_{R}}$$

where $\Delta{H}_{R}$ denotes the enthalpy difference for a _reversible_ process. Using the efficiency expression gives us a method to compute
$\Delta{H}_{I}$, the _irreversible_ enthalpy difference:

$$\Delta{H}_{I} = \eta_{T}\Delta{H}_{R}$$

which we can substitute into the expression for work to give:

$$\dot{W}^{*}_{sh} = \dot{m}\eta_{T}\left(H_{3} - H_{2}\right)_{R}$$

__Why__? Because we can follow along the line of constant entropy on a PH-diagram until we hit the outlet pressure (or temperature), which gives $H_{3}$.

"""

# ╔═╡ 75376d0d-c4eb-4bc5-a58b-752187aae593
md"""

##### Estimate $\Delta{H}_{R}$ from the Pressure-Enthalpy (PH) diagram

"""

# ╔═╡ c0d0a66d-4651-4317-9e08-c16ffa290a68
PlutoUI.LocalResource("./figs/PH-diagram-Water-OPs.png")

# ╔═╡ e449301c-c198-4506-9512-02da59cac3f1
md"""
##### Estimate $\Delta{H}_{I}$ and visualize enthalpy difference on a PH-diagram
"""

# ╔═╡ b7ddbd1a-de34-4cd1-ae89-e0c70f3e6c7e
begin

	# Data from the problem and the PH-diagram -
	mdot = 4.0 	# units: kg/s
	H₂ = 3400  	# units: kJ/kg
	H₃= 2400  	# units: kJ/kg
	ηₜ = 0.86;	# units: dimensionless

	# return -
	nothing
end

# ╔═╡ 571973ba-b1a3-4698-a02a-c8f84337d512
# What is the work recovered from the turbine?
Wdot = ηₜ*(mdot*(H₃ - H₂))

# ╔═╡ 5f4dff5a-5a4c-4948-82df-161e27456578
PlutoUI.LocalResource("./figs/PH-diagram-Water-OPs-IRR.png")

# ╔═╡ db4037dd-d4a3-42a3-a357-2540972294bd
# What is the irreversible enthlapy change?
ΔH_I = ηₜ*(H₃-H₂)

# ╔═╡ c285ccd0-bf5c-423f-9577-361ff21f29ad
md"""
##### What is the rate of entropy generation $\dot{S}_{G}$?

The open second law (entropy balance) at steady-state (with a single heat-transfer surface) is given by:

$$-\left(\sum_{s}^{\mathcal{S}}v_{s}\dot{m}_{s}S_{s} + \frac{\dot{Q}}{T_{\sigma}}\right) = \dot{S}_{G}\geq{0}$$

The turbine is adiabatic ($\dot{Q}$ = 0), is operating at steady-state, and has a single input and output ( $\dot{m}_{2}=\dot{m}_{3}\equiv{\dot{m}}$). Putting these ideas together, and throwing out terms in the open entropy balance that we don't need, gives an expression for the rate of entropy generation:

$$-\dot{m}\left(S_{2} - S_{3}\right) = \dot{S}_{G}$$

"""

# ╔═╡ 988ab1b3-a5e5-40df-abc5-c36c9d6cf9b9
begin
	# lookup entropy values from the PH diagram -
	S₂ = 6.5 	# units: kJ/kg-K
	S₃ = 7.0 	# units: kJ/kg-K

	# return -
	nothing
end

# ╔═╡ c787d47d-7a06-4603-8a12-cdaec354374e
# compute SG -
SG_dot = -1*mdot*(S₂-S₃)

# ╔═╡ bb9d6c95-1391-489c-8c00-d5f27ac832d2
md"""
##### What is the rate of lost work $\dot{W}_{lost}$?
From the drivation in the lectire notes, we know that the rate of lost work $\dot{W}_{lost}$ is given by:

$$\dot{W}_{lost} = \left(\eta_{T} - 1\right)\dot{W}_{sh}$$

However, $\dot{W}_{sh}$ is given by:

$$\dot{W}_{sh} = \dot{m}\Delta{H}_{R}$$

Substituting $\dot{W}_{sh}$ into the lost work expression gives:

$$\dot{W}_{lost} = \left(\eta_{T} - 1\right) \dot{m}\Delta{H}_{R}$$

"""

# ╔═╡ db18c3bf-15b3-47c9-be0b-1f543e07fc19
# compute the reversible (ideal) enthalpy change -
ΔH_R = H₃ - H₂

# ╔═╡ 3fc1d446-16dd-479e-bc80-03df6804f557
W_dot_lost = (ηₜ - 1)*ΔH_R

# ╔═╡ c6b5e534-d843-4637-aeb8-c26528605f7c
md"""
#### What is the surrounding temperature $T_{\sigma}$?

We know from the class notes that:

$$\dot{W}_{lost} = T_{\sigma}\dot{S}_{G}$$

thus we can solve for $T_{\sigma}$:

$$T_{\sigma} = \frac{\dot{W}_{lost}}{\dot{S}_{G}}$$

"""

# ╔═╡ e1738875-0394-486a-94c2-0a20647a6857
# Compute T_σ -
T_σ = W_dot_lost/SG_dot

# ╔═╡ e4681733-6926-44a9-9eea-ace0f74b797d
html"""<style>
main {
    max-width: 1200px;
  	width: 80%;
	margin: auto;
	font-family: "Roboto, monospace";
}

a {
	color: blue;
	text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.12"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "f35ae11e070dbf123d5a6f54cbda45818d765ad2"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.12"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═bb6a9f88-2134-11ec-3337-2bdd2c6747a9
# ╟─7a46f79a-c4f4-4cea-8187-839337a9c962
# ╟─88b9e9f5-f1ed-4146-a63d-f941be18c559
# ╟─5cf4b776-e31a-48c3-be81-1d7517e9db8b
# ╟─75376d0d-c4eb-4bc5-a58b-752187aae593
# ╟─c0d0a66d-4651-4317-9e08-c16ffa290a68
# ╟─e449301c-c198-4506-9512-02da59cac3f1
# ╠═b7ddbd1a-de34-4cd1-ae89-e0c70f3e6c7e
# ╠═571973ba-b1a3-4698-a02a-c8f84337d512
# ╟─5f4dff5a-5a4c-4948-82df-161e27456578
# ╠═db4037dd-d4a3-42a3-a357-2540972294bd
# ╟─c285ccd0-bf5c-423f-9577-361ff21f29ad
# ╠═988ab1b3-a5e5-40df-abc5-c36c9d6cf9b9
# ╠═c787d47d-7a06-4603-8a12-cdaec354374e
# ╟─bb9d6c95-1391-489c-8c00-d5f27ac832d2
# ╠═db18c3bf-15b3-47c9-be0b-1f543e07fc19
# ╠═3fc1d446-16dd-479e-bc80-03df6804f557
# ╟─c6b5e534-d843-4637-aeb8-c26528605f7c
# ╠═e1738875-0394-486a-94c2-0a20647a6857
# ╟─e4681733-6926-44a9-9eea-ace0f74b797d
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module CommonEq

import Base.Order.Lt  # Base.Order.Lt(lt::T) is signature

export  Lt, ≪, Le, ≦, Eq, ⩵, Ne, ≶, ≷, Ge, ≫, Gt, ≧

"Lt"
Lt(a::Any, b::Any) = nothing

"Le"
Le(a::Any, b::Any) = nothing

"Eq"
Eq(a::Any, b::Any) = nothing

"Ne"
Ne(a::Any, b::Any) = nothing

"Ge"
Ge(a::Any, b::Any) = nothing

"Gt"
Gt(a::Any, b::Any) = nothing

"≪, `\\ll[tab]` or `Lt`"
const ≪  = Lt

"≦, `\\leqq[tab]` or `Le`"
const ≦ = Le

"⩵, `\\Equal[tab]` or `Eq`"
const ⩵ = Eq

"≶ `\\lessgtr[tab]` or `Ne`"
const ≶ = Ne

"≷ `\\gtrless[tab]` or `Ne`"
const ≷ = Ne

"≧, \\geqq[tab]` or `Ge`"
const ≧ = Ge

"≫, \\gg[tab] or `Gt`"
const ≫ = Gt

end


text:
mutable struct CachedBoolVariate <: AbstractBoolVariate
    cache::Union{BitArray{1}, Nothing}
    basevariate::AbstractBoolVariate
    lockobj::Threads.TatasLock
end

Base.length(var::CachedBoolVariate) = length(var.basevariate)

getname(var::CachedBoolVariate) = getname(var.basevariate)

function CachedBoolVariate(basevariate::AbstractBoolVariate)
    CachedBoolVariate(nothing, basevariate, Threads.TatasLock())
end

function cache(basevariate::AbstractBoolVariate)
    CachedBoolVariate(nothing, basevariate, Threads.TatasLock())
end

function slice(boolvariate::CachedBoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    basevar = boolvariate.basevariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if isa(basevar, CachedBoolVariate) || isa(basevar, BoolVariate)
        slice(basevar, fromobs, toobs, slicelength)
    else
        lockobj = boolvariate.lockobj
        lock(lockobj)
        try
            if boolvariate.cache === nothing
                slices = slice(basevar, 1, length(basevar), SLICELENGTH)
                cachedata = BitArray{1}(undef, length(basevar))
                fold(0, slices) do offset, slice
                    n = length(slice)
                    view(cachedata, (1 + offset):(n + offset)) .= slice
                    offset += n
                end
                boolvariate.cache = cachedata
            end
        finally
            unlock(lockobj)
        end
        slice(boolvariate.cache, fromobs, toobs, slicelength)
    end
end

function Base.map(boolvariate::CachedBoolVariate, dataframe::AbstractDataFrame)
    CachedBoolVariate(map(boolvariate.basevariate, dataframe))
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
export Transfer, transfer#, colstop, resizedata!

@compat abstract type AbstractTransfer{T} <: Operator{T} end

ApproxFunBase.israggedbelow(L::AbstractTransfer) = true
#Base.issymmetric(L::AbstractTransfer) = false

#for OP in (:(ApproxFun.domain),:rangedomain)
#    @eval $OP(L::AbstractTransfer) = $OP(markovmap(L))
#end

@compat struct ConcreteTransfer{T,D<:Space,R<:Space,M<:AbstractMarkovMap} <: AbstractTransfer{T}
  m::M
  domainspace::D #Domain space
  rangespace::R
  colstops::Array{Int,1}
  #    f::Array{Vector{T},1}
  # function ConcreteTransfer{T,D<:Space,R<:Space,M<:AbstractMarkovMap}(m::AbstractMarkovMap,domainspace::Space,rangespace::Space,colstops::Array{Int,1}=Int[])
  #   @assert domain(m)==domain(domainspace)
  #   @assert rangedomain(m)==domain(rangespace)
  #   nneutral(m) != 0 && error("Neutral fixed points not supported for transfer operator")
  #   new{T,typeof(domainspace),typeof(rangespace),typeof(m)}(m,domainspace,rangespace,colstops)
  # end
end

"""
    Transfer(m::AbstractMarkovMap)

Create a `CachedOperator` of a `ConcreteTransfer` operator encoding the transfer operator of `m`.

Caching is used for speed, as entries of the transfer operator are most efficiently calculated whole columns at a time.
"""
Transfer(stuff...;padding=false) = cache(ConcreteTransfer(stuff...),padding=padding)

# ConcreteTransfer(::Type{T},m::AbstractMarkovMap,dom::Space=Space(domain(m)),
#                     ran::Space=(domain(m)==rangedomain(m) ? dom : Space(rangedomain(m))),
#                     colstops::Array{Int,1}=Int[]) where T =
#   ConcreteTransfer{T,typeof(dom),typeof(ran),typeof(m)}(m,dom,ran,colstops)#,domainspace(m),rangespace(m))
ConcreteTransfer(m,dom::Space=Space(domain(m)),ran=(domain(m)==rangedomain(m) ? dom : Space(rangedomain(m))),colstops::Array{Int,1}=Int[]) =
  ConcreteTransfer{prectype(rangedomain(m)),typeof(dom),typeof(ran),typeof(m)}(m,dom,ran,colstops)

for OP in (:domainspace,:rangespace)
  @eval ApproxFun.$OP(L::ConcreteTransfer) = L.$OP
end

function Base.convert(::Type{Operator{T}},D::ConcreteTransfer) where T
    if T==eltype(D)
        D
    else
        ConcreteTransfer{T,typeof(D.domainspace),typeof(D.rangespace),typeof(D.m)}(D.m,D.domainspace,D.rangespace,D.colstops)
    end
end

markovmap(L::ConcreteTransfer) = L.m
markovmap(L::CachedOperator) = markovmap(L.op)

function resizecolstops!(L::ConcreteTransfer,n)
  colstopslength = length(L.colstops)
  if n > colstopslength
    append!(L.colstops,-ones(Int,n-colstopslength)) #-1 is clunky but w/e
  end
end

function setcolstops!(L::ConcreteTransfer,inds,vals) #should be generic range/collection I guess?
  resizecolstops!(L,maximum(inds))
  L.colstops[inds] = vals
end

# Transfer a fun - TODO

# TODO: structs for transferfunction
struct TransferFunction{M,F} <: Function
    m::M
    fn::F
end
(tf::TransferFunction)(x) = transferfunction(x,tf.m,tf.fn)
function transfer(m::AbstractIntervalMap,fn, sp=Space(rangedomain(m)))
  Fun(TransferFunction(m,fn),sp)
end
transfer(L::AbstractTransfer,fn) = transfer(markovmap(L),fn,rangespace(L))
# transfer(m::AbstractIntervalMap,fn,x) = transferfunction(x,m,fn)

# Indexing

transferfunction_nodes(L::ConcreteTransfer{TT,D,R,M},n::Integer,kk) where {TT,D,R,M<:AbstractIntervalMap} =
  prectype(domainspace(L))[transferfunction(p,markovmap(L),BasisFun(domainspace(L),kk)) for p in points(rangespace(L),n)]
# transferfunction_nodes{TT,D,R,M<:MarkovInverseCache}(L::ConcreteTransfer{TT,D,R,M},n::Integer,kk) =
#   T[transferfunction(InterpolationNode(rangespace(markovmap(L)),k,n),markovmap(L),BasisFun(domainspace(L),kk)) for k = 1:n]

function transfer_getindex(L::ConcreteTransfer{T},
    jdat::Tuple{Integer,Integer,Union{Integer,Infinity}},
    k::OrdinalRange,padding::Bool=false) where T
  Tr = real(T)
  @compat dat = Array{T}(undef,0)
  @compat cols = Array{eltype(k)}(undef,Base.length(k)+1)
  cols[1] = 1
  K = isfinite(jdat[3]) ? length(jdat[1]:jdat[2]:jdat[3]) : 0 # largest colstop
  padding && start(k) > 1 && (K = max(maximum([colstop(L,i) for i =1:start(k)-1]),K))
  rs = rangespace(L)

  resizecolstops!(L,maximum(k))
  mc = maximum(L.colstops[1:first(k)])

  for (kind,kk) in enumerate(k)
    kind > 1 && (mc = max(mc,maximum(L.colstops[k[kind-1]+1:kk])))

    tol = (Tr==Any) ? 200eps() : 200eps(Tr)

    if L.colstops[kk] >= 1
      @compat coeffs = ApproxFun.transform(rs,transferfunction_nodes(L,max(16,nextpow(2,L.colstops[kk])),kk))[1:L.colstops[kk]]
      maxabsc = max(maximum(abs.(coeffs)),one(Tr))
      chop!(coeffs,tol*maxabsc*log2(length(coeffs)))
    elseif L.colstops[kk]  == 0
      coeffs = zeros(T,1)
    else

      #       if mc ≤ 2^4
      #         coeffs = Fun(f,rs).coefficients
      #       else

      # code from ApproxFun v0.5 (src/Fun/constructors.jl) under BSD license (same as this package)
      # ( the difference is we start at n \approx mc )

      r=ApproxFun.checkpoints(rs)
      fr=[transferfunction(rr,markovmap(L),BasisFun(domainspace(L),kk)) for rr in r]
      maxabsfr=norm(fr,Inf)

      logn = min(round(Int,log2(max(mc,16)),RoundUp),20)
      while logn < 21
        coeffs = ApproxFun.transform(rs,transferfunction_nodes(L,2^logn,kk))

        maxabsc = max(one(Tr),maximum(abs.(coeffs)))
        if maxabsc == 0 && maxabsfr == 0
          coeffs = zeros(T,1)
          break
        else
          maxabsfr = max(maxabsfr,one(Tr))
          b = ApproxFun.block(rs,length(coeffs))
          bs = ApproxFun.blockstart(rs,max(div(2Integer(b),3),1))
          if length(coeffs) > 8 && maximum(abs.(coeffs[bs:end])) < 20tol*maxabsc*logn &&
              all(kkk->norm(Fun(rs,coeffs)(r[kkk])-fr[kkk],1)<tol*length(coeffs)*maxabsfr*500,1:length(r))
            chop!(coeffs,tol*maxabsc*logn)
            break
          end
        end
        logn += 1
      end

      if logn == 21
        warn("Maximum number of coefficients "*string(2^20+1)*" reached in constructing $(kk)th column.")
        coeffs = ApproxFun.transform(rs,transferfunction_nodes(L,2^logn,kk))
      end

    end

    #    tol = 4ceil(log2(length(tk.coefficients)))*eps(T) # 4 is dummy should be like 1 in Fourier, 2C+1 in Cheby
    #    chop!(tk,tol)
    lcfc = length(coeffs)
    maxabsc = lcfc > 0 ? maximum(abs.(coeffs)) : zero(eltype(coeffs))

    for i = 1:lcfc # set small enough coefficients (in middle of vector) to zero
      abs(coeffs[i])<tol*maxabsc*log2(lcfc)/10 && (coeffs[i] = 0)
    end

    cutcfc = coeffs[(jdat[1]:jdat[2]:min(lcfc,jdat[3]))::OrdinalRange]

    K = max(K,length(cutcfc))
    padding && ApproxFun.pad!(cutcfc,K)
    append!(dat,cutcfc)
    cols[kind+1] = cols[kind]+length(cutcfc)

    setcolstops!(L,kk,lcfc)
    #   kk == 1 && display(stacktrace())
  end
  RaggedMatrix(dat,cols,K)
end

Base.getindex(L::ConcreteTransfer,j::OrdinalRange,k::OrdinalRange) = transfer_getindex(L,(start(j),step(j),last(j)),k)
Base.getindex(L::ConcreteTransfer,j::Colon,k::OrdinalRange) = transfer_getindex(L,(1,1,InfiniteArrays.∞),k)
Base.getindex(L::ConcreteTransfer,j::Integer,k::Integer) = Base.getindex(L,j:j,k:k)[1,1]
Base.getindex(L::ConcreteTransfer,j::Integer,k::OrdinalRange) = Base.getindex(L,j:j,k).data
Base.getindex(L::ConcreteTransfer,j::OrdinalRange,k::Integer) = Base.getindex(L,j,k:k).data

Base.convert(::Type{RaggedMatrix},S::ApproxFun.SubOperator{T,LL,Tuple{R1,R2}}) where
        {T,LL<:AbstractTransfer,R1<:OrdinalRange,R2<:OrdinalRange} =
  Base.getindex(parent(S),parentindexes(S)[1],parentindexes(S)[2])

function BandedMatrices.colstop(L::ConcreteTransfer,k::Integer)
  k <= length(L.colstops) && L.colstops[k] != -1 && return L.colstops[k]
  transfer_getindex(L,(1,1,Infinity),k:k)
  L.colstops[k]
end

function ApproxFun.resizedata!(co::CachedOperator{T,RaggedMatrix{T},CT},
    ::Colon,n::Integer) where {T<:Number,CT<:ConcreteTransfer}
  if n > co.datasize[2]
    RO = transfer_getindex(co.op,(1,1,InfiniteArrays.∞),(co.datasize[2]+1):n,co.padding)
    if co.datasize[2] == 0
      co.data = RO
      co.datasize = (co.data.m,n)
    else
      append!(co.data.data,RO.data)
      append!(co.data.cols,RO.cols[2:end].+(co.data.cols[end]-1))
      co.data.m = max(co.data.m,RO.m)
      co.datasize = (co.data.m,n)
    end
  end

  co
end

function BandedMatrices.colstop(co::CachedOperator{T,DM,CT},n::Integer) where
    {T<:Number,DM<:AbstractMatrix,CT<:ConcreteTransfer}
  ApproxFun.resizedata!(co,:,n)
  BandedMatrices.colstop(co.data,n)
end

# # fast colstop for cached operators
# function BandedMatrices.colstop{T,M<:ConcreteTransfer,DS<:Space,RS<:Space}(B::CachedOperator{T,RaggedMatrix{T},M,DS,RS,Tuple{Infinity},k::Integer)
#   ApproxFun.resizedata!(B,:,k)
#   BandedMatrices.colstop(co.op,k)
# end


text:
@doc raw"""
Symmetric quadratic equilibrium in (x,y,z) coordinates with covariant components of
the vector potential given by
```math
A (x,y,z) = \frac{B_0}{4} \big( - y \, (2 + x^2 + y^2) , \, x \, (2 + x^2 + y^2) , \, 0 \big)^T
```
resulting in the magnetic field with covariant components
```math
B(x,y,z) = B_0 \, \begin{pmatrix}
0 \\
0 \\
1 + x^2 + y^2 \\
\end{pmatrix}
```

Parameters: `B₀`
"""
struct SymmetricQuadratic{T <: Number} <: CartesianEquilibrium
    name::String
    B₀::T
    SymmetricQuadratic{T}(B₀::T) where T <: Number = new("SymmetricQuadraticEquilibrium", B₀)
end

SymmetricQuadratic(B₀::T) where T <: Number = SymmetricQuadratic{T}(B₀)

function Base.show(io::IO, equ::SymmetricQuadratic)
    print(io, "Symmetric Quadratic Equilibrium")
end

function R(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    r(x,equ)
end

function r(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    sqrt(r²(x,equ))
end

function r²(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    X(x,equ)^2 + Y(x,equ)^2
end

function θ(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    atan2(Y(x,equ), X(x,equ))
end

function ϕ(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    θ(x,equ)
end

function A₁(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    - equ.B₀ * x[2] * (2 + x[1]^2 + x[2]^2) / 4
end

function A₂(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    + equ.B₀ * x[1] * (2 + x[1]^2 + x[2]^2) / 4
end

function A₃(x::AbstractArray{T,1}, equ::SymmetricQuadratic) where {T <: Number}
    zero(eltype(x))
end

macro symmetric_quadratic_equilibrium()
    generate_equilibrium_code(SymmetricQuadratic(1.); output=false)
end


text:
using LinearAlgebra, SphericalHarmonicExpansions, Test

include("sphericalHarmonics.jl")
include("sphericalHarmonicsExpansion.jl")
include("sphericalQuadrature.jl")
include("fastfunc.jl")
include("translation.jl")
include("errorPropagation.jl")


text:
#=
function evolve(prob, ctrl_amps, timeslices)
    Δt = prob.transfer_time/prob.n_steps
    U = Matrix(I, size(prob.drift_Ham))
    n_pulses = length(prob.ctrl_Hams)

    for i = 1:timeslices
        H_tot = prob.drift_Ham

        for k = 1:n_pulses
            H_tot = H_tot + ctrl_amps[i, k]*prob.ctrl_Hams[k]
        end
        U = exp(-1.0im * Δt * H_tot) * U
    end
    U
end
=#

function evolve(drift_Ham, ctrl_Hams, n_steps, transfer_time, ctrl_amps, timeslices)
    Δt = transfer_time/n_steps
    U = Matrix(I, size(drift_Ham))
    n_pulses = length(ctrl_Hams)

    for i = 1:timeslices
        H_tot = drift_Ham

        for k = 1:n_pulses
            H_tot = H_tot + ctrl_amps[i, k]*ctrl_Hams[k]
        end
        U = exp(-1.0im * Δt * H_tot) * U
    end
    U
end

function evolve_back(prob, ctrl_amps, timeslices)
    U = Matrix(I, size(prob.drift_Ham))
    n_pulses = length(prob.ctrl_Hams)

    for timeslice in prob.transfer_time:-1:timeslices+1
        U = evolve(prob, ctrl_amps, timeslices)' * U
    end
    U
end

function density_op(prob, ctrl_amps, timeslices)
    evo = evolve(prob, ctrl_amps, timeslices)
    return evo*prob.init_density_op*evo'
end

function back_target_op(prob, ctrl_amps, timeslices)
    back_evo = evolve_back(prob, ctrl_amps, timeslices)
    return back_evo*prob.target_density_op*back_evo'
end


text:
##### deal with nutrients uptake
function gpu_calc_consume_kernel!(ctsdic, plank, ac, x, y, z, ΔT)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:size(ac,1)
        @inbounds CUDA.@atomic ctsdic[x[i], y[i], z[i]] += (plank.resp[i] - plank.PS[i]) * ΔT * ac[i]
    end
    return nothing
end
function calc_consume!(ctsdic, plank, ac, x, y, z, ΔT, ::GPU)
    @cuda threads=256 blocks=ceil(Int, size(ac,1)/256) gpu_calc_consume_kernel!(ctsdic, plank, ac, x, y, z, ΔT)
    return nothing
end
function calc_consume!(ctsdic, plank, ac, x, y, z, ΔT, ::CPU)
    for i in 1:size(ac,1)
        @inbounds ctsdic[x[i], y[i], z[i]] += (plank.resp[i] - plank.PS[i]) * ΔT * ac[i]
    end
    return nothing
end

##### deal with grazed or dead individuals
function gpu_calc_loss_kernel!(ctsdoc, ctspoc, plank, ac, x, y, z, loss, lossFracC)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:size(ac,1)
        @inbounds CUDA.@atomic ctsdoc[x[i], y[i], z[i]] += plank.Bm[i] * lossFracC * ac[i] * loss[i]
        @inbounds CUDA.@atomic ctspoc[x[i], y[i], z[i]] += plank.Bm[i] * (1.0-lossFracC) * ac[i] * loss[i]
    end
    return nothing
end
function calc_loss!(ctsdoc, ctspoc, plank, ac, x, y, z,
                    loss, lossFracC, ::GPU)
    @cuda threads=256 blocks=ceil(Int, size(ac,1)/256) gpu_calc_loss_kernel!(ctsdoc, ctspoc, plank, ac, x, y, z, loss, lossFracC)
    return nothing
end
function calc_loss!(ctsdoc, ctspoc, plank, ac, x, y, z, loss, lossFracC, ::CPU)
    for i in 1:size(ac,1)
        @inbounds ctsdoc[x[i], y[i], z[i]] += plank.Bm[i] * lossFracC * ac[i] * loss[i]
        @inbounds ctspoc[x[i], y[i], z[i]] += plank.Bm[i] * (1.0-lossFracC) * ac[i] * loss[i]
    end
    return nothing
end