text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libarpack"], :libarpack),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaLinearAlgebra/ArpackBuilder/releases/download/v3.5.0-3"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.aarch64-linux-gnu-gcc4.tar.gz", "d19fd4f1aa62c189f93462cbe25d6765c2ad42f8762d6c374e9e4e916f573aba"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.aarch64-linux-gnu-gcc7.tar.gz", "66242daa33cfeeb2085ffb1c374e93642bca3f2a5dc61f864d214c041441d25d"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.aarch64-linux-gnu-gcc8.tar.gz", "33245806067e5f96f4ad0b5070b52296e0f92b474fbfc9e4c0790253013ff127"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.aarch64-linux-musl-gcc4.tar.gz", "f8858dd20c7fd83dca71616be8adf1b2e3e86a7ece7e3ec52e05867da2290364"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.aarch64-linux-musl-gcc7.tar.gz", "a23b510a9fe2a78469c8d942f498ff1729bd9d7da9ceafda8318acf6a0057da0"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.aarch64-linux-musl-gcc8.tar.gz", "75f08c85072c4ae6d2d004efb99a6f9202e139746b120fa59e9c82fcdc207427"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.arm-linux-gnueabihf-gcc4.tar.gz", "d0601f37d571dff33f810aa3c2ba724c4cd163c17cbc57eb001057992cf01973"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.arm-linux-gnueabihf-gcc7.tar.gz", "c366fe80bdb3f705ebb5ddbd25267376e34085633931441ff50e0f0793c53dd5"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.arm-linux-gnueabihf-gcc8.tar.gz", "ea83ae4112c39cfec4fc317520827d0a89c61dfe01d1ecbf21b7349aa135035a"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.arm-linux-musleabihf-gcc4.tar.gz", "164150c704b4db50662987a05e23d0581e72157ca56cc9389a68ebfd242baef9"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.arm-linux-musleabihf-gcc7.tar.gz", "5b83bf6dac70d9ad492d69a4ca08a6b2efa0e80ddd59b4c44ec0db6169a01384"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.arm-linux-musleabihf-gcc8.tar.gz", "398d5b42ae1a631264f680c7133e27ff17396e5a7cf9a4830b54e3f88b06052f"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-linux-gnu-gcc4.tar.gz", "624c7b4e9ff830c8ebd88ead20250f02b9a2cbcd78b83b2e350d1d53abb3bdc3"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-linux-gnu-gcc7.tar.gz", "be5022190cb3365e5e3e43354c13c6adb245f4529d18b97b4efd6dea5fdb36e8"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-linux-gnu-gcc8.tar.gz", "b110c5ccb564f0617efb57cefbd46459a4830fbaa390f225167561b73b694081"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-linux-musl-gcc4.tar.gz", "8835f810f27661813e1fdb0fa02540215b02b5a6205a45dad6d02d00de82ca4c"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-linux-musl-gcc7.tar.gz", "9e3d16c23eb539f2345d5ebc7dc5df2a8a3236ae02492d13de50e3759cbbb6be"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-linux-musl-gcc8.tar.gz", "4db9a56c70c4cd5d55101c13250e800476b9cd6160a53249821cb00618b00b34"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-w64-mingw32-gcc4.tar.gz", "0a6b238b843491fa2dd783ddcc440d82f4d295ab33f91ce1a0c80835fa7805e2"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-w64-mingw32-gcc7.tar.gz", "17471e4289c7393e499ca3066e32f1a8776d1aad1fa994c9ec4398275be6e3cd"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.i686-w64-mingw32-gcc8.tar.gz", "c2da0cc187bdca8edee8bc3528c29f80c8a70eef165b449bcefe62b9affe51c5"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.powerpc64le-linux-gnu-gcc4.tar.gz", "dbb601c9c34d27fdbd74829874620e3dcf325bb05a3584d1be22b45bae3c1376"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.powerpc64le-linux-gnu-gcc7.tar.gz", "762601ab49f9d5ca5702968f167ce1a04b14ba8dc02d6a61bf7650c006c5bb42"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.powerpc64le-linux-gnu-gcc8.tar.gz", "efc2337f3a911377cf018d1827d9cd65f2314e1a81e000f92cb21b8b7fd6aa57"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-apple-darwin14-gcc4.tar.gz", "7ef977fd7214bbb96b57eb852da6e9227a89ea168403a23d0feb1602e77ade43"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-apple-darwin14-gcc7.tar.gz", "282bc463663328dd7f613b43b84459a7b87f59e61ceff780e8a606026613b67f"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-apple-darwin14-gcc8.tar.gz", "78295ca3721b22db9e76760676a7157c13ec3125e276c7da571a0d5fc5858929"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-linux-gnu-gcc4.tar.gz", "137a1617ff1c9bb64080fa0a36103294c2670568a00ae6a953307fb0625cdb46"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-linux-gnu-gcc7.tar.gz", "35e2f96333bfde0fb1749bf2f7d9a84fa8dfac6119d60a84a6b274d42dd85037"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-linux-gnu-gcc8.tar.gz", "11ae0abf2e42fa19364391bd9994c321aceb6604b204a62b062766f68ff7ef60"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-linux-musl-gcc4.tar.gz", "0a01691c184c6251521a54d24d8095210fa79d6ea4cc5d6777e4b45d1d4390e7"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-linux-musl-gcc7.tar.gz", "a40c92a2799b42605e0ccb17806b66bb1e1c848934a6d7fb7498df4aeaa29d9d"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-linux-musl-gcc8.tar.gz", "235e3816d930ef1552456586cb93cff798ee240f05561a42ac9e79eb0dc84a9a"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-unknown-freebsd11.1-gcc4.tar.gz", "8b0f3e6219c5fdcd58a5b0c7bfac6141b613c5846f0536e6d621b312d88204ee"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-unknown-freebsd11.1-gcc7.tar.gz", "f69db95295c172cd4281339f9727f1abc4491171ef0f08632000ad87f9d90cf3"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-unknown-freebsd11.1-gcc8.tar.gz", "f8ea43eb0624f77ab18582cea9e21dbab6981320b14625eaf461b152be4ea1b9"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-w64-mingw32-gcc4.tar.gz", "017832f856354f156910bfd9db9c04671863d1ea1c761435917fff3e5863c0b8"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-w64-mingw32-gcc7.tar.gz", "28beddf5a3e178012f910e406897d2abdf70ff3373c4b946ac39eae61070d042"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/Arpack.v3.5.0-3.x86_64-w64-mingw32-gcc8.tar.gz", "1fc0035e5bd62448243f8ed09041dee6a2480ea23e648c47b9d56d1f467a0ca0"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
#===============================================================================

            The main macro of this package `@diffusion_process`
            used for defining diffusion processes

===============================================================================#
# define some names for the document
_DIMENSIONS = [:dimensions, :dim, :dims, :dimension]
_PARAMETERS = [:parameters, :param, :params]
_CONSTPARAMETERS = [
    :const_parameters, :const_param, :const_params,
    :constparameters, :constparam, :constparams,
    :constant_parameters, :constant_param, :constant_params,
    :constantparameters, :constantparam, :constantparams
]
#_CONJUGATE = [:conjugate]
_ENDPOINTINFO = [:aux_info, :auxiliary_info, :end_points, :end_point_info]
_ENDPOINTEXTRA = [
    [:t0, :t_0],
    [:T],
    [:v0, :obs0, :v_0, :obs_0],
    [:vT, :obsT,  :v_T, :obs_T],
]
_EXTRA = [:additional, :extra]
_WIENER = [:wiener, :noise, :gaus, :gaussian]
_PROCESS = [:proc, :process, :state, :statespace]
_STATESPACE = [:statespace, :state, :domain, :domains]
_CONSTDIFF = [
    :constdiff,
    :constvola,
    :constdiffusivity,
    :constvolatility,
    :constσ,
    :constantdiff,
    :constantvola,
    :constantdiffusivity,
    :constantvolatility,
    :constantσ,
]
_DIAGONALDIFF = [
    :diagonaldiff,
    :diagonalvola,
    :diagonaldiffusivity,
    :diagonalvolatility,
    :diagonalσ,
    :diagdiff,
    :diagvola,
    :diagdiffusivity,
    :diagvolatility,
    :diagσ,
]
_SPARSEDIFF = [
    :sparsediff,
    :sparsevola,
    :sparsediffusivity,
    :sparsevolatility,
    :sparseσ,
]
_DIAGONALBMAT = [
    :diagonalbmat,
    :diagonalb,
    :diagonalbmatrix,
]
_SPARSEBMAT = [
    :sparsebmat,
    :sparseb,
    :sparsebmatrix,
]
_LINEAR = [:linear, :lineardiffusion]
_ELTYPE = [:eltype]
#_NUMNONHYPO = [:num_non_hypo, :numnonhypo]
#_PHI = [:phi, :ϕ, :φ]

_TEMP_DIMENSION_PROCESS = -1
_TEMP_DIMENSION_WIENER = -1

"""
    Base.lowercase(s::Symbol)

Lowercase all letters in a symbol
"""
Base.lowercase(s::Symbol) = Symbol(lowercase(string(s)))

"""
    _symbol_in(s::Symbol, symbols)

Check if symbol `s` is listed in a list of symbols `symbols`
"""
_symbol_in(s::Symbol, symbols) = lowercase(s) in symbols

"""
    _symbol_in(s::QuoteNode, symbols)

Check if quote of a symbol `s` is listed in a list of symbols `symbols`
"""
_symbol_in(s::QuoteNode, symbols) = _symbol_in(eval(s), symbols)

"""
    _symbol_in(::Any, ::Any)

Return `false` by default i.e. if the first argument is not a symbol nor its
quote
"""
_symbol_in(::Any, ::Any) = false

"""
    diffusion_process(name, ex::Expr, p...)

Defines a diffusion process according to a template described in the
documentation of the github repository:
https://github.com/mmider/DiffusionDefinition.jl
"""
macro diffusion_process(name, ex::Expr, p...)
    parse_process(name, MacroTools.striplines(ex), p)
end

"""
    @load_diffusion

Displays available choices of predefined diffusion processes that can be loaded
"""
macro load_diffusion()
    println("The following diffusions have been defined inside the package ",
            "DiffusionDefinition.jl and can be loaded by calling a macro ",
            "@load_diffusion name-of-a-diffusion:")
    for name in keys(_NAMES_TO_ACCEPTEDNAMES)
        println("- $name")
    end
    nothing
end

"""
    @load_diffusion name

Loads the predefined diffusion process.
"""
macro load_diffusion(name)
    found_me, importname = _load_diffusion(name)
    found_me && return Meta.parse("import DiffusionDefinition.$importname")
    nothing
end

function _load_diffusion(name)
    if _symbol_in(name, _ADMISSIBLENAMES)
        if typeof(name) <: QuoteNode
            name = eval(name)
        end
        homedir = joinpath(@__DIR__, "..", "examples")
        importname = _ACCEPTEDNAMES_TO_NAMES[lowercase(name)]
        filepath = _NAMES_TO_PATH[importname]
        path = joinpath(homedir, filepath)
        isfile(path) && include(path)
        !isfile(path) && println("Error, diffusion $name is supposed to be ",
                                "defined but it seems the file does not exist")
        return true, importname
    end
    println("Diffusion $name does not seem to be defined...")
    false, nothing
end

macro load_variable_diffusion(name, dim_proc, args...)
    global _TEMP_DIMENSION_PROCESS = eval(dim_proc)
    global _TEMP_DIMENSION_WIENER = (
        length(args) > 0 ?
        eval(args[1]) :
        _TEMP_DIMENSION_PROCESS
    )
    found_me, importname = _load_diffusion(name)
    found_me && return Meta.parse("import DiffusionDefinition.$importname")
    nothing
end

"""
    parse_process(name , ex::Expr, ::Any)

Parse a template defining a diffusion process, create a corresponding struct and
specified functions, evaluate them in the environment of a package and then
import the struct name to `Main` scope, in which the package has been imported
to.
"""
function parse_process(name , ex::Expr, ::Any)
    name, curly_name, template_args = process_name(name)
    p = (
        name = name,
        curly_name = curly_name,
        template_args = template_args,
        parameters = Vector{Tuple{Symbol,Union{Symbol,Expr},Symbol,Int64}}(
            undef,
            0
        ),
        dims = Dict{Symbol,Int64}(),
        extras = Dict{Symbol, Any}(),
        fns = Vector{Expr}(undef, 0),
        constp = Symbol[],
        nonconstp = Symbol[],
    )

    # parse lines defining parameters to get all parameter names first
    parse_lines!(ex, p, x->( x in [_PARAMETERS[1], _CONSTPARAMETERS[1]] ))
    # parse all other lines
    parse_lines!(ex, p, x-> !( x in [_PARAMETERS[1], _CONSTPARAMETERS[1]] ))
    fill_unspecified_with_defaults(p)

    abstract_type = prepare_abstract_type(
        p.extras[_LINEAR[1]] ? :LinearDiffusion : :DiffusionProcess,
        p.dims,
        p.extras[_ELTYPE[1]],
        p.extras[_STATESPACE[1]],
    )
    #struct_def, struct_body, struct_const = createstruct(abstract_type, p)
    struct_def = createstruct(abstract_type, p)
    #println(struct_def)
    #struct_def = Expr(
    #    :macrocall,
    #    :(Core.@doc),
    #    join(repr(struct_body), "\n\n\t", repr(struct_const)),
    #    struct_def,
    #)

    add_diff_function!(p.fns, p)
    add_parameter_names_function!(p.fns, p)
    add_parameter_names_function!(p.fns, p, :const_parameter_names, p.constp, identity)
    #add_parameter_names_function!(p.fns, p, :var_parameter_names, p.nonconstp, identity)
    #add_end_point_info_function!(p.fns, p)
    add_end_point_info_names_function!(p.fns, p)

    eval(struct_def)
    for fn in p.fns
        eval(fn)
    end
    println("A new struct `$curly_name` has been defined.")
    diffusion_message(Val(name))
    Meta.parse("import DiffusionDefinition.$name")
end

function diffusion_message(::Val{T}) where T
    #println("To learn more about how to define an instance of this struct")
    #println("please type in `?$T` and hit `ENTER`.")
end

"""
    process_name

Process the name of a struct by returning the pure name (without template
arguments), the name with template arguments and a list of template arguments.
"""
function process_name end
process_name(name::Symbol) = name, name, Any[]
process_name(name::Expr) = name.args[1], name, name.args[2:end]

"""
    parse_lines!(ex::Expr, p, condition)

Parse all lines of the expression `ex`, but process only those which satisfy
`condition`. `p` is a passed-around structure that accumulates processed
information.
"""
function parse_lines!(ex::Expr, p, condition)
    current_label = nothing
    for line in ex.args
        current_label = update_label(line, current_label)
        typeof(line) <: Union{Symbol, QuoteNode} && (continue)
        @assert line.head == :-->
        condition(current_label) && parse_line!(Val{current_label}(), line, p)
    end
end

"""
    update_label(line, current_label)

Update the label, which signifies what type of information a given line in a
template is supposed to be encoding.
"""
function update_label(line, current_label)
    _symbol_in(line, _DIMENSIONS) && return _DIMENSIONS[1]
    _symbol_in(line, _PARAMETERS) && return _PARAMETERS[1]
    _symbol_in(line, _CONSTPARAMETERS) && return _CONSTPARAMETERS[1]
    _symbol_in(line, _ENDPOINTINFO) && return _ENDPOINTINFO[1]
    #_symbol_in(line, _CONJUGATE) && return _CONJUGATE[1]
    _symbol_in(line, _EXTRA) && return _EXTRA[1]
    current_label
end

#------------------------------------------------------------------------------#
#
#                          For parsing: PARAMETERS
#
#------------------------------------------------------------------------------#
"""
    parse_line!(::Val{:parameters}, line, p)

Parse a line that defines parameters of the diffusion. The line must be in a
format:
    name --> parameter-description
"""
function parse_line!(
        ::K, line, p
    ) where K <: Union{Val{:parameters}, Val{:const_parameters}}
    param_name = line.args[1]

    if typeof(param_name) <: Symbol
        parse_param_single_name(line, p, K==Val{:const_parameters})
    else
        @assert param_name.head == :tuple
        parse_param_multi_names(line, p, K==Val{:const_parameters})
    end
end

"""
    parse_param_single_name(line, p)

Parse a line that defines parameters of the diffusion. The line must be in one
of the formats:
    parameter_name --> (number_of_parameters, datatype)
    parameter_name --> datatype
In the former case defines `number_of_parameters`-many parameters, with names
`parameter_name`i and of `datatype` type. In the latter case defines a
single parameters with name `parameter_name` and of type `datatype`.
"""
function parse_param_single_name(line, p, const_param=false)
    name_stem, disambig_idx, generic_name = get_name_stem(
        line.args[1], p.parameters
    )

    data_type = line.args[2]

    if _is_datatype(data_type, p)
        num_params = 1
    else
        @assert data_type.head == :tuple
        num_params, data_type = data_type.args[1], data_type.args[2]
        generic_name = true
    end

    for i in 1:num_params
        name = ( generic_name ? Symbol(name_stem, disambig_idx+i) : name_stem )
        append!(p.parameters, [(name, data_type, name_stem, disambig_idx+i)])
        push!( (const_param ? p.constp : p.nonconstp), name)
    end
end

"""
    _is_datatype(sym, p)

Utility function that checks whether `sym` is a datatype. It returns true if
`sym` is either an in-built datatype (say float, StaticArray etc.) or if it is a
template argument.
"""
function _is_datatype(sym, p)
    (sym in p.template_args) && return true
    typeof(sym) <: Expr && sym.head == :tuple && return false
    typeof(eval(sym)) <: DataType
end

"""
    get_name_stem(name_stem::Symbol, parameters)

Get the stem of a name for a parameters and then add a disambiguation index.
Underscore `_` used in place of name is defaulted to `p`.
"""
function get_name_stem(name_stem::Symbol, parameters)
    if name_stem == :_
        name_stem = :p
    end
    disambig_idx = highest_idx_used(name_stem, parameters)
    name_stem, disambig_idx, (name_stem == :p)
end

"""
    highest_idx_used(name_stem, params)

Find the highest disambiguation index that has been used thus far for a given
`name_stem`.
"""
function highest_idx_used(name_stem, params)
    idxs = [p[4] for p in params if name_stem == p[3]]
    length(idxs) == 0 && return 0
    maximum(idxs)
end

"""
    parse_param_multi_names(line, p)

Parse a line that defines parameters of the diffusion. The line must be in one
of the formats:
    (p_name1, p_name2, ...) --> (number_of_parameters, datatype)
    (p_name1, p_name2, ...) --> datatype
    (p_name1, p_name2, ...) --> (datatype1, datatype2, ...)
In the former two cases defines `number_of_parameters`-many parameters, with
names `p_name1`, `p_name2`, ... and of `datatype` type. In the last case the
datatypes differ from parameter to parameter.
"""
function parse_param_multi_names(line, p, const_param=false)
    names = line.args[1].args
    data_types = line.args[2]

    if _is_datatype(data_types, p)
        data_types = repeat([data_types], length(names))
    else
        @assert data_types.head == :tuple
        if _is_datatype(data_types.args[1], p)
            @assert length(data_types.args) == length(names)
            @assert all( map(x->_is_datatype(x, p), data_types.args) )
            data_types = data_types.args
        else
            @assert typeof(eval(data_types.args[1])) <: Number
            @assert length(data_types.args) == 2
            @assert eval(data_types.args[1]) == length(names)
            @assert _is_datatype(data_types.args[2], p)
            data_types = repeat([data_types.args[2]], length(names))
        end
    end

    current_param_names = [param[1] for param in p.parameters]
    for (n,d) in zip(names, data_types)
        @assert !(n in current_param_names)
        append!(p.parameters, [(n,d,n,1)])
        push!( (const_param ? p.constp : p.nonconstp), n)
    end
end

function add_parameter_names_function!(
        fns, p, f_name=:parameter_names, coll=p.parameters, λ=(x->x[1])
    )

    fn_def = Expr(
        :call,
        f_name,
        Expr(
            :(::),
            Expr(
                :curly,
                :Type,
                Expr(
                    :(<:),
                    p.name,
                )
            )
        )
    )
    _, used_names = organize_end_point_info(p)

    fn_body = Expr(
        :tuple,
        [QuoteNode(λ(x)) for x in coll]...,
        [QuoteNode(name) for name in used_names]...,
    )
    push!(fns, Expr(:(=), fn_def, fn_body))
end

#------------------------------------------------------------------------------#
#
#                      For parsing: END POINT INFORMATION
#
#------------------------------------------------------------------------------#
"""
    parse_line!(::Val{:aux_info}, line, p)

Parse a line that defines parameters of the diffusion. The line must be in a
format:
    name --> parameter-description
"""
function parse_line!(::Val{:aux_info}, line, p)
    name = line.args[1]
    for i in eachindex(_ENDPOINTEXTRA)
        _symbol_in(name, _ENDPOINTEXTRA[i]) && (name = _ENDPOINTEXTRA[i][1])
    end

    p.extras[name] = line.args[2]
end

function add_end_point_info_names_function!(fns, p)
    fn_def = Expr(
        :call,
        :end_point_info_names,
        Expr(
            :(::),
            Expr(
                :curly,
                :Type,
                Expr(
                    :(<:),
                    p.name,
                )
            )
        )
    )

    _, used_names = organize_end_point_info(p)

    fn_body = Expr(
        :tuple,
        [QuoteNode(name) for name in used_names]...,
    )

    push!(fns, Expr(:(=), fn_def, fn_body))
end

#=
#------------------------------------------------------------------------------#
#
#                       For parsing: CONJUGATE UPDATES
#
#------------------------------------------------------------------------------#

"""
    parse_line!(::Val{:conjugate}, line, p)

Parse a line that defines conjugate updates.
"""
function parse_line!(::Val{:conjugate}, line, p)
    if typeof(line.args[1]) <: Expr
        @assert line.args[1].head == :call
        function_name = line.args[1].args[1]
        fns = Vector{Expr}(undef, 0)
        if _symbol_in(function_name, _PHI)
            t_sym = line.args[1].args[2]
            x_sym = line.args[1].args[3]
            @assert line.args[2].head == :tuple
            num_funs = length(line.args[2].args)
            for expr in line.args[2].args
                add_phi_function!(expr, fns, t_sym, x_sym, p)
            end
        elseif function_name == :nonhypo
            add_nonhypo_function!(fns, line, p)
        end
        append!(p.fns, fns)
    elseif _symbol_in(line.args[1], _NUMNONHYPO)
        num = line.args[2]
        @assert typeof(num) <: Integer
        p.extras[_NUMNONHYPO[1]] = num
    end
end

"""
    add_nonhypo_function!(fns, line, p)

Add a definition of a function `nonhypo` that takes a current diffusion state
`x` and returns the elliptic coordinates.
"""
function add_nonhypo_function!(fns, line, p)
    fn_def = Expr(:call,
        :nonhypo,
        Expr(:(::),
            :P,
            p.name
        ),
        line.args[1].args[2]
    )
    push!(fns, Expr(:(=), fn_def, line.args[2]))
end

"""
    add_phi_function!(expr::Expr, fns::Vector{Expr}, t::Symbol, x::Symbol, p)

Construct a `phi` function used for conjugate updates.
"""
function add_phi_function!(
        expr::Expr,
        fns::Vector{Expr},
        t::Symbol,
        x::Symbol,
        p
    )
    idx = length(fns)
    name = p.name
    phi_def = Expr(:call,
        _PHI[1],
        Expr(:(::),
            Val{idx}
        ),
        :t,
        :x,
        Expr(:(::), :P, p.name)
    )

    @assert expr.head == :tuple
    @assert length(expr.args) == p.dims[_PROCESS[1]]
    phi_body = copy(expr)
    for i in 1:p.dims[_PROCESS[1]]
        cleanup_param_names!(phi_body.args[i], t, x, map(x->x[1], p.parameters))
    end

    phi = Expr(:(=), phi_def, phi_body)
    append!(fns, [phi])
end

"""
    cleanup_param_names!(expr, t::Symbol, x::Symbol, params)

A hepler function for cleaning up the names of variables in the definition of a
`phi` function.
"""
function cleanup_param_names!(expr, t::Symbol, x::Symbol, params)
    for (i,arg) in enumerate(expr.args)
        if typeof(arg) <: Symbol
            arg == t && (expr.args[i] = :t)
            arg == x && (expr.args[i] = :x)
            arg in params && (expr.args[i] = Expr(:., :P, QuoteNode(arg)))
        elseif typeof(arg) == Expr
            cleanup_param_names!(arg, t, x, params)
        end
    end
end

cleanup_param_names!(num::Number, t::Symbol, x::Symbol, params) = num
=#
#------------------------------------------------------------------------------#
#
#                       For parsing: ADDITIONAL INFO
#
#------------------------------------------------------------------------------#

"""
    parse_line!(::Val{:additional}, line, p)

Parse a line that defines additional information about a diffusion process.
"""
function parse_line!(::Val{:additional}, line, p)
    name = line.args[1]
    _symbol_in(name, _STATESPACE) && (name = _STATESPACE[1])
    _symbol_in(name, _CONSTDIFF) && (name = _CONSTDIFF[1])
    _symbol_in(name, _SPARSEDIFF) && (name = _SPARSEDIFF[1])
    _symbol_in(name, _DIAGONALDIFF) && (name = _DIAGONALDIFF[1])
    _symbol_in(name, _SPARSEBMAT) && (name = _SPARSEBMAT[1])
    _symbol_in(name, _DIAGONALBMAT) && (name = _DIAGONALBMAT[1])
    _symbol_in(name, _LINEAR) && (name = _LINEAR[1])
    _symbol_in(name, _ELTYPE) && (name = _ELTYPE[1])

    p.extras[name] = line.args[2]
end

"""
    add_diff_function!(fns, p)

Add a definition of a function `consdiff` that indicates if the diffusion
coefficient is constant
"""
function add_diff_function!(fns, p)
    to_add = [
        (:constdiff, _CONSTDIFF[1]),
        (:sparsediff, _SPARSEDIFF[1]),
        (:diagonaldiff, _DIAGONALDIFF[1]),
        (:sparseBmat, _SPARSEBMAT[1]),
        (:diagonalBmat, _DIAGONALBMAT[1])
    ]
    for (name, _NAME) in to_add
        fn_def = Expr(:call,
            name,
            Expr(:(::),
                :P,
                p.name
            )
        )
        push!(fns, Expr(:(=), fn_def, p.extras[_NAME]))
    end
end

#------------------------------------------------------------------------------#
#
#                   For parsing: DIFFUSION DIMENSION
#
#------------------------------------------------------------------------------#

"""
    parse_line!(::Val{:dimensions}, line, p)

Parse a line that defines the dimension of a diffusion process and the driving
Brownian motion.
"""
function parse_line!(::Val{:dimensions}, line, p)
    dim_of_what, dim = line.args[1], line.args[2]
    _symbol_in(dim_of_what, _WIENER) && (p.dims[_WIENER[1]] = eval(dim))
    _symbol_in(dim_of_what, _PROCESS) && (p.dims[_PROCESS[1]] = eval(dim))
end

#------------------------------------------------------------------------------#
#
#                            final routines
#
#------------------------------------------------------------------------------#

"""
    fill_unspecified_with_defaults(p)

Fill all unspecified variables with default values
"""
function fill_unspecified_with_defaults(p)
    fill_unspecified_with_defaults(Val{:dimensions}(), p)
    fill_unspecified_with_defaults(Val{:additional}(), p)
end

"""
    fill_unspecified_with_defaults(::Val{:dimensions}, p)

If unspecified, the dimension of the stochastic process and the driving Brownian
motion is set to 1.
"""
function fill_unspecified_with_defaults(::Val{:dimensions}, p)
    for key in [_WIENER[1], _PROCESS[1]]
        !haskey(p.dims, key) && (p.dims[key] = 1)
    end
end

"""
    fill_unspecified_with_defaults(::Val{:additional}, p)

If unspecified, there are no restriction on a state space, the volatility
coefficient is assumed constant, the diffusion is not linear and the datatype
of each coordinate is set to `Float64`.
"""
function fill_unspecified_with_defaults(::Val{:additional}, p)
    defaults = [
        (_STATESPACE[1], UnboundedStateSpace()),
        (_CONSTDIFF[1], false),
        (_SPARSEDIFF[1], false),
        (_DIAGONALDIFF[1], false),
        (_SPARSEBMAT[1], false),
        (_DIAGONALBMAT[1], false),
        (_LINEAR[1], false),
        (_ELTYPE[1], Float64),
        #(_NUMNONHYPO[1], (
        #    haskey(p.extras, _STATESPACE[1]) ?
        #    p.extras[_STATESPACE[1]] :
        #    1
        #)),
    ]
    for (key, default_val) in defaults
        !haskey(p.extras, key) && (p.extras[key] = default_val)
    end
end

"""
    prepare_abstract_type(stem, dims, data_type, state_restr)

Create a string defining a parent, abstract type from its `stem`, the dimensions
`dims` of the process and the driving Brownian motion, the datatype `data_type`
of each coordinate and the restrictions on the state space `state_restr`.
"""
function prepare_abstract_type(stem, dims, data_type, state_restr)
    Expr(:curly,
        stem,
        data_type,
        dims[_PROCESS[1]],
        dims[_WIENER[1]],
        state_restr,
    )
end

"""
    createstruct(abstract_type, name, params)

Create code that defines a struct defining a diffusion process.
"""
function createstruct(abstract_type, p)
    param_vec = map(p.parameters) do (name, data_type, _, _)
        Expr(:(::), name, data_type)
    end
    end_point_info_vec, end_point_info_used_names = organize_end_point_info(p)

    _new = (
        p.template_args == Any[] ?
        :new :
        Expr(:curly, :new, p.template_args...,)
    )
    #=
    xT_default = (
        length(end_point_info_vec)>0 ?
        [
            Expr(:kw,
                end_point_info_vec[end],
                Expr(
                    :call,
                    :zero,
                    end_point_info_vec[end].args[2],
                    p.dims[_PROCESS[1]],
                    Expr(
                        :call,
                        Expr(
                            :.,
                            :DiffusionDefinition,
                            :(:ismutable),
                        ),
                        end_point_info_vec[end].args[2],
                    )
                ),
            )
        ] :
        []
    )
    =#

    constructor_def = add_where_decorator(
        Expr(
            :call,
            p.name,
            param_vec...,
            end_point_info_vec...,
        ),
        p
    )

    named_constructor_def = add_where_decorator(
        Expr(
            :call,
            p.name,
            Expr(
                :parameters,
                param_vec...,
                end_point_info_vec...,
            )
        ),
        p
    )

    partial_clone_constructor_def = add_where_decorator(
        Expr(
            :call,
            p.name,
            Expr(
                :parameters,
                param_vec...,
            ),
            end_point_info_vec...,
        ),
        p
    )

    build_costructor(x)= Expr(
        :function, # constructor
        x,
        Expr(:block,
            Expr(:call,
                _new,
                map(x->x[1], p.parameters)...,
                end_point_info_used_names...,
            ),
        ), # body of the constructor
    ) # constructor

    struct_def = Expr(:struct,
        true, # immutable
        Expr(:<:,
            p.curly_name,
            abstract_type,
        ), # first line defining the struct
        Expr(:block,
            param_vec..., # parameters
            end_point_info_vec..., # extra info about end-points
            build_costructor(constructor_def),
            build_costructor(named_constructor_def),
            build_costructor(partial_clone_constructor_def),
        ), # body defining the struct
    )
    struct_def#, Expr(:block, param_vec..., end_point_info_vec...,), constructor_def
end

function add_where_decorator(expr, p)
    (
        p.template_args == Any[] ?
        expr :
        Expr(:where,
            expr,
            p.template_args...,
        )
    )
end

function organize_end_point_info(p)
    end_point_info_names = [:t0, :T, :v0, :vT]
    end_point_info_vec = Any[]
    end_point_info_used_names = Symbol[]
    #if any([haskey(p.extras, eio) for eio in end_point_info_names])
    for name in end_point_info_names
        haskey(p.extras, name) && (
            push!(end_point_info_vec, Expr(:(::), name, p.extras[name]));
            push!(end_point_info_used_names, name)
        )
    end
    #=
        # the space for the exact end-point must always be there due to blocking
        default_datatype = find_state_datatype(
            p.extras[_ELTYPE[1]], # eltype
            p.dims[_PROCESS[1]], # dim of the process
        )
        push!(
            end_point_info_vec,
            Expr(:(::), :xT, get(p.extras, :xT, default_datatype))
        )
        push!(end_point_info_used_names, :xT)
        =#
    #end
    end_point_info_vec, end_point_info_used_names
end

#=
function find_state_datatype(_eltype, dim_process)
    dim_process < 9 && return SVector{dim_process,eval(_eltype)}
    return Vector{eval(_eltype)}
end
=#

ext: jl


text:
using RemoteGtkREPL

module GtkREPLWorker

    using Sockets
    using RemoteGtkREPL

    gtkrepl_port = parse(Int,ARGS[1])
    global const  id = parse(Int,ARGS[2])
    port, server = RemoteGtkREPL.start_server()

    global const gtkrepl = connect(gtkrepl_port)

end

GtkREPLWorker.RemoteGtkREPL.remotecall_fetch(
    include_string, GtkREPLWorker.gtkrepl,
    Main, "client_start_cb($(GtkREPLWorker.port))"
)


text:
"""

Plotting of symbolic objects.

The `Plots` package provide a uniform interface to many of `Julia`'s
plotting packages. `SymPy` plugs into `Plots`' "recipes."

The basic goal is that when `Plots` provides an interface for function
objects, this package extends the interface to symbolic expressions.

In particular:

* `plot(ex::Sym, a, b; kwargs...)` will plot a function evaluating `ex` over [a,b]

Example. Here we use the default backend for `Plots` to make a plot:

```
using Plots
@vars x
plot(x^2 - 2x, 0, 4)
```

* `plot(exs::Vector{Sym}, a, b; kwargs...)` will plot the functions evaluating `exs` over [a,b]

Example:

```
@vars x
plot([sin(x), cos(x)], 0, 2pi)
```

* `plot(ex1, ex2, a, b; kwargs...)` will plot the two expressions in a parametric plot over the interval `[a,b]`.

Example:

```
@vars x
plot(sin(2x), cos(3x), 0, 4pi) ## also
```

For a few backends (those that support `:path3d`) a third symbolic
expression may be added to have a 3d parametric plot rendered:

```
plot(sin(x), cos(x), x, 0, 4pi) # helix in 3d
```

* `plot(xs, ys, expression)` will make a contour plot (for many backends).

```
@vars x y
plot(linspace(0,5), linspace(0,5), x*y)
```

* To plot the surface  `z=ex(x,y)` over a region we have `Plots.surface`. For example,

```
@vars x y
surface(-5:5, -5:5, 25 - x^2 - y^2)
```

* a vectorfield plot can (inefficiently but directly) be produced following this example:

```
function vfieldplot(fx, fy; xlim=(-5,5), ylim=(-5,5), n=7)
    xs = linspace(xlim..., n)
    ys = linspace(ylim..., n)

    us = vec([x for x in xs, y in ys])
    vs = vec([y for x in xs, y in ys])
    fxs = vec([fx(x,y) for x in xs, y in ys])
    fys = vec([fy(x,y) for x in xs, y in ys])

    quiver(us, vs, quiver=(fxs, fys))
end
fx = (x + y) / sqrt(x^2 + y^2)
fy = (x - y) / sqrt(x^2 + y^2)
vfieldplot(fx, fy)
```

----

Some graphics provided by `SymPy` are available if `PyPlot` is installed.

* `plot_parametric_surface(exs::Tuple, (uvar,a0,b0), (vvar,a1,b1);
  kwargs)` will make parametrically defined surface plots.

Plot the parametrically defined surface `[exs[1](u,v), exs[2](u,v), exs[3](u,v)]` over `[a0,a1] x
[b0,b1]`. The specification of the variables uses a tuple of the form
`(Sym, Real, Real)` following the style of SymPy in `integrate`, say,
where disambiguation of variable names is needed.

```
@vars theta, phi
r = 1
plot_parametric_surface((r*sin(theta)*sin(phi), r*sin(theta)*cos(phi), r*cos(theta)),
                        (theta, 0, pi), (phi, 0, pi/2))
```

(The SymPy name for this function is `plot3d_parametric_surface`, we have dropped the "`3d`" part.)

* `plot_implicit(equation, (xvar, x0, x1), (yvar, y0, y1))` will plot implicitly the equation.

```
@syms x y
plot_implicit(Eq(x^2+ y^2,3), (x, -2, 2), (y, -2, 2))  # draw a circle
```

"""
sympy_plotting = nothing
export sympy_plotting

## Recipes for hooking into Plots

using RecipesBase

##
@recipe f{T<:Sym}(::Type{T}, v::T) = lambdify(v)

## for vectors of expressions
@recipe f{S<:AbstractVector{Sym}}(::Type{S}, ss::S) = Function[lambdify(s) for s in ss]

## vectorfieldplot is an addon as we avoid quiver
## Not ready yet!
## @recipe function f(::Type{Val{:vectorfieldplot}}, fx, fy, xlim=(-5,5), ylim=(-5,5); n=15)

##     x₀, x₁ = xlim
##     y₀, y₁ = ylim
##     Δx = (x₁ - x₀) / (n-1)
##     Δy = (y₁ - y₀) / (n-1)

##     xs = x₀:Δx:x₁
##     ys = y₀:Δy:y₁

##     mx, my = 0.0, 0.0

##     ## two loops, first to identify scaling factor so
##     ## lines stay within box of width Δx by Δy
##     for x in xs, y in ys
##         mx = max(mx, abs(fx(x,y)))
##         my = max(my, abs(fy(x,y)))
##     end

##     ## we want all lines to be in the box, so we scale
##     λ = .95 *  min(Δx/mx, Δy/my)

##     XS = Float64[]
##     YS = Float64[]

##     for x in xs, y in ys
##         append!(XS, [x, x + λ * fx(x,y), NaN])
##         append!(YS, [y, y + λ * fy(x,y), NaN])
##     end
##     pop!(XS); pop!(YS)

##     seriestype := :path
##     xlims := xlim
##     ylims := ylim
##     legend := false

##     x := XS
##     y := YS

##     ()
## end

## function vectorfieldplot(fx::Sym, fy::Sym, xlim=(-5,5), ylim=(-5,5); n::Int=15, kwargs...)
##     xs = free_symbols([fx, fy])
##     vectorfieldplot(lambdify(fx, xs) , lambdify(fy, xs), xlim, ylim; n=n, kwargs...)
## end

## function vectorfieldplot!(fx::Sym, fy::Sym, xlim=(-5,5), ylim=(-5,5); n::Int=15, kwargs...)
##     xs = free_symbols([fx, fy])
##     vectorfieldplot!(lambdify(fx, xs), lambdify(fy, xs), xlim, ylim; n=n, kwargs...)
## end

## function vectorfieldplot(fx, fy, xlim=(-5,5), ylim=(-5,5); n::Int=15, kwargs...)
##     try
##         eval(Expr(:using, :Plots)) # XXX Super hacky way to do this
##     catch err
##         warn("Can't plot without the `Plots` package.")
##     end
##     plot(Val{:vectorfieldplot}, fx, fy, xlim, ylim; n=n, kwargs...)
## end

## function vectorfieldplot!(fx, fy, xlim=(-5,5), ylim=(-5,5); n::Int=15, kwargs...)
##     try
##         eval(Expr(:using, :Plots)) # XXX Super hacky way to do this
##     catch err
##         warn("Can't plot without the `Plots` package.")
##     end
##     plot!(Val{:vectorfieldplot}, fx, fy, xlim, ylim; n=n, kwargs...)
## end
## export vectorfieldplot, vectorfieldplot!

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

## These functions give acces to SymPy's plotting module. They will work if PyPlot is installed, but may otherwise cause an error

## surface plot xvar = Tuple(Sym, Real, Real)
##
"""

Render a parametrically defined surface plot.

Example:
```
@vars u, v
plot_parametric_surface((u*v,u-v,u+v), (u,0,1), (v,0,1))
```

This uses `PyPlot`, not `Plots` for now.
"""
function plot_parametric_surface(exs::(@compat Tuple{Sym,Sym,Sym}),
                                 xvar=(-5.0, 5.0),
                                 yvar=(-5.0, 5.0),
                                 args...;
                                 kwargs...)

    SymPy.call_sympy_fun(sympy[:plotting][:plot3d_parametric_surface], exs..., args...; kwargs...)

end
export plot_parametric_surface

"""
Plot an implicit equation

```
@syms x y
plot_implicit(Eq(x^2+ y^2,3), (x, -2, 2), (y, -2, 2))
```

"""
plot_implicit(ex, args...; kwargs...) = SymPy.call_sympy_fun(sympy[:plotting][:plot_implicit], ex, args...; kwargs...)
export plot_implicit


text:
# This file is a part of NeighborhoodApproximationIndex.jl

using SearchModels, Random, StatsBase
import SearchModels: combine, mutate
import SimilaritySearch: optimize!, MinRecall, ParetoRecall, ParetoRadius, ErrorFunction, setconfig!, runconfig, optimization_space
export optimize!, KnrOptSpace

@with_kw struct KnrOptSpace <: AbstractSolutionSpace
    ksearch = 1:3:21
    ksearch_scale = (s=1.1, p1=0.8, p2=0.8, lower=1, upper=128)
end

Base.hash(c::KnrOpt) = hash(c.ksearch)
Base.isequal(a::KnrOpt, b::KnrOpt) = a.ksearch == b.ksearch
Base.eltype(::KnrOptSpace) = KnrOpt
Base.rand(space::KnrOptSpace) = KnrOpt(rand(space.ksearch))

combine(a::KnrOpt, b::KnrOpt) = KnrOpt((a.ksearch + b.ksearch) ÷ 2)
mutate(sp::KnrOptSpace, a::KnrOpt, iter) = KnrOpt(SearchModels.scale(a.ksearch; sp.ksearch_scale...))

optimization_space(index::KnrIndex) = KnrOptSpace()

function setconfig!(c::KnrOpt, index::KnrIndex, perf)
    index.opt.ksearch = c.ksearch
end

function runconfig(c::KnrOpt, index::KnrIndex, q, res::KnnResult, pools)
    search(index, q, res; pools, ksearch=c.ksearch)
end

max_stars_repo_name: sadit/NeighborhoodApproximationIndex

max_issues_repo_name: sadit/NeighborhoodApproximationIndex

max_forks_repo_name: sadit/NeighborhoodApproximationIndex


text:
# if VERSION <= v"0.5.0-"
#     import Base: isprime
# end

if isdefined(Base, :isprime)
    import Base: isprime
else
    import Primes: isprime
end

# Gaussian primes
function isprime{T<:Integer}(z::Complex{T})
    (x,y) = (real(z),imag(z))
    x == 0 && isprime(y) && y % 4 == 3 && return true
    y == 0 && isprime(x) && x % 4 == 3 && return true
    isprime(abs2(z))
end

#Base.@vectorize_1arg Complex isprime


text:
module NNBasicFunction

  export hardlim , hardlims, purelin, satlin, satlins, logsig, tansig, poslin;

   function hardlim(x)
		if (x<0)
	     return 0
	    else
	    return  +1
	   end
	end

    function hardlims(x)
		if (x<0)
	     return -1
	    else
	    return +1
	   end
    end

    function purelin(x)
	   x
   	end

   function satlin(x)
	   if(x<0)
		   return 0
	   elseif(0<=x<=1)
		   return x
	   else
		   return 1
	   end
   end

	function satlins(x)
	   if(x<-1)
		   return -1
	   elseif(-1<=x<=1)
		   return x
	   else
		   return 1
	   end
   end

	function logsig(x)
		let top=1,bottom=1+ℯ^-x
			return top/bottom
		end
	end

	function tansig(x)
		let first= ℯ^x, last=ℯ^-x
			top= first-last
			bottom=first+last
		return top/bottom
		end
	end
	function poslin(x)
		if(x<0)
		  return 0
		else
		  return x
		end
	end

end

max_stars_repo_name: phpsmarter/plutonote

max_issues_repo_name: phpsmarter/plutonote

max_forks_repo_name: phpsmarter/plutonote


text:
module TypeSortedCollections

export
    TypeSortedCollection,
    num_types,
    eltypes,
    vectortypes

const TupleOfVectors = Tuple{Vararg{Vector{T} where T}}

struct TypeSortedCollection{D<:TupleOfVectors, N}
    data::D
    indices::NTuple{N, Vector{Int}}

    function TypeSortedCollection{D, N}(data::D, indices::NTuple{N, Vector{Int}}) where {D<:TupleOfVectors, N}
        fieldcount(D) == N || error()
        l = mapreduce(length, +, data, init=0)
        l == mapreduce(length, +, indices, init=0) || error()
        if N > 0
            allindices = Base.Iterators.flatten(indices)
            allunique(allindices) || error()
            isempty(allindices) || extrema(allindices) == (1, l) || error()
        end
        new{D, N}(data, indices)
    end

    function TypeSortedCollection{D, N}(indices::NTuple{N, Vector{Int}}) where {D<:TupleOfVectors, N}
        lengths = map(length, indices)
        data = ntuple(i -> D.parameters[i](undef, lengths[i]), N)
        TypeSortedCollection{D, N}(data, indices)
    end

    function TypeSortedCollection{D, N}() where {D<:TupleOfVectors, N}
        indices = ntuple(_ -> Int[], N)
        TypeSortedCollection{D, N}(indices)
    end

    TypeSortedCollection{D}() where {D<:TupleOfVectors} = TypeSortedCollection{D, length(D.parameters)}()
    TypeSortedCollection{D, N}(A) where {D<:TupleOfVectors, N} = append!(TypeSortedCollection{D, N}(), A)
    TypeSortedCollection{D}(A) where {D<:TupleOfVectors} = append!(TypeSortedCollection{D}(), A)
end

function TypeSortedCollection(data::D, indices::NTuple{N, Vector{Int}}) where {D<:TupleOfVectors, N}
    TypeSortedCollection{D, N}(data, indices)
end

function TypeSortedCollection(A, preserve_order::Bool = false)
    if preserve_order
        data = Vector[]
        indices = Vector{Vector{Int}}()
        for (i, x) in enumerate(A)
            T = typeof(x)
            if isempty(data) || T != eltype(last(data))
                push!(data, T[])
                push!(indices, Int[])
            end
            push!(last(data), x)
            push!(last(indices), i)
        end
        TypeSortedCollection(tuple(data...), tuple(indices...))
    else
        types = unique(typeof.(A))
        D = Tuple{[Vector{T} for T in types]...}
        TypeSortedCollection{D}(A)
    end
end

function TypeSortedCollection(A, indices::NTuple{N, Vector{Int}} where {N})
    @assert length(A) == mapreduce(length, +, indices, init=0)
    data = []
    for indicesvec in indices
        T = length(indicesvec) > 0 ? typeof(A[indicesvec[1]]) : Nothing
        Tdata = Vector{T}()
        sizehint!(Tdata, length(indicesvec))
        push!(data, Tdata)
        for i in indicesvec
            A[i]::T
            push!(Tdata, A[i])
        end
    end
    TypeSortedCollection(tuple(data...), indices)
end

@inline Base.eltype(A::TypeSortedCollection) = Union{map(eltype, A.data)...}

eltypes(::Type{TypeSortedCollection{D, N}}) where {D, N} = eltypes(D)
function eltypes(::Type{T}) where {T <: TupleOfVectors}
    Base.tuple_type_cons(eltype(Base.tuple_type_head(T)), eltypes(Base.tuple_type_tail(T)))
end
eltypes(::Type{Tuple{}}) = Tuple{}

function vectortypes(::Type{T}) where {T <: Tuple}
    Base.tuple_type_cons(Vector{Base.tuple_type_head(T)}, vectortypes(Base.tuple_type_tail(T)))
end
vectortypes(::Type{Tuple{}}) = Tuple{}

@generated function Base.push!(dest::TypeSortedCollection{D}, x::X) where {D, X}
    i = 0
    for j = 1 : length(D.parameters)
        Vector{X} == D.parameters[j] && (i = j; break)
    end
    i == 0 && return :(throw(ArgumentError("Destination cannot store arguments of type $(typeof(x)).")))
    quote
        Base.@_inline_meta
        index = length(dest) + 1
        push!(dest.data[$i], x)
        push!(dest.indices[$i], index)
        return dest
    end
end

function Base.append!(dest::TypeSortedCollection, A)
    # TODO: consider resizing first
    for x in A
        push!(dest, x)
    end
    dest
end

num_types(::Type{<:TypeSortedCollection{<:Any, N}}) where {N} = N
num_types(x::TypeSortedCollection) = num_types(typeof(x))

const TSCOrAbstractVector{N} = Union{<:TypeSortedCollection{<:Any, N}, AbstractVector}

Base.isempty(x::TypeSortedCollection) = all(isempty, x.data)
Base.empty!(x::TypeSortedCollection) = foreach(empty!, x.data)
@inline Base.length(x::TypeSortedCollection) = mapreduce(length, +, x.data, init=0)
indices(x::TypeSortedCollection) = x.indices

# Trick from StaticArrays:
@inline first_tsc(a1::TypeSortedCollection, as...) = a1
@inline first_tsc(a1, as...) = first_tsc(as...)

@inline first_tsc_type(a1::Type{<:TypeSortedCollection}, as::Type...) = a1
@inline first_tsc_type(a1::Type, as::Type...) = first_tsc_type(as...)

# inspired by Base.ith_all
@inline _getindex_all(::Val, j, vecindex) = ()
Base.@propagate_inbounds @inline _getindex_all(vali::Val{i}, j, vecindex, a1, as...) where {i} = (_getindex(vali, j, vecindex, a1), _getindex_all(vali, j, vecindex, as...)...)
@inline _getindex(::Val, j, vecindex, a) = a # for anything that's not an AbstractVector or TypeSortedCollection, don't index (for use in broadcast!)
@inline _getindex(::Val, j, vecindex, a::AbstractVector) = a[vecindex]
@inline _getindex(::Val, j, vecindex, a::Ref) = a[]
@inline _getindex(::Val{i}, j, vecindex, a::TypeSortedCollection) where {i} = a.data[i][j]
@inline _setindex!(::Val, j, vecindex, a::AbstractVector, val) = a[vecindex] = val
@inline _setindex!(::Val{i}, j, vecindex, a::TypeSortedCollection, val) where {i} = a.data[i][j] = val

@inline lengths_match(a1) = true
@inline lengths_match(a1::TSCOrAbstractVector, a2::TSCOrAbstractVector, as...) = length(a1) == length(a2) && lengths_match(a2, as...)
@inline lengths_match(a1::TSCOrAbstractVector, a2, as...) = lengths_match(a1, as...) # case: a2 is not indexable: skip it
@noinline lengths_match_fail() = throw(DimensionMismatch("Lengths of input collections do not match."))

@inline indices_match(::Val, indices::Vector{Int}, ::Any) = true
@inline function indices_match(::Val{i}, indices::Vector{Int}, tsc::TypeSortedCollection) where {i}
    tsc_indices = tsc.indices[i]
    length(indices) == length(tsc_indices) || return false
    @inbounds for j in eachindex(indices)
        indices[j] == tsc_indices[j] || return false
    end
    true
end
@inline indices_match(vali::Val, indices::Vector{Int}, a1, as...) = indices_match(vali, indices, a1) && indices_match(vali, indices, as...)
@noinline indices_match_fail() = throw(ArgumentError("Indices of TypeSortedCollections do not match."))

@generated function Base.map!(f::F, dest::TSCOrAbstractVector{N}, src1::TypeSortedCollection{<:Any, N}, srcs::TSCOrAbstractVector{N}...) where {F, N}
    expr = Expr(:block)
    for i = 1 : N
        vali = Val(i)
        push!(expr.args, quote
            let inds = leading_tsc.indices[$i]
                @boundscheck indices_match($vali, inds, dest, src1, srcs...) || indices_match_fail()
                @inbounds for j in LinearIndices(inds)
                    vecindex = inds[j]
                    _setindex!($vali, j, vecindex, dest, f(_getindex_all($vali, j, vecindex, src1, srcs...)...))
                end
            end
        end)
    end
    quote
        Base.@_inline_meta
        leading_tsc = first_tsc(dest, src1, srcs...)
        @boundscheck lengths_match(dest, src1, srcs...) || lengths_match_fail()
        $expr
        dest
    end
end

@generated function Base.foreach(f::F, A1::TypeSortedCollection{<:Any, N}, As::TSCOrAbstractVector{N}...) where {F, N}
    expr = Expr(:block)
    for i = 1 : N
        vali = Val(i)
        push!(expr.args, quote
            let inds = leading_tsc.indices[$i]
                @boundscheck indices_match($vali, inds, A1, As...) || indices_match_fail()
                @inbounds for j in LinearIndices(inds)
                    vecindex = inds[j]
                    f(_getindex_all($vali, j, vecindex, A1, As...)...)
                end
            end
        end)
    end
    quote
        Base.@_inline_meta
        leading_tsc = first_tsc(A1, As...)
        @boundscheck lengths_match(A1, As...) || lengths_match_fail()
        $expr
        nothing
    end
end

@generated function Base.mapreduce(f::F, op::O, tsc::TypeSortedCollection{<:Any, N}; init) where {F, O, N}
    quote
        Base.@_inline_meta
        ret = init
        $([:(ret = mapreduce(f, op, tsc.data[$i], init=ret)) for i = 1 : N]...)
        return ret
    end
end

@generated function Base.any(f::F, tsc::TypeSortedCollection{<:Any, N}) where {F, N}
    expr = Expr(:block)
    for i = 1 : N
        push!(expr.args, quote
            v = any(f, tsc.data[$i])
            if ismissing(v)
                anymissing = true
            elseif v
                return true
            end
        end)
    end
    quote
        Base.@_inline_meta
        anymissing = false
        $expr
        anymissing ? missing : false
    end
end

@generated function Base.all(f::F, tsc::TypeSortedCollection{<:Any, N}) where {F, N}
    expr = Expr(:block)
    for i = 1 : N
        push!(expr.args, quote
            v = all(f, tsc.data[$i])
            if ismissing(v)
                anymissing = true
            elseif !v
                return false
            end
        end)
    end
    quote
        Base.@_inline_meta
        anymissing = false
        $expr
        anymissing ? missing : true
    end
end

## broadcast!
@generated function _broadcast!(f::F, dest::D, A, Bs...) where {F, D}
    T = first_tsc_type(dest, A, Bs...)
    N = num_types(T)
    expr = Expr(:block)
    for i = 1 : N
        vali = Val(i)
        push!(expr.args, quote
            let inds = leading_tsc.indices[$i]
                @boundscheck indices_match($vali, inds, A, Bs...) || indices_match_fail()
                @inbounds for j in LinearIndices(inds)
                    vecindex = inds[j]
                    _setindex!($vali, j, vecindex, dest, f(_getindex_all($vali, j, vecindex, A, Bs...)...))
                end
            end
        end)
    end
    quote
        Base.@_inline_meta
        leading_tsc = first_tsc(dest, A, Bs...)
        @boundscheck lengths_match(dest, A, Bs...) || lengths_match_fail()
        $expr
        dest
    end
end

import Base.Broadcast: Broadcasted, BroadcastStyle, broadcastable, axes, instantiate
struct TypeSortedStyle <: Broadcast.BroadcastStyle end
BroadcastStyle(::Type{<:TypeSortedCollection}) = TypeSortedStyle()
BroadcastStyle(::TypeSortedStyle, ::Broadcast.DefaultArrayStyle{1}) = TypeSortedStyle()
BroadcastStyle(::TypeSortedStyle, ::Broadcast.DefaultArrayStyle{0}) = TypeSortedStyle()
broadcastable(x::TypeSortedCollection) = x
axes(tsc::TypeSortedCollection) = nothing
instantiate(bc::Broadcasted{TypeSortedStyle}) = bc
@inline Base.copyto!(dest::AbstractArray, bc::Broadcasted{TypeSortedStyle}) = _copy!(dest, bc)
@inline Base.copyto!(dest::TypeSortedCollection, bc::Broadcasted) = _copy!(dest, bc)
@inline function _copy!(dest, bc::Broadcasted)
    flat = Broadcast.flatten(bc)
    _broadcast!(flat.f, dest, flat.args...)
end

end # module


text:
using Distributed
using DistributedArrays

DA = DistributedArrays

"""
Split the work into `Task`s that are spawned on all workers available.
"""
struct DistributedTasks <: DistributionStrategy end

function step!(new::DArray, old::DArray, ::DistributedTasks)
    m, n = size(old)

    Indices = UnitRange{Int}
    Tile = Tuple{Indices,Indices}
    channels = [RemoteChannel(()->Channel{Tile}(Sys.CPU_THREADS)) for _ in procs(old)]

    indices = old.indices
    @assert size(indices) == size(channels)

    for (chan,indices) in zip(channels,indices)
        @async begin
            t_rows, t_cols = indices
            rend = min(m-1, last(t_rows))
            cend = min(n-1, last(t_cols))
            for c in max(2,first(t_cols)):TILESIZE:cend
                for r in max(2,first(t_rows)):TILESIZE:rend
                    rows = r:min(r+TILESIZE-1, rend)
                    cols = c:min(c+TILESIZE-1, cend)
                    I = (rows,cols)
                    put!(chan, I)
                end
            end
            close(chan)
        end
    end

    prepare!(old)

    @sync for (p,chan) in zip(procs(old),channels)
        # One could even use threads or local tasks on worker p.
        # In this scenario it's not useful though.
        @spawnat p begin
            for I in chan
                oldindices = I
                # Transform tile indices onto local indices:
                rows, cols = I
                offset = map(first,localindices(old))
                newindices = (first(rows)-offset[1]+1:last(rows)-offset[1]+1,
                              first(cols)-offset[2]+1:last(cols)-offset[2]+1)
                life_rule!(localpart(new), old, newindices, oldindices)
            end
        end
    end
end


text:
function bkkcrypt_encrypt(s)
  return s
end


text:
# Econometric Notes - Matlab code
#  Matteo Courthoud
#  01/01/2017

#import Pkg
#Pkg.add("LinearAlgebra")
#Pkg.add("Random")
#Pkg.add("Statistics")
#Pkg.add("Distributions")
#Pkg.add("DataFrames")
#Pkg.add("CSV")
#Pkg.add("Optim")

using LinearAlgebra
using Random
using Statistics
using Distributions
using DataFrames
using CSV
using Optim
using Plots

## Lecture 1

# Set seed
Random.seed!(123);

# Set the number of observations
n = 100;

# Set the dimension of X
k = 2;

# Draw a sample of explanatory variables
X = rand(Uniform(0,1), n, k);

# Draw the error term
σ = 1;
ε = rand(Normal(0,1), n, 1) * sqrt(σ);

# Set the parameters
β = [2; -1];

# Calculate the dependent variable
y = X*β + ε;

# Estimate beta
β_hat = inv(X'*X)*(X'*y)

# Equivalent but faster formulation
β_hat = (X'*X)\(X'*y);

# Even faster (but less intuitive) formulation
β_hat = X\y;

# Note that is generally not equivalent to Var(X)^-1 * Cov(X,y)...
β_wrong = inv(cov(X)) * cov(X, y)

# ...unless you include a constant
α = 3;
y1 = α .+ X*β + ε;
β_hat1 = [ones(n,1) X] \ y1
β_correct1 = inv(cov(X)) * cov(X, y1)

# Predicted y
y_hat = X*β_hat;

# Residuals
ε_hat = y - X*β_hat;

# Projection matrix
P = X * inv(X'*X) * X';

# Annihilator matrix
M = I - P;

# Leverage
h = diag(P);

# Biased variance estimator
σ_hat = ε_hat'*ε_hat / n;

# Unbiased estimator 1
σ_hat_2 = ε_hat'*ε_hat / (n-k);

# Unbiased estimator 2
σ_hat_3 = mean( ε_hat.^2 ./ (1 .- h) );

# R squared - uncentered
R2_uc = (y_hat'*y_hat)/ (y'*y);

# R squared
y_bar = mean(y);
R2 = ((y_hat .- y_bar)'*(y_hat .- y_bar))/ ((y .- y_bar)'*(y .- y_bar));

# Ideal variance of the OLS estimator
var_β = σ * inv(X'*X);

# Standard errors
std_β = sqrt.(diag(var_β));

## Lecture 2

# Set seed
Random.seed!(123);

# Homoskedastic standard errors
std_h = var(ε_hat) * inv(X'*X);

# HC0 variance and standard errors
Ω_hc0 = X' * (I(n) .* ε_hat.^2) * X;
std_hc0 = sqrt.(diag(inv(X'*X) * Ω_hc0 * inv(X'*X)))

# HC1 variance and standard errors
Ω_hc1 = n/(n-k) * X' * (I(n) .* ε_hat.^2) * X;
std_hc1 = sqrt.(diag(inv(X'*X) * Ω_hc1 * inv(X'*X)))

# HC2 variance and standard errors
Ω_hc2 = X' * (I(n) .* ε_hat.^2 ./ (1 .- h)) * X;
std_hc2 = sqrt.(diag(inv(X'*X) * Ω_hc2 * inv(X'*X)))

# HC3 variance and standard errors
Ω_hc3 = X' * (I(n) .* ε_hat.^2 ./ (1 .- h).^2) * X;
std_hc3 = sqrt.(diag(inv(X'*X) * Ω_hc3 * inv(X'*X)))

# Note what happens if you allow for full autocorrelation
omega_full = X'*ε_hat*ε_hat'*X;

# t-test for beta=0
t = abs.(β_hat ./ (std_hc1));

# p-value
p_val = 1 .- cdf.(Normal(0,1), t);

# F statistic of joint significance
SSR_u = ε_hat'*ε_hat;
SSR_r = y'*y;
F = (SSR_r - SSR_u)/k / (SSR_u/(n-k));

# 95# confidente intervals
conf_int = [β_hat - 1.96*std_hc1, β_hat + 1.96*std_hc1];

## Lecture 3

# Set seed
Random.seed!(123);

# Set the dimension of Z
l = 3;

# Draw instruments
Z = rand(Normal(0,1), n,l);

# Correlation matrix for error terms
S = [1 0.8; 0.8 1];

# Endogenous X
γ = [2 0; 0 -1; -1 3];
ε = rand(Normal(0,1), n, 2) * cholesky(S).U;
X = Z*γ .+ ε[:,1];

# Calculate y
y = X*β .+ ε[:,2];

# Estimate beta OLS
β_OLS = (X'*X)\(X'*y)

# IV: l=k=2 instruments
Z_IV = Z[:,1:k];
β_IV = (Z_IV'*X)\(Z_IV'*y)

# Calculate standard errors
ε_hat = y - X*β_IV;
V_NHC_IV = var(ε_hat) * inv(Z_IV'*X)*Z_IV'*Z_IV*inv(Z_IV'*X);
V_HC0_IV = inv(Z_IV'*X)*Z_IV' * (I(n) .* ε_hat.^2) * Z_IV*inv(Z_IV'*X);

# 2SLS: l=3 instruments
Pz = Z*inv(Z'*Z)*Z';
β_2SLS = (X'*Pz*X)\(X'*Pz*y)

# Calculate standard errors
ε_hat = y - X*β_2SLS;
V_NCH_2SLS = var(ε_hat) * inv(X'*Pz*X);
V_HC0_2SLS = inv(X'*Pz*X)*X'*Pz * (I(n) .* ε_hat.^2) *Pz*X*inv(X'*Pz*X);

# GMM 1-step: inefficient weighting matrix
W_1 = I(l);

# Objective function
gmm_1(b) = ( y - X*b )' * Z * W_1 *  Z' * ( y - X*b );

# Estimate GMM
β_gmm_1 = optimize(gmm_1, β_OLS).minimizer
ε_hat = y - X*β_gmm_1;

# Standard errors GMM
S_hat = Z' * (I(n) .* ε_hat.^2) * Z;
d_hat = -X'*Z;
V_gmm_1 = inv(d_hat * inv(S_hat) * d_hat');

# GMM 2-step: efficient weighting matrix
W_2 = inv(S_hat);

# Objective function
gmm_2(b) = ( y - X*b )' * Z * W_2 *  Z' * ( y - X*b );

# Estimate GMM
β_gmm_2 = optimize(gmm_2, β_OLS).minimizer

# Standard errors GMM
ε_hat = y - X*β_gmm_2;
S_hat = Z' * (I(n) .* ε_hat.^2) * Z;
d_hat = -X'*Z;
V_gmm_2 = inv(d_hat * inv(S_hat) * d_hat');

## Lecture 5

# Set seed
Random.seed!(123);

# Define Model Parameters
n = 1000;
α = 1;
β = .1;
γ = -4;

# Set number of simulations
n_simulations = 10000;

# Preallocate simulation results
α_short = zeros(n_simulations,1);
α_long = zeros(n_simulations,1);
α_pretest = zeros(n_simulations,1);

# Loop over simulations
for sim = 1:n_simulations

    # Generate Data
    Z = randn(n,1);
    X = γ*Z + randn(n,1);
    ε = randn(n,1);
    y = α*X + β*Z + ε;

    # Alpha estimate from long regression
    α_short[sim] = ((X'*X)\(X'*y))[1,1];

    # Alpha estimate from long regression
    estimates_temp = ([X Z]'*[X Z]) \ ([X Z]'*y);
    α_long[sim] = estimates_temp[1];

    # Compute test statistic
    e = y - [X Z]*estimates_temp;
    t = estimates_temp[2]/sqrt(var(ε)*([0 1]*inv([X Z]'*[X Z])*[0;1])[1]);

    # Pick estimate depending on test statistic
    if abs(t)>1.96
        α_pretest[sim] = α_long[sim];
    else
        α_pretest[sim] = α_short[sim];
    end
end

# Plot results
p1 = histogram(α_long, bins=30, title = "alpha long")
plot!([1], seriestype="vline", linewidth = 2)
p2 = histogram(α_short, bins=30, title = "alpha short")
plot!([1], seriestype="vline", linewidth = 2)
p3 = histogram(α_pretest, bins=30, title = "alpha pretest")
plot!([1], seriestype="vline", linewidth = 2)
p = plot(p1, p2, p3, layout=(1, 3), legend=false, size=(800, 300))
savefig(p, "../img/Fig_621.png")

# Sequence of sample sizes
n_sequence = [100,1000,10000,50000,100000];

# Initialize figure
figure
set(gcf,'position',[0,0,2000,1000])
n_fig = 1;

# Loop over data generating processes
for dgp=0:1
    for n=n_sequence
        for sim=1:n_simulations

            # Select beta according to dgp
            if dgp==0
                beta_n = beta;
            else
                beta_n = beta/sqrt(n)*sqrt(100);
            end

            # Generate Data
            Z = randn(n,1);
            X = γ*Z + randn(n,1);
            e = randn(n,1);
            y = alpha*X+ beta_n*Z + e;

            # Alpha estimate from long regression
            α_short[sim] = inv(X'*X)*X'*y;

            # Alpha estimate from long regression
            estimates_temp = inv([X,Z]'*[X,Z])*[X,Z]'*y;
            α_long[sim] = estimates_temp[1];

            # Compute test statistic
            e = y - [X,Z]*estimates_temp;
            t = estimates_temp[2]/sqrt(var(e)*([0,1]*inv([X,Z]'*[X,Z])*[0;1]));

            # Pick estimate depending on test statistic
            if abs(t)>1.96
                α_pretest[sim] = α_long[sim];
            else
                α_pretest[sim] = α_short[sim];
            end
        end

        # Add subplot
        subplot(2,length(n_sequence),n_fig)
        hist((α_pretest-alpha)*sqrt(n/100), 12)
        title(['beta=',num2str(beta_n)],'fontsize',10)
        xlabel(['n=',num2str(n)], 'fontsize',10)
        n_fig = n_fig+1;

    end
end
saveas(gcf,'figures/Fig_624.png')

max_stars_repo_name: matteocourthoud/Graduate-Econometrics

max_issues_repo_name: matteocourthoud/Graduate-Econometrics

max_forks_repo_name: matteocourthoud/Graduate-Econometrics


text:
using Swizzles.Properties
using Swizzles.Virtuals
using Swizzles.WrapperArrays
using Swizzles.ArrayifiedArrays
using Swizzles.StylishArrays
using Swizzles.ExtrudedArrays
using Swizzles.ValArrays
using Swizzles.NamedArrays
using Swizzles.ScalarArrays

using Base: checkbounds_indices, throw_boundserror, tail, dataids, unaliascopy, unalias
using Base.Iterators: reverse, repeated, countfrom, flatten, product, take, peel, EltypeUnknown
using Base.Broadcast: Broadcasted, BroadcastStyle, Style, DefaultArrayStyle, AbstractArrayStyle, Unknown, ArrayConflict
using Base.Broadcast: materialize, materialize!, instantiate, broadcastable, _broadcast_getindex, combine_eltypes, broadcast_shape, broadcast_unalias
using Base.FastMath: add_fast, mul_fast, min_fast, max_fast
using StaticArrays
using Base.Cartesian

struct SwizzledArray{T, N, Op, mask, Init<:AbstractArray, Arg<:AbstractArray} <: StylishArray{T, N}
    op::Op
    init::Init
    arg::Arg
    Base.@propagate_inbounds function SwizzledArray{T, N, Op, mask, Init, Arg}(op::Op, init::Init, arg::Arg) where {T, N, Op, mask, Init, Arg}
        @assert T isa Type
        @assert max(0, mask...) <= ndims(arg)
        @assert length(mask) == N
        #TODO assert mask is unique
        if op === nothing
            @boundscheck begin
                arg_keeps = keeps(arg)
                if any(imasktuple(d->kept(arg_keeps[d]), d->false, Val(mask), Val(ndims(arg))))
                    throw(DimensionMismatch("TODO"))
                end
            end
        end
        new(op, init, arg)
    end
end

@inline function SwizzledArray{T, N, Op, mask}(op::Op, init::Init, arg::Arg) where {T, N, Op, mask, Init, Arg}
    SwizzledArray{T, N, Op, mask, Init, Arg}(op, init, arg)
end

@inline function Base.convert(::Type{SwizzledArray{T}}, arr::SwizzledArray{S, N, Op, mask, Init, Arg}) where {T, S, N, Op, mask, Init, Arg}
    return SwizzledArray{T, N, Op, mask, Init, Arg}(arr.op, arr.init, arr.arg)
end

@inline function Properties.eltype_bound(arr::SwizzledArray)
    S = Properties.eltype_bound(arr.arg)
    if arr.op === nothing
        return S
    end
    T = Properties.eltype_bound(arr.init)
    T! = Union{T, Properties.return_type(arr.op, T, S)}
    if T! <: T
        return T!
    end
    arg_keeps = keeps(arr.arg)
    arr_mask = mask(arr)
    if all(imasktuple(d->arg_keeps[d] isa Extrude, d->true, Val(mask(arr)), Val(ndims(arr.arg))))
        return T!
    end
    T = T!
    T! = Union{T, Properties.return_type(arr.op, T, S)}
    if T! <: T
        return T!
    end
    return Any
end

@inline mask(::Type{<:SwizzledArray{<:Any, <:Any, <:Any, _mask}}) where {_mask} = _mask
@inline mask(::SwizzledArray{<:Any, <:Any, <:Any, _mask}) where {_mask} = _mask

function Base.show(io::IO, arr::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    print(io, SwizzledArray)
    print(io, "{$T, $N, $Op, $mask}($(arr.op), $(arr.init), $(arr.arg))")
    nothing
end

Base.parent(arr::SwizzledArray) = arr.arg
Base.parent(::Type{<:SwizzledArray{T, N, Op, mask, Init, Arg}}) where {T, N, Op, mask, Init, Arg} = Arg
WrapperArrays.iswrapper(arr::SwizzledArray) = true
Base.@propagate_inbounds function WrapperArrays.adopt(arr::SwizzledArray{T, N, Op, mask, Init}, arg::Arg) where {T, N, Op, mask, Init, Arg}
    SwizzledArray{T, N, Op, mask, Init, Arg}(arr.op, arr.init, arg)
end

Base.dataids(src::SwizzledArray) = (dataids(src.op), dataids(src.init), dataids(src.arg))
function Base.unaliascopy(src::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    init = unaliascopy(src.init)
    arg = unaliascopy(src.arg)
    SwizzledArray{T, N, Op, mask, typeof(init), typeof(arg)}(src.op, init, arg)
end
function Base.unalias(dst, src::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    init = unalias(dst, src.init)
    arg = unalias(dst, src.arg)
    SwizzledArray{T, N, Op, mask, typeof(init), typeof(arg)}(src.op, init, arg)
end
function Broadcast.broadcast_unalias(dst, src::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    init = broadcast_unalias(dst, src.init)
    if ndims(src.arg) != 0
        arg = unalias(dst, src.arg)
    else
        arg = src.arg
    end
    SwizzledArray{T, N, Op, mask, typeof(init), typeof(arg)}(src.op, init, arg)
end
function Broadcast.broadcast_unalias(::Nothing, src::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    return src
end
function ArrayifiedArrays.preprocess(dst, src::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    init = src.init #preprocess that later
    if ndims(src.arg) != 0
        #arg = preprocess(nothing, unalias(dst, src.arg)) #FIXME for some reason, calling mightalias screws up vectorization. I do not understand.
        arg = preprocess(nothing, src.arg)
    else
        arg = src.arg
    end
    SwizzledArray{T, N, Op, mask, typeof(init), typeof(arg)}(src.op, init, arg)
end

@inline function Base.size(arr::SwizzledArray)
    arg_size = size(arr.arg)
    masktuple(d->1, d->arg_size[d], Val(mask(arr)))
end

@inline function Base.axes(arr::SwizzledArray)
    arg_axes = axes(arr.arg)
    masktuple(d->Base.OneTo(1), d->arg_axes[d], Val(mask(arr)))
end

Base.@propagate_inbounds function Base.convert(::Type{SwizzledArray}, src::SubArray{T, M, Arr, <:Tuple{Vararg{Any, N}}}) where {T, N, M, Op, Arr <: SwizzledArray{T, N, Op}}
    arr = parent(src)
    inds = Base.parentindices(src)
    arg = arr.arg
    init = arr.init
    init′ = SubArray(init, ntuple(n -> (Base.@_inline_meta; size(init, n) == 1 ? firstindex(init, n) : inds[n]), Val(ndims(init))))
    arg_axes = axes(arr.arg)
    inds′ = imasktuple(d->Base.Slice(arg_axes[d]), d->inds[d], Val(mask(arr)), Val(ndims(arr.arg)))
    arg′ = SubArray(arg, inds′)
    mask′ = remask(inds, inds′, mask(arr))
    return SwizzledArray{eltype(src), M, Op, mask′}(arr.op, init′, arg′)
end

@inline function remask(inds, inds′, mask::Tuple{})
    ()
end
@inline function remask(inds, inds′, mask)
    counts = _remask_counts((), inds′)
    _remask_mask(counts, inds, mask)
end
@inline function _remask_mask(counts, inds, mask)
    rest = _remask_mask(counts, Base.tail(inds), Base.tail(mask))
    if Base.index_dimsum(first(inds)) isa Tuple{}
        return rest
    elseif first(mask) === nil
        return (nil, rest...)
    else
        return (counts[first(mask)], rest...)
    end
end
@inline _remask_mask(counts, ::Tuple{}, ::Tuple{}) = ()
@inline function _remask_counts(firsts, inds′)
    firsts = (firsts..., first(inds′))
    return (length(Base.index_dimsum(firsts...)), _remask_counts(firsts, Base.tail(inds′))...)
end
@inline _remask_counts(counts, ::Tuple{}) = ()

Base.similar(::Broadcasted{DefaultArrayStyle{0}, <:Any, typeof(identity), <:Tuple{<:SwizzledArray{T}}}) where {T} = ScalarArray{T}()

Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, sty::Styled{<:AbstractArrayStyle, <:SubArray{T, <:Any, <:SwizzledArray{T, N}, <:Tuple{Vararg{Any, N}}}}) where {T, N}
    #A view of a Swizzle can be computed as a swizzle of a view (hiding the
    #complexity of dropping view indices). Therefore, we convert first.
    sub = sty.arg
    src = convert(SwizzledArray, sub)
    return Base.copyto!(dst, src)
end

Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, sty::Styled{<:AbstractArrayStyle, <:SubArray{T, 0, <:SwizzledArray{T, 0}, <:Tuple{}}}) where {T}
    #Zero-dimensional reswizzling case.
    src = parent(sty.arg)
    return Base.copyto!(dst, src)
end

Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, sty::Styled{<:AbstractArrayStyle, <:SwizzledArray})
    @boundscheck axes(dst) == axes(sty.arg) || error("TODO")
    if eltype(sty.arg) <: eltype(dst)
        src = preprocess(dst, sty.arg)
        arg = src.arg
        op = src.op
        init = src.init
        @inbounds begin
            index = swizzleindex(dst, src)
            drive = eachindex(arg, index)
            if op === nothing
                assign!(dst, index, arg, drive)
            else
                dst .= init
                increment!(op, dst, index, arg, drive)
            end
        end
        return dst
    else
        #if the destination is unsuitable for directly accumulating the swizzle, we just give up and copy it.
        copyto!(dst, copy(sty))
    end
end

#Broadcast Optimizations
for Style = (AbstractArrayStyle{0}, AbstractArrayStyle, DefaultArrayStyle{0}, DefaultArrayStyle)
    for Identity = (typeof(identity), typeof(myidentity))
        @eval begin
            Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, bc::Broadcasted{S, <:Any, $Identity, <:Tuple{SubArray{T, <:Any, <:SwizzledArray{T, N}, <:Tuple{Vararg{Any, N}}}}}) where {S<:$Style, T, N}
                if axes(dst) == axes(bc.args[1]) && eltype(bc.args[1]) <: eltype(dst)
                    copyto!(dst, Styled{S}(bc.args[1]))
                else
                    dst .= copy(bc.args[1])
                end
            end

            Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, bc::Broadcasted{S, <:Any, $Identity, <:Tuple{SwizzledArray}}) where {S<:$Style}
                if axes(dst) == axes(bc.args[1]) && eltype(bc.args[1]) <: eltype(dst)
                    copyto!(dst, Styled{S}(bc.args[1]))
                else
                    dst .= copy(bc.args[1])
                end
            end
        end
    end
end

is_nil_mask(mask) = mask == ntuple(n->nil, length(mask))
@generated function is_nil_mask(::Val{mask}) where {mask}
    return is_nil_mask(mask)
end

is_oneto_mask(mask) = mask == 1:length(mask)
@generated function is_oneto_mask(::Val{mask}) where {mask}
    return is_oneto_mask(mask)
end

Base.@propagate_inbounds indices(arr) = indices(IndexStyle(arr), arr)
Base.@propagate_inbounds indices(::IndexLinear, arr) = LinearIndices(arr)
Base.@propagate_inbounds indices(::IndexCartesian, arr) = CartesianIndices(arr)

Base.@propagate_inbounds swizzleindex(dst, arr) = swizzleindex(IndexStyle(dst), dst, arr)

Base.@propagate_inbounds function swizzleindex(::IndexLinear, dst, arr)
    if is_nil_mask(Val(mask(arr)))
        return ConstantIndices(1, indices(arr.arg))
    elseif is_oneto_mask(Val(mask(arr)))
        return LinearIndices(dst)
    else
        return SwizzledIndices(arr)
    end
end

Base.@propagate_inbounds function swizzleindex(::IndexCartesian, dst, arr)
    if is_nil_mask(Val(mask(arr)))
        i = CartesianIndex(ntuple(n->1, Val(ndims(dst))))
        return ConstantIndices(i, indices(arr.arg))
    elseif is_oneto_mask(Val(mask(arr)))
        return CartesianIndices(dst)
    else
        return SwizzledIndices(arr)
    end
end

#FIXME This is a separate method for the scalar case so that swizzles can avoid
#recursion depth limiting.
Base.@propagate_inbounds function assign!(dst, index, src, drive::CartesianIndices{0})
    i = CartesianIndex()
    dst[index[i]] = src[i]
end

#FIXME This is a separate method for the scalar case so that swizzles can avoid
#recursion depth limiting.
Base.@propagate_inbounds function increment!(op::Op, dst, index, src, drive::CartesianIndices{0}) where {Op}
    i = CartesianIndex()
    i′ = index[i]
    dst[i′] = op(dst[i′], src[i])
end

Base.@propagate_inbounds function assign!(dst, index, src, drive)
    for i in drive
        dst[index[i]] = src[i]
    end
end

Base.@propagate_inbounds function increment!(op::Op, dst, index, src, drive) where {Op}
    for i in drive
        i′ = index[i]
        dst[i′] = op(dst[i′], src[i])
    end
end

@generated function assign!(dst, index, src, drive::CartesianIndices{N}) where {N}
    return quote
        Base.@_propagate_inbounds_meta
        @nloops $N i n -> drive.indices[n] begin
            i = CartesianIndex(@ntuple $N i)
            dst[index[i]] = src[i]
        end
    end
end

@generated function increment!(op::Op, dst, index, src, drive::CartesianIndices{N}) where {Op, N}
    return quote
        Base.@_propagate_inbounds_meta
        @nloops $N i n -> drive.indices[n] begin
            i = CartesianIndex(@ntuple $N i)
            i′ = index[i]
            dst[i′] = op(dst[i′], src[i])
        end
    end
end

struct SwizzledIndices{T, N, mask, Inds <: CartesianIndices{N}} <: AbstractArray{T, N}
    inds::Inds
    function SwizzledIndices{T, N, mask, Inds}(inds::Inds) where {T, N, mask, Inds}
        @assert T <: CartesianIndex{length(mask)}
        @assert max(0, mask...) <= ndims(inds)
        #TODO assert mask is unique
        return new{T, N, mask, Inds}(inds)
    end
end

Base.@propagate_inbounds SwizzledIndices(arr::SwizzledArray) = SwizzledIndices(Val(mask(arr)), CartesianIndices(arr.arg))
Base.@propagate_inbounds SwizzledIndices(mask, inds) = SwizzledIndices(mask, CartesianIndices(inds))
Base.@propagate_inbounds SwizzledIndices(mask, inds::CartesianIndices) = SwizzledIndices(Val(mask), inds)
Base.@propagate_inbounds function SwizzledIndices(::Val{mask}, inds::CartesianIndices) where {mask}
    return SwizzledIndices{CartesianIndex{length(mask)}, ndims(inds), mask, typeof(inds)}(inds)
end

@inline Base.size(arr::SwizzledIndices) = size(arr.inds)

@inline Base.axes(arr::SwizzledIndices) = axes(arr.inds)

@inline mask(arr::SwizzledIndices{<:Any, <:Any, _mask}) where {_mask} = _mask
@inline mask(::Type{<:SwizzledIndices{<:Any, <:Any, _mask}}) where {_mask} = _mask

Base.@propagate_inbounds function Base.getindex(arr::SwizzledIndices{T, N}, i::Vararg{Any, N})::T where {T, N}
    return CartesianIndex(masktuple(d->1, d->i[d], Val(mask(arr))))
end

struct ConstantIndices{T, N, i, Inds <: AbstractArray{<:Any, N}} <: AbstractArray{T, N}
    inds::Inds
    @inline function ConstantIndices{T, N, i, Inds}(inds::Inds) where {T, N, i, Inds}
        @assert i isa T
        return new{T, N, i, Inds}(inds)
    end
end

Base.@propagate_inbounds ConstantIndices(i, inds) = ConstantIndices(Val(i), inds)
Base.@propagate_inbounds function ConstantIndices(::Val{i}, inds) where {i}
    return ConstantIndices{typeof(i), ndims(inds), i, typeof(inds)}(inds)
end

@inline Base.size(arr::ConstantIndices) = size(arr.inds)
@inline Base.axes(arr::ConstantIndices) = axes(arr.inds)
@inline Base.IndexStyle(arr::ConstantIndices) = IndexStyle(arr.inds)
Base.@propagate_inbounds (Base.getindex(arr::ConstantIndices{T, N, i}, ::Vararg{Any, N})::T) where {T, N, i} = i
Base.@propagate_inbounds (Base.getindex(arr::ConstantIndices{T, 1, i}, ::Integer)::T) where {T, i} = i
Base.@propagate_inbounds (Base.getindex(arr::ConstantIndices{T, <:Any, i}, ::Integer)::T) where {T, i} = i

@inline childstyle(Arr::Type{<:SwizzledArray}, ::DefaultArrayStyle) = DefaultArrayStyle{ndims(Arr)}()
@inline childstyle(Arr::Type{<:SwizzledArray}, ::BroadcastStyle) = DefaultArrayStyle{ndims(Arr)}()
@inline childstyle(::Type{<:SwizzledArray}, ::ArrayConflict) = ArrayConflict()
@inline childstyle(Arr::Type{<:SwizzledArray}, ::Style{Tuple}) = mask(Arr) == (1,) ? Style{Tuple}() : DefaultArrayStyle{ndims(Arr)}()

@inline function Broadcast.BroadcastStyle(Arr::Type{<:SwizzledArray})
    childstyle(Arr, BroadcastStyle(parent(Arr)))
end

@inline function ExtrudedArrays.keeps(arr::SwizzledArray)
    arg_keeps = keeps(arr.arg)
    arr_keeps = masktuple(d->extrude, d->arg_keeps[d], Val(mask(arr)))
    init_keeps = keeps(arr.init)
    return combinetuple(|, arr_keeps, init_keeps)
end

function ExtrudedArrays.freeze_extrudes(arr::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    return SwizzledArray{T, N, Op, mask}(arr.op, freeze_extrudes(arr.init), freeze_extrudes(arr.arg)) #FIXME this doesn't stabilize the swizzle though
end

function ExtrudedArrays.lift_keeps(arr::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    return SwizzledArray{T, N, Op, mask}(arr.op, lift_keeps(arr.init), lift_keeps(arr.arg))
end

function ValArrays.lift_vals(arr::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    return SwizzledArray{T, N, Op, mask}(arr.op, lift_vals(arr.init), lift_vals(arr.arg))
end

function NamedArrays.lift_names(arr::SwizzledArray{T, N, Op, mask}) where {T, N, Op, mask}
    return SwizzledArray{T, N, Op, mask}(arr.op, lift_names(arr.init), lift_names(arr.arg))
end

function Virtuals.virtualize(root, ::Type{<:SwizzledArray{T, N, Op, mask, Init, Arg}}) where {T, N, Op, mask, Init, Arg}
    init = virtualize(:($root.init), Init)
    arg = virtualize(:($root.arg), Arg)
    op = virtualize(:($root.op), Op)
    return SwizzledArray{T, N, Op, mask}(op, init, arg)
end


text:
function dict_add_reduce(local_res)
  res = reduce(local_res) do  rr, rr1
    for k = keys(rr1)
      rr[k] = get(rr, k, 0) + rr1[k]
    end
    return rr
  end
  return res
end

function dict_mean_reduce(local_res)
  res = reduce(local_res) do rr, rr1
    for k in keys(rr1)
      a = get(rr,k,(0,0))
      b = rr1[k]
      c = (a[1] + b[1], a[2] + b[2])
      rr[k] = c
    end
    rr
  end
  res
end

function pgroupreduce{T}(byfn, map, local_reduce, master_reduce, by::SharedArray{T,1}, val::SharedArray{T,1})
  if nprocs()==1
    throw(ErrorException("Only 1 worker; use groupreduce instead"))
  end
  l = length(by)
  ii = sort(collect(Set([1:Int64(round(l/nprocs())):l...,l])))
  res = pmap(2:length(ii)) do i
    iii = ii[i-1]:ii[i]
    groupreduce(byfn, map, local_reduce, zip(by[iii],val[iii]))
  end
  return master_reduce(res)
end


text:
module Revise354

test() = println("Original")

end # module


text:
#####Convenience functions for filtering arrays of Imagine commands by channel name, channel type######

findname(coms::AbstractVector{T}, nm::AbstractString) where{T<:ImagineSignal} =
    findfirst(x->name(x) == nm, coms)
function getname(coms::AbstractVector{T}, nm::AbstractString) where{T<:ImagineSignal}
    namei = findname(coms, nm)
    if namei == nothing
        error("Name $nm not found")
    else
        coms[namei]
    end
end
finddaqchan(coms::AbstractVector{T}, nm::AbstractString) where{T<:ImagineSignal} =
    findfirst(x->daq_channel(x) == nm, coms)
getdaqchan(coms::AbstractVector{T}, nm::AbstractString) where{T<:ImagineSignal} =
    coms[finddaqchan(coms, nm)]

isdigital(com::T) where{T<:ImagineSignal} =
    isdigital(daq_channel(com), rig_name(com))
finddigital(coms::AbstractVector{T}) where{T<:ImagineSignal} =
    findall(isdigital, coms)
getdigital(coms::AbstractVector{T}) where{T<:ImagineSignal} =
    view(coms, finddigital(coms))
isanalog(daq_chan::AbstractString, rig::AbstractString) where{T<:ImagineSignal} =
    !isdigital(daq_chan, rig)
isanalog(com::T) where{T<:ImagineSignal} =
    isanalog(daq_channel(com), rig_name(com))
findanalog(coms::AbstractVector{T}) where{T<:ImagineSignal} =
    findall(x->!isdigital(x), coms)
getanalog(coms::AbstractVector{T}) where{T<:ImagineSignal} =
    view(coms, findanalog(coms))

isoutput(com::T) where{T<:ImagineSignal} = isoutput(daq_channel(com), rig_name(com))
isinput(com::T) where{T<:ImagineSignal} = !isoutput(com)
findoutputs(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(isoutput, coms)
getoutputs(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findoutputs(coms))
findinputs(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(isinput, coms)
getinputs(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findinputs(coms))

isfree(com::T) where{T<:ImagineSignal} = isfree(daq_channel(com), rig_name(com))
isfixed(com::T) where{T<:ImagineSignal} = !isfree(com)
findfree(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(isfree, coms)
getfree(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findfree(coms))
findfixed(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(x->!isfree(x), coms)
getfixed(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findfixed(coms))

ispos(com::T) where{T<:ImagineSignal} = ispos(daq_channel(com), rig_name(com))
findpositioners(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(ispos, coms)
getpositioners(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findpositioners(coms))

isposmonitor(com::T) where{T<:ImagineSignal} = isposmonitor(daq_channel(com), rig_name(com))
findpositionermonitors(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(isposmonitor, coms)
getpositionermonitors(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findpositionermonitors(coms))

iscam(com::T) where{T<:ImagineSignal} = iscam(daq_channel(com), rig_name(com))
findcameras(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(iscam, coms)
getcameras(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findcameras(coms))

iscammonitor(com::T) where{T<:ImagineSignal} = iscammonitor(daq_channel(com), rig_name(com))
findcameramonitors(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(iscammonitor, coms)
getcameramonitors(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findcameramonitors(coms))

islas(com::T) where{T<:ImagineSignal} = islas(daq_channel(com), rig_name(com))
findlasers(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(islas, coms)
getlasers(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findlasers(coms))

isstim(com::T) where{T<:ImagineSignal} = isstim(daq_channel(com), rig_name(com))
findstimuli(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(isstim, coms)
getstimuli(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findstimuli(coms))

isgalvo(com::T) where{T<:ImagineSignal} = iscam(daq_channel(com), rig_name(com))
findgalvos(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(iscam, coms)
getgalvos(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findcameras(coms))

isgalvomonitor(com::T) where{T<:ImagineSignal} = isgalvomonitor(daq_channel(com), rig_name(com))
findgalvomonitors(coms::AbstractVector{T}) where{T<:ImagineSignal} = findall(isgalvomonitor, coms)
getgalvomonitors(coms::AbstractVector{T}) where{T<:ImagineSignal} = view(coms, findgalvomonitors(coms))

hasmonitor(com::T) where{T<:ImagineSignal} = iscam(com) || ispos(com) || isgalvo(com)
hasactuator(com::T) where{T<:ImagineSignal} = iscammonitor(com) || isposmonitor(com) || isgalvomonitor(com)

function monitor_name(com::T) where{T<:ImagineSignal}
    if !hasmonitor(com)
        error("There is no monitor (input) corresponding to this channel")
    end
    if iscam(com)
        return name(com) * " frame monitor"
    else
        return name(com) * " monitor"
    end
end
function getmonitor(com::T) where{T<:ImagineSignal}
    mon_nm = monitor_name(com)
    temp = rigtemplate(rig_name(com); sample_rate = samprate(com))
    return getname(temp, mon_nm)
end

function actuator_name(com::T) where{T<:ImagineSignal}
    if !hasactuator(com)
        error("There is no actuator (output) corresponding to this channel")
    end
    if iscammonitor(com)
        return String(split(name(com), " frame monitor")[1])
    else
        return String(split(name(com), " monitor")[1])
    end
end
function getactuator(com::T) where{T<:ImagineSignal}
    act_nm = actuator_name(com)
    temp = rigtemplate(rig_name(com); sample_rate = samprate(com))
    return getname(temp, act_nm)
end

max_forks_repo_name: HolyLab/ImagineInterface


text:
using Printf
using Turing
using CairoMakie
using ElectronDisplay

# Data generation

N = 100
x = range(0, 1, length=N)
y = @. 0.16x^2 + 4.5x + π
y = @. y * (1 + 0.1 * randn())

scatter(x, y)

# Linear model (y = mx + c)

@model function linear_model(x, y)
    m ~ Normal(0, 1)
    c ~ Normal(0, 1)
    σ ~ Exponential(1)

    N = length(x)
    for n in 1:N
        y[n] ~ Normal(m*x[n] + c, σ)
    end
end

chain_linear = sample(linear_model(x, y), NUTS(0.65), 10000, progress=true)

# Quadratic model (y = ax² + bx + c)

@model function quadratic_model(x, y)
    a ~ Normal(0, 1)
    b ~ Normal(0, 1)
    c ~ Normal(0, 1)
    σ ~ Exponential(1)

    N = length(x)
    for n in 1:N
        y[n] ~ Normal(a*x[n]^2 + b*x[n] + c, σ)
    end
end

chain_quadratic = sample(quadratic_model(x, y), NUTS(0.65), 10000, progress=true)

# Exponential model [y = α + β * exp(γ * x)]

@model function exponential_model(x, y)
    α ~ Normal(0, 1)
    β ~ Normal(0, 1)
    γ ~ Normal(0, 1)
    σ ~ Exponential(1)

    N = length(x)
    for n in 1:N
        y[n] ~ Normal(α + β * exp(γ*x[n]), σ)
    end
end

chain_exponential = sample(exponential_model(x, y), NUTS(0.65), 10000, progress=true)

# Approximate Bayesian computation

loss(y₁, y₂) = mean((y₂ .- y₁).^2) / length(y₁)

function rand_model(chain, θs)
    N = length(chain)
    n = rand(1:N)
    return NamedTuple(θ => chain[θ][n] for θ in θs)
end

function rand_y_linear(chain, x)
    m, c = rand_model(chain, (:m, :c))
    return @. m*x + c
end

function rand_y_quadratic(chain, x)
    a, b, c = rand_model(chain, [:a, :b, :c])
    return @. a*x^2 + b*x + c
end

function rand_y_exponential(chain, x)
    α, β, γ = rand_model(chain, [:α, :β, :γ])
    return @. α + β*exp(γ*x)
end

@enum LineModel linear=1 quadratic=2 exponential=3

function abc_line_model(x, y, M₀, chain_linear, chain_quadratic, chain_exponential; N_samples, max_loss)
    M = zeros(Int, N_samples)
    samples = 0
    trials = 0

    while samples < N_samples
        model = rand(M₀) |> LineModel

        if model == linear
            l = loss(y, rand_y_linear(chain_linear, x))
        elseif model == quadratic
            l = loss(y, rand_y_quadratic(chain_quadratic, x))
        elseif model == exponential
            l = loss(y, rand_y_exponential(chain_exponential, x))
        end

        if l <= max_loss
            samples += 1
            M[samples] = Int(model)
        end

        trials += 1
        @info @sprintf("Trial %d, sample %d: %s, loss=%.3e (%s)", trials, samples, model, l, l <= max_loss ? "accepted!" : "rejected")
    end

    @info *([@sprintf("P(m=%s)=%.4f%s", LineModel(m), count(M .== m) / N_samples, m < 3 ? ", " : "") for m in 1:3]...)

    return M
end

losses_linear = [loss(y, rand_y_linear(chain_linear, x)) for _ in 1:1000]
losses_quadratic = [loss(y, rand_y_quadratic(chain_quadratic, x)) for _ in 1:1000]
losses_exponential = [loss(y, rand_y_exponential(chain_exponential, x)) for _ in 1:1000]

M₀ = Categorical(3)  # Prior distribution over models?
M = abc_line_model(x, y, M₀, chain_linear, chain_quadratic, chain_exponential, N_samples=10000, max_loss=0.0034)


text:
import Base: rand
using Random: AbstractRNG, Sampler, SamplerType

rand(rng::AbstractRNG, ::SamplerType{Leaf{N}}) where {N} = Leaf{N}(rand(rng, 1:N))

"""
    BoltzmannSampler{N}(m, n)
Configuration for Boltzmann sampling `DecisionTree{N}`s with minimum size `m` and maximum size `n`.
"""
struct BoltzmannSampler{N} <: Sampler{Tree{N}}
    minsize::Int
    maxsize::Int
end

function boltzmann_ub(rng::AbstractRNG, sampler::BoltzmannSampler{N}, cursize) where {N}
    minsize = sampler.minsize
    maxsize = sampler.maxsize

    if cursize ≥ maxsize        # maximum bound exceeded
        return nothing, cursize
    elseif rand(rng) ≤ 0.5      # generate random leaf
        return rand(rng, Leaf{N}), cursize
    else                        # try generating a branch -- if its size is small enough
        # try generating children, or immediately fail
        left_children, cursize = boltzmann_ub(rng, sampler, cursize + 1)
        left_children === nothing && return nothing, cursize

        right_children, cursize = boltzmann_ub(rng, sampler, cursize + 1)
        right_children === nothing && return nothing, cursize

        # `cursize` is now the size of the whole thing, and less than maxsize
        return Branch{N}(left_children, right_children), cursize
    end
end

function rand(rng::AbstractRNG, sampler::BoltzmannSampler{N}) where {N}
    candidate, size = boltzmann_ub(rng, sampler, 0)
    while candidate === nothing || size ≤ sampler.minsize
        candidate, size = boltzmann_ub(rng, sampler, 0)
    end

    return candidate
end


text:
# This file was generated, do not modify it. # hide
url = "https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data";
feature_names = ["ID", "Class", "mean radius", "mean texture", "mean perimeter", "mean area", "mean smoothness", "mean compactness", "mean concavity", "mean concave points", "mean symmetry", "mean fractal dimension", "radius error", "texture error", "perimeter error", "area error", "smoothness error", "compactness error", "concavity error", "concave points error", "symmetry error", "fractal dimension error", "worst radius", "worst texture", "worst perimeter", "worst area", "worst smoothness", "worst compactness", "worst concavity", "worst concave points", "worst symmetry", "worst fractal dimension"]
data = urldownload(url, true, format = :CSV, header = feature_names);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function import_edcd_data(filepath = "", url = "https://opendata.ecdc.europa.eu/covid19/casedistribution/csv")
    if filepath == "" && startswith(url, "https://")
        filepath = string(dirname(@__DIR__)) * "/data/raw/latest_edcd_covid.csv"
        download(url, filepath)
    end

    raw_dat = DataFrame(CSV.File(filepath))

    raw_dat.dateRep = Date.(raw_dat[:, :dateRep], "dd/mm/yyyy")
    sort!(raw_dat, (:countriesAndTerritories, :dateRep));
    filter!(row -> !ismissing(row.popData2018), raw_dat);

    function add_cumsum_time(df)
        cum_case_prop = cumsum(df.cases ./ df.popData2018)
        t = df.dateRep - minimum(df.dateRep)
        t = map(x -> convert(Float64, x.value), t)
        return (cum_case_prop = cum_case_prop, time_since_first_case = t, popData2018 = df.popData2018)
    end

    dat = @pipe (groupby(raw_dat, :countriesAndTerritories) |>
        map(add_cumsum_time, _) |>
        DataFrame(_));
    filter!(row -> row.cum_case_prop > 0.0, dat);

    return dat
end

function import_hopkins_data(filepath, value_col_name)

    function aggregate_across_regions(df)
        colnames = (Symbol(value_col_name), :date)
        colvalues = (sum(df[:, Symbol(value_col_name)]), df.date[1])
        res = namedtuple(colnames, colvalues)
        return res
    end

    dat = @pipe (DataFrame(CSV.File(filepath)) |>
        stack(_, 5:size(_)[2]) |>
        rename(_, Dict( "Country/Region" => "Country_Region",
            "Province/State" => "Province_State",
            "variable" => "date",
            "value" => value_col_name)) |>
            groupby(_, [:Country_Region, :date]) |>
            map(aggregate_across_regions, _) |>
            DataFrame(_))
    return dat

end

function import_data(drop_threshold = 100; download_data = false)

    confirmed_fn = string(dirname(@__DIR__)) * "/data/raw/time_series_covid19_confirmed_global.csv"
    deaths_fn = string(dirname(@__DIR__)) * "/data/raw/time_series_covid19_deaths_global.csv"
    recovered_fn = string(dirname(@__DIR__)) * "/data/raw/time_series_covid19_recovered_global.csv"

    if download_data
        download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv",
            confirmed_fn)
        download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv",
            deaths_fn)
        download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv",
            recovered_fn)
    end

    hopkins_case_dat = import_hopkins_data(confirmed_fn, "confirmed_cases")
    hopkins_death_dat = import_hopkins_data(deaths_fn, "deaths")
    hopkins_recov_dat = import_hopkins_data(recovered_fn, "recovered")

    hopkins_dat = @pipe (join(hopkins_case_dat, hopkins_recov_dat, on = [:date, :Country_Region]) |>
        join(_, hopkins_death_dat, on = [:date, :Country_Region]))

    edcd_dat = @pipe(import_edcd_data() |>
        select(_, :countriesAndTerritories, :popData2018) |>
        groupby(_, :countriesAndTerritories) |>
        map(df -> (countriesAndTerritories = df.countriesAndTerritories[1], popData2018 = df.popData2018[1]), _) |>
        DataFrame(_))

    replace!(edcd_dat.countriesAndTerritories, "United_States_of_America" => "US")

    hopkins_ecdc = join(hopkins_dat, edcd_dat, on = [:Country_Region => :countriesAndTerritories])

    # drop observations bewlow 100 cases
    filter!(row -> row.confirmed_cases > drop_threshold, hopkins_ecdc);

    function add_SIR_time(df)
        case_prop = df.confirmed_cases ./ df.popData2018
        cum_case_prop = cumsum(df.confirmed_cases ./ df.popData2018)
        death_prop = df.deaths ./ df.popData2018
        recov_prop = df.recovered ./ df.popData2018
        susc_prop = ((df.popData2018 - df.confirmed_cases - df.deaths - df.recovered) ./ df.popData2018)

        dates = Date.(string.(df.date), "mm/dd/yy") + Dates.Year("2000")
        t = dates - minimum(dates)
        res = (time_since_case_thresh = map(x -> convert(Float64, x.value), t),
            case_prop = case_prop,
            cum_case_prop = cum_case_prop,
            death_prop = death_prop,
            recov_prop = recov_prop,
            susc_prop = susc_prop,
            pop = df.popData2018)
        return res
    end

    dat = @pipe (groupby(hopkins_ecdc, :Country_Region) |>
        map(add_SIR_time, _) |>
        DataFrame(_));

    return dat
end


text: fCheck2(x) = x^2

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# Autogenerated wrapper script for diffutils_jll for x86_64-linux-gnu
export _cmp, _diff, diff3, sdiff

JLLWrappers.@generate_wrapper_header("diffutils")
JLLWrappers.@declare_executable_product(_cmp)
JLLWrappers.@declare_executable_product(_diff)
JLLWrappers.@declare_executable_product(diff3)
JLLWrappers.@declare_executable_product(sdiff)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        _cmp,
        "bin/cmp",
    )

    JLLWrappers.@init_executable_product(
        _diff,
        "bin/diff",
    )

    JLLWrappers.@init_executable_product(
        diff3,
        "bin/diff3",
    )

    JLLWrappers.@init_executable_product(
        sdiff,
        "bin/sdiff",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()