text:
type Scene
    vehicles::Vector{Vehicle} # this is a pre-allocated array that is at least as large as the maximum number of vehicles in a Trajdata frame
    n_vehicles::Int

    function Scene(n_vehicles::Int=500)
        vehicles = Array(Vehicle, n_vehicles)
        for i in 1 : length(vehicles)
            vehicles[i] = Vehicle()
        end
        new(vehicles, 0)
    end
    function Scene(
        vehicles::Vector{Vehicle},
        n_vehicles::Int=length(vehicles);
        )

        new(vehicles, n_vehicles)
    end
end

# iteration
Base.start(scene::Scene) = 1
Base.done(scene::Scene, i::Int) = i > length(scene)
Base.next(scene::Scene, i::Int) = (scene.vehicles[i], i+1)

# copying
function Base.copy!(dest::Scene, src::Scene)
    for i in 1 : src.n_vehicles
        copy!(dest.vehicles[i], src.vehicles[i])
    end
    dest.n_vehicles = src.n_vehicles
    dest
end

function Base.push!(scene::Scene, veh::Vehicle)
    scene.n_vehicles += 1
    scene.vehicles[scene.n_vehicles] = veh
    scene
end

function get_first_available_id(scene::Scene)
    ids = Set{Int}([veh.def.id for veh in scene])
    id = 1
    while in(id, ids)
        id += 1
    end
    return id
end
Base.push!(scene::Scene, state::VehicleState) = push!(scene, Vehicle(state, VehicleDef(get_first_available_id(scene))))

Base.length(scene::Scene) = scene.n_vehicles
Base.getindex(scene::Scene, i::Int) = scene.vehicles[i]
function Base.setindex!(scene::Scene, veh::Vehicle, i::Int)
    scene.vehicles[i] = veh
    scene
end
function Base.empty!(scene::Scene)
    scene.n_vehicles = 0
    scene
end
function Base.get!(scene::Scene, trajdata::Trajdata, frame::Int)

    scene.n_vehicles = 0

    if frame_inbounds(trajdata, frame)
        tdframe = trajdata.frames[frame]
        for i in tdframe.lo : tdframe.hi
            scene.n_vehicles += 1
            veh = scene.vehicles[scene.n_vehicles]
            s = trajdata.states[i]
            veh.state = s.state
            veh.def = get_vehicledef(trajdata, s.id)
        end
    end

    scene
end
function Base.get!(scene::Scene, trajdata::Trajdata, frame_lo::Int, frame_hi::Int, γ::Float64)
    # add all vehicles that exist in both the lower and upper frames

    if frame_hi != frame_lo+1 && frame_hi != frame_lo
        # NOTE: frame_hi should pretty much always be frame_lo + 1
        #       the most frequent exception would be frame_lo = frame_hi
        #       you otherwise pretty much never want another case
        warn("get!(scene, trajdata, frame_lo, frame_hi, γ) used with frame_hi = frame_lo + $(frame_hi - frame_lo)!")
    end

    scene.n_vehicles = 0

    tdframe_lo = trajdata.frames[frame_lo]
    tdframe_hi = trajdata.frames[frame_hi]

    for i in tdframe_lo.lo : tdframe_lo.hi
        s = trajdata.states[i]
        for j in tdframe_hi.lo : tdframe_hi.hi
            s2 = trajdata.states[j]
            if s2.id == s.id # we are good
                # add it
                scene.n_vehicles += 1

                state_lo = s.state
                state_hi = s2.state

                veh = scene.vehicles[scene.n_vehicles]
                veh.state = lerp(state_lo, state_hi, γ, trajdata.roadway)
                veh.def = get_vehicledef(trajdata, s.id)
                break
            end
        end
    end

    scene
end
function Base.get!(scene::Scene, trajdata::Trajdata, time::Float64)

    # performs linear interpolation between frames

    i = searchsortedfirst(trajdata.frames, time, by=frame_or_time->isa(frame_or_time, TrajdataFrame) ? frame_or_time.t : frame_or_time) # index of first frame ≥ t_current

    if i == 1
        lo, hi = 1, 2
    elseif i ≤ nframes(trajdata)
        lo, hi = i-1, i
    else # t_current > all times in trajdata.frames
        lo, hi = i-2, i-1
    end

    t_lo = get_time(trajdata, lo)
    t_hi = get_time(trajdata, hi)
    γ = (time - t_lo) / (t_hi - t_lo)

    if isapprox(γ, 0.0)
        get!(scene, trajdata, lo)
    elseif isapprox(γ, 1.0)
        get!(scene, trajdata, hi)
    else
        get!(scene, trajdata, lo, hi, γ)
    end

    scene
end
function Base.deleteat!(scene::Scene, vehicle_index::Int)
    for i in vehicle_index : scene.n_vehicles - 1
        copy!(scene.vehicles[i], scene.vehicles[i+1])
    end
    scene.n_vehicles -= 1
    scene
end
Base.delete!(scene::Scene, veh::Vehicle) = deleteat!(scene, get_index_of_first_vehicle_with_id(scene, veh.def.id))

function get_by_id(scene::Scene, id::Int)
    for i in 1 : scene.n_vehicles
        if scene.vehicles[i].def.id == id
            return scene.vehicles[i]
        end
    end
    error("vehicle with id $id not found in scene")
end
function get_index_of_first_vehicle_with_id(scene::Scene, id::Int)
    retval = 0
    for i in 1 : scene.n_vehicles
        if scene.vehicles[i].def.id == id
            retval = i
            break
        end
    end
    retval
end
iscarinframe(scene::Scene, id::Int) = get_index_of_first_vehicle_with_id(scene, id) != 0

function get_vehiclestate(scene::Scene, id::Int)
    for i in 1 : scene.n_vehicles
        if scene.vehicles[i].def.id == id
            return scene.vehicles[i].state
        end
    end
    error("get_vehiclestate: vehicle with id $id not found")
end
function get_vehicle!(veh::Vehicle, scene::Scene, id::Int)
    for i in 1 : scene.n_vehicles
        if scene.vehicles[i].def.id == id
            copy!(veh, scene.vehicles[i])
            return veh
        end
    end
    error("get_vehicle!: vehicle with id $id not found")
end
get_vehicle(scene::Scene, id::Int) = get_vehicle!(Vehicle(), scene, id)

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

"""
    get_neighbor_index_fore(scene::Scene, vehicle_index::Int, roadway::Roadway)
Return the index of the vehicle that is in the same lane as scene[vehicle_index] and
in front of it with the smallest distance along the lane

    The method will search on the current lane first, and if no vehicle is found it
    will continue to travel along the lane following next_lane(lane, roadway).
    If no vehicle is found within `max_distance_fore,` a value of 0 is returned instead.
"""
immutable NeighborLongitudinalResult
    ind::Int # index in scene of the neighbor
    Δs::Float64 # positive distance along lane between vehicles' positions
end

export
    VehicleTargetPoint,
    VehicleTargetPointFront,
    VehicleTargetPointCenter,
    VehicleTargetPointRear,
    get_targetpoint_delta

abstract VehicleTargetPoint
immutable VehicleTargetPointFront <: VehicleTargetPoint end
get_targetpoint_delta(::VehicleTargetPointFront, veh::Vehicle) = veh.def.length/2*cos(veh.state.posF.ϕ)
immutable VehicleTargetPointCenter <: VehicleTargetPoint end
get_targetpoint_delta(::VehicleTargetPointCenter, veh::Vehicle) = 0.0
immutable VehicleTargetPointRear <: VehicleTargetPoint end
get_targetpoint_delta(::VehicleTargetPointRear, veh::Vehicle) = -veh.def.length/2*cos(veh.state.posF.ϕ)

const VEHICLE_TARGET_POINT_CENTER = VehicleTargetPointCenter()

function get_neighbor_fore_along_lane(
    scene::Scene,
    roadway::Roadway,
    tag_start::LaneTag,
    s_base::Float64,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
    index_to_ignore::Int=-1,
    )

    best_ind = 0
    best_dist = max_distance_fore
    tag_target = tag_start

    dist_searched = 0.0
    while dist_searched < max_distance_fore

        lane = roadway[tag_target]

        for (i,veh) in enumerate(scene)
            if i != index_to_ignore

                s_adjust = NaN
                if veh.state.posF.roadind.tag == tag_target
                    s_adjust = 0.0

                elseif is_between_segments_hi(veh.state.posF.roadind.ind, lane.curve) &&
                       is_in_entrances(roadway[tag_target], veh.state.posF.roadind.tag)

                    distance_between_lanes = abs(roadway[tag_target].curve[1].pos - roadway[veh.state.posF.roadind.tag].curve[end].pos)
                    s_adjust = -(roadway[veh.state.posF.roadind.tag].curve[end].s + distance_between_lanes)

                elseif is_between_segments_lo(veh.state.posF.roadind.ind) &&
                       is_in_exits(roadway[tag_target], veh.state.posF.roadind.tag)

                    distance_between_lanes = abs(roadway[tag_target].curve[end].pos - roadway[veh.state.posF.roadind.tag].curve[1].pos)
                    s_adjust = roadway[tag_target].curve[end].s + distance_between_lanes
                end

                if !isnan(s_adjust)
                    s_valid = veh.state.posF.s + get_targetpoint_delta(targetpoint_valid, veh) + s_adjust
                    dist_valid = s_valid - s_base + dist_searched
                    if dist_valid ≥ 0.0
                        s_primary = veh.state.posF.s + get_targetpoint_delta(targetpoint_primary, veh) + s_adjust
                        dist = s_primary - s_base + dist_searched
                        if dist < best_dist
                            best_dist = dist
                            best_ind = i
                        end
                    end
                end
            end
        end

        if best_ind != 0
            break
        end

        if !has_next(lane) ||
           (tag_target == tag_start && dist_searched != 0.0) # exit after visiting this lane a 2nd time
            break
        end

        dist_searched += (lane.curve[end].s - s_base)
        s_base = -abs(lane.curve[end].pos - next_lane_point(lane, roadway).pos) # negative distance between lanes
        tag_target = next_lane(lane, roadway).tag
    end

    NeighborLongitudinalResult(best_ind, best_dist)
end
function get_neighbor_fore_along_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway,
    targetpoint_ego::VehicleTargetPoint,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_fore::Float64 = 250.0 # max distance to search forward [m]
    )

    veh_ego = scene[vehicle_index]
    tag_start = veh_ego.state.posF.roadind.tag
    s_base = veh_ego.state.posF.s + get_targetpoint_delta(targetpoint_ego, veh_ego)

    get_neighbor_fore_along_lane(scene, roadway, tag_start, s_base,
        targetpoint_primary, targetpoint_valid,
        max_distance_fore=max_distance_fore, index_to_ignore=vehicle_index)
end
function get_neighbor_fore_along_left_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway,
    targetpoint_ego::VehicleTargetPoint,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_fore::Float64 = 250.0 # max distance to search forward [m]
    )

    retval = NeighborLongitudinalResult(0, max_distance_fore)

    veh_ego = scene[vehicle_index]
    lane = roadway[veh_ego.state.posF.roadind.tag]
    if n_lanes_left(lane, roadway) > 0
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        roadproj = proj(veh_ego.state.posG, lane_left, roadway)
        tag_start = roadproj.tag
        s_base = lane_left[roadproj.curveproj.ind, roadway].s + get_targetpoint_delta(targetpoint_ego, veh_ego)

        retval = get_neighbor_fore_along_lane(scene, roadway, tag_start, s_base,
                                              targetpoint_primary, targetpoint_valid,
                                              index_to_ignore=vehicle_index,
                                              max_distance_fore=max_distance_fore)
    end

    retval
end
function get_neighbor_fore_along_right_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway,
    targetpoint_ego::VehicleTargetPoint,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_fore::Float64 = 250.0 # max distance to search forward [m]
    )

    retval = NeighborLongitudinalResult(0, max_distance_fore)

    veh_ego = scene[vehicle_index]
    lane = roadway[veh_ego.state.posF.roadind.tag]
    if n_lanes_right(lane, roadway) > 0
        lane_right = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        roadproj = proj(veh_ego.state.posG, lane_right, roadway)
        tag_start = roadproj.tag
        s_base = lane_right[roadproj.curveproj.ind, roadway].s + get_targetpoint_delta(targetpoint_ego, veh_ego)

        retval = get_neighbor_fore_along_lane(scene, roadway, tag_start, s_base,
                                              targetpoint_primary, targetpoint_valid,
                                              index_to_ignore=vehicle_index,
                                              max_distance_fore=max_distance_fore)
    end

    retval
end

function get_neighbor_fore_along_lane(
    scene::Scene,
    roadway::Roadway,
    tag_start::LaneTag,
    s_base::Float64;
    max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
    index_to_ignore::Int=-1,
    )

    get_neighbor_fore_along_lane(scene, roadway, tag_start, s_base,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        max_distance_fore=max_distance_fore,
        index_to_ignore=index_to_ignore)
end
function get_neighbor_fore_along_lane(scene::Scene, vehicle_index::Int, roadway::Roadway;
    max_distance_fore::Float64 = 250.0 # max distance to search forward [m]
    )

    get_neighbor_fore_along_lane(scene, vehicle_index, roadway,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        VEHICLE_TARGET_POINT_CENTER, max_distance_fore=max_distance_fore)
end
function get_neighbor_fore_along_left_lane(scene::Scene, vehicle_index::Int, roadway::Roadway;
    max_distance_fore::Float64 = 250.0 # max distance to search forward [m]
    )

    get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        VEHICLE_TARGET_POINT_CENTER, max_distance_fore=max_distance_fore)
end
function get_neighbor_fore_along_right_lane(scene::Scene, vehicle_index::Int, roadway::Roadway;
    max_distance_fore::Float64 = 250.0 # max distance to search forward [m]
    )

    get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        VEHICLE_TARGET_POINT_CENTER, max_distance_fore=max_distance_fore)
end

function get_neighbor_rear_along_lane(
    scene::Scene,
    roadway::Roadway,
    tag_start::LaneTag,
    s_base::Float64,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_rear::Float64 = 250.0, # max distance to search rearward [m]
    index_to_ignore::Int=-1,
    )

    best_ind = 0
    best_dist = max_distance_rear
    tag_target = tag_start

    dist_searched = 0.0
    while dist_searched < max_distance_rear

        lane = roadway[tag_target]

        for (i,veh) in enumerate(scene)
            if i != index_to_ignore

                s_adjust = NaN

                if veh.state.posF.roadind.tag == tag_target
                    s_adjust = 0.0

                elseif is_between_segments_hi(veh.state.posF.roadind.ind, lane.curve) &&
                       is_in_entrances(roadway[tag_target], veh.state.posF.roadind.tag)

                    distance_between_lanes = abs(roadway[tag_target].curve[1].pos - roadway[veh.state.posF.roadind.tag].curve[end].pos)
                    s_adjust = -(roadway[veh.state.posF.roadind.tag].curve[end].s + distance_between_lanes)

                elseif is_between_segments_lo(veh.state.posF.roadind.ind) &&
                       is_in_exits(roadway[tag_target], veh.state.posF.roadind.tag)

                    distance_between_lanes = abs(roadway[tag_target].curve[end].pos - roadway[veh.state.posF.roadind.tag].curve[1].pos)
                    s_adjust = roadway[tag_target].curve[end].s + distance_between_lanes
                end

                if !isnan(s_adjust)
                    s_valid = veh.state.posF.s + get_targetpoint_delta(targetpoint_valid, veh) + s_adjust
                    dist_valid = s_base - s_valid + dist_searched
                    if dist_valid ≥ 0.0
                        s_primary = veh.state.posF.s + get_targetpoint_delta(targetpoint_primary, veh) + s_adjust
                        dist = s_base - s_primary + dist_searched
                        if dist < best_dist
                            best_dist = dist
                            best_ind = i
                        end
                    end
                end
            end
        end

        if best_ind != 0
            break
        end

        if !has_prev(lane) ||
           (tag_target == tag_start && dist_searched != 0.0) # exit after visiting this lane a 2nd time
            break
        end

        dist_searched += s_base
        s_base = lane.curve[end].s + abs(lane.curve[end].pos - prev_lane_point(lane, roadway).pos) # end of prev lane plus crossover
        tag_target = prev_lane(lane, roadway).tag
    end

    NeighborLongitudinalResult(best_ind, best_dist)
end
function get_neighbor_rear_along_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway,
    targetpoint_ego::VehicleTargetPoint,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_rear::Float64 = 250.0 # max distance to search rearward [m]
    )

    veh_ego = scene[vehicle_index]
    tag_start = veh_ego.state.posF.roadind.tag
    s_base = veh_ego.state.posF.s + get_targetpoint_delta(targetpoint_ego, veh_ego)

    get_neighbor_rear_along_lane(scene, roadway, tag_start, s_base,
        targetpoint_primary, targetpoint_valid,
        max_distance_rear=max_distance_rear, index_to_ignore=vehicle_index)
end
function get_neighbor_rear_along_left_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway,
    targetpoint_ego::VehicleTargetPoint,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_rear::Float64 = 250.0 # max distance to search forward [m]
    )

    retval = NeighborLongitudinalResult(0, max_distance_rear)

    veh_ego = scene[vehicle_index]
    lane = roadway[veh_ego.state.posF.roadind.tag]
    if n_lanes_left(lane, roadway) > 0
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        roadproj = proj(veh_ego.state.posG, lane_left, roadway)
        tag_start = roadproj.tag
        s_base = lane_left[roadproj.curveproj.ind, roadway].s + get_targetpoint_delta(targetpoint_ego, veh_ego)

        retval = get_neighbor_rear_along_lane(scene, roadway, tag_start, s_base,
                                              targetpoint_primary, targetpoint_valid,
                                              index_to_ignore=vehicle_index,
                                              max_distance_rear=max_distance_rear)
    end

    retval
end
function get_neighbor_rear_along_right_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway,
    targetpoint_ego::VehicleTargetPoint,
    targetpoint_primary::VehicleTargetPoint, # the reference point whose distance we want to minimize
    targetpoint_valid::VehicleTargetPoint; # the reference point, which if distance to is positive, we include the vehicle
    max_distance_rear::Float64 = 250.0 # max distance to search forward [m]
    )

    retval = NeighborLongitudinalResult(0, max_distance_rear)

    veh_ego = scene[vehicle_index]
    lane = roadway[veh_ego.state.posF.roadind.tag]
    if n_lanes_right(lane, roadway) > 0
        lane_right = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        roadproj = proj(veh_ego.state.posG, lane_right, roadway)
        tag_start = roadproj.tag
        s_base = lane_right[roadproj.curveproj.ind, roadway].s + get_targetpoint_delta(targetpoint_ego, veh_ego)

        retval = get_neighbor_rear_along_lane(scene, roadway, tag_start, s_base,
                                              targetpoint_primary, targetpoint_valid,
                                              index_to_ignore=vehicle_index,
                                              max_distance_rear=max_distance_rear)
    end

    retval
end

function get_neighbor_rear_along_lane(
    scene::Scene,
    roadway::Roadway,
    tag_start::LaneTag,
    s_base::Float64;
    max_distance_rear::Float64 = 250.0, # max distance to search rearward [m]
    index_to_ignore::Int=-1,
    )

    get_neighbor_rear_along_lane(scene, roadway, tag_start, s_base,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        max_distance_rear=max_distance_rear,
        index_to_ignore=index_to_ignore)
end
function get_neighbor_rear_along_lane(scene::Scene, vehicle_index::Int, roadway::Roadway;
    max_distance_rear::Float64 = 250.0 # max distance to search forward [m]
    )

    get_neighbor_rear_along_lane(scene, vehicle_index, roadway,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        VEHICLE_TARGET_POINT_CENTER, max_distance_rear=max_distance_rear)
end
function get_neighbor_rear_along_left_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway;
    max_distance_rear::Float64 = 250.0 # max distance to search forward [m]
    )

    get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        VEHICLE_TARGET_POINT_CENTER, max_distance_rear=max_distance_rear)
end
function get_neighbor_rear_along_right_lane(
    scene::Scene,
    vehicle_index::Int,
    roadway::Roadway;
    max_distance_rear::Float64 = 250.0 # max distance to search forward [m]
    )

    get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway,
        VEHICLE_TARGET_POINT_CENTER, VEHICLE_TARGET_POINT_CENTER,
        VEHICLE_TARGET_POINT_CENTER, max_distance_rear=max_distance_rear)
end

"""
    Project the given point to the same lane as the given RoadIndex.
This will return the projection of the point, along with the Δs along the lane from the RoadIndex.
"""
immutable FrenetRelativePosition
    origin::RoadIndex
    target::RoadIndex
    Δs::Float64
    t::Float64
    ϕ::Float64
end
function get_frenet_relative_position(posG::VecSE2, roadind::RoadIndex, roadway::Roadway;
    max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
    max_distance_rear::Float64 = 250.0, # max distance to search backward [m]
    improvement_threshold::Float64 = 1e-4,
    )

    # project to current lane first
    tag_start = roadind.tag
    lane_start = roadway[tag_start]
    curveproj_start = proj(posG, lane_start, roadway, move_along_curves=false).curveproj
    curvept_start = lane_start[curveproj_start.ind, roadway]
    s_base = lane_start[roadind.ind, roadway].s
    s_proj = curvept_start.s
    Δs = s_proj - s_base
    sq_dist_to_curve = abs2(posG - curvept_start.pos)
    retval = FrenetRelativePosition(roadind,
                RoadIndex(curveproj_start.ind, tag_start), Δs, curveproj_start.t, curveproj_start.ϕ)

    # search downstream
    if has_next(lane_start)
        dist_searched = lane_start.curve[end].s - s_base
        s_base = -abs(lane_start.curve[end].pos - next_lane_point(lane_start, roadway).pos) # negative distance between lanes
        tag_target = next_lane(lane_start, roadway).tag
        while dist_searched < max_distance_fore

            lane = roadway[tag_target]
            curveproj = proj(posG, lane, roadway, move_along_curves=false).curveproj
            sq_dist_to_curve2 = abs2(posG - lane[curveproj.ind, roadway].pos)

            if sq_dist_to_curve2 < sq_dist_to_curve - improvement_threshold

                sq_dist_to_curve = sq_dist_to_curve2
                s_proj = lane[curveproj.ind, roadway].s
                Δs = s_proj - s_base + dist_searched
                retval = FrenetRelativePosition(
                            roadind, RoadIndex(curveproj.ind, tag_target),
                            Δs, curveproj.t, curveproj.ϕ)
            end

            if !has_next(lane)
                break
            end

            dist_searched += (lane.curve[end].s - s_base)
            s_base = -abs(lane.curve[end].pos - next_lane_point(lane, roadway).pos) # negative distance between lanes
            tag_target = next_lane(lane, roadway).tag

            if tag_target == tag_start
                break
            end
        end
    end

    # search upstream
    if has_prev(lane_start)
        dist_searched = s_base
        tag_target = prev_lane(lane_start, roadway).tag
        s_base = roadway[tag_target].curve[end].s + abs(lane_start.curve[1].pos - prev_lane_point(lane_start, roadway).pos) # end of the lane

        while dist_searched < max_distance_fore

            lane = roadway[tag_target]
            curveproj = proj(posG, lane, roadway, move_along_curves=false).curveproj
            sq_dist_to_curve2 = abs2(posG - lane[curveproj.ind, roadway].pos)

            if sq_dist_to_curve2 < sq_dist_to_curve - improvement_threshold

                sq_dist_to_curve = sq_dist_to_curve2
                s_proj = lane[curveproj.ind, roadway].s
                dist = s_base - s_proj  + dist_searched
                Δs = -dist
                retval = FrenetRelativePosition(
                            roadind, RoadIndex(curveproj.ind, tag_target),
                            Δs, curveproj.t, curveproj.ϕ)
            end

            if !has_prev(lane)
                break
            end

            dist_searched += s_base
            s_base = lane.curve[end].s + abs(lane.curve[1].pos - prev_lane_point(lane, roadway).pos) # length of this lane plus crossover
            tag_target = prev_lane(lane, roadway).tag

            if tag_target == tag_start
                break
            end
        end
    end

    retval
end
get_frenet_relative_position(veh_fore::Vehicle, veh_rear::Vehicle, roadway::Roadway) = get_frenet_relative_position(veh_fore.state.posG, veh_rear.state.posF.roadind, roadway)


text:
using AbInitioTopo
using Test

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

ext: jl


text:
module KG_3_1

using Jecco
using LinearAlgebra
using OrdinaryDiffEq

import Base.Threads.@threads
import Base.Threads.@spawn

# abstract types and structs used throughout
include("types.jl")

include("system.jl")

include("potential.jl")

include("initial_data.jl")

# nested system
include("equation_coeff.jl")
include("solve_nested.jl")

# evolution equations
include("compute_bulkevolved_t.jl")

# time marching orders
include("setup_rhs.jl")

# input/output
include("IO.jl")

# run the model
include("run.jl")

# always set the number of BLAS threads to 1 upon loading the module. by default
# it uses a bunch of them and we don't want that since they trample over each
# other when solving the nested systems equations. it's much better to thread
# over the loop. see also the discussion here:
# https://github.com/JuliaLang/julia/issues/33409
#
# this saves us the tedious task of always setting OMP_NUM_THREADS=1 before
# launching julia.
function __init__()
    LinearAlgebra.BLAS.set_num_threads(1)
    nothing
end

export SpecCartGrid3D
export Potential, ConstPotential, SquarePotential
export System, SystemPartition
export BulkEvolvedPartition, BulkConstrainedPartition, Boundary
export Uniform2D, Sine2D
export AffineNull
export Integration, InOut, run_model

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ResourceMetricStatus indicates the current value of a resource metric known to Kubernetes, as specified in requests and limits, describing each pod in the current scale target (e.g. CPU or memory).  Such metrics are built in to Kubernetes, and have special scaling options on top of those available to normal per-pod metrics using the \"pods\" source.

    IoK8sApiAutoscalingV2beta1ResourceMetricStatus(;
        currentAverageUtilization=nothing,
        currentAverageValue=nothing,
        name=nothing,
    )

    - currentAverageUtilization::Int32 : currentAverageUtilization is the current value of the average of the resource metric across all relevant pods, represented as a percentage of the requested value of the resource for the pods.  It will only be present if `targetAverageValue` was set in the corresponding metric specification.
    - currentAverageValue::IoK8sApimachineryPkgApiResourceQuantity : currentAverageValue is the current value of the average of the resource metric across all relevant pods, as a raw value (instead of as a percentage of the request), similar to the \"pods\" metric source type. It will always be set, regardless of the corresponding metric specification.
    - name::String : name is the name of the resource in question.
"""
mutable struct IoK8sApiAutoscalingV2beta1ResourceMetricStatus <: SwaggerModel
    currentAverageUtilization::Any # spec type: Union{ Nothing, Int32 } # spec name: currentAverageUtilization
    currentAverageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: currentAverageValue
    name::Any # spec type: Union{ Nothing, String } # spec name: name

    function IoK8sApiAutoscalingV2beta1ResourceMetricStatus(;currentAverageUtilization=nothing, currentAverageValue=nothing, name=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricStatus, Symbol("currentAverageUtilization"), currentAverageUtilization)
        setfield!(o, Symbol("currentAverageUtilization"), currentAverageUtilization)
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricStatus, Symbol("currentAverageValue"), currentAverageValue)
        setfield!(o, Symbol("currentAverageValue"), currentAverageValue)
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricStatus, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        o
    end
end # type IoK8sApiAutoscalingV2beta1ResourceMetricStatus

const _property_map_IoK8sApiAutoscalingV2beta1ResourceMetricStatus = Dict{Symbol,Symbol}(Symbol("currentAverageUtilization")=>Symbol("currentAverageUtilization"), Symbol("currentAverageValue")=>Symbol("currentAverageValue"), Symbol("name")=>Symbol("name"))
const _property_types_IoK8sApiAutoscalingV2beta1ResourceMetricStatus = Dict{Symbol,String}(Symbol("currentAverageUtilization")=>"Int32", Symbol("currentAverageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("name")=>"String")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1ResourceMetricStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1ResourceMetricStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1ResourceMetricStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1ResourceMetricStatus)
    (getproperty(o, Symbol("currentAverageValue")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricStatus }, name::Symbol, val)
end


text:
using Distributions

struct FiniteSimulation{F<:Number,W<:World, P<:Plotter} <: Simulation
    time::Base.RefValue{F}
    t_sheep_boredom::Base.RefValue{F}
    dt::F
    end_time::F
    world::W
    plotter::P
end

FiniteSimulation(;
    time = 0.0,
    dt::F,
    end_time = 10.0,
    world::W,
    t_sheep_boredom = rand(Exponential(world.meanSheepDiffusionTime / world.sheeps.current_sheep[])),
    plotter::P = VoidPlotter()
    ) where {F<:Number,W<:World,P<:Plotter} = FiniteSimulation{F,W,P}( Ref(time), Ref(t_sheep_boredom), dt, end_time, world, plotter )

function runsim( sim::FiniteSimulation )
    p = plot( sim.plotter, sim.world, sim.time[] )
    io = IOBuffer()
    print(io, "\e[25l")
    # try
        for _ in 0:sim.end_time÷sim.dt
            iterate!( sim )

            plot!( io, p, sim.plotter, sim.world, sim.time[] )
        end # while
    # catch err
    #     @error "Caught $err !"
    #     rethrow( err )
    # finally
        print(io, "\e[25h")
        endtime = sim.time[]
        sim.time[] = zero(sim.time[])
        return endtime
    # end # try
end # function


text:
# TODO: introduce a way for solvers to indicate that they support some features. For now, only the barest FZN file is generated, with extremely little structure.

# A few generic helpers for parsing the FZN output.

mutable struct _FznResults
    raw_status_string::String
    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    objective_value::Real # TODO: implement this. Compute it from the solutions? Fzn doesn't really output it (unless the FZN file is tweaked).
    primal_solutions::Vector{Dict{MOI.VariableIndex, Real}} # Solution ID -> (variable index -> value)
end

function _FznResults()
    return _FznResults(
        "Optimize not called.",
        MOI.OPTIMIZE_NOT_CALLED,
        MOI.NO_SOLUTION,
        NaN,
        Dict{MOI.VariableIndex, Float64}[],
    )
end

function _parse_fzn_value(str::AbstractString)
    # Heuristically guess the type of the output value: either integer or
    # float.
    if '.' in str
        return parse(Float64, str)
    else
        return parse(Int, str)
    end
end

"""
    _parse_to_assignments(str::String)::Vector{Dict{String, Vector{Number}}}

Parses the output of a FlatZinc-compatible solver into a list of dictionaries
mapping the name of the variables to their values (either a scalar or a vector
of numbers). The values are automatically transformed into the closest type
(integer or float).
"""
function _parse_to_assignments(str::String)::Vector{Dict{String, Vector{Number}}}
    results = Dict{String, Vector{Number}}[]

    # If the infeasibility marker is discovered, return an empty list.
    if occursin("=====UNSATISFIABLE=====", str)
        return results
    end

    # Remove comments from the output (starting with % and ending with a new line).
    str = replace(str, r"%(.*)(\\r|\\n|\\r\\n)" => s"")

    # There may be several results returned by the solver. Each solution is
    # separated from the others by `'-' ^ 10`.
    str_split = split(str, '-' ^ 10)[1:(end - 1)]
    n_results = length(str_split)
    sizehint!(results, n_results)

    for i in 1:n_results
        push!(results, Dict{String, Vector{Number}}())

        # Each value is indicated in its own statement, separated by a
        # semi-colon.
        for part in split(strip(str_split[i]), ';')
            if isempty(part)
                continue
            end

            var, val = split(part, '=')
            var = strip(var)
            val = strip(val)

            # Either an array or a scalar. Always return an array for
            # simplicity. A scalar is simply an array with one element.
            if !occursin("array", val)
                # Scalar. Just a value: "1", "1.0".
                results[i][var] = [_parse_fzn_value(val)]
            else
                # Array. Several arguments: "array1d(1..2, [1, 2])",
                # "array2d(1..2, 1..2, [1, 2, 3, 4])".
                # TODO: should dimensions be preserved? (First argument[s] of arrayNd.)
                val = split(split(val, '[')[2], ']')[1]
                results[i][var] = map(_parse_fzn_value, map(strip, split(val, ',')))
            end
        end
    end

    return results
end

# MOI wrapper.
# Based on AmplNLWriter.jl's _NLResults and Optimizer.
# The main difference is that typical solutions do not have a Float64 type,
# but rather Int. However, it all depends on the actual FZN solver that is
# used below (some of them can still deal with floats).

mutable struct Optimizer <: MOI.AbstractOptimizer
    # Solver to call and options.
    inner::CP.FlatZinc.Model
    solver_command::Cmd
    time_limit_ms::Float64
    verboseness::Bool # True: verbose (default).
    options::Vector{String} # Custom-set options.

    # Results once solver called.
    results::_FznResults
    solve_time::Float64 # For the call to optimize!
    fzn_time::Float64 # For writing the FZN file
    parse_time::Float64 # For reading the FZN solution
end

"""
    Optimizer(
        solver_command::Union{String, Function},
        solver_args::Vector{String},
    )

Create a new FlatZinc-backed Optimizer object. `solver_command` is the path
to the FlatZinc-compatible CLI of the solver or a `Cmd` object to call the
solver (typically, this is useful to set environment variables).

`solver_args` is a vector of `String` arguments passed solver executable.
However, prefer passing `key=value` options via `MOI.RawOptimizerAttribute`.
These arguments are passed to `Base.pipeline`.
[See the Julia documentation for more details](https://docs.julialang.org/en/v1/base/base/#Base.pipeline-Tuple{Base.AbstractCmd}).

## Examples

A string to an executable, no required argument:

```julia
Optimizer("/path/to/fzn.exe")

A `Cmd` object with an environment variable:

```julia
Optimizer(Cmd(`/path/to/fzn.exe`, env=["PATH=/usr/bin"]))
```
"""
function Optimizer(
    solver_command::Union{String, Cmd}="",
    solver_args::Vector{String}=String[],
)
    return Optimizer(
        CP.FlatZinc.Model(),
        `$(solver_command)`,
        0.0,
        true,
        solver_args,

        _FznResults(),
        NaN,
        NaN,
        NaN,
    )
end

Base.show(io::IO, ::Optimizer) = print(io, "A FlatZinc (flattened MiniZinc) model")

MOI.get(model::Optimizer, ::MOI.SolverName) = "FlatZincWriter"

function MOI.supports(model::Optimizer, attr::MOI.AnyAttribute, x...)
    return MOI.supports(model.inner, attr, x...)::Bool
end

function MOI.supports_add_constrained_variable(model::Optimizer, x::Type{S}) where {S <: MOI.AbstractScalarSet}
    return MOI.supports_add_constrained_variable(model.inner, x)::Bool
end

function MOI.supports_add_constrained_variables(model::Optimizer, x::Type{S}) where {S <: MOI.AbstractScalarSet}
    return MOI.supports_add_constrained_variables(model.inner, x)::Bool
end

function MOI.supports_constraint(model::Optimizer, f::Type{F}, s::Type{S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return MOI.supports_constraint(model.inner, f, s)::Bool
end

function MOI.get(model::Optimizer, attr::MOI.AnyAttribute, x...)
    return MOI.get(model.inner, attr, x...)
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    return MOI.get(model.inner, MOI.VariableIndex, name)
end

function MOI.set(model::Optimizer, attr::MOI.AnyAttribute, x...)
    MOI.set(model.inner, attr, x...)
    return
end

function MOI.add_variable(model::Optimizer)
    return MOI.add_variable(model.inner)
end

function MOI.add_constrained_variable(model::Optimizer, x::MOI.AbstractScalarSet)
    return MOI.add_constrained_variable(model.inner, x)
end

function MOI.add_constraint(model::Optimizer, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    return MOI.add_constraint(model.inner, f, s)
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    return MOIU.default_copy_to(dest, src)
end

function MOI.supports_incremental_interface(::Optimizer)
    return true
end

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.inner)
    # Only two attributes to empty, the other ones link the actual solver.
    model.results = _FznResults()
    model.solve_time = NaN
    return
end

MOI.is_empty(model::Optimizer) = MOI.is_empty(model.inner)

function MOI.optimize!(model::Optimizer)
    start_time = time()

    # Generate the FZN file.
    temp_dir = mktempdir()
    fzn_file = joinpath(temp_dir, "model.fzn")
    open(io -> write(io, model.inner), fzn_file, "w")

    model.fzn_time = time() - start_time

    # Generate the list of options. Always put the user-defined options at
    # the end.
    opts = copy(model.options)
    if !iszero(model.time_limit_ms)
        prepend!(opts, ["-t $(model.time_limit_ms)"])
    end
    if model.verboseness
        # Interpret as statistics to stdout and log to stderr.
        prepend!(opts, ["-s", "-v"])
    end

    # Call the FZN solver and gather the results in a string.
    try
        io = IOBuffer()
        ret = run(
            pipeline(
                `$(model.solver_command) $(opts) $(fzn_file)`,
                stdout=io,
            ),
            wait=true
        )

        if ret.exitcode != 0
            error("Nonzero exit code: $(ret.exitcode)")
        end

        start_parse_time = time()

        seekstart(io)
        sols_str = String(take!(io))
        sols_parsed = _parse_to_assignments(sols_str)

        model.parse_time = start_parse_time - time()

        _parse_to_moi_solutions(sols_parsed, model)
    catch err
        model.results = _FznResults(
            "Error calling the solver. Failed with: $(err)",
            MOI.OTHER_ERROR,
            MOI.NO_SOLUTION,
            NaN,
            Dict{MOI.VariableIndex, Float64}[],
        )
        # throw(err)
    end

    # Compute the total time of solving the MOI model.
    model.solve_time = time() - start_time
    return
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return model.results.termination_status
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if length(model.results.primal_solutions) >= attr.result_index
        return model.results.primal_solutions[attr.result_index][vi]
    else
        # No solution with this number. In particular, with infeasibility,
        # there is no solution, and this always returns nothing.
        return nothing
    end
end

MOI.supports(::Optimizer, ::MOI.SolveTimeSec) = true
MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.get(model::Optimizer, ::MOI.TimeLimitSec) = model.time_limit_ms / 1_000.0
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, limit::Real)
    model.time_limit_ms = limit * 1_000.0
end
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, limit::Nothing)
    model.time_limit_ms = 0.0
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.get(model::Optimizer, ::MOI.Silent) = !model.verboseness
function MOI.set(model::Optimizer, ::MOI.Silent, silentness::Bool)
    model.verboseness = !silentness
end

MOI.supports(::Optimizer, ::MOI.ResultCount) = true
MOI.get(model::Optimizer, ::MOI.ResultCount) = length(model.results.primal_solutions)

# Specific case of dual solution: getting it must be supported (MOI
# requirement), but few CP solvers have it accessible (none?).
# https://github.com/jump-dev/MathOptInterface.jl/pull/1561#pullrequestreview-740032701

MOI.supports(::Optimizer, ::MOI.DualStatus) = true
MOI.get(::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION

# Final helper to parse FZN output. It needs access to the model above...

"""
    _parse_to_assignments(sols::Vector{Dict{String, Vector{Number}}}, model::CP.FlatZinc.Model)

Parses the output of `_parse_to_assignments` and stores the solutions into
`model`. This function is responsible for filling `model.results`.
"""
function _parse_to_moi_solutions(sols::Vector{Dict{String, Vector{Number}}}, model::Optimizer)
    @assert length(model.results.primal_solutions) == 0

    for i in 1:length(sols)
        sol = sols[i]
        push!(model.results.primal_solutions, Dict{MOI.VariableIndex, Real}())

        for (var_name, val) in sol
            # At least for now: no vector at the FlatZinc layer. Just ensure
            # that this is consistent.
            @assert length(val) == 1

            var_idx = MOI.get(model, MOI.VariableIndex, var_name)
            model.results.primal_solutions[i][var_idx] = val[1]
        end
    end

    model.results.termination_status = (length(sols) == 0) ? MOI.INFEASIBLE : MOI.OPTIMAL
    model.results.primal_status = (length(sols) == 0) ? MOI.NO_SOLUTION : MOI.FEASIBLE_POINT

    @assert length(model.results.primal_solutions) == length(sols)
    return
end


text:
# # Cook panel under plane stress

# Source code: [`Cook-plane-stress_tut.jl`](Cook-plane-stress_tut.jl)

# ## Description

# In this example we investigate the well-known benchmark of a tapered panel
# under plane stress conditions known under the name of Cook.  The problem has
# been solved many times with a variety of finite element models  and hence the
# solution is well-known.

# ## Goals

# - Show how to generate the mesh by creating a rectangular block and reshaping it.
# - Execute the simulation with a static-equilibrium algorithm (solver).

##
# ## Definitions

# The problem is solved in a script.  We begin  by `using` the top-level module `FinEtools`.
# Further, we use the linear-deformation application package.
using FinEtools
using FinEtoolsDeforLinear

# With the algorithm modules, the problem can be set up (the materials, boundary
# conditions, and mesh are defined) and handed off to an algorithm (in this
# case linear static solution).  Then for postprocessing another set of
# algorithms can be invoked.
using FinEtoolsDeforLinear.AlgoDeforLinearModule

# A few  input parameters are defined: the material parameters. Note: the units
# are consistent, but unnamed.
E = 1.0;
nu = 1.0/3;

# The geometry of the tapered panel.
width = 48.0; height = 44.0; thickness  = 1.0;
free_height  = 16.0;
# Location of tracked  deflection is the midpoint of the loaded edge.
Mid_edge  = [48.0, 52.0];

# The tapered panel is loaded along the free edge with a unit force, which is
# here converted to loading per unit area.
magn = 1.0/free_height/thickness;# Magnitude of applied load

# For the above input parameters the converged displacement of the tip  of the
# tapered panel in the direction of the applied shear load is
convutip = 23.97;

# The mesh is generated as a rectangular block to begin with, and then the
# coordinates of the nodes are tweaked into the tapered panel shape. In this
# case we are using quadratic triangles (T6).
n = 10; # number of elements per side
fens, fes = T6block(width, height, n, n)

# Reshape the rectangle into a trapezoidal panel:
for i in 1:count(fens)
    fens.xyz[i,2] += (fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
end

# The  boundary conditions  are applied to selected finite element nodes.   The
# selection is based on the inclusion in a selection "box".
tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance
# Clamped edge of the membrane
l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance);

# The list of the selected nodes is then used twice,  to fix the degree of
# freedom  in the direction 1 and  in the direction 2. The essential-boundary
# condition data is stored in  dictionaries: `ess1` and `ess2 `.  These
# dictionaries  are used below to compose the computational model.
ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1);
ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1);

# The traction boundary condition is applied to the finite elements on the boundary of the panel. First we generate the three-node "curve" elements on the entire boundary of the panel.

boundaryfes =  meshboundary(fes);

# Then from these finite elements we choose the ones that are inside the box
# that captures the edge of the geometry to which the traction should be
# applied.
Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);

# To apply the traction we create a finite element model machine (FEMM). For the
# evaluation of the traction it is sufficient to create a  "base" FEMM.  It
# consists of the geometry data `IntegDomain` (connectivity,  integration rule,
# evaluation  of the basis functions  and basis function gradients with respect
# to the parametric coordinates). This object is composed of the list of the
# finite elements and  an appropriate quadrature rule (Gauss rule here).
el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 3), thickness));

# The traction boundary condition is specified with a constant traction vector and the FEMM that will be used to evaluate  the load vector.
flux1 = FDataDict("traction_vector"=>[0.0,+magn],
    "femm"=>el1femm
    );

# We make the dictionary for the region (the interior of the domain).  The FEMM
# for the evaluation of the integrals over the interior of the domain (that is
# the stiffness matrix) and the material are needed. The geometry data  now is
# equipped with the  triangular  three-point rule. Note the model-reduction
# type which is used to dispatch to appropriate specializations of the material
# routines and the FEMM which needs to execute different code for different
# reduced-dimension models. Here the model reduction is "plane stress".
MR = DeforModelRed2DStress
material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)
region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, TriRule(3), thickness), material));

# The model data is a dictionary.   In the present example it consists of the
# node set, the array of dictionaries for the regions, and arrays of
# dictionaries for each essential and natural boundary condition.
modeldata = FDataDict("fens"=>fens,
 "regions"=>[region1],
 "essential_bcs"=>[ess1, ess2],
 "traction_bcs"=>[flux1]
 );

# When the model data is defined, we simply pass it to the algorithm.
modeldata = AlgoDeforLinearModule.linearstatics(modeldata);

# The model data is augmented in the algorithm by the nodal field representing
# the geometry and the displacement field  computed by solving the system of
# linear algebraic equations of equilibrium.
u = modeldata["u"];
geom = modeldata["geom"];

# The complete information returned from the algorithm  is
@show keys(modeldata)

# Now we can extract the displacement at the mid-edge node and compare to the
# converged (reference) value. The code below selects the node inside a very
# small box of the size `tolerance` which presumably contains only a single
# node, the one at the midpoint of the edge.
nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
          inflate=tolerance);
theutip = u.values[nl,:]
println("displacement =$(theutip[2]) as compared to converged $convutip")

# For postprocessing  we will export a VTK file  with the displacement field
# (vectors)  and  one scalar field ($\sigma_{xy}$).
modeldata["postprocessing"] = FDataDict("file"=>"cookstress",
   "quantity"=>:Cauchy, "component"=>:xy);
modeldata = AlgoDeforLinearModule.exportstress(modeldata);

# The  attribute `"postprocessing"` holds additional data computed and returned
# by the algorithm:
@show keys(modeldata["postprocessing"])

# The exported data can be digested as follows: `modeldata["postprocessing"]
# ["exported"]` is an array of exported items.
display(keys(modeldata["postprocessing"]["exported"]))

# Each entry of the array is a dictionary:
display(keys(modeldata["postprocessing"]["exported"][1]))

# Provided we have  `paraview` in the PATH, we can bring it up  to display the
# exported data.
File = modeldata["postprocessing"]["exported"][1]["file"]
@async run(`"paraview.exe" $File`);

# We can also extract the minimum and maximum value of the shear stress
# (-0.06, and 0.12).
display(modeldata["postprocessing"]["exported"][1]["quantity"])
display(modeldata["postprocessing"]["exported"][1]["component"])
fld = modeldata["postprocessing"]["exported"][1]["field"]
println("$(minimum(fld.values)) $(maximum(fld.values))")

true


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
include("scenarios.jl")
println("Finished loading libraries."); flush(stdout)

flatten = Iterators.flatten

if length(ARGS) > 1
  #=
  SCRATCH_DIR = ENV["SCRATCH_DIR"]
  SIMNAME = "sim5"
  RESULTS_DIR = "$(SCRATCH_DIR)/cde/simstudy/$(SIMNAME)/results"
  AWS_BUCKET = "s3://cytof-density-estimation/simstudy/$(SIMNAME)"
  resultsdir = RESULTS_DIR
  awsbucket = AWS_BUCKET
  =#
  resultsdir = ARGS[1]
  awsbucket = ARGS[2]
  println("ARGS: ", ARGS)
  istest = false
  Ni = 1000
  nburn = 10000
  nsamps = 5000
  thin=2
  Ks = [2,4,6]
  p=0.5
else
  resultsdir = "results/test/"
  awsbucket = nothing
  istest = true
  Ni = 1000
  nburn = 300
  nsamps = 300
  thin=1
  Ks = [2]
  p=0.5
end
flush(stdout)

using Distributed
println("Load libraries on workers ..."); flush(stdout)
if !istest
  rmprocs(workers())
  addprocs(12)
  @everywhere include("imports.jl")
end
println("Finished loading libraries on workers."); flush(stdout)

# Helpers.
simname(K, snum, beta) = "K$(K)/scenario$(snum)/m$(beta)"
function make_bucket(K, snum, beta, bucket)
  return (bucket == nothing) ? bucket : "$(bucket)/$(simname(K, snum, beta))"
end
make_resdir(K, snum, beta) = "$(resultsdir)/$(simname(K, snum, beta))"
function make_imgdir(K, snum, beta)
  resdir = make_resdir(K, snum, beta)
  imgdir = "$(resdir)/img"
  mkpath(imgdir)
  return imgdir
end

configs = [[let
  # NOTE
  simdata = scenarios(snum, seed=1, Ni=Ni)
  bucket = make_bucket(K, snum, beta, awsbucket)
  resdir = make_resdir(K, snum, beta)
  imgdir = make_imgdir(K, snum, beta)
  (awsbucket=bucket, simdata=simdata, resultsdir=resdir, imgdir=imgdir,
   snum=snum, beta=beta, K=K, nsamps=nsamps, thin=thin, nburn=nburn, p=p)
end for beta in 0:1] for K in Ks for snum in 1:4]

# Parallel run.
println("Starting runs ..."); flush(stdout)
res = pmap(run, istest ? configs[1] : flatten(configs), on_error=identity)

println("Status of runs:"); flush(stdout)
foreach(z -> println(z[2], " => ", z[1]),
        zip(res, getfield.(istest ? configs[1] : flatten(configs), :resultsdir)))

# Bayes factor
println("Compute BF:"); flush(stdout)
res = pmap(c -> let
  # NOTE
  @assert length(unique(getfield.(c, :snum))) == 1
  resd0 = c[1][:resultsdir]
  resd1 = c[2][:resultsdir]
  out0 = BSON.load("$(resd0)/results.bson")
  out1 = BSON.load("$(resd1)/results.bson")
  imgdir = "$(c[1][:imgdir])/../../img"; mkpath(imgdir)
  bucket = "$(c[1][:awsbucket])/../img"
  postprocess(out0[:chain], out1[:chain], out0[:data],
              imgdir, bucket, simdata=c[1][:simdata],
              density_legend_pos=:topleft, bw_postpred=.3, binsC=50, binsT=100,
              p=c[1][:p])
end, istest ? [configs[1]] : configs, on_error=identity)

println("Status of BF computation:"); flush(stdout)
foreach(z -> println(z[2], " => ", z[1]),
        zip(res, [c[1][:snum] for c in (istest ? [configs[1]] : configs)]))

println("DONE WITH ALL!"); flush(stdout)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
"""
测试ult时的计算
"""

using Plots
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Kagome
using MARY_fRG.FlowEquation
using Profile

"""
获取极低温的时候的bubble
"""
function get_ult_bubb(Γ4; maxnum=10000)
    Γ4_ult = fliter_away_surface(Γ4, 2.0)
    Γ4_ult = refine_to_surface(Γ4_ult)
    while length(Γ4_ult.ltris) < maxnum
        Γ4_ult = refine_to_surface(Γ4_ult)
    end
    #plt = draw_points([tri.center for tri in Γ4_ult.ltris])
    #@savepng plt "ult_surface"
    maxi, mini = engpeak_to_surface(Γ4_ult)
    bubb_pp_u, bubb_fs_u, bubb_ex_u = @profile all_bubble_tf_ult_mt(Γ4_ult, maxi)
    return bubb_pp_u, bubb_fs_u, bubb_ex_u
end

"""
运行
"""
function run()
    model = upperband_kagome_lattice(0.)
    Γ4 = TFGamma4(
        model, 6.0, 24, 200
    )
    bpu, bfu, beu = get_ult_bubb(Γ4)
    #
    println(all(isfinite.(bpu.V)))
    println(all(isfinite.(bfu.V)))
    println(all(isfinite.(beu.V)))
    #
    umat = bpu.V[1, 1, 1, 1, :, :, :]
    maxv, maxi = findmin(isfinite.(umat))
    println(maxi, maxv)
    println(umat[maxi])
end

run()

run()

outf = open("kag_profile", "w")
Profile.print(outf, format = :flat)
close(outf)


text:
using Test
using QuantumControlBase

@testset "chain_infohooks" begin

    function hook1(wrk, a1, a2, args...)
        return (1, )
    end

    function hook2(wrk, a1, a2, args...)
        return (2, )
    end

    function hook3(wrk, a1, a2, args...)
        @test length(args) == 2
        return (3.1, 3.2)
    end

    function hook4(wrk, a1, a2, args...)
        return nothing
    end

    chained = chain_infohooks(hook1, hook2, hook3, hook4)

    res = chained(nothing, 0, 0)

    @test res == (1, 2, 3.1, 3.2)
    @test res[1] isa Int64
    @test res[3] isa Float64

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the general functions using quaternions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/quaternion.jl
# =========================

# Functions: conj
# ---------------

@testset "General functions of quaternions: conj" begin
    q = Quaternion{Int}(1, 2, 3, 4)
    qc = conj(q)
    @test qc.q0 == +1
    @test qc.q1 == -2
    @test qc.q2 == -3
    @test qc.q3 == -4
    @test eltype(qc) === Int

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    qc = conj(q)
    @test qc.q0 ≈ +q.q0
    @test qc.q1 ≈ -q.q1
    @test qc.q2 ≈ -q.q2
    @test qc.q3 ≈ -q.q3
    @test eltype(qc) === eltype(q)

    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    qc = conj(q)
    @test qc.q0 ≈ +q.q0
    @test qc.q1 ≈ -q.q1
    @test qc.q2 ≈ -q.q2
    @test qc.q3 ≈ -q.q3
    @test eltype(qc) === eltype(q)
end

# Functions: copy
# ---------------

@testset "General functions of quaternions: copy" begin
    q = Quaternion{Int}(1, 2, 3, 4)
    qc = copy(q)
    @test q === qc
    @test eltype(qc) === Int

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    qc = copy(q)
    @test q === qc
    @test eltype(qc) === eltype(q)

    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    qc = copy(q)
    @test q === qc
    @test eltype(qc) === eltype(q)
end

# Functions: imag
# ---------------

@testset "General functions of quaternions: imag" begin
    q = Quaternion{Int}(1, 2, 3, 4)
    v = imag(q)
    @test v[1] == 2
    @test v[2] == 3
    @test v[3] == 4
    @test eltype(v) === Int

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    v = imag(q)
    @test v[1] == q.q1
    @test v[2] == q.q2
    @test v[3] == q.q3
    @test eltype(v) === eltype(q)

    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    v = imag(q)
    @test v[1] == q.q1
    @test v[2] == q.q2
    @test v[3] == q.q3
    @test eltype(v) === eltype(q)
end

# Function: inv
# -------------

@testset "General functions of quaternions: inv" begin
    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    qi = inv(q)
    norm_q² = q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3
    @test qi.q0 == +q.q0 / norm_q²
    @test qi.q1 == -q.q1 / norm_q²
    @test qi.q2 == -q.q2 / norm_q²
    @test qi.q3 == -q.q3 / norm_q²
    @test eltype(qi) === Float64

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    qi = inv(q)
    norm_q² = q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3
    @test qi.q0 == +q.q0 / norm_q²
    @test qi.q1 == -q.q1 / norm_q²
    @test qi.q2 == -q.q2 / norm_q²
    @test qi.q3 == -q.q3 / norm_q²
    @test eltype(qi) === Float32
end

# Function: norm
# --------------

@testset "General functions of quaternions: norm" begin
    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    norm_q = norm(q)
    @test norm(q) ≈ sqrt(q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3)
    @test norm_q isa Float64

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    norm_q = norm(q)
    @test norm(q) ≈ sqrt(q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3)
    @test norm_q isa Float32
end

# Function: real
# --------------

@testset "General functions of quaternions: real" begin
    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    real_q = real(q)
    @test real_q == q.q0
    @test real_q isa Float64

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    real_q = real(q)
    @test real_q == q.q0
    @test real_q isa Float32

    q = Quaternion{Int}(1, 0, 0, 0)
    real_q = real(q)
    @test real_q == 1
    @test real_q isa Int
end

# Functions: vect
# ---------------

@testset "General functions of quaternions: vect" begin
    q = Quaternion{Int}(1, 2, 3, 4)
    v = vect(q)
    @test v[1] == 2
    @test v[2] == 3
    @test v[3] == 4
    @test eltype(v) === Int

    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    v = vect(q)
    @test v[1] == q.q1
    @test v[2] == q.q2
    @test v[3] == q.q3
    @test eltype(v) === eltype(q)

    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    v = vect(q)
    @test v[1] == q.q1
    @test v[2] == q.q2
    @test v[3] == q.q3
    @test eltype(v) === eltype(q)
end

# Functions: one
# --------------

@testset "General functions of quaternions: one" begin
    q = one(Quaternion)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Float64

    q = one(Quaternion{Float32})
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Float32

    qi = Quaternion(1, 0, 0, 0)
    q = one(qi)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Int
end

# Functions: zero
# ---------------

@testset "General functions of quaternions: zero" begin
    q = zero(Quaternion)
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Float64

    q = zero(Quaternion{Float32})
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Float32

    qi = Quaternion(1, 0, 0, 0)
    q = zero(qi)
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Int
end

# Functions: getindex, length, setindex!
# --------------------------------------

@testset "General functions of quaternions: iterable object API" begin
    q = Quaternion{Float64}(randn(), randn(), randn(), randn())
    @test length(q) == 4
    @test q[1] == q.q0
    @test q[2] == q.q1
    @test q[3] == q.q2
    @test q[4] == q.q3

    @test_throws BoundsError q[0]
    @test_throws BoundsError q[5]

    v = zeros(10)

    v[4:7] = q

    @test v[1]  == 0
    @test v[2]  == 0
    @test v[3]  == 0
    @test v[4]  == q.q0
    @test v[5]  == q.q1
    @test v[6]  == q.q2
    @test v[7]  == q.q3
    @test v[8]  == 0
    @test v[9]  == 0
    @test v[10] == 0

    @test firstindex(q) === 1
    @test lastindex(q) === 4
end

# Functions: show
# ---------------

@testset "General functions of quaternions: show" begin
    buf = IOBuffer()
    io = IOContext(buf)
    q = Quaternion{Float64}(1, 2, exp(1), π)

    # Extended printing.
    show(io, MIME"text/plain"(), q)
    expected = """
        Quaternion{Float64}:
          + 1.0 + 2.0⋅i + 2.71828⋅j + 3.14159⋅k"""
    @test String(take!(io.io)) == expected

    # Comapct printing.
    show(io, q)
    expected = """
        Quaternion{Float64}: + 1.0 + 2.0⋅i + 2.71828⋅j + 3.14159⋅k"""
    @test String(take!(io.io)) == expected

    # Colors.
    io = IOContext(buf, :color => true)
    show(io, MIME"text/plain"(), q)
    expected = """
        Quaternion{Float64}:
          + 1.0 + 2.0⋅\e[1mi\e[0m + 2.71828⋅\e[1mj\e[0m + 3.14159⋅\e[1mk\e[0m"""
    @test String(take!(io.io)) == expected
end


text:
using KernelAbstractions
using KernelGradients
using Enzyme
using CUDA
using CUDAKernels
using Test

include(joinpath(dirname(pathof(KernelAbstractions)), "..", "test", "testsuite.jl"))
include(joinpath(dirname(pathof(KernelGradients)), "..", "test", "testsuite.jl"))

CI = parse(Bool, get(ENV, "CI", "false"))
if CI
    default = "CPU"
else
    default = "CUDA"
end

backend = get(ENV, "KERNELABSTRACTIONS_TEST_BACKEND", default)

if backend != "CUDA"
    @info "CUDA backend not selected"
    exit()
end

# if CUDA.functional(true)
if CUDA.functional()
    CUDA.versioninfo()
    CUDA.allowscalar(false)
    Testsuite.testsuite(CUDADevice, backend, CUDA, CuArray, CUDA.CuDeviceArray)
    # GradientsTestsuite.testsuite(CUDADevice, backend, CUDA, CuArray, CUDA.CuDeviceArray)
elseif !CI
    error("No CUDA GPUs available!")
end


text:
"""
    AbstractScheduler

Documentation placeholder for AbstractScheduler

See also: [`scheduled_subscription!`](@ref), [`scheduled_next!`](@ref), [`scheduled_error!`](@ref), [`scheduled_complete!`](@ref)
"""
abstract type AbstractScheduler end

"""
    getscheduler

Documentation placeholder for getscheduler

See also: [`scheduled_subscription!`](@ref), [`scheduled_next!`](@ref), [`scheduled_error!`](@ref), [`scheduled_complete!`](@ref)
"""
function getscheduler end

"""
    scheduled_subscription!

Documentation placeholder for scheduled_subscription!

See also: [`getscheduler`](@ref), [`scheduled_next!`](@ref), [`scheduled_error!`](@ref), [`scheduled_complete!`](@ref)
"""
function scheduled_subscription! end

"""
    scheduled_next!

Documentation placeholder for scheduled_next!

See also: [`getscheduler`](@ref), [`scheduled_subscription!`](@ref), [`scheduled_error!`](@ref), [`scheduled_complete!`](@ref)
"""
function scheduled_next! end

"""
    scheduled_error!

Documentation placeholder for scheduled_error!

See also: [`getscheduler`](@ref), [`scheduled_subscription!`](@ref), [`scheduled_next!`](@ref), [`scheduled_complete!`](@ref)
"""
function scheduled_error! end

"""
    scheduled_complete!

Documentation placeholder for scheduled_complete!

See also: [`getscheduler`](@ref), [`scheduled_subscription!`](@ref), [`scheduled_next!`](@ref), [`scheduled_error!`](@ref)
"""
function scheduled_complete! end

"""
    makeinstance

Documentation placeholder for makeinstance

See also: [`getscheduler`](@ref), [`scheduled_subscription!`](@ref), [`scheduled_next!`](@ref), [`scheduled_error!`](@ref)
"""
function makeinstance end

"""
    instancetype

Documentation placeholder for makeinstance

See also: [`getscheduler`](@ref), [`scheduled_subscription!`](@ref), [`scheduled_next!`](@ref), [`scheduled_error!`](@ref)
"""
function instancetype end


text:
#=
    JumpProblem(rn::ReactionSystemMod, prob::DiscreteProblem, aggregator, args...; kwargs...)

Given a `ReactionSystemMod`, DifferentialEquations [`DiscreteProblem`]
(https://diffeq.sciml.ai/stable/types/discrete_types/#Discrete-Problems) and [the method
for aggregating constant jumps](https://diffeq.sciml.ai/stable/types/jump_types/#Jump-Problems),
return a corresponding  [`JumpProblem`](https://diffeq.sciml.ai/stable/types/jump_types/#Jump-Problems)
which can be solved using standard DifferentialEquations methods.

# Notes
- This function essentially massages a `ReactionSystemMod` into a form which could be used
  with DifferentialEquations SSA solvers. Namely, a `ReactionSystemMod` is converted into a
  system of [ConstantRateJumps]
  (https://diffeq.sciml.ai/stable/types/jump_types/#Types-of-Jumps:-Regular,-Variable,-Constant-Rate-and-Mass-Action)
  taking into account the geometrically distributed reaction products.
- The implementation here is very restricted: time-dependent propensities are currently not
  allowed; [MassActionJumps]
  (https://diffeq.sciml.ai/stable/types/jump_types/#Types-of-Jumps:-Regular,-Variable,-Constant-Rate-and-Mass-Action)
  cannot be used due to limitations of the `ReactionSystemMod` API, hence the SSA performance is subpar.
=#
function DiffEqJump.JumpProblem(rn::ReactionSystemMod, prob::DiscreteProblem, aggregator, args...; kwargs...)

    # TODO:
    # - add random seed for sampling from geometric
    # - allow VariableRateJumps
    # - using MassActionJumps where possible would be significantly faster but they require
    #  the reactant stoichometry matrix and knowing which reactions are mass-action (not possible with current API)

    # construct the mapping from each paramater to their value if the DiscreteProblem
    # contains a vector of parameter values
    if prob.p isa Array{<:Real, 1}
        pmap = Dict(Pair.(params(rn), prob.p / 1.))
    else
        # otherwise we assume `prob` already contains the required mapping
        pmap = Dict(prob.p)
    end

    # build rate functions
    as = [substitute(a, pmap) for a in propensities(rn)]
    # in case of VariableRateJumps or symbolic params this should be amended
    af = [eval(build_function(a, species(rn)...)) for a in as]

    S_mat = get_S_mat(rn)

    is_geometric = fill(false, size(S_mat))
    rngs = similar(S_mat, Any)

    for (i, j) in Tuple.(findall(x->isa(x, Sym), S_mat))
        is_geometric[i, j] = true
        m = pmap[S_mat[i, j]]
        # note that Geometric considers success probability p
        # the mean is (1-p)/p, hence p = 1/(mean+1)
        rngs[i, j] = Geometric(1/(m+1))
    end

    jumps = []

    for r in 1:numreactions(rn)

        rate(u, p, t) = af[r](u...)

        if sum(is_geometric[:, r]) > 0

            inds_n = findall(s->!s, is_geometric[:, r])
            inds_g = findall(is_geometric[:, r])
            function affect1!(integrator)
                integrator.u[inds_n] .+= S_mat[inds_n, r]
                for ind in inds_g
                    integrator.u[ind] += rand(rngs[ind, r])
                end
            end
            push!(jumps, ConstantRateJump(rate, affect1!))
        else

            function affect2!(integrator)
                integrator.u .+= S_mat[:, r]
            end
            push!(jumps, ConstantRateJump(rate, affect2!))
        end

    end

    return JumpProblem(prob, aggregator, jumps..., args...; kwargs...)

end


text:
function Worksheet(xf::XLSXFile, sheet_element::EzXML.Node)
    @assert EzXML.nodename(sheet_element) == "sheet"
    sheetId = parse(Int, sheet_element["sheetId"])
    relationship_id = sheet_element["r:id"]
    name = sheet_element["name"]
    dim = read_worksheet_dimension(xf, relationship_id, name)

    return Worksheet(xf, sheetId, relationship_id, name, dim)
end

# 18.3.1.35 - dimension (Worksheet Dimensions). This is optional, and not required.
function read_worksheet_dimension(xf::XLSXFile, relationship_id, name) :: Union{Nothing, CellRange}
    local result::Union{Nothing, CellRange} = nothing

    wb = get_workbook(xf)
    target_file = get_relationship_target_by_id("xl", wb, relationship_id)
    zip_io, reader = open_internal_file_stream(xf, target_file)

    try
        # read Worksheet dimension
        while EzXML.iterate(reader) != nothing
            if EzXML.nodetype(reader) == EzXML.READER_ELEMENT && EzXML.nodename(reader) == "dimension"
                @assert EzXML.nodedepth(reader) == 1 "Malformed Worksheet \"$(ws.name)\": unexpected node depth for dimension node: $(EzXML.nodedepth(reader))."
                ref_str = reader["ref"]
                if is_valid_cellname(ref_str)
                    result = CellRange("$(ref_str):$(ref_str)")
                else
                    result = CellRange(ref_str)
                end

                break
            end
        end
    finally
        close(reader)
        close(zip_io)
    end

    return result
end

@inline isdate1904(ws::Worksheet) = isdate1904(get_workbook(ws))

# Retuns the dimension of this worksheet as a CellRange.
# Returns `nothing` if the dimension is unknown.
@inline get_dimension(ws::Worksheet) :: Union{Nothing, CellRange} = ws.dimension

function set_dimension!(ws::Worksheet, rng::CellRange)
    ws.dimension = rng
    nothing
end

"""
    getdata(sheet, ref)
    getdata(sheet, row, column)

Returns a escalar or a matrix with values from a spreadsheet.
`ref` can be a cell reference or a range.

Indexing in a `Worksheet` will dispatch to `getdata` method.

# Example

```julia
julia> f = XLSX.readxlsx("myfile.xlsx")

julia> sheet = f["mysheet"]

julia> matrix = sheet["A1:B4"]

julia> single_value = sheet[2, 2] # B2
```

See also [`XLSX.readdata`](@ref).
"""
getdata(ws::Worksheet, single::CellRef) = getdata(ws, getcell(ws, single))
getdata(ws::Worksheet, row::Integer, col::Integer) = getdata(ws, CellRef(row, col))

function getdata(ws::Worksheet, rng::CellRange) :: Array{Any,2}
    result = Array{Any, 2}(undef, size(rng))
    fill!(result, missing)

    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    left = column_number(rng.start)
    right = column_number(rng.stop)

    for sheetrow in eachrow(ws)
        if top <= sheetrow.row && sheetrow.row <= bottom
            for column in left:right
                cell = getcell(sheetrow, column)
                if !isempty(cell)
                    (r, c) = relative_cell_position(cell, rng)
                    result[r, c] = getdata(ws, cell)
                end
            end
        end

        # don't need to read new rows
        if sheetrow.row > bottom
            break
        end
    end

    return result
end

function getdata(ws::Worksheet, rng::ColumnRange) :: Array{Any,2}
    columns_count = length(rng)
    columns = Vector{Vector{Any}}(undef, columns_count)
    for i in 1:columns_count
        columns[i] = Vector{Any}()
    end

    left, right = column_bounds(rng)

    for sheetrow in eachrow(ws)
        for column in left:right
            cell = getcell(sheetrow, column)
            c = relative_column_position(cell, rng) # r will be ignored
            push!(columns[c], getdata(ws, cell))
        end
    end

    rows = length(columns[1])
    for i in 1:columns_count
        @assert length(columns[i]) == rows "Inconsistent state: Each column should have the same number of rows."
    end

    return hcat(columns...)
end

function getdata(ws::Worksheet, ref::AbstractString) :: Union{Array{Any,2}, Any}
    if is_valid_cellname(ref)
        return getdata(ws, CellRef(ref))
    elseif is_valid_cellrange(ref)
        return getdata(ws, CellRange(ref))
    elseif is_valid_column_range(ref)
        return getdata(ws, ColumnRange(ref))
    elseif is_worksheet_defined_name(ws, ref)
        v = get_defined_name_value(ws, ref)
        if is_defined_name_value_a_constant(v)
            return v
        elseif is_defined_name_value_a_reference(v)
            return getdata(ws, v)
        else
            error("Unexpected defined name value: $v.")
        end
    elseif is_workbook_defined_name(get_workbook(ws), ref)
        wb = get_workbook(ws)
        v = get_defined_name_value(wb, ref)
        if is_defined_name_value_a_constant(v)
            return v
        elseif is_defined_name_value_a_reference(v)
            return getdata(get_xlsxfile(ws), v)
        else
            error("Unexpected defined name value: $v.")
        end
    else
        error("$ref is not a valid cell or range reference.")
    end
end

getdata(ws::Worksheet, rng::SheetCellRange) = getdata(get_xlsxfile(ws), rng)

function getdata(ws::Worksheet)
    if ws.dimension != nothing
        return getdata(ws, get_dimension(ws))
    else
        error("Worksheet dimension is unknown.")
    end
end

Base.getindex(ws::Worksheet, r) = getdata(ws, r)
Base.getindex(ws::Worksheet, r, c) = getdata(ws, r, c)
Base.getindex(ws::Worksheet, ::Colon) = getdata(ws)

function Base.show(io::IO, ws::Worksheet)
    if get_dimension(ws) != nothing
        rg = get_dimension(ws)
        nrow, ncol = size(rg)
        @printf(io, "%d×%d %s: [\"%s\"](%s)", nrow, ncol, typeof(ws), ws.name, rg)
    else
        @printf(io, "%s: [\"%s\"]", typeof(ws), ws.name)
    end
end

"""
    getcell(sheet, ref)

Returns an `AbstractCell` that represents a cell in the spreadsheet.

Example:

```julia
julia> xf = XLSX.readxlsx("myfile.xlsx")

julia> sheet = xf["mysheet"]

julia> cell = XLSX.getcell(sheet, "A1")
```
"""
function getcell(ws::Worksheet, single::CellRef) :: AbstractCell

    for sheetrow in eachrow(ws)
        if row_number(sheetrow) == row_number(single)
            return getcell(sheetrow, column_number(single))
        end
    end

    return EmptyCell(single)
end

function getcell(ws::Worksheet, ref::AbstractString)
    if is_valid_cellname(ref)
        return getcell(ws, CellRef(ref))
    else
        error("$ref is not a valid cell reference.")
    end
end

getcell(ws::Worksheet, row::Integer, col::Integer) = getcell(ws, CellRef(row, col))

"""
    getcellrange(sheet, rng)

Returns a matrix with cells as `Array{AbstractCell, 2}`.
`rng` must be a valid cell range, as in `"A1:B2"`.
"""
function getcellrange(ws::Worksheet, rng::CellRange) :: Array{AbstractCell,2}
    result = Array{AbstractCell, 2}(undef, size(rng))
    for cellref in rng
        (r, c) = relative_cell_position(cellref, rng)
        result[r, c] = EmptyCell(cellref)
    end

    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    left = column_number(rng.start)
    right = column_number(rng.stop)

    for sheetrow in eachrow(ws)
        if top <= sheetrow.row && sheetrow.row <= bottom
            for column in left:right
                cell = getcell(sheetrow, column)
                if !isempty(cell)
                    (r, c) = relative_cell_position(cell, rng)
                    result[r, c] = cell
                end
            end
        end

        # don't need to read new rows
        if sheetrow.row > bottom
            break
        end
    end

    return result
end

function getcellrange(ws::Worksheet, rng::ColumnRange) :: Array{AbstractCell,2}
    columns_count = length(rng)
    columns = Vector{Vector{AbstractCell}}(undef, columns_count)
    for i in 1:columns_count
        columns[i] = Vector{AbstractCell}()
    end

    let
        left, right = column_bounds(rng)

        for sheetrow in eachrow(ws)
            for column in left:right
                cell = getcell(sheetrow, column)
                c = relative_column_position(cell, rng) # r will be ignored
                push!(columns[c], cell)
            end
        end
    end

    rows = length(columns[1])
    for i in 1:columns_count
        @assert length(columns[i]) == rows "Inconsistent state: Each column should have the same number of rows."
    end

    return hcat(columns...)
end

function getcellrange(ws::Worksheet, rng::AbstractString)
    if is_valid_cellrange(rng)
        return getcellrange(ws, CellRange(rng))
    elseif is_valid_column_range(rng)
        return getcellrange(ws, ColumnRange(rng))
    else
        error("$rng is not a valid cell range.")
    end
end


text:
module TaylorBigF

    include("calc_Vandermonde.jl")
    include("calc_Vandermonde_inv_formula.jl")
    include("calc_Vandermonde_dot.jl")
    include("calc_Vandermonde_dot_dot.jl")
    include("calc_Vandermonde_dot_dot_dot.jl")
    include("calc_Vandermonde_dot_dot_dot_dot.jl")
    include("calc_Vandermonde_int.jl")
    include("calc_Vandermonde_int_int.jl")
    include("calc_Vandermonde_2D.jl")
    include("calc_Vandermonde_xxxx_2D.jl")
    include("calc_Vandermonde_yyyy_2D.jl")
    include("calc_Vandermonde_xxyy_2D.jl")
    include("calc_Vandermonde_xxx_2D.jl")
    include("calc_Vandermonde_yyy_2D.jl")
    include("calc_Vandermonde_xyy_2D.jl")
    include("calc_Vandermonde_yxx_2D.jl")

    export calc_Vandermonde,calc_Vandermonde_inv_formula,calc_Vandermonde_dot,
    calc_Vandermonde_dot_dot,calc_Vandermonde_dot_dot_dot,calc_Vandermonde_dot_dot_dot_dot,
    calc_Vandermonde_int,calc_Vandermonde_int_int,calc_Vandermonde_2D,
    calc_Vandermonde_xxxx_2D,calc_Vandermonde_yyyy_2D,calc_Vandermonde_xxyy_2D,
    calc_Vandermonde_xxx_2D,calc_Vandermonde_yyy_2D,calc_Vandermonde_xyy_2D,calc_Vandermonde_yxx_2D

end


text:
# Using Base modules.
using Random
# Load the distributions library.
using Distributions

mutable struct UGMM_CAVI <: InferenceAlgorithm
    X # Data
    K::Int       # number of iterations
    N::Int
    sigma2
    m
    phi
    s2
end

function UGMM_CAVI(X, num_components::Int, n_samples::Int, sigma2=1)
    K=num_components
    N=n_samples
    rng = RandomDevice()
    phi = rand(Dirichlet(K*N, (rand(rng, Float32,(1))*rand(rng, 1:9))[1]),(K));
    m = rand(rng,floor(minimum(X)):floor(maximum(X)), (1,K))
    m += maximum(X)*rand(rng, Float32, (1,K))
    s2 = rand(rng, Float32, (1, K))
    print("\nInit mean\n")
    println(m)
    print("\nInit s2\n")
    println(s2)
    UGMM_CAVI(X, num_components, n_samples, sigma2, m, phi, s2)

end

function get_elbo(alg::UGMM_CAVI)
    t1 = log.(alg.s2) - alg.m./(alg.sigma2)
    t1 = sum(t1)
    t2 = -0.5.*(reshape(alg.X, (alg.N*alg.K, 1)).^2 * ones(1, alg.K) + ones(alg.N*alg.K, 1)*((alg.s2)+(alg.m).^2 ))
    t2 .+= reshape(alg.X, (alg.N*alg.K, 1))*(alg.m)
    t2 .-= log.(alg.phi)
    t2 .+= alg.phi
    t2 = sum(t2)
    return t1 + t2
end

function _update_phi(alg::UGMM_CAVI)
    t1 = reshape(alg.X, (alg.N*alg.K, 1))*alg.m
    t2 = -(0.5*alg.m.^2 + 0.5*alg.s2)
    exponent = t1 .+ t2
    phi = exp.(exponent)
    alg.phi = phi ./ sum(phi,dims=2)
end
function _update_mu(alg::UGMM_CAVI)
    alg.m = sum((alg.phi.*reshape(alg.X, (alg.N*alg.K, 1))),dims=1) .* (1 ./alg.sigma2 .+ sum(alg.phi,dims=1)).^(-1)
    @assert size(alg.m)[2] == alg.K
    #print(self.m)
    alg.s2 = (1 ./alg.sigma2 .+ sum(alg.phi,dims=1)).^(-1)
    @assert size(alg.s2)[2] == alg.K
end

function _cavi_step(alg::UGMM_CAVI)
    _update_phi(alg)
    _update_mu(alg)
end

function fit(alg::UGMM_CAVI; max_iter=100, tol=1e-20, print_interval=50)
    elbo_values = [get_elbo(alg)]
    m_history = [alg.m]
    s2_history = [alg.s2]
    for iter_ in 1:(max_iter+1)
        _cavi_step(alg)
        append!(m_history, [alg.m])
        append!(s2_history, [alg.s2])
        append!(elbo_values, [get_elbo(alg)])

        if iter_ % print_interval == 0
            println(iter_, " Means: ", m_history[iter_], " ELBO: ", elbo_values[end])
        end
        if abs(elbo_values[end-1] - elbo_values[end]) < tol
            i=iter_
            elbo_val = elbo_values[end]
            println("ELBO converged with $elbo_val at iteration $i")
            break
        end
        if iter_ == max_iter
            elbo_val = elbo_values[end]
            print("ELBO ended with ll $elbo_val")
        end
    end
    elbo_values, m_history ,s2_history
end


text:
## -------------- Rank is a weighted average of endowment percentiles

"""
	$(SIGNATURES)

Switches for linear endowment percentile weights.

`eNameV` contains the endowments used for ranking. It must be possible to call `retrieve_draws(EndowmentDraws, eName)`.

`wtV` are the weights to be used if not calibrated. This omits the first weight, which is fixed at 1 (or -1).
This is empty when there is only one endowment to rank on.
Weights are bounded in the interval `lbV` to `ubV`. Weights may be negative.

`doCal` determines whether weights are calibrated or fixed.
"""
mutable struct EndowPctRankingSwitches{F1} <: AbstractRankingSwitches{F1}
    # Endowments to rank on
    eNameV :: Vector{Symbol}
    # Weights on those endowments, if fixed. First is omitted as fixed.
    wtV :: Vector{F1}
    # Bounds on the weights (first again omitted)
    lbV :: Vector{F1}
    ubV :: Vector{F1}
    # Calibrate weights? The first is always fixed (normalization).
    doCal :: Bool
    # High first endowments => rank first?
    highDrawsFirst :: Bool
end

"""
	$(SIGNATURES)

Rank students by a linear combination of endowment percentiles.
"""
mutable struct EndowPctRanking{F1} <: AbstractRanking{F1}
    objId :: ObjectId
    pvec :: ParamVector
    # Weights, excluding the first, normalized one.
    wtV :: Vector{F1}
    switches :: EndowPctRankingSwitches{F1}
end

# ------  Access and show

Base.show(io :: IO,  e :: EndowPctRankingSwitches{F1}) where F1 =
    print(io, typeof(e), " with endowments ", endow_names(e));

Base.show(io :: IO,  e :: EndowPctRanking{F1}) where F1 =
    print(io, typeof(e), " with endowments ", endow_names(e),
        " and weights ",  round.(weights(e), digits = 2));

function StructLH.describe(e :: EndowPctRankingSwitches{F1}) where F1
    endowNameV = CollegeEntry.endow_names(e);
    calStr = calibrate_weights(e)  ?  "calibrated"  :  "fixed";
    return [
        "Ranking for sequential entry"  " ";
        "Ranking based on"  "$endowNameV";
        "Weights on endowments:"  "$calStr"
    ]
end

# StructLH.describe(e :: EndowPctRanking) = StructLH.describe(e.switches);

# high_draws_first(e :: EndowPctRanking{F1}) where F1 =
#     high_draws_first(e.switches);

high_draws_first(switches :: EndowPctRankingSwitches{F1}) where F1 =
    switches.highDrawsFirst;

first_weight(switches :: EndowPctRankingSwitches{F1}) where F1 =
    high_draws_first(switches) ? one(F1) : -one(F1);

fixed_weights(switches :: EndowPctRankingSwitches{F1}) where F1 =
    [first_weight(switches), switches.wtV...];

weights(e :: EndowPctRanking{F1}) where F1 =
    [first_weight(e.switches), e.wtV...];

calibrate_weights(switches :: EndowPctRankingSwitches{F1}) where F1 =
    switches.doCal;
# calibrate_weights(e :: EndowPctRanking{F1}) where F1 =
#     calibrate_weights(e.switches);

lower_bounds(e :: EndowPctRanking{F1}) where F1 = zeros(F1, length(endow_names(e)));
upper_bounds(e :: EndowPctRanking{F1}) where F1 = ones(F1, length(endow_names(e)));

"""
	$(SIGNATURES)

Set bounds for one endowment.

# Example
```
set_bounds!(switches, :ability, -1.0, 1.0);
```
"""
function set_bounds!(switches :: EndowPctRankingSwitches{F1}, eName :: Symbol,
    lb :: F1, ub :: F1) where F1

    idx = findfirst(endow_names(switches) .== eName);
    @assert idx > 1  "Cannot set bounds for $eName in $switches"
    switches.lbV[idx-1] = lb;
    switches.ubV[idx-1] = ub;
    switches.wtV[idx-1] = 0.5 * (lb + ub);
    return nothing
end

# ------  Constructing

"""
	$(SIGNATURES)

Constructor with keyword arguments.
Properly handles the case of a single endowment.
"""
function EndowPctRankingSwitches(eNameV :: Vector{Symbol};
    wtInV = nothing, lbInV = nothing, ubInV = nothing,
    doCalIn :: Bool = true, highDrawsFirst :: Bool = true)

    if length(eNameV) == 1
        wtV = Vector{Float64}();
        lbV = Vector{Float64}();
        ubV = Vector{Float64}();
        doCal = false;
    else
        n = length(eNameV) - 1;
        wtV = isnothing(wtInV) ? fill(1.0, n) : wtInV;
        lbV = isnothing(lbInV) ? zeros(n) : lbInV;
        ubV = isnothing(ubInV) ? fill(3.0, n) : ubInV;
        doCal = doCalIn;
    end
    return EndowPctRankingSwitches(eNameV, wtV, lbV, ubV, doCal, highDrawsFirst);
end

"""
	$(SIGNATURES)

Rank on a single endowment.
"""
EndowPctRankingSwitches(eName :: Symbol; highDrawsFirst :: Bool = true) =
    EndowPctRankingSwitches([eName]; highDrawsFirst = highDrawsFirst);

function make_test_endowpct_switches(n :: Integer, highDrawsFirst :: Bool)
    eNameV = [:abilPct, :parentalPct, :hsGpaPct, :h0Pct];
    lbV = [-1.0, 0.0, 0.2];
    ubV = [0.0, 3.0, 2.0];
    wtV = 0.6 .* lbV .+ 0.4 .* ubV;
    # This produces an [] wtV when n == 1
    e = EndowPctRankingSwitches(eNameV[1 : n],
        wtV[1 : (n-1)], lbV[1 : (n-1)], ubV[1 : (n-1)], true, highDrawsFirst);
    @assert validate_ranking_switches(e);
    return e
end

function make_student_ranking(objId :: ObjectId,
    switches :: EndowPctRankingSwitches{F1},
    st :: SymbolTable) where F1
    # switches = make_test_endowpct_switches(n);
    # objId = make_child_id(parentId, :studentRanking);
    pWtV = init_endow_pct_weights(switches, st);
    pvec = ParamVector(objId, [pWtV]);
    return EndowPctRanking(objId, pvec, ModelParams.value(pWtV), switches)
end

function init_endow_pct_weights(switches :: EndowPctRankingSwitches{F1},
    st :: SymbolTable) where F1
    wtV = switches.wtV;
    n = length(wtV);
    p = Param(:wtV,
        LatexLH.description(st, :rankWt), latex(st, :rankWt),
        wtV, wtV, switches.lbV, switches.ubV, switches.doCal);
    return p
end

function validate_ranking_switches(e :: EndowPctRankingSwitches{F1}) where F1
    isValid = true
    if !(length(endow_names(e)) == 1 + length(e.wtV))
        @warn "Weights have wrong length: $e"
        isValid = false;
    end
    if !isempty(e.wtV)
        if !all_greater(e.wtV, e.lbV; atol = 1e-5)  ||  any(isinf.(e.wtV))
            @warn "Invalid weights: $(e.wtV)"
            isValid = false;
        end
        if !all(e.ubV .> e.lbV)
            @warn "Invalid bounds: $(e.lbV), $(e.ubV)"
            isValid = false;
        end
    end
    if !(length(e.lbV) == length(e.ubV) == length(e.wtV))
        @warn "Invalid sizes: $(e.lbV), $(e.wtV), $(e.ubV)"
        isValid = false;
    end
    return isValid
end

function validate_ranking(e :: EndowPctRanking{F1}) where F1
    isValid = true;
    isValid = isValid && (size(endow_names(e)) == size(weights(e)))
    return isValid
end

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


text:
"""
Parse fe variables
"""
function parse_fe_variables(fe_vars, expvars; depvar=nothing, is_pair=false)
    valid_vars = copy(expvars)

    if depvar != nothing
        append!(valid_vars, [depvar])
    end

    selected_vars = []

    if is_pair
        vars = []
        if isa(fe_vars, Pair)
            fe_vars = [fe_vars]
        end
        for var in fe_vars
            vars = vcat(vars, [Symbol(var[1])=>var[2]])
            selected_vars = append!(selected_vars, [Symbol(var[1])])
        end
        fe_vars = vars
    else
        if !isa(fe_vars, Array)
            fe_vars = [fe_vars]
        end
        fe_vars = [Symbol(s) for s in fe_vars ]
        selected_vars = fe_vars
    end

    if !ModelSelection.in_vector(selected_vars, valid_vars)
        error(SOME_VARIABLES_NOT_FOUND)
    end

    return fe_vars
end

"""
Adds square feature extraction to data
"""
function data_add_fe_sqr(data, fe_vars)
    postfix = "_sqrt"
    new_vars = []
    existent_vars = []
    for fe_var in fe_vars
        if !(Symbol(string(fe_var, postfix)) in data.expvars)
            append!(new_vars, [fe_var])
        else
            append!(existent_vars, [fe_var])
        end
    end

    if size(new_vars, 1) > 0
        data.expvars_data = hcat(data.expvars_data, (data.expvars_data[:, [ ModelSelection.get_column_index(var, data.expvars) for var in new_vars ]]).^2)
        data.expvars = vcat(data.expvars, [Symbol(string(var, postfix)) for var in new_vars ])
    end

    if size(existent_vars, 1) > 0
        data.expvars_data[:, [ ModelSelection.get_column_index(Symbol(string(var, postfix)), data.expvars) for var in existent_vars ]] = data.expvars_data[:, [ ModelSelection.get_column_index(var, data.expvars) for var in existent_vars ]].^2
    end

    return data
end

"""
Adds log feature extraction to data
"""
function data_add_fe_log(data, fe_vars)
    try
        postfix = "_log"
        new_vars = []
        existent_vars = []
        for fe_var in fe_vars
            if !(Symbol(string(fe_var, postfix)) in data.expvars)
                append!(new_vars, [fe_var])
            else
                append!(existent_vars, [fe_var])
            end
        end

        if size(new_vars, 1) > 0
            data.expvars_data = hcat(data.expvars_data, log.(data.expvars_data[:, [ ModelSelection.get_column_index(var, data.expvars) for var in new_vars ]]))
            data.expvars = vcat(data.expvars, [Symbol(string(var, postfix)) for var in new_vars ])
        end

        if size(existent_vars, 1) > 0
            data.expvars_data[:, [ ModelSelection.get_column_index(Symbol(string(var, postfix)), data.expvars) for var in existent_vars ]] = log.(data.expvars_data[:, [ ModelSelection.get_column_index(var, data.expvars) for var in existent_vars ]])
        end
        return data
    catch
        error(LOG_FUNCTION_ERROR)
    end
end

"""
Adds inverse feature extraction to data
"""
function data_add_fe_inv(data, fe_vars)
    postfix = "_inv"
    new_vars = []
    existent_vars = []
    for fe_var in fe_vars
        if !(Symbol(string(fe_var, postfix)) in data.expvars)
            append!(new_vars, [fe_var])
        else
            append!(existent_vars, [fe_var])
        end
    end

    if size(new_vars, 1) > 0
        data.expvars_data = hcat(data.expvars_data, 1 ./ (data.expvars_data[:, [ ModelSelection.get_column_index(var, data.expvars) for var in new_vars ]]))
        data.expvars = vcat(data.expvars, [Symbol(string(var, postfix)) for var in new_vars ])
    end

    if size(existent_vars, 1) > 0
        data.expvars_data[:, [ ModelSelection.get_column_index(Symbol(string(var, postfix)), data.expvars) for var in existent_vars ]] = 1 ./ (data.expvars_data[:, [ ModelSelection.get_column_index(var, data.expvars) for var in existent_vars ]])
    end
    return data
end

"""
Adds lag feature extraction to data
"""
function data_add_fe_lag(data, fe_vars)
    nobs = size(data.expvars_data, 1)
    postfix = "_l"
    csis = (data.panel != nothing) ? unique(data.panel_data) : [nothing]
    depvar_enabled = false

    for var in fe_vars
        if var[1] in data.expvars
            col = ModelSelection.get_column_index(var[1], data.expvars)
            num_cols = 0
            for i = 1:var[2]
                if !(Symbol(string(var[1], postfix, i)) in data.expvars)
                    num_cols = num_cols + 1
                end
            end
            var_data = Array{Union{Missing, data.datatype}}(missing, nobs, num_cols)
            m = 1
            for i = 1:var[2]
                expvar = Symbol(string(var[1], postfix, i))

                col_added = false
                for csi in csis
                    rows = (csi != nothing) ? findall(x->x == csi, data.panel_data) : collect(1:1:nobs)
                    num_rows = size(rows, 1)
                    if !(expvar in data.expvars)
                        var_data[rows[1]:rows[1]+num_rows-1, m] .= lag(data.expvars_data[rows[1]:rows[1]+num_rows-1, col], i)
                        col_added = true
                    else
                        n = ModelSelection.get_column_index(expvar, data.expvars)
                        data.expvars_data[rows[1]:rows[1]+num_rows-1, n].= lag(data.expvars_data[rows[1]:rows[1]+num_rows-1, col], i)
                    end
                end
                if col_added
                    m = m + 1
                end

                if !(expvar in data.expvars)
                    data.expvars = vcat(data.expvars, [expvar])
                end
            end
        else
            num_cols = 0
            for i = 1:var[2]
                if !(Symbol(string(var[1], postfix, i)) in data.expvars)
                    num_cols = num_cols + 1
                end
            end
            var_data = Array{Union{Missing, data.datatype}}(missing, nobs, num_cols)
            m = 1
            for i = 1:var[2]
                expvar = Symbol(string(var[1], postfix, i))

                col_added = false
                for csi in csis
                    rows = (csi != nothing) ? findall(x->x == csi, data.panel_data) : collect(1:1:nobs)
                    num_rows = size(rows, 1)
                    if !(expvar in data.expvars)
                        var_data[rows[1]:rows[1]+num_rows-1, m] .= lag(data.depvar_data[rows[1]:rows[1]+num_rows-1], i)
                        col_added = true
                    else
                        n = ModelSelection.get_column_index(expvar, data.expvars)
                        data.expvars_data[rows[1]:rows[1]+num_rows-1, n].= lag(data.depvar_data[rows[1]:rows[1]+num_rows-1], i)
                    end
                end
                if col_added
                    m = m + 1
                end

                if !(expvar in data.expvars)
                    data.expvars = vcat(data.expvars, [expvar])
                end
            end
        end
        if size(var_data, 2) > 0
            data.expvars_data = hcat(data.expvars_data, var_data)
        end
    end

    return data
end

"""
Adds interaction between variables
"""
function data_add_interaction(data, interaction)
    if ModelSelection.get_column_index(data.depvar, interaction) != nothing
        error(INTERACTION_DEPVAR_ERROR)
    end

    if !ModelSelection.in_vector(interaction, data.expvars)
        error(INTERACTION_EQUATION_ERROR)
    end

    num_variables = size(interaction, 1)
    infix = "_"

    for i = 1:num_variables-1
        for j = i+1:num_variables
            col = Symbol(string(interaction[i], infix, interaction[j]))
            var_1 = data.expvars_data[:, ModelSelection.get_column_index(interaction[i], data.expvars)]
            var_2 = data.expvars_data[:, ModelSelection.get_column_index(interaction[j], data.expvars)]
            res = var_1 .* var_2
            if !(col in data.expvars)
                data.expvars_data = hcat(data.expvars_data, res)
                data.expvars = vcat(data.expvars, [Symbol(string(interaction[i], infix, interaction[j]))])
            else
                data.expvars_data[:, ModelSelection.get_column_index(col, data.expvars)] = res
            end
        end
    end

    return data
end

"""
Add values to extras
"""
function addextras(data, fe_sqr, fe_log, fe_inv, fe_lag, interaction, removemissings)
    data.extras[ModelSelection.generate_extra_key(FEATUREEXTRACTION_EXTRAKEY, data.extras)] = Dict(
        :fe_sqr => fe_sqr,
        :fe_log => fe_log,
        :fe_inv => fe_inv,
        :fe_lag => fe_lag,
        :interaction => interaction,
        :removemissings => removemissings
    )
    return data
end