text:
abstract type Camera{T} end
const Q = Quaternions # save some writing!

mutable struct OrthographicCamera{T} <: Camera{T}
    window_size     ::Signal{SimpleRectangle{Int}}
    view            ::Signal{Mat4{T}}
    projection      ::Signal{Mat4{T}}
    projectionview  ::Signal{Mat4{T}}
end

mutable struct PerspectiveCamera{T} <: Camera{T}
    window_size     ::Signal{SimpleRectangle{Int}}
    nearclip        ::Signal{T}
    farclip         ::Signal{T}
    fov             ::Signal{T}
    view            ::Signal{Mat4{T}}
    projection      ::Signal{Mat4{T}}
    projectionview  ::Signal{Mat4{T}}
    eyeposition     ::Signal{Vec{3, T}}
    lookat          ::Signal{Vec{3, T}}
    up              ::Signal{Vec{3, T}}
    trans           ::Signal{Vec{3, T}}
    theta           ::Signal{Vec{3, T}}
    projectiontype  ::Signal{Projection}
end

mutable struct DummyCamera{T, IT} <: Camera{T}
    window_size     ::Signal{SimpleRectangle{IT}}
    view            ::Signal{Mat4{T}}
    projection      ::Signal{Mat4{T}}
    projectionview  ::Signal{Mat4{T}}
end

function DummyCamera(;
        window_size    = Signal(SimpleRectangle(-1, -1, 1, 1)),
        view           = Signal(eye(Mat4f0)),
        nearclip       = Signal(-10_000f0),
        farclip        = Signal(10_000f0),
        projection     = const_lift(orthographicprojection, window_size, nearclip, farclip),
        projectionview = const_lift(*, projection, view)
    )
    DummyCamera(window_size, view, projection, projectionview)
end

function Base.collect(camera::Camera, collected = Dict{Symbol, Any}())
    names = fieldnames(camera)
    for name in (:view, :projection, :projectionview, :eyeposition)
        if name in names
            collected[name] = getfield(camera, name)
        end
    end
    return collected
end

"""
Creates an orthographic camera with the pixel perfect plane in z == 0
Signals needed:
Dict(
    :window_size           => Signal(SimpleRectangle{Int}),
    :buttons_pressed       => Signal(Int[]),
    :mouse_buttons_pressed => Signal(Int[]),
    :mouseposition         => mouseposition, -> Panning
    :scroll_y              => Signal(0) -> Zoomig
)
"""
function OrthographicPixelCamera(
        inputs;
        fov=41f0, near=0.01f0, up=Vec3f0(0,1,0),
        translation_speed = Signal(1), theta = Signal(Vec3f0(0)), keep = Signal(true)
    )
    @materialize mouseposition, mouse_buttons_pressed, buttons_pressed, scroll = inputs
    left_ctrl     = Set([GLFW.KEY_LEFT_CONTROL])
    use_cam       = map(AND, const_lift(==, buttons_pressed, left_ctrl), keep)

    mouseposition = droprepeats(map(Vec2f0, mouseposition))
    left_pressed  = const_lift(pressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    xytranslate   = dragged_diff(mouseposition, left_pressed, use_cam)

    ztranslate    = filterwhen(use_cam, 0f0, map(x->Float32(x[2]), scroll))

    trans = map(translationlift, xytranslate, ztranslate, translation_speed)
    OrthographicPixelCamera(
        theta, trans, Signal(up), Signal(fov), Signal(near),
        inputs[:window_area],
    )
end
function OrthographicPixelCamera(
        theta, trans, up, fov_s, near_s, area_s
    )
    fov, near = Reactive.value(fov_s), Reactive.value(near_s)

    # lets calculate how we need to adjust the camera, so that it mapps to
    # the pixel of the window (area)
    area = Reactive.value(area_s)
    h = Float32(tan(fov / 360.0 * pi) * near)
    w_, h_ = area.w / 2f0, area.h / 2f0
    zoom = min(h_, w_) / h
    x, y = w_, h_
    eyeposition = Signal(Vec3f0(x, y, zoom))
    lookatvec   = Signal(Vec3f0(x, y, 0))
    far         = Signal(zoom * 10f0) # this should probably be not calculated
    # since there is no scene independant, well working far clip

    PerspectiveCamera(
        theta,
        trans,
        lookatvec,
        eyeposition,
        up,
        area_s,
        fov_s, # Field of View
        near_s,  # Min distance (clip distance)
        far, # Max distance (clip distance)
        Signal(GLAbstraction.ORTHOGRAPHIC)
    )
end

pressed(keys, key) = key in keys
singlepressed(keys, key) = length(keys) == 1 && first(keys) == key

mouse_dragg(v0, args) = mouse_dragg(v0..., args...)
function mouse_dragg(
        started::Bool, waspressed::Bool, startpoint, diff,
        ispressed::Bool, position, start_condition::Bool
    )
    started && ispressed && return (true, ispressed, startpoint, position-startpoint)
    if !started && !waspressed && ispressed && start_condition
        return (true, true, position, Vec2f0(0))
    end
    (false, ispressed, Vec2f0(0), Vec2f0(0))
end
mouse_dragg_diff(v0, args) = mouse_dragg_diff(v0..., args...)
function mouse_dragg_diff(
        started::Bool, waspressed, position0, diff,
        ispressed::Bool, position, start_condition::Bool
    )
    if !started && ispressed && (start_condition && !waspressed)
        return (true, ispressed, position, Vec2f0(0))
    end
    started && ispressed && return (true, ispressed, position, position0-position)
    (false, ispressed, Vec2f0(0), Vec2f0(0))
end

function dragged(mouseposition, key_pressed, start_condition = true)
    v0 = (false, false, Vec2f0(0), Vec2f0(0))
    args = const_lift(tuple, key_pressed, mouseposition, start_condition)
    dragg_sig = foldp(mouse_dragg, v0, args)
    is_dragg = map(first, dragg_sig)
    dragg = map(last, dragg_sig)
    dragg_diff = filterwhen(is_dragg, Reactive.value(dragg), dragg)
    dragg_diff
end
function dragged_diff(mouseposition, key_pressed, start_condition=true)
    v0 = (false, false, Vec2f0(0), Vec2f0(0))
    args = const_lift(tuple, key_pressed, mouseposition, start_condition)
    dragg_sig = foldp(mouse_dragg_diff, v0, args)
    is_dragg = map(first, dragg_sig)
    dragg_diff = map(last, dragg_sig)
    dragg_diff
end

"""
Transforms a mouse drag into a selection from drag start to drag end
"""
function drag2selectionrange(v0, selection)
    mousediff, id_start, current_id = selection
    if mousediff != Vec2f0(0) # Mouse Moved
        if current_id[1] == id_start[1]
            return min(id_start[2],current_id[2]):max(id_start[2],current_id[2])
        end
    else # if mouse did not move while dragging, make a single point selection
        if current_id.id == id_start.id
            return current_id.index:0 # this is the type stable way of indicating, that the selection is between currend_index
        end
    end
    v0
end

"""
Returns two signals, one boolean signal if clicked over `robj` and another
one that consists of the object clicked on and another argument indicating that it's the first click
"""
function clicked(robj::RenderObject, button::MouseButton, window)
    @materialize mouse_hover, mousebuttonspressed = window.inputs
    clicked_on = const_lift(mouse_hover, mousebuttonspressed) do mh, mbp
        mh.id == robj.id && in(button, mbp)
    end
    clicked_on_obj = keepwhen(clicked_on, false, clicked_on)
    clicked_on_obj = const_lift((mh, x)->(x,robj,mh), mouse_hover, clicked_on)
    clicked_on, clicked_on_obj
end
export is_same_id
is_same_id(id_index, id::Int) = id_index.id == id
is_same_id(id_index, robj::RenderObject) = id_index.id == robj.id
is_same_id(id_index, ids::Tuple) = id_index.id in ids
"""
Returns a boolean signal indicating if the mouse hovers over `robj`
"""
function is_hovering(robj::RenderObject, window)
    droprepeats(const_lift(is_same_id, window.inputs[:mouse_hover], robj))
end

"""
returns a signal which becomes true whenever there is a doublecklick
"""
function doubleclick(mouseclick, threshold::Real)
    ddclick = foldp((time(), Reactive.value(mouseclick), false), mouseclick) do v0, mclicked
        t0, lastc, _ = v0
        t1 = time()
        isclicked = (length(mclicked) == 1 &&
            length(lastc) == 1 &&
            first(lastc) == first(mclicked) &&
            t1-t0 < threshold
        )
        return (t1, mclicked, isclicked)
    end
    dd = const_lift(last, ddclick)
    return dd
end

export doubleclick
function default_camera_control(
        inputs, rotation_speed, translation_speed, keep = Signal(true)
    )
    @materialize mouseposition, mouse_buttons_pressed, scroll = inputs

    mouseposition = droprepeats(map(Vec2f0, mouseposition))
    left_pressed  = const_lift(pressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    right_pressed = const_lift(pressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_RIGHT)
    xytheta       = dragged_diff(mouseposition, left_pressed, keep)
    xytranslate   = dragged_diff(mouseposition, right_pressed, keep)

    ztranslate    = filterwhen(keep, 0f0,
        map(last, scroll)
    )
    translate_theta(
        xytranslate, ztranslate, xytheta,
        rotation_speed, translation_speed
    )
end

function thetalift(xy, speed)
    Vec3f0(xy[1], -xy[2], 0f0).*speed
end
function translationlift(up_left, zoom, speed)
    Vec3f0(zoom, up_left[1], up_left[2]).*speed
end
function diff_vector(v0, p1)
    p0, diff = v0
    p1, p0-p1
end
function translate_theta(
        xytranslate, ztranslate, xytheta,
        rotation_speed, translation_speed
    )
    theta = map(thetalift, xytheta, rotation_speed)
    trans = map(translationlift, xytranslate, ztranslate, translation_speed)
    theta, trans
end

"""
Creates a perspective camera from a dict of signals
Args:

inputs: Dict of signals, looking like this:
[
    :window_size            => Signal(Vec{2, Int}),
    :buttons_pressed        => Signal(Int[]),
    :mouse_buttons_pressed  => Signal(Int[]),
    :mouseposition          => mouseposition, -> Panning + Rotation
    :scroll_y               => Signal(0) -> Zoomig
]
eyeposition: Position of the camera
lookatvec: Point the camera looks at
"""
function PerspectiveCamera(
        inputs::Dict{Symbol,Any},
        eyeposition::Vec{3, T}, lookatvec::Vec{3, T};
        upvector = Vec3f0(0, 0, 1),
        keep=Signal(true), theta=nothing, trans=nothing
    ) where T
    lookat, eyepos = Signal(lookatvec), Signal(eyeposition)
    # TODO make this more elegant!
    _theta, _trans = default_camera_control(
        inputs, Signal(0.1f0), Signal(1f0), keep
    )
    theta = theta == nothing ? _theta : theta
    trans = trans == nothing ? _trans : trans
    farclip = map(eyepos, lookat) do a,b
        max(norm(b-a) * 5f0, 30f0)
    end
    minclip = map(eyepos, lookat) do a,b
        norm(b-a) * 0.007f0
    end
    PerspectiveCamera(
        theta,
        trans,
        lookat,
        eyepos,
        Signal(upvector),
        inputs[:window_area],
        Signal(41f0), # Field of View
        Signal(0.001f0),  # Min distance (clip distance)
        farclip # Max distance (clip distance)
    )
end
function PerspectiveCamera(
        area,
        eyeposition::Signal{Vec{3, T}}, lookatvec::Signal{Vec{3, T}}, upvector
    ) where T
    PerspectiveCamera(
        Signal(Vec3f0(0)),
        Signal(Vec3f0(0)),
        lookatvec,
        eyeposition,
        upvector,
        area,
        Signal(41f0), # Field of View
        Signal(0.1f0),  # Min distance (clip distance)
        Signal(50f0) # Max distance (clip distance)
    )
end

function projection_switch(
        wh::SimpleRectangle,
        fov::T, near::T, far::T,
        projection::Projection, zoom::T
    ) where T<:Real
    aspect = T(wh.w/wh.h)
    h      = T(tan(fov / 360.0 * pi) * near)
    w      = T(h * aspect)
    projection == PERSPECTIVE && return frustum(-w, w, -h, h, near, far)
    h, w   = h*zoom, w*zoom
    orthographicprojection(-w, w, -h, h, near, far)
end

w_component(::Point{N, T}) where {N, T} = T(1)
w_component(::Vec{N, T}) where {N, T} = T(0)

function to_worldspace(point::T, cam) where T <: StaticVector
    to_worldspace(
        point,
        Reactive.value(cam.projection) * Reactive.value(cam.view),
        T(widths(Reactive.value(cam.window_size)))
    )
end
function to_worldspace(
        p::StaticVector{N, T},
        projectionview::Mat4,
        cam_res::StaticVector
    ) where {N, T}
    VT = typeof(p)
    prj_view_inv = inv(projectionview)
    clip_space = T(4) * (VT(p) ./ VT(cam_res))
    pix_space = Vec{4, T}(
        clip_space[1],
        clip_space[2],
        T(0), w_component(p)
    )
    ws = prj_view_inv * pix_space
    ws # ./ ws[4]
end

"""
Takes a point and a camera and transforms it from mouse (imagespace) to world space
"""
function imagespace(pos, camera)
    # Setup transformation matrix
    pv = Reactive.value(camera.projection) * Reactive.value(camera.view)
    inv_pv = inv(pv)
    width, height = widths(Reactive.value(camera.window_size)) # get pixel resolution
    x, y = pos
    # transform to normalized device coordinates [-1, 1]
    device_space = Vec4f0(
        2 * (x / width)  - 1,
        2 * (y / height) - 1,
        0.0,
        1.0
    )
    pos = inv_pv * device_space
    Point2f0(pos[1], pos[2]) / pos[4]
end

function translate_cam(
        translate, proj, view, window_size, prj_type,
        eyepos_s, lookat_s, up_s,
    )
    translate == Vec3f0(0) && return nothing # nothing to do

    lookat, eyepos, up, prjt = map(value, (lookat_s, eyepos_s, up_s, prj_type))
    dir = eyepos - lookat
    dir_len = norm(dir)
    cam_res = Vec2f0(widths(Reactive.value(window_size)))

    zoom, x, y = translate
    zoom *= 0.1f0 * dir_len
    if prjt != PERSPECTIVE
        x, y = to_worldspace(Vec2f0(x, y), Reactive.value(proj) * Reactive.value(view), cam_res)
    else
        x, y = (Vec2f0(x, y) ./ cam_res) .* dir_len
    end
    dir_norm = normalize(dir)
    right = normalize(cross(dir_norm, up))
    zoom_trans = dir_norm*zoom
    side_trans = right * (-x) + normalize(up) * y
    push!(eyepos_s, eyepos + side_trans + zoom_trans)
    push!(lookat_s, lookat + side_trans)
    nothing
end

function rotate_cam(
        theta::Vec{3, T},
        cam_right::Vec{3,T}, cam_up::Vec{3,T}, cam_dir::Vec{3, T}
    ) where T
    rotation = one(Q.Quaternion{T})
    # first the rotation around up axis, since the other rotation should be relative to that rotation
    if theta[1] != 0
        rotation *= Q.qrotation(cam_up, theta[1])
    end
    # then right rotation
    if theta[2] != 0
        rotation *= Q.qrotation(cam_right, theta[2])
    end
    # last rotation around camera axis
    if theta[3] != 0
        rotation *= Q.qrotation(cam_dir, theta[3])
    end
    rotation
end
"""
Creates a perspective camera from signals, controlling the camera
Args:

`window_size`: Size of the window

fov: Field of View
nearclip: Near clip plane
farclip: Far clip plane
`theta`: rotation around camera axis
`trans`: translation in camera space (xyz are the camera axes)
`lookatposition`: point the camera looks at
`eyeposition`: the actual position of the camera (the lense, the \"eye\")
"""
function PerspectiveCamera(
        theta,
        trans::Signal{T},
        lookatposition::Signal{T},
        eyeposition::Signal{T},
        upvector::Signal{T},
        window_size,
        fov,
        nearclip,
        farclip,
        projectiontype = Signal(PERSPECTIVE)
    ) where T<:Vec3
    # we have three ways to manipulate the camera: rotation, lookat/eyeposition and translation
    positions = (eyeposition, lookatposition, upvector)

    zoomlen = map(norm, map(-, lookatposition, eyeposition))
    projectionmatrix = map(projection_switch,
        window_size, fov, nearclip,
        farclip, projectiontype, zoomlen
    )

    # create the vievmatrix with the help of the lookat function
    viewmatrix = map(lookat, eyeposition, lookatposition, upvector)
    projectionview = map(*, projectionmatrix, viewmatrix)

    preserve(map(translate_cam,
       trans, Signal(projectionmatrix), Signal(viewmatrix), Signal(window_size),
       Signal(projectiontype),  Signal(eyeposition), Signal(lookatposition),
       Signal(upvector)
    ))

    preserve(map(theta) do theta_v
        theta_v == Vec3f0(0) && return nothing #nothing to do!
        eyepos_v, lookat_v, up_v = map(value, positions)

        dir = normalize(eyepos_v - lookat_v)
        right_v = normalize(cross(up_v, dir))
        up_v  = normalize(cross(dir, right_v))

        rotation = rotate_cam(theta_v, right_v, Vec3f0(0, 0, sign(up_v[3])), dir)
        r_eyepos = lookat_v + rotation*(eyepos_v - lookat_v)
        r_up = normalize(rotation*up_v)
        push!(eyeposition, r_eyepos)
        push!(upvector, r_up)
    end)

    PerspectiveCamera{eltype(T)}(
        window_size,
        nearclip,
        farclip,
        fov,
        viewmatrix,
        projectionmatrix,
        projectionview,
        eyeposition, lookatposition, upvector,
        trans,
        theta,
        projectiontype
    )
end

"""
get's the boundingbox of a render object.
needs value, because boundingbox will always return a boundingbox signal
"""
signal_boundingbox(robj) = Reactive.value(boundingbox(robj))

"""
Calculates union boundingbox of all elements in renderlist
(Can't do ::Vector{RenderObject{T}}, because t is not always the same)
"""
function renderlist_boundingbox(renderlist::Vector)
    renderlist = filter(x-> x != nothing, renderlist)
    isempty(renderlist) && return AABB(Vec3f0(NaN), Vec3f0(0)) # nothing to do here
    robj1 = first(renderlist)
    bb = Reactive.value(robj1[:model])*signal_boundingbox(robj1)
    for elem in renderlist[2:end]
        bb = union(Reactive.value(elem[:model])*signal_boundingbox(elem), bb)
    end
    bb
end

"""
Centers the camera on a list of render objects
"""
function center!(camera::PerspectiveCamera, renderlist::Vector; border = 0)
    bb = renderlist_boundingbox(renderlist)
    bb = AABB(minimum(bb) .- border, widths(bb) .+ 2border)
    center!(camera, bb)
end
"""
Centers a camera onto a boundingbox
"""
function center!(camera::PerspectiveCamera, bb::AABB)
    width        = widths(bb)
    half_width   = width/2f0
    lower_corner = minimum(bb)
    middle       = maximum(bb) - half_width
    if Reactive.value(camera.projectiontype) == ORTHOGRAPHIC
        area, fov, near, far = map(value,
            (camera.window_size, camera.fov, camera.nearclip, camera.farclip)
        )
        aspect = Float32(area.w/area.h)
        h = Float32(tan(fov / 360.0 * pi) * near)
        w = h * aspect
        w_, h_, _ = half_width
        zoom = Vec2f0(w_, h_)./Vec2f0(w,h)
        x, y, _ = middle
        push!(camera.eyeposition, Vec3f0(x, y, maximum(zoom)))
        push!(camera.lookat, Vec3f0(x, y, 0))
        push!(camera.up, Vec3f0(0,1,0))
    else
        push!(camera.lookat, middle)
        neweyepos = middle + (width*1.2f0)
        push!(camera.eyeposition, neweyepos)
        push!(camera.up, Vec3f0(0,0,1))
        push!(camera.nearclip, 0.1f0 * norm(widths(bb)))
        push!(camera.farclip, 3f0 * norm(widths(bb)))
    end
    return
end
function robj_from_camera() end
function renderlist() end
export renderlist, robj_from_camera

"""
Centers the camera(=:perspective) on all render objects in `window`
"""
function center!(window, camera::Symbol=:perspective; border = 0)
    rl = robj_from_camera(window, camera)
    center!(window.cameras[camera], rl, border = border)
end


text:
using CanvasWebIO, WebIO, Mux

try
    global port += 1
catch
    global port = 8000
end

canvas = Canvas()
bg = dom"svg:rect[height=800, width=800, fill=blue]"()
box1 = dom"svg:rect[id=box1, height=50, width=50, x=50, y=50]"()
box2 = dom"svg:rect[id=box2, height=100, width=25, x=250, y=250]"()
box3 = dom"svg:rect[id=box3, height=100, width=25, x=550, y=350]"()
circ1 = dom"svg:circle[id=circ1, cx=200, cy=25, r=100]"()
addmovable!(canvas, box1)
addmovable!(canvas, box2)
addmovable!(canvas, circ1)
addclickable!(canvas, box3)
addstatic!(canvas, bg)

webio_serve(page("/", req -> canvas()), port)

ext: jl


text:
# input one single character
print("enter any character: ")
c = read(stdin,Char)
print("\n")
# selection ladder
if c in '0':'9'
   println("digit")
elseif c in  'A':'Z'
   println("uppercase letter")
elseif c in  'a':'z'
   println("lowercase letter")
else
   println("special character")
end

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

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

# =========================================================================================
"""
    get_json_safe(
        ::Type{T},
        json::AbstractDict,
        key::AbstractString
    )::Union{T, nothing}

Return a value from a JSON dictionary of the specified type.  Returns `false` if the item
doesn't exist or is not of the specified type.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | ::Type{T} |         | The data type expected for the value.
| R | json      |         | The JSON dictionary to pull the value from.
| R | key       |         | Key to look for in the dictionary.
"""
function get_json_safe(
            ::Type{T},
            json::AbstractDict,
            key::AbstractString
        )::Union{T, Nothing} where T

    if haskey(json, key) == false
        return nothing
    end

    local value = json[key]

    if isa(value, T) == false
        return nothing
    end

    return value
end

# =========================================================================================
"""
    create_data_dir(
        dir::AbstractString
    )::Bool

Return a value from a JSON dictionary of the specified type.  Returns `false` if the item
doesn't exist or is not of the specified type.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | ::Type{T} |         | The data type expected for the value.
| R | json      |         | The JSON dictionary to pull the value from.
| R | key       |         | Key to look for in the dictionary.
"""
function create_data_dir(
            dir::AbstractString
        )::Bool

    # -------------------------------------------------------------------------------------
    # Directory exists, our work is done.
    if isdir(dir) == true
        return true
    end

    # -------------------------------------------------------------------------------------
    # Make sure the path isn't something else
    if ispath(dir) == true
        @warn("Data directory is not a directory: $(abspath(dir))")
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the directory.
    local retval = false

    try
        mkpath(dir)
        retval = true
    catch ex
        if isa(ex, Base.IOError) == false
            rethrow(ex)
        end
        @warn("Could not create data directory: $(abspath(dir))")
    end

    return retval
end

# =========================================================================================
"""
    calc_sha1(
        fullPath::AbstractString
    )::String

Calculate the GitHub SHA1 for a file on disk.  Returns the SHA1 code for the file.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | fullpath |         | The full path to the file we want to read.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| FileError    | There was an issue reading from the file.

__Details:__

GitHub's SHA1 value is slightly different from a sha1sum calculated locally.  First the
text: "blob $size\0" is hashed followed by the contents of the file.  Size is replaced by
the file size in decimal.
"""
function calc_sha1(
            fullpath::AbstractString
        )::String
    local sha1 = SHA1_CTX()
    local size = filesize(fullpath)
    update!(sha1, Vector{UInt8}("blob $size\0"))

    try
        open(fullpath; read = true) do io
            local buffer = Vector{UInt8}(undef, 8192)

            while eof(io) == false
                local bytes = readbytes!(io, buffer)
                update!(sha1, buffer, bytes)
            end
        end
    catch ex
        if isa(ex, SystemError)
            @warn "Could not read file: $(fullpath)"
            throw(FileError(fullpath))
        end
        rethrow(ex)
    end

    return bytes2hex(digest!(sha1))
end


text:
function pca(X::Array, k::Int)

    m=size(X)
    n=size(X)

    μ=mean(X,dims=1)
    σ=std(X,dims=1)
    Xn=(X.-μ)./σ

    CovM=(1/m).*(Xn'*Xn)
    U,S,V=svd(CovM)
    Ur=U[:,1:k]
    Z=Xn*Ur

    return Z
end


text:
"""
    apply_title!(g, t, p, parent_font)

Internal function to apply a `Title` object `t` in a GLE context.
The argument `p` specifies the prefix.
"""
function apply_title!(g::GLE, t::Title, p::String, parent_font::String)
    # [x]title ...
    "\n\t$(p)title \"$(t.text)\""     |> g
    isdef(t.dist) && "dist $(t.dist)" |> g
    apply_textstyle!(g, t.textstyle, parent_font)
    return nothing
end

"""
    apply_ticks!(g, t, p)

Internal function to apply a `Ticks` object `t` in a GLE context for an axis
prefixed by `p`.
"""
function apply_ticks!(g::GLE, t::Ticks, prefix::String, parent_font::String)
    # [x]ticks ...
    "\n\t$(prefix)ticks" |> g
    if t.off
        "off" |> g
        # also discard the labels and don't call anything else
        "\n\t$(prefix)labels off" |> g
        return nothing
    end
    # - style
    isdef(t.length) && "length $(t.length)" |> g
    apply_linestyle!(g, t.linestyle)
    # [x]places pos1 pos2 ...
    isempty(t.places) || "\n\t$(prefix)places $(vec2str(t.places))" |> g
    # [x]xaxis symticks
    t.symticks && "\n\t$(prefix)axis symticks" |> g
    apply_tickslabels!(g, t.labels, prefix, parent_font)
    return nothing
end

"""
    apply_tickslabels!(g, t, p)

Internal function to apply a `TicksLabels` object `t` in a GLE context.
The prefix `p` indicates which axis we're on.
"""
function apply_tickslabels!(g::GLE, t::TicksLabels, prefix::String, parent_font::String)
    # [x]names "names1" ...
    isempty(t.names) || "\n\t$(prefix)names $(vec2str(t.names))" |> g
    # [x]labels ...
    "\n\t$(prefix)labels"                     |> g
    ifelse(t.off, "off", "")                  |> g
    isdef(t.dist) && "dist $(t.dist)"         |> g
    apply_textstyle!(g, t.textstyle, parent_font)
    # [x]axis ...
    if any(isdef, (t.angle, t.format))
        "\n\t$(prefix)axis" |> g
        isdef(t.angle)  && "angle $(t.angle)"   |> g
        isdef(t.format) && "format $(t.format)" |> g
        isdef(t.shift)  && "shift $(t.shift)"   |> g
    end
    return nothing
end


text:
# http://web.mit.edu/puzzle/www/2011/puzzles/world1/squared_key/answer/

@testset "squared key" begin
    f = best_feature(wordlist("""
        EIGHTH NOTE, CHAINS, JUNO, PLUTO, COMET, RAIN, FLOWER, SEXTILE, PICK, NEUTER, PISCES, FLEUR-DE-LIS, ALEMBIC, ANCHOR, HOT SPRINGS, SHAMROCK, EARTH, YIN YANG, WHITE FLAG, SATURN, HOT BEVERAGE"""))
    @test f.description == "is the name of a unicode character"
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdrate

"""
    'feffery'_antdrate(;kwargs...)

An AntdRate component.

Keyword arguments:
- `id` (String; optional)
- `allowClear` (Bool; optional)
- `allowHalf` (Bool; optional)
- `className` (String; optional)
- `count` (Real; optional)
- `defaultValue` (Real; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `style` (Dict; optional)
- `tooltips` (Array of Strings; optional)
- `value` (Real; optional)
"""
function 'feffery'_antdrate(; kwargs...)
        available_props = Symbol[:id, :allowClear, :allowHalf, :className, :count, :defaultValue, :disabled, :loading_state, :style, :tooltips, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdrate", "AntdRate", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
using Quantica
using Documenter

DocMeta.setdocmeta!(Quantica, :DocTestSetup, :(using Quantica); recursive=true)

makedocs(;
    modules=[Quantica],
    authors="Pablo San-Jose",
    repo="https://github.com/pablosanjose/Quantica.jl/blob/{commit}{path}#L{line}",
    sitename="Quantica.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pablosanjose.github.io/Quantica.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Examples" => "examples.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/pablosanjose/Quantica.jl",
)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FXdiv_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FXdiv")
JLLWrappers.@generate_main_file("FXdiv", UUID("04d6f918-1a1e-5862-bbec-f5b7bf9ea9e9"))
end  # module FXdiv_jll


text:
using LinearAlgebra

export translate, scale, rotate

import Base: position

"""
    translate(v::Vec2D) -> Matrix
Creates a 3x3 matrix that would translate points
in the `v` direction.
"""
function translate(v::Union{Vec2D, Point})
    A = zeros(3, 3) + I
	A[1, 3] = v.x
	A[2, 3] = v.y

    A
end

"""
    scale(x, y) -> Matrix
Create 3x3 matrix for scaling in object. If you multiply
this with a collection of points making up a shape, it
has the effect of scaling that shape.
"""
function scale(x::Number, y::Number)
    A = zeros(3, 3)
    A[1, 1] = x
    A[2, 2] = y
    A[3, 3] = 0.0
    A
end

"""
    rotate(angle) -> Matrix
Create a 3x3 matrix that can be used to rotate a shape
made up of multiple points along its axis. Each point
will have to be a column in matrix.
"""
function rotate(angle::Number)
    A = zeros(3, 3) + I
    sine   = sin(angle)
    cosine = cos(angle)

	A[1, 1] = cosine
	A[1, 2] = -sine
	A[2, 1] = sine
	A[2, 2] = cosine

    A
end


text:
module hughes_cantilever_examples
using FinEtools
using FinEtools.AlgoBaseModule: evalconvergencestudy
using FinEtools.AlgoDeforLinearModule: linearstatics, exportstresselementwise, exportstress

function evaluateerrors(filebase, modeldatasequence)
    println("")
    println("Stress RMS error")
    for md = modeldatasequence
        md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
    end
    elementsizes, errornorms, p = evalconvergencestudy(modeldatasequence)

    println("Normalized Approximate Error = $(errornorms)")

    f = log.(vec(errornorms))
    A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
    p = A \ f
    println("Linear log-log fit: p = $(p)")

    csvFile = filebase * "_Stress" * ".CSV"
    savecsv(csvFile,
    elementsizes=vec(elementsizes[1:end-1]),
    elementsizes2=vec(elementsizes[1:end-1].^2),
    elementsizes3=vec(elementsizes[1:end-1].^3),
    errornorms=vec(errornorms)
    )
    println("Wrote $csvFile")

    println("")
    println("Displacement RMS error")
    for md = modeldatasequence
        md["targetfields"] = [md["u"] for r in md["regions"]]
    end
    elementsizes, errornorms, p = evalconvergencestudy(modeldatasequence)

    println("Normalized Approximate Error = $(errornorms)")

    f = log.(vec(errornorms))
    A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
    p = A \ f
    println("Linear log-log fit: p = $(p)")

    csvFile = filebase * "_Displ" * ".CSV"
    savecsv(csvFile,
    elementsizes=vec(elementsizes[1:end-1]),
    elementsizes2=vec(elementsizes[1:end-1].^2),
    elementsizes3=vec(elementsizes[1:end-1].^3),
    errornorms=vec(errornorms)
    )
    println("Wrote $csvFile")
end

function hughes_cantilever_stresses_MST10()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)

    # Isotropic material
    E=1.0;
    nu=0.3;
    # nu=0.499999999;
    P=1.0;
    L=16.0;
    c=2.0;
    h=1.0;
    smult=0;
    E1=E/(1-nu^2); nu1=nu/(1-nu); I=(2*c)^3/12;
    uexx(x,y) = (P/(6*E1*I)/2*(-y).*(3*(L^2-(L-x).^2)+(2+nu1)*(y.^2-c^2)));
    uexy(x,y) = (P/(6*E1*I)/2*(((L-x).^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x).*y.^2));

    CTE = 0.0

    tolerance = 0.00001*c

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linspace(0.0, L, nL+1))
        ys = collect(linspace(0.0, h, nh+1))
        zs = collect(linspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copy!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
        IntegData(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
        end

        function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
        end

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip)")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_MST10

function hughes_cantilever_stresses_MST10_incompressible()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)

    # Isotropic material
    E=1.0;
    # nu=0.3;
    nu=0.499999;
    P=1.0;
    L=16.0;
    c=2.0;
    h=1.0;
    smult=0;
    E1=E/(1-nu^2); nu1=nu/(1-nu); I=(2*c)^3/12;
    uexx(x,y) = (P/(6*E1*I)/2*(-y).*(3*(L^2-(L-x).^2)+(2+nu1)*(y.^2-c^2)));
    uexy(x,y) = (P/(6*E1*I)/2*(((L-x).^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x).*y.^2));

    CTE = 0.0

    tolerance = 0.00001*c

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linspace(0.0, L, nL+1))
        ys = collect(linspace(0.0, h, nh+1))
        zs = collect(linspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,  0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copy!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
        IntegData(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
        end

        function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
        end

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip)")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_incompressible_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_MST10_incompressible

function hughes_cantilever_stresses_nodal_MST10()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)

    # Isotropic material
    E=1.0;
    nu=0.3;
    # nu=0.499999999;
    P=1.0;
    L=16.0;
    c=2.0;
    h=1.0;
    smult=0;
    E1=E/(1-nu^2); nu1=nu/(1-nu); I=(2*c)^3/12;
    uexx(x,y) = (P/(6*E1*I)/2*(-y).*(3*(L^2-(L-x).^2)+(2+nu1)*(y.^2-c^2)));
    uexy(x,y) = (P/(6*E1*I)/2*(((L-x).^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x).*y.^2));

    CTE = 0.0

    tolerance = 0.00001*c

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linspace(0.0, L, nL+1))
        ys = collect(linspace(0.0, h, nh+1))
        zs = collect(linspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,   0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copy!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4) # rule for tetrahedral meshes

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
        IntegData(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
        end

        function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
        end

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip)")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_nodal_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
        modeldata = exportstress(modeldata)
        stressfields = NodalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_nodal_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_nodal_MST10

function hughes_cantilever_stresses_nodal_T10()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)

    # Isotropic material
    E=1.0;
    nu=0.3;
    # nu=0.499999999;
    P=1.0;
    L=16.0;
    c=2.0;
    h=1.0;
    smult=0;
    E1=E/(1-nu^2); nu1=nu/(1-nu); I=(2*c)^3/12;
    uexx(x,y) = (P/(6*E1*I)/2*(-y).*(3*(L^2-(L-x).^2)+(2+nu1)*(y.^2-c^2)));
    uexy(x,y) = (P/(6*E1*I)/2*(((L-x).^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x).*y.^2));

    CTE = 0.0

    tolerance = 0.00001*c

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linspace(0.0, L, nL+1))
        ys = collect(linspace(0.0, h, nh+1))
        zs = collect(linspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j? 1.0: 0.0 for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copy!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR,
        IntegData(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
        end

        function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
        end

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip)")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_nodal_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstress(modeldata)
        stressfields = NodalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_nodal_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_nodal_T10

function hughes_cantilever_stresses_T10()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)

    # Isotropic material
    E=1.0;
    nu=0.3;
    # nu=0.499999999;
    P=1.0;
    L=16.0;
    c=2.0;
    h=1.0;
    smult=0;
    E1=E/(1-nu^2); nu1=nu/(1-nu); I=(2*c)^3/12;
    uexx(x,y) = (P/(6*E1*I)/2*(-y).*(3*(L^2-(L-x).^2)+(2+nu1)*(y.^2-c^2)));
    uexy(x,y) = (P/(6*E1*I)/2*(((L-x).^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x).*y.^2));

    CTE = 0.0

    tolerance = 0.00001*c

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linspace(0.0, L, nL+1))
        ys = collect(linspace(0.0, h, nh+1))
        zs = collect(linspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copy!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR,
        IntegData(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
        end

        function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
        end

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip)")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_T10

function hughes_cantilever_stresses_T10_incompressible()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)

    # Isotropic material
    E=1.0;
    # nu=0.3;
    nu=0.499999;
    P=1.0;
    L=16.0;
    c=2.0;
    h=1.0;
    smult=0;
    E1=E/(1-nu^2); nu1=nu/(1-nu); I=(2*c)^3/12;
    uexx(x,y) = (P/(6*E1*I)/2*(-y).*(3*(L^2-(L-x).^2)+(2+nu1)*(y.^2-c^2)));
    uexy(x,y) = (P/(6*E1*I)/2*(((L-x).^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x).*y.^2));

    CTE = 0.0

    tolerance = 0.00001*c

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linspace(0.0, L, nL+1))
        ys = collect(linspace(0.0, h, nh+1))
        zs = collect(linspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,  0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copy!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR,
        IntegData(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
        end

        function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            W=2/3*c^3;
            copy!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
        end

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegData(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip)")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_incompressible_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_T10_incompressible

function allrun()
    println("#####################################################")
    println("# hughes_cantilever_stresses_MST10 ")
    hughes_cantilever_stresses_MST10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_MST10_incompressible ")
    hughes_cantilever_stresses_MST10_incompressible()
    println("#####################################################")
    println("# hughes_cantilever_stresses_nodal_MST10 ")
    hughes_cantilever_stresses_nodal_MST10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_nodal_T10 ")
    hughes_cantilever_stresses_nodal_T10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_T10 ")
    hughes_cantilever_stresses_T10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_T10_incompressible ")
    hughes_cantilever_stresses_T10_incompressible()
    return true
end # function allrun

end # module hughes_cantilever_examples


text:
#!/usr/bin/julia
using DataFrames

df = readtable("data.csv")

df1 = by(df, :pid) do df
      times = df[:time]
      dt = Array(Int64, size(times, 1)-1)
      for i in 1:size(dt, 1)
        t0 = times[i]
        t1 = times[i+1]

        dt[i] = t1-t0
      end
      DataFrame(μ = mean(dt), var = var(dt), min = minimum(times))
end

stimes = df1[:min]
mdts = Array(Int64, size(stimes, 1)-1)

for i in 1:size(mdts, 1)
  t0 = stimes[i]
  t1 = stimes[i+1]

  mdts[i] = t1-t0
end

println("Mean ΔT: ", mean(mdts))
println("Var ΔT: ", var(mdts))
println("Mean of mean δt: ", mean(df1[:μ]))
println("Var of mean δt: ", var(df1[:μ]))
println("Mean of var δt: ", mean(df1[:var]))
println("Var of var δt: ", var(df1[:var]))

max_stars_repo_name: Hack-a-thingie/Wall-Scriptor

max_issues_repo_name: Hack-a-thingie/Wall-Scriptor

max_forks_repo_name: Hack-a-thingie/Wall-Scriptor


text:
abstract type YTDLDeps end
struct YTDLbin <: YTDLDeps end
struct FFMPEGbin <: YTDLDeps end

# test availability of COMMANDS (copied from Franklin.jl)
isexcutable(cmd)::Bool = try success(cmd); catch; false; end

isavailable(::Type{YTDLbin}) = isexcutable(`youtube-dl --help`)

isavailable(::Type{FFMPEGbin}) = isexcutable(`ffmpeg -h`)

function deps()
    yt = try
        `youtube-dl --version` |> readlines |> first
    catch
        "unavailable"
    end
    ff = try
        `ffmpeg -version` |> readlines |> first |> x->split(x," ",limit=4)[3]
    catch
        "unavailable"
    end
    @info("Dependencies' version",
        ytdl = yt,
        ffmpeg = ff
    )
    return nothing
end


text:
# Snippet 09/clip-03.jl

using StatisticalRethinking

ProjDir = @__DIR__
cd(ProjDir)

# ### snippet 9.3

# Construct the logdensity problem

struct clip_9_3_model{TY <: AbstractVector, TX <: AbstractVector}
    "Observations."
    y::TY
    "Covariate"
    x::TX
end

# Make the type callable with the parameters *as a single argument*.

function (problem:: clip_9_3_model)(θ)
    @unpack y, x, = problem    # extract the data
    @unpack muy, mux = θ     # works on the named tuple too
    ll = 0.0
    ll += loglikelihood(Normal(mux, 1), x)
    ll += loglikelihood(Normal(muy, 1), y)
    ll += logpdf(Normal(0, 1), mux)
    ll += logpdf(Normal(0, 1), muy)
    ll
end

# Instantiate the model with data and inits.

Random.seed!(1234591)

N = 100
x = rand(Normal(0, 1), N)
y = rand(Normal(0, 1), N)

p = clip_9_3_model(y, x)
θ = (muy = 0.0, mux=0.0)
p(θ)

# Write a function to return properly dimensioned transformation.

problem_transformation(p::clip_9_3_model) =
    as((muy = asℝ, mux = asℝ))

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P);

# Generate and show fig 9.3

fig = StatisticalRethinking.generate_n_samples(p, ∇P);

# End of `09/clip-03d.jl`


text:
# Thanks @tlienart

using Markdown

function hfun_doc(params)
    fname = params[1]
    head = length(params) > 1 ? params[2] : fname
    type = length(params) == 3 ? params[3] : ""
    doc = eval(Meta.parse("using SymbolicUtils; @doc SymbolicUtils.$fname"))
    txt = Markdown.plain(doc)
    # possibly further processing here
    body = Franklin.fd2html(txt, internal=true)
    return """

            $head

$type

$body

    """
end


text:
@testset "$(format("Concept: MaybeMissing[Conversion][VecOrMat]"))" begin
    let
        tc1     = [1,2,3,4]
        output1 = convert(MaybeMissing{Number},tc1)
        @test isa(output1,Array{MaybeMissing{Number}}) == true
        @test output1 == tc1
        #see(output1)
        output2 = convert(MaybeMissing{Float64},tc1)
        @test isa(output2, Array{MaybeMissing{Float64}}) == true
        #display(output2)
        # test array of missing
        tc2 = [missing, missing, missing]
        output3 = convert(MaybeMissing{Float64},tc2)
        @test isa(output3,Array{MaybeMissing{Float64}}) == true
        output3[1] = 1.0
        @test output3[1] == 1.0 && ismissing(output3[2]) && ismissing(output3[3])
        @test isa(convert(MaybeMissing{Float64},[missing missing;missing missing]),
            Array{MaybeMissing{Float64}}) == true
    end
end


text:
using LinkedLists
using Test

ll = LinkedList(["apple", "pear", "banana"])
state = start(ll)
insert!(ll, state, "kiwi")

@test first(ll) == "kiwi"


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

SEQUENCES = [
#I000002,L000002,
V000005,
V000010,
G000032,L000032,
L000041,V000041,
F000045,G000045,I000045,L000045,R000045,V000045,is000045,
G000073,L000073,
L000085,
G000108,L000108,
L000110,V000110,
V000111,
L000118,
V000120,
L000122,
L000132,
L000141,
V000142,
L000143,
L000144,
L000145,
V000165,
V000203,
V000217,
V000255,
G000257,L000257,
V000262,
V000265,
V000290,
V000292,
V000326,
V000330,
T000369,
V000384,
V000407,
V000566,
V000567,
#L000587,V000587,
L000594,
L000698,
L000712,
L000716,
L000720,
L000727,
L000728,
L000729,
L000730,
L000731,
L000735,
L000739,
V000796,
G000957,L000957,
L000961,
G001003,L001003,
G001006,L001006,
G001045,L001045,
V001106,
V001107,
V001147,
V001157,
V001221,
V001222,
V001316,
V001339,
L001464,
V001477,
T001497,
V001511,
L001710,
V001783,
V001790,
V001803,
V001813,
V001855,
L001934,
F002093,I002093,L002093,V002093,
L002107,
V002110,
#L002115,
V002411,
V002412,
V002413,
V002414,
G002426,L002426,
L002445,V002445,
L002448,
L002476,
L002808,
L002858,
F002997,I002997,L002997,
V003188,
I003277,L003277,
V003314,
L003319,
L003723,
V003817,
L004018,
L004402,
L004403,
L004404,
L004405,
L004406,
L004407,
L004408,
L004409,
L004410,
L004411,
L004412,
L004413,
L004414,
L004415,
L004416,
L004417,
L004418,
L004419,
L004420,
L004421,
L004422,
L004423,
L004424,
L004425,
T004747,
G005043,L005043,
#L005046,V005046,
F005101,I005101,L005101,V005101,is005101,
L005117,
V005187,
L005411,
L005412,
L005425,
V005564,
L005651,
L005758,
L005773,
L005875,
L006012,
V006053,
L006171,V006171,
L006229,
V006257,
G006318,L006318,
V006519,
is206864,
L006922,
is206942,
V007060,
T007318,
L007522,
V007559,
V007584,
V007585,
V007696,
V007814,
V007917,
L008279,M008279,T008279,V008279,
L008451,
L008452,
L008453,
V008544,
V008545,
V008683,
L008784,
L010815,
L010816,
L010817,
L010819,
L010820,
L010821,
L010822,
L010823,
L010824,
L010825,
L010826,
L010827,
L010828,
L010829,
L010830,
L010831,
L010832,
L010833,
L010834,
L010835,
L010836,
L010837,
L010838,
L010839,
L010840,
L010841,
T011117,
V011371,
L013928,
V014107,
#L014606,
L015128,
L023003,
L023004,
L023005,
L023006,
L023007,
L023008,
L023009,
L023010,
L023011,
L023012,
L023013,
L023014,
L023015,
L023016,
L023017,
L023018,
L023019,
L023020,
L023021,
is325143,
is325145,
L025528,
V026791,
V027641,
V027642,
is327988,
#L028296,
L031363,
#T131689,
V032031,
V033156,
V033815,
L034017,
V034386,
V034444,
F034885,I034885,L034885,V034885,
L035016,
L035251,
V035327,
L036038,
#L036040,
T037027,
V038554,
V038712,
L038872,
L038873,
T039599,
T039683,
L042965,
V046161,
V047053,
L047891,
V048724,
V048735,
L048993,M048993,T048993,V048993,
T049218,
T049310,
I050384,L050384,
T051141,
T051142,
V051933,
L052847,V052847,
T053117,
T053120,
T053121,
V054248,
V055634,
T055883,
V056982,
L057126,
L057127,
V058373,
T059419,
T060693,
V060818,
V061142,
V061168,
L061256,V061256,
L062980,
T064189,
L065515,
L065619,V065619,
L065855,
T066325,
V066570,
V067998,
L068228,
G068875,L068875,
L069637,
I072233,L072233,
L078760,
V080079,
V080577,
V080940,
V080956,
V081125,
L082298,
L082301,
L082302,
L082556,
L082557,
L082558,
L082559,
F083207,I083207,L083207,V083207,is083207,
V083652,
L084916,
T084938,
V086799,
#L088218,
T088969,
L089270,
T090238,
T090981,
is293654,
#L094088,
L094587,M094587,T094587,V094587,
T094665,
V095794,
is296095,
L096727,
V097383,
#L097805,M097805,P097805,T097805,V097805,
V099022,
T099174,
is299498,
is299733,
is299928,
is299929,
is299930,
V102376,
T104556,
T104562,
L104794,
L105523,
T106566,
L107716,
L107895,V107895,
L108624,
T111062,
L111529,
L111530,
L111531,
L111532,
L111533,
T111593,
T116392,
T119274,
V120738,
V123753,
V124320,
V124399,
V124441,
V124442,
V129760,
#T131689,P131689,
T132062,
L132393,M132393,T132393,V132393,
V135481,
V135521,
T137286,
T137338,
T137452,
T137513,
L141158,
V142151,
L146559,
V151800,
#L156289,
T159834,
V163590,
V163617,
L167872,
T172094,
V178963,
V178964,
V181830,
V181936,
T184962,
V188021,
L190905,V190905,
F193338,I193338,L193338,V193338,
F193339,I193339,L193339,V193339,
V195441,
T202327,
T203412,
F206864,I206864,L206864,
F206942,I206942,L206942,
#L210029,
#L211212,
V213370,
L213384,
L214551,
T216916,
V216919,
T217537,
L217924,
T225478,
F225498,I225498,L225498,
V231181,
V232980,
V232981,
V232982,
#T241171,P241171,
L242660,
L243168,
#L243655,
#L243664,
#L243665,
L244291,
L244779,
L244780,
L244819,
L246547,
L246655,
V250283,
V251568,
V254749,
#L257490,
V257993,
#L260884,V260884,
L262071,
T264428,
T265606,
V265609,
L271703,M271703,T271703,V271703,
L275585,V275585,
L276285,
L276286,
L276287,
V277233,
#L278073,P278073,
#L278074,P278074,
L280172,
#L281478,
#L281479,
#L281480,
V281586,
V281588,
L290351,V290351,
#L291451,
#L291452,
#L291973,V291973,
#V291974,
#V291975,
#V291976,
L292222,
F293654,I293654,L293654,
V295513,
F296095,I296095,L296095,
L299214,V299214,
F299498,I299498,L299498,
L299733,
F299928,I299928,L299928,
F299929,I299929,L299929,
F299930,I299930,L299930,
V309896,
L319284,
L319514,V319514,
T321620,
T321621,
T321623,
T321624,
T322942,
F325143,I325143,L325143,
F325145,I325145,L325145,
#L327003,
#L327004,
V327491,
V327492,
V327493,
V327494,
V327495,
V327987,
L327988,
isAbundant,
isCarmichael,
isComposite,
isCyclic,
isEven,
isNonnegative,
isOdd,
isPerfectPower,
isPositive,
isPowerOfPrimes,
isPrime,
isPrimePower,
#isPrimeTo,
isSquare,
isSquareFree,
isStrongCyclic,
#isStrongPrimeTo,
isTriangular,
#isUlam,
isZumkeller,
isweakCarmichael,
Ω,
σ,
σ2,
τ,
ω,
ϕ,
↑,
↓,
∏,
⊥,
⍊
]

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# optimized parsing functions, copied from julia/base/parse.jl

import Base: tryparse, tryparse_internal

function tryparse(::Type{Float64}, s::DenseStringViewAndSub)
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
                          (Ptr{UInt8},Csize_t,Csize_t), s, 0, sizeof(s))
    hasvalue ? val : nothing
end
function tryparse_internal(::Type{Float64}, s::DenseStringViewAndSub, startpos::Int, endpos::Int)
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
                          (Ptr{UInt8},Csize_t,Csize_t), s, startpos-1, endpos-startpos+1)
    hasvalue ? val : nothing
end
function tryparse(::Type{Float32}, s::DenseStringViewAndSub)
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
                          (Ptr{UInt8},Csize_t,Csize_t), s, 0, sizeof(s))
    hasvalue ? val : nothing
end
function tryparse_internal(::Type{Float32}, s::DenseStringViewAndSub, startpos::Int, endpos::Int)
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
                          (Ptr{UInt8},Csize_t,Csize_t), s, startpos-1, endpos-startpos+1)
    hasvalue ? val : nothing
end

function tryparse_internal(::Type{Complex{T}}, s::DenseStringViewAndSub, i::Int, e::Int, raise::Bool) where {T<:Real}
    # skip initial whitespace
    while i ≤ e && isspace(s[i])
        i = nextind(s, i)
    end
    if i > e
        raise && throw(ArgumentError("input string is empty or only contains whitespace"))
        return nothing
    end

    # find index of ± separating real/imaginary parts (if any)
    i₊ = something(findnext(in(('+','-')), s, i), 0)
    if i₊ == i # leading ± sign
        i₊ = something(findnext(in(('+','-')), s, i₊+1), 0)
    end
    if i₊ != 0 && s[i₊-1] in ('e','E') # exponent sign
        i₊ = something(findnext(in(('+','-')), s, i₊+1), 0)
    end

    # find trailing im/i/j
    iᵢ = something(findprev(in(('m','i','j')), s, e), 0)
    if iᵢ > 0 && s[iᵢ] == 'm' # im
        iᵢ -= 1
        if s[iᵢ] != 'i'
            raise && throw(ArgumentError("expected trailing \"im\", found only \"m\""))
            return nothing
        end
    end

    if i₊ == 0 # purely real or imaginary value
        if iᵢ > i && !(iᵢ == i+1 && s[i] in ('+','-')) # purely imaginary (not "±inf")
            x = tryparse_internal(T, s, i, iᵢ-1, raise)
            x === nothing && return nothing
            return Complex{T}(zero(x),x)
        else # purely real
            x = tryparse_internal(T, s, i, e, raise)
            x === nothing && return nothing
            return Complex{T}(x)
        end
    end

    if iᵢ < i₊
        raise && throw(ArgumentError("missing imaginary unit"))
        return nothing # no imaginary part
    end

    # parse real part
    re = tryparse_internal(T, s, i, i₊-1, raise)
    re === nothing && return nothing

    # parse imaginary part
    im = tryparse_internal(T, s, i₊+1, iᵢ-1, raise)
    im === nothing && return nothing

    return Complex{T}(re, s[i₊]=='-' ? -im : im)
end

function tryparse_internal(::Type{Bool}, sbuff::DenseStringViewAndSub,
        startpos::Int, endpos::Int, base::Integer, raise::Bool)
    if isempty(sbuff)
        raise && throw(ArgumentError("input string is empty"))
        return nothing
    end

    if isnumeric(sbuff[1])
        intres = tryparse_internal(UInt8, sbuff, startpos, endpos, base, false)
        (intres == 1) && return true
        (intres == 0) && return false
        raise && throw(ArgumentError("invalid Bool representation: $(repr(sbuff))"))
    end

    orig_start = startpos
    orig_end   = endpos

    # Ignore leading and trailing whitespace
    while isspace(sbuff[startpos]) && startpos <= endpos
        startpos = nextind(sbuff, startpos)
    end
    while isspace(sbuff[endpos]) && endpos >= startpos
        endpos = prevind(sbuff, endpos)
    end

    len = endpos - startpos + 1
    p   = pointer(sbuff) + startpos - 1
    GC.@preserve sbuff begin
        (len == 4) && (0 == _memcmp(p, "true", 4)) && (return true)
        (len == 5) && (0 == _memcmp(p, "false", 5)) && (return false)
    end

    if raise
        substr = SubString(sbuff, orig_start, orig_end) # show input string in the error to avoid confusion
        if all(isspace, substr)
            throw(ArgumentError("input string only contains whitespace"))
        else
            throw(ArgumentError("invalid Bool representation: $(repr(substr))"))
        end
    end
    return nothing
end


text:
using AdventOfCode, Test, Mocking, HTTP, Dates
Mocking.activate()
@testset "All tests" begin
    AOC = AdventOfCode
    get_patch = @patch HTTP.get(a...; kw...) = (status = 200, body = "TESTDATA")
    get_error_patch = @patch HTTP.get(a...; kw...) = (status = 404, body = "ERROR")

    @testset "Basic Operations" begin
        @testset "_is_unlocked" begin
            @test AOC._is_unlocked(2018, 1)
            @test (@test_logs (:warn, "Advent of Code for year 3000 and day 1 hasn't unlocked yet.") !AOC._is_unlocked(3000, 1))
        end
        @testset "_get_cookies" begin
            cleanup = false
            if "AOC_SESSION" ∈ keys(ENV)
                temp = pop!(ENV, "AOC_SESSION")
                cleanup = true
            end
            @test_throws ErrorException("Session cookie in ENV[\"AOC_SESSION\"] needed to download data.") AOC._get_cookies()
            ENV["AOC_SESSION"] = "COOKIE"
            @test AOC._get_cookies() == Dict("session" => "COOKIE")

            if cleanup
                ENV["AOC_SESSION"] = temp
            end
        end
        @testset "_base_url" begin
            @test AOC._base_url(2019, 1) == "https://adventofcode.com/2019/day/1"
        end
        @testset "_download_data" begin
            apply(get_patch) do
                @test AOC._download_data(2019, 1) == "TESTDATA"
            end
            apply(get_error_patch) do
                @test_throws ErrorException("Unable to download data") AOC._download_data(2019, 1)
            end
        end
    end
    for include_year in (true, false)
        include_year ? base_path = joinpath(@__DIR__, "2019") : base_path = @__DIR__
        data_path = joinpath(base_path, "data")
        src_path = joinpath(base_path, "src")

        @testset "Include Year = $include_year" begin
            @testset "_template" begin
                @test AOC._template(2019, 1; include_year = include_year) == """
                # https://adventofcode.com/2019/day/1
                using AdventOfCode

                input = readlines("$(relpath(joinpath(data_path, "day_1.txt")))")

                function part_1(input)
                    nothing
                end
                @info part_1(input)

                function part_2(input)
                    nothing
                end
                @info part_2(input)
                """
            end
            @testset "setup_files" begin
                @testset "_setup_data_file" begin
                    apply(get_patch) do
                        try
                            AOC._setup_data_file(2019, 1, include_year = include_year)
                            @test isdir(data_path)
                            @test isfile(joinpath(data_path, "day_1.txt"))
                            @test readlines(joinpath(data_path, "day_1.txt")) == ["TESTDATA"]
                            @test_logs (:warn, r"will not redownload") AOC._setup_data_file(2019, 1; include_year = include_year)
                        finally
                            rm(data_path, recursive = true, force = true)
                        end
                    end
                end

                @testset "works for full method" begin
                    apply(get_patch) do
                        try
                            path = AOC.setup_files(2019, 1, force = true, include_year = include_year)
                            joinpath(src_path, "day_1.jl")
                            @test isdir(data_path)
                            @test isfile(joinpath(data_path, "day_1.txt"))
                            @test readlines(joinpath(data_path, "day_1.txt")) == ["TESTDATA"]
                            @test isdir(src_path)
                            @test isfile(joinpath(src_path, "day_1.jl"))
                            @test String(read(joinpath(src_path, "day_1.jl"))) == AOC._template(2019, 1,include_year = include_year)
                            @test_logs (:warn, Regex("day_1.txt already exists.")) (:warn, Regex("day_1.jl already exists.")) AOC.setup_files(2019, 1, force = false, include_year = include_year)
                        finally
                            if include_year
                                rm(base_path, recursive = true, force = true)
                            else
                                rm(data_path, recursive = true, force = true)
                                rm(src_path, recursive = true, force = true)
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
"""
`module  JAC.NonresonantTwoColorPhotonIon`  ... a submodel of JAC that contains all methods for computing nonresonant two-color
                                                (two-photon & single-electron) ionization cross sections; it is using JAC.
"""
module NonresonantTwoColorPhotonIon

    using JAC, JAC.ManyElectron

    """
    `struct  NonresonantTwoColorPhotonIon.Settings`  ... defines a type for the settings in estimating resonant two-color (two-photon &
                                                         single-electron) ionization cross sections.

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + photonEnergies          ::Array{Float64,1}             ... List of photon energies.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).

    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        photonEnergies            ::Array{Float64,1}
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `JAC.NonresonantTwoColorPhotonIon.Settings()`  ... constructor for the default values of resonant two-color (two-photon &
                                                       single-electron) ionization estimates.
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], Float64[], false, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::NonresonantTwoColorPhotonIon.Settings)`
    #		... prepares a proper printout of the variable settings::NonresonantTwoColorPhotonIon.Settings.
    function Base.show(io::IO, settings::NonresonantTwoColorPhotonIon.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "photonEnergies:           $(settings.photonEnergies)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

    """
    `struct  NonresonantTwoColorPhotonIon.Channel`  ... defines a type for a two-color channel to help characterize resonant two-color
                                                        (two-photon & single-electron) ionization with well-defined multipolarities.

        + multipoles     ::Array{EmMultipole,1}   ... Multipoles of all N incoming photons.
        + gauge          ::EmGauge                ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                  ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry          ... total angular momentum and parity of the scattering state
        + phase          ::Float64                ... phase of the partial wave
        + amplitude      ::Complex{Float64}       ... multi-photon ionization amplitude associated with the given channel.
   """
    struct  Channel
        multipoles       ::Array{EmMultipole,1}
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  NonresonantTwoColorPhotonIon.Line`  ... defines a type for a multi-photon ionization line that may include the definition of channels.

        + initialLevel     ::Level                  ... initial-(state) level
        + finalLevel       ::Level                  ... final-(state) level
        + photonEnergy     ::Float64                ... Energy of the incoming photons; all photons are assumed to have equal energy.
        + crossSection     ::EmProperty             ... Cross section for this multi-photon ionization.
        + hasChannels      ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                        multipolarities, etc., or not.
        + channels         ::Array{NonresonantTwoColorPhotonIon.Channel,1}  ... List of NonresonantTwoColorPhotonIon.Channels of this line.
    """
    struct  Line
        initialLevel       ::Level
        finalLevel         ::Level
        photonEnergy       ::Float64
        crossSection       ::EmProperty
        hasChannels        ::Bool
        channels           ::Array{NonresonantTwoColorPhotonIon.Channel,1}
    end

    # `Base.show(io::IO, line::NonresonantTwoColorPhotonIon.Line)`
    #		... prepares a proper printout of the variable  line::NonresonantTwoColorPhotonIon.Line.
    function Base.show(io::IO, line::NonresonantTwoColorPhotonIon.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "photonEnergy:      $(line.photonEnergy)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
        println(io, "hasChannels:       $(line.hasChannels)  ")
        println(io, "channels:          $(line.channels)  ")
    end

end # module


text:
using Catalyst, Graphviz_jll

rn = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
    S^2, R --> 0
end α β

# check can make a graph
gr = Graph(rn)

# check can save a graph
fname = Base.Filesystem.tempname()
savegraph(gr, "$fname.svg", "svg")

# this are broken in the jll, see
# https://github.com/JuliaPackaging/Yggdrasil/issues/1428
# savegraph(gr, "$fname.pdf", "pdf")
# savegraph(gr, "$fname.png", "png")


text:
"""
Reverse function for `sinh`.
"""
function sinh_rev(y::Interval, x::Interval)
    x = x ∩ asinh(y)

    return y, x
end

"""
Reverse function for `cosh`.
"""
function cosh_rev(y::Interval,x::Interval)
    y_new = y ∩ Interval(1.,∞)
    x = x ∩ acosh(y)

    return y_new, x
end

"""
Reverse function for `tanh`.
"""
function tanh_rev(y::Interval,x::Interval)
    y_new = y ∩ Interval(-1.,1.)
    x = x ∩ atanh(y)

    return y_new, x
end


text:
using DataGraphs
using NamedGraphs
using Dictionaries
using Graphs

g = NamedDimGraph(grid((2, 2)); dims=(2, 2))
dg = NamedDimDataGraph{String,String}(g)

dg[1, 1] = "X11"

@show dg[1, 1] == "X11"

dg[(1, 1) => (1, 2)] = "X11↔X12"

@show dg[(1, 1) => (1, 2)] == "X11↔X12"
@show dg[(1, 2) => (1, 1)] == "X11↔X12"

# XXX: Broken
#@show isassigned(dg, (1, 1))
#@show isassigned(dg, 1, 1)

@show isassigned(dg, NamedDimEdge((1, 1), (1, 2)))
@show !isassigned(dg, NamedDimEdge((1, 1), (2, 2)))
@show isassigned(dg, (1, 1) => (1, 2))
@show !isassigned(dg, (1, 1) => (2, 2))

## @show has_edge(dg, 1, 2)
## @show has_edge(dg, 1 => 2)
## @show !has_edge(dg, 1, 3)
## @show !has_edge(dg, 1 => 3)
## @show has_vertex(dg, 1)
## @show has_vertex(dg, 4)
## @show !has_vertex(dg, 0)
## @show !has_vertex(dg, 5)
##
## dg[1] = "V1"
## dg[2] = "V2"
## dg[3] = "V3"
## dg[4] = "V4"
## @show isassigned(dg, 1)
## @show dg[1] == "V1"
## @show dg[2] == "V2"
## @show dg[3] == "V3"
## @show dg[4] == "V4"
##
## dg[1 => 2] = :E12
## dg[2 => 3] = :E23
## dg[Edge(3, 4)] = :E34
## #@show isassigned(dg, (1, 2))
## @show isassigned(dg, Edge(2, 3))
## @show isassigned(dg, 3 => 4)
## @show dg[Edge(1, 2)] == :E12
## @show dg[2 => 3] == :E23
## @show dg[3 => 4] == :E34


text:
using Test

import MathOptSetDistances
const MOD = MathOptSetDistances
const MOI = MathOptSetDistances.MOI

using LinearAlgebra
using Random
import BlockDiagonals: BlockDiagonal

include("distances.jl")
include("projections.jl")
include("projection_gradients.jl")
include("chainrules.jl")


text:
using KNITRO
using Compat.Test

 ## Solve test problem HS15 from the Hock & Schittkowski collection.
 #
 #  min   100 (x2 - x1^2)^2 + (1 - x1)^2
 #  s.t.  x1 x2 >= 1
 #        x1 + x2^2 >= 0
 #        x1 <= 0.5
 #
 #  The standard start point (-2, 1) usually converges to the standard
 #  minimum at (0.5, 2.0), with final objective = 306.5.
 #  Sometimes the solver converges to another local minimum
 #  at (-0.79212, -1.26243), with final objective = 360.4.
 ##

function eval_f(x::Vector{Float64})
    tmp = x[2] - x[1]^2
    100.0*tmp*tmp + (1.0-x[1])^2
end

function eval_g(x::Vector{Float64}, cons::Vector{Float64})
    cons[1] = x[1] * x[2]
    cons[2] = x[1] + x[2]^2
end

function eval_grad_f(x::Vector{Float64}, grad::Vector{Float64})
    tmp = x[2] - x[1]^2
    grad[1] = (-400.0 * tmp * x[1]) - (2.0 * (1.0 - x[1]))
    grad[2] = 200.0 * tmp
end

function eval_jac_g(x::Vector{Float64}, jac::Vector{Float64})
    jac[1] = x[2]
    jac[2] = x[1]
    jac[3] = 1.0
    jac[4] = 2.0 * x[2]
end

function eval_h(x::Vector{Float64}, lambda::Vector{Float64},
                sigma::Float64, hess::Vector{Float64})
    hess[1] = sigma * ( (-400.0 * x[2]) + (1200.0 * x[1]*x[1]) + 2.0)
    hess[2] = (sigma * (-400.0 * x[1])) + lambda[1]
    hess[3] = (sigma * 200.0) + (lambda[2] * 2.0)
end

function eval_hv(x::Vector{Float64}, lambda::Vector{Float64},
                 sigma::Float64, hv::Vector{Float64})
    # H[0,0]*v[0] + H[0,1]*v[1]
    hv1 = (sigma*(((-400.0*x[2])+(1200.0*x[1]*x[1])+2.0)))*hv[1]+(sigma*(-400.0*x[1])+lambda[1])*hv[2]
    # H[1,0]*v[0] + H[1,1]*v[1]
    hv2 = (sigma*(-400.0*x[1])+lambda[1])*hv[1]+(sigma*200.0+(lambda[2]*2.0))*hv[2]
    hv[1] = hv1
    hv[2] = hv2
end

objGoal = KTR_OBJGOAL_MINIMIZE
objType = KTR_OBJTYPE_GENERAL

n = 2
x_L = [-KTR_INFBOUND, -KTR_INFBOUND]
x_U = [0.5, KTR_INFBOUND]
x = [-2.0,1.0] # initial guess

m = 2
c_Type = [KTR_CONTYPE_QUADRATIC, KTR_CONTYPE_QUADRATIC]
c_L = [1.0, 0.0]
c_U = [KTR_INFBOUND, KTR_INFBOUND]

jac_con = Int32[0,0,1,1]
jac_var = Int32[0,1,0,1]
hess_row = Int32[0,0,1]
hess_col = Int32[0,1,1]

kp = createProblem()
@test applicationReturnStatus(kp) == :Uninitialized
loadOptionsFile(kp, joinpath(dirname(@__FILE__),"tuner-fixed.opt"))
setOption(kp, KTR_PARAM_TUNER, KTR_TUNER_ON)
loadTunerFile(kp, joinpath(dirname(@__FILE__), "tuner-explore.opt"))

initializeProblem(kp, objGoal, objType, x_L, x_U, c_Type, c_L, c_U,
                  jac_var, jac_con, hess_row, hess_col; initial_x = x)
@test applicationReturnStatus(kp) == :Initialized
setCallbacks(kp, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, eval_hv)
solveProblem(kp)

# --- test optimal solutions ---
@testset "Test optimal solutions" begin
    @test applicationReturnStatus(kp) == :Optimal
    @test (abs(kp.x[1]-0.5) < 1e-4 || abs(kp.x[1]+0.79212) < 1e-4)
    @test (abs(kp.x[2]-2.0) < 1e-4 || abs(kp.x[2]+1.26243) < 1e-4)
    @test (abs(kp.obj_val[1]-306.5) < 0.025 || abs(kp.obj_val[1]-360.4) < 0.025)
end

freeProblem(kp)