text:
module TestDeprecated

using Test, DataFrames, CategoricalArrays

const ≅ = isequal

@testset "DataFrame!" begin
    x = [1, 2, 3]
    y = [4, 5, 6]
    @test DataFrame!(x=x, y=y, copycols=true) == DataFrame(x=x, y=y)
    df1 = DataFrame(x=x, y=y)
    df2 = DataFrame!(df1)
    @test df1 == df2
    @test df1.x === df2.x
    @test df1.y === df2.y

    a=[1, 2, 3]
    df = DataFrame!(:a=>a, :b=>1, :c=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!("a"=>a, "b"=>1, "c"=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!(Dict(:a=>a, :b=>1, :c=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!(Dict("a"=>a, "b"=>1, "c"=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!((x, y))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), (:x1, :x2))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), ("x1", "x2"))
    @test names(df) == ["x1", "x2"]
    @test df."x1" === x
    @test df."x2" === y

    @test_throws MethodError DataFrame!([Union{Int, Missing}, Union{Float64, Missing}],
                                        [:x1, :x2], 2)
end

@testset "test categorical" begin
    df = DataFrame(x=["a", "b", "c"],
                   y=["a", "b", missing],
                   z=[1, 2, 3])
    for x in [df, view(df, :, :)]
        y = categorical(x)
        @test y isa DataFrame
        @test x ≅ y
        @test x.x !== y.x
        @test x.y !== y.y
        @test x.z !== y.z
        @test y.x isa CategoricalVector{String}
        @test y.y isa CategoricalVector{Union{Missing, String}}
        @test y.z isa Vector{Int}

        y = categorical(x, Int)
        @test y isa DataFrame
        @test x ≅ y
        @test x.x !== y.x
        @test x.y !== y.y
        @test x.z !== y.z
        @test y.x isa Vector{String}
        @test y.y isa Vector{Union{Missing, String}}
        @test y.z isa CategoricalVector{Int}

        for colsel in [:, names(x), [1, 2, 3], [true, true, true], r"", Not(r"a")]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa CategoricalVector{String}
            @test y.y isa CategoricalVector{Union{Missing, String}}
            @test y.z isa CategoricalVector{Int}
        end

        for colsel in [:x, "x", 1, [:x], ["x"], [1], [true, false, false], r"x", Not(2:3)]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa CategoricalVector{String}
            @test y.y isa Vector{Union{Missing, String}}
            @test y.z isa Vector{Int}
        end

        for colsel in [:z, "z", 3, [:z], ["z"], [3], [false, false, true], r"z", Not(1:2)]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa Vector{String}
            @test y.y isa Vector{Union{Missing, String}}
            @test y.z isa CategoricalVector{Int}
        end

        for colsel in [Int[], Symbol[], [false, false, false], r"a", Not(:)]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa Vector{String}
            @test y.y isa Vector{Union{Missing, String}}
            @test y.z isa Vector{Int}
        end
    end
end

@testset "test categorical!" begin
    using DataFrames: _columns
    df = DataFrame(A = Vector{Union{Int, Missing}}(1:3), B = Vector{Union{Int, Missing}}(4:6))
    DRT = CategoricalArrays.DefaultRefType
    @test all(c -> isa(c, Vector{Union{Int, Missing}}), eachcol(categorical!(deepcopy(df))))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), [1, 2])))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), [:A, :B])))
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), [:A]))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), :A))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), [1]))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), 1))) == 1

    @test all(c -> isa(c, Vector{Union{Int, Missing}}), eachcol(categorical!(deepcopy(df))))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), Not(Not([1, 2])))))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), Not(Not([:A, :B])))))
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not([:A]))))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not(:A))))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not([1]))))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not(1))))) == 1
end

@testset "categorical!" begin
    df = DataFrame([["a", "b"], ['a', 'b'], [true, false], 1:2, ["x", "y"]], :auto)
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df)))),
                  [CategoricalArrays.CategoricalValue{String, UInt32},
                   Char, Bool, Int,
                   CategoricalArrays.CategoricalValue{String, UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), :))),
                  [CategoricalArrays.CategoricalValue{String, UInt32},
                   CategoricalArrays.CategoricalValue{Char, UInt32},
                   CategoricalArrays.CategoricalValue{Bool, UInt32},
                   CategoricalArrays.CategoricalValue{Int, UInt32},
                   CategoricalArrays.CategoricalValue{String, UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), compress=true))),
                  [CategoricalArrays.CategoricalValue{String, UInt8},
                   Char, Bool, Int,
                   CategoricalArrays.CategoricalValue{String, UInt8}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), names(df)))),
                  [CategoricalArrays.CategoricalValue{String, UInt32},
                   CategoricalArrays.CategoricalValue{Char, UInt32},
                   CategoricalArrays.CategoricalValue{Bool, UInt32},
                   CategoricalArrays.CategoricalValue{Int, UInt32},
                   CategoricalArrays.CategoricalValue{String, UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), names(df), compress=true))),
                  [CategoricalArrays.CategoricalValue{String, UInt8},
                   CategoricalArrays.CategoricalValue{Char, UInt8},
                   CategoricalArrays.CategoricalValue{Bool, UInt8},
                   CategoricalArrays.CategoricalValue{Int, UInt8},
                   CategoricalArrays.CategoricalValue{String, UInt8}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), Not(1:0)))),
                  [CategoricalArrays.CategoricalValue{String, UInt32},
                   CategoricalArrays.CategoricalValue{Char, UInt32},
                   CategoricalArrays.CategoricalValue{Bool, UInt32},
                   CategoricalArrays.CategoricalValue{Int, UInt32},
                   CategoricalArrays.CategoricalValue{String, UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), Not(1:0), compress=true))),
                  [CategoricalArrays.CategoricalValue{String, UInt8},
                   CategoricalArrays.CategoricalValue{Char, UInt8},
                   CategoricalArrays.CategoricalValue{Bool, UInt8},
                   CategoricalArrays.CategoricalValue{Int, UInt8},
                   CategoricalArrays.CategoricalValue{String, UInt8}]))

    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), Integer))),
                  [String, Char,
                   CategoricalArrays.CategoricalValue{Bool, UInt32},
                   CategoricalArrays.CategoricalValue{Int, UInt32},
                   String]))

    df = DataFrame([["a", missing]], :auto)
    categorical!(df)
    @test df.x1 isa CategoricalVector{Union{Missing, String}}

    df = DataFrame(x1=[1, 2])
    categorical!(df)
    @test df.x1 isa Vector{Int}
    categorical!(df, :)
    @test df.x1 isa CategoricalVector{Int}
end

@testset "categorical with Cols, All and Between" begin
    df = DataFrame(x1=["a", "b"], y=[2, 3])
    categorical(df, All())
    categorical(df, Cols())
    categorical(df, Between(1, 2))
    categorical!(df, All())
    categorical!(df, Cols())
    categorical!(df, Between(1, 2))
end

@testset "deprecated describe syntax" begin
    @test describe(DataFrame(a=[1, 2]), cols = :a, :min, :min2 => minimum, "max2" => maximum, :max) ==
          DataFrame(variable=:a, min=1, min2=1, max2=2, max=2)
end

@testset "All indexing" begin
    df = DataFrame(a=1, b=2, c=3)

    @test select(df, All(1, 2)) == df[:, 1:2]
    @test select(df, All(1, :b)) == df[:, 1:2]
    @test select(df, All(:a, 2)) == df[:, 1:2]
    @test select(df, All(:a, :b)) == df[:, 1:2]
    @test select(df, All(2, 1)) == df[:, [2, 1]]
    @test select(df, All(:b, 1)) == df[:, [2, 1]]
    @test select(df, All(2, :a)) == df[:, [2, 1]]
    @test select(df, All(:b, :a)) == df[:, [2, 1]]

    @test df[:, All(1, 2)] == df[:, 1:2]
    @test df[:, All(1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2)] == df[:, 1:2]
    @test df[:, All(:a, :b)] == df[:, 1:2]
    @test df[:, All(2, 1)] == df[:, [2, 1]]
    @test df[:, All(:b, 1)] == df[:, [2, 1]]
    @test df[:, All(2, :a)] == df[:, [2, 1]]
    @test df[:, All(:b, :a)] == df[:, [2, 1]]

    @test df[:, All(1, 1, 2)] == df[:, 1:2]
    @test df[:, All(:a, 1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2, :b)] == df[:, 1:2]
    @test df[:, All(:a, :b, 2)] == df[:, 1:2]
    @test df[:, All(2, 1, :a)] == df[:, [2, 1]]

    @test select(df, All(1, "b")) == df[:, 1:2]
    @test select(df, All("a", 2)) == df[:, 1:2]
    @test select(df, All("a", "b")) == df[:, 1:2]
    @test select(df, All("b", 1)) == df[:, [2, 1]]
    @test select(df, All(2, "a")) == df[:, [2, 1]]
    @test select(df, All("b", "a")) == df[:, [2, 1]]

    @test df[:, All(1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2)] == df[:, 1:2]
    @test df[:, All("a", "b")] == df[:, 1:2]
    @test df[:, All("b", 1)] == df[:, [2, 1]]
    @test df[:, All(2, "a")] == df[:, [2, 1]]
    @test df[:, All("b", "a")] == df[:, [2, 1]]

    @test df[:, All("a", 1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2, "b")] == df[:, 1:2]
    @test df[:, All("a", "b", 2)] == df[:, 1:2]
    @test df[:, All(2, 1, "a")] == df[:, [2, 1]]

    df = DataFrame(a1=1, a2=2, b1=3, b2=4)
    @test df[:, All(r"a", Not(r"1"))] == df[:, [1, 2, 4]]
    @test df[:, All(Not(r"1"), r"a")] == df[:, [2, 4, 1]]
end

@testset "deprecated DataFrame constructors" begin
    @test DataFrame(([1, 2], [3, 4])) == DataFrame([[1, 2], [3, 4]], :auto)
    @test DataFrame((categorical([1, 2]), categorical([3, 4]))) ==
          DataFrame([categorical([1, 2]), categorical([3, 4])], :auto)
    @test DataFrame(([1, 2], [3, 4]), ("a", "b")) == DataFrame([[1, 2], [3, 4]], ["a", "b"])
    @test DataFrame(([1, 2], [3, 4]), (:a, :b)) == DataFrame([[1, 2], [3, 4]], [:a, :b])
    @test DataFrame(([1, 2, 3], [1, 2, 3])) == DataFrame((1:3, 1:3)) == DataFrame((1:3, [1, 2, 3]))
    @test DataFrame(("x1"=>1:3, "x2"=>[1, 2, 3])) == DataFrame(["x1"=>1:3, "x2"=>[1, 2, 3]])
    @test DataFrame((:x1=>1:3, :x2=>[1, 2, 3])) == DataFrame([:x1=>1:3, :x2=>[1, 2, 3]])
    @inferred DataFrame((1:3, 1:3))
    @inferred DataFrame((1:3, 1:3), (:a, :b))
    @inferred DataFrame((1:3, 1:3), ("a", "b"))
    @inferred DataFrame((:x1=>1:3, :x2=>[1, 2, 3]))
    @inferred DataFrame(("x1"=>1:3, "x2"=>[1, 2, 3]))
    @test DataFrame(Union{Float64, Missing}[0.0 1.0;
                                            0.0 1.0;
                                            0.0 1.0]) ==
          convert(DataFrame, Union{Float64, Missing}[0.0 1.0;
                                                     0.0 1.0;
                                                     0.0 1.0])
    @test names(DataFrame([0.0 1.0;
                           0.0 1.0;
                           0.0 1.0], ["a", "b"])) == ["a", "b"]
    @test names(DataFrame([0.0 1.0;
                           0.0 1.0;
                           0.0 1.0], [:a, :b])) == ["a", "b"]

    x = [1, 2, 3]
    y = [1, 2, 3]

    df = DataFrame((x, y))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), copycols=true)
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), copycols=false)
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame((x, y), (:x1, :x2))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), (:x1, :x2), copycols=true)
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), (:x1, :x2), copycols=false)
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame((x, y), ("x1", "x2"))
    @test names(df) == ["x1", "x2"]
    @test df."x1" == x
    @test df."x2" == y
    @test df."x1" !== x
    @test df."x2" !== y
    df = DataFrame((x, y), ("x1", "x2"), copycols=true)
    @test names(df) == ["x1", "x2"]
    @test df."x1" == x
    @test df."x2" == y
    @test df."x1" !== x
    @test df."x2" !== y
    df = DataFrame((x, y), ("x1", "x2"), copycols=false)
    @test names(df) == ["x1", "x2"]
    @test df."x1" === x
    @test df."x2" === y

    df = DataFrame([1 2; 3 4], :auto)
    @test size(df) == (2, 2)
    @test df.x1 == [1, 3]
    @test df.x2 == [2, 4]

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                   [:A, :B, :C], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[!, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[!, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[!, 3]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[!, 1])
    @test all(ismissing, df[!, 2])
    @test all(ismissing, df[!, 3])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                   ["A", "B", "C"], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[!, "A"]) == Vector{Union{Int, Missing}}
    @test typeof(df[!, "B"]) == Vector{Union{Float64, Missing}}
    @test typeof(df[!, "C"]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[!, "A"])
    @test all(ismissing, df[!, "B"])
    @test all(ismissing, df[!, "C"])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}], [:x1, :x2], 2)
    @test size(df) == (2, 2)
    @test eltype.(eachcol(df)) == [Union{Int, Missing}, Union{Float64, Missing}]

    @test_throws MethodError DataFrame([Union{Int, Missing}, Union{Float64, Missing}],
                                       [:x1, :x2], 2, copycols=false)
    @test size(df) == (2, 2)
    @test eltype.(eachcol(df)) == [Union{Int, Missing}, Union{Float64, Missing}]

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                   [:A, :B, :C])
    @test size(df, 1) == 0
    @test size(df, 2) == 3
    @test typeof(df[!, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[!, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[!, 3]) == Vector{Union{String, Missing}}
    @test propertynames(df) == [:A, :B, :C]

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                   ["A", "B", "C"])
    @test size(df, 1) == 0
    @test size(df, 2) == 3
    @test typeof(df[!, "A"]) == Vector{Union{Int, Missing}}
    @test typeof(df[!, "B"]) == Vector{Union{Float64, Missing}}
    @test typeof(df[!, "C"]) == Vector{Union{String, Missing}}
    @test names(df) == ["A", "B", "C"]

    df = convert(DataFrame, zeros(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[!, 1]) == Vector{Float64}
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, ones(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[!, 1]) == Vector{Float64}
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, Matrix{Float64}(undef, 10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[!, 1]) == Vector{Float64}
    @test typeof(df[:, 1]) == Vector{Float64}
end

end # module


text:
#!/home/kimura/julia/julia
#
# メカ女子将棋システム (C) 2013 メカ女子将棋部☆
#
# メカ女子将棋部：
#
# 竹部　さゆり（女流三段）様
#   メカさゆりん
# 渡辺　弥生（女流一級）様
#   メカみおたん
# T.R.　（女子大学院生）様
#   メカりえぽん
# 木村　健（メカウーサーメカ担当、実装責任者）
#   メカきむりん、プロジェクトリーダー、小五女子w
#
#

const MECHA_JYOSHI_SHOGI = 1
const MECHAJYO_VERSION = "0.6"

srand(1234)

#require("Profile")
#using IProfile

require("BoardConsts.jl") #remain when registration
require("Board.jl")       #remain when registration
require("Move.jl")        #remain when registration
require("GameStatus.jl")  #remain when registration
require("OldGenMove.jl")  #comment out when registration
require("BitBoard.jl")    #remain when registration
require("fvbin.jl")       #comment out when registration
require("Eval.jl")        #comment out when registration
require("EvalBonanza.jl") #comment out when registration
require("GenMove.jl")     #comment out when registration
require("Search.jl")      #comment out when registration
require("PVS.jl")         #comment out when registration
require("SimpleSearch.jl")#comment out when registration

# おまじない
function setupIO()
    stdinx::Ptr{Uint8} = 0
    stdoutx::Ptr{Uint8} = 0
    ret = ccall((:setunbuffering, "../lib/libMJ.so.1"), Int32, ())

    gs = InitGS()
    InitTables(gs)
    #BBTest2(gs)
    #BBTestForMoveGeneration(gs)
    return gs
end

function producer()
    while true
        s = readline(STDIN)
        s = chomp(s)
        produce(s)
        if s == "quit" || s == "exit"
            break
        end
    end
end

function main()
    sfenHirate = "lnsgkgsnl/1r5b1/ppppppppp/9/9/9/PPPPPPPPP/1B5R1/LNSGKGSNL w - 1"
    li::Array{ASCIIString,1} = ["" for x = 1:100]::Array{ASCIIString,1}
    li2::Array{ASCIIString,1} = ["" for x = 1:100]::Array{ASCIIString,1}
    lilen = 0
    lilen2 = 0
    side::Int = SENTE # default
    count::Int = 0
    count2::Int = 0
    gs = setupIO()
    #GenMoveTest(gs)
    t = Task(producer)
    while true
        st = consume(t)
        #println("consumed $st")
        if st == "quit" || st == "exit"
            break
        elseif st == "usi"
            println("id name Mecha Jyoshi Shogi ",MECHAJYO_VERSION)
            println("id author Sayuri TAKEBE, Mio WATANABE, Rieko TSUJI and Takeshi KIMURA")
            println("option name BookFile type string default $(gs.bookfile)")
            println("option name UseBook type check default $(gs.usebook)")
            println("usiok");
        elseif st == "isready"
            srand(time_ns())
            println("readyok")
        elseif beginswith(st,"setoption")
            if beginswith(st,"setoption name USI_Ponder value true")
                gs.canponder = true
            elseif beginswith(st,"setoption name USI_Ponder value false")
                gs.canponder = false
            elseif beginswith(st,"setoption name USI_Hash value ")
                gs.hashsize = uint32(st[length("setoption name USI_Hash value "):])
            elseif beginswith(st,"setoption name BookFile value ")
                gs.bookfile = st[length("setoption name BookFile value "):]
            else
            end
        elseif beginswith(st,"usinewgame")
            # do nothing
            out = [Move(0,0,0,0,0,0) for n = 1:30000]
            history = [0 for x = 1:NumSQ, y = 1:NumSQ]
        elseif beginswith(st,"position startpos")
            li = split(st)
            count = 0
            count2 = 0
            if li[2] == "startpos" # 平手
                if size(li,1) == 2 # 初手
                    gs.board = InitSFEN(sfenHirate, gs.board)
                    count = 0
                    side = SENTE
                else
                    if size(li,1) == 4
                        side = GOTE
                    end
                    gs.board = InitSFEN(sfenHirate, gs.board)
                    count = 0

                    for x = 4:size(li,1)
                        count = 0
                        if iseven(x) # 先手
                            count2 = generateMoves(gs.board, out, SENTE, count, gs)
                            idx = findIndex(out,li[x],SENTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println("Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                println("Warning: invalid moves:", li[x])
                            else
                                #println("index = ", idx)
                                makeMove(gs.board,
		                         idx,
		                         out,
		                         SENTE)
                                #println("move=",move2USIString(out[idx]))
                            end
                        else # 後手
                            count2 = generateMoves(gs.board, out, GOTE, count, gs)
                            idx = findIndex(out,li[x],GOTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println("Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                println("Warning: invalid moves:", li[x])
                            else
                                #println("index = ", idx)
                                makeMove(gs.board,
		                         idx,
		                         out,
		                         GOTE)
                                #println("move=",move2USIString(out[idx]))
                            end
                        end

                    end
                end
                #DisplayBoard(gs.board)
            end
        elseif beginswith(st,"position sfen")
            li = split(st)
            count = 0
            count2 = 0
            sfen = join(li[3:6]," ")
            #println("sfen=",sfen)
            if li[2] == "sfen" # 平手
                if size(li,1) == 6 # 初手
                    gs.board = InitSFEN(sfen, gs.board)
                    count = 0
                    side = SENTE
                else
                    if size(li,1) == 8
                        side = GOTE
                    end
                    gs.board = InitSFEN(sfen, gs.board)
                    count = 0

                    for x = 8:size(li,1)
                        count = 0
                        if iseven(x) # 先手
                            count2 = generateMoves(gs.board, out, SENTE, count, gs)
                            idx = findIndex(out,li[x],SENTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println("Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                printlpn("Warning: invalid moves:", li[x])
                            else
                                #println("index = ", idx)
                                makeMove(gs.board,
		                         idx,
		                         out,
		                         SENTE)
                                #println("move=",move2USIString(out[idx]))
                            end
                        else # 後手
                            count2 = generateMoves(gs.board, out, GOTE, count, gs)
                            idx = findIndex(out,li[x],GOTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println("Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                println("Warning: invalid moves:", li[x])
                            else
                                #println("index = ", idx)
                                makeMove(gs.board,
		                         idx,
		                         out,
		                         GOTE)
                                #println("move=",move2USIString(out[idx]))
                            end
                        end
                    end
                end
            end
            #DisplayBoard(gs.board)
        elseif beginswith(st,"go")
            li2 = split(st)
            btime::Int = parseint(li2[3],10)
            wtime::Int = parseint(li2[5],10)
            byoyomi::Int = parseint(li2[7],10)
            #if in_check( side, gs.board)
            #    println("check!")
            #end
            count2 = generateMoves(gs.board, out, side, 0, gs)
            count3 = generateBB(gs.board, out, side, 0, gs)

            if count2 == 0
                println("bestmove resign")
            else
                # chose random moves
                #randomIndex::Int = rand(Uint32) % (count2)
                Index::Int = -1
                #m::Move = think(side,gs)
                m::Move = thinkASP(side,gs)
                #println("move=",m)
                for q = 1:count2
                    if out[q].move == m.move
                        Index = q
                        break
                    end
                end
                if Index == -1
                    println("bestmove resign")
                else
                    makeMove(gs.board,Index,out,side)
                    if in_check( side, gs.board)
                        #println("check!")
                    end
                    println("bestmove ",move2USIString(out[Index]))
                    #@iprofile report
                end
            end
        elseif beginswith(st,"gameover")
            # do nothing
        else
            println("COMMAND NOT FOUND($st)")
        end
        #println("COMMAND: $st")
    end
end

main()

ext: jl

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess


text:
__precompile__()

module POSIXClock

const librt = Libdl.find_library(["librt.so"])
const juliasrc = "julia"
include(joinpath(juliasrc,"clock.jl"))

end


text:
module ArrayInterfaces

import Base: connect, close, getindex, start, next, done, length, isopen, isempty
importall ..DBAPIBase
import Compat.view

### Underlying column data structures

abstract AbstractColumn{T}

immutable Column{T} <: AbstractColumn{T}
    name::Symbol
    data::Vector{T}
end

immutable SubColumn{T} <: AbstractColumn{T}
    name::Symbol
    data::SubArray{T}
end

SubColumn(column::Column, indices...) = SubColumn(column.name, view(column.data, indices...))

start(column::AbstractColumn) = start(column.data)
next{T}(column::AbstractColumn{T}, state) = next(column.data, state)::Tuple{T, Any}
done(column::AbstractColumn, state) = done(column.data, state)

getindex(column::AbstractColumn, indexes...) = getindex(column.data, indexes...)
length(column::AbstractColumn) = length(column.data)

function allequal{T}(things::AbstractArray{T})
    thing_set = Set{T}()

    for thing in things
        push!(thing_set, thing)

        if length(thing_set) > 1
            return false
        end
    end

    return true
end

### Interface

type ColumnarArrayInterface <: DatabaseInterface end

### Connections

type ColumnarArrayConnection <: DatabaseConnection{ColumnarArrayInterface}
    columns::Vector{Column}
    closed::Bool
end

function connect(
        ::Type{ColumnarArrayInterface},
        names::AbstractArray{Symbol},
        columns::AbstractArray{Vector}
    )

    if length(names) != length(columns)
        throw(ArrayInterfaceError("Arrays of names and columns must be the same length"))
    end

    if !allequal(map(length, columns))
        throw(ArrayInterfaceError("Columns must be the same length"))
    end

    ColumnarArrayConnection(
        map(names, columns) do name, column
            Column(name, column)
        end,
        false,
    )
end

function close(connection::ColumnarArrayConnection)
    connection.closed = true
    empty!(connection.columns)

    return nothing
end

# does not support transactions
commit(connection::ColumnarArrayConnection) = nothing

isopen(connection::ColumnarArrayConnection) = !connection.closed

### Cursors

type ColumnarArrayCursor <: FixedLengthDatabaseCursor{ColumnarArrayInterface}
    connection::ColumnarArrayConnection
    columns::Vector{SubColumn}

    ColumnarArrayCursor(connection) = new(connection)
end

cursor(connection::ColumnarArrayConnection) = ColumnarArrayCursor(connection)

connection(cursor::ColumnarArrayCursor) = cursor.connection

isempty(cursor::ColumnarArrayCursor) = isempty(cursor.columns)

length(cursor::ColumnarArrayCursor) = isempty(cursor) ? 0 : length(first(cursor.columns))

### Queries

immutable ColumnarArrayQuery{T<:OrdinalRange} <: DatabaseQuery
    columns::Vector{Symbol}
    rows::T
end

function execute!(cursor::ColumnarArrayCursor, query::ColumnarArrayQuery)
    try
        # handle case where there are rows requested from no columns
        if isempty(query.columns) && !isempty(query.rows)
            throw(BoundsError())
        end

        remaining_columns = Set(query.columns)
        cursor.columns = map(filter(cursor.connection.columns) do col
                is_queried = col.name in query.columns
                if is_queried
                    pop!(remaining_columns, col.name)
                end
                is_queried
            end) do col
            SubColumn(col, query.rows)
        end

        # don't allow queries for nonexistent columns
        if !isempty(remaining_columns)
            throw(BoundsError())
        end
    catch error
        if isa(error, BoundsError) || isa(error, KeyError)
            rethrow(DatabaseQueryError(interface(cursor), query))
        else
            rethrow(error)
        end
    end

    return nothing
end

### Results

immutable ColumnarArrayRowIterator
    cursor::ColumnarArrayCursor
end

function rows(cursor::ColumnarArrayCursor)
    if !isdefined(cursor, :columns)
        throw(ArrayInterfaceError("No query has been run on $cursor"))
    end

    return ColumnarArrayRowIterator(cursor)
end

function start(iter::ColumnarArrayRowIterator)
    if isempty(iter.cursor.columns)
        # needed because the method zip() does not exist
        zip_iterator = zip(())
    else
        zip_iterator = zip(iter.cursor.columns...)
    end

    return (zip_iterator, start(zip_iterator))
end

function next(iter::ColumnarArrayRowIterator, state)
    (zip_iterator, current_state) = state
    (next_set, next_state) = next(zip_iterator, current_state)
    return (next_set, (zip_iterator, next_state))
end

function done(iter::ColumnarArrayRowIterator, state)
    return done(state...)
end

immutable ColumnarArrayColumnIterator
    cursor::ColumnarArrayCursor
end

function columns(cursor::ColumnarArrayCursor)
    if !isdefined(cursor, :columns)
        throw(ArrayInterfaceError("No query has been run on $cursor"))
    end

    return ColumnarArrayColumnIterator(cursor)
end

start(iter::ColumnarArrayColumnIterator) = start(iter.cursor.columns)

function next(iter::ColumnarArrayColumnIterator, state)
    (next_col, next_state) = next(iter.cursor.columns, state)
    return (collect(next_col), next_state)
end

done(iter::ColumnarArrayColumnIterator, state) = done(iter.cursor.columns, state)

function getindex(cursor::ColumnarArrayCursor, row_ind::Int, column_name::Symbol)
    if_failed = BoundsError(cursor, (row_ind, column_name))

    for column in cursor.columns
        if column.name == column_name
            try
                return column[row_ind]
            catch error
                rethrow(if_failed)
            end
        end
    end

    throw(if_failed)
end

function getindex(cursor::ColumnarArrayCursor, row_ind::Int, column_ind::Int)
    if_failed = BoundsError(cursor, (row_ind, column_ind))

    try
        return cursor.columns[column_ind][row_ind]
    catch error
        rethrow(if_failed)
    end
end

### Errors

immutable ArrayInterfaceError{T<:AbstractString} <: DatabaseError{ColumnarArrayInterface}
    message::T
end

if VERSION >= v"0.5-"
    Base.iteratorsize(::ColumnarArrayRowIterator) = Base.SizeUnknown()
    Base.iteratoreltype(::ColumnarArrayRowIterator) = Base.EltypeUnknown()
    Base.iteratorsize(::ColumnarArrayColumnIterator) = Base.SizeUnknown()
    Base.iteratoreltype(::ColumnarArrayColumnIterator) = Base.EltypeUnknown()
end

end  # module


text:
using PyPlot, DataFrames, Query
pygui(true)

function sosagg(csvraw, agglabel)
	csvform = tempname()
	if !(ispath(csvform))
		run(pipeline(`./sos.awk`, stdin=csvraw, stdout=csvform))
	end
	sosmeta = split(readlines(open(csvraw))[1], [','; '/'])
	measure = lowercase(strip(sosmeta[2]))
	reg = strip(sosmeta[end-1])
	if reg == "Riket"
		reg = "Sverige"
	end

	sos = readtable(csvform)
	sosa = @from i in sos begin
			@group i by @NT(Ålder=i.Ålder, Kön=i.Kön) into g
			@select {Ålder=g.key.Ålder,
				Kön=g.key.Kön,
				Incidens=mean(g..Incidens)}
			@collect DataFrame
		end
	return Dict(:sosa => sosa,
		:agglabel=>agglabel,
		:minyr=>minimum(sos[:År]),
		:maxyr=>maximum(sos[:År]),
		:measure=>measure,
		:reg => reg)
end

function sexplot(sosad)
	sosa = sosad[:sosa]
	lab = sosad[:agglabel]
	minyr = sosad[:minyr]
	maxyr = sosad[:maxyr]
	measure = sosad[:measure]
	reg = sosad[:reg]
	for sex in ["Kvinnor"; "Män"]
		sosas = sosa[sosa[:Kön].==sex, :]
		plot(sosas[:Ålder], log(sosas[:Incidens]./10^5),
			"-o", label = "$sex")
	end
	grid(1)
	xlabel("Ålder")
	ylabel("log(incidens, $measure)")
	title("Sjukhusvårdade $lab $reg $minyr\u2013$maxyr")
	legend()
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
using Test
include("test_1.jl")
include("test_2.jl")
include("test_3.jl")
include("test_4.jl")
include("test_5.jl")


text:
using Test
using MPIMapReduce
using MPI
MPI.Init()

const comm = MPI.COMM_WORLD
const np = MPI.Comm_size(comm)
const rank = MPI.Comm_rank(comm)

function reductiontest(op)
	arraysize = (10, 1)
	for n in 1:5, (fmap, itrs) in Any[
		(x -> ones(arraysize) * x, (1:n,)),
		((x, y) -> ones(arraysize) * (x + y), (1:n, 1:n)),
		]

		yexp = mapreduce(fmap, (x,y) -> op.(x,y), itrs...)

		for root in 0:MPI.Comm_size(comm) - 1
			y = MPIMapReduce.pmapreduce(fmap, op, itrs..., root = root)
			if MPI.Comm_rank(comm) == root
				@test y == yexp
			else
				@test y == nothing
			end
		end
	end

	@test_throws ArgumentError MPIMapReduce.pmapreduce(x -> ones(2), op, 1:0)
	@test_throws ArgumentError MPIMapReduce.pmapreduce(x -> ones(2), op, 1:1, root = -1)
	@test_throws ArgumentError MPIMapReduce.pmapreduce(x -> ones(2), op, 1:1, root = np + 1)
end

function reductiontest(op, MPIop)
    arraysize = (10, 1)
    for (fmap, itrs) in Any[
        (x -> ones(arraysize) * x, (1:np,)),
        ((x, y) -> ones(arraysize) * (x + y), (1:np, 1:np)),
        ]

        yexp = mapreduce(fmap, (x,y) -> op.(x,y), itrs...)

        for root in 0:MPI.Comm_size(comm) - 1
            y = MPIMapReduce.pmapreduce(fmap, MPIop, itrs..., root = root)
            if MPI.Comm_rank(comm) == root
                @test y == yexp
            else
                @test y == nothing
            end
        end
    end
end

@testset "Reduction" begin
	@testset "elementwise sum" begin
		reductiontest(+)
        reductiontest(+, MPI.SUM)
	end
	@testset "elementwise product" begin
		reductiontest(*)
        reductiontest(*, MPI.PROD)
	end
	@testset "elementwise max" begin
		reductiontest(max)
        reductiontest(max, MPI.MAX)
	end
	@testset "elementwise min" begin
		reductiontest(min)
        reductiontest(min, MPI.MIN)
	end
end

function concattest(op)
	for n in 1:10, (ind, (fmap, itrs)) in  enumerate(Any[
		(x -> x^2, (1:n,)),
		((x, y) -> (x^2 + y^2), (1:n, 1:n)),
		(x -> ones(2) * x^2, (1:n,)),
		((x, y) -> ones(2) * (x^2 + y^2), (1:n, 1:n)),
		(x -> ones(1,2) * x^2, (1:n,)),
		((x, y) -> ones(1,2) * (x^2 + y^2), (1:n, 1:n)),
		(x -> ones(2,2,2) * x^2, (1:n,)),
		((x, y) -> ones(2,2,2) * (x^2 + y^2), (1:n, 1:n)),
		])

		yexp = mapreduce(fmap, op, itrs...)

		for root in 0:MPI.Comm_size(comm) - 1
			y = MPIMapReduce.pmapgatherv(fmap, op, itrs..., root = root)
			if MPI.Comm_rank(comm) == root
				@test begin
					res = y == yexp
					if !res
						@show n, np, ind, itrs, root, y, yexp
					end
					res
				end
			else
				@test y === nothing
			end
		end
	end

	@test_throws ArgumentError MPIMapReduce.pmapgatherv(x -> ones(2), op, 1:0)
	@test_throws ArgumentError MPIMapReduce.pmapgatherv(x -> ones(2), op, 1:1, root = -1)
	@test_throws ArgumentError MPIMapReduce.pmapgatherv(x -> ones(2), op, 1:1, root = np + 1)
end

@testset "Concatenation" begin
	@testset "vcat" begin
		rank == 0 && @info "vcat"
	    concattest(vcat)
	end

	@testset "hcat" begin
		rank == 0 && @info "hcat"
		concattest(hcat)
	end

	@testset "Cat" begin
	    @testset "single dimension" begin
	    	@testset "Cat(1) is vcat" begin
	    		rank == 0 && @info "Cat(1) == vcat"
	    		fmap = x -> x^2
	    		for itrs in [(1:1,), (1:4,)]
	    			yexp = mapreduce(fmap, vcat, itrs...)
	    			for root in 1:MPI.Comm_size(comm) - 1
			    	    y1 = MPIMapReduce.pmapgatherv(fmap, MPIMapReduce.Cat(1), itrs..., root = root)
			    	    y2 = MPIMapReduce.pmapgatherv(fmap, vcat, itrs..., root = root)
			    	    if MPI.Comm_rank(comm) == root
					    	@test y1 == y2 == yexp
					    else
					    	@test y1 === y2 === nothing
					    end
			    	end
		    	end
	    	end
	    	@testset "Cat(2) is hcat" begin
	    		rank == 0 && @info "Cat(2) == hcat"
	    		fmap = x -> x^2
	    		for itrs in [(1:1,), (1:4,)]
	    			yexp = mapreduce(fmap, hcat, itrs...)
	    			for root in 1:MPI.Comm_size(comm) - 1
					    y1 = MPIMapReduce.pmapgatherv(fmap, MPIMapReduce.Cat(2), itrs..., root = root)
					    y2 = MPIMapReduce.pmapgatherv(fmap, hcat, itrs..., root = root)
					    if MPI.Comm_rank(comm) == root
					    	@test y1 == y2 == yexp
					    else
					    	@test y1 === y2 === nothing
					    end
					end
				end
	    	end
	    	rank == 0 && @info "Cat single dimension"
	        concattest(MPIMapReduce.Cat(3))
	        concattest(MPIMapReduce.Cat(Val(3)))
	        concattest(MPIMapReduce.Cat(4))
	        concattest(MPIMapReduce.Cat(Val(4)))
	    end
	    @testset "multiple dimensions" begin
	    	rank == 0 && @info "Cat multiple dimensions"
	        concattest(MPIMapReduce.Cat(3:4))
	        concattest(MPIMapReduce.Cat([3,4]))
	        concattest(MPIMapReduce.Cat([1,4]))
	    end
	end
end


text:
using Printf
using Plots

using Oceananigans
using Oceananigans.Advection
using JULES

using JULES: IdealGas
using Oceananigans.Grids: Cell, xnodes

ENV["GKSwstype"] = "100"

@inline x′(x, t, L, U) = mod(x + L/2 - U * t, L) - L/2
@inline ϕ_Square(x, t; L, U, w=0.15) = -w <= x′(x, t, L, U) <= w ? 1.0 : 0.0

function every(n)
      0 < n <= 128 && return 1
    128 < n <= 256 && return 2
    256 < n <= 512 && return 4
    512 < n        && return 8
end

function periodic_advection_verification(N, L, T, U, CFL, advection, solution)
    ideal_gas = IdealGas(Float64, JULES.R⁰, 1; T₀=1, p₀=1, s₀=1, u₀=0)

    topo = (Periodic, Periodic, Periodic)
    domain = (x=(-L/2, L/2), y=(0, 1), z=(0, 1))
    grid = RegularCartesianGrid(topology=topo, size=(N, 1, 1), halo=(4, 4, 4); domain...)

    Δt = CFL * grid.Δx / abs(U)
    Nt = ceil(Int, T/Δt)

    model = CompressibleModel(
                      grid = grid,
                 advection = advection,
                     gases = (ρ=ideal_gas,),
    thermodynamic_variable = Entropy(),
                   closure = IsotropicDiffusivity(ν=0, κ=0),
                   gravity = 0.0,
               tracernames = (:ρ, :ρs, :ρc)
    )

    gas = model.gases.ρ
    R, cₚ, cᵥ = gas.R, gas.cₚ, gas.cᵥ
    u₀₀, T₀₀, ρ₀₀, s₀₀ = gas.u₀, gas.T₀, gas.ρ₀, gas.s₀

    ρ₀(x, y, z) = 1
    p₀(x, y, z) = 1
    T₀(x, y, z) = p₀(x, y, z) / (R*ρ₀(x, y, z))
    ρs₀(x, y, z) = ρ₀(x, y, z) * (s₀₀ + cᵥ * log(T₀(x, y, z)/T₀₀) - R * log(ρ₀(x, y, z)/ρ₀₀))

    set!(model.tracers.ρ, ρ₀)
    set!(model.momenta.ρu, U)
    set!(model.tracers.ρs, ρs₀)
    update_total_density!(model)

    initial_condition(x, y, z) = solution(x, 0)
    set!(model.tracers.ρc, initial_condition)
    set!(model.momenta.ρv, initial_condition)
    set!(model.momenta.ρw, initial_condition)

    anim = @animate for n in 1:Nt
        @info "Running periodic advection [N=$N, CFL=$CFL, $(typeof(advection)), U=$U]... iteration $n/$Nt"

        time_step!(model, Δt)

        x = xnodes(Cell, grid)
        analytic_solution = solution.(x, model.clock.time)
        ρc = interior(model.tracers.ρc)[:]
        ρv = interior(model.momenta.ρv)[:]
        ρw = interior(model.momenta.ρv)[:]

        title = @sprintf("N=%d, CFL=%.2f %s", N, CFL, typeof(advection))
        plot(x, analytic_solution, lw=2, label="analytic", title=title, xlims=(-L/2, L/2), ylims=(-0.2, 1.2), dpi=200)
        plot!(x, ρc, lw=2, ls=:solid, label="Oceananigans ρc")
        plot!(x, ρv, lw=2, ls=:dash,  label="Oceananigans ρv")
        plot!(x, ρw, lw=2, ls=:dot,   label="Oceananigans ρw")

    end every every(Nt)

    filename = @sprintf("periodic_advection_N%d_CFL%.2f_%s_U%+d.mp4", N, CFL, typeof(advection), U)
    mp4(anim, filename, fps=15)
end

T = 2
L = 1

Ns = [64]
Us = [+1, -1]
CFLs = [0.1]
advection_schemes = [CenteredSecondOrder(), CenteredFourthOrder(), UpwindBiasedThirdOrder(), UpwindBiasedFifthOrder(), WENO5()]

for N in Ns, CFL in CFLs, scheme in advection_schemes, U in Us
    @info "Running periodic advection [N=$N, CFL=$CFL, $(typeof(scheme)), U=$U]..."
    solution(x, t) = ϕ_Square(x, t, L=L, U=U)
    periodic_advection_verification(N, L, T, U, CFL, scheme, solution)
end


text:
module voxel_bracket_examples
using FinEtools

function voxel_bracket_mesh()
    V = VoxelBoxVolume(Int, 6*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    h1 = solidcylinder((2.0, 2.5, 2.5), (1.0, 0.0, 0.0), 00.75)
    fillsolid!(V, differenceop(unionop(b1, b2), h1), 1)

    fens, fes = H8voximg(V.data, vec([voxeldims(V)...]), [1])
    fens = meshsmoothing(fens, fes; npass = 5)
    println("count(fes) = $(count(fes))")

    File = "voxel_bracket_mesh.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)

end # voxel_bracket_mesh

function voxel_bracket_mesh_tet()
    V = VoxelBoxVolume(Int, 6*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    h1 = solidcylinder((2.0, 2.5, 2.5), (1.0, 0.0, 0.0), 00.75)
    fillsolid!(V, differenceop(unionop(b1, b2), h1), 1)

    fens, fes = T4voximg(V.data, vec([voxeldims(V)...]), [1])
    fens = meshsmoothing(fens, fes; method = :laplace, npass = 5)
    println("count(fes) = $(count(fes))")

    File = "voxel_bracket_mesh_tet.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)

end # voxel_bracket_mesh_tet

function voxel_bracket_mesh_tet_remesh()
    V = VoxelBoxVolume(Int, 8*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    h1 = solidcylinder((2.0, 2.5, 2.5), (1.0, 0.0, 0.0), 0.75)
    fillsolid!(V, differenceop(unionop(b1, b2), h1), 1)

    function checkvolumes(im)
        for i = 1:size(im.t, 1)
            if FinEtools.MeshTetrahedronModule.tetv1times6(im.v, im.t[i,1], im.t[i,2], im.t[i,3], im.t[i,4]) < 0.0
                println("Negative volume = $([im.t[i,1], im.t[i,2], im.t[i,3], im.t[i,4]])")
            end
        end
    end

    im = ImageMesher(V, zero(eltype(V.data)), eltype(V.data)[1])
    mesh!(im)
    checkvolumes(im)
    println("Mesh size: initial = $(size(im.t,1))")

    im.elementsizeweightfunctions = [ElementSizeWeightFunction(20.0, vec([0.0, 2.5, 2.5]), 1.0), ElementSizeWeightFunction(1.0, vec([0.0, 2.5, 2.5]), 3.5)]
    for i = 1:12
        println("Phase $i");
        checkvolumes(im)
        mesh!(im, 1.1)
        checkvolumes(im)
        println("Mesh size: final = $(size(im.t,1))")
        V = volumes(im)
        println("length(find(x -> x <= 0.0, V)) = $(length(find(x -> x <= 0.0, V)))")
        # open("im$(i)" * ".jls", "w") do file
        #     serialize(file, im)
        # end
    end

    fens = FENodeSet(im.v)
    fes = FESetT4(im.t)
    setlabel!(fes, im.tmid)

    # bfes = meshboundary(fes)
    # list = selectelem(fens, fes; overlappingbox = boundingbox([0.2018 2.1537 3.9064]), inflate = 0.01, allin = false)
    # File = "voxel_bracket_mesh_tet.vtk"
    # vtkexportmesh(File, fens, subset(fes, list))
    # @async run(`"paraview.exe" $File`)
    File = "voxel_bracket_mesh_tet.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)

    ne = NASTRANExporter("voxel_bracket_mesh_tet.nas")
    BEGIN_BULK(ne)
    for i = 1:count(fens)
        GRID(ne, i, fens.xyz[i, :])
    end
    for i = 1:count(fes)
        CTETRA(ne, i, 1, [fes.conn[i]...])
    end
    PSOLID(ne, 1, 1)
    MAT1(ne, 1, 20.0e9, 0.3, 1800.0)
    ENDDATA(ne)
    close(ne)

    stle = STLExporter("voxel_bracket_mesh_tet.stl")
    solid(stle)
    bfes = meshboundary(fes)
    for i = 1:count(bfes)
        facet(stle, fens.xyz[bfes.conn[i][1], :], fens.xyz[bfes.conn[i][2], :], fens.xyz[bfes.conn[i][3], :])
    end
    endsolid(stle)
    close(stle)
end # voxel_bracket_mesh_tet_remesh

function allrun()
    println("#####################################################")
    println("# voxel_bracket_mesh ")
    voxel_bracket_mesh()
    println("#####################################################")
    println("# voxel_bracket_mesh_tet ")
    voxel_bracket_mesh_tet()
    println("#####################################################")
    println("# voxel_bracket_mesh_tet_remesh ")
    voxel_bracket_mesh_tet_remesh()
    return true
end # function allrun

end # module voxel_bracket_examples


text:
# ------------------------------------------------------------------------------------------------------------------------------------------------
mutable struct Coords
    x::Integer
    y::Integer

    Coords() = new(0,0) # - такой конструктор определен для удобства использования
    Coords(x::Integer,y::Integer) = new(x,y)
end
# ------------------------------------------------------------------------------------------------------------------------------------------------

function move_Coords!(coords::Coords, side::HorizonSide)
    if side==Nord
        coords.y += 1
    elseif side==Sud
        coords.y -= 1
    elseif side==Ost
        coords.x += 1
    else #if side==West
        coords.x -= 1
    end
end

get_coords(coord::Coords) = (coord.x, coord.y)

#=
function get_coords(coord::Coords)
    return (coord.x, cord.y)
end
=#

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
# [[file:../fresa.org::gamsfmo][gamsfmo]]
function fmo(x::Array{Float64,1})
    open("gamsexample.gms", "w") do f
        write(f, "\$include gamsdeclarations.gms\n")
        write(f, "X1.fx = $(x[1]); \n")
        write(f, "X2.fx = $(x[2]); \n")
        write(f, "X3.fx = $(x[3]); \n")
        write(f, "solve TEST using NLP minimizing Z; \n")
        write(f, "file fresa /'gamsoutput.txt'/ ;\n")
        write(f, "put fresa ;\n")
        write(f, "put z.l /;\n")
        write(f, "put res.l /;\n")
        write(f, "put TEST.modelstat /;\n")
    end
    # execute GAMS
    run( `/opt/gams/latest/gams gamsexample.gms` )
    # read in results
    z = [0.0; 0.0]
    g = 0.0;
    open("gamsoutput.txt", "r") do f
        lines = readlines(f)
        z[1] = parse(Float64, lines[1])
        z[2] = abs(parse(Float64, lines[2]))
        modelstat = parse(Float64, lines[3])
        if modelstat != 1 && modelstat != 5
            g = 1
        end
    end
    # return results
    ( z, g )
end
# gamsfmo ends here

# [[file:../fresa.org::gamsfsingle][gamsfsingle]]
function fsingle(x::Array{Float64,1})
    open("gamsexample.gms", "w") do f
        write(f, "\$include gamsdeclarations.gms\n")
        write(f, "X1.fx = $(x[1]); \n")
        write(f, "X2.fx = $(x[2]); \n")
        write(f, "X3.fx = $(x[3]); \n")
        write(f, "solve TEST using NLP minimizing Z; \n")
        write(f, "file fresa /'gamsoutput.txt'/ ;\n")
        write(f, "put fresa ;\n")
        write(f, "put z.l /;\n")
        write(f, "put res.l /;\n")
        write(f, "put TEST.modelstat /;\n")
    end
    # execute GAMS
    run( `/opt/gams/latest/gams gamsexample.gms` )
    # read in results
    z = 0.0
    g = 0.0
    open("gamsoutput.txt", "r") do f
        lines = readlines(f)
        z = parse(Float64, lines[1])
        g = abs(parse(Float64, lines[2]))
        modelstat = parse(Float64, lines[3])
        if modelstat != 1 && modelstat != 5
            g = 10 # penalty function
        end
    end
    # return results
    ( z, g )
end
# gamsfsingle ends here

# [[file:../fresa.org::*solve the multi-objective problem using Fresa][solve the multi-objective problem using Fresa:1]]
using Fresa
domain = Fresa.Domain(x -> [0.0;0.0;0.0], x -> [5.0;3.0;3.0])
x0 = [4.0;2.0;2.0]
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,fmo)]
# now invoke Fresa to solve the problem
pareto, population = Fresa.solve(fmo, p0, domain;
                                 fitnesstype = :borda,
                                 ngen = 100)
println("Pareto front:")
println(population[pareto])
# solve the multi-objective problem using Fresa:1 ends here

# [[file:../fresa.org::*solve the multi-objective problem using Fresa][solve the multi-objective problem using Fresa:2]]
using PyPlot
z = [population[pareto[i]].z for i in 1:length(pareto)];
PyPlot.plot([z[i][1] for i=1:length(z)],
            [z[i][2] for i=1:length(z)],
            "ro")
PyPlot.savefig("gamsmo.pdf")
# solve the multi-objective problem using Fresa:2 ends here

# [[file:../fresa.org::*solve the single objective version][solve the single objective version:1]]
best, pop = Fresa.solve(fsingle, p0, domain; ngen = 100)
println("Population: $pop")
println("Best: f($(best.x)) = $(best.z), $( best.g )")
# solve the single objective version:1 ends here


text:
using Flux
using Statistics
using Plots
using DifferentialEquations

# Solve the ode given by $$ u' = \cos 2\pi t $$ by approximation of $u$ with a Neural Network
NNODE = Chain(x -> [x], # Take in a scalar and transform it into an array
              Dense(1,32,tanh),
              Dense(32,1),
              first) # Take first value, i.e. return a scalar
NNODE(1.0)

# construct function $g(t)$ to enforce initial condition $$ g(t) = t NNODE(t) + u(0)
u₀ = 1.0
g(t) = t*NNODE(t) + u₀

#Construct loss function to numericaly compute derivative and check against R.H.S.
ϵ = sqrt(eps(Float32))
loss() = mean(abs2(((g(t+ϵ)-g(t))/ϵ) - cos(2π*t)) for t in 0:1f-2:1f0)

opt = Flux.Descent(0.01)
data = Iterators.repeated((), 5000)
iter = 0
cb = function () #callback function to observe training
    global iter += 1
    if iter % 500 == 0
        display(loss())
    end
end
display(loss())
Flux.train!(loss, Flux.params(NNODE), data, opt; cb=cb)

# compare against true solution
t = 0:1e-3:1.0

plot(t, g.(t), label="NN")
plot!(t, u₀ .+ sin.(2π.*t)/2π, label="True Solution")


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03sd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    LDA = NMAX
    LDQG = NMAX
    LDWORK = NMAX*( NMAX+1 )
    A = Array{Float64,2}(undef, LDA,NMAX)
    QG = Array{Float64,2}(undef, LDQG,NMAX+1)
    WR = Array{Float64,1}(undef, NMAX)
    WI = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    JOBSCL =  vs[2][1]
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*(_jsz+1)÷2
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    _i0 = 0
    for j in 1:_jsz
       QG[j,j+1:_jsz+1] .= parsex.(Float64, vs[_i0+1:_i0+_jsz-j+1])
       _i0 += _jsz-j+1
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*(_jsz+1)÷2
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    _i0 = 0
    for j in 1:_jsz
       QG[j:_jsz,j] .= parsex.(Float64, vs[_i0+1:_i0+_jsz-j+1])
       _i0 += _jsz-j+1
    end
    close(f)
    # interp call 1

    INFO = SLICOT.mb03sd!(JOBSCL, N, A, QG, WR, WI, LDWORK)
    @test INFO == 0
    INFO != 0 && return

    println(io, "eigvals:")
    show(io, "text/plain", vcat(WR[1:N] + im*WI[1:N],
                 -WR[N:-1:1] - im*WI[N:-1:1]))
    println(io)

end # run_mb03sd()


text:
using MetabolicFluxAnalysis
using OpenCL

fluxi_1_f = open("testdata/cmfa_01.fluxi")
fluxi_1_s = readstring(fluxi_1_f)
fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

fluxi_target_semus = reduce( (x,y)->[x;y] , map( x->x.sass , fluxi1.measurements["Data_Exp_01"] ) )
emudecomp_1 = MetabolicFluxAnalysis.compute_emu_decomp( fluxi1.net , fluxi_target_semus )
emu_graph_1 = MetabolicFluxAnalysis.create_emu_graph(emudecomp_1)

fluxi_ast_1 = MetabolicFluxAnalysis.create_ast_from_emu_graph( emu_graph_1 )

target_smrs = MetabolicFluxAnalysis.construct_all_mass_ratios_for_measurement(fluxi1,"Data_Exp_01")
julia_fw_sim = MetabolicFluxAnalysis.ast_to_juliacode_01(fluxi1.net,fluxi_ast_1,target_smrs)

## OpenCL:

ocl_devices = OpenCL.cl.devices()
ocl_device  = ocl_devices[2]
ocl_context =  MetabolicFluxAnalysis.OCLContext(ocl_device)
print("OpenCL Device: $(ocl_device)")

#semu_buff = cl.Buffer(Float32, ocl_context.ocl_ctx, (:r, :copy), hostbuf=semu_data)
#flux_buff = cl.Buffer(Float32, ocl_context.ocl_ctx, (:r, :copy), hostbuf=flux_data)
#c_buff = cl.Buffer(Float32, ocl_context.ocl_ctx, :w, length(a))
ocl_fw_sim = MetabolicFluxAnalysis.ast_to_opencl_01(fluxi1.net,fluxi_ast_1,target_smrs)

p = OpenCL.cl.Program(ocl_context.ocl_ctx, source=ocl_fw_sim.code)
OpenCL.cl.build!( p , options="-cl-opt-disable",raise = true)
k = cl.Kernel(p, "fwsim")
##

constraints_config_01 = ["Constraints_FixedRatios","Constraints_BiomassRatios","Constraints_Measurements_01","Constraints_IrreversibilityPublication"]
constraints_01 = reduce( (x,y)->[x;y] , map( x -> fluxi1.linear_constraints[x] , constraints_config_01 ) )
fp1 = MetabolicFluxAnalysis.get_flux_polyhedron( fluxi1.net , constraints_01 , max_val=2000 )
# parametrize the flux polytope
rfd_parametrization = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( fp1 ; remove_red_constr=1)
@assert MetabolicFluxAnalysis.isinside( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt[:,:] )
print("Flux Polytop has dimension $(size(rfd_parametrization.p_opt.G,2))")
xx = MetabolicFluxAnalysis.sample_unif_hr( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt , 5000 , 1000 );
xx_flux = MetabolicFluxAnalysis.eval_opt2raw(rfd_parametrization,xx)

ilee_jul = MetabolicFluxAnalysis.create_julia_ile_evaluator( fluxi1.net , fluxi1.substrate_configs["SC_Conf_C1"] ,fluxi1.measurements["Data_Exp_01"] )

tic()
@profile sse_xx = MetabolicFluxAnalysis.eval_flux2sse( ilee_jul , xx_flux )
toc()

tic()
semu_data = convert( Array{Float32,1} , ilee_jul.substrate_mds[:] )
flux_data = convert( Array{Float32,1} , xx_flux[:] )
data_out  = convert( Array{Float32,1} , zeros( size(xx_flux,2)*length(ilee_jul.fwsim.target_mass_ratios) ) )

#device, ctx, queue = OpenCL.cl.create_compute_context()
semu_buff = cl.Buffer(Float32, ocl_context.ocl_ctx, (:r, :copy), hostbuf=semu_data)
flux_buff = cl.Buffer(Float32, ocl_context.ocl_ctx, (:r, :copy), hostbuf=flux_data)
out_buff  = cl.Buffer(Float32, ocl_context.ocl_ctx, :w, length(data_out))
#k = ocl_fwsim.k.value #cl.Kernel(p, "fwsim")

req_sx = 256
req_sy = 256

println("length data_out = $(length(data_out))")
println("nfluxes = $(size(xx_flux,2))")
println("run config: $(req_sx),$(req_sy*size(xx_flux,2)) ; $(req_sx),$(req_sy)")

#queue(k, [ size(flux_data,2) * req_sx ; req_sy ], [req_sx;req_sy], semu_buff, flux_buff, out_buff)
ocl_context.ocl_queue(k, [ req_sx ; req_sy * size(xx_flux,2) ], [req_sx;req_sy], semu_buff, flux_buff, out_buff)

r = cl.read(queue, out_buff)

#sse_xx = MetabolicFluxAnalysis.eval_flux2sse( ilee_1 , xx_flux )
toc()


text:
function WeightedExplore(A,b,C,d,XLB,XUB,VariableDef,YR,TimeLim = 1e10)
        const EPS = 1e-5;
        t0 = time_ns();
        nObj = length(YR)

        MMModel= WeightedMinMaxModel(A,b,C,d,XLB,XUB,VariableDef,YR)

        # Time commands ---------------------------------------------------
        if (TimeLim-(time_ns() - t0)/1e9) < 0.001;
            return [],[],[],false
        end
        # -----------------------------------------------------------------

        setsolver(MMModel,CplexSolver(CPXPARAM_Threads = 1,CPXPARAM_ScreenOutput = 0,CPXPARAM_TimeLimit = TimeLim-(time_ns() - t0)/1e9));
        status = solve(MMModel);

        # Time commands ---------------------------------------------------
        if (status != :Optimal) | (TimeLim-(time_ns() - t0)/1e9 < 0.001)
            return [],[],[],false
        end
        # -----------------------------------------------------------------

        Var = getvalue(getVar(MMModel,:X))
        Y_bar1 = zeros(size(YR))
        for j = 1:nObj
            Y_bar1[j] = sum(Var.*C[j,:]) + d[j];
        end

        # println(Y_bar2)

        Y = copy(Y_bar1)
        if any(abs(YR-Y_bar1) .< EPS)
            Y = YR + maximum(abs(Y_bar1-YR))
        end

        return Y_bar1,Var,Y,true
end


text:
using Test

using SingleColumnModels.FiniteDifferenceGrids
using SingleColumnModels.StateVecs
using LinearAlgebra

output_root = joinpath("..", "output", "tests", "PDEs")

include("HeatEquation.jl")
include("AdvectionEquation.jl")


text:
@testset "Normalize LUT" begin
    lut = SVector(complex.(ones(4,4), ones(4,4)) * 3, complex.(ones(4,4), ones(4,4)) * 2)
    lut_array = [lut[i][j,k] for i = 1:length(lut), j = 1:size(lut[1],1), k = 1:size(lut[1],2)]
    normalized_lut = @inferred PhasedArray.norm_manifold(lut_array)
    lut_test = SVector(complex.(ones(4,4), ones(4,4)) * 3 / sqrt(26 / 2), complex.(ones(4,4), ones(4,4)) * 2 / sqrt(26 / 2))
    lut_array_test = [lut_test[i][j,k] for i = 1:length(lut_test), j = 1:size(lut_test[1],1), k = 1:size(lut_test[1],2)]
    @test normalized_lut ≈ lut_array_test
end

@testset "LUT expansion" begin
    lut = SVector(randn(ComplexF64, 10, 10), randn(ComplexF64, 10, 10), randn(ComplexF64, 10, 10), randn(ComplexF64, 10, 10))
    lut_array = [lut[i][j,k] for i = 1:length(lut), j = 1:size(lut[1],1), k = 1:size(lut[1],2)]
    lut_expanded = @inferred PhasedArray.expand(lut_array, 2)
    num_ants = length(lut)
    @test size(lut_expanded) == (num_ants, 14, 14)
    for i = 1:num_ants
        @test lut_expanded[i,:,1] == lut_expanded[i,:,11]
        @test lut_expanded[i,:,2] == lut_expanded[i,:,12]
        @test lut_expanded[i,:,3] == lut_expanded[i,:,13]
        @test lut_expanded[i,:,4] == lut_expanded[i,:,14]
        @test lut_expanded[i,1,3:12] == circshift(lut_expanded[i,5,3:12], 5)
        @test lut_expanded[i,2,3:12] == circshift(lut_expanded[i,4,3:12], 5)
        @test lut_expanded[i,14,3:12] == circshift(lut_expanded[i,10,3:12], 5)
        @test lut_expanded[i,13,3:12] == circshift(lut_expanded[i,11,3:12], 5)
    end

    @test PhasedArray.calc_expansion_length(Constant) == 0
    @test PhasedArray.calc_expansion_length(Constant) == 0
    @test PhasedArray.calc_expansion_length(Quadratic) == 18
end

@testset "LUT manifold" begin
    λ = 0.1904
    ant_positions = SVector(λ / 4 * SVector(1, 1, 0), λ / 4 * SVector(-1, 1, 0), λ / 4 * SVector(1, -1, 0), λ / 4 * SVector(-1, -1 , 0))
    num_ants = length(ant_positions)
    # Simulated measured pattern
    antenna_gain = 2
    sph2cart = CartesianFromSpherical()
    lut = map(ant_position -> [antenna_gain * cis(2 * π / λ * sph2cart(Spherical(1.0, θ * π / 180, π / 2 - ϕ * π / 180))' * ant_position) for ϕ=0:90, θ=0:359], ant_positions)
    manifold = @inferred RealManifold(lut, max_elevation = π / 2)
    @test typeof(manifold) <: AbstractManifold{4}
    @test norm(@inferred(get_steer_vec(manifold, SVector(1, 1, 1)))) ≈ sqrt(num_ants)
    @test @inferred(get_steer_vec(manifold, SVector(0,0,1))) ≈ [1, 1, 1, 1]
    @test @inferred(get_steer_vec(manifold, SVector(0,1,0))) ≈ [1im, 1im, -1im, -1im]
    @test @inferred(get_steer_vec(manifold, SVector(1,0,0))) ≈ [1im, -1im, 1im, -1im]

    @test @inferred(get_steer_vec(manifold, Spherical(1,0,0))) ≈ [1im, -1im, 1im, -1im]
    @test @inferred(get_steer_vec(manifold, Spherical(1,0,0), RotXYZ(0.0,0.0,π/2))) ≈ [1im, 1im, -1im, -1im]
    @test @inferred(get_steer_vec(manifold, SVector(1,0,0), RotXYZ(0.0,0.0,π/2))) ≈ [1im, 1im, -1im, -1im]

    manifold = @inferred RealManifold(lut[1], lut[2], lut[3], lut[4], max_elevation = π / 2)
    @test typeof(manifold) <: AbstractManifold{4}
    @test norm(@inferred(get_steer_vec(manifold, SVector(1, 1, 1)))) ≈ sqrt(num_ants)
    @test @inferred(get_steer_vec(manifold, SVector(0,0,1))) ≈ [1, 1, 1, 1]
    @test @inferred(get_steer_vec(manifold, SVector(0,1,0))) ≈ [1im, 1im, -1im, -1im]
    @test @inferred(get_steer_vec(manifold, SVector(1,0,0))) ≈ [1im, -1im, 1im, -1im]

    @test @inferred(get_steer_vec(manifold, Spherical(1,0,0))) ≈ [1im, -1im, 1im, -1im]
    @test @inferred(get_steer_vec(manifold, Spherical(1,0,0), RotXYZ(0.0,0.0,π/2))) ≈ [1im, 1im, -1im, -1im]
    @test @inferred(get_steer_vec(manifold, SVector(1,0,0), RotXYZ(0.0,0.0,π/2))) ≈ [1im, 1im, -1im, -1im]
end

@testset "LUT manifold interpolated" begin
    λ = 0.1904
    ant_positions = SVector(λ / 4 * SVector(1, 1, 0), λ / 4 * SVector(-1, 1, 0), λ / 4 * SVector(1, -1, 0), λ / 4 * SVector(-1, -1 , 0))
    num_ants = length(ant_positions)
    # Simulated measured pattern
    antenna_gain = 2;
    sph2cart = CartesianFromSpherical()
    lut = map(ant_position -> [antenna_gain * cis(2 * π / λ * sph2cart(Spherical(1.0, θ * π / 180, π / 2 - ϕ * π / 180))' * ant_position) for ϕ=0:90, θ=0:359], ant_positions)
    manifold = @inferred RealManifold(lut, Linear, max_elevation = π / 2)
    @test typeof(manifold) <: AbstractManifold{4}

    @test isapprox(norm(@inferred(get_steer_vec(manifold, SVector(1, 1, 1)))), sqrt(num_ants); rtol = 1e-4)
    @test @inferred(get_steer_vec(manifold, SVector(0,0,1))) ≈ [1, 1, 1, 1]
    @test @inferred(get_steer_vec(manifold, SVector(0,1,0))) ≈ [1im, 1im, -1im, -1im]
    @test @inferred(get_steer_vec(manifold, SVector(1,0,0))) ≈ [1im, -1im, 1im, -1im]
end

@testset "Example LUT manifold" begin
    num_ants = length(EXAMPLE_LUT)
    manifold = RealManifold(EXAMPLE_LUT)
    @test typeof(manifold) <: AbstractManifold{2}
    @test norm(@inferred(get_steer_vec(manifold, SVector(0,0,1)))) < sqrt(num_ants)
end


text:
# An implementation of LSQR for the solution of the
# over-determined linear least-squares problem
#
#  minimize ‖Ax - b‖
#
# equivalently, of the normal equations
#
#  A'Ax = A'b.
#
# LSQR is formally equivalent to applying the conjugate gradient method
# to the normal equations but should be more stable. It is also formally
# equivalent to CGLS though LSQR should be expected to be more stable on
# ill-conditioned or poorly scaled problems.
#
# This implementation follows the original implementation by
# Michael Saunders described in
#
# C. C. Paige and M. A. Saunders, LSQR: An Algorithm for Sparse Linear
# Equations and Sparse Least Squares, ACM Transactions on Mathematical
# Software, 8(1), pp. 43--71, 1982.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Montreal, QC, May 2015.

export lsqr

"""Solve the regularized linear least-squares problem

  minimize ‖b - Ax‖₂² + λ² ‖x‖₂²

using the LSQR method, where λ ≥ 0 is a regularization parameter.
LSQR is formally equivalent to applying CG to the normal equations

  (A'A + λ² I) x = A'b

(and therefore to CGLS) but is more stable.

LSQR produces monotonic residuals ‖r‖₂ but not optimality residuals ‖A'r‖₂.
It is formally equivalent to CGLS, though can be slightly more accurate.

Preconditioners M and N may be provided in the form of linear operators and are
assumed to be symmetric and positive definite. If `sqd` is set to `true`,
we solve the symmetric and quasi-definite system

  [ E   A' ] [ r ]   [ b ]
  [ A  -F  ] [ x ] = [ 0 ],

where E = M⁻¹  and F = N⁻¹.

If `sqd` is set to `false` (the default), we solve the symmetric and
indefinite system

  [ E   A' ] [ r ]   [ b ]
  [ A   0  ] [ x ] = [ 0 ].

In this case, `N` can still be specified and indicates the norm
in which `x` should be measured.
"""
function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(A,1)),
              N :: AbstractLinearOperator=opEye(size(A,2)),
              sqd :: Bool=false,
              λ :: Float64=0.0, axtol :: Float64=1.0e-8, btol :: Float64=1.0e-8,
              atol :: Float64=0.0, rtol :: Float64=0.0,
              etol :: Float64=1.0e-8, window :: Int=5,
              itmax :: Int=0, conlim :: Float64=1.0e+8,
              radius :: Float64=0.0, verbose :: Bool=false) where T <: Number

  m, n = size(A)
  size(b, 1) == m || error("Inconsistent problem size")
  verbose && @printf("LSQR: system of %d equations in %d variables\n", m, n)

  # If solving an SQD system, set regularization to 1.
  sqd && (λ = 1.0)
  λ² = λ * λ
  ctol = conlim > 0.0 ? 1/conlim : 0.0
  x = zeros(T, n)

  # Initialize Golub-Kahan process.
  # β₁ M u₁ = b.
  Mu = copy(b)
  u = M * Mu
  β₁ = sqrt(@kdot(m, u, Mu))
  β₁ == 0.0 && return (x, SimpleStats(true, false, [0.0], [0.0], "x = 0 is a zero-residual solution"))
  β = β₁

  @kscal!(m, 1.0/β₁, u)
  @kscal!(m, 1.0/β₁, Mu)
  Nv = copy(A' * u)
  v = N * Nv
  α = sqrt(@kdot(n, v, Nv))
  Anorm² = α * α
  Anorm = α
  Acond  = 0.0
  xNorm  = 0.0
  xNorm² = 0.0
  dNorm² = 0.0
  var = zeros(T, n)
  c2 = -1.0
  s2 =  0.0
  z  =  0.0

  xENorm² = 0.0
  err_lbnd = 0.0
  err_vec = zeros(T, window)

  verbose && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s\n",
                     "Aprod", "‖r‖", "‖A'r‖", "β", "α", "cos", "sin", "‖A‖²")
  verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e\n",
                     1, β₁, α, β₁, α, 0, 1, Anorm²)

  # A'b = 0 so x = 0 is a minimum least-squares solution
  α == 0.0 && return (x, SimpleStats(true, false, [β₁], [0.0], "x = 0 is a minimum least-squares solution"))
  @kscal!(n, 1.0/α, v)
  @kscal!(n, 1.0/α, Nv)
  w = copy(v)

  # Initialize other constants.
  ϕbar = β₁
  ρbar = α
  # θ = 0.0
  rNorm = ϕbar
  r1Norm = rNorm
  r2Norm = rNorm
  res2   = 0.0
  rNorms = [r2Norm]
  ArNorm = ArNorm0 = α * β
  ArNorms = [ArNorm]

  iter = 0
  itmax == 0 && (itmax = m + n)

  status = "unknown"
  on_boundary = false
  solved_lim = ArNorm / (Anorm * rNorm) <= axtol
  solved_mach = 1.0 + ArNorm / (Anorm * rNorm) <= 1.0
  solved = solved_mach | solved_lim
  tired  = iter >= itmax
  ill_cond = ill_cond_mach = ill_cond_lim = false
  zero_resid_lim = rNorm / β₁ <= axtol
  zero_resid_mach = 1.0 + rNorm / β₁ <= 1.0
  zero_resid = zero_resid_mach | zero_resid_lim
  fwd_err = false

  while ! (solved || tired || ill_cond)
    iter = iter + 1

    # Generate next Golub-Kahan vectors.
    # 1. βu = Av - αu
    @kscal!(m, -α, Mu)
    @kaxpy!(m, 1.0, A * v, Mu)
    u = M * Mu
    β = sqrt(@kdot(m, u, Mu))
    if β != 0.0
      @kscal!(m, 1.0/β, u)
      @kscal!(m, 1.0/β, Mu)
      Anorm² = Anorm² + α * α + β * β;  # = ‖B_{k-1}‖²
      λ > 0.0 && (Anorm² += λ²)

      # 2. αv = A'u - βv
      @kscal!(n, -β, Nv)
      @kaxpy!(n, 1.0, A' * u, Nv)
      v = N * Nv
      α = sqrt(@kdot(n, v, Nv))
      if α != 0.0
        @kscal!(n, 1.0/α, v)
        @kscal!(n, 1.0/α, Nv)
      end
    end

    # Continue QR factorization
    # 1. Eliminate the regularization parameter.
    (c1, s1, ρbar1) = sym_givens(ρbar, λ)
    ψ = s1 * ϕbar
    ϕbar = c1 * ϕbar

    # 2. Eliminate β.
    # Q [ Lₖ  β₁ e₁ ] = [ Rₖ   zₖ  ] :
    #   [ β    0    ]   [ 0   ζbar ]
    #
    #       k  k+1    k    k+1      k  k+1
    # k   [ c   s ] [ ρbar    ] = [ ρ  θ⁺    ]
    # k+1 [ s  -c ] [ β    α⁺ ]   [    ρbar⁺ ]
    #
    # so that we obtain
    #
    # [ c  s ] [ ζbar ] = [ ζ     ]
    # [ s -c ] [  0   ]   [ ζbar⁺ ]
    (c, s, ρ) = sym_givens(ρbar1, β)
    ϕ = c * ϕbar
    ϕbar = s * ϕbar

    xENorm² = xENorm² + ϕ * ϕ
    err_vec[mod(iter, window) + 1] = ϕ
    iter >= window && (err_lbnd = norm(err_vec))

    τ = s * ϕ;
    θ = s * α
    ρbar = -c * α

    d =  w / ρ;  # TODO: Use BLAS call.
    dNorm² += @kdot(n, d, d)
    var += d .* d

    # if a trust-region constraint is give, compute step to the boundary
    # the step ϕ/ρ is not necessarily positive
    σ = ϕ / ρ
    if radius > 0.0
      t1, t2 = to_boundary(x, w, radius)
      tmax, tmin = max(t1, t2), min(t1, t2)
      on_boundary = σ > tmax || σ < tmin
      σ = σ > 0 ? min(σ, tmax) : max(σ, tmin)
    end

    @kaxpy!(n,  σ, w, x)  # x = x + ϕ / ρ * w
    @kscal!(n, -θ/ρ, w)
    @kaxpy!(n,  1.0, v, w)  # w = v - θ / ρ * w

    # Use a plane rotation on the right to eliminate the super-diagonal
    # element (θ) of the upper-bidiagonal matrix.
    # Use the result to estimate norm(x).
    δ = s2 * ρ
    γbar = -c2 * ρ
    rhs = ϕ - δ * z
    zbar = rhs / γbar
    xNorm = sqrt(xNorm² + zbar * zbar)
    (c2, s2, γ) = sym_givens(γbar, θ)
    z = rhs / γ
    xNorm² += z * z

    Anorm = sqrt(Anorm²)
    Acond = Anorm * sqrt(dNorm²)
    res1  = ϕbar * ϕbar
    res2 += ψ * ψ
    rNorm = sqrt(res1 + res2)

    ArNorm = α * abs(τ)
    push!(ArNorms, ArNorm)

    r1sq = rNorm * rNorm - λ² * xNorm²
    r1Norm = sqrt(abs(r1sq))
    r1sq < 0.0 && (r1Norm = -r1Norm)
    r2Norm = rNorm
    push!(rNorms, r2Norm)

    test1 = rNorm / β₁
    test2 = ArNorm / (Anorm * rNorm)
    test3 = 1 / Acond
    t1    = test1 / (1.0 + Anorm * xNorm / β₁)
    rNormtol = btol + axtol * Anorm * xNorm / β₁

    verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e\n",
                       1 + 2 * iter, rNorm, ArNorm, β, α, c, s, Anorm²)

    # Stopping conditions that do not depend on user input.
    # This is to guard against tolerances that are unreasonably small.
    ill_cond_mach = (1.0 + test3 <= 1.0)
    solved_mach = (1.0 + test2 <= 1.0)
    zero_resid_mach = (1.0 + t1 <= 1.0)

    # Stopping conditions based on user-provided tolerances.
    tired  = iter >= itmax
    ill_cond_lim = (test3 <= ctol)
    solved_lim = (test2 <= axtol)
    solved_opt = ArNorm <= atol + rtol * ArNorm0
    zero_resid_lim = (test1 <= rNormtol)
    iter >= window && (fwd_err = err_lbnd <= etol * sqrt(xENorm²))

    ill_cond = ill_cond_mach | ill_cond_lim
    solved = solved_mach | solved_lim | solved_opt | zero_resid_mach | zero_resid_lim | fwd_err | on_boundary
  end

  tired         && (status = "maximum number of iterations exceeded")
  ill_cond_mach && (status = "condition number seems too large for this machine")
  ill_cond_lim  && (status = "condition number exceeds tolerance")
  solved        && (status = "found approximate minimum least-squares solution")
  zero_resid    && (status = "found approximate zero-residual solution")
  fwd_err       && (status = "truncated forward error small enough")
  on_boundary   && (status = "on trust-region boundary")

  stats = SimpleStats(solved, !zero_resid, rNorms, ArNorms, status)
  return (x, stats)
end


text:
# Shannon entropy of discrete distribution with probabilities given by p
entropy(p::Real) = iszero(p) ? p : -p*log2(p)
entropy(p::AbstractVector) = sum(entropy, p)
perplexity(p::AbstractVector) = 2^entropy(p)

# Kullback-Leibler divergence
kl(p::Real, q::Real) = iszero(p) ? zero(p) : p*log(p/q) # what about q = 0?
kl(pq::NTuple{2}) = kl(pq...) # unpack tuple argument
# Z is assumed to be normalization constant of Q
function kl(P::AbstractMatOrFac, Q::AbstractMatOrFac, Z::Real = 1)
    val = zero(promote_type(eltype(P), eltype(Q)))
    for i in 1:size(P, 1), j in 1:size(Q, 2)
        val += kl(P[i, j], Q[i, j] / Z)
    end
    return val
end
# specialization for sparse matrix P
# Z is assumed to be normalization constant of Q
function kl(P::SparseMatrixCSC, Q::AbstractMatOrFac, Z::Real = 1)
    I, J, _ = findnz(P)
    val = zero(promote_type(eltype(P), eltype(Q)))
    for (i, j) in zip(I, J) # IDEA: parallelize, simd?
        val += kl(P[i, j], Q[i, j] / Z)
    end
    return val
end

# euclidean distance
_d2(x::Real, y::Real) = (x-y)^2
_d2(x::NTuple{2}) = _d2(x...)
_euclidean(x, y) = sqrt(sum(_d2, zip(x, y)))
euclidean(x::AbstractVector, y::AbstractVector) = _euclidean(x, y) # avoids temporary array
euclidean(x::NTuple, y::NTuple) = _euclidean(x, y) # avoids temporary array

# kernels
gaussian(x::Real) = exp(-x^2/2)
cauchy(x::Real) = inv(1+x^2)

# linear algebra
function hadamard_product(P::AbstractMatrix, Q::AbstractMatOrFac)
    n, m = size(P)
    H = zeros(eltype(P), n, m)
    for j in 1:n, i in 1:m
        H[i, j] = P[i, j] * Q[i, j]
    end
    return H
end

function hadamard_product(P::SparseMatrixCSC, Q::AbstractMatOrFac)
    PQ = similar(P) # this copies the sparsity pattern of P, therefore Q is nevery instantiated densely
    I, J, _ = findnz(P)
     for (i, j) in zip(I, J)
          PQ[i, j] = P[i, j] * Q[i, j]
     end
     return PQ
end

# synthetic data generators
function two_bump_data(d::Int, n::Int)
    X = randn(d, n)
    r = rem(n, 2)
    n = div(n, 2)
    @. X[:, 1:n] = .1X[:, 1:n] + 1
    @. X[:, n+1:end] = .1X[:, n+1:end] - 1
    return X
end

############################# optimization #####################################
function newtons_method(objective, x; tol::Real = 1e-6, max_iter::Int = 32)
    g(x) = ForwardDiff.derivative(objective, x)
    for t in 1:max_iter # Newton's method
        ft, gt = value_derivative(objective, x)
        gt = sign(gt) * maximum(abs, (ft, gt)) # results in abs(ft/gt) ≤ 1 which leads to more robust convergence
        abs(ft) > tol || break
        x -= ft / gt
    end
    return x
end

# auto-diff helper functions to calculate value and gradient with one call
function value_derivative(f, x::Real, diffresults::Val{true} = Val(true))
    r = DiffResults.DiffResult(zero(x), zero(x))
    r = ForwardDiff.derivative!(r, f, x)
    r.value, r.derivs[1]
end
function value_derivative(f, x::Real, diffresults::Val{false})
    f(x), ForwardDiff.derivative(f, x)
end

############################### neighborhoods ##################################
# defaults to euclidean distance
function k_nearest_neighbors(X::AbstractMatrix, k::Int)
    metric = Euclidean()
    tree = BallTree(X, metric; leafsize = 32, reorder = true)
    neighbors, distances = knn(tree, X, k + 1) # plus one because we are putting in the original points as well
    for i in 1:size(X, 2) # delete self-references
         filter!(!=(i), neighbors[i])
    end
    return neighbors
end

# ε = sqrt(-2log(δ)) # based on Gaussian kernel, minimum distance for kernel to be smaller than δ
function in_range_neighbors(X::AbstractMatrix, r::Real)
    reorder = true # reorder copies the data and puts nearby points close in memory, think about using Distance.jl for the rest of the code too
    leafsize = 8 # to optimize
    metric = Euclidean()
    tree = BallTree(X, metric; leafsize = 8, reorder = true)
    neighbors = inrange(tree, X, r) # plus one because we are putting in the original points as well
    for i in 1:size(X, 2) # delete self-references
         filter!(!=(i), neighbors[i])
    end
    return neighbors
end


text:
using HTTPClient.HTTPC
#using Dates
using Gumbo

# =====================================================================
function fetchHistoricalData(sym::String, fromDate::Dates.Date,
                             toDate=Dates.today()::Dates.Date, freq="d")

    fromY = Dates.year(fromDate);
    fromM = Dates.month(fromDate);
    fromD = Dates.day(fromDate);

    toY = Dates.year(toDate);
    toM = Dates.month(toDate);
    toD = Dates.day(toDate);

    url1 = "http://ichart.yahoo.com/table.csv?s=$sym&a=$(fromM-1)&b=$(fromD)"
    url2 = "&c=$(fromY)&d=$(toM-1)&e=$(toD)&f=$(toY)&g=$freq&ignore=.csv"
    url = "$url1$url2"

    page = HTTPC.get(url);
    text = String(page.body.data)

    if (text[1] == '<')
        error("Bad historical data pull for $sym")
    end

    table = readdlm(IOBuffer(text), ',', Any)

    HD = HistoricalData(sym, [Dates.Date(s) for s in table[end:-1:2,1]],
                        float64(table[end:-1:2,2]),
                        float64(table[end:-1:2,3]),
                        float64(table[end:-1:2,4]),
                        float64(table[end:-1:2,5]),
                        float64(table[end:-1:2,6]),
                        float64(table[end:-1:2,7]))

    return HD
end

# =====================================================================
function fetchQuotes(syms, props)
    urlbase = "http://download.finance.yahoo.com/d/quotes.csv?s="
    if (isa(syms, Array))
        symspart = "";
        for i=1:length(syms)-1
            symspart = "$symspart$(syms[i]),"
        end
        symspart = "$symspart$(syms[length(syms)])"
    else
        symspart = syms;
    end

    if (isa(props, Array))
        propspart = "";
        for i=1:length(props)
            propspart = "$propspart$(props[i])"
        end
    else
        propspart = props
    end

    url = "$urlbase$symspart&f=$propspart&e=.csv"

    page = HTTPC.get(url)
    text = String(page.body.data)
    if (text[1] == '<')
        error("Bad quotes pull")
    end

    table = readdlm(IOBuffer(text), ',', Any)

    return table;
end

# =====================================================================
function fetchSectors(sortBy::String = "coname", sortDir::String = "u")

    url = "http://biz.yahoo.com/p/csv/s_$sortBy$sortDir.csv"

    page = HTTPC.get(url)
    text = String(page.body.data)
    if (text[1] == '<')
        error("Bad sector pull")
    end

    table = readdlm(IOBuffer(text), ',', Any);

    # table from Yahoo is "null terminated", so return all but last line
    return table[1:end-1,:];
end

# =====================================================================
function fetchIndustries(sectorNum, sortBy::String = "coname", sortDir::String = "u")

    url = "http://biz.yahoo.com/p/csv/$(sectorNum)$sortBy$sortDir.csv";
    println(url)
    page = HTTPC.get(url)
    text = String(page.body.data)
    if (text[1] == '<')
        error("Bad Industries pull")
    end

    table = readdlm(IOBuffer(text), ',', Any);

    return table[1:end-1,:];
end

# =====================================================================
function fetchCompanies(industryNum, sortBy::String = "coname", sortDir::String = "u")

    url = "http://biz.yahoo.com/p/csv/$(industryNum)$sortBy$sortDir.csv";
    println(url)
    page = HTTPC.get(url)
    text = String(page.body.data)
    if (text[1] == '<')
        error("Bad Companies pull")
    end

    table = readdlm(IOBuffer(text), ',', Any);

    return table[1:end-1,:];
end

# =====================================================================
function fetchRSS(symbol::String, locale::String = "en-us")

    url = "http://feeds.finance.yahoo.com/rss/2.0/headline?s=$symbol&region=$(locale[4:end])&lang=$locale"

    page = HTTPC.get(url)
    text = String(page.body.data)

    # parse html with gumbo package
    doc = parsehtml(text)

    # prepare empty Array{Any,1} for return
    articles = []

    # Parse returned html document
    chan = doc.root.children[2].children[1].children[1]
    for ch in chan.children
        if(isa(ch, Gumbo.HTMLElement{:item}))
            # This condition handles the event where no page is found
            if (isa(ch.children[1],Gumbo.HTMLElement{:title}) && isa(ch.children[3],HTMLText))
                title = ch.children[1].children[1].text
                link = ch.children[3].text
                push!(articles, (title, link))
            end
        end
    end

    return articles;
end


text:
function f(b,q,n,nl,po,pf,c,s) println(s,nl,po,q,b,b,q,c,q,b,q,q,c,q,n,q,c,q,b,n,q,c,q,po,q,c,q,pf,q,c,q,c,q,c,q,s,q,pf) end
f("\\","\"","n","\n","f(",")",",","function f(b,q,n,nl,po,pf,c,s) println(s,nl,po,q,b,b,q,c,q,b,q,q,c,q,n,q,c,q,b,n,q,c,q,po,q,c,q,pf,q,c,q,c,q,c,q,s,q,pf) end")

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
__precompile__()
module Ising
using ..Random
using DataStructures
using Erdos

export  ground_state_mincut, ground_state_ϵgreedy, ground_state_τeo,
        random_couplings

getJ(J::T, i, j, k) where {T} = J
getJ(J::Vector{Vector{T}}, i, j, k) where {T} = J[i][k]
getJ(J::AbstractMatrix{T}, i, j, k) where {T} = J[i,j]

Jtype(J::T) where {T} = T
Jtype(J::Vector{Vector{T}}) where {T} = T
Jtype(J::AbstractMatrix{T}) where {T} = T

function energy(g::AGraph, σ::Vector, h::Vector{TH}, Js) where TH
    E = zero(promote_type(Jtype(Js), TH))
    for i=1:nv(g)
        E -= h[i]*σ[i]
        for (k,j) in enumerate(neighbors(g,i))
            if i < j
                E -= getJ(Js, i, j, k) * σ[i] * σ[j]
            end
        end
    end
    return E
end

function random_couplings(g::AGraph)
    adj = adjacency_list(g)
    coupls = [ones(Int, degree(g, i)) for i=1:nv(g)]
    for i=1:nv(g)
        for (k,j) in enumerate(neighbors(g,i))
            if i < j
                J = rand([-1,1])
                coupls[i][k] = J
                ki = findfirst(adj[j], i)
                coupls[j][ki] = J
            end
        end
    end
    return coupls
end

function graph_and_couplings(Js::Vector{Vector{Int}})
    n = length(Js)
    g = Graph(n)
    for i=1:n
        neigs = find(i->i!=0, Js[i])
        for j in neigs
            if i < j
                add_edge!(g, i, j)
            end
        end
    end

    myJ = [zeros(Int,degree(g,i)) for i=1:nv(g)]
    for i=1:n
        for (k,j) in enumerate(neighbors(g,i))
            myJ[i][k] = Js[i][j]
        end
    end
    return g, myJ
end

include("mincut.jl")
include("tau-eo.jl")

end # submodule

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


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

@testset ExtendedTestSet "DatasetRegistry" begin
    reg = Datasets.DatasetRegistry()
    @testset ExtendedTestSet "registerdataset!" begin
        @test_nowarn Datasets.registerdataset!(
            reg, "mnist_var_size_tiny", () -> datasetpath("mnist_var_size_tiny"))

        # Reregistering should error
        @test_throws ErrorException Datasets.registerdataset!(
            reg, "mnist_var_size_tiny", () -> datasetpath("mnist_var_size_tiny"))
    end

    @testset ExtendedTestSet "listdatasources" begin
        @test length(Datasets.listdatasources(reg)) == 1
    end

    @testset ExtendedTestSet "datasetpath" begin
        @test_nowarn datasetpath(reg, "mnist_var_size_tiny")
    end

    @testset ExtendedTestSet "registerrecipe!" begin
        @test_nowarn Datasets.registerrecipe!(
            reg, "mnist_var_size_tiny", Datasets.ImageFolders())
    end

    @testset ExtendedTestSet "finddatasets" begin
        @test finddatasets(reg) |> length == 1
        @test finddatasets(reg, name="mnist_var_size_tiny") |> length == 1
        @test finddatasets(reg, blocks=Tuple{Image, Label}) |> length == 1
        @test finddatasets(reg, blocks=Tuple{Image, LabelMulti}) |> length == 0
        @test finddatasets(reg, name="mnist_var_size_tiny", blocks=Tuple{Image, Label}) |> length == 1
        @test finddatasets(reg, name="mnist", blocks=Tuple{Image, Label}) |> length == 0
    end
end


text:
module DustExtinction

using Unitful, UnitfulAstro, DataDeps

export redden,
       deredden,
       ccm89,
       cal00,
       od94,
       SFD98Map,
       ebv_galactic

include("color_laws.jl")
include("dust_maps.jl")

@deprecate ccm89(x::AbstractArray, r_v::Real = 3.1) ccm89.(x, r_v)
@deprecate od94(x::AbstractArray, r_v::Real = 3.1) od94.(x, r_v)

#--------------------------------------------------------------------------------

# reddening functions
"""
    redden(f::Real, λ::Real, Av; Rv=3.1, law=ccm89)
    redden(f::Quantity, λ::Quantity, Av; Rv=3.1, law=ccm89)

Redden the value `f` by the value calculated via the given law and total
extinction value `Av`. By default we use `Rv=3.1` which is the Milky Way
average selective attenuation. Note that λ should be in Angstrom if it is not
a `Quantity`.
"""
redden(f::Real, λ::Real, Av::Real; Rv = 3.1, law = ccm89) = f * 10^(-0.4 * Av * law(λ, Rv))
redden(f::Quantity, λ::Quantity, Av::Real; Rv = 3.1, law = ccm89) = f * (Av * law(λ, Rv))

"""
    deredden(f::Real, λ::Real, Av; Rv=3.1, law=ccm89)
    deredden(f::Quantity, λ::Quantity, Av; Rv=3.1, law=ccm89)

Deredden the value `f` by the value calculated via the given law and total
extinction value `Av`. By default we use `Rv=3.1` which is the Milky Way
average selective attenuation. Note that λ should be in Angstrom if it is not
a `Quantity`.
"""
deredden(f::Real, λ::Real, Av::Real; Rv = 3.1, law = ccm89) = f / 10^(-0.4 * Av * law(λ, Rv))
deredden(f::Quantity, λ::Quantity, Av::Real; Rv = 3.1, law = ccm89) = f / (Av * law(λ, Rv))

#--------------------------------------------------------------------------------

function __init__()
    # register our data dependencies
    register(DataDep("sfd98_map",
    """
    SFD98 Galactic Dust Maps
    Website: https://sncosmo.github.io
    """,
    ["https://sncosmo.github.io/data/dust/SFD_dust_4096_ngp.fits",
    "https://sncosmo.github.io/data/dust/SFD_dust_4096_sgp.fits"],
    ["50b6aaad0b880762d0fd081177802dcc17c39d7044a410dd5649e2dfd0503e97",
    "84891a59054adab44a7be54051e4dcf0e66e3f13eee0d845ce3739242f553b83"]))
end

end # module


text:
using SemiseparableMatrices, BandedMatrices, LinearAlgebra, Test

include("test_semiseparable.jl")
include("test_almostbanded.jl")
include("test_invbanded.jl")


text:
function multidwellplans(stockdata)
    # Find other dwellings within same footprint
    undersameroof(uprns,osgbs,osgb) = uprns[findall(osgbs.==osgb)]
    uprns_undersameroof = undersameroof.((stockdata.UPRN,),(stockdata.osgb,),stockdata.osgb) # <- WARNING: SLOW OPERATION

    # Exclude UPRNs with no OSGB link (failed geo-ref)
    i_nolink = findall(stockdata.iGml .== 0)

    # Find multi-dwelling buildings only
    n_uprns_undersameroof = length.(uprns_undersameroof)
    n_uprns_undersameroof[i_nolink] .= 0
    i_flats = findall(n_uprns_undersameroof.>1)

    # Create unique_id from combinations of all UPRNs (including self ref)
    unique_id = join.(uprns_undersameroof[i_flats],"+")

    # Create a template missing/string array
    uprns_undersameroof = fill!(Array{String}(undef,nrow(stockdata)),"")
    uprns_undersameroof[i_flats] .= unique_id
    return n_uprns_undersameroof, uprns_undersameroof
end


text:
"""
`Ratios` provides `SimpleRatio`, a faster variant of `Rational`.
Speed comes at the cost of greater vulnerability to overflow.

API summary:

- `r = SimpleRatio(num, den)` is equivalent to `num // den`
- `common_denominator` standardizes a collection of `SimpleRatio`s to have the same denominator,
  making some arithmetic operations among them less likely to overflow.
"""
module Ratios

import Base: convert, promote_rule, *, /, +, -, ^, ==, decompose

using Requires

export SimpleRatio, common_denominator

struct SimpleRatio{T<:Integer} <: Real
    num::T
    den::T
end

"""
    SimpleRatio(num::Integer, den::Integer)

Construct the equivalent of the rational number `num // den`.

Operations with `SimpleRatio` are faster, but also more vulnerable to integer overflow,
than with `Rational`. Arithmetic with `SimpleRatio` does not perform any simplification of the resulting ratios.
The best defense against overflow is to use ratios with the same denominator: in such cases,
`+` and `-` will skip forming the product of denominators. See [`common_denominator`](@ref).

If overflow is a risk, consider constructing them using SaferIntegers.jl.

# Examples

```julia
julia> x, y, z = SimpleRatio(1, 8), SimpleRatio(1, 4), SimpleRatio(2, 8)
(SimpleRatio{$Int}(1, 8), SimpleRatio{$Int}(1, 4), SimpleRatio{$Int}(2, 8))

julia> x+y
SimpleRatio{$Int}(12, 32)

julia> x+z
SimpleRatio{$Int}(3, 8)
```
"""
SimpleRatio(num::Integer, den::Integer) = SimpleRatio(promote(num, den)...)

convert(::Type{BigFloat}, r::SimpleRatio{S}) where {S} = BigFloat(r.num)/r.den
function convert(::Type{T}, r::SimpleRatio{S}) where {T<:AbstractFloat,S}
    P = promote_type(T,S)
    convert(T, convert(P, r.num)/convert(P, r.den))
end
SimpleRatio{T}(i::Integer) where {T<:Integer} = SimpleRatio{T}(convert(T, i), oneunit(T))
SimpleRatio{T}(r::Rational{S}) where {T<:Integer, S<:Integer} = SimpleRatio(convert(T, r.num), convert(T, r.den))
Rational{T}(r::SimpleRatio{S}) where {T<:Integer, S<:Integer} = convert(T, r.num) // convert(T, r.den)

*(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.num, x.den*y.den)
*(x::SimpleRatio, y::Bool) = SimpleRatio(x.num*y, x.den)
*(x::SimpleRatio, y::Integer) = SimpleRatio(x.num*y, x.den)
*(x::Bool, y::SimpleRatio) = SimpleRatio(x*y.num, y.den)
*(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.num, y.den)
/(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den, x.den*y.num)
/(x::SimpleRatio, y::Integer) = SimpleRatio(x.num, x.den*y)
/(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den, y.num)
+(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den + y.num, y.den)
-(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den - y.num, y.den)
+(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num + y.num, x.den) :
                                                     SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)
-(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num - y.num, x.den) :
                                                     SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)
^(x::SimpleRatio, y::Integer) = SimpleRatio(x.num^y, x.den^y)

-(x::SimpleRatio{T}) where {T<:Signed} = SimpleRatio(-x.num, x.den)
-(x::SimpleRatio{T}) where {T<:Unsigned} = throw(VERSION < v"0.7.0-DEV.1269" ? OverflowError() : OverflowError("cannot negate unsigned number"))

promote_rule(::Type{SimpleRatio{T}}, ::Type{S}) where {T<:Integer,S<:Integer} = SimpleRatio{promote_type(T,S)}
promote_rule(::Type{SimpleRatio{T}}, ::Type{SimpleRatio{S}}) where {T<:Integer,S<:Integer} = SimpleRatio{promote_type(T,S)}
promote_rule(::Type{SimpleRatio{T}}, ::Type{S}) where {T<:Integer,S<:AbstractFloat} = promote_type(T,S)
promote_rule(::Type{SimpleRatio{T}}, ::Type{Rational{S}}) where {T<:Integer,S<:Integer} = Rational{promote_type(T,S)}

==(x::SimpleRatio, y::SimpleRatio) = x.num*y.den == x.den*y.num

==(x::SimpleRatio, y::Integer) = x.num == x.den*y
==(x::Integer, y::SimpleRatio) = x*y.den == y.num

function ==(x::AbstractFloat, q::SimpleRatio)
    if isfinite(x)
        (count_ones(q.den) == 1) & (x*q.den == q.num)
    else
        x == q.num/q.den
    end
end

==(q::SimpleRatio, x::AbstractFloat) = x == q

decompose(x::SimpleRatio) = x.num, 0, x.den

"""
    common_denominator(x::SimpleRatio, ys::SimpleRatio...)

Return the equivalent of `(x, ys...)` but using a common denominator.
This can be useful to avoid overflow.

!!! info
    This function is quite slow.  In performance-sensitive contexts where the
    ratios are constructed with literal constants, it is better to ensure a common
    denominator at the time of original construction.
"""
function common_denominator(x::SimpleRatio, ys::SimpleRatio...)
    all(y -> y.den == x.den, ys) && return (x, ys...)
    cd = gcd(x.den, map(y -> y.den, ys)...)        # common divisor
    # product of the denominators
    pd = Base.Checked.checked_mul(cd, mapreduce(z -> z.den ÷ cd, Base.Checked.checked_mul, (x, ys...)))
    return map((x, ys...)) do z
        SimpleRatio(z.num * (pd ÷ z.den), pd)
    end
end

function __init__()
    @require FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" begin
        using .FixedPointNumbers: FixedPoint, Fixed, Normed, rawone
        rawone_noerr(::Type{Fixed{T,f}}) where {T,f} = widen(oneunit(T)) << f
        rawone_noerr(::Type{N}) where N<:Normed = rawone(N)
        rawone_noerr(x::FixedPoint) = rawone_noerr(typeof(x))
        Base.promote_rule(::Type{SimpleRatio{S}}, ::Type{<:FixedPoint{T}}) where {S<:Integer,T<:Integer} = SimpleRatio{promote_type(S, T)}
        SimpleRatio{S}(x::FixedPoint) where S<:Integer = SimpleRatio{S}(reinterpret(x), rawone_noerr(x))
        SimpleRatio(x::FixedPoint) = SimpleRatio(reinterpret(x), rawone_noerr(x))
        Base.convert(::Type{S}, x::FixedPoint) where S<:SimpleRatio = S(x)
    end
end

end


text:
using HOHQMesh
using Test
using AbaqusReader

# Start with a clean environment: remove HOHQMesh output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "HOHQMesh.jl" begin

  @testset "examples_dir()" begin
    @test occursin("examples", HOHQMesh.examples_dir())
  end

  @testset "generate_mesh()" begin
    control_file = joinpath(HOHQMesh.examples_dir(), "GingerbreadMan.control")
    @test generate_mesh(control_file) isa String
  end

  @testset "generate_mesh(; verbose=true)" begin
    control_file = joinpath(HOHQMesh.examples_dir(), "GingerbreadMan.control")
    @test generate_mesh(control_file, verbose=true) isa String
  end

  @testset "generate_mesh() with ABAQUS output" begin
    control_file = joinpath(HOHQMesh.examples_dir(), "IceCreamCone_Abaqus.control")
    generate_mesh(control_file)
    parse_mesh = abaqus_read_mesh(joinpath(outdir, "IceCreamCone_Abaqus.inp"))
    # set some reference values for comparison. These are the corner IDs for element 114
    reference_ids = [140, 141, 154, 153]
    @test parse_mesh["elements"][114] == reference_ids
  end

  @testset "generate_mesh() with invalid format" begin
    # Create temporary control file option that is invalid
    mktemp() do path, io
      write(io, "mesh file format = ABBAKISS")
      flush(io)
      @test_throws ErrorException generate_mesh(path)
    end
  end

end # testset "HOHQMesh.jl"

# Clean up afterwards: delete HOHQMesh output directory
@test_nowarn rm(outdir, recursive=true)


text:
using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots; gr()
H(q,p) = norm(p)^2/2 - inv(norm(q))
L(q,p) = q[1]*p[2] - p[1]*q[2]

pdot(dp,p,q,params,t) = ForwardDiff.gradient!(dp, q->-H(q, p), q)
qdot(dq,p,q,params,t) = ForwardDiff.gradient!(dq, p-> H(q, p), p)

initial_position = [.4, 0]
initial_velocity = [0., 2.]
initial_cond = (initial_position, initial_velocity)
initial_first_integrals = (H(initial_cond...), L(initial_cond...))
tspan = (0,20.)
prob = DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)
sol = solve(prob, KahanLi6(), dt=1//10);

plot_orbit(sol) = plot(sol,vars=(3,4), lab="Orbit", title="Kepler Problem Solution")

function plot_first_integrals(sol, H, L)
    plot(initial_first_integrals[1].-map(u->H(u[2,:], u[1,:]), sol.u), lab="Energy variation", title="First Integrals")
    plot!(initial_first_integrals[2].-map(u->L(u[2,:], u[1,:]), sol.u), lab="Angular momentum variation")
end
analysis_plot(sol, H, L) = plot(plot_orbit(sol), plot_first_integrals(sol, H, L))

analysis_plot(sol, H, L)

sol2 = solve(prob, DPRKN6())  # dt is not necessary, because unlike symplectic
                              # integrators DPRKN6 is adaptive
@show sol2.u |> length
analysis_plot(sol2, H, L)

sol3 = solve(prob, ERKN4()) # dt is not necessary, because unlike symplectic
                            # integrators ERKN4 is adaptive
@show sol3.u |> length
analysis_plot(sol3, H, L)

sol4 = solve(prob, Tsit5())
@show sol4.u |> length
analysis_plot(sol4, H, L)

using DiffEqCallbacks

plot_orbit2(sol) = plot(sol,vars=(1,2), lab="Orbit", title="Kepler Problem Solution")

function plot_first_integrals2(sol, H, L)
    plot(initial_first_integrals[1].-map(u->H(u[1:2],u[3:4]), sol.u), lab="Energy variation", title="First Integrals")
    plot!(initial_first_integrals[2].-map(u->L(u[1:2],u[3:4]), sol.u), lab="Angular momentum variation")
end

analysis_plot2(sol, H, L) = plot(plot_orbit2(sol), plot_first_integrals2(sol, H, L))

function hamiltonian(du,u,params,t)
    q, p = u[1:2], u[3:4]
    qdot(@view(du[1:2]), p, q, params, t)
    pdot(@view(du[3:4]), p, q, params, t)
end

prob2 = ODEProblem(hamiltonian, [initial_position; initial_velocity], tspan)
sol_ = solve(prob2, RK4(), dt=1//5, adaptive=false)
analysis_plot2(sol_, H, L)

function first_integrals_manifold(residual,u)
    residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])
    residual[3:4] .= initial_first_integrals[2] - L(u[1:2], u[3:4])
end

cb = ManifoldProjection(first_integrals_manifold)
sol5 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=cb)
analysis_plot2(sol5, H, L)

function energy_manifold(residual,u)
    residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])
    residual[3:4] .= 0
end
energy_cb = ManifoldProjection(energy_manifold)
sol6 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=energy_cb)
analysis_plot2(sol6, H, L)

function angular_manifold(residual,u)
    residual[1:2] .= initial_first_integrals[2] - L(u[1:2], u[3:4])
    residual[3:4] .= 0
end
angular_cb = ManifoldProjection(angular_manifold)
sol7 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=angular_cb)
analysis_plot2(sol7, H, L)

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])