text:
function compute_routing(rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::MinimumTotal, ::PathFormulation, ::Val{true}, ::Automatic, ::NoUncertaintyHandling, ::NoUncertainty)
    _ensure_objective_compatibility(edge_obj, agg_obj)

    start = time_ns()

    # Create the problem.
    m = Model(rd.solver)
    set_silent(m)
    rm = basic_routing_model_unitary(m, rd)

    # Constraints.
    capacity_constraints(rm, rd.traffic_matrix)

    # Objective function.
    obj = sum(objective_edge_expression(rm, edge_obj, e) for e in edges(rd))
    @objective(m, Min, obj)

    time_create_master_model_ms = (time_ns() - start) / 1_000_000.

    # Solve this by column generation.
    routings = Routing[]
    objectives = Float64[]
    times_ms = Float64[]
    n_new_paths = 0
    n_iter = 0

    while true
        start = time_ns()

        # Optimise for the current iteration.
        optimize!(m)

        push!(routings, Routing(rd, value.(rm.routing)))
        push!(objectives, objective_value(m))

        # First, the pricing.
        dual_value_convexity = dual.(rm.constraints_convexity) # Demand -> dual value
        dual_values_capacity = Dict(e => dual(rm.constraints_capacity[e]) for e in edges(rd)) # Edge -> dual value

        new_paths = _mintotload_solve_pricing_problem_master(rd, dual_values_capacity, dual_value_convexity)
        n_new_paths += length(new_paths)

        # Add the new columns to the formulation.
        for (d, path) in new_paths
            # Add the new columns to the data object and retrieve a path ID.
            path_id = add_path!(rm, d, path)

            # Get your hands dirty in the formulation.
            # - Create a new variable for this path, add it into the
            #   convexity and capacity constraints.
            add_routing_var!(rm, demand, path_id, constraint_capacity=true)
        end

        rd.logmessage("[$(n_iter)] Added $(length(new_paths)) new paths.")
        n_iter += 1
        push!(times_ms, (time_ns() - start) / 1_000_000.)

        if length(new_paths) == 0
            break
        end
    end

    return RoutingSolution(rd,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_create_master_model_ms=time_create_master_model_ms,
                           time_solve_ms=times_ms,
                           objectives=objective_value(m),
                           n_columns=n_new_paths,
                           matrices=rd.traffic_matrix,
                           routings=Routing(rd, value.(rm.routing)),
                           master_model=rm)
end

function _mintotload_solve_pricing_problem_master(rd::RoutingData,
                                                  dual_values_capacity,
                                                  dual_value_convexity)
    paths = Dict{Edge{Int}, Vector{Edge}}()

    for demand in demands(rd)
        # Determine the weight matrix to use for the computations:
        #     for e: \sum_{D\in\Delta} dual_{e, d_k} d_k
        rd.logmessage("Pricing for $demand")
        weight_matrix = spzeros(Float64, n_nodes(rd), n_nodes(rd))
        for e in edges(rd)
            value = 1.0 / capacity(rd, e) - rd.traffic_matrix[demand] * dual_values_capacity[e]
            if value >= - CPLEX_REDUCED_COST_TOLERANCE
                weight_matrix[src(e), dst(e)] = value
            end
        end

        # Check that at least some weights are nonzero.
        if nnz(weight_matrix) == 0
            # No need for the potentially expensive _check_weight_matrix,
            # as no value in weight_matrix can be below the tolerance by
            # construction (no addition, unlike iterative version of oblivious).
            return paths
        end

        # Compute the shortest path for this demand.
        state = desopo_pape_shortest_paths(rd.g, src(demand), weight_matrix)

        # Is there a solution? Has this path a negative reduced cost?
        if all(state.parents .== 0)
           continue
        end
        reduced_cost = state.dists[dst(demand)] - dual_value_convexity[demand]
        if reduced_cost >= - CPLEX_REDUCED_COST_TOLERANCE
           continue
        end
        rd.logmessage(reduced_cost)

        # Build the path.
        path = _build_path(state, demand)
        if path in rd.paths_edges
           rd.logmessage("New path $(path) is already in the formulation!")
        else
           paths[demand] = path
        end
    end

    return paths
end

function compute_routing(rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::MinimumTotal, ::PathFormulation, ::Val{true}, ::Automatic, ::NoUncertaintyHandling, ::NoUncertainty)
    start = time_ns()

    # Create the problem.
    m = _create_model(rd)
    rm = basic_routing_model_unitary(m, rd)

    # Objective function.
    @variable(m, mu >= 0)
    rm.mu = mu
    @objective(m, Min, mu)

    # Constraints.
    mu_capacity_constraints(rm, rd.traffic_matrix)

    time_create_master_model_ms = (time_ns() - start) / 1_000_000.

    # Solve this by column generation.
    routings = Routing[]
    objectives = Float64[]
    times_ms = Float64[]
    n_new_paths = 0
    n_iter = 0

    while true
        start = time_ns()

        # Optimise for the current iteration.
        optimize!(m)

        push!(routings, Routing(rd, value.(rm.routing)))
        push!(objectives, objective_value(m))

        # First, the pricing.
        dual_value_convexity = dual.(rm.constraints_convexity) # Demand -> dual value
        dual_values_capacity = Dict(e => dual(rm.constraints_capacity[e]) for e in edges(rd)) # Edge -> dual value

        new_paths = _minmaxload_solve_pricing_problem_master(rd, dual_values_capacity, dual_value_convexity)
        n_new_paths += length(new_paths)

        # Add the new columns to the formulation.
        for (d, path) in new_paths
            # Add the new columns to the data object and retrieve a path ID.
            path_id = add_path!(rm, d, path)

            # Get your hands dirty in the formulation.
            # - Create a new variable for this path, add it into the
            #   convexity and capacity constraints.
            add_routing_var!(rm, demand, path_id, constraint_capacity=true)
        end

        rd.logmessage("[$(n_iter)] Added $(length(new_paths)) new paths.")
        n_iter += 1
        push!(times_ms, (time_ns() - start) / 1_000_000.)

        if length(new_paths) == 0
            break
        end
    end

    return RoutingSolution(rd,
                           time_create_master_model_ms=time_create_master_model_ms,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_solve_ms=(stop - start) / 1_000_000.,
                           objectives=objective_value(m),
                           n_columns=n_new_paths,
                           matrices=rd.traffic_matrix,
                           routings=Routing(rd, value.(rm.routing)),
                           master_model=rm)
end

function _minmaxload_solve_pricing_problem_master(rd::RoutingData,
                                                  dual_values_capacity,
                                                  dual_value_convexity)
    paths = Dict{Edge{Int}, Vector{Edge}}()

    for demand in demands(rd)
        # Determine the weight matrix to use for the computations:
        #     for e: \sum_{D\in\Delta} dual_{e, d_k} d_k
        rd.logmessage("Pricing for $demand")
        weight_matrix = spzeros(Float64, n_nodes(rd), n_nodes(rd))
        for e in edges(rd)
            value = 1.0 / capacity(rd, e) - rd.traffic_matrix[demand] * dual_values_capacity[e]
            if value >= - CPLEX_REDUCED_COST_TOLERANCE
                weight_matrix[src(e), dst(e)] = value
            end
        end

        # Check that at least some weights are nonzero.
        if nnz(weight_matrix) == 0
            # No need for the potentially expensive _check_weight_matrix,
            # as no value in weight_matrix can be below the tolerance by
            # construction (no addition, unlike iterative version of oblivious).
            return paths
        end

        # Compute the shortest path for this demand.
        state = desopo_pape_shortest_paths(rd.g, src(demand), weight_matrix)

        # Is there a solution? Has this path a negative reduced cost?
        if all(state.parents .== 0)
           continue
        end
        reduced_cost = state.dists[dst(demand)] - dual_value_convexity[demand]
        if reduced_cost >= - CPLEX_REDUCED_COST_TOLERANCE
           continue
        end
        rd.logmessage(reduced_cost)

        # Build the path.
        path = _build_path(state, demand)
        if path in rd.paths_edges
           rd.logmessage("New path $(path) is already in the formulation!")
        else
           paths[demand] = path
        end
    end

    return paths
end


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

using AWS.Compat
using AWS.UUIDs
"""
    AssociateMemberAccount()

Associates a specified AWS account with Amazon Macie Classic as a member account.

# Required Parameters
- `memberAccountId`: The ID of the AWS account that you want to associate with Amazon Macie Classic as a member account.

"""

associate_member_account(memberAccountId; aws_config::AWSConfig=global_aws_config()) = macie("AssociateMemberAccount", Dict{String, Any}("memberAccountId"=>memberAccountId); aws_config=aws_config)
associate_member_account(memberAccountId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("AssociateMemberAccount", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("memberAccountId"=>memberAccountId), args)); aws_config=aws_config)

"""
    AssociateS3Resources()

Associates specified S3 resources with Amazon Macie Classic for monitoring and data classification. If memberAccountId isn't specified, the action associates specified S3 resources with Macie Classic for the current master account. If memberAccountId is specified, the action associates specified S3 resources with Macie Classic for the specified member account.

# Required Parameters
- `s3Resources`: The S3 resources that you want to associate with Amazon Macie Classic for monitoring and data classification.

# Optional Parameters
- `memberAccountId`: The ID of the Amazon Macie Classic member account whose resources you want to associate with Macie Classic.
"""

associate_s3_resources(s3Resources; aws_config::AWSConfig=global_aws_config()) = macie("AssociateS3Resources", Dict{String, Any}("s3Resources"=>s3Resources); aws_config=aws_config)
associate_s3_resources(s3Resources, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("AssociateS3Resources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("s3Resources"=>s3Resources), args)); aws_config=aws_config)

"""
    DisassociateMemberAccount()

Removes the specified member account from Amazon Macie Classic.

# Required Parameters
- `memberAccountId`: The ID of the member account that you want to remove from Amazon Macie Classic.

"""

disassociate_member_account(memberAccountId; aws_config::AWSConfig=global_aws_config()) = macie("DisassociateMemberAccount", Dict{String, Any}("memberAccountId"=>memberAccountId); aws_config=aws_config)
disassociate_member_account(memberAccountId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("DisassociateMemberAccount", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("memberAccountId"=>memberAccountId), args)); aws_config=aws_config)

"""
    DisassociateS3Resources()

Removes specified S3 resources from being monitored by Amazon Macie Classic. If memberAccountId isn't specified, the action removes specified S3 resources from Macie Classic for the current master account. If memberAccountId is specified, the action removes specified S3 resources from Macie Classic for the specified member account.

# Required Parameters
- `associatedS3Resources`: The S3 resources (buckets or prefixes) that you want to remove from being monitored and classified by Amazon Macie Classic.

# Optional Parameters
- `memberAccountId`: The ID of the Amazon Macie Classic member account whose resources you want to remove from being monitored by Amazon Macie Classic.
"""

disassociate_s3_resources(associatedS3Resources; aws_config::AWSConfig=global_aws_config()) = macie("DisassociateS3Resources", Dict{String, Any}("associatedS3Resources"=>associatedS3Resources); aws_config=aws_config)
disassociate_s3_resources(associatedS3Resources, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("DisassociateS3Resources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("associatedS3Resources"=>associatedS3Resources), args)); aws_config=aws_config)

"""
    ListMemberAccounts()

Lists all Amazon Macie Classic member accounts for the current Amazon Macie Classic master account.

# Optional Parameters
- `maxResults`: Use this parameter to indicate the maximum number of items that you want in the response. The default value is 250.
- `nextToken`: Use this parameter when paginating results. Set the value of this parameter to null on your first call to the ListMemberAccounts action. Subsequent calls to the action fill nextToken in the request with the value of nextToken from the previous response to continue listing data.
"""

list_member_accounts(; aws_config::AWSConfig=global_aws_config()) = macie("ListMemberAccounts"; aws_config=aws_config)
list_member_accounts(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("ListMemberAccounts", args; aws_config=aws_config)

"""
    ListS3Resources()

Lists all the S3 resources associated with Amazon Macie Classic. If memberAccountId isn't specified, the action lists the S3 resources associated with Amazon Macie Classic for the current master account. If memberAccountId is specified, the action lists the S3 resources associated with Amazon Macie Classic for the specified member account.

# Optional Parameters
- `maxResults`: Use this parameter to indicate the maximum number of items that you want in the response. The default value is 250.
- `memberAccountId`: The Amazon Macie Classic member account ID whose associated S3 resources you want to list.
- `nextToken`: Use this parameter when paginating results. Set its value to null on your first call to the ListS3Resources action. Subsequent calls to the action fill nextToken in the request with the value of nextToken from the previous response to continue listing data.
"""

list_s3_resources(; aws_config::AWSConfig=global_aws_config()) = macie("ListS3Resources"; aws_config=aws_config)
list_s3_resources(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("ListS3Resources", args; aws_config=aws_config)

"""
    UpdateS3Resources()

Updates the classification types for the specified S3 resources. If memberAccountId isn't specified, the action updates the classification types of the S3 resources associated with Amazon Macie Classic for the current master account. If memberAccountId is specified, the action updates the classification types of the S3 resources associated with Amazon Macie Classic for the specified member account.

# Required Parameters
- `s3ResourcesUpdate`: The S3 resources whose classification types you want to update.

# Optional Parameters
- `memberAccountId`: The AWS ID of the Amazon Macie Classic member account whose S3 resources' classification types you want to update.
"""

update_s3_resources(s3ResourcesUpdate; aws_config::AWSConfig=global_aws_config()) = macie("UpdateS3Resources", Dict{String, Any}("s3ResourcesUpdate"=>s3ResourcesUpdate); aws_config=aws_config)
update_s3_resources(s3ResourcesUpdate, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = macie("UpdateS3Resources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("s3ResourcesUpdate"=>s3ResourcesUpdate), args)); aws_config=aws_config)

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Shapes (Modia3D/Shapes/_module.jl)
#
# Content:
#   This file contains elements to show text in the visualization window

const fontFamilyPalette = Dict{String,String}([("FreeSans"     , "FreeSans"),
                                               ("FreeSerif"    , "FreeSerif"),
                                               ("Arial"        , "ARIAL"),
                                               ("ArialNarrow"  , "ARIALN"),
                                               ("CourierNew"   , "COUR"),
                                               ("TimesNewRoman", "Times"),
                                               ("Verdana"      , "verdana")
                                               ])

function getFontFileName(fontFamily::String, bold::Bool, italic::Bool)::String
    # Handle special cases first

    if fontFamily == "ArialNarrow" && bold && !italic
        fontFileName = "ARIALNB"
    elseif fontFamily == "Verdana"
        fontFileName = bold ? ( italic ? "verdanaz" : "verdanab") :
                                ( italic ? "verdanai" : "verdana" )
    else
        if !haskey(fontFamilyPalette, fontFamily)
            error("\nError from Modia3D.visualizer: fontFamily = \"", fontFamily, "\" is not supported.\n",
                "Supported font families: FreeSans, FreeSerif, Arial, ArialNarrow, CourierNew, TimesNewRoman, Verdana.")
        end
        name = fontFamilyPalette[fontFamily]
        fontFileName = bold ? ( italic ? name * "BI" : name * "BD" ) :
                                ( italic ? name * "I"  : name )
    end
    fontFileName = fontFileName * ".TTF"
end

"""
    Font(; fontFamily="FreeSans", bold=false, italic=false,
           charSize=0.1, color=defaultColor(), transparency=0.0)

Generate a Font.

# Arguments
- `fontFamily`: name of the font family
    - "FreeSans", "FreeSerif", "Arial", "ArialNarrow", "CourierNew", "TimesNewRoman", or "Verdana"
![fontFamily](../../resources/images/TextFonts.png)
Courtesy to DLR Visualization Library
- `bold`: if true, the text is displayed in bold font.
- `italic`: if true, the text is displayed in italic font.
- `charSize`: character size in [m].
- `color` defines the material color as a String or RGB values. In case of Strings the color
          is constructed using [Colors.jl](https://github.com/JuliaGraphics/Colors.jl). So
          [named colors](https://juliagraphics.github.io/Colors.jl/stable/namedcolors/) and
          many other kinds of [color specifications](https://juliagraphics.github.io/Colors.jl/stable/constructionandconversion/#Color-Parsing)
          are available. In addition, colors can be directly defined by a vector of RGB values.
          E.g. `color="MediumVioletRed"` or `color=[255, 0, 0]`.
- `transparency`: transparency of the fon , 0.0 (opaque) ... 1.0 (transparent)

# Notes
- `fontFamily` like "Arial", "Verdana" are only supported under Windows.

# Examples
```julia
import Modia3D
font1 = Font()
font2 = Font(fontFamily="Arial", bold=true, charSize=0.2,
             color="LightBlue", transparency=0.5)
```
"""
struct Font
    fontFamily::String
    bold::Bool
    italic::Bool
    charSize::Float64
    color::MVector{3,Cint}
    transparency::Float64

    # Constructed by constructor
    fontFileName::String  # File name of font

    function Font(; fontFamily::String = "FreeSans",
        bold::Bool            = false,
        italic::Bool          = false,
        charSize::Number      = 0.1,
        color                 = defaultColor() ,
        transparency::Number  = 0.0)
        @assert(charSize >= 0.0)
        @assert(transparency >= 0.0 && transparency <= 1.0)
        fontFileName = getFontFileName(fontFamily, bold, italic)

        new(fontFamily, bold, italic, charSize, rgb(color), transparency, fontFileName)
    end
end

@enum AxisAlignment::Cint Screen=0 XY_Plane=1 XZ_Plane=2 YZ_Plane=3
@enum Alignment::Cint     Left=1 Right=2 Center=3

"""
    TextShape(; text="Hello world", font=Font(), offset=[0.0,0.0,0.0],
                axisAlignment=Screen, alignment=Center)

Generate a new visual shape representing a text shape.

# Arguments
- `text::AbstractString`: String of the text
- `font::`[`Font`](@ref): [`Font`](@ref) of the text
- `offset`: a 3D vector as offset from origin to text alignment point
- `axisAlignment::Modia3D.AxisAlignment`: alignment of text
   - `Modia3D.Screen`: parallel to screen
   - `Modia3D.XY_Plane`: in xy-planes of Object3D
   - `Modia3D.XZ_Plane`: in xz-planes of Object3D
   - `Modia3D.YZ_Plane`: in yz-planes of Object3D
- `axisAlignment::Modia3D.Alignment`: defines the direction the text is displayed, relative to its origin
  - `Modia3D.Left`, `Modia3D.Right` or `Modia3D.Center`

# Notes
- TextShape is not supported by the community edition of SimVis.
- TextShape is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).

# Examples
```julia
import Modia3D
font  = Font(fontFamily="Arial", charSize=0.4, color="MediumVioletRed")
text1 = TextShape("This is a box")
text2 = TextShape("This is the xy plane";
                   font = font, axisAlignment=Modia3D.XY_Plane,
                   alignment = Modia3D.Left)
```
"""
mutable struct TextShape <: Modia3D.AbstractShape
    text::String                   # String to be displayed
    font::Font                     # Font
    offset::MVector{3,Float64}     # Offset from origin to text alignment point
    axisAlignment::AxisAlignment   # Alignment of Text (parallel to screen or in planes of frame)
    alignment::Alignment           # Alignment of Text relative to its origin

    function TextShape(; text::AbstractString,
        font::Font=Font(),
        offset::AbstractVector = Basics.zeroMVector(),
        axisAlignment::AxisAlignment=Screen,
        alignment::Alignment=Center)

        new(text,font,offset,axisAlignment,alignment)
    end
end


text:
struct ArmijoWork
	xtilde::Vector{Float64}
end

ArmijoWork(n::Int) = ArmijoWork(Array{Float64}(undef, n))

struct ExactLinesearchWork
	tmp_n1::Vector{Float64}
	tmp_n2::Vector{Float64}
	tmp_n3::Vector{Float64}
	tmp_n4::Vector{Float64}
end

ExactLinesearchWork(n::Int) = ExactLinesearchWork([Array{Float64}(undef, n) for i in 1:4]...)

#=
Performs Armijo backtracking linesearch

OUTPUT
flag - any flags from the retraction or Armijo error
iter1 - from retraction
iter2 - from retraction
newf - new function value at xnew
f_diff - difference in function value from previous iterate
step_diff - norm of x difference from pervious iterate

OVERWRITES
xnew - with new x

=#
function armijo!(xnew, x, n, d, g, f, fval, retract_method, cval, c!, param::LFPSQPParams, work::ArmijoWork)
	f_diff = Inf
	step_diff = Inf
	α = param.α
	flag = 0
	tot_iter1 = 0
	tot_iter2 = 0
	newf = 0.0

	ar_dot = dot(d, g)

	# unpack
	xtilde = work.xtilde
	step = xtilde				# reuse xtilde for step size calculation at end

	while step_diff > param.ϵ_x
		xtilde .= x .+ α.*d

		# perform retraction
		flag, iter1, iter2 = retract!(cval, xnew, c!, xtilde, x, retract_method)
		tot_iter1 += iter1
		tot_iter2 += iter2

		# if there's an issue, shrink the step size
		if flag > 0
			α *= param.s
			continue
		end

		# calculate new function values
		step .= xnew .- x
		newf = f(xnew)

		step_diff = norm(view(step, 1:n))	# only use first n values for step_diff
		f_diff = abs(newf - fval)

		# break conditions
		if param.disable_linesearch
			break
		end

		# Armijo-Goldstein condition
		if (newf - fval) <= param.σ * α * ar_dot
			break
		end

		α *= param.s

		# to prevent infinite loop
		if α < 1e-100
			flag = 99
			break
		end
	end

	return flag, tot_iter1, tot_iter2, newf, f_diff, step_diff, α
end

#=
Performs "exact" golden-section bisection with upper bounding from McCormick

OUTPUT
flag - any flags from the retraction or Armijo error
iter1 - from retraction
iter2 - from retraction
newf - new function value at xnew
f_diff - difference in function value from previous iterate
step_diff - norm of x difference from pervious iterate

OVERWRITES
xnew - with new x

=#
function exact_linesearch!(xnew, x, n, d, f, fval, retract_method, cval, c!, param::LFPSQPParams, work::ExactLinesearchWork)

	ϕ1 = (3 - sqrt(5))/2
	ϕ2 = (sqrt(5) - 1)/2
	ϕ3 = (sqrt(5) + 1)/2
	Δ = param.α     		# use previous α as step length guess

	flag = 0
	tot_iter1 = 0
	tot_iter2 = 0
	newf = 0.0

	f_a = 0.0
	f_b = 0.0
	f_c = 0.0
	f_d = 0.0

	α_a = 0.0
	α_b = 0.0
	α_c = 0.0
	α_d = 0.0

	# unpack work and assign to x_{a,b,c,d} values
	x_a = work.tmp_n1
	x_b = work.tmp_n2
	x_c = work.tmp_n3
	x_d = work.tmp_n4
	step = work.tmp_n1		# reuse tmp_n1 for step size calculation at end

	do_shrinking = true
	flag = 0
	exact_iter_count = 0

	# find the upper bound and maintain rotating list of points
	x_d .= x
	f_d = fval

	# growing
	while true
		exact_iter_count += 1

		swp = x_b	# rotate all of the position vectors
		x_b = x_c
		x_c = x_d
		x_d = swp

		f_b = f_c
		f_c = f_d
		α_b = α_c
		α_c = α_d

		x_d .= x .+ (α_d + Δ).*d

		# perform retraction
		flag, iter1, iter2 = retract!(cval, xnew, c!, x_d, x, retract_method)
		tot_iter1 += iter1
		tot_iter2 += iter2

		x_d .= xnew 	# update x_d with retracted step

		α_d += Δ

		# should break if projection failed or if α > 1.0
		if flag > 0 || α_d > 1.0
			f_d = Inf 		# Inf indicates projection failed
			break
		end

		f_d = f(x_d)

		if f_d > f_c
			break
		end

		do_shrinking = false
		Δ *= ϕ3
	end

	# no feasible upper bound found, so perform shrinking
	if do_shrinking
		# shrinking

		# assign b with α=0 so that point a is assigned correctly below
		f_b = fval
		α_b = 0.0
		x_b .= x

		f_c = Inf
		α_c = Δ

		swp = x_d	# swap with d since d contains the point for α=Δ already
		x_d = x_c
		x_c = swp

		while true
			exact_iter_count += 1

			swp = x_d	# swapping is now "up" toward d since Δ is getting shrunk
			x_d = x_c
			x_c = swp

			f_d = f_c
			α_d = α_c

			x_c .= x .+ (ϕ1*α_c).*d

			# perform retraction
			flag, iter1, iter2 = retract!(cval, xnew, c!, x_c, x, retract_method)
			tot_iter1 += iter1
			tot_iter2 += iter2

			x_c .= xnew 	# update x_c with retracted step

			α_c *= ϕ1

			# should break if projection failed or if α > 1.0
			if flag > 0 || α_c > 1.0
				f_c = Inf
			else
				f_c = f(x_c)
			end

			if f_c <= fval || α_c < 1e-100
				break
			end
		end
	end

	# assign values from upper bounding procedure and calculate point c
	exact_iter_count += 1
	f_a = f_b
	f_b = f_c
	α_a = α_b
	α_b = α_c

	swp = x_a
	x_a = x_b
	x_b = x_c
	x_c = swp

	α_c = α_a + ϕ2*(α_d - α_a)
	x_c .= x .+ α_c.*d

	# perform retraction
	flag, iter1, iter2 = retract!(cval, xnew, c!, x_c, x, retract_method)
	tot_iter1 += iter1
	tot_iter2 += iter2

	x_c .= xnew 	# update x_c with retracted step

	if flag > 0 || α_c > 1.0
		f_c = Inf
	else
		f_c = f(x_c)
	end

	# ----------- main loop ------------------

	# do golden ratio bisection
	nd = norm(d)
	while (α_c - α_b) > 1e-6*nd
		exact_iter_count += 1

		if f_b < f_c || isinf(f_c)	# shrink interval to the left
			swp = x_d
			x_d = x_c
			x_c = x_b
			x_b = swp

			f_d = f_c
			f_c = f_b
			α_d = α_c
			α_c = α_b

			# calculate point b (which can never be infinite)
			α_b = α_a + ϕ1*(α_d - α_a)
			x_b .= x .+ α_b.*d

			# perform retraction
			flag, iter1, iter2 = retract!(cval, xnew, c!, x_b, x, retract_method)
			tot_iter1 += iter1
			tot_iter2 += iter2

			x_b .= xnew 	# update x_b with retracted step

			f_b = f(x_b)
		else 						# shrink interval to the right
			swp = x_a
			x_a = x_b
			x_b = x_c
			x_c = swp

			f_a = f_b
			f_b = f_c
			α_a = α_b
			α_b = α_c

			# calcaulte point c
			α_c = α_a + ϕ2*(α_d - α_a)
			x_c .= x .+ α_c.*d

			# perform retraction
			flag, iter1, iter2 = retract!(cval, xnew, c!, x_c, x, retract_method)
			tot_iter1 += iter1
			tot_iter2 += iter2

			x_c .= xnew 	# update x_c with retracted step

			if flag > 0 || α_c > 1.0
				f_c = Inf
			else
				f_c = f(x_c)
			end
		end
	end

	# assign final function, x, and α value
	if f_b < f_c
		xnew .= x_b
		newf = f_b
		α = α_b
	else
		xnew .= x_c
		newf = f_c
		α = α_c
	end

	step .= xnew .- x
	step_diff = norm(view(step, 1:n))	# only use first n values for step_diff
	f_diff = abs(newf - fval)

	return flag, tot_iter1, tot_iter2, newf, f_diff, step_diff, α
end


text:
"""
    mes = meshmeasures(p, t)

compute areas of every element
"""
function meshareas(p, t)
    d12 = p[t[:, 2], :] - p[t[:, 1], :]
    d13 = p[t[:, 3], :] - p[t[:, 1], :]
    mes = (d12[:, 1] .* d13[:, 2] .- d12[:, 2] .* d13[:, 1]) / 2.
    return reshape(abs.(mes), length(mes))
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
   IK, JK, SK, IF, SF = assembKM_vemKsource(pv, cellsb, rhs, boundary_condition,
                                            meshboundary)

Assemble stiffness and source terms of virtual element approach
"""
function assembKM_vemKsource(pv, cellsb, rhs, boundary_condition,
                             meshboundary)
    # type of input (for overloaded differentiation calls)
    typein = typeof(pv[1])
    n_dofs, n_polys = size(pv, 1), 3 # method has 1 degree of freedom per vertex
    # stiffness coefficients and source coefficients
    IK, JK, IF = Int64[], Int64[], Int64[]
    SK, SF = if typein == Float64
        Float64[], Float64[]
    else
        Any[], Any[]
    end
    # impose an ordering on the linear polynomials
    linear_polynomials = [[0,0], [1,0], [0,1]]
    # a utility function for wrapping around a vector
    mod_wrap(x, a) = mod(x - 1, a) + 1
    for el_id = 1:length(cellsb)
        #println(cellsb[el_id])
        vert_ids = cellsb[el_id]   # IDs of the vertices of this element
        verts = pv[vert_ids, :]    # coordinates of the vertices of this element
        n_sides = length(vert_ids) # start computing the geometric information
        area_components = verts[:, 1] .* [verts[2:end, 2]..., verts[1, 2]] .-
                         [verts[2:end, 1]..., verts[1, 1]] .* verts[:, 2]
        area = 0.5 * abs(sum(area_components))
        centroid = sum((verts .+ vcat(verts[2:end, :], verts[1, :]')) .*
                        repeat(area_components, 1, 2), dims = 1) / (6 * area)
        diameter = 0 # compute the diameter by looking at every pair of vertices
        for i = 1:(n_sides - 1)
            for j = (i + 1):n_sides
                diameter = max(diameter, norm(verts[i, :] .- verts[j, :]))
            end
        end
        D = zeros(typein, n_sides, n_polys)
        D[:, 1] .= 1
        B = zeros(typein, n_polys, n_sides)
        B[1, :] .= 1 / n_sides
        for vertex_id = 1:n_sides
            vert = verts[vertex_id, :] # this vertex and its neighbours
            prev = verts[mod_wrap(vertex_id - 1, n_sides), :]
            next = verts[mod_wrap(vertex_id + 1, n_sides), :]
            # average of normals on edges
            vertex_normal = [next[2] - prev[2], prev[1] - next[1]]
            # only need to loop over non-constant polynomials
            for poly_id = 2:n_polys
                poly_degree = linear_polynomials[poly_id]
                # gradient of a linear polynomial is constant
                monomial_grad = poly_degree / diameter
                D[vertex_id, poly_id] = dot(vert .- centroid[:],
                                            poly_degree) / diameter
                B[poly_id, vertex_id] = 0.5 * dot(monomial_grad, vertex_normal)
            end
        end
        # compute the local Ritz projector to polynomials
        projector = (B * D) \ B
        stabilising_term = (eye(n_sides) .- D * projector)' *
                           (eye(n_sides) .- D * projector)
        G = B * D
        G[1, :] .= 0
        local_stiffness = projector' * G * projector + stabilising_term
        # copy local to global
        # K[vert_ids, vert_ids] += local_stiffness
        # F[vert_ids] .+= rhs(centroid) * area / n_sides
        append!(IK, repeat(vert_ids, length(vert_ids)))
        append!(JK, repeat(vert_ids', length(vert_ids))[:])
        append!(SK, local_stiffness[:])
        append!(IF, vert_ids)
        append!(SF, repeat(rhs(centroid) * area / n_sides, length(vert_ids)))
    end
    # add boundary contribution in the source term
    # old term F -= K[:, meshboundary] * boundary_vals
    boundary_vals = boundary_condition(pv[meshboundary, :])
    Jinmesh = indexin(JK, meshboundary)
    Jboundary = findall(Jinmesh .!= nothing)
    for l ∈ Jboundary
        push!(IF, IK[l])
        push!(SF, - SK[l] * boundary_vals[Jinmesh[l]])
    end
    return IK, JK, SK, IF, SF
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
   I, J = indKM_sparse(t)

non zero indices in stiffness and mass P1 sparse matrices
"""
function indKM_sparse(t)
    it1, it2, it3 = t[:, 1], t[:, 2], t[:, 3]
    I = vcat(it1, it2, it3, it2, it3, it1, it1, it2, it3)
    J = vcat(it2, it3, it1, it1, it2, it3, it1, it2, it3)
    return I, J
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
     K, M = assembKM_P12D(p, t)

assemble the P1 stifness and mass matrices.
"""
function assembKM_P12D(p, t)
    np, dim, nt = size(p, 1), size(p, 2), size(t, 1)
    # compute ar add define arK and arM (not used ?)
    ar = meshareas(p, t)
    arK, arM = ar, ar
    it1, it2, it3 = t[:, 1], t[:, 2], t[:, 3]
    g1x, g1y, g2x, g2y, g3x, g3y = pdetrg(p, t, ar)
    c3 = (g1x .* g2x .+ g1y .* g2y) .* arK
    c1 = (g2x .* g3x .+ g2y .* g3y) .* arK
    c2 = (g3x .* g1x .+ g3y .* g1y) .* arK
    SK = vcat(c3, c1, c2, c3, c1, c2, -c2 - c3, -c3 - c1, -c1 - c2)
    # mass matrix
    aod = arM / 12. # off diagonal element
    ad = 2 * aod    # diagonal element
    SM = vcat(aod, aod, aod, aod, aod, aod, ad, ad, ad)
    # matrices
    I, J = indKM_sparse(t)
    K = sparse(I, J, SK, np, np)
    M = sparse(I, J, SM, np, np)
    return K, M
end
#-------------------------------------------------------------------------------
"""
   g1x, g1y, g2x, g2y, g3x, g3y = pdetrg(p, t, ar)

returns the  the gradient components of the triangle base functions
"""
function pdetrg(p, t, ar)
    # corner point indices
    a1, a2, a3 =  t[:, 1], t[:, 2], t[:, 3]
    # triangle sides
    r23x =  p[a3, 1] .- p[a2, 1]
    r23y =  p[a3, 2] .- p[a2, 2]
    r31x =  p[a1, 1] .- p[a3, 1]
    r31y =  p[a1, 2] .- p[a3, 2]
    r12x =  p[a2, 1] .- p[a1, 1]
    r12y =  p[a2, 2] .- p[a1, 2]
    g1x =  - 0.5 .* r23y ./ ar
    g1y =    0.5 .* r23x ./ ar
    g2x =  - 0.5 .* r31y ./ ar
    g2y =    0.5 .* r31x ./ ar
    g3x =  - 0.5 .* r12y ./ ar
    g3y =    0.5 .* r12x ./ ar
    return g1x, g1y, g2x, g2y, g3x, g3y
end
#-------------------------------------------------------------------------------
"""
    be = btri(T::Array{Int64, 2})

finds boundary edges
"""
function btri(T::Array{Int64, 2})
    # form all faces, non - duplicates are surface triangles
    faces = vcat(T[:, [1, 2]], T[:, [1, 3]], T[:, [2, 3]])
    faces = sort(faces, dims = 2)
    _, ix, jx = uniquerows_int(faces)
    ee = collect(1:(maximum(jx) + 1)) .- 1e-8
    result = fit(Histogram, jx, (1:(maximum(jx) + 1)) .- 1e-8,
                 closed = :right)
    vec = result.weights
    qx = findall(x -> x != 0, vec .- 2)
    return faces[ix[qx], :]
end
#-------------------------------------------------------------------------------
function uniquerows_int(x::Matrix{Int64})
    Vx = [x[k, :] for k = 1:size(x, 1)]
    jj = groupslices(Vx)
    ix, jx = firstinds(jj), similar(jj)
    for k = 1:length(ix)
        jx[jj .== ix[k]] .= k
    end
    Vy = Vx[ix]
    y = [Vy[k][l] for k = 1:length(Vy), l = 1:length(Vy[1])]
    return y, ix, jx
end


text:
@inline apply(o::Sum, diags::Vector{Diagram{W}}) where {W<:Number} = sum(d.weight for d in diags)
@inline apply(o::Prod, diags::Vector{Diagram{W}}) where {W<:Number} = prod(d.weight for d in diags)
@inline apply(o::Sum, diag::Diagram{W}) where {W<:Number} = diag.weight
@inline apply(o::Prod, diag::Diagram{W}) where {W<:Number} = diag.weight

@inline eval(d::DiagramId) = error("eval for $d has not yet implemented!")

function evalDiagNode!(diag::Diagram, varK, varT, evalBare::Function)
    if length(diag.subdiagram) == 0
        if hasproperty(diag.id, :extK)
            K = varK * diag.id.extK
            diag.weight = evalBare(diag.id, K, diag.id.extT, varT) * diag.factor
        else
            diag.weight = evalBare(diag.id, nothing, diag.id.extT, varT) * diag.factor
        end
    else
        diag.weight = apply(diag.operator, diag.subdiagram) * diag.factor
    end
    return diag.weight
end

function evalDiagTree!(diag::Diagram, varK, varT, evalBare::Function)
    for d in PostOrderDFS(diag)
        evalDiagNode!(d, varK, varT, evalBare)
    end
    return diag.weight
end

function evalDiagTree!(diags::Vector{Diagram{W}}, varK, varT, evalBare::Function) where {W}
    for d in diags
        evalDiagTree!(d, varK, varT, evalBare)
    end
    # return W[d.weight for d in diags]
end

function evalDiagTree!(df::DataFrame, varK, varT, evalBare::Function) where {W}
    for d in df[!, :diagram]
        evalDiagTree!(d, varK, varT, evalBare)
    end
    # return W[d.weight for d in df[!, :Diagram]]
end

function evalDiagNode!(diag::Diagram, evalBare::Function, vargs...)
    if length(diag.subdiagram) == 0
        diag.weight = evalBare(diag.id, vargs...) * diag.factor
    else
        diag.weight = apply(diag.operator, diag.subdiagram) * diag.factor
    end
    return diag.weight
end

function evalDiagTree!(diag::Diagram, evalBare::Function, vargs...)
    for d in PostOrderDFS(diag)
        evalDiagNode!(d, evalBare, vargs...)
    end
    return diag.weight
end

function evalDiagTree!(diags::Vector{Diagram{W}}, evalBare::Function, vargs...) where {W}
    for d in diags
        evalDiagTree!(d, evalBare, vargs...)
    end
    return W[d.weight for d in diags]
end

function evalDiagTree!(df::DataFrame, evalBare::Function, vargs...) where {W}
    for d in df[!, :diagram]
        evalDiagTree!(d, evalBare, vargs...)
    end
    # return W[d.weight for d in df[!, :Diagram]]
end


text:
using GRASS
using Test

# run the tests
include("test_geometry.jl")
include("test_physics.jl")
include("test_ccf.jl")
include("test_interpolations.jl")
include("test_input.jl")
include("test_synthesis.jl")

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
module BetaVQE
using Yao, LinearAlgebra
using Flux.Optimise: ADAM, update!
using StatsBase
using VAN

export train, exact

include("hamiltonian.jl")
include("circuit.jl")
include("free_energy.jl")
include("zygote_patch.jl")

# exact value
function exact(β, h)
    H = mat(h)
    w, _ = eigen(Matrix(H))
    Z = sum(e->exp(-β*e), w)
    F = -log(Z)/β
    E = sum(e->exp(-β*e)*e, w)/Z
    E2 = sum(e->exp(-β*e)*e*e, w)/Z
    γ =  sum(e->exp(-2*β*e), w)/Z^2
    S = (E-F)*β
    Cv = (E2-E^2)*β^2
    F, E, S, Cv, γ
end

function loss(β, H, sampler, circuit, nbatch::Int)
    samples = gen_samples(sampler, nbatch)
    free_energy(β, H, sampler, circuit, samples)
end

function train(β::Real, H::AbstractBlock{N}, sampler::AbstractSampler, circuit::AbstractBlock{N}, logfile=nothing;
                optimizer=ADAM(0.1), nbatch::Int=1000, niter::Int=100) where N
    ϕ = model_parameters(sampler)
    θ = parameters(circuit)
    for i = 1:niter
        _, _, gϕ, gθ, _ = gradient(loss, β, H, sampler, circuit, nbatch)
        update!.(Ref(optimizer), ϕ, gϕ)
        update!(optimizer, θ, gθ)
        model_dispatch!(sampler, ϕ)
        dispatch!(circuit, θ)
        message = "$i $(loss(β, H, sampler, circuit, nbatch))"
        println(message)
        flush(stdout)
        if !(logfile === nothing)
            write(logfile, message*"\n")
            flush(logfile)
        end
    end
    ϕ, θ
end

end # module


text:
# Command: Set Timeouts
"""
    timeouts!(session::Session, timeouts::Timeouts)::Nothing
Sets the timeouts.
"""
function timeouts!(session::Session, timeouts::Timeouts)::Nothing
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/timeouts",
        [("Content-Type" => "application/json")],
        JSON3.write(timeouts),
    )
    @assert response.status == 200
    JSON3.read(response.body).value
    nothing
end


text:
using Distributions
using PDMats
using Clustering
using JLD
using Compat

include("compat.jl")
if !isdefined(:GMM)
    include("gmmtypes.jl")
    include("bayestypes.jl")
end
include("gmms.jl")
include("train.jl")
include("io.jl")
include("stats.jl")
include("rand.jl")
include("data.jl")
include("recognizer.jl")

include("bayes.jl")

include("distributions.jl")

## some init code.  Turn off subnormal computation, as it is slow.  This is a global setting...
set_zero_subnormals(true)


text:
function checkaxes(ax_exp, ax_recv)
	ax_exp == ax_recv || throw_axismismatcherror(ax_exp, ax_recv)
end
function throw_axismismatcherror(ax_exp, ax_recv)
	throw(ArgumentError("expected an axis $ax_exp, received $ax_recv"))
end

struct ModeMismatchError{A,B} <: Exception
	mode1 :: A
	mode2 :: B
end

function Base.showerror(io::IO, e::ModeMismatchError{A,A}) where {A}
	print(io,"Modes $(e.mode1) and $(e.mode2) differ")
end
function Base.showerror(io::IO, e::ModeMismatchError{A,B}) where {A,B}
	print(io,"Modes are of different types: $A and $B")
end


text:
# Import all package functions
using Test
using Random
using CpelAsm

# Test set for estimation related functions
@testset "Estimation" begin
    # Initialize all variables
    n=[2,2]; a=[-1.0,1.0]; b=1.0; θ=vcat([a,b]...); ∇logZ=CpelAsm.get_grad_logZ(n,θ);
    ex=CpelAsm.comp_ex(n,a,b); exx=CpelAsm.comp_exx(n,a,b);
    Random.seed!(1234); xobs=[CpelAsm.gen_x_mc(n,a,b) for i=1:50];
    θhat=CpelAsm.est_theta_sa(n,xobs);
    # Check proper computation partition function
    @test CpelAsm.comp_Z(n,a,b)≈235.912 atol=1e-3
    # Check proper computation scaling factor
    @test CpelAsm.comp_g(n,a,b,1.0,1.0)≈143.330 atol=1e-3
    # Check proper likelihood
    @test CpelAsm.comp_lkhd([1,1,0,1,1],n,a,b)≈0.232 atol=1e-3
    # Check proper estimation of θ
    @test sum((θhat-θ).^2)<0.25
end

# Test set for information theory related functions
@testset "Output" begin
    # Initialize all variables
    n=[2,2]; a=[-1.0,1.0]; b=1.0; θ=vcat([a,b]...); ∇logZ=CpelAsm.get_grad_logZ(n,θ);
    ex=CpelAsm.comp_ex(n,a,b); exx=CpelAsm.comp_exx(n,a,b);
    # Check proper E[X] computation
    @test sum(ex)≈0.0 atol=1e-3
    # Check proper E[XX] computation
    @test sum(exx)≈1.285 atol=1e-3
    # Check proper MML computation
    @test comp_mml(trues(4),ex)≈0.5 atol=1e-3
    # Check proper MML computation (2)
    @test CpelAsm.comp_mml_∇(n,∇logZ)≈0.5 atol=1e-3
    # Check proper NME computation
    @test CpelAsm.comp_nme(trues(4),n,a,b,ex,exx)≈0.462 atol=1e-3
    # Check proper NME computation (2)
    h = CpelAsm.comp_nme_∇(n,θ,∇logZ)
    @test h≈0.463 atol=1e-3
    # Check proper NME computation (3)
    @test CpelAsm.comp_nme_mix_exact(trues(4),trues(4),n,n,θ,θ)≈0.463 atol=1e-3
    # Check proper computation of UC
    @test CpelAsm.comp_uc(trues(4),trues(4),n,n,θ,θ,h,h)≈0.0 atol=1e-3
end


text:
################################################################################
# Model
################################################################################

mutable struct WallPendulum{T} <: Model
	n::Int # state dim
	m::Int # control dim
	mp::T # mass
	l::T # length
	g::T # gravity
	k::T # spring
	d::T # distance to wall
end

function dynamics_model(model::WallPendulum{T}, dt::T; mode::Symbol = :none) where {T}
	mp = model.mp
	l = model.l
	g = model.g
	k = model.k
	d = model.d
	B = dt * [0 1 / (mp*l^2)]'
	if mode == :none
		A = I + dt * [0   1;
			          g/l 0]
		c = dt * [0, 0]
	elseif mode == :left
		A = I + dt * [0          1;
			          g/l - k/mp 0]
		c = dt * [0, k*d / (mp*l)]
	elseif mode == :right
		A = I + dt * [0          1;
			          g/l - k/mp 0]
		c = dt * [0, -k*d / (mp*l)]
	else
		error("Unknown contact mode.")
	end
	return A, B, c
end

function get_mode(x::AbstractVector{T}, u::AbstractVector{T}, model::WallPendulum{T}) where {T}
	d = model.d
	l = model.l
	if -d/l <= x[1] <= d/l # none
		return 1
	elseif x[1] > d/l # left
		return 2
	elseif x[1] < -d/l # right
		return 3
	end
end

function dynamics(x0::AbstractVector{T}, u0::AbstractVector{T},
		A::Vector{Matrix{T}}, B::Vector{Matrix{T}}, c::Vector{Vector{T}},
		model::WallPendulum{T}) where {T}
	i = get_mode(x0, u0, model)
	x1 = A[i] * x0 + B[i] * u0 + c[i]
	return x1
end

function kinematics(x::AbstractVector{T}) where {T}
	θ = x[1] + π/2
	θd = x[2]
	p = l * [cos(θ), sin(θ)]
	return p
end

################################################################################
# Domain
################################################################################

"""
	Box domain for the state-control vector (x, u).
	Domain C = {x_min <= x <= x_max} × {u_min <= u <= u_max}.
	C = {S x + R u <= T}
"""
mutable struct BoxDomain{T}
	n::Int # state dim
	m::Int # control dim
	x_min::Vector{T} # box boundary
	x_max::Vector{T} # box boundary
	u_min::Vector{T} # box boundary
	u_max::Vector{T} # box boundary
	S::Matrix{T} # affine inequality representation
	R::Matrix{T} # affine inequality representation
	T::Vector{T} # affine inequality representation
end

function BoxDomain(x_min, x_max, u_min, u_max)
	n = length(x_min)
	m = length(u_min)
	S = Matrix([-I(n); I(n); zeros(m, n); zeros(m, n)])
	R = Matrix([zeros(n, m); zeros(n, m); -I(m); I(m)])
	T = [-x_min; x_max; -u_min; u_max]
	return BoxDomain(n, m, x_min, x_max, u_min, u_max, S, R, T)
end

function inclusion(C::BoxDomain, x, u)
	all(C.S * x + C.R * u .<= C.T)
end

function domain(model::WallPendulum{T}; mode::Symbol = :none) where{T}
	m = model.m
	l = model.l
	g = model.g
	k = model.k
	d = model.d

	u_min = [-4.0]
	u_max = [ 4.0]

	if mode == :none
		# @warn "need to revise when adding right wall."
		x_min = [-d/l, -1.5]
		# x_min = [-d/l*2, -1.5]
		x_max = [ d/l,    1.5]
	elseif mode == :left
		x_min = [ d/l,   -1.5]
		x_max = [ d/l*2,  1.5]
	elseif mode == :right
		x_min = [-d/l*2,   -1.5]
		x_max = [-d/l,  1.5]
	else
		error("Unknown contact mode.")
	end
	return BoxDomain(x_min, x_max, u_min, u_max)
end

# @testset "Domain" begin
# 	# Test
# 	x_min = [-d/l*2, -1.5]
# 	x_max = [ d/l*2,  1.5]
#
# 	x_min_1 = [-d/l*2, -1.5]
# 	x_max_1 = [ d/l,    1.5]
#
# 	x_min_2 = [ d/l,   -1.5]
# 	x_max_2 = [ d/l*2,  1.5]
#
# 	u_min = [-4.0]
# 	u_max = [ 4.0]
#
# 	C0 = BoxDomain(x_min, x_max, u_min, u_max)
# 	C1 = BoxDomain(x_min_1, x_max_1, u_min, u_max)
# 	C2 = BoxDomain(x_min_2, x_max_2, u_min, u_max)
#
# 	xu1 = [[d/l/2,   -1.0], [-2.0]]
# 	xu2 = [[1.5*d/l,  1.0], [ 2.0]]
# 	xu3 = [[2.1d/l,    0.0], [ 0.0]]
#
# 	@test inclusion(C0, xu1...)
# 	@test inclusion(C0, xu2...)
# 	@test !inclusion(C0, xu3...)
# 	@test inclusion(C1, xu1...)
# 	@test !inclusion(C1, xu2...)
# 	@test !inclusion(C1, xu3...)
# 	@test !inclusion(C2, xu1...)
# 	@test inclusion(C2, xu2...)
# 	@test !inclusion(C2, xu3...)
# end

################################################################################
# Problem
################################################################################

mutable struct WallProblem{F}
	model::WallPendulum{F}
	T::Int
	x0::Vector{F}
	Q::F
	Qf::F
	R::F
	β::F
	A::Vector{Matrix{F}}
	B::Vector{Matrix{F}}
	c::Vector{Vector{F}}
	C::Vector{BoxDomain{F}}
end

function build_optimizer(prob::WallProblem{F}) where {F}
	n = prob.model.n
	m = prob.model.m
	nd = length(C)
	T = prob.T
	x0 = prob.x0
	Q = prob.Q
	Qf = prob.Qf
	R = prob.R
	β = prob.β

	Mstar = [β * ones(2n+2m) for i = 1:nd]
	MU = β * ones(n) # upper bound
	ML = - β * ones(n) # lower bound

	optimizer = Model(Gurobi.Optimizer)
	set_silent(optimizer)

	# variables
	@variable(optimizer, x[1:T+1, 1:n])
	@variable(optimizer, u[1:T, 1:m])
	@variable(optimizer, δ[1:nd, 1:T], Bin)
	@variable(optimizer, z[1:nd, 1:T, 1:n])

	# constraints
	# intitial state
	@constraint(optimizer, initial_state, x[1,:] .== x0);
	# 16.a
	@constraint(optimizer, c16a[i in 1:nd, t in 1:T], C[i].S * x[t,:] + C[i].R * u[t,:] - C[i].T .<= Mstar[i] * (1 - δ[i,t]));
	# 16.b
	@constraint(optimizer, c16b[t in 1:T], sum(δ[:, t]) == 1)
	# 18
	@constraint(optimizer, c18[t in 1:T], x[t+1, :] .== sum([z[i,t,:] for i = 1:nd]))
	# 22.a
	@constraint(optimizer, c22a[i in 1:nd, t in 1:T], z[i,t,:] .<= MU * δ[i,t])
	# 22.b
	@constraint(optimizer, c22b[i in 1:nd, t in 1:T], z[i,t,:] .>= ML * δ[i,t])
	# 22.c
	@constraint(optimizer, c22c[i in 1:nd, t in 1:T], z[i,t,:] .<= A[i] * x[t,:] + B[i] * u[t,:] + c[i] - ML * (1 - δ[i,t]))
	# 22.d
	@constraint(optimizer, c22d[i in 1:nd, t in 1:T], z[i,t,:] .>= A[i] * x[t,:] + B[i] * u[t,:] + c[i] - MU * (1 - δ[i,t]))

	# objective
	@objective(optimizer, Min, Q * sum(x[1:T,:].^2) + Qf * sum(x[T+1,:].^2) + R * sum(u.^2))

	return optimizer, x, u, δ, z
end

function get_control(prob::WallProblem{T}) where {T}
	optimizer, x, u, δ, z = build_optimizer(prob)
	optimize!(optimizer)
	return value.(u)[1,:]
end

function simulate!(prob::WallProblem{T}, x0::Vector{T}, H::Int;
		w::Vector{Vector{T}} = Vector{Vector{T}}(), iw::Vector{Int} = Vector{Int}()) where {T}
	x = Vector{Vector{T}}([x0])
	u = Vector{Vector{T}}()
	s = Vector{T}()
	for h = 1:H
		println(h, "/", H)
		prob.x0 = deepcopy(x0)
		s0 = @elapsed begin
			u0 = get_control(prob)
		end
		i = findfirst(x -> x == h, iw)
		(i != nothing) && (u0 += w[i])
		x0 = dynamics(x0, u0, prob.A, prob.B, prob.c, prob.model)
		push!(x, deepcopy(x0))
		push!(u, deepcopy(u0))
		push!(s, s0)
	end
	return x, u, s
end

abstract type ControlPolicy{T}
end

mutable struct OpenLoopPolicy16{T} <: ControlPolicy{T}
	u::Vector{Vector{T}}
	N::Int
	N_sample::Int
	count::Int
	dt::T
end

function OpenLoopPolicy16(u::AbstractVector{V}; N_sample::Int = 1, dt::T = 0.01) where {V,T}
	OpenLoopPolicy16(u, length(u), N_sample, 0, dt)
end

mutable struct Simulator12{T}
	p::ControlPolicy{T}
	x::Vector{Vector{T}}
	u::Vector{Vector{T}}
	N_sample::Int
	dt::T
end

function Simulator12(p::ControlPolicy{T}; dt::T = 0.01) where {T}
	x = Vector{Vector{T}}()
	u = Vector{Vector{T}}()
	Simulator12(p, x, u, p.N_sample, p.dt)
end

function simulate!(sim::Simulator12{T}, x0::Vector{T}, H::Int) where {T}
	N_sample = sim.N_sample
	cnt = 0
	x = deepcopy(x0)
	push!(sim.x, deepcopy(x0))
	for t = 1:H
		x, u = step!(p, x, t)
		push!(sim.x, deepcopy(x))
		push!(sim.u, deepcopy(u))
	end
	return nothing
end

function step!(p::OpenLoopPolicy16{T}, x0::AbstractVector{T}, t::Int) where{T}
	@show t
	@show t%p.N_sample
	((t-1) % p.N_sample) == 0 && (p.count += 1)
	u0 = p.u[p.count]
	x1 = dyna(x0, u0)
	return x1, u0
end


text:
using StochasticDiffEq, Test
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear
prob = prob_sde_2Dlinear
choice_function(integrator) = (Int(integrator.t<0.5) + 1)
alg_double = StochasticCompositeAlgorithm((SRIW1(),SRIW1()),choice_function)
alg_double2 = StochasticCompositeAlgorithm((SRI(),SRI()),choice_function)
alg_switch = StochasticCompositeAlgorithm((EM(),RKMil()),choice_function)

Random.seed!(100)
@time sol1 = solve(prob_sde_linear,alg_double)
Random.seed!(100)
@time sol2 = solve(prob_sde_linear,SRIW1())
@test sol1.t == sol2.t
@test sol1(0.8) == sol2(0.8)

Random.seed!(100)
integrator1 = init(prob,alg_double2)
Random.seed!(100)
integrator2 = init(prob,SRI())
Random.seed!(100)
solve!(integrator1)
Random.seed!(100)
solve!(integrator2)

@test integrator1.sol.t == integrator2.sol.t

sol = solve(prob,alg_switch,dt=1/8)


text:
import .Compose

# This is necessary as without `false` compose tries to draw a png by themself when calling show
# instead of writing it to `io`
show_custom(io::IO, ::MIME"image/png", ctx::Compose.Context) =
    Compose.draw(Compose.PNG(io, Compose.default_graphic_width, Compose.default_graphic_height, false), ctx)


text:
@testset "E" begin
    @testset "edlimb" begin
        viewpt = [2.0, 0.0, 0.0]
        limb = edlimb(sqrt(2), 2.0 * sqrt(2), sqrt(2), viewpt)
        expected_sminor = [0.0, 0.0, -1.0]
        expected_smajor = [0.0, 2.0, 0.0]
        expected_center = [1.0, 0.0, 0.0]
        @test center(limb) ≈ expected_center
        @test semi_major(limb) ≈ expected_smajor
        @test semi_minor(limb) ≈ expected_sminor
        @test_throws ArgumentError edlimb(sqrt(2), 2.0 * sqrt(2), sqrt(2), viewpt[1:2])
        @test_throws SpiceError edlimb(sqrt(2), 2.0 * sqrt(2), -sqrt(2), viewpt)
    end
    @testset "edterm" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            et = str2et("2007 FEB 3 00:00:00.000")
            trgepc, obspos, trmpts = edterm("UMBRAL", "SUN", "MOON", et, "IAU_MOON", "LT+S", "EARTH", 3)
            expected_trgepc = 223732863.86351674795
            expected_obspos = [394721.1024056578753516078,
                               27265.11780063395417528227,
                               -19069.08478859506431035697]
            expected_trmpts0 = [-1.53978381936825627463e+02,
                                -1.73056331949840728157e+03,
                                1.22893325627419600088e-01]
            expected_trmpts1 = [ 87.37506200891714058798,
                                864.40670594653545322217,
                                1504.56817899807947469526]
            expected_trmpts2 = [ 42.21324376177891224415,
                                868.21134635239388899208,
                                -1504.3223923468244720425]
            expected_trmpts = [ -1.53978381936825627463e+02,
                               -1.73056331949840728157e+03,
                               1.22893325627419600088e-01,
                               87.37506200891714058798,
                               864.40670594653545322217,
                               1504.56817899807947469526,
                               42.21324376177891224415,
                               868.21134635239388899208,
                               -1504.3223923468244720425]

            @test trgepc ≈ expected_trgepc
            @test obspos[1] ≈ expected_obspos[1]
            @test obspos[2] ≈ expected_obspos[2]
            @test obspos[3] ≈ expected_obspos[3]

            @testset for i in eachindex(trmpts, expected_trmpts)
            @test obspos[1] ≈ expected_obspos[1]
                @test trmpts[i] ≈ expected_trmpts[i]
            end
            iluet0, srfvec0, phase0, solar0, emissn0 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,1], abcorr="LT+S")
            @test rad2deg(solar0) ≈ 90.269765819
            iluet1, srfvec1, phase1, solar1, emissn1 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,2], abcorr="LT+S")
            @test rad2deg(solar1) ≈ 90.269765706
            iluet2, srfvec2, phase2, solar2, emissn2 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,3], abcorr="LT+S")
            @test rad2deg(solar2) ≈ 90.269765730
            # penumbral
            trgepc, obspos, trmpts = edterm("PENUMBRAL", "SUN", "MOON", et, "IAU_MOON", "LT+S", "EARTH", 3)
            expected_trmpts = [1.54019056755619715204e+02,
                               1.73055969989532059117e+03,
                               -1.23508409498995316844e-01,
                               -87.33436047798454637814,
                               -864.41003834758112134296,
                               -1504.56862757530461749411,
                               -42.17254722919552278881,
                               -868.21467833235510624945,
                               1504.32161075630597224517]
            @test trgepc ≈ expected_trgepc
            @test obspos[1] ≈ expected_obspos[1]
            @test obspos[2] ≈ expected_obspos[2]
            @test obspos[3] ≈ expected_obspos[3]
            @testset for i in eachindex(trmpts, expected_trmpts)
                @test trmpts[i] ≈ expected_trmpts[i]
            end
            iluet0, srfvec0, phase0, solar0, emissn0 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,1], abcorr="LT+S")
            @test rad2deg(solar0) ≈ 89.730234406
            iluet1, srfvec1, phase1, solar1, emissn1 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,2], abcorr="LT+S")
            @test rad2deg(solar1) ≈ 89.730234298
            iluet2, srfvec2, phase2, solar2, emissn2 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,3], abcorr="LT+S")
            @test rad2deg(solar2) ≈ 89.730234322
        finally
            kclear()
        end
    end
    @testset "ek" begin
        try
            ekpath = tempname()
            tablename = "TEST_TABLE_EK"
            handle = ekopn(ekpath, ekpath, 0)
            decls = ["DATATYPE = CHARACTER*(10),   NULLS_OK = TRUE, SIZE = VARIABLE",
                     "DATATYPE = DOUBLE PRECISION, NULLS_OK = TRUE, SIZE = VARIABLE",
                     "DATATYPE = INTEGER,          NULLS_OK = TRUE, SIZE = VARIABLE"]
            segno = ekbseg(handle, tablename, ["c1", "d1", "i1"], decls)
            c_data = [["100"], ["101", "101"], ["102", "102", "102"]]
            d_data = [[100.0], [101.0, 101.0], [102.0, 102.0, 102.0]]
            i_data = [[100], [101, 101], [102, 102, 102]]
            for r in 1:3
                ekinsr(handle, segno, r)
                ekacec(handle, segno, r, "c1", c_data[r], false)
                ekaced(handle, segno, r, "d1", d_data[r], false)
                ekacei(handle, segno, r, "i1", i_data[r], false)
            end
            ekinsr(handle, segno, 4)
            ekacec(handle, segno, 4, "c1", c_data[1], true)
            ekaced(handle, segno, 4, "d1", d_data[1], true)
            ekacei(handle, segno, 4, "i1", i_data[1], true)
            @test_throws SpiceError ekinsr(handle, segno, 6)
            ekcls(handle)
            kclear()

            handle = eklef(ekpath)
            query = "SELECT c1, d1, i1 from $tablename"
            xbegs, xends, xtypes, xclass, tabs, cols = ekpsel(query)
            @test xtypes == collect(0:2)
            @test xclass == zeros(Int, 3)
            @test tabs == fill("TEST_TABLE_EK", 3)
            @test cols == ["C1", "D1", "I1"]
            nmrows = ekfind(query)
            @test nmrows == 4
            for r in 1:nmrows - 1
                n_elm = eknelt(1, r)
                @test n_elm == r
                for e in 1:n_elm
                    i_datum = ekgi(3, r, e)
                    @test i_datum == i_data[r][e]
                    d_datum = ekgd(2, r, e)
                    @test d_datum == d_data[r][e]
                    c_datum = ekgc(1, r, e)
                    @test c_datum == c_data[r][e]
                end
            end
            for r in 1:nmrows - 1
                ri_data = ekrcei(handle, segno, r, "i1")
                @test ri_data == i_data[r]
                rd_data = ekrced(handle, segno, r, "d1")
                @test rd_data == d_data[r]
                rc_data = ekrcec(handle, segno, r, "c1")
                @test rc_data == c_data[r]

                @test_throws SpiceError ekrcei(handle, segno, 5, "i1")
                @test_throws SpiceError ekrced(handle, segno, 5, "d1")
                # This segfaults on some platforms
                @test_skip ekrcec(handle, segno, 5, "c1")
            end
            @test ekrcei(handle, segno, 4, "i1") === missing
            @test ekrced(handle, segno, 4, "d1") === missing
            @test ekrcec(handle, segno, 4, "c1") === missing

            ekuef(handle)
            handle = ekopw(ekpath)
            c_data = [["200"], ["201", "201"], ["202", "202", "202"]]
            d_data = [[200.0], [201.0, 201.0], [202.0, 202.0, 202.0]]
            i_data = [[200], [201, 201], [202, 202, 202]]
            for r in 1:3
                ekucec(handle, segno, r, "c1", c_data[r], false)
                ekuced(handle, segno, r, "d1", d_data[r], false)
                ekucei(handle, segno, r, "i1", i_data[r], false)
            end
            ekcls(handle)
            handle = ekopr(ekpath)
            for r in 1:nmrows - 1
                ri_data = ekrcei(handle, segno, r, "i1")
                @test ri_data == i_data[r]
                rd_data = ekrced(handle, segno, r, "d1")
                @test rd_data == d_data[r]
                rc_data = ekrcec(handle, segno, r, "c1")
                @test rc_data == c_data[r]
            end
            ekuef(handle)
        finally
            kclear()
        end
    end
    @testset "ekntab/ektnam/ekccnt" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno = ekbseg(handle, "TEST_TABLE_EKCCNT", ["c1"], ["DATATYPE  = INTEGER, NULLS_OK = TRUE"])
            recno = ekappr(handle, segno)
            ekacei(handle, segno, recno, "c1", [1, 2], false)
            ekcls(handle)
            kclear()
            furnsh(ekpath)
            @test ekntab() == 1
            @test ektnam(1) == "TEST_TABLE_EKCCNT"
            @test ekccnt("TEST_TABLE_EKCCNT") == 1
        finally
            kclear()
        end
    end
    @testset "ekcii" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno = ekbseg(handle, "TEST_TABLE_EKCII", ["c1"], ["DATATYPE = INTEGER, NULLS_OK = TRUE"])
            recno = ekappr(handle, segno)
            ekacei(handle, segno, recno, "c1", [1, 2], false)
            ekcls(handle)
            kclear()
            furnsh(ekpath)
            @test ekntab() == 1
            @test ektnam(1) == "TEST_TABLE_EKCII"
            @test ekccnt("TEST_TABLE_EKCII") == 1
            column, attdsc = ekcii("TEST_TABLE_EKCII", 1)
            @test column == "C1"
            @test attdsc.cclass == 1
            @test attdsc.dtype == 2
            @test attdsc.size == 1
            @test attdsc.strlen == 1
            @test attdsc.indexd == 0
            @test attdsc.nullok == 1
        finally
            kclear()
        end
    end
    @testset "ekdelr" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno, rcptrs = ekifld(handle, "TEST_TABLE_EKDELR", 2, ["c1"],
                                   ["DATATYPE = INTEGER, NULLS_OK = TRUE"])
            ekacli(handle, segno, "c1", [[1], [2]], [false, false], rcptrs)
            ekffld(handle, segno, rcptrs)
            ekdelr(handle, segno, 1)
            ekcls(handle)
        finally
            kclear()
        end
    end
    @testset "ekfind" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno, rcptrs = ekifld(handle, "test_table_ekfind", 2, ["cc1"],
                                   ["DATATYPE = INTEGER, NULLS_OK = TRUE"])
            ekacli(handle, segno, "cc1", [[1], [2]], [false, false], rcptrs)
            ekffld(handle, segno, rcptrs)
            ekcls(handle)
            kclear()
            furnsh(ekpath)
            nmrows = ekfind("SELECT CC1 FROM TEST_TABLE_EKFIND WHERE CC1 > 0", 100)
            @test nmrows == 2
        finally
            kclear()
        end
    end
    @testset "ekclc/ekgc" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno, rcptrs = ekifld(handle, "TEST_TABLE_EKGC", 3, ["c1"],
                                   ["DATATYPE = CHARACTER*(*), INDEXED = TRUE, NULLS_OK = TRUE"])
            ekaclc(handle, segno, "c1", ["1.0", "2.0", ""], [false, false, true], rcptrs)
            ekffld(handle, segno, rcptrs)
            ekcls(handle)
            kclear()
            furnsh(ekpath)
            nmrows = ekfind("SELECT C1 FROM TEST_TABLE_EKGC", 100)
            @test nmrows == 3
            c = ekgc(1, 1, 1)
            @test c == "1.0"
            c = ekgc(1, 2, 1)
            @test c == "2.0"
            c = ekgc(1, 3, 1)
            @test c === missing
        finally
            kclear()
        end
    end
    @testset "ekcld/ekgd" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno, rcptrs = ekifld(handle, "TEST_TABLE_EKGD", 3, ["c1"],
                                   ["DATATYPE = DOUBLE PRECISION, NULLS_OK = TRUE"])
            ekacld(handle, segno, "c1", [[1.0], [2.0], [1.0]], [false, false, true], rcptrs)
            ekffld(handle, segno, rcptrs)
            ekcls(handle)
            kclear()
            furnsh(ekpath)
            nmrows = ekfind("SELECT C1 FROM TEST_TABLE_EKGD", 100)
            @test nmrows == 3
            d = ekgd(1, 1, 1)
            @test d == 1.0
            d = ekgd(1, 2, 1)
            @test d == 2.0
            d = ekgd(1, 3, 1)
            @test d === missing
        finally
            kclear()
        end
    end
    @testset "ekcli/ekgi" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno, rcptrs = ekifld(handle, "TEST_TABLE_EKGI", 3, ["c1"],
                                   ["DATATYPE = INTEGER, NULLS_OK = TRUE"])
            ekacli(handle, segno, "c1", [[1, 2], [2, 3], [1, 1]], [false, false, true], rcptrs)
            ekffld(handle, segno, rcptrs)
            ekcls(handle)
            kclear()
            furnsh(ekpath)
            nmrows = ekfind("SELECT C1 FROM TEST_TABLE_EKGI", 100)
            @test nmrows == 3
            d = ekgi(1, 1, 1)
            @test d == 1
            d = ekgi(1, 2, 1)
            @test d == 2
            d = ekgi(1, 3, 1)
            @test d === missing
        finally
            kclear()
        end
    end
    @testset "ekinsr" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno = ekbseg(handle, "TEST_TABLE_EKINSR", ["c1"],
                           ["DATATYPE = INTEGER, NULLS_OK = TRUE"])
            ekinsr(handle, segno, 1)
            ekacei(handle, segno, 1, "c1", [1, 2], false)
            ekcls(handle)
            furnsh(ekpath)
            nmrows = ekfind("SELECT C1 FROM TEST_TABLE_EKINSR", 100)
            @test nmrows == 1
        finally
            kclear()
        end
    end
    @testset "eklef/eknelt/ekuef" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno = ekbseg(handle, "TEST_TABLE_EKLEF", ["C1"],
                           ["DATATYPE  = INTEGER, NULLS_OK = TRUE"])
            recno = ekappr(handle, segno)
            ekacei(handle, segno, recno, "c1", [1, 2], false)
            ekcls(handle)
            handle = eklef(ekpath)
            nmrows = ekfind("SELECT C1 FROM TEST_TABLE_EKLEF", 100)
            @test eknelt(1, 1) == 1
            @test nmrows == 1
            ekuef(handle)
        finally
            kclear()
        end
    end
    @testset "eknseg" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno = ekbseg(handle, "TEST_TABLE_EKNSEG", ["C1"],
                           ["DATATYPE  = INTEGER, NULLS_OK = TRUE"])
            recno = ekappr(handle, segno)
            ekacei(handle, segno, recno, "c1", [1, 2], false)
            ekcls(handle)
            kclear()
            handle = ekopr(ekpath)
            @test eknseg(handle) == 1
            ekcls(handle)
        finally
            kclear()
        end
    end
    @testset "ekops" begin
        try
            handle = ekops()
            @test handle != -1
            ekcls(handle)
        finally
            kclear()
        end
    end
    @testset "ekssum" begin
        try
            ekpath = tempname()
            handle = ekopn(ekpath, ekpath, 0)
            segno, rcptrs = ekifld(handle, "test_table_ekssum", 2, ["c1"],
                                   ["DATATYPE = INTEGER, NULLS_OK = TRUE"])
            ekacli(handle, segno, "c1", [[1], [2]], [false, false], rcptrs)
            ekffld(handle, segno, rcptrs)
            segsum = ekssum(handle, segno)
            @test segsum.ncols == 1
            @test segsum.nrows == 2
            @test SPICE.cnames(segsum) == ["C1"]
            @test SPICE.tabnam(segsum) == "TEST_TABLE_EKSSUM"
            c1descr = segsum.cdescrs[1]
            @test c1descr.dtype == 2
            @test !Bool(c1descr.indexd)
            @test Bool(c1descr.nullok)
            ekcls(handle)
        finally
            kclear()
        end
    end
    @testset "el2cgv" begin
        vec1 = [1.0, 1.0, 1.0]
        vec2 = [1.0, -1.0, 1.0]
        center = [1.0, 1.0, 1.0]
        smajor, sminor = saelgv(vec1, vec2)
        ellipse = cgv2el(center, smajor, sminor)
        out_center, out_smajor, out_sminor = el2cgv(ellipse)
        expected_center = [1.0, 1.0, 1.0]
        expected_smajor = [sqrt(2.0), 0.0, sqrt(2.0)]
        expected_sminor = [0.0, sqrt(2.0), 0.0]
        @test out_center ≈ expected_center
        @test out_smajor ≈ expected_smajor
        @test out_sminor ≈ expected_sminor
    end
    @testset "eqncpv" begin
        p = 10000.0
        gm = 398600.436
        ecc = 0.1
        a = p / (1.0 - ecc)
        n = sqrt(gm / a) / a
        argp = deg2rad(30.0)
        node = deg2rad(15.0)
        inc = deg2rad(10.0)
        m0 = deg2rad(45.0)
        t0 = -100000000.0
        eqel = [a, ecc * sin(argp + node), ecc * cos(argp + node), m0 + argp + node,
                tan(inc / 2.0) * sin(node), tan(inc / 2.0) * cos(node), 0.0, n, 0.0]
        state = eqncpv(t0 - 9750.0, t0, eqel,  -π/2, π/2)
        expected = [-10732.167433285387, 3902.505790600528, 1154.4516152766892,
                    -2.540766899262123, -5.15226920298345, -0.7615758062877463]
        @test expected ≈ state rtol=1e-6
    end
    @testset "eqstr" begin
        @test eqstr("A short string    ", "ashortstring")
        @test eqstr("Embedded        blanks", "Em be dd ed bl an ks")
        @test !eqstr("One word left out", "WORD LEFT OUT")
    end
    @testset "esrchc" begin
        array = ["This", "is", "a", "test"]
        @test esrchc("This", array) == 1
        @test esrchc("is", array) == 2
        @test esrchc("a", array) == 3
        @test esrchc("test", array) == 4
        @test esrchc("fail", array) == -1
    end
    @testset "et2lst" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            et = str2et("2004 may 17 16:30:00")
            hr, mn, sc, time, ampm = et2lst(et, 399, deg2rad(281.49521300000004), "planetocentric", 51, 51)
            @test hr == 11
            @test mn == 19
            @test sc == 22
            @test time == "11:19:22"
            @test ampm == "11:19:22 A.M."
        finally
            kclear()
        end
    end
    @testset "et2utc" begin
        try
            furnsh(path(CORE, :lsk))
            et = -527644192.5403653
            output = et2utc(et, :J, 6)
            @test output == "JD 2445438.006415"
            @test_throws SpiceError et2utc(et, :NORBERT, 6)
        finally
            kclear()
        end
    end
    @testset "etcal" begin
        et = 0.0
        cal = etcal(et)
        @test cal == "2000 JAN 01 12:00:00.000"
    end
    @testset "eul2m" begin
        act = eul2m(3, 2, 1, 3, 2, 1)
        @test size(act) == (3, 3)
        exp = [0.411982245665683 -0.6812427202564033 0.6051272472413688;
               0.05872664492762098 -0.642872836134547 -0.7637183366502791;
               0.9092974268256817 0.35017548837401463  -0.2248450953661529]
        @testset for i in eachindex(act, exp)
            @test act[i] ≈ exp[i]
        end
    end
    @testset "eul2xf" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck))
            et = str2et("Jan 1, 2009")
            expected = sxform("IAU_EARTH", "J2000", et)
            eul = [1.571803284049681, 0.0008750002978301174, 2.9555269829740034,
                   3.5458495690569166e-12, 3.080552365717176e-12, -7.292115373266558e-05]
            out = eul2xf(eul, 3, 1, 3)
            @test out ≈ expected
        finally
            kclear()
        end
    end
    @testset "expool" begin
        try
            textbuf = ["DELTET/K = 1.657D-3", "DELTET/EB = 1.671D-2"]
            lmpool(textbuf)
            @test expool("DELTET/K")
            @test expool("DELTET/EB")
        finally
            kclear()
        end
    end
end


text:
using .Contexts
const CallbackContextItems = Union{Nothing, Vector{NamedTuple}}
const TriggeredParam = NamedTuple{(:prop_id, :value)}
mutable struct CallbackContext
    response::HTTP.Response
    inputs::Dict{String, Any}
    states::Dict{String, Any}
    outputs_list::Vector{Any}
    inputs_list::Vector{Any}
    states_list::Vector{Any}
    triggered::Vector{TriggeredParam}
    function CallbackContext(response, outputs_list, inputs_list, states_list, changed_props)
        input_values = inputs_list_to_dict(inputs_list)
        state_values = inputs_list_to_dict(states_list)
        triggered = TriggeredParam[(prop_id = id, value = input_values[id]) for id in changed_props]
        return new(response, input_values, state_values, outputs_list, inputs_list, states_list, triggered)
    end
end

const _callback_context_storage = TaskContextStorage()

function with_callback_context(f, context::CallbackContext)
    return with_context(f, _callback_context_storage, context)
end

"""
    callback_context()::CallbackContext

    Get context of current callback, available only inside callback processing function
"""
function callback_context()
    !has_context(_callback_context_storage) && error("callback_context() is only available from a callback processing function")
    return get_context(_callback_context_storage)
end

function inputs_list_to_dict(list::Vector{Any})
    result = Dict{String, Any}()
    _item_to_dict!.(Ref(result), list)
    return result
end

dep_id_string(id::NamedTuple) = sorted_json(id)
dep_id_string(id::AbstractString) = String(id)
function _item_to_dict!(target::Dict{String, Any}, item)
    target["$(dep_id_string(item.id)).$(item.property)"] = get(item, :value, nothing)
end

_item_to_dict!(target::Dict{String, Any}, item::Vector) =  _item_to_dict!.(Ref(target), item)


text:
immutable PersistentSet{T}
    dict::PersistentHashMap{T, (@compat Void)}

    PersistentSet(d::PersistentHashMap{T, (@compat Void)}) = new(d)
    PersistentSet() = new(PersistentHashMap{T, (@compat Void)}())
    PersistentSet(itr) = _union(new(PersistentHashMap{T, (@compat Void)}()), itr)
end
PersistentSet() = PersistentSet{Any}()
PersistentSet(itr) = PersistentSet{eltype(itr)}(itr)

# TODO: Depricate this invocation
PersistentSet{T}(x1::T, x2::T, xs::T...) =
    PersistentSet{T}(PersistentHashMap((x1, nothing), (x2, nothing),
                                        [(x, nothing) for x in xs]...))

Base.hash(s::PersistentSet,h::UInt) =
    hash(s.dict, h+(@compat UInt(0xf7dca1a5fd7090be)))

Base.conj{T}(s::PersistentSet{T}, val) =
    PersistentSet{T}(assoc(s.dict, val, nothing))

push(s::PersistentSet, val) = conj(s, val)

disj{T}(s::PersistentSet{T}, val) =
    PersistentSet{T}(dissoc(s.dict, val))

Base.length(s::PersistentSet) = length(s.dict)
Base.eltype{T}(s::PersistentSet{T}) = T

Base.isequal(s1::PersistentSet, s2::PersistentSet) = isequal(s1.dict, s2.dict)
==(s1::PersistentSet, s2::PersistentSet) = s1.dict == s2.dict

Base.in(val, s::PersistentSet) = haskey(s.dict, val)

Base.start(s::PersistentSet) = start(s.dict)
Base.done(s::PersistentSet, state) = done(s.dict, state)
function Base.next(s::PersistentSet, state)
    kv, state = next(s.dict, state)
    (kv[1], state)
end

function Base.filter{T}(f::Function, s::PersistentSet{T})
    filtered = filter((kv) -> f(kv[1]), s.dict)
    PersistentSet{T}(keys(filtered))
end

function Base.setdiff(l::PersistentSet, r::(@compat Union{PersistentSet, Set}))
    notinr(el) = !(el in r)
    filter(notinr, l)
end

import Base.-
-(l::PersistentSet, r::(@compat Union{PersistentSet, Set})) = setdiff(l, r)

Base.isempty(s::PersistentSet) = length(s.dict) == 0

function _union(s::PersistentSet, xs)
    for x in xs
        s = conj(s, x)
    end
    s
end

join_eltype() = @compat Union{}
join_eltype(v1, vs...) = typejoin(eltype(v1), join_eltype(vs...))

Base.union(s::PersistentSet) = s
Base.union(x::PersistentSet, xs::PersistentSet...) =
    reduce(_union, PersistentSet{join_eltype(x, xs...)}(), [x, xs...])

function Base.isless(s1::PersistentSet, s2::PersistentSet)
    length(s1) < length(s2) &&
    all(el -> el in s2, s1)
end

import Base.<=
<=(s1::PersistentSet, s2::PersistentSet) = all(el -> el in s2, s1)

function Base.show{T}(io::IO, s::PersistentSet{T})
    print(io, "PersistentSet{$T}(")
    print(io, join([k for (k, v) in s.dict], ", "))
    print(io, ")")
end


text:
@nospecialize
export TagfulPattern,
    And,
    Or,
    Literal,
    Wildcard,
    Deconstrucution,
    Guard,
    Effect,
    untagless,
    TagfulPattern,
    PatternInfo

abstract type TagfulPattern end

struct PatternInfo
    pattern::TagfulPattern
    typetag::TypeObject
end

struct And <: TagfulPattern
    ps::Vector{PatternInfo}
end

struct Or <: TagfulPattern
    ps::Vector{PatternInfo}
end

struct Literal{T} <: TagfulPattern
    val::T
end

struct Wildcard <: TagfulPattern end

struct Deconstrucution <: TagfulPattern
    comp::PComp
    extract::Function
    params::Vector{PatternInfo}
end

struct Guard <: TagfulPattern
    predicate::Any
end

struct Effect <: TagfulPattern
    perform::Any
end

@specialize
function _uncurry_call_argtail(f)
    function (_, args...)
        f(args...)
    end
end
@nospecialize

function untagless(points_of_view::Dict{Function, Int})
    myviewpoint = points_of_view[untagless]
    typetag_viewpoint::Int = points_of_view[tag_extract]
    mk_info(all_info)::PatternInfo =
        PatternInfo(all_info[myviewpoint], all_info[typetag_viewpoint])
    ! = mk_info
    function decons(comp::PComp, extract::Function, ps)
        Deconstrucution(comp, extract, PatternInfo[!p for p in ps])
    end
    and(ps::Vector{Vector{Any}}) = And(PatternInfo[!e for e in ps])
    or(ps::Vector{Vector{Any}}) = Or(PatternInfo[!e for e in ps])

    (
        and = and,
        or = or,
        literal = Literal,
        wildcard = Wildcard(),
        decons = decons,
        guard = Guard,
        effect = Effect,
    )
end
@specialize


text:
module MCMCChain

import Showoff: showoff
import StatsBase: autocor, autocov, countmap, counts, describe, predict,
       quantile, sample, sem, summarystats
import LinearAlgebra: diag

using RecipesBase
import RecipesBase: plot

using Distributions
using SpecialFunctions

export Chains, getindex, setindex!
export plot, traceplot, meanplot, densityplot, histogramplot, mixeddensityplot, autcorplot
export describe

# export diagnostics functions
export discretediag, gelmandiag, gewekediag, heideldiag, rafterydiag

abstract type AbstractChains end

"""
    Chains type

Parameters:

- `value`: `iterations × variables × chains` Data array
- `range`: Range describing the iterations (considering thinning)
- `names`: List of variable names (strings)
- `chains`: List of chain ids
"""
struct Chains{T<:Real} <: AbstractChains
    logevidence::T
    value::Array{Union{Missing, T}, 3}
    range::AbstractRange{Int}
    names::Vector
    uniquenames::Dict{Symbol, Int}
    chains::Vector{Int}
end

# imports
include("utils.jl")

include("chains.jl")
include("chainsummary.jl")
include("discretediag.jl")
include("fileio.jl")
include("gelmandiag.jl")
include("gewekediag.jl")
include("heideldiag.jl")
include("mcse.jl")
#include("modelchains.jl")
#include("modelstats.jl")
include("rafterydiag.jl")
include("stats.jl")
include("plot.jl")
#include("plot2.jl")

end # module


text:
struct LazyModel{Name,T} <: AbstractMillModel
    m::T
end

const LazyModel{Name} = LazyModel{Name, T} where {T}
LazyModel{Name}(m::M) where {Name, M} = LazyModel{Name,M}(m)

Flux.@functor LazyModel

function (m::LazyModel{Name})(x::LazyNode{Name}) where {Name}
    ds = unpack2mill(x)
    m.m(ds)
end

function _reflectinmodel(ds::LazyNode{Name}, db, da, b, a, s) where {Name}
    pm, d = Mill._reflectinmodel(unpack2mill(ds), db, da, b, a, s * Mill.encode(1, 1))
    LazyModel{Name}(pm), d
end

function Mill.HiddenLayerModel(m::LazyModel{N}, ds::LazyNode{N}, n) where {N}
    hm, o = Mill.HiddenLayerModel(m.m, unpack2mill(ds), n)
    return(LazyModel{N}(hm), o )
end

function Mill.mapactivations(hm::LazyModel{N}, x::LazyNode{N}, m::LazyModel{N}) where {N}
    ho, o = Mill.mapactivations(hm.m, unpack2mill(x), m.m)
end

function unpack2mill end


text:
# Our own symbolic graph

using Debugger

# include("utils.jl")
# using Main.utils
include("aimajulia.jl");
# include("planning.jl");
using Main.aimajulia;
# using Main.planning;

function define()
    knowledge_base = [
        expr("Door(D1)"),
        expr("Door(D2)"),
        expr("Spot(S1)"),
        expr("Spot(S2)"),
        expr("Cabinet(C1)"),
        expr("Cabinet(C2)"),
        expr("CabinetDoor(D1, C1)"),
        expr("CabinetDoor(D2, C2)"),
        expr("Arm(A1)"),

        expr("State(S1, Empty)"),
        expr("State(S2, Empty)"),
        expr("State(C1, Occupied)"),
        expr("State(C2, Occupied)"),
        expr("State(A1, Empty)"),
        expr("State(D1, Closed)"),
        expr("State(D2, Closed)"),
    ]

    # Open door action
    precond_pos = [
        expr("State(d, Closed)"),
        expr("State(a, Empty)"),
        expr("Arm(a)"),
        expr("Door(d)"),
    ]
    precond_neg = []
    effect_add = [expr("State(d, Open)")]
    effect_rem = [expr("State(d, Closed)")] # arm is still empty after the action
    open_door = PlanningAction(
        expr("OpenDoor(d, a)"),
        (precond_pos, precond_neg),
        (effect_add, effect_rem),
    )

    # Pick from cabinet action (Pick from cabinet c that's occupied, it's door d is open, arm a should be empty)
    precond_pos = [
        expr("Cabinet(c)"),
        expr("Door(d)"),
        expr("CabinetDoor(d,c)"),
        expr("State(d, Open)"),
        expr("State(c, Occupied)"),
        expr("Arm(a)"),
        expr("State(a, Empty)"),
    ]
    precond_neg = []
    effect_add = [expr("State(c, Empty)"), expr("State(a, Occupied)")]
    effect_rem = [expr("State(c, Occupied)"), expr("State(a, Empty)")] # arm is still empty after the action
    pick_from_cabinet = PlanningAction(
        expr("PickFromCabinet(c, d, a)"),
        (precond_pos, precond_neg),
        (effect_add, effect_rem),
    )

    # Place into spot s (Arm a should be occupied, spot s should be empty)
    precond_pos = [
        expr("Spot(s)"),
        expr("State(s, Empty)"),
        expr("Arm(a)"),
        expr("State(a, Occupied)"),
    ]
    precond_neg = []
    effect_add = [expr("State(s, Occupied)"), expr("State(a, Empty)")]
    effect_rem = [expr("State(s, Empty)"), expr("State(a, Occupied)")]
    place_on_spot = PlanningAction(
        expr("PlaceOnSpot(s,a)"),
        (precond_pos, precond_neg),
        (effect_add, effect_rem),
    )

    # function goal_test(kb::PDDL)
    #     return ask(kb, expr("State(D1, Open)"));
    # end

    # prob = PDDL(knowledge_base, [open_door, pick_from_cabinet, place_on_spot], goal_test)
    prob = PDDL(
        knowledge_base,
        [open_door, pick_from_cabinet, place_on_spot],
        goal_test,
    )

    fol_kb = FirstOrderLogicKnowledgeBase(knowledge_base)
    gpp = GraphPlanProblem(prob, fol_kb)

    # @bp
    graphplan(gpp, ([expr("State(D1, Open)"), expr("State(C1, Empty)")], []))
    println(fol_kb)
    println("STARTING STATE")

    execute_action(open_door, fol_kb, (expr("D1"),expr("A1")))
    println(fol_kb)
    println("EXECUTED DOOR OPEN")

    execute_action(pick_from_cabinet, fol_kb, (expr("C1"),expr("D1"),expr("A1")))
    println(fol_kb)
    println("EXECUTED PICK")

    execute_action(open_door, fol_kb, (expr("D2"),expr("A1")))
    println(fol_kb)
    println("EXECUTED DOOR OPEN")

    for soln in gpp.solution
        println(soln)
    end
    println("Finito!!")
end

# @enter define()

max_stars_repo_name: uscresl/aima-julia

max_issues_repo_name: uscresl/aima-julia

max_forks_repo_name: uscresl/aima-julia


text:
# test 4: manual help/version, import another parser

using ArgParse
using Base.Test

function ap_test4(args)

    s0 = ArgParseSettings()  # a "parent" structure e.g. one with some useful set of rules
                             # which we want to extend

    # So we just add a simple table
    @add_arg_table s0 begin
        "--parent-flag", "-o"
            action = "store_true"
            help = "parent flag"
        "--flag"
            action = "store_true"
            help = "another parent flag"
        "parent-argument"
            help = "parent argument"
    end

    s = ArgParseSettings("Test 4 for ArgParse.jl",
                         add_help = false,           # disable auto-add of --help option
                         version = "Version 1.0",    # we set the version info, but --version won't be added
                         error_on_conflict = false)  # do not error-out when trying to override an option

    import_settings(s, s0)       # now s has all of s0 arguments (except help/version)

    @add_arg_table s begin
        "-o"                       # this will partially override s0's --parent-flag
            action = :store_true
            help = "child flag"
        "--flag"                   # this will fully override s0's --flag
            action = :store_true
            help = "another child flag"
        "-?", "--HELP", "--¡ḧëļṕ"                # (almost) all characters allowed
            action = :show_help                  # will invoke the help generator
            help = "this will help you"
        "-v", "--VERSION"
            action = :show_version               # will show version information
            help = "show version information" *
                   "and exit"
    end

    s.exc_handler = (settings, err)->error(err.text)

    parsed_args = parse_args(args, s)
end

function ap_test4_fails(args)

    s0 = ArgParseSettings()  # a "parent" structure e.g. one with some useful set of rules
                             # which we want to extend

    # So we just add a simple table
    @add_arg_table s0 begin
        "--parent-flag", "-o"
            action = "store_true"
            help = "parent flag"
        "--flag"
            action = "store_true"
            help = "another parent flag"
        "parent-argument"
            help = "parent argument"
    end

    s = ArgParseSettings("Test 4 for ArgParse.jl",
                         add_help = false,
                         version = "Version 1.0")

    import_settings(s, s0)       # now s has all of s0 arguments (except help/version)

    @add_arg_table s begin
        "-o"                       # this will partially override s0's --parent-flag
            action = :store_true
            help = "child flag"
        "--flag"                   # this will fully override s0's --flag
            action = :store_true
            help = "another child flag"
    end

    s.exc_handler = (settings, err)->error(err.text)

    parsed_args = parse_args(args, s)
end

@test ap_test4([]) == (String=>Any)["parent-flag"=>false, "o"=>false, "flag"=>false, "parent-argument"=>nothing]
@test ap_test4(["-o", "X"]) == (String=>Any)["parent-flag"=>false, "o"=>true, "flag"=>false, "parent-argument"=>"X"]
@test_throws ap_test4(["-h"])
@test_throws ap_test4_fails([])


text:
"""
    KT_residuals

Compute Kuhn-Turck residuals

# Arguments
- `df`: gradient
- `lambda`: Lagrangian multipliers with respect to the constraints
- `mult_x_U`: reduced cost with respect to the upper bounds
- `mult_x_L`: reduced cost with respect to the lower bounds
- `Jac`: Jacobian matrix
- `norm`: whether the residual is normalized or not
"""
function KT_residuals(
    df::Tv, lambda::Tv, mult_x_U::Tv, mult_x_L::Tv, Jac::Tm
) where {T, Tv<:AbstractArray{T}, Tm<:AbstractMatrix{T}}
    KT_res = norm(df - Jac' * lambda - mult_x_U - mult_x_L)
    scalar = max(1.0, norm(df))
    for i = 1:size(Jac,1)
        scalar = max(scalar, abs(lambda[i]) * norm(Jac[i,:]))
    end
    return KT_res / scalar
end

"""
    norm_complementarity

Compute the normalized complementeraity
"""
function norm_complementarity(
    E::Tv, g_L::Tv, g_U::Tv, x::Tv, x_L::Tv, x_U::Tv,
    lambda::Tv, mult_x_U::Tv, mult_x_L::Tv,
    p = Inf
) where {T, Tv <: AbstractArray{T}}
    m = length(E)
    compl = Tv(undef, m)
    denom = 0.0
    for i = 1:m
        if g_L[i] == g_U[i]
            compl[i] = 0.0
        else
            compl[i] = min(E[i] - g_L[i], g_U[i] - E[i]) * lambda[i]
            denom += lambda[i]^2
        end
    end
    return norm(compl, p) / (1 + sqrt(denom))
end

"""
    norm_violations

Compute the normalized constraint violation
"""
function norm_violations(
    E::Tv, g_L::Tv, g_U::Tv, x::Tv, x_L::Tv, x_U::Tv, p = 1
) where {T, Tv <: AbstractArray{T}}

    m = length(E)
    n = length(x)
    viol = Tv(undef, m+n)
    fill!(viol, 0.0)
    for i = 1:m
        if E[i] > g_U[i]
            viol[i] = E[i] - g_U[i]
        elseif E[i] < g_L[i]
            viol[i] = g_L[i] - E[i]
        end
    end
    for j = 1:n
        if x[j] > x_U[j]
            viol[m+j] = x[j] - x_U[j]
        elseif x[j] < x_L[j]
            viol[m+j] = x_L[j] - x[j]
        end
    end
    return norm(viol, p)
end


text:
export combinePRS,
    combinePRSbyType,
    makeModelFrameCC,
    estimate_h2,
    estimate_neff,
    get_chi_markers,
    postfix_names!,
    snp_join,
    read_chi_markers

using DataFrames
using CSV

"add up PRS scores from chromosome files in dir and write to outfile"
function combinePRS(dir::AbstractString, outfile)

    columns = []
    fh = open(joinpath(dir, "prstypes.txt"))
    for line in eachline(fh)
        push!(columns, Symbol(chomp(line)))
    end
    close(fh)

    prs = Dict{String, Array{Float64}}()
    files = [joinpath(dir, f) for f in readdir(dir) if endswith(f, ".prs")]
    for f in files
        for line in eachline(f)
            cols = split(chomp(line), "\t")
            pn = cols[1]
            vals = parse.(Float64, cols[5:end])
            haskey(prs, pn) || (prs[pn] = zeros(Float64, length(vals)))
            prs[pn] .+= vals
        end
    end

    fh = open(outfile, "w")
    columns == [] || println(fh, join(vcat([:PN], columns), "\t"))
    for (k,vs) in sort(collect(prs))
        println(fh, k, "\t", join([string(round(v; sigdigits=4)) for v in vs], "\t"))
    end
    close(fh)
end

"for each type of prs weight in columns, create file containing weights by chromosome in columns"
function combinePRSbyType(dir::AbstractString, outdir)

    columns = []
    fh = open(joinpath(dir, "prstypes.txt"))
    for line in eachline(fh)
        push!(columns, Symbol(chomp(line)))
    end
    close(fh)

    chroms = ["chr" * string(i) for i = 1:22]

    prs = Dict{String, Dict{String, Array{Float64}}}()
    for chrom in chroms
        prs[chrom] = Dict{String, Array{Float64}}()  # prs[chrom][pn][type] = weight
        prs_chrom = prs[chrom]
        path = joinpath(dir, chrom * ".prs")
        for line in eachline(path)
            cols = split(chomp(line), "\t")
            pn = cols[1]
            vals = parse.(Float64, cols[5:end])
            @assert length(vals) == length(columns) "Wrong number of columns for $pn in $path"
            prs_chrom[pn] = vals
        end
    end

    idx = 0
    for column in columns
        idx += 1
        path = joinpath(outdir, string(column) * ".prs")
        fh = open(path, "w")
        println(fh, "PN\t", join([string(column) * "_" *string(chrom) for chrom in chroms], "\t"))
        pns = sort(collect(keys(prs[chroms[1]])))
        for pn in pns
            println(fh, pn, "\t", join([string(round(prs[chrom][pn][idx]; sigdigits=4)) for chrom in chroms], "\t"))
        end
        close(fh)
    end
end

function makeModelFrameCC(prs_path, pheno_path; header = 1)
    prs = CSV.read(prs_path; delim = "\t", header = header)

    fitted_path = joinpath(pheno_path, "results/casecontrol", basename(pheno_path), "fitted.auto.tsv")
    fitted = CSV.read(fitted_path; delim = "\t", header = [:PN, :fitted])

    case_path = joinpath(pheno_path, "results/casecontrol", basename(pheno_path), "aff.auto.pnlist")
    case = CSV.read(case_path; delim = "\t", header = [:PN])

    pns = Set(case.PN)

    m = join(prs, fitted, on = :PN)
    m[!,:affected] .= 0
    m.affected[in(pns).(m.PN)] .= 1

    m
end

"""
    estimate_h2(z, lds, n)

Estimate total heritability ``h^2`` for markers using LD score
regression with intercept forced to 1.

# Arguments

- `lds::Vector`: LD score as computed by `ldscore(LDMatrix)`
- `n::Integer`: effective number of samples (``4 n_0*n_1/(n_0+n_1)`` for case-control)
"""
function estimate_h2(z, lds, n)
    M = length(lds)
    meanlds = sum(lds) / M

    h2 = (mean(z .^2) - 1) / n / meanlds * M
end

"""
    estimate_neff(beta, sebeta, freq)

Estimate effective number of samples in case-control study.

Effective number of samples is total number of samples for a cohort
with same number of samples and controls, that would result in the
observed `sebeta`.

`freq` can be minor or major allele freq,
as only `freq*(1-freq)` is being used.
"""
function estimate_neff(beta, sebeta, freq)
    neff = 2. ./ (sebeta .^2 .* freq .* (1.0 .- freq))
    median(neff[isfinite.(neff)])
end

"""
    snp_join(df1, df2; on = [:Chrom, :Pos],
             alleles1 = [:A1_1, :A2_1], alleles2 = [:A1_2, :A2_2],
             matchcol = :sign)

Join data frames `df1` and `df2` and match variants.

Data frames `df1` and `df2` are first matched by columns specified in
`on`. Then column `matchcol` is set to

- `matchcol == 1` if values for `alleles1` and `alleles2` match
- `matchcol == -1` if values for `alleles1` and `alleles2` are swapped
- `matchcol == 0` otherwise

Note: Swapped means simply swapping of alleles, not applying reverse
complement as done in other implementations.

"""
function snp_join(df1, df2; on = [:Chrom, :Pos], alleles1 = [:A1_1, :A2_1],
                  alleles2 = [:A1_2, :A2_2], matchcol = :sign)

    m = innerjoin(df1, df2, on = on)
    m[!, matchcol] .= 0
    m[(m[!, alleles1[1]] .== m[!, alleles2[1]]) .& (m[!, alleles1[2]] .== m[!, alleles2[2]]), matchcol] .= 1
    m[(m[!, alleles1[1]] .== m[!, alleles2[2]]) .& (m[!, alleles1[2]] .== m[!, alleles2[1]]), matchcol] .= -1

    n_match = sum(m[!, matchcol] .== 1)
    n_mismatch = sum(m[!, matchcol] .== 0)
    n_flip = sum(m[!, matchcol] .== -1)

    println("Found $n_match matching, $n_flip flipped, and $n_mismatch mismatching SNPs")
    m
end

"""
read chi markers

return data frame with columns
Chrom, Name, Pos, A1, A2
"""
function read_chi_markers(chidir::AbstractString, chrom)
    path = joinpath(chidir, chrom * ".markers")
    tmp = CSV.read(path; delim="\t")[!,1:5]
    rename!(tmp, [:Chrom_, :Name, :Pos, :A1, :A2])
    tmp[!,:Chrom] .= "chr" .* string.(tmp.Chrom_)
    tmp[!, [:Chrom, :Name, :Pos, :A1, :A2]]
end

function get_chi_markers(dir::AbstractString, markers::Vector)
    requested = in(Set(markers))
    files = sort([joinpath(dir, f) for f in readdir(dir) if endswith(f, ".markers")])
    tmp = CSV.read(files[1]; delim="\t")
    tmp = tmp[requested.(tmp.name), 1:5]
    rename!(tmp, [:Chrom, :Name, :Pos, :A1, :A2])
    res = tmp
    for f in files[2:end]
        tmp = CSV.read(f; delim="\t")
        tmp = tmp[requested.(tmp.name),1:5]
        rename!(tmp, [:Chrom, :Name, :Pos, :A1, :A2])
        res = vcat(res, tmp)
    end
    res
end

function postfix_names!(df::DataFrame, postfix::String)
    oldnames = [string(x) for x in names(df)]
    newnames = [ x * postfix for x in oldnames]
    rename!(df, newnames)
end

# do not export this
# should be called as PRS.pkgpath
pkgpath(paths...) = normpath(joinpath(@__DIR__, "..", paths...))


text:
abstract ReferenceSystem
abstract IndependentReferenceSystem <: ReferenceSystem
abstract DependentReferenceSystem <: ReferenceSystem

abstract Perturbation
abstract NFEPerturbation <: Perturbation
abstract TBPerturbation <: Perturbation

abstract PseudoPotential

typealias Pairpotential Function


text:
# Visualize 2 sessions in Arena using local DistributedFactorGraphs types
using UUIDs
using IncrementalInference
using RoME
using DistributedFactorGraphs
using Amphitheatre

##
# create a local fg hexslam
# start with an empty factor graph object

fg1 = LightDFG{SolverParams}(solverParams=SolverParams())
# Add the first pose :x0
v = addVariable!(fg1, :x0, Pose2, tags=[:POSE])
# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg1, [:x0], Prior( MvNormal([0; 0; 0], Matrix(Diagonal([0.1;0.1;0.05].^2)) )))
romeVis1 = BasicFactorGraphPose(fg1, meanmax=:mean, poseProp = plDrawProp(0.3, 0.1, RGBA(0,1,1,0.5)))

fg2 = LightDFG{SolverParams}(description="LG DFG",userId="a_user",robotId="a_robot",solverParams=SolverParams())
# Add the first pose :x0
addVariable!(fg2, :x0, Pose2, tags=[:POSE])
# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg2, [:x0], Prior( MvNormal([4; 0; -pi], Matrix(Diagonal([0.1;0.1;0.05].^2)) )))
romeVis2 = BasicFactorGraphPose(fg2, meanmax=:mean)

# Create AbstractVarsVis container to hold different visualizers
visdatasets = AbstractAmphitheatre[romeVis1, romeVis2]

vis, vistask = visualize(visdatasets, start_browser=true)

##
for session in visdatasets
    fg = session.fg

    # Drive around in a hexagon
    for i in 0:5
        psym = Symbol("x$i")
        nsym = Symbol("x$(i+1)")
        addVariable!(fg, nsym, RoME.Pose2, tags=[:POSE])
        pp = Pose2Pose2(MvNormal([1.0;0;pi/3+0.02], Matrix(Diagonal([0.01;0.01;0.1].^2))))
        addFactor!(fg, [psym;nsym], pp )
        sleep(0.5)
    end

    # Add landmarks with Bearing range measurements
    addVariable!(fg, :l1, RoME.Point2, tags=[:LANDMARK])
    p2br = Pose2Point2BearingRange(IIF.Normal(0,0.01),IIF.Normal(2.0,0.1))
    addFactor!(fg, [:x6; :l1], p2br)

    # Add landmarks with Bearing range measurements
    p2br2 = Pose2Point2BearingRange(IIF.Normal(0,0.01),IIF.Normal(2.0,0.1))
    addFactor!(fg, [:x6; :l1], p2br2)
end

addFactor!(fg1, [:l1], PriorPoint2(MvNormal([2.,0], Matrix(Diagonal([0.001, 0.001].^2))) ))
addFactor!(fg2, [:l1], PriorPoint2(MvNormal([2.,0], Matrix(Diagonal([0.001, 0.001].^2))) ))

## solve
@async solveTree!(fg1)
@async solveTree!(fg2)

##
@info "To stop call stopAmphiVis!()"
#
stopAmphiVis!()
##
# To restart:
vistask = Amphitheatre.restartVisualize(vis, visdatasets)

max_stars_repo_name: Affie/Amphitheatre

max_issues_repo_name: Affie/Amphitheatre