text:
!isempty(ARGS) || error("No config supplied.")
isfile(ARGS[1]) || error("Cannot read '$(ARGS[1])'")
isabspath(ARGS[1]) || error("Please use an absolute path for the config.")
println("Config supplied: '$(ARGS[1])'")
config_file = ARGS[1]

using SVDD, OneClassActiveLearning
using MLKernels
using Memento
using Serialization, Logging, Random, StatsBase

Random.seed!(0)
include(config_file)

function generate_experiments(data_files, data_split_strategy, initial_pool_strategy, solver, num_resamples_initial_pool=5, oracle_param=Dict{Symbol, Any}())
    exp_name = "$data_info"
    exp_dir = joinpath(data_output_root, exp_name)
    println(exp_dir)
    if isdir(exp_dir)
        print("Type 'yes' or 'y' to delete and overwrite experiment $(exp_dir): ")
        argin = readline()
        if argin == "yes" || argin == "y"
            rm(exp_dir, recursive=true)
        else
            error("Overwriting anyways... Just kidding, nothing happened.")
        end
    else
        !isdir(exp_dir) || error("The experiment directory $(exp_dir) already exists.")
    end

    isdir(exp_dir) || mkpath(exp_dir)
    mkpath(joinpath(exp_dir, "log", "experiment"))
    mkpath(joinpath(exp_dir, "log", "worker"))
    mkpath(joinpath(exp_dir, "results"))
    @info "Generating experiment directory with name: $exp_dir and config: $(config_file). ($num_resamples_initial_pool resamples of the initial pool)"

    experiments = []

    for data_file in data_files
        data, labels = load_data(data_file)
        Random.seed!(0)
        for sg in subspace_generators
            for ns in 1:num_resamples_subspaces
                Random.seed!(ns)
                if isa(sg, Function)
                    subspaces = sg(data)
                end
                @info "Learning oracle."
                oracle = OneClassActiveLearning.initialize_oracle(oracle_param[:type], data, labels, oracle_param[:param])
                @info "Learning oracle done."
                for ss in data_split_strategy
                    for ip in initial_pool_strategy
                        for n in 1:num_resamples_initial_pool
                            Random.seed!(n)
                            split_strategy, initial_pools = get_splits_and_init_pools(data, labels, ss[1], ip[1];
                                                                                ss[2]..., ip[2]...)
                            param = Dict(:num_al_iterations => num_al_iterations,
                                         :solver => solver,
                                         :initial_pools => initial_pools,
                                         :adjust_K => true,
                                         :initial_pool_resample_version => n,
                                         :subspace_resample_version => ns,
                                         :classify_precision => classify_precision)

                            for model in models
                                for init_strategy in init_strategies
                                    for query_strategy in query_strategies
                                        out_dir = split(data_file, '/')[end-1]
                                        output_path = joinpath(exp_dir, "results", out_dir)
                                        isdir(output_path) || mkdir(output_path)
                                        experiment = Dict{Symbol, Any}(
                                                :data_file => data_file,
                                                :data_set_name => out_dir,
                                                :split_strategy_name => ss,
                                                :initial_pool_strategy_name => ip,
                                                :model => merge(model, Dict(:init_strategy => init_strategy)),
                                                :query_strategy => Dict(:type => query_strategy[:type],
                                                                        :param => query_strategy[:param]),
                                                :split_strategy => split_strategy,
                                                :oracle => oracle,
                                                :param => param)

                                        if isa(sg, Function)
                                            experiment[:model][:param][:subspaces] = deepcopy(subspaces)
                                            experiment[:query_strategy][:param][:subspaces] = deepcopy(subspaces)
                                        end

                                        exp_hash = hash(sprint(print, experiment))
                                        @show data_file, ss, ip, n, exp_hash
                                        @assert exp_hash == hash(sprint(print, deepcopy(experiment)))
                                        experiment[:hash] = "$exp_hash"

                                        out_name = splitext(splitdir(data_file)[2])[1]
                                        out_name = joinpath(output_path, "$(out_name)_$(query_strategy[:type])_$(model[:type])_$(exp_hash).json")
                                        experiment[:output_file] = out_name
                                        experiment[:log_dir] = joinpath(exp_dir, "log")
                                        push!(experiments, deepcopy(experiment))
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    # save the experimental setup
    cp(@__FILE__, joinpath(exp_dir, splitdir(@__FILE__)[2]), follow_symlinks=true)
    cp(config_file, joinpath(exp_dir, splitdir(config_file)[2]), follow_symlinks=true)

    @info "Created $exp_dir with $(length(experiments)) instances."

    open(joinpath(exp_dir, "experiment_hashes"), "a") do f
        for e in experiments
            write(f, "$(e[:hash])\n")
        end
    end
    # subspaces = map(e -> e[:model][:param][:subspaces], experiments)
    # println("Checksum: final number of unique subspaces: $(length(unique(subspaces)))")
    serialize(open(joinpath(exp_dir, "experiments.jser"), "w"), experiments)
    return experiments
end

all(isdir.(joinpath.(data_input_root, data_dirs))) || error("Not all data dirs are valid.")
data_files = vcat(([joinpath.(data_input_root, x, readdir(joinpath(data_input_root, x))) for x in data_dirs])...)
println("Found $(length(data_files)) data files.")

generate_experiments(data_files, split_strategy, initial_pool_strategy, solver, num_resamples_initial_pool, oracle_param)

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


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

# What is included in the file `FMI3_fmu_to_md.jl` (FMU to model description)?
# - wrappers to call the model description functions from a FMU-instance [exported]

function fmi3GetModelName(fmu::FMU3)
    fmi3GetModelName(fmu.modelDescription)
end

function fmi3GetInstantiationToken(fmu::FMU3)
    fmi3GetInstantiationToken(fmu.modelDescription)
end

function fmi3GetGenerationTool(fmu::FMU3)
    fmi3GetGenerationTool(fmu.modelDescription)
end

function fmi3GetGenerationDateAndTime(fmu::FMU3)
    fmi3GetGenerationDateAndTime(fmu.modelDescription)
end

function fmi3GetVariableNamingConvention(fmu::FMU3)
    fmi3GetVariableNamingConvention(fmu.modelDescription)
end

function fmi3GetNumberOfEventIndicators(fmu::FMU3)
    fmi3GetNumberOfEventIndicators(fmu.modelDescription)
end

function fmi3CanGetSetState(fmu::FMU3)
    fmi3CanGetSetState(fmu.modelDescription)
end

function fmi3CanSerializeFMUState(fmu::FMU3)
    fmi3CanSerializeFMUState(fmu.modelDescription)
end

function fmi3ProvidesDirectionalDerivatives(fmu::FMU3)
    fmi3ProvidesDirectionalDerivatives(fmu.modelDescription)
end

function fmi3ProvidesAdjointDerivatives(fmu::FMU3)
    fmi3ProvidesAdjointDerivatives(fmu.modelDescription)
end

function fmi3IsCoSimulation(fmu::FMU3)
    fmi3IsCoSimulation(fmu.modelDescription)
end

function fmi3IsModelExchange(fmu::FMU3)
    fmi3IsModelExchange(fmu.modelDescription)
end

function fmi3IsScheduledExecution(fmu::FMU3)
    fmi3IsScheduledExecution(fmu.modelDescription)
end

ext: jl


text:
println("Hello, World")
println()
println("Bye, World!")


text:
module StructuralDDP

	using BasisMatrices: Basis, SplineParams, nodes, BasisMatrix, Expanded
	using QuantEcon: qnwnorm, qnwunif, qnwlogn, gridmake, DiscreteDP, VFI
	using Interpolations, Distributions, TreeViews, DataFrames
	using SparseArrays, LinearAlgebra
	using DocStringExtensions
	using Parameters

	import QuantEcon: solve
	import DataFrames: DataFrame
	import Base: Array, isapprox

	export DiscreteDynamicProblem, DDP
	export createmodel
	export solve, transitionmatrix, value, policy
	export drawshocks, simulate, policy, states, value, DataFrame

	include("utils.jl")

	# type constructor
	include("problem/Transition.jl")
	include("problem/DiscreteDynamicProblem.jl")
	include("problem/DiscreteDynamicProblem_interface.jl")

	# Solver
	include("solve/constructors.jl")
	include("solve/main.jl")
	include("solve/rewardmatrix.jl")
	include("solve/transitionmatrix.jl")
	include("solve/singlechoicevar.jl")
	include("solve/twochoicevar.jl")
	include("solve/state_action.jl")
	include("solve/initialendogstatevars.jl")

	# simulator
	include("simulate/constructor.jl")
	include("simulate/shocks.jl")
	include("simulate/simulate.jl")

	# Models for developing: requires NLsolve
	# include("../test/models/CapitalAdjustModel.jl")
	# include("../test/models/CapitalAdjustModel2.jl")
	# createmodel(model::Symbol; kwargs...) = eval(model)(; kwargs...)
	# export CapitalAdjustModel, CapitalAdjustModel2, createmodel

end


text:
# add example for if run case requires a square mesh (e.g. for testing)

function initializeMeshes(p::precipParam)
    sMesh,fMesh = loadMesh(p)

    mesh = precipMesh(sMesh,fMesh)
    if p.MeshFile == "sin_x2.msh"
        sin_x2_transform!(mesh,1.0)
    elseif p.MeshFile == "sin_wall.msh"
        sin_wall_transform!(mesh,1.0)
    end

    return precipMesh(sMesh,fMesh)
  end

function loadMesh(p::precipParam)
    fileDir = string(p.MeshFolder,"/",p.MeshFile)
    N = 50; square = [-2.0,2.0,-1.0,1.0]; order=2

    if p.runID == 1
        sMesh = Mesh(fileDir)
        fMesh = FluidMesh(fileDir)
    elseif p.runID == 2
        sMesh = Mesh(fileDir)
        fMesh = FluidMesh(fileDir)
    elseif p.runID == 3
        error("no mesh defined. see mesh_functions.jl")
    elseif p.runID == 4
        error("no mesh defined. see mesh_functions.jl")
    elseif p.runID == 5
        error("no mesh defined. see mesh_functions.jl")
    elseif p.runID == 6
        sMesh = squareMesh(square,N,order)
        fMesh = squareMeshFluid(square,N)
    elseif p.runID == 7
        sMesh = Mesh(fileDir)
        fMesh = FluidMesh(fileDir)
    else
        error("you shouldn't get this error")
    end

    return sMesh,fMesh
end

function sin_x2_transform!(mesh::precipMesh,ω::T) where T<:Real
    Ns = length(mesh.sMesh.xy)
    Nfu = length(mesh.fMesh.xy)
    Nfp = length(mesh.fMesh.xyp)
    amp = 0.1

    for ti=1:Ns
        xval = mesh.sMesh.xy[ti].x
        yval = mesh.sMesh.xy[ti].y

        # transform
        if xval >= 1.0
            mesh.sMesh.xy[ti].y += sin_x2_transform_func(xval,yval,ω,amp)
        end
    end

    for ti=1:Nfu
        xval = mesh.fMesh.xy[ti].x
        yval = mesh.fMesh.xy[ti].y

        # transform
        if xval >= 1.0
            mesh.fMesh.xy[ti].y += sin_x2_transform_func(xval,yval,ω,amp)
        end
    end

    for ti=1:Nfp
        xval = mesh.fMesh.xyp[ti].x
        yval = mesh.fMesh.xyp[ti].y

        # transform
        if xval >= 1.0
            mesh.fMesh.xyp[ti].y += sin_x2_transform_func(xval,yval,ω,amp)
        end
    end

    nothing
end

function sin_wall_transform!(mesh::precipMesh,ω::T) where T<:Real
    Ns = length(mesh.sMesh.xy)
    Nfu = length(mesh.fMesh.xy)
    Nfp = length(mesh.fMesh.xyp)
    amp = 0.1

    for ti=1:Ns
        xval = mesh.sMesh.xy[ti].x
        yval = mesh.sMesh.xy[ti].y

        # transform
        if xval >= 1.0
            mesh.sMesh.xy[ti].y += sin_wall_transform_func(xval,yval,ω,amp)
        end
    end

    for ti=1:Nfu
        xval = mesh.fMesh.xy[ti].x
        yval = mesh.fMesh.xy[ti].y

        # transform
        if xval >= 1.0
            mesh.fMesh.xy[ti].y += sin_wall_transform_func(xval,yval,ω,amp)
        end
    end

    for ti=1:Nfp
        xval = mesh.fMesh.xyp[ti].x
        yval = mesh.fMesh.xyp[ti].y

        # transform
        if xval >= 1.0
            mesh.fMesh.xyp[ti].y += sin_wall_transform_func(xval,yval,ω,amp)
        end
    end

    nothing
end

sin_x2_transform_func(x,y,ω,amp) = amp*(cos(2*pi*ω*(x - 1.0)) - 1.0)

function sin_wall_transform_func(x,y,ω,amp)
    valTop = amp*(cos(2*pi*ω*(x - 1.0)) - 1.0)
    valBottom = -3/8

    ytop = 3/8
    ybottom = -3/8

    return (y-ybottom)*valTop/(3/4)
end


text:
# This test stuff in the TOML README at https://github.com/toml-lang/toml
@testset "README" begin

@testset "Example" begin
str = """
# This is a TOML document.

title = "TOML Example"

[owner]
name = "Tom Preston-Werner"

[database]
server = "192.168.1.1"
ports = [ 8001, 8001, 8002 ]
connection_max = 5000
enabled = true

[servers]

  # Indentation (tabs and/or spaces) is allowed but not required
  [servers.alpha]
  ip = "10.0.0.1"
  dc = "eqdc10"

  [servers.beta]
  ip = "10.0.0.2"
  dc = "eqdc10"

[clients]
data = [ ["gamma", "delta"], [1, 2] ]

# Line breaks are OK when inside arrays
hosts = [
  "alpha",
  "omega"
]
"""
d = parsestring(str)
@test d["title"] == "TOML Example"
@test d["owner"]["name"] == "Tom Preston-Werner"
@test d["database"] == Dict(
  "server" => "192.168.1.1",
  "ports" => [ 8001, 8001, 8002 ],
  "connection_max" => 5000,
  "enabled" => true,
)
@test d["servers"] == Dict(
  "alpha" => Dict(
  "ip" => "10.0.0.1",
  "dc" => "eqdc10",
  ),
  "beta" => Dict(
    "ip" => "10.0.0.2",
    "dc" => "eqdc10",
  )
)

@test d["clients"]["data"] == [ ["gamma", "delta"], [1, 2] ]
@test d["clients"]["hosts"] == ["alpha", "omega"]

end # testset

@testset "Comment" begin
str = """
# This is a full-line comment
key = "value"  # This is a comment at the end of a line
another = "# This is not a comment"
"""
@test parsestring(str) == Dict("key" => "value", "another" => "# This is not a comment")

end # testset

@testset "Key/Value Pair" begin

str = """
key = # INVALID
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrUnexpectedStartOfValue

str = """
first = "Tom" last = "Preston-Werner" # INVALID
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrExpectedNewLineKeyValue

end

@testset "Keys" begin

str = """
key = "value"
bare_key = "value"
bare-key = "value"
1234 = "value"
"""
@test parsestring(str) == Dict(
  "key" => "value",
  "bare_key" => "value",
  "bare-key" => "value",
  "1234" => "value",
)

str = """
"127.0.0.1" = "value"
"character encoding" = "value"
"ʎǝʞ" = "value"
'key2' = "value"
'quoted "value"' = "value"
"""
@test parsestring(str) == Dict(
  "127.0.0.1" => "value",
  "character encoding" => "value",
  "ʎǝʞ" => "value",
  "key2" => "value",
  "quoted \"value\"" => "value",
)

str = """
= "no key name"  # INVALID
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrEmptyBareKey

str = """
"" = "blank"     # VALID but discouraged
"""
@test parsestring(str) == Dict("" => "blank")

str = """
'' = 'blank'     # VALID but discouraged
"""
@test parsestring(str) == Dict("" => "blank")

str = """
name = "Orange"
physical.color = "orange"
physical.shape = "round"
site."google.com" = true
"""

@test parsestring(str) == Dict(
  "name" => "Orange",
  "physical" => Dict("color" => "orange", "shape" => "round"),
  "site" => Dict("google.com" => true)
)

str = """
# DO NOT DO THIS
name = "Tom"
name = "Pradyun"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrKeyAlreadyHasValue

str = """
# THIS WILL NOT WORK
spelling = "favorite"
"spelling" = "favourite"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrKeyAlreadyHasValue

str = """
3.14159 = "pi"
"""
@test parsestring(str) == Dict("3" => Dict("14159" => "pi"))

str = """
# This makes the key "fruit" into a table.
fruit.apple.smooth = true

# So then you can add to the table "fruit" like so:
fruit.orange = 2
"""
@test parsestring(str) == Dict(
  "fruit" => Dict("orange" => 2,
                  "apple" => Dict("smooth" => true)))

str = """
# THE FOLLOWING IS INVALID

# This defines the value of fruit.apple to be an integer.
fruit.apple = 1

# But then this treats fruit.apple like it's a table.
# You can't turn an integer into a table.
fruit.apple.smooth = true
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrKeyAlreadyHasValue

str = """
# VALID BUT DISCOURAGED

apple.type = "fruit"
orange.type = "fruit"

apple.skin = "thin"
orange.skin = "thick"

apple.color = "red"
orange.color = "orange"
"""
@test parsestring(str) == Dict(
  "apple" => Dict("type" => "fruit", "skin" => "thin", "color" => "red"),
  "orange" => Dict("type" => "fruit", "skin" => "thick", "color" => "orange"),
)

str = """
# RECOMMENDED

apple.type = "fruit"
apple.skin = "thin"
apple.color = "red"

orange.type = "fruit"
orange.skin = "thick"
orange.color = "orange"
"""
@test parsestring(str) == Dict(
  "apple" => Dict("type" => "fruit", "skin" => "thin", "color" => "red"),
  "orange" => Dict("type" => "fruit", "skin" => "thick", "color" => "orange"),
)

end #testset

@testset "String" begin

str = """str = "I'm a string. \\"You can quote me\\". Name\\tJos\\u00E9\\nLocation\\tSF." """
@test parsestring(str) == Dict("str" => "I'm a string. \"You can quote me\". Name\tJos\u00E9\nLocation\tSF.")

str = """str = "I'm a string. \\"You can quote me\\". Name\\tJos\\u00E9\\nLocation\\tSF." """

str = """str1 = \"\"\"
Roses are red
Violets are blue
\"\"\"
"""
@test parsestring(str) == Dict("str1" => """
Roses are red
Violets are blue
""")

str = """
# The following strings are byte-for-byte equivalent:
str1 = "The quick brown fox jumps over the lazy dog."

str2 = \"\"\"
The quick brown \\

  fox jumps over \\
    the lazy dog.\"\"\"

str3 = \"\"\"\\
       The quick brown \\
       fox jumps over \\
       the lazy dog.\\
       \"\"\"
"""
d = @test parsestring(str)["str1"] == parsestring(str)["str2"] == parsestring(str)["str3"]

str = """
str4 = \"\"\"Here are two quotation marks: \"\". Simple enough.\"\"\"
str5 = \"\"\"Here are three quotation marks: \"\"\\\".\"\"\"
str6 = \"\"\"Here are fifteen quotation marks: \"\"\\\"\"\"\\\"\"\"\\\"\"\"\\\"\"\"\\\".\"\"\"
"""
d = parsestring(str)
@test d["str4"] == "Here are two quotation marks: \"\". Simple enough."
@test d["str5"] == "Here are three quotation marks: \"\"\"."
@test d["str6"] == "Here are fifteen quotation marks: \"\"\"\"\"\"\"\"\"\"\"\"\"\"\"."

str = """
# "This," she said, "is just a pointless statement."
str7 = \"\"\"\"This,\" she said, \"is just a pointless statement.\"\"\"\"
"""
@test_broken parsestring(str)

"""
str5 = \"\"\"Here are three quotation marks: \"\"\".\"\"\"  # INVALID
"""

str = raw"""
# What you see is what you get.
winpath  = 'C:\Users\nodejs\templates'
winpath2 = '\\ServerX\admin$\system32\'
quoted   = 'Tom "Dubs" Preston-Werner'
regex    = '<\i\c*\s*>'
"""
d = parsestring(str)
@test d["winpath"] == raw"C:\Users\nodejs\templates"
@test d["winpath2"] == raw"\\ServerX\admin$\system32\\"
@test d["quoted"] == raw"""Tom "Dubs" Preston-Werner"""
@test d["regex"] == raw"<\i\c*\s*>"

str = raw"""
regex2 = '''I [dw]on't need \d{2} apples'''
lines  = '''
The first newline is
trimmed in raw strings.
   All other whitespace
   is preserved.
'''
"""
d = parsestring(str)
@test d["regex2"] == raw"I [dw]on't need \d{2} apples"
@test d["lines"] == raw"""
The first newline is
trimmed in raw strings.
   All other whitespace
   is preserved.
"""

str = """
quot15 = '''Here are fifteen quotation marks: \"\"\"\"\"\"\"\"\"\"\"\"\"\"\"'''

apos15 = "Here are fifteen apostrophes: '''''''''''''''"

# 'That,' she said, 'is still pointless.'
str = ''''That,' she said, 'is still pointless.''''
"""
@test_broken parsestring(str)

str = """
apos15 = '''Here are fifteen apostrophes: ''''''''''''''''''  # INVALID
"""
@test tryparsestring(str) isa ParserError

end # testset

@testset "Integer" begin

str = """
int1 = +99
int2 = 42
int3 = 0
int4 = -17
int5 = -0
int6 = +0
"""
d = parsestring(str)
@test d["int1"] === 99
@test d["int2"] === 42
@test d["int3"] === 0
@test d["int4"] === -17
@test d["int5"] === 0
@test d["int6"] === 0

str = """
int5 = 1_000
int6 = 5_349_221
int7 = 53_49_221  # Indian number system grouping
int8 = 1_2_3_4_5  # VALID but discouraged
"""
d = parsestring(str)
@test d["int5"] == 1_000
@test d["int6"] == 5_349_221
@test d["int7"] == 53_49_221
@test d["int8"] == 1_2_3_4_5

str = """
# hexadecimal with prefix `0x`
hex1 = 0xDEADBEEF
hex2 = 0xdeadbeef
hex3 = 0xdead_beef

# octal with prefix `0o`
oct1 = 0o01234567
oct2 = 0o755 # useful for Unix file permissions

# binary with prefix `0b`
bin1 = 0b11010110
"""
d = parsestring(str)
@test d["hex1"] == 0xDEADBEEF
@test d["hex2"] == 0xdeadbeef
@test d["hex3"] == 0xdead_beef
@test d["oct1"] == 0o01234567
@test d["oct2"] == 0o755
@test d["bin1"] == 0b11010110

#Arbitrary 64-bit signed integers (from −2^63 to 2^63−1) should be accepted and
#handled losslessly. If an integer cannot be represented losslessly, an error
#must be thrown.
str = """
low = -9_223_372_036_854_775_808
high = 9_223_372_036_854_775_807
"""
d = parsestring(str)
@test d["low"] == -9_223_372_036_854_775_808
@test d["high"] == 9_223_372_036_854_775_807

str = """
toolow = -9_223_372_036_854_775_809
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrOverflowError

str = """
toohigh = 9_223_372_036_854_775_808
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrOverflowError

end

@testset "Float" begin

str = """
# fractional
flt1 = +1.0
flt2 = 3.1415
flt3 = -0.01

# exponent
flt4 = 5e+22
flt5 = 1e06
flt6 = -2E-2

# both
flt7 = 6.626e-34
flt8 = 224_617.445_991_228
"""
d = parsestring(str)
@test d["flt1"] == +1.0
@test d["flt2"] == 3.1415
@test d["flt3"] == -0.01
@test d["flt4"] == 5e+22
@test d["flt5"] == 1e+6
@test d["flt6"] == -2E-2
@test d["flt7"] == 6.626e-34
@test d["flt8"] == 224_617.445_991_228

#The decimal point, if used, must be surrounded by at least one digit on each side.
str = """
# INVALID FLOATS
invalid_float_1 = .7
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrLeadingDot

str = """
# INVALID FLOATS
invalid_float_2 = 7.
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrNoTrailingDigitAfterDot

str = """
# INVALID FLOATS
invalid_float_3 = 3.e+20
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrNoTrailingDigitAfterDot

str = """
# infinity
sf1 = inf  # positive infinity
sf2 = +inf # positive infinity
sf3 = -inf # negative infinity

# not a number
sf4 = nan  # actual sNaN/qNaN encoding is implementation-specific
sf5 = +nan # same as `nan`
sf6 = -nan # valid, actual encoding is implementation-specific
"""
d = parsestring(str)
@test d["sf1"] == Inf
@test d["sf2"] == Inf
@test d["sf3"] == -Inf
@test isnan(d["sf4"])
@test isnan(d["sf5"])
@test isnan(d["sf6"])

end

@testset "Boolean" begin
str = """
bool1 = true
bool2 = false
"""
d =  parsestring(str)
@test d["bool1"] === true
@test d["bool2"] === false

@testset "Offset Date-Time" begin

# Added some fractional digits here
str = "odt1 = 1979-05-27T07:32:00.99999Z"
# Truncated milliseconds
@test parsestring(str)["odt1"] == parse(DateTime, "1979-05-27T07:32:00.999Z", dateformat"yyyy-mm-ddTHH:MM:SS.sZ")

# Julia doesn't support offset datetimes
str = "odt2 = 1979-05-27T00:32:00-07:00"
err = tryparsestring(str)
@test_broken err isa Dict
@test err isa Internals.ParserError
@test err.type == Internals.ErrOffsetDateNotSupported

str = "odt3 = 1979-05-27T00:32:00.999999-07:00"
err = tryparsestring(str)
@test_broken err isa Dict
@test err isa Internals.ParserError
@test err.type == Internals.ErrOffsetDateNotSupported

str = "odt4 = 1979-05-27 07:32:00Z"
d = parsestring(str)
@test d["odt4"] == parse(DateTime, "1979-05-27 07:32:00Z", dateformat"yyyy-mm-dd HH:MM:SSZ")

end

@testset "Local Date-Time" begin

str = """
ldt1 = 1979-05-27T07:32:00
ldt2 = 1979-05-27T00:32:00.999999
"""
d = parsestring(str)
@test d["ldt1"] == parse(DateTime, "1979-05-27T07:32:00")
@test d["ldt2"] == parse(DateTime, "1979-05-27T00:32:00.999")

end

@testset "Local Date" begin

str = """
ld1 = 1979-05-27
"""
d = parsestring(str)
@test d["ld1"] == parse(Date, "1979-05-27")

end

@testset "Local Time" begin

str = """
lt1 = 07:32:00
lt2 = 00:32:00.999999
"""
d = parsestring(str)
@test d["lt1"] == Time(07, 32, 00)
@test d["lt2"] == Time(00, 32, 00, 999)

end

end # testset

@testset "Array" begin

str = """
integers = [ 1, 2, 3 ]
colors = [ \"red\", \"yellow\", \"green\" ]
nested_array_of_int = [ [ 1, 2 ], [3, 4, 5] ]
nested_mixed_array = [ [ 1, 2 ], [\"a\", \"b\", \"c\"] ]
string_array = [ \"all\", 'strings', \"\"\"are the same\"\"\", '''type''' ]

# Mixed-type arrays are allowed
numbers = [ 0.1, 0.2, 0.5, 1, 2, 5 ]
contributors = [
  \"Foo Bar <foo@example.com>\",
  { name = \"Baz Qux\", email = \"bazqux@example.com\", url = \"https://example.com/bazqux\" }
]
"""
d =  parsestring(str)
@test d["integers"] == [1,2,3]
@test d["colors"] == ["red", "yellow", "green"]
@test d["nested_array_of_int"] == [[1,2], [3,4,5]]
@test d["nested_mixed_array"] == [[1,2], ["a", "b", "c"]]
@test d["string_array"] == ["all", "strings", "are the same", "type"]
@test all(d["numbers"] .=== Any[0.1, 0.2, 0.5, 1, 2, 5])
@test d["contributors"] == [
    "Foo Bar <foo@example.com>",
    Dict("name" => "Baz Qux",
         "email" => "bazqux@example.com",
         "url" => "https://example.com/bazqux")
]

str = """
integers2 = [
  1, 2, 3
]

integers3 = [
  1,
  2, # this is ok
]
"""
d = parsestring(str)
@test d["integers3"] == [1, 2]
@test d["integers2"] == [1, 2, 3]

end # testset

@testset "Table" begin

s = """
[table]
key1 = "some string"
key2 = "some other string"
"""
d = parsestring(s)
@test d == Dict("table" => Dict("key2" => "some other string", "key1" => "some string"))

str = """
[dog."tater.man"]
type.name = "pug"
"""
d = parsestring(str)
@test d == Dict("dog" => Dict("tater.man" => Dict("type" => Dict("name" => "pug"))))

str = """
[a.b.c]            # this is best practice
[ d.e.f ]          # same as [d.e.f]
[ g .  h  . i ]    # same as [g.h.i]
[ j . "ʞ" . 'l' ]  # same as [j."ʞ".'l']
"""
d = parsestring(str)
@test d == Dict(
  "a" => Dict("b" => Dict("c" => Dict())),
  "d" => Dict("e" => Dict("f" => Dict())),
  "g" => Dict("h" => Dict("i" => Dict())),
  "j" => Dict("ʞ" => Dict("l" => Dict())),
)

str = """
# [x] you
# [x.y] don't
# [x.y.z] need these
[x.y.z.w] # for this to work

[x] # defining a super-table afterward is ok
"""
d =  parsestring(str)
@test d == Dict("x" => Dict("y" => Dict("z" => Dict("w" => Dict()))))

str = """
# [x] you
# [x.y] don't
# [x.y.z] need these
[x.y.z.w] # for this to work
a = 3

[x] # defining a super-table afterward is ok
b = 2
"""
d =  parsestring(str)
@test d == Dict(
  "x" => Dict(
      "b" => 2,
      "y" => Dict("z" => Dict("w" => Dict("a" => 3)))
  )
)

str = """
# DO NOT DO THIS

[fruit]
apple = "red"

[fruit]
orange = "orange"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrDuplicatedKey

str = """
# DO NOT DO THIS EITHER

[fruit]
apple = "red"

[fruit.apple]
texture = "smooth"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrKeyAlreadyHasValue

str = """
# VALID BUT DISCOURAGED
[fruit.apple]
[animal]
[fruit.orange]
"""
d = parsestring(str)
@test d == Dict(
  "fruit" => Dict("apple" => Dict(), "orange" => Dict()),
  "animal" => Dict()
)

str = """
# RECOMMENDED
[fruit.apple]
[fruit.orange]
[animal]
"""
@test d == Dict(
  "fruit" => Dict("apple" => Dict(), "orange" => Dict()),
  "animal" => Dict()
)

str = """
[fruit]
apple.color = "red"
apple.taste.sweet = true

# [fruit.apple]  # INVALID
# [fruit.apple.taste]  # INVALID

[fruit.apple.texture]  # you can add sub-tables
smooth = true
"""
d = parsestring(str)
@test d == Dict(
"fruit" => Dict("apple" => Dict("color" => "red",
                                "taste" => Dict("sweet" => true),
                                "texture" => Dict("smooth" => true)))
)

str = """
[fruit]
apple.color = "red"
apple.taste.sweet = true

[fruit.apple]  # INVALID
"""

str = """
[fruit]
apple.color = "red"
apple.taste.sweet = true

[fruit.apple.taste]  # INVALID
"""

end # testset

@testset "Inline table" begin

str = """
name = { first = "Tom", last = "Preston-Werner" }
point = { x = 1, y = 2 }
animal = { type.name = "pug" }
"""

str2 = """
[name]
first = "Tom"
last = "Preston-Werner"

[point]
x = 1
y = 2

[animal]
type.name = "pug"
"""

@test parsestring(str) == parsestring(str2)

str = """
[product]
type = { name = "Nail" }
type.edible = false  # INVALID
"""
err =  tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrAddKeyToInlineTable

str = """
[product]
type.name = "Nail"
type = { edible = false }  # INVALID
"""
tryparsestring(str) isa ParserError
@test err isa ParserError
@test err.type == Internals.ErrAddKeyToInlineTable

@testset "Array of Tables" begin

str = """
[[products]]
name = "Hammer"
sku = 738594937

[[products]]

[[products]]
name = "Nail"
sku = 284758393

color = "gray"
"""
@test parsestring(str) == Dict(
  "products" => [
    Dict("name" => "Hammer",
         "sku" => 738594937),
    Dict(),
    Dict("name" => "Nail",
         "sku" => 284758393,
         "color" => "gray")
  ]
)

str = """
[[fruit]]
  name = "apple"

  [fruit.physical]  # subtable
    color = "red"
    shape = "round"

  [[fruit.variety]]  # nested array of tables
    name = "red delicious"

  [[fruit.variety]]
    name = "granny smith"

[[fruit]]
  name = "banana"

  [[fruit.variety]]
    name = "plantain"
"""
d = parsestring(str)
@test d == Dict(
  "fruit" => [
      Dict("name" => "apple",
           "physical" => Dict("color" => "red", "shape" => "round"),
           "variety" => [
        Dict("name" => "red delicious"),
        Dict("name" => "granny smith"),
      ]),
      Dict("name" => "banana",
           "variety" => [
             Dict("name" => "plantain"),
          ]
      )
])

str = """
# INVALID TOML DOC
[fruit.physical]  # subtable, but to which parent element should it belong?
  color = "red"
  shape = "round"

[[fruit]]  # parser must throw an error upon discovering that "fruit" is
           # an array rather than a table
  name = "apple"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrArrayTreatedAsDictionary

str = """
# INVALID TOML DOC
fruit = []

[[fruit]] # Not allowed
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrAddArrayToStaticArray

str = """
# INVALID TOML DOC
[[fruit]]
  name = "apple"

  [[fruit.variety]]
    name = "red delicious"

  # INVALID: This table conflicts with the previous array of tables
  [fruit.variety]
    name = "granny smith"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrDuplicatedKey

str = """
# INVALID TOML DOC
[[fruit]]
  name = "apple"

  [fruit.physical]
    color = "red"
    shape = "round"

  # INVALID: This array of tables conflicts with the previous table
  [[fruit.physical]]
    color = "green"
"""
err = tryparsestring(str)
@test err isa ParserError
@test err.type == Internals.ErrArrayTreatedAsDictionary

str = """
points = [ { x = 1, y = 2, z = 3 },
           { x = 7, y = 8, z = 9 },
           { x = 2, y = 4, z = 8 } ]
"""
d = tryparsestring(str)
@test d == Dict(
  "points" => [
    Dict("x" => 1, "y" => 2, "z" => 3) ,
    Dict("x" => 7, "y" => 8, "z" => 9) ,
    Dict("x" => 2, "y" => 4, "z" => 8)
  ]
)

end

end

end # testset README


text:
# Autogenerated wrapper script for IntelOpenMP_jll for i686-linux-gnu
export libiomp

JLLWrappers.@generate_wrapper_header("IntelOpenMP")
JLLWrappers.@declare_library_product(libiomp, "libiomp5.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libiomp,
        "lib/libiomp5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct QGramIterator{S <: Union{AbstractString, AbstractVector}}
	s::S   # Collection
	q::Int # Length of Qgram
	function QGramIterator{S}(s, q) where {S <: Union{AbstractString, AbstractVector}}
		q > 0 || throw(ArgumentError("The qgram length must be higher than zero"))
		new(s, q)
	end
end
function QGramIterator(s::Union{AbstractString, AbstractVector}, q::Integer)
	QGramIterator{typeof(s)}(s, q)
end
Base.length(qgram::QGramIterator) = max(length(qgram.s) - qgram.q + 1, 0)

# q-grams of AbstractString
function Base.iterate(qgram::QGramIterator{<: AbstractString},
	state = (1, nextind(qgram.s, 0, qgram.q)))
	istart, iend = state
	iend > ncodeunits(qgram.s) && return nothing
	element = SubString(qgram.s, istart, iend)
	nextstate = nextind(qgram.s, istart), nextind(qgram.s, iend)
	element, nextstate
end
Base.eltype(qgram::QGramIterator{SubString{S}}) where {S} = SubString{S}
Base.eltype(qgram::QGramIterator{S}) where {S <: AbstractString} = SubString{S}
qgrams(s::AbstractString, q::Integer) = QGramIterator(s, q)

#q-grams of AbstractVector
# Alternatively, I could also use partition in IterTools but it creates a vector for each iteration
# so it does not seem to be worth it.
function Base.iterate(qgram::QGramIterator{<: AbstractVector}, state = firstindex(qgram.s))
	state + qgram.q - 1 > lastindex(qgram.s) && return nothing
	view(qgram.s, state:(state + qgram.q - 1)), state + 1
end
Base.eltype(qgram::QGramIterator{<: AbstractVector}) = typeof(first(qgram))
qgrams(s::AbstractVector, q::Integer) = QGramIterator(s, q)
qgrams(s, q::Integer) = QGramIterator(collect(s), q)

@doc """
Return an iterator corresponding to the the q-gram of an iterator.
When the iterator is a String, qgrams are SubStrings.

### Arguments
* `s` iterator
* `q::Integer`: length of q-gram

## Examples
```julia
for x in qgrams("hello", 2)
	println(x)
end
```
"""
qgrams

# For two iterators s1 and s2, that define a length and eltype method,
# this returns an iterator that,
# for each element in s1 ∪ s2, returns (numbers of times it appears in s1, numbers of times it appears in s2)
function _count(qgrams1, qgrams2)
	K = promote_type(eltype(qgrams1), eltype(qgrams2))
	d = Dict{K, Tuple{Int, Int}}()
	sizehint!(d, length(qgrams1) + length(qgrams2))
	# I use a faster way to change a dictionary key
	# see setindex! in https://github.com/JuliaLang/julia/blob/master/base/dict.jl#L380
	for x1 in qgrams1
		index = Base.ht_keyindex2!(d, x1)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x1
			@inbounds d.vals[index] = (d.vals[index][1] + 1, 0)
		else
			@inbounds Base._setindex!(d, (1, 0), x1, -index)
		end
	end
	for x2 in qgrams2
		index = Base.ht_keyindex2!(d, x2)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x2
			@inbounds d.vals[index] = (d.vals[index][1], d.vals[index][2] + 1)
		else
			@inbounds Base._setindex!(d, (0, 1), x2, -index)
		end
	end
	return values(d)
end

abstract type AbstractQGramMatchCounter end

abstract type AbstractQGramDistance <: SemiMetric end

function (dist::AbstractQGramDistance)(s1, s2)
	((s1 === missing) | (s2 === missing)) && return missing
	counter = newcounter(dist)
	for (n1, n2) in _count(qgrams(s1, dist.q), qgrams(s2, dist.q))
		count!(dist, counter, n1, n2)
	end
	calculate(dist, counter)
end

"""
	QGram(q::Int)

Creates a QGram distance.

The distance corresponds to

``||v(s1, q) - v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q,
that contains the number of times a q-gram appears for the string s
"""
struct QGram <: AbstractQGramDistance
	q::Int
end

mutable struct SingleCounter <: AbstractQGramMatchCounter
	shared::Int
end

newcounter(::QGram) = SingleCounter(0)
@inline function count!(::QGram, c::SingleCounter, n1::Integer, n2::Integer)
	c.shared += abs(n1 - n2)
end
calculate(::QGram, c::SingleCounter) = c.shared

"""
	Cosine(q::Int)

Creates a Cosine distance.

The distance corresponds to

`` 1 - v(s1, q).v(s2, q)  / ||v(s1, q)|| * ||v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q,
that contains the  number of times a q-gram appears for the string s
"""
struct Cosine <: AbstractQGramDistance
	q::Int
end

mutable struct ThreeCounters <: AbstractQGramMatchCounter
	left::Int
	right::Int
	shared::Int
end

newcounter(::Cosine) = ThreeCounters(0, 0, 0)
@inline function count!(::Cosine, c::ThreeCounters, n1::Integer, n2::Integer)
	c.left += n1^2
	c.right += n2^2
	c.shared += n1 * n2
end
calculate(::Cosine, c::ThreeCounters) =
	1.0 - c.shared / (sqrt(c.left) * sqrt(c.right))

"""
	Jaccard(q::Int)

Creates a Jaccard distance.

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)| / |Q(s1, q) ∪ Q(s2, q))|``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Jaccard <: AbstractQGramDistance
	q::Int
end
newcounter(::Jaccard) = ThreeCounters(0, 0, 0)
@inline function count!(::Jaccard, c::ThreeCounters, n1::Integer, n2::Integer)
	c.left += n1 > 0
	c.right += n2 > 0
	c.shared += (n1 > 0) & (n2 > 0)
end
calculate(::Jaccard, c::ThreeCounters) =
	1.0 - c.shared / (c.left + c.right - c.shared)

"""
	SorensenDice(q::Int)

Creates a SorensenDice distance.

The distance corresponds to

``1 - 2 * |Q(s1, q) ∩ Q(s2, q)|  / (|Q(s1, q)| + |Q(s2, q))|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct SorensenDice <: AbstractQGramDistance
	q::Int
end
newcounter(::SorensenDice) = ThreeCounters(0, 0, 0)
@inline function count!(::SorensenDice, c::ThreeCounters, n1::Integer, n2::Integer)
	c.left += n1 > 0
	c.right += n2 > 0
	c.shared += (n1 > 0) & (n2 > 0)
end
calculate(::SorensenDice, c::ThreeCounters) =
	1.0 - 2.0 * c.shared / (c.left + c.right)

"""
	Overlap(q::Int)

Creates a Overlap distance.

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)|  / min(|Q(s1, q)|, |Q(s2, q)|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Overlap <: AbstractQGramDistance
	q::Int
end
newcounter(::Overlap) = ThreeCounters(0, 0, 0)
@inline function count!(::Overlap, c::ThreeCounters, n1::Integer, n2::Integer)
	c.left += n1 > 0
	c.right += n2 > 0
	c.shared += (n1 > 0) & (n2 > 0)
end
calculate(::Overlap, c::ThreeCounters) =
	1.0 - c.shared / min(c.left, c.right)

"""
	NMD(q::Int)
	NMD(q::Int)

Creates a NMD (Normalized Multiset Distance) as introduced by Besiris and
Zigouris 2013. The goal with this distance is to behave similarly to a normalized
compression distance without having to do any actual compression (and thus being
faster to compute).

The distance corresponds to

``(sum(max.(m(s1), m(s2)) - min(M(s1), M(s2))) / max(M(s1), M(s2))``

where ``m(s)`` is the vector of q-gram counts for string ``s`` and ``M(s)`` is the
sum of those counts.

For details see:
https://www.sciencedirect.com/science/article/pii/S1047320313001417
"""
struct NMD <: AbstractQGramDistance
	q::Int
end

newcounter(::NMD) = ThreeCounters(0, 0, 0)
@inline function count!(::NMD, c::ThreeCounters, n1::Integer, n2::Integer)
	c.left += n1
	c.right += n2
	c.shared += max(n1, n2)
end
calculate(::NMD, c::ThreeCounters) =
	(c.shared - min(c.left, c.right)) / max(c.left, c.right)

"""
	MorisitaOverlap(q::Int)

Creates a MorisitaOverlap distance, a general, statistical measure of
dispersion which can also be used on dictionaries such as created
from q-grams. See https://en.wikipedia.org/wiki/Morisita%27s_overlap_index
This is more fine-grained than many of the other QGramDistances since
it is based on the counts per q-gram rather than only which q-grams are
in the strings.

The distance corresponds to

``(2 * sum(m(s1) .* m(s2)) / (sum(m(s1).^2)*M(s2)/M(s1) + sum(m(s2).^2)*M(s1)/M(s2))``

where ``m(s)`` is the vector of q-gram counts for string ``s`` and ``M(s)`` is the
sum of those counts.
"""
struct MorisitaOverlap <: AbstractQGramDistance
	q::Int
end

mutable struct FiveCounters <: AbstractQGramMatchCounter
	leftsum::Int    # sum(m(s1))
	rightsum::Int   # sum(m(s2))
	leftsq::Int     # sum(m(s1).^2)
	rightsq::Int    # sum(m(s2).^2)
	shared::Int     # sum(m(s1) .* m(s2))
end

newcounter(::MorisitaOverlap) = FiveCounters(0, 0, 0, 0, 0)
@inline function count!(::MorisitaOverlap, c::FiveCounters, n1::Integer, n2::Integer)
	c.leftsum += n1
	c.rightsum += n2
	c.leftsq += n1^2
	c.rightsq += n2^2
	c.shared += n1 * n2
end
calculate(::MorisitaOverlap, c::FiveCounters) =
	1.0 - ((2 * c.shared) / (c.leftsq*c.rightsum/c.leftsum + c.rightsq*c.leftsum/c.rightsum))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SLICOT_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SLICOT")
JLLWrappers.@generate_main_file("SLICOT", UUID("545525a2-e20e-568b-b87f-b40a06098995"))
end  # module SLICOT_jll


text:
@testset "Filter_test" begin

    mt_a = Mod("a")
    mt_b = Mod("b")
    mt_c = Mod("c")
    mt_d = Mod("d")
    mt_e = Mod("e")

    # some test constructs
    oc1 = OrderedConstruct([mt_a mt_d mt_b])
    oc2 = OrderedConstruct([mt_a mt_b mt_e])
    oc3 = OrderedConstruct([mt_e mt_d mt_c])
    oc4 = OrderedConstruct([mt_e mt_a mt_d])
    # constrains
    con1 = UnOrderedConstraint([Mod("a") , Mod("b")])
    con2 = OrderedConstraint([1 3],[Mod("a") Mod("b")])
    con3 = OrderedConstraint([1 3],[Mod("e") Mod("c")])
    con4 = con1 + con2 + con3

    @testset "UnCon" begin
        @test filterconstraint(oc1,con1) == true
        @test filterconstraint(oc2,con1) == true
        @test filterconstraint(oc3,con1) == false
    end

    @testset "OrCon" begin
        @test filterconstraint(oc1,con2) == true
        @test filterconstraint(oc2,con2) == false
        @test filterconstraint(oc3,con2) == false
    end

    @testset "compose" begin
        @test filterconstraint(oc1,con4) == true
        @test filterconstraint(oc2,con4) == true
        @test filterconstraint(oc3,con4) == true
        @test filterconstraint(oc4,con4) == false
    end
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
module LIFSimple
export LIFSimpleAgent, accept, LIFSimpleParams, LIFSimpleStates
using Printf
using ...Types
using ...AgentParts.LIFNeuron
using ...AgentParts.DC
using ...Network
import ...Network: act, state_dict, accept
using ...Signals
import ...Signals: add_acceptor, add_donor, can_add_acceptor, can_add_donor, accept

struct LIFSimpleParams{T <: AbstractFloat}
    lif::LIFParams{T}
    fire_delta_v::T
end

mutable struct LIFSimpleStates{T <: AbstractFloat}
    lif::LIFStates{T}
    accepted_delta_v::T
end

function LIFSimpleStates(v::T,
                         v_steady::T) where {T <: AbstractFloat}
    LIFSimpleStates(LIFStates(v, zero(T), zero(T), v_steady), zero(T))
end

mutable struct LIFSimpleAgent{T <: AbstractFloat, U <: Unsigned} <: Agent
    states::LIFSimpleStates{T}
    params::LIFSimpleParams{T}
    acceptors_delta_v::Vector{Address{U}} # agents that accept from self.
    donors_delta_v::Vector{Address{U}} # agents that donate to self.
    donors_dc::Vector{Address{U}}
end

function LIFSimpleAgent{T, U}(states::LIFSimpleStates{T},
                              params::LIFSimpleParams{T}) where {T <: AbstractFloat, U <: Unsigned}
    LIFSimpleAgent{T, U}(states, params, [], [], [])
    # LIFSimpleAgent{T, U}(states,
    #                      params,
    #                      Vector{Address{U}}(undef, 0),
    #                      Vector{Address{U}}(undef, 0),
    #                      Vector{TimedDeltaV{T}}(undef, 0),
    #                      [])
end

# function update(agent::LIFSimpleAgent{T, U},
#                 updates::LIFSimpleStates{T}) where {T <: AbstractFloat, U <: Unsigned}
#     agent.states = updates
# end

function act(address::Address, # self address.
             agent::LIFSimpleAgent{T, U},
             dt::T) where{T <: AbstractFloat, U <: Unsigned}
    fired, triggered, new_lif_states = evolve(dt,
                                              agent.states.lif,
                                              agent.params.lif,
                                              agent.states.accepted_delta_v)
    triggered_agents = Set{Address{U}}([])
    signals_acceptors = Vector{Tuple{Signal, Vector{Address{U}}}}([])
    if fired
        signal = DeltaV(agent.params.fire_delta_v)
        push!(signals_acceptors, (signal, agent.acceptors_delta_v))
        union!(triggered_agents, Set(agent.acceptors_delta_v))
    end
    if triggered
        push!(triggered_agents, address)
    end

    agent.states = LIFSimpleStates(new_lif_states, zero(T))

    (triggered_agents, signals_acceptors)
end

function accept(agent::LIFSimpleAgent{T, U}, signal::DeltaV{T}) where{T <: AbstractFloat, U <: Unsigned}
    agent.states.accepted_delta_v += signal.delta_v
end

function accept(agent::LIFSimpleAgent{T, U},
                signal::DCInstruction{T}) where{T <: AbstractFloat, U <: Unsigned}
    update_dc(agent.states.lif, agent.params.lif, signal)
end

function can_add_donor(agent::LIFSimpleAgent{T, U},
                       signal_name::String) where{T <: AbstractFloat, U <: Unsigned}
    signal_name == name_delta_v || signal_name == name_dc_instruction
end

function add_donor(agent::LIFSimpleAgent{T, U},
                   signal_name::String,
                   address::Address{U}) where{T <: AbstractFloat, U <: Unsigned}
    if can_add_donor(agent, signal_name)
        if signal_name == name_delta_v
            push!(agent.donors_delta_v, address)
        elseif signal_name == name_dc_instruction
            push!(agent.donors_dc, address)
        else
            error("Got unhandled signal_name on add_donor.")
        end
    else
        error("LIFSimpleAgent cannot add $signal_name for donor of address: $(address.population)-$(address.num)!")
    end
end

function can_add_acceptor(agent::LIFSimpleAgent{T, U},
                          signal_name::String) where{T <: AbstractFloat, U <: Unsigned}
    #println("LIFSimple can_add_acceptor: $(signal_name) vs $(name_t_delta_v), $(signal_name == name_t_delta_v)")
    signal_name == name_delta_v
end

function add_acceptor(agent::LIFSimpleAgent{T, U},
                      signal_name::String,
                      address::Address{U}) where{T <: AbstractFloat, U <: Unsigned}
    if can_add_acceptor(agent, signal_name)
        push!(agent.acceptors_delta_v, address)
    else
        error("LIFSimpleAgent cannot add $signal_name acceptors of address: $(address.population)-$(address.num)!")
    end
end

function state_dict(agent::LIFSimpleAgent{T, U}) where{T <: AbstractFloat, U <: Unsigned}
    Dict(["v" => agent.states.lif.v,
          "t_refractory" => agent.states.lif.t_refractory,
          "dc" => agent.states.lif.dc,
          "v_equilibrium" => agent.states.lif.v_equilibrium])
end

end # module end


text:
using Distributions
using Random

function logistic_population(r, K, p0, t)
  K / (1 + ((K - p0) / p0) * exp(-r * t))
end

function logistic_difeq(x)
    P, r, K = x
    r * P * (1 - P / K)
end

@doc raw"""
    beverton_holt(x, K, nu, mu)

Given a population of size x, generate a population at the next time step.
`K` is the limiting population. `nu` is the growth rate. The Beverton-Holt
model normally has just `nu` as a growth rate. This one adds `mu` as
a rate at which the system decreases to the limiting population.
"The Beverton-Holt dynamic equation" by Bohner and Warth, 2007
"""
function beverton_holt(x, K, nu, mu)
    @assert nu > 1
    @assert mu > 1
    if x <= K
        rate = nu
    else
        rate = mu
    end
    (rate * K * x) / (K + (rate - 1) * x)
end

@doc raw"""
    beverton_holt_constant_food(young, step_cnt, params)

Predict `step_cnt` time steps of a future population, starting with `young`
individuals and using parameters in the `params` dictionary. In this version,
it reads the amoung of food from the parameter called `:max_food`.
"""
function beverton_holt_constant_food(young, step_cnt, params)
    x = zeros(Float64, step_cnt)
    x[1] = young
    for step_idx = 1:(step_cnt - 1)
        x[step_idx + 1] = beverton_holt(
            x[step_idx], params[:max_food], params[:growth_rate], params[:death_rate]
            )
    end
    x
end

function beverton_holt_carrying(young, carrying, growth_rate, step_cnt)
    x = zeros(Float64, step_cnt)
    x[1] = young
    for step_idx = 1:(step_cnt - 1)
        x[step_idx + 1] = ((growth_rate * carrying[step_idx] * x[step_idx]) /
            (carrying[step_idx] + (growth_rate - 1) * x[step_idx]))
    end
    x
end

@doc raw"""
    forecast_population(young, step_cnt, params, rng)

Predict `step_cnt` time steps of a future population, starting with `young`
individuals and using parameters in the `params` dictionary. `rng` is
a random number generator.
"""
function forecast_population(young, step_cnt, params, rng)
    food = zeros(Float64, step_cnt)
    rand!(rng, Distributions.Beta(params[:food_alpha], params[:food_beta]), food)
    food .*= params[:max_food]

    x = zeros(Float64, step_cnt)
    x[1] = young
    for step_idx = 1:(step_cnt - 1)
        x[step_idx + 1] = beverton_holt(x[step_idx], food[step_idx], params[:growth_rate], params[:death_rate])
    end
    x
end

function run_pop()
    rng = MersenneTwister(9237429)
    parameters = Dict(
        :max_food => 400,
        :food_alpha => 2,
        :food_beta => 10,
        :growth_rate => 1.02,
        :death_rate => 1.02
    )
    forecast_population(2, 5000, parameters, rng)
end

function paperfig(fn)
    savefig(fn * ".pdf")
    savefig(fn * ".png")
end

function make_plots()
    y = run_pop()

    histogram(y[1000:5000])
    paperfig("insects_histogram_4k")

    plot(y[1:500])
    paperfig("insects_rise")

    plot(y[1:5000])
    paperfig("insects_long")
end

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
using PkgBenchmark

results = benchmarkpkg("McCormick")
show(results)

# specify tag and uncommit to benchmark versus prior tagged version
#tag = "9b5ba7d756e3feac2a9ab5d42af464f1139b8900" # v0.4.1
#results = judge("McCormick", tag)
#show(results)

export_markdown("results.md", results)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using PyPlot
#assumption the price each wants to sell/buy is determined by Uniform Distribution

Num_suppliers = 20
Num_consumers = 20

type Consumer
  id::Int
  u::Float64
  p::Float64
  buy::Function
  bought::Bool
  function Consumer(id)
    C = new()
    C.id = id
    C.u = 0
    C.p = rand()
    C.bought = false
    return C
  end
end

whether_to_buy(c, p) = c.p >= p
function buy(c, p)
  if whether_to_buy(c, p)
    #c.bought = true
    c.u += c.p - p
  end
end

all_u_of_eles(eles) = sum(map(e -> e.u, eles))

type Supplier
  id::Int
  p::Float64
  u::Float64
  sell::Function
  sold::Bool
  function Supplier(id)
    S = new()
    S.id = id
    S.u = 0
    S.p = rand()
    S.sold = false
    return S
  end
end

function interact(c, s)
  c.p > s.p ? true : false
end

whether_to_sell(s, p) = s.p <= p
function sell(s, p)
  if whether_to_sell(s, p)
    #s.sold = true
    s.u += p - s.p
  end
end

function calc_num(whether, eles, ps)
  e4p = []
  for p in ps
    push!(e4p, sum(map(e -> whether(e, p), eles)))
  end
  return e4p
end

function find_equilibrium(ps, supplies, demands)
  i = indmin(map(abs, supplies .- demands))
  return ps[i], min(supplies[i], demands[i])
end

function trade(p_s, p_c, ss, cs)
  map(s -> sell(s, p_s), ss)
  map(c -> buy(c, p_c), cs)
end

function show_data(x, p_s, p_c, ss, cs)
  println("under free trade")
  println("suppliers u: $(all_u_of_eles(ss))")
  println("consumers u: $(all_u_of_eles(cs))")
  println("all u: $(all_u_of_eles(ss) + all_u_of_eles(cs))")
  println("traded amount: $x")
  println("price for consumers: $p_c, price for suppliers: $p_s")
end

calc_supply(ss, ps) = calc_num(whether_to_sell, ss, ps)
calc_demand(cs, ps) = calc_num(whether_to_buy, cs, ps)

cs = [Consumer(i) for i in 1:Num_consumers]
ss = [Supplier(i) for i in 1:Num_consumers]

println(ss)
println(cs)

ps = linspace(0, 1, 100)
supplies = calc_supply(ss, ps)
demands = calc_demand(cs, ps)
plot(supplies, ps)
plot(demands, ps)
p, x = find_equilibrium(ps, supplies, demands)

trade(p, p, ss, cs)
show_data(x, p, p, ss, cs)
"""
println("under taxation")
println("suppliers u: $(all_u_of_eles(ss))")
println("consumers u: $(all_u_of_eles(cs))")
println("all u: $(all_u_of_eles(ss) + all_u_of_eles(cs))")
"""


text:
# Solves diffusion equation for zoned crystal using Crank-Nicolson scheme

using LinearAlgebra

Nt = 100
t_end = 1.0
dt = t_end/(Nt-1)
t1 = collect(LinRange(70.0,50.0,ceil(Int,Nt/2)).*3.1558e7*1e6)
t2 = collect(LinRange(49.9,0.01,floor(Int,Nt/2)).*3.1558e7*1e6)
t = vcat(t1,t2)

T1 = collect(LinRange(900.0,600.0,ceil(Int,Nt/2)).+273.15)
T2 = collect(LinRange(600.0,400.0,floor(Int,Nt/2)).+273.15)
T = vcat(T1,T2)
Ea = 250.0*1e3
D0 = 3.9*1e-10
#Rjoules = 8.314
L = 1e-4
Nx = 513
U238 = 1.0
# Nt,t_end
function zonation_forward(Ea,R,D0,U235,Th232,L,Nt,Nx,tTrU...)

tTcopy = collect(tTrU)
t = tTcopy[1:Nt]#[1:ceil(Int,length(tT)/2)] # type unstable but fastest so far
T = tTcopy[(Nt+1):(2*Nt)]#[ceil(Int,length(tT)/2)+1:end]
r = tTcopy[(2*Nt+1):(2*Nt+Nx)]
U238 = tTcopy[(2*Nt+Nx+1):end]
U238_0 = U238*(exp(COAST.lambda_38*t[1])) # U238 is measured at present
U235_0 = U235*(exp(COAST.lambda_35*t[1])) # U238 is measured at present
Th232_0 = Th232*(exp(COAST.lambda_32*t[1])) # U238 is measured at present
D_diffn = [0.0]

Pb = 0.0*ones(Nx)
Pb_solns = zeros(Nx,Nt) # includes time t0
Pb_solns[:,1] = Pb
Pb_rad = 0.0

rhs = zeros(Nx)
U238_cur = U238_0
for i in 1:Nt-1

    dt = t[i] - t[i+1]
    U238_cur = U238_0*(exp(-COAST.lambda_38*(t[1]-t[i]+dt/2))) # Crank-Nicolson is staggered in time, hence dt/2
    D_diffn[1] = D0*exp(-Ea/(R*T[i]))
    A = D_diffn[1]*dt/(dr^2)
    B = D_diffn[1]*dt/(2*dr)
    lhs_center = [1.0 + 3*A; fill(1.0 + A, Nx-2);1.0]
    lhs_lo = [-(A/2).+ B./r[2:Nx-1];0.0]
    lhs_hi = [-3*A; -A/2 .- B./r[2:Nx-1]]
    lhs = LinearAlgebra.Tridiagonal(lhs_lo,lhs_center,lhs_hi)
    rhs[1] = (1.0-3*A)*Pb[1] + 3*A*Pb[2] + U238_cur[1]*COAST.lambda_38*dt
    for j in 2:(Nx-1)
        rhs[j] = (A/2 - B/r[j])*Pb[j-1] + (1.0 - A)*Pb[j] + (A/2 + B/r[j])*Pb[j+1] + U238_cur[j]*COAST.lambda_38*dt
    end
    rhs[end] = Pb_rad

    Pb[:] = lhs\rhs
    Pb_solns[:,i+1] = Pb

end

return Pb_solns[:,end]
end


text:
using Test

function sum(arr)
    if isempty(arr)
        return 0
    end
    return arr[1] + sum(arr[2:end])
end

@test sum([1,2,3,4]) == 10


text:
import IterativeSolvers: gmres!, bicgstabl!
import LinearAlgebra: ldiv!

abstract type AbstractLinearSolver end
abstract type DirectSolver <: AbstractLinearSolver end
abstract type IterativeSolver <: AbstractLinearSolver end

struct LU  <:DirectSolver end
struct GMRES  <: IterativeSolver end
struct BICGSTAB  <: IterativeSolver end

linsolve(A,b ; args...) = linsolve(A,b,LU() ; args...)

function linsolve(A,b,solver::AbstractLinearSolver ; args...)
    x = similar(b)
    linsolve!(x,A,b,solver; args...)
    return x
end

linsolve!(x,A,b ; args...) = linsolve!(x,A,b,LU() ; args...)

function linsolve!(x,A,b , ::LU ; args...)
  ldiv!(x, lu(A), b; args...)
end

function linsolve!(x,A,b,::GMRES ; args...)
    gmres!(x,A,b; args...)
end

function linsolve!(x,A,b,::BICGSTAB ; args...)
    bicgstabl!(x,A,b; args...)
end

# methods for sensitivty analysis
function linsolve!(x,x_adj,A,b,b_adj,::LU ; args...)
    F = lu(A)
    ldiv!(x, F, b; args...)
    ldiv!(x_adj, transpose(F), b_adj; args...)
end

function linsolve(A,b,b_adj,::LU ;  args...)
    x = similar(b) ; x_adj = similar(b_adj)
    linsolve!(x,x_adj,A,b,b_adj,LU(); args...)
    return x , x_adj
end


text:
using Documenter, DocumenterTools, Juqbox

const ROOT = joinpath(@__DIR__, "..")

makedocs(
    modules = [Juqbox],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    #format = Documenter.LaTeX(platform = "docker"), # errors with no such file or directory???
    #clean = false,
    sitename="Juqbox.jl",
    authors = "Anders Petersson, Fortino Garcia, and contributors.",
    pages = [
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Workflow" => "workflow.md",
        "Examples" => "examples.md",
        "Types" => "types.md",
        "Methods" => "methods.md",
        "Index" => "function-index.md",
    ],
    doctest = false, # set to true once initial implementation is done
    source = "src",
)

deploydocs(
    repo = "github.com/LLNL/Juqbox.jl",
    target = "build",
    devurl="docs",
)


text:
struct DijkstraHeapEntry{T, U<:Integer}
    vertex::U
    dist::T
end

isless(e1::DijkstraHeapEntry, e2::DijkstraHeapEntry) = e1.dist < e2.dist

"""
    struct DijkstraState{T, U}

An [`AbstractPathState`](@ref) designed for Dijkstra shortest-paths calculations.
"""
struct DijkstraState{T, U<:Integer}<: AbstractPathState
    parents::Vector{U}
    dists::Vector{T}
    predecessors::Vector{Vector{U}}
    pathcounts::Vector{U}
end

"""
    dijkstra_shortest_paths(g, srcs, distmx=DefaultDistance());

Perform [Dijkstra's algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)
on a graph, computing shortest distances between `srcs` and all other vertices.
Return a [`DijkstraState`](@ref) that contains various traversal information.

### Optional Arguments
- `allpaths=false`: If true, returns a [`DijkstraState`](@ref) that keeps track of all
predecessors of a given vertex.
"""
function dijkstra_shortest_paths(
    g::AbstractGraph,
    srcs::Vector{U},
    distmx::AbstractMatrix{T}=DefaultDistance();
    allpaths=false
    ) where T where U<:Integer
    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    preds = fill(Vector{U}(),nvg)
    visited = zeros(Bool, nvg)
    pathcounts = zeros(Int, nvg)
    H = Vector{DijkstraHeapEntry{T, U}}()  # this should be Vector{T}() in 0.4, I think.
    dists[srcs] = zero(T)
    pathcounts[srcs] = 1

    sizehint!(H, nvg)

    for v in srcs
        heappush!(H, DijkstraHeapEntry{T, U}(v, dists[v]))
        visited[v] = true
    end

    while !isempty(H)
        hentry = heappop!(H)
            # info("Popped H - got $(hentry.vertex)")
        u = hentry.vertex
        for v in out_neighbors(g,u)
            alt = (dists[u] == typemax(T))? typemax(T) : dists[u] + distmx[u,v]

            if !visited[v]
                dists[v] = alt
                parents[v] = u
                pathcounts[v] += pathcounts[u]
                visited[v] = true
                if allpaths
                    preds[v] = [u;]
                end
                heappush!(H, DijkstraHeapEntry{T, U}(v, alt))
                # info("Pushed $v")
            else
                if alt < dists[v]
                    dists[v] = alt
                    parents[v] = u
                    heappush!(H, DijkstraHeapEntry{T, U}(v, alt))
                end
                if alt == dists[v]
                    pathcounts[v] += pathcounts[u]
                    if allpaths
                        push!(preds[v], u)
                    end
                end
            end
        end
    end

    pathcounts[srcs] = 1
    parents[srcs] = 0
    for src in srcs
        preds[src] = []
    end

    return DijkstraState{T, U}(parents, dists, preds, pathcounts)
end

dijkstra_shortest_paths(g::AbstractGraph, src::Integer, distmx::AbstractMatrix = DefaultDistance(); allpaths=false) =
dijkstra_shortest_paths(g, [src;], distmx; allpaths=allpaths)

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
export SACPolicy

mutable struct SACPolicy{
    BA<:NeuralNetworkApproximator,
    BC1<:NeuralNetworkApproximator,
    BC2<:NeuralNetworkApproximator,
    P,
    R<:AbstractRNG,
} <: AbstractPolicy
    policy::BA
    qnetwork1::BC1
    qnetwork2::BC2
    target_qnetwork1::BC1
    target_qnetwork2::BC2
    γ::Float32
    τ::Float32
    α::Float32
    batch_size::Int
    start_steps::Int
    start_policy::P
    update_after::Int
    update_freq::Int
    automatic_entropy_tuning::Bool
    lr_alpha::Float32
    target_entropy::Float32
    update_step::Int
    rng::R
    # Logging
    reward_term::Float32
    entropy_term::Float32
end

"""
    SACPolicy(;kwargs...)

# Keyword arguments

- `policy`, used to get action.
- `qnetwork1`, used to get Q-values.
- `qnetwork2`, used to get Q-values.
- `target_qnetwork1`, used to estimate the target Q-values.
- `target_qnetwork2`, used to estimate the target Q-values.
- `start_policy`,
- `γ::Float32 = 0.99f0`, reward discount rate.
- `τ::Float32 = 0.005f0`, the speed at which the target network is updated.
- `α::Float32 = 0.2f0`, entropy term.
- `batch_size = 32`,
- `start_steps = 10000`,
- `update_after = 1000`,
- `update_freq = 50`,
- `automatic_entropy_tuning::Bool = false`, whether to automatically tune the entropy.
- `lr_alpha::Float32 = 0.003f0`, learning rate of tuning entropy.
- `action_dims = 0`, the dimensionality of the action. if `automatic_entropy_tuning = true`, must enter this parameter.
- `update_step = 0`,
- `rng = Random.GLOBAL_RNG`,

`policy` is expected to output a tuple `(μ, logσ)` of mean and
log standard deviations for the desired action distributions, this
can be implemented using a `GaussianNetwork` in a `NeuralNetworkApproximator`.

Implemented based on http://arxiv.org/abs/1812.05905
"""
function SACPolicy(;
    policy,
    qnetwork1,
    qnetwork2,
    target_qnetwork1,
    target_qnetwork2,
    start_policy,
    γ = 0.99f0,
    τ = 0.005f0,
    α = 0.2f0,
    batch_size = 32,
    start_steps = 10000,
    update_after = 1000,
    update_freq = 50,
    automatic_entropy_tuning = true,
    lr_alpha = 0.003f0,
    action_dims = 0,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(qnetwork1, target_qnetwork1)  # force sync
    copyto!(qnetwork2, target_qnetwork2)  # force sync
    if automatic_entropy_tuning
        @assert action_dims != 0
    end
    SACPolicy(
        policy,
        qnetwork1,
        qnetwork2,
        target_qnetwork1,
        target_qnetwork2,
        γ,
        τ,
        α,
        batch_size,
        start_steps,
        start_policy,
        update_after,
        update_freq,
        automatic_entropy_tuning,
        lr_alpha,
        Float32(-action_dims),
        update_step,
        rng,
        0f0,
        0f0,
    )
end

# TODO: handle Training/Testing mode
function (p::SACPolicy)(env)
    p.update_step += 1

    if p.update_step <= p.start_steps
        p.start_policy(env)
    else
        D = device(p.policy)
        s = state(env)
        s = Flux.unsqueeze(s, ndims(s) + 1)
        # trainmode:
        action = dropdims(p.policy(p.rng, s; is_sampling=true), dims=2) # Single action vec, drop second dim

        # testmode:
        # if testing dont sample an action, but act deterministically by
        # taking the "mean" action
        # action = dropdims(p.policy(s)[1], dims=2)
    end
end

function RLBase.update!(
    p::SACPolicy,
    traj::CircularArraySARTTrajectory,
    ::AbstractEnv,
    ::PreActStage,
)
    length(traj) > p.update_after || return
    p.update_step % p.update_freq == 0 || return
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(p.batch_size))
    update!(p, batch)
end

function RLBase.update!(p::SACPolicy, batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(p.qnetwork1), batch)

    γ, τ, α = p.γ, p.τ, p.α

    a′, log_π = p.policy(p.rng, s′; is_sampling=true, is_return_log_prob=true)
    q′_input = vcat(s′, a′)
    q′ = min.(p.target_qnetwork1(q′_input), p.target_qnetwork2(q′_input))

    y = r .+ γ .* (1 .- t) .* vec(q′ .- α .* log_π)

    # Train Q Networks
    q_input = vcat(s, a)

    q_grad_1 = gradient(Flux.params(p.qnetwork1)) do
        q1 = p.qnetwork1(q_input) |> vec
        mse(q1, y)
    end
    update!(p.qnetwork1, q_grad_1)
    q_grad_2 = gradient(Flux.params(p.qnetwork2)) do
        q2 = p.qnetwork2(q_input) |> vec
        mse(q2, y)
    end
    update!(p.qnetwork2, q_grad_2)

    # Train Policy
    p_grad = gradient(Flux.params(p.policy)) do
        a, log_π = p.policy(p.rng, s; is_sampling=true, is_return_log_prob=true)
        q_input = vcat(s, a)
        q = min.(p.qnetwork1(q_input), p.qnetwork2(q_input))
        reward = mean(q)
        entropy = mean(log_π)
        ignore() do
            p.reward_term = reward
            p.entropy_term = entropy
        end
        α * entropy - reward
    end
    update!(p.policy, p_grad)

    # Tune entropy automatically
    if p.automatic_entropy_tuning
        p.α -= p.lr_alpha * mean(-log_π .- p.target_entropy)
    end

    # polyak averaging
    for (dest, src) in zip(
        Flux.params([p.target_qnetwork1, p.target_qnetwork2]),
        Flux.params([p.qnetwork1, p.qnetwork2]),
    )
        dest .= (1 - τ) .* dest .+ τ .* src
    end
end


text:
using Cytof5, RCall
import Cytof5.Util.similarity_FM
@rimport graphics as plt
@rimport grDevices as dev

# Get path to results
if length(ARGS) == 0
  PATH_TO_RESULTS = "results/cb-paper/"
else
  PATH_TO_RESULTS = ARGS[1]
end
results_dirs = filter(d -> occursin("K_MCMC", d), readdir(PATH_TO_RESULTS))

struct Estimate
  Z::Vector{Matrix{Int}}
  W::Vector{Vector{Float64}}
  I::Int
  J::Int
  K::Int
end

function Estimate(Z::Vector{Matrix{Int}}, W::Vector{Vector{Float64}})
  @assert length(Z) == length(W)
  I = length(Z)
  J, K = size(Z[1])
  for i in 1:I
    @assert size(Z[i], 2) == length(W[i])
    @assert size(Z[i], 1) == J
  end
  return Estimate(Z, W, I, J, K)
end

function read_Zhat(path_to_zhat)
  rows = Vector{Int}[]
  open(path_to_zhat, "r") do f
    lines = readlines(f)
    for line in lines
      row = parse.(Int, split(line, ","))
      append!(rows, [row])
    end
  end
  Z = Matrix(hcat(rows...)')
  return Z
end

function read_What(path_to_what)
  what = nothing
  open(path_to_what, "r") do f
    lines = readlines(f)
    what = parse.(Float64, lines)
  end
  return what
end

function get_est(d)
  imgdir = "$(PATH_TO_RESULTS)/$(d)/img/"
  Z_hat_files = filter(d -> occursin(r"Z\d+_hat\.txt", d), readdir(imgdir))
  W_hat_files = filter(d -> occursin(r"W_\d+_hat\.txt", d), readdir(imgdir))
  I = length(Z_hat_files)
  @assert I == length(W_hat_files)

  Z_hats = [read_Zhat("$imgdir/$d") for d in Z_hat_files]
  W_hats = [read_What("$imgdir/$d") for d in W_hat_files]

  Estimate(Z_hats, W_hats)
end

# Read estimates
estimates = [get_est(d) for d in results_dirs]
ord = sortperm(getfield.(estimates, :K))
estimates = estimates[ord]
Ks = [est.K for est in estimates[2:end]]

z_metric = Float64[]
for r in 2:length(estimates)
  W = Matrix(hcat(estimates[r].W...)')
  zm = similarity_FM(estimates[r].Z[1], W, estimates[r-1].Z[1])
  append!(z_metric, zm)
end

dev.pdf("$(PATH_TO_RESULTS)/metrics/L0_MCMC5/z-metric.pdf")
plt.plot(Ks, z_metric, typ="o", xlab="K-MCMC", ylab="Z-metric", xaxt="n")
plt.axis(1, at=Ks, label=Ks)
dev.dev_off()

z_metric_scaled = Float64[]
for r in 2:length(estimates)
  W = Matrix(hcat(estimates[r].W...)')
  zm = similarity_FM(estimates[r].Z[1], W, estimates[r-1].Z[1], z_diff=Cytof5.Util.meanabsdiff)
  append!(z_metric_scaled, zm)
end

dev.pdf("$(PATH_TO_RESULTS)/metrics/L0_MCMC5/z-metric-scaled.pdf")
plt.plot(Ks, z_metric_scaled, typ="o", xlab="K-MCMC", ylab="Z-metric", xaxt="n")
plt.axis(1, at=Ks, label=Ks)
dev.dev_off()


text:
import AndTraits
export Variable, AbstractVariable, isvariable, traitvartype, AbstractClass, AbstractVariableOrClass, TraitIsVariableOrClass

abstract type AbstractVariableOrClass end

abstract type AbstractVariable <: AbstractVariableOrClass end

abstract type AbstractClass <: AbstractVariableOrClass end

# Traits
const TraitIsVariable = AndTraits.Trait{:TraitIsVariable}
const TraitIsClass = AndTraits.Trait{:TraitIsClass}
const TraitIsVariableOrClass = AndTraits.Trait{:TraitIsVariableOrClass}
const TraitUnknownVariableType = AndTraits.Trait{:TraitUnknownVariableType}

"Is `v` a variable?"
function isvariable(v::Function)
  mthds = methods(v)
  any(mthds.ms) do m
    # Core.println(m.sig)
    !(m.sig isa UnionAll) && length(m.sig.types) == 2 && (m.sig.types[2] == AbstractΩ)
  end
end

function isclass(v::Function)
  mthds = methods(v)
  any(mthds.ms) do m
    # Core.println(m.sig)
    !(m.sig isa UnionAll) && length(m.sig.types) == 3 && (m.sig.types[3] == AbstractΩ)
  end
end

function traitvartype_(f::Function)
  # Core.println(f)
  if isvariable(f)
    TraitIsVariable
  elseif isclass(f)
    TraitIsClass
  else
    TraitUnknownVariableType
  end
end

# @generated function traitvartype(f::Function)
#   traitvartype(f)
# end

@generated function traitvartype(f::Type{<:Function})
  functype = f.parameters[1]
  traitvartype_(functype.instance)
end

# By default we don't know the variable type
traitvartype(T) = TraitUnknownVariableType

traitvartype(::AbstractVariable) = TraitIsVariable
traitvartype(::Type{<:AbstractVariable}) = TraitIsVariable

# # Variable
# A Variable is a parametric or random variable, which is just any function of
# ω::AbstractΩ.  We have this data structure because we need to intercept inner
# calls to `f(ω)`, which is difficult if it is not a special type and just a
# normal function

"A variable is just a function of ω."
struct Variable{F} <: AbstractVariable
  f::F
end
recurse(f::Variable, ω) = f.f(ω)

Base.show(io::IO, f::Variable) = Base.print(io, "ᵛ", f.f)

struct Class{C} <: AbstractClass
  c::C
end
Var.traitvartype(::Class) = TraitIsClass
Var.traitvartype(::Type{<:Class}) = TraitIsClass

@inline (c::Class)(i, ω) = c.c(i, ω)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(
        inradius,
        outradius,
        center,
        ls,
        ms,
        lc,
        mc,
        theta0,
    )
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function shell_radial_stress(ls, ms, theta0, A1, A2, r)
    return (ls + 2ms) * (A1 - A2 / r^2) + ls * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_circumferential_stress(ls, ms, theta0, A1, A2, r)
    return ls * (A1 - A2 / r^2) + (ls + 2ms) * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_out_of_plane_stress(ls, ms, A1, theta0)
    return 2 * ls * A1 - (ls + 2ms / 3) * theta0
end

function core_in_plane_stress(lc, mc, A1)
    return (lc + 2mc) * A1 + lc * A1
end

function core_out_of_plane_stress(lc, A1)
    return 2 * lc * A1
end

function rotation_matrix(x, r)
    costheta = x[1] / r
    sintheta = x[2] / r
    Q = [
        costheta -sintheta
        sintheta costheta
    ]
    return Q
end

function shell_radial_stress(A::AnalyticalSolution, r)
    return shell_radial_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
end

function shell_circumferential_stress(A::AnalyticalSolution, r)
    return shell_circumferential_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
end

function shell_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    Q = rotation_matrix(relpos, r)

    srr = shell_radial_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
    stt = shell_circumferential_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)

    cylstress = [
        srr 0.0
        0.0 stt
    ]

    cartstress = Q * cylstress * Q'
    s11 = cartstress[1, 1]
    s22 = cartstress[2, 2]
    s12 = cartstress[1, 2]
    s33 = shell_out_of_plane_stress(A.ls, A.ms, A.A1s, A.theta0)

    return [s11, s22, s12, s33]
end

function core_in_plane_stress(A::AnalyticalSolution)
    return core_in_plane_stress(A.lc, A.mc, A.A1c)
end

function core_stress(A::AnalyticalSolution)
    s11 = core_in_plane_stress(A.lc, A.mc, A.A1c)
    s33 = core_out_of_plane_stress(A.lc, A.A1c)
    return [s11, s11, 0.0, s33]
end

function exact_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    if r < A.inradius
        return core_stress(A)
    else
        return shell_stress(A, x)
    end
end

function nodal_displacement(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    boundarydisplacement,
    penalty;
    eta = 1,
)

    L, W = CutCellDG.widths(mesh)
    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        basis,
        cellquads,
        stiffness,
        mesh,
    )
    CutCellDG.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        mesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_interelement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
    )

    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        basis,
        interfacequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        mesh,
    )

    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        boundarydisplacement,
        basis,
        facequads,
        stiffness,
        mesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
        x -> onboundary(x, L, W),
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mesh)

    solution = matrix \ rhs

    return solution
end

function construct_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    interfacecenter,
    interfaceradius,
    numqp;
    tinyratio = 0.2,
)
    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)

    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCellDG.levelset_coefficients(
        x -> -circle_distance_function(x, interfacecenter, interfaceradius),
        mesh,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset, levelsetcoeffs)
    cellquads =
        CutCellDG.CellQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    interfacequads =
        CutCellDG.InterfaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    facequads =
        CutCellDG.FaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)

    mergedwithcell, hasmergedcells = CutCellDG.merge_tiny_cells_in_mesh!(
        cutmesh,
        cellquads,
        facequads,
        interfacequads,
        tinyratio = tinyratio,
    )
    mergedmesh = CutCellDG.MergedMesh(cutmesh, mergedwithcell)

    return mergedmesh, cellquads, facequads, interfacequads
end

function construct_unmerged_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    interfacecenter,
    interfaceradius,
    numqp,
)
    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)

    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCellDG.levelset_coefficients(
        x -> -circle_distance_function(x, interfacecenter, interfaceradius),
        mesh,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset, levelsetcoeffs)
    cellquads =
        CutCellDG.CellQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    interfacequads =
        CutCellDG.InterfaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    facequads =
        CutCellDG.FaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)

    return cutmesh, cellquads, facequads, interfacequads
end

function product_stress(
    celldisp,
    basis,
    stiffness,
    transfstress,
    theta0,
    point,
    jac,
    vectosymmconverter,
)

    dim = length(vectosymmconverter)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, +1)

    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp

    inplanestress = (stiffness[+1] * symmdispgrad) - transfstress
    s33 =
        lambda * (symmdispgrad[1] + symmdispgrad[2]) -
        (lambda + 2mu / 3) * theta0

    stress = vcat(inplanestress, s33)

    return stress
end

function product_stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    transfstress,
    theta0,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    nphase, dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert size(referencecellids) == (nphase, numpts)

    productstress = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    row = CutCellDG.cell_sign_to_row(+1)

    for i = 1:numpts
        cellid = referencecellids[row, i]
        nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[row, :, i]

        productstress[:, i] .= product_stress(
            celldisp,
            basis,
            stiffness,
            transfstress,
            theta0,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return productstress
end

function parent_stress(
    celldisp,
    basis,
    stiffness,
    point,
    jac,
    vectosymmconverter,
)
    dim = length(vectosymmconverter)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, -1)

    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp

    inplanestress = stiffness[-1] * symmdispgrad
    s33 = lambda * (symmdispgrad[1] + symmdispgrad[2])

    stress = vcat(inplanestress, s33)

    return stress
end

function parent_stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    nphase, dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert size(referencecellids) == (nphase, numpts)

    parentstress = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    row = CutCellDG.cell_sign_to_row(-1)

    for i = 1:numpts
        cellid = referencecellids[row, i]
        nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[row, :, i]

        parentstress[:, i] .= parent_stress(
            celldisp,
            basis,
            stiffness,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return parentstress
end


text:
# named tests
let
  @pytest function test_name()
    @test 1 == 1
  end
end

# nesting @pytest in @testset
module basetest_integration  # for some reason cant use let here
  using PyTest
  using Base.Test
  using ResumableFunctions
  import Base.Test.DefaultTestSet
  @fixture params=(100, 200) function h(request) request.param end
  @fixture function f(h) h end
  @fixture function g(f) return (f, :g, ) end
  tests_performed = Dict("key" => 0)
  const number_v = 5
  @testset DefaultTestSet "testset description $v, $w" for w in [:a, :b], v in 1:number_v
    @pytest function pytest_description(g, h)
      @test g == (h, :g, )
      tests_performed["key"] = tests_performed["key"] + 1
    end
  end
  @test tests_performed["key"] == number_v * 2 * 2
end


text:
module wiring

export manhattan_distance, signal_delay

struct Point
    x::Int64
    y::Int64
end

function manhattan_distance(left, right)
    intersections = intersect(keys(trace(left)), keys(trace(right)))
    return minimum(map(p -> abs(p.x) + abs(p.y), collect(intersections)))
end

function signal_delay(left, right)
    l = trace(left)
    r = trace(right)
    intersections = intersect(keys(l), keys(r))
    return minimum(
        map(location -> l[location] + r[location],
            collect(intersections)
        )
    )
end

function trace(ops)
    x = 0
    y = 0
    acc::Dict{Point, Int64} = Dict()
    counter = 0
    for op in ops
        direction = SubString(op, 1, 1)
        distance = parse(Int, SubString(op, 2))
        for i = 1:distance
            location = move(direction, x, y)
            x = location.x
            y = location.y
            counter += 1
            if !haskey(acc, location)
                acc[location] = counter
            end
        end
    end
    return acc
end

function move(direction, x, y)
    if direction == "R"
        return Point(x + 1, y)
    elseif direction == "L"
        return Point(x - 1, y)
    elseif direction == "U"
        return Point(x, y + 1)
    elseif direction == "D"
        return Point(x, y - 1)
    end
end

end


text:
include("testdata.jl")

@testset "SGcoef" begin
    c = coef(Linear(sgtabX, sgtabY))
    @test c.slope ≈  1.08862150
    @test c.intercept ≈ -2.63197077
end

@testset "NEATcoef" begin
    c = coef(ChainedLinear(neattabX, neattabY))
    @test c.slope ≈ 1.0212716885774726
    @test c.intercept ≈ 0.39367983874242185
end


text:
using Plasma

TD = 30000 # eV
Te = 10000 # eV

D = species.D
e = species.e

D_D = Distribution(Maxwellian(TD, D.m), D)
D_e = Distribution(Maxwellian(Te, e.m), e)
G = Geometry()

plasma = ElectrostaticPlasma([D_D, D_e], G)

sol = Plasma.solve(plasma, dim=1, GPU=false)

Plasma.plot(sol)

## 2D with custom P and species
using Plasma

Tα = 70000 # eV

α = Species(1.602176634e-19, 6.6446562e-27)

function HotCarrier(T)
    Kb = 8.617333262145e-5
    P(x,v) = exp(-v/(Kb*T))
end

Dα = Distribution(HotCarrier(Tα), α)
G = Geometry() # TODO define a custom geometry

plasma = ElectrostaticPlasma([Dα], G)

Plasma.solve(plasma, dim=2) # with GPU

Plasma.plot(sol)

## 3D CollisionlessPlasma
using Plasma

TD = 15000 # eV
TT = 15000 # eV
Te = 13000 # eV

e = species.e
T = species.T
D = species.D

De = Distribution(Maxwellian(Te, e.m), e)
DT = Distribution(Maxwellian(TT, T.m), T)
DD = Distribution(Maxwellian(TD, D.m), D)
G = Geometry()

plasma = CollisionlessPlasma([De,DT,DD], G)

Plasma.solve(plasma)

Plasma.plot(sol)


text:
using Revise
using AdFem
using PyPlot
using SparseArrays

cmesh = CrackMesh(50,50,1/50,15)
mmesh = cmesh.mesh

NT = 500
Δt = 2/NT

Ca = spzeros(2cmesh.nslip, 2mmesh.ndof)
ls = zeros(2cmesh.nslip)
for i = 1:cmesh.nslip
    u = cmesh.upper[i+1]
    l = cmesh.lower[i+1]
    Ca[i, u] = 1.0
    Ca[i, l] = -1.0
    ls[i] = 0.0

    u = cmesh.upper[i+1] + mmesh.ndof
    l = cmesh.lower[i+1] + mmesh.ndof
    Ca[i+cmesh.nslip, u] = 1.0
    Ca[i+cmesh.nslip, l] = -1.0
    ls[i+cmesh.nslip] = 0.0
end
Ca = constant(Ca)

upbd = bcnode((x,y)->y>1.0-0.001, mmesh)
downbd = bcnode((x,y)->y<0.001, mmesh)

ν = 0.1
E = 1.0
D = compute_plane_strain_matrix(E, ν)
K = constant(compute_fem_stiffness_matrix(D, mmesh))

rhs = zeros(2mmesh.ndof)
K, rhs = impose_Dirichlet_boundary_conditions(
    K, rhs,
    [upbd; downbd; upbd .+ mmesh.ndof; downbd .+ mmesh.ndof],
    [0.001*ones(length(upbd)); zeros(length(downbd)+length(upbd)+length(downbd))]
)

A = [K Ca';Ca spzero(2cmesh.nslip)]
rhs = [rhs; constant(ls)]
u = A\rhs

u = u[1:2mmesh.ndof]

sess = Session(); init(sess)

U = run(sess, u)

# close("all")
# visualize_displacement(100*U, mmesh)
# savefig("mesh.png")

Se = compute_von_mises_stress_term(D, U, mmesh)
close("all")
visualize_scalar_on_gauss_points(Se, mmesh)
savefig("stress.png")


text:
module CloudFormation

using AWS
using AWS: AWSExceptions.AWSException
using EzXML
using MbedTLS: MbedException
using Memento
using Mocking
using OrderedCollections: OrderedDict
using XMLDict

# Improper casing to avoid issues with Module name and AWS.AWSService
@service cloudFormation

const logger = getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(MyModule)`
# NOTE: If this line is not included then the precompiled `MyModule.logger` won't be registered at runtime.
__init__() = Memento.register(logger)

export raw_stack_description, stack_output

function describe_stacks(config::AWSConfig, params::AbstractDict)
    dep_msg = """
    `describe_stacks(config::AWSConfig, params::AbstractDict)` is deprecated and will be removed.
    Use the AWS @service CloudFormation.describe_stacks() functionality instead.
    """
Base.depwarn(dep_msg, :get_authorization_token)
    return cloudFormation.describe_stacks(params; aws_config=config)
end

const _NRETRIES = 5

# for APIs we don't want to hammer
function cautious_delays(; kwargs...)
    return ExponentialBackOff(; n=_NRETRIES, first_delay=5, max_delay=300, kwargs...)
end

function minimal_delays(; kwargs...)
    return ExponentialBackOff(; n=_NRETRIES, first_delay=0.1, max_delay=60, kwargs...)
end

"""
    raw_stack_description(stack_name::AbstractString) -> String

Returns the description for the specified stack. Can optionally pass in the aws `config`
as a keyword argument.
"""
function raw_stack_description(
    stack_name::AbstractString; config::AWSConfig=global_aws_config()
)
    function retry_cond(s, e)
        if e isa AWSException
            if 500 <= e.cause.status <= 504
                debug(logger, "CloudFormation request encountered $(e.code); retrying")
                return (s, true)
            elseif e.cause.status == 429 || (e.cause.status == 400 && e.code == "Throttling")
                debug(logger, "CloudFormation request encountered $(e.code); retrying")
                return (s, true)
            end
        elseif e isa MbedException
            debug(logger, "CloudFormation request encountered $e; retrying")
            return (s, true)
        end

        return (s, false)
    end

    f = retry(; delays=cautious_delays(; jitter=0.2), check=retry_cond) do
        @mock describe_stacks(config, Dict("StackName" => stack_name, "return_raw" => true))
    end

    response = String(f())

    return response
end

"""
    stack_output(stack_name::AbstractString) -> OrderedDict

The stack's OutputKey and OutputValue values as a dictionary. Can pass in the aws `config`
as a keyword argument.
"""
function stack_output(stack_name::AbstractString; config::AWSConfig=global_aws_config())
    outputs = OrderedDict{String,String}()
    description = raw_stack_description(stack_name; config=config)

    xml = root(parsexml(description))
    ns = ["ns" => namespace(xml)]
    output_elements = findall("//ns:Outputs/ns:member", xml, ns)

    for el in output_elements
        key = nodecontent(findfirst("ns:OutputKey", el, ns))
        val = nodecontent(findfirst("ns:OutputValue", el, ns))
        outputs[key] = val
    end

    return outputs
end

function output_pair(item::AbstractDict)
    key = item["OutputKey"]::String
    value = if isa(item["OutputValue"], String)
        item["OutputValue"]::String
    elseif isa(item["OutputValue"], AbstractDict) && isempty(item["OutputValue"])
        ""
    else
        throw(ArgumentError("Unhandled output value: $(repr(item["OutputValue"]))"))
    end

    return key => value
end

end  # CloudFormation


text:
function model()
    T = Float32
    x = Var()
    y = Embedding(T,10000,100)(x)
end

function getscore(m, s::State)
    fs = features(s)
    score = 0.0

end


text:
"""
Layouts are used to organize the `canvas` and other widgets
within the `window`.

They are simple function or callable structures.

Example of layout implementation:

```julia
struct MyLayout <: AbstractLayout
    fields...
    function MyLayout(args...)
        ...
    end
end

function (layout::MyLayout)(window, canvas, widgets...)
    ...
end
```

The `call` implementation of a layout should return a widget
that contains (or not) the canvas and the widgets.
The returned widget will then be added to the window.
"""
abstract type AbstractLayout <: Function end

"""
The `canvas` is the only widget and fills the entire window.
"""
function CanvasOnly(::Widget, canvas::Canvas, widgets...)
    for widget in widgets
        if widget isa AbstractWidget{T} where {T}
            @warn "Attempt to add widget of type $T, but the current layout is CanvasOnly."
        end
    end
    return canvas
end

"""
The `SideBar` layout divides the window between the `canvas`
and a `sidebar` containing the other widgets.
"""
struct SideBar <: AbstractLayout
    position::Symbol
    size::Real
    props::NamedTuple
end
"""
        SideBar(position = :left, size = .25, props = NamedTuple()) -> SideBar

Creates a [`SideBar`](@ref) layout. The `position` determines where the sidebar is placed relative
to the canvas (can be either `:left` or `:right`). The `size` controls the width of the
sidebar relative to the width of the window (must be a value between 0 and 1). You can
also pass aditional properties for the sidebar widget using the `props` parameter.
"""
function SideBar(position = :right; size = .25, props = NamedTuple())
    @assert position in (:left, :right) "Invalid position '$position'. Valid positions are :left or :right"
    @assert size > 0 "Size cannot be negative"
    return SideBar(position, size, props)
end

function (this::SideBar)(window::Widget, canvas::Canvas, widgets...)
    props = this.props
    if this.position == :left
        pos = floor(Int, width(window) * this.size)
        return Paned(:h, position = pos,
            Box(:v, widgets...; props...),
            canvas,
        ) # Paned
    else
        pos = floor(Int, width(window) - width(window) * this.size)
        return Paned(:h, position = pos,
            canvas,
            Box(:v, widgets...; props...),
        ) # Paned
    end
end


text:
"""
Module with general utilities and missing standard library features
that could be useful in any Julia project
"""
module Utils

include("misc.jl")
include("cubitvector.jl")
include("data.jl")

end #module

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
@doc raw"""
    postproc(model::XY, param::Parameter, obs::MCObservableSet)

# Observables to be calculated
In the following, $m$ is total magnetization per site and $\epsilon$ is total energy per site.

- `"Binder Ratio x"`
    - $\frac{\left \langle m_x^4 \right \rangle}{\left \langle m_x^2 \right\rangle^2}$
- `"Binder Ratio y"`
    - $\frac{\left \langle m_y^4 \right \rangle}{\left \langle m_y^2 \right\rangle^2}$
- `"Binder Ratio"`
    - $\frac{\left \langle |m|^4 \right \rangle}{\left \langle |m|^2 \right\rangle^2}$
- `"Susceptibility x"`
    - $\frac{N}{T}\left(\left\langle m_x^2\right\rangle\right)$
- `"Susceptibility y"`
    - $\frac{N}{T}\left(\left\langle m_y^2\right\rangle\right)$
- `"Susceptibility y"`
    - $\frac{N}{T}\left(\left\langle |m|^2\right\rangle\right)$
- `"Connected Susceptibility x"`
    - $\frac{N}{T}\left(\left\langle m_x^2\right\rangle - \left\langle |m_x| \right\rangle^2\right)$
- `"Connected Susceptibility y"`
    - $\frac{N}{T}\left(\left\langle m_y^2\right\rangle - \left\langle |m_y| \right\rangle^2\right)$
- `"Connected Susceptibility"`
    - $\frac{N}{T}\left(\left\langle |m|^2\right\rangle - \left\langle |m| \right\rangle^2\right)$
- `"Specific Heat"`
    - $\frac{N}{T^2}\left(\left\langle \epsilon^2\right\rangle - \left\langle \epsilon \right\rangle^2\right)$
"""
function postproc(model::XY, param::Parameter, obs::MCObservableSet)
    nsites = numsites(model)
    T = convert(Float64, param["T"])
    beta = 1.0/T

    jk = jackknife(obs)
    jk["Binder Ratio x"] = jk["Magnetization x^4"] / (jk["Magnetization x^2"]^2)
    jk["Binder Ratio y"] = jk["Magnetization y^4"] / (jk["Magnetization y^2"]^2)
    jk["Binder Ratio"] = jk["|Magnetization|^4"] / (jk["|Magnetization|^2"]^2)
    jk["Susceptibility x"] = (nsites*beta)*jk["Magnetization x^2"]
    jk["Susceptibility y"] = (nsites*beta)*jk["Magnetization y^2"]
    jk["Susceptibility"] = (nsites*beta)*jk["|Magnetization|^2"]
    jk["Connected Susceptibility x"] = (nsites*beta)*(jk["Magnetization x^2"] - jk["|Magnetization x|"]^2)
    jk["Connected Susceptibility y"] = (nsites*beta)*(jk["Magnetization y^2"] - jk["|Magnetization y|"]^2)
    jk["Connected Susceptibility"] = (nsites*beta)*(jk["|Magnetization|^2"] - jk["|Magnetization|"]^2)
    jk["Specific Heat"] = (nsites*beta*beta)*(jk["Energy^2"] - jk["Energy"]^2)
    return jk
end