text:
module Impute

using AutoHashEquals
using IterTools
using Statistics
using StatsBase
using Tables: Tables, materializer, istable

import Base.Iterators: drop

export impute, impute!, chain, chain!, drop, drop!, interp, interp!, ImputeError

function __init__()
    sym = join(["chain", "chain!", "drop", "drop!", "interp", "interp!"], ", ", " and ")

    @warn(
        """
        The following symbols will not be exported in future releases: $sym.
        Please qualify your calls with `Impute.<method>(...)` or explicitly import the symbol.
        """
    )

    @warn(
        """
        The default limit for all impute functions will be 1.0 going forward.
        If you depend on a specific threshold please pass in an appropriate `AbstractContext`.
        """
    )
end

"""
    ImputeError{T} <: Exception

Is thrown by `impute` methods when the limit of imputable values has been exceeded.

# Fields
* msg::T - the message to print.
"""
struct ImputeError{T} <: Exception
    msg::T
end

Base.showerror(io::IO, err::ImputeError) = println(io, "ImputeError: $(err.msg)")

include("context.jl")
include("imputors.jl")

const global imputation_methods = (
    drop = DropObs,
    dropobs = DropObs,
    dropvars = DropVars,
    interp = Interpolate,
    fill = Fill,
    locf = LOCF,
    nocb = NOCB,
)

include("deprecated.jl")

for (f, v) in pairs(imputation_methods)
    typename = nameof(v)
    f! = Symbol(f, :!)

    @eval begin
        $f(data; kwargs...) = impute(data, $typename(; _extract_context_kwargs(kwargs...)...))
        $f!(data; kwargs...) = impute!(data, $typename(; _extract_context_kwargs(kwargs...)...))
        $f(; kwargs...) = data -> impute(data, $typename(; _extract_context_kwargs(kwargs...)...))
        $f!(; kwargs...) = data -> impute!(data, $typename(; _extract_context_kwargs(kwargs...)...))
    end
end

@doc """
    Impute.dropobs(data; vardim=2, context=Context())

Removes missing observations from the `AbstractArray` or `Tables.table` provided.
See [DropObs](@ref) for details.

# Example
```
julia> using DataFrames; using Impute: Impute, Context

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.dropobs(df; vardim=1, context=Context(; limit=1.0))
3×2 DataFrames.DataFrame
│ Row │ a       │ b       │
│     │ Float64 │ Float64 │
├─────┼─────────┼─────────┤
│ 1   │ 1.0     │ 1.1     │
│ 2   │ 2.0     │ 2.2     │
│ 3   │ 5.0     │ 5.5     │
```
""" dropobs

@doc """
    Impute.dropvars(data; vardim=2, context=Context())

Finds variables with too many missing values in a `AbstractMatrix` or `Tables.table` and
removes them from the input data. See [DropVars](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute, Context

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.dropvars(df; vardim=1, context=Context(; limit=0.2))
5×1 DataFrames.DataFrame
│ Row │ b        │
│     │ Float64⍰ │
├─────┼──────────┤
│ 1   │ 1.1      │
│ 2   │ 2.2      │
│ 3   │ 3.3      │
│ 4   │ missing  │
│ 5   │ 5.5      │
```
""" dropvars

@doc """
    Impute.interp(data; vardim=2, context=Context())

Performs linear interpolation between the nearest values in an vector.
See [Interpolate](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute, Context

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.interp(df; vardim=1, context=Context(; limit=1.0))
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ 3.0      │ 3.3      │
│ 4   │ 4.0      │ 4.4      │
│ 5   │ 5.0      │ 5.5      │
```
""" interp

@doc """
    Impute.fill(data; value=mean, vardim=2, context=Context())

Fills in the missing data with a specific value. See [Fill](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute, Context

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.fill(df; value=-1.0, vardim=1, context=Context(; limit=1.0))
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ -1.0     │ 3.3      │
│ 4   │ -1.0     │ -1.0     │
│ 5   │ 5.0      │ 5.5      │
```
""" fill

@doc """
    Impute.locf(data; vardim=2, context=Context())

Iterates forwards through the `data` and fills missing data with the last existing
observation. See [LOCF](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute, Context

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.locf(df; vardim=1, context=Context(; limit=1.0))
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ 2.0      │ 3.3      │
│ 4   │ 2.0      │ 3.3      │
│ 5   │ 5.0      │ 5.5      │
```
""" locf

@doc """
    Impute.nocb(data; vardim=2, context=Context())

Iterates backwards through the `data` and fills missing data with the next existing
observation. See [LOCF](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute, Context

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.nocb(df; vardim=1, context=Context(; limit=1.0))
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64⍰ │ Float64⍰ │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ 5.0      │ 3.3      │
│ 4   │ 5.0      │ 5.5      │
│ 5   │ 5.0      │ 5.5      │
```
""" nocb

end  # module


text:
# module for drawing moving robots
# add random noise against robot moving

module AnimeMoveRandomNoise
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function animate_per_time(time_interval, world, noise_robot,
                            no_noise_robot)
    draw(world)
    annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
    draw!(noise_robot)
    draw!(no_noise_robot)
  end

  function main(time_interval=0.1; is_test=false)
    w = World(-5.0, 5.0, -5.0, 5.0)

    init_pose = [0.0, 0.0, 0.0]
    circling = Agent(0.2, 10.0/180*pi)

    nr = DifferentialWheeledRobot(init_pose, 0.2, "gray",
                                  circling, time_interval,
                                  noise_per_meter=5,
                                  noise_std=pi/30)

    nnr = DifferentialWheeledRobot(init_pose, 0.2, "red",
                                   circling, time_interval)

    if is_test == false
      anime = @animate for t in 0:time_interval:30
        animate_per_time(t, w, nr, nnr)
      end
      path = "src/model/uncertainty/movement/random_noise/anime_move_random_noise.gif"
      gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], path))
    else
      for t in 0:time_interval:10
        animate_per_time(t, w, nr, nnr)
      end
    end
  end
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
##########
##########
# Set equations to track biomass in each population
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-06-25
##########
##########

# set up =======================================================================

using DifferentialEquations
include("./src/eqn_component_functions.jl")
include("./src/species_list_creation.jl")
include("./src/network_metrics.jl")

"""
    check_extinctions!(x)

Check if any species have low enough abundance to warrent extinction - if so,
remove them from the array of abundances right away.
"""
function check_extinctions!(x::Array)
    for i in 1:length(x)
        if x[i] < 1.0e-30
            x[i] = 0
        end
    end
    return x
end

"""
biomass!(dx, x, p, t)

Calling biomass functions, create the actual equation for each species, to be
passed to the simulation.
"""
function biomass!(dx, x, p, t)
    for i in 1:size(species_list, 1)
        dx[i] = growth_rate(i, species_list, x) +
                biomass_gain(species_list, network, x, 2, i) -
                biomass_loss(species_list, network, x, 2, i) -
                mortality_change(i, x)
    end
end

function starting_eq_sims(
    species_list::Array,
    network::Array
    )

    #  allocate memory
    dx = zeros(50)

    # save all initial parameters (maybe a seaprate object)

    # sims

function perturbation_simulation(
    species_list:::Array,
    network::Array,
    type_of_perturbation::String,
    x::Array
    )

# pre-allocate memory for the objects we want to save from the time series
# save initial parameters
# sim from the starting biomass out to 5000 - save that point
# ^ maybe every 10 points stop it and check for extinction in any of them
# ^ also, add a component to the species_list to keep track of alive (1) or
#    extinct (0)
# check the shape of the web at that point & save
# sim 10 at a time, check for extinctions, check for perturbation ,
#    get the parameters again

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
using Documenter, Simpsons

makedocs(sitename = "Simpsons Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
"""
    FGMSolver
"""
mutable struct FGMSolver
    step::FGM.AbstractStep
    cond::Vector{FGM.AbstractStopCondition}

    function FGMSolver( step::FGM.AbstractStep = FGM.ConstantStep{T}(), sc::Vector{FGM.AbstractStopCondition} = [FGM.Best(1e-4)] )
        solver = new( step, sc )

        return solver
    end
end

function createsolver!( ::FGMSolver )

end

function updatesolver!( ::FGMSolver )

end

function solve!( ::FGMSolver  )

end


text:
const ASSET_REGISTRY = Dict{String, String}()

function unique_file_key(path::String)
    return bytes2hex(sha1(abspath(path))) * "-" * basename(path)
end

function register_local_file(file::String)
    target = normpath(abspath(expanduser(file)))
    key = "/assetserver/" * unique_file_key(target)
    get!(()-> target, ASSET_REGISTRY, key)
    return key
end

function is_key_registered(key::String)
    return haskey(ASSET_REGISTRY, key)
end

function assetserver_to_localfile(key::String)
    path = get(ASSET_REGISTRY, key, nothing)
    if path === nothing
        error("Key does not map to a local path (is not registered via `register_local_file`): $(key)")
    end
    return path
end

include("mimetypes.jl")

function file_server(context)
    path = context.request.target
    if is_key_registered(path)
        filepath = assetserver_to_localfile(path)
        if isfile(filepath)
            header = ["Access-Control-Allow-Origin" => "*",
                      "Content-Type" => file_mimetype(filepath)]
            return HTTP.Response(200, header, body = read(filepath))
        end
    end
    return HTTP.Response(404)
end

"""
    dependency_path(paths...)

Path to serve downloaded dependencies
"""
dependency_path(paths...) = joinpath(@__DIR__, "..", "js_dependencies", paths...)

mediatype(asset::Asset) = asset.media_type

function Base.show(io::IO, asset::Asset)
    print(io, url(asset))
end

function Asset(online_path::String, onload::Union{Nothing, JSCode}=nothing; check_isfile=false)
    local_path = ""; real_online_path = ""
    if is_online(online_path)
        local_path = ""
        real_online_path = online_path
    else
        local_path = normpath(abspath(expanduser(online_path)))
        if check_isfile && !isfile(local_path)
            error("File $(local_path) does not exist!")
        end
    end
    return Asset(Symbol(getextension(online_path)), real_online_path, local_path, onload)
end

"""
    getextension(path)
Get the file extension of the path.
The extension is defined to be the bit after the last dot, excluding any query
string.
# Examples
```julia-repl
julia> JSServe.getextension("foo.bar.js")
"js"
julia> JSServe.getextension("https://my-cdn.net/foo.bar.css?version=1")
"css"
```
Taken from WebIO.jl
"""
getextension(path) = lowercase(last(split(first(split(path, "?")), ".")))

"""
    is_online(path)

Determine whether or not the specified path is a local filesystem path (and not
a remote resource that is hosted on, for example, a CDN).
"""
is_online(path) = any(startswith.(path, ("//", "https://", "http://", "ftp://")))

function Dependency(name::Symbol, urls::AbstractVector)
    return Dependency(name, Asset.(urls), Dict{Symbol, JSCode}())
end

"""
    construct_arguments(args, keyword_arguments)
Constructs the arguments for a JS call.
Can only use either keyword arguments or positional arguments.
"""
function construct_arguments(args, keyword_arguments)
    if isempty(keyword_arguments)
        return args
    elseif isempty(args)
        # tojs isn't recursive bug:
        return keyword_arguments
    else
        # TODO: I'm not actually sure about this :D
        error("""
        Javascript only supports keyword arguments OR arguments.
        Found posititional arguments and keyword arguments
        """)
    end
end

"""
Implement to call functions in a Dependency:
const Three = Dependency(...)
Three.camera(session, args...; kw...)
"""
function Base.getproperty(dependency::Dependency, func_name::Symbol)
    func_name in (:name, :assets, :function) && return getfield(dependency, func_name)
    func_name_js = JSString(string(func_name))
    return function (session::Session, args...; kw...)
        args = construct_arguments(args, kw)
        evaljs(session, js"""
            $(dependency).$(func_name_js)(...$(args))
        """)
    end
end

# With this, one can just put a dependency anywhere in the dom to get loaded
function jsrender(session::Session, x::Dependency)
    push!(session, x)
    return nothing
end

function jsrender(session::Session, asset::Asset)
    register_resource!(session, asset)
    return nothing
end

function include_asset(asset::Asset, serializer::UrlSerializer=UrlSerializer())
    function to_url(mime)
        # In case we don't rely on a running server,
        # we base64 encode the content in a data url
        # This isn't needed, if the asset is already hosted online (!isempty(asset.online_path))
        if serializer.inline_all && isempty(asset.online_path)
            data_str = base64encode(read(asset.local_path))
            return "data:$(mime);base64," * data_str
        else
            return url(asset, serializer)
        end
    end
    if mediatype(asset) == :js
        return DOM.script(src=to_url("text/javascript"))
    elseif mediatype(asset) == :css
        return DOM.link(href=to_url("text/css"), rel="stylesheet", type="text/css")
    else
        error("Unrecognized asset media type: $(mediatype(asset))")
    end
end

function include_asset(dep::Dependency, serializer::UrlSerializer=UrlSerializer())
    if length(dep.assets) == 1
        include_asset(dep.assets[1], serializer)
    else
        return include_asset(dep.assets, serializer)
    end
end

function include_asset(assets::Union{Vector{Asset}, Set{Asset}}, serializer::UrlSerializer=UrlSerializer())
    return DOM.div(include_asset.(assets, (serializer,))...)
end

function url(str::String, serializer::UrlSerializer=UrlSerializer())
    return url(Asset(str), serializer)
end

function unique_name_in_folder(folder, dream_name)
    taken_names = readdir(folder)
    try_name = dream_name
    unique_id = 1
    # Large number of tries - I would make this a while true
    # but I think its better to have a fixed termination!
    for unique_id in 1:1000000
        if !(try_name in taken_names)
            return try_name
        end
        try_name = string(unique_id, "_", dream_name)
    end
    error("No unique names found after 1000000 tries for name $(dream_name) in $(folder)")
end

function url(asset::Asset, serializer::UrlSerializer=UrlSerializer())
    if isempty(asset.online_path)
        path = asset.local_path
        relative_path = if serializer.assetserver
            # we use assetserver, so we register the local file with the server
            register_local_file(path)
        else
            # we don't use assetserver, so we copy the asset to asset_folder
            # for someone else to serve them!
            if serializer.asset_folder === nothing
                error("Not using assetserver requires to set `asset_folder` to a valid local folder")
            end
            if !isdir(serializer.asset_folder)
                error("`asset_folder` doesn't exist: $(serializer.asset_folder)")
            end
            # make sure we don't get name clashes in asset_folder
            path_base = dirname(path)
            file_name = basename(path)
            unique_file_name = unique_name_in_folder(serializer.asset_folder, file_name)
            unique_path = joinpath(serializer.asset_folder, unique_file_name)
            cp(path, unique_path, force=true)
            unique_file_name
        end
        if serializer.absolute
            # we serve from an absolute URL, so we append the content_delivery_url!
            return serializer.content_delivery_url * relative_path
        else
            return relative_path
        end
    else
        # We dont touch online urls for now
        # we may to decide to download them and add them to the asset_folder
        # if desired
        return asset.online_path
    end
end

const MsgPackLib = Asset(dependency_path("msgpack.min.js"))
const PakoLib = Asset(dependency_path("pako_inflate.min.js"))
const JSServeLib = Dependency(:JSServe, [dependency_path("JSServe.js")])
const Base64Lib = Dependency(:Base64, [dependency_path("Base64.js")])

const MarkdownCSS = Asset(dependency_path("markdown.css"))
const TailwindCSS = Asset(dependency_path("tailwind.min.css"))
const Styling = Asset(dependency_path("styled.css"))


text:
#    Problem : GROUP A
#    *********
#    The chained Rosenbrock function as given by Fletcher.
#
#    Original SIF Source:
#    The second problem given by
#    R. Fletcher,
#    "An optimal positive definite update for sparse Hessian matrices"
#    Numerical Analysis report NA/145, University of Dundee, 1992.
#
#    FLETCHR.SIF classification OUR2-AN-V-0
#
#    The Number of variables is variable.
#
# Daniel Henderson, 08/2021

f = (x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)-1
        γ = x[i+1] - x[i] + 1 - x[i]^2
        fx += 100*γ^2
    end
    return fx
end

g! = (g, x) -> begin
    for i in firstindex(x):lastindex(x)-1
        g[i] += 200*(1 + 2x[i])(-1+x[i]+x[i]^2-x[i+1])
        g[i+1] += 200*(1 - x[i] - x[i]^2 + x[i+1])
    end
    return g
end

fg! = (g, x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)-1
        γ = x[i+1] - x[i] + 1 - x[i]^2
        fx += 100*γ^2
        g[i] += 200*(-1 - 2x[i])*(1 - x[i] - x[i]^2 + x[1+i])
        g[i+1] += 200*(1 - x[i] - x[i]^2 + x[1+i])
    end
    return fx, g
end

init = (n::Int=1000) -> begin
	n < 2 && @warn("FLETCHCR: number of variables must be ≥ 2")
	n = max(n, 2)

    return n, zeros(n)
end

@warn "TODO: FLETCHCR, determine starting iterate - see whiteboard"
# TestSet["FLETCHCR"] = UncProgram("FLETCHCR", f, g!, fg!, init)


text:
module DynamicPPLBenchmarks

using DynamicPPL
using BenchmarkTools

using Weave: Weave
using Markdown: Markdown

using LibGit2: LibGit2
using Pkg: Pkg

export weave_benchmarks

function time_model_def(model_def, args...)
    return @time model_def(args...)
end

function benchmark_untyped_varinfo!(suite, m)
    vi = VarInfo()
    # Populate.
    m(vi)
    # Evaluate.
    suite["evaluation_untyped"] = @benchmarkable $m($vi, $(DefaultContext()))
    return suite
end

function benchmark_typed_varinfo!(suite, m)
    # Populate.
    vi = VarInfo(m)
    # Evaluate.
    suite["evaluation_typed"] = @benchmarkable $m($vi, $(DefaultContext()))
    return suite
end

function typed_code(m, vi=VarInfo(m))
    rng = DynamicPPL.Random.MersenneTwister(42)
    spl = DynamicPPL.SampleFromPrior()
    ctx = DynamicPPL.SamplingContext(rng, spl, DynamicPPL.DefaultContext())

    results = code_typed(m.f, Base.typesof(m, vi, ctx, m.args...))
    return first(results)
end

"""
    make_suite(model)

Create default benchmark suite for `model`.
"""
function make_suite(model)
    suite = BenchmarkGroup()
    benchmark_untyped_varinfo!(suite, model)
    benchmark_typed_varinfo!(suite, model)

    return suite
end

"""
    weave_child(indoc; mod, args, kwargs...)

Weave `indoc` with scope of `mod` into markdown.

Useful for weaving within weaving, e.g.
```julia
weave_child(child_jmd_path, mod = @__MODULE__, args = WEAVE_ARGS)
```
together with `results="markup"` and `echo=false` will simply insert
the weaved version of `indoc`.

# Notes
- Currently only supports `doctype == "github"`. Other outputs are "supported"
  in the sense that it works but you might lose niceties such as syntax highlighting.
"""
function weave_child(indoc; mod, args, kwargs...)
    # FIXME: Make this work for other output formats than just `github`.
    doc = Weave.WeaveDoc(indoc, nothing)
    doc = Weave.run_doc(doc; doctype="github", mod=mod, args=args, kwargs...)
    rendered = Weave.render_doc(doc)
    return display(Markdown.parse(rendered))
end

"""
    pkgversion(m::Module)

Return version of module `m` as listed in its Project.toml.
"""
function pkgversion(m::Module)
    projecttoml_path = joinpath(dirname(pathof(m)), "..", "Project.toml")
    return Pkg.TOML.parsefile(projecttoml_path)["version"]
end

"""
    default_name(; include_commit_id=false)

Construct a name from either repo information or package version
of `DynamicPPL`.

If the path of `DynamicPPL` is a git-repo, return name of current branch,
joined with the commit id if `include_commit_id` is `true`.

If path of `DynamicPPL` is _not_ a git-repo, it is assumed to be a release,
resulting in a name of the form `release-VERSION`.
"""
function default_name(; include_commit_id=false)
    dppl_path = abspath(joinpath(dirname(pathof(DynamicPPL)), ".."))

    # Extract branch name and commit id
    local name
    try
        githead = LibGit2.head(LibGit2.GitRepo(dppl_path))
        branchname = LibGit2.shortname(githead)

        name = replace(branchname, "/" => "_")
        if include_commit_id
            gitcommit = LibGit2.peel(LibGit2.GitCommit, githead)
            commitid = string(LibGit2.GitHash(gitcommit))
            name *= "-$(commitid)"
        end
    catch e
        if e isa LibGit2.GitError
            @info "No git repo found for $(dppl_path); extracting name from package version."
            name = "release-$(pkgversion(DynamicPPL))"
        else
            rethrow(e)
        end
    end

    return name
end

"""
    weave_benchmarks(input="benchmarks.jmd"; kwargs...)

Weave benchmarks present in `benchmarks.jmd` into a single file.

# Keyword arguments
- `benchmarkbody`: JMD-file to be rendered for each model.
- `include_commit_id=false`: specify whether to include commit-id in the default name.
- `name`: the name of directory in `results/` to use as output directory.
- `name_old=nothing`: if specified, comparisons of current run vs. the run pinted to
  by `name_old` will be included in the generated document.
- `include_typed_code=false`: if `true`, output of `code_typed` for the evaluator
  of the model will be included in the weaved document.
- Rest of the passed `kwargs` will be passed on to `Weave.weave`.
"""
function weave_benchmarks(
    input=joinpath(dirname(pathof(DynamicPPLBenchmarks)), "..", "benchmarks.jmd");
    benchmarkbody=joinpath(
        dirname(pathof(DynamicPPLBenchmarks)), "..", "benchmark_body.jmd"
    ),
    include_commit_id=false,
    name=default_name(; include_commit_id=include_commit_id),
    name_old=nothing,
    include_typed_code=false,
    doctype="github",
    outpath="results/$(name)/",
    kwargs...,
)
    args = Dict(
        :benchmarkbody => benchmarkbody,
        :name => name,
        :include_typed_code => include_typed_code,
    )
    if !isnothing(name_old)
        args[:name_old] = name_old
    end
    @info "Storing output in $(outpath)"
    mkpath(outpath)
    return Weave.weave(input, doctype; out_path=outpath, args=args, kwargs...)
end

end # module

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
## FALL-BACK METHODS FOR COMPOSITE MODELS (EXPORTED LEARNING NETWORKS)

# *Note.* Be sure to read Note 4 in src/operations.jl to see see how
# fallbacks are provided for operations acting on Composite models.

caches_data_by_default(::Type{<:Composite}) = true

fitted_params(::Union{Composite,Surrogate},
              fitresult::NamedTuple) =
                  fitted_params(glb(values(fitresult)...))

# legacy code (fitresult a Node):
fitted_params(::Union{Composite,Surrogate},
              fitresult::Node) =
                  fitted_params(fitresult)

function update(model::M,
                verbosity::Integer,
                fitresult::NamedTuple,
                cache,
                args...) where M <: Composite

    # This method falls back to `fit` to force rebuilding the
    # underlying learning network if, since the last fit:
    #
    # (i) Any field associated with a model in the learning network
    #     has been replaced with a new model instance (and not merely
    #     mutated), OR

    # (ii) Any OTHER field has changed it's value (in the sense
    # of `==`).

    # Otherwise, a "smart" fit is carried out by calling `fit!` on a
    # greatest lower bound node for nodes in the signature of the
    # underlying learning network machine. For this it is necessary to
    # temporarily "de-anonymize" the source nodes.

    network_model_fields = cache.network_model_fields
    old_model = cache.old_model

    glb_node = glb(values(fitresult)...) # greatest lower bound

    if fallback(model, old_model, network_model_fields, glb_node)
        return fit(model, verbosity, args...)
    end

    # return data to source nodes for fitting:
    sources, data = cache.sources, cache.data
    for k in eachindex(sources)
        rebind!(sources[k], data[k])
    end

    fit!(glb_node; verbosity=verbosity)

    # anonymize data again:
    for s in sources
        rebind!(s, nothing)
    end

    # record current model state:
    cache = (sources=cache.sources,
             data=cache.data,
             network_model_fields=cache.network_model_fields,
             old_model = deepcopy(model))

    return fitresult, cache, report(glb_node)

end

# helper for preceding method (where logic is explained):
function fallback(model::M, old_model, network_model_fields, glb_node) where M
    network_model_ids = objectid.(models(glb_node))
    for name in fieldnames(M)
        value = getproperty(model, name)
        if name in network_model_fields
            objectid(value) in network_model_ids || return true
        else
            old_value = getproperty(old_model, name)
            value == old_value || return true
        end
    end
    return false
end

# legacy version of above (fitresult a Node):
function update(model::Composite,
                verbosity::Integer,
                yhat::Node,
                cache,
                args...)

    # If any `model` field has been replaced (and not just mutated)
    # then we actually need to fit rather than update (which will
    # force build of a new learning network). If `model` has been
    # created using a learning network export macro, the test used
    # below is perfect. In any other case it is at least conservative:
    network_model_ids = objectid.(models(yhat))
    fields = [getproperty(model, name) for name in fieldnames(typeof(model))]
    submodels    = filter(f->f isa Model, fields)
    submodel_ids = objectid.(submodels)
    if !issubset(submodel_ids, network_model_ids)
        return fit(model, verbosity, args...)
    end

    is_anonymized = cache isa NamedTuple{(:sources, :data)}

    if is_anonymized
        sources, data = cache.sources, cache.data
        for k in eachindex(sources)
            rebind!(sources[k], data[k])
        end
    end

    fit!(yhat; verbosity=verbosity)
    if is_anonymized
        for s in sources
            rebind!(s, nothing)
        end
    end

    return yhat, cache, report(yhat)
end

# legacy method (replacements defined in operations.jl):
predict(::SupervisedComposite, fitresult::Node, Xnew)     = fitresult(Xnew)

# legacy method (replacements defined in operations.jl):
transform(::UnsupervisedComposite, fitresult::Node, Xnew) = fitresult(Xnew)


text:
function build_kalman_model(obj::IntegratedPropagatorBuilder, arg0::List, arg1::List, arg2::ParameterDriversList)
    return jcall(obj, "buildKalmanModel", KalmanODModel, (List, List, ParameterDriversList), arg0, arg1, arg2)
end

function build_ls_model(obj::IntegratedPropagatorBuilder, arg0::Vector{IntegratedPropagatorBuilder}, arg1::List, arg2::ParameterDriversList, arg3::ModelObserver)
    return jcall(obj, "buildLSModel", BatchLSODModel, (Vector{IntegratedPropagatorBuilder}, List, ParameterDriversList, ModelObserver), arg0, arg1, arg2, arg3)
end

function get_frame(obj::PropagatorBuilder)
    return jcall(obj, "getFrame", Frame, ())
end

function get_initial_orbit_date(obj::PropagatorBuilder)
    return jcall(obj, "getInitialOrbitDate", AbsoluteDate, ())
end

function get_orbit_type(obj::PropagatorBuilder)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_orbital_parameters_drivers(obj::PropagatorBuilder)
    return jcall(obj, "getOrbitalParametersDrivers", ParameterDriversList, ())
end

function get_position_angle(obj::PropagatorBuilder)
    return jcall(obj, "getPositionAngle", PositionAngle, ())
end

function get_propagation_parameters_drivers(obj::PropagatorBuilder)
    return jcall(obj, "getPropagationParametersDrivers", ParameterDriversList, ())
end

function get_selected_normalized_parameters(obj::PropagatorBuilder)
    return jcall(obj, "getSelectedNormalizedParameters", Vector{jdouble}, ())
end

function reset_orbit(obj::IntegratedPropagatorBuilder, arg0::Orbit)
    return jcall(obj, "resetOrbit", void, (Orbit,), arg0)
end


text:
using Raylib

const screenWidth = 800
const screenHeight = 450

function main()
    Raylib.InitWindow(
        screenWidth, screenHeight,
        "raylib [core] example - keyboard input"
    )

    ballColor = Raylib.DARKBLUE

    Raylib.SetTargetFPS(60)
    while !Raylib.WindowShouldClose()    # Detect window close button or ESC key
        ballPosition = Raylib.GetMousePosition()

        # update
        if Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_LEFT))
            ballColor = Raylib.MAROON
        elseif Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_MIDDLE))
            ballColor = Raylib.LIME
        elseif Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_RIGHT))
            ballColor = Raylib.DARKBLUE
        elseif Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_SIDE))
            ballColor = Raylib.PURPLE
        elseif Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_EXTRA))
            ballColor = Raylib.YELLOW
        elseif Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_FORWARD))
            ballColor = Raylib.ORANGE
        elseif Raylib.IsMouseButtonPressed(Int(Raylib.MOUSE_BUTTON_BACK))
            ballColor = Raylib.BEIGE
        end

        # draw
        Raylib.BeginDrawing()
        Raylib.ClearBackground(Raylib.RAYWHITE)
        Raylib.DrawText("move the ball with arrow keys", 10, 10, 20, Raylib.DARKGRAY)
        Raylib.DrawCircleV(ballPosition, 50, ballColor)
        Raylib.EndDrawing()
    end

    Raylib.CloseWindow()
end


text:
using TropicalGEMM, Octavian, LoopVectorization
using TropicalNumbers, LinearAlgebra
using VectorizationBase: VecUnroll, Vec, StaticInt
using TropicalGEMM: naive_mul!
using Test

function distance(a::AbstractArray{<:Tropical}, b::AbstractArray{<:Tropical})
    sum(abs.(content.(a) .- content.(b)))
end

_eps(::Type{Tropical{T}}) where T = eps(T)
_eps(::Type{Tropical{T}}) where T<:Integer = 0
_rand(::Type{T}, args...) where T = T <: Tropical{<:Integer} ? T.(rand(1:10, args...)) : T.(randn(args...))

macro test_close(a, b, atol)
    esc(
        :(@test isapprox(distance($a, $b), 0; atol=$atol))
    )
end

@testset "mul with static ints" begin
    @test StaticInt{0}() * Tropical(3.0) == Tropical(-Inf)
    @test StaticInt{1}() * Tropical(3.0) == Tropical(3.0)
    @test Tropical(3.0) * StaticInt{0}() == Tropical(-Inf)
    @test Tropical(3.0) * StaticInt{1}() == Tropical(3.0)
end

@testset "mydot" begin
    function mydot(a,b)
        s = zero(promote_type(eltype(a),eltype(b)))
        @avx for i in eachindex(a,b)
            s += a[i]*b[i]
        end
        s
    end

    a = Tropical.(randn(10))
    b = Tropical.(randn(10))

    @test mydot(a, b) ≈ transpose(a) * b
    @test LoopVectorization.check_args(TropicalF64, TropicalF64)
end

@testset "matmul" begin
    for (T1,T2) in [[TropicalF64, TropicalF64], [TropicalF32, TropicalF32], [TropicalF64, Tropical{Int64}], [Tropical{Int64}, Tropical{Int64}]]
        To = promote_type(T1, T2)
        atol = sqrt(max(_eps(T1), _eps(T2)))
        for n in [0, 1, 4, 40]  # 0 does not work yet!
            A = _rand(T1, n, n)
            B = _rand(T2, n, n)
            for (f, finplace) in [(Octavian.matmul_serial, Octavian.matmul_serial!),
                    (Octavian.matmul, Octavian.matmul!)]
                for tA in [true, false]
                    a = tA ? transpose(A) : A
                    @show T1,T2,n,f,tA
                    @test_close f(a, a) naive_mul!(similar(a), a, a) atol
                    for tB in [true, false]
                        b = tB ? transpose(B) : B
                        @show T1,T2,n,f,tA,tB
                        @test_close f(a, b) naive_mul!(similar(a), a, b) atol
                        α, β = _rand(To,2)
                        c = _rand(To, n, n) .|> T1
                        @test_close finplace(copy(c), a, b, α, β) naive_mul!(copy(c), a, b, α, β) atol
                        sa = view(a, 1:min(n, 2), :)
                        sb = view(b, :, 1:min(n,2))
                        c = _rand(To, min(n,2), min(n,2))
                        @test_close finplace(copy(c), sa, sb, α, β) naive_mul!(copy(c), sa, sb, α, β) atol
                    end
                end
            end
        end
    end
end

@testset "*" begin
    a = CountingTropical{Float64}.(randn(5,5))
    b = CountingTropical{Float64}.(randn(5,50))
    @test a * b == naive_mul!(similar(b), a, b)
    a = Tropical{Float64}.(randn(5,5))
    b = Tropical{Float64}.(randn(5,50))
    @test_close a * b naive_mul!(similar(b), a, b) 1e-12
end

@testset "fix julia-1.5" begin
    x=Tropical(Vec(1.0, 2.0))
    @test VecUnroll((x, x)) === Tropical(VecUnroll((Vec(1.0, 2.0), Vec(1.0, 2.0))))
end

@testset "fix nan bug" begin
    res = LinearAlgebra.mul!(Tropical.(fill(NaN, 2, 2)), transpose(Tropical.(randn(2,2))), Tropical.(randn(2,2)), 1, 0)
    @test !any(isnan, res)
end


text:
function F_dimerdecay(x, parms)
  S1, S2, S3 = x
  c1, c2, c3, c4 = parms

  R1 = c1 * S1
  R2 = c2 * S1 * (S1 - 1) / 2
  R3 = c3 * S2
  R4 = c4 * S2

  return [R1, R2, R3, R4]
end

function gillespie_dimerdecay(;
  x = [1_000, 0, 0], c = [1.0, 0.002, 0.5, 0.04]
  )
  # initial conditions
  x0 = copy(x)
  # stoichiometries
  nu = [
    [-1  0  0]; # S1 --> 0
    [-2  1  0]; # 2 * S1 --> S2
    [ 2 -1  0]; # S2 --> 2 * S1
    [ 0 -1  1]  # S2 --> S3
  ]
  # parameter vector
  parms = copy(c)
  return x0, nu, parms
end

x0, nu, parms = gillespie_dimerdecay()
F = F_dimerdecay


text:
# Number of iterations Euclidean Algorithm needs for gcd(a, b)
function gcd_iters(a, b)
    iters = 0
    while b != 0
        iters += 1
        a, b = b, a % b
    end
    return iters
end

# Number of bits needed to represent some non-negative n
function num_bits(n)
    ct = 1
    n >>= 1
    while n != 0
        n >>= 1
        ct += 1
    end
    return ct
end

# Test num_bits
num_bits_tests = [(0, 1), (1, 1), (2, 2), (3, 2), (4, 3), (5, 3), (6, 3), (7, 3), (8, 4)]
for (num, ans) in num_bits_tests
    @assert num_bits(num) == ans
end

small_numbers = collect(0:2^12-1)
fibs = [0, 1] # Fibonacci numbers
for i in 1:90
    push!(fibs, fibs[length(fibs) - 1] + fibs[length(fibs)])
end
two_powers = [2 ^ x for x in 1 : 60]
mersenne_numbers = [x - 1 for x in two_powers]

operands_of_interest = Set(vcat(small_numbers, fibs, two_powers, mersenne_numbers))

# Verify that for all operands of interest, num iters is <= 1 + num_bits(y) ÷ log2_of_φ (note floor divison)
log2_of_φ = log2(MathConstants.golden)
for x in operands_of_interest, y in operands_of_interest
# for x in 0:2^16-1, y in 0:2^16-1
    iters = gcd_iters(x, y)
    if 1 + num_bits(y) ÷ log2_of_φ < iters
        println("gcd($(x), $(y)) takes $(iters) iters (max bits: $(max_bits), ratio: $(iters / max_bits))")
    end
    @assert iters <= 1 + num_bits(y) ÷ log2_of_φ
end
println("Done")

#==
Operands where (# of iterations)/(max_bits) = 1.5 (none have a higher ratio, up through 14 bits)
Note that all are sequential fibonacci #s
gcd(2, 3) takes 3 iters (max bits: 2)
gcd(8, 13) takes 6 iters (max bits: 4)
gcd(34, 55) takes 9 iters (max bits: 6)
gcd(144, 233) takes 12 iters (max bits: 8)
gcd(610, 987) takes 15 iters (max bits: 10)

https://stackoverflow.com/questions/9060816/running-time-of-euclids-gcd-algorithm
y = second arg
b = # of bits in second arg
bound is:
1 + log(φ, y)
1 + log(φ, 2^b - 1)
1 + log(2, 2^b - 1)/log(2, φ)
1 + b/log(2, φ)

experimentally verified for all numbers representable in 16 bits

==#


text:
module ins

using LightGraphs
using JuliaDB
using JuliaDBMeta
using TensorDecompositions: CANDECOMP
using IndexedTables: IndexedTable

using StatsBase: countmap

import SparseArrays
import Dates
import Plots

#= struct D =#
#=     indegree::Any, =#
#=     outdegree::Any, =#
#=     time::Any, =#
#= end =#

mutable struct DGraph
    meta::IndexedTable
    adj_mats::Array{SparseArrays.SparseMatrixCSC{Int64,Int64}}
    D::CANDECOMP
    g::LightGraphs.AbstractGraph
    n::Int
    r::Int
    nlookup::Vector{Int}
    projmeta

    DGraph(meta, adj_mats, D) = begin
        n, r = size(D.factors[1])
        g = _digraph(adj_mats[end])
        nlookup = (@groupby meta :Project_Node {ID=:Project_ID[1]}) |>
            @select :ID
        new(meta, adj_mats, D, g, n, r, nlookup, missing)
    end
end

"I made the graphs backwards in ProcessLibrariesIO.jl so reverse that."
_digraph(adj_mat) = DiGraph(adj_mat |> permutedims)

months(dg) = begin
    months = select(dg.meta, :Month)
    minimum(months):Dates.Month(1):maximum(months)
end

activity(dg, component, factor=2) =
    dg.D.factors[factor][1:end,component]

node_viewer(dg, component, nodes) = begin
    # TODO: Look at activities across all components in factor 2 instead.
    nodes = sort(nodes)
    activities = [activity(dg, component, f)[nodes] for f in 1:2]
    res = select(dg.projmeta[nodes], (:Name, :SourceRank))
    #= res = @transform_vec res ( =#
    #=     Activity_k=activities[2], =#
    #=     Activity_h=activities[1], =#
    #=     Node=nodes) =#
    res = setcol(res, :Activity_k, activities[2])
    res = setcol(res, :Activity_h, activities[1])
    res = reindex(res, :SourceRank)
    sort(res, rev=true)
end

contributing_nodes(dg, component, pct=2) = begin
    nodes = top_percentile(dg, component, pct)
    node_viewer(dg, component, nodes)
end

function rank(dg::DGraph, component)
    fs = dg.D.factors
    ranks = [sortperm(f[1:end,component], rev=true) for f in fs[1:2]]
end

function top_n(dg::DGraph, component, n)
    ranks = rank(dg, component)
    ranks = ranks[2:2] # ignore factor 1
    vcat([rank[1:min(n, length(rank))] for rank in ranks]...) |> unique
end

function top_percentile(dg::DGraph, component, pct)
    top_n(dg, component, Int(dg.n ÷ (100/pct)))
end

function links_over_time(dg)
    table((Month=ins.months(dg), nlinks=[SparseArrays.nnz(a) for a in dg.adj_mats]))
end

function get_projmeta(dg, platform::AbstractString = platform(dg))
    projmeta = loadtable("../data/sample-1.4/$(platform)_projects_with_repository_fields-1.4.0-2018-12-22.csv")
    get_projmeta(dg, projmeta)
end

function get_projmeta(dg, projmeta)
    project_ids = vcat(
        select(dg.meta, :Project_ID),
        select(dg.meta, :Dependency_Project_ID)) |>
        unique
    node_numbers = Dict(project_ids[i] => i for i in 1:length(project_ids))
    projmeta = setcol(projmeta, :Node, :ID => pid->getkey(node_numbers, pid, missing))
    dg.projmeta = reindex(projmeta, :Node)
end

function first_appearance(dg, node)
    (@filter dg.meta :Project_Node == node) |>
        @with minimum(:Month)
end

function first_appearances(dg)
    (@groupby dg.meta :Project_Node minimum(:Month)) |> @select _[2]
end

function platform(dg)
    platform_freq = dg.meta.columns.Platform |> countmap |> collect
    sort(platform_freq, by = x -> x[2])[end][1]
end

function name_pkg(dg::DGraph, node::Int)
    if ismissing(dg.projmeta) || node > length(dg.projmeta)
        (@filter dg.meta node in (:Project_Node, :Dependency_Project_Node)) |>
            @with ((node == :Project_Node[1] ? :Project_Name : :Dependency_Name) |> unique)[1]
    else
        dg.projmeta[node].Name
    end
end

Base.display(dg::DGraph) = begin
    println("""DGraph(platform = $(platform(dg)), n = $(dg.n), r = $(dg.r))

    Top packages degree:

        $(join(
            map(n -> name_pkg(dg ,n) * "  --  $n", sortperm(indegree(dg.g), rev=true)[1:5]),
            "\n        "))""")
end

"Component activity over time"
plot_time(dg, nox=false) = begin
    if nox
        #= cs = dg.D.factors[3] |> eachcol |> collect =#
        #= plt = UnicodePlots.lineplot(cs[1]) =#
        #= [UnicodePlots.lineplot!(plt, cs[c]) for c in 2:length(cs)] =#
        #= plt =#
        Plots.unicodeplots()
        Plots.plot(
            dg.D.factors[3],
            legend=:topleft,
            xlabel="Time",
            size=(600,400),
        )
    else
        Plots.plot(
            dg.D.factors[3],
            legend=:topleft,
            xlabel="Time",
            ylabel="Activity",
        )
    end
end

"Component heatmaps by degree"
plot_degree(dg, factor, nox=false) = begin
    nox ? Plots.gr() : missing
    plt = Plots.heatmap(dg.D.factors[factor], xlabel="Component", ylabel="Project", size=(1366,768))
    if nox
        filename = "/tmp/heatmap-$(rand()).png"
        Plots.savefig(plt, filename)
        run(`fake_browser_copy_to $filename`)
    else
        plt
    end
end

end


text:
export logψ!

# Definitions for batched evaluation of networks
# When the networrks are not cached and therefore allocate
# the result structure
@inline logψ!(out::AbstractMatrix, net::NeuralNetwork, σ::ADoubleStateBatch) where N =
    out .= net(σ...)
@inline logψ!(out::AbstractMatrix, net::NeuralNetwork, σ::AStateBatch) =
    out .= net(σ)

"""
    NNCache{N}

The base abstract type holding a cache for the network type `N<:NeuralNetwork`.
"""
abstract type NNBatchedCache{N} <: NNCache{N}
    #valid::Bool
end

cached(net::NeuralNetwork, batch_sz::Int) =
    CachedNet(net, cache(net, batch_sz))
cached(net::CachedNet, batch_sz::Int) =
    CachedNet(net.net, cache(net.net, batch_sz))

batch_size(net::CachedNet) = batch_size(net.cache)
batch_size(cache::NNBatchedCache) = throw("Not Implemented")

out_similar(net::CachedNet{N,C}) where{N,C<:NNBatchedCache} =
    similar(trainable_first(net), out_type(net), 1, batch_size(net))
out_similar(net::CachedNet{N,C}, dims::Vararg{T,I}) where{T<:Integer, I, N,C<:NNBatchedCache} =
    similar(trainable_first(net), out_type(net), 1, dims...)

# Definition for inplace evaluation of batched cached networks
@inline logψ!(out::AbstractMatrix, net::CachedNet, σ::Vararg{T,N}) where {N,T<:AbstractVecOrMat} =
    logψ!(out, net.net, net.cache, σ...)
@inline logψ!(out::AbstractMatrix, net::CachedNet, σ::ADoubleStateBatch) =
    logψ!(out, net.net, net.cache, row(σ), col(σ))
@inline logψ!(out::AbstractMatrix, cnet::CachedNet, σ::AStateBatch) =
    logψ!(out, cnet.net, cnet.cache, σ)

# Definition for allocating evaluation of batched cached networks
# Shadowing things at ~80 of base_cached_networks.jl
@inline logψ(net::CachedNet{NN,NC}, σr::AStateBatch, σc::AStateBatch) where {N,T,NN,NC<:NNBatchedCache} = begin
    logψ(net, (σr, σc))
end
@inline logψ(net::CachedNet{NN,NC}, σ::ADoubleStateBatch) where {N,T,NN,NC<:NNBatchedCache} = begin
    b_sz = last(size(first(σ)))
    out = similar(trainable_first(net), out_type(net), 1, b_sz)
    logψ!(out, net.net, net.cache, σ...)
end
@inline logψ(cnet::CachedNet{NN,NC}, σ::AStateBatch) where {NN,NC<:NNBatchedCache} = begin
    b_sz = last(size(σ))
    out = similar(trainable_first(cnet), out_type(cnet), 1, b_sz)
    logψ!(out, cnet.net, cnet.cache, σ)
end

# Declare the two functions, even if config(blabla)=blabla, because of a shitty
# Julia's performance bug #32761
# see https://github.com/JuliaLang/julia/issues/32761
@inline logψ_and_∇logψ!(der, n::CachedNet{NN,NC}, σ::AbstractArray) where {NN,NC<:NNBatchedCache} = begin
    b_sz = last(size(σ))
    out = similar(trainable_first(n), out_type(n), 1, b_sz)
    logψ_and_∇logψ!(der, out, n.net, n.cache, σ)
    return out, der
end
@inline logψ_and_∇logψ!(der, n::CachedNet{NN,NC}, σ::NTuple{N,<:AbstractArray}) where {N,NN,NC<:NNBatchedCache} = begin
    b_sz = last(size(first(σ)))
    out = similar(trainable_first(n), out_type(n), 1, b_sz)
    logψ_and_∇logψ!(der, out, n.net, n.cache, σ...)
    return out, der
end
@inline logψ_and_∇logψ!(der, n::CachedNet{NN,NC}, σ::Vararg{<:AbstractArray,N}) where {N,T,NN,NC<:NNBatchedCache} = begin
    b_sz = last(size(first(σ)))
    out = similar(trainable_first(n), out_type(n), 1, b_sz)
    logψ_and_∇logψ!(der, out, n.net, n.cache, σ...)
    return out, der
end

@inline function logψ_and_∇logψ!(der::AbstractDerivative, out, n::CachedNet, σ::NTuple{N,<:AbstractArray}) where N
    logψ_and_∇logψ!(der, out, n.net, n.cache, σ...)
    return (out, der)
end
@inline function logψ_and_∇logψ!(der::AbstractDerivative, out, n::CachedNet, σ::Vararg{<:AbstractArray,N}) where N
    logψ_and_∇logψ!(der, out, n.net, n.cache, σ...)
    return (out, der)
end
@inline function logψ_and_∇logψ!(der::AbstractDerivative, out, n::CachedNet, σ::AbstractArray) where N
    logψ_and_∇logψ!(der, out, n.net, n.cache, σ);
    return (out, der)
end

#
grad_cache(net::NeuralNetwork, batch_sz) =
    grad_cache(out_type(net), net, batch_sz)
function grad_cache(T::Type{<:Number}, net::NeuralNetwork, batch_sz)
    is_analytic(net) && return RealDerivative(T, net, batch_sz)
    return WirtingerDerivative(T, net, batch_sz)
end

#
grad_cache(T::Type{<:Number}, net::CachedNet{N,C}) where {N,C<:NNBatchedCache} =
    grad_cache(T, net, batch_size(net))

function RealDerivative(T::Type{<:Number}, net::NeuralNetwork, batch_sz::Int)
    pars = trainable(net)

    vec       = similar(trainable_first(pars), T, _tlen(pars), batch_sz)
    i, fields = batched_weight_tuple(net, vec)
    return RealDerivative(fields, [vec])
end

#
"""
    grad_cache_vec(net, batch_sz, vec_len)

Builds a vector of `grad_cache` batching `batch_sz` evaluations.
All share the same underlying vector.
"""
grad_cache(net::NeuralNetwork, batch_sz::Integer, vec_len::Integer) =
    grad_cache(out_type(net), net, batch_sz, vec_len)

function grad_cache(T::Type{<:Number}, net::NeuralNetwork, batch_sz, vec_len)
    pars  = trainable(net)
    npars = _tlen(pars)

    if is_analytic(net)
        vec       = similar(trainable_first(pars), T, npars, batch_sz, vec_len)
        # compute type
        j, fields = batched_weight_tuple(net, view(vec, :, :, 1))
        der       = RealDerivative(fields, (view(vec, :, :, 1),))
        ders      = Vector{typeof(der)}()
        for i=1:vec_len
            vec_rsp   = view(vec, :, :, i)
            j, fields = batched_weight_tuple(net, vec_rsp)
            der       = RealDerivative(fields, (vec_rsp,))
            push!(ders, der)
        end
        vecs = (vec,)
    else
        throw("To implement")
    end
    return ders, vecs
end

@inline store_state!(cache::Array,
             v::AbstractVector,
             i::Integer) = begin
    #@uviews cache v begin
        uview(cache, :, i) .= v
    #end
    return cache
end

@inline store_state!(cache::NTuple{2,<:Matrix},
             (vl, vr)::NTuple{2,<:AbstractVector},
             i::Integer) = begin
    cache_l, cache_r = cache
    #@uviews cache_l cache_r vl vr begin
        uview(cache_l, :,i) .= vl
        uview(cache_r, :,i) .= vr
    #end
    return cache
end

@inline @inbounds unsafe_get_batch(σ::AStateBatch, i::Integer) = uview(σ, :, i)
@inline @inbounds unsafe_get_batch(σ::ADoubleStateBatch, i::Integer) = (uview(row(σ), :, i), uview(col(σ), :, i))
@inline num_batches(σ::AStateBatch) = size(σ, 2)
@inline num_batches(σ::ADoubleStateBatch) = size(row(σ), 2)

# Automatic evaluation of batches
logψ!(out::AbstractArray{T,3}, cnet::CachedNet, σ::NTuple{2,AbstractArray{T2,3}}) where{T,T2} =
    logψ!(out, cnet, first(σ), last(σ))
function logψ!(out::AbstractArray{T,3}, cnet::Union{NeuralNetwork, CachedNet}, σr::AT, σc::AT) where {T, T2, AT<:AbstractArray{T2, 3}}
    n_batches = size(σr, 3)
    for i=1:n_batches
        σr_i  = state_uview(σr,  i)
        σc_i  = state_uview(σc,  i)
        out_i = state_uview(out, i)
        logψ!(out_i, cnet, σr_i, σc_i)
    end
    return out
end

function logψ!(out::AbstractArray{T,3}, cnet::NeuralNetwork, σ::AT) where {T, T2, AT<:AbstractArray{T2, 3}}
    n_batches = size(σ, 3)
    for i=1:n_batches
        σ_i   = state_uview(σ,  i)
        out_i = state_uview(out, i)
        logψ!(out_i, cnet, σ_i)
    end
    return out
end

# Compute batches of those things
function logψ_and_∇logψ!(der::Vector{<:AbstractDerivative}, out::AbstractArray{T,3},
                         net::Union{NeuralNetwork, CachedNet},
                         σ::Union{AStateBatchVec, ADoubleStateBatchVec}) where {T}
    n_batches  = length(der)
    for i=1:n_batches
        σ_i   = state_uview(σ,  i)
        out_i = state_uview(out, i)
        logψ_and_∇logψ!(der[i], out_i, net, σ_i)
    end
    return out
end

function logψ!(out::AbstractArray{T,3}, net::Union{NeuralNetwork, CachedNet},
               σ::Union{AStateBatchVec, ADoubleStateBatchVec}) where {T}
    n_batches  = length(der)
    for i=1:n_batches
        σ_i   = state_uview(σ,  i)
        out_i = state_uview(out, i)
        logψ!(out_i, net, σ_i)
    end
    return out
end

function logψ_and_∇logψ(net::Union{NeuralNetwork, CachedNet},
                        σ::Union{AStateBatchVec, ADoubleStateBatchVec}) where {T}
    ∇vals, ∇vec    = grad_cache(bnet, 1, batch_size(σ), chain_length(σ))
    out = similar(out_similar(net), batch_size(σ), chain_length(σ))
    logψ_and_∇logψ!(∇vals, out, net, σ)
    return ∇vals, out
end

function logψ(net::Union{NeuralNetwork, CachedNet},
               σ::Union{AStateBatch, ADoubleStateBatch}) where {T}
    out = out_similar(net, batch_size(σ))
    logψ!(out, net, σ)
end

function logψ(net::Union{NeuralNetwork, CachedNet},
               σ::Union{AStateBatchVec, ADoubleStateBatchVec}) where {T}
    out = similar(out_similar(net), 1, batch_size(σ), chain_length(σ))
    logψ!(out, net, σ)
end

# diagonal
logψ!(out::AbstractArray, net::MatrixNeuralNetwork, cache::NNCache, σ::AStateBatch) =
    logψ!(out, net, cache, σ, σ)

log_prob_ψ!(prob::AbstractArray, net::NeuralNetwork, σ) =
    log_prob_ψ!(prob, prob, net, σ)

function log_prob_ψ!(prob::AbstractArray, out_tmp::AbstractArray, net::MatrixNet, σ::AStateOrBatchOrVec)
    logψ!(out_tmp, net, σ, σ)
    prob .= abs.(out_tmp)
    return prob
end

function log_prob_ψ!(prob::AbstractArray, out_tmp::AbstractArray, net::NeuralNetwork, σ)
    logψ!(out_tmp, net, σ)
    prob .= 2 .* real.(out_tmp)
    return prob
end


text:
using Test
using StaticArrays
using LinearAlgebra
using BenchmarkTools

using iLQGames:
    LinearSystem,
    AffineStrategy,
    QuadraticPlayerCost,
    LTVSystem,
    LQGame,
    SystemTrajectory,
    solve_lq_game!,
    n_states,
    n_controls,
    n_players,
    uindex,
    horizon,
    samplingtime,
    dynamics,
    strategytype,
    player_costs,
    trajectory!

function generate_1D_pointmass_game()
    # Testing the solver at a simple example: A two-player point mass 1D system.
    # The state composes of position and oritation. Therefore, the system dynamics
    # are a pure integrator.
    ΔT = 0.1
    H = 10.0
    N_STEPS = Int(H / ΔT)

    # dynamical system
    A = SMatrix{2, 2}([1. ΔT; 0. 1.])
    B = SMatrix{2, 2}([0.5*ΔT^2 ΔT; 0.32*ΔT^2 0.11*ΔT]')
    dyn = LinearSystem{ΔT}(A, B)
    # costs for each player
    c1 = QuadraticPlayerCost(@SVector([0., 0.]),       # l
                             @SMatrix([1. 0.; 0. 1.]), # Q
                             @SVector([0., 0.]),       # r
                             @SMatrix([1. 0.; 0. 0.])) # R
    c2 = QuadraticPlayerCost(-c1.l,                    # l
                             -c1.Q,                    # Q
                             @SVector([0., 0.]),       # r
                             @SMatrix([0. 0.; 0. 1.])) # R

    costs = @SVector [c1, c2]
    # the lq game (player one has control input 1 and 2; player 2 has control input 3
    ltv_dyn = LTVSystem(SizedVector{N_STEPS}(repeat([dyn], N_STEPS)))
    qtv_costs = SizedVector{N_STEPS}(repeat([costs], N_STEPS))
    uids = (SVector(1), SVector(2))
    lqGame = LQGame(uids, ltv_dyn, qtv_costs)

    # test all the function calls:
    @test n_players(lqGame) == length(costs)
    @test horizon(lqGame) == N_STEPS
    uindex(lqGame)

    return lqGame
end

function generate_2D_pointmass_game()
    ΔT = 0.1
    HORIZON = 10.0
    N_STEPS = Int(HORIZON / ΔT)

    # dynamical system
    # state layout: x = (x, y, vx, vy)
    A = I + ΔT * @SMatrix [0. 0. 1. 0.;
                           0. 0. 0. 1.;
                           0. 0. 0. 0.;
                           0. 0. 0. 0.]
    B = ΔT * @SMatrix [0. 0. 1. 0.
                       0. 0. 0. 1.
                       1. 0. 0. 0.
                       0. 1. 0. 0.]
    dyn = LinearSystem{ΔT}(A, B)

    # cost
    # player 1 want's the position to be zero
    c1 = QuadraticPlayerCost(@SVector(zeros(4)),      # l
                             @SMatrix([1. 0. 0. 0.;
                                       0. 2. 0. 0.;
                                       0. 0. 0. 0.;
                                       0. 0. 0. 0.]), # Q
                             @SVector(zeros(4)),      # r
                             @SMatrix([1. 0. 0. 0.;
                                       0. 1. 0. 0.;
                                       0. 0. 0. 0.;
                                       0. 0. 0. 0.]))# R)
    # cost
    # player 2 is adversarial
    c2 = QuadraticPlayerCost(-c1.l,                   # Q
                             -c1.Q ,                  # l
                             @SVector(zeros(4)),      # r
                             @SMatrix([0. 0. 0. 0.;
                                       0. 0. 0. 0.;
                                       0. 0. 1. 0.;
                                       0. 0. 0. 1.]))# R)
    costs = SVector(c1, c2)

    ltv_dyn = LTVSystem(SizedVector{N_STEPS}(repeat([dyn], N_STEPS)))
    qtv_costs = SizedVector{N_STEPS}(repeat([costs], N_STEPS))
    uids =  (SVector(1, 2), SVector(3, 4))
    lqGame = LQGame(uids, ltv_dyn, qtv_costs)

    return lqGame
end

# solve two-player inifinite horizion (time-invariant) LQ game by Lyapunov
# iterations
function solve_lyapunov_iterations(dyn::LinearSystem, c1::QuadraticPlayerCost, c2::QuadraticPlayerCost, ud1, ud2, n_iter=100)
    A = dyn.A
    B1 = dyn.B[:, ud1]
    B2 = dyn.B[:, ud2]
    Q1 = c1.Q
    Q2 = c2.Q
    R11 = c1.R[ud1, ud1]
    R12 = c1.R[ud2, ud2]
    R21 = c2.R[ud1, ud1]
    R22 = c2.R[ud2, ud2]

    # Initialize cost-to-go with terminal state cost
    Z1 = Q1
    Z2 = Q2

    # Initialize controls
    P1 = (R11 + B1' * Z1 * B1) \ (B1' * Z1 * A)
    P2 = (R22 + B2' * Z2 * B2) \ (B2' * Z2 * A)

    for ii in 1:n_iter
        P1_old = P1
        P2_old = P2

        P1 = (R11 + B1' * Z1 * B1) \ (B1' * Z1 * (A - B2 * P2_old))
        P2 = (R22 + B2' * Z2 * B2) \ (B2' * Z2 * (A - B1 * P1_old))

        Z1 = (A - B1 * P1 - B2 * P2)' * Z1 * (A - B1 * P1 - B2 * P2) + P1' * R11 * P1 + P2' * R12 * P2 + Q1
        Z2 = (A - B1 * P1 - B2 * P2)' * Z2 * (A - B1 * P1 - B2 * P2) + P1' * R21 * P1 + P2' * R22 * P2 + Q2
    end

    return (P1, P2)
end

function test_lyapunov(g::LQGame, strategies)

    # 1. Lyapunov solution for the inifnite horizion problem
    P1_lyap, P2_lyap = solve_lyapunov_iterations(dynamics(g)[1],
                                                 player_costs(g)[1][1],
                                                 player_costs(g)[1][2],
                                                 uindex(g)[1],
                                                 uindex(g)[2])
    P_lyap = [P1_lyap; P2_lyap]

    # 2. LQ game solution
    solve_lq_game!(strategies, g)
    γ0 = first(strategies)
    P_lqg = γ0.P

    print("""
          P_lyap: $(round.(P_lyap; digits=2))
          P_lqg:  $(round.(P_lqg; digits=2))
          """)

    @testset "lyapunov" begin
        @test isapprox(P_lyap, P_lqg; norm=v->norm(v, Inf), atol=1e-4)
    end;

    return P_lyap, P_lqg
end

@testset "solve_lq_game" begin
    g1D = generate_1D_pointmass_game()
    strategies1D = strategytype(g1D)(undef)
    test_lyapunov(g1D, strategies1D)
    (@benchmark solve_lq_game!($strategies1D, $g1D)) |> display

    g2D = generate_2D_pointmass_game()
    strategies2D = strategytype(g2D)(undef)
    test_lyapunov(g2D, strategies2D)
    (@benchmark solve_lq_game!($strategies2D, $g2D)) |> display
end;

"--------------- manual test to run and check for sanity ---------------"

function sim()
    g = generate_2D_pointmass_game()
    strategy = solve_lq_game(g)

    nx = n_states(g)
    nu = n_controls(g)
    ΔT = samplingtime(g)
    h = horizon(g)

    # dummy operating point
    last_op = zero(SystemTrajectory{h, ΔT, nx, nu})
    traj = deepcopy(last_op)
    x0 = @SVector randn(nx)
    trajectory!(traj, dynamics(g), strategy, last_op, x0)
end


text:
#---------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Model for tank volume calculation
#*----------------------------------------------------------------------
#*
#* 	Routine to calculate tank volume and level tank volume from
#* different geometries and orientations.
#*
#*	Geometry:
#*		* Flat (no head) (default)
#*		* Spherical
#*
#*	Orientation:
#*		* Vertical (default)
#*		* Horizontal
#*
#*----------------------------------------------------------------------
#* Author: Rodolfo Rodrigues
#* $Id$
#*--------------------------------------------------------------------
type vol_tank
	vol_tank()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"Pi value",
				:Default=>3.141593,
				:Symbol=>"\\pi"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tank head type",
				:Valid=>["flat","spherical"],
				:Default=>"flat"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tank orientation",
				:Valid=>["vertical","horizontal"],
				:Default=>"vertical"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Level tank volume"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Tank volume"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tank length"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tank level"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tank diameter"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Tank cross section area"
			)),
			[
				:(V = Across*Level),
				:(Vt = pi/4*(D^2)*L),
				:(V = pi/4*(D^2)*Level),
				:(Vt = pi/6*D^3),
				:(V = pi/3*Level^2*(3*D/2 - Level)),
				:(Vt = pi/4*(D^2)*L),
				:(V = ((D^2)*acos((D - 2*Level)/D)/4/"rad" - (D - 2*Level)*sqrt(D*Level - Level^2)/2)*L),
				:(Vt = pi/6*D^3),
				:(V = pi/3*Level^2*(3*D/2 - Level)),
			],
			[
				"Content volume","Tank volume","Level tank volume","Tank volume","Level tank volume","Tank volume","Level tank volume","Tank volume","Level tank volume",
			],
			[:pi,:Geometry,:Orientation,],
			[:V,:Vt,:L,:Level,:D,:Across,]
		)
	end
	pi::positive
	Geometry::DanaSwitcher
	Orientation::DanaSwitcher
	V::volume
	Vt::volume
	L::length
	Level::length
	D::length
	Across::area
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export vol_tank
function setEquationFlow(in::vol_tank)
	addEquation(1)
	let switch=Orientation
		if switch=="vertical"
			let switch=Geometry
				if switch=="flat"
					addEquation(2)
					addEquation(3)
				elseif switch=="spherical"
					addEquation(4)
					addEquation(5)
				end
			end
		elseif switch=="horizontal"
			let switch=Geometry
				if switch=="flat"
					addEquation(6)
					addEquation(7)
				elseif switch=="spherical"
					addEquation(8)
					addEquation(9)
				end
			end
		end
	end
end
function atributes(in::vol_tank,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Routine to calculate tank volume"
	fields[:Info]="
Based in 2 geometry and 2 configurations/orientations

== Geometry ==
* Flat (no head) (default)
* Spherical

== Orientation ==
* Vertical (default)
* Horizontal
"
	drive!(fields,_)
	return fields
end
vol_tank(_::Dict{Symbol,Any})=begin
	newModel=vol_tank()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
    map_set(::Type{BT}, set) where {BT}

Return the image of `set` through the linear map `A` defined in
[`Variable.SetMapBridge`](@ref) and [`Constraint.SetMapBridge`](@ref). This is
used for bridging the constraint and setting
the [`MathOptInterface.ConstraintSet`](@ref).
"""
function map_set end

"""
    inverse_map_set(::Type{BT}, set) where {BT}

Return the preimage of `set` through the linear map `A` defined in
[`Variable.SetMapBridge`](@ref) and [`Constraint.SetMapBridge`](@ref). This is
used for getting the [`MathOptInterface.ConstraintSet`](@ref).
"""
function inverse_map_set end

"""
    map_function(::Type{BT}, func) where {BT}

Return the image of `func` through the linear map `A` defined in
[`Variable.SetMapBridge`](@ref) and [`Constraint.SetMapBridge`](@ref). This is
used for getting the [`MathOptInterface.ConstraintPrimal`](@ref) of variable
bridges. For constraint bridges, this is used for bridging the constraint,
setting the [`MathOptInterface.ConstraintFunction`](@ref) and
[`MathOptInterface.ConstraintPrimalStart`](@ref) and
modifying the function with [`MathOptInterface.modify`](@ref).

    map_function(::Type{BT}, func, i::IndexInVector) where {BT}

Return the scalar function at the `i`th index of the vector function that
would be returned by `map_function(BT, func)` except that it may compute the
`i`th element. This is used by [`bridged_function`](@ref) and for getting
the [`MathOptInterface.VariablePrimal`](@ref) and
[`MathOptInterface.VariablePrimalStart`](@ref) of variable bridges.
"""
function map_function end

function map_function(::Type{BT}, func, i::IndexInVector) where {BT}
    return MOIU.eachscalar(map_function(BT, func))[i.value]
end

"""
    inverse_map_function(::Type{BT}, func) where {BT}

Return the image of `func` through the inverse of the linear map `A` defined in
[`Variable.SetMapBridge`](@ref) and [`Constraint.SetMapBridge`](@ref). This is
used by [`Variable.unbridged_map`](@ref) and for setting the
[`MathOptInterface.VariablePrimalStart`](@ref) of variable bridges
and for getting the [`MathOptInterface.ConstraintFunction`](@ref),
the [`MathOptInterface.ConstraintPrimal`](@ref) and the
[`MathOptInterface.ConstraintPrimalStart`](@ref) of constraint bridges.
"""
function inverse_map_function end

"""
    adjoint_map_function(::Type{BT}, func) where {BT}

Return the image of `func` through the adjoint of the linear map `A` defined in
[`Variable.SetMapBridge`](@ref) and [`Constraint.SetMapBridge`](@ref). This is
used for getting the [`MathOptInterface.ConstraintDual`](@ref) and
[`MathOptInterface.ConstraintDualStart`](@ref) of constraint bridges.
"""
function adjoint_map_function end

"""
    inverse_adjoint_map_function(::Type{BT}, func) where {BT}

Return the image of `func` through the inverse of the adjoint of the linear map
`A` defined in [`Variable.SetMapBridge`](@ref) and
[`Constraint.SetMapBridge`](@ref). This is used for getting the
[`MathOptInterface.ConstraintDual`](@ref) of variable bridges and setting the
[`MathOptInterface.ConstraintDualStart`](@ref) of constraint bridges.
"""
function inverse_adjoint_map_function end


text:
using ImageMagick, ColorVectorSpace, ImageMetadata, ImageTransformations

using Random: bitrand
using Base.CoreLogging: SimpleLogger, with_logger
using Pkg

function is_ci()
    get(ENV, "TRAVIS", "") == "true" ||
    get(ENV, "APPVEYOR", "") in ("true", "True") ||
    get(ENV, "CI", "") in ("true", "True")
end

@info "ImageMagick version ", ImageMagick.libversion()

include("constructed_images.jl")
include("readremote.jl")
include("badimages.jl")

workdir = joinpath(tempdir(), "Images")
try
    rm(workdir, recursive=true)
catch
end


text:
## ----------------  Individual moments

## Mean time to dropout (conditional on dropping out)
# Some small cells. 4Y only.
function time_to_drop_4y_qual_gpa(ds :: DataSettings)
    load_fct =
        mt -> read_matrix_by_xy(ds, :timeToDrop4y_qgM, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert size(m) == (n_4year(ds), n_parental(ds));
    @assert all(m .< 6.0)  &&  all(m .>= 0.0)
    return m, ses, cnts
end

## Fraction in each quality, conditional on entry, by gpa
function frac_qual_by_gpa(ds :: DataSettings)
    dataV = read_matrix_by_xy(raw_entry_qual_gpa(ds));
    @assert check_float_array(dataV, 0.001, 1.0);
    @assert isapprox(sum(dataV), 1.0, atol = 1e-6);

    # Make conditional on entry (columns sum to 1)
    dataV = dataV ./ sum(dataV, dims = 1);
    @assert all(isapprox.(sum(dataV, dims = 1), 1.0))
    nc = size(dataV, 1);

    # Counts by GPA
    cntV = read_row_totals(raw_entry_qual_gpa(ds; momentType = :count));
    @assert length(cntV) == size(dataV, 2)
    cnts = repeat(cntV', outer = (nc, 1));
    @assert size(cnts) == size(dataV)

    ses = ses_from_choice_probs(dataV, cnts);
    cnts = round.(Int, cnts);
    return dataV, ses, cnts
end

## Mass of freshmen by quality / gpa. Sums to 1.
function mass_entry_qual_gpa(ds :: DataSettings)
    m = read_matrix_by_xy(ds, :massEntry_qgM, MtMean);
    @assert size(m) == (n_colleges(ds), n_gpa(ds));
    # SES treats this as a discrete choice problem with total count.
    cnts = read_matrix_by_xy(ds, :massEntry_qgM, MtCount);
    cnts = fill(sum(cnts), size(cnts));
    @assert all(m .< 1.0)  &&  all(m .> 0.0)
    @assert isapprox(sum(m), 1.0,  atol = 0.0001)
    ses = ses_from_choice_probs(m, cnts);
    cnts = round.(Int, cnts);
    return m, ses, cnts
end

"""
	$(SIGNATURES)

Graduation rates by [quality, gpa]. All colleges.
"""
function grad_rate_qual_gpa(ds :: DataSettings)
    load_fct =
        mt -> read_matrix_by_xy(ds, :fracGrad_qgM, mt);
    m, ses, cnts = choice_prob_from_xy(load_fct);
    @assert all(m .<= 1.0)  &&  all(m .> 0.0);
    @assert size(m) == (n_colleges(ds), n_gpa(ds));

    # Set to 0 for colleges that do not produce graduates
    m[no_grad_idx(ds), :] .= 0.0;
    cnts[no_grad_idx(ds), :] .= 0;
    # wtM = 1.0 ./ max.(0.1, m);
    return m, ses, cnts
end

## --------------  Study time

"""
	$(SIGNATURES)

Study time from NLSY79. Class and study time combined.
Files are copied BY HAND from Dropbox to DataCollegeStrat.
But need to adjust b/c grand mean should be lower for NLSY97.
"""
function study_time_qual(ds :: DataSettings)
    studyV = [read_row_totals(study_time79_path(ds; momentType = mt))
        for mt in (:mean, :std, :count)];
    classV = [read_row_totals(class_time79_path(ds; momentType = mt))
        for mt in (:mean, :std, :count)];

    m, ses, cnts = total_study_time(studyV, classV);
    return m, ses, cnts
end

function study_time79_path(ds; momentType = :mean)
    fPath = data_file(raw_study_time_qual_gpa(ds; momentType = momentType));
    fPath = nlsy79_path(fPath);
    return fPath
end

function class_time79_path(ds; momentType = :mean)
    fPath = data_file(raw_class_time_qual_gpa(ds; momentType = momentType));
    fPath = nlsy79_path(fPath);
    return fPath
end

nlsy79_path(fPath) = replace(fPath, "97" => "79");

function total_study_time(studyV, classV)
    @assert size(studyV) == size(classV) == (3, );
    for j = 1 : 3
        @assert size(studyV[j]) == size(studyV[1]);
        @assert size(classV[j]) == size(studyV[1]);
    end

    # Multiply by mean of NSSE relative to NLSY79
    # Babcock/Marks 2011 tb 2
    studyFactor = 13.28 / 19.75;
    # Adjust to match HERI 2004 mean (NSSE does not have class time)
    classFactor = 13.01 / 15.84;

    cnts = round.(Int, 0.5 .* (studyV[3] .+ classV[3]));
    m = studyV[1] .* studyFactor .+ classV[1] .* classFactor;
    ses = (studyV[2] .* studyFactor .+ classV[2] .* classFactor) ./
        sqrt.(cnts);
    return m, ses, cnts
end

function study_time_qual_gpa(ds :: DataSettings)
    studyV = [read_matrix_by_xy(study_time79_path(ds; momentType = mt))
        for mt in (:mean, :std, :count)];
    classV = [read_matrix_by_xy(class_time79_path(ds; momentType = mt))
        for mt in (:mean, :std, :count)];

    dataM, ses, cnts = total_study_time(studyV, classV);

    # fPath = data_file(raw_study_time_qual_gpa(ds; momentType = :mean));
    # fPath = replace(fPath, "97" => "79");
    # m = read_row_totals(fPath);

    # # Study time only, by quality, gpa
    # studyM = [10.7	12.9	12.5	16.2;
    #         19.3	21.6	18.2	15.2;
    #         19.3	13.1	19.5	20.2;
    #         15.1	24.2	28.5	24.8];
    # # Multiply by mean of NSSE relative to NLSY79
    # # Babcock/Marks 2011 tb 2
    # dFactor = 13.28 / 19.75;
    # studyM = studyM .* dFactor;

    # # Class time
    # classM = [13.5	16.0	13.4	17.3;
    #         17.7	18.3	18.3	13.9;
    #         20.1	17.5	15.2	15.9;
    #         32.6	15.9	16.3	16.1];
    # # Adjust to match HERI 2004 mean (NSSE does not have class time)
    # dFactor = 13.01 / 15.84;
    # classM = classM .* dFactor;

    # dataM = studyM .+ classM;
    # # if modelUnits
    # #     dataM = hours_per_week_to_mtu(dataM);
    # #     validate_mtu(dataM);
    # # end

    # # These are made up
    # ses = ones(size(dataM));
    # cnts = fill(100, size(dataM)...);

    @check size(dataM) == (n_colleges(ds), n_gpa(ds))

    return dataM, ses, cnts
end

function study_time(ds :: DataSettings)
    studyV = [read_total(study_time79_path(ds; momentType = mt))
        for mt in (:mean, :std, :count)];
    classV = [read_total(class_time79_path(ds; momentType = mt))
        for mt in (:mean, :std, :count)];

    m, ses, cnts = total_study_time(studyV, classV);
    return m, ses, cnts
end

# function nlsy79_dir(ds)
#     @assert occursin(data_dir(ds), "97")  "Expecting 97 NLSY";
#     d2 = replace(data_dir(ds), "97" => "79");
#     return d2
# end

# study_time_dir79(ds) = joinpath(nlsy79_dir(ds), "dat_files");

# function study_time_qual_gpa(ds :: DataSettings)
#     dataM = load_study_time_qual_gpa(ds; modelUnits = true);
#     return Deviation{Double}(name = :studyTime_qgM, dataV = dataM, modelV = dataM,
#         scalarWt = 1.0 / sum(dataM),
#         shortStr = "studyTimeQualGpa",
#         longStr = "Average study time by (quality, gpa)",
#         showPath = "studyTimeQualGpa");
# end

# In hours per week
# function mean_study_times_by_qual(ds :: DataSettings)
#     st_qgM, _ = study_time_qual_gpa(ds);
#     st_qV = [mean_by_gpa(vec(st_qgM[ic,:]), ds)  for ic = 1 : n_colleges(ds)];
#     @check size(st_qV) == (n_colleges(ds), )
#     @assert all(st_qV .> 7.0)
#     ses = ones(size(st_qV));
#     cnts = fill(100, size(st_qV));
#     return st_qV, ses, cnts
# end

# In hours per week
# function mean_study_time(ds :: DataSettings)
#     st_qV, _ = mean_study_times_by_qual(ds);
#     st = mean_by_qual(st_qV, ds);
#     return st, 1.0, 100
# end

# This is not a useful target. High quality / low AFQT cells are almost empty.
# Instead, use marginals (by quality, AFQT) computed for 4y colleges only.
function time_to_grad_4y_qual_gpa(ds :: DataSettings)
    load_fct =
        mt -> read_matrix_by_xy(ds, :timeToGrad4y_qgM, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert check_float_array(m, 3.0, 7.0);
    # 4y only
    @assert size(m) == (n_4year(ds), n_gpa(ds))
    # # zero out 2 year colleges
    # m[two_year_colleges(ds), :] .= 0.0;
    # cnts[two_year_colleges(ds), :] .= 0;
    return m, ses, cnts
end

# -----------------


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

using Markdown
using InteractiveUtils

# ╔═╡ 29b26c10-fd63-11ea-19d9-51bbbdf5e6f6
using LinearAlgebra

# ╔═╡ 15626ee0-fd68-11ea-2c83-2dc85ba373e7
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 973991a0-fd68-11ea-27c8-61246753b4b3
ingredients("../AcousticPropagation.jl")

# ╔═╡ ce78547e-fd68-11ea-3d71-61ec6c934761
t_inc = [0.0014855627054164151, 0.003713906763541037]

# ╔═╡ e64e2ad0-fd68-11ea-040d-6907d34eca35
t_bnd = [1.0, 0.6770329614269008]

# ╔═╡ ef543200-fd68-11ea-1c37-c158193f71be
n_bnd = [-t_bnd[2], t_bnd[1]]

# ╔═╡ 0bf32920-fd69-11ea-293c-4949b8fee869
t_rfl = t_inc - 2LinearAlgebra.dot(t_inc, n_bnd)*n_bnd

# ╔═╡ Cell order:
# ╠═29b26c10-fd63-11ea-19d9-51bbbdf5e6f6
# ╠═15626ee0-fd68-11ea-2c83-2dc85ba373e7
# ╠═973991a0-fd68-11ea-27c8-61246753b4b3
# ╠═ce78547e-fd68-11ea-3d71-61ec6c934761
# ╠═e64e2ad0-fd68-11ea-040d-6907d34eca35
# ╠═ef543200-fd68-11ea-1c37-c158193f71be
# ╠═0bf32920-fd69-11ea-293c-4949b8fee869


text:
# 创建  周期性的执行命令
#	类似 cron 的效果
# 本文件的执行命令:
# 	nohup:输入输出重定向,忽略sigup信号,来创建进程
# >>> nohup julia nohup_order.jl > log.txt
import Dates.dayofweek
import Dates.Hour
import Dates.Minute
import Dates.now

function 时间对正(hour::Int,min::Int)
	整点=Hour(now())
	分钟=Minute(now())
	if (整点==hour)&&(分钟==min)
		return true  #到达指定时间
	else
		return false
	end
end
#=========================
	每周1-5执行: 股市信息采集
		上午,9:30
		下午,1:00
	110m的内存占用,得不偿失!
=======================#
function 股市信息采集()
	星期=dayofweek(now())
	if 0< 星期< 6
		while true
			if 时间对正(9,30)||时间对正(13,00)
				#run(`julia julia/股市/analyse.jl`)
				@info "执行;时间:$(now())"
			else
				@info "不执行;时间:$(now())"
			end
			sleep(60)
		end
	end
end
股市信息采集()

max_stars_repo_name: Jonson-Sun/C-plus-plus

max_issues_repo_name: Jonson-Sun/C-plus-plus

max_forks_repo_name: Jonson-Sun/C-plus-plus


text:
@testset "core" begin

@testset "Tape" begin
    # Simple tests for `Tape`.
    @test getindex(setindex!(Tape(5), "hi", 5), 5) == "hi"
    @test lastindex(Tape(50)) == 50
    @test eachindex(Tape(50)) == Base.OneTo(50)
    @test length(Tape()) == 0
    @test length(Tape(50)) == 50
    @test !isassigned(Tape(1), 1)
    @test !isassigned(Tape(26), 26)

    # Check that printing works as expected.
    @testset "printing" begin
        buffer = IOBuffer()
        show(buffer, Tape())
        @test String(take!(buffer)) == "Tape with 0 elements"

        buffer = IOBuffer()
        show(buffer, Tape(1))
        @test String(take!(buffer)) == "Tape with 1 element:\n  [1]: #undef"

        buffer = IOBuffer()
        show(buffer, Tape(2))
        @test String(take!(buffer)) == "Tape with 2 elements:\n  [1]: #undef\n  [2]: #undef"

        buffer = IOBuffer()
        tape_ = Tape(1)
        tape_[1] = 5
        show(buffer, tape_)
        @test String(take!(buffer)) == "Tape with 1 element:\n  [1]: 5"
    end

    # Check isassigned consistency.
    leaf = Leaf(Tape(), 5)
    @test isassigned(tape(tape(leaf)), pos(leaf)) == isassigned(tape(leaf), leaf)

    # Simple tests for `Leaf`.
    tp1, tp2 = Tape(), Tape(50)
    @test tape(Leaf(tp1, 5.0)) == tp1
    @test unbox(Leaf(tp1, 5.0)) == 5.0
    @test pos(Leaf(tp1, 5.0)) == 3
    @test pos(Leaf(tp2, 5)) == 51
    @test unbox(Leaf(tp2, 5)) == 5
    @test tape(Leaf(tp2, 5)) == tp2

    # Simple tests for `Branch`.
    foo_coeff = 10
    foo(x::Real) = foo_coeff * x
    foo(x::Node{T} where T<:Real) = Branch(foo, (x,), tape(x))
    Nabla.∇(::typeof(foo), ::Type{Arg{1}}, p, y, ȳ, x) = ȳ * foo_coeff
    function get_new_branch()
        leaf = Leaf(Tape(), 5)
        return foo(leaf)
    end
    @test isa(get_new_branch(), Branch)
    @test unbox(get_new_branch()) == 50
    @test getfield(get_new_branch(), :f) == foo
    @test length(getfield(get_new_branch(), :args)) == 1
    @test isa(getfield(get_new_branch(), :args)[1], Leaf)
    @test unbox(getfield(get_new_branch(), :args)[1]) == 5
    @test pos(getfield(get_new_branch(), :args)[1]) == 1
    @test pos(get_new_branch()) == 2
    @test tape(getfield(get_new_branch(), :args)[1]) != tape(get_new_branch())
    branch = get_new_branch()
    @test tape(branch) == tape(getfield(branch, :args)[1])

    # Simple test for `pos`.
    @test Nabla.pos(1) == -1
    @test Nabla.pos("hi") == -1
    @test Nabla.pos(Leaf(Tape(), 5)) == 1
    @test Nabla.pos(Leaf(Tape(50), 5.0)) == 51

    # Simple tests for `unbox`.
    @test Nabla.unbox(1) == 1
    @test Nabla.unbox("hi") == "hi"
    @test Nabla.unbox(Leaf(Tape(), 5)) == 5
    @test Nabla.unbox(Leaf(Tape(), 5.0)) == 5.0

    # Simple tests for reverse_tape.
    @test isa(Nabla.reverse_tape(Leaf(Tape(), 5), 5), Tape)
    @test length(Nabla.reverse_tape(Leaf(Tape(), 5), 5)) == 1
    @test Nabla.reverse_tape(Leaf(Tape(), 5), 4)[end] == 4
    @test isassigned(Nabla.reverse_tape(Leaf(Tape(), 5), 5), 1)
    @test isassigned(Nabla.reverse_tape(get_new_branch(), 5), 2)
    @test !isassigned(Nabla.reverse_tape(get_new_branch(), 5), 1)

    # Simple tests for `propagate`.
    @test Nabla.propagate(Leaf(Tape(), 5), Tape()) == nothing
    br = get_new_branch()
    rvs_tape = Nabla.reverse_tape(br, 5)
    @test Nabla.propagate(br, rvs_tape) == nothing
    @test rvs_tape[br] == 5
    @test rvs_tape[br.args[1]] == foo_coeff * 5

    # Simple integration tests for ∇.
    br = get_new_branch()
    @test ∇(br, 3)[br] == 3
    @test ∇(br, 2)[br.args[1]] == 2 * foo_coeff

end # testset Tape

@testset "Check that functions involving `isapprox` can be differentiated" begin
    @testset "Test Case 1" begin
        f(x) = x ≈ 5.0 ? 1.0 : 3.0 * x
        g(x) = 5.0 * x
        h(x) = g(x) ≈ 25.0 ? x : f(x) + g(x)
        ∇f = ∇(f)
        ∇h = ∇(h)
        @test ∇f(5.0) == (0.0,)
        @test ∇f(6.0) == (3.0,)
        @test ∇h(5.0) == (1.0,)
        @test ∇h(6.0) == (8.0,)
    end

    @testset "Test Case 2" begin
        f(x) = x ≈ [5.0] ? 1.0 : 3.0 * sum(x)
        ∇f = ∇(f)
        @test ∇f([5.0]) == ([0.0],)
        @test ∇f([6.0]) == ([3.0],)
        f(x, y) = x ≈ y ? 2y : 3x
        ∇f = ∇(f)
        @test ∇f(5.0, 5.0) == (0.0, 2.0)
        @test ∇f(6.0, 5.0) == (3.0, 0.0)
    end
end

@testset "Check that functions with extra, unused variables can be differentiated" begin
    f(a,b,c,d) = a*c
    ∇f = ∇(f)
    g(a,b) = 12
    ∇g = ∇(g)

    @test ∇f(1,2,3,4) == (3, 0, 1, 0)
    @test ∇f(1,[2.0],3,4.0) == (3, [0.0], 1, 0.0)
    @test ∇g(1,2) == (0,0)
end

@testset "Check that functions with `zero` and `one` can be differentiated" begin
    f(a) = zero(a)
    g(a) = one(a)
    h(a) = zero(3 * a) + one(4 * a)
    ∇f = ∇(f)
    ∇g = ∇(g)
    ∇h = ∇(h)

    @test ∇f(1) == (0,)
    @test ∇f([1]) == ([0],)
    @test ∇g(4) == (0,)
    @test ∇h(8) == (0,)
end

@testset "Check that the convenience implementation of ∇ works as intended." begin
    f(x, y) = 2x + y
    ∇f = ∇(f)
    ∇f_out = ∇(f; get_output=true)

    @test_throws MethodError ∇f(randn(5), randn(5))
    x, y = randn(), randn()
    ∇z = ∇(f(Leaf.(Tape(), (x, y))...))
    @test ∇f(x, y) == (∇z[1], ∇z[2])
    z, (∇x, ∇y) = ∇f_out(x, y)
    @test unbox(z) == f(x, y)
    @test (∇x, ∇y) == ∇f(x, y)
end

@testset "get_output" begin
    y = ∇(-, get_output=true)(2)
    @test y isa Tuple{Branch{Int}, Tuple{Real}}
    @test last(y) == (-1,)
    @test ∇(unbox, get_output=true)(2) == (2, (0,))
end

@testset "Tests for zero'd and one'd containers." begin
    import Nabla: zerod_container, oned_container
    @test zerod_container(1.0) == 0.0
    @test zerod_container(1) == 0
    @test oned_container(1.0) == 1.0
    @test oned_container(5) == 1

    @test zerod_container(randn(5)) == zeros(5)
    @test oned_container(randn(5)) == ones(5)
    @test zerod_container(randn(5, 4, 3, 2, 1)) == zeros(5, 4, 3, 2, 1)
    @test oned_container(randn(5, 4, 3, 2, 1)) == ones(5, 4, 3, 2, 1)

    @test zerod_container((1.0, 1)) == (0.0, 0)
    @test oned_container((0, 0.0)) == (1, 1.0)
    @test zerod_container((randn(), randn(5))) == (0.0, zeros(5))
    @test oned_container((randn(5), randn())) == (ones(5), 1.0)
    @test zerod_container((1.0, (randn(5), randn(5)))) == (0.0, (zeros(5), zeros(5)))
    @test oned_container((randn(), (randn(5), randn(5)))) == (1.0, (ones(5), ones(5)))

    @test zerod_container([[1.0], [1.0]]) == [[0.0], [0.0]]
    @test oned_container([[5.0], [4.0]]) == [[1.0], [1.0]]

    @test zerod_container(Dict("a"=>5.0, "b"=>randn(3))) == Dict("a"=>0.0, "b"=>zeros(3))
    @test oned_container(Dict("a"=>5.0, "b"=>randn(3))) == Dict("a"=>1.0, "b"=>ones(3))

    @test ref_equal(zerod_container(Ref(4)), Ref(0))
    @test ref_equal(oned_container(Ref(4)), Ref(1))
end

# To ensure we end up using the fallback machinery for ∇(x̄, f, ...) we'll define a new
# function and setup for it to use in the testset below
quad(A::Matrix, B::Matrix) = B'A*B
@explicit_intercepts quad Tuple{Matrix, Matrix}
Nabla.∇(::typeof(quad), ::Type{Arg{1}}, p, Y, Ȳ, A::Matrix, B::Matrix) = B*Ȳ*B'
Nabla.∇(::typeof(quad), ::Type{Arg{2}}, p, Y, Ȳ, A::Matrix, B::Matrix) = A*B*Ȳ' + A'B*Ȳ

@testset "Mutating values in the tape" begin
    rng = MersenneTwister(123456)
    n = 5
    A = Leaf(Tape(), randn(rng, n, n))
    B = randn(rng, n, n)
    Q = quad(A, B)
    QQ = quad(Q, B)
    rt = ∇(QQ, Matrix(1.0I, n, n))
    oldvals = map(deepcopy∘unbox, getfield(rt, :tape))
    Nabla.propagate(Q, rt)  # This triggers a mutating addition
    newvals = map(unbox, getfield(rt, :tape))
    @test !(oldvals[1] ≈ newvals[1])
    @test oldvals[2:end] ≈ newvals[2:end]
end

end # testset "core"


text:
# define our `AbstractInterpreter` (not fancier one, which set ups another abstract layer and allows yet more composablility)
# and satisfy its interface requirements

"""
    XInterpreter <: AbstractInterpreter

Basically, this `AbstractInterpreter` just traces the original `NativeInterpreter`'s
abstract interpretation methods, but with additional tweaks to enable extensible inspections:
- `XInterpreter` setups its own code cache (both global and local), which can be referred later
- Entry points of the `XInterpreter`'s compilation pipeline can easily configure its parameters,
  and each different configuration will yield a different code cache.
"""
mutable struct XInterpreter <: AbstractInterpreter
    # NativeInterpreter
    native::NativeInterpreter
    optimize::Bool
    compress::Bool
    discard_trees::Bool
    verbose_stmt_info::Bool

    # customized cache
    cache_key::UInt

    # debug
    depth::Int

    function XInterpreter(world             = get_world_counter();
                          inf_params        = nothing,
                          opt_params        = nothing,
                          optimize          = true,
                          compress          = false,
                          discard_trees     = false,
                          verbose_stmt_info = false,
                          xconfigs...)
        isnothing(inf_params) && (inf_params = XInferenceParams(;    xconfigs...))
        isnothing(opt_params) && (opt_params = XOptimizationParams(; xconfigs...))

        # generate cache key, and initialize the cache partition
        h = @static UInt === UInt64 ? 0xacd4508fa540c2c9 : 0x43e699d6
        h = hash(inf_params, h)
        h = hash(opt_params, h)
        h = hash(optimize, h)
        h = hash(compress, h)
        cache_key = hash(discard_trees, h)
        __init_cache!(cache_key)

        native = NativeInterpreter(world; inf_params, opt_params)
        return new(native,
                   optimize,
                   compress,
                   discard_trees,
                   verbose_stmt_info,
                   cache_key,
                   0,
                   )
    end
end

# define these functions just to make them able to accept other JET configrations
XInferenceParams(; ipo_constant_propagation::Bool        = true,
                   aggressive_constant_propagation::Bool = false,
                   unoptimize_throw_blocks::Bool         = true,
                   max_methods::Int                      = 3,
                   union_splitting::Int                  = 4,
                   apply_union_enum::Int                 = 8,
                   tupletype_depth::Int                  = 3,
                   tuple_splat::Int                      = 32,
                   __xconfigs...) =
    return InferenceParams(; ipo_constant_propagation,
                             aggressive_constant_propagation,
                             unoptimize_throw_blocks,
                             max_methods,
                             union_splitting,
                             apply_union_enum,
                             tupletype_depth,
                             tuple_splat,
                             )
XOptimizationParams(; inlining::Bool                = inlining_enabled(),
                      inline_cost_threshold::Int    = 100,
                      inline_nonleaf_penalty::Int   = 1000,
                      inline_tupleret_bonus::Int    = 250,
                      inline_error_path_cost::Int   = 20,
                      max_methods::Int              = 3,
                      tuple_splat::Int              = 32,
                      union_splitting::Int          = 4,
                      unoptimize_throw_blocks::Bool = true,
                      __xconfigs...) =
    return OptimizationParams(; inlining,
                                inline_cost_threshold,
                                inline_nonleaf_penalty,
                                inline_tupleret_bonus,
                                inline_error_path_cost,
                                max_methods,
                                tuple_splat,
                                union_splitting,
                                unoptimize_throw_blocks,
                                )
# # assert here that they create same objects as the original constructors
@assert XInferenceParams()    == InferenceParams()
@assert XOptimizationParams() == OptimizationParams()

CC.InferenceParams(interp::XInterpreter)    = InferenceParams(interp.native)
CC.OptimizationParams(interp::XInterpreter) = OptimizationParams(interp.native)
CC.get_world_counter(interp::XInterpreter)  = get_world_counter(interp.native)

CC.lock_mi_inference(::XInterpreter,   ::MethodInstance) = nothing
CC.unlock_mi_inference(::XInterpreter, ::MethodInstance) = nothing

CC.add_remark!(interp::XInterpreter, sv, s) = add_remark!(interp.native, sv, s)

CC.may_optimize(interp::XInterpreter)      = interp.optimize
CC.may_compress(interp::XInterpreter)      = interp.compress
CC.may_discard_trees(interp::XInterpreter) = interp.discard_trees
CC.verbose_stmt_info(interp::XInterpreter) = interp.verbose_stmt_info


text:
# Rational Quadratic ARD Covariance Function

@doc """
# Description
Constructor for the ARD Rational Quadratic kernel (covariance)

k(x,x') = σ²(1+(x-x')ᵀL⁻²(x-x')/2α)^{-α}, where L = diag(ℓ₁,ℓ₂,...)
# Arguments:
* `ll::Vector{Float64}`: Log of length scale ℓ
* `lσ::Float64`        : Log of the signal standard deviation σ
* `lα::Float64`        : Log of shape parameter α
""" ->
type RQArd <: StationaryARD
    ℓ2::Vector{Float64}      # Length scale
    σ2::Float64              # Signal std
    α::Float64               # Shape parameter
    RQArd(ll::Vector{Float64}, lσ::Float64, lα::Float64) = new(exp(2.0*ll), exp(2.0*lσ), exp(lα))
end

function set_params!(rq::RQArd, hyp::Vector{Float64})
    length(hyp) == num_params(rq) || throw(ArgumentError("RQArd kernel has $(num_params(rq_ard)) parameters"))
    d = length(rq.ℓ2)
    rq.ℓ2 = exp(2.0*hyp[1:d])
    rq.σ2 = exp(2.0*hyp[d+1])
    rq.α = exp(hyp[d+2])
end

get_params(rq::RQArd) = [log(rq.ℓ2)/2.0; log(rq.σ2)/2.0; log(rq.α)]
get_param_names(rq::RQArd) = [get_param_names(rq.ℓ2, :ll); :lσ; :lα]
num_params(rq::RQArd) = length(rq.ℓ2) + 2

metric(rq::RQArd) = WeightedSqEuclidean(1.0./(rq.ℓ2))
cov(rq::RQArd,r::Float64) = rq.σ2*(1+0.5*r/rq.α)^(-rq.α)

function grad_kern(rq::RQArd, x::Vector{Float64}, y::Vector{Float64})
    wdiff = ((x-y).^2)./rq.ℓ2
    r  = sum(wdiff)
    part  = (1+0.5*r/rq.α)

    g1   = rq.σ2*wdiff*part^(-rq.α-1)
    g2   = 2.0*rq.σ2*part^(-rq.α)
    g3   = rq.σ2*part^(-rq.α)*(0.5*r/part-rq.α*log(part))
    return [g1; g2; g3]
end

function grad_stack!(stack::AbstractArray, rq::RQArd, X::Matrix{Float64}, data::StationaryARDData)
    d = size(X,1)
    R = distance(rq, X, data)
    part  = (1+0.5*R/rq.α)

    stack[:,:,d+2] = cov(rq, X)
    ck = view(stack, :, :, d+2)
    part2 = ck./part

    broadcast!(*, view(stack, :, :, 1:d), rq.σ2, data.dist_stack, reshape(1.0./rq.ℓ2, (1,1,d)), part.^(-rq.α - 1))

    stack[:,:, d+1] = 2.0 * ck
    stack[:,:, d+2] = ck.*(0.5*R./part-rq.α*log(part))
    return stack
end


text:
#
#	apply Heff
#
function applyHeff(x::AbstractArray{ComplexF64}, localMPO::MPOsparseT, L::Array{<:Any,1}, R::Array{<:Any,1})

	res = zeros(ComplexF64, size(L[1])[1], size(x)[2], size(R[1])[1])
	tmp = zeros(ComplexF64, size(L[1])[2], size(R[1])[1])

	@inbounds for op in eachoperation(localMPO) # potentially possible to run in parallel

		@inbounds r = @view R[op[3][2]][:, :]
		@inbounds l = @view L[op[3][1]][:, :]
		c = ComplexF64(op[1])

		m = @view x[:, op[2][1], :]
		new_m = @view res[:, op[2][2], :]

		LinearAlgebra.BLAS.gemm!('N','T', c, m, r, ComplexF64(0.0), tmp)
		LinearAlgebra.BLAS.gemm!('N', 'N', ComplexF64(1.0), l, tmp, ComplexF64(1.0), new_m)
	 end

	return res #(α_bra, d_bra, β_bra)
end

function applyHeff2(x::AbstractArray{<:Number}, localMPO1::MPO, localMPO2::MPO, L, R)

	#=
	res = zeros(ComplexF64, size(L[1])[1], size(x)[2:3]..., size(R[1])[1])
	tmp = zeros(ComplexF64, size(L[1])[2], size(R[1])[1])

	@inbounds for op2 in eachoperation(localMPO2) # potentially possible to run in parallel
		@inbounds for op in eachoperation(localMPO1) # potentially possible to run in parallel

			if op2[3][1] == op[3][2]
				@inbounds r = @view R[op2[3][2]][:, :]
				@inbounds l = @view L[op[3][1]][:, :]
				@inbounds m = @view x[:,op[2][1],op2[2][1],:]
				@inbounds new_m = @view res[:,op[2][2],op2[2][2],:]
				c1 = ComplexF64(op2[1])
				c2 = ComplexF64(op[1])

				LinearAlgebra.BLAS.gemm!('N','T', c1*c2, m, r, ComplexF64(0.0), tmp)
				LinearAlgebra.BLAS.gemm!('N', 'N', ComplexF64(1.0), l, tmp, ComplexF64(1.0), new_m)

			else
				continue
			end
		end

	end
	=#
	res = zeros(ComplexF64, size(x))

	mpo1 = get_MPOTensor(localMPO1)
	mpo2 = get_MPOTensor(localMPO2)

	@tensor bondOp[α, d1, d2, d1', d2', β] := mpo1[d1, α, γ, d2]*mpo2[d1', γ, β, d2']

	for (i, l) in enumerate(L)

		@tensor tmp[α, d1, d2, β] := l[α, γ]*x[γ, d1, d2, β]

		for (j, r) in enumerate(R)
			op = @view bondOp[i,:,:,:,:,j]
			@tensor tmp2[α, d1, d2, β] := tmp[α, d1, d2, γ]*r[β, γ]

			@tensor tmp2[α, d2, d2', β] := op[d1, d2, d1', d2']*tmp2[α, d1, d1', β]

			res+=tmp2
		end
	end

	return res

end

function applyHCeff(x::AbstractArray{<:Any}, L::Array{<:Any,1}, R::Array{<:Any,1})
	res = zeros(ComplexF64, size(L[1])[1], size(R[1])[1])
	tmp = similar(res)
	for diag_iter = 1:length(L)

		@inbounds l = @view L[diag_iter][:,:]
		@inbounds r = @view R[diag_iter][:,:]

		LinearAlgebra.BLAS.gemm!('N', 'N', ComplexF64(1.0), l, x, ComplexF64(0.0), tmp)
		LinearAlgebra.BLAS.gemm!('N', 'T', ComplexF64(1.0), tmp, r, ComplexF64(1.0), res)

 	end

	return res
end

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS


text:
using YaoDensityMatrixRegisters
using Test

@testset "YaoDensityMatrixRegisters.jl" begin
    # Write your tests here.
end


text:
#plot solution for 2d
## Transform to cartesian grid
@unpack ωx,ωy,ωz,Γ̄,M̄,g,t0,μ,ti,tf,Nt,t,dt = siminfo
R = sqrt(2μ/ωx^2)
xMax=1.5R
Nx = 500
x = collect(linspace(-xMax,xMax,Nx))
Tx = eigmat(M,x)

#Plot
f=figure(figsize=(9,2))
@manipulate for i=1:length(t) withfig(f,clear=true) do
        ψ = Tx*sol[i];
        θ = unwrap(angle(ψ));
        subplot(1,2,1);plot(x/R,g*abs.(ψ).^2,x/R,ones(x)*μ,"g:")
        ylim(-1,1.2*μ);xlim(-xMax/R,xMax/R)
        xlabel(L"x/R");ylabel(L"g|ψ|^2")
        subplot(1,2,2);plot(x/R,θ)
        xlim(-xMax/R,xMax/R);#ylim(-pi*1.1,pi*1.1);
        xlabel(L"x/R");ylabel(L"angle(ψ)")
    end
end


text:
using NewtonsMethod
using Test

@testset "NewtonsMethod.jl" begin

    # Define some functions to be tested

    # a. root = 1.0
    f(x) = (x - 1)^3
    df(x) = 3*(x - 1)^2

    # b. root = 0.6
    g(x) = (5*x - 3)^2
    dg(x) = 10*(5*x - 3)

    # c. no root (root = nothing)
    h(x) = 2 + x^2
    dh(x) = 2*x

    # 1. Test function f(x)
    @test NewtonsMethod.newtonroot(f, df)[1] ≈ 1.0
    @test NewtonsMethod.newtonroot(f)[1] ≈ 1.0

    # 2. Test function g(x)
    @test NewtonsMethod.newtonroot(g, dg)[1] ≈ 0.6
    @test NewtonsMethod.newtonroot(g)[1] ≈ 0.6

    # 3. Test using BigFloat precision instead of Float
    @test BigFloat(NewtonsMethod.newtonroot(f, df)[1]) ≈ 1.0
    @test BigFloat(NewtonsMethod.newtonroot(f)[1]) ≈ 1.0

    x0 = BigFloat(0.3)
    @test NewtonsMethod.newtonroot(f, df; iv = x0)[1] ≈ 1.0
    @test NewtonsMethod.newtonroot(f; iv = x0)[1] ≈ 1.0

    # 4. Test non-convergence of h(x)
    @test NewtonsMethod.newtonroot(h, dh)[1] === nothing
    @test NewtonsMethod.newtonroot(h)[1] === nothing

    # 5. Test maxiter
    # From PS1, root of f(x) found in ≈ 20 iterations
    @test NewtonsMethod.newtonroot(f, df; maxiter = 1)[1] === nothing
    @test NewtonsMethod.newtonroot(f; maxiter = 1)[1] === nothing
    @test NewtonsMethod.newtonroot(f, df; maxiter = 200)[1] ≈ 1.0
    @test NewtonsMethod.newtonroot(f; maxiter = 200)[1] ≈ 1.0

    # 6. Test tolerance
    @test !(NewtonsMethod.newtonroot(f, df; tol = 0.1)[1] ≈ 1.0)
    @test !(NewtonsMethod.newtonroot(f; tol = 0.1)[1] ≈ 1.0)

    # 7. Test initial value
    @test NewtonsMethod.newtonroot(f, df; iv = 0.5)[1] ≈ 1.0
    @test NewtonsMethod.newtonroot(f; iv = 0.5)[1] ≈ 1.0

    # 8. Break a test
    # @test NewtonsMethod.newtonroot(f, df)[1] ≈ 0.2
    # @test NewtonsMethod.newtonroot(f)[1] ≈ 0.2

end

max_stars_repo_name: loforteg/NewtonsMethod

max_issues_repo_name: loforteg/NewtonsMethod


text:
using RecurrenceAnalysis
using Test

include("dynamicalsystems.jl")
include("smallmatrix.jl")