text:
function nk_glfw3_device_create()
    vertex_shader = read(joinpath(@__DIR__, "vertex.glsl"), String)
    fragment_shader = read(joinpath(@__DIR__, "fragment.glsl"), String)
    glfw.ogl = nk_glfw_device()
    dev = glfw.ogl

    glfw.ogl.tex_ids = Vector{Cuint}()
    glfw.ogl.tex_count = 0

    dev.cmds = nk_create_buffer()
    nk_buffer_init_default(dev.cmds)
    dev.null = nk_draw_null_texture(nk_handle(0), nk_vec2(0,0))
    dev.prog = glCreateProgram()
    dev.vert_shdr = glCreateShader(GL_VERTEX_SHADER)
    dev.frag_shdr = glCreateShader(GL_FRAGMENT_SHADER)
    glShaderSource(dev.vert_shdr, 1, Ptr{GLchar}[pointer(vertex_shader)], C_NULL)
    glShaderSource(dev.frag_shdr, 1, Ptr{GLchar}[pointer(fragment_shader)], C_NULL)
    glCompileShader(dev.vert_shdr)
    glCompileShader(dev.frag_shdr)

    status_ref = Ref{GLint}(-1)
    glGetShaderiv(dev.vert_shdr, GL_COMPILE_STATUS, status_ref)
    if status_ref[] != GL_TRUE
        max_length_ref = Ref{GLsizei}(0)
        glGetShaderiv(dev.vert_shdr, GL_INFO_LOG_LENGTH, max_length_ref)
        actual_length_ref = Ref{GLsizei}(0)
        log = Vector{GLchar}(undef, max_length_ref[])
        glGetShaderInfoLog(dev.vert_shdr, max_length_ref[], actual_length_ref, log)
        @error "[GL]: failed to compile vertex shader: $(dev.vert_shdr): $(String(log))"
    end

    status_ref = Ref{GLint}(-1)
    glGetShaderiv(dev.frag_shdr, GL_COMPILE_STATUS, status_ref)
    if status_ref[] != GL_TRUE
        max_length_ref = Ref{GLsizei}(0)
        glGetShaderiv(dev.frag_shdr, GL_INFO_LOG_LENGTH, max_length_ref)
        actual_length_ref = Ref{GLsizei}(0)
        log = Vector{GLchar}(undef, max_length_ref[])
        glGetShaderInfoLog(dev.frag_shdr, max_length_ref[], actual_length_ref, log)
        @error "[GL]: failed to compile fragment shader: $(dev.vert_shdr): $(String(log))"
    end

    glAttachShader(dev.prog, dev.vert_shdr)
    glAttachShader(dev.prog, dev.frag_shdr)
    glLinkProgram(dev.prog)
    status_ref = Ref{GLint}(-1)
    glGetProgramiv(dev.prog, GL_LINK_STATUS, status_ref)
    @assert status_ref[] == GL_TRUE
    glDetachShader(dev.prog, dev.vert_shdr)
    glDetachShader(dev.prog, dev.frag_shdr)
    glDeleteShader(dev.vert_shdr)
    glDeleteShader(dev.frag_shdr)

    dev.uniform_tex = glGetUniformLocation(dev.prog, "tex")
    dev.uniform_proj = glGetUniformLocation(dev.prog, "proj_mat")
    dev.attrib_pos = glGetAttribLocation(dev.prog, "position")
    dev.attrib_uv = glGetAttribLocation(dev.prog, "texcoord")
    dev.attrib_col = glGetAttribLocation(dev.prog, "color")

    # buffer setup
    vs = sizeof(nk_glfw_vertex)
    vp = fieldoffset(nk_glfw_vertex, 1)
    vt = fieldoffset(nk_glfw_vertex, 2)
    vc = fieldoffset(nk_glfw_vertex, 3)

    vbo_ref = Ref{GLuint}(0)
    ebo_ref = Ref{GLuint}(0)
    vao_ref = Ref{GLuint}(0)
    glGenBuffers(1, vbo_ref)
    glGenBuffers(1, ebo_ref)
    glGenVertexArrays(1, vao_ref)
    dev.vao = vao_ref[]
    dev.vbo = vbo_ref[]
    dev.ebo = ebo_ref[]

    glBindVertexArray(dev.vao)
    glBindBuffer(GL_ARRAY_BUFFER, dev.vbo)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, dev.ebo)

    glEnableVertexAttribArray(dev.attrib_pos)
    glEnableVertexAttribArray(dev.attrib_uv)
    glEnableVertexAttribArray(dev.attrib_col)

    glVertexAttribPointer(dev.attrib_pos, 2, GL_FLOAT, GL_FALSE, vs, Ptr{Cvoid}(vp))
    glVertexAttribPointer(dev.attrib_uv, 2, GL_FLOAT, GL_FALSE, vs, Ptr{Cvoid}(vt))
    glVertexAttribPointer(dev.attrib_col, 4, GL_UNSIGNED_BYTE, GL_TRUE, vs, Ptr{Cvoid}(vc))

    glBindTexture(GL_TEXTURE_2D, 0)
    glBindBuffer(GL_ARRAY_BUFFER, 0)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0)
    glBindVertexArray(0)

    return dev
end

function nk_glfw3_init(win::GLFW.Window, init_state, max_vertex_buffer, max_element_buffer)
    glfw.win = win
    if init_state == NK_GLFW3_INSTALL_CALLBACKS
        GLFW.SetScrollCallback(win, nk_gflw3_scroll_callback)
        GLFW.SetCharCallback(win, nk_glfw3_char_callback)
        GLFW.SetMouseButtonCallback(win, nk_glfw3_mouse_button_callback)
    end
    glfw.ctx = nk_create_context()
    nk_init_default(glfw.ctx, C_NULL)
    clipboard_copy = @cfunction(nk_glfw3_clipboard_copy, Cvoid, (nk_handle, Ptr{UInt8}, Cint))
    clipboard_paste = @cfunction(nk_glfw3_clipboard_paste, Cvoid, (nk_handle, Ptr{nk_text_edit}))
    nk_set_clip(glfw.ctx, clipboard_copy, clipboard_paste)
    glfw.last_button_click = 0
    nk_glfw3_device_create()

    glfw.is_double_click_down = false
    glfw.double_click_pos = nk_vec2(0, 0)
    glfw.text_len = 0
    glfw.text = Vector{Char}(undef, NK_GLFW_TEXT_MAX)
    return glfw.ctx
end

function nk_glfw3_device_upload_atlas(image, width, height)
    dev = glfw.ogl
    tex = Ref{GLuint}(0)
    glGenTextures(1, tex)
    dev.font_tex = tex[]
    glBindTexture(GL_TEXTURE_2D, dev.font_tex)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, GLsizei(width), GLsizei(height), 0,
                 GL_RGBA, GL_UNSIGNED_BYTE, image)
end

function nk_glfw3_font_stash_begin()
    atlas = nk_create_font_atlas()
    GC.@preserve atlas begin
        nk_font_atlas_init_default(atlas)
        nk_font_atlas_begin(atlas)
        glfw.atlas = atlas
    end
    return nothing
end

function nk_glfw3_font_stash_end()
    w = Ref{Cint}(0)
    h = Ref{Cint}(0)
    image = nk_font_atlas_bake(glfw.atlas, w, h, NK_FONT_ATLAS_RGBA32)
    nk_glfw3_device_upload_atlas(image, w[], h[])
    null_ref = Ref{nk_draw_null_texture}(glfw.ogl.null)
    nk_font_atlas_end(glfw.atlas, nk_handle_id(glfw.ogl.font_tex), null_ref)
    glfw.ogl.null = null_ref[]
    nk_set_default_font(glfw.ctx, glfw.atlas)
    return nothing
end

function nk_glfw3_render(anti_alias, max_vertex_buffer::Integer, max_element_buffer::Integer)
    dev = glfw.ogl
    orthomatrix = GLfloat[2.0, 0.0, 0.0, 0.0,
                          0.0,-2.0, 0.0, 0.0,
                          0.0, 0.0,-1.0, 0.0,
                         -1.0, 1.0, 0.0, 1.0]
    orthomatrix[1] /= GLfloat(glfw.width)
    orthomatrix[6] /= GLfloat(glfw.height)

    glEnable(GL_BLEND)
    glBlendEquation(GL_FUNC_ADD)
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
    glDisable(GL_CULL_FACE)
    glDisable(GL_DEPTH_TEST)
    glEnable(GL_SCISSOR_TEST)
    glActiveTexture(GL_TEXTURE0)

    glUseProgram(dev.prog)
    glUniform1i(dev.uniform_tex, 0)
    glUniformMatrix4fv(dev.uniform_proj, 1, GL_FALSE, orthomatrix)
    glViewport(0, 0, glfw.display_width, glfw.display_height)

    # convert from command queue into draw list and draw to screen
    # allocate vertex and element buffer
    glBindVertexArray(dev.vao)
    glBindBuffer(GL_ARRAY_BUFFER, dev.vbo)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, dev.ebo)

    glBufferData(GL_ARRAY_BUFFER, max_vertex_buffer, C_NULL, GL_STREAM_DRAW)
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, max_element_buffer, C_NULL, GL_STREAM_DRAW)

    # load draw vertices & elements directly into vertex + element buffer
    vertices = glMapBuffer(GL_ARRAY_BUFFER, GL_WRITE_ONLY)
    elements = glMapBuffer(GL_ELEMENT_ARRAY_BUFFER, GL_WRITE_ONLY)

    # fill convert configuration
    cfg = nk_set_convert_config(anti_alias, dev.null)

    # setup buffers to load vertices and elements
    vbo_buffer = nk_create_buffer()
    ebo_buffer = nk_create_buffer()
    GC.@preserve vbo_buffer ebo_buffer vertices elements cfg begin
        nk_buffer_init_fixed(vbo_buffer, vertices, max_vertex_buffer)
        nk_buffer_init_fixed(ebo_buffer, elements, max_element_buffer)
        res = nk_convert(glfw.ctx, dev.cmds, vbo_buffer, ebo_buffer, cfg)
        @assert res == NK_CONVERT_SUCCESS "convert failed!"
    end

    glUnmapBuffer(GL_ARRAY_BUFFER)
    glUnmapBuffer(GL_ELEMENT_ARRAY_BUFFER)

    # iterate over and execute each draw command
    offset = Ptr{nk_draw_index}(0)
    draw_cmd = nk__draw_begin(glfw.ctx, dev.cmds)
    while draw_cmd != C_NULL
        elem_count_ref = Ref{UInt32}()
        clip_rect_ref = Ref{nk_rect}()
        texture_ref = Ref{nk_handle}()
        nk_unpack_draw_command(draw_cmd, elem_count_ref, clip_rect_ref, texture_ref)
        elem_count = elem_count_ref[]
        clip_rect = clip_rect_ref[]
        texture = texture_ref[]
        if elem_count == 0
            draw_cmd = nk__draw_next(draw_cmd, dev.cmds, glfw.ctx)
            continue
        end
        glBindTexture(GL_TEXTURE_2D, texture.id)
        glScissor(
            GLint(clip_rect.x * glfw.fb_scale.x),
            GLint((glfw.height - (clip_rect.y + clip_rect.h)) * glfw.fb_scale.y),
            GLint(clip_rect.w * glfw.fb_scale.x),
            GLint(clip_rect.h * glfw.fb_scale.y))
        glDrawElements(GL_TRIANGLES, GLsizei(elem_count), GL_UNSIGNED_SHORT, offset)
        offset += elem_count * Core.sizeof(nk_draw_index)
        draw_cmd = nk__draw_next(draw_cmd, dev.cmds, glfw.ctx)
    end
    nk_clear(glfw.ctx)
    nk_buffer_free(vbo_buffer)
    nk_buffer_free(ebo_buffer)
    # reset OpenGL state
    glUseProgram(0)
    glBindBuffer(GL_ARRAY_BUFFER, 0)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0)
    glBindVertexArray(0)
    glDisable(GL_BLEND)
    glDisable(GL_SCISSOR_TEST)
end

function nk_glfw3_new_frame()
    ctx = glfw.ctx
    win = glfw.win

    glfw.width, glfw.height = GLFW.GetWindowSize(win)
    glfw.display_width, glfw.display_height = GLFW.GetFramebufferSize(win)
    glfw.fb_scale = nk_vec2(glfw.display_width/glfw.width, glfw.display_height/glfw.height)

    nk_input_begin(ctx)
    i = 1
    while i <= glfw.text_len
        nk_input_unicode(ctx, glfw.text[i])
        i += 1
    end

    nk_input_key(ctx, NK_KEY_DEL, GLFW.GetKey(win, GLFW.KEY_DELETE))
    nk_input_key(ctx, NK_KEY_ENTER, GLFW.GetKey(win, GLFW.KEY_ENTER))
    nk_input_key(ctx, NK_KEY_TAB, GLFW.GetKey(win, GLFW.KEY_TAB))
    nk_input_key(ctx, NK_KEY_BACKSPACE, GLFW.GetKey(win, GLFW.KEY_BACKSPACE))
    nk_input_key(ctx, NK_KEY_UP, GLFW.GetKey(win, GLFW.KEY_UP))
    nk_input_key(ctx, NK_KEY_DOWN, GLFW.GetKey(win, GLFW.KEY_DOWN))
    nk_input_key(ctx, NK_KEY_TEXT_START, GLFW.GetKey(win, GLFW.KEY_HOME))
    nk_input_key(ctx, NK_KEY_TEXT_END, GLFW.GetKey(win, GLFW.KEY_END))
    nk_input_key(ctx, NK_KEY_SCROLL_START, GLFW.GetKey(win, GLFW.KEY_HOME))
    nk_input_key(ctx, NK_KEY_SCROLL_END, GLFW.GetKey(win, GLFW.KEY_END))
    nk_input_key(ctx, NK_KEY_SCROLL_DOWN, GLFW.GetKey(win, GLFW.KEY_PAGE_DOWN))
    nk_input_key(ctx, NK_KEY_SCROLL_UP, GLFW.GetKey(win, GLFW.KEY_PAGE_UP))
    nk_input_key(ctx, NK_KEY_SHIFT, GLFW.GetKey(win, GLFW.KEY_LEFT_SHIFT)  ||
                                    GLFW.GetKey(win, GLFW.KEY_RIGHT_SHIFT))

    if GLFW.GetKey(win, GLFW.KEY_LEFT_CONTROL) || GLFW.GetKey(win, GLFW.KEY_RIGHT_CONTROL)
        nk_input_key(ctx, NK_KEY_COPY, GLFW.GetKey(win, GLFW.KEY_C))
        nk_input_key(ctx, NK_KEY_PASTE, GLFW.GetKey(win, GLFW.KEY_V))
        nk_input_key(ctx, NK_KEY_CUT, GLFW.GetKey(win, GLFW.KEY_X))
        nk_input_key(ctx, NK_KEY_TEXT_UNDO, GLFW.GetKey(win, GLFW.KEY_Z))
        nk_input_key(ctx, NK_KEY_TEXT_REDO, GLFW.GetKey(win, GLFW.KEY_R))
        nk_input_key(ctx, NK_KEY_TEXT_WORD_LEFT, GLFW.GetKey(win, GLFW.KEY_LEFT))
        nk_input_key(ctx, NK_KEY_TEXT_WORD_RIGHT, GLFW.GetKey(win, GLFW.KEY_RIGHT))
        nk_input_key(ctx, NK_KEY_TEXT_LINE_START, GLFW.GetKey(win, GLFW.KEY_B))
        nk_input_key(ctx, NK_KEY_TEXT_LINE_END, GLFW.GetKey(win, GLFW.KEY_E))
    else
        nk_input_key(ctx, NK_KEY_LEFT, GLFW.GetKey(win, GLFW.KEY_LEFT))
        nk_input_key(ctx, NK_KEY_RIGHT, GLFW.GetKey(win, GLFW.KEY_RIGHT))
        nk_input_key(ctx, NK_KEY_COPY, 0)
        nk_input_key(ctx, NK_KEY_PASTE, 0)
        nk_input_key(ctx, NK_KEY_CUT, 0)
        nk_input_key(ctx, NK_KEY_SHIFT, 0)
    end

    x, y = GLFW.GetCursorPos(win)
    x, y = trunc(Int, x), trunc(Int, y)
    nk_input_motion(ctx, x, y)
    nk_input_button(ctx, NK_BUTTON_LEFT, x, y, GLFW.GetMouseButton(win, GLFW.MOUSE_BUTTON_LEFT))
    nk_input_button(ctx, NK_BUTTON_MIDDLE, x, y, GLFW.GetMouseButton(win, GLFW.MOUSE_BUTTON_MIDDLE))
    nk_input_button(ctx, NK_BUTTON_RIGHT, x, y, GLFW.GetMouseButton(win, GLFW.MOUSE_BUTTON_RIGHT))
    dx, dy = trunc(Int, glfw.double_click_pos.x), trunc(Int, glfw.double_click_pos.y)
    nk_input_button(ctx, NK_BUTTON_DOUBLE, dx, dy, glfw.is_double_click_down)
    nk_input_scroll(ctx, glfw.scroll)
    nk_input_end(ctx)

    glfw.text_len = 0
    glfw.scroll = nk_vec2(0,0)
end

function nk_glfw3_device_destroy()
    dev = glfw.ogl
    glDeleteProgram(dev.prog)
    font_tex_ref = Ref{UInt32}(dev.font_tex)
    vbo_ref = Ref{UInt32}(dev.vbo)
    ebo_ref = Ref{UInt32}(dev.ebo)
    glDeleteTextures(1, font_tex_ref)
    for tid in dev.tex_ids
        tid_ref = Ref{UInt32}(tid)
        glDeleteTextures(1, tid_ref)
    end
    glDeleteBuffers(1, vbo_ref)
    glDeleteBuffers(1, ebo_ref)
    nk_buffer_free(dev.cmds)
end

function nk_glfw3_shutdown()
    nk_font_atlas_clear(glfw.atlas)
    nk_free(glfw.ctx)
    nk_glfw3_device_destroy()
end

function nk_glfw3_create_texture(width, height; format=GL_RGBA, type=GL_UNSIGNED_BYTE)
    tex = Ref{GLuint}(0)
    glGenTextures(1, tex)
    push!(glfw.ogl.tex_ids, tex[])
    glfw.ogl.tex_count += 1
    glBindTexture(GL_TEXTURE_2D, tex[])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, format, GLsizei(width), GLsizei(height), 0, format, type, C_NULL)
    return glfw.ogl.tex_count
end

function nk_glfw3_update_texture(idx, image, width, height; format=GL_RGBA, type=GL_UNSIGNED_BYTE)
    glBindTexture(GL_TEXTURE_2D, glfw.ogl.tex_ids[idx])
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, GLsizei(width), GLsizei(height), format, type, image)
end

function nk_glfw3_delete_texture(idx)
    tid = glfw.ogl.tex_ids[idx]
    tid_ref = Ref{UInt32}(tid)
    glDeleteTextures(1, tid_ref)
    deleteat!(glfw.ogl.tex_ids, idx)
    glfw.ogl.tex_count -= 1
end

function create_nk_image(idx)
    tid = glfw.ogl.tex_ids[idx]
    return nk_image(nk_handle(tid), 0, 0, (0,0,0,0))
end


text:
using Documenter
using Documenter.Anchors
using Documenter.Documents
using Documenter.Expanders
using Documenter.Selectors

using Bibliography: xnames, xyear, xlink, xtitle, xin

abstract type BibliographyBlock <: Expanders.ExpanderPipeline end

Selectors.order(::Type{BibliographyBlock}) = 12.0  # Expand bibliography last
Selectors.matcher(::Type{BibliographyBlock}, node, page, doc) = Expanders.iscode(node, r"^@bibliography")

function Selectors.runner(::Type{BibliographyBlock}, x, page, doc)
    @info "Expanding bibliography."
    raw_bib = "
"
    for (id, entry) in sort(BIBLIOGRAPHY)
        @info "Expanding bibliography entry: $id."

        # Add anchor that citations can link to from anywhere in the docs.
        Anchors.add!(doc.internal.headers, entry, entry.id, page.build)

        entry_text = """

$id

          $(xnames(entry)) ($(xyear(entry))), $(xtitle(entry)), $(xin(entry))

"""
        raw_bib *= entry_text
    end
    raw_bib *= "\n
"

    page.mapping[x] = Documents.RawNode(:html, raw_bib)
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PublicIPAddressDnsSettings <: SwaggerModel
    domainNameLabel::Any # spec type: Union{ Nothing, String } # spec name: domainNameLabel
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn
    reverseFqdn::Any # spec type: Union{ Nothing, String } # spec name: reverseFqdn

    function PublicIPAddressDnsSettings(;domainNameLabel=nothing, fqdn=nothing, reverseFqdn=nothing)
        o = new()
        validate_property(PublicIPAddressDnsSettings, Symbol("domainNameLabel"), domainNameLabel)
        setfield!(o, Symbol("domainNameLabel"), domainNameLabel)
        validate_property(PublicIPAddressDnsSettings, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        validate_property(PublicIPAddressDnsSettings, Symbol("reverseFqdn"), reverseFqdn)
        setfield!(o, Symbol("reverseFqdn"), reverseFqdn)
        o
    end
end # type PublicIPAddressDnsSettings

const _property_map_PublicIPAddressDnsSettings = Dict{Symbol,Symbol}(Symbol("domainNameLabel")=>Symbol("domainNameLabel"), Symbol("fqdn")=>Symbol("fqdn"), Symbol("reverseFqdn")=>Symbol("reverseFqdn"))
const _property_types_PublicIPAddressDnsSettings = Dict{Symbol,String}(Symbol("domainNameLabel")=>"String", Symbol("fqdn")=>"String", Symbol("reverseFqdn")=>"String")
Base.propertynames(::Type{ PublicIPAddressDnsSettings }) = collect(keys(_property_map_PublicIPAddressDnsSettings))
Swagger.property_type(::Type{ PublicIPAddressDnsSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PublicIPAddressDnsSettings[name]))}
Swagger.field_name(::Type{ PublicIPAddressDnsSettings }, property_name::Symbol) =  _property_map_PublicIPAddressDnsSettings[property_name]

function check_required(o::PublicIPAddressDnsSettings)
    true
end

function validate_property(::Type{ PublicIPAddressDnsSettings }, name::Symbol, val)
end


text:
function raster_one_to_all(T, V, cfg)::Matrix{T}

    # Load the data
    rasterdata = load_raster_data(T, V, cfg)

    # Get flags
    flags = get_raster_flags(cfg)

    # Send to main kernel
    onetoall_kernel(rasterdata, flags, cfg)
end

function onetoall_kernel(data::RasData{T,V}, flags, cfg)::Matrix{T} where {T,V}

    # Data
    strengths = data.strengths
    included_pairs = data.included_pairs
    points_rc = data.points_rc
    gmap = data.cellmap
    polymap = data.polymap
    hbmeta = data.hbmeta
    source_map = data.source_map

    # Flags
    use_variable_strengths = !isempty(strengths)
    use_included_pairs = !isempty(included_pairs)
    mode = included_pairs.mode == :include ? 0 : 1
    one_to_all = flags.is_onetoall
    avg_res = flags.avg_res
    four_neighbors = flags.four_neighbors

    if use_included_pairs
        points_unique = included_pairs.point_ids
        prune_points!(points_rc, included_pairs.point_ids)
        if use_variable_strengths
            prune_strengths!(strengths, included_pairs.point_ids)
        end
    end

    # Construct point map
    point_map = zeros(V, size(gmap))
    f(i, x) = points_rc[i][x]
    for x = 1:size(points_rc[1], 1)
        point_map[f(1,x), f(2,x)] = f(3, x)
    end

    points_unique = unique(points_rc[3])

    newpoly = create_new_polymap(gmap, polymap, points_rc, 0, 0, point_map)

    nodemap = construct_node_map(gmap, newpoly)

    a = construct_graph(gmap, nodemap, avg_res, four_neighbors)
    cc = connected_components(SimpleWeightedGraph(a))
    G = laplacian(a)
    csinfo("There are $(size(a, 1)) points and $(length(cc)) connected components")

    # source_map = Matrix{eltype(a)}(0, 0)
    # ground_map = Matrix{eltype(a)}(0, 0)
    s = zeros(eltype(a), size(point_map))
    z = deepcopy(s)
    cum = initialize_cum_maps(gmap, flags.outputflags.write_max_cur_maps)

    point_ids = included_pairs.point_ids
    res = zeros(eltype(a), size(points_unique, 1)) |> SharedArray
    num_points_to_solve = size(points_unique, 1)
    original_point_map = copy(point_map)
    unique_point_map = zeros(V, size(gmap))

    for i in points_unique
        ind = findfirst(x -> x == i, points_rc[3])
        unique_point_map[f(1,ind), f(2,ind)] = f(3,ind)
    end

    # @distributed for i = 1:num_points_to_solve
    function f(i)
        # copyto!(point_map, original_point_map)
        point_map = copy(original_point_map)
        str = use_variable_strengths ? strengths[i,2] : 1
        csinfo("Solving point $i of $num_points_to_solve")
        # copyto!(s, z)
        s = copy(z)
        n = points_unique[i]
        if use_included_pairs
            for j = 1:size(point_ids,1)
                if i != j && included_pairs.include_pairs[i,j] == mode
                    exclude = point_ids[j]
                    map!(x -> x == exclude ? 0 : x, point_map, point_map)
                end
            end
            # polymap = create_new_polymap(gmap, Polymap(polymap), points_rc, point_map = point_map)
            newpoly = create_new_polymap(gmap, polymap, points_rc, 0, 0, point_map)
            nodemap = construct_node_map(gmap, polymap)
            a = construct_graph(gmap, nodemap, avg_res, four_neighbors)
        end
        # T = eltype(a)
        if one_to_all
            #source_map = map(x -> x == n ? str : 0, point_map)
            source_map = map(x -> x == n ? T(str) : T(0), unique_point_map)
            ground_map = map(x -> x == n ? T(0) : T(x), point_map)
            map!(x -> x > 0 ? Inf : x, ground_map, ground_map)
        else
            source_map = map(x -> x != 0 ? T(x) : T(0), point_map)
            map!(x -> x == n ? 0 : x, source_map, source_map)
            map!(x -> x != 0 ? 1 : x, source_map, source_map)
            ground_map = map(x -> x == n ? Inf : T(0), point_map)
        end

        check_node = nodemap[points_rc[1][i], points_rc[2][i]]

        policy = one_to_all ? :rmvgnd : :rmvsrc
        sources, grounds, finite_grounds =
                    _get_sources_and_grounds(source_map, ground_map,
                            flags, G, nodemap, policy)

        advanced_data = AdvancedData(G, cc, nodemap, newpoly, hbmeta,
                        sources, grounds, source_map, finite_grounds,
                        check_node, n, gmap)

        if one_to_all
            # v = advanced(cfg, a, source_map, ground_map; nodemap = nodemap, policy = :rmvgnd,
            #                check_node = check_node, src = n, polymap = Polymap(newpoly), hbmeta = hbmeta)
            v, curr = advanced_kernel(advanced_data, flags, cfg)
        else
            # v = advanced(cfg, a, source_map, ground_map; nodemap = nodemap, policy = :rmvsrc,
            #                check_node = check_node, src = n, polymap = Polymap(newpoly), hbmeta = hbmeta)
            v, curr = advanced_kernel(advanced_data, flags, cfg)
        end
        res[i] = v[1]

        cum.cum_curr[mycsid()] .+= curr
        flags.outputflags.write_max_cur_maps && (cum.max_curr[mycsid()] .= max.(cum.max_curr[mycsid()], curr))
    end

    pmap(x -> f(x), 1:num_points_to_solve)

    if flags.outputflags.write_cur_maps
        write_cum_maps(cum, gmap, cfg, hbmeta,
                       flags.outputflags.write_max_cur_maps,
                       flags.outputflags.write_cum_cur_map_only)
    end

    hcat(points_unique, res)
end

function prune_points!(points_rc, point_ids::Vector{V}) where V
    rmv = V[]
    for (i,p) in enumerate(points_rc[3])
        if p in point_ids
            continue
        else
            #for it in 1:3 deleteat!(points_rc[it], i) end
            push!(rmv, i)
        end
    end
    for i in 1:3 deleteat!(points_rc[i], rmv) end
end

function prune_strengths!(strengths, point_ids::Vector{V}) where V
    pts = strengths[:,1]
    l = length(pts)
    rmv = V[]
    for (i,p) in enumerate(pts)
        if !(p in point_ids)
           push!(rmv, i)
       end
    end
    rng = collect(1:l)
    deleteat!(rng, rmv)
    strengths[rng,:]
end


text:
using Soss
using AbstractGPs
using KernelFunctions

x = randn(100)
y = sinpi.(x) + 0.3*randn(100)

import TransformVariables
const TV = TransformVariables
Soss.as(gp::AbstractGPs.FiniteGP, _data=NamedTuple()) = as(Array, asℝ,size(gp)...)

sqexpkernel(alpha::Real, rho::Real) =
alpha^2 * transform(SEKernel(), 1/(rho*sqrt(2)))
using LinearAlgebra
using PositiveFactorizations

# Define model.
m = @model x begin
    # Priors.
    alpha ~ LogNormal(0.0, 0.1)
    rho ~ LogNormal(0.0, 1.0)
    sigma ~ LogNormal(0.0, 1.0)

    # Covariance function.
    kernel = sqexpkernel(alpha, rho)

    # GP (implicit zero-mean).
    gp = GP(kernel)

    y ~ gp(x, 0.6+sigma^2)
    # Sampling Distribution (MvNormal likelihood).
    # y ~ gp(x, sigma^2 + 1e-6)  # add 1e-6 for numerical stability.
end;

PositiveFactorizations.floattype(::Type{T}) where {T} = T

post = dynamicHMC(m(x=x), (;y=y); ad_backend=Val(:Zygote));
particles(post)

gp1 = predict(m(x=x), post[1]).gp(x, post[1].sigma^2 + 1e-6)

gp1post = posterior(gp1, y)

using Plots
plt = scatter(x, y; label = "Data")
plot!(plt, gp1post, range(extrema(x)...; length=100); label="Posterior")

plot(gp1post)


text:
include("../src/NMRSpectraSimulator.jl")
import .NMRSpectraSimulator

using LinearAlgebra
using FFTW
import PyPlot
#import BSON
#import JLD

#import Clustering
import Statistics

import Random
Random.seed!(25)

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

### user inputs.

tol_coherence = 1e-2
α_relative_threshold = 0.05
λ0 = 3.4
Δcs_max = 0.2
κ_λ_lb = 0.5
κ_λ_ub = 2.5

#molecule_names = ["L-Serine"; "L-Phenylalanine"; "DSS";]
molecule_names = ["D-(+)-Glucose"; "DSS"]

# machine values taken from the BMRB 700 MHz 20 mM glucose experiment.
fs = 14005.602240896402
SW = 20.0041938620844
ν_0ppm = 10656.011933076665

# # machine values for the BMRB 500 MHz glucose experiment.
# ν_0ppm = 6752.490995937095
# SW = 16.0196917451925
# fs = 9615.38461538462

# path to the GISSMO Julia storage folder.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"

### end inputs.

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

Δcs_max_mixture = collect( Δcs_max for i = 1:length(molecule_names))

dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0)
mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
    H_params_path, ppm2hzfunc, fs, SW, λ0,
    ν_0ppm, dummy_SSFID;
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold)
As = mixture_params

u_min = ppm2hzfunc(-0.5)
u_max = ppm2hzfunc(4.0)

NMRSpectraSimulator.fitproxies!(As;
    κ_λ_lb = κ_λ_lb,
    κ_λ_ub = κ_λ_ub,
    u_min = u_min,
    u_max = u_max,
    Δr = 1.0,
    Δκ_λ = 0.05)

### plot.

function convertΔcstoΔω0(x::T, fs::T, SW::T)::T where T
    return x*2*π*fs/SW
end

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)
U_rad = U .* (2*π)

# purposely distort the spectra by setting non-autophased FID values.
Ag = As[end]
Ag.ss_params.d[1] = convertΔcstoΔω0(1.0*0.05, fs, SW) # increase by 0.05 ppm
Ag.ss_params.κs_λ = rand(length(Ag.ss_params.κs_λ)) .+ 1
Ag.ss_params.κs_β[:] = collect( rand(length(Ag.ss_params.κs_β[i])) .* (2*π) for i = 1:length(Ag.ss_params.κs_β) )

f = uu->NMRSpectraSimulator.evalmixture(uu, mixture_params)
f_U = f.(U_rad)

# test params.
ΩS_ppm = collect( hz2ppmfunc.( NMRSpectraSimulator.combinevectors(A.Ωs) ./ (2*π) ) for A in mixture_params )
#ΩS_ppm_flat = NMRSpectraSimulator.combinevectors(ΩS_ppm)

As2 = collect( NMRSpectraSimulator.κCompoundFIDType(As[i]) for i = 1:length(As) )

# purposely modify As2, DSS element.
Eg = As2[end]
#Eg.ss_params.κ = collect( rand(length(Eg.κ[i])) for i = 1:length(Eg.κ) )
#Eg.κs_α_singlets = rand(length(Eg.κs_α_singlets))
Eg.κs_α[1][1] = 0.3
Eg.κs_α[1][2] = 0.7
Eg.κs_α[1][3] = 0.1
Eg.κs_α_singlets[1] = 0.4

fill!(Ag.ss_params.d, 0.0) # shift back.

#@assert 1==4

## parameters that affect qs.
# A.d, A.κs_λ, A.κs_β
# A.d_singlets, A.αs_singlets, A.Ωs_singlets, A.β_singlets, A.λ0, A.κs_λ_singlets

#q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, mixture_params)
q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As2)

q_U = q.(U_rad)

## visualize.
PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P, real.(f_U), label = "f")
PyPlot.plot(P, real.(q_U), label = "q")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("f vs q")


text:
__precompile__(false)

import Base.print, Base.println

using TimeManager
using LoggerManager
using ThreadManager

using Distributed

######################################################
# Overwrite Base functions

PrintMutex = Threads.Mutex()
#ThreadManager.setPrintMutex(PrintMutex)

print(args...) = thread_call(() -> LoggerManager.print(args...);mutex=PrintMutex)
println(args...) = thread_call(() -> LoggerManager.println(args...);mutex=PrintMutex)

PushMutex = Threads.Mutex()
thread_push!(args...) = thread_call(() -> Base.push!(args...);mutex=PushMutex)

sleep(sec) = thread_sleep(sec)

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

info(args...) = thread_call(() -> LoggerManager.info(args...);mutex=PrintMutex)
debug(args...) = thread_call(() -> LoggerManager.debug(args...);mutex=PrintMutex)
warn(args...) = thread_call(() -> LoggerManager.warn(args...);mutex=PrintMutex)
error(args...) = thread_call(() -> LoggerManager.error(args...);mutex=PrintMutex)

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

Messages = String[] #shared object among threads
function message(this::Thread, args... ;mode="", title="", lineBreak=true)
  global Messages
  thread_push!(Messages, LoggerManager.msg(args... ;time=true, name=thread_id()*":"*this.name, mode=mode, title=title, lineBreak=lineBreak))
end

function showMessages()
  global Messages
  msgs = String[]

  # copy & clear messages
  thread_call(() -> begin msgs = deepcopy(Messages); Messages = Array{String,1}() end ;mutex=PushMutex)

  thread_call(function()
    if length(msgs) > 0
      LoggerManager.println("---[Messages]---")
      for msg in msgs LoggerManager.print(msg) end
      LoggerManager.println("----------------")
    end
  end ;mutex=PrintMutex)
end

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

print(this::Thread, args... ;title="") = message(this, args... ;title=title, lineBreak=false)
println(this::Thread, args... ;title="") = message(this, args... ;title=title)

info(this::Thread, args... ;title="") = message(this, args... ;mode=:Info, title=title)
debug(this::Thread, args... ;title="") = message(this, args... ;mode=:Debug, title=title)
warn(this::Thread, args... ;title="") = message(this, args... ;mode=:Warning, title=title)
error(this::Thread, args... ;title="") = message(this, args... ;mode=:Error, title=title)

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Make input file for a subject, block, and epi
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_data_to_input(Exp_data_frame::DataFrames.DataFrame;
                            Sub=1, Epi=1, Block=2)
    selected_data = @where(Exp_data_frame,
                           :env .== Block+2, :subID .== Sub, :epi .== Epi)

    obs = selected_data[:state] .- 1
    act = selected_data[:action] .- 1

    if obs[end] != 0
        obs = cat(obs, 0; dims=1)
        act = cat(act, -1; dims=1)
    end

    return Str_Input(; obs=obs, act=act)
end

function Func_data_to_input(Exp_data_frame::Array{Str_Input,3};
                            Sub=1, Epi=1, Block=2)
    return Exp_data_frame[Sub,Block,Epi]
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Initialize an Agent for a specific input
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_initial_Agent(; Sub=1, Epi=1, Block=1,
                             Param=Str_Param(),
                             HyperParam=Str_HyperParam(),
                             Input=Str_Input(),
                             Agent_Prev=Str_Agent())

    First_Epi = (Epi == 1)

    if First_Epi&&(Block==1)
        # Belief
        Π = ones(Param.num_obs, Param.num_obs*Param.num_act, 1)
        Π *= Param.ϵ

        # Counter
        Counter = zeros(Param.num_obs,1)
        Counter[ Input.obs[1] + 1 ] += 1

        # Novelty
        if HyperParam.novelty_form == "log-count"
            P_obs = (Counter .+ 1) ./ (sum( Counter .+ 1) )
            Novelty = - log.(P_obs)

        elseif HyperParam.novelty_form == "inv-count"
            P_obs = (Counter .+ 1) ./ (sum( Counter .+ 1) )
            Novelty = 1 ./ P_obs

        elseif HyperParam.novelty_form == "trap-det"
            Novelty = -1 .* (Counter .> Param.C_thresh)

        elseif HyperParam.novelty_form == "trap-det-2"
            Count = Counter
            Novelty = zeros(size(Count))
            ind_Count = sortperm(Count[:],rev=true)[1:Int(Param.C_thresh)]
            Novelty[ind_Count[Count[ind_Count] .> 0]] .= -1

        else
            error("Novelty format is not known!")
        end

        Q_R_MF = ones(Param.num_obs*Param.num_act, 1) .* Param.Q_R_0
        Q_N_MF = ones(Param.num_obs*Param.num_act, 1) .* Param.Q_N_0
        E_R = zeros(Param.num_obs*Param.num_act, 1)
        E_N = zeros(Param.num_obs*Param.num_act, 1)

        # QR and UR
        if HyperParam.First_Epi_reward_update
            Ur0 = 1 / (1 - Param.γ_RL_R)
        else
            Ur0 = 0
        end
        Q_R_MB = ones(Param.num_obs*Param.num_act, 1) .* Ur0
        U_R = ones(Param.num_obs, 1) .* Ur0
        # QN and UN
        if HyperParam.exploration_form == "Novelty"
            Un00 = mean(Novelty)
        elseif HyperParam.exploration_form == "MB_Curiosity"
            Un00 = log(Param.num_obs)
        end
        if Param.γ_RL_N_MB == -1
            Un0 = Un00 / (1 - Param.γ_RL_N)
        else
            Un0 = Un00 / (1 - Param.γ_RL_N_MB)
        end
        Q_N_MB = ones(Param.num_obs*Param.num_act, 1) .* Un0
        U_N = ones(Param.num_obs, 1) .* Un0

        State = ones(1)
        State *= Input.obs[1]*Param.num_act + Input.act[1]
        P_Action = ones(Param.num_act,1) ./ Param.num_act

        t_since_change = ones(1)
    else
        # Belief
        Π = Agent_Prev.Π[:,:,end:end]

        # Counter
        Counter = Agent_Prev.Counter[:,end:end] .* Param.count_decay
        Counter[ Input.obs[1] + 1 , 1] += 1

        # Novelty
        if HyperParam.novelty_form == "log-count"
            P_obs = (Counter .+ 1) ./ (sum( Counter .+ 1) )
            Novelty = - log.(P_obs)

        elseif HyperParam.novelty_form == "inv-count"
            P_obs = (Counter .+ 1) ./ (sum( Counter .+ 1) )
            Novelty = 1 ./ P_obs

        elseif HyperParam.novelty_form == "trap-det"
            Novelty = -1 .* (Counter .> Param.C_thresh)

        elseif HyperParam.novelty_form == "trap-det-2"
            Count = Counter
            Novelty = zeros(size(Count))
            ind_Count = sortperm(Count[:],rev=true)[1:Int(Param.C_thresh)]
            Novelty[ind_Count[Count[ind_Count] .> 0]] .= -1

        else
            error("Novelty format is not known!")
        end

        Q_R_MF = Agent_Prev.Q_R_MF[:,end:end]
        Q_N_MF = Agent_Prev.Q_N_MF[:,end:end]
        E_R = zeros(Param.num_obs*Param.num_act, 1)
        E_N = zeros(Param.num_obs*Param.num_act, 1)

        Q_R_MB = Agent_Prev.Q_R_MB[:,end:end]
        Q_N_MB = Agent_Prev.Q_N_MB[:,end:end]
        U_R = Agent_Prev.U_R[:,end:end]
        U_N = Agent_Prev.U_N[:,end:end]

        State = [Input.obs[1]*Param.num_act + Input.act[1]]

        # Model based updating
        if First_Epi
            Q_temp_MB = Q_R_MB .+ ( Param.β_ratio_N2R .* Q_N_MB )
            Q_temp_MF = Q_R_MF .+ ( Param.β_ratio_N2R .* Q_N_MF )

            Q_temp = (Param.w_MF_to_MB_1 * Param.w_MF_to_MB_Scaler) .* Q_temp_MF .+
                     (1-Param.w_MF_to_MB_1) .* Q_temp_MB
        else
            Q_temp_MB = Q_R_MB
            Q_temp_MF = Q_R_MF

            Q_temp = (Param.w_MF_to_MB_2 * Param.w_MF_to_MB_Scaler) .* Q_temp_MF .+
                     (1-Param.w_MF_to_MB_2) .* Q_temp_MB
        end

        Q_temp = Q_temp[ (Input.obs[1]*Param.num_act) .+ (1:Param.num_act), 1:1]
        Q_temp = (Q_temp .- findmin(Q_temp)[1]) .* Param.β_R

        P_Action = exp.(Q_temp)
        P_Action = P_Action ./ sum(P_Action)

        t_since_change = ones(1) .* Agent_Prev.t_since_change[end]
    end

    Str_Agent(; Sub=Sub, Epi=Epi, Block=Block,
             Param=Param, HyperParam=HyperParam, Input=Input,
             Π=Π, Counter=Counter, Novelty=Novelty,
             Q_N_MF=Q_N_MF, Q_R_MF=Q_R_MF, E_N=E_N, E_R=E_R,
             Q_N_MB=Q_N_MB, Q_R_MB=Q_R_MB, U_N=U_N, U_R=U_R,
             State=State, P_Action=P_Action,
             t_since_change = t_since_change,
             First_Epi=First_Epi)
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Training an agent
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_train_Agent(; Sub=1, Epi=1, Block=1,
                           Param=Str_Param(),
                           HyperParam=Str_HyperParam(),
                           Input=Str_Input(),
                           Agent_Prev=Str_Agent())

    Agent = Func_initial_Agent(;Sub=Sub, Epi=Epi, Block=Block,
                                Param=Param, HyperParam=HyperParam,
                                Input=Input, Agent_Prev=Agent_Prev)
    T = length(Input.obs)
    for t = 2:T
        S = Func_compute_surprise(Agent, t)
        Agent.Surprise = cat(Agent.Surprise, S; dims=1)

        gamma = Func_compute_gamma(Agent,t)
        Agent.Gamma = cat(Agent.Gamma, gamma; dims=1)

        t_since_change = Func_compute_t_since_change(Agent,t)
        Agent.t_since_change = cat(Agent.t_since_change, t_since_change; dims=1)

        counter = Func_update_counter(Agent, t)
        Agent.Counter = cat(Agent.Counter, counter; dims=2)

        Π = Func_update_belief(Agent, t)
        Agent.Π = cat(Agent.Π, Π; dims=3)

        Novelty = Func_update_novelty(Agent, t)
        Agent.Novelty = cat(Agent.Novelty, Novelty; dims=2)

        # Model free updating
        δ_R = Func_compute_δ(Agent, t, "Reward")
        Agent.δ_R_Seq = cat(Agent.δ_R_Seq, δ_R; dims=1)
        δ_N = Func_compute_δ(Agent, t, HyperParam.exploration_form)
        Agent.δ_N_Seq = cat(Agent.δ_N_Seq, δ_N; dims=1)

        E_R = Func_compute_ElTrace(Agent, t, "Reward")
        Agent.E_R = cat(Agent.E_R, E_R; dims=2)
        E_N = Func_compute_ElTrace(Agent, t, HyperParam.exploration_form)
        Agent.E_N = cat(Agent.E_N, E_N; dims=2)

        Q_R_MF = Func_update_Q_MF(Agent, t, "Reward")
        Agent.Q_R_MF = cat(Agent.Q_R_MF, Q_R_MF; dims=2)
        Q_N_MF = Func_update_Q_MF(Agent, t, HyperParam.exploration_form)
        Agent.Q_N_MF = cat(Agent.Q_N_MF, Q_N_MF; dims=2)

        # Model based updating
        if (!Agent.First_Epi)||(t==T)||(Block==2)||(HyperParam.First_Epi_reward_update)
            Q_R_MB, U_R = Func_update_prior_sweep(Agent, t, "Reward")
        else
            Q_R_MB = Agent.Q_R_MB[:,1]
            U_R = Agent.U_R[:,1]
        end
        Agent.Q_R_MB = cat(Agent.Q_R_MB, Q_R_MB; dims=2)
        Agent.U_R = cat(Agent.U_R, U_R; dims=2)

        Q_N_MB, U_N = Func_update_prior_sweep(Agent, t, HyperParam.exploration_form)
        Agent.Q_N_MB = cat(Agent.Q_N_MB, Q_N_MB; dims=2)
        Agent.U_N = cat(Agent.U_N, U_N; dims=2)

        State = Func_update_State(Agent, t)
        Agent.State = cat(Agent.State, State; dims=1)

        if t<T
            if Agent.First_Epi
                P_Action = Func_compute_P_action(Agent, t; Only_Reward=false)
            else
                P_Action = Func_compute_P_action(Agent, t; Only_Reward=true)
            end

            Agent.P_Action = cat(Agent.P_Action, P_Action; dims=2)
        end

        Agent.Novelty_Seq = cat(Agent.Novelty_Seq,
                                Agent.Novelty[Input.obs[t]+1,t-1];
                                dims=1)
    end

    # --------------------------------------------------------------------------
    # Action Analysis - LogL
    for t = 1:(T-1)
        Agent.LogL += log(Agent.P_Action[Input.act[t]+1 ,t])

        (p, a) = findmax(Agent.P_Action[: ,t])
        if t==1
            Agent.Est_Act[1] = a-1
            Agent.Est_Act_P[1] = p
        else
            Agent.Est_Act = cat(Agent.Est_Act, a-1; dims=1)
            Agent.Est_Act_P = cat(Agent.Est_Act_P, p; dims=1)
        end
    end

    Agent.LogL_Unif = (T-1) * log( 1/Param.num_act )
    # --------------------------------------------------------------------------
    # Action Analysis - Acc Rate
    Agent.Sort_Dist_Data = zeros(size(Agent.P_Action))
    Agent.Sort_Dist_Model = zeros(size(Agent.P_Action))

    for t = 1:(T-1)
        temp = sortperm(Agent.P_Action[:,t])
        Agent.Sort_Dist_Model[:,t] = Agent.P_Action[temp,t]

        Agent.Sort_Dist_Data[Input.act[t] + 1 ,t] = 1
        Agent.Sort_Dist_Data[:,t] = Agent.Sort_Dist_Data[temp,t]

        if (Agent.Sort_Dist_Model[findmax(Agent.Sort_Dist_Data[:,t])[2],t] -
            Agent.Sort_Dist_Model[end,t]) == 0

            Agent.Sort_Dist_Data[:,t] = zeros(size(Agent.Sort_Dist_Data[:,t]))

            Agent.Sort_Dist_Data[Agent.Sort_Dist_Model[:,t] .==
                              Agent.Sort_Dist_Model[end,t],t] .= 1

            Agent.Sort_Dist_Data[:,t] = Agent.Sort_Dist_Data[:,t] ./
                                     sum(Agent.Sort_Dist_Data[:,t])
        end
    end

    Agent.LogL_Unif = (T-1) * log( 1/Param.num_act )
    Agent.Acc_Rate = mean(Agent.Sort_Dist_Data[4,:])

    return Agent
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Training a set of agents
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_train_set_of_Agents(BehavData;
                                  Sub_Set = [1], Epi_Set=[1], Block_Set=1:2,
                                  Param=[Str_Param(),Str_Param()],
                                  HyperParam=Str_HyperParam())

    LogL = 0
    N_sub = length(Sub_Set)
    N_blo = length(Block_Set)
    N_epi = length(Epi_Set)
    Agents = Array{Str_Agent,3}(undef,N_sub,N_blo,N_epi)
    for Sub_Ind = 1:length(Sub_Set)
        Sub = Sub_Set[Sub_Ind]
        for Block = Block_Set
            for Epi = Epi_Set
                Input = Func_data_to_input(BehavData; Sub = Sub,
                                           Epi = Epi, Block = Block)
                if (Block==1)&&(Epi==1)
                    Agent_Epi = Func_train_Agent(; Sub=Sub, Epi=Epi,
                                                   Block=Block,
                                                   Param=Param[Block],
                                                   HyperParam=HyperParam,
                                                   Input=Input)
                else
                    if Epi==1
                        Agent_Prev = Agents[Sub_Ind,Block-1,end]
                    else
                        Agent_Prev = Agents[Sub_Ind,Block,Epi-1]
                    end

                    Agent_Epi = Func_train_Agent(; Sub=Sub, Epi=Epi,
                                                   Block=Block,
                                                   Param=Param[Block],
                                                   HyperParam=HyperParam,
                                                   Input=Input,
                                                   Agent_Prev=Agent_Prev)
                end

                Agents[Sub_Ind,Block,Epi] = Agent_Epi
                LogL += Agent_Epi.LogL
            end
        end
    end
    return Agents, LogL
end


text:
# Compute the finite difference approximation of ∂ₓy and ∂ₓₓy at node (i) on an even/uneven grid
function fntdiff1d(x::StaticVector{N,T}, y::AbstractArray{T,1}, i::Int64) where {T<:Union{Float64, BigFloat}, N}
    length(y) == N || throw(DimensionMismatch("x and y in fntdiff1d"))

    if 1 < i < N
        ∂ₓy = /( y[i+1] - y[i-1] , x[i+1] - x[i-1] )

        ∂ₓₓy = /(+ *(x[i]   - x[i-1] , y[i+1])
                 - *(x[i+1] - x[i-1] , y[i]  )
                 + *(x[i+1] - x[i]   , y[i-1])
                ,
                 *(1/2            ,
                   x[i+1] - x[i-1],
                   x[i+1] - x[i]  ,
                   x[i]   - x[i-1])
                )
    elseif 1 == i
        ∂ₓy = /( y[i+1] - y[i] , x[i+1] - x[i] )

        ∂ₓₓy = /(+ *(x[i+1] - x[i]   , y[i+2])
                 - *(x[i+2] - x[i]   , y[i+1])
                 + *(x[i+2] - x[i+1] , y[i]  )
                ,
                 *(1/2            ,
                   x[i+2] - x[i]  ,
                   x[i+2] - x[i+1],
                   x[i+1] - x[i]  )
                )
    elseif N == i
        ∂ₓy = /( y[i] - y[i-1] , x[i] - x[i-1] )

        ∂ₓₓy = /(+ *(x[i-1] - x[i-2] , y[i]  )
                 - *(x[i]   - x[i-2] , y[i-1])
                 + *(x[i]   - x[i-1] , y[i-2])
                ,
                 *(1/2            ,
                   x[i]   - x[i-2],
                   x[i]   - x[i-1],
                   x[i-1] - x[i-2])
                )
    else
        throw(BoundsError(x, i))
    end

    return ∂ₓy, ∂ₓₓy
end

# Compute the finite difference approximation of ∂ₓy and ∂ₓₓy at all nodes on an even/uneven grid
function fntdiff1d!(∂ₓy::AbstractArray{T,1}, ∂ₓₓy::AbstractArray{T,1}, x::StaticVector{N,T}, y::AbstractArray{T,1}) where {T<:Union{Float64, BigFloat}, N}
    length(y) == length(∂ₓy) == length(∂ₓₓy) == N || throw(DimensionMismatch("x, y, ∂ₓy, and ∂ₓₓy in fntdiff1d!"))

    for i in 1:N
        ∂ₓy[i], ∂ₓₓy[i] = fntdiff1d(x, y, i)
    end
    return nothing
end

max_stars_repo_name: cziFinEcon/HPCT


text:
"""
Created on Fri 15 May 2020

@author: Yoann Pradat

Performs NMF rank selection by the method of consensus matrices and cophenetic correlation on the
the first 200 genes of the 5000 x 38 leukemia dataset used in the article of Brunet et al.
"""

using CSV
using Distributed
using DelimitedFiles
using Distributions

@everywhere include("./src/VNMF.jl")
@everywhere using .VNMF

df_pos = CSV.read("../../../data/all_aml_38/golub_nmf_R_200.txt", delim='\t')
mt_pos = convert(Matrix{Float64}, df_pos[:, 2:end])

#### # 2. PARAMETERS ALGORITHM
#### # ##############################################################################################################

#### params for rank selection
rc_params = RCParams(
    K_min = 2,
    K_max = 4,
    n_iter = 10
)

#### nmf β-divergence, multiplicative updates
nmf_global_params = NMFParams(
    init          = :random,
    dist          = truncated(Normal(1, 1), 0, Inf),
    max_iter      = convert(Int, 1e4),
    stopping_crit = :conn,
    stopping_iter = 100,
    verbose       = false,
)

nmf_local_params = FIParams(
    β         = 1,
    l₁ratio_H = 0,
    α_H       = 1,
    alg       = :mm
)

#### NMF struct
nmf = NMF(
    solver        = nmf_FI,
    global_params = nmf_global_params,
    local_params  = nmf_local_params
)

#### get consensus metrics for each rank
@time rc_results = rank_by_consensus(mt_pos, rc_params, nmf)

#### # 3. PLOT WITH PYCALL
#### # ################################################################################################################

#using Pkg
#ENV["Python"] = "/usr/local/anaconda3/bin/python"
#Pkg.build("PyCall")

using PyCall

py"""
import os
import sys
sys.path.append(".")

from demo.src.plot_consensus import *

def make_plot(X, rank):
    folder   = "demo/plot/simulation/rank_cophenetic"
    os.makedirs(folder, exist_ok=True)
    filename = "rank_brunet_consensus_%s.png"  % rank
    filepath = os.path.join(folder, filename)

    config   = ConsensusConfig()

    fig = plot_consensus(X, config)
    fig.savefig(filepath, bbox_inches='tight')
    print("file saved at %s" % filepath)
"""

for (cons, rank) in zip(rc_results.list_cons, rc_results.list_rank)
    py"make_plot"(cons, rank)
end


text:
# First compile all functions in RT_tracking_functions_list.jl
#

## Section 1: import packages

using Images, ImageView, ImageDraw, ImageFeatures, Gtk.ShortNames, VideoIO, ImageFiltering
using SensorFeatureTracking

#######
## Section 2: Setup structs and use Harris corner detection to create a list of features

const TrackingType_Forward = "Forward"
const TrackingType_Pyramid = "Pyramid"
const TrackingType_Inversed = "Inversed"
const TrackingType_InversedPyramid = "InversedPyramid"

const windowSize_25 = 25            # this is half the size of the window to look for a matching feature. Actual window size is windowSize * 2 + 1
const windowSize_20 = 20            # windowSize = windowSize_20 will result in a 41x41 window
const windowSize_15 = 15
const windowSize_10 = 10

frame_counter_from = 2             # in an image sequence you can select a range of images from frame_counter_from + 1 to frame_counter_to
frame_counter_to = 213
nfeatures=10

frame_counter = 2
window_counter = 0
update_reference_frame_count = 10

number_frames = frame_counter_to - frame_counter_from +1

img = load(joinpath(dirname(@__FILE__),"../Data/testSequence/image_$(frame_counter_from - 1).jpg"))
img_end = load(joinpath(dirname(@__FILE__),"../Data/testSequence/image_$(frame_counter_to).jpg"))

I = Gray.(img);
orgI_setup = deepcopy(img)
orgI_setup = Gray.(orgI_setup)
orgI_setup = convert(Array{Float64}, orgI_setup)

#harris corner detection=============

corners = getApproxBestHarrisInWindow(I, nfeatures=nfeatures)   # function in RT_tracking_functions_list.jl
number_features_setup = length(corners[:,1])

#harris corner detection=============

# images used to draw and display the tracked feature's path
feature_path_img_start = deepcopy(img)
feature_path_img_start = Gray.(feature_path_img_start);

feature_path_img_end = deepcopy(img_end)
feature_path_img_end = Gray.(feature_path_img_end);

#constructor function for structs ITVar and ITConst
ITVar, ITConst = ImageTrackerSetup(orgI_setup, corners, windowSize = windowSize_15, TrackingType_setup = TrackingType_Inversed)
# fillNewImageTemplates!(ITVar, ITConst)
#########
# Section 3: Run KLT tracker on the selected image sequence
frame_counter = 2
window_counter = 0
tic()
while frame_counter <= number_frames
    @show frame_counter
    window_counter += 1
    # Load next image
    NextFrame = load(joinpath(dirname(@__FILE__),"../Data/testSequence/image_$(frame_counter + frame_counter_from - 1).jpg"))
    NextFrame = Gray.(NextFrame)
    NextFrame = convert(Array{Float64}, NextFrame)
    I_nextFrame = convert(Array{Float64}, NextFrame)

    ITVar.I_nextFrame = I_nextFrame
    if (ITConst.TrackingType == "Pyramid" || ITConst.TrackingType == "InversedPyramid")
        ITVar.I_nextFrame_downsample = imresize(I_nextFrame, (Int(length(I_nextFrame[:,1])/ITConst.downsampleFactor),Int(length(I_nextFrame[1,:])/ITConst.downsampleFactor)));
    end

        # Main KLT tracker function ==============
        # @time KTL_Tracker!(ITVar, ITConst)
        KTL_Tracker!(ITVar, ITConst)
        # Main KLT tracker function ==============

    # Draw the path of the feature on both the first and last image in the sequence
    # This code draws lines between the reference frame and the tracked feature and not sequential frames.
    # If you need a frame to frame feature path set update_reference_frame_count = 1, however this will decrease the accuracy of the tracker
        for columnCount = 1:length(ITVar.p_reference[1,:])
            if (ITVar.p_reference[2, columnCount][1] <= 0   ||   ITVar.p_reference[2, columnCount][2] <= 0)
                @show fail=1
            else
                draw!(feature_path_img_start, LineSegment(ITVar.p_reference[1, columnCount],ITVar.p_reference[2, columnCount]))
            end
        end
        for columnCount = 1:length(ITVar.p_reference[1,:])
            if (ITVar.p_reference[2, columnCount][1] <= 0   ||   ITVar.p_reference[2, columnCount][2] <= 0)
                @show fail=1
            else
                draw!(feature_path_img_end, LineSegment(ITVar.p_reference[1, columnCount],ITVar.p_reference[2, columnCount]))
            end
        end

    # Update the frame used as reference for tracking
    # Making update_reference_frame_count larger reduces feature drift but can make the tracker lose the feature if the feature change to much
    if (window_counter == update_reference_frame_count)
        # if frame_counter == 96
        #     @show ITVar
        # end
        fillNewImageTemplates!(ITVar, ITConst)
        window_counter = 0
    end

    frame_counter += 1
end
toc()
ImageView.imshow(feature_path_img_start)
ImageView.imshow(feature_path_img_end)
####

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
export TAMessage, TAArc, TA, TTS, MSG_DIRECTION, MSG_IN, MSG_OUT, MSG_SILENT

@enum MSG_DIRECTION MSG_IN=1 MSG_OUT=0 MSG_SILENT=-1

struct TAMessage
    direction::MSG_DIRECTION
    symbol::Symbol
end
TAMessage() = TAMessage(MSG_SILENT, Symbol("—"))
function Base.show(io::IO, m::TAMessage)
    dir = if m.direction==MSG_IN
        "?"
    elseif m.direction==MSG_OUT
        "!"
    else
        ""
    end
    print(io, "$(m.symbol)$(dir)")
end

Base.iterate(x::TAMessage) = (x, nothing)
Base.iterate(::TAMessage, ::Any) = nothing
Base.length(x::TAMessage) = 1
Base.Broadcast.broadcastable(c::TAMessage) = Ref(c)
Base.adjoint(c::TAMessage) = c
Base.copy(m::TAMessage) = TAMessage(m.direction,m.symbol)
Base.deepcopy(m::TAMessage) = TAMessage(deepcopy(m.direction),deepcopy(m.symbol))

mutable struct TAArc{State,T}
    source::State
    target::State
    guard::TAConstraint{T}
    message::TAMessage
    resets::Set{Symbol}
end

Base.iterate(x::TAArc) = (x, nothing)
Base.iterate(::TAArc, ::Any) = nothing
Base.length(x::TAArc) = 1
Base.Broadcast.broadcastable(c::TAArc) = Ref(c)
Base.adjoint(c::TAArc) = c
Base.copy(a::T) where T<:TAArc = T(a.source,a.target,a.guard,a.message,a.resets)
Base.deepcopy(a::T) where T<:TAArc = T(deepcopy(a.source),deepcopy(a.target),deepcopy(a.guard),deepcopy(a.message),deepcopy(a.resets))
function Base.show(io::IO, arc::TAArc)
    m=repr(arc.message)
    r=join(arc.resets, ",")
    r=isempty(r) ? "" : "{$(r)}"
    g=repr(arc.guard)
    i=join(filter(!isempty, (m,g,r)),",")
    i = isempty(i) ? i : "($(i))"
    print(io, "$(repr(arc.source))——$(i)—⟶$(repr(arc.target))")
end

function set_clock_list!(arc::TAArc{State,T}, clocks::Vector{Symbol}; mapping=identity) where {State,T}
    set_clock_list!(arc.guard, clocks; mapping)
    new_resets = mapping.(arc.resets)
    empty!(arc.resets)
    union!(arc.resets, new_resets)
    return arc
end

struct TAState{CombinedState}
    components::CombinedState
    function TAState(arg::T, args...) where T
        s = (arg,args...)
        new{typeof(s)}(s)
    end
end
Base.:(==)(s::TAState, o::TAState) = s.components==o.components
Base.hash(s::TAState) = hash(s.components)
Base.iterate(x::TAState) = iterate(x.components)
Base.iterate(x::TAState, s) = iterate(x.components, s)
Base.length(x::TAState) = length(x.components)
Base.lastindex(x::TAState) = Base.lastindex(x.components)
Base.getindex(x::TAState, s) = x.components[s]
Base.Broadcast.broadcastable(x::TAState) = x.components
Base.adjoint(x::TAState) = x
Base.copy(x::T) where T<:TAState = T(x.components)
Base.deepcopy(x::T) where T<:TAState = TAState(deepcopy.(x.components)...)
Base.:|(args::TAState...) = TAState(flatten(x.components for x in args)...)
function Base.show(io::IO, x::TAState)
    s=join(repr.(x.components),"|")
    print(io, "⟨$(s)⟩")
end

struct TA{State, T}
    states::Vector{State}                    # set of states
    initial_state::State                           # initial state
    clocks::Vector{Symbol}                   # set of clocks - the numeric index is used internally
    symbols::Set{Symbol}                         # alphabet of symbols
    arcs::Vector{TAArc{State,T}}       # arcs of the TA
    arcs_by_state::DefaultDict{State, Vector{TAArc{State,T}}}
    invariants::DefaultDict{State, TAConstraint{T}}   # invariants of the TA's states
    function TA(states::Vector{State},initial_state::State,clocks::Vector{Symbol},symbols::Set{Symbol},arcs::Vector{TAArc{State,T}},invariants::Dict{State, TAConstraint{T}}) where {State<:TAState, T}
        # make sure all arcs have the same clock-set
        arcs_by_state = DefaultDict{State, Vector{TAArc{State,T}}}(()->Vector{TAArc{State,T}}())

        for arc ∈ arcs
            source_id = findfirst(==(arc.source), states)
            @assert !isnothing(source_id) "Source state '$(arc.source)' of arc '$(arc)' not in set of states."
            target_id = findfirst(==(arc.target), states)
            @assert !isnothing(target_id) "Target state '$(arc.target)' of arc '$(arc)' not in set of states."

            # make sure the source/target actually refer to the exact object in the TA's state-vector
            arc.source = states[source_id]
            arc.target = states[target_id]

            set_clock_list!(arc.guard, clocks)
            for reset ∈ arc.resets
                @assert reset ∈ clocks "Clock '$(reset)' should be reset by arc '$(arc)', but no such clock is known! Clock list: '$(clocks)'."
            end
            push!(arcs_by_state[arc.source], arc)
        end

        # make sure all guards have the same clock-set, make sure there is an (empty?) invariance for each state
        default_guard = @constraint( true, T )
        set_clock_list!(default_guard, clocks)

        new_invariants = DefaultDict(default_guard, invariants)
        for (state,I) ∈ pairs(invariants)
            set_clock_list!(I, clocks)
            state_id = findfirst(==(state), states)
            @assert !isnothing(state_id) "State '$(state)' of guard expression '$(I)' not in the set states."
            new_invariants[states[state_id]] = I
        end

        new{State, T}(states,initial_state,clocks,symbols,arcs,arcs_by_state,new_invariants)
    end
end
TA(states::Vector,initial_state::State,clocks::Vector{Symbol},symbols::Set{Symbol},arcs::Vector{TAArc{State,T}},invariants::Dict{State, TAConstraint{T}}) where {State,T} =
    TA(TAState.(states),TAState(initial_state),clocks, symbols, [TAArc(TAState(a.source),TAState(a.target),a.guard,a.message,a.resets) for a in arcs], Dict(TAState(k)=>v for (k,v) in pairs(invariants)))

Base.copy(ta::T) where T<:TA = TA(ta.states,ta.initial_state,ta.clocks,ta.symbols,ta.arcs,Dict(ta.invariants))
Base.deepcopy(ta::T) where T<:TA = TA(deepcopy(ta.states),deepcopy(ta.initial_state),deepcopy(ta.clocks),deepcopy(ta.symbols),deepcopy(ta.arcs),Dict(deepcopy(ta.invariants)))
Base.iterate(x::TA) = (x, nothing)
Base.iterate(::TA, ::Any) = nothing
Base.length(x::TA) = 1
Base.Broadcast.broadcastable(c::TA) = Ref(c)
Base.adjoint(c::TA) = c

function set_clock_list!(ta::TA{State,T}, clocks::Vector{Symbol}; mapping=identity) where {State,T}
    clocks=copy(clocks)
    # update clock for each state invariant
    for state in ta.states
        set_clock_list!(ta.invariants[state], clocks; mapping)
    end

    # update clock for each arc guard & reset
    #should also apply changes to arcs_by_state, because they share the same objects
    for arc in ta.arcs
        set_clock_list!(arc, clocks; mapping)
    end

    empty!(ta.clocks)
    append!(ta.clocks,clocks)

    return ta
end

# Timed Transition System
struct TTS{State, T}
    states::Vector{Tuple{State,TAConstraint{T}}}
    transitions::Vector{Pair{Tuple{State,TAConstraint{T}}, Tuple{TAMessage, Tuple{State,TAConstraint{T}}}}}
    function TTS(ta::TA{State,T}, args...; kwargs...) where {State<:TAState, T}
        new{State,T}(TimedAutomata.zone_graph_fwd(ta, args...; kwargs...)...)
    end
end

Base.iterate(x::TTS) = (x, nothing)
Base.iterate(::TTS, ::Any) = nothing
Base.length(x::TTS) = 1
Base.Broadcast.broadcastable(c::TTS) = Ref(c)
Base.adjoint(c::TTS) = c
Base.copy(tts::T) where T<:TTS = T(tts.states,tts.transitions)
Base.deepcopy(tts::T) where T<:TTS = T(deepcopy(tts.states),deepcopy(tts.transitions))


text:
function get_module_3_excersise_3_author_name()
    return "YOUR NAME HERE"
end

function problem_3(A)
    dimension_of_A = size(A)[1]
    is_symmetric = true
    #implement your loops here
    return is_symmetric #return the x vector
end


text:
module ESS

# These methods have been deprecated / moved
macro current_module()
    return VERSION >= v"0.7-" ? :(@__MODULE__) : :(current_module())
end

parse = VERSION >= v"0.7-" ? Base.Meta.parse : Base.parse
function_module = VERSION >= v"0.7-" ? Base.parentmodule : Base.function_module

function all_help_topics()
    ## There are not clear topics anymore. Approximate those with a very general
    ## apropos(" ")
    Base.Docs.apropos(" ")
end

function help(topic::AbstractString)
    if (VERSION >= v"1.0-")
        Core.eval(parentmodule(ESS), parse("@doc $topic"))
    elseif  (VERSION >= v"0.4-")
        Core.eval(@current_module(), parse("@doc $topic"))
    else
        Base.Help.help(topic)
    end
end

## modified version of function show(io::IO, m::Method)
function fun_args(m::Method)
    tv, decls, file, line = Base.arg_decl_parts(m)
    io = VERSION >= v"0.7-" ? Base.stdout : STDOUT::IO # STDOUT is no longer in 1.0
    if !isempty(tv)
        Base.show_delim_array(io, tv, '{', ',', '}', false)
    end
    print(io, "(")
    join(io, [escape_string(isempty(d[2]) ? d[1] : d[1]*"::"*d[2]) for d in decls],
         ",", ",")
    Base.print(io, ")")
end

## modified versionof show(io::IO, mt::MethodTable)
function fun_args(f::Function)
    mt = Base.MethodList(methods(f).mt)
    mod = function_module(f)   # Base.function_module deprecated in 0.7
    if mod == Main
        mod = "nil"
    end
    print("(list \"$mod\" nil '(")
    for d in mt
        print("\"")
        ## method
        fun_args(d)
        print("\" ")
    end
    print("))")
end

function fun_args(s::AbstractString)
    try
        mod = VERSION >= v"1.0-" ?  parentmodule(ESS) : @current_module()
        m = Core.eval(mod, parse(s))
        if ! isa(m, String)
            fun_args(m)
        end
    catch
        print("(list nil nil nil)")
    end
end

function fun_args(t::DataType)
    print("(list nil nil '(")
    for d = fieldnames(t)
        print("\"$d\" ")
    end
    print("))")
end

### OBJECT COMPLETION
# Must print an output of the form:
#
# Cache                         Module
# Write                         Module
# add                           Function
# free                          Function
function components(m::Module)
    for v in sort(names(m, all=true, imported=true))
        s = string(v)
        if !startswith(s, "#") && isdefined(m,v)
            println(rpad(s, 30), summary(Core.eval(m,v)))
        end
    end
end

function components(t::DataType)
    for v in sort(fieldnames(t))
        println(rpad(string(v), 30), "field")
    end
end

function components(v)
    t = typeof(v)
    if isa(t, DataType)
        return components(t)
    end
end

### MISC
function main_modules(m::Module)
    for nm in names(m)
        if isdefined(m, nm)
            mod = Core.eval(m, nm)
            if isa(mod, Module)
                print("\"$nm\" ")
            end
        end
    end
end

if VERSION >= v"0.7-"
    main_modules() = main_modules(Base.parentmodule(@current_module()))
else
    main_modules() = main_modules(@current_module())
end

end


text:
"""
    @reset_parallel_stencil()

Reset the ParallelStencil module.

See also: [`init_parallel_stencil`](@ref)
"""
macro reset_parallel_stencil() esc(reset_parallel_stencil(__module__)) end

function reset_parallel_stencil(caller::Module)
    ParallelKernel.reset_parallel_kernel(caller)
    set_initialized(false)
    set_package(PKG_NONE)
    set_numbertype(NUMBERTYPE_NONE)
    set_ndims(NDIMS_NONE)
    return nothing
end


text:
"""
    Constraints Data
"""
struct ConstraintsData{T,C,CX,CU}
    constraints::Constraints{T}
    violations::Vector{C}
    jacobian_state::Vector{CX}
    jacobian_action::Vector{CU}
end

function constraint_data(model::Model, constraints::Constraints)
    H = length(constraints)
    c = [zeros(constraints[t].num_constraint) for t = 1:H]
    cx = [zeros(constraints[t].num_constraint, t < H ? model[t].num_state : model[H-1].num_next_state) for t = 1:H]
    cu = [zeros(constraints[t].num_constraint, model[t].num_action) for t = 1:H-1]
    ConstraintsData(constraints, c, cx, cu)
end

function constraint!(constraint_data::ConstraintsData, x, u, w)
    constraint!(constraint_data.violations, constraint_data.constraints, x, u, w)
end

function constraint_violation(constraint_data::ConstraintsData;
    norm_type=Inf)

    constraints = constraint_data.constraints
    H = length(constraints)
    max_violation = 0.0
    for t = 1:H
        num_constraint = constraints[t].num_constraint
        ineq = constraints[t].indices_inequality
        for i = 1:num_constraint
            c = constraint_data.violations[t][i]
            cti = (i in ineq) ? max(0.0, c) : abs(c)
            max_violation = max(max_violation, cti)
        end
    end
    return max_violation
end

function constraint_violation(constraint_data::ConstraintsData, x, u, w;
    norm_type=Inf)
    constraint!(constraint_data, x, u, w)
    constraint_violation(constraint_data,
        norm_type=norm_type)
end


text:
function solve_step!(
  x1::NTuple{3,AbstractVector},
  solver::Newmark,
  op::AffineODEOperator,
  x0::NTuple{3,AbstractVector},
  t0::Real,
  cache) # -> (uF,tF)

  dt = solver.dt
  γ = solver.γ
  β = solver.β
  t1 = t0+dt
  u0, v0, a0 = x0
  u1, v1, a1 = x1
  newmatrix = true

  if cache === nothing
    # Allocate caches
    newmark_cache = allocate_cache(op,v0,a0)
    (v,a, ode_cache) = newmark_cache

    # Allocate matrices and vectors
    A, b = _allocate_matrix_and_vector(op,x0,ode_cache)

    # Create affine operator cache
    affOp_cache = (A,b,nothing)
  else
    newmark_cache, affOp_cache = cache
  end

  # Unpack and update caches
  (v,a, ode_cache) = newmark_cache
  ode_cache = update_cache!(ode_cache,op,t1)
  A,b,l_cache = affOp_cache

  # Define Newmark operator
  newmark_affOp = NewmarkAffineOperator(op,t1,dt,γ,β,(u0,v0,a0),newmark_cache)

  # Fill matrix and vector
  _matrix_and_vector!(A,b,newmark_affOp,u1)

  # Create affine operator with updated RHS
  affOp = AffineOperator(A,b)
  l_cache = solve!(u1,solver.nls,affOp,l_cache,newmatrix)

  # Update auxiliar variables
  u1 = u1 + u0
  v1 = γ/(β*dt)*(u1-u0) + (1-γ/β)*v0 + dt*(1-γ/(2*β))*a0
  a1 = 1.0/(β*dt^2)*(u1-u0) - 1.0/(β*dt)*v0 - (1-2*β)/(2*β)*a0

  # Pack caches
  affOp_cache = A,b,l_cache
  cache = (newmark_cache, affOp_cache)
  x1 = (u1,v1,a1)

  return (x1,t1,cache)

end

"""
Affine operator that represents the Newmark Affine operator at a
given time step, i.e., M(t)(u_n+1-u_n)/dt + K(t)u_n+1 + b(t)
"""
struct NewmarkAffineOperator <: NonlinearOperator
  odeop::AffineODEOperator
  t1::Float64
  dt::Float64
  γ::Float64
  β::Float64
  x0::NTuple{3,AbstractVector}
  ode_cache
end

function _matrix_and_vector!(
  A::AbstractMatrix,
  b::AbstractVector,
  affOp::NewmarkAffineOperator,
  x::AbstractVector)
  jacobian!(A,affOp,x)
  residual!(b,affOp,x)
end

function residual!(b::AbstractVector,op::NewmarkAffineOperator,x::AbstractVector)
  u1 = x
  u0, v0, a0 = op.x0
  v1, a1, cache = op.ode_cache
  a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
  v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
  residual!(b,op.odeop,op.t1,(u1,v1,a1),cache)
  b .*= -1.0
end

function jacobian!(A::AbstractMatrix,op::NewmarkAffineOperator,x::AbstractVector)
  u1 = x
  u0, v0, a0 = op.x0
  v1, a1, cache = op.ode_cache
  a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
  v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
  z = zero(eltype(A))
  fillstored!(A,z)
  jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
end

function _allocate_matrix(odeop::ODEOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
  A = allocate_jacobian(odeop,x[1],ode_cache)
  return A
end

function _allocate_matrix_and_vector(odeop::ODEOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
  b = allocate_residual(odeop,x[1],ode_cache)
  A = allocate_jacobian(odeop,x[1],ode_cache)
  return A, b
end

# # function allocate_residual(op::NewmarkAffineOperator,x::AbstractVector)
# #   v1, a1, cache = op.ode_cache
# #   allocate_residual(op.odeop,x,cache)
# # end

# # function allocate_jacobian(op::NewmarkAffineOperator,x::AbstractVector)
# #   v1, a1, cache = op.ode_cache
# #   allocate_jacobian(op.odeop,x,cache)
# # end

# function _allocate_matrix_and_vector(op::NewmarkAffineOperator,x::AbstractVector)
#   A = allocate_jacobian(op,x)
#   b = allocate_residual(op,x)
#   return A,b
# end


text:
using SafeTestsets

@safetestset "Sims 2010" begin include("Sims2010.jl") end
@safetestset "Mackowiack and Wiederholt 2009" begin include("MW2009.jl") end
@safetestset "Endogenous Feedback" begin include("EndFeedback.jl") end


text:
@doc raw"""
    HeteroscedasticLikelihood(λ::T=1.0)->HeteroscedasticGaussianLikelihood

## Arguments
- `λ::Real` : The maximum precision possible (this is optimized during training)

---

Gaussian with heteroscedastic noise given by another gp:
```math
    p(y|f,g) = \mathcal{N}(y|f,(\lambda \sigma(g))^{-1})
```
Where ``\sigma`` is the logistic function

The augmentation is not trivial and will be described in a future paper
"""
HeteroscedasticLikelihood(λ::Real) = HeteroscedasticGaussianLikelihood(InvScaledLogistic(λ))

struct InvScaledLogistic{T} <: AbstractLink
    λ::Vector{T}
end

InvScaledLogistic(λ::Real) = InvScaledLogistic([λ])

(l::InvScaledLogistic)(f::Real) = inv(l.λ[1] * logistic(f))

implemented(::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic}, ::AnalyticVI) = true

function (l::HeteroscedasticGaussianLikelihood)(f, y::Real)
    return pdf(l(f), y)
end

function Distributions.loglikelihood(
    l::HeteroscedasticGaussianLikelihood, f::AbstractVector, y
)
    return logpdf(l(f), y)
end

function Base.show(io::IO, ::HeteroscedasticGaussianLikelihood)
    return print(io, "Gaussian likelihood with heteroscedastic noise")
end

n_latent(::HeteroscedasticGaussianLikelihood) = 2

function init_local_vars(
    ::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},
    batchsize::Int,
    T::DataType=Float64,
)
    return (;
        c=ones(T, batchsize),
        ϕ=ones(T, batchsize),
        γ=ones(T, batchsize),
        θ=ones(T, batchsize),
        σg=ones(T, batchsize),
    )
end

function compute_proba(
    l::HeteroscedasticGaussianLikelihood,
    μ::AbstractVector{<:AbstractVector},
    Σ::AbstractVector{<:AbstractVector},
)
    return μ[1], Σ[1] .+ l.invlink.(μ[2])
end

function local_updates!(
    local_vars,
    l::HeteroscedasticGaussianLikelihood,
    y::AbstractVector,
    μ::NTuple{2,<:AbstractVector},
    diagΣ::NTuple{2,<:AbstractVector},
)
    # gp[1] is f and gp[2] is g (for approximating the noise)
    @. local_vars.ϕ = 0.5 * (abs2(μ[1] - y) + diagΣ[1])
    @. local_vars.c = sqrt(abs2(μ[2]) + diagΣ[2])
    @. local_vars.γ =
        0.5 * l.invlink.λ[1] * local_vars.ϕ * safe_expcosh(-0.5 * μ[2], 0.5 * local_vars.c)
    @. local_vars.θ = 0.5 * (0.5 + local_vars.γ) / local_vars.c * tanh(0.5 * local_vars.c)
    @. local_vars.σg = expectation(logistic, μ[2], diagΣ[2])
    l.invlink.λ .= max(
        0.5 * length(local_vars.ϕ) / dot(local_vars.ϕ, local_vars.σg), l.invlink.λ[1]
    )
    return local_vars
end

function variational_updates!(
    m::AbstractGPModel{
        T,<:HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},<:AnalyticVI
    },
) where {T}
    local_vars = local_updates!(
        state.local_vars, likelihood(m), yview(m), mean_f(m), var_f(m)
    )
    natural_gradient!(
        m.f[2],
        ∇E_μ(likelihood(m), opt(inference(m)), yview(m), local_vars)[2],
        ∇E_Σ(likelihood(m), opt(inference(m)), yview(m), local_vars)[2],
        ρ(m),
        opt(inference(m)),
        last(Zviews(m)),
        state.kernel_matrices[2],
        state.opt_state[2],
    )
    vi_opt_state_2 = global_update!(
        m.f[2], opt(inference(m)), inference(m), state.opt_state[2]
    )
    local_vars = heteroscedastic_expectations!(
        local_vars, likelihood(m), mean_f(m.f[2]), var_f(m.f[2])
    )
    natural_gradient!(
        m.f[1],
        ∇E_μ(likelihood(m), opt(inference(m)), yview(m), local_vars)[1],
        ∇E_Σ(likelihood(m), opt(inference(m)), yview(m), local_vars)[1],
        ρ(m),
        opt(inference(m)),
        first(Zviews(m)),
        state.kernel_matrices[1],
        state.opt_state[1],
    )
    vi_opt_state_1 = global_update!(
        m.f[1], opt(inference(m)), inference(m), state.opt_state[1]
    )
    opt_state = (vi_opt_state_1, vi_opt_state_2)
    return merge(state, (; opt_state))
end

function heteroscedastic_expectations!(
    local_vars,
    l::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},
    μ::AbstractVector,
    Σ::AbstractVector,
)
    @. local_vars.σg = expectation(logistic, μ, Σ)
    l.invlink.λ .= max(
        0.5 * length(local_vars.ϕ) / dot(local_vars.ϕ, local_vars.σg), l.invlink.λ[1]
    )
    return local_vars
end

@inline function ∇E_μ(
    l::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},
    ::AOptimizer,
    y::AbstractVector,
    state,
)
    return (0.5 * y .* l.invlink.λ[1] .* state.σg, 0.5 * (0.5 .- state.γ))
end

@inline function ∇E_Σ(
    l::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},
    ::AOptimizer,
    ::AbstractVector,
    state,
)
    return (0.5 * l.invlink.λ[1] .* state.σg, 0.5 * state.θ)
end

function compute_proba(
    l::HeteroscedasticGaussianLikelihood,
    μs::Tuple{<:AbstractVector,<:AbstractVector},
    σs::Tuple{<:AbstractVector,<:AbstractVector},
) where {T<:Real}
    return μs[1], σs[1] + expectation.(Ref(l.invlink), μs[2], σs[2])
end

function predict_y(
    ::HeteroscedasticGaussianLikelihood, μs::Tuple{<:AbstractVector,<:AbstractVector}
)
    return first(μs) # For predict_y the variance is ignored
end

function expec_loglikelihood(
    l::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},
    ::AnalyticVI,
    y::AbstractVector,
    μ,
    diag_cov,
    state,
)
    tot = length(y) * (0.5 * log(l.invlink.λ[1]) - log(2 * sqrt(twoπ)))
    tot +=
        0.5 * (
            dot(μ[2], (0.5 .- state.γ)) - dot(abs2.(μ[2]), state.θ) -
            dot(diag_cov[2], state.θ)
        )
    tot -= PoissonKL(l, y, μ[1], diag_cov[1], state)
    return tot
end

function AugmentedKL(
    l::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic}, state, ::Any
)
    return PolyaGammaKL(l, state)
end

function PoissonKL(
    l::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic},
    y::AbstractVector,
    μ::AbstractVector,
    Σ::AbstractVector,
    state,
)
    return PoissonKL(
        state.γ,
        0.5 * l.invlink.λ[1] * (abs2.(y - μ) + Σ),
        log.(0.5 * l.invlink.λ[1] * (abs2.(μ - y) + Σ)),
    )
end

function PolyaGammaKL(::HeteroscedasticGaussianLikelihood{<:InvScaledLogistic}, state)
    return PolyaGammaKL(0.5 .+ state.γ, state.c, state.θ)
end


text:
#SignalProcessing show functions
#-------------------------------------------------------------------------------

#==Generate friendly show functions
===============================================================================#

function Base.show(io::IO, p::Pole{T, AbastractFrequency{U}}) where {T,U}
	print(io, "Pole{$T,$U}($(p.v))")
end

function Base.show(io::IO, ::Type{DS{T}}) where T
	print(io, "DS{:$T}")
end

function Base.show(io::IO, ::DS{T}) where T
	print(io, "DispatchableSymbol:$T")
end

function Base.show(io::IO, d::DataTime)
	strtype = d.data.tperiodic ? "Periodic" : "Finite"
	print(io, "DataTime($strtype, t=", d.data.t, ", x(t)=", d.data.xt, ")")
end

function Base.show(io::IO, d::DataFreq)
	strtype = d.data.tperiodic ? "Series" : "Spectrum"
	print(io, "DataFreq($strtype, f=", d.data.f, ", X(f)=", d.data.Xf, ")")
end

#Last line


text:
# A list of members of the population, with easy constructors,
#  which allow for random generation of new populations
mutable struct Population{T<:Real}
    members::Array{PopMember{T}, 1}
    n::Int
end

Population(pop::Array{PopMember{T}, 1}) where {T<:Real} = Population{T}(pop, size(pop)[1])
Population(pop::Array{PopMember{T}, 1}, npop::Int) where {T<:Real} = Population{T}(pop, npop)
Population(dataset::Dataset{T}, baseline::T;
           npop::Int, nlength::Int=3,
           options::Options,
           nfeatures::Int) where {T<:Real} = Population([PopMember(dataset, baseline, genRandomTree(nlength, options, nfeatures), options) for i=1:npop], npop)
Population(X::AbstractMatrix{T}, y::AbstractVector{T}, baseline::T;
           npop::Int, nlength::Int=3,
           options::Options,
           nfeatures::Int) where {T<:Real} = Population(Dataset(X, y), baseline, npop=npop, options=options, nfeatures=nfeatures)

# Sample 10 random members of the population, and make a new one
function samplePop(pop::Population, options::Options)::Population
    idx = rand(1:pop.n, options.ns)
    return Population(pop.members[idx])
end

# Sample the population, and get the best member from that sample
function bestOfSample(pop::Population, options::Options)::PopMember
    sample = samplePop(pop, options)
    best_idx = argmin([sample.members[member].score for member=1:sample.n])
    return sample.members[best_idx]
end

function finalizeScores(dataset::Dataset{T},
                        baseline::T, pop::Population,
                        options::Options)::Population where {T<:Real}
    need_recalculate = options.batching
    if need_recalculate
        @inbounds @simd for member=1:pop.n
            pop.members[member].score = scoreFunc(dataset, baseline,
                                                  pop.members[member].tree,
                                                  options)
        end
    end
    return pop
end

# Return best 10 examples
function bestSubPop(pop::Population; topn::Int=10)::Population
    best_idx = sortperm([pop.members[member].score for member=1:pop.n])
    return Population(pop.members[best_idx[1:topn]])
end

# Return best 10 examples size that are pareto dominating
function bestSubPopParetoDominating(pop::Population{T}; topn::Int=10)::Population where {T<:Real}
    scores = [pop.members[member].score for member=1:pop.n]
    best_idx = sortperm(scores)

    sorted = pop.members[best_idx]
    score_sorted = scores[best_idx]
    sizes = [countNodes(sorted[i].tree) for i=1:pop.n]

    best = [1]
    for i=2:pop.n
        better_than_all_smaller = true
        for j=i-1:1
            if sizes[j] < sizes[i] #Another model is smaller AND better
                better_than_all_smaller = false
                break
            end
        end
        if better_than_all_smaller
            best = vcat(best, [i])
        end
        if length(best) == topn
            break
        end
    end
    return Population(sorted[best])
end


text:
@testset "Robots: biped" begin
    q0 = [0.0; 0.0; 0.5 * π * ones(7)]
    v0 = zeros(9)
    @test norm(lagrangian(biped, q0, v0)) < 1.0e-8

    # visualize
    vis = RoboDojo.Visualizer();
    @test visualize!(vis, biped, [q0], Δt=0.1);
end


text:
using Documenter, RobRecSolver, RobRecSolver.Experiments

makedocs(
    modules = [
        RobRecSolver,
        RobRecSolver.Experiments
    ]
)
run(`mkdocs build`)


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    envelope(f::PiecewiseQuadratic)

Computes the greatest convex lower bound of a `PiecewiseQuadratic` `f`.

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 0., 0., 1., 0.),
                               BoundedQuadratic(0., 3., 0., 0., 0.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [-Inf, 0.00000]
BoundedQuadratic: f(x) = 0, ∀x ∈ [0.00000, 3.00000]

julia> envelope(f)
Piecewise quadratic function:
BoundedQuadratic: f(x) = + 1.00000 x - 3.00000, ∀x ∈ [-Inf, 3.00000]

```
"""
function envelope(f::PiecewiseQuadratic)
    h = FixedMemoryPwq(2 * length(f))
    env = FixedMemoryPwq(3)  # Pre-allocated workspace.
    for fi in f
        append_envelope!(h, fi, env)
    end
    return PiecewiseQuadratic(h)
end

function append_envelope!(f::FixedMemoryPwq, g::BoundedQuadratic, env::FixedMemoryPwq)
    if length(f) == 0
        push!(f, g)
    else
        intersection_at_left = true
        empty!(env)

        while (length(f) > 0) & intersection_at_left
            empty!(env)
            fi = pop!(f)
            intersection_at_left = append_envelope!(fi, g, env)
        end
        append!(f, env)
    end
end

function append_envelope!(f::BoundedQuadratic, g::BoundedQuadratic, env::FixedMemoryPwq)
    @assert f.ub ≤ g.lb
    @assert length(env) == 0

    _, _ = _envelope_midpt_midpt!(f, g, env)
    if length(env) > 0
        return false
    end

    _, _ = _envelope_midpt_endpt!(f, g, env)
    if length(env) > 0
        return false
    end

    _, iar = _envelope_midpt_endpt!(reverse(g), reverse(f), env)
    if length(env) > 0
        reverse!(env)
        return iar
    end

    ial, _ = _envelope_endpt_endpt!(f, g, env, true)
    if length(env) > 0
        return ial
    end

    _, iar = _envelope_endpt_endpt!(reverse(g), reverse(f), env, false)
    if length(env) > 0
        reverse!(env)
        return iar
    end

    error("Error in convex envelope computation.")

    return false, 0
end

function _envelope_midpt_midpt!(f::BoundedQuadratic, g::BoundedQuadratic,
                                env::FixedMemoryPwq)
    if (f.p > 0) & (g.p > 0)
        a = (f.p^2 / g.p - f.p)
        b = (f.p / g.p) * (f.q - g.q)
        c = f.r - g.r + (f.q - g.q)^2 / (4 * g.p)
        x1, x2 = solve_quad(a, b, c)
        for xf in [x1, x2]
            xg = f.p / g.p * xf + (f.q - g.q) / (2 * g.p)
            if (xf ≳ f.lb) && (xf ≲ f.ub) && (xg ≳ g.lb) && (xg ≲ g.ub)
                xf = clip(xf, f.lb, f.ub)
                xg = clip(xg, g.lb, g.ub)
                m = 2 * f.p * xf + f.q
                env[1] = BoundedQuadratic(f.lb, xf, f.p, f.q, f.r)
                env[2] = BoundedQuadratic(xf, xg, 0.0, m, f(xf) - m * xf)
                env[3] = BoundedQuadratic(xg, g.ub, g.p, g.q, g.r)
                env.len = 3
                return false, false
            end
        end
    end
    return false, false, 0
end

function _envelope_midpt_endpt!(f::BoundedQuadratic, g::BoundedQuadratic,
                                env::FixedMemoryPwq)
    if f.p > 0

        # Midpoint to lower
        xg = g.lb
        xf = _envelope_point_quad(xg, g(xg), f.p, f.q, f.r)
        if !isnan(xf)
            if (xf ≳ f.lb) && (xf ≲ f.ub)
                xf = clip(xf, f.lb, f.ub)
                h = get_tangent(f, xf)
                if h ≲ g
                    clip(xf, f.lb, f.ub)
                    env[1] = restrict_dom(f, -Inf, xf)
                    env[2] = restrict_dom(h, xf, xg)
                    env[3] = restrict_dom(g, xg, Inf)
                    env.len = 3
                    return false, false
                end
            end
        end

        # Midpoint to (finite) upper
        xg = g.ub
        xf = _envelope_point_quad(xg, g(xg), f.p, f.q, f.r)
        if !isnan(xf)
            if (xf ≳ f.lb) && (xf ≲ f.ub)
                xf = clip(xf, f.lb, f.ub)
                h = get_tangent(f, xf)
                if h ≲ g
                    env[1] = restrict_dom(f, -Inf, xf)
                    env[2] = restrict_dom(h, xf, xg)
                    env.len = 2
                    return false, true
                end
            end
        end

        # Midpoint to (infinite) upper
        if isinf(g.ub) & (g.p == 0.0)
            xf = (g.q - f.q) / (2 * f.p)
            if (xf ≳ f.lb) && (xf ≲ f.ub)
                xf = clip(xf, f.lb, f.ub)
                h = BoundedQuadratic(0.0, g.q, extend_dom(f)(xf) - g.q * xf)
                if h ≲ g
                    env[1] = restrict_dom(f, -Inf, xf)
                    env[2] = restrict_dom(h, xf, Inf)
                    env.len = 2
                    return false, false
                end
            end
        end
    end
    return false, false, 0
end

function _envelope_endpt_endpt!(f::BoundedQuadratic, g::BoundedQuadratic,
                                env::FixedMemoryPwq, do_symmetric::Bool)
    if do_symmetric
        # Upper to lower (no gap)
        if f.ub == g.lb
            if is_point(f) && (f(f.ub) ≳ g(g.lb))
                env[1] = g
                env.len = 1
                return true, is_point(g)
            elseif is_point(g) && (g(g.lb) ≳ f(f.ub))
                env[1] = f
                env.len = 1
                return true, is_point(g)
            elseif (f(f.ub) ≈ g(g.lb)) && (2 * f.p * f.ub + f.q ≲ 2 * g.p * g.lb + g.q)
                env[1] = f
                env[2] = g
                env.len = 2
                return is_point(f), is_point(g)
            end
        end

        # Upper to lower (with gap)
        if f.ub != g.lb
            h = get_line(f.ub, f(f.ub), g.lb, g(g.lb))
            if (h ≲ f) && (h ≲ g)
                env[1] = f
                env[2] = restrict_dom(h, f.ub, g.lb)
                env[3] = g
                env.len = 3
                return is_point(f), is_point(g)
            end
        end

        # (Finite) lower to (finite) upper
        if (f.lb > -Inf) && (g.ub < Inf) && (f.lb != g.ub)
            h = get_line(f.lb, f(f.lb), g.ub, g(g.ub))
            if (h ≲ f) && (h ≲ g)
                env[1] = restrict_dom(h, f.lb, g.ub)
                env.len = 1
                return true, true
            end
        end
    end

    # Upper to (finite) upper
    if (g.ub < Inf) && (f.ub != g.ub)
        h = get_line(f.ub, f(f.ub), g.ub, g(g.ub))
        if (h ≲ f) && (h ≲ g)
            env[1] = f
            env[2] = restrict_dom(h, f.ub, g.ub)
            env.len = 2
            return is_point(f), true
        end
    end

    # (Finite) lower to (infinite) upper
    if (f.lb > -Inf) && (g.ub == Inf) && (g.p == 0)
        h = BoundedQuadratic(0.0, g.q, f(f.lb) - g.q * f.lb)
        if (h ≲ f) && (h ≲ g)
            env[1] = restrict_dom(h, f.lb, g.ub)
            env.len = 1
            return true, false
        end
    end

    # Upper to (infinite) upper
    if g.ub == Inf && (g.p == 0)
        h = BoundedQuadratic(0.0, g.q, f(f.ub) - g.q * f.ub)
        if (h ≲ f) && (h ≲ g)
            env[1] = f
            env[2] = restrict_dom(h, f.ub, Inf)
            env.len = 2
            return is_point(f), false
        end
    end

    return false, false
end

function _envelope_point_quad(x1::Float64, y1::Float64, p::Float64, q::Float64, r::Float64)
    a = p
    b = -2 * p * x1
    c = -q * x1 - r + y1

    x1, x2 = solve_quad(a, b, c)
    x_tan = minimum(solve_quad(a, b, c))
    return x_tan
end


text:
module CalibrationErrors

using Reexport

using DataStructures
@reexport using Distances
@reexport using KernelFunctions
using StatsBase: StatsBase
using UnPack: @unpack

using LinearAlgebra
using Statistics

# estimation
export calibrationerror

# estimators
export ECE, BiasedSKCE, UnbiasedSKCE, BlockUnbiasedSKCE, UCME

# binning algorithms
export UniformBinning, MedianVarianceBinning

include("generic.jl")
include("binning/generic.jl")
include("binning/uniform.jl")
include("binning/medianvariance.jl")
include("ece.jl")

include("skce/generic.jl")
include("skce/biased.jl")
include("skce/unbiased.jl")

include("ucme.jl")

include("deprecated.jl")

end # module


text:
# Test specific data for one network:
println("- number/case9 check")
mpc = loadcase("case9")
gencost = [
        2.0  1500.0  0.0  3.0  0.11    5.0  150.0
        2.0  2000.0  0.0  3.0  0.085   1.2  600.0
        2.0  3000.0  0.0  3.0  0.1225  1.0  335.0
        ]
@test mpc["gencost"] == gencost

# Ensure these networks exist in the `data`:
println("- file check")
@test issubset(["case96","case1354pegase"], casenames())

# Ensure every network's mpc Dict contains these keys:
println("- key check")
# key_list = ["gen";"branch" ;"docstring";"bus"]
key_list = ["gen";"branch";"bus"]
key_test = Bool[]
for name in casenames()
    mpc = loadcase(name,describe=false)
    push!(key_test,key_list ⊆ [key for key in keys(mpc)])
end
@test all(key_test)

println("- graceful failure test")
@test_throws ErrorException loadcase("invalidcasename")
@test_throws ErrorException loadcase("doesnotexist.mat")


text:
#SIMULATION RESULTS AND CHARTS FOR MC HAMMER
#by Eric Torkia, April 2019

"""
MCH_Charts contains standard simulation charts for sensitivity, density, trends (time series with confidence bands) for simulation arrays generated by mc_hammer.
"""

# Density Chart
"""
    density_chrt(Data, x_label="Sim. Values")

Data is your array, either simulated or historical.
x_label [optional] allows you to customize your X axis label.
"""
function density_chrt(Data, x_label="Sim. Values")
      if x_label ==""
            x_label="x"
      end
chart = plot(Data, x=Data, Guide.xlabel(x_label), Guide.ylabel("Frequency"), Geom.density)
describe(Data)
println("")
print("Mean: ", mean(Data),"\n")
print("Std.Dev: ", std(Data),"\n")
print("Prob. of Neg.: ", GetCertainty(Data,0,0),"\n")
println("")
print("p10, p50, p90 : ", quantile(collect(Float64, Data),[0.1,0.5,0.9]),"\n")
return chart
end

# Histogram Chart
"""
    histogram_chrt(Data, x_label="Sim. Values")

Data is your array, either simulated or historical.
x_label [optional] allows you to customize your X axis label.
"""
function histogram_chrt(Data, x_label="Sim. Values")
      if x_label ==""
            x_label="x"
      end
chart = plot(Data, x=Data, Guide.xlabel(x_label), Guide.ylabel("Frequency"), Geom.histogram)
describe(Data)
println("")
print("Mean: ", mean(Data),"\n")
print("Std.Dev: ", std(Data),"\n")
print("Prob. of Neg.: ", GetCertainty(Data,0,0),"\n")
println("")
print("p10, p50, p90 : ", quantile(collect(Float64, Data),[0.1,0.5,0.9]),"\n")
return chart
end

#SENSITIVITY CHART
"""
    sensitivity_chrt(ArrayName, TargetCol, Chrt_Type=1)

**TargetCol**: used to select the output against which the other variables are analyzed for influence.

**Chrt_Type**: allows to change the chart metric: Spearman (1), PPMC (2) and  Contribution to Variance % (3)
"""
function sensitivity_chrt(ArrayName, TargetCol, Chrt_Type=1)
      #Chrt_Type =="" || "please select a chart"
      # if TargetCol == nothing
      #       error("please select your target output")
      # end
      # if Chrt_Type == nothing
      #       Chrt_Type = 3
      # end
M_Size = size(ArrayName,2)
ArrayName = DataFrame(ArrayName)

#Calculate Spearman
cor_mat_s = []
for i=TargetCol
cor_vector = []
      for i2=1:M_Size
       cor_i = cor(tiedrank(ArrayName[!,i]),tiedrank(ArrayName[!,i2]));
       push!(cor_vector,cor_i)
       #print(cov_i)
      end
 push!(cor_mat_s,cor_vector)
end
correl_vals_s = hcat(cor_mat_s...)

#Calculate PPMC
cor_mat_p = []
for i=TargetCol
cor_vector = []
      for i2=1:M_Size
       cor_i = cor(ArrayName[!,i], ArrayName[!,i2]);
       push!(cor_vector,cor_i)
      end
push!(cor_mat_p,cor_vector)
end
correl_vals_p = hcat(cor_mat_p...)

#This removes all entries of 1 (self)
#correl_vals = correl_vals[(correl_vals[:].!=1)]

ystr = names(ArrayName)
if Chrt_Type == 1
      color_code = correl_vals_s .<0
elseif Chrt_Type == 2
      color_code = correl_vals_p .<0
elseif Chrt_Type == 3
      color_code = correl_vals_s .<0
else
      Chrt_Type = 3
      color_code = correl_vals_s .<0
end

impact=[]

for i in 1:size(correl_vals_s,1)
      if color_code[i] == true
            value="Negative"
      else
            value="Positive"
      end
      push!(impact, value)
end

var_sign=[]
for i in 1:size(correl_vals_s,1)
      if color_code[i] == true
            value = -1
      else
            value = 1
      end
      push!(var_sign, value)
end

#Contribution to Variance
cont_var = correl_vals_s .^2
cont_var = cont_var ./ (sum(cont_var)-1)
cont_var = cont_var .* var_sign

graph_tbl = DataFrame(hcat(ystr,correl_vals_s, abs.(correl_vals_s), correl_vals_p, cont_var, impact))
names!(graph_tbl, [:name, :correlation, :abs_cor, :PPMC, :cont_var, :impact])

#This removes all entries of 1 (self) in DataTable
graph_tbl = graph_tbl[graph_tbl[:correlation].!=1,:]

#Sort
graph_tbl = sort(graph_tbl,3,rev=false)

println(graph_tbl)

if Chrt_Type ==1
      return plot(graph_tbl, y=:name, x=:correlation, Guide.Title("Variables with Biggest Impact"), Guide.xlabel("Rank Correlation"), Guide.ylabel("Input"), Geom.bar(orientation=:horizontal), Theme(bar_spacing=10pt), color=:impact, Scale.color_discrete_manual(colorant"red", colorant"deepskyblue"), Scale.x_continuous(minvalue=-1, maxvalue=1))
elseif Chrt_Type ==2
      return plot(graph_tbl, y=:name, x=:PPMC, Guide.Title("Variables with Biggest Impact"), Guide.xlabel("Pearson Correlation"), Guide.ylabel("Input"), Geom.bar(orientation=:horizontal), Theme(bar_spacing=10pt), color=:impact, Scale.color_discrete_manual(colorant"red", colorant"deepskyblue"), Scale.x_continuous(minvalue=-1, maxvalue=1))
elseif Chrt_Type ==3
      return plot(graph_tbl, y=:name, x=:cont_var, Guide.Title("Variables with Biggest Impact"), Guide.xlabel("% Contribution to Variance"), Guide.ylabel("Input"), Geom.bar(orientation=:horizontal), Theme(bar_spacing=10pt), color=:impact, Scale.color_discrete_manual(colorant"red", colorant"deepskyblue"),Scale.x_continuous(minvalue=-1, maxvalue=1))
else
      return plot(graph_tbl, y=:name, x=:cont_var, Guide.Title("Variables with Biggest Impact"), Guide.xlabel("% Contribution to Variance"), Guide.ylabel("Input"), Geom.bar(orientation=:horizontal), Theme(bar_spacing=10pt), color=:impact, Scale.color_discrete_manual(colorant"red", colorant"deepskyblue"), Scale.x_continuous(minvalue=-1, maxvalue=1))

end
end

#TREND CHART
"""
    trend_chrt(SimTimeArray, PeriodRange, quantiles=[0.05,0.5,0.95])

**trend_chrt** allows the visualization of a simulated time series. These can be generated using the GBMM function.

**PeriodRange** must constructed using the `Dates` package and use the following syntax :

      dr = collect(Date(2019,1,01):Dates.Year(1):Date(2023,01,01))
"""
function trend_chrt(SimTimeArray, PeriodRange, quantiles=[0.05,0.5,0.95])
#In order to stack vector entries, use the mapping function ****

if typeof(SimTimeArray) == DataFrame
    AC_DF = SimTimeArray
else
AC_DF = vcat(map(x->x',SimTimeArray)...)
AC_DF = DataFrame(AC_DF)
end

trend_chart = []

for i = 1:size(AC_DF,2)
    push!(trend_chart, quantile!(AC_DF[!,i], quantiles))
end

#this allows us to join arrays generated by the quantile function
trend_chart = vcat(map(x->x',trend_chart)...)

#convert to time array to join dates and reconvert to DataFrame
trend_chart = TimeArray(PeriodRange, trend_chart)
trend_chart = DataFrame(trend_chart)

#Build Charts
names!(trend_chart, [:timestamp, :LowerBound, :p50, :UpperBound])
return plot(stack(trend_chart), y=:value, x=:timestamp, color=:variable, Geom.line)

end

#

#ADD Error Bar Example

# x = sort(rand(10)); xmin = 0.9x; xmax = 1.1x;
# ystr = map(i -> string(i) * " as a string", 1:10)
# plot(y = ystr, x = ex, xmin = xmin, xmax = xmax, Geom.errorbar)

#Generate an emprical s_curve
function s_curve(results)
    n=size(results,1)
    p = plot(y=collect((1:n)./n), x=sort(results), Geom.line)
    return p
end

# #My cumulative function for my Pareto magic
# test = sort(rand(Normal(10,1),1000))
# cum = cumsum(test)
# y_plot =cum ./ sum(test)
# plot(y=y_plot, x = test)


text:
module MyLib

include("../deps/deps.jl")
export delete_tms, get_tms_cstime, get_tms_cutime, get_tms_stime, get_tms_utime, libc_times, new_tms

function delete_tms(ptms)
    return ccall((:delete_tms, mylib), Void, (Ptr{Void},), ptms)
end

function get_tms_cstime(ptms)
    return ccall((:get_tms_cstime, mylib), Int32, (Ptr{Void},), ptms)
end

function get_tms_cutime(ptms)
    return ccall((:get_tms_cutime, mylib), Int32, (Ptr{Void},), ptms)
end

function get_tms_stime(ptms)
    return ccall((:get_tms_stime, mylib), Int32, (Ptr{Void},), ptms)
end

function get_tms_utime(ptms)
    return ccall((:get_tms_utime, mylib), Int32, (Ptr{Void},), ptms)
end

function libc_times(ptms)
    return ccall((:times, "libc"), Int32, (Ptr{Void},), ptms)
end

function new_tms()
    return ccall((:new_tms, mylib), Ptr{Void}, ())
end

end #module