text:
module DataFramesMeta

using DataFrames, Tables

# Basics:
export @with, @where, @orderby, @transform, @by, @based_on, @select

include("linqmacro.jl")
include("byrow.jl")

##############################################################################
##
## @with
##
##############################################################################

function addkey!(membernames, nam)
    if !haskey(membernames, nam)
        membernames[nam] = gensym()
    end
    membernames[nam]
end

onearg(e, f) = e.head == :call && length(e.args) == 2 && e.args[1] == f
mapexpr(f, e) = Expr(e.head, map(f, e.args)...)

replace_syms!(x, membernames) = x
replace_syms!(q::QuoteNode, membernames) =
    replace_syms!(Meta.quot(q.value), membernames)
replace_syms!(e::Expr, membernames) =
    if onearg(e, :^)
        e.args[2]
    elseif onearg(e, :_I_)
        @warn "_I_() for escaping variables is deprecated, use cols() instead"
        addkey!(membernames, :($(e.args[2])))
    elseif onearg(e, :cols)
        addkey!(membernames, :($(e.args[2])))
    elseif e.head == :quote
        addkey!(membernames, Meta.quot(e.args[1]) )
    elseif e.head == :.
        replace_dotted!(e, membernames)
    else
        mapexpr(x -> replace_syms!(x, membernames), e)
    end

protect_replace_syms!(e, membernames) = e
function protect_replace_syms!(e::Expr, membernames)
    if e.head == :quote
        e
    else
        replace_syms!(e, membernames)
    end
end

function replace_dotted!(e, membernames)
    x_new = replace_syms!(e.args[1], membernames)
    y_new = protect_replace_syms!(e.args[2], membernames)
    Expr(:., x_new, y_new)
end

function with_helper(d, body)
    membernames = Dict{Any, Symbol}()
    funname = gensym()
    body = replace_syms!(body, membernames)
    if isempty(membernames)
        body
    else
        quote
            function $funname($(values(membernames)...))
                $body
            end
            $funname($((:($d[!, $key]) for key in keys(membernames))...))
        end
    end
end

function with_anonymous(body)
    d = gensym()
    :($d -> $(with_helper(d, body)))
end

"""
    @with(d, expr)

`@with` allows DataFrame columns keys to be referenced as symbols.

### Arguments

* `d` : an AbstractDataFrame type
* `expr` : the expression to evaluate in `d`

### Details

`@with` works by parsing the expression body for all columns indicated
by symbols (e.g. `:colA`). Then, a function is created that wraps the
body and passes the columns as function arguments. This function is
then called. Operations are efficient because:

- A pseudo-anonymous function is defined, so types are stable.
- Columns are passed as references, eliminating DataFrame indexing.

The following

```julia
@with(d, :a .+ :b .+ 1)
```

becomes

```julia
tempfun(a, b) = a .+ b .+ 1
tempfun(d[!, :a], d[!, :b])
```

All of the other DataFramesMeta macros are based on `@with`.

If an expression is wrapped in `^(expr)`, `expr` gets passed through untouched.
If an expression is wrapped in  `cols(expr)`, the column is referenced by the
variable `expr` rather than a symbol.

### Examples

```jldoctest
julia> using DataFramesMeta

julia> y = 3;

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> x = [2, 1, 0];

julia> @with(df, :y .+ 1)
3-element Array{Int64,1}:
 3
 2
 3

julia> @with(df, :x + x)
3-element Array{Int64,1}:
 3
 3
 3

julia> @with df begin
            res = 0.0
            for i in 1:length(:x)
                res += :x[i] * :y[i]
            end
            res
        end
10.0

julia> @with(df, df[:x .> 1, ^(:y)]) # The ^ means leave the :y alone
2-element Array{Int64,1}:
 1
 2

julia> colref = :x;

julia> @with(df, :y + cols(colref)) # Equivalent to df[!, :y] + df[!, colref]
 3
 3
 5
```

`@with` creates a function, so scope within `@with` is a local scope.
Variables in the parent can be read. Writing to variables in the parent scope
differs depending on the type of scope of the parent. If the parent scope is a
global scope, then a variable cannot be assigned without using the `global` keyword.
If the parent scope is a local scope (inside a function or let block for example),
the `global` keyword is not needed to assign to that parent scope.

"""
macro with(d, body)
    esc(with_helper(d, body))
end

##############################################################################
##
## @where - select row subsets
##
##############################################################################

where(d::AbstractDataFrame, arg) = d[arg, :]
where(d::AbstractDataFrame, f::Function) = d[f(d), :]
where(g::GroupedDataFrame, f::Function) = g[Bool[f(x) for x in g]]

and(x, y) = :($x .& $y)

function where_helper(d, args...)
    :($where($d, $(with_anonymous(reduce(and, args)))))
end

"""
    @where(d, i...)

Select row subsets in AbstractDataFrames or groups in GroupedDataFrames.

### Arguments

* `d` : an AbstractDataFrame or GroupedDataFrame
* `i...` : expression for selecting rows

Multiple `i` expressions are "and-ed" together.

### Examples

```jldoctest
julia> using DataFramesMeta, DataFrames

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> x = [2, 1, 0];

julia> @where(df, :x .> 1)
2×2 DataFrame
│ Row │ x │ y │
├─────┼───┼───┤
│ 1   │ 2 │ 1 │
│ 2   │ 3 │ 2 │

julia> @where(df, :x .> x)
2×2 DataFrame
│ Row │ x │ y │
├─────┼───┼───┤
│ 1   │ 2 │ 1 │
│ 2   │ 3 │ 2 │

julia> @where(df, :x .> x, :y .== 3)
0×2 DataFrame

julia> d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1,
                                    2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @where(d, :x .== 3)
5×2 DataFrame
│ Row │ n  │ x │
├─────┼────┼───┤
│ 1   │ 1  │ 3 │
│ 2   │ 2  │ 3 │
│ 3   │ 3  │ 3 │
│ 4   │ 4  │ 3 │
│ 5   │ 17 │ 3 │

julia> @where(g, length(:x) > 5)   # pick out some groups
GroupedDataFrame  2 groups with keys: Symbol[:x]
First Group:
9×2 SubDataFrame{Array{Int64,1}}
│ Row │ n  │ x │
├─────┼────┼───┤
│ 1   │ 5  │ 1 │
│ 2   │ 6  │ 1 │
│ 3   │ 7  │ 1 │
│ 4   │ 9  │ 1 │
│ 5   │ 10 │ 1 │
│ 6   │ 12 │ 1 │
│ 7   │ 13 │ 1 │
│ 8   │ 18 │ 1 │
│ 9   │ 19 │ 1 │
⋮
Last Group:
6×2 SubDataFrame{Array{Int64,1}}
│ Row │ n  │ x │
├─────┼────┼───┤
│ 1   │ 8  │ 2 │
│ 2   │ 11 │ 2 │
│ 3   │ 14 │ 2 │
│ 4   │ 15 │ 2 │
│ 5   │ 16 │ 2 │
│ 6   │ 20 │ 2 │
```
"""
macro where(d, args...)
    esc(where_helper(d, args...))
end

##############################################################################
##
## select - select columns
##
##############################################################################

select(d::AbstractDataFrame, arg) = d[!, arg]

##############################################################################
##
## @orderby
##
##############################################################################

# needed on Julia 1.0 till #1489 in DataFrames is merged
orderby(d::DataFrame, arg::DataFrame) = d[sortperm(arg), :]

function orderby(d::AbstractDataFrame, args...)
    D = typeof(d)(args...)
    d[sortperm(D), :]
end

orderby(d::AbstractDataFrame, f::Function) = d[sortperm(f(d)), :]
orderby(g::GroupedDataFrame, f::Function) = g[sortperm([f(x) for x in g])]

orderbyconstructor(d::AbstractDataFrame) = (x...) -> DataFrame(Any[x...], Symbol.(1:length(x)))
orderbyconstructor(d) = x -> x

function orderby_helper(d, args...)
    _D = gensym()
    quote
        let $_D = $d
            $orderby($_D, $(with_anonymous(:($orderbyconstructor($_D)($(args...))))))
        end
    end
end

"""
    @orderby(d, i...)

Sort by criteria. Normally used to sort groups in GroupedDataFrames.

### Arguments

* `d` : an AbstractDataFrame or GroupedDataFrame
* `i...` : expression for sorting

### Examples

```jldoctest
julia> using DataFrames, DataFramesMeta, Statistics

julia> d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1,
                                    2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @orderby(g, mean(:n))
GroupedDataFrame  3 groups with keys: Symbol[:x]
First Group:
5×2 SubDataFrame{Array{Int64,1}}
│ Row │ n  │ x │
├─────┼────┼───┤
│ 1   │ 1  │ 3 │
│ 2   │ 2  │ 3 │
│ 3   │ 3  │ 3 │
│ 4   │ 4  │ 3 │
│ 5   │ 17 │ 3 │
⋮
Last Group:
6×2 SubDataFrame{Array{Int64,1}}
│ Row │ n  │ x │
├─────┼────┼───┤
│ 1   │ 8  │ 2 │
│ 2   │ 11 │ 2 │
│ 3   │ 14 │ 2 │
│ 4   │ 15 │ 2 │
│ 5   │ 16 │ 2 │
│ 6   │ 20 │ 2 │
```

"""
macro orderby(d, args...)
    # I don't esc just the input because I want _DF to be visible to the user
    esc(orderby_helper(d, args...))
end

##############################################################################
##
## transform & @transform
##
##############################################################################

function transform(d::AbstractDataFrame; kwargs...)
    result = copy(d)
    for (k, v) in kwargs
        result[!, k] = isa(v, Function) ? v(d) : v
    end
    return result
end

function transform(g::GroupedDataFrame; kwargs...)
    result = DataFrame(g)
    ends = cumsum(Int[size(g[i],1) for i in 1:length(g)])
    starts = [1; 1 .+ ends[1:end-1]]
    lengths = [ends[i] - starts[i] + 1 for i in 1:length(starts)]
    for (k, v) in kwargs
        first = v(g[1])
        if first isa AbstractVector
            t = _transform!(Tables.allocatecolumn(eltype(first), size(result, 1)),
                            first, 1, g, v, starts, ends)
        else
            t = _transform!(Tables.allocatecolumn(typeof(first), size(result, 1)),
                            first, 1, g, v, starts, ends)
        end
        result[!, k] = t
    end
    return result
end

function _transform!(t::AbstractVector, first::AbstractVector, start::Int,
                     g::GroupedDataFrame, v::Function, starts::Vector, ends::Vector)
    @inline function fill_column!(t::AbstractVector, out, startpoint::Int, endpoint::Int,
                                      len::Int)
        if !(out isa AbstractVector)
            throw(ArgumentError("Return value must be an `AbstractVector` for all groups or" *
                                "for none of them"))
        elseif length(out) != len
            throw(ArgumentError("If a function returns a vector, the result " *
                                "must have the same length as the groups it operates on"))
        end
        eltypout = eltype(out)
        T = eltype(t)
        if eltypout <: T || (newtype = promote_type(eltypout, T)) <: T
           t[startpoint:endpoint] = out
            return nothing
        else
            return newtype
        end
        return nothing
    end

    # handle the first case
    newtype_first = fill_column!(t, first, starts[start], ends[start], size(g[start], 1))
    @assert newtype_first === nothing
    @inbounds for i in (start+1):length(g)
        out = v(g[i])
        newtype = fill_column!(t, out, starts[i], ends[i], size(g[i], 1))
        if newtype !== nothing
             t = copyto!(Tables.allocatecolumn(newtype, length(t)),
                         1, t, 1, ends[i-1])
             _transform!(t, out, i, g, v, starts, ends)
         end
    end
    return t
end

function _transform!(t::AbstractVector, first::Any, start::Int,
                     g::GroupedDataFrame, v::Function, starts::Vector, ends::Vector)
    @inline function fill_column!(t::AbstractVector, out, startpoint::Int, endpoint::Int)
        if out isa AbstractVector
            throw(ArgumentError("Return value must be an `AbstractVector` for all groups or" *
                                 "for none of them"))
        end
        typout = typeof(out)
        T = eltype(t)
        if typout <: T || (newtype = promote_type(typout, T)) <: T
            t[startpoint:endpoint] .= Ref(out)
            return nothing
        else
            return newtype
        end
    end
    # handle the first case
    newtype_first = fill_column!(t, first, starts[start], ends[start])
    @assert newtype_first === nothing
    @inbounds for i in (start+1):length(g)
        out = v(g[i])
        newtype = fill_column!(t, out, starts[i], ends[i])
        if newtype !== nothing
             t = copyto!(Tables.allocatecolumn(newtype, length(t)),
                         1, t, 1, ends[i-1])
             _transform!(t, out, i, g, v, starts, ends)
         end
    end
    return t
end

function transform_helper(x, args...)
    quote
        $transform($x, $(map(args) do kw
            Expr(:kw, kw.args[1], with_anonymous(kw.args[2]))
        end...) )
    end
end

"""
    @transform(d, i...)

Add additional columns or keys based on keyword arguments.

### Arguments

* `d` : an `AbstractDataFrame`, or `GroupedDataFrame`
* `i...` : keyword arguments defining new columns or keys

### Returns

* `::AbstractDataFrame` or `::GroupedDataFrame`

### Examples

```jldoctest
julia> using DataFramesMeta, DataFrames

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> @transform(df, a = 2 * :A, x = :A .+ :B)
3×4 DataFrame
│ Row │ A │ B │ a │ x │
├─────┼───┼───┼───┼───┤
│ 1   │ 1 │ 2 │ 2 │ 3 │
│ 2   │ 2 │ 1 │ 4 │ 3 │
│ 3   │ 3 │ 2 │ 6 │ 5 │
```

"""
macro transform(x, args...)
    esc(transform_helper(x, args...))
end

##############################################################################
##
## @based_on - summarize a grouping operation
##
##############################################################################

function based_on_helper(x, args...)
    with_args =
        with_anonymous(:($DataFrame($(map(replace_equals_with_kw, args)...))))
    :(DataFrames.combine($with_args, $x))
end

"""
    @based_on(g, i...)

Summarize a grouping operation

### Arguments

* `g` : a GroupedDataFrame
* `i...` : keyword arguments defining new columns

### Examples

```jldoctest
julia> using DataFramesMeta, DataFrames

julia> d = DataFrame(
            n = 1:20,
            x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @based_on(g, nsum = sum(:n))
3×2 DataFrame
│ Row │ x │ nsum │
├─────┼───┼──────┤
│ 1   │ 1 │ 99   │
│ 2   │ 2 │ 84   │
│ 3   │ 3 │ 27   │

julia> @based_on(g, x2 = 2 * :x, nsum = sum(:n))
20×3 DataFrame
│ Row │ x │ x2 │ nsum │
├─────┼───┼────┼──────┤
│ 1   │ 1 │ 2  │ 99   │
│ 2   │ 1 │ 2  │ 99   │
│ 3   │ 1 │ 2  │ 99   │
│ 4   │ 1 │ 2  │ 99   │
│ 5   │ 1 │ 2  │ 99   │
│ 6   │ 1 │ 2  │ 99   │
│ 7   │ 1 │ 2  │ 99   │
│ 8   │ 1 │ 2  │ 99   │
⋮
│ 12  │ 2 │ 4  │ 84   │
│ 13  │ 2 │ 4  │ 84   │
│ 14  │ 2 │ 4  │ 84   │
│ 15  │ 2 │ 4  │ 84   │
│ 16  │ 3 │ 6  │ 27   │
│ 17  │ 3 │ 6  │ 27   │
│ 18  │ 3 │ 6  │ 27   │
│ 19  │ 3 │ 6  │ 27   │
│ 20  │ 3 │ 6  │ 27   │
```
"""
macro based_on(x, args...)
    esc(based_on_helper(x, args...))
end

##############################################################################
##
## @by - grouping
##
##############################################################################

function by_helper(x, what, args...)
    :(DataFrames.combine($(with_anonymous(:($DataFrame($(map(replace_equals_with_kw, args)...))))),
              DataFrames.groupby($x, $what)))
end

"""
    @by(d::AbstractDataFrame, cols, e...)

Split-apply-combine in one step.

### Arguments

* `d` : an AbstractDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `e` :  keyword arguments specifying new columns in terms of column groupings

### Returns

* `::DataFrame`

### Examples

```jldoctest
julia> using DataFrames, DataFramesMeta, Statistics

julia> df = DataFrame(
            a = repeat(1:4, outer = 2),
            b = repeat(2:-1:1, outer = 4),
            c = randn(8));

julia> @by(df, :a, d = sum(:c))
4×2 DataFrame
│ Row │ a │ d        │
├─────┼───┼──────────┤
│ 1   │ 1 │ 1.27638  │
│ 2   │ 2 │ 1.00951  │
│ 3   │ 3 │ 1.48328  │
│ 4   │ 4 │ -2.42621 │

julia> @by(df, :a, d = 2 * :c)
8×2 DataFrame
│ Row │ a │ d         │
├─────┼───┼───────────┤
│ 1   │ 1 │ 1.22982   │
│ 2   │ 1 │ 1.32294   │
│ 3   │ 2 │ 1.93664   │
│ 4   │ 2 │ 0.0823819 │
│ 5   │ 3 │ -0.670512 │
│ 6   │ 3 │ 3.63708   │
│ 7   │ 4 │ -3.06436  │
│ 8   │ 4 │ -1.78806  │

julia> @by(df, :a, c_sum = sum(:c), c_mean = mean(:c))
4×3 DataFrame
│ Row │ a │ c_sum    │ c_mean   │
├─────┼───┼──────────┼──────────┤
│ 1   │ 1 │ 1.27638  │ 0.63819  │
│ 2   │ 2 │ 1.00951  │ 0.504755 │
│ 3   │ 3 │ 1.48328  │ 0.741642 │
│ 4   │ 4 │ -2.42621 │ -1.2131  │

julia> @by(df, :a, c = :c, c_mean = mean(:c))
8×3 DataFrame
│ Row │ a │ c         │ c_mean   │
├─────┼───┼───────────┼──────────┤
│ 1   │ 1 │ 0.61491   │ 0.63819  │
│ 2   │ 1 │ 0.66147   │ 0.63819  │
│ 3   │ 2 │ 0.968319  │ 0.504755 │
│ 4   │ 2 │ 0.041191  │ 0.504755 │
│ 5   │ 3 │ -0.335256 │ 0.741642 │
│ 6   │ 3 │ 1.81854   │ 0.741642 │
│ 7   │ 4 │ -1.53218  │ -1.2131  │
│ 8   │ 4 │ -0.894029 │ -1.2131  │
```
"""
macro by(x, what, args...)
    esc(by_helper(x, what, args...))
end

##############################################################################
##
## @select - select and transform columns
##
##############################################################################

function select(d::AbstractDataFrame; kwargs...)
    result = typeof(d)()
    for (k, v) in kwargs
        result[!, k] = v
    end
    return result
end

function replace_equals_with_kw(e)
    if e.head == :(=)
        Expr(:kw, e.args[1], e.args[2])
    else
        e
    end
end

expandargs(x) = x
expandargs(q::QuoteNode) = Expr(:kw, q.value, q)
function expandargs(e::Expr)
    if e.head == :quote
        Expr(:kw, e.args[1], e)
    else
        replace_equals_with_kw(e)
    end
end

function select_helper(x, args...)
    DF = gensym()
    select_args = with_helper(DF, :($select($DF, $(map(expandargs, args)...))))
    quote
        let $DF = $x
            $(with_helper(DF, :($select($DF, $(map(expandargs, args)...)))))
        end
    end
end

"""
    @select(d, e...)

Select and transform columns.

### Arguments

* `d` : an AbstractDataFrame
* `e` :  keyword arguments specifying new columns in terms of existing columns
  or symbols to specify existing columns

### Returns

* `::AbstractDataFrame`

### Examples

```jldoctest
julia> using DataFrames, DataFramesMeta

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = randn(8))
8×3 DataFrame
│ Row │ a │ b │ c         │
├─────┼───┼───┼───────────┤
│ 1   │ 1 │ 2 │ -0.354685 │
│ 2   │ 2 │ 1 │ 0.287631  │
│ 3   │ 3 │ 2 │ -0.918007 │
│ 4   │ 4 │ 1 │ -0.352519 │
│ 5   │ 1 │ 2 │ 0.743501  │
│ 6   │ 2 │ 1 │ -1.27415  │
│ 7   │ 3 │ 2 │ 0.258456  │
│ 8   │ 4 │ 1 │ -0.460486 │

julia> @select(df, :c, :a)
8×2 DataFrame
│ Row │ c         │ a │
├─────┼───────────┼───┤
│ 1   │ -0.354685 │ 1 │
│ 2   │ 0.287631  │ 2 │
│ 3   │ -0.918007 │ 3 │
│ 4   │ -0.352519 │ 4 │
│ 5   │ 0.743501  │ 1 │
│ 6   │ -1.27415  │ 2 │
│ 7   │ 0.258456  │ 3 │
│ 8   │ -0.460486 │ 4 │

julia> @select(df, :c, x = :b + :c)
8×2 DataFrame
│ Row │ c         │ x         │
├─────┼───────────┼───────────┤
│ 1   │ -0.354685 │ 1.64531   │
│ 2   │ 0.287631  │ 1.28763   │
│ 3   │ -0.918007 │ 1.08199   │
│ 4   │ -0.352519 │ 0.647481  │
│ 5   │ 0.743501  │ 2.7435    │
│ 6   │ -1.27415  │ -0.274145 │
│ 7   │ 0.258456  │ 2.25846   │
│ 8   │ -0.460486 │ 0.539514  │
```
"""
macro select(x, args...)
    esc(select_helper(x, args...))
end

end # module


text:
using Test

@testset "/src/BellScenario/DichotomicLocalPolytope.jl" begin

using BellScenario: DichotomicLocalPolytope

@testset "bipartite_no_signaling_strategies()" begin

    @testset "16 strategies are returned" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies()
        @test size(s) == (16,)

        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies(subspace="p12")
        @test size(s) == (16,)
    end

    @testset "strategies are 4x4" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies()
        for i in 1:size(s)[1]
            @test size(s[i]) == (4,4)
        end
    end

    @testset "strategies are 3x4 in p12 subspace" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies(subspace="p12")
        for i in 1:size(s)[1]
            @test size(s[i]) == (3,4)
        end
    end

    @testset "returned strategies are unique" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies()
        for i in 1:size(s)[1]
            strategy = pop!(s)
            non_unique_count = 0
            for j in 1:size(s)[1]
                if strategy == s[j]
                    non_unique_count += 1
                end
            end
            @test non_unique_count == 0
        end
    end

    @testset "returned p12 strategies are unique" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies(subspace="p12")
        for i in 1:size(s)[1]
            strategy = pop!(s)
            non_unique_count = 0
            for j in 1:size(s)[1]
                if strategy == s[j]
                    non_unique_count += 1
                end
            end
            @test non_unique_count == 0
        end
    end

    @testset "spot checking returned values" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies()
        @test s[1] == [
            1 0 0 0;
            0 1 0 0;
            0 0 1 0;
            0 0 0 1;
        ]
        @test s[2] == [
            1 1 0 0;
            0 0 0 0;
            0 0 1 1;
            0 0 0 0;
        ]
        @test s[15] == [
            0 0 0 0;
            0 0 1 1;
            0 0 0 0;
            1 1 0 0;
        ]
    end

    @testset "spot checking returned values for p12" begin
        s = DichotomicLocalPolytope.bipartite_no_signaling_strategies(subspace="p12")
        @test s[1] == [
            1 0 0 0;
            0 1 0 0;
            0 0 1 0;
        ]
        @test s[2] == [
            1 1 0 0;
            0 0 0 0;
            0 0 1 1;
        ]
        @test s[15] == [
            0 0 0 0;
            0 0 1 1;
            0 0 0 0;
        ]
    end
end

@testset "deterministic_fixed_direction_signaling_behaviors()" begin
    test_cases = [
        ("p16",16), # (subspace, dimension of returned behaviors)
        ("p12",12),
        ("p10",10)
    ]

    @testset "subspace=$(case[1])" for case in test_cases
        vertices = DichotomicLocalPolytope.deterministic_fixed_direction_signaling_behaviors(subspace=case[1])

        @test size(vertices) == (64,)

        @testset "behaviors are unique" for v_test in vertices
            @test size(v_test) == (case[2],)

            duplicate_count = 0
            for v_match in vertices
                if v_test == v_match
                    duplicate_count += 1
                end
            end
            @test duplicate_count == 1
        end
    end

    @testset "verify equivalence between subspaces" begin
        p16_vertices = DichotomicLocalPolytope.deterministic_fixed_direction_signaling_behaviors(subspace="p16")
        p12_vertices = DichotomicLocalPolytope.deterministic_fixed_direction_signaling_behaviors(subspace="p12")
        p10_vertices = DichotomicLocalPolytope.deterministic_fixed_direction_signaling_behaviors(subspace="p10")
        for i in 1:64
            catv_16 = cat([1], p16_vertices[i], dims = 1)
            catv_12 = cat([1], p12_vertices[i], dims = 1)
            catv_10 = cat([1], p10_vertices[i], dims = 1)

            @test catv_10 == Behavior.gen_to_fd_proj((2,2),(2,2))*catv_16
            @test catv_12 == Behavior.gen_to_norm_proj((2,2),(2,2))*catv_16
            @test catv_16 == Behavior.fd_to_gen_proj((2,2),(2,2))*catv_10
            @test catv_16 == Behavior.norm_to_gen_proj((2,2),(2,2))*catv_12
        end
    end

    @testset "verify that p16_vertices are normalized" begin
        p16_vertices = DichotomicLocalPolytope.deterministic_fixed_direction_signaling_behaviors(subspace="p16")

        for v in p16_vertices
            for i in 1:4
                @test 1 == v[i] + v[i+4] + v[i+8] + v[i+12]
            end
        end
    end
end

@testset "deterministic_one_way_signaling_behaviors()" begin
    vertices = DichotomicLocalPolytope.deterministic_one_way_signaling_behaviors()

    @test size(vertices) == (112,)

    p12_vertices = DichotomicLocalPolytope.deterministic_one_way_signaling_behaviors(subspace="p12")
    @test size(p12_vertices) == (112,)

    @testset "behaviors are normalized" begin
        for vertex in vertices
            @test size(vertex) == (16,)
            # there are 4 normalization constraints per behavior
            for i in 1:4
                @test 1 == vertex[i] + vertex[4+i] + vertex[8+i] + vertex[12+i]
            end
        end

        for vertex in p12_vertices
            @test size(vertex) == (12,)
            for i in 1:4
                @test 1 >= vertex[i] + vertex[4+i] + vertex[8+i] >= 0
            end
        end
    end

    @testset "behaviors are unique" begin
        for test_vertex in vertices
            duplicate_count = 0
            for match_vertex in vertices
                if test_vertex == match_vertex
                    duplicate_count += 1
                end
            end
            @test duplicate_count == 1
        end

        for test_vertex in p12_vertices
            duplicate_count = 0
            for match_vertex in p12_vertices
                if test_vertex == match_vertex
                    duplicate_count += 1
                end
            end
            @test duplicate_count == 1
        end
    end
end

@testset "deterministic_no_signaling_behaviors()" begin

    @testset "16 vertices are returned for bipartite case" begin
        vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors()
        @test size(vertices) == (16,)

        p12_vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors(subspace="p12")
        @test size(p12_vertices) == (16,)
    end

    @testset "vertices are normalized with respect to unique inputs" begin
        vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors()
        for vertex in vertices
            # there are 4 normalization constraints per behavior
            for i in 1:4
                @test 1 == vertex[i] + vertex[4+i] + vertex[8+i] + vertex[12+i]
            end
        end

        p12_vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors(subspace="p12")
        for vertex in p12_vertices
            for i in 1:4
                @test 1 >= vertex[i] + vertex[4+i] + vertex[8+i] >= 0
            end
        end
    end

    @testset "vertices are unique" begin
        vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors()

        for i in 1:size(vertices)[1]
            test_vertex = pop!(vertices)
            non_unique_count = 0
            for vertex in vertices
                if test_vertex == vertex
                    non_unique_count += 1
                end
            end
            @test non_unique_count == 0
        end

        p12_vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors(subspace="p12")

        for test_vertex in p12_vertices
            duplicate_count = 0
            for match_vertex in p12_vertices
                if test_vertex == match_vertex
                    duplicate_count += 1
                end
            end
            @test duplicate_count == 1
        end
    end

    @testset "p16 vertice match p12 vertices" begin
        p16_vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors(subspace="p16")
        p12_vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors(subspace="p12")

        p12_match = []
        for vertex in p16_vertices
            catv_16 = cat([1],vertex,dims=1)
            push!(p12_match, (Behavior.gen_to_norm_proj((2,2),(2,2))*catv_16)[2:end] )
        end
        @test p12_vertices == p12_match

        p16_match = []
        for vertex in p12_vertices
            catv_12 = cat([1],vertex,dims=1)
            push!(p16_match, (Behavior.norm_to_gen_proj((2,2),(2,2))*catv_12)[2:end])
        end
        @test p16_vertices == p16_match
    end

    @testset "vertices match the no-signaling polytope" begin
        vertices = DichotomicLocalPolytope.deterministic_no_signaling_behaviors()

        no_signaling_vertices_match = [
            [0,0,0,0,0,0,0,0],
            [0,0,0,1,0,0,0,0],
            [0,0,1,0,0,0,0,0],
            [0,0,1,1,0,0,0,0],
            [0,1,0,0,0,0,0,0],
            [0,1,0,1,0,0,0,1],
            [0,1,1,0,0,0,1,0],
            [0,1,1,1,0,0,1,1],
            [1,0,0,0,0,0,0,0],
            [1,0,0,1,0,1,0,0],
            [1,0,1,0,1,0,0,0],
            [1,0,1,1,1,1,0,0],
            [1,1,0,0,0,0,0,0],
            [1,1,0,1,0,1,0,1],
            [1,1,1,0,1,0,1,0],
            [1,1,1,1,1,1,1,1],
        ]

        # loop through each computed vertex and make sure that it is in the
        # match set only once.
        for vertex in vertices
            catv = cat([1],vertex,dims=1)
            no_signaling_vertex = (Behavior.gen_to_ns_proj((2,2),(2,2))*catv)[2:end]

            match_found = 0
            for match_vertex in no_signaling_vertices_match
                if no_signaling_vertex == match_vertex
                    match_found += 1
                end
            end
            @test match_found == 1
        end
    end
end

@testset "two_input_strategies()" begin
    strategies = DichotomicLocalPolytope.two_input_strategies()

    @test size(strategies) == (16,)

    @testset "two_input_strategies are unique amongst themselves" begin
        for test_strategy in strategies
            duplicate_count = 0
            for match_strategy in strategies
                if test_strategy == match_strategy
                    duplicate_count += 1
                end
            end
            # test_strategy should only match with itself
            @test duplicate_count == 1
        end
    end
end

@testset "dichotomic_strategies()" begin
    s = DichotomicLocalPolytope.dichotomic_strategies()
    @test s[1] == [1 0; 0 1]
    @test s[2] == [1 1; 0 0]
    @test s[3] == [0 0; 1 1]
    @test s[4] == [0 1; 1 0]
end

@testset "Permutations of communication protocols create redundant strategies" begin
    # consider a two-input strategy α(x,σ(y)) → α*(x⊗σ(y)) = α*(I⊗σ)*(x⊗y)
    # in the below tests verify that { α } = { α*(I⊗σ) } i.e. the set of
    # permuted strategies contains all members of {α} and that {α} contains the
    # set of all permuted strategies.

    two_input_strategies = DichotomicLocalPolytope.two_input_strategies()
    dichotomic_strategies = DichotomicLocalPolytope.dichotomic_strategies()

    permuted_strategies = []
    for ds in dichotomic_strategies[[2 3 4]]
        input_matrix = kron([1 0; 0 1], ds)
        for ts in two_input_strategies
            push!(permuted_strategies,ts*input_matrix)
        end
    end

    @testset "each permuted strategy appears once in the set of two input strategies" begin
        for ps_test in permuted_strategies
            duplicate_count = 0
            for ts_match in two_input_strategies
                if ps_test == ts_match
                    duplicate_count += 1
                end
            end

            @test duplicate_count == 1
        end
    end

    @testset "each two input strategy appears at least once in the set of permuted strategies" begin
        for ts_test in two_input_strategies
            duplicate_count = 0
            for ps_match in permuted_strategies
                if ts_test == ps_match
                    duplicate_count += 1
                end
            end
            @test duplicate_count > 0
        end
    end
end

end

ext: jl


text:
using UUIDs, Random

u1 = uuid1()
u4 = uuid4()
@test uuid_version(u1) == 1
@test uuid_version(u4) == 4
@test u1 == UUID(string(u1)) == UUID(GenericString(string(u1)))
@test u4 == UUID(string(u4)) == UUID(GenericString(string(u4)))
@test u1 == UUID(UInt128(u1))
@test u4 == UUID(UInt128(u4))
@test uuid4(MersenneTwister(0)) == uuid4(MersenneTwister(0))
@test_throws ArgumentError UUID("550e8400e29b-41d4-a716-446655440000")
@test_throws ArgumentError UUID("550e8400e29b-41d4-a716-44665544000098")
@test_throws ArgumentError UUID("z50e8400-e29b-41d4-a716-446655440000")

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
# LLVM IR generation

function irgen(@nospecialize(job::CompilerJob), method_instance::Core.MethodInstance;
               ctx::Context)
    mod, compiled = @timeit_debug to "emission" compile_method_instance(job, method_instance; ctx)
    if job.entry_abi === :specfunc
        entry_fn = compiled[method_instance].specfunc
    else
        entry_fn = compiled[method_instance].func
    end

    # clean up incompatibilities
    @timeit_debug to "clean-up" begin
        for llvmf in functions(mod)
            # only occurs in debug builds
            delete!(function_attributes(llvmf), EnumAttribute("sspstrong", 0; ctx))

            if Sys.iswindows()
                personality!(llvmf, nothing)
            end

            # remove the non-specialized jfptr functions
            # TODO: Do we need to remove these?
            if job.entry_abi === :specfunc
                if startswith(LLVM.name(llvmf), "jfptr_")
                    unsafe_delete!(mod, llvmf)
                end
            end
        end

        # remove the exception-handling personality function
        if Sys.iswindows() && "__julia_personality" in functions(mod)
            llvmf = functions(mod)["__julia_personality"]
            @compiler_assert isempty(uses(llvmf)) job
            unsafe_delete!(mod, llvmf)
        end
    end

    # target-specific processing
    process_module!(job, mod)
    entry = functions(mod)[entry_fn]

    # sanitize function names
    # FIXME: Julia should do this, but apparently fails (see maleadt/LLVM.jl#201)
    for f in functions(mod)
        LLVM.isintrinsic(f) && continue
        llvmfn = LLVM.name(f)
        startswith(llvmfn, "julia.") && continue # Julia intrinsics
        startswith(llvmfn, "llvm.") && continue # unofficial LLVM intrinsics
        llvmfn′ = safe_name(llvmfn)
        if llvmfn != llvmfn′
            @assert !haskey(functions(mod), llvmfn′)
            LLVM.name!(f, llvmfn′)
        end
    end

    # rename and process the entry point
    if job.source.name !== nothing
        LLVM.name!(entry, safe_name(string("julia_", job.source.name)))
    end
    if job.source.kernel
        LLVM.name!(entry, mangle_call(entry, job.source.tt))
    end
    entry = process_entry!(job, mod, entry)
    if job.entry_abi === :specfunc
        func = compiled[method_instance].func
        specfunc = LLVM.name(entry)
    else
        func = LLVM.name(entry)
        specfunc = compiled[method_instance].specfunc
    end

    compiled[method_instance] =
        (; compiled[method_instance].ci, func, specfunc)

    # minimal required optimization
    @timeit_debug to "rewrite" ModulePassManager() do pm
        global current_job
        current_job = job

        linkage!(entry, LLVM.API.LLVMExternalLinkage)

        # internalize all functions, but keep exported global variables
        exports = String[LLVM.name(entry)]
        for gvar in globals(mod)
            push!(exports, LLVM.name(gvar))
        end
        internalize!(pm, exports)

        # inline llvmcall bodies
        always_inliner!(pm)

        can_throw(job) || add!(pm, ModulePass("LowerThrow", lower_throw!))

        add_lowering_passes!(job, pm)

        run!(pm, mod)
    end

    return mod, compiled
end

## name mangling

# we generate function names that look like C++ functions, because many NVIDIA tools
# support them, e.g., grouping different instantiations of the same kernel together.

function mangle_param(t, substitutions)
    t == Nothing && return "v"

    if isa(t, DataType) && t <: Ptr
        tn = mangle_param(eltype(t), substitutions)
        "P$tn"
    elseif isa(t, DataType) || isa(t, Core.Function)
        tn = safe_name(t)

        # handle substitutions
        sub = findfirst(isequal(tn), substitutions)
        if sub === nothing
            str = "$(length(tn))$tn"
            push!(substitutions, tn)
        elseif sub == 1
            str = "S_"
        else
            str = "S$(sub-2)_"
        end

        # encode typevars as template parameters
        if !isempty(t.parameters)
            str *= "I"
            for t in t.parameters
                str *= mangle_param(t, substitutions)
            end
            str *= "E"
        end

        str
    elseif isa(t, Integer)
        t > 0 ? "Li$(t)E" : "Lin$(abs(t))E"
    else
        tn = safe_name(t)   # TODO: actually does support digits...
        if startswith(tn, r"\d")
            # C++ classes cannot start with a digit, so mangling doesn't support it
            tn = "_$(tn)"
        end
        "$(length(tn))$tn"
    end
end

function mangle_call(f, tt)
    fn = safe_name(f)
    str = "_Z$(length(fn))$fn"

    substitutions = String[]
    for t in tt.parameters
        str *= mangle_param(t, substitutions)
    end

    return str
end

# make names safe for ptxas
safe_name(fn::String) = replace(fn, r"[^A-Za-z0-9_]"=>"_")
safe_name(f::Union{Core.Function,DataType}) = safe_name(String(nameof(f)))
safe_name(f::LLVM.Function) = safe_name(LLVM.name(f))
safe_name(x) = safe_name(repr(x))

## exception handling

# this pass lowers `jl_throw` and friends to GPU-compatible exceptions.
# this isn't strictly necessary, but has a couple of advantages:
# - we can kill off unused exception arguments that otherwise would allocate or invoke
# - we can fake debug information (lacking a stack unwinder)
#
# once we have thorough inference (ie. discarding `@nospecialize` and thus supporting
# exception arguments) and proper debug info to unwind the stack, this pass can go.
function lower_throw!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "lower throw" begin

    throw_functions = [
        # unsupported runtime functions that are used to throw specific exceptions
        "jl_throw"                      => "exception",
        "jl_error"                      => "error",
        "jl_too_few_args"               => "too few arguments exception",
        "jl_too_many_args"              => "too many arguments exception",
        "jl_type_error"                 => "type error",
        "jl_type_error_rt"              => "type error",
        "jl_undefined_var_error"        => "undefined variable error",
        "jl_bounds_error"               => "bounds error",
        "jl_bounds_error_v"             => "bounds error",
        "jl_bounds_error_int"           => "bounds error",
        "jl_bounds_error_tuple_int"     => "bounds error",
        "jl_bounds_error_unboxed_int"   => "bounds error",
        "jl_bounds_error_ints"          => "bounds error",
        "jl_eof_error"                  => "EOF error",
    ]

    for f in functions(mod)
        fn = LLVM.name(f)
        for (throw_fn, name) in throw_functions
            occursin(throw_fn, fn) || continue

            for use in uses(f)
                call = user(use)::LLVM.CallInst

                # replace the throw with a PTX-compatible exception
                let builder = Builder(ctx)
                    position!(builder, call)
                    emit_exception!(builder, name, call)
                    dispose(builder)
                end

                # remove the call
                call_args = arguments(call)
                unsafe_delete!(LLVM.parent(call), call)

                # HACK: kill the exceptions' unused arguments
                #       this is needed for throwing objects with @nospecialize constructors.
                for arg in call_args
                    # peek through casts
                    if isa(arg, LLVM.AddrSpaceCastInst)
                        cast = arg
                        arg = first(operands(cast))
                        isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                    end

                    if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                        unsafe_delete!(LLVM.parent(arg), arg)
                    end
                end

                changed = true
            end

            @compiler_assert isempty(uses(f)) job
            break
         end
     end

    end
    return changed
end

# report an exception in a GPU-compatible manner
#
# the exact behavior depends on the debug level. in all cases, a `trap` will be emitted, On
# debug level 1, the exception name will be printed, and on debug level 2 the individual
# stack frames (as recovered from the LLVM debug information) will be printed as well.
function emit_exception!(builder, name, inst)
    job = current_job::CompilerJob
    bb = position(builder)
    fun = LLVM.parent(bb)
    mod = LLVM.parent(fun)
    ctx = context(mod)

    # report the exception
    if Base.JLOptions().debug_level >= 1
        name = globalstring_ptr!(builder, name, "exception")
        if Base.JLOptions().debug_level == 1
            call!(builder, Runtime.get(:report_exception), [name])
        else
            call!(builder, Runtime.get(:report_exception_name), [name])
        end
    end

    # report each frame
    if Base.JLOptions().debug_level >= 2
        rt = Runtime.get(:report_exception_frame)
        ft = convert(LLVM.FunctionType, rt; ctx)
        bt = backtrace(inst)
        for (i,frame) in enumerate(bt)
            idx = ConstantInt(parameters(ft)[1], i)
            func = globalstring_ptr!(builder, String(frame.func), "di_func")
            file = globalstring_ptr!(builder, String(frame.file), "di_file")
            line = ConstantInt(parameters(ft)[4], frame.line)
            call!(builder, rt, [idx, func, file, line])
        end
    end

    # signal the exception
    call!(builder, Runtime.get(:signal_exception))

    emit_trap!(job, builder, mod, inst)
end

function emit_trap!(@nospecialize(job::CompilerJob), builder, mod, inst)
    ctx = context(mod)
    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    call!(builder, trap)
end

## kernel promotion

@enum ArgumentCC begin
    BITS_VALUE  # bitstype, passed as value
    BITS_REF    # bitstype, passed as pointer
    MUT_REF     # jl_value_t*, or the anonymous equivalent
    GHOST       # not passed
end

function classify_arguments(@nospecialize(job::CompilerJob), codegen_ft::LLVM.FunctionType)
    source_sig = typed_signature(job)

    source_types = [source_sig.parameters...]

    codegen_types = parameters(codegen_ft)

    args = []
    codegen_i = 1
    for (source_i, source_typ) in enumerate(source_types)
        if isghosttype(source_typ) || Core.Compiler.isconstType(source_typ)
            push!(args, (cc=GHOST, typ=source_typ))
            continue
        end

        codegen_typ = codegen_types[codegen_i]
        if codegen_typ isa LLVM.PointerType && !issized(eltype(codegen_typ))
            push!(args, (cc=MUT_REF, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        elseif codegen_typ isa LLVM.PointerType && issized(eltype(codegen_typ)) &&
               !(source_typ <: Ptr) && !(source_typ <: Core.LLVMPtr)
            push!(args, (cc=BITS_REF, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        else
            push!(args, (cc=BITS_VALUE, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        end
        codegen_i += 1
    end

    return args
end

if VERSION >= v"1.7.0-DEV.204"
    function is_immutable_datatype(T::Type)
        isa(T,DataType) && !Base.ismutabletype(T)
    end
else
    function is_immutable_datatype(T::Type)
        isa(T,DataType) && !T.mutable
    end
end

if VERSION >= v"1.7.0-DEV.204"
    function is_inlinealloc(T::Type)
        mayinlinealloc = (T.name.flags >> 2) & 1 == true
        # FIXME: To simple
        return mayinlinealloc
    end
else
    function is_inlinealloc(T::Type)
        return T.isinlinealloc
    end
end

function is_concrete_immutable(T::Type)
    is_immutable_datatype(T) && T.layout !== C_NULL
end

function is_pointerfree(T::Type)
    if !is_immutable_datatype(T)
        return false
    end
    return Base.datatype_pointerfree(T)
end

function deserves_stack(@nospecialize(T))
    if !is_concrete_immutable(T)
        return false
    end
    return is_inlinealloc(T)
end

deserves_argbox(T) = !deserves_stack(T)
deserves_retbox(T) = deserves_argbox(T)
function deserves_sret(T, llvmT)
    @assert isa(T,DataType)
    sizeof(T) > sizeof(Ptr{Cvoid}) && !isa(llvmT, LLVM.FloatingPointType) && !isa(llvmT, LLVM.VectorType)
end

# byval lowering
#
# some back-ends don't support byval, or support it badly, so lower it eagerly ourselves
# https://reviews.llvm.org/D79744
function lower_byval(@nospecialize(job::CompilerJob), mod::LLVM.Module, f::LLVM.Function)
    ctx = context(mod)
    ft = eltype(llvmtype(f))
    @compiler_assert LLVM.return_type(ft) == LLVM.VoidType(ctx) job

    # find the byval parameters
    byval = BitVector(undef, length(parameters(ft)))
    if LLVM.version() >= v"12"
        for i in 1:length(byval)
            attrs = collect(parameter_attributes(f, i))
            byval[i] = any(attrs) do attr
                kind(attr) == kind(EnumAttribute("byval", 0; ctx))
            end
        end
    else
        # XXX: byval is not round-trippable on LLVM < 12 (see maleadt/LLVM.jl#186)
        #      so we need to re-classify the Julia arguments.
        #      remove this once we only support 1.7.
        args = classify_arguments(job, ft)
        filter!(args) do arg
            arg.cc != GHOST
        end
        for arg in args
            if arg.cc == BITS_REF
                # NOTE: +1 since this pass runs after introducing the kernel state
                byval[arg.codegen.i] = true
            end
        end
    end

    # fixup metadata
    #
    # Julia emits invariant.load and const TBAA metadta on loads from pointer args,
    # which is invalid now that we have materialized the byval.
    for (i, param) in enumerate(parameters(f))
        if byval[i]
            # collect all uses of the argument
            worklist = Vector{LLVM.Instruction}(user.(collect(uses(param))))
            while !isempty(worklist)
                value = popfirst!(worklist)

                # remove the invariant.load attribute
                md = metadata(value)
                if haskey(md, LLVM.MD_invariant_load)
                    delete!(md, LLVM.MD_invariant_load)
                end
                if haskey(md, LLVM.MD_tbaa)
                    delete!(md, LLVM.MD_tbaa)
                end

                # recurse on the output of some instructions
                if isa(value, LLVM.BitCastInst) ||
                   isa(value, LLVM.GetElementPtrInst) ||
                   isa(value, LLVM.AddrSpaceCastInst)
                    append!(worklist, user.(collect(uses(value))))
                end
            end
        end
    end

    # generate the new function type & definition
    new_types = LLVM.LLVMType[]
    for (i, param) in enumerate(parameters(ft))
        if byval[i]
            push!(new_types, eltype(param::LLVM.PointerType))
        else
            push!(new_types, param)
        end
    end
    new_ft = LLVM.FunctionType(LLVM.return_type(ft), new_types)
    new_f = LLVM.Function(mod, "", new_ft)
    linkage!(new_f, linkage(f))
    for (arg, new_arg) in zip(parameters(f), parameters(new_f))
        LLVM.name!(new_arg, LLVM.name(arg))
    end

    # emit IR performing the "conversions"
    new_args = LLVM.Value[]
    Builder(ctx) do builder
        entry = BasicBlock(new_f, "entry"; ctx)
        position!(builder, entry)

        # perform argument conversions
        for (i, param) in enumerate(parameters(ft))
            if byval[i]
                # copy the argument value to a stack slot, and reference it.
                ptr = alloca!(builder, eltype(param))
                if LLVM.addrspace(param) != 0
                    ptr = addrspacecast!(builder, ptr, param)
                end
                store!(builder, parameters(new_f)[i], ptr)
                push!(new_args, ptr)
            else
                push!(new_args, parameters(new_f)[i])
                for attr in collect(parameter_attributes(f, i))
                    push!(parameter_attributes(new_f, i), attr)
                end
            end
        end

        # map the arguments
        value_map = Dict{LLVM.Value, LLVM.Value}(
            param => new_args[i] for (i,param) in enumerate(parameters(f))
        )

        value_map[f] = new_f
        clone_into!(new_f, f; value_map,
                    changes=LLVM.API.LLVMCloneFunctionChangeTypeGlobalChanges)

        # fall through
        br!(builder, blocks(new_f)[2])
    end

    # remove the old function
    # NOTE: if we ever have legitimate uses of the old function, create a shim instead
    fn = LLVM.name(f)
    @assert isempty(uses(f))
    replace_metadata_uses!(f, new_f)
    unsafe_delete!(mod, f)
    LLVM.name!(new_f, fn)

    return new_f
end

# kernel state arguments
#
# to facilitate passing stateful information to kernels without having to recompile, e.g.,
# the storage location for exception flags, or the location of a I/O buffer, we enable the
# back-end to specify a Julia object that will be passed to the kernel by-value, and to
# every called function by-reference. Access to this object is done using the
# `julia.gpu.state_getter` intrinsic. after optimization, these intrinsics will be lowered
# to refer to the state argument.
#
# note that we deviate from the typical Julia calling convention, by always passing the
# state objects by value instead of by reference, this to ensure that the state object
# is not copied to the stack (because LLVM doesn't see that all uses are read-only).
# in principle, `readonly byval` should be equivalent, but LLVM doesn't realize that.
# also see https://github.com/JuliaGPU/CUDA.jl/pull/1167 and the comments in that PR.
# once LLVM supports this pattern, consider going back to passing the state by reference,
# so that the julia.gpu.state_getter` can be simplified to return an opaque pointer.

# add a state argument to every function in the module, starting from the kernel entry point
function add_kernel_state!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)

    # check if we even need a kernel state argument
    state = kernel_state_type(job)
    @assert job.source.kernel
    if state === Nothing
        return false
    end
    T_state = convert(LLVMType, state; ctx)

    # intrinsic returning an opaque pointer to the kernel state.
    # this is both for extern uses, and to make this transformation a two-step process.
    state_intr = kernel_state_intr(mod, T_state)

    kernels = []
    kernels_md = metadata(mod)["julia.kernel"]
    for kernel_md in operands(kernels_md)
        push!(kernels, Value(operands(kernel_md)[1]; ctx))
    end

    # determine which functions need a kernel state argument
    #
    # previously, we add the argument to every function and relied on unused arg elim to
    # clean-up the IR. however, some libraries do Funny Stuff, e.g., libdevice bitcasting
    # function pointers. such IR is hard to rewrite, so instead be more conservative.
    worklist = Set{LLVM.Function}([state_intr, kernels...])
    worklist_length = 0
    while worklist_length != length(worklist)
        # iteratively discover functions that use the intrinsic or any function calling it
        worklist_length = length(worklist)
        additions = LLVM.Function[]
        function check_user(val)
            if val isa Instruction
                bb = LLVM.parent(val)
                new_f = LLVM.parent(bb)
                in(new_f, worklist) || push!(additions, new_f)
            elseif val isa ConstantExpr
                # constant expressions don't have a parent; we need to look up their uses
                for use in uses(val)
                    check_user(user(use))
                end
            else
                error("Don't know how to check uses of $val. Please file an issue.")
            end
        end
        for f in worklist, use in uses(f)
            check_user(user(use))
        end
        for f in additions
            push!(worklist, f)
        end
    end
    delete!(worklist, state_intr)

    # add a state argument
    workmap = Dict{LLVM.Function, LLVM.Function}()
    for f in worklist
        fn = LLVM.name(f)
        ft = eltype(llvmtype(f))
        LLVM.name!(f, fn * ".stateless")

        # create a new function
        new_param_types = [T_state, parameters(ft)...]
        new_ft = LLVM.FunctionType(LLVM.return_type(ft), new_param_types)
        new_f = LLVM.Function(mod, fn, new_ft)
        LLVM.name!(parameters(new_f)[1], "state")
        linkage!(new_f, linkage(f))
        for (arg, new_arg) in zip(parameters(f), parameters(new_f)[2:end])
            LLVM.name!(new_arg, LLVM.name(arg))
        end

        workmap[f] = new_f
    end

    # clone and rewrite the function bodies, replacing uses of the old stateless function
    # with the newly created definition that includes the state argument.
    #
    # most uses are rewritten by LLVM by putting the functions in the value map.
    # a separate value materializer is used to recreate constant expressions.
    #
    # note that this only _replaces_ the uses of these functions, we'll still need to
    # _correct_ the uses (i.e. actually add the state argument) afterwards.
    function materializer(val)
        if val isa ConstantExpr
            if opcode(val) == LLVM.API.LLVMBitCast
                target = operands(val)[1]
                if target isa LLVM.Function && haskey(workmap, target)
                    # the function is being bitcasted to a different function type.
                    # we need to mutate that function type to include the state argument,
                    # or we'd be invoking the original function in an invalid way.
                    #
                    # XXX: ptrtoint/inttoptr pairs can also lose the state argument...
                    #      is all this even sound?
                    typ = llvmtype(val)::LLVM.PointerType
                    ft = eltype(typ)::LLVM.FunctionType
                    new_ft = LLVM.FunctionType(LLVM.return_type(ft), [T_state, parameters(ft)...])
                    return const_bitcast(workmap[target], LLVM.PointerType(new_ft, addrspace(typ)))
                end
            elseif opcode(val) == LLVM.API.LLVMPtrToInt
                target = operands(val)[1]
                if target isa LLVM.Function && haskey(workmap, target)
                    return const_ptrtoint(workmap[target], llvmtype(val))
                end
            end
        end
        return val
    end
    for (f, new_f) in workmap
        # use a value mapper for rewriting function arguments
        value_map = Dict{LLVM.Value, LLVM.Value}()
        for (param, new_param) in zip(parameters(f), parameters(new_f)[2:end])
            LLVM.name!(new_param, LLVM.name(param))
            value_map[param] = new_param
        end

        # rewrite references to the old function
        merge!(value_map, workmap)

        clone_into!(new_f, f; value_map, materializer,
                    changes=LLVM.API.LLVMCloneFunctionChangeTypeGlobalChanges)

        # remove the function IR so that we won't have any uses left after this pass.
        empty!(f)
    end

    # ensure the old (stateless) functions don't have uses anymore, and remove them
    for f in keys(workmap)
        for use in uses(f)
            val = user(use)
            if val isa ConstantExpr
                # XXX: shouldn't clone_into! remove unused CEs?
                isempty(uses(val)) || error("old function still has uses (via a constant expr)")
                LLVM.unsafe_destroy!(val)
            else
                error("old function still has uses")
            end
        end
        replace_metadata_uses!(f, workmap[f])
        unsafe_delete!(mod, f)
    end

    # update uses of the new function, modifying call sites to include the kernel state
    function rewrite_uses!(f)
        # update uses
        Builder(ctx) do builder
            for use in uses(f)
                val = user(use)
                if val isa LLVM.CallBase && called_value(val) == f
                    # NOTE: we don't rewrite calls using Julia's jlcall calling convention,
                    #       as those have a fixed argument list, passing actual arguments
                    #       in an array of objects. that doesn't matter, for now, since
                    #       GPU back-ends don't support such calls anyhow. but if we ever
                    #       want to support kernel state passing on more capable back-ends,
                    #       we'll need to update the argument array instead.
                    if callconv(val) == 37 || callconv(val) == 38
                        # TODO: update for LLVM 15 when JuliaLang/julia#45088 is merged.
                        continue
                    end

                    # forward the state argument
                    position!(builder, val)
                    state = call!(builder, state_intr, Value[], "state")
                    new_val = if val isa LLVM.CallInst
                        call!(builder, f, [state, arguments(val)...], operand_bundles(val))
                    else
                        # TODO: invoke and callbr
                        error("Rewrite of $(typeof(val))-based calls is not implemented: $val")
                    end
                    callconv!(new_val, callconv(val))

                    replace_uses!(val, new_val)
                    @assert isempty(uses(val))
                    unsafe_delete!(LLVM.parent(val), val)
                elseif val isa LLVM.CallBase
                    # the function is being passed as an argument, which we'll just permit,
                    # because we expect to have rewritten the call down the line separately.
                elseif val isa LLVM.StoreInst
                    # the function is being stored, which again we'll permit like before.
                elseif val isa ConstantExpr
                    rewrite_uses!(val)
                else
                    error("Cannot rewrite $(typeof(val)) use of function: $val")
                end
            end
        end
    end
    for f in values(workmap)
        rewrite_uses!(f)
    end

    return true
end

# lower calls to the state getter intrinsic. this is a two-step process, so that the state
# argument can be added before optimization, and that optimization can introduce new uses
# before the intrinsic getting lowered late during optimization.
function lower_kernel_state!(fun::LLVM.Function)
    job = current_job::CompilerJob
    mod = LLVM.parent(fun)
    ctx = context(fun)
    changed = false

    # check if we even need a kernel state argument
    state = kernel_state_type(job)
    if state === Nothing
        return false
    end

    # fixup all uses of the state getter to use the newly introduced function state argument
    if haskey(functions(mod), "julia.gpu.state_getter")
        state_intr = functions(mod)["julia.gpu.state_getter"]
        state_arg = nothing # only look-up when needed

        Builder(ctx) do builder
            for use in uses(state_intr)
                inst = user(use)
                @assert inst isa LLVM.CallInst
                bb = LLVM.parent(inst)
                LLVM.parent(bb) == fun || continue

                position!(builder, inst)
                bb = LLVM.parent(inst)
                f = LLVM.parent(bb)

                if state_arg === nothing
                    # find the kernel state argument. this should be the first argument of
                    # the function, but only when this function needs the state!
                    state_arg = parameters(fun)[1]
                    T_state = convert(LLVMType, state; ctx)
                    @assert llvmtype(state_arg) == T_state
                end

                replace_uses!(inst, state_arg)

                @assert isempty(uses(inst))
                unsafe_delete!(LLVM.parent(inst), inst)

                changed = true
            end
        end
    end

    return changed
end

function cleanup_kernel_state!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false

    # remove the getter intrinsic
    if haskey(functions(mod), "julia.gpu.state_getter")
        intr = functions(mod)["julia.gpu.state_getter"]
        if isempty(uses(intr))
            # if we're not emitting a kernel, we can't resolve the intrinsic to an argument.
            unsafe_delete!(mod, intr)
            changed = true
        end
    end

    return changed
end

function kernel_state_intr(mod::LLVM.Module, T_state)
    ctx = context(mod)

    state_intr = if haskey(functions(mod), "julia.gpu.state_getter")
        functions(mod)["julia.gpu.state_getter"]
    else
        LLVM.Function(mod, "julia.gpu.state_getter", LLVM.FunctionType(T_state))
    end
    push!(function_attributes(state_intr), EnumAttribute("readnone", 0; ctx))

    return state_intr
end

# run-time equivalent
function kernel_state_value(state)
    Context() do ctx
        T_state = convert(LLVMType, state; ctx)

        # create function
        llvm_f, _ = create_function(T_state)
        mod = LLVM.parent(llvm_f)

        # get intrinsic
        state_intr = kernel_state_intr(mod, T_state)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            val = call!(builder, state_intr, Value[], "state")

            ret!(builder, val)
        end

        call_function(llvm_f, state)
    end
end


text:
using BHAtp, CSV, Profile
const TP = BHAtp
const profile = false

ProjDir = dirname(@__FILE__)
ProjName = split(ProjDir, "/")[end]
bhaj = BHAJ(ProjName, ProjDir)
bhaj.ratio = 1.7

segs = [
# Element type,  Material,    Length,     OD,         ID,    fc
  (:bit,                      :steel,        0.00,   2.75,        9.0,  0.0),
  (:collar,                  :steel,     45.00,   2.75,        7.0,  0.0),
  (:stabilizer,            :steel,        0.00,   2.75,       9.0,  0.0),
  (:pipe,                    :steel,     30.00,   2.75,       7.0,  0.0),
  (:stabilizer,             :steel,       0.00,   2.75,       9.0,  0.0),
  (:pipe,                     :steel,    90.00,   2.75,       7.0,  0.0),
  (:stabilizer,             :steel,       0.00,   2.75,       9.0,  0.0),
  (:pipe,                    :steel,   100.00,   2.75,       7.0,  0.0)
]

traj = [
#   Heading,      Diameter
  ( 60.0,      9.0)
]

wobs = 20:10:40
incls = 20:10:40             # Or e.g. incls = [5 10 20 30 40 45 50]

if profile
  Profile.clear()
  Profile.@profile bhaj(segs, traj, wobs, incls)
  Profile.print()
  println()
else
  @time bhaj(segs, traj, wobs, incls)
  println()
end

display("Fetch tp=false, wob=40, incl@bit=20 solution:")
df,df_tp = show_solution(ProjDir, 40, 20, show=false, tp=false)
df[:,[2; 5:6; 9:12]] |> display
println()

println("Fetch tp=true, wob=40, incl@bit=20 solution:")
df,df_tp = show_solution(ProjDir, 40, 20, show=false);
df[:,[2; 5:6; 9:12]] |> display
println()

show_tp(ProjDir, wobs, incls)


text:
"""
    LinearRegression()
Classic regression model. This struct has no parameter.
If you want to use polynomial model, use `Regression.make_design_matrix()`.

see also: [`make_design_matrix`](@ref)

# Example
```jldoctest regression
julia> x = [
    16.862463771320925 68.10823385851712
    15.382965696961577 65.4313485700859
    8.916228406218375 53.92034559524475
    10.560285659132695 59.17305391117168
    12.142253214135884 62.28708207525656
    5.362107221163482 43.604947901567414
    13.893239446341777 62.44348617377496
    11.871357065173395 60.28433066289655
    29.83792267802442 69.22281924803998
    21.327107214235483 70.15810991597944
    23.852372696012498 69.81780163668844
    26.269031430914108 67.61037566099782
    22.78907104644012 67.78105545358633
    26.73342178134947 68.59263965946904
    9.107259141706415 56.565383817343495
    29.38551885863976 68.1005579469209
    7.935966787763017 53.76264777936664
    29.01677894379809 68.69484161138638
    6.839609488194577 49.69794758177567
    13.95215840314148 62.058116579899085]; #These data are also used to explanations of other functions.

julia> t = [169.80980778351542, 167.9081124078835, 152.30845618985222, 160.3110300206261, 161.96826472170756, 136.02842285615077, 163.98131131382686, 160.117817321485, 172.22758529098235, 172.21342437006865, 171.8939175591617, 169.83018083884602, 171.3878062674257, 170.52487535026015, 156.40282783981309, 170.6488327896672, 151.69267899906185, 172.32478221316322, 145.14365314788827, 163.79383292080666];

julia> model = LinearRegression()
LinearRegression(Float64[])

julia> fit!(model, x, t)
3-element Vector{Float64}:
 -0.04772448076255398
  1.395963968616736
 76.7817095600793

julia> model(x)
20-element Vector{Float64}:
 171.05359766482795
 167.38737053144575
 151.62704681535598
 158.88117658330424
 163.15274911747872
 137.3967419011542
 163.28751869479999
 160.3699086857777
 171.99027166957023
 173.70207799107243
 173.10650291105486
 169.9096820022986
 170.31402414534642
 171.25872436348817
 155.31030802635905
 170.44522606721017
 151.45368882321284
 171.29242257091374
 145.83183688699864
 162.74674475052848
```
"""
mutable struct LinearRegression
    w::Array
    LinearRegression() = new(Array{Float64}(undef, 0))
end

"""
    fit!(model, x, t)
`x` must be the number of features in the first dimension and the second dimension must be the number of data.
"""
function fit!(model::LinearRegression, x, t)
    check_size(x, t)
    x = expand(x)
    model.w = inv(x * x') * x * t
end

(model::LinearRegression)(x) = expand(x)' * model.w


text:
# Test of the infeasibility detection feature of the COSMO solver on problems from the SDPLibrary
# SDP problems are given in the following format
# problem: max  tr( F0*Y)
#          s.t. tr(Fi Y) = ci, for i=1,...,m
#               Y ⪴ 0

# RESULTS SO FAR:
# infp1 primal form --> detected as prim inf, dual form --> detected as dual inf
# infp2 primal form --> detected as prim inf, dual form --> detected as dual inf
# infd1 primal form --> not detected, dual form --> not detected (both forms dual inf according to MOSEK)
# infd2 primal form --> not detected, dual form --> not detected (both forms dual inf according to MOSEK)

# workspace()
# include("../../../src/Solver.jl")
# using JLD
# using Base.Test
# using COSMO
#  using JuMP, Mosek
using JLD

problems = ["infp1.jld";"infp2.jld";"infd1.jld";"infd2.jld"]

@testset "Infeasible problems from SDP Library" begin
  for kkk = 1:length(problems)
    data = JLD.load("../sdplib/"*problems[kkk])
    F = data["F"]
    c = data["c"]
    m = data["m"]
    n = data["n"]
    optVal = data["optVal"]
    kkk <= 2 ? (problem_type = :Primal_infeasible) : (problem_type=:Dual_infeasible)

    # Rewrite problem to COSMO compatible format:
    # min   1/2 x'Px + q'x
    # s.t.  Ax + s = b
    #       s in K
    # -----------------------
    # primal form
    # -----------------------
    P = zeros(m,m)
    q = vec(c)
    A = zeros(n^2,m)
    for iii = 1:m
      A[:,iii] = -vec(F[iii+1])
    end
    b = -vec(F[1])
    Kf = 0
    Kl = 0
    Kq = []
    Ks = [n^2]
    # model = Model(solver=MosekSolver())
    # @variable(model, x[1:m])
    # @variable(model, S[1:n,1:n],SDP)
    # @objective(model, Min, q'*x)
    # s = vec(S)
    # @constraint(model, A*x+s .== b)
    # status = JuMP.solve(model)

    # -----------------------
    # dual form
    # -----------------------
    # P = zeros(n^2,n^2)
    # q = -vec(F[1])
    # A = zeros(m,n^2)
    # for iii = 1:m
    #   A[iii,:] = vec(F[iii+1])'
    # end
    # A = [A; -speye(n^2)]
    # b = [c;zeros(n^2)]
    # Kf = m
    # Kl = 0
    # Kq = []
    # Ks = [n^2]

    # model = Model(solver=MosekSolver())
    # @variable(model, x[1:n^2])
    # @variable(model, s1[1:m])
    # @variable(model, S2[1:n,1:n],SDP)
    # s = [s1;vec(S2)]
    # @objective(model, Min, q'*x)
    # @constraint(model, A*x+s .== b)
    # status = JuMP.solve(model)
    K = Cone(Kf,Kl,Kq,Ks)
    settings = Settings(max_iter=3000,verbose=false,check_termination=1,checkInfeasibility=50,scaling = 10 ,scaleFunc=2,adaptive_rho=true,eps_abs=1e-4,eps_rel=1e-4)
    res = COSMO.solve(P,q,A,b,K,settings);
    @test res.status == problem_type
  end
end


text:
# # Algebraic networks
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/experimental/algebraic_nets.ipynb)
#

using Catlab, Catlab.Doctrines, Catlab.Graphics
using Catlab.Experimental.AlgebraicNets
using Catlab.Experimental.MathFormulas

import TikzPictures
#-
R = Ob(AlgebraicNet, "\\mathbb{R}")
f_sin = Hom(:sin, R, R)
f_cos = Hom(:cos, R, R)
display(f_sin)
display(f_cos)
#-
f = compose(mcopy(R),otimes(f_sin,f_cos),mmerge(R))
#-
to_tikz(f, labels=true)
#-
compile_expr(f; args=[:x])
#-
to_formula(f, [:x])
#-
f = compose(mcopy(R,3), otimes(f_sin, f_cos, f_sin))
#-
to_tikz(f)
#-
compile_expr(f; args=[:x])
#-
f = compose(linear(2,R,R), f_sin, linear(2,R,R))
#-
to_tikz(f)
#-
compile_expr(f; args=[:x])
#-
to_formula(f, [:x])
#-
f = compose(mcopy(R), otimes(id(R),Hom(:cos,R,R)), Hom(:*,otimes(R,R),R))
#-
to_tikz(f)
#-
compile_expr(f; args=[:x])
#-
to_formula(f, [:x])


text:
using NESSie.Rjasanow
using NESSie.Rjasanow: laplacepot

@test_skip laplacecoll!

@testset "laplacepot" begin
    for T in testtypes
        # 1. Single layer potentials
        # 1.1 ξ is a triangle vertex
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        ξ = elem.v1     # h = √2, φ1 = -π/4, φ2=π/4
        res = laplacepot(SingleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.099189377627951192068415700486459476011275034113019436460)
        ξ = elem.v2     # h = 1, φ1 = -π/4, φ2 = 0
        res = laplacepot(SingleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.070137481542397515891964948616904961601109053078512490737)
        ξ = elem.v3     # h = 1, φ1 = 0, φ2 = π/4
        res = laplacepot(SingleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.070137481542397515891964948616904961601109053078512490737)

        # 1.2 ξ lies in the same plane as the triangle
        ξ = T[0, -1, 0]
        res = laplacepot(SingleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.029051896085553676176450751869554514410165981034506945723)
        ξ = T[0, -1, -1]
        res = laplacepot(SingleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ T(4π * (0.110553267370012295141653085597608171948437903548411687443
                                - 2 * 0.044743379406364904901413439333009046862112643066826198636))

        # 1.3 ξ lies on a line through a vertex, perpendicular
        #     to the triangle plane. In the following code, ξ
        #     is already proceted onto the plane!
        ξ = elem.v1
        res = laplacepot(SingleLayer, ξ, elem, one(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.034775298546314300210123433855429220291493953452986787115)
        ξ = elem.v2
        res = laplacepot(SingleLayer, ξ, elem, one(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.03156675645153950833618064292205638527751593336549)
        ξ = elem.v3
        res = laplacepot(SingleLayer, ξ, elem, one(T))
        @test typeof(res) == T
        @test res ≈ T(4π * 0.03156675645153950833618064292205638527751593336549)

        # 1.4 ξ lies somewhere else in space. Again, ξ is
        #     already projected onto the plane!
        ξ = T[0, -1, 0]
        res = laplacepot(SingleLayer, ξ, elem, one(T))
        @test typeof(res) == T
        @test res ≈ T(4π * (0.05486307812003175845441110089494128884416545308420259271
                                     -0.03156675645153950833618064292205638527751593336549))
        ξ = T[0, -1, -1]
        res = laplacepot(SingleLayer, ξ, elem, one(T))
        @test typeof(res) == T
        @test res ≈ T(4π * (2 * 0.035315594609376722705763636889615906447125307241104235442
                                    - 2 * 0.0260303288327139057352736693314968773933071931295084))

        # 2. Double layer potentials
        # 1.1 ξ is a triangle vertex
        ξ = elem.v1
        res = laplacepot(DoubleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ zero(T)
        ξ = elem.v2
        res = laplacepot(DoubleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ zero(T)
        ξ = elem.v3
        res = laplacepot(DoubleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ zero(T)

        # 1.2 ξ lies in the same plane as the triangle
        ξ = T[0, -1, 0]
        res = laplacepot(DoubleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ zero(T)
        ξ = T[0, -1, -1]
        res = laplacepot(DoubleLayer, ξ, elem, zero(T))
        @test typeof(res) == T
        @test res ≈ zero(T)

        @test_skip dl_ξ_above_vertex
        @test_skip dl_ξ_somewhere_in_space
    end
end


text:
using YaoBase, YaoArrayRegister, Random
export Measure, AllLocs, ComputationalBasis, chmeasureoperator

"""
    Measure{N, K, OT, RNG, IT<:Integer} <: PrimitiveBlock{N}
    Measure(n::Int; rng=Random.GLOBAL_RNG, operator=ComputationalBasis(), locs=1:n, collapseto=nothing, remove=false)

Measure operator.
"""
mutable struct Measure{N,K,OT,RNG,IT<:Integer} <: PrimitiveBlock{N}
    rng::RNG
    operator::OT
    locations::Union{NTuple{K,Int},AllLocs}
    collapseto::Union{BitStr64{N},Nothing}
    remove::Bool
    results::Vector{IT}
    function Measure{N,K,OT,RNG,IT}(rng::RNG, operator, locations, collapseto, remove) where {RNG, N, K, OT, IT}
        locations isa AllLocs || @assert_locs_safe N locations
        if collapseto !== nothing && remove == true
            error("invalid keyword combination, expect collapseto or remove, got (collapseto=$collapseto, remove=true)")
        end
        new{N,K,OT,RNG,IT}(rng, operator, locations, collapseto, remove)
    end
end

"""
    chmeasureoperator(m::Measure, op::AbstractBlock)

change the measuring `operator`. It will also discard existing measuring results.
"""
function chmeasureoperator(m::Measure{N}, op::AbstractBlock) where {N}
    Measure(
        N;
        rng = m.rng,
        operator = op,
        locs = m.locations,
        collapseto = m.collapseto,
        remove = m.remove,
    )
end

function Base.:(==)(m1::Measure, m2::Measure)
    res = m1.rng == m2.rng &&
          m1.operator == m2.operator &&
          m1.locations == m2.locations && m1.collapseto == m2.collapseto && m1.remove == m2.remove
    res = res && isdefined(m1, :results) == isdefined(m2, :results)
    res && (!isdefined(m1, :results) || m1.results == m2.results)
end

@interface nqubits_measured(::Measure{N,K}) where {N,K} = K

"""
    Measure(n::Int; rng=Random.GLOBAL_RNG, operator=ComputationalBasis(), locs=AllLocs(), collapseto=nothing, remove=false)

Create a `Measure` block with number of qubits `n`.

# Example

You can create a `Measure` block on given basis (default is the computational basis).

```jldoctest; setup=:(using YaoBlocks)
julia> Measure(4)
Measure(4)
```

Or you could specify which qubits you are going to measure

```jldoctest; setup=:(using YaoBlocks)
julia> Measure(4; locs=1:3)
Measure(4;locs=(1, 2, 3))
```

by default this will collapse the current register to measure results.

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister, BitBasis, Random; Random.seed!(2))
julia> r = rand_state(3)
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 3/3

julia> state(r)
8×1 Array{Complex{Float64},2}:
    0.21633342515406265 - 0.21776267239802458im
   -0.17798384008375148 - 0.5040979387214165im
   -0.19761243345925425 + 0.16281482444784728im
   -0.25200691415025867 + 0.15153595884416518im
     0.3650977378140692 + 0.3419566592091794im
  -0.027207023333497483 - 0.3780181361735894im
 -0.0034728372576413743 + 0.1693915490059622im
   -0.19898587237095824 - 0.07607057769761456im

julia> r |> Measure(3)
Measure(3)

julia> state(r)
8×1 Array{Complex{Float64},2}:
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7298587746534583 + 0.6835979586433478im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
```

But you can also specify the target bit configuration you want to collapse to with keyword `collapseto`.

```jldoctest; setup=:(using YaoBlocks; using BitBasis)
julia> m = Measure(4; collapseto=bit"0101")
Measure(4;collapseto=0101 ₍₂₎)

julia> m.collapseto
0101 ₍₂₎
```
"""
function Measure(n::Int; rng::RNG=Random.GLOBAL_RNG, operator::OT=ComputationalBasis(), locs=AllLocs(), collapseto=nothing, remove=false, result_dtype=BitStr64{n}) where {OT, RNG}
    if locs isa AllLocs
        Measure{n, n, OT, RNG, result_dtype}(rng, operator, locs, collapseto, remove)
    else
        Measure{n, length(locs), OT, RNG, result_dtype}(rng, operator, tuple(locs...), collapseto, remove)
    end
end

Measure(
    ;
    rng = Random.GLOBAL_RNG,
    locs = AllLocs(),
    operator = ComputationalBasis(),
    collapseto = nothing,
    remove = false,
) where {K} = @λ(
    n -> Measure(
        n;
        rng = rng,
        locs = locs,
        operator = operator,
        collapseto = collapseto,
        remove = remove,
    )
)
mat(x::Measure) = error("use BlockMap to get its matrix.")

function apply!(r::AbstractRegister, m::Measure{N}) where {N}
    _check_size(r, m)
    if m.collapseto !== nothing
        m.results = measure_collapseto!(m.rng, m.operator, r, m.locations; config = m.collapseto)
    elseif m.remove
        m.results = measure_remove!(m.rng, m.operator, r, m.locations)
    else
        m.results = measure!(m.rng, m.operator, r, m.locations)
    end
    return m
end


text:
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using CMSSPs

inhor_file = ARGS[1]
outhor_file = ARGS[2]
params_fn = ARGS[3]
trials = 5

rng = MersenneTwister(5542)

modal_horizon_policy = load_modal_horizon_policy_localapprox(inhor_file, outhor_file)

# Create parameters and then CMSSP
params = toy2d_parse_params(params_fn)
cmssp = create_toy2d_cmssp(params)

# Set up the modal mdp
# Need to do it for just one mode
modal_mdp = Toy2DModalMDPType(1, params, cmssp.control_actions, 0.8, params.horizon_limit)

for i = 1:trials

    tot_cost = 0.0

    target_state = sample_toy2d(rng)
    curr_state = Toy2DContState(0.0,0.0)

    @show target_state
    inf_hor_rch = 0

    while true
        (a,inf_hor_rch) = get_best_intramodal_action_infhor(modal_horizon_policy, inf_hor_rch, curr_state, target_state)
        (new_state, reward) = generate_sr(modal_mdp, curr_state, a.action, rng)
        curr_state = new_state
        tot_cost += -1.0*reward
        if POMDPs.isterminal(modal_mdp,get_relative_state(modal_mdp, curr_state,target_state))
            println("Success!")
            break
        end
    end
end


text:
module Cinta
export LineRetrieverFromString, LineRetrieverFromFile, Parameters, get_next_line, is_header, is_detail, get_output

abstract type LineRetriever end

struct Parameters
		base::Float32
		variation::Float32
end

function sorted_by_schedule(m::Array{Any})
	schedule = Dict([
			(1, 1), (7, 2),
			(2, 3), (8, 4),
			(3, 5), (9, 6),
			(4, 7), (10, 8),
			(5, 9), (11, 10),
			(6, 11), (12, 12),
			])
	return sort!(m, by=r->schedule[r[1]])
end

struct LineRetrieverFromString <: LineRetriever
	buffer::IOBuffer
	function LineRetrieverFromString(received_text::String)
			new(IOBuffer(received_text))
	end
end

struct LineRetrieverFromFile <: LineRetriever
	buffer::IOStream
	function LineRetrieverFromFile(file_name::String)
			new(open(file_name))
	end
end

get_next_line(t::LineRetrieverFromString) = eof(t.buffer) ? nothing : readline(t.buffer)
get_next_line(t::LineRetrieverFromFile) = eof(t.buffer) ? nothing : readline(t.buffer)

function is_header(line::String)
		header_re = r"cinta ruta [0-9]+ [0-9]+$"
		occursin(header_re, line)
end

function is_detail(line::String)
		detail_re = r"cinta [0-9]+$"
		occursin(detail_re, line)
end

function get_output(p::Parameters, t::LineRetriever)
		step = 0
		ruta = 0
		cab = r"cinta ruta [0-9]+ [0-9]+$"
		det = r"cinta [0-9]+$"
		ref = 3
		button = 9
		up = Char(0xFE3D)
		down = Char(0xFE3E)
		eq = Char(0x1d117)
		dont_touch = Char(0xa555)
		output = []

		line = get_next_line(t)
		while !isnothing(line)
				if is_header(line)
						step = 1
						m = match(cab, line).match
						ruta, vel_inc = parse(Int, split(m)[3]), parse(Int, split(m)[4])
						vel = Int(trunc(vel_inc / 100))
						inc = vel_inc % 100
						new_vel = p.base + (ref - inc) * p.variation
						button = abs(new_vel - 9) < abs(new_vel - 12) ? 9 : 12
						closest = abs(new_vel - vel) < abs(new_vel - button) ? vel : button
						append!(output, [[ruta, inc, string(43 - 4step) * ':', vel, '→', vel == closest ? dont_touch : button, new_vel == closest ? eq : new_vel > closest ? up : down, new_vel]])
						step += 1
				elseif is_detail(line)
						m = match(det, line).match
						vel_inc = parse(Int, split(m)[2])
						vel = Int(trunc(vel_inc / 100))
						inc = vel_inc % 100
						new_vel = p.base + (ref - inc) * p.variation
						button = abs(new_vel - 9) < abs(new_vel - 12) ? 9 : 12
						closest = abs(new_vel - vel) < abs(new_vel - button) ? vel : button
						append!(output, [[ruta, inc, string(43 - 4step) * ':', vel, '→', vel == closest ? dont_touch : button, new_vel == closest ? eq : new_vel > closest ? up : down, new_vel]])
						step += 1
				end
				line = get_next_line(t)
		end
		output = sorted_by_schedule(output)
		str_out = ""
		for r in output
				str_out *= join(r, '\t') * '\n'
		end
		return str_out
end

end

max_stars_repo_name: AndresParraSilva/julia-start

max_issues_repo_name: AndresParraSilva/julia-start

max_forks_repo_name: AndresParraSilva/julia-start


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@doc raw"""
    EulerAcousticsCouplingCallback

!!! warning "Experimental code"
    This callback is experimental and may change in any future release.

A callback that couples the acoustic perturbation equations and compressible Euler equations. Must
be used in conjunction with [`SemidiscretizationEulerAcoustics`](@ref).
This callback manages the flow solver - which is always one time step ahead of the
acoustics solver - and calculates the acoustic source term after each time step. The linearized
Lamb vector is used as the source term, i.e.
```math
\mathbf{s} = -(\mathbf{\omega'} \times \bar{\mathbf{v}}
  + \bar{\mathbf{\omega}} \times \mathbf{v'}),
```
where ``\mathbf{v}`` denotes the velocity, ``\mathbf{\omega}`` denotes the vorticity, the bar
``\bar{(\cdot)}`` indicates time-averaged quantities (see [`AveragingCallback`](@ref)) and prime
``(\cdot)'`` denotes perturbed quantities defined by ``\phi' = \phi - \bar{\phi}``. Note that
the perturbed quantities here are based entirely on the pure flow solution and should not be
confused with the state variables of the acoustic perturbation equations.

In addition, this callback manages the time step size for both solvers
and initializes the mean values of the acoustic perturbation equations using results obtained with
the [`AveragingCallback`](@ref).

- Michael Schlottke-Lakemper (2017)
  A direct-hybrid method for aeroacoustic analysis
  [DOI: 10.18154/RWTH-2017-04082](https://doi.org/10.18154/RWTH-2017-04082)
"""
mutable struct EulerAcousticsCouplingCallback{RealT<:Real, MeanValues, IntegratorEuler}
  stepsize_callback_acoustics::StepsizeCallback{RealT}
  stepsize_callback_euler::StepsizeCallback{RealT}
  mean_values::MeanValues
  integrator_euler::IntegratorEuler
end

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EulerAcousticsCouplingCallback})
  @nospecialize cb # reduce precompilation time
  euler_acoustics_coupling = cb.affect!

  print(io, "EulerAcousticsCouplingCallback(")
  print(io,       euler_acoustics_coupling.stepsize_callback_acoustics)
  print(io, ", ", euler_acoustics_coupling.stepsize_callback_euler, ")")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:EulerAcousticsCouplingCallback})
  @nospecialize cb # reduce precompilation time
  euler_acoustics_coupling = cb.affect!

  summary_header(io, "EulerAcousticsCouplingCallback")
  summary_line(io, "acoustics StepsizeCallback", euler_acoustics_coupling.stepsize_callback_acoustics)
  summary_line(io, "Euler StepsizeCallback", euler_acoustics_coupling.stepsize_callback_euler)
  summary_footer(io)
end

"""
    EulerAcousticsCouplingCallback(ode_euler,
                                   averaging_callback::DiscreteCallback{<:Any, <:AveragingCallback},
                                   alg, cfl_acoustics::Real, cfl_euler::Real; kwargs...)

!!! warning "Experimental code"
    This callback is experimental and may change in any future release.

Creates an [`EulerAcousticsCouplingCallback`](@ref) based on the pure flow `ODEProblem` given by
`ode_euler`. Creates an integrator using the time integration method `alg` and the keyword arguments
to solve `ode_euler` (consult the [OrdinaryDiffEq documentation](https://diffeq.sciml.ai/stable/)
for further information).
Manages the step size for both solvers by using the minimum of the maximum step size obtained with
CFL numbers `cfl_acoustics` for the acoustics solver and `cfl_euler` for and flow solver,
respectively.
The mean values for the acoustic perturbation equations are read from `averaging_callback`
(see [`AveragingCallback`](@ref)).
"""
function EulerAcousticsCouplingCallback(ode_euler,
                                        averaging_callback::DiscreteCallback{<:Any, <:AveragingCallback},
                                        alg, cfl_acoustics::Real, cfl_euler::Real; kwargs...)
  @unpack mean_values = averaging_callback.affect!

  return EulerAcousticsCouplingCallback(ode_euler, mean_values, alg, cfl_acoustics, cfl_euler;
                                        kwargs...)
end

"""
    EulerAcousticsCouplingCallback(ode_euler, averaging_file::AbstractString, alg,
                                   cfl_acoustics::Real, cfl_euler::Real; kwargs...)

!!! warning "Experimental code"
    This callback is experimental and may change in any future release.

Creates an [`EulerAcousticsCouplingCallback`](@ref) based on the pure flow `ODEProblem` given by
`ode_euler`. Creates an integrator using the time integration method `alg` and the keyword arguments
to solve `ode_euler` (consult the [OrdinaryDiffEq documentation](https://diffeq.sciml.ai/stable/)
for further information).
Manages the step size for both solvers by using the minimum of the maximum step size obtained with
CFL numbers `cfl_acoustics` for the acoustics solver and `cfl_euler` for and flow solver,
respectively.
The mean values for the acoustic perturbation equations are read from `averaging_file`
(see [`AveragingCallback`](@ref)).
"""
function EulerAcousticsCouplingCallback(ode_euler, averaging_file::AbstractString, alg,
                                        cfl_acoustics::Real, cfl_euler::Real; kwargs...)
  semi_euler = ode_euler.p
  mean_values = load_averaging_file(averaging_file, semi_euler)

  return EulerAcousticsCouplingCallback(ode_euler, mean_values, alg, cfl_acoustics, cfl_euler;
                                        kwargs...)
end

function EulerAcousticsCouplingCallback(ode_euler, mean_values, alg, cfl_acoustics, cfl_euler;
                                        kwargs...)
  # Set up ODE Integrator for Euler equations
  integrator_euler = init(ode_euler, alg, save_everystep=false, dt=1.0; kwargs...) # dt will be overwritten

  euler_acoustics_coupling = EulerAcousticsCouplingCallback{typeof(cfl_acoustics),
                                                            typeof(mean_values),
                                                            typeof(integrator_euler)}(
    StepsizeCallback(cfl_acoustics), StepsizeCallback(cfl_euler), mean_values, integrator_euler)
  condition = (u, t, integrator) -> true

  return DiscreteCallback(condition, euler_acoustics_coupling, save_positions=(false, false),
                          initialize=initialize!)
end

# This is called before the main loop and initializes the mean values in u_ode
function initialize!(cb::DiscreteCallback{Condition,Affect!}, u_ode, t, integrator_acoustics) where {Condition, Affect!<:EulerAcousticsCouplingCallback}
  euler_acoustics_coupling = cb.affect!
  semi = integrator_acoustics.p
  @unpack semi_acoustics = semi

  # Initialize mean values in u_ode
  u_acoustics = wrap_array(u_ode, semi_acoustics)
  @unpack mean_values = euler_acoustics_coupling
  @views @. u_acoustics[4:5, .., :] = mean_values.v_mean
  @views @. u_acoustics[6, .., :] = mean_values.c_mean
  @views @. u_acoustics[7, .., :] = mean_values.rho_mean

  # Adjust stepsize, advance the flow solver by one time step
  cb.affect!(integrator_acoustics)

  return nothing
end

# This function is called at the end of every time step and advances the Euler solution by one
# time step, manages the time stepsize for both the acoustics and Euler solvers and calculates the
# acoustic sources for the next acoustics time step
function (euler_acoustics_coupling::EulerAcousticsCouplingCallback)(integrator_acoustics)
  @unpack stepsize_callback_acoustics, stepsize_callback_euler, integrator_euler = euler_acoustics_coupling

  @assert integrator_acoustics.t == integrator_euler.t

  # Use the minimum of the acoustics and Euler stepsizes for both solvers
  stepsize_callback_acoustics(integrator_acoustics)
  stepsize_callback_euler(integrator_euler)
  dt = min(get_proposed_dt(integrator_acoustics), get_proposed_dt(integrator_euler))

  set_proposed_dt!(integrator_acoustics, dt)
  integrator_acoustics.opts.dtmax = dt
  integrator_acoustics.dtcache = dt

  set_proposed_dt!(integrator_euler, dt)
  integrator_euler.opts.dtmax = dt
  integrator_euler.dtcache = dt

  # Advance the Euler solution by one step and check for errors
  if !isfinished(integrator_euler)
    @trixi_timeit timer() "Euler solver" step!(integrator_euler)
    return_code = check_error(integrator_euler)
    if return_code !== :Success && return_code !== :Default
      error("Error during compressible Euler time integration. Received return code $(return_code)")
    end
  end

  # Calculate acoustic sources based on linearized lamb vector
  semi = integrator_acoustics.p
  semi_euler = integrator_euler.p
  u_acoustics = wrap_array(integrator_acoustics.u, semi)
  u_euler = wrap_array(integrator_euler.u, semi_euler)
  @unpack acoustic_source_terms, coupled_element_ids = semi.cache
  @unpack vorticity_mean = euler_acoustics_coupling.mean_values

  @trixi_timeit timer() "calc acoustic source terms" calc_acoustic_sources!(
    acoustic_source_terms, u_euler, u_acoustics, vorticity_mean, coupled_element_ids,
    mesh_equations_solver_cache(semi_euler)...)

  # avoid re-evaluation possible FSAL stages
  u_modified!(integrator_acoustics, false)
  u_modified!(integrator_euler, false)

  return nothing
end

include("euler_acoustics_coupling_dg2d.jl")

end # @muladd


text:
using AbstractDAQs
using Documenter

DocMeta.setdocmeta!(AbstractDAQs, :DocTestSetup, :(using AbstractDAQs); recursive=true)

makedocs(;
    modules=[AbstractDAQ],
    authors="Paulo Jabardo <pjabardo@ipt.br>",
    repo="https://github.com/pjsjipt/AbstractDAQs.jl/blob/{commit}{path}#{line}",
    sitename="AbstractDAQs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
# ## Non-zero inflated, continuous. e.g. LogNormal
# function expert_tn_bar_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:NonZIContinuousExpert}

#     expert_tn_bar = (tl == tu) ? 0.0 : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu)))

#     return expert_tn_bar
# end

# ## Zero inflated, continuous. e.g. ZILogNormal
# function expert_tn_bar_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:ZIContinuousExpert}

#     expert_tn_bar = (tl == tu) ? 0.0 : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu)))

#     return expert_tn_bar
# end

# ## Non-zero inflated, discrete. e.g. Poisson
# function expert_tn_bar_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:NonZIDiscreteExpert}

#     expert_tn_bar = (tl == tu) ? log1mexp.(logpdf.(d, tl)) : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu)))

#     return expert_tn_bar
# end

# ## Zero inflated, discrete. e.g. ZIPoisson
# function expert_tn_bar_pos(d::e, tl::Real, yl::Real, yu::Real, tu::Real) where {e<:ZIDiscreteExpert}

#     expert_tn_bar = (tl == tu) ? log1mexp.(logpdf.(d, tl)) : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu)))

#     return expert_tn_bar
# end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
using Base.Meta, Pkg
if Pkg.installed()["JuMP"] >= v"0.18.4"
    import JuMP.undef
end

"""
    @randomvariable(m, x, mean)
Defines a random variable using [`RandomVariable(m::Model, mean::Number, name::AbstractString)`](@ref)
and requires that a mean for the variable be provided. This can later be overwritten with [`setmean`](@ref).

**Arguments**
- `m::Model` The flexibility model.
- `x::Symbol` The variable name.
- `mean::Number` The variable mean

```julia
julia> @randomvariable(m2, w, mean = 42)

julia> @randomvariable(m2, ws[i = 1:4], mean = 42)
4-element Array{FlexJuMP.RandomVariable,1}:
 FlexJuMP.RandomVariable(Feasibility problem with:
 * 9 linear constraints
 * 89 variables

 julia> @randomvariable(m2, ws[i = 1:4], mean = [1; 2; 3; 4][i])
 4-element Array{FlexJuMP.RandomVariable,1}:
  FlexJuMP.RandomVariable(Feasibility problem with:
  * 9 linear constraints
  * 96 variables
```
"""
macro randomvariable(m, x, mean)
    m = esc(m)
    kwsymbol = VERSION < v"0.6.0-dev.1934" ? :kw : :(=) # changed by julia PR #19868
    if isexpr(mean,kwsymbol)
        @assert isexpr(mean,kwsymbol) && mean.args[1] == :mean
        mean = esc(mean.args[2])
    end

    if isa(x,Symbol)
        # easy case
        return quote
            $(esc(x)) = RandomVariable($m,$mean,$(string(x)))
            nothing
        end
    else
        if !isexpr(x,:ref)
            error("Syntax error: Expected $var to be of form var[...]")
        end

        variable = gensym()
        refcall, idxvars, idxsets, idxpairs, condition = JuMP.buildrefsets(x, variable)
        varname = JuMP.getname(x)
        escvarname = esc(varname)

        varstr = :(string($(string(varname)),"["))
        for idxvar in idxvars
            push!(varstr.args,:(string($(esc(idxvar)))))
            push!(varstr.args,",")
        end
        deleteat!(varstr.args,length(varstr.args))
        push!(varstr.args,"]")

        code = :( $(refcall) = RandomVariable($m, $mean, $varstr ) )
        looped = JuMP.getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :RandomVariable)
        return quote
            $looped
            $escvarname = $variable
        end
    end
end

"""
    @recoursevariable(m, x)
Defines a recourse variable using [`RecourseVariable(m::Model, name::AbstractString)`](@ref).

**Arguments**
- `m::Model` The flexibility model.
- `x::Symbol` The variable name.

```julia
julia> @recoursevariable(m2, d)

julia> @recoursevariable(m2, ds[1:4])
4-element Array{FlexJuMP.RecourseVariable,1}:
 FlexJuMP.RecourseVariable(Feasibility problem with:
 * 9 linear constraints
 * 101 variables
```
"""
macro recoursevariable(m, x)
    m = esc(m)
    kwsymbol = VERSION < v"0.6.0-dev.1934" ? :kw : :(=) # changed by julia PR #19868

    if isa(x,Symbol)
        # easy case
        return quote
            $(esc(x)) = RecourseVariable($m,$(string(x)))
            nothing
        end
    else
        if !isexpr(x,:ref)
            error("Syntax error: Expected $var to be of form var[...]")
        end

        variable = gensym()
        refcall, idxvars, idxsets, idxpairs, condition = JuMP.buildrefsets(x, variable)
        varname = JuMP.getname(x)
        escvarname = esc(varname)

        varstr = :(string($(string(varname)),"["))
        for idxvar in idxvars
            push!(varstr.args,:(string($(esc(idxvar)))))
            push!(varstr.args,",")
        end
        deleteat!(varstr.args,length(varstr.args))
        push!(varstr.args,"]")

        code = :( $(refcall) = RecourseVariable($m, $varstr ) )
        looped = JuMP.getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :RecourseVariable)
        return quote
            $looped
            $escvarname = $variable
        end
    end
end

"""
    JuMP.constructconstraint!(flex_aff::FlexibilityExpr, sense::Symbol)
Extends `JuMP.constructconstraint!` for `FlexibilityExpr` types.
"""
function JuMP.constructconstraint!(flex_aff::FlexibilityExpr, sense::Symbol)
    if sense == :(<=) || sense == :≤
        return FlexibilityConstraint(flex_aff, :(<=))
    elseif sense == :(>=) || sense == :≥
        return FlexibilityConstraint(flex_aff, :(>=))
    elseif sense == :(==)
        return FlexibilityConstraint(flex_aff, :(==))
    end
    error("Unrecognized constraint type $sense")
end