text:
# Greetings! from path ~/.julia/config/startup.jl
using SeisIO, SeisNoise, Dates, Glob, Base.Threads, Distributed
export FFT_CC

function get_alltime_path_fftcc(FFTDIR::String, select_time::Bool, TIMETXT::String)
    if select_time == true
        alltime_path = readlines(TIMETXT)
    else
        alltime_path = sort(glob("*/*/*", FFTDIR))
    end
    return alltime_path
end

function get_SR_pairs_fftcc(select_sta::String, STATXT::String, SRTXT::String)
    All_NET_STA = Vector{Vector{String}}(undef, 0)
    SR_pairs = Vector{Vector{String}}(undef, 0)
    if select_sta == "no"
        nothing
    elseif select_sta == "sta"
        # a.read STATXT --> station_txt
        station_txt = readlines(STATXT)  # read sta.txt
        # b.make all stations --> rAll_NET_STA (for reading data[sac/mseed] from database and auto-corr)
        N_sta = size(station_txt, 1)
        rAll_NET_STA = Vector{Vector{String}}(undef, N_sta)
        for ll = 1:N_sta
            D_net, D_sta = split(station_txt[ll])
            rAll_NET_STA[ll] = [string(D_net), string(D_sta)]
        end
        # c.make s-r pairs --> rSR_pairs (for cross-corr)
        N_SR = Int(N_sta * (N_sta - 1) / 2)
        rSR_pairs = Vector{Vector{String}}(undef, N_SR)
        num = 0
        for ss = 1:N_sta-1
            for rr = ss+1:N_sta
                num += 1
                S_net, S_sta = split(station_txt[ss])
                R_net, R_sta = split(station_txt[rr])
                rSR_pairs[num] =
                    [string(S_net), string(S_sta), string(R_net), string(R_sta)]
            end
        end
        # d.return
        All_NET_STA = rAll_NET_STA
        SR_pairs = rSR_pairs
    elseif select_sta == "sr"
        # a.read SRTXT --> station_txt
        sr_txt = readlines(SRTXT)
        # b.generate pSR_pairs with repeated sta-pairs --> pSR_pairs
        source_mark = 1   # "1" mark the source of one read-line
        S_net = ""
        S_sta = ""
        pSR_pairs = Vector{Vector{String}}(undef, 0)
        for i = 1:size(sr_txt, 1)
            if sr_txt[i] != ""
                if source_mark == 1
                    s_net, s_sta = split(sr_txt[i])
                    S_net = s_net
                    S_sta = s_sta
                    source_mark = -1   # "-1" mark next line is not the source
                else
                    R_net, R_sta = split(sr_txt[i])
                    if !((S_net == R_net) && (S_sta == R_sta))
                        push!(
                            pSR_pairs,
                            [string(S_net), string(S_sta), string(R_net), string(R_sta)],
                        )
                    end
                end
            else
                source_mark = 1   # mark next line is the source
            end
        end
        # c.delete the repeated pSR_pairs --> SR_pairs
        rSR_pairs = Vector{Vector{String}}(undef, 0)
        repeat_mark = -1   # "-1" means not repeat
        for i = 1:size(pSR_pairs, 1)-1
            S_net_1, S_sta_1, R_net_1, R_sta_1 =
                pSR_pairs[i][1], pSR_pairs[i][2], pSR_pairs[i][3], pSR_pairs[i][4]
            for j = i+1:size(pSR_pairs, 1)
                S_net_2, S_sta_2, R_net_2, R_sta_2 =
                    pSR_pairs[j][1], pSR_pairs[j][2], pSR_pairs[j][3], pSR_pairs[j][4]
                condition_1 =
                    (S_net_1 == S_net_2) &&
                    (S_sta_1 == S_sta_2) &&
                    (R_net_1 == R_net_2) &&
                    (R_sta_1 == R_sta_2)   # Completely corresponding
                condition_2 =
                    (S_net_1 == R_net_2) &&
                    (S_sta_1 == R_sta_2) &&
                    (R_net_1 == S_net_2) &&
                    (R_sta_1 == S_sta_2)   # Reverse S and R, then corresponde
                if condition_1 || condition_2
                    repeat_mark = 1
                    break
                else
                    repeat_mark = -1
                end
            end
            if repeat_mark != 1
                push!(
                    rSR_pairs,
                    [string(S_net_1), string(S_sta_1), string(R_net_1), string(R_sta_1)],
                )
            end
        end
        # The last group will never repeat
        S_net_1, S_sta_1, R_net_1, R_sta_1 =
            pSR_pairs[end][1], pSR_pairs[end][2], pSR_pairs[end][3], pSR_pairs[end][4]
        push!(
            rSR_pairs,
            [string(S_net_1), string(S_sta_1), string(R_net_1), string(R_sta_1)],
        )

        # d.make all stations --> All_NET_STA (for read data[sac/mseed] from FFT-database and auto-corr)
        double_SR_pairs = Vector{Vector{String}}(undef, 0)
        for i = 1:size(rSR_pairs, 1)
            S_net, S_sta, R_net, R_sta =
                rSR_pairs[i][1], rSR_pairs[i][2], rSR_pairs[i][3], rSR_pairs[i][4]
            push!(double_SR_pairs, [S_net, S_sta])
            push!(double_SR_pairs, [R_net, R_sta])
        end
        rAll_NET_STA = Vector{Vector{String}}(undef, 0)
        NET_1 = ""
        STA_1 = ""
        repeat_mark = -1
        for i = 1:size(double_SR_pairs, 1)-1
            net_1, sta_1 = double_SR_pairs[i][1], double_SR_pairs[i][2]
            NET_1, STA_1 = net_1, sta_1
            for j = i+1:size(double_SR_pairs, 1)
                net_2, sta_2 = double_SR_pairs[j][1], double_SR_pairs[j][2]
                if NET_1 == net_2 && STA_1 == sta_2
                    repeat_mark = 1
                    break
                else
                    repeat_mark = -1
                end
            end
            if repeat_mark != 1
                push!(rAll_NET_STA, [string(NET_1), string(STA_1)])
            end
        end
        NET_1, STA_1 = double_SR_pairs[end][1], double_SR_pairs[end][2]
        push!(rAll_NET_STA, [string(NET_1), string(STA_1)])
        # e.conclusion --> All_NET_STA,SR_pairs
        All_NET_STA = rAll_NET_STA
        SR_pairs = rSR_pairs
    end

    ### modify sequence of SR_pairs
    net_1,sta_1,net_2,sta_2 = "","","",""
    for ii = 1:size(SR_pairs,1)
        net_1,sta_1,net_2,sta_2 = SR_pairs[ii][1],SR_pairs[ii][2],SR_pairs[ii][3],SR_pairs[ii][4]
        if net_1*sta_1 > net_2*sta_2
            SR_pairs[ii][1] = net_2
            SR_pairs[ii][2] = sta_2
            SR_pairs[ii][3] = net_1
            SR_pairs[ii][4] = sta_1
        end
    end

    return sort(All_NET_STA), sort(SR_pairs)
end

function get_allfile_list_fftcc(Onetime_path::String, channel_regular::String)
    files = Vector{String}(undef, 0)
    files_cha = glob(channel_regular, Onetime_path)
    append!(files, files_cha)
    sort!(files)
    return files
end

function get_S_dict_fftcc(S::SeisData)
    S_dict = Dict{String,SeisChannel}()
    for i = 1:S.n
        net, sta, kel, cha = split(S[i].id, '.')
        key = string(net, '.', sta)
        setindex!(S_dict, S[i], key)
    end
    return S_dict
end

function get_SeisData_fftcc(
    Onetime_path::String,
    channel_regular::String,
    select_sta::String,
    STATXT::String,
    SRTXT::String,
    input_fmt::String,
)
    # 1.read one time-chunck paths
    files = get_allfile_list(Onetime_path, channel_regular)
    # 2.read data
    N_files = size(files, 1)
    All_NET_STA = Vector{Vector{String}}(undef, 0)
    SR_pairs = Vector{Vector{String}}(undef, 0)
    S = SeisData()
    if select_sta == "no"
        # a.read all data --> S::SeisData
        for j = 1:N_files
            file = files[j]
            S_one = ""
            try
                s_one = read_data(input_fmt, file)
                S_one = deepcopy(s_one)
            catch e
                println("Errorpass in function 'get_SeisData':", e)
                continue
            end
            append!(S, S_one)
        end
        if S.n == 0   # one time-chunck include 0 files
            return S, 0, 0
        end
        # b.make all stations --> All_NET_STA (for read data[sac/mseed] from database and auto-corr)
        N_sta = S.n
        rAll_NET_STA = Vector{Vector{String}}(undef, 0)
        for j = 1:N_sta
            net, sta, kno, cha = split(S[j].id, '.')
            push!(rAll_NET_STA, [string(net), string(sta)])
        end
        All_NET_STA = deepcopy(rAll_NET_STA)
        # c.make s-r pairs --> SR_pairs (for cross-corr)
        N_SR = Int(N_sta * (N_sta - 1) / 2)
        rSR_pairs = Vector{Vector{String}}(undef, N_SR)
        num = 0
        for ss = 1:N_sta-1
            for rr = ss+1:N_sta
                num += 1
                S_net, S_sta = rAll_NET_STA[ss]
                R_net, R_sta = rAll_NET_STA[rr]
                rSR_pairs[num] = [S_net, S_sta, R_net, R_sta]
            end
        end
        SR_pairs = deepcopy(rSR_pairs)
    else
        # a.read sta-pairs --> rAll_NET_STA, rSR_pairs
        rAll_NET_STA, rSR_pairs = get_SR_pairs(select_sta, STATXT, SRTXT)
        # b.read data from path with selected 'All_NET_STA' --> error_net_sta, S
        error_net_sta = Vector{Vector{String}}(undef, 0)
        read_mark = 1
        for ii = 1:size(rAll_NET_STA, 1)
            for j = 1:N_files
                file = files[j]
                if occursin(rAll_NET_STA[ii][1], basename(file)) &&
                   occursin(rAll_NET_STA[ii][2], basename(file))
                    read_mark = 1
                    try
                        S_one = read_data(input_fmt, file)
                        append!(S, S_one)
                        break
                    catch e
                        push!(error_net_sta, [rAll_NET_STA[ii][1], rAll_NET_STA[ii][2]])   # read error
                        println("Errorpass in function 'get_SeisData':", e)
                        continue
                    end
                else
                    read_mark = -1
                end
            end
            if read_mark == -1   # could not find file with net-sta provided by sta.txt
                push!(error_net_sta, [rAll_NET_STA[ii][1], rAll_NET_STA[ii][2]])
            end
        end
        if S.n == 0   # one time-chunck include 0 files
            return S, 0, 0
        end
        # c.modify rAll_NET_STA, rSR_pairs --> All_NET_STA, SR_pairs
        error_num = length(error_net_sta)
        if error_num == 0
            All_NET_STA = rAll_NET_STA
            SR_pairs = rSR_pairs
        else
            # modify All_NET_STA
            for i = 1:size(rAll_NET_STA, 1)
                net_1, sta_1 = rAll_NET_STA[i][1], rAll_NET_STA[i][2]
                mark_sta = -1
                for j = 1:error_num
                    net_2, sta_2 = error_net_sta[j][1], error_net_sta[j][2]
                    if net_1 == net_2 && sta_1 == sta_2
                        mark_sta = 1
                        break
                    end
                end
                if mark_sta == -1
                    push!(All_NET_STA, [net_1, sta_1])
                end
            end
            # modify SR_pairs
            for i = 1:size(rSR_pairs, 1)
                s_net_1, s_sta_1, r_net_1, r_sta_1 =
                    rSR_pairs[i][1], rSR_pairs[i][2], rSR_pairs[i][3], rSR_pairs[i][4]
                mark_sta = -1
                for j = 1:error_num
                    net_2, sta_2 = error_net_sta[j][1], error_net_sta[j][2]
                    if (s_net_1 == net_2 && s_sta_1 == sta_2) ||
                        (r_net_1 == net_2 && r_sta_1 == sta_2)
                        mark_sta = 1
                        break
                    end
                end
                if mark_sta == -1
                    push!(SR_pairs, [s_net_1, s_sta_1, r_net_1, r_sta_1])
                end
            end
        end
    end
    # 3.convert S --> S_dict
    S_dict = get_S_dict(S)

    return S_dict, All_NET_STA, SR_pairs
end

function get_SeisData2RawData_fftcc(S_dict::Dict{String,SeisChannel}, cc_len, cc_step)
    R_dict = Dict{String,RawData}()
    for (k, v) in S_dict
        R = RawData(v, cc_len, cc_step)
        demean!(R)
        taper!(R)
        setindex!(R_dict, R, k)
    end
    return R_dict
end

function get_time_norm_fftcc(
    R_dict::Dict{String,RawData},
    time_norm::String,
    factor_clip_std::Float64,
    factor_mute::Float64,
    time_half_win::Int64,
)
    if time_norm == "no"
        nothing
    elseif time_norm == "one_bit"
        for (k, _) in R_dict
            onebit!(R_dict[k])
        end
    elseif time_norm == "clip"
        for (k, _) in R_dict
            clip!(R_dict[k], factor_clip_std)
        end
    elseif time_norm == "mute"
        for (k, _) in R_dict
            mute!(R_dict[k], factor_mute)
        end
    elseif time_norm == "smooth"
        for (k, _) in R_dict
            smooth!(R_dict[k].x, time_half_win)
        end
    end
    return R_dict
end

function get_rfft_fftcc(R_dict::Dict{String,RawData}, cc_method::String)
    # N_station = R_dict.count
    Fs_dict = Dict{String,FFTData}()
    for (k, _) in R_dict
        if cc_method == "CC"
            Fs = rfft(R_dict[k])
        else
            Fs = phase(R_dict[k])
        end
        setindex!(Fs_dict, Fs, k)
    end
    return Fs_dict
end

function get_freq_norm_fftcc(
    Fs_dict::Dict{String,FFTData},
    freq_norm::String,
    freqmin::Float64,
    freqmax::Float64,
    freq_half_win::Int64,
)
    if freq_norm == "no"
        nothing
    elseif freq_norm == "whiten"
        for (k, _) in Fs_dict
            whiten!(Fs_dict[k], freqmin, freqmax)
        end
    elseif freq_norm == "coherence"
        for (k, _) in Fs_dict
            coherence!(Fs_dict[k], freq_half_win)
        end
    elseif freq_norm == "deconvolution"
        for (k, _) in Fs_dict
            deconvolution!(Fs_dict[k], freq_half_win)
        end
    end
    return Fs_dict
end

function get_cc_substack_fftcc(
    Onetime_path::String,
    Fs_dict::Dict{String,FFTData},
    All_NET_STA::Vector{Vector{String}},
    SR_pairs::Vector{Vector{String}},
    auto_corr::Bool,
    cross_corr::Bool,
    cc_method::String,
    CCFDIR::String,
    maxlag::Float64,
    substack::Bool,
    substack_method::Vector{String},
)
    year = splitpath(Onetime_path)[end-2]
    yearday = splitpath(Onetime_path)[end-1]
    daytime = splitpath(Onetime_path)[end]

    if auto_corr == true
        Threads.@threads for i = 1:size(All_NET_STA, 1)
            # debug:
            println("thread num = ", Threads.threadid())
            #
            net, sta = All_NET_STA[i][1], All_NET_STA[i][2]
            key = string(net, '.', sta)
            C = correlate(Fs_dict[key], Fs_dict[key], maxlag, corr_type = cc_method)
            if substack == true
                for jj in substack_method
                    if jj == "mean"
                        C_mean = stack(C, allstack = true)
                        SPLITOUT_mean =
                            joinpath(CCFDIR, "substack_mean", year, yearday, daytime, key)
                        if !isdir(SPLITOUT_mean)
                            mkpath(SPLITOUT_mean)
                        end
                        save_corr(C_mean, SPLITOUT_mean)
                    end
                    if jj == "pws"
                        C_pws = pws(C)
                        SPLITOUT_pws =
                            joinpath(CCFDIR, "substack_pws", year, yearday, daytime, key)
                        if !isdir(SPLITOUT_pws)
                            mkpath(SPLITOUT_pws)
                        end
                        save_corr(C_pws, SPLITOUT_pws)
                    end
                    if jj == "robust"
                        C_robust = deepcopy(C)
                        C_robust.corr = robustpws(C.corr)
                        SPLITOUT_robust =
                            joinpath(CCFDIR, "substack_robust", year, yearday, daytime, key)
                        if !isdir(SPLITOUT_robust)
                            mkpath(SPLITOUT_robust)
                        end
                        save_corr(C_robust, SPLITOUT_robust)
                    end
                end
            else
                SPLITOUT_no = joinpath(CCFDIR, "substack_no", year, yearday, daytime, key)
                if !isdir(SPLITOUT_no)
                    mkpath(SPLITOUT_no)
                end
                save_corr(C, SPLITOUT_no)
            end
        end
    end

    if cross_corr == true
        Threads.@threads for i = 1:size(SR_pairs, 1)
            s_net, s_sta, r_net, r_sta =
                SR_pairs[i][1], SR_pairs[i][2], SR_pairs[i][3], SR_pairs[i][4]
            key_1 = string(s_net, '.', s_sta)
            key_2 = string(r_net, '.', r_sta)
            C = correlate(Fs_dict[key_1], Fs_dict[key_2], maxlag, corr_type = cc_method)
            if substack == true
                for jj in substack_method
                    if jj == "mean"
                        C_mean = stack(C, allstack = true)
                        SPLITOUT_mean =
                            joinpath(CCFDIR, "substack_mean", year, yearday, daytime, key_1)
                        if !isdir(SPLITOUT_mean)
                            mkpath(SPLITOUT_mean)
                        end
                        save_corr(C_mean, SPLITOUT_mean)
                    end
                    if jj == "pws"
                        C_pws = pws(C)
                        SPLITOUT_pws =
                            joinpath(CCFDIR, "substack_pws", year, yearday, daytime, key_1)
                        if !isdir(SPLITOUT_pws)
                            mkpath(SPLITOUT_pws)
                        end
                        save_corr(C_pws, SPLITOUT_pws)
                    end
                    if jj == "robust"
                        C_robust = deepcopy(C)
                        C_robust.corr = robustpws(C.corr)
                        SPLITOUT_robust = joinpath(
                            CCFDIR,
                            "substack_robust",
                            year,
                            yearday,
                            daytime,
                            key_1,
                        )
                        if !isdir(SPLITOUT_robust)
                            mkpath(SPLITOUT_robust)
                        end
                        save_corr(C_robust, SPLITOUT_robust)
                    end
                end
            else
                SPLITOUT_no = joinpath(CCFDIR, "substack_no", year, yearday, daytime, key_1)
                if !isdir(SPLITOUT_no)
                    mkpath(SPLITOUT_no)
                end
                save_corr(C, SPLITOUT_no)
            end
        end
    end

    return nothing
end

function main_fft_cc(Onetime_path::String, Parameters::Dict{String,Any})
    flag = Parameters["flag"]
    channel_regular = Parameters["channel_regular"]
    select_sta = Parameters["select_sta"]
    STATXT = Parameters["STATXT"]
    SRTXT = Parameters["SRTXT"]
    input_fmt = Parameters["input_fmt"]
    rm_resp = Parameters["rm_resp"]
    respdir = Parameters["respdir"]
    cc_len = Parameters["cc_len"]
    cc_step = Parameters["cc_step"]
    time_norm = Parameters["time_norm"]
    factor_clip_std = Parameters["factor_clip_std"]
    factor_mute = Parameters["factor_mute"]
    time_half_win = Parameters["time_half_win"]
    cc_method = Parameters["cc_method"]
    freq_norm = Parameters["freq_norm"]
    freqmin = Parameters["freqmin"]
    freqmax = Parameters["freqmax"]
    freq_half_win = Parameters["freq_half_win"]
    auto_corr = Parameters["auto_corr"]
    cross_corr = Parameters["cross_corr"]
    CCFDIR = Parameters["CCFDIR"]
    maxlag = Parameters["maxlag"]
    substack = Parameters["substack"]
    substack_method = Parameters["substack_method"]

    println("Cross-correlating Onetime_day: $(Onetime_path) || now:", now())

    if flag == true
        ###### 1.read sac/mseed -> SeisData() ######
        tALL = @elapsed begin
            S_dict, All_NET_STA, SR_pairs = get_SeisData_fftcc(
                Onetime_path,
                channel_regular,
                select_sta,
                STATXT,
                SRTXT,
                input_fmt,
            )
        end
        println("'Func:get_SeisData' took $tALL seconds")

        ###### 2.remove response ######
        #...

        ###### 3.SeisData() -> RawData() ######
        tALL = @elapsed begin
            R_dict = get_SeisData2RawData_fftcc(S_dict, cc_len, cc_step)
        end
        println("'Func:get_SeisData2RawData' took $tALL seconds")

        ###### 4.time_norm ######
        tALL = @elapsed begin
            R_dict = get_time_norm_fftcc(
                R_dict,
                time_norm,
                factor_clip_std,
                factor_mute,
                time_half_win,
            )
        end
        println("'Func:get_time_norm' took $tALL seconds")

        ###### 5.fft ######
        tALL = @elapsed begin
            Fs_dict = get_rfft_fftcc(R_dict, cc_method)
        end
        println("'Func:get_rfft' took $tALL seconds")

        ###### 6.freq_norm ######
        tALL = @elapsed begin
            Fs_dict = get_freq_norm_fftcc(Fs_dict, freq_norm, freqmin, freqmax, freq_half_win)
        end
        println("'Func:get_freq_norm' took $tALL seconds")

        ###### 7.cc_substack ######
        tALL = @elapsed begin
            get_cc_substack_fftcc(
                Onetime_path,
                Fs_dict,
                All_NET_STA,
                SR_pairs,
                auto_corr,
                cross_corr,
                cc_method,
                CCFDIR,
                maxlag,
                substack,
                substack_method,
            )
        end
        println("'Func:get_cc_substack' took $tALL seconds")
    else
        ###### 1.read sac/mseed -> SeisData() ######
        S_dict, All_NET_STA, SR_pairs = get_SeisData_fftcc(
            Onetime_path,
            channel_regular,
            select_sta,
            STATXT,
            SRTXT,
            input_fmt,
        )

        ###### 2.remove response ######
        #...

        ###### 3.SeisData() -> RawData() ######
        R_dict = get_SeisData2RawData_fftcc(S_dict, cc_len, cc_step)

        ###### 4.time_norm ######
        R_dict =
            get_time_norm_fftcc(R_dict, time_norm, factor_clip_std, factor_mute, time_half_win)

        ###### 5.fft ######
        Fs_dict = get_rfft_fftcc(R_dict, cc_method)

        ###### 6.freq_norm ######
        Fs_dict = get_freq_norm_fftcc(Fs_dict, freq_norm, freqmin, freqmax, freq_half_win)

        ###### 7.cc_substack ######
        get_cc_substack_fftcc(
            Onetime_path,
            Fs_dict,
            All_NET_STA,
            SR_pairs,
            auto_corr,
            cross_corr,
            cc_method,
            CCFDIR,
            maxlag,
            substack,
            substack_method,
        )

    end
    return nothing
end

function FFT_CC(
    cores::Int64,
    threads::Int64,
    flag::Bool,
    precompile::Bool,
    DATADIR::String,
    CCFDIR::String,
    select_time::Bool,
    TIMETXT::String,
    select_sta::String,
    STATXT::String,
    SRTXT::String,
    channel_regular::String,
    rm_resp::String,
    respdir::String,
    input_fmt::String,
    time_norm::String,
    freq_norm::String,
    cc_method::String,
    auto_corr::Bool,
    cross_corr::Bool,
    freqmin::Float64,
    freqmax::Float64,
    cc_len::Int64,
    cc_step::Int64,
    maxlag::Float64,
    substack::Bool,
    substack_method::Vector{String},
    factor_clip_std::Float64,
    factor_mute::Float64,
    time_half_win::Int64,
    freq_half_win::Int64,
)

    ### a. parameters --> Dict-Parameters
    Parameters = Dict(
        "cores" => cores,
        "threads" => threads,
        "flag" => flag,
        "DATADIR" => DATADIR,
        "CCFDIR" => CCFDIR,
        "select_time" => select_time,
        "TIMETXT" => TIMETXT,
        "select_sta" => select_sta,
        "STATXT" => STATXT,
        "SRTXT" => SRTXT,
        "rm_resp" => rm_resp,
        "respdir" => respdir,
        "input_fmt" => input_fmt,
        "time_norm" => time_norm,
        "freq_norm" => freq_norm,
        "cc_method" => cc_method,
        "auto_corr" => auto_corr,
        "cross_corr" => cross_corr,
        "channel_regular" => channel_regular,
        "cc_len" => cc_len,
        "cc_step" => cc_step,
        "maxlag" => maxlag,
        "freqmin" => freqmin,
        "freqmax" => freqmax,
        "substack" => substack,
        "substack_method" => substack_method,
        "factor_clip_std" => factor_clip_std,
        "factor_mute" => factor_mute,
        "time_half_win" => time_half_win,
        "freq_half_win" => freq_half_win,
    )

    ### b. get_alltime_path
    alltime_path = get_alltime_path_fftcc(DATADIR, select_time, TIMETXT)
    N_time_chunks = size(alltime_path, 1)

    ### c. precompile
    if precompile == true
        println("precompile begin...")
        tALL_pre = @elapsed begin
            pmap(main_fft_cc, [alltime_path[1]], [Parameters])
            toremove = glob("*/*/*/*/*/*", Parameters["CCFDIR"])
            rm.(toremove)
        end
        println("      precompile took $tALL_pre seconds")
        println("precompile end...")
    end

    ### d. run
    println("\n\n******************* START ******************")
    println(
        "Beginning XCORR with $(cores) processors and each core with $(threads) threads... || N_time_chunks: $(N_time_chunks)",
    )
    tALL = @elapsed begin
        pmap(main_fft_cc, alltime_path, fill(Parameters, N_time_chunks))
    end
    println("XCORR took $tALL seconds")
    println("******************* END ******************\n\n")
end

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


text:
module ccall

using Logging
@Logging.configure(level=DEBUG)
# package code goes here

end # module

ext: jl


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Round 1A 2008
# Problem B: Milkshakes

# Determines whether customer c is satisfied by the milkshake malts ms.
function satisfied(c, ms)
	for (f, m) in c
		if ms[f] == m
			return true
		end
	end
	false
end

# Returns the first malted flavour in the preferences of customer c. If all
# flavours are unmalted, returns zero.
function firstmalted(c)
	i = findfirst(p -> p[2] == 1, c)
	i == 0 ? 0 : c[i][1]
end

# Returns the list of nf malts (0 or 1) such that the number of malted flavours
# is a minimum while satisfying the preferences of all customers in cs.
function minmalts(nf, cs)
	ms = zeros(Int, nf)
	done = false
	while !done
		done = true
		for c in cs
			if !satisfied(c, ms)
				f = firstmalted(c)
				if f == 0
					return nothing
				end
				ms[f] = 1
				done = false
			end
		end
	end
	ms
end

# Parse a string as a customer's milkshake preferences. The first integer is
# ignored (the number of preferences), and subsequent pairs of integers are
# parsed as (f, m) tuples where f is the flavour and m is the malt.
function parse_customer(str)
	xs = int(split(str))
	[(xs[i], xs[i+1]) for i in 2:2:length(xs)-1]
end

function main()
	n = int(readline())
	for i = 1:n
		nf = int(readline()) # number of flavours
		nc = int(readline()) # number of customers
		cs = map(ci -> parse_customer(readline()), 1:nc)
		ms = minmalts(nf, cs)
		str = (ms === nothing) ? "IMPOSSIBLE" : join(ms, " ")
		println("Case #$i: ", str)
	end
end

main()


text:
# ------------------------------------------------------------
# Extended Flux.Chain
# ------------------------------------------------------------
function Flux.Chain(D::Integer, N::Integer, σ::Function)
    t = ()
    for i = 1:N
        t = (t..., Dense(D, D, σ))
    end

    return Chain(t...)
end

"""
    dense_layer(
        in::T,
        out::T,
        σ = identity::Function;
        fw = randn::Function,
        fb = zeros::Function,
        isBias = true::Bool,
        precision = Float32,
    ) where {T<:Integer}
Create dense layer with meticulous settings
"""
function dense_layer(
    in::T,
    out::T,
    σ = identity::Function;
    fw = randn::Function,
    fb = zeros::Function,
    isBias = true::Bool,
    precision = Float32,
) where {T<:Integer}
    if isBias
        return Dense(fw(precision, out, in), fb(precision, out), σ)
    else
        return Dense(fw(precision, out, in), Flux.Zeros(precision, out), σ)
    end
end


text:
using EulerLagrange
using Test

@testset "EulerLagrange.jl" begin
    # Write your tests here.
end


text: @non_differentiable has_fermion_string(::AbstractString, ::Index)


text:
include("splinesParallel.jl")

Su0p = plarBezier(S1, [[0,0,0],[10,0,0]])
Su1p = plarBezier(S1, [[0,10,0],[2.5,10,3],[5,10,-3],[7.5,10,3],[10,10,0]])
Sv0p = plarBezier(S2, [[0,0,0],[0,0,3],[0,10,3],[0,10,0]])
Sv1p = plarBezier(S2, [[10,0,0],[10,5,3],[10,10,0]])

coons = plarCoonsPatch([Su0p,Su1p,Sv0p,Sv1p])
point = [0.2, 0.7]
coons(point)


text:
@enum Direction Left=1 Right=-1
@enum Bearing East=0 North=1 West=2 South=3

const PatchDict = Dict{Tuple{Int,Int},PyCall.PyObject}

const COLORS = Dict{NTuple{3,Float64}, String}(
    (1.0, 0.0, 0.0) => "red",
    (0.0, 1.0, 0.0) => "green",
    (0.0, 0.0, 1.0) => "blue",
    (0.5, 0.5, 0.5) => "gray"
)

const SPEED = 2.0


text:
"""
    Coulomb(; cutoff, nl_only, weight_14, coulomb_const, force_unit, energy_unit)

The Coulomb electrostatic interaction.
"""
struct Coulomb{C, W, T, F, E} <: GeneralInteraction
    cutoff::C
    nl_only::Bool
    weight_14::W
    coulomb_const::T
    force_unit::F
    energy_unit::E
end

const coulombconst = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0

function Coulomb(;
                    cutoff=NoCutoff(),
                    nl_only=false,
                    weight_14=1.0,
                    coulomb_const=coulombconst,
                    force_unit=u"kJ * mol^-1 * nm^-1",
                    energy_unit=u"kJ * mol^-1")
    return Coulomb{typeof(cutoff), typeof(weight_14), typeof(coulomb_const), typeof(force_unit), typeof(energy_unit)}(
        cutoff, nl_only, weight_14, coulomb_const, force_unit, energy_unit)
end

@inline @inbounds function force(inter::Coulomb{C},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size,
                                    weight_14::Bool=false) where C
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    cutoff = inter.cutoff
    coulomb_const = inter.coulomb_const
    qi, qj = atom_i.charge, atom_j.charge

    params = (coulomb_const, qi, qj)

    if cutoff_points(C) == 0
        f = force_divr_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_unit

        f = force_divr_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_unit

        if r2 < cutoff.activation_dist
            f = force_divr_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_divr_cutoff(cutoff, r2, inter, params)
        end
    end

    if weight_14
        return f * dr * inter.weight_14
    else
        return f * dr
    end
end

@fastmath function force_divr_nocutoff(::Coulomb, r2, invr2, (coulomb_const, qi, qj))
    (coulomb_const * qi * qj) / √(r2 ^ 3)
end

@inline @inbounds function potential_energy(inter::Coulomb{C},
                                    s::System,
                                    i::Integer,
                                    j::Integer,
                                    weight_14::Bool=false) where C
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)

    cutoff = inter.cutoff
    coulomb_const = inter.coulomb_const
    qi, qj = s.atoms[i].charge, s.atoms[j].charge
    params = (coulomb_const, qi, qj)

    if cutoff_points(C) == 0
        pe = potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(s.box_size[1])) * inter.energy_unit

        pe = potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(s.box_size[1])) * inter.energy_unit

        if r2 < cutoff.activation_dist
            pe = potential(inter, r2, inv(r2), params)
        else
            pe = potential_cutoff(cutoff, r2, inter, params)
        end
    end

    if weight_14
        return pe * inter.weight_14
    else
        return pe
    end
end

@fastmath function potential(::Coulomb, r2, invr2, (coulomb_const, qi, qj))
    (coulomb_const * qi * qj) * √invr2
end


text:
using MPI

function do_sendrecv()
    comm  = MPI.COMM_WORLD
    MPI.Barrier(comm)

    rk = MPI.Comm_rank(comm)
    sz = MPI.Comm_size(comm)

    dst = mod(rk+1, sz)
    src = mod(rk-1, sz)

    N = 4

    send_msg = Array{Float64}(undef,N)
    fill!(send_msg, Float64(rk))
    recv_msg = Array{Float64}(undef,N)

    tag = 1
    r_req = MPI.Irecv!(recv_msg, src, tag, comm)
    println("$rk: Sending, $rk -> $dst = $send_msg")

    s_req  = MPI.Isend(send_msg,dst,tag, comm)

    stats = MPI.Waitall!([r_req,s_req])
    println("$rk: Receiving $src -> $rk = $recv_msg")

    MPI.Barrier(comm)

end

function main()
    MPI.Init()
    do_sendrecv()
    MPI.Finalize()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Random
#Implementing the general Advantage Estimation according to https://arxiv.org/abs/1506.02438

#Return the Generalized Advantage Estimation:
@everywhere function gae!(params::Params,traj::Trajectory)
    l_traj = traj.len
    γ      = params.γ_disc
    λ      = params.λ_actor

    Val = traj.val_t
    r_t = traj.r_t

    #Calculate the local TD-error:
    delta_V_i = zeros(l_traj)
    for i_t = 1:l_traj-1
        delta_V_i[i_t] = -Val[i_t] + r_t[i_t] + γ*Val[i_t+1]
    end

    #Discount the TD-Error as the Advantage
    A_t = vec(discount(delta_V_i,λ*γ))

    #Write result in trajectory
    traj.A_t = A_t
end

#Parallel PPO update
#-------------------------------------------------------------------------------
"""
Function training the Value function on the last batch of data in parallel
"""
function refit_VF(params::Params,θv::Any,all_traj::Traj_Data, opt_para_v)
        #Take all partial traj and shuffle their order!
        n_traj = all_traj.n_traj
        shuffle_arr = Random.shuffle(1:n_traj)

        active_workers = workers()
        n_workers = length(active_workers)
        ĝ_v = copy(θv)

        # Create Channels storing results and starting values
        input_channel  = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks and traj
        output_channel = RemoteChannel(()->Channel{Array{Any,1}}(length(workers())+1)) # will contain networks

        # Create problem on every worker
        begin
            for p in workers()
                remote_do(distributed_vf_training, p, params, input_channel, output_channel)
            end
        end

        # Start optimizations
        n_parallel = Int(ceil(n_traj / n_workers))
        prog = Progress(n_traj,1)

        for i_parallel = 1 : n_parallel

            #1. Send the work tasks:
            for i_w = 1 : n_workers
                i_true = ((i_parallel-1)*n_workers + i_w) #get true order
                if i_true <= n_traj
                    i_traj = shuffle_arr[i_true] #shuffle order
                    traj_local = all_traj.traj[i_traj]
                    @async put!(input_channel, (traj_local, θv))
                end
            end

            #2. Take the result
            ĝ_v *= 0
            for i_w = 1: n_workers
                i_true = ((i_parallel-1)*n_workers + i_w)
                if i_true <= n_traj
                    #Take from Channel
                    ĝ_v += take!(output_channel)
                end
            end

            #Update the Network
            Knet.update!(θv, ĝ_v, opt_para_v)
        end

        # Tell workers that its done
        NaN_mat = zeros(length(params.batch_size))
        fill!(NaN_mat,NaN)
        ending_traj = init_traj(params,params.batch_size)
        ending_traj.r_t = NaN_mat

        for j = 1 : n_workers
                @async put!(input_channel, (ending_traj, θv))
        end

        #Close Channels
        sleep(0.1)
        close(input_channel)
        close(output_channel)

        #Return the updated Network
        return θv
end

# Function running on the worker
@everywhere function distributed_vf_training(params::Params, input_channel::RemoteChannel, output_channel::RemoteChannel)

    # Start loop where work is done
    while true
        # break loop if NaN state received
        try
            #Take the new Input
            input = take!(input_channel)

            #Reorder data:
            traj_local = input[1]
            θv         = input[2]

            #Check if u reached final
            if (sum(isnan, traj_local.r_t) > 0)
                      break
            end

            #DO THE WORK:
            #-----------------------------------------------------------------------
            #Calculate the Value Function gradient using TD(1)
            ĝ_v_loc = grad_value(θv,params,traj_local)
            #-----------------------------------------------------------------------

            put!(output_channel, ĝ_v_loc)

        catch
            println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
            println(traj_local)
            break
        end
    end

    return nothing
end

#Function applying general advantage estimation for each trajectory
function gae_all!(all_traj::Traj_Data,θv::Any,params::Params)
        #Take all partial traj and shuffle their order!
        n_traj = all_traj.n_traj

        active_workers = workers()
        n_workers = length(active_workers)

        # Create Channels storing results and starting values
        input_channel  = RemoteChannel(()->Channel{Trajectory}(length(workers())+1)) # will contain networks and traj
        output_channel = RemoteChannel(()->Channel{Trajectory}(length(workers())+1)) # will contain networks

        # Create problem on every worker
        begin
            for p in workers()
                remote_do(distributed_gae, p, params, input_channel, output_channel, θv)
            end
        end

        # Start optimizations
        n_parallel = Int(ceil(n_traj / n_workers))
        prog       = Progress(n_traj,1)

        # Start optimizations
        for i_traj = 1:n_traj
                traj_local = all_traj.traj[i_traj]
                @async put!(input_channel, traj_local)
        end

        for i_traj = 1:n_traj
                #Take from Channel:
                all_traj.traj[i_traj] = copy_traj(take!(output_channel),params)
        end

        # Tell workers that its done
        NaN_mat = zeros(length(params.batch_size))
        fill!(NaN_mat,NaN)
        ending_traj = init_traj(params,params.batch_size)
        ending_traj.r_t = NaN_mat

        for j = 1 : n_workers
                @async put!(input_channel, ending_traj)
        end

        #Close Channels
        sleep(0.1)
        close(input_channel)
        close(output_channel)
end

#Function running on the worker
@everywhere function distributed_gae(params::Params, input_channel::RemoteChannel, output_channel::RemoteChannel, θv)

   traj_local = init_traj(params,params.l_traj)
   # Start loop where work is done
   while true
       # break loop if NaN state received
       try
           #Take the new Input
           traj_local = copy_traj(take!(input_channel),params)

           #Check if u reached final
           if (sum(isnan, traj_local.r_t) > 0)
                     break
           end

           #DO THE WORK:
           #-----------------------------------------------------------------------
           #Recalculate the values for the states
           generate_values!(traj_local,params,θv)
           #recalculate the gae
           gae!(params,traj_local)
           #-----------------------------------------------------------------------

           put!(output_channel,copy_traj(traj_local,params))

       catch
           println(traj_local)
           println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
           break
       end
   end

   return nothing
end

#Function calculating the values for each state in a trajectory
@everywhere function generate_values!(traj::Trajectory,params::Params,θv::Any)
    for i_t = 1:traj.len
        s_t = traj.s_t[i_t,:]

        #Only for Segway
        s_mod = s_t
        #=---------------------------------
        s_mod        = zeros(8)
        s_mod[1:2]   = s_t[1:2]
        s_mod[3]     = sin(s_t[3])
        s_mod[4]     = cos(s_t[3])
        s_mod[5:end] = s_t[4:end]
        #-------------------------------=#

        val_t = mlp_nlayer_value_norm(θv,s_mod,params)[1]
        traj.val_t[i_t] = val_t
    end
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
mutable struct LUState{S<:State, L} <: FiniteBasisState
    state::S
    lookup::L
    valid::Bool
end

LUState(s, lu) = LUState(s, lu, false)
export LUState

# accessors
state(s::LUState) = s.state
lut(s::LUState)   = s.lookup
isvalid(s::LUState) = s.valid

raw_state(s::LUState) = raw_state(state(s))
invalidate!(s::LUState) = s.valid = false
validate!(s::LUState) = s.valid = true
# Modifiers
init_lut!(s::LUState, net::NeuralNetwork, force::Bool=false) = begin
    if isvalid(s) && !force
        return update_lut!(s, net)
    end
    apply_chages!(state(s))
    prepare_lut!(s, net)
    validate!(s)
    return s
end

init_lut!(s::State, net) = nothing
update_lut!(s::State, net) = nothing

update_lut!(s::LUState, net::NeuralNetwork) = begin
    apply_lut_updates!(s, net)
    apply_chages!(state(s))
    return s
end

clear_changes!(s::LUState) = clear_changes!(state(s))

## State Interface : custom accessors
@inline config(s::LUState) = (invalidate!(s); config(state(s)) )

## State Interface : Property Accesors
@inline spacedimension(s::LUState) = spacedimension(state(s))
@inline nsites(s::LUState) = nsites(state(s))
@inline local_dimension(s::LUState) = local_dimension(state(s))
@inline eltype(s::LUState) = eltype(state(s))

toint(s::LUState) = (invalidate!(s); toint(state(s)))
index(s::LUState) = (invalidate!(s); return index(state(s)))
index_to_int(s::LUState) = (invalidate!(s); index_to_int(state(s)))

## State Interface : Checks
same_basis(s1::LUState, s2::LUState) =
    same_basis(state(s1), state(s2))

## State Interface : flip operations
# This is the standard method, that applies all exhisting transformations
function flipat!(rng::AbstractRNG, s::LUState, args...)
    invalidate!(s)
    return flipat!(rng, state(s), args...)
end

function flipat_fast!(rng::AbstractRNG, s::LUState, args...)
    return flipat_fast!(rng, state(s), args...)
end

setat!(s::LUState, args...) = (invalidate!(s);
    setat!(state(s), args...); s)

set_index!(s::LUState, val) = (invalidate!(s);
    set_index!(state(s), val); s)

set!(s::LUState, val) = begin
    invalidate!(s)
    set!(state(s), val)
    s
end
add!(s::LUState, val) = (invalidate!(s);
    add!(state(state), val); s)
rand!(rng, s::LUState) = (invalidate!(s);
    rand!(rng, state(s)); s)

Base.show(io::IO, ::MIME"text/plain", s::LUState) =
    print(io, "LUState: $(state(s)) valid: $(s.valid)")


text:
#=
day5:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-05
=#

using AoC, Test
using Base.Iterators: filter

parseLine(line) = split(line, r",| -> ") |> cc -> parse.(Int, cc) |> cc -> reshape(cc, 2,2)
@test parseLine("3,5 -> 6,8") == [[3, 5] [6,8]]

parseLines(lines) = cat((parseLine.(lines))...; dims=3)
@test parseLines(["3,5 -> 6,8", "4,6 -> 2,3"]) == [[3,5] [6,8] ;;; [4,6] [2,3]]

maxX(lines) = max(lines[1,:,:]...)
maxY(lines) = max(lines[2,:,:]...)
maxCoord(lines) = [maxX(lines), maxY(lines)]
@test maxCoord([[1,2] [3,4] ;;; [5,3] [1,2]]) == [5,4]

minX(lines) = min(lines[1,:,:]...)
minY(lines) = min(lines[2,:,:]...)
minCoord(lines) = [minX(lines), minY(lines)]
@test minCoord([[0,2] [3,4] ;;; [5,3] [1,2]]) == [0,2]

mapOffset(lines) = [1,1] - minCoord(lines)

mapSize(lines) = maxCoord(lines) + mapOffset(lines)

direction(line) = direction(line[:,1],line[:,2])
direction(start, finish) = sign.(finish - start)
@test direction([[1,1] [7,7]]) == [1,1]
@test direction([[1,7] [7,7]]) == [1,0]
@test direction([[7,1] [1,7]]) == [-1,1]

function createMap(lines)
    map = zeros(Int, mapSize(lines)...)
    offset = mapOffset(lines)

    for line in eachslice(lines; dims=3)
        point = line[:,1] + offset
        endPoint = line[:,2] + offset
        dir = direction(line)
        while true
            map[point...]+=1
            if point == endPoint
                break
            end
            point += dir
        end
    end

    return map
end

selectPart1(lines) = cat(filter(eachslice(lines, dims=3)) do line
    (start, finish) = (line[:,1], line[:,2])
    return start[1] == finish[1] || start[2] == finish[2]
end...; dims=3)
@test selectPart1([[2,2] [4,4] ;;; [2,2] [2,5]]) == [[2,2] [2,5] ;;;]

part1(lines) = parseLines(lines) |> selectPart1 |> createMap |> map -> count(x -> x>=2, map)
part2(lines) = parseLines(lines) |> createMap |> map -> count(x -> x>=2, map)

@test part1(exampleLines(5,1)) == 5
@test part2(exampleLines(5,1)) == 12

lines(5) |> ll -> @time part1(ll) |> show
lines(5) |> ll -> @time part2(ll) |> show


text:
thermal_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch)
load_model = DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad)
line_model = DeviceModel(PSY.Line, PSI.StaticLine)
transformer_model = DeviceModel(PSY.Transformer2W, PSI.StaticTransformer)
ttransformer_model = DeviceModel(PSY.TapTransformer, PSI.StaticTransformer)
dc_line = DeviceModel(PSY.HVDCLine, PSI.HVDCDispatch)

@testset "Network Copper Plate" begin
    network = CopperPlatePowerModel
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters = [true, false]
    test_results = Dict{PSY.System, Vector{Int64}}(c_sys5 => [120, 120, 0, 0, 24],
                                                   c_sys14 => [120, 120, 0, 0, 24],
                                                   c_sys14_dc => [120, 120, 0, 0, 24])
    constraint_names = [:CopperPlateBalance]
    objfuncs = [GAEVF, GQEVF, GQEVF]

    for (ix, sys) in enumerate(systems), p in parameters
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = OSQP_optimizer, parameters = p)
        construct_device!(ps_model, :Thermal, thermal_model; parameters = p);
        construct_device!(ps_model, :Load, load_model; parameters = p);
        construct_network!(ps_model, network; parameters = p);
        construct_device!(ps_model, :Line, line_model; parameters = p);
        construct_device!(ps_model, :Tf, transformer_model; parameters = p);
        construct_device!(ps_model, :TTf, ttransformer_model; parameters = p);
        construct_device!(ps_model, :DCLine, dc_line; parameters = p);

        moi_tests(ps_model, p, test_results[sys][1],
                                test_results[sys][2],
                                test_results[sys][3],
                                test_results[sys][4],
                                test_results[sys][5], false)
        psi_constraint_test(ps_model, constraint_names)
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(ps_model, [MOI.OPTIMAL])
    end
end

@testset "Network DC-PF with PTDF formulation" begin
    network = StandardPTDFModel
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names = [:RateLimit_Line, :nodal_balance, :network_flow]
    parameters = [true, false]
    PTDF_ref = Dict{PSY.System, PSY.PTDF}(c_sys5 => PTDF5, c_sys14 => PTDF14, c_sys14_dc => PTDF14_dc);
    test_results = Dict{PSY.System, Vector{Int64}}(c_sys5 => [264, 264, 0, 0, 264],
                                                    c_sys14 => [600, 600, 0, 0, 816],
                                                    c_sys14_dc => [600, 600, 0, 0, 768])

    for (ix, sys) in enumerate(systems), p in parameters
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = OSQP_optimizer, parameters = p)
        construct_device!(ps_model, :Thermal, thermal_model; parameters = p);
        construct_device!(ps_model, :Load, load_model; parameters = p);
        construct_network!(ps_model, network; PTDF = PTDF_ref[sys], parameters = p);
        construct_device!(ps_model, :Line, line_model; parameters = p);
        construct_device!(ps_model, :Tf, transformer_model; parameters = p);
        construct_device!(ps_model, :TTf, ttransformer_model; parameters = p);
        construct_device!(ps_model, :DCLine, dc_line; parameters = p);

        moi_tests(ps_model, p, test_results[sys][1],
                                test_results[sys][2],
                                test_results[sys][3],
                                test_results[sys][4],
                                test_results[sys][5], false)
        psi_constraint_test(ps_model, constraint_names)
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(ps_model, [MOI.OPTIMAL])
    end

    #PTDF input Error testing
    ps_model = OperationModel(TestOptModel, network, c_sys5; optimizer = GLPK_optimizer)
    construct_device!(ps_model, :Thermal, thermal_model);
    construct_device!(ps_model, :Load, load_model);
    @test_throws ArgumentError construct_network!(ps_model, network)

end

@testset "Network DC lossless -PF network with PowerModels DCPlosslessForm" begin
    network = DCPPowerModel
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names = [:RateLimit_Line]
    parameters = [true, false]
    test_results = Dict{PSY.System, Vector{Int64}}(c_sys5 => [384, 264, 144, 144, 288],
                                                    c_sys14 => [936, 600, 480, 480, 840],
                                                    c_sys14_dc => [984, 600, 432, 432, 840])

    for (ix, sys) in enumerate(systems), p in parameters
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = OSQP_optimizer, parameters = p)
        construct_device!(ps_model, :Thermal, thermal_model; parameters = p);
        construct_device!(ps_model, :Load, load_model; parameters = p);
        construct_network!(ps_model, network; parameters = p);
        construct_device!(ps_model, :Line, line_model; parameters = p);
        construct_device!(ps_model, :Tf, transformer_model; parameters = p);
        construct_device!(ps_model, :TTf, ttransformer_model; parameters = p);
        construct_device!(ps_model, :DCLine, dc_line; parameters = p);

        moi_tests(ps_model, p, test_results[sys][1],
                                test_results[sys][2],
                                test_results[sys][3],
                                test_results[sys][4],
                                test_results[sys][5], false)
        psi_constraint_test(ps_model, constraint_names)
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(ps_model, [MOI.OPTIMAL])
    end

end

@testset  "Network Solve AC-PF PowerModels StandardACPModel" begin
    network = ACPPowerModel
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names = [:RateLimitFT_Line, :RateLimitTF_Line]
    parameters = [true, false]
    test_results = Dict{PSY.System, Vector{Int64}}(c_sys5 => [1056, 240, 144, 144, 264],
                                                    c_sys14 => [2832, 240, 480, 480, 696],
                                                    c_sys14_dc => [2832, 336, 432, 432, 744]) # TODO: changed the interval constraint number to 336 from 240. double check

    for (ix, sys) in enumerate(systems), p in parameters
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = ipopt_optimizer, parameters = p)
        construct_device!(ps_model, :Thermal, thermal_model; parameters = p);
        construct_device!(ps_model, :Load, load_model; parameters = p);
        construct_network!(ps_model, network; parameters = p);
        construct_device!(ps_model, :Line, line_model; parameters = p);
        construct_device!(ps_model, :Tf, transformer_model; parameters = p);
        construct_device!(ps_model, :TTf, ttransformer_model; parameters = p);
        construct_device!(ps_model, :DCLine, dc_line; parameters = p);

        moi_tests(ps_model, p, test_results[sys][1],
                                test_results[sys][2],
                                test_results[sys][3],
                                test_results[sys][4],
                                test_results[sys][5], false)
        psi_constraint_test(ps_model, constraint_names)
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(ps_model, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
    end
end

@testset  "Network Solve AC-PF PowerModels linear approximation models" begin
    networks = [DCPPowerModel, NFAPowerModel]
    systems = [c_sys5, c_sys14, c_sys14_dc]
    p = true
    for network in networks, sys in systems
        @info "Testing construction of a $(network) network"
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = OSQP_optimizer, parameters = p)
        construct_device!(ps_model, :Thermal, thermal_model; parameters = p);
        construct_device!(ps_model, :Load, load_model; parameters = p);
        construct_network!(ps_model, network; parameters = p);
        construct_device!(ps_model, :Line, line_model; parameters = p);
        construct_device!(ps_model, :DCLine, dc_line; parameters = p);
        psi_checksolve_test(ps_model, [MOI.OPTIMAL])

    end

end

@testset  "Network AC-PF PowerModels non-convex models" begin
    networks = [#ACPPowerModel, Already tested
                ACRPowerModel,
                ACTPowerModel
                ]
    systems = [c_sys5, c_sys14, c_sys14_dc]

    for network in networks, sys in systems
        @info "Testing construction of a $(network) network"
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = ipopt_optimizer)
        construct_device!(ps_model, :Thermal, thermal_model);
        construct_device!(ps_model, :Load, load_model);
        construct_network!(ps_model, network);
        construct_device!(ps_model, :Line, line_model);
        @test !isnothing(ps_model.canonical.pm_model)
    end

end

@testset  "Network AC-PF PowerModels quadratic loss approximations models" begin
    networks = [DCPLLPowerModel, LPACCPowerModel]
    systems = [c_sys5, c_sys14, c_sys14_dc]

    for network in networks, sys in systems
        @info "Testing construction of a $(network) network"
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = ipopt_optimizer)
        construct_device!(ps_model, :Thermal, thermal_model);
        construct_device!(ps_model, :Load, load_model);
        construct_network!(ps_model, network);
        construct_device!(ps_model, :Line, line_model);
        construct_device!(ps_model, :DCLine, dc_line);
        @test !isnothing(ps_model.canonical.pm_model)
    end

end

@testset  "Network AC-PF PowerModels quadratic relaxations models" begin
    networks = [ SOCWRPowerModel,
                 QCRMPowerModel,
                 QCLSPowerModel,
                 ]
    systems = [c_sys5, c_sys14, c_sys14_dc]

    for network in networks, sys in systems
        @info "Testing construction of a $(network) network"
        ps_model = OperationModel(TestOptModel, network, sys; optimizer = ipopt_optimizer)
        construct_device!(ps_model, :Thermal, thermal_model);
        construct_device!(ps_model, :Load, load_model);
        construct_network!(ps_model, network);
        construct_device!(ps_model, :Line, line_model);
        construct_device!(ps_model, :DCLine, dc_line);
        @test !isnothing(ps_model.canonical.pm_model)
    end

end

@testset  "Network Unsupported Power Model Formulations" begin
incompat_list = [PM.SDPWRMPowerModel,
                 PM.SparseSDPWRMPowerModel,
                 PM.SOCBFPowerModel,
                 PM.SOCBFConicPowerModel]

    for network in incompat_list
        ps_model = OperationModel(TestOptModel, network, c_sys5; optimizer = ipopt_optimizer)
        construct_device!(ps_model, :Thermal, thermal_model);
        construct_device!(ps_model, :Load, load_model);

        @test_throws ArgumentError construct_network!(ps_model, network);
    end

end


text:
module ChainRulesCore
using Base.Broadcast: materialize, materialize!, broadcasted, Broadcasted, broadcastable

export frule, rrule
export wirtinger_conjugate, wirtinger_primal, refine_differential
export @scalar_rule, @thunk
export extern, store!
export Wirtinger, Zero, One, DoesNotExist, Thunk, InplaceableThunk
export NO_FIELDS

include("differentials.jl")
include("differential_arithmetic.jl")
include("operations.jl")
include("rules.jl")
include("rule_definition_tools.jl")

Base.@deprecate_binding DNE DoesNotExist

end # module


text:
function run_condition!(ex, stimuli; parms...)
    for trial in 1:ex.n_trials
        run_trial!(ex, stimuli[trial]...; parms...)
    end
    return nothing
end

function run_condition!(ex; parms...)
    for trial in 1:ex.n_trials
        run_trial!(ex; parms...)
    end
    return nothing
end

function run_trial!(ex; parms...)
    stimuli = generate_stimuli(ex)
    return run_trial!(ex, stimuli...; parms...)
end

function run_trial!(ex, target, present, visual_objects; parms...)
    ex.trial_fixations = Fixation[]
    actr = initialize_model(ex, target, visual_objects; parms...)
    compute_angular_size!(actr)
    orient!(actr, ex)
    ex.visible ? draw_cross!(actr, ex) : nothing
    ex.trace ? println("\n", get_time(actr), " start search sequence") : nothing
    search!(actr, ex)
    return actr
end

function generate_stimuli(ex)
    target,present = initialize_trial!(ex)
    visual_objects = ex.populate_visicon(ex, target..., present)
    return target, present, visual_objects
end

function initialize_model(ex, target, visual_objects; parms...)
    T = eltype(visual_objects)
    visicon=visual_objects
    visual_location = VisualLocation(buffer=T[])
    visual_location.iconic_memory = visual_objects
    target_chunk = Chunk(;target...)
    goal = Goal(buffer=target_chunk)
    visual = Visual(buffer=T[])
    return ACTRV(;goal, visual_location, visual, visicon, parms...)
end

function initialize_trial!(ex::Experiment)
    target = sample_target(ex)
    present = rand() < ex.base_rate ? :present : :absent
    ex.trial_fixations = Fixation[]
    data = Data()
    data.target_present = present
    data.target_shape = target.shape
    data.target_color = target.color
    ex.trial_data = data
    return target,present
end

function set_trial_type!(data)
    if data.target_present == :present
        if data.response == :present
            return data.trial_type = :hit
        else
            return data.trial_type = :miss
        end
    end
    if data.target_present == :absent
        if data.response == :present
            return data.trial_type = :fa
        else
            return data.trial_type = :cr
        end
    end
end

get_width(ex) = ex.object_width

function conjunctive_ratio(ex, target_color, target_shape, present)
    distractor_color = setdiff(ex.colors, [target_color])[1]
    distractor_shape = setdiff(ex.shapes, [target_shape])[1]
    color_fun() = populate_features((:color,:shape), [distractor_color,target_shape])
    shape_fun() = populate_features((:color,:shape), [target_color,distractor_shape])
    visicon = [VisualObject(features=color_fun(), width=get_width(ex)) for _ in 1:ex.n_color_distractors]
    temp = [VisualObject(features=shape_fun(), width=get_width(ex)) for _ in 1:ex.n_shape_distractors]
    push!(visicon, temp...)
    if present == :present
        push!(visicon, VisualObject(features=populate_features((:color,:shape),
        [target_color,target_shape]), width=get_width(ex), target=true))
    end
    set_locations!(ex, visicon)
    return visicon
end

function conjunctive_set(ex, target_color, target_shape, present)
    distractor_color = setdiff(ex.colors, [target_color])[1]
    distractor_shape = setdiff(ex.shapes, [target_shape])[1]
    n = round(Int, ex.set_size/2)
    n = max(n, 1)
    color_fun() = populate_features((:color,:shape), [distractor_color,target_shape])
    shape_fun() = populate_features((:color,:shape), [target_color,distractor_shape])
    visicon = [VisualObject(features=color_fun(), width=get_width(ex)) for _ in 1:n]
    temp = [VisualObject(features=shape_fun(), width=get_width(ex)) for _ in 1:n]
    push!(visicon, temp...)
    n == 1 ? (visicon = [rand(visicon)]) : nothing
    if present == :present
        vo = rand(visicon)
        vo.target = true
        vo.features = populate_features((:color,:shape),
        [target_color,target_shape])
    end
    set_locations!(ex, visicon)
    return visicon
end

function feature_set(ex, target_color, target_shape, present)
    distractor_color = setdiff(ex.colors, [target_color])[1]
    color_fun() = populate_features((:color,:shape), [distractor_color,target_shape])
    visicon = [VisualObject(features=color_fun(), width=get_width(ex)) for _ in 1:ex.set_size]
    if present == :present
        vo = rand(visicon)
        vo.target = true
        vo.features = populate_features((:color,:shape),
        [target_color,target_shape])
    end
    set_locations!(ex, visicon)
    return visicon
end

function set_locations!(ex::Experiment, visicon)
    n = length(visicon)
    indices = [[r,c] for r in 1:ex.n_cells for c in 1:ex.n_cells]
    locations = sample(indices, n; replace=false)
    for (l,v) in zip(locations, visicon)
        set_location!(ex, v, l)
    end
end

"""
set location of visual object at center of object
"""
function set_location!(ex::Experiment, visual_object, indices)
    center = visual_object.width/2
    lb = @. ex.cell_width*(indices-1) + center
    ub = @. ex.cell_width*indices - center
    location = @. rand(Uniform(lb, ub))
    visual_object.location = location
    return nothing
end

function sample_target(ex)
    target_color = rand(ex.colors)
    target_shape = rand(ex.shapes)
    return (color=target_color,shape=target_shape)
end


text:
## Claudia: commented these tests out because they were not part of runtests.jl
## when I added the test for readNexusTrees. So, I was not sure if they were
## not included for a reason (slow test?)

## @testset "test: read CF data" begin
## # d=readInputData("1.ms"); #tableCF0.txt
## d = (@test_nowarn PhyloNetworks.readInputData("1.ms",:all,0,["1","2","3","4","5","6"],false,"none",false,false);) #tableCF0.txt
## d = PhyloNetworks.readInputData("1.ms",:all,0,false); # writes summaryTreesQuartets.txt
## d=readInputData("1.ms",:rand,10); #tableCF3.txt
## d=readInputData("1.ms",[1,2,3,4,5]); #tableCF4.txt

## d=readInputData("1.ms","allQuartets.txt"); #tableCF1.txt
## d=readInputData("1.ms","allQuartets.txt",:rand,10); #tableCF2.txt
## d=readInputData("1.ms","allQuartets.txt",true,"try4.txt"); #try4.txt

## descData(d)

## d = readTrees2CF("1.ms");
## d = readTrees2CF("1.ms",filename="try.txt");
## d = readTrees2CF("1.ms","allQuartets.txt",filename="hola.txt");
## d = readTrees2CF("1.ms","allQuartets.txt",whichQ=:rand,numQ=10,filename="hola.txt");

## d= readTableCF("tableCF.txt");
## d2= readTableCF("tableCFbad.txt")
## end

nexusfile = joinpath(@__DIR__, "..", "examples", "test.nex")
@testset "test: reading nexus file" begin
    vnet = readNexusTrees(nexusfile)
    @test length(vnet) == 10
end


text:
module SymSpellChecker

using StringDistances: matching_blocks

export update!, SymSpell, lookup, term, set_options!

include("utils.jl")
include("symspell.jl")
include("distance.jl")
include("lookup.jl")

end # module


text:
#=display.jl - will replay the game
from start to finish, displaying the current game state.
Accepts 1 command line argument,<filename> => database
=#
#include("square.jl")
include("dParse.jl")
include("board.jl")
using BM

#include("start.jl")
using SQLite
pathToDatabase = ARGS[1]
db = SQLite.DB(pathToDatabase) #Opens the database gamefile
#global board = ST.loadBoard()
#calculates all black piece positions
res = SQLite.query(db,"""SELECT "value" FROM meta WHERE "key" = 'type';""")
gameType = get(res[1,1])
#=----Gets last move number----=#
res = SQLite.query(db,"""SELECT MAX("move_number") FROM moves;""") #Finds the last played move (maximum move_number)
lastMoveIDNullable = res[1,1] #SQL query with max move_number (POSSIBLY "NULL" if no moves have been made)

if (!isnull(lastMoveIDNullable)) #Checks that lastMoveID was not NULL
lastMoveID = get(res[1,1]) #Parses the last move move_number as an Int

else #lastMoveID is NULL
lastMoveID = 0 #move_number "0" is unsused. This implies no moves have been made

end
currentMoveID=1

#=Helper functions=#
function displaymini()
  board = BM.startGame(gameType)
  global dboard = Array(Char,9,9)
  for (x,y) in keys(board)
    dboard[(x,y)]=board[(x,y)][1]
  end

  for x_index in (1:11)
    for y_index in (1:11)
      if y_index==1
        if x_index==1
          print("┌")
        elseif x_index==11
          print("└")
        elseif rem(x_index,2)==0
          print("|")
        else
          print("├")
        end
      elseif y_index==11
        if x_index==1
          print("┐")
        elseif x_index==11
          print("┘")
        elseif rem(x_index,2)==0
          print("|")
        else
          print("┤")
        end
      elseif rem(y_index,2)==1
        if x_index==1
          print("┬")
        elseif x_index==11
          print("┴")
        elseif rem(x_index,2)==1
          print("┼")
        else
          print("|")
        end
      else
        if rem(x_index,2)==1
          print("-")
        else
          if dboard[(div(x_index,2),(6-div(y_index,2)))]=='k' || dboard[(div(x_index,2),(6-div(y_index,2)))]=='K'
            print_with_color(:yellow, "$(dboard[div(x_index,2),(6-div(y_index,2))])")
            continue
          end
          print(dboard[(div(x_index,2),(6-div(y_index,2)))])
        end
      end
    end
    print("\n")
  end
end

function displaystandard()

  board = BM.startGame(gameType)
  global dboard = Array(Char,9,9)
  for (x,y) in keys(board)
    dboard[(x,y)]=board[(x,y)][1]
  end

for x_index in (1:19)
  for y_index in (1:19)
    if y_index==1
      if x_index==1
        print("┌")
      elseif x_index==19
        print("└")
      elseif rem(x_index,2)==0
        print("|")
      else
        print("├")
      end
    elseif y_index==19
      if x_index==1
        print("┐")
      elseif x_index==19
        print("┘")
      elseif rem(x_index,2)==0
        print("|")
      else
        print("┤")
      end
    elseif rem(y_index,2)==1
      if x_index==1
        print("┬")
      elseif x_index==19
        print("┴")
      elseif rem(x_index,2)==1
        print("┼")
      else
        print("|")
      end
    else
      if rem(x_index,2)==1
        print("-")
      else
        if dboard[(div(x_index,2),(10-div(y_index,2)))]=='k' || dboard[(div(x_index,2),(10-div(y_index,2)))]=='K'
          print_with_color(:yellow, "$(dboard[(div(x_index,2),(10-div(y_index,2)))])")
          continue
        end
        print(dboard[(div(x_index,2),(10-div(y_index,2)))])
      end
    end
  end
  print("\n")
end

function Setupboard()
  res = SQLite.query(db,"""SELECT sourcex,sourcey,targetx,targety,move_type,option FROM moves WHERE "move_number" = $(currentMoveID);""")
  sourcexNullable = res[1,1]
  sourceyNullable = res[1,2]
  targetxNullable = res[1,3]
  targetyNullable = res[1,4]
  targetx = get(targetxNullable)
  targety = get(targetyNullable)
  sourcex = get(sourcexNullable)
  sourcey = get(sourceyNullable)
  board[targetx,targety].piece = board[sourcex,sourcey].piece
  board[targetx,targety].team = board[sourcex,sourcey].team
  ST.clear!(board[sourcex,sourcey])
  ST.saveBoard(board)
end

end#Func
#=----=#
#=----Replays the game until move_id = lastMoveID----=#
while currentMoveID<=lastMoveID
  global board
  res = SQLite.query(db,"""SELECT sourcex,sourcey,targetx,targety,move_type,option FROM moves WHERE "move_number" = $(currentMoveID);""")
  sourcexNullable = res[1,1]
  sourceyNullable = res[1,2]
  targetxNullable = res[1,3]
  targetyNullable = res[1,4]
  move_type = get(res[1,5])
  optionNullable = res[1,6]

  if move_type == "move" #Regular Move
    targetx = get(targetxNullable)
    targety = get(targetyNullable)
    sourcex = get(sourcexNullable)
    sourcey = get(sourceyNullable)
    if !(ST.isEmpty(board[targetx,targety]))# capture
      #push(captures,board[targetx][targety])
    end
    board[targetx,targety].piece = board[sourcex,sourcey].piece
    board[targetx,targety].team = board[sourcex,sourcey].team
    ST.clear!(board[sourcex,sourcey])

  elseif move_type == "drop"
    option = get(optionNullable)
    try
      #deleteat!(captures,findfirst(x->x.piece==option))
    end
    targetx = get(targetxNullable)
    targety = get(targetyNullable)
    board[targetx,targety].piece = board[sourcex,sourcey].piece
    board[targetx,targety].team = board[sourcex,sourcey].team
  elseif move_type == "resign"
    #Do nothing
  end
  ST.saveBoard(board)

  if gameType == "standard"
    displaystandard()
  else
    displaymini()
  end
  sleep(1) #Waits half a second before running displaying a new board
  currentMoveID+=1
end

println("Enter b to go backward, or enter f to go forward\nTo exit viewing the display, enter quit")

instruction=readline(STDIN)

function backward()#function to go backward
  currentMoveID-=1
  if currentMoveID<1
    println("This is the first move")
    continue
  end
  Setupboard()
end

while instruction!="quit\n"
  global board
  if instruction=="b\n"
    currentMoveID-=1
    if currentMoveID<1
      println("This is the first move")
      continue
    end
    #Set up the team and piece of the board#
    Setupboard()
    if gameType == "standard"
      displaystandard()
    else
      displaymini()
    end
    sleep(1)
  elseif instruction=="f\n"
    currentMoveID+=1
    if currentMoveID>lastMoveID
      println("Exceed the last move")
      continue
    end
    Setupboard()
    if gameType == "standard"
      displaystandard()
    else
      displaymini()
    end
    sleep(1)
  end
  instruction=readline(STDIN)
end
export backward, Setupboard

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

"""
    guyan_reduction(K, m, s)

Reduce the stiffness matrix by Guyan Reduction. K = original stiffness matrix, m = master nodes, s= slave nodes.

"""
function guyan_reduction(K, m, s)
    Kred = K[m,m] - K[m,s]*inv(K[s,s])*K[s,m]
    return Kred
end


text:
## Samplers using Distributions package

using Distributions

# single sample
categorical_sample(p::Vector{Float64}) = categorical_sample(p, 1)[:,1]
poisson_sample(lambda::Float64) = poisson_sample(lambda, 1)[1]
gamma_sample(a::Float64, b::Float64) = gamma_sample(a, b, 1)[1]
dirichlet_sample(alpha::Vector{Float64}) = dirichlet_sample(alpha, 1)[:,1]

# samplers
function categorical_sample(p::Vector{Float64}, N::Int)
    K = length(p)
    Z = zeros(K, N)
    Z_tmp = rand(Categorical(p), N)
    for k in 1 : K
        Z[k,find(Z_tmp.==k)] = 1
    end
    return Z
end

function gw_sample(beta::Float64, m::Vector{Float64}, nu::Float64, W::Matrix{Float64})
    # sample Lambda
    Lambda = rand(Wishart(nu, W))
    # sample mu
    mu =  rand(MvNormal(m, inv(beta*Lambda)))
    return mu, Lambda
end

function poisson_sample(lambda::Float64, N::Int)
    X = rand(Poisson(lambda), N)
    return X
end

function gamma_sample(a::Float64, b::Float64, N::Int)
    # a: shape parameter
    # b: rate parameter
    lambda = rand(Gamma(a, 1 / b), N)
    return lambda
end

function dirichlet_sample(alpha::Vector{Float64}, N::Int)
    phi = zeros(length(alpha), N)
    for n in 1 : N
        phi[:,n] = rand(Dirichlet(alpha))
    end
    return phi
end

max_stars_repo_name: sammy-suyama/mlblog

max_issues_repo_name: sammy-suyama/mlblog

max_forks_repo_name: sammy-suyama/mlblog


text:
function rotate(value, wp_x, wp_y)
    if value == 90
        wp_x, wp_y = -wp_y, wp_x
    elseif value == 180
        wp_x, wp_y = -wp_x, -wp_y
    elseif value == 270
        wp_x, wp_y = wp_y, - wp_x
    else
        error("unknown value")
    end
    return wp_x, wp_y
end

function run(s)
    # Your code here
    x::Int, y::Int = 0, 0
    wp_x::Int, wp_y::Int = 10, 1
    for line in readlines(IOBuffer(s))
        action = line[1]
        value = parse(Int, line[2:end])
        if action == 'N'
            wp_y += value
        elseif action == 'E'
            wp_x += value
        elseif action == 'S'
            wp_y -= value
        elseif action == 'W'
            wp_x -= value
        elseif action == 'L'
            wp_x, wp_y = rotate(value % 360, wp_x, wp_y)
        elseif action == 'R'
            wp_x, wp_y = rotate((-value + 360) % 360, wp_x, wp_y)
        elseif action == 'F'
            x += value * wp_x
            y += value * wp_y
        end
    end
    return abs(y) + abs(x)
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
module RemoteSensing

#types
export Landsat7Image
export Landsat8Image
export Sentinel2AImage
export Sentinel2BImage

#spectral indices
export NDVI, TasseledCap

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ b1ed3812-8cce-11eb-14ae-2f44e7ea24a2
begin
	using Lale
	using InteractiveUtils
    using Random
	using Statistics
	using Test
	using DataFrames: DataFrame
	using AutoMLPipeline: Utils
end

# ╔═╡ d27e91ae-8cce-11eb-1e68-cd8f37cb9a8b
begin
	iris = getiris()
	Xreg = iris[:,1:3] |> DataFrame
	Yreg = iris[:,4]   |> Vector
	Xcl  = iris[:,1:4] |> DataFrame
	Ycl  = iris[:,5]   |> Vector
end

# ╔═╡ c7756074-8cce-11eb-2698-91e4400a0b92
begin
	pca     = LalePreprocessor("PCA")
	rb      = LalePreprocessor("RobustScaler")
	noop    = LalePreprocessor("NoOp")
	rfr     = LaleLearner("RandomForestRegressor")
	rfc     = LaleLearner("RandomForestClassifier")
	treereg = LaleLearner("DecisionTreeRegressor")
	ohe  = OneHotEncoder()
	catf = CatFeatureSelector()
	numf = NumFeatureSelector()
end

# ╔═╡ 20cd74ea-8ccf-11eb-2f55-55aa24161f1e
begin
	# regression
	lalepipe =  (pca + noop) >>  (rfr | treereg )
	lale_hopt = LaleOptimizer(lalepipe,"Hyperopt",max_evals=10,cv=3)
	lalepred = fit_transform!(lale_hopt,Xreg,Yreg)
	lalermse=score(:rmse,lalepred,Yreg)
end

# ╔═╡ 3fa9bd06-8ccf-11eb-2d9d-159d5a7cc2b4
begin
	amlpipe = @pipeline  (pca + noop) |> (rfr * treereg)
	amlpred = fit_transform!(amlpipe,Xreg,Yreg)
	crossvalidate(amlpipe,Xreg,Yreg,"mean_squared_error")
	amlprmse=score(:rmse,amlpred,Yreg)
end

# ╔═╡ 52c60cfa-8ccf-11eb-3eaf-3dd8d8f59be1
begin
	# classification
	lalepipec =  (pca) |> rfc
	lale_hoptc = LaleOptimizer(lalepipec,"Hyperopt",max_evals = 10,cv = 3)
	lalepredc  = fit_transform!(lale_hoptc,Xcl,Ycl)
	laleaccc   = score(:accuracy,lalepredc,Ycl)
end

# ╔═╡ 5d79e374-8ccf-11eb-0c48-03cb707018c6
begin
	amlpipec = @pipeline  (pca + rb) |> rfc
	amlpredc = fit_transform!(amlpipec,Xcl,Ycl)
	crossvalidate(amlpipec,Xcl,Ycl,"accuracy_score")
	amlpaccc = score(:accuracy,amlpredc,Ycl)
end

# ╔═╡ Cell order:
# ╠═b1ed3812-8cce-11eb-14ae-2f44e7ea24a2
# ╠═d27e91ae-8cce-11eb-1e68-cd8f37cb9a8b
# ╠═c7756074-8cce-11eb-2698-91e4400a0b92
# ╠═20cd74ea-8ccf-11eb-2f55-55aa24161f1e
# ╠═3fa9bd06-8ccf-11eb-2d9d-159d5a7cc2b4
# ╠═52c60cfa-8ccf-11eb-3eaf-3dd8d8f59be1
# ╠═5d79e374-8ccf-11eb-0c48-03cb707018c6


text:
# struct Arpack_IRAM{Tv} <: AbstractEigenMethod{Tv}
#     precond::Vector{Tv}
# end

# (::Type{<:Arpack_IRAM})(h::AbstractMatrix{Tv}) where {Tv} = Arpack_IRAM(rand(Tv, size(h, 2)))

# function (d::Diagonalizer{<:Arpack_IRAM,Tv})(nev::Integer; precond = true, kw...) where {Tv}
#     if isfinite(d.point)
#         which = :LM
#         # sigma = real(Tv) === Tv ? d.point : d.point + 1.0im
#         sigma = d.point
#     elseif point > 0
#         which = :LR
#         sigma = nothing
#     else
#         which = :SR
#         sigma = nothing
#     end
#     λs, ϕs, _ = Arpack.eigs(d.matrix; nev = nev, sigma = sigma, which = which,
#                              v0 = d.method.precond, kw...)
#     if precond
#         d.method.precond .= zero(Tv)
#         foreach(ϕ -> (d.method.precond .+= ϕ), eachcol(ϕs))
#     end
#     return Eigen(real(λs), ϕs)
# end

ArpackPackage(; sigma = 1im, kw...) = ArpackPackage((sigma = sigma, values(kw)...))


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

using Markdown
using InteractiveUtils

# ╔═╡ 736da374-bfee-11eb-23a7-f1ad290a242e
using Speckles

# ╔═╡ feeaa93e-1a23-4cfc-b219-701531191d10
begin
	# specify all parameters
	νHα1 = [456808] #GHz
	νHα2 = [456808,456811] #GHz
	νHα = [456812, 456808,456811, 456802] #GHz

	paramDict = Dict(
					:n    => [10], # number of atoms
					:νm   => [νHα1,νHα2], # line frequencies in GHz
					:Em   => ["ones"], # relative line magnitudes
					:σ    => [20.0], # Doppler broadening in GHz
					:fγ   => [2*1e4/0.01], # mean photon count rate in GHz
					:deadtime   => [0.00],#,0.01,0.025,0.05,0.1], # detector deadtime in nanoseconds
					:resolution => [0.010], # detector resolution in nanoseconds
					:jitter     => [0.015], # detector timing jitter in nanoseconds
					:efficiency => [0.9], # detector efficiency
					:darkcounts => [1.0e-8], # detector dark count rate in GHz
					:duration   => [20.0], # duration of each correlation measurement in nanoseconds
					:repeat     => [10], # number of times to repeat correlation measurement
					:reinstance => [true], # control whether or not frequencies and phases should be reinstanced between measurements
					:timeint    => ["halfwindow"] # time over which to average correlations in nanoseconds
					# :directory  => [] # defaults to main package directory
					)

end

# ╔═╡ 5e100134-d913-44a5-ad2b-55d29e889304
# split into vector of dictionaries: one for each run
iterParams = paramVector(paramDict);

# ╔═╡ 4d9ab06d-6bae-4066-9333-4287a948112e
# define beamsplitter
bs = Beamsplitter(0.5,0.5)

# ╔═╡ adc32c75-b46f-4818-8de8-50d9990b26b0
for i=1:length(iterParams)
	iterParams[i][:Em] = ones(length(iterParams[i][:νm]))
	iterParams[i][:timeint] = iterParams[i][:duration]/2
end

# ╔═╡ d237c492-4bc3-45a7-a3b5-63a33988a94f
prefixes = map(x->makeName(x),iterParams)

# ╔═╡ e893f820-4be5-4810-84a9-b5629119fe1f
function makeSource(params::Dict)
	return LightSource(
				params[:n],
				params[:Em],
				params[:νm],
				params[:σ],
				params[:fγ]
			)
end

# ╔═╡ 4a19ac54-9617-4a29-a9e6-04cb91450470
source = makeSource(iterParams[1])

# ╔═╡ f4f1c6d6-2b3a-41b0-88e5-d0a59debed04
function makeDetector(params::Dict)
	return Detector(
            params[:deadtime],
            params[:resolution],
            params[:jitter],
            params[:efficiency],
            params[:darkcounts]
        )
end

# ╔═╡ dac6bbdc-8996-425b-8410-f789b71ea3d9
detect = makeDetector(iterParams[1])

# ╔═╡ 65bf052e-43c6-4743-b264-b04ae9843a2c
γint = γIntensity(nbar(iterParams[1][:duration],source)*bs.t,
	iterParams[1][:duration],
	detect.resolution,
	source)

# ╔═╡ b10ddbe6-90eb-4e65-a972-d20f52ed7547

# ╔═╡ Cell order:
# ╠═736da374-bfee-11eb-23a7-f1ad290a242e
# ╠═feeaa93e-1a23-4cfc-b219-701531191d10
# ╠═5e100134-d913-44a5-ad2b-55d29e889304
# ╠═4d9ab06d-6bae-4066-9333-4287a948112e
# ╠═adc32c75-b46f-4818-8de8-50d9990b26b0
# ╠═d237c492-4bc3-45a7-a3b5-63a33988a94f
# ╠═e893f820-4be5-4810-84a9-b5629119fe1f
# ╠═4a19ac54-9617-4a29-a9e6-04cb91450470
# ╠═f4f1c6d6-2b3a-41b0-88e5-d0a59debed04
# ╠═dac6bbdc-8996-425b-8410-f789b71ea3d9
# ╠═65bf052e-43c6-4743-b264-b04ae9843a2c
# ╠═b10ddbe6-90eb-4e65-a972-d20f52ed7547


text:
#=
Metaprogramming: Julia Cookbook
=#

function new_struct(fields::Vector{Tuple{String,DataType}})
    name = "A" * string(hash(fields), base = 16)
    code = "begin\nstruct $name\n"
    for f in fields
        code *= f[1] * "::" * string(f[2]) * "\n"
    end
    code *= "end\n$name\nend"
    eval(Meta.parse(code))
end

function parse_data(data::AbstractString)
    lines = filter(x -> length(x) > 0, strip.(split(data, ('\n', '\r'))))
    colnames = string.(split(lines[1], ','))
    row1 = split.(lines[2], ',')
    coltypes = [occursin(r"^-?\d+$", val) ? Int64 : String for val in row1]
    (lines[2:end], new_struct(collect(zip(colnames, coltypes))))
end

function parse_text(data::AbstractString)
    lines, MyStruct = parse_data(data)
    res = MyStruct[]
    for line in lines
        colvals = split(line, ',')
        f = (t, v) -> t <: Int ? parse(Int, v) : string(v)
        vals = f.(MyStruct.types, colvals)
        push!(res, Base.invokelatest(MyStruct, vals...))
    end
    res
end

MyS = new_struct([("a", Int), ("b", String), ("c", Int)])
#repr(dump(MyS))
#repr(dump(parse_data("col1,col2,col3\nabc,123,123.5")[2]))
data = """
id,val,class
1,4,A
2,39,B
3,44,C
"""
#print(parse_text(data))

code = "begin\n
      struct Goma\n
          x::Int\n
          y::Int\n
      end\n
      Goma\n
  end"

g = eval(Meta.parse(code))(1, 2)
print(g) # should be Goma(1,2)


text:
export TraceTape, recover_ast

struct TraceTape{B} <: AbstractRegister{B}
    commands::NTuple{B,Vector{Any}}
end

TraceTape() = TraceTape{1}((Any[],))

function trace!(tape::TraceTape{1}, stmt)
    push!(tape.commands[1], stmt)
end

function Base.show(io::IO, tape::TraceTape{1})
    for ex in tape.commands[1]
        printstyled(io, ex.args[1]; color = :light_blue, bold = true)
        print(io, " ")
        join(io, string.(ex.args[2:end]), " ")
        println(io)
    end
end

function recover_ast(tape::TraceTape{1})
    ex = Expr(:block)
    for each in tape.commands[1]
        if each.args[1] === :gate
            push!(ex.args, Expr(:call, :(=>), each.args[3], each.args[2]))
        elseif each.args[1] === :ctrl
            push!(
                ex.args,
                Expr(
                    :macrocall,
                    Symbol("@ctrl"),
                    LineNumberNode(@__LINE__, @__FILE__),
                    each.args[3],
                    Expr(:call, :(=>), each.args[3], each.args[2]),
                ),
            )
        elseif each.args[1] === :measure
            # push!(ex.args,
            #     Expr(:macrocall, Symbol("@measure"),
            #         LineNumberNode(@__LINE__, @__FILE__),
            #         ex.args...
            #     )
            # )
        end
    end
    return ex
end

function YaoIR(m::Module, tape::TraceTape{1}, name = gensym())
    ast = recover_ast(tape)
    lowered_ast = Meta.lower(m, to_function(m, ast))
    if lowered_ast === nothing
        body = IR()
    else
        body = IR(lowered_ast.args[], 0)
    end

    ir = YaoIR(m, name, Any[], Any[], mark_quantum(body), true, true)
    update_slots!(ir)
    return ir
end

function quantum_m(m::Module, n::Int, ex)
    if !((ex isa Expr) && (ex.head === :call))
        throw(ParseError("expect a call statement"))
    end

    name = gensym(Symbol(ex.args[1]))
    quote
        old_ir = $(esc(IRTools.xcall(YaoLang, :code_yao, ex.args...)))
        if $(GlobalRef(YaoLang, :is_pure_quantum))(old_ir)
            $(esc(ex.args[1]))
        else
            tape = $TraceTape()
            $(esc(ex))(tape, $Locations(1:$n))
            ir = $YaoIR($m, tape, $(QuoteNode(name)))
            $(GlobalRef(Base, :eval))($m, codegen(ir))
            $(esc(name))
        end
    end
end

export @quantum

"""
    @quantum nqubits::Int <generic circuit call>

Partially evaluate all classical parts of the generic circuit.

!!! note
    The classical part of the program could result in
    different circuits if it is not deterministic.
"""
macro quantum(n::Int, ex)
    return quantum_m(__module__, n, ex)
end


text:
# This file was generated, do not modify it. # hide
ytest, Xtest = unpack(df_test, ==(:Churn), !=(:customerID));

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# using Images, ImageView
using LCMCore
using FixedPointNumbers
using ColorTypes
using ImageView

mutable struct image_metadata_t <: LCMType
    key::String
    n::Int32
    value::Vector{UInt8}
end

@lcmtypesetup(image_metadata_t,
    value => (n,)
)

mutable struct image_t <: LCMType
    utime::Int64
    width::Int32
    height::Int32
    row_stride::Int32
    pixelformat::Int32
    size::Int32
    data::Vector{UInt8}
    nmetadata::Int32
    metadata::Vector{image_metadata_t}
end

@lcmtypesetup(image_t,
    data => (size,),
    metadata => (nmetadata,)
)

# these values are to be kept in sync with the CamPixelFormat enumeration
# in the Camunits library.  See camunits/pixels.h in Camunits
const PIXEL_FORMAT_UYVY             = Int32(1498831189)
const PIXEL_FORMAT_YUYV             = Int32(1448695129)
const PIXEL_FORMAT_IYU1             = Int32(827677001)
const PIXEL_FORMAT_IYU2             = Int32(844454217)
const PIXEL_FORMAT_YUV420           = Int32(842093913)
const PIXEL_FORMAT_YUV411P          = Int32(1345401140)
const PIXEL_FORMAT_I420             = Int32(808596553)
const PIXEL_FORMAT_NV12             = Int32(842094158)
const PIXEL_FORMAT_GRAY             = Int32(1497715271)
const PIXEL_FORMAT_RGB              = Int32(859981650)
const PIXEL_FORMAT_BGR              = Int32(861030210)
const PIXEL_FORMAT_RGBA             = Int32(876758866)
const PIXEL_FORMAT_BGRA             = Int32(877807426)
const PIXEL_FORMAT_BAYER_BGGR       = Int32(825770306)
const PIXEL_FORMAT_BAYER_GBRG       = Int32(844650584)
const PIXEL_FORMAT_BAYER_GRBG       = Int32(861427800)
const PIXEL_FORMAT_BAYER_RGGB       = Int32(878205016)
const PIXEL_FORMAT_BE_BAYER16_BGGR  = Int32(826360386)
const PIXEL_FORMAT_BE_BAYER16_GBRG  = Int32(843137602)
const PIXEL_FORMAT_BE_BAYER16_GRBG  = Int32(859914818)
const PIXEL_FORMAT_BE_BAYER16_RGGB  = Int32(876692034)
const PIXEL_FORMAT_LE_BAYER16_BGGR  = Int32(826360396)
const PIXEL_FORMAT_LE_BAYER16_GBRG  = Int32(843137612)
const PIXEL_FORMAT_LE_BAYER16_GRBG  = Int32(859914828)
const PIXEL_FORMAT_LE_BAYER16_RGGB  = Int32(876692044)
const PIXEL_FORMAT_MJPEG            = Int32(1196444237)
const PIXEL_FORMAT_BE_GRAY16        = Int32(357)
const PIXEL_FORMAT_LE_GRAY16        = Int32(909199180)
const PIXEL_FORMAT_BE_RGB16         = Int32(358)
const PIXEL_FORMAT_LE_RGB16         = Int32(1279412050)
const PIXEL_FORMAT_BE_SIGNED_GRAY16 = Int32(359)
const PIXEL_FORMAT_BE_SIGNED_RGB16  = Int32(360)
const PIXEL_FORMAT_FLOAT_GRAY32     = Int32(842221382)
const PIXEL_FORMAT_INVALID          = Int32(-2)
const PIXEL_FORMAT_ANY              = Int32(-1)

#
# def decode(self, data):
#     msg = image_t.decode(data)
#     if msg.pixelformat == image_t.PIXEL_FORMAT_GRAY:
#         return im_resize(np.asarray(bytearray(msg.data),
#                                     dtype=np.uint8).reshape(msg.height,
#                                                             msg.width),
#                          scale=self.scale)
#     elif msg.pixelformat == image_t.PIXEL_FORMAT_MJPEG:
#         im = cv2.imdecode(np.asarray(bytearray(msg.data),
#                                      dtype=np.uint8),
#                           -1)
#         return im_resize(im, scale=self.scale)
#     else:
#     raise RuntimeError('Unknown pixelformat for ImageDecoder')

# global imgt = image_t()

global canvas = Dict{String, Any}()

function updateCamCallback!(channel::String, msg::image_t)
    global canvas
    imc = Array{RGB{N0f8}}(undef, msg.height, msg.width)
    if (msg.pixelformat == PIXEL_FORMAT_RGB)
        im1 = reshape(msg.data, (3, msg.width, msg.height))
        im2 = reinterpret(N0f8, im1)
        for i in 1:msg.height, j in 1:msg.width
          imc[i, j] = RGB(im2[:, j,i]...)
        end
    end
    if length(canvas) == 0
      canvas = imshow(imc)
    else
      imshow(canvas["gui"]["canvas"], imc)
    end
end

# This colour mapping may be wrong for the PS3eye.
# ycrcb = Video4Linux.YUYV(640,480)
# vidchan = Channel((c::Channel) -> videoproducer(c, ycrcb, devicename = "/dev/video0", iomethod = Video4Linux.IO_METHOD_MMAP ))
#
# ##
# #capture one frame to create im1 and canvas needed to diplay
# A =  take!(vidchan)
# im1 = RGB.(YCbCr.(view(A,:,:,1), view(A,:,:,2), view(A,:,:,3)))

global const lcm_handle = LCM("udpm://239.255.76.67:7667?ttl=1")

subscribe(lcm_handle, "IMAGE_STREAM", updateCamCallback!, image_t)

# canvas = imshow(im1)

while handle(lcm_handle)
end


text:
module WCSUtils

import WCS: WCSTransform, world_to_pix, pix_to_world

export pixel_deriv_to_world_deriv, wcs_id

"""
Convert a world location to a 1-indexed pixel location using a linear
approximation to the true transform.

Args:
    - wcs_jacobian: The jacobian of the transform pixel_loc = F(world_loc)
    - world_offset: The world location at which the jacobian is evaluated
    - pix_offset: The pixel location of the world offset
    - worldcoords: The world location to be translated to pixel coordinates.

Returns:
    - The 1-indexed pixel locations in the same shape as the input.
"""
world_to_pix{T <: Number}(wcs_jacobian::Matrix{Float64},
                          world_offset::Vector{Float64},
                          pix_offset::Vector{Float64},
                          worldcoords::VecOrMat{T}) =
    wcs_jacobian * (worldcoords .- world_offset) + pix_offset

"""
A finite-differences approximation to the jacobian of the transform
from pixel to world coordinates.

Args:
 - wcs: The world coordinate system
 - pix_loc: The location at which to evaluate the jacobian
 - pixel_delt: The step size for the finite difference
               approximation in pixel coordinates

Returns:
 - The jacobian of the transform pixel_coord = F(world_coord).  Following the
   standard definition, the pixel coordinates vary across rows and the world
   coordinates across columns.
"""
function pixel_world_jacobian(wcs::WCSTransform, pix_loc::Array{Float64, 1};
	                      pixel_delt=0.5)

    # Choose a step size.
    # Assume that about a half a pixel is a reasonable step size and the
    # directions are about the same.
    world_loc = pix_to_world(wcs, pix_loc)
    world_delt = maximum(abs(
      pix_to_world(wcs, pix_loc + Float64[pixel_delt, pixel_delt]) - world_loc))

    world_loc_1 = world_loc + world_delt * Float64[1, 0]
    world_loc_2 = world_loc + world_delt * Float64[0, 1]

    hcat((world_to_pix(wcs, world_loc_1) - pix_loc) ./ world_delt,
         (world_to_pix(wcs, world_loc_2) - pix_loc) ./ world_delt)
end

"""
Transform a derivative of a scalar function with respect to pixel
coordinates into a derivatve with respect to world coordinates.

Args:
 - wcs: The world coordinate system object
 - df_dpix: The derivative of a scalar function with respect to pixel coordinates
 - pix_loc: The pixel location at which the derivative was taken.
 - pixel_delt: The step size for the finite difference
               approximation in pixel coordinates

Returns:
 - The derivative with respect to world coordinates.
"""
function pixel_deriv_to_world_deriv(
  wcs::WCSTransform, df_dpix::Array{Float64, 1},
  pix_loc::Array{Float64, 1}; pixel_delt=1e-3)
    @assert length(pix_loc) == length(df_dpix) == 2
    trans = pixel_world_jacobian(wcs, pix_loc, pixel_delt=pixel_delt)
    trans' * df_dpix
end

# A world coordinate system where the world and pixel coordinates are the same.
const wcs_id = WCSTransform(2,
                            cd = Float64[1 0; 0 1],
                            ctype = ["none", "none"],
                            crpix = Float64[1, 1],
                            crval = Float64[1, 1]);

end


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

using Markdown
using InteractiveUtils

# ╔═╡ ac89ffd2-f717-11eb-3fb2-d77d7f818725
begin
    import Pkg
    Pkg.activate(Base.current_project())

	using QuantumOptics
end

# ╔═╡ 65dd2c78-3317-4fc4-867b-ae5a1043f3f3
begin
	using PyPlot
	ion()
end

# ╔═╡ ccd42b7b-0c0c-4bb8-b3f0-fa424bbf127c
begin
	const dim = 50

	fbasis = FockBasis(dim)
end

# ╔═╡ e0fbeeb5-02ef-4c4f-a8cf-bb8143af4660
begin
	a = destroy(fbasis)
	n = number(fbasis)
end;

# ╔═╡ 97e730c1-c09e-44dd-8074-283b993394f0
H(ω, λ) = ω * dagger(a) * a + λ * (a + dagger(a))

# ╔═╡ 05531381-00fa-424e-87cd-fab529b72d3d
tlist = [0:0.05:10;];

# ╔═╡ 1bed70a4-27d2-4a7b-bd49-ab565ea128c0
ψ₀(α) = coherentstate(fbasis, α);

# ╔═╡ b7cb690e-3e0e-4855-97af-5015ac354fa2
tout, ψₜ = timeevolution.schroedinger(tlist, ψ₀(4), H(4, 0.7));

# ╔═╡ cc571d2a-3abe-4c07-a64d-9799b5aef844
begin
	ev1 = expect(n, ψₜ)
	ev2 = expect(dagger(a)+a, ψₜ)
end;

# ╔═╡ 45e61eba-9dbe-4322-8924-155058f90a02
begin
	fig, (ax1, ax2) = subplots(2, 1, figsize = (6,10))
	ax1.plot(tout, ev1)
	ax2.plot(tout, ev2)
end

# ╔═╡ Cell order:
# ╠═ac89ffd2-f717-11eb-3fb2-d77d7f818725
# ╠═65dd2c78-3317-4fc4-867b-ae5a1043f3f3
# ╠═ccd42b7b-0c0c-4bb8-b3f0-fa424bbf127c
# ╠═e0fbeeb5-02ef-4c4f-a8cf-bb8143af4660
# ╠═97e730c1-c09e-44dd-8074-283b993394f0
# ╠═05531381-00fa-424e-87cd-fab529b72d3d
# ╠═1bed70a4-27d2-4a7b-bd49-ab565ea128c0
# ╠═b7cb690e-3e0e-4855-97af-5015ac354fa2
# ╠═cc571d2a-3abe-4c07-a64d-9799b5aef844
# ╠═45e61eba-9dbe-4322-8924-155058f90a02

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
struct Queue <: LavaAbstraction
    handle::Vk.Queue
    capabilities::Vk.QueueFlag
    index::Int
    family::Int
end

vk_handle_type(::Type{Queue}) = Vk.Queue

struct QueueDispatch
    queues::Dictionary{Int,Vector{Queue}}
    present_queue::Optional{Queue}
    """
    Build a `QueueDispatch` structure from a given device and configuration.
    If a surface is provided, then a queue that supports presentation on this surface will be filled in `present_queue`.

    !!! warning
        `device` must have been created with a consistent number of queues as requested in the provided queue configuration.
        It is highly recommended to have created the device with the result of `queue_infos(QueueDispatch, physical_device, config)`.
    """
    function QueueDispatch(device, infos; surface = nothing)
        pdevice = physical_device(device)
        families = dictionary(map(infos) do info
            info.queue_family_index => length(info.queue_priorities)
        end)
        props = Vk.get_physical_device_queue_family_properties(pdevice)
        queues = dictionary(map(pairs(families)) do (family, count)
            family => map(1:count) do index
                info = Vk.DeviceQueueInfo2(family, index - 1)
                Queue(Vk.get_device_queue_2(device, info), props[family + 1].queue_flags, index, family + 1)
            end
        end)
        present_queue = if !isnothing(surface)
            idx = findfirst(families) do family
                unwrap(Vk.get_physical_device_surface_support_khr(pdevice, family, surface))
            end
            if isnothing(idx)
                error("Could not find a queue that supports presentation on the provided surface.")
            else
                first(values(queues)[idx])
            end
        else
            nothing
        end
        new(queues, present_queue)
    end
end

function queue_infos(::Type{QueueDispatch}, physical_device::Vk.PhysicalDevice, config)
    # queue family index => count
    families = Dictionary{Int,Int}()
    props::Vector{Union{Nothing,Vk.QueueFamilyProperties}} = Vk.get_physical_device_queue_family_properties(physical_device)
    config = deepcopy(config)

    # resolve exact matches first
    for (i, prop) in enumerate(props); isnothing(prop) && continue
        for (flags, count) in pairs(config)
            # exact match
            if flags == prop.queue_flags
                remaining = prop.queue_count - get(families, i-1, 0)
                if remaining ≥ count
                    delete!(config, flags)
                    set!(families, i-1, get(families, i-1, 0) + count)
                else
                    config[flags] -= remaining
                    set!(families, i-1, get(families, i-1, 0) + remaining)
                    props[i-1] = nothing
                    break
                end
            end
        end
    end

    for (i, prop) in enumerate(props); isnothing(prop) && continue
        for (flags, count) in pairs(config)
            # match
            if flags in prop.queue_flags
                remaining = prop.queue_count - get(families, i-1, 0)
                if remaining ≥ count
                    delete!(config, flags)
                    set!(families, i-1, get(families, i-1, 0) + count)
                else
                    config[flags] -= remaining
                    set!(families, i-1, get(families, i-1, 0) + remaining)
                    props[i-1] = nothing
                    break
                end
            end
        end
    end

    map(pairs(families)) do (index, count)
        Vk.DeviceQueueCreateInfo(index, ones(count))
    end |> collect
end

function get_queue_family(dispatch::QueueDispatch, properties::Vk.QueueFlag)
    # try exact match
    for (family, queues) in pairs(dispatch.queues)
        properties == first(queues).capabilities && return family
    end

    # try queues that contain the required properties
    for (family, queues) in pairs(dispatch.queues)
        properties in first(queues).capabilities && return family
    end

    # panic
    error("Could not find a queue matching with the required properties $properties")
end

function queue(dispatch::QueueDispatch, family_index)
    haskey(dispatch.queues, family_index) && return first(dispatch.queues[family_index])
    !isnothing(dispatch.present_queue) && dispatch.present_queue.family == family_index && return dispatch.present_queue
    error("Could not find queue with family index $family_index")
end

function present(dispatch::QueueDispatch, present_info::Vk.PresentInfoKHR)
    queue = dispatch.present_queue
    if isnothing(queue)
        error("No presentation queue was specified for $dispatch")
    else
        Vk.queue_present_khr(queue, present_info)
    end
end

function queue_family_indices(dispatch::QueueDispatch; include_present = true)
    indices = collect(keys(dispatch.queues))
    !isnothing(dispatch.present_queue) && include_present && union!(indices, dispatch.present_queue.family)
    indices
end