text:
# Methods related to online control of dynamic systems

"""
Objects relevant to online control problem
"""
struct OnlineControlProblem
    initial_state
    control_modes
    # TODO: We don't need this IMDP here, get rid of it.
    imdp::IMDP
    pimdp::PIMDP
    dfa::DFA
    extents
    post_image_data
    control_policy
    local_gps::Dict
    global_gps::Dict # stores the paths to the full GPs generated offline
    reward_type::Tuple
    # TODO: There is probably a better place for this somewhere else!
    primary_reward_history
    secondary_reward_history
    online_data # Stores the data collected online
    global_gp_infos # Stores the persistent info dicts
end

"""
Setup the online control problem.
"""
function setup_online_control_problem(x0, modes, imdp, pimdp, dfa, extents, posts, policy, local_gps, global_gps, reward_type)
    reset_pimdp(x0, imdp, dfa, pimdp, extents, policy)

    # Setup local gp data
    online_data = Dict()
    for i=1:length(modes)
        online_data[i] = Dict("x_online" => [], "y_online" => Dict(i=>[] for i in keys(extents[1])))
    end
    # online_data["x_online"] = []
    # online_data["y_online"] = Dict(i=>[] for i in keys(extents[1]))

    global_gp_infos = Dict()

    return OnlineControlProblem(x0, modes, imdp, pimdp, dfa, extents, posts, policy, local_gps, global_gps, reward_type, [], [], online_data, global_gp_infos)
end

"""
One step of a control loop with an offline policy.
"""
function online_control_loop(ocp::OnlineControlProblem)
    res = check_pimdp_exit_conditions(ocp.pimdp, length(ocp.dfa.states))
    if !isnothing(res)
        return res
    end
    x_c = ocp.pimdp.trajectory[end]
    mode = ocp.pimdp.policy[map_pimdp_state_to_index(ocp.pimdp.state_history[end], length(ocp.dfa.states) )]
    x_n = ocp.control_modes[Int(mode)](x_c)
    propogate_pimdp_trajectory(ocp.pimdp, ocp.dfa, ocp.extents, x_n)
    return nothing
end

"""
Concatonate Online and Offline Data
"""
function concatonate_data_points(x_data_offline, y_data_offline, x_data_online, y_data_online)
    x_cat = [x_data_offline x_data_online...]
    y_cat = [y_data_offline; y_data_online]
    return x_cat, y_cat
end

using NearestNeighbors
"""
Create local GP for an extent.
"""
function get_local_data_knn(center, x_data, y_data; num_neighbors = 50)
    num_neighbors = minimum([num_neighbors, size(x_data,2)])
    kdtree = KDTree(x_data);
    sub_idx, _ = knn(kdtree, center, num_neighbors, true)

    # WHY THE HELL IS IT DOING THIS
    if typeof(sub_idx) == Array{Array{Int64,1},1}
        sub_idx = sub_idx[1]
    end

    sub_x_data = x_data[:, sub_idx]
    sub_y_data = y_data[sub_idx]
    return sub_x_data, sub_y_data
end

"""
One step optimal control loop with an offline policy and GP update.
"""
function online_control_loop_optimal(ocp::OnlineControlProblem;
                                     update_gp_flag=false, gp_use_key="global", update_pimdp_ring_states=false,
                                     pimdp_reward_dict=nothing)
    etime = @elapsed begin
        # Check if an exit condition has been reached
        res = check_pimdp_exit_conditions(ocp.pimdp, length(ocp.dfa.states))
        if !isnothing(res)

            return res
        end

        x_c = ocp.pimdp.trajectory[end]     # Current real state
        q_c = ocp.pimdp.state_history[end]
        pimdp_state = map_pimdp_state_to_index(q_c, length(ocp.dfa.states))  # Current PIMDP state index

        best_mode, best_reward, gp_info_dict = get_best_action(ocp, gp_use_key=gp_use_key, pimdp_reward_dict=pimdp_reward_dict)
        x_n = ocp.control_modes[best_mode](x_c)

        propogate_pimdp_trajectory(ocp.pimdp, ocp.dfa, ocp.extents, x_n)

        @debug "BEST ACTION: " best_mode
        @debug "BEST REWARD: " best_reward

        push!(ocp.primary_reward_history, best_reward[1])
        push!(ocp.secondary_reward_history, best_reward[2])

        # TODO: Move this to its own fcn / subfcns
        if update_gp_flag && isnothing(check_pimdp_exit_conditions(ocp.pimdp, length(ocp.dfa.states)))
            dim_keys = keys(gp_info_dict)
            push!(ocp.online_data[best_mode]["x_online"], x_c)
            # Update the GP with a new datapoint
            for (i, dim_key) in enumerate(dim_keys)
                y_meas = (x_n[i] - x_c[i])
                push!(ocp.online_data[best_mode]["y_online"][dim_key], y_meas)
                add_datapoint_to_gp(gp_info_dict[dim_key].gp, (x_c, y_meas))
                update_gp_info(gp_info_dict[dim_key])
            end

            pimdp_states_to_update = [q_c, ]
            if update_pimdp_ring_states
                new_states = get_pimdp_ring_states(q_c)

                for new_state in new_states
                    push!(pimdp_states_to_update, new_state)
                end
            end

            for pimdp_state in pimdp_states_to_update

                # Update transition probability to some states
                extent_id = pimdp_state[1]
                regions_to_update = []
                # pimdp_row = map_pimdp_state_to_Pmat_row(q_c, length(ocp.dfa.states), , best_mode)
                pimdp_rows = get_pimdp_row_idx_same_imdp(pimdp_state, ocp.dfa, length(ocp.pimdp.actions), best_mode)
                @assert length(pimdp_rows) > 0

                # og_post = ocp.post_image_data[best_mode][pimdp_state[1]]

                # Get the new post-image of the GP
                current_extent = ocp.extents[extent_id]
                lb = [current_extent[dim_key][1] for dim_key in dim_keys]
                ub = [current_extent[dim_key][2] for dim_key in dim_keys]

                ## TODO: THIS IS MESSY
                region_post_dict_new = bound_extent(current_extent, lb, ub, gp_info_dict, dim_keys, Dict("x1" => [1., 1.], "x2" => [1., 1.]); known_part_flag=true)
                dummy_region_post_dict = Dict(extent_id => region_post_dict_new)

                for pimdp_row_idx in pimdp_rows
                    prow_min = ocp.pimdp.Pmin[pimdp_row_idx, :]
                    prow_max = ocp.pimdp.Pmax[pimdp_row_idx, :]

                    # Find all states with interval bounds
                    regions_to_update = findall(x->x>0, prow_min .!= prow_max)
                    pimdp_row_vals_min = prow_min[regions_to_update]
                    pimdp_row_vals_max = prow_max[regions_to_update]

                    noise_proc = Truncated(Normal(0, 0.01), -0.01, 0.01)
                    η = 0.01
                    ϵ = -1

                    for target in regions_to_update
                        target_id = pimdp_col_to_imdp_state(target, length(ocp.dfa.states))
                        region_pair = (extent_id, target_id)
                        frame_row = region_pair_transitions(region_pair, ocp.extents, dummy_region_post_dict, ϵ, gp_info_dict, noise_proc, η)

                        if target_id == 1025.
                            minPrn = 1 - frame_row[4]
                            maxPrn = 1 - frame_row[3]

                            if minPrn > prow_min[target]
                                prow_min[target] = minPrn
                            end
                            if maxPrn < prow_max[target]
                                prow_max[target] = maxPrn
                            end
                        else
                            if frame_row[4] < prow_max[target]
                                prow_max[target] = frame_row[4]
                            end
                            if frame_row[3] > prow_min[target]
                                prow_min[target] = frame_row[3]
                            end
                        end
                    end

                    ocp.pimdp.Pmin[pimdp_row_idx, :] = prow_min
                    ocp.pimdp.Pmax[pimdp_row_idx, :] = prow_max
                end
            end
        end
    end
    @debug "Total time: $etime seconds"
end

"""
Update PIMDP row
"""
function update_pimdp_transition_bound_row()
    nothing
end

"""
Get the relevant GP info
"""
function get_gp_info(ocp, mode, current_state; gp_use_case="global")
    # Begin all of the fancy local GP processing
    global_gp_key = "global-$mode"

    # gp_use_case = "use-local-fresh"
    if gp_use_case == "global"
        if global_gp_key in keys(ocp.global_gp_infos)
            global_gp_info_dict = ocp.global_gp_infos[global_gp_key]
        else
            gp_filename = ocp.global_gps[mode]
            gp_info = load_gp("$gp_filename")
            global_gp_info_dict = gp_info[:gp_info]
            ocp.global_gp_infos[global_gp_key] = global_gp_info_dict
        end
        gp_info_dict = global_gp_info_dict
    elseif gp_use_case == "local-persistent"
        nothing
        # gp_key = "$imdp_state-$best_mode"
        # if gp_key in keys(ocp.local_gps)
        #     gp_info_dict = ocp.local_gps[gp_key]
        #     # gp_info_dict
        # else
        #     #gp_filename = filter(x -> occursin("-$imdp_state.bin", x), readdir(gp_dir))[1]
        #     # Load the global GP this could be optimized
        #     gp_filename = ocp.global_gps[best_mode]
        #     gp_info = load_gp("$gp_filename")
        #     # gp_info = load_gp("$gp_dir/$gp_filename")
        #     # gp_local = create_local_gp_knn(x_c, gp_info; num_neighbors = 100)
        #     gp_info_dict = gp_info[:gp_info]
        # end

    elseif gp_use_case == "local-fresh"
        if global_gp_key in keys(ocp.global_gp_infos)
            global_gp_info_dict = ocp.global_gp_infos[global_gp_key]
        else
            gp_filename = ocp.global_gps[mode]
            gp_info = load_gp("$gp_filename")
            global_gp_info_dict = gp_info[:gp_info]
            ocp.global_gp_infos[global_gp_key] = global_gp_info_dict
        end

        # Create a fresh local GP from the global GP
        gp_info_dict = copy(global_gp_info_dict)
        x_offline = global_gp_info_dict["x1"].gp.x
        x_online = ocp.online_data[mode]["x_online"]

        for dim_key in keys(gp_info_dict)
            y_offline = global_gp_info_dict[dim_key].gp.y
            y_online = ocp.online_data[mode]["y_online"][dim_key]
            x_cat, y_cat = concatonate_data_points(x_offline, y_offline, x_online, y_online)
            x_nn, y_nn = get_local_data_knn(current_state, x_cat, y_cat, num_neighbors=75)

            gp_info_dict[dim_key].gp = train_gp_1dim(x_nn, y_nn; se_params=[0., 0.65])
            update_gp_info(gp_info_dict[dim_key])
        end
    else
        @error "Unknown GP use case: " + gp_use_case
    end

    return gp_info_dict
end

"""
Get the policy from a verification result.
"""
function get_policy_from_offline_synthesis(filename)
    policy = Dict()
    open(filename, "r") do f
        while !eof(f)
            row_split = split(readline(f))
            if length(row_split) > 0
                q0 = parse(Int, row_split[1]) + 1
                u = parse(Int, row_split[2]) + 1
                policy[q0] = u
            end
        end
    end
    return policy
end

"""
Gets the best action at each state.
"""
function get_best_action(ocp::OnlineControlProblem; gp_use_key="global", pimdp_reward_dict=nothing)
    # Get the current state and PIMDP state
    x_c = ocp.pimdp.trajectory[end]     # Current real state
    q_c = ocp.pimdp.state_history[end]  # Current IMDP (extent ID) / DFA state
    pimdp_state = map_pimdp_state_to_index(q_c, length(ocp.dfa.states))  # Current PIMDP state index

    # Initialize the rewards and best control mode
    best_reward = [-Inf, -Inf]
    best_ub = -Inf
    best_mode = nothing
    best_gp_info_dict = nothing

    # Check if the current policy gives us a shortcut
    if ocp.control_policy[pimdp_state, 2] == 1.0
        best_mode = Int(ocp.control_policy[pimdp_state, 1])
        best_reward[1] = 2
        best_gp_info_dict = get_gp_info(ocp, best_mode, x_c, gp_use_case=gp_use_key)
    else
        for mode = 1:length(ocp.control_modes)

            gp_info_dict = get_gp_info(ocp, mode, x_c, gp_use_case=gp_use_key)
            # Extract the relevant row of the PIMDP transition matrix
            P_row_idx = map_pimdp_state_to_Pmat_row(q_c, length(ocp.dfa.states), length(ocp.pimdp.actions), mode)
            prow_max = ocp.pimdp.Pmax[P_row_idx, :]

            # Find all states with interval bounds with non-zero UB
            regions_with_nonzero_ub = findall(x->x>0., prow_max)
            @debug regions_with_nonzero_ub

            ##### Test out the point region transition fcn
            ϵ = -1.
            μ, σ = predict_one_step_full(gp_info_dict, x_c)
            # TODO: get rid of this noise dist implicit
            noise_sigma = 0.01
            dist = Truncated(Normal(0., noise_sigma), -noise_sigma, noise_sigma)

            non_trivial_regions = []
            prob_intervals = Dict()
            imdp_extent_ids = []
            nt_pimdp_states = []

            for reg in regions_with_nonzero_ub
                region_idx = pimdp_col_to_imdp_state(reg, length(ocp.dfa.states))
                push!(nt_pimdp_states, (region_idx, q_c[2]))
                push!(imdp_extent_ids, region_idx)
                res = point_region_transitions(μ, σ, region_idx, ocp.extents, ϵ, gp_info_dict; process_noise_dist=dist)
                if res[2] > 0.
                    push!(non_trivial_regions, map_pimdp_state_to_index((Int(region_idx), q_c[2]), length(ocp.dfa.states) ))
                    prob_intervals[map_pimdp_state_to_index((Int(region_idx), q_c[2]), length(ocp.dfa.states) )] = [res[1], res[2]]
                end
            end

            # Extract the potential actions based on the potential regions we could transition to.
            potential_actions = zeros(length(non_trivial_regions), 6)
            total_reward = 0.

            max_reward = 0.
            for (i,region) in enumerate(non_trivial_regions)
                potential_actions[i, 1] = region
                potential_actions[i, 2:4] = ocp.control_policy[region,:]
                potential_actions[i, 5] = prob_intervals[region][1]
                potential_actions[i, 6] = prob_intervals[region][2]
                # total_reward += ocp.control_policy[region,2]
                # total_reward += ocp.control_policy[region,3]
            end
            # total_reward *= 1. / length(non_trivial_regions)

            # Calculate the worst expected lower bound of satisfaction
            sorted_actions = potential_actions[sortperm(potential_actions[:,3]), :]
            region_idx = sorted_actions[:,1]
            minPrs = sorted_actions[:,5]
            maxPrs = sorted_actions[:,6]
            worst_probs = get_true_transition_probabilities(minPrs, maxPrs, collect(1:1:length(region_idx)))

            worst_exp_lb = sum(worst_probs .* sorted_actions[:, 3]')

            sorted_actions2 = potential_actions[sortperm(potential_actions[:,4]), :]
            region_idx = sorted_actions2[:,1]
            minPrs = sorted_actions2[:,5]
            maxPrs = sorted_actions2[:,6]
            worst_probs2 = get_true_transition_probabilities(minPrs, maxPrs, collect(1:1:length(region_idx)))
            # worst_exp_ub = sum(worst_probs2 .* sorted_actions[:, 4]')
            worst_exp_ub = sum(sorted_actions[:,4]')/length(sorted_actions[:,4]')

            if ocp.reward_type[1] == "worst-exp-lb"
                total_reward = worst_exp_lb
            elseif ocp.reward_type[1] == "worst-exp-sum"
                sorted_actions2 = potential_actions[sortperm(potential_actions[:,4]), :]
                region_idx = sorted_actions2[:,1]
                minPrs = sorted_actions2[:,5]
                maxPrs = sorted_actions2[:,6]
                worst_probs2 = get_true_transition_probabilities(minPrs, maxPrs, collect(1:1:length(region_idx)))
                total_reward = 0.5*(sum(worst_probs .* sorted_actions[:, 3]') + sum(worst_probs2 .* sorted_actions2[:, 4]'))
            elseif ocp.reward_type[1] == "worst-exp-lb-ub"
                # TOdo: Move all logic for this reward here
                total_reward = worst_exp_lb
            else
                @error "Unknown primary reward type!!!"
            end

            ###### Defining secondary reward just on DFA can result in weird behavior!!
            if ocp.reward_type[2] == "expected-dfa"
                x_new = [μ["x1"], μ["x2"]]
                x_new = reshape(x_new, 2, 1)
                secondary_reward = propogate_pimdp_trajectory_test(ocp.pimdp, ocp.dfa, ocp.extents, x_new)
            elseif ocp.reward_type[2] == "worst-dfa"
                secondary_reward = calculate_worst_dfa_reward(ocp.pimdp, ocp.dfa, imdp_extent_ids)
            elseif ocp.reward_type[2] == "best-dfa"
                secondary_reward = calculate_best_dfa_reward(ocp.pimdp, ocp.dfa, imdp_extent_ids)
            elseif ocp.reward_type[2] == "exp-probs"
                secondary_reward = (sum(minPrs) + sum(maxPrs))/length(maxPrs)
            elseif ocp.reward_type[2] == "combo"
                dfa_metric = calculate_worst_dfa_reward(ocp.pimdp, ocp.dfa, imdp_extent_ids)
                if dfa_metric > -Inf
                    secondary_reward = sum(minPrs) + sum(maxPrs)
                else
                    secondary_reward = dfa_metric
                end
            elseif ocp.reward_type[2] == "combo2"
                x_new = [μ["x1"], μ["x2"]]
                x_new = reshape(x_new, 2, 1)
                dfa_metric = propogate_pimdp_trajectory_test(ocp.pimdp, ocp.dfa, ocp.extents, x_new)
                if dfa_metric > -Inf
                    if isnothing(pimdp_reward_dict)
                        secondary_reward = calculate_pimdp_metric(ocp, mode)
                    else
                        secondary_reward = -pimdp_reward_dict[q_c]
                    end
                else
                    secondary_reward = dfa_metric
                end
            elseif ocp.reward_type[2] == "worst-exp-ub"
                sorted_actions = potential_actions[sortperm(potential_actions[:,4]), :]
                region_idx = sorted_actions[:,1]
                minPrs = sorted_actions[:,5]
                maxPrs = sorted_actions[:,6]
                worst_probs = get_true_transition_probabilities(minPrs, maxPrs, collect(1:1:length(region_idx)))
                secondary_reward = sum(worst_probs .* sorted_actions[:, 4]')
            elseif ocp.reward_type[2] == "worst-exp-ub-same"
                sorted_actions = potential_actions[sortperm(potential_actions[:,4]), :]
                secondary_reward = sum(worst_probs .* sorted_actions[:, 4]')
            elseif ocp.reward_type[2] == "none"
                secondary_reward = -Inf
            elseif ocp.reward_type[2] == "exp-pimdp-metric1"
                # x_new = [μ["x1"], μ["x2"]]
                # x_new = reshape(x_new, 2, 1)
                # pimdp_state_ex = propogate_pimdp_trajectory_test(ocp.pimdp, ocp.dfa, ocp.extents, x_new)

                secondary_reward = 0
                for ps in nt_pimdp_states
                    val = pimdp_reward_dict[ps]
                    if val < Inf
                        secondary_reward -= pimdp_reward_dict[ps]
                    end
                end
                secondary_reward *= 1/length(nt_pimdp_states)

                # if isnothing(pimdp_reward_dict)
                #     secondary_reward = calculate_pimdp_metric(ocp, mode)
                # else
                #     secondary_reward = -pimdp_reward_dict[pimdp_state_ex]
                # end
            elseif ocp.reward_type[2] == "worst-dfa-exp-pimdp-metric1"
                dfa_metric = calculate_worst_dfa_reward(ocp.pimdp, ocp.dfa, imdp_extent_ids)
                if dfa_metric > -Inf
                    # x_new = [μ["x1"], μ["x2"]]
                    # x_new = reshape(x_new, 2, 1)
                    # pimdp_state_ex = propogate_pimdp_trajectory_test(ocp.pimdp, ocp.dfa, ocp.extents, x_new)
                    secondary_reward = 0
                    for ps in nt_pimdp_states
                        val = pimdp_reward_dict[ps]
                        if val < Inf
                            secondary_reward -= pimdp_reward_dict[ps]
                        end
                    end
                    secondary_reward *= 1/length(nt_pimdp_states)
                    # mean([-pimdp_reward_dict[ps] for ps in nt_pimdp_states])
                    # if isnothing(pimdp_reward_dict)
                    #     secondary_reward = calculate_pimdp_metric(ocp, mode)
                    # else
                    #     secondary_reward = -pimdp_reward_dict[pimdp_state_ex]
                    # end
                else
                    secondary_reward = dfa_metric
                end
            else
                @error "Unknown sec reward type!!!"
            end

            if total_reward > best_reward[1]
                best_mode = mode
                best_reward[1] = total_reward
                best_reward[2] = secondary_reward
                best_gp_info_dict = gp_info_dict
                best_ub = worst_exp_ub
            elseif total_reward == best_reward[1]
                if worst_exp_ub > best_ub && ocp.reward_type[1] == "worst-exp-lb-ub"
                    best_mode = mode
                    best_ub = worst_exp_ub
                    best_reward[2] = secondary_reward
                    best_gp_info_dict = gp_info_dict
                elseif secondary_reward > best_reward[2]
                    best_reward[2] = secondary_reward
                    best_mode = mode
                    best_gp_info_dict = gp_info_dict
                elseif secondary_reward == best_reward[2]
                    # This chooses a random action among a final tie
                    if rand(1)[1] > 0.5
                        best_reward[2] = secondary_reward
                        best_mode = mode
                        best_gp_info_dict = gp_info_dict
                    end
                end
            end
        end
    end
    return best_mode, best_reward, best_gp_info_dict
end

"""
Get pimdp states with non-trivial transition bounds
"""
function get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_state, mode, foo)
    P_row_idx = map_pimdp_state_to_Pmat_row(pimdp_state, length(dfa.states), length(pimdp.actions), mode)
    # @info pimdp_state
    prow_max = pimdp.Pmax[P_row_idx, :]

    # Find all states with interval bounds with non-zero UB
    regions_with_nonzero_ub = findall(x->x>0., prow_max)
    return regions_with_nonzero_ub
end

"""
Get pimdp states with non-trivial transition bounds for all actions
"""
function get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_state, num_actions)
    pimdp_states = []
    for mode=1:num_actions
        pimdp_states = pimdp_states ∪ get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_state, mode, "foo")
    end
    # Get unique states
    unique!(pimdp_states)
    return pimdp_states
end

"""
PIMDP Distance Metric
"""
function calculate_pimdp_metric(pimdp, pimdp_state, dfa, num_actions; max_depth=Inf, helper_dict=nothing)
    # Get the current pimdp state
    # pimdp_state = pimdp.state_history[end]

    # Get the labels that will induce a transition
    target_labels = get_next_dfa_labels(dfa, pimdp_state[2])
    sink_labels = get_dfa_sink_labels(dfa)

    distance = 0
    # From the current PIMDP state, get the next possible states
    pimdp_states = get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_state, 1, "foo")

    pimdp_states_tried = []
    searching_flag = true

    best_dist = Inf
    while !isempty(pimdp_states) && searching_flag && distance < max_depth
        distance += 1
        next_states = []

        for state in pimdp_states
            imdp_state_idx = Int(pimdp_col_to_imdp_state(state, length(dfa.states)))
            pimdp_statep =  (imdp_state_idx, pimdp_state[2])
            if pimdp_statep ∈ keys(helper_dict)
                searching_flag = false
                subdistance = helper_dict[pimdp_statep] + distance
                if subdistance < best_dist
                    best_dist = subdistance
                end
            elseif pimdp.imdp_label_dict[imdp_state_idx] ∈ target_labels
                searching_flag = false
                if distance < best_dist
                    best_dist = distance
                end
            end

            push!(pimdp_states_tried, state)
            if searching_flag
                next_states = next_states ∪ get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_statep, num_actions)
            end
        end

        pimdp_states = setdiff(next_states, pimdp_states_tried)
        @assert isempty(pimdp_states ∩ pimdp_states_tried)
    end
        # For each next possible state, check the imdp label
    return best_dist
end

"""
Calculate Min Distance to Accepting state
"""
function calculate_pimdp_metric_all_offline(pimdp, dfa, num_actions; max_depth=Inf)
    # Get the current pimdp state
    pimdp_dist_dict = Dict()

    # Get the labels that will induce a transition
    states_finished = []
    sink_labels = get_dfa_sink_labels(dfa)
    for imdp_idx in keys(pimdp.imdp_label_dict), dfa_state=1:length(dfa.states)
    # for imdp_idx in [1], dfa_state=1:length(dfa.states)
        pimdp_state = (imdp_idx, dfa_state)
        imdp_label = pimdp.imdp_label_dict[imdp_idx]

        if pimdp_state[2] == dfa.sink_state || imdp_label ∈ sink_labels
            pimdp_dist_dict[pimdp_state] = Inf
        elseif pimdp_state[2] == dfa.accepting_state || possible_accept_transition(dfa, dfa_state, imdp_label)
            pimdp_dist_dict[pimdp_state] = 0
        elseif pimdp_state ∈ keys(pimdp_dist_dict)
            continue
        else
            val = calculate_pimdp_metric(pimdp, pimdp_state, dfa, num_actions, helper_dict=pimdp_dist_dict)
            pimdp_dist_dict[pimdp_state] = val
        end
        #     # From the current PIMDP state, get the next possible states
        #     pimdp_states = get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_state, num_actions)
        #     @info pimdp_states
        #     pimdp_states_tried = []
        #     searching_flag = true
        #     distance = 0
        #     while !isempty(pimdp_states) && searching_flag && distance < max_depth
        #         distance += 1
        #         @info distance
        #         next_states = []
        #         for state in pimdp_states

        #             imdp_idxp = pimdp_col_to_imdp_state(state, length(dfa.states))
        #             imdp_labelp = pimdp.imdp_label_dict[imdp_idxp]
        #             dfa_statep = state % length(dfa.states)
        #             pimdp_statep = (Int(imdp_idxp), dfa_statep)

        #             if (pimdp_statep[2] == 0 && imdp_labelp ∉ sink_labels) || possible_accept_transition(dfa, dfa_statep, imdp_labelp)
        #                 @info "STATE FOUND --------"
        #                 @info pimdp_statep
        #                 searching_flag = false
        #                 pimdp_dist_dict[pimdp_state] = distance
        #             elseif pimdp_statep in keys(pimdp_dist_dict)
        #                 searching_flag = false
        #                 pimdp_dist_dict[pimdp_state] = pimdp_dist_dict[pimdp_statep] + distance
        #             end
        #             push!(pimdp_states_tried, state)
        #             next_states = next_states ∪ get_nontrivial_next_pimdp_states(pimdp, dfa, pimdp_statep, num_actions)
        #             @info next_states
        #         end
        #         pimdp_states = setdiff(next_states, pimdp_states_tried)
        #         # Check for states already procced here
        #         @assert isempty(pimdp_states ∩ pimdp_states_tried)
        #     end
        # end
        @assert pimdp_state ∈ keys(pimdp_dist_dict)
    end
    return pimdp_dist_dict
end

function get_distance_to_dfa_label(ocp, mode)
    # TODO: Finish this
   # Get the current pimdp state
   pimdp_state = ocp.pimdp.state_history[end]

   # Get the labels that will induce a transition
   target_labels = get_next_dfa_labels(ocp.dfa, pimdp_state[2])

   # From the current PIMDP state, get the next possible states
   pimdp_states = get_nontrivial_next_pimdp_states(ocp, pimdp_state, mode)

   distance = Inf
   for state in pimdp_states
        imdp_state_idx = Int(pimdp_col_to_imdp_state(state, length(ocp.dfa.states)))

   end
end

"""
Get the policy from a synthesis result saved as a matrix.
"""
function get_policy_from_offline_synthesis_mat(filename)
    policy = Dict()
    synth_dict = matread(filename)

    policy = synth_dict["policy"]
    minPr = synth_dict["indVmin"]
    maxPr = synth_dict["indVmax"]

    for (i, action) in enumerate(synth_dict["policy"])
        policy[i] = Int(action)
    end
    return policy, minPr, maxPr
end


text:
export TfRanker

using PyCall: python

@with_kw mutable struct TfRanker <: BaseEstimator
    tf::Dict{String, Vector{UInt8}} = Dict()
    nstep::Int = 1000
    lr::Float32 = 1f-2
    pdrop::Float32 = 0
    hidden::String = "10"
    loss::String = "pairwise_logistic_loss"
end

is_classifier(::TfRanker) = false
is_ranker(::TfRanker) = true

function gridparams(m::TfRanker)
    grid = [
        "horizon" => [2, 5, 10, 20, 50],
        "nlevel" => [5, 10, 20],
        "nstep" => [1000, 2000, 5000],
        "lr" => [1f-3, 1f-2, 1f-1],
        "pdrop" => [0f0, 0.3f0, 0.5f0],
        "hidden" => ["10", "20", "50", "100"],
        "loss" =>
        ["pairwise_logistic_loss",
         "pairwise_hinge_loss",
         "pairwise_soft_zero_one_loss",
         "softmax_loss",
         "sigmoid_cross_entropy_loss",
         "mean_squared_loss",
         "list_mle_loss",
         "approx_ndcg_loss"]
    ]
    params = gridparams(grid)
end

function fit!(m::TfRanker, x, y, w = nothing; group = nothing, columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    if isnothing(group) && ndims(x) == 3
        group = repeat([size(x, 2)], size(x, 3))
    end
    ENV["nfea"] = size(x, 1)
    ENV["nlist"] = min(500, minimum(group))
    dst = to_svm(x, y, w, group = group)
    out = mktempdir()
    run(rankcmd(m, dst, out))
    m.tf = readtf(out)
    foreach(rmdir, [dst, out])
    BSON.bson("model.bson", model = m)
    return m
end

function predict(m::TfRanker, x)
    @unpack tf = m
    ENV["nfea"] = size(x, 1)
    ENV["nlist"] = prod(size(x)[2:end])
    dst = to_svm(x)
    out = mktempdir()
    writetf(out, tf)
    run(rankcmd(m, dst, out, predonly = true))
    ŷ = readdlm("preds", Float32)
    foreach(rmdir, ["preds", dst, out])
    return ŷ
end

function rankcmd(m::TfRanker, dst, out; predonly = false)
    @unpack lr, pdrop, hidden, loss, nstep = m
    @unpack nfea, nlist = ENV
    nstep = predonly ? 0 : nstep
    if !usegpu()
        ENV["CUDA_VISIBLE_DEVICES"] = -1
        ENV["MKL_NUM_THREADS"] = 20
        ENV["OPENBLAS_NUM_THREADS"] = 20
        ENV["OMP_NUM_THREADS"] = 20
    end
    tfrank = joinpath(@__DIR__, "tfrank.py")
    `$python $tfrank
    --train_path=$dst
    --vali_path=$dst
    --test_path=$dst
    --output_dir=$out
    --train_batch_size=32
    --num_features=$nfea
    --list_size=$nlist
    --num_train_steps=$nstep
    --learning_rate=$lr
    --dropout_rate=$pdrop
    --hidden_layer_dims=$hidden
    --loss=$loss`
end

readtf(dir) = Dict(src => read(joinpath(dir, src)) for src in readdir(dir))

function writetf(dir, dict)
    mkpath(dir)
    for (src, bytes) in dict
        write(joinpath(dir, src), bytes)
    end
end

ext: jl


text:
__precompile__()

module BioTools

include("blast/BLAST.jl")

end # module Tools


text:
"""
    initdb(name)

Initialize database in file `name`.
Accept only empty database or database with SQLiteFs specific tables.
All not-existing tables and other DB objects are created.
Returns the db-session handle needed in all database statements.
"""
function initdb(name::AbstractString)
    db = SQLite.DB(name)
    for command in CREATE_SQL
        DBInterface.execute(db, command)
    end
    SQLite.register(db, modestring)
    SQLite.register(db, updatemode)

    st = FStatus(db) # initialize home directory name and inode
    create_rootnode(st)
    st
end

# const TODO
CREATE_SQL = ["""
    PRAGMA FOREIGN_KEYS = TRUE;
    """,
    """
    CREATE TABLE IF NOT EXISTS inode(
        ino INTEGER PRIMARY KEY,
        nlinks INTEGER DEFAULT 0,
        mode INTEGER,
        size INTEGER,
        atime INTEGER,
        btime INTEGER,
        ctime INTEGER,
        mtime INTEGER
    );
    """,
    """
    CREATE TABLE IF NOT EXISTS data(
        ino INTEGER,
        data BLOB,
        FOREIGN KEY (ino) REFERENCES inode ON DELETE CASCADE
    );
    """,
    """
    CREATE TABLE IF NOT EXISTS direntry(
        dino INTEGER NOT NULL,
        name TEXT NOT NULL,
        ino INTEGER NOT NULL,
        PRIMARY KEY (dino, name ASC),
        FOREIGN KEY (ino) REFERENCES inode ON DELETE CASCADE,
        FOREIGN KEY (dino) REFERENCES inode (ino) ON DELETE CASCADE
    );
    """,
    """
    CREATE TRIGGER IF NOT EXISTS nlinks_ins INSERT ON direntry
        BEGIN
            UPDATE inode SET nlinks = nlinks + 1 WHERE ino = new.ino;
        END
    ;
    """,
    """
    CREATE TRIGGER IF NOT EXISTS nlinks_del DELETE ON direntry
        BEGIN
            UPDATE inode SET nlinks = nlinks - 1 WHERE ino = old.ino;
        END
    ;
    """,
    """
    CREATE TRIGGER IF NOT EXISTS nlinks_up UPDATE of ino ON direntry
        BEGIN
            UPDATE inode SET nlinks = nlinks - 1 WHERE ino = old.ino;
            UPDATE inode SET nlinks = nlinks + 1 WHERE ino = new.ino;
        END
    ;
    """,
    """
    CREATE TRIGGER IF NOT EXISTS nlinks_zero UPDATE of nlinks ON inode
        WHEN new.nlinks = 0
        BEGIN
            DELETE FROM inode WHERE ino = old.ino;
        END
    ;
    """,
    """
    CREATE TRIGGER IF NOT EXISTS direntry_del DELETE ON direntry
        WHEN old.ino IN (SELECT dino from direntry)
        BEGIN
            SELECT RAISE(FAIL, 'directory is not empty');
        END
    ;
    """,
    """
    CREATE VIEW IF NOT EXISTS pathes AS
    WITH RECURSIVE pathnames (path, ino) AS
    (
        SELECT name, ino from direntry where dino <= 1
            UNION ALL
        SELECT path || '/' || name, direntry.ino FROM pathnames, direntry
            WHERE pathnames.ino = dino AND dino > 1
    )
    SELECT * FROM pathnames;
    """,
    """
    CREATE VIEW IF NOT EXISTS ls_data AS
    SELECT mode, mtime, dino, inode.ino, nlinks, name FROM inode, direntry
    WHERE direntry.ino = inode.ino
    ORDER BY dino, name;
    """,
    """
    CREATE VIEW IF NOT EXISTS ls_data_r AS
    WITH RECURSIVE under_root(ino, name, level) AS (
        SELECT ino, name, 0 FROM direntry where ino = 1
        UNION ALL
        SELECT direntry.ino, direntry.name, under_root.level + 1
        FROM direntry, under_root
        WHERE dino = under_root.ino AND dino != direntry.ino
        ORDER BY 3 DESC
        )
    SELECT printf('%*s', level, '') || name
    FROM under_root;
    """,
    ]


text:
@testset "Composite types" begin
    c = CompositeObject((1,2,3,4))
    @test iscomposite(c)
    @test components(c) == (1,2,3,4)
    @test ncomponents(c) == length(components(c))
    @test component(c, 2) == 2
    @test component(c, 4) == 4
    @test_throws BoundsError component(c, 5)

    @test !iscomposite(5)
    @test components(5) == ()
    @test ncomponents(5) == 0
    @test_throws BoundsError component(5, 1)

    c2 = CompositeObject([1,2,3,4])
    @test setcomponent!(c2, 4, 2) == [1,4,3,4]
    @test component(c2, 2) == 4
end

Base.getindex(x::CompositeObject, I::ComponentIndex...) =
    component(x, map(Indexing.to_index, I)...)
Base.setindex!(x::CompositeObject, val, I::ComponentIndex...) =
    setcomponent!(x, val, map(Indexing.to_index, I)...)

@testset "Indexing" begin
    c1 = CompositeObject([1, 2, 3, 4])
    @test c1[Component(1)] == 1
    @test c1[Component(3)] == 3
    @test c1[Component(1):Component(3)] == 1:3
    @test_throws BoundsError c1[Component(5)]

    A = rand(10,10)
    c2 = CompositeObject(A)
    @test c2[Component(1),Component(2)] == A[1,2]
    @test c2[Component(1):Component(4),Component(2)] == A[1:4,2]
    @test c2[Component(:),Component(2)] == A[:,2]
    @test c2[Component(:)] == A[:]
    @test c2[Component(CartesianIndex(2,3))] == A[2,3]
    @test c2[Component(CartesianIndices((2,2)))] == A[CartesianIndices((2,2))]
    c2[Component(1),Component(2)] = 2
    @test component(c2, 1, 2) == 2
end


text:
#####
##### Profiling backpropagation
#####

using AlphaZero
using Setfield
import CUDA

# nsys launch julia --project
# using Revise; Revise.includet("scripts/profile/backprop.jl")
# profile_backprop()
# nsight-sys ...

# Takeaway: time spent on GC skyrockets when you approach the memory limit

function profile_backprop(
    exp::Experiment = Examples.experiments["connect-four"];
    profile=false,
    num_games=ConstSchedule(5000),
    num_batches=100,
    batch_size=ConstSchedule(1024),
    num_filters=64)

  exp = @set exp.params.learning.batch_size = batch_size
  exp = @set exp.netparams.num_filters = num_filters
  session = Session(exp, autosave=false, dir="sessions/profile-backprop-$(exp.name)")
  env = session.env
  for i in 1:num_games[session.env.itc]
    trace = play_game(exp.gspec, RandomPlayer())
    AlphaZero.push_trace!(session.env.memory, trace, 1.0)
  end
  nparams = Network.num_parameters(env.curnn)
  UI.Log.section(session.logger, 1, "Profiling backprop")
  experience = AlphaZero.get_experience(env.memory)
  tr = AlphaZero.Trainer(env.gspec, env.curnn, experience, env.params.learning)
  AlphaZero.batch_updates!(tr, 1)  # to compile everything
  if profile
    CUDA.@profile AlphaZero.batch_updates!(tr, num_batches)
  else
    CUDA.@time AlphaZero.batch_updates!(tr, num_batches)
  end
  return
end

profile_backprop()


text:
import MeasureTheory.Kernel, MeasureTheory.kernel

MeasureTheory.mapcall(t) = map(func -> func(), t)
(k::Kernel{M,<:NamedTuple})() where {M} = M(; MeasureTheory.mapcall(k.ops)...)

const GaussKernel = MeasureTheory.Kernel{<:Gaussian}
#const Copy = MeasureTheory.Kernel{<:Dirac}

"""
    AffineMap(B, β)

Represents a function `f = AffineMap(B, β)`
such that `f(x) == B*x + β`.
"""
struct AffineMap{S,T}
    B::S
    β::T
end
(a::AffineMap)(x) = a.B*x + a.β
(a::AffineMap)(p::Gaussian) = Gaussian(μ = a.B*mean(p) + a.β, Σ = a.B*cov(p)*a.B')
(a::AffineMap)(p::WGaussian) = WGaussian(μ = a.B*mean(p) + a.β, Σ = a.B*cov(p)*a.B', c=p.c)

"""
    LinearMap(B)

Represents a function `f = LinearMap(B)`
such that `f(x) == B*x`.
"""
struct LinearMap{T}
    B::T
end
(a::LinearMap)(x) = a.B*x
(a::LinearMap)(p::Gaussian) = Gaussian(μ = a.B*mean(p), Σ = a.B*cov(p)*a.B')
(a::LinearMap)(p::WGaussian) = WGaussian(μ = a.B*mean(p), Σ = a.B*cov(p)*a.B', c=p.c)

"""
    ConstantMap(β)

Represents a function `f = ConstantMap(β)`
such that `f(x) == β`.
"""
struct ConstantMap{T}
    x::T
end
(a::ConstantMap)(x) = a.x
(a::ConstantMap)() = a.x

const LinearGaussianKernel = Kernel{T,NamedTuple{(:μ, :Σ),Tuple{A, C}}} where {T<:Gaussian, A<:LinearMap, C<:ConstantMap}
params(k::LinearGaussianKernel) = k.ops.μ.B, Zero(), k.ops.Σ.x
const AffineGaussianKernel = Kernel{T,NamedTuple{(:μ, :Σ),Tuple{A, C}}} where {T<:Gaussian, A<:AffineMap, C<:ConstantMap}
params(k::AffineGaussianKernel) = k.ops.μ.B, k.ops.μ.β, k.ops.Σ.x
const ConstantGaussianKernel =  Kernel{T,NamedTuple{(:μ, :Σ),Tuple{A, C}}} where {T<:Gaussian, A<:ConstantMap, C<:ConstantMap}
params(k::ConstantGaussianKernel) = k.ops.μ.x, k.ops.Σ.x

"""
    correct(prior, obskernel, obs) = u, yres, S

Joseph form correction step of a Kalman filter with `prior` state
and `obs` the observation with observation kernel
`obskernel = kernel(Gaussian; μ=LinearMap(H), Σ=ConstantMap(R))`
`H` is the observation operator and `R` the observation covariance. Returns corrected/conditional
distribution `u`, the residual and the innovation covariance.
See https://en.wikipedia.org/wiki/Kalman_filter#Update.
"""
function correct(u::Gaussian{T}, k::LinearGaussianKernel, y) where {T}
    x, Ppred = meancov(u)
    H = k.ops.μ.B
    R = k.ops.Σ.x
    (x, P), yres, S = correct_joseph((x, Ppred), H, R, y)
    Gaussian{T}(μ=x, Σ=P), yres, S
end

function correct_joseph((x, Ppred), H, R, y)
    yres = y - H*x # innovation residual
    S = (H*Ppred*H' + R) # innovation covariance

    K = Ppred*H'/S # Kalman gain
    x = x + K*yres
    P = (I - K*H)*Ppred*(I - K*H)' + K*R*K' #  Joseph form
    (x, P), yres, S
end

function correct_kalman((x, Ppred), H, R, y)
    yres = y - H*x # innovation residual
    S = (H*Ppred*H' + R) # innovation covariance

    K = Ppred*H'/S # Kalman gain
    x = x + K*yres
    P = (I - K*H)*Ppred # Kalman form
    (x, P), yres, S
end

function smooth(Gs::T, Gf, Gpred, Φ) where {T}
    xs, Ps = meancov(Gs)
    xf, Pf = meancov(Gf)
    xpred, Ppred = meancov(Gpred)

    J = Pf*Φ'/Ppred # C/(C+w)
    xs = xf + J*(xs - xpred)
    Ps = Pf + J*(Ps - Ppred)*J'

    (xs, Ps)
end


text:
@broadcastref struct Block
    id::ID
    insts::Vector{Instruction}
end

struct ControlFlowGraph
    blocks::Vector{Block}
    graph::SimpleDiGraph
end

ControlFlowGraph() = ControlFlowGraph([], SimpleDiGraph())

function get_vertex!(cfg::ControlFlowGraph, block::Block)
    i = findfirst(==(block), cfg.blocks)
    if isnothing(i)
        push!(cfg.blocks, block)
        add_vertex!(cfg.graph)
        nv(cfg.graph)
    else
        i
    end
end

Graphs.add_edge!(cfg::ControlFlowGraph, src::Block, dst::Block) = add_edge!(cfg.graph, get_vertex!(cfg, src), get_vertex!(cfg, dst))

struct FunctionDefinition
    type::ID
    control::FunctionControl
    args::Vector{ID}
    cfg::ControlFlowGraph
end

body(fdef::FunctionDefinition) = foldl(append!, map(x -> x.insts, fdef.cfg.blocks); init=Instruction[])


text:
module MultiPartParsing

import ..access_threaded, ..Request, ..Multipart, ..payload
using ..Parsers

export parse_multipart_form

const CR_BYTE = 0x0d # \r
const LF_BYTE = 0x0a # \n
const DASH_BYTE = 0x2d # -
const HTAB_BYTE = 0x09 # \t
const SPACE_BYTE = 0x20
const SEMICOLON_BYTE = UInt8(';')
const CRLFCRLF = (CR_BYTE, LF_BYTE, CR_BYTE, LF_BYTE)

"compare byte buffer `a` from index `i` to index `j` with `b` and check if they are byte-equal"
function byte_buffers_eq(a, i, j, b)
    l = 1
    @inbounds for k = i:j
        a[k] == b[l] || return false
        l += 1
    end
    return true
end

"""
    find_multipart_boundary(bytes, boundaryDelimiter; start::Int=1)

Find the first and last index of the next boundary delimiting a part, and if
the discovered boundary is the terminating boundary.
"""
function find_multipart_boundary(bytes::AbstractVector{UInt8}, boundaryDelimiter::AbstractVector{UInt8}; start::Int=1)
    # The boundary delimiter line is prepended with two '-' characters
    # The boundary delimiter line starts on a new line, so must be preceded by a \r\n.
    # The boundary delimiter line ends with \r\n, and can have "optional linear whitespace" between
    # the end of the boundary delimiter, and the \r\n.
    # The last boundary delimiter line has an additional '--' at the end of the boundary delimiter
    # [RFC2046 5.1.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)

    i = start
    end_index = i + length(boundaryDelimiter) + 1
    while end_index <= length(bytes)
        if bytes[i] == DASH_BYTE && bytes[i + 1] == DASH_BYTE && byte_buffers_eq(bytes, i + 2, end_index, boundaryDelimiter)
            # boundary delimiter line start on a new line ...
            if i > 1
                (i == 2 || bytes[i-2] != CR_BYTE || bytes[i-1] != LF_BYTE) && error("boundary delimiter found, but it was not the start of a line")
                # the CRLF preceding the boundary delimiter is "conceptually attached
                # to the boundary", so account for this with the index
                i -= 2
            end

            # need to check if there are enough characters for the CRLF or for two dashes
            end_index < length(bytes)-1 || error("boundary delimiter found, but did not end with new line")

            is_terminating_delimiter = bytes[end_index+1] == DASH_BYTE && bytes[end_index+2] == DASH_BYTE
            is_terminating_delimiter && (end_index += 2)

            # ... there can be arbitrary SP and HTAB space between the boundary delimiter ...
            while end_index < length(bytes) && (bytes[end_index+1] in (HTAB_BYTE, SPACE_BYTE))
                end_index += 1
            end
            # ... and ends with a new line
            newlineEnd = end_index < length(bytes)-1 &&
                         bytes[end_index+1] == CR_BYTE &&
                         bytes[end_index+2] == LF_BYTE
            if !newlineEnd
                error("boundary delimiter found, but did not end with new line")
            end

            end_index += 2

            return (is_terminating_delimiter, i, end_index)
        end

        i += 1
        end_index += 1
    end

    error("boundary delimiter not found")
end

"""
    find_multipart_boundaries(bytes, boundary; start=1)

Find the start and end indexes of all the parts of the multipart object.  Ultimately this method is
looking for the data between the boundary delimiters in the byte array.  A vector containing all
the start/end pairs is returned.
"""
function find_multipart_boundaries(bytes::AbstractVector{UInt8}, boundary::AbstractVector{UInt8}; start=1)
    idxs = Tuple{Int, Int}[]
    while true
        (is_terminating_delimiter, i, end_index) = find_multipart_boundary(bytes, boundary; start = start)
        push!(idxs, (i, end_index))
        is_terminating_delimiter && break
        start = end_index + 1
    end
    return idxs
end

"""
    find_header_boundary(bytes)

Find the end of the multipart header in the byte array.  Returns a Tuple with the
start index(1) and the end index. Headers are separated from the body by CRLFCRLF.

[RFC2046 5.1](https://tools.ietf.org/html/rfc2046#section-5.1)
[RFC822 3.1](https://tools.ietf.org/html/rfc822#section-3.1)
"""
function find_header_boundary(bytes::AbstractVector{UInt8})
    length(CRLFCRLF) > length(bytes) && return nothing

    l = length(bytes) - length(CRLFCRLF) + 1
    i = 1
    end_index = length(CRLFCRLF)
    while (i <= l)
        byte_buffers_eq(bytes, i, end_index, CRLFCRLF) && return (1, end_index)
        i += 1
        end_index += 1
    end
    error("no delimiter found separating header from multipart body")
end

const content_disposition_regex = Parsers.RegexAndMatchData[]
function content_disposition_regex_f()
    r = Parsers.RegexAndMatchData(r"^Content-Disposition:[ \t]*form-data;[ \t]*(.*)\r\n"x)
    Parsers.init!(r)
end

const content_disposition_flag_regex = Parsers.RegexAndMatchData[]
function content_disposition_flag_regex_f()
    r = Parsers.RegexAndMatchData(r"""^
    [ \t]*([!#$%&'*+\-.^_`|~[:alnum:]]+);?
    """x)
    Parsers.init!(r)
end

const content_disposition_pair_regex = Parsers.RegexAndMatchData[]
function content_disposition_pair_regex_f()
    r = Parsers.RegexAndMatchData(r"""^
    [ \t]*([!#$%&'*+\-.^_`|~[:alnum:]]+)[ \t]*=[ \t]*"(.*?)";?
    """x)
    Parsers.init!(r)
end

const content_type_regex = Parsers.RegexAndMatchData[]
function content_type_regex_f()
    r = Parsers.RegexAndMatchData(r"(?i)Content-Type: (\S*[^;\s])"x)
    Parsers.init!(r)
end

"""
    parse_multipart_chunk(chunk)

Parse a single multi-part chunk into a Multipart object.  This will decode
the header and extract the contents from the byte array.
"""
function parse_multipart_chunk(chunk)
    startIndex, end_index = find_header_boundary(chunk)
    header = SubString(unsafe_string(pointer(chunk, startIndex), end_index - startIndex + 1))
    content = view(chunk, end_index+1:lastindex(chunk))

    # find content disposition
    re = access_threaded(content_disposition_regex_f, content_disposition_regex)
    if !Parsers.exec(re, header)
        @warn "Content disposition is not specified dropping the chunk." String(chunk)
        return nothing # Specifying content disposition is mandatory
    end
    content_disposition = Parsers.group(1, re, header)

    re_flag = access_threaded(content_disposition_flag_regex_f, content_disposition_flag_regex)
    re_pair = access_threaded(content_disposition_pair_regex_f, content_disposition_pair_regex)
    name = nothing
    filename = nothing
    while !isempty(content_disposition)
        if Parsers.exec(re_pair, content_disposition)
            key = Parsers.group(1, re_pair, content_disposition)
            value = Parsers.group(2, re_pair, content_disposition)
            if key == "name"
                name = value
            elseif key == "filename"
                filename = value
            else
                # do stuff with other content disposition key-value pairs
            end
            content_disposition = Parsers.nextbytes(re_pair, content_disposition)
        elseif Parsers.exec(re_flag, content_disposition)
            # do stuff with content disposition flags
            content_disposition = Parsers.nextbytes(re_flag, content_disposition)
        else
            break
        end
    end

    name === nothing && return

    re_ct = access_threaded(content_type_regex_f, content_type_regex)
    contenttype = Parsers.exec(re_ct, header) ? Parsers.group(1, re_ct, header) : "text/plain"

    return Multipart(filename, IOBuffer(content), contenttype, "", name)
end

"""
    parse_multipart_body(body, boundary)::Vector{Multipart}

Parse the multipart body received from the client breaking it into the various
chunks which are returned as an array of Multipart objects.
"""
function parse_multipart_body(body::AbstractVector{UInt8}, boundary::AbstractString)::Vector{Multipart}
    multiparts = Multipart[]
    idxs = find_multipart_boundaries(body, codeunits(boundary))
    length(idxs) > 1 || (return multiparts)

    for i in 1:length(idxs)-1
        chunk = view(body, idxs[i][2]+1:idxs[i+1][1]-1)
        push!(multiparts, parse_multipart_chunk(chunk))
    end
    return multiparts
end

"""
    parse_multipart_form(req::Request)::Vector{Multipart}

Parse the full mutipart form submission from the client returning and
array of Multipart objects containing all the data.

The order of the multipart form data in the request should be preserved.
[RFC7578 5.2](https://tools.ietf.org/html/rfc7578#section-5.2).

The boundary delimiter MUST NOT appear inside any of the encapsulated parts. Note
that the boundary delimiter does not need to have '-' characters, but a line using
the boundary delimiter will start with '--' and end in \r\n.
[RFC2046 5.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)
"""
function parse_multipart_form(req::Request)::Union{Vector{Multipart}, Nothing}
    # parse boundary from Content-Type
    m = match(r"multipart/form-data; boundary=(.*)$", req["Content-Type"])
    m === nothing && return nothing

    boundary_delimiter = m[1]

    # [RFC2046 5.1.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)
    length(boundary_delimiter) > 70 && error("boundary delimiter must not be greater than 70 characters")

    return parse_multipart_body(payload(req), boundary_delimiter)
end

function __init__()
    resize!(empty!(content_disposition_regex), Threads.nthreads())
    resize!(empty!(content_disposition_flag_regex), Threads.nthreads())
    resize!(empty!(content_disposition_pair_regex), Threads.nthreads())
    resize!(empty!(content_type_regex), Threads.nthreads())
    return
end

end # module MultiPartParsing


text:
using Memoization
using Test

##

@testset "Memoization" begin

global n = 0

@eval @memoize f(x::T, y=nothing, args... ; z::Union{Int,Nothing}, w=nothing, kwargs...) where {T<:Int} = (global n+=1; ((x,y,args...), (;z=z,w=w,kwargs...)))

# args
n=0; @test (f(1,    z=1) == ((1,nothing),(z=1,w=nothing)) && n==1)
n=0; @test (f(1,    z=1) == ((1,nothing),(z=1,w=nothing)) && n==0)

n=0; @test (f(1,1,  z=1) == ((1,1),      (z=1,w=nothing)) && n==1)
n=0; @test (f(1,1,  z=1) == ((1,1),      (z=1,w=nothing)) && n==0)

n=0; @test (f(1,1,1,z=1) == ((1,1,1),    (z=1,w=nothing)) && n==1)
n=0; @test (f(1,1,1,z=1) == ((1,1,1),    (z=1,w=nothing)) && n==0)

# kwargs
n=0; @test (f(1,z=1,w=1)     == ((1,nothing),(z=1,w=1))     && n==1)
n=0; @test (f(1,z=1,w=1)     == ((1,nothing),(z=1,w=1))     && n==0)

n=0; @test (f(1,z=1,w=1,p=1) == ((1,nothing),(z=1,w=1,p=1)) && n==1)
n=0; @test (f(1,z=1,w=1,p=1) == ((1,nothing),(z=1,w=1,p=1)) && n==0)

# cache clear
Memoization.empty_cache!(f)
n=0; @test (f(1,z=1) == ((1,nothing),(z=1,w=nothing)) && n==1)
Memoization.empty_cache!(f)
n=0; @test (f(1,z=1) == ((1,nothing),(z=1,w=nothing)) && n==1)

# Dict vs. IdDict cache vs. custom cache
@eval @memoize IdDict g1(x) = (global n+=1; x)
n=0; @test (g1([1,2])==[1,2] && n==1)
n=0; @test (g1([1,2])==[1,2] && n==1)
@eval @memoize   Dict g2(x) = (global n+=1; x)
n=0; @test (g2([1,2])==[1,2] && n==1)
n=0; @test (g2([1,2])==[1,2] && n==0)
@eval @memoize Dict() g3(x) = (global n+=1; x)
n=0; @test (g3([1,2])==[1,2] && n==1)
n=0; @test (g3([1,2])==[1,2] && n==0)

# redefinition
# should only clear cache in toplevel, hence @eval
@eval @memoize h(x) = x
@test h(2)==2 && h(2)==2
@eval @memoize h(x) = 2x
@test h(2)==4 && h(2)==4

# inference
@eval @memoize foo(x) = x
@eval @memoize foo(;x) = x
@test @inferred(foo(2)) == 2
@test @inferred(foo(x=2)) == 2

# closures
function make_func(x)
    @memoize func(y) = (global n+=1; (x,y))
end
k = make_func(1)
n=0; @test (k(2)==(1,2) && n==1)
n=0; @test (k(2)==(1,2) && n==0)
k′ = make_func(2)
n=0; @test (k′(2)==(2,2) && n==1)
n=0; @test (k′(2)==(2,2) && n==0)
n=0; @test (k(2)==(1,2) && n==0)

# https://github.com/JuliaLang/julia/issues/40576
let
    @memoize foo(x) = (n+=1; x)
    global n=0; @test foo(1)==1 && n==1
    global n=1; @test foo(1)==1 && n==1
end

# callables
@eval struct Baz{T} end
@eval @memoize (::Baz{T})(x::X) where {X, T<:Int} = (global n+=1; Int)
@eval @memoize (::Baz{T})(x::X) where {X, T<:String} = (global n+=1; String)
bazint = Baz{Int}()
bazstring = Baz{String}()
n=0; @test (bazint(2)==Int && n==1)
n=0; @test (bazint(2)==Int && n==0)
n=0; @test (bazstring(2)==String && n==1)
n=0; @test (bazstring(2)==String && n==0)
Memoization.empty_cache!(Baz{Int})
n=0; @test (bazint(2)==Int && n==1)
Memoization.empty_cache!(Baz{String})
n=0; @test (bazstring(2)==String && n==1)
Memoization.empty_cache!(Baz)
n=0; @test (bazint(2)==Int && n==1)
n=0; @test (bazstring(2)==String && n==1)

# unnamed args
@eval @memoize uarg(::Type{T}) where {T} = (global n+=1; T)
n=0; @test (uarg(Int)==Int && n==1)
n=0; @test (uarg(Int)==Int && n==0)
n=0; @test (uarg(Float64)==Float64 && n==1)
n=0; @test (uarg(Float64)==Float64 && n==0)

# redefining cache type
@eval @memoize IdDict redef_cache_toplevel(x) = x
@test_throws Exception @eval @memoize Dict redef_cache_toplevel(x) = x

# works inside generated functions
@memoize gen1(x) = x
@generated gen2() = gen1(1)
@test gen2() == 1

# works when called during precompilation of a package
@test read(`$(Base.julia_cmd()[1]) --project=$(joinpath(@__DIR__,"TestPackage")) -e "using TestPackage; print(TestPackage.foo(2))"`,String) == "2"

end


text:
module MOITests

using BARON
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const optimizer = BARON.Optimizer(PrLevel=0)
const caching_optimizer = MOIU.CachingOptimizer(
    MOIU.UniversalFallback(MOIU.Model{Float64}()), BARON.Optimizer(PrLevel=0));

function test_runtests()
    model = caching_optimizer#MOI.instantiate(BARON.Optimizer, with_bridge_type = Float64)
    MOI.set(model, MOI.RawOptimizerAttribute("PrLevel"), 0)
    # MOI.set(model, MOI.Silent(), true) # todo
    MOI.Test.runtests(model,
        MOI.Test.Config(
            atol = 1e-3,
            rtol = 1e-3,
            # optimal_status = MOI.LOCALLY_SOLVED,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.DualObjectiveValue,
                MOI.ObjectiveBound,
                MOI.DualStatus,
                MOI.ConstraintDual,
            ],
        );
        exclude = [
            "test_attribute_SolverVersion",      # unavailable
            "test_nonlinear_hs071_NLPBlockDual", # MathOptInterface.NLPBlockDual(1)
            "test_nonlinear_invalid",            # see below
            # returns NaN in expression and solver has to responde with:
            # MOI.get(model, MOI.TerminationStatus()) == MOI.INVALID_MODEL
            # this code will error when NaN is found (better than waiting to knoe about bad stuff)
            "test_variable_solve_ZeroOne_with_upper_bound",# fail is upstream
            "test_objective_ObjectiveFunction_blank", # fail is upstream
            "test_objective_FEASIBILITY_SENSE_clears_objective", # fail is upstream
            "test_linear_integer_solve_twice", # simply fails in the first solve
            # objective fails
            # BARON will set the same large number
            # for both abj and variables in case of unbounded
            "test_unbounded_MIN_SENSE_offset",
            "test_unbounded_MIN_SENSE",
            "test_unbounded_MAX_SENSE_offset",
            "test_unbounded_MAX_SENSE",
        ]
    )
    return
end

@testset "MOI Unit" begin
    test_runtests()
end

end # module


text:
using Plots
#using LsqFit
import Distributions as di
using Random
using Turing
using Statistics
using StatsPlots
using FillArrays

n = di.Normal(1,2)
sample_n = rand(n,20)

@model function normal_fit(data)
	μ ~ Uniform(-10,10)
	σ ~ Uniform(0,20)
    data ~ MvNormal(Fill(μ,length(data)),σ)
end

model1 = normal_fit(sample_n)
chain = Turing.sample(model1,NUTS(0.65),1000)
plot(chain)

rwd = Normal(0,1)
rw_steps = rand(rwd,999)

x = [0.0]
for step in rw_steps
    push!(x,x[end]+step)
end
plot(x)

steps = []
for i in 1:length(x)-1
    push!(steps,x[i+1]-x[i])
end

rw_steps
steps

@model function diffusion_fit(data)
	μ = 0
	σ ~ Uniform(0,20)
	for i in 1:length(data)
		data[i] ~ Normal(μ,σ)
	end
end

model2 = diffusion_fit(steps)
chain2 = Turing.sample(model2,NUTS(0.65),1000)


text:
"""

The `LibSerialPort` module provides access to [serial
ports](https://en.wikipedia.org/wiki/Serial_port)
([UARTs](https://en.wikipedia.org/wiki/Universal_asynchronous_receiver-transmitter)
or USB/Bluetooth devices that emulate their operating-system
interface) using the portable C library
[libserialport](http://sigrok.org/wiki/Libserialport).

It defines the `SerialPort` type, which is returned by
`LibSerialPort.open`. This is a subtype of `Base.IO` and can
therefore be used like a file handle, using many of the same `read`,
`write`, `print`, etc. methods defined for `Base.IO`.

# Example
```
using LibSerialPort

list_ports()
ports = get_port_list()

sp = LibSerialPort.open("/dev/ttyUSB0", 115200)
set_flow_control(sp)
sp_flush(sp, SP_BUF_BOTH)
write(sp, "hello\\n")
println(readline(sp))
close(sp)
```
"""
module LibSerialPort

include("wrap.jl")
using .Lib   # low-level wrapper for the libserialport C API
import .Lib: sp_flush, sp_drain, sp_output_waiting

export
    SerialPort,
    list_ports,
    get_port_list,
    get_port_settings,
    print_port_metadata,
    print_port_settings,
    set_speed,
    set_frame,
    set_flow_control,
    set_read_timeout,
    set_write_timeout,
    clear_read_timeout,
    clear_write_timeout,
    seteof,
    reseteof,
    nonblocking_read,

    # Re-export the following symbol definitions from LibSerialPort.Lib.
    # (We may want to phase out some of these later.)

    # Low-level functions currently extended by the high-level API
    sp_flush,
    sp_drain,
    sp_output_waiting,

    # Enum types
    SPMode,
    SPBuffer,
    SPParity,
    SPrts,
    SPcts,
    SPdtr,
    SPdsr,
    SPXonXoff,

    # Enum values
    SP_MODE_READ,
    SP_MODE_WRITE,
    SP_MODE_READ_WRITE,
    SP_BUF_INPUT,
    SP_BUF_OUTPUT,
    SP_BUF_BOTH,
    SP_PARITY_NONE,
    SP_PARITY_ODD,
    SP_PARITY_EVEN,
    SP_PARITY_MARK,
    SP_PARITY_SPACE,
    SP_RTS_OFF,
    SP_RTS_ON,
    SP_RTS_FLOW_CONTROL,
    SP_CTS_IGNORE,
    SP_CTS_FLOW_CONTROL,
    SP_DTR_OFF,
    SP_DTR_ON,
    SP_DTR_FLOW_CONTROL,
    SP_DSR_IGNORE,
    SP_DSR_FLOW_CONTROL,
    SP_XONXOFF_DISABLED,
    SP_XONXOFF_IN,
    SP_XONXOFF_OUT,
    SP_XONXOFF_INOUT

import Base: isopen, open, close, write, unsafe_write, flush,
    read, unsafe_read, bytesavailable, eof

mutable struct SerialPort <: IO
    ref::Port
    is_eof::Bool
    is_open::Bool

    # cumulative timeout limits, passed on to blocking libserialport functions
    # 0 means wait indefinitely, as per sigrok libserialport interface
    read_timeout_ms::Cuint
    write_timeout_ms::Cuint

    function SerialPort(ref, is_eof, is_open)
        sp = new(ref, is_eof, is_open, 0, 0)
        finalizer(destroy!, sp)
        return sp
    end
end

"""
A `Timeout` exception is thrown by blocking read or write functions
when the cumulative blocking time since that last call to
`set_read_timeout(t)` or `set_write_timeout(t)` has exceeded `t`
seconds.
"""
struct Timeout <: Exception
end

"""
    SerialPort(portname::AbstractString)

Constructs and returns a `SerialPort` object.
"""
SerialPort(portname::AbstractString) = SerialPort(sp_get_port_by_name(portname), false, false)

"""
    destroy!(sp::SerialPort)

Close the `SerialPort` object and deallocate associated memory
objects.
"""
function destroy!(sp::SerialPort)
    close(sp)
    sp_free_port(sp.ref)
end

"""
    isopen(sp::SerialPort) -> Bool

Determine whether a `SerialPort` object is open.
"""
isopen(sp::SerialPort) = sp.is_open

"""
    eof(sp::SerialPort) -> Bool

Return the “end-of-file” state (`true` or `false`). Since serial ports
do not have any standard mechanism for signalling the end of a
transmitted file, this is just a dummy function that returns whatever
Boolean value `eof` was previously set with `seteof(sp, eof)`.
Returns `false` by default.
"""
eof(sp::SerialPort) = sp.is_eof

"""
    set_speed(sp::SerialPort, baudrate::Integer)

Set the data signalling rate, or the reciprocal duration of one data
bit, in bits/s. The set of values supported depends on the UART
hardware, but typically includes e.g. 9600, 19200 and 115200.

Raise an `ErrorException` if `bps` is not a value supported by the
driver or hardware.
"""
function set_speed(sp::SerialPort, baudrate::Integer)
     sp_set_baudrate(sp.ref, baudrate)
     return nothing
end

"""
    set_frame(sp::SerialPort;
              ndatabits::Integer = 8,
              parity::SPParity   = SP_PARITY_NONE,
              nstopbits::Integer = 1)

Configure the [data
framing](https://en.wikipedia.org/wiki/Universal_asynchronous_receiver/transmitter#Data_framing)
parameters. Defaults to the very common “8N1” scheme, which consists
of a start bit followed by eight data bits, no parity bit, one stop
bit.

for more details.

* `ndatabits` is the number of data bits, which is `8` in the common "8N1" scheme

* `parity` controls the presence and value of a party bit and can take
  the values `SP_PARITY_NONE` (default), `SP_PARITY_ODD`,
  `SP_PARITY_EVEN`, `SP_PARITY_MARK` and `SP_PARITY_SPACE`

* `nstopbits` sets the number of stop bits, typically `1` (default) or `2`
"""
function set_frame(sp::SerialPort;
    ndatabits::Integer=8,
    parity::SPParity=SP_PARITY_NONE,
    nstopbits::Integer=1)

    sp_set_bits(sp.ref, ndatabits)
    sp_set_parity(sp.ref, parity)
    sp_set_stopbits(sp.ref, nstopbits)
    return nothing
end

"""
    set_flow_control(sp::SerialPort;
                     rts::SPrts = SP_RTS_OFF,
                     cts::SPcts = SP_CTS_IGNORE,
                     dtr::SPdtr = SP_DTR_OFF,
                     dst::SPdsr = SP_DSR_IGNORE,
                     xonxoff::SPXonXoff = SP_XONXOFF_DISABLED)`

Configure the flow-control lines and method. Many systems don't support all options.
If an unsupported option is requested, the library will return SP_ERR_SUPP.

* `rts` controls the output line _Ready To Send (RTS)_ and can take
  the values `SP_RTS_OFF` (default), `SP_RTS_ON` and
  `SP_RTS_FLOW_CONTROL`.

* `cts` controls the input line _Clear To Send (CTS)_ and can take the
  values `SP_CTS_IGNORE` (default) and `SP_CTS_FLOW_CONTROL`

* `dtr` controls the output line _Data Terminal Ready (DTR)_ and can
  take the values `SP_DTR_OFF` (default), `SP_DTR_ON`, and
  `SP_DTR_FLOW_CONTROL`

* `dsr` controls the input line _Data Set Ready (DSR)_ and can take
  the values `SP_DSR_IGNORE` (default) and `SP_DSR_FLOW_CONTROL`

* `xonxoff` controls whether [software flow
  control](https://en.wikipedia.org/wiki/Software_flow_control) via the
  control bytes XOFF (pause transmission, 0x13, Ctrl-S) and XON (resume
  transmission, 0x11, Ctrl-Q) is active, and in which direction; it can
  take the values: `SP_XONXOFF_DISABLED` (default), `SP_XONXOFF_IN`,
  `SP_XONXOFF_OUT`, and `SP_XONXOFF_INOUT`
"""
function set_flow_control(sp::SerialPort;
    rts::SPrts=SP_RTS_OFF,
    cts::SPcts=SP_CTS_IGNORE,
    dtr::SPdtr=SP_DTR_OFF,
    dsr::SPdsr=SP_DSR_IGNORE,
    xonxoff::SPXonXoff=SP_XONXOFF_DISABLED)

    # CTS and RTS must be enabled or disabled as a pair
    sp_set_rts(sp.ref, rts)
    sp_set_cts(sp.ref, cts)

    sp_set_dtr(sp.ref, dtr)
    sp_set_dsr(sp.ref, dsr)

    sp_set_xon_xoff(sp.ref, xonxoff)
    return nothing
end

"""
    list_ports([nports_guess::Integer])`

Print a list of currently visible ports, along with some basic info.

`nports_guess` provides the number of ports guessed. Its default is `64`.
"""
function list_ports(;nports_guess::Integer=64)
    ports = sp_list_ports()

    for port in unsafe_wrap(Array, ports, nports_guess, own=false)
        port == C_NULL && return

        println(sp_get_port_name(port))
        println("\tDescription:\t",    sp_get_port_description(port))
        println("\tTransport type:\t", sp_get_port_transport(port))
    end

    sp_free_port_list(ports)
    return nothing
end

"""
    get_port_list([nports_guess::Integer])

Return a vector of currently visible ports.

`nports_guess` provides the number of ports guessed. Its default is `64`.
"""
function get_port_list(;nports_guess::Integer=64)
    ports = sp_list_ports()
    port_list = String[]
    for port in unsafe_wrap(Array, ports, nports_guess, own=false)
        port == C_NULL && return port_list
        push!(port_list, sp_get_port_name(port))
    end
    sp_free_port_list(ports)
    return port_list
end

"""
    print_port_metadata(sp::SerialPort; show_config::Bool = true)

Print info found for this port.
Note: port should be open to obtain a valid FD/handle before accessing fields.

`show_config` is `true` by default and prints out the current port settings.
"""
function print_port_metadata(sp::SerialPort; show_config::Bool=true)
    print_port_metadata(sp.ref, show_config=show_config)
    return nothing
end

function print_port_metadata(port::LibSerialPort.Port; show_config::Bool=true)
    println("\nPort name:\t",       sp_get_port_name(port))
    transport = sp_get_port_transport(port)
    print("\nPort transport:\t");
    if transport == SP_TRANSPORT_NATIVE
        println("native serial port")
    elseif transport == SP_TRANSPORT_USB
        println("USB")
        println("Manufacturer:\t",      sp_get_port_usb_manufacturer(port))
        println("Product:\t",           sp_get_port_usb_product(port))
        println("USB serial number:\t", sp_get_port_usb_serial(port))
        bus, addr = sp_get_port_usb_bus_address(port)
        println("USB bus #:\t", bus)
        println("Address on bus:\t", addr)
        vid, pid = sp_get_port_usb_vid_pid(port)
        println("Vendor ID:\t", vid)
        println("Product ID:\t", pid)
    elseif transport == SP_TRANSPORT_BLUETOOTH
        println("Bluetooth")
        println("Bluetooth address:\t", sp_get_port_bluetooth_address(port))
    end
    println("File descriptor:\t",   sp_get_port_handle(port))

    if show_config
        print_port_settings(port)
    end
    return nothing
end

function get_port_settings(config::LibSerialPort.Config)
    return Dict(
        "baudrate" => sp_get_config_baudrate(config),
        "bits"     => sp_get_config_bits(config),
        "parity"   => sp_get_config_parity(config),
        "stopbits" => sp_get_config_stopbits(config),
        "RTS"      => sp_get_config_rts(config),
        "CTS"      => sp_get_config_cts(config),
        "DTR"      => sp_get_config_dtr(config),
        "DSR"      => sp_get_config_dsr(config),
        "XonXoff"  => sp_get_config_xon_xoff(config),
        )
end

function get_port_settings(port::LibSerialPort.Port)
    config = sp_get_config(port)
    settings = get_port_settings(config)
    sp_free_config(config)
    return settings
end

"""
    get_port_settings(sp::SerialPort)

Return port settings for `sp` as a dictionary.
"""
get_port_settings(sp::SerialPort) = get_port_settings(sp.ref)

"""
    print_port_settings(sp::SerialPort)

Print port settings for `sp`.
"""
function print_port_settings(sp::SerialPort)
    print_port_settings(sp.ref)
    println("Read timeout (ms, forever if 0): ", sp.read_timeout_ms)
    return
end

function print_port_settings(port::LibSerialPort.Port)
    println("Configuration for ", sp_get_port_name(port), ":")
    config = sp_get_config(port)
    print_port_settings(config)
    sp_free_config(config)
    return
end

function print_port_settings(config::LibSerialPort.Config)
    s = get_port_settings(config)
    for (k, v) in s
        println("\t$k\t", v)
    end
    return
end

"""
    open(sp::SerialPort; mode::SPMode = SP_MODE_READ_WRITE)

Open the serial port `sp`.

`mode` selects in which direction of transmission access is requested
and can take the values: `SP_MODE_READ`, `SP_MODE_WRITE`, and
`SP_MODE_READ_WRITE` (default).
"""
function open(sp::SerialPort; mode::SPMode=SP_MODE_READ_WRITE)
    sp_open(sp.ref, mode)
    sp.is_open = true
    return sp
end

"""
    open(portname::AbstractString, baudrate::Integer;
         mode::SPMode, ndatabits::Integer,
         parity::SPParity, nstopbits::Integer)

Construct, configure and open a `SerialPort` object.

* `portname` is the name of the serial port to open. Typical port
  names available depend on the operating system. Some valid names are
  listed by [`get_port_list()`](@ref), but there are can be others and aliases:
* * Linux: `"/dev/ttyS0"`, `"/dev/ttyUSB0"`, `"/dev/serial/by-id/..."`
* * macOS: `"/dev/cu.usbserial-0001"`, `"/dev/cu.Bluetooth-Incoming-Port"`
* * Windows: `"COM1"`, `"COM2"`,`"COM3"`
* `baudrate` is the data signalling rate, or the reciprocal duration
  of one data bit, in bits/s. The set of values supported depends on
  the UART hardware, but typically includes e.g. 9600, 19200 and 115200.
* `mode` selects in which direction of transmission access is
  requested and can take the values: `SP_MODE_READ`, `SP_MODE_WRITE`,
  and `SP_MODE_READ_WRITE` (default).

The parameters `ndatabits`, `parity` and `nstopbits` have the same
meaning as in [`set_frame`](@ref) and default to the common “8N1”
frame format (8 data bits, no parity, one stop bit).

To set the flow-control method, use [`set_flow_control`](@ref).
"""
function open(portname::AbstractString,
              baudrate::Integer;
              mode::SPMode=SP_MODE_READ_WRITE,
              ndatabits::Integer=8,
              parity::SPParity=SP_PARITY_NONE,
              nstopbits::Integer=1)
    sp = SerialPort(portname)
    sp_open(sp.ref, mode)
    sp.is_open = true
    set_speed(sp, baudrate)
    set_frame(sp, ndatabits=ndatabits, parity=parity, nstopbits=nstopbits)
    return sp
end

"""
    close(sp::SerialPort)

Close the serial port `sp`.
"""
function close(sp::SerialPort)
    if isopen(sp)
        sp_close(sp.ref)
        sp.is_open = false
    end
    return
end

# fixes https://github.com/JuliaIO/LibSerialPort.jl/issues/53
@deprecate flush(sp::SerialPort, buffer::SPBuffer=SP_BUF_BOTH) sp_flush(sp, buffer)

# pass through some methods from the low-level interface
sp_output_waiting(sp::SerialPort) = sp_output_waiting(sp.ref)
sp_flush(sp::SerialPort, args...) = sp_flush(sp.ref, args...)
sp_drain(sp::SerialPort)          = sp_drain(sp.ref)

"""
    set_read_timeout(sp::SerialPort, seconds::Real)

Set a read timeout limit of `t` > 0 seconds for the total (cumulative)
time that subsequently called blocking read functions can wait before
a `Timeout` exception is thrown.

# Example

```
sp=LibSerialPort.open("/dev/ttyUSB0", 115200)
# wait until either two lines have been received
# or 10 seconds have elapsed
set_read_timeout(sp, 10)
try
    line1 = readuntil(sp, '\\n')
    line2 = readuntil(sp, '\\n')
catch e
    if isa(e, LibSerialPort.Timeout)
        println("Too late!")
    else
        rethrow()
    end
end
clear_read_timeout(sp)
```

See also: [`clear_read_timeout`](@ref), [`set_write_timeout`](@ref)
"""
function set_read_timeout(sp::SerialPort, seconds::Real)
    @assert seconds > 0
    sp.read_timeout_ms =  seconds * 1000;
    return
end

"""
    set_write_timeout(sp::SerialPort, seconds::Real)

Set a write timeout limit of `t` > 0 seconds for the total (cumulative)
time that subsequently called blocking read functions can wait before
a `Timeout` exception is thrown.

# Example

```
sp=LibSerialPort.open("/dev/ttyUSB0", 300)
# wait until either 4000 periods have been
# passed on to the serial-port driver or
# 10 seconds have elapsed
set_write_timeout(sp, 10)
try
    for i=1:50 ; write(sp, '.' ^ 80); end
catch e
    if isa(e, LibSerialPort.Timeout)
        println("This took too long!")
    else
        rethrow()
    end
end
clear_write_timeout(sp)
```

See also: [`clear_write_timeout`](@ref), [`set_read_timeout`](@ref)
"""
function set_write_timeout(sp::SerialPort, seconds::Real)
    @assert seconds > 0
    sp.write_timeout_ms =  seconds * 1000;
    return
end

"""
    clear_read_timeout(sp::SerialPort)

Cancel any previous read timeout, such that blocking read operations
will now wait without any time limit.
"""
function clear_read_timeout(sp::SerialPort)
    sp.read_timeout_ms = 0;
    return
end

"""
    clear_write_timeout(sp::SerialPort)

Cancel any previous write timeout, such that blocking write
operations will block without any time limit.
"""
function clear_write_timeout(sp::SerialPort)
    sp.write_timeout_ms = 0;
    return
end

# We define here only the basic methods for reading and writing
# bytes. All other methods for reading/writing the canonical binary
# representation of any type, and print() methods for writing
# its text representation, are inherited from the IO supertype
# (see julia/base/io.jl), i.e. work just like for files.

function read(sp::SerialPort, ::Type{UInt8})
    byte_ref = Ref{UInt8}(0)
    unsafe_read(sp, byte_ref, 1)
    return byte_ref[]
end

function unsafe_read(sp::SerialPort, p::Ptr{UInt8}, nb::UInt)
    if sp.read_timeout_ms > 0
        starttime = time_ns()
    end
    nbytes = sp_blocking_read(sp.ref, p, nb, sp.read_timeout_ms)
    if nbytes < nb
        @assert sp.read_timeout_ms > 0
        sp.read_timeout_ms = 0
        throw(Timeout())
    elseif sp.read_timeout_ms > 0
        # how much of the timeout have we used up already?
        elapsed_ms = (time_ns() - starttime + 999_999) ÷ 1_000_000
        if sp.read_timeout_ms > elapsed_ms
            sp.read_timeout_ms -= elapsed_ms
        else
            sp.read_timeout_ms = 0
            throw(Timeout())
        end
    end
    nothing
end

function write(sp::SerialPort, b::UInt8)
    unsafe_write(sp, Ref(b), 1)
end

function unsafe_write(sp::SerialPort, p::Ptr{UInt8}, nb::UInt)
    if sp.write_timeout_ms > 0
        starttime = time_ns()
    end
    nbytes = sp_blocking_write(sp.ref, p, nb, sp.write_timeout_ms)
    if nbytes < nb
        @assert sp.write_timeout_ms > 0
        sp.write_timeout_ms = 0
        throw(Timeout())
    elseif sp.write_timeout_ms > 0
        # how much of the timeout have we used up already?
        elapsed_ms = (time_ns() - starttime + 999_999) ÷ 1_000_000
        if sp.write_timeout_ms > elapsed_ms
            sp.write_timeout_ms -= elapsed_ms
        else
            sp.write_timeout_ms = 0
            throw(Timeout())
        end
    end
    nbytes
end

"""
    seteof(sp::SerialPort, state::Bool)

Set the return value of `eof(sp)` to `state`.
"""
function seteof(sp::SerialPort, state::Bool)
    sp.is_eof = state
    return nothing
end

"""
    reseteof(sp::SerialPort)

Set the return value of `eof(sp)` to its default of `false`
"""
reseteof(sp::SerialPort) = seteof(sp, false)

"""
    read(sp::SerialPort)

Return all incoming bytes currently available in the UART as a
Vector{UInt8}.

!!! note

    Julia's Base module defines `read(s::IO, nb::Integer = typemax(Int))`.
    This method overrides the default value `nb` to `bytesavailable(sp)`,
    which is useful for this context.
"""
read(sp::SerialPort) = read(sp, bytesavailable(sp))

"""
    nonblocking_read(sp::SerialPort)

Read everything from the specified serial ports `sp` input buffer, one byte at
a time, until it is empty. Returns a `String`.
"""
function nonblocking_read(sp::SerialPort)
    result = UInt8[]
    byte_ref = Ref{UInt8}(0)
    while bytesavailable(sp) > 0
        sp_nonblocking_read(sp.ref, byte_ref, 1)
        push!(result, byte_ref.x)
    end
    return result
end

"""
    bytesavailable(sp::SerialPort)

Return the number of bytes waiting in the input buffer.
"""
bytesavailable(sp::SerialPort) = sp_input_waiting(sp.ref)

end # LibSerialPort


text:
struct PointEvalHandler{dim, C, T, U}
    dh::DofHandler{dim, C, T}
    cells::Vector{Int}
    cellvalues::Vector{U}
end

# function PointEvalHandler(grid, interpolation, points)
#     """
#     This is what you would call in reality to set up the PointEvalHandler
#     """
#     kdtree = KDTree([gn.x for gn in grid.nodes])
#     """
#     Assume that we have a conformal mesh for now so that the nearest node will
#     be a part of the cell that the point is in
#     """
#     nnodes, _ = knn(kdtree, points, 1)
#     nnodes = reduce(vcat, nnodes)
# end

function PointEvalHandler(dh::DofHandler{dim, C, T}, interpolation::Interpolation{dim, S}, points::AbstractVector{Vec{dim, T}}) where {dim, C, S, T<:Real}
    grid = dh.grid
    dummy_cell_itp_qr = QuadratureRule{dim, S}(1)
    dummy_face_qr = QuadratureRule{dim-1, S}(1)
    cellvalues = CellScalarValues(dummy_cell_itp_qr, interpolation)
    facevalues = FaceScalarValues(dummy_face_qr, interpolation)
    points_array = [(i, p) for (i, p) in enumerate(points)]
    cells_of_points = Vector{Int}(undef, length(points))
    cellvalues_of_points = Vector{typeof(cellvalues)}(undef, length(points))

    for cell in 1:length(grid.cells)
        cell_coords = getcoordinates(grid, cell)
        face_nodes = faces(grid.cells[cell])
        nodes_on_faces = [grid.nodes[fn[1]] for fn in face_nodes]
        deletions = []
        for (j, ipoint) in enumerate(points_array)
            i, point = ipoint
            if is_point_inside_cell(cell_coords, point, facevalues, nodes_on_faces)
                push!(deletions, j)
                local_coordinate = find_local_coordinate(interpolation, cell_coords, point)
                point_qr = QuadratureRule{2, S, T}([1], [local_coordinate])
                cells_of_points[i] = cell
                # since these are unique for each point to evaluate, we need to create a cellvalue for each point anyway, so we might as well do it all at once
                cellvalues = CellScalarValues(point_qr, interpolation)
                reinit!(cellvalues, cell_coords)
                cellvalues_of_points[i] = cellvalues
            end
        end
        deleteat!(points_array, deletions)
    end
    if length(points_array) > 0
        println(points_array)
        error("Did not find cells for all points")
    end
    return PointEvalHandler(dh, cells_of_points, cellvalues_of_points)
end

function (peh::PointEvalHandler)(dof_values)
    [function_value(cellvalue, 1,
                    dof_values[celldofs(peh.dh, cell)]) for (cell, cellvalue) in zip(peh.cells,
                                                                                peh.cellvalues)]
end

function create_eval_matrix(peh::PointEvalHandler)
    I = reduce(vcat, [[j for i in 1:length(celldofs(peh.dh, cell))] for (j, cell) in enumerate(peh.cells)])
    J = reduce(vcat, [celldofs(peh.dh, cell) for cell in peh.cells])
    V = reduce(vcat, [[shape_value(cellvalue, 1, i) for i in 1:length(celldofs(peh.dh, cell))] for (cell, cellvalue) in zip(peh.cells, peh.cellvalues)])
    m = length(peh.cells)
    n = peh.dh.ndofs.x
    return sparse(I,J,V,m,n)
end

point_eval(peh::PointEvalHandler, dof_values) = peh(dof_values)

function is_point_inside_cell(cell_coords, point, facevalues, nodes_on_faces)
    is_on_wrong_side_of_a_plane = false

    for face in 1:length(nodes_on_faces)
        reinit!(facevalues, cell_coords, face)
        normal = facevalues.normals[1]
        face_node = nodes_on_faces[face]
        test = dot(normal, point-face_node.x) > 0
        is_on_wrong_side_of_a_plane = is_on_wrong_side_of_a_plane || test
    end
    return !is_on_wrong_side_of_a_plane
end

function find_local_coordinate(interpolation, cell_coordinates, global_coordinate)
    """
    currently copied verbatim from https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2
    other than to make J dim x dim rather than 2x2
    """
    dim = length(global_coordinate)
    local_guess = zero(Vec{dim})
    n_basefuncs = getnbasefunctions(interpolation)
    max_iters = 10
    tol_norm = 1e-10
    for iter in 1:10
        if iter == max_iters
            error("did not find a local coordinate")
        end
        N = JuAFEM.value(interpolation, local_guess)

        global_guess = zero(Vec{dim})
        for j in 1:n_basefuncs
            global_guess += N[j] * cell_coordinates[j]
        end
        residual = global_guess - global_coordinate
        if norm(residual) <= tol_norm
            break
        end
        dNdξ = JuAFEM.derivative(interpolation, local_guess)
        J = zero(Tensor{dim, dim})
        for j in 1:n_basefuncs
            J += cell_coordinates[j] ⊗ dNdξ[j]
        end
        local_guess -= inv(J) ⋅ residual
    end
    return local_guess
end


text:
function ping(gl::GitLink;
        tout = 30.0,
        verbose = false,
        onping::Function = () -> nothing
    )

    # send signal
    @info("Sending ping signal")
    ok_sync = writewdir(gl; verbose) do wdir_
        _write_ping_signal(gl, tout)
    end
    !ok_sync && (@error("Sync fail"); return false)
    @info("Ping signal sended")

    rhash0 = _remote_HEAD_hash(gl)

    # globals
    tot_t = time()
    ping_t = time()
    ping_count = 0

    client_tout = tout + min(tout * 0.5, 15.0)

    try
        @info("Waiting for response...", )

        while true
            for _ in 1:3
                if (time() - tot_t) > client_tout
                    tot_time = round(time() - tot_t; sigdigits = 3)
                    @info("Time out, total time: $(tot_time)(s)")
                    return ping_count > 0
                end

                chash = _HEAD_hash(gl)
                rhash = _remote_HEAD_hash(gl)

                # ping!
                new_remote_hash = rhash0 != rhash
                remote_ahead = rhash != chash
                if new_remote_hash && remote_ahead
                    time_ = round(time() - ping_t; sigdigits = 3)
                    ping_count += 1
                    @info("Ping $(ping_count) succeded, time: $(time_)(s)")
                    onping()
                    ping_t = time()
                    rhash0 = rhash
                end

                sleep(1.0)
            end
        end

    catch err
        (err isa InterruptException) && return
        rethrow(err)
    end

    return ping_count > 0
end


text:
using Base.Test
using AGRN

config = AGRN.Config()

@testset "Evolution functions" begin
    @testset "tournament" begin
        @test AGRN.tournament([rand()]) == 1
        for i in 2:10
            fits = rand(i)
            winner = AGRN.tournament(fits)
            @test winner >= 1
            @test winner <= i
            @test fits[winner] > minimum(fits)
        end
    end
    @testset "speciation" begin
        distance(a::Float64, b::Float64, config::AGRN.Config) = abs(a-b)
        population = rand(rand(5:15))
        reprs = rand(rand(2:5))
        species = AGRN.speciation(population, reprs, config; dfunc=distance)
        @test length(species) == length(population)
        @test all(species .<= length(reprs))
    end
    @testset "species sizes" begin
        fits = rand(config.ga_population)
        species = rand(1:5, config.ga_population)
        spec_sizes = AGRN.species_sizes(fits, species, config)
        @test length(spec_sizes) == 5
        @test sum(spec_sizes) == config.ga_population
    end
end

@testset "Evolutionary runs" begin
    @testset "Length fit" begin
        function length_fit(ind::GRN)
            length(ind)
        end
        fitness = Fitness(1, 1, length_fit)
        config = AGRN.Config()
        max_fit, best = evolve(fitness, config)
        @test max_fit > (1.0 + 1.0 + 1.0)
        @test length(best) == max_fit
    end
    @testset "sin fit" begin
        function sin_fit(ind::GRN)
            inputs = (sin.(0.0:0.1:5.0) + 1.0) / 2.0
            outs = zeros(length(inputs))
            for i in eachindex(inputs)
                outs[i] = step!(ind, [inputs[i]])[1]
            end
            -sqrt(sum((inputs - outs).^2)/length(inputs))
        end
        fitness = Fitness(1, 1, sin_fit)
        config = AGRN.Config()
        max_fit, best = evolve(fitness, config)
        @test max_fit < 0.0
        rand_grn = GRN(1, 1, 1, config)
        @test sin_fit(best) >= sin_fit(rand_grn)
    end
end


text:
module Metrics

using ..Photon: MetricStore, stop, Workout, Tensor, Meter, getmetricvalue, saveWorkout, getContext
using Printf
using Statistics
import Knet

include("core.jl")
export SmartReducer, update!, history, BinaryAccuracy, OneHotBinaryAccuracy

include("meters.jl")
export ConsoleMeter, TensorBoardMeter, FileMeter, PlotMeter

include("callbacks.jl")
export AutoSave, EpochSave, EarlyStop

end


text:
using ForceImport
VERSION > v"0.7-" ? (using Test) : (using Base.Test)

module Foo
    export +
    +() = 7
end

module Bar
    using ForceImport
    @force using Main.Foo
end

# write your own tests here
@test Bar.:+() == 7
@test 1+1 == 2


text:
using LineSearches
using LinearAlgebra

function ls(uk::Vector{Float64}, grad_L::Vector{Float64}, Lu::Function, grad_Lu::Function)
    s = copy(-grad_L)
    Lalpha(alpha) = Lu(uk .+ alpha.*s)
    function grad_Lalpha(alpha)
        return dot(grad_Lu(uk .+ alpha .* s), s)
    end
    function Lgrad_Lalpha(alpha)
        gvec = grad_Lu(uk .+ alpha .* s)
        phi = Lu(uk .+ alpha .*s)
        dphi = dot(gvec, s)
        return (phi, dphi)
    end
    dL_0 = dot(s, grad_L)
    obj = Lu(uk)
    ls_ = BackTracking()
    alpha, obj = ls_(Lalpha, grad_Lalpha, Lgrad_Lalpha, 1.0, obj, dL_0)
    return alpha
end


text:
using PowerModels
using Ipopt
using JuMP
using DataFrames
using CSV

include("post_flex_model.jl")

filename = ARGS[1]
ADC_filename = ARGS[2]

# read data files
data = PowerModels.parse_file(filename)
flex_data = CSV.read(ADC_filename)

# number of loads
num_load = size(flex_data,1)
# dictionary contining the A and B matrices
A = Dict{Int64,Any}()
B = Dict{Int64,Any}()

# populating the A and B dictionaries
for i=1:num_load
    ca = split(flex_data[:A][i][2:end-1], '|')
    cb = split(flex_data[:b][i][2:end-1], '|')
    a = zeros(length(ca),2)
    b = zeros(length(ca))
    for j=1:length(ca)
        za = split(ca[j],' ')
        a[j,:] = parse.(Float64,za)
    end
    b = parse.(Float64,cb)
    A[i] = a
    B[i] = b
end

# build base optimization model with loads as variables
m = JuMP.Model(solver=IpoptSolver())
model, var_refs = post_ac_opf_withload(data, m, A, B)

# solver the model
status = solve(model)

# extract solution and write output.csv
gen_nums = [parse(Int64,i) for (i,gen) in data["gen"]]
pg = [getvalue(var_refs["pg"][i]) for i in gen_nums]
qg = [getvalue(var_refs["qg"][i]) for i in gen_nums]

df = DataFrame()
df[:gen_num] = gen_nums
df[:active_power] = pg
df[:reactive_power] = qg

CSV.write("output.csv",df)

max_stars_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_issues_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_forks_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed


text:
# Precompose with diagonal scaling and translation

"""
  Precompose(f::ProximableFunction, a::AbstractArray, b::AbstractArray)

Returns the function `g(x) = f(diag(a)*x) + b`. Function `f` must be separable, or `a` must be a scalar, for the `prox` of `g` to be computable.
"""

immutable Precompose{T <: ProximableFunction, R <: Union{Real, AbstractArray}, S <: Union{Real, AbstractArray}} <: ProximableFunction
  f::T
  a::R
  b::S
  function Precompose(f::T, a::R, b::S)
    if !(eltype(a) <: Real && eltype(b) <: Real)
      error("a and b must be real")
    end
    if any(a == 0.0)
      error("elements of a must be nonzero")
    else
      new(f, a, b)
    end
  end
end

Precompose{T <: ProximableFunction, S <: Real}(f::T, a::S=1.0, b::S=0.0) = Precompose{T, S, S}(f, a, b)

Precompose{T <: ProximableFunction, R <: AbstractArray, S <: Real}(f::T, a::R, b::S=1.0) = Precompose{T, R, S}(f, a, b)

Precompose{T <: ProximableFunction, R <: Union{AbstractArray, Real}, S <: AbstractArray}(f::T, a::R, b::S) = Precompose{T, R, S}(f, a, b)

is_prox_accurate(f::Precompose) = is_prox_accurate(f.f)

function (g::Precompose{T, S, V}){T <: ProximableFunction, S <: Union{Real, AbstractArray}, V <: Union{Real, AbstractArray}, R <: RealOrComplex}(x::AbstractArray{R})
  return g.f((g.a).*x .+ g.b)
end

function prox!{T <: RealOrComplex, R <: Real}(g::Precompose, x::AbstractArray{T}, y::AbstractArray{T}, gamma::R=1.0)
  y .= g.a .* x .+ g.b
  v = prox!(g.f, y, y, (g.a .* g.a) .* gamma)
  y .-= g.b
  y ./= g.a
  return v
end

function prox_naive{T <: RealOrComplex, R <: Real}(g::Precompose, x::AbstractArray{T}, gamma::R=1.0)
  z = g.a .* x + g.b
  y, fy = prox_naive(g.f, z, (g.a .* g.a) * gamma)
  return (y - g.b)./g.a, fy
end

fun_name(f::Precompose) = string("Precomposition by affine diagonal mapping of ", fun_name(f.f))
fun_dom(f::Precompose) = fun_dom(f.f)
fun_expr(f::Precompose) = "x ↦ f(diag(a)*x + b)"
fun_params(f::Precompose) = string("f(x) = ", fun_expr(f.f), ", a = ", length(f.a) == 1 ? string(f.a[1]) : string(typeof(f.a)), ", b = ", length(f.b) == 1 ? string(f.b[1]) : string(typeof(f.b)))