text:
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Sample an Action
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_sample_action(Agent, t::Int; Greedy = false)
    P_Action = Agent.P_Action[:,t]
    if Greedy
        (p, act) = findmax(P_Action)
        act = rand((1:4)[P_Action.==p])
        act = act - 1
        return act
    else
        P_cum = cumsum(P_Action)
        p = rand()
        act = (P_cum[1]<p) + (P_cum[2]<p) + (P_cum[3]<p)
        return act
    end

end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Make Observation
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_make_obs(Agent, t::Int, EnvGraph::Str_EnvGraph)
    obs = Agent.Input.obs[t]
    act = Agent.Input.act[t]
    next_obs = EnvGraph.TranMat[obs+1,act+1]
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Simulating an agent
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_simulate_Agent(EnvGraph;
                             Epi=1, Block=1,
                             Param=Str_Param(),
                             HyperParam=Str_HyperParam(),
                             Agent_Prev=Str_Agent(),
                             Greedy = false,
                             T_stop = 1000)

    Input = Str_Input(; obs = [EnvGraph.InitObs[Epi]], act = [-1])
    Agent = Func_initial_Agent(; Sub=1, Epi=Epi, Block=Block,
                                Param=Param, HyperParam=HyperParam,
                                Input=Input, Agent_Prev=Agent_Prev)

    t = 1
    act = Func_sample_action(Agent,t; Greedy=Greedy)
    Agent.Input.act[t] = act

    State = Func_update_State(Agent, t)
    Agent.State[t] = State

    obs = Func_make_obs(Agent,t,EnvGraph)
    Agent.Input.obs = cat(Agent.Input.obs, obs; dims=1)

    Goal = false
    Stop = false
    while (!Goal)&&(!Stop)
        t += 1

        Stop = (t>=T_stop)
        Goal = (Agent.Input.obs[t] == EnvGraph.Goal)

        S = Func_compute_surprise(Agent, t)
        Agent.Surprise = cat(Agent.Surprise, S; dims=1)

        gamma = Func_compute_gamma(Agent,t)
        Agent.Gamma = cat(Agent.Gamma, gamma; dims=1)

        t_since_change = Func_compute_t_since_change(Agent,t)
        Agent.t_since_change = cat(Agent.t_since_change, t_since_change; dims=1)

        counter = Func_update_counter(Agent, t)
        Agent.Counter = cat(Agent.Counter, counter; dims=2)

        Π = Func_update_belief(Agent, t)
        Agent.Π = cat(Agent.Π, Π; dims=3)

        Novelty = Func_update_novelty(Agent, t)
        Agent.Novelty = cat(Agent.Novelty, Novelty; dims=2)

        # Model free updating
        δ_R = Func_compute_δ(Agent, t, "Reward")
        Agent.δ_R_Seq = cat(Agent.δ_R_Seq, δ_R; dims=1)
        δ_N = Func_compute_δ(Agent, t, HyperParam.exploration_form)
        Agent.δ_N_Seq = cat(Agent.δ_N_Seq, δ_N; dims=1)

        E_R = Func_compute_ElTrace(Agent, t, "Reward")
        Agent.E_R = cat(Agent.E_R, E_R; dims=2)
        E_N = Func_compute_ElTrace(Agent, t, HyperParam.exploration_form)
        Agent.E_N = cat(Agent.E_N, E_N; dims=2)

        Q_R_MF = Func_update_Q_MF(Agent, t, "Reward")
        Agent.Q_R_MF = cat(Agent.Q_R_MF, Q_R_MF; dims=2)
        Q_N_MF = Func_update_Q_MF(Agent, t, HyperParam.exploration_form)
        Agent.Q_N_MF = cat(Agent.Q_N_MF, Q_N_MF; dims=2)

        # Model based updating
        if (!Agent.First_Epi)||(Goal)||(Block==2)||(HyperParam.First_Epi_reward_update)
            Q_R_MB, U_R = Func_update_prior_sweep(Agent, t, "Reward")
        else
            Q_R_MB = Agent.Q_R_MB[:,1]
            U_R = Agent.U_R[:,1]
        end
        Agent.Q_R_MB = cat(Agent.Q_R_MB, Q_R_MB; dims=2)
        Agent.U_R = cat(Agent.U_R, U_R; dims=2)

        Q_N_MB, U_N = Func_update_prior_sweep(Agent, t, HyperParam.exploration_form)
        Agent.Q_N_MB = cat(Agent.Q_N_MB, Q_N_MB; dims=2)
        Agent.U_N = cat(Agent.U_N, U_N; dims=2)

        if !Goal
            if Agent.First_Epi
                P_Action = Func_compute_P_action(Agent, t; Only_Reward=false)
            else
                P_Action = Func_compute_P_action(Agent, t; Only_Reward=true)
            end
            Agent.P_Action = cat(Agent.P_Action, P_Action; dims=2)

            act = Func_sample_action(Agent,t; Greedy=Greedy)
            Agent.Input.act = cat(Agent.Input.act, act; dims=1)

            obs = Func_make_obs(Agent,t,EnvGraph)
            Agent.Input.obs = cat(Agent.Input.obs, obs; dims=1)
        else
            act = -1
            Agent.Input.act = cat(Agent.Input.act, act; dims=1)
        end

        State = Func_update_State(Agent, t)
        Agent.State = cat(Agent.State, State; dims=1)

        Agent.Novelty_Seq = cat(Agent.Novelty_Seq,
                                  Agent.Novelty[Input.obs[t]+1,t-1];
                                  dims=1)
    end

    Input = Agent.Input
    T = length(Input.obs)

    # --------------------------------------------------------------------------
    # Action Analysis - LogL
    for t = 1:(T-1)
        Agent.LogL += log(Agent.P_Action[Input.act[t]+1 ,t])

        (p, a) = findmax(Agent.P_Action[: ,t])
        if t==1
            Agent.Est_Act[1] = a-1
            Agent.Est_Act_P[1] = p
        else
            Agent.Est_Act = cat(Agent.Est_Act, a-1; dims=1)
            Agent.Est_Act_P = cat(Agent.Est_Act_P, p; dims=1)
        end
    end

    Agent.LogL_Unif = (T-1) * log( 1/Param.num_act )
    # --------------------------------------------------------------------------
    # Action Analysis - Acc Rate
    Agent.Sort_Dist_Data = zeros(size(Agent.P_Action))
    Agent.Sort_Dist_Model = zeros(size(Agent.P_Action))

    for t = 1:(T-1)
        temp = sortperm(Agent.P_Action[:,t])
        Agent.Sort_Dist_Model[:,t] = Agent.P_Action[temp,t]

        Agent.Sort_Dist_Data[Input.act[t] + 1 ,t] = 1
        Agent.Sort_Dist_Data[:,t] = Agent.Sort_Dist_Data[temp,t]

        if (Agent.Sort_Dist_Model[findmax(Agent.Sort_Dist_Data[:,t])[2],t] -
            Agent.Sort_Dist_Model[end,t]) == 0

            Agent.Sort_Dist_Data[:,t] = zeros(size(Agent.Sort_Dist_Data[:,t]))

            Agent.Sort_Dist_Data[Agent.Sort_Dist_Model[:,t] .==
                              Agent.Sort_Dist_Model[end,t],t] .= 1

            Agent.Sort_Dist_Data[:,t] = Agent.Sort_Dist_Data[:,t] ./
                                     sum(Agent.Sort_Dist_Data[:,t])
        end
    end

    Agent.LogL_Unif = (T-1) * log( 1/Param.num_act )
    Agent.Acc_Rate = mean(Agent.Sort_Dist_Data[4,:])

    return Agent, Stop
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Simulating blocks
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_simulate_Agent_Blocks(EnvGraph_Array;
                                   Epi_Set=1:5, Block_Set=1:2,
                                   Param=[Str_Param(),Str_Param()],
                                   HyperParam=Str_HyperParam(),
                                   Greedy = false,
                                   Seed = 0,
                                   T_stop = 1000)
    Random.seed!(Seed)

    N_blo = length(Block_Set)
    N_epi = length(Epi_Set)
    Agents = Array{Str_Agent,2}(undef,N_blo,N_epi)
    for Block = Block_Set
        EnvGraph = EnvGraph_Array[Block]
        for Epi = Epi_Set
            if (Block==1)&&(Epi==1)
                Agent_Epi, Stop = Func_simulate_Agent(EnvGraph;
                                                      Epi=Epi, Block=Block,
                                                      Param=Param[Block],
                                                      HyperParam=HyperParam,
                                                      Greedy = Greedy,
                                                      T_stop = T_stop)
            else
                if Epi==1
                    Agent_Prev = Agents[Block-1,end]
                else
                    Agent_Prev = Agents[Block,Epi-1]
                end

                Agent_Epi, Stop = Func_simulate_Agent(EnvGraph;
                                                      Epi=Epi, Block=Block,
                                                      Param=Param[Block],
                                                      HyperParam=HyperParam,
                                                      Greedy = Greedy,
                                                      Agent_Prev=Agent_Prev,
                                                      T_stop = T_stop)
            end

            if Stop
                return [Agent_Epi]
            end
            Agents[Block,Epi] = Agent_Epi
        end
    end

    return Agents
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# η to params SurNoR
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function Func_param_generator(η0;
                              σ = 0.1*ones(18),
                              η_low = [0,1e-12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                              η_up = [0.999,1,Inf,Inf,Inf,0.999,0.999,0.999,0.999,
                                      1,1,Inf,Inf,1,1,1,1,Inf,Inf],
                              Seed = 1,
                              noise_free=true)
    Random.seed!(Seed)
    n_par = length(η0)
    if noise_free
        η = η0
    else
        η = η0 .+ (randn(n_par).*σ)

        for i = 1:n_par
            if η[i] < η_low[i]
                η[i] = η_low[i]
            elseif η[i] > η_up[i]
                η[i] = η_up[i]
            end
        end
    end

    Param1 = Str_Param(; num_act=4, num_obs=11,
                        count_decay=1, Q_R_0 = 0)

    Param1.p_change = η[1]
    Param1.ϵ = η[2]
    Param1.β_R = η[3]
    Param1.β_ratio_N2R = η[4]
    Param1.γ_RL_R = η[6]
    Param1.γ_RL_N = η[7]
    Param1.λ_elig_R = η[8]
    Param1.λ_elig_N = η[9]
    Param1.α = η[10]
    Param1.Δα = η[11]
    Param1.Q_N_0 = η[12]
    Param1.N_prior_sweep = Int64(round(η[13]))
    Param1.w_MF_to_MB_2 = η[14]
    Param1.w_MF_to_MB_1 = η[15]
    Param1.w_MF_to_MB_Scaler = η[17]

    Param2 = deepcopy(Param1)
    Param2.β_ratio_N2R = η[5]
    Param2.w_MF_to_MB_1 = η[16]
    Param2.β_R = η[18]

    Param = [Param1, Param2]
end


text:
# # Routines for extracting segment/face data for tracks (rays) overlaid on a mesh
#
# # Return the HRPM/face indices in which each segment resides
# # Type-stable, other than warning
# function find_segment_faces(segment_points::Vector{Vector{Vector{Point_2D{F}}}},
#                             HRPM::HierarchicalRectangularlyPartitionedMesh{F, U},
#                             template_vec::MVector{N, U}
#                            ) where {F <: AbstractFloat, U <: Unsigned, N}
#
#     @info "    - Finding faces for each segment"
#     if !has_materialized_faces(HRPM)
#         @warn "Faces are not materialized for this mesh. This will be VERY slow"
#     end
#     nγ = length(segment_points)
#     bools = fill(false, nγ)
#     # Preallocate indices in the most frustrating way
#     indices =   [
#                     [
#                         [
#                             MVector{N, U}(zeros(U, N))
#                                 for i = 1:length(segment_points[iγ][it])-1 # Segments
#                         ] for it = 1:length(segment_points[iγ]) # Tracks
#                     ] for iγ = 1:nγ # Angles
#                 ]
#     Threads.@threads for iγ = 1:nγ
#         bools[iγ] = find_segment_faces_in_angle!(segment_points[iγ], indices[iγ], HRPM)
#     end
#     if !all(bools)
#         iγ_bad = findall(x->!x, bools)
#         @error "Failed to find indices for some points in segment_points$iγ_bad"
#     end
#     return indices
# end
#
# Return the face id in which each segment resides
function find_segment_faces(segment_points::Vector{Vector{Vector{Point_2D}}}, mesh::UnstructuredMesh_2D)
    @info "    - Finding faces for each segment"
    if 0 == length(mesh.materialized_faces)
        @warn "Faces are not materialized for this mesh. This will be VERY slow"
    end
    nγ = length(segment_points)
    # Preallocate indices in the most frustrating way
    segment_faces =   [
                          [
                            fill(UInt32(0), length(segment_points[iγ][it])-1) # Segments
                            for it = 1:length(segment_points[iγ]) # Tracks
                          ] for iγ = 1:nγ # Angles
                      ]
    Threads.@threads for iγ = 1:nγ
        find_segment_faces_in_angle!(segment_points[iγ], segment_faces[iγ], mesh)
        for it = 1:length(segment_points[iγ])
            if any(x->x == 0x00000000, segment_faces[iγ][it])
                @warn "Segment face not found for face in [$iγ][$it]"
            end
        end
    end
    return segment_faces
end

# # Ray trace an HRPM given the ray spacing and angular quadrature
# # Not type-stable
# function ray_trace(tₛ::F,
#                    ang_quad::ProductAngularQuadrature{nᵧ, nₚ, F},
#                    HRPM::HierarchicalRectangularlyPartitionedMesh{F, U}
#                    ) where {nᵧ, nₚ, F <: AbstractFloat, U <: Unsigned}
#     @info "Ray tracing"
#     tracks = generate_tracks(tₛ, ang_quad, HRPM)
#     segment_points = segmentize(tracks, HRPM)
#     ind_size = node_height(HRPM) + 1
#     @info "  - Using the naive segmentize + find face algorithm"
#     template_vec = MVector{ind_size, U}(zeros(U, ind_size))
#     segment_faces = find_segment_faces(segment_points, HRPM, template_vec)
#     return segment_points, segment_faces
# end
#
# Ray trace a mesh given the ray spacing and angular quadrature
function ray_trace(tₛ::Float64,
                   ang_quad::ProductAngularQuadrature{nᵧ, nₚ},
                   mesh::UnstructuredMesh_2D) where {nᵧ, nₚ}
    @info "Ray tracing"
    tracks = generate_tracks(tₛ, ang_quad, mesh, boundary_shape = "Rectangle")
    if use_E2E_raytracing(mesh)
        @info "  - Using the edge-to-edge algorithm"
        segment_points, segment_faces = ray_trace_edge_to_edge(tracks, mesh)
        validate_ray_tracing_data(segment_points, segment_faces, mesh, plot = enable_visualization)
        return segment_points, segment_faces
    else
        @info "  - Using the naive segmentize + find face algorithm"
        segment_points = segmentize(tracks, mesh)
        segment_faces = find_segment_faces(segment_points, mesh)
        validate_ray_tracing_data(segment_points, segment_faces, mesh, plot = enable_visualization)
        return segment_points, segment_faces
    end
end

# Get the segment points and the face which the segment lies in for all segments,
# in all tracks, in all angles, using the edge-to-edge ray tracing method.
# Assumes a rectangular boundary
function ray_trace_edge_to_edge(tracks::Vector{Vector{LineSegment_2D}}, mesh::UnstructuredMesh_2D)
    if visualize_ray_tracing
        @error "visualize_ray_tracing = true. Please reset to false in constants.jl"
    end
    if length(mesh.boundary_edges) != 4
        @error "Mesh does not have 4 boundary edges needed for edge-to-edge ray tracing!"
    end
    # index 1 = γ
    # index 2 = track
    # index 3 = point/segment
    nγ = length(tracks)
    segment_points =[
                        [
                            Point_2D[] for it = 1:length(tracks[iγ]) # Tracks
                        ] for iγ = 1:nγ # Angles
                    ]
    segment_faces = [
                        [
                            UInt32[] for it = 1:length(tracks[iγ]) # Tracks
                        ] for iγ = 1:nγ # Angles
                    ]
    # For each angle, get the segments and segment faces for each track
    Threads.@threads for iγ = 1:nγ
        ray_trace_angle_edge_to_edge!(tracks[iγ],
                                      segment_points[iγ],
                                      segment_faces[iγ],
                                      mesh)
    end
    return segment_points, segment_faces
end

# # Return the points of intersection between the tracks and the mesh edges.
# # Returns a Vector{Vector{Vector{Point_2D{F}}}}.
# #   index 1 = γ
# #   index 2 = track
# #   index 3 = point/segment
# # Type-stable other than the error message
# function segmentize(tracks::Vector{Vector{LineSegment_2D{F}}},
#                     HRPM::HierarchicalRectangularlyPartitionedMesh{F, U}
#                     ) where {F <: AbstractFloat, U <: Unsigned}
#
#     # Give info about intersection algorithm being used
#     int_alg = get_intersection_algorithm(HRPM)
#     @info "    - Segmentizing using the '$int_alg' algorithm"
#     # index 1 = γ
#     # index 2 = track
#     # index 3 = point/segment
#     nγ = length(tracks)
#     segment_points = Vector{Vector{Vector{Point_2D{F}}}}(undef, nγ)
#     Threads.@threads for iγ = 1:nγ
#         # for each track, intersect the track with the mesh
#         segment_points[iγ] = tracks[iγ] .∩ HRPM
#     end
#     return segment_points
# end
#
# Return the points of intersection between the tracks and the mesh edges.
# Returns a Vector{Vector{Vector{Point_2D{F}}}}.
#   index 1 = γ
#   index 2 = track
#   index 3 = point/segment
function segmentize(tracks::Vector{Vector{LineSegment_2D}}, mesh::UnstructuredMesh_2D)
    # Give info about intersection algorithm being used
    # index 1 = γ
    # index 2 = track
    # index 3 = point/segment
    nγ = length(tracks)
    segment_points = Vector{Vector{Vector{Point_2D}}}(undef, nγ)
    if length(mesh.materialized_edges) !== 0
        @info "    - Segmentizing using the 'Edges - Explicit' algorithm"
        Threads.@threads for iγ = 1:nγ
            segment_points[iγ] = intersect_edges_explicit.(tracks[iγ], Ref(mesh.materialized_edges))
        end
    elseif length(mesh.edges) !== 0
        @info "    - Segmentizing using the 'Edges - Implicit' algorithm"
        Threads.@threads for iγ = 1:nγ
            segment_points[iγ] = intersect_edges_implicit.(tracks[iγ], Ref(mesh.edges), Ref(mesh.points))
        end
    elseif length(mesh.materialized_faces) !== 0
        @info "    - Segmentizing using the 'Faces - Explicit' algorithm"
        Threads.@threads for iγ = 1:nγ
            segment_points[iγ] = intersect_faces_explicit.(tracks[iγ], Ref(mesh.materialized_faces))
        end
    else
        @info "    - Segmentizing using the 'Faces - Implicit' algorithm"
        Threads.@threads for iγ = 1:nγ
            segment_points[iγ] = intersect_faces_implicit.(tracks[iγ], Ref(mesh.faces), Ref(mesh.points))
        end
    end
    return segment_points
end

# Assumes materialized faces
function validate_ray_tracing_data(segment_points::Vector{Vector{Vector{Point_2D}}},
                                   segment_faces::Vector{Vector{Vector{UInt32}}},
                                   mesh::UnstructuredMesh_2D;
                                   plot::Bool = false)
    @info "  - Validating ray tracing data"
    nthreads = Threads.nthreads()
    nsegs = zeros(Int64, nthreads)
    nsegs_problem = zeros(Int64, nthreads) # Problem segment if the face doesn't match, and it's over 10 μm
    problem_indices = [ SVector{3, Int64}[] for i = 1:nthreads ]
    plot_segs_face = [ LineSegment_2D[] for i = 1:nthreads ]
    plot_points_face = [ Point_2D[] for i = 1:nthreads ]

    # Validate faces
    Threads.@threads for iγ = 1:length(segment_faces)
        for it = 1:length(segment_faces[iγ])
            nsegs[Threads.threadid()] += length(segment_faces[iγ][it])
            for iseg = 1:length(segment_faces[iγ][it])
                p_midpoint = midpoint(segment_points[iγ][it][iseg], segment_points[iγ][it][iseg+1])
                face = segment_faces[iγ][it][iseg]
                if p_midpoint ∉  mesh.materialized_faces[face]
                    p1 = segment_points[iγ][it][iseg]
                    p2 = segment_points[iγ][it][iseg + 1]
                    l = LineSegment_2D(p1, p2)
                    problem_length = 1e-3 < arclength(l)
                    if problem_length
                        nsegs_problem[Threads.threadid()] += 1
                        push!(problem_indices[Threads.threadid()], SVector(iγ, it, iseg))
                    end
                    # Append the points and line if we want to plot them.
                    # We only want to plot if actually a problem.
                    if enable_visualization && plot && problem_length
                        append!(plot_points_face[Threads.threadid()], [p1, p2])
                        push!(plot_segs_face[Threads.threadid()], l)
                    end
                end
            end
        end
    end

    # Attempt to fix problem segments if used E2E ray tracing
    if use_E2E_raytracing(mesh)
        @info "    - Attempting to fix $(sum(nsegs_problem)) segments with incorrect faces"
        fixed = [ Int64[] for i = 1:nthreads ]
        Threads.@threads for i = 1:nthreads
            for (iprob, problem_index) in enumerate(problem_indices[i])
                # Ray trace the track in reverse, then test to see if that fixed things
                new_info_ok = true
                iγ, it, iseg = problem_index
                p_midpoint = midpoint(segment_points[iγ][it][iseg], segment_points[iγ][it][iseg+1])
                face = segment_faces[iγ][it][iseg]
                # Check to see if still a problem. May have been fixed by a previous segment updating the
                # points/faces
                if p_midpoint ∈  mesh.materialized_faces[face]
                    nsegs_problem[i] -= 1
                    continue
                end
                npoints = length(segment_points[iγ][it])
                problem_line = LineSegment_2D(segment_points[iγ][it][1], segment_points[iγ][it][npoints])
                line_reversed = LineSegment_2D(problem_line.points[2], problem_line.points[1])
                reversed_points, reversed_faces = ray_trace_track_edge_to_edge(line_reversed,  mesh)
                # Check
                for iseg = 1:length(reversed_faces)
                    p_midpoint = midpoint(reversed_points[iseg], reversed_points[iseg+1])
                    face = reversed_faces[iseg]
                    if p_midpoint ∉  mesh.materialized_faces[face]
                        l = LineSegment_2D(reversed_points[iseg], reversed_points[iseg + 1])
                        problem_length = 1e-3 < arclength(l)
                        if problem_length
                            new_info_ok = false
                        end
                    end
                end
                if new_info_ok
                    push!(fixed[i], iprob)
                    nsegs_problem[i] -= 1
                    segment_points[iγ][it] = reverse(reversed_points)
                    segment_faces[iγ][it] = reverse(reversed_faces)
                else
                    new_face = findface(p_midpoint, mesh)
                    if new_face == 0
                        @warn "Could not find new segment face for segment $(problem_index[i])"
                    else
                        segment_faces[iγ][it][iseg] = new_face
                    end
                end
            end
        end
    end

    # Visualize
    problem_segs = sum(nsegs_problem)
    if enable_visualization && plot && 0 < problem_segs
        f = Figure()
        ax = Axis(f[1, 1], aspect = 1)
        display(f)
        mesh!(mesh.materialized_faces, color = (:black, 0.15))
        if 0 < length(mesh.materialized_edges)
            linesegments!(mesh.materialized_edges, color = :blue)
        else
            linesegments!(materialize_edges(mesh), color = :blue)
        end
        for i = 1:nthreads
            if 0 < length(fixed[i])
                deleteat!(plot_segs_face[i], fixed[i])
            end
            if 0 < length(plot_segs_face[i])
                linesegments!(plot_segs_face[i], color = :red)
            end
            if 0 < length(fixed[i])
                deleteat!(plot_points_face[i], sort!(vcat(2 .* fixed[i], 2 .* fixed[i] .- 1)))
            end
            if 0 < length(plot_points_face[i])
                scatter!(plot_points_face[i], color = :red)
            end
        end
    end
    nsegs_total = sum(nsegs)
    prob_percent = 100*problem_segs/nsegs_total
    @info "    - Problem segments: $problem_segs, Total segments: $nsegs_total"
    if 0 < problem_segs
        @info "      - Check for mesh overlap (darker gray) around areas of many problem segments."
    end
    return problem_segs == 0
end

# Plot
# -------------------------------------------------------------------------------------------------
# Plot ray tracing data one angle at a time.
if enable_visualization
#    function linesegments!(segment_points::Vector{Vector{Vector{Point_2D}}},
#                           segment_faces::Vector{Vector{Vector{UInt32}}})
#        println("Press enter to plot the segments in the next angle")
#        colormap = ColorSchemes.tab20.colors
#        lines_by_color = Vector{Vector{LineSegment_2D{T}}}(undef, 20)
#        nγ = length(segment_points)
#        for iγ = 1:nγ
#            f = Figure()
#            ax = Axis(f[1, 1], aspect = 1)
#            display(f)
#            for icolor = 1:20
#                lines_by_color[icolor] = LineSegment_2D{T}[]
#            end
#            for it = 1:length(segment_points[iγ])
#                for iseg = 1:length(segment_points[iγ][it])-1
#                    l = LineSegment_2D(segment_points[iγ][it][iseg], segment_points[iγ][it][iseg+1])
#                    face = segment_faces[iγ][it][iseg]
#                    if face == 0
#                        @error "Segment [$iγ][$it][$iseg] has a face id of 0"
#                    end
#                    push!(lines_by_color[face % 20 + 1], l)
#                end
#            end
#            for icolor = 1:20
#                linesegments!(lines_by_color[icolor], color = colormap[icolor])
#            end
#            s = readline()
#            println(iγ)
#        end
#    end

    # Set visualize_ray_tracing = true in constants.jl to get this to work.
    function plot_track_edge_to_edge(track::LineSegment_2D, mesh::UnstructuredMesh_2D)
        @info "Plotting ray tracing of track. Press enter to advance the ray"
        f = Figure()
        ax = Axis(f[1, 1], aspect = 1)
        display(f)
        linesegments!(mesh.materialized_edges, color = :blue)
        linesegments!(track, color = :orange)
        segment_points, segment_faces = ray_trace_track_edge_to_edge(track, mesh)
        return segment_points, segment_faces
    end
end

ext: jl

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
function truncated(c, k, x)
    if x ≥ c
        (x - c)^k
    else
        zero(x)
    end
end


text:
# Automatically generated using Clang.jl

const HIGHSINT_FORMAT = "d"
const HighsInt = Cint
const HighsUInt = UInt32


text:
"""
Check if IJulia is available (anywhere).
"""
function check_IJulia()
    if _check_IJulia()
        @info "IJulia seems to be installed correctly."
    else
        @warn "Couldn't find IJulia. You can call `install_IJulia()` to install it.`"
    end
    return nothing
end

function _check_IJulia()
    try
        @eval import IJulia
        return true
    catch
        return false
    end
end

"""
    install_IJulia(; python = nothing)

Install IJulia to the to the global `v#.#` environment,
or to the environment that has been given to `--project` when Julia was started.

IJulia should automatically install python and jupyter if necessary.
The keyword `python` can be used to manually point to a python executable (to give it a hint).
"""
function install_IJulia(; python = nothing)
    if !isnothing(python)
        !isfile(python) &&
        isdir(python) &&
        throw(ArgumentError("Please specify the path to the python executable " *
                            "including the executable name itself!"))
        @info "Setting ENV[\"PYTHON\"] = $python"
        ENV["PYTHON"] = python
    end

    with_pkg_env("", globalenv = true) do
        if "IJulia" in pkgs_in_env()
            @info "IJulia already seems to be installed."
        else
            @info "Installing IJulia to global environment"
            pkg"add IJulia"
        end
    end
    return nothing
end

"""
    uninstall_IJulia()

Remove IJulia from the global environment.
"""
function uninstall_IJulia()
    with_pkg_env("", globalenv = true) do
        if "IJulia" in pkgs_in_env()
            @info "Uninstalling IJulia"
            pkg"rm IJulia"
        else
            @info "IJulia couldn't be found in the global environment"
        end
    end
    return nothing
end


text:
export TransverseIsing1D, getHamiltonian, getObservables, getLocalOperators, getCorrelationFunctions, getGatesEven, getGatesOdd

struct TransverseIsing1D <: SpinHalf
    parameters::Dict{String,Any}
    sites::Vector{Index{Int64}}

    TransverseIsing1D(L::Int64) = new(
        Dict{String,Any}(
            "L" => L,
            "J" => 1.0,
            "hx" => 1.0,
            "hz" => 0.5,
            "pbc" => false,
        ),
        siteinds("S=1/2", L)
    )
end

function getHamiltonian(model::TransverseIsing1D)::MPO
    ampo = OpSum()

    L = model.parameters["L"]
    J = model.parameters["J"]
    hx = model.parameters["hx"]
    hz = model.parameters["hz"]

    for j = 1:L-1
        ampo += -4.0 * J, "Sz", j, "Sz", j + 1
    end

    if model.parameters["pbc"]
        ampo += -4.0 * J, "Sz", L, "Sz", 1
    end

    for j = 1:L
        ampo += -2.0 * hx, "Sx", j
    end

    for j = 1:L
        ampo += -2.0 * hz, "Sz", j
    end

    return MPO(ampo, model.sites)
end

function getObservables(model::TransverseIsing1D)::Vector{Observable}
    return [
        Observable("H", getHamiltonian(model)),
        Observable("Sx", getTotalSx(model)),
        Observable("Sz", getTotalSz(model))
    ]
end

function getLocalOperators(model::TransverseIsing1D)::Vector{LocalOperator}
    return [
        LocalOperator("sx", 2.0, "Sx"),
        LocalOperator("sz", 2.0, "Sz"),
    ]
end

function getCorrelationFunctions(model::TransverseIsing1D)::Vector{CorrelationFunction}
    return [
        CorrelationFunction("sz_sz", 4.0, "Sz", "Sz"),
    ]
end

function getGatesEven(model::TransverseIsing1D, dt::Float64)::Vector{ITensor}
    if model.parameters["pbc"]
        throw(DomainError(model.parameters["pbc"], "not implemented for periodic boundary conditions"))
    end

    gates = ITensor[]

    L = model.parameters["L"]
    J = model.parameters["J"]
    hx = model.parameters["hx"]
    hz = model.parameters["hz"]

    for j = 2:2:L-1
        s1 = model.sites[j]
        s2 = model.sites[j+1]

        hj = -4.0 * J * op("Sz", s1) * op("Sz", s2)
        hj -= 2.0 * hx * op("Sx", s1) * op("Id", s2)
        hj -= 2.0 * hz * op("Sz", s1) * op("Id", s2)

        push!(gates, exp(-1.0im * dt * hj))
    end

    if L % 2 == 0
        hj = -2.0 * hx * op("Sx", model.sites[L])
        hj -= 2.0 * hz * op("Sz", model.sites[L])
        push!(gates, exp(-1.0im * dt * hj))
    end

    return gates
end

function getGatesOdd(model::TransverseIsing1D, dt::Float64)::Vector{ITensor}
    if model.parameters["pbc"]
        throw(DomainError(model.parameters["pbc"], "not implemented for periodic boundary conditions"))
    end

    gates = ITensor[]

    L = model.parameters["L"]
    J = model.parameters["J"]
    hx = model.parameters["hx"]
    hz = model.parameters["hz"]

    for j = 1:2:L-1
        s1 = model.sites[j]
        s2 = model.sites[j+1]

        hj = -4.0 * J * op("Sz", s1) * op("Sz", s2)
        hj -= 2.0 * hx * op("Sx", s1) * op("Id", s2)
        hj -= 2.0 * hz * op("Sz", s1) * op("Id", s2)

        push!(gates, exp(-1.0im * dt * hj))
    end

    if L % 2 == 1
        hj = -2.0 * hx * op("Sx", model.sites[L])
        hj -= 2.0 * hz * op("Sz", model.sites[L])
        push!(gates, exp(-1.0im * dt * hj))
    end

    return gates
end


text:
#airyeigmax.jl
#Algorithm <21.1> of Random Eigenvalues by Alan Edelman

#Experiment:    Largest eigenvalue of a Stochastic Airy Operator
#Plot:          Histogram of the largest eigenvalues
#Theory:        The Tracy-Widom Law

## Parameters   # Include the most interesting configurable values at the top, with brief descriptions
t = 10000       # trials
n = 1e9         # level of discretization
beta = 2
h = n^(-1/3)    # h serves as dx

function airyeigmax_experiment(t,n,beta,h) # wrap the experimental logic in a function to enable faster JIT
    ## Experiment
    v = zeros(t)    # samples
    x = 0:h:10      # discretization of x
    N = length(x)

    #Parallel experiment
    b = (1 / h^2) * ones(N-1)
    v = pmap(1:t) do i
        ### The kernel is equivalent to evolve(AiryProcess(beta), h, 10)
        ## discretize stochastic airy operator
        # discretize airy operator
        a = -(2 / h ^ 2) * ones(N);  # differential operator: d^2 / dx^2
        a -= x                       # d^2 / dx^2 - x
        # add the stochastic part
        dW = randn(N) * sqrt(h)
        a += (2 / sqrt(beta)) * dW / h
        ## calculate the largest eigenvalue of tridiagonal matrix T
        # diagonal of T: a
        # subdiagonal of T: b
        return eigmax(SymTridiagonal(a,b))
    end

    binsize = 1/6
    return hist(v, -6:binsize:6)
end
grid, count = airyeigmax_experiment(t,n,beta,h) #run the experiment, making the global variables local for speed

## Plot
include("../1/tracywidom.jl")
h = Histogram(count/(sum(count)*step(grid)), step(grid))    # add the histogram of data
h.x0 = first(grid)
add(p, h)

if isinteractive()
    Winston.display(p)
else
    file(p, "<filename>.png")
end


text:
using PyPlot
A,B,C=1,2,3 # give variables a number
Nlevels=10

function demoSkill()
    A,B,C=1,2,3 # give variables a number
    Nlevels=10

    sigmoid(x)=1./(1+exp(-x))

    pLeftPlayerWins=zeros(Nlevels,Nlevels)
    for levelleft=1:10
        for levelright=1:10
            pLeftPlayerWins[levelleft,levelright]=sigmoid(levelleft-levelright)
        end
    end

    #prior_pot=exp.(-0.05*(Nlevels/2-linspace(0,Nlevels,Nlevels)).^2) # peaked around skill 5
    prior_pot=ones(Nlevels,1) # uniform prior on skill
    prior_pot=prior_pot./sum(prior_pot)
    priorA=PotArray(A,prior_pot)
    priorB=PotArray(B,prior_pot)
    priorC=PotArray(C,prior_pot)

    fAbeatsB=PotArray([A B],pLeftPlayerWins)
    fBbeatsC=PotArray([B C],pLeftPlayerWins)
    fAbeatsC=PotArray([A C],pLeftPlayerWins)

    Data = "A>B, B>C"
    pABC=fAbeatsB*fBbeatsC*priorA*priorB*priorC # p(A,B,C,Data)
    plotstuff(Data,pABC,fAbeatsC)

    Data = "A>B, A>B"
    pABC=fAbeatsB*fAbeatsB*priorA*priorB*priorC # p(A,B,C,Data)
    plotstuff(Data,pABC,fAbeatsC)

    Data = "(A>B)^6, (B>C)^6"
    pABC=fAbeatsB^6*fBbeatsC^6*priorA*priorB*priorC # p(A,B,C,Data)
    plotstuff(Data,pABC,fAbeatsC)

    Data = "A>B, C>B"
    pABC=fAbeatsB*(1-fBbeatsC)*priorA*priorB*priorC # p(A,B,C,Data)
    plotstuff(Data,pABC,fAbeatsC)

end

function plotstuff(Data,pABC,fAbeatsC)
    pABC=pABC./sum(pABC) # normalise to get p(A,B,C|Data)
    pAC=sum(pABC,B)
    pAbeatsC=sum(fAbeatsC*pAC)
    # What is the probability that A will beat C?
    println(Data*" : probability that A beats C = $(pAbeatsC.content)")
    pA=sum(pABC,[B C]); pB=sum(pABC,[A C]); pC=sum(pABC,[B A])
    figure(); title(Data)
    plot(pA.content,label="A skill")
    plot(pB.content,label="B skill")
    plot(pC.content,label="C skill")
    legend()
    # Can calculate mean skill levels using:
    println("mean skill for A=",sum(pA.content.*(1:Nlevels))) # mean skill level for A
    println("mean skill for B=",sum(pB.content.*(1:Nlevels))) # mean skill level for B
    println("mean skill for C=",sum(pC.content.*(1:Nlevels))) # mean skill level for C
end


text:
@testset "User defined model" begin
    y = [1 2; 1 2.0] # y definition

    Z = [1 2.0; 2.0 1.0]
    T = [1 2.0]
    R = [1 2.0]
    @test_throws ErrorException StateSpaceModel(y, Z, T, R)

    Z = Array{Float64, 3}(undef, 2, 2, 2)
    Z[:, :, 1] = [1 2.0; 2.0 1.0]
    Z[:, :, 2] = [1 2.0; 2.0 1.0]
    T = [1 2.0; 2.0 1.0]
    R = [1 2.0; 2.0 1.0]
    ss = StateSpaceModel(y, Z, T, R)
    @test ss.mode == "time-variant"

    Z = Array{Float64, 3}(undef, 2, 2, 2)
    Z[:, :, 1] = [1 2.0; 2.0 1.0]
    Z[:, :, 2] = [1 2.0; 2.0 1.0]
    T = [1 2.0]
    R = [1 2.0; 2.0 1.0]
    @test_throws ErrorException StateSpaceModel(y, Z, T, R)
end


text:
module BuzzHash

export sprand_fd, buzzhash, inverse

""" sprand_fd(KC, PN, p)

    Returns a KC by PN sparse 1/0 (technically Int8) matrix in which each entry is 1 with probability p. KC and PN abbreviate Kenyon Cells and Projection Neurons, respectively, to acknowledge the algorithm's biological source. Generally KC is much larger than PN (e.g., 2000 vs 50) and p is relatively small (e.g., 0.12).

    (I.e., sprand_fd is essentially a wrapper for Base.sprand.)
    """
function sprand_fd{T1<:Integer, T2<:Integer, R<:Real}(KC::T1, PN::T2, p::R)::SparseMatrixCSC{Int8,Int}
    tmp = sprand(Bool,KC,PN,p) # Bool to make it binary
    i,j = findnz(tmp)
    return sparse(i,j,Int8(1)) # Convert to Int8
end

""" sprand_fd(KC, PN, PN_per_KC)

    If the last argument is an integer sprand_fd returns a KC by PN sparse 1/0 (technically Int8) matrix in which each row has exactly PN_per_KC randomly placed 1's. This version is much slower.
    """
function sprand_fd{T1<:Integer, T2<:Integer, T3<:Integer}(KC::T1, PN::T2, PN_per_KC::T3)::SparseMatrixCSC{Int8,Int}
    KC > 0 && PN > 0 || error("KC and PN must be positive.")
    tmp = collect(1:PN)
    idx = 1:PN_per_KC
    ans = spzeros(Int8,KC,PN)
    for i in 1:KC
        shuffle!(tmp)
        ans[i,tmp[idx]]=Int8(1)
    end
    return ans
end

""" buzzhash(A, x, topN; clip=true)

    Return a sparse vector formed by applying A to x-mean(x), zeroizing all but the topN largest components of the result, and by default setting the remaining components to 1 ("clipping" them.) To retain the component values, call with `clip=false`. A is a random, sparse, binary (technically Int8) matrix implementing the expansion stage of the algorithm.
"""
function buzzhash{R<:Real,T<:Integer}(A::SparseMatrixCSC{Int8,Int}, x::Vector{R}, topN::T; clip::Bool=true)::SparseVector{Float64,Int}
    # center x and apply A
    tmp = A*(x-mean(x))
    # find the permutation which sorts tmp in descending order
    perm = sortperm(tmp, rev=true)
    # return a sparse vector containing the largest topN elements of the result.
    return sparsevec(perm[1:topN], clip ? 1 : tmp[perm[1:topN]], size(A,1))
end

""" inverse(A)

    Returns a closure (function + data) which, when applied to Ax, returns x.

    When KC >> PN the PN columns of a sparse random matrix formed by sprand_fd are linearly independent with high probability. This implies that the map, y = Ax, is invertable--that x can be recovered from y. The inverse mapping is inv(A'A)A' where A' is the transpose of A and inv indicates inverse. The inversion will be at best approximate, of course, after the final hashing step in which all but the largest values are zeroized.
    """
function inverse(A::SparseMatrixCSC{Int8,Int})
    B = Matrix{Float64}(A)
    S = inv(B'*B)
    (x) -> S*(B'*x)
end

end # module

max_stars_repo_name: WilCrofter/BuzzHash

max_issues_repo_name: WilCrofter/BuzzHash

max_forks_repo_name: WilCrofter/BuzzHash


text:
module tinygrad

import Base: +, -, *, /, ^, inv, exp, log, sin, cos, abs

export  +, -, *, /, ^, inv, exp, log, sin, cos, abs,
        Diffable, backwards!, clear!

include("core.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ f354c7da-63b5-11eb-154d-0191edfd3c93
begin
	using Pkg
	Pkg.add(["DataFrames", "RDatasets"])

	using DataFrames
	using RDatasets
end

# ╔═╡ 23ae525c-63bb-11eb-31d1-cdc26c21546b
begin
	Pkg.add(["DataFramesMeta", "Query"])

	using DataFramesMeta
	using Query
end

# ╔═╡ b9b81dfc-63be-11eb-2982-29f73ad29cf2
using Statistics

# ╔═╡ 9e6e4bce-63b5-11eb-3e5a-0912e0975882
md"# Работа с датафреймами на Julia"

# ╔═╡ c24bea42-63b5-11eb-1aad-cd94dc0c4679
md"## Создание датафрейма"

# ╔═╡ d070c570-63b5-11eb-1946-bd81d539224e
md"Для начала скачаем датасеты из R и пакет DataFrames для работы с ними:"

# ╔═╡ 1a7d0f34-63b6-11eb-2d21-735f0df9ad0f
md"Сохраняем датасет diamonds как DataFrame:"

# ╔═╡ 57bbbd5a-63b6-11eb-177a-edd02c67a162
diamonds = dataset("ggplot2", "diamonds")

# ╔═╡ 27e1ebc6-63b7-11eb-1d1d-9b2c6c337410
md"Функции first и last позволяют вывести на экран первые и последние строки датафрейма соответственно:"

# ╔═╡ 3b03aa00-63b7-11eb-1412-1b4ed91ec1ff
first(diamonds, 5)

# ╔═╡ 4e6f3802-63b7-11eb-3046-a15ea8996d7a
last(diamonds, 5)

# ╔═╡ 5977749e-63b7-11eb-314c-0f803af2384d
md"## Выбор отдельных столбцов"

# ╔═╡ 6a54d3ec-63b7-11eb-3ae3-b5c85d118271
md"По названию стобцы могут быть выбраны одним из следующих способов:"

# ╔═╡ a642aa28-63b7-11eb-1087-05b49af3a8f9
diamonds.Carat

# ╔═╡ d03a4eba-63b7-11eb-3d33-4383c57dcbbb
diamonds."Carat"

# ╔═╡ d8239908-63b7-11eb-1c7d-0d1a01f91907
diamonds[!, :Carat]

# ╔═╡ e22ebc02-63b7-11eb-0eb0-91077a72127d
diamonds[!, "Carat"]

# ╔═╡ f518b20a-63b7-11eb-3471-b336d3d7e887
md"Поскольку df[!, :col] не создает копию, изменение элементов получаемого таким способом вектора приводит к изменению значений в исходном df."

# ╔═╡ 0a0fba6e-63b8-11eb-2deb-d97aa7162d23
begin
	vec = diamonds[!, "Carat"]
	vec[1] = 111 # изменяем первое значение вектора на 111
	vec
end

# ╔═╡ 209a82d0-63b8-11eb-07cc-eb324848c7af
 first(diamonds, 5) # значение в df изменилось на 111

# ╔═╡ 2bfebcd8-63b8-11eb-1ba5-7bae39541004
md"Для создания копии столбца надо использовать df[:, :col] , изменение полученного вектора не будет влиять на исходный df:"

# ╔═╡ 3a2c2502-63b8-11eb-0269-4f17ae7948d1
begin
	vec2 = diamonds[:, "Carat"]
	vec2[1] = 222 # изменяем первое значение вектора на 222
	vec2
end

# ╔═╡ 4e022a7a-63b8-11eb-1d81-b99f33e6de2e
first(diamonds, 5) # значение в df не изменилось на 222

# ╔═╡ 61578892-63b8-11eb-060f-bffb99a079f6
diamonds[1, 1] = 0.23 # вернем исходное значение

# ╔═╡ 6e273bb2-63b8-11eb-0859-25605875e8bd
md"## Выбор подмножества датафрейма"

# ╔═╡ 7ca723d2-63b8-11eb-3482-df3ad4448ca7
md"Можно выбрать подмножество датафрейма с помощью индексов. Двоеточие : означает, что
нужно выбрать все элементы (строки или столбца, в зависимости от позиции)."

# ╔═╡ 87f55de4-63b8-11eb-394a-5fde7ade3192
diamonds[2:10, :]

# ╔═╡ 912f32c0-63b8-11eb-2bdd-f540663f6e66
diamonds[[1, 5, 10], :]

# ╔═╡ a4f9c902-63b8-11eb-20c6-f5b161e12ae3
diamonds[3:6, ["Carat", "Color", "Cut"]]

# ╔═╡ b7c21206-63b8-11eb-1b7d-39ed5f2a1ae1
diamonds[10:15, [:Carat, :Color, :Cut]]

# ╔═╡ c902c736-63b8-11eb-1427-f7168b69de43
md"Стоит отметить, что df[:, [:col]] и df[!, [:col]] возвращают объект DataFrame, а df[:, :col] и df[!, :col] - вектор:"

# ╔═╡ db91a390-63b8-11eb-13ab-2b7680a2114b
diamonds[:, [:Carat]]

# ╔═╡ ecdd9186-63b8-11eb-331c-cf9af704178c
diamonds[:, :Carat]

# ╔═╡ f7882602-63b8-11eb-0c90-f14fc95e323f
md"Селектор Not позволяет выбрать все столбцы, кроме обозначенного подмножества:"

# ╔═╡ 04e3bfda-63b9-11eb-1411-917898a9416b
diamonds[1:5, Not("Carat")]

# ╔═╡ 16882e1a-63b9-11eb-08b4-673fd8d60167
diamonds[1:5, Not(["Carat", "Cut"])]

# ╔═╡ 293aa204-63b9-11eb-3377-df1c4e6b9aa3
md"Селектор Between позволяет выбрать все столбцы между указанными:"

# ╔═╡ 3458fef6-63b9-11eb-043e-7950b0d1ce4b
diamonds[1:5, Between("X", "Z")]

# ╔═╡ 41bfbfce-63b9-11eb-3315-f1ad0a828c67
md"С помощью индексов можно осуществлять выбор наблюдений (строк), удовлетворяющих
условиям на значения переменных (столбцов):"

# ╔═╡ 4aceb928-63b9-11eb-38a9-0f4ec57d95f6
diamonds[diamonds.Carat .> 0.5, :] # все наблюдения со значением "Carat" строго больше 0.5

# ╔═╡ 5a022d9e-63b9-11eb-0b0a-3dc16ed2b1ce
# все наблюдения, где "Carat" строго больше 0.5 и "Price" строго меньше 2500
diamonds[(diamonds.Carat.> 0.5) .& (diamonds.Price .<2500) , :]

# ╔═╡ 734f54de-63b9-11eb-27c5-a135601ba3ad
# все наблюдения, где "Carat" строго больше 0.5 и значение "Cut" равно "Ideal"
diamonds[(diamonds.Carat.> 0.5) .& (diamonds.Cut .== "Ideal") , :]

# ╔═╡ 8be7b8c4-63b9-11eb-0e5d-effb5856e5ce
# все наблюдения, где "Carat" строго больше 0.5 или "Price" строго больше 2750
diamonds[(diamonds.Carat .> 0.5) .| (diamonds.Price .> 2750), :]

# ╔═╡ 9e4b1012-63b9-11eb-33bb-d766e60e61e2
md"## Работа со столбцами"

# ╔═╡ d19b6690-63b9-11eb-362b-550e7e0bb190
md"### select, select!, transform и transform!"

# ╔═╡ dcd6f150-63b9-11eb-19a4-0304f869203a
md"С помощью функций select и select! можно осуществлять выбор, переименование и
транформацию столбцов датафрейма."

# ╔═╡ e6fa7cba-63b9-11eb-0770-6d2cdd08cf6d
md"Функция select создает новый объект DataFrame. Она всегда возвращает объект DataFrame, даже если выбран только один столбец:"

# ╔═╡ fb3a1a52-63b9-11eb-0bc6-e1fba2ada82f
select(diamonds, :Carat) # датафрейм с одним столбцом

# ╔═╡ 088dea6a-63ba-11eb-2aea-8bd4410071b3
select(diamonds, ["Carat", "Cut", "Price"]) # выбор нескольких указанных столбцов

# ╔═╡ 1584a62a-63ba-11eb-1ccf-1b8651770f75
select(diamonds, Not(["Clarity", "Color"])) # выбор всех столбцов, кроме указанных

# ╔═╡ 21b772fe-63ba-11eb-287a-b36c72b79084
select(diamonds, Between("Carat", "Price")) # выбор диапазона столбцов

# ╔═╡ 2f0f250a-63ba-11eb-187f-b1c93806addb
# переименование столбцов
select(diamonds, "Carat" => "Карат", "Price" => "Цена", "Cut" => "Огранка")

# ╔═╡ 50de6d94-63ba-11eb-0b43-ed16faff9dd8
# трансформация столбца
select(diamonds, "Price", "Price" => (x -> log.(x)) => "LogPrice")

# ╔═╡ 68064a0c-63ba-11eb-3bbf-17fa3851f382
# альтернативный способ сделать то же самое
select(diamonds, "Price", "Price" => ByRow(log))

# ╔═╡ 7659aca0-63ba-11eb-2efe-3d952047ee7d
md"Чтобы не создавать копию датафрейма, можно использовать функцию select!"

# ╔═╡ 7faeaeb6-63ba-11eb-103c-cbb4fbf4f24b
md"(для изменений inplace)"

# ╔═╡ 89a7418c-63ba-11eb-2e20-1fac8245672d
md"Функции transform и transform! работают аналогично select и select!. Единственное
отличие - сохраняются все исходные столбцы:"

# ╔═╡ 99a30cb0-63ba-11eb-04b2-f9e62b4bfb98
# добавляем новый столбец (логарифм цены) в исходный датафрейм
transform!(diamonds, "Price" => ByRow(log))

# ╔═╡ b179754a-63ba-11eb-3c7b-55e009489380
diamonds # в датафрейм добавился столбец

# ╔═╡ bc086efa-63ba-11eb-065f-cd5ff0803f53
md"# Фреймворки для работы с данными и запросами"

# ╔═╡ c94842d2-63ba-11eb-1249-dd3341cdcd1a
md"## DataFramesMeta"

# ╔═╡ d193b4bc-63ba-11eb-2bc4-63ec49ef1d71
md"Догрузим и запустим необходимые пакеты:"

# ╔═╡ 63cfff56-63bc-11eb-322b-8329cb04fded
md"Рассматриваемый пакет предоставляет короткий способ вызова столбцов DataFrame."

# ╔═╡ 7635a9d4-63bc-11eb-2b6a-2f2e1ffb6b6a
md"Вызов столбца тут осуществляется только по его названию. Помимо этого можно создавать цепи вызовом макроса @linq: оператор |> подает результат вычисления предыдущего шага на вход в следующее выражение. Например:"

# ╔═╡ b6547496-63bc-11eb-3f77-034e29bdd11d
@linq diamonds |>
	where(:Price .> 500) |>
	select(Mass=:Carat, :Price)

# ╔═╡ c13e21ae-63bc-11eb-0323-f58fa9e7783a
md"Данный фреймворк также позволяет делить данные из исходного датасета на группы вычислять функции от групп, комбинировать полученные результаты."

# ╔═╡ cf752a42-63bc-11eb-0450-51a28f855353
md"Например, для группировки качественных данных (типа CategoricalValue) можно использовать функцию sort(), а для создания нового столбца на основе имеющихся – функцию transform(). Чтобы применить функцию ко всему столбцу (поточечно) непосредственно после ее названия нужно поставить знак точки:"

# ╔═╡ dcfd4b72-63bc-11eb-29ba-e5a7fa7ee50d
@linq diamonds |>
	sort(:Cut) |>
	transform(LogRatio = log.(:Price ./ :Carat))

# ╔═╡ e7a9393c-63bc-11eb-0da9-ad6ab03f30f3
md"## Query"

# ╔═╡ f9441f54-63bc-11eb-1273-affc5e83fde8
md"Другой способ извлекать необходимые паттерны из датасета – пакет для обработки запросов Query."

# ╔═╡ 05d1599e-63bd-11eb-1899-d387dfe903cf
md"Запрос начинается с макроса @from. Первый аргумент – это переменная, пробегающая строки таблицы, второй – непосредственно запрашиваемая таблица."

# ╔═╡ 107934c0-63bd-11eb-309d-1f3c517deb3c
md"Команда @where отфильтровывает строки, для которых значение в указанном столбце не
удовлетворяет заданному условию."

# ╔═╡ 1dc4a2ae-63bd-11eb-235c-1fd93d071330
md"Команда @select сохраняет требуемые данные в новый столбец. Чтобы выход имел вид таблицы следует использовать {}."

# ╔═╡ 23a47848-63bd-11eb-0519-b9d191b54902
md"При помощи команды @collect можно задать тип объекта, желаемый получить в качестве ответа на запрос (по умолчанию возвращается iterator)."

# ╔═╡ 2ea3cf80-63bd-11eb-10bf-910fd48c9ba0
md"Тело запроса (часть кода после @from) начинается с команды begin и заканчивается командой end."

# ╔═╡ 3985089e-63bd-11eb-38e8-a7dab8fa51d6
md"Выполним те же операции, что и в предыдущем разделе, но с ипользованием фреймворка
Query:"

# ╔═╡ 44604774-63bd-11eb-0228-81c4346b1d16
q1 = @from i in diamonds begin
	@where i.Price > 500
	@select {Mass=i.Carat, i.Price}
	@collect DataFrame
end

# ╔═╡ 53069b3e-63bd-11eb-2b7f-5baf50ae000d
md"Вызов команды @collect без указания типа данных вернет ответ на запрос в формате array:"

# ╔═╡ 5be87394-63bd-11eb-33ec-cda3cfcae091
q2 = @from i in diamonds begin
	@where i.Price > 500 && i.Carat > 0.5
	@select i.Price
	@collect
end

# ╔═╡ 6c6f26ac-63bd-11eb-0350-57481a29969f
md"# Summarizing"

# ╔═╡ 6916c56a-63be-11eb-3fb4-fbf373f00a61
md"Команда describe показывает статистическую информацию о датафрейме, среди которой такие значения: среднее, минимальное, максимальное значения, медиана, количество пропусков, тип переменной:"

# ╔═╡ 7e4b52f0-63be-11eb-35b0-5b05e05e6136
describe(diamonds)

# ╔═╡ 86ea35e0-63be-11eb-2fe2-3dcd50b8e279
md"Данную команду можно применить также не только к целому датафрейму, но и к его отдельной части, например:"

# ╔═╡ 9d718ab6-63be-11eb-1e39-350dd83b1cf3
describe(diamonds[!, [:Cut]])

# ╔═╡ a8b3c574-63be-11eb-27de-3515fec4c748
md"Также с помощью пакета Statistics можно найти статистические значения для отдельного столбца. Например, среднее:"

# ╔═╡ c5be9950-63be-11eb-0404-0f06dd409a46
mean(diamonds[!, "Carat"]) #среднее

# ╔═╡ cdf71a48-63be-11eb-2c94-db3981fcfe15
md"Аналогично можно найти другие значения:"

# ╔═╡ d64eb566-63be-11eb-36df-a16db49233ac
median(diamonds[!, "Price"]) #медиана

# ╔═╡ e07ea2b2-63be-11eb-1a8f-55af47d59a88
maximum(diamonds[!, "Price"]) #максимум

# ╔═╡ e8143922-63be-11eb-3a93-4bf933abe948
minimum(diamonds[!, "Price"]) #минимум

# ╔═╡ ee11aece-63be-11eb-24d9-1d4b90a7a6b9
md"Функция combine позволяет применять функцию к отдельному столбцу. Например, следующим образом:"

# ╔═╡ fe101e5a-63be-11eb-20d7-7147d3d6afdb
combine(diamonds, :Price .=> sum) #сумма по столбцу

# ╔═╡ 05c3c7dc-63bf-11eb-1620-a90dbedf4157
combine(diamonds, :Price => maximum) #максимальное значение по столбцу

# ╔═╡ 0c46b7cc-63bf-11eb-061a-31586565022d
md"Select позволяет получить тот же результат, только на выходе дает то же число строк, которое изначально в датафрейме:"

# ╔═╡ 1677f7e2-63bf-11eb-059f-83ea78ead39f
select(diamonds, :Price => maximum) #максимальное значение по столбцу

# ╔═╡ 226fbada-63bf-11eb-3cd0-5d07d6c8c1e8
md"# Работа со столбцами"

# ╔═╡ c1d18022-63bf-11eb-21a9-876544dd3553
md"Все функции, которые меняют датафрейм по умолчанию копируют его столбцы. Так, например:"

# ╔═╡ cd273b24-63bf-11eb-3126-1f37ce4cdc6e
diamonds2 = copy(diamonds)

# ╔═╡ df2ce4de-63bf-11eb-1f00-c78758bb98e3
diamonds2.Price === diamonds.Price

# ╔═╡ ee638900-63bf-11eb-250a-53ed8a5523cd
md"Тем не менее, будьте осторожны! Функции, которые заканчиваются символом ! могут изменить столбец, который ей скормили в качестве аргумента. Обратите внимание, к примеру, на такую ситуацию:"

# ╔═╡ 17a3d89c-63c0-11eb-2ae1-ebd583514f0e
price = diamonds.Price

# ╔═╡ 203bc4e2-63c0-11eb-33fb-9bc93e0fd74f
df = DataFrame(price=price)

# ╔═╡ 285682fc-63c0-11eb-064e-2b34fafa81ce
sort!(df)

# ╔═╡ 3e0eebb4-63c0-11eb-3168-edbe95b535cb
price

# ╔═╡ 4a7e4a16-63c0-11eb-169a-fbcb1fc5c874
df.price[1] = 100

# ╔═╡ 5215f69a-63c0-11eb-0947-bf7b827cb4cd
df

# ╔═╡ 5a8d0dae-63c0-11eb-2a3b-83519757fae9
price

# ╔═╡ 60fbe23a-63c0-11eb-0090-1b7c7c60af19
md"Как Вы могли заметить, в указанном выше примере вектор цен price, взятый из таблицы diamonds не был изменён, так как конструктор DataFrame(price=price) по умолчанию создаёт копию."

# ╔═╡ 6b3d3160-63c0-11eb-0ac0-8388c4d47603
md"Тем не менее, если Вам вдруг понадобится обратиться напрямую к некоторому столбцу датафрейма, это возможно сделать при помощи функции eachcol."

# ╔═╡ 7875630c-63c0-11eb-241f-cb8e76e78df1
df2 = DataFrame(price=price)

# ╔═╡ 8181233c-63c0-11eb-2828-87b66090b047
df2.price == price

# ╔═╡ 896d9526-63c0-11eb-2f0e-f790d633ce92
df2[!, 1] !== price

# ╔═╡ 912a93a4-63c0-11eb-0e3b-bfc907b230cc
eachcol(df2)[1] === df2.price

# ╔═╡ 985256f8-63c0-11eb-2856-b975188c3fa1
md"Будьте осторожны, когда меняете столбцы из датафрейма таким способом!"

# ╔═╡ Cell order:
# ╟─9e6e4bce-63b5-11eb-3e5a-0912e0975882
# ╟─c24bea42-63b5-11eb-1aad-cd94dc0c4679
# ╟─d070c570-63b5-11eb-1946-bd81d539224e
# ╠═f354c7da-63b5-11eb-154d-0191edfd3c93
# ╟─1a7d0f34-63b6-11eb-2d21-735f0df9ad0f
# ╠═57bbbd5a-63b6-11eb-177a-edd02c67a162
# ╟─27e1ebc6-63b7-11eb-1d1d-9b2c6c337410
# ╠═3b03aa00-63b7-11eb-1412-1b4ed91ec1ff
# ╠═4e6f3802-63b7-11eb-3046-a15ea8996d7a
# ╟─5977749e-63b7-11eb-314c-0f803af2384d
# ╟─6a54d3ec-63b7-11eb-3ae3-b5c85d118271
# ╠═a642aa28-63b7-11eb-1087-05b49af3a8f9
# ╠═d03a4eba-63b7-11eb-3d33-4383c57dcbbb
# ╠═d8239908-63b7-11eb-1c7d-0d1a01f91907
# ╠═e22ebc02-63b7-11eb-0eb0-91077a72127d
# ╟─f518b20a-63b7-11eb-3471-b336d3d7e887
# ╠═0a0fba6e-63b8-11eb-2deb-d97aa7162d23
# ╠═209a82d0-63b8-11eb-07cc-eb324848c7af
# ╟─2bfebcd8-63b8-11eb-1ba5-7bae39541004
# ╠═3a2c2502-63b8-11eb-0269-4f17ae7948d1
# ╠═4e022a7a-63b8-11eb-1d81-b99f33e6de2e
# ╠═61578892-63b8-11eb-060f-bffb99a079f6
# ╟─6e273bb2-63b8-11eb-0859-25605875e8bd
# ╟─7ca723d2-63b8-11eb-3482-df3ad4448ca7
# ╠═87f55de4-63b8-11eb-394a-5fde7ade3192
# ╠═912f32c0-63b8-11eb-2bdd-f540663f6e66
# ╠═a4f9c902-63b8-11eb-20c6-f5b161e12ae3
# ╠═b7c21206-63b8-11eb-1b7d-39ed5f2a1ae1
# ╟─c902c736-63b8-11eb-1427-f7168b69de43
# ╠═db91a390-63b8-11eb-13ab-2b7680a2114b
# ╠═ecdd9186-63b8-11eb-331c-cf9af704178c
# ╟─f7882602-63b8-11eb-0c90-f14fc95e323f
# ╠═04e3bfda-63b9-11eb-1411-917898a9416b
# ╠═16882e1a-63b9-11eb-08b4-673fd8d60167
# ╟─293aa204-63b9-11eb-3377-df1c4e6b9aa3
# ╠═3458fef6-63b9-11eb-043e-7950b0d1ce4b
# ╟─41bfbfce-63b9-11eb-3315-f1ad0a828c67
# ╠═4aceb928-63b9-11eb-38a9-0f4ec57d95f6
# ╠═5a022d9e-63b9-11eb-0b0a-3dc16ed2b1ce
# ╠═734f54de-63b9-11eb-27c5-a135601ba3ad
# ╠═8be7b8c4-63b9-11eb-0e5d-effb5856e5ce
# ╟─9e4b1012-63b9-11eb-33bb-d766e60e61e2
# ╟─d19b6690-63b9-11eb-362b-550e7e0bb190
# ╟─dcd6f150-63b9-11eb-19a4-0304f869203a
# ╟─e6fa7cba-63b9-11eb-0770-6d2cdd08cf6d
# ╠═fb3a1a52-63b9-11eb-0bc6-e1fba2ada82f
# ╠═088dea6a-63ba-11eb-2aea-8bd4410071b3
# ╠═1584a62a-63ba-11eb-1ccf-1b8651770f75
# ╠═21b772fe-63ba-11eb-287a-b36c72b79084
# ╠═2f0f250a-63ba-11eb-187f-b1c93806addb
# ╠═50de6d94-63ba-11eb-0b43-ed16faff9dd8
# ╠═68064a0c-63ba-11eb-3bbf-17fa3851f382
# ╟─7659aca0-63ba-11eb-2efe-3d952047ee7d
# ╟─7faeaeb6-63ba-11eb-103c-cbb4fbf4f24b
# ╟─89a7418c-63ba-11eb-2e20-1fac8245672d
# ╠═99a30cb0-63ba-11eb-04b2-f9e62b4bfb98
# ╠═b179754a-63ba-11eb-3c7b-55e009489380
# ╟─bc086efa-63ba-11eb-065f-cd5ff0803f53
# ╟─c94842d2-63ba-11eb-1249-dd3341cdcd1a
# ╟─d193b4bc-63ba-11eb-2bc4-63ec49ef1d71
# ╠═23ae525c-63bb-11eb-31d1-cdc26c21546b
# ╟─63cfff56-63bc-11eb-322b-8329cb04fded
# ╟─7635a9d4-63bc-11eb-2b6a-2f2e1ffb6b6a
# ╠═b6547496-63bc-11eb-3f77-034e29bdd11d
# ╟─c13e21ae-63bc-11eb-0323-f58fa9e7783a
# ╟─cf752a42-63bc-11eb-0450-51a28f855353
# ╠═dcfd4b72-63bc-11eb-29ba-e5a7fa7ee50d
# ╟─e7a9393c-63bc-11eb-0da9-ad6ab03f30f3
# ╟─f9441f54-63bc-11eb-1273-affc5e83fde8
# ╟─05d1599e-63bd-11eb-1899-d387dfe903cf
# ╟─107934c0-63bd-11eb-309d-1f3c517deb3c
# ╟─1dc4a2ae-63bd-11eb-235c-1fd93d071330
# ╟─23a47848-63bd-11eb-0519-b9d191b54902
# ╟─2ea3cf80-63bd-11eb-10bf-910fd48c9ba0
# ╟─3985089e-63bd-11eb-38e8-a7dab8fa51d6
# ╠═44604774-63bd-11eb-0228-81c4346b1d16
# ╟─53069b3e-63bd-11eb-2b7f-5baf50ae000d
# ╠═5be87394-63bd-11eb-33ec-cda3cfcae091
# ╟─6c6f26ac-63bd-11eb-0350-57481a29969f
# ╟─6916c56a-63be-11eb-3fb4-fbf373f00a61
# ╠═7e4b52f0-63be-11eb-35b0-5b05e05e6136
# ╟─86ea35e0-63be-11eb-2fe2-3dcd50b8e279
# ╠═9d718ab6-63be-11eb-1e39-350dd83b1cf3
# ╟─a8b3c574-63be-11eb-27de-3515fec4c748
# ╠═b9b81dfc-63be-11eb-2982-29f73ad29cf2
# ╠═c5be9950-63be-11eb-0404-0f06dd409a46
# ╟─cdf71a48-63be-11eb-2c94-db3981fcfe15
# ╠═d64eb566-63be-11eb-36df-a16db49233ac
# ╠═e07ea2b2-63be-11eb-1a8f-55af47d59a88
# ╠═e8143922-63be-11eb-3a93-4bf933abe948
# ╟─ee11aece-63be-11eb-24d9-1d4b90a7a6b9
# ╠═fe101e5a-63be-11eb-20d7-7147d3d6afdb
# ╠═05c3c7dc-63bf-11eb-1620-a90dbedf4157
# ╟─0c46b7cc-63bf-11eb-061a-31586565022d
# ╠═1677f7e2-63bf-11eb-059f-83ea78ead39f
# ╟─226fbada-63bf-11eb-3cd0-5d07d6c8c1e8
# ╟─c1d18022-63bf-11eb-21a9-876544dd3553
# ╟─cd273b24-63bf-11eb-3126-1f37ce4cdc6e
# ╠═df2ce4de-63bf-11eb-1f00-c78758bb98e3
# ╟─ee638900-63bf-11eb-250a-53ed8a5523cd
# ╠═17a3d89c-63c0-11eb-2ae1-ebd583514f0e
# ╠═203bc4e2-63c0-11eb-33fb-9bc93e0fd74f
# ╠═285682fc-63c0-11eb-064e-2b34fafa81ce
# ╠═3e0eebb4-63c0-11eb-3168-edbe95b535cb
# ╠═4a7e4a16-63c0-11eb-169a-fbcb1fc5c874
# ╠═5215f69a-63c0-11eb-0947-bf7b827cb4cd
# ╠═5a8d0dae-63c0-11eb-2a3b-83519757fae9
# ╟─60fbe23a-63c0-11eb-0090-1b7c7c60af19
# ╟─6b3d3160-63c0-11eb-0ac0-8388c4d47603
# ╠═7875630c-63c0-11eb-241f-cb8e76e78df1
# ╠═8181233c-63c0-11eb-2828-87b66090b047
# ╠═896d9526-63c0-11eb-2f0e-f790d633ce92
# ╠═912a93a4-63c0-11eb-0e3b-bfc907b230cc
# ╟─985256f8-63c0-11eb-2856-b975188c3fa1


text:
#####
##### Some multivariate test problems
#####

####
#### generic code
####

"""
A type for test functions with uniform bounds.

All subtypes (test functions) support:

1. being called with a vector of arbitrary length

2. `minimum_location`,

3. `lower_bounds`, `upper_bounds` (defined via `lu_bounds`),

4. having a global minimum *value* of `0`.

See [`TEST_FUNCTIONS`](@ref).
"""
abstract type UniformBounds end

lower_bounds(f::UniformBounds, n::Integer) = fill(Float64(lu_bounds(f)[1]), n)

upper_bounds(f::UniformBounds, n::Integer) = fill(Float64(lu_bounds(f)[2]), n)

####
#### shifted quadratic
####

"A test function with global minimum at [a, …]. For sanity checks."
struct ShiftedQuadratic{T <: Real} <: UniformBounds
    a::T
end

(f::ShiftedQuadratic)(x) = (a = f.a; sum(x -> abs2(x - a), x))

minimum_location(f::ShiftedQuadratic, n::Integer) = fill(f.a, n)

lu_bounds(f::ShiftedQuadratic) = (f.a - 50, f.a + 100)

const SHIFTED_QUADRATIC = ShiftedQuadratic(0.5)

####
#### Griewank
####

"""
The Griewank problem.

From Ali, Khompatraporn, and Zabinsy (2005, p 559).
"""
struct Griewank{T <: Real} <: UniformBounds
    a::T
end

function (f::Griewank)(x)
    sum(abs2, x) / f.a - prod(((i, x),) -> cos(x / √i), enumerate(x)) + 1
end

minimum_location(::Griewank, n::Integer) = zeros(n)

lu_bounds(::Griewank) = (-50, 100)

const GRIEWANK = Griewank(200.0)

####
#### Levy Montalvo 2.
####

"""
Levi and Montalvo 2 problem.

From Ali, Khompatraporn, and Zabinsy (2005, p 662).
"""
struct LevyMontalvo2 <: UniformBounds end

function (::LevyMontalvo2)(x)
    xn = last(x)
    0.1 * abs2(sinpi(3 * first(x))) + abs2(xn -  1) * (1 + abs2(sinpi(2 * xn))) +
        sum(@. abs2($(x[1:(end - 1)]) - 1) * (1 + abs2(sinpi(3 * $(x[2:end])))))
end

minimum_location(::LevyMontalvo2, n::Integer) = ones(n)

lu_bounds(::LevyMontalvo2) = (-10, 15)

const LEVY_MONTALVO2 = LevyMontalvo2()

####
#### Rastrigin
####

"""
Rastrigin problem.

From Ali, Khompatraporn, and Zabinsy (2005, p 665).
"""
struct Rastrigin <: UniformBounds end

function (::Rastrigin)(x)
    10 * length(x) + sum(@. abs2(x) - 10 * cospi(2 * x))
end

minimum_location(::Rastrigin, n::Integer) = zeros(n)

lu_bounds(::Rastrigin) = (-5.12, 5.12)

const RASTRIGIN = Rastrigin()

####
#### Rosenbrock
####

"""
Rosenbrock problem.

From Ali, Khompatraporn, and Zabinsy (2005, p 666).
"""
struct Rosenbrock <: UniformBounds end

function (::Rosenbrock)(x)
    x1 = x[1:(end - 1)]
    x2 = x[2:end]
    sum(@. 100 * abs2(x2 - abs2(x1)) + abs2(x1 - 1))
end

minimum_location(::Rosenbrock, n::Integer) = ones(n)

lu_bounds(::Rosenbrock) = (-30, 30)

const ROSENBROCK = Rosenbrock()

####
#### helper code
####

"A tuple of all test functions."
const TEST_FUNCTIONS = (SHIFTED_QUADRATIC, GRIEWANK, LEVY_MONTALVO2, RASTRIGIN, ROSENBROCK)


text:
#
# Factors-DataFrames conversion
#
# Bridge between the two

"""
Convert a Factor to a DataFrame
"""
function Base.convert(::Type{DataFrame}, ϕ::Factor)
    df = DataFrames.DataFrame(pattern(ϕ))
    DataFrames.rename!(df, [f => t for (f, t) = zip(names(df), names(ϕ))] )
    df[:potential] = ϕ.potential[:]

    return df
end

function DataFrame(ϕ::Factor)
    return convert(DataFrame, ϕ)
end


text:
# Autogenerated wrapper script for GNUMake_jll for i686-w64-mingw32
export make

JLLWrappers.@generate_wrapper_header("GNUMake")
JLLWrappers.@declare_executable_product(make)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        make,
        "bin\\make.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Practical Machine learning
# Regression Analysis  - Logistic Regression example
# Chapter 10

## Initialization
using PyPlot
using NLopt

function costFunction(theta, X, y)

  # Initialize some useful values
  m = length(y) # number of training examples

  # You need to return the following variables correctly
  J = COSTFUNCTION(theta, X, y)
  grad = zeros(size(theta)[1])
  return (J, grad)
end

function plotData(X, y)

  # Create New Figure
  figure()
  hold(true)
  PLOTDATA(x,y)
  hold(false)
end

function plotDecisionBoundary(theta, X, y)
  # Plot Data
  plotData(X[:, 2:3], y)
  hold(true)

  if size(X, 2) <= 3
    # Only need 2 points to define a line, so choose two endpoints
    plot_x = [minimum(X[:, 2])-2,  maximum(X[:, 2])+2]

    # Calculate the decision boundary line
    plot_y = (-1./theta[3]).*(theta[2].*plot_x + theta[1])

    # Plot, and adjust axes for better viewing
    plot(plot_x, plot_y)

    # Legend, specific for the exercise
    # legend("Admitted", "Not admitted", "Decision Boundary")
    legend("Admitted", "Not admitted")
    axis([30, 100, 30, 100])
  else
    # Here is the grid range
    u = linspace(-1, 1.5, 50)
    v = linspace(-1, 1.5, 50)

    z = zeros(length(u), length(v))
    # Evaluate z = theta*x over the grid
    for j in 1:length(v), i in 1:length(u)
      z[i, j] = (mapFeature([u[i]], [v[j]]) * theta)[1]
    end
    z = z' # important to transpose z before calling contour

    # Plot z = 0
    # Notice you need to specify the range [0, 0]
    contour(u, v, z, [0, 0], linewidth=2)
  end
  hold(false)
end

function sigmoid(z)
  J = SIGMOID(z)
  return J
end

function  predict(theta, X)
  m = size(X, 1) # Number of training examples
  p = PREDICT(theta, X)
 return p
end

## Load Data
#  The first two columns contains the exam scores and the third column
#  contains the label.

data = readcsv("dataset1.txt")
X = data[:, 1:2]
y = data[:, 3]

## ==================== Plotting ====================

@printf("Plotting data with + indicating (y = 1) examples and o indicating (y = 0) examples.\n")

plotData(X, y)

# Put some labels
hold(true)
# Labels and Legend
xlabel("Exam 1 score")
ylabel("Exam 2 score")

# Specified in plot order
legend(["Admitted", "Not admitted"])
hold(false)

@printf("Program paused. Press enter to continue.\n")
readline()

## ============ Compute Cost and Gradient ============

#  Setup the data matrix appropriately, and add ones for the intercept term
m, n = size(X)

# Add intercept term to x and X_test
X = [ones(m, 1) X]

# Initialize fitting parameters
initial_theta = zeros(n + 1)

# Compute and display initial cost and gradient
cost, grad = costFunction(initial_theta, X, y)

@printf("Cost at initial theta (zeros): %f\n", cost)
@printf("Gradient at initial theta (zeros): \n")
@printf("%s", join(map(x -> @sprintf(" %f ", x), grad), "\n"))

@printf("\nProgram paused. Press enter to continue.\n")
readline()

## ============= Part 3: Optimizing using fminunc  =============
#  In this exercise, you will use a built-in function (fminunc) to find the
#  optimal parameters theta.

#  Set options for fminunc
#options = optimset("GradObj", "on", "MaxIter", 400)
options = Opt(:LN_NELDERMEAD, n+1)
min_objective!(options, (theta, grad) -> costFunction(theta, X, y)[1])
maxeval!(options, 400)
#  Run fminunc to obtain the optimal theta
#  This function will return theta and the cost
(cost, theta, _) = optimize(options, initial_theta)

# Print theta to screen
@printf("Cost at theta found by optimize: %f\n", cost)
@printf("theta: \n")
@printf("%s\n", join(map(x -> @sprintf(" %f ", x), theta), "\n"))

# Plot Boundary
plotDecisionBoundary(theta, X, y)

# Put some labels
hold(true)
# Labels and Legend
xlabel("Exam 1 score")
ylabel("Exam 2 score")

# Specified in plot order
legend(["Regression", "Admitted", "Not admitted"])
hold(false)

@printf("\nProgram paused. Press enter to continue.\n")
#readline()

prob = sigmoid([1 45 85] * theta)
@printf("For a student with scores 45 and 85, we predict an admission probability of %f\n\n", prob[1])

# Compute accuracy on our training set
p = predict(theta, X)

@printf("Train Accuracy: %f\n", mean(p .== y) * 100)


text:
# # Perfect convective adjustment calibration with Unscented Ensemble Kalman Inversion
#
# This example calibrates a convective adjustment model in the "perfect model context".
# In this context, synthetic observations are generated by a convective adjustment model
# with "true" parameters. The true parameters are then "rediscovered" by calibrating the model
# to match the synthetic observations.
#
# We use the discrepency between observed and modeled buoyancy ``b`` to calibrate
# the convective adjustment model.
# The calibration problem is solved by Ensemble Kalman Inversion.
# For more information about Ensemble Kalman Inversion, see the
# [EnsembleKalmanProcesses.jl documentation](https://clima.github.io/EnsembleKalmanProcesses.jl/stable/ensemble_kalman_inversion/).
#
# The calibration is done here using Unscented Kalman Inversion. For more information about the
# algorithm refer to [EnsembleKalmanProcesses.jl documentation](https://clima.github.io/EnsembleKalmanProcesses.jl/stable/unscented_kalman_inversion/).

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add OceanTurbulenceParameterEstimation, Oceananigans, Distributions, CairoMakie"
# ```

# We load some packages,

using OceanTurbulenceParameterEstimation
using LinearAlgebra
using EnsembleKalmanProcesses.ParameterDistributionStorage

# and reuse some some code from a previous example to generate observations,

examples_path = joinpath(pathof(OceanTurbulenceParameterEstimation), "..", "..", "examples")
include(joinpath(examples_path, "intro_to_inverse_problems.jl"))

data_path = generate_synthetic_observations()
observations = SyntheticObservations(data_path, field_names=:b, normalize=ZScore)

# ### Free parameters
#
# We construct some prior distributions for our free parameters. We found that it often helps to
# constrain the prior distributions so that neither very high nor very low values for diffusivities
# can be drawn out of the distribution.

priors = (
    convective_κz = ConstrainedNormal(0.0, 1.0, 0.0, 4 * θ★.convective_κz),
    background_κz = ConstrainedNormal(0.0, 1.0, 0.0, 4 * θ★.background_κz)
)

free_parameters = FreeParameters(priors)

# and an ensemble_simulation,

Nparameters = length(priors)
Nensemble = 2 * Nparameters + 1

ensemble_simulation, closure★ = build_ensemble_simulation(observations; Nensemble)

# To visualize the prior distributions we randomly sample out values from then and plot the p.d.f.

using CairoMakie
using OceanTurbulenceParameterEstimation.EnsembleKalmanInversions: convert_prior, inverse_parameter_transform

samples(prior) = [inverse_parameter_transform(prior, x) for x in rand(convert_prior(prior), 50000000)]

samples_convective_κz = samples(priors.convective_κz)
samples_background_κz = samples(priors.background_κz)

f = Figure()

axtop = Axis(f[1, 1],
             xlabel = "convective_κz [m² s⁻¹]",
             ylabel = "p.d.f.")

axbottom = Axis(f[2, 1],
                xlabel = "background_κz [m² s⁻¹]",
                ylabel = "p.d.f.")

densities = []

push!(densities, density!(axtop, samples_convective_κz))
push!(densities, density!(axbottom, samples_background_κz))

save("visualize_prior_diffusivities_convective_adjustment_uki.svg", f); nothing #hide

# ![](visualize_prior_diffusivities_convective_adjustment_uki.svg)

# ### The inverse problem

# We can construct the inverse problem ``y = G(θ) + η``. Here, ``y`` are the `observations` and ``G`` is the
# `ensemble_model`.
calibration = InverseProblem(observations, ensemble_simulation, free_parameters)

# Next, we construct an `UnscentedKalmanInversion` (UKI) object,

prior_mean = fill(0.0, Nparameters)
prior_cov = Matrix(Diagonal(fill(1.0, Nparameters)))
α_reg = 1.0   # regularization parameter
update_freq = 1
noise_covariance = 0.05^2  # error is about 5%

uki = UnscentedKalmanInversion(calibration, prior_mean, prior_cov;
                               noise_covariance = noise_covariance, α_reg = α_reg, update_freq = update_freq)

# and perform few iterations to see if we can converge to the true parameter values.

iterate!(uki; iterations = 10)

# Last, we visualize the outputs of UKI calibration:

θ_mean, θθ_cov, θθ_std_arr, error =  UnscentedKalmanInversionPostprocess(uki)

N_iter = size(θ_mean, 2)

f = Figure(resolution = (800, 600))
ax1 = Axis(f[1, 1],
           xlabel = "iterations",
           xticks = 1:N_iter,
           ylabel = "convective_κz [m² s⁻¹]")

ax2 = Axis(f[2, 1],
           xlabel = "iterations",
           xticks = 1:N_iter,
           ylabel = "background_κz [m² s⁻¹]")

band!(ax1, 1:N_iter, θ_mean[1, :] .+ θθ_std_arr[1, :], θ_mean[1, :] .- θθ_std_arr[1, :])
lines!(ax1, Float64.(1:N_iter), θ_mean[1, :])
hlines!(ax1, [θ★.convective_κz], color=:red)

band!(ax2, 1:N_iter, θ_mean[2, :] .+ θθ_std_arr[2, :], θ_mean[2, :] .- θθ_std_arr[2, :])
lines!(ax2, 1:N_iter, θ_mean[2, :])
hlines!(ax2, [θ★.background_κz], color=:red)

xlims!(ax1, 0.5, N_iter+0.5)
xlims!(ax2, 0.5, N_iter+0.5)

save("uki_results.svg", f); nothing #hide

# ![](uki_results.svg)

# and the error:

f = Figure(resolution = (800, 300))

lines(f[1, 1], 2:N_iter, error, color = :red, linewidth = 2,
               axis = (xlabel = "iterations",
                       xticks = 1:N_iter,
                       ylabel = "error"))

xlims!(0.5, N_iter+0.5)

save("uki_error.svg", f); nothing #hide

# ![](uki_error.svg)


text:
function compare(a, b)
    println("\n$a is of type $(typeof(a)) and $b is of type $(typeof(b))")
    if a <  b println("$a is strictly less than $b") end
    if a <= b println("$a is less than or equal to $b") end
    if a >  b println("$a is strictly greater than $b") end
    if a >= b println("$a is greater than or equal to $b") end
    if a == b println("$a is equal to $b") end
    if a != b println("$a is not equal to $b") end
    if a === b println("$a has object identity with $b") end
    if a !== b println("$a has negated object identity with $b") end
end

compare("YUP", "YUP")
compare('a', 'z')
compare("24", "123")
compare(24, 123)
compare(5.0, 5)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module ModalC0BasesTests

using Test
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials
# using BenchmarkTools

import Gridap.Fields: Broadcasting

# Real-valued Q space with isotropic order

x1 = Point(0.0)
x2 = Point(0.5)
x3 = Point(1.0)

V = Float64
G = gradient_type(V,x1)
H = gradient_type(G,x1)
order = 12

order = 3
a = fill(Point(-0.5),order+1)
b = fill(Point(2.5),order+1)
b1 = ModalC0Basis{1}(V,order,a,b)
∇b1 = Broadcasting(∇)(b1)
∇∇b1 = Broadcasting(∇)(∇b1)

@test evaluate(b1,[x1,x2,x3,]) ≈ [1.0 0.0 -0.0 0.0;
                                  0.5 0.5 -0.4330127018922193 0.18633899812498247;
                                  0.0 1.0 -0.0 -0.0]
@test evaluate(∇b1,[x1,x2,x3,]) ≈ G[(-1.0,) (1.0,) (-1.7320508075688772,) (1.4907119849998598,);
                                    (-1.0,) (1.0,) (-0.0,) (-0.37267799624996495,);
                                    (-1.0,) (1.0,) (1.7320508075688772,) (-0.0,)]
@test evaluate(∇∇b1,[x1,x2,x3,]) ≈ H[(0.0,) (0.0,) (3.4641016151377544,) (-5.962847939999439,);
                                     (0.0,) (0.0,) (3.4641016151377544,) (-1.4907119849998598,);
                                     (0.0,) (0.0,) (3.4641016151377544,) (2.9814239699997196,)]

x1 = Point(0.0,0.0)
x2 = Point(0.5,0.5)
x3 = Point(1.0,1.0)
a = [ Point(0.0,0.0),Point(0.0,0.0),Point(0.0,0.0),Point(0.0,0.0),
      Point(-0.5,2.5),Point(-0.5,2.5),Point(0.0,1.5),Point(0.0,1.5),
      Point(-1.0,-1.0),Point(-1.0,-1.0),Point(-1.0,-1.0),Point(-1.0,-1.0),
      Point(0.0,0.0),Point(0.0,0.0),Point(0.0,0.0),Point(0.0,0.0) ]
b = [ Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0),
      Point(2.5,2.5),Point(2.5,2.5),Point(1.5,1.5),Point(1.5,1.5),
      Point(-1.0,1.0),Point(-1.0,1.0),Point(-1.0,1.25),Point(-1.0,1.25),
      Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0) ]
b2 = ModalC0Basis{2}(V,order,a,b)
∇b2 = Broadcasting(∇)(b2)
∇∇b2 = Broadcasting(∇)(∇b2)

G = gradient_type(V,x1)
H = gradient_type(G,x1)

@test evaluate(b2,[x1,x2,x3,]) ≈ [ 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
                                   0.25 0.25 0.25 0.25 -0.21650635094610965 0.09316949906249124 -0.21650635094610965 0.09316949906249124 -0.21650635094610965 -0.13975424859373686 -0.21650635094610965 -0.09316949906249122 0.18749999999999997 0.0 0.0 0.0;
                                   0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
@test evaluate(∇b2,[x1,x2,x3,])[:,10] ≈ G[ (0.0, 0.0,); (0.2795084971874737, -0.2795084971874737,); (0.0, 0.0,) ]
@test evaluate(∇∇b2,[x1,x2,x3,])[:,10] ≈ H[ (0.0, 0.0, 0.0, -4.47213595499958);
                                            (0.0, 0.5590169943749475, 0.5590169943749475, 1.118033988749895);
                                            (0.0, -2.23606797749979, -2.23606797749979, 0.0) ]

end # module


text:
col_negs(A::Matrix) = [sum(A[:,j].<0) for j=1:size(A)[2]]
row_negs(A::Matrix) = [sum(A[i,:].<0) for i=1:size(A)[1]]
negs(A) = (row_negs(A), col_negs(A))

straight_cols(A) = all(sort(col_negs(A)).==Vector(0:size(A)[2]-1))
straight_rows(A) = all(sort(row_negs(A)).==Vector(0:size(A)[1]-1))
is_straight(A) = straight_cols(A) && straight_rows(A)

crossing_cols(A) = any(sort(col_negs(A)).>Vector(0:size(A)[2]-1))
crossing_rows(A) = any(sort(row_negs(A)).>Vector(0:size(A)[1]-1))
is_crossing(A) = crossing_cols(A) || crossing_rows(A)

⊖(a::Vector, b::Vector, ⊟₁= >, ⊟₂= -) = (a.⊟₁b)⋅(a.⊟₂b)

Δᶜ⁺(A::Matrix) = sort(col_negs(A)) ⊖ Vector(0:size(A)[2]-1)
Δʳ⁺(A::Matrix) = reverse(sort(row_negs(A))) ⊖ Vector(size(A)[1]-1:-1:0)

Δ⁺(A) = max(Δᶜ⁺(A), Δʳ⁺(A))

Δᶜ⁻(A::Matrix) = Vector(0:size(A)[2]-1) ⊖ sort(col_negs(A))
Δʳ⁻(A::Matrix) = Vector(size(A)[1]-1:-1:0) ⊖ reverse(sort(row_negs(A)))

col_negs_sums(A::Matrix) = [sum(zeros(eltype(A),size(A)[2])⊖A[:,j]) for j=1:size(A)[2]]
row_negs_sums(A::Matrix) = [sum(zeros(eltype(A),size(A)[1])⊖A[i,:]) for i=1:size(A)[1]]
negs_sums(A) = (row_negs_sums(A), col_negs_sums(A))


text:
# Tests for EXP3 agents

@testset "EXP" begin

noOfArms = 5

# EXP3 agnet
@testset "EXP3" begin

    srand(1729)
    agent = EXP3( noOfArms, 0.25 )
    @test agent.noOfArms == noOfArms
    for t = 1:10
        armIdx  = get_arm_index( agent )
        update_reward!( agent, rand()  )
    end
    @test agent.noOfSteps == 10

    @test info_str( agent ) == "EXP3 (γ = 0.250)"
    @test info_str( agent, false ) == "EXP3 (γ = 0.250)"
    @test info_str( agent, true ) == "EXP3 (\$\\gamma = 0.250\$)"

    reset!( agent )
    @test ( agent.noOfArms==noOfArms &&
            agent.noOfSteps==0 &&
            agent.γ==0.25 &&
            agent.wVec==ones(noOfArms)  &&
            agent.pDist.p==1/noOfArms*ones(noOfArms) )

    # With γ = 1.0, EXP3 will be like uniform strategy
    srand(1729)
    agent = EXP3( noOfArms, 1.00 )
    for t = 1:10
        armIdx = get_arm_index( agent )
        update_reward!( agent, rand() )
    end
    @test( agent.pDist.p == 1/noOfArms*ones(noOfArms) )
end

# EXP3.1 agent
@testset "EXP3.1" begin
    srand(1729)
    agent = EXP31( noOfArms )
    # No need to test sub agent EXP3 here again
    @test info_str( agent ) == "EXP3.1"
    @test info_str( agent, false ) == "EXP3.1"
    @test info_str( agent, true ) == "EXP3.1"

    reset!( agent )
    @test ( agent.G_hat == zeros(noOfArms) &&
            agent.r == 0 &&
            agent.g_r == (noOfArms*log(noOfArms))/(e-1) )

end

# REXP3 agent
@testset "REXP3" begin
    srand(1729)
    agent = REXP3( noOfArms, 0.25, 10 )

    @test info_str( agent ) == "REXP3 (γ = 0.250, Δ = 10)"
    @test info_str( agent, false ) == "REXP3 (γ = 0.250, Δ = 10)"
    @test info_str( agent, true ) == "REXP3 (\$\\gamma = 0.250, \\Delta = 10\$)"

    # Take 9 steps and agent shouldn't be reset
    for t = 1:9
        get_arm_index( agent )
        update_reward!( agent, rand() )
    end
    @test agent._EXP3.noOfSteps == 9
    @test agent.noOfSteps == 9
    @test agent.j == 1
    # Take one more step and agent should reset
    get_arm_index( agent ); update_reward!( agent, rand() )
    @test agent._EXP3.noOfSteps == 0
    @test agent.noOfSteps == 10
    @test agent.j == 2
    # Take one more step and it should not reset
    get_arm_index( agent ); update_reward!( agent, rand() )
    @test agent._EXP3.noOfSteps == 1
    @test agent.noOfSteps == 11
    @test agent.j == 2

end

# EXP3-IX agent
@testset "EXP3-IX" begin

    srand(1729)
    agent = EXP3IX( noOfArms, 0.20, 0.15 )

    @test info_str( agent ) == "EXP3-IX (η = 0.200, γ = 0.150)"
    @test info_str( agent, false ) == "EXP3-IX (η = 0.200, γ = 0.150)"
    @test info_str( agent, true ) == "EXP3-IX (\$\\eta = 0.200, \\gamma = 0.150\$)"

    for t = 1:10
        get_arm_index( agent ); update_reward!( agent, rand() )
    end
    @test agent.noOfSteps == 10

    reset!( agent )
    @test agent.noOfSteps == 0 &&
            agent.lastPlayedArm == 0 &&
            agent.wVec == ones( noOfArms ) &&
            agent.pDist.p == 1/noOfArms*ones(noOfArms) &&
            agent.η == 0.200 &&
            agent.γ == 0.150

    # With η = 0, weights shouldn't change
    agent = EXP3IX( noOfArms, 0.00, 0.10 )
    for t = 1:5
        get_arm_index( agent ); update_reward!( agent, rand() )
    end
    @test agent.wVec == ones( noOfArms )
    @test agent.pDist.p == 1/noOfArms * ones(noOfArms)
end

end


text:
"""
    cl4(x::Float64)::Float64

Returns the value of the Clausen function ``\\operatorname{Cl}_4(x)``
for a real angle ``x`` of type `Float64`.  This function is defined as

```math
\\operatorname{Cl}_4(x) = \\Im[\\operatorname{Li}_4(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\sin(kx)}{k^4}
```

Author: Alexander Voigt

License: MIT

# Example
```jldoctest; setup = :(using ClausenFunctions), output = false
julia> cl4(1.0)
0.8958052386793799
```
"""
function cl4(x::Float64)::Float64
    zeta3 = 1.2020569031595943
    pi28 = pi*pi/8.0

    (x, sgn) = range_reduce_even(x)

    if x == zero(x) || x == pi
        return zero(x)
    end

    if x < 0.5*pi
        P = (-3.0555555555555556e-01,  6.0521392328447206e-03,
             -1.9587493942041528e-05, -3.1137343767030358e-08)
        Q = (1.0000000000000000e+00, -2.2079728398400851e-02,
             1.0887447112236682e-04, -6.1847621370547954e-08)
        y = x*x
        y2 = y*y
        p = P[1] + y * P[2] + y2 * (P[3] + y * P[4])
        q = Q[1] + y * Q[2] + y2 * (Q[3] + y * Q[4])
        sgn*x*(zeta3 + y*(p/q + 1/6*log(x)))
    else
        P = (7.6223911686491336e-01, -2.4339587368267260e-01,
             2.8715364937979943e-02, -1.5368612510964667e-03,
             3.6261044225761673e-05, -2.8557977333851308e-07)
        Q = (1.0000000000000000e+00, -1.7465715261403233e-01,
             9.5439417991615653e-03, -1.7325070821666274e-04,
             5.9283675098376635e-07,  9.4127575773361230e-10)
        y = pi - x
        z = y*y - pi28
        z2 = z*z
        z4 = z2*z2
        p = P[1] + z * P[2] + z2 * (P[3] + z * P[4]) + z4 * (P[5] + z * P[6])
        q = Q[1] + z * Q[2] + z2 * (Q[3] + z * Q[4]) + z4 * (Q[5] + z * Q[6])
        sgn*y*p/q
    end
end


text:
# Trace utilities

function get_scene_graph(trace)
    get_retval(trace).scene_graph
end
function structure_addr()
    (:scene_graph => :structure)
end
function obs_addr()
    :obs
end
function get_structure(trace)
    trace[structure_addr()]
end
function get_gen_cloud(trace)
    rendered_clouds = get_retval(trace).rendered_clouds
    length(rendered_clouds) == 1 || error(
        "rendered_clouds does not have length 1." *
        " Did you mean to call get_gen_cloud(trace, i)?")
    only(rendered_clouds)
end
function get_gen_cloud(trace, i)
    get_retval(trace).rendered_clouds[i]
end
function get_obs_cloud(trace)
    trace[:obs]
end
function get_object_boxes(trace)
    get_args(trace)[1].boxes
end
function get_object_meshes(trace)
    get_args(trace)[1].meshes
end
function get_num_objects(trace)
    length(get_object_boxes(trace))
end
function get_cam_pose(trace)::Pose
    get_args(trace)[1].camera_pose
end
function get_poses(trace)
    floatingPosesOf(get_scene_graph(trace))
end
function obj_name_from_idx(idx)
    Symbol("obj_", idx)
end

function get_edges(trace::Gen.Trace)
    collect(LG.edges(get_structure(trace)))
end

export scene, get_scene_graph, structure_addr, obs_addr, get_structure, get_gen_cloud, get_obs_cloud, get_object_names,
        get_object_boxes, get_num_objects, get_poses, get_cam_pose, get_edges, obj_name_from_idx, get_object_meshes

function get_unexplained_obs_cloud(trace, radius=1.0)
    obs = get_obs_cloud(trace)
    obs_tree = NearestNeighbors.KDTree(obs)

    valid = fill(true, size(obs)[2])

    idxs = NearestNeighbors.inrange(obs_tree, get_gen_cloud(trace), radius)
    valid[unique(vcat(idxs...))] .= false
    valid
end

export get_unexplained_obs_cloud

floating_pose_addr(i::Int) = (:scene_graph => :params => i => :floating => :pose)
floating_addr(i::Int) = (:scene_graph => :params => i => :floating)
contact_addr(i::Int) = (:scene_graph => :params => i => :contact)
contact_addr(i::Int, addr::Symbol) = (:scene_graph => :params => i => :contact => addr)

export floating_pose_addr, contact_addr


text:
using RoME
using Test
# using RoMEPlotting, Distributions

import IncrementalInference: getSample

mutable struct NorthSouthPartial{T} <: AbstractPrior
  Z::T
  partial::Tuple{Int}
  NorthSouthPartial{D}() where D = new{D}()
  NorthSouthPartial{D}(Z::D) where {D <: IIF.SamplableBelief} = new{D}(Z, (2,))
end
NorthSouthPartial(Z::D) where {D <: IIF.SamplableBelief} = NorthSouthPartial{D}(Z)

getSample(cfo::CalcFactor{<:NorthSouthPartial}, N::Int=1) = (reshape(rand(cfo.factor.Z, N),1,N),)

@testset "test standard multimodal range bearing factor setup..." begin

# Start with an empty graph
global N = 100
global fg = initfg(sessionname="MULTIMODAL_2D_TUTORIAL")

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
addVariable!(fg, :l2, Point2, tags=[:LANDMARK;])

addFactor!(fg, [:l1], Prior(MvNormal([10.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) )
addFactor!(fg, [:l2], Prior(MvNormal([30.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) )

val, = predictbelief(fg, :l1, [:l1f1;])
setVal!(fg, :l1, val)

val, = predictbelief(fg, :l2, [:l2f1;])
setVal!(fg, :l2, val)

addVariable!(fg, :x0, Pose2)
# addFactor!(fg, [:x0], Prior(MvNormal([0.0;0.0;0], Matrix(Diagonal([1.0;1.0;0.01].^2)))) )
addFactor!(fg, [:x0;], NorthSouthPartial(Normal(0,1.0)))

global p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1; :l2], p2br, multihypo=[1.0; 0.5; 0.5])

predictbelief(fg, :x0, ls(fg, :x0))

# writeGraphPdf(fg)
end

@testset "test multimodal bearing range factors calculate pose position properly..." begin

ensureAllInitialized!(fg)

tree, smt, hist = solveTree!(fg)
# global tree = batchSolve!(fg, N=N, recursive=true)
# global tree = wipeBuildNewTree!(fg)
# global infv = [inferOverTreeR!(fg, tree, N=N) for i in 1:4]
# @test length(infv) == 4

# X0 = getKDE(getVariable(fg, :x0))
# X0pts = getPoints(X0)
global X0pts = getKDE(getVariable(fg, :x0)) |> getPoints

@test size(X0pts,1) == 3
@test size(X0pts,2) == N

# find at least two of the four possible modes
@test sum([0 < sum(-20 .< X0pts[1,:] .< 0);
           0 < sum(  0 .< X0pts[1,:] .< 20);
           0 < sum( 20 .< X0pts[1,:] .< 40);
           0 < sum( 40 .< X0pts[1,:] .< 60)]) >= 2

end

@testset "test multimodal landmark locations are computed correclty..." begin

# Start with an empty graph
global fg = initfg(sessionname="MULTIMODAL_2D_TUTORIAL")

addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal([0.0;0.0;0], Matrix(Diagonal([1.0;1.0;0.01].^2)))) ) # TODO IIF.Prior with IIF 0.3.9

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
addFactor!(fg, [:l1], PriorPose2(MvNormal([40.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) ) # TODO IIF.Prior with IIF 0.3.9

addVariable!(fg, :l2, Point2, tags=[:LANDMARK;])
addFactor!(fg, [:l2;], NorthSouthPartial(Normal(0,1.0)))
# addFactor!(fg, [:l2], PriorPose2(MvNormal([30.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) ) # TODO IIF.Prior with IIF 0.3.9

global p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1; :l2], p2br, multihypo=[1.0; 0.5; 0.5])

# @test true

tree, smt, hist = solveTree!(fg)
# setVal!(getVert(fg, :l2), zeros(2,N))
# writeGraphPdf(fg)
# global tree = wipeBuildNewTree!(fg)
# global infv = [inferOverTreeR!(fg, tree, N=N) for i in 1:4]
# @test length(infv) == 4

# L1 = getKDE(getVariable(fg, :l1))
global L2 = getKDE(getVariable(fg, :l2))
# X0 = getKDE(getVariable(fg, :x0))

global L2pts = getPoints(L2)

@test size(L2pts,1) == 2
@test size(L2pts,2) == N

# plot(L2, dims=[1;2])

numM1 = sum(10 .< L2pts[1,:] .< 30)
numM2 = sum(30 .< L2pts[1,:] .< 50)

@test 20 < numM1 < 70
@test 20 < numM2 < 70
@test 0.6*N <= numM1 + numM2

end

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,20cm)
#
# drawPosesLandms(fg, spscale = 1.0)
#
# plotLocalProduct(fg, :l1, levels=10)
# plotLocalProduct(fg, :l2, levels=10)