text:
using Base64

"""
    run(doc::WeaveDoc; kwargs...)

Run code chunks and capture output from the parsed document.

## Keyword options

- `doctype::Union{Symbol,AbstractString} = :auto`: Output document format. `:auto` will set it automatically based on file extension. You can also manually specify it; see [`list_out_formats()`](@ref) for the supported formats
- `out_path::Union{Symbol,AbstractString} = :doc`: Path where the output is generated can be either of:
  * `:doc`: Path of the source document (default)
  * `:pwd`: Julia working directory
  * `"somepath"`: `String` of output directory e.g. `"~/outdir"`, or of filename e.g. `"~/outdir/outfile.tex"`
- `args::Dict = Dict()`: Arguments to be passed to the weaved document; will be available as `WEAVE_ARGS` in the document
- `mod::Union{Module,Symbol} = :sandbox`: Module where Weave `eval`s code. Defaults to `:sandbox` to create new sandbox module. You also can also pass a `Module` e.g. `Main`
- `fig_path::AbstractString = "figures"`: Where figures will be generated, relative to `out_path`
- `fig_ext::Union{Nothing,AbstractString} = nothing`: Extension for saved figures e.g. `".pdf"`, `".png"`. Default setting depends on `doctype`
- `cache_path::AbstractString = "cache"`: Where of cached output will be saved
- `cache::Symbol = :off`: Controls caching of code:
  * `:off` means no caching (default)
  * `:all` caches everything
  * `:user` caches based on chunk options
  * `:refresh` runs all code chunks and save new cache
- `throw_errors::Bool = false`: If `false` errors are included in output document and the whole document is executed. If `true` errors are thrown when they occur
- `latex_keep_unicode::Bool = false`: If `true`, do not convert unicode characters to their respective latex representation. This is especially useful if a font and tex-engine with support for unicode characters are used

!!! note
    Run Weave from terminal and try to avoid weaving from IJulia or ESS; they tend to mess with capturing output.
"""
function Base.run(
    doc::WeaveDoc;
    doctype::Union{Symbol,AbstractString} = :auto,
    out_path::Union{Symbol,AbstractString} = :doc,
    args::Dict = Dict(),
    mod::Union{Module,Symbol} = :sandbox,
    fig_path::AbstractString = "figures",
    fig_ext::Union{Nothing,AbstractString} = nothing,
    cache_path::AbstractString = "cache",
    cache::Symbol = :off,
    throw_errors::Bool = false,
    latex_keep_unicode::Bool = false
)
    #cache :all, :user, :off, :refresh

    doc.cwd = get_cwd(doc, out_path)
    doctype == :auto && (doctype = detect_doctype(doc.source))
    doc.doctype = doctype
    doc.format = formats[doctype]

    if (haskey(doc.format.formatdict, :keep_unicode))
        doc.format.formatdict[:keep_unicode] = latex_keep_unicode
    end

    isdir(doc.cwd) || mkpath(doc.cwd)

    if occursin("2pdf", doctype) && cache == :off
        fig_path = mktempdir(abspath(doc.cwd))
    elseif occursin("2html", doctype)
        fig_path = mktempdir(abspath(doc.cwd))
    end

    cache == :off || @eval import Serialization

    #This is needed for latex and should work on all output formats
    Sys.iswindows() && (fig_path = replace(fig_path, "\\" => "/"))

    doc.fig_path = fig_path
    set_rc_params(doc, fig_path, fig_ext)

    #New sandbox for each document with args exposed
    if mod == :sandbox
        sandbox = "WeaveSandBox$(rcParams[:doc_number])"
        mod = Core.eval(Main, Meta.parse("module $sandbox\nend"))
    end
    @eval mod WEAVE_ARGS = Dict()
    merge!(mod.WEAVE_ARGS, args)

    rcParams[:doc_number] += 1

    if haskey(doc.format.formatdict, :mimetypes)
      mimetypes = doc.format.formatdict[:mimetypes]
    else
      mimetypes = default_mime_types
    end

    report = Report(doc.cwd, doc.basename, doc.format.formatdict, mimetypes, throw_errors)
    pushdisplay(report)

    if cache != :off && cache != :refresh
        cached = read_cache(doc, cache_path)
        cached == nothing && @info("No cached results found, running code")
    else
        cached = nothing
    end

    executed = Any[]
    n = length(doc.chunks)

    for i = 1:n
        chunk = doc.chunks[i]

        if isa(chunk, CodeChunk)
            options = merge(doc.chunk_defaults, chunk.options)
            merge!(chunk.options, options)
        end

        restore = (cache ==:user && typeof(chunk) == CodeChunk && chunk.options[:cache])

        if cached != nothing && (cache == :all || restore)
            result_chunks = restore_chunk(chunk, cached)
        else
            result_chunks = run_chunk(chunk, doc, report, mod)
        end

        executed = [executed; result_chunks]
    end

    doc.header_script = report.header_script

    popdisplay(report)

    #Clear variables from used sandbox
    mod == :sandbox && clear_sandbox(SandBox)
    doc.chunks = executed

    if cache != :off
        write_cache(doc, cache_path)
    end

    return doc
end

"""
    detect_doctype(path::AbstractString)

Detect the output format based on file extension.
"""
function detect_doctype(path::AbstractString)
    _, ext = lowercase.(splitext(path))

    match(r"^\.(jl|.?md|ipynb)", ext) !== nothing && return "md2html"
    ext == ".rst" && return "rst"
    ext == ".tex" && return "texminted"
    ext == ".txt"  && return "asciidoc"

    return "pandoc"
end

function run_chunk(chunk::CodeChunk, doc::WeaveDoc, report::Report, SandBox::Module)
    @info("Weaving chunk $(chunk.number) from line $(chunk.start_line)")
    result_chunks = eval_chunk(chunk, report, SandBox)
    occursin("2html", report.formatdict[:doctype]) && (result_chunks = embed_figures(result_chunks, report.cwd))
    return result_chunks
end

function embed_figures(chunk::CodeChunk, cwd)
    chunk.figures = [img2base64(fig, cwd) for fig in chunk.figures]
    return chunk
end

function embed_figures(result_chunks, cwd)
    for i in 1:length(result_chunks)
        figs = result_chunks[i].figures
        if !isempty(figs)
            result_chunks[i].figures = [img2base64(fig, cwd) for fig in figs]
        end
    end
    return result_chunks
end

function img2base64(fig, cwd)
  ext = splitext(fig)[2]
  f = open(joinpath(cwd, fig), "r")
    raw = read(f)
  close(f)
  if ext == ".png"
    return "data:image/png;base64," * stringmime(MIME("image/png"), raw)
  elseif ext == ".svg"
    return "data:image/svg+xml;base64," * stringmime(MIME("image/svg"), raw)
  elseif ext == ".gif"
    return "data:image/gif;base64," * stringmime(MIME("image/gif"), raw)
  else
    return(fig)
  end
end

function run_chunk(chunk::DocChunk, doc::WeaveDoc, report::Report, SandBox::Module)
    chunk.content =  [run_inline(c, doc, report, SandBox) for c in chunk.content]
    return chunk
end

function run_inline(inline::InlineText, doc::WeaveDoc, report::Report, SandBox::Module)
    return inline
end

function run_inline(inline::InlineCode, doc::WeaveDoc, report::Report, SandBox::Module)
    #Make a temporary CodeChunk for running code. Collect results and don't wrap
    chunk = CodeChunk(inline.content, 0, 0, "", Dict(:hold => true, :wrap => false))
    options = merge(doc.chunk_defaults, chunk.options)
    merge!(chunk.options, options)

    chunks = eval_chunk(chunk, report, SandBox)
    occursin("2html", report.formatdict[:doctype]) && (chunks = embed_figures(chunks, report.cwd))

    output = chunks[1].output
    endswith(output, "\n") && (output = output[1:end-1])
    inline.output = output
    inline.rich_output = chunks[1].rich_output
    inline.figures = chunks[1].figures
    return inline
end

function reset_report(report::Report)
    report.cur_result = ""
    report.figures = AbstractString[]
    report.term_state = :text
end

function run_code(chunk::CodeChunk, report::Report, SandBox::Module)
    expressions = parse_input(chunk.content)
    N = length(expressions)
    #@show expressions
    result_no = 1
    results = ChunkOutput[ ]

    for (str_expr, expr) = expressions
        reset_report(report)
        lastline = (result_no == N)
        (obj, out) = capture_output(expr, SandBox, chunk.options[:term],
                      chunk.options[:display], lastline, report.throw_errors)
        figures = report.figures #Captured figures
        result = ChunkOutput(str_expr, out, report.cur_result, report.rich_output, figures)
        report.rich_output = ""
        push!(results, result)
        result_no += 1
    end
    return results
end

getstdout() = stdout

function capture_output(expr, SandBox::Module, term, disp,
                        lastline, throw_errors=false)
    #oldSTDOUT = STDOUT
    oldSTDOUT = getstdout()
    out = nothing
    obj = nothing
    rw, wr = redirect_stdout()
    reader = @async read(rw, String)
    try
        obj = Core.eval(SandBox, expr)
        if (term || disp) && (typeof(expr) != Expr || expr.head != :toplevel)
            obj != nothing && display(obj)
        #This shows images and lone variables, result can
        #Handle last line sepately
        elseif lastline && obj != nothing
            (typeof(expr) != Expr || expr.head != :toplevel) && display(obj)
        end
    catch E
        throw_errors && throw(E)
        display(E)
        @warn("ERROR: $(typeof(E)) occurred, including output in Weaved document")
    finally
        redirect_stdout(oldSTDOUT)
        close(wr)
        out = fetch(reader)
        close(rw)
    end
    out = replace(out, r"\u001b\[.*?m" => "") #Remove ANSI color codes
    return (obj, out)
end

#Parse chunk input to array of expressions
function parse_input(input::AbstractString)
    parsed = Tuple{AbstractString, Any}[]
    input = lstrip(input)
    n = sizeof(input)
    pos = 1 #The first character is extra line end
    while pos ≤ n
        oldpos = pos
        code,  pos = Meta.parse(input, pos)
        push!(parsed, (input[oldpos:pos-1] , code ))
    end
    parsed
end

function eval_chunk(chunk::CodeChunk, report::Report, SandBox::Module)
    if !chunk.options[:eval]
        chunk.output = ""
        chunk.options[:fig] = false
        return chunk
    end

    #Run preexecute_hooks
    for hook in preexecute_hooks
      chunk = Base.invokelatest(hook, chunk)
    end

    report.fignum = 1
    report.cur_chunk = chunk

    if haskey(report.formatdict, :out_width) && chunk.options[:out_width] == nothing
        chunk.options[:out_width] = report.formatdict[:out_width]
    end

    chunk.result = run_code(chunk, report, SandBox)

    #Run post_execute chunks
    for hook in postexecute_hooks
      chunk = Base.invokelatest(hook, chunk)
    end

    if chunk.options[:term]
        chunks = collect_results(chunk, TermResult())
    elseif chunk.options[:hold]
        chunks = collect_results(chunk, CollectResult())
    else
        chunks = collect_results(chunk, ScriptResult())
    end

      #else
     #   chunk.options[:fig] && (chunk.figures = copy(report.figures))
    #end

    chunks
end

#function eval_chunk(chunk::DocChunk, report::Report, SandBox)
#    chunk
#end

#Set all variables to nothing
function clear_sandbox(SandBox::Module)
    for name = names(SandBox, all=true)
        if name != :eval && name != names(SandBox)[1]
            try eval(SandBox, Meta.parse(AbstractString(AbstractString(name), "=nothing"))) catch; end
        end
    end
end

function get_figname(report::Report, chunk; fignum = nothing, ext = nothing)
    figpath = joinpath(report.cwd, chunk.options[:fig_path])
    isdir(figpath) || mkpath(figpath)
    ext == nothing && (ext = chunk.options[:fig_ext])
    fignum == nothing && (fignum = report.fignum)

    chunkid = (chunk.options[:label] == nothing) ? chunk.number : chunk.options[:label]
    full_name = joinpath(report.cwd, chunk.options[:fig_path],
        "$(report.basename)_$(chunkid)_$(fignum)$ext")
    rel_name = "$(chunk.options[:fig_path])/$(report.basename)_$(chunkid)_$(fignum)$ext" #Relative path is used in output
    return full_name, rel_name
end

function get_cwd(doc::WeaveDoc, out_path)
    #Set the output directory
    if out_path == :doc
        cwd = doc.path
    elseif out_path == :pwd
        cwd = pwd()
    else
        #If there is no extension, use as path
        splitted = splitext(out_path)
        if splitted[2] == ""
            cwd = expanduser(out_path)
        else
            cwd = splitdir(expanduser(out_path))[1]
        end
    end
    return cwd
end

"""Get output file name based on out_path"""
function get_outname(out_path::Symbol, doc::WeaveDoc; ext = nothing)
    ext == nothing && (ext = doc.format.formatdict[:extension])
    outname = "$(doc.cwd)/$(doc.basename).$ext"
end

"""Get output file name based on out_path"""
function get_outname(out_path::AbstractString, doc::WeaveDoc; ext = nothing)
    ext == nothing && (ext = doc.format.formatdict[:extension])
    splitted = splitext(out_path)
    if (splitted[2]) == ""
        outname = "$(doc.cwd)/$(doc.basename).$ext"
    else
        outname = expanduser(out_path)
    end
end

function set_rc_params(doc::WeaveDoc, fig_path, fig_ext)
    formatdict = doc.format.formatdict
    if fig_ext == nothing
        doc.chunk_defaults[:fig_ext] = formatdict[:fig_ext]
    else
        doc.chunk_defaults[:fig_ext] = fig_ext
    end
        doc.chunk_defaults[:fig_path] = fig_path
    return nothing
end

function collect_results(chunk::CodeChunk, fmt::ScriptResult)
    content = ""
    result_no = 1
    result_chunks = CodeChunk[ ]
    for r = chunk.result
        #Check if there is any output from chunk
        if strip(r.stdout) == "" && isempty(r.figures)  && strip(r.rich_output) == ""
            content *= r.code
        else
            content = "\n" * content * r.code
            rchunk = CodeChunk(content, chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
            content = ""
            rchunk.result_no = result_no
            result_no *=1
            rchunk.figures = r.figures
            rchunk.output = r.stdout * r.displayed
            rchunk.rich_output = r.rich_output
            push!(result_chunks, rchunk)
        end
    end
    if content != ""
        startswith(content, "\n") || (content = "\n" * content)
        rchunk = CodeChunk(content, chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_results(chunk::CodeChunk, fmt::TermResult)
    output = ""
    prompt = chunk.options[:prompt]
    result_no = 1
    result_chunks = CodeChunk[ ]
    for r = chunk.result
        output *= prompt * r.code
        output *= r.displayed * r.stdout
        if !isempty(r.figures)
            rchunk = CodeChunk("", chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
            rchunk.output = output
            output = ""
            rchunk.figures = r.figures
            push!(result_chunks, rchunk)
        end
    end
    if output != ""
         rchunk = CodeChunk("", chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
         rchunk.output = output
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_results(chunk::CodeChunk, fmt::CollectResult)
    result_no = 1
    for r =chunk.result
        chunk.output *=  r.stdout
        chunk.rich_output *= r.rich_output
        chunk.figures = [chunk.figures; r.figures]
    end
    return [chunk]
end

ext: jl


text:
function download_layer(w::WorldClim, layer::Integer)
    1 ≤ layer ≤ 19 || throw(ArgumentError("The layer must be between 1 and 19"))

    path = SimpleSDMLayers.assets_path()

    res = Dict(2.5 => "2.5", 5.0 => "5", 10.0 => "10")

    output_file = joinpath(path, "wc2.1_$(res[w.resolution])m_bio_$(layer).tif")
    zip_file = joinpath(path, "bioclim_2.1_$(res[w.resolution])m.zip")

    if !isfile(path)
        if !isfile(zip_file)
            root = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/"
            stem = "wc2.1_$(res[w.resolution])m_bio.zip"
            r = HTTP.request("GET", root * stem)
            open(zip_file, "w") do f
                write(f, String(r.body))
            end
        end
        zf = ZipFile.Reader(zip_file)
        file_to_read =
            first(filter(f -> joinpath(path, f.name) == output_file, zf.files))

        if !isfile(joinpath(path, file_to_read.name))
            write(joinpath(path, file_to_read.name), read(file_to_read))
        end
        close(zf)
    end

    return joinpath(path, file_to_read.name)
end

function download_layer(l::EarthEnv, layer::Integer)
    1 ≤ layer ≤ 12 || throw(ArgumentError("The layer must be between 1 and 12"))

    path = SimpleSDMLayers.assets_path()

    root = "https://data.earthenv.org/consensus_landcover/"
    stem = l.full ? "with_DISCover/consensus_full_class_$(layer).tif" :
        "without_DISCover/Consensus_reduced_class_$(layer).tif"
    filetype = l.full ? "complete" : "partial"
    filename = "landcover_$(filetype)_$(layer).tif"

    if !isfile(joinpath(path, filename))
        layerrequest = HTTP.request("GET", root * stem)
        open(joinpath(path, filename), "w") do layerfile
            write(layerfile, String(layerrequest.body))
        end
    end

    return joinpath(path, filename)
end

function download_layer(::BioClim, layer::Integer)
    1 ≤ layer ≤ 19 || throw(ArgumentError("The layer must be between 1 and 19"))
    path = SimpleSDMLayers.assets_path()
    layer = lpad(layer, 2, "0")
    filename = "CHELSA_bio10_$(layer).tif"
    url_root = "https://envicloud.os.zhdk.cloud.switch.ch/chelsa/chelsa_V1/climatologies/bio/"

    filepath = joinpath(path, filename)
    if !(isfile(filepath))
        res = HTTP.request("GET", url_root * filename)
        open(filepath, "w") do f
            write(f, String(res.body))
        end
    end
    return filepath
end


text:
# ------------------------------------------------------------------------------------------
# # Basics: Arithmetic, Assignment, and Strings
#
# Topics:
# 1. Infix arithmetic operators
# 2. Numeric literals
# 3. Comparisons
# 4. Assignments
# 5. Strings
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Arithmetic should largely be familiar
# ------------------------------------------------------------------------------------------

3 + 7 # addition

10 - 3 # subtraction

20 * 5 # multiplication

100 / 10 # division

10 ^ 2 # exponentiation

101 % 2 # remainder (modulus)

sqrt(2) # square root

√2 # Unicode to the rescue: \sqrt + TAB

# ------------------------------------------------------------------------------------------
# Note that dividing two integers yields a floating point number. There are two additional
# operators that may be helpful here:
# ------------------------------------------------------------------------------------------

10 / 6

10 ÷ 6 # \div TAB or the `div` function

div(10, 6)

10 // 6

# ------------------------------------------------------------------------------------------
# ### Numbers: Many different ways to write the number forty-two
# ------------------------------------------------------------------------------------------

fortytwos = (42, 42.0, 4.20e1, 4.20f1, 84//2, 0x2a)

for x in fortytwos
    show(x)
    println("\tisa $(typeof(x))")
end

# ------------------------------------------------------------------------------------------
# ### Bitwise arithmetic
# ------------------------------------------------------------------------------------------

0x2a & 0x70 # AND

0x2a | 0x70 # OR

42 & 112

0b0010 << 2 # == 0b1000

# ------------------------------------------------------------------------------------------
# ### Logical operators
# ------------------------------------------------------------------------------------------

false && true # AND

false || true # OR

# ------------------------------------------------------------------------------------------
# Note that they "short-circuit!"
# ------------------------------------------------------------------------------------------

x = -42
x > 0 || error("x must be positive")

# ------------------------------------------------------------------------------------------
# ### Comparisons
# ------------------------------------------------------------------------------------------

1 == 1.0 # Equality

1 === 1.0 # Programmatically identical

3 < π

1 <= 1

.1 + .2

.1 + .2 ≈ .3 # \approx + TAB

# ------------------------------------------------------------------------------------------
# Comparisons "chain"
#
# Try inserting parentheses around one of these comparisons
# ------------------------------------------------------------------------------------------

2 == 2.0 == 0x02

x = 42
0 < x < 100 || error("x must be between 0 and 100")

# ------------------------------------------------------------------------------------------
# # Assignment
#
# Assignment in Julia is done with the single `=`. All it does is associates a name (on the
# left) to a value (on the right).
# ------------------------------------------------------------------------------------------

x = 1 # Use the name `x` for the value `1`

y = x # Use the name `y` for the value that `x` refers to

x = "hello!" # Decide you have a better use for the name `x`

y # Is still the value 1

ϵ = eps(1.0) # You can make your own unicode names

5ϵ # Juxtaposition is multiplication

# ------------------------------------------------------------------------------------------
# ## Updating operators
#
# All the infix arithmetic operators above can be used as "updating" operators in
# conjunction with an assignment:
# ------------------------------------------------------------------------------------------

y += 1

# This is exactly the same as:
y = y + 1

# ------------------------------------------------------------------------------------------
# Note that it's just re-purposing the _same name_ for a new value. This means that the type
# might even change:
# ------------------------------------------------------------------------------------------

y /= 2

# ------------------------------------------------------------------------------------------
# # Strings
# ------------------------------------------------------------------------------------------

s1 = "I am a string."

s2 = """I am also a string. """

"Here, we get an "error" because it's ambiguous where this string ends "

"""Look, Mom, no "errors"!!! """

    println("""The other nice thing about triple-quoted
               string literals is that they ignore leading
               indentation, which is nice for long strings
               in real code. Try changing these quotes!""")

# ------------------------------------------------------------------------------------------
# Strings are not written with single `'`s — that's used for a single character:
# ------------------------------------------------------------------------------------------

first(s1)

'⊂'

'If you try writing a string in single-quotes, you will get an error'

# ------------------------------------------------------------------------------------------
# ## String interpolation
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# You can use the dollar sign inside a string to evaluate a Julia expression inside a string
# — either a single variable or a more complicated expression:
# ------------------------------------------------------------------------------------------

name = "Jane"
num_fingers = 10
num_toes = 10
println("Hello, my name is $name.")
println("I have $num_fingers fingers and $num_toes toes.")

 println("That is $(num_fingers + num_toes) digits in all!!")

@assert days == 365
@assert days_float == 365.0

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia/long-version/010 Arithmetic, Assignment, and Strings.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# We output a message after loading each package to work around a
# ten-minute timeout limit on Travis. Travis assumes the tests have hung
# if the interval between printing something to stdio exceeds ten minutes.
using Images
@info "Finished loading Images package."
using StaticArrays
@info "Finished loading StaticArrays package."

@testset "Error Evaluation" begin
    @info "Running flow evaluation test."

    img1 = load("../demo/car2.jpg")
    img2 = load("../demo/car1.jpg")

    algorithm = Farneback(50, estimation_window = 11,
                         σ_estimation_window = 9.0,
                         expansion_window = 6,
                         σ_expansion_window = 5.0)

    @time flow = optical_flow(Gray{Float32}.(img1), Gray{Float32}.(img2), algorithm)

    # Endpoint error statistics
    @time endpoint_error = evaluate_flow_error(flow, flow, EndpointError())
    @time ep_mean, ep_sd, ep_RX, ep_AX = calculate_statistics(endpoint_error)

    # Angular error statistics
    @time angular_error = evaluate_flow_error(flow, flow, AngularError())
    @time ang_mean, ang_sd, ang_RX, ang_AX = calculate_statistics(angular_error)

    @test ep_mean == 0
    @test ep_sd == 0
    @test ang_mean == 0
    @test ang_sd == 0
    @test ep_RX[0.5] == 0
    @test ep_AX[0.5] < 1e-5
    @test ang_RX[0.5] == 0
    @test ang_AX[0.5] < 1e-5
end;


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

using LinearAlgebra
using FFTW

import PlotlyJS
using Plots; plotly()

import OffsetArrays
import Interpolations

import BSON
import JSON

import Random
Random.seed!(25)

# simulation compounds. # paths.
plot_title = "resonance groupings"

# name of molecules to simulate.
# Limit to only one compound in this tutorial since we just want to visualize the different resonance groups for one compound.
#molecule_names = ["L-Histidine";]

tag = "L-Isoleucine"
molecule_names = [tag;]

# get mapping from molecule names to their spin system info json files.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"
dict_compound_to_filename = JSON.parsefile("/home/roy/Documents/repo/NMRData/input/compound_mapping/select_compounds.json")

# where the bson file is located.
root_folder = "/home/roy/MEGAsync/outputs/NMR/experiments/BMRB-500-0.5mM"
project_path = joinpath(root_folder, tag)
load_path = joinpath(project_path, "experiment.bson")

# where to save the resultant plot. Warning: This script will create the following path if it doesn't it exist.
save_folder = joinpath(project_path, "plots")

# spin-Hamiltonian-related.
tol_coherence = 1e-2 # resonances are pairs of eigenvalues of the Hamiltonian that have quantum numbers that differ by -1. This is the tolerance away from -1 that is allowed.
α_relative_threshold = 0.05 # resonances with relative amplitude less than this factor compared to the maximum resonance in the spin group will be removed. Set to 0.0 to see every single resonance component.
Δc_partition_radius = 0.17 # determines how many resonances get grouped together. Larger number means less groups and thus more resonances per group.
SH_config_path = "/home/roy/Documents/repo/NMRData/input/SH_configs/select_compounds_SH_configs.json"

# surrogate-related.
# default values. The alternative is to load from a config file.
Δr_default = 1.0 # the samples used to build the surrogate is taken every `Δr` radian on the frequency axis. Decrease for improved accuracy at the expense of computation resources.
Δκ_λ_default = 0.05 # the samples used to build thes urrogate for κ_λ are taken at this sampling spacing. Decrease for improved accuracy at the expense of computation resources.
Δcs_max_scalar_default = 0.2 # In units of ppm. interpolation border that is added to the lowest and highest resonance frequency component of the mixture being simulated.
κ_λ_lb_default = 0.5 # interpolation lower limit for κ_λ.
κ_λ_ub_default = 2.5 # interpolation upper limit for κ_λ.
λ0 = 3.4
surrogate_config_path = "/home/roy/Documents/repo/NMRData/input/surrogate_configs/select_compounds_surrogate_configs.json"

dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0) # level 2 model.

##### end user inputs.

### load block.
dic = BSON.load(load_path)
fs = dic[:fs]
SW = dic[:SW]
ν_0ppm = dic[:ν_0ppm]
λ_0ppm = dic[:λ_0ppm]

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

# get a surrogate where K_{n,i} is encouraged to be no larger than `early_exit_part_size`.
#println("Timing: setupmixtureproxies()")
#@time mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
mixture_params = NMRSpectraSimulator.setupmixtureSH(molecule_names,
    H_params_path, dict_compound_to_filename, fs, SW,
    ν_0ppm;
    config_path = SH_config_path,
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold,
    Δc_partition_radius = Δc_partition_radius)
As = mixture_params

# We only work with a single compound in this tutorial. Assign a new object for this compound to reduce clutter.
A = As[1];

## frequency locations. For plotting.
ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(mixture_params, hz2ppmfunc)
ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))

u_offset = 0.5 # in units ppm.
u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

# This is the frequency range that we shall work with.
P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)
U_rad = U .* (2*π)

# fit the surrogate.
Bs = NMRSpectraSimulator.fitproxies(As, dummy_SSFID, λ0;
    names = molecule_names,
    config_path = surrogate_config_path,
    Δcs_max_scalar_default = Δcs_max_scalar_default,
    κ_λ_lb_default = κ_λ_lb_default,
    κ_λ_ub_default = κ_λ_ub_default,
    u_min = u_min,
    u_max = u_max,
    Δr_default = Δr_default,
    Δκ_λ_default = Δκ_λ_default)
B = Bs[1] # looking at the first (and only) molecule in this example script.

# create the functions for each resonance group.
qs = collect( collect( ωω->B.qs[i][k](ωω-B.ss_params.d[i], B.ss_params.κs_λ[i]) for k = 1:length(B.qs[i]) ) for i = 1:length(B.qs) )
q_singlets = ωω->NMRSpectraSimulator.evalsinglets(ωω, B.d_singlets, A.αs_singlets, A.Ωs_singlets, B.β_singlets, B.λ0, B.κs_λ_singlets)

# create the function for the entire compound.
q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As[1:1], Bs[1:1])

# evaluate at the plotting positions.
q_U = q.(U_rad)

qs_U = collect( collect( qs[i][k].(U_rad) for k = 1:length(qs[i]) ) for i = 1:length(qs) )
q_singlets_U = q_singlets.(U_rad)

# sanity check.
q_check_U = q_singlets_U
if !isempty(qs) # some compounds only have singlets.
    q_check_U += sum( sum( qs[i][k].(U_rad) for k = 1:length(qs[i]) ) for i = 1:length(qs) )
end
discrepancy = norm(q_check_U- q_U)
println("sanity check. This should be numerically zero: ", discrepancy)

"""
save resonance groupings of a compound.
Choices for `f`` are: `real()`, `imag()`, or `abs()`. These corresponds to real part, imaginary part, and magnitude spectrum, respectively.
"""
function plotgroups(title_string::String,
    P,
    U,
    q,
    qs,
    q_singlets,
    f::Function,
    return_val_type::T;
    canvas_size::Tuple{Int,Int} = (1000, 400)) where T

    U_rad = U .* (2*π)

    q_U = q.(U_rad)
    plot_obj = Plots.plot( P,
        f.(q_U),
        title = title_string,
        label = "sum of sub-models",
        seriestype = :line,
        ticks = :native,
        xlims = (P[1],P[end]),
        hover = P,
        linewidth = 4,
        xflip = true,
        size = canvas_size)

    qs_U = Vector{Vector{Vector{Complex{T}}}}(undef, length(qs))
    for i = 1:length(qs)

        qs_U[i] = Vector{Vector{Complex{T}}}(undef, length(qs[i]))
        for k = 1:length(qs[i])

            qs_U[i][k] = qs[i][k].(U_rad)

            Plots.plot!(plot_obj, P, f.(qs_U[i][k]), label = "sub-group ($(i),$(k))",
                seriestype = :line,
                linestyle = :dot,
                xflip = true,
                linewidth = 4)
        end
    end

    q_singlets_U = q_singlets.(U_rad)
    Plots.plot!(plot_obj, P, f.(q_singlets_U), label = "group of all singlets",
        seriestype = :line,
        linestyle = :dot,
        xflip = true,
        linewidth = 4)
    #
    #Plots.plot!(plot_obj, P, f.(q_U),
    #markershape = :circle,
    #seriestype = :scatter,
    #xflip = true)

    return plot_obj, q_U, qs_U, q_singlets_U
end

# reduce the plotting positions for low signal regions. Otherwise the plot store size will be too large, and the time to load the plot will be long.
display_reduction_factor = 100
display_threshold_factor =  α_relative_threshold/10

inds, _ = NMRSpectraSimulator.prunelowsignalentries(q_U, display_threshold_factor, display_reduction_factor)
P_display = P[inds]
U_display = U[inds]

# plot.
canvas_size = (1000, 400)

plots_save_path = joinpath(save_folder, "groups_real.html")
title_string = "$(plot_title), real"
plot_obj, q_U, qs_U, q_singlets_U = plotgroups(title_string, P_display, U_display, q, qs, q_singlets, real, P[1]; canvas_size = canvas_size)
Plots.savefig(plot_obj, plots_save_path)
display(plot_obj)


text:
struct MbplsrSo
    fm
    T::Matrix{Float64}
    fit::Matrix{Float64}
    b
end

"""
    mbplsr_so(X, Y, weights = ones(size(X, 1)); nlv)
Multiblock sequentially orthogonalized PLSR (SO-PLSR).
* `X` : List (vector) of blocks (matrices) of X-data.
    Each component of the list is a block.
* `Y` : Y-data.
* `weights` : Weights of the observations (rows).
* `nlv` : Nb. latent variables (LVs) to consider
    for each block. Vector that must have a length equal to the nb. blocks.

Vector `weights` (row-weighting) is internally normalized to sum to 1.
See the help of `plskern` for details.

## References

- Biancolillo et al. , 2015. Combining SO-PLS and linear discriminant analysis
    for multi-block classification. Chemometrics and Intelligent Laboratory Systems,
    141, 58-67.

- Biancolillo, A. 2016. Method development in the area of multi-block analysis focused on
    food analysis. PhD. University of copenhagen.

- Menichelli et al., 2014. SO-PLS as an exploratory tool for path modelling.
    Food Quality and Preference, 36, 122-134.
"""
function mbplsr_so(X, Y, weights = ones(size(X[1], 1)); nlv)
    Y = ensure_mat(Y)
    n = size(X[1], 1)
    q = size(Y, 2)
    nbl = length(X)
    weights = mweights(weights)
    D = Diagonal(weights)
    fm = list(nbl)
    pred = similar(X[1], n, q)
    zX = copy(X)
    b = list(nbl)
    # First block
    fm[1] = plskern(zX[1], Y, weights; nlv = nlv[1])
    T = fm[1].T
    pred .= predict(fm[1], zX[1]).pred
    b[1] = nothing
    # Other blocks
    if nbl > 1
        for i = 2:nbl
            b[i] = inv(T' * (D * T)) * T' * (D * X[i])
            zX = X[i] - T * b[i]
            fm[i] = plskern(zX, Y - pred, weights; nlv = nlv[i])
            T = hcat(T, fm[i].T)
            pred .+= predict(fm[i], zX).pred
        end
    end
    MbplsrSo(fm, T, pred, b)
end

"""
    transform(object::MbplsrSo)
Compute LVs ("scores" T) from a fitted model.
* `object` : The maximal fitted model.
* `X` : A list (vector) of blocks (matrices) of X-data for which LVs are computed.
"""
function transform(object::MbplsrSo, X)
    nbl = length(object.fm)
    T = transform(object.fm[1], X[1])
    if nbl > 1
        @inbounds for i = 2:nbl
            zX = X[i] - T * object.b[i]
            T = hcat(T, transform(object.fm[i], zX))
        end
    end
    T
end

"""
    predict(object::MbplsrSo, X)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : A list (vector) of X-data for which predictions are computed.
"""
function predict(object::MbplsrSo, X)
    nbl = length(object.fm)
    T = transform(object.fm[1], X[1])
    pred =  object.fm[1].ymeans' .+ T * object.fm[1].C'
    if nbl > 1
        @inbounds for i = 2:nbl
            zX = X[i] - T * object.b[i]
            zT = transform(object.fm[i], zX)
            pred .+= object.fm[i].ymeans' .+ zT * object.fm[i].C'
            T = hcat(T, transform(object.fm[i], zX))
        end
    end
    (pred = pred,)
end


text:
using Test
using SHIPs, JuLIP, JuLIP.Testing, QuadGK, LinearAlgebra, SHIPs.JacobiPolys,
      BenchmarkTools
using SHIPs: TransformedJacobi, transform, transform_d, eval_basis!,
             alloc_B, alloc_temp, alloc_temp_d, alloc_dB, IntS,
             eval_basis_d!

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N), 0
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

##

Nmax = 4
rl, ru = 0.5, 3.0
fcut =  PolyCutoff2s(2, rl, ru)
trans = PolyTransform(2, 1.0)
spec = SparseSHIP(Nmax, 10)
shpB = SHIPBasis(spec, trans, fcut)

bgrp = shpB.bgrps[1]

##

Nr = 50
Rs, Zs = randR(Nr)
tmp = alloc_temp(shpB, Nr)
B = SHIPs.alloc_B(shpB)
tmp2 = alloc_temp(shpB2, Nr)
B2 = SHIPs.alloc_B(shpB2)

SHIPs.eval_basis!(B, tmp, shpB, Rs, Zs, 0)
SHIPs.eval_basis!(B2, tmp2, shpB2, Rs, Zs, 0)
@show B ≈ B2

@btime SHIPs.eval_basis!($B, $tmp, $shpB, $Rs, $Zs, 0);
@btime SHIPs.eval_basis!($B2, $tmp2, $shpB2, $Rs, $Zs, 0);

##

tmpd = alloc_temp_d(shpB, Nr)
dB = alloc_dB(shpB, Nr)
eval_basis_d!(B, dB, tmpd, shpB, Rs, Zs, 0)

tmpd2 = alloc_temp_d(shpB2, Nr)
dB2 = alloc_dB(shpB2, Nr)
eval_basis_d!(dB2, tmpd2, shpB2, Rs, Zs, 0)

@show dB2 ≈ dB

@btime eval_basis_d!($B, $dB, $tmpd, $shpB, $Rs, $Zs, 0)
@btime eval_basis_d!($dB2, $tmpd2, $shpB2, $Rs, $Zs, 0)

##

# ##
#
# f = let dB2 = dB2, tmpd2 = tmpd2, shpB2 = shpB2, Rs = Rs, Zs = Zs
#    () -> eval_basis_d!(dB2, tmpd2, shpB2, Rs, Zs, 0)
# end
# runn(f, N) = (for n=1:N; f(); end)
# runn(f, 1)
# ##
#
# using Profile
# Profile.clear()
# @profile runn(f, 100)
# Profile.print()
#
# ##
#
#
# # debugging -> shows that SHIPs.grad_AA_Rj! is correct...
#
# _AA(Rs) = (
#    SHIPs.precompute_A!(tmp2, shpB2, Rs, Zs, 1);
#    SHIPs.precompute_AA!(tmp2, shpB2, 1);
#    return copy(tmp2.AA[1])
#    )
#
# _dAA(Rs, j) = (
#    SHIPs.precompute_dA!(tmpd2, shpB2, Rs, Zs, 1);
#    SHIPs.precompute_AA!(tmpd2, shpB2, 1);
#    SHIPs.grad_AA_Rj!(tmpd2, shpB2, j, Rs, Zs, 1);
#    return copy(tmpd2.dAAj[1])
#    )
#
#
# u1 = rand(JVecF); u1 /= norm(u1)
# jj = 4
#
# for p = 2:12
#    h = 0.1^p
#
#    A = _AA(Rs)
#    Rs_h = copy(Rs); Rs_h[jj] += h * u1
#    A_h = _AA(Rs_h)
#    dA_h = (A_h - A) / h
#
#    dAA = _dAA(Rs, jj)
#    dAA_u1 = [ dot(u1, da) for da in dAA ]
#
#    @show norm(real.(dAA_u1 - dA_h), Inf)
# end
#
#

# h = 1e-5
#
# A = _AA(Rs)
# Rs_h = copy(Rs); Rs_h[1] += h * u1
# A_h = _AA(Rs_h)
# dA_h = (A_h - A) / h
#
# dAA = _dAA1(Rs)
# dAA_u1 = [ dot(da, u1) for da in dAA ]
#
# @show norm(real.(dAA_u1 - dA_h), Inf)
#
# @show imag.(dAA_u1 - dA_h)
#
#
#
# _A(Rs) = (
#    SHIPs.precompute_A!(tmp2, shpB2, Rs, Zs, 1);
#    return copy(tmp2.A[1])
#    )
#
# _dA1(Rs) = begin
#       SHIPs.precompute_dA!(tmpd2, shpB2, Rs, Zs, 1);
#       alist = shpB2.alists[1]
#       dA = zeros(JVec{ComplexF64}, length(alist))
#       for n = 1:length(alist)
#          zklm = alist[n]
#          dA[n] = SHIPs.grad_phi_Rj(Rs[1], 1, zklm, tmpd2)
#       end
#       return dA
#    end
#
# _dA1(Rs)
#
# for p = 2:12
#    h = 0.1^p
#
#    A = _A(Rs)
#    Rs_h = copy(Rs); Rs_h[1] += h * u1
#    A_h = _A(Rs_h)
#    dA_h = (A_h - A) / h
#
#    dA = _dA1(Rs)
#    dA_u1 = [ dot(u1, da) for da in dA ]
#
#    @show norm(imag.(dA_u1 - dA_h), Inf)
# end
#
#
# R = Rs[1]
# _phi(R) = _A([R])
# _dphi(R) = begin
#       SHIPs.precompute_dA!(tmpd2, shpB2, [R], [0], 1);
#       alist = shpB2.alists[1]
#       dphi = zeros(JVec{ComplexF64}, length(alist))
#       for n = 1:length(alist)
#          zklm = alist[n]
#          dphi[n] = SHIPs.grad_phi_Rj(R, 1, zklm, tmpd2)
#       end
#       return dphi
#    end
#
# for p = 2:12
#    h = 0.1^p
#
#    phi = _phi(R)
#    phi_h = _phi(R + h * u1)
#    dphi_h = (phi_h - phi) / h
#
#    dphi = dot.(Ref(u1), _dphi(R))
#
#    @show norm(imag.(dphi - dphi_h), Inf)
# end
#
# _phi(R)
# dot.(_dphi(R), Ref(u1))
#
#
#
# tmpd.J ≈ tmpd2.JJ
# tmpd.dJ ≈ tmpd2.dJJ
# tmpd.Y ≈ tmpd2.YY
# tmpd.dY ≈ tmpd2.dYY
#
#
# SH = shpB2.SH
# Y = copy(tmpd2.Y)
# dY = copy(tmpd2.dY)
#
# _Y(R) = SHIPs.eval_basis(SH, R)
# _dY(R, u) = dot.(Ref(u), SHIPs.eval_basis_d(SH, R)[2])
#
# h = 1e-5
# dY0 = _dY(R, u1)
# dYh = (_Y(R+h*u1) - _Y(R)) / h
# norm(imag.(dY0 - dYh), Inf)
# imag.(dY0 - dYh)

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
#=

In the following, ``uhat`` and ``tauhat`` are what the planner would
choose if he could reset at time t, ``uhatdif`` and ``tauhatdif`` are
the difference between those and what the planner is constrained to
choose.  The variable ``mu`` is the Lagrange multiplier associated with
the constraint at time t.

For more complete description of inputs and outputs see the website.

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

@date: 2014-08-21

References
----------

Simple port of the file examples/evans_sargent.py

http://quant-econ.net/hist_dep_policies.html

=#
using QuantEcon
using Roots
using Plots
using LaTeXStrings
pyplot()

struct HistDepRamsey{TF<:AbstractFloat}
    # These are the parameters of the economy
    A0::AbstractFloat
    A1::AbstractFloat
    d::AbstractFloat
    Q0::AbstractFloat
    tau0::AbstractFloat
    mu0::AbstractFloat
    bet::AbstractFloat

    # These are the LQ fields and stationary values
    R::Matrix{TF}
    A::Matrix{TF}
    B::Vector{TF}
    Q::TF
    P::Matrix{TF}
    F::Matrix{TF}
    lq::LQ
end

struct RamseyPath{TF<:AbstractFloat}
    y::Matrix{TF}
    uhat::Vector{TF}
    uhatdif::Vector{TF}
    tauhat::Vector{TF}
    tauhatdif::Vector{TF}
    mu::Vector{TF}
    G::Vector{TF}
    GPay::Vector{TF}
end

function HistDepRamsey(A0::AbstractFloat,
                       A1::AbstractFloat,
                       d::AbstractFloat,
                       Q0::AbstractFloat,
                       tau0::AbstractFloat,
                       mu::AbstractFloat,
                       bet::AbstractFloat)
    # Create Matrices for solving Ramsey problem
    R = [0.0  -A0/2  0.0    0.0
        -A0/2 A1/2   -mu/2  0.0
        0.0   -mu/2  0.0    0.0
        0.0   0.0    0.0    d/2]

    A = [1.0   0.0  0.0  0.0
         0.0   1.0  0.0  1.0
         0.0   0.0  0.0  0.0
         -A0/d A1/d 0.0  A1/d+1.0/bet]

    B = [0.0; 0.0; 1.0; 1.0/d]

    Q = 0.0

    # Use LQ to solve the Ramsey Problem.
    lq = LQ(Q, -R, A, B, bet=bet)

    P, F, _d = stationary_values(lq)

    HistDepRamsey(A0, A1, d, Q0, tau0, mu0, bet, R, A, B, Q, P, Array(F), lq)
end

function compute_G(hdr::HistDepRamsey, mu::AbstractFloat)
    # simplify notation
    Q0, tau0, A, B, Q = hdr.Q0, hdr.tau0, hdr.A, hdr.B, hdr.Q
    bet = hdr.bet

    R = hdr.R
    R[2, 3] = R[3, 2] = -mu/2
    lq = LQ(Q, -R, A, B, bet=bet)

    P, F, _d = stationary_values(lq)

    # Need y_0 to compute government tax revenue.
    u0 = compute_u0(hdr, P)
    y0 = vcat([1.0 Q0 tau0]', u0)

    # Define A_F and S matricies
    AF = A - B * F
    S = [0.0 1.0 0.0 0]' * [0.0 0.0 1.0 0]

    # Solves equation (25)
    Omega = solve_discrete_lyapunov(sqrt(bet) .* AF', bet * AF' * S * AF)
    T0 = y0' * Omega * y0

    return T0[1], A, B, F, P
end

function compute_u0(hdr::HistDepRamsey, P::Matrix)
    # simplify notation
    Q0, tau0 = hdr.Q0, hdr.tau0

    P21 = P[4, 1:3]
    P22 = P[4, 4]
    z0 = [1.0 Q0 tau0]'
    u0 = reshape(-P22^(-1).* P21,1,3)*(z0)

    return u0[1]
end

function init_path{TF<:AbstractFloat}(hdr::HistDepRamsey{TF}, mu0::TF, T::Integer=20)
    # Construct starting values for the path of the Ramsey economy
    G0, A, B, F, P = compute_G(hdr, mu0)

    # Compute the optimal u0
    u0 = compute_u0(hdr, P)

    # Initialize vectors
    y = Array{TF}(4, T)
    uhat       = Vector{TF}(T)
    uhatdif    = Vector{TF}(T)
    tauhat     = Vector{TF}(T)
    tauhatdif  = Vector{TF}(T)
    mu         = Vector{TF}(T)
    G          = Vector{TF}(T)
    GPay       = Vector{TF}(T)

    # Initial conditions
    G[1] = G0
    mu[1] = mu0
    uhatdif[1] = 0
    uhat[1] = u0
    tauhatdif[1] = 0
    y[:, 1] = vcat([1.0 hdr.Q0 hdr.tau0]', u0)

    return RamseyPath(y, uhat, uhatdif, tauhat, tauhatdif, mu, G, GPay)
end

function compute_ramsey_path!(hdr::HistDepRamsey, rp::RamseyPath)
    # simplify notation
    y, uhat, uhatdif, tauhat, = rp.y, rp.uhat, rp.uhatdif, rp.tauhat
    tauhatdif, mu, G, GPay = rp.tauhatdif, rp.mu, rp.G, rp.GPay
    bet = hdr.bet

    _, A, B, F, P = compute_G(hdr, mu[1])

    for t=2:T
        # iterate government policy
        y[:, t] = (A - B * F) * y[:, t-1]

        # update G
        G[t] = (G[t-1] - bet*y[2, t]*y[3, t])/bet
        GPay[t] = bet.*y[2, t]*y[3, t]

        #=
        Compute the mu if the government were able to reset its plan
        ff is the tax revenues the government would receive if they reset the
        plan with Lagrange multiplier mu minus current G
        =#
        ff(mu) = compute_G(hdr, mu)[1]-G[t]

        # find ff = 0
        mu[t] = fzero(ff, mu[t-1])
        _, Atemp, Btemp, Ftemp, Ptemp = compute_G(hdr, mu[t])

        # Compute alternative decisions
        P21temp = Ptemp[4, 1:3]
        P22temp = P[4, 4]
        uhat[t] = dot(-P22temp^(-1) .* P21temp, y[1:3, t])

        yhat = (Atemp-Btemp * Ftemp) * [y[1:3, t-1]; uhat[t-1]]
        tauhat[t] = yhat[4]
        tauhatdif[t] = tauhat[t] - y[4, t]
        uhatdif[t] = uhat[t] - y[4, t]
    end

    return rp
end

function plot1(rp::RamseyPath)
    tt = 0:length(rp.mu)  # tt is used to make the plot time index correct.
    y = rp.y

    ylabels = [L"$Q$" L"$\tau$" L"$u$"]
    #y_vals = [squeeze(y[2, :], 1) squeeze(y[3, :], 1) squeeze(y[4, :], 1)]
    y_vals = [y[2, :] y[3, :] y[4, :]]
    p = plot(tt, y_vals, color=:blue,
            label=["output" "tax rate" "first difference in output"],
            lw=2, alpha=0.7, ylabel=ylabels, layout=(3,1),
            xlims=(0, 15), xlabel=["" "" "time"], legend=:topright,
            xticks=0:5:15)
     return p
end

function plot2(rp::RamseyPath)
    y, uhatdif, tauhatdif, mu = rp.y, rp.uhatdif, rp.tauhatdif, rp.mu
    G, GPay = rp.G, rp.GPay
    T = length(rp.mu)
    tt = 0:T  # tt is used to make the plot time index correct.

    y_vals = [tauhatdif uhatdif mu G]
    ylabels = [L"$\Delta\tau$" L"$\Delta u$" L"$\mu$" L"$G$"]
    labels = hcat("time inconsistency differential for tax rate",
              L"time inconsistency differential for $u$",
              "Lagrange multiplier", "government revenue")
    p = plot(tt, y_vals, ylabel=ylabels, label=labels,
             layout=(4, 1), xlims=(-0.5, 15), lw=2, alpha=0.7,
             legend=:topright, color=:blue, xlabel=["" "" "" "time"])
    return p
end

# Primitives
T    = 20
A0   = 100.0
A1   = 0.05
d    = 0.20
bet = 0.95

# Initial conditions
mu0  = 0.0025
Q0   = 1000.0
tau0 = 0.0

# Solve Ramsey problem and compute path
hdr = HistDepRamsey(A0, A1, d, Q0, tau0, mu0, bet)
rp = init_path(hdr, mu0, T)
compute_ramsey_path!(hdr, rp)  # updates rp in place
p1=plot1(rp)
display(p1)
p2=plot2(rp)
display(p2)


text:
sz = (4,5)
filename = tempname()

# The mode "c" stands for creating a new file (clobber)

NCDatasets.NCDataset(filename,"c") do ds
    # define the dimension "lon" and "lat"
    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    v = NCDatasets.defVar(ds,"var1",UInt8,("lon","lat"))
    v[:,:] = fill(123,size(v))
end

# The "a" stands for appending to an existing file

NCDatasets.NCDataset(filename,"a") do ds
    v = NCDatasets.defVar(ds,"var2",UInt16,("lon","lat"))
    v[:,:] = fill(1234,size(v))

    # check if there are two variables
    @test length(keys(ds)) == 2
end

rm(filename)


text:
# DESCRIPTION
# Simulates open-loop dynamics using a dynamic inverse controller.
#
# INPUT
# - acfun: name of system dynamics function
# - finp: name of control input function
# - Tsim: total simulation time
# - dt: time step
# - x0: initial state vector
# - xdot0: initial state derivative vector
# - u0: initial control vector
#
# OUTPUT
# - state: state vector history
# - time: time vector history
#
# ------------------------------------------------------------------------------

function simulate(acfun,finp,Tsim,dt,x0,xdot0,u0)
    # set up simulation
    # vector of times
    time=collect(0:dt:Tsim)
    # initialize state time history
    state=zeros(NSTATES,length(time))
    dstate=zeros(NSTATES,length(time))
    # set initial condition
    state[:,1]=x0
    dstate[:,1]=xdot0
    # simulate open-loop system
    for i=1:length(time)-1
        # integrate
        sol=rk4(H60!,finp,time[i],dt,dstate[:,i],state[:,i],u0)
        state[:,i+1]=sol[1]
        dstate[:,i+1]=sol[2]
    end
    #
    return state, time
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
# Import `libgobject` and whatnot
using Glib_jll
import Base.eltype
#using Gtk.ShortNames # ShortNames to be fetched #meaningless #defined in GTK.jl's footer
#TODO: add _LList
# GType
abstract type GObject end
abstract type GInterface <: GObject end
abstract type GBoxed end
#= Infer: misleading ref
mutable struct GBoxedUnkown <: GBoxed #Gboxed
    handle::Ptr{GBoxed}
end
=#
const GEnum = Int32
const GType = Csize_t

struct GParamSpec
    g_type_instance::Ptr{Nothing}
    name::Ptr{UInt8}
    flags::Cint
    value_type::GType
    owner_type::GType
end

const fundamental_types = (
    #(:name,      Ctype,            JuliaType,      g_value_fn)
    (:invalid, Nothing, Union{}, :error),
    (:void, Nothing, Nothing, :error),
    (:GInterface, Ptr{Nothing}, GInterface, :error),
    (:gchar, Int8, Int8, :schar),
    (:guchar, UInt8, UInt8, :uchar),
    (:gboolean, Cint, Bool, :boolean),
    (:gint, Cint, Union{}, :int),
    (:guint, Cuint, Union{}, :uint),
    (:glong, Clong, Union{}, :long),
    (:gulong, Culong, Union{}, :ulong),
    (:gint64, Int64, Signed, :int64),
    (:guint64, UInt64, Unsigned, :uint64),
    (:GEnum, GEnum, Union{}, :enum),
    (:GFlags, GEnum, Union{}, :flags),
    (:gfloat, Float32, Float32, :float),
    (:gdouble, Float64, AbstractFloat, :double),
    (:gchararray, Ptr{UInt8}, AbstractString, :string),
    (:gpointer, Ptr{Nothing}, Ptr, :pointer),
    (:GBoxed, Ptr{GBoxed}, GBoxed, :boxed),
    (:GParam, Ptr{GParamSpec}, Ptr{GParamSpec}, :param),
    (:GObject, Ptr{GObject}, GObject, :object),
    #(:GVariant,  Ptr{GVariant},    GVariant,       :variant),
)
# NOTE: in general do not cache ids, except for these fundamental values
g_type_from_name(name::Symbol) = ccall((:g_type_from_name, libgobject), GType, (Ptr{UInt8},), name)
const fundamental_ids = tuple(GType[g_type_from_name(name) for (name, c, j, f) in fundamental_types]...)
#ERROR: libgobject is not defined
g_type(gtyp::GType) = gtyp
let jtypes = Expr(:block, :(g_type(::Type{Nothing}) = $(g_type_from_name(:void))))
    for i = 1:length(fundamental_types)
        (name, ctype, juliatype, g_value_fn) = fundamental_types[i]
        if juliatype != Union{}
            push!(jtypes.args, :(g_type(::Type{T}) where {T<:$juliatype} = convert(GType, $(fundamental_ids[i]))))
        end
    end
    Core.eval(GLib, jtypes) #GLib.jl requirement (GLib.jl or Glib_jll) #infer: another misleading trace
end
# NOTE: in general do not cache ids, except for these fundamental values
g_type_from_name(name::Symbol) = ccall((:g_type_from_name, libgobject), GType, (Ptr{UInt8},), name)
#=const=# fundamental_ids = tuple(GType[g_type_from_name(name) for (name, c, j, f) in fundamental_types]...)# TODO:define GType Info:Invalid redefinition of Constant #
#either remove constant, or remove redefinition() function ; redefinition proves to be costly, thus, action taken is to remove constant #TODO
#pre-supposes all fundamentals compile (without errors ) #requires error (file)

G_TYPE_FROM_CLASS(w::Ptr{Nothing}) = unsafe_load(convert(Ptr{GType}, w))
G_OBJECT_GET_CLASS(w::GObject) = G_OBJECT_GET_CLASS(w.handle)
G_OBJECT_GET_CLASS(hnd::Ptr{GObject}) = unsafe_load(convert(Ptr{Ptr{Nothing}}, hnd))
G_OBJECT_CLASS_TYPE(w) = G_TYPE_FROM_CLASS(G_OBJECT_GET_CLASS(w))

g_isa(gtyp::GType, is_a_type::GType) = ccall((:g_type_is_a, libgobject), Cint, (GType, GType), gtyp, is_a_type) != 0
g_isa(gtyp, is_a_type) = g_isa(g_type(gtyp), g_type(is_a_type))
g_type_parent(child::GType) = ccall((:g_type_parent, libgobject), GType, (GType,), child)
g_type_name(g_type::GType) = Symbol(bytestring(ccall((:g_type_name, libgobject), Ptr{UInt8}, (GType,), g_type)))

g_type_test_flags(g_type::GType, flag) = ccall((:g_type_test_flags, libgobject), Bool, (GType, GEnum), g_type, flag)
const G_TYPE_FLAG_CLASSED = 1 << 0
const G_TYPE_FLAG_INSTANTIATABLE = 1 << 1
const G_TYPE_FLAG_DERIVABLE = 1 << 2
const G_TYPE_FLAG_DEEP_DERIVABLE = 1 << 3
mutable struct GObjectLeaf <: GObject
    handle::Ptr{GObject}
    function GObjectLeaf(handle::Ptr{GObject})
        if handle == C_NULL
            error("Cannot construct $gname with a NULL pointer")
        end
        return gobject_ref(new(handle))
    end
end

g_type(obj::GObject) = g_type(typeof(obj))

gtypes(types...) = GType[g_type(t) for t in types]

const gtype_abstracts = Dict{Symbol,Type}()
const gtype_wrappers = Dict{Symbol,Type}()
const gtype_ifaces = Dict{Symbol,Type}()

gtype_abstracts[:GObject] = GObject # GObject #Definition  #abstract
gtype_wrappers[:GObject] = GObjectLeaf #GObject # Definition #leaf # concrete

let libs = Dict{AbstractString,Any}()
    global get_fn_ptr
    function get_fn_ptr(fnname, lib, cm)
        if !isa(lib, AbstractString)
            lib = Core.eval(cm, lib)
        end
        libptr = get(libs, lib, C_NULL)::Ptr{Nothing}
        if libptr == C_NULL
            libs[lib] = libptr = dlopen(lib)
        end
        fnptr = dlsym_e(libptr, fnname)
    end
end
function g_type(name::Symbol, lib, symname::Symbol, cm)
    if name in keys(gtype_wrappers)
        return g_type(gtype_wrappers[name])
    end
    fnptr = get_fn_ptr(string(symname, "_get_type"), lib, cm)
    if fnptr != C_NULL
        ccall(fnptr, GType, ())
    else
        convert(GType, 0)
    end
end
g_type(name::Symbol, lib, symname::Expr, cm) = Core.eval(cm, symname)
g_type(name::Expr, lib::Expr, symname::Expr, cm) = info((name, lib, symname))

function get_interface_decl(iname::Symbol, gtyp::GType, gtyp_decl, cm)
    if isdefined(cm, iname)
        return nothing
    end
    parent = g_type_parent(gtyp)
    @assert parent != 0
    piname = g_type_name(parent)
    quote
        if $(QuoteNode(iname)) in keys(gtype_ifaces)
            $(esc(iname)) = gtype_abstracts[$(Meta.quot(iname))]
        else
            struct $(esc(iname)) <: GInterface
                handle::Ptr{GObject}
                gc::Any
                $(esc(iname))(x::GObject) = new(unsafe_convert(Ptr{GObject}, x), x)
                # Gtk does an interface type check when calling methods. So, it's
                # not worth repeating it here. Plus, we might as well just allow
                # the user to lie, since we aren't using this for dispatch
                # (like C & unlike most other languages), the user may be able
                # to write more generic code
            end
            gtype_ifaces[$(QuoteNode(iname))] = $(esc(iname))
            local T #to prevent Julia-0.2 from name-mangling <: T
            $gtyp_decl
        end
        nothing
    end
end

function get_itype_decl(iname::Symbol, gtyp::GType, cm)
    if isdefined(cm, iname)
        return nothing
    end
    if iname === :GObject
        return :($(esc(iname)) = gtype_abstracts[:GObject])
    end
    #ntypes = mutable(Cuint)
    #interfaces = ccall((:g_type_interfaces, libgobject), Ptr{GType}, (GType, Ptr{Cuint}), gtyp, ntypes)
    #for i = 1:ntypes[]
    #    interface = unsafe_load(interfaces, i)
    #    # what do we care to do here?!
    #end
    #g_free(interfaces)
    parent = g_type_parent(gtyp)
    @assert parent != 0
    piname = g_type_name(parent)
    piface_decl = get_itype_decl(piname, parent, cm)
    quote
        if $(QuoteNode(iname)) in keys(gtype_abstracts)
            $(esc(iname)) = gtype_abstracts[$(QuoteNode(iname))]
        else
            $piface_decl
            abstract type $(esc(iname)) <: $(esc(piname)) end
            gtype_abstracts[$(QuoteNode(iname))] = $(esc(iname))
        end
        nothing
    end
end

get_gtype_decl(name::Symbol, lib, symname::Expr) =
    :(GLib.g_type(::Type{T}) where {T<:$(esc(name))} = $(esc(symname)))
let handled = Set()
    global get_gtype_decl
    function get_gtype_decl(name::Symbol, lib, symname::Symbol)
        if !(name in handled)
            push!(handled, name)
            return :(GLib.g_type(::Type{T}) where {T<:$(esc(name))} =
                ccall(($(QuoteNode(Symbol(string(symname, "_get_type")))), $(esc(lib))), GType, ()))
        end
        nothing
    end
end #let

function get_type_decl(name, iname, gtyp, gtype_decl, cm)
    ename = esc(name)
    einame = esc(iname)
    quote
        if $(QuoteNode(iname)) in keys(gtype_wrappers)
            $einame = gtype_abstracts[$(QuoteNode(iname))]
        else
            $(get_itype_decl(iname, gtyp, cm))
        end
        mutable struct $ename <: $einame
            handle::Ptr{GObject}
            function $ename(handle::Ptr{GObject})
                if handle == C_NULL
                    error($("Cannot construct $name with a NULL pointer"))
                end
                return gobject_ref(new(handle))
            end
        end
        local kwargs, T #to prevent Julia-0.2 from name-mangling kwargs, <: T
        function $ename(args...; kwargs...)
            if isempty(kwargs)
                error(MethodError($ename, args))
            end
            w = $ename(args...)
            for (kw, val) in kwargs
                set_gtk_property!(w, kw, val)
            end
            w
        end
        gtype_wrappers[$(QuoteNode(iname))] = $ename
        function $einame(args...; kwargs...)
            $ename(args...; kwargs...)
        end
        $gtype_decl
        nothing
    end
end

macro Gtype_decl(name, gtyp, gtype_decl)
    get_type_decl(name, Symbol(string(name, __module__.suffix)), gtyp, gtype_decl, __module__)
end

macro Gtype(iname, lib, symname)
    gtyp = g_type(iname, lib, symname, __module__)
    if gtyp == 0
        return Expr(:call, :error, string("Could not find ", symname, " in ", lib,
            ". This is likely a issue with a missing Gtk.jl version check."))
    end
    @assert iname === g_type_name(gtyp)
    if !g_type_test_flags(gtyp, G_TYPE_FLAG_CLASSED)
        error("GType is currently only implemented for G_TYPE_FLAG_CLASSED")
    end
    gtype_decl = get_gtype_decl(iname, lib, symname)
    name = Symbol(string(iname, __module__.suffix))
    get_type_decl(name, iname, gtyp, gtype_decl, __module__)
end
#=Error: no ref
macro Gabstract(iname, lib, symname)
    gtyp = g_type(iname, lib, symname) #Error: no ref
    if gtyp == 0
        return Expr(:call, :error, string("Could not find ", symname, " in ", lib, ". This is likely a issue with a missing Gtk.jl version check."))
    end
    @assert iname === g_type_name(gtyp)
    Expr(:block,
        get_itype_decl(iname, gtyp, __module__),
        get_gtype_decl(iname, lib, symname))
end
=#

macro Giface(iname, lib, symname)
    gtyp = g_type(iname, lib, symname, __module__)
    if gtyp == 0
        return Expr(:call, :error, string("Could not find ", symname, " in ", lib, ". This is likely a issue with a missing Gtk.jl version check."))
    end
    @assert iname === g_type_name(gtyp)
    gtype_decl = get_gtype_decl(iname, lib, symname)
    get_interface_decl(iname::Symbol, gtyp::GType, gtype_decl, __module__)
end

macro quark_str(q)
    :(ccall((:g_quark_from_string, libglib), UInt32, (Ptr{UInt8},), bytestring($q)))
end

unsafe_convert(::Type{Ptr{T}}, box::T) where {T<:GBoxed} = convert(Ptr{T}, box.handle)
convert(::Type{GBoxed}, boxed::GBoxed) = boxed
convert(::Type{GBoxedUnkown}, boxed::GBoxedUnkown) = boxed #TODO:Info an imported function has been exended without using module defined typed arguments

convert(::Type{T}, boxed::T) where {T<:GBoxed} = boxed
convert(::Type{T}, boxed::GBoxed) where {T<:GBoxed} = convert(T, boxed.handle)
convert(::Type{GBoxed}, unbox::Ptr{GBoxed}) = GBoxedUnkown(unbox)
convert(::Type{GBoxed}, unbox::Ptr{T}) where {T<:GBoxed} = GBoxedUnkown(unbox)
convert(::Type{T}, unbox::Ptr{GBoxed}) where {T<:GBoxed} = convert(T, convert(Ptr{T}, unbox))
convert(::Type{T}, unbox::Ptr{T}) where {T<:GBoxed} = T(unbox)

cconvert(::Type{Ptr{GObject}}, @nospecialize(x::GObject)) = x

# All GObjects are expected to have a 'handle' field
# of type Ptr{GObject} corresponding to the GLib object
# or to override this method (e.g. GtkNullContainer, AbstractString)
unsafe_convert(::Type{Ptr{GObject}}, w::GObject) = getfield(w, :handle)

# this method should be used by gtk methods returning widgets of unknown type
# and/or that might have been wrapped by julia before,
# instead of a direct call to the constructor
convert(::Type{T}, w::Ptr{GObject}) where {T<:GObject} = convert_(T, convert(Ptr{T}, w)) # this definition must be first due to a 0.2 dispatch bug
convert(::Type{T}, ptr::Ptr{T}) where {T<:GObject} = convert_(T, ptr)

# need to introduce convert_ since otherwise there was a StackOverFlow error
function convert_(::Type{T}, ptr::Ptr{T}) where {T<:GObject}
    hnd = convert(Ptr{GObject}, ptr)
    if hnd == C_NULL
        throw(UndefRefError())
    end
    x = ccall((:g_object_get_qdata, libgobject), Ptr{GObject}, (Ptr{GObject}, UInt32), hnd, jlref_quark::UInt32)
    if x != C_NULL
        return gobject_ref(unsafe_pointer_to_objref(x)::T)
    end
    wrap_gobject(hnd)::T
end

function wrap_gobject(hnd::Ptr{GObject})
    gtyp = G_OBJECT_CLASS_TYPE(hnd)
    typname = g_type_name(gtyp)
    while !(typname in keys(gtype_wrappers))
        gtyp = g_type_parent(gtyp)
        @assert gtyp != 0
        typname = g_type_name(gtyp)
    end
    T = gtype_wrappers[typname]
    return T(hnd)
end

eltype(::Type{_LList{T}}) where {T<:GObject} = T #_LList is not defined
ref_to(::Type{T}, x) where {T<:GObject} = gobject_ref(unsafe_convert(Ptr{GObject}, x))
deref_to(::Type{T}, x::Ptr) where {T<:GObject} = convert(T, x)
empty!(li::Ptr{_LList{Ptr{T}}}) where {T<:GObject} = gc_unref(unsafe_load(li).data)
#_LList not defined
### Miscellaneous types
baremodule GConnectFlags
const AFTER = 1
const SWAPPED = 2
get(s::Symbol) =
    if s === :after
        AFTER
    elseif s === :swapped
        SWAPPED
    else
        Main.Base.error(Main.Base.string("invalid GConnectFlag ", s))
    end
end

### Garbage collection [prevention]
const gc_preserve = IdDict{Any,Any}() # reference counted closures
function gc_ref(@nospecialize(x))
    global gc_preserve
    local ref::Ref{Any}, cnt::Int
    if x in keys(gc_preserve)
        ref, cnt = gc_preserve[x]::Tuple{Ref{Any},Int}
    else
        ref = Ref{Any}(x)
        cnt = 0
    end
    gc_preserve[x] = (ref, cnt + 1)
    return unsafe_load(convert(Ptr{Ptr{Nothing}}, unsafe_convert(Ptr{Any}, ref)))
end
function gc_unref(@nospecialize(x))
    global gc_preserve
    ref, cnt = gc_preserve[x]::Tuple{Ref{Any},Int}
    @assert cnt > 0
    if cnt == 1
        delete!(gc_preserve, x)
    else
        gc_preserve[x] = (ref, cnt - 1)
    end
    nothing
end
_gc_unref(@nospecialize(x), ::Ptr{Nothing}) = gc_unref(x)
gc_ref_closure(@nospecialize(cb::Function)) = (invoke(gc_ref, Tuple{Any}, cb), @cfunction(_gc_unref, Nothing, (Any, Ptr{Nothing})))
gc_ref_closure(x::T) where {T} = (gc_ref(x), @cfunction(_gc_unref, Nothing, (Any, Ptr{Nothing})))

# generally, you shouldn't be calling gc_ref(::Ptr{GObject})
gc_ref(x::Ptr{GObject}) = ccall((:g_object_ref, libgobject), Nothing, (Ptr{GObject},), x)
gc_unref(x::Ptr{GObject}) = ccall((:g_object_unref, libgobject), Nothing, (Ptr{GObject},), x)

const gc_preserve_glib = Dict{Union{WeakRef,GObject},Bool}() # glib objects
const gc_preserve_glib_lock = Ref(false) # to satisfy this lock, must never decrement a ref counter while it is held
const topfinalizer = Ref(true) # keep recursion to a minimum by only iterating from the top
const await_finalize = Set{Any}()

Base.isequal(x::GObject, w::WeakRef) = x === w.value   # cuts the number of MethodInstances from O(N^2) to O(N)

function finalize_gc_unref(@nospecialize(x::GObject))
    # this records that the are no user references left to the object from Julia
    # and notifies GLib that it can free the object (if no reference exist from C)
    # it is intended to be called by GC, not in user code function
    istop = topfinalizer[]
    topfinalizer[] = false
    gc_preserve_glib_lock[] = true
    delete!(gc_preserve_glib, x)
    if getfield(x, :handle) != C_NULL
        gc_preserve_glib[x] = true # convert to a strong-reference
        gc_preserve_glib_lock[] = false
        gc_unref(unsafe_convert(Ptr{GObject}, x)) # may clear the strong reference
    else
        gc_preserve_glib_lock[] = false
    end
    topfinalizer[] = istop
    istop && run_delayed_finalizers()
    nothing
end

function delref(@nospecialize(x::GObject))
    # internal helper function
    exiting[] && return # unnecessary to cleanup if we are about to die anyways
    if gc_preserve_glib_lock[] || g_yielded[]
        push!(await_finalize, x)
        return # avoid running finalizers at random times
    end
    finalize_gc_unref(x)
    nothing
end
function addref(@nospecialize(x::GObject))
    # internal helper function
    ccall((:g_object_ref_sink, libgobject), Ptr{GObject}, (Ptr{GObject},), x)
    finalizer(delref, x)
    delete!(gc_preserve_glib, x) # in v0.2, the WeakRef assignment below wouldn't update the key
    gc_preserve_glib[WeakRef(x)] = false # record the existence of the object, but allow the finalizer
    nothing
end
function gobject_ref(x::T) where {T<:GObject}
    gc_preserve_glib_lock[] = true
    strong = get(gc_preserve_glib, x, nothing)
    if strong === nothing
        if ccall((:g_object_get_qdata, libgobject), Ptr{Cvoid},
            (Ptr{GObject}, UInt32), x, jlref_quark::UInt32) != C_NULL
            # have set up metadata for this before, but its weakref has been cleared. restore the ref.
            delete!(await_finalize, x)
            finalizer(delref, x)
            gc_preserve_glib[WeakRef(x)] = false # record the existence of the object, but allow the finalizer
        else
            # we haven't seen this before, setup the metadata
            deref = @cfunction(gc_unref, Nothing, (Ref{T},))
            ccall((:g_object_set_qdata_full, libgobject), Nothing,
                (Ptr{GObject}, UInt32, Any, Ptr{Nothing}), x, jlref_quark::UInt32, x,
                deref) # add a circular reference to the Julia object in the GObject
            addref(Ref{GObject}(x)[])
        end
    elseif strong
        # oops, we previously deleted the link, but now it's back
        addref(Ref{GObject}(x)[])
    else
        # already gc-protected, nothing to do
    end
    gc_preserve_glib_lock[] = false
    run_delayed_finalizers()
    return x
end

gc_ref(x::GObject) = pointer_from_objref(gobject_ref(x))

function run_delayed_finalizers()
    exiting[] && return # unnecessary to cleanup if we are about to die anyways
    g_yielded[] && return # can't run them right now
    topfinalizer[] = false
    while !isempty(await_finalize)
        x = pop!(await_finalize)
        finalize_gc_unref(x)
    end
    topfinalizer[] = true
end

function gc_unref_weak(x::GObject)
    # this strongly destroys and invalidates the object
    # it is intended to be called by GLib, not in user code function
    # note: this may be called multiple times by GLib
    x.handle = C_NULL
    gc_preserve_glib_lock[] = true
    delete!(gc_preserve_glib, x)
    gc_preserve_glib_lock[] = false
    nothing
end

function gc_unref(x::GObject)
    # this strongly destroys and invalidates the object
    # it is intended to be called by GLib, not in user code function
    ref = ccall((:g_object_get_qdata, libgobject), Ptr{Nothing}, (Ptr{GObject}, UInt32), x, jlref_quark::UInt32)
    if ref != C_NULL && x !== unsafe_pointer_to_objref(ref)
        # We got called because we are no longer the default object for this handle, but we are still alive
        @warn("Duplicate Julia object creation detected for GObject")
        deref = cfunction_(gc_unref_weak, Nothing, (Ref{typeof(x)},))
        ccall((:g_object_weak_ref, libgobject), Nothing, (Ptr{GObject}, Ptr{Nothing}, Any), x, deref, x)
    else
        ccall((:g_object_steal_qdata, libgobject), Any, (Ptr{GObject}, UInt32), x, jlref_quark::UInt32)
        gc_unref_weak(x)
    end
    nothing
end

gc_unref(::Ptr{GObject}, x::GObject) = gc_unref(x)
gc_ref_closure(x::GObject) = (gc_ref(x), C_NULL)

function gc_force_floating(x::GObject)
    ccall((:g_object_force_floating, libgobject), Nothing, (Ptr{GObject},), x)
end
function gobject_move_ref(new::GObject, old::GObject)
    h = unsafe_convert(Ptr{GObject}, new)
    @assert h == unsafe_convert(Ptr{GObject}, old) != C_NULL
    gc_ref(h)
    gc_unref(old)
    gc_ref(new)
    gc_unref(h)
    new
end

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
using LightXML

struct JointDamping{T}
    coeffs::Dict{JointID, Vector{T}}
end

function JointDamping{T}(mechanism::Mechanism, urdf::AbstractString) where T
    xdoc = parse_file(urdf)
    xroot = LightXML.root(xdoc)
    @assert LightXML.name(xroot) == "robot"
    coeffs = Vector{Pair{JointID, Vector{T}}}()
    xml_joints = get_elements_by_tagname(xroot, "joint")
    for xml_joint in xml_joints
        name = attribute(xml_joint, "name")
        try
            joint = findjoint(mechanism, name)
            xml_dynamics = find_element(xml_joint, "dynamics")
            if xml_dynamics != nothing
                xml_damping = attribute(xml_dynamics, "damping")
                if xml_damping != nothing
                    push!(coeffs, JointID(joint) => [parse(T, xml_damping)])
                end
            end
        catch # TODO: be more specific
        end
    end
    return JointDamping(Dict(coeffs))
end

function (damping::JointDamping)(τ::AbstractVector, t, state::MechanismState)
    τ .= 0
    v = velocity(state)
    @inbounds for (jointid, coeffs) in damping.coeffs
        τjoint = τ[jointid]
        vjoint = v[jointid]
        τjoint .= .-coeffs .* vjoint
    end
    return τ
end


text:
module FrostHelperCapDashOnGroundTriggerModule

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/CapDashOnGroundTrigger" FrostHelperCapDashOnGroundTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16)

const placements = Ahorn.PlacementDict(
	"Cap Dash Count On Ground (Frost Helper)" => Ahorn.EntityPlacement(
		FrostHelperCapDashOnGroundTrigger,
		"rectangle"
	)
)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


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

mutable struct OutboundRulePropertiesFormat <: SwaggerModel
    allocatedOutboundPorts::Any # spec type: Union{ Nothing, Int32 } # spec name: allocatedOutboundPorts
    frontendIPConfigurations::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: frontendIPConfigurations
    backendAddressPool::Any # spec type: Union{ Nothing, SubResource } # spec name: backendAddressPool
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    enableTcpReset::Any # spec type: Union{ Nothing, Bool } # spec name: enableTcpReset
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes

    function OutboundRulePropertiesFormat(;allocatedOutboundPorts=nothing, frontendIPConfigurations=nothing, backendAddressPool=nothing, provisioningState=nothing, protocol=nothing, enableTcpReset=nothing, idleTimeoutInMinutes=nothing)
        o = new()
        validate_property(OutboundRulePropertiesFormat, Symbol("allocatedOutboundPorts"), allocatedOutboundPorts)
        setfield!(o, Symbol("allocatedOutboundPorts"), allocatedOutboundPorts)
        validate_property(OutboundRulePropertiesFormat, Symbol("frontendIPConfigurations"), frontendIPConfigurations)
        setfield!(o, Symbol("frontendIPConfigurations"), frontendIPConfigurations)
        validate_property(OutboundRulePropertiesFormat, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(OutboundRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(OutboundRulePropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(OutboundRulePropertiesFormat, Symbol("enableTcpReset"), enableTcpReset)
        setfield!(o, Symbol("enableTcpReset"), enableTcpReset)
        validate_property(OutboundRulePropertiesFormat, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        o
    end
end # type OutboundRulePropertiesFormat

const _property_map_OutboundRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("allocatedOutboundPorts")=>Symbol("allocatedOutboundPorts"), Symbol("frontendIPConfigurations")=>Symbol("frontendIPConfigurations"), Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("protocol")=>Symbol("protocol"), Symbol("enableTcpReset")=>Symbol("enableTcpReset"), Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"))
const _property_types_OutboundRulePropertiesFormat = Dict{Symbol,String}(Symbol("allocatedOutboundPorts")=>"Int32", Symbol("frontendIPConfigurations")=>"Vector{SubResource}", Symbol("backendAddressPool")=>"SubResource", Symbol("provisioningState")=>"ProvisioningState", Symbol("protocol")=>"String", Symbol("enableTcpReset")=>"Bool", Symbol("idleTimeoutInMinutes")=>"Int32")
Base.propertynames(::Type{ OutboundRulePropertiesFormat }) = collect(keys(_property_map_OutboundRulePropertiesFormat))
Swagger.property_type(::Type{ OutboundRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OutboundRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ OutboundRulePropertiesFormat }, property_name::Symbol) =  _property_map_OutboundRulePropertiesFormat[property_name]

const _allowed_OutboundRulePropertiesFormat_protocol = ["Tcp", "Udp", "All"]

function check_required(o::OutboundRulePropertiesFormat)
    (getproperty(o, Symbol("frontendIPConfigurations")) === nothing) && (return false)
    (getproperty(o, Symbol("backendAddressPool")) === nothing) && (return false)
    true
end

function validate_property(::Type{ OutboundRulePropertiesFormat }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "OutboundRulePropertiesFormat", :enum, val, _allowed_OutboundRulePropertiesFormat_protocol)
    end
end


text:
using Divergences
using Base.Test

# write your own tests here
@test 1 == 1

cr = CressieRead(2)

srand(1)
v = rand(10)
v = 10*v./sum(v)

x = rand(10)
x = 10*x./sum(x)

evaluate(cr, v)
evaluate(cr, v, x)

gradient(cr, 1.1)
gradient(cr, 0.9)

hessian(cr, 1.0)
hessian(cr, 1.0)

################################################################################
##
## Function - Cressie Read
##
################################################################################

seq =0:.1:3
cr = CressieRead(2.)

ch =[0.333333, 0.2835, 0.234667, 0.187833, 0.144, 0.104167, 0.0693333, 0.0405,
	 0.0186667, 0.00483333, 0., 0.00516667, 0.0213333, 0.0495, 0.0906667,
	 0.145833, 0.216, 0.302167, 0.405333, 0.5265, 0.666667, 0.826833, 1.008,
	 1.21117, 1.43733, 1.6875, 1.96267, 2.26383, 2.592, 2.94817, 3.33333]

[@test_approx_eq_eps ch[index] evaluate(cr, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

cr = CressieRead(.5)

ch = [0.666667, 0.50883, 0.385924, 0.285756, 0.203976, 0.138071, 0.086344, 0.0475494,
	 0.0207223, 0.00508662, 0., 0.00491964, 0.0193789, 0.0429707, 0.0753365, 0.116156,
	 0.165144, 0.222038, 0.286605, 0.358626, 0.437903, 0.524252, 0.617503, 0.717497,
	 0.824085, 0.937129, 1.0565, 1.18207, 1.31373, 1.45136, 1.59487];

[@test_approx_eq_eps ch[index] evaluate(cr, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

cr = CressieRead(-.5)

ch = [2., 0.935089, 0.611146, 0.40911, 0.270178, 0.171573, 0.101613, 0.0533599,
	  0.0222912, 0.00526681, 0., 0.00476461, 0.0182195, 0.0392983, 0.0671362,
	  0.101021, 0.140356, 0.184638, 0.233437, 0.28638, 0.343146, 0.403449,
	  0.467041, 0.5337, 0.603227, 0.675445, 0.750194, 0.827329, 0.90672, 0.988245, 1.0718]

[@test_approx_eq_eps ch[index] evaluate(cr, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

################################################################################
##
## Gradient - Cressie Read
##
################################################################################

cr = CressieRead(2.)

ch = [-0.5, -0.495, -0.48, -0.455, -0.42, -0.375, -0.32, -0.255, -0.18,
	  -0.095, 0., 0.105, 0.22, 0.345, 0.48, 0.625, 0.78, 0.945, 1.12,
	  1.305, 1.5, 1.705, 1.92, 2.145, 2.38, 2.625, 2.88, 3.145, 3.42, 3.705, 4.]

[@test_approx_eq_eps ch[index] gradient(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

cr = CressieRead(.5)

ch = [-2., -1.36754, -1.10557, -0.904555, -0.735089, -0.585786, -0.450807,
-0.32668, -0.211146, -0.102633, 0., 0.0976177, 0.19089, 0.280351, 0.366432,
0.44949, 0.529822, 0.607681, 0.683282, 0.75681, 0.828427, 0.898275,
0.966479, 1.03315, 1.09839, 1.16228, 1.2249, 1.28634,  1.34664, 1.40588, 1.4641]

[@test_approx_eq_eps ch[index] gradient(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

cr = CressieRead(-.5)

ch = [-Inf, -4.32456, -2.47214, -1.65148, -1.16228, -0.828427, -0.581989, -0.390457,
	  -0.236068, -0.108185, 0., 0.0930748, 0.174258, 0.245884, 0.309691,
	  0.367007, 0.418861, 0.46607, 0.509288, 0.549047, 0.585786, 0.619869,
	  0.6516, 0.681239, 0.709006, 0.735089, 0.759653, 0.782839, 0.804771,
	  0.82556, 0.84529]

[@test_approx_eq_eps ch[index] gradient(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

###############################################################################
##
## Hessian - Cressie Read
##
#############################################################################
cr = CressieRead(2.)
ch = [0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4,
   1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.]
[@test_approx_eq_eps ch[index] hessian(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

cr = CressieRead(-.5)
ch = [+Inf, 31.6228, 11.1803, 6.08581, 3.95285, 2.82843, 2.15166, 1.70747,
      1.39754, 1.17121, 1., 0.866784, 0.760726, 0.67466, 0.603682, 0.544331,
      0.494106, 0.451156, 0.414087, 0.38183, 0.353553, 0.328603, 0.306454,
      0.286687, 0.268957, 0.252982, 0.238528, 0.2254, 0.213434, 0.20249, 0.19245]
[@test_approx_eq_eps ch[index] hessian(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

###############################################################################
##
## Function - Modified Cressie Read
##
#############################################################################

cr = ModifiedCressieRead(2, .2)

ch = [0.333333, 0.2835, 0.234667, 0.187833, 0.144, 0.104167, 0.0693333, 0.0405,
	  0.0186667, 0.00483333, 0., 0.00516667, 0.0213333, 0.0493333, 0.0893333,
	  0.141333, 0.205333, 0.281333, 0.369333, 0.469333, 0.581333, 0.705333,
	  0.841333, 0.989333, 1.14933, 1.32133, 1.50533, 1.70133, 1.90933, 2.12933, 2.36133]

[@test_approx_eq_eps ch[index] evaluate(cr, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

cr = ModifiedCressieRead(-.5, .2)

ch = [2., 0.935089, 0.611146, 0.40911, 0.270178, 0.171573, 0.101613, .0533599,
      0.0222912, 0.00526681, 0., 0.00476461, 0.0182195, 0.039449, 0.0682857,
      0.10473, 0.148781, 0.200439, 0.259705, 0.326578, 0.401058, 0.483146,
      0.572841, 0.670143, 0.775052, 0.887568, 1.00769, 1.13542, 1.27076,
      1.41371, 1.56426]
[@test_approx_eq_eps ch[index] evaluate(cr, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

################################################################################
##
## Gradient - Modified Cressie Read
##
################################################################################

seq =0:.101:3
cr = ModifiedCressieRead(2., .2)

	ch = [-0.5, -0.494899, -0.479598, -0.454096, -0.418392, -0.372487, -0.316382,
	      -0.250075, -0.173568, -0.0868595, 0.01005, 0.117161, 0.2344, 0.3556,
	      0.4768, 0.598, 0.7192, 0.8404, 0.9616, 1.0828, 1.204, 1.3252, 1.4464,
	      1.5676, 1.6888, 1.81, 1.9312, 2.0524, 2.1736, 2.2948]
[@test_approx_eq_eps ch[index] gradient(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

cr = ModifiedCressieRead(-.5, .2)
ch = [-Inf, -4.29317, -2.44994, -1.63336, -1.14658, -0.81439, -0.569175,
      -0.378594, -0.224971, -0.0977226, 0.00992562, 0.102539, 0.183387,
      0.26022, 0.337053, 0.413887, 0.49072, 0.567553, 0.644387, 0.72122,
      0.798053, 0.874887, 0.95172, 1.02855, 1.10539, 1.18222, 1.25905,
      1.33589, 1.41272, 1.48955]
[@test_approx_eq_eps ch[index] gradient(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

################################################################################
##
## Hessian - Modified Cressie Read
##
################################################################################

cr = ModifiedCressieRead(3., .2)
ch = [0., 0.010201, 0.040804, 0.091809, 0.163216, 0.255025, 0.367236, 0.499849,
      0.652864, 0.826281, 1.0201, 1.23432, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44,
      1.44, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44, 1.44]
[@test_approx_eq_eps ch[index] hessian(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

cr = ModifiedCressieRead(-1/3., .2)
ch = [+Inf, 21.2604, 8.4372, 4.91371, 3.3483, 2.48663, 1.95001, 1.58772,
      1.32878, 1.13566, 0.986821, 0.869056, 0.784197, 0.784197, 0.784197,
      0.784197, 0.784197, 0.784197, 0.784197, 0.784197, 0.784197, 0.784197,
      0.784197, 0.784197, 0.784197, 0.784197, 0.784197, 0.784197, 0.784197, 0.784197]
[@test_approx_eq_eps ch[index] hessian(cr, value) 1.e-4
									for (index, value) in enumerate(seq)];

#############################################################################
##
## Kullback Leibler
##
#############################################################################
seq =0:.1:3
kl = KullbackLeibler()
ch = [1.0, 0.669741, 0.478112, 0.338808, 0.233484, 0.153426, .0935046, .0503275,
	  .0214852, 0.00517554, 0., 0.0048412, 0.0187859, 0.0410735, 0.0710611, 0.108198,
	  0.152006, 0.202068, 0.258016, 0.319522, 0.386294, 0.458068, 0.534606, 0.615691,
	  0.701125, 0.790727, 0.88433, 0.98178, 1.08293, 1.18766, 1.29584]
[@test_approx_eq_eps ch[index] evaluate(kl, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

ch = [-Inf, -2.30259, -1.60944, -1.20397, -0.916291, -0.693147, -0.510826,
	  -0.356675, -0.223144, -0.105361, 0., 0.0953102, 0.182322, 0.262364,
	  0.336472, 0.405465, 0.470004, 0.530628, 0.587787, 0.641854, 0.693147,
	  0.741937, 0.788457, 0.832909, 0.875469, 0.916291, 0.955511, 0.993252,
	  1.02962, 1.06471, 1.09861]

[@test_approx_eq_eps ch[index] gradient(kl, value) 1.e-4
									for (index, value) in enumerate(seq)];

ch = [+Inf, 10., 5., 3.33333, 2.5, 2., 1.66667, 1.42857, 1.25,
	  1.11111, 1., 0.909091, 0.833333, 0.769231, 0.714286, 0.666667, 0.625,
	  0.588235, 0.555556, 0.526316, 0.5, 0.47619, 0.454545, 0.434783,
	  0.416667, 0.4, 0.384615, 0.37037, 0.357143, 0.344828, 0.333333]

[@test_approx_eq_eps ch[index] hessian(kl, value) 1.e-4
									for (index, value) in enumerate(seq)];

mkl = ModifiedKullbackLeibler(.2)
ch = [+Inf, 0.669741, 0.478112, 0.338808, 0.233484, 0.153426, 0.0935046,
	  0.0503275, 0.0214852, 0.00517554, 0., 0.0048412, 0.0187859, 0.0411847,
	  0.0719168, 0.110982, 0.158381, 0.214113, 0.278179, 0.350578, 0.43131,
	  0.520375, 0.617774, 0.723506, 0.837572, 0.959971, 1.0907, 1.22977, 1.37717,
	  1.5329, 1.69696]

[@test_approx_eq_eps ch[index] evaluate(mkl, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

seq =0:.11:3
ch = [-Inf,-2.20727,-1.51413,-1.10866,-0.820981,-0.597837,-0.415515,-0.261365,-0.127833,-0.0100503,0.0953102,0.190655,0.282322,0.373988,0.465655,0.557322,0.648988,0.740655,0.832322,0.923988,1.01565,1.10732,1.19899,1.29065,1.38232,1.47399,1.56565,1.65732]

[@test_approx_eq_eps ch[index] gradient(mkl, value) 1.e-4
									for (index, value) in enumerate(seq)];

ch = [+Inf, 9.09091,4.54545,3.0303,2.27273,1.81818,1.51515,1.2987,1.13636,1.0101,0.909091,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333,0.833333]

[@test_approx_eq_eps ch[index] hessian(mkl, value) 1.e-4
									for (index, value) in enumerate(seq)];

#############################################################################
##
## Reverse Kullback Leibler
##
#############################################################################

seq =0:.1:3
mrkl = ModifiedReverseKullbackLeibler(.2)

ch = [+Inf,1.40259,0.809438,0.503973,0.316291,0.193147,0.110826,0.0566749,0.0231436,0.00536052,0.,0.00468982,0.0176784,0.0378173,0.0649007,0.0989284,0.139901,0.187817,0.242678,0.304484,0.373234,0.448928,0.531567,0.621151,0.717678,0.821151,0.931567,1.04893,1.17323,1.30448,1.44268]

[@test_approx_eq_eps ch[index] evaluate(mrkl, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

ch = [+Inf,-8.09091,-3.54545,-2.0303,-1.27273,-0.818182,-0.515152,-0.298701,-0.136364,-0.010101,0.0909091,0.173611,0.25,0.326389,0.402778,0.479167,0.555556,0.631944,0.708333,0.784722,0.861111,0.9375,1.01389,1.09028,1.16667,1.24306,1.31944,1.39583]

ch = [+Inf,82.6446,20.6612,9.18274,5.16529,3.30579,2.29568,1.68663,1.29132,1.0203,0.826446,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444]

#############################################################################
##
## Fully Modified Kullback Leibler
##
#############################################################################

seq =-3:.1:3
mrkl = FullyModifiedReverseKullbackLeibler(.1, 1)

# \begin{array}{cc}
#  \{ &
# \begin{array}{cc}
#  50. (z-0.1)^2-9. (z-0.1)+1.40259 & z\leq 0.1 \\
#  z-\log (z)-1 & z>0.1\land z<2 \\
#  \frac{1}{8} (z-2)^2+\frac{z-2}{2}+1-\log (2) & z\geq 2 \\
# \end{array}
#  \\
# \end{array}

ch =    [509.80258509299404,478.402585092994,448.00258509299397,418.60258509299405,
         390.202585092994,362.80258509299404,336.402585092994,311.00258509299397,
         286.60258509299405,263.202585092994,240.802585092994,219.402585092994,
         199.002585092994,179.60258509299405,161.202585092994,143.80258509299404,
         127.40258509299403,112.00258509299402,97.60258509299405,84.20258509299404,
         71.80258509299404,60.402585092994016,50.00258509299401,40.60258509299402,
         32.20258509299401,24.802585092994043,18.402585092994038,13.002585092994035,
         8.602585092994033,5.202585092994035,2.8025850929940455,1.4025850929940447,
         0.8094379124340996,0.5039728043259353,0.3162907318741546,0.1931471805599453,
         0.11082562376599059,0.05667494393873229,0.023143551314209698,0.005360515657826262,
         0.,0.004689820195675182,0.01767844320604539,0.03763573553250893,0.06352776337878718,
         0.09453489189183562,0.12999637075426462,0.1693717489378297,0.21221333509788132,
         0.25814611382760544,0.3068528194400547,0.35810281944005495,0.4118528194400547,
         0.46810281944005505,0.5268528194400549,0.5881028194400547,0.6518528194400551,
         0.7181028194400548,0.7868528194400551,0.858102819440055,0.9318528194400547]

[@test_approx_eq_eps ch[index] evaluate(mrkl, [value]) 1.e-4
									for (index, value) in enumerate(seq)];

ch = [-319., -309., -299., -289., -279., -269., -259., -249., -239., -229.,
      -219., -209., -199., -189., -179., -169., -159., -149., -139., -129.,
      -119., -109., -99., -89., -79., -69., -59., -49., -39., -29., -19.,
      -9., -4., -2.33333, -1.5, -1., -0.666667, -0.428571, -0.25,
      -0.111111, 0., 0.0909091, 0.166667, 0.230769, 0.285714, 0.333333,
      0.375, 0.411765, 0.444444, 0.473684, 0.5, 0.525, 0.55, 0.575, 0.6,
      0.625, 0.65, 0.675, 0.7, 0.725, 0.75]

[@test_approx_eq_eps ch[index] gradient(mrkl, value) 1.e-4
									for (index, value) in enumerate(seq)];

#############################################################################
##
## Both a and b
##
#############################################################################

a = .1
b = .2

@test evaluate(KullbackLeibler(), [a], [b]) == 0.030685281944005494
@test evaluate(ReverseKullbackLeibler(), [a], [b]) == 0.038629436111989046
@test evaluate(CressieRead(3), a, b) == 0.017708333333333333
@test evaluate(CressieRead(-3), a, b) == 0.06666666666666667
@test evaluate(Divergences.HD(), a, b) == 0.03431457505076194
@test evaluate(Divergences.ChiSquared(), a, b) == (a-b)^2/(2*b)

@test gradient(KullbackLeibler(), a, b) == log(a/b)
@test gradient(ReverseKullbackLeibler(), a, b) == 1.0-b/a
@test gradient(CressieRead(3), a, b) == ((a/b)^(3)-1)/3
@test gradient(CressieRead(-3), a, b) == ((a/b)^(-3)-1)/(-3)
@test gradient(Divergences.HD(), a, b) == ((a/b)^(-1/2)-1)/(-1/2)
@test gradient(Divergences.ChiSquared(), a, b) == a/b-1.0

@test hessian(Divergences.ChiSquared(), [a], [b]) == [1/b]
@test hessian(KullbackLeibler(), [a], [b]) == [1/a]
@test hessian(ReverseKullbackLeibler(), [a], [b]) == [b/a^2]

@test hessian(CressieRead(3), [a], [b]) ==    [( (a/b)^3 )/a]
@test hessian(CressieRead(1), [a], [b]) ==    [( (a/b)   )/a]
@test hessian(CressieRead(-1/2), [a], [b]) == [( (a/b)^(-.5) )/a]

@test_approx_eq evaluate(CressieRead(1), [a], [b])[1] evaluate(ChiSquared(), [a], [b])[1]
@test hessian(CressieRead(1), [a], [b])[1] == hessian(ChiSquared(), [a], [b])[1]
@test gradient(CressieRead(1), [a], [b])[1] == gradient(ChiSquared(), [a], [b])[1]

## Additional tests

@test CressieRead(1) === CressieRead(1.0)
@test CressieRead(-1/2) === HD()
@test_throws(AssertionError, CressieRead(-1))
@test_throws(AssertionError, ModifiedCressieRead(-1, -1))
@test_throws(AssertionError, ModifiedKullbackLeibler(-1))
@test_throws(AssertionError, ModifiedReverseKullbackLeibler(-1))

@test_throws(AssertionError, FullyModifiedCressieRead(1, -1, 2))
@test_throws(AssertionError, FullyModifiedKullbackLeibler(-1, 2))
@test_throws(AssertionError, FullyModifiedReverseKullbackLeibler(-1, 2))

@test_throws(AssertionError, FullyModifiedCressieRead(1, .5, -2))
@test_throws(AssertionError, FullyModifiedKullbackLeibler(.5, -2))
@test_throws(AssertionError, FullyModifiedReverseKullbackLeibler(.5, -2))

@test_throws(DimensionMismatch, evaluate(KL(), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(RKL(), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(CR(1), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(ChiSquared(), rand(10), rand(11)))

@test_throws(DimensionMismatch, evaluate(MKL(1), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(MRKL(1), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(MCR(1, 1), rand(10), rand(11)))

@test_throws(DimensionMismatch, evaluate(FMKL(.5, .5), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(FMRKL(.5, .5), rand(10), rand(11)))
@test_throws(DimensionMismatch, evaluate(FMCR(1, .5, .5), rand(10), rand(11)))

## Test normalization

@test evaluate(KL(), [1.], [1.]) == 0.0
@test evaluate(RKL(), [1.], [1.]) == 0.0
@test evaluate(CR(2), [1.], [1.]) == 0.0
@test evaluate(CR(-3/2), [1.], [1.]) == 0.0

@test gradient(KL(), [1.], [1.]) == [0.0]
@test gradient(RKL(), [1.], [1.]) == [0.0]
@test gradient(CR(2), [1.], [1.]) == [0.0]
@test gradient(CR(-3/2), [1.], [1.]) == [0.0]

@test hessian(KL(), [1.], [1.]) == [1.0]
@test hessian(RKL(), [1.], [1.]) == [1.0]
@test hessian(CR(2), [1.], [1.]) == [1.0]
@test hessian(CR(-3/2), [1.], [1.]) == [1.0]

## Modified CR with both a, b

@test evaluate(MCR(2, .2), [.2], [.1]) == 0.05813333333333335
@test evaluate(MCR(2, .2), [.3], [.1]) == 0.23613333333333322
@test evaluate(MCR(2, .2), [.4], [.1]) == 0.5341333333333332

@test evaluate(MCR(-2, .2), [.2], [.1]) == 0.032407407407407406
@test evaluate(MCR(-2, .2), [.3], [.1]) == 0.12291666666666663
@test evaluate(MCR(-2, .2), [.4], [.1]) == 0.2712962962962963

@test evaluate(FMCR(2, .1, .2), [.2], [.1]) == 0.05813333333333335
@test evaluate(FMCR(2, .1, .2), [.3], [.1]) == 0.23613333333333322
@test evaluate(FMCR(2, .1, .2), [.4], [.1]) == 0.5341333333333332

@test evaluate(FMCR(-2, .1, .2), [.2], [.1]) == 0.032407407407407406
@test evaluate(FMCR(-2, .1, .2), [.3], [.1]) == 0.12291666666666663
@test evaluate(FMCR(-2, .1, .2), [.4], [.1]) == 0.2712962962962963

@test evaluate(FMCR(-2, .1, .2), [-.2], [.1]) == 231.29999999999995
@test evaluate(FMCR(-2, .1, .2), [-.3], [.1]) == 496.24999999999983
@test evaluate(FMCR(2, .1, .2), [-.2], [.1]) == 0.15435000000000001

@test gradient(FMCR(2, .1, .2), [.25], [.2]) == gradient(MCR(2, .2), [.25], [.2])
@test gradient(FMCR(-2, .1, .2), [.25], [.2]) == gradient(MCR(-2, .2), [.25], [.2])

@test_approx_eq gradient(FMCR(2, .1, .2), [-2.], [.2]) -1.505
@test_approx_eq gradient(FMCR(2, .1, .2), [.21], [.2]) 0.05125
@test_approx_eq gradient(FMCR(2, .1, .2), [2.], [.2]) 10.78

@test_approx_eq hessian(FMCR(2, .1, .2), [-2.], [.2]) 0.5
@test_approx_eq hessian(FMCR(2, .1, .2), [.21], [.2]) 5.25
@test_approx_eq hessian(FMCR(2, .1, .2), [2.], [.2]) 6

@test evaluate(MCR(-2, .1), [-.2], [.1]) == Inf
@test evaluate(KL(), [-.2], [.1]) == Inf
@test evaluate(RKL(), [-.2], [.1]) == Inf
@test evaluate(CR(1), [-.2], [.1]) == Inf
@test evaluate(ChiSquared(), [-.2], [.1]) == 0.4500000000000001

## Chi-squared
@test evaluate(ChiSquared(), .2, .1) == evaluate(ChiSquared(), [.2], [.1])
@test_approx_eq evaluate(ChiSquared(), .2, .1) evaluate(CR(1), [.2], [.1])
@test_approx_eq gradient(ChiSquared(), .2, .1) gradient(CR(1), [.2], [.1])
@test_approx_eq hessian(ChiSquared(), .2, .1) hessian(CR(1), [.2], [.1])
@test_approx_eq gradient(ChiSquared(), .2) gradient(CR(1), [.2])
@test_approx_eq hessian(ChiSquared(), .2) hessian(CR(1), [.2])
@test_approx_eq evaluate(ChiSquared(), [.2]) evaluate(CR(1), [.2])
@test_approx_eq gradient(ChiSquared(), [.2]) gradient(CR(1), [.2])
@test_approx_eq hessian(ChiSquared(), [.2]) hessian(CR(1), [.2])

## KL
@test evaluate(KL(), [.2], [1.0]) == evaluate(KL(), [.2], [1.0])
@test gradient(KL(), [.2], [1.0]) == gradient(KL(), [.2], [1.0])
@test hessian(KL(), [.2], [1.0]) == hessian(KL(), [.2], [1.0])

@test gradient(KL(), [0.], [1.0]) == gradient(KL(), [0.], [1.0])
@test hessian(KL(), [0.], [1.0]) == hessian(KL(), [0.], [1.0])

## RKL
@test evaluate(RKL(), [.2], [1.0]) == evaluate(RKL(), [.2], [1.0])
@test gradient(RKL(), [.2], [1.0]) == gradient(RKL(), [.2], [1.0])
@test hessian(RKL(), [.2], [1.0]) == hessian(RKL(), [.2], [1.0])

@test evaluate(RKL(), [0.], [1.0]) == evaluate(RKL(), [0.], [1.0])
@test gradient(RKL(), [0.], [1.0]) == gradient(RKL(), [0.], [1.0])
@test hessian(RKL(), [0.], [1.0]) == hessian(RKL(), [0.], [1.0])


text:
# batched routines
@inline indexpos(iAs, i)::Int = findfirst(==(i), iAs)

# can be used in either static or dynamic invoke
function analyse_batched(iAs, sA, iBs, sB, iOuts)
    iABs = iAs ∩ iBs
    pres   = iABs ∩ iOuts
    broad  = setdiff((iAs ∩ iOuts) ∪ (iBs ∩ iOuts), pres)
    summed = setdiff(iABs, pres)

    iAps, iAbs, iAss = pres ∩ iAs, broad ∩ iAs, summed ∩ iAs
    iBps, iBbs, iBss = pres ∩ iBs, broad ∩ iBs, summed ∩ iBs

    pA   = indexpos.(Ref(iAs), vcat(iAbs, iAss, iAps))
    pB   = indexpos.(Ref(iBs), vcat(iBss, iBbs, iBps))
    iABs = vcat(iAbs, iBbs, iAps)
    pOut = indexpos.(Ref(iABs), iOuts)

    sAbs = getindex.(Ref(sA), indexpos.(Ref(iAs), iAbs))
    sAb = reduce(*, sAbs, init=1)
    sAs = mapreduce(i->sA[indexpos(iAs,i)], *, iAss, init=1)
    sAps = getindex.(Ref(sA), indexpos.(Ref(iAs),iAps))
    sAp = reduce(*, sAps, init=1)

    sBbs = getindex.(Ref(sB), indexpos.(Ref(iBs),iBbs))
    sBb = reduce(*, sBbs, init=1)
    sBs = mapreduce(i->sB[indexpos(iBs, i)], *, iBss, init=1)
    sBp = mapreduce(i->sB[indexpos(iBs, i)], *, iBps, init=1)
    sAB = vcat(sAbs, sBbs, sAps)
    return pA, sAb, sAs, sAp, pB, sBs, sBb, sBp, sAB, pOut
end

# batched, dynamic version
function batched_contract(iAs, A::AbstractArray, iBs, B::AbstractArray, iOuts::NTuple{NO,T}) where {NO,T}
    pA, sAb, sAs, sAp, pB, sBs, sBb, sBp, sAB, pOut = analyse_batched(iAs, size(A), iBs, size(B), iOuts)

    A, B = align_eltypes(A, B)
    Apr = reshape(conditioned_permutedims(A, pA), sAb, sAs, sAp)
    Bpr = reshape(conditioned_permutedims(B, pB), sBs, sBb, sBp)
    AB = _batched_gemm('N','N', Apr, Bpr)
    AB = conditioned_permutedims(reshape(AB, sAB...), [pOut...])
end

# reload this function for GPU support!
function _batched_gemm(C1::Char, C2::Char, A::StridedArray{T,3}, B::StridedArray{T2,3}) where {T<:Number, T2<:Number}
    batched_gemm(C1, C2, A, B)
end