text:
using ObjectFile.ELF
import ObjectFile: rpaths, canonical_rpaths

"""
    platform_for_object(oh::ObjectHandle)

Returns the platform the given `ObjectHandle` should run on.  E.g.
if the given `ObjectHandle` is an x86_64 Linux ELF object, this function
will return `Linux(:x86_64)`.  This function does not yet distinguish
between different libc's such as `:glibc` and `:musl`.
"""
function platform_for_object(oh::ObjectHandle)
    if oh isa ELFHandle
        mach_to_arch = Dict(
            ELF.EM_386 => :i686,
            ELF.EM_X86_64 => :x86_64,
            ELF.EM_AARCH64 => :aarch64,
            ELF.EM_PPC64 => :powerpc64le,
            ELF.EM_ARM => :armv7l,
        )
        mach = oh.header.e_machine
        if !haskey(mach_to_arch, mach)
            error("Unknown ELF architecture $(mach)")
        end

        arch = mach_to_arch[mach]

        if oh.ei.osabi == ELF.ELFOSABI_LINUX || oh.ei.osabi == ELF.ELFOSABI_NONE
            return Linux(arch)
        elseif oh.ei.osabi == ELF.ELFOSABI_FREEBSD
            return FreeBSD(arch)
        else
            error("Unknown ELF OSABI $(oh.ei.osabi)")
        end
    elseif oh isa MachOHandle
        return MacOS()
    elseif oh isa COFFHandle
        if is64bit(oh)
            return Windows(:x86_64)
        else
            return Windows(:i686)
        end
    else
        error("Unknown ObjectHandle type $(typeof(oh))")
    end
end

function rpaths(file::AbstractString)
    readmeta(file) do oh
        rpaths(RPath(oh))
    end
end

function canonical_rpaths(file::AbstractString)
    readmeta(file) do oh
        canonical_rpaths(RPath(oh))
    end
end

"""
    is_for_platform(h::ObjectHandle, platform::Platform)

Returns `true` if the given `ObjectHandle` refers to an object of the given
`platform`; E.g. if the given `platform` is for AArch64 Linux, then `h` must
be an `ELFHandle` with `h.header.e_machine` set to `ELF.EM_AARCH64`.

In particular, this method and `platform_for_object()` both exist because
the latter is not smart enough to deal with `:glibc` and `:musl` yet.
"""
function is_for_platform(h::ObjectHandle, platform::Platform)
    if platform isa Linux || platform isa FreeBSD
        # First off, if h isn't an ELF object, quit out
        if !(h isa ELFHandle)
            return false
        end
        # If the ELF object has an OSABI, check it matches platform
        if h.ei.osabi != ELF.ELF.ELFOSABI_NONE
            if platform isa Linux
                if h.ei.osabi != ELF.ELFOSABI_LINUX
                    return false
                end
            elseif platform isa FreeBSD
                if h.ei.osabi != ELF.ELFOSABI_FREEBSD
                    return false
                end
            else
                error("Unknown OS ABI type $(typeof(platform))")
            end
        end
        # Check that the ELF arch matches our own
        m = h.header.e_machine
        if platform.arch == :i686
            return m == ELF.EM_386
        elseif platform.arch == :x86_64
            # Allow i686 on x86_64, because that's technically ok
            return m == ELF.EM_386 || m == ELF.EM_X86_64
        elseif platform.arch == :aarch64
            return m == ELF.EM_AARCH64
        elseif platform.arch == :powerpc64le
            return m == ELF.EM_PPC64
        elseif platform.arch == :armv7l
            return m == ELF.EM_ARM
        else
            error("Unknown $(typeof(platform)) architecture $(platform.arch)")
        end
    elseif platform isa Windows
        if !(h isa COFFHandle)
            return false
        end

        if platform.arch == :x86_64
            return true
        elseif platform.arch == :i686
            return !is64bit(h)
        else
            error("Unknown $(typeof(platform)) architecture $(platform.arch)")
        end
    elseif platform isa MacOS
        # We'll take any old Mach-O handle
        if !(h isa MachOHandle)
            return false
        end
        return true
    else
        error("Unkown platform $(typeof(platform))")
    end
end

# These are libraries we should straight-up ignore, like libsystem on OSX
function should_ignore_lib(lib, ::ELFHandle)
    ignore_libs = [
        # Basic runtimes for both Linux and FreeBSD
        "libc.so",
        "libc.so.6",
        "libc.so.7",
        "libstdc++.so.6",
        "libdl.so.2",
        "librt.so.1",
        "libgcc_s.1.so",
        "libgcc_s.so.1",
        "libm.so.6",
        "libgfortran.so.3",
        "libgfortran.so.4",
        "libthr.so.3",
        # libpthread and libgomp are pretty safe bets
        "libpthread.so.0",
        "libgomp.so.1",
        # dynamic loaders from our alpine environment
        "ld-linux-x86-64.so.2",
    ]
    return lowercase(basename(lib)) in ignore_libs
end
function should_ignore_lib(lib, ::MachOHandle)
    ignore_libs = [
        "libsystem.b.dylib",
        "libstdc++.6.dylib",
        "libc++.1.dylib",
    ]
    return lowercase(basename(lib)) in ignore_libs
end
function should_ignore_lib(lib, ::COFFHandle)
    ignore_libs = [
        "msvcrt.dll",
        "kernel32.dll",
        "user32.dll",
        "libgcc_s_seh-1.dll",
        "libgcc_s_sjlj-1.dll",
        "libgfortran-3.dll",
        "libgfortran-4.dll",
        "libstdc++-6.dll",
    ]
    return lowercase(basename(lib)) in ignore_libs
end

# Determine whether a library is a "default" library or not, if it is we need
# to map it to `@rpath/$libname` on OSX.
is_default_lib(lib, oh) = false
function is_default_lib(lib, ::MachOHandle)
    default_libs = [
        "libgcc_s.1.dylib",
        "libgfortran.3.dylib",
        "libgfortran.4.dylib",
        "libquadmath.0.dylib",
    ]
    return lowercase(basename(lib)) in default_libs
end

function relink_to_rpath(prefix::Prefix, platform::Platform, path::AbstractString,
                         old_libpath::AbstractString; verbose::Bool = false)
    ur = preferred_runner()(prefix.path; cwd="/workspace/", platform=platform)
    rel_path = relpath(path, prefix.path)
    libname = basename(old_libpath)
    relink_cmd = ``

    if Compat.Sys.isapple(platform)
        install_name_tool = "/opt/x86_64-apple-darwin14/bin/install_name_tool"
        relink_cmd = `$install_name_tool -change $(old_libpath) @rpath/$(libname) $(rel_path)`
    elseif Compat.Sys.islinux(platform) || Compat.Sys.isbsd(platform)
        patchelf = "/usr/local/bin/patchelf"
        relink_cmd = `$patchelf --replace-needed $(old_libpath) $(libname) $(rel_path)`
    end

    # Create a new linkage that looks like @rpath/$lib on OSX,
    logpath = joinpath(logdir(prefix), "relink_to_rpath_$(basename(rel_path))_$(libname).log")
    run(ur, relink_cmd, logpath; verbose=verbose)
end

"""
    update_linkage(prefix::Prefix, platform::Platform, path::AbstractString,
                   old_libpath, new_libpath; verbose::Bool = false)

Given a binary object located at `path` within `prefix`, update its dynamic
linkage to point to `new_libpath` instead of `old_libpath`.  This is done using
a tool within the cross-compilation environment such as `install_name_tool` on
MacOS or `patchelf` on Linux.  Windows platforms are completely skipped, as
they do not encode paths or RPaths within their executables.
"""
function update_linkage(prefix::Prefix, platform::Platform, path::AbstractString,
                        old_libpath, new_libpath; verbose::Bool = false)
    # Windows doesn't do updating of linkage
    if Compat.Sys.iswindows(platform)
        return
    end

    ur = preferred_runner()(prefix.path; cwd="/workspace/", platform=platform)
    rel_path = relpath(path, prefix.path)

    add_rpath = x -> ``
    relink = (x, y) -> ``
    patchelf = "/usr/local/bin/patchelf"
    install_name_tool = "/opt/x86_64-apple-darwin14/bin/install_name_tool"
    if Compat.Sys.isapple(platform)
        add_rpath = rp -> `$install_name_tool -add_rpath @loader_path/$(rp) $(rel_path)`
        relink = (op, np) -> `$install_name_tool -change $(op) $(np) $(rel_path)`
    elseif Compat.Sys.islinux(platform) || Compat.Sys.isbsd(platform)
        current_rpaths = [r for r in rpaths(path) if !isempty(r)]
        add_rpath = rp -> begin
            # Join together RPaths to set new one
            rpaths = unique(vcat(current_rpaths, rp))

            # If any rpaths are `.`, map that to `$ORIGIN`
            remap_to_origin = path -> begin
                if path == "."
                    return "\$ORIGIN"
                end
                return path
            end
            rpaths = remap_to_origin.(rpaths)

            # I don't like strings ending in '/.', like '$ORIGIN/.'.  I don't think
            # it semantically makes a difference, but why not be correct AND beautiful?
            chomp_slashdot = path -> begin
                if length(path) > 2 && path[end-1:end] == "/."
                    return path[1:end-2]
                end
                return path
            end
            rpaths = chomp_slashdot.(rpaths)

            rpath_str = join(rpaths, ':')
            return `$patchelf --set-rpath $(rpath_str) $(rel_path)`
        end
        relink = (op, np) -> `$patchelf --replace-needed $(op) $(np) $(rel_path)`
    end

    # If the relative directory doesn't already exist within the RPATH of this
    # binary, then add it in.
    new_libdir = abspath(dirname(new_libpath) * "/")
    if !(new_libdir in canonical_rpaths(path))
        libname = basename(old_libpath)
        logpath = joinpath(logdir(prefix), "update_rpath_$(libname).log")
        cmd = add_rpath(relpath(new_libdir, dirname(path)))
        run(ur, cmd, logpath; verbose=verbose)
    end

    # Create a new linkage that depends on the RPATH to find things.
    # This allows us to split things up into multiple packages, and as long as the
    # libraries that this guy is interested in have been `dlopen()`'ed previously,
    # (and have the appropriate SONAME) things should "just work".
    logpath = joinpath(logdir(prefix), "update_linkage_$(basename(path))_$(basename(old_libpath)).log")
    new_libpath = relpath(new_libpath, dirname(path))
    if Compat.Sys.isapple(platform)
        # On MacOS, we need to explicitly add `@rpath/` before our library linkage path.
        new_libpath = joinpath("@rpath", new_libpath)
    end
    cmd = relink(old_libpath, new_libpath)
    run(ur, cmd, logpath; verbose=verbose)

    return new_libpath
end


text:
using Test, StaticArrays, Random
using LogisticRegression: optimize!

Random.seed!(0)
X = [rand(2) for _ in 1:10]
Y = [x[1] + x[2] > 1.0 ? 1.0 : 0.0 for x in X]
J = [0.0]
w = [1.0 -1.0]
b = [0.0]
optimize!(J, w, b, X, Y, num_iterations=50, learning_rate=3.0)
@test isapprox(J, [0.15912527602988893])
@test isapprox(w, [4.9046541360009215 4.340606046548434])
@test isapprox(b, [-3.867136029168926])

ext: jl


text:
module CircuitConstants

# Note: only import the constants that you want to use (to not pollute the namespace)
# For example:
#using CircuitConstants: T, G, M, k, m, u, n, p, f, a
#using CircuitConstants: V, mV, uV, nV, uV, pV
#using CircuitConstants: A, mA, uA, nA, uA, pA, fA, aA
#using CircuitConstants: s, ms, us, ns, us, ps, fs
#using CircuitConstants: GHz, MHz, kHz, Hz
#using CircuitConstants: TOhm, GOhm, MOhm, kOhm, Ohm, mOhm
#using CircuitConstants: TΩ, GΩ, MΩ, kΩ, Ω, mΩ

# Common SI prefixis:
const T = 1e12
const G = 1e9
const M = 1e6
const k = 1e3
const m = 1e-3
const u = 1e-6
const n = 1e-9
const p = 1e-12
const f = 1e-15
const a = 1e-18

# Other duplicates:
const μ = 1e-6

# Volts:
const TV = 1e12
const GV = 1e9
const MV = 1e6
const kV = 1e3
const mV = 1e-3
const V =  1
const uV = 1e-6
const μV = 1e-6
const nV = 1e-9
const pV = 1e-12
const fV = 1e-15
const aV = 1e-18

# Amps:
const TA = 1e12
const GA = 1e9
const MA = 1e6
const kA = 1e3
const A = 1
const mA = 1e-3
const uA = 1e-6
const μA = 1e-6
const nA = 1e-9
const pA = 1e-12
const fA = 1e-15
const aA = 1e-18

# Seconds:
const Ts = 1e12
const Gs = 1e9
const Ms = 1e6
const ks = 1e3
const s = 1
const ms = 1e-3
const us = 1e-6
const μs = 1e-6
const ns = 1e-9
const ps = 1e-12
const fs = 1e-15
const as = 1e-18

const Tsec = 1e12
const Gsec = 1e9
const Msec = 1e6
const ksec = 1e3
const sec = 1
const msec = 1e-3
const usec = 1e-6
const μsec = 1e-6
const nsec = 1e-9
const psec = 1e-12
const fsec = 1e-15
const asec = 1e-18

# Hertz:
const THz = 1e12
const GHz = 1e9
const MHz = 1e6
const kHz = 1e3
const Hz = 1
const mHz = 1e-3
const uHz = 1e-6
const μHz = 1e-6
const nHz = 1e-9
const pHz = 1e-12
const fHz = 1e-15
const aHz = 1e-18

# Farads:
const TF = 1e12
const GF = 1e9
const MF = 1e6
const kF = 1e3
const F = 1
const mF = 1e-3
const uF = 1e-6
const μF = 1e-6
const nF = 1e-9
const pF = 1e-12
const fF = 1e-15
const aF = 1e-18

# Henrys:
const TH = 1e12
const GH = 1e9
const MH = 1e6
const kH = 1e3
const H = 1
const mH = 1e-3
const uH = 1e-6
const μH = 1e-6
const nH = 1e-9
const pH = 1e-12
const fH = 1e-15
const aH = 1e-18

# Ohms:
const TOhms = 1e12
const GOhms = 1e9
const MOhms = 1e6
const kOhms = 1e3
const Ohms = 1
const mOhms = 1e-3
const uOhms = 1e-6
const μOhms = 1e-6
const nOhms = 1e-9
const pOhms = 1e-12
const fOhms = 1e-15
const aOhms = 1e-18

const TΩ = 1e12
const GΩ = 1e9
const MΩ = 1e6
const kΩ = 1e3
const Ω = 1
const mΩ = 1e-3
const uΩ = 1e-6
const μΩ = 1e-6
const nΩ = 1e-9
const pΩ = 1e-12
const fΩ = 1e-15
const aΩ = 1e-18

end # module

max_stars_repo_name: GlenHertz/CircuitConstants

max_issues_repo_name: GlenHertz/CircuitConstants

max_forks_repo_name: GlenHertz/CircuitConstants


text:
using Plots

"""Plot """
function canvas(; kwargs...)
    plot(; xticks = :none, yticks = :none, legend = false, axis = false, kwargs...)
end

"""Create unique color for all line types in the grid."""
function color_map(grid::Grid)
    t = LineTypes(grid)
    qs = collect(zip(t.x, t.y))
    g = cgrad(:inferno)
    zs = LinRange(0, 1, length(t))
    return Dict(d => RGB(g[z]) for (d, z) in zip(qs, zs))
end

function plot_grid!(plt, grid::Grid; kwargs...)
    scatter!(plt, grid.x, grid.y; kwargs...)
end

function plot_lines!(plt, lines::Vector{NTuple{2,Vector{Int}}}, colors::Vector{RGB}; kwargs...)
    for ((x, y), color) in zip(lines, colors)
        plot!(plt, x, y; linecolor=color, kwargs...)
    end
end

function plot_lock_pattern!(plt, grid::Grid, p::Vector{Int})
    c_map = color_map(grid)
    lines = NTuple{2,Vector{Int}}[]
    colors = RGB[]
    for i = 1:length(p)-1
        j = p[[i, i + 1]]
        x = grid.x[j]
        y = grid.y[j]
        t = line_type(x[1], x[2], y[1], y[2])
        push!(lines, (x, y))
        push!(colors, c_map[t])
    end
    plot_lines!(plt, lines, colors; linewidth = 3)
    plot_grid!(plt, grid; marker = :circle, color = :black, markersize = 5)
    return plt
end

function create_plots(n::Int, path::AbstractString, directory::AbstractString; limit::Int=0)
    grid = Grid(n)
    solutions = extract_solutions(path)
    for (distance, patterns) in solutions
        if distance ≥ limit
            root = mkpath(joinpath(directory, "$(join(grid.n, "x"))", "$(distance)"))
            for pattern in patterns
                plt = canvas(; size=(400, 400), backgroud=:lightgray)
                plot_lock_pattern!(plt, grid, pattern)
                savefig(plt, joinpath(root, "$(hash(pattern)).svg"))
            end
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

function prefix(pts, i)
    inds = eachindex(pts)
    index = first(inds) : i
    view(pts, index)
end

function move_to_front!(pts, i)
    @assert i in eachindex(pts)
    pt = pts[i]
    for j in eachindex(pts)
        qt = pts[j]
        pts[j] = pt
        pt = qt
        j == i && break
    end
    pts
end

function leq_approx(x,y;kw...)
    x < y || isapprox(x,y;kw...)
end


text:
for n in names(Actions, all=true)
    v = eval(Actions,n)
    if typeof(v) == Action

        docont = false
        for i in names(GdkKeySyms, all=true)
            vk = eval(GdkKeySyms,i)
            if typeof(vk) <: Integer
                if v.keyval == vk
                    docont = true
                    break
                end
            end
        end
        docont && continue

        k = gdk_keyval_name(v.keyval)
        s = ""
        if v.state == PrimaryModifier
            s = "Ctrl"
        elseif v.state == PrimaryModifier + GdkModifierType.SHIFT
            s = "Ctrl+Shift"
        elseif v.state == SecondaryModifer
            s = "Alt"
        end
        if s != ""
            d = v.description
            println( "- `$s+$k` $d.")
        end
    end
end


text:
using Qsosed
using Test

include("parameters_tests.jl")
include("black_hole_tests.jl")
include("disk_tests.jl")
include("corona_tests.jl")
include("warm_tests.jl")
include("xspec_routines_tests.jl")
include("fraction_tests.jl")
include("compton_tests.jl")


text:
using Distributions
type SmallFiniteHorizonMDP
    P::Array{Float64, 4}
    atbls::Array{Distributions.AliasTable,3}
    R::Array{Float64, 3}
    p0::Categorical
    s::UInt64
    t::UInt64
end
function SmallFiniteHorizonMDP(P, R, p0, s, t)
    S, A, H = size(R)
    atbls = Array{Distributions.AliasTable, 3}(S, A, H)
    for t=1:H, a=1:A, s=1:S
        atbls[s,a,t] = Distributions.AliasTable(P[:, s, a, t])
    end
    SmallFiniteHorizonMDP(P, atbls, R, p0, s, t)
end

horizon(mdp::SmallFiniteHorizonMDP) = size(mdp.R, 3)
nS(mdp::SmallFiniteHorizonMDP) = size(mdp.R, 1)
nA(mdp::SmallFiniteHorizonMDP) = size(mdp.R, 2)
rewards(mdp::SmallFiniteHorizonMDP) = mdp.R
function restart!(mdp::SmallFiniteHorizonMDP)
    s = rand(mdp.p0)
    mdp.s = s
    mdp.t = 1
    s
end
function step!(mdp::SmallFiniteHorizonMDP, a)
    r = mdp.R[mdp.s, a, mdp.t]
    sn = rand(mdp.atbls[mdp.s, a, mdp.t])
    mdp.s = sn
    mdp.t += 1
    stop = mdp.t > horizon(mdp)
    sn, r, stop
end

function expected_return(mdp::SmallFiniteHorizonMDP, policy)
    H = horizon(mdp)
    S = nS(mdp)
    A = nA(mdp)

    Vold = zeros(S)
    Vnew = zeros(S)

    for t=H:-1:1
        for s=1:S
            a = policy(s, t)
            Vnew[s] = mdp.R[s, a, t] + mdp.P[:, s, a, t] ⋅ Vold
        end
        Vnew, Vold = Vold, Vnew
    end
    mdp.p0.p ⋅ Vold
end

function optimal_policy(mdp::SmallFiniteHorizonMDP)
    H = horizon(mdp)
    S = nS(mdp)
    A = nA(mdp)
    policy = zeros(UInt64, S, H)
    Vold = zeros(S)
    Qnew = zeros(S, A)

    for t=H:-1:1
        Qnew = mdp.R[:, :, t] + squeeze(sum(Vold .* mdp.P[:, :, :, t], 1), 1)
        for s=1:S
            policy[s, t] = indmax(Qnew[s, :])
            Vold[s] = Qnew[s, policy[s, t]]
        end
    end
    (s, t) -> policy[s, t]
end

function value_fun(mdp::SmallFiniteHorizonMDP, policy)
    H = horizon(mdp)
    S = nS(mdp)
    A = nA(mdp)

    V = zeros(S, H+1)

    for t=H:-1:1
        for s=1:S
            a = policy(s, t)
            V[s, t] = mdp.R[s, a, t] + mdp.P[:, s, a, t] ⋅ V[:, t+1]
        end
    end
    V
end

gammarand(a) = if a >= 1
            rand(Distributions.GammaGDSampler(Gamma(a)))
        else
            rand(Distributions.GammaGSSampler(Gamma(a)))
        end

function dirrand!(res, αs)
    ss = 0.
    for (i, α) in enumerate(αs)
        b = gammarand(α)
        ss += b
        res[i] = b
    end
    res ./= ss
end
function dirrand(αs)
    res = zeros(length(αs))
    dirrand!(res, αs)
    res
end
function randomMDP(S, A, H, α =.1)
    dir = Dirichlet(S, α)
    P = zeros(S, S, A, H)
    R = rand(S, A, H)
    R[rand(S, A, H) .<= 0.85] = 0.
    p0 = Categorical(dirrand(repeated(α, S)))
    for s=1:S, a=1:A, t=1:H
        P[:, s, a, t] = dirrand(repeated(α, S))
    end
    SmallFiniteHorizonMDP(P, R, p0, 1, 1)
end
function stateDepTestMDP(N, ϵ=1e-2)
    S = 2N + 2
    A = 2
    H = 2
    pp0 = zeros(S)
    pp0[1] = 1
    p0 = Categorical(pp0)
    R = zeros(S, A, H)
    R[(S - N + 1):end, :, 2] = 1
    R[2, :, 2] = 0.5 + ϵ
    P = zeros(S, S, A, H)
    for s=1:S, a=1:A
        P[s, s, a, 2] = 1.
    end
    P[3:end, 1, 1, 1] = 1 / 2 / N
    P[2, 1, 2, 1] = 1.
    SmallFiniteHorizonMDP(P, R, p0, 1, 1)
end

type ChainMDP
    P::Array{Float64, 3}
    R::Vector{Float64}
    N::Int64
    s::UInt64
    t::UInt64
    r0::Float64
end
nS(m::ChainMDP) = m.N
nA(m::ChainMDP) = 2
horizon(m::ChainMDP) = m.N
function rewards(m::ChainMDP)
    R = zeros(nS(m), nA(m), horizon(m))
    R .= m.R
    R
end

function ChainMDP(N, r0=0.)
    P = zeros(N, N, 2)
    for s=1:N
        sn = max(s - 1, 1)
        P[sn, s, 1] = 1.
        P[sn, s, 2] = 1./N
        sn = min(s+1, N)
        P[sn, s, 2] = 1 - 1./N
    end

    R = zeros(N)
    R[end] = 1.
    R[1] = r0
    ChainMDP(P, R, N, 1, 1, r0)
end
function restart!(mdp::ChainMDP)
    mdp.s = 1
    mdp.t = 1
    1
end
function step!(mdp::ChainMDP, a)
    r = 0.
    if mdp.s == mdp.N
        r = 1.
    elseif mdp.s == 1
        r = mdp.r0
    end
    sn = max(mdp.s-1, 1)
    if a == 2 && rand() > 1/mdp.N
        # go right
        sn = min(mdp.s+1,mdp.N)
    end
    mdp.s = sn
    mdp.t += 1
    stop = mdp.t > horizon(mdp)
    sn, r, stop
end

function expected_return(mdp::ChainMDP, policy)
    H = horizon(mdp)
    S = nS(mdp)
    A = nA(mdp)

    Vold = zeros(S)
    Vnew = zeros(S)

    for t=H:-1:1
        for s=1:S
            a = policy(s, t)
            Vnew[s] = mdp.R[s] + mdp.P[:, s, a] ⋅ Vold
        end
        Vnew, Vold = Vold, Vnew
    end
    Vold[1]
end

function optimal_policy(mdp::ChainMDP)
    H = horizon(mdp)
    S = nS(mdp)
    A = nA(mdp)
    policy = zeros(UInt64, S, H)
    Vold = zeros(S)
    Qnew = zeros(S, A)

    for t=H:-1:1
        Qnew = reshape(mdp.R, S, 1) .+ squeeze(sum(Vold .* mdp.P, 1), 1)
        for s=1:S
            policy[s, t] = indmax(Qnew[s, :])
            Vold[s] = Qnew[s, policy[s, t]]
        end
    end
    (s, t) -> policy[s, t]
end


text:
using InterpretableModels: getfactor, parambounds, getintercept, featurename
using InterpretableModels: setfeaturenames!, modelasstring, setparams!, hasparams
using InterpretableModels: minlag, actualminlag, maxlag, actualmaxlag, params, horizon
using InterpretableModels: predict, predict!, fit!
using DataFrames

@testset "LaggedFactorLM" begin

lm = LaggedFactorLITRM(2, 4;
        minlag = 3, maxlag = 8,
        hasintercept = true, mincoef = -10.0, maxcoef = 97.5)

@test !isclassificationmodel(lm)
@test  isregressionmodel(lm)
@test  istimeseriesmodel(lm)
@test !hasparams(lm)
@test minlag(lm) == 3
@test actualminlag(lm) == 3
@test horizon(lm) == 3
@test maxlag(lm) == 8
@test actualmaxlag(lm) == 8

setparams!(lm, [1.0, 2.24, 5.5, 4.1, 1.2, 8.0, 0.5])
@test hasparams(lm)
@test minlag(lm) == 3
@test actualminlag(lm) == 5 # Since 5.5 < 8.0
@test actualmaxlag(lm) == 8 # Since 8.0 > 5.5
@test horizon(lm) == 5
@test params(lm) == [1.0, 2.24, 5.5, 4.1, 1.2, 8.0, 0.5]

c, f, l = getfactor(lm, 1)
@test c == 1.0
@test isa(f, Int)
@test f == 2
@test isa(l, Int)
@test l == 5

c, f, l = getfactor(lm, 2)
@test c == 4.1
@test isa(f, Int)
@test f == 1
@test isa(l, Int)
@test l == 8

@test getintercept(lm) == 0.5

@testset "parambounds" begin
    bs = parambounds(lm)
    @test length(bs) == (3*2+1) # 3 values per factor and 1 for intercept
    for fi in 1:2
        o = (fi-1)*3
        @test bs[1+o] == (-10.0, 97.5)
        @test bs[2+o][1] == 1.0  # feature 1 is lowest possible feature number
        @test isapprox(bs[2+o][2], 1.0+4.0, atol=1e-4)  # Should almost be 1+nfeatures
        @test bs[3+o][1] == 3.0                         # minlag
        @test isapprox(bs[3+o][2], 1.0+8.0, atol=1e-4)  # maxlag
    end
    @test bs[end] == (-10.0, 97.5)
end # @testset "parambounds" begin

@testset "featurenames and setfeaturenames!" begin
    @test featurename(lm, 1) == "feature_1"
    @test featurename(lm, 2) == "feature_2"
    setfeaturenames!(lm, ["a", "bfg"])
    @test featurename(lm, 1) == "a"
    @test featurename(lm, 2) == "bfg"
end # @testset "featurenames and setfeaturenames!" begin

@testset "modelasstring" begin
    setparams!(lm, [1.2, 1.6, 5.8, 6.987, 2.3, 7.3, 0.762])
    @test modelasstring(lm) == "0.76 + 1.2 * a[T-5] + 6.99 * bfg[T-7]"
    df = DataFrame(MyFactor1 = Int[], AnotherFac = Float64[])
    @test modelasstring(lm, df) == "0.76 + 1.2 * MyFactor1[T-5] + 6.99 * AnotherFac[T-7]"
end

@testset "predict, predict! and fit!" begin
    lm = LaggedFactorLITRM(2, 3;
            minlag = 1, maxlag = 3,
            hasintercept = true, mincoef = -10.0, maxcoef = 97.5)

    df = DataFrame(
        FacA = Int[    1,    2,   3,    4,    5],
        Fac2 = Float64[3.14, 5.4, 9.8, 10.0, 25.0],
        Fac3 = Float64[4.00, 6.0, 8.0, 10.0, 12.0])

    X = Array{Float64}(df)

    # Since we have a lag range of 1-2 we can only get three predictions
    # from this model, i.e. for rows 3-5.
    # True expression: y = 1.34 .+ X[T-1, 1] .* (-2.5) .+ X[T-2, 3] .* 4.0
    trueparams = [-2.5, 1.0, 1.0, 4.0, 3.0, 2.0, 1.34]
    y = [Inf, Inf,
            1.34 + (-2.5)*X[3-1, 1] + 4.0*X[3-2, 3],
            1.34 + (-2.5)*X[4-1, 1] + 4.0*X[4-2, 3],
            1.34 + (-2.5)*X[5-1, 1] + 4.0*X[5-2, 3]]
    yhat = predict(lm, X, trueparams)
    @test length(yhat) == 3
    @test yhat[1] == y[3]
    @test yhat[2] == y[4]
    @test yhat[3] == y[5]

    yhatpreds = zeros(Float64, 3)
    predict!(lm, yhatpreds, X, trueparams)
    @test yhatpreds[1] == y[3]
    @test yhatpreds[2] == y[4]
    @test yhatpreds[3] == y[5]

    # If we don't have the true params we can fit the model to data to recover them.
    dummyparams = [3.5, 1.5, 1.5, 2.1, 2.0, 1.7, 1.34]
    setparams!(lm, dummyparams) # Set a random param vector to ensure a better one is set after fitting
    fit!(lm, X, y; MaxTime = 3.0)
    @test length(params(lm)) == (1+2*3)
    @test params(lm) != dummyparams

    # Since the problem is underconstrained we cannot expect to find the "true"
    # params since there are many that works...
    yhat = predict(lm, X)
    for i in 1:length(yhat)
        @test isapprox(yhat[i], y[end-(length(yhat)-i)], atol=1e-1)
    end
end

end # @testset "LaggedFactorLM" begin


text:
using Agents
using Agents.CA2D

# 0. Define the rule
# Rules of Conway's game of life: DSR (Death, Survival, Reproduction). Cells die if the number of their living neighbors are $<D$, survive if the number of their living neighbors are $<=S$, come to life if their living neighbors are as many as $R$.
rules = (2,3,3)

# 1. Build the model
model = CA2D.build_model(rules=rules, dims=(100, 100), Moore=true)  # creates a model where all cells are "0"
# make some random cells alive
for i in 1:gridsize(model.space.dimensions)
  if rand() < 0.1
    model.agents[i].status="1"
  end
end

# 2. Run the model, collect data, and visualize it
runs = 50
CA2D.ca_run(model, runs)

# convert -delay 10 -loop 1 *.png animation.gif


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using EzXML

"""
Extract the FMU variables and meta data from the ModelDescription
"""
function fmi2ReadModelDescription(pathToModellDescription::String)
    md = fmi2ModelDescription()

    md.stringValueReferences = Dict{String, fmi2ValueReference}()
    md.outputValueReferences = Array{fmi2ValueReference}(undef, 0)
    md.inputValueReferences = Array{fmi2ValueReference}(undef, 0)
    md.stateValueReferences = Array{fmi2ValueReference}(undef, 0)
    md.derivativeValueReferences = Array{fmi2ValueReference}(undef, 0)

    # CS specific entries
    md.CSmodelIdentifier = ""
    md.CScanHandleVariableCommunicationStepSize = false
    md.CScanInterpolateInputs = false
    md.CSmaxOutputDerivativeOrder = -1
    md.CScanGetAndSetFMUstate = false
    md.CScanSerializeFMUstate = false
    md.CSprovidesDirectionalDerivative = false

    # ME specific entries
    md.MEmodelIdentifier = ""
    md.MEcanGetAndSetFMUstate = false
    md.MEcanSerializeFMUstate = false
    md.MEprovidesDirectionalDerivative = false

    md.enumerations = []
    typedefinitions = nothing
    modelvariables = nothing
    modelstructure = nothing

    doc = readxml(pathToModellDescription)

    root = doc.root

    md.fmiVersion = root["fmiVersion"]
    md.modelName = root["modelName"]
    md.guid = root["guid"]
    md.isModelExchange = false
    md.isCoSimulation = false
    md.generationTool = parseNodeString(root, "generationTool"; onfail="[Unknown generation tool]")
    md.generationDateAndTime = parseNodeString(root, "generationDateAndTime"; onfail="[Unknown generation date and time]")
    md.variableNamingConvention = parseNodeString(root, "variableNamingConvention"; onfail="[Unknown variable naming convention]")
    md.numberOfEventIndicators = parseNodeInteger(root, "numberOfEventIndicators"; onfail=0)
    md.description = parseNodeString(root, "description"; onfail="[Unknown Description]")

    for node in eachelement(root)

        if node.name == "CoSimulation"
            md.isCoSimulation = true
            md.CSmodelIdentifier                        = node["modelIdentifier"]
            md.CScanHandleVariableCommunicationStepSize = parseNodeBoolean(node, "canHandleVariableCommunicationStepSize"   ; onfail=false)
            md.CScanInterpolateInputs                   = parseNodeBoolean(node, "canInterpolateInputs"                     ; onfail=false)
            md.CSmaxOutputDerivativeOrder               = parseNodeInteger(node, "maxOutputDerivativeOrder"                 ; onfail=-1)
            md.CScanGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUstate"                     ; onfail=false)
            md.CScanSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUstate"                     ; onfail=false)
            md.CSprovidesDirectionalDerivative          = parseNodeBoolean(node, "providesDirectionalDerivative"            ; onfail=false)

        elseif node.name == "ModelExchange"
            md.isModelExchange = true
            md.MEmodelIdentifier                        = node["modelIdentifier"]
            md.MEcanGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUstate"                     ; onfail=false)
            md.MEcanSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUstate"                     ; onfail=false)
            md.MEprovidesDirectionalDerivative          = parseNodeBoolean(node, "providesDirectionalDerivative"            ; onfail=false)

        elseif node.name == "TypeDefinitions"
            typedefinitions = node

        elseif node.name == "ModelVariables"
            modelvariables = node

        elseif node.name == "ModelStructure"
            modelstructure = node
        end
    end

    if typedefinitions == nothing
        md.enumerations = []
    else
        md.enumerations = createEnum(typedefinitions)
    end

    md.valueReferences = []
    md.valueReferenceIndicies = Dict{Integer,Integer}()

    derivativeindices = getDerivativeIndices(modelstructure)
    md.modelVariables = parseModelVariables(modelvariables, md, derivativeindices)

    # parse model dependencies (if available)
    for element in eachelement(modelstructure)
        if element.name == "Derivatives" || element.name == "InitialUnknowns"
            parseDependencies(element, md)
        elseif element.name == "Outputs"

        else
            @warn "Unknown tag `$(element.name)` for node `ModelStructure`."
        end
    end

    # creating an index for value references (fast look-up for dependencies)
    for i in 1:length(md.valueReferences)
        md.valueReferenceIndicies[md.valueReferences[i]] = i
    end

    md
end

"""
Returns the tag 'modelName' from the model description.
"""
function fmi2GetModelName(md::fmi2ModelDescription)#, escape::Bool = true)
    md.modelName
end

"""
Returns the tag 'guid' from the model description.
"""
function fmi2GetGUID(md::fmi2ModelDescription)
    md.guid
end

"""
Returns the tag 'generationtool' from the model description.
"""
function fmi2GetGenerationTool(md::fmi2ModelDescription)
    md.generationTool
end

"""
Returns the tag 'generationdateandtime' from the model description.
"""
function fmi2GetGenerationDateAndTime(md::fmi2ModelDescription)
    md.generationDateAndTime
end

"""
Returns the tag 'varaiblenamingconvention' from the model description.
"""
function fmi2GetVariableNamingConvention(md::fmi2ModelDescription)
    md.variableNamingConvention
end

"""
Returns the tag 'numberOfEventIndicators' from the model description.
"""
function fmi2GetNumberOfEventIndicators(md::fmi2ModelDescription)
    md.numberOfEventIndicators
end

"""
Returns if the FMU model description contains `dependency` information.
"""
function fmi2DependenciesSupported(md::fmi2ModelDescription)
    for mv in md.modelVariables
        if mv.dependencies != nothing && length(mv.dependencies) > 0
            return true
        end
    end

    return false
end

"""
Returns the tag 'modelIdentifier' from CS or ME section.
"""
function fmi2GetModelIdentifier(md::fmi2ModelDescription)
    if fmi2IsCoSimulation(md)
        return md.CSmodelIdentifier
    elseif fmi2IsModelExchange(md)
        return md.MEmodelIdentifier
    else
        @assert false "fmi2GetModelName(...): FMU does not support ME or CS!"
    end
end

"""
Returns true, if the FMU supports the getting/setting of states
"""
function fmi2CanGetSetState(md::fmi2ModelDescription)
    if md.CScanGetAndSetFMUstate || md.MEcanGetAndSetFMUstate
        return true
    else
        return false
    end
end

"""
Returns true, if the FMU state can be serialized
"""
function fmi2CanSerializeFMUstate(md::fmi2ModelDescription)
    if md.CScanSerializeFMUstate || md.MEcanSerializeFMUstate
        return true
    else
        return false
    end
end

"""
Returns true, if the FMU provides directional derivatives
"""
function fmi2ProvidesDirectionalDerivative(md::fmi2ModelDescription)
    if md.CSprovidesDirectionalDerivative || md.MEprovidesDirectionalDerivative
        return true
    else
        return false
    end
end

"""
Returns true, if the FMU supports co simulation
"""
function fmi2IsCoSimulation(md::fmi2ModelDescription)
    md.isCoSimulation
end

"""
Returns true, if the FMU supports model exchange
"""
function fmi2IsModelExchange(md::fmi2ModelDescription)
    md.isModelExchange
end

# Returns the indices of the state derivatives.
function getDerivativeIndices(node::EzXML.Node)
    indices = []
    for element in eachelement(node)
        if element.name == "Derivatives"
            for derivative in eachelement(element)
                ind = parse(Int, derivative["index"])
                der = nothing
                derKind = nothing

                if haskey(derivative, "dependencies")
                    der = split(derivative["dependencies"], " ")

                    if der[1] == ""
                        der = fmi2Integer[]
                    else
                        der = collect(parse(fmi2Integer, e) for e in der)
                    end
                end

                if haskey(derivative, "dependenciesKind")
                    derKind = split(derivative["dependenciesKind"], " ")
                end

                push!(indices, (ind, der, derKind))
            end
        end
    end
    sort!(indices, rev=true)
end

# Parses the model variables of the FMU model description.
function parseModelVariables(nodes::EzXML.Node, md::fmi2ModelDescription, derivativeIndices)
    lastValueReference = fmi2ValueReference(0)
    derivativeIndex = nothing
    if derivativeIndices != []
        derivativeIndex = pop!(derivativeIndices)
    end
    numberOfVariables = 0
    for node in eachelement(nodes)
        numberOfVariables += 1
    end
    scalarVariables = Array{fmi2ScalarVariable}(undef, numberOfVariables)
    index = 1

    for node in eachelement(nodes)
        name = node["name"]
        ValueReference = parse(fmi2ValueReference, (node["valueReference"]))
        description = ""
        causality = ""
        variability = ""
        initial = ""

        if !(ValueReference in md.valueReferences)
            push!(md.valueReferences, ValueReference)
        end

        if haskey(node, "description")
            description = node["description"]
        end
        if haskey(node, "causality")
            causality = node["causality"]
        end
        if haskey(node, "variability")
            variability = node["variability"]
        end
        if haskey(node, "initial")
            initial = node["initial"]
        end
        datatype = setDatatypeVariables(node, md)

        dependencies = []
        dependenciesKind = []

        if derivativeIndex != nothing
            if index == derivativeIndex[1]
                push!(md.stateValueReferences, lastValueReference)
                push!(md.derivativeValueReferences, ValueReference)

                if derivativeIndices != []
                    derivativeIndex = pop!(derivativeIndices)
                end
            end
        end

        scalarVariables[index] = fmi2ScalarVariable(name, ValueReference, datatype, description, causality, variability, initial, dependencies, dependenciesKind)

        if causality == "output"
            push!(md.outputValueReferences, ValueReference)
        elseif causality == "input"
            push!(md.inputValueReferences, ValueReference)
        end
        md.stringValueReferences[name] = ValueReference

        lastValueReference = ValueReference
        index += 1
    end
    md.numberOfContinuousStates = length(md.stateValueReferences)
    scalarVariables
end

# Parses the model variables of the FMU model description.
function parseDependencies(nodes::EzXML.Node, md::fmi2ModelDescription)
    for node in eachelement(nodes)

        if node.name == "Unknown"

            index = 0
            dependencies = nothing
            dependenciesKind = nothing

            if haskey(node, "index") && haskey(node, "dependencies") && haskey(node, "dependenciesKind")
                index = parseInteger(node["index"])
                dependencies = node["dependencies"]
                dependenciesKind = node["dependenciesKind"]

                if length(dependencies) > 0 && length(dependenciesKind) > 0
                    dependenciesSplit = split(dependencies, " ")
                    dependenciesKindSplit = split(dependenciesKind, " ")

                    if length(dependenciesSplit) != length(dependenciesKindSplit)
                        @warn "Length of field dependencies ($(length(dependenciesSplit))) doesn't match length of dependenciesKind ($(length(dependenciesKindSplit)))."
                    else
                        md.modelVariables[index].dependencies = vcat(md.modelVariables[index].dependencies, collect(parseInteger(s) for s in dependenciesSplit))
                        md.modelVariables[index].dependenciesKind = vcat(md.modelVariables[index].dependenciesKind,  dependenciesKindSplit)
                    end
                else
                    md.modelVariables[index].dependencies = []
                    md.modelVariables[index].dependenciesKind = []
                end
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`."
            end
        else
            @warn "Unknown entry in `ModelStructure` named `$(node.name)`."
        end
    end
end

"""
Returns the model variable(s) fitting the value reference.
"""
function fmi2ModelVariablesForValueReference(md::fmi2ModelDescription, vr::fmi2ValueReference)
    ar = []
    for modelVariable in md.modelVariables
        if modelVariable.valueReference == vr
            push!(ar, modelVariable)
        end
    end
    ar
end

# Parses a Bool value represented by a string.
function parseBoolean(s::Union{String, SubString{String}}; onfail=nothing)
    if s == "true"
        return true
    elseif s == "false"
        return false
    else
        @assert onfail != nothing ["parseBoolean(...) unknown boolean value '$s'."]
        return onfail
    end
end

function parseNodeBoolean(node, key; onfail=nothing)
    if haskey(node, key)
        return parseBoolean(node[key]; onfail=onfail)
    else
        return onfail
    end
end

# Parses an Integer value represented by a string.
function parseInteger(s::Union{String, SubString{String}}; onfail=nothing)
    if onfail == nothing
        return parse(Int, s)
    else
        try
            return parse(Int, s)
        catch
            return onfail
        end
    end
end

function parseNodeInteger(node, key; onfail=nothing)
    if haskey(node, key)
        return parseInteger(node[key]; onfail=onfail)
    else
        return onfail
    end
end

function parseNodeString(node, key; onfail=nothing)
    if haskey(node, key)
        return node[key]
    else
        return onfail
    end
end

# Parses a fmi2Boolean value represented by a string.
function parseFMI2Boolean(s::Union{String, SubString{String}})
    if parseBoolean(s)
        return fmi2True
    else
        return fmi2False
    end
end

# set the datatype and attributes of an model variable
function setDatatypeVariables(node::EzXML.Node, md::fmi2ModelDescription)
    type = datatypeVariable()
    typenode = node.firstelement
    typename = typenode.name
    type.start = nothing
    type.min = nothing
    type.max = nothing
    type.quantity = nothing
    type.unit = nothing
    type.displayUnit = nothing
    type.relativeQuantity = nothing
    type.nominal = nothing
    type.unbounded = nothing
    type.derivative = nothing
    type.reinit = nothing
    if typename == "Real"
        type.datatype = fmi2Real
    elseif typename == "String"
        type.datatype = fmi2String
    elseif typename == "Boolean"
        type.datatype = fmi2Boolean
    elseif typename == "Integer"
        type.datatype = fmi2Integer
    elseif typename == "Enumeration"
        type.datatype = fmi2Enum
    else
        @warn "Unknown data type `$(type.datatype)`."
    end

    if haskey(typenode, "declaredType")
        type.declaredType = typenode["declaredType"]
    end

    if haskey(typenode, "start")
        if typename == "Real"
            type.start = parse(fmi2Real, typenode["start"])
        elseif typename == "Integer"
            type.start = parse(fmi2Integer, typenode["start"])
        elseif typename == "Boolean"
            type.start = parseFMI2Boolean(typenode["start"])
        elseif typename == "Enumeration"
            for i in 1:length(md.enumerations)
                if type.declaredType == md.enumerations[i][1] # identify the enum by the name
                    type.start = md.enumerations[i][1 + parse(Int, typenode["start"])] # find the enum value and set it
                end
            end
        elseif typename == "String"
            type.start = typenode["start"]
        else
            @warn "setDatatypeVariables(...) unimplemented start value type $typename"
            type.start = typenode["start"]
        end
    end

    if haskey(typenode, "min") && (type.datatype == fmi2Real || type.datatype == fmi2Integer || type.datatype == fmi2Enum)
        if type.datatype == fmi2Real
            type.min = parse(fmi2Real, typenode["min"])
        else
            type.min = parse(fmi2Integer, typenode["min"])
        end
    end
    if haskey(typenode, "max") && (type.datatype == fmi2Real || type.datatype == fmi2Integer || type.datatype == fmi2Enum)
        if type.datatype == fmi2Real
            type.max = parse(fmi2Real, typenode["max"])
        elseif type.datatype == fmi2Integer
            type.max = parse(fmi2Integer, typenode["max"])
        end
    end
    if haskey(typenode, "quantity") && (type.datatype == fmi2Real || type.datatype == fmi2Integer || type.datatype == fmi2Enum)
        type.quantity = typenode["quantity"]
    end
    if haskey(typenode, "unit") && type.datatype == fmi2Real
        type.unit = typenode["unit"]
    end
    if haskey(typenode, "displayUnit") && type.datatype == fmi2Real
        type.displayUnit = typenode["displayUnit"]
    end
    if haskey(typenode, "relativeQuantity") && type.datatype == fmi2Real
        type.relativeQuantity = convert(fmi2Boolean, parse(Bool, typenode["relativeQuantity"]))
    end
    if haskey(typenode, "nominal") && type.datatype == fmi2Real
        type.nominal = parse(fmi2Real, typenode["nominal"])
    end
    if haskey(typenode, "unbounded") && type.datatype == fmi2Real
        type.unbounded = parse(fmi2Boolean, typenode["unbounded"])
    end
    if haskey(typenode, "derivative") && type.datatype == fmi2Real
        type.derivative = parse(fmi2Integer, typenode["derivative"])
    end
    if haskey(typenode, "reinit") && type.datatype == fmi2Real
        type.reinit = parseFMI2Boolean(typenode["reinit"])
    end
    type
end

#=
Read all enumerations from the modeldescription and store them in a matrix. First entries are the enum names
-------------------------------------------
Example:
"enum1name" "value1"    "value2"
"enum2name" "value1"    "value2"
=#
function createEnum(node::EzXML.Node)
    enum = 1
    idx = 1
    enumerations = []
    for simpleType in eachelement(node)
        name = simpleType["name"]
        for type in eachelement(simpleType)
            if type.name == "Enumeration"
                enum = []
                push!(enum, name)
                for item in eachelement(type)
                    push!(enum, item["name"])
                end
                push!(enumerations, enum)
            end
        end
    end
    enumerations
end


text:
module book

end

max_stars_repo_name: johnfoster-pge-utexas/hpc-book

max_issues_repo_name: johnfoster-pge-utexas/hpc-book

max_forks_repo_name: johnfoster-pge-utexas/hpc-book


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
nq = s.model.dim.q
nu = s.model.dim.u
nc = s.model.dim.c
nb = nc * friction_dim(s.env)
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd

# get trajectory
ref_traj = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_in_place.jld2"),
    load_type=:joint_traj)

H = ref_traj.H
h = ref_traj.h
κ = 1.0e-8

# Cost function
obj = TrackingObjective(s.model, s.env, H,
    q = [[Diagonal(1.0e-10 * ones(nq)) for t = 1:H-1]; [Diagonal(1.0e9 * ones(nq))]],
    u = [Diagonal(1.0e-5 * [0.1, 1.0]) for t = 1:H],
    γ = [Diagonal(1.0e-100 * ones(nc)) for t = 1:H],
    b = [Diagonal(1.0e-100 * ones(nb)) for t = 1:H])

n_opts = NewtonOptions(r_tol=3e-8, max_iter=100, verbose = true)

im_traj_ = ImplicitTraj(ref_traj, s, κ=κ)

core = Newton(s, H, h, ref_traj, im_traj_, obj = obj, opts = n_opts)

q0_dist = deepcopy(ref_traj.q[1] + [-0.1,0.0,0,0])
q1_dist = deepcopy(ref_traj.q[2] + [-0.1,0.0,0,0])
newton_solve!(core, s, im_traj_, ref_traj, q0=q0_dist, q1=q1_dist)

l = 1
lu = 1
plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.([q[1:end] for q in core.traj.q])...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][lu:lu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[3,1], hcat(Vector.(vcat([fill(ref_traj.γ[i][1:nc], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[1:end] for u in core.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)

visualize_robot!(vis, s.model, core.traj.q)

plot(hcat(Vector.(ref_traj.u)...)')
plot(hcat(Vector.(core.traj.u)...)')
# Check ~perfect loop
@show round.(core.traj.q[end-2] - ref_traj.q[1], digits=3)
@show round.(core.traj.q[end-1] - ref_traj.q[2], digits=3)

# Save trajectory
traj = deepcopy(core.traj)
gait_path = joinpath(@__DIR__, "..", "dynamics", "hopper_2D", "gaits", "gait_forward.jld2")
@save gait_path traj

# Reload trajectory
res = JLD2.jldopen(gait_path)
loaded_traj = res["traj"]

traj = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_forward.jld2"),
    load_type=:joint_traj)


text:
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of Intel Corporation nor the names of its contributors
#       may be used to endorse or promote products derived from this software
#       without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# FIXME: This should be moved into the Latte core as it uses special internal
# handling.
export ReshapeLayer

function ReshapeEnsemble(name::Symbol, input::AbstractEnsemble, shape)
    neurons = Array(ReshapeNeuron, shape...)
    for i in 1:length(neurons)
        neurons[i] = ReshapeNeuron(0.0, 0.0)
    end
    ReshapeEnsemble(name, neurons, [], TrainTest, 1, input, shape)
end

function init(ensemble::ReshapeEnsemble, net::Net)
    resized = (ensemble.shape..., batch_size(net))
    for t in 1:net.time_steps
        input = get_buffer(net, symbol(ensemble.input.name, :value), t)
        ∇ = get_buffer(net, symbol(ensemble.input.name, :∇), t)
        set_buffer(net, symbol(ensemble.name, :value), reshape(input, resized), t)
        set_buffer(net, symbol(ensemble.name, :∇), reshape(∇, resized), t)
    end
end

function ReshapeLayer(name::Symbol, net::Net, input_ensemble, shape)
    ens = ReshapeEnsemble(name, input_ensemble, shape)
    add_ensemble(net, ens)
    ens
end


text:
@testset "for-basic" begin
    s = """
        @def v1 = [1, 2, 3]
        ~~~
        {{for v in v1}}
            v: {{fill v}}
        {{end}}
        ~~~
        """ |> fd2html_td
    @test isapproxstr(s, """
        v: 1
        v: 2
        v: 3
        """)

    s = """
        @def v1 = [[1,2], [3,4]]
        ~~~
        {{for (a,b) in v1}}
            a: {{fill a}} b: {{fill b}}
        {{end}}
        ~~~
        """ |> fd2html_td
    @test isapproxstr(s, """
        a: 1 b: 2
        a: 3 b: 4
        """)

    s = """
        @def v_1 = ("a"=>1, "b"=>2, "c"=>3)
        ~~~
        {{for (a,b) in v_1}}
            a: {{fill a}} b: {{fill b}}
        {{end}}
        ~~~
        """ |> fd2html_td
    @test isapproxstr(s, """
        a: a b: 1
        a: b b: 2
        a: c b: 3
        """)
end


text:
"""
    SplitByProportion(X::Array, Proportion::Float64 = 0.5)

Splits `X` Array along the observations dimension into a 2-Tuple based on the `Proportion`.
The form of the output is the following: ( X1, X2 )
"""
function SplitByProportion(X::Array, Proportion::Float64 = 0.5)
    FirstChunk = Int(floor(size(X)[1] * Proportion))
    return ( (X[ 1:FirstChunk, : ] ), (X[ ( FirstChunk + 1 ):end, : ] ) )
end

"""
    SplitByProportion(X::Array, Y::Array,Proportion::Float64 = 0.5)

Splits an `X` and Associated `Y` Array along the observations dimension into a 2-Tuple of 2-Tuples based on the `Proportion`.
The form of the output is the following: ( (X1, Y1), (X2, Y2) )
"""
function SplitByProportion(X::Array, Y::Array, Proportion::Float64 = 0.5)
    Y = forceMatrix(Y)
    FirstChunk = Int(floor(size(X)[1] * Proportion))
    return ((X[1:FirstChunk,:], Y[1:FirstChunk,:]), (X[(FirstChunk+1):end,:], Y[(FirstChunk+1):end,:])   )
end

"""
    VenetianBlinds(X)

Splits an `X` Array along the observations dimension into a 2-Tuple of 2-Tuples based on the whether it is even or odd.
The form of the output is the following: ( X1, X2 )
"""
VenetianBlinds(X) = ( X[1:2:end,:], X[2:2:end,:] )

"""
    VenetianBlinds(X,Y)

Splits an `X` and associated `Y` Array along the observation dimension into a 2-Tuple of 2-Tuples based on the whether it is even or odd.
The form of the output is the following: ( (X1,Y1), (X2, Y2) )
"""
VenetianBlinds(X,Y) = (( X[1:2:end,:], Y[1:2:end,:] ), ( X[2:2:end,:], Y[2:2:end,:] ) )

CICol(a::CartesianIndex) = collect(Tuple(a))[2]

CIRow(a::CartesianIndex) = collect(Tuple(a))[1]

"""
    KennardStone(X, TrainSamples; distance = "euclidean")

Returns the indices of the Kennard-Stone sampled exemplars (E), and those not sampled (O) as a 2-Tuple (E, O).

Note: there is an indexing bug here...

R. W. Kennard & L. A. Stone (1969) Computer Aided Design of Experiments, Technometrics, 111, 137-148, DOI: 10.1080/00401706.1969.10490666
"""
function KennardStone(X, TrainSamples; distance = "euclidean")
    Obs = 1:size(X)[1]
    FullSet = collect( Obs )
    DistMat = []
    if distance == "euclidean"
        DistMat = SquareEuclideanDistance( X )
    elseif distance == "manhattan"
        DistMat = ManhattanDistance( X )
    end
    @assert DistMat != []
    DistMat = hcat([[0] ; Obs], vcat( Obs', DistMat ))
    Candidates = Array{Int64, 1}(zeros( TrainSamples ))
    #Grab the 2 most seperated points
    Candidates[ 1 : 2 ] =  collect(Tuple(argmax( DistMat[ 2 : end, 2 : end ] ))) .|> Int64
    Candidates[ 1 : 2 ] .+= 1
    Available = setdiff( FullSet, Candidates[1:2])
    for sample in 3 : TrainSamples
        CurrentView = DistMat[ [ [ 1 ];   Available ],  [ [ 1 ];  Candidates[ 1:(sample - 1)  ]  ] ]
        #Find nearest observation(column) to each of the points that were selected(row)
        NearestNeighbors = argmin( CurrentView[2:end, 2:end], dims = 2)
        Values = CurrentView[2:end, 2:end][NearestNeighbors]#Get those values...
        FarNeighbor = argmax( Values )#Find the maximum
        #Get the column these values were associated with
        PointRelativeIndex = CurrentView[ CIRow(NearestNeighbors[ FarNeighbor ]) + 1, 1 ] |> Int64
        Candidates[ sample ] = PointRelativeIndex + 1
        Available = setdiff( FullSet, Candidates[ 1 : sample ]  ) #.- 1 )
    end
    return ( Candidates, Available )
end


text:
hascas(::Computed) = true  # maybe

function tryreact!(actr::Reactor{<:Computed}, a, rx::Reaction, offer::Union{Offer,Nothing})
    (; f) = actr.reagent
    r = f(a)::Union{Nothing,Failure,Reagent}
    if r === nothing
        return Block()
    elseif r isa Reagent
        return tryreact!(then(r, actr.continuation), a, rx, offer)
    else
        return r
    end
end

hascas(::Return) = false

function tryreact!(actr::Reactor{<:Return}, _, rx::Reaction, offer::Union{Offer,Nothing})
    (; value) = actr.reagent
    value isa Failure && return value
    return tryreact!(actr.continuation, value, rx, offer)
end

hascas(::Map) = false

function tryreact!(actr::Reactor{<:Map}, a, rx::Reaction, offer::Union{Offer,Nothing})
    (; f) = actr.reagent
    b = f(a)
    b isa Failure && return b
    return tryreact!(actr.continuation, b, rx, offer)
end

Reagents.Map(::Type{T}) where {T} = Reagents.Map{Type{T}}(T)


text:
include("main.jl")

REWRITE_DATA = true

BUSES, LINES, GENERATORS, rPTDF = load_case_data()
ERROR_VARIANCES_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

# Test dependency on variance estimation confidence
casename = "true_out_of_sample"
case_path = "results/$(casename)"

if REWRITE_DATA && isdir(case_path)
    rm(case_path, recursive=true)
end

if !(isdir(case_path)) || REWRITE_DATA
    # Only run if there are no saved results
    mkdir(case_path)

    # testvalues = [0.15, 0.12, 0.1, 0.08, 0.05, 0.03]
    testvalues = [0.1, 0.08, 0.05, 0.03, 0.01, 0.005]
    # var_conf = [0.1, 0.05, 0.01, 0.005, 0.001, 0.0005, 0.0001]
    var_conf = [0.1, 0.01, 0.005, 0.0001]
    # var_conf = [0.01]

    sample_size = 1500
    # Create global set of true observations to make tests comparable
    OBSERVATIONS = create_data(ERROR_VARIANCES_WC.*0.7, sample_size)

    cctest_dict = DataFrame()
    ddtest_dict = DataFrame()
    ldtest_dict = DataFrame()

    solvetimes_dr = []

    objective, bus_results, line_results = Solve_LinDist(BUSES, LINES, GENERATORS)
    det_alpha = vcat([1], zeros(length(BUSES)-1))
    # In Deterministic Case all imbalance is controlled by the substation
    bus_results = hcat(bus_results, DataFrame(α=det_alpha))
    ldtest = test_results(bus_results, line_results, objective; samples = sample_size, message="Standard CC", print_on=false)
    ldtest_dict = vcat(ldtest_dict, ldtest)

    for v in testvalues
        # Set Confidence Intervals
        η_v = v # 1 - Confidence for voltage limit cc
        η_g = v # 1 - Confidence for generation limit cc
        # SND Quantiles for global use
        snd = Normal(0,1)
        z_v = quantile(snd, 1-η_v)
        z_g = quantile(snd, 1-η_g)
        z_g2 = quantile(snd, 1-(η_g/2))

        buses_ad = change_load_same_pf(3.5)
        objective, bus_results, line_results, solvetimecc = Solve_DD_DR_CCLinDist(buses_ad, LINES, GENERATORS, ERROR_VARIANCES_SAMPLE)
        cctest = test_results(bus_results, line_results, objective; samples = sample_size, message="Standard CC", print_on=false)

        for vc in var_conf
            info("Runnig dr for η = $v and ξ=$vc")

            ξ = vc
            ERROR_VARIANCES_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity(create_new_data=false)

            ddobjective, ddbus_results, ddline_results, solvetimedr = Solve_DD_DR_CCLinDist(buses_ad, LINES, GENERATORS, ERROR_VARIANCES_WC)
            ddtest = test_results(ddbus_results, ddline_results, ddobjective; samples = sample_size, message="DD CC", print_on=false)

            param = DataFrame(xi=vc, eta=v)
            cctest_dict = vcat(cctest_dict, hcat(param, cctest)) # Redundant Data Saving for easier plotting
            ddtest_dict = vcat(ddtest_dict, hcat(param, ddtest))

        end
    end

    writetable("$(case_path)/cctest_dict.csv", cctest_dict)
    writetable("$(case_path)/ddtest_dict.csv", ddtest_dict)
    writetable("$(case_path)/ldtest_dict.csv", ldtest_dict)

else
    # If there is saved data load it
    info("Data already exists at $case_path")
    cctest_dict = readtable("$(case_path)/cctest_dict.csv")
    ddtest_dict = readtable("$(case_path)/ddtest_dict.csv")
end


text:
include("../src/BlueMesh7.jl")
using .BlueMesh7

function run()
    positions = generate_positions(dims=(30, 30), n=64)

    node_roles = ones(Int, length(positions))
    mesh = initialize_mesh(positions, node_roles)

    println("Running a simulation for 7 minutes (420000 steps, 1 step = 1ms)...")
    stats = start(mesh, minutes = 7)

    println("Packet delivery rate: $(stats.pdr), Worst delivery rate: $(stats.worstpdr), Average delay: $(stats.delay)ms, Centrality: $(stats.centrality), Average packet loss: $(stats.packetloss)")
end

run()


text:
# Adapted from Knet's src/data.jl (author: Deniz Yuret)
module Data

export
    DataLoader,
    VariableBatchSizeUniformSampleWithReplacement

using
    Random

using Base: @propagate_inbounds
using Random: GLOBAL_RNG
struct DataLoader{S,D,R<:AbstractRNG}
    data::D
    batchsize::Int
    nobs::Int
    partial::Bool
    imax::Int
    indices::Vector{Int}
    shuffle::Bool
    rng::R
    datasampler::S
end

function UniformSampleWithoutReplacement(d::DataLoader, i::Integer)
    if d.shuffle && i == 0
        shuffle!(d.rng, d.indices)
        @info "Using UniformSampleWithoutReplacement"
    end
    nexti = min(i + d.batchsize, d.nobs)
    ids = d.indices[i + 1:nexti]
    return (ids, nexti)
end

function UniformSampleWithReplacement(d::DataLoader, i::Integer)
    if d.shuffle && i == 0
        # only thing different from above... should probably combine
        rand!(d.rng, d.indices, [1:d.nobs;])
        @info "Using UniformSampleWithReplacement"
    end
    nexti = min(i + d.batchsize, d.nobs)
    ids = d.indices[i + 1:nexti]
    return (ids, nexti)
end

function VariableBatchSizeUniformSampleWithReplacement(d::DataLoader, i::Integer)
    if i == 0
        @info "Using VariableBatchSizeUniformSampleWithReplacement"
    end

    total_batches = d.nobs / d.batchsize
    sample_rate = float(d.batchsize) / float(d.nobs)

    ids = 0 # hack for scope, probably not Julian
    while true
        mask = rand(d.rng, Float32, size(d.indices)) .< sample_rate
        ids = d.indices[mask]
        length(ids) > 0 && break # ensure we have at least one observation per batch
    end
    # @show length(ids)
    nexti = i + 1

    if nexti >= total_batches
        nexti = d.imax # hacky
    end

    return (ids, nexti)
end

"""
    DataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG)

An object that iterates over mini-batches of `data`, each mini-batch containing
`batchsize` observations (except possibly the last one).

Takes as input a single data tensor, or a tuple (or a named tuple) of tensors.
The last dimension in each tensor is considered to be the observation dimension.

If `shuffle=true`, shuffles the observations each time iterations are re-started.
If `partial=false`, drops the last mini-batch if it is smaller than the batchsize.

The original data is preserved in the `data` field of the DataLoader.

Usage example:

    Xtrain = rand(10, 100)
    train_loader = DataLoader(Xtrain, batchsize=2)
    # iterate over 50 mini-batches of size 2
    for x in train_loader
        @assert size(x) == (10, 2)
        ...
    end

    train_loader.data   # original dataset

    # similar, but yielding tuples
    train_loader = DataLoader((Xtrain,), batchsize=2)
    for (x,) in train_loader
        @assert size(x) == (10, 2)
        ...
    end

    Xtrain = rand(10, 100)
    Ytrain = rand(100)
    train_loader = DataLoader((Xtrain, Ytrain), batchsize=2, shuffle=true)
    for epoch in 1:100
        for (x, y) in train_loader
            @assert size(x) == (10, 2)
            @assert size(y) == (2,)
            ...
        end
    end

    # train for 10 epochs
    using IterTools: ncycle
    Flux.train!(loss, ps, ncycle(train_loader, 10), opt)

    # can use NamedTuple to name tensors
    train_loader = DataLoader((images=Xtrain, labels=Ytrain), batchsize=2, shuffle=true)
    for datum in train_loader
        @assert size(datum.images) == (10, 2)
        @assert size(datum.labels) == (2,)
    end
"""
function DataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG,
                    datasampler=UniformSampleWithoutReplacement)
    batchsize > 0 || throw(ArgumentError("Need positive batchsize"))

    n = _nobs(data)
    if n < batchsize
        @warn "Number of observations less than batchsize, decreasing the batchsize to $n"
        batchsize = n
    end
    imax = partial ? n : n - batchsize + 1
    DataLoader(data, batchsize, n, partial, imax, [1:n;], shuffle, rng, datasampler)
end

# Return data in d.indices[i+1:i+batchsize]
@propagate_inbounds function Base.iterate(d::DataLoader, i=0)
    i >= d.imax && return nothing
    # if d.shuffle && i == 0
    #     shuffle!(d.rng, d.indices)
    # end
    # nexti = min(i + d.batchsize, d.nobs)
    # ids = d.indices[i + 1:nexti]
    ids, nexti = d.datasampler(d, i)
    batch = _getobs(d.data, ids)
    return (batch, nexti)
end

function Base.length(d::DataLoader)
    n = d.nobs / d.batchsize
    d.partial ? ceil(Int, n) : floor(Int, n)
end

_nobs(data::AbstractArray) = size(data)[end]

function _nobs(data::Union{Tuple,NamedTuple})
    length(data) > 0 || throw(ArgumentError("Need at least one data input"))
    n = _nobs(data[1])
    if !all(x -> _nobs(x) == n, Base.tail(data))
        throw(DimensionMismatch("All data should contain same number of observations"))
    end
    return n
end

_getobs(data::AbstractArray, i) = data[ntuple(i -> Colon(), Val(ndims(data) - 1))..., i]
_getobs(data::Union{Tuple,NamedTuple}, i) = map(Base.Fix2(_getobs, i), data)

Base.eltype(::DataLoader{D}) where D = D

end


text:
using Images
using ImageView
using HDF5, JLD
using LIBSVM
using DecisionTree
using MLBase
using Boltzmann
using NaiveBayes

include("constants.jl")
include("interp.jl")
include("shape.jl")
include("data.jl")
include("analysis.jl")
include("fit.jl")

function test_raw()
    img_dat, lab_dat, nzs = read_aligned(labeled_only=true)

    X = img_dat
    y = lab_dat

    p, n = size(X)
    train_frac = 0.9
    k = int(floor(train_frac * n))
    idxs = randperm(n)
    train = idxs[1:k]
    test = idxs[k+1:end]

    m = svmtrain(y[train], X[:, train], kernel_type=LIBSVM.Linear)

    accuracy = countnz(svmpredict(m, X[:,test])[1] .== y[test]) / countnz(test)
    return accuracy
end

function test_compressed()
    img_dat, lab_dat, nzs = read_aligned(labeled_only=true)
    rbm, nzs2 = load_model()

    X = transform(rbm, img_dat)
    y = lab_dat

    p, n = size(X)
    train_frac = 0.9
    k = int(floor(train_frac * n))
    idxs = randperm(n)
    train = idxs[1:k]
    test = idxs[k+1:end]

    m = svmtrain(y[train], X[:, train], kernel_type=LIBSVM.Linear)

    accuracy = countnz(svmpredict(m, X[:,test])[1] .== y[test]) / countnz(test)
    return accuracy
end

function test_tree_raw()
    img_dat, lab_dat, nzs = read_aligned(labeled_only=true)
    rbm, nzs2 = load_model()

    X = img_dat
    y = lab_dat

    accuracy = nfoldCV_tree(y, X', 0.9, 10)
    return mean(accuracy)
end

function test_tree_compressed()
    img_dat, lab_dat, nzs = read_aligned(labeled_only=true)
    rbm, nzs2 = load_model()

    X = transform(rbm, img_dat)
    y = lab_dat

    accuracy = nfoldCV_tree(y, X', 0.9, 10)
    return mean(accuracy)
end

function dothis()
    X = rand(1000, 100) * 100
    y = rand(1:5, 100)
    for i=1:100
        println("$(i)th iteration")
        m = svmtrain(y, X)
        yhat = svmpredict(m, X)
    end
    println("Done.")
end

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
"""
Abstract parent type of `MeshCat.Visualizer`.
"""
abstract type AbstractVisualizer end

# Interface methods
function setobject! end
function settransform! end
function delete! end
function setprop! end
function setcontrol! end

# Convenient shortcuts
"""
$(TYPEDSIGNATURES)

This will construct an appropriate three.js object from the given geometry and
a default material.
"""
setobject!(vis::AbstractVisualizer, geom::GeometryLike) = setobject!(vis, Object(geom))

"""
$(TYPEDSIGNATURES)

This will construct an appropriate three.js object from the given geometry and
the given material.
"""
setobject!(vis::AbstractVisualizer, geom::GeometryLike, material::AbstractMaterial) = setobject!(vis, Object(geom, material))

"""
Toggle visibility of the visualizer at the current path.

$(TYPEDSIGNATURES)
"""
setvisible!(vis::AbstractVisualizer, visible::Bool) = setprop!(vis, "visible", visible)


text:
using LinearAlgebra
@testset "Polynomial evaluation" begin
    A = [3 4; 5 6.6]

    coeff = [3.0; -1.0; 2.0; 0.1; 1.2]

    P1 =
        coeff[1] * I +
        coeff[2] * A +
        coeff[3] * A^2 +
        coeff[4] * A^3 +
        coeff[5] * A^4
    P2 = coeff[1] * I + coeff[2] * A
    P3 = coeff[1] * I

    poly_gens = Dict(
        "Monomial" => :graph_monomial,
        "Horner"   => :graph_horner,
        "PS"       => :graph_ps,
    )

    for key in keys(poly_gens)
        @testset "$key" begin
            (graph, cref) = eval(poly_gens[key])(coeff)
            @test eval_graph(graph, A) ≈ P1

            (graph, cref) = eval(poly_gens[key])(coeff[1:2])
            @test eval_graph(graph, A) ≈ P2

            (graph, cref) = eval(poly_gens[key])([coeff[1]])
            @test eval_graph(graph, A) ≈ P3

            if key == "Horner" # Add test here if scaling is available
                s = -0.89
                (graph_1, cref_1) =
                    eval(poly_gens[key])(coeff .* s .^ (0:length(coeff)-1))
                (graph_2, cref_2) = eval(poly_gens[key])(coeff, scaling = s)
                @test eval_graph(graph_1, A) ≈ eval_graph(graph_2, A)
            end
        end
    end
end


text:
function Base.show(io::IO, m::Clp.ClpCInterface.ClpModel)

    @printf "%0s \n" "Clp Model"
    @printf "   %-15s : %10s \n" "sense " Clp.ClpCInterface.optimization_direction(m) == -1 ? "maximize" : "minimize"
    @printf "   %-15s = %10s \n" "variables" Clp.ClpCInterface.get_num_cols(m)
    @printf "   %-15s = %10s \n" "constraints" Clp.ClpCInterface.get_num_rows(m)
    @printf "   %-15s = %10s \n" "coefficients" Clp.ClpCInterface.get_num_elements(m)
end

function settings_clp()
    lp = Clp.ClpCInterface.ClpModel()

    Clp.ClpCInterface.set_log_level(lp,0)
    return lp
end

function change_objective_coef(lp::ClpModel, objective::Number, value::Number = 1.0)
    obj = zeros(Float64, number_cols(lp))
    obj[objective] = value
    chg_obj_coefficients(lp, obj)
end

function change_objective_coef{T <: Number}(lp::ClpModel, objective::Array{T}, value::Array = [])
    if value == []
        value = ones(Float64, length(objective))
    end

    obj = zeros(Float64, number_cols(lp))
    obj[objective] = value
    chg_obj_coefficients(lp, obj)
end

function get_variable_bounds(lp::ClpModel)
    return get_col_lower(lp), get_col_upper(lp)
end

function get_variable_bounds(lp::ClpModel, index::Number)
    lb,ub = get_variable_bounds(lp)
    return lb[index], ub[index]
end

function get_constraint_bounds(lp::ClpModel)
    return get_row_lower(lp), get_row_upper(lp)
end

function get_constraint_matrix(lp::ClpModel)
    ClpCInterface.get_constraint_matrix(lp)
end

function get_direction(lp::ClpModel)
   optimization_direction(lp) == -1 ? "maximize" : "minimize"
end

function get_objective_value(lp::ClpModel)
    objective_value(lp)
end

function get_reduced_costs(lp::ClpModel)
    get_reduced_cost(lp)
end

function get_solution(lp::ClpModel)
    get_col_solution(lp)
end

function get_solution_status(lp::ClpModel)
    code = get_solution_status_code(lp)
    if code == 0
        return "Optimal"
    else
        return "Not optimal"
    end
end

function get_solution_status_code(lp::ClpModel)
    status(lp)
end

function get_slack(lp::ClpModel)
    dual_column_solution(lp)
end

function get_objective(lp::ClpModel)
    get_obj_coefficients(lp)
end

function set_direction(lp::ClpModel, direction::AbstractString)
    set_optimization_direction(lp, lowercase(direction) == "max" ? -1.0 : 1.0)
end

function set_col_bounds(lp::ClpModel, index::Number, lb::Number, ub::Number)
    lower, upper = get_variable_bounds(lp)
    lower[index] = lb
    upper[index] = ub
    chg_column_lower(lp, lower)
    chg_column_upper(lp, upper)
end

function set_col_bounds(lp::ClpModel, index::Number, fixed_value::Number)
    set_col_bounds(lp, index, fixed_value, fixed_value)
end

function set_objective(lp::ClpModel, objective::Number, value::Number = 1.0)
    obj = Clp.ClpCInterface.objective(lp)
    obj[objective] = value
    chg_obj_coefficients(lp, obj)
end

function setup_clp(model::Model, objective::AbstractString = "max")
    lp = settings_clp()

    load_problem(lp,model.S, model.lb, model.ub, model.c, model.b,model.b)
    set_direction(lp, objective)

    return lp
end

#function setup_clp{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b::Array{T}, S::SparseMatrixCSC, objective::AbstractString = "max")
#end

#function setup_clp{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, beq::Array{T}, Seq::SparseMatrixCSC, bineq::Array{T}, Sineq::SparseMatrixCSC, objective::AbstractString = "max")
#end

function solve(lp::ClpModel)
    primal(lp,0)
end

function solver_status(lp::ClpModel)
    get_solution_status(lp)
end

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
using Serialization
using JSON
using HTTP
using ZipFile

const DataDir = joinpath(dirname(@__FILE__()), "..", "data")
if !isdir(DataDir)
    mkdir(DataDir)
end

function http_get_binary_file(url, destfile)
    r = HTTP.request("GET", url)
    open(destfile, "w") do fh
        print(fh, String(r.body))
    end
end

function unzip(filename, writefiles = [])
    r = ZipFile.Reader(filename)
    for f in r.files
        if length(writefiles) == 0 || in(f.name, writefiles)
            open(f.name, "w") do fh
                write(fh, read(f))
            end
        end
    end
end

# To map Swedish postal numbers to Long and Lat we can use the export dumps
# of Geonames.org available here:
const GeonamesPostalDumpBaseUrl = "https://download.geonames.org/export/zip/"
# However, the resolution of this data is poor and most postal numbers are just
# mapped to the city that is closest to them.

function get_latest_geonames_postal_country_dump(country = "SE", redownload = false)
    filename = country * ".zip"
    url = GeonamesPostalDumpBaseUrl * filename
    destdir = DataDir
    if !isdir(destdir)
        mkpath(destdir)
    end
    destfile = joinpath(destdir, filename)
    resultfile = joinpath(destdir, country * ".txt")
    if redownload || !isfile(resultfile)
        http_get_binary_file(url, destfile)
        if isfile(destfile)
            cd(destdir) do
                unzip(destfile)
            end
            if isfile(resultfile)
                cd(destdir) do
                    isfile(destfile) && rm(destfile)
                    isfile("readme.txt") && rm("readme.txt")
                end
                return resultfile
            else
                error("Couldn't unpack file $destfile")
            end
        else
            error("Couldn't download from $url")
        end
    end
    return resultfile
end

function parse_geonames_country_dump(rawdatafile)
    d = Dict{String,Tuple{Float64, Float64, String, String, String}}()
    for line in readlines(rawdatafile)
        parts = split(line, r"\s+")
        postalnumber = replace(strip(parts[2] * parts[3]), r"\s+" => "")
        # Index from the back for lat and long since the number of names
        # in between sometimes varies
        lat = parse(Float64, parts[end-2])
        long = parse(Float64, parts[end-1])
        d[postalnumber] = (lat, long, parts[4], parts[5], parts[6])
    end
    d
end

GeonamesSwedishPostnrDictCache = nothing # This is parsed on first access below

function swedish_postnr_dict(redownload = false)
    rawdatafile = get_latest_geonames_postal_country_dump("SE", redownload)
    jdfile = rawdatafile * ".juliadata"
    global GeonamesSwedishPostnrDictCache
    if !isfile(jdfile)
        GeonamesSwedishPostnrDictCache = parse_geonames_country_dump(rawdatafile)
        serialize(jdfile, GeonamesSwedishPostnrDictCache)
        return GeonamesSwedishPostnrDictCache
    else
        try
            return deserialize(jdfile)
        catch _err
            dict = parse_geonames_country_dump(rawdatafile)
            serialize(dict, jdfile)
            return dict
        end
    end
end

function geonames_swedish_postnr2position(postnr::String, redownload = false)
    postnr = replace(strip(postnr), r"\s+" => "")
    return swedish_postnr_dict(redownload)[postnr]
end

# Google maps api string:
# https://maps.googleapis.com/maps/api/geocode/json?address=417+56+G%C3%B6teborg&key=YOUR_API_KEY
mutable struct GoogleMapsApi
    key::String
    mindelay::Float64
    lasttime::Float64
    GoogleMapsApi(key::String, mindelay::Float64 = 0.2) =
        new(key, mindelay, 0.0)
end

function addapikey(g::GoogleMapsApi, prefix, postfix = "")
    prefix * "&key=$(g.key)" * postfix
end

function get(g::GoogleMapsApi, url::String)
    t = time()
    if t - g.lasttime < g.mindelay
        sleep(g.mindelay - (t - g.lasttime))
    end
    g.lasttime = time()
    HTTP.get(url)
end

function getjson(g::GoogleMapsApi, url::String)
    searchurl = addapikey(g, url)
    JSON.parse(String(get(g, searchurl).body))
end

DefaultGoogleMapsApi = nothing

function setgooglemapsapikey!(key::String, defaultmindelay = 0.5)
    HTTP.setuseragent!("Mozilla/5.0 (Windows NT 10.0; rv:68.0) Gecko/20100101 Firefox/68.0")
    global DefaultGoogleMapsApi
    DefaultGoogleMapsApi = GoogleMapsApi(key, defaultmindelay)
end

removews(s) = replace(s, r"\s+" => "")
replacews(s, c="+") = replace(strip(s), r"\s+" => c)

function format_postnummer_for_search(nr::String)
    nr = removews(nr)
    occursin(r"\d{5}", nr) ? (nr[1:3] * "+" * nr[4:5]) : nr
end

function geocode(g::GoogleMapsApi, searchterms...)
    sts = map(replacews, searchterms)
    s = join(sts, "+")
    url = "https://maps.googleapis.com/maps/api/geocode/json?address=$(s)"
    js = getjson(g, url)
    if !isnothing(js) && haskey(js, "results") && (length(js["results"]) >= 1)
        return first(js["results"])
    else
        return nothing
    end
end

function geocode(searchterms...)
    global DefaultGoogleMapsApi
    if isnothing(DefaultGoogleMapsApi)
        error("You must first set the Google Maps API key with setgooglemapsapikey!")
    end
    geocode(DefaultGoogleMapsApi, searchterms...)
end

const SwePostnrDictCache = joinpath(DataDir, "SwePostnrDict.juliadata")
const SwePostnrDict = CachedDict(SwePostnrDictCache) do
    Dict{Union{String,Int},Tuple{Float64, Float64, String}}()
end

function swedish_postnr2position(postnr::String, addr::String = "";
    redownload = false, usegeonames = false)
    shortform = replace(strip(postnr), r"\s+" => "")
    global SwePostnrDict
    if !redownload && haskey(SwePostnrDict, shortform)
        return SwePostnrDict[shortform]
    end
    searchform = format_postnummer_for_search(shortform)
    res = if usegeonames
        lat, long, c1, c2, c3 = geonames_swedish_postnr2position(shortform)
        geocode(searchform, c1, c2, c3)
    else
        geocode(searchform, addr)
    end
    isnothing(res) && return nothing
    loc = res["geometry"]["location"]
    lat = loc["lat"]
    lng = loc["lng"]
    entry = (lat, lng, res["formatted_address"])
    update!(SwePostnrDict, entry, postnr)
    update!(SwePostnrDict, entry, shortform)
    if length(shortform) == 5
        longform = shortform[1:3] * " " * shortform[4:5]
        update!(SwePostnrDict, entry, longform)
    end
    asint = parse(Int, shortform)
    update!(SwePostnrDict, entry, asint, true)
    return entry
end

const AddressDictCache = joinpath(DataDir, "AddressDict.juliadata")
const AddressDict = CachedDict(AddressDictCache) do
    Dict{Union{String,Int},Tuple{Float64, Float64, String}}()
end

address2position(a::String) = address2position([a])

normalizeaddress(elems) =
    join(map(s -> replace(lowercase(strip(s)), r"\s{2}" => " "), elems), ", ")

function address2position(addressParts::Vector{String}; redownload = false)
    key = normalizeaddress(addressParts)
    global AddressDict
    if !redownload && haskey(AddressDict, key)
        return AddressDict[key]
    end
    res = geocode(addressParts...)
    isnothing(res) && return nothing
    loc = res["geometry"]["location"]
    lat = loc["lat"]
    lng = loc["lng"]
    entry = (lat, lng, res["formatted_address"])
    update!(AddressDict, entry, key)
    return entry
end


text:
# List of possible test functions to use
using Combinatorics

# The Sobol G Function

function sobolG(x::Array{Float64,1}; a = [0.0, 0.0, 0.0, 0.4, 0.4, 5.0])
    if length(x) != length(a)
        error("Inputs are of different lenths!")
    end

    return prod((abs.(4*x .- 2) .+ a)./(1 .+ a));
end

function sobolG_analytic(; a = [0.0, 0.0, 0.0, 0.4, 0.4, 5.0])
    y = zeros(2,length(a))

    # Compute Vi, the first order indicies
    Vi = (1/3) .* (1 .+ a) .^ (-2)
    y[1,:] = Vi

    # Use the Vi to compute the total order index, but first determine all
    # possible combinations of variables:
    I = monomialDegrees(length(a),length(a))
    I[(LinearIndices(I))[findall(x->x!=0,I)]] .= 1
    I = unique(I,dims=1)          # All unique combinations
    J = I[(LinearIndices(I[:,1]))[findall(x->x==1,I[:,1])],:]    # First element is non-zero

    #println(J)
    # Find the total variance
    Vn = 0;
    for k = 2:size(I,1) # First element is all zero, so skip
        Vn += prod(Vi[(LinearIndices(I[k,:]))[findall(x->x==1,I[k,:])]]);
    end

    for j = 1:length(a)
        # move the ith element to the start
        tempVi = [Vi[j]; Vi]
        deleteat!(tempVi,j+1)
        #println(tempVi)
        for k = 1:size(J,1)
            y[2,j] += prod(tempVi[J[k,:].==1])
        end
    end

    # Devide all the computed quantities by the total variance and return
    y = y./Vn
end

function monomialDegrees(numVars::Integer, maxDegree::Integer)
## Construct the Graded polynomial index
#               I = monomialDegrees(3,2)
#
#               I =
#
#                    0     0     0
#                    1     0     0
#                    0     1     0
#                    0     0     1
#                    2     0     0
#                    1     1     0
#                    1     0     1
#                    0     2     0
#                    0     1     1
#                    0     0     2
#

# Check for the trivial case
    if numVars==1
        	degrees = collect(0:maxDegree);
    else
    	# Define needed variables
    	degrees = zeros(Int64,1,numVars);
    	k = numVars;

    	# For each degree (starting at 1) up to maxDegree...
    	for n = 1:maxDegree
    	    # Find the integers that divide up the RV permutations
    	    dividers = reverse(transpose(hcat(collect(combinations(collect(1:(n+k-1)), k-1))...)),dims=1);
    	    div_diff = zeros(Int64,size(dividers,1),size(dividers,2)-1);
    	    for i = 1:size(dividers,1)
    		          div_diff[i,:] = diff(vec(dividers[i,:]));
    	    end

    	    # Form the index
            degrees = vcat(degrees, hcat(dividers[:,1] .- 1, div_diff .- 1, (n+k) .-dividers[:,end] .- 1));
    	end
    	# Output result
    	degrees
    end
end

max_stars_repo_name: mstobb/VarianceBasedSensitivity


text:
#!/usr/bin/env -S JULIA_VERSION={{ julia_version }} julia
{% if contest_id and problem_id and problem_alphabet %}
# contest: {{ contest_id }}, problem: {{ problem_id}}, alphabet: {{ problem_alphabet }}
{% endif %}
{% if prediction_success %}
{% endif %}
{% if mod or yes_str or no_str %}
{% endif %}
{% if mod %}
const MOD = {{ mod }}
{% endif %}
{% if yes_str %}
const YES = "{{ yes_str }}"
macro Y() :( println(YES); return ) end
{% endif %}
{% if no_str %}
const NO = "{{ no_str }}"
macro N() :( println(NO); return ) end
{% endif %}
{% if prediction_success %}

function solve({{ formal_arguments }})

end
{% endif %}

function main()
    {% if prediction_success %}
    tokens = Channel{String}(32)
    Task() do
        for line in eachline(@static VERSION < v"0.6" ? STDIN : stdin)
            startswith(line, '\0') && break
            for token in split(chomp(line))
                put!(tokens, token)
            end
        end
        close(tokens)
    end |> schedule
    {{ input_part }}
    solve({{ actual_arguments }})
    {% else %}
    # Failed to predict input format
    {% endif %}
end

if isempty(get(ENV, "ATCODER_LOCAL", ""))
    isempty(ARGS) && main()
else
    @eval begin
        {{ samples }}
        function test(sampleids...)
            isempty(sampleids) && return test(collect(1:length(samples))...)
            ostdin, ostdout = @static VERSION < v"0.6" ? (STDIN, STDOUT) : (stdin, stdout)
            rd, wr = first(redirect_stdout()), last(redirect_stdin())
            try
                map(sampleids) do sampleid
                    input, output = samples[sampleid]
                    print(ostdout, "Testing sample #$(sampleid)..."); flush(ostdout)
                    println(wr, input)
                    println(wr, "\0")
                    t = @elapsed main()
                    println()
                    result = strip(String(readavailable(rd)))
                    if result == output
                        println(ostdout, "OK ($t sec).")
                        true
                    else
                        println(ostdout, "ERROR ($t sec)")
                        println(ostdout, "== expected ==\n$output")
                        println(ostdout, "== actual ==\n$result\n")
                        false
                    end
                end |> all
            finally
                redirect_stdin(ostdin)
                redirect_stdout(ostdout)
            end
        end
        submit() = test() && (prog = @__FILE__; run(`atcoder-submit $prog`))
        !isinteractive() && submit()
    end
end


text:
using SymEngine
using Compat
using Compat.Test
import Compat: MathConstants.γ, MathConstants.e, MathConstants.φ, MathConstants.catalan

include("test-dense-matrix.jl")

x = symbols("x")
y = symbols(:y)
@vars z

# Check Basic conversions
@test eltype([Basic(u) for u in [1, 1/2, 1//2, pi, e]]) == Basic

# make sure @vars defines in a local scope
let
    @vars w
end
@test_throws UndefVarError isdefined(w)
@test_throws Exception show(Basic())

a = x^2 + x/2 - x*y*5
b = diff(a, x)
@test b == 2*x + 1//2 - 5*y

c = x + Rational(1, 5)
c = expand(c * 5)
@test c == 5*x + 1

c = x ^ 5
@test diff(c, x) == 5 * x ^ 4

c = x ^ y
@test c != y^x

c = Basic(-5)
@test abs(c) == 5
@test abs(c) != 4

show(a)
println()
show(b)
println()

## mathfuns
@test abs(Basic(-1)) == 1
@test sin(Basic(1)) == subs(sin(x), x, 1)
@test sin(PI) == 0
@test subs(sin(x), x, pi) == 0
@test sind(Basic(30)) == 1 // 2

## calculus
x,y = symbols("x y")
n = Basic(2)
ex = sin(x*y)
@test diff(log(x),x) == 1/x
@test diff(ex, x) == y * cos(x*y)
@test diff(ex, x, 2) == diff(diff(ex,x), x)
@test diff(ex, x, n) == diff(diff(ex,x), x)
@test diff(ex, x, y) == diff(diff(ex,x), y)
@test diff(ex, x, y,x) == diff(diff(diff(ex,x), y), x)
@test series(sin(x), x, 0, 2) == x
@test series(sin(x), x, 0, 3) == x - x^3/6

## ntheory
@test mod(Basic(10), Basic(4)) == 2
for j in [-3, 3], p in [-5,5]
    @test mod(Basic(j), Basic(p)) == mod(j, p)
end
@test mod(Basic(10), 4) == 2               # mod(::Basic, ::Number)
@test_throws MethodError mod(10, Basic(4)) # no mod(::Number, ::Basic)
@test gcd(Basic(10), Basic(4)) == 2
@test lcm(Basic(10), Basic(4)) == 20
@test binomial(Basic(5), 2) == 10

## type information
a = Basic(1)
b = Basic(1//2)
c = Basic(0.125)
@test isa(SymEngine.BasicType(a+a), SymEngine.BasicType{Val{:Integer}})
@test isa(SymEngine.BasicType(a+b), SymEngine.BasicType{Val{:Rational}})
@test isa(SymEngine.BasicType(a+c), SymEngine.BasicType{Val{:RealDouble}})
@test isa(SymEngine.BasicType(b+c), SymEngine.BasicType{Val{:RealDouble}})
@test isa(SymEngine.BasicType(c+c), SymEngine.BasicType{Val{:RealDouble}})

## can we do math with items of BasicType?
a1 = SymEngine.BasicType(a)
tot = a1
for i in 1:100
    global tot
    tot = tot + a1
end
@test tot == 101
sin(a1)

# samples of different types:
# (Int, Rational{Int}, Complex{Int}, Float64, Complex{Float64})
samples = (1, 1//2, (1 + 2im), 1.0, (1.0 + 0im))
## subs - check all different syntaxes and types
ex = x^2 + y^2
for val in samples
    @test subs(ex, x, val) == val^2 + y^2
    @test subs(ex, (x, val)) == val^2 + y^2
    @test subs(ex, x => val) == val^2 + y^2
    @test subs(ex, SymEngine.CMapBasicBasic(Dict(x=>val))) == val^2 + y^2
    @test subs(ex, Dict(x=>val)) == val^2 + y^2
end
# This probably results in a number of redundant tests (operator order).
for val1 in samples, val2 in samples
    @test subs(ex, (x, val1), (y, val2)) == val1^2 + val2^2
    @test subs(ex, x => val1, y => val2) == val1^2 + val2^2
end

## lambidfy
@test abs(lambdify(sin(Basic(1))) - sin(1)) <= 1e-14
fn = lambdify(exp(PI/2*x))
@test abs(fn(1) - exp(pi/2)) <= 1e-14
for val in samples
    ex2 = sin(x + val)
    fn2 = lambdify(ex2)
    @test abs(fn2(val) - sin(2*val)) <= 1e-14
end
@test lambdify(x^2)(3) == 9

A = [x 2; x 1]
@test lambdify(A, [x])(0) == [0 2; 0 1]
@test lambdify(A)(0) == [0 2; 0 1]
A = [x 2]
@test lambdify(A, [x])(1) == [1 2]
@test lambdify(A)(1) == [1 2]
@test isa(convert.(Expr, [0 x x+1]), Array{Expr})

## N
for val in samples
    @test N(Basic(val)) == val
end

for val in [π, γ, e, φ, catalan]
    @test N(Basic(val)) == val
end

@test !isnan(x)
@test isnan(Basic(0)/0)

## generic linear algebra
x = symbols("x")
A = [x 2; x 1]
#@test det(A) == -x
#@test det(inv(A)) == - 1/x
#(A \ [1,2])[1] == 3/x

## check that unique work (hash)
x,y,z = symbols("x y z")
@test length(SymEngine.free_symbols([x*y, y,z])) == 3

## check that callable symengine expressions can be used as functions for duck-typed functions
@vars x
function simple_newton(f, fp, x0)
    x = float(x0)
    while abs(f(x)) >= 1e-14
        x = x - f(x)/fp(x)
    end
    x
end
@test abs(simple_newton(sin(x), diff(sin(x), x), 3) - pi) <= 1e-14

## Check conversions SymEngine -> Julia
z,flt, rat, ima, cplx = btypes = [Basic(1), Basic(1.23), Basic(3//5), Basic(2im), Basic(1 + 2im)]

@test Int(z) == 1
@test BigInt(z) == 1
@test Float64(flt) ≈ 1.23
@test Real(flt) ≈ 1.23
@test convert(Rational{Int}, rat) == 3//5
@test convert(Complex{Int}, ima) == 2im
@test convert(Complex{Int}, cplx) == 1 + 2im

@test_throws InexactError convert(Int, flt)
@test_throws InexactError convert(Int, rat)

x = symbols("x")
Number[1 2 3 x]
@test_throws ArgumentError Int[1 2 3 x]

t = BigFloat(1.23)
@test !SymEngine.have_component("mpfr") || t == convert(BigFloat, convert(Basic, t))

@test typeof(N(Basic(-1))) != BigInt

# Check that libversion works. VersionNumber should always be >= 0.2.0
# since 0.2.0 is the first public release
@test SymEngine.libversion >= VersionNumber("0.2.0")

# Check that constructing Basic from Expr works
@vars x y
@test Basic(:(-2*x)) == -2*x
@test Basic(:(-3*x*y)) == -3*x*y
@test Basic(:((x-y)*-3)) == (x-y)*(-3)
@test Basic(:(-y)) == -y
@test Basic(:(-2*(x-2*y))) == -2*(x-2*y)

@test Basic(0)/0 == NAN
@test subs(1/x, x, 0) == Basic(1)/0

d = Dict(x=>y, y=>x)
@test subs(x + 2*y, d) == y + 2*x

@test sin(x+PI/4) != sin(x)
@test sin(PI/2-x) == cos(x)

f = SymFunction("f")
@test string(f(x, y)) == "f(x, y)"
@test string(f([x, y])) == "f(x, y)"
@test string(f(2*x)) == "f(2*x)"

@funs g, h
@test string(g(x, y)) == "g(x, y)"
@test string(h(x, y)) == "h(x, y)"

if SymEngine.libversion >= VersionNumber("0.4.0")
    # Function symbols
    @funs f, g, h, s
    @vars n, m
    expr = Basic("f(n) + f(n+1) + g(m)^2 - 3*h(n)*s(n)^3")
    @test function_symbols(expr) == [h(n), f(n), s(n), g(m), f(n+1)]
    @test get_name(f(n+1)) == "f"
    @test get_args(f(n+1)) == [n+1]
    @test get_args(f(n, m^2)) == [n, m^2]

    # Coefficients
    @vars x y
    expr = x^3 + 3*x^2*y + 3*x*y^2 + y^3 + 1
    @test coeff(expr, x, Basic(3)) == 1
    @test coeff(expr, x, Basic(2)) == 3*y
    @test coeff(expr, x, Basic(1)) == 3*y^2
    if SymEngine.libversion >= VersionNumber("0.5.0")
        @test coeff(expr, x, Basic(0)) == y^3 + 1
    end
end