text:
    using SACFiles
using SACFiles.parsetext

@testset "IO" begin
    @testset "Binary Data Decoding" begin
        @testset "Float32" begin
            let testnums = Float32[-12345:12345;]
                decnums = SACFiles.decodesacbytes(Float32, reinterpret(UInt8, testnums))
                @test isapprox(testnums, decnums)

                decnums = SACFiles.decodesacbytes(Float32, reinterpret(UInt8, map(bswap, testnums)), true)
                @test isapprox(testnums, decnums)
            end
        end

        @testset "Int32" begin
            let testnums = Int32[-12345:12345;]
                decnums = SACFiles.decodesacbytes(Int32, reinterpret(UInt8, testnums))
                @test testnums == decnums

                decnums = SACFiles.decodesacbytes(Int32, reinterpret(UInt8, map(bswap, testnums)), true)
                @test testnums == decnums
            end
        end

        @testset "Bool" begin
            # Testing Bools which are just encoded as Int32s
            let testnums = rand(Int32[1, 0], 1000)
                decnums = SACFiles.decodesacbytes(Bool, reinterpret(UInt8, testnums))
                @test testnums == decnums

                decnums = SACFiles.decodesacbytes(Bool, reinterpret(UInt8, map(bswap, testnums)), true)
                @test testnums == decnums

                # Test that -12345 -> true, or anything else for that matter
                # (because that's what SAC does)
                testnums = Int32[1, 0, -12345]
                decnums = SACFiles.decodesacbytes(Bool, reinterpret(UInt8, testnums))
                @test [true, false, true] == decnums

                decnums = SACFiles.decodesacbytes(Bool, reinterpret(UInt8, map(bswap, testnums)), true)
                @test [true, false, true] == decnums
            end
        end

        @testset "Enums" begin
            let testnums = collect(instances(HeaderEnum))
                decnums = SACFiles.decodesacbytes(HeaderEnum, reinterpret(UInt8, testnums))
                @test testnums == decnums

                decnums = SACFiles.decodesacbytes(HeaderEnum,
                                                  reinterpret(UInt8, map((e)->bswap(Int32(e)), testnums)),
                                                  true)
                @test testnums == decnums
            end
        end

        @testset "Strings" begin
            let teststrings = String["hello", "world"]
                # Test 2 word strings
                map!((s) -> rpad(s, 8, ' '), teststrings)
                bs = convert(Array{UInt8}, string(teststrings...))
                decstrings = SACFiles.decodesacbytes(String, bs, 8)
                @test teststrings == decstrings

                # Test 4 word strings
                map!((s) -> rpad(s, 16, ' '), teststrings)
                bs = convert(Array{UInt8}, string(teststrings...))
                decstrings = SACFiles.decodesacbytes(String, bs, 16)
                @test teststrings == decstrings

                # Test decoding into a single string
                decstrings = SACFiles.decodesacbytes(String, bs, 32)
                @test string(teststrings...) == decstrings[1]
            end
        end

        @testset "Endianness Tests" begin
            let endianness = ENDIAN_BOM == 0x04030201 ? :le : :be
                if endianness == :le
                    open("./test-files/test-seismo.sac") do f
                        @test SACFiles.isalienend(f) == false
                    end
                    open("./test-files/be-test-seismo.sac") do f
                        @test SACFiles.isalienend(f) == true
                    end
                else
                    open("./test-files/test-seismo.sac") do f
                        @test SACFiles.isalienend(f) == true
                    end
                    open("./test-files/be-test-seismo.sac") do f
                        @test SACFiles.isalienend(f) == false
                    end
                end
            end
            open("./test-files/test-seismo.sac") do f
                seek(f, 42)
                SACFiles.isalienend(f)
                @test position(f) == 42
            end
        end
    end

    @testset "ASCII Text Parsing" begin
        @testset "Float32" begin
            testnums = Float32[-12345, 12345, 3.14159, 2.71]
            teststring = @sprintf("%15.7f%15.7f\n%15.7f%15.7f", testnums...)
            decnums = parsetext(Float32, teststring)
            @test isapprox(testnums, decnums)
        end

        @testset "Int32" begin
            testnums = Int32[-12345, 12345, 42, 1337]
            teststring = @sprintf("%10d%10d\n%10d%10d", testnums...)
            decnums = parsetext(Int32, teststring)
            @test testnums == decnums
        end

        @testset "Bool" begin
            testnums = [true, false, true, false]
            teststring = @sprintf("%10d%10d\n%10d%10d", testnums...)
            decnums = parsetext(Bool, teststring)
            @test testnums == decnums

            testnums = Int32[1, 0, 1, -12345]
            teststring = @sprintf("%10d%10d\n%10d%10d", testnums...)
            decnums = parsetext(Bool, teststring)
            @test [true, false, true, true] == decnums
        end

        @testset "Enums" begin
            testnums = [SACFiles.itime, SACFiles.iday, SACFiles.icaltech,
                        SACFiles.undefined, SACFiles.ios]
            teststring = @sprintf("%10d%10d%10d\n%10d%10d", reinterpret(Int32, testnums)...)
            decnums = parsetext(HeaderEnum, teststring)
            @test testnums == decnums
        end

        @testset "Strings" begin
            teststrings = ["hello", "world", "this", "is", "a", "test!"]
            teststring = @sprintf("%8s%8s\n%8s%8s%8s%8s", teststrings...)
            decstrings = parsetext(String, teststring, 8)
            @test teststrings == decstrings
        end
    end

@testset "Header File Reading" begin
    @testset "Binary Files" begin
        # Make sure we're reading the whole header
        open("./test-files/test-seismo.sac", "r") do f
            _ = readsachdr(f)
            @test position(f) == 632
        end

        # Test header reading on a hexed file.
        open("./test-files/test-hexed-header.sac", "r") do f
            inhdr = readsachdr(f)
            hdr_hexedvalues = Dict{Type, Any}(Float32     => Float32(1337.0),
                                              Int32       => Int32(1337),
                                              HeaderEnum  => SACFiles.inucl,
                                              String => ascii("BLEEPBLO"),
                                              Bool        => true)
            testhdr = Header()

            map(fieldnames(Header), Header.types) do field, T
                if field == :kevnm
                    setfield!(testhdr, field, ascii("BLEEPBLOOPBLEEPS"))
                elseif field == :nvhdr
                    setfield!(testhdr, field, Int32(6))
                else
                    setfield!(testhdr, field, hdr_hexedvalues[T])
                end
                @test getfield(inhdr, field) == getfield(testhdr, field)
            end
        end

        # Test header reading functions on a synthetic SAC generated seismogram
        testhdr = SACTestUtilities.make_test_seismo_hdr()
        open("./test-files/test-seismo.sac", "r") do f
            inhdr = readsachdr(f)
            SACTestUtilities.testhdrequal(inhdr, testhdr)
        end
        SACTestUtilities.testhdrequal(readsachdr("./test-files/test-seismo.sac"), testhdr)
        SACTestUtilities.testhdrequal(readsachdr("./test-files/be-test-seismo.sac"), testhdr)
    end

    @testset "ASCII Files" begin
        testhdr = SACTestUtilities.make_test_seismo_hdr()
        open("./test-files/test-seismo.txt") do f
            inhdr = readsachdr(f, ascii=true)
            @test position(f) == 1672
            SACTestUtilities.testhdrequal(inhdr, testhdr)
        end
        SACTestUtilities.testhdrequal(readsachdr("./test-files/test-seismo.txt", ascii=true), testhdr)
    end
end

@testset "Data File Reading" begin
    field_mappings = Dict(EvenTimeSeries => Dict("data1" => :data),
                          UnevenTimeSeries => Dict("data1" => :ddata, "data2" => :idata),
                          AmplitudeSpectrum => Dict("data1" => :ampdata, "data2" => :phasedata),
                          ComplexSpectrum => Dict("data1" => :data),
                          GeneralXY => Dict("data1" => :y, "data2" => :x))

    for (T, files) in testfile_specs
        data1field = field_mappings[T]["data1"]
        data2field = get(field_mappings[T], "data2", :none)

        for file in files
            # Test reading from file name
            readdata = readsac(file["fname"], ascii=file["ascii"])
            @test isapprox(getfield(readdata, data1field), file["data1"])
            if data2field != :none && haskey(file, "data2")
                @test isapprox(getfield(readdata, data2field), file["data2"])
            end

            # Test reading from IOStream
            open(file["fname"]) do f
                # Non-type stable
                readdata = readsac(f, ascii=file["ascii"])
                @test isapprox(getfield(readdata, data1field), file["data1"])
                if data2field != :none && haskey(file, "data2")
                    @test isapprox(getfield(readdata, data2field), file["data2"])
                end

                # Type stable
                readdata = readsac(T, f, ascii=file["ascii"])
                @test isapprox(getfield(readdata, data1field), file["data1"])
                if data2field != :none && haskey(file, "data2")
                    @test isapprox(getfield(readdata, data2field), file["data2"])
                end
            end
        end
    end
end
end


text:
using Pathfinder
using Random
using Test

Random.seed!(0)

@testset "Pathfinder.jl" begin
    include("transducers.jl")
    include("woodbury.jl")
    include("optimize.jl")
    include("inverse_hessian.jl")
    include("mvnormal.jl")
    include("elbo.jl")
    include("resample.jl")
    include("singlepath.jl")
    include("multipath.jl")
end

ext: jl


text:
@testset "system" begin
	dims = (32, 1, 1)
	cell = UnitCell(1, :lamellar, 10.0)

	amon = Monomer(; id = 1, vol = 1.0, charge = 0.0)
	bmon = Monomer(; id = 2, vol = 1.0, charge = 0.0)
	cmon = Monomer(; id = 3, vol = 1.0, charge = 0.0)
	point = Point(cmon)
	chain = Diblock(amon, bmon, 100, 0.5, 1.0, 1.0, 100)

	itx = FloryInteraction()
	set_interaction!(itx, amon.id, bmon.id, 0.2)

	sys = FieldSystem(dims, cell; monomers = [amon, bmon, cmon], ensemble = Canonical)
	add_species!(sys, point, 0.75)
	add_species!(sys, chain, 0.25)
	add_interaction!(sys, itx)

	# Checks many initialization fields
	@test_throws ErrorException validate(sys)

	# Check actual system usage
	fieldinit!(sys; scale = 0.1, seed = 12345)

	# Species density operators
	density!(point)
	density!(chain)

	@test point.Q ≈ 1.004177455788938
	@test chain.Q ≈ 8.24300700515971

	@test free_energy(sys) ≈ -0.1681580628744596
	@test free_energy_bulk(sys) ≈ -0.9686022902416352
	@test scfstress(sys)[1] ≈ -0.0004272750367130805
end


text:
## SCENARIO PARTITION

function memoryPart(Vec_old, duals, d; memory = 0.)
    # If memory = 1: Vec_old = duals_old
    # If memory = Inf: Vec_old = Ind_old
    if memory == "0"
        Sim = calcSim(duals, d, S)
        Ind_new = AffProp(Sim)

    else
        if memory == "1"
            # Calculating partitions for the duals on the previos iteration
            Sim_old = calcSim(Vec_old, d, S)
            # Calculating new subpartition indices
            Vec_old = AffProp(Sim_old)
        end
        DualPart = Array{Float64}[]
            locsPart = Array{Float64}[]
            N = size(duals,1)
            NT = size(duals,2)
            NP = maximum(Vec_old)
            Ind_new = zeros(N)
            NsubParts = 1

        Scenarios = [duals Vec_old collect(1:N)]

        for p in 1:NP
            DualPart = Scenarios[Scenarios[:,NT+1].==p,:] # Accessing the Vec_old column
                sizePart = size(DualPart,1)
                locsPart = Integer.(DualPart[:,end])# indices of scenarios in subpart

            if sizePart > 1
                # Calculating Distance within partition
                Sim = calcSim(DualPart[:,1:NT], d, 1:sizePart)
                # Calculating new subpartition indices
                clustInds = AffProp(Sim)

                # Assigning new subpart indices pre-clustering
                Ind_new[locsPart] .= NsubParts
                # Reference value of subpartition/ first value
                ref = clustInds[1]

                for n in 2:sizePart
                    # Comparing sub-assignments
                    if clustInds[n] != ref
                        NsubParts += 1
                        ref = clustInds[n] # Reference (1st) value for the partition
                    end
                    Ind_new[locsPart[n]] = NsubParts
                end
            else
                Ind_new[locsPart] .= NsubParts
            end
            NsubParts += 1
        end
    end
    return Integer.(Ind_new)
end

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
const INT_MAX = typemax(Cint)

ZSTD_EXPAND_AND_QUOTE(expr) = string(expr)


text:
using Distributed
using Test
using ParallelUtilities
import ParallelUtilities: ProductSplit, ProductSection, ZipSplit, zipsplit,
minimumelement, maximumelement, extremaelement, nelements, dropleading, indexinproduct,
extremadims, localindex, extrema_commonlastdim, whichproc, procrange_recast, whichproc_localindex,
getiterators, _niterators
using SplittablesBase

macro testsetwithinfo(str, ex)
    quote
        @info "Testing "*$str
        @testset $str begin $(esc(ex)); end;
    end
end

@testsetwithinfo "AbstractConstrainedProduct" begin

    various_iters = Any[(1:10,), (1:1:10,), (1:10, 4:6), (1:1:10, 4:6), (1:10, 4:6, 1:4), (1:2:9,), (1:2:9, 4:1:6),
                    (1:2, Base.OneTo(4), 1:3:10), (1:0.5:3, 2:4)]

    @testsetwithinfo "ProductSplit" begin

        function split_across_processors_iterators(arr::Iterators.ProductIterator, num_procs, proc_id)

            num_tasks = length(arr);

            num_tasks_per_process, num_tasks_leftover = divrem(num_tasks, num_procs)

            num_tasks_on_proc = num_tasks_per_process + (proc_id <= mod(num_tasks, num_procs) ? 1 : 0 );
            task_start = num_tasks_per_process*(proc_id-1) + min(num_tasks_leftover, proc_id-1) + 1;

            Iterators.take(Iterators.drop(arr, task_start-1), num_tasks_on_proc)
        end

        function split_product_across_processors_iterators(arrs_tuple, num_procs, proc_id)
            split_across_processors_iterators(Iterators.product(arrs_tuple...), num_procs, proc_id)
        end

        @testset "Constructor" begin

            function checkPSconstructor(iters, npmax = 10)
                ntasks_total = prod(length, iters)
                for np = 1:npmax, p = 1:np
                    ps = ProductSplit(iters, np, p)
                    @test eltype(ps) == Tuple{map(eltype, iters)...}
                    @test _niterators(ps) == length(iters)
                    if !isempty(ps)
                        @test collect(ps) == collect(split_product_across_processors_iterators(iters, np, p))
                    end
                    @test prod(length, getiterators(ps)) == ntasks_total
                    @test ParallelUtilities.workerrank(ps) == p
                    @test nworkers(ps) == np
                end

                @test_throws ArgumentError ProductSplit(iters, npmax, npmax + 1)
            end

            @testset "0D" begin
                @test_throws ArgumentError ProductSplit((), 2, 1)
            end

            @testset "cumprod" begin
                @test ParallelUtilities._cumprod(1,()) == ()
                @test ParallelUtilities._cumprod(1,(2,)) == (1,)
                @test ParallelUtilities._cumprod(1,(2, 3)) == (1, 2)
                @test ParallelUtilities._cumprod(1,(2, 3, 4)) == (1, 2, 6)
            end

            @testset "1D" begin
                iters = (1:10,)
                checkPSconstructor(iters)
            end
            @testset "2D" begin
                iters = (1:10, 4:6)
                checkPSconstructor(iters)
            end
            @testset "3D" begin
                iters = (1:10, 4:6, 1:4)
                checkPSconstructor(iters)
            end
            @testset "steps" begin
                iters = (1:2:10, 4:1:6)
                checkPSconstructor(iters)
            end
            @testset "mixed" begin
                for iters in [(1:2, 4:2:6), (1:2, Base.OneTo(4), 1:3:10)]
                    checkPSconstructor(iters)
                end
            end

            @testset "empty" begin
                iters = (1:1,)
                ps = ProductSplit(iters, 10, 2)
                @test isempty(ps)
                @test length(ps) == 0
            end

            @testset "first and last ind" begin
                for iters in Any[(1:10,), (1:2, Base.OneTo(4), 1:3:10)]
                    ps = ProductSplit(iters, 2, 1)
                    @test firstindex(ps) == 1
                    @test ParallelUtilities.firstindexglobal(ps) == 1
                    @test ParallelUtilities.lastindexglobal(ps) == div(prod(length, iters), 2)
                    @test lastindex(ps) == div(prod(length, iters), 2)
                    @test lastindex(ps) == length(ps)
                    ps = ProductSplit(iters, 2, 2)
                    @test ParallelUtilities.firstindexglobal(ps) == div(prod(length, iters), 2) + 1
                    @test firstindex(ps) == 1
                    @test ParallelUtilities.lastindexglobal(ps) == prod(length, iters)
                    @test lastindex(ps) == length(ps)

                    for np in prod(length, iters) + 1:prod(length, iters) + 10,
                        p in prod(length, iters) + 1:np

                        ps = ProductSplit(iters, np, p)
                        @test ParallelUtilities.firstindexglobal(ps) == prod(length, iters) + 1
                        @test ParallelUtilities.lastindexglobal(ps) == prod(length, iters)
                    end
                end
            end
        end

        @testset "firstlast" begin
            @testset "first" begin

                @test ParallelUtilities._first(()) == ()

                for iters in various_iters, np = 1:prod(length, iters)

                    ps = ProductSplit(iters, np, 1)
                    @test first(ps) == map(first, iters)
                end
            end
            @testset "last" begin

                @test ParallelUtilities._last(()) == ()

                for iters in various_iters, np = 1:prod(length, iters)

                    ps = ProductSplit(iters, np, np)
                    @test last(ps) == map(last, iters)
                end
            end
        end

        @testset "extrema" begin

            @testset "min max extrema" begin
                function checkPSextrema(iters, (fn_el, fn), npmax = 10)
                    for np = 1:npmax, p = 1:np
                        ps = ProductSplit(iters, np, p)
                        if isempty(ps)
                            continue
                        end
                        pcol = collect(ps)
                        for dims in 1:length(iters)
                            @test begin
                                res = fn_el(ps, dims = dims) == fn(x[dims] for x in pcol)
                                if !res
                                    show(ps)
                                end
                                res
                            end
                        end
                        if _niterators(ps) == 1
                            @test begin
                                res = fn_el(ps) == fn(x[1] for x in pcol)
                                if !res
                                    show(ps)
                                end
                                res
                            end
                        end
                    end
                end

                for iters in various_iters,
                    fntup in [(maximumelement, maximum), (minimumelement, minimum), (extremaelement, extrema)]
                    checkPSextrema(iters, fntup)
                end

                @test minimumelement(ProductSplit((1:5,), 2, 1)) == 1
                @test maximumelement(ProductSplit((1:5,), 2, 1)) == 3
                @test extremaelement(ProductSplit((1:5,), 2, 1)) == (1, 3)

                @test minimumelement(ProductSplit((1:5,), 2, 2)) == 4
                @test maximumelement(ProductSplit((1:5,), 2, 2)) == 5
                @test extremaelement(ProductSplit((1:5,), 2, 2)) == (4, 5)
            end

            @testset "extremadims" begin
                ps = ProductSplit((1:10,), 2, 1)
                @test ParallelUtilities._extremadims(ps, 1,()) == ()
                for iters in various_iters
                    dims = length(iters)
                    for np = 1:prod(length, iters) + 1, proc_id = 1:np
                        ps = ProductSplit(iters, np, proc_id)
                        if isempty(ps)
                            @test_throws ArgumentError extremadims(ps)
                        else
                            ext = Tuple(map(extrema, zip(collect(ps)...)))
                            @test extremadims(ps) == ext
                        end
                    end
                end
            end

            @testset "extrema_commonlastdim" begin
                iters = (1:10, 4:6, 1:4)
                ps = ProductSplit(iters, 37, 8)
                @test extrema_commonlastdim(ps) == ([(9, 1), (6, 1)], [(2, 2), (4, 2)])
                ps = ProductSplit(iters, prod(length, iters) + 1, prod(length, iters) + 1)
                @test extrema_commonlastdim(ps) === nothing
            end
        end

        @testset "in" begin

            function checkifpresent(iters, npmax = 10)
                for np = 1:npmax, p = 1:np
                    ps = ProductSplit(iters, np, p)
                    if isempty(ps)
                        continue
                    end
                    pcol = collect(ps)

                    for el in pcol
                        # It should be contained in this iterator
                        @test el in ps
                        for p2 in 1:np
                            # It should not be contained anywhere else
                            p2 == p && continue
                            ps2 = ProductSplit(iters, np, p2)
                            @test !(el in ps2)
                        end
                    end
                end
            end

            for iters in various_iters
                checkifpresent(iters)
            end

            @test ParallelUtilities._infullrange((), ())
        end

        @testset "whichproc + procrange_recast" begin
            np, proc_id = 5, 5
            iters = (1:10, 4:6, 1:4)
            ps = ProductSplit(iters, np, proc_id)
            @test whichproc(iters, first(ps), 1) == 1
            @test whichproc(ps, first(ps)) == proc_id
            @test whichproc(ps, last(ps)) == proc_id
            @test whichproc(iters,(100, 100, 100), 1) === nothing
            @test procrange_recast(iters, ps, 1) == 1:1
            @test procrange_recast(ps, 1) == 1:1

            smalleriter = (1:1, 1:1, 1:1)
            err = ParallelUtilities.TaskNotPresentError(smalleriter, first(ps))
            @test_throws err procrange_recast(smalleriter, ps, 1)
            smalleriter = (7:9, 4:6, 1:4)
            err = ParallelUtilities.TaskNotPresentError(smalleriter, last(ps))
            @test_throws err procrange_recast(smalleriter, ps, 1)

            iters = (1:1, 2:2)
            ps = ProductSplit(iters, np, proc_id)
            @test procrange_recast(iters, ps, 2) == nothing
            @test procrange_recast(ps, 2) == nothing

            iters = (1:1, 2:2)
            ps = ProductSplit(iters, 1, 1)
            @test procrange_recast(iters, ps, 2) == 1:1
            @test procrange_recast(ps, 2) == 1:1

            iters = (Base.OneTo(2), 2:4)
            ps = ProductSplit(iters, 2, 1)
            @test procrange_recast(iters, ps, 1) == 1:1
            @test procrange_recast(iters, ps, 2) == 1:1
            @test procrange_recast(iters, ps, prod(length, iters)) == 1:length(ps)

            for np_new in 1:prod(length, iters)
                for proc_id_new = 1:np_new
                    ps_new = ProductSplit(iters, np_new, proc_id_new)

                    for val in ps_new
                        # Should loop only if ps_new is non-empty
                        @test whichproc(iters, val, np_new) == proc_id_new
                    end
                end
                @test procrange_recast(iters, ps, np_new) == (isempty(ps) ? nothing : (whichproc(iters, first(ps), np_new):whichproc(iters, last(ps), np_new)))
                @test procrange_recast(ps, np_new) == (isempty(ps) ? nothing : (whichproc(iters, first(ps), np_new):whichproc(iters, last(ps), np_new)))
            end

            @testset "different set" begin
                iters = (1:100, 1:4000)
                ps = ProductSplit((20:30, 1:1), 2, 1)
                @test procrange_recast(iters, ps, 700) == 1:1
                ps = ProductSplit((20:30, 1:1), 2, 2)
                @test procrange_recast(iters, ps, 700) == 1:1

                iters = (1:1, 2:2)
                ps = ProductSplit((20:30, 2:2), 2, 1)
                @test_throws ParallelUtilities.TaskNotPresentError procrange_recast(iters, ps, 3)
                ps = ProductSplit((1:30, 2:2), 2, 1)
                @test_throws ParallelUtilities.TaskNotPresentError procrange_recast(iters, ps, 3)
            end
        end

        @testset "indexinproduct" begin
            @test indexinproduct((1:4, 2:3:8), (3, 5)) == 7
            @test indexinproduct((1:4, 2:3:8), (3, 6)) === nothing
            @test_throws ArgumentError indexinproduct((), ())
        end

        @testset "localindex" begin
            for iters in various_iters
                for np = 1:prod(length, iters), proc_id = 1:np
                    ps = ProductSplit(iters, np, proc_id)
                    for (ind, val) in enumerate(ps)
                        @test localindex(ps, val) == ind
                    end
                end
            end
        end

        @testset "whichproc_localindex" begin
            for iters in various_iters
                iters isa Tuple{AbstractUnitRange, Vararg{AbstractUnitRange}} || continue
                for np = 1:prod(length, iters), proc_id = 1:np
                    ps_col = collect(ProductSplit(iters, np, proc_id))
                    ps_col_rev = [reverse(t) for t in ps_col]
                    for val in ps_col
                        p, ind = whichproc_localindex(iters, val, np)
                        @test p == proc_id
                        ind_in_arr = searchsortedfirst(ps_col_rev, reverse(val))
                        @test ind == ind_in_arr
                    end
                end
            end
            @test whichproc_localindex((1:1,1:1), (1,2), 1) === nothing
        end

        @testset "getindex" begin

            @test ParallelUtilities._getindex((), 1) == ()
            @test ParallelUtilities._getindex((), 1, 2) == ()

            @test ParallelUtilities.childindex((), 1) == (1,)

            for iters in various_iters
                for np = 1:prod(length, iters), p = 1:np
                    ps = ProductSplit(iters, np, p)
                    ps_col = collect(ps)
                    for i in 1:length(ps)
                        @test ps[i] == ps_col[i]
                    end
                    @test ps[end] == ps[length(ps)]
                    for ind in [0, length(ps) + 1]
                       @test_throws ParallelUtilities.BoundsError(ps, ind) ps[ind]
                    end
                end
            end
        end
    end
    @testsetwithinfo "ProductSection" begin
        @testset "Constructor" begin
            function testPS(iterators)
                itp = collect(Iterators.product(iterators...))
                l = length(itp)
                for startind in 1:l, endind in startind:l
                    ps = ProductSection(iterators, startind:endind)
                    @test eltype(ps) == Tuple{map(eltype, iterators)...}
                    for (psind, ind) in enumerate(startind:endind)
                        @test ps[psind] == itp[ind]
                    end
                end
            end

            for iter in various_iters
                testPS(iter)
            end

            @test_throws ArgumentError ProductSection((), 2:3)
        end
    end
    @testset "dropleading" begin
        ps = ProductSplit((1:5, 2:4, 1:3), 7, 3);
        @test dropleading(ps) isa ProductSection
        @test collect(dropleading(ps)) == [(4, 1), (2, 2), (3, 2)]
        @test collect(dropleading(dropleading(ps))) == [(1,), (2,)]

        ps = ProductSection((1:5, 2:4, 1:3), 5:8);
        @test dropleading(ps) isa ProductSection
        @test collect(dropleading(ps)) == [(2, 1), (3, 1)]
        @test collect(dropleading(dropleading(ps))) == [(1,)]
    end
    @testset "nelements" begin
        ps = ProductSplit((1:5, 2:4, 1:3), 7, 3);
        @test nelements(ps, dims = 1) == 5
        @test nelements(ps, dims = 2) == 3
        @test nelements(ps, dims = 3) == 2
        @test_throws ArgumentError nelements(ps, dims = 0)
        @test_throws ArgumentError nelements(ps, dims = 4)

        ps = ProductSection((1:5, 2:4, 1:3), 5:8);
        @test nelements(ps, dims =1) == 4
        @test nelements(ps, dims =2) == 2
        @test nelements(ps, dims =3) == 1

        ps = ProductSection((1:5, 2:4, 1:3), 5:11);
        @test nelements(ps, dims = 1) == 5
        @test nelements(ps, dims = 2) == 3
        @test nelements(ps, dims = 3) == 1

        ps = ProductSection((1:5, 2:4, 1:3), 4:8);
        @test nelements(ps, dims = 1) == 5
        @test nelements(ps, dims = 2) == 2
        @test nelements(ps, dims = 3) == 1

        ps = ProductSection((1:5, 2:4, 1:3), 4:9);
        @test nelements(ps, dims = 1) == 5
        @test nelements(ps, dims = 2) == 2
        @test nelements(ps, dims = 3) == 1
    end

    @testset "SplittablesBase" begin
        for iters in [(1:4, 1:3), (1:4, 1:4)]
            for ps = Any[ProductSplit(iters, 3, 2), ProductSection(iters, 3:8)]
                l, r = SplittablesBase.halve(ps)
                lc, rc = SplittablesBase.halve(collect(ps))
                @test collect(l) == lc
                @test collect(r) == rc
            end
        end
    end

    @test ParallelUtilities._checknorollover((), (), ())
end;

@testset "ReverseLexicographicTuple" begin
    @testset "isless" begin
        a = ParallelUtilities.ReverseLexicographicTuple((1, 2, 3))
        b = ParallelUtilities.ReverseLexicographicTuple((2, 2, 3))
        @test a < b
        @test a <= b
        b = ParallelUtilities.ReverseLexicographicTuple((1, 1, 3))
        @test b < a
        @test b <= a
        b = ParallelUtilities.ReverseLexicographicTuple((2, 1, 3))
        @test b < a
        @test b <= a
        b = ParallelUtilities.ReverseLexicographicTuple((2, 1, 4))
        @test a < b
        @test a <= b
    end
    @testset "equal" begin
        a = ParallelUtilities.ReverseLexicographicTuple((1, 2, 3))
        @test a == a
        @test isequal(a, a)
        @test a <= a
        b = ParallelUtilities.ReverseLexicographicTuple(a.t)
        @test a == b
        @test isequal(a, b)
        @test a <= b
    end
end;

@testset "ZipSplit" begin
    @testset "SplittablesBase" begin
        for ps in [zipsplit((1:4, 1:4), 3, 2), zipsplit((1:5, 1:5), 3, 2)]
            l, r = SplittablesBase.halve(ps)
            lc, rc = SplittablesBase.halve(collect(ps))
            @test collect(l) == lc
            @test collect(r) == rc
        end
    end
end


text:
export getBoundaryCondition

"""
function getBoundaryCondition(gd::Function, M::AbstractMesh)

evaluate function gd on boundary faces of mesh M

"""
function getBoundaryCondition(gd::Function, M::AbstractMesh)
	if M.dim==2
		xf1,xf2 = getFaceGrids(M)

		gf1  = reshape(gd(xf1),tuple(M.n+[1,0]...))
		gdl  = gf1[1,:]
		gdr  = gf1[end,:]

		gf2  = reshape(gd(xf2),tuple(M.n+[0,1]...))
		gdb  = gf2[:,1]
		gdt  = gf2[:,end]

		return (gdl,gdb,gdr,gdt)
	elseif M.dim==3
		xf1,xf2,xf3 = getFaceGrids(M)

		gf1  = reshape(gd(xf1),tuple(M.n+[1,0,0]...))
		gx1   = gf1[1,:,:]
		gx2   = gf1[end,:,:]

		gf2  = reshape(gd(xf2),tuple(M.n+[0,1,0]...))
		gy1   = gf2[:,1,:]
		gy2   = gf2[:,end,:]

		gf3  = reshape(gd(xf3),tuple(M.n+[0,0,1]...))
		gz1   = gf3[:,:,1]
		gz2   = gf3[:,:,end]

		return (gx1,gy1,gz1,gx2,gy2,gz2)
	end
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
"""
    do_ump2

performs unrestricted Moller-Plesset perturbation theory to second order
on a UHF reference.
## paramters
refWfn::Wfn -> wavefunction to which the MP2 equations will be applied.

## outputs
dmp2::Float -> UMP2 energy correction
"""
function do_ump2(refWfn::Wfn)
    dmp2 = 0.0
    unr = refWfn.unrestricted
    nocca = refWfn.nalpha
    noccb = refWfn.nbeta
    rocca = 1:refWfn.nalpha
    roccb = 1:refWfn.nbeta
    rvira = refWfn.nalpha+1:refWfn.nmo
    rvirb = refWfn.nbeta+1:refWfn.nmo
    epsa = refWfn.epsa
    epsb = refWfn.epsb
    #moeri = refWfn.ijab
    #TODO: add MP2 singles energy s.t. non-HF reference can be used.
    Cao = refWfn.Cao
    Cbo = refWfn.Cbo
    Cav = refWfn.Cav
    Cbv = refWfn.Cbv
    moeri = tei_transform(refWfn.uvsr, Cao, Cav, Cao, Cav, "oovv")
    for b in rvira
        for a = b+1:refWfn.nmo
            aa = a - nocca
            bb = b - nocca
            cache = moeri[:, aa:aa, :, bb:bb] - moeri[:, bb:bb, :, aa:aa]
            cache = squeeze(cache)
            for j in rocca
                for i = j+1:refWfn.nalpha
                    moint = cache[i, j]
                    dmp2 += (moint * moint) / (epsa[i] + epsa[j] - epsa[a] - epsa[b])
                end
            end
        end
    end
    #spin case ABAB
    moeri = tei_transform(refWfn.uvsr, Cao, Cav, Cbo, Cbv, "oOvV")
    for b in rvirb
        for a in rvira
            aa = a - nocca
            bb = b - noccb
            cache = moeri[:, aa:aa, :, bb:bb]
            cache = squeeze(cache)
            for j in roccb
                for i in rocca
                    dmp2 +=
                        cache[i, j] * cache[i, j] / (epsa[i] + epsb[j] - epsa[a] - epsb[b])
                end
            end
        end
    end
    #spin case BBBB
    moeri = tei_transform(refWfn.uvsr, Cbo, Cbv, Cbo, Cbv, "OOVV")
    cache = zeros(noccb, 1, noccb, 1)
    for b in rvirb
        for a = b+1:refWfn.nmo
            aa = a - noccb
            bb = b - noccb
            cache = moeri[:, aa:aa, :, bb:bb] .- moeri[:, bb:bb, :, aa:aa]
            cache = squeeze(cache)
            for i in roccb
                for j = i+1:refWfn.nbeta
                    dmp2 +=
                        (cache[i, j] * cache[i, j]) /
                        (epsb[i] + epsb[j] - epsb[a] - epsb[b])
                end
            end
        end
    end
    return dmp2
end


text:
include("sec_6_2.jl")
using Printf

r = 3
M = 8
τ = FractionalTimeDG.reconstruction_pts(T, r)
ppI = M
store = Store(α, r+1, M, ppI)

nrows = 6
N0 = 8 # should be a multiple of 4
Nvals = Integer[ N0 * 2^i for i = 0:nrows-1 ]
err = OffsetArray{T}(undef, 0:r, 1:nrows)
for row = 1:nrows
    N = Nvals[row]
    @printf("%5d ", N)
    t, U = FODEdG!(λ, tmax, f, u0, N, r, M, store)
    pcwise_t, pcwise_U = evaluate_pcwise_poly!(U, t, τ[0:r], store)
    pcwise_u = T[ u(t, λ, u0, f, 20, store) for t in pcwise_t ]
    weight = pcwise_t.^(r-α)
    pcwise_err = OffsetArray(weight.*(pcwise_U - pcwise_u), 0:r, 1:N)
    for j = 0:r
        err[j,row] = maximum(abs.(pcwise_err[j,:]))
    end
    if row == 1
        for j = 0:r
            @printf("&%8.1e&%5s ", err[j,row], "")
        end
    else
        N_ratio = Nvals[row] / Nvals[row-1]
        for j = 0:r
            err_ratio = err[j,row-1] / err[j,row]
            rate = log(err_ratio) / log(N_ratio)
            @printf("&%8.1e&%5.2f ", err[j,row], rate)
        end
    end
    @printf("\\\\\n")
end


text:
@indextype JointID

"""
$(TYPEDEF)

The abstract supertype of all concrete joint types.
"""
abstract type JointType{T} end

Base.eltype(::Type{<:JointType{T}}) where {T} = T
isfloating(::T) where {T<:JointType} = isfloating(T)
num_velocities(::T) where {T<:JointType} = num_velocities(T)
num_positions(::T) where {T<:JointType} = num_positions(T)
num_constraints(::Type{T}) where {T<:JointType} = 6 - num_velocities(T)

"""
$(TYPEDEF)

A joint represents a kinematic restriction of the relative twist between two
rigid bodies to a linear subspace of dimension ``k``.

A joint has a direction. The rigid body before the joint is called the
joint's predecessor, and the rigid body after the joint is its successor.

The state related to the joint is parameterized by two sets of variables, namely

* a vector ``q \\in \\mathcal{Q}``, parameterizing the relative homogeneous transform.
* a vector ``v \\in \\mathbb{R}^k``, parameterizing the relative twist.

The twist of the successor with respect to the predecessor is a linear function
of ``v``.

For some joint types (notably those using a redundant representation of relative
orientation, such as a unit quaternion), ``\\dot{q}``, the time derivative of
``q``, may not be the same as ``v``.
However, an invertible linear transformation exists between
``\\dot{q}`` and ``v``.

See also:
* Definition 2.9 in Duindam, "Port-Based Modeling and Control for Efficient Bipedal Walking Robots", 2006.
* Section 4.4 of Featherstone, "Rigid Body Dynamics Algorithms", 2008.
"""
struct Joint{T, JT<:JointType{T}}
    name::String
    frame_before::CartesianFrame3D
    frame_after::CartesianFrame3D
    joint_type::JT
    id::Base.RefValue{JointID}
    joint_to_predecessor::Base.RefValue{Transform3D{T}}
    joint_to_successor::Base.RefValue{Transform3D{T}}
    position_bounds::Vector{Bounds{T}}
    velocity_bounds::Vector{Bounds{T}}
    effort_bounds::Vector{Bounds{T}}

    function Joint(name::String, frame_before::CartesianFrame3D, frame_after::CartesianFrame3D, joint_type::JointType{T};
            position_bounds::Vector{Bounds{T}}=fill(Bounds{T}(), num_positions(joint_type)),
            velocity_bounds::Vector{Bounds{T}}=fill(Bounds{T}(), num_velocities(joint_type)),
            effort_bounds::Vector{Bounds{T}}=fill(Bounds{T}(), num_velocities(joint_type))) where {T}
        JT = typeof(joint_type)
        id = Ref(JointID(-1))
        joint_to_predecessor = Ref(one(Transform3D{T}, frame_before))
        joint_to_successor = Ref(one(Transform3D{T}, frame_after))
        new{T, JT}(name, frame_before, frame_after, joint_type, id,
            joint_to_predecessor, joint_to_successor,
            position_bounds, velocity_bounds, effort_bounds)
    end
end

function Joint(name::String, jtype::JointType; kw...)
    Joint(name, CartesianFrame3D(string("before_", name)), CartesianFrame3D(string("after_", name)), jtype; kw...)
end

Base.print(io::IO, joint::Joint) = print(io, joint.name)
frame_before(joint::Joint) = joint.frame_before
frame_after(joint::Joint) = joint.frame_after
joint_type(joint::Joint) = joint.joint_type
joint_to_predecessor(joint::Joint) = joint.joint_to_predecessor[]
joint_to_successor(joint::Joint) = joint.joint_to_successor[]

"""
$(SIGNATURES)

Return a `Vector{Bounds{T}}` giving the upper and lower bounds of the
configuration for `joint`
"""
position_bounds(joint::Joint) = joint.position_bounds

"""
$(SIGNATURES)

Return a `Vector{Bounds{T}}` giving the upper and lower bounds of the
velocity for `joint`
"""
velocity_bounds(joint::Joint) = joint.velocity_bounds

"""
$(SIGNATURES)

Return a `Vector{Bounds{T}}` giving the upper and lower bounds of the
effort for `joint`
"""
effort_bounds(joint::Joint) = joint.effort_bounds

JointID(joint::Joint) = joint.id[]
Base.convert(::Type{JointID}, joint::Joint) = JointID(joint)
@inline RigidBodyDynamics.Graphs.edge_id_type(::Type{<:Joint}) = JointID
@inline RigidBodyDynamics.Graphs.edge_id(joint::Joint) = convert(JointID, joint)
@inline RigidBodyDynamics.Graphs.set_edge_id!(joint::Joint, id::JointID) = (joint.id[] = id)
function RigidBodyDynamics.Graphs.flip_direction(joint::Joint)
    jtype = RigidBodyDynamics.flip_direction(joint_type(joint))
    Joint(string(joint), frame_after(joint), frame_before(joint), jtype;
        position_bounds = .-joint.position_bounds,
        velocity_bounds = .-joint.velocity_bounds,
        effort_bounds = .-joint.effort_bounds)
end

function set_joint_to_predecessor!(joint::Joint, tf::Transform3D)
    @framecheck tf.from frame_before(joint)
    joint.joint_to_predecessor[] = tf
    joint
end

function set_joint_to_successor!(joint::Joint, tf::Transform3D)
    @framecheck tf.from frame_after(joint)
    joint.joint_to_successor[] = tf
    joint
end

function Base.show(io::IO, joint::Joint)
    if get(io, :compact, false)
        print(io, joint)
    else
        print(io, "Joint \"$(string(joint))\": $(joint.joint_type)")
    end
end

num_positions(itr) = reduce((val, joint) -> val + num_positions(joint), 0, itr)
num_velocities(itr) = reduce((val, joint) -> val + num_velocities(joint), 0, itr)

@inline function check_num_positions(joint::Joint, vec::AbstractVector)
    length(vec) == num_positions(joint) || error("wrong size")
    nothing
end

@inline function check_num_velocities(joint::Joint, vec::AbstractVector)
    length(vec) == num_velocities(joint) || error("wrong size")
    nothing
end

@propagate_inbounds function set_configuration!(q::AbstractVector, joint::Joint, config::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    copyto!(q, config)
    q
end

@propagate_inbounds function set_velocity!(v::AbstractVector, joint::Joint, vel::AbstractVector)
    @boundscheck check_num_velocities(joint, v)
    copyto!(v, vel)
    v
end

"""
$(SIGNATURES)

Return the length of the configuration vector of `joint`.
"""
@inline num_positions(joint::Joint) = num_positions(joint.joint_type)

"""
$(SIGNATURES)

Return the length of the velocity vector of `joint`.
"""
@inline num_velocities(joint::Joint) = num_velocities(joint.joint_type)

"""
$(SIGNATURES)

Return the number of constraints imposed on the relative twist between the joint's predecessor and successor
"""
@inline num_constraints(joint::Joint) = 6 - num_velocities(joint)

"""
$(SIGNATURES)

Return a `Transform3D` representing the homogeneous transform from the frame
after the joint to the frame before the joint for joint configuration vector ``q``.
"""
@propagate_inbounds function joint_transform(joint::Joint, q::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @inbounds return joint_transform(joint.joint_type, frame_after(joint), frame_before(joint), q)
end

"""
$(SIGNATURES)

Return a basis for the motion subspace of the joint in configuration ``q``.

The motion subspace basis is a ``6 \\times  k`` matrix, where ``k`` is the dimension of
the velocity vector ``v``, that maps ``v`` to the twist of the joint's successor
with respect to its predecessor. The returned motion subspace is expressed in
the frame after the joint, which is attached to the joint's successor.
"""
@propagate_inbounds function motion_subspace(joint::Joint, q::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @inbounds return motion_subspace(joint.joint_type, frame_after(joint), joint_to_predecessor(joint).to, q)
end

"""
$(SIGNATURES)

Return a basis for the constraint wrench subspace of the joint, where
`joint_transform` is the transform from the frame after the joint to the frame
before the joint.

The constraint wrench subspace is a ``6 \\times (6 - k)`` matrix, where ``k``
is the dimension of the velocity vector ``v``, that maps a vector of Lagrange
multipliers ``\\lambda`` to the constraint wrench exerted across the joint onto
its successor.

The constraint wrench subspace is orthogonal to the motion subspace.
"""
@propagate_inbounds function constraint_wrench_subspace(joint::Joint, joint_transform::Transform3D)
    @framecheck joint_transform.from frame_after(joint)
    @framecheck joint_transform.to frame_before(joint)
    @inbounds return constraint_wrench_subspace(joint.joint_type, joint_transform)
end

"""
$(SIGNATURES)

Return the acceleration of the joint's successor with respect to its predecessor
in configuration ``q`` and at velocity ``v``, when the joint acceleration
``\\dot{v}`` is zero.
"""
@propagate_inbounds function bias_acceleration(joint::Joint, q::AbstractVector, v::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_velocities(joint, v)
    @inbounds return bias_acceleration(joint.joint_type, frame_after(joint), joint_to_predecessor(joint).to, q, v)
end

"""
$(SIGNATURES)

Whether the joint's motion subspace and constraint wrench subspace depend on
``q``.
"""
has_fixed_subspaces(joint::Joint) = has_fixed_subspaces(joint.joint_type)

"""
$(SIGNATURES)

Compute joint velocity vector ``v`` given the joint configuration vector ``q``
and its time derivative ``\\dot{q}`` (in place).

Note that this mapping is linear.

See also [`velocity_to_configuration_derivative!`](@ref), the inverse mapping.
"""
@propagate_inbounds function configuration_derivative_to_velocity!(v::AbstractVector, joint::Joint, q::AbstractVector, q̇::AbstractVector)
    @boundscheck check_num_velocities(joint, v)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_positions(joint, q̇)
    @inbounds return configuration_derivative_to_velocity!(v, joint.joint_type, q, q̇)
end

"""
$(SIGNATURES)

Given  a linear function

```math
f(v) = \\langle f_v, v \\rangle
```
where ``v`` is the joint velocity vector, return a vector ``f_q`` such that

```math
\\langle f_v, v \\rangle = \\langle f_q, \\dot{q}(v) \\rangle.
```

Note: since ``v`` is a linear function of ``\\dot{q}`` (see [`configuration_derivative_to_velocity!`](@ref)),
we can write ``v = J_{\\dot{q} \\rightarrow v} \\dot{q}``, so
```math
\\langle f_v, v \\rangle = \\langle f_v, J_{\\dot{q} \\rightarrow v} \\dot{q} \\rangle = \\langle J_{\\dot{q} \\rightarrow v}^{*} f_v, \\dot{q} \\rangle
```
so ``f_q = J_{\\dot{q} \\rightarrow v}^{*} f_v``.

To compute ``J_{\\dot{q} \\rightarrow v}`` see [`configuration_derivative_to_velocity_jacobian`](@ref).
"""
@propagate_inbounds function configuration_derivative_to_velocity_adjoint!(fq, joint::Joint, q::AbstractVector, fv)
    @boundscheck check_num_positions(joint, fq)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_velocities(joint, fv)
    @inbounds return configuration_derivative_to_velocity_adjoint!(fq, joint_type(joint), q, fv)
end

"""
$(SIGNATURES)

Compute the time derivative ``\\dot{q}`` of the joint configuration vector ``q``
given ``q`` and the joint velocity vector ``v`` (in place).

Note that this mapping is linear.

See also [`configuration_derivative_to_velocity!`](@ref), the inverse mapping.
"""
@propagate_inbounds function velocity_to_configuration_derivative!(q̇::AbstractVector, joint::Joint, q::AbstractVector, v::AbstractVector)
    @boundscheck check_num_positions(joint, q̇)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_velocities(joint, v)
    @inbounds return velocity_to_configuration_derivative!(q̇, joint.joint_type, q, v)
end

"""
$(SIGNATURES)

Compute the jacobian ``J_{v \\rightarrow \\dot{q}}`` which maps joint velocity to configuration
derivative for the given joint:

```math
\\dot{q} = J_{v \\rightarrow \\dot{q}} v
```
"""
@propagate_inbounds function velocity_to_configuration_derivative_jacobian(joint::Joint, q::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @inbounds return velocity_to_configuration_derivative_jacobian(joint.joint_type, q)
end

"""
$(SIGNATURES)

Compute the jacobian ``J_{\\dot{q} \\rightarrow v}`` which maps joint
configuration derivative to velocity for the given joint:

```math
v = J_{\\dot{q} \\rightarrow v} \\dot{q}
```
"""
@propagate_inbounds function configuration_derivative_to_velocity_jacobian(joint::Joint, q::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @inbounds return configuration_derivative_to_velocity_jacobian(joint.joint_type, q)
end

"""
$(SIGNATURES)

Set ``q`` to the 'zero' configuration, corresponding to an identity joint
transform.
"""
@propagate_inbounds function zero_configuration!(q::AbstractVector, joint::Joint)
    @boundscheck check_num_positions(joint, q)
    @inbounds return zero_configuration!(q, joint.joint_type)
end

"""
$(SIGNATURES)

Set ``q`` to a random configuration. The distribution used depends on the
joint type.
"""
@propagate_inbounds function rand_configuration!(q::AbstractVector, joint::Joint)
    @boundscheck check_num_positions(joint, q)
    @inbounds return rand_configuration!(q, joint.joint_type)
end

"""
$(SIGNATURES)

Return the twist of `joint`'s  successor with respect to its predecessor,
expressed in the frame after the joint.

Note that this is the same as `Twist(motion_subspace(joint, q), v)`.
"""
@propagate_inbounds function joint_twist(joint::Joint, q::AbstractVector, v::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_velocities(joint, v)
    @inbounds return joint_twist(joint.joint_type, frame_after(joint), joint_to_predecessor(joint).to, q, v)
end

"""
$(SIGNATURES)

Return the spatial acceleration of `joint`'s  successor with respect to its predecessor,
expressed in the frame after the joint.
"""
@propagate_inbounds function joint_spatial_acceleration(joint::Joint, q::AbstractVector, v::AbstractVector, vd::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_velocities(joint, v)
    @boundscheck check_num_velocities(joint, vd)
    @inbounds return joint_spatial_acceleration(joint.joint_type, frame_after(joint), joint_to_predecessor(joint).to, q, v, vd)
end

"""
$(SIGNATURES)

Given the wrench exerted across the joint on the joint's successor, compute the
vector of joint torques ``\\tau`` (in place), in configuration `q`.
"""
@propagate_inbounds function joint_torque!(τ::AbstractVector, joint::Joint, q::AbstractVector, joint_wrench::Wrench)
    @boundscheck check_num_velocities(joint, τ)
    @boundscheck check_num_positions(joint, q)
    @framecheck(joint_wrench.frame, frame_after(joint))
    @inbounds return joint_torque!(τ, joint.joint_type, q, joint_wrench)
end

"""
$(SIGNATURES)

Compute a vector of local coordinates ``\\phi`` around configuration ``q_0``
corresponding to configuration ``q`` (in place). Also compute the time
derivative ``\\dot{\\phi}`` of ``\\phi`` given the joint velocity vector ``v``.

The local coordinate vector ``\\phi`` must be zero if and only if ``q = q_0``.

For revolute or prismatic joint types, the local coordinates can just be
``\\phi = q - q_0``, but for joint types with configuration vectors that are
restricted to a manifold (e.g. when unit quaternions are used to represent
orientation), elementwise subtraction may not make sense. For such joints,
exponential coordinates could be used as the local coordinate vector ``\\phi``.

See also [`global_coordinates!`](@ref).
"""
@propagate_inbounds function local_coordinates!(ϕ::AbstractVector, ϕ̇::AbstractVector,
        joint::Joint, q0::AbstractVector, q::AbstractVector, v::AbstractVector)
    @boundscheck check_num_velocities(joint, ϕ)
    @boundscheck check_num_velocities(joint, ϕ̇)
    @boundscheck check_num_positions(joint, q0)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_velocities(joint, v)
    @inbounds return local_coordinates!(ϕ, ϕ̇, joint.joint_type, q0, q, v)
end

"""
$(SIGNATURES)

Compute the global parameterization of the joint's configuration, ``q``, given
a 'base' orientation ``q_0`` and a vector of local coordinates ``ϕ`` centered
around ``q_0``.

See also [`local_coordinates!`](@ref).
"""
@propagate_inbounds function global_coordinates!(q::AbstractVector, joint::Joint, q0::AbstractVector, ϕ::AbstractVector)
    @boundscheck check_num_positions(joint, q)
    @boundscheck check_num_positions(joint, q0)
    @boundscheck check_num_velocities(joint, ϕ)
    @inbounds return global_coordinates!(q, joint.joint_type, q0, ϕ)
end

"""
$(SIGNATURES)

Whether the joint is a floating joint, i.e., whether it imposes no constraints
on the relative motions of its successor and predecessor bodies.
"""
isfloating(joint::Joint) = isfloating(joint.joint_type)

"""
$(SIGNATURES)

Renormalize the configuration vector ``q`` associated with `joint` so that it
lies on the joint's configuration manifold.
"""
@propagate_inbounds function normalize_configuration!(q::AbstractVector, joint::Joint)
    @boundscheck check_num_positions(joint, q)
    @inbounds return normalize_configuration!(q, joint.joint_type)
end

@propagate_inbounds function is_configuration_normalized(joint::Joint, q::AbstractVector{X}; rtol::Real = sqrt(eps(X)), atol::Real = zero(X)) where {X}
    @boundscheck check_num_positions(joint, q)
    @inbounds return is_configuration_normalized(joint.joint_type, q, rtol, atol)
end

"""
$(SIGNATURES)

Applies the principal_value functions from [Rotations.jl](https://github.com/FugroRoames/Rotations.jl/blob/d080990517f89b56c37962ad53a7fd24bd94b9f7/src/principal_value.jl)
to joint angles. This currently only applies to `SPQuatFloating` joints.
"""
@propagate_inbounds function principal_value!(q::AbstractVector, joint::Joint)
    @boundscheck check_num_positions(joint, q)
    @inbounds return principal_value!(q, joint.joint_type)
end


text:
# ---
# title: 233. Number of Digit One
# id: problem233
# author: zhwang
# date: 2022-03-12
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/number-of-digit-one/description/>
# hidden: true
# ---
#
# Given an integer n, count the total number of digit 1 appearing in all non-
# negative integers less than or equal to n.
#
# **Example:**
#
#
#
#     Input: 13
#     Output: 6
#     Explanation: Digit 1 occurred in the following numbers: 1, 10, 11, 12, 13.
#
#
#
## @lc code=start
using LeetCode

function count_digit_one(n::Int)::Int
    nums = reverse!([parse(Int, i) for i in string(n + 1)])
    res, ones = 0, count(==(1), nums)
    for (i, num) in enumerate(nums)
        num == 1 && (ones -= 1)
        res += 10^(i - 1) * (num * ones + (num > 1))
        res += i > 1 && (10^(i - 2) * (i - 1) * num)
    end
    return res
end

## example: reduce 21413 to 0
## num: step => reduced num
## 2141|4|-1: 4*2*1 + 4*0 + (4>1)*10 => 21410-1
## 214|1|0-1: 1*1*10 + 1*1*1 + (1>1)*10 => 21400-1
## 21|4|00-1: 4*1*100 + 4*2*10 + (4>1)*100 => 21000-1
## 2|1|000-1: 1*0*1000 + 1*3*100 + (1>1)*1000 => 20000-1
## |2|0000-1: 2*0*10000 + 2*4*10000 + (2>1)*10000 => 0

## @lc code=end


text:
# udat.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files udat.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
"""
module UDAT
const NONE     = Int32(-1)
const FULL     = Int32(0)
const LONG     = Int32(1)
const MEDIUM   = Int32(2)
const SHORT    = Int32(3)
const RELATIVE = Int32(128)
end # module UDAT

export UDAT, UDateFormat

macro libdat(s)     ; _libicu(s, iculibi18n, "udat_")     ; end

type UDateFormat
    ptr::Ptr{Void}

    UDateFormat(tstyle::Integer, dstyle::Integer, tz::UniStr) =
        UDateFormat(tstyle, dstyle, Vector{UInt16}(tz))
    function UDateFormat(tstyle::Integer, dstyle::Integer, tz::Vector{UInt16})
        err = UErrorCode[0]
        p = ccall(@libdat(open), Ptr{Void},
                  (Int32, Int32, Ptr{UInt8}, Ptr{UChar}, Int32,
                   Ptr{UChar}, Int32, Ptr{UErrorCode}),
                  tstyle, dstyle, locale[], tz, length(tz)-1, C_NULL, 0, err)
        FAILURE(err[1]) && error("bad date format")
        self = new(p)
        finalizer(self, close)
        self
    end

    UDateFormat(pattern::UniStr, tz::UniStr) =
        UDateFormat(Vector{UInt16}(pattern), Vector{UInt16}(tz))
    function UDateFormat(pattern::Vector{UInt16}, tz::Vector{UInt16})
        err = UErrorCode[0]
        p = ccall(@libdat(open), Ptr{Void},
                  (Int32, Int32, Ptr{UInt8}, Ptr{UChar}, Int32,
                   Ptr{UChar}, Int32, Ptr{UErrorCode}),
                  -2, -2, locale[], tz, length(tz)-1,
                  pattern, length(pattern)-1, err)
        FAILURE(err[1]) && error("bad date format")
        self = new(p)
        finalizer(self, close)
        self
    end
end

UDateFormat(pattern::AbstractString, tz::AbstractString) =
    UDateFormat(cvt_utf16(pattern), cvt_utf16(tz))

UDateFormat(tstyle::Integer, dstyle::Integer, tz::AbstractString) =
    UDateFormat(tstyle, dstyle, cvt_utf16(tz))

close(df::UDateFormat) =
    df.ptr == C_NULL || (ccall(@libdat(close), Void, (Ptr{Void},), df.ptr) ; df.ptr = C_NULL)

function format(df::UDateFormat, millis::Float64)
    err = UErrorCode[0]
    buflen = 64
    buf = zeros(UChar, buflen)
    len = ccall(@libdat(format), Int32,
                (Ptr{Void}, Float64, Ptr{UChar}, Int32, Ptr{Void}, Ptr{UErrorCode}),
                df.ptr, millis, buf, buflen, C_NULL, err)
    FAILURE(err[1]) && error("failed to format time")
    UniStr(buf[1:len+1])
end

parse(df::UDateFormat, s::AbstractString) = parse(df, cvt_utf16(s))
parse(df::UDateFormat, s::UniStr) = parse(df, Vector{UInt16}(s))
function parse(df::UDateFormat, s16::Vector{UInt16})
    err = UErrorCode[0]
    ret = ccall(@libdat(parse), Float64,
                (Ptr{Void}, Ptr{UChar}, Int32, Ptr{Int32}, Ptr{UErrorCode}),
                df.ptr, s16, length(s16)-1, C_NULL, err)
    FAILURE(err[1]) && error("failed to parse string")
    ret
end


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

# FileSystem Utils
const UNSAFE_CHARS = Dict(
    '/' => '／',
    '\\' => '＼',
    '?' => '？',
    '<' => '＜',
    '>' => '＞',
    ':' => '：',
    '*' => '＊',
    '|' => '｜',
    '"' => '”'
)

function safename(str, opts)
    mode = get(opts, :unsafe_chars, :delete)
    m = Symbol(mode)
    if m === :delete
        filter(x -> x ∉ keys(UNSAFE_CHARS), str)
    elseif m === :replace
        io = IOBuffer(str)
        out = IOBuffer()
        while !eof(io)
            c = read(io, Char)
            if haskey(UNSAFE_CHARS, c)
                write(out, UNSAFE_CHARS[c])
            else
                write(out, c)
            end
        end
        String(take!(out))
    else
        filter(x -> x ∉ keys(UNSAFE_CHARS), str)
    end
end

# URL Utils
# return dict of queries
queries(url) = url |> HTTP.URI |> HTTP.queryparams

# HTTP Utils
function print_global_useragent()
    url = "https://httpbin.org/user-agent"
    response = HTTP.get(url).body |> String
    println(response)
    return nothing
end

# Print Utils
function available_codecs(url; color=true)
    HTTP.setuseragent!(nothing)
    v_id = queries(url)["v"]
    metadata = get_metadata(v_id)
    if isnothing(metadata)
        return @error("""

            Failed to extract video info.
            Use `youtube-dl -F $url` instead!
        """)
    end
    codecs = metadata.codecs |> keys
    df_a = DataFrame(ID=Int[], Codec=String[], Bitrate=Int[])
    df_v = DataFrame(ID=Int[], Resolution=Int[], FPS=String[], Codec=String[])
    for c ∈ codecs
        if isaudio(c)
            push!(df_a, (c, name(c), bitrate(c)))
        else
            push!(df_v, (c, resolution(c), fps(c), name(c)))
        end
    end
    # sorting table
    sort!(df_a, [order(:Codec, rev=true), order(:Bitrate, rev=true)])
    sort!(df_v, [order(:Resolution, rev=true), order(:FPS, rev=true), order(:Codec, rev=true)])
    header_a = [
        "ID" "Codec" "Bitrate";
        "" "" "[Kbps]"
    ]
    header_v = [
        "ID" "Resolution" "Framerate" "Codec";
        "" "[pixels]" "[fps]" ""
    ]
    fmt = []
    if color
        push!(fmt, :header_crayon => crayon"green bold")
        #push!(fmt, :border_crayon => crayon"yellow")
        #push!(fmt, :highlighters => (hl_lt(2), hl_gt(4)))
    end
    pretty_table(df_a, header_a; title="  Audio Codecs", fmt...)
    pretty_table(df_v, header_v; title="  Video Codecs", fmt...)
end


text:
struct AffiliatedKeywordsWrapper{E <: Element} <: Element
    element::E
    keywords::Vector{<:AffiliatedKeyword}
end

AffiliatedKeywordsWrapper(element::Element, keywords::Vector{<:Pair{String, <:Any}}) =
    AffiliatedKeywordsWrapper(element, AffiliatedKeyword.(keywords))


text:
##### Beginning of file

module CommandLine # Begin submodule Snapshots.CommandLine

__precompile__(true)

import ArgParse
import Dates
import HTTP
import Pkg
import TimeZones

import ..Types
import ..Utils
import ..Common
import ..Run

import ..delayederror
import ..process_delayed_error_list

function run_snapshots_command_line!!(
        ;
        arguments::Vector{String} = String[],
        git_user_name,
        git_user_email,
        src_provider,
        dst_provider,
        include_branches,
        exclude_branches,
        do_not_push_to_these_destinations::Vector{String} = String[],
        do_not_try_url_list::Vector{String} = String[],
        try_but_allow_failures_url_list::Vector{String} =
            String[],
        time_zone::TimeZones.TimeZone =
            TimeZones.TimeZone("America/New_York"),
        )::Nothing
    @info(
        string(
            "Running Snapshots.CommandLine.",
            "run_snapshots_command_line!!"
            )
        )
    @info("parsing command line arguments...")
    parsed_arguments::Dict = _parse_arguments(
        arguments
        )
    @info("processing parsed command line arguments...")
    processed_arguments::Dict = _process_parsed_arguments(
        parsed_arguments
        )
    task::String = processed_arguments[:task]
    is_dry_run::Bool = processed_arguments[:is_dry_run]
    Run.run_snapshots!!(
        ;
        include_branches = include_branches,
        exclude_branches = exclude_branches,
        src_provider = src_provider,
        dst_provider = dst_provider,
        git_user_name = git_user_name,
        git_user_email = git_user_email,
        task = task,
        is_dry_run = is_dry_run,
        time_zone = time_zone,
        do_not_push_to_these_destinations =
            do_not_push_to_these_destinations,
        do_not_try_url_list =
            do_not_try_url_list,
        try_but_allow_failures_url_list =
            try_but_allow_failures_url_list,
        )
    process_delayed_error_list()
    return nothing
end

function _parse_arguments(arguments::Vector{String})::Dict
    s = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table s begin
        "--task"
            help = "which task to run"
            arg_type = String
            default = ""
        "--dry-run"
            help = "do everything except actually pushing the repos"
            action = :store_true
    end
    result::Dict = ArgParse.parse_args(arguments, s)
    return result
end

function _process_parsed_arguments(parsed_arguments::Dict)::Dict{Symbol, Any}
    task_argument::String = strip(
        convert(String, parsed_arguments["task"])
        )
    if length(task_argument) > 0
        task = task_argument
    else
        task = "all"
    end

    is_dry_run::Bool = parsed_arguments["dry-run"]

    result::Dict{Symbol, Any} = Dict{Symbol, Any}()
    result[:task] = task
    result[:is_dry_run]=is_dry_run
    return result
end

end # End submodule Snapshots.CommandLine

##### End of file


text:
#Test code
#-------------------------------------------------------------------------------

using SpiceData
include("importCppSimData.jl")

#No real test code yet... just demonstrate use:
stdout_ct = IOContext(stdout, :compact=>true)

testpath(testfile::String) = joinpath(CppSimData.rootpath, "core/data", testfile)

testfile = "test.tr0"
filepath = testpath(testfile)
println("\nLoad $filepath:")
reader = SpiceData._open(filepath)
@show(reader)
println("\nCompact output:")
show(stdout_ct, reader)
println("\n")

println("\nRead in list of signal names:")
@show names(reader)

signame = reader.sweepname
println("\nRead in sweep vector \"$signame\":")
t = reader.sweep
@show t[1], t[end]

signame = "vin"
println("\nRead in \"$signame\" vector:")
v = read(reader, signame)
@show v[1], v[end]

:Test_Complete


text:
using DataDeps
using ExpectationStubs

@testset "Checksum" begin
    @testset "fail" begin
        @stub dummyhash
        @expect(dummyhash(::Any) = [0x12, 0x34])

        register(
            DataDep(
                "TestFailChecksumPreupload",
                "dummy message",
                "http://example.void",
                (dummyhash, "0000")
            )
        )

        local_filepath = @__FILE__ # Need some real file to test on
        @test preupload_check("TestFailChecksumPreupload",  local_filepath) == false
        @test @used dummyhash(@__FILE__)
    end

    @testset "pass" begin
        @stub dummyhash
        @expect(dummyhash(::Any) = [0x12, 0x34])

        register(
            DataDep(
                "TestPassChecksumPreupload",
                "dummy message",
                "http://example.void",
                (dummyhash, "1234")
            )
        )

        local_filepath = @__FILE__ # Need some real file to test on
        @test preupload_check("TestPassChecksumPreupload",  local_filepath) == true
        @test @used dummyhash(@__FILE__)
    end
end

@testset "PostFetch" begin
    @testset "fail" begin
        register(
            DataDep(
                "TestFailPostFetchPreupload",
                "dummy message",
                "http://example.void",
                Any;
                post_fetch_method = error
            )
        )

        local_filepath = @__FILE__ # Need some real file to test on
        @test preupload_check("TestFailPostFetchPreupload",  local_filepath) == false
    end

    @testset "pass" begin
        @stub dummy_postfetch
        @expect(dummy_postfetch(::Any) = [0x12, 0x34])

        register(
            DataDep(
                "TestPassPostFetchPreupload",
                "dummy message",
                "http://example.void",
                Any,
                post_fetch_method = dummy_postfetch,
            )
        )

        local_filepath = @__FILE__ # Need some real file to test on
        @test preupload_check("TestPassPostFetchPreupload",  local_filepath) == true
        @test !(@used dummy_postfetch(@__FILE__)) # should not run on file given, but on a copy.
        @test @used dummy_postfetch(::Any) # should run on some file
    end
end


text:
using GLMakie, Random
Random.seed!(13)
x = -6:0.5:6
y = -6:0.5:6
z = 6exp.( -(x.^2 .+ y' .^ 2)./4)

box = Rect3(Point3f(-0.5), Vec3f(1))
colors = to_colormap(:linear_worb_100_25_c53_n256)
n = length(colors)
g(x) = x^(1/10)
alphas = [g(x) for x in LinRange(0,1,n)]
cmap_alpha = RGBAf0.(colors, alphas)

set_theme!(theme_dark())
fig, ax, = meshscatter(x, y, z; marker=box, markersize = 0.5,
    color = vec(z), colormap = cmap_alpha, colorrange = (0,6),
    axis = (; type = Axis3, aspect = :data, azimuth = 7.3,
    elevation = 0.189, perspectiveness = 0.5),
    figure = (;resolution =(1200,800)))
meshscatter!(ax, x .+ 7, y, z./2; markersize = 0.25, color = vec(z./2),
    colormap = cmap_alpha, colorrange = (0, 6),
    ambient = Vec3f0(0.85, 0.85, 0.85), backlight = 1.5f0)
xlims!(-5.5,10)
ylims!(-5.5,5.5)
hidedecorations!(ax; grid = false)
hidespines!(ax)
save("meshscatter.png", fig)
display(fig)


text:
using JuMP, Gurobi, Cbc
import XLSX
using BenchmarkTools

#### read data ####
#xf = XLSX.readxlsx("../data/part1/data_numP10_numC10_09_05_13_29.xlsx")
xf = XLSX.readxlsx("../data/data1.xlsx")

l = xf["General!B1"]     # work day duration
N = xf["General!B3"] + 1  # number of clients + 1
M = xf["General!B2"]  # number of providers
# + 1 is for including the headquarter (HS14)

# Clients
s = xf["Clients!A2:A$N"]   # client service start times
d = vcat([0], xf["Clients!B2:B$N"])   # durations for clients, set d[1] = 0 for HS14

# Providers
w = xf["Providers!B2:B$(M+1)"]  # start time of providers
f = xf["Providers!A2:A$(M+1)"]  # hiring cost of providers

# make separate s vector for each provider, inserting into the first element the provider's work start time
s_test = ones(M, N)
for k in 1:M
    global s_test[k,:] = vcat([w[k]], s)
end
s = s_test

# Distances
t = xf["Distances"][:]  # 2-d array of travel times
##########

model = Model();
set_optimizer(model, Gurobi.Optimizer);
set_time_limit_sec(model, 300.0);

#set_optimizer(model, Gurobi.Optimizer);
# set 5 mintues time limit

#(13)
@variable(model, x[1:N, 1:N, 1:M], Bin)

# subtour constraint variables using u
@variable(model, 1 <= u[1:N, 1:M] <= N-1)

######################## routing constraints ######################

# (2) each client served once by one provider
@constraint(model, served[j in 2:N], sum(x[i,j,k] for i in 1:N for k in 1:M) == 1)

# (3) route continuity
@constraint(model, continuity[p in 2:N, k in 1:M], sum(x[i,p,k] for i in 1:N) == sum(x[p,j,k] for j in 1:N) )

# (4) not all providers need to leave 14HS (and start working)
@constraint(model, leave_depot[k in 1:M], sum(x[1,j,k] for j in 2:N) <= 1)

# (7) at least 1 provider must leave 14HS and start working to serve clients
@constraint(model, sum(x[1,j,k] for j in 2:N for k in 1:M) >= 1)

#(9) subtour elimination for each provider
for k in 1:M
    for i in 2:N
        for j in 2:N
            @constraint(model, u[i,k] - u[j,k] + (N-1) * x[i,j,k] <= (N-1) - 1 )
        end
    end
end

######################## end routing constraints ######################

######################## time constraints ######################

# (10) initial client service times
@constraint(model, start[j in 2:N, k in 1:M], x[1,j,k]*(s[k,1] + t[1,j] + d[1]) <= s[k,j])

# (11) in between start and end of work, service times
@constraint(model, bet[i in 2:N, j in 2:N, k in 1:M], x[i,j,k]*(s[k,i] + t[i,j] + d[i]) <= s[k,j])

# (12) come back to 14HS by hour l.
@constraint(model, endofwork[i in 2:N, k in 1:M], x[i,1,k]*(s[k,i] + t[i,1] + d[i]) <= l)

######################## end time constraints ######################

# newer objective to minimise total cost based on hours of work.
### Note. the max() function inside the objective can only be solved by Gurobi solver (Cbc cannot solve)
@objective(model, Min, sum(x[1,j,k]*f[k]*(sum(x[p,q,k]*(t[p,q] + d[p] + max(0, s[k,q] - s[k,p] - t[p,q] - d[p]) ) for p in 1:N for q in 1:N)) for j in 2:N for k in 1:M))

# original non-hourly objective for part 1
#@objective(model, Min, sum(f[k]*x[1,j,k] for j in 1:N for k in 1:M))

optimize!(model)

#@btime optimize!(model)

# see which providers are chosen to provide service.
for k in 1:M
    if sum(value.(x[1,:,k])) == 1
        println("Provider $k provides service.")
        println("Provider $k costs $(f[k]) per hour.")
        # come back hr
        for i in 2:N
            if value.(x[i,1,k]) == 1
                println("Provider $k comes back to 14HS at $(t[i,1] + s[k,i] + d[i])hr. \n")
                break # added
            end
        end
    end
end

# Given N: number of clients+1, and provider number k, find the path
### So if you found that providers 2, 3, 4 and 5 are chosen to work, run get_path(N, 2), get_path(N, 3) etc.
function get_path(N, k)
    path = []
    for j in 2:N
        if value.(x[1,j,k]) == 1
            println("\nstart of path")
            println(1, " => ")
            path = vcat(path, [1])
            println(j, " => ")
            path = vcat(path, [j])
            i = j
            while value.(x[i,1,k]) == 0
                for m in 2:N
                    if value.(x[i,m,k]) == 1
                        println(m, " => ")
                        path = vcat(path, [m])
                        i = m
                        break
                    end
                end
            end
            println(1)
            path = vcat(path, [1])
            println("end of path.... or in client index below (not location)")
        end
    end

    for j in 2:N
        if value.(x[1,j,k]) == 1
            println("\nstart of path")
            println(0, " => ")
            println(j-1, " => ")
            i = j
            while value.(x[i,1,k]) == 0
                for m in 2:N
                    if value.(x[i,m,k]) == 1
                        println(m-1, " => ")
                        i = m
                        break
                    end
                end
            end
            println(0)
            println("end of path (in client index)")
        end
    end
    return path
end

k = 2
path = get_path(N, k)

function overview(path, k)
    println("$(s[k,1])hr: Becomes available, leaves to location $(path[2]).")
    for i in 2:length(path)-1
        println("$(s[k,path[i-1]]+t[path[i-1], path[i]] + d[path[i-1]])hr: Arrived at location $(path[i]).")
        if s[k,path[i-1]]+t[path[i-1], path[i]] + d[path[i-1]] != s[k,path[i]]
            println("   Arrived early... Wait until service start time.")
        end
        println("$(s[k,path[i]])hr: Serve for $(d[path[i]]) hrs.")
        println("$(s[k,path[i]]+d[path[i]])hr: Finished service. Leave to location $(path[i+1]).")
    end
    println("$(s[k,path[length(path)-1]] + d[path[length(path)-1]] + t[path[length(path)-1], 1])hr: Arrived at location 1 (14HS).\n")
end

# print the overview of the path for provider k, where 'path' is found above using get_path(N,k)
path = [1 7 1]
k = 2
overview(path, k)

# find the number of worked hours = provider come back time (to 14HS) -  provider start time (w_k)
worked_hr = s[k,path[length(path)-1]] + d[path[length(path)-1]] + t[path[length(path)-1], 1] - s[k,1]
println("Provider $k worked for $worked_hr hours.")


text:
using Documenter, ViscousFlow

ENV["GKSwstype"] = "nul" # removes GKS warnings during plotting

makedocs(
    sitename = "ViscousFlow.jl",
    doctest = true,
    clean = true,
    pages = [
        "Home" => "index.md",
        "Manual" => ["manual/1.-Basic-viscous-flow.md",
                     "manual/2.-Basic-flow-with-a-stationary-body.md",
                     "manual/3.-Applying-pulse-forcing-to-a-flow.md",
                     "manual/4.-Multiple-stationary-bodies.md",
                     "manual/5.-Viscous-flow-about-a-moving-body.md",
                     "manual/6.-Variable-free-stream.md",
                     "manual/7.-Lid-driven-cavity.md",
                     "manual/functions.md"
                     ]
        #"Internals" => [ "internals/properties.md"]
    ],
    #format = Documenter.HTML(assets = ["assets/custom.css"])
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict()
            )
        ))
    ),
    #assets = ["assets/custom.css"],
    #strict = true
)

#if "DOCUMENTER_KEY" in keys(ENV)
deploydocs(
     repo = "github.com/JuliaIBPM/ViscousFlow.jl.git",
     target = "build",
     deps = nothing,
     make = nothing
     #versions = "v^"
)
#end


text:
module SymbolicTensors2
using LinearAlgebra
using SymbolicUtils
using SymbolicUtils: istree, operation, arguments, similarterm

abstract type AbstractNVector end
abstract type AbstractTensor <: AbstractNVector end
struct NVector <: AbstractNVector
    nvec::Array
    idxs::Tuple{Int}
    metric::Array

end

abstract type AbstractContraction end
#struct Contraction <: AbstractContraction
    #A::NVector
    #B::NVector
#end

struct Contraction <: AbstractContraction
    ex :: Expr
end

SymbolicUtils.istree(x::Contraction) = x.ex.head == :call
SymbolicUtils.operation(x::Contraction) = x.ex.args[1]
SymbolicUtils.arguments(x::Contraction) = x.ex.args[2:end]
SymbolicUtils.symtype(x::Contraction) = Real # ?
SymbolicUtils.similarterm(x::Contraction, f, xs) = Contraction(:(I * I))

export Contraction
end


text:
N = 3
E = SNN.IF(!;N = N)
EE = SNN.SpikingSynapse(E, E, :ge; σ=0.5, p=0.8)
for n in 1:(N - 1)
    SNN.connect!(EE, n, n + 1, 50)
end
E.I[1] = 30

SNN.monitor(E, [(:v, [1, N])])
SNN.sim!([E], [EE]; duration = 100ms)
#SNN.vecplot(E, :v) |> display


text:
__precompile__()

"""
    import Eglob: eglob, eglobt

`Eglob` provides extended globbing functions that support a recursive `**`
wildcard pattern, similar to the `**` supported by the
[Bash](https://www.gnu.org/software/bash/) (version 4+) `globstar` option
and by [Zsh](http://www.zsh.org/). This package is built on top of
[Glob.lj](https://github.com/vtjnash/Glob.jl).

* `eglob()`: An extended globbing function that supports leading `~` expansion
  and a recursive `**` wildcard pattern
* `eglobt()`: A Task- and Channel-based version of `eglob()`.

See the docs for the individual functions for more information.
"""
module Eglob

export eglob, eglobt

using Glob: glob

pathsep = Base.Filesystem.path_separator

#     _pathsplit(path)
#
# Split a path into an array of path components, using the file separator
# ('/' on POSIX systems, '\' on Windows) that's appropriate for the
# underlying operating system. Does not take drive letters into account.
# If there's a Windows drive letter in the path, it'll end up with the
# first component.
#
# **Parameters**
#
# * `path`: path to split. Can be relative or absolute
#
# **Returns** an array of path components
function _pathsplit(path::AbstractString)
    split(path, Base.Filesystem.path_separator_re, keep=false)
end

"""
    eglob(pattern)

An extended globbing function that supports all the wildcards of the `glob()`
function, plus:

- a leading `~`, signifying the user's home directory (`~user` isn't supported)
- a special `**` wildcard that recursively matches any directory.

`eglob()` is just a convenience front-end. It is exactly equivalent to:

    [match for match in eglobt(pattern)]

**Parameters:**

* `pattern`: The wildcard pattern

**Returns** an array of expanded path names, which might be empty. If the
original path contains no wildcards, the result will be a 1-element array
containing the original path.

**Note:** This function expands the entire pattern before returning. If you're
expanding a pattern in a deeply-nested directory tree, `eglob()` can take
awhile to return. If you want to process the results in a way that allows
you to stop before seeing all of them, use `eglobt()`.
"""
function eglob(pattern::AbstractString)
    [match for match in eglobt(pattern)]
end

"""
    eglobt(pattern)

An extended, Task- and Channel-based globbing function that supports all the
wildcards of the `glob()` function, plus:

- a leading `~`, signifying the user's home directory (`~user` isn't supported)
- a special `**` wildcard that recursively matches any directory.

**Parameters:**

* `pattern`: The wildcard pattern

**Returns** the channel to which results will be written.

**Examples**

Example 1:

```julia-repl
julia> import Glob: eglobt
julia> c = eglobt("**/*.csv") # find all CSV files below the current directory
julia> match1 = take!(c)
julia> match2 = take!(c)
```

Example 2:

```julia
# We only want three of them, but there might not be that many, so handle
# a prematurely empty channel.
import Glob: eglobt
c = eglobt("src/**/*.scala")
matches = []
for (i, match) in enumerate(c)
    if i > 3
        break
    end
    append!(matches, [match])
end
```

Example 3:

```julia
# Consume all matches. This is exactly what eglob() does.
matches = [match for match in eglobt("test/**/*.py")]
```

"""
function eglobt(pattern::AbstractString)

    dotslash = ".$(pathsep)"

    function putmatch!(channel::Channel, match::String)
        # Matching artifact: The code adds "./" to the beginning of matches
        # for things like "**/*.c". Remove it, if necessary.
        m = if startswith(match, dotslash) && !startswith(pattern, dotslash)
            match[3:end]
        else
            match
        end
        put!(channel, m)
    end

    function handle_recursive_glob(channel, directory, pieces)
        last = length(pieces) == 1
        remaining_pieces = last ? [] : @views pieces[2:end]
        for (root, dirs, files) in walkdir(directory, topdown=true)
            if last
                # At the end of a pattern, "**" just recursively matches
                # directories.
                putmatch!(channel, root)
            else
                # Recurse downward, trying to match the rest of the
                # pattern.
                if length(remaining_pieces) > 0
                    do_glob(channel, remaining_pieces, root)
                end
            end
        end
    end

    function handle_non_recursive_glob(channel, directory, piece, pieces)
        last = length(pieces) == 1
        path = joinpath(directory, piece)
        matches = glob(path)
        if length(matches) > 0
            if last
                # Save the last set of matches, and stop.
                for m in matches
                    putmatch!(channel, m)
                end
            else
                # Have to continue recursing.
                remaining_pieces = @views pieces[2:end]
                for m in filter(f -> isdir(f), matches)
                    do_glob(channel, remaining_pieces, m)
                end
            end
        end
    end

    function do_glob(channel, pieces, directory)
        if length(pieces) == 0
            return
        end

        piece = pieces[1]

        if piece == "**"
            handle_recursive_glob(channel, directory, pieces)
        else
            # It's just a regular (non-recursive) glob pattern.
            handle_non_recursive_glob(channel, directory, piece, pieces)
        end
    end

    function split_windows_eglob_pattern(pattern::AbstractString)
        drive, path = splitdrive(pattern)
        if (drive, path) == ("", "")
            (".", ".")
        elseif drive == ""
            (path, ".")
        elseif path == ""
            (".", drive)
        else
            # Hack: Can't handle non-absolute paths in a drive. Pretend a
            # drive letter means "absolute". Note that "drive" can be empty
            # here, which is fine.
            path2 = path[1] == '\\' ? path[2:end] : path
            (path2, drive * '\\')
        end
    end

    function split_posix_eglob_pattern(pattern::AbstractString)
        if length(pattern) == 0
            (".", ".")
        elseif pattern[1] == '/'
            (pattern[2:end], "/")
        else
            (pattern, ".")
        end
    end

    function main(channel::Channel)
        eglob_pattern_splitter = (is_windows() ? split_windows_eglob_pattern
                                               : split_posix_eglob_pattern)

        # Main eglob() logic.

        # First, account for leading "~"
        adjusted_pattern =
            if length(pattern) == 0
                "."
            elseif pattern[1] == '~'
                normpath(expanduser(pattern))
            else
                pattern
            end

        # Determine leading directory, which is different per OS (because of
        # Windows' stupid drive letters).
        relative_pattern, directory = eglob_pattern_splitter(adjusted_pattern)

        # Do the actual globbing.
        pieces = _pathsplit(relative_pattern)
        do_glob(channel, pieces, directory)
    end

    # Kick it off
    Channel(main)
end

end


text:
module TestARMA

using QuantEcon
using Base.Test
using FactCheck

# set up
phi = [.95, -.4, -.4]
theta = zeros(3)
sigma = .15
lp = ARMA(phi, theta, sigma)

facts("Testing arma.jl") do
    # test simulate
    sim = simulation(lp, ts_length=250)
    @fact length(sim) --> 250

    # test impulse response
    imp_resp = impulse_response(lp, impulse_length=75)
    @fact length(imp_resp) --> 75

    context("test constructors") do
        phi = 0.5
        theta = 0.4
        sigma = 0.15

        a1 = ARMA(phi, theta, sigma)
        a2 = ARMA([phi;], theta, sigma)
        a3 = ARMA(phi, [theta;], sigma)

        for nm in fieldnames(a1)
            @fact getfield(a1, nm) --> getfield(a2, nm)
            @fact getfield(a1, nm) --> getfield(a3, nm)
        end
    end

    context("test autocovariance") do
        θ = 0.5
        σ = 0.15
        ma1 = ARMA(Float64[], [θ], σ)
        ac = autocovariance(ma1; num_autocov=5)

        # first is the variance. equal to (1 + θ^2) sigma^2
        @fact ac[1] --> roughly((1+θ^2)*σ^2; atol=1e-3)

        # second should be θ σ^2
        @fact ac[2] --> roughly(θ*σ^2; atol=1e-3)

        # all others should be 0
        @fact ac[3:end] --> roughly(zeros(ac[3:end]); atol=1e-3)
    end

end  # facts
end  # module


text:
include("../imports.jl")

@testset ExtendedTestSet "ImageKeypointRegression" begin
    method = ImageKeypointRegression((16, 16), 10)
    DLPipelines.checkmethod_core(method)
    FastAI.checkmethod_plot(method)
end


text:
module words_extractor_jl

using Distributed
using YAML

const folder = "words"

function worker(yaml_path)
    path = get_filepath(yaml_path)
    words = get_words(yaml_path)
    write(path, join(words, "\n"))
    println(string("Saved...", path))
end

function get_words(yaml_path)
    text_path = replace(yaml_path, ".yml" => ".txt")
    text = read(text_path, String) |> lowercase
    split(text, r"[\W\d]+") |> Set |> collect
end

function get_filepath(path)
    meta = YAML.load_file(path)
    string(folder, "/extracted-words-for-", meta["label"], ".txt")
end

function walk(path, file_ext)
    res = []
    for (root, _, files) in walkdir(path, topdown = true)
        for file in files
            if endswith(file, file_ext)
                filepath = joinpath(root, file)
                push!(res, filepath)
            end
        end
    end
    res
end

function main()
    if ispath(folder)
        rm(folder, recursive = true)
    end
    mkdir(folder)
    Threads.@threads for path in walk("../data/pl/", ".yml")
        println("Spawn $path")
        worker(path)
    end
end

# addprocs()
# println(string("Workers ", nworkers()))
println(string("Processing... using ", Threads.nthreads(), " threads"))
@time main()
end # module


text:
# Modified from Reexport.jl, with one difference: it doesn't reexport the module names themselves
#
# Copyright (c) 2014: Simon Kornblith.
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
# associated documentation files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial
# portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT
# LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
macro reexport(ex)
    isa(ex, Expr) && (ex.head == :module ||
        ex.head == :using ||
        ex.head == :importall ||
        (ex.head == :toplevel &&
        all(e->isa(e, Expr) && (e.head == :using || e.head == :importall), ex.args))) ||
    error("@reexport: syntax error")

    if ex.head == :module
    modules = Any[ex.args[2]]
    ex = Expr(:toplevel, ex, :(using .$(ex.args[2])))
    elseif (ex.head == :using || ex.head == :importall) && all(e->isa(e, Symbol), ex.args)
    modules = Any[ex.args[end]]
    else
    modules = Any[e.args[end] for e in ex.args]
    end

    Expr(:toplevel,
        ex,
        [:(eval(Expr(:export, setdiff(names($mod), [Compat.nameof($mod)])...))) for mod in modules]...
    ) |> esc
end


text:
using LinAlgExtensions
using Test
using LinearAlgebra
using SparseArrays
using SuiteSparse

@testset "utils"     begin include("utils.jl") end
@testset "randutils" begin include("randutils.jl") end
@testset "pinvfact"  begin include("pinv.jl") end


text:
module SthenoCompat

# Write your package code here.

using Stheno
using Zygote
import Stheno: AV, ew, pw

include("transformed_linear_kernel.jl")

end


text:
mutable struct DSPFilter
    a::Vector{Float64}
    b::Vector{Float64}
    si::Vector{Float64}
    need_restart::Bool

    function DSPFilter()
        new(Float64[], Float64[], Float64[], true)
    end
end

# MovingMean

# MovingDiff

# MovingSpectra

#
function init(obj::DSPFilter, ftype::FilterType, coefs::FilterCoefficients)
    d = digitalfilter(ftype, coefs) |> DF2TFilter
    obj.a = coefa(d.coef)
    obj.b = coefb(d.coef)

    order = length(obj.a) - 1
    resize!(obj.si, order)
    fill!(obj.si, 0)

    obj
end

function run_filter(obj::DSPFilter, y::Vector, x::Vector)
    a, b, si = obj.a, obj.b, obj.si
    order = length(si)

    for i in 1:length(x)
        xi = x[i]
        yi = xi * b[1] + si[1]
        for j in 2:order
            si[j-1] = si[j] + b[j] * xi - a[j] * yi
        end
        si[order] = b[order+1] * xi - a[order+1] * yi
        y[i] = yi
    end
    return obj
end


text:
#
# Factors Join Tests
#

@testset "Factors Join" begin

@testset "main" begin
    ϕ1 = Factor(:X =>3, :C =>2, :Y =>3, :A=>2, :Z=>3)
    ϕ2 = Factor(:A=>2, :B=>3, :C=>2)
    ϕ1.potential[:] = collect(1.0:length(ϕ1))
    ϕ2.potential[:] = collect(1.0:length(ϕ2))

    # why did I do this?
    true_pot = [1.0,2.0,3.0,28.0,35.0,42.0,7.0,8.0,9.0,70.0,77.0,84.0,13.0,14.0,
            15.0,112.0,119.0,126.0,38.0,40.0,42.0,176.0,184.0,192.0,50.0,52.0,
            54.0,224.0,232.0,240.0,62.0,64.0,66.0,272.0,280.0,288.0,37.0,38.0,
            39.0,280.0,287.0,294.0,43.0,44.0,45.0,322.0,329.0,336.0,49.0,50.0,
            51.0,364.0,371.0,378.0,110.0,112.0,114.0,464.0,472.0,480.0,122.0,
            124.0,126.0,512.0,520.0,528.0,134.0,136.0,138.0,560.0,568.0,576.0,
            73.0,74.0,75.0,532.0,539.0,546.0,79.0,80.0,81.0,574.0,581.0,588.0,
            85.0,86.0,87.0,616.0,623.0,630.0,182.0,184.0,186.0,752.0,760.0,768.0,
            194.0,196.0,198.0,800.0,808.0,816.0,206.0,208.0,210.0,848.0,856.0,
            864.0,3.0,6.0,9.0,36.0,45.0,54.0,21.0,24.0,27.0,90.0,99.0,108.0,39.0,
            42.0,45.0,144.0,153.0,162.0,76.0,80.0,84.0,220.0,230.0,240.0,100.0,
            104.0,108.0,280.0,290.0,300.0,124.0,128.0,132.0,340.0,350.0,360.0,
            111.0,114.0,117.0,360.0,369.0,378.0,129.0,132.0,135.0,414.0,423.0,
            432.0,147.0,150.0,153.0,468.0,477.0,486.0,220.0,224.0,228.0,580.0,
            590.0,600.0,244.0,248.0,252.0,640.0,650.0,660.0,268.0,272.0,276.0,
            700.0,710.0,720.0,219.0,222.0,225.0,684.0,693.0,702.0,237.0,240.0,
            243.0,738.0,747.0,756.0,255.0,258.0,261.0,792.0,801.0,810.0,364.0,
            368.0,372.0,940.0,950.0,960.0,388.0,392.0,396.0,1000.0,1010.0,1020.0,
            412.0,416.0,420.0,1060.0,1070.0,1080.0,5.0,10.0,15.0,44.0,55.0,66.0,
            35.0,40.0,45.0,110.0,121.0,132.0,65.0,70.0,75.0,176.0,187.0,198.0,
            114.0,120.0,126.0,264.0,276.0,288.0,150.0,156.0,162.0,336.0,348.0,
            360.0,186.0,192.0,198.0,408.0,420.0,432.0,185.0,190.0,195.0,440.0,
            451.0,462.0,215.0,220.0,225.0,506.0,517.0,528.0,245.0,250.0,255.0,
            572.0,583.0,594.0,330.0,336.0,342.0,696.0,708.0,720.0,366.0,372.0,
            378.0,768.0,780.0,792.0,402.0,408.0,414.0,840.0,852.0,864.0,365.0,
            370.0,375.0,836.0,847.0,858.0,395.0,400.0,405.0,902.0,913.0,924.0,
            425.0,430.0,435.0,968.0,979.0,990.0,546.0,552.0,558.0,1128.0,1140.0,
            1152.0,582.0,588.0,594.0,1200.0,1212.0,1224.0,618.0,624.0,630.0,
            1272.0,1284.0,1296.0]

    ϕ12 = ϕ1 * ϕ2
    @test ϕ12.potential[:] == true_pot
end

@testset "2" begin
    # dimensions don't have the same lengths
    ϕ1 = Factor(:X => 3, :Y =>3, :Z=>3)
    ϕ2 = Factor(:A=>3, :X=>31, :Y=>2)

    @test_throws DimensionMismatch ϕ1 * ϕ2
end

@testset "swapping" begin
    # current implementation swaps the smaller one to the end
    ϕ1 = Factor(:Y, [2, 3, 4])
    ϕ2 = Factor([:X, :Y, :Z], reshape(collect(1:27), (3, 3, 3)))
    ϕ3 = ϕ1 - ϕ2

    p1 = ϕ1.potential
    p2 = ϕ2.potential
    p3 = reshape(p1, 1, 3, 1) .- p2

    @test_approx_eq p3 ϕ3.potential
end

@testset "singleton" begin
    ϕ2 = Factor(:A=>2, :B=>3, :C=>2)
    ϕ2.potential[:] = collect(1:length(ϕ2))
    ϕ1 = Factor(0.1)
    ϕ3 = ϕ1 * ϕ2

    @test_approx_eq ϕ3.potential ϕ2.potential * 0.1
end
@testset "singleton 2" begin
    ϕ1 = Factor(:A=>2, :B=>3, :C=>2)
    ϕ1.potential[:] = collect(1:length(ϕ1))
    ϕ2 = Factor(0.1)
    ϕ3 = ϕ1 * ϕ2

    @test_approx_eq ϕ3.potential ϕ1.potential * 0.1
end

@testset "both singleton" begin
    ϕ1 = Factor(3)
    ϕ2 = Factor(-2)
    ϕ3 = join(atan2, ϕ2, ϕ1)
    @test_approx_eq first(ϕ3.potential) atan2(-2, 3)
end

@testset "none common" begin
    ϕ1 = Factor(:X=>2, :Y=>3)
    ϕ2 = Factor(:A=>3, :B=>2)
    rand!(ϕ1)
    rand!(ϕ2)
    ϕ3 = ϕ1 * ϕ2

    p1 = ϕ1.potential
    p2 = ϕ2.potential
    p3 = Factors.duplicate(p1, 3, 2)
    p3 .*= reshape(p2, (1, 1, 3, 2))

    @test_approx_eq ϕ3.potential p3
end

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of cstr reactor
#*--------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type cstr_basic
	cstr_basic()=begin
		NComp=outers.NComp
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number Of Components"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number Of Reactions"
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Stoichiometric Matrix"
			)),(NComp, NReac)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Stream",
				:PosX=>0,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Stream",
				:PosX=>1,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Heat"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Reacting Volume"
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar total amount"
			)),
			fill(conc_mol(Dict{Symbol,Any}(
				:Brief=>"Components concentration"
			)),(NComp)),
			[
				:(Outlet.z * M = C * Vr),
				:(sum(Outlet.z) = 1),
			],
			[
				"Molar Concentration","Molar fraction",
			],
			[:NComp,:NReac,:stoic,],
			[:Inlet,:Outlet,:q,:Vr,:M,:C,]
		)
	end
	NComp::DanaInteger
	NReac::DanaInteger
	stoic::Array{DanaReal}
	Inlet::stream
	Outlet::streamPH
	q::heat_rate
	Vr::volume
	M::mol
	C::Array{conc_mol}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export cstr_basic
function setEquationFlow(in::cstr_basic)
	addEquation(1)
	addEquation(2)
end
function atributes(in::cstr_basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Basic model for a CSTR reactor"
	drive!(fields,_)
	return fields
end
cstr_basic(_::Dict{Symbol,Any})=begin
	newModel=cstr_basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
__precompile__(true)

module GraphLayout
    using Compose   # for plotting features

    typealias AdjList{T} Vector{Vector{T}}

    # Spring-based force layout algorithm
    export layout_spring_adj
    include("spring.jl")

    # Stress majorization layout algorithm
    export layout_stressmajorize_adj
    include("stress.jl")

    # Tree layout algorithms
    export layout_tree
    include("tree.jl")
    # Heuristic algortihms for tree layout
    include("tree_heur.jl")
    # Optimal algorithms for tree layout, that require JuMP
    include("tree_opt.jl")

    # Drawing utilities
    export draw_layout_adj
    include("draw.jl")
end