text:
"""
   the header for regularly sampled data, the length of each string field
is 80 and support maximum of 9 dimensions.
"""
struct RegularSampleHeader{Ti<:Int64, Tv<:Float64}
    #
    n1 :: Ti
    n2 :: Ti
    n3 :: Ti
    n4 :: Ti
    n5 :: Ti
    n6 :: Ti
    n7 :: Ti
    n8 :: Ti
    n9 :: Ti
    #
    o1 :: Tv
    o2 :: Tv
    o3 :: Tv
    o4 :: Tv
    o5 :: Tv
    o6 :: Tv
    o7 :: Tv
    o8 :: Tv
    o9 :: Tv
    #
    d1 :: Tv
    d2 :: Tv
    d3 :: Tv
    d4 :: Tv
    d5 :: Tv
    d6 :: Tv
    d7 :: Tv
    d8 :: Tv
    d9 :: Tv
    #
    label1 :: String
    label2 :: String
    label3 :: String
    label4 :: String
    label5 :: String
    label6 :: String
    label7 :: String
    label8 :: String
    label9 :: String
    #
    unit1  :: String
    unit2  :: String
    unit3  :: String
    unit4  :: String
    unit5  :: String
    unit6  :: String
    unit7  :: String
    unit8  :: String
    unit9  :: String
    #
    title  :: String
    data_format :: DataType
end

"""
   the byte location of each field of the header of Uniform sampled data
"""
const field_location = Dict(
    # domension
    :n1 => 0,
    :n2 => 8,
    :n3 => 16,
    :n4 => 24,
    :n5 => 32,
    :n6 => 40,
    :n7 => 48,
    :n8 => 56,
    :n9 => 64,
    # origin
    :o1 => 72,
    :o2 => 80,
    :o3 => 88,
    :o4 => 96,
    :o5 => 104,
    :o6 => 112,
    :o7 => 120,
    :o8 => 128,
    :o9 => 136,
    # interval
    :d1 => 144,
    :d2 => 152,
    :d3 => 160,
    :d4 => 168,
    :d5 => 176,
    :d6 => 184,
    :d7 => 192,
    :d8 => 200,
    :d9 => 208,
    # label
    :label1 => 216,
    :label2 => 296,
    :label3 => 376,
    :label4 => 456,
    :label5 => 536,
    :label6 => 616,
    :label7 => 696,
    :label8 => 776,
    :label9 => 856,
    # unit
    :unit1 => 936,
    :unit2 => 1016,
    :unit3 => 1096,
    :unit4 => 1176,
    :unit5 => 1256,
    :unit6 => 1336,
    :unit7 => 1416,
    :unit8 => 1496,
    :unit9 => 1576,
    # title
    :title => 1656,
    :data_format => 1736,
    # data_format is saved as UInt8 (1 byte)
    :data  => 1737
);

"""
   constructor for the header of regularly sampled data
"""
function RegularSampleHeader(;n1=0, n2=1, n3=1, n4=1, n5=1, n6=1, n7=1, n8=1, n9=1,
         o1=0., o2=0., o3=0., o4=0., o5=0., o6=0., o7=0., o8=0., o9=0.,
         d1=0., d2=0., d3=0., d4=0., d5=0., d6=0., d7=0., d8=0., d9=0.,
         label1="", label2="", label3="", label4="", label5="", label6="", label7="", label8="", label9="",
         unit1="" , unit2="" , unit3="" , unit4="" , unit5="" , unit6="" , unit7="" , unit8="" , unit9="" ,
         title="" , data_format=Float32)

    if !(data_format <: Number)
       error("data_format must be a legal subtype of number")
    end

    hdr = RegularSampleHeader(Int64(n1)  , Int64(n2)  , Int64(n3)  , Int64(n4)  , Int64(n5)  , Int64(n6)  , Int64(n7)  , Int64(n8)  , Int64(n9)  ,
                              Float64(o1), Float64(o2), Float64(o3), Float64(o4), Float64(o5), Float64(o6), Float64(o7), Float64(o8), Float64(o9),
                              Float64(d1), Float64(d2), Float64(d3), Float64(d4), Float64(d5), Float64(d6), Float64(d7), Float64(d8), Float64(d9),
                              label1, label2, label3, label4, label5, label6, label7, label8, label9,
                              unit1 , unit2 , unit3 , unit4 , unit5 , unit6 , unit7 , unit8 , unit9 ,
                              title , data_format)

    return hdr

end

"""
   constructor for the header of regularly sampled data
"""
function RegularSampleHeader(d::Array{Tv}; n1=0, n2=1, n3=1, n4=1, n5=1, n6=1, n7=1, n8=1, n9=1,
         o1=0., o2=0., o3=0., o4=0., o5=0., o6=0., o7=0., o8=0., o9=0.,
         d1=0., d2=0., d3=0., d4=0., d5=0., d6=0., d7=0., d8=0., d9=0.,
         label1="", label2="", label3="", label4="", label5="", label6="", label7="", label8="", label9="",
         unit1="" , unit2="" , unit3="" , unit4="" , unit5="" , unit6="" , unit7="" , unit8="" , unit9="" ,
         title="" ) where {Tv<:Number}

    # get the dimensions of data
    dim  = size(d)
    ndim = length(dim)
    if ndim == 1
       n1 = dim[1]
    elseif ndim == 2
       n1 = dim[1]; n2 = dim[2];
    elseif ndim == 3
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3];
    elseif ndim == 4
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3]; n4 = dim[4];
    elseif ndim == 5
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3]; n4 = dim[4]; n5 = dim[5];
    elseif ndim == 6
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3]; n4 = dim[4]; n5 = dim[5]; n6 = dim[6];
    elseif ndim == 7
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3]; n4 = dim[4]; n5 = dim[5]; n6 = dim[6]; n7 = dim[7];
    elseif ndim == 8
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3]; n4 = dim[4]; n5 = dim[5]; n6 = dim[6]; n7 = dim[7]; n8 = dim[8];
    elseif ndim == 9
       n1 = dim[1]; n2 = dim[2]; n3 = dim[3]; n4 = dim[4]; n5 = dim[5]; n6 = dim[6]; n7 = dim[7]; n8 = dim[8]; n9 = dim[9];
    end

    # get data format
    data_format = eltype(d)

    hdr = RegularSampleHeader(n1, n2, n3, n4, n5, n6, n7, n8, n9,
                              o1, o2, o3, o4, o5, o6, o7, o8, o9,
                              d1, d2, d3, d4, d5, d6, d7, d8, d9,
                              label1, label2, label3, label4, label5, label6, label7, label8, label9,
                              unit1 , unit2 , unit3 , unit4 , unit5 , unit6 , unit7 , unit8 , unit9 ,
                              title , data_format)
    return hdr

end

"""
   construct data based on the header information
"""
function zeros(hdr::RegularSampleHeader)

    # get the dimensions of data
    if hdr.n9 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3, hdr.n4, hdr.n5, hdr.n6, hdr.n7, hdr.n8, hdr.n9)
    elseif hdr.n8 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3, hdr.n4, hdr.n5, hdr.n6, hdr.n7, hdr.n8)
    elseif hdr.n7 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3, hdr.n4, hdr.n5, hdr.n6, hdr.n7)
    elseif hdr.n6 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3, hdr.n4, hdr.n5, hdr.n6)
    elseif hdr.n5 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3, hdr.n4, hdr.n5)
    elseif hdr.n4 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3, hdr.n4)
    elseif hdr.n3 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2, hdr.n3)
    elseif hdr.n2 > 1
       return zeros(hdr.data_format, hdr.n1, hdr.n2)
    elseif hdr.n1 > 1
       return zeros(hdr.data_format, hdr.n1)
    else
       error("empty header")
    end
end

"""
   overloading the show function for the header of regularly sampled data
"""
function show(io::IO, hdr::RegularSampleHeader)

    @printf("n1= %6d, o1=%6.2f, d1=%6.2f, label1=%10s, unit1=%10s\n", hdr.n1, hdr.o1, hdr.d1, hdr.label1, hdr.unit1)
    @printf("n2= %6d, o2=%6.2f, d2=%6.2f, label2=%10s, unit2=%10s\n", hdr.n2, hdr.o2, hdr.d2, hdr.label2, hdr.unit2)
    @printf("n3= %6d, o3=%6.2f, d3=%6.2f, label3=%10s, unit3=%10s\n", hdr.n3, hdr.o3, hdr.d3, hdr.label3, hdr.unit3)
    @printf("n4= %6d, o4=%6.2f, d4=%6.2f, label4=%10s, unit4=%10s\n", hdr.n4, hdr.o4, hdr.d4, hdr.label4, hdr.unit4)
    @printf("n5= %6d, o5=%6.2f, d5=%6.2f, label5=%10s, unit5=%10s\n", hdr.n5, hdr.o5, hdr.d5, hdr.label5, hdr.unit5)
    @printf("n6= %6d, o6=%6.2f, d6=%6.2f, label6=%10s, unit6=%10s\n", hdr.n6, hdr.o6, hdr.d6, hdr.label6, hdr.unit6)
    @printf("n7= %6d, o7=%6.2f, d7=%6.2f, label7=%10s, unit7=%10s\n", hdr.n7, hdr.o7, hdr.d7, hdr.label7, hdr.unit7)
    @printf("n8= %6d, o8=%6.2f, d8=%6.2f, label8=%10s, unit8=%10s\n", hdr.n8, hdr.o8, hdr.d8, hdr.label8, hdr.unit8)
    @printf("n9= %6d, o9=%6.2f, d9=%6.2f, label9=%10s, unit9=%10s\n", hdr.n9, hdr.o9, hdr.d9, hdr.label9, hdr.unit9)

    @printf("title=%20s\n", hdr.title)
    @printf("data_format=%16s", hdr.data_format)

end

"""
   write regularly sampled data to binary file
"""
function write_RSheader(path::String, hdr::RegularSampleHeader)

    #for each field, Int64, Float64,
    fid = open(path, "w")
    for n = 1 : 27
        a = getfield(hdr, n)
        write(fid, a)
    end

    # 80 byte string
    for n = 28 : 46
        str = getfield(hdr, n)
        tmp = convert(UInt8, '\n') * ones(UInt8, 80)
        for i = 1 : length(str)
            tmp[i] = convert(UInt8, str[i])
        end
        write(fid, tmp)
    end

    #data_format: 1->Float32
                # 2->Float64
                # 3->Complex64
                # 4->Complex128
                # 5->Int32
                # 6->Int64
    if hdr.data_format == Float32
       write(fid, UInt8(1))
    elseif hdr.data_format == Float64
       write(fid, UInt8(2))
    elseif hdr.data_format == Complex64
       write(fid, UInt8(3))
    elseif hdr.data_format == Complex128
       write(fid, UInt8(4))
    elseif hdr.data_format == Int32
       write(fid, UInt8(5))
    elseif hdr.data_format == Int64
       write(fid, UInt8(6))
    end
    # the size of header is 1001 byte

    flush(fid)
    return fid

end

function write_RSdata(path::String, hdr::RegularSampleHeader, d::Array{Tv}) where {Tv<:Number}

    # check the data_format
    if hdr.data_format != eltype(d)
       error("the data format is wrong")
    end

    # check the length of data
    N = hdr.n1 * hdr.n2 * hdr.n3 * hdr.n4 * hdr.n5 * hdr.n6 * hdr.n7 * hdr.n8 * hdr.n9
    if N != length(d)
       error("the length of data is wrong")
    end

    #for each field, Int64, Float64,
    fid = open(path, "w")
    for n = 1 : 27
        a = getfield(hdr, n)
        write(fid, a)
    end

    # 80 byte string
    for n = 28 : 46
        str = getfield(hdr, n)
        tmp = convert(UInt8, '\n') * ones(UInt8, 80)
        for i = 1 : length(str)
            tmp[i] = convert(UInt8, str[i])
        end
        write(fid, tmp)
    end

    #data_format: 1->Float32
                # 2->Float64
                # 3->Complex64
                # 4->Complex128
                # 5->Int32
                # 6->Int64
    if hdr.data_format == Float32
       write(fid, UInt8(1))
    elseif hdr.data_format == Float64
       write(fid, UInt8(2))
    elseif hdr.data_format == Complex{Float32}
       write(fid, UInt8(3))
    elseif hdr.data_format == Complex{Float64}
       write(fid, UInt8(4))
    elseif hdr.data_format == Int32
       write(fid, UInt8(5))
    elseif hdr.data_format == Int64
       write(fid, UInt8(6))
    end
    # the size of header is 1001 byte

    # write the data
    write(fid, vec(d))
    close(fid)
end

"""
   read the header of regularly sampled data file
"""
function read_RSheader(path::String)

    fid = open(path, "r")

    # the size of each dimension
    n1  = read(fid, Int64)
    n2  = read(fid, Int64)
    n3  = read(fid, Int64)
    n4  = read(fid, Int64)
    n5  = read(fid, Int64)
    n6  = read(fid, Int64)
    n7  = read(fid, Int64)
    n8  = read(fid, Int64)
    n9  = read(fid, Int64)

    # start index
    o1  = read(fid, Float64)
    o2  = read(fid, Float64)
    o3  = read(fid, Float64)
    o4  = read(fid, Float64)
    o5  = read(fid, Float64)
    o6  = read(fid, Float64)
    o7  = read(fid, Float64)
    o8  = read(fid, Float64)
    o9  = read(fid, Float64)

    # interval aling each dimension
    d1  = read(fid, Float64)
    d2  = read(fid, Float64)
    d3  = read(fid, Float64)
    d4  = read(fid, Float64)
    d5  = read(fid, Float64)
    d6  = read(fid, Float64)
    d7  = read(fid, Float64)
    d8  = read(fid, Float64)
    d9  = read(fid, Float64)

    # read 11 string
    indicator = convert(UInt8, '\n')
    str = Vector{String}(undef, 19)
    a = zeros(UInt8, 80 * 19)
    read!(fid, a)
    a = reshape(a, 80, 19)
    for i = 1 : 19
        k = 1
        while a[k,i] != indicator
              k = k + 1
        end
        k = k-1
        str[i] = String(a[1:k,i])
    end

    # data format
    indicator = read(fid, UInt8)
    if indicator == 1
       data_format = Float32
    elseif indicator == 2
       data_format = Float64
    elseif indicator == 3
       data_format = Complex{Float32}
    elseif indicator == 4
       data_format = Complex{Float64}
    elseif indicator == 5
       data_format = Int32
    elseif indicator == 6
       data_format = Int64
    end

    # construct a header
    return RegularSampleHeader(n1, n2, n3, n4, n5, n6, n7, n8, n9,
                               o1, o2, o3, o4, o5, o6, o7, o8, o9,
                               d1, d2, d3, d4, d5, d6, d7, d8, d9,
                               str[1] , str[2] , str[3] , str[4] , str[5] , str[6] , str[7] , str[8] , str[9] ,
                               str[10], str[11], str[12], str[13], str[14], str[15], str[16], str[17], str[18],
                               str[19], data_format)
end

"""
   read the regularly sampled data, produce two output, one is header and the other
is the data
"""
function read_RSdata(path::String)

    fid = open(path, "r")

    # the size of each dimension
    n1  = read(fid, Int64)
    n2  = read(fid, Int64)
    n3  = read(fid, Int64)
    n4  = read(fid, Int64)
    n5  = read(fid, Int64)
    n6  = read(fid, Int64)
    n7  = read(fid, Int64)
    n8  = read(fid, Int64)
    n9  = read(fid, Int64)

    # start index
    o1  = read(fid, Float64)
    o2  = read(fid, Float64)
    o3  = read(fid, Float64)
    o4  = read(fid, Float64)
    o5  = read(fid, Float64)
    o6  = read(fid, Float64)
    o7  = read(fid, Float64)
    o8  = read(fid, Float64)
    o9  = read(fid, Float64)

    # interval aling each dimension
    d1  = read(fid, Float64)
    d2  = read(fid, Float64)
    d3  = read(fid, Float64)
    d4  = read(fid, Float64)
    d5  = read(fid, Float64)
    d6  = read(fid, Float64)
    d7  = read(fid, Float64)
    d8  = read(fid, Float64)
    d9  = read(fid, Float64)

    # read 11 string
    indicator = convert(UInt8, '\n')
    str = Vector{String}(undef, 19)
    a = zeros(UInt8, 80 * 19)
    read!(fid, a)
    a = reshape(a, 80, 19)
    for i = 1 : 19
        k = 1
        while a[k,i] != indicator
              k = k + 1
        end
        k = k-1
        str[i] = String(a[1:k,i])
    end

    # data format
    indicator = read(fid, UInt8)
    if indicator == 1
       data_format = Float32
    elseif indicator == 2
       data_format = Float64
    elseif indicator == 3
       data_format = Complex{Float32}
    elseif indicator == 4
       data_format = Complex{Float64}
    elseif indicator == 5
       data_format = Int32
    elseif indicator == 6
       data_format = Int64
    end

    # construct a header
    hdr =  RegularSampleHeader(n1, n2, n3, n4, n5, n6, n7, n8, n9,
                               o1, o2, o3, o4, o5, o6, o7, o8, o9,
                               d1, d2, d3, d4, d5, d6, d7, d8, d9,
                               str[1] , str[2] , str[3] , str[4] , str[5] , str[6] , str[7] , str[8] , str[9] ,
                               str[10], str[11], str[12], str[13], str[14], str[15], str[16], str[17], str[18],
                               str[19], data_format)

    N = n1 * n2 * n3 * n4 * n5 * n6 * n7 * n8 * n9
    data = Vector{data_format}(undef, N)

    # read data
    read!(fid, data)
    if n9 != 1
       data = reshape(data, n1, n2, n3, n4, n5, n6, n7, n8, n9)
    elseif n8 != 1
       data = reshape(data, n1, n2, n3, n4, n5, n6, n7, n8)
    elseif n7 != 1
       data = reshape(data, n1, n2, n3, n4, n5, n6, n7)
    elseif n6 != 1
       data = reshape(data, n1, n2, n3, n4, n5, n6)
    elseif n5 != 1
       data = reshape(data, n1, n2, n3, n4, n5)
    elseif n4 != 1
       data = reshape(data, n1, n2, n3, n4)
    elseif n3 != 1
       data = reshape(data, n1, n2, n3)
    elseif n2 != 1
       data = reshape(data, n1, n2)
    end
    close(fid)

    return hdr, data
end

function l2norm_rsf(path::String)
    (hdr, d) = read_RSdata(path)
    return norm(d)
end


text:
export
    AbstractCBSModel,
    DefaultCBSModel,
    PrioritizedDFSPlanner,
    AbstractISPSModel,
    DefaultISPSModel,
    PrioritizedPlannerISPSModel,
    AbstractPathFinderModel,
    AStarPathFinderModel,
    PrioritizedAStarModel,
    DFS_PathFinder,
    PC_TAPF_Solver

abstract type AbstractCBSModel end
struct DefaultCBSModel <: AbstractCBSModel end
@with_kw struct PrioritizedDFSPlanner <: AbstractCBSModel
    max_iters::Int = 3000
end
abstract type AbstractISPSModel end
@with_kw struct DefaultISPSModel <: AbstractISPSModel
    n_repair_iters::Int = 2
end
abstract type AbstractPathFinderModel end
@with_kw struct AStarPathFinderModel <: AbstractPathFinderModel
    search_history::Vector{State}   = State[]
    replan::Bool                    = false # Flag for replanning with empty conflict table after timeout
end
@with_kw struct PrioritizedAStarModel <: AbstractPathFinderModel
    search_history::Vector{State}   = State[]
    replan::Bool                    = false # Flag for replanning with empty conflict table after timeout
end
struct DFS_PathFinder <: AbstractPathFinderModel end

@with_kw mutable struct PC_TAPF_Solver{M,C,I,A,T}
    # TODO parameterize by MILP Solver, CBS solver, ISPS solver, A_star solver
    nbs_model                   ::M = SparseAdjacencyMILP()
    cbs_model                   ::C = DefaultCBSModel()
    isps_model                  ::I = DefaultISPSModel()
    astar_model                 ::A = AStarPathFinderModel()
    LIMIT_assignment_iterations ::Int = 50
    LIMIT_CBS_iterations        ::Int = 100
    LIMIT_A_star_iterations     ::Int = 5000

    num_assignment_iterations   ::Int = 0
    num_CBS_iterations          ::Int = 0
    num_A_star_iterations       ::Int = 0
    best_cost                   ::T   = get_infeasible_cost(astar_model)

    total_assignment_iterations ::Int = 0
    total_CBS_iterations        ::Int = 0
    total_A_star_iterations     ::Int = 0

    max_CBS_iterations          ::Int = 0
    max_A_star_iterations       ::Int = 0

    start_time                  ::Float64 = time()
    time_limit                  ::Float64 = 200.0
    nbs_time_limit              ::Float64 = 190.0

    verbosity                   ::Int = 0
    l1_verbosity                ::Int = 0
    l2_verbosity                ::Int = 0
    l3_verbosity                ::Int = 0
    l4_verbosity                ::Int = 0
    DEBUG                       ::Bool = false
end

function TOML.parse(solver::S) where {S<:PC_TAPF_Solver}
    toml_dict = Dict()

    toml_dict["LIMIT_assignment_iterations"] = solver.LIMIT_assignment_iterations
    toml_dict["LIMIT_CBS_iterations"] = solver.LIMIT_CBS_iterations
    toml_dict["LIMIT_A_star_iterations"] = solver.LIMIT_A_star_iterations
    toml_dict["num_assignment_iterations"] = solver.num_assignment_iterations
    toml_dict["num_CBS_iterations"] = solver.num_CBS_iterations
    toml_dict["num_A_star_iterations"] = solver.num_A_star_iterations
    if any(i->i==Inf, solver.best_cost)
        toml_dict["best_cost"] = map(i->-1,1:length(solver.best_cost))
    else
        toml_dict["best_cost"] = collect(solver.best_cost)
    end
    toml_dict["total_assignment_iterations"] = solver.total_assignment_iterations
    toml_dict["total_CBS_iterations"] = solver.total_CBS_iterations
    toml_dict["total_A_star_iterations"] = solver.total_A_star_iterations
    toml_dict["max_CBS_iterations"] = solver.max_CBS_iterations
    toml_dict["max_A_star_iterations"] = solver.max_A_star_iterations
    toml_dict["time_limit"] = solver.time_limit
    toml_dict["verbosity"] = solver.verbosity
    toml_dict["l1_verbosity"] = solver.l1_verbosity
    toml_dict["l2_verbosity"] = solver.l2_verbosity
    toml_dict["l3_verbosity"] = solver.l3_verbosity
    toml_dict["l4_verbosity"] = solver.l4_verbosity
    toml_dict["DEBUG"]        = solver.DEBUG

    toml_dict
end
function read_solver(toml_dict::Dict)
    solver = PC_TAPF_Solver(
        LIMIT_assignment_iterations = toml_dict["LIMIT_assignment_iterations"],
        LIMIT_CBS_iterations = toml_dict["LIMIT_CBS_iterations"],
        LIMIT_A_star_iterations = toml_dict["LIMIT_A_star_iterations"],

        num_assignment_iterations = toml_dict["num_assignment_iterations"],
        num_CBS_iterations = toml_dict["num_CBS_iterations"],
        num_A_star_iterations = toml_dict["num_A_star_iterations"],
        best_cost = Tuple(toml_dict["best_cost"]),

        total_assignment_iterations = toml_dict["total_assignment_iterations"],
        total_CBS_iterations = toml_dict["total_CBS_iterations"],
        total_A_star_iterations = toml_dict["total_A_star_iterations"],

        max_CBS_iterations = toml_dict["max_CBS_iterations"],
        max_A_star_iterations = toml_dict["max_A_star_iterations"],

        time_limit = toml_dict["time_limit"],

        verbosity = toml_dict["verbosity"],
        l1_verbosity = toml_dict["l1_verbosity"],
        l2_verbosity = toml_dict["l2_verbosity"],
        l3_verbosity = toml_dict["l3_verbosity"],
        l4_verbosity = toml_dict["l4_verbosity"],
        DEBUG = toml_dict["DEBUG"]
    )
end
function read_solver(io)
    read_solver(TOML.parsefile(io))
end

function CRCBS.check_termination_criteria(solver::S,env::E,cost_so_far,s) where {S<:PC_TAPF_Solver,E<:AbstractLowLevelEnv}
    # solver.num_A_star_iterations += 1
    if solver.num_A_star_iterations > solver.LIMIT_A_star_iterations
        # throw(SolverAstarMaxOutException(string("# MAX OUT: A* limit of ",solver.LIMIT_A_star_iterations," exceeded.")))
        return true
    end
    return false
end

ext: jl


text:
"""
    randomIntegerGenerator(rng::Random.AbstractRNG, genotypeSize::Integer,
                           min::Integer, max::Integer)::IntegerGenotype

Generates an integer individual with genes generated randomly within **min** and **max**, both included.

`Selfprovided Arguments` are already provided by the library, so, only `User Arguments` must be provided.

# Selfprovided Arguments
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `genotypeSize::Integer`: number of genes.
- `min::Integer`: lower limit for random generation.
- `max::Integer`: upper limit for random generation.
"""
function randomIntegerGenerator(rng::Random.AbstractRNG, genotypeSize::Integer,
                                min::Integer, max::Integer)::IntegerGenotype

    # Checks if user introduce correctly the bounds
    if min > max
        min, max = max, min
    end
    # Memory reserve
    genotype = Array{Integer}(undef, genotypeSize)
    # Range for remainder
    range = max - min + 1
    # For each gen
    for i = eachindex(genotype)
        # Generates a number within range
        genotype[i] = (rand(rng, UInt)%range)+min
    end

    return IntegerGenotype{Int64}(genotype)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
module RocketMapToOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: map_to()" begin

    println("Testing: operator map_to()")

    run_testset([
        (
            source      = from(1:5) |> map_to('a'),
            values      = @ts([ 'a', 'a', 'a', 'a', 'a', c ]),
            source_type = Char
        ),
        (
            source = timer(0, 10) |> take(3) |> map_to(1),
            values = @ts([ 1 ] ~ 10 ~ [ 1 ] ~ 10 ~ [ 1, c ])
        ),
        (
            source      = completed() |> map_to(1),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted(String, "e") |> map_to(1),
            values      = @ts(e("e")),
            source_type = Int
        ),
        (
            source = never() |> map_to(1),
            values = @ts()
        )
    ])

end

end


text:
path = @get_scratch!("planarpush")

nq = planarpush.nq
nu = planarpush.nu
nc = planarpush.nc
nz = num_var(planarpush)
nθ = 2 * planarpush.nq + planarpush.nu + planarpush.nw + 1

# Declare variables
@variables z[1:nz]
@variables θ[1:nθ]
@variables κ[1:1]

# Residual
r_pp = residual(planarpush, z, θ, κ)
rz_pp = Symbolics.jacobian(r_pp, z)
rθ_pp = Symbolics.jacobian(r_pp, θ)

# Build function
r_pp_func = build_function(r_pp, z, θ, κ)[2]
rz_pp_func = build_function(rz_pp, z, θ)[2]
rθ_pp_func = build_function(rθ_pp, z, θ)[2]
rz_pp_array = similar(rz_pp, Float64)
rθ_pp_array = similar(rθ_pp, Float64)

@save joinpath(path, "residual.jld2") r_pp_func rz_pp_func rθ_pp_func rz_pp_array rθ_pp_array
@load joinpath(path, "residual.jld2") r_pp_func rz_pp_func rθ_pp_func rz_pp_array rθ_pp_array


text:
using Pkg
Pkg.activate(".")
println("Loading SpineBasedRecordLinkage.jl")
using SpineBasedRecordLinkage
configfile = ARGS[1]
run_linkage(configfile)


text:
struct Root
    function Root()
        local_boil(@pot DefaultLogger() name=LOGGER logger=current_logger())
        local_boil(@pot Scheduler() name=SCHEDULER)
        new()
    end
end

function banner(io::IO=stdout;color=true)
    c = Base.text_colors
    tx = c[:normal] # text
    d1 = c[:bold] * c[:blue]    # first dot
    d2 = c[:bold] * c[:red]     # second dot
    d3 = c[:bold] * c[:green]   # third dot
    d4 = c[:bold] * c[:magenta] # fourth dot

    if color
        print(io,
        """
          ____        _                     |  > 是非成败转头空
         / $(d1)__$(tx) \\      | |                    |  > Success or failure,
        | $(d1)|  |$(tx) | ___ | | ___  _ __   __ _   |  > right or wrong,
        | $(d1)|  |$(tx) |/ $(d2)_$(tx) \\| |/ $(d3)_$(tx) \\| '_ \\ / $(d4)_$(tx)` |  |  > all turn out vain.
        | $(d1)|__|$(tx) | $(d2)(_)$(tx) | | $(d3)(_)$(tx) | | | | $(d4)(_)$(tx) |  |
         \\____/ \\___/|_|\\___/|_| |_|\\__, |  |  The Immortals by the River
                                     __/ |  |  -- Yang Shen
                                    |___/   |  (Translated by Xu Yuanchong)
        """)
    else
        print(io,
        """
          ____        _                     |  > 是非成败转头空
         / __ \\      | |                    |  > Success or failure,
        | |  | | ___ | | ___  _ __   __ _   |  > right or wrong,
        | |  | |/ _ \\| |/ _ \\| '_ \\ / _` |  |  > all turn out vain.
        | |__| | (_) | | (_) | | | | (_) |  |
         \\____/ \\___/|_|\\___/|_| |_|\\__, |  |  The Immortals by the River
                                     __/ |  |  -- Yang Shen
                                    |___/   |  (Translated by Xu Yuanchong)
        """)
    end
end

function start(config_file::String="Oolong.yaml";kw...)
    config = nothing
    if isfile(config_file)
        @info "Found $config_file. Loading configs..."
        config = Configurations.from_dict(Config, YAML.load_file(config_file; dicttype=Dict{String, Any});kw...)
    else
        @info "$config_file not found. Using default configs."
        config = Config(;kw...)
    end
    start(config)
end

function start(config::Config)
    config.banner && banner(color=config.color)

    @info "$(@__MODULE__) starting..."
    if myid() == 1
        local_boil(@pot Root() name=ROOT logger=current_logger())
    end

    if myid() in workers()
        local_boil(@pot LocalScheduler() name=local_scheduler())
    end
end

function stop()
end


text:
#=
LCM Server: an LCM interface to Caesar.jl

=#

using Caesar, RoME
using TransformUtils, Rotations, CoordinateTransformations
using Distributions
using PyCall, PyLCM
# using LibBSON
# using CloudGraphs # for sorryPedro

function gen_bindings()
    @show lcmtpath = joinpath(dirname(@__FILE__),"lcmtypes")
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_point_cloud_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_node_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_pose_nh_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_pose_xyh_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_prior_zpr_t.lcm`)
    println("Adding lcmtypes dir to Python path: $(lcmtpath)")
    unshift!(PyVector(pyimport("sys")["path"]),lcmtpath)
end

println("[Caesar.jl] (re)generating LCM bindings")
gen_bindings()

println("[Caesar.jl] Importing LCM message types")
@pyimport rome

## ========================================================================== ##

wTx = Vector{AffineMap}(5)
wTx[1] = Translation(8.7481, 7.51, 4.5662) ∘ LinearMap(Quat(0.9939, -0.0280019, -0.00824962, -0.106355))
wTx[2] = Translation(9.43, 11.0879, 4.9) ∘ LinearMap(Quat(0.992387, 0.0575732, -0.00698004, -0.108645))
wTx[3] = Translation(9.94, 15.1484, 4.8792) ∘ LinearMap(Quat(0.997597, 0.0563353, -0.0115003, -0.0386678))
wTx[4] = Translation(9.93104, 14.5831, 5.8) ∘ LinearMap(Quat(0.997626, 0.0336346, -0.011705, -0.0589371))
wTx[5] = Translation(9.12, 10.1615, 5.84927) ∘ LinearMap(Quat(0.998351, 0.0362801, -0.00973871, -0.0434122))

wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
wTx3 = convert(SE3, wTx[3])
wTx4 = convert(SE3, wTx[4])
wTx5 = convert(SE3, wTx[5])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R);  wEx1.Y = 0.0;
wRlx1 = SE3(zeros(3), wEx1)
wEx2 = convert(Euler, wTx2.R);  wEx2.Y = 0.0;
wRlx2 = SE3(zeros(3), wEx2)
wEx3 = convert(Euler, wTx3.R);  wEx3.Y = 0.0;
wRlx3 = SE3(zeros(3), wEx3)
wEx4 = convert(Euler, wTx4.R);  wEx4.Y = 0.0;
wRlx4 = SE3(zeros(3), wEx4)
wEx5 = convert(Euler, wTx5.R);  wEx5.Y = 0.0;
wRlx5 = SE3(zeros(3), wEx5)
# Odometries
wRlx1Tx2 = wRlx1 * (wTx1\wTx2)
wRlx2Tx3 = wRlx2 * (wTx2\wTx3)
wRlx3Tx4 = wRlx3 * (wTx3\wTx4)
wRlx4Tx5 = wRlx4 * (wTx4\wTx5)

vEx1_2 = veeEuler(wRlx1Tx2)
vEx2_3 = veeEuler(wRlx2Tx3)
vEx3_4 = veeEuler(wRlx3Tx4)
vEx4_5 = veeEuler(wRlx4Tx5)

XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
XYH2_3 = [vEx2_3[1], vEx2_3[2], vEx2_3[6]]
XYH3_4 = [vEx3_4[1], vEx3_4[2], vEx3_4[6]]
XYH4_5 = [vEx4_5[1], vEx4_5[2], vEx4_5[6]]

zrp2 = veeEuler(wTx2)[[3;4;5]]
zrp3 = veeEuler(wTx3)[[3;4;5]]
zrp4 = veeEuler(wTx4)[[3;4;5]]
zrp5 = veeEuler(wTx5)[[3;4;5]]

lcm_node = LCM()

# send the first pose
msg = rome.pose_node_t()
msg[:utime] = 0
msg[:id] = 1
msg[:mean_dim] = 7
q1 = convert(Quaternion, wTx1.R)
msg[:mean] = Float64[wTx1.t..., q1.s, q1.v...]
msg[:covar_dim] = 6
msg[:covar] = Float64[0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
publish(lcm_node, "ROME_POSES", msg)
sleep(0.01)

# send the second pose
msg = rome.pose_node_t()
msg[:utime] = 0
msg[:id] = 2
msg[:mean_dim] = 7
q2 = convert(Quaternion, wTx2.R)
msg[:mean] = Float64[wTx2.t..., q2.s, q2.v...]
msg[:covar_dim] = 6
msg[:covar] = Float64[0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
publish(lcm_node, "ROME_POSES", msg)
sleep(0.01)

msg = rome.prior_zpr_t()
msg[:utime] = 0
msg[:id] = 2
msg[:z] = zrp2[1]
msg[:roll] = zrp2[2]
msg[:pitch] = zrp2[3]
msg[:var_z] = 0.01
msg[:var_pitch] = 0.001
msg[:var_roll] = 0.001
publish(lcm_node, "ROME_PARTIAL_ZPR", msg)
sleep(0.01)

msg = rome.pose_pose_xyh_t()
msg[:utime] = 0
msg[:node_1_utime] = 0
msg[:node_1_id] = 1
msg[:node_2_utime] = 0
msg[:node_2_id] = 2
msg[:delta_x] = XYH1_2[1]
msg[:delta_y] = XYH1_2[2]
msg[:delta_yaw] = XYH1_2[3]
msg[:var_x] = 0.001
msg[:var_y] = 0.001
msg[:var_yaw] = 0.001
publish(lcm_node, "ROME_PARTIAL_XYH", msg)

# send the third pose
msg = rome.pose_node_t()
msg[:utime] = 0
msg[:id] = 3
msg[:mean_dim] = 7
q3 = convert(Quaternion, wTx3.R)
msg[:mean] = Float64[wTx3.t..., q3.s, q3.v...]
msg[:covar_dim] = 6
msg[:covar] = Float64[0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
publish(lcm_node, "ROME_POSES", msg)
sleep(0.01)

msg = rome.prior_zpr_t()
msg[:utime] = 0
msg[:id] = 3
msg[:z] = zrp3[1]
msg[:roll] = zrp3[2]
msg[:pitch] = zrp3[3]
msg[:var_z] = 0.01
msg[:var_pitch] = 0.001
msg[:var_roll] = 0.001
publish(lcm_node, "ROME_PARTIAL_ZPR", msg)
sleep(0.01)

msg = rome.pose_pose_xyh_t()
msg[:utime] = 0
msg[:node_1_utime] = 0
msg[:node_1_id] = 2
msg[:node_2_utime] = 0
msg[:node_2_id] = 3
msg[:delta_x] = XYH2_3[1]
msg[:delta_y] = XYH2_3[2]
msg[:delta_yaw] = XYH2_3[3]
msg[:var_x] = 0.001
msg[:var_y] = 0.001
msg[:var_yaw] = 0.001
publish(lcm_node, "ROME_PARTIAL_XYH", msg)

# send the forth pose
msg = rome.pose_node_t()
msg[:utime] = 0
msg[:id] = 4
msg[:mean_dim] = 7
q4 = convert(Quaternion, wTx4.R)
msg[:mean] = Float64[wTx4.t..., q4.s, q4.v...]
msg[:covar_dim] = 6
msg[:covar] = Float64[0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
publish(lcm_node, "ROME_POSES", msg)
sleep(0.01)

msg = rome.prior_zpr_t()
msg[:utime] = 0
msg[:id] = 4
msg[:z] = zrp4[1]
msg[:roll] = zrp4[2]
msg[:pitch] = zrp4[3]
msg[:var_z] = 0.01
msg[:var_pitch] = 0.001
msg[:var_roll] = 0.001
publish(lcm_node, "ROME_PARTIAL_ZPR", msg)
sleep(0.01)

msg = rome.pose_pose_xyh_t()
msg[:utime] = 0
msg[:node_1_utime] = 0
msg[:node_1_id] = 3
msg[:node_2_utime] = 0
msg[:node_2_id] = 4
msg[:delta_x] = XYH3_4[1]
msg[:delta_y] = XYH3_4[2]
msg[:delta_yaw] = XYH3_4[3]
msg[:var_x] = 0.001
msg[:var_y] = 0.001
msg[:var_yaw] = 0.001
publish(lcm_node, "ROME_PARTIAL_XYH", msg)

# send the fifth pose
msg = rome.pose_node_t()
msg[:utime] = 0
msg[:id] = 5
msg[:mean_dim] = 7
q5 = convert(Quaternion, wTx4.R)
msg[:mean] = Float64[wTx4.t..., q5.s, q5.v...]
msg[:covar_dim] = 6
msg[:covar] = Float64[0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
publish(lcm_node, "ROME_POSES", msg)
sleep(0.01)

msg = rome.prior_zpr_t()
msg[:utime] = 0
msg[:id] = 5
msg[:z] = zrp5[1]
msg[:roll] = zrp5[2]
msg[:pitch] = zrp5[3]
msg[:var_z] = 0.01
msg[:var_pitch] = 0.001
msg[:var_roll] = 0.001
publish(lcm_node, "ROME_PARTIAL_ZPR", msg)
sleep(0.01)

msg = rome.pose_pose_xyh_t()
msg[:utime] = 0
msg[:node_1_utime] = 0
msg[:node_1_id] = 4
msg[:node_2_utime] = 0
msg[:node_2_id] = 5
msg[:delta_x] = XYH4_5[1]
msg[:delta_y] = XYH4_5[2]
msg[:delta_yaw] = XYH4_5[3]
msg[:var_x] = 0.001
msg[:var_y] = 0.001
msg[:var_yaw] = 0.001
publish(lcm_node, "ROME_PARTIAL_XYH", msg)

#
# function handle_poses!(slam::SLAMWrapper,
#                        message_data)
#
#     message = rome.pose_node_t[:decode](message_data)
#     id = message[:id]
#     println("[Caesar.jl] Received pose message for x$(id)")
#
#     mean = message[:mean]
#     covar = message[:covar]
#     t = [mean[1], mean[2], mean[3]]
#     qw = mean[4]
#     qxyz = [mean[5], mean[6], mean[7]]
#     q = Quaternion(qw,qxyz) # why not a (w,x,y,z) constructor?
#     pose = SE3(t,q)
#     euler = Euler(q)
#
#     node_label = Symbol("x$(id)")
#     nothing
# end
#
# function handle_priors!(slam::SLAMWrapper,
#                          message_data)
#
#
#     message = rome.prior_zpr_t[:decode](message_data)
#     id = message[:id]
#     println("[Caesar.jl] Adding prior on RPZ to x$(id)")
#
#     @show node_label = Symbol("x$(id)")
#
#     @show message[:z]
#     @show message[:pitch]
#     @show message[:roll]
#
#     @show message[:var_z]
#     @show message[:var_pitch]
#     @show message[:var_roll]
#     nothing
# end
#
#
# function handle_partials!(slam::SLAMWrapper,
#                          message_data)
#     # add XYH factor
#
#     message = rome.pose_pose_xyh_t[:decode](message_data)
#
#     origin_id = message[:node_1_id]
#     destination_id = message[:node_2_id]
#     origin_label = Symbol("x$(origin_id)")
#     destination_label = Symbol("x$(destination_id)")
#
#     println("[Caesar.jl] Adding XYH odometry constraint betwee (x$(origin_id), x$(destination_id))")
#
#     @show message[:delta_x]
#     @show message[:delta_y]
#     @show message[:delta_yaw]
#
#     @show message[:var_x]
#     @show message[:var_y]
#     @show message[:var_yaw]
#
#     @show origin_label, destination_label
#     nothing
# end
#
#
#
# # Now test and evaluate if the lcm log is in fact correct
#
# lcm_pose_handler = (channel, message_data) -> handle_poses!(slam_client, message_data )
# lcm_odom_handler = (channel, message_data) -> handle_partials!(slam_client, message_data )
# lcm_prior_handler = (channel, message_data) -> handle_priors!(slam_client, message_data )

#


text:
#=
test_cyk:
- Julia version:
- Author: bramb
- Date: 2019-02-15
=#
using Test

@testset "CYK" begin
    using HyperGraphTools.CYK

    g = Grammar()
    r1 = Rule("S", ("A","B")) # S -> A B
    @test CYK.non_terminals(r1) == ["S", "A", "B"]
    add_rule!(g,r1)

    r2 = Rule("A", "a") # A -> a
    @test CYK.non_terminals(r2) == ["A"]
    add_rule!(g,r2)

    r3 = Rule("B", "b") # B -> b
    @test CYK.non_terminals(r3) == ["B"]
    add_rule!(g,r3)

    r4 = Rule("A", ("AA","C")) # A -> AA C
    @test CYK.non_terminals(r4) == ["A", "AA", "C"]
    add_rule!(g,r4)

    r5 = Rule("C", "c") # C -> c
    @test CYK.non_terminals(r5) == ["C"]
    add_rule!(g,r5)

    r6 = Rule("AA", "ä") # C -> c
    @test CYK.non_terminals(r6) == ["AA"]
    add_rule!(g,r6)

    it = CYK.index_triples(g)
    @test it == [(1, 2, 3), (2, 4, 5)] # (S->A B, A -> AA C)

    @test validate(["ä", "c", "b"],g)
    @test !validate(["x", "y", "z"],g)
    @test validate(["a", "b"],g)

    g = Grammar()
    r1 = Rule("S",("NP","VP"))  # 1
    r2a = Rule("NP",("DET","N"))# 2
    r2b = Rule("NP","she")      # 3
    r3a = Rule("VP",("VP","PP"))# 4
    r3b = Rule("VP",("V","NP")) # 5
    r3c = Rule("VP","eats")     # 6
    r4 = Rule("DET","a")        # 7
    r5a = Rule("N","fish")      # 8
    r5b = Rule("N","fork")      # 9
    r6 = Rule("PP",("P","NP"))  # 10
    r7 = Rule("V","eats")       # 11
    r8 = Rule("P","with")       # 12
    add_rule!(g,r1)
    add_rule!(g,r2a)
    add_rule!(g,r2b)
    add_rule!(g,r3a)
    add_rule!(g,r3b)
    add_rule!(g,r3c)
    add_rule!(g,r4)
    add_rule!(g,r5a)
    add_rule!(g,r5b)
    add_rule!(g,r6)
    add_rule!(g,r7)
    add_rule!(g,r8)
    @test g.nonTerminals == ["S", "NP", "VP", "DET", "N", "PP", "V", "P"]
    @test g.ruleIndex == Dict(
        1 => [1], # S is lhs in rule 1 (r1)
        2 => [2, 3], # NP is lhs in rules 2,3 (r2a,r2b)
        3 => [4, 5, 6], # VP is lhs in rules 4,5,6 (r3a,r3b,r3c)
        4 => [7], # DET is lhs in rule 7 (r4)
        5 => [8, 9], # N is lhs in rules 8,9 (r5a,r5b)
        6 => [10], # PP is lhs in rule 10 (r6)
        7 => [11], # V is lhs in rule 11 (r7)
        8 => [12] # P is lhs in rule 12 (r8)
    )
    @test CYK.index_triples(g) == [(1, 2, 3), (2, 4, 5), (3, 3, 6), (3, 7, 2), (6, 8, 2)]
#     println(CYK.lhs_index(g,"she"))
#     println(CYK.lhs_index(g,"eats"))
#     println(CYK.lhs_index(g,"a"))
#     println(CYK.lhs_index(g,"fish"))
#     println(CYK.lhs_index(g,"with"))
#     println(CYK.lhs_index(g,"a"))
#     println(CYK.lhs_index(g,"fork"))
    @test validate(["she", "eats", "a", "fish", "with", "a", "fork"],g)
end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: brambg/HyperGraphTools

max_forks_repo_name: brambg/HyperGraphTools


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

function knl_fillsendbuf!(::Val{Np}, ::Val{nvar}, sendbuf, buf,
                          vmapsend, nvmapsend) where {Np, nvar}

  @inbounds @loop for i in (1:nvmapsend;
                            threadIdx().x + blockDim().x * (blockIdx().x-1))
    e, n = fldmod1(vmapsend[i], Np)
    @unroll for s = 1:nvar
      sendbuf[s, i] = buf[n, s, e]
    end
  end
  nothing
end

function knl_transferrecvbuf!(::Val{Np}, ::Val{nvar}, buf, recvbuf,
                              vmaprecv, nvmaprecv) where {Np, nvar}

  @inbounds @loop for i in (1:nvmaprecv;
                            threadIdx().x + blockDim().x * (blockIdx().x-1))
    e, n = fldmod1(vmaprecv[i], Np)
    @unroll for s = 1:nvar
      buf[n, s, e] = recvbuf[s, i]
    end
  end
  nothing
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
module spectra
using FFTW

get_freq(t, fs) = fftshift(fftfreq(round(Int,fs*(maximum(t)-minimum(t)))+1, fs))
get_freq(t) = get_freq(t, 1.0/(t[2]-t[1]))

get_spectrum(v::Vector) = fftshift(fft((v.*hanning(length(v)))))
get_spectrum(v::Matrix) = fftshift(fft((v'.*hanning(size(v,2)))',2),2)
function get_spectrum(v::Matrix, tlim, ndt)
    j = tlim*ndt
    j = j[1]:j[2]
    get_spectrum(v[:,j])
end
function get_spectrum(v::Vector, tlim, ndt)
    j = tlim*ndt
    j = j[1]:j[2]
    get_spectrum(v[j])
end

export get_freq, get_spectrum
end


text:
"""
    mz_integrals_to_csv(spectra, compound_name, compounds)
Determine mz integrals for specific compound and output to csv
"""
function mz_integrals_to_csv(pathin, compound_name)
    json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
	settings_json = JSON3.read(json_string)
	main_compound_name = settings_json[:main_settings]["main_compound"]
    IS_name = settings_json[:main_settings]["internal_standard"]

    spectra, metadata_headers = batch_import(pathin, settings_json)

    # Import RT and mz info of valid compounds into DataFrame
	compounds_csv = CSV.read(joinpath(COMPOUNDS_CSV_LOCATION, "compounds.csv"), DataFrame)
    compound_to_analyse = filter(row -> row.compound == compound_name, compounds_csv) # DataFrame
    main_compound = filter(row -> row.compound == main_compound_name, compounds_csv)[1, :] # DataFrameRow
    internal_standard = filter(row -> row.compound == IS_name, compounds_csv)[1, :] # DataFrameRow

    # Create DataFrame for storing mz integrals
    mz_values = parse_data(compound_to_analyse[1, :])[1]
    mz_integrals_dataframe = DataFrame()
	for header in metadata_headers
		insertcols!(mz_integrals_dataframe, Symbol(header) => String[])
	end
	for mz in mz_values
		insertcols!(mz_integrals_dataframe, Symbol(mz) => Float32[])
	end

    # Analyse all spectra
	for i=1:length(spectra)
		println("Analysing spectrum $i...")
        spectrum = spectra[i]["MS1"]

        # Determine mz integrals, retrieve metadata, and store in DataFrame
        mz_integrals = analyse_spectrum(spectrum, compounds_csv, main_compound, internal_standard)
        sample_metadata = Vector{Any}([spectrum[header] for header in metadata_headers])
        if mz_integrals == 0
            push!(mz_integrals_dataframe, append!(sample_metadata, zeros(length(mz_values))))
        else
            mz_integrals = mz_integrals[Symbol(compound_name)]
            push!(mz_integrals_dataframe, append!(sample_metadata, [mz_integrals[Symbol(mz)] for mz in mz_values]))
        end

    end

    CSV.write(joinpath(@__DIR__, "$(compound_name)_mz integrals.csv"), mz_integrals_dataframe)
end

"""
    manual_integration(spectrum, RT_range, mz_values)
Manually integrate spectrum between RT range at chosen mz values
"""
function manual_integration(spectrum, RT_range, mz_values)
    index_left, index_right = RT_to_scans(spectrum, RT_range)
    spectrum_XIC = filter_XIC(spectrum, mz_values)

    return sum(spectrum_XIC[index_left : index_right])
end

"""
    visualize_peak_range(spectrum, compound_name, mz; secondary_ylims=0, predicted_RT=0, hidetitle=false, transparent=false)
Shows plot of peak with determined bounds and baseline. Will plot secondary zoomed in plot if peak has high intensity.
Set secondary_ylims to enforce a specific y max for secondary plot.
Set predicted_RT to enforce specific retention time, instead of using the retention time from compounds.csv
Setting hidetitle and transparent to true will not show a title and set the background to be transparent, respectively.
"""
function visualize_peak_range(spectrum, compound_name, mz; secondary_ylims=0, predicted_RT=0, hidetitle=false, transparent=false)
    PLOT_EXTRA_SCANS = 60
    H_V_LINEWIDTH = 3

    if transparent
        default(background_color=:transparent, foreground_color=:black)
    end

    json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
	settings_json = JSON3.read(json_string)
	main_compound_name = settings_json[:main_settings]["main_compound"]
    IS_name = settings_json[:main_settings]["internal_standard"]

    # Read compounds.csv
	compounds_csv = CSV.read(joinpath(COMPOUNDS_CSV_LOCATION, "compounds.csv"), DataFrame)
	filter!(row -> !(any(ismissing, (row.RT, row.mz)) || any((row.RT, row.mz) .== 0)), compounds_csv)

    # Retrieve row of compound and main compound from compounds.csv
    compound = filter(row -> row.compound == compound_name, compounds_csv)[1, :]
    main_compound = filter(row -> row.compound == main_compound_name, compounds_csv)[1, :]
    internal_standard = filter(row -> row.compound == IS_name, compounds_csv)[1, :]

    RT_modifier = determine_RT_modifier(spectrum, main_compound, internal_standard)

    # Cocaine not found
    if RT_modifier == -1
        return @error "Cocaine peak not found"
    end
    # Determine predicted RT
    if predicted_RT == 0
        predicted_RT = compound.RT * RT_modifier
    end

    spectrum_XIC = filter_XIC(spectrum, mz)
    peak_exists, max_scan = search_peak(spectrum_XIC, predicted_RT, spectrum)

    if peak_exists == false
        _x_lims = [i for i in RT_to_scans(spectrum, predicted_RT)]
        _x_lims[1] -= PLOT_EXTRA_SCANS; _x_lims[2] += PLOT_EXTRA_SCANS
        _plot = plot(spectrum_XIC, xlims=_x_lims, ylims=(0,
                maximum(spectrum_XIC[_x_lims[1]:_x_lims[2]]) * 1.1),
                        xlabel="scan number", ylabel="intensity", label=spectrum)
        default(background_color=:white, background_color_outside=:white)
        return _plot
    end

    # Determine baseline and bounds
    pre_bounds = max_scan - MAX_SCANS_PEAK_LEFT : max_scan + MAX_SCANS_PEAK_RIGHT
    overlap_RT_vals = parse_data(compound)[2] .* RT_modifier
    overlap_max_left, overlap_max_right = determine_overlap(spectrum_XIC, max_scan, overlap_RT_vals, spectrum)
    baseline = determine_baseline(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, pre_bounds)
    bounds = determine_bounds(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, baseline, pre_bounds)
    left_bound, right_bound = bounds[begin], bounds[end]

    if secondary_ylims == 0
        secondary_ylims = max(spectrum_XIC[left_bound] * 5, spectrum_XIC[right_bound] * 5, 4000)
    end
    _xlims = (left_bound - PLOT_EXTRA_SCANS, right_bound + PLOT_EXTRA_SCANS)

    _plot = plot(_xlims[1]:_xlims[2], spectrum_XIC[_xlims[1]:_xlims[2]], xlims=_xlims, ylims=(0, maximum(spectrum_XIC[left_bound : right_bound]) * 1.1),
               xlabel="scan number", ylabel="intensity", label="spectrum", legend=:topleft)
    if !hidetitle
        try
            plot!(title="Filename: $(spectrum["Filename"]) \nItem Number: $(spectrum["item_number"]) \nmz: $mz", titlefontsize=7)
        catch KeyError
            @warn "No metadata for title"
        end
    end

    vline!([left_bound, right_bound], linestyle=:dash, label="peak cutoffs", linewidth=H_V_LINEWIDTH)
    plot!(pre_bounds, baseline, linestyle=:dash, label="baseline ($(round(Int, mean(baseline))))", seriescolor=:green, linewidth=H_V_LINEWIDTH)

    if maximum(spectrum_XIC[left_bound : right_bound]) / secondary_ylims > 2.5

        plot!(twinx(), pre_bounds, baseline, xlims=_xlims, linestyle=:dot, ylims=(0, secondary_ylims), label="",
              seriescolor=:green, xaxis=:false, linewidth=1, right_margin=35mm)

        plot!(twinx(), _xlims[1]:_xlims[2], spectrum_XIC[_xlims[1]:_xlims[2]], xlims=_xlims,
        ylims=(0, secondary_ylims), linestyle=:dot, ylabel="intensity (zoomed)", label="", grid=:off,
        linewidth=1)
        # plot!(linestyle=:dot)
    end
    default(background_color=:white, background_color_outside=:white)
    return _plot
end

"""
    visualize_peak_determination_steps(spectrum, compound_name, mz; pathout="", hidetitle=false, transparent=false)
Visualizes different steps in determining bounds and baseline with plots.
Setting a pathout will save the plots to this path.
Setting hidetitle and transparent to true will not show a title and set the background to be transparent, respectively.
"""
function visualize_peak_determination_steps(spectrum, compound_name, mz; transparent=false, pathout="", hidetitle=false)
    if transparent
        default(background_color=:transparent, foreground_color=:black)
    end
    # default(linewidth=4, labelfontsize=20, tickfontsize=18)
    H_V_LINEWIDTH = 3
    plots = []

	# Read settings.json
	json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
	settings_json = JSON3.read(json_string)
	main_compound_name = settings_json[:main_settings]["main_compound"]
	IS_name = settings_json[:main_settings]["internal_standard"]

	# Import RT and mz info of valid compounds into DataFrame
	compounds_csv = CSV.read(joinpath(COMPOUNDS_CSV_LOCATION, "compounds.csv"), DataFrame)
	filter!(row -> !(any(ismissing, (row.RT, row.mz)) || any((row.RT, row.mz) .== 0)), compounds_csv)
	main_compound = filter(row -> row.compound == main_compound_name, compounds_csv)[1, :]
	internal_standard = filter(row -> row.compound == IS_name, compounds_csv)[1, :]
    compound_row = filter(row -> row.compound == compound_name, compounds_csv)[1, :]

    RT_modifier = determine_RT_modifier(spectrum, main_compound, internal_standard)
    @info "" RT_modifier
    if RT_modifier == -1
        @error "Main compound and/or internal standard not detected"
    end

    spectrum_XIC = filter_XIC(spectrum, mz)
    RT = compound_row.RT * RT_modifier
    overlap_RT_vals = parse_data(compound_row)[2]
    overlap_RT_vals .*= RT_modifier

    # --- Function search_peak --- #
    RT_range = (RT - MAX_RT_SHIFT, RT + MAX_RT_SHIFT)
	scan_range = collect(RT_to_scans(spectrum, RT_range))

    push!(plots, plot(1:length(spectrum_XIC), spectrum_XIC, legend=false, xlabel="Scan Number", ylabel="Intensity", title="First"))
    push!(plots, plot(deepcopy(plots[1]), title="Range to search max"))
    vline!(scan_range, linestyle=:dot, linewidth=H_V_LINEWIDTH)

    max_intensity, max_scan_number = findmax(spectrum_XIC[scan_range[1] : scan_range[2]])
	max_scan_number += scan_range[1] - 1
    if max_intensity <= 0
        @warn "$compound_name at m/z: $mz not detected"
        return plots
    end

    push!(plots, plot(deepcopy(plots[2]), title="Found max"))
    vline!([max_scan_number], linestyle=:dot, linewidth=H_V_LINEWIDTH, seriescolor=:black)

	# Bounds used for determining whether peak is noise or not
	left_scan = max_scan_number - round(Int, MAX_SCANS_PEAK_SEARCH / 2)
	right_scan = max_scan_number + round(Int, MAX_SCANS_PEAK_SEARCH / 2)
    push!(plots, plot(deepcopy(plots[1]), title="Range to determine noise"))
    vline!([left_scan, right_scan], linestyle=:dot, linewidth=H_V_LINEWIDTH)
    push!(plots, plot(deepcopy(plots[4]), title="Median"))
    hline!([median(spectrum_XIC[left_scan:right_scan])], linestyle=:dot, linewidth=H_V_LINEWIDTH)

	# Count number of median crosses
	spectrum_part_median = median(spectrum_XIC[left_scan:right_scan])
	median_crossings_count = 0
	for i in (left_scan + 1):right_scan
		if (spectrum_XIC[i] - spectrum_part_median) * (spectrum_XIC[i - 1] - spectrum_part_median) <= 0
			median_crossings_count += 1
		end
	end

	# Count number of mean crosses
	spectrum_part_mean = mean(spectrum_XIC[left_scan:right_scan])
	mean_crossings_count = 0
	for i in (left_scan + 1):right_scan
		if (spectrum_XIC[i] - spectrum_part_mean) * (spectrum_XIC[i - 1] - spectrum_part_mean) <= 0
			mean_crossings_count += 1
		end
	end

	# More crossings than allowed, noise instead of peak
	if median_crossings_count / length(spectrum_XIC[left_scan:right_scan]) > NOISE_MEDIAN_CROSSINGS_FRACTION ||
	   		mean_crossings_count / length(spectrum_XIC[left_scan:right_scan]) > NOISE_MEAN_CROSSINGS_FRACTION
        @warn "Noise detected instead of peak"
		return plots
    end
    # --- End of function --- #

    max_scan = max_scan_number
    pre_bounds = max_scan - MAX_SCANS_PEAK_LEFT : max_scan + MAX_SCANS_PEAK_RIGHT
    overlap_max_left, overlap_max_right = determine_overlap(spectrum_XIC, max_scan, overlap_RT_vals, spectrum)
	push!(plots, plot(deepcopy(plots[1]), title="Show overlap"))
    vline!([overlap_max_left, overlap_max_right], linestyle=:dot, linewidth=H_V_LINEWIDTH, seriescolor=:black)

    # --- Function determine_baseline --- #
    # Noise end as minimum closest to peak
	noise_end = max_scan - findmin(spectrum_XIC[max_scan:-1:max_scan - 30])[2] + 1
	spectrum_part = spectrum_XIC[noise_end - BASELINE_SCANS:noise_end]

	# Calculate median and set median as baseline
	spectrum_part_median = round(Int, median(spectrum_part))
	baseline = fill(spectrum_part_median, length(pre_bounds))

    push!(plots, plot(deepcopy(plots[1]), title="Show baseline"))
    vline!([noise_end - BASELINE_SCANS, noise_end], linestyle=:dot, linewidth=H_V_LINEWIDTH)
    plot!(pre_bounds, baseline, linestyle=:dot, linewidth=H_V_LINEWIDTH)
    # --- End of function --- #

    push!(plots, plot(deepcopy(plots[6]), title="Baseline with overlap values"))
    plot!(pre_bounds, baseline, linestyle=:dot, linewidth=H_V_LINEWIDTH)

    # --- Function determine_bounds --- #
    peak_bounds = determine_bounds(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, baseline, pre_bounds)
	pre_left_bound = overlap_max_left > 0 && overlap_max_left > pre_bounds[begin] ? overlap_max_left : pre_bounds[begin]
	pre_right_bound = overlap_max_right > 0 && overlap_max_right < pre_bounds[end] ? overlap_max_right : pre_bounds[end]
    # --- End of function --- #

    push!(plots, plot(deepcopy(plots[1]), title="Determining bounds"))
    vline!([peak_bounds[begin], peak_bounds[end]], linestyle=:dashdot, linewidth=H_V_LINEWIDTH)
    plot!([pre_bounds[begin]:pre_bounds[end]], baseline, linestyle=:dot, linewidth=H_V_LINEWIDTH)

    # Fill integral
    fill_bounds = peak_bounds[begin]:peak_bounds[end]
    baseline_bounds = (peak_bounds[1] - pre_bounds[begin] + 1):((peak_bounds[1] - pre_bounds[begin] + 1) + length(fill_bounds) - 1)
    plot_ = plot(deepcopy(plots[9]))
    plot!(plot_, fill_bounds, spectrum_XIC[fill_bounds], fillrange=baseline[baseline_bounds], fillalpha = 0.5, fillcolor = :grey, linewidth=0)
    vline!([pre_left_bound, pre_right_bound], seriescolor=:grey, linewidth=H_V_LINEWIDTH, linestyle=:dot)
    push!(plots, plot_)

    for (i, plot) in enumerate(plots)
        xlims!(plot, max_scan_number - (MAX_SCANS_PEAK_LEFT + 10), max_scan_number + (MAX_SCANS_PEAK_RIGHT + 10))
        ylims!(plot, 0, max_intensity*1.1)
        plot!(right_margin=15mm)
        if hidetitle
            plot!(title="")
        end
        if pathout != ""
            savefig(plot, joinpath(pathout, "plot_$i"))
        end
    end
    default(background_color=:white, background_color_outside=:white)
    return plots
end

"""
    function sort_spectra(spectra, sort_by)
Sorts spectra by chosen header
"""
function sort_spectra(spectra, sort_by)
    sort!(spectra, by = x -> x["MS1"][sort_by])
end

"""
    retrieve_spectrum(spectra, header, value)
Returns first spectrum found that has chosen value under chosen header
"""
function retrieve_spectrum(spectra, header, value)
    for spectrum in spectra
        if spectrum["MS1"][header] == value
            return spectrum["MS1"]
        end
    end

    println("Spectrum Not Found")
    return 0
end


text:
module js

import GR

pxwidth = 640
pxheight = 480

id_count = 0
js_running = false
checking_js = false
conditions = Dict()
l = ReentrantLock()

mutable struct JSTermWidget
    identifier::Int
    width::Int
    height::Int
    disp::Int
    visible::Bool
end

function inject_js()
  global wss, port

  wss = nothing
  _gr_js = if isfile(joinpath(ENV["GRDIR"], "lib", "gr.js"))
    _gr_js = try
      _gr_js = open(joinpath(ENV["GRDIR"], "lib", "gr.js")) do f
        _gr_js = read(f, String)
        _gr_js
      end
    catch e
      nothing
    end
    _gr_js
  end
  if _gr_js === nothing
      error(string("Unable to open '", joinpath(ENV["GRDIR"], "lib", "gr.js"), "'."))
  else
      display(HTML(string("""

      """)))
  end
end

function check_js()
    global js_running, checking_js, port
    if !checking_js
        checking_js = true
        d = Dict("text/html"=>string("""

        """))
        transient = Dict("display_id"=>"jsterm_check")
        Main.IJulia.send_ipython(Main.IJulia.publish[], Main.IJulia.msg_pub(Main.IJulia.execute_msg, "display_data", Dict("data"=>d, "metadata"=>Dict(), "transient"=>transient)))
    end
end

function JSTermWidget(id::Int, width::Int, height::Int, disp::Int)
  global id_count, js_running
  if GR.isijulia()
    id_count += 1
    if !js_running
      @async check_js()
    end
    JSTermWidget(id, width, height, disp, false)
  else
    error("JSTermWidget is only available in IJulia environments")
  end
end

wss = nothing

evthandler = Dict{Int32,Any}()
global_evthandler = nothing

function register_evthandler(f::Function, id)
  global evthandler, l
  if GR.isijulia()
    send_command(Dict("command" => "enable_events"), "cmd", id)
    lock(l) do
      evthandler[id] = f
    end
  else
    error("register_evthandler is only available in IJulia environments")
  end
end

function unregister_evthandler(id)
  global evthandler, l
  if GR.isijulia()
    if global_evthandler === nothing
      send_command(Dict("command" => "disable_events"), "cmd", id)
    end
    lock(l) do
      evthandler[id] = nothing
    end
  else
    error("unregister_evthandler is only available in IJulia environments")
  end
end

function register_evthandler(f::Function)
  global global_evthandler, l
  if GR.isijulia()
    send_command(Dict("command" => "enable_events"), "cmd", nothing)
    lock(l) do
      evthandler[id] = f
    end
  else
    error("register_evthandler is only available in IJulia environments")
  end
end

function unregister_evthandler()
  global global_evthandler, evthandler, l
  if GR.isijulia()
    send_command(Dict("command" => "disable_events"), "cmd", nothing)
    lock(l) do
      for key in keys(evthandler)
        if evthandler[key] !== nothing
          send_command(Dict("command" => "enable_events"), "cmd", key)
        end
      end
      evthandler[id] = nothing
    end
  else
    error("unregister_evthandler is only available in IJulia environments")
  end
end

function send_command(msg, msgtype, id=nothing)
  global wss, ws
  if GR.isijulia()
    if id !== nothing
      m = merge(msg, Dict("type" => msgtype, "id" => id))
    else
      m = merge(msg, Dict("type" => msgtype))
    end
    if !js_running
      conditions["sendonconnect"] = Condition()
      @async check_js()
      wait(conditions["sendonconnect"])
    end
    if ws != nothing
      try
        HTTP.write(ws, Array{UInt8}(JSON.json(m)))
      catch e
        ws = nothing
        js_running = false
      end
    end
  else
    error("send_command is only available in IJulia environments")
  end
end

function send_evt(msg, id)
  if GR.isijulia()
    send_command(msg, "evt", id)
  else
    error("send_evt is only available in IJulia environments")
  end
end

function disable_jseventhandling(id)
  if GR.isijulia()
    send_command(Dict("command" => "disable_jseventhandling"), "cmd", id)
  else
    error("disable_jseventhandling is only available in IJulia environments")
  end
end

function enable_jseventhandling(id)
  if GR.isijulia()
    send_command(Dict("command" => "enable_jseventhandling"), "cmd", id)
  else
    error("enable_jseventhandling is only available in IJulia environments")
  end
end

function disable_jseventhandling()
  if GR.isijulia()
    send_command(Dict("command" => "disable_jseventhandling"), "cmd", nothing)
  else
    error("disable_jseventhandling is only available in IJulia environments")
  end
end

function enable_jseventhandling()
  global ws
  if GR.isijulia()
    send_command(Dict("command" => "enable_jseventhandling"), "cmd", nothing)
  else
    error("enable_jseventhandling is only available in IJulia environments")
  end
end

function comm_msg_callback(msg)
  global conditions, evthandler, l
  data = msg
  if haskey(data, "type")
    if data["type"] == "save"
      d = Dict("text/html"=>string(""))
      transient = Dict("display_id"=>string("save_display_", data["display_id"]))
      Main.IJulia.send_ipython(Main.IJulia.publish[], Main.IJulia.msg_pub(Main.IJulia.execute_msg, "update_display_data", Dict("data"=>d, "metadata"=>Dict(), "transient"=>transient)))
    elseif data["type"] == "evt"
      lock(l) do
        if haskey(evthandler, data["id"]) && evthandler[data["id"]] !== nothing
          evthandler[data["id"]](data["content"])
        end
      end
    elseif data["type"] == "value"

    elseif data["type"] == "createDisplay"
      d = Dict("text/html"=>"")
      transient = Dict("display_id"=>string("save_display_", data["dispid"]))
      Main.IJulia.send_ipython(Main.IJulia.publish[], Main.IJulia.msg_pub(Main.IJulia.execute_msg, "display_data", Dict("data"=>d, "metadata"=>Dict(), "transient"=>transient)))
      display(HTML(string("
")))
    elseif data["type"] == "ack"
      notify(conditions[data["dispid"]])
    end
  end
end

function settooltip(tthtml, ttdata, id)
  send_command(Dict("command" => "settooltip", "html" => tthtml, "data" => ttdata), "cmd", id)
  return nothing
end

pluto_data = Dict()
pluto_disp = ""

function get_pluto_html()
  global pluto_data, pluto_disp, plutoisinit
  str = JSON.json(pluto_data)
  # remove leading and trailing '"'
  str = str[2:lastindex(str)]
  if plutoisinit
    return HTML(string("""

    """))
  else
    return "JSTerm (GR) not initialized. Run `GR.js.init_pluto()` at the end of a codecell"
  end
end

function jsterm_send(data::String, disp)
  global js_running, draw_end_condition, ws, conditions, pluto_data, pluto_disp
  if GR.isijulia()
    if !js_running
      conditions["sendonconnect"] = Condition()
      @async check_js()
      wait(conditions["sendonconnect"])
    end
      if ws != nothing
        try
          conditions[disp] = Condition()
          HTTP.write(ws, Array{UInt8}(JSON.json(Dict("json" => data, "type"=>"draw", "display"=>disp))))
          wait(conditions[disp])
        catch e
          ws = nothing
          js_running = false
        end
      end
  elseif GR.displayname() == "pluto"
    pluto_data = data
    pluto_disp = disp
  else
    error("jsterm_send is only available in IJulia environments and Pluto.jl notebooks")
  end
end

function get_prev_plot_id()
  send_command(Dict("value"=>"prev_id"), "inq")
end

function set_ref_id(id::Int)
  send_command(Dict("id"=>id), "set_ref_id")
end

function recv(name::Cstring, id::Int32, msg::Cstring)
    # receives string from C and sends it to JS via Comm
    global js_running, draw_end_condition

    id = string(UUIDs.uuid4());
    jsterm_send(unsafe_string(msg), id)
    return convert(Int32, 1)
end

function send(name::Cstring, id::Int32)
    # Dummy function, not in use
    return convert(Cstring, "String")
end

send_c = nothing
recv_c = nothing
init = false

function ws_cb(webs)
  global comm_msg_callback, ws, js_running, conditions, checking_js
  check = false
  while !eof(webs)
    data = readavailable(webs)
    if length(data) > 0
      data = String(data)
      if data == "inject-js"
        check = true
        inject_js()
        d = Dict("text/html"=>string(""))
        transient = Dict("display_id"=>"jsterm_check")
        Main.IJulia.send_ipython(Main.IJulia.publish[], Main.IJulia.msg_pub(Main.IJulia.execute_msg, "update_display_data", Dict("data"=>d, "metadata"=>Dict(), "transient"=>transient)))
        checking_js = false
      elseif data == "js-running"
        d = Dict("text/html"=>string(""))
        transient = Dict("display_id"=>"jsterm_check")
        Main.IJulia.send_ipython(Main.IJulia.publish[], Main.IJulia.msg_pub(Main.IJulia.execute_msg, "update_display_data", Dict("data"=>d, "metadata"=>Dict(), "transient"=>transient)))
        ws = webs
        js_running = true
        if haskey(conditions, "sendonconnect")
          notify(conditions["sendonconnect"])
        end
      else
        comm_msg_callback(JSON.parse(data))
      end
    end
  end
  if !check
    ws = nothing
    js_running = false
    @async check_js()
  end
end

plutoisinit = false

function init_pluto(jssource="https://gr-framework.org/downloads/gr-latest.js")
  global plutoisinit
  plutoisinit = true
  return HTML(string("""

  """))
end

function initjs()
    global send_c, recv_c, init, checking_js, port, connect_cond, connected
    if !init
      init = true
      send_c = @cfunction(send, Cstring, (Cstring, Int32))
      recv_c = @cfunction(recv, Int32, (Cstring, Int32, Cstring))
      @eval js begin
        import UUIDs
        import HTTP
        import Sockets
        import JSON
      end
      if GR.isijulia()
        connect_cond = Condition()
        connected = false
        ws_server_task = @async begin
          global port, connect_cond, connected
          #port, server = Sockets.listenany(8081)
          port = 8081
          #server = Sockets.listen(8081)
          #@async HTTP.listen(server=server, verbose=true) do webs
          @async HTTP.WebSockets.listen("0.0.0.0", UInt16(8081)) do webs
            ws_cb(webs)
            #HTTP.WebSockets.upgrade(ws_cb, webs)
          end
          connected = true
          notify(connect_cond)
        end
        if !connected
          wait(connect_cond)
        end
      end
    end
    send_c, recv_c
end

end # module


text:
"""Validates the contents of a Bus."""
function validate(bus::Bus)::Bool
    is_valid = true

    @debug "Bus validation" bus.name bus.number is_valid
    return is_valid
end


text:
#=
  Euler #15 in Julia.

  Problem 15
  """
  Starting in the top left corner of a 2×2 grid, there are 6 routes
  (without backtracking) to the bottom right corner.

  How many routes are there through a 20×20 grid?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

function prodlist(from,to)
    return prod(from:to)
end

# 0.01229874s
function euler15a()
    return prodlist(BigInt(21),40) / prodlist(2,20);
end

run_euler(euler15a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
"""
Abstract type for struct (functor) that can
block the current task for a specified delay.
"""
abstract type DelayFunc end

"""
    SleepFuncStruct()

A struct that can block (when called)
the current task for a specified number of seconds.

`SleepFuncStruct` implements `DelayFunc` abstract type.

The minimum sleep time is 1 millisecond or input of 0.001.

`SleepFunc` is an instance of `SleepFuncStruct`

    SleepFunc(duration)

is same as `sleep(duration)`
"""
struct SleepFuncStruct <: DelayFunc
    func
    SleepFuncStruct() = new(sleep)
end
function (delayfunc::SleepFuncStruct)(args...)
    delayfunc.func(args...)
end
SleepFunc = SleepFuncStruct()

"""
    NoSleepStruct()

`NoSleepStruct` implements `DelayFunc` abstract type.

`NoSleep` is an instance of `NoSleepStruct`

When called `NoSleep`, doesn't do anything.

This is simply a workaround for simulating a scheduler waiting
for next job to be processed.
"""
struct NoSleepStruct <: DelayFunc
end
function (delayfunc::NoSleepStruct)(args...)
end
NoSleep = NoSleepStruct()

"""
Default sleep function
"""
_sleep = SleepFunc


text:
### 5 point calculation

function SparseArray(t::Type,f::Function,indicies::Tuple{Int,Int})
    n,m = indicies

    V = Array{t}(0)
    I = Array{Int}(0)
    J = Array{Int}(0)

    for i in 1:n
        for j in 1:m
            val = f(i,j)
            if val!=0
                push!(I,i)
                push!(J,j)
                push!(V,val)
            end
        end
    end

    return sparse(I,J,V)
end

function d1(i,j,y)
    x(p) = y(i+p) - y(i)

    j==i-2 ? -((x(-1)*x(1)*x(2))/(x(-2)*(x(-2)-x(-1))*(x(-2)-x(1))*(x(-2)-x(2)))) :
        j==i-1 ? (x(-2)*x(1)*x(2))/((x(-2)-x(-1))*x(-1)*(x(-1)-x(1))*(x(-1)-x(2))) :
        j==i ? -(1/x(-1))-1/x(1)-1/x(2)-1/x(-2) :
        j==i+1 ? (x(-2)*x(-1)*x(2))/((x(-2)-x(1))*x(1)*(x(1)-x(-1))*(x(1)-x(2))) :
        j==i+2 ? (x(-2)*x(-1)*x(1))/((x(-2)-x(2))*x(2)*(x(2)-x(-1))*(x(2)-x(1))) :
        0
end

function d2(i,j,y)
    x(p) = y(i+p) - y(i)

    j==i-2 ? (2*(x(1)*x(2)+x(-1)*(x(1)+x(2))))/(x(-2)*(x(-2)-x(-1))*(x(-2)-x(1))*(x(-2)-x(2))) :
        j==i-1 ? -((2*(x(1)*x(2)+x(-2)*(x(1)+x(2))))/((x(-2)-x(-1))*x(-1)*(x(-1)-x(1))*(x(-1)-x(2)))) :
        j==i ? (2*(x(1)*x(2)+x(-1)*(x(1)+x(2))+x(-2)*(x(-1)+x(1)+x(2))))/(x(-2)*x(-1)*x(1)*x(2)) :
        j==i+1 ? -((2*(x(-1)*x(2)+x(-2)*(x(-1)+x(2))))/((x(-2)-x(1))*x(1)*(x(1)-x(-1))*(x(1)-x(2)))) :
        j==i+2 ? -((2*(x(-1)*x(1)+x(-2)*(x(-1)+x(1))))/((x(-2)-x(2))*x(2)*(x(2)-x(-1))*(x(2)-x(1)))) :
        0
end

### There is a bug if I do put d1,d2 inside a function

function Hamiltonian(α,II,G,p0,L,nn)

    sigmax = [0 1;1 0]
    sigmay = [0 -1im;1im 0]
    sigmaz = [1 0;0 -1]

    Ix,Iy = II
    Gx,Gy = G
    p0x,p0y = p0

    ### One could make theese nicer
    xx(i) = i<= nn//2 ? -(L+1)^((nn/2-i+1)/(nn/2+1)) + 1 : (L+1)^((i-nn/2)/(nn/2+1)) - 1
    x(i) = i<= nn//2 ? -(L+1)^((nn/2-i+1)/(nn/2+1)) + 1 - xx(nn/2)/2 : (L+1)^((i-nn/2)/(nn/2+1)) - 1 - xx(nn/2+1)/2

    ### Laplacian operator

    D2 = SparseArray(Float64,(i,j)->d2(i,j,x),(nn,nn))
    Lap = Gx*kron(D2,speye(nn)) + Gy*kron(speye(nn),D2)

    ### Potential

    V = Float64[]
    for i in 1:nn
        for j in 1:nn
            x_ = x(i)
            y_ = x(j)
            r = sqrt(x_^2+y_^2)
            push!(V,-α/r)
        end
    end
    V = spdiagm(V)

    ### Weyl coupling

    D1 = SparseArray(Float64,(i,j)->d1(i,j,x),(nn,nn))

    Hm = kron(-1/2*Lap + V,eye(2)) + Ix*kron(1im*kron(D1,speye(nn)) + p0x*speye(nn^2),sigmax) + Iy*kron(1im*kron(speye(nn),D1) + p0y*speye(nn^2),sigmay)
end

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc


text:
using MPOmodule2
using finiteMPS
using LinearAlgebra
using TensorOperations

Id(n) = Array{ComplexF64}(I, n, n)

sz = localOp([1.0 0.0; 0.0 -1.0])
sx = localOp([0.0 1.0; 1.0 0.0])
sy = localOp([0.0 -1.0im; 1.0im 0.0])
id = localOp(Id(2))

cr = localOp([0.0 1.0; 0.0 0.0])
ca = localOp([0.0 0.0; 1.0 0.0])
n = localOp([0.0 0.0; 0.0 1.0])

function startMPS(N::Int, d::Int, maxD::Int)
	MPSvec = []
	Cvec = []

	N_half = round(N/2, RoundDown)
	for i = 1:N_half

		α_dim = d^(i-1) > maxD ? maxD : Int(d^(i-1))
		β_dim = d^i > maxD ? maxD : Int(d^i)
		push!(MPSvec, randn(ComplexF64, (α_dim, d, β_dim)))
		push!(Cvec, randn(ComplexF64, (β_dim, β_dim)))
	end

	if N-2*N_half != 0
		bDim = d^N_half > maxD ? maxD : Int(d^N_half)
		push!(MPSvec, randn(ComplexF64, (bDim, d, bDim)))
		push!(Cvec, randn(ComplexF64, (bDim, bDim)))
	end

	for i = N_half:-1:1

		β_dim = d^(i-1) > maxD ? maxD : Int(d^(i-1))
		α_dim = d^i > maxD ? maxD : Int(d^i)

		push!(MPSvec, randn(ComplexF64, (α_dim, d, β_dim)))
		push!(Cvec, randn(ComplexF64, (β_dim, β_dim)))

	end

	return MPSvec, Cvec
end

function TDVP(MPSvec::Vector{<:Any}, MPOvec::Vector{MPO}, Tsteps::Int, dτ::Union{Float64, ComplexF64})
	#the inital
	leftCanMPS(MPSvec)
	Cpre, R, Cvec = rightCanMPS(MPSvec)

	Lenv = [Array{ComplexF64,2}(I,1,1)]
	@show  Lenv = applyTM_MPO(MPSvec[1:end], MPOvec[1:end], Lenv; left = false)

	@show size(Cvec)
	#
	# build up right environment blocks
	#
	Renv = [Array{ComplexF64,2}(I,1,1)]
	RBlocks = [Renv]
	for site = 1:length(MPSvec)-1
		Renv = applyTM_MPO([MPSvec[end-(site-1)]],[MPOvec[end-(site-1)]], RBlocks[site]; left = false)
		push!(RBlocks, Renv)#[init,N,.....,2]
	end

	measureMPS = deepcopy(MPSvec)
	Cnext, Remainder, __ = leftCanMPS(measureMPS, 1, Int(length(MPSvec)/2)-1)

	measurelist = []

	for i in 1:Tsteps
		println("\n\n")
		@show i
		Cpre, RBlocks, Cvec = evo_sweep(MPSvec, Cpre, MPOvec, RBlocks, dτ)

		Lenv = [Array{ComplexF64,2}(I,1,1)]
		expSz= applyTM_MPO(deepcopy(MPSvec[1:end]), MPOvecSz[1:end], Lenv; left = false)
		@show expSz./2
		push!(measurelist, (expSz./2)[1,1])

	end

	@show measurelist
	return measurelist

end

""" transver Ising """

# creat MPO representation for transver ising
#
#	 1  	0	0
#	-J*sx	λ	0
#	-h*sz	sx	1
#

J = 1.0
h = -0.5
N = 10

# for now λ = 0
Op_vec_TIsing(J, h) = Vector{localOp}([id, sx, h*sz, J*sx, id])
Op_idx_TIsing = Vector{Int}([1, 2, 3, 6, 9])
Op_vec_TIsingStart(J, h) = Vector{localOp}([h*sz, J*sx, id])
Op_vec_TIsingEnd(J, h) = Vector{localOp}([id, sx, h*sz])

TIsing = MPO(3,Op_vec_TIsing(J, h),Op_idx_TIsing)

MPOvecTIsing = [MPO((1,3),Op_vec_TIsingStart(J,h),[1,2,3]), [TIsing for i = 1:N-2]...,MPO((3,1), Op_vec_TIsingEnd(J,h),[1,2,3])]

Op_vec_Id = Vector{localOp}([id, id, id])
Op_vec_Sz = Vector{localOp}([id, sz, id])
Op_idx_Sz = Vector{Int}([1,2,4])

SzMPO = MPO(2, Op_vec_Sz, Op_idx_Sz)

MPOvecSz = [MPO((1,2),Vector{localOp}([sz, id]), Vector{Int}([1,2])),[SzMPO for i = 1:N-2]..., MPO((2,1),Vector{localOp}([id, sz]),Vector{Int}([1,2]))]

"""****************************** init MPS ********************************"""
MPSvec, C = startMPS(N,2,4)

A = zeros(ComplexF64, 20, 2, 20)
A[1, 1, 1] = 1.0
B = zeros(ComplexF64, 1, 2, 20)
B[1, 1, 1] = 1.0
D = zeros(ComplexF64, 20, 2, 1)
D[1, 1, 1] = 1.0

MPSvec[1] = B
C[1] = Id(2)
for i = 2:1:N-1
	MPSvec[i] = A
	C[i] = Id(2)
end
MPSvec[end] = D
C[end] = Id(1)

rightCanMPS(MPSvec)
leftCanMPS(MPSvec)

Lenv = [Array{ComplexF64,2}(I,1,1)]
Lenv = applyTM_MPO(MPSvec[2:end], MPOvecSz[2:end], Lenv; left = false)
@show  "exp right" , applyTM_MPO([MPSvec[1]], [MPOvecSz[1]], Lenv; left = false)

Lenv = [Array{ComplexF64,2}(I,1,1)]
Lenv = applyTM_MPO(MPSvec[1:end-1], MPOvecSz[1:end-1], Lenv; left = true)
@show  "exp left", applyTM_MPO([MPSvec[end]], [MPOvecSz[end]], Lenv; left = true)

rightCanMPS(MPSvec)
measureMPS = deepcopy(MPSvec)
Cnext, Remainder, __ = leftCanMPS(measureMPS, 1, Int(length(MPSvec)/2)-1)
@tensor C[α, β] := Cnext[α, γ]*Remainder[γ, β]

@tensor MPS[α, d, β] := C[α, γ]*MPSvec[Int(length(MPSvec)/2)][γ, d, β]

@tensor expvalue[] := MPS[α, d, β]*Array(sz)[d, d']*conj(MPS)[α, d', β]
@show "here ",expvalue

"""****************************** init Ham ********************************"""

measurelist = TDVP(MPSvec, MPOvecTIsing, 1000, 0.01im)

measurelist =[x[1,1] for x in measurelist]

#using PyCall, PyPlot
#pygui(:tk)
#PyPlot.plot(measurelist)
#PyPlot.show()


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{Any, Any}(
"physics" => "Euler",
"operator_type" => "SBPOmega",
"dimensions" => 2,
"run_type" => 5,
"jac_method" => 2,
"jac_type" => 2,
"order" => 1,
"IC_name" => "ICIsentropicVortex",
"use_DG" => true,
"Flux_name" => "RoeFlux",
"numBC" => 3,
"BC1" => [0],
"BC1_name" => "isentropicVortexBC",  # outlet
"BC2" => [2],
"BC2_name" => "isentropicVortexBC", # inlet
"BC3" => [1, 3],
"BC3_name" => "noPenetrationBC",  # was noPenetrationBC
"aoa" => 0.0,
#"SRCname" => "SRCp3",
#"delta_t" => 0.005,
"smb_name" => "SRCMESHES/vortex_3x3_.smb",
#"smb_name" => "small.smb",
"dmg_name" => ".null",
"itermax" => 20,
#"print_cond" => true,
#"write_eigs" => true,
"res_abstol" => 1e-9,
"res_reltol" => 1e-9,
#"step_tol" => 1e-10,
"write_vis" => true,
"write_adjoint_vis" => true,
"output_freq" => 1,
"do_postproc" => true,
"exact_soln_func" => "ICIsentropicVortex",
"force_solution_complex" => true,
"force_mesh_complex" => true,
"need_adjoint" => true,
)


text:
module DivConformingFESpacesTests

using Test
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Io

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

order = 1

u(x) = x

reffe = ReferenceFE(raviart_thomas,order)

V = TestFESpace(model,reffe,dirichlet_tags = [1,6])
test_single_field_fe_space(V)

U = TrialFESpace(V,u)

uh = interpolate(u,U)

e = u - uh

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

#using Gridap.Visualization
#
#writevtk(Ω,"trian",nsubcells=10,cellfields=["uh"=>uh])

order = 1

reffe = ReferenceFE(TET,raviart_thomas,order)

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = simplexify(CartesianDiscreteModel(domain,partition))

labels = get_face_labeling(model)
dir_tags = Array{Integer}(undef,0)

V = FESpace(model,reffe,conformity=DivConformity())

v(x) = VectorValue(-0.5*x[1]+1.0,-0.5*x[2],-0.5*x[3])
vh = interpolate(v,V)

e = v - vh

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

#using Gridap.Visualization
#
#writevtk(trian,"test",order=3,cellfields=["vh"=>vh])

# Tests on manifold

# Create domain
domain = (0,1,0,1,0,1)
cells  = (2,2,2)
model  = CartesianDiscreteModel(domain,cells)

# Restrict model to cube surface (using new BoundaryDiscreteModel)
labels = get_face_labeling(model)
bgface_to_mask = get_face_mask(labels,"boundary",2)
Γface_to_bgface = findall(bgface_to_mask)
Dc2Dp3model = BoundaryDiscreteModel(Polytope{2},model,Γface_to_bgface)

order  = 0
degree = 1

reffe_rt = ReferenceFE(raviart_thomas,Float64,order)
V  = FESpace(Dc2Dp3model, reffe_rt ; conformity=:HDiv)
uh = FEFunction(V,rand(num_free_dofs(V)))
vh = interpolate_everywhere(uh,V)

Tₕ_K = Triangulation(Dc2Dp3model)
Qₕ_K = CellQuadrature(Tₕ_K,degree)
e=sqrt(sum(∫((uh-vh)⋅(uh-vh))Qₕ_K))
@test e < 1.0e-12

end # module


text:
cd("INSERT_PATH")

using Distributions
using Polynomials
using DataFrames

include("toolBox.jl")
include("fqBackPredGSB.jl")

#########################################################
# Parameter Declaration
#########################################################

# map parameters
theta, x0, y0 = [1.38, 0., 0., 0., -1., 0.211], 1.5, 0.5

# sample size
ss = 250
# past prediction horizon
T = 2

# noise f_{2,l} precisions
w1, lam1 = 0.9, 1e07
lam2 = 1e02

xdet = zeros(ss+T)
xdet[1] = henon(theta, x0, y0)
xdet[2] = henon(theta, xdet[1], x0)
for i in 3:(ss+T)
  xdet[i] = henon(theta, xdet[i-1], xdet[i-2])
end

dataSeed = 50
fulldata = copy(xdet) # no noise
# fulldata = genDataf1(ss + T, theta, x0, y0, lam2, dataSeed) #f_{1} noise
# fulldata = genData(ss + T, theta, x0, y0, w1, lam1, lam2, dataSeed) #f_{2,l} noise
predValues = fulldata[1:T]
data = fulldata[(T+1):end]

# GSBR reconstruction - bPrediction
gibbsIter, burnIn, gEps1, gEps2, thLow, thUp, zpLow, zpUp, zLow, zUp, papr, pbpr, thin, samplerSeed =
 50000, 25000, 1e-03, 1e-03, -10.0, 10.0, -2.0, 2.0, -10.0, 10.0, 10., 10., 5, 12345;

filename = "/Results"
savelocation = string(pwd(), filename, "/seed$samplerSeed")
mkpath(savelocation)
writedlm(string(savelocation,"/predValues.txt"), predValues)
writedlm(string(savelocation,"/data.txt"), data, '\n')

@time f = open("PATH/noninv.txt", "a") do f
	for i in 1:25
		println("#######################################")
	    println("Simulation: $i/25")
		println("#######################################")
		x0s = fqBackPredGSB(data[i:end], T, gibbsIter, burnIn, gEps1, gEps2, thLow, thUp, zpLow, zpUp, zLow, zUp, papr, pbpr, thin, samplerSeed);
		writedlm(f, x0s)
		end
	x0s = 0
	gc()
end

max_stars_repo_name: KosKal/SM-GSBR

max_issues_repo_name: KosKal/SM-GSBR

max_forks_repo_name: KosKal/SM-GSBR


text:
using TimeSeries, HDF5, JLD, Requests

include(Pkg.dir("Quandl/src/timearray.jl"))

md  = timearray(load(Pkg.dir("Quandl/test/response.jld"))["missing"])
taa = timearray(load(Pkg.dir("Quandl/test/response.jld"))["asc"])
tad = timearray(load(Pkg.dir("Quandl/test/response.jld"))["des"])

facts("timearray works on Request object") do

  context("there are three rows") do
      @fact  length(taa) => 3
      @fact  length(tad) => 3
  end

  context("oldest date first regardless of ordering argument") do
      @fact  taa[1].timestamp[1] < taa[3].timestamp[1] => true
      @fact  tad[1].timestamp[1] < tad[3].timestamp[1] => true
  end

  context("correct value at first row, first column, regardless of ordering argument") do
      @fact  taa[1][1].values[1] => 103.02
      @fact  tad[1][1].values[1] => 103.02
  end
end

facts("floats and NaNs present") do

  context("NaN fills in missing values slot") do
      @fact  isnan(sum(md.values[:,1])) => true
  end

  context("existing values remain floats") do
      @fact  sum(md.values[:,2]) => 66.86
  end
end


text:
using Publish
using ParameterSchedulers

p = Publish.Project(ParameterSchedulers)

# needed to prevent error when overwriting
rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)

# build documentation
deploy(ParameterSchedulers; root = "/ParameterSchedulers.jl", force = true, label = "dev")


text:
function regression_plot(;max=10)
    @printf("...regression plot ...\n")
    l_idt = zeros(minimum([max, n_exp]), 3)
    for i_exp in 1:minimum([max, n_exp])
        local T0 = conds[i_exp, 1]
        local P = conds[i_exp, 2] * one_atm
        local phi = conds[i_exp, 3]
        local idt0 = conds[i_exp, 4]
        idt = f_idt(T0, P, phi, p; dT=dT)
        if i_exp < max
            check_sol(T0, P, phi, p; i_exp=i_exp)
        end
        @printf("%d idt %.2e idt0 %.2e \n", i_exp, idt, idt0)
        _idt0 = f_idt(T0, P, phi, p .* 0.0; dT=dT)
        l_idt[i_exp, :] .= [idt0, idt, _idt0]
    end
    plt = Plots.scatter(l_idt[:, 1], l_idt[:, 2], xscale=:log10, yscale=:log10, label="train-ref")
    Plots.scatter!(plt, l_idt[:, 1], l_idt[:, 3], xscale=:log10, yscale=:log10, label="sk-ref")
    plot!(plt, [minimum(l_idt), maximum(l_idt)], [minimum(l_idt), maximum(l_idt)], label="y=x")
    xlabel!(plt, "ref [log-s]")
    ylabel!(plt, "train [log-s]")
    plot!(plt, legend=:topleft)
    png(plt, string(fig_path, "/regression"))
end

function check_sol(T0, P, phi, p; i_exp=0)
    l_plt = []
    ts, pred = get_idt(T0, P, phi, p; dT=dT)
    ts .+= 1.e-8
    plt = plot(ts, pred[end, :], lw=2, label="Train")
    ylabel!(plt, "Temperature [K]")
    xlabel!(plt, "Time [s]")
    title!(plt, @sprintf("%s, IDT=%.2e [s] \n @%.1f K, %.1f atm, phi=%.1f",
                fuel, ts[end], T0, P / one_atm, phi))
    push!(l_plt, plt)
    for s in [fuel, oxygen]
        plt =
            plot(ts,
            pred[species_index(gas, "$s"), :],
            lw=2, label="Train")
        ylabel!(plt, "Y $s")
        xlabel!(plt, "Time [s]")
        push!(l_plt, plt)
    end
    pltsum = plot(l_plt..., legend=false, framestyle=:box, xscale=:log10)

    png(pltsum, string(fig_path, "/conditions/sol_$i_exp"))
end

function plot_loss()
    plt_loss = plot(l_loss, yscale=:log10, label="train");
    plot!(plt_loss, l_loss_val, yscale=:log10, label="val");
    plt_grad = plot(l_grad, yscale=:log10, label="grad_norm");
    plt_pnorm = plot(l_pnorm .+ 1.e-8, yscale=:identity, label="p_norm");
    xlabel!(plt_loss, "Epoch");
    xlabel!(plt_grad, "Epoch");
    xlabel!(plt_pnorm, "Epoch");
    ylabel!(plt_loss, "Loss");
    ylabel!(plt_grad, "Grad Norm");
    ylabel!(plt_pnorm, "p Norm");
    plt_all = plot([plt_loss, plt_grad, plt_pnorm]..., legend=:bottomleft);
    png(plt_all, string(fig_path, "/loss_grad"));
end

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction


text:
module Reinforcement

	using MLJ: categorical, mcc
	using Flux

	"""
		generate_score_fitness(x::Vector, y::Vector)

	Generates a fitness function that increases the fitness score depending on whether the
	model finds the right answer.

	Example: y = 0, predicted_y = 0, model_output = model(example), then
				score += model_output[predicted_y+1]
			#takes the value of the respective output node

			 y = 0, predicted_y = 1, model_output = model(example), then
			 	score = score
			#does nothing
	"""
	function generate_score_fitness(x::Vector, y::Vector)
		function fitness(chain::Chain)
			score = 0.0
			for (xi, yi) in zip(x, y)
				out = chain(xi)
				predicted = argmax(out) - 1
				if predicted == yi
					if yi == 1
						score += out[2]
					else
						score += out[1]
					end
				end
			end
			score
		end
	end

	"""
		generate_mcc_fitness(x::Vector, y::Vector)

	Generates a fitness function that increases the fitness score depending on the
	Matthews Correlation Coefficient (MCC).
	"""
	function generate_mcc_fitness(x::Vector, y::Vector)
		predicted = categorical(zeros(Int, length(y)))
		caty = categorical(y)
		function fitness(chain::Chain)
			for (xi, i) in zip(x, eachindex(y))
				out = chain(xi)
				predicted[i] = argmax(out) - 1
			end
			mcc(caty, predicted)
		end
	end

	export generate_score_fitness, generate_mcc_fitness

end


text:
#testing seed is preserved
Fs = 4e4                    # Sampling frequency
c = 343                     # create new acoustic env with default values
Nt = round(Int64,0.1*Fs)    # Number of time samples, 0.25 secs
L  = (4., 4.,  4.)          # Room dimensions
Lx,Ly,Lz = L
xs = (2., 1.5, 1.)      # Source position
xr = (1., 2.,  2.)   # Receiver position

T60 = 0.1           # Reverberation Time
beta = ImageMethodReverb.revTime2beta(L, T60, c)
Tw = 20              # samples of Low pass filter
Fc = 0.9             # cut-off frequency

println("testing β for a given T60 gives correct T60")
h, seed  = rim(xs,xr,L,beta,Nt,Fs; N = (4,4,4))
rt, edc = revTime(h, Fs)
@test norm(rt.-T60)/norm(T60)< 0.25 #testing rev time

println("testing seed is preserved")
h2,  = rim(xs,xr,L,beta,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(h-h2) < 1e-8

println("testing calls")
# with T60
h2,  = rim(xs,xr,L,T60,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(h-h2) < 1e-8

## original image source method
println("test with no fractional delay")
xs = (2., 1.5, 1.)      # Source position
xr = (1., 2.,  2.)   # Receiver position
h2, = rim(xs,xr,L,beta,Nt,Fs; Tw = 0, Fc = 0., N = (4,4,4), seed = seed)

println("test line of sight")
xs = (1., 1., 1.)
for i = 1:10
  d = (0.2 .*randn(3)...,) # dandom dispacement
  h_rand, = rim(xs, xs .+ d, L, 0.001, Nt,Fs; Rd = 1e-1, N = (0,0,0))
  @test argmax(h_rand)*1/Fs - norm(xr .- xs)/343 < 1e-4
end

# testing errors
L, xr, xs = (1,1,1), (0.5,0.5,0.5),(0.7,0.7,0.7)
@test_throws ErrorException rim((0.1,), xr, L, beta, Nt,Fs)
@test_throws ErrorException rim((1.1,0,0), xr, L, beta, Nt,Fs)
@test_throws ErrorException rim(xs,(0.1,), L, beta, Nt,Fs)
@test_throws ErrorException rim(xs,(1.1,0,0), L, beta, Nt,Fs)
@test_throws ErrorException rim(xs,xr,(0,0), beta, Nt,Fs)
@test_throws ErrorException rim(xs,xr,L, beta, Nt,Fs; N=(-1,1,1))
@test_throws ErrorException rim(xs,xr,L, beta, Nt,Fs; N=(1,1))
@test_throws ErrorException rim(xs,xr,L, beta, Nt,Fs; N=(1,1))
@test_throws ErrorException rim(xs,xr,L, beta, Nt,Fs; N=(1.5,1,1))
@test_throws ErrorException rim(xs,xr,L, (1,1), Nt,Fs)

max_stars_repo_name: nantonel/JuRIM

max_issues_repo_name: nantonel/JuRIM

max_forks_repo_name: nantonel/JuRIM