text:
@testset "Low-level" begin
    @testset "Definitions" begin
        test(StructDefinition{false}, struct_by_name, :VkPhysicalDeviceProperties, :(
            struct _PhysicalDeviceProperties <: VulkanStruct{false}
                vks::VkPhysicalDeviceProperties
            end))

        test(StructDefinition{false}, struct_by_name, :VkApplicationInfo, :(
            struct _ApplicationInfo <: VulkanStruct{true}
                vks::VkApplicationInfo
                deps::Vector{Any}
            end))

        test(StructDefinition{false}, struct_by_name, :VkExtent2D, :(
            struct _Extent2D <: VulkanStruct{false}
                vks::VkExtent2D
            end))

        test(StructDefinition{false}, struct_by_name, :VkExternalBufferProperties, :(
            struct _ExternalBufferProperties <: VulkanStruct{true}
                vks::VkExternalBufferProperties
                deps::Vector{Any}
            end
        ))

        test(StructDefinition{false}, struct_by_name, :VkPipelineExecutableInternalRepresentationKHR, :(
            struct _PipelineExecutableInternalRepresentationKHR <: VulkanStruct{true}
                vks::VkPipelineExecutableInternalRepresentationKHR
                deps::Vector{Any}
            end
        ))

        test(StructDefinition{false}, struct_by_name, :VkDescriptorSetAllocateInfo, :(
            struct _DescriptorSetAllocateInfo <: VulkanStruct{true}
                vks::VkDescriptorSetAllocateInfo
                deps::Vector{Any}
                descriptor_pool::DescriptorPool
            end
        ))

        test(StructDefinition{false}, struct_by_name, :VkAccelerationStructureBuildGeometryInfoKHR, :(
            struct _AccelerationStructureBuildGeometryInfoKHR <: VulkanStruct{true}
                vks::VkAccelerationStructureBuildGeometryInfoKHR
                deps::Vector{Any}
                src_acceleration_structure::OptionalPtr{AccelerationStructureKHR}
                dst_acceleration_structure::OptionalPtr{AccelerationStructureKHR}
            end
        ))
    end

    @testset "Friendly constructors" begin
        test(Constructor, StructDefinition{false}(struct_by_name(:VkPhysicalDeviceProperties)), :(
            function _PhysicalDeviceProperties(api_version::VersionNumber, driver_version::VersionNumber, vendor_id::Integer, device_id::Integer, device_type::PhysicalDeviceType, device_name::AbstractString, pipeline_cache_uuid::NTuple{Int(VK_UUID_SIZE), UInt8}, limits::_PhysicalDeviceLimits, sparse_properties::_PhysicalDeviceSparseProperties)
                _PhysicalDeviceProperties(VkPhysicalDeviceProperties(to_vk(UInt32, api_version), to_vk(UInt32, driver_version), vendor_id, device_id, device_type, device_name, pipeline_cache_uuid, limits.vks, sparse_properties.vks))
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkInstanceCreateInfo)), :(
            function _InstanceCreateInfo(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; next=C_NULL, flags=0, application_info=C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                application_info = cconvert(Ptr{VkApplicationInfo}, application_info)
                enabled_layer_names = cconvert(Ptr{Cstring}, enabled_layer_names)
                enabled_extension_names = cconvert(Ptr{Cstring}, enabled_extension_names)
                deps = Any[
                    next,
                    application_info,
                    enabled_layer_names,
                    enabled_extension_names
                ]
                vks = VkInstanceCreateInfo(structure_type(VkInstanceCreateInfo),
                                        unsafe_convert(Ptr{Cvoid}, next),
                                        flags,
                                        unsafe_convert(Ptr{VkApplicationInfo}, application_info),
                                        pointer_length(enabled_layer_names),
                                        unsafe_convert(Ptr{Cstring}, enabled_layer_names),
                                        pointer_length(enabled_extension_names),
                                        unsafe_convert(Ptr{Cstring}, enabled_extension_names),
                                        )
                _InstanceCreateInfo(vks, deps)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkSubpassSampleLocationsEXT)), :(
            function _SubpassSampleLocationsEXT(subpass_index::Integer, sample_locations_info::_SampleLocationsInfoEXT)
                _SubpassSampleLocationsEXT(VkSubpassSampleLocationsEXT(subpass_index, sample_locations_info.vks))
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDebugUtilsMessengerCreateInfoEXT)), :(
            function _DebugUtilsMessengerCreateInfoEXT(message_severity::DebugUtilsMessageSeverityFlagEXT, message_type::DebugUtilsMessageTypeFlagEXT, pfn_user_callback::FunctionPtr; next = C_NULL, flags = 0, user_data = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                user_data = cconvert(Ptr{Cvoid}, user_data)
                deps = Any[next, user_data]
                vks = VkDebugUtilsMessengerCreateInfoEXT(structure_type(VkDebugUtilsMessengerCreateInfoEXT), unsafe_convert(Ptr{Cvoid}, next), flags, message_severity, message_type, pfn_user_callback, unsafe_convert(Ptr{Cvoid}, user_data))
                _DebugUtilsMessengerCreateInfoEXT(vks, deps)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDescriptorSetAllocateInfo)), :(
            function _DescriptorSetAllocateInfo(descriptor_pool, set_layouts::AbstractArray; next = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                set_layouts = cconvert(Ptr{VkDescriptorSetLayout}, set_layouts)
                deps = Any[next, set_layouts]
                vks = VkDescriptorSetAllocateInfo(structure_type(VkDescriptorSetAllocateInfo), unsafe_convert(Ptr{Cvoid}, next), descriptor_pool, pointer_length(set_layouts), unsafe_convert(Ptr{VkDescriptorSetLayout}, set_layouts))
                _DescriptorSetAllocateInfo(vks, deps, descriptor_pool)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkPipelineShaderStageCreateInfo)), :(
            function _PipelineShaderStageCreateInfo(stage::ShaderStageFlag, _module, name::AbstractString; next = C_NULL, flags = 0, specialization_info = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                name = cconvert(Cstring, name)
                specialization_info = cconvert(Ptr{VkSpecializationInfo}, specialization_info)
                deps = Any[next, name, specialization_info]
                vks = VkPipelineShaderStageCreateInfo(structure_type(VkPipelineShaderStageCreateInfo), unsafe_convert(Ptr{Cvoid}, next), flags, VkShaderStageFlagBits(stage.val), _module, unsafe_convert(Cstring, name), unsafe_convert(Ptr{VkSpecializationInfo}, specialization_info))
                _PipelineShaderStageCreateInfo(vks, deps, _module)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDescriptorImageInfo)), :(
            function _DescriptorImageInfo(sampler, image_view, image_layout::ImageLayout)
                _DescriptorImageInfo(VkDescriptorImageInfo(sampler, image_view, image_layout), sampler, image_view)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDescriptorSetLayoutBinding)), :(
            function _DescriptorSetLayoutBinding(binding::Integer, descriptor_type::DescriptorType, stage_flags::ShaderStageFlag; descriptor_count = 0, immutable_samplers = C_NULL)
                immutable_samplers = cconvert(Ptr{VkSampler}, immutable_samplers)
                deps = Any[immutable_samplers]
                vks = VkDescriptorSetLayoutBinding(binding, descriptor_type, descriptor_count, stage_flags, unsafe_convert(Ptr{VkSampler}, immutable_samplers))
                _DescriptorSetLayoutBinding(vks, deps)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkXcbSurfaceCreateInfoKHR)), :(
            function _XcbSurfaceCreateInfoKHR(connection::Ptr{vk.xcb_connection_t}, window::vk.xcb_window_t; next = C_NULL, flags = 0)
                next = cconvert(Ptr{Cvoid}, next)
                connection = cconvert(Ptr{vk.xcb_connection_t}, connection)
                deps = Any[next, connection]
                vks = VkXcbSurfaceCreateInfoKHR(structure_type(VkXcbSurfaceCreateInfoKHR), unsafe_convert(Ptr{Cvoid}, next), flags, unsafe_convert(Ptr{vk.xcb_connection_t}, connection), window)
                _XcbSurfaceCreateInfoKHR(vks, deps)
            end
        ))
    end

    @testset "Manual tweaks" begin
        test(Constructor, StructDefinition{false}(struct_by_name(:VkWriteDescriptorSet)), :(
            function _WriteDescriptorSet(dst_set, dst_binding::Integer, dst_array_element::Integer, descriptor_type::DescriptorType, image_info::AbstractArray, buffer_info::AbstractArray, texel_buffer_view::AbstractArray; next = C_NULL, descriptor_count = max(pointer_length(image_info), pointer_length(buffer_info), pointer_length(texel_buffer_view)))
                next = cconvert(Ptr{Cvoid}, next)
                image_info = cconvert(Ptr{VkDescriptorImageInfo}, image_info)
                buffer_info = cconvert(Ptr{VkDescriptorBufferInfo}, buffer_info)
                texel_buffer_view = cconvert(Ptr{VkBufferView}, texel_buffer_view)
                deps = Any[next, image_info, buffer_info, texel_buffer_view]
                vks = VkWriteDescriptorSet(structure_type(VkWriteDescriptorSet), unsafe_convert(Ptr{Cvoid}, next), dst_set, dst_binding, dst_array_element, descriptor_count, descriptor_type, unsafe_convert(Ptr{VkDescriptorImageInfo}, image_info), unsafe_convert(Ptr{VkDescriptorBufferInfo}, buffer_info), unsafe_convert(Ptr{VkBufferView}, texel_buffer_view))
                _WriteDescriptorSet(vks, deps, dst_set)
            end
        ))
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
		Pkg.activate(".")
end
atreplinit() do repl
    try
        @eval using OhMyREPL
    catch e
        @warn "error while importing OhMyREPL" e
    end
    try
        @eval using Revise
    catch e
        @warn "error while importing Revise" e
    end
end

ext: jl

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# tests generate function

@generator GNStartFuncGen() begin # prefix GN (for GeNerate) to avoid type name clashes
    start() = 'a'
    other() = 'b'
end

@generator GNChoiceModelGen() begin # prefix GN (for GeNerate) to avoid type name clashes
    start() = plus(a)
    a() = 'a'
end

@generator GNIntGen begin
    start() = choose(Int,5,9)
end

@generator GNMainGen(intGen) begin
    start() = plus(choose(intGen))
end

@generator GNInfRecursionGen() begin
    start() = [plus(start)]
end

@generator GN100ChoicesGen() begin
    start() = map(i->a(),1:50) # each invocation of a requires two choice points
    a() = choose(Int,1,10)
    a() = plus(b)
    b() = 'b'
end

@generator GNMultGen() begin
    start() = mult(a)
    a() = :a
end

@testset "generate" begin

	@testset "generate using different start rules" begin

		gn = GNStartFuncGen()

		@testset "default start rule" begin
		    td = first(generate(gn))
		    @test td == 'a'
		end

		@testset "non-default start rule" begin
		    td = first(generate(gn, startrule=:other))
		    @test td == 'b'
		end

	end

	@testset "default choice model is a sampler choice model" begin

		gn = GNChoiceModelGen()
		@test isa(choicemodel(gn), DataGenerators.SamplerChoiceModel)

	end

	@testset "generate using default and non-default choice models" begin

		gn = GNChoiceModelGen()

		@test isa(choicemodel(gn), DataGenerators.SamplerChoiceModel)

		@mtestset "default choice model" reps=Main.REPS alpha=Main.ALPHA begin
	    	td = choose(gn)
			@mtest_values_include [1,2,3] length(td)
		end

		setminimumvaluechoicemodel!(gn)

		@mtestset "non-default choice model" reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
			@mtest_values_are [1,] length(td)
		end

	end

	@testset "generate a generator with sub-generators using default and non-default choice models" begin

		ign = GNIntGen()
		gn = GNMainGen(ign)

		@mtestset "default choice model" reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @mtest_values_include [5,6,7] first(td)
		end

		setminimumvaluechoicemodel!(gn)
		setminimumvaluechoicemodel!(ign)

		@mtestset "non-default choice model" reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @mtest_values_are [5,] first(td)
		end

	end

	@testset "limit on number of choice points (maxchoices)" begin

		infgn = GNInfRecursionGen()

		@testset "infinite recursion throws GenerationTerminatedException" begin
		    exc = nothing
		    try
		        td = choose(infgn)
		    catch e
		        if isa(e,GenerationTerminatedException)
		            exc = e
		        else
		            throw(e)
		        end
		    end
		    # TODO use @check_throws instead
		    @test typeof(exc) == GenerationTerminatedException
		    # check for sensible reason that specifies the limit
		    @test match(r"choices", exc.reason) != nothing
		    @test match(r"exceeded", exc.reason) != nothing
		    @test match(Regex("$(DataGenerators.defaultgenerateparams[:maxchoices])"), exc.reason) != nothing
		end

		c100gn = GN100ChoicesGen()

		@testset "number of choices <= maxchoices parameter to generator" begin
		    exc = nothing
		    try
		        td = choose(c100gn, maxchoices = 99)
		    catch e
		        if isa(e,GenerationTerminatedException)
		            exc = e
		        else
		            throw(e)
		        end
		    end
		    # TODO use @check_throws instead
		    @test typeof(exc) == GenerationTerminatedException
		end

		@testset "number of choices > maxchoices parameter to generator" begin
		    exc = nothing
		    try
		        td = choose(c100gn, maxchoices = 100)
		    catch e
		        if isa(e,GenerationTerminatedException)
		            exc = e
		        else
		            throw(e)
		        end
		    end
		    # TODO use @check_throws instead
		    @test typeof(exc) != GenerationTerminatedException
		end

		@testset "robustgen with number of choices <= maxchoices parameter to generator" begin
		    @test robustchoose(c100gn, maxchoices = 99) == nothing
		end

		@testset "robustgen with number of choices > maxchoices parameter to generator" begin
		    @test robustchoose(c100gn, maxchoices = 100) != nothing
		end

	end

	@testset "limit on sequence choice point reps (maxseqreps)" begin

		gn = GNMultGen()

		# note that to force sufficiently long sequences to reach the default limit
		# we need to use the sampler choice model rather than a default choice model
		# (but sampler choice model itself is only tested later)
		setchoicemodel!(gn, SamplerChoiceModel(gn))

		# sampler choice model should have a single parameter that controls the geometric distribution
		# set this to a value close to zero so that median length is much larger than the default maxseqreps
		setparams!(choicemodel(gn), [0.000001])

		@mtestset "sequence are no longer than the default limit" reps=Main.REPS alpha=Main.ALPHA begin
		    exc = nothing
		    td = nothing
		    try
		        td = choose(gn)
		    catch e
		        if isa(e,GenerationTerminatedException)
		            exc = e
		        else
		            throw(e)
		        end
		    end
		    @test (exc != nothing) || (length(td) <= DataGenerators.defaultgenerateparams[:maxseqreps]) # check that limit applies
		    @mtest_that_sometimes typeof(exc) == GenerationTerminatedException
		    if (exc != nothing)
		        @test match(r"repetitions", exc.reason) != nothing
		        @test match(r"exceeded", exc.reason) != nothing
		        @test match(Regex("$(DataGenerators.defaultgenerateparams[:maxseqreps])"), exc.reason) != nothing
		    end
		end

		@mtestset "sequence are no longer than the specified limit" reps=Main.REPS alpha=Main.ALPHA begin
		    exc = nothing
		    td = nothing
		    try
		        td = choose(gn, maxseqreps = 87)
		    catch e
		        if isa(e,GenerationTerminatedException)
		            exc = e
		        else
		            throw(e)
		        end
		    end
		    @test (exc != nothing) || (length(td) <= DataGenerators.defaultgenerateparams[:maxseqreps]) # check that limit applies
		    @mtest_that_sometimes typeof(exc) == GenerationTerminatedException
		    if (exc != nothing)
		        @test match(r"repetitions", exc.reason) != nothing
		        @test match(r"exceeded", exc.reason) != nothing
		        @test match(Regex("87"), exc.reason) != nothing
		    end
		end

	end

end


text:
# ClimateMachine driver configurations
#
# Contains helper functions to establish simulation configurations to be
# used with the ClimateMachine driver. Currently:
# - AtmosLESConfiguration
# - AtmosGCMConfiguration
# - OceanBoxGCMConfiguration
# - SingleStackConfiguration
# - MultiColumnLandModel
#
# User-customized configurations can use these as templates.

using CLIMAParameters
using CLIMAParameters.Planet: planet_radius

abstract type ConfigSpecificInfo end
struct AtmosLESSpecificInfo <: ConfigSpecificInfo end
struct AtmosGCMSpecificInfo{FT} <: ConfigSpecificInfo
    domain_height::FT
    nelem_vert::Int
    nelem_horz::Int
    vert_range::Array{FT, 1}
end
struct OceanBoxGCMSpecificInfo <: ConfigSpecificInfo end
struct SingleStackSpecificInfo{FT} <: ConfigSpecificInfo
    zmax::FT
    hmax::FT
end
struct MultiColumnLandSpecificInfo <: ConfigSpecificInfo end

include("SolverTypes/SolverTypes.jl")

"""
    ArgParse.parse_item

Parses custom command line option for tuples of three or fewer integers.
"""
function ArgParse.parse_item(::Type{NTuple{2, Int}}, s::AbstractString)

    str_array = split(s, ",")
    int_1 = length(str_array) > 0 ? parse(Int, str_array[1]) : 0
    int_2 = length(str_array) > 1 ? parse(Int, str_array[2]) : 0

    return (int_1, int_2)
end

"""
    ArgParse.parse_item

Parses custom command line option for tuples of three or fewer integers.
"""
function ArgParse.parse_item(::Type{NTuple{3, Int}}, s::AbstractString)

    str_array = split(s, ",")
    int_1 = length(str_array) > 0 ? parse(Int, str_array[1]) : 0
    int_2 = length(str_array) > 1 ? parse(Int, str_array[2]) : 0
    int_3 = length(str_array) > 2 ? parse(Int, str_array[3]) : 0

    return (int_1, int_2, int_3)
end

"""
    ArgParse.parse_item

Parses custom command line option for tuples of three integers.
"""
function ArgParse.parse_item(::Type{NTuple{3, Float64}}, s::AbstractString)

    str_array = split(s, ",")
    ft_1 = length(str_array) > 0 ? parse(Float64, str_array[1]) : 0
    ft_2 = length(str_array) > 1 ? parse(Float64, str_array[2]) : 0
    ft_3 = length(str_array) > 2 ? parse(Float64, str_array[3]) : 0

    return (ft_1, ft_2, ft_3)
end

"""
    ArgParse.parse_item

Parses custom command line option for `DateTime`s.
"""
function ArgParse.parse_item(::Type{DateTime}, s::AbstractString)
    DateTime(s, "yyyy-mm-ddTHH:MM:SS")
end

"""
    get_polyorders

Utility function that gets the polynomial orders for the given configuration
either passed from command line or as default values
"""
function get_polyorders(N, setting = ClimateMachine.Settings.degree)

    # Check if polynomial degree was passed as a CL option
    if setting != (-1, -1)
        setting
    elseif N isa Int
        (N, N)
    else
        N
    end
end

"""
    ClimateMachine.DriverConfiguration

Collects all parameters necessary to set up a ClimateMachine simulation.
"""
struct DriverConfiguration{FT}
    config_type::ClimateMachineConfigType

    name::String
    # Polynomial order tuple (polyorder_horz, polyorder_vert)
    polyorders::NTuple{2, Int}
    array_type::Any
    #
    # Model details
    param_set::AbstractParameterSet
    bl::BalanceLaw
    #
    # Execution details
    mpicomm::MPI.Comm
    #
    # Mesh details
    grid::DiscontinuousSpectralElementGrid
    #
    # DGModel details
    numerical_flux_first_order::NumericalFluxFirstOrder
    numerical_flux_second_order::NumericalFluxSecondOrder
    numerical_flux_gradient::NumericalFluxGradient
    # DGFVModel details, used when polyorder_vert = 0
    fv_reconstruction::Union{Nothing, AbstractReconstruction}
    # Cutoff filter to emulate overintegration
    filter::Any
    #
    # Configuration-specific info
    config_info::ConfigSpecificInfo

    function DriverConfiguration(
        config_type,
        name::String,
        polyorders::NTuple{2, Int},
        FT,
        array_type,
        param_set::AbstractParameterSet,
        bl::BalanceLaw,
        mpicomm::MPI.Comm,
        grid::DiscontinuousSpectralElementGrid,
        numerical_flux_first_order::NumericalFluxFirstOrder,
        numerical_flux_second_order::NumericalFluxSecondOrder,
        numerical_flux_gradient::NumericalFluxGradient,
        fv_reconstruction::Union{Nothing, AbstractReconstruction},
        filter,
        config_info::ConfigSpecificInfo,
    )
        return new{FT}(
            config_type,
            name,
            polyorders,
            array_type,
            param_set,
            bl,
            mpicomm,
            grid,
            numerical_flux_first_order,
            numerical_flux_second_order,
            numerical_flux_gradient,
            fv_reconstruction,
            filter,
            config_info,
        )
    end
end

function print_model_info(model, mpicomm)
    mpirank = MPI.Comm_rank(mpicomm)
    if mpirank == 0
        @show ClimateMachine.array_type()
        msg = "Model composition\n"
        for key in fieldnames(typeof(model))
            msg =
                msg * @sprintf(
                    "    %s = %s\n",
                    string(key),
                    string((getproperty(model, key)))
                )
        end
        @info msg
        show_tendencies(model; table_complete = model isa AtmosModel)
    end
end

function AtmosLESConfiguration(
    name::String,
    N::Union{Int, NTuple{2, Int}},
    (Δx, Δy, Δz)::NTuple{3, FT},
    xmax::FT,
    ymax::FT,
    zmax::FT,
    param_set::AbstractParameterSet,
    init_LES!;
    xmin = zero(FT),
    ymin = zero(FT),
    zmin = zero(FT),
    array_type = ClimateMachine.array_type(),
    physics = AtmosPhysics{FT}(param_set;),
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        init_state_prognostic = init_LES!,
    ),
    mpicomm = MPI.COMM_WORLD,
    boundary = ((0, 0), (0, 0), (1, 2)),
    periodicity = (true, true, false),
    meshwarp = (x...) -> identity(x),
    grid_stretching = (nothing, nothing, nothing),
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
    fv_reconstruction = nothing,
    Ncutoff = get_polyorders(N),
) where {FT <: AbstractFloat}

    (polyorder_horz, polyorder_vert) = get_polyorders(N)
    (cutofforder_horz, cutofforder_vert) =
        get_polyorders(Ncutoff, ClimateMachine.Settings.cutoff_degree)

    # Check if the element resolution was passed as a CL option
    if ClimateMachine.Settings.resolution != (-1, -1, -1)
        (Δx, Δy, Δz) = ClimateMachine.Settings.resolution
    end

    # Check if the min of the domain was passed as a CL option
    if ClimateMachine.Settings.domain_min != (-1, -1, -1)
        (xmin, ymin, zmin) = ClimateMachine.Settings.domain_min
    end

    # Check if the max of the domain was passed as a CL option
    if ClimateMachine.Settings.domain_max != (-1, -1, -1)
        (xmax, ymax, zmax) = ClimateMachine.Settings.domain_max
    end

    print_model_info(model, mpicomm)

    brickrange = (
        grid1d(
            xmin,
            xmax,
            grid_stretching[1],
            elemsize = Δx * max(polyorder_horz, 1),
        ),
        grid1d(
            ymin,
            ymax,
            grid_stretching[2],
            elemsize = Δy * max(polyorder_horz, 1),
        ),
        grid1d(
            zmin,
            zmax,
            grid_stretching[3],
            elemsize = Δz * max(polyorder_vert, 1),
        ),
    )
    topology = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = (polyorder_horz, polyorder_vert),
        meshwarp = meshwarp,
    )

    if cutofforder_horz < polyorder_horz && cutofforder_vert < polyorder_vert
        filter = CutoffFilter(
            grid,
            (cutofforder_horz + 1, cutofforder_horz + 1, cutofforder_vert + 1),
        )
    else
        filter = nothing
        cutofforder_horz = cutofforder_vert = nothing
    end

    @info @sprintf(
        """
Establishing Atmos LES configuration for %s
    precision               = %s
    horiz polynomial order  = %d
    vert polynomial order   = %d
    horiz cutoff order      = %s
    vert cutoff order       = %s
    domain_min              = %.2f m, %.2f m, %.2f m
    domain_max              = %.2f m, %.2f m, %.2f m
    resolution              = %dx%dx%d
    MPI ranks               = %d
    min(Δ_horz)             = %.2f m
    min(Δ_vert)             = %.2f m""",
        name,
        FT,
        polyorder_horz,
        polyorder_vert,
        cutofforder_horz,
        cutofforder_vert,
        xmin,
        ymin,
        zmin,
        xmax,
        ymax,
        zmax,
        Δx,
        Δy,
        Δz,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        AtmosLESConfigType(),
        name,
        (polyorder_horz, polyorder_vert),
        FT,
        array_type,
        param_set,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        fv_reconstruction,
        filter,
        AtmosLESSpecificInfo(),
    )
end

function AtmosGCMConfiguration(
    name::String,
    N::Union{Int, NTuple{2, Int}},
    (nelem_horz, nelem_vert)::NTuple{2, Int},
    domain_height::FT,
    param_set::AbstractParameterSet,
    init_GCM!;
    array_type = ClimateMachine.array_type(),
    physics = AtmosPhysics{FT}(param_set),
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        physics;
        init_state_prognostic = init_GCM!,
    ),
    mpicomm = MPI.COMM_WORLD,
    grid_stretching = (nothing,),
    meshwarp::Function = equiangular_cubed_sphere_warp,
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
    fv_reconstruction = nothing,
    Ncutoff = get_polyorders(N),
) where {FT <: AbstractFloat}

    (polyorder_horz, polyorder_vert) = get_polyorders(N)
    (cutofforder_horz, cutofforder_vert) =
        get_polyorders(Ncutoff, ClimateMachine.Settings.cutoff_degree)

    # Check if the number of elements was passed as a CL option
    if ClimateMachine.Settings.nelems != (-1, -1, -1)
        (nelem_horz, nelem_vert) = ClimateMachine.Settings.nelems
    end

    # Check if the domain height was passed as a CL option
    if ClimateMachine.Settings.domain_height != -1
        domain_height = ClimateMachine.Settings.domain_height
    end

    print_model_info(model, mpicomm)

    _planet_radius::FT = planet_radius(param_set)
    vert_range = grid1d(
        _planet_radius,
        FT(_planet_radius + domain_height),
        grid_stretching[1],
        nelem = nelem_vert,
    )

    topology = StackedCubedSphereTopology(
        mpicomm,
        nelem_horz,
        vert_range;
        boundary = (1, 2),
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = (polyorder_horz, polyorder_vert),
        meshwarp = meshwarp,
    )

    if cutofforder_horz < polyorder_horz && cutofforder_vert < polyorder_vert
        filter = CutoffFilter(
            grid,
            (cutofforder_horz + 1, cutofforder_horz + 1, cutofforder_vert + 1),
        )
    else
        filter = nothing
        cutofforder_horz = cutofforder_vert = nothing
    end

    @info @sprintf(
        """
Establishing Atmos GCM configuration for %s
    precision               = %s
    horiz polynomial order  = %d
    vert polynomial order   = %d
    horiz cutoff order      = %s
    vert cutoff order       = %s
    # horiz elem            = %d
    # vert elems            = %d
    domain height           = %.2e m
    MPI ranks               = %d
    min(Δ_horz)             = %.2f m
    min(Δ_vert)             = %.2f m""",
        name,
        FT,
        polyorder_horz,
        polyorder_vert,
        cutofforder_horz,
        cutofforder_vert,
        nelem_horz,
        nelem_vert,
        domain_height,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        AtmosGCMConfigType(),
        name,
        (polyorder_horz, polyorder_vert),
        FT,
        array_type,
        param_set,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        fv_reconstruction,
        filter,
        AtmosGCMSpecificInfo(
            domain_height,
            nelem_vert,
            nelem_horz,
            collect(vert_range),
        ),
    )
end

function OceanBoxGCMConfiguration(
    name::String,
    N::Union{Int, NTuple{2, Int}},
    (Nˣ, Nʸ, Nᶻ)::NTuple{3, Int},
    param_set::AbstractParameterSet,
    model;
    FT = Float64,
    array_type = ClimateMachine.array_type(),
    mpicomm = MPI.COMM_WORLD,
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
    fv_reconstruction = nothing,
    periodicity = (false, false, false),
    boundary = ((1, 1), (1, 1), (2, 3)),
)

    (polyorder_horz, polyorder_vert) = get_polyorders(N)

    # Check if the number of elements was passed as a CL option
    if ClimateMachine.Settings.nelems != (-1, -1, -1)
        (Nˣ, Nʸ, Nᶻ) = ClimateMachine.Settings.nelems
    end

    brickrange = (
        range(FT(0); length = Nˣ + 1, stop = model.problem.Lˣ),
        range(FT(0); length = Nʸ + 1, stop = model.problem.Lʸ),
        range(FT(-model.problem.H); length = Nᶻ + 1, stop = 0),
    )

    topology = StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = (polyorder_horz, polyorder_vert),
    )

    @info @sprintf(
        """
Establishing Ocean Box GCM configuration for %s
    precision               = %s
    horiz polynomial order  = %d
    vert polynomial order   = %d
    Nˣ # elems              = %d
    Nʸ # elems              = %d
    Nᶻ # elems              = %d
    domain depth            = %.2e m
    MPI ranks               = %d""",
        name,
        FT,
        polyorder_horz,
        polyorder_vert,
        Nˣ,
        Nʸ,
        Nᶻ,
        -model.problem.H,
        MPI.Comm_size(mpicomm)
    )

    return DriverConfiguration(
        OceanBoxGCMConfigType(),
        name,
        (polyorder_horz, polyorder_vert),
        FT,
        array_type,
        param_set,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        fv_reconstruction,
        nothing, # filter
        OceanBoxGCMSpecificInfo(),
    )
end

function SingleStackConfiguration(
    name::String,
    N::Union{Int, NTuple{2, Int}},
    nelem_vert::Int,
    zmax::FT,
    param_set::AbstractParameterSet,
    model::BalanceLaw;
    zmin = zero(FT),
    hmax = one(FT),
    array_type = ClimateMachine.array_type(),
    mpicomm = MPI.COMM_WORLD,
    boundary = ((0, 0), (0, 0), (1, 2)),
    periodicity = (true, true, false),
    meshwarp = (x...) -> identity(x),
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
    fv_reconstruction = nothing,
) where {FT <: AbstractFloat}

    (polyorder_horz, polyorder_vert) = get_polyorders(N)

    # Check if the number of elements was passed as a CL option
    if ClimateMachine.Settings.nelems != (-1, -1, -1)
        nE = ClimateMachine.Settings.nelems
        nelem_vert = nE[1]
    end

    # Check if the domain height was passed as a CL option
    if ClimateMachine.Settings.domain_height != -1
        zmax = ClimateMachine.Settings.domain_height
    end

    print_model_info(model, mpicomm)

    xmin, xmax = zero(FT), hmax
    ymin, ymax = zero(FT), hmax
    brickrange = (
        grid1d(xmin, xmax, nelem = 1),
        grid1d(ymin, ymax, nelem = 1),
        grid1d(zmin, zmax, nelem = nelem_vert),
    )
    topology = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = (polyorder_horz, polyorder_vert),
        meshwarp = meshwarp,
    )

    @info @sprintf(
        """
Establishing single stack configuration for %s
    precision               = %s
    horiz polynomial order  = %d
    vert polynomial order   = %d
    domain_min              = %.2f m, %.2f m, %.2f m
    domain_max              = %.2f m, %.2f m, %.2f m
    # vert elems            = %d
    MPI ranks               = %d
    min(Δ_horz)             = %.2f m
    min(Δ_vert)             = %.2f m""",
        name,
        FT,
        polyorder_horz,
        polyorder_vert,
        xmin,
        ymin,
        zmin,
        xmax,
        ymax,
        zmax,
        nelem_vert,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        SingleStackConfigType(),
        name,
        (polyorder_horz, polyorder_vert),
        FT,
        array_type,
        param_set,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        fv_reconstruction,
        nothing, # filter
        SingleStackSpecificInfo(zmax, hmax),
    )
end

function MultiColumnLandModel(
    name::String,
    N::Union{Int, NTuple{2, Int}},
    (Δx, Δy, Δz)::NTuple{3, FT},
    xmax::FT,
    ymax::FT,
    zmax::FT,
    param_set::AbstractParameterSet,
    model::BalanceLaw;
    xmin = zero(FT),
    ymin = zero(FT),
    zmin = zero(FT),
    array_type = ClimateMachine.array_type(),
    mpicomm = MPI.COMM_WORLD,
    boundary = ((3, 4), (5, 6), (1, 2)),
    solver_type = ExplicitSolverType(),
    periodicity = (false, false, false),
    meshwarp = (x...) -> identity(x),
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
    fv_reconstruction = nothing,
) where {FT <: AbstractFloat}

    (polyorder_horz, polyorder_vert) = get_polyorders(N)

    # Check if the element resolution was passed as a CL option
    if ClimateMachine.Settings.resolution != (-1, -1, -1)
        (Δx, Δy, Δz) = ClimateMachine.Settings.resolution
    end

    # Check if the min of the domain was passed as a CL option
    if ClimateMachine.Settings.domain_min != (-1, -1, -1)
        (xmin, ymin, zmin) = ClimateMachine.Settings.domain_min
    end

    # Check if the max of the domain was passed as a CL option
    if ClimateMachine.Settings.domain_max != (-1, -1, -1)
        (xmax, ymax, zmax) = ClimateMachine.Settings.domain_max
    end

    print_model_info(model, mpicomm)

    brickrange = (
        grid1d(xmin, xmax, elemsize = Δx * max(polyorder_horz, 1)),
        grid1d(ymin, ymax, elemsize = Δy * max(polyorder_horz, 1)),
        grid1d(zmin, zmax, elemsize = Δz * max(polyorder_vert, 1)),
    )

    topology = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = (polyorder_horz, polyorder_vert),
        meshwarp = meshwarp,
    )

    @info @sprintf(
        """
Establishing MultiColumnLandModel configuration for %s
    precision        = %s
    vert polyn order = %d
    horz polyn order = %d
    domain_min       = %.2f m, %.2f m, %.2f m
    domain_max       = %.2f m, %.2f m, %.2f m
    resolution       = %dx%dx%d
    MPI ranks        = %d
    min(Δ_horz)      = %.2f m
    min(Δ_vert)      = %.2f m""",
        name,
        FT,
        polyorder_vert,
        polyorder_horz,
        xmin,
        ymin,
        zmin,
        xmax,
        ymax,
        zmax,
        Δx,
        Δy,
        Δz,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        MultiColumnLandConfigType(),
        name,
        (polyorder_horz, polyorder_vert),
        FT,
        array_type,
        param_set,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        fv_reconstruction,
        nothing, # filter
        MultiColumnLandSpecificInfo(),
    )
end

import ..DGMethods: DGModel

"""
    DGModel(driver_config; kwargs...)

Initialize a [`DGModel`](@ref) given a
[`DriverConfiguration`](@ref) and keyword
arguments supported by [`DGModel`](@ref).
"""
DGModel(driver_config; kwargs...) = DGModel(
    driver_config.bl,
    driver_config.grid,
    driver_config.numerical_flux_first_order,
    driver_config.numerical_flux_second_order,
    driver_config.numerical_flux_gradient;
    gradient_filter = driver_config.filter,
    tendency_filter = driver_config.filter,
    kwargs...,
)

"""
    SpaceDiscretization(driver_config; kwargs...)

Initialize a [`SpaceDiscretization`](@ref) given a
[`DriverConfiguration`](@ref) and keyword
arguments supported by [`SpaceDiscretization`](@ref).
"""
SpaceDiscretization(driver_config; kwargs...) =
    (driver_config.polyorders[2] == 0) ?
    DGFVModel(
        driver_config.bl,
        driver_config.grid,
        driver_config.fv_reconstruction,
        driver_config.numerical_flux_first_order,
        driver_config.numerical_flux_second_order,
        driver_config.numerical_flux_gradient;
        gradient_filter = driver_config.filter,
        tendency_filter = driver_config.filter,
        kwargs...,
    ) :
    DGModel(
        driver_config.bl,
        driver_config.grid,
        driver_config.numerical_flux_first_order,
        driver_config.numerical_flux_second_order,
        driver_config.numerical_flux_gradient;
        gradient_filter = driver_config.filter,
        tendency_filter = driver_config.filter,
        kwargs...,
    )


text:
#
# Selection Sort
#

function selectionSort(A)
#=
Choose max value of A[1:i]
Insert it to the correct position
=#
  i = length(A)
  while i > 0
    max_pos = selectMax(A, 1,  i)
    if max_pos != i
      tmp = A[i]
      A[i] = A[max_pos]
      A[max_pos] = tmp
    end
    i = i - 1
  end

end

function selectMax(A, left, right)
#=
Select a max value index of A[left:right]
=#

  max_pos = left
  i = left
  while i <= right
    if A[i] > A[max_pos]
      max_pos = i
    end
    i = i + 1
  end

   return max_pos

end

function main()

  # list to be sorted
  rng = MersenneTwister(1234);
  A = shuffle(rng, Vector(1:100))
  Aorigin = copy(A)

  # sorting
  @time selectionSort(A)

  print("The original list: ")
  println(Aorigin)

  print("The sorted list: ")
  println(A)

end

if contains(@__FILE__, PROGRAM_FILE)
    main()
end

max_stars_repo_name: szmlb/algorithmsInJulia

max_issues_repo_name: szmlb/algorithmsInJulia

max_forks_repo_name: szmlb/algorithmsInJulia


text:
@testset "Recorder" begin
    learner = test_learner()

    tlr = Recorder.TrainLoss()
    vlr = Recorder.ValidateLoss()
    tslr = Recorder.SmoothTrainLoss()
    vslr = Recorder.SmoothValidateLoss()

    add_cb!(learner,tlr)
    add_cb!(learner,vlr)
    add_cb!(learner,tslr)
    add_cb!(learner,vslr)

    run_learner(learner)

    @test tlr[1,1] > tlr[16,1]
    @test vlr[1,1] > vlr[16,1]
    @test tslr[1,1] > tslr[16,1]
    @test vslr[1,1] > vslr[16,1]
end


text:
typeof(2) #> Int64
typeof(Int64) #> DataType
typeof(Complex{Int64}) #> DataType
typeof(Any) #> DataType
typeof(DataType) #> DataType

supertype(Int64) #> Signed
supertype(Signed) #> Integer
supertype(Integer) #> Real
supertype(Real) #> Number
supertype(Number) #> Any
supertype(Any) #> Any

using InteractiveUtils
subtypes(Integer) #> 3-element Array{Any,1}:
# Bool
# Signed
# Unsigned
subtypes(Signed) #> 6-element Array{Any,1}:
# BigInt
# Int128
# Int16
# Int32
# Int64
# Int8
subtypes(Int64) #> 0-element Array{Any,1}
typeof(Nothing) #> DataType

Bool <: Integer #> true
Bool <: Any #> true
Bool <: Char #> false
Bool <: Integer #> true
Float64 <: Integer #> false

typeof(5) #> Int64
isa(5, Number) #> true


text: waitforkey(prompt) = (print(stdout, prompt); read(stdin, 1); nothing)


text:
# best practice: Ctrl-D out of REPL first, then run test
using Test
using LightGraphs
using LinearAlgebra
include("FFGmain.jl")
include("lattices.jl")

ffgmathtest = @testset "test-FFGMATH/LATTICES" begin
    ############### testing lattice routines #################################
    dims = (5,5);
    lg = Lattices.get_SquareLattice(dims,pbc=true).graph
    @test sort(unique(neighborhood(lg, 3,2))) ==
    sort([3,4,5,2,1,9,8,7,13,23,18,24,22])

    dims = (4,4)
    lattice = Lattices.get_MixedBoundarySquareLattice(dims)
    F = Lattices.get_Tightbinding_Wavefunctions(lattice)
    @test F.values ≈ [-3.414213562373095,-3.414213562373095,
    -1.4142135623730951,-1.4142135623730951,-1.4142135623730951,-1.4142135623730951,
    -0.5857864376269051,-0.5857864376269051,0.5857864376269051,0.5857864376269051,
    1.4142135623730951,1.4142135623730951,1.4142135623730951,1.4142135623730951,
    3.414213562373095,3.414213562373095]

    dims = (2,3);
    lattice = Lattices.get_SquareLattice(dims,pbc=true);
    @test Matrix(adjacency_matrix(lattice.graph)) == [0 1 1 0 1 0;
                                                     1 0 0 1 0 1;
                                                     1 0 0 1 1 0;
                                                     0 1 1 0 0 1;
                                                     1 0 1 0 0 1;
                                                     0 1 0 1 1 0]
    F = Lattices.get_Tightbinding_Wavefunctions(lattice);
    evals = F.values;
    @test evals ≈ [-3, -1, 0, 0, 2, 2]
    states = F.vectors[:,1:3]
    r = [1,2,3]
    ############### testing math routines #####################################
    @test FFG.get_Determinant(r,states) ≈ 0.3333333333
    tmat = [0 -1.2247448713915894 -1.2247448713915894;
            -1.2247448713915894 1.2247448713915894 0;
            1 0 -1]
    @test FFG.get_Inverse_Matrix(r,states) ≈ tmat

    ############### set up small matrix to test det ratio factor & inverse #####
    m1 = [1 2 4;
          4 8 12;
          8 12 16];
    m1inv = inv(m1);
    m2 = [1 4 4;
           4 12 12;
           8 15 16];

    u = [2 ;4 ;3 ];
    v = [0 ; 1 ; 0 ];

    @test FFG.det_ratio_factor(m1inv,u,v) ≈ det(m2)/det(m1)

    @test FFG.sherman_morrison_inverse_update(m1inv,u,v) ≈ inv(m2)
end

@testset "test-FFGMAIN/HELPERS" begin
    ############# test chain initialization #######################
    model = Dict(
        "dims" => (2,3),
        "lattice" => Lattices.get_SquareLattice((2,3),pbc=true),
        "filling" => 1,
        "hamiltonian" => Lattices.get_Tightbinding_Wavefunctions,
        "fermi_energy" => 0,
        )
    policy = FFG.SwapNeighborsPolicy()
    mc_spec = Dict("mc_steps" => 1,
                        "mc_warmup_steps" => 1)

    # just tests that the function doesn't error out
    c1 = FFG.GutzwillerChain()
    @test FFG.init_Chain!(c1,model=model,observable=Nothing,policy=policy,
    mc_spec=mc_spec) == Nothing

    c2 = FFG.DoubleGutzwillerChain()
    @test FFG.init_Chain!(c2,model=model,observable=Nothing,policy=policy,
    mc_spec=mc_spec) == Nothing

    ############ make a test chain $$$$$$$$$$$$$$$$$$$$$$$$

    model = Dict(
        "dims" => (2,3),
        "lattice" => Lattices.get_SquareLattice((2,3),pbc=true),
        "filling" => 1,
        "hamiltonian" => Lattices.get_Tightbinding_Wavefunctions,
        "fermi_energy" => 0,
        )

    function observable(chain::FFG.GutzwillerChain)
        state = FFG.get_State(chain)
        Sz1 = state.spin_config.sc[1]
        return Dict("Sz1" => Sz1)
    end

    policy = FFG.SwapNeighborsPolicy()
    mc_spec = Dict("mc_steps" => 1,
                    "mc_warmup_steps" => 1)

    testchain = FFG.GutzwillerChain();
    testchain.basechain = FFG.MarkovChain()

    r_up = [3,1,5]
    r_down = [2,6,4]

    # STATE CREATION TESTS
    bd = [FFG.Electron(2,"up"), FFG.Electron(1,"down"), FFG.Electron(1,"up"),
        FFG.Electron(3,"down"), FFG.Electron(3,"up"), FFG.Electron(2,"down")]
    @test FFG.make_business_directory(r_up,r_down) == bd

    spinc = FFG.SpinConfiguration([1,-1,1,-1,1,-1])
    @test FFG.make_spin_config(r_up,r_down).sc == spinc.sc

    # testing construction of initial bond graph
    using LightGraphs: SimpleGraph, add_edge!
    ib = SimpleGraph(6);
    add_edge!(ib,1,2);
    add_edge!(ib,3,4);
    add_edge!(ib,5,6);

    @test FFG.get_init_bonds(model["lattice"].graph,spinc) ==
    ib

    ## STATE UPDATE TESTS
    teststate = FFG.GutzwillerState(
    r_up,r_down, bd, spinc,ib,
    [], 0, 0, [],[]
    )

    testchain.basechain.state = teststate
    testchain.basechain.model = model

    ## test a simple move
    move = FFG.ExchangeMove((1,2))

    @test FFG.get_proposal_factor_ratio(testchain,move) == 3/7
    FFG.update_Rs!(teststate,move)
    @test teststate.r_up == [3,2,5]
    @test teststate.r_down == [1,6,4]

    FFG.update_Bonds!(testchain,move)

    ub = SimpleGraph(6);
    add_edge!(ub,1,2);
    add_edge!(ub,1,3);
    add_edge!(ub,1,5);
    add_edge!(ub,2,4);
    add_edge!(ub,2,6);
    add_edge!(ub,3,4);
    add_edge!(ub,5,6);
    @test teststate.bonds == ub

    FFG.update_Spin_config!(teststate,move)
    @test teststate.spin_config.sc == [-1,1,1,-1,1,-1]

    bd = [FFG.Electron(1,"down"),FFG.Electron(2,"up"),FFG.Electron(1,"up"),
        FFG.Electron(3,"down"),FFG.Electron(3,"up"),FFG.Electron(2,"down")]
    FFG.update_Business_directory!(teststate,move)
    @test teststate.business_directory == bd

end

@testset "maintest" begin
    include("maintest.jl")
    @test size(maintest()) == (3,)
end

@testset "swaproutines" begin
    """ Test suite for swap routines. """
    dims = (6,) # dimension of the lattice
    lattice = Lattices.get_SquareLattice(dims,pbc=true) # make a square lattice
    filling = 1 # setting filling ≢ 1 means there are holes.
    model = Dict(
        "dims" => dims,
        "lattice" => lattice,
        "filling" => filling,
        "hamiltonian" => Lattices.get_Tightbinding_Wavefunctions,
        "fermi_energy" => 0,
        )
    function observable(chain::FFG.GutzwillerChain)
        state = FFG.get_State(chain)
        Sz1 = state.spin_config.sc[1]
        data = Dict("Sz1" => Sz1)
        return data
    end
    policy = FFG.SwapNeighborsPolicy()
    mc_spec = Dict("mc_steps" => 1E3,
                    "mc_warmup_steps" => 1E3,
                    "sample_interval" => 1E2)

    dg = FFG.DoubleGutzwillerChain()
    c1 = dg.replicas[1]
    c2 = dg.replicas[2]
    c1.basechain = FFG.MarkovChain()
    c2.basechain = FFG.MarkovChain()
    c1.basechain.model = model
    c2.basechain.model = model

    r_up1 = [5,6,4]
    r_down1 = [2,1,3]
    r_up2 = [2,5,1]
    r_down2 = [6,4,3]

    teststate1 = FFG.get_test_state(c1,r_up1,r_down1)
    teststate2 = FFG.get_test_state(c2,r_up2,r_down2)

    dg.replicas[1].basechain.state = teststate1
    dg.replicas[2].basechain.state = teststate2

    # first lets test a region that has a single exchange
    swapsites = [1,4,5]

    (spinc1,spinc2) = (teststate1.spin_config.sc,teststate2.spin_config.sc)

    Tinder1 = FFG.ElectronTinder([],[])
    Tinder2 = FFG.ElectronTinder([],[])
    for site = swapsites

        if spinc1[site] != spinc2[site]
            FFG.addUser!(Tinder1,site,spinc1[site])
            FFG.addUser!(Tinder2,site,spinc2[site])
        end
    end
    @test FFG.SwipeRight(Tinder1) == FFG.SwipeRight(FFG.ElectronTinder([4],[1]))
    @test FFG.SwipeRight(Tinder2) == FFG.SwipeRight(FFG.ElectronTinder([1],[4]))

    # now comes the difficult test...

    #
    r_up1prime = [5,6,1]
    r_down1prime = [2,4,3]
    r_up2prime = [2,5,4]
    r_down2prime = [6,1,3]
    #
    # ts1p = FFG.get_test_state(c1,r_up1prime,r_down1prime)
    # ts2p = FFG.get_test_state(c2,r_up2prime,r_down2prime)

    # print(ts1p.det_A_up * ts1p.det_A_down * ts2p.det_A_up * ts2p.det_A_down,"\n")
    # print(teststate1.det_A_up * teststate1.det_A_down * teststate2.det_A_up * teststate2.det_A_down,"\n")

    # manual_eval = (ts1p.det_A_up * ts1p.det_A_down * ts2p.det_A_up * ts2p.det_A_down)/(teststate1.det_A_up * teststate1.det_A_down * teststate2.det_A_up * teststate2.det_A_down)
    @test FFG.compute_swapregion(dg,swapsites) ≈ 1

end


text:
#! /usr/bin/julia

# Rosetta Code, Bitmap/Read a PPM file

using Color, Images, FixedPointNumbers

const M_RGB_Y = reshape(Color.M_RGB_XYZ[2,:], 3)

function rgb2gray(img::Image)
    g = red(img)*M_RGB_Y[1] + green(img)*M_RGB_Y[2] + blue(img)*M_RGB_Y[3]
    g = clamp(g, 0.0, 1.0)
    return grayim(g)
end

ima = imread("bitmap_read_ppm_in.ppm")
imb = convert(Image{Gray{Ufixed8}}, ima)
imwrite(imb, "bitmap_read_ppm_out.png")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# The four adjacent digits in the 1000-digit number that have the greatest
# product are 9 × 9 × 8 × 9 = 5832.
#
# 73167176531330624919225119674426574742355349194934
# 96983520312774506326239578318016984801869478851843
# 85861560789112949495459501737958331952853208805511
# 12540698747158523863050715693290963295227443043557
# 66896648950445244523161731856403098711121722383113
# 62229893423380308135336276614282806444486645238749
# 30358907296290491560440772390713810515859307960866
# 70172427121883998797908792274921901699720888093776
# 65727333001053367881220235421809751254540594752243
# 52584907711670556013604839586446706324415722155397
# 53697817977846174064955149290862569321978468622482
# 83972241375657056057490261407972968652414535100474
# 82166370484403199890008895243450658541227588666881
# 16427171479924442928230863465674813919123162824586
# 17866458359124566529476545682848912883142607690042
# 24219022671055626321111109370544217506941658960408
# 07198403850962455444362981230987879927244284909188
# 84580156166097919133875499200524063689912560717606
# 05886116467109405077541002256983155200055935729725
# 71636269561882670428252483600823257530420752963450
#
# Find the thirteen adjacent digits in the 1000-digit number that have the
# greatest product. What is the value of this product?

using ProjectEulerSolutions

# Use Julia built-ins to convert string into individual numbers then multiply
# every adjacent chunck.
function p008solution(number::String="23", adjacent::Integer=1)::Integer
    array = map(x->parse(Int8, x), split(number, ""))
    # Below, another way to split, using digits function, not quite as fast
    # array = digits(Int8, parse(BigInt, input))
    products = map(x -> prod(array[x:x+adjacent-1]), 1:(length(array)-adjacent+1))
    return maximum(products)
end

p008 = Problems.Problem(p008solution)

input = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"

Problems.benchmark(p008, input, 13)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
"""
    rand(N::ProbabilisticNetwork)

Converts a probabilistic network into a deterministic one by performing random
draws. All interactions are treated as independent Bernoulli events. Note that
this network is *not* checked for degeneracy, *i.e.* species can end up with no
interactions.

#### References

- Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
  D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
  Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function Base.rand(N::ProbabilisticNetwork)
    # Get the correct network type
    newtype = typeof(N) <: AbstractUnipartiteNetwork ? UnipartiteNetwork : BipartiteNetwork
    int = interactions(N)
    I = zeros(Int64, length(int))
    J = zeros(Int64, length(int))
    V = zeros(Bool, length(int))
    for (i,idata) in enumerate(int)
        I[i] = findfirst(isequal(idata.from), species(N; dims=1))
        J[i] = findfirst(isequal(idata.to), species(N; dims=2))
        V[i] = rand() <= idata.probability
    end
    edges = sparse(I, J, V, richness(N; dims=1), richness(N; dims=2))
    dropzeros!(edges)
    K = newtype(edges, EcologicalNetworks._species_objects(N)...)
    return K
end

"""
    rand(N::ProbabilisticNetwork, n::T) where {T<:Integer}

Generates a number of random deterministic networks based on a probabilistic
network.

#### References

- Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
  D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
  Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function Base.rand(N::ProbabilisticNetwork, n::T) where {T<:Integer}
    @assert n > 0
    newtype = typeof(N) <: AbstractUnipartiteNetwork ? UnipartiteNetwork : BipartiteNetwork
    R = Vector{newtype}(undef, n)
    for i in 1:n
        R[i] = rand(N)
    end
    return R
end

"""
    rand(N::ProbabilisticNetwork, S::Tuple{Int64,Int64})

Generates a number of random deterministic networks based on a probabilistic
network, and returns them as a matrix.

#### References

- Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
  D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
  Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function Base.rand(N::ProbabilisticNetwork, S::Tuple{Int64,Int64})
    @assert minimum(S) > 0
    return reshape(rand(N, prod(S)), S)
end


text:
# Format is defined by
# https://www.ibm.com/support/knowledgecenter/en/SSQ2R2_14.2.0/com.ibm.rsar.analysis.codereview.cobol.doc/topics/cac_useresults_junit.html
# http://help.catchsoftware.com/display/ET/JUnit+Format

"""
    set_attribute!(node, attr, val)

Add the attritube with name `attr` and value `val` to `node`.
"""
set_attribute!(node, attr, val) = setindex!(node, string(val), attr)

"""
    testsets_xml(name, id, ntests, nfails, nerrors, x_children)

Create the testsets element of a JUnit XML.
"""
function testsets_xml(name, id, ntests, nfails, nerrors, x_children)
    x_testset = ElementNode("testsets")
    link!.(Ref(x_testset), x_children)
    set_attribute!(x_testset, "name", name)
    set_attribute!(x_testset, "id", id)
    set_attribute!(x_testset, "tests", ntests)
    set_attribute!(x_testset, "failures", nfails)
    set_attribute!(x_testset, "errors", nerrors)
    x_testset
end

"""
    testset_xml(name, id, ntests, nfails, nerrors, x_children)

Create a testset element of a JUnit XML.
"""
function testset_xml(name, id, ntests, nfails, nerrors, x_children)
    x_testset = ElementNode("testset")
    link!.(Ref(x_testset), x_children)
    set_attribute!(x_testset, "name", name)
    set_attribute!(x_testset, "id", id)
    set_attribute!(x_testset, "tests", ntests)
    set_attribute!(x_testset, "failures", nfails)
    set_attribute!(x_testset, "errors", nerrors)
    x_testset
end

"""
    testcase_xml(name, id, x_children)

Create a testcase element of a JUnit XML.

This is the generic form (with name, id and children) that is used by other methods.
"""
function testcase_xml(name, id, x_children)
    x_tc = ElementNode("testcase")
    link!.(Ref(x_tc), x_children)
    set_attribute!(x_tc, "name", name)
    set_attribute!(x_tc, "id", id)
    x_tc
end

"""
    testcase_xml(v::Result, childs)

Create a testcase element of a JUnit XML for the result `v`.

The original expression of the test is used as the name, whilst the id is defaulted to
_testcase_id_.
"""
testcase_xml(v::Result, childs) = testcase_xml(string(v.orig_expr), "_testcase_id_", childs)

"""
    failure_xml(message, test_type, content)

Create a failure node (which will be the child of a testcase).
"""
function failure_xml(message, test_type, content)
    x_fail = ElementNode("failure")
    set_attribute!(x_fail, "message", message)
    set_attribute!(x_fail, "type", test_type)
    link!(x_fail, TextNode(content))
    x_fail
end

"""
    skip_xml()

Create a skip node (which will be the child of a testcase).
"""
function skip_xml()
    ElementNode("skip")
end

"""
    failure_xml(message, test_type, content)

Create an error node (which will be the child of a testcase).
"""
function error_xml(message, ex_type, content)
    x_fail = ElementNode("error")
    set_attribute!(x_fail, "message", message)
    set_attribute!(x_fail, "type",  ex_type)
    link!(x_fail, TextNode(content))
    x_fail
end

#####################

"""
    report(ts::AbstractTestSet)

Will produce an `XMLDocument` that contains a report about the `TestSet`'s results.
This report will follow the JUnit XML schema.

In theory this works on many kinds of `TestSet`s, but it is primarily intended for use on
`ReportingTestSet`s. To report correctly, the `TestSet` must have the following structure:

    AbstractTestSet
      └─ AbstractTestSet
           └─ Result

That is, the results of the top level `TestSet` must all be `AbstractTestSet`s,
and the results of those `TestSet`s must all be `Results`.
"""
function TestReports.report(ts_rich::RichReportingTestSet)
    ts = ts_rich.flattened_reporting_test_set[]
    total_ntests = 0
    total_nfails = 0
    total_nerrors = 0
    x_testsets = map(ts.results) do result
        x_testset, ntests, nfails, nerrors = to_xml(result)
        total_ntests += ntests
        total_nfails += nfails
        total_nerrors += nerrors;
        x_testset
    end

    xdoc = XMLDocument()
    root = setroot!(xdoc, testsets_xml(
        get_description(ts_rich),
        "_id_",
        total_ntests,
        total_nfails,
        total_nerrors,
        x_testsets))

    xdoc
end

"""
    to_xml(ts::AbstractTestSet)

Create a testset node from an `AbstractTestSet`, by creating nodes for each result
in `ts.results`. For creating a JUnit XML, all results must be `Result`s, that is
they cannot be `AbstractTestSet`s, as the XML cannot have one testset nested inside
another.
"""
function to_xml(ts_rich::RichReportingTestSet)
    ts = ts_rich.reporting_test_set[]
    total_ntests = 0
    total_nfails = 0
    total_nerrors = 0

    x_testcases = map(ts.results) do result
        x_testcase, ntests, nfails, nerrors = to_xml(result)
        total_ntests += ntests
        total_nfails += nfails
        total_nerrors += nerrors
        set_attribute!(x_testcase, "classname", get_description(ts))
        result isa Pass && set_attribute!(x_testcase, "name", x_testcase["name"] * " (Test $total_ntests)")
        x_testcase
    end

    x_testset = testset_xml(get_description(ts_rich), "_id_", total_ntests, total_nfails, total_nerrors, x_testcases)
    ts isa ReportingTestSet && add_testset_properties!(x_testset, ts)

    include_test_outputs(ts_rich, x_testset)

    x_testset, total_ntests, total_nfails, total_nerrors
end

function include_test_outputs(ts_rich, x_testset)
    output_children = []
    seekstart(ts_rich.out_buff)
    out_data = read(ts_rich.out_buff, String)
    if !isempty(out_data)
        x_output = ElementNode("std_out_data")
        link!(x_output, TextNode(out_data))
        push!(output_children, x_output)
    end
    err_data = read(ts_rich.err_buff, String)
    if !isempty(err_data)
        x_err = ElementNode("std_err_data")
        link!(x_output, TextNode(err_data))
        push!(output_children, x_err)
    end
    if !isempty(output_children)
        link!.(Ref(x_testset), output_children)
    end
end

function to_xml(ts::ReportingTestSet)
    return TestReports.to_xml(ts)
end

"""
    to_xml(res::Pass)
    to_xml(res::Fail)
    to_xml(res::Broken)
    to_xml(res::Error)

Create a testcase node from the result and return it along with
information on number of tests.
"""
function to_xml(res::Pass)
    x_testcase = testcase_xml("pass (info lost)", "_testcase_id_", [])
    x_testcase, 1, 0, 0  # Increment number of tests by 1
end

function to_xml(v::Fail)
    data = v.data === nothing ? "" : v.data  # Needed for V1.0
    x_failure = failure_xml(string(data), string(v.test_type), string(v))
    x_testcase = testcase_xml(v, [x_failure])
    x_testcase, 1, 1, 0  # Increment number of tests and number of failures by 1
end

function to_xml(v::Broken)
    x_testcase = testcase_xml(v, [skip_xml()]) # it is not actually skipped
    x_testcase, 0, 0, 0
end

function to_xml(v::Error)
    x_testcase = error_xml(string(v.value), typeof(v.value), v.backtrace)
    x_testcase, 0, 0, 1  # Increment number of errors by 1
end

"""
    add_testset_properties!(x_testset, ts::ReportingTestSet)

Add all key value pairs in the `properties` field of a `ReportingTestSet` to the
corresponding testset xml element.
"""
function add_testset_properties!(x_testset, ts::ReportingTestSet)
    if !isempty(keys(ts.properties))
        x_properties = ElementNode("properties")
        for (name, value) in ts.properties
            x_property = ElementNode("property")
            set_attribute!(x_property, "name", name)
            set_attribute!(x_property, "value", value)
            link!(x_properties, x_property)
        end
        link!(x_testset, x_properties)
    end
    return x_testset
end


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@info("Basic test of PIPotential construction and evaluation")
maxdeg = 10
r0 = 1.0
rcut = 3.0
trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
D = ACE1.SparsePSHDegree()
P1 = ACE1.BasicPSH1pBasis(Pr; species = :X, D = D)
basis = ACE1.PIBasis(P1, 2, D, maxdeg)
c = ACE1.Random.randcoeffs(basis)
V = combine(basis, c)
Nat = 15
Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, :X)
val_basis = real(sum(c .* evaluate(basis, Rs, Zs, z0)))
val_V = evaluate(V, Rs, Zs, z0)
println(@test(val_basis ≈ val_V))
J = evaluate_d(basis, Rs, Zs, z0)
grad_basis = real(sum(c[i] * J[i,:] for i = 1:length(c)))[:]
grad_V = evaluate_d(V, Rs, Zs, z0)
println(@test(grad_basis ≈ grad_V))

#---
# using BenchmarkTools
# Vgr = ACE1.GraphPIPot(V)
# tmp = ACE1.alloc_temp(V, Nat)
# tmpgr = ACE1.alloc_temp(Vgr, Nat)
# @btime ACE1.evaluate!($tmp, $V, $Rs, $Zs, $z0)
# @btime ACE1.evaluate!($tmpgr, $Vgr, $Rs, $Zs, $z0)
# @btime ACE1.evaluate_old!($tmpgr, $Vgr, $Rs, $Zs, $z0)
#
# #---
#
# @code_warntype ACE1.evaluate!(tmpgr, Vgr, Rs, Zs, z0)


text:
@testset "TestCompat" begin

#---

using JuLIP, ACE
import ZipFile, JSON

include("../artifacts.jl")
testdatadir = joinpath(datadir, "tests")

#---

@info("Test auto-generated testset `v07_compat`")
zipname = joinpath(testdatadir, "v07_compat.zip")
zipdir = ZipFile.Reader(zipname)

for fptr in zipdir.files
   print("    $(fptr.name): ")
   D = JSON.parse(fptr)
   calc = read_dict(D)
   ACE.Testing.runtests(calc, D["_tests"])
   if fptr.name[1:8] == "rpibasis"
      basis = ACE.Testing.test_basis(D["_params"])
      print(" / ")
      ACE.Testing.runtests(basis, D["_tests"])
   end
   println()
end

close(zipdir)

#---

@info("Test auto-generated testset `v08_compat`")
zipname = joinpath(testdatadir, "v08_compat.zip")
zipdir = ZipFile.Reader(zipname)

for fptr in zipdir.files
   print("    $(fptr.name): ")
   D = JSON.parse(fptr)
   calc = read_dict(D)
   ACE.Testing.runtests(calc, D["_tests"])
   if fptr.name[1:8] == "rpibasis"
      basis = ACE.Testing.test_basis(D["_params"])
      print(" / ")
      ACE.Testing.runtests(basis, D["_tests"])
   end
   println()
end

close(zipdir)

end


text:
module Sobol

import NLopt.libnlopt
import Base: ndims, length, skip, next, start, done, show, writemime
export SobolSeq, ScaledSobolSeq, next!

type SobolSeq{N}
    p::Ptr{Void} # nlopt_sobol pointer
    function SobolSeq()
        (N < 1 || N > 1111) && error("invalid Sobol dimension")
        p = ccall((:nlopt_sobol_create,libnlopt), Ptr{Void}, (Cuint,), N)
        p == C_NULL && error("error initializing Sobol sequence")
        s = new(p)
        finalizer(s, destroy)
        s
    end
end
SobolSeq(N::Integer) = SobolSeq{N}()

destroy(s::SobolSeq) =
    ccall((:nlopt_sobol_destroy,libnlopt), Void, (Ptr{Void},), s.p)

ndims{N}(s::SobolSeq{N}) = N

function show{N}(io::IO, s::SobolSeq{N})
    print(io, "$N-dimensional Sobol sequence on [0,1]^$N")
end
function writemime{N}(io::IO, ::MIME"text/html", s::SobolSeq{N})
    print(io, "$N-dimensional Sobol sequence on [0,1]$N")
end

function skip!{N}(s::SobolSeq{N}, n::Integer, x::Vector{Float64})
    (n < 0 || length(x) < N) && throw(BoundsError())
    ccall((:nlopt_sobol_skip,libnlopt), Void, (Ptr{Void}, Cuint, Ptr{Float64}),
          s.p, n, x)
    x
end

function skip{N}(s::SobolSeq{N}, n::Integer)
    skip!(s, n, Array(Float64, N))
    nothing
end

function next!{N}(s::SobolSeq{N}, x::Vector{Float64})
    length(x) < N && throw(BoundsError())
    ccall((:nlopt_sobol_next01,libnlopt), Void,
          (Ptr{Void}, Ptr{Float64}), s.p, x)
    x
end
next{N}(s::SobolSeq{N}) = next!(s, Array(Float64,N))

function next!{N}(s::SobolSeq{N}, x::Vector{Float64},
                  lb::Vector, ub::Vector)
    length(x) < N && throw(BoundsError())
    ccall((:nlopt_sobol_next,libnlopt), Void,
          (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
          s.p, x, float64(lb), float64(ub))
    x
end
next{N}(s::SobolSeq{N}, lb::Vector, ub::Vector) =
    next!(s, Array(Float64,N), lb, ub)

# Make an iterator so that we can do "for x in SobolSeq(...)".
# Technically, the Sobol sequence ends after 2^32-1 points, but NLopt
# falls back on pseudorandom numbers after this.  In practice, one is
# unlikely to reach that point.
start(s::SobolSeq) = s
next(s::SobolSeq, s_::SobolSeq) = (next(s), s_)
done(s::SobolSeq, s_::SobolSeq) = false

# Convenience wrapper for scaled Sobol sequences

type ScaledSobolSeq{N}
    s::SobolSeq{N}
    lb::Vector{Float64}
    ub::Vector{Float64}
    ScaledSobolSeq(lb::Vector{Float64}, ub::Vector{Float64}) =
        new(SobolSeq(N), lb, ub)
end
SobolSeq(N::Integer, lb, ub) =
    ScaledSobolSeq{N}(copy!(Array(Float64,N),lb), copy!(Array(Float64,N),ub))

next!(s::ScaledSobolSeq, x::Vector{Float64}) = next!(s.s, x, s.lb, s.ub)
next{N}(s::ScaledSobolSeq{N}) = next!(s, Array(Float64,N))

start(s::ScaledSobolSeq) = s
next(s::ScaledSobolSeq, s_::ScaledSobolSeq) = (next(s), s_)
done(s::ScaledSobolSeq, s_::ScaledSobolSeq) = false

skip(s::ScaledSobolSeq, n) = skip(s.s, n)

function show{N}(io::IO, s::ScaledSobolSeq{N})
    lb = s.lb; ub = s.ub
    print(io, "$N-dimensional scaled Sobol sequence on [$(lb[1]),$(ub[1])]")
    cnt = 1
    for i = 2:N
        if lb[i] == lb[i-1] && ub[i] == ub[i-1]
            cnt += 1
        else
            if cnt > 1
                print(io, "^", cnt)
            end
            print(io, " x [", lb[i], ",", ub[i], "]")
            cnt = 1
        end
    end
    if cnt > 1
        print(io, "^", cnt)
    end
end

function writemime{N}(io::IO, ::MIME"text/html", s::ScaledSobolSeq{N})
    lb = s.lb; ub = s.ub
    print(io, "$N-dimensional scaled Sobol sequence on [$(lb[1]),$(ub[1])]")
    cnt = 1
    for i = 2:N
        if lb[i] == lb[i-1] && ub[i] == ub[i-1]
            cnt += 1
        else
            if cnt > 1
                print(io, "", cnt, "")
            end
            print(io, " × [", lb[i], ",", ub[i], "]")
            cnt = 1
        end
    end
    if cnt > 1
        print(io, "", cnt, "")
    end
end

end # module


text:
using Documenter
using MPIMeasurements

makedocs(
    sitename = "MPIMeasurements",
    authors = "Tobias Knopp et al.",
    format = Documenter.HTML(prettyurls = false),
    modules = [MPIMeasurements],
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "Devices" => "man/devices.md",
            "Protocols" => "man/protocols.md",
            "Sequences" => "man/sequences.md",
            "Examples" => "man/examples.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
            "Internals" => map(
                s -> "lib/internals/$(s)",
                sort(readdir(normpath(@__DIR__, "src/lib/internals")))
            ),
        ],
    ],
)

# deploydocs(repo   = "github.com/MagneticParticleImaging/MPIMeasurements.jl.git",
#           target = "build")


text:
using DataFrames
using CSV
using Plots

function main()

   theme(:ggplot2)
   data = CSV.read("../data/clusters.csv", DataFrame, header=false)

   plot(data[!, 1], data[!, 2], data[!, 3], color=data[!, 4], st=:scatter, label="")
   title!("Data set")

   savefig("../images/clusters.png")
end

main()

max_stars_repo_name: axmat/Parallel-Kmeans

max_issues_repo_name: axmat/ParallelKMeans

max_forks_repo_name: axmat/ParallelKMeans


text:
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    DFMStruct(factorlags::Int64, errorlags::Int64, ndraws::Int64, burnin::Int64)

Description:
1-level DFM lag structure specification and MCMC sample size for Bayesian estimation.

Inputs:
- factorlags = Number of lags in the autoregressive specification of the latent factors.
- errorlags = Number of lags in the autoregressive specification of the observable variable idiosyncratic errors.
- ndraws = Number of MCMC draws used for posterior distributions.
- burnin = Number of initial MCMC draws discarded.
"""
@with_kw mutable struct DFMStruct
    factorlags::Int64
    errorlags::Int64
    ndraws::Int64
    burnin::Int64
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    HDFMStruct(nlevels::Int64, nfactors::Array{Int64,1}, factorassign::Array{Int64,2}, factorlags::Array{Int64,1}, errorlags::Array{Int64,1}, ndraws::Int64, burnin::Int64)

Description:
Multi-level/hierarchical DFM (HDFM) level, factor assignment, and lag structure specification, and MCMC sample size for Bayesian estimation.

Inputs:
- nlevels = Number of levels in the HDFM specification.
- nvars = Number of observable variables in the HDFM specification.
- nfactors = Number of factor per level in the HDFM specification.
- factorassign = Factors assigned to each variable across all levels.
- factorlags = Number of lags in the autoregressive specification of the latent factors.
- errorlags = Number of lags in the autoregressive specification of the observable variable idiosyncratic errors.
- ndraws = Number of MCMC draws used for posterior distributions.
- burnin = Number of initial MCMC draws discarded.
"""
@with_kw mutable struct HDFMStruct
    nlevels::Int64                  # number of levels in the multi-level model structure
    nfactors::Array{Int64,1}        # number of factors for each level (vector of length `nlevels`)
    factorassign::Array{Int64,2}         # integer matrix of size `nvar` × `nlevels`
    factorlags::Array{Int64,1}           # number of autoregressive lags for each factor level (vector of length `nlevels`)
    errorlags::Array{Int64,1}         # number of obs. variable error autoregressive lags (vector of length `nvar`)
    ndraws::Int64
    burnin::Int64
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    DFMMeans(F::Array{Float64}, B::Array{Float64}, S::Array{Float64}, P::Array{Float64}, P2::Array{Float64})

Description:
HDFM Bayesian estimator-generated latent factor and hyperparameter sample means (expected values).

Inputs:
- F = MCMC-generated latent factor sample mean.
- B = MCMC-generated observation equation regression coefficient sample means.
- S = MCMC-generated observable variable idiosyncratic error disturbance variance sample means.
- P = MCMC-generated latent factor autoregressive coefficient sample means.
- P2 = MCMC-generated idiosyncratic error autoregressive coefficient sample means.
"""
@with_kw mutable struct DFMMeans
    F::Array{Float64}   # Factor means
    B::Array{Float64}   # Obs. equation coefficient means
    S::Array{Float64}   # Idiosyncratic disturbance variance means
    P::Array{Float64}   # Factor autoregressive coefficient means
    P2::Array{Float64}  # Idiosyncratic disturbance autoregressive means
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    DFMResults(F::Array{Float64}, B::Array{Float64}, S::Array{Float64}, P::Array{Float64}, P2::Array{Float64}, means::DFMMeans)

Description:
HDMF Bayesian estimator-generated MCMC posterior distribution samples and their means for latent factors and hyperparameters.

Inputs:
- F = MCMC-generated latent factor sample.
- B = MCMC-generated observation equation regression coefficient sample.
- S = MCMC-generated observable variable idiosyncratic error disturbance variance sample.
- P = MCMC-generated latent factor autoregressive coefficient sample.
- P2 = MCMC-generated idiosyncratic error autoregressive coefficient sample.
- means = HDFM Bayesian estimator-generated latent factor and hyperparameter sample means (expected values).
"""
@with_kw mutable struct DFMResults
    F::Array{Float64}   # Factor sample
    B::Array{Float64}   # Obs. equation coefficient sample
    S::Array{Float64}   # Idiosyncratic disturbance variance sample
    P::Array{Float64}   # Factor autoregressive coefficient sample
    P2::Array{Float64}  # Idiosyncratic disturbance autoregressive sample
    means::DFMMeans     # Factor and hyperparameter means
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################


text:
export getNodalMassMatrix

function getNodalMassMatrix(mesh::TetraMesh, sigma::Array{Float64,1})

	nc = mesh.nc
  nn = mesh.nn

	M = zeros(4,4,nc)

	for i=1:nc
		Ti = mesh.Tetras[i,:]
	  	Xi = mesh.Points[Ti,:]

		M[:,:,i] = getLocalNodalMassMatrix(sigma[i],Xi)

	end
	nodalInd = mesh.Tetras

	ii = zeros(Int64,16*nc)
	jj = zeros(Int64,16*nc)
	# mm = zeros(Float64,16*nc)

	cnt = 1
	for i=1:nc
		for i2=1:4
			for i1=1:4
				ii[cnt] = nodalInd[i,i1]
				jj[cnt] = nodalInd[i,i2]
				# mm[cnt] = M[i1,i2,i]
				cnt    += 1
			end
		end
	end
  mm = vec(M)

	return sparse(ii,jj,mm,nn,nn)

end

# Nodal mass matrix is independent of geometry. Store mass matrix for unit coefficient and unit volume.
const M = (1/20) * [
  2 1 1 1
  1 2 1 1
  1 1 2 1
  1 1 1 2]

function getLocalNodalMassMatrix(sigma::Float64,X::Array{Float64})
	v =  getCellVolume(X[1,:], X[2,:], X[3,:], X[4,:])
	return (v * sigma) * M
end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
using DynamicallyLoadedEmbedding
using Test
using Libdl
# Piggyback on PackageCompiler internals. May have to replicate that
# code if the internals change.
using PackageCompiler: get_compiler, run_with_env, default_sysimg_path

# Compile the C code.
srcdir = joinpath(dirname(@__DIR__), "c")
destdir = mktempdir()
compiler = get_compiler()
if Sys.iswindows()
    binary_path = "$(destdir)/julia_embedding.exe"
    cmd = `$(compiler) $(srcdir)/julia_embedding.c $(srcdir)/julia_cfunctions.c $(srcdir)/main.c -o $(binary_path)`
else
    binary_path = "$(destdir)/julia_embedding"
    cmd = `$(compiler) $(srcdir)/julia_embedding.c $(srcdir)/julia_cfunctions.c $(srcdir)/main.c -ldl -o $(binary_path)`
end
run_with_env(cmd, compiler)
if !isfile(binary_path)
    error("Failed to compile binary.")
end

if Sys.iswindows()
    # Place `julia` in PATH. It is assumed that this makes `libjulia`
    # also available from PATH.
    ENV["PATH"] = string(Sys.BINDIR, ";", ENV["PATH"])
else
    # Make `libjulia` available in LD_LIBRARY_PATH.
    libdir = dirname(abspath(Libdl.dlpath("libjulia")))
    ENV["LD_LIBRARY_PATH"] = string(libdir, ":", get(ENV, "LD_LIBRARY_PATH", ""))
end

cfunctions_path = abspath(dirname(@__DIR__), "src", "julia_cfunctions.jl")

# `VERSION` is usually not the same as this for nightly builds.
# TODO: Find out whether there's some more convenient way to get
#       exactly the same string as this `ccall`.
version_string = unsafe_string(ccall(:jl_ver_string, Ptr{Cchar}, ()))
expected_stdout =
    """
    Julia version $(version_string)
    0.5
    fibonacci(10) = 55
    xi = 2 3 4 4
    xf = 1.0 2.5 4.0
    """

expected_stderr =
    """
    Exception test: MethodError: no method matching this_function_has_no_methods()
    """

@testset "without sysimage" begin
    out = Pipe()
    err = Pipe()
    p = run(pipeline(`$(binary_path) "" $(cfunctions_path)`,
                     stdin=devnull, stdout=out, stderr=err))
    close(out.in)
    close(err.in)
    @test replace(read(out, String), "\r\n" => "\n") == expected_stdout
    @test replace(read(err, String), "\r\n" => "\n") == expected_stderr
end

@testset "with sysimage" begin
    out = Pipe()
    err = Pipe()
    # Use the default system image to test the custom system image
    # code paths.
    p = run(pipeline(`$(binary_path) $(default_sysimg_path()) $(cfunctions_path)`,
                     stdin=devnull, stdout=out, stderr=err))
    close(out.in)
    close(err.in)
    @test replace(read(out, String), "\r\n" => "\n") == expected_stdout
    @test replace(read(err, String), "\r\n" => "\n") == expected_stderr
end


text:
module Day2
export solvePart1, solvePart2

import DelimitedFiles
import Combinatorics

function getInput()
    open("day2input.txt") do f
        return DelimitedFiles.readdlm(f, String)
    end
end

function getLetterCounts(s)
    counts = Dict()
    for c in s
        counts[c] = get(counts, c, 0) + 1
    end
    return counts
end

function numCodesWithNRepeats(n, inputs)
    letters = Char[]
    num = 0
    for code in inputs
        counts = getLetterCounts(code)
        if n in values(counts)
            num +=1
        end
    end
    return num
end

function solvePart1()
    inputs = getInput()
    repeat2 = numCodesWithNRepeats(2, inputs)
    repeat3 = numCodesWithNRepeats(3, inputs)
    return repeat2 * repeat3
end

function numChanges(a, b)
    n = 0
    for (c1, c2) in zip(a, b)
        if c1 != c2
            n +=1
        end
    end
    return n
end

function pairsWithOneChange(inputs)
    pairs = []
    for (s1, s2) in Combinatorics.combinations(inputs, 2)
        if numChanges(s1, s2) == 1
            push!(pairs, (s1, s2))
        end
    end
    return pairs
end

function commonLetters(s1, s2)
    output = []
    for (c1, c2) in zip(s1, s2)
        if c1 == c2
            push!(output, c1)
        end
    end
    return join(output)
end

function solvePart2()
    inputs = getInput()
    oneChange = pairsWithOneChange(inputs)
    commonLetters(oneChange[1]...)
end
end


text:
line = read("input.txt", String)
numbers = parse.(Int, split(line, ','))
counts = zeros(Int, 9)

for n = numbers
    counts[n+1] += 1
end

days = 256
for _ = 1:days
    newfish = counts[1]
    for i = 1:8
        counts[i] = counts[i+1]
    end
    counts[7] += newfish
    counts[9] = newfish
end
println(sum(counts))


text:
# Complementarity constraints (M x + b  complements x)
# As Knitro supports only complementarity constraints between variables,
# we reformulate it by adding auxiliary variables:
#
# x_aux = Mx + b
# (x_aux complements x)
#
function MOI.add_constraint(
    model::Optimizer,
    func::MOI.VectorAffineFunction,
    set::MOI.Complements,
)
    (model.number_solved >= 1) && throw(AddConstraintError())

    # Number of complementarity in Knitro is half the dimension of the MOI set
    n_comp = div(set.dimension, 2)

    # Add auxiliary variables `x_aux`.
    x_aux = KN_add_vars(model.inner, n_comp)

    # Add new constraints in Knitro to formulate x_aux = Mx + b
    n_cons = KN_add_cons(model.inner, n_comp)
    offset = n_cons[1]

    # Parse VectorAffineFunction defining the complementarity constraints
    indexcoords, indexvars, coefs = canonical_vector_affine_reduction(func)

    # Convert VectorAffineFunction in Knitro format
    x_comp = Cint[]
    jac_cons = Cint[]
    jac_vars = Cint[]
    jac_coefs = Cdouble[]
    for (c, v, f) in zip(indexcoords, indexvars, coefs)
        # By convention, if c is greater than the dimension of
        # the complementarity constraint set, then it specifies
        # the complementarity variable `x`.
        if c >= n_comp
            push!(x_comp, v)
        else
            # Otherwise, the index corresponds to one of the matrix `M`.
            # Add it to the Jacobian structure of the problem.
            push!(jac_cons, c + offset)
            push!(jac_vars, v)
            push!(jac_coefs, f)
        end
    end
    # Add the index corresponding to `x_aux` in the Jacobian.
    for (icons, ivar) in zip(n_cons, x_aux)
        push!(jac_cons, icons)
        push!(jac_vars, ivar)
        push!(jac_coefs, -1.0)
    end
    # Get constant structure
    q = func.constants[1:n_comp]

    # Add structure of new constaints to Knitro.
    KN_add_con_linear_struct(model.inner, jac_cons, jac_vars, jac_coefs)
    KN_set_con_eqbnds(model.inner, n_comp, n_cons, -q)

    # Currently, only complementarity constraints between two variables
    # are supported.
    comp_type = fill(KN_CCTYPE_VARVAR, n_comp)

    # Number of complementarity constraint previously added
    n_comp_cons = model.complementarity_cache.n
    _add_complementarity_constraint!(model.complementarity_cache, x_comp, x_aux, comp_type)

    return MOI.ConstraintIndex{typeof(func),typeof(set)}(n_comp_cons)
end

# Complementarity constraints (x_1 complements x_2), with x_1 and x_2
# being two variables of the problem.
function MOI.add_constraint(
    model::Optimizer,
    func::MOI.VectorOfVariables,
    set::MOI.Complements,
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    indv = Cint[v.value - 1 for v in func.variables]

    # Number of complementarity in Knitro is half the dimension of the MOI set
    n_comp = div(set.dimension, 2)
    # Currently, only complementarity constraints between two variables
    # are supported.
    comp_type = fill(KN_CCTYPE_VARVAR, n_comp)

    # Number of complementarity constraint previously added
    n_comp_cons = model.complementarity_cache.n

    _add_complementarity_constraint!(
        model.complementarity_cache,
        indv[1:n_comp],
        indv[n_comp+1:end],
        comp_type,
    )

    return MOI.ConstraintIndex{typeof(func),typeof(set)}(n_comp_cons)
end


text:
export mess2var
function mess2var(vars,FG;ShowFactors=false)
    #MESS2VAR Returns the message numbers that connect into variables v
    # [messnum fact]=mess2var(vars,FG)
    # vars is a vector of variables
    #
    # messnum(i) is the message number from factor fact(i) to vars(i)
    # FG is a Factor Graph
    # see also mess2fact.m, FactorConnectingVariable.m, FactorGraph.m

    V=minimum(find(FG[1,:]))-1; # variables are first in the order
    messnum=[]; fact=[];
    for v=vars
        tmp=vec(full(FG[:,v]));
        tmp2=tmp.>0;
        messnum = vcat(messnum,tmp[tmp2]);
        fact=vcat(fact, find(tmp2)-V);
    end
    if ShowFactors
        return messnum, fact
    else
        return messnum
    end
end

export sumprod
function sumprod(pot,A;InitialMessages=[])

    variables=potvariables(pot)
    V=length(variables); N=size(A,1)
    fnodes=zeros(Int64,1,N); fnodes[V+1:end]=1:N-V; # factor nodes
    vnodes=zeros(Int64,1,N); vnodes[1:V]=1:V; # variable nodes
    nmess=maxarray(A);
    mess=Array(Potential,nmess)
    InitialMessages=[]
    if !isempty(InitialMessages); mess=InitialMessages; end
    if isempty(InitialMessages) # message initialisation
        for mcount=1:nmess
            mess[mcount]=Const(1)
            FGnodeA, FGnodeB=findone(A,mcount);
            if fnodes[FGnodeA].>0 # factor to variable message:
                # if the factor is at the edge (simplical), need to set message to the factor potentials
                if length(find(A[FGnodeA,:]))==1
                    mess[mcount]=pot[fnodes[FGnodeA]];
                end
            end
        end
    end

 # Do the message passing:
    for mcount=1:length(mess)
        FGnodeA, FGnodeB=findone(A,mcount);
        FGparents=setdiff(find(A[FGnodeA,:]),FGnodeB) # FG parent nodes of FGnodeA
        if !isempty(FGparents)
#            tmpmess = prod(mess[A[FGparents,FGnodeA]]) # product of incoming messages
            # At the moment (0.3.0 rc2) Julia assumes prod(single element) is an error. I don't know how to dispatch on the size of an input, so unfortunately need to define my own prod function:
            tmpmess = multpots(mess[A[FGparents,FGnodeA]]) # product of incoming messages
            mfactor=fnodes[FGnodeA]
            if mfactor==0 # variable to factor message:
                mess[mcount]=tmpmess
            else # factor to variable message:
                tmpmess = multpots([tmpmess pot[mfactor]])
                mess[mcount] = sum(tmpmess,FGnodeB,SumOver=false)
            end
        end
    end
    # Get all the marginals: variable nodes are first in the ordering, so
    marg=Array(Potential,V)
    for i=1:V
        dum1, dum2, incoming=findnz(A[:,i]');
        tmpmess = multpots(mess[incoming]);
        marg[i]=tmpmess;
    end
    normconstpot=sum(multpots(mess[mess2var(1,A)]));
    return marg, mess, normconstpot
end

export maxprod
function maxprod(pot,A;InitialMessages=[])

    variables=potvariables(pot)
    V=length(variables); N=size(A,1)
    fnodes=zeros(Int64,1,N); fnodes[V+1:end]=1:N-V; # factor nodes
    vnodes=zeros(Int64,1,N); vnodes[1:V]=1:V; # variable nodes
    nmess=maxarray(A);
    mess=Array(Potential,nmess)
    InitialMessages=[]
    if !isempty(InitialMessages); mess=InitialMessages; end
    if isempty(InitialMessages) # message initialisation
        for mcount=1:nmess
            mess[mcount]=Const(1)
            FGnodeA, FGnodeB=findone(A,mcount);
            if fnodes[FGnodeA].>0 # factor to variable message:
                # if the factor is at the edge (simplical), need to set message to the factor potentials
                if length(find(A[FGnodeA,:]))==1
                    mess[mcount]=pot[fnodes[FGnodeA]];
                end
            end
        end
    end

 # Do the message passing:
    for mcount=1:length(mess)
        FGnodeA, FGnodeB=findone(A,mcount);
        FGparents=setdiff(find(A[FGnodeA,:]),FGnodeB) # FG parent nodes of FGnodeA
        if !isempty(FGparents)
            tmpmess = multpots(mess[A[FGparents,FGnodeA]]) # product of incoming messages
            mfactor=fnodes[FGnodeA]
            if mfactor==0 # variable to factor message:
                mess[mcount]=tmpmess
            else # factor to variable message:
                tmpmess = multpots([tmpmess pot[mfactor]])
                mess[mcount] = max(tmpmess,FGnodeB,MaxOver=false)
            end
        end
    end

    # now find the maximum states: variable nodes are first in the ordering:
    maxstate=Dict{Integer,Any}()
    maxval=[] # ensure that is in scope
    for i=1:V
        #dum1, dum2, incoming=findnz(A[:,i])
        dum1, incoming=findnz(A[:,i])
        tmpmess = multpots(mess[incoming])
        tmpmess,maxstate[i]=max(tmpmess,i,ReturnState=true)
        maxval=tmpmess.content
    end

    return maxstate, maxval, mess
end

#function [f fact2varmess var2factmess] = FactorConnectingVariable(vars,A)
export FactorConnectingVariable
function FactorConnectingVariable(vars,A)
    #%FACTORCONNECTINGVARIABLE Factor nodes connecting to a set of variables
    #% [f fact2varmess var2factmess] = FactorConnectingVariable(vars,A)
    #% find the intersection of factor indices that connect to variables vars.
    #% fact2varmess are the message indices connecting the factors f to the variable
    #% var2factmess are the message indices from the variables to the factors f
    #% A is a Factor Graph
    #% vars are the variables
    #% see also FactorGraph.m, demoSumProd.m, mess2fact.m, fact2mess.m
    c=1
    for v=vars
        if c==1
            f = find(A[v,:]);
        else
            f = intersect(f,find(A[v,:]));
        end
        c=c+1
    end
    fact2varmess=[];var2factmess=[];
    for fn=f
        fact2varmess=full(union(fact2varmess,A[fn,vars]));
        var2factmess=full(union(var2factmess,A[vars,fn]));
    end
    V=minimum(find(A[1,:]))-1; # variables are first in the order
    f=full(f-V)
    return f, fact2varmess, var2factmess
end

export FactorGraph
function FactorGraph(pot::DictOrArray)
    #FACTORGRAPH Returns a Factor Graph adjacency matrix based on a cell of potentials.
    # A = FactorGraph(pot)
    #
    # pot is a set of potentials and the routine returns the Factor Graph
    # (sparse) adjacency matrix A. The size of A is equal to
    # (V+F)*(V+F) where V are the total number of variables and F the total
    # number of factors. A(1:V,1:V) and A(V+1:end,V+1:end) are empty;
    # A(1:V,V+1:end) contains the variable to factor message indices and
    # A(V+1:end,1:V) contains the factor to variable message indices
    # If the set of potentials is not singly-connected, all message indices are -1
    #
    #
    # A(i,j)=k means that message number k is from FGnodei -> FGnodej
    # Going through the messages in sequence corresponds to a valid
    # forward-backward procedure over all variable nodes in the Factor Graph.
    # See also FactorConnectingVariable.m, VariableConnectingFactor.m
    # Note that the variables in pot must be numbered 1,2,3,...
    #
    # See also demoSumProd.jl

    F=length(pot) # number of factors (potentials in distribution)
    variables=potvariables(pot)
    if variables[end]!=length(variables);
        error("potential variables used are not numbered 1:end. Use standarisevariables.jl. See demoSumProdStandariseVariables"); end

    V=length(variables)
    N=V+F # all nodes in factor graph

    vnodes=zeros(Int64,1,N); vnodes[1:V]=1:V # variable nodes
    fnodes=zeros(Int64,1,N); fnodes[V+1:end]=1:F # factor nodes
    A = SparseIntMatrix(N,N)
    for f=1:length(pot)
        FGnodeA=find(fnodes.==f)
        FGnodeB=vec(memberinds(pot[f].variables,variables))
        fillmatrix!(A,FGnodeA,FGnodeB[:],1)
        fillmatrix!(A,FGnodeB[:],FGnodeA,1)
    end
    # get a message passing sequence and initialise the messages
    #tree, elimseq, forwardschedule=istree(full(A),ReturnElimSeq=true); ## change to sparse
    tree, elimseq, forwardschedule=istree(A,ReturnElimSeq=true); ## change to sparse
    #reverseschedule=flipud(fliplr(forwardschedule));
    reverseschedule=flipdim(flipdim(forwardschedule,2),1);
    schedule=vcat(forwardschedule,reverseschedule);

    if tree
        for count=1:size(schedule,1)
            # setup the structure for a message from FGnodeA -> FGnodeB
            FGnodeA, FGnodeB = schedule[count,:]
            A[FGnodeA,FGnodeB]=count; # store the FG adjacency matrix with mess number on edge
        end
    else
        A = replace(A,1,-1);
    end

    return A

end


text:
using Revise

using BenchmarkTools
using Clustering
using ISODATA
using MLDataUtils
using MLLabelUtils
using StatsBase
using AssignmentSolver

X, Y, fnames = load_iris()
Y = convertlabel(MLLabelUtils.LabelEnc.Indices{Int64,3}(), Y, labelenc(Y))
# X = zscore(X)
# res = base(X, 3, θn=0.25, θe=0.12, θc=1.5, iter=15)
# println(unique(res.assignments))
# println(randindex(res.assignments, Y))

function cluster_acc(y_true, y_pred)
    w = counts(y_true, y_pred)
    cost = maximum(w) .- w
    reward = Float64.(reward2cost(cost))
    auctionsol, lambda = auction_assignment(reward)
    compute_objective(auctionsol, cost)
    ind = auctionsol.r2c

    n = size(ind, 1)
    x = 0
    for i in 1:n
        x += w[i, ind[i]]
    end

    return (x * 1.0/(size(y_true, 1))) * 100.0
end
a = [0, 0, 1, 0, 2, 2]
b = [1, 1, 2, 2, 0, 0]
print(cluster_acc(a, b))
#  print(cluster_acc(Y, res.assignments))
# runs = 1000
# results = []
# for i in 1:runs
#     res = base(X, 3, θn=0.25, θe=0.12, θc=1.5, iter=15)
#     acc = cluster_acc(Y, res.assignments)
#     push!(results, acc)
#     #push!(results, randindex(res.assignments, Y)[1])
# end

# avg_accuracy = mean(results)
# min_accuracy = minimum(results)
# max_accuracy = maximum(results)
# println("ISODATA:")
# println(avg_accuracy)
# println(min_accuracy)
# println(max_accuracy)
# println("Average number of clusters: ")

# results = []
# for i in 1:runs
#     res = kmeans(X, 3, init=:rand)
#     # res = kmeans(X, 3)
#     acc = cluster_acc(Y, res.assignments)
#     push!(results, acc)
#     # push!(results, randindex(res.assignments, Y)[2])
# end
# avg_accuracy = mean(results)
# min_accuracy = minimum(results)
# max_accuracy = maximum(results)
# println("Kmeans:")
# println(avg_accuracy)
# println(min_accuracy)
# println(max_accuracy)