text:
##########
#  TRef  #
##########

"""
    TRef(x)

Implementation of an abstract data structure that
automatically performs deepcopy during task copying.

Atomic (single-valued) TRef objects must be set or updated
by indexing. For example, to access `val = TRef(1)`, you
must use `val[]`.

Usage:

```julia
TRef(x)
```

Example:

```julia
# Initialize an atomic value
z = TRef(19.2)
z[] += 31

# Initialize a multi-index object
x = TRef([1 2 3; 4 5 6])
x[1, 3] = 999

# Initialize a TRef holding a dictionary.
d = TRef(Dict("A" => 1, 5 => "B"))
d["A"] = 10
```
"""
mutable struct TRef{T}
    val::T
end

Base.get(r::TRef) = r.val
Base.show(io::IO, r::TRef) = Base.show(io::IO, get(r))
Base.size(r::TRef) = Base.size(get(r))
Base.ndims(r::TRef) = Base.ndims(get(r))
tape_copy(r::TRef) = TRef(deepcopy(get(r)))

function Base.getindex(r::TRef, I::Vararg{Any,N}) where {N}
    return get(r)[I...]
end

function Base.setindex!(r::TRef, x, I::Vararg{Any,N}) where {N}
    d = get(r)
    if isa(d, Real)
        r.val = x
    else
        setindex!(d, x, I...)
    end
    return d
end

# Implements eltype, firstindex, lastindex, and iterate
# functions.
for F in (:eltype, :firstindex, :lastindex, :iterate)
    @eval Base.$F(a::TRef, args...) = $F(get(a), args...)
end

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
module Display

using UUIDs
import LibGit2

using ..Types

const colors = Dict(
    ' ' => :white,
    '+' => :light_green,
    '-' => :light_red,
    '↑' => :light_yellow,
    '~' => :light_yellow,
    '↓' => :light_magenta,
    '?' => :red,
)
const color_dark = :light_black

function git_file_stream(repo::LibGit2.GitRepo, spec::String; fakeit::Bool=false)::IO
    blob = try LibGit2.GitBlob(repo, spec)
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
        fakeit && return devnull
    end
    iob = IOBuffer(LibGit2.content(blob))
    close(blob)
    return iob
end

function status(ctx::Context, mode::PackageMode, use_as_api=false)
    env = ctx.env
    project₀ = project₁ = env.project
    manifest₀ = manifest₁ = env.manifest
    diff = nothing
    if !use_as_api
        pkg = ctx.env.pkg
        if pkg !== nothing
           printstyled("Project "; color=Base.info_color(), bold=true)
           println(pkg.name, " v", pkg.version)
        end
    end
    if env.git != nothing
        git_path = LibGit2.path(env.git)
        project_path = relpath(env.project_file, git_path)
        manifest_path = relpath(env.manifest_file, git_path)
        project₀ = read_project(git_file_stream(env.git, "HEAD:$project_path", fakeit=true))
        manifest₀ = read_manifest(git_file_stream(env.git, "HEAD:$manifest_path", fakeit=true))
    end
    if mode == PKGMODE_PROJECT || mode == PKGMODE_COMBINED
        # TODO: handle project deps missing from manifest
        m₀ = filter_manifest(in_project(project₀["deps"]), manifest₀)
        m₁ = filter_manifest(in_project(project₁["deps"]), manifest₁)
        diff = manifest_diff(ctx, m₀, m₁)
        if !use_as_api
            printpkgstyle(ctx, :Status, pathrepr(ctx, env.project_file); ignore_indent=true)
            print_diff(ctx, diff)
        end
    end
    if mode == PKGMODE_MANIFEST
        diff = manifest_diff(ctx, manifest₀, manifest₁)
        if !use_as_api
            printpkgstyle(ctx, :Status, pathrepr(ctx, env.manifest_file); ignore_indent=true)
            print_diff(ctx, diff)
        end
    elseif mode == PKGMODE_COMBINED
        p = not_in_project(merge(project₀["deps"], project₁["deps"]))
        m₀ = filter_manifest(p, manifest₀)
        m₁ = filter_manifest(p, manifest₁)
        c_diff = filter!(x->x.old != x.new, manifest_diff(ctx, m₀, m₁))
        if !isempty(c_diff)
            if !use_as_api
                printpkgstyle(ctx, :Status, pathrepr(ctx, env.manifest_file); ignore_indent=true)
                print_diff(ctx, c_diff)
            end
            diff = Base.vcat(c_diff, diff)
        end
    end
    return diff
end

function print_project_diff(ctx::Context, env₀::EnvCache, env₁::EnvCache)
    pm₀ = filter_manifest(in_project(env₀.project["deps"]), env₀.manifest)
    pm₁ = filter_manifest(in_project(env₁.project["deps"]), env₁.manifest)
    diff = filter!(x->x.old != x.new, manifest_diff(ctx, pm₀, pm₁))
    if isempty(diff)
        printstyled(color = color_dark, " [no changes]\n")
    else
        print_diff(ctx, diff)
    end
end

function print_manifest_diff(ctx::Context, env₀::EnvCache, env₁::EnvCache)
    diff = manifest_diff(ctx, env₀.manifest, env₁.manifest)
    diff = filter!(x->x.old != x.new, diff)
    if isempty(diff)
        printstyled(color = color_dark, " [no changes]\n")
    else
        print_diff(ctx, diff)
    end
end

struct VerInfo
    hash::Union{SHA1,Nothing}
    path::Union{String,Nothing}
    ver::Union{VersionNumber,Nothing}
    pinned::Bool
    repo::Union{Types.GitRepo, Nothing}
end

revstring(str::String) = occursin(r"\b([a-f0-9]{40})\b", str) ? str[1:7] : str

vstring(ctx::Context, a::VerInfo) =
    string((a.ver == nothing && a.hash != nothing) ? "[$(string(a.hash)[1:16])]" : "",
           a.ver != nothing ? "v$(a.ver)" : "",
           a.path != nothing ? " [$(pathrepr(ctx, a.path))]" : "",
           a.repo != nothing ? " #$(revstring(a.repo.rev))" : "",
           a.pinned == true ? " ⚲" : "",
           )

Base.:(==)(a::VerInfo, b::VerInfo) =
    a.hash == b.hash && a.ver == b.ver && a.pinned == b.pinned

≈(a::VerInfo, b::VerInfo) = a.hash == b.hash &&
    (a.ver == nothing || b.ver == nothing || a.ver == b.ver) &&
    (a.pinned == b.pinned)

struct DiffEntry
    uuid::UUID
    name::String
    old::Union{VerInfo,Nothing}
    new::Union{VerInfo,Nothing}
end

function print_diff(io::IO, ctx::Context, diff::Vector{DiffEntry})
    same = all(x.old == x.new for x in diff)
    for x in diff
        warnings = String[]
        if x.old != nothing && x.new != nothing
            if x.old ≈ x.new
                verb = ' '
                vstr = vstring(ctx, x.new)
            else
                if x.old.hash != x.new.hash && x.old.ver != x.new.ver
                    verb = x.old.ver == nothing || x.new.ver == nothing ||
                           x.old.ver == x.new.ver ? '~' :
                           x.old.ver < x.new.ver  ? '↑' : '↓'
                elseif x.old.ver == x.new.ver && x.old.pinned != x.new.pinned ||
                        x.old.repo != nothing || x.new.repo != nothing
                    verb = '~'
                else
                    verb = '?'
                    msg = x.old.hash == x.new.hash ?
                        "hashes match but versions don't: $(x.old.ver) ≠ $(x.new.ver)" :
                        "versions match but hashes don't: $(x.old.hash) ≠ $(x.new.hash)"
                    push!(warnings, msg)
                end
                vstr = (x.old.ver == x.new.ver && x.old.pinned == x.new.pinned) ?
                      vstring(ctx, x.new) :
                      vstring(ctx, x.old) * " ⇒ " * vstring(ctx, x.new)
            end
        elseif x.new != nothing
            verb = '+'
            vstr = vstring(ctx, x.new)
        elseif x.old != nothing
            verb = '-'
            vstr = vstring(ctx, x.old)
        else
            verb = '?'
            vstr = "[unknown]"
        end
        v = same ? "" : " $verb"
        printstyled(io, " [$(string(x.uuid)[1:8])]"; color = color_dark)
        printstyled(io, "$v $(x.name) $vstr\n"; color = colors[verb])
    end
end
# TODO: Use the Context stream
print_diff(ctx::Context, diff::Vector{DiffEntry}) = print_diff(stdout, ctx, diff)

function manifest_by_uuid(manifest::Dict)
    entries = Dict{UUID,Dict}()
    for (name, infos) in manifest, info in infos
        uuid = UUID(info["uuid"])
        haskey(entries, uuid) && @warn("Duplicate UUID in manifest: $uuid")
        entries[uuid] = merge(info, Dict("name" => name))
    end
    return entries
end

function name_ver_info(info::Dict)
    name = info["name"]
    hash = haskey(info, "git-tree-sha1") ? SHA1(info["git-tree-sha1"])    : nothing
    ver  = haskey(info, "version")       ? VersionNumber(info["version"]) : nothing
    path =  get(info, "path", nothing)
    pin  =  get(info, "pinned", false)
    if haskey(info, "repo-url")
        repo = Types.GitRepo(info["repo-url"], info["repo-rev"])
    else
        repo = nothing
    end
    name, VerInfo(hash, path, ver, pin, repo)
end

function manifest_diff(ctx::Context, manifest₀::Dict, manifest₁::Dict)
    diff = DiffEntry[]
    entries₀ = manifest_by_uuid(manifest₀)
    entries₁ = manifest_by_uuid(manifest₁)
    for uuid in union(keys(entries₀), keys(entries₁))
        name₀ = name₁ = v₀ = v₁ = nothing
        haskey(entries₀, uuid) && ((name₀, v₀) = name_ver_info(entries₀[uuid]))
        haskey(entries₁, uuid) && ((name₁, v₁) = name_ver_info(entries₁[uuid]))
        name₀ == nothing && (name₀ = name₁)
        name₁ == nothing && (name₁ = name₀)
        if name₀ == name₁
            push!(diff, DiffEntry(uuid, name₀, v₀, v₁))
        else
            push!(diff, DiffEntry(uuid, name₀, v₀, nothing))
            push!(diff, DiffEntry(uuid, name₁, nothing, v₁))
        end
    end
    sort!(diff, by=x->(x.uuid in keys(ctx.stdlibs), x.name, x.uuid))
end

function filter_manifest!(predicate, manifest::Dict)
    empty = String[]
    for (name, infos) in manifest
        filter!(infos) do info
            predicate(name, info)
        end
        isempty(infos) && push!(empty, name)
    end
    for name in empty
        pop!(manifest, name)
    end
    return manifest
end
filter_manifest(predicate, manifest::Dict) =
    filter_manifest!(predicate, deepcopy(manifest))

# This is precompilable, an anonymous function is not.
struct InProject{D <: Dict}
    deps::D
    neg::Bool
end
function (ip::InProject)(name::String, info::Dict)
    v = haskey(ip.deps, name) && haskey(info, "uuid") && ip.deps[name] == info["uuid"]
    return ip.neg ? !v : v
end
in_project(deps::Dict) = InProject(deps, false)
not_in_project(deps::Dict) = InProject(deps, true)

end # module

ext: jl

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
module lorenz_model

using DifferentialEquations, Flux, Optim, DiffEqFlux, Plots

model_params=[10.0,28.0,8.0/3]
all_ic=[[1.0,0.0,0.0],[0.5,1.0,0.0], [3,0.2,0.1]];
tend=10.0
tstart=0.0
sampling=0.05

function model(du,u,p,t)
    x,y,z = u
    σ,ρ,β = p
    du[1] = dx = σ*(y-x)
    du[2] = dy = x*(ρ-z) - y
    du[3] = dz = x*y - β*z
end

function InitPlot()

    mycolor=[:blue, :red, :green, :purple, :black, :cyan, :orange]
    t=tstart:sampling:100.0
    odedata=[ solve(ODEProblem(model,η,(tstart,100.0), model_params), Tsit5(), saveat=t, abstol=1e-8,reltol=1e-6) for η in all_ic];
    pl = plot(t,odedata[1][1,:],label="data 1", color=mycolor[1], lw=2, legend=false)

    for k in 2:length(all_ic)
        plot!(pl,t,odedata[k][1,:],label="data $k", color=mycolor[k], lw=2)
    end
    title!(pl, "Lorenz Chaos")
    xlabel!(pl, "Time")
    ylabel!(pl, "Variable")
    savefig("Lorenz_init.png")
    display(pl)

end

function predict_adjoint(param) # Our 1-layer neural network
    prob=ODEProblem(model,[1.0,0.0,0.0],(tstart,tend), model_params)
    Array(concrete_solve(prob,Tsit5(),param[1:3],param[4:end],saveat=tstart:sampling:tend,abstol=1e-8,reltol=1e-6))
end

# Generate some data to fit, and add some noise to it
data=predict_adjoint([1.0,0.0,0.0,10.0,28.0,8.0/3])
σN=0.05
data+=σN*randn(size(data))

#Init guess for the parameters
pguess=[0.8,0.1,0.1,10.5,28.5,9.0/3]

function loss_adjoint(param)
  prediction = predict_adjoint(param)
  loss = sum(abs2,prediction - data)
  loss
end

function train_model(;pguess=[1.1,0.05,-0.05,10.2,28.2,9.0/3])
    println("The initial loss function is $(loss_adjoint(pguess)[1])")
    resinit = DiffEqFlux.sciml_train(loss_adjoint,pguess,ADAM(), maxiters=4000)
    res = DiffEqFlux.sciml_train(loss_adjoint,resinit.minimizer,BFGS(initial_stepnorm=0.0001), maxiters=4000)
    println("The parameters are $(res.minimizer) with final loss value $(res.minimum)")
    return(res)

end

function plotFit(param; vars=(1,2,3), tend=100.0)

    validationPlot(param, param[1:3], vars=vars, tend=tend)

end

function validationPlot(param, ic; vars=(1,2,3), tend=100.0)

    sol_fit=solve(ODEProblem(model,ic,(0.0,tend),param[4:end]), Tsit5())
    sol_actual=solve(ODEProblem(model,ic,(0.0,tend),model_params), Tsit5())
    pl=plot(sol_fit, lw=2, legend=false, vars=vars, color=:green)
    plot!(pl,sol_actual,vars=vars,color=:blue)
    title!(pl,"Model Parameter Fits")
    savefig("valid_plot_lorenz2.png")
    display(pl)

end

end #module

lorenz_model.InitPlot()

resL=lorenz_model.train_model()

lorenz_model.plotFit(resL.minimizer; vars=(1,2,3), tend=30.0)

lorenz_model.plotFit(resL.minimizer; vars=(0,1), tend=30.0)

max_stars_repo_name: khannay/FittingParamsDiffEqFlux

max_issues_repo_name: khannay/FittingParamsDiffEqFlux

max_forks_repo_name: khannay/FittingParamsDiffEqFlux


text:
@testset "Unit tests" begin
    include("TableLogger.jl")
    include("utils.jl")
    include("header.jl")
    include("line.jl")
    include("format.jl")
end


text:
using MIDI
using StatsBase

# export scale_identification, SCALES
export note_to_fundamental, pitch_frequency

const SCALES = Dict{String, Vector{String}}()

SCALES["C Major/A minor"] = ["C", "D", "E", "F", "G", "A", "B"]
SCALES["C♯ Major/A♯ minor"] = ["C", "C♯", "D♯", "F", "F♯", "G♯", "A♯"]
SCALES["D Major/B minor"] = ["C♯", "D", "E", "F♯", "G", "A", "B"]
SCALES["D♯ Major/C minor"] = ["C", "D", "D♯", "F", "G", "G♯", "A♯"]
SCALES["E Major/C♯ minor"] = ["C♯", "D♯", "E", "F♯", "G♯", "A", "B"]
SCALES["F Major/D minor"] = ["C", "D", "E", "F", "G", "A", "A♯"]
SCALES["F♯ Major/D♯ minor"] = ["C♯", "D♯", "F", "F♯", "G♯", "A♯", "B"]
SCALES["G Major/E minor"] = ["C", "D", "E", "F♯", "G", "A", "B"]
SCALES["G♯ Major/F minor"] = ["C","C♯", "D♯", "F", "G", "G♯", "A♯"]
SCALES["A Major/F♯ minor"] = ["C♯", "D", "E", "F♯", "G♯", "A", "B"]
SCALES["A♯ Major/G minor"] = ["C", "D", "D♯", "F", "G", "A", "A♯"]
SCALES["B Major/G♯ minor"] = ["C♯", "D♯", "E", "F♯", "G♯", "A♯", "B"]
SCALES["C minor harmonic"] = ["C", "D", "D♯", "F", "G", "G♯", "B"]
SCALES["C♯ minor harmonic"] = ["C", "C♯","D♯", "E", "F♯", "G♯", "A"]
SCALES["D minor harmonic"] = ["C♯", "D", "E", "F", "G", "A", "A♯"]
SCALES["D♯ minor harmonic"] = ["B", "D", "D♯", "F", "F♯", "G♯", "A♯"]
SCALES["E minor harmonic"] = ["C", "D♯", "E", "F♯", "G", "A", "B"]
SCALES["F minor harmonic"] = ["C", "C♯", "E", "F", "G", "G♯", "A♯"]
SCALES["F♯ minor harmonic"] =  ["C♯", "D", "F", "F♯", "G♯", "A", "B"]
SCALES["G minor harmonic"] = ["C", "D", "D♯", "F♯", "G", "A", "A♯"]
SCALES["G♯ minor harmonic"] = ["C♯", "D♯", "E", "G", "G♯", "A♯", "B"]
SCALES["A minor harmonic"] = ["C", "D", "E", "F", "G♯", "A", "B"]
SCALES["A♯ minor harmonic"] = ["C", "C♯", "D♯", "F", "F♯", "A", "A♯"]
SCALES["B minor harmonic"] = ["C♯", "D", "E", "F♯", "G", "A♯", "B"]
SCALES["C minor melodic"] = ["C", "D", "D♯", "F", "G", "A", "B"]
SCALES["C♯ minor melodic"] = ["C", "C♯","D♯", "E", "F♯", "G♯", "A♯"]
SCALES["D minor melodic"] = ["C♯", "D", "E", "F", "G", "A", "B"]
SCALES["D♯ minor harmonic"] = ["C", "D", "D♯", "F", "F♯", "G♯", "A♯"]
SCALES["E minor melodic"] = ["C♯", "D♯", "E", "F♯", "G", "A", "B"]
SCALES["F minor melodic"] = ["C", "D", "E", "F", "G", "G♯", "A♯"]
SCALES["F♯ minor melodic"] =  ["C♯", "D♯", "F", "F♯", "G♯", "A", "B"]
SCALES["G minor melodic"] = ["C", "D", "E", "F♯", "G", "A", "A♯"]
SCALES["G♯ minor melodic"] = ["C♯", "D♯", "F", "G", "G♯", "A♯", "B"]
SCALES["A minor melodic"] = ["C", "D", "E", "F♯", "G♯", "A", "B"]
SCALES["A♯ minor melodic"] = ["C", "C♯", "D♯", "F", "G", "A", "A♯"]
SCALES["B minor melodic"] = ["C♯", "D", "E", "F♯", "G♯", "A♯", "B"]

"""
    note_to_fundamental(note(s))

Return a `String` or `Vector{String}` with the fundamental pitch of the notes
(i.e. without the octave information).
"""
note_to_fundamental(n::AbstractNote) = MIDI.PITCH_TO_NAME[mod(n.pitch, 12)]
note_to_fundamental(notes) = note_to_fundamental.(notes)

"""
    pitch_frequency(notes)

Return the pitches (without octaves) in the given `notes`,
sorted by most to least frequent.
"""
pitch_frequency(notes) = pitch_frequency(note_to_fundamental(notes))
function pitch_frequency(fundamentals::Vector{String})
    occurency = sort(Dict(value => key for (key, value) in countmap(fundamentals)), rev = true)
    most_frequent = collect(values(occurency))
    return most_frequent
end

"""
    scale_identification(notes)

Return the most probable scale that represents the given `notes`.
"""
function scale_identification(notes)
    fund = note_to_fundamental(notes)
    mostfreq = pitch_frequency(fund)[1:7] # only 7 most frequent
    for (k,v) in SCALES
        tester = true
        for n in mostfreq
            if n ∉ v
                tester = false
            end
        end
        if tester == true
            return k
        end
    end
    if tester == false
        error("Atonal/modulating musical key center or Unregistered exotic scale")
    end
end


text:
@testset "Constraint List" begin
    model = Cartpole()
    n,m = size(model)
    N = 11
    x,u = rand(model)
    z = KnotPoint(x,u,0.1)

    #--- Generate some constraints
    # Circle Constraint
    xc = SA[1,1,1]
    yc = SA[1,2,3]
    r  = SA[1,1,1]
    cir = CircleConstraint(n, xc, yc, r)

    # Goal Constraint
    xf = @SVector rand(n)
    goal = GoalConstraint(xf)

    # Linear Constraint
    p = 5
    A = @SMatrix rand(p,n+m)
    b = @SVector rand(p)
    lin = LinearConstraint(n,m,A,b, Inequality())

    # Bound Constraint
    xmin = -@SVector rand(n)
    xmax = +@SVector rand(n)
    umin = -@SVector rand(m)
    umax = +@SVector rand(m)
    bnd = BoundConstraint(n,m, x_min=xmin, x_max=xmax, u_min=umin, u_max=umax)

    # Dynamics Constraint
    dyn = TO.DynamicsConstraint(model, N)

    #--- Create a List
    cons = ConstraintList(n,m,N)
    add_constraint!(cons, cir, 1:N)
    @test cons.constraints[1] === cir
    @test cons[1] === cir
    @test cons.inds[1] == 1:N
    @test cons.p == fill(length(cir), N)

    add_constraint!(cons, goal, N)
    @test cons.constraints[2] === goal
    @test cons[2] === goal
    @test cons.inds[2] == N:N
    @test cons.p[1:N-1] == fill(length(cir), N-1)
    @test cons.p[end] == length(cir) + length(goal)

    add_constraint!(cons, lin, 1:4, 1)
    @test cons[1] === lin
    @test cons[2] === cir
    @test cons[end] === goal
    @test cons.inds[1] === 1:4
    @test cons.p[1:4] == fill(length(cir)+length(lin), 4)
    @test cons.p[5:N-1] == fill(length(cir), N-1-4)
    @test length(cons) == 3

    cons2 = copy(cons)
    add_constraint!(cons, bnd, 1:N-1)
    @test length(cons) == 4
    @test length(cons2) == 3

    add_constraint!(cons, dyn, 1:N-1, 2)
    @test cons[2] == dyn
    @test collect(zip(cons))[2] == (1:N-1, dyn)
    @test cons.p[1] == length(cir) + length(bnd) + length(lin) + n
    @test length(cons) == 5

    @test TO.has_dynamics_constraint(cons) == true
    @test TO.has_dynamics_constraint(cons2) == false
    @test TO.num_constraints(cons) === cons.p
    @test TO.num_constraints(cons2) !== cons.p
    @test cons[end] == bnd
    sort!(cons)
    @test cons[end] == dyn
    @test cons[end-1] == bnd
    @test cons[1] == lin
    @test length(cons) == 5

    # Try adding a constraint with incorrect dimensions
    lin2 = LinearConstraint(2, 1, rand(3,2), rand(3), Inequality(), 1:2)
    @test_throws AssertionError add_constraint!(cons, lin2, 1:4)
    @test_throws AssertionError add_constraint!(cons2, dyn, 1:N+1)

    # Test iteration
    conlist = [lin, cir, goal, bnd, dyn]
    @test all(conlist .=== [con for con in cons])
    @test length.(cons) == length.(conlist)
    @test eltype(cons) == TO.AbstractConstraint
end


text:
# https://github.com/JuliaLabs/brutus/blob/9d87d7aa543090ab6e8b98dc7413476c107a7075/Brutus/src/Brutus.jl#L98-L146

function get_methodinstance(@nospecialize(sig);
                            world=Base.get_world_counter(),
                            interp=Core.Compiler.NativeInterpreter(world))
    ms = Base._methods_by_ftype(sig, 1, Base.get_world_counter())
    @assert length(ms) == 1
    m = ms[1]
    mi = ccall(:jl_specializations_get_linfo,
               Ref{Core.MethodInstance}, (Any, Any, Any),
               m[3], m[1], m[2])
    return mi
end

function code_ircode_by_signature(@nospecialize(sig);
                                  world=Base.get_world_counter(),
                                  interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(ccall(:jl_specializations_get_linfo,
                              Ref{Core.MethodInstance},
                              (Any, Any, Any),
                              data[3], data[1], data[2]);
                        world=world, interp=interp)
            for data in Base._methods_by_ftype(sig, -1, world)]
end

function code_ircode(@nospecialize(f), @nospecialize(types=Tuple);
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(mi; world=world, interp=interp)
            for mi in Base.method_instances(f, types, world)]
end

function code_ircode(mi::Core.Compiler.MethodInstance;
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = @static if hasmethod(
        Core.Compiler.InferenceState,
        (Core.Compiler.InferenceResult, Bool, Core.Compiler.AbstractInterpreter),
    )
        Core.Compiler.InferenceState(result, #=cached=# false, interp)
    else
        Core.Compiler.InferenceState(result, #=cache=# :no, interp)
    end
    frame === nothing && return nothing
    nargs = Int((mi.def::Method).nargs)
    if Core.Compiler.typeinf(interp, frame)
        opt_params = Core.Compiler.OptimizationParams(interp)
        opt = Core.Compiler.OptimizationState(frame, opt_params, interp)
        ir = @static if hasmethod(
            Core.Compiler.run_passes,
            (Core.CodeInfo, Int, Core.Compiler.OptimizationState),
        )
            Core.Compiler.run_passes(opt.src, nargs - 1, opt)
        else
            Core.Compiler.run_passes(opt.src, opt, result)
        end
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return nothing
    # TODO(yhls): Fix this upstream
    resize!(ir.argtypes, nargs)
    return ir => Core.Compiler.widenconst(result.result)
end


text:
module Kannaadi

include("utf8.jl")

end # module


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

module IndexController

using Genie.Renderer
using Genie.Router

function index()
    html!(:index, :index)
end

end


text:
module APIKeyHandlers

using ConfParser
using Printf

export AbstractKeyHandler
export get_name, get_key, get_api_url
export set_name!, set_key!, set_api_url!
export load_key!, load!, save!

export GenericKeyHandler

include("abstractkeyhandler.jl")
include("generickeyhandler.jl")

end # module


text:
## Common code for CategoricalArray and NullableCategoricalArray

import Base: convert, copy, getindex, setindex!, similar, size, linearindexing, vcat

# Used for keyword argument default value
_ordered(x::AbstractCategoricalArray) = ordered(x)
_ordered(x::AbstractNullableCategoricalArray) = ordered(x)
_ordered(x::Any) = false

function reftype(sz::Int)
    if sz <= typemax(UInt8)
        return UInt8
    elseif sz <= typemax(UInt16)
        return UInt16
    elseif sz <= typemax(UInt32)
        return UInt32
    else
        return UInt64
    end
end

for (A, V, M) in ((:CategoricalArray, :CategoricalVector, :CategoricalMatrix),
                  (:NullableCategoricalArray, :NullableCategoricalVector, :NullableCategoricalMatrix))
    @eval begin
        # Uninitialized array constructors

        $A{T, N}(::Type{T}, dims::NTuple{N,Int}; ordered=false) =
            $A(zeros(DefaultRefType, dims), CategoricalPool{T}(ordered))
        function $A{T}(::Type{T}, dims::Int...; ordered=false)
            A = $A{T}(dims)
            ordered!(A, ordered)
            A
        end
        $A(dims::Int...; ordered=false) = $A{String}(dims, ordered=ordered)

        $A{T, N, R}(::Type{CategoricalValue{T, R}}, dims::NTuple{N,Int}) = $A{T, N, R}(dims)
        $A{T, N}(::Type{CategoricalValue{T}}, dims::NTuple{N,Int}) = $A{T, N}(dims)
#        $A{N}(::Type{CategoricalValue}, dims::NTuple{N,Int}) = $A{String, N}(dims)

        @compat (::Type{$A{T, N, R}}){T, N, R}(dims::NTuple{N,Int}; ordered=false) =
            $A{T, N, R}(zeros(R, dims), CategoricalPool{T, R}(ordered))
        @compat (::Type{$A{T, N}}){T, N}(dims::NTuple{N,Int}; ordered=false) =
            $A{T, N, DefaultRefType}(dims, ordered=ordered)
        @compat (::Type{$A{T}}){T, N}(dims::NTuple{N,Int}; ordered=false) =
            $A{T, N}(dims, ordered=ordered)
        @compat (::Type{$A{T, 1, R}}){T, R}(m::Int; ordered=false) =
            $A{T, 1, R}((m,), ordered=ordered)
        # R <: Integer is required to prevent default constructor from being called instead
        @compat (::Type{$A{T, 2, R}}){T, R <: Integer}(m::Int, n::Int; ordered=false) =
            $A{T, 2, R}((m, n), ordered=ordered)
        @compat (::Type{$A{T, 3, R}}){T, R}(m::Int, n::Int, o::Int; ordered=false) =
            $A{T, 3, R}((m, n, o), ordered=ordered)
        @compat (::Type{$A{T}}){T}(m::Int; ordered=false) =
            $A{T}((m,), ordered=ordered)
        @compat (::Type{$A{T}}){T}(m::Int, n::Int; ordered=false) =
            $A{T}((m, n), ordered=ordered)
        @compat (::Type{$A{T}}){T}(m::Int, n::Int, o::Int; ordered=false) =
            $A{T}((m, n, o), ordered=ordered)

        @compat (::Type{$A{CategoricalValue{T, R}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                    ordered=false) =
            $A{T, N, R}(dims, ordered=ordered)
        @compat (::Type{$A{CategoricalValue{T}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                 ordered=false) =
            $A{T, N, R}(dims, ordered=ordered)
        @compat (::Type{$A{CategoricalValue{T, R}, N}}){T, N, R}(dims::NTuple{N,Int};
                                                                 ordered=false) =
            $A{T, N, R}(dims, ordered=ordered)
        @compat (::Type{$A{CategoricalValue{T}, N}}){T, N}(dims::NTuple{N,Int};
                                                           ordered=false) =
            $A{T, N}(dims, ordered=ordered)
#        @compat (::Type{$A{CategoricalValue, N}}){N}(dims::NTuple{N,Int};
#                                                     ordered=false) =
#            $A{String, N}(dims, ordered=ordered)
#        @compat (::Type{$A{CategoricalValue}}){N}(dims::NTuple{N,Int};
#                                                  ordered=false) =
#            $A{String, N}(dims, ordered=ordered)

if VERSION >= v"0.5.0-dev"
        $V{T}(::Type{T}, m::Integer; ordered=false) = $A{T}((m,), ordered=ordered)
        $V(m::Integer; ordered=false) = $A(m, ordered=ordered)
end
        @compat (::Type{$V{T}}){T}(m::Int; ordered=false) = $A{T}((m,), ordered=ordered)

if VERSION >= v"0.5.0-dev"
        $M{T}(::Type{T}, m::Int, n::Int; ordered=false) = $A{T}((m, n), ordered=ordered)
        $M(m::Int, n::Int; ordered=false) = $A(m, n, ordered=ordered)
end
        @compat (::Type{$M{T}}){T}(m::Int, n::Int; ordered=false) = $A{T}((m, n), ordered=ordered)

        ## Constructors from arrays

        # This method is needed to ensure ordered!() only mutates a copy of A
        @compat function (::Type{$A{T, N, R}}){T, N, R}(A::$A{T, N, R};
                                                        ordered=_ordered(A))
            ret = copy(A)
            ordered!(ret, ordered)
            ret
        end

        # Note this method is also used for CategoricalArrays when T, N or R don't match
        @compat function (::Type{$A{T, N, R}}){T, N, R}(A::AbstractArray;
                                                        ordered=_ordered(A))
            ret = convert($A{T, N, R}, A)
            ordered!(ret, ordered)
            ret
        end

        @compat (::Type{$A{T, N, R}}){T<:CategoricalValue, N, R}(A::AbstractArray;
                                                                 ordered=_ordered(A)) =
            $A{T.parameters[1], N, R}(A, ordered=ordered)

        # From AbstractArray
        @compat (::Type{$A{T, N}}){S, T, N}(A::AbstractArray{S, N}; ordered=_ordered(A)) =
            $A{T, N, DefaultRefType}(A, ordered=ordered)
        @compat (::Type{$A{T}}){S, T, N}(A::AbstractArray{S, N}; ordered=_ordered(A)) =
            $A{T, N}(A, ordered=ordered)
        @compat (::Type{$A}){T, N}(A::AbstractArray{T, N}; ordered=_ordered(A)) =
            $A{T, N}(A, ordered=ordered)

        @compat (::Type{$V{T}}){S, T}(A::AbstractVector{S}; ordered=_ordered(A)) =
            $A{T, 1}(A, ordered=ordered)
        @compat (::Type{$V}){T}(A::AbstractVector{T}; ordered=_ordered(A)) =
            $A{T, 1}(A, ordered=ordered)

        @compat (::Type{$M{T}}){S, T}(A::AbstractMatrix{S}; ordered=_ordered(A)) =
            $A{T, 2}(A, ordered=ordered)
        @compat (::Type{$M}){T}(A::AbstractMatrix{T}; ordered=_ordered(A)) =
            $A{T, 2}(A, ordered=ordered)

        # From CategoricalArray (preserve R)
        @compat (::Type{$A{T, N}}){S, T, N, R}(A::$A{S, N, R}; ordered=_ordered(A)) =
            $A{T, N, R}(A, ordered=ordered)
        @compat (::Type{$A{T}}){S, T, N, R}(A::$A{S, N, R}; ordered=_ordered(A)) =
            $A{T, N, R}(A, ordered=ordered)
        @compat (::Type{$A}){T, N, R}(A::$A{T, N, R}; ordered=_ordered(A)) =
            $A{T, N, R}(A, ordered=ordered)

        @compat (::Type{$V{T}}){S, T, R}(A::$V{S, R}; ordered=_ordered(A)) =
            $A{T, 1, R}(A, ordered=ordered)
        @compat (::Type{$V}){T, R}(A::$V{T, R}; ordered=_ordered(A)) =
            $A{T, 1, R}(A, ordered=ordered)

        @compat (::Type{$M{T}}){S, T, R}(A::$M{S, R}; ordered=_ordered(A)) =
            $A{T, 2, R}(A, ordered=ordered)
        @compat (::Type{$M}){T, R}(A::$M{T, R}; ordered=_ordered(A)) =
            $A{T, 2, R}(A, ordered=ordered)

        ## Conversion methods

        # From AbstractArray
        convert{T, N}(::Type{$A{T, N}}, A::AbstractArray{T, N}) =
            convert($A{T, N, DefaultRefType}, A)
        convert{T, N}(::Type{$A{T}}, A::AbstractArray{T, N}) = convert($A{T, N}, A)
        convert{S, T, N}(::Type{$A{T}}, A::AbstractArray{S, N}) = convert($A{T, N}, A)
        convert{T, N}(::Type{$A}, A::AbstractArray{T, N}) = convert($A{T, N}, A)

        convert{T, N, R}(::Type{$A{CategoricalValue{T, R}, N}}, A::AbstractArray{T, N}) =
            convert($A{T, N, R}, A)
        convert{T, N}(::Type{$A{CategoricalValue{T}, N}}, A::AbstractArray{T, N}) =
            convert($A{T, N}, A)

        convert{T}(::Type{$V{T}}, A::AbstractVector) = convert($V{T, DefaultRefType}, A)
        convert{T}(::Type{$V}, A::AbstractVector{T}) = convert($V{T}, A)
        convert{T}(::Type{$V{T}}, A::$V{T}) = A
        convert(::Type{$V}, A::$V) = A

        convert{T}(::Type{$M{T}}, A::AbstractMatrix) = convert($M{T, DefaultRefType}, A)
        convert{T}(::Type{$M}, A::AbstractMatrix{T}) = convert($M{T}, A)
        convert{T}(::Type{$M{T}}, A::$M{T}) = A
        convert(::Type{$M}, A::$M) = A

        function convert{S, T, N, R}(::Type{$A{T, N, R}}, A::AbstractArray{S, N})
            res = $A{T, N, R}(size(A))
            copy!(res, A)

            if method_exists(isless, (T, T))
                levels!(res, sort(levels(res)))
            end

            res
        end

        # From CategoricalArray (preserve R)
        function convert{T, N, R}(::Type{$A{T, N, R}}, A::$A)
            refs = convert(Array{R, N}, A.refs)
            pool = convert(CategoricalPool{T, R}, A.pool)
            ret = $A(refs, pool)
            ordered!(ret, ordered(A))
            ret
        end
        convert{S, T, N, R}(::Type{$A{T, N}}, A::$A{S, N, R}) =
            convert($A{T, N, R}, A)
        convert{S, T, N, R}(::Type{$A{T}}, A::$A{S, N, R}) =
            convert($A{T, N, R}, A)
        convert{T, N, R}(::Type{$A}, A::$A{T, N, R}) =
            convert($A{T, N, R}, A)

        # R<:Integer is needed for this method to be considered more specific
        # than the generic one above (JuliaLang/julia#18443)
        convert{T, N, R<:Integer}(::Type{$A{T, N, R}}, A::$A{T, N, R}) = A
        convert{T, N}(::Type{$A{T, N}}, A::$A{T, N}) = A
        convert{T}(::Type{$A{T}}, A::$A{T}) = A
        convert(::Type{$A}, A::$A) = A

        ## Other methods

        similar{S, T, M, N, R}(A::$A{S, M, R}, ::Type{T}, dims::NTuple{N, Int}) =
            $A{T, N, R}(dims; ordered=ordered(A))

        function compact{T, N}(A::$A{T, N})
            R = reftype(length(index(A.pool)))
            convert($A{T, N, R}, A)
        end

        function vcat(A::$A...)
            newlevels, isordered = mergelevels(map(levels, A)...)

            refs = [indexin(index(a.pool), newlevels)[a.refs] for a in A]
            $A(DefaultRefType[refs...;],
               CategoricalPool(newlevels, isordered && all(ordered, A)))
        end
    end
end

function @compat(Base.:(==))(A::CatArray, B::CatArray)
    if size(A) != size(B)
        return false
    end
    if A.pool === B.pool
        for (a, b) in zip(A.refs, B.refs)
            if a != b
                return false
            end
        end
    else
        for (a, b) in zip(A, B)
            if a != b
                return false
            end
        end
    end
    return true
end

size(A::CatArray) = size(A.refs)
linearindexing{T <: CatArray}(::Type{T}) = Base.LinearFast()

setindex!(A::CatArray, v::Any, i::Int) = A.refs[i] = get!(A.pool, v)
setindex!{T}(A::CatArray, v::CategoricalValue{T}, i::Int) =
    A.refs[i] = get!(A.pool, convert(T, v))

# Method preserving levels and more efficient than AbstractArray one
copy(A::CatArray) = deepcopy(A)

## Categorical-specific methods

levels(A::CatArray) = levels(A.pool)

function _levels!(A::CatArray, newlevels::Vector; nullok=false)
    if !allunique(newlevels)
        throw(ArgumentError(string("duplicated levels found: ",
                                   join(unique(filter(x->sum(newlevels.==x)>1, newlevels)), ", "))))
    end

    # first pass to check whether changes can be applied without error
    # TODO: save original levels and undo changes in case of error to skip this step
    if !all(l->l in newlevels, index(A.pool))
        deleted = [!(l in newlevels) for l in index(A.pool)]
        @inbounds for (i, x) in enumerate(A.refs)
            if isa(A, CategoricalArray) && deleted[x]
                throw(ArgumentError("cannot remove level $(repr(index(A.pool)[x])) as it is used at position $i. Convert array to a Nullable$(typeof(A).name.name) if you want to transform some levels to missing values."))
            elseif isa(A, NullableCategoricalArray) && !nullok && deleted[x]
                throw(ArgumentError("cannot remove level $(repr(index(A.pool)[x])) as it is used at position $i and nullok=false."))
            end
        end
    end

    # actually apply changes
    oldindex = copy(index(A.pool))
    levels!(A.pool, newlevels)

    if index(A.pool) != oldindex
        # indexin returns 0 when not found, which maps to a missing value
        levelsmap = indexin(oldindex, index(A.pool))

        @inbounds for (i, x) in enumerate(A.refs)
            x > 0 && (A.refs[i] = levelsmap[x])
        end
    end

    levels(A.pool)
end

function droplevels!(A::CatArray)
    found = fill(false, length(index(A.pool)))
    @inbounds for i in A.refs
        i > 0 && (found[i] = true)
    end
    levels!(A, intersect(levels(A.pool), index(A.pool)[found]))
end

ordered(A::CatArray) = ordered(A.pool)
ordered!(A::CatArray, ordered) = ordered!(A.pool, ordered)

function Base.push!(A::CatArray, item)
    resize!(A.refs, length(A.refs) + 1)
    A[end] = item
    return A
end

function Base.append!(A::CatArray, B::CatArray)
    levels!(A, union(levels(A), levels(B)))
    len = length(A.refs)
    len2 = length(B.refs)
    resize!(A.refs, len + length(B.refs))
    for i = 1:len2
        A[len + i] = B[i]
    end
    return A
end

Base.empty!(A::CatArray) = (empty!(A.refs); return A)

categorical(A::AbstractArray; ordered=false) = CategoricalArray(A, ordered=ordered)
categorical{T<:Nullable}(A::AbstractArray{T}; ordered=false) =
    NullableCategoricalArray(A, ordered=ordered)

# Type-unstable methods
function categorical{T, N}(A::AbstractArray{T, N}, compact; ordered=false)
    RefType = compact ? reftype(length(unique(A))) : DefaultRefType
    CategoricalArray{T, N, RefType}(A, ordered=ordered)
end
function categorical{T<:Nullable, N}(A::AbstractArray{T, N}, compact; ordered=false)
    RefType = compact ? reftype(length(unique(A))) : DefaultRefType
    NullableCategoricalArray{T, N, RefType}(A, ordered=ordered)
end
function categorical{T, N, R}(A::CategoricalArray{T, N, R}, compact; ordered=false)
    RefType = compact ? reftype(length(levels(A))) : R
    CategoricalArray{T, N, RefType}(A, ordered=ordered)
end
function categorical{T, N, R}(A::NullableCategoricalArray{T, N, R}, compact; ordered=false)
    RefType = compact ? reftype(length(levels(A))) : R
    NullableCategoricalArray{T, N, RefType}(A, ordered=ordered)
end

## Code specific to CategoricalArray

function getindex(A::CategoricalArray, i::Int)
    j = A.refs[i]
    j > 0 || throw(UndefRefError())
    A.pool[j]
end

levels!(A::CategoricalArray, newlevels::Vector) = _levels!(A, newlevels)

function mergelevels(levels...)
    T = Base.promote_eltype(levels...)
    res = Array{T}(0)
    isordered = true

    for l in levels
        levelsmap = indexin(l, res)

        isordered &= issorted(levelsmap[levelsmap.!=0])
        if !isordered
            # Give up attempt to order res
            append!(res, l[levelsmap.==0])
        else
            i = length(res)+1
            for j = length(l):-1:1
                if levelsmap[j] == 0
                    insert!(res, i, l[j])
                else
                    i = levelsmap[j]
                end
            end
        end
    end

    res, isordered
end


text:
using EcologicalNetworks
using Test
using LinearAlgebra

my_tests = [
    "types/declaration.jl",
    "types/utilities.jl",
    "types/conversions.jl",
    "data.jl",
    "links/degree.jl",
    "links/specificity.jl",
    "links/connectance.jl",
    "links/degenerate.jl",
    "rand/draws.jl",
    "rand/shuffle.jl",
    "rand/null.jl",
    "rand/sample.jl",
    "rand/cascademodel.jl",
    "rand/nestedhierarchymodel.jl",
    "rand/nichemodel.jl",
    "rand/mpnmodel.jl",
    "community/nestedness.jl",
    "community/paths.jl",
    "community/overlap.jl",
    "community/centrality.jl",
    "community/motifs.jl",
    "community/foodwebs.jl",
    "community/resilience.jl",
    "community/spectralradius.jl",
    "betadiversity/operations.jl",
    "betadiversity/partitions.jl",
    "modularity/utilities.jl",
    "modularity/starters.jl",
    "modularity/labelpropagation.jl",
    "modularity/brim.jl",
    "modularity/roles.jl",
    "information/entropy.jl",
    "information/otsin.jl"
]
# "mangal.jl",

global test_n
global anyerrors

test_n = 1
anyerrors = false

for my_test in my_tests
  try
    include(my_test)
    println("[TEST $(lpad(test_n,2))] \033[1m\033[32mPASS\033[0m $(my_test)")
  catch e
    global anyerrors = true
    println("[TEST $(lpad(test_n,2))] \033[1m\033[31mFAIL\033[0m $(my_test)")
    showerror(stdout, e, backtrace())
    println()
    throw("TEST FAILED")
  end
  global test_n += 1
end

if anyerrors
  throw("Tests failed")
end


text:
# common functions for ESPPRC
# pulse.jl and labeling.jl

# Main ESPPRC function
function solveESPPRC(org_pg::ESPPRC_Instance; method="pulse", DSSR=false)

    if method == "pulse"
        best_label, _ = solveESPPRCpulse(org_pg)
    elseif method == "monodirectional"
        best_label, _ = monodirectional(org_pg; DSSR=DSSR)
    elseif method == "bidirectional"
        best_label, _ = bidirectional(org_pg; DSSR=DSSR)
    else
        best_label, _ = solveESPPRCpulse(org_pg)
    end

    return best_label
end

function solveESPPRC_vrp(org_pg::ESPPRC_Instance; max_neg_routes=MAX_INT, method="pulse", DSSR=false)

    if DSSR == true && max_neg_routes < MAX_INT
        @warn("DSSR option cannot be used with max_neg_routes. Non-DSSR algorithm is used.")
        DSSR = false
    end

    if method == "pulse"
        return solveESPPRCpulse(org_pg; max_neg_routes=max_neg_routes)
    elseif method == "monodirectional"
        return monodirectional(org_pg; max_neg_routes=max_neg_routes, DSSR=DSSR)
    elseif method == "bidirectional"
        return bidirectional(org_pg; max_neg_routes=max_neg_routes, DSSR=DSSR)
    else
        return solveESPPRCpulse(org_pg; max_neg_routes=max_neg_routes)
    end
end


text:
using JPkg
using Test

@testset "JPkg.jl" begin

    @test 1==1
    # Write your own tests here.
end


text:
@with_kw struct ContinuousLightDark <: POMDPs.POMDP{Float64,Int,Float64}
    discount_factor::Float64        = 0.95
    correct_r::Float64              = 100.0
    incorrect_r::Float64            = -100.0
    light_loc::Float64              = 10.0
    big_step::Int                   = 3
    movement_cost::Float64          = -1.0
    transition_std::Float64         = 1.0
    observation_std::Float64        = sqrt(2)
    observation_std_base::Float64   = 0.5
    goal_radius::Float64            = 1.0
    initialstate_distribution      = Normal(0.0, 10.0)
end
discount(p::ContinuousLightDark) = p.discount_factor
isterminal_act(::ContinuousLightDark, a::Int) = a == 0

actions(p::ContinuousLightDark) = [-p.big_step, -1, 0, 1, p.big_step]
n_actions(p::ContinuousLightDark) = length(actions(p))
actionindex(p::ContinuousLightDark, a::Int) = findfirst(actions(p), a)

transition(p::ContinuousLightDark, s::Float64, a::Int) = Normal(s+a, p.transition_std)

function observation(p::ContinuousLightDark, sp::Float64)
    Normal(sp, p.observation_std*abs(sp-p.light_loc) + p.observation_std_base)
end

function reward(p::ContinuousLightDark, s::Float64, a::Int, sp::Float64)
    if a == 0
        if abs(s) < p.goal_radius
            return p.correct_r
        elseif abs(s) < 1.5 * p.goal_radius
            return p.correct_r + (p.correct_r-p.incorrect_r)/2 * (p.goal_radius - abs(s))
        elseif abs(s) < 2.0 * p.goal_radius
            return p.incorrect_r + (p.correct_r-p.incorrect_r)*3/2 * (2.0 * p.goal_radius - abs(s))
        else
            return p.incorrect_r
        end
    else
        return p.movement_cost
    end
end

function initialstate_distribution(pomdp::ContinuousLightDark)
    return pomdp.initialstate_distribution
end

# Added to use PFT algorithm

sampletype(d::Normal{T}) where {T} = T

function max_possible_weight(pomdp::ContinuousLightDark, a, o)
    return obs_weight(pomdp, o, o)
end

function new_particle(p::ContinuousLightDark, b, a, o, rng)
    return rand(rng, observation(p, o))
end

max_entropy_distribution(pomdp::ContinuousLightDark) = Normal(0., 100.)

# For LocalApproximationValueIteration
function POMDPs.convert_s(::Type{V} where V <: AbstractVector{Float64}, s::Float64, pomdp::ContinuousLightDark)
    v = SVector{1,Float64}(s)
end

function POMDPs.convert_s(::Type{Float64}, v::AbstractVector{Float64}, pomdp::ContinuousLightDark)
    s = v[1]
end

# Walk towards a goal with big and small steps

mutable struct CLDGoalPolicy <: Policy
    pomdp::ContinuousLightDark
    goal::Float64
end

function action(p::CLDGoalPolicy, b)
    m = mean(b)
    small_step = 1
    if abs(m - p.goal) >= (p.pomdp.big_step + 1)/2
        return -p.pomdp.big_step*Int(sign(m-p.goal))
    else
        return -small_step*Int(sign(m-p.goal))
    end
end

# Simple Policy

@with_kw mutable struct CLDPingPong <: Solver
    turnpoint::Float64  = 10.
    direction::Int      = 1
end

struct CLDPingPongPolicy <: Policy
    s::CLDPingPong
end

solve(solver::CLDPingPong, pomdp::ContinuousLightDark) = CLDPingPongPolicy(solver)
Random.seed!(p::CLDPingPongPolicy, s) = p

function action(policy::CLDPingPongPolicy, b)
    if mean(b) > policy.s.turnpoint
        policy.s.direction = -1
    elseif mean(b) < -policy.s.turnpoint
        policy.s.direction = 1
    end
    policy.s.direction
end

# Information greedy Policy

mutable struct CLDInformationGain <: Solver end

solve(solver::CLDInformationGain, pomdp::ContinuousLightDark) = CLDGoalPolicy(pomdp, pomdp.light_loc)
Random.seed!(p::CLDGoalPolicy, s) = p

# Heuristic policy based on standard deviation

mutable struct CLDHeuristic <: Policy
    p::ContinuousLightDark
    policy::CLDGoalPolicy
    std_thresh::Float64
    localized::Bool
end

solve(sol::LDHSolver, pomdp::ContinuousLightDark) = CLDHeuristic(pomdp, CLDGoalPolicy(pomdp, 0.0), sol.std_thresh, false)

function action(p::CLDHeuristic, b::AbstractParticleBelief)
    big_step = p.p.big_step
    small_step = 1
    s = std(b)
    m = mean(b)

    if !p.localized && s <= p.std_thresh
        p.localized = true
    end

    if p.localized
        # If localized set origin as goal
        p.policy.goal = 0.0
        # If close enough tag goal with action 0
        if abs(m - p.policy.goal) < small_step
            return 0
        end
    else
        # Else set light source as goal
        p.policy.goal = p.p.light_loc
    end
    # Move to goal if not close enough
    return action(p.policy, b)
end

action(p::CLDHeuristic, s::Float64) = (p.policy.goal = 0.0; action(p.policy, ParticleCollection([s])))

Base.copy(p::CLDHeuristic) = CLDHeuristic(p.p, p.policy, p.std_thresh, p.localized)


text:
using Documenter, ParametricOptInterface

makedocs(sitename="My Documentation")


text:
module EventL1

    include("../RawFileofflineProc/BaseLib.jl")
    include( "LogLstLib.jl")
    using ..SegParaAlarmProc
    function ReportEventL1(s_1::UInt,s_2::UInt, curTime::UInt,curFrame ,filterinfo="NoWind")
        posStr = @sprintf("%04u:%04u,%04u", s_1,s_2,s_2-s_1)
        a = ("intrusion",posStr, curTime , UInt2DTStr( curTime ),curFrame,filterinfo )
        a1= join(a,",")*"\r\n"
        LogEvent(a1)
    end

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
"""
The convolution of a one galaxy component with one PSF component.
It also contains the derivatives of sigma with respect to the shape parameters.
It does not contain the derivatives with respect to other parameters
(pos and gal_frac_dev) because they have easy expressions in terms of other known
quantities.

Args:
 - gal_frac_dev_dir: "Theta direction": this is 1 or -1, depending on whether
     increasing gal_frac_dev increases the weight of this GalaxyCacheComponent
     (1) or decreases it (-1).
 - gal_frac_dev_i: The weight given to this type of galaxy for this celestial object.
     This is either gal_frac_dev or (1 - gal_frac_dev).
 - gc: The galaxy component to be convolved
 - pc: The psf component to be convolved
 - pos: The location of the celestial object in pixel coordinates as a 2x1 vector
 - gal_axis_ratio: The ratio of the galaxy minor axis to major axis (0 < gal_axis_ratio <= 1)
 - gal_radius_px: The scale of the galaxy major axis

Attributes:
 - gal_frac_dev_dir: Same as input
 - gal_frac_dev_i: Same as input
 - bmc: A BvnComponent with the convolution.
 - dSigma: A 3x3 matrix containing the derivates of
     [Sigma11, Sigma12, Sigma22] (in the rows) with respect to
     [gal_axis_ratio, gal_angle, gal_radius_px] (in the columns)
"""
struct GalaxyCacheComponent{T<:Number}
    gal_frac_dev_dir::Float64
    gal_frac_dev_i::T
    bmc::BvnComponent{T}
    sig_sf::GalaxySigmaDerivs{T}
    # [Sigma11, Sigma12, Sigma22] x [gal_axis_ratio, gal_angle, gal_radius_px]
end

function GalaxyCacheComponent(
        gal_frac_dev_dir::Float64,
        gal_frac_dev_i::T,
        gc::GalaxyComponent,
        pc::PsfComponent,
        pos::Vector{T},
        gal_axis_ratio::T,
        gal_angle::T,
        gal_radius_px::T,
        calculate_gradient::Bool,
        calculate_hessian::Bool) where {T<:Number}

    XiXi = get_bvn_cov(gal_axis_ratio, gal_angle, gal_radius_px)
    mean_s = @SVector T[pc.xiBar[1] + pos[1], pc.xiBar[2] + pos[2]]
    var_s = pc.tauBar + gc.nuBar * XiXi
    weight = pc.alphaBar * gc.etaBar  # excludes gal_frac_dev

    # d siginv / dsigma is only necessary for the Hessian.
    bmc = BvnComponent(mean_s, var_s, weight, calculate_gradient && calculate_hessian)

    if calculate_gradient
        sig_sf = GalaxySigmaDerivs(gal_angle, gal_axis_ratio, gal_radius_px, XiXi,
                                   gc.nuBar, calculate_hessian)
    else
        sig_sf = GalaxySigmaDerivs(T)
    end

    GalaxyCacheComponent(gal_frac_dev_dir, gal_frac_dev_i, bmc, sig_sf)
end

struct BvnBundle{T<:Real}
    bvn_derivs::BivariateNormalDerivatives{T}
    star_mcs::Matrix{BvnComponent{T}}
    gal_mcs::Array{GalaxyCacheComponent{T},4}
    function (::Type{BvnBundle{T}}){T}(psf_K::Int, S::Int)
        return new{T}(BivariateNormalDerivatives{T}(),
                      Matrix{BvnComponent{T}}(psf_K, S),
                      Array{GalaxyCacheComponent{T}}(psf_K, 8, 2, S))
    end
end

function zero!(bvn_bundle::BvnBundle)
    BivariateNormals.zero!(bvn_bundle.bvn_derivs)
    return bvn_bundle
end

"""
Convolve the current locations and galaxy shapes with the PSF.  If
calculate_gradient is true, also calculate derivatives and hessians for
active sources.

Args:
 - S (formerly from ElboArgs)
 - patches (formerly from ElboArgs)
 - vp: (formerly from ElboArgs)
 - active_sources (formerly from ElboArgs)
 - psf_K: Number of psf components (psf from patches object)
 - n: the image id (not the band)
 - calculate_gradient: Whether to calculate derivatives for active sources.
 - calculate_hessian

Returns:
 - star_mcs: An array of BvnComponents with indices
    - PSF component
    - Source (index within active_sources)
 - gal_mcs: An array of BvnComponents with indices
    - PSF component
    - Galaxy component
    - Galaxy type
    - Source (index within active_sources)
  Hessians are only populated for s in ea.active_sources.
"""
function load_bvn_mixtures!(
                    #outputs
                    star_mcs::Matrix{BvnComponent{T}},
                    gal_mcs::Array{GalaxyCacheComponent{T},4},
                    #inputs
                    S::Int64,
                    patches::Matrix{ImagePatch},
                    source_params::Vector{Vector{T}},
                    active_sources::Vector{Int},
                    psf_K::Int64,
                    n::Int,
                    calculate_gradient::Bool=true,
                    calculate_hessian::Bool=true) where {T<:Number}
    @assert size(star_mcs, 1) == psf_K
    @assert size(star_mcs, 2) == S
    # TODO: do not keep any derviative information if the sources are not in
    # active_sources.

    for s in 1:S
        psf = patches[s, n].psf
        sp  = source_params[s]

        # TODO: it's a lucky coincidence that lidx.pos = ids.pos.
        # That's why this works when called from both log_prob.jl with `sp`
        # and elbo_objective.jl with `vp`.
        # We need a safer way to let both methods call this method.
        world_loc = sp[lidx.pos]
        m_pos = Model.linear_world_to_pix(patches[s, n].wcs_jacobian,
                                          patches[s, n].world_center,
                                          patches[s, n].pixel_center, world_loc)

        # Convolve the star locations with the PSF.
        for k in 1:psf_K
            pc = psf[k]
            mean_s = @SVector T[pc.xiBar[1] + m_pos[1], pc.xiBar[2] + m_pos[2]]
            star_mcs[k, s] = BvnComponent(mean_s, pc.tauBar, pc.alphaBar, false)
        end

        # Convolve the galaxy representations with the PSF.
        for i = 1:2 # i indexes dev vs exp galaxy types.
            gal_frac_dev_dir = (i == 1) ? 1. : -1.
            gal_frac_dev_i = (i == 1) ? sp[lidx.gal_frac_dev] : 1. - sp[lidx.gal_frac_dev]

            # Galaxies of type 1 have 8 components, and type 2 have 6 components.
            for j in 1:ifelse(i == 1, 8, 6)
                for k = 1:psf_K
                    gal_mcs[k, j, i, s] = GalaxyCacheComponent(
                        gal_frac_dev_dir, gal_frac_dev_i, galaxy_prototypes[i][j], psf[k],
                        m_pos,
                        sp[lidx.gal_axis_ratio], sp[lidx.gal_angle], sp[lidx.gal_radius_px],
                        calculate_gradient && (s in active_sources),
                        calculate_hessian)
                end
            end
        end
    end

    star_mcs, gal_mcs
end

function load_bvn_mixtures(
        S::Int64,
        patches::Matrix{ImagePatch},
        source_params::Vector{Vector{T}},
        active_sources::Vector{Int},
        psf_K::Int64,
        n::Int;
        calculate_gradient::Bool=true,
        calculate_hessian::Bool=true) where {T<:Number}
    star_mcs = Matrix{BvnComponent{T}}(psf_K, S)
    gal_mcs  = Array{GalaxyCacheComponent{T}}(psf_K, 8, 2, S)
    load_bvn_mixtures!(star_mcs, gal_mcs, S, patches, source_params,
                       active_sources, psf_K, n, calculate_gradient,
                       calculate_hessian)
end

"""
Populate fs0m and fs1m for source s in the a given pixel.

Non-standard args:
    - x: The pixel location in the image
"""
function populate_gal_fsm!(
        fs1m::SensitiveFloat{T},
        bvn_derivs::BivariateNormalDerivatives{T},
        s::Int,
        h::Int64,
        w::Int64,
        is_active_source::Bool,
        wcs_jacobian::Matrix{Float64},
        gal_mcs::Array{GalaxyCacheComponent{T}, 4}) where {T<:Number}
    SensitiveFloats.zero!(fs1m)
    x = SVector{2,Float64}(h, w)

    for i = 1:2 # Galaxy types
        for j in 1:8 # Galaxy component
            # If i == 2 then there are only six galaxy components.
            if (i == 1) || (j <= 6)
                for k = 1:size(gal_mcs, 1) # PSF component
                    accum_galaxy_pos!(
                        fs1m, bvn_derivs,
                        gal_mcs[k, j, i, s], x, wcs_jacobian,
                        is_active_source)
                end
            end
        end
    end
end

softpluslike(x) = 1000x > 1 ? 1000x - 1 : log(1000x)
softpluslikeinv(y) = y < 0 ? 1e-3exp(y) : 1e-3(y + 1)

function star_light_density!(fs0m::SensitiveFloat{T},
                             p::ImagePatch,
                             h::Int,
                             w::Int,
                             pos::Vector,
                             is_active_source::Bool) where {T <: Number}
     f(local_pos::Vector) = begin
         m_pos = Model.linear_world_to_pix(p.wcs_jacobian,
                                           p.world_center,
                                           p.pixel_center,
                                           local_pos)
         softpluslikeinv(p.itp_psf[h - m_pos[1] + 26, w - m_pos[2] + 26])
     end

     fs0m.v[] = f(pos)

     if fs0m.has_gradient && is_active_source
         ForwardDiff.gradient!(fs0m.d, f, pos)
     end

     if fs0m.has_hessian && is_active_source
         ForwardDiff.hessian!(fs0m.h, f, pos)
     end
end

"""
Add the contributions of a galaxy component term to the ELBO by
updating fs1m in place.
"""
function accum_galaxy_pos!(fs1m::SensitiveFloat{T},
                           bvn_derivs::BivariateNormalDerivatives{T},
                           gcc::GalaxyCacheComponent{T},
                           x::SVector{2,Float64},
                           wcs_jacobian::Array{Float64, 2},
                           is_active_source::Bool) where {T<:Number}
    eval_bvn_pdf!(bvn_derivs, gcc.bmc, x)
    f = bvn_derivs.f_pre[1] * gcc.gal_frac_dev_i
    fs1m.v[] += f

    if fs1m.has_gradient && is_active_source
        get_bvn_derivs!(bvn_derivs, gcc.bmc,
                        fs1m.has_gradient, fs1m.has_hessian)
        transform_bvn_derivs!(
            bvn_derivs, gcc.sig_sf, wcs_jacobian, fs1m.has_hessian)

        @aliasscope begin
            bvn_u_d = Const(bvn_derivs.bvn_u_d)
            bvn_s_d = Const(bvn_derivs.bvn_s_d)
            bvn_derivs_f_pre = Const(bvn_derivs.f_pre)

            # Accumulate the derivatives.
            @inbounds for u_id in 1:2
                fs1m.d[gal_ids.pos[u_id]] += f * bvn_u_d[u_id]
            end

            @inbounds for gal_id in 1:length(gal_shape_ids)
                fs1m.d[gal_shape_alignment[gal_id]] += f * bvn_s_d[gal_id]
            end

            # The gal_frac_dev derivative. gal_frac_dev just scales the entire component.
            # The direction is positive or negative depending on whether this
            # is an exp or dev component.
            @inbounds fs1m.d[gal_ids.gal_frac_dev] += gcc.gal_frac_dev_dir * bvn_derivs_f_pre[1]

            if fs1m.has_hessian
                # The Hessians:
                bvn_uu_h = Const(bvn_derivs.bvn_uu_h)
                bvn_ss_h = Const(bvn_derivs.bvn_ss_h)
                bvn_us_h = Const(bvn_derivs.bvn_us_h)
                gal_ids_u = Const(gal_ids.pos)

                # Hessian terms involving only the shape parameters.
                @inbounds for shape_id1 in 1:length(gal_shape_ids)
                    @inbounds for shape_id2 in 1:length(gal_shape_ids)
                        s1 = gal_shape_alignment[shape_id1]
                        s2 = gal_shape_alignment[shape_id2]
                        fs1m.h[s1, s2] +=
                            f * (bvn_ss_h[shape_id1, shape_id2] +
                                     bvn_s_d[shape_id1] * bvn_s_d[shape_id2])
                    end
                end

                # Hessian terms involving only the location parameters.
                @inbounds for u_id1 in 1:2
                    @inbounds for u_id2 in 1:2
                        u1 = gal_ids_u[u_id1]
                        u2 = gal_ids_u[u_id2]
                        fs1m.h[u1, u2] +=
                            f * (bvn_uu_h[u_id1, u_id2] + bvn_u_d[u_id1] * bvn_u_d[u_id2])
                    end
                end

                # Hessian terms involving both the shape and location parameters.
                for u_id in 1:2
                  @inbounds for shape_id in 1:length(gal_shape_ids)
                    ui = gal_ids_u[u_id]
                    si = gal_shape_alignment[shape_id]
                    fs1m.h[ui, si] +=
                        f * (bvn_us_h[u_id, shape_id] + bvn_u_d[u_id] * bvn_s_d[shape_id])
                    fs1m.h[si, ui] = fs1m.h[ui, si]
                  end
                end

                # Do the gal_frac_dev hessian terms.
                devi = gal_ids.gal_frac_dev
                @inbounds for u_id in 1:2
                    ui = gal_ids.pos[u_id]
                    fs1m.h[ui, devi] +=
                        bvn_derivs_f_pre[1] * gcc.gal_frac_dev_dir * bvn_u_d[u_id]
                    fs1m.h[devi, ui] = fs1m.h[ui, devi]
                end
                @inbounds for shape_id in 1:length(gal_shape_ids)
                    si = gal_shape_alignment[shape_id]
                    fs1m.h[si, devi] +=
                        bvn_derivs_f_pre[1] * gcc.gal_frac_dev_dir * bvn_s_d[shape_id]
                    fs1m.h[devi, si] = fs1m.h[si, devi]
                end
            end # if calculate hessian
        end
    end # if is_active_source
end

function write_star_nmgy!(world_pos::Array{Float64,1},
                          flux::Float64,
                          patch::ImagePatch,
                          pixels::Matrix{Float32};
                          write_to_patch::Bool=false)
    fs0m = SensitiveFloat{Float64}(length(StarPosParams), 1, false, false)
    H2, W2 = size(patch.active_pixel_bitmap)
    for w2 in 1:W2, h2 in 1:H2
        h = patch.bitmap_offset[1] + h2
        w = patch.bitmap_offset[2] + w2
        Model.star_light_density!(fs0m, patch, h, w, world_pos, false)
        if write_to_patch
            pixels[h2, w2] += fs0m.v[] * flux
        else
            pixels[h, w] += fs0m.v[] * flux
        end
    end
end

function write_galaxy_nmgy!(world_pos::Array{Float64,1},
                            flux::Float64,
                            gal_frac_dev::Float64,
                            gal_axis_ratio::Float64,
                            gal_angle::Float64,
                            gal_radius_px::Float64,
                            psf::Array{Model.PsfComponent,1},
                            patches::Array{ImagePatch, 2},
                            pixels::Matrix{Float32};
                            write_to_patch::Bool=false)
    bvn_derivs = Model.BivariateNormalDerivatives{Float64}()
    fs1m = SensitiveFloat{Float64}(length(GalaxyPosParams), 1, false, false)
    source_params = [[world_pos[1], world_pos[2], gal_frac_dev, gal_axis_ratio,
                     gal_angle, gal_radius_px],]
    star_mcs, gal_mcs = Model.load_bvn_mixtures(1, patches,
                          source_params, [1,], length(psf), 1,
                          calculate_gradient=false,
                          calculate_hessian=false)
    p = patches[1]
    H2, W2 = size(p.active_pixel_bitmap)
    for w2 in 1:W2, h2 in 1:H2
        # (h2, w2) index the local patch, while (h, w) index the image
        h = p.bitmap_offset[1] + h2
        w = p.bitmap_offset[2] + w2
        Model.populate_gal_fsm!(fs1m, bvn_derivs, 1, h, w, false, p.wcs_jacobian, gal_mcs)
        if write_to_patch
            pixels[h2, w2] += fs1m.v[] * flux
        else
            pixels[h, w] += fs1m.v[] * flux
        end
    end
end


text:
d = dataeye(3, 3)
d[1, 1] = NA

svd(d)


text:
# this is just a helper module to test cross-module fixtures
module different_module
  using PyTest

  # FIXME: can we get rid of this?
  using ResumableFunctions

  export different_module_f, different_module_g,
         some_function

  @fixture function different_module_f() :dmf_result end
  @fixture function different_module_g(different_module_f)
    return [different_module_f, :dmg_result]
  end

  function some_function() :some_function_result end
end


text:
"""
    InitialValueProblem

Parametric composite type for initial value problems. It is parameterized in the
system's type.

### Fields

- `s`  -- system
- `x0` -- initial state

### Examples

The linear system ``x' = -x`` with initial condition ``x_0 = [-1/2, 1/2]``:

```jldoctest
julia> p = InitialValueProblem(LinearContinuousSystem([-1.0 0.0; 0.0 -1.0]), [-1/2, 1/2]);

julia> p.x0
2-element Array{Float64,1}:
 -0.5
  0.5

julia> statedim(p)
2

julia> inputdim(p)
0
```
"""
struct InitialValueProblem{S <: AbstractSystem, XT} <: AbstractSystem
    s::S
    x0::XT
end
statedim(ivp::InitialValueProblem) = statedim(ivp.s)
inputdim(ivp::InitialValueProblem) = inputdim(ivp.s)
inputset(ivp::InitialValueProblem) = inputset(ivp.s)
islinear(ivp::InitialValueProblem) = islinear(ivp.s)
isaffine(ivp::InitialValueProblem) = isaffine(ivp.s)

"""
    IVP

`IVP` is an alias for `InitialValueProblem`.
"""
const IVP = InitialValueProblem


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7]);F = [2824*x1^4+7500*x1^3*x3+288*x1^3*x5-1100*x1^2*x2*x5+5625*x1^2*x3^2-1332*x1^2*x3*x6-7744*x1^2*x4^2+352*x1^2*x4*x6-1400*x1^2*x4*x7+3520*x1^2*x5^2+2952*x1^2*x5*x6+125*x1^2*x6^2+5600*x1^2*x6*x7-1650*x1*x2*x3*x5-16*x1*x2*x3*x6-3008*x1*x2*x4*x6-2048*x1*x2*x6^2-160*x1*x3^2*x4-110*x1*x3^2*x6-12000*x1*x3*x4^2-2100*x1*x3*x4*x7-592*x1*x3*x5*x6+8400*x1*x3*x6*x7-2080*x1*x4^3-1430*x1*x4^2*x6+1536*x1*x5^3+1312*x1*x5^2*x6+16*x2^2*x3^2-40*x2^2*x3*x6+121*x2^2*x5^2+8861*x2^2*x6^2+940*x2*x3^2*x6+1760*x2*x4^2*x5+12220*x2*x4^2*x6+308*x2*x4*x5*x7-1232*x2*x5*x6*x7+25*x3^4+650*x3^2*x4^2+1369*x3^2*x6^2-7104*x3*x5^2*x6-6068*x3*x5*x6^2+10625*x4^4+2240*x4^3*x7-8960*x4^2*x6*x7+196*x4^2*x7^2-1568*x4*x6*x7^2+9216*x5^4+15744*x5^3*x6+6724*x5^2*x6^2+3136*x6^2*x7^2-864*x1^2*x4+28*x1^2*x6-56*x1*x2*x3+26*x1*x2*x6+2304*x1*x3*x4+1584*x1*x3*x6-384*x1*x4*x5+88*x2^2*x3-110*x2^2*x6-13536*x2*x3*x6-720*x3^3-9360*x3*x4^2+1776*x3*x4*x6-4608*x4*x5^2-3936*x4*x5*x6+49*x1^2-154*x1*x2+344*x1*x6+121*x2^2-688*x2*x3+860*x2*x6+5184*x3^2+576*x4^2+1204*x1-1892*x2+7396, -1100*x1^2*x5-1650*x1*x3*x5-16*x1*x3*x6-3008*x1*x4*x6-2048*x1*x6^2+32*x2*x3^2-80*x2*x3*x6+242*x2*x5^2+17722*x2*x6^2+940*x3^2*x6+1760*x4^2*x5+12220*x4^2*x6+308*x4*x5*x7-1232*x5*x6*x7-56*x1*x3+26*x1*x6+176*x2*x3-220*x2*x6-13536*x3*x6-154*x1+242*x2-688*x3+860*x6-1892, 7500*x1^3+11250*x1^2*x3-1332*x1^2*x6-1650*x1*x2*x5-16*x1*x2*x6-320*x1*x3*x4-220*x1*x3*x6-12000*x1*x4^2-2100*x1*x4*x7-592*x1*x5*x6+8400*x1*x6*x7+32*x2^2*x3-40*x2^2*x6+1880*x2*x3*x6+100*x3^3+1300*x3*x4^2+2738*x3*x6^2-7104*x5^2*x6-6068*x5*x6^2-56*x1*x2+2304*x1*x4+1584*x1*x6+88*x2^2-13536*x2*x6-2160*x3^2-9360*x4^2+1776*x4*x6-688*x2+10368*x3, -15488*x1^2*x4+352*x1^2*x6-1400*x1^2*x7-3008*x1*x2*x6-160*x1*x3^2-24000*x1*x3*x4-2100*x1*x3*x7-6240*x1*x4^2-2860*x1*x4*x6+3520*x2*x4*x5+24440*x2*x4*x6+308*x2*x5*x7+1300*x3^2*x4+42500*x4^3+6720*x4^2*x7-17920*x4*x6*x7+392*x4*x7^2-1568*x6*x7^2-864*x1^2+2304*x1*x3-384*x1*x5-18720*x3*x4+1776*x3*x6-4608*x5^2-3936*x5*x6+1152*x4, 288*x1^3-1100*x1^2*x2+7040*x1^2*x5+2952*x1^2*x6-1650*x1*x2*x3-592*x1*x3*x6+4608*x1*x5^2+2624*x1*x5*x6+242*x2^2*x5+1760*x2*x4^2+308*x2*x4*x7-1232*x2*x6*x7-14208*x3*x5*x6-6068*x3*x6^2+36864*x5^3+47232*x5^2*x6+13448*x5*x6^2-384*x1*x4-9216*x4*x5-3936*x4*x6, -1332*x1^2*x3+352*x1^2*x4+2952*x1^2*x5+250*x1^2*x6+5600*x1^2*x7-16*x1*x2*x3-3008*x1*x2*x4-4096*x1*x2*x6-110*x1*x3^2-592*x1*x3*x5+8400*x1*x3*x7-1430*x1*x4^2+1312*x1*x5^2-40*x2^2*x3+17722*x2^2*x6+940*x2*x3^2+12220*x2*x4^2-1232*x2*x5*x7+2738*x3^2*x6-7104*x3*x5^2-12136*x3*x5*x6-8960*x4^2*x7-1568*x4*x7^2+15744*x5^3+13448*x5^2*x6+6272*x6*x7^2+28*x1^2+26*x1*x2+1584*x1*x3-110*x2^2-13536*x2*x3+1776*x3*x4-3936*x4*x5+344*x1+860*x2, -1400*x1^2*x4+5600*x1^2*x6-2100*x1*x3*x4+8400*x1*x3*x6+308*x2*x4*x5-1232*x2*x5*x6+2240*x4^3-8960*x4^2*x6+392*x4^2*x7-3136*x4*x6*x7+6272*x6^2*x7];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeWatchedFiles")}}, params) = DidChangeWatchedFilesParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeWatchedFiles")},DidChangeWatchedFilesParams}, server)
    for change in r.params.changes
        uri = change.uri
        !haskey(server.documents, URI2(uri)) && continue
        if change.type == FileChangeTypes["Created"] || (change.type == FileChangeTypes["Changed"] && !get_open_in_editor(server.documents[URI2(uri)]))
            doc = server.documents[URI2(uri)]
            filepath = uri2filepath(uri)
            content = String(read(filepath))
            content == get_text(doc) && return
            server.documents[URI2(uri)] = Document(uri, content, true, server)
            parse_all(server.documents[URI2(uri)], server)

        elseif change.type == FileChangeTypes["Deleted"] && !get_open_in_editor(server.documents[URI2(uri)])
            delete!(server.documents, URI2(uri))

            response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(uri, Diagnostic[]))
            send(response, server)
        end
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeConfiguration")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeConfiguration")},Dict{String,Any}}, server::LanguageServerInstance)
    send(JSONRPC.Request{Val{Symbol("workspace/configuration")},ConfigurationParams}(-100, ConfigurationParams([
        (ConfigurationItem(missing, "julia.format.$opt") for opt in fieldnames(DocumentFormat.FormatOptions))...;
        ConfigurationItem(missing, "julia.runLinter")
        ])), server)
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeWorkspaceFolders")}}, params) = didChangeWorkspaceFoldersParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeWorkspaceFolders")}}, server)
    for wksp in r.params.event.added
        push!(server.workspaceFolders, uri2filepath(wksp.uri))
        load_folder(wksp, server)
    end
    for wksp in r.params.event.removed
        delete!(server.workspaceFolders, uri2filepath(wksp.uri))
        remove_workspace_files(wksp, server)
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/symbol")}}, params) = WorkspaceSymbolParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/symbol")},WorkspaceSymbolParams}, server)
    syms = SymbolInformation[]
    for (uri,doc) in server.documents
        bs = collect_toplevel_bindings_w_loc(getcst(doc), query = r.params.query)
        for x in bs
            p, b = x[1], x[2]
            push!(syms, SymbolInformation(valof(b.name), 1, false, Location(doc._uri, Range(doc, p)), nothing))
        end
    end

    response = JSONRPC.Response(r.id, syms)
    send(response, server)
end


text:
using Documenter, Mera

makedocs(modules = [Mera],
         sitename = "Mera.jl",
         doctest = false,
         clean = true,
         checkdocs = :all,
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true", sidebar_sitename = false),
         authors = "Manuel Behrendt",
		 pages = Any[ "Home"                  => "index.md",
		              "First Steps"           => "00_multi_FirstSteps.md",
                      "1-Data Inspection"     => Any[ "Hydro"      =>  "01_hydro_First_Inspection.md",
                                                      "Particles"  =>  "01_particles_First_Inspection.md",
                                                      "Clumps"     =>  "01_clumps_First_Inspection.md"],

                      "2-Load by Selection"   => Any[ "Hydro"      =>  "02_hydro_Load_Selections.md",
                                                      "Particles"  =>  "02_particles_Load_Selections.md",
                                                      "Clumps"     =>  "02_clumps_Load_Selections.md"],

                       "3-Get Subregions"     => Any[ "Hydro"      => "03_hydro_Get_Subregions.md",
                                                      "Particles"  => "03_particles_Get_Subregions.md",
                                                      "Clumps"     => "03_clumps_Get_Subregions.md"],

                      "4-Basic Calculations"  => "04_multi_Basic_Calculations.md",

                      "5-Mask/Filter/Meta"    => "05_multi_Masking_Filtering.md",

                      "6-Projection"          => Any[ "Hydro"      =>  "06_hydro_Projection.md",
                                                      "Particles"  =>  "06_particles_Projection.md"],

                      "Examples"              => "examples.md",

                      "API Documentation"     => "api.md"
                      ]
		 )

deploydocs(repo = "github.com/ManuelBehrendt/Mera.jl.git")


text:
using Distributions, StatsBase, Random
Random.seed!(1)

dist1 = TriangularDist(0,10,5)
dist2 = DiscreteUniform(1,5)
theorMean1, theorMean2 = mean(dist1), mean(dist2)

N = 10^6
data1 = rand(dist1,N)
data2 = rand(dist2,N)
estMean1, estMean2 = mean(data1), mean(data2)

println("Symmetric Triangular Distiribution on [0,10] has mean $theorMean1
	(estimated: $estMean1)")
println("Discrete Uniform Distiribution on {1,2,3,4,5} has mean $theorMean2
	(estimated: $estMean2)")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Currencies
using Test

currencies = ((:USD, 2, 840, "US Dollar"),
              (:EUR, 2, 978, "Euro"),
              (:JPY, 0, 392, "Yen"),
              (:JOD, 3, 400, "Jordanian Dinar"),
              (:CNY, 2, 156, "Yuan Renminbi"))

# This just makes sure that the data was loaded and at least some basic values are set as expected
@testset "Basic currencies" begin
    for (s, u, c, n) in currencies
        ccy = Currency{s}
        @test currency(s) == ccy
        @test symbol(ccy) == s
        @test unit(ccy) == u
        @test name(ccy) == n
        @test code(ccy) == c
    end
end

# This makes sure that the values are within expected ranges
@testset "Validation" begin
    @test length(Currencies.allpairs()) >= 155
    for (s, (ct,u,c,n)) in Currencies.allpairs()
        ccy = Currency(s)
        @test symbol(ct) == s
        @test length(string(s)) == 3
        @test u >= 0
        @test 0 < c < 2000
        @test length(n) < 40
    end
end

max_stars_repo_name: TotalVerb/Currencies

max_issues_repo_name: TotalVerb/Currencies

max_forks_repo_name: TotalVerb/Currencies


text:
module SeisCore

# packages
using Base, Core, SeisIO, SeisNoise, Dates, CSV, DataFrames, SCEDC, AWS, Distributed, Parallelism, AWSS3, Glob,
        Statistics, JLD2, HDF5, AbstractFFTs, Plots

# functions
include("download.jl")
include("process.jl")
include("correlate.jl")
include("big_correlate.jl")
end


text:
## extend logical operators so that we can have a notation like

#f(x,y) < 0 or f(x,y) == 0

import Base: <, <=, ==, !=, !==, >=, >, &, |, !

abstract Predicate

"""

A predicate is defined in terms of a function of two variables, an
inquality, and either another function or a real number.  For example,
`f < 0` or `f >= g`. The case `f==g` and `f != g` are not defined, as doing so crosses
up `Gadfly` and other code that compares functions for equality. Use
`eq(f,g)`  or `f \Equal<tab> g` for equality and `neq(f,g)` or `f \gtrless<tab> g` for not equal.

Available operations to produce predicates:

* `<`
* `<=`, `≤` (`\le<tab>`)
* `==` (`Function == Real`), `eq(f,g)`, `⩵` (`\Equal<tab>`)
* `!==` (`Function != Real`), `neq(f,g)`, `≷` (`\gtrless<tab>`), `≶` (`\lessgtr<tab>`)
* `>=`, `≥` (`\ge<tab>`)
* `>`

To combine predicates, `&` and `|` can be used.

To negate a predicate, `!` is used.

"""
type Pred <: Predicate
    f::Function
    op
    val
end

## meta these
<(f::Function, x::Real) = Pred(f, < , x)
<(f::Function, g::Function) = Pred((x,y) -> f(x,y) - g(x,y), < , 0)

<=(f::Function, x::Real) = Pred(f, <= , x)
<=(f::Function, g::Function) = Pred((x,y) -> f(x,y) - g(x,y), <= , 0)

==(f::Function, x::Real) = Pred(f, == , x)
## ==(f::Function, g::Function) this crosses up Gadfly and others so...
eq(f::Function, g::Function) = Pred((x,y) -> f(x,y) - g(x,y), == , 0)
## unicode variants
⩵(f::Function, x::Real) =  f == x
⩵(f::Function, g::Function) = eq(f,g)

!=(f::Function, x::Real) = Pred(f, != , x)
neq(f::Function, g::Function) = Pred((x,y) -> f(x,y) - g(x,y), != , 0)

≶(x::Real, y::Real) = (x != y)
≶(f::Function, x::Real) = (f != x)
≶(f::Function, g::Function) = neq(f, g)

≷(x::Real, y::Real) = (x != y)
≷(f::Function, x::Real) = (f != x)
≷(f::Function, g::Function) = neq(f, g)

>=(f::Function, x::Real) = Pred(f, >= , x)
>=(f::Function, g::Function) = Pred((x,y) -> f(x,y) - g(x,y), >= , 0)

>(f::Function, x::Real) = Pred(f, > , x)
>(f::Function, g::Function) = Pred((x,y) -> f(x,y) - g(x,y), > , 0)

Base.isless(x::Real, f::Function) = (f >= x)
Base.isless(f::Function, x::Real) = (f < x)

"""

Predicates can be joined together with either `&` or `|`. Individual
predicates can be negated with `!`. The parsing rules require the
individual predicates to be enclosed with parentheses, as in `(f==0) | (g==0)`.

"""
type Preds <: Predicate
    ps
    ops
end

## Some algebra for Pred and Preds
(&)(r1::Pred, r2::Pred) = Preds([r1,r2], Any[&])
(|)(r1::Pred, r2::Pred) = Preds([r1,r2], Any[|])

(&)(ps::Preds, r1::Pred) = Preds([ps.ps, r1], [ps.ops, &])
(&)(r1::Pred, ps::Preds) = ps & r1
(|)(ps::Preds, r1::Pred) = Preds(vcat(ps.ps, r1), vcat(ps.ops, |))
(|)(r1::Pred, ps::Preds) = ps | r1


text:
import Base.size, Base.endof, Base.ndims, Base.convert
export Constant, AffineExpr, AffineConstant, Variable, SumSquaresExpr
export Value, Numeric, AffineOrConstant
export evaluate
export convert
export sum_squares
export endof, size, ndims

Value = Union{Float64, SparseMatrixCSC{Float64, Int64}}
Numeric = Union{Number, AbstractArray}

ValueOrVoid = Union{Value, Void}

abstract AbstractExpr
abstract AffineOrConstant <: AbstractExpr

function evaluate(x::AbstractExpr)
  if x.size == (1, 1)
    return (x.evaluate())[1]
  elseif x.size[2] == 1
    return vec(full(x.evaluate()))
  else
    return full(x.evaluate())
  end
end

type Constant <: AffineOrConstant
  head::Symbol
  value::ValueOrVoid
  size::Tuple{Int64, Int64}
  uid::UInt64
  evaluate::Function

  function Constant(value::Value)
    this = new(:constant, value, (size(value, 1), size(value, 2)))
    this.uid = object_id(this)
    this.evaluate = ()->this.value
    return this
  end
end

function convert(::Type{Constant}, value::Number)
  return Constant(convert(Float64, value))
end

function convert(::Type{Constant}, value::AbstractArray)
  sparse_value = sparse(value)
  if size(sparse_value) == (1, 1)
    return Constant(convert(Float64, value[1]))
  end
  return Constant(convert(SparseMatrixCSC{Float64, Int64}, sparse(value)))
end

Constant(value) = convert(Constant, value)

type AffineExpr <: AffineOrConstant
  head::Symbol
  value::ValueOrVoid
  children::Tuple
  vars_to_coeffs_map::Dict{UInt64, Constant}
  constant::Constant
  size::Tuple{Int64, Int64}
  uid::UInt64
  evaluate::Function

  function AffineExpr(head::Symbol, children::Tuple, vars_to_coeffs_map::Dict{UInt64, Constant}, constant::Constant, size::Tuple{Int64, Int64})
    this = new(head, nothing, children, vars_to_coeffs_map, constant, size)
    this.uid = object_id(this)
    return this
  end
end

function AffineConstant(value::Value)
  constant = Constant(value)
  this = AffineExpr(:constant, (), Dict{UInt64, Constant}(), constant, constant.size)
  this.value = this.constant.value
  this.evaluate = ()->this.value
  return this
end

function Variable(size::Tuple{Int64, Int64})
  if (size[1] < 1 || size[2] < 1)
    error("invalid variable size")
  end
  vec_sz = size[1]*size[2]
  this = AffineExpr(:variable, (), Dict{UInt64, Constant}(), Constant(spzeros(vec_sz, 1)), size)
  this.vars_to_coeffs_map[this.uid] = Constant(speye(vec_sz))
  this.evaluate = ()->begin
    if this.value == nothing
      error("value of the variable is yet to be calculated")
    elseif this.size == (1, 1)
      return this.value[1]
    else
      return this.value
    end
  end
  return this
end

Variable(size...) = Variable(size)
Variable() = Variable((1, 1))
Variable(size::Integer) = Variable((size, 1))

type SumSquaresExpr <: AbstractExpr
  head::Symbol
  value::ValueOrVoid
  affines::Array{AffineExpr}
  multipliers::Array{Float64}
  scalar::Float64
  uid::UInt64
  evaluate::Function

  function SumSquaresExpr(head::Symbol, affines::Array{AffineExpr}, multipliers::Array{Float64}, scalar::Float64)
    this = new(head, nothing, affines, multipliers, scalar)
    this.uid = object_id(this)
    this.evaluate = ()->begin
      value = 0
      for i in 1:length(this.affines)
        value += multipliers[i] * sum((this.affines[i].evaluate()).^2)
      end
      return value + this.scalar
    end
    return this
  end
end

function sum_squares()
  this = SumSquaresExpr(:sum_squares, AffineExpr[], Float64[], 0.0)
end

function sum_squares(num::Number)
  try
    num = convert(Float64, num)
  catch
    error("Sum Squares expressions can only be constructed from real numbers")
  end
  if num < 0
    error("Sum Squares expressions can only be constructed from nonnegative scalars")
  end
  this = SumSquaresExpr(:sum_squares, AffineExpr[], Float64[], num^2)
end

function sum_squares(affine::AffineExpr)
  affines = AffineExpr[]
  push!(affines, affine)
  this = SumSquaresExpr(:sum_squares, affines, ones(length(affines)), 0.0)
end

function evaluate(x::SumSquaresExpr)
  return (x.evaluate())[1]
end

endof(x::AffineOrConstant) = x.size[1] * x.size[2]
function size(x::AffineOrConstant, dim::Integer)
  if dim < 1
    error("dimension out of range")
  elseif dim > 2
    return 1
  else
    return x.size[dim]
  end
end
ndims(x::AffineOrConstant) = 2


text:
import PowerDynamics: rhs, get_current
using DifferentialEquations: DAEFunction
using NetworkDynamics: network_dynamics, StaticEdgeFunction

#function rhs(pg::PowerGrid)
#    network_dynamics(map(construct_vertex, pg.nodes), map(construct_edge, pg.lines), pg.graph; parallel=true)
#end

get_current1(state, n) = begin
    vertices = map(construct_vertex, state.grid.nodes)
    edges = map(construct_edge, state.grid.lines)
    sef = StaticEdgeFunction(vertices,edges,state.grid.graph, parallel=true)
    (e_s, e_d) = sef(state.vec, Nothing, 0)
    total_current(e_s[n], e_d[n])
end

function rhs_dae(pg)
    rpg_ode = rhs(pg)

    du_temp = zeros(systemsize(pg))

    function dae_form(res, du, u, p, t)
        rpg_ode.f(du_temp, u, p, t)
        @. res = du - du_temp
    end

    return DAEFunction{true, true}(dae_form, syms=rpg_ode.syms)
end

using LinearAlgebra: diag
using OrdinaryDiffEq: ODEFunction

function differential_vars(pg)
    rpg_ode = rhs(pg)
    return diag(Array(rpg_ode.mass_matrix)) #(diag(Array(rpg_ode.mass_matrix)) .== 1)
end

function differential_vars(rpg_ode::ODEFunction)
    return (diag(Array(rpg_ode.mass_matrix)) .== 1)
end

import PowerDynamics: find_operationpoint, State
using NLsolve: nlsolve, converged

function find_operationpoint(pg::PowerGrid; ic_guess = nothing, tol=1E-9)
    if SlackAlgebraic ∉ pg.nodes .|> typeof
        @warn "There is no slack bus in the system to balance powers. Currently not making any checks concerning assumptions of whether its possible to find the fixed point"
    end
    if SwingEq ∈ pg.nodes .|> typeof
        throw(OperationPointError("found SwingEq node but these should be SwingEqLVS (just SwingEq is not yet supported for operation point search)"))
    end

    if ic_guess === nothing
        ic_guess = initial_guess(pg)
    end

    rr = RootRhs(rhs(pg))
    nl_res = nlsolve(rr, ic_guess; xtol = tol)
    if converged(nl_res) == true
        return State(pg, nl_res.zero)
    else
        throw(OperationPointError("Failed to find initial conditions on the constraint manifold!"))
    end
end

using ForwardDiff: jacobian
using SparseArrays: sparse
using LinearAlgebra: eigvals
using NLsolve: OnceDifferentiable

"""
Ist das so richtig beschrieben?
Calculates and returns the eigenvalues λ of the jacobian of the right hand side of
the powergrid at the state s.
The sign of the real part of the eigenvalues will decide whether a fixed point
will be attracting (λ_max < 0) or repelling (λ_max > 0).
Inputs:
       pg: Power grid, a graph containg nodes and lines
       s: State of the powergrid
Outputs:
        λ: The biggest and the smallest eigenvalue
        stable: Determines the stability by checking if λ_max is positive
"""
function check_eigenvalues(pg::PowerGrid, s::State)
    rpg = rhs(pg)
    M = Array(rpg.mass_matrix)
    f!(dx, x) = rpg(dx, x, nothing, 0.)
    j(x) = (dx = similar(x); jacobian(f!, dx, x))
    λ = eigvals(j(s.vec) * pinv(M) * M) .|> real |> extrema
    stable = isapprox(last(λ), 0, atol=1e-8)
    return λ, stable
end

function find_operationpoint_sparse(pg::PowerGrid; ic_guess = nothing, tol=1E-9)
    if SlackAlgebraic ∉ pg.nodes .|> typeof
        @warn "There is no slack bus in the system to balance powers. Currently not making any checks concerning assumptions of whether its possible to find the fixed point"
    end
    if SwingEq ∈ pg.nodes .|> typeof
        throw(OperationPointError("found SwingEq node but these should be SwingEqLVS (just SwingEq is not yet supported for operation point search)"))
    end

    if ic_guess === nothing
        system_size = systemsize(pg)
        ic_guess = ones(system_size)
    end

    # we can specify the Jacobian to be sparse
    rpg = rhs(pg)
    f!(dx, x) =  rpg(dx, x, nothing, 0.)
    f(x) = (dx = similar(x); f!(dx, x); dx) # this is the same as RootRhs
    j!(J, x) = (J .= jacobian(f, x))

    n = systemsize(pg)

    dx0 = similar(ic_guess)
    F0 = similar(ic_guess)
    J0 = sparse(zeros(n, n))
    df = OnceDifferentiable(f!, j!, dx0, F0, J0)
    nl_res = nlsolve(df, ic_guess; xtol = tol)

    if converged(nl_res) == true
        return State(pg, nl_res.zero)
    else
        throw(OperationPointError("Failed to find initial conditions on the constraint manifold!"))
    end
end

using OrdinaryDiffEq: ODEFunction
using NLsolve: nlsolve, converged
using LinearAlgebra: pinv

struct RootRhs_ic
    rhs
    mpm
end
function (rr::RootRhs_ic)(x)
    dx = similar(x)
    rr.rhs(dx, x, nothing, 0.)
    rr.mpm * dx .- dx
end

function RootRhs_ic(of::ODEFunction)
    mm = of.mass_matrix
    @assert mm != nothing
    mpm = pinv(mm) * mm
    RootRhs_ic(of.f, mpm)
end

function find_valid_initial_condition(pg, ic_guess)
    rr = RootRhs_ic(rhs(pg))
    nl_res = nlsolve(rr, ic_guess)
    if converged(nl_res) == true
        return State(pg, nl_res.zero)
    else
        println("Failed to find initial conditions on the constraint manifold!")
        println("Try running nlsolve with other options.")
    end
end

using PowerDynamics: AbstractNode
"""
Makes an type specific initial guess to help the operation point search.
The voltage is of all nodes is fixed to the voltage of the first SlackAlgebraic
in the system. The other symbols are set to zero.
Inputs:
    pg: Power grid, a graph containg nodes and lines
Outputs:
    guess: Type specific initial guess
"""
function initial_guess(pg)
    if SlackAlgebraic ∉ pg.nodes .|> typeof
        @warn "There is no slack bus in the system to balance powers."
    end

    sl = findfirst(SlackAlgebraic  ∈  pg.nodes .|> typeof)
    slack = pg.nodes[sl]

    type_guesses = guess.(pg.nodes, slack.U)
    return vcat(type_guesses...)
end

"""
guess(::Type{SlackAlgebraic}) = [slack.U, 0.]         #[:u_r, :u_i]
guess(::Type{PQAlgebraic}) = [slack.U, 0.]            #[:u_r, :u_i]
guess(::Type{ThirdOrderEq}) = [slack.U, 0., 0., 0.]   #[:u_r, :u_i, :θ, :ω]
"""
function guess(n::AbstractNode, slack_voltage)
    voltage = zeros(dimension(n))
    voltage[1] = real(slack_voltage)  #[:u_r, :u_i]
    voltage[2] = imag(slack_voltage)
    return voltage
end

function guess(n::SlackAlgebraic, slack_voltage)
    @assert n.U == slack_voltage
    return [real(n.U), imag(n.U)]
end

function find_steady_state(pg, op_guess=nothing)
    if op_guess isa Nothing
        op_guess = initial_guess(pg)
    end
    ode = rhs(pg)
    op_prob = ODEProblem(ode, op_guess, (0.0, 1.0))
    sol = solve(
        SteadyStateProblem(op_prob),
        DynamicSS(Rodas5(); abstol = 1e-8, reltol = 1e-6, tspan = Inf),
    )
    return State(pg, sol.u)
end

using OrdinaryDiffEq: ODEFunction
using LightGraphs: AbstractGraph
# stack dynamics with higher level controller
struct ControlledPowerGrid
    graph:: G where G <: AbstractGraph
    nodes
    lines
    controller # oop function
end

# Constructor
function ControlledPowerGrid(control, pg::PowerGrid, args...)
    ControlledPowerGrid(pg.graph, pg.nodes, pg.lines, (u, p, t) -> control(u, p, t, args...))
end

import PowerDynamics: PowerGrid

PowerGrid(pg::PowerGrid) = pg
PowerGrid(pg::ControlledPowerGrid) = PowerGrid(pg.graph, pg.nodes, pg.lines)

import PowerDynamics: State, PowerGridSolution
State(pg::ControlledPowerGrid, vec) = State(PowerGrid(pg), vec)
PowerGridSolution(vec, pg::ControlledPowerGrid) = PowerGridSolution(vec, PowerGrid(pg))

import PowerDynamics: systemsize

systemsize(cpg::ControlledPowerGrid) = systemsize(PowerGrid(cpg.graph, cpg.nodes, cpg.lines))

function rhs(cpg::ControlledPowerGrid)
    open_loop_dyn = PowerGrid(cpg.graph, cpg.nodes, cpg.lines) |> rhs
    function cpg_rhs!(du, u, p, t)
        # Get controller state
        p_cont = cpg.controller(u, p, t)
        # Calculate the derivatives, passing the controller state through
        open_loop_dyn(du, u, p_cont, t)
    end
    ODEFunction{true}(cpg_rhs!, mass_matrix=open_loop_dyn.mass_matrix, syms=open_loop_dyn.syms)
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC