text:
module SCIP

include("scip_objsense.jl")
include("scip_retcode.jl")
include("scip_stage.jl")
include("scip_status.jl")
include("scip_vartype.jl")

export run_test

################################################################################
# Mappings to SCIP types
################################################################################
typealias SCIP_Real Float64
typealias SCIP_Bool Bool

typealias PtrSCIP Ptr{Void}
typealias PtrPtrSCIP Ptr{Ptr{Void}}

typealias PtrSCIP_Var Ptr{Void}
typealias PtrPtrSCIP_Var Ptr{Ptr{Void}}

typealias PtrSCIP_Cons Ptr{Void}
typealias PtrPtrSCIP_Cons Ptr{Ptr{Void}}

typealias PtrSCIP_Sol Ptr{Void}
typealias PtrPtrSCIP_Sol Ptr{Ptr{Void}}

################################################################################
# Macros for interacting with SCIP
################################################################################
# Macro for calling SCIP functions that return misc. types
macro scip_ccall(func, args...)
    f = "SCIP$(func)"
    return quote
        ccall(($f, "libscipopt"), $(args...))
    end
end

# Macro for calling SCIP functions that have checked return codes
macro scip_ccall_check(func, args...)
    f = "SCIP$(func)"
    return quote
        ret = ccall(($f, "libscipopt"), Int8, $(args...))
        if ret != SCIP_OKAY
            error(SCIP_RETCODE[ret])
        end
    end
end

################################################################################
# Miscellaneous methods
################################################################################
version() = @scip_ccall("version", SCIP_Real, ())
major_version() = @scip_ccall("majorVersion", Int, ())
minor_version() = @scip_ccall("minorVersion", Int, ())
tech_version() = @scip_ccall("techVersion", Int, ())
subversion() = @scip_ccall("subversion", Int, ())
# TODO: SCIPprintVersion
print_error(e) = @scip_ccall("printError", Void, (SCIP_Retcode,), e)
store_solution_gap(scip) = @scip_ccall("storeSolutionGap", Void, (PtrSCIP,), scip)

################################################################################
# General SCIP methodds
################################################################################
create(scip) = @scip_ccall_check("create", (PtrPtrSCIP,), scip) # TODO: Add !?
free(scip) = @scip_ccall_check("free", (PtrPtrSCIP,), scip)     # TODO: Add !?
get_stage(scip) = @scip_ccall("getStage", SCIP_Stage, (PtrSCIP,), scip[1])
# TODO: SCIPprintStage
get_status(scip) = @scip_ccall("getStatus", SCIP_Status, (PtrSCIP,), scip[1])
# TODO: SCIPprintStatus
is_transformed(scip) = @scip_ccall("isTransformed", SCIP_Bool, (PtrSCIP,), scip[1])
is_exact_solve(scip) = @scip_ccall("isExactSolve", SCIP_Bool, (PtrSCIP,), scip[1])
is_presolve_finished(scip) = @scip_ccall("isPresolveFinished", SCIP_Bool, (PtrSCIP,), scip[1])
pressed_ctrl_c(scip) = @scip_ccall("pressedCtrlC", SCIP_Bool, (PtrSCIP,), scip[1])
is_stopped(scip) = @scip_ccall("isStopped", SCIP_Bool, (PtrSCIP,), scip[1])

################################################################################
# More methods to get a basic example working
################################################################################
include_default_plugins!(scip) = @scip_ccall_check("includeDefaultPlugins", (PtrSCIP,), scip[1])

create_prob_basic!(scip, name) = @scip_ccall_check("createProbBasic", (PtrSCIP, String), scip[1], name)
free_prob!(scip) = @scip_ccall_check("freeProb", (PtrSCIP,), scip[1])

create_var_basic!(scip, var, name, lb, ub, obj, vartype) = @scip_ccall_check(
    "createVarBasic",
    (PtrSCIP, PtrPtrSCIP_Var, String, SCIP_Real, SCIP_Real, SCIP_Real, SCIP_Vartype),
    scip[1], var, name, lb, ub, obj, vartype
)
add_var!(scip, var) = @scip_ccall_check("addVar", (PtrSCIP, PtrSCIP_Var), scip[1], var[1])
release_var!(scip, var) = @scip_ccall_check("releaseVar", (PtrSCIP, PtrPtrSCIP_Var), scip[1], var)

# TODO: handle vars & vals conversion
include_conshdlr_linear!(scip) = @scip_ccall_check("includeConshdlrLinear", (PtrSCIP,), scip[1])
create_cons_basic_linear!(scip, cons, name, nvars, vars, vals, lhs, rhs) = @scip_ccall_check(
    "createConsBasicLinear",
    (PtrSCIP, PtrPtrSCIP_Cons, String, Int, PtrPtrSCIP_Var, Ptr{SCIP_Real}, SCIP_Real, SCIP_Real),
    scip[1], cons, name, nvars, vars, vals, lhs, rhs
)
add_cons!(scip, cons) = @scip_ccall_check("addCons", (PtrSCIP, PtrSCIP_Cons), scip[1], cons[1])
del_cons!(scip, cons) = @scip_ccall_check("delCons", (PtrSCIP, PtrSCIP_Cons), scip[1], cons[1])
release_cons!(scip, cons) = @scip_ccall_check("releaseCons", (PtrSCIP, PtrPtrSCIP_Cons), scip[1], cons)

set_objsense!(scip, objsense) = @scip_ccall_check("setObjsense", (PtrSCIP, SCIP_Objsense), scip[1], objsense)
presolve!(scip) = @scip_ccall_check("presolve", (PtrSCIP,), scip)
solve!(scip) = @scip_ccall_check("solve", (PtrSCIP,), scip[1])

get_best_sol(scip) = @scip_ccall("getBestSol", PtrSCIP_Sol, (PtrSCIP,), scip[1])
get_sol_orig_obj(scip, sol) = @scip_ccall("getSolOrigObj", SCIP_Real, (PtrSCIP, PtrSCIP_Sol), scip[1], sol)
get_sol_val(scip, sol, var) = @scip_ccall("getSolVal", SCIP_Real, (PtrSCIP, PtrSCIP_Sol, PtrSCIP_Var), scip[1], sol, var[1])

################################################################################
# Basic SCIP example - this will be wrapped up in an API
################################################################################
function run_test()
    scip = Array(Ptr{Void}, 1)
    create(scip)
    include_default_plugins!(scip)

    # Create a problem
    create_prob_basic!(scip, "Test problem")

    # Create a couple binary variables.
    x1 = Array(Ptr{Void}, 1)
    create_var_basic!(scip, x1, "x1", 0, 1, 4, SCIP_VARTYPE_BINARY)
    add_var!(scip, x1)

    x2 = Array(Ptr{Void}, 1)
    create_var_basic!(scip, x2, "x2", 0, 1, 3, SCIP_VARTYPE_BINARY)
    add_var!(scip, x2)

    # Add a constraint: 0 <= x1 + x2 <= 1
    cons = Array(Ptr{Void}, 1)
    vars = Array(Ptr{Void}, 2)
    vars[1] = x1[1]
    vars[2] = x2[1]
    vals = Array(SCIP_Real, 2)
    vals[1] = 1
    vals[2] = 1
    #include_conshdlr_linear!(scip)
    create_cons_basic_linear!(scip, cons, "c1", 2, vars, vals, 0, 1)
    add_cons!(scip, cons)

    # Objective sense = maximize
    set_objsense!(scip, SCIP_OBJSENSE_MAXIMIZE)
    solve!(scip)

    sol = get_best_sol(scip)
    println("\n========================================")
    println("MODEL SOLUTION:")
    println("OBJ = ", get_sol_orig_obj(scip, sol))
    println("X1 = ", get_sol_val(scip, sol, x1))
    println("X2 = ", get_sol_val(scip, sol, x2))
    println("========================================")

    # Cleaning up
    release_cons!(scip, cons)
    release_var!(scip, x1)
    release_var!(scip, x2)
    free_prob!(scip)
    free(scip)
end

end


text:
"""
    read_model_document(path_to_file::String;
            strip_comments::Bool = true)::Array{String,1}

Read each line of the Network.vff input file and extract contents of the file. Comments in the file (beginning by `//` symbol) are excluded by default, which can be overridden by the user.

Input arguments:
`path_to_file::String` - user-specified path where the Network.vff input file can be found and parsed
`strip_comments::Bool` - if true, comment lines (beginning by ‘//’ symbol) are excluded by the parser (optional).

Output arguments:
`buffer::Array{String,1}` - single array holding information extracted from the Network.vff file.

"""
function read_model_document(path_to_file::String;
    strip_comments::Bool = true)::Array{String,1}

    # initialize -
    buffer = String[]

    # Read in the file -
    open("$(path_to_file)", "r") do file
        for line in eachline(file)

            # exclude comments -
            if (occursin("//",line) == false && isempty(line) == false && strip_comments == true)
                +(buffer,line)
            end
        end
    end

    # return -
    return buffer
end

function transfer_distribution_file(path_to_distribution_files::String,
                                      input_file_name_with_extension::String,
                                      path_to_output_files::String,
                                      output_file_name_with_extension::String)

    # Load the specific file -
    # create src_buffer -
    src_buffer::Array{String} = String[]

    # check - do we have the file path?
    if (isdir(path_to_output_files) == false)
        mkpath(path_to_output_files)
    end

    # path to distrubtion -
    path_to_src_file = path_to_distribution_files*"/"*input_file_name_with_extension
    open(path_to_src_file,"r") do src_file
        for line in eachline(src_file)

            # need to add a new line for some reason in Julia 0.6
            new_line_with_line_ending = line*"\n"
            push!(src_buffer,new_line_with_line_ending)
        end
    end

    # Write the file to the output -
    path_to_program_file = path_to_output_files*"/"*output_file_name_with_extension
    open(path_to_program_file, "w") do f
        for line in src_buffer
            write(f,line)
        end
    end
end

function transfer_distribution_files(path_to_distribution_files::String,
                                      path_to_output_files::String,
                                      file_extension::String)

    # Search the directory for src files -
    # load the files -
    searchdir(path,key) = filter(x->contains(x,key),readdir(path))

    # build src file list -
    list_of_src_files = searchdir(path_to_distribution_files,file_extension)

    # check - do we have the file path?
    if (isdir(path_to_output_files) == false)
        mkpath(path_to_output_files)
    end

    # go thru the src file list, and copy the files to the output path -
    for src_file in list_of_src_files

        # create src_buffer -
        src_buffer::Array{String,1} = String[]

        # path to distrubtion -
        path_to_src_file = path_to_distribution_files*"/"*src_file
        open(path_to_src_file,"r") do src_file
            for line in eachline(src_file)

                # need to add a new line for some reason in Julia 0.6
                new_line_with_line_ending = line*"\n"
                push!(src_buffer,new_line_with_line_ending)
            end
        end

        # Write the file to the output -
        path_to_program_file = path_to_output_files*"/"*src_file
        open(path_to_program_file, "w") do f
            for line in src_buffer
                write(f,line)
            end
        end
    end
end

function write_program_components_to_disk(file_path::String, set_of_program_components::Set{NamedTuple})

    # check - do we have teh file path?
    if (isdir(file_path) == false)
      mkpath(file_path)
    end

    # go through each component, and dump the buffer to disk -
    for program_component in set_of_program_components

      # We switch on type -
      filename = program_component.filename
      component_type = program_component.component_type
      if (component_type == :buffer)

          # get the data -
          program_buffer = program_component.buffer

          # build the path -
          path_to_program_file = file_path*"/"*filename

          # Write the file -
          outfile = open(path_to_program_file, "w")
          write(outfile,program_buffer);
          close(outfile);

      elseif (component_type == :matrix || component_type == :vector)

          # get the matrix -
          program_matrix = program_component.matrix

          # build the path -
          path_to_program_file = file_path*"/"*filename

          # write the file -
          writedlm(path_to_program_file, program_matrix)

      else
          error("unsupported program component type: $(component_type)")
      end
    end
end

function move_existing_project_at_path(path_to_existing_project::String)::Bool

    # we are getting called *if* we already know there is a dir conflict -
    # if this is getting called, we have an existing dir where the user wants to write code.
    # we need then create a new dir called *.0, and mv the offending dir to this location?
    # return true if this worked, otherwise false -

    # parent and child dir that we are generating into -
    parent_dir = dirname(path_to_existing_project)
    child_dir = basename(path_to_existing_project)
    destination_path = ""

    # current backup index  -
    current_backup_index = 0

    # do we already have the destination?
    loop_flag = true
    while loop_flag

         # make a destination path -
        destination_path = joinpath(parent_dir,"$(child_dir).$(current_backup_index)")

        # we don't have this dir, we are done -
        if (isdir(destination_path) == false)
            loop_flag = false
        end

        # ok, looks like we already have this dir, update the counter -
        current_backup_index = current_backup_index + 1
    end

    # mv -
    mv(path_to_existing_project, destination_path)

    # check -
    if (isdir(destination_path) == false)
        return false
    end

    return true
end

"""
    generate_default_project(path_to_project_dir::String)::VLResult

Generates a default project structure which contains an empty model file and Defaults.toml file.

Input arguments:
`path_to_project_file::String` - path to where you want model code to be generated

Output arguments:
None
"""
function generate_default_project_file(path_to_project_file::String)::VLResult

    # ok, if we get here, then we have a clean place to generate the default project structure -
    # We need to two things for a project, the defaults file, and a blank network file with all the sections -
    filename = basename(path_to_project_file)
    path_to_output_files = dirname(path_to_project_file)

    # Transfer distrubtion files to the output -> these files are shared between model types -
    # TODO: replace w/joinpath?
    transfer_distribution_file("$(_PATH_TO_SRC)/distribution/dynamic", filename, path_to_output_files, filename)

    # ok, so return nothing if everything worked ...
    return VLResult(nothing)
end

ext: jl


text:
include("../src/casingcalculations.jl")

@testset "Casing pressure segment drop" begin
    sg_gas = 0.7
    P_pc = HankinsonWithWichertPseudoCriticalPressure(sg_gas, 0, 0)
    _, T_pc, _ = HankinsonWithWichertPseudoCriticalTemp(sg_gas, 0, 0)
    ΔP = calculate_casing_pressuresegment(300, 10, (100+180)/2, 4500, sg_gas, KareemEtAlZFactor, P_pc, T_pc)
    @test ΔP ≈ 332.7-300 atol = 1
end

@testset "Casing pressure traverse" begin
    md = [0, 2250, 4500]
    tvd = md
    inc = [0, 0, 0]
    id = 2.441

    temps = [100., 140., 180.]

    testwell = Wellbore(md, inc, tvd, id)

    pressures = casing_traverse_topdown(wellbore = testwell, temperatureprofile = temps,
                                        CHP = 300 - pressure_atmospheric, sg_gas = 0.7, dp_est = 10)

    @test pressures[end] ≈ (332.7 - pressure_atmospheric) atol = 5

    #model = WellModel(wellbore = testwell, temperatureprofile = temps, CHP = 300 - pressure_atmospheric, sg_gas_inj = 0.7, dp_est_inj = 10)
    #pressures2 = casing_traverse_topdown(model)
    #@test pressures2[end] ≈ (332.7 - pressure_atmospheric) atol = 5
end


text:
export NAU

"""
    NAU(in::Int, out::Int; init=glorot_uniform)

Neural addition unit.

As suggested in https://openreview.net/pdf?id=H1gNOeHKPS
"""
struct NAU{T<:AbstractMatrix}
    W::T
end

NAU(in::Int, out::Int; init=glorot_uniform) = NAU(init(out,in))

(m::NAU)(x) = m.W * x

Flux.@functor NAU

Base.show(io::IO, l::NAU) = print(io,"NAU(in=$(size(l.W,2)), out=$(size(l.W,1)))")


text:
"An `AbstractMappedGrid` represents the lazy application of a map to a grid."
abstract type AbstractMappedGrid{T,N} <: SimpleLazyGrid{T,N} end

const AbstractMappedGrid1d{T<:Number,N} = AbstractMappedGrid{T,N}

for op in (:length, :size, :eachindex, :indextype, :isperiodic)
	@eval $op(g::AbstractMappedGrid) = $op(supergrid(g))
end

covering(g::AbstractMappedGrid) = map_domain(forward_map(g), covering(supergrid(g)))

for op in (:minimum, :maximum)
	# TODO: the map could flip minimum and maximum...
	@eval $op(g::AbstractMappedGrid1d) = forward_map(g).($op(supergrid(g)))
end

resize(g::AbstractMappedGrid, n::Int) = map_grid(forward_map(g), resize(supergrid(g), n))

function rescale(g::AbstractGrid1d, a, b)
	m = mapto(covering(g), a..b)
	map_grid(m, g)
end

(→)(g::AbstractGrid1d, d::AbstractInterval) = rescale(g, infimum(d), supremum(d))

map_grid(grid::AbstractGrid, map) = (@warn "Interchange arguments of map_grid here"; map_grid(map, grid))
map_grid(map, grid::AbstractGrid) = map_grid1(map, grid)
map_grid1(map, grid::AbstractGrid) = map_grid2(map, grid)
map_grid2(map, grid) = MappedGrid(map, grid)

# some simplifications
map_grid1(map, g::AbstractMappedGrid) = map_grid(map∘forward_map(g), supergrid(g))
map_grid2(map::IdentityMap, grid) = grid

# Convenience function, similar to apply_map for Dictionary's
apply_map(grid::AbstractGrid, map::AbstractMap) = map_grid(map, grid)

forward_map(g::AbstractGrid, x...) = forward_map(g)(x...)
inverse_map(g::AbstractGrid, x...) = inverse_map(g)(x...)

unsafe_grid_getindex(g::AbstractMappedGrid, I...) = forward_map(g, supergrid(g, I...))

# Preserve tensor product structure
function rescale(g::ProductGrid, a::SVector{N}, b::SVector{N}) where {N}
	scaled_grids = [ rescale(component(g, i), a[i], b[i]) for i in 1:N]
	ProductGrid(scaled_grids...)
end

# Base.show(io::IO, mime::MIME"text/plain", g::AbstractMappedGrid) =
# 	composite_show(io, mime, g)
# Display.displaystencil(g::AbstractMappedGrid) =
#     DomainSets.map_stencil_broadcast(forward_map(g), supergrid(g))
# Display.object_parentheses(g::AbstractMappedGrid) =
#     Display.object_parentheses(forward_map(g))
# Display.stencil_parentheses(g::AbstractMappedGrid) =
#     Display.stencil_parentheses(forward_map(g))

"""
A MappedGrid consists of a grid and a map. Each grid point of the mapped grid
is the map of the corresponding point of the underlying grid.
"""
struct MappedGrid{T,N,M,G} <: AbstractMappedGrid{T,N}
	map		::	M
	grid	::	G
end

const MappedGrid1d{T<:Number,N,M,G} = MappedGrid{T,N,M,G}

MappedGrid(map, grid::AbstractGrid{T,N}) where {T,N} =
	MappedGrid{T,N}(map, grid)
MappedGrid(map::Map{S}, grid::AbstractGrid{T,N}) where {S,T,N} =
	MappedGrid{codomaintype(map),N}(map, grid)
MappedGrid{T,N}(map, grid) where {T,N} =
	MappedGrid{T,N,typeof(map),typeof(grid)}(map, grid)

MappedGrid(grid::AbstractGrid, map) = MappedGrid(map, grid)

forward_map(g::MappedGrid) = g.map
inverse_map(g::MappedGrid) = inverse(forward_map(g))
inverse_map(g::MappedGrid, x) = inverse(forward_map(g), x)


text:
# FusionTreeIterator:
# iterate over fusion trees for fixed incoming and outgoing sector labels
#==============================================================================#
struct FusionTreeIterator{G<:Sector,N}
    outgoing::NTuple{N,G}
    incoming::G
end

fusiontrees(outgoing::NTuple{N,G}, incoming::G = one(G)) where {N,G<:Sector} = FusionTreeIterator{G,N}(outgoing, incoming)

Base.iteratorsize(::FusionTreeIterator) = Base.SizeUnknown()
Base.iteratoreltype(::FusionTreeIterator) = Base.HasEltype()
Base.eltype(T::Type{FusionTreeIterator{G,N}}) where {G<:Sector, N} = _eltype(T, valsub(Val(N),Val(2)), valsub(Val(N),Val(1)), vertex_labeltype(G))
_eltype(::Type{FusionTreeIterator{G,N}}, ::Val{M}, ::Val{L}, ::Type{T}) where {G<:Sector,N,M,L,T} = FusionTree{G,N,M,L,T}

# * Iterator methods: start, next, done
#   Start with special cases:
Base.start(it::FusionTreeIterator{G,0} where {G<:Sector}) = false
function Base.next(it::FusionTreeIterator{G,0}, state) where {G<:Sector}
    T = vertex_labeltype(G)
    tree = FusionTree{G,0,0,0,T}((), one(G), (), ())
    return tree, true
end
Base.done(it::FusionTreeIterator{G,0} where {G<:Sector}, state) = state

Base.start(it::FusionTreeIterator{G,1}) where {G<:Sector} = it.outgoing[1] != it.incoming
function Base.next(it::FusionTreeIterator{G,1}, state) where {G<:Sector}
    T = vertex_labeltype(G)
    tree = FusionTree{G,1,0,0,T}(it.outgoing, it.incoming, (), ())
    return tree, true
end
Base.done(it::FusionTreeIterator{G,1} where {G<:Sector}, state) = state

#   General case:
Base.start(it::FusionTreeIterator{G,N} where {N}) where {G<:Sector} = _start(it.outgoing, it.incoming)
function Base.next(it::FusionTreeIterator{G,N} where {N}, state) where {G<:Sector}
    lines, vertices = _nextval(it.outgoing, it.incoming, state)
    f = FusionTree(it.outgoing, it.incoming, lines, vertices)
    return f, _nextstate(it.outgoing, it.incoming, state)
end
Base.done(it::FusionTreeIterator{G,N} where {N}, state) where {G<:Sector} = _done(it.outgoing, it.incoming, state)

# NOTE: Alternative code below specializes to optimize Abelian case, but in practice the
# speedup seems negligable
#
# Base.start(it::FusionTreeIterator{G,N₁,N₂} where {N₁,N₂}) where {G<:Sector} = _start(it, fusiontype(G))
# Base.next(it::FusionTreeIterator{G,N₁,N₂} where {N₁,N₂}, state) where {G<:Sector} = _next(it, fusiontype(G), state)
# Base.done(it::FusionTreeIterator{G,N₁,N₂} where {N₁,N₂}, state) where {G<:Sector} = _done(it, fusiontype(G), state)
# #   - Optimized version for Abelian case
# function _start(it::FusionTreeIterator,::Type{Abelian})
#     slines = _abeliantree(it.outgoing)
#     flines = _abeliantree(it.incoming)
#     match = last(slines) == last(flines)
#     return (!match, slines, flines)
# end
# function _next(it::FusionTreeIterator,::Type{Abelian}, state)
#     _, slines, flines = state
#     svertices = map(l->nothing, slines)
#     fvertices = map(l->nothing, flines)
#     tree = FusionTree(it.outgoing, it.incoming, slines, svertices, flines, fvertices)
#     return tree, (true, slines, flines)
# end
# _done(it::FusionTreeIterator,::Type{Abelian}, state) = state[1]
#
# @inline _abeliantree(incoming::NTuple{2,G}) where {G<:Sector} = (first(incoming[1] × incoming[2]),)
# @inline function _abeliantree(incoming::NTuple{N,G}) where {N,G<:Sector}
#     a,b, = incoming
#     c = first(a × b)
#     return (c, _abeliantree(tuple(c, tail2(incoming)...))...)
# end
# #   - General non-abelian case
# function _start(it::FusionTreeIterator,::Type{<:Fusion})
#     si = _start(it.incoming, nothing)
#     c = last(si) # last entry of state should be the outcoming line label
#     so = _start(it.outgoing, c)
#     while _done(it.outgoing, c, so)
#         si = _nextstate(it.incoming, nothing, si)
#         _done(it.incoming, nothing, si) && break
#         c = last(si)
#         so = _start(it.outgoing, c)
#     end
#     return (si, so)
# end
# function _next(it::FusionTreeIterator,::Type{<:Fusion}, state)
#     si, so = state
#     c = last(si)
#     # create tree
#     flines, fvertices = _nextval(it.incoming, nothing, si)
#     slines, svertices = _nextval(it.outgoing, c, so)
#     tree = FusionTree(it.outgoing, it.incoming, slines, svertices, flines, fvertices)
#     # create next state
#     c = last(si)
#     so = _nextstate(it.outgoing, c, so)
#     while _done(it.outgoing, c, so)
#         si = _nextstate(it.incoming, nothing, si)
#         _done(it.incoming, nothing, si) && break
#         c = last(si)
#         so = _start(it.outgoing, c)
#     end
#     return tree, (si, so)
# end
# _done(it::FusionTreeIterator,::Type{<:Fusion}, state) = _done(it.incoming, nothing, state[1])

# Actual implementation
function _start(outgoing::NTuple{2,G}, incoming::G) where {G<:Sector}
    a, b = outgoing
    return ((1,),)
end
_nextstate(outgoing::NTuple{2,G}, incoming::G, state) where {G<:Sector} = ((state[1][1]+1,),)
_nextval(outgoing::NTuple{2,G}, incoming::G, state) where {G<:Sector}  = (), (vertex_ind2label(state[1][1], outgoing[1], outgoing[2], incoming),)
_done(outgoing::NTuple{2,G}, incoming::G, state) where {G<:Sector} = state[1][1] > Nsymbol(outgoing[1],outgoing[2], incoming)

function _start(outgoing::NTuple{N,G}, incoming::G) where {N, G<:Sector}
    a, b, = outgoing
    it = a ⊗ b
    s = start(it) # done(it1,s1) == false: there should always be at least one fusion output
    c, snext = next(it, s)
    resttree = tuple(c, tail2(outgoing)...)
    reststate = _start(resttree, incoming)
    while _done(resttree, incoming, reststate)
        s = snext
        done(it, s) && break
        c, snext = next(it, s)
        resttree = tuple(c, tail2(outgoing)...)
        reststate = _start(resttree, incoming)
    end
    return tuple((1, s), reststate...)
end
function _nextstate(outgoing::NTuple{N,G}, incoming::G, state) where {N, G<:Sector}
    a,b, = outgoing
    it = a ⊗ b
    n, s = state[1]
    c, snext = next(it, s)
    resttree = tuple(c, tail2(outgoing)...)
    reststate = _nextstate(resttree, incoming, tail(state))
    while _done(resttree, incoming, reststate)
        if n < Nsymbol(a, b, c)
            n += 1
        else
            s = snext
            done(it, s) && break
            c, snext = next(it, s)
            resttree = tuple(c, tail2(outgoing)...)
        end
        reststate = _start(resttree, incoming)
    end
    return ((n, s), reststate...)
end
function _nextval(outgoing::NTuple{N,G}, incoming::G, state) where {N,G<:Sector}
    a, b, = outgoing
    it = a ⊗ b
    n, s = state[1]
    c, = next(it, s)
    resttree = tuple(c, tail2(outgoing)...)
    reststate = tail(state)
    lines, vertices = _nextval(resttree, incoming, reststate)
    return tuple(c, lines...), tuple(vertex_ind2label(n, outgoing[1], outgoing[2],c), vertices...)
end
_done(outgoing::NTuple{N,G}, incoming::G, state) where {N, G<:Sector} = done(outgoing[1] ⊗ outgoing[2], state[1][2])


text:
push!(LOAD_PATH, joinpath(pwd(), "src"))

using Neapolitan
using Base.Test

#

source = """
---
# Example
"""

result = Neapolitan.render(source)

@test result == "
Example
\n"

#

source = """
---
title: Example Title
---
# {{title}}
"""

result = Neapolitan.render(source)

@test result == "
Example Title
\n"


text:
using Flux.Optimise
using Flux.Optimise: runall
using Flux.Tracker
using Test
@testset "Optimise" begin
  w = randn(10, 10)
  @testset for opt in [ADAMW(), ADAGrad(0.1), AdaMax(), ADADelta(0.9), AMSGrad(),
                       NADAM(), Descent(0.1), ADAM(), Nesterov(), RMSProp(),
                       Momentum()]
    w′ = param(randn(10, 10))
    loss(x) = Flux.mse(w*x, w′*x)
    for t = 1: 10^5
      θ = Params([w′])
      θ̄ = gradient(() -> loss(rand(10)), θ)
      Optimise.update!(opt, θ, θ̄)
    end
    @test Flux.mse(w, w′) < 0.01
  end
end

@testset "Optimiser" begin
  w = randn(10, 10)
  @testset for Opt in [InvDecay, WeightDecay, ExpDecay]
    w′ = param(randn(10, 10))
    loss(x) = Flux.mse(w*x, w′*x)
    opt = Optimiser(Opt(), ADAM(0.001))
    for t = 1:10^5
      l = loss(rand(10))
      back!(l)
      delta = Optimise.apply!(opt, w′.data, w′.grad)
      w′.data .-= delta
    end
    @test Flux.mse(w, w′) < 0.01
  end
end

@testset "Training Loop" begin
  i = 0
  l = param(1)

  Flux.train!(() -> (sleep(0.1); i += 1; l),
              (),
              Iterators.repeated((), 100),
              Descent(),
              cb = Flux.throttle(() -> (i > 3 && Flux.stop()), 1))

  @test 3 < i < 50

  # Test multiple callbacks
  x = 0
  fs = [() -> (), () -> x = 1]
  cbs = runall(fs)
  cbs()
  @test x == 1
end

@testset "ExpDecay" begin
    w = randn(10, 10)
    o = ExpDecay(0.1, 0.1, 1000, 1e-4)
    w1 = param(randn(10,10))
    loss(x) = Flux.mse(w*x, w1*x)
    flag = 1
    decay_steps = []
    for t = 1:10^5
      l = loss(rand(10))
      back!(l)
      prev_eta = o.eta
      prev_grad = collect(w1.grad)
      delta = Optimise.apply!(o, w1.data, w1.grad)
      w1.data .-= delta
      new_eta = o.eta
      if new_eta != prev_eta
        push!(decay_steps, t)
      end
      array = fill(o.eta, size(prev_grad))
      if array .* prev_grad != delta
        flag = 0
      end
    end
    @test flag == 1
    # Test to check if decay happens at decay steps. Eta reaches clip value eventually.
    ground_truth = []
    for i in 1:11
      push!(ground_truth, 1000*i)  # Expected decay steps for this example.
    end
    @test decay_steps == ground_truth
    @test o.eta == o.clip
end


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

include("model_BGPCommunity.jl")
include("model_BgpServiceCommunityListResult.jl")
include("model_BgpServiceCommunityPropertiesFormat.jl")
include("model_Resource.jl")
include("model_BgpServiceCommunity.jl")


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

module ContainerServiceClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_ContainerServiceApi.jl")
include("api_ContainerServicesApi.jl")

# export models
export convert, ContainerService
export convert, ContainerServiceAgentPoolProfile
export convert, ContainerServiceCustomProfile
export convert, ContainerServiceDiagnosticsProfile
export convert, ContainerServiceLinuxProfile
export convert, ContainerServiceListResult
export convert, ContainerServiceMasterProfile
export convert, ContainerServiceOrchestratorProfile
export convert, ContainerServiceProperties
export convert, ContainerServiceServicePrincipalProfile
export convert, ContainerServiceSshConfiguration
export convert, ContainerServiceSshPublicKey
export convert, ContainerServiceVMDiagnostics
export convert, ContainerServiceWindowsProfile
export convert, Resource

# export operations
export convert, ContainerServiceApi, ContainerServicesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
"""
    sphere_mesh(ne, radius::FT, wtype::AbstractSphereWarp) where {FT<:AbstractFloat}

Generates a sphere mesh of a given `radius`, warp type and number of elements along each side
of a cubed surface. The sphere mesh is generated by applying the warping function to a cube panel mesh.
"""
function sphere_mesh(ne, radius::FT, wtype::AbstractSphereWarp) where {FT<:AbstractFloat}
    mesh = cube_panel_mesh(ne, FT)
    mesh.coordinates .-= FT(0.5)
    mesh.coordinates .*= (FT(2) * radius)
    for i = 1:mesh.nverts
        mesh.coordinates[i, :] .= cubed_sphere_warp(wtype, mesh.coordinates[i, :]...)
    end
    return mesh
end

"""
    cubed_sphere_warp(::EquiangularSphereWarp, a, b, c)

Given points `(a, b, c)` on the surface of a cube, warp the points out to a
spherical shell of radius `R` based on the equiangular gnomonic grid proposed by
[Ronchi1996]

Returns a tuple of (x1, x2, x3), where x1, x2, x3 are the corresponding coordinates on a sphere
of radius R.

The "Cubed Sphere": A New Method for the Solution of Partial Differential Equations in Spherical Geometry
C. RONCHI,*,1 R. IACONO,* AND P. S. PAOLUCCI†
JOURNAL OF COMPUTATIONAL PHYSICS 124, 93–114 (1996)
https://github.com/CliMA/ClimateMachine.jl/blob/master/src/Numerics/Mesh/Topologies.jl
"""
function cubed_sphere_warp(wtype::EquiangularSphereWarp, a, b, c)
    R = max(abs(a), abs(b), abs(c))
    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        ζ1 = sR / sqrt(X^2 + Y^2 + 1)
        ζ2, ζ3 = X * ζ1, Y * ζ1
        ζ1, ζ2, ζ3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1
        if a < 0
            # (-R, *, *) : formulas for Face I from Ronchi, Iacono, Paolucci (1996)
            #              but for us face IV of the developed net of the cube
            x1, x2, x3 = f(-R, b / a, c / a)
        elseif a > 0
            # ( R, *, *) : formulas for Face III from Ronchi, Iacono, Paolucci (1996)
            #              but for us face II of the developed net of the cube
            x1, x2, x3 = f(R, b / a, c / a)
        end
    elseif fdim == 2
        if b < 0
            # ( *,-R, *) : formulas for Face II from Ronchi, Iacono, Paolucci (1996)
            #              but for us face V of the developed net of the cube
            x2, x1, x3 = f(-R, a / b, c / b)
        elseif b > 0
            # ( *, R, *) : formulas for Face IV from Ronchi, Iacono, Paolucci (1996)
            #              but for us face III of the developed net of the cube
            x2, x1, x3 = f(R, a / b, c / b)
        end
    elseif fdim == 3
        if c > 0
            # ( *, *, R) : formulas for Face V from Ronchi, Iacono, Paolucci (1996)
            #              but for us face VI of the developed net of the cube
            x3, x2, x1 = f(R, b / c, a / c)
        elseif c < 0
            # ( *, *,-R) : formulas for Face VI from Ronchi, Iacono, Paolucci (1996)
            #              but for us face I of the developed net of the cube
            x3, x2, x1 = f(-R, b / c, a / c)
        end
    end

    return x1, x2, x3
end

"""
    cubed_sphere_warp(::EquidistantSphereWarp, a, b, c)

Returns a tuple of (x1, x2, x3), where x1, x2, x3 are the corresponding coordinates on a sphere
of radius R.

https://github.com/CliMA/ClimateMachine.jl/blob/master/src/Numerics/Mesh/Topologies.jl
"""
function cubed_sphere_warp(wtype::EquidistantSphereWarp, a, b, c)
    R = max(abs(a), abs(b), abs(c))
    r = hypot(a, b, c)

    return R * a / r, R * b / r, R * c / r
end

"""
    cube_panel_mesh(ne)

This function builds a cube panel mesh with a resolution of `ne` elements along each edge.

               v8 (xs,xe,xe)          v7 (xe,xe,xe)
                 o--------e11---------o
                /|                   /|
               / |                  / |
              /  |                 /  |
            e12  e8               e10 e7
            /    |               /    |
           /     |            v6/     |
       v5 o--------e9----------o      |
          |    v4o------e3-----|------o v3 (xe,xe,xs)
          |     /   o------->  |     /
          |    /   /           |    /
          e5  e4  /           e6   e2
          |  /   /             |  /
          | /   /              | /
          |/   o               |/
          o--------e1----------o
         v1                    v2
       (xs,xs,xs)               (xe,xs,xs)

       panel 1 => 1 4 3 2
       panel 2 => 2 3 7 6
       panel 3 => 3 4 8 7
       panel 4 => 1 5 8 4
       panel 5 => 1 2 6 5
       panel 6 => 5 6 7 8

       edge  1 => 1 2
       edge  2 => 2 3
       edge  3 => 3 4
       edge  4 => 4 1
       edge  5 => 1 5
       edge  6 => 2 6
       edge  7 => 3 7
       edge  8 => 4 8
       edge  9 => 5 6
       edge 10 => 6 7
       edge 11 => 7 8
       edge 12 => 8 5

            v8 +---e11---+ v7
               | ^       |
               | |       |
              e8 ^   3   e7
               | |       |
            v4 | o--<-o  | v3      v7        v8
  v8 +---e8----+---e3----+----e7---+---e11---+
     | ^       | o-->->  | o-->--> | ^       |
     | |       | |       | |       | |       |
   e12 ^  4   e4 ^  1   e2 ^  2   e10^  6   e12
     | |       | |       | |       | |       |
     | o-<-<-o | o       | o       | o<-<--o |
     +---e5----+---e1----+----e6---+----e9---+
    v5       v1| o-->--o | v2      v6        v5
               |       | |
              e5    5  | e6
               |       ↓ |
               |         |
               +---e9----+
              v5         v6
"""
function cube_panel_mesh(ne::I, ::Type{FT}) where {FT<:AbstractFloat,I<:Integer}
    xs, xe = FT(0), FT(1)
    nverts = (ne + 1)^3 - (ne - 1)^3
    nfaces = 12 * ne + 6 * (2 * ne * (ne - 1))
    nelems = 6 * ne * ne
    nbndry = 0

    nx = ne + 1

    nfaces_edg = 12 * ne

    emat = reshape(1:nelems, ne, ne, 6)
    ndmat = zeros(I, ne + 1, ne + 1)

    panel_verts = [
        1 2 3 1 1 5
        4 3 4 5 2 6
        3 7 8 8 6 7
        2 6 7 4 5 8
    ]
    panel_edges = [
        4 2 3 5 1 9
        3 7 8 12 6 10
        2 10 11 8 9 11
        1 6 7 4 5 12
    ]
    panel_edges_rev = Bool.([
        1 0 0 0 0 0
        1 0 0 1 0 0
        1 1 1 1 1 0
        1 1 1 0 1 0
    ])
    # node coordinates
    xc = range(xs, xe; step = FT(1 / ne)) # [xs, xs+Δ, xs+2Δ, ..., xe]
    xci = view(xc, 2:ne)                   # [xs+Δ, xs+2Δ, ..., xe-Δ]
    xcir = view(xc, ne:-1:2)                # reverse(xci)
    sc = ones(FT, ne - 1) * xs              # [xs, ...., xs]
    ec = ones(FT, ne - 1) * xe              # [xe, ...., xe]
    # where Δ = (xe-xs)/ne
    sc2 = ones(FT, (ne - 1) * (ne - 1)) * xs
    ec2 = ones(FT, (ne - 1) * (ne - 1)) * xe

    xci12 = repeat(xci, outer = ne - 1)
    xci1r2 = repeat(xcir, outer = ne - 1)
    xci21 = repeat(xci, inner = ne - 1)

    edge_nodes = reshape(1:((ne-1)*12), ne - 1, 12) .+ 8
    # coordinates
    coordinates = vcat(
        hcat(
            [xs, xe, xe, xs, xs, xe, xe, xs], # x1,
            [xs, xs, xe, xe, xs, xs, xe, xe], # x2,
            [xs, xs, xs, xs, xe, xe, xe, xe], # x3 vertex coordinates
        ),
        vcat(
            hcat(xci, sc, sc), # edge 1
            hcat(ec, xci, sc), # edge 2
            hcat(xcir, ec, sc), # edge 3
            hcat(sc, xcir, sc), # edge 4
            hcat(sc, sc, xci), # edge 5
            hcat(ec, sc, xci), # edge 6
            hcat(ec, ec, xci), # edge 7
            hcat(sc, ec, xci), # edge 8
            hcat(xci, sc, ec), # edge 9
            hcat(ec, xci, ec), # edge 10
            hcat(xcir, ec, ec), # edge 11
            hcat(sc, xcir, ec), # edge 12
        ),
        hcat(xci21, xci12, sc2), # panel 1
        hcat(ec2, xci12, xci21), # panel 2
        hcat(xci1r2, ec2, xci21), # panel 3
        hcat(sc2, xci21, xci12), # panel 4
        hcat(xci12, sc2, xci21), # panel 5
        hcat(xci12, xci21, ec2), # panel 6
    )
    face_interior = reshape(1:((ne-1)*(ne-1)), ne - 1, ne - 1)

    nfc1i = (nx - 2) * (nx - 1) # panels with normals along local first direction
    nfc2i = (nx - 1) * (nx - 2) # panels with normals along local second direction
    nfci = nfc1i + nfc2i
    fci1 = reshape(1:nfc1i, nx - 2, nx - 1)
    fci2 = reshape(1:nfc2i, nx - 1, nx - 2)

    fcmat1 = zeros(I, nx, nx - 1) # face numbering
    fcmat2 = zeros(I, nx - 1, nx) # for each panel

    edge_faces = reshape(1:(12*ne), ne, 12)

    face_verts = zeros(I, nfaces, 2)
    face_neighbors = zeros(I, nfaces, 5)
    face_boundary = Vector{I}(zeros(nfaces)) # all interior nodes (no boundaries)
    elem_verts = zeros(I, nelems, 4)
    elem_faces = zeros(I, nelems, 4)

    for sfc = 1:6
        ndmat[1, 1],
        ndmat[ne+1, 1],  # panel vertices
        ndmat[ne+1, ne+1],
        ndmat[1, ne+1] = panel_verts[:, sfc]

        ndmat[2:ne, 1] .=
            panel_edges_rev[1, sfc] ? reverse(edge_nodes[:, panel_edges[1, sfc]]) :
            edge_nodes[:, panel_edges[1, sfc]]  # panel edges

        ndmat[end, 2:ne] .=
            panel_edges_rev[2, sfc] ? reverse(edge_nodes[:, panel_edges[2, sfc]]) :
            edge_nodes[:, panel_edges[2, sfc]]

        ndmat[ne:-1:2, end] .=
            panel_edges_rev[3, sfc] ? reverse(edge_nodes[:, panel_edges[3, sfc]]) :
            edge_nodes[:, panel_edges[3, sfc]]

        ndmat[1, ne:-1:2, end] .=
            panel_edges_rev[4, sfc] ? reverse(edge_nodes[:, panel_edges[4, sfc]]) :
            edge_nodes[:, panel_edges[4, sfc]]

        offset = 8 + 12 * (ne - 1) + (sfc - 1) * (ne - 1) * (ne - 1) # interior
        ndmat[2:ne, 2:ne] .= face_interior .+ offset

        fcmat1[1, end:-1:1] .=
            panel_edges_rev[4, sfc] ? reverse(edge_faces[:, panel_edges[4, sfc]]) :
            edge_faces[:, panel_edges[4, sfc]]

        fcmat1[end, :] .=
            panel_edges_rev[2, sfc] ? reverse(edge_faces[:, panel_edges[2, sfc]]) :
            edge_faces[:, panel_edges[2, sfc]]

        fcmat2[:, 1] .=
            panel_edges_rev[1, sfc] ? reverse(edge_faces[:, panel_edges[1, sfc]]) :
            edge_faces[:, panel_edges[1, sfc]]

        fcmat2[end:-1:1, end] .=
            panel_edges_rev[3, sfc] ? reverse(edge_faces[:, panel_edges[3, sfc]]) :
            edge_faces[:, panel_edges[3, sfc]]

        off = ne * 12 + (sfc - 1) * nfci

        fcmat1[2:(end-1), :] .= fci1 .+ off
        fcmat2[:, 2:(end-1)] .= fci2 .+ (off + nfc1i)

        face_verts[fcmat1[:], 1] .= ndmat[:, 1:ne][:] # face nodes
        face_verts[fcmat1[:], 2] .= ndmat[:, 2:(ne+1)][:]
        face_verts[fcmat2[:], 1] .= ndmat[1:ne, :][:]
        face_verts[fcmat2[:], 2] .= ndmat[2:(ne+1), :][:]

        if sfc == 1
            bdy1 = emat[:, 1:1, 5]'
            bdy2 = emat[ne:-1:1, 1:1, 3]'
            bdy3 = emat[1:1, 1:ne, 4]'
            bdy4 = emat[:, 1:1, 2]
        elseif sfc == 2
            bdy1 = emat[ne:ne, :, 5]
            bdy2 = emat[1:1, :, 3]
            bdy3 = emat[:, ne:ne, 1]
            bdy4 = emat[ne:ne, :, 6]'
        elseif sfc == 3
            bdy1 = emat[ne:ne, :, 2]
            bdy2 = emat[:, ne:ne, 4]'
            bdy3 = emat[ne:ne, ne:-1:1, 1]'
            bdy4 = emat[ne:-1:1, ne:ne, 6]
        elseif sfc == 4
            bdy1 = emat[:, 1:1, 1]'
            bdy2 = emat[1:1, :, 6]
            bdy3 = emat[1:1, :, 5]'
            bdy4 = emat[ne:ne, :, 3]'
        elseif sfc == 5
            bdy1 = emat[:, 1:1, 4]'
            bdy2 = emat[1:1, :, 2]
            bdy3 = emat[1:1, :, 1]'
            bdy4 = emat[:, 1:1, 6]
        else # sfc == 6
            bdy1 = emat[ne:ne, :, 4]
            bdy2 = emat[:, ne:ne, 2]'
            bdy3 = emat[:, ne:ne, 5]
            bdy4 = emat[ne:-1:1, ne:ne, 3]
        end
        face_neighbors[fcmat1[:], 1] .= vcat(bdy1, emat[:, :, sfc])[:]
        face_neighbors[fcmat1[:], 3] .= vcat(emat[:, :, sfc], bdy2)[:]
        face_neighbors[fcmat2[:], 1] .= hcat(bdy3, emat[:, :, sfc])[:]
        face_neighbors[fcmat2[:], 3] .= hcat(emat[:, :, sfc], bdy4)[:]

        elem_verts[emat[:, :, sfc][:], :] .= hcat(
            ndmat[1:ne, 1:ne][:], # node numbers (local node 1)
            ndmat[2:(ne+1), 1:ne][:], # for each element (local node 2)
            ndmat[2:(ne+1), 2:(ne+1)][:], #(local node 3)
            ndmat[1:ne, 2:(ne+1)][:], # (local node 4)
        )

        elem_faces[emat[:, :, sfc][:], :] .= hcat(
            fcmat2[:, 1:(nx-1)][:], # (local face 1)
            fcmat1[2:nx, :][:],   # each element (local face 2)
            fcmat2[:, 2:nx][:], # (local face 3)
            fcmat1[1:(nx-1), :][:], # face numbers for (local face 4)
        )
    end

    ref_fc_verts = [
        1 2 3 4
        2 3 4 1
    ]

    for fc = 1:nfaces
        elems = (face_neighbors[fc, 1], face_neighbors[fc, 3])
        for e = 1:2
            el = elems[e]
            localface = findfirst(elem_faces[el, :] .== fc)
            face_neighbors[fc, 2+(e-1)*2] = localface
        end
        # setting up relative orientation
        or1 =
            elem_verts[elems[1], ref_fc_verts[:, face_neighbors[fc, 2]]] ==
            face_verts[fc, :] ? 1 : -1
        or2 =
            elem_verts[elems[2], ref_fc_verts[:, face_neighbors[fc, 4]]] ==
            face_verts[fc, :] ? 1 : -1
        face_neighbors[fc, 5] = or1 * or2
    end

    boundary_tags = sort(unique(face_boundary))
    face_boundary_tags = sort(face_boundary)
    face_boundary = sortperm(face_boundary)

    face_boundary_offset = I[1, length(face_boundary)+1]

    boundary_tag_names = (:interior,)

    # add unique vertex iterator information
    vtconn = map(i -> zeros(I, i), zeros(I, nverts))

    for el = 1:nelems
        for lv = 1:4
            vt = elem_verts[el, lv]
            push!(vtconn[vt], el, lv)
        end
    end
    unique_verts = I[]
    uverts_conn = I[]
    uverts_offset = I.([1])
    for vt = 1:nverts
        lconn = length(vtconn[vt])
        if lconn > 0
            push!(unique_verts, vt)
            push!(uverts_conn, vtconn[vt]...)
            push!(uverts_offset, uverts_offset[end] + lconn)
        end
    end

    return Mesh2D(
        nverts,
        nfaces,
        nelems,
        nbndry,
        coordinates,
        unique_verts,
        uverts_conn,
        uverts_offset,
        face_verts,
        face_neighbors,
        face_boundary,
        boundary_tags,
        boundary_tag_names,
        face_boundary_offset,
        elem_verts,
        elem_faces,
    )
end

function unfold_cube_panel_to_plane(
    x::FTA1D,
    y::FTA1D,
    z::FTA1D,
    panelno::Int,
) where {FT<:AbstractFloat,FTA1D<:AbstractVector{FT}}
    @assert 1 ≤ panelno ≤ 6
    n = length(x)
    xx = similar(x)
    yy = similar(x)

    if panelno == 1
        xx = x
        yy = y
    elseif panelno == 2
        xx = z .+ FT(1)
        yy = y
    elseif panelno == 3
        xx = x
        yy = z .+ FT(1)
    elseif panelno == 4
        xx = -z
        yy = y
    elseif panelno == 5
        xx = x
        yy = -z
    elseif panelno == 6
        xx = 3 .- x
        yy = y
    end

    return xx, yy
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
module LRP_CONFIG
using ExplainableAI
using ExplainableAI: LRPSupportedLayer, LRPSupportedActivation
"""
    LRP_CONFIG.supports_layer(layer)

Check whether LRP can be used on a layer or a Chain.
To extend LRP to your own layers, define:
```julia
LRP_CONFIG.supports_layer(::MyLayer) = true          # for structs
LRP_CONFIG.supports_layer(::typeof(mylayer)) = true  # for functions
```
"""
supports_layer(l) = false
supports_layer(::LRPSupportedLayer) = true
"""
    LRP_CONFIG.supports_activation(σ)

Check whether LRP can be used on a given activation function.
To extend LRP to your own activation functions, define:
```julia
LRP_CONFIG.supports_activation(::typeof(myactivation)) = true  # for functions
LRP_CONFIG.supports_activation(::MyActivation) = true          # for structs
```
"""
supports_activation(σ) = false
supports_activation(::LRPSupportedActivation) = true
end # LRP_CONFIG module

_check_layer(::Val{:LRP}, layer) = LRP_CONFIG.supports_layer(layer)
_check_layer(::Val{:LRP}, c::Chain) = all(_check_layer(Val(:LRP), l) for l in c)

function _check_activation(::Val{:LRP}, layer)
    hasproperty(layer, :σ) && return LRP_CONFIG.supports_activation(layer.σ)
    return true
end
_check_activation(::Val{:LRP}, c::Chain) = all(_check_activation(Val(:LRP), l) for l in c)

"""
    check_model(method::Symbol, model; verbose=true)

Check whether the given method can be used on the model.
Currently, model checks are only implemented for LRP, using the symbol `:LRP`.

# Example
julia> check_model(:LRP, model)
"""
check_model(method::Symbol, model; kwargs...) = check_model(Val(method), model; kwargs...)
function check_model(::Val{:LRP}, c::Chain; verbose=true)
    layer_checks = collect(_check_layer.(Val(:LRP), c.layers))
    activation_checks = collect(_check_activation.(Val(:LRP), c.layers))
    passed_layer_checks = all(layer_checks)
    passed_activation_checks = all(activation_checks)

    passed_layer_checks && passed_activation_checks && return true

    layer_names = [_print_name(l) for l in c]
    activation_names = [_print_activation(l) for l in c]

    verbose && _show_check_summary(
        c, layer_names, layer_checks, activation_names, activation_checks
    )
    if !passed_layer_checks
        verbose && display(
            Markdown.parse(
                """# Layers failed model check
                Found unknown layers `$(join(unique(layer_names[.!layer_checks]), ", "))`
                that are not supported by ExplainableAI's LRP implementation yet.

                If you think the missing layer should be supported by default, please [submit an issue](https://github.com/adrhill/ExplainableAI.jl/issues).

                These model checks can be skipped at your own risk by setting the LRP-analyzer keyword argument `skip_checks=true`.

                ## Using custom layers
                If you implemented custom layers, register them via
                ```julia
                LRP_CONFIG.supports_layer(::MyLayer) = true               # for structs
                LRP_CONFIG.supports_layer(::typeof(mylayer)) = true  # for functions
                ```
                The default fallback for this layer will use Automatic Differentiation according to "Layer-Wise Relevance Propagation: An Overview".
                """,
            ),
        )
        throw(ArgumentError("Unknown layers found in model"))
    end
    if !passed_activation_checks
        verbose && display(
            Markdown.parse(
                """ # Activations failed model check
                Found layers with unknown  or unsupported activation functions
                `$(join(unique(activation_names[.!activation_checks]), ", "))`.
                LRP assumes that the model is a "deep rectifier network" that only contains ReLU-like activation functions.

                If you think the missing activation function should be supported by default, please [submit an issue](https://github.com/adrhill/ExplainableAI.jl/issues).

                These model checks can be skipped at your own risk by setting the LRP-analyzer keyword argument `skip_checks=true`.

                ## Using custom activation functions
                If you use custom ReLU-like activation functions, register them via
                ```julia
                LRP_CONFIG.supports_activation(::typeof(myfunction)) = true  # for functions
                LRP_CONFIG.supports_activation(::MyActivation) = true        # for structs
                ```
                """,
            ),
        )
        throw(ArgumentError("Unknown or unsupported activation functions found in model"))
    end
end


text:
function makeit()
    n1 = SVector(1,0,0.0);
    n2 = SVector(0,1,0.0);
    n3 = SVector(0,0,1.0);

    pl1 = ImplicitPlane(n1, n1)
    pl2 = ImplicitPlane(n2, n2)
    pl3 = ImplicitPlane(n3, n3)
    pl4 = ImplicitPlane(-n1, -n1)
    pl5 = ImplicitPlane(-n2, -n2)
    pl6= ImplicitPlane(-n3, -n3)
    mpl = ImplicitPlane([4,4,0.25], -n3)
    cyl = ImplicitCylinder([0,0,1], [0,0,0], .25)
    return [pl1, pl2, pl3, pl4, pl5, pl6, cyl, mpl]
end

surfac = makeit();
Random.seed!(1234);
pontok = [rand(3)*20*rand() for i in 1:100];
cnodes, cvals, cnorms = cachenodes(surfac, pontok)

Random.seed!(1234);
trc1 = randomcachedtree(cnodes, 12);
trc2 = randomcachedtree(cnodes, 1);
randomcachedtree(cnodes, 1);
trc3 = randomcachedtree(cnodes, 1);
trc4 = randomcachedtree(cnodes, 2);
randomcachedtree(cnodes, 4);
trc5 = randomcachedtree(cnodes, 4);
trc6 = CachedCSGNode(:complement, [cnodes[1]], gensym());
trc7 = cnodes[7];
trc8 = CachedCSGNode(:complement, [cnodes[7]], gensym());
trc9 = CachedCSGNode(:subtraction, [cnodes[1], cnodes[7]], gensym());
trc10 = randomcachedtree(cnodes, 4);
trc11 = CachedCSGNode(:subtraction, [cnodes[1], cnodes[1]], gensym());

tr1 = cached2normaltree(trc1, surfac);
tr2 = cached2normaltree(trc2, surfac);
tr3 = cached2normaltree(trc3, surfac);
tr4 = cached2normaltree(trc4, surfac);
tr5 = cached2normaltree(trc5, surfac);
tr6 = cached2normaltree(trc6, surfac);
tr7 = cached2normaltree(trc7, surfac);
tr8 = cached2normaltree(trc8, surfac);
tr9 = cached2normaltree(trc9, surfac);
tr10 = cached2normaltree(trc10, surfac);
tr11 = cached2normaltree(trc11, surfac);

f1 = tree2func(trc1)
f2 = tree2func(trc2)
f3 = tree2func(trc3)
f4 = tree2func(trc4)
f5 = tree2func(trc5)
f6 = tree2func(trc6)
f7 = tree2func(trc7)
f8 = tree2func(trc8)
f9 = tree2func(trc9)
f10 = tree2func(trc10)
f11 = tree2func(trc11)

@testset "t1" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc1, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr1, p)

        fv1 = Base.invokelatest(f1, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        if vc1 != vf1
            @test fv1.index == 7
        end

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t2" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc2, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr2, p)

        fv1 = Base.invokelatest(f2, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t3" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc3, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr3, p)

        fv1 = Base.invokelatest(f3, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t4" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc4, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr4, p)

        fv1 = Base.invokelatest(f4, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t5" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc5, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr5, p)

        fv1 = Base.invokelatest(f5, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t6" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc6, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr6, p)

        fv1 = Base.invokelatest(f6, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t7" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc7, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr7, p)

        fv1 = Base.invokelatest(f7, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t8" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc8, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr8, p)

        fv1 = Base.invokelatest(f8, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t9" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc9, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr9, p)

        fv1 = Base.invokelatest(f9, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t10" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc10, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr10, p)

        fv1 = Base.invokelatest(f10, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end

@testset "t11" begin
    for i in eachindex(pontok)
        p = pontok[i]
        vc1, nc1 = valueandnormal(trc11, cvals, cnorms, i)
        v1, n1 = valueandnormal(tr11, p)

        fv1 = Base.invokelatest(f11, cvals, i)::CachedResult
        vf1 = value(fv1)
        nf1 = normal(fv1, cnorms, i)

        @test vc1 == v1 == vf1
        @test nc1 == n1 == nf1
    end
end


text:
"""
    printsubtypes(T::Type, level::Int = 0, indent::Int = 0)

Prints tree of subtypes of given type T.
"""
function printsubtypes(T::Type, level::Int = 0, indent::Int = 0)
    colors = [:light_magenta; :light_blue; :light_green; :light_yellow; :white]
    name = String(nameof(T))
    printstyled(repeat("   ", indent), name,"\n"; color=colors[level+1])
    printsubtypes.(subtypes(T), level + 1, indent + 1)
    return
end

"""
    printproblemtypes(P::Type)

Prints names of tabular dataset that are of given type (Classification or Regression).
"""
function printproblemtypes(P::Type)
    datasets = [subtypes(Tabular); subtypes(GrayImage); subtypes(ColorImage)]
    printstyled(String(nameof(P))," datatsets\n", bold=true)

    for i in datasets
        if problem(i()) == P
            printstyled(String(nameof(i))," "; color=:light_yellow)
            println("(",(nameof(supertype(typeof(i())))),")")
        end
    end
end

"""
    listdatasets()

Print datatsets in a tree structure. Prints all datasets of the PreprocessData package.
"""
listdatasets() = printsubtypes(DatasetName)

"""
    listdatasets(type)

Print datatsets in a tree structure. Prints all datasets of given type known to PreprocessData
package. Type can be any subtype of `DatasetName` or `Classifacation` or `Regression`.
"""
listdatasets(type::Type{T} where T <: DatasetName) = printsubtypes(type,  1)

listdatasets(type::Type{T} where T <: Problem) = printproblemtypes(type)

"""
    info(dataset::DatasetName)

Prints information about given dataset.

Dataset must have a registration file.
"""
function info(dataset::Tabular)
    printstyled("\n$dataset\n"; bold=true, color=:light_yellow)
    text = """
        Target column:  $(target(dataset))
    """
    println(infotext(dataset),text)
end

function info(dataset::Image{MLImage})
    printstyled("\n$dataset\n"; bold=true, color=:light_yellow)
    println(infotext(dataset))
end

"""
    infotext(dataset::DatasetName)

Returns String with basic info about given dataset.
"""
function infotext(dataset::DatasetName)
    text = """
        Name:           $(name(dataset))
        Type:           $(nameof(supertype(typeof(dataset))))
        Downloaded:     $(isdownloaded(dataset) ? "Yes $(DataDeps.determine_save_path((name(dataset))))" : "No")
        Source:         $(url(dataset))
        Size:           $(PreprocessData.size(dataset)[1]) (train data)
                        $(PreprocessData.size(dataset)[2]) (valid data)
                        $(PreprocessData.size(dataset)[3]) (test data)
        Problem type:   $(nameof(problem(dataset)))
    """
    msg = message(dataset)
    if !isempty(msg)
        text = text*"""
            Message:        $msg
        """
    end
    return text
end

"""
    remove(dataset::DatasetName)

Removes dataset directory.
"""
function remove(dataset::DatasetName)
    rm(@datadep_str name(dataset); recursive=true)
end

"""
    isdownloaded(dataset::DatasetName)

Returns true if file data-train.csv is present in dataset directory.
"""
function isdownloaded(dataset::DatasetName)
    path = joinpath(download_dir(),name(dataset))
    return isdir(path)
end

# Variation on function determine_save_path(name, rel=nothing) from DataDeps.jl (file
# locations.jl). Here are no arguments, returns path to directory where DataDeps files are
# downloaded (usually something like "/home/user/.julia/datadeps")
"""
    download_dir()

Determines the location to save all datadeps. Same as function from `DataDeps.jl` (see the
original function [`DataDeps.determine_save_path`](@ref).)
"""
function download_dir()
    rel=nothing;
    cands = DataDeps.preferred_paths(rel; use_package_dir=false)
    path_ind = findfirst(cands) do path
        0 == first(DataDeps.uv_access(path, DataDeps.W_OK))
    end
    if path_ind === nothing
            @error """
            No writable path exists to save the data. Make sure there exists as writable path in your DataDeps Load Path.
            See https://www.oxinabox.net/DataDeps.jl/stable/z10-for-end-users/#The-Load-Path-1
            The current load path contains:
            """ cands
        throw(DataDeps.NoValidPathError("No writable path exists to save the data."))
    end
    return String(cands[path_ind])
end

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
module HullGMP

using LazySets, LazySets.Approximations
using Polyhedra, CDDLib

using LinearAlgebra
using Parameters
using Interpolations # only for PiecewiseLinear

import LazySets: dim, HalfSpace # necessary to avoid conflict with Polyhedra

using Requires

using Base.Threads

include("activation.jl")
include("network.jl")
include("problem.jl")
include("util.jl")

function __init__()
  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("flux.jl")
end

export
    Solver,
    Network,
    AbstractActivation,
    #PolytopeComplement,
    #complement,
    # NOTE: not sure if exporting these is a good idea as far as namespace conflicts go:
    # ReLU,
    # Max,
    # Id,
    GeneralAct,
    PiecewiseLinear,
    LinearPieces,
    PiecewiseLinearActivation,
    Problem,
    Result,
    BasicResult,
    CounterExampleResult,
    AdversarialResult,
    ReachabilityResult,
    read_nnet,
    solve,
    forward_network,
    check_inclusion

export solve

include("reachability.jl")

include("maxSens.jl")
include("maxSens1.jl")
include("maxSens2.jl")

include("hullReach.jl")
include("hullReach1.jl")
include("hullReach2.jl")
include("hullReach3.jl")

include("speGuid.jl")
include("speGuid1.jl")
include("speGuid2.jl")

include("hullSearch.jl")
include("hullSearch1.jl")
include("hullSearch2.jl")

export MaxSens, MaxSens1, MaxSens2
export HullReach, HullReach1, HullReach2, HullReach3
export SpeGuid, SpeGuid1, SpeGuid2
export HullSearch, HullSearch1, HullSearch2

end # module


text:
module UnPack

export @unpack, @pack!

###########################
# Packing and unpacking @unpack, @pack!
##########################
# Below code slightly adapted from Simon Danisch's GLVisualize via PR
# https://github.com/mauro3/Parameters.jl/pull/13
#
# Note, this used to be part of the package Parameters.jl

"""
This function is invoked to unpack one field/entry of some DataType
`dt` and has signature:

`unpack(dt::Any, ::Val{sym}) -> value of sym`

The `sym` is the symbol of the assigned variable.

Three definitions are included in the package to unpack a composite type
or a dictionary with Symbol or string keys:
```
@inline unpack(x, ::Val{f}) where {f} = getproperty(x, f)
@inline unpack(x::AbstractDict{Symbol}, ::Val{k}) where {k} = x[k]
@inline unpack(x::AbstractDict{S}, ::Val{k}) where {S<:AbstractString,k} = x[string(k)]
```

More methods can be added to allow for specialized unpacking of other datatypes.

See also `pack!`.
"""
function unpack end
@inline unpack(x, ::Val{f}) where {f} = getproperty(x, f)
@inline unpack(x::AbstractDict{Symbol}, ::Val{k}) where {k} = x[k]
@inline unpack(x::AbstractDict{<:AbstractString}, ::Val{k}) where {k} = x[string(k)]

"""
This function is invoked to pack one entity into some DataType and has
signature:

`pack!(dt::Any, ::Val{property}, value) -> value`

Two definitions are included in the package to pack into a composite
type or into a dictionary with Symbol or string keys:

```
@inline pack!(x, ::Val{f}, val) where {f} = setproperty!(x, f, val)
@inline pack!(x::AbstractDict{Symbol}, ::Val{k}, val) where {k} = x[k]=val
@inline pack!(x::AbstractDict{S}, ::Val{k}, val) where {S<:AbstractString,k} = x[string(k)]=val
```

More methods can be added to allow for specialized packing of other
datatypes.

See also `unpack`.

To "pack" immutable datatypes, see the package Setfield.jl
"""
function pack! end
@inline pack!(x, ::Val{f}, val) where {f} = setproperty!(x, f, val)
@inline pack!(x::AbstractDict{Symbol}, ::Val{k}, val) where {k} = x[k]=val
@inline pack!(x::AbstractDict{<:AbstractString}, ::Val{k}, val) where {k} = x[string(k)]=val

"""
Unpacks fields/properties/keys from a composite type, a `Dict{Symbol}`, a `Dict{String}`,
or a module into variables
```julia_skip
@unpack a, b, c = dict_or_typeinstance
```

Example with dict:
```julia
d = Dict{Symbol,Any}(:a=>5.0,:b=>2,:c=>"Hi!")
@unpack a, c = d
a == 5.0 #true
c == "Hi!" #true
```

Example with type:
```julia
struct A; a; b; c; end
d = A(4,7.0,"Hi")
@unpack a, c = d
a == 4 #true
c == "Hi" #true
```

Note that its functionality can be extende by adding methods to the
`UnPack.unpack` function.
"""
macro unpack(args)
    args.head!=:(=) && error("Expression needs to be of form `a, b = c`")
    items, suitecase = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitecase_instance = gensym()
    kd = [:( $key = $UnPack.unpack($suitecase_instance, Val{$(Expr(:quote, key))}()) ) for key in items]
    kdblock = Expr(:block, kd...)
    expr = quote
        $suitecase_instance = $suitecase # handles if suitecase is not a variable but an expression
        $kdblock
        $suitecase_instance # return RHS of `=` as standard in Julia
    end
    esc(expr)
end

"""
Packs variables into a mutable, composite type, a `Dict{Symbol}`, or a `Dict{String}`
```julia_skip
@pack! dict_or_typeinstance = a, b, c
```

Example with dict:
```julia
a = 5.0
c = "Hi!"
d = Dict{Symbol,Any}()
@pack! d = a, c
d # Dict{Symbol,Any}(:a=>5.0,:c=>"Hi!")
```

Example with type:
```julia
a = 99
c = "HaHa"
mutable struct A; a; b; c; end
d = A(4,7.0,"Hi")
@pack! d = a, c
d.a == 99 #true
d.c == "HaHa" #true
```

Note that its functionality can be extende by adding methods to the
`UnPack.pack!` function.

To "pack" immutables use the package Setfield.jl.
"""
macro pack!(args)
    esc(_pack_bang(args))
end

function _pack_bang(args)
    args.head!=:(=) && error("Expression needs to be in the form of an assignment.")
    suitecase, items = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitecase_instance = gensym()
    kd = [:( $UnPack.pack!($suitecase_instance, Val{$(Expr(:quote, key))}(), $key) ) for key in items]
    kdblock = Expr(:block, kd...)
    return quote
        $suitecase_instance = $suitecase # handles if suitecase is not a variable but an expression
        $kdblock
        ($(items...),)
    end
end

end # module


text:
using BibUtils
using Test

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


text:
@testset "Viewing" begin
  @testset "Domain" begin
    g = CartesianGrid{T}(10,10)
    b = Box(P2(1,1), P2(5,5))
    v = view(g, b)
    @test v == CartesianGrid(P2(1,1), P2(5,5), dims=(4,4))

    p = PointSet(collect(vertices(g)))
    v = view(p, b)
    @test coordinates(v, 1) == T[1,1]
    @test coordinates(v, nelements(v)) == T[5,5]
  end

  @testset "Data" begin
    # dummy type implementing the Data trait
    struct DummyData{𝒟,𝒯} <: Data
      domain::𝒟
      table::𝒯
    end
    Meshes.domain(data::DummyData) = data.domain
    Meshes.values(data::DummyData) = data.table
    Meshes.constructor(::Type{D}) where {D<:DummyData} = DummyData

    g = CartesianGrid{T}(10,10)
    t = (a=1:100, b=1:100)
    d = DummyData(g, t)
    b = Box(P2(1,1), P2(5,5))
    v = view(d, b)
    @test domain(v) == CartesianGrid(P2(1,1), P2(5,5), dims=(4,4))
    @test values(v) == (a=[12,13,14,15,22,23,24,25,32,33,34,35,42,43,44,45], b=[12,13,14,15,22,23,24,25,32,33,34,35,42,43,44,45])

    p = PointSet(collect(vertices(g)))
    d = DummyData(p, t)
    v = view(d, b)
    dd = domain(v)
    @test coordinates(dd, 1) == T[1,1]
    @test coordinates(dd, nelements(dd)) == T[5,5]
    tt = values(v)
    @test tt == (a=[13,14,15,16,17,24,25,26,27,28,35,36,37,38,39,46,47,48,49,50,57,58,59,60,61],
                 b=[13,14,15,16,17,24,25,26,27,28,35,36,37,38,39,46,47,48,49,50,57,58,59,60,61])
  end

  @testset "Utilities" begin
    d = CartesianGrid{T}(10,10)
    s = slice(d, T(0.5):T(11), T(0.5):T(11))
    @test s == CartesianGrid(P2(1,1), P2(10,10), dims=(9,9))
  end
end


text:
function baz1()::String
    return "foo"
end


text:
#!/usr/bin/env julia
#
# Note: Nothing changed in this script. Only broken URLs removed.
#
function printsum(a)
    # generates a summary of an object
    println(summary(a), ": ", repr(a))
end

# arrays can be initialised directly:
a1 = [1,2,3]
printsum(a1)
#> 3-element Array{Int64,1}: [1,2,3]

# or initialised empty:
a2 = []
printsum(a2)
#> 0-element Array{None,1}: None[]

# since this array has no type, functions like push! (see below) don't work
# instead arrays can be initialised with a type:
a3 = Int64[]
printsum(a3)
#> 0-element Array{Int64,1}: []

# ranges are different from arrays:
a4 = 1:20
printsum(a4)
#> 20-element UnitRange{Int64}: 1:20

# however they can be used to create arrays thus:
a4 = collect(1:20)
printsum(a4)
#> 20-element Array{Int64,1}: [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]

# arrays can also be generated from comprehensions
a5 = [2^i for i = 1:10]
printsum(a5)
#> 10-element Array{Int64,1}: [2,4,8,16,32,64,128,256,512,1024]

# arrays can be any type, so arrays of arrays can be created:
a6 = (Array{Int64, 1})[]
printsum(a6)
#> 0-element Array{Array{Int64,1},1}: []
# (note this is a "jagged array" (i.e., an array of arrays),
# not a multidimensional array, these are not covered here)

# Julia provided a number of dequeues functions, the most common for appending
# to the end of arrays is push!
# ! at the end of a function name indicates that the first argument is updated.

push!(a1, 4)
printsum(a1)
#> 4-element Array{Int64,1}: [1,2,3,4]

# push!(a2, 1) would cause error:

push!(a3, 1)
printsum(a3) #> 1-element Array{Int64,1}: [1]
#> 1-element Array{Int64,1}: [1]

push!(a6, [1,2,3])
printsum(a6)
#> 1-element Array{Array{Int64,1},1}: [[1,2,3]]

# using repeat() to create arrays
# you must use the keywords "inner" and "outer"
# all arguments must be arrays (not ranges)
a7 = repeat(a1,inner=[2],outer=[1])
printsum(a7)
#> 8-element Array{Int64,1}: [1,1,2,2,3,3,4,4]
a8 = repeat(collect(4:-1:1),inner=[1],outer=[2])
printsum(a8)
#> 8-element Array{Int64,1}: [4,3,2,1,4,3,2,1]

max_stars_repo_name: wbeebe/julia

max_issues_repo_name: wbeebe/julia

max_forks_repo_name: wbeebe/julia


text:
using TypeDBORM
using Test

using TypeDBORM: _resulting_fieldtype

include(joinpath(dirname(@__FILE__),"types_structs.jl"))

@testset "util_functions" begin
    rel_type = _resulting_fieldtype(Employee, :address_rel)
    @info rel_type
    @test rel_type == Address
end

max_issues_repo_name: Humans-of-Julia/TypeDBORM

max_forks_repo_name: Humans-of-Julia/TypeDBORM


text:
const SBit = @NamedTuple{pos::Bool, neg::Bool}
Base.convert(::Type{SBit}, b::Tuple{Bool, Bool}) = (pos = b[1], neg = b[2])
Base.promote_rule(::Type{SBit}, ::Type{Tuple{Bool, Bool}}) = SBit

pos(b::SBit) = b.pos
neg(b::SBit) = b.neg

Base.zero(::Type{SBit}) = SBit((false, false))
Base.one(::Type{SBit}) = SBit((true, false))
SBit(value::Number) = isone(value) ? one(SBit) :
                       iszero(value) ? zero(SBit) :
                       error("Cannot create an SBit($value). Use pop!(SBitstream($value)) instead.")

Base.show(io::IO, b::SBit) = print(io, "($(b.pos), $(b.neg))")
Base.show(io::IO, ::MIME"text/plain", b::SBit) = print(io, "SBit(pos = $(b.pos), neg = $(b.neg))")

"""
    SBitstream(value::Real)
    SBitstream{T<:Real}(value)

A stochastic bitstream that represents a real (floating-point) number
between [-1, 1].

Fields:
- `bits::Vector{SBit}`: the underlying bitstream
- `value::Float64`: the underlying floating-point number being represented
"""
struct SBitstream{T<:Real} <: AbstractBitstream
    bits::Vector{SBit}
    value::T

    function SBitstream{T}(bits::Vector{SBit}, value::T) where {T<:Real}
        if value > 1 || value < -1
            @warn "SBitstream can only be ∈ [-1, 1] (saturation occurring)."
        end

        new{T}(bits, min(max(value, -1), 1))
    end
end
SBitstream(value::T) where {T<:Real} = SBitstream{T}(Vector{SBit}(), value)
SBitstream{T}(value::Real) where {T<:Real} = SBitstream(convert(T, value))
Base.convert(::Type{SBitstream{T}}, s::SBitstream) where {T<:Real} =
    SBitstream{T}(s.bits, convert(T, s.value))

const SBitstreamLike = Union{<:SBitstream, AbstractArray{<:SBitstream}}

Base.float(s::SBitstream) = s.value
Base.zero(::SBitstream{T}) where T = SBitstream(zero(T))
Base.one(::SBitstream{T}) where T = SBitstream(one(T))

Base.typemin(::Type{T}) where {T<:SBitstream} = zero(T)
Base.typemax(::Type{T}) where {T<:SBitstream} = one(T)

Base.isless(x::SBitstream, y::SBitstream) = Base.isless(float(x), float(y))
Base.:(==)(x::SBitstream, y::SBitstream) = (float(x) == float(y))
Base.hash(x::SBitstream, h::UInt) = hash(float(x), h)

Base.show(io::IO, s::SBitstream) = print(io, "SBitstream($(s.value))")
Base.show(io::IO, ::MIME"text/plain", s::SBitstream{T}) where T =
    print(io, "SBitstream{$T}(value = $(s.value))\n    with $(length(s)) bits.")

include("./soperators.jl")

Base.:(+)(x::SBitstream, y::SBitstream) = SBitstream(float(x) + float(y))

Base.:(-)(x::SBitstream, y::SBitstream) = SBitstream(float(x) - float(y))

Base.:(*)(x::SBitstream, y::SBitstream) = SBitstream(float(x) * float(y))

function Base.:(/)(x::SBitstream, y::SBitstream)
    if float(y) <= 0
        error("SBitstream only supports divisors > 0 (y == $y).")
    end

    SBitstream(float(x) / y.value)
end

for (op, sim) in ((:+, :SSignedAdder),
                  (:-, :SSignedSubtractor),
                  (:*, :SSignedMultiplier),
                  (:/, :SSignedDivider))
    @eval begin
        is_trace_primitive(::Type{typeof($op)}, ::Type{<:SBitstream}, ::Type{<:SBitstream}) = true
        is_trace_primitive(::Type{typeof(Base.broadcasted)},
                           ::Type{typeof($op)},
                           ::Type{<:SBitstreamLike},
                           ::Type{<:SBitstreamLike}) = true
        getsimulator(::typeof($op), x::SBitstream, y::SBitstream) = $(sim)()
        getsimulator(::typeof(Base.broadcasted),
                     ::typeof($op),
                     x::SBitstreamLike,
                     y::SBitstreamLike) = getsimulator.($op, x, y)
    end
end

function Base.:(÷)(x::SBitstream, y::Real)
    if y < 1
        error("SBitstream only supports fixed-gain divisors >= 1 (y == $y).")
    end

    SBitstream(float(x) / y)
end
is_trace_primitive(::Type{typeof(÷)}, ::Type{<:SBitstreamLike}, ::Type{<:Real}) = true
is_trace_primitive(::Type{typeof(Base.broadcasted)},
                   ::Type{typeof(÷)},
                   ::Type{<:SBitstreamLike},
                   ::Type{<:Real}) = true
getsimulator(::typeof(÷), x::SBitstream, y::Real) = SSignedFixedGainDivider()
getsimulator(::typeof(Base.broadcasted), ::typeof(÷), x::SBitstreamLike, y::Real) =
    getsimulator.(÷, x, y)

Base.sqrt(x::SBitstream) = SBitstream(sqrt(float(x)))
is_trace_primitive(::Type{typeof(sqrt)}, ::Type{<:SBitstreamLike}) = true
is_trace_primitive(::Type{typeof(Base.broadcasted)},
                   ::Type{typeof(sqrt)},
                   ::Type{<:SBitstreamLike}) = true
getsimulator(::typeof(sqrt), x::SBitstream) = SSquareRoot()
getsimulator(::typeof(Base.broadcasted), ::typeof(sqrt), x::SBitstream) = getsimulator.(sqrt, x)

decorrelate(x::SBitstream) = SBitstream(float(x))
is_trace_primitive(::Type{typeof(decorrelate)}, ::Type{<:SBitstream}) = true
is_trace_primitive(::Type{typeof(Base.broadcasted)},
                   ::Type{typeof(decorrelate)},
                   ::Type{<:SBitstreamLike}) = true
getsimulator(::typeof(decorrelate), x::SBitstream) = SSignedDecorrelator()
getsimulator(::typeof(Base.broadcasted), ::typeof(decorrelate), x::SBitstreamLike) =
    getsimulator.(decorrelate, x)

Base.:(*)(x::AbstractVecOrMat{<:SBitstream}, y::AbstractVecOrMat{<:SBitstream}) =
    SBitstream.(float.(x) * float.(y))
Base.:(*)(x::LinearAlgebra.Transpose{<:SBitstream, <:AbstractVector}, y::AbstractMatrix{<:SBitstream}) =
    SBitstream.(float.(x) * float.(y))
is_trace_primitive(::Type{typeof(*)},
                   ::Type{<:AbstractVecOrMat{<:SBitstream}},
                   ::Type{<:AbstractVecOrMat{<:SBitstream}}) = true
getsimulator(::typeof(*), x::AbstractVecOrMat{<:SBitstream}, y::AbstractVecOrMat{<:SBitstream}) =
    SSignedMatMultiplier(size(x, 1), size(y, 2))

LinearAlgebra.norm(x::AbstractVector{<:SBitstream}) = SBitstream(norm(float.(x)))
is_trace_primitive(::Type{typeof(LinearAlgebra.norm)}, ::Type{<:AbstractVector{<:SBitstream}}) = true
getsimulator(::typeof(LinearAlgebra.norm), x::AbstractVector{<:SBitstream}) = SL2Normer()

Base.max(x::SBitstream, y::SBitstream, zs::SBitstream...) = SBitstream(max(float(x), float(y), float.(zs)...))
is_trace_primitive(::Type{typeof(max)}, ::Type{<:SBitstream}, ::Type{<:SBitstream}, ::Type{<:SBitstream}...) = true
is_trace_primitive(::Type{typeof(Base.broadcasted)},
                   ::Type{typeof(max)},
                   ::Type{<:SBitstreamLike},
                   ::Type{<:SBitstreamLike},
                   ::Type{<:SBitstreamLike}...) = true
getsimulator(::typeof(max), x::SBitstream, y::SBitstream, zs::SBitstream...) = SSignedNMaxer(length(zs) + 2)
getsimulator(::typeof(max),
             x::SBitstreamLike,
             y::SBitstreamLike,
             zs::SBitstreamLike...) = getsimulator.(max, x, y, zs...)

bits(s::SBitstream) = s.bits

"""
    generate(s::SBitstream, T::Integer = 1)
    generate!(s::SBitstream, T::Integer = 1)

Generate `T` samples of the bitstream, `s`.
Add them to the queue for `generate!`, otherwise return the vector of bits.
"""
function generate(s::SBitstream, T::Integer = 1)
    bits = rand(T) .< abs(s.value)
    sbits = map(b -> (s.value >= 0) ? (pos = b, neg = false) : (pos = false, neg = b), bits)

    return sbits
end
generate!(s::SBitstream, T::Integer = 1) = push!(s, generate(s, T))

Base.push!(s::SBitstream, b) = push!(s.bits, b)
Base.pop!(s::SBitstream) = isempty(s.bits) ? generate(s)[1] : popfirst!(s.bits)
observe(s::SBitstream) = last(s.bits)

"""
    estimate(buffer::AbstractVector)
    estimate(s::SBitstream)

Get the empirical mean of the bits in `buffer`/`s`
"""
estimate(buffer::AbstractVector) = sum(buffer) / length(buffer)
estimate(s::SBitstream) = mapreduce(x -> x.pos - x.neg, +, s.bits) / length(s)

"""
    estimate!(buffer::AbstractVector, b::SBit)
    estimate!(buffer::AbstractVector, b::VecOrMat{SBit})
    estimate(buffer::AbstractVector)
    estimate(s::SBitstream)

Push `b` into the `buffer` and return the current [`estimate`](#).
"""
function estimate!(buffer::AbstractVector, s::SBitstream)
    b = observe(s)
    push!(buffer, pos(b) - neg(b))

    return estimate(buffer)
end
function estimate!(buffer::AbstractVector, s::VecOrMat{<:SBitstream})
    bs = observe(s)
    push!(buffer, pos.(bs) - neg.(bs))

    return estimate(buffer)
end


text:
using Immerse

y = rand(10^4)
p = plot(x=1:length(y),y=y,Geom.line)

figure()
hfig = display(p)


text:
import Pkg; Pkg.activate(".")
using CSV
using CmdStan
using ArviZ
using DelimitedFiles

# include functions
include("common_functions.jl")

# define cmdstan location
set_cmdstan_home!(homedir() * "/Desktop/cmdstan/")

# get the data and filter for predation only
d = CSV.read(joinpath(pwd(), "data", "network_data.dat"))
d = d[d.predation .> 0 , :]

### LSSL

const lssl = """
data{
    int W;
    int L[W];
    int S[W];
    int cf;
}
parameters{
    real a;
    real phi;
}
model{
    vector[W] mu;
    phi ~ normal( 2 , 1 );
    a ~ normal( 0.7 , 0.02 );
    for ( i in 1:W ) {
        mu[i] = exp(a) * S[i];
    }
    L ~ neg_binomial_2( mu , exp(phi) );
}
generated quantities{
    vector[W] log_lik;
    vector[W] mu;
    vector[W] y_hat;
    vector[cf] counterfactual_links;
    for ( i in 1:W ) {
        mu[i] = exp(a) * S[i];
        log_lik[i] = neg_binomial_2_lpmf( L[i] | mu[i] , exp(phi) );
        y_hat[i] = neg_binomial_2_rng(mu[i], exp(phi));
    }
    for (j in 1:cf){
     counterfactual_links[j] = neg_binomial_2_rng(exp(a) * j, exp(phi));
    }
}
"""

data_dict = Dict("W" => length(d.id),
    "L" => d.links,
    "S" => d.nodes,
    "cf" => 750)

lssl_stan_model = Stanmodel(
    model =lssl,
    nchains = 4,
    num_warmup = 1000,
    num_samples = 1000,
    name = "lssl"
)
_, lssl_stan_chns, _ = stan(lssl_stan_model, data_dict, summary = true);

### Constant connectance

const constant_connect = """
data{
    int W;
    int L[W];
    int S[W];
    int cf;
}
parameters{
    real<lower=0,upper=1> a;
    real phi;
}
model{
    vector[W] mu;
    phi ~ normal( 2 , 1 );
    a ~ beta( 3 , 7 );
    for ( i in 1:W ) {
        mu[i] = a * S[i]^2;
    }
    L ~ neg_binomial_2( mu , exp(phi) );
}
generated quantities{
    vector[W] log_lik;
    vector[W] mu;
    vector[W] y_hat;
    vector[cf] counterfactual_links;
    for ( i in 1:W ) {
        mu[i] = a * S[i] ^ 2;
        log_lik[i] = neg_binomial_2_lpmf( L[i] | mu[i] , exp(phi) );
        y_hat[i] = neg_binomial_2_rng(mu[i], exp(phi));
    }
    for (j in 1:cf){
     counterfactual_links[j] = neg_binomial_2_rng(a * j^2, exp(phi));
    }
}
"""

data_dict = Dict("W" => length(d.id),
    "L" => d.links,
    "S" => d.nodes,
    "cf" => 750)

const_connect_stan_model = Stanmodel(
    model =constant_connect,
    nchains = 4,
    num_warmup = 1000,
    num_samples = 1000,
    name = "constant_connect"
)
_, constant_connect_stan_chns, _ = stan(const_connect_stan_model, data_dict, summary = true);

##### power law connectance

const pwrlaw_connectance = """
data{
    int W;
    int L[W];
    int S[W];
    int cf;
}
parameters{
    real a;
    real b;
    real phi;
}
model{
    vector[W] mu;
    phi ~ normal( 2 , 1 );
    a ~ normal( -3 , 1 );
    b ~ normal(2, 0.6);
    for ( i in 1:W ) {
        mu[i] = exp(a) * S[i] ^ b;
    }
    L ~ neg_binomial_2( mu , exp(phi) );
}
generated quantities{
    vector[W] log_lik;
    vector[W] mu;
    vector[W] y_hat;
    vector[cf] counterfactual_links;
    for ( i in 1:W ) {
        mu[i] = exp(a) * S[i] ^ b;
        log_lik[i] = neg_binomial_2_lpmf( L[i] | mu[i] , exp(phi) );
        y_hat[i] = neg_binomial_2_rng(mu[i], exp(phi));
    }
    for (j in 1:cf){
     counterfactual_links[j] = neg_binomial_2_rng(exp(a) * j^b, exp(phi));
    }
}
"""

pwrlaw_conn_stan_model = Stanmodel(
    model = pwrlaw_connectance,
    nchains = 4,
    num_warmup = 1000,
    num_samples = 1000,
    name = "powerlaw_connectance"
)

_, pwrlaw_stan_chns, _ = stan(pwrlaw_conn_stan_model, data_dict, summary = false);

### beta-binomial connectance
const betabin_connectance = """
data{
    int W;
    int L[W];
    int S[W];
    int cf;
}
transformed data{
    int F[W];
    int R[W];
    int M[W];
    for ( i in 1:W ) {
        M[i] = S[i] - 1;
        F[i] = S[i] * S[i] - M[i];
        R[i] = L[i]        - M[i];
    }
}
parameters{
    real<lower=0,upper=1> mu;
    real phi;
}
model{
    phi ~ normal( 3,0.5 );
    mu ~ beta( 3 , 7 );
    for (i in 1:W){
       target += beta_binomial_lpmf(  R[i] | F[i] ,  mu * exp(phi) , (1 - mu) * exp(phi));
    }
}
generated quantities{
    vector[W] log_lik;
    vector[W] y_hat;
    vector[cf] counterfactual_links;
    for ( i in 1:W ) {
        log_lik[i] = beta_binomial_lpmf( R[i] | F[i] , mu * exp(phi), (1 - mu) * exp(phi)  );
        y_hat[i] = beta_binomial_rng(F[i] , mu * exp(phi), (1 - mu) * exp(phi) ) + M[i];
    }
    for (j in 1:cf){
        counterfactual_links[j] = beta_binomial_rng(j * j - j + 1 , mu * exp(phi), (1 - mu) * exp(phi) ) + j - 1;
    }
}
"""

data_dict = Dict(
    "W" => length(d.id),
    "L" => d.links,
    "S" => d.nodes,
    "cf" => 750)

bb_model = Stanmodel(
    model = betabin_connectance,
    nchains = 2,
    num_warmup = 1000,
    num_samples = 1000,
    name = "betabin_connectance"
)

_, bb_chains , _ = stan(bb_model, data_dict, summary = true);

##### write out posterior samples as csvs

size(lssl_stan_chns)

lssl_df = DataFrame(lssl_stan_chns)
CSV.write("data/posterior_distributions/lssl.csv", lssl_df, delim=',')

constant_connect_array = DataFrame(constant_connect_stan_chns)
CSV.write("data/posterior_distributions/const_posterior.csv", constant_connect_array, delim=',')

pwrlaw_df = DataFrame(pwrlaw_stan_chns)
CSV.write("data/posterior_distributions/powerlaw_posterior.csv", pwrlaw_df, delim=',')

bb_df = DataFrame(bb_chains)
CSV.write("data/posterior_distributions/beta_binomial_posterior.csv", bb_df, delim=',')

### more posterior samples for Flexible links model

data_dict_2 = Dict(
    "W" => length(d.id),
    "L" => d.links,
    "S" => d.nodes,
    "cf" => 1500)

bb_model = Stanmodel(
        model = betabin_connectance,
        nchains = 2,
        num_warmup = 1000,
        num_samples = 1000,
        name = "betabin_connectance_bigger"
    )

_, bb_chains_bigger , _ = stan(bb_model, data_dict_2, summary = true);

bb_df = DataFrame(bb_chains_bigger)
CSV.write("data/posterior_distributions/beta_binomial_posterior_bigger.csv", bb_df, delim=',')

## could be used to plot a "ribbon"
bb_hpd = MCMCChains.hpd(bb_chains, alpha = 0.89)

bb_hpd.df

###  calculation of LOO

bb_chains_infdata = foodweb_model_output(bb_chains)

lssl_stan_infdata = foodweb_model_output(lssl_stan_chns)

const_stan_infdata = foodweb_model_output(constant_connect_stan_chns)

pwrlaw_stan_infdata = foodweb_model_output(pwrlaw_stan_chns)

### get parameter estimates
summary(const_stan_infdata)
summary(pwrlaw_stan_infdata)
summary(bb_chains_infdata)
summary(lssl_stan_infdata)

loo(lssl_stan_infdata)
# 2940.68 +- 118.63
## calculate loo
loo(const_stan_infdata) #2798 +- 104
# 2940 +- 118 with a narrower prior

loo(pwrlaw_stan_infdata) # 2595 +- 49

loo(bb_chains_infdata) # 2543 +- 46

### this should produce pointwise loo calculations but I don't know how to get the numbers out
loo_pw = loo(pwrlaw_stan_infdata, pointwise = true)


text:
sets = ["Fam2a","Fam4a","Fam8a","Group2a","Group4a","Group5a","Wed2a","Wed3a","Wed5a"]

#http://chenlab.ece.cornell.edu/people/Andy/ImagesOfGroups.html

# 1		Female
# 2		Male
for i = 1:length(sets)

f = readlines("../original-data/Group_Pic_Data/$(sets[i])/PersonData.txt")
Men = Vector{Int64}()
People = Vector{Int64}()
r = 10
N = zeros(r,r+1)
for j = 1:length(f)
    global m, p
    if j == 1
        p = 0
        m = 0
        continue
    end
    line = split(f[j],"\t")
    if length(line) == 1
        @assert(occursin("@",string(line)))
        @assert(occursin(".jpg",string(line)))
        # Whenever you get to a new picture
        push!(Men,m)
        push!(People,p)
        if p < 11
            N[p,p-m+1] += 1
        end
        m = 0
        p = 0
    else
        p += 1
        m += parse(Int64,line[end])-1
    end
end
matwrite("N_$(sets[i]).mat",Dict("N"=>N))

end

## Combine them together

wed = ["Wed2a","Wed3a","Wed5a"]
Nw = zeros(10,11)
for i = 1:3
    global Nw
    M = matread("N_$(wed[i]).mat")
    Nw += M["N"]
end

group = ["Group2a","Group4a","Group5a"]
Ng = zeros(10,11)
for i = 1:3
    global Ng
    M = matread("N_$(group[i]).mat")
    Ng += M["N"]
end

fam = ["Fam2a","Fam4a"]
Nf = zeros(10,11)
for i = 1:2
    global Nf
    M = matread("N_$(fam[i]).mat")
    Nf += M["N"]
end

matwrite("Total_Ns.mat",Dict("Nf"=>Nf,"Nw"=>Nw,"Ng"=>Ng))

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
using TurbulenceConvection, Documenter
using DocumenterCitations
import Glob

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))

#! format: off
pages = Any[
    "Home" => "index.md",
    "Reference states" => "ReferenceStates.md",
    "References" => "References.md",
]

mathengine = MathJax(Dict(
    :TeX => Dict(
        :equationNumbers => Dict(:autoNumber => "AMS"),
        :Macros => Dict(),
    ),
))

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
)
#! format: on

makedocs(
    bib,
    sitename = "TurbulenceConvection.jl",
    strict = true,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [TurbulenceConvection],
    pages = pages,
)

# Clean up
build_dir = joinpath(@__DIR__, "build")
if isdir(build_dir)
    cd(build_dir) do
        for filename in Glob.glob("Output.*")
            rm(filename; force = true, recursive = true)
        end
        for filename in Glob.glob("*.in")
            rm(filename; force = true)
        end
    end
end

deploydocs(
    repo = "github.com/CliMA/TurbulenceConvection.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
a = Turing.getargnames(divorce_AM)
m5_3 = divorce_AM(eachcol(d[!, [a...]])...)
b = typeof(m5_3)
c = Turing.VarInfo(m5_3)
d = Turing.tonamedtuple(c) |> keys


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct TGFixed <: TurbineGov
        efficiency::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of a fixed Turbine Governor that returns a fixed mechanical torque
 given by the product of P_ref*efficiency

# Arguments
- `efficiency::Float64`:  Efficiency factor that multiplies P_ref, validation range: `(0, nothing)`
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int`: TGFixed has no states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct TGFixed <: TurbineGov
    " Efficiency factor that multiplies P_ref"
    efficiency::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    "TGFixed has no states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function TGFixed(efficiency, P_ref=1.0, ext=Dict{String, Any}(), )
    TGFixed(efficiency, P_ref, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function TGFixed(; efficiency, P_ref=1.0, ext=Dict{String, Any}(), states=Vector{Symbol}(), n_states=0, internal=InfrastructureSystemsInternal(), )
    TGFixed(efficiency, P_ref, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function TGFixed(::Nothing)
    TGFixed(;
        efficiency=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`TGFixed`](@ref) `efficiency`."""
get_efficiency(value::TGFixed) = value.efficiency
"""Get [`TGFixed`](@ref) `P_ref`."""
get_P_ref(value::TGFixed) = value.P_ref
"""Get [`TGFixed`](@ref) `ext`."""
get_ext(value::TGFixed) = value.ext
"""Get [`TGFixed`](@ref) `states`."""
get_states(value::TGFixed) = value.states
"""Get [`TGFixed`](@ref) `n_states`."""
get_n_states(value::TGFixed) = value.n_states
"""Get [`TGFixed`](@ref) `internal`."""
get_internal(value::TGFixed) = value.internal

"""Set [`TGFixed`](@ref) `efficiency`."""
set_efficiency!(value::TGFixed, val) = value.efficiency = val
"""Set [`TGFixed`](@ref) `P_ref`."""
set_P_ref!(value::TGFixed, val) = value.P_ref = val
"""Set [`TGFixed`](@ref) `ext`."""
set_ext!(value::TGFixed, val) = value.ext = val


text:
using StatsBase
#
# Auxiliary globals
#
const _empty        = 0
const _nonpolarized = 1

#
#  Auxiliary functions
#

function update(k,cell,lattice,_rates,prob,rng)
    # k is 1-based index

    iflag  = 0

    ipos   = cell[k]
    istate = lattice[ipos]
    # NOTE comment the double-check below
    if (istate < 1)
        return iflag,prob
    end
    #
    # time interval
    #
    dt = 1e0/prod(size(lattice)) # detailed scale
    #    dt = 1e0/size(cell)[1] # only for particle conservation
    #
    #  ONE  BODY transitios
    #
    itmp = Int( istate == _nonpolarized )
    # using total polarization rate! i.e., 4 times the indiv. polarz
    prob += _rates["polarization"]*dt*itmp
    if (rng < prob)
        # 4 is the coordination number. NOTE Use a variable
        lattice[ipos] = rand((_nonpolarized+1):(_nonpolarized+4))
        return 1,prob
    end
    itmp = Int(istate > _nonpolarized)
    prob += _rates["depolarization"]*dt*itmp
    if (rng < prob)
        lattice[ipos] = _nonpolarized
        return 2,prob
    end
    # get_neighbours works for nonpolarized as well
    L = [i for i in size(lattice)]
    neigh = get_neighbours(ipos,L)[istate]
    # cell is polarized and the target site is empty
    itmp = itmp * Int( lattice[neigh] == _empty  )
    prob += _rates["displacement"]*itmp*dt
    if (rng < prob)
        cell[k]        = neigh
        lattice[neigh] = istate
        lattice[ipos]  = _empty
        return 3,prob
    end
    return iflag,prob
end

function init(L,N0 = 1,inum_states=5)
    L_tot = prod(L)
    lattice = zeros(Int,L...)
    cell = StatsBase.sample(1:L_tot,N0,replace=false)
    lattice[cell] = rand(1:inum_states,N0)
    return cell,lattice
end

function get_neighbours(k,L)
    kshift  = k - 1 # because julia uses 1-based indices
    idims   = size(L)[1]
    icoords = get_coords(kshift,L)
    neighs  = zeros(Int,2*idims+1)
    neighs[1]=kshift
    # forward neighs
    for i=1:idims
        ioutro = icoords[:]
        ioutro[i] = mod(icoords[i]+1+L[i],L[i])
        neighs[i+1] = get_index(ioutro,L)
    end
    # back neighs
    for i=1:idims
        ioutro = icoords[:]
        ioutro[i] = mod(icoords[i]-1+L[i],L[i])
        neighs[i+idims+1] = get_index(ioutro,L)
    end
    return neighs + ones(Int,2*idims+1) # 1-based indices
end

function get_coords(k,L)
    # k is a 0-based index here
    idims      = size(L)[1]
    icoords    = zeros(Int,idims)
    icoords[1] = mod(k,L[1])
    ifrac      = L[1]
    for j = 2:idims
        icoords[j] = div(k,ifrac)
        ifrac = ifrac*L[j]
    end
    return icoords # 0-based grid index
end

function get_index(icoords,L)
    idims  = size(L)[1]
    index = icoords[1]
    ifrac = L[1]
    for i=2:idims
        index += icoords[i]*ifrac
        ifrac = ifrac*L[i]
    end
    return index # 0-based index
end