text:
module Randoms

using ScottishTaxBenefitModel
using .Utils:isordered

#=
Some functions to make a huge, 60 digit random number, and extract
groups of digits from it.

One such number is attached to each person and household
when the dataset is created. This seems a simple way to create
repeatable but kinda-sort random numbers.
=#

export
    DLA_TO_PIP,
    UC_TRANSITION,
    R_EMPLOYERS_PENSION,
    mybigrand,
    mybigrandstr,
    randchunk,
    pickfirst,
    strtobi,
    testp

const BR_DIGITS = 60

const DEFAULT_CHUNK_SIZE = 8
const DLA_TO_PIP = 1
const R_EMPLOYERS_PENSION = 7
const UC_TRANSITION = 9

"""
extract the chars from start:len-1 and convert to a probability
"""
function randchunk( b :: String, start::Int, len :: Int = DEFAULT_CHUNK_SIZE ) :: Float64
    div = Float64( 10^(len) )
    s =b[start:start+len-1]
    p = parse( Float64, s )/div

    @assert 0 <= p <= 1 "p out of range $p"
    return p
end

#
# if thresh is 0.2 we want this to return true 20% of the time
#
function testp( b :: String, thresh :: Real, start :: Int ) :: Bool
    p = randchunk( b, start )
    # println( "thresh=|$thresh| p=$p")
    return thresh > p
end

"""
given a set of thresholds e.g 0.1, 0.2 .. 1 pick the first less
than or equal to the random from b,start
"""
function pickfirst( b :: String, start :: Int, thresholds :: Vector ) :: Int
    @assert thresholds[end] ≈ 1.0
    @assert( isordered(thresholds))
    r = randchunk( b, start)
    n = size( thresholds )[1]
    for i in 1:n
        if r <= thresholds[i]
            return i
        end
    end
    # FAIL here on purpose
end

"""
The random number is mangled by spreadsheets even when
it's enclosed in quotes, so one crude fix is to prefix it
with a letter and extract from beyond that.
"""
function strtobi( s :: String ) :: String
   return s[2:end]
end

"""
A random number which always has BR_DIGITS decimal digits.
"""
function mybigrand()::String
   return String(rand( '0':'9', BR_DIGITS  ))
end

"""
one of our a `BR_DIGITS` random numbers, prefixed with a character so spreadsheets
don't mangle it.
"""
function mybigrandstr()::String
   return "X"*mybigrand()
end

# print( "randoms is loaded")

end # module Randoms

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module NeuroAdmin
    using NeuroJulia
    using DataFrames

    function getactivenotebooksessions()
        service = "notebookmanagementservice"
        method = "getdetailedsessionlist"
        requestbody=nothing
        response = NeuroJulia.neurocall(service,method,requestbody)
        if response["Error"] != nothing
            error("Neuroverse Error: " * response["Error"])
        end
        df = DataFrame()
        df[:Host] = response["SessionList"]["Host"]
        df[:Session_User] = response["SessionList"]["Session User"]
        df[:Session_Id] = response["SessionList"]["Session Id"]
        df[:Session_CPU] = response["SessionList"]["Session CPU"]
        df[:Session_Memory] = response["SessionList"]["Session Memory"]
        df[:Session_TmpStorage] = response["SessionList"]["Session Temporary Storage"]
        return df
    end
    function getendpointlog(endpointName,startDate,endDate)
        endPointResultEnvelope = NeuroJulia.neurocall("80","endpointmanagement","GetEndpoints",nothing)
        endpointId = ""
        for endpoint in endPointResultEnvelope["EndPointInfo"]
            if endpoint["Name"] == endpointName
                endpointId = endpoint["EndPointId"]
            end
        end
        request = Dict("EndpointId"=>endpointId,"MessageTypeId"=> "", "EndDate"=>endDate, "StartDate"=>startDate)
        resultEnvelope = NeuroJulia.neurocall("80","endpointmanagement","GetMonitorLogEntries",request)
        return resultEnvelope["DataIngestionLogEntries"]
    end
end

ext: jl

max_stars_repo_name: SnowdenNeuroverse/NeuroJulia

max_issues_repo_name: SnowdenNeuroverse/NeuroJulia

max_forks_repo_name: SnowdenNeuroverse/NeuroJulia


text:
using MakieGallery, AbstractPlotting, GLMakie, Test
using MakieGallery: @block, @cell

push!(MakieGallery.plotting_backends, "GLMakie")
database = MakieGallery.load_database()
tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end
recordings = MakieGallery.record_examples(test_record_path)
@test length(recordings) == length(database)
MakieGallery.run_comparison(test_record_path, tested_diff_path)

empty!(database) # remove other examples
include("glmakie_tests.jl") # include GLMakie specific tests
# THese examples download additional data - don't want to deal with that!
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end
examples = MakieGallery.record_examples(test_record_path)
MakieGallery.run_comparison(test_record_path, tested_diff_path, maxdiff = 0.00001)
# MakieGallery.generate_preview(test_record_path)
# repo = joinpath(homedir(), "ReferenceImages", "gallery")
# cp(test_record_path, repo, force = true)


text:
using Distances: Euclidean, SqEuclidean

function dist_tests(
    Dist,
    ȳ::AbstractVector{<:Real},
    Ȳ::AbstractMatrix{<:Real},
    Ȳsq::AbstractMatrix{<:Real},
    x1::AbstractVector,
    x2::AbstractVector,
    x3::AbstractVector;
    rtol=_rtol,
    atol=_atol,
)
    # Check inputs have valid lengths.
    @assert length(x1) == length(x2)
    @assert length(x1) != length(x3)

    @assert length(ȳ) == length(x1)
    @assert size(Ȳ, 1) == length(x1)
    @assert size(Ȳ, 2) == length(x3)
    @assert size(Ȳsq, 1) == size(Ȳsq, 2)
    @assert size(Ȳsq, 1) == length(x1)

    # Binary == Unary
    @test ew(Dist, x1, x1) ≈ ew(Dist, x1)
    @test pw(Dist, x1, x1) ≈ pw(Dist, x1)

    # Symmetry
    @test ew(Dist, x1, x2) ≈ ew(Dist, x2, x1)
    @test pw(Dist, x1, x3) ≈ collect(pw(Dist, x3, x1)')

    # Diagonal of pw is ew
    @test ew(Dist, x1, x2) ≈ diag(pw(Dist, x1, x2))
    @test ew(Dist, x1) ≈ diag(pw(Dist, x1))

    # adjoints
    adjoint_test((x, x′)->ew(Dist, x, x′), ȳ, x1, x2; rtol=rtol, atol=atol)
    adjoint_test((x, x′)->pw(Dist, x, x′), Ȳ, x1, x3; rtol=rtol, atol=atol)

    adjoint_test(x->ew(Dist, x), ȳ, x1; rtol=rtol, atol=atol)
    adjoint_test(x->pw(Dist, x), Ȳsq, x1; rtol=rtol, atol=atol)
end

@testset "distances" begin
    @testset "$Dist" for Dist in [Euclidean(), SqEuclidean()]
        @testset "Vector" begin
            rng = MersenneTwister(123456)
            P, Q = 5, 3
            dist_tests(
                Dist,
                randn(rng, P),
                randn(rng, P, Q),
                randn(rng, P, P),
                collect(range(-2.0, 2.0; length=P)),
                randn(rng, P),
                randn(rng, Q),
            )
        end
        @testset "ColVecs" begin
            rng = MersenneTwister(123456)
            P, Q, D = 5, 3, 2
            dist_tests(
                Dist,
                randn(rng, P),
                randn(rng, P, Q),
                randn(rng, P, P),
                ColVecs(randn(rng, D, P)),
                ColVecs(randn(rng, D, P)),
                ColVecs(randn(rng, D, Q)),
            )
        end
    end

    @testset "Euclidean with repeated inputs" begin
        N = 5
        x_ = randn(N)
        x = vcat(x_, x_)
        ΔD = randn(2N, 2N)
        adjoint_test(x -> pairwise(Euclidean(), x), ΔD, x; atol=1e-6, rtol=1e-6, print_results=true)
    end
end


text:
using CubeSatSim
using Test

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


text:
# Install dependencies needed to build the documentation.
Base.ACTIVE_PROJECT[] = nothing
empty!(LOAD_PATH)
push!(LOAD_PATH, @__DIR__, "@stdlib")
empty!(DEPOT_PATH)
pushfirst!(DEPOT_PATH, joinpath(@__DIR__, "deps"))
using Pkg
Pkg.instantiate()

using Documenter

baremodule GenStdLib end

# Documenter Setup.

symlink_q(tgt, link) = isfile(link) || symlink(tgt, link)
cp_q(src, dest) = isfile(dest) || cp(src, dest)

# make links for stdlib package docs, this is needed until #552 in Documenter.jl is finished
const STDLIB_DOCS = []
const STDLIB_DIR = Sys.STDLIB
const EXT_STDLIB_DOCS = ["Pkg"]
cd(joinpath(@__DIR__, "src")) do
    Base.rm("stdlib"; recursive=true, force=true)
    mkdir("stdlib")
    for dir in readdir(STDLIB_DIR)
        sourcefile = joinpath(STDLIB_DIR, dir, "docs", "src")
        if dir in EXT_STDLIB_DOCS
            sourcefile = joinpath(sourcefile, "basedocs.md")
        else
            sourcefile = joinpath(sourcefile, "index.md")
        end
        if isfile(sourcefile)
            targetfile = joinpath("stdlib", dir * ".md")
            push!(STDLIB_DOCS, (stdlib = Symbol(dir), targetfile = targetfile))
            if Sys.iswindows()
                cp_q(sourcefile, targetfile)
            else
                symlink_q(sourcefile, targetfile)
            end
        end
    end
end

# Check if we are building a PDF
const render_pdf = "pdf" in ARGS

# Generate a suitable markdown file from NEWS.md and put it in src
str = read(joinpath(@__DIR__, "..", "NEWS.md"), String)
splitted = split(str, "<!--- generated by NEWS-update.jl: -->")
@assert length(splitted) == 2
replaced_links = replace(splitted[1], r"\[\#([0-9]*?)\]" => s"[#\g<1>](https://github.com/JuliaLang/julia/issues/\g<1>)")
write(
    joinpath(@__DIR__, "src", "NEWS.md"),
    """
    ```@meta
    EditURL = "https://github.com/JuliaLang/julia/blob/master/NEWS.md"
    ```

    """ * replaced_links)

Manual = [
    "manual/getting-started.md",
    "manual/variables.md",
    "manual/integers-and-floating-point-numbers.md",
    "manual/mathematical-operations.md",
    "manual/complex-and-rational-numbers.md",
    "manual/strings.md",
    "manual/functions.md",
    "manual/control-flow.md",
    "manual/variables-and-scoping.md",
    "manual/types.md",
    "manual/methods.md",
    "manual/constructors.md",
    "manual/conversion-and-promotion.md",
    "manual/interfaces.md",
    "manual/modules.md",
    "manual/documentation.md",
    "manual/metaprogramming.md",
    "manual/arrays.md",
    "manual/missing.md",
    "manual/networking-and-streams.md",
    "manual/parallel-computing.md",
    "manual/asynchronous-programming.md",
    "manual/multi-threading.md",
    "manual/distributed-computing.md",
    "manual/running-external-programs.md",
    "manual/calling-c-and-fortran-code.md",
    "manual/handling-operating-system-variation.md",
    "manual/environment-variables.md",
    "manual/embedding.md",
    "manual/code-loading.md",
    "manual/profile.md",
    "manual/stacktraces.md",
    "manual/performance-tips.md",
    "manual/workflow-tips.md",
    "manual/style-guide.md",
    "manual/faq.md",
    "manual/noteworthy-differences.md",
    "manual/unicode-input.md",
    "manual/command-line-options.md",
]

BaseDocs = [
    "base/base.md",
    "base/collections.md",
    "base/math.md",
    "base/numbers.md",
    "base/strings.md",
    "base/arrays.md",
    "base/parallel.md",
    "base/multi-threading.md",
    "base/constants.md",
    "base/file.md",
    "base/io-network.md",
    "base/punctuation.md",
    "base/sort.md",
    "base/iterators.md",
    "base/c.md",
    "base/libc.md",
    "base/stacktraces.md",
    "base/simd-types.md",
]

StdlibDocs = [stdlib.targetfile for stdlib in STDLIB_DOCS]

DevDocs = [
    "devdocs/reflection.md",
    "Documentation of Julia's Internals" => [
        "devdocs/init.md",
        "devdocs/ast.md",
        "devdocs/types.md",
        "devdocs/object.md",
        "devdocs/eval.md",
        "devdocs/callconv.md",
        "devdocs/compiler.md",
        "devdocs/functions.md",
        "devdocs/cartesian.md",
        "devdocs/meta.md",
        "devdocs/subarrays.md",
        "devdocs/isbitsunionarrays.md",
        "devdocs/sysimg.md",
        "devdocs/llvm.md",
        "devdocs/stdio.md",
        "devdocs/boundscheck.md",
        "devdocs/locks.md",
        "devdocs/offset-arrays.md",
        "devdocs/require.md",
        "devdocs/inference.md",
        "devdocs/ssair.md",
        "devdocs/gc-sa.md",
    ],
    "Developing/debugging Julia's C code" => [
        "devdocs/backtraces.md",
        "devdocs/debuggingtips.md",
        "devdocs/valgrind.md",
        "devdocs/sanitizers.md",
    ]
]

if render_pdf
const PAGES = [
    "Manual" => ["index.md", Manual...],
    "Base" => BaseDocs,
    "Standard Library" => StdlibDocs,
    "Developer Documentation" => DevDocs,
    hide("NEWS.md"),
]
else
const PAGES = [
    "Julia Documentation" => "index.md",
    hide("NEWS.md"),
    "Manual" => Manual,
    "Base" => BaseDocs,
    "Standard Library" => StdlibDocs,
    "Developer Documentation" => DevDocs,
]
end

const use_revise = "revise=true" in ARGS
if use_revise
    let revise_env = joinpath(@__DIR__, "deps", "revise")
        Pkg.activate(revise_env)
        Pkg.add("Revise"; preserve=Pkg.PRESERVE_NONE)
        Base.ACTIVE_PROJECT[] = nothing
        pushfirst!(LOAD_PATH, revise_env)
    end
end
function maybe_revise(ex)
    use_revise || return ex
    STDLIB_DIR = Sys.STDLIB
    STDLIBS = filter!(x -> isfile(joinpath(STDLIB_DIR, x, "src", "$(x).jl")), readdir(STDLIB_DIR))
    return quote
        $ex
        using Revise
        const STDLIBS = $STDLIBS
        union!(Revise.stdlib_names, Symbol.(STDLIBS))
        Revise.track(Core.Compiler)
        Revise.track(Base)
        for (id, mod) in Base.loaded_modules
            if id.name in STDLIBS
                Revise.track(mod)
            end
        end
        Revise.revise()
    end
end

for stdlib in STDLIB_DOCS
    @eval using $(stdlib.stdlib)
    # All standard library modules get `using $STDLIB` as their global
    DocMeta.setdocmeta!(
        Base.root_module(Base, stdlib.stdlib),
        :DocTestSetup,
        maybe_revise(:(using $(stdlib.stdlib)));
        recursive=true,
    )
end
# A few standard libraries need more than just the module itself in the DocTestSetup.
# This overwrites the existing ones from above though, hence the warn=false.
DocMeta.setdocmeta!(
    SparseArrays,
    :DocTestSetup,
    maybe_revise(:(using SparseArrays, LinearAlgebra));
    recursive=true, warn=false,
)
DocMeta.setdocmeta!(
    SuiteSparse,
    :DocTestSetup,
    maybe_revise(:(using SparseArrays, LinearAlgebra, SuiteSparse));
    recursive=true, warn=false,
)
DocMeta.setdocmeta!(
    UUIDs,
    :DocTestSetup,
    maybe_revise(:(using UUIDs, Random));
    recursive=true, warn=false,
)
DocMeta.setdocmeta!(
    Pkg,
    :DocTestSetup,
    maybe_revise(:(using Pkg, Pkg.Artifacts));
    recursive=true, warn=false,
)
DocMeta.setdocmeta!(
    Base,
    :DocTestSetup,
    maybe_revise(:(;;));
    recursive=true,
)
DocMeta.setdocmeta!(
    Base.BinaryPlatforms,
    :DocTestSetup,
    maybe_revise(:(using Base.BinaryPlatforms));
    recursive=true, warn=false,
)
DocMeta.setdocmeta!(
    Pkg.LazilyInitializedFields,
    :DocTestSetup,
    maybe_revise(:(using Pkg.LazilyInitializedFields));
    recursive=true, warn=false,
)

let r = r"buildroot=(.+)", i = findfirst(x -> occursin(r, x), ARGS)
    global const buildroot = i === nothing ? (@__DIR__) : first(match(r, ARGS[i]).captures)
end

const format = if render_pdf
    Documenter.LaTeX(
        platform = "texplatform=docker" in ARGS ? "docker" : "native"
    )
else
    Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
        canonical = ("deploy" in ARGS) ? "https://docs.julialang.org/en/v1/" : nothing,
        assets = [
            "assets/julia-manual.css",
            "assets/julia.ico",
        ],
        analytics = "UA-28835595-6",
        collapselevel = 1,
        sidebar_sitename = false,
    )
end

makedocs(
    build     = joinpath(buildroot, "doc", "_build", (render_pdf ? "pdf" : "html"), "en"),
    modules   = [Main, Base, Core, [Base.root_module(Base, stdlib.stdlib) for stdlib in STDLIB_DOCS]...],
    clean     = true,
    doctest   = ("doctest=fix" in ARGS) ? (:fix) : ("doctest=only" in ARGS) ? (:only) : ("doctest=true" in ARGS) ? true : false,
    linkcheck = "linkcheck=true" in ARGS,
    linkcheck_ignore = ["https://bugs.kde.org/show_bug.cgi?id=136779"], # fails to load from nanosoldier?
    strict    = true,
    checkdocs = :none,
    format    = format,
    sitename  = "The Julia Language",
    authors   = "The Julia Project",
    pages     = PAGES,
)

# Define our own DeployConfig
struct BuildBotConfig <: Documenter.DeployConfig end
function Documenter.deploy_folder(::BuildBotConfig; devurl, repo, branch, kwargs...)
    haskey(ENV, "DOCUMENTER_KEY") || return Documenter.DeployDecision(; all_ok=false)
    if Base.GIT_VERSION_INFO.tagged_commit
        # Strip extra pre-release info (1.5.0-rc2.0 -> 1.5.0-rc2)
        ver = VersionNumber(VERSION.major, VERSION.minor, VERSION.patch,
            isempty(VERSION.prerelease) ? () : (VERSION.prerelease[1],))
        subfolder = "v$(ver)"
        return Documenter.DeployDecision(; all_ok=true, repo, branch, subfolder)
    elseif Base.GIT_VERSION_INFO.branch == "master"
        return Documenter.DeployDecision(; all_ok=true, repo, branch, subfolder=devurl)
    end
    return Documenter.DeployDecision(; all_ok=false)
end

const devurl = "v$(VERSION.major).$(VERSION.minor)-dev"

# Hack to make rc docs visible in the version selector
struct Versions versions end
function Documenter.Writers.HTMLWriter.expand_versions(dir::String, v::Versions)
    # Find all available docs
    available_folders = readdir(dir)
    cd(() -> filter!(!islink, available_folders), dir)
    filter!(x -> occursin(Base.VERSION_REGEX, x), available_folders)

    # Look for docs for an "active" release candidate and insert it
    vnums = [VersionNumber(x) for x in available_folders]
    master_version = maximum(vnums)
    filter!(x -> x.major == 1 && x.minor == master_version.minor-1, vnums)
    rc = maximum(vnums)
    if !isempty(rc.prerelease) && occursin(r"^rc", rc.prerelease[1])
        src = "v$(rc)"
        @assert src ∈ available_folders
        push!(v.versions, src => src, pop!(v.versions))
    end

    return Documenter.Writers.HTMLWriter.expand_versions(dir, v.versions)
end

deploydocs(
    repo = "github.com/JuliaLang/docs.julialang.org.git",
    deploy_config = BuildBotConfig(),
    target = joinpath(buildroot, "doc", "_build", "html", "en"),
    dirname = "en",
    devurl = devurl,
    versions = Versions(["v#.#", devurl => devurl]),
)


text:
using LightGraphs
import JSON

model = JSON.parsefile("meta.json")
code = "# This file has been generated automatically by 'generate.jl'\n\n"
code = code * "__precompile__()\n\n"
code = code * "module Traits\n\n"
code = code * "export VertexTrait, EdgeTrait, refines\n\n"
# code = code * "import ForSyDeIO.Models as Models\n\n"
# open("package_template.py", 'r') do template_file
#     template = jinja2.Template(template_file.read())
# end

# build the trait vertex graph
vertex_traits_idx = [m for m in keys(model["vertexTraits"])]
vertexTraitGraph = SimpleDiGraph(length(model["vertexTraits"]))
# for (vidx, (vname, _)) in enumerate(model["vertexTraits"])
#     vertexTraitGraph.add_node(vname)
# end
for (vidx, (vname, vdata)) in enumerate(model["vertexTraits"])
    if !isnothing(vdata) && "superTraits" in keys(vdata)
        for superTrait in vdata["superTraits"]
            super_idx = findfirst(superTrait .== vertex_traits_idx)
            add_edge!(vertexTraitGraph, vidx, super_idx)
        end
    end
end

code = code * "@enum VertexTrait begin\n"
for (vidx, trait_name) in enumerate(vertex_traits_idx)
    global code
    code = code * "  $trait_name\n"
end
code = code * "end\n\n"
# code = code * "  if t1 == t2\n"
# code = code * "\n"
# for (vidx, trait_name) in enumerate(vertex_traits_idx)
#     global code
#     instance_name = trait_name * "Instance"
#     code = code * "const $instance_name = $trait_name()\n"
# end

code = code * "function refines(t1::VertexTrait, t2::VertexTrait)\n"
for (vidx, trait_name) in enumerate(vertex_traits_idx)
    global code
    code = code * "  if t1 == $trait_name\n"
    for (oidx, trait_name_other) in enumerate(vertex_traits_idx)
        if has_path(vertexTraitGraph, vidx, oidx)
            global code
            code = code * "    if t2 == $trait_name_other\n"#"refines(t1::$trait_name, t2::$trait_name_other) = true\n"
            code = code * "      return true\n"
            #  elseif t1 == $trait_name && t2 == $trait_name_other\n"
            code = code * "    end\n"
        end
    end
    code = code * "  end\n"
end
code = code * "  return false\n"
code = code * "end\n\n"

# build the trait edge graph
edge_traits_idx = [m for m in keys(model["edgeTraits"])]
edgeTraitGraph = SimpleDiGraph(length(edge_traits_idx))
for (vidx, (vname, vdata)) in enumerate(model["edgeTraits"])
    if !isnothing(vdata) && "superTraits" in keys(vdata)
        for superTrait in vdata["superTraits"]
            super_idx = findfirst(superTrait .== edge_traits_idx)
            add_edge!(edgeTraitGraph, vidx, super_idx)
        end
    end
end

code = code * "@enum EdgeTrait begin\n"
for (vidx, trait_name) in enumerate(edge_traits_idx)
    global code
    code = code * "  $trait_name\n"
end
code = code * "end\n\n"
    # code = code * "  if t1 == t2\n"
    # code = code * "\n"
    # for (vidx, trait_name) in enumerate(vertex_traits_idx)
    #     global code
    #     instance_name = trait_name * "Instance"
    #     code = code * "const $instance_name = $trait_name()\n"
    # end

code = code * "function refines(t1::EdgeTrait, t2::EdgeTrait)\n"
for (vidx, trait_name) in enumerate(edge_traits_idx)
    global code
    code = code * "  if t1 == $trait_name\n"
    for (oidx, trait_name_other) in enumerate(edge_traits_idx)
        if has_path(edgeTraitGraph, vidx, oidx)
            global code
            code = code * "    if t2 == $trait_name_other\n"#"refines(t1::$trait_name, t2::$trait_name_other) = true\n"
            code = code * "      return true\n"
            #  elseif t1 == $trait_name && t2 == $trait_name_other\n"
            code = code * "    end\n"
        end
    end
    code = code * "  end\n"
end
code = code * "  return false\n"
code = code * "end\n\n"

#code = code * "@enum EdgeTrait <: ForSyDeIO.Model.Trait begin\n"
# for (vidx, trait_name) in enumerate(edge_traits_idx)
#     global code
#     code = code * "struct $trait_name <: Models.EdgeTrait end\n"
# end

# code = code * "\n"
# for (vidx, trait_name) in enumerate(edge_traits_idx)
#     global code
#     instance_name = trait_name * "Instance"
#     code = code * "const $instance_name = $trait_name()\n"
# end
# code = code * "end\n\nfunction refines(t1::EdgeTrait, t2::EdgeTrait)\n"
# code = code * "  if t1 == t2\n"
# for (vidx, trait_name) in enumerate(edge_traits_idx)
#     global code
#     code = code * "\n"
#     for (oidx, trait_name_other) in enumerate(edge_traits_idx)
#         if has_path(edgeTraitGraph, vidx, oidx)
#             global code
#             code = code * "refines(t1::$trait_name, t2::$trait_name_other) = true\n"
#             #code = code * "  elseif t1 == $trait_name && t2 == $trait_name_other\n"
#             #code = code * "    return true\n"
#         end
#     end
# end
code = code * "\n"
code = code * "const vertex_trait_lut = Dict(\n"
code = code * Base.join(map((t) -> "  \"$t\" => $t", vertex_traits_idx), ",\n")
code = code * ")\n\n"

code = code * "\n"
code = code * "const edge_trait_lut = Dict(\n"
code = code * Base.join(map((t) -> "  \"$t\" => $t", edge_traits_idx), ",\n")
code = code * ")\n\n"

#code = code * "@enum EdgeTrait <: ForSyDeIO.Model.Trait begin\n"
code = code * "make_trait_vertex(label::AbstractString) = vertex_trait_lut[label]\n"

code = code * "make_trait_edge(label::AbstractString) = edge_trait_lut[label]\n"

# code = code * "function make_trait(label::AbstractString)::M where M <: Models.Trait\n"
# code = code * "  global vertex_trait_lut\n"
# code = code * "  global edge_trait_lut\n"
# code = code * "  if haskey(vertex_trait_lut, label)\n"
# code = code * "    return vertex_trait_lut[label]\n"
# code = code * "  else\n"
# code = code * "    return edge_trait_lut[label]\n"
# code = code * "end\n\n"

code = code * "is_trait(label::AbstractString)::Bool = haskey(vertex_trait_lut, label) || haskey(edge_trait_lut, label)\n"
# code = code * "end\n\n"

code = code * "end # module"

open("src/ForSyDeIO/Models/Traits.jl", "w") do out_file
    write(out_file, code)
end

max_stars_repo_name: forsyde/forsyde-io

max_issues_repo_name: forsyde/forsyde-io

max_forks_repo_name: forsyde/forsyde-io


text:
# # Iterators

# `for` loops "lower" to `while` loops plus calls to the `iterate` function:

# ```
# for i in iter
#     # body
# end
# ```

# is internally converted to:

# ```
# next = iterate(iter)
# while next !== nothing
#     (i, state) = next
#     # body
#     next = iterate(iter, state)
# end
# ```

# The same applies to comprehensions and generators.

# Note `nothing` is a singleton value (the only value of its type `Nothing`)
# used by convention when there is no value to return.

typeof(print("hello"))

A = ['a','b','c'];

iterate(A)
iterate(A, 2)
iterate(A, 3)
iterate(A, 4)

# Iteration is also used by "destructuring" assignment:

x, y = A
x
y

# Yet another user of this "iteration protocol" is argument "splatting":

string(A)
string('a','b','c')
string(A...)

# ## Iteration utilities

# `collect` gives you elements of an iterator as an array.
# Comprehensions work by calling `collect` on a generator.

collect(pairs(A))
summary(collect(pairs(A)))
collect(zip(100:102,A))

# There are several other highly popular iteration utilities in the
# built-in module `Iterators`:
# - `enumerate`
# - `rest`
# - `take`
# - `drop`
# - `product`
# - `flatten`
# - `partition`

#-

# Some iterators are infinite!
# - `countfrom`
# - `repeated`
# - `cycle`

I = zip(Iterators.cycle(0:1), Iterators.flatten([[2,3],[4,5]]))

collect(I)

collect(Iterators.product(I,A))

string(I...)

# ## Generator syntax

# Generators provide convenient syntax for composing map, filter, product, and
# flatten.

collect( (2i,j,k) for i = 1:3, j = 1:3 for k = 1:j if i==k)

change(x) = ( (p,n÷5,d÷10,q÷25) for q = 0:25:x for d = 0:10:x-q for n = 0:5:x-q-d for p = x-q-d-n )

collect(change(200))

# ## Defining iterators
#
# Defining iterable types is not necessary that often, since you can often
# get what you need by composing existing iterators.

struct SimpleRange
    lo::Int
    hi::Int
end

Base.iterate(r::SimpleRange, state = r.lo) = state > r.hi ? nothing : (state, state+1)
Base.length(r::SimpleRange) = r.hi - r.lo + 1

collect(SimpleRange(2,8))

# ## Iterator traits
#
# For many algorithms, it's useful to know certain properties of an iterator
# up front.
# The most useful is whether it has a known length.

Base.IteratorSize([1])
Base.IteratorSize(Iterators.repeated(1))
Base.IteratorSize(eachline(open("/dev/null")))

# ## Exercise
#
# Define an iterator giving the first N fibonacci numbers.

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# The file test.mseed comes from an older IRIS libmseed, found by anowacki
# It has a more complicated structure than the test.mseed file in more recent
# versions of libmseed, which reads with no issues
printstyled("SEED submodule\n", color=:light_green)
using SeisIO.SEED

printstyled("  info dump\n", color=:light_green)
redirect_stdout(out) do
  dataless_support()
  mseed_support()
  seed_support()
  resp_wont_read()
end

printstyled("  internals\n", color=:light_green)

printstyled("    seed_time\n", color=:light_green)
u16 = ones(UInt16, 3)
u16[3] = 0x0000
@test u2d(1.0e-6*SEED.seed_time(u16, 0x00, 0x00, 0x00, 0)) == DateTime("0001-01-01T00:00:00")

u16[1] = 0x0640 # 1600
u16[2] = 0x003c # 60
@test u2d(1.0e-6*SEED.seed_time(u16, 0x00, 0x00, 0x00, 0)) == DateTime("1600-02-29T00:00:00")

u16[1] = 0x076c # 1900
@test u2d(1.0e-6*SEED.seed_time(u16, 0x00, 0x00, 0x00, 0)) == DateTime("1900-03-01T00:00:00")

u16[1] = 0x07d0 # 2000
@test u2d(1.0e-6*SEED.seed_time(u16, 0x00, 0x00, 0x00, 0)) == DateTime("2000-02-29T00:00:00")

# 23, 59, 59
@test u2d(1.0e-6*SEED.seed_time(u16, 0x17, 0x3b, 0x3b, 0)) == DateTime("2000-02-29T23:59:59")
@test u2d(1.0e-6*SEED.seed_time(u16, 0x17, 0x3b, 0x3b, -110000000)) == DateTime("2000-02-29T23:58:09")
@test u2d(1.0e-6*SEED.seed_time(u16, 0x17, 0x3b, 0x3b, -115900000)) == DateTime("2000-02-29T23:58:03.1")
@test u2d(1.0e-6*SEED.seed_time(u16, 0x17, 0x3b, 0x3b, Int64(typemax(Int32)))) == DateTime("2000-03-01T00:35:46.484")


text:
module DesignDriveCoilFilt
using FFTW,PyPlot,Optim
using LinearAlgebra
using PyPlot
using CSV

# using Interpolations

include("Load_LTSpice_Net.jl")
include("RunAnalysis.jl")
include("SPICE2Matrix.jl")
include("ToroidOptimizer.jl")
include("Filter_Designer.jl")
include("ThermalModeling.jl")
include("WriteLTSPICEFile.jl")
include("ImpedanceTransformations.jl")
include("SaveToroidSVG.jl")
export DesignDriveFilter,
       ToroidOptimizer,
       PipeFlow,
       findResPair,
       findEquivLC,
       Z_Cap,
       Z_Ind

# Write your package code here.

end


text:
using ITensors,
      Test

@testset "Fermions" begin

  @testset "parity_sign function" begin

    # Full permutations
    p1 = [1,2,3]
    @test ITensors.parity_sign(p1) == +1
    p2 = [2,1,3]
    @test ITensors.parity_sign(p2) == -1
    p3 = [2,3,1]
    @test ITensors.parity_sign(p3) == +1
    p4 = [3,2,1]
    @test ITensors.parity_sign(p4) == -1

    ## Partial permutations
    p5 = [2,7]
    @test ITensors.parity_sign(p5) == +1
    p6 = [5,3]
    @test ITensors.parity_sign(p6) == -1
    p7 = [1,9,3,10]
    @test ITensors.parity_sign(p7) == -1
    p8 = [1,12,9,3,11]
    @test ITensors.parity_sign(p8) == +1

  end

end


text:
module ProtoBuf

import Base: setproperty!, getproperty, hasproperty, propertynames, show, copy!, deepcopy, hash, isequal, ==

export writeproto, readproto, ProtoMeta, ProtoMetaAttribs, meta
export isfilled, which_oneof
export setproperty!, getproperty, hasproperty, show, copy!, deepcopy, clear, isinitialized
export hash, isequal, ==
export ProtoType, lookup, enumstr
export ProtoServiceException, ProtoRpcChannel, ProtoRpcController, MethodDescriptor, ServiceDescriptor, ProtoService,
       AbstractProtoServiceStub, GenericProtoServiceStub, ProtoServiceStub, ProtoServiceBlockingStub,
       find_method, get_request_type, get_response_type, get_descriptor_for_type, call_method

include("codec.jl")
include("svc.jl")

include("google/google.jl")

include("gen.jl")
include("utils.jl")

# Include Google ProtoBuf well known types (https://developers.google.com/protocol-buffers/docs/reference/google.protobuf).
# These are part of the `google.protobuf` sub-module and are included automatically by the code generator.
# For hand coded modules, include them with: `using ProtoBuf; using ProtoBuf.google.protobuf`.
include("google/wellknown.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    message(c::GitCommit, raw::Bool=false)

Return the commit message describing the changes made in commit `c`. If
`raw` is `false`, return a slightly "cleaned up" message (which has any
leading newlines removed). If `raw` is `true`, the message is not stripped
of any such newlines.
"""
function message(c::GitCommit, raw::Bool=false)
    local msg_ptr::Cstring
    msg_ptr = raw ? ccall((:git_commit_message_raw, :libgit2), Cstring, (Ptr{Void},), c.ptr) :
                    ccall((:git_commit_message, :libgit2), Cstring, (Ptr{Void},), c.ptr)
    if msg_ptr == C_NULL
        return nothing
    end
    return unsafe_string(msg_ptr)
end

"""
    author(c::GitCommit)

Return the `Signature` of the author of the commit `c`. The author is
the person who made changes to the relevant file(s). See also [`committer`](@ref).
"""
function author(c::GitCommit)
    ptr = ccall((:git_commit_author, :libgit2), Ptr{SignatureStruct}, (Ptr{Void},), c.ptr)
    @assert ptr != C_NULL
    return Signature(ptr)
end

"""
    committer(c::GitCommit)

Return the `Signature` of the committer of the commit `c`. The committer is
the person who committed the changes originally authored by the [`author`](@ref), but
need not be the same as the `author`, for example, if the `author` emailed a patch to
a `committer` who committed it.
"""
function committer(c::GitCommit)
    ptr = ccall((:git_commit_committer, :libgit2), Ptr{SignatureStruct}, (Ptr{Void},), c.ptr)
    @assert ptr != C_NULL
    return Signature(ptr)
end

function Base.show(io::IO, c::GitCommit)
    authstr = sprint(show, author(c))
    cmtrstr = sprint(show, committer(c))
    print(io, "Git Commit:\nCommit Author: $authstr\nCommitter: $cmtrstr\nSHA: $(GitHash(c))\nMessage:\n$(message(c))")
end

""" Wrapper around `git_commit_create` """
function commit(repo::GitRepo,
                refname::AbstractString,
                msg::AbstractString,
                author::GitSignature,
                committer::GitSignature,
                tree::GitTree,
                parents::GitCommit...)
    commit_id_ptr = Ref(GitHash())
    nparents = length(parents)
    parentptrs = Ptr{Void}[c.ptr for c in parents]
    @check ccall((:git_commit_create, :libgit2), Cint,
                 (Ptr{GitHash}, Ptr{Void}, Ptr{UInt8},
                  Ptr{SignatureStruct}, Ptr{SignatureStruct},
                  Ptr{UInt8}, Ptr{UInt8}, Ptr{Void},
                  Csize_t, Ptr{Ptr{Void}}),
                 commit_id_ptr, repo.ptr, isempty(refname) ? C_NULL : refname,
                 author.ptr, committer.ptr,
                 C_NULL, msg, tree.ptr,
                 nparents, nparents > 0 ? parentptrs : C_NULL)
    return commit_id_ptr[]
end

"""Commit changes to repository"""
function commit(repo::GitRepo, msg::AbstractString;
                refname::AbstractString=Consts.HEAD_FILE,
                author::Signature = Signature(repo),
                committer::Signature = Signature(repo),
                tree_id::GitHash = GitHash(),
                parent_ids::Vector{GitHash}=GitHash[])
    # Retrieve tree identifier
    if iszero(tree_id)
        tree_id = with(GitIndex, repo) do idx; write_tree!(idx) end
    end

    # Retrieve parents from HEAD
    if isempty(parent_ids)
        try # if throws then HEAD not found -> empty repo
            push!(parent_ids, GitHash(repo, refname))
        end
    end

    # return commit id
    commit_id  = GitHash()

    # get necessary objects
    tree = GitTree(repo, tree_id)
    auth_sig = convert(GitSignature, author)
    comm_sig = convert(GitSignature, committer)
    parents = GitCommit[]
    try
        for id in parent_ids
            push!(parents, GitCommit(repo, id))
        end
        commit_id = commit(repo, refname, msg, auth_sig, comm_sig, tree, parents...)
    finally
        for parent in parents
            close(parent)
        end
        close(tree)
        close(auth_sig)
        close(comm_sig)
    end
    return commit_id
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
println(.√[1,2,3,4])  # vectorized dot operators
const ⩾ = >=
println(1 ⩾ 2)

println(NaN == NaN, NaN ≠ NaN)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using DataDeps

register(DataDep(
  "construction-maintenance",
  "http://www.stat.ufl.edu/~winner/data/const_maint.txt",
  "http://www.stat.ufl.edu/~winner/data/const_maint.dat",
  "31fe1ce069e1b50c0fabb0e19270df58ecfaed0cf37be74b26a42251000c6c84",
  post_fetch_method=(path -> begin
    # First char of each row is whitespace, remove
    write(path, join(strip.(readlines(path)), "\n"))

    UCIData.process_dataset(path,
        target_index=2,
        feature_indices=3:6,
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# ------------------------------------------------------------------
# Includes para os outros arquivos do resto da turma
#include("funcoes_relax.jl")
include("feasible_solution.jl")

# Código original do Raphael
mutable struct node
  level::Int
  model::JuMP.Model
end

## Checks if solution is binary
function isBinary(model::JuMP.Model, binaryIndices::Vector{Int64})
  isZero = model.colVal[binaryIndices] .== 0
  isOne  = model.colVal[binaryIndices] .== 1
  if sum(isZero + isOne) == length(binaryIndices)
    return true
  end
  return false
end

## Checks if model is max: if it is, converts to min
function convertSense!(m::JuMP.Model)
  if m.objSense == :Max
    m.objSense = :Min
    m.obj = -m.obj
  else
    m.objSense = :Max
    m.obj = -m.obj
  end
end

## Strong branching -- if amount of branches == Inf --> full strong branching
function strong(currentNode::node, binaryIndices::Vector{Int64}, amountOfBranches::Int64,
                pscMatrix::Matrix{Float64})

  isZero = currentNode.model.colVal[binaryIndices] .== 0
  isOne  = currentNode.model.colVal[binaryIndices] .== 1
  indFracIndices = find(isZero .+ isOne .== 0)
  fracIndices = binaryIndices[indFracIndices]

  n = min(amountOfBranches, length(fracIndices))

  indCandidatesToBranch = randperm(length(fracIndices))[1:n]
  candidatesToBranch = fracIndices[indCandidatesToBranch]
  bounds = Array{Float64}(2,n)
  for i = 1:n
    leftModel = deepcopy(currentNode.model)
    leftModel.colUpper[candidatesToBranch[i]] = 0
    leftModel.colLower[candidatesToBranch[i]] = 0

    rightModel = deepcopy(currentNode.model)
    rightModel.colUpper[candidatesToBranch[i]] = 1
    rightModel.colLower[candidatesToBranch[i]] = 1

    solve(leftModel)
    bounds[1,i] = leftModel.objVal
    solve(rightModel)
    bounds[2,i] = rightModel.objVal
  end

  indBestFrac = ceil(Int,indmax(bounds)/2)
  indToBranch = candidatesToBranch[indBestFrac]

  leftModel = deepcopy(currentNode.model)
  leftModel.colUpper[indToBranch] = 0
  leftModel.colLower[indToBranch] = 0

  rightModel = deepcopy(currentNode.model)
  rightModel.colUpper[indToBranch] = 1
  rightModel.colLower[indToBranch] = 1

  leftDelta = abs(currentNode.model.objVal - leftModel.objVal)/currentNode.model.colVal[indToBranch]
  rightDelta = abs(currentNode.model.objVal - rightModel.objVal)/(1-currentNode.model.colVal[indToBranch])

  pscMatrix[indBestFrac, 2]+=1
  pscMatrix[indBestFrac, 1] = (pscMatrix[indBestFrac,1] + ((leftDelta + rightDelta)/2))/pscMatrix[indBestFrac,2]

  leftChild = node(currentNode.level+1, leftModel)
  rightChild = node(currentNode.level+1, rightModel)

  return leftChild, rightChild, pscMatrix
end

## Most fractional branching
function fractional(currentNode::node, binaryIndices::Vector{Int64})

  distance = abs.(currentNode.model.colVal[binaryIndices] - 0.5)
  indFrac = indmin(distance)
  indToSet = binaryIndices[indFrac]

  leftModel = deepcopy(currentNode.model)
  leftModel.colUpper[indToSet] = 0
  leftModel.colLower[indToSet] = 0

  rightModel = deepcopy(currentNode.model)
  rightModel.colUpper[indToSet] = 1
  rightModel.colLower[indToSet] = 1

  leftChild = node(currentNode.level+1, leftModel)
  rightChild = node(currentNode.level+1, rightModel)

  return leftChild, rightChild
end

## Pseudo-cost branching
function pseudocost(currentNode::node, binaryIndices::Vector{Int64}, pscMatrix::Matrix{Float64})

    isZero = currentNode.model.colVal[binaryIndices] .== 0
    isOne  = currentNode.model.colVal[binaryIndices] .== 1
    indFracIndices = find(isZero .+ isOne .== 0)
    fracIndices = binaryIndices[indFracIndices]

    indMaxPSC = indmax(pscMatrix[indFracIndices,1].*(1-currentNode.model.colVal[indFracIndices]))
    indToBranch = fracIndices[indMaxPSC]

    leftModel = deepcopy(currentNode.model)
    leftModel.colUpper[indToBranch] = 0
    leftModel.colLower[indToBranch] = 0

    rightModel = deepcopy(currentNode.model)
    rightModel.colUpper[indToBranch] = 1
    rightModel.colLower[indToBranch] = 1

    leftDelta = abs(currentNode.model.objVal - leftModel.objVal)/currentNode.model.colVal[indToBranch]
    rightDelta = abs(currentNode.model.objVal - rightModel.objVal)/(1-currentNode.model.colVal[indToBranch])

    pscMatrix[indMaxPSC, 2]+=1
    pscMatrix[indMaxPSC, 1] = (pscMatrix[indMaxPSC,1] + ((leftDelta + rightDelta)/2))/pscMatrix[indMaxPSC, 2]

    leftChild = node(currentNode.level+1, leftModel)
    rightChild = node(currentNode.level+1, rightModel)

    return leftChild, rightChild, pscMatrix
end

## Receives node and creates two children by setting a variable to 0 and 1 respectively
function branch(currentNode::node, binaryIndices::Vector{Int64}, method::Symbol;
                amountOfBranches::Int64 = 10^10, pscMatrix::Matrix{Float64} = ones(1,1))

  if method == :fractional
    leftChild, rightChild = fractional(currentNode, binaryIndices)
  elseif method == :strong
    leftChild, rightChild, pscMatrix = strong(currentNode, binaryIndices, amountOfBranches, pscMatrix)
  elseif method == :pseudocost
    leftChild, rightChild, pscMatrix = pseudocost(currentNode, binaryIndices, pscMatrix)
  elseif method == :hybrid
    if currentNode.level <= 3
      leftChild, rightChild, pscMatrix = strong(currentNode, binaryIndices, 10^10, pscMatrix)
    else
      leftChild, rightChild, pscMatrix = pseudocost(currentNode, binaryIndices, pscMatrix)
    end
  else
    println("Error on branching method defition")
  end

  return leftChild, rightChild, pscMatrix
end

function obtainBoundList(nodeList)
  boundList = Array{Float64}(length(nodeList))
  for i = 1 : length(nodeList)
    boundList[i] = nodeList[i].model.objVal
  end
  return boundList
end

## Receives a mixed binary linear JuMP model
function solveMIP(m::JuMP.Model; branchMethod = :strong, traverseMethod = :mixed, boolgrasptsp::Bool = false)

  tic()

  # Check if model is max; if it is, converts to min
  flagConverted = 0
  if m.objSense == :Max
    convertSense!(m)
    flagConverted = 1
  end

  # Best bounds: start as Inf
  bestBound = -1e200
  bestVal = 1e200

  # Create vector of indices of the binary variables
  binaryIndices = find(m.colCat .== :Bin)
  binarySolutions = 0

  flagOpt = 0 # flag that indicates if a viable solution has been found

  # Solve using initial heuristic (GRASP)
  if boolgrasptsp
    m_aux = deepcopy(m)
    status = feasible_solution.grasp(m_aux,true,200,0.03,10000000)
    if status == :SubOptimal && m_aux.objVal < bestVal
      bestVal = m_aux.objVal
      m.colVal = m_aux.colVal
      binarySolutions += 1
      flagOpt+=1
    end
  end

  # Solve linear relaxation
  m.colCat[:] = :Cont
  status = solve(m)
  nodes = Vector{node}(0)
  if status == :Optimal && isBinary(m, binaryIndices)
    # Solution of the relaxed problem is binary: optimal solution
    bestBound = m.objVal
    bestVal = m.objVal
    binarySolutions = 1
  else
    push!(nodes, node(0, m)) # Add root to branch and bound tree
    lastNodeLevel = 0
  end

  iter = 1 # number of visited nodes
  branched = false # flag that indicates if a branch has occurred in this iteration
  tol = 0.01 # tolerance (%)

  if traverseMethod == :depth
    traverse = -1
  else
    traverse = 1
  end

  # Initializing pseudo-cost matrix
  pscMatrix = [ones(length(binaryIndices)) zeros(length(binaryIndices))]

  time0 = time_ns()
  timeLimit = 1800

  while !isempty(nodes) && abs((bestVal - bestBound)/bestVal) > tol && (time_ns()-time0)/1e9 < timeLimit

    # Change traverse method every 10 iterations for better bound discovery
    if traverseMethod == :mixed && iter%10 == 0
      traverse = (-1)*traverse
    end

    branched = false
    # Check node lower bound. If greater than current best UB, prune by limit
    if nodes[1].model.objVal <= bestVal
      status = solve(nodes[1].model)
      if status == :Optimal
        boundList = obtainBoundList(nodes)
        bestBound = minimum(boundList)
        if isBinary(nodes[1].model, binaryIndices)
          # Relaxed solution is binary: optimal solution -- don't branch
          if nodes[1].model.objVal < bestVal
            bestVal = nodes[1].model.objVal
            m.colVal = copy(nodes[1].model.colVal)
            flagOpt = 1
            binarySolutions+=1
          end
        elseif nodes[1].model.objVal <= bestVal
          # Relaxed solution is not binary and should not be pruned by limit -- branch
          leftChild, rightChild, pscMatrix = branch(nodes[1], binaryIndices, branchMethod; pscMatrix = pscMatrix)
          branched = true
        end
      end
    end

    # Check if branched -- if true, add children to the list
    if branched == true
      if traverse == 1 # breadth -- insert children at the end of the list
        push!(nodes, leftChild, rightChild)
        lastNodeLevel = nodes[1].level
        deleteat!(nodes, 1)
      else # depth -- insert children at the beginning of the list
        lastNodeLevel = nodes[1].level
        deleteat!(nodes, 1)
        unshift!(nodes, leftChild, rightChild)
      end
    else
      deleteat!(nodes, 1)
    end
    iter+=1

    if iter == 1 || iter%10 == 0
      println("UB: $bestVal")
      println("LB: $bestBound")
    end
  end

  if flagConverted == 1
    convertSense!(m)
    m.objVal = -bestVal
    m.objBound = -bestBound
  else
    m.objVal = bestVal
    m.objBound = bestBound
  end

  # Return binary variables to original state
  m.colCat[binaryIndices] = :Bin

  # Outputs
  m.ext[:status] = status
  if flagOpt == 1
    m.ext[:status] = :Optimal
  end
  m.ext[:nodes] = iter
  m.ext[:solutions] = binarySolutions
  t = toc()
  m.ext[:time] = t

  # Retorna status :userlimit caso tenha extrapolado o tempo
  if m.ext[:time] >= timeLimit
    m.ext[:status] = :userlimit
  end

  return m.ext[:status]
end


text:
include("./fixed_parameters.jl")
using .fixed_parameters
include("./louisenlund.jl")
using .louisenlund

num_cohorts = 4
regulation = true
move_intensity = 1 # times per day (two-way trip)
#student_break_rules = 0.05 # percent
#teacher_break_rules = 0.01 # percent
free_from_quarantine = 0.05 # percent
time_of_lessons = 0.625 # percent
num_lessons = 10
time_in_house = 0.375 # percent
emergency_regulation = true
intern_random_infection_rate = 0.002
extern_random_infection_rate = 0.01

function simulation(student_break_rules, teacher_break_rules)
    LUND = Louisenlund(num_cohorts, regulation, intern_random_infection_rate, extern_random_infection_rate, teacher_break_rules, student_break_rules)
    LUND.start_infection(LUND)
    data = Array{Float32, 2}(undef, (6, sim_time*day_tick))
    end_infection = false
    while !end_infection
        healthy, dead, infected, immunized = 0, 0, 0, 0

        LUND.ticking(LUND, move_intensity, free_from_quarantine, time_of_lessons, num_lessons, time_in_house)

        for i in LUND.cohorts
            for j in i.crew
                for k in j
                    if k.infected && !k.in_quarantine
                        infected += 1
                    elseif !k.alive
                        dead += 1
                    elseif k.immunized
                        immunized += 1
                    elseif !k.infected
                        healthy += 1
                    end
                end
            end
        end

        data[1,LUND.tick] = Float32(LUND.tick/day_tick)
        data[2,LUND.tick] = Float32(healthy)
        data[3,LUND.tick] = Float32(infected)
        data[4,LUND.tick] = Float32(immunized)
        data[5,LUND.tick] = Float32(dead)
        data[6,LUND.tick] = Float32(LUND.quarantine)

        if LUND.tick%day_tick==0
            println("Day ", data[1,LUND.tick], " Healthy ", data[2,LUND.tick], " Infected ", data[3,LUND.tick], " Immunized ", data[4,LUND.tick], " Dead ", data[5,LUND.tick], " In quarantine ", data[6,LUND.tick])
        end

        if (length(LUND.infected_list)==0 && LUND.quarantine==0) || length(LUND.healthy_list)==0 || LUND.tick÷day_tick==sim_time || data[3,LUND.tick]==0
            end_infection = true
        end
    end

    return data[:,1:LUND.tick]
end

using HDF5
function save_Data(data, path::String)
    h5open(path, "w") do file
        write(file, "data", data)
    end
end

using .Threads
Threads.@threads for break_rules in 0.0:0.2:1.0
    save_Data(simulation(break_rules, break_rules), "./Data/default_c=4_5/break_rules="*string(break_rules)*".h5")
end


text:
using DrWatson
@quickactivate :Continuation # activate package +  using

using Plots
hmap(A::Matrix) = heatmap(1:size(A, 1), 1:size(A, 2), abs.(A); c = :viridis)

N = 128
σ = 0.4
α = 0.8832
γ = 0.01

params = SpiralParams(N, σ, α, γ)
#S = read_spiral_state(params)
S = read_spiral_x(params, sym=true)
U2 = extend_symmetric(S.u, dim=1)

hmap(U2)

S₁ = State(gp=S.gp, u=U2)
S₀_conj = conj(S₁)
g_ft = get_g_ft(S.gp)

V, V1, V2, V3, (g, h, i) = RHS_Spiral(S₁, S₀_conj, g_ft)

hmap(V1)
hmap(V2)
hmap(V3)
hmap(abs.(V))

extend_symmetric([1 2; 3 4; 5 6; 7 8], dim=2)

U2 = vcat(S.u, reverse(S.u, dims=1))

X = vec(S.u)
#X2 = append!(X, reverse(X))

xs = range(-π, π, length=S.gp.N)

hmap(U2)

import PyPlot

PyPlot.imshow(abs.(U2))
PyPlot.savefig(plotsdir("spiral", "mods.png"))


text:
using ParserCombinator

include("ast.jl")
include("ast2cnf.jl")

function create_ast(io::IO)::AstNode
  ws = Drop(Star(Space()))
  mws = Drop(Plus(Space()))
  variable_pattern = p"[a-zA-Z][a-zA-Z0-9]*"

  n_variable = variable_pattern > VariableNode

  n_formula = Delayed()

  left_parenthesis = Alt!(E"(", Error("Expected '('"))
  right_parenthesis = Alt!(E")", Error("Expected ')'"))

  n_compound_begin = Seq!(left_parenthesis, ws)
  n_compound_end = Seq!(ws, right_parenthesis)

  n_and = Seq!(E"and", mws, n_formula, mws, n_formula) > AndNode
  n_or = Seq!(E"or", mws, n_formula, mws, n_formula) > OrNode
  n_not = Seq!(E"not", mws, variable_pattern) > NotNode

  n_compound_inner = Alt!(n_and, n_or, n_not)

  n_compound = Seq!(n_compound_begin, n_compound_inner, n_compound_end)

  e_formula = Error("Expected an identifier or '('")
  n_formula.matcher = Alt!(n_variable, n_compound, e_formula)

  parse_try(io, n_formula)[1]
end

"""
    parse_nnf_formula(io::IO; <keyword arguments>)::Ast2CnfResult

Reads a formula in negation normal form from `io` and returns the same formula in CNF.
The conversion is implemented using Tseytin encoding.

# Arguments
- `only_left_to_right:Bool=false`: only generate left to right implications in Tseytin encoding.

See also [`Ast2CnfResult`](@ref).
"""
function parse_nnf_formula(io::IO; only_left_to_right=false)::Ast2CnfResult
  ast = create_ast(io)
  ast2cnf(ast, only_left_to_right=only_left_to_right)
end


text:
defalg(::typeof(mean), ::RandVar) = SampleMean
defalg(::typeof(mean), ::Prim.PrimRandVar) = Analytic # Fix ambiguities
defalg(::typeof(mean), ::ReplaceRandVar{A,R2} where R2<:RandVar where A<:Prim.PrimRandVar) = Analytic

defalg(dop::DistOp, ::Prim.PrimRandVar) = Analytic
defalg(dop::DistOp, ::ReplaceRandVar{<:Prim.PrimRandVar}) = Analytic
# defalg(dop::DistOp, ::ReplaceRandVar{<:RandVar, RandVar}) = Analytic

for dop in distops_names
  expr = :($dop(x::RandVar) = $dop(x, defalg($dop, x)))
  eval(expr)
end


text:
isdefined(Base, :__precompile__) && __precompile__()
module PowerSimulationsDynamics

####################################### Structs Exports ####################################

# Base Exports
export Simulation
export Simulation!
export execute!
export NetworkSwitch
export ControlReferenceChange
export BranchTrip
export SourceBusVoltageChange
# export BusTrip

export ImplicitModel
export MassMatrixModel

# Export for routines
export small_signal_analysis
export get_state_series
export get_voltage_magnitude_series
export get_voltage_angle_series
export print_device_states
export get_initial_conditions
export get_activepower_series
export get_reactivepower_series

####################################### Package Imports ####################################
import Logging
import InfrastructureSystems
import SciMLBase
import DiffEqBase
import ForwardDiff
import SparseArrays
import LinearAlgebra
import Base.to_index
import NLsolve
import Base.ImmutableDict
import PowerSystems
const PSY = PowerSystems
const IS = InfrastructureSystems
const PSID = PowerSimulationsDynamics

using DocStringExtensions

#@template (FUNCTIONS, METHODS) = """
#                                 $(TYPEDSIGNATURES)
#                                 $(DOCSTRING)
#                                 """

#Structs for General Devices and System
include("base/definitions.jl")
include("base/ports.jl")
include("base/perturbations.jl")
include("base/file_system.jl")
include("base/simulation_model.jl")
include("base/simulation_inputs.jl")
include("base/device_indexing.jl")
include("base/mass_matrix.jl")
include("base/frequency_reference.jl")
include("base/simulation.jl")
include("base/supplemental_accesors.jl")
include("base/simulation_initialization.jl")
include("base/small_signal.jl")

#Common Models
include("models/branch.jl")
include("models/device.jl")
include("models/kirchoff_laws.jl")
include("models/dynline_model.jl")
include("models/ref_transformations.jl")

#Generator Component Models
include("models/generator_models/machine_models.jl")
include("models/generator_models/pss_models.jl")
include("models/generator_models/avr_models.jl")
include("models/generator_models/tg_models.jl")
include("models/generator_models/shaft_models.jl")

#Inverter Component Models
include("models/inverter_models/DCside_models.jl")
include("models/inverter_models/filter_models.jl")
include("models/inverter_models/frequency_estimator_models.jl")
include("models/inverter_models/outer_control_models.jl")
include("models/inverter_models/inner_control_models.jl")
include("models/inverter_models/converter_models.jl")

#Injection Models
include("models/load_models.jl")
include("models/source_models.jl")

#Saturation Models
include("models/saturation_models.jl")

#Initialization Parameters
include("initialization/init_device.jl")

#Initialization Generator
include("initialization/generator_components/init_machine.jl")
include("initialization/generator_components/init_shaft.jl")
include("initialization/generator_components/init_avr.jl")
include("initialization/generator_components/init_tg.jl")
include("initialization/generator_components/init_pss.jl")

#Initialization Inverter
include("initialization/inverter_components/init_filter.jl")
include("initialization/inverter_components/init_DCside.jl")
include("initialization/inverter_components/init_converter.jl")
include("initialization/inverter_components/init_frequency_estimator.jl")
include("initialization/inverter_components/init_inner.jl")
include("initialization/inverter_components/init_outer.jl")

#System Model
include("models/system.jl")

#PostProcessing
include("post_processing/post_proc_common.jl")
include("post_processing/post_proc_generator.jl")
include("post_processing/post_proc_inverter.jl")

#Utils
include("utils/plot_utils.jl")
include("utils/immutable_dicts.jl")
include("utils/print.jl")
include("utils/kwargs_check.jl")
include("utils/logging.jl")

end # module


text:
import MathProgBase

"""
    solvehook(m::Model; [suppress_warnings::Bool = false, U::Number = 10000, diag::Bool = false, active_constr::Bool = false, real_recourse_dim::Int = -1, conic_δ::Bool = false, inactives::Vector = []])
Returns the solution status to solving the flexibility model `m`. This solvehook what `JuMP.solve(::Model)` for flexibility models.
This solves the flexibility index problem using the variables and constraints specified in `m`.

**Arguments**
- `m::Model` The flexibility model.

**Keyword Arguments**
- `suppress_warnings::Bool = false` Indicates if solver warnings should be suppressed.
- `U::Number = 10000` The slack variable upper bound.
- `diag::Bool = false` Indicates if the ellipsoidal uncertainty set is diagonalized (this is only active when an ellipsoidal set is used).
- `active_constr::Bool = false` Indicates if the optional active constraint should be used which enforces how many inequalities are active at the solution, this must be set to `true` for systems without control variables and/or contain state variables.
- `real_recourse_dim::Int = -1` The actual number of recourse variables in case state variables are included as recourse variables. This is mandatory if `active_constr = true` and no state variables are provided.
- `conic_δ::Bool = false` This should be set to `true` if a conic solver is used such as Pajarito.jl.
- `inactives::Vector = []` The indexes of inequality constraints that should be turned off.
"""
function solvehook(m::Model; suppress_warnings::Bool = false, U::Number = 10000, diag::Bool = false, active_constr::Bool = false,
                   real_recourse_dim::Int = -1, conic_δ::Bool = false, inactives::Vector = [])

    # Pull in the constraint data from the JuMP model
    flex_data = getflexibilitydata(m)

    # Read in the data and extract it from the dictionary
    input_dict = MakeInputDict(m)
    fConsts = input_dict["fConsts"]
    fControls = input_dict["fControls"]
    fRandoms = input_dict["fRandoms"]
    fStates = input_dict["fStates"]
    hConsts = input_dict["hConsts"]
    hControls = input_dict["hControls"]
    hRandoms = input_dict["hRandoms"]
    hStates = input_dict["hStates"]
    θ_nom = input_dict["theta_nom"]
    n_z = input_dict["dims"][1]
    n_θ = input_dict["dims"][2]
    n_x = input_dict["dims"][3]
    n_f = input_dict["dims"][4]
    n_h = input_dict["dims"][5]
    if n_x != 0
        state_cols = input_dict["state_cols"]
    end
    inequal_inds = input_dict["inequal_inds"]

    # Check that U is sufficiently large
    if maximum(abs.(fConsts)) >= U
        U_new = 10 ^ ceil(log10(maximum(abs.(fConsts))))
        @warn "The slack upper bound (U) value of $U is too low for this problem. Thus, U is being set to $U_new !"
        U = U_new
    end

    # Inilialize the JuMP model to solve the MICP Problem
    m_solve = Model(solver = m.solver)

    # Initialize the variables
    @variable(m_solve, θ[1:n_θ])
    @variable(m_solve, δ >= 0)
    @variable(m_solve, s[1:n_f] >= 0)
    @variable(m_solve, y[1:n_f], Bin)
    if n_z != 0
        @variable(m_solve, z[1:n_z])
        @variable(m_solve, λ[1:n_f] >= 0)
        if n_h != 0
            @variable(m_solve, μ[1:n_h])
        end
    end
    if n_x != 0
        @variable(m_solve, x[1:n_x])
    end

    # Set objective function
    @objective(m_solve, Min, δ)

    # Turn off constraints that are set to be inactive
    if length(inactives) != 0
        for i = 1:length(inactives)
            JuMP.fix(y[inactives[i]], 0)
        end
    end

    # Add slack variable constraints
    @constraint(m_solve, sUyconst[j = 1:n_f], s[j] - U * (1 - y[j]) <= 0)

    # Add active set contraint
    if real_recourse_dim != -1
        active_constr = true
    end
    if active_constr || n_x != 0 || n_z == 0
        if !active_constr
            @warn "Problem requires use of active set constraint. Setting active_constr = true"
        end
        if real_recourse_dim == -1
            real_recourse_dim = n_z
            if n_h != 0 && n_x == 0
                @warn "real_recourse_dim not specified. Setting real_recourse_dim to $n_z"
            end
        end
        @constraint(m_solve, sum(y[j] for j = 1:n_f) == real_recourse_dim + 1)
    end

    # Add system constraints
    exprs = AddSystemExpressions(m_solve, input_dict)
    if length(inactives) != 0
        actives = deleteat!(Vector(1:n_f), sort(inactives))
        @constraint(m_solve, f[j = actives], exprs[1][j] + s[j] == 0)
    else
        @constraint(m_solve, f[j = 1:n_f], exprs[1][j] + s[j] == 0)
    end
    if n_h != 0
        @constraint(m_solve, h[i = 1:n_h], exprs[2][i] == 0)
    end

    # Add stationary conditions
    if n_z != 0
        @constraint(m_solve, sum(λ[j] for j = 1:n_f) == 1)
        if n_h != 0
            @constraint(m_solve, deriv[i = 1:n_z], sum(λ[j] * fControls[j, i] for j = 1:n_f) + sum(μ[j] * hControls[j, i] for j = 1:n_h) == 0)
        else
            @constraint(m_solve, deriv[i = 1:n_z], sum(λ[j] * fControls[j, i] for j = 1:n_f) == 0)
        end
    end

    # Add binary variable constraints
    if n_z != 0
        @constraint(m_solve, lyconst[j = 1:n_f], λ[j] - y[j] <= 0)
    end

    ## Setup the uncertainty set
    # Add hyperbox contraints if specified
    sqrt_δ = false
    if flex_data.uncertainty_set.name == :Hyperbox
        if length(flex_data.uncertainty_set.pos_dev) != n_θ
            error("The dimensions of random variables and the uncertainty set do not match.")
        end
        @constraint(m_solve, upper[i = 1:n_θ], θ[i] <= θ_nom[i] + flex_data.uncertainty_set.pos_dev[i] * δ)
        @constraint(m_solve, lower[i = 1:n_θ], θ[i] >= θ_nom[i] - flex_data.uncertainty_set.neg_dev[i] * δ)

    # Add p-norm constraints if specified
    elseif flex_data.uncertainty_set.name == :PNorm
        p = flex_data.uncertainty_set.p
        if p == 1
            @variable(m_solve, aux[1:n_θ] >= 0)
            @constraint(m_solve, abs_upper[i = 1:n_θ], aux[i] >= θ[i] - θ_nom[i])
            @constraint(m_solve, abs_lower[i = 1:n_θ], aux[i] >= -(θ[i] - θ_nom[i]))
            @constraint(m_solve, sum(aux[i] for i = 1:n_θ) <= δ)
        elseif p == 2
            @constraint(m_solve, sum((θ - θ_nom)[i]^2 for i = 1:n_θ) <= δ)
            sqrt_δ = true
        else
            @variable(m_solve, aux[1:n_θ] >= 0)
            @variable(m_solve, max_aux >= 0)
            @constraint(m_solve, abs_upper[i = 1:n_θ], aux[i] >= θ[i] - θ_nom[i])
            @constraint(m_solve, abs_lower[i = 1:n_θ], aux[i] >= -(θ[i] - θ_nom[i]))
            @constraint(m_solve, maxs[i = 1:n_θ], max_aux >= aux[i])
            @constraint(m_solve, max_aux <= δ)
        end

    # Add constraints for diagonalized region
    elseif diag && n_θ > 1
        covar = getcovariance(m)
        if n_θ != size(covar)[1]
            error("The dimensions of the covariance matrix and the random variables do not match.")
        end
        D, V = eigen(inv(covar))
        @variable(m_solve, w[1:n_θ])
        @constraint(m_solve, ws[i = 1:n_θ], w[i] == sqrt(D[i]) * sum(V[j, i] * (θ - θ_nom)[j] for j = 1:n_θ))
        if conic_δ
            @constraint(m_solve, sum(w[i] * w[i] for i = 1:n_θ) <= δ^2)
        else
            @constraint(m_solve, sum(w[i] * w[i] for i = 1:n_θ) <= δ)
        end

    # Add constraint for normal undiagonalized region
    elseif n_θ > 1
        covar = getcovariance(m)
        if n_θ != size(covar)[1]
            error("The dimensions of the covariance matrix and the random variables do not match.")
        end
        inv_covar = inv(covar)
        if conic_δ
            @constraint(m_solve, sum(sum((θ - θ_nom)[j] * inv_covar[j, i] * (θ - θ_nom)[i] for j = 1:n_θ) for i = 1:n_θ) <= δ^2)
        else
            @constraint(m_solve, sum(sum((θ - θ_nom)[j] * inv_covar[j, i] * (θ - θ_nom)[i] for j = 1:n_θ) for i = 1:n_θ) <= δ)
        end

    # Add constraint for speical case that there is only 1 random variable
    else
        covar = getcovariance(m)
        if n_θ != size(covar)[1]
            error("The dimensions of the covariance matrix and the random variables do not match.")
        end
        if isa(covar, Array)
            if conic_δ
                @constraint(m_solve, (θ[1] - θ_nom[1]) * (1 / covar[1]) * (θ[1] - θ_nom[1]) <= δ^2)
            else
                @constraint(m_solve, (θ[1] - θ_nom[1]) * (1 / covar[1]) * (θ[1] - θ_nom[1]) <= δ)
            end
        else
            if conic_δ
                @constraint(m_solve, (θ[1] - θ_nom) * (1 / covar) * (θ[1] - θ_nom) <= δ^2)
            else
                @constraint(m_solve, (θ[1] - θ_nom) * (1 / covar) * (θ[1] - θ_nom) <= δ)
            end
        end
    end

    # Enforce the intersection of the traditional set with the positive θ set
    if flex_data.uncertainty_set.only_positive
        @constraint(m_solve, pos_theta[i = 1:n_θ], θ[i] >= 0)
    end

    # Solve the model
    tic_time = @elapsed status = solve(m_solve, suppress_warnings = suppress_warnings, ignore_solve_hook = true)

    if status == :Optimal
        # Parse the optimized random variable values
        theta = getvalue(θ)
        for i = 1:n_θ
            setvalue(Variable(m, flex_data.RVcols[i]), theta[i])
        end

        # Parse the optimized recourse variable values
        if n_z != 0
            recourse = getvalue(z)
            for i = 1:n_z
                setvalue(Variable(m, flex_data.recourse_cols[i]), recourse[i])
            end
        end

        # Parse the optimized recourse variable values
        if n_x != 0
            state = getvalue(x)
            for i = 1:n_x
                setvalue(Variable(m, state_cols[i]), state[i])
            end
        end

        # Parse the active active constraints
        active_inds = findall(abs.(getvalue(y) - ones(n_f)) .<= 1e-4)
        flex_data.active_constraints = inequal_inds[active_inds]

        # Save the flexibility index
        if sqrt_δ
            flex_data.flexibility_index = sqrt(getobjectivevalue(m_solve))
        elseif flex_data.uncertainty_set.name == :Ellipsoid && conic_δ
            flex_data.flexibility_index = getobjectivevalue(m_solve)^2
        else
            flex_data.flexibility_index = getobjectivevalue(m_solve)
        end

        # Save the solution time
        try
            flex_data.solution_time = getsolvetime(m_solve)
        catch
            flex_data.solution_time = tic_time
        end
    end

    return status
end

"""
    MakeInputDict(m::Model)
Returns `input_dict::Dict` which contains the state space representation of the system equations stored in
the flexibility model. The resulting `input_dict` is used by most of the flexibility analysis functions.

**Arguments**
- `m::Model` The flexibility model.
"""
function MakeInputDict(m::Model)
    # Pull in the constraint data from the JuMP model
    flex_data = getflexibilitydata(m)
    constr_coeffs = JuMP.prepConstrMatrix(m)  #A coefficients
    constr_bounds = JuMP.prepConstrBounds(m)  #A lower and A upper

    # Determine which constraints are inequalities and equalities
    equal_inds = findall(constr_bounds[1] .== constr_bounds[2])
    inequal_inds = findall(constr_bounds[1] .!= constr_bounds[2])

    # Parse the problem dimensions
    n_f = length(inequal_inds)
    n_h = length(equal_inds)
    n_θ = flex_data.numRVs
    n_z = flex_data.num_recourse_vars
    n_x = size(constr_coeffs)[2] - n_θ - n_z

    # Check inequality directions and adjust as necessary
    reversed_inds = findall(constr_bounds[2] .== Inf)
    constr_coeffs[reversed_inds, :] *= -1
    constr_bounds[2][reversed_inds] = -constr_bounds[1][reversed_inds]
    constr_bounds[1][reversed_inds] .= -Inf

    # Parse the coefficent matrices needed for the flexibility index problem
    fConsts = -constr_bounds[2][inequal_inds]
    hConsts = -constr_bounds[2][equal_inds]
    fRandoms = constr_coeffs[inequal_inds, flex_data.RVcols]
    hRandoms = constr_coeffs[equal_inds, flex_data.RVcols]
    if n_z != 0
        fControls = constr_coeffs[inequal_inds, flex_data.recourse_cols]
        hControls = constr_coeffs[equal_inds, flex_data.recourse_cols]
    else
        fControls = 0
        hControls = 0
    end
    if n_x != 0
        num_vars = n_x + n_z + n_θ
        var_cols = collect(1:num_vars)
        state_cols = var_cols[BitArray(!any((y-> ==(y, var_cols[i])), flex_data.RVcols) && !any((y-> ==(y, var_cols[i])), flex_data.recourse_cols) for i = 1:num_vars)]
        fStates = constr_coeffs[inequal_inds, state_cols]
        hStates = constr_coeffs[equal_inds, state_cols]
    else
        fStates = 0
        hStates = 0
    end

    # Get the mean
    θ_nom = flex_data.RVmeans

    # Make the input dictionary
    input_dict = Dict("fConsts" => fConsts, "fControls" => fControls, "fRandoms" => fRandoms, "fStates" => fStates,
                      "hConsts" => hConsts, "hControls" => hControls, "hRandoms" => hRandoms, "hStates" => hStates,
                      "theta_nom" => θ_nom, "dims" => [n_z; n_θ; n_x; n_f; n_h], "inequal_inds" => inequal_inds)
    if n_x != 0
        input_dict["state_cols"] = state_cols
    end
    return input_dict
end

"""
    AddSystemExpressions(m::Model, input_dict::Dict, [num_scenarios::Int = 0])
Returns a vector of vectors where the first contains all the inequality expressions corresponding
to the inequalities defined in `input_dict` and the second contains all of the equality expressions
corresponding to the equalities defined in `input_dict`.

**Arguments**
- `m::Model` The flexibility model.
- `input_dict::Dict` Input dictionary as defined by [`MakeInputDict`](@ref).
- `num_scenarios::Int = 0` The number of scenerio subproblems, 0 turns off this feature.
"""
function AddSystemExpressions(m::Model, input_dict::Dict, num_scenarios::Int = 0)
    # Extract information
    fConsts = input_dict["fConsts"]
    fControls = input_dict["fControls"]
    fRandoms = input_dict["fRandoms"]
    fStates = input_dict["fStates"]
    hConsts = input_dict["hConsts"]
    hControls = input_dict["hControls"]
    hRandoms = input_dict["hRandoms"]
    hStates = input_dict["hStates"]
    n_z = input_dict["dims"][1]
    n_θ = input_dict["dims"][2]
    n_x = input_dict["dims"][3]
    n_f = input_dict["dims"][4]
    n_h = input_dict["dims"][5]

    # Make the expressions
    if num_scenarios == 0
        if n_z != 0
            if n_x != 0
                @expression(m, fexpr[j = 1:n_f], fConsts[j] + sum(fControls[j, i] * m[:z][i] for i = 1:n_z) + sum(fRandoms[j, i] * m[:θ][i] for i = 1:n_θ) +
                            sum(fStates[j, i] * m[:x][i] for i = 1:n_x))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h], hConsts[j] + sum(hControls[j, i] * m[:z][i] for i = 1:n_z) + sum(hRandoms[j, i] * m[:θ][i] for i = 1:n_θ) +
                                sum(hStates[j, i] * m[:x][i] for i = 1:n_x))
                end
            else
                @expression(m, fexpr[j = 1:n_f], fConsts[j] + sum(fControls[j, i] * m[:z][i] for i = 1:n_z) + sum(fRandoms[j, i] * m[:θ][i] for i = 1:n_θ))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h], hConsts[j] + sum(hControls[j, i] * m[:z][i] for i = 1:n_z) + sum(hRandoms[j, i] * m[:θ][i] for i = 1:n_θ))
                end
            end
        else
            if n_x != 0
                @expression(m, fexpr[j = 1:n_f], fConsts[j] + sum(fRandoms[j, i] * m[:θ][i] for i = 1:n_θ) + sum(fStates[j, i] * m[:x][i] for i = 1:n_x))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h], hConsts[j] + sum(hRandoms[j, i] * m[:θ][i] for i = 1:n_θ) + sum(hStates[j, i] * m[:x][i] for i = 1:n_x))
                end
            else
                @expression(m, fexpr[j = 1:n_f], fConsts[j] + sum(fRandoms[j, i] * m[:θ][i] for i = 1:n_θ))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h], hConsts[j] + sum(hRandoms[j, i] * m[:θ][i] for i = 1:n_θ))
                end
            end
        end
    else
        if n_z != 0
            if n_x != 0
                @expression(m, fexpr[j = 1:n_f, k = 1:num_scenarios], fConsts[j] + sum(fControls[j, i] * m[:z][i, k] for i = 1:n_z) + sum(fRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ) +
                            sum(fStates[j, i] * m[:x][i, k] for i = 1:n_x))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h, k = 1:num_scenarios], hConsts[j] + sum(hControls[j, i] * m[:z][i, k] for i = 1:n_z) + sum(hRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ) +
                                sum(hStates[j, i] * m[:x][i, k] for i = 1:n_x))
                end
            else
                @expression(m, fexpr[j = 1:n_f, k = 1:num_scenarios], fConsts[j] + sum(fControls[j, i] * m[:z][i, k] for i = 1:n_z) + sum(fRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h, k = 1:num_scenarios], hConsts[j] + sum(hControls[j, i] * m[:z][i, k] for i = 1:n_z) + sum(hRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ))
                end
            end
        else
            if n_x != 0
                @expression(m, fexpr[j = 1:n_f, k = 1:num_scenarios], fConsts[j] + sum(fRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ) + sum(fStates[j, i] * m[:x][i, k] for i = 1:n_x))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h, k = 1:num_scenarios], hConsts[j] + sum(hRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ) + sum(hStates[j, i] * m[:x][i, k] for i = 1:n_x))
                end
            else
                @expression(m, fexpr[j = 1:n_f, k = 1:num_scenarios], fConsts[j] + sum(fRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ))
                if n_h != 0
                    @expression(m, hexpr[j = 1:n_h, k = 1:num_scenarios], hConsts[j] + sum(hRandoms[j, i] * m[:θ][i, k] for i = 1:n_θ))
                end
            end
        end
    end

    # Return expressions
    if n_h != 0
        return [fexpr, hexpr]
    else
        return [fexpr, []]
    end
end


text:
using Caesar
using IncrementalInference
using RoME
using CloudGraphs

@show dbaddress = ARGS[1]
configuration = CloudGraphs.CloudGraphConfiguration(dbaddress, 7474, "", "", dbaddress, 27017, false, "", "");
cloudGraph = connect(configuration, encodePackedType, getpackedtype, decodePackedType);

# register types of interest in CloudGraphs
registerGeneralVariableTypes!(cloudGraph)
IncrementalInference.setCloudDataLayerAPI!()

nprocs() > 1 ? thxl = 2 : nothing

type SLAMWrapper
  fg::FactorGraph
  tree
  lndmidx::Int
end

USEREADY=0
include("SlamConvenienceFunctions.jl")

function tcpStringSLAMServer(;slamdata=Union{},port::Int=60001)
  println("Empty slam object created")
  if slamdata == Union{}
    # this is being replaced by cloudGraph, added here for development period
    fg = emptyFactorGraph()
    fg.cg = cloudGraph
    slamdata = SLAMWrapper(fg, Union{}, 0)
  end

  println("Listenting on $(port)")
  server = listen(port)
  loop = true
  while loop
   sock = accept(server)
     while isopen(sock)
         loop, retstr = parseTCP!(slamdata, readline(sock)[1:(end-1)])
         loop ? (isopen(sock) ? println(sock, retstr) : nothing) : close(sock)
     end
     println("connection lost")
  end
  !loop ? close(server) : nothing
  return slamdata
end

tcpStringSLAMServer()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -8.232182974914569 <= q <= 6.070970235703473)

                     add_NL_constraint(m, :((-0.17705543456567296 + 0.5634748900178774*(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3]))/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + 0.6106029608527352*(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3]))/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + 0.5384388086098779*(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3])))))/(1 + exp(-(-0.17705543456567296 + 0.5634748900178774*(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3]))/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + 0.6106029608527352*(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3]))/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + 0.5384388086098779*(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))))))) + (0.08053767290922664 + 0.911477397877237*(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3]))/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + 0.7134526118957916*(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3]))/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + -0.34097520379571256*(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3])))))/(1 + exp(-(0.08053767290922664 + 0.911477397877237*(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3]))/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + 0.7134526118957916*(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3]))/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + -0.34097520379571256*(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))))))) + (0.5992621446321471 + 0.15304126986204158*(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3]))/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + -0.37134199071135*(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3]))/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + 0.8747239793986465*(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3])))))/(1 + exp(-(0.5992621446321471 + 0.15304126986204158*(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3]))/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + -0.37134199071135*(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3]))/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + 0.8747239793986465*(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "Generators" begin

  @safetestset "New resource" begin
    using SearchLight

    SearchLight.Generator.newresource("Foo")

    @test isfile(joinpath("app", "resources", "foos", "Foos.jl")) == true

    @test isdir(joinpath("db", "migrations")) == true
    # @test occursin("_create_table_foos.jl", Base.Filesystem.readdir(joinpath("db", "migrations"))[1]) == true
    # @test length(SearchLight.Migrations.downed_migrations()) == 1

    @test isfile(joinpath("app", "resources", "foos", "FoosValidator.jl")) == true

    @test isfile(joinpath("test", "foos_test.jl")) == true

  end;

end;


text:
module Test_ode_pendulum_with_Measurement

#=
This test model demonstrates how to simulate a simple pendulum
with DifferentialEquations.jl using different types:

- model1: states and parameters are Float64
- model2: states and some parameters are Measurement{Float64}
          and the remaining parameters are Float64.

DifferentialEquations.jl stores by default t, x, der_x in its result data structure
at communication points (this can be switched off)
In this test model, additionally the complete model (all variables)
are stored at communication points at every time instant as well (= Vector{Model})
=#

using DifferentialEquations, Measurements, StaticArrays, PyPlot

mutable struct Model{R<:Real, F<:Real}
    m::R  # mass
    L::R  # length
    d::R  # damping
    g::F  # gravity (not used as Measurment)
    ϕ::R  # angle
    ω::R  # der(angle)
    α::R  # der(der(angle))
    r::SVector{2,R}   # absolute position of end point

    function Model{R,F}(m::R, L::R, d::R, g::F) where {R<:Real, F<:Real}
        new(m, L, d, g, 0, 0, 0, SVector{2,R}(0,0))
    end
end

Model(;m::R=0, L::R=0, d::R=0, g::F=0) where {R<:Real, F<:Real} = Model{R,F}(m,L,d,g)

#Model(m::R, L::R, d::R, g::F) where {R<:Real, F<:Real} = Model{R,F}(m,L,d,g)
#Model(;m::Real=0, L::Real=0, d::Real=0, g::Real=0) = Model(promote(m,L,d),...)

"simplePendulum: Compute m::Model at actual time instant"
function simplePendulum!(m::Model, t, x)::Nothing
    m.ϕ = x[1]
    m.ω = x[2]
    m.α = (-m.m*m.g*m.L*sin(m.ϕ) - m.d*m.ω) / (m.m*m.L^2)
    m.r = @SVector [m.L*sin(m.ϕ), -m.L*cos(m.ϕ)]
    return nothing
end

"SimplePendulum1: Called by integrator"
function simplePendulum1!(der_x, x, m::Model, t)
    simplePendulum!(m, t, x)
    der_x[1] = m.ω
    der_x[2] = m.α
end

"SimplePendulum2: Called at every communication point"
function simplePendulum2(x, t, integrator)
    m::Model = integrator.p
    simplePendulum!(m, t, x)
    return deepcopy(m)
end

t_start = 0.0
t_inc   = 0.01
t_end   = 7.0
tspan   = (t_start, t_end)
tspan2  = t_start:t_inc:t_end

# Problem 1 with Float64
model1 = Model(m=1.0, L=1.0, d=0.2, g=9.81)
# println("model1 = $model1")
x1₀ = [π/2, 0.0]

saved_values1=SavedValues(typeof(x1₀[1]), Model)
cb1 = SavingCallback(simplePendulum2, saved_values1, saveat=tspan2)

prob1 = ODEProblem(simplePendulum1!, x1₀, tspan, model1)
sol1 = solve(prob1, Tsit5(), reltol = 1e-6, saveat=t_inc, callback=cb1)
#prob1b = deepcopy(prob1)
#sol3 = solve(prob1, Tsit5(), reltol = 1e-6, saveat=t_inc)

figure(1)
clf()
plot(sol1.t, getindex.(sol1.u, 1), label="\$\\varphi\$")
grid(true)
legend()
title("SimplePendulum with Float64")

figure(2)
clf()
t = saved_values1.t
m = saved_values1.saveval
r = getfield.(m, :r)
r1 = getindex.(r,1)
r2 = getindex.(r,2)
plot(t, r1, label="r1")
plot(t, r2, label="r2")
grid(true)
legend()
title("SimplePendulum with Float64")

# Problem 2 with Measurement
model2 = Model(m=1.0 ± 0.1, L=1.00 ± 0.01, d=0.2± 0.05, g=9.81)
# println("model2 = $model2")
# println("typeof(model2.r) = ", typeof(model2.r))
x2₀ = [π/2 ± 0.1, 0.0 ± 0]

saved_values2=SavedValues(typeof(x2₀[1]), Model)
cb2 = SavingCallback(simplePendulum2, saved_values2, saveat=tspan2)

prob2 = ODEProblem(simplePendulum1!, x2₀, tspan, model2)
sol2  = solve(prob2, Tsit5(), reltol = 1e-6, saveat=t_inc, callback=cb2)

ϕ = getindex.(sol2.u, 1)
ϕm = Measurements.value.(ϕ)
ϕu = Measurements.uncertainty.(ϕ)
ϕmax = ϕm + ϕu
ϕmin = ϕm - ϕu

figure(3)
clf()
plot(sol2.t, ϕmax, label="\$\\varphi_{max}\$")
plot(sol2.t, ϕmin, label="\$\\varphi_{min}\$" )
plot(sol2.t, ϕm  , label="\$\\varphi_{mean}\$" )
grid(true)
legend()
title("SimplePendulum with Measurement{Float64}")

end


text:
using Logging
using Qaintessent
using Qaintellect
using Flux
using IterTools: ncycle
using JLD
using Plots
using LaTeXStrings
using LinearAlgebra
push!(LOAD_PATH, "./density_demo/")
using density_demo

global debug = false
global visualize = false

folder_name = "demo"
if !isdir(folder_name)
    mkdir(folder_name)
end
Base.Filesystem.cd(folder_name)

global Nx = 2
global Ny = 2
global isboundary = false

@assert Nx > 0 "Dimension of lattice must be positive and greater than 0"
@assert Ny > 0 "Dimension of lattice must be positive and greater than 0"

"""
Identity map as sparse matrix with real-valued entries.
"""
sparse_identity(n) = sparse(1.0*I, n, n)

function train(θlist, σtherm, ftarget, circ, cgc, steps)
    paras1 = Flux.params(θlist)
    Flux.params!(paras1, circ)
    delete!(paras1, circ.meas[1].operator)
    # there is not actually any input data for training
    data = ncycle([()], steps)

    # define optimizer
    opt = AdaMax(0.005)

    # define evaluation function
    result_loss = []
    result_trdistance = []
    result_fidelity = []

    evalcb() = begin
        push!(result_loss, ftarget(θlist))
        push!(result_trdistance, trdistance(θlist, σtherm, cgc))
        push!(result_fidelity, fidelity(θlist, σtherm, cgc))
        if debug
            @info "ftarget(θlist): $(ftarget(θlist)), reference: $(thermal_logZ(β*H)); trdistance(θlist, σtherm, cgc): $(trdistance(θlist, σtherm, cgc)); fidelity(θlist, σtherm, cgc): $(fidelity(θlist, σtherm, cgc))"
        end
    end

    # perform optimization
    Flux.train!(() -> ftarget(θlist), paras1, data, opt, cb=evalcb)

    result_loss, result_trdistance, result_fidelity
end

function main(H::AbstractMatrix)
    if debug
        @debug "Constructed Hamiltonian is equivalent to sample $(H ≈ H')"
        @assert(H ≈ H')
    end

    if visualize
        # show eigenvalues
        scatter(eigvals(Matrix(H)), xlabel=L"i", ylabel=L"\lambda_i", legend=false)
    end

    reference_energy = thermal_logZ(β*H)

    # Ground truth
    σtherm = thermal_state(β*H)

    # Consistency check
    @assert tr(σtherm)  ≈ 1 "Trace of thermal state should be 1"

    if visualize
        heatmap(log.(max.(σtherm, 1e-4)), title=L"\sigma_{\beta}", yflip=true, aspect_ratio=:equal)
    end

    # parameterized quantum circuit gate chain

    if visualize
        #visualize binary entropy
        plot(0:0.05:1, binary_entropy.(0:0.05:1), xlabel="p", ylabel="S(p)", legend=false)
    end

    losses = []
    trdistances = []
    fidelities = []

    # create circuit_gate

    samples = 50
    steps = 500
    for i in 1:samples
        cgc = create_cgc(L; isboundary=isboundary)

        # parameterized quantum circuit
        circ = Circuit{L}(cgc, [Hmeas])

        # randomize inputs
        θopt = π/2 .* (1 .+ 0.2*randn(L))

        # target function: β tr[ρ H] - S(ρ), with ρ = U ρlatent U†
        ftarget(θlist) = β * apply(latent_density(θlist), circ)[1] - sum(binary_entropy.((1 .+ cos.(θlist)) ./ 2))

        result_loss, result_trdistance, result_fidelity = train(θopt, σtherm, ftarget, circ, cgc, steps)
        push!(losses, result_loss .- reference_energy)
        push!(trdistances, result_trdistance)
        push!(fidelities, result_fidelity)
    end
    @info "Processing data"
    results_string = "beta"*string(β)*".jld"
    save(results_string, "losses", losses,
             "trdistances", trdistances,
             "fidelities", fidelities)
end

Js = 0
Jy = 0
Jx = 1
Jz = 1

# number of lattice sites
#L = size(adj, 1)
L = Nx*Ny

@info "Construct Hamiltonian"
# construct 2d hamiltonian
H = construct_2d_hamiltonian(Nx, Ny; Jx=Jx, Jy=Jy, Jz=Jz)
# construct 1d hamiltonian
#H = construct_1d_hamiltonian(Nx, Ny; Jx=Jx, Jy=Jy, Jz=Jz, Js=Js)

# use Hamiltonian as measurement operator; note that H gets converted to a complex matrix here
global Hmeas = MeasurementOperator(H, Tuple(1:L));

for b in 0:0.1:0.2
    @info "Running fidelity: $b"
    #try
        global β = b
        main(H)
    #catch e
    #    showerror(stdout, e)
    #    @info "Did not run complete $(b)"
    #end
end

max_stars_repo_name: cmendl/density-matrix-bloch-qcircuit

max_issues_repo_name: cmendl/density-matrix-bloch-qcircuit

max_forks_repo_name: cmendl/density-matrix-bloch-qcircuit


text:
"""
    particle 3D
        -fully actuated
"""
struct Particle{T} <: Model{T}
	# dimensions
	nq::Int # generalized coordinates
	nu::Int # controls
	nw::Int # parameters
	nc::Int # contact points

    mass::T
    gravity::T

    friction_body_world::Vector{T}
    friction_joint::Vector{T}
end

# Methods
function lagrangian(model::Particle, q, q̇)
    L = 0.0
    L += 0.5 * model.mass * dot(q̇[1:3], q̇[1:3])
    L -= model.mass * model.gravity * q[3]

    return L
end

function kinematics(model::Particle, q)
    q[1:3]
end

function signed_distance(model::Particle, q)
    q[3:3]
end

function input_jacobian(model::Particle, q)
	[
	 1.0 0.0 0.0;
	 0.0 1.0 0.0;
     0.0 0.0 1.0;
    ]
end

function contact_jacobian(model::Particle, q)
	[
	 1.0 0.0 0.0;
	 0.0 1.0 0.0;
     0.0 0.0 1.0;
    ]
end

# nominal configuration
function nominal_configuration(model::Particle)
	[0.0; 0.0; 1.0]
end

# friction coefficients
friction_coefficients(model::Particle) = model.friction_body_world

# Dimensions
nq = 3 # configuration dimension
nu = 3 # control dimension
nw = 0 # parameters
nc = 1 # number of contact points

# Parameters
gravity = 10.0#9.81
mass = 1.0 # mass

# Contact corner

friction_body_world = [1.0]  # coefficient of friction

# Model
particle = Particle(nq, nu, nw, nc,
		  mass, gravity,
		  friction_body_world, zeros(0))

particle_contact_kinematics = [q -> kinematics(particle, q)]

particle_contact_kinematics_jacobians = [q -> contact_jacobian(particle, q)]

name(::Particle) = :particle
floating_base_dim(::Particle) = 3


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

    melting_curve.jl

    parse data from request for melting curve analysis

    Author: Tom Price
    Date:   July 2019

===============================================================================#

import DataStructures.OrderedDict
import DataFrames.DataFrame
import Memento: debug, info

#===============================================================================
    field names >>
===============================================================================#

const MC_RAW_FIELDS = Dict(
    :temperature            => TEMPERATURE_KEY,
    :fluorescence           => FLUORESCENCE_VALUE_KEY,
    :well                   => WELL_NUM_KEY,
    :channel                => CHANNEL_KEY)

#===============================================================================
    function definitions >>
===============================================================================#

## called by dispatch()
function act(
    ::Type{Val{meltcurve}},
    req         ::Associative;
    out_format  ::OutputFormat = pre_json
)
    debug(logger, "at act(::Type{Val{meltcurve}})")
    #
    ## required fields
    @get_calibration_data_from_req(meltcurve)
    @parse_raw_data_from_req(meltcurve)
    #
    ## keyword arguments
    kwargs = MC_FIELD_DEFS |>
        sift(req_key ∘ field(:key)) |>
        mold() do x
			if isa(req[x.key], Int64) && Int_T != Int64
				x.name => Int_T(req[x.key])
			else
				x.name => req[x.key]
			end
        end
    #
    ## create container for data and parameter values
    interface = McInput(
            parsed_raw_data...,
            calibration_data;
            out_format = out_format,
            kwargs...)
    #
    ## pass data and parameter values to mc_analysis()
    ## which will perform the analysis for the entire dataset
    response =
        mc_analysis(interface)

    return response |> out(out_format)
end ## act(::Type{Val{meltcurve}})

#==============================================================================#

"Extract dimensions of raw melting curve data and format as a DataFrame."
function parse_raw_data(
    ::Union{Type{Val{meltcurve}},Type{Val{thermal_consistency}}},
    raw_dict ::Associative
)
    mc_raw_df = DataFrame()
    foreach(keys(MC_RAW_FIELDS)) do key
        try
            mc_raw_df[key] = raw_dict[MC_RAW_FIELDS[key]]
        catch err
            if isa(err, ArgumentError)
                throw(ArgumentError("the format of the raw data is incorrect:" *
                    "each data field should have the same length"))
            else
                rethrow()
            end ## if
        end ## try
    end ## next key
    mc_raw_df[:well] = mc_raw_df[:well] |> mold(Symbol ∘ Int_T)
    (wells, channels) =
        map([WELL_NUM_KEY, CHANNEL_KEY]) do key
            raw_dict[key] |> unique |> sort
        end
    (num_wells, num_channels) =
        map(length, (wells, channels))
    return (
        mc_raw_df,
        num_wells,
        num_channels,
        wells |> mold(Symbol ∘ Int_T) |> SVector{num_wells,Symbol},
        channels |> SVector{num_channels,Int_T})
end ## parse_raw_data(::Type{Val{meltcurve}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
function integrate(
    f::Function,
    a::Real,
    b::Real;
    method::Symbol = :quadrature,
)
    if method == :quadrature
        res, err = quadgk(f, a, b)
        return res
    elseif method == :simpsons
        return adaptivesimpsons(f, a, b)
    elseif method == :monte_carlo
        return montecarlo(f, a, b)
    else
        throw(
            ArgumentError(
                "method must be :quadrature, :simpsons or :monte_carlo"
            )
        )
    end
end


text:
using Test

@testset "Utilities" begin
    using WeightsAndBiasLogger: string_dict

    test1 = (x=1, y="x")
    @test string_dict("", test1) == Dict("x" => 1, "y" => "x")
    @test string_dict("m1", test1) == Dict("m1/x" => 1, "m1/y" => "x")
    test2 = Dict(:x => 1, :y => "x")
    @test string_dict("", test2; ignores=[:y]) == Dict("x" => 1)
end

@testset "Demo" begin include("demo.jl") end


text:
@testset "compatibility" begin
    p = 2
    k = Knots(rand(10)) + p * Knots(0,1)
    P = FastBSplineSpace(p, k)
    n = dim(P)
    i = n÷2
    t = 0.5

    @test bsplinebasis₊₀(i::Integer,P::AbstractBSplineSpace,t::Real) == bsplinebasis(P,i,t)
    @test bsplinebasis₋₀(i::Integer,P::AbstractBSplineSpace,t::Real) == bsplinebasis(P,i,t)
    @test bsplinebasis(i::Integer,P::AbstractBSplineSpace,t::Real) == bsplinebasis(P,i,t)
    @test bsplinebasis′₊₀(i::Integer,P::AbstractBSplineSpace,t::Real) == bsplinebasis(P,i,t)
    @test bsplinebasis′₋₀(i::Integer,P::AbstractBSplineSpace,t::Real) == bsplinebasis(P,i,t)
    @test bsplinebasis′(i::Integer,P::AbstractBSplineSpace,t::Real) == bsplinebasis(P,i,t)
    @test bsplinesupport(i::Integer,P::AbstractBSplineSpace) == bsplinesupport(P,i)
end