text:
"""
    grdimage(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; kwargs...)

Produces a gray-shaded (or colored) map by plotting rectangles centered on each grid node and assigning them a gray-shade (or color) based on the z-value.

Full option list at [`grdimage`]($(GMTdoc)grdimage.html)

Parameters
----------

- **A** | **img_out** | **image_out** :: [Type => Str]

    Save an image in a raster format instead of PostScript.
    ($(GMTdoc)grdimage.html#a)
- $(GMT.opt_J)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** | **img_in** | **image_in** :: [Type => Str]

    Specifies that the grid supplied is an image file to be read via GDAL.
    ($(GMTdoc)grdimage.html#d)
- **E** | **dpi** :: [Type => Int]

    Sets the resolution of the projected grid that will be created.
    ($(GMTdoc)grdimage.html#e)
- **G** | **bit_color** :: [Type => Int]

    ($(GMTdoc)grdimage.html#g)
- **I** | **shade** | **intensity** :: [Type => Bool | Str | GMTgrid]

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    ($(GMTdoc)grdimage.html#i)
- **M** | **monochrome** :: [Type => Bool]

    Force conversion to monochrome image using the (television) YIQ transformation.
    ($(GMTdoc)grdimage.html#m)
- **N** | **noclip** :: [Type => Bool]

    Do not clip the image at the map boundary.
    ($(GMTdoc)grdimage.html#n)
- $(GMT.opt_P)
- **Q** | **nan_t** | **nan_alpha** :: [Type => Bool]

    Make grid nodes with z = NaN transparent, using the colormasking feature in PostScript Level 3.
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdimage(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdimage", cmd0, arg1, arg2, arg3)
	arg4 = nothing		# For the r,g,b + intensity case

	d = KW(kwargs)
    K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd = parse_common_opts(d, cmd, [:UVXY :params :c :f :n :p :t], first)
	cmd = parse_these_opts(cmd, d, [[:A :img_out :image_out], [:D :img_in :image_in], [:E :dpi], [:G :bit_color],
				[:M :monochrome], [:N :noclip], [:Q :nan_t :nan_alpha], ["," :mem :mem_layout]])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)		# Find how data was transmitted
	if (got_fname == 0 && isa(arg1, Tuple))			# Then it must be using the three r,g,b grids
		cmd, got_fname, arg1, arg2, arg3 = find_data(d, cmd0, cmd, arg1, arg2, arg3)
	end

	if (isa(arg1, Array{<:Number}))
		arg1 = mat2grid(arg1)
		if (isa(arg2, Array{<:Number}))  arg2 = mat2grid(arg2)  end
		if (isa(arg3, Array{<:Number}))  arg3 = mat2grid(arg3)  end
	end

	cmd, N_used, arg1, arg2, arg3 = get_cpt_set_R(d, cmd0, cmd, opt_R, got_fname, arg1, arg2, arg3, "grdimage")
	cmd, arg1, arg2, arg3, arg4 = common_shade(d, cmd, arg1, arg2, arg3, arg4, "grdimage")

	if (isa(arg1, GMTimage) && !occursin("-D", cmd))  cmd *= " -D"  end	# GMT bug. It says not need but it is.
	cmd = "grdimage " * cmd				# In any case we need this
	do_finish = false
	if (!occursin("-A", cmd))			# -A means that we are requesting the image directly
		cmd, K = finish_PS_nested(d, cmd, "", K, O, [:coast :colorbar :basemap])
		do_finish = true
	end
	return finish_PS_module(d, cmd, "", K, O, do_finish, arg1, arg2, arg3, arg4)
end

# ---------------------------------------------------------------------------------------------------
function common_shade(d::Dict, cmd::String, arg1, arg2, arg3, arg4, prog)
	# Used both by grdimage and grdview
	if ((val = find_in_dict(d, [:I :shade :intensity])[1]) !== nothing)
		if (!isa(val, GMTgrid))			# Uff, simple. Either a file name or a -A type modifier
			if (isa(val, String) || isa(val, Symbol) || isa(val, Bool))
				val = arg2str(val)
				(val == "" || val == "default" || val == "auto") ? cmd *= " -I+a-45+nt1" : cmd *= " -I" * val
			else
				cmd = add_opt(cmd, 'I', d, [:I :shade :intensity],
				              (auto="_+", azim="+a", azimuth="+a", norm="+n", default="_+d+a-45+nt1"))
			end
		else
			if (prog == "grdimage")  cmd, N_used = put_in_slot(cmd, val, 'I', [arg1, arg2, arg3, arg4])
			else                     cmd, N_used = put_in_slot(cmd, val, 'I', [arg1, arg2, arg3])
			end
			if     (N_used == 1)  arg1 = val
			elseif (N_used == 2)  arg2 = val
			elseif (N_used == 3)  arg3 = val
			elseif (N_used == 4)  arg4 = val	# grdview doesn't have this case but no harm to not test for that
			end
		end
	end
	return cmd, arg1, arg2, arg3, arg4
end

# ---------------------------------------------------------------------------------------------------
grdimage!(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; first=false, kw...) =
	grdimage(cmd0, arg1, arg2, arg3; first=false, kw...)

grdimage(arg1,  arg2=nothing, arg3=nothing; first=true, kw...)  = grdimage("", arg1, arg2, arg3; first=first, kw...)
grdimage!(arg1, arg2=nothing, arg3=nothing; first=false, kw...) = grdimage("", arg1, arg2, arg3; first=first, kw...)


text:
help_str(x; color = true) = sprint(print_cmd, x; context = :color => color)
# printing in expression
emit_help(x, ptr::Int=1; color = true) = quote
    if !isnothing(findnext(isequal("-h"), ARGS, $ptr)) ||
        !isnothing(findnext(isequal("--help"), ARGS, $ptr))
        print($(help_str(x; color = color)))
        return 0
    end
end

function emit_error(cmd, msg::String; color::Bool=true)
    msg = "Error: $msg, use -h or --help to check more detailed help info"
    return quote
        printstyled($msg; color = :red, bold = true)
        println()
        print($(help_str(cmd; color = color)))
        println()
        return 1
    end
end

function emit_error(cmd, msg::Expr; color::Bool=true)
    msg.head === :string || throw(Meta.ParseError("expect string expression, got $msg"))

    ex = Expr(:string, "Error: ")
    for each in msg.args
        if each isa Number
            push!(ex.args, string(each))
        else
            push!(ex.args, each)
        end
    end
    push!(ex.args, ", use -h or --help to check more detailed help info")
    return quote
        printstyled($msg; color = :red, bold = true)
        println()
        print($(help_str(cmd; color = color)))
        println()
        return 1
    end
end

"""
    emit_expr(cmd::Entry[, ptr::Int=1])

Emit `Expr` from a `Entry`.
"""
function emit_expr(cmd::Entry, ptr::Int=1)
    jlfn = JLFunction(;
        name=:command_main,
        args=[Expr(:kw, :(ARGS::Vector{String}), :ARGS)],
        body=quote
            $(emit_scan_version(cmd))
            $(emit_body(cmd.root, ptr))
        end
    )
    return codegen_ast(jlfn)
end

function emit_scan_version(cmd::Entry)
    quote
        if "-V" in ARGS || "--version" in ARGS
            print($(cmd.version))
            return 0
        end
    end
end

function emit_body(cmd::NodeCommand, ptr::Int=1)
    nargs_assert = quote
        if length(ARGS) < $ptr
            $(emit_error(cmd, "valid sub-commands for command $(cmd.name) are: $(join(keys(cmd.subcmds), ", "))"))
        end
    end

    jl = JLIfElse()
    for (name, subcmd) in cmd.subcmds
        jl[:(ARGS[$ptr] == $name)] = emit_body(subcmd, ptr+1)
    end
    jl.otherwise = emit_error(cmd, :("Error: unknown command $(ARGS[$ptr])"))
    return quote
        if length(ARGS) == $ptr && (ARGS[$(ptr)] == "-h" || ARGS[$(ptr)] == "--help")
            print($(help_str(cmd)))
            return 0
        end

        $nargs_assert
        $(codegen_ast(jl))
    end
end

function emit_body(cmd::LeafCommand, ptr::Int=1)
    @gensym idx
    quote
        $(emit_help(cmd, ptr))

        $idx = findnext(isequal("--"), ARGS, $ptr)
        if isnothing($idx) # no dash
            $(emit_norm_body(cmd, ptr))
        else # dash
            $(emit_dash_body(cmd, idx, ptr))
        end
    end
end

function emit_dash_body(cmd::LeafCommand, idx::Symbol, ptr::Int=1)
    @gensym token token_ptr args kwargs

    quote # parse option/flag
        $kwargs = []
        $token_ptr = $ptr
        while $token_ptr ≤ $idx-1
            $token = ARGS[$token_ptr]
            if startswith($token, "-")
                $(emit_kwarg(cmd, token, kwargs, token_ptr))
            else
                $(emit_error(cmd, :("unknown command: $($token)")))
            end
            $token_ptr += 1
        end

        $args = ARGS[$idx+1:end]
        $(emit_leaf_call(cmd, args, kwargs))
    end
end

function emit_norm_body(cmd::LeafCommand, ptr::Int=1)
    @gensym token_ptr args kwargs token

    quote
        $args = []; $kwargs = [];
        sizehint!($args, $(cmd.nrequire))
        $token_ptr = $ptr
        while $token_ptr ≤ length(ARGS)
            $token = ARGS[$token_ptr]
            if startswith($token, "-")
                $(emit_kwarg(cmd, token, kwargs, token_ptr))
            else # argument
                push!($args, $token)
            end
            $token_ptr += 1
        end

        $(emit_leaf_call(cmd, args, kwargs))
    end
end

function emit_leaf_call(cmd::LeafCommand, args::Symbol, kwargs::Symbol)
    @gensym nargs
    ret = quote
        $nargs = length($args)
    end

    if cmd.nrequire > 0
        push!(ret.args, quote
            if $(cmd.nrequire) > $nargs
                $(emit_error(cmd, "expect $(cmd.nrequire) positional arguments"))
            end
        end)
    end

    # check maximum number of arguments
    if isnothing(cmd.vararg)
        push!(ret.args, quote
            if $(length(cmd.args)) < $nargs
                $(emit_error(cmd, "expect at most $(length(cmd.args)) positional arguments"))
            end
        end)
    end

    call = Expr(:call, cmd.fn)
    if !isempty(cmd.flags) || !isempty(cmd.options)
        push!(call.args, Expr(:parameters, :($kwargs...)))
    end

    for (i, arg) in enumerate(cmd.args)
        arg.require || break
        push!(call.args, emit_parse_value(cmd, arg.type, :($args[$i])))
    end

    ifelse = JLIfElse()
    ifelse[:($nargs == $(cmd.nrequire))] = quote
        $call
        return 0
    end

    for i in cmd.nrequire+1:length(cmd.args)
        call = copy(call)
        type = cmd.args[i].type
        push!(call.args, emit_parse_value(cmd, type, :($args[$i])))
        ifelse[:($nargs == $i)] = quote
            $call
            return 0
        end
    end

    if isnothing(cmd.vararg)
        ifelse.otherwise = emit_error(cmd, "expect at most $(length(cmd.args)) positional arguments")
    else
        @gensym varargs
        call = copy(call)
        push!(call.args, :($varargs...))
        type = cmd.vararg.type
        if type === Any || type === String || type === AbstractString
            ifelse.otherwise = quote
                $varargs = $args[$(length(cmd.args))+1:end]
                $call
                return 0
            end
        else
            ifelse.otherwise = quote
                $varargs = map($args[$(length(cmd.args)+1):end]) do value
                    $(emit_parse_value(cmd, type, :value))
                end
                $call
                return 0
            end
        end
    end

    push!(ret.args, codegen_ast(ifelse))
    return ret
end

function emit_kwarg(cmd::LeafCommand, token::Symbol, kwargs::Symbol, token_ptr)
    if isempty(cmd.flags) && isempty(cmd.options)
        return emit_error(cmd, :("do not have $($token)"))
    end

    @gensym sym key value

    ifelse = JLIfElse()
    # short flag
    ifelse[:(length($token) == 2)] = quote
        $key = $token[2:2]
        $(emit_short_flag(cmd, token, sym, key, value))
    end

    # long option/flag
    ifelse[:(startswith($token, "--"))] = quote
        $key = lstrip(split($token, '=')[1], '-')
        $(emit_long_option_or_flag(cmd, token, sym, key, value, token_ptr))
    end

    # short option
    ifelse.otherwise = quote
        $key = $token[2:2]
        $(emit_short_option(cmd, token, sym, key, value, token_ptr))
    end

    return quote
        $(codegen_ast(ifelse))
        push!($kwargs, $sym=>$value)
    end
end

function emit_short_flag(cmd::LeafCommand, token::Symbol, sym::Symbol, key::Symbol, value::Symbol)
    ifelse = JLIfElse()
    for (_, flag) in cmd.flags
        if flag.short
            name = flag.name[1:1]
            ifelse[:($key == $name)] = quote
                $sym = $(QuoteNode(flag.sym))
                $value = true
            end
        end
    end
    ifelse.otherwise = emit_error(cmd, :("cannot find flag: $($token)"))
    return codegen_ast(ifelse)
end

function emit_short_option(cmd::LeafCommand, token::Symbol, sym::Symbol, key::Symbol, value::Symbol, token_ptr)
    ifelse = JLIfElse()
    for (_, option) in cmd.options
        if option.short
            name = option.name[1:1]
            ifelse[:($key == $name)] = quote
                if occursin('=', $token)
                    _, $value = split($token, '=')
                elseif length($token) == 2 # read next
                    $token_ptr += 1
                    if $token_ptr > length(ARGS)
                        $(emit_error(option, "expect a value"))
                    end
                    $value = ARGS[$token_ptr]
                else # -o<value>
                    $value = $token[3:end]
                end
                $sym = $(QuoteNode(option.sym))
                $value = $(emit_parse_value(option, option.type, value))
            end
        end
    end
    ifelse.otherwise = emit_error(cmd, :("cannot find $($token)"))
    return codegen_ast(ifelse)
end

function emit_long_option_or_flag(cmd::LeafCommand, token::Symbol, sym::Symbol, key::Symbol, value::Symbol, token_ptr)
    ifelse = JLIfElse()
    for (name, flag) in cmd.flags
        ifelse[:($key == $name)] = quote
            $sym = $(QuoteNode(flag.sym))
            $value = true
        end
    end

    for (name, option) in cmd.options
        ifelse[:($key == $name)] = quote
            $sym = $(QuoteNode(option.sym))
            $(emit_option(option, token, value, token_ptr))
        end
    end

    ifelse.otherwise = emit_error(cmd, :("cannot find $($token)"))
    return codegen_ast(ifelse)
end

function emit_parse_value(cmd, type, value)
    if type === Any || type === String || type === AbstractString
        return value
    else
        @gensym ret
        return quote
            $ret = tryparse($type, $value)
            if isnothing($ret)
                $(emit_error(cmd, "expect value of type: $(type)"))
            end
            $ret
        end
    end
end

function emit_option(option::Option, token::Symbol, value::Symbol, token_ptr::Symbol)
    return quote
        if occursin('=', $token)
            _, $value = split($token, '=')
        else # read next token
            $token_ptr += 1
            if $token_ptr > length(ARGS)
                $(emit_error(option, "expect a value"))
            end
            $value = ARGS[$token_ptr]
        end
        $value = $(emit_parse_value(option, option.type, value))
    end
end

ext: jl


text:
struct CPAParam <: EoSParam
    a::PairParam{Float64}
    b::PairParam{Float64}
    c1::SingleParam{Float64}
    Tc::SingleParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type CPAModel <: EoSModel end

struct CPA{T <: IdealModel,c <: CubicModel} <: CPAModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    cubicmodel::c
    params::CPAParam
    sites::SiteParam
    idealmodel::T
    assoc_options::AssocOptions
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel CPA
export CPA
function CPA(components;
    idealmodel=BasicIdeal,
    cubicmodel=RK,
    alpha=CPAAlpha,
    mixing=vdW1fRule,
    activity=nothing,
    translation=NoTranslation,
    userlocations=String[],
    ideal_userlocations=String[],
    alpha_userlocations=String[],
    activity_userlocations=String[],
    mixing_userlocations=String[],
    translation_userlocations=String[],
    verbose=false,
    assoc_options = AssocOptions())

    icomponents = 1:length(components)

    params,sites = getparams(components, ["SAFT/CPA", "properties/molarmass.csv","properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    Mw  = params["Mw"]
    k  = params["k"]
    Tc = params["Tc"]
    c1 = params["c1"]
    params["a"].values .*= 1E-1
    params["b"].values .*= 1E-3
    a  = epsilon_LorentzBerthelot(params["a"], k)
    b  = sigma_LorentzBerthelot(params["b"])
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]
    packagedparams = CPAParam(a, b, c1, Tc, epsilon_assoc, bondvol,Mw)

    init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
    init_alpha = init_model(alpha,components,alpha_userlocations,verbose)
    init_mixing = init_model(mixing,components,activity,mixing_userlocations,activity_userlocations,verbose)
    init_translation = init_model(translation,components,translation_userlocations,verbose)

    if occursin("RK",string(cubicmodel))
        cubicparams = RKParam(a, b, params["Tc"],params["pc"],Mw)
    elseif occursin("PR",string(cubicmodel))
        cubicparams = PRParam(a, b, params["Tc"],params["pc"],Mw)
    end

    init_cubicmodel = cubicmodel(components,icomponents,init_alpha,init_mixing,init_translation,cubicparams,init_idealmodel,1e-12,String[])

    references = ["10.1021/ie051305v"]

    model = CPA(components, icomponents, init_cubicmodel, packagedparams, sites, init_idealmodel, assoc_options, 1e-12, references)
    return model
end

lb_volume(model::CPAModel,z = SA[1.0]) = lb_volume(model.cubicmodel,z)
T_scale(model::CPAModel,z=SA[1.0]) = T_scale(model.cubicmodel,z)
p_scale(model::CPAModel,z=SA[1.0]) = p_scale(model.cubicmodel,z)

function x0_crit_pure(model::CPAModel)
    lb_v = lb_volume(model)
    if isempty(model.params.epsilon_assoc.values[1,1])
        [2.0, log10(lb_v/0.3)]
    else
        [2.5, log10(lb_v/0.3)]
    end
end

function a_res(model::CPAModel, V, T, z)
    n = sum(z)
    ā,b̄,c̄ = cubic_ab(model.cubicmodel,V,T,z,n)
    return a_res(model.cubicmodel,V,T,z) + a_assoc(model,V+c̄*n,T,z)
end

ab_consts(model::CPAModel) = ab_consts(model.cubicmodel)

function Δ(model::CPAModel, V, T, z, i, j, a, b)
    ϵ_associjab = model.params.epsilon_assoc.values[i,j][a,b] * 1e2/R̄
    βijab = model.params.bondvol.values[i,j][a,b] * 1e-3
    Σz = sum(z)
    b = model.params.b.values
    b̄ = dot(z,b,z)
    η = b̄/(4*V*Σz)
    g = (1-η/2)/(1-η)^3
    bij = (b[i,i]+b[j,j])/2
    return g*(exp(ϵ_associjab/T)-1)*βijab*bij/N_A
end


text:
using DataStructures
using JLD2

"""
    $(SIGNATURES)

Saves `object` into `.jld2` file `filename` (the suffix is added automatically if not entered).
The resulting file contains a dictionary with a single entry.

"""
function save(filename, object)
    JLD2.save(_jld2_name(filename), Dict("object" => object))
end

function save(filename, x::Result)
    x_nofunc = deepcopy(x)

    # compiled functions cause problems in saving: ignore J now, compile when loading
    null_J() = 0
    x_nofunc.jacobian = null_J
    JLD2.save(_jld2_name(filename), Dict("object" => x_nofunc))
end

_jld2_name(filename) = filename[end-4:end] == ".jld2" ? filename : filename * ".jld2"

"""
    $(SIGNATURES)

Loads an object from `filename`. For objects containing symbolic expressions such as `HarmonicEquation`, the symbolic variables are
reinstated in the `HarmonicBalance` namespace.

"""
function load(filename)
    loaded = JLD2.load(filename)

    if haskey(loaded,"object") #otherwise save data is from a plot
        loaded = loaded["object"]

        # we need the symbols in our namespace to parse strings with `transform_solutions`
        _parse_symbol_names(loaded)

        # automatic saving fails for some objects: reconstruct these manually
        _parse_loaded(loaded)
    else
        return loaded
    end
end

function _parse_loaded(x::Problem)
    # reconstruct the HomotopyContinuation System
    system = HC_wrapper.System(x.eom)
    x.system = system
    return x
end

function _parse_loaded(x::Result)
    # reconstruct System and the compiled Jacobian
    x.problem.system = HC_wrapper.System(x.problem.eom)
    x.jacobian = _compile_Jacobian(x.problem, x.swept_parameters, x.fixed_parameters)
    return x
end

_parse_loaded(x) = x

"Retrieve names for all symbols and declare them in this namespace."
function _parse_symbol_names(x::Problem)
    all_symbols= cat(x.parameters, get_variables(x.variables),
    get_independent_variables(x.eom), get_independent_variables(x.eom.natural_equation), dims=1)
    return declare_variable.(string.(all_symbols))
end

_parse_symbol_names(x::Result) = _parse_symbol_names(x.problem)
_parse_symbol_names(x) = nothing


text:
@testset "Refine triangulation with images" begin
    @testset "E = $E" for E in 3:3
        @testset "k = $k" for k in 3:3
            n_points = 10
            tau = 1
            #println(E, " ", k)
            splitting_rules = simplicial_subdivision(k, E)
            embedding = embedding_example(n_points+E*2, E, tau)

            # Triangulate all but the last point
            triang_points, simplex_indices = triangulate(embedding[1:size(embedding, 1), :])

            # Take the images of the simplices in the triangulation as the forward projections
            # in time of the original simplices.
            image_indices = simplex_indices + 1
            image_points = vcat(triang_points[2:end, :], embedding[1, :].')
            #@show E, k, triang_points, simplex_indices, image_points

            # How many simplices did we get?
            n_orig = size(simplex_indices, 1)

            simplices = [triang_points[simplex_indices[i, :], :] for i in 1:n_orig]

            split_indices = collect(1:size(embedding, 2))

            refined_triang = refine_triangulation_images(
                triang_points, simplex_indices, image_points, split_indices, k)

            # How many new simplices did we get?
            n_refined = size(refined_triang[2], 1)
            refined_simplices = [refined_triang[1][refined_triang[2][i, :], :]
                                    for i in 1:n_refined]

            refined_vertices = refined_triang[1]
            refined_simplex_indices = refined_triang[2]

            @testset "Arrays are of right dimensions and types" begin
                @test isa(refined_vertices, Array{Float64, 2})
                @test isa(refined_simplex_indices, Array{Int, 2})
            end
            # Compare volume of original triangulation with the refined triangulation

            @testset "Volume is preserved" begin
                triang_vols = [abs(det(hcat(simplices[i], ones(E + 1)))) for i in 1:n_orig]
                refined_triang_vols = [abs(det(hcat(refined_simplices[i], ones(E + 1)))) for i in 1:n_refined]
                @test sum(triang_vols) ≈ sum(refined_triang_vols) atol = 1/10^7
            end
        end
    end
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
using Dictionaries
include("loggingprofiler.jl")

function retrieve_code_info(sigtypes, world = Base.get_world_counter())
    S = Tuple{map(s -> Core.Compiler.has_free_typevars(s) ? typeof(s.parameters[1]) : s, sigtypes)...}
    _methods = Base._methods_by_ftype(S, -1, world)
    if isempty(_methods)
        @info("method $(sigtypes) does not exist")
        return(nothing)
    end
    type_signature, raw_static_params, method = _methods[1]
    mi = Core.Compiler.specialize_method(method, type_signature, raw_static_params, false)
    ci = Base.isgenerated(mi) ? Core.Compiler.get_staged(mi) : Base.uncompressed_ast(method)
    Base.Meta.partially_inline!(ci.code, [], method.sig, Any[raw_static_params...], 0, 0, :propagate)
    ci
end

function overdubbable(ex::Expr)
    ex.head != :call && return(false)
    length(ex.args) < 2 && return(false)
    return(overdubbable(ex.args[1]))
end
overdubbable(gr::Core.GlobalRef) = gr.name ∉ [:overdub, :record_start, :record_end, :promote, :convert, :tuple]
# overdubbable(gr::Symbol) =
overdubbable(ex) = false
timable(ex) = overdubbable(ex)

#
remap(ex::Expr, maps) = Expr(ex.head, remap(ex.args, maps)...)
remap(args::AbstractArray, maps) = map(a -> remap(a, maps), args)
remap(c::Core.GotoNode, maps) = Core.GotoNode(maps.goto[c.label])
remap(c::Core.GotoIfNot, maps) = Core.GotoIfNot(remap(c.cond, maps), maps.goto[c.dest])
remap(r::Core.ReturnNode, maps) = Core.ReturnNode(remap(r.val, maps))
remap(a::Core.SlotNumber, maps) = maps.slots[a.id]
remap(a::Core.SSAValue, maps) = Core.SSAValue(maps.ssa[a.id])
remap(a::Core.NewvarNode, maps) = Core.NewvarNode(maps.slots[a.slot.id])
remap(a::GlobalRef, maps) = a
remap(a::QuoteNode, maps) = a
remap(ex, maps) = ex

exportname(ex::GlobalRef) = QuoteNode(ex.name)
exportname(ex::Symbol) = QuoteNode(ex)
exportname(ex::Expr) = exportname(ex.args[1])
exportname(i::Int) = QuoteNode(Symbol("Int(",i,")"))

dummy() = return
function empty_codeinfo()
    new_ci = code_lowered(dummy, Tuple{})[1]
    empty!(new_ci.code)
    empty!(new_ci.slotnames)
    empty!(new_ci.linetable)
    empty!(new_ci.codelocs)
    new_ci
end

overdub(f::Core.IntrinsicFunction, args...) = f(args...)
overdub(f::Core.Builtin, args...) = f(args...)

@generated function overdub(f::F, args...) where {F}
    ci = retrieve_code_info((F, args...))
    if ci === nothing
        return(Expr(:call, :f, [:(args[$(i)]) for i in 1:length(args)]...))
    end

    new_ci = empty_codeinfo()
    new_ci.slotnames = vcat([Symbol("#self#"), :f, :args], ci.slotnames[2:end])
    new_ci.slotflags = vcat([0x00, 0x00, 0x00], ci.slotflags[2:end])
    foreach(s -> push!(new_ci.linetable, s), ci.linetable)

    maps = (
        ssa = Dict{Int, Int}(),
        slots = Dict{Int, Any}(),
        goto = Dict{Int,Int}(),
        )

    #we need to map indexes of slot-variables from ci to their new values.
    # except the first one, we just remap them
    maps.slots[1] = Core.SlotNumber(1)
    foreach(i -> maps.slots[i] = Core.SlotNumber(i + 2), 2:length(ci.slotnames)) # they are shifted by 2 accomondating inserted `f` and `args`

    #if somewhere the original parameters of the functions will be used
    #they needs to be remapped to an SSAValue from here, since the overdubbed
    # function has signatures overdub(f, args...) instead of f(x,y,z...)
    newci_no = 0
    for i in 1:length(args)
        newci_no +=1
        push!(new_ci.code, Expr(:call, Base.getindex, Core.SlotNumber(3), i))
        maps.slots[i+1] = Core.SSAValue(newci_no)
        push!(new_ci.codelocs, ci.codelocs[1])
    end

    for (ci_no, ex) in enumerate(ci.code)
        if timable(ex)
            fname = exportname(ex)
            push!(new_ci.code, Expr(:call, GlobalRef(LoggingProfiler, :record_start), fname))
            push!(new_ci.codelocs, ci.codelocs[ci_no])
            newci_no += 1
            maps.goto[ci_no] = newci_no
            # if overdubbable(ex)
            #     ex = Expr(:call, GlobalRef(Main, :overdub), ex.args...)
            # end
            push!(new_ci.code, ex)
            push!(new_ci.codelocs, ci.codelocs[ci_no])
            newci_no += 1
            maps.ssa[ci_no] = newci_no
            push!(new_ci.code, Expr(:call, GlobalRef(LoggingProfiler, :record_end), fname))
            push!(new_ci.codelocs, ci.codelocs[ci_no])
            newci_no += 1
        else
            push!(new_ci.code, ex)
            push!(new_ci.codelocs, ci.codelocs[ci_no])
            newci_no += 1
            maps.ssa[ci_no] = newci_no
        end
    end

    for i in length(args)+1:length(new_ci.code)
        new_ci.code[i] = remap(new_ci.code[i], maps)
    end
    new_ci
    new_ci.inferred = false
    new_ci.ssavaluetypes = length(new_ci.code)
    # new_ci
    return(new_ci)
end

LoggingProfiler.reset!()
new_ci = overdub(sin, 1.0)
LoggingProfiler.to

function foo(x, y)
   z =  x * y
   z + sin(y)
end

LoggingProfiler.reset!()
overdub(foo, 1.0, 1.0)
LoggingProfiler.to

macro record(ex)
    Expr(:call, :overdub, ex.args...)
end

LoggingProfiler.reset!()
@record foo(1.0, 1.0)
LoggingProfiler.to


text:
function load!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    println("loading wind data")
    loadwind!(options,sizes,data,output)
    println("wind data loaded")

    println("loading turbine data")
    loadturbine!(options,sizes,data,output)
    println("turbine data loaded")

    println("initialize! data")
    initialize!(options,sizes,data,output)
    println("initialize! loaded")

    nothing
end

function loadwind!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    windPath = joinpath(options.PATHCASE,"wind.csv")
    windFile = CSV.File(windPath,header=0) |> Tables.matrix

    data.windnames = windFile[1,2:end]
    data.onshoreDistance = parse.(Float64,windFile[2,2:end])
    data.offshoreDistance = parse.(Float64,windFile[3,2:end])
    data.surveysCost = parse.(Float64,windFile[4,2:end])
    data.vesselTimeWT = parse.(Float64,windFile[5,2:end])
    data.vesselTimeSS = parse.(Float64,windFile[6,2:end])
    data.WD = parse.(Float64,windFile[7,2:end])
    data.Ntln = parse.(Int,windFile[8,2:end])
    data.Ctln = parse.(Float64,windFile[9,2:end])
    data.Ntlf = parse.(Float64,windFile[10,2:end])
    data.Ctlf = parse.(Float64,windFile[11,2:end])
    data.Cprot = parse.(Float64,windFile[12,2:end])
    data.Cport = parse.(Float64,windFile[13,2:end])
    data.Lr = parse.(Float64,windFile[14,2:end])
    data.Cins = parse.(Float64,windFile[15,2:end])

    data.l = parse.(Float64,windFile[16,2:end])
    data.Pe = parse.(Float64,windFile[17,2:end])
    data.ComIns = parse.(Float64,windFile[18,2:end])
    data.Ctransunit = parse.(Float64,windFile[19,2:end])
    data.d = parse.(Float64,windFile[20,2:end])
    data.tc = parse.(Float64,windFile[21,2:end])
    data.Cindport = parse.(Float64,windFile[22,2:end])
    data.Cindves = parse.(Float64,windFile[23,2:end])
    data.Cindlab = parse.(Float64,windFile[24,2:end])

    data.Cddport = parse.(Float64,windFile[25,2:end])
    data.Nddport = parse.(Float64,windFile[26,2:end])
    data.Lddport = parse.(Float64,windFile[27,2:end])
    data.Nremove = parse.(Float64,windFile[28,2:end])
    data.Vremove = parse.(Float64,windFile[29,2:end])
    data.Asc = parse.(Float64,windFile[30,2:end])
    data.Cscunit = parse.(Float64,windFile[31,2:end])

    data.windspeed = parse.(Float64,windFile[32:end,1])
    data.RFwindspeed = parse.(Float64,windFile[32:end,2:end])
    sizes.winds = length(data.windnames)

    nothing
end

function loadturbine!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    turbinePath = joinpath(options.PATHCASE,"turbines.csv")
    turbineFile = CSV.File(turbinePath,header=0) |> Tables.matrix

    data.turbinenames = turbineFile[1,2:end]
    data.turbineQuantity = parse.(Int,turbineFile[2,2:end])
    data.turbinePR = parse.(Float64,turbineFile[3,2:end])
    data.height = parse.(Float64,turbineFile[4,2:end])
    data.diameter = parse.(Float64,turbineFile[5,2:end])

    data.N1d = parse.(Float64,turbineFile[6,2:end])
    data.Ccomp = parse.(Float64,turbineFile[7,2:end])
    data.Ncom = parse.(Float64,turbineFile[8,2:end])
    data.Lcom = parse.(Float64,turbineFile[9,2:end])

    data.N1dOM = parse.(Float64,turbineFile[10,2:end])
    data.LrOM = parse.(Float64,turbineFile[11,2:end])
    data.Pd = parse.(Float64,turbineFile[12,2:end])
    data.lambda = parse.(Float64,turbineFile[13,2:end])
    data.Csm = parse.(Float64,turbineFile[14,2:end])

    data.Wjproc = parse.(Float64,turbineFile[15,2:end])
    data.Cjproc = parse.(Float64,turbineFile[16,2:end])
    data.Wjtrans = parse.(Float64,turbineFile[17,2:end])
    data.Wttrans = parse.(Float64,turbineFile[18,2:end])
    data.Cjtrnas = parse.(Float64,turbineFile[19,2:end])
    data.Wnr = parse.(Float64,turbineFile[20,2:end])
    data.Cnr = parse.(Float64,turbineFile[21,2:end])
    data.Wr = parse.(Float64,turbineFile[22,2:end])
    data.SV = parse.(Float64,turbineFile[23,2:end])

    windspeed = turbineFile[24:end,1]
    if parse.(Float64,windspeed) != data.windspeed
        error("wind speed is inconsistenty")
    end

    powercurve = turbineFile[24:end,2:end]
    data.powercurve = parse.(Float64,powercurve)
    sizes.turbines = length(data.turbinenames)

    nothing
end

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
function get_fluent_range_indices(graph, level, max_num_fluents)
    lvl = 5(level-1) #+ 1
    row = length(graph.props[level-1])
    relevant = (1:row, lvl)
    nonrelevant = (row+1:max_num_fluents, lvl+1:lvl+5)
    return relevant, nonrelevant
end

function get_fluent_terminal_indices(graph; loc=:init)
    indices = []
    not_ind = []
    T = graph.num_levels
    if loc == :init
        for i=1:length(graph.initprops)
            ind = findall(x->x == graph.initprops[i], graph.props[0])[1]
            push!(indices, ind)
        end
    else
        for i=1:length(graph.goalprops)
            ind = findall(x->x == graph.goalprops[i], graph.props[T])[1]
            push!(indices, ind)
        end
    end
    not_ind = [i for i=1:length(graph.props[0]) if !(i in indices)]
    return indices, not_ind
end

function get_action_indices(graph, fluent, level, set, except_set)
    # actions = graph.acts[level]
    actions = [act for act in graph.acts[level] if act.name != :NoOp]
    action_indices = []
    if set == :add && except_set == :pre
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.pos_eff && !(fluent in act.pos_prec)
                push!(action_indices, i)
            end
        end
    elseif set == :del && except_set == :pre
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.neg_eff && !(fluent in act.pos_prec)
                push!(action_indices, i)
            end
        end
    elseif set == :pre && except_set == :del
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.pos_prec && !(fluent in act.neg_eff)
                push!(action_indices, i)
            end
        end
    else
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.pos_prec && fluent in act.neg_eff
                push!(action_indices, i)
            end
        end
    end
    return action_indices
end

#=
state change variables xₜ are encoded as
            1     2      3    4     5
        preadd, predel, add, del, maintain
            pa,     pd,  ad,  dl,   mn
=#
function encode_optiplan(domain_name, problem_name; max_levels=10)
    domain = load_domain(domain_name)
    problem = load_problem(problem_name)
    graph = create_graph(domain, problem; max_levels=max_levels)

    model = Model(Gurobi.Optimizer)
    set_silent(model)
    T = graph.num_levels #+ 1
    max_nf = max([length(graph.props[i]) for i=1:T]...)
    max_na = max([length([act for act in graph.acts[i] if act.name != :NoOp]) for i=1:T]...)

    pa, pd, ad, dl, mn = 1,2,3,4,5
    #=
        VARIABLES
    =#
    @variable(model, x[i=1:max_nf, j=1:5*T], Bin)
    @variable(model, y[i=1:max_na, j=1:T], Bin)

    #=
        CONSTRAINTS
    =#

    #init constraints
    # Eq. 2
    init_idx, _ = get_fluent_range_indices(graph, 1, max_nf)
    inits, notinits = get_fluent_terminal_indices(graph; loc=:init)
    for i in inits
        @constraint(model, x[i, init_idx[2]+ad] == 1)
    end
    #Eq. 3
    for i in notinits
        @constraint(model, x[i, init_idx[2]+ad] == 0)
        @constraint(model, x[i, init_idx[2]+mn] == 0)
        @constraint(model, x[i, init_idx[2]+pa] == 0)
    end
    # @constraint(model, x[init_idx[1], init_idx[2]+mn] .== 0)
    # @constraint(model, x[init_idx[1], init_idx[2]+pa] .== 0)
    # @constraint(model, x[init_idx[1], init_idx[2]+pd] .== 0)
    # @constraint(model, x[init_idx[1], init_idx[2]+dl] .== 0)

    # @constraint(model, x[nr_init_idx[1], nr_init_idx[2]] .== 0)

    #goal constraints
    #Eq. 4
    goals, _ = get_fluent_terminal_indices(graph; loc=:goal)
    gidx, _ = get_fluent_range_indices(graph, T, max_nf)
    for g in goals
        @constraint(model, x[g,gidx[2]+ad] + x[g,gidx[2]+mn] + x[g, gidx[2]+pa] >= 1)
    end

    # one action at first time step
    num_actions = length([act for act in graph.acts[1] if act.name != :NoOp])
    @constraint(model, sum(y[1:num_actions, 1]) == 1)

    for t=2:T
        #effect implication constraints
        fidx,_ = get_fluent_range_indices(graph, t, max_nf)
        # Eq. 5
        at = t-1
        for f = 1:length(graph.props[t-1])
            fluent = graph.props[t-1][f]
            act_inds = get_action_indices(graph, fluent, at, :add, :pre)
            if !isempty(act_inds)
                @constraint(model, sum([y[i,t] for i in act_inds]) >= x[f, fidx[2]+ad])
                # Eq. 6
                [@constraint(model, y[i,t]  <= x[f, fidx[2]+ad]) for i in act_inds]
            end

        # # Eq. 7
        # for f = 1:length(graph.props[t])
            # fluent = graph.props[t][f]
            act_inds = get_action_indices(graph, fluent, at, :del, :pre)
            if !isempty(act_inds)
                @constraint(model, sum([y[i,t] for i in act_inds]) >= x[f, fidx[2]+dl])
                #Eq. 8
                [@constraint(model, y[i,t]  <= x[f, fidx[2]+dl]) for i in act_inds]
            end

        # # Eq. 9
        # for f = 1:length(graph.props[t])
            # fluent = graph.props[t][f]
            act_inds = get_action_indices(graph, fluent, at, :pre, :del)
            if !isempty(act_inds)
                @constraint(model, sum([y[i,t] for i in act_inds]) >= x[f, fidx[2]+pa])
                #Eq. 10
                [@constraint(model, y[i,t]  <= x[f, fidx[2]+pa]) for i in act_inds]
            end

        # # Eq. 11
        # for f = 1:length(graph.props[t])
            # fluent = graph.props[t][f]
            act_inds = get_action_indices(graph, fluent, at, :pre, :anddel)
            if !isempty(act_inds)
                @constraint(model, sum([y[i,t] for i in act_inds]) == x[f, fidx[2]+pd])
            end
            # Eq. 12
            @constraint(model, x[f,fidx[2]+ad] + x[f,fidx[2]+mn] + x[f,fidx[2]+dl]
                            + x[f, fidx[2]+pd] <= 1)

            # Eq. 13
            @constraint(model, x[f,fidx[2]+pa] + x[f,fidx[2]+mn] + x[f,fidx[2]+dl]
                            + x[f, fidx[2]+pd] <= 1)

            # Eq. 14
            fidx_prev,_ = get_fluent_range_indices(graph, at, max_nf)
            @constraint(model, x[f,fidx[2]+pa] + x[f,fidx[2]+mn] + x[f,fidx[2]+pd]
                            <= x[f,fidx_prev[2]+pa] + x[f,fidx_prev[2]+ad] + x[f,fidx_prev[2]+mn])

            # one action at each time step
            num_actions = length([act for act in graph.acts[t] if act.name != :NoOp])
            @constraint(model, sum(y[1:num_actions, t]) == 1)
        end

    end

    #objective
    @objective(model, Min, sum(y))
    optimize!(model)
    sol = value.(y)
    render_action(sol, graph)
    # graph
    # value.(x)
    # return model
end

function render_action(y, graph)
    T = graph.num_levels
    actions = []
    for t=1:T
        ind = findall(x->x == 1.0, y[1:end,t])[1]
        act  = graph.acts[t][ind]
        push!(actions, (act.name, act.args))
    end
    println(actions)
    return actions

end


text:
using Lava: IFeature
mutable struct TestFeature <: IFeature
    mOnInstanceCreatedCalled
    mBeforeInstanceDestruction

    mOnLogicalDeviceCreated
    mOnPhysicalDeviceSelected
    mBeforeDeviceDestructionCalled

    function TestFeature()
        this = new()
        this.mOnInstanceCreatedCalled = false
        this.mBeforeInstanceDestruction = false

        this.mOnLogicalDeviceCreated = false
        this.mOnPhysicalDeviceSelected = false
        this.mBeforeDeviceDestructionCalled = false
        return this
    end
end

@class TestFeature

function create(::Type{TestFeature})
    return TestFeature()
end

function lava.:onInstanceCreated(this::TestFeature, vkInstance::vk.VkInstance)
    this.mOnInstanceCreatedCalled = true
    this.mBeforeInstanceDestruction = false
end

function lava.:beforeInstanceDestruction(this::TestFeature)
    this.mBeforeInstanceDestruction = true
    this.mOnInstanceCreatedCalled = false
end

function lava.:onLogicalDeviceCreated(this::TestFeature, device)
    this.mOnLogicalDeviceCreated = true
    this.mBeforeDeviceDestructionCalled = false
end

function lava.:onPhysicalDeviceSelected(this::TestFeature, phy::vk.VkPhysicalDevice)
    this.mOnPhysicalDeviceSelected = true
    this.mBeforeDeviceDestructionCalled = false
end

function lava.:beforeDeviceDestruction(this::TestFeature)
    this.mBeforeDeviceDestructionCalled = true
    this.mOnPhysicalDeviceSelected = false
    this.mOnLogicalDeviceCreated = false
end


text:
using DH32Utils
using Base.Test
using GZip

# write your own tests here
x = vec(readcsv(gzopen("testdist.csv.gz")))
w = vec(readcsv(gzopen("testweights.csv.gz")))
q1,q2 = block_average_spread(x, dimin=1000, ddi=1000, w=w)
@test abs(q1[1]-1.7432) < 0.001
@test abs(q2[1]-2.0877) < 0.001


text:
export NEATInd
import Base: print

mutable struct Connection
    in_node::Int
    out_node::Int
    weight::Float64
    enabled::Bool
    innovation::Int
end

mutable struct Neuron
    input::Float64
    output::Float64
    processed::Bool
end

struct NEATInd <: Cambrian.Individual
    n_in::Int
    n_out::Int
    n_hidden::Int
    connections::Array{Connection}
    neurons::Array{Neuron}
    fitness::Array{Float64}
end

function rand_weight()
    rand() * 2.0 - 1.0 # initial weight is uniformly distributed between -1 and 1
end

function NEATInd(cfg::Dict)
    innovation = cfg["innovation_max"]
    n_in = cfg["n_in"]
    n_out = cfg["n_out"]
    connections = Array{Connection}(undef, 0)
    for i in 1:n_in
        for j in 1:n_out
            innovation += 1
            c = Connection(i, n_in + j, rand_weight(), true, innovation)
            push!(connections, c)
        end
    end
    cfg["innovation_max"]= innovation
    neurons = Array{Neuron}(undef, 0)
    for i in 1:(n_in+n_out)
        push!(neurons, Neuron(0.0, 0.0, false))
    end
    fitness = -Inf .* ones(cfg["d_fitness"])
    NEATInd(n_in, n_out, 0, connections, neurons, fitness)
end

"""clone an individual and give new fitness values"""
function NEATInd(ind::NEATInd)
    ind2 = deepcopy(ind)
    for n in ind2.neurons
        n.input = 0.0
        n.output = 0.0
    end
    ind2.fitness .= -Inf
    ind2
end

# TODO: remove, already in abstracting branch of Cambrian.jl
function print(io::IO, ind::NEATInd)
    print(io, JSON.json(ind))
end

function String(ind::NEATInd)
    string(ind)
end

function NEATInd(cfg::Dict, ind_s::String)
    d = JSON.parse(ind_s)
    n_in = d["n_in"]
    n_out = d["n_out"]
    n_hidden = d["n_hidden"]
    connections = Array{Connection}(undef, 0)
    for ci in d["connections"]
        c = Connection(ci["in_node"], ci["out_node"],
                       ci["weight"], ci["enabled"], ci["innovation"])
        push!(connections, c)
    end
    neurons = Array{Neuron}(undef, 0)
    for ni in d["neurons"]
        push!(neurons, Neuron(ni["input"], ni["output"], ni["processed"]))
    end
    fitness = zeros(length(d["fitness"]))
    for i in eachindex(d["fitness"])
        if d["fitness"][i] == nothing
            fitness[i] = -Inf
        else
            fitness[i] = Float64(d["fitness"][i])
        end
    end
    NEATInd(n_in, n_out, n_hidden, connections, neurons, fitness)
end


text:
module VolterraGP
	using Base
	using LinearAlgebra
	using Random
	using Distributions
	using DistributionsAD
	using Plots
	using QuadGK
	using Flux
	using Flux:@adjoint
	using ForwardDiff
	using TensorCast
	using Tullio

	export Data
	export GaussianProcess, DiffableParameters
	export negloglikelihood, posterior, threeEQs, scaledEQs, plotgp, generate_toy_data, fit!
	export NMSE

	include("gp.jl")
	include("GGu.jl")
	include("volterra.jl")
	include("plot.jl")
	include("data/load.jl")
	include("fit.jl")
	include("utils.jl")
	include("metrics.jl")
end


text:
export LogisticBinomial, HalfCauchy, HalfNormal
import Distributions.logpdf
# using SymPy

struct HalfCauchy{T} <: Distribution{Univariate,Continuous}
    σ::T

    function HalfCauchy{T}(σ::T) where T
        new{T}(σ)
    end
end

HalfCauchy(σ::T) where {T<:Real} = HalfCauchy{T}(σ)
HalfCauchy(σ::Integer) = HalfCauchy(Float64(σ))
HalfCauchy() = HalfCauchy(1.0)

Distributions.logpdf(d::HalfCauchy{T} ,x::Real) where {T}  = log(2.0) + logpdf(Cauchy(0.0,d.σ),x)

Distributions.pdf(d::HalfCauchy,x) = 2 * pdf(Cauchy(0,d.σ),x)

Distributions.rand(d::HalfCauchy) = abs(rand(Cauchy(0,d.σ)))

Distributions.quantile(d::HalfCauchy, p) = quantile(Cauchy(0, d.σ), (p+1)/2)

Distributions.support(::HalfCauchy{T} where T) = RealInterval(0.0, Inf)

struct HalfNormal
    σ
end

HalfNormal() = HalfNormal(1)

Distributions.logpdf(d::HalfNormal,x::Real) = log(2) + logpdf(Normal(0,d.σ),x)

Distributions.pdf(d::HalfNormal,x) = 2 * pdf(Normal(0,d.σ),x)

Distributions.rand(d::HalfNormal) = abs(rand(Normal(0,d.σ)))

Distributions.support(::HalfNormal) = RealInterval(0.0, Inf)

# HalfNormal(s) = Truncated(Normal(0,s),0,Inf)
# HalfNormal() = Truncated(Normal(0,1.0),0,Inf)

# Binomial distribution, parameterized by logit(p)
LogisticBinomial(n,x)=Binomial(n,logistic(x))


text:
using ModelingToolkit
using Base.Test

# `Expr`, `Number` -> `Operation`
@IVar a
@Param b
@DVar x(t)
@Var y
@test parse(Operation, 2) == 2
expr = :(-inv(2sqrt(+($a, $b))))
op   = Operation(-, [Operation(inv,
                [Operation(*, [2, Operation(sqrt,
                              [Operation(+, [a, b])])])])])
@test parse(Operation, expr) == op
expr1 = :($x^($y-1))
op1   = Operation(^, [x, Operation(-, [y, 1])])
@test parse(Operation, expr1) == op1


text:
using FractionalDiffEq, Plots, SpecialFunctions

phi(x)=0
α=0.9; τ = 0.1
f(t, y, ϕ)=2/gamma(3-α)*t^(2-α)-1/gamma(2-α)*t^(1-α)+2*τ*t-τ^2-τ-y+ϕ

prob=FDDEProblem(f, phi, α, τ, (2, 5))
sol=solve(prob, 0.01, DelayPI())

plot(collect(2:0.01:5), sol)

plot!(t->t^2-t, ls=:dash)


text:
function check_folder_integrity(folder::String)
    folder_files = readdir(folder)
    alien_files = filter(!∈(KNOWN_SIMULATION_PATHS), folder_files)
    if isempty(alien_files)
        return true
    end
    if "data_store" ∉ folder_files
        error("The file path doesn't contain any data_store folder")
    end
    return false
end

function _fill_result_value_container(fields)
    return FieldResultsByTime(x => ResultsByTime() for x in fields)
end

struct SimulationResults
    path::String
    params::SimulationStoreParams
    decision_problem_results::Dict{
        String,
        SimulationProblemResults{DecisionModelSimulationResults},
    }
    emulation_problem_results::SimulationProblemResults{EmulationModelSimulationResults}
    store::Union{Nothing, SimulationStore}
end

"""
Construct SimulationResults from a simulation output directory.

# Arguments

  - `path::AbstractString`: Simulation output directory
  - `execution::AbstractString`: Execution number. Default is the most recent.
  - `ignore_status::Bool`: If true, return results even if the simulation failed.
"""
function SimulationResults(path::AbstractString, execution=nothing; ignore_status=false)
    # path will be either the execution_path or the directory containing all executions.
    contents = readdir(path)
    if "data_store" in contents
        execution_path = path
    else
        if execution === nothing
            executions = [parse(Int, f) for f in contents if occursin(r"^\d+$", f)]
            if isempty(executions)
                error("There are no simulation results in the path")
            end
            execution = maximum(executions)
        end
        execution_path = joinpath(path, string(execution))
        if !isdir(execution_path)
            error("Execution $execution not in the simulations results")
        end
    end

    status = deserialize_status(joinpath(execution_path, RESULTS_DIR))
    _check_status(status, ignore_status)

    if !check_folder_integrity(execution_path)
        @warn "The results folder $(execution_path) is not consistent with the default folder structure. " *
              "This can lead to errors or unwanted results."
    end

    simulation_store_path = joinpath(execution_path, "data_store")
    check_file_integrity(simulation_store_path)

    return open_store(HdfSimulationStore, simulation_store_path, "r") do store
        decision_problem_results =
            Dict{String, SimulationProblemResults{DecisionModelSimulationResults}}()
        sim_params = get_params(store)
        container_key_lookup = get_container_key_lookup(store)
        for (name, problem_params) in sim_params.decision_models_params
            name = string(name)
            problem_result = SimulationProblemResults(
                DecisionModel,
                store,
                name,
                problem_params,
                sim_params,
                execution_path,
                container_key_lookup,
            )
            decision_problem_results[name] = problem_result
        end

        emulation_params = get_emulation_model_params(sim_params)
        emulation_result = SimulationProblemResults(
            EmulationModel,
            store,
            string(first(keys(sim_params.emulation_model_params))),
            first(values(sim_params.emulation_model_params)),
            sim_params,
            execution_path,
            container_key_lookup,
        )

        return SimulationResults(
            execution_path,
            sim_params,
            decision_problem_results,
            emulation_result,
            nothing,
        )
    end
end

"""
Construct SimulationResults from a simulation.
"""
function SimulationResults(sim::Simulation; ignore_status=false, kwargs...)
    _check_status(get_simulation_status(sim), ignore_status)
    store = get_simulation_store(sim)
    execution_path = get_simulation_dir(sim)
    decision_problem_results =
        Dict{String, SimulationProblemResults{DecisionModelSimulationResults}}()
    sim_params = get_params(store)
    models = get_models(sim)
    container_key_lookup = get_container_key_lookup(store)
    for (name, problem_params) in sim_params.decision_models_params
        model = get_simulation_model(models, name)
        name = string(name)
        problem_result = SimulationProblemResults(
            DecisionModel,
            store,
            name,
            problem_params,
            sim_params,
            execution_path,
            container_key_lookup,
            system=get_system(model),
        )
        decision_problem_results[name] = problem_result
    end

    emulation_model = get_emulation_model(models)
    emulation_results = SimulationProblemResults(
        EmulationModel,
        store,
        string(first(keys(sim_params.emulation_model_params))),
        first(values(sim_params.emulation_model_params)),
        sim_params,
        execution_path,
        container_key_lookup,
        system=isnothing(emulation_model) ? nothing : get_system(emulation_model),
    )

    return SimulationResults(
        execution_path,
        sim_params,
        decision_problem_results,
        emulation_results,
        store,
    )
end

function Base.empty!(res::SimulationResults)
    foreach(empty!, values(res.decision_problem_results))
    empty!(res.emulation_problem_results)
end

Base.isempty(res::SimulationResults) = all(isempty, values(res.decision_problem_results))
Base.length(res::SimulationResults) =
    mapreduce(length, +, values(res.decision_problem_results))
get_exports_folder(x::SimulationResults) = joinpath(x.path, "exports")

function get_decision_problem_results(results::SimulationResults, problem)
    if !haskey(results.decision_problem_results, problem)
        throw(IS.InvalidValue("$problem is not stored"))
    end

    return results.decision_problem_results[problem]
end

function get_emulation_problem_results(results::SimulationResults)
    return results.emulation_problem_results
end

"""
Return the problem names in the simulation.
"""
list_decision_problems(results::SimulationResults) =
    collect(keys(results.decision_problem_results))

"""
Export results to files in the results directory.

# Arguments

  - `results::SimulationResults`: simulation results
  - `exports`: SimulationResultsExport or anything that can be passed to its constructor.
    (such as Dict or path to JSON file)

An example JSON file demonstrating possible options is below. Note that `start_time`,
`end_time`, `path`, and `format` are optional.

```
{
  "decision_models": [
    {
      "name": "ED",
      "variables": [
        "P__ThermalStandard",
        "E__HydroEnergyReservoir"
      ],
      "parameters": [
        "all"
      ]
    },
    {
      "name": "UC",
      "variables": [
        "On__ThermalStandard"
      ],
      "parameters": [
        "all"
      ],
      "duals": [
        "all"
      ]
    }
  ],
  "start_time": "2020-01-01T04:00:00",
  "end_time": null,
  "path": null,
  "format": "csv"
}

```
"""
function export_results(results::SimulationResults, exports)
    if results.store isa InMemorySimulationStore
        export_results(results, exports, results.store)
    else
        simulation_store_path = joinpath(results.path, "data_store")
        open_store(HdfSimulationStore, simulation_store_path, "r") do store
            export_results(results, exports, store)
        end
    end
    return
end

function export_results(results::SimulationResults, exports, store::SimulationStore)
    if !(exports isa SimulationResultsExport)
        exports = SimulationResultsExport(exports, results.params)
    end

    file_type = get_export_file_type(exports)

    for problem_results in values(results.decision_problem_results)
        problem_exports = get_problem_exports(exports, problem_results.problem)
        path =
            exports.path === nothing ? problem_results.results_output_folder : exports.path
        for timestamp in get_timestamps(problem_results)
            !should_export(exports, timestamp) && continue

            export_path = mkpath(joinpath(path, problem_results.problem, "variables"))
            for name in list_variable_names(problem_results)
                if should_export_variable(problem_exports, name)
                    dfs = read_variable(
                        problem_results,
                        name;
                        initial_time=timestamp,
                        count=1,
                        store=store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end

            export_path = mkpath(joinpath(path, problem_results.problem, "aux_variables"))
            for name in list_aux_variable_names(problem_results)
                if should_export_aux_variable(problem_exports, name)
                    dfs = read_aux_variable(
                        problem_results,
                        name;
                        initial_time=timestamp,
                        count=1,
                        store=store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end

            export_path = mkpath(joinpath(path, problem_results.problem, "parameters"))
            for name in list_parameter_names(problem_results)
                if should_export_parameter(problem_exports, name)
                    dfs = read_parameter(
                        problem_results,
                        name;
                        initial_time=timestamp,
                        count=1,
                        store=store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end

            export_path = mkpath(joinpath(path, problem_results.problem, "duals"))
            for name in list_dual_names(problem_results)
                if should_export_dual(problem_exports, name)
                    dfs = read_dual(
                        problem_results,
                        name;
                        initial_time=timestamp,
                        count=1,
                        store=store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end
        end

        export_path = mkpath(joinpath(path, problem_results.problem, "expression"))
        for name in list_expression_names(problem_results)
            if should_export_expression(problem_exports, name)
                dfs = read_expression(
                    problem_results,
                    name;
                    initial_time=timestamp,
                    count=1,
                    store=store,
                )
                export_result(file_type, export_path, name, timestamp, dfs[timestamp])
            end
        end

        if problem_exports.optimizer_stats
            export_path = joinpath(path, problem_results.problem, "optimizer_stats.csv")
            df = read_optimizer_stats(problem_results, store=store)
            export_result(file_type, export_path, df)
        end
    end
    return
end

function export_result(
    ::Type{CSV.File},
    path,
    key::OptimizationContainerKey,
    timestamp::Dates.DateTime,
    df::DataFrames.DataFrame,
)
    name = encode_key_as_string(key)
    export_result(CSV.File, path, name, timestamp, df)
    return
end

function export_result(
    ::Type{CSV.File},
    path,
    name::AbstractString,
    timestamp::Dates.DateTime,
    df::DataFrames.DataFrame,
)
    filename = joinpath(path, name * "_" * convert_for_path(timestamp) * ".csv")
    export_result(CSV.File, filename, df)
    return
end

function export_result(
    ::Type{CSV.File},
    path,
    key::OptimizationContainerKey,
    df::DataFrames.DataFrame,
)
    name = encode_key_as_string(key)
    export_result(CSV.File, path, name, df)
    return
end

function export_result(
    ::Type{CSV.File},
    path,
    name::AbstractString,
    df::DataFrames.DataFrame,
)
    filename = joinpath(path, name * ".csv")
    export_result(CSV.File, filename, df)
    return
end

function export_result(::Type{CSV.File}, filename, df::DataFrames.DataFrame)
    open(filename, "w") do io
        CSV.write(io, df)
    end

    @debug "Exported $filename"
    return
end

function _check_status(status::RunStatus, ignore_status)
    status == RunStatus.SUCCESSFUL && return

    if ignore_status
        @warn "Simulation was not successful: $status. Results may not be valid."
    else
        error(
            "Simulation was not successful: status = $status. Set ignore_status = true to override.",
        )
    end
    return
end


text:
using JFSM
using Base.Test

# Load driving data

metdata = readdlm(joinpath(dirname(@__FILE__), "..", "data", "met_CdP_0506.txt"));

# Load data from fortran run

valdata = readdlm(joinpath(dirname(@__FILE__), "..", "data", "out_CdP_0506.txt"));

# Model data

md = FsmType(0, 0, 0, 0, 0);

# Run model

hs_new = run_fsm(md, metdata);

# From fortran run

hs_old = valdata[:, 6];

# Tests

@test maximum(abs(hs_new-hs_old)) < 0.01


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayBackendHealthOnDemand <: SwaggerModel
    backendAddressPool::Any # spec type: Union{ Nothing, ApplicationGatewayBackendAddressPool } # spec name: backendAddressPool
    backendHealthHttpSettings::Any # spec type: Union{ Nothing, ApplicationGatewayBackendHealthHttpSettings } # spec name: backendHealthHttpSettings

    function ApplicationGatewayBackendHealthOnDemand(;backendAddressPool=nothing, backendHealthHttpSettings=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendHealthOnDemand, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(ApplicationGatewayBackendHealthOnDemand, Symbol("backendHealthHttpSettings"), backendHealthHttpSettings)
        setfield!(o, Symbol("backendHealthHttpSettings"), backendHealthHttpSettings)
        o
    end
end # type ApplicationGatewayBackendHealthOnDemand

const _property_map_ApplicationGatewayBackendHealthOnDemand = Dict{Symbol,Symbol}(Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("backendHealthHttpSettings")=>Symbol("backendHealthHttpSettings"))
const _property_types_ApplicationGatewayBackendHealthOnDemand = Dict{Symbol,String}(Symbol("backendAddressPool")=>"ApplicationGatewayBackendAddressPool", Symbol("backendHealthHttpSettings")=>"ApplicationGatewayBackendHealthHttpSettings")
Base.propertynames(::Type{ ApplicationGatewayBackendHealthOnDemand }) = collect(keys(_property_map_ApplicationGatewayBackendHealthOnDemand))
Swagger.property_type(::Type{ ApplicationGatewayBackendHealthOnDemand }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendHealthOnDemand[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendHealthOnDemand }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendHealthOnDemand[property_name]

function check_required(o::ApplicationGatewayBackendHealthOnDemand)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendHealthOnDemand }, name::Symbol, val)
end


text:
"""
    Physical constants
    ==================
"""
# speed of light
const c = 299792458
# elementary charge
const q = 1.6021766208e-19
# atomic mass unit
const u = 1.66053904e-27
# electron mass
const m_e = 9.10938356e-31
# Plank constant
const h = 6.62607004e-34
# Bohr radius
const a0 = 5.2917721067e-11
# Rydberg constant
const Ry = 10973731.6
# Boltzman constant
const kB = 1.38064852e-23


text:
const SUPPORTED_MOSD_INPUT_SETS = Union{
    MOI.Zeros,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    MOI.PositiveSemidefiniteConeTriangle,
}

const SUPPORTED_NON_MOSD_INPUT_SETS = Union{
    MOI.ExponentialCone,
    MOI.DualExponentialCone,
    MOI.PowerCone,
    MOI.DualPowerCone
}

const SUPPORTED_MOSD_SETS = Union{
    SUPPORTED_MOSD_INPUT_SETS,
    MOI.Reals
}
const SUPPORTED_NON_MOSD_SETS = SUPPORTED_NON_MOSD_INPUT_SETS

const SUPPORTED_PROJ_SETS = Union{
    SUPPORTED_MOSD_SETS,
    SUPPORTED_NON_MOSD_SETS
}

# Same as cones supported by SCS
const SUPPORTED_INPUT_SETS = Union{
    SUPPORTED_MOSD_INPUT_SETS,
    SUPPORTED_NON_MOSD_INPUT_SETS
}

const EXP_CONE_THRESH = 1e-8
const POW_CONE_THRESH = 1e-8

function project_onto_cone(x, cone_prod)
    # TODO: fix to be undefined.
    pi_x = zeros(length(x))
    # pi_x = Vector{typeof(x)}(undef, length(x))
    curr = 1
    for cone in cone_prod
        x_curr = @view(x[curr:curr+MOI.dimension(cone)-1])
        pi_x[curr:curr+MOI.dimension(cone)-1] .= _proj(x_curr, cone)
        curr += MOI.dimension(cone)
    end
    return pi_x
end

function _proj(x, cone)
    if typeof(cone) <: SUPPORTED_MOSD_SETS
        return MOSD.projection_on_set(MOSD.DefaultDistance(), x, cone)
    elseif typeof(cone) <: SUPPORTED_NON_MOSD_SETS
        return _proj_custom(x, cone)
    else
        throw(ArgumentError("Projections for $(typeof(cone)) are not supported"))
    end
end

function _proj_custom(x, cone)
    if typeof(cone) <: MOI.ExponentialCone
        _proj_exp_cone(x; dual=false)
    elseif typeof(cone) <: MOI.DualExponentialCone
        _proj_exp_cone(x; dual=true)

    # TODO: Power cone
    elseif typeof(cone) <: MOI.PowerCone
        _proj_pow_cone(x, cone.exponent; dual=false)
    elseif typeof(cone) <: MOI.DualPowerCone
        _proj_pow_cone(x, cone.exponent; dual=true)
    end
end

function _proj_exp_cone(v; dual=false)
    # v has dimension 3
    # Moreau’s decomposition theorem:
    #   v = P_K(v) + P_Ko(v) where K is primal cone and Ko is polar cone = -K*
    #   Dual cone proj: P_K*(v) = -P_Ko(-v)
    # See section 6.3 of https://web.stanford.edu/~boyd/papers/pdf/prox_algs.pdf
    v = dual ? -v : v

    if in_exp_cone(v)
        return dual ? zeros(3) : v
    elseif in_exp_cone_dual(v)
        return dual ? -v : zeros(3)
    elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
        return dual ? [0.0; -v[2]; max(-v[3],0)] : [v[1]; 0.0; max(v[3],0)]
    else
        return get_exp_proj_case4(v; dual=dual)
    end
end

function in_exp_cone(v)
    return (
        (v[1] <= 0 && abs(v[2]) <= EXP_CONE_THRESH && v[3] >= 0) ||
        (v[2] > 0 && v[2] * exp(v[1] / v[2]) - v[3] <= EXP_CONE_THRESH)
    )
end

function in_exp_cone_dual(v)
    return (
        (abs(v[1]) <= EXP_CONE_THRESH && v[2] <= 0 && v[3] <= 0) ||
        (v[1] > 0 && v[1]*exp(v[2]/v[1]) + ℯ*v[3] <= EXP_CONE_THRESH)
    )
end

function get_exp_proj_case4(v; dual=false)
    # Bisection method for root finding of h(x) from
    # https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf
    # Thm: h(x) is smooth, strictly increasing, and changes sign
    @views r, s, t = v[1], v[2], v[3]

    h(x) = (((x-1)*r + s) * exp(x) - (r - x*s)*exp(-x))/(x^2 - x + 1) - t

    # Note: these won't both be infinity (by case 3)
    # TODO: figure out a better replacement for Inf
    if r > 0
        lb = 1 - s/r
    else
        lb = -1
        for _ in 1:10
            h(lb) < 0 && break
            lb *= 10
        end
    end
    if s > 0
        ub = r/s
    else
        ub = 1
        for _ in 1:10
            h(ub) > 0 && break
            ub *= 10
        end
    end

    (h(lb) > 0 || h(ub)) < 0 && error("Numerical error in Proj_exp_cone")

    # ∇h(x) = (exp(x)*(r*(x - 1) + s) + exp(-x)*(r - s*x) + r*exp(x) +
    #         s*exp(-x))/(x^2 - x + 1) - ((2x - 1)*(exp(x)*(r*(x - 1) + s) -
    #         exp(-x)*(r - s*x)))/(x^2 - x + 1)^2

    # TODO: anything more efficient than bisection search?
    x = find_zero(h, (lb,ub), verbose=false)
    vp = ((x - 1)*r + s)/(x^2 - x + 1) * [x; 1.0; exp(x)]
    # vd = (r - x*s)/(x^2 - x + 1) * [1.0; 1.0-x; -exp(-x)]

    return dual ? -v + vp : vp
end

# K = {(x,y,z) | x^a * y^(1-a) >= |z|, x>=0, y>=0}
# K* = {(u,v,w) | (u/a)^a * (v/(1-a))^(1-a) >= |w|, u>=0, v>=0}
function _proj_pow_cone(v, α; dual=false)
    # Prop 2.2 of https://link.springer.com/article/10.1007/s00186-015-0514-0
    # v = dual ? -v : v
    x, y, z = dual ? -v : v

    if in_pow_cone(x, y, z, α)
        # if in power cone
        ret = [x; y; z]
    elseif in_pow_cone_polar(x, y, z, α)
        # if in polar cone Ko = -K*
        ret = zeros(3)
    elseif abs(z) <= 1e-6#POW_CONE_THRESH
        # if not in K, Ko and z = 0
        ret = [max(x,0); max(y,0); 0.0]
    else
        ret = _proj_pow_cone_case4(x, y, z, α)
    end

    return dual ? v + ret : ret
end

function in_pow_cone(x, y, z, α)
    return x >= 0 && y >= 0 && POW_CONE_THRESH + x^α * y^(1-α) >= abs(z)
end

function in_pow_cone_polar(x, y, z, α)
    # -v = -[x;y;z] in K*(α)
    return (x <= 0 && y <=0 &&
        POW_CONE_THRESH + (-x)^α * (-y)^(1-α) >= α^α * (1-α)^(1-α) * abs(z))
end

function _proj_pow_cone_case4(x, y, z, α)
    Phi(x,y,z,α,r) = 0.5*(Phi_prod(x,α,z,r)^α * Phi_prod(y,1-α,z,r)^(1-α)) - r
    Phi_prod(xi,αi,z,r) = (xi + sqrt(xi^2 + 4*αi*r*(abs(z) - r)))

    lb_ub = (0.0+POW_CONE_THRESH, abs(z)-POW_CONE_THRESH)
    # println(Phi(x,y,z,α,lb_ub[1]), Phi(x,y,z,α,lb_ub[2]))
    r = find_zero(r -> Phi(x,y,z,α,r), lb_ub, verbose=false)

    return [0.5*Phi_prod(x,α,z,r); 0.5*Phi_prod(y,1-α,z,r); sign(z)*r]
end

# TODO: can probably write this code more efficiently? lots of repeats
function d_project_onto_cone(x, cone_prod)
    # TODO: fix to be undefined.
    ret = Vector{SparseMatrixCSC}(undef, length(cone_prod))
    # pi_x = Vector{typeof(x)}(undef, length(x))
    curr = 1
    for (ind, cone) in enumerate(cone_prod)
        ret[ind] = sparse(_d_proj(@view(x[curr:curr+MOI.dimension(cone)-1]), cone))
        curr += MOI.dimension(cone)
    end
    return blockdiag(ret...)
end

function _d_proj(x, cone)
    if typeof(cone) <: SUPPORTED_MOSD_SETS
        return MOSD.projection_gradient_on_set(MOSD.DefaultDistance(), x, cone)
    elseif typeof(cone) <: SUPPORTED_NON_MOSD_SETS
        return _d_proj_custom(x, cone)
    else
        throw(ArgumentError("Cone type $(typeof(cone)) is not supported."))
    end
end

function _d_proj_custom(x, cone)
    if typeof(cone) <: MOI.ExponentialCone
        _d_proj_exp_cone(x; dual=false)
    elseif typeof(cone) <: MOI.DualExponentialCone
        _d_proj_exp_cone(x; dual=true)

    # TODO: Power cone
    elseif typeof(cone) <: MOI.PowerCone
        _d_proj_pow_cone(x, cone.exponent; dual=false)
    elseif typeof(cone) <: MOI.DualPowerCone
        _d_proj_pow_cone(x, cone.exponent; dual=false)
    end
end

function _d_proj_exp_cone(v; dual=false)
    # Lemma 3.6 https://arxiv.org/pdf/1705.00772.pdf
    # J_PK*(v) = I - J_PK(-v)
    v = dual ? -v : v
    Ip(z) = z >= 0 ? 1 : 0

    if in_exp_cone(v)
        return dual ? zeros(3,3) : Matrix{Float64}(I, 3, 3)
    elseif in_exp_cone_dual(v)
        return dual ? Matrix{Float64}(I, 3, 3) : zeros(3,3)
    elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
        vp = diagm([1; Ip(v[2]); Ip(v[3])])
        return dual ? I - vp : vp
    else
        # Have to look at KKT conditions of projection problem
        return get_exp_d_proj_case4(v; dual=dual)
    end
end

function get_exp_d_proj_case4(v; dual=false)
    Pv = get_exp_proj_case4(v, dual=false)
    nu = Pv[3] - v[3]
    z1,z2,z3 = Pv
    abs(z2) <= 100eps() && error("Numerical error in projection onto exp cone")
    rs = z1/z2
    exp_rs = exp(rs)

    mat = inv([
        1+nu*exp_rs/z2     -nu*exp_rs*rs/z2       0     exp_rs;
        -nu*exp_rs*rs/z2   1+nu*exp_rs*rs^2/z2    0     (1-rs)*exp_rs;
        0                  0                      1     -1
        exp_rs             (1-rs)*exp_rs          -1    0
    ])
    return @view(mat[1:3,1:3])
end

function _d_proj_pow_cone(v, α; dual=false)
    # Thm 3.1 of https://link.springer.com/article/10.1007/s00186-015-0514-0
    x, y, z = dual ? -v : v

    if in_pow_cone(x, y, z, α)
        # if in power cone
        ret = Matrix{Float64}(I, 3, 3)
    elseif in_pow_cone_polar(x, y, z, α)
        # if in polar cone Ko = -K*
        ret = zeros(3,3)
    elseif abs(z) <= POW_CONE_THRESH
        # if not in K, Ko and z = 0
        ret = _d_proj_pow_cone_case3(x, y, z, α)
    else
        ret = _d_proj_pow_cone_case4(x, y, z, α)
    end

    return dual ? Matrix{Float64}(I, 3, 3) - ret : ret
end

function _d_proj_pow_cone_case3(x, y, z, α)
    v = [x; y]
    I(t) = t > 0 ? 1 : 0
    αs = [α; 1-α]

    if sum(αs[v .> 0]) > sum(αs[v .< 0])
        d = 1
    elseif sum(αs[v .> 0]) < sum(αs[v .< 0])
        d = 0
    else
        num = reduce(*, (-v[v .< 0]).^αs[v .< 0])
        denom = reduce(*, v[v .> 0].^αs[v .> 0]) * reduce(*, αs[v .< 0].^αs[v .< 0])
        d = 1/((num/denom)^2 + 1)
    end
    return diagm([I(x), I(y), d])
end

function _d_proj_pow_cone_case4(x, y, z, α)
    Phi(x,y,z,α,r) = 0.5*(Phi_prod(x,α,z,r)^α * Phi_prod(y,1-α,z,r)^(1-α)) - r
    Phi_prod(xi,αi,z,r) = (xi + sqrt(xi^2 + 4*αi*r*(abs(z) - r)))

    lb_ub = (0.0+POW_CONE_THRESH, abs(z)-POW_CONE_THRESH)
    r = find_zero(r -> Phi(x,y,z,α,r), lb_ub, verbose=false)

    za = abs(z)
    gx = sqrt(x^2 + 4*α*r*(za - r))
    gy = sqrt(y^2 + 4*(1-α)*r*(za - r))
    fx = 0.5*(x + gx)
    fy = 0.5*(y + gy)

    β = 1-α
    T_neg = (α*x/gx + β*y/gy)
    L = 2*(za - r) / (za + (za - 2r) * T_neg)
    T = -T_neg
    J_ii(w, γ, g) = 0.5 + w/2g + γ^2*(za - 2r)*r*L/g^2
    J_ij = α*β*(za - 2r)*r*L/(gx*gy)
    J = [
        J_ii(x, α, gx)      J_ij                sign(z)*α*r*L/gx;
        J_ij                J_ii(y, β, gy)      sign(z)*β*r*L/gy;
        sign(z)*α*r*L/gx    sign(z)*β*r*L/gy    r/za*(1+T*L)
    ]
    return J
end

function MOSD.projection_gradient_on_set(::MOSD.DefaultDistance, v::AbstractVector{T}, ::MOI.Reals) where {T}
    return Matrix{Float64}(I, (length(v), length(v)))
end

function MOSD.projection_gradient_on_set(::MOSD.DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}
    return zeros(T, (length(v), length(v)))
end


text:
using FixedSizeDictionaries
using BaseTestNext

function test_fixedkeyvalue(kvdict)

    @test kvdict[Val{:a}()] == 22
    @test @get(kvdict.b) == 3f0
    @test kvdict[Val{:c}()] == 3f0

    @test haskey(kvdict, Val{:a}()) == true
    @test haskey(kvdict, Val{:b}()) == true
    @test haskey(kvdict, Val{:c}()) == true
    @test haskey(kvdict, Val{:Y}()) == false

    val = get(kvdict, Val{:a}()) do
        "oh nose"
    end
    @test val == 22
    val = get(kvdict, Val{:Y}()) do
        "u do this"
    end
    @test val == "u do this"

end

@testset "FixedKeyValueDict" begin
    kvdict = FixedKeyValueDict((:a => 22, :b => 3f0, :c => 3f0))
    @test kvdict == FixedKeyValueDict(((:a, 22), (:b, 3f0), (:c, 3f0)))
    @test kvdict == FixedKeyValueDict(:a => 22, :b => 3f0, :c => 3f0)
    @test kvdict == FixedKeyValueDict((:a, :b, :c), (22, 3f0, 3f0))
    @test kvdict == FixedKeyValueDict(Dict(:a => 22, :b => 3f0, :c => 3f0))
    @test kvdict != FixedKeyDict(Dict(:a => 33, :b => 3f0, :c => 3f0))

    @test FixedSizeDictionaries.key2index(kvdict, Val{:a}()) == 1
    @test FixedSizeDictionaries.key2index(kvdict, Val{:b}()) == 2
    @test FixedSizeDictionaries.key2index(kvdict, Val{:c}()) == 3
    @test keys(kvdict) == (Val{:a}(), Val{:b}(), Val{:c}())
    @test values(kvdict) == (22, 3f0, 3f0)

    # test 2 times to make sure nothing gets compiled and different dicts
    # don't influence each other (method table is global state)
    test_fixedkeyvalue(kvdict)
    test_fixedkeyvalue(kvdict)

    #try the same with different order
    kvdict2 = FixedKeyValueDict((:b => 3f0, :c => 3f0, :a => 22))

    @test FixedSizeDictionaries.key2index(kvdict2, Val{:b}()) == 1
    @test FixedSizeDictionaries.key2index(kvdict2, Val{:c}()) == 2
    @test FixedSizeDictionaries.key2index(kvdict2, Val{:a}()) == 3
    @test keys(kvdict2) == (Val{:b}(), Val{:c}(), Val{:a}())
    @test values(kvdict2) == (3f0, 3f0, 22)

    test_fixedkeyvalue(kvdict2)
    test_fixedkeyvalue(kvdict2)

end

FixedSizeDictionaries.FixedKeyDict{Tuple{Val{:a},Val{:b},Val{:c}},Array{Any,1}}((Val{:a}(),Val{:b}(),Val{:c}()),Any[22,3.0,3.0]) != FixedSizeDictionaries.FixedKeyDict{Tuple{Val{:c},Val{:a},Val{:b}},Array{Any,1}}((Val{:c}(),Val{:a}(),Val{:b}()),Any[3.0,33,3.0])

@testset "FixedKeyDict" begin

    kvdict = FixedKeyDict((:a => 22, :b => 3f0, :c => 3f0))
    @test kvdict == FixedKeyDict(((:a, 22), (:b, 3f0), (:c, 3f0)))
    @test kvdict == FixedKeyDict(:a => 22, :b => 3f0, :c => 3f0)
    @test kvdict == FixedKeyDict((:a, :b, :c), [22, 3f0, 3f0])
    @test kvdict == FixedKeyDict(Dict(:a => 22, :b => 3f0, :c => 3f0))

    @test kvdict != FixedKeyDict(Dict(:a => 33, :b => 3f0, :c => 3f0))

    @test FixedSizeDictionaries.key2index(kvdict, Val{:a}()) == 1
    @test FixedSizeDictionaries.key2index(kvdict, Val{:b}()) == 2
    @test FixedSizeDictionaries.key2index(kvdict, Val{:c}()) == 3
    @test keys(kvdict) == (Val{:a}(), Val{:b}(), Val{:c}())
    @test keys(typeof(kvdict)) == (Val{:a}(), Val{:b}(), Val{:c}())
    @test values(kvdict) == [22, 3f0, 3f0]

    test_fixedkeyvalue(kvdict)
    test_fixedkeyvalue(kvdict)

    @get kvdict.a = 10
    @get kvdict.b = 77f0
    @get kvdict.c = 88f0

    @test kvdict[Val{:a}()] == 10
    @test kvdict[Val{:b}()] == 77f0
    @test kvdict[Val{:c}()] == 88f0

    @test get(kvdict, Val{:a}(), 999) == 10
    @test get(kvdict, Val{:b}(), 999) == 77f0
    @test get(kvdict, Val{:c}(), 999) == 88f0
    @test get(kvdict, Val{:a}(), 999) == 10

    @test get(kvdict, Val{:Y}(), 999) == 999

    val = get(kvdict, Val{:a}()) do
        "oh nose"
    end
    @test val == 10
    val = get(kvdict, Val{:Y}()) do
        "u do this"
    end
    @test val == "u do this"

    @test_throws ArgumentError @get(kvdict[x])
    # try
    #     @get(kvdict[x])
    # catch e
    #     @test isa(e, ArgumentError)
    # end

end


text:
using SparseDiffTools
using SparseArrays
using VertexSafeGraphs
using LightGraphs
using Random

const LG = LightGraphs

n = rand(10:100)
matrices = Array{SparseMatrixCSC, 1}(undef, 0)

for _ in 1:20
    n = rand(10:100)
    matrix = sprand(Int64, n, n, 0.5)
    push!(matrices, copy(matrix))
end

for i in 1:20
    matrix = matrices[i]
    g = matrix2graph(matrix, false)
    for e in edges(g)
        src = LG.src(e)
        dst = LG.dst(e)
        col1 = abs.(matrix[:, src])
        col2 = abs.(matrix[:, dst])
        pr = col1' * col2
        @test pr != 0
    end
end

for i in 1:20
    matrix = matrices[i]
    g = matrix2graph(matrix, true)
    for e in edges(g)
        src = LG.src(e)
        dst = LG.dst(e)
        row1 = abs.(matrix[src, :])
        row2 = abs.(matrix[dst, :])
        pr = row1' * row2
        @test pr != 0
    end
end

@info "stored zeros"
for i in 1:20
    matrix = matrices[i]
    g = matrix2graph(matrix)
    # recalculate graph with stored zeros
    matrix_sz = copy(matrix)
    fill!(matrix_sz, 0.0)
    g_sz = matrix2graph(matrix_sz)
    # check that graphs are the same
    @test nv(g) == nv(g_sz)
    @test ne(g) == ne(g_sz)
    for e in edges(g)
        @test has_edge(g_sz, e)
    end
    for e in edges(g_sz)
        @test has_edge(g, e)
    end
end


text:
const Optimals = Dict{Symbol,Float64}(
               [:a280 => 2579,
                :ali535 => 202339,
                :att48 => 10628,
                :att532 => 27686,
                :bayg29 => 1610,
                :bays29 => 2020,
                :berlin52 => 7542,
                :bier127 => 118282,
                :brazil58 => 25395,
                :brd14051 => 469385,
                :brg180 => 1950,
                :burma14 => 3323,
                :ch130 => 6110,
                :ch150 => 6528,
                :d198 => 15780,
                :d493 => 35002,
                :d657 => 48912,
                :d1291 => 50801,
                :d1655 => 62128,
                :d2103 => 80450,
                :d15112 => 1573084,
                :d18512 => 645238,
                :dantzig42 => 699,
                :dsj1000 => 18659688,
                :eil51 => 426,
                :eil76 => 538,
                :eil101 => 629,
                :fl417 => 11861,
                :fl1400 => 20127,
                :fl1577 => 22249,
                :fl3795 => 28772,
                :fnl4461 => 182566,
                :fri26 => 937,
                :gil262 => 2378,
                :gr17 => 2085,
                :gr21 => 2707,
                :gr24 => 1272,
                :gr48 => 5046,
                :gr96 => 55209,
                :gr120 => 6942,
                :gr137 => 69853,
                :gr202 => 40160,
                :gr229 => 134602,
                :gr431 => 171414,
                :gr666 => 294358,
                :hk48 => 11461,
                :kroA100 => 21282,
                :kroB100 => 22141,
                :kroC100 => 20749,
                :kroD100 => 21294,
                :kroE100 => 22068,
                :kroA150 => 26524,
                :kroB150 => 26130,
                :kroA200 => 29368,
                :kroB200 => 29437,
                :lin105 => 14379,
                :lin318 => 42029,
                :linhp318 => 41345,
                :nrw1379 => 56638,
                :p654 => 34643,
                :pa561 => 2763,
                :pcb442 => 50778,
                :pcb1173 => 56892,
                :pcb3038 => 137694,
                :pla7397 => 23260728,
                :pla33810 => 66048945,
                :pla85900 => 142382641,
                :pr76 => 108159,
                :pr107 => 44303,
                :pr124 => 59030,
                :pr136 => 96772,
                :pr144 => 58537,
                :pr152 => 73682,
                :pr226 => 80369,
                :pr264 => 49135,
                :pr299 => 48191,
                :pr439 => 107217,
                :pr1002 => 259045,
                :pr2392 => 378032,
                :rat99 => 1211,
                :rat195 => 2323,
                :rat575 => 6773,
                :rat783 => 8806,
                :rd100 => 7910,
                :rd400 => 15281,
                :rl1304 => 252948,
                :rl1323 => 270199,
                :rl1889 => 316536,
                :rl5915 => 565530,
                :rl5934 => 556045,
                :rl11849 => 923288,
                :si175 => 21407,
                :si535 => 48450,
                :si1032 => 92650,
                :st70 => 675,
                :swiss42 => 1273,
                :ts225 => 126643,
                :tsp225 => 3916,
                :u159 => 42080,
                :u574 => 36905,
                :u724 => 41910,
                :u1060 => 224094,
                :u1432 => 152970,
                :u1817 => 57201,
                :u2152 => 64253,
                :u2319 => 234256,
                :ulysses16 => 6859,
                :ulysses22 => 7013,
                :usa13509 => 19982859,
                :vm1084 => 239297,
                :vm1748 => 336556])


text:
module TORA

# Workaround for: https://github.com/JuliaRobotics/RigidBodyDynamics.jl/issues/500
using LinearAlgebra; BLAS.set_num_threads(1)

using Colors
using ForwardDiff
using GeometryBasics
using Ipopt
using MeshCat
using MeshCatMechanisms
using NPZ
using Plots
using Requires
using RigidBodyDynamics
using SparseArrays
using SparseDiffTools

using Plots.PlotMeasures: px
using Random: rand!

import Base: length

"""
    solve_with_knitro(problem, robot;
                      initial_guess=Float64[],
                      use_inv_dyn=false,
                      minimise_τ=false,
                      user_options=Dict())

Solve the nonlinear optimization problem with Knitro.

Further options can be set using the keyword arguments. See [Solver Interfaces](@ref).

# Keyword arguments
- `initial_guess::Array{Float64}=Float64[]`: the starting point for the solver.
- `use_inv_dyn::Bool=false`: if true, enables the use of inverse dynamics instead of forward dynamics.
- `minimise_τ::Bool=false`: if true, activates a cost function to minimize the joint torques.
- `user_options::Dict=Dict()`: the user options for Knitro.

See also: [`solve_with_ipopt`](@ref)
"""
function solve_with_knitro end

function __init__()
    @require KNITRO = "67920dd8-b58e-52a8-8622-53c4cffbe346" include("./transcription/knitro.jl")
end

greet() = print("Hello World!")

# `Struct`-defining files
include("./jacobian_data.jl")
include("./robot.jl")
include("./problem.jl")
include("./solverlog.jl")

# Regular source code
include("./constraints/dynamics.jl")
include("./constraints/end_effector.jl")
include("./transcription/ipopt.jl")
include("./plots.jl")

export solve_with_knitro

end # module


text:
# algorithm/z-test.jl

using Test

@test ncg(:test)
@test ogm_ls(:test)
@test pogm_restart(:test)


text:
module KoLesky

# Write your package code here.

include("./DOF.jl")
include("./Sorting.jl")
include("./SuperNodes.jl")

export mat2points
export maximin
export partition_into_supernodes

end


text:
using Test

include("util.jl")

@testset "hypergraph_external_nodes_count" begin
    using HyperGraphTools.HRGVisualization

    hg = HyperGraph{String}(HyperEdge[
        HyperEdge("John", ["_"], ["3"]),
        HyperEdge("LOVES",["3"], ["_"])
    ])
    external_node_count = HRGVisualization._external_node_count(hg)
    @test external_node_count == 2

    hg = HyperGraph{String}(HyperEdge[
        HyperEdge("Cookie Monster", ["_"], ["3"]),
        HyperEdge("EATS",["3"], ["_", "_"])
    ])
    external_node_count = HRGVisualization._external_node_count(hg)
    @test external_node_count == 3
end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: rhdekker/HyperGraphTools

max_forks_repo_name: rhdekker/HyperGraphTools


text:
"""
    space(::SiteType"S=1/2";
          conserve_qns = false,
          conserve_sz = conserve_qns,
          conserve_szparity = false,
          qnname_sz = "Sz",
          qnname_szparity = "SzParity")

Create the Hilbert space for a site of type "S=1/2".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function space(::SiteType"S=1/2";
               conserve_qns = false,
               conserve_sz = conserve_qns,
               conserve_szparity = false,
               qnname_sz = "Sz",
               qnname_szparity = "SzParity")
  if conserve_sz && conserve_szparity
    return [QN((qnname_sz, +1), (qnname_szparity, 1, 2)) => 1,
            QN((qnname_sz, -1), (qnname_szparity, 0, 2)) => 1]
  elseif conserve_sz
    return [QN(qnname_sz, +1) => 1,
            QN(qnname_sz, -1) => 1]
  elseif conserve_szparity
    return [QN(qnname_szparity, 1, 2) => 1,
            QN(qnname_szparity, 0, 2) => 1]
  end
  return 2
end

state(::SiteType"S=1/2", ::StateName"Up") = 1
state(::SiteType"S=1/2", ::StateName"Dn") = 2

state(st::SiteType"S=1/2", ::StateName"↑") =
  state(st, StateName("Up"))
state(st::SiteType"S=1/2", ::StateName"↓") =
  state(st, StateName("Dn"))

op(::OpName"Z",::SiteType"S=1/2") =
  [1  0
   0 -1]

op(::OpName"Sz",::SiteType"S=1/2") =
  [0.5  0.0
   0.0 -0.5]

op(::OpName"Sᶻ",t::SiteType"S=1/2") =
  op(OpName("Sz"),t)

op(::OpName"S+",::SiteType"S=1/2") =
  [0  1
   0  0]

op(::OpName"S⁺",t::SiteType"S=1/2") =
  op(OpName("S+"),t)

op(::OpName"Splus",t::SiteType"S=1/2") =
  op(OpName("S+"),t)

op(::OpName"S-",::SiteType"S=1/2") =
  [0  0
   1  0]

op(::OpName"S⁻",t::SiteType"S=1/2") =
  op(OpName("S-"),t)

op(::OpName"Sminus",t::SiteType"S=1/2") =
  op(OpName("S-"),t)

op(::OpName"X",::SiteType"S=1/2") =
  [0  1
   1  0]

op(::OpName"Sx",::SiteType"S=1/2") =
  [0.0  0.5
   0.5  0.0]

op(::OpName"Sˣ",t::SiteType"S=1/2") =
  op(OpName("Sx"),t)

op(::OpName"iY",::SiteType"S=1/2") =
  [ 0 1
   -1 0]

op(::OpName"iSy",::SiteType"S=1/2") =
  [ 0.0  0.5
   -0.5  0.0]

op(::OpName"iSʸ",t::SiteType"S=1/2") =
  op(OpName("iSy"),t)

op(::OpName"Y",::SiteType"S=1/2") =
  [0.0   -1.0im
   1.0im  0.0  ]

op(::OpName"Sy",::SiteType"S=1/2") =
  [0.0   -0.5im
   0.5im  0.0  ]

op(::OpName"Sʸ",t::SiteType"S=1/2") =
  op(OpName("Sy"),t)

op(::OpName"S2",::SiteType"S=1/2") =
  [0.75  0.0
   0.0   0.75]

op(::OpName"S²",t::SiteType"S=1/2") =
  op(OpName("S2"),t)

op(::OpName"ProjUp",::SiteType"S=1/2") =
  [1  0
   0  0]

op(::OpName"projUp",t::SiteType"S=1/2") =
  op(OpName("ProjUp"),t)

op(::OpName"ProjDn",::SiteType"S=1/2") =
  [0  0
   0  1]

op(::OpName"projDn",t::SiteType"S=1/2") =
  op(OpName("ProjDn"),t)

# Support the tag "SpinHalf" as equivalent to "S=1/2"

space(::SiteType"SpinHalf"; kwargs...) =
  space(SiteType("S=1/2"); kwargs...)

state(::SiteType"SpinHalf", n::StateName) =
  state(SiteType("S=1/2"), n)

op(o::OpName, ::SiteType"SpinHalf";kwargs...) =
  op(o, SiteType("S=1/2");kwargs...)

# Support the tag "S=½" as equivalent to "S=1/2"

space(::SiteType"S=½"; kwargs...) =
  space(SiteType("S=1/2"); kwargs...)

state(::SiteType"S=½", n::StateName) =
  state(SiteType("S=1/2"), n)

op(o::OpName, ::SiteType"S=½";kwargs...) =
  op(o, SiteType("S=1/2");kwargs...)


text:
using Distributed
@everywhere using Mamba

## Data
stacks = Dict{Symbol, Any}(
  :y => [42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9,
         15, 15],
  :x =>
    [80 27 89
     80 27 88
     75 25 90
     62 24 87
     62 22 87
     62 23 87
     62 24 93
     62 24 93
     58 23 87
     58 18 80
     58 18 89
     58 17 88
     58 18 82
     58 19 93
     50 18 89
     50 18 86
     50 19 72
     50 19 79
     50 20 80
     56 20 82
     70 20 91]
)
stacks[:N] = size(stacks[:x], 1)
stacks[:p] = size(stacks[:x], 2)

stacks[:meanx] = map(j -> mean(stacks[:x][:, j]), 1:stacks[:p])
stacks[:sdx] = map(j -> std(stacks[:x][:, j]), 1:stacks[:p])
stacks[:z] = Float64[
  (stacks[:x][i, j] - stacks[:meanx][j]) / stacks[:sdx][j]
  for i in 1:stacks[:N], j in 1:stacks[:p]
]

## Model Specification
model = Model(

  y = Stochastic(1,
    (mu, s2, N) ->
      UnivariateDistribution[Laplace(mu[i], s2) for i in 1:N],
    false
  ),

  beta0 = Stochastic(
    () -> Normal(0, 1000),
    false
  ),

  beta = Stochastic(1,
    () -> Normal(0, 1000),
    false
  ),

  mu = Logical(1,
    (beta0, z, beta) -> beta0 .+ z * beta,
    false
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001),
    false
  ),

  sigma = Logical(
    s2 -> sqrt(2.0) * s2
  ),

  b0 = Logical(
    (beta0, b, meanx) -> beta0 - dot(b, meanx)
  ),

  b = Logical(1,
    (beta, sdx) -> beta ./ sdx
  ),

  outlier = Logical(1,
    (y, mu, sigma, N) ->
      Float64[abs((y[i] - mu[i]) / sigma) > 2.5 for i in 1:N],
    [1, 3, 4, 21]
  )

)

## Initial Values
inits = [
  Dict(:y => stacks[:y], :beta0 => 10, :beta => [0, 0, 0], :s2 => 10),
  Dict(:y => stacks[:y], :beta0 => 1, :beta => [1, 1, 1], :s2 => 1)
]

## Sampling Scheme
scheme = [NUTS([:beta0, :beta]),
          Slice(:s2, 1.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, stacks, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
function main()
    set_of_terms = Set{BigInt}()
    for a in 2:100
        for b in 2:100
            push!(set_of_terms, big(a) ^ big(b))
        end
    end
    @printf("%d distinct terms\n", length(set_of_terms))
end

main()


text:
using SafeTestsets

@time @safetestset "optimisation" begin include("optimisation/optimisation.jl") end
@time @safetestset "ouput" begin include("optimisation/output.jl") end
@time @safetestset "integration" begin include("integration.jl") end
@time @safetestset "downsampling" begin include("downsampling.jl") end
@time @safetestset "layers" begin include("layers.jl") end
@time @safetestset "human dispersal" begin include("rules/human.jl") end
@time @safetestset "jump dispersal" begin include("rules/jump.jl") end
@time @safetestset "kernels" begin include("rules/kernels.jl") end
@time @safetestset "allee effects" begin include("rules/allee.jl") end
@time @safetestset "growth" begin include("rules/growth.jl") end


text:
mutable struct Config{T <: AbstractFloat}
    folds::Integer
    epochs::Integer
    bs::Integer
    lr::T
end

show(io::IO,c::Config) = print(io, "A trainer config.\n(epochs $(c.epochs), batch size $(c.bs), learn rate $(c.lr)).\nThere are $(c.folds) folds.")

mutable struct Status
    epoch::Integer
    loss::Dict{String,Array{Float64}}
end

function reset_model!(model;init=Flux.glorot_uniform)
    Flux.loadparams!(model, map(x -> init(size(x)), Flux.params(model)))
    model
end

struct Trainer
    model
    loss::Function
    data::Dataset
    folds::Int
    epochs::Int
    bs::Int
end

Trainer(model,loss,data;folds=5,epochs=100,batch_size=16) = Trainer(model, loss, data, folds, epochs, batch_size)

function training(trainer::Trainer, opt;log_dir::String="./result",use_gpu=true)
    RD = trainer.data.associations
    knowns = findall(RD .== 1)
    unknowns = findall(RD .== 0)
    shuffle!(knowns)
    shuffle!(unknowns)
    selected = eachindex(knowns)
    cv = CrossValidation(trainer.folds, selected)
    ps = params(trainer.model)
    plot_train = plot()
    plot_validate = plot()
    isdir("./result") || mkdir("./result")
    writedlm("./result/RD.csv", RD, ',')
    writedlm("./result/drugs.csv", trainer.data.drugs, ',')
    writedlm("./result/diseases.csv", trainer.data.diseases, ',')
    for f in 1:trainer.folds
        @info "fold $f"
        rd = copy(RD)
        loss_train = []
        loss_val = []
        train, test = cv(f)
        rd[knowns[test]] .= false
        mask = falses(size(RD))
        mask[knowns[test]] .= true
        useful = map(>=(1), sum(mask, dims=2) |> vec)
        mask = trues(size(RD))
        mask[knowns[train]] .= false
        mask[unknowns[train]] .= false
        reset_model!(trainer.model)
        samples = Samples(trainer.data, rd, RD, findall(!, mask);use_gpu=use_gpu,mode=:both)
        @info "training..."
        for e in 1:trainer.epochs
            loss_epoch = []
            val_epoch = []
            change_indices!(samples, findall(!, mask))
            loader = Flux.Data.DataLoader(samples, batchsize=trainer.bs)
            for (i, (x1, x2), y) in loader
                local loss
                gs = gradient(ps) do
                    ŷ = trainer.model(x1, x2)
                    loss = trainer.loss(ŷ, y)
                end
                Flux.update!(opt, ps, gs)
                push!(loss_epoch, loss)
            end
            push!(loss_train, mean(loss_epoch))
            validate = map(>=(4), sum(mask, dims=2) |> vec)
            change_indices!(samples, findall(mask .& validate))
            loader = Flux.Data.DataLoader(samples, batchsize=512)
            for (i, (x1, x2), y) in loader
                ŷ = trainer.model(x1, x2)
                loss = trainer.loss(ŷ, y)
                push!(val_epoch, loss)
            end
            push!(loss_val, mean(val_epoch))
        end
        plot!(plot_train, loss_train, label="fold $f")
        plot!(plot_validate, loss_val, label="fold $f")
        plot(plot_train,plot_validate,layout=(2,1))|>display
        @info "testing..."
        @info "$(sum(useful)) drugs to test."
        mask = mask .& useful
        @info "total $(sum(mask)) drug-disease pairs to test."
        change_indices!(samples, findall(mask))
        loader = Flux.Data.DataLoader(samples, batchsize=1024)
        result = zeros(size(RD))
        for (i, x, y) in loader
            ŷ = trainer.model(x...)
            ȳ = cpu(ŷ)
            result[i] .= softmax(ȳ)[2,:]
        end
        h = heatmap(result[useful,:])
        yticks!(h, [1:length(useful);], trainer.data.drugs[useful])
        display(h)
        isdir("./result/fold_$f") || mkdir("./result/fold_$f")
        writedlm("./result/fold_$f/result.csv", result, ',')
        writedlm("./result/fold_$f/mask.csv", mask, ',')
    end
end

export Config,Status,reset_model!,Trainer,training

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
using Test
using Manifolds

struct NotImplementedNumbers <: Manifolds.AbstractNumbers end

@testset "Number systems" begin
    @test_throws ErrorException real_dimension(NotImplementedNumbers())

    @test ℝ isa Manifolds.RealNumbers
    @test Manifolds.RealNumbers() === ℝ
    @test real_dimension(ℝ) == 1
    @test repr(ℝ) == "ℝ"
    @test ℝ^2 === Euclidean(2)
    @test ℝ^(2,3) === Euclidean(2, 3)

    @test ℂ isa Manifolds.ComplexNumbers
    @test Manifolds.ComplexNumbers() === ℂ
    @test real_dimension(ℂ) == 2
    @test repr(ℂ) == "ℂ"
    @test ℂ^2 === Euclidean(2; field = ℂ)
    @test ℂ^(2,3) === Euclidean(2, 3; field = ℂ)

    @test ℍ isa Manifolds.QuaternionNumbers
    @test Manifolds.QuaternionNumbers() === ℍ
    @test real_dimension(ℍ) == 4
    @test repr(ℍ) == "ℍ"
    @test ℍ^2 === Euclidean(2; field = ℍ)
    @test ℍ^(2,3) === Euclidean(2, 3; field = ℍ)
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

MyGame.Example.eval(quote

#=
#  Composite components of Monster color.
=#
@enum Color::UInt8 begin
    ColorRed = 1
#=
#  \brief color Green
#  Green is bit_flag with value (1u << 1)
=#
    ColorGreen = 2
#=
#  \brief color Blue (1u << 3)
=#
    ColorBlue = 8
end

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia