text:
## Wrapping hybrd routine
const _hybrd_messages = Dict{Int, String}(
    0 => "improper input parameters",
    1 => string("algorithm estimates that the relative error between x and the ",
                "solution is at most tol"),
    2 => "maximum iterations has been exceeded",
    3 => "tol is too small, no further improvement in x is possible",
    4 => string("iteration is not making good progress, measured by ",
                "improvement from last 5 jacobian evaluations"),
    5 => string("iteration is not making good progress, measured by ",
                "improvement from last 10 iterations"),
    -1 => "exceeded user imposed number of iterations",
    -2 => "user terminated iterations with code "
)

function hybrd(f!::Function, x0::Vector{Float64}, tol::Float64,
               show_trace::Bool, tracing::Bool, maxit::Int;
               _n::Int=length(x0), ml::Int=_n-1, mu::Int=_n-1,
               epsfcn::Float64=0.0, diag::Vector{Float64}=fill(1.0, _n),
               mode::Int=2, factor::Float64=100.0, nprint::Int=0,
               lr::Cint=ceil(Cint, _n*(_n+1)/2))
    n = length(x0)
    x = copy(x0)
    fvec = similar(x)
    fjac = Array{Float64}(undef, n, n)
    lwa = ceil(Int, (n*(3*n+13))/2)
    wa = ones(lwa)
    trace = AlgoTrace(x0, show_trace, tracing, maxit)

    r = Array{Float64}(undef, lr)
    qtf = Array{Float64}(undef, n)
    wa1 = Array{Float64}(undef, n)
    wa2 = Array{Float64}(undef, n)
    wa3 = Array{Float64}(undef, n)
    wa4 = Array{Float64}(undef, n)

    if show_trace
        show(trace)
    end

    function _hybrd_func_wrapper(_p::Ptr{Cvoid}, n::Cint, _x::Ptr{Cdouble},
                                 _fvec::Ptr{Cdouble}, iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, n)
        local x = unsafe_wrap(Array, _x, n)
        if iflag < 0
            print(fvec)
            return Cint(0)
        end
        f!(fvec, x)

        local trace = unsafe_pointer_to_objref(_p)::AlgoTrace
        push!(trace, x, fvec, iflag)

        trace.f_calls > trace.maxit ? Cint(-1) : Cint(0)
    end
    _hybrd_cfunc = @cfunction(
        $_hybrd_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint)
    )

    return_code = ccall(
        (:hybrd, cminpack),
        Cint,
        (
            Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble,
            Cint, Cint, Cint, Cdouble, Ptr{Cdouble}, Cint, Cdouble, Cint,
            Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble},
            Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}
        ),
        _hybrd_cfunc, pointer_from_objref(trace), n, x, fvec, tol, typemax(Cint),
        ml, mu, epsfcn, diag, mode, factor, nprint, [0], fjac, n, r, lr, qtf,
        wa1, wa2, wa3, wa4
    )

    msg = _hybrd_messages[max(-2, return_code)]
    if return_code < -1
        msg = msg * string(return_code)
    end

    coverged = return_code == 1
    trace.tot_time = time() - trace.start_time

    SolverResults("Modified Powell (Expert)", x0, x, fvec, return_code, coverged, msg, trace)
end

## Wrapping hybrj routine
const _hybrj_messages = Dict{Int, String}(
    0 => "improper input parameters",
    1 => string("algorithm estimates that the relative error between x and the ",
                "solution is at most tol"),
    2 => "maximum iterations has been exceeded",
    3 => "tol is too small, no further improvement in x is possible",
    4 => string("iteration is not making good progress, measured by ",
                "improvement from last 5 jacobian evaluations"),
    5 => string("iteration is not making good progress, measured by ",
                "improvement from last 10 iterations"),
    -1 => "exceeded user imposed number of iterations",
    -2 => "user terminated iterations with code "
)

function hybrj(f!::Function, g!::Function, x0::Vector{Float64}, xtol::Float64,
               show_trace::Bool, tracing::Bool, maxit::Int;
               _n::Int=length(x0), ml::Int=_n-1, mu::Int=_n-1,
               maxfev::Int=Int(typemax(Cint)),
               epsfcn::Float64=0.0, diag::Vector{Float64}=fill(1.0, _n),
               mode::Int=2, factor::Float64=100.0, nprint::Int=0,
               lr::Cint=ceil(Cint, _n*(_n+1)/2))
    n = length(x0)
    ldfjac = n
    x = copy(x0)
    fvec = similar(x)
    fjac = Array{Float64}(undef, n, n)
    lwa = ceil(Int, (n*(3*n+13))/2)
    wa = ones(lwa)
    trace = AlgoTrace(x0, show_trace, tracing, maxit)

    r = Array{Float64}(undef, lr)
    qtf = Array{Float64}(undef, n)
    wa1 = Array{Float64}(undef, n)
    wa2 = Array{Float64}(undef, n)
    wa3 = Array{Float64}(undef, n)
    wa4 = Array{Float64}(undef, n)

    if show_trace
        show(trace)
    end

    function _hybrj_func_wrapper(_p::Ptr{Cvoid}, n::Cint, _x::Ptr{Cdouble},
                                 _fvec::Ptr{Cdouble}, _fjac::Ptr{Cdouble},
                                 ldfjac::Cint, iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, n)
        fjac_flat = unsafe_wrap(Array, _fjac, ldfjac*n)
        fjac = reshape(fjac_flat, Int(ldfjac), Int(n))
        local x = unsafe_wrap(Array, _x, n)
        if iflag < 0
            print(fvec)
            return Cint(0)
        end

        if iflag == 1
            f!(fvec, x)
        elseif iflag == 2
            g!(fjac, x)
        end

        local trace = unsafe_pointer_to_objref(_p)::AlgoTrace
        push!(trace, x, fvec, iflag)

        trace.f_calls > trace.maxit ? Cint(-1) : Cint(0)
    end
    _hybrj_cfunc = @cfunction(
        $_hybrj_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint)
    )

    return_code = ccall(
        (:hybrj, cminpack),
        Cint,
        (
            Ptr{Cvoid}, Ptr{Cvoid}, Cint,
            Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Cint,
            Ptr{Cdouble}, Cint, Cdouble, Cint, Ptr{Cint}, Ptr{Cint},
            Ptr{Cdouble}, Cint,
            Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}
        ),
        _hybrj_cfunc, pointer_from_objref(trace), n,
        x, fvec, fjac, ldfjac, xtol, Cint(maxfev),
        diag, mode, factor, nprint, [0], [0],
        r, lr,
        qtf, wa1, wa2, wa3, wa4
    )

    msg = _hybrj_messages[max(-2, return_code)]
    if return_code < -1
        msg = msg * string(return_code)
    end

    coverged = return_code == 1 || norm(fvec, Inf) <= xtol
    trace.tot_time = time() - trace.start_time

    SolverResults("Modified Powell (User Jac, Expert)", x0, x, fvec, return_code, coverged, msg, trace)
end

## Wrapping hybrd1 routine
const _hybr_messages = Dict{Int, String}(
    0 => "improper input parameters",
    1 => string("algorithm estimates that the relative error between x and the ",
                "solution is at most tol"),
    2 => "maximum iterations has been exceeded",
    3 => "tol is too small, no further improvement in x is possible",
    4 => "iteration is not making good progress",
    -1 => "exceeded user imposed number of iterations",
    -2 => "user terminated iterations with code "
)

function hybrd1(f!::Function, x0::Vector{Float64}, tol::Float64,
                show_trace::Bool, tracing::Bool, maxit::Int)
    x = copy(x0)
    fvec = similar(x)
    n = length(x)
    lwa = ceil(Int, (n*(3*n+13))/2)
    wa = ones(lwa)
    trace = AlgoTrace(x0, show_trace, tracing, maxit)

    function _hybrd1_func_wrapper(_p::Ptr{Cvoid}, n::Cint, _x::Ptr{Cdouble},
                                 _fvec::Ptr{Cdouble}, iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, n)
        local x = unsafe_wrap(Array, _x, n)
        if iflag < 0
            print(fvec)
            return Cint(0)
        end
        f!(fvec, x)

        local trace = unsafe_pointer_to_objref(_p)::AlgoTrace
        push!(trace, x, fvec, iflag)

        trace.f_calls > trace.maxit ? Cint(-1) : Cint(0)
    end

    _hybrd1_cfunc = @cfunction(
        $_hybrd1_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint)
    )

    if show_trace
        show(trace)
    end
    return_code = ccall(
        (:hybrd1, cminpack),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Cint),
        _hybrd1_cfunc, pointer_from_objref(trace), n, x, fvec, tol, wa, lwa
    )

    msg = _hybr_messages[max(-2, return_code)]
    if return_code < -1
        msg = msg * string(return_code)
    end

    coverged = return_code == 1 || norm(fvec, Inf) <= tol
    trace.tot_time = time() - trace.start_time

    SolverResults("Modified Powell", x0, x, fvec, return_code, coverged, msg, trace)
end

## Wrapping lmdif1 routine
const _lmdif1_messages = Dict{Int, String}(
    0 => "improper input parameters",
    1 => string("algorithm estimates that the relative error between x and the ",
                "solution is at most tol"),
    2 => string("algorithm estimates that the relative error between x and the ",
                "solution is at most tol"),
    3 => string("algorithm estimates that the relative error in the sum of ",
                "squares and the relative error between x and the solution is ",
                "at most tol"),
    4 =>  "fvec is orthogonal to the columns of the jacobian to machine precision.",
    5 => "maximum iterations has been exceeded",
    6 => "tol is too small, no further reduction of sum of squares is possible",
    7 => "tol is too small, no further improvement in x is possible",
    -1 => "exceeded user imposed number of iterations",
    -2 => "user terminated iterations with code "
)

# NOTE: default doesn't always hold
function lmdif1(f!::Function, x0::Vector{Float64}, m::Int, tol::Float64,
                show_trace::Bool, tracing::Bool, maxit::Int)
    x = copy(x0)
    n = length(x)
    if n > m
        msg = "Must have at least as many variables as equations"
        throw(ArgumentError(msg))
    end

    fvec = Array{Float64}(undef, m)
    lwa = m*n+5*n+m
    iwa = Array{Int}(undef, n)
    wa = Array{Float64}(undef, lwa)
    trace = AlgoTrace(x0, show_trace, tracing, maxit)

    if show_trace
        show(trace)
    end

    function _lmdif1_func_wrapper(_p::Ptr{Cvoid}, m::Cint, n::Cint, _x::Ptr{Cdouble},
                                 _fvec::Ptr{Cdouble}, iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, m)
        local x = unsafe_wrap(Array, _x, n)
        if iflag < 0
            print(fvec)
            return Cint(0)
        end
        f!(fvec, x)

        local trace = unsafe_pointer_to_objref(_p)::AlgoTrace
        push!(trace, x, fvec, iflag)
        trace.f_calls > trace.maxit ? Cint(-1) : Cint(0)
    end
    _lmdif1_cfunc = @cfunction(
        $_lmdif1_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint)
    )

    return_code = ccall(
        (:lmdif1, cminpack),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Ptr{Cint}, Ptr{Cdouble}, Cint),
        _lmdif1_cfunc, pointer_from_objref(trace), m, n, x, fvec, tol, iwa, wa, lwa
    )

    msg = _lmdif1_messages[max(-2, return_code)]
    if return_code < -1
        msg = msg * string(return_code)
    end
    converged = return_code in [1, 2, 3] || norm(fvec, Inf) <= tol
    trace.tot_time = time() - trace.start_time

    SolverResults("Levenberg-Marquardt", x0, x, fvec, return_code, converged, msg, trace)
end

## Wrapping lmdif routine
const _lmdif_messages = Dict{Int, String}(
    0 => "improper input parameters",
    1 => string("Both actual and predicted relative errors in the sum of ",
                "squares are at most ftol"),
    2 => "relative error between two consecutive iterates is at most xtol",
    3 => string("Both actual and predicted relative errors in the sum of ",
                "squares are at most ftol",
                "\nAND relative error between two consecutive iterates is ",
                "at most xtol"),
    4 =>  string("the cosine of the angle between fvec and any column of the ",
                 "jacobian is at most gtol in absolute value"),
    5 => "number of calls to fcn has reached or exceeded maxfev",
    6 => "ftol is too small, no further reduction of sum of squares is possible",
    7 => "xtol is too small, no further improvement in x is possible",
    8 => string("gtol is too small. fvec is orthogonal to the columns of the ",
                "jacobian to machine precision"),
    -1 => "exceeded user imposed number of iterations",
    -2 => "user terminated iterations with code "
)

function lmdif(f!::Function, x0::Vector{Float64}, m::Int, tol::Float64,
                show_trace::Bool, tracing::Bool, maxit::Int;
                _n::Int = length(x0),
                gtol::Float64=0.0, ftol::Float64=tol, xtol::Float64=tol,
                epsfcn::Float64=0.0, mode::Int=1, nprint::Int=0,
                maxfev::Int=(_n+1)*200, factor::Float64=100.0,
                ldfjac::Int=m)
    x = copy(x0)
    n = length(x)
    if n > m
        msg = "Must have at least as many variables as equations"
        throw(ArgumentError(msg))
    end

    fvec = Array{Float64}(undef, m)
    lwa = m*n+5*n+m
    iwa = Array{Int}(undef, n)
    wa = Array{Float64}(undef, lwa)
    trace = AlgoTrace(x0, show_trace, tracing, maxit)

    diag = Array{Float64}(undef, n)
    nfev = [0]
    fjac = Array{Float64}(undef, m, n)
    ipvt = Array{Cint}(undef, n)
    qtf = Array{Float64}(undef, n)
    wa1 = Array{Float64}(undef, n)
    wa2 = Array{Float64}(undef, n)
    wa3 = Array{Float64}(undef, n)
    wa4 = Array{Float64}(undef, m)

    if show_trace
        show(trace)
    end

    function _lmdif_func_wrapper(_p::Ptr{Cvoid}, m::Cint, n::Cint, _x::Ptr{Cdouble},
                                 _fvec::Ptr{Cdouble}, iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, m)
        local x = unsafe_wrap(Array, _x, n)
        if iflag < 0
            print(fvec)
            return Cint(0)
        end
        f!(fvec, x)

        local trace = unsafe_pointer_to_objref(_p)::AlgoTrace
        push!(trace, x, fvec, iflag)
        trace.f_calls > trace.maxit ? Cint(-1) : Cint(0)
    end
    _lmdif_cfunc = @cfunction(
        $_lmdif_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint)
    )

    return_code = ccall(
        (:lmdif, cminpack),
        Cint,
        # func         p        m     n        x             fvec       ftol
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble,
        # xtol     gtol   maxfev  epsfcn    diag        mode   factor
        Cdouble, Cdouble, Cint, Cdouble, Ptr{Cdouble}, Cint, Cdouble,
        # nprint nfev       fjac       ldfjac   ipvt     qtf
        Cint, Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cdouble},
        # wa1           wa2           wa3            wa4
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
        _lmdif_cfunc, pointer_from_objref(trace), m, n, x, fvec, ftol, xtol,
        gtol, maxfev, epsfcn, diag, mode, factor, nprint, nfev, fjac, ldfjac,
        ipvt, qtf, wa1, wa2, wa3, wa4
    )

    msg = _lmdif_messages[max(-2, return_code)]
    if return_code < -1
        msg = msg * string(return_code)
    end
    converged = return_code in [1, 2, 3, 4] || norm(fvec, Inf) <= ftol
    trace.tot_time = time() - trace.start_time

    SolverResults("Levenberg-Marquardt (expert)", x0, x, fvec, return_code,
        converged, msg, trace)
end

## Wrapping lmder routine
const _lmder_messages = Dict{Int, String}(
    0 => "improper input parameters",
    1 => string("Both actual and predicted relative errors in the sum of ",
                "squares are at most ftol"),
    2 => "relative error between two consecutive iterates is at most xtol",
    3 => string("Both actual and predicted relative errors in the sum of ",
                "squares are at most ftol",
                "\nAND relative error between two consecutive iterates is ",
                "at most xtol"),
    4 =>  string("the cosine of the angle between fvec and any column of the ",
                 "jacobian is at most gtol in absolute value"),
    5 => "number of calls to fcn has reached or exceeded maxfev",
    6 => "ftol is too small, no further reduction of sum of squares is possible",
    7 => "xtol is too small, no further improvement in x is possible",
    8 => string("gtol is too small. fvec is orthogonal to the columns of the ",
                "jacobian to machine precision"),
    -1 => "exceeded user imposed number of iterations",
    -2 => "user terminated iterations with code "
)

function lmder(f!::Function, g!::Function, x0::Vector{Float64}, m::Int,
               tol::Float64, show_trace::Bool, tracing::Bool, maxit::Int;
               _n::Int=length(x0),
               gtol::Float64=0.0, ftol::Float64=tol, xtol::Float64=tol,
               epsfcn::Float64=0.0, mode::Int=1, nprint::Int=0,
               maxfev::Int=(_n+1)*200, factor::Float64=100.0,
               ldfjac::Int=m)
    x = copy(x0)
    n = length(x)
    if n > m
        msg = "Must have at least as many variables as equations"
        throw(ArgumentError(msg))
    end

    fvec = Array{Float64}(undef, m)
    lwa = m*n+5*n+m
    iwa = Array{Int}(undef, n)
    wa = Array{Float64}(undef, lwa)
    trace = AlgoTrace(x0, show_trace, tracing, maxit)

    diag = Array{Float64}(undef, n)
    nfev = [0]
    njev = [0]
    fjac = Array{Float64}(undef, m, n)
    ipvt = Array{Cint}(undef, n)
    qtf = Array{Float64}(undef, n)
    wa1 = Array{Float64}(undef, n)
    wa2 = Array{Float64}(undef, n)
    wa3 = Array{Float64}(undef, n)
    wa4 = Array{Float64}(undef, m)

    if show_trace
        show(trace)
    end

    function _lmder_func_wrapper(_p::Ptr{Cvoid}, m::Cint, n::Cint,
                                 _x::Ptr{Cdouble}, _fvec::Ptr{Cdouble},
                                 _fjac::Ptr{Cdouble}, ldfjac::Cint,
                                 iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, m)
        fjac_flat = unsafe_wrap(Array, _fjac, ldfjac*n)
        fjac = reshape(fjac_flat, Int(ldfjac), Int(n))
        local x = unsafe_wrap(Array, _x, n)
        if iflag < 0
            print(fvec)
            return Cint(0)
        end

        if iflag == 1
            f!(fvec, x)
        elseif iflag == 2
            g!(fjac, x)
        end

        local trace = unsafe_pointer_to_objref(_p)::AlgoTrace
        push!(trace, x, fvec, iflag)

        trace.f_calls > trace.maxit ? Cint(-1) : Cint(0)
    end

    _lmder_cfunc = @cfunction(
        $_lmder_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint)
    )

    return_code = ccall(
        (:lmder, cminpack),
        Cint,
        # func         p        m     n        x             fvec       fjac
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        # ldfjac  ftol  xtol     gtol   maxfev   diag        mode   factor
        Cint, Cdouble, Cdouble, Cdouble, Cint, Ptr{Cdouble}, Cint, Cdouble,
        # nprint  nfev    njev        ipvt        qtf         wa1
        Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
        # wa2            wa3          wa4
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),

        _lmder_cfunc, pointer_from_objref(trace), m, n, x, fvec, fjac,
        ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, nfev,
        njev, ipvt, qtf, wa1, wa2, wa3, wa4
    )

    msg = _lmder_messages[max(-2, return_code)]
    if return_code < -1
        msg = msg * string(return_code)
    end
    converged = return_code in [1, 2, 3, 4, 8] || norm(fvec, Inf) <= ftol
    trace.tot_time = time() - trace.start_time

    SolverResults("Levenberg-Marquardt (expert)", x0, x, fvec, return_code,
        converged, msg, trace)
end

## Wrapping fdjac1 routine
function fdjac1(f!::Function, x0::Vector{Float64};
               _n::Int=length(x0), ml::Int=_n-1, mu::Int=_n-1,
               epsfcn::Float64=0.0)
    n = length(x0)
    x = copy(x0)
    fvec = similar(x)
    fjac = Array{Float64}(undef, n, n)
    ldjfac = n

    f!(x, fvec)

    wa1 = Array{Float64}(undef, n)
    wa2 = Array{Float64}(undef, n)

    function _fdjac1_func_wrapper(_p::Ptr{Cvoid}, n::Cint, _x::Ptr{Cdouble},
                                  _fvec::Ptr{Cdouble}, iflag::Cint)
        local fvec = unsafe_wrap(Array, _fvec, n)
        local x = unsafe_wrap(Array, _x, n)
        f!(fvec, x)

        Cint(0)
    end

    _fdjac1_cfunc = @cfunction(
        $_fdjac1_func_wrapper,
        Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint)
    )

    return_code = ccall(
        (:fdjac1, cminpack),
        Cint,
        (
            Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
            Cint, Cint, Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}

        ),
        _fdjac1_cfunc, pointer_from_objref(nothing), n, x, fvec, fjac,
        ldjfac, ml, mu, epsfcn, wa1, wa2
    )

    return fjac
end


text:
#!/usr/bin/env julia

using JLD2, Distributions, Serialization
import GEAS:sim_pt_QTL, breeding_program, top_QTL

"""
    function run_2021_06_01(jld, nqtl, rt, nsib, nrpt; rst = "")
---
## Description
- 20 repeats
- 3 comparisons
  - Normal genomic selection, no known QTL
  - Known 5 QTL as fixed effects
  - Edit the biggest of known 5 QTL in generation 2-6

## parameters
1. jld: the base
2. nqtl: e.g., 500
3. rt: correlation of QTL effects between 2 traits
4. nsib: number of sibs challenged and in production. total 2sib
5. nrpt: number of repeats
6. rst: result path
"""
function run_2021_06_01(jld, nqtl, rt, nsib, nrpt; rst = "")
    @load jld base
    if length(rst) > 1
        isdir(rst) || mkpath(rst)
    end

    w = ndigits(nrpt)
    for r in 1:nrpt
        d = MvNormal(zeros(2), [1 rt; rt 1])
        qtl = sim_pt_QTL(base, nqtl, d)

        suffix = lpad(string(r), w, '0')
        serialize(joinpath(rst, "qtl-$suffix.ser"), qtl)
        par = begin
            Parameters = Dict(
                :nSire=> 100,
                :nDam => 200,       # male:female = 1:2
                :nSib => 2nsib,
                :nC7e => nsib,      # number of sib to be challenged
                :nG8n => 6,         # number of generations
                :nQTL => [nqtl, nqtl],
                :h²   => [.5, .5],
                :p8e  => .5,        # death rate after challenge
                :e19e => 1.,        # edit_successsful_rate
                :w4t  => 1,         # weight on the binary trait EBV
                :nk3n => 5          # number of known QTL on binary trait
            )
            (; Parameters...)       # named tuple.  contents as above
        end

        # genomic selection with no known QTL
        prd, snp = breeding_program(base, par, qtl)
        serialize(joinpath(rst, "prd-a-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-a-$suffix.ser"), snp)

        # genomic selection with known QTL as fixed effects
        prd, snp = breeding_program(base, par, qtl, fixed=true)
        serialize(joinpath(rst, "prd-b-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-b-$suffix.ser"), snp)

        # genomic selection, edit biggest known QTL each generation
        prd, snp = breeding_program(base, par, qtl, edit=true)
        serialize(joinpath(rst, "prd-c-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-c-$suffix.ser"), snp)
    end
end

# run_2021_06_01("dat/run/base.jld2", 500, -0.2, 40, 20, rst="dat")

#=
Note:
- this is a task discription to be submitted on slurm system
- below is slurm task bash script.
- copy below to `example.sh`
- uncomment the function call
- `sbatch example.sh`

#!/usr/bin/env bash
############################################################
# Environment setup
############################################################
#SBATCH --job-name=GEAS
#SBATCH --account=nn9891k
#SBATCH --job-name=GEAS
#SBATCH --time=0-10:00:00       # d-hh:mm:ss
#SBATCH --mem-per-cpu=2G        # 10x2 => 20G memory
#SBATCH --cpus-per-task=10      # 80 threads available
##sbatch --ntasks=6             # 190/20 => 9 independent tasks
##sbatch --ntasks-per-node=2    # not sure for the moment

# other job-controls may be needed:
#  - nodes=;
#  - ntasks-per-node=;
#  - cpus-per-task=;

log=results.geas

############################################################
# copy to workdirectory on the computation node
# i always write the working scripts in Julia
############################################################
cp $SUBMITDIR/geas.{sh,jl} $SCRATCH/
cp $SUBMITDIR/../dat/base.jld2 $SCRATCH

############################################################
# run the task
############################################################
cd $SCRATCH
echo "Time start: $(date)" > $log
./geas.jl
echo "Time stop: $(date)" >> $log

############################################################
# collect results
############################################################
cp *ser $log $SUBMITDIR/
=#

ext: jl


text:
######### Stan program example  ###########

using StanSample, MCMCChains

  binorm_model = "
  transformed data {
      matrix[2,2] Sigma;
      vector[2] mu;

      mu[1] = 0.0;
      mu[2] = 0.0;
      Sigma[1,1] = 1.0;
      Sigma[2,2] = 1.0;
      Sigma[1,2] = 0.10;
      Sigma[2,1] = 0.10;
  }
  parameters {
      vector[2] y;
  }
  model {
        y ~ multi_normal(mu, Sigma);
  }
  "

sm = SampleModel("binormal", binorm_model);

rc = stan_sample(sm)

if success(rc)
  chn = read_samples(sm, :mcmcchains)
  chn |> display
end


text:
using Deconvolution
#using PyPlot
#####################################
print("Program begins...... ")
tic()
println("Time in...")
########parameters setting###########
λ = 632.8E-9
k = 2*pi/λ
f = 200E-3
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Rav1 = 2100E-6
dR1 = 50E-6
R1 = Rav1-dR1/2
R2 = Rav1+dR1/2
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Rav2 = 4300E-6
dR2 = 200E-6
R3 = Rav2-dR2/2
R4 = Rav2+dR2/2
println("General parameters have been set...")
########range setting & ring mask###########
XYrange = linspace(-1E-0,1E-0,1000)         # 5000 shows the resolution of the mask
n = length(XYrange)
x = reshape(XYrange,1,n)
y = reshape(XYrange,n,1)
(X,Y) = (broadcast(+,x,zeros(n,1)),broadcast(+,y,zeros(1,n)))
ρ,ta,res2,res4 = zeros(n,n),zeros(n,n),zeros(n,n),zeros(n,n)
tb = ones(n,n)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
print("Arranging double-ring mask... ")
for ii = 1:1:n
  for jj = 1:1:n
    ρ[jj,ii] = sqrt(x[ii]^2 + y[jj]^2)
    if (ρ[jj,ii]>R1 && ρ[jj,ii]<R2)||(ρ[jj,ii]>R3 && ρ[jj,ii]<R4)
      ta[jj,ii]=1
    elseif ρ[jj,ii]==R1 || ρ[jj,ii]==R2 || ρ[jj,ii]==R3 || ρ[jj,ii]==R4
      ta[jj,ii]=0.5
    end
  end
end
println("Completed!")
########random points setting###########
pn = 2
polypx,polypy = [rand(pn) rand(pn)],[rand(pn) rand(pn)]
println("Creating random points...NumOfPoints = ",pn," ...")
print("Build the central Gaussian and Bessel... ")
W0 = 3E-3;
gaussian0 = exp(-(ρ.^2)/(W0^2))
res = (fftshift(abs(fft(gaussian0)))).^2
res = res./maximum(res)
res3 = (fftshift(abs(fft(ta.*gaussian0)))).^2
res3 = res3./maximum(res3)
println("Completed!")
########Fourier Transformation###########
println("Loop starts... Total turns = ",pn," .")
for cont = 1:pn
  gaussian0d = exp(-(ρ.^2)/(W0^2)).*exp(-1*im*k*Y*polypx[cont]-1*im*k*X*polypy[cont])
  res2 = res2+(fftshift(abs(fft(gaussian0d)))).^2
  res4 = res4+(fftshift(abs(fft(ta.*gaussian0d)))).^2
  println("The ",cont,"th turn has completed...")
end
print("Loop finishes... Organizing data... ")
res2 = res2./maximum(res2)
restot = res + res2
res4 = res4./maximum(res4)
restot2 = res3 + res4
println("Completed!")
print("Deconvolution step starts... Setting PSF... ")
PSF = res3
print("Completed!")
imobj = wiener(restot2,PSF,0)

toc()

max_stars_repo_name: xwghua/DRiP-public

max_issues_repo_name: xwghua/DRiP-public

max_forks_repo_name: xwghua/DRiP-public


text:
# Example which demonstrates removing an edge

using LightGraphs
using NetworkViz
using ThreeJS
using Escher

g = CompleteGraph(10)

@show typeof(g)

main(window) = begin
  push!(window.assets, "widgets")
  push!(window.assets,("ThreeJS","threejs"))

  inp = Signal(Dict{Any, Any}(:node1=>"1",:node2=>"3"))

  s = Escher.sampler()
  sampler = Escher.sampler()
  form = vbox(
      h1("Remove Edge"),
      watch!(s, :node1, textinput("", label="Source Vertex")),
      watch!(s, :node2, textinput("", label="Destination Vertex")),
      trigger!(s, :submit, button("Remove Edge"))
  ) |> maxwidth(400px)

  map(inp) do dict
      node1 = parse(Int,dict[:node1])
      node2 = parse(Int,dict[:node2])
      vbox(
          intent(s, form) >>> inp,
          vskip(2em),
          removeEdge(g, node1, node2, 0)
      ) |> Escher.pad(2em)
  end
end


text:
__precompile__()
module ReSHOP

import MathOptInterface
const MOI = MathOptInterface

import MathProgBase
const MPB = MathProgBase

using Compat
using LinearAlgebra, SparseArrays

export ReSHOPSolver, getsolvername, getsolveresult, getsolveresultnum, getsolvemessage, getsolveexitcode, LinearQuadraticModel

if VERSION <= v"1.3.0"
   # Load in `deps.jl`, complaining if it does not exist
   const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
   if !isfile(depsjl_path)
      error("ReSHOP not installed properly, run Pkg.build(\"ReSHOP\"), restart Julia and try again")
   end
   include(depsjl_path)
else
   using ReSHOP_jll
end

reshop_valid_index_max = Csize_t(0)

# This is will be called immediately after the module is loaded
function __init__()
   global reshop_valid_index_max = ccall((:reshop_get_index_max, libreshop), Csize_t, (),)
end

const CONFIG = Dict(
:debug => false,
:export_gms => false,
:solver_log => false
)

const solverdata_dir = joinpath(dirname(@__DIR__), ".solverdata")

include("reshop_data.jl")
include("reshop_linearity.jl")
include("reshop_expr.jl")
include("reshop_params.jl")
include("reshop_convert.jl")
include("reshop_fun.jl")
include("reshop_nlp.jl")
include("reshop_gams.jl")

mutable struct ReSHOPSolver <: MPB.AbstractMathProgSolver
    solver_name::String
    options::Dict{String,Any}
    emp::Union{Nothing,Function}
end

# TODO(xhub) write a better struct/enum here
@enum MODEL_TYPE qcp=7 nlp=2 miqcp=6 minlp=5 emp=10

"change the debug state"
setdebug(b::Bool) = global CONFIG[:debug] = b
"export the problem to a GAMS Model file (.gms)"
setexport(b::Bool) = global CONFIG[:export_gms] = b
"printout the log from the solver"
setsolverlog(b::Bool) = global CONFIG[:solver_log] = b

"""
Create a ReSHOPSolver Solver for MPB. The optional arguments are:

# Optional Arguments
- `solver_name::String=""`: solver used for this problem
- `options::Dict{String,Any}=Dict{String,Any}()`: the ReSHOP options

"""
function ReSHOPSolver(solver_name::String="",
                     options::Dict{String,Any}=Dict{String,Any}())
    ReSHOPSolver(solver_name, options, nothing)
end

include("reshop_utils.jl")

include("MBP_wrapper.jl")
include("MOI_wrapper.jl")

include("reshop_mathprgm.jl")
include("reshop_ovf.jl")
include("reshop_solve.jl")

end


text:
# Locally Linear Embedding (LLE)
# ------------------------
# Nonlinear dimensionality reduction by locally linear embedding,
# Roweis, S. & Saul, L., Science 290:2323 (2000)

#### LLE type
immutable LLE{T <: Real} <: SpectralResult
    k::Int
    λ::AbstractVector{T}
    proj::Projection{T}

    LLE{T}(k::Int, λ::AbstractVector{T}, proj::Projection{T}) = new(k, λ, proj)
end

## properties
outdim(M::LLE) = size(M.proj, 1)
projection(M::LLE) = M.proj

eigvals(M::LLE) = M.λ
neighbors(M::LLE) = M.k

## show & dump
function show(io::IO, M::LLE)
    print(io, "LLE(indim = $(indim(M)), outdim = $(outdim(M)), neighbors = $(neighbors(M)))")
end

function dump(io::IO, M::LLE)
    show(io, M)
    println(io, "eigenvalues: ")
    Base.showarray(io, M.λ', header=false, repr=false)
    println(io)
    println(io, "projection:")
    Base.showarray(io, M.proj, header=false, repr=false)
end

## interface functions
function transform{T<:Real}(::Type{LLE}, X::DenseMatrix{T}; d::Int=2, k::Int=12)
    n = size(X, 2)

    # Construct NN graph
    D, E = find_nn(X, k)

    # Select largest connected component
    cc = components(E)
    c = cc[indmax(map(size, cc))]
    if length(cc) == 1
        c = cc[1]
    else
        c = cc[indmax(map(size, cc))]
        # renumber edges
        R = Dict(c, 1:length(c))
        Ec = zeros(Int,k,length(c))
        for i = 1 : length(c)
            Ec[:,i] = map(i->get(R,i,0), E[:,c[i]])
        end
        E = Ec
        X = X[:,c]
    end

    if k > d
        warn("K>D: regularization will be used")
        tol = 1e-5
    else
        tol = 0
    end

    # Reconstruct weights
    W = zeros(k, n)
    for i = 1 : n
        Z = X[:, E[:,i]] .- X[:,i]
        C = Z'*Z
        C = C + eye(k, k) * tol * trace(C)
        wi = vec(C\ones(k, 1))
        W[:, i] = wi./sum(wi)
    end

    # Compute embedding
    M = speye(n,n)
    for i = 1 : n
        w = W[:, i]
        jj = E[:, i]
        M[i,jj] = M[i,jj] - w'
        M[jj,i] = M[jj,i] - w
        M[jj,jj] = M[jj,jj] + w*w'
    end

    λ, V = decompose(M, d)
    return LLE{T}(k, λ, scale!(V', sqrt(n)))
end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ c4381b3c-6b8f-11eb-34c1-9760c2fc6d87
begin
	using Test
end

# ╔═╡ 968a3130-6b8e-11eb-1d6c-1317dbf89bc8
md"""
## \#7
"""

# ╔═╡ b4d8b1a2-6b8e-11eb-3273-bd06e5e16a9b
begin
	b = [
		4
		1
		0
		1
	]

	x = [
		0
		1
		2
		3
	]

	A = [
		ones(4, 1) x
	]

	u = round.(Int, inv(A' * A) * A' * b)

	@test u == [
		3
		-1
	]
end

# ╔═╡ d2278240-6b90-11eb-0bc3-df340513edd4
md"""
## \#8
"""

# ╔═╡ e6c51852-6b90-11eb-1e64-1b29b65808fd
begin
	p = A * u

	e = b - p

	@test A' * e == [
		0
		0
	]
end

# ╔═╡ 00970928-6b93-11eb-1707-f157198becdd
md"""
## \#12
"""

# ╔═╡ 09012d82-6b93-11eb-00ab-41d5b7a8615e
begin
	A₂ = [
		A x.^2
	]

	u₂ = round.(Int, inv(A₂' * A₂) * A₂' * b)

	@test u₂ == [
		4
		-4
		1
	]

	p₂ = A₂ * u₂

	e₂ = b - p₂

	# Exact solution
	e₂ == [
		0.0
		0.0
		0.0
		0.0
	]

	A₃ = [
		A₂ x.^3
	]

	u₃ = round.(Int, inv(A₃' * A₃) * A₃' * b)

	p₃ = A₃ * u₃

	e₃ = b - p₃

	# Exact solution for cubic
	e₃ == [
		0.0
		0.0
		0.0
		0.0
	]
end

# ╔═╡ a50a2640-6ba0-11eb-1247-2388f1b38906
md"""
## \#22
"""

# ╔═╡ aa83f89c-6ba0-11eb-32ae-8ba9f30015eb
begin
	A₄ = [
		1 -1
		1 1
		1 2
	]

	b₄ = [
		5
		13
		17
	]

	u₄ = round.(Int, inv(A₄' * A₄) * A₄' * b₄)

	@test u₄ == [
		9
		4
	]

	e₄ = A₄ * u₄ - b₄

	@test e₄ == [
		0
		0
		0
	]
end

# ╔═╡ Cell order:
# ╠═c4381b3c-6b8f-11eb-34c1-9760c2fc6d87
# ╟─968a3130-6b8e-11eb-1d6c-1317dbf89bc8
# ╠═b4d8b1a2-6b8e-11eb-3273-bd06e5e16a9b
# ╟─d2278240-6b90-11eb-0bc3-df340513edd4
# ╠═e6c51852-6b90-11eb-1e64-1b29b65808fd
# ╟─00970928-6b93-11eb-1707-f157198becdd
# ╠═09012d82-6b93-11eb-00ab-41d5b7a8615e
# ╠═a50a2640-6ba0-11eb-1247-2388f1b38906
# ╠═aa83f89c-6ba0-11eb-32ae-8ba9f30015eb


text:
l_bottom_vertex = 2
l_top_vertex = 5
r_top_vertex = 7
r_bottom_vertex = 11

mf = TrapezoidalMF(l_bottom_vertex, l_top_vertex, r_top_vertex, r_bottom_vertex)

# Evaluation tests
@assert mf.eval(l_top_vertex) == mf.eval(r_top_vertex) == 1
@assert mf.eval(l_bottom_vertex) == mf.eval(r_bottom_vertex) == 0
@assert mf.eval((l_bottom_vertex + l_top_vertex) / 2) == mf.eval((r_top_vertex + r_bottom_vertex) / 2) == 0.5

# Mean finding tests
@assert mf.mean_at(1) == (l_top_vertex + r_top_vertex) / 2
@assert mf.mean_at(0) == (l_bottom_vertex + r_bottom_vertex) / 2


text:
using DSP, Plots

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

for LENGTH=3:5
    for i = 1:10
        x = rand((-1,0,1), LENGTH)
        h = rand((-1,0,1), LENGTH)
        y = conv(x, h)

        p1 = plot(1:LENGTH, x,
                    seriestype=:bar,
                    title="x[n]",
                    linealpha=0,
                    color=:blue,
                    alpha=0.7,
                    legend=:none,
                    showaxis=true,
                    ylims=(-1,1))

        p2 = plot(1:LENGTH, h,
                    seriestype=:bar,
                    title="h[n]",
                    linealpha=0,
                    color=:grey,
                    alpha=0.7,
                    legend=:none,
                    showaxis=true,
                    ylims=(-1,1))

        p3 = plot(1:LENGTH*2-1,
                    y,
                    seriestype=:bar,
                    linealpha=0,
                    alpha=0.7,
                    color=:red,
                    title="y[n] = x[n] * h[n] = h[n] * x[n]",
                    ylims=(-5,5),
                    showaxis=true,
                    minorgrid=true,
                    legend=:none)
        plot!(1:LENGTH*2-1,
                    [x;zeros(LENGTH-1)],
                    seriestype=:bar,
                    linealpha=0,
                    alpha=0.2,
                    color=:blue)
        plot!(1:LENGTH*2-1,
                    [h;zeros(LENGTH-1)],
                    seriestype=:bar,
                    linealpha=0,
                    alpha=0.1,
                    color=:grey)

        question_layout = @layout[a;b]

        if (!isdir("images"))
            mkdir("images")
        end

        plot(p1, p2, layout = question_layout, title="Convolve these two series.")
        png("images/$(LENGTH)_$(i)_side-a_question")

        plot(p3, title="Convolution in red! Did you get it?")
        png("images/$(LENGTH)_$(i)_side-b_answer")
    end
end


text:
module MatLang

export clcM
include("Language_Fundamentals/Entering_Commands.jl")

using LinearAlgebra

export zerosM, onesM, randM, trueM, falseM, eyeM, diagM, blkdiagM, catM, horzcatM, vertcatM, repelemM, repmatM, linspaceM, logspaceM, meshgridM, ndgridM, freqspaceM, lengthM, sizeM, ndimsM, numelM, isscalarM, isvectorM, ismatrixM, isrowM, iscolumnM, sortM, flipM, fliplrM, flipudM, flipdimM, rot90M, transposeM, squeezeM
include("Language_Fundamentals/Matrices_and_Arrays.jl")

include("Language_Fundamentals/Data_Types/Numeric_Types.jl")
################################################################
# Conditional REQUIRE - https://github.com/MikeInnes/Requires.jl
using Requires
function __init__()
    # clcM() in Juno
    @require Juno = "e5e0dc1b-0480-54bc-9374-aad01c23163d" begin
        function clcM()
            println("\33[2J")
            Juno.clearconsole()
        end
    end
end
################################################################

include("precompile_includer.jl")
end # module

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
__precompile__(true)

module SortMerge

using Printf, SparseArrays

import Base.show
import Base.sortperm
import Base.zip

import Base.iterate, Base.length, Base.size, Base.getindex,
       Base.firstindex, Base.lastindex, Base.IndexStyle

export sortmerge, sortperm, nmatch, countmatch, multimatch

# --------------------------------------------------------------------
# Source structure
#
struct Source
    size::Int
    sortperm::Vector{Int}
    cmatch::SparseVector{Int,Int}
end

# --------------------------------------------------------------------
# Matched structure
#
struct Matched <: AbstractVector{Vector{Int}}
    nsrc::Int
    nrow::Int
    sources::Vector{Source}
    matched::Matrix{Int}
end

# Methods
length(mm::Matched) = mm.nsrc
size(mm::Matched) = (mm.nsrc,)
getindex(mm::Matched, inds) = mm.matched[:, inds]
firstindex(mm::Matched) = 1
lastindex(mm::Matched) = mm.nsrc
IndexStyle(::Type{Matched}) = IndexLinear()

nmatch(mm::Matched) = mm.nrow
countmatch(mm::Matched, source::Int) = mm.sources[source].cmatch
sortperm(  mm::Matched, source::Int) = mm.sources[source].sortperm

zip(mm::Matched) = zip(map(i -> mm.matched[:,i], 1:mm.nsrc)...)
#unzip(mm::Matched)  = ntuple(i -> mm.matched[:,i],  length(mm.sources))

function multimatch(mm::Matched, source::Int, multi::Int; group=false)
    @assert multi >= 1
    index = findall(countmatch(mm, source) .== multi)
    index_out = sortmerge(index, mm.matched[:,source], quiet=true)[2]
    matrix = mm.matched[index_out,:]

    out = Vector{Matched}()
    ngroups = (group  ?  length(index)  :  1)
    for igroup in 1:ngroups
        jj = (group  ?
              findall(matrix[:,source] .== index[igroup])  :
              collect(1:length(index_out)))

        sources = Vector{Source}()
        for i in 1:mm.nsrc
            push!(sources, Source(mm.sources[i].size, Vector{Int}(),
                                  sparsevec(unique(mm.matched[index_out[jj], i]),
                                            multi, length(mm.sources[i].cmatch))))
        end
        push!(out, Matched(mm.nsrc, length(jj), sources, mm.matched[index_out[jj],:]))
    end

    (group)  &&  (return out)
    return out[1]
end

function default_lt(v, i, j)
    return v[i] < v[j]
end
function default_sd(A, B, i, j)
    return sign(A[i] - B[j])
end

# sortmerge(j::NTuple{2, Matched}, A, B, args...; sd=default_sd, quiet=false) =
#   sortmerge(A, B, args..., sort1=j[1].sortperm, sort2=j[2].sortperm, sd=sd, quiet=quiet)
function sortmerge(A, B, args...;
                   sd=default_sd,
                   quiet=false,
                   sort1=Vector{Int}(),
                   sort2=Vector{Int}(),
                   lt1=default_lt,
                   lt2=default_lt,
                   sorted=false)

    elapsed = (Base.time_ns)()
    size1 = size(A)[1]
    size2 = size(B)[1]

    elapsed_sorting = fill(0., 2)
    elapsed_sorting[1] = (Base.time_ns)()
    if length(sort1) == 0
        sort1 = collect(range(1, length=size1))
        (sorted)  ||  (sort1 = sortperm(sort1, lt=(i, j) -> (lt1(A, i, j))))
    end
    elapsed_sorting[1] = ((Base.time_ns)() - elapsed_sorting[1]) / 1.e9
    elapsed_sorting[2] = (Base.time_ns)()
    if length(sort2) == 0
        sort2 = collect(range(1, length=size2))
        (sorted)  ||  (sort2 = sortperm(sort2, lt=(i, j) -> (lt2(B, i, j))))
    end
    elapsed_sorting[2] = ((Base.time_ns)() - elapsed_sorting[2]) / 1.e9

    i2a = 1
    match1 = Array{Int}(undef, 0)
    match2 = Array{Int}(undef, 0)
    cm1 = fill(0, size1)
    cm2 = fill(0, size2)

    lastlog = -1.
    progress = false
    for i1 in 1:size1
        for i2 in i2a:size2
            if !quiet  &&  !progress  &&  (((Base.time_ns)() - elapsed) / 1.e9 > 1)
                progress = true
            end
            if !quiet  &&  progress
                completed = ceil(1000. * ((i1-1) * size2 + i2) / (size1 * size2))
                if completed > lastlog
                    @printf("Completed: %5.1f%%, matched: %d \r", completed/10., length(match1))
                    lastlog = completed
                end
            end

            j1 = sort1[i1]
            j2 = sort2[i2]
            if length(args) > 0  # This improves performances
                dd = Int(sd(A, B, j1, j2, args...))
            else
                dd = Int(sd(A, B, j1, j2))
            end

            if     dd == -1; break
            elseif dd ==  1; i2a += 1
            elseif dd ==  0
                push!(match1, j1)
                push!(match2, j2)
                cm1[j1] += 1
                cm2[j2] += 1
            end
        end
    end
    if !quiet  &&  progress
        @printf("Completed: %5.1f%%, matched: %d \n", 100., length(match1))
    end

    ii = findall(cm1 .> 0); side1 = Source(size1, sort1, sparsevec(ii, cm1[ii], length(cm1)))
    ii = findall(cm2 .> 0); side2 = Source(size2, sort2, sparsevec(ii, cm2[ii], length(cm2)))
    mm = Matched(2, length(match1), [side1, side2], [match1 match2])
    elapsed = ((Base.time_ns)() - elapsed) / 1.e9

    if !quiet
        for i in 1:length(mm.sources)
            ss = mm.sources[i]
            uu = length(unique(mm.matched[:,i]))
            @printf("Input %2d: %10d / %10d  (%6.2f%%)  -  max mult. %d | sort : %.3gs\n",
                    i, uu, ss.size, 100. * uu/float(ss.size), maximum(ss.cmatch), elapsed_sorting[i])
        end
        @printf("Output  : %10d                                         | total: %.3gs\n",
                size(mm.matched)[1], elapsed)
    end
    return mm
end

end # module


text:
#=
Functions used by the CGSolver for matrix free solutions.
=#

# Note: This uses the conjugate gradient iterative method,
# which assumes an SPD matrix.
function solve_matrix_free_sym(var, bilinear, linear, stepper=nothing; assemble_func=nothing)
    start_time = time_ns();
    tol = config.linalg_matfree_tol;
    maxiters = config.linalg_matfree_max;

    # If more than one variable
    if typeof(var) <: Array
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        dofs_per_loop = 0;
        for vi=1:length(var)
            dofs_per_loop += length(var[vi].symvar);
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_loop = length(var.symvar);
        dofs_per_node = var.total_components;
    end

    N1 = size(grid_data.allnodes,2);
    Nn = dofs_per_node * N1;
    Np = refel.Np;
    nel = size(grid_data.loc2glb,2);

    # For partitioned meshes, keep these numbers handy
    (b_order, b_sizes) = get_partitioned_ordering(N1, dofs_per_node);

    # Allocate arrays that will be used by assemble
    rhsvec = zeros(Nn);
    allocated_vecs = [rhsvec];

    if prob.time_dependent && !(stepper === nothing)
        #TODO
    else
        # Use regular rhs assembly
        b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, 0, 0; rhs_only = true, assemble_loops=assemble_func)
        b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes, rescatter_b=true);

        normb = norm(b, Inf);

        x = zeros(Nn);

        # Do initial matvec
        #Ax = elem_matvec(x, bilinear, dofs_per_node, var);
        # But this will just be zeros

        r0 = copy(b); # = b - Ax;
        r1 = copy(b);
        p = copy(r0);

        iter = 0;
        err = 1;
        while iter < maxiters && err > tol
            iter = iter+1;

            # Parts that require the global communication###########
            Ap = elem_matvec(p,bilinear, dofs_per_node, var, 0, 0);
            Ap = gather_system(nothing, Ap, N1, dofs_per_node, b_order, b_sizes, rescatter_b=true);

            r0sq = owned_dot(r0,r0,dofs_per_node);
            pAp = owned_dot(p,Ap,dofs_per_node);
            tmp = combine_values([r0sq, pAp]); # tmp[1] = dot(r0,r0), tmp[2] = dot(p,Ap)
            alpha = tmp[1] / tmp[2];

            for i=1:Nn
                r1[i] -= alpha * Ap[i]; # r1 = r0 .- alpha.*Ap;
            end
            beta = combine_values(owned_dot(r1,r1,dofs_per_node)) / tmp[1];
            ########################################################

            for i=1:Nn
                x[i] += alpha * p[i];   # x = x .+ alpha.*p;
                p[i] = r1[i] + beta * p[i]; # p = r1 .+ beta.*p;
                r0[i] = r1[i]; # r0 = copy(r1);
            end

            err = norm(r0, Inf)/normb;

            err = combine_values(err) / config.num_procs;

            if iter%50 == 0
                log_entry("iteration "*string(iter)*": res = "*string(err));
            end
        end

        #println("Converged to "*string(err)*" in "*string(iter)*" iterations");
        log_entry("Converged to "*string(err)*" in "*string(iter)*" iterations");

        total_time = time_ns() - start_time;
        log_entry("Matrix-free solve took "*string(total_time/1e9)*" seconds");

        return x;
    end
end

# Note: This uses the stabilized biconjugate gradient method which works for nonsymmetric matrices
function solve_matrix_free_asym(var, bilinear, linear, stepper=nothing, t=0, dt=0)
    start_time = time_ns();
    tol = config.linalg_matfree_tol;
    maxiters = config.linalg_matfree_max;
    Np = refel.Np;
    nel = mesh_data.nel;
    N1 = size(grid_data.allnodes,2);
    multivar = typeof(var) <: Array;
    if multivar  # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        for vi=1:length(var)
            tmp = dofs_per_node;
            dofs_per_node += length(var[vi].symvar);
        end
    else
        # one variable
        dofs_per_node = length(var.symvar);
    end
    Nn = dofs_per_node * N1;

    if prob.time_dependent && !(stepper === nothing)
        #TODO
    else
        # Use regular rhs assembly
        b = assemble_rhs_only(var, linear, t, dt);
        normb = norm(b, Inf);

        x = zeros(Nn);

        r0 = copy(b); # = b - Ax; but x is initially zeros
        rs = copy(b);
        rho0 = 1;
        alpha = 1;
        w0 = 1;
        p = zeros(Nn);
        Ap = zeros(Nn);
        As = zeros(Nn);

        iter = 0;
        err = 1;
        while iter < maxiters && err > tol
            iter = iter+1;

            rho1 = dot(rs,r0);
            beta = (rho1*alpha) / (rho0*w0);
            p = r0 + beta*(p - w0*Ap);

            Ap = elem_matvec(p,bilinear, dofs_per_node, var, t, dt);

            alpha = rho1/dot(rs,Ap);
            s = r0 - alpha*Ap;

            As = elem_matvec(s,bilinear, dofs_per_node, var, t, dt);

            w1 = dot(s,As) / dot(As,As);
            x = x + alpha*p + w1*s;
            r0 = s - w1*As;

            err = norm(r0, Inf)/normb;

            if iter%50 == 0
                log_entry("iteration "*string(iter)*": res = "*string(err));
            end
        end

        log_entry("Converged to "*string(err)*" in "*string(iter)*" iterations");

        total_time = time_ns() - start_time;
        log_entry("Matrix free solve took "*string(total_time/1e9)*" seconds");

        return x;
    end
end

# Does the elemental matvec Ax=b
# x is input, b is output, A is from bilinear
function elem_matvec(x, bilinear, dofs_per_node, var, t = 0.0, dt = 0.0)
    Np = refel.Np;
    nel = mesh_data.nel;
    multivar = typeof(var) <: Array;
    maxvarindex = 0;
    if multivar # multiple variables being solved for simultaneously
        for vi=1:length(var)
            maxvarindex = max(maxvarindex,var[vi].index);

            # check for neumann bcs
            for bi=1:length(prob.bc_type[vi,:])
                if prob.bc_type[vi,bi] == NEUMANN
                    printerr("Neumann boundary conditions are not yet supported for matrix free. Sorry.", fatal=true)
                end
            end
        end
    else
        # one variable
        maxvarindex = var.index;
    end

    Ax = zeros(size(x));

    # Stiffness and mass are precomputed for uniform grid meshes
    if config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
        glb = grid_data.loc2glb[:,1];
        xe = grid_data.allnodes[:,glb[:]];
        (detJ, J) = geometric_factors(refel, xe);
        wgdetj = refel.wg .* detJ;
        if config.dimension == 1
            (RQ1, RD1) = build_deriv_matrix(refel, J);
            TRQ1 = RQ1';
            stiffness = [(TRQ1 * diagm(wgdetj) * RQ1)];
        elseif config.dimension == 2
            (RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);
            (TRQ1, TRQ2) = (RQ1', RQ2');
            stiffness = [(TRQ1 * diagm(wgdetj) * RQ1) , (TRQ2 * diagm(wgdetj) * RQ2)];
        else
            (RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);
            (TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');
            stiffness = [(TRQ1 * diagm(wgdetj) * RQ1) , (TRQ2 * diagm(wgdetj) * RQ2) , (TRQ3 * diagm(wgdetj) * RQ3)];
        end
        mass = (refel.Q)' * diagm(wgdetj) * refel.Q;
    else
        stiffness = 0;
        mass = 0;
    end

    #Elemental loop follows elemental ordering
    for e=elemental_order;
        eid = elemental_order[e];
        loc2glb = grid_data.loc2glb[:,eid]; # global indices of this element's nodes

        subx = extract_linear(x, loc2glb, dofs_per_node);

        volargs = (var, eid, 0, grid_data, geo_factors, refel, t, dt, stiffness, mass);
        bilinchunk = bilinear.func(volargs); # the elemental bilinear part
        # Neumann bcs need to be applied to this
        # TODO

        if dofs_per_node == 1
            Ax[loc2glb] = Ax[loc2glb] + bilinchunk * subx;
        else
            insert_linear_matfree!(Ax, bilinchunk*subx, loc2glb, 1:dofs_per_node, dofs_per_node);
        end
    end

    # Apply Dirichlet boudary conditions
    bidcount = length(grid_data.bids); # the number of BIDs
    if multivar
        d = 0;
        for vi=1:length(var)
            for compo=1:length(var[vi].symvar)
                d = d + 1;
                for bid=1:bidcount
                    if prob.bc_type[var[vi].index, bid] == NO_BC
                        # do nothing
                    else
                        if config.num_partitions > 1
                            owned_bdry = similar(grid_data.bdry[bid]);
                            borrowed_bdry = similar(grid_data.bdry[bid]);
                            next_o_ind=1;
                            next_b_ind = 1;
                            for ni=1:length(owned_bdry)
                                if grid_data.node_owner[grid_data.bdry[bid][ni]] == config.partition_index
                                    owned_bdry[next_o_ind] = grid_data.bdry[bid][ni];
                                    next_o_ind += 1;
                                else
                                    borrowed_bdry[next_b_ind] = grid_data.bdry[bid][ni];
                                    next_b_ind += 1;
                                end
                            end
                            owned_bdry = owned_bdry[1:next_o_ind-1];
                            borrowed_bdry = borrowed_bdry[1:next_b_ind-1];
                            Ax = zero_bc_matfree(Ax, x, borrowed_bdry, d, dofs_per_node);
                        else
                            owned_bdry = grid_data.bdry[bid];
                        end
                        Ax = dirichlet_bc_matfree(Ax, x, owned_bdry, d, dofs_per_node);
                    end
                end
            end
        end
    else
        for bid=1:bidcount
            if prob.bc_type[var.index, bid] == NO_BC
                # do nothing
            else
                if config.num_partitions > 1
                    owned_bdry = similar(grid_data.bdry[bid]);
                    borrowed_bdry = similar(grid_data.bdry[bid]);
                    next_o_ind=1;
                    next_b_ind = 1;
                    for ni=1:length(owned_bdry)
                        if grid_data.node_owner[grid_data.bdry[bid][ni]] == config.partition_index
                            owned_bdry[next_o_ind] = grid_data.bdry[bid][ni];
                            next_o_ind += 1;
                        else
                            borrowed_bdry[next_b_ind] = grid_data.bdry[bid][ni];
                            next_b_ind += 1;
                        end
                    end
                    owned_bdry = owned_bdry[1:next_o_ind-1];
                    borrowed_bdry = borrowed_bdry[1:next_b_ind-1];
                    Ax = zero_bc_matfree(Ax, x, borrowed_bdry, 1:dofs_per_node, dofs_per_node);
                else
                    owned_bdry = grid_data.bdry[bid];
                end
                Ax = dirichlet_bc_matfree(Ax, x, owned_bdry, 1:dofs_per_node, dofs_per_node);
            end
        end
    end

    # Reference points
    if size(prob.ref_point,1) >= maxvarindex
        if multivar
            posind = zeros(Int,0);
            vals = zeros(0);
            dof_offset = 1;
            for vi=1:length(var)
                if prob.ref_point[var[vi].index,1]
                    eii = prob.ref_point[var[vi].index, 2];
                    tmp = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + dof_offset;
                    if length(prob.ref_point[var[vi].index, 3]) > 1
                        tmp = tmp:(tmp+length(prob.ref_point[var[vi].index, 3])-1);
                    end
                    posind = [posind; tmp];
                    vals = [vals; prob.ref_point[var[vi].index, 3]];
                end
                dof_offset += length(var[vi].symvar);
            end
            if length(vals) > 0
                Ax[posind] = vals;
            end

        else
            if prob.ref_point[var.index,1]
                eii = prob.ref_point[var.index, 2];
                posind = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + 1;
                if length(prob.ref_point[var.index, 3]) > 1
                    posind = posind:(posind+length(prob.ref_point[var[vi].index, 3])-1);
                end
                Ax[posind] = prob.ref_point[var.index, 3];
            end
        end
    end

    return Ax;
end

# does a dot product, but only adds components corresponding to owned nodes
function owned_dot(a, b, dofs_per_node)
    if config.num_partitions > 1
        val = 0.0;
        if dofs_per_node > 1
            for ni=1:size(grid_data.allnodes,2)
                if grid_data.node_owner[ni] == config.partition_index
                    offset = (ni-1)*dofs_per_node;
                    for di=1:dofs_per_node
                        val += a[offset+di]*b[offset+di];
                    end
                end
            end
        else
            for i=1:length(a)
                if grid_data.node_owner[i] == config.partition_index
                    val += a[i]*b[i];
                end
            end
        end

        return val;
    else
        return dot(a, b);
    end
end

# sets b[bdry] = x[bdry]
function dirichlet_bc_matfree(b, x, bdryind, dofind=1, totaldofs=1)
    if totaldofs > 1
        for d=dofind
            ind2 = totaldofs.*(bdryind .- 1) .+ d;
            b[ind2] = x[ind2];
        end
    else
        b[bdryind] = x[bdryind];
    end

    return b;
end

# sets b[bdry] = 0
function zero_bc_matfree(b, x, bdryind, dofind=1, totaldofs=1)
    if totaldofs > 1
        for d=dofind
            ind2 = totaldofs.*(bdryind .- 1) .+ d;
            b[ind2] .= 0;
        end
    else
        b[bdryind] .= 0;
    end

    return b;
end

function extract_linear(b, glb, dofs)
    if dofs == 1
        return b[glb];
    else
        np = length(glb);
        part = zeros(np*dofs);

        for d=1:dofs
            ind1 = ((d-1)*np+1):(d*np);
            ind2 = dofs.*(glb .- 1) .+ d;
            part[ind1] = b[ind2];
        end

        return part;
    end
end

# Inset the single dof into the greater construct
function insert_linear_matfree!(b, bel, glb, dof, Ndofs)
    # group nodal dofs
    np = length(glb);
    for d=1:length(dof)
        ind1 = ((d-1)*np+1):(d*np);
        ind2 = Ndofs.*(glb .- 1) .+ d;

        b[ind2] = b[ind2] + bel[ind1];
    end
end

function block_to_interlace(x, Np, dofs)
    newx = similar(x);
    for d=1:dofs
        newx[d:dofs:end] = x[((d-1)*Np+1):(d*Np)];
    end
    return newx;
end
function interlace_to_block(x, Np, dofs)
    newx = similar(x);
    for d=1:dofs
        newx[((d-1)*Np+1):(d*Np)] = x[d:dofs:end];
    end
    return newx;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using ProgressView
using Test

function inner()
    sleep(1.1)
    return 2
end

function fun()
    sleep(1.1)
    x = @monitor "inner" inner()
    @test x == 2
    sleep(1.1)
    return nothing
end

@testset "ProgressView" begin
    @monitor "testset" begin
        @monitor "fun" fun()
    end
    sleep(1.1)
    @monitor "fun" fun()
    monitor_flush()
end


text:
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(dirname(dirname(@__DIR__)))
using ArgParse

@everywhere begin
    import Dates
    import Random
    using ArgParse
    using CalibrateEDMF
    using CalibrateEDMF.HelperFuncs
    using CalibrateEDMF.DistributionUtils
    using CalibrateEDMF.TurbulenceConvectionUtils
    using Combinatorics
    using TurbulenceConvection
    tc = pkgdir(TurbulenceConvection)
    include(joinpath(tc, "driver", "generate_namelist.jl"))

    """
        get_parameter_pairs(namelist, p1, p2, v1, v2)

    Get dictionary of parameter name-value pairs given parameters `p1` and `p2` with values
    `v1` and `v2`, respectively. If either p1 or p2 are components of a vector parameter,
    the returned dictionary also contains name-value pairs for all components of that
    vector parameter. The other components are defined in the `namelist`.
    """
    function get_parameter_pairs(namelist::Dict, p1::S, p2::S, v1::Number, v2::Number) where {S <: AbstractString}
        turbconv_params = namelist["turbulence"]["EDMF_PrognosticTKE"]  # TODO: Generalize for vector parameters from other parts of the namelist

        # If any parameters are vectors, fetch relevant defaults from namelist
        vec_param_names = getfield.(filter(!isnothing, match.(r".*(?=_{\d+})", [p1, p2])), :match)
        namelist_vec_params = filter(p -> p.first ∈ vec_param_names, turbconv_params)
        # expand vector parameters to index form
        flatten_vector_parameter(p) = "$(p.first)_{" .* string.(1:length(p.second)) .* "}" .=> p.second
        params = Dict((map(flatten_vector_parameter, collect(namelist_vec_params))...)...)
        # update params with custom parameters
        params[p1] = v1
        params[p2] = v2
        return params
    end

    struct SimConfig
        "Name of first parameter"
        parameter1::String
        "Name of second parameter"
        parameter2::String
        "Value of first parameter"
        value1::Number
        "Value of second parameter"
        value2::Number
        "Ensemble index"
        ens_i::Integer
        "Case name"
        case_name::String
        "Case name index"
        case_id::Integer
        "Simulation output root directory"
        output_root::String
        "Additional namelist arguments"
        namelist_args::Union{AbstractVector, Nothing}
    end

    """
        run_sims(t) = run_sims(t...)
        run_sims(param_values::Tuple{Number, Number}, case::String, ens_ind::Integer, nt::NamedTuple)

    Run one forward simulation of the SCM given a case and pair of parameter values to be modified.

    Parameters:
    - param_values  :: Tuple of two values that are to be modified when running the model
    - case          :: Name of the simulation case to be run (e.g. Bomex)
    - ens_ind       :: Integer value for which ensemble member is to be run
    - nt            :: Named tuple that contains information that is constant across all simulations.
        Has the entires `output_dir`, `param_pair`, `namelist_args`, denoting the simulation output
        directory, names of parameters corresponding to `param_values`, and additional namelist
        arguments to be modified with respect to the default namelist, respectively.
    """
    run_sims(s::SimConfig) = run_sims(
        s.parameter1,
        s.parameter2,
        s.value1,
        s.value2,
        s.ens_i,
        s.case_name,
        s.case_id,
        s.output_root,
        s.namelist_args,
    )
    function run_sims(
        param1::S,
        param2::S,
        value1::Number,
        value2::Number,
        ens_i::Integer,
        case::S,
        case_id::Integer,
        output_root::S,
        namelist_args::Union{AbstractVector, Nothing},
    ) where {S <: AbstractString}
        # Create path to store forward model output
        case_dir = joinpath(output_root, "$param1.$param2/$(value1)_$(value2)/$case.$case_id")  # e.g. output/220101_abc/param1.param2/0.1_0.2/Bomex.1

        # If the output simulation directory already exists, don't run a simulation for this configuration.
        output_dir = joinpath(case_dir, "Output.$case.$ens_i")
        if isdir(output_dir)
            return nothing
        end

        # Otherwise, create the case directory and run a simulation
        mkpath(case_dir)

        # Get namelist for case
        namelist = NameList.default_namelist(case, write = false, set_seed = false)

        params = get_parameter_pairs(namelist, param1, param2, value1, value2)

        # Run forward model
        run_SCM_handler(
            case,
            case_dir;
            u = collect(Float64, values(params)),
            u_names = collect(String, keys(params)),
            namelist = namelist,
            namelist_args = namelist_args,
            uuid = "$ens_i",
            les = get(namelist["meta"], "lesfile", nothing),
        )
    end
end  # end @everywhere

function grid_search(config::Dict, config_path::String)
    root = get_entry(config["grid_search"], "output_root_dir", pwd())
    now = Dates.format(Dates.now(), "YYmmdd")
    suffix = Random.randstring(3) # ensure output folder is unique
    out_dir = joinpath(root, "output", "$(now)_$(suffix)")
    grid_search(config, config_path, out_dir)
end
function grid_search(config::Dict, config_path::String, out_dir::String)
    sim_type = get_entry(config["grid_search"], "sim_type", "reference")
    @assert sim_type in ("reference", "validation")

    # Make output folder
    mkpath(out_dir)
    if config_path != joinpath(out_dir, "config.jl")
        cp(config_path, joinpath(out_dir, "config.jl"), force = true)
    end

    # get config entries
    parameters = config["grid_search"]["parameters"]
    n_ens = config["grid_search"]["ensemble_size"]

    # Get cases and count repeat instances of same case names
    cases = config[sim_type]["case_name"]
    case_name_id = Tuple[]
    for (i, case) in enumerate(cases)
        case_id = length(cases[1:i][(cases .== case)[1:i]])
        push!(case_name_id, (case, case_id))
    end
    namelist_args = get_entry(config["scm"], "namelist_args", nothing)

    # Construct simulation configs
    param_names = collect(keys(parameters))
    param_pairs = combinations(param_names, 2)
    sim_configs = SimConfig[]
    for (param1, param2) in param_pairs
        for (case, case_id) in case_name_id
            config_product = Iterators.product(parameters[param1], parameters[param2], 1:n_ens)
            append!(
                sim_configs,
                SimConfig.(
                    param1,
                    param2,
                    getfield.(config_product, 1),  # param1 value
                    getfield.(config_product, 2),  # param2 value
                    getfield.(config_product, 3),  # ensemble id
                    case,
                    case_id,
                    out_dir,
                    Ref(namelist_args),
                ),
            )
        end
    end
    pmap(run_sims, sim_configs, on_error = e -> @warn "Worker failure" exception = (e, catch_backtrace()))  #, on_error = e -> NaN)
end

function parse_commandline_gs()
    s = ArgParseSettings(; description = "Run config input")

    @add_arg_table s begin
        "--config"
        help = "config file"
        arg_type = String
        default = "config.jl"
        "--outdir"
        help = "output directory"
        arg_type = String
        default = ""
        "--mode"
        help = "grid search mode (new or restart)"
        arg_type = String
        default = "new"
    end

    return ArgParse.parse_args(s)
end

if abspath(PROGRAM_FILE) == @__FILE__
    args = parse_commandline_gs()

    sim_mode = args["mode"]
    if sim_mode == "new"
        @info "Starting new grid search"
        config_path = args["config"]
        include(config_path)
        config = get_config()
        grid_search(config, config_path)
    elseif sim_mode == "restart"
        @info "Restarting / continuing existing grid search"
        outdir = args["outdir"]
        config_path = joinpath(outdir, "config.jl")
        include(config_path)
        config = get_config()
        grid_search(config, config_path, outdir)
    else
        throw(ArgumentError("Invalid mode: $sim_mode"))
    end
end


text:
using Metaprogramming
using BaseTestNext

include("primitives.jl")
include("mint.jl")
include("vals.jl")


text:
# Benchmark on heaps

using DataStructures

# benchmark function

function benchmark_heap(title::ASCIIString, h::AbstractHeap, xs::Vector{Float64})
    @assert isempty(h)

    # warming
    push!(h, 0.5)
    pop!(h)

    # bench
    n = length(xs)

    t1 = @elapsed for i = 1 : n
        push!(h, xs[i])
    end
    t2 = @elapsed for i = 1 : n
        pop!(h)
    end

    @printf("   On %-24s:  add.elapsed = %7.4fs  pop.elapsed = %7.4fs\n", title, t1, t2)
end

# Benchmark on add! and pop!

xs = rand(10^6)

h_bin  = binary_minheap(Float64)
h_mbin = mutable_binary_minheap(Float64)

benchmark_heap("BinaryHeap", h_bin, xs)
benchmark_heap("MutableBinaryHeap", h_mbin, xs)


text:
"""
    memory_bandwidth_saxpy_scaling() -> sizes, bandwidths
Measures the memory bandwidth (via [`memory_bandwidth_saxpy`](@ref)) as a function of vector length.
If `verbose=true` (default), displays a unicode plot. Returns the considered lengths and GiB/s.
For further options, see [`memory_bandwidth_saxpy`](@ref).
"""
function memory_bandwidth_saxpy_scaling(;
    device=CUDA.device(), sizes=[2^20 * i for i in 10:10:300], verbose=true, io::IO=stdout, kwargs...
)
    # sizes = [2^20 * i for i in 8:128] # V100
    bandwidths = zeros(length(sizes))
    for (i, s) in enumerate(sizes)
        bandwidths[i] = memory_bandwidth_saxpy(;
            device=device, size=s, verbose=false, kwargs...
        )
        clear_gpu_memory(device)
    end
    if verbose
        peak_val, idx = findmax(bandwidths)
        peak_size = sizes[idx]
        p = UnicodePlots.lineplot(
            sizes,
            bandwidths;
            xlabel="vector length",
            ylabel="GiB/s",
            title=string(
                "Peak: ", round(peak_val; digits=2), " GiB/s (size = $(bytes(peak_size)))"
            ),
            xscale=:log2,
        )
        UnicodePlots.lineplot!(p, [peak_size, peak_size], [0.0, peak_val]; color=:red)
        println(io) # top margin
        println(io,p)
        println(io) # bottom margin
    end
    return sizes, bandwidths
end

"""
Tries to estimate the peak memory bandwidth of a GPU in GiB/s by measuring the time
it takes to perform a SAXPY, i.e. `a * x[i] + y[i]`.

**Keyword arguments:**
* `device` (default: `CUDA.device()`): CUDA device to be used.
* `dtype` (default: `Float32`): element type of the vectors.
* `size` (default: `2^20 * 10`): length of the vectors.
* `nbench` (default: `5`): number of measurements to be performed the best of which is used for the GiB/s computation.
* `verbose` (default: `true`): toggle printing.
* `cublas` (default: `true`): toggle between `CUDA.axpy!` and a custom `saxpy_gpu_kernel!`.
* `io` (default: `stdout`): set the stream where the results should be printed.

See also: [`memory_bandwidth_saxpy_scaling`](@ref).
"""
function memory_bandwidth_saxpy(;
    device=CUDA.device(),
    size=2^20 * 10,
    nbench=10,
    dtype=Float32,
    cublas=true,
    verbose=true,
    io::IO=stdout
)
    device!(device) do
        a = dtype(pi)
        x = CUDA.rand(dtype, size)
        y = CUDA.rand(dtype, size)
        z = CUDA.zeros(dtype, size)

        nthreads = CUDA.attribute(device, CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
        nblocks = cld(size, nthreads)
        t = Inf
        for _ in 1:nbench
            if cublas
                Δt = CUDA.@elapsed CUBLAS.axpy!(size, a, x, y)
            else
                Δt = CUDA.@elapsed @cuda(
                    threads = nthreads, blocks = nblocks, saxpy_gpu_kernel!(z, a, x, y)
                )
            end
            t = min(t, Δt)
        end

        bandwidth = 3.0 * sizeof(dtype) * size * (1024)^(-3) / t
        if verbose
            printstyled(io,"Memory Bandwidth (GiB/s):\n"; bold=true)
            print(io," └ max: ")
            printstyled(io,round(bandwidth; digits=2), "\n"; color=:green, bold=true)
        end
        return bandwidth
    end
end

function saxpy_gpu_kernel!(z, a, x, y)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(z)
        @inbounds z[i] = a * x[i] + y[i]
    end
    return nothing
end


text:
########## Dumbo ##########

mutable struct Dumbo
    energy::Int
    neighbors::Vector{Dumbo}
    has_flashed::Bool
end
energy(dumbo::Dumbo) = dumbo.energy
has_flashed(dumbo::Dumbo) = dumbo.has_flashed
neighbors(dumbo::Dumbo) = dumbo.neighbors

Dumbo(energy::Int) = Dumbo(energy, Vector{Dumbo}[], false)

is_ready(dumbo::Dumbo) = !has_flashed(dumbo) && energy(dumbo) > 9

energize!(dumbo::Dumbo) = dumbo.energy += 1

function reset!(dumbo::Dumbo)
    if has_flashed(dumbo)
        dumbo.energy = 0
        dumbo.has_flashed = false
    end
end

function flash!(dumbo::Dumbo)
    dumbo.has_flashed = true
    energize!.(neighbors(dumbo))
    for neighbor in neighbors(dumbo)
        is_ready(neighbor) && flash!(neighbor)
    end
end

function set_neighbors!(dumbos::Matrix{Dumbo})
    indices = Tuple.(CartesianIndices(dumbos))
    for (x, y) in indices
        # 3 by 3 matrix centered on the current position but excluding said position
        # as well as positions outside the reference matrix
        neighbors_indices = indices ∩
            [(x2, y2) for x2 = x .+ (-1:1), y2 = y .+ (-1:1) if (x2, y2) != (x, y)]

        dumbos[x, y].neighbors = [dumbos[x2, y2] for (x2, y2) in neighbors_indices]
    end
end

function step!(dumbos::Matrix{Dumbo})::Tuple{Int, Bool}
    energize!.(dumbos)
    for dumbo in dumbos
        is_ready(dumbo) && flash!(dumbo)
    end
    flashes = has_flashed.(dumbos)
    n_flashes = sum(flashes)
    synchronized = all(flashes)
    reset!.(dumbos)
    return n_flashes, synchronized
end

########## Script ##########

function main(input_name::String)

    # open file and store content into an input string
    input = open(join(["inputs//", input_name]), "r") do io
        read(io, String)
    end

    # split input into a data vector of vector of integers
    data = [parse.(Int, collect(row))
        for row in String.(split(input, "\r\n"))]

    # parse data into a matrix of dumbos and set neighbors
    dumbos = Dumbo.([data[i][j]
        for i = 1:length(data), j = 1:length(data[1])])
    set_neighbors!(dumbos)

    # a bit of complexity to make sure part 2 does'nt end the loop
    # before part 1 is done and vice-versa (unecessary for the input
    # I got but I pretend I don't know that)
    total_flashes = 0
    synchronized = false
    has_synchronized = false
    step_sync = Inf
    step = 0
    while !synchronized || step <= 100
        step += 1
        n_flashes, synchronized = step!(dumbos)

        # part 1
        step <= 100 && (total_flashes += n_flashes)

        # part 2
        if synchronized && !has_synchronized
            step_sync = step
            has_synchronized = true
        end
    end
    println("Part 1: ", total_flashes)
    println("Part 2: ", step_sync)

end

main("day11.txt");


text:
"""
Treat WAV files as lazy arrays stored on disk.

Documentation: https://github.com/baggepinnen/LazyWAVFiles.jl

## Quick start
### LazyWAVFile
```
f1 = LazyWAVFile(joinpath(d,"f1.wav"))
f1[1]
f1[1:5]
size(f1)
f1.fs # Get the sample rate

[f1; f2] # Creates a `DistributedWAVFile`
```

### DistributedWAVFile
```
df = DistributedWAVFile(folder_path)
df[1]
df[1:12]
df[:]

size(df) # Other array functions are defined as well
length(df)
df.fs    # Get the sample rate
```
"""
module LazyWAVFiles
using WAV, BlockArrays
export LazyWAVFile, DistributedWAVFile, path

struct LazyWAVFile{T,N,FS,S<:Tuple} <: AbstractArray{T,N}
    path::String
    size::S
    fs::FS
end
function LazyWAVFile(path)
    r,fs = wavread(path, format="native", subrange=1)
    s = wavread(path, format="size")
    T = eltype(r)
    dim = s[2]
    if dim == 1
        s = (s[1],)
    end
    LazyWAVFile{T,length(s),typeof(fs),typeof(s)}(path, s, fs)
end

Base.size(f::LazyWAVFile) = f.size
Base.size(f::LazyWAVFile{T,N},i) where {T,N} = i > N ? 1 : f.size[i]
Base.length(f::LazyWAVFile) = prod(f.size)
path(f::LazyWAVFile) = f.path

function Base.copyto!(dst::SubArray{T,N,Array{T,N},Tuple{UnitRange{Int64}}}, sa::SubArray{T,N,<:LazyWAVFile{T,N}}) where {T,N}
    copyto!(dst, sa.parent[sa.indices...])
end

function Base.copyto!(dst::AbstractArray, rd::Int, src::LazyWAVFile, rs::Int, N::Int)
    dst[rd:rd+N-1] .= src[rs:rs+N-1]
end

Base.getindex(f::LazyWAVFile{T,N}, i::Integer) where {T,N} = wavread(f.path, format="native", subrange=i:i)[1][1]::T

Base.getindex(f::LazyWAVFile{T,N}, i::Integer,j::Integer) where {T,N} = wavread(f.path, format="native", subrange=i:i)[1][j]::T
Base.getindex(f::LazyWAVFile{T,N}, i::Integer,j::Union{Colon,AbstractRange}) where {T,N} = wavread(f.path, format="native", subrange=i:i)[1][1,j]::Array{T,1}
Base.getindex(f::LazyWAVFile{T,N}, i::Union{Colon,AbstractRange},j::Integer) where {T,N} = wavread(f.path, format="native", subrange=i)[1][:,j]::Array{T,1}
Base.getindex(f::LazyWAVFile{T,N}, i,j) where {T,N} = wavread(f.path, format="native", subrange=i)[1][:,j]
Base.getindex(f::LazyWAVFile{T,1}, i::AbstractRange) where {T} = vec(wavread(f.path, format="native", subrange=i)[1])::Array{T,1}
Base.getindex(f::LazyWAVFile{T,N}, i::AbstractRange) where {T,N} = wavread(f.path, format="native", subrange=i)[1]::Array{T,2}

Base.getindex(f::LazyWAVFile{T,1}, ::Colon) where {T} = vec(wavread(f.path, format="native")[1])::Array{T,1}
Base.getindex(f::LazyWAVFile{T,N}, ::Colon, ::Colon) where {T,N} = wavread(f.path, format="native")[1]::Array{T,2}

Base.eltype(f::LazyWAVFile{T}) where T = T
Base.ndims(f::LazyWAVFile{T,N}) where {T,N} = N
Base.show(io::IO, ::MIME"text/plain", f::LazyWAVFile{T,N,S}) where {T,N,S} = println(io, "LazyWAV{$T, $N, $(f.size), fs=$(f.fs)}: ", f.path)

Base.view(f::LazyWAVFile, i...) = f[i...]

struct DistributedWAVFile{T,N,L,FS} <: AbstractArray{T,N}
    files::Vector{LazyWAVFile{T,N,FS}}
    blockarray::L
    fs::FS
end
function DistributedWAVFile(wavpaths::Vector{String})
    isempty(wavpaths) && error("Found no wav paths.")
    files = LazyWAVFile.(wavpaths)
    fs0 = files[1].fs
    all(x->x.fs==fs0, files) || error("WAV files have different sample rates.")
    DistributedWAVFile(files, fs0)
end
function DistributedWAVFile(folder::String)
    files = filter(readdir(folder)) do file
        lowercase(file[end-2:end]) == "wav"
    end
    isempty(files) && error("Found no wav files in the specified path: $folder")
    files = sort(files)
    files = LazyWAVFile.(joinpath.(Ref(folder), files))
    fs0 = files[1].fs
    if !all(x->x.fs==fs0, files)
        error("WAV files in $folder have different sample rates.")
    end
    DistributedWAVFile(files,fs0)
end
function DistributedWAVFile(files::AbstractVector{L},fs) where L <: LazyWAVFile{T,N} where {T,N}
    blockarray = N == 1 ? mortar(files) : mortar(reshape(files, (length(files), 1)))
    try
        return DistributedWAVFile{eltype(files[1]), ndims(files[1]), typeof(blockarray), typeof(fs)}(files, blockarray, fs)
    catch e
        @error "Creating distributed WAV file failed. This can happen if the wav-files in the folder have different number of channels."
        rethrow(e)
    end
end
Base.length(f::DistributedWAVFile) = sum(length, f.files)
Base.size(f::DistributedWAVFile{T,N}) where {T,N} = size(f.blockarray)

Base.show(io::IO, ::MIME"text/plain", f::DistributedWAVFile{T,N}) where {T,N} = println(io, "DistributedWAVFile{$T, $N} with $(length(f.files)) files, $(length(f)) total datapoints and samplerate $(f.fs)")

function blockindexranges(bindices::Vector{BlockIndex{1}})
    indices = first.((x -> x.I).(bindices))
    blks = unique(indices)
    biranges = BlockArrays.BlockIndexRange{1,Tuple{UnitRange{Int}}}[]
    for blk in blks
        bindicesindices = (indices .== blk)
        push!(
            biranges,
            Block(blk)[bindices[bindicesindices][1].α[1]:bindices[bindicesindices][end].α[1]],
        )
    end
    biranges
end

function Base.getindex(blockarray::BlockVector{T}, indices::UnitRange{Int}) where T
    bindices = findblockindex.(axes(blockarray), indices)
    biranges = blockindexranges(bindices)
    x = Vector{eltype(blockarray)}(undef, length(indices))
    startind = 1
    for birange in biranges
        stopind = startind+length(birange)-1
        x[startind:stopind] = blockarray[birange]
        startind = stopind+1
    end
    x
end

Base.getindex(df::DistributedWAVFile, i...) = getindex(df.blockarray, i...)

function Base.vcat(lfs::LazyWAVFile...)
    fs = lfs[1].fs
    any(x->x.fs != fs, lfs) && error("Concatenated WAV files have different sample rates.")
    DistributedWAVFile([lfs...], fs)
end

function Base.vcat(dfs::DistributedWAVFile...)
    fs = dfs[1].fs
    any(x->x.fs != fs, dfs) && error("Distributed WAV files have different sample rates.")
    DistributedWAVFile(reduce(vcat, getfield.(dfs, :files)), fs)
end

Base.view(df::DistributedWAVFile, i...) = df[i...]

end


text:
using Plots
result_triangles=[[[0.0;0.0], [0.5;sqrt(3)/2],[1.0;0.0],[0.0;0.0]]]
for step in 1:5
    newOnes=[]
    for triangle in result_triangles
        first_triangle=[]
        second_triangle=[]
        third_triangle=[]
        for i in 1:4
            vertix=(triangle[i]-triangle[1])/2+triangle[1]
            push!(first_triangle, vertix)
        end
        for i in 1:4
            push!(second_triangle, first_triangle[i]+(first_triangle[2]-first_triangle[1]))
            push!(third_triangle, first_triangle[i]+(first_triangle[3]-first_triangle[1]))
        end
        push!(newOnes, first_triangle)
        push!(newOnes, second_triangle)
        push!(newOnes, third_triangle)
        #println(result_triangles)
    end
    result_triangles=newOnes
    #println(result_triangles)
end
plot()
for little_triangles in result_triangles
    #println(little_triangles)
    plot!(hcat(little_triangles...)[1,:] , hcat(little_triangles...)[2,:], fill=(0,0.5,"orange"), legend=false)
end
plot!()
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter2\\Fig\\SierpinskiTriangle.png")
# plot(hcat(first_triangle...)[1,:],hcat(first_triangle...)[2,:])
# plot!(hcat(second_triangle...)[1,:],hcat(second_triangle...)[2,:])
# plot!(hcat(third_triangle...)[1,:],hcat(third_triangle...)[2,:])


text:
module SubModule
using Test

RUN = []
@testset "A" begin
    push!(RUN, 1)
end

include("subsubmodule.jl")

end # SubModule

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
"""
   get the observations for multiple single source, parallel computation is implemented,
crop the physical model adaptively according to acquisition geometry.
"""
function get_shotgather_adaptive(dir_obs::Ts, isz, isx, wlet, irz::Ti, irx::Ti,
                                 vel::Matrix{Tv}, rho::Matrix{Tv},  dz, dx, dt, tmax;
                                 location_flag="index",  data_format=Float32, order=2, free_surface=true, npad=100) where {Ts<:String, Ti<:Vector, Tv<:Real}

    # define function inside
    function wrap_internal(params::NamedTuple)

        # get the original model size
        (nz, nx) = size(params.vel)

        # determine model range in spatial direction
        ix_lower = minimum(vcat(params.source_x, params.receiver_x))
        if ix_lower > params.npad
           ix_lower = ix_lower - params.npad
        else
           ix_lower = 1
        end

        # upper end
        ix_upper = maximum(vcat(params.source_x, params.receiver_x))
        if ix_upper + params.npad > nx
           ix_upper = nx
        else
           ix_upper = ix_upper + params.npad
        end

        # crop part of velocity model
        vel = params.vel[:, ix_lower:ix_upper]
        rho = params.rho[:, ix_lower:ix_upper]
        fidiff = TdParams(rho, vel, params.free_surface, params.dz, params.dx, params.dt, params.tmax;
                          data_format=params.data_format, order=params.order)

        # adjust receiver's location
        receiver_x = copy(params.receiver_x)
        for i = 1 : length(receiver_x)
            receiver_x[i] = receiver_x[i] - ix_lower + 1
        end

        # adjust source's location
        source_x = params.source_x - ix_lower + 1
        src = Source(params.source_z, source_x, fidiff; p=params.wlet)

        # do simulation
        # rec = multi_step_forward!(params.receiver_z, receiver_x, src, fidiff; location_flag=params.location_flag)
        multi_step_forward!(params.receiver_z, receiver_x, src, fidiff; location_flag=params.location_flag, path_shot=params.path_obs)

        # save the result to disk
        # write_recordings(params.path_obs, rec)

        return nothing
    end

    # create a fold save observations
    if isdir(dir_obs)
       rm(dir_obs, force=true, recursive=true)
    end
    mkdir(dir_obs)
    if !isdir(dir_obs)
       error("can't create directory for reflections")
    end

    # determine number of shot
    ns =  length(isx)
    ns == length(isz) || error("source coordinate doesn't match")

    # wrap parameters into vector of named tuple
    argument_collection = Vector{NamedTuple}(undef, ns)
    for i = 1 : ns

        file_name = join(["recordings_" "$i" ".bin"])
        path_obs  = joinpath(dir_obs, file_name)

        # OBN acquisition geometry
        if eltype(irz) <: Real
           argument_collection[i] = (path_obs=path_obs, receiver_z=irz, receiver_x=irx, location_flag=location_flag, source_z=isz[i], source_x=isx[i], wlet=wlet,
                                     free_surface=free_surface, dz=dz, dx=dx, dt=dt, tmax=tmax, data_format=data_format, order=order, vel=vel, rho=rho, npad=npad)

        # towed streamer
        elseif eltype(irz) <: Vector
           argument_collection[i] = (path_obs=path_obs, receiver_z=irz[i], receiver_x=irx[i], location_flag=location_flag, source_z=isz[i], source_x=isx[i], wlet=wlet,
                                     free_surface=free_surface, dz=dz, dx=dx, dt=dt, tmax=tmax, data_format=data_format, order=order, vel=vel, rho=rho, npad=npad)
        else
           error("wrong type receiver locations")
        end
    end

    # do simulation parallel
    if nprocs() == 1
       for i = 1 : ns
           wrap_internal(argument_collection[i])
       end
    else
       pmap(wrap_internal, argument_collection)
    end

    return nothing
end

"""
   get the observations for multiple single source, parallel computation is implemented
"""
function get_shotgather(dir_obs::Ts, irz::Ti, irx::Ti, src::T, fidiff::P;
                         location_flag="index") where {Ts<:String, Ti<:Vector, T<:Union{Source, Vector{Source}}, P<:TdParams}

    # define function accept named tuple
    function wrap_internal(params::NamedTuple)

        # do simulation
        multi_step_forward(params.source, params.fidiff; rz=params.receiver_z, rx=params.receiver_x,
                           location_flag=params.location_flag, path_shot=params.path_obs)
    end

    # create a fold save observations
    if isdir(dir_obs)
       rm(dir_obs, force=true, recursive=true)
    end
    mkdir(dir_obs)
    if !isdir(dir_obs)
       error("can't create directory for reflections")
    end

    # determine number of shot
    if typeof(src) <: Source
       ns  = 1
       src = [src]
    else
       ns = length(src)
    end

    # wrap parameters into vector of named tuple
    argument_collection = Vector{NamedTuple}(undef, ns)
    for i = 1 : ns

        file_name = join(["recordings_" "$i" ".bin"])
        path_obs  = joinpath(dir_obs, file_name)

        # OBN acquisition geometry
        if eltype(irz) <: Real
           argument_collection[i] = (path_obs=path_obs, receiver_z=irz, receiver_x=irx, location_flag=location_flag,
                                     source=src[i], fidiff=fidiff)

        # towed streamer
        elseif eltype(irz) <: Vector
           argument_collection[i] = (path_obs=path_obs, receiver_z=irz[i], receiver_x=irx[i], location_flag=location_flag,
                                     source=src[i], fidiff=fidiff)
        else
           error("wrong type receiver locations")
        end
    end

    # do simulation parallel
    if nprocs() == 1
       for i = 1 : ns
           wrap_internal(argument_collection[i])
       end
    else
       pmap(wrap_internal, argument_collection)
    end

    return nothing
end

"""
   get the reflections for multiple single source, parallel computation is implemented
"""
function get_reflections(dir_obs::Ts, irz::Ti, irx::Ti, src::T, params_hete::P, params_homo::P;
                         location_flag="index") where {Ts<:String, Ti<:Vector, T<:Union{Source, Vector{Source}}, P<:TdParams}

    # define function accept named tuple
    function wrap_get_reflections(params::NamedTuple)

        # do simulation
        dobs = multi_step_forward(params.source, params.fidiff_hete; rz=params.receiver_z, rx=params.receiver_x, location_flag=params.location_flag)
        dire = multi_step_forward(params.source, params.fidiff_homo; rz=params.receiver_z, rx=params.receiver_x, location_flag=params.location_flag)

        # remove direct wave
        dobs.p .= dobs.p .- dire.p

        # save the reflections to disk
        write_recordings(params.path_obs, dobs)
    end

    # create a fold save observations
    rm(dir_obs, force=true, recursive=true)
    mkdir(dir_obs)
    if !isdir(dir_obs)
       error("can't create directory for reflections")
    end

    # determine number of shot
    if typeof(src) <: Source
       ns  = 1
       src = [src]
    else
       ns = length(src)
    end

    # wrap parameters into vector of named tuple
    argument_collection = Vector{NamedTuple}(undef, ns)
    for i = 1 : ns

        file_name = join(["recordings_" "$i" ".bin"])
        path_obs  = joinpath(dir_obs, file_name)

        # OBN acquisition geometry
        if eltype(irz) <: Real
           argument_collection[i] = (path_obs=path_obs, receiver_z=irz, receiver_x=irx, location_flag=location_flag,
                                     source=src[i], fidiff_hete=params_hete, fidiff_homo=params_homo)

        # towed streamer
        elseif eltype(irz) <: Vector
           argument_collection[i] = (path_obs=path_obs, receiver_z=irz[i], receiver_x=irx[i], location_flag=location_flag,
                                     source=src[i], fidiff_hete=params_hete, fidiff_homo=params_homo)
        else
           error("wrong type receiver locations")
        end
    end

    # do simulation parallel
    if nprocs() == 1
       for i = 1 : ns
           wrap_get_reflections(argument_collection[i])
       end
    else
       pmap(wrap_get_reflections, argument_collection)
    end

    return nothing
end

"""
   get wavefield boundary, last wavefield and sourceside wavefield strength
"""
function get_wavefield_bound(dir_sourceside::Ts, src::T, fidiff::P; remove_flag=true) where {Ts<:String, T<:Union{Source, Vector{Source}}, P<:TdParams}

    # define function accept named tuple
    function wrap_get_wavefield_bound(params::NamedTuple)

        # do simulation
        multi_step_forward(params.source, params.fidiff;
                           path_bnd=params.path_bnd, path_lwfd=params.path_wfd, path_sws=params.path_sws)

        # save the source wavelet
        write_source(params.path_src, params.source)

        return nothing
    end

    # create a folder for sourceside wavefield
    rm(dir_sourceside, force=true, recursive=true)
    mkdir(dir_sourceside)
    if !isdir(dir_sourceside) # check the directory is created
       error("can't create directory for source-side wavefield")
    end

    # create folders to save the file
    dir_bnd = joinpath(dir_sourceside, "boundary")
    dir_wfd = joinpath(dir_sourceside, "wavefield")
    dir_sws = joinpath(dir_sourceside, "strength")
    dir_src = joinpath(dir_sourceside, "source")

    # if this folder doesn't exist, create one
    mkdir(dir_bnd)
    if !isdir(dir_bnd) # check the directory is created
       error("can't create directory for boundary")
    end

    mkdir(dir_wfd)
    if !isdir(dir_wfd) # check the directory is created
       error("can't create directory for wavefield")
    end

    mkdir(dir_sws)
    if !isdir(dir_sws) # check the directory is created
       error("can't create directory for source strength")
    end

    mkdir(dir_src)
    if !isdir(dir_src) # check the directory is created
       error("can't create directory for source wavelet")
    end

    # determine number of shot
    if typeof(src) <: Source
       ns  = 1
       src = [src]
    else
       ns = length(src)
    end

    # wrap parameters into vector of named tuple
    argument_collection = Vector{NamedTuple}(undef, ns)
    for i = 1 : ns

        file_name = join(["boundary_" "$i" ".rsf"])
        path_bnd  = joinpath(dir_bnd, file_name)

        file_name = join(["wavefield_" "$i" ".rsf"])
        path_wfd  = joinpath(dir_wfd, file_name)

        file_name = join(["strength_" "$i" ".rsf"])
        path_sws  = joinpath(dir_sws, file_name)

        file_name = join(["source_" "$i" ".bin"])
        path_src  = joinpath(dir_src, file_name)

        # pack parameter for source-side wavefield
        argument_collection[i] = (source=src[i], fidiff=fidiff, path_bnd=path_bnd,
                                  path_wfd=path_wfd, path_sws=path_sws, path_src=path_src)
    end

    # do simulation parallel
    if nprocs() == 1
       for i = 1 : ns
           wrap_get_wavefield_bound(argument_collection[i])
       end
    else
       pmap(wrap_get_wavefield_bound, argument_collection)
    end

    # compute normalization factor
    nal = zeros(fidiff.data_format, fidiff.nz, fidiff.nx)
    for i = 1 : ns
        (hdr, sws) = read_RSdata(argument_collection[i].path_sws)
        nal .= nal .+ sws
    end
    nal .= 1.0 ./ nal

    # remove the source strength file
    remove_flag && rm(dir_sws, force=true, recursive=true)

    # save the normalization factor
    hdr = RegularSampleHeader(nal; title="source normalization")
    path_nal = joinpath(dir_sourceside, "normalization.rsf")
    write_RSdata(path_nal, hdr, nal)

    return nothing
end

"""
   create a project folder for LSRTM
"""
function initialize_lsrtm(work_dir::String)

    # check physical model
    dir_physical = joinpath(work_dir, "physical_model")
    isdir(dir_physical) || error("no physical model")

    dir_born = joinpath(work_dir, "data_space")
    rm(dir_born, force=true, recursive=true)
    mkdir(dir_born)
    if !isdir(dir_born) # check the directory is created
       error("can't create directory for variables in data space")
    end

    dir_model = joinpath(work_dir, "model_space")
    rm(dir_model, force=true, recursive=true)
    mkdir(dir_model)
    if !isdir(dir_model) # check the directory is created
       error("can't create directory for variables in model space")
    end

    # dir_sourceside = joinpath(work_dir, "sourceside")
    # rm(dir_sourceside, force=true, recursive=true)
    # mkdir(dir_sourceside)
    # if !isdir(dir_sourceside) # check the directory is created
    #    error("can't create directory for boundary values")
    # end

    return nothing
end