text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

# Utility function that should not be directly called (only to be called from attach(..)
function attachAndReverseParents(newParent::Object3D, obj::Object3D)::NOTHING
   @assert(!(newParent ≡ obj))

   # Save elements of obj
   r_rel  = obj.r_rel
   R_rel  = obj.R_rel
   parent = obj.parent

   # Reverse obj, so that newParent is the new parent
   obj.parent = newParent
   push!(newParent.children, obj)

   # Reverse all parents of obj
   obj1::Object3D = obj
   obj2::Object3D = parent
   while !(parent === obj1)
      parent_r_rel = parent.r_rel
      parent_R_rel = parent.R_rel

      if r_rel ≡ ModiaMath.ZeroVector3D
         parent.r_rel = ModiaMath.ZeroVector3D
         parent.r_abs = obj1.r_abs
      else
         parent.r_rel = -R_rel*r_rel
      end

      if R_rel ≡ ModiaMath.NullRotation
         parent.R_rel = ModiaMath.NullRotation
         parent.R_abs = obj1.R_abs
      else
         parent.R_rel = R_rel'
      end

      obj2 = parent
      parent = obj2.parent
      r_rel  = parent_r_rel
      R_rel  = parent_R_rel
      obj2.parent = obj1
      removeChild!(obj2, obj1)
      push!(obj1.children, obj2)
      obj1 = obj2
   end

   # If scene != nothing, visit all objs and set scene
#=
   if typeof(scene) != NOTHING
      stack = scene.stack
      empty!(stack)
      append!(stack, obj.children)
      while length(stack) > 0
         obj1 = pop!(stack)
         obj1.scene = scene
         append!(stack, obj1.children)
      end
   end
=#
   return nothing
end

"""
    (obj1, obj2, cutJoint) = attach(frame_a, frame_b)

If frame_a and frame_b do not have the same root, they are attached to each other
and cutJoint = false is returned.

If they have the same root, the tree is not modified and cutJoint=true is returned.
"""
function attach(obj1::Object3D, obj2::Object3D)
   root1 = rootObject3D(obj1)
   root2 = rootObject3D(obj2)
   #println("attach: obj1 = ", ModiaMath.fullName(obj1), ", root = ", ModiaMath.fullName(root1))
   #println("attach: obj2 = ", ModiaMath.fullName(obj2), ", root = ", ModiaMath.fullName(root2))

   if root1 ≡ root2
      # Compute absolute positions
      updatePosition!(root1)
      return (obj1,obj2,true)
   end

   if isWorld(root2)
      attachAndReverseParents(obj2, obj1)
      return (obj2,obj1,false)
   else
      attachAndReverseParents(obj1, obj2)
      return (obj1,obj2,false)
   end
end

function connect(obj1::Object3D, obj2::Object3D;
                 r::AbstractVector = ModiaMath.ZeroVector3D,
                 R::Union{ModiaMath.RotationMatrix,NOTHING} = nothing,
                 q::Union{ModiaMath.Quaternion,NOTHING} = nothing,
                 fixed::Bool = true)::Nothing
   if typeof(R) != NOTHING && typeof(q) != NOTHING
      error("Modia3D.connect(...): either R or q must be nothing but both have a value.")
   end
   if typeof(R) != NOTHING
      ModiaMath.assertRotationMatrix(R)
   elseif typeof(q) != NOTHING
      ModiaMath.assertQuaternion(q)
   end

   (parentObject3D, obj, cutJoint) = attach(obj1, obj2)
  # println("... connect, fixed = ", fixed, ", obj1=",ModiaMath.fullName(obj1), ", obj2 = ", ModiaMath.fullName(obj2), ", obj = ", ModiaMath.fullName(obj))

   if cutJoint
      error("Error from Modia3D.Composition.connect(", obj1.name, ",", obj2.name, "):\n",
            "Not yet supported to rigidly connect two objs that have the same root.")
   end

   r_rel = SVector{3,Float64}(r)
   R_rel = typeof(R) != NOTHING ? R                   :
           typeof(q) != NOTHING ? ModiaMath.from_q(q) : ModiaMath.NullRotation

   obj.r_rel = obj===obj2 ? r_rel : -R_rel*r_rel
   obj.R_rel = obj===obj2 ? R_rel :  R_rel'

   #r_abs = parent.r_abs + r_rel
   #R_abs = R_rel*parent.R_abs

   if fixed
      obj.joint = fixedJoint
   else
      obj1.hasChildJoint = true
      q_start = typeof(R) != NOTHING ? ModiaMath.from_R(R) :
                typeof(q) != NOTHING ? q                   : ModiaMath.NullQuaternion
      q_start = obj===obj2 ? q_start : ModiaMath.inverseRotation(q_start)

      obj.joint = FreeMotion(obj, r_start = obj.r_rel, q_start = q_start)
   end
   return nothing
end

updatePosition!(assembly::Modia3D.AbstractAssembly) = updatePosition!(assembly._internal.referenceObject3D)
function updatePosition!(frame::Object3D)::NOTHING
   stack = Object3D[]
   # Push initial children on stack
   append!(stack, frame.children)
   while length(stack) > 0
      frame = pop!(stack)

      computeKinematics!(frame.joint, frame, ModiaMath.KinematicAnalysis, 0.0)

#=
      parent = frame.parent

      if frame.r_rel ≡ ModiaMath.ZeroVector3D
         frame.r_abs = parent.r_abs
      else
         frame.r_abs = parent.r_abs + parent.R_abs'*frame.r_rel
      end

      if frame.R_rel ≡ ModiaMath.NullRotation
         frame.R_abs = parent.R_abs
      else
         frame.R_abs = frame.R_rel*parent.R_abs
      end
=#
      if typeof(frame.visualizationFrame) != NOTHING
         frame.visualizationFrame.r_abs = frame.r_abs
         frame.visualizationFrame.R_abs = frame.R_abs
      end

      append!(stack, frame.children)
   end
   return nothing
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 5/3
equations = CompressibleEulerEquations3D(gamma)

initial_condition = initial_condition_blob

surface_flux = flux_hllc
volume_flux  = flux_ranocha
solver = DGSEM(3, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-20, -20, -20)
coordinates_max = ( 20,  20,  20)

refinement_patches = (
  (type="box", coordinates_min=(-20, -10, -10), coordinates_max=(-10, 10, 10)),
  (type="box", coordinates_min=(-20,  -5,  -5), coordinates_max=(-10,  5,  5)),
  (type="box", coordinates_min=(-17,  -2,  -2), coordinates_max=(-13,  2,  2)),
  (type="box", coordinates_min=(-17,  -2,  -2), coordinates_max=(-13,  2,  2)),
)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                refinement_patches=refinement_patches,
                n_cells_max=100_000,)

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

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

tspan = (0.0, 2.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 200
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

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

amr_indicator = IndicatorLöhner(semi,
                                variable=Trixi.density)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=1,
                                      med_level =0, med_threshold=0.1, # med_level = current level
                                      max_level =6, max_threshold=0.3)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=3,
                           adapt_initial_condition=false,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.7)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        amr_callback, stepsize_callback)

limiter! = PositivityPreservingLimiterZhangShu(thresholds=(1.0e-4, 1.0e-4),
                                               variables=(Trixi.density, pressure))
stage_limiter! = limiter!
step_limiter!  = limiter!

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

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, step_limiter!, 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

ext: jl


text:
function GetGlobalPosition(R::Route, Dis::Float64)
    GIndex, GDis = CurrGeo_Index_Dis(R, Dis)
    return Length2XY(R.Geos[GIndex], GDis)
end

function GetGlobalPosAngle(R::Route, Dis::Float64)
    GIndex, GDis = CurrGeo_Index_Dis(R, Dis)
    return Length2XYθ(R.Geos[GIndex], GDis)
end

function CurrGeo_Index_Dis(R::Route, distance::Float64)
    if distance > R.Length
        error("Input distance $distance is longer than route length $(R.Length)")
    end
    dis = distance
    index = 1
    for i in 1:length(R.Geos)
        if dis > round(R.Geos[i].Length, 1.0e-8) + 1.0e-8 # plus precision lose.
            dis -= round(R.Geos[i].Length, 1.0e-8) + 1.0e-8
            index += 1
        else
            break
        end
    end
    return (index, dis)
end

function Dist2Ip(R::Route, GDist::Float64, GIndex::UInt16) # total distance from start point to intersect point
    Dist::Float64 = 0.0
    if GIndex > 1
        for i in 1 : GIndex-1
            Dist += R.Geos[i].Length
        end
    end
    Dist += GDist
    return Dist
end

function Routes_Update!(Rs::Vector{Route}, RIndex::UInt16, newGeo::Geometry)
    push!(Rs[RIndex].Geos, newGeo)
    empty!(Rs.intersect_Infos)
    RoutesInit!(Rs)
end

function RoutesIntersect(R1::Route, R2::Route)
    IntersectInfos = []
    for i in 1:length(R1.Geos)  # what if there are more than one intersect points? what if the car has already passed this point?
        for j in 1:length(R2.Geos)
            Result = GeosIntersect(R1.Geos[i], R2.Geos[j])
            if Result.intersect
                push!(IntersectInfos, (Result..., GIndex1=UInt16(i), GIndex2=UInt16(j))) # if intersect, return the intersection point, distance in Geos and Geo indexes.
            end
        end
    end
    if !isempty(IntersectInfos)
        return IntersectInfos
    else
        return [(intersect=false,)] # if no intersection, just return (intersect=false,) #message="xxx")
    end
end

function RoutesInit!(Rs::Vector{Route})
    for i = 2:length(Rs)
        IntersectInfos = RoutesIntersect(Rs[1], Rs[i])
        PtDistVec1 = Vector{Tuple{String, Point2D{Float64}, Float64}}()
        PtDistVeck = Vector{Tuple{String, Point2D{Float64}, Float64}}()
        for Re in IntersectInfos
            if Re.intersect
                Dist2Ip_1 = Dist2Ip(Rs[1], Re.Dist1, Re.GIndex1)
                Dist2Ip_k = Dist2Ip(Rs[i], Re.Dist2, Re.GIndex2)
                push!(PtDistVec1, (Re.type, Re.Ip, Dist2Ip_1))
                push!(PtDistVeck, (Re.type, Re.Ip, Dist2Ip_k))
            end
        end
        if !isempty(PtDistVec1)
            push!(Rs[1].intersect_Infos, UInt16(i) => PtDistVec1) # used to calculate ΔTTC; ego car should not use this infomation
        end
        if !isempty(PtDistVeck)
            push!(Rs[i].intersect_Infos, UInt16(1) => PtDistVeck)
        end
    end
    return nothing
end

function RouteAVrefset!(R::Route, Sstep::Float64,
    Aset::NamedTuple{(:min, :max, :comfort),Tuple{Float64, Float64, Float64}})
    # first write the funciton of Vref related to s, then sample it with step = Sstep.
    Vvec = Vector{Float64}()
    Avec = Vector{Float64}()
    FVvec = Vector{Function}()
    FAvec = Vector{Function}()
    Geo_num = length(R.Geos)
    if Geo_num < 1
        error("""
              Route has empty Geos.
              Route: $R
              """)
    end
    if Geo_num == 1
        vec_length = 1+UInt16(floor(R.Geos[1].Length/Sstep))
        Vvec = [R.Geos[1].Vstd for i in 1:vec_length]
        Avec = [0.0 for i in 1:vec_length]
    else # Geo_num > 1
        for Gi in 1:Geo_num # one Geo, one Expression
            if R.Geos[Gi].GType == "line" # Vref changes only in straight line; otherwise keeps constant.
                if Gi == 1 # if no preceding Geo exists, calculate Vref with the Vstd of current and next Geo.
                    if  R.Geos[Gi].Vstd == R.Geos[Gi+1].Vstd # a == 0
                        push!(FVvec, s->R.Geos[Gi].Vstd)
                        push!(FAvec, s->0.0)
                    else # a != 0
                        rel_v, rel_a = Exp_approach(R.Geos[Gi].Vstd, R.Geos[Gi+1].Vstd, R.Geos[Gi].Length, Aset)
                        push!(FVvec, rel_v)
                        push!(FAvec, rel_a)
                    end
                else # Gi > 1;
                    # if preceding Geo exists, and has the same Vstd, then calculate Vref with the Vstd of current and next Geo.
                    if Gi != length(R.Geos) && R.Geos[Gi-1].Vstd == R.Geos[Gi].Vstd
                        if  R.Geos[Gi].Vstd == R.Geos[Gi+1].Vstd # a == 0
                            push!(FVvec, s->R.Geos[Gi].Vstd)
                            push!(FAvec, s->0.0)
                        else # a != 0
                            rel_v, rel_a = Exp_approach(R.Geos[Gi].Vstd, R.Geos[Gi+1].Vstd, R.Geos[Gi].Length, Aset)
                            push!(FVvec, rel_v)
                            push!(FAvec, rel_a)
                        end
                    else # Gi == length(R.Geos) || R.Geos[Gi-1].Vstd != R.Geos[Gi].Vstd)
                        # if preceding Geo exists, but has a different Vstd (if pre_Geo != line, Vend == Vstd),
                        # then calculate Vref with the Vend of preceding Geo and Vstd of current Geo // and Vstd of next Geo.
                        if R.Geos[Gi-1].Vstd != R.Geos[Gi].Vstd
                            rel_v, rel_a = Exp_leave(R.Geos[Gi-1].Vstd, R.Geos[Gi].Vstd, R.Geos[Gi].Length, Aset)
                            push!(FVvec, rel_v)
                            push!(FAvec, rel_a)
                        else # R.Geos[Gi-1].Vstd == R.Geos[Gi].Vstd && Gi == length(R.Geos)
                            push!(FVvec, s->R.Geos[Gi].Vstd)
                            push!(FAvec, s->0.0)
                        end
                    end
                end
            else # GType != "line"
                push!(FVvec, s->R.Geos[Gi].Vstd)
                push!(FAvec, s->0.0)
            end
        end
        if length(FVvec) == length(FAvec) == Geo_num
            for i in 0:UInt16(floor(R.Length/Sstep))
                Gi, Δs = CurrGeo_Index_Dis(R, i*Sstep)
                Vstdi = FVvec[Gi](Δs)
                Astdi = FAvec[Gi](Δs)
                push!(Vvec, Vstdi)
                push!(Avec, Astdi)
            end
        else
            error("""
            Expression number is not the same with Geometry number.
            Exp_num: $(length(FVvec))
            Geo_num: $Geo_num
            Related Route: $R
            """)
        end
    end
    R.Vref = Vvec
    R.Aref = Avec
end

function Exp_approach(Vnow::Float64, Vnext::Float64, Length,
    Aset::NamedTuple{(:min, :max, :comfort),Tuple{Float64, Float64, Float64}}) # entering

    a_req = (Vnext^2 - Vnow^2)/(2*Length) # a_req has sign
    if a_req > Aset.max || a_req < Aset.min
        error("""
              This Route is too short to accelerate/brake the car.
              Required a : $a_req
              Length available: $Length.
              """)
    end
    # Vehicle should accelerate at least with a comfort a.
    a = a_req
    if abs(a_req) < Aset.comfort
        a = a_req > 0 ? Aset.comfort : -Aset.comfort
    end
    s_acc = (Vnext^2 - Vnow^2)/(2*a)
    Δs = Length - s_acc
    rel_v(s::Float64) = begin
                            if s <= Δs
                                return Vnow
                            else
                                 return sqrt(Vnow^2 + 2*a*(s-Δs))
                            end
                       end

    rel_a(s::Float64) = begin
                            if s <= Δs
                                return 0.0
                            else
                                 return a
                            end
                       end

    return rel_v, rel_a
end

function Exp_leave(Vpre::Float64, Vnow::Float64, Length,
    Aset::NamedTuple{(:min, :max, :comfort),Tuple{Float64, Float64, Float64}})  # leaving
    a_req = (Vnow^2 - Vpre^2)/(2*Length)
    if a_req > Aset.max || a_req < Aset.min
        error("""
              This Route is too short to accelerate/brake the car.
              Vpre: $Vpre
              Vcurrent: $Vnow
              Required a : $a_req
              Length available: $Length.
              """)
    end
    # Vehicle should accelerate itself as fast as it can.
    a = a_req > 0 ? Aset.max : Aset.min
    s_acc = (Vnow^2 - Vpre^2)/(2*a)
    rel_v(s::Float64) = begin
                            if s < s_acc
                                return sqrt(Vpre^2 + 2*a*s)
                            else
                                 return Vnow
                            end
                       end

    rel_a(s::Float64) = begin
                           if s < s_acc
                               return a
                           else
                                return 0.0
                           end
                      end
    return rel_v, rel_a
end

function RoutesGenerator(Vgeos::Matrix{Float64}, Rlength::Float64, gauge::Float64, θ::Float64, GRatio_Other::Float64, GRatio_Ego::Float64; style::Symbol)

    if style == :Crossroad
        if size(Vgeos) != (4, 3)
            error("The size of Start Velocity Matrix is not correct! Input size is $(size(Vgeos)), while required size is 4x3.")
        end
    elseif style == :TJunction
        if size(Vgeos) != (3, 3)
            error("The size of Start Velocity Matrix is not correct! Input size is $(size(Vgeos)), while required size is 3x3.")
        end
    else
        error("The input style $style is not supported. Please choose either :Crossroad or :TJunciton.")
    end

    Linepts_shared = [0.0 0.0; gauge/sin(θ) gauge/sin(θ)+GRatio_Other*Rlength]
    Linepts_sharedR = Rmat(π/2-θ)*Linepts_shared
    Linepts_sharedR[1,:] .-= gauge/(2*sin(θ))
    Geo_shared = Geometry(Point2D{Float64}(Linepts_sharedR[1,2], Linepts_sharedR[2,2]), Point2D{Float64}(Linepts_sharedR[1,1], Linepts_sharedR[2,1]), Vgeos[2, 1])

    r3 = 1.5*gauge/(1-cos(θ))
    c3_endpt = (x=gauge*(1/sin(θ)+1/(2*tan(θ))), y=-gauge/2)
    Geo32 = Geometry(Point4D{Float64}(Linepts_sharedR[1,1], Linepts_sharedR[2,1], 1/r3, -θ), Point4D{Float64}(c3_endpt.x, c3_endpt.y, 1/r3, 0.0), Vgeos[3, 2])
    if Geo_shared.Length + Geo32.Length >= Rlength
        error("GRatio_Other $GRatio_Other is too big. Please Reduce it and retry.")
    end
    Geo33 = Geometry(Point2D{Float64}(c3_endpt.x, c3_endpt.y), Point2D{Float64}(c3_endpt.x+Rlength-Geo_shared.Length-Geo32.Length, c3_endpt.y), Vgeos[3, 3])
    R3 = Route([Geo_shared, Geo32, Geo33])

    r2 = 0.5*gauge/(1+cos(θ))
    Geo22 = Geometry(Point4D{Float64}(Linepts_sharedR[1,1], Linepts_sharedR[2,1], -1/r2, -θ), Point4D{Float64}(-c3_endpt.x, -c3_endpt.y, -1/r2, -π), Vgeos[2, 2])
    Geo23 = Geometry(Point2D{Float64}(-c3_endpt.x, -c3_endpt.y), Point2D{Float64}(-c3_endpt.x-Rlength+Geo_shared.Length+Geo22.Length, -c3_endpt.y), Vgeos[2, 3])
    R2 = Route([Geo_shared, Geo22, Geo23])

    if style == :Crossroad
        Line1pts = [0.0 0.0; -GRatio_Ego*Rlength (1.0-GRatio_Ego)*Rlength]
        Line1ptsR = Rmat(π/2-θ)*Line1pts
        Line1ptsR[1,:] .+= gauge/(2*sin(θ))
        Geo1 = Geometry(Point2D{Float64}(Line1ptsR[1,1], Line1ptsR[2,1]), Point2D{Float64}(Line1ptsR[1,2], Line1ptsR[2,2]), Vgeos[1, 1])
        R1 = Route([Geo1])

        line4endpt = Rmat(π/2-θ)*[0.0; gauge/sin(θ)-(1.0-GRatio_Other)*Rlength] - [gauge/(2*sin(θ)); 0.0]
        Geo42 = Geometry(Point2D{Float64}(Linepts_sharedR[1,1], Linepts_sharedR[2,1]), Point2D{Float64}(line4endpt[1], line4endpt[2]), Vgeos[4, 2])
        R4 = Route([Geo_shared, Geo42])

        Routes = [R1, R2, R3, R4]
        return Routes
    else # style == :TJunction
        Geo1 = Geometry(Point2D{Float64}(GRatio_Ego*Rlength, gauge/2), Point2D{Float64}((GRatio_Ego-1.0)*Rlength, gauge/2), Vgeos[1, 1])
        R1 = Route([Geo1])
        Routes = [R1, R2, R3]
        return Routes
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: s3outposts
using AWS.Compat
using AWS.UUIDs

"""
    CreateEndpoint()

S3 on Outposts access points simplify managing data access at scale for shared datasets in Amazon S3 on Outposts. S3 on Outposts uses endpoints to connect to Outposts buckets so that you can perform actions within your virtual private cloud (VPC).  This action creates an endpoint and associates it with the specified Outpost.   Related actions include:    DeleteEndpoint     ListEndpoints

# Required Parameters
- `OutpostId`: The ID of the AWS Outpost.
- `SecurityGroupId`: The ID of the security group to use with the endpoint.
- `SubnetId`: The ID of the subnet in the selected VPC.

"""
create_endpoint(OutpostId, SecurityGroupId, SubnetId; aws_config::AWSConfig=global_aws_config()) = s3outposts("POST", "/S3Outposts/CreateEndpoint", Dict{String, Any}("OutpostId"=>OutpostId, "SecurityGroupId"=>SecurityGroupId, "SubnetId"=>SubnetId); aws_config=aws_config)
create_endpoint(OutpostId, SecurityGroupId, SubnetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = s3outposts("POST", "/S3Outposts/CreateEndpoint", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("OutpostId"=>OutpostId, "SecurityGroupId"=>SecurityGroupId, "SubnetId"=>SubnetId), args)); aws_config=aws_config)

"""
    DeleteEndpoint()

S3 on Outposts access points simplify managing data access at scale for shared datasets in Amazon S3 on Outposts. S3 on Outposts uses endpoints to connect to Outposts buckets so that you can perform actions within your virtual private cloud (VPC).  This action deletes an endpoint.   Related actions include:    CreateEndpoint     ListEndpoints

# Required Parameters
- `endpointId`: The ID of the end point.
- `outpostId`: The ID of the AWS Outpost.

"""
delete_endpoint(endpointId, outpostId; aws_config::AWSConfig=global_aws_config()) = s3outposts("DELETE", "/S3Outposts/DeleteEndpoint", Dict{String, Any}("endpointId"=>endpointId, "outpostId"=>outpostId); aws_config=aws_config)
delete_endpoint(endpointId, outpostId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = s3outposts("DELETE", "/S3Outposts/DeleteEndpoint", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endpointId"=>endpointId, "outpostId"=>outpostId), args)); aws_config=aws_config)

"""
    ListEndpoints()

S3 on Outposts access points simplify managing data access at scale for shared datasets in Amazon S3 on Outposts. S3 on Outposts uses endpoints to connect to Outposts buckets so that you can perform actions within your virtual private cloud (VPC).  This action lists endpoints associated with the Outpost.   Related actions include:    CreateEndpoint     DeleteEndpoint

# Optional Parameters
- `maxResults`: The max number of endpoints that can be returned on the request.
- `nextToken`: The next endpoint requested in the list.
"""
list_endpoints(; aws_config::AWSConfig=global_aws_config()) = s3outposts("GET", "/S3Outposts/ListEndpoints"; aws_config=aws_config)
list_endpoints(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = s3outposts("GET", "/S3Outposts/ListEndpoints", args; aws_config=aws_config)


text:
import PowerGraphics
using PowerSystems
const PG = PowerGraphics
include("../../src/get_test_data.jl")

# 3.0

path = mkdir(joinpath(pwd(), "plots-01"));
PG.stack_plot(re_results; save=path, display=false, title="Example GR Plot");

# ## To make an interactive PlotlyJS plot, reset the backend
#Plots.plotlyjs()
path = mkdir(joinpath(pwd(), "plots-02"));
PG.stack_plot(re_results; save=path, display=false, title="Example PlotlyJS Plot");

# 3.1

path = mkdir(joinpath(pwd(), "plots-10"));
PG.stack_plot(re_results; save=path, title="Example Stack Plot");
path = mkdir(joinpath(pwd(), "plots-11"));
PG.stack_plot(re_results; save=path, title="Example saved Stack Plot");
path = mkdir(joinpath(pwd(), "plots-12"));
PG.stack_plot(
    re_results;
    reserves=true,
    save=path,
    format="png",
    title="Example Stack Plot with Reserves",
);

path = mkdir(joinpath(pwd(), "plots-13"));
colors = [:pink :green :blue :magenta :black]
PG.stack_plot(
    re_results;
    seriescolor=colors,
    save=path,
    format="png",
    title="Example Stack Plot with Other Colors",
);

path = mkdir(joinpath(pwd(), "plots-14"));
PG.stack_plot(re_results; stair=true, save=path, title="Stair Plot");
path = mkdir(joinpath(pwd(), "plots-15"));
title = "Example of a Title"
PG.stack_plot(re_results; save=path, title=title);
# 3.2
path = mkdir(joinpath(pwd(), "plots-20"));
PG.bar_plot(re_results; save=path, title="Example Bar Plot");

path = mkdir(joinpath(pwd(), "plots-21"));
PG.bar_plot(re_results; save=path, title="Example saved Bar Plot");

path = mkdir(joinpath(pwd(), "plots-22"))
PG.bar_plot(re_results; reserves=true, save=path, title="Example Bar Plot with Reserves");

path = mkdir(joinpath(pwd(), "plots-23"))
PG.bar_plot(
    re_results;
    seriescolor=colors,
    save=path,
    title="Example Bar Plot with Other Colors",
);

title = "Example of a Title"

path = mkdir(joinpath(pwd(), "plots-24"));
PG.bar_plot(re_results; save=path, title=title);

# 3.3

path = mkdir(joinpath(pwd(), "plots-3"));
PG.fuel_plot(re_results, system; save=path, title="Example Fuel Plot");
PG.fuel_plot(re_results, system; save=path, load=true, title="Fuel Plot with Load");
PG.fuel_plot(
    re_results,
    system;
    save=path,
    curtailment=true,
    title="Fuel Plot with Curtailment",
);
PG.fuel_plot(
    op_results,
    system;
    reserves=true,
    save=path,
    title="Example Fuel Plot with Reserves",
);

colors = [:pink :green :blue :magenta :black]
PG.fuel_plot(
    re_results,
    system;
    seriescolor=colors,
    save=path,
    title="Example Fuel Plot with Other Colors",
);

title = "Example of a Title";
PG.fuel_plot(re_results, system; save=path, title=title);
PG.fuel_plot(re_results, system; save=path, stair=true);

# 3.4 FORECAST PLOTS
path = mkdir(joinpath(pwd(), "plots-4"));
PG.plot_reserves(op_results; save=path);
PG.plot_demand(re_results; save=path, title="Example Demand Plot");
PG.plot_demand(system; save=path, title="Example Demand Plot From System");
initial_time = Dates.DateTime(2024, 01, 01, 02, 0, 0);
horizon = 6;

PG.plot_demand(
    system;
    horizon=horizon,
    initial_time=initial_time,
    save=path,
    title="Example Demand Plot Subsection",
);

PG.plot_demand(system; aggregate=PSY.System, save=path, title="Example Demand Plot by Type");

PG.plot_demand(system; stair=true, title="Example Stair Demand Plot", save=path);
#
colors = [:orange :pink :blue :red :grey]
PG.plot_demand(
    system;
    seriescolor=colors,
    save=path,
    title="Example Demand Plot with Different Colors",
);

PG.plot_demand(system; title="Example Demand Plot with Title", save=path);

# 3.5

path = mkdir(joinpath(pwd(), "plots-50"));

PG.plot_variable(op_results, "P__ThermalStandard"; save=path);
path = mkdir(joinpath(pwd(), "plots-53"));

p = PG.plot_variable(re_results, "P__ThermalStandard");
PG.plot_variable(p, re_results, "P__RenewableDispatch"; title="overlay", save=path);
path = mkdir(joinpath(pwd(), "plots-54"));
PG.plot_dataframe(
    re_results.variable_values[:P__ThermalStandard],
    re_results.timestamp;
    save=path,
);
path = mkdir(joinpath(pwd(), "plots-55"));

p2 =
    PG.plot_dataframe(re_results.variable_values[:P__ThermalStandard], re_results.timestamp);
PG.plot_dataframe(
    p2,
    re_results.variable_values[:P__RenewableDispatch],
    re_results.timestamp;
    title="overlay",
    save=path,
    format="png",
);

variables = [Symbol("P__ThermalStandard")]
path = mkdir(joinpath(pwd(), "plots-51"));

PG.stack_plot(op_results, variables; save=path, title="Plot with Fewer Variables");
#
path = mkdir(joinpath(pwd(), "plots-52"));

selected_variables = Dict(Symbol("P__ThermalStandard") => [:Brighton, :Solitude]);
results_subset = PG.sort_data(op_results; Variables=selected_variables);

PG.stack_plot(results_subset; save=path, title="Selected Variables Plot");

# 3.6
results_one = PSI.run_economic_dispatch(system; optimizer=solver, use_parameters=true);
results_two = PSI.run_economic_dispatch(system; optimizer=solver, use_parameters=true);

path = mkdir(joinpath(pwd(), "plots-6"));
PG.stack_plot([results_one, results_two]; save=path, title="Comparison");
#
PG.fuel_plot([results_one, results_two], system; save=path, title="Comparison");
#
variables = [Symbol("P__ThermalStandard")]
PG.stack_plot(
    [results_one, results_two],
    variables;
    save=path,
    title="Comparison with fewer variables",
);


text:
include("utils.jl")
mathy_makedocs()
deploydocs(;
    repo="github.com/tkf/Mathy.jl",
)


text:
module Notations # Musial

export A, B, C, D, E, F, G

struct Notation
    name::String
end

const A = Notation("A")
const B = Notation("B")
const C = Notation("C")
const D = Notation("D")
const E = Notation("E")
const F = Notation("F")
const G = Notation("G")

end # module Musical.Notations


text:
"""
For educational purposes to demonstrate how to create composite types in Julia
and combine them in different ways. Show how to use composition, define interfaces and more.
"""
module ToyRockets

include("interfaces.jl")
include("physics.jl")

include("tank.jl")
include("tanks.jl")
include("flexi-tank.jl")

include("engines.jl")
include("custom-engine.jl")
include("cluster.jl")
include("cluster-index.jl")

include("payloads.jl")
include("empty-rocket.jl")
include("staged-rocket.jl")
include("iterate-rocket.jl")
include("operations.jl")

include("simulate.jl")

end


text:
"""
    ARDTransform(v::AbstractVector)

Transformation that multiplies the input elementwise by `v`.

# Examples

```jldoctest
julia> v = rand(10); t = ARDTransform(v); X = rand(10, 100);

julia> map(t, ColVecs(X)) == ColVecs(v .* X)
true
```
"""
struct ARDTransform{Tv<:AbstractVector{<:Real}} <: Transform
    v::Tv
end

"""
    ARDTransform(s::Real, dims::Integer)

Create an [`ARDTransform`](@ref) with vector `fill(s, dims)`.
"""
ARDTransform(s::Real, dims::Integer) = ARDTransform(fill(s, dims))

@functor ARDTransform

function set!(t::ARDTransform{<:AbstractVector{T}}, ρ::AbstractVector{T}) where {T<:Real}
    @assert length(ρ) == dim(t) "Trying to set a vector of size $(length(ρ)) to ARDTransform of dimension $(dim(t))"
    return t.v .= ρ
end

dim(t::ARDTransform) = length(t.v)

(t::ARDTransform)(x::Real) = first(t.v) * x
(t::ARDTransform)(x) = t.v .* x

_map(t::ARDTransform, x::AbstractVector{<:Real}) = t.v' .* x
_map(t::ARDTransform, x::ColVecs) = ColVecs(t.v .* x.X)
_map(t::ARDTransform, x::RowVecs) = RowVecs(t.v' .* x.X)

Base.isequal(t::ARDTransform, t2::ARDTransform) = isequal(t.v, t2.v)

Base.show(io::IO, t::ARDTransform) = print(io, "ARD Transform (dims: ", dim(t), ")")


text:
# Support Data for color maps
include("maps_data.jl")

# color scale generation
# ----------------------

"""
    distinguishable_colors(n, [seed]; [transform, lchoices, cchoices, hchoices])

Generate n maximally distinguishable colors.

This uses a greedy brute-force approach to choose n colors that are maximally
distinguishable. Given seed color(s), and a set of possible hue, chroma, and
lightness values (in LCHab space), it repeatedly chooses the next color as the
one that maximizes the minimum pairwise distance to any of the colors already
in the palette.

Args:

- `n`: Number of colors to generate.
- `seed`: Initial color(s) included in the palette.  Default is `Vector{RGB{N0f8}}(0)`.

Keyword arguments:

- `transform`: Transform applied to colors before measuring distance. Default is the identity; other choices include `deuteranopic` to simulate color-blindness.
- `lchoices`: Possible lightness values (default `linspace(0,100,15)`)
- `cchoices`: Possible chroma values (default `linspace(0,100,15)`)
- `hchoices`: Possible hue values (default `linspace(0,340,20)`)

Returns:
  A `Vector` of colors of length `n`, of the type specified in `seed`.
"""
function distinguishable_colors{T<:Color}(n::Integer,
                            seed::AbstractVector{T};
                            transform::Function = identity,
                            lchoices::AbstractVector = linspace(0, 100, 15),
                            cchoices::AbstractVector = linspace(0, 100, 15),
                            hchoices::AbstractVector = linspace(0, 340, 20))
    if n <= length(seed)
        return seed[1:n]
    end

    # Candidate colors
    N = length(lchoices)*length(cchoices)*length(hchoices)
    candidate = Vector{Lab{Float64}}(N)
    j = 0
    for h in hchoices, c in cchoices, l in lchoices
        rgb = convert(RGB, LCHab(l, c, h))
        candidate[j+=1] = convert(LCHab, rgb)
    end

    # Transformed colors
    candidate_t = Vector{Lab{Float64}}(N)
    for i = 1:N
        candidate_t[i] = transform(candidate[i])
    end

    # Start with the seed colors
    colors = Vector{T}(n)
    copy!(colors, seed)

    # Minimum distances of the current color to each previously selected color.
    ds = fill(Inf, N)
    for i = 1:length(seed)
        ts = convert(Lab{Float64}, transform(seed[i]))::Lab{Float64}
        for k = 1:N
            ds[k] = min(ds[k], colordiff(ts, candidate_t[k]))
        end
    end

    for i in length(seed)+1:n
        j = indmax(ds)
        colors[i] = candidate[j]
        tc = candidate_t[j]
        for k = 1:N
            d = colordiff(tc, candidate_t[k])
            ds[k] = min(ds[k], d)
        end
    end

    colors
end

distinguishable_colors(n::Integer, seed::Color; kwargs...) = distinguishable_colors(n, [seed]; kwargs...)
distinguishable_colors(n::Integer; kwargs...) = distinguishable_colors(n, Vector{RGB{N0f8}}(0); kwargs...)

@deprecate distinguishable_colors(n::Integer,
                                transform::Function,
                                seed::Color,
                                ls::Vector{Float64},
                                cs::Vector{Float64},
                                hs::Vector{Float64})    distinguishable_colors(n, [seed], transform = transform, lchoices = ls, cchoices = cs, hchoices = hs)

# Sequential palette
# ----------------------
# Sequential_palette implements the color palette creation technique by
# Wijffelaars, M., et al. (2008)
# http://magnaview.nl/documents/MagnaView-M_Wijffelaars-Generating_color_palettes_using_intuitive_parameters.pdf
#
# Colormaps are formed using Bezier curves in LCHuv colorspace
# with some constant hue. In addition, start and end points can be given
# that are then blended to the original hue smoothly.
#
# The arguments are:
# N        - number of colors
# h        - the main hue [0,360]
# c        - the overall lightness contrast [0,1]
# s        - saturation [0,1]
# b        - brightness [0,1]
# w        - cold/warm parameter, i.e. the strength of the starting color [0,1]
# d        - depth of the ending color [0,1]
# wcolor   - starting color (warmness)
# dcolor   - ending color (depth)
# logscale - true/false for toggling logspacing
#
function sequential_palette(h,
                            N::Int=100;
                            c=0.88,
                            s=0.6,
                            b=0.75,
                            w=0.15,
                            d=0.0,
                            wcolor=RGB(1,1,0),
                            dcolor=RGB(0,0,1),
                            logscale=false)

    function MixHue(a,h0,h1)
        M=mod(180.0+h1-h0, 360)-180.0
        mod(h0+a*M, 360)
    end
    function mix_linearly{C<:Color}(a::C, b::C, s)
        base_color_type(C)((1-s)*comp1(a)+s*comp1(b), (1-s)*comp2(a)+s*comp2(b), (1-s)*comp3(a)+s*comp3(b))
    end

    pstart=convert(LCHuv, wcolor)
    p1=MSC(h)
    #p0=LCHuv(0,0,h) #original end point
    pend=convert(LCHuv, dcolor)

    #multi-hue start point
    p2l=100*(1.-w)+w*pstart.l
    p2h=MixHue(w,h,pstart.h)
    p2c=min(MSC(p2h,p2l), w*s*pstart.c)
    p2=LCHuv(p2l,p2c,p2h)

    #multi-hue ending point
    p0l=20.0*d
    p0h=MixHue(d,h,pend.h)
    p0c=min(MSC(p0h,p0l), d*s*pend.c)
    p0=LCHuv(p0l,p0c,p0h)

    q0=mix_linearly(p0,p1,s)
    q2=mix_linearly(p2,p1,s)
    q1=mix_linearly(q0,q2,0.5)

    pal = RGB{Float64}[]

    if logscale
        absc = logspace(-2.,0.,N)
    else
        absc = linspace(0.,1.,N)
    end

    for t in absc
        u=1.0-t

        #Change grid to favor light colors and to be uniform along the curve
        u = (125.0-125.0*0.2^((1.0-c)*b+u*c))
        u = invBezier(u,p0.l,p2.l,q0.l,q1.l,q2.l)

        #Get color components from Bezier curves
        ll = Bezier(u, p0.l, p2.l, q0.l, q1.l, q2.l)
        cc = Bezier(u, p0.c, p2.c, q0.c, q1.c, q2.c)
        hh = Bezier(u, p0.h, p2.h, q0.h, q1.h, q2.h)

        push!(pal, convert(RGB, LCHuv(ll,cc,hh)))
    end

    pal
end

# Diverging palettes
# ----------------------
# Create diverging palettes by combining 2 sequential palettes
#
# The arguments are:
# N        - number of colors
# h1       - the main hue of the left side [0,360]
# h2       - the main hue of the right side [0,360]
# c        - the overall lightness contrast [0,1]
# s        - saturation [0,1]
# b        - brightness [0,1]
# w        - cold/warm parameter, i.e. the strength of the starting color [0,1]
# d1       - depth of the ending color in the left side [0,1]
# d2       - depth of the ending color in the right side [0,1]
# wcolor   - starting color (warmness)
# dcolor1  - ending color of the left side (depth)
# dcolor2  - ending color of the right side (depth)
# logscale - true/false for toggling logspacing
#
function diverging_palette(h1,
                           h2,
                           N::Int=100;
                           mid=0.5,
                           c=0.88,
                           s=0.6,
                           b=0.75,
                           w=0.15,
                           d1=0.0,
                           d2=0.0,
                           wcolor=RGB(1,1,0),
                           dcolor1=RGB(1,0,0),
                           dcolor2=RGB(0,0,1),
                           logscale=false)

    if isodd(N)
        n=N-1
    else
        n=N
    end
    N1 = max(ceil(Int, mid*n), 1)
    N2 = Int(max(n-N1, 1))

    pal1 = sequential_palette(h1, N1+1, w=w, d=d1, c=c, s=s, b=b, wcolor=wcolor, dcolor=dcolor1, logscale=logscale)
    pal1 = flipdim(pal1, 1)

    pal2 = sequential_palette(h2, N2+1, w=w, d=d2, c=c, s=s, b=b, wcolor=wcolor, dcolor=dcolor2, logscale=logscale)

    if isodd(N)
        midcol = weighted_color_mean(0.5, pal1[end], pal2[1])
        return [pal1[1:end-1]; midcol; pal2[2:end]]
    else
        return [pal1[1:end-1]; pal2[2:end]]
    end
end

# Colormap
# ----------------------
# Main function to handle different predefined colormaps
#
"""
    colormap(cname, [N; mid, logscale, kvs...])

Returns a predefined sequential or diverging colormap computed using
the algorithm by Wijffelaars, M., et al. (2008).  Sequential colormaps
`cname` choices are `Blues`, `Greens`, `Grays`, `Oranges`, `Purples`,
and `Reds`.  Diverging colormap choices are `RdBu`.  Optionally, you
can specify the number of colors `N` (default 100). Keyword arguments
include the position of the middle point `mid` (default 0.5) and the
possibility to switch to log scaling with `logscale` (default false).
"""
function colormap(cname::AbstractString, N::Int=100; mid=0.5, logscale=false, kvs...)

    cname = lowercase(cname)
    if haskey(colormaps_sequential, cname)
        p = copy(colormaps_sequential[cname][1:8])

        const allowedkeys = [:h, :w, :d, :c, :s, :b, :wcolor, :dcolor]
        for (k,v) in kvs
            k in allowedkeys || throw(ArgumentError("Unknown keyword argument $k"))
            ind = findfirst(allowedkeys, k)
            if ind > 0
                p[ind] = v
            end
        end

        return sequential_palette(p[1], N, w=p[2], d=p[3], c=p[4], s=p[5], b=p[6], wcolor=p[7], dcolor=p[8], logscale=logscale)

    elseif haskey(colormaps_diverging, cname)
        p = copy(colormaps_diverging[cname][1:11])

        for (k,v) in kvs
            k in allowedkeys || throw(ArgumentError("Unknown keyword argument $k"))
            ind = findfirst(allowedkeys, k)
            if ind > 0
                p[ind] = v
            end
        end

        return diverging_palette(p[1], p[2], N, w=p[3], d1=p[4], d2=p[5], c=p[6], s=p[7], b=p[8], wcolor=p[9], dcolor1=p[10], dcolor2=p[11], mid=mid, logscale=logscale)

    else
        throw(ArgumentError(string("Unknown colormap: ", cname)))
    end
end


text:
# wrap_silo.jl: Generate wrapping code. See Sundials.jl (thanks!)
# This file is not an active part of the package. This is the code
# that uses the Clang.jl package to wrap Silo using the headers.

cd(dirname(@__FILE__))

using Clang.wrap_c

# silo_names = ["silo.h",]

#Clang include path for system clang: /usr/include/clang/3.4/include

JULIAHOME = "/home/jgoldfar/Documents/work/projects/julia"

clang_includes = map(x->joinpath(JULIAHOME, x),[
                       "usr/include/clang/",
                       "usr/include",
                       "usr/include/llvm",
                       "usr/include/llvm-c",
                       ]
                     )

path = "/home/jgoldfar/Documents/work/projects/Silo.jl/deps/usr/include"
headerlist = ["silo.h",]
check_use_header(path) = true

clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS", "-v"]
context = wrap_c.init(output_file = "libsilo.jl",
                      headers = [joinpath(path, x) for x in headerlist],
                      header_library=x->"libsilo",
                      clang_args = clang_extraargs,
                      clang_includes = clang_includes,
                      common_file = "libsilo_h.jl",
#                       rewriter = x->begin
#                         println(x)
#                         replace(x, "DB", "")
#                       end
                      )

run(context)


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

using Markdown
using InteractiveUtils

# ╔═╡ 844417d4-247b-11eb-19ed-6f83ab3e8a3d
using Plots

# ╔═╡ 26cbdbe6-3001-11eb-1302-5d0412fe34df
md"A paragraph containing a link to [Julia](http://www.julialang.org)."

# ╔═╡ 8b5121d8-2479-11eb-1160-8506756e005d
begin
	import Pkg
	Pkg.add("Plots")
end

# ╔═╡ c395c30a-2486-11eb-16ec-9968033a9d5a

# ╔═╡ Cell order:
# ╟─26cbdbe6-3001-11eb-1302-5d0412fe34df
# ╠═8b5121d8-2479-11eb-1160-8506756e005d
# ╠═844417d4-247b-11eb-19ed-6f83ab3e8a3d
# ╠═c395c30a-2486-11eb-16ec-9968033a9d5a


text:
export matching_poly

"""
`matching_poly(G)` returns the matching polynomial of the
`SimpleGraph` `G`.
"""
function matching_poly(G::SimpleGraph, cache_flag::Bool=true)
    if cache_flag && cache_check(G,:matching_poly)
        return cache_recall(G,:matching_poly)
    end
    if NE(G)==0
        p = Poly([0,1])^NV(G)
        if cache_flag
            SimpleGraphs.cache_save(G,:matching_poly,p)
        end
        return p
    end

    if is_connected(G)
        e = first(G.E)  # get any edge
        v,w = e

        G1 = deepcopy(G)
        SimpleGraphs.delete!(G1,v,w)
        p1 = matching_poly(G1,false)

        SimpleGraphs.delete!(G1,v)
        SimpleGraphs.delete!(G1,w)
        p2 = matching_poly(G1,false)

        p = p1 - p2
        if cache_flag
            cache_save(G,:matching_poly,p)
        end
        return p
    end

    comps = parts(components(G))
    p = Poly([1])
    for S in comps
        H = induce(G,S)
        pH = matching_poly(H,false)
        p *= pH
    end
    if cache_flag
        cache_save(G,:matching_poly,p)
    end
    return p
end


text:
using CUDAdrv, CUDAnative
using CUDAnativelib
using .CURANDkernel # notice the dot prefix

function fillRandom(out)
  i = (blockIdx().x-1) * blockDim().x + threadIdx().x
  if i <= length(out)
    # Initialise state
    state = new(curandState_t)
    curand_init(0, i, 0, state) # one should use a decent seed here ;)
    @inbounds out[i] = curand_uniform(eltype(out), state) # @inbounds is optional
    delete(state)
  end
  return nothing
end

# code_warntype(STDOUT, fillRandom, (CuDeviceArray{Float32,1},))
# code_warntype(STDOUT, fillRandom, (CuDeviceArray{Float64,1},))

# CUDAnative.code_llvm(STDOUT, fillRandom, (CuDeviceArray{Float32,1},))
# CUDAnative.code_llvm(STDOUT, fillRandom, (CuDeviceArray{Float64,1},))

dev = CuDevice(0)
ctx = CuContext(dev)
N = 100
out = CuArray{Float32}(N)
@cuda (N,1) fillRandom(out)
c = Array(out) # force sync


text:
include("crossover.jl")
include("mutation.jl")
include("BinaryGenotype/includes.jl")
include("IntegerGenotype/includes.jl")

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
function grad(Q::MPIStateArray, _1)
    Q_x1 = view(Q, (:,:,_1))

    v_∇Q_x1 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
    s_∇Q_x1 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)

    event = launch_volume_gradient!(grid, v_∇Q, Q_x1)
    wait(event)

    event = launch_interface_gradient!(grid, s_∇Q, Q_x1)
    wait(event)

    grad = v_∇Q_x1 + s_∇Q_x1

    return grad
end

function grad(Q::MPIStateArray, _1, _2, _3)
    Q_x1 = view(Q, (:,:,_1))
    Q_x2 = view(Q, (:,:,_2))
    Q_x3 = view(Q, (:,:,_3))

    v_∇Q_x1 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
    v_∇Q_x2 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
    v_∇Q_x3 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)

    s_∇Q_x1 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
    s_∇Q_x2 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
    s_∇Q_x3 = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)

    ## strong form gradients here
    event_x1 = launch_volume_gradient!(grid, v_∇Q_x1, Q_x1)
    event_x2 = launch_volume_gradient!(grid, v_∇Q_x2, Q_x2)
    event_x3 = launch_volume_gradient!(grid, v_∇Q_x3, Q_x3)
    wait(event_x1)
    wait(event_x2)
    wait(event_x3)

    ## strong form gradients here
    event_x1 = launch_interface_gradient!(grid, s_∇Q_x1, Q_x1)
    event_x2 = launch_interface_gradient!(grid, s_∇Q_x2, Q_x2)
    event_x3 = launch_interface_gradient!(grid, s_∇Q_x3, Q_x3)
    wait(event_x1)
    wait(event_x2)
    wait(event_x3)

    ∇Q_x1 = v_∇Q_x1 + s_∇Q_x1
    ∇Q_x2 = v_∇Q_x2 + s_∇Q_x2
    ∇Q_x3 = v_∇Q_x3 + s_∇Q_x3

    return (; ∇Q_x1, ∇Q_x2, ∇Q_x3)
end

function div(Q::MPIStateArray, _1, _2, _3)
    ∇Q_x1, ∇Q_x2, ∇Q_x3 = grad(Q, _1, _2, _3)

    div = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 1)
    div[:,:,1] = ∇Q_x1[:,:,1] + ∇Q_x2[:,:,2] + ∇Q_x3[:,:,3]

    return div
end

function curl(Q::MPIStateArray, _1, _2, _3)
    ∇Q_x1, ∇Q_x2, ∇Q_x3 = grad(Q, _1, _2, _3)

    curl = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)

    curl[:,:,1] = ∇Q_x3[:,:,2] - ∇Q_x2[:,:,3]
    curl[:,:,2] = ∇Q_x1[:,:,3] - ∇Q_x3[:,:,1]
    curl[:,:,3] = ∇Q_x2[:,:,1] - ∇Q_x1[:,:,2]

    return curl
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using Distributions, OnlineStats

function queue1(until::Real, burnin::Real,
                ad::Distribution, sd::Distribution)
    now, nextArrival, nextDeparture = 0.0, rand(ad), Inf
    queue, waits = Float64[], Mean()
    while now < until
        if nextArrival < nextDeparture
            now = nextArrival
            if isempty(queue)
                nextDeparture = nextArrival + rand(sd)
            end
            push!(queue, nextArrival)
            nextArrival += rand(ad)
        else
            now = nextDeparture
            insystem = nextDeparture - popfirst!(queue)
            burnin < now < until && fit!(waits, insystem)
            nextDeparture += isempty(queue) ? Inf : rand(sd)
        end
    end
    value(waits)
end

mm1_exact(ad::Exponential, sd::Exponential) = 1/(1/mean(sd)-1/mean(ad))

function bootCI(data, stat::Function, CI::Float64, reps::Integer)
    boot = [stat(rand(data, length(data))) for i in 1:reps]
    low, high = quantile(boot, [(1-CI)/2, (1+CI)/2])
    (value=stat(data), low=low, high=high)
end

max_stars_repo_path: Chapter04/06. Analyzing a queuing system/mm1.jl


text:
"""
"""
abstract type AssemblyStrategy end

"""
"""
function row_map(a::AssemblyStrategy,row)
  @abstractmethod
end

"""
"""
function col_map(a::AssemblyStrategy,col)
  @abstractmethod
end

"""
"""
function row_mask(a::AssemblyStrategy,row)
  @abstractmethod
end

"""
"""
function col_mask(a::AssemblyStrategy,col)
  @abstractmethod
end

struct DefaultAssemblyStrategy <: AssemblyStrategy end

row_map(a::DefaultAssemblyStrategy,row) = row

col_map(a::DefaultAssemblyStrategy,col) = col

row_mask(a::DefaultAssemblyStrategy,row) = true

col_mask(a::DefaultAssemblyStrategy,col) = true

"""
"""
abstract type Assembler <: GridapType end

"""
"""
function get_test(a::Assembler)
  @abstractmethod
end

"""
"""
function get_trial(a::Assembler)
  @abstractmethod
end

"""
"""
function get_assembly_strategy(a::Assembler)
  @abstractmethod
end

"""
"""
function allocate_matrix(a::Assembler,matdata)
  @abstractmethod
end

"""
"""
function allocate_vector(a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function allocate_matrix_and_vector(a::Assembler,data)
  @abstractmethod
end

"""
"""
function assemble_matrix!(A,a::Assembler,matdata)
  @abstractmethod
end

"""
"""
function assemble_matrix_add!(mat,a::Assembler, matdata)
  @abstractmethod
end

"""
"""
function assemble_vector!(b,a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function assemble_vector_add!(b,a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function assemble_matrix_and_vector!(A,b,a::Assembler, data)
  @abstractmethod
end

function assemble_matrix_and_vector_add!(A,b,a::Assembler, data)
  @abstractmethod
end

"""
"""
function assemble_matrix(a::Assembler,matdata)
  A = allocate_matrix(a,matdata)
  assemble_matrix!(A,a,matdata)
  A
end

"""
"""
function assemble_vector(a::Assembler,vecdata)
  b = allocate_vector(a,vecdata)
  assemble_vector!(b,a,vecdata)
  b
end

"""
"""
function assemble_matrix_and_vector(a::Assembler,data)
  A, b = allocate_matrix_and_vector(a,data)
  assemble_matrix_and_vector!(A,b,a,data)
  (A, b)
end

"""
"""
function test_assembler(a::Assembler,matdata,vecdata,data)
  trial_fesp = get_trial(a)
  test_fesp = get_test(a)
  A = allocate_matrix(a,matdata)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  assemble_matrix!(A,a,matdata)
  assemble_matrix_add!(A,a,matdata)
  A = assemble_matrix(a,matdata)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  b = allocate_vector(a,vecdata)
  @test num_free_dofs(test_fesp) == length(b)
  assemble_vector!(b,a,vecdata)
  assemble_vector_add!(b,a,vecdata)
  b = assemble_vector(a,vecdata)
  @test num_free_dofs(test_fesp) == length(b)
  A, b = allocate_matrix_and_vector(a,data)
  assemble_matrix_and_vector!(A,b,a,data)
  assemble_matrix_and_vector_add!(A,b,a,data)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  @test num_free_dofs(test_fesp) == length(b)
  A, b = assemble_matrix_and_vector(a,data)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  @test num_free_dofs(test_fesp) == length(b)
end

# Some syntactic sugar for assembling from anonymous functions
# and objects from which one can collect cell matrices/vectors

function assemble_matrix(f::Function,a::Assembler)
  v = get_cell_shapefuns(get_test(a))
  u = get_cell_shapefuns_trial(get_trial(a))
  assemble_matrix(a,collect_cell_matrix(f(u,v)))
end

function assemble_vector(f::Function,a::Assembler)
  v = get_cell_shapefuns(get_test(a))
  assemble_vector(a,collect_cell_vector(f(v)))
end

function assemble_matrix_and_vector(f::Function,b::Function,a::Assembler)
  v = get_cell_shapefuns(get_test(a))
  u = get_cell_shapefuns_trial(get_trial(a))
  assemble_matrix_and_vector(a,collect_cell_matrix_and_vector(f(u,v),b(v)))
end

function assemble_matrix(f,a::Assembler)
  assemble_matrix(a,collect_cell_matrix(f))
end

function assemble_vector(f,a::Assembler)
  assemble_vector(a,collect_cell_vector(f))
end

function assemble_matrix_and_vector(f,b,a::Assembler)
  assemble_matrix_and_vector(a,collect_cell_matrix_and_vector(f,b))
end

function assemble_matrix(f,U::FESpace,V::FESpace)
  a = SparseMatrixAssembler(U,V)
  assemble_matrix(f,a)
end

function assemble_vector(f,V::FESpace)
  a = SparseMatrixAssembler(V,V)
  assemble_vector(f,a)
end

function assemble_matrix_and_vector(f,b,U::FESpace,V::FESpace)
  a = SparseMatrixAssembler(U,V)
  assemble_matrix_and_vector(f,b,a)
end

# Abstract interface for computing the data to be sent to the assembler

function collect_cell_matrix(mat_contributions)
  @abstractmethod
end

function collect_cell_vector(vec_contributions)
  @abstractmethod
end

function collect_cell_matrix_and_vector(mat_contributions,vec_contributions)
  @abstractmethod
end

function collect_cell_matrix_and_vector(mat_contributions,vec_contributions,uhd::FEFunction)
  @abstractmethod
end

# Implementation of this interface for DomainContribution

function collect_cell_matrix(a::DomainContribution)
  w = []
  r = []
  for trian in get_domains(a)
    cell_mat = get_contribution(a,trian)
    @assert eltype(cell_mat) <: AbstractMatrix
    push!(w,cell_mat)
    push!(r,get_cell_to_bgcell(trian))
  end
  (w,r,r)
end

function collect_cell_vector(a::DomainContribution)
  w = []
  r = []
  for trian in get_domains(a)
    cell_vec = get_contribution(a,trian)
    @assert eltype(cell_vec) <: AbstractVector
    push!(w,cell_vec)
    push!(r,get_cell_to_bgcell(trian))
  end
  (w,r)
end

function _collect_cell_matvec(a::DomainContribution)
  w = []
  r = []
  for trian in get_domains(a)
    cell_mat = get_contribution(a,trian)
    @assert eltype(cell_mat) <: Tuple
    push!(w,cell_mat)
    push!(r,get_cell_to_bgcell(trian))
  end
  (w,r,r)
end

function collect_cell_matrix_and_vector(
  biform::DomainContribution,liform::DomainContribution)

  matvec, mat, vec = _pair_contribution_when_possible(biform,liform)
  matvecdata = _collect_cell_matvec(matvec)
  matdata = collect_cell_matrix(mat)
  vecdata = collect_cell_vector(vec)
  (matvecdata, matdata, vecdata)
end

function collect_cell_matrix_and_vector(
  biform::DomainContribution,liform::DomainContribution,uhd::FEFunction)

  matvec, mat, vec = _pair_contribution_when_possible(biform,liform,uhd)

  matvecdata = _collect_cell_matvec(matvec)
  matdata = collect_cell_matrix(mat)
  vecdata = collect_cell_vector(vec)
  (matvecdata, matdata, vecdata)
end

function _pair_contribution_when_possible(biform,liform)
  matvec = DomainContribution()
  mat = DomainContribution()
  vec = DomainContribution()
  for (trian,t) in biform.dict
    if haskey(liform.dict,trian)
      matvec.dict[trian] = pair_arrays(t,liform.dict[trian])
    else
      mat.dict[trian] = t
    end
  end
  for (trian,t) in liform.dict
    if ! haskey(biform.dict,trian)
      vec.dict[trian] = t
    end
  end
  matvec, mat, vec
end

function _pair_contribution_when_possible(biform,liform,uhd)
  _matvec, _mat, _vec = _pair_contribution_when_possible(biform,liform)
  matvec = DomainContribution()
  mat = DomainContribution()
  for (trian,t) in _matvec.dict
    cellvals = get_cell_dof_values(uhd,get_cell_to_bgcell(trian))
    matvec.dict[trian] = attach_dirichlet(t,cellvals)
  end
  for (trian,t) in _mat.dict
    cellvals = get_cell_dof_values(uhd,get_cell_to_bgcell(trian))
    matvec.dict[trian] = attach_dirichlet(t,cellvals)
  end
  matvec, mat, _vec
end

# allow linear forms like `l(v) = 0

function collect_cell_vector(l::Number)
  @notimplementedif l != 0
  w = []
  r = []
  (w,r)
end

function collect_cell_matrix_and_vector(mat_contributions::DomainContribution,l::Number)
  @notimplementedif l != 0
  vec_contributions = DomainContribution()
  collect_cell_matrix_and_vector(mat_contributions,vec_contributions)
end

function collect_cell_matrix_and_vector(
  biform::DomainContribution,l::Number,uhd::FEFunction)
  liform = DomainContribution()
  collect_cell_matrix_and_vector(biform,liform,uhd)
end


text:
"""
A Box Truncated normal distribution with a naive implementation and representation in the state of the mean and covariance. Works well for very low dimensions (e.g. 2,3,4).
"""
const BasicBoxTruncatedMvNormal = TruncatedMvDistribution{MvNormal,BoxTruncationRegion,TruncatedMvDistributionSecondOrderState}

function BasicBoxTruncatedMvNormal( μₑ::Vector{Float64},
                                    Σₑ::PDMat,
                                    a::Vector{Float64},
                                    b::Vector{Float64})
    d = MvNormal(μₑ,Σₑ)
    r = BoxTruncationRegion(a,b)
    TruncatedMvDistribution{MvNormal,BoxTruncationRegion,TruncatedMvDistributionSecondOrderState}(d,r)
end

"""
A Box Truncated normal distribution with a recursive moment computation implementation.
"""
const RecursiveMomentsBoxTruncatedMvNormal = TruncatedMvDistribution{MvNormal,BoxTruncationRegion,BoxTruncatedMvNormalRecursiveMomentsState}

function RecursiveMomentsBoxTruncatedMvNormal(  μₑ::Vector{Float64},
                                                Σₑ::PDMat,
                                                a::Vector{Float64},
                                                b::Vector{Float64};
                                                max_moment_levels::Int = 2)
    d = MvNormal(μₑ,Σₑ)
    r = BoxTruncationRegion(a,b)
    s = BoxTruncatedMvNormalRecursiveMomentsState(d,r,max_moment_levels)
    TruncatedMvDistribution{MvNormal,BoxTruncationRegion,BoxTruncatedMvNormalRecursiveMomentsState}(d,r,s)
end

# const EllipsoidTruncatedMvNormal = TruncatedMvDistribution{MvNormal,EllipticalTruncationRegion}
# const PolytopeTruncatedMvNormal = TruncatedMvDistribution{MvNormal,PolytopeTruncationRegion}


text:
function ProperFractionFormat()
    return ProperFractionFormat(())
end

function ProperFractionFormat(arg0::NumberFormat)
    return ProperFractionFormat((NumberFormat,), arg0)
end

function ProperFractionFormat(arg0::NumberFormat, arg1::NumberFormat, arg2::NumberFormat)
    return ProperFractionFormat((NumberFormat, NumberFormat, NumberFormat), arg0, arg1, arg2)
end

function format(obj::ProperFractionFormat, arg0::Fraction, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (Fraction, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function get_whole_format(obj::ProperFractionFormat)
    return jcall(obj, "getWholeFormat", NumberFormat, ())
end

function parse(obj::ProperFractionFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", Fraction, (JString, ParsePosition), arg0, arg1)
end


text:
################################################################################
#                      TEST: DCM <=> Euler Angle and Axis
################################################################################

for k = 1:samples

    # Euler angle and axis => DCM
    # ==========================================================================

    # The test will use the conversion from DCM to Quaternions as the
    # ground-truth value, since it was extensively tested.
    v  = randn(3)
    v /= norm(v)
    θ  = -π + 2π*rand()

    ea = EulerAngleAxis(θ,v)
    q  = Quaternion( cos(θ/2), sin(θ/2)*v )
    D1 = quat_to_dcm(q)
    D2 = angleaxis_to_dcm(ea)

    @test norm(D1-D2) ≈ 0 atol = 1e-13

    # Test rotations along the Euler axis.
    r     = v*randn()
    r_rot = D2*r

    @test norm(r-r_rot) ≈ 0 atol = 1e-10

    # DCM => Euler angle and axis
    # ==========================================================================

    eac1 = dcm_to_angleaxis(D1)
    eac2 = dcm_to_angleaxis(D2)

    # Get the angle in the interval [0,π] with the corresponding vector.
    θ = mod(θ,2π)

    s = 1
    if θ > π
        s = -1
        θ = 2π - θ
    end

    @test    eac1.a ≈ θ      atol = 5e-7
    @test eac1.v[1] ≈ s*v[1] atol = 5e-7
    @test eac1.v[2] ≈ s*v[2] atol = 5e-7
    @test eac1.v[3] ≈ s*v[3] atol = 5e-7

    @test    eac2.a ≈ θ      atol = 5e-7
    @test eac2.v[1] ≈ s*v[1] atol = 5e-7
    @test eac2.v[2] ≈ s*v[2] atol = 5e-7
    @test eac2.v[3] ≈ s*v[3] atol = 5e-7
end

# Test special cases
# ==================

dcm_z = angle_to_dcm(π, 0, 0, :ZYX)
dcm_y = angle_to_dcm(0, π, 0, :ZYX)
dcm_x = angle_to_dcm(0, 0, π, :ZYX)

@test EulerAngleAxis(π, [0, 0, 1]) ≈ dcm_to_angleaxis(dcm_z)  atol = 1e-14
@test EulerAngleAxis(π, [0, 1, 0]) ≈ dcm_to_angleaxis(dcm_y)  atol = 1e-14
@test EulerAngleAxis(π, [1, 0, 0]) ≈ dcm_to_angleaxis(dcm_x)  atol = 1e-14
@test EulerAngleAxis(0, [0, 0, 0]) ≈ dcm_to_angleaxis(DCM(I)) atol = 1e-14


text:
nlevels = 2

nvar = 9

nfactors = [1, 2]

fassign = [1 1
    1 1
    1 1
    1 1
    1 2
    1 2
    1 2
    1 2
    1 2]

flags = [2, 2]

varlags = [2, 2, 2, 2, 2, 2, 2, 2, 2]

varcoefs = [0.0 1.0 1.0
    0.0 0.5 0.2
    0.0 0.7 0.4
    0.0 0.3 0.5
    0.0 0.5 1.0
    0.0 0.5 0.7
    0.0 0.4 0.5
    0.0 0.5 0.2
    0.0 0.5 0.2]

varlagcoefs = [0.5 0.25
    0.5 0.25
    0.5 0.25
    0.5 0.25
    0.5 0.25
    0.5 0.25
    0.5 0.25
    0.5 0.25
    0.5 0.25]

fcoefs = Any[]
fmat = [0.85 -0.3][:, :]
push!(fcoefs, fmat)
fmat = [0.5 0.05
    0.2 -0.1]
push!(fcoefs, fmat)

fvars = Any[]
fmat = [1.0]
push!(fvars, fmat)
fmat = [1.0, 1.0]
push!(fvars, fmat)

varvars = 0.5 * ones(nvar);

hdfm = HDFM(nlevels = nlevels,
    nvar = nvar,
    nfactors = nfactors,
    fassign = fassign,
    flags = flags,
    varlags = varlags,
    varcoefs = varcoefs,
    varlagcoefs = varlagcoefs,
    fcoefs = fcoefs,
    fvars = fvars,
    varvars = varvars)

ssmodel = convertHDFMtoSS(hdfm)

num_obs = 100
data_y, data_z, data_β = simulateSSModel(num_obs, ssmodel::SSModel)

#flags = [3, 3]

#varlags = [3, 3, 3, 3, 3, 3, 3, 3, 3]

hdfmpriors = HDFMStruct(nlevels = nlevels,
    nfactors = nfactors,
    factorassign = fassign,
    factorlags = flags,
    errorlags = varlags,
    ndraws = 1000,
    burnin = 50)

results = PCA2LevelEstimator(data_y, hdfmpriors)

stds = Any[]
quant33 = Any[]
quant66 = Any[]
medians = Any[]

j = 1
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, j, :]))
    push!(quant33, quantile(results.F[i, j, :], 0.05))
    push!(quant66, quantile(results.F[i, j, :], 0.95))
    push!(medians, median(results.F[i, j, :]))
end

plot(data_β[:, 1+j])
plot!(results.means.F[:, j])
#plot!(medians)
#plot!(quant33)
#plot!(quant66)


text:
using ArgParse
using Pkg

Pkg.instantiate()

# Parse command line arguments.
parser = ArgParseSettings()
@add_arg_table! parser begin
    "--data"
        help =
            "Data set: eq-small, eq, matern52, eq-mixture, noisy-mixture, " *
            "weakly-periodic, sawtooth, or mixture. " *
            "Append \"-noisy\" to a data set to make it noisy."
        arg_type = String
        required = true
    "--model"
        help =
            "Model: conv[c]np, corconvcnp, a[c]np, or [c]np. " *
            "Append \"-global-{mean,sum}\" to introduce a global latent variable. " *
            "Append \"-amortised-{mean,sum}\" to use amortised observation noise. " *
            "Append \"-het\" to use heterogeneous observation noise."
        arg_type = String
        required = true
    "--num-samples"
        help =
            "Number of samples to estimate the training loss. Defaults to 20 for " *
            "\"loglik\" and 5 for \"elbo\"."
        arg_type = Int
    "--batch-size"
        help = "Batch size."
        arg_type = Int
        default = 16
    "--loss"
        help = "Loss: loglik, loglik-iw, or elbo."
        arg_type = String
        required = true
    "--starting-epoch"
        help = "Set to a number greater than one to continue training."
        arg_type = Int
        default = 1
    "--epochs"
        help = "Number of epochs to training for."
        arg_type = Int
        default = 20
    "--evaluate"
        help = "Evaluate model."
        action = :store_true
    "--evaluate-iw"
        help = "Force to use importance weighting for the evaluation objective."
        action = :store_true
    "--evaluate-no-iw"
        help = "Force to NOT use importance weighting for the evaluation objective."
        action = :store_true
    "--evaluate-num-samples"
        help = "Number of samples to estimate the evaluation loss."
        arg_type = Int
        default = 4096
    "--evaluate-only-within"
        help = "Evaluate with only the task of interpolation within training range."
        action = :store_true
    "--models-dir"
        help = "Directory to store models in."
        arg_type = String
        default = "models"
    "--bson"
        help = "Directly specify the file to save the model to and load it from."
        arg_type = String
end
args = parse_args(parser)

using BSON
using NeuralProcesses
using NeuralProcesses.Experiment
using Distributions
using Flux
using Stheno
using Tracker

# Determine the noise level.
if endswith(args["data"], "-noisy")
    trimmed_data_name = args["data"][1:end - length("-noisy")]
    noise = 0.05^2  # Use a total error of 0.05.
else
    trimmed_data_name = args["data"]
    noise = 1e-8  # Use very little noise, but still some for regularisation.
end

# Set up experiment.
if trimmed_data_name == "eq-small"
    process = GP(stretch(EQ(), 1 / 0.25), GPC())
    receptive_field = 1f0
    points_per_unit = 32f0
    num_context = DiscreteUniform(0, 50)
    num_target = DiscreteUniform(50, 50)
    num_channels = 16
    dim_embedding = 32
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "eq"
    process = GP(stretch(EQ(), 1 / 0.25), GPC())
    receptive_field = 2f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 50)
    num_target = DiscreteUniform(50, 50)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "matern52"
    process = GP(stretch(Matern52(), 1 / 0.25), GPC())
    receptive_field = 2f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 50)
    num_target = DiscreteUniform(50, 50)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "noisy-mixture"
    process = GP(stretch(EQ(), 1 / 0.25) + EQ() + 1e-3 * Stheno.Noise(), GPC())
    receptive_field = 4f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 50)
    num_target = DiscreteUniform(50, 50)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "eq-mixture"
    process = GP(stretch(EQ(), 1 / 0.25) + EQ(), GPC())
    receptive_field = 4f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 50)
    num_target = DiscreteUniform(50, 50)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "weakly-periodic"
    process = GP(stretch(EQ(), 1 / 0.5) * stretch(Stheno.PerEQ(1), 1 / 0.25), GPC())
    receptive_field = 4f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 50)
    num_target = DiscreteUniform(50, 50)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "sawtooth"
    process = Sawtooth()
    receptive_field = 16f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 100)
    num_target = DiscreteUniform(100, 100)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
elseif trimmed_data_name == "mixture"
    process = Mixture(
        GP(stretch(EQ(), 1 / 0.25), GPC()),
        GP(stretch(Matern52(), 1 / 0.25), GPC()),
        GP(stretch(EQ(), 1 / 0.25) + EQ() + 1e-3 * Stheno.Noise(), GPC()),
        GP(stretch(EQ(), 1 / 0.5) * stretch(Stheno.PerEQ(1), 1 / 0.25), GPC()),
        Sawtooth()
    )
    receptive_field = 16f0
    points_per_unit = 64f0
    num_context = DiscreteUniform(0, 100)
    num_target = DiscreteUniform(100, 100)
    num_channels = 64
    dim_embedding = 128
    num_context_eval = DiscreteUniform(0, 10)
else
    error("Unknown data \"" * args["data"] * "\".")
end

# Set the loss.
if args["model"] in ["convcnp", "corconvcnp", "acnp", "cnp"]
    # Determine training loss.
    if args["loss"] == "loglik"
        # Use a single sample: there is nothing random.
        loss(xs...) = NeuralProcesses.loglik(xs..., num_samples=1)
    elseif args["loss"] in ["elbo", "loglik-iw"]
        error("Losses \"elbo\" and \"loglik-iw\" not applicable conditional models.")
    else
        error("Unknown loss \"" * args["loss"] * "\".")
    end

    # Use the train loss for evaluation.
    eval_loss = loss
elseif args["model"] in [
    "convnp",
    "convnp-global-mean", "convnp-global-sum",
    "convnp-amortised-mean", "convnp-amortised-sum",
    "convnp-het",
    "anp",
    "anp-amortised-mean", "anp-amortised-sum",
    "anp-het",
    "np",
    "np-amortised-mean", "np-amortised-sum",
    "np-het"
]
    # Determine training loss.
    if args["loss"] == "loglik"
        if !isnothing(args["num-samples"])
            num_samples = args["num-samples"]
            args["loss"] *= "-$num_samples"  # Incorporate number of samples in the loss.
        else
            num_samples = 20
        end
        loss(xs...) = NeuralProcesses.loglik(
            xs...,
            num_samples=num_samples,
            importance_weighted=false,
            fixed_σ_epochs=3
        )
        eval_importance_weighted = false  # Encoder is not suited for IW.
    elseif args["loss"] == "loglik-iw"
        if !isnothing(args["num-samples"])
            num_samples = args["num-samples"]
            args["loss"] *= "-$num_samples"  # Incorporate number of samples in the loss.
        else
            num_samples = 20
        end
        loss(xs...) = NeuralProcesses.loglik(
            xs...,
            num_samples=num_samples,
            importance_weighted=true,
            fixed_σ_epochs=3
        )
        eval_importance_weighted = true  # Encoder is suited for IW!
    elseif args["loss"] == "elbo"
        if !isnothing(args["num-samples"])
            num_samples = args["num-samples"]
            args["loss"] *= "-$num_samples"  # Incorporate number of samples in the loss.
        else
            num_samples = 5
        end
        loss(xs...) = NeuralProcesses.elbo(
            xs...,
            num_samples=num_samples,
            fixed_σ_epochs=3
        )
        eval_importance_weighted = true  # Encoder is suited for IW!
    else
        error("Unknown loss \"" * args["loss"] * "\".")
    end

    # Check if `eval_importance_weighted` needs to be forced.
    if args["evaluate-iw"] && args["evaluate-no-iw"]
        error("Cannot set both \"--evaluate-iw\" and \"--evaluate-no-iw\".")
    elseif args["evaluate-iw"]
        println("Force using importance weighting for evaluation objective.")
        eval_importance_weighted = true
    elseif args["evaluate-no-iw"]
        println("Force NOT using importance weighting for evaluation objective.")
        eval_importance_weighted = false
    end

    # Use a high-sample log-EL for the eval loss.
    eval_loss(xs...) = NeuralProcesses.loglik(
        xs...,
        num_samples=args["evaluate-num-samples"],
        importance_weighted=eval_importance_weighted,
        fixed_σ_epochs=0
    )
else
    error("Unknown model \"" * args["model"] * "\".")
end

# Determine name of file to write model to.
if !isnothing(args["bson"])
    bson = args["bson"]
else
    bson =
        args["models-dir"] * "/" *
        args["model"] * "/" *
        args["loss"] * "/" *
        args["data"] * ".bson"
    mkpath(args["models-dir"] * "/" * args["model"] * "/" * args["loss"])
end

# Determine folder to output images.
path = "output/" * args["model"] * "/" * args["loss"] * "/" * args["data"]
mkpath("output/" * args["model"] * "/" * args["loss"] * "/" * args["data"])

function build_data_gen(; x_context, x_target, num_context, num_target, batch_size)
    return DataGenerator(
        process,
        batch_size=batch_size,
        x_context=x_context,
        x_target=x_target,
        num_context=num_context,
        num_target=num_target,
        σ²=noise
    )
end

if args["evaluate"]
    # Use the best model for evaluation.
    model = best_model(bson) |> gpu
    report_num_params(model)

    # Use a batch size of one to support a high number of samples. We alleviate the increase
    # in variance of the objective by using a high number of batches.
    batch_size = 1
    num_batches = 5000

    # Determine which evaluation tasks to perform.
    tasks = [(
        "interpolation on training range",
        build_data_gen(
            x_context=Uniform(-2, 2),
            x_target=Uniform(-2, 2),
            num_context=num_context_eval,
            num_target=num_target,
            batch_size=batch_size
        )
    )]
    if !args["evaluate-only-within"]
        push!(tasks, (
            "interpolation beyond training range",
            build_data_gen(
                x_context=Uniform(2, 6),
                x_target=Uniform(2, 6),
                num_context=num_context_eval,
                num_target=num_target,
                batch_size=batch_size
            )
        ))
        push!(tasks, (
            "extrapolation beyond training range",
            build_data_gen(
                x_context=Uniform(-2, 2),
                x_target=UniformUnion(Uniform(-4, -2), Uniform(2, 4)),
                num_context=num_context_eval,
                num_target=num_target,
                batch_size=batch_size
            )
        ))
    end

    # Perform evaluation tasks with `epoch` set to 1000.
    for (name, data_gen) in tasks
        println("Evaluation task: $name")
        eval_model(model, eval_loss, data_gen, 1000, num_batches=num_batches)
    end
else
    # Construct data generator for training.
    data_gen = build_data_gen(
        x_context=Uniform(-2, 2),
        x_target=Uniform(-2, 2),
        num_context=num_context,
        num_target=num_target,
        batch_size=args["batch-size"]
    )

    if args["starting-epoch"] > 1
        # Continue training from most recent model.
        model = recent_model(bson) |> gpu
    else
        # Instantiate a new model to start training. Ideally, the margin should be the
        # receptive field size, but that creates large memory requirements for models with
        # large receptive field.
        margin = 0f1
        if args["model"] == "convcnp"
            model = convcnp_1d(
                receptive_field=receptive_field,
                num_layers=8,
                num_channels=num_channels,
                points_per_unit=points_per_unit,
                margin=margin
            ) |> gpu
        elseif args["model"] == "corconvcnp"
            model = corconvcnp_1d(
                receptive_field_μ=receptive_field,
                receptive_field_Σ=min(receptive_field, 8f0),
                num_layers=8,
                num_channels=num_channels,
                points_per_unit_μ=points_per_unit,
                points_per_unit_Σ=20f0,
                margin=margin
            ) |> gpu
        elseif args["model"] in [
            "convnp",
            "convnp-global-mean", "convnp-global-sum",
            "convnp-amortised-mean", "convnp-amortised-sum",
            "convnp-het"
        ]
            if args["model"] == "convnp"
                num_global_channels = 0
                noise_type = "fixed"
                pooling_type = "mean"  # This doesn't matter, but must be set to something.
            elseif args["model"] == "convnp-global-mean"
                num_global_channels = 16
                noise_type = "fixed"
                pooling_type = "mean"
            elseif args["model"] == "convnp-global-sum"
                num_global_channels = 16
                noise_type = "fixed"
                pooling_type = "sum"
            elseif args["model"] == "convnp-amortised-mean"
                num_global_channels = 0
                noise_type = "amortised"
                pooling_type = "mean"
            elseif args["model"] == "convnp-amortised-sum"
                num_global_channels = 0
                noise_type = "amortised"
                pooling_type = "sum"
            elseif args["model"] == "convnp-het"
                num_global_channels = 0
                noise_type = "het"
                pooling_type = "mean"  # This doesn't matter, but must be set to something.
            else
                error("Unknown model \"" * args["model"] * "\".")
            end

            model = convnp_1d(
                receptive_field=receptive_field,
                num_encoder_layers=8,
                num_decoder_layers=8,
                num_encoder_channels=num_channels,
                num_decoder_channels=num_channels,
                num_latent_channels=16,
                num_global_channels=num_global_channels,
                points_per_unit=points_per_unit,
                margin=1f0,
                noise_type=noise_type,
                pooling_type=pooling_type,
                σ=5f-2,
                learn_σ=false
            ) |> gpu
        elseif args["model"] == "acnp"
            model = acnp_1d(
                dim_embedding=dim_embedding,
                num_encoder_heads=8,
                num_encoder_layers=6,
                num_decoder_layers=6
            ) |> gpu
        elseif args["model"] in [
            "anp",
            "anp-amortised-mean", "anp-amortised-sum",
            "anp-het"
        ]
            if args["model"] == "anp"
                noise_type = "fixed"
                pooling_type = "mean"  # This doesn't matter, but must be set to something.
            elseif args["model"] == "anp-amortised-mean"
                noise_type = "amortised"
                pooling_type = "mean"
            elseif args["model"] == "anp-amortised-sum"
                noise_type = "amortised"
                pooling_type = "sum"
            elseif args["model"] == "anp-het"
                noise_type = "het"
                pooling_type = "mean"  # This doesn't matter, but must be set to something.
            else
                error("Unknown model \"" * args["model"] * "\".")
            end

            model = anp_1d(
                dim_embedding=dim_embedding,
                num_encoder_heads=8,
                num_encoder_layers=6,
                num_decoder_layers=6,
                noise_type=noise_type,
                pooling_type=pooling_type,
                σ=5f-2,
                learn_σ=false
            ) |> gpu
        elseif args["model"] == "cnp"
            model = cnp_1d(
                dim_embedding=dim_embedding,
                num_encoder_layers=6,
                num_decoder_layers=6
            ) |> gpu
        elseif args["model"] in [
            "np",
            "np-amortised-mean", "np-amortised-sum",
            "np-het"
        ]
            if args["model"] == "np"
                noise_type = "fixed"
                pooling_type = "mean"  # This doesn't matter, but must be set to something.
            elseif args["model"] == "np-amortised-mean"
                noise_type = "amortised"
                pooling_type = "mean"
            elseif args["model"] == "np-amortised-sum"
                noise_type = "amortised"
                pooling_type = "sum"
            elseif args["model"] == "np-het"
                noise_type = "het"
                pooling_type = "mean"  # This doesn't matter, but must be set to something.
            else
                error("Unknown model \"" * args["model"] * "\".")
            end

            model = np_1d(
                dim_embedding=dim_embedding,
                num_encoder_layers=6,
                num_decoder_layers=6,
                noise_type=noise_type,
                pooling_type=pooling_type,
                σ=5f-2,
                learn_σ=false
            ) |> gpu
        else
            error("Unknown model \"" * args["model"] * "\".")
        end
    end

    report_num_params(model)

    train!(
        model,
        loss,
        data_gen,
        ADAM(5e-4),
        bson=bson,
        starting_epoch=args["starting-epoch"],
        tasks_per_epoch=2^14,
        epochs=args["epochs"],
        path=path
    )
end


text:
module PDControl

using RigidBodyDynamics.Spatial
using StaticArrays
using Rotations

# functions
export
    pd

# gain types
export
    PDGains,
    FramePDGains,
    SE3PDGains

# pd control methods
export
    SE3PDMethod

abstract type AbstractPDGains end

group_error(x, xdes) = x - xdes
group_error(x::Rotation, xdes::Rotation) = inv(xdes) * x
group_error(x::Transform3D, xdes::Transform3D) = inv(xdes) * x

pd(gains::AbstractPDGains, x, xdes, ẋ, ẋdes) = pd(gains, group_error(x, xdes), -ẋdes + ẋ)
pd(gains::AbstractPDGains, x, xdes, ẋ, ẋdes, method) = pd(gains, group_error(x, xdes), -ẋdes + ẋ, method)

struct PDGains{K, D} <: AbstractPDGains
    k::K
    d::D
end

pd(gains::PDGains, e, ė) = -gains.k * e - gains.d * ė
pd(gains::PDGains, e::RodriguesVec, ė::AbstractVector) = pd(gains, SVector(e.sx, e.sy, e.sz), ė)
pd(gains::PDGains, e::Rotation{3}, ė::AbstractVector) = pd(gains, RodriguesVec(e), ė)

rotategain(gain::Number, R::Rotation) = gain
rotategain(gain::AbstractMatrix, R::Rotation) = R * gain * R'

function Spatial.transform(gains::PDGains, tf::Transform3D)
    R = rotation(tf)
    PDGains(rotategain(gains.k, R), rotategain(gains.d, R))
end

# Gains with a frame annotation
struct FramePDGains{K, D} <: AbstractPDGains
    frame::CartesianFrame3D
    gains::PDGains{K, D}
end

function pd(fgains::FramePDGains, e::FreeVector3D, ė::FreeVector3D)
    @framecheck fgains.frame e.frame
    @framecheck fgains.frame ė.frame
    FreeVector3D(fgains.frame, pd(fgains.gains, e.v, ė.v))
end

function Spatial.transform(fgains::FramePDGains, tf::Transform3D)
    @framecheck tf.from fgains.frame
    FramePDGains(tf.to, transform(fgains.gains, tf))
end

# PD control on the special Euclidean group
struct SE3PDGains{A<:Union{PDGains, FramePDGains}, L<:Union{PDGains, FramePDGains}} <: AbstractPDGains
    angular::A
    linear::L
end

function SE3PDGains(frame::CartesianFrame3D, angular::PDGains, linear::PDGains)
    SE3PDGains(FramePDGains(frame, angular), FramePDGains(frame, linear))
end

Spatial.angular(gains::SE3PDGains) = gains.angular
Spatial.linear(gains::SE3PDGains) = gains.linear

function Spatial.transform(gains::SE3PDGains, t::Transform3D)
    SE3PDGains(transform(gains.angular, t), transform(gains.linear, t))
end

struct SE3PDMethod{T} end

pd(gains::SE3PDGains, e::Transform3D, ė::Twist) = pd(gains, e, ė, SE3PDMethod{:DoubleGeodesic}())

function pd(gains::SE3PDGains, e::Transform3D, ė::Twist, ::SE3PDMethod{:DoubleGeodesic})
    # Theorem 12 in Bullo, Murray, "Proportional derivative (PD) control on the Euclidean group", 1995 (4.6 in the Technical Report).
    # Note: in Theorem 12, even though the twist and spatial acceleration are expressed in body frame, the gain Kv is expressed in base (or desired) frame,
    # since it acts on the translation part of the transform from body frame to base frame (i.e. the definition of body frame expressed in base frame),
    # and force Kv * p needs to be rotated back to body frame to match the spatial acceleration in body frame.
    # Instead, we express Kv in body frame by performing a similarity transform on Kv: R * Kv * R', where R is the rotation from actual body frame to desired body frame.
    # This turns the linear and proportional part of the PD law into R' * R * Kv * R' * p = Kv * R' * p
    # Note also that in Theorem 12, the frame in which the orientation gain Kω must be expressed is ambiguous, since R' * log(R) = log(R).
    # This explains why the Kω used here is the same as the Kω in Theorem 12.
    # Gains should be expressed in actual body frame.
    bodyframe = ė.body
    @framecheck e.from bodyframe
    @framecheck ė.base e.to
    @framecheck ė.frame bodyframe

    R = rotation(e)
    p = translation(e)
    ψ = RodriguesVec(R)

    ang = pd(angular(gains), FreeVector3D(bodyframe, ψ.sx, ψ.sy, ψ.sz), FreeVector3D(bodyframe, angular(ė)))
    lin = pd(linear(gains), FreeVector3D(bodyframe, R' * p), FreeVector3D(bodyframe, linear(ė)))
    SpatialAcceleration(ė.body, ė.base, ang, lin)
end

function pd(gains::SE3PDGains, e::Transform3D, ė::Twist, ::SE3PDMethod{:Linearized})
    bodyframe = ė.body
    @framecheck e.from bodyframe
    @framecheck ė.base e.to
    @framecheck ė.frame bodyframe

    R = rotation(e)
    p = translation(e)
    ψ = linearized_rodrigues_vec(R)

    ang = pd(angular(gains), FreeVector3D(bodyframe, ψ.sx, ψ.sy, ψ.sz), FreeVector3D(bodyframe, angular(ė)))
    lin = pd(linear(gains), FreeVector3D(bodyframe, R' * p), FreeVector3D(bodyframe, linear(ė)))
    SpatialAcceleration(ė.body, ė.base, ang, lin)
end

end # module


text:
using Sundials, Test

println("Test error handling")

f_error(u,p,t) = u/t
u0 = 1.0
prob = ODEProblem(f_error,u0,(0.0,1.0))
sol = solve(prob,CVODE_BDF())
sol = solve(prob,CVODE_BDF(),verbose=false)

f_error2(du,u,p,t) = u/t-1
u0 = 1.0; du0 = 1.0
prob = DAEProblem(f_error2,u0,du0,(0.0,1.0),differential_vars=[1])
sol = solve(prob,IDA())
sol = solve(prob,IDA(),verbose=false)

@test sol.retcode == :InitialFailure


text:
using LSystems
using Base.Test

# Algae-growth from https://en.wikipedia.org/wiki/L-system#Example_1:_Algae
algae_start = "A"
algae_transition = Dict([('A', "AB"), ('B', "A")])

s = ""
for c in @task lindenmayer(algae_start, algae_transition, 0)
    s *= string(c)
end
@test "A" == s

s = ""
for c in @task lindenmayer(algae_start, algae_transition, 1)
    s *= string(c)
end
@test "AB" == s

s = ""
for c in @task lindenmayer(algae_start, algae_transition, 5)
    s *= string(c)
end
@test "ABAABABAABAAB" == s

s = ""
for c in @task lindenmayer(algae_start, algae_transition, 7)
    s *= string(c)
end
@test "ABAABABAABAABABAABABAABAABABAABAAB" == s

# Pythagoras tree from https://en.wikipedia.org/wiki/L-system#Example_2:_Pythagoras_tree
pythagoras_start = "0"
pythagoras_transition = Dict([('1', "11"), ('0', "1[0]0")])

s = ""
for c in @task lindenmayer(pythagoras_start, pythagoras_transition, 3)
    s *= string(c)
end
@test "1111[11[1[0]0]1[0]0]11[1[0]0]1[0]0" == s

# Koch curve from https://en.wikipedia.org/wiki/L-system#Example_4:_Koch_curve
koch_curve_start = "F"
koch_curve_transition = Dict([('F', "F+F-F-F+F")])

s = ""
for c in @task lindenmayer(koch_curve_start, koch_curve_transition, 2)
    s *= string(c)
end
@test "F+F-F-F+F+F+F-F-F+F-F+F-F-F+F-F+F-F-F+F+F+F-F-F+F" == s


text:
# https://discourse.julialang.org/t/composition-and-inheritance-the-julian-way/11231

# https://github.com/lwabeke/Unmarshal.jl

struct IrisData
    class::String
    petallength::Float64
    petalwidth::Float64
    sepallength::Float64
    sepalwidth::Float64
end


text:
using Test, Sockets, HTTP, JSON
using GraphQL

function before()
  server = GraphQLServer()
  client = GraphQLClient("http://$(server.address)")
  return Dict("server"=>server, "client"=>client)
end

function testPing(client)
  data = query(client)
  @test get(data,"TODO",nothing) === "TODO"
end

function after(server)
  close(server.server)
end

@testset "Server" begin
  @info "Running server tests..."
  fixture = before()
  testPing(fixture["client"])
  after(fixture["server"])
end


text:
function hello()
    "Hello, World!"
end

max_stars_repo_name: cdrlos/exercism

max_issues_repo_name: cdrlos/exercism

max_forks_repo_name: cdrlos/exercism


text:
# Nesting calls to apply should take the appropriate union of patches.
@testset "nested apply calls" begin
    multiply(x::Number) = 2x
    multiply(x::Int) = 2x - 1
    add(x::Number) = x + 2
    add(x::Int) = x + 1

    @testset "simple" begin
        patches = Patch[
            @patch multiply(x::Integer) = 3x
            @patch multiply(x::Int) = 4x
            @patch add(x::Int) = x + 4
        ]

        apply(patches) do
            @test (@mock multiply(2)) == 8
            @test (@mock multiply(0x2)) == 0x6
            @test (@mock multiply(2//1)) == 4//1
            @test (@mock add(2//1)) == 4//1
            @test (@mock add(2)) == 6
        end

        apply(patches[1]) do
            @test (@mock multiply(2)) == 6
            @test (@mock multiply(0x2)) == 0x6
            @test (@mock multiply(2//1)) == 4//1
            @test (@mock add(2//1)) == 4//1
            @test (@mock add(2)) == 3

            apply(patches[2]) do
                @test (@mock multiply(2)) == 8
                @test (@mock multiply(0x2)) == 0x6
                @test (@mock multiply(2//1)) == 4//1
                @test (@mock add(2//1)) == 4//1
                @test (@mock add(2)) == 3

                apply(patches[3]) do
                    @test (@mock multiply(2)) == 8
                    @test (@mock multiply(0x2)) == 0x6
                    @test (@mock multiply(2//1)) == 4//1
                    @test (@mock add(2//1)) == 4//1
                    @test (@mock add(2)) == 6
                end

                @test (@mock multiply(2)) == 8
                @test (@mock multiply(0x2)) == 0x6
                @test (@mock multiply(2//1)) == 4//1
                @test (@mock add(2//1)) == 4//1
                @test (@mock add(2)) == 3
            end

            @test (@mock multiply(2)) == 6
            @test (@mock multiply(0x2)) == 0x6
            @test (@mock multiply(2//1)) == 4//1
            @test (@mock add(2//1)) == 4//1
            @test (@mock add(2)) == 3
        end
    end

    @testset "repeated patch" begin
        patches = Patch[
            @patch multiply(x::Integer) = 3x
            @patch multiply(x::Integer) = 4x
        ]

        apply(patches) do
            @test (@mock multiply(2)) == 8
        end

        apply(patches[1]) do
            @test (@mock multiply(2)) == 6
            apply(patches[2]) do
                @test (@mock multiply(2)) == 8
            end
            @test (@mock multiply(2)) == 6
        end
    end
end


text:
function RBMs.log_pseudolikelihood(rbm::ConvRBM, v::AbstractArray; β::Real=true)
    return log_pseudolikelihood_stoch(rbm, v; β)
end

function log_pseudolikelihood_stoch(rbm::ConvRBM, v::AbstractArray; β::Real=true)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    sites = [
        rand(CartesianIndices(sitesize(rbm)))
        for _ in CartesianIndices(batch_size(rbm, v))
    ]
    return log_pseudolikelihood_sites(rbm, v, sites; β)
end

sitesize(rbm::ConvRBM) = size(rbm.w)[1:(end - ndims(rbm.hidden))]
sitesize(rbm::ConvRBM{<:RBMs.Potts}) = size(rbm.w)[2:(end - ndims(rbm.hidden))]

function log_pseudolikelihood_sites(
    rbm::ConvRBM{<:Potts},
    v::AbstractArray,
    sites::AbstractArray{<:CartesianIndex};
    β::Real=true
)
    ΔE = substitution_matrix_sites(rbm, v, sites; β)
    lPL = -logsumexp(-β * ΔE; dims=1)
    return reshape(lPL, batch_size(rbm, v))
end

function substitution_matrix_sites(
    rbm::ConvRBM{<:Potts},
    v::AbstractArray,
    sites::AbstractArray{<:CartesianIndex};
    β::Real = true
)
    E_ = zeros(RBMs.colors(visible(rbm)), batch_size(rbm, v)...)
    for x in 1:RBMs.colors(visible(rbm))
        v_ = copy(v)
        for (b, i) in pairs(sites)
            v_[:, i, b] .= false
            v_[x, i, b] = true
        end
        selectdim(E_, 1, x) .= free_energy(rbm, v_; β)
    end
    c = RBMs.onehot_decode(v)
    E = [E_[c[i, b], b] for (b, i) in pairs(sites)]
    return E_ .- reshape(E, 1, batch_size(rbm, v)...)
end

function log_pseudolikelihood_sites(
    rbm::ConvRBM{<:RBMs.Binary},
    v::AbstractArray,
    sites::AbstractArray{<:CartesianIndex};
    β::Real = true
)
    v_ = copy(v)
    for (b, i) in pairs(sites)
        v_[i, b] = 1 - v_[i, b]
    end
    F = free_energy(rbm, v; β)
    F_ = free_energy(rbm, v_; β)
    return -log1pexp.(β * (F - F_))
end

function log_pseudolikelihood_sites(
    rbm::ConvRBM{<:RBMs.Spin},
    v::AbstractArray,
    sites::AbstractVector{<:CartesianIndex};
    β::Real = true
)
    v_ = copy(v)
    for (b, i) in pairs(sites)
        v_[i, b] = -v_[i, b]
    end
    F = free_energy(rbm, v; β)
    F_ = free_energy(rbm, v_; β)
    return -log1pexp.(β * (F - F_))
end


text:
struct FilterHierarchicalBatchQs <: MultiObjectiveBatchQs
    model::SVDD.OCClassifier
    inf_measure::SequentialPQs
    rep_measure::Function
    div_measure::Function
    k::Int

    function FilterHierarchicalBatchQs(model::SVDD.OCClassifier, informativeness::SequentialPQs; representativeness::Symbol, diversity::Symbol,
        k::Int)::FilterHierarchicalBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))

        representativeness_measure = get_rep_measure(representativeness)
        diversity_measure = get_incremental_div_measure(diversity)
        return new(model, informativeness, representativeness_measure, diversity_measure, k)
    end
end

"""
Idea: Find 4*k elements with largest representativeness,
    select 2*k elements with largest informativeness,
    iteratively generate batch by greedyly adding sample with maximum diversity
        to already selected batch samples
"""
function select_batch(qs::FilterHierarchicalBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    # representativeness
    rep_scores = qs.rep_measure(qs.model, x, labels, candidate_indices)
    most_representative_indices = candidate_indices[sortperm(rep_scores, rev=true)[1:min(4 * qs.k, num_observations)]]

    #informativeness
    inf_scores = qs_score(qs.inf_measure, x, labels)[most_representative_indices]
    batch_candidate_indices = most_representative_indices[sortperm(inf_scores, rev=true)[1:min(2 * qs.k, num_observations)]]

    batch_samples = batch_candidate_indices[1:1]
    div_scores = Float64[]
    for iteration in 2:qs.k
        div_scores = qs.div_measure(qs.model, x, batch_samples[end], candidate_indices, div_scores)
        best_sample_index = candidate_indices[argmax(div_scores)]

        push!(batch_samples, best_sample_index)
    end
    return batch_samples
end