text:
#=
TokenizeTests:
- Julia version: 1.0
- Author: singularitti
- Date: 2019-08-02
=#
module TokenizeTests

using Test

using Fortran90Namelists.Tokenize

@testset "Test bcast" begin
    benchmark = [["&", "bcast_nml"],
    ["    ", "x", " ", "=", " ", "2", "*", "2.0"],
    ["    ", "y", " ", "=", " ", "3", "*"],
    ["    ", "z", " ", "=", " ", "4", "*", ".true."],
    ["/"],
    [],
    ["&", "bcast_endnull_nml"],
    ["    ", "x", " ", "=", " ", "2", "*", "2.0"],
    ["    ", "y", " ", "=", " ", "3", "*"],
    ["/"],
    [],
    ["&", "bcast_mixed_nml"],
    ["    ", "x", " ", "=", " ", "3", "*", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "y", " ", "=", " ", "3", "*", "1", ",", " ", "2", "*", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/bcast.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test bcast target" begin
    benchmark = [["&", "bcast_nml"],
    ["    ", "x", " ", "=", " ", "2.0", ",", " ", "2.0"],
    ["    ", "y", " ", "=", " ", ",", " ", ",", " ", ","],
    ["    ", "z", " ", "=", " ", ".true.", ",", " ", ".true.", ",", " ", ".true.", ",", " ", ".true."],
    ["/"],
    [],
    ["&", "bcast_endnull_nml"],
    ["    ", "x", " ", "=", " ", "2.0", ",", " ", "2.0"],
    ["    ", "y", " ", "=", " ", ",", " ", ",", " ", ","],
    ["/"],
    [],
    ["&", "bcast_mixed_nml"],
    ["    ", "x", " ", "=", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "y", " ", "=", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "2", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/bcast_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test string" begin
    benchmark = [["&", "string_nml"],
    ["    ", "str_basic", " ", "=", " ", "'hello'"],
    ["    ", "str_no_delim", " ", "=", " ", "hello"],
    ["    ", "str_no_delim_no_esc", " ", "=", " ", "a''b"],
    ["    ", "single_esc_delim", " ", "=", " ", "'a ''single'' delimiter'"],
    ["    ", "double_esc_delim", " ", "=", " ", "\"a \"\"double\"\" delimiter\""],
    ["    ", "double_nested", " ", "=", " ", "\"''x'' \"\"y\"\"\""],
    ["    ", "str_list", " ", "=", " ", "'a'", ",", " ", "'b'", ",", " ", "'c'"],
    ["    ", "slist_no_space", " ", "=", " ", "'a'", ",", "'b'", ",", "'c'"],
    ["    ", "slist_no_quote", " ", "=", " ", "a", ",", "b", ",", "c"],
    ["    ", "slash", " ", "=", " ", "'back\\slash'"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/string.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test dollar" begin
    benchmark = [[raw"$", "dollar_nml"],
    ["    ", "v", " ", "=", " ", "1.00"],
    [raw"$"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/dollar.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test dollar target" begin
    benchmark = [["&", "dollar_nml"],
    ["    ", "v", " ", "=", " ", "1.0"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/dollar_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test comment argument" begin
    benchmark = [["&", "comment_alt_nml"],
    ["    ", "x", " ", "=", " ", "1"],
    ["    ", "#y = 2"],
    ["    ", "z", " ", "=", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/comment_alt.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test comment patch" begin
    benchmark = [["! This is an external comment"],
    ["&", "comment_nml"],
    ["    ", "v_cmt_inline", " ", "=", " ", "456", "  ", "! This is an inline comment"],
    ["    ", "! This is a separate comment"],
    ["    ", "v_cmt_in_str", " ", "=", " ", "'This token ! is not a comment'"],
    ["    ", "v_cmt_after_str", " ", "=", " ", "'This ! is not a comment'", " ", "! But this is"],
    ["/"],
    ["! This is a post-namelist comment"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/comment_patch.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test comment target" begin
    benchmark = [["&", "comment_nml"],
    ["    ", "v_cmt_inline", " ", "=", " ", "123"],
    ["    ", "v_cmt_in_str", " ", "=", " ", "'This token ! is not a comment'"],
    ["    ", "v_cmt_after_str", " ", "=", " ", "'This ! is not a comment'"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/comment_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test comment" begin
    benchmark = [["! This is an external comment"],
    ["&", "comment_nml"],
    ["    ", "v_cmt_inline", " ", "=", " ", "123", "  ", "! This is an inline comment"],
    ["    ", "! This is a separate comment"],
    ["    ", "v_cmt_in_str", " ", "=", " ", "'This token ! is not a comment'"],
    ["    ", "v_cmt_after_str", " ", "=", " ", "'This ! is not a comment'", " ", "! But this is"],
    ["/"],
    ["! This is a post-namelist comment"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/comment.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test default index" begin
    benchmark = [["&", "default_index_nml"],
    ["    ", "v", "(", "3", ":", "5", ")", " ", "=", " ", "3", ",", " ", "4", ",", " ", "5"],
    ["    ", "v", " ", "=", " ", "1", ",", " ", "2"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/default_index.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test empty namelist" begin
    benchmark = [["&", "empty_nml"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/empty.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test external token" begin
    benchmark = [["a"],
    ["123"],
    ["&", "ext_token_nml"],
    ["    ", "x", " ", "=", " ", "1"],
    ["/"],
    ["456"],
    ["z"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/ext_token.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

# FIXME:
# @testset "Test external comment" begin
#     benchmark = [[],
#     ["&", "efitin"],
#     ["abc", " ", "=", " ", "0"],
#     ["/"]]
#     tk = Tokenizer()
#     open(joinpath(dirname(@__FILE__), "data/extern_cmt.nml"), "r") do io
#         for (i, line) in enumerate(eachline(io))
#             @test lex(tk, line) == benchmark[i]
#         end
#     end
# end # testset

@testset "Test float" begin
    benchmark = [["&", "float_nml"],
    ["    ", "v_float", " ", "=", " ", "1.0"],
    ["    ", "v_decimal_start", " ", "=", " ", ".1"],
    ["    ", "v_decimal_end", " ", "=", " ", "1."],
    ["    ", "v_negative", " ", "=", " ", "-1."],
    [],
    ["    ", "v_single", " ", "=", " ", "1.0e0"],
    ["    ", "v_double", " ", "=", " ", "1.0d0"],
    [],
    ["    ", "v_single_upper", " ", "=", " ", "1.0E0"],
    ["    ", "v_double_upper", " ", "=", " ", "1.0D0"],
    [],
    ["    ", "v_positive_index", " ", "=", " ", "1.0e+01"],
    ["    ", "v_negative_index", " ", "=", " ", "1.0e-01"],
    [],
    ["    ", "v_no_exp_pos", " ", "=", " ", "1+0"],
    ["    ", "v_no_exp_neg", " ", "=", " ", "1-0"],
    ["    ", "v_no_exp_pos_dot", " ", "=", " ", "1.+0"],
    ["    ", "v_no_exp_neg_dot", " ", "=", " ", "1.-0"],
    ["    ", "v_neg_no_exp_pos", " ", "=", " ", "-1+0"],
    ["    ", "v_neg_no_exp_neg", " ", "=", " ", "-1-0"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/float.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test float target" begin
    benchmark = [["&", "float_nml"],
    ["    ", "v_float", " ", "=", " ", "1.0"],
    ["    ", "v_decimal_start", " ", "=", " ", "0.1"],
    ["    ", "v_decimal_end", " ", "=", " ", "1.0"],
    ["    ", "v_negative", " ", "=", " ", "-1.0"],
    ["    ", "v_single", " ", "=", " ", "1.0"],
    ["    ", "v_double", " ", "=", " ", "1.0"],
    ["    ", "v_single_upper", " ", "=", " ", "1.0"],
    ["    ", "v_double_upper", " ", "=", " ", "1.0"],
    ["    ", "v_positive_index", " ", "=", " ", "10.0"],
    ["    ", "v_negative_index", " ", "=", " ", "0.1"],
    ["    ", "v_no_exp_pos", " ", "=", " ", "1.0"],
    ["    ", "v_no_exp_neg", " ", "=", " ", "1.0"],
    ["    ", "v_no_exp_pos_dot", " ", "=", " ", "1.0"],
    ["    ", "v_no_exp_neg_dot", " ", "=", " ", "1.0"],
    ["    ", "v_neg_no_exp_pos", " ", "=", " ", "-1.0"],
    ["    ", "v_neg_no_exp_neg", " ", "=", " ", "-1.0"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/float_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test float format" begin
    benchmark = [["&", "float_nml"],
    ["    ", "v_float", " ", "=", " ", "1.000"],
    ["    ", "v_decimal_start", " ", "=", " ", "0.100"],
    ["    ", "v_decimal_end", " ", "=", " ", "1.000"],
    ["    ", "v_negative", " ", "=", " ", "-1.000"],
    ["    ", "v_single", " ", "=", " ", "1.000"],
    ["    ", "v_double", " ", "=", " ", "1.000"],
    ["    ", "v_single_upper", " ", "=", " ", "1.000"],
    ["    ", "v_double_upper", " ", "=", " ", "1.000"],
    ["    ", "v_positive_index", " ", "=", " ", "10.000"],
    ["    ", "v_negative_index", " ", "=", " ", "0.100"],
    ["    ", "v_no_exp_pos", " ", "=", " ", "1.000"],
    ["    ", "v_no_exp_neg", " ", "=", " ", "1.000"],
    ["    ", "v_no_exp_pos_dot", " ", "=", " ", "1.000"],
    ["    ", "v_no_exp_neg_dot", " ", "=", " ", "1.000"],
    ["    ", "v_neg_no_exp_pos", " ", "=", " ", "-1.000"],
    ["    ", "v_neg_no_exp_neg", " ", "=", " ", "-1.000"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/float_format.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test global index" begin
    benchmark = [["&", "global_index_nml"],
    ["    ", "v_zero", "(", "0", ":", "3", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "v_neg", "(", "-2", ":", "1", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "v_pos", "(", "2", ":", "5", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/global_index.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test group repeat" begin
    benchmark = [["&", "grp_repeat_nml"],
    ["    ", "x", " ", "=", " ", "1"],
    ["/"],
    ["&", "grp_repeat_nml"],
    ["    ", "x", " ", "=", " ", "2"],
    ["/"],
    ["&", "CASE_CHECK_nml"],
    ["    ", "y", " ", "=", " ", "1"],
    ["/"],
    ["&", "CASE_CHECK_nml"],
    ["    ", "y", " ", "=", " ", "2"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/grp_repeat.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test group repeat target" begin
    benchmark = [["&", "grp_repeat_nml"],
    ["    ", "x", " ", "=", " ", "1"],
    ["/"],
    [],
    ["&", "grp_repeat_nml"],
    ["    ", "x", " ", "=", " ", "2"],
    ["/"],
    [],
    ["&", "case_check_nml"],
    ["    ", "y", " ", "=", " ", "1"],
    ["/"],
    [],
    ["&", "case_check_nml"],
    ["    ", "y", " ", "=", " ", "2"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/grp_repeat_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index bad end" begin
    benchmark = [["&", "bad_index_nml"],
    ["    ", "y", "(", "1", ":", "~", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_bad_end.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index bad start" begin
    benchmark = [["&", "bad_index_nml"],
    ["    ", "y", "(", "1", "~", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_bad_start.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index bad stride" begin
    benchmark = [["&", "bad_index_nml"],
    ["    ", "y", "(", "1", ":", "3", ":", "~", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_bad_stride.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index bad" begin
    benchmark = [["&", "bad_index_nml"],
    ["    ", "y", "(", "~", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_bad.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index empty end" begin
    benchmark = [["&", "empty_index_nml"],
    ["    ", "x", "(", "1", ":", ":", "2", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_empty_end.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index empty stride" begin
    benchmark = [["&", "empty_index_nml"],
    ["    ", "x", "(", "1", ":", "3", ":", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_empty_stride.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index empty" begin
    benchmark = [["&", "empty_index_nml"],
    ["    ", "x", "(", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_empty.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test index zero stride" begin
    benchmark = [["&", "bad_index_nml"],
    ["    ", "y", "(", "1", ":", "3", ":", "0", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/index_zero_stride.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test logical representation" begin
    benchmark = [["&", "logical_nml"],
    ["    ", "a", " ", "=", " ", "T"],
    ["    ", "b", " ", "=", " ", "F"],
    ["    ", "c", " ", "=", " ", "T"],
    ["    ", "d", " ", "=", " ", "F"],
    ["    ", "e", " ", "=", " ", "T"],
    ["    ", "f", " ", "=", " ", "F"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/logical_repr.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test logical" begin
    benchmark = [["&", "logical_nml"],
    ["    ", "a", " ", "=", " ", ".true."],
    ["    ", "b", " ", "=", " ", ".false."],
    ["    ", "c", " ", "=", " ", "t"],
    ["    ", "d", " ", "=", " ", "f"],
    ["    ", "e", " ", "=", " ", ".t"],
    ["    ", "f", " ", "=", " ", ".f"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/logical.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multidimensional target" begin
    benchmark = [["&", "multidim_ooo_nml"],
    ["    ", "a", "(", "1", ",", "1", ")", " ", "=", " ", "1"],
    ["    ", "a", "(", "1", ":", "2", ",", "2", ")", " ", "=", " ", ",", " ", "2"],
    ["    ", "b", "(", "1", ",", "1", ")", " ", "=", " ", "1"],
    ["    ", "b", "(", "1", ":", "3", ",", "2", ")", " ", "=", " ", ",", " ", ",", " ", "3"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multidim_ooo_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multidimensional" begin
    benchmark = [["&", "multidim_ooo_nml"],
    ["    ", "a", "(", "2", ",", "2", ")", " ", "=", " ", "2"],
    ["    ", "a", "(", "1", ",", "1", ")", " ", "=", " ", "1"],
    [],
    ["    ", "b", "(", "3", ",", "2", ")", " ", "=", " ", "3"],
    ["    ", "b", "(", "1", ",", "1", ")", " ", "=", " ", "1"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multidim_ooo.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multidimensional space" begin
    benchmark = [["&", "multidim_nml"],
    ["    ", "v2d", "(", "1", ":", "2", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d", "(", "1", ":", "2", ",", " ", "2", ")", " ", "=", " ", "3", ",", " ", "4"],
    ["    ", "v3d", "(", "1", ":", "2", ",", " ", "1", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v3d", "(", "1", ":", "2", ",", " ", "2", ",", " ", "1", ")", " ", "=", " ", "3", ",", " ", "4"],
    ["    ", "v3d", "(", "1", ":", "2", ",", " ", "1", ",", " ", "2", ")", " ", "=", " ", "5", ",", " ", "6"],
    ["    ", "v3d", "(", "1", ":", "2", ",", " ", "2", ",", " ", "2", ")", " ", "=", " ", "7", ",", " ", "8"],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "1", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "2", ",", " ", "1", ")", " ", "=", " ", "5", ",", " ", "6", ",", " ", "7", ",", " ", "8"],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "3", ",", " ", "1", ")", " ", "=", " ", "9", ",", " ", "10", ",", " ", "11", ",", " ", "12"],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "1", ",", " ", "2", ")", " ", "=", " ", "13", ",", " ", "14", ",", " ", "15", ",", " ", "16"],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "2", ",", " ", "2", ")", " ", "=", " ", "17", ",", " ", "18", ",", " ", "19", ",", " ", "20"],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "3", ",", " ", "2", ")", " ", "=", " ", "21", ",", " ", "22", ",", " ", "23", ",", " ", "24"],
    ["    ", "v2d_explicit", "(", "1", ":", "2", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d_explicit", "(", "1", ":", "2", ",", " ", "2", ")", " ", "=", " ", "3", ",", " ", "4"],
    ["    ", "v2d_outer", "(", "1", ",", " ", "1", ")", " ", "=", " ", "1"],
    ["    ", "v2d_outer", "(", "1", ",", " ", "2", ")", " ", "=", " ", "2"],
    ["    ", "v2d_outer", "(", "1", ",", " ", "3", ")", " ", "=", " ", "3"],
    ["    ", "v2d_outer", "(", "1", ",", " ", "4", ")", " ", "=", " ", "4"],
    ["    ", "v2d_inner", "(", ":", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "v2d_sparse", "(", ":", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d_sparse", "(", ":", ",", " ", "2", ")", " ", "=", " ", ",", " ", ","],
    ["    ", "v2d_sparse", "(", ":", ",", " ", "3", ")", " ", "=", " ", "5", ",", " ", "6"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multidim_space.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multidimensional target" begin
    benchmark = [["&", "multidim_nml"],
    ["    ", "v2d", "(", "1", ":", "2", ",", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d", "(", "1", ":", "2", ",", "2", ")", " ", "=", " ", "3", ",", " ", "4"],
    ["    ", "v3d", "(", "1", ":", "2", ",", "1", ",", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v3d", "(", "1", ":", "2", ",", "2", ",", "1", ")", " ", "=", " ", "3", ",", " ", "4"],
    ["    ", "v3d", "(", "1", ":", "2", ",", "1", ",", "2", ")", " ", "=", " ", "5", ",", " ", "6"],
    ["    ", "v3d", "(", "1", ":", "2", ",", "2", ",", "2", ")", " ", "=", " ", "7", ",", " ", "8"],
    ["    ", "w3d", "(", "1", ":", "4", ",", "1", ",", "1", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "w3d", "(", "1", ":", "4", ",", "2", ",", "1", ")", " ", "=", " ", "5", ",", " ", "6", ",", " ", "7", ",", " ", "8"],
    ["    ", "w3d", "(", "1", ":", "4", ",", "3", ",", "1", ")", " ", "=", " ", "9", ",", " ", "10", ",", " ", "11", ",", " ", "12"],
    ["    ", "w3d", "(", "1", ":", "4", ",", "1", ",", "2", ")", " ", "=", " ", "13", ",", " ", "14", ",", " ", "15", ",", " ", "16"],
    ["    ", "w3d", "(", "1", ":", "4", ",", "2", ",", "2", ")", " ", "=", " ", "17", ",", " ", "18", ",", " ", "19", ",", " ", "20"],
    ["    ", "w3d", "(", "1", ":", "4", ",", "3", ",", "2", ")", " ", "=", " ", "21", ",", " ", "22", ",", " ", "23", ",", " ", "24"],
    ["    ", "v2d_explicit", "(", "1", ":", "2", ",", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d_explicit", "(", "1", ":", "2", ",", "2", ")", " ", "=", " ", "3", ",", " ", "4"],
    ["    ", "v2d_outer", "(", "1", ",", "1", ")", " ", "=", " ", "1"],
    ["    ", "v2d_outer", "(", "1", ",", "2", ")", " ", "=", " ", "2"],
    ["    ", "v2d_outer", "(", "1", ",", "3", ")", " ", "=", " ", "3"],
    ["    ", "v2d_outer", "(", "1", ",", "4", ")", " ", "=", " ", "4"],
    ["    ", "v2d_inner", "(", ":", ",", "1", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "v2d_sparse", "(", ":", ",", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d_sparse", "(", ":", ",", "2", ")", " ", "=", " ", ",", " ", ","],
    ["    ", "v2d_sparse", "(", ":", ",", "3", ")", " ", "=", " ", "5", ",", " ", "6"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multidim_target.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multidimensional" begin
    benchmark = [["&", "multidim_nml"],
    ["    ", "v2d", "(", "1", ":", "2", ",", " ", "1", ":", "2", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    [],
    ["    ", "v3d", "(", "1", ":", "2", ",", " ", "1", ":", "2", ",", " ", "1", ":", "2", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4", ",", " ", "5", ",", " ", "6", ",", " ", "7", ",", " ", "8"],
    [],
    ["    ", "w3d", "(", "1", ":", "4", ",", " ", "1", ":", "3", ",", " ", "1", ":", "2", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4", ",", " ", "5", ",", " ", "6", ",", " ", "7", ",", " ", "8", ",", " ", "9", ",", " ", "10", ",", " ", "11", ",", " ", "12", ","],
    ["                         ", "13", ",", " ", "14", ",", " ", "15", ",", " ", "16", ",", " ", "17", ",", " ", "18", ",", " ", "19", ",", " ", "20", ",", " ", "21", ",", " ", "22", ",", " ", "23", ",", " ", "24"],
    [],
    ["    ", "v2d_explicit", "(", "1", ",", " ", "1", ")", " ", "=", " ", "1"],
    ["    ", "v2d_explicit", "(", "2", ",", " ", "1", ")", " ", "=", " ", "2"],
    ["    ", "v2d_explicit", "(", "1", ",", " ", "2", ")", " ", "=", " ", "3"],
    ["    ", "v2d_explicit", "(", "2", ",", " ", "2", ")", " ", "=", " ", "4"],
    [],
    ["    ", "v2d_outer", "(", "1", ",", " ", ":", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    ["    ", "v2d_inner", "(", ":", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2", ",", " ", "3", ",", " ", "4"],
    [],
    ["    ", "v2d_sparse", "(", ":", ",", " ", "1", ")", " ", "=", " ", "1", ",", " ", "2"],
    ["    ", "v2d_sparse", "(", ":", ",", " ", "3", ")", " ", "=", " ", "5", ",", " ", "6"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multidim.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multidimensional column width" begin
    benchmark = [["&", "multiline_nml"],
    ["    ", "x", " ", "=", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["        ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["        ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["        ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["        ", "1", ",", " ", "1", ",", " ", "1"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multiline_colwidth.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multiline index" begin
    benchmark = [["&", "multiline_nml"],
    ["    ", "x", "(", "1", ":", "47", ")", " ", "=", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["              ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["              ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multiline_index.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

@testset "Test multiline" begin
    benchmark = [["&", "multiline_nml"],
    ["    ", "x", " ", "=", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["        ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ",", " ", "1", ","],
    ["        ", "1", ",", " ", "1", ",", " ", "1"],
    ["/"]]
    tk = Tokenizer()
    open(joinpath(dirname(@__FILE__), "data/multiline.nml"), "r") do io
        for (i, line) in enumerate(eachline(io))
            @test lex(tk, line) == benchmark[i]
        end
    end
end # testset

end


text:
using Sobol, Compat
using Compat.Test

# compare results with results from C++ code sobol.cc published on
# http://web.maths.unsw.edu.au/~fkuo/sobol/
# with new-joe-kuo-6.21201 file used as input.
#
# Command line used to generate output is
#  ./sobol $N 1024 new-joe-kuo-6.21201 > exp_results_$N
#
# For each of the dimensions below
# (except 21201 where only 64 samples are generated)

dimensions = [5, 10, 20, 50, 100, 500, 21201]

for dim in dimensions
    println("Testing dimension $(dim)")
    open(joinpath(dirname(@__FILE__), "results", "exp_results_$(dim)")) do exp_results_file
        s = SobolSeq(dim)
        x = zeros(dim)
        for line in eachline(exp_results_file)
            values = [parse(Float64, item) for item in split(line)]
            if length(values) > 0
                @test x == values
                next!(s, x)
            end
        end
    end
end

# issue #8
using Base.Iterators: take
@test [x[1] for x in collect(take(Sobol.SobolSeq(1),5))] == [0.5,0.75,0.25,0.375,0.875]

ext: jl


text:
export isbase64

function isbase64(str::AbstractString, urlsafe::Bool=false)::Bool
    notBase64 = r"[^A-Z0-9+\/=]"i
    urlSafeBase64 = r"^[A-Z0-9_\-]*$"i

    if str === ""
        return true
    end

    if urlsafe
        return contains(str, urlSafeBase64)
    end

    if length(str) % 4 !== 0 || contains(str, notBase64)
        return false
    end

    firstPadding = findfirst("=", str)
    if firstPadding !== nothing
        firstPadding = first(firstPadding)
    end

    return firstPadding === nothing ||
        firstPadding === length(str) ||
        (firstPadding === length(str) - 1 && str[end] === '=')
end


text:
module ConfusionMatrix

using NamedArrays
#using NamedArrays: NamedMatrix
using Flux: OneHotArray, OneHotMatrix, onecold
#import Base.push!
#import Base.summary
#import Base.show

export ConfMat, update!, metrics, metrics!, toidx

include("confmat.jl")
include("metrics.jl")

#=export ConfMatrix

using NamedArrays
using Flux: OneHotMatrix, onecold

struct ConfMatrix{N}
    labs
    cmat
end

ConfMatrix(l, m) = ConfMatrix{length(l)}(l, m)

function ConfMatrix(l)
    n = length(l)
    m = NamedArray(zeros(Int, n, n), (l, l), ("Actual", "Predicted"))
    ConfMatrix(l, m)
end

function ConfMatrix(l, x, y)
    cm = ConfMatrix(l)
    push!(cm, x, y)
end

function ConfMatrix(l, x, y; thrd)
    cm = ConfMatrix(l)
    push!(cm, x, y, thrd = thrd)
end

Base.show(io::IO, cm::ConfMatrix) = show(cm.cmat)

function labidx(v, l)
    i = indexin(v, l)
    all(ismissing.(i)) && throw(ArgumentError(
        "Unknown labels in supplied data."))
    Vector{Int}(i)
end

function Base.push!(cm::ConfMatrix, x, y::Vector{T}; thrd) where T <: AbstractFloat
    x = labidx(x, cm.labs)
    y = @. Int(y > thrd) + 1
    push!(cm, x, y)
end

function Base.push!(cm::ConfMatrix, x::T, y::T) where T <: OneHotMatrix
    x = onecold(x)
    y = onecold(y)
    push!(cm, x, y)
end

function Base.push!(cm::ConfMatrix, x::Vector{T}, y::Vector{T}) where T
    x = labidx(x, cm.labs)
    y = labidx(y, cm.labs)
    push!(cm, x, y)
end

function Base.push!(cm::ConfMatrix, x::Vector{Int}, y::Vector{Int})
    length(y) == length(x) || throw(DimensionMismatch(
        "Number of actual values and number of predicted values mismatch."))
    for i in eachindex(y)
        @inbounds r = x[i]
        @inbounds c = y[i]
        cm.cmat[r, c] += 1
    end
    cm
end
=#
end


text:
## User should just use TransferFunction
struct SisoRational{T} <: SisoTf{T}
    num::Poly{T}
    den::Poly{T}
    function SisoRational{T}(num::Poly{T}, den::Poly{T}) where T <: Number
        if all(den == zero(den))
            error("Cannot create SisoRational with zero denominator")
        elseif all(num == zero(num))
            # The numerator is zero, make the denominator 1
            den = one(den)
        end
        new{T}(num, den)
    end
end
function SisoRational(num::Poly{T1}, den::Poly{T2}) where T1 <: Number where T2 <: Number
    T = promote_type(T1,T2)
    SisoRational{T}(Poly{T}(num.a), Poly{T}(den.a))
end
function SisoRational{T}(num::AbstractVector, den::AbstractVector) where T <: Number # NOTE: Typearguemnts on the parameters?
    SisoRational{T}(Poly{T}(reverse(num)), Poly{T}(reverse(den)))
end
function SisoRational(num::AbstractVector{T1}, den::AbstractVector{T2}) where T1 <: Number where T2 <: Number
    T = promote_type(T1,T2)
    SisoRational{T}(num, den)
end
# NOTE: How many of these above are actually needed?
# TODO: Add method for scalar inputs

Base.zero(::Type{SisoRational{T}}) where T = SisoRational{T}([zero(T)], [one(T)])
Base.one(::Type{SisoRational{T}}) where T = SisoRational{T}([one(T)], [one(T)])

Base.one(f::SisoRational) = one(typeof(f))
Base.zero(f::SisoRational) = zero(typeof(f))

isproper(f::SisoRational) = (length(f.num) <= length(f.den))

function minreal(sys::SisoRational{T}, eps::Real=sqrt(eps())) where T
    T2 = typeof(zero(T)/one(T))
    return convert(SisoRational{T2}, minreal(convert(SisoZpk,sys), eps))
end

function print_siso(io::IO, f::SisoRational, var=:s)
    # Convert the numerator and denominator to strings
    numstr = sprint(printpolyfun(var), f.num)
    denstr = sprint(printpolyfun(var), f.den)

    # Figure out the length of the separating line
    len_num = length(numstr)
    len_den = length(denstr)
    dashcount = max(len_num, len_den)

    # Center the numerator or denominator
    if len_num < dashcount
        numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr"
    else
        denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr"
    end
    println(io, numstr)
    println(io, repeat("-", dashcount))
    println(io, denstr)
end

# TODO is this working?
function print_compact(io::Base.IO, f::SisoRational, var)
    numstr = sprint(print_poly, f.num)
    denstr = sprint(print_poly, f.den)
    println(io, "($numstr)/($denstr)")
end

numvec(f::SisoRational) = reverse(coeffs(f.num))
denvec(f::SisoRational) = reverse(coeffs(f.den))

denpoly(f::SisoRational) = f.den
numpoly(f::SisoRational) = f.num

tzero(f::SisoRational) = roots(f.num)
pole(f::SisoRational) = roots(f.den)

function evalfr(f::SisoRational{T}, s::Number) where T
    S = promote_op(/, promote_type(T, typeof(s)), promote_type(T, typeof(s)))
    den = polyval(f.den, s)
    if den == zero(S)
        convert(S,Inf)
    else
        polyval(f.num, s)/den
    end
end

==(f1::SisoRational, f2::SisoRational) = (f1.num * f2.den[1] == f2.num * f1.den[1] && f1.den * f2.den[1] == f2.den * f1.den[1]) # NOTE: Not in analogy with how it's done for SisoZpk

# We might want to consider alowing scaled num and den as equal
function isapprox(f1::SisoRational, f2::SisoRational; rtol::Real=sqrt(eps()), atol::Real=0)
    isapprox(f1.num * f2.den[1], f2.num * f1.den[1], rtol=rtol, atol=atol) && isapprox(f1.den * f2.den[1], f2.den * f1.den[1], rtol=rtol, atol=atol)
end

+(f1::SisoRational, f2::SisoRational) = SisoRational(f1.num*f2.den + f2.num*f1.den, f1.den*f2.den)
+(f::SisoRational, n::Number) = SisoRational(f.num + n*f.den, f.den)
+(n::Number, f::SisoRational) = f + n
#.+(f::SisoRational, n::Number) = t + n
#.+(n::Number, f::SisoRational) = t + n

-(f1::SisoRational, f2::SisoRational) = SisoRational(f1.num*f2.den - f2.num*f1.den, f1.den*f2.den)
-(n::Number, f::SisoRational) = SisoRational(n*f.den - f.num, f.den)
-(f::SisoRational, n::Number) = +(f, -n)
#.-(f::SisoRational, n::Number) = -(t, n)
#.-(n::Number, f::SisoRational) = -(n, t)

-(f::SisoRational) = SisoRational(-f.num, f.den)

*(f1::SisoRational, f2::SisoRational) = SisoRational(f1.num*f2.num, f1.den*f2.den)
*(f::SisoRational, n::Number) = SisoRational(f.num*n, f.den)
*(n::Number, f::SisoRational) = *(f, n)
#.*(f1::SisoRational, f2::SisoRational) = *(f1, f2)
#.*(f::SisoRational, n::Number) = *(f, n)
#.*(n::Number, f::SisoRational) = *(f, n)

/(n::Number, f::SisoRational) = SisoRational(n*f.den, f.num)
/(f::SisoRational, n::Number) = f*(1/n)
/(f1::SisoRational, f2::SisoRational) = f1*(1/f2)


text:
Base.include(@__MODULE__, "../../get_indat.jl")
#=
function xyz_to_geometry(xyzfile)
  input_file = []
  #println(joinpath(@__DIR__, xyzfile))
  open(joinpath(@__DIR__, xyzfile)) do file
    input_file = readlines(file)
  end

  natoms = parse(Int, input_file[1])
  #println(natoms)
  geometry = input_file[3:end]
  #display(geometry); println()

  coords_vector = Vector{Float64}([])
  symbols = Vector{String}([])
  for geometry_line in geometry
    coords_match = eachmatch(r"([-]?[\d]{1,}\.[\d]{1,})",geometry_line)
    atom_coords = collect(coords_match)
    atom_coords_match = map(x -> x.match, atom_coords)
    append!(coords_vector, parse.(Float64,atom_coords_match[1:3]))

    symbol_match = match(r"[A-Za-z]{1,2}",geometry_line)
    push!(symbols, symbol_match.match)
  end
  #coords = reshape(coords_vector, (3,length(geometry)))
  #coords = transpose(coords)
  coords = coords_vector

  #@assert natoms == size(coords)[1]
  #@assert natoms == length(symbols)

  #display(coords); println()
  #display(symbols); println()

  charge = 0

  return coords, symbols, charge
end
=#
function create_input_mbe()
  #== read in fragment files ==#
  directory = joinpath(@__DIR__) * "/"
  all_files = readdir(directory)
  all_files .= directory .* all_files

  #input_match = match.(r"frag[0-9]+.xyz",all_files)
  #input_match_filter = filter!(x -> x != nothing, input_match)
  #inputs = map(x->x.match, input_match)

  num_frags = begin
    nfrags = 0
    Nfrag = open("../Nfrags", "r") do file
      argv = readlines(file)
      nfrags = parse(Int64, argv[2])
    end
    nfrags
  end

  fragments = [ "$i" for i in 1:num_frags ]
  inputs = "frag" .* fragments .* ".xyz"
  display(inputs)

  fmobnd::Vector{String} = []
  for frag in 1:num_frags
    bnd = open("../BROKENBONDS", "r") do file
      fmobnd = readlines(file)
      #push!(fmobnd, bonds[2:end])
    end
  end

  coords, symbols, indat = get_indat("../molecule.xyz", num_frags)
  charge = 0

  #== write input json file ==#
  output = open("input.json", "w") do file
    write(file, "{\n")
    #== write in molecule information ==#
    write(file, "  \"molecule\": {\n")
    #== write fragment geometry ==#
    write(file, "    \"geometry\" : [\n")
    for icoord in 1:length(coords)
      write(file, "      ")

      value = coords[icoord]
      if icoord == length(coords)
        write(file, "$value\n")
      elseif icoord%3 == 0
        write(file, "$value,\n")
      else
        write(file, "$value,")
      end
    end
    write(file, "    ],\n")

    #== write fragment symbols ==#
    write(file, "    \"symbols\" : [\n")
    for iatom in 1:length(symbols)
      write(file, "      ")
      symbol = symbols[iatom]
      if iatom == length(symbols)
        write(file, "\"$symbol\"\n")
      elseif iatom%5 == 0
        write(file, "\"$symbol\",\n")
      else
        write(file, "\"$symbol\",")
      end
    end
    write(file, "    ],\n")

    #== assign atoms to fragments ==#
    write(file, "    \"fragments\" : {\n")
    write(file, "      \"nfrag\": $num_frags,\n")
    write(file, "      \"fragid\": [\n")
    for ifrag in 1:length(indat)
      write(file, "      ")
      fragment = indat[ifrag]
      if ifrag == length(indat)
        write(file, "$fragment\n")
      elseif ifrag%5 == 0
        write(file, "$fragment,\n")
      else
        write(file, "$fragment,")
      end
    end
    write(file, "      ],\n")

    #== write broken bonds ==#
    write(file, "      \"broken_bonds\" : [\n")
    for ibnd in fmobnd[2:end]
      bond_match = eachmatch(r"[0-9]+",ibnd)
      broken_bonds = collect(bond_match)
      broken_bonds_match = map(x -> x.match, broken_bonds)
      if ibnd == fmobnd[end]
        write(file,"        $(broken_bonds_match[1]), $(broken_bonds_match[2])\n")
      else
        write(file,"        $(broken_bonds_match[1]), $(broken_bonds_match[2]),\n")
      end
    end
    write(file, "      ],\n")

    #== write fragment charge ==#
    #write(file, "    \"molecular_charge\": $charge\n")
    #write(file, "  },\n")

    #== write fragment charges ==#
    write(file, "      \"fragment_charges\" : [\n")
    for ifrag in 1:num_frags
      write(file, "        ")
      fragment = indat[ifrag]
      if ifrag == num_frags
        write(file, "0\n")
      elseif ifrag%5 == 0
        write(file, "0,\n")
      else
        write(file, "0,")
      end
    end
    write(file, "      ]\n")
    write(file, "    }\n")
    write(file, "  },\n")

    #== write calculation driver and model information ==#
    write(file,"  \"driver\": \"energy\",\n")
    write(file,"  \"model\": {\n")
    write(file,"    \"method\": \"MBE-RIMP2\",\n")
    write(file,"    \"basis\": \"6-31G\",\n")
    write(file,"    \"aux_basis\": \"cc-pVDZ\"\n")
    write(file,"  },\n")

    #== write keywords ==#
    write(file,"  \"keywords\": {\n")
    write(file,"    \"scf\": {\n")
    write(file,"      \"niter\": 200\n")
    write(file,"    },\n")
    write(file,"    \"rimp2\": {\n")
    write(file,"      \"box_dim\": 15\n")
    write(file,"    }\n")
    write(file,"  }\n")
    write(file,"}")
  end
end

function create_input_rhf(input)
  #== write input json file ==#
  filename = input[1:(end-4)]
  #display(filename)
  output = open("$filename.json", "w") do file
    write(file, "{\n")
    #== write in molecule information ==#
    write(file, "  \"molecule\": {\n")
    #== write each fragment's molecule information ==#

    coords, symbols, charge = xyz_to_geometry(input)
    #== write fragment geometry ==#
    write(file, "    \"geometry\" : [\n")
    for icoord in 1:length(coords)
      write(file, "      ")

      value = coords[icoord]
      #println(value)
      if icoord == length(coords)
        write(file, "$value\n")
      elseif icoord%3 == 0
        write(file, "$value,\n")
      else
        write(file, "$value,")
      end
    end
    write(file, "    ],\n")

    #== write fragment symbols ==#
    write(file, "    \"symbols\" : [\n")
    for iatom in 1:length(symbols)
      write(file, "      ")
      symbol = symbols[iatom]
      if iatom == length(symbols)
        write(file, "\"$symbol\"\n")
      elseif iatom%5 == 0
        write(file, "\"$symbol\",\n")
      else
        write(file, "\"$symbol\",")
      end
    end
    write(file, "    ],\n")
    #== write fragment charge ==#
    write(file, "    \"molecular_charge\": $charge\n")
    write(file, "  },\n")

    #== write calculation driver and model information ==#
    write(file,"  \"driver\": \"energy\",\n")
    write(file,"  \"model\": {\n")
    write(file,"    \"method\": \"RHF\",\n")
    write(file,"    \"basis\": \"6-31+G*\"\n")
    write(file,"  },\n")

    #== write keywords ==#
    write(file,"  \"keywords\": {\n")
    write(file,"    \"scf\": {\n")
    write(file,"      \"guess\":\"hcore\",\n")
    write(file,"      \"dele\":1.0, \n")
    write(file,"      \"rmsd\":2.0E-5 \n")
    write(file,"    }\n")
    write(file,"  }\n")
    write(file,"}")
  end
end

create_input_mbe()


text:
export print

###############################################################################
####  Relaxation context
###############################################################################
function print(io::IO, relctx::RelaxationContext)
    print(io, "RelaxationContext:\n")
    print(io, "ismultiordered         : $(relctx.ismultiordered)\n")
    print(io, "issparse               : $(relctx.issparse)\n")
    print(io, "symmetries             : $(relctx.symmetries)\n")
    print(io, "hierarchykind          : $(relctx.hierarchykind)\n")
    print(io, "renamevars             : $(relctx.renamevars)\n")
    for cstrname in sort(collect(keys(relctx.di)))
        di = relctx.di[cstrname]
        print(io, "di                     : $cstrname  \t=> $di\n")
    end
    for cstrname in sort(collect(keys(relctx.ki)))
        ki = relctx.ki[cstrname]
        print(io, "ki                     : $cstrname  \t=> $ki\n")
    end
    for cstrname in sort(collect(keys(relctx.cstrtypes)))
        cstrtype = relctx.cstrtypes[cstrname]
        print(io, "bar var types          : $cstrname  \t=> $(string(cstrtype))\n")
    end
end

###############################################################################
####  SDPDual
###############################################################################
function print(io::IO, momentrelax::SDPDual{T}) where T
    println(io, "Moment Relaxation Problem:")
    println(io, "→ Objective: ")
    momentlen = maximum(x->length(string(x)), keys(momentrelax.objective))
    for moment in sort(collect(keys(momentrelax.objective)))
        coeff = momentrelax.objective[moment]
        print_string(io, string(moment), momentlen)
        println(io, " $coeff")
    end

    println(io, "→ Constraints:")
    for (cstrname, blocname) in sort(collect(keys(momentrelax.constraints)))
        mmtmat = momentrelax.constraints[(cstrname, blocname)]
        println(io, " → $cstrname, $blocname")
        println(io, mmtmat)
    end

    println(io, "→ Moments clique overlap:")
    if length(momentrelax.moments_overlap) > 0
        mmtlength = maximum(x->length(string(x)), keys(momentrelax.moments_overlap))
        for moment in sort(collect(keys(momentrelax.moments_overlap)))
            cliquenames = momentrelax.moments_overlap[moment]
            print(io, " → ")
            print_string(io, string(moment), mmtlength)
            for clique in sort(collect(cliquenames)) print(io, "$clique, ") end
            @printf(io, "\b\b \n")
        end
    else
        print(io, "  None")
    end
end

###############################################################################
####  SDPPrimal
###############################################################################
function print(io::IO, sdpinst::SDPPrimal)
    println(io, " -- SDP Blocks:")
    print(io, sdpinst.blocks)
    println(io, " -- linear part:")
    print(io, sdpinst.lin, sdpinst.linsym)
    println(io, " -- const part:")
    print(io, sdpinst.cst)
    println(io, " -- mat var types:")
    for blockname in sort(collect(keys(sdpinst.block_to_vartype)))
        blocktype = sdpinst.block_to_vartype[blockname]
        println(io, "   $blockname  \t $blocktype")
    end
end

function print(io::IO, sdpblocks::DictType{Tuple{Moment, String, Exponent, Exponent}, T}; indentedprint=true) where T
    print_blocksfile(io, sdpblocks; indentedprint=indentedprint, print_header=false)
end

function print(io::IO, sdplin::DictType{Tuple{Moment, Exponent}, T}, sdplinsym::DictType{Tuple{Moment, String, Exponent}, T}; indentedprint=true) where T
    print_linfile(io, sdplin, sdplinsym; indentedprint=indentedprint, print_header=false)
end

function print(io::IO, sdpcst::DictType{Moment, T}; indentedprint=true) where T
    print_cstfile(io, sdpcst; indentedprint=indentedprint, print_header=false)
end

###############################################################################
####  SDP_Problem
###############################################################################
function print(io::IO, sdp::SDP_Problem)
  for cstr in keys(sdp.name_to_sdpblock)
    block = sdp.name_to_sdpblock[cstr]
    println(io, "  sdp   : $cstr -> $block")
  end

  println(io, "  objk  : $(sdp.obj_keys)")

  for name in keys(sdp.name_to_ctr)
    ctr = sdp.name_to_ctr[name]
    println(io, "  ctr   : $name \t $ctr")
  end

  for name in keys(sdp.matrices)
    mat = sdp.matrices[name]
    println(io, "  matrix: $name \t $mat")
  end

  for name in keys(sdp.linear)
    lin = sdp.linear[name]
    println(io, "  lin   : $name \t $lin")
  end

  for name in keys(sdp.cst_ctr)
    cst = sdp.cst_ctr[name]
    println(io, "  cst   : $name \t $cst")
  end
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using SparseArrays
Random.seed!(233)

function push_matrices(A,B)
    ii1, jj1, vv1 = find(A)
    ii2, jj2, vv2 = find(B)
    d = size(A, 1)
    push_matrices_ = load_op_and_grad("$(@__DIR__)/build/libPushMatrices","push_matrices")
    ii1,jj1,vv1,ii2,jj2,vv2,d = convert_to_tensor([ii1,jj1,vv1,ii2,jj2,vv2,d], [Int64,Int64,Float64,Int64,Int64,Float64,Int64])
    push_matrices_(ii1,jj1,vv1,ii2,jj2,vv2,d)
end

function visco_solve(rhs,vv,op)
    visco_solve_ = load_op_and_grad("$(@__DIR__)/build/libViscoSolve","visco_solve")
    rhs,vv, op = convert_to_tensor([rhs,vv, op], [Float64,Float64, Int64])
    visco_solve_(rhs,vv, op)
end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
export generateGraph_Kaess, generateGraph_TestSymbolic, generateGraph_CaesarRing1D, generateGraph_LineStep
export calcHelix_T
export generateGraph_EuclidDistance

"""
    $SIGNATURES

Canonical example from literature, Kaess, et al.: ISAM2, IJRR, 2011.

Notes
- Paper variable ordering: p = [:l1;:l2;:x1;:x2;:x3]
"""
function generateGraph_Kaess(;graphinit::Bool=false)
  fg = initfg()

  addVariable!(fg,:x1, ContinuousScalar)
  addFactor!(fg, [:x1;], Prior(Normal()), graphinit=graphinit)

  addVariable!(fg,:x2, ContinuousScalar)
  addFactor!(fg,[:x1, :x2], LinearRelative(Normal()), graphinit=graphinit)

  addVariable!(fg, :x3, ContinuousScalar)
  addFactor!(fg,[:x2,:x3],LinearRelative(Normal()), graphinit=graphinit)

  addVariable!(fg, :l1, ContinuousScalar)
  addFactor!(fg, [:x1,:l1], LinearRelative(Normal()) , graphinit=graphinit)
  addFactor!(fg, [:x2,:l1], LinearRelative(Normal()) , graphinit=graphinit)

  addVariable!(fg, :l2, ContinuousScalar)
  addFactor!(fg, [:x3,:l2], LinearRelative(Normal()), graphinit=graphinit)

  return fg
end

"""
    $SIGNATURES

Canonical example introduced by Borglab.

Notes
- Known variable ordering: p = [:x1; :l3; :l1; :x5; :x2; :l2; :x4; :x3]
"""
function generateGraph_TestSymbolic(;graphinit::Bool=false)
  fg = initfg()

  addVariable!(fg, :x1, ContinuousScalar)
  addVariable!(fg, :x2, ContinuousScalar)
  addVariable!(fg, :x3, ContinuousScalar)
  addVariable!(fg, :x4, ContinuousScalar)
  addVariable!(fg, :x5, ContinuousScalar)
  addVariable!(fg, :l1, ContinuousScalar)
  addVariable!(fg, :l2, ContinuousScalar)
  addVariable!(fg, :l3, ContinuousScalar)

  addFactor!(fg, [:x1;:l1], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x1;:x2], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x2;:l1], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x2;:x3], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x3;:x4], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x4;:l2], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x4;:x5], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:l2;:x5], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x4;:l3], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x5;:l3], LinearRelative(Normal()), graphinit=graphinit)

  return fg
end

"""
    $SIGNATURES

Generate generalized helix parameterized by a curve along "t-axis" (i.e. z-axis, assuming z(t)=t).

Notes
- Returns vectors for (`t`, `x,y`, and `yaw` angle).
- Offset to t_start at origin and facing direction along +y-axis.
- Use callbacks `xr_t(t)` and `yr_t(t)` to skew the helix with any desired curve, examples include
  - `xr_t = (t) -> (1/3)t` to generate helix pattern along x-axis,
  - or make spiral along t using xr_t, yr_t to generate a rose pattern on xy,
  - use `spine_t(t)=xr_t(t) + im*yr_t(t)` as shortcut for more complicated patterns,
  - note `xr_t` and `yr_t` are scaled by a factor `radius`, unscale the input by division if desired.
- Use the function twice for simulated and noisy trajectories (i.e. easier Gauss-Markov processes)
- Gradient (i.e. angle) calculations are on the order of 1e-8.

Related

[`RoME.generateGraph_Helix2D!`](@ref)
"""
function calcHelix_T( t_start::Real=0,
                      t_stop::Real=1,
                      pointsperturn=20;
                      direction::Int=-1,
                      T::AbstractVector{<:Real}=(t_start:(t_stop*pointsperturn))./pointsperturn,
                      radius::Real = 0.5,
                      spine_t=(t)->0 + im*0,
                      xr_t::Function=(t)->real(spine_t(t)),
                      yr_t::Function=(t)->imag(spine_t(t)),
                      h::Real=1e-8  )
  #
  # calc the position
  f(t, x=xr_t(t), y=yr_t(t)) = radius*( cis(pi + direction*2pi*t) + 1 + x + im*y)
  vals = f.(T)

  # calc the gradient
  g(t) = (f(t+h)-f(t))/h
  grad = g.(T)

  return T, hcat(real.(vals), imag.(vals)), angle.(grad)
end

"""
    $SIGNATURES

Canonical example introduced originally as Caesar Hex Example.

Notes
- Paper variable ordering: p = [:x0;:x2;:x4;:x6;:x1;:l1;:x5;:x3;]
"""
function generateGraph_CaesarRing1D(;graphinit::Bool=false)

  fg = initfg()

  addVariable!(fg, :x0, ContinuousScalar)
  addVariable!(fg, :x1, ContinuousScalar)
  addVariable!(fg, :x2, ContinuousScalar)
  addVariable!(fg, :x3, ContinuousScalar)
  addVariable!(fg, :x4, ContinuousScalar)
  addVariable!(fg, :x5, ContinuousScalar)
  addVariable!(fg, :x6, ContinuousScalar)

  addFactor!(fg, [:x0], Prior(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x0;:x1], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x1;:x2], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x2;:x3], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x3;:x4], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x4;:x5], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x5;:x6], LinearRelative(Normal()), graphinit=graphinit)

  addVariable!(fg, :l1, ContinuousScalar)
  addFactor!(fg, [:x0;:l1], LinearRelative(Normal()), graphinit=graphinit)
  addFactor!(fg, [:x6;:l1], LinearRelative(Normal()), graphinit=graphinit)

  return fg
end

"""
    $SIGNATURES
Continuous, linear scalar and multivariate test graph generation. Follows a line
with the pose id equal to the ground truth.
"""
function generateGraph_LineStep(  lineLength::Int;
                                        poseEvery::Int=2,
                                        landmarkEvery::Int=4,
                                        posePriorsAt = Int[0],
                                        landmarkPriorsAt = Int[],
                                        sightDistance::Int=4,
                                        vardims = 1,
                                        noisy = false,
                                        graphinit = false,
                                        σ_pose_prior = 0.1,
                                        σ_lm_prior = 0.1,
                                        σ_pose_pose = 0.1,
                                        σ_pose_lm = 0.1,
                                        solverParams=SolverParams())
                                        # solverParams=SolverParams(algorithms=[:default, :parametric]))

  vtype = (vardims == 1) ? ContinuousScalar() : ContinuousEuclid(vardims)

  fg = LocalDFG{SolverParams}( solverParams=solverParams)

  function xNoise(i::Int, σ::Float64=1.0)
    if (vardims == 1)
      return noisy ? Normal(σ*randn() + i, σ) : Normal(0.0*randn() + i, σ)
    else
      return noisy ? MvNormal(σ*randn(vardims) .+ i, σ) : MvNormal(0.0*randn(vardims) .+ i, σ)
    end
  end

  x = Int[]
  lm = Int[]
  for i=0:lineLength
    if mod(i,poseEvery) == 0
      push!(x, i)
      addVariable!(fg, Symbol("x",i), vtype) #, autoinit = graphinit)
      (i in posePriorsAt) && addFactor!(fg, [Symbol("x",i)], Prior(xNoise(i, σ_pose_prior)), graphinit=graphinit)
      # "odo" type
      (i > 0) && addFactor!(fg, [Symbol("x",i-poseEvery); Symbol("x",i)], LinearRelative(xNoise(poseEvery, σ_pose_pose)), graphinit=graphinit)
    end

    if landmarkEvery != 0 && mod(i,landmarkEvery) == 0
      push!(lm, i)
      addVariable!(fg, Symbol("lm",i), vtype) #, autoinit = graphinit)
      (i in landmarkPriorsAt) && addFactor!(fg, [Symbol("lm",i)], Prior(xNoise(i, σ_lm_prior)), graphinit=graphinit)
    end
  end

  #add landmarks sightings
  for xi in x, lmi in lm
    dist = lmi - xi
    if abs(dist) < sightDistance
      # @info "adding landmark lm$lmi to x$xi with dist $dist"
      addFactor!(fg, [Symbol("x",xi); Symbol("lm",lmi)], LinearRelative(xNoise(dist, σ_pose_lm)), graphinit=graphinit)
    end
  end

  return fg
end

"""
    $SIGNATURES
Generate a EuclidDistance test graph where 1 landmark position is unknown.
"""
function generateGraph_EuclidDistance(points::Vector{Vector{Float64}} = [[100.,0],[0.,100]];
                                            dist = 100.0,
                                            σ_prior=1.0, σ_dist=1.0,
                                            N=100,
                                            graphinit=false)
  #
  dims = length(points[1])
  fg = initfg()
  fg.solverParams.N = N
  fg.solverParams.graphinit = graphinit

  for (i,p) in enumerate(points)
  xlbl = Symbol("x",i)
  addVariable!(fg, xlbl, ContinuousEuclid{dims})
  addFactor!(fg, [xlbl], Prior(MvNormal(p, σ_prior*ones(dims))))
  end

  addVariable!(fg, :l1, ContinuousEuclid{dims})

  for i=1:length(points)
  xlbl = Symbol("x",i)
  addFactor!(fg, [xlbl;:l1], EuclidDistance(Normal(dist, σ_dist)))
  end

  return fg
end


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

mutable struct VpnPacketCaptureStopParameters <: SwaggerModel
    sasUrl::Any # spec type: Union{ Nothing, String } # spec name: sasUrl

    function VpnPacketCaptureStopParameters(;sasUrl=nothing)
        o = new()
        validate_property(VpnPacketCaptureStopParameters, Symbol("sasUrl"), sasUrl)
        setfield!(o, Symbol("sasUrl"), sasUrl)
        o
    end
end # type VpnPacketCaptureStopParameters

const _property_map_VpnPacketCaptureStopParameters = Dict{Symbol,Symbol}(Symbol("sasUrl")=>Symbol("sasUrl"))
const _property_types_VpnPacketCaptureStopParameters = Dict{Symbol,String}(Symbol("sasUrl")=>"String")
Base.propertynames(::Type{ VpnPacketCaptureStopParameters }) = collect(keys(_property_map_VpnPacketCaptureStopParameters))
Swagger.property_type(::Type{ VpnPacketCaptureStopParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VpnPacketCaptureStopParameters[name]))}
Swagger.field_name(::Type{ VpnPacketCaptureStopParameters }, property_name::Symbol) =  _property_map_VpnPacketCaptureStopParameters[property_name]

function check_required(o::VpnPacketCaptureStopParameters)
    true
end

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


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type DoublePipe_HeatTransfer
	DoublePipe_HeatTransfer()=begin
		new(
			area(Dict{Symbol,Any}(
				:Brief=>"Cross Sectional Area for Flow",
				:Default=>0.05,
				:Lower=>1e-8
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Hydraulic Diameter of Pipe for Heat Transfer",
				:Lower=>1e-8
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Reynolds Number",
				:Default=>100,
				:Lower=>1
			)),
			heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Film Coefficient",
				:Default=>1,
				:Lower=>1e-12,
				:Upper=>1e6
			)),
			fricfactor(Dict{Symbol,Any}(
				:Brief=>"Friction Factor",
				:Default=>0.05,
				:Lower=>1e-10,
				:Upper=>2000
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Nusselt Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Prandtl Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Phi Correction",
				:Default=>1,
				:Lower=>1e-3
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Tube Velocity",
				:Lower=>1e-8
			)),
			[:As,:Dh,],
			[:Re,:hcoeff,:fi,:Nu,:PR,:Phi,:Vmean,]
		)
	end
	As::area
	Dh::length
	Re::positive
	hcoeff::heat_trans_coeff
	fi::fricfactor
	Nu::positive
	PR::positive
	Phi::positive
	Vmean::velocity
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export DoublePipe_HeatTransfer
function atributes(in::DoublePipe_HeatTransfer,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
DoublePipe_HeatTransfer(_::Dict{Symbol,Any})=begin
	newModel=DoublePipe_HeatTransfer()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
include("k_n.jl")
include("k_br.jl")
include("k_kappa.jl")
include("k_f.jl")

include("k_tb.jl")
include("k_kf.jl")
include("k_wl.jl")
include("k_dv.jl")
include("k_pc.jl")

include("k_p.jl")
include("k_pb.jl")

include("k_dr.jl")
include("k_do.jl")

include("k_vi.jl")
include("k_vt.jl")

include("k_rp.jl")
include("k_lp.jl")

include("k_cs.jl")

include("k_bs.jl")
include("k_cd.jl")

include("k_w.jl")

include("k_omega.jl")

include("k_wm.jl")

include("k_nt.jl")


text:
using SafeTestsets

@safetestset "Testing Finsts" begin
    using  VisualSearchACTR, Test
    ex = Experiment()
    target,present = initialize_trial!(ex)
    visual_objects = ex.populate_visicon(ex, target..., present)
    T = typeof(visual_objects)
    visual_location = VisualLocation(buffer=T)
    visicon = visual_objects
    visual_location.iconic_memory = visual_objects
    target_chunk = Chunk(;target...)
    goal = Goal(buffer=target_chunk)
    visual = Visual(buffer=T)
    actr = ACTRV(;T=Parm, goal, visual_location, visual, visicon)
    iconic_memory = get_iconic_memory(actr)
    map(x->x.attended=true, iconic_memory)
    map(x->x.attend_time = 0.0, iconic_memory)
    iconic_memory[1].attend_time = 1.1
    iconic_memory[2].attend_time = 1.1
    iconic_memory[3].attend_time = 1.1
    iconic_memory[4].attend_time = 1.1
    iconic_memory[6].attend_time = 2.0
    VisualSearchACTR.update_finst!(actr)
    @test iconic_memory[6].attend_time == 2.0
    @test iconic_memory[6].attended
    @test sum(x->x.attended, iconic_memory) == actr.parms.n_finst
    @test maximum(x->x.attend_time, iconic_memory) <= actr.parms.finst_span
end

@safetestset "Testing Iconic Decay" begin
    using  VisualSearchACTR, Test
    ex = Experiment()
    target,present = initialize_trial!(ex)
    visual_objects = ex.populate_visicon(ex, target..., present)
    T = typeof(visual_objects)
    visual_location = VisualLocation(buffer=T)
    visicon = visual_objects
    visual_location.iconic_memory = visual_objects
    target_chunk = Chunk(;target...)
    goal = Goal(buffer=target_chunk)
    visual = Visual(buffer=T)
    actr = ACTRV(;T=Parm, goal, visual_location, visual, persistence=1.0,
        visicon, time=4.0)
    iconic_memory = get_iconic_memory(actr)
    visible_objects = iconic_memory[1:5]
    map(x->x.features.color.visible=true, visible_objects[1:2])
    map(x->x.features.color.fixation_time=5.0, visible_objects[1:2])
    map(x->x.features.shape.visible=true, visible_objects[3:end])
    map(x->x.features.shape.fixation_time=6.0, visible_objects[3:end])
    map(x->VisualSearchACTR.object_visibility!(x), iconic_memory)
    @test sum(x->x.visible, iconic_memory) == 5
    VisualSearchACTR.update_decay!(actr)
    map(x->VisualSearchACTR.object_visibility!(x), iconic_memory)
    @test sum(x->x.visible, iconic_memory) == 5
    actr.time = 6.5
    VisualSearchACTR.update_decay!(actr)
    map(x->VisualSearchACTR.object_visibility!(x), iconic_memory)
    @test sum(x->x.visible, iconic_memory) == 3
end

@safetestset "Testing Activation" begin
    using  VisualSearchACTR, Test

    vals = [(color = :gray, shape = :p, x = 338, y =515),
            (color = :black, shape = :q, x = 351 , y = 193),
            (color = :gray, shape = :q, x = 511 , y = 462),
            (color = :black, shape = :p, x = 511 , y = 400)]

    experiment = Experiment(n_trials=10^1, trace=true, visible=false)

    target = (color=:gray,shape=:p)
    visual_objects = [VisualObject(features=(color=Feature(;value=x.color), shape=Feature(;value=x.shape)),
        width=32.0, location=[x.x-0.,x.y-0.]) for x in vals]

    T = typeof(visual_objects)
    visual_location = VisualLocation(buffer=T)
    visicon = visual_objects
    visual_location.iconic_memory = visual_objects
    target_chunk = Chunk(;target...)
    goal = Goal(buffer=target_chunk)
    visual = Visual(buffer=T)
    actr = ACTRV(;T=Parm, goal, visual_location, visual, visicon, σ=0.0)
    iconic_memory = get_iconic_memory(actr)

    # model = Model(;target=target, iconic_memory=visicon, noise=0.0)
    VisualSearchACTR.compute_angular_size!(actr)
    actr.visual.focus = [324.0,324.0]
    VisualSearchACTR.update_decay!(actr)
    VisualSearchACTR.update_finst!(actr)
    VisualSearchACTR.update_visibility!(actr)
    VisualSearchACTR.compute_activations!(actr)

    correct_top_down = [2,0,1,1]
    @test all(x->x[2].topdown_activation == x[1], zip(correct_top_down, iconic_memory))

    vo1 = iconic_memory[1]
    vo2 = iconic_memory[2]
    distance = VisualSearchACTR.compute_distance(vo1, vo2)
    bottomup_activation1 = VisualSearchACTR.bottomup_activation(vo1.features, vo2.features, distance)
    @test bottomup_activation1 == 1/sqrt(distance)*2

    vo1 = iconic_memory[1]
    vo3 = iconic_memory[3]
    distance = VisualSearchACTR.compute_distance(vo1, vo3)
    bottomup_activation2 = VisualSearchACTR.bottomup_activation(vo1.features, vo3.features, distance)
    @test bottomup_activation2 == 1/sqrt(distance)

    vo1 = iconic_memory[1]
    vo4 = iconic_memory[4]
    distance = VisualSearchACTR.compute_distance(vo1, vo4)
    bottomup_activation3 = VisualSearchACTR.bottomup_activation(vo1.features, vo4.features, distance)
    @test bottomup_activation3 == 1/sqrt(distance)

    @test vo1.bottomup_activation == (bottomup_activation1 + bottomup_activation2 + bottomup_activation3)
end

@safetestset "Testing Feature Search" begin
    using VisualSearchACTR, Test, DataFrames, GLM, Statistics
    include("simulation.jl")
    Random.seed!(95025181)
    set_sizes = [1,2,5,10,20,30]
    Δτ = .8
    topdown_weight = .60
    noise = true
    rnd_time = true
    results = run_simulation(set_sizes; fun=feature_set, Δτ, topdown_weight, noise, rnd_time)

    df_present = filter(x->x[:target_present] ==:present && x[:response] ==:present, results)
    ols_present = lm(@eval(@formula(rt_mean ~ distractors)), df_present)
    β0,β1 = coef(ols_present)
    @test β0 ≈ 0.36 rtol = .05
    @test β1 ≈ 0.0 atol = .01
    df_absent = filter(x->x[:target_present] ==:absent && x[:response] ==:absent, results)
    ols_absent = lm(@eval(@formula(rt_mean ~ distractors)), df_absent)
    β0,β1 = coef(ols_absent)
    @test β0 ≈ 0.45 rtol = .05
    @test β1 ≈ 0.0 atol = .01
end

@safetestset "Testing Conjunctive Search" begin
    using VisualSearchACTR, Test, DataFrames, GLM, Statistics
    include("simulation.jl")
    Random.seed!(52484)
    set_sizes = [1,2,5,10,20,30]
    Δτ = .4
    topdown_weight = .60
    noise = true
    rnd_time = true
    results = run_simulation(set_sizes; fun=conjunctive_set, Δτ, topdown_weight, noise, rnd_time)

    df_present = filter(x->x[:target_present] ==:present && x[:response] ==:present, results)
    ols_present = lm(@eval(@formula(rt_mean ~ distractors)), df_present)
    β0,β1 = coef(ols_present)
    @test β0 ≈ 0.34 rtol = .05
    @test β1 ≈ 0.020 atol = .01
    df_absent = filter(x->x[:target_present] ==:absent && x[:response] ==:absent, results)
    ols_absent = lm(@eval(@formula(rt_mean ~ distractors)), df_absent)
    β0,β1 = coef(ols_absent)
    @test β0 ≈ 0.38 rtol = .05
    @test β1 ≈ 0.067 atol = .01
end

@safetestset "Gamma Parms" begin
    using VisualSearchACTR, Test, Distributions
    import VisualSearchACTR: gamma_parms

    θ = gamma_parms(3, 2)
    @test mean(Gamma(θ...)) ≈ 3 rtol = .0005
    @test std(Gamma(θ...)) ≈ 2 rtol = .0005

    θ = gamma_parms(.2, 1)
    @test mean(Gamma(θ...)) ≈ .2 rtol = .0005
    @test std(Gamma(θ...)) ≈ 1 rtol = .0005
end

@safetestset "Run model" begin
    using VisualSearchACTR, Test
    experiment = Experiment(set_size=10,  n_trials=2,
        trace=true, speed =.5)
    run_condition!(experiment)
    @test true
end


text:
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

code = @dice begin
    # all tails
    probs = [1/i for i=2:20]
    mapreduce(p -> !flip(p), &, probs)
end

# BDD analysis
bdd = compile(code)
num_flips(bdd), num_nodes(bdd)
infer(code, :bdd)

# IR analysis
println(to_dice_ir(code))
has_dice_binary() && rundice(code)
has_dice_binary() && infer(code, :ocaml)


text:
#unum-comparison.jl

import Base:  ==, <, >, <=, >=  #this statement is necessary to redefine these functions directly

@universal function ==(a::Unum, b::Unum)
  #first compare the ubits.... These must be the same or else they aren't equal.
  ((a.flags & UNUM_UBIT_MASK) != (b.flags & UNUM_UBIT_MASK)) && return false
  #next, compare the sign bits...  The only one that spans is zero.
  ((a.flags & UNUM_SIGN_MASK) != (b.flags & UNUM_SIGN_MASK)) && (return (is_zero(a) && is_zero(b)))
  #check if either is nan.  Note that NaN != NaN.
  (isnan(a) || isnan(b)) && return false

  #resolve degenerate forms.
  resolve_degenerates!(a)
  resolve_degenerates!(b)

  #first calculate the exponents.
  _aexp::Int64 = decode_exp(a)
  _bexp::Int64 = decode_exp(b)

  (is_exp_zero(a) != is_exp_zero(b)) && return false

  (is_exp_zero(a) && is_exp_zero(b) && is_frac_zero(a) && is_frac_zero(b)) && return true

  (_aexp != _bexp) && return false
  #now that we know that the exponents are the same,
  #the fractions must also be identical
  (a.fraction != b.fraction) && return false

  #the ubit on case is simple - the fsizes must be equal, except if one is mmr.
  if ((a.flags & UNUM_UBIT_MASK) == UNUM_UBIT_MASK)
    is_inf_ulp(a) || return (a.fsize == b.fsize)
    (a.fsize == b.fsize) && return true
    #mmr is categorically going to be equal to the inf_ulp with fsize one less
    #than max_fsize
    _mfsize = max_fsize(FSS)
    (a.fsize == _mfsize) && (b.fsize == _mfsize - o16) && return true
    (b.fsize == _mfsize) && (a.fsize == _mfsize - o16) && return true
    return false
  end
  #so we are left with exact floats at any resolution, ok, or uncertain floats
  #with the same resolution.  These must be equal.
  return true
end

#note that unlike IEEE floating points, 0 === -0 for Unums.  So the only isequal
#exception should be the NaN exception.  In order to achieve this, we take negative
#zero and make it the degenerate zero.  We also collapse fsize and esize to the
#most reasonable value.
@universal function Base.hash(a::Unum, h::UInt)
  #use the full_decode function to extract
  resolve_degenerates!(a)
  is_zero(a) && is_negative(a) && (a.flags &= ~UNUM_SIGN_MASK)

  #now generate the hash.
  h = hash(a.esize << 32 | a.fsize << 16 | a.flags, h)
  h = hash(a.fraction, h)
  h = hash(a.exponent, h)
  h
end

#the corresponding isequal function.
@universal Base.isequal(a::Unum, b::Unum) = (hash(a) == hash(b))

doc"""
  `Unums.is_inward(a::Unum, b::Unum)` tells if a is inward of b.  If
  a and b are both positive,` is_inward(a, b) === a < b`, and if
  a and b are both negative, `is_inward(a, b) === b < a`.
"""
@universal function is_inward(a::Unum, b::Unum)
  #resolve degeneracies to make comparison much much easier..

  resolve_degenerates!(a)
  resolve_degenerates!(b)

  _as = is_subnormal(a)
  _bs = is_subnormal(b)

  _aexp = decode_exp(a)
  _bexp = decode_exp(b)

  #so now we know that these two have the same sign.
  (_aexp < _bexp) && return true

  #check for the strange subnormal ulps.
  if is_ulp(a) && is_strange_subnormal(a)
    is_ulp(b) && is_strange_subnormal(b) && return false
    #redo _aexp
    _aexp -= a.fsize + 1
    return _aexp < _bexp
  end

  (_aexp > _bexp) && return false

  is_exact(a) & is_exact(b) && return (a.fraction < b.fraction)

  #in the case that b is an ulp and a is exact, then exact equality means a is
  #less, otherwise exact equality does not.
  orequal = is_ulp(b) & is_exact(a)

  lessthanwithubit(a.fraction, b.fraction, a.fsize, orequal)
end

#TODO:  reimplement this as isless.
@universal function >(a::Unum, b::Unum)
  (is_nan(a) || is_nan(b)) && return false
  _b_pos::Bool = (is_positive(b))
  _a_pos::Bool = (is_positive(a))
  _b_zer::Bool = is_zero(b)
  _a_zer::Bool = is_zero(a)
  _a_inf::Bool = is_inf(a)

  (_b_pos) & (!_a_pos) && return false
  (!_b_pos) & (_a_pos) && return (!(_a_zer & _b_zer))

  #zero can cause problems down the line.
  _a_zer && return !(_b_pos | _b_zer)
  _b_zer && return _a_pos

  #infinity has strange interactions with some inf_ulp values.
  (_a_inf & _a_pos) && return !(is_inf(b) & _b_pos)  #if a is positive infinity, greater than only if b is not.
  (_a_inf & !_a_pos) && return false  #if a is negative infinity, never greater than

  #use this as a trampoline for is_inward.
  _a_pos && return is_inward(b, a)

  return is_inward(a, b)
end
#hopefully the julia compiler knows what to do here.
@universal <(a::Unum, b::Unum) = (b > a)
@universal <=(a::Unum, b::Unum) = !(a > b)
@universal >=(a::Unum, b::Unum) = !(b > a)

@universal function inner(a::Unum, b::Unum)
  #adjust them in case they are subnormal
  resolve_degenerates!(a)
  resolve_degenerates!(b)

  #next criterion, are the exponents different
  _aexp = decode_exp(a)
  _bexp = decode_exp(b)
  if (_aexp != _bexp)
    return (_aexp < _bexp) ? a : b
  end

  #next criteria, check the fractions
  if (a.fraction != b.fraction)
    return (a.fraction < b.fraction) ? a : b
  end

  #finally, check the ubit
  return is_ulp(a) ? a : b
end

@universal function outer(a::Unum, b::Unum)
  #adjust them in case they are subnormal
  resolve_degenerates!(a)
  resolve_degenerates!(b)

  #next criterion, are the exponents different
  _aexp = decode_exp(a)
  _bexp = decode_exp(b)
  if (_aexp != _bexp)
    return (_aexp > _bexp) ? a : b
  end

  #next criteria, check the fractions
  lessthanwithubit(a.fraction, b.fraction, min(a.fsize, b.fsize)) ? b : a
end

@universal function Base.min(a::Unum, b::Unum)
  #first, fastest criterion:  Are they not the same sign?
  if (a.flags $ b.flags) & UNUM_SIGN_MASK != 0
    return is_negative(a) ? a : b
  end

  return is_negative(a) ? outer(a, b) : inner(a, b)
end

@universal function Base.max(a::Unum, b::Unum)
  #first, fastest criterion:  Are they not the same sign?
  if (a.flags $ b.flags) & UNUM_SIGN_MASK != 0
    return is_negative(a) ? b : a
  end

  return is_negative(a) ? inner(a, b) : outer(a, b)
end

doc"""
  `Unums.mag_greater_than_one(a::Unum)`
  outputs true if the magnitude of the passed unum is greater than one.
"""
@universal mag_greater_than_one(a::Unum) = decode_exp(a) >= 0
mag_greater_than_one{FSS}(a::UnumSmall{0,FSS}) = (a.exponent != z64) | frac_top_bit(a)
mag_greater_than_one{FSS}(a::UnumLarge{0,FSS}) = (a.exponent != z64) | frac_top_bit(a)

@universal function ≸(a::Unum, b::Unum)
  !(a < b) && !(b < a)
end

≹(a::Unum, b::Unum) = ≸(a, b)

#simless
@universal function notless(a::Unum, b::Unum)
  !(a < b)
end
≮ = notless
#simgtr
@universal function notgtr(a::Unum, b::Unum)
  !(a > b)
end
≯ = notgtr

export ≸,≹,≯,≮

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
  import JSON

  type Emitter
    io::IO
    cache::Cache
  end

  function make_emitter(io, verbose::Bool)
    let cache = verbose ? NoopCache() : RollingCache()
      Emitter(io, cache)
    end
  end

  function emit_raw(e::Emitter, s::AbstractString)
    print(e.io, s)
  end

  function emit_tag(e::Emitter, x::AbstractString)
    emit(e, "~$x", true)
  end

  function emit(e::Emitter, x::AbstractString, cacheable::Bool)
    let value = iscacheable(x, cacheable) ? write!(e.cache, x) : x
        print(e.io, JSON.json(value))
    end
  end

  function emit(e::Emitter, x::Integer)
    print(e.io, JSON.json(x))
  end

  function emit(e::Emitter, x::Bool)
    print(e.io, JSON.json(x))
  end

  function emit_null(e::Emitter, askey::Bool)
    askey ? emit_tag(e, "_") : emit_raw(e, "null")
  end

  function emit_array_start(e::Emitter)
    print(e.io, "[")
  end

  function emit_array_end(e::Emitter)
    print(e.io, "]")
  end

  function emit_map_start(e::Emitter)
    print(e.io, "{")
  end

  function emit_map_sep(e::Emitter, i=2)
    emit_array_sep(e, i)
  end

  function emit_key_sep(e::Emitter)
    print(e.io, ":")
  end

  function emit_map_end(e::Emitter)
    print(e.io, "} ")
  end

  function emit_array_sep(e::Emitter, i=2)
    if i != 1
      print(e.io, ",")
    end
  end


text:
Transducers.transduce(xf, rf, init, xs, ex::TaskPoolEx) =
    transduce_taskpool(xf, rf, init, xs; ex.kwargs...)

function transduce_taskpool(
    xf::Transducer,
    rf,
    init,
    xs;
    simd = Val(false),
    basesize = nothing,
    ntasks::Int = Threads.nthreads(),
    background::Bool = false,
)
    rf0 = _reducingfunction(xf, rf; init = init, simd = simd)
    rf1, xs0 = retransform(rf0, xs)
    if basesize === nothing
        basesize = amount(xs0) ÷ Threads.nthreads()
    end
    xs1 = SizedReducible(xs0, basesize)
    return transduce_dac(
        TaskPoolScheduler(ntasks, background),
        CancellableDACContext(),
        rf1,
        init,
        xs1,
    )
end

struct TaskPoolScheduler <: AbstractScheduler
    ntasks::Int
    semaphore::Union{Base.Semaphore,Nothing}
    request::Channel{Any}
end

TaskPoolScheduler(ntasks::Int, background::Bool) = TaskPoolScheduler(
    ntasks,
    ntasks > Threads.nthreads() ? nothing : Base.Semaphore(ntasks),
    background ? get_background_taskpool_request() : get_taskpool_request(),
)

function spawn(@nospecialize(f), sch::TaskPoolScheduler)
    semaphore = sch.semaphore
    if semaphore === nothing
        put!(sch.request, f)
    else
        Base.acquire(semaphore)
        function wrapper()
            try
                return f()
            finally
                Base.release(semaphore)
            end
        end
        put!(sch.request, wrapper)
    end
    return
end

const TASKPOOL_REQUEST = Ref{typeof(Channel{Any}(Inf))}()
const BACKGROUND_TASKPOOL_REQUEST = Ref{typeof(Channel{Any}(Inf))}()

get_taskpool_request() = _get_taskpool_request(TASKPOOL_REQUEST, init_taskpool)

function get_background_taskpool_request()
    if Threads.nthreads() == 1
        error(
            "background thread not usable when Julia is not started with multiple threads.",
        )
    end
    return _get_taskpool_request(BACKGROUND_TASKPOOL_REQUEST, init_background_taskpool)
end

function _get_taskpool_request(ref, init)
    ch = ref[]
    isopen(ch) && return ch
    @error "taskpool closed; trying to recover..."
    # TODO: lock
    init()
    ch = ref[]
    isopen(ch) && return ch
    error("failed to create a root task pool")
end

function init_taskpool()
    @debug "`init_taskpool`: Initializing `TASKPOOL_REQUEST` channels..."
    TASKPOOL_REQUEST[] = request = request_channel()
    start_taskpool(() -> true, request)
    @debug "`init_taskpool`: Initializing `TASKPOOL_REQUEST` channels... DONE"
end

function init_background_taskpool()
    @debug "`init_background_taskpool`: Initializing `BACKGROUND_TASKPOOL_REQUEST` channels..."
    BACKGROUND_TASKPOOL_REQUEST[] = request = request_channel()
    start_taskpool(() -> Threads.threadid() != 1, request)
    @debug "`init_background_taskpool`: Initializing `BACKGROUND_TASKPOOL_REQUEST` channels... DONE"
end

function start_taskpool(p, request)
    ok = foreach_thread() do
        p() || return
        @async request_handler(request)
    end
    if !ok
        @warn "failed to start task pool workers on all threads"
    end
end


text:
using ExamplePackage
using Test

@test foo() == 0.11388071406436832
@test foo(1, 1.5) == 0.2731856314283442
@test_broken foo(1, 0) # tells us this is broken


text:
# Arvind Prasadan
# 2019 January
# DMF Package
# Tests for functions in package

push!(LOAD_PATH, "../src/")
using DMF
using LinearAlgebra
using Statistics
using StatsBase
using Test

# gen_arma_sequence, no arma
tol = 1e-2
n = 10000
ar = []
ma = []
arma_std = 1.0
trials = 200
running_sum = 0.0
for tr = 1:1:trials
    x = gen_arma_sequence(n, ar, ma, arma_std)
    global running_sum += std(x) / trials
end
@test isapprox(arma_std, running_sum; atol = tol)

# gen_arma_sequence, AR(1)
tol = 1e-2
n = 10000
ar = 0.5
ma = []
arma_std = 1.0
trials = 200
running_sum = 0.0
for tr = 1:1:trials
    x = gen_arma_sequence(n, ar, ma, arma_std)
    global running_sum += std(x) / trials
end
@test isapprox(arma_std / sqrt(1.0 - ar^2.0), running_sum; atol = tol)

# gen_cos_sequence, element
tol = 1e-8
n = 1000
f = 2.0
fs = 1.0
phase = 0.1
x, t = gen_cos_sequence(n, f, fs, phase)
@test isapprox(x[2], cos((f / fs) + phase); atol = tol)
@test isapprox(t[2], 1.0 / fs; atol = tol)

# autocorrelation, standard
tol = 1e-8
x = [1; 1; 1; 1]
gamma = autocorrelation(x; circular = false, bias_correct = false)
@test isapprox(gamma[1], 0.0; atol = tol)
@test isapprox(gamma[2], 0.0; atol = tol)

# autocorrelation, bias correction
tol = 1e-8
x = [1; 1; 1; 1]
gamma = autocorrelation(x; circular = false, bias_correct = true)
@test isapprox(gamma[1], 0.0; atol = tol)
@test isapprox(gamma[2], 0.0; atol = tol)

# autocorrelation, circular
tol = 1e-8
x = [1; 1; 1; 1]
gamma = autocorrelation(x; circular = true, bias_correct = true)
@test isapprox(gamma[1], 0.0; atol = tol)
@test isapprox(gamma[2], 0.0; atol = tol)

# autocorr_mat
tol = 1e-2
p = 5
n = 10000
trials = 200
running_sum = zeros(1, p, p)
for tr = 1:1:trials
    global running_sum += autocorr_mat(randn(p, n), 1)[1] / trials
end
@test isapprox(running_sum[1, :, :], zeros(p, p); atol = tol)

# eigenvector_error
tol = 1e-8
p = 10
Q = randn(p, 2)
Q_hat = [Q[:, 2] Q[:, 1] randn(p)]
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)

# For Source Separation Tests

p = 100
k = 2
n = 1000
Q = qr(randn(p, k)).Q
Q = Q[:, 1:k]
S = [gen_cos_sequence(n, 2.0)[1] gen_cos_sequence(n, 0.25)[1]]
S = mapslices(normalize, S; dims = 1)
D = Diagonal(1:k)
X = Q * D * S'

# SOBI_Wrapper
tol = 1e-2
Q_hat, S_hat = SOBI_Wrapper(X, 1)
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)
@test isapprox(eigenvector_error(S, S_hat), 0.0; atol = tol)

# DMF: Standard
tol = 1e-3
w, Q_hat, C_hat, A = dmf(X; tsvd = false, os = false, hil = false, nsv = k, C_nsv = k, demean = false, lag = 1)
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)
@test isapprox(eigenvector_error(S, C_hat), 0.0; atol = tol)
@test size(C_hat, 2) == k

# DMF: Hilbert
tol = 1e-3
w, Q_hat, C_hat, A = dmf(X; tsvd = false, os = false, hil = true, nsv = k, C_nsv = k, demean = false, lag = 1)
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)
@test isapprox(eigenvector_error(S, C_hat), 0.0; atol = tol)
@test size(C_hat, 2) == k

# DMF: tSVD
tol = 1e-3
w, Q_hat, C_hat, A = dmf(X; tsvd = true, os = false, hil = false, nsv = k, C_nsv = k, demean = false, lag = 1)
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)
@test isapprox(eigenvector_error(S, C_hat), 0.0; atol = tol)
@test size(C_hat, 2) == k

# DMF: OptShrink
tol = 1e-3
w, Q_hat, C_hat, A = dmf(X; tsvd = false, os = true, hil = false, nsv = k, C_nsv = k, demean = false, lag = 1)
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)
@test isapprox(eigenvector_error(S, C_hat), 0.0; atol = tol)
@test size(C_hat, 2) == k

# DMF: De-Mean
tol = 1e-3
w, Q_hat, C_hat, A = dmf(X; tsvd = false, os = false, hil = false, nsv = k, C_nsv = k, demean = true, lag = 1)
@test isapprox(eigenvector_error(Q, Q_hat), 0.0; atol = tol)
@test isapprox(eigenvector_error(S, C_hat), 0.0; atol = tol)
@test size(C_hat, 2) == k

# OptShrink, rank = p or n
tol = 1e-8
p = 10
n = 100
X = randn(p, n)
Xh, w = optshrink(X, p)
@test isapprox(X, Xh; atol = tol)
Xh, w = optshrink(X, n)
@test isapprox(X, Xh; atol = tol)
Xh, w = optshrink(X, n + p)
@test isapprox(X, Xh; atol = tol)

# OptShrink
tol = 1e-8
X = ones(p, n)
Xh, w = optshrink(X, 1)
@test isapprox(w[1], norm(ones(p * n), 2.0); atol = tol)
@test isapprox(X, Xh; atol = tol)


text:
@testset "Text" begin

@testset "sentences" begin
    @test [string(x) for x = sentences("Hello, world!")] ==
            ["Hello, world!"]
    @test [string(x) for x = sentences("Hello. This is a test.")] ==
            ["Hello.", "This is a test."]
    @test [string(x) for x = sentences(
            "Do question marks work? Yes, they do.")] ==
            ["Do question marks work?", "Yes, they do."]

    mobydick = replace(
        read(joinpath(dirname(@__DIR__), "text", "mobydick.txt"), String),
        "\n" => " ")
    result = [string(x) for x = sentences(mobydick)]
    @test result[1] == "Call me Ishmael."
    @test result[2] == string(
        "Some years ago — never mind how long precisely — having little or ",
        "no money in my purse, and nothing particular to interest me on ",
        "shore, I thought I would sail about a little and see the watery ",
        "part of the world.")
    @test result[end] == string(
        "If they but knew it, almost all men in their degree, some time or ",
        "other, cherish very nearly the same feelings towards the ocean ",
        "with me.")

    # test ending without punctuation
    @test [string(x) for x = sentences("Look! It's Iron Man")] ==
        ["Look!", "It's Iron Man"]
end

end


text:
#####
##### Building blocks for traversing a Hamiltonian deterministically, using the leapfrog
##### integrator.
#####

####
#### kinetic energy
####

"""
$(TYPEDEF)

Kinetic energy specifications. Implements the methods

- `Base.size`

- [`kinetic_energy`](@ref)

- [`calculate_p♯`](@ref)

- [`∇kinetic_energy`](@ref)

- [`rand_p`](@ref)

For all subtypes, it is implicitly assumed that kinetic energy is symmetric in
the momentum `p`,

```julia
kinetic_energy(κ, p, q) == kinetic_energy(κ, .-p, q)
```

When the above is violated, the consequences are undefined.
"""
abstract type KineticEnergy end

"""
$(TYPEDEF)

Euclidean kinetic energies (position independent).
"""
abstract type EuclideanKineticEnergy <: KineticEnergy end

"""
$(TYPEDEF)

Gaussian kinetic energy, with ``K(q,p) = p ∣ q ∼ 1/2 pᵀ⋅M⁻¹⋅p + log|M|`` (without constant),
which is independent of ``q``.

The inverse covariance ``M⁻¹`` is stored.

!!! note
    Making ``M⁻¹`` approximate the posterior variance is a reasonable starting point.
"""
struct GaussianKineticEnergy{T <: AbstractMatrix,
                             S <: AbstractMatrix} <: EuclideanKineticEnergy
    "M⁻¹"
    M⁻¹::T
    "W such that W*W'=M. Used for generating random draws."
    W::S
    function GaussianKineticEnergy(M⁻¹::T, W::S) where {T, S}
        @argcheck checksquare(M⁻¹) == checksquare(W)
        new{T,S}(M⁻¹, W)
    end
end

"""
$(SIGNATURES)

Gaussian kinetic energy with the given inverse covariance matrix `M⁻¹`.
"""
GaussianKineticEnergy(M⁻¹::AbstractMatrix) = GaussianKineticEnergy(M⁻¹, cholesky(inv(M⁻¹)).L)

"""
$(SIGNATURES)

Gaussian kinetic energy with the given inverse covariance matrix `M⁻¹`.
"""
GaussianKineticEnergy(M⁻¹::Diagonal) = GaussianKineticEnergy(M⁻¹, Diagonal(.√inv.(diag(M⁻¹))))

"""
$(SIGNATURES)

Gaussian kinetic energy with a diagonal inverse covariance matrix `M⁻¹=m⁻¹*I`.
"""
GaussianKineticEnergy(N::Integer, m⁻¹ = 1.0) = GaussianKineticEnergy(Diagonal(fill(m⁻¹, N)))

function Base.show(io::IO, κ::GaussianKineticEnergy{T}) where {T}
    print(io::IO, "Gaussian kinetic energy ($(nameof(T))), √diag(M⁻¹): $(.√(diag(κ.M⁻¹)))")
end

## NOTE about implementation: the 3 methods are callable without a third argument (`q`)
## because they are defined for Gaussian (Euclidean) kinetic energies.

Base.size(κ::GaussianKineticEnergy, args...) = size(κ.M⁻¹, args...)

"""
$(SIGNATURES)

Return kinetic energy `κ`, at momentum `p`.
"""
kinetic_energy(κ::GaussianKineticEnergy, p, q = nothing) = dot(p, κ.M⁻¹ * p) / 2

"""
$(SIGNATURES)

Return ``p♯ = M⁻¹⋅p``, used for turn diagnostics.
"""
calculate_p♯(κ::GaussianKineticEnergy, p, q = nothing) = κ.M⁻¹ * p

"""
$(SIGNATURES)

Calculate the gradient of the logarithm of kinetic energy in momentum `p`.
"""
∇kinetic_energy(κ::GaussianKineticEnergy, p, q = nothing) = calculate_p♯(κ, p)

"""
$(SIGNATURES)

Generate a random momentum from a kinetic energy at position `q`.
"""
rand_p(rng::AbstractRNG, κ::GaussianKineticEnergy, q = nothing) = κ.W * randn(rng, size(κ.W, 1))

####
#### Hamiltonian
####

struct Hamiltonian{K,L}
    "The kinetic energy specification."
    κ::K
    """
    The (log) density we are sampling from. Supports the `LogDensityProblem` API.
    Technically, it is the negative of the potential energy.
    """
    ℓ::L
    """
    $(SIGNATURES)

    Construct a Hamiltonian from the log density `ℓ`, and the kinetic energy specification
    `κ`. `ℓ` with a vector are expected to support the `LogDensityProblems` API, with
    gradients.
    """
    function Hamiltonian(κ::K, ℓ::L) where {K <: KineticEnergy,L}
        @argcheck capabilities(ℓ) ≥ LogDensityOrder(1)
        @argcheck dimension(ℓ) == size(κ, 1)
        new{K,L}(κ, ℓ)
    end
end

Base.show(io::IO, H::Hamiltonian) = print(io, "Hamiltonian with $(H.κ)")

"""
$(TYPEDEF)

A log density evaluated at position `q`. The log densities and gradient are saved, so that
they are not calculated twice for every leapfrog step (both as start- and endpoints).

Because of caching, a `EvaluatedLogDensity` should only be used with a specific Hamiltonian,
preferably constructed with the `evaluate_ℓ` constructor.

In composite types and arguments, `Q` is usually used for this type.
"""
struct EvaluatedLogDensity{T,S}
    "Position."
    q::T
    "ℓ(q). Saved for reuse in sampling."
    ℓq::S
    "∇ℓ(q). Cached for reuse in sampling."
    ∇ℓq::T
    function EvaluatedLogDensity(q::T, ℓq::S, ∇ℓq::T) where {T <: AbstractVector,S <: Real}
        @argcheck length(q) == length(∇ℓq)
        new{T,S}(q, ℓq, ∇ℓq)
    end
end

# general constructors below are necessary to sanitize input from eg Diagnostics, or an
# initial position given as integers, etc

function EvaluatedLogDensity(q::AbstractVector, ℓq::Real, ∇ℓq::AbstractVector)
    q, ∇ℓq = promote(q, ∇ℓq)
    EvaluatedLogDensity(q, ℓq, ∇ℓq)
end

EvaluatedLogDensity(q, ℓq::Real, ∇ℓq) = EvaluatedLogDensity(collect(q), ℓq, collect(∇ℓq))

"""
$(SIGNATURES)

Evaluate log density and gradient and save with the position. Preferred interface for
creating `EvaluatedLogDensity` instances.
"""
function evaluate_ℓ(ℓ, q)
    ℓq, ∇ℓq = logdensity_and_gradient(ℓ, q)
    if isfinite(ℓq)
        EvaluatedLogDensity(q, ℓq, ∇ℓq)
    else
        EvaluatedLogDensity(q, oftype(ℓq, -Inf), q) # second q used just as a placeholder
    end
end

"""
$(TYPEDEF)

A point in phase space, consists of a position (in the form of an evaluated log density `ℓ`
at `q`) and a momentum.
"""
struct PhasePoint{T <: EvaluatedLogDensity,S}
    "Evaluated log density."
    Q::T
    "Momentum."
    p::S
    function PhasePoint(Q::T, p::S) where {T,S}
        @argcheck length(p) == length(Q.q)
        new{T,S}(Q, p)
    end
end

"""
$(SIGNATURES)

Log density for Hamiltonian `H` at point `z`.

If `ℓ(q) == -Inf` (rejected), skips the kinetic energy calculation.

Non-finite values (incl `NaN`, `Inf`) are automatically converted to `-Inf`. This can happen
if

1. the log density is not a finite value,

2. the kinetic energy is not a finite value (which usually happens when `NaN` or `Inf` got
mixed in the leapfrog step, leading to an invalid position).
"""
function logdensity(H::Hamiltonian{<:EuclideanKineticEnergy}, z::PhasePoint)
    @unpack ℓq = z.Q
    isfinite(ℓq) || return oftype(ℓq, -Inf)
    K = kinetic_energy(H.κ, z.p)
    ℓq - (isfinite(K) ? K : oftype(K, Inf))
end

function calculate_p♯(H::Hamiltonian{<:EuclideanKineticEnergy}, z::PhasePoint)
    calculate_p♯(H.κ, z.p)
end

"""
    leapfrog(H, z, ϵ)

Take a leapfrog step of length `ϵ` from `z` along the Hamiltonian `H`.

Return the new phase point.

The leapfrog algorithm uses the gradient of the next position to evolve the momentum. If
this is not finite, the momentum won't be either, `logdensity` above will catch this and
return an `-Inf`, making the point divergent.
"""
function leapfrog(H::Hamiltonian{<: EuclideanKineticEnergy}, z::PhasePoint, ϵ)
    @unpack ℓ, κ = H
    @unpack p, Q = z
    @argcheck isfinite(Q.ℓq) "Internal error: leapfrog called from non-finite log density"
    pₘ = p + ϵ/2 * Q.∇ℓq
    q′ = Q.q + ϵ * ∇kinetic_energy(κ, pₘ)
    Q′ = evaluate_ℓ(H.ℓ, q′)
    p′ = pₘ + ϵ/2 * Q′.∇ℓq
    PhasePoint(Q′, p′)
end


text:
function AdVIL_π_loss(π, 𝒫, 𝒟; info=Dict())
    π_a = action(π, 𝒟[:s])
    mean(value(π, 𝒟[:s], π_a)) + 𝒫[:λ_BC]*Flux.mse(π_a, 𝒟[:a])
end

function AdVIL_D_loss(π, 𝒫, 𝒟, y; info=Dict())
    π_sa = vcat(𝒟[:s], action(π, 𝒟[:s]))
    expert_sa = vcat(𝒟[:s], 𝒟[:a])
    mean(value(π, expert_sa)) - mean(value(π, π_sa)) + 𝒫[:λ_GP]*gradient_penalty(critic(π), expert_sa, π_sa, target=0.4f0)
end

function AdVIL(;π,
                S,
                𝒟_demo,
                normalize_demo::Bool=true,
                λ_GP::Float32=10f0,
                λ_orth::Float32=1f-4,
                λ_BC::Float32=2f-1,
                a_opt::NamedTuple=(;),
                c_opt::NamedTuple=(;),
                log::NamedTuple=(;),
                kwargs...)

    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π)))
    𝒟_demo = 𝒟_demo |> device(π)

    BatchSolver(;agent=PolicyParams(π),
                 S=S,
                 𝒫=(λ_GP=λ_GP, λ_BC=λ_BC,),
                 𝒟_train = 𝒟_demo,
                 a_opt=TrainingParams(;name="actor_", loss=AdVIL_π_loss, regularizer=OrthogonalRegularizer(λ_orth), a_opt...),
                 c_opt=TrainingParams(;name="discriminator_", loss=AdVIL_D_loss, c_opt...),
                 log=LoggerParams(;dir="log/AdVIL", period=1, log...),
                 kwargs...)
end


text:
mutable struct six_axis_force_torque_array_t <: LCMType
    utime::Int64
    num_sensors::Int32
    names::Vector{String}
    sensors::Vector{six_axis_force_torque_t}
end

@lcmtypesetup(six_axis_force_torque_array_t,
    names => (num_sensors,),
    sensors => (num_sensors,)
)


text:
struct WeightedBagNode{T <: Union{Missing, AbstractNode}, B <: AbstractBags, W, C} <: AbstractBagNode
    data::T
    bags::B
    weights::Vector{W}
    metadata::C

    function WeightedBagNode(d::T, b::B, w::Vector{W}, m::C=nothing) where {T <: Union{Missing, AbstractNode}, B <: AbstractBags, W, C}
        ismissing(d) && any(_len.(b.bags) .!= 0) && error("WeightedBagNode with nothing in data cannot have a non-empty bag")
        new{T, B, W, C}(d, b, w, m)
    end
end

WeightedBagNode(x::T, b::Vector, weights::Vector{W}, metadata::C=nothing) where {T, W, C} =
WeightedBagNode(x, bags(b), weights, metadata)

mapdata(f, x::WeightedBagNode) = WeightedBagNode(mapdata(f, x.data), x.bags, x.weights, x.metadata)

function Base.getindex(x::WeightedBagNode{T, B, W}, i::VecOrRange) where {T, B, W}
    nb, ii = remapbag(x.bags, i)
    _emptyismissing[] && isempty(ii) && return(WeightedBagNode(missing, nb, W[], nothing))
    WeightedBagNode(subset(x.data,ii), nb, subset(x.weights, ii), subset(x.metadata, ii))
end

function reduce(::typeof(catobs), as::Vector{T}) where {T <: WeightedBagNode}
    data = filter(!ismissing, [x.data for x in as])
    metadata = filter(!isnothing, [x.metadata for x in as])
    bags = vcat((d.bags for d in as)...)
    WeightedBagNode(isempty(data) ? missing : reduce(catobs, data),
                    bags, vcat((a.weights for a in as)...),
            isempty(metadata) ? nothing : reduce(catobs, metadata))
end

removeinstances(a::WeightedBagNode, mask) = WeightedBagNode(subset(a.data, findall(mask)), adjustbags(a.bags, mask), subset(a.weights, findall(mask)), a.metadata)

Base.hash(e::WeightedBagNode{T,B,W,C}, h::UInt) where {T,B,W,C} = hash((T,B,W,C, e.data, e.bags, e.weights, e.metadata), h)
Base.:(==)(e1::WeightedBagNode{T,B,W,C}, e2::WeightedBagNode{T,B,W,C}) where {T,B,W,C} = e1.data == e2.data && e1.bags == e2.bags && e1.weights == e2.weights && e1.metadata == e2.metadata


text:
module MRFingerprintingRecon

using Polyester
using LinearAlgebra
using FFTW
using NFFT
using NFFTTools
import MRIReco
using LinearOperators

export NFFTNormalOpBasisFunc, NFFTNormalOpBasisFuncLO, calcCoilMaps, calculateBackProjection, kooshball, kooshballGA

function __init__()
    if Threads.nthreads() > 1
      BLAS.set_num_threads(1)
    end
    FFTW.set_num_threads(Threads.nthreads())
end

include("NFFTNormalOpBasisFunc.jl")
include("CoilMaps.jl")
include("BackProjection.jl")
include("Trajectories.jl")

end # module