text:
typealias Scene Frame{Vehicle}
Scene(n::Int=100) = Frame(Vehicle, n)
Scene(arr::Vector{Vehicle}) = Frame{Vehicle}(arr, length(arr))

Base.show(io, scene::Scene) = print(io, "Scene(with $(length(scene)) cars)")

"""
    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

    ignore = Set{Int}()

    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 && !in(veh.id, 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
                    else
                        push!(ignore, veh.id)
                    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:
#*****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

#supplement to TensorFlow.jl -- remove after TensorFlow.jl update
using PyCall

@pyimport tensorflow as tf
using TensorFlow.CoreTypes

one_hot(indices::Union{AbstractTensor,Void}, depth::Union{AbstractTensor,Void}, on_value::Union{AbstractTensor,Void}=nothing, off_value::Union{AbstractTensor,Void}=nothing, axis::Union{AbstractTensor,Void}=nothing, dtype::Dtype=DT_FLOAT32, name::Union{AbstractString,Void}=nothing) = Tensor(tf.one_hot(;Dict(:indices=>indices, :depth=>depth, :on_value=>on_value, :off_value=>off_value, :axis=>axis, :dtype=>dtype, :name=>name)...))

export one_hot

function cond_(pred::Union{AbstractTensor,Void}, fn1::Function,
    fn2::Function, name::Union{AbstractString,Void}=nothing)

    fun = pyeval("""lambda f: lambda: f()""")
    Tensor(tf.cond(;Dict(:pred=>pred, :fn1=>fun(fn1), :fn2=>fun(fn2), :name=>name)...))
end

export cond_

ext: jl


text:
"""
A CategoricalVector is an AbstractVector which is treated as a categorical axis regardless
of the element type. Duplicate values are not allowed but are not filtered out.

A CategoricalVector axis can be indexed with an ClosedInterval, with a value, or with a
vector of values. Use of a CategoricalVector{Tuple} axis allows indexing similar to the
hierarchical index of the Python Pandas package or the R data.table package.

In general, indexing into a CategoricalVector will be much slower than the corresponding
SortedVector or another sorted axis type, as linear search is required.

### Constructors

```julia
CategoricalVector(x::AbstractVector)
```

### Arguments

* `x::AbstractVector` : the wrapped vector

### Examples

```julia
v = CategoricalVector(collect([1; 8; 10:15]))
A = AxisArray(reshape(1:16, 8, 2), v, [:a, :b])
A[Axis{:row}(1), :]
A[Axis{:row}(10), :]
A[Axis{:row}([1, 10]), :]

## Hierarchical index example with three key levels

data = reshape(1.:40., 20, 2)
v = collect(zip([:a, :b, :c][rand(1:3,20)], [:x,:y][rand(1:2,20)], [:x,:y][rand(1:2,20)]))
A = AxisArray(data, CategoricalVector(v), [:a, :b])
A[:b, :]
A[[:a,:c], :]
A[(:a,:x), :]
A[(:a,:x,:x), :]
```
"""
struct CategoricalVector{T, A<:AbstractVector{T}} <: AbstractVector{T}
    data::A
end

Base.getindex(v::CategoricalVector, idx::Int) = v.data[idx]
Base.getindex(v::CategoricalVector, idx::AbstractVector) = CategoricalVector(v.data[idx])

Base.length(v::CategoricalVector) = length(v.data)
Base.size(v::CategoricalVector) = size(v.data)
Base.size(v::CategoricalVector, i) = size(v.data, i)
Base.axes(v::CategoricalVector) = Base.axes(v.data)

axistrait(::Type{CategoricalVector{T,A}}) where {T,A} = Categorical
checkaxis(::CategoricalVector) = nothing

## Add some special indexing for CategoricalVector{Tuple}'s to achieve something like
## Panda's hierarchical indexing

axisindexes(ax::Axis{S,CategoricalVector{T,A}}, idx) where {T<:Tuple,S,A} = axisindexes(ax, (idx,))

function axisindexes(ax::Axis{S,CategoricalVector{T,A}}, idx::Tuple) where {T<:Tuple,S,A}
    collect(filter(ax_idx->_tuple_matches(ax.val[ax_idx], idx), Base.axes(ax.val)...))
end

function _tuple_matches(element::Tuple, idx::Tuple)
    length(idx) <= length(element) || return false

    for (x, y) in zip(element, idx)
        x == y || return false
    end

    return true
end

axisindexes(ax::Axis{S,CategoricalVector{T,A}}, idx::AbstractArray) where {T<:Tuple,S,A} =
    vcat([axisindexes(ax, i) for i in idx]...)


text:
function time(phifunc, x, y)
    dx = x[2:end] - x[1:(end-1)]
    dy = y[2:end] - y[1:(end-1)]
    dpath = sqrt(dx.*dx + dy.*dy)
    phi = zeros(length(dx))
    C2 = C*C
    @inbounds for i in 1:length(dx)
        xc = x[i] + 0.5*dx[i]
        yc = y[i] + 0.5*dy[i]
        phi[i] = phifunc(xc,yc)./C2
    end
    time_dilation = 1.0-2.*phi
    sum(dpath./time_dilation) ./ C
end

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
# assume 25 Celcius. From PubChem.
# returns solubility in grams / litre.
function solubilityinwater(molecule_name)::Float64
    #

    if Unicode.normalize(molecule_name, casefold = true) == Unicode.normalize("L-Tyrosine", casefold=true)
        # https://pubchem.ncbi.nlm.nih.gov/compound/Tyrosine#section=Melting-Point
        #return 0.453

        # https://www.sigmaaldrich.com/content/dam/sigma-aldrich/docs/Sigma-Aldrich/Product_Information_Sheet/t3754pis.pdf
        return 0.45

    elseif Unicode.normalize(molecule_name, casefold = true) == Unicode.normalize("L-Tryptophan", casefold=true)
        # https://pubchem.ncbi.nlm.nih.gov/compound/Tryptophan#section=Solubility
        return 11.4

    elseif Unicode.normalize(molecule_name, casefold = true) == Unicode.normalize("L-Glutamic acid", casefold=true)
        # https://pubchem.ncbi.nlm.nih.gov/compound/Glutamic-acid#section=Solubility
        return 8570 * 1e-3

    elseif Unicode.normalize(molecule_name, casefold = true) == Unicode.normalize("Caffeine", casefold=true)
        # https://pubchem.ncbi.nlm.nih.gov/compound/Caffeine#section=Solubility
        return (2.16*10+4) * 1e-3
    end

    # TODO
    # 30: Epinephrine
    # solute_concentration = NaN
    # reference_concentration = 0.5
    #
    # 31: (-)-Norepinephrine
    # solute_concentration = NaN
    # reference_concentration = 0.5
    #
    # 32: 3,4-Dihydroxy-L-phenylalanine
    # solute_concentration = NaN
    # reference_concentration = 0.5

    return NaN
end

########### fetch metadata from GISSMO or BMRB.

# look up from JLD.
function getmolarmass(base_dir, target_name::String)
    load_path = joinpath(base_dir, "molar_masses.jld")
    dic = JLD.load(load_path)

    molecule_names = dic["molecule_names"]
    molar_masses = dic["molar_masses"]

    ind = findfirst(xx->xx==target_name, molecule_names)

    if typeof(ind) == Nothing
        return NaN
    end

    return molar_masses[ind]
end

function fetchentrymetadata(entry::String)

    # get molecular weight (natural abundance).
    b = "http://api.bmrb.io/v2/entry/$(entry)"
    HTTP_result = HTTP.request("GET", b)
    json_data = JSON3.read(HTTP_result.body)
    frame_data = json_data[Symbol(entry)][:saveframes]

    ind = findfirst(xx->(xx[:category] == "chem_comp"), frame_data)
    ind2 = findfirst(xx->(xx[1] == "Formula_weight"), frame_data[ind][:tags])
    molar_mass_solute = parse(Float64, frame_data[ind][:tags][ind2][2]) # grams per mol.

    return molar_mass_solute, json_data
end

# was called parseconcentrationfromNMReDATA
function convertsampleconcentration(compound_names,
                                    concentration_value_strings,
                                    concentration_units,
                                    compound_types,
                                    solute_concentration_string,
                                    solute_concentration_unit,
                                    molar_mass_solute::Float64,
                                    solubility_gL_solute::Float64)

    ### case if solute is saturated.
    solute_mM = convertstringtomM(solute_concentration_string,
                                solute_concentration_unit,
                                solubility_gL_solute,
                                molar_mass_solute)

    ### find reference.
    ref_mM = NaN
    ref_name = ""
    molar_mass_ref = NaN
    solubility_gL_ref = NaN

    ind = findfirst(xx->(occursin("Reference", xx) || occursin("reference", xx)), compound_types)
    if length(ind) < 1
        println("Warning, cannot find reference entry in provided sample composition metadata. Assigning NaN as reference concentration.")
        ref_mM = NaN

    else

        # convert reference.
        ref_name = compound_names[ind]
        molar_mass_ref = getmolarmass(ref_name)
        solubility_gL_ref = solubilityinwater(ref_name)

        ref_mM = convertstringtomM(concentration_value_strings[ind],
                                        concentration_units[ind],
                                        solubility_gL_ref,
                                        molar_mass_ref)
    end

    return solute_mM, ref_mM, molar_mass_ref, solubility_gL_ref, ref_name
end

function convertstringtomM( value_string,
                            unit_string,
                            solubility_gL::T,
                            molar_mass::T) where T

    #
    x_mM = NaN

    ret_val = tryparse(T, value_string)

    if typeof(ret_val) != T
        # fall back method: use saturation value.
        x_mM = convertsaturationtomM(solubility_gL, molar_mass)
    else
        x_mM = ret_val
    end

    # if saturated, set unit to a known case of convertconcentrationtomM().
    if occursin("saturated", value_string) || occursin("Saturated", value_string)
        unit_string = "mM"
    end

    ### solute conversion.
    x_mM = convertconcentrationtomM(x_mM, unit_string, molar_mass)

    return x_mM
end

function convertsaturationtomM(solubility_gL::T, molar_mass::T)::T where T
    # # old
    # solubility_g_100mL = solubility_gL/10
    # x_M = solubility_g_100mL/molar_mass # molar concentration.
    # x_mM = x_M*1000 # millimolar concentration.

    x_M = solubility_gL/molar_mass # molar concentration.
    x_mM = x_M*1000 # millimolar concentration.

    return x_mM
end

function convertconcentrationtomM(x::T, unit, molar_mass::T)::T where T <: Real

    if unit == "uM"
        return x * 1e-3

    elseif unit == "%"
        # see http://abacus.bates.edu/~ganderso/biology/resources/dilutions.html#percenttomolarity

        x_M = x*10/molar_mass # molar concentration.
        return x_M*1000 # millimolar concentration.

    elseif unit == "mM"

        return x
    end

    return NaN # unknown case.
end

# assumes temperature is always in Kelvins.
function parseconditionsfromNMReDATA(file_strings::Vector{String})
    # pH.
    header_string = "> <NMREDATA_PH>"

    ind = findfirst(xx->xx==header_string, file_strings)
    @assert typeof(ind) == Int # must not be empty.

    pH = NaN
    val = tryparse(Float64, file_strings[ind+1])
    if typeof(val) == Float64
        pH = val
    end

    # temperature.
    header_string = "> <NMREDATA_TEMPERATURE>"

    ind = findfirst(xx->xx==header_string, file_strings)
    @assert typeof(ind) == Int # must not be empty.

    s = file_strings[ind+1]
    tokens = split(s, ' ')

    temperature = NaN
    val = tryparse(Float64, tokens[1])
    if typeof(val) == Float64
        temperature = val
    end

    return pH, temperature
end

function parsecompositionfromNMReDATA(file_strings::Vector{String})

    ### find the solute concentration.
    header_string = "> <NMREDATA_CONCENTRATION>"

    ind = findfirst(xx->xx==header_string, file_strings)
    @assert typeof(ind) == Int # must not be empty.

    # split and parse.
    s = file_strings[ind+1]
    tokens = split(s, ' ')

    ### solute concentration. Possibly numerical (string) value or "saturated".
    # current strategy: when in doubt, use 100%.
    solute_concentration_unit = tokens[2]
    solute_concentration_string = tokens[1]

    ### find the mixture composotion info
    header_string = "> <NMREDATA_SOLVENT>"

    ind = findfirst(xx->xx==header_string, file_strings)
    @assert typeof(ind) == Int # must not be empty.

    # split and parse.
    s = file_strings[ind+1]
    tokens = split(s, '/')
    tmp = tokens[1:end-1]

    tokens = split(tokens[end], ' ')

    compound_names = [tmp; tokens[1]]

    concentration_values = split(tokens[2], ':')

    concentration_units = split(tokens[3], ':')
    compound_types = split(tokens[4], ':')
    @assert length(compound_names) == length(concentration_values) == length(concentration_units) == length(compound_types)

    return compound_names, concentration_values, concentration_units, compound_types, solute_concentration_string, solute_concentration_unit
end

### convert concentration values to millimolar.
function converttomM(   compound_types,
                        concentration_units,
                        concentration_values,
                        molar_mass_solute::T,
                        molar_mass_solvent::T,
                        molar_mass_ref::T) where T

    concentration_values_mM = collect( parse(T, concentration_values[i]) for i = 1:length(concentration_values))

    for i = 1:length(compound_types)
        cv = concentration_values_mM[i]

        molar_mass = molar_mass_solute
        if compound_types[i] == "Solvent"
            molar_mass = molar_mass_solvent

        elseif compound_types[i] == "Reference"
            molar_mass = molar_mass_ref
        end

        if compound_types[i] == "Reference" || compound_types[i] == "Solvent"

            if concentration_units[i] == "uM"
                concentration_values_mM[i] = cv/1000
            elseif concentration_units[i] == "%"
                # see http://abacus.bates.edu/~ganderso/biology/resources/dilutions.html#percenttomolarity

                concentration_values_mM[i] = cv*10/molar_mass # molar concentration.
                concentration_values_mM[i] = concentration_values_mM[i]*1000 # millimolar concentration.
            end
        else
            concentration_values_mM[i] = NaN
        end
    end

    return concentration_values_mM
end


text:
using Ipopt
using ForwardDiff

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

function eval_f(x)

  VM3 = x[1]
  VA3 = x[2]
  P1  = x[4]

  VM1 = x[5]
  P2 = x[6]

  VA1 = p[1]

  # we fix generation weights inside the
  # function to simplify the script and
  # follow the paper closely.
  w1 = 1.0
  w2 = 1.0

  cost = abs(P1 + P2)

  return cost
end

function eval_g(x, g)
  # Get Float64 of Vectors{Float64}, that is the first parameter

  # retrieve variables
  VM3 = x[1]
  VA3 = x[2]
  VA2 = x[3]
  P1  = x[4]

  VM1 = x[5]
  P2 = x[6]
  VM2 = x[7]

  VA1 = p[1]
  P3 = p[2]
  Q3 = p[3]

  # intermediate quantities
  VA23 = VA2 - VA3
  VA31 = VA3 - VA1
  VA32 = VA3 - VA2
  VA13 = VA1 - VA3

  g[1] = 4.0*VM2*VM2 + VM2*VM3*(-4*cos(VA23) + 10*sin(VA23)) - P2
  g[2] = (8.0*VM3*VM3 + VM3*VM1*(-4*cos(VA31) + 5*sin(VA31))
          + VM3*VM2*(-4*cos(VA32) + 10*sin(VA32)) + P3)
  g[3] = (15.0*VM3*VM3 + VM3*VM1*(-4*sin(VA31) - 5*cos(VA31))
          + VM3*VM2*(-4*sin(VA32) - 10*cos(VA32)) + Q3)
  g[4] = 4.0*VM1*VM1 + VM1*VM3*(-4*cos(VA13) + 5*sin(VA13)) - P1
  return g
end

function eval_g(x)
  T = typeof(x)
  g = zeros(T.parameters[1], m)
  eval_g(x,g)
  return g
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
  # Bad: grad_f    = zeros(4)  # Allocates new array
  # OK:  grad_f[:] = zeros(4)  # Modifies 'in place'
  ForwardDiff.gradient!(grad_f, eval_f, x)
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
  if mode == :Structure
    idx = 1
    for c in 1:m #number of constraints
        for i in 1:n # number of variables
            rows[idx] = c ; cols[idx] = i
            idx += 1
        end
    end
  else
    jac = ForwardDiff.jacobian(eval_g, x)
    k = 1
    for i in 1:m
        for j in 1:n
            values[k] = jac[i,j]
            k += 1
        end
    end
  end
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
  if mode == :Structure
    # Symmetric matrix, fill the lower left triangle only
    idx = 1
    for row = 1:6
      for col = 1:row
        rows[idx] = row
        cols[idx] = col
        idx += 1
      end
    end
  else
    hf = ForwardDiff.hessian(eval_f, x)
    k = 1
    for i in 1:n
        for j in 1:i
            values[k] = obj_factor * hf[i,j]
            k += 1
        end
    end
    # hg = ForwardDiff.hessian(eval_g, x)
    h = x -> ForwardDiff.jacobian(x -> ForwardDiff.jacobian(eval_g, x), x)
    hg = h(x)
    hg = reshape(hg, m, n, n)
    for l in 1:m
        k = 1
        for i in 1:n
            for j in 1:i
            values[k] += lambda[l] * hg[l,i,j]
            k += 1
            end
        end
    end
  end
end

n = 7
x_L = [-1e18,-1e18,-1e18,-1e18,-1e18,-1e18,-1e18]
x_U = [1e18,1e18,1e18,1e18,1e18,1e18,1e18]

m = 4
g_L = [0.0, 0.0, 0.0,0.0]
g_U = [0.0, 0.0, 0.0,0.0]

p = zeros(3)
p[1] = 0.0 #VA1, slack angle
p[2] = 2.0 #P3
p[3] = 1.0 #Q3

# Number of nonzeros in upper triangular Hessian
hnnz = Int(n*(n+1)/2)
prob = createProblem(n, x_L, x_U, m, g_L, g_U, m*n, hnnz,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)

prob.x = [1.0, 0.0, 0.0, 1.5, 1.0, 1.7, 1.0]

# This tests callbacks.
function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)
  return iter_count < 100  # Interrupts after one iteration.
end

setIntermediateCallback(prob, intermediate)

solvestat = solveProblem(prob)


text:
### RoadGraphs.jl

"""
The following Julia code consists of separate files that are to be bundled together and installed as a module (`]pkg dev .` in REPL)

Project.toml
src/
    RoadGraphs.jl
    road_model.jl
    road_constraints.jl
    road_utils.jl
"""
module RoadGraphs

include("road_model.jl")

end


text:
function equispaced_rectangular_mesh(
    domain::RectangleDomain{CT},
    n1,
    n2,
) where {FT <: AbstractFloat, CT <: Geometry.Abstract2DPoint{FT}}
    x1x2min, x1x2max = domain.x1x2min, domain.x1x2max
    x1max = Geometry.component(x1x2max, 1)
    x2max = Geometry.component(x1x2max, 2)
    x1min = Geometry.component(x1x2min, 1)
    x2min = Geometry.component(x1x2min, 2)
    x1c = range(x1min, x1max; length = n1 + 1)
    x2c = range(x2min, x2max; length = n2 + 1)
    return Mesh2D(domain, x1c, x2c)
end

"""
    rectangular_mesh(x1c, x2c, per = (false, false))

This function builds a 2D rectangular mesh with points located at x1c and x2c in x1 and x2 directions respectively.
"""
function Mesh2D(
    domain::RectangleDomain{CT},
    x1c,
    x2c,
) where {FT <: AbstractFloat, CT <: Geometry.Abstract2DPoint{FT}}
    per = tuple(domain.x1boundary === nothing, domain.x2boundary === nothing)
    I = Int
    nx1, nx2 = length(x1c), length(x2c)
    nel1, nel2 = nx1 - 1, nx2 - 1

    nverts = nx1 * nx2
    nfaces = nx1 * (nx2 - 1) + (nx1 - 1) * nx2
    nelems = (nx1 - 1) * (nx2 - 1)
    nbndry = 4

    nfc1 = nx1 * (nx2 - 1) # faces with normals along x1 direction

    ndmat = reshape(1:nverts, nx1, nx2)

    fcmat1 = reshape(1:nfc1, nx1, nx2 - 1)        # faces with normals along x1 direction
    fcmat2 = reshape((nfc1 + 1):nfaces, nx1 - 1, nx2) # faces with normals along x2 direction

    emat = reshape(1:nelems, nx1 - 1, nx2 - 1)

    coordinates = [repeat(x1c, 1, nx2)[:] repeat(x2c', nx1, 1)[:]] # node coordinates

    face_verts = hcat(
        vcat(ndmat[:, 1:(nx2 - 1)][:], ndmat[1:(nx1 - 1), :][:]), # face nodes
        vcat(ndmat[:, 2:nx2][:], ndmat[2:nx1, :][:]),
    )

    face_boundary = Vector{I}(zeros(nfaces))
    # accounting for periodic boundaries, wherever applicable
    periodic = per[1] || per[2]
    periodic_faces = Dict{Int, Int}()
    global_vert = Vector{I}(1:nverts)
    if per[1]
        bdy1 = emat[end:end, :]
        bdy2 = emat[1:1, :]
        nbndry -= 2
        for pfc in 1:(nx2 - 1)
            pfc1, pfc2 = fcmat1[1, pfc], fcmat1[nx1, pfc]
            periodic_faces[pfc1] = fcmat1[pfc2]
            periodic_faces[pfc2] = fcmat1[pfc1]
        end
        for i in 1:nx2
            global_vert[ndmat[nx1, i]] = ndmat[1, i]
        end
    else
        bdy1 = bdy2 = zeros(I, 1, nx2 - 1)
        face_boundary[fcmat1[1, :]] .= 1  # left boundary
        face_boundary[fcmat1[end, :]] .= 2  # right boundary
    end
    if per[2]
        bdy3 = emat[:, end:end]
        bdy4 = emat[:, 1:1]
        nbndry -= 2
        for pfc in 1:(nx1 - 1)
            periodic_faces[fcmat2[pfc, 1]] = fcmat2[pfc, nx2]
            periodic_faces[fcmat2[pfc, nx2]] = fcmat2[pfc, 1]
        end
        for i in 1:nx1
            global_vert[ndmat[i, nx2]] = global_vert[ndmat[i, 1]]
        end
    else
        bdy3 = bdy4 = zeros(I, nx1 - 1, 1)
        face_boundary[fcmat2[:, 1]] .= 3  # bottom boundary
        face_boundary[fcmat2[:, end]] .= 4  # top boundary
    end

    face_neighbors = hcat(
        vcat(vcat(bdy1, emat)[:], hcat(bdy3, emat)[:]),
        zeros(I, nfaces),
        vcat(vcat(emat, bdy2)[:], hcat(emat, bdy4)[:]),
        zeros(I, nfaces),
        ones(I, nfaces),
    )

    elem_verts = hcat(
        ndmat[1:(nx1 - 1), 1:(nx2 - 1)][:], # node numbers
        ndmat[2:nx1, 1:(nx2 - 1)][:], # for each element
        ndmat[1:(nx1 - 1), 2:nx2][:],
        ndmat[2:nx1, 2:nx2][:],
    )
    elem_faces = hcat(
        fcmat1[1:(nx1 - 1), :][:], # face numbers for
        fcmat1[2:nx1, :][:],   # each element
        fcmat2[:, 1:(nx2 - 1)][:],
        fcmat2[:, 2:nx2][:],
    )

    for fc in 1:nfaces
        elems = (face_neighbors[fc, 1], face_neighbors[fc, 3])

        if periodic && haskey(periodic_faces, fc)
            shadow = periodic_faces[fc]
            localface1 = findfirst(elem_faces[elems[1], :] .== fc)
            if !isnothing(localface1)
                localface2 = findfirst(elem_faces[elems[2], :] .== shadow)
                if !isnothing(localface2)
                    face_neighbors[fc, 2] = localface1
                    face_neighbors[fc, 4] = localface2
                else
                    error(
                        "rectangular_mesh: Fatal error, face could not be located in neighboring element",
                    )
                end
            else
                localface1 = findfirst(elem_faces[elems[1], :] .== shadow)
                localface2 = findfirst(elem_faces[elems[2], :] .== fc)
                if !isnothing(localface1) && !isnothing(localface2)
                    face_neighbors[fc, 2] = localface1
                    face_neighbors[fc, 4] = localface2
                else
                    error(
                        "rectangular_mesh: Fatal error, face could not be located in neighboring element",
                    )
                end
            end
        else # no periodic boundaries
            for e in 1:2
                el = elems[e]
                if el ≠ 0
                    localface = findfirst(elem_faces[el, :] .== fc)
                    if isnothing(localface)
                        error(
                            "rectangular_mesh: Fatal error, face could not be located in neighboring element",
                        )
                    else
                        face_neighbors[fc, 2 + (e - 1) * 2] = localface
                    end
                end
            end
        end # note: no orientation mismatch in this case
    end

    if periodic # marking shadow faces with -1
        for (per, shd) in periodic_faces
            face_boundary[max(per, shd)] = -1
        end
    end
    boundary_tags = sort(unique(face_boundary))
    face_boundary_offset = ones(I, length(boundary_tags) + 1)
    face_boundary_tags = sort(face_boundary)
    face_boundary = sortperm(face_boundary)

    if periodic # remove redundant shadow faces
        loc = findlast(face_boundary_tags .== -1) + 1
        boundary_tags = boundary_tags[2:end]
        face_boundary_offset = face_boundary_offset[2:end]
        face_boundary_tags = face_boundary_tags[loc:end]
        face_boundary = face_boundary[loc:end]
    end

    for i in 1:(length(boundary_tags) - 1)
        tag = boundary_tags[i + 1]
        face_boundary_offset[i + 1] = findfirst(face_boundary_tags .== tag)
    end
    face_boundary_offset[end] = length(face_boundary) + 1
    tag_names = (:interior, :west, :east, :south, :north)

    boundary_tag_names = tuple(tag_names[boundary_tags .+ 1]...)
    NB = length(boundary_tag_names)

    # add unique vertex iterator information
    vtconn = map(i -> zeros(I, i), zeros(I, nverts))

    for el in 1:nelems
        for lv in 1:4
            vt = global_vert[elem_verts[el, lv]]
            push!(vtconn[vt], el, lv)
        end
    end
    unique_verts = I[]
    uverts_conn = I[]
    uverts_offset = I.([1])
    for vt in 1:nverts
        lconn = length(vtconn[vt])
        if lconn > 0
            push!(unique_verts, vt)
            push!(uverts_conn, vtconn[vt]...)
            push!(uverts_offset, uverts_offset[end] + lconn)
        end
    end

    return Mesh2D(
        domain,
        NoWarp(),
        nverts,
        nfaces,
        nelems,
        nbndry,
        coordinates,
        unique_verts,
        uverts_conn,
        uverts_offset,
        face_verts,
        face_neighbors,
        face_boundary,
        boundary_tags,
        boundary_tag_names,
        face_boundary_offset,
        elem_verts,
        elem_faces,
    )
end


text:
## extensions/FrameTools.jl --- tools for describing and manipulating frames

export getchaindict
export addprotons!

"""
    getchaindict(frame::StructureFrame)

Get a dictionary of chains in the `frame` containing `UnitRange`s of standard
protein residues.
"""
function getchaindict(fr::StructureFrame)
    dct = Dict{Symbol,UnitRange{UInt}}()

    current_ch = :o
    start = stop = one(UInt)
    stopped = true

    (; res_list) = fr
    for i in eachindex(res_list)
        res = @inbounds res_list[i]

        chid_ = res.chainid
        pdb_  = res.standard_pdb

        @assert !isnothing(pdb_) && !ismissing(chid_) "cannot define chaindict"

        ## start condition: when previously stopped, but encounter standard residue
        if stopped
            !pdb_ && continue
            stop = i
            current_ch = Symbol(chid_)
            stopped = false
        end

        ## stop condition: when on diff chain or not standard residue (PDB)
        if !(pdb_ && Symbol(chid_) == current_ch)
            stopped = true
            stop = i-1
            push!(dct, current_ch=>(start:stop))
            ## restart start index
            start = i
        end
    end

    return dct
end

"""
    addprotons!(fr::StructureFrame[; has_na])

Add backbone protons to residues in the Frame `fr`. Estimate the position of the
missing proton using the method used by DSSP (within its source code).
Optional keyword argument `has_na` indicates whether the frame `fr` has nucleic
acids added in parsing. Defaults to `false`.

**Note**: Setting `has_na` is important, when it can be shown that the frame has
nucleic acids.
"""
function addprotons!(fr::StructureFrame; has_na=false)
    (; at_pos, at_list, res_list) = fr                ## unpack these
    len = length(at_pos)

    ## drop first element in iteration
    for i in Iterators.drop(eachindex(res_list), 1)
        res = @inbounds res_list[i]
        ## don't do anything if it's not a standard pdb residue (protein)
        # N terminal end or proline residue; skip
        nucleic_acid = has_na && res.name ∉ kAminoAcidNames
        (!res.standard_pdb || res.name == :PRO || nucleic_acid) && continue

        ## calculate new position and other shit for proton
        Nid = getatom(res, :N)          ## amide nitrogen

        ## hydrogen assignment 1) position; 2) name; 3) charge; 4) type/mass
        Hpos = at_pos[Nid]  # initially set position to amide N

        # carbonyl of previous residue; estimate H from that
        # H must be opposite O and peptide bond is estimated planar
        res_prev = @inbounds res_list[i-1]
        res_prev.standard_pdb || continue

        Oprevid = getatom(res_prev, :O); Cprevid = getatom(res_prev, :C)
        _opos = at_pos[Oprevid]; _cpos = at_pos[Cprevid]
        codist = distance(_opos, _cpos)
        Hpos += (_cpos - _opos) / codist

        ## make the actual atom; i have no idea about charge tbh set to 0
        ## DSSP has its own constants for calculating Hbond energy
        Hatom = JAtom(:H, "H", 1.008, float(0))
        push!(at_list, Hatom); push!(at_pos, Hpos)
        Hidx = length(at_pos)
        setatom!(res, Hatom, Hidx)
    end

    nothing
end


text:
using TimeZones.TZData: parse_components
using TimeZones: Transition

dt = DateTime(1942,12,25,1,23,45)
custom_dt = DateTime(1800,1,1)

utc = FixedTimeZone("UTC")
gmt = FixedTimeZone("GMT", 0)
foo = FixedTimeZone("FOO", 0)
null = FixedTimeZone("", 10800)
fixed = FixedTimeZone("UTC+01:00")
est = FixedTimeZone("EST", -18000)
warsaw = first(compile("Europe/Warsaw", tzdata["europe"]))
apia = cache_tz(compile("Pacific/Apia", tzdata["australasia"]))
honolulu = cache_tz(compile("Pacific/Honolulu", tzdata["northamerica"]))  # Uses cutoff
ulyanovsk = first(compile("Europe/Ulyanovsk", tzdata["europe"]))  # No named abbreviations
new_york = first(compile("America/New_York", tzdata["northamerica"]))  # Underscore in name
custom = VariableTimeZone("Test/Custom", [Transition(custom_dt, utc)])  # Non-cached variable time zone

@test sprint(print, utc) == "UTC"
@test sprint(print, gmt) == "GMT"
@test sprint(print, foo) == "FOO"
@test sprint(print, null) == "UTC+03:00"
@test sprint(print, fixed) == "UTC+01:00"
@test sprint(print, est) == "EST"
@test sprint(print, warsaw) == "Europe/Warsaw"
@test sprint(print, apia) == "Pacific/Apia"
@test sprint(print, honolulu) == "Pacific/Honolulu"
@test sprint(print, ulyanovsk) == "Europe/Ulyanovsk"
@test sprint(print, custom) == "Test/Custom"

@test sprint(show_compact, utc) == "tz\"UTC\""
@test sprint(show_compact, gmt) == "tz\"GMT\""
@test sprint(show_compact, foo) == "FixedTimeZone(\"FOO\", 0)"
@test sprint(show_compact, null) == "FixedTimeZone(\"\", 10800)"
@test sprint(show_compact, fixed) == "tz\"UTC+01:00\""
@test sprint(show_compact, est) == "tz\"EST\""
@test sprint(show_compact, warsaw) == "tz\"Europe/Warsaw\""
@test sprint(show_compact, apia) == "tz\"Pacific/Apia\""
@test sprint(show_compact, honolulu) == "tz\"Pacific/Honolulu\""
@test sprint(show_compact, ulyanovsk) == "tz\"Europe/Ulyanovsk\""
@test sprint(show_compact, custom) == "VariableTimeZone(\"Test/Custom\", ...)"

@test sprint(show, utc) == "tz\"UTC\""
@test sprint(show, gmt) == "tz\"GMT\""
@test sprint(show, foo) == "FixedTimeZone(\"FOO\", 0)"
@test sprint(show, null) == "FixedTimeZone(\"\", 10800)"
@test sprint(show, fixed) == "tz\"UTC+01:00\""
@test sprint(show, est) == "tz\"EST\""
@test sprint(show, warsaw) == "tz\"Europe/Warsaw\""
@test sprint(show, apia) == "tz\"Pacific/Apia\""
@test sprint(show, honolulu) == "tz\"Pacific/Honolulu\""
@test sprint(show, ulyanovsk) == "tz\"Europe/Ulyanovsk\""
@test sprint(show, custom) == "VariableTimeZone(\"Test/Custom\", Transition[Transition($(repr(custom_dt)), tz\"UTC\")], nothing)"

@test sprint(show, MIME("text/plain"), utc) == "UTC"
@test sprint(show, MIME("text/plain"), gmt) == "GMT"
@test sprint(show, MIME("text/plain"), foo) == "FOO (UTC+0)"
@test sprint(show, MIME("text/plain"), null) == "UTC+03:00"
@test sprint(show, MIME("text/plain"), fixed) == "UTC+01:00"
@test sprint(show, MIME("text/plain"), est) == "EST (UTC-5)"
@test sprint(show, MIME("text/plain"), warsaw) == "Europe/Warsaw (UTC+1/UTC+2)"
@test sprint(show, MIME("text/plain"), apia) == "Pacific/Apia (UTC+13/UTC+14)"
@test sprint(show, MIME("text/plain"), honolulu) == "Pacific/Honolulu (UTC-10)"
@test sprint(show, MIME("text/plain"), ulyanovsk) == "Europe/Ulyanovsk (UTC+4)"
@test sprint(show, MIME("text/plain"), custom) == "Test/Custom (UTC+0)"

# ZonedDateTime as a string
zdt = ZonedDateTime(dt, warsaw)
@test string(zdt) == "1942-12-25T01:23:45+01:00"
@test sprint(show_compact, zdt) == "1942-12-25T01:23:45+01:00"
@test sprint(show, zdt) == "ZonedDateTime(1942, 12, 25, 1, 23, 45, tz\"Europe/Warsaw\")"
@test sprint(show, MIME("text/plain"), zdt) == "1942-12-25T01:23:45+01:00"

# TimeZone parsing

# Make sure that TimeZone conversion specifiers are set (issue #24)
CONVERSION_SPECIFIERS = isdefined(Dates, :SLOT_RULE) ? Dates.keys(Dates.SLOT_RULE) : keys(Dates.CONVERSION_SPECIFIERS)
@test 'z' in CONVERSION_SPECIFIERS
@test 'Z' in CONVERSION_SPECIFIERS

df = Dates.DateFormat("z")
@test parse_components("+0100", df) == Any[FixedTimeZone("+01:00")]
@test_throws ArgumentError parse_components("EST", df)
@test_throws ArgumentError parse_components("UTC", df)
@test_throws ArgumentError parse_components("Europe/Warsaw", df)

df = Dates.DateFormat("Z")
@test_throws ArgumentError parse_components("+0100", df)
@test_throws ArgumentError parse_components("EST", df)
@test parse_components("UTC", df) == Any[FixedTimeZone("UTC")]
@test parse_components("Europe/Warsaw", df) == Any[warsaw]

# ZonedDateTime parsing.
# Note: uses compiled time zone information. If these tests are failing try to rebuild
# the TimeZones package.
@test ZonedDateTime("1942-12-25T01:23:45.000+01:00") == ZonedDateTime(dt, fixed)
@test ZonedDateTime("1942-12-25T01:23:45+0100", "yyyy-mm-ddTHH:MM:SSzzz") == ZonedDateTime(dt, fixed)
@test ZonedDateTime("1942-12-25T01:23:45 Europe/Warsaw", "yyyy-mm-ddTHH:MM:SS ZZZ") == ZonedDateTime(dt, warsaw)
@test ZonedDateTime("1942-12-25T01:23:45 America/New_York", "yyyy-mm-ddTHH:MM:SS ZZZ") == ZonedDateTime(dt, new_york)

x = "1942-12-25T01:23:45.123+01:00"
@test string(ZonedDateTime(x)) == x

# Note: CET here represents the FixedTimeZone used in Europe/Warsaw and not the
# VariableTimeZone CET.
@test_throws ArgumentError ZonedDateTime("1942-12-25T01:23:45 CET", "yyyy-mm-ddTHH:MM:SS ZZZ")

# Creating a ZonedDateTime requires a TimeZone to be present.
@test_throws ArgumentError ZonedDateTime("1942-12-25T01:23:45", "yyyy-mm-ddTHH:MM:SSzzz")

# ZonedDateTime formatting
f = "yyyy/m/d H:M:S ZZZ"
@test Dates.format(ZonedDateTime(dt, fixed), f) == "1942/12/25 1:23:45 UTC+01:00"
@test Dates.format(ZonedDateTime(dt, warsaw), f) == "1942/12/25 1:23:45 CET"
@test Dates.format(ZonedDateTime(dt, ulyanovsk), f) == "1942/12/25 1:23:45 UTC+04:00"

f = "yyyy/m/d H:M:S zzz"
@test Dates.format(ZonedDateTime(dt, fixed), f) == "1942/12/25 1:23:45 +01:00"
@test Dates.format(ZonedDateTime(dt, warsaw), f) == "1942/12/25 1:23:45 +01:00"
@test Dates.format(ZonedDateTime(dt, ulyanovsk), f) == "1942/12/25 1:23:45 +04:00"

# The "Z" slot displays the time zone abbreviation for VariableTimeZones. It is fine to use
# the abbreviation for display purposes but not fine for parsing. This means that we
# currently cannot parse all strings produced by format.
df = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS ZZZ")
zdt = ZonedDateTime(dt, warsaw)
@test_throws ArgumentError parse(ZonedDateTime, Dates.format(zdt, df), df)

# Displaying VariableTimeZone's vector of transitions on the REPL has a special printing
@testset "Transitions I/O" begin
    transitions = honolulu.transitions  # Only contains a few transitions
    t = transitions[2]  # 1896-01-13T22:31:26 UTC-10:30/+0 (HST)

    @testset "basic" begin
        @test sprint(print, t) == "1896-01-13T22:31:26 UTC-10:30/+0 (HST)"
        @test sprint(show_compact, t) == "1896-01-13T22:31:26 UTC-10:30/+0 (HST)"
        @test sprint(show, t) == "Transition($(repr(t.utc_datetime)), FixedTimeZone(\"HST\", -37800))"
        @test sprint(show, MIME("text/plain"), t) == "1896-01-13T22:31:26 UTC-10:30/+0 (HST)"
    end

    @testset "REPL vector" begin
        expected_full = string(
            Transition,
            "[",
            join(map(t -> sprint(show, t, context=:compact => false), transitions), ", "),
            "]",
        )

        # Note: The output here is different from the interactive REPL but is representative
        # of the output.
        expected_repl = string(
            Transition,
            "[",
            join(map(t -> sprint(show, t, context=:compact => true), transitions), ", "),
            "]",
        )

        @test sprint(show, transitions, context=:compact => false) == expected_full
        @test sprint(show, transitions; context=:compact_el => true) == expected_repl
    end
end


text:
using NestedTuples

export lazymerge

"""
    lazymerge(x::NamedTuple, y::NamedTuple)

Create a `LazyMerge` struct that behaves like a recursively merged NamedTuple,
but is much faster to construct.

In addition to the usual `getproperty`, this structure supports `get`ting by
`Val` type:
```
(lm::LazyMerge).a == get(lm, Val(:a))
```

In some situatiuons, the latter can be much faster.

The original use case for this is to support using named tuples as namespaces,
especially in the context of probabilistic programming. There, it's common to
have one (possibly nested) named tuple for observed data, and another for a
proposal in an MCMC algorithm. The merge is therefore in the body of a loop
that's executed many times.
"""
function lazymerge(x, y)
    return LazyMerge(x,y)
end

struct LazyMerge{Nx,Ny,Tx,Ty}
    x::NamedTuple{Nx,Tx}
    y::NamedTuple{Ny,Ty}
end

_getx(m::LazyMerge) = getfield(m, :x)
_gety(m::LazyMerge) = getfield(m, :y)

import Base

function Base.show(io::IO, m::LazyMerge)
    io = IOContext(io, :compact => true)
    print(io, "LazyMerge(")
    print(io, _getx(m))
    print(io, ", ")
    print(io, _gety(m))
    print(io, ")")
end

function Base.getproperty(m::LazyMerge, k::Symbol)
    return get(m, Val(k))
end

function Base.get(m::LazyMerge, Val_k::Val{k}) where {k}
    return _get(m, Val_k)
end

Base.propertynames(m::LazyMerge) = union(propertynames(_getx(m)), propertynames(_gety(m)))

_get_code(tx::Missing, ty::Missing, k) = :(error("type LazyMerge has no field ", k))

_get_code(tx, ty::Missing, k) = :(getproperty(_getx(m), k))
_get_code(tx::Missing, ty, k) = :(getproperty(_gety(m), k))

function _get_code(tx::NamedTuple, ty::NamedTuple, k)
    quote
        xk = getproperty(_getx(m), k)
        yk = getproperty(_gety(m), k)
        LazyMerge(xk, yk)
    end
end

# Other than for NamedTuples, `y` gets priority if both match
# Note that this may need to be updated, e.g. for merging arrays with missing values
_get_code(tx, ty, k) = :(getproperty(_gety(m), k))

@generated function _get(m::LazyMerge{Nx,Ny,Tx,Ty}, ::Val{k}) where {Nx,Ny,Tx,Ty,k}
    x = schema(NamedTuple{Nx,Tx})
    y = schema(NamedTuple{Ny,Ty})

    tx = if (k ∈ propertynames(x)) getproperty(x, k) else missing end
    ty = if (k ∈ propertynames(y)) getproperty(y, k) else missing end

    _get_code(tx, ty, k)
end

Base.convert(::Type{NamedTuple}, lm::LazyMerge{Nx,Ny}) where {Nx, Ny} = _convert(NamedTuple, lm)

@gg function _convert(::Type{NamedTuple}, lm::LazyMerge{Nx,Ny}) where {Nx, Ny}
    A = tuple(setdiff(Nx, Ny)...)
    B = tuple(setdiff(Ny, Nx)...)
    dupvars = tuple((Nx ∩ Ny)...)

    quote
        x = _getx(lm)
        y = _gety(lm)
        nt = NamedTuple{$A}(x)
        nt = merge(nt, NamedTuple{$B}(y))
        dup = NamedTuple{$dupvars}(map(v -> getproperty(lm, v), $dupvars))
        merge(nt, dup)
    end
end

# x = (a = (b = 1, c = 2), f = (g = 3, h = 4), g = 2)
# y = (a = (d = 5,), e = (g = 6, j = 7), g = 3)
# m = LazyMerge(x,y)

# m.a
# m.e
# m.g
# m.f
# m.d


text:
remove_tracker(x) = x

_match_container(y, y_el::Number) = first(y)
_match_container(y, y_el) = y


text:
using Random
using Test
using MPI
using ClimateMachine.Mesh.BrickMesh

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    Random.seed!(1234)
    d = 4
    A = rand(1:10, d, 3rank)
    B = BrickMesh.parallelsortcolumns(comm, A, rev = true)

    root = 0
    Acounts = MPI.Gather(Cint(length(A)), root, comm)
    A_all = MPI.Gatherv!(
        A,
        MPI.Comm_rank(comm) == root ?
        VBuffer(similar(A, sum(Acounts)), Acounts) : nothing,
        root,
        comm,
    )

    Bcounts = MPI.Gather(Cint(length(B)), root, comm)
    B_all = MPI.Gatherv!(
        B,
        MPI.Comm_rank(comm) == root ?
        VBuffer(similar(B, sum(Bcounts)), Bcounts) : nothing,
        root,
        comm,
    )

    if MPI.Comm_rank(comm) == root
        A_all = reshape(A_all, d, div(length(A_all), d))
        B_all = reshape(B_all, d, div(length(B_all), d))

        A_all = sortslices(A_all, dims = 2, rev = true)

        @test A_all == B_all
    end
end

main()


text:
#######################################################################
# EffectiveMatrix
#######################################################################
"""
    EffectiveMatrix

A dense matrix of the form `matrix = [ZaL(ω) L' 0; L*ZaR(ω) I*ω-H R*ZrL(ω); 0 R' ZrR(ω)]`,
where H, L and R are given at construction, while ω, φL, φR are given at runtime.
"""
struct EffectiveMatrix{T}
    matrix::Matrix{T}
    L::Matrix{T}
    R::Matrix{T}
    n::Int
    r::Int
end

function EffectiveMatrix(H, L, R)
    size(L) == size(R) || throw("L and R in Effective matrix should have same `size`, but are $(size(L)) and $(size(R))")
    Hd = Matrix(H)
    Ld = Matrix(L)
    Rd = Matrix(R)
    matrix = [0I Ld' 0I; 0Ld -Hd 0Rd; 0I Rd' 0I]
    n, r = size(R)
    return EffectiveMatrix(matrix, L, R, n, r)
end

function effective_matrix!(mat, e::EffectiveMatrix, ω, (ZrL, ZrR, ZaL, ZaR))
    copy!(mat, e.matrix)
    for n in axes(mat, 1)
        mat[n, n] += ω
    end
    r, n = e.r, e.n
    i1, i2, i3 = 1:r, r+1:r+n, r+n+1:n+2r
    copy!(view(mat, i1, i1), ZaL)
    copy!(view(mat, i3, i3), ZrR)
    mul!(view(mat, i2, i1), e.L, ZaR)
    mul!(view(mat, i2, i3), e.R, ZrL)
    return mat
end

pad_effective(m::StridedMatrix{T}, e::EffectiveMatrix, ::Val{:L}) where {T} =
    [zeros(T, e.r, size(m, 2)); m]
pad_effective(m::StridedMatrix{T}, e::EffectiveMatrix, ::Val{:R}) where {T} =
    [m; zeros(T, e.r, size(m, 2))]
pad_effective(m::StridedMatrix{T}, e::EffectiveMatrix, ::Val{:LR}) where {T} =
    [zeros(T, e.r, size(m, 2)); m; zeros(T, e.r, size(m, 2))]

unpad_effective(m::StridedMatrix, e::EffectiveMatrix, ::Val{:L}) = view(m, e.r+1:e.r+e.n, :)
unpad_effective(m::StridedMatrix, e::EffectiveMatrix, ::Val{:R}) = view(m, 1:e.n, :)
unpad_effective(m::StridedMatrix, e::EffectiveMatrix, ::Val{:LR}) = view(m, e.r+1:e.r+e.n, :)

padded_lu!(invmat, e::EffectiveMatrix, ::Val{:L}) = lu!(view(invmat, 1:e.r+e.n, 1:e.r+e.n))
padded_lu!(invmat, e::EffectiveMatrix, ::Val{:R}) = lu!(view(invmat, e.r+1:2e.r+e.n, e.r+1:2e.r+e.n))
padded_lu!(invmat, e::EffectiveMatrix, ::Val{:LR}) = lu!(invmat)

padded_ldiv(invmat::AbstractMatrix, mat::AbstractMatrix, e::EffectiveMatrix, side) =
    padded_ldiv(padded_lu!(invmat, e, side), mat, e, side)
padded_ldiv(fact::Factorization, mat, e::EffectiveMatrix, side) =
    unpad_effective(ldiv!(fact, pad_effective(mat, e, side)), e, side)

Base.eltype(::EffectiveMatrix{T}) where {T} = T

#######################################################################
# EffectiveHamiltonian
#######################################################################


text:
#----------------------------------------------------------------------------------------------#
#                                 Dimensions-based Conversions                                 #
#----------------------------------------------------------------------------------------------#

"""
`function AMT(x::Number)`\n
Generates the default `AMOUNTS` from `a`, based on its unit dimensions.  The eltype-undecorated
`Quantity` constructors are evoked, so that the resulting type precision is taken from the `x`
argument. This function is extensively used in operations that result in a unit change.
"""
function AMT(x::Number)
    X, D = float(real(x)), dimension(x)
    # --- GenerAmt default
    if      D == dimension(1);              _Amt(X)
    # --- WholeAmt
    elseif  D == dimension(u"K");           sysT(X)
    elseif  D == dimension(u"kPa");         sysP(X)
    elseif  D == dimension(u"m/s");         VELO(X)
    elseif  D == dimension(u"s");           TIME(X)
    elseif  D == dimension(u"m/s^2");       grav(X)
    elseif  D == dimension(u"m");           alti(X)
    # --- BasedAmt
    elseif  D == dimension(u"kg");          mAmt(X)
    elseif  D == dimension(u"kg/s");        mAmt(X)
    elseif  D == dimension(u"kg/kmol");     mAmt(X)
    elseif  D == dimension(u"kmol");        nAmt(X)
    elseif  D == dimension(u"kmol/s");      nAmt(X)
    elseif  D == dimension(u"kmol/kg");     nAmt(X)
    elseif  D == dimension(u"m^3");         vAmt(X)
    elseif  D == dimension(u"m^3/s");       vAmt(X)
    elseif  D == dimension(u"m^3/kg");      vAmt(X)
    elseif  D == dimension(u"m^3/kmol");    vAmt(X)
    elseif  D == dimension(u"kJ");          ΔeAmt(X)    # energy fallback
    elseif  D == dimension(u"kJ/s");        ΔeAmt(X)
    elseif  D == dimension(u"kJ/kg");       ΔeAmt(X)
    elseif  D == dimension(u"kJ/kmol");     ΔeAmt(X)
    elseif  D == dimension(u"kJ/K");        ΔsAmt(X)    # ntropy fallback
    elseif  D == dimension(u"kJ/K/s");      ΔsAmt(X)
    elseif  D == dimension(u"kJ/K/kg");     ΔsAmt(X)
    elseif  D == dimension(u"kJ/K/kmol");   ΔsAmt(X)
    # --- GenerAmt fallback
    else                                    _Amt(X)
    end
end

export AMT

#----------------------------------------------------------------------------------------------#
#                                     Same-Type Operations                                     #
#----------------------------------------------------------------------------------------------#

+(x::AMOUNTS) = x
-(x::AMOUNTS) = (typeof(x).name.wrapper)(-amt(x))

#----------------------------------------------------------------------------------------------#
#                               Same-Unit (Same-Base) Operations                               #
#----------------------------------------------------------------------------------------------#

# Diff-{type,parameters} converting/promoting sum,sub of same-base energies
+(x::ENERGYA{𝗽,𝘅,𝗯}, y::ENERGYA{𝘀,𝘆,𝗯}) where {𝗽,𝘀,𝘅,𝘆,𝗯} = begin
    +(promote(map(x -> ΔeAmt(amt(x)), (x, y))...)...)
end
-(x::ENERGYA{𝗽,𝘅,𝗯}, y::ENERGYA{𝘀,𝘆,𝗯}) where {𝗽,𝘀,𝘅,𝘆,𝗯} = begin
    -(promote(map(x -> ΔeAmt(amt(x)), (x, y))...)...)
end

# Diff-{type,parameters} converting/promoting sum,sub of same-base entropies
+(x::NTROPYA{𝗽,𝘅,𝗯}, y::NTROPYA{𝘀,𝘆,𝗯}) where {𝗽,𝘀,𝘅,𝘆,𝗯} = begin
    +(promote(map(x -> ΔsAmt(amt(x)), (x, y))...)...)
end
-(x::NTROPYA{𝗽,𝘅,𝗯}, y::NTROPYA{𝘀,𝘆,𝗯}) where {𝗽,𝘀,𝘅,𝘆,𝗯} = begin
    -(promote(map(x -> ΔsAmt(amt(x)), (x, y))...)...)
end

# Diff-{type,parameters} converting/promoting sum,sub of velocities
+(x::VELOCYP{𝗽,𝘅}, y::VELOCYP{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    +(promote(map(x -> VELO(amt(x)), (x, y))...)...)
end
-(x::VELOCYP{𝗽,𝘅}, y::VELOCYP{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    -(promote(map(x -> VELO(amt(x)), (x, y))...)...)
end

## # Diff-{type,parameters} converting/promoting sum,sub of GenerAmt's
## +(x::GenerAmt{𝗽,𝘅}, y::GenerAmt{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
##     +(promote(map(x -> _Amt(amt(x)), (x, y))...)...)
## end
## -(x::GenerAmt{𝗽,𝘅}, y::GenerAmt{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
##     -(promote(map(x -> _Amt(amt(x)), (x, y))...)...)
## end

# Diff-{type,parameters} converting/promoting sum,sub of AMOUNTS'
+(x::AMOUNTS{𝗽,𝘅}, y::AMOUNTS{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    AMT(amt(+(promote(map(x -> _Amt(amt(x)), (x, y))...)...)))
end
-(x::AMOUNTS{𝗽,𝘅}, y::AMOUNTS{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    AMT(amt(-(promote(map(x -> _Amt(amt(x)), (x, y))...)...)))
end

#----------------------------------------------------------------------------------------------#
#                              Known-type Products and Divisions                               #
#----------------------------------------------------------------------------------------------#

# MA-based * mass => SY-based; with Unitful promotion
*(x::BasedAmt{𝗽,𝘅,MA}, y::mAmt{𝘀,𝘆,SY}) where {𝗽,𝘀,𝘅,𝘆} = begin
    (typeof(x).name.wrapper)(*(amt(x), amt(y)))
end
*(y::mAmt{𝘀,𝘆,SY}, x::BasedAmt{𝗽,𝘅,MA}) where {𝗽,𝘀,𝘅,𝘆} = x * y     # as to fallback

# MO-based * mole => SY-based; with Unitful promotion
*(x::BasedAmt{𝗽,𝘅,MO}, y::nAmt{𝘀,𝘆,SY}) where {𝗽,𝘀,𝘅,𝘆} = begin
    (typeof(x).name.wrapper)(*(amt(x), amt(y)))
end
*(y::nAmt{𝘀,𝘆,SY}, x::BasedAmt{𝗽,𝘅,MO}) where {𝗽,𝘀,𝘅,𝘆} = x * y     # as to fallback

# DT-based * TIME => SY-based; with Unitful promotion
*(x::BasedAmt{𝗽,𝘅,DT}, y::TIME{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    (typeof(x).name.wrapper)(*(amt(x), amt(y)))
end
*(y::TIME{𝘀,𝘆}, x::BasedAmt{𝗽,𝘅,DT}) where {𝗽,𝘀,𝘅,𝘆} = x * y        # as to fallback

# SY-based / mass => MA-based; with Unitful promotion
# SY-based / mole => MO-based; with Unitful promotion
# SY-based / TIME => DT-based; with Unitful promotion
/(x::BasedAmt{𝗽,𝘅,SY}, y::Union{mAmt{𝘀,𝘆,SY},nAmt{𝘀,𝘆,SY},TIME{𝘀,𝘆}}) where {𝗽,𝘀,𝘅,𝘆} = begin
    (typeof(x).name.wrapper)(/(amt(x), amt(y)))
end

# MA-based * MO-based mass => MO-based; with Unitful promotion
*(x::BasedAmt{𝗽,𝘅,MA}, y::mAmt{𝘀,𝘆,MO}) where {𝗽,𝘀,𝘅,𝘆} = begin
    (typeof(x).name.wrapper)(*(amt(x), amt(y)))
end
*(y::mAmt{𝘀,𝘆,MO}, x::BasedAmt{𝗽,𝘅,MA}) where {𝗽,𝘀,𝘅,𝘆} = x * y     # as to fallback

# MO-based * MA-based mole => MA-based; with Unitful promotion
*(x::BasedAmt{𝗽,𝘅,MO}, y::nAmt{𝘀,𝘆,MA}) where {𝗽,𝘀,𝘅,𝘆} = begin
    (typeof(x).name.wrapper)(*(amt(x), amt(y)))
end
*(y::nAmt{𝘀,𝘆,MA}, x::BasedAmt{𝗽,𝘅,MO}) where {𝗽,𝘀,𝘅,𝘆} = x * y     # as to fallback

#----------------------------------------------------------------------------------------------#
#                          Generic (fallback) Products and Divisions                           #
#----------------------------------------------------------------------------------------------#

*(x::AMOUNTS{𝗽,𝘅}, y::AMOUNTS{𝗽,𝘅}) where {𝗽,𝘅} = AMT(*(amt(x), amt(y)))
/(x::AMOUNTS{𝗽,𝘅}, y::AMOUNTS{𝗽,𝘅}) where {𝗽,𝘅} = AMT(/(amt(x), amt(y)))

*(x::AMOUNTS{𝗽,𝘅}, y::AMOUNTS{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    *(promote(map(x -> _Amt(amt(x)), (x, y))...)...)
end
/(x::AMOUNTS{𝗽,𝘅}, y::AMOUNTS{𝘀,𝘆}) where {𝗽,𝘀,𝘅,𝘆} = begin
    /(promote(map(x -> _Amt(amt(x)), (x, y))...)...)
end

#----------------------------------------------------------------------------------------------#
#                                    Other Base operations                                     #
#----------------------------------------------------------------------------------------------#

import Base: inv

inv(x::AMOUNTS) = AMT(inv(amt(x)))

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

import Base: ^, sqrt, cbrt

^(x::AMOUNTS, y::Real) = AMT(^(amt(x), y))
sqrt(x::AMOUNTS) = AMT(sqrt(amt(x)))
cbrt(x::AMOUNTS) = AMT(cbrt(amt(x)))

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

import Base: log, log2, log10

log(x::AMOUNTS) = _Amt(log(amt(x).val))
log2(x::AMOUNTS) = _Amt(log2(amt(x).val))
log10(x::AMOUNTS) = _Amt(log10(amt(x).val))

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

import Base: real, float, abs, abs2, min, max

real(x::AMOUNTS) = x
float(x::AMOUNTS) = x
abs(x::𝗧) where 𝗧<:AMOUNTS = 𝗧(abs(amt(x)))
abs2(x::AMOUNTS) = x^2

min(x::𝗧...) where 𝗧<:AMOUNTS = 𝗧(min((amt(i) for i in x)...))
max(x::𝗧...) where 𝗧<:AMOUNTS = 𝗧(max((amt(i) for i in x)...))

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

import Base: widen, eps

widen(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧.name.wrapper{widen(𝗽)}
widen(x::AMOUNTS) = widen(typeof(x))(x)

eps(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = (𝗧.name.wrapper)(eps(𝗽))
eps(x::𝗧) where 𝗧<:AMOUNTS = (𝗧.name.wrapper)(eps(amt(x)))

"""
`precof(::Type{𝗧} | x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗽`\n
Returns the precision of the `AMOUNTS` subtype or instance as a `DataType`.
"""
precof(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗽
precof(x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗽

"""
`exacof(::Type{𝗧} | x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗽`\n
Returns the exactness of the `AMOUNTS` subtype or instance as a `DataType`.
"""
exacof(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽,𝘅} where {𝗽,𝘅} = 𝘅
exacof(x::𝗧) where 𝗧<:AMOUNTS{𝗽,𝘅} where {𝗽,𝘅} = 𝘅

"""
`baseof(::Type{𝗧} | x::𝗧) where 𝗧<:BasedAmt{𝗽,𝘅,𝗯} where {𝗽,𝘅,𝗯} = 𝗯`\n
Returns the thermodynamic base of the `AMOUNTS` subtype or instance as a `DataType`.
"""
baseof(::Type{𝗧}) where 𝗧<:BasedAmt{𝗽,𝘅,𝗯} where {𝗽,𝘅,𝗯} = 𝗯
baseof(x::𝗧) where 𝗧<:BasedAmt{𝗽,𝘅,𝗯} where {𝗽,𝘅,𝗯} = 𝗯

export precof, exacof, baseof

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

import Base: prevfloat, nextfloat, zero, one, typemin, typemax

prevfloat(x::𝗧) where 𝗧<:AMOUNTS = (𝗧.name.wrapper)(prevfloat(amt(x)))
nextfloat(x::𝗧) where 𝗧<:AMOUNTS = (𝗧.name.wrapper)(nextfloat(amt(x)))

zero(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(zero(𝗽))
zero(x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(zero(𝗽))

one(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(one(𝗽))
one(x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(one(𝗽))

typemin(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(typemin(𝗽))
typemin(x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(typemin(𝗽))

typemax(::Type{𝗧}) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(typemax(𝗽))
typemax(x::𝗧) where 𝗧<:AMOUNTS{𝗽} where 𝗽 = 𝗧(typemax(𝗽))

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

import Base: floor, ceil, trunc, round, sign, signbit

for FUN in (:floor, :ceil, :trunc, :round)
    @eval $FUN(x::𝗧) where 𝗧<:AMOUNTS = (𝗧.name.wrapper)(($FUN)(amt(x).val) * unit(amt(x)))
end

round(x::𝗧, r::RoundingMode; digits, sigdigits, base) where 𝗧<:AMOUNTS = begin
    (𝗧.name.wrapper)(round(amt(x).val, r,
                           digits=digits,
                           sigdigits=sigdigits,
                           base=base) * unit(amt(x)))
end

for FUN in (:sign, :signbit)
    @eval $FUN(x::𝗧) where 𝗧<:AMOUNTS = ($FUN)(amt(x))
end

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

import Base: isfinite, isnan, isinf

for FUN in (:isfinite, :isnan, :isinf)
    @eval $FUN(x::𝗧) where 𝗧<:AMOUNTS = ($FUN)(amt(x))
end

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

import Base: ==, >, <, isequal, isless, isapprox

==(x::𝗧, y::𝗧) where 𝗧<:AMOUNTS = begin
    # We don't care about the 3 least significant bits of the wider type
    RTOL = (1<<3) * Base.rtoldefault(amt(x), amt(y), 0)
    isapprox(amt(x), amt(y), rtol=RTOL)
end

for FUN in (:>, :<, :isequal, :isless)
    @eval ($FUN)(x::𝗧, y::𝗧) where 𝗧<:AMOUNTS = ($FUN)(amt(x),amt(y))
end

function isapprox(x::𝗧, y::𝗧; atol::Real=0,
                  rtol::Real = Base.rtoldefault(amt(x), amt(y), atol),
                  nans::Bool=false) where 𝗧<:AMOUNTS
    isapprox(amt(x), amt(y), atol=atol, rtol=rtol, nans=nans)
end


text:
######
##FEM Poisson Nonlinear System Tests
######
using FiniteElementDiffEq

prob = prob_poisson_birthdeathsystem

sol = solve(prob,FEMDiffEqPoisson())

#plot(sol,plot_analytic=false,zlim=(0,2))

#Returns true if computed solution is homogenous near 2
bool1 = maximum(abs.(sol[end] .- [2 1]))< 1e-8

### Harder system

prob = prob_poisson_birthdeathinteractingsystem

sol = solve(prob,FEMDiffEqPoisson())

#plot(sol,plot_analytic=false,zlim=(0,2),cbar=false)

bool2 = maximum(abs.(sol[end] .- [2 1]))< 1e-8

bool1 && bool2


text:
using GraffSDK
using GraffSDK.DataHelpers
using ProgressMeter
using UUIDs
using Caesar

# 1a. Create a Configuration
config = loadGraffConfig()
# config.sessionId = "HexDemoSample1_"*replace(string(uuid4())[1:6], "-" => "")
config.sessionId = "HexDemoSample1_3529fd"

println(getGraffConfig())

# 1b. Check the credentials and the service status
printStatus()
# 1c. Check the session queue length
@info "Session backlog (queue length) = $(getSessionBacklog())"

# Variable data
GraffSDK.ls()
x0 = getVariable(:x0)
x0.packed
x0.type

## Getting all data entries for the whole session.
sessionDataEntries = getSessionDataEntries()

## Getting individual elements
# By variable
dataEntries = getDataEntries(getVariable(:x1))
@info dataEntries
dataEntry = GraffSDK.getData(getVariable("x1"), dataEntries[1])
# This is raw bytestream now, but we can convert it if we know it's JSON.
if dataEntry.mimeType == "application/json"
    @info String(dataEntry.data)
end
dataEntry = GraffSDK.getData(getVariable("x1"), dataEntries[2])
# TODO: Show as PNG

## Making a new element
# Now lets make a new element
newElement = BigDataElementRequest(
    "NewElement"*replace(string(uuid4())[1:6], "-" => ""),
    "Mongo", "A simple description", rand(UInt8, 100), "application/octet-stream")
GraffSDK.setData(getVariable("x1"), newElement)
dataEntries = getDataEntries(getVariable(:x1))
# Now let's get it back...
dataEntry = GraffSDK.getData(getVariable("x1"), newElement.id)

# Kablam!!! :D


text:
@testset "reported issues" begin
    @testset "plot hang (tan → ∞) (#129)" begin
        p = lineplot([cos, sin, tan], -π/2, 2π)
        test_ref("references/issues/cos_sin_tan.txt", @show_col(p))
    end

    @testset "integer edges (#139)" begin
        p_int = histogram(fit(Histogram, rand(10), 0:3))
        p_float = histogram(fit(Histogram, rand(10), 0.0:1.0:3.0))
        @test @print_col(p_int) == @print_col(p_float)
    end
end


text:
#!/usr/bin/env julia
using PkgServer, Sockets

# This stub script only really used for testing, but feel free to adapt to your needs
host = "127.0.0.1"
port = 8001
PkgServer.start(;
    listen_addr=Sockets.InetAddr(host, port),
    storage_root=mktempdir(),
    dotflavors=[""],
)

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
module LowLevelParticleFilters

export KalmanFilter, UnscentedKalmanFilter, ParticleFilter, AuxiliaryParticleFilter, AdvancedParticleFilter, SigmaFilter, PFstate, index, state, covariance, num_particles, effective_particles, weights, expweights, particles, particletype, smooth, sample_measurement, simulate, loglik, log_likelihood_fun, forward_trajectory, mean_trajectory, mode_trajectory, update!, predict!, correct!, reset!, metropolis, shouldresample, TupleProduct

export densityplot, debugplot, commandplot, trajectorydensity, dimensiondensity

using StatsBase, Parameters, Lazy, Yeppp, Random, LinearAlgebra, Printf
import PDMats # To extend some methods on static arrays
using StaticArrays
using Distributions
using Plots

abstract type ResamplingStrategy end
struct ResampleSystematic <: ResamplingStrategy end

abstract type AbstractFilter end

include("PFtypes.jl")
include("kalman.jl")
include("filtering.jl")
include("resample.jl")
include("utils.jl")
include("smoothing.jl")
include("plotting.jl")

index(f::AbstractFilter)               = f.t[]

end # module


text:
include("init.jl")

@testset "Iteration execution_iteration_test" begin

    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")
    # register the publications

    pubid = h.helicsFederateRegisterGlobalPublication(vFed1, "pub1", h.HELICS_DATA_TYPE_DOUBLE, "")

    subid = h.helicsFederateRegisterSubscription(vFed1, "pub1", "")
    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_DELTA, 1.0)

    h.helicsFederateEnterInitializingMode(vFed1)
    h.helicsPublicationPublishDouble(pubid, 27.0)

    comp = h.helicsFederateEnterExecutingModeIterative(vFed1, h.HELICS_ITERATION_REQUEST_ITERATE_IF_NEEDED)
    @test comp == h.HELICS_ITERATION_RESULT_ITERATING
    val = h.helicsInputGetDouble(subid)
    @test val == 27.0

    comp = h.helicsFederateEnterExecutingModeIterative(vFed1, h.HELICS_ITERATION_REQUEST_ITERATE_IF_NEEDED)

    @test comp == h.HELICS_ITERATION_RESULT_NEXT_STEP

    val2 = h.helicsInputGetDouble(subid)

    @test val2 == val

    h.helicsFederateFinalize(vFed1)
    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)
end

@testset "Iteration async test" begin

    broker = createBroker(1)
    vFed1, fedinfo1 = createValueFederate(1, "fed0")
    vFed2, fedinfo2 = createValueFederate(1, "fed1")

    # register the publications
    pub1 = h.helicsFederateRegisterGlobalPublication(vFed1, "pub1", h.HELICS_DATA_TYPE_INT)

    sub1 = h.helicsFederateRegisterSubscription(vFed2, "pub1")
    pub2 = h.helicsFederateRegisterGlobalPublication(vFed2, "pub2", h.HELICS_DATA_TYPE_INT)

    sub2 = h.helicsFederateRegisterSubscription(vFed1, "pub2")
    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)
    h.helicsFederateSetTimeProperty(vFed2, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)
    # vFed1->setLoggingLevel(5)
    # vFed2->setLoggingLevel(5)

    h.helicsFederateEnterInitializingModeAsync(vFed1)
    h.helicsFederateEnterInitializingMode(vFed2)
    h.helicsFederateEnterInitializingModeComplete(vFed1)

    c1 = 0
    c2 = 0

    h.helicsPublicationPublishInteger(pub1, c1)
    h.helicsPublicationPublishInteger(pub2, c2)

    h.helicsFederateEnterExecutingModeAsync(vFed1)
    h.helicsFederateEnterExecutingMode(vFed2)
    h.helicsFederateEnterExecutingModeComplete(vFed1)

    while c1 <= 10
        h.helicsInputGetInteger(sub1), c1
        h.helicsInputGetInteger(sub2), c2
        c1 += 1
        c2 += 1

        if c1 <= 10
            h.helicsPublicationPublishInteger(pub1, c1)
            h.helicsPublicationPublishInteger(pub2, c2)
        end

        h.helicsFederateRequestTimeIterativeAsync(vFed1, 1.0, h.HELICS_ITERATION_REQUEST_ITERATE_IF_NEEDED)

        grantedTime, state = h.helicsFederateRequestTimeIterative(vFed2, 1.0, h.HELICS_ITERATION_REQUEST_ITERATE_IF_NEEDED)
        if c1 <= 10
            @test state == h.HELICS_ITERATION_RESULT_ITERATING
            @test grantedTime == 0.0
        else
            @test state == h.HELICS_ITERATION_RESULT_NEXT_STEP
            @test grantedTime == 1.0
        end

        grantedTime, state = h.helicsFederateRequestTimeIterativeComplete(vFed1)
        if c1 <= 10
            @test state == h.HELICS_ITERATION_RESULT_ITERATING
            @test grantedTime == 0.0
        else
            @test state == h.HELICS_ITERATION_RESULT_NEXT_STEP
            @test grantedTime == 1.0
        end
    end

end


text:
"""
Struct for ShiftedLogNormalMixture

		lnmModel=ShiftedLogNormalMixture(η::Array{num1},λ::Array{num2},num3) where {num1,num2,num3<: Number}

Where:\n
		η  =	Array of volatilities.
		λ  = 	Array of weights.
		α  = 	shift.
"""
mutable struct ShiftedLogNormalMixture{num <: Number, num2 <: Number, num3 <: Number, abstrUnderlying <: AbstractUnderlying, numtype <: Number} <: ItoProcess{numtype}
    η::Array{num, 1}
    λ::Array{num2, 1}
    α::num3
    underlying::abstrUnderlying
    function ShiftedLogNormalMixture(η::Array{num, 1}, λ::Array{num2, 1}, α::num3, underlying::abstrUnderlying) where {num <: Number, num2 <: Number, num3 <: Number, abstrUnderlying <: AbstractUnderlying}
        if minimum(η) <= 0
            error("Volatilities must be positive")
        elseif minimum(λ) <= 0
            error("weights must be positive")
        elseif sum(λ) > 1.0
            error("λs must be weights")
        elseif length(λ) != length(η) - 1
            error("Check vector lengths")
        else
            zero_typed = zero(num) + zero(num2) + zero(num3)
            return new{num, num2, num3, abstrUnderlying, typeof(zero_typed)}(η, λ, α, underlying)
        end
    end
end

export ShiftedLogNormalMixture;

function simulate!(X, mcProcess::ShiftedLogNormalMixture, rfCurve::AbstractZeroRateCurve, mcBaseData::AbstractMonteCarloConfiguration, T::Number)
    @assert T > 0
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nstep = mcBaseData.Nstep
    η_gbm = copy(mcProcess.η)
    λ_gmb = copy(mcProcess.λ)
    mu_gbm = r - d
    dt = T / Nstep
    A0 = S0 * (1 - mcProcess.α)
    simulate!(X, LogNormalMixture(η_gbm, λ_gmb, Underlying(A0, d)), rfCurve, mcBaseData, T)
    tt = collect(0:dt:T)
    X .= X .+ mcProcess.α .* S0 .* exp.(integral(mu_gbm, t_) for t_ in tt)'
    return nothing
end


text:
module LJM
using CEnum

export LJDeviceType, LJConnectionType
export dtANY, T7, TSERIES, ctANY, USB, ETHERNET
import Base.cconvert
include("liblabjackm_common.jl")

Base.cconvert(::Type{Ref{Int32}}, a::Array{LJMDataType}) = Base.cconvert(Ref{Int32}, Cint.(a))

include("liblabjackm_api.jl")
end # module


text:
#========== GPU etc ==========#

using .GPUArrays
# https://github.com/JuliaGPU/GPUArrays.jl/blob/master/src/host/broadcast.jl

using Base.Broadcast
import Base.Broadcast: BroadcastStyle, Broadcasted, ArrayStyle

TransmuteGPU{AT} = TransmutedDimsArray{T,N,P,Q,AT} where {T,N,P,Q}

BroadcastStyle(::Type{<:TransmuteGPU{AT}}) where {AT<:AbstractGPUArray} =
    BroadcastStyle(AT)

GPUArrays.backend(::Type{<:TransmuteGPU{AT}}) where {AT<:AbstractGPUArray} =
    GPUArrays.backend(AT)

@inline function Base.copyto!(dest::TransmuteGPU{AT}, bc::Broadcasted{Nothing}) where {AT<:AbstractGPUArray}
    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
    isempty(dest) && return dest
    bc′ = Broadcast.preprocess(dest, bc)

    function broadcast_kernel(ctx, dest, bc′, nelem)
        for i in 1:nelem
            I = @cartesianidx(dest, i)
            @inbounds dest[I] = bc′[I]
        end
        return
    end
    heuristic = launch_heuristic(backend(dest), broadcast_kernel, dest, bc′, 1)
    config = launch_configuration(backend(dest), heuristic, length(dest), typemax(Int))
    gpu_call(broadcast_kernel, dest, bc′, config.elements_per_thread;
             threads=config.threads, blocks=config.blocks)

    return dest
end

@inline function Base.copyto!(
        dest::TransmuteGPU{AT},
        bc::Broadcasted{<:Broadcast.AbstractArrayStyle{0}}
    )  where {AT<:AbstractGPUArray}
    copyto!(dest, convert(Broadcasted{Nothing}, bc))
end

# https://github.com/JuliaGPU/GPUArrays.jl/blob/master/src/host/abstractarray.jl#L49
# display
Base.print_array(io::IO, X::TransmuteGPU{AT} where {AT <: AbstractGPUArray}) =
    Base.print_array(io, GPUArrays.convert_to_cpu(X))

Base._show_nonempty(io::IO, X::TransmuteGPU{AT} where {AT <: AbstractGPUArray}, prefix::String) =
    Base._show_nonempty(io, GPUArrays.convert_to_cpu(X), prefix)
Base._show_empty(io::IO, X::TransmuteGPU{AT} where {AT <: AbstractGPUArray}) =
    Base._show_empty(io, GPUArrays.convert_to_cpu(X))
Base.show_vector(io::IO, v::TransmuteGPU{AT} where {AT <: AbstractGPUArray}, args...) =
    Base.show_vector(io, GPUArrays.convert_to_cpu(X), args...)

using Adapt  # is this still needed?
# https://github.com/JuliaGPU/Adapt.jl/blob/master/src/base.jl

function Adapt.adapt_structure(to, A::TransmutedDimsArray{T,N,P,Q,AT}) where {T,N,P,Q,AT}
    data = adapt(to, A.parent)
    TransmutedDimsArray{eltype(data),N,P,Q,typeof(data)}(data)
end


text:
using Test
using LinearAlgebra
using DFTK: compute_fft_size, compute_fft_size_precise
using DFTK

include("testcases.jl")

@testset "Test compute_fft_size on Silicon" begin
    @test compute_fft_size(silicon.lattice,  3, supersampling=2) == [15, 15, 15]
    @test compute_fft_size(silicon.lattice,  4, supersampling=2) == [15, 15, 15]
    @test compute_fft_size(silicon.lattice,  5, supersampling=2) == [18, 18, 18]
    @test compute_fft_size(silicon.lattice, 15, supersampling=2) == [27, 27, 27]
    @test compute_fft_size(silicon.lattice, 25, supersampling=2) == [36, 36, 36]
    @test compute_fft_size(silicon.lattice, 30, supersampling=2) == [40, 40, 40]

    # Test the model interface as well
    model = Model(silicon.lattice; n_electrons=silicon.n_electrons)
    @test compute_fft_size(model, 30) == [40, 40, 40]
    @test compute_fft_size(model, 30, supersampling=1.8) == [36, 36, 36]
end

@testset "Test compute_fft_size on skewed lattice" begin
    lattice = Diagonal([1, 1e-12, 1e-12])
    @test compute_fft_size(lattice, 15, supersampling=2) == [5, 1, 1]
    @test compute_fft_size(lattice, 300, supersampling=2) == [18, 1, 1]
end

@testset "Test compute_fft_size_precise" begin
    atoms = [ElementPsp(:Si, psp=load_psp(silicon.psp)) => silicon.positions]
    model = Model(silicon.lattice; n_electrons=silicon.n_electrons, atoms=atoms)

    function fft_precise(model, kcoords, Ecut; supersampling=2)
        fft_size_fast = compute_fft_size(model.lattice, Ecut, supersampling=supersampling)
        kpoints = DFTK.build_kpoints(model, fft_size_fast, kcoords, Ecut; variational=true)
        fft_size = compute_fft_size_precise(model.lattice, Ecut, kpoints; supersampling=supersampling)
        @test all(fft_size .<= fft_size_fast)
        fft_size
    end

    @test fft_precise(model, silicon.kcoords, 20, supersampling=2)   == [30, 30, 30]
    @test fft_precise(model, silicon.kcoords, 20, supersampling=1.6) == [24, 24, 24]
end


text:
using OrdinaryDiffEq
using Trixi

##########  EINSTELLUNGEN     #########
# BLZ: we adjust the CFL and the terminal time to test the different methods for stability
# maximal setting for N = 3, Nq = 16: Chandrashekar without dissipation, CFL = 2.37, t = 100

initial_condition = initial_condition_weak_blast_wave
N = 3
cells_per_dimension = (16, 16)

CFL = 2.37
tspan = (0.0, 100)

# Standard DGSEM Entropy STability
# polydeg = N
# surface_flux = flux_lax_friedrichs # flux_chandrashekar # flux_central #
# volume_integral = VolumeIntegralWeakForm()
# solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral )

# CHandrashekar DGSEM Entropy STability
surface_flux = flux_chandrashekar # FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive)) #
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
basis = LobattoLegendreBasis(N)
solver = DGSEM(basis, surface_flux, volume_integral)

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

coordinates_min = (-2.0, -2.0)
coordinates_max = ( 2.0,  2.0)
mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)
analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                                 extra_analysis_errors=(:conservation_error,),
                                                 extra_analysis_integrals=(entropy,)# energy_total,
                                                                           #energy_kinetic, energy_internal)

                                                )

callbacks = CallbackSet(summary_callback, stepsize_callback, analysis_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

max_stars_repo_path: src/WR2/MS2 Stability.jl

max_issues_repo_path: src/WR2/MS2 Stability.jl

max_forks_repo_path: src/WR2/MS2 Stability.jl