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

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

Runs a batch SQL statement over an array of data. You can run bulk update and insert operations for multiple records using a DML statement with different parameter sets. Bulk operations can provide a significant performance improvement over individual insert and update operations.  If a call isn't part of a transaction because it doesn't include the transactionID parameter, changes that result from the call are committed automatically.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secretArn`: The name or ARN of the secret that enables access to the DB cluster.
- `sql`: The SQL statement to run.

# Optional Parameters
- `database`: The name of the database.
- `parameterSets`: The parameter set for the batch operation. The SQL statement is executed as many times as the number of parameter sets provided. To execute a SQL statement with no parameters, use one of the following options:   Specify one or more empty parameter sets.   Use the ExecuteStatement operation instead of the BatchExecuteStatement operation.    Array parameters are not supported.
- `schema`: The name of the database schema.
- `transactionId`: The identifier of a transaction that was started by using the BeginTransaction operation. Specify the transaction ID of the transaction that you want to include the SQL statement in. If the SQL statement is not part of a transaction, don't set this parameter.
"""

batch_execute_statement(resourceArn, secretArn, sql; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/BatchExecute", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql); aws_config=aws_config)
batch_execute_statement(resourceArn, secretArn, sql, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/BatchExecute", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql), args)); aws_config=aws_config)

"""
    BeginTransaction()

Starts a SQL transaction.  <important>
A transaction can run for a maximum of 24 hours. A transaction is terminated and rolled back automatically after 24 hours.

A transaction times out if no calls use its transaction ID in three minutes. If a transaction times out before it's committed, it's rolled back automatically.

DDL statements inside a transaction cause an implicit commit. We recommend that you run each DDL statement in a separate ExecuteStatement call with continueAfterTimeout enabled.
 </important>

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secretArn`: The name or ARN of the secret that enables access to the DB cluster.

# Optional Parameters
- `database`: The name of the database.
- `schema`: The name of the database schema.
"""

begin_transaction(resourceArn, secretArn; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/BeginTransaction", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn); aws_config=aws_config)
begin_transaction(resourceArn, secretArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/BeginTransaction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn), args)); aws_config=aws_config)

"""
    CommitTransaction()

Ends a SQL transaction started with the BeginTransaction operation and commits the changes.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secretArn`: The name or ARN of the secret that enables access to the DB cluster.
- `transactionId`: The identifier of the transaction to end and commit.

"""

commit_transaction(resourceArn, secretArn, transactionId; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/CommitTransaction", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId); aws_config=aws_config)
commit_transaction(resourceArn, secretArn, transactionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/CommitTransaction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId), args)); aws_config=aws_config)

"""
    ExecuteSql()

Runs one or more SQL statements.  This operation is deprecated. Use the BatchExecuteStatement or ExecuteStatement operation.

# Required Parameters
- `awsSecretStoreArn`: The Amazon Resource Name (ARN) of the secret that enables access to the DB cluster.
- `dbClusterOrInstanceArn`: The ARN of the Aurora Serverless DB cluster.
- `sqlStatements`: One or more SQL statements to run on the DB cluster. You can separate SQL statements from each other with a semicolon (;). Any valid SQL statement is permitted, including data definition, data manipulation, and commit statements.

# Optional Parameters
- `database`: The name of the database.
- `schema`: The name of the database schema.
"""

execute_sql(awsSecretStoreArn, dbClusterOrInstanceArn, sqlStatements; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/ExecuteSql", Dict{String, Any}("awsSecretStoreArn"=>awsSecretStoreArn, "dbClusterOrInstanceArn"=>dbClusterOrInstanceArn, "sqlStatements"=>sqlStatements); aws_config=aws_config)
execute_sql(awsSecretStoreArn, dbClusterOrInstanceArn, sqlStatements, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/ExecuteSql", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("awsSecretStoreArn"=>awsSecretStoreArn, "dbClusterOrInstanceArn"=>dbClusterOrInstanceArn, "sqlStatements"=>sqlStatements), args)); aws_config=aws_config)

"""
    ExecuteStatement()

Runs a SQL statement against a database.  If a call isn't part of a transaction because it doesn't include the transactionID parameter, changes that result from the call are committed automatically.  The response size limit is 1 MB. If the call returns more than 1 MB of response data, the call is terminated.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secretArn`: The name or ARN of the secret that enables access to the DB cluster.
- `sql`: The SQL statement to run.

# Optional Parameters
- `continueAfterTimeout`: A value that indicates whether to continue running the statement after the call times out. By default, the statement stops running when the call times out.  For DDL statements, we recommend continuing to run the statement after the call times out. When a DDL statement terminates before it is finished running, it can result in errors and possibly corrupted data structures.
- `database`: The name of the database.
- `includeResultMetadata`: A value that indicates whether to include metadata in the results.
- `parameters`: The parameters for the SQL statement.  Array parameters are not supported.
- `resultSetOptions`: Options that control how the result set is returned.
- `schema`: The name of the database schema.
- `transactionId`: The identifier of a transaction that was started by using the BeginTransaction operation. Specify the transaction ID of the transaction that you want to include the SQL statement in. If the SQL statement is not part of a transaction, don't set this parameter.
"""

execute_statement(resourceArn, secretArn, sql; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/Execute", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql); aws_config=aws_config)
execute_statement(resourceArn, secretArn, sql, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/Execute", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql), args)); aws_config=aws_config)

"""
    RollbackTransaction()

Performs a rollback of a transaction. Rolling back a transaction cancels its changes.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secretArn`: The name or ARN of the secret that enables access to the DB cluster.
- `transactionId`: The identifier of the transaction to roll back.

"""

rollback_transaction(resourceArn, secretArn, transactionId; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/RollbackTransaction", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId); aws_config=aws_config)
rollback_transaction(resourceArn, secretArn, transactionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = rds_data("POST", "/RollbackTransaction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId), args)); aws_config=aws_config)


text:
function codegen end

include("utils.jl")
include("emulation.jl")
include("qasm.jl")
include("qobj.jl")

ext: jl


text:
@doc_str """
    push_relabel(residual_graph, source, target, capacity_matrix)

Return the maximum flow of `residual_graph` from `source` to `target` using the
FIFO push relabel algorithm with gap heuristic.

### Performance
Takes approximately ``\\mathcal{O}(|V|^{3})`` time.
"""
function push_relabel end
@traitfn function push_relabel(
        residual_graph::::lg.IsDirected,       # the input graph
        source::Integer,                       # the source vertex
        target::Integer,                       # the target vertex
        capacity_matrix::AbstractMatrix{T}     # edge flow capacities
    ) where {T}

    n = lg.nv(residual_graph)
    flow_matrix = zeros(T, n, n)

    height = zeros(Int, n)
    height[source] = n

    count = zeros(Int, 2 * n + 1)
    count[0 + 1] = n - 1
    count[n + 1] = 1

    excess = zeros(T, n)
    excess[source] = typemax(T)

    active = falses(n)
    active[source] = true
    active[target] = true

    Q = Array{Int,1}()
    sizehint!(Q, n)

    for v in lg.outneighbors(residual_graph, source)
        push_flow!(residual_graph, source, v, capacity_matrix, flow_matrix, excess, height, active, Q)
    end

    while length(Q) > 0
        v = pop!(Q)
        active[v] = false
        discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
    end

    return sum([flow_matrix[v, target] for v in lg.inneighbors(residual_graph, target)]), flow_matrix
end

"""
    enqueue_vertex!(Q, v, active, excess)

Push inactive node `v` into queue `Q` and activates it. Requires preallocated
`active` and `excess` vectors.
"""
function enqueue_vertex!(
        Q::AbstractVector,
        v::Integer,                   # input vertex
        active::AbstractVector{Bool},
        excess::AbstractVector
    )
    if !active[v] && excess[v] > 0
        active[v] = true
        pushfirst!(Q, v)
    end
    return nothing
end

"""
    push_flow!(residual_graph, u, v, capacity_matrix, flow_matrix, excess, height, active, Q)

Using `residual_graph` with capacities in `capacity_matrix`, push as much flow
as possible through the given edge(`u`, `v`). Requires preallocated `flow_matrix`
matrix, and `excess`, `height, `active`, and `Q` vectors.
"""
function push_flow! end
@traitfn function push_flow!(
        residual_graph::::lg.IsDirected,   # the input graph
        u::Integer,                        # input from-vertex
        v::Integer,                        # input to-vetex
        capacity_matrix::AbstractMatrix,
        flow_matrix::AbstractMatrix,
        excess::AbstractVector,
        height::AbstractVector{Int},
        active::AbstractVector{Bool},
        Q::AbstractVector
    )
    flow = min(excess[u], capacity_matrix[u, v] - flow_matrix[u, v])

    flow == 0 && return nothing
    height[u] <= height[v] && return nothing

    flow_matrix[u, v] += flow
    flow_matrix[v, u] -= flow

    excess[u] -= flow
    excess[v] += flow

    enqueue_vertex!(Q, v, active, excess)
    nothing
end

"""
    gap!(residual_graph, h, excess, height, active, count, Q)

Implement the push-relabel gap heuristic. Relabel all vertices above a cutoff height.
Reduce the number of relabels required.

Requires arguments:

- residual_graph::DiGraph                # the input graph
- h::Int                                 # cutoff height
- excess::AbstractVector
- height::AbstractVector{Int}
- active::AbstractVector{Bool}
- count::AbstractVector{Int}
- Q::AbstractVector
"""
function gap! end
@traitfn function gap!(
        residual_graph::::lg.IsDirected,       # the input graph
        h::Int,                                # cutoff height
        excess::AbstractVector,
        height::AbstractVector{Int},
        active::AbstractVector{Bool},
        count::AbstractVector{Int},
        Q::AbstractVector                      # FIFO queue
    )
    n = lg.nv(residual_graph)
    for v in lg.vertices(residual_graph)
        height[v] < h && continue
        count[height[v] + 1] -= 1
        height[v] = max(height[v], n + 1)
        count[height[v] + 1] += 1
        enqueue_vertex!(Q, v, active, excess)
    end
    nothing
end

"""
    relabel!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)

Relabel a node `v` with respect to its neighbors to produce an admissable edge.
"""
function relabel! end
@traitfn function relabel!(
        residual_graph::::lg.IsDirected,   # the input graph
        v::Integer,                        # input vertex to be relabeled
        capacity_matrix::AbstractMatrix,
        flow_matrix::AbstractMatrix,
        excess::AbstractVector,
        height::AbstractVector{Int},
        active::AbstractVector{Bool},
        count::AbstractVector{Int},
        Q::AbstractVector
    )
    n = lg.nv(residual_graph)
    count[height[v] + 1] -= 1
    height[v] = 2 * n
    for to in lg.outneighbors(residual_graph, v)
        if capacity_matrix[v, to] > flow_matrix[v, to]
            height[v] = min(height[v], height[to] + 1)
        end
    end
    count[height[v] + 1] += 1
    enqueue_vertex!(Q, v, active, excess)
    nothing
end

"""
    discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)

Drain the excess flow out of node `v`. Run the gap heuristic or relabel the
vertex if the excess remains non-zero.
"""
function discharge! end
@traitfn function discharge!(
        residual_graph::::lg.IsDirected,    # the input graph
        v::Integer,                         # vertex to be discharged
        capacity_matrix::AbstractMatrix,
        flow_matrix::AbstractMatrix,
        excess::AbstractVector,
        height::AbstractVector{Int},
        active::AbstractVector{Bool},
        count::AbstractVector{Int},
        Q::AbstractVector                   # FIFO queue
    )
    for to in lg.outneighbors(residual_graph, v)
        excess[v] == 0 && break
        push_flow!(residual_graph, v, to, capacity_matrix, flow_matrix, excess, height, active, Q)
    end

    if excess[v] > 0
        if count[height[v] + 1] == 1
            gap!(residual_graph, height[v], excess, height, active, count, Q)
        else
            relabel!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
        end
    end
    nothing
end


text:
macro inc(a)
    :(@testset $a begin let n, p, i; include($a) end end)
end

@testset "MatrixDepot generator tests" begin

    @inc("test_magic.jl")
    @inc("test_cauchy.jl")
    @inc("test_circul.jl")
    @inc("test_hadamard.jl")
    @inc("test_hilb.jl")
    @inc("test_dingdong.jl")
    @inc("test_frank.jl")
    @inc("test_invhilb.jl")
    @inc("test_forsythe.jl")
    @inc("test_grcar.jl")
    @inc("test_triw.jl")
    @inc("test_moler.jl")
    @inc("test_pascal.jl")
    @inc("test_kahan.jl")
    @inc("test_pei.jl")
    @inc("test_vand.jl")
    @inc("test_invol.jl")
    @inc("test_chebspec.jl")
    @inc("test_lotkin.jl")
    @inc("test_clement.jl")
    @inc("test_fiedler.jl")
    @inc("test_minij.jl")
    @inc("test_binomial.jl")
    @inc("test_tridiag.jl")
    @inc("test_lehmer.jl")
    @inc("test_parter.jl")
    @inc("test_chow.jl")
    @inc("test_randcorr.jl")
    @inc("test_poisson.jl")
    @inc("test_neumann.jl")
    @inc("test_rosser.jl")
    @inc("test_sampling.jl")
    @inc("test_wilkinson.jl")
    @inc("test_rando.jl")
    @inc("test_randsvd.jl")
    @inc("test_rohess.jl")
    @inc("test_kms.jl")
    @inc("test_wathen.jl")
    @inc("test_oscillate.jl")
    @inc("test_toeplitz.jl")
    @inc("test_hankel.jl")
    @inc("test_golub.jl")
    @inc("test_companion.jl")
    @inc("test_erdrey.jl")
    @inc("test_gilbert.jl")
    @inc("test_smallworld.jl")
end

@testset "regu and include tests" begin
    include("regu.jl")
    include("include_generator.jl")
end


text:
macro define_traitfn(trait,traitfndef)
    !(@capture(trait,Trait_Symbol)) && error("Please provide a trait to dispatch on")
    (fndef, impls) = @match traitfndef begin
        fn_(args__) => (:($traitfndef = begin end),Expr(:block))
        fn_(args__) where wargs__ => (:($traitfndef = begin end),Expr(:block))
        (fndef_ = body_) => (:($fndef = begin end),block(prettify(body)))
    end
    # Match function definition
    traitfn_def_split = try
        splitdef(fndef)
    catch er
        error("Invalid trait function syntax. Parser says: ", er.msg)
    end
    isempty(traitfn_def_split[:args]) && error("No arguments to (trait)dispatch on")

    # Almost the same definition will be used to implement the NullTrait- and NotImplemented behaviour
    nulltrait_split = deepcopy(traitfn_def_split)
    noimplement_split = deepcopy(traitfn_def_split)

    impltraits = Symbol[]
    if !isempty(impls.args)
        for (i,fnimpl) in enumerate(impls.args)
            (fn,args,wparams,impltrait) = @match fnimpl begin
                (fn_(args__,trait_) = body_) => (fn,args,[],trait)
                (fn_(args__,trait_) where wparams__ = body_) => (fn,args,wparams,trait)
                (function fn_(args__,trait_)
                    body_
                 end) => (fn,args,[],trait)
                (function fn_(args__,trait_) where wparams__
                    body_
                end) => (fn,args,wparams,trait)
                _ => error("Syntax error in implementation list. Provide function implementations with the trait as last argument.")
            end
            if isa(impltrait,Expr)
                if impltrait.head == :call && impltrait.args[1] == :! && impltrait.args[2] == trait
                    impltrait = :NullTrait
                    if !isempty(wparams)
                        impls.args[i].args[1].args[1].args[end] = :NullTrait
                    else
                        impls.args[i].args[1].args[end] = :NullTrait
                    end
                else
                    error("Syntax error in last argument of implementation list. Use !Trait to specify that the function is implemented for types that do not implement Trait.")
                end
            end
            isempty(args) && error("Syntax error in implementation list. Provide function implementations with the trait as last argument.")
            length(args) < 1 && error("Syntax error in implementation list. Provide function implementations with the trait as last argument.")
            !isa(impltrait,Symbol) && error("Syntax error in last argument of implementation list. Specify for which trait the function is being implemented as last argument.")
            fn != traitfn_def_split[:name] && error("Name of implemented function does not match defined function.")
            !(all(args .== traitfn_def_split[:args])) && error("Inconsistent arguments in definition and implementation, in function ", traitfn_def_split[:name])
            if !isempty(wparams)
                !(all(wparams .== traitfn_def_split[:whereparams])) && error("Inconsistent where parameters in definition and implementation, in function ",traitfn_def_split[:name])
            end
            push!(impltraits,impltrait)
        end
    end

    # Prepare the lhs of the traitfn def.
    parname = prepare_dispatchparameter!(traitfn_def_split)

    # Change rhs of traitfn definition to trait dispatch call
    rhs = Expr(:block)
    push!(rhs.args,:(Base.@_inline_meta()))
    call = Expr(:call)
    push!(call.args,traitfn_def_split[:name])
    append!(call.args,[isa(arg,Symbol) ? arg : arg.args[1] for arg in traitfn_def_split[:args]])
    traitpick = @q begin
        TraitDispatch.traitdispatch($trait{$parname})
    end
    push!(call.args,traitpick)
    push!(rhs.args,call)
    traitfn_def_split[:body] = rhs

    # Construct the traitfn definition
    traitfn_def = combinecall(traitfn_def_split)

    # Finish NullTrait implementation with NullTrait dispatch
    push!(nulltrait_split[:args],:NullTrait)
    nulltrait_split[:body] = :(error($(nulltrait_split[:args][1])," has no applicable traits"))
    nulltrait_impl = combinecall(nulltrait_split)

    # Finish NotImplemented implementation with dispatch on parent trait
    push!(noimplement_split[:args],:($trait))
    noimplement_split[:body] = :(error("Trait function ",$(traitfn_def_split[:name])," not implemented for trait ",Trait.name))
    noimplement_impl = combinecall(noimplement_split)

    # Construct the definition code, with sanity checks
    code = @q begin
        # Sanity checks
        $(esc(trait)) == AbstractTrait && error("Cannot define trait function for the AbstractTrait")
        $(esc(trait)) == NullTrait && error("Cannot define trait function for the NullTrait")
        !($(esc(trait)) <: AbstractTrait) && error($(esc(trait))," is not a trait")
        # Definition
        $(esc(traitfn_def))
    end
    # Provided implementations
    for traitfn_impl in impls.args
        push!(code.args,esc(:(@implement_traitfn $traitfn_impl)))
    end
    # Default implementations
    if !(:NullTrait in impltraits)
        push!(code.args,esc(:(@implement_traitfn $nulltrait_impl)))
    end
    if !(trait in impltraits)
        push!(code.args,esc(@q begin
            if !isempty(subtraits($trait))
                @implement_traitfn $noimplement_impl
            end
        end))
    end
    prettify(code)
end

macro implement_traitfn(fndef)
    # Match function definition
    traitfn_impl_split = try
        splitdef(fndef)
    catch er
        error("Invalid trait function syntax. Parser says: ", er.msg)
    end
    isempty(traitfn_impl_split[:args]) && error("No arguments to (trait)dispatch on")
    length(traitfn_impl_split[:args]) == 1 && error("Provide at least one argument to trait dispatch on, and one trait")
    trait = pop!(traitfn_impl_split[:args])

    # Prepare the lhs of the traitfn impl.
    parname = prepare_dispatchparameter!(traitfn_impl_split)

    # Finish the lhs of the traitfn implementation, by adding TraitDispatch
    push!(traitfn_impl_split[:args],:(::Type{Trait}))
    traitfn_impl_split[:whereparams] = (traitfn_impl_split[:whereparams]...,:(Trait <: $trait{$parname}))

    # Construct the traitfn definition
    traitfn_impl = combinecall(traitfn_impl_split)

    # Construct the definition code, with sanity checks
    code = @q begin
        # Sanity checks
        $(esc(trait)) == AbstractTrait && error("Cannot implement trait function for the AbstractTrait")
        $(esc(trait)) != NullTrait && !($(esc(trait)) <: AbstractTrait) && error($(esc(trait))," is not a trait. Provide a trait as the final argument")

        # Implementation
        $(esc(traitfn_impl))
    end
    prettify(code)
end

function prepare_dispatchparameter!(fn_split::Dict)
    # Prepare the argument x that will be subjected to trait dispatch
    first_arg = fn_split[:args][1]
    x,x_type = @capture(first_arg,x_::T_) ? (x,T) : (first_arg,:none)

    # Make sure parameter restriction matches type restriction on x
    parname = :T
    if x_type == :none
        fn_split[:whereparams] = tuple(fn_split[:whereparams]...,:T)
    else
        if isempty(fn_split[:whereparams])
            fn_split[:whereparams] = tuple(fn_split[:whereparams]...,:(T <: $x_type))
        else
            if !(x_type ∈ fn_split[:whereparams])
                if !any(map(wparam -> @capture(wparam,T_ <: types_) && x_type == T,fn_split[:whereparams]))
                    fn_split[:whereparams] = tuple(fn_split[:whereparams]...,:(VarType <: $x_type))
                    parname = :VarType
                end
            end
        end
    end
    fn_split[:args][1] = :($x::$parname)

    parname
end

function combinecall(dict::Dict)
    params = get(dict, :params, [])
    wparams = get(dict, :whereparams, [])
    name = dict[:name]
    name_param = isempty(params) ? name : :($name{$(params...)})
    if isempty(wparams)
        :(function $name_param($(dict[:args]...);$(dict[:kwargs]...))
          $(dict[:body])
          end)
    else
        :(function $name_param($(dict[:args]...);$(dict[:kwargs]...)) where {$(wparams...)}
          $(dict[:body])
          end)
    end
end


text:
function asymmetry(df, id, x, y)
    result = DataFrame(TrackID = Int[], len = Int[], asymmetry = Float64[])
    @inbounds for n = minimum(df[!, id]):maximum(df[!, id])
        m = extract(df, n, id, [x, y])
        T = Matrix{Float64}(undef, 2, 2)
        mean_x = @views mean(m[:, 1])
        mean_y = @views mean(m[:, 2])
        t1 = 0.0
        t2 = 0.0
        t3 = 0.0
        N = size(m, 1)
        for t = 1:N
            _x = m[t, 1] - mean_x
            _y = m[t, 2] - mean_y
            t1 += abs2(_x)
            t2 += _x * _y
            t3 += abs2(_y)
        end
        t1 /= N
        t2 /= N
        t3 /= N
        T[1, 1] = t1
        T[1, 2] = t2
        T[2, 1] = t2
        T[2, 2] = t3
        λ1, λ2 = eigvals(T)
        A = -log(1 - abs2(λ1 - λ2)/(2*abs2(λ1 + λ2)))
        push!(result, [Int(n), Int(N), A])
    end
    result
end

function non_gaussian_parameter(df, x, frame=:FRAME2)
    Δ = diff(df[!, x])
    prepend!(Δ, NaN)
    Δ[findall(x -> x == 1, df[!, frame])] .= NaN
    filter!(!isnan, Δ)
    a = mean(Δ.^4)
    b = mean(abs2.(Δ))^2
    a/3b - 1
end

function efficiency(df, id, x, y)
    result = DataFrame(TrackID = Int[], value = Float64[])
    @inbounds for n in 1:maximum(df[!, id])
        m = extract(df, n, id, [x, y])
        c = 0.0
        T = size(m, 1)
        for t in 1:T-1
            c += squared_displacement(m, t, 1)
        end
        c /= (T-1)
        c₀ = squared_displacement(m, 1, T-1)
        push!(result, [n, c₀/c])
    end
    result
end

function _largest_distance(m)
    T = size(m, 1)
    d = displacement(m, 1)
    for δ = 1:T-1
        for t = 1:T-δ
            dₜ = displacement(m, t, δ)
            if dₜ > d
                d = dₜ
            end
        end
    end
    d
end

function _total_length(m)
    T = size(m, 1)
    c = 0.0
    for t = 1:T-1
        c += displacement(m, t, 1)
    end
    c
end

function fractal_dimension(df, id, x, y)
    result = DataFrame(TrackID = Int[], value = Float64[])
    @inbounds @simd for n in minimum(df[!, id]):maximum(df[!, id])
        m = extract(df, n, id, [x, y])
        N = size(m, 1)-1
        d = _largest_distance(m)
        L = _total_length(m)
        push!(result, [n, log(N)/log(N*d/L)])
    end
    result
end

function gaussianity(
    df::DataFrame,
    id::Symbol = :TrackID,
    x::Symbol = :POSITION_X,
    y::Symbol = :POSITION_Y,
)
    gaussianity = DataFrame(TrackID = Int64[], gaussianity = Float64[], delta_t = Int64[], n = Int64[])
    @inbounds for n in sort(collect(Set(df[!, id])))
        m = extract(df, Int(n), id, [x, y])
        T = size(m, 1)
        for δ = 1:T-1
            r⁴ = 0.0
            r² = 0.0
            @simd for t = 1:(T-δ)
                r⁴ += abs2(squared_displacement(m, t, δ))
                r² += squared_displacement(m, t, δ)
            end
            r⁴ /= (T - δ)
            r² /= (T - δ)
            push!(gaussianity, [n, r⁴/(2*abs2(r²)) - 1, δ, T - δ])
        end
    end
    gaussianity
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
"""
    GenericRoutine{name} <: Function

A `GenericRoutine` cannot be directly execute on a quantum
device. It is a Julia `Function` that returns `Operation`,
and `Operation` can be executed on quanutm device.

!!! note
    An instance of `GenericRoutine` should be treated like
    `Function`.
"""
abstract type GenericRoutine{name} <: Function end

"""
    abstract type Routine end

This defines operations that one can execute on a quantum device.
"""
abstract type Routine end

"""
    abstract type IntrinsicRoutine <: Routine end

`IntrinsicRoutine` are the routine that can be executed
by the compile target that is pre-defined in the compiler.
"""
abstract type IntrinsicRoutine <: Routine end

"""
    Operation{P, Args} <: Routine

An `Operation` is a user defined composite routine
that can be called in other `Operation` or execute
on target device.
"""
struct Operation{P,Args} <: Routine
    parent::P
    args::Args
end

"""
    AdjointOperation{P} <: Routine

An `AdjointOperation` is the adjoint of another
`Routine`.
"""
struct AdjointOperation{P} <: Routine
    parent::P
end

Base.adjoint(x::Routine) = AdjointOperation(x)
Base.adjoint(x::AdjointOperation) = x.parent

routine_name(::Type) = nothing
routine_name(x) = routine_name(typeof(x))
routine_name(::Type{<:GenericRoutine{name}}) where {name} = name
routine_name(::Type{T}) where {T<:IntrinsicRoutine} = nameof(T)
routine_name(::Type{<:Operation{P}}) where {P} = routine_name(P)
routine_name(::Type{<:AdjointOperation{P}}) where {P} = Symbol("adjoint_", routine_name(P))

# this is only a place holder
# we are not using struct because
# singleton types will get const prop
struct MeasureResult{T}
    result::T
end

@inline function Base.:(==)(lhs::MeasureResult, rhs)
    measure_cmp(lhs, rhs)
end

@inline function Base.:(==)(lhs, rhs::MeasureResult)
    measure_cmp(lhs, rhs)
end

@inline function Base.:(==)(lhs::MeasureResult, rhs::MeasureResult)
    measure_cmp(lhs, rhs)
end

@noinline function measure_cmp(lhs, rhs)
    throw(IntrinsicError("cannot compare measurement result outside @device"))
end

struct Chain <: Routine
    args::Vector{Any}
    Chain(args...) = new(collect(Any, args))
end

# prop adjoint to leaves
Base.adjoint(c::Chain) = Chain(map(adjoint, c.args))

struct Typed
    obj
    type
end

struct Gate <: Routine
    operation # SSAValue/Routine
    locations # SSAValue/Locations
end

Base.adjoint(x::Gate) = Gate(adjoint(x.operation), x.locations)

struct Ctrl <: Routine
    gate::Gate
    ctrl # SSAValue/CtrlLocation
end

Base.adjoint(x::Ctrl) = Ctrl(adjoint(x.gate), x.ctrl)

struct BlockIR
    parent::IRCode
    nqubits::Int
    circuit::Chain
end

@as_record Chain
@as_record Gate
@as_record Ctrl
@as_record Typed
@as_record Operation
@as_record AdjointOperation

Base.length(c::Chain) = length(c.args)

function leaves(root::Chain)
    nodes = []
    sizehint!(nodes, length(root))
    for each in root.args
        if each isa Chain
            append!(nodes, leaves(each::Chain))
        else
            push!(nodes, each)
        end
    end
    return nodes
end

function Base.:(==)(lhs::Chain, rhs::Chain)
    mapreduce(==, &, lhs.args, rhs.args)
end

function Base.:(==)(lhs::Gate, rhs::Gate)
    lhs.operation == rhs.operation && lhs.locations == rhs.locations
end

function Base.:(==)(lhs::Ctrl, rhs::Ctrl)
    lhs.gate == rhs.gate && lhs.ctrl == rhs.ctrl
end


text:
@testsuite "mapreduce essentials" AT->begin
    @testset "mapreducedim! $ET" for ET in supported_eltypes() begin
        T = AT{ET}
        range = ET <: Real ? (ET(1):ET(10)) : ET
        for (sz,red) in [(10,)=>(1,), (10,10)=>(1,1), (10,10,10)=>(1,1,1), (10,10,10)=>(10,10,10),
                         (10,10,10)=>(1,10,10), (10,10,10)=>(10,1,10), (10,10,10)=>(10,10,1)]
            @test compare((A,R)->Base.mapreducedim!(identity, +, R, A), AT, rand(range, sz), zeros(ET, red))
            @test compare((A,R)->Base.mapreducedim!(identity, *, R, A), AT, rand(range, sz), ones(ET, red))
            @test compare((A,R)->Base.mapreducedim!(x->x+x, +, R, A), AT, rand(range, sz), zeros(ET, red))
        end

        # implicit singleton dimensions
        @test compare((A,R)->Base.mapreducedim!(identity, +, R, A), AT, rand(range, (2,2)), zeros(ET, (2,)))
        @test compare((A,R)->Base.mapreducedim!(identity, +, R, A), AT, rand(range, (2,3)), zeros(ET, (2,)))
    end
    end

    @testset "reducedim! $ET" for ET in supported_eltypes() begin
        T = AT{ET}
        range = ET <: Real ? (ET(1):ET(10)) : ET
        for (sz,red) in [(10,)=>(1,), (10,10)=>(1,1), (10,10,10)=>(1,1,1), (10,10,10)=>(10,10,10),
                         (10,10,10)=>(1,10,10), (10,10,10)=>(10,1,10), (10,10,10)=>(10,10,1)]
            @test compare((A,R)->Base.reducedim!(+, R, A), AT, rand(range, sz), zeros(ET, red))
            @test compare((A,R)->Base.reducedim!(*, R, A), AT, rand(range, sz), ones(ET, red))
        end
    end
    end

    @testset "mapreduce $ET" for ET in supported_eltypes() begin
        T = AT{ET}
        range = ET <: Real ? (ET(1):ET(10)) : ET
        for (sz,dims) in [(10,)=>[1], (10,10)=>[1,2], (10,10,10)=>[1,2,3], (10,10,10)=>[],
                          (10,)=>:, (10,10)=>:, (10,10,10)=>:,
                          (10,10,10)=>[1], (10,10,10)=>[2], (10,10,10)=>[3]]
            @test compare(A->mapreduce(identity, +, A; dims=dims, init=zero(ET)), AT, rand(range, sz))
            @test compare(A->mapreduce(identity, *, A; dims=dims, init=one(ET)), AT, rand(range, sz))
            @test compare(A->mapreduce(x->x+x, +, A; dims=dims, init=zero(ET)), AT, rand(range, sz))
        end
    end
    end

    @testset "reduce $ET" for ET in supported_eltypes() begin
        T = AT{ET}
        range = ET <: Real ? (ET(1):ET(10)) : ET
        for (sz,dims) in [(10,)=>[1], (10,10)=>[1,2], (10,10,10)=>[1,2,3], (10,10,10)=>[],
                          (10,)=>:, (10,10)=>:, (10,10,10)=>:,
                          (10,10,10)=>[1], (10,10,10)=>[2], (10,10,10)=>[3]]
            @test compare(A->reduce(+, A; dims=dims, init=zero(ET)), AT, rand(range, sz))
            @test compare(A->reduce(*, A; dims=dims, init=one(ET)), AT, rand(range, sz))
        end
    end
    end
end

@testsuite "mapreduce derivatives" AT->begin
    @testset "sum prod minimum maximum $ET" for ET in supported_eltypes() begin
        T = AT{ET}
        range = ET <: Real ? (ET(1):ET(10)) : ET
        for (sz,dims) in [(10,)=>[1], (10,10)=>[1,2], (10,10,10)=>[1,2,3], (10,10,10)=>[],
                          (10,)=>:, (10,10)=>:, (10,10,10)=>:,
                          (10,10,10)=>[1], (10,10,10)=>[2], (10,10,10)=>[3]]
            @test compare(A->sum(A), AT, rand(range, sz))
            @test compare(A->sum(abs, A), AT, rand(range, sz))
            @test compare(A->sum(A; dims=dims), AT, rand(range, sz))
            @test compare(A->prod(A), AT, rand(range, sz))
            @test compare(A->prod(abs, A), AT, rand(range, sz))
            @test compare(A->prod(A; dims=dims), AT, rand(range, sz))
            if !(ET <: Complex)
                @test compare(A->minimum(A), AT, rand(range, sz))
                @test compare(A->minimum(x->x*x, A), AT, rand(range, sz))
                @test compare(A->minimum(A; dims=dims), AT, rand(range, sz))
                @test compare(A->maximum(A), AT, rand(range, sz))
                @test compare(A->maximum(x->x*x, A), AT, rand(range, sz))
                @test compare(A->maximum(A; dims=dims), AT, rand(range, sz))
            end
        end
        OT = isbitstype(widen(ET)) ? widen(ET) : ET
        for (sz,red) in [(10,)=>(1,), (10,10)=>(1,1), (10,10,10)=>(1,1,1), (10,10,10)=>(10,10,10),
                         (10,10,10)=>(1,10,10), (10,10,10)=>(10,1,10), (10,10,10)=>(10,10,1)]
            if !(ET <: Complex)
                @test compare((A,R)->minimum!(R, A), AT, rand(range, sz), fill(typemax(ET), red))
                @test compare((A,R)->maximum!(R, A), AT, rand(range, sz), fill(typemin(ET), red))
            end
        end
        # smaller-scale test to avoid very large values and roundoff issues
        for (sz,red) in [(2,)=>(1,), (2,2)=>(1,1), (2,2,2)=>(1,1,1), (2,2,2)=>(2,2,2),
                         (2,2,2)=>(1,2,2), (2,2,2)=>(2,1,2), (2,2,2)=>(2,2,1)]
            @test compare((A,R)->sum!(R, A), AT, rand(range, sz), zeros(OT, red))
            @test compare((A,R)->prod!(R, A), AT, rand(range, sz), ones(OT, red))
        end
    end
    end

    @testset "any all count ==" begin
        for Ac in ([false, false], [false, true], [true, true],
                   [false false; false false], [false true; false false],
                   [true true; false false], [true true; true true])
            @test compare(A->any(A), AT, Ac)
            @test compare(A->all(A), AT, Ac)
        end
        for Ac in ([1, 1], [1, 2], [2, 2],
                   [1 1; 1 1], [1 2; 1 1],
                   [2 2; 1 1], [2 2; 2 2])
            @test compare(A->any(iseven, A), AT, Ac)
            @test compare(A->all(iseven, A), AT, Ac)
            @test compare(A->count(iseven, A), AT, Ac)

            A = AT(Ac)
            @test A == copy(A)
            @test A !== copy(A)
            @test A == deepcopy(A)
            @test A !== deepcopy(A)

            B = similar(A)
            @allowscalar B[1] = 3
            @test A != B
        end
    end
end

@testsuite "mapreduce (old tests)" AT->begin
    # old tests: can be removed, but left in here for a while to ensure the new impl works
    @testset "mapreduce" begin
        for ET in supported_eltypes()
            T = AT{ET}
            @testset "$ET" begin
                range = ET <: Integer ? (ET(-2):ET(2)) : ET
                @testset "mapreducedim" begin
                    for N in (2, 10)
                        @test compare(x -> sum(x, dims=2),      AT, rand(range, N, N))
                        @test compare(x -> sum(x, dims=1),      AT, rand(range, N, N))
                        @test compare(x -> sum(x, dims=(1, 2)), AT, rand(range, N, N))

                        @test compare(x -> sum(x, dims=2),      AT, rand(range, N, 10))
                        @test compare(x -> sum(x, dims=1),      AT, rand(range, N, 10))

                        @test compare(x -> sum(x, dims=2),      AT, rand(range, 10, N))
                        @test compare(x -> sum(x, dims=1),      AT, rand(range, 10, N))

                        _zero = zero(ET)
                        _addone(z) = z + one(z)
                        @test compare(x->mapreduce(_addone, +, x; dims = 2),
                                      AT, rand(range, N, N))
                        @test compare(x->mapreduce(_addone, +, x; dims = 2, init = _zero),
                                      AT, rand(range, N, N))

                        @test compare((x,y)->mapreduce(+, +, x, y; dims = 2),
                                      AT, rand(range, N, N), rand(range, N, N))
                        @test compare((x,y)->mapreduce(+, +, x, y; dims = 2, init = _zero),
                                      AT, rand(range, N, N), rand(range, N, N))
                    end
                end
                @testset "sum maximum minimum prod" begin
                    for dims in ((4048,), (77,), (1923,209))
                        @test compare(sum,  AT, rand(range, dims))
                        @test compare(prod, AT, rand(range, dims))
                        @test compare(x -> sum(abs, x),  AT, rand(range, dims))
                        @test compare(x -> prod(abs, x), AT, rand(range, dims))
                        @test compare(x -> sum(abs2, x),  AT, rand(range, dims))
                        @test compare(x -> prod(abs2, x), AT, rand(range, dims))
                        ET <: Complex || @test compare(maximum, AT,rand(range, dims))
                        ET <: Complex || @test compare(minimum, AT,rand(range, dims))
                    end
                end
                @testset "broadcasting behavior" begin
                        @test compare((x,y)->mapreduce(+, +, x, y), AT,
                                      rand(range, 1), rand(range, 2, 2))
                        @test compare(AT, rand(range, 1), rand(range, 2, 2)) do x, y
                            bc = Broadcast.instantiate(Broadcast.broadcasted(*, x, y))
                            reduce(+, bc)
                        end

                end
            end
        end
        @testset "any all ==" begin
            for Ac in ([false, false], [false, true], [true, true])
                A = AT(Ac)
                @test A isa AT{Bool,1}
                @test any(A) == any(Ac)
                @test all(A) == all(Ac)
                @test A == copy(A)
                @test A !== copy(A)
                @test A == deepcopy(A)
                @test A !== deepcopy(A)
            end
        end
    end
end


text:
@polly function kernel_trisolv(L, x, b)
    n = size(L,1)
    for i = 1:n
        x[i] = b[i]
        for j = 1:(i-1)
            x[i] -= L[i,j] * x[j]
        end
        x[i] = x[i] / L[i,i]
    end
end

let
    n = 2000

    L = zeros(Float32, n, n)
    x = zeros(Float32, n)
    b = zeros(Float32, n)

    for i = 1:n
        x[i] = -999
        b[i] = i
        for j = 1:i
            L[i,j] = (i+n-j+1)*2/n
        end
    end

    SUITE["trisolv"] = @benchmarkable kernel_trisolv(L, x, b) setup = (L = copy($L); x = copy($x); b = copy($b))
end


text:
#md # Needed packages
using DataFrames
using GraphRecipes
using Plots
using AbstractTrees
gr()

#md # this file holds the needed function to build a desition tree base on
#md # [Iterative Dichotomiser 3 method](http://hunch.net/~coms-4771/quinlan.pdf)
#md # this file documentation is base on [Weave.jl](https://github.com/JunoLab/Weave.jl) and [Literate](https://github.com/fredrikekre/Literate.jl)

#md # # ENTROPY

#md # this function calculate the entropy, at Base 2, of the input data for a given
#md # category.
#md #+
#md # base 2 entropy is good for binary class that has only two values.
#md # for more information see: [Entropy](https://en.wikipedia.org/wiki/Entropy)

#src """
#src    entropy2(data::DataFrame, categorycol::Symbol)
#src
#src this function calculate the category entropy using base 2.
#src this is good for binary class that can have only two values.
#src for more information see: [Entropy](https://en.wikipedia.org/wiki/Entropy)
#src the function output the category entropy
#src """

function entropy2(data::DataFrame, categorycol::Symbol)
    #src # find the total number of samples
    totalnumofsamples = length(data[!,categorycol])
    #src # get the category sets (there should be only 2)
    categorysets = Set(data[!,categorycol])
    #src # define the Dictionary for the Probability of each dataset to have each categoryset
    Pcategory = Dict{String,Real}()
    for categoryset in categorysets
        #src # find the number of samples of the choosen dataset
        Pcategory[string(categoryset)] = count(i->i == categoryset, data[!,categorycol])/totalnumofsamples
    end
    return -sum([v == 0 ? 0 : v*log2(v) for v in values(Pcategory)])
end # function entropy2

#md # # Best feature to split
#md # this function run on all the features in the DataFrame and find the best
#md # feature to split about, base on 'entropy2' function

#src """
#src    best_feature_to_split(data::DataFrame, features::Tuple{Symbol}, categorycol::Symbol)
#src
#src this function run on all the features in the DataFrame and find the best
#src feature to split about, base on 'entropy2' function
#src """
function best_feature_to_split(data::DataFrame, features::Tuple, categorycol::Symbol)
    baseline = entropy2(data,categorycol) # the base line is the entropy of the total data
    feature_entropy =  Dict{String,Real}()
    #src # this for loop run on each feature column in the data and calculate his entropy
    for feature in features
        feature_entropy[string(feature)] = 0.0
        # categorysets = Set(data[!,categorycol])
        featuresets = Set(data[!,feature])
        # for categoryset in categorysets
        # partitioned_data = filter(row -> row[categorycol] == categoryset,data)
        # porportion = length(partitioned_data[!,feature])/length(data[!,feature])
        # feature_entropy[string(feature)] += porportion * entropy2(partitioned_data, feature)
        for featureset in featuresets
            partitioned_data = filter(row -> row[feature] == featureset,data)
            porportion = length(partitioned_data[!,feature])/length(data[!,feature])
            feature_entropy[string(feature)] += porportion * entropy2(partitioned_data, categorycol)
        end
    end
    information_gain = Dict{String,Real}()
    for feature in features
        information_gain[string(feature)] = baseline - feature_entropy[string(feature)]
    end
    @info "information_gain: $information_gain"
    return findmax(information_gain)
end # function best_feature_to_split

#md # Potential leaf node
#md # this function find the most common category

#src """
#src    potential_leaf_node(data::DataFrame, categorycol::Symbol)
#src
#src this function find the most common category
#src """
function potential_leaf_node(data::DataFrame, categorycol::Symbol)
    counter = Dict{String,Real}()
    categorysets = Set(data[!,categorycol])
    for categoryset in categorysets
        counter[string(categoryset)] = count(i->i == categoryset, data[!,categorycol])
    end

    return findmax(counter)
end # function potential_leaf_node

#md # # Create Tree
#md # the most common is define as the category which has larger percentage then the 'persent_threshold'
#md # the default value is set to 100%

#src """
#src    create_tree(data::DataFrame, features::Tuple, categirycol::Symbol, persent_threshold::Int = 100)
#src
#src the most common is define as the category which has larger percentage then the 'persent_threshold'
#src the default value is set to 100%
#src """
function create_tree(data::DataFrame, features::Tuple, categorycol::Symbol, persent_threshold::Int = 100)
    counter , category = potential_leaf_node(data, categorycol)
    counter == length(data[!,categorycol]) && return category
    node = Dict{String,Dict}()
    feature_gain, feature_lable = best_feature_to_split(data, features, categorycol)
    node[feature_lable] = Dict()
    classes = Set(data[!,Symbol(feature_lable)])
    for c in classes
        partitioned_data = data[data[!,Symbol(feature_lable)] .== c,:]
        node[feature_lable][string(c)] = create_tree(partitioned_data, features, categorycol)
    end
    return node
end # function create_tree

# AbstractTrees.children(d::Dict) = [p for p in d]
# AbstractTrees.children(p::Pair) = AbstractTrees.children(p[2])
# function AbstractTrees.printnode(io::IO, p::Pair)
#     str = isempty(AbstractTrees.children(p[2])) ? string(p[1], ": ", p[2]) : string(p[1], ": ")
#     print(io, str)
# end
"""
    plottree(tree)

this function plot the tree
"""
function plottree(tree)
    default(size=(1000, 1000))
    plot(TreePlot(tree), method=:tree, fontsize=10, nodeshape=:rectangle, nodesize = 0.05)
end # function plottree

"""
    classify(tree::Dict, data)

this function classify data using the a tree
"""
function classify(tree::Dict,data::DataFrame)
    root = collect(keys(tree))[1]
    node = tree[root]
    res = nothing
    for k in keys(node)
        if string(data[!,Symbol(root)][1]) == k
            if node[k] isa Dict
                res = classify(node[k],data)
            else
                return node[k]
            end
        end
    end
    return res
end # function classify

#====== exmaples for building desition tree and calssify new data =============#
# desition_tree_data_weather = DataFrame( outlook = ["sunny","sunny","overcast","rain","rain","rain","overcast","sunny","sunny","rain","sunny","overcast","overcast","rain"],
#                                         temperature = ["hot","hot","hot","mild","cool","cool","cool","mild","cool","mild","mild","mild","hot","mild"],
#                                         humidity = ["high","high","high","high","normal","normal","normal","high","normal","normal","normal","high","normal","high"],
#                                         windy = [0,1,0,0,0,1,1,0,0,0,1,1,0,1],
#                                         class = [0,0,1,1,1,0,1,0,1,1,1,1,1,0])
# desition_tree_data_weather_features = (:outlook,:temperature,:humidity,:windy)
# desition_tree_data_weather_categorycol = :class
#
# tree = create_tree(desition_tree_data_weather,desition_tree_data_weather_features,desition_tree_data_weather_categorycol)
#
# plottree(tree)
#
# data = DataFrame( outlook = "rain", temperature = "hot", humidity = "normal", windy = 0)
#
# @show classify(tree,data)
#
# desition_tree_data_box2 = DataFrame(x = [0,-1,1,0,0], y = [0,0,0,-1,1], out = [0,1,1,1,1])
# desition_tree_data_box_features = (:x,:y)
# desition_tree_data_box_categorycol = :out
#
# tree = create_tree(desition_tree_data_box2,desition_tree_data_box_features,desition_tree_data_box_categorycol)
#
# plottree(tree)
#
# data = DataFrame( y = 1, x = 1)
#
# @show classify(tree,data)


text:
OS_NAME == :Windows && Pkg.add("FactCheck") # Hack for appveyor.

module LexiconTests

using Docile
@document

using Docile.Interface, Lexicon, FactCheck

import Lexicon: Query

include(joinpath(dirname(@__FILE__), "testcases.jl"))
Docile.Cache.register_module(joinpath(dirname(@__FILE__), "testcases.jl"))

import .TestCases: f, g, A, T, @m, K, S

include("facts/query-parsing.jl")
include("facts/querying-results.jl")
include("facts/rendering.jl")
include("facts/filtering.jl")

end


text:
module Murmur3

include("x86.jl")
include("x64.jl")

end


text:
import Gtk

function gtkwindow(name, w, h, closecb=nothing)
    c = Gtk.@Canvas()
    win = Gtk.@Window(c, name, w, h)
    if closecb !== nothing
        Gtk.on_signal_destroy(closecb, win)
    end
    show(c)
end

function display(c::Gtk.Canvas, pc::PlotContainer)
    c.draw = let bad=false
        function (_)
            bad && return
            ctx = getgc(c)
            set_source_rgb(ctx, 1, 1, 1)
            paint(ctx)
            try
                Winston.page_compose(pc, Gtk.cairo_surface(c))
            catch e
                bad = true
                isa(e, WinstonException) || rethrow(e)
                println("Winston: ", e.msg)
            end
        end
    end
    Gtk.draw(c)
end

gtkdestroy(c::Gtk.Canvas) = Gtk.destroy(Gtk.toplevel(c))

# JWN: copied the following from tk.jl, but I don't know
# what it does, so I can't make sure it works
#
#function get_context(c::Gtk.Canvas, pc::PlotContainer)
#    device = CairoRenderer(Gtk.cairo_surface(c))
#    ext_bbox = BoundingBox(0,width(c),0,height(c))
#    _get_context(device, ext_bbox, pc)
#end
#
#get_context(pc::PlotContainer) = get_context(curfig(_display), pc)


text:
# States.jl
# Define electronic and nuclear states

"""
Function to define diabatic states as gaussian functions of width σ
--------------
R_n = bond length
σ = width of states
--------------
phi_l = gaussian centered on left site
phi_r = gaussian centered on right site
"""
function diabatic_state(R_n::Float64,σ::Float64=0.998)
    gaussian = function (R) return exp(-R^2/(2*σ^2)) end
    gaussian2(R) = gaussian(R)^2
    A = riemann(gaussian2, -10,10,1000)
    A = 1/sqrt(A)
    phi(R) = A*gaussian(R-R_n)
    return phi
end

"""
Overlap of diabatic states to determine coupling strength
--------------
phi_l = gaussian centered on left site
phi_r = gaussian centered on right site
--------------
J_lr = coupling strength to be used in diabatic Hamiltonian
"""
function overlap_phi(phi_l::Function, phi_r::Function)
    overlap_phi(R) = conj(phi_l(R))*phi_r(R)
    J_lr = riemann(overlap_phi, -10,10,1000)
    return J_lr
end

"""
Function for the derivative of a gaussian function
    (to be used in NACV when the derivative of phi w.r.t R is needed)
--------------
R_n = bond length
σ = width of states
--------------
d_phi = gaussian derivative
"""
function diabatic_derivative(σ::Float64=0.374)
    dphi = function (R) return (R)*exp(-(R)^2/(2*σ^2))/σ^2 end
end

"""
Function to define the adiabatic states from the diabatic states and
transformation Matrix
------------
U_nk = Unitary transformation matrix
phi_l = diabatic state on left site
phi_r = diabatic state on right site
------------
psi_g = adiabatic lower state
phi_e = adiabatic upper state
"""
function adiabatic_states(phi_l::Function, phi_r::Function, U_nk::Array)
    psi_g = function (R) return U_nk[1]*phi_l(R) + U_nk[2]*phi_r(R) end
    psi_e = function (R) return U_nk[3]*phi_l(R) + U_nk[4]*phi_r(R) end
    return psi_g, psi_e
end

"""
Function to define the non-adiabatic coupling vectors for all possible transisions
--------------
phi_l = diabatic state on left site
phi_r = diabatic state on right site
dphi = generic derivative of gaussian function
U_nk = Unitary transformation matrix
R_n = bond length
dR_l = change of nuclear coordinate on the left site
dR_r = change of nuclear coordinate on the right site
dt = time step this change has occured in
--------------
d_ge = Non-adiabatic coupling vector from lower to upper adiabatic states
d_eg = Non-adiabatic coupling vector from upper to lower adiabatic states
d_gg = Non-adiabatic coupling vector from lower to lower adiabatic states
d_ee = Non-adiabatic coupling vector from upper to upper adiabatic states
"""
function NACV(phi_l::Function, phi_r::Function, dphi::Function, U_nk::Array, R_n::Float64, dR_l, dR_r, dt)
    dphi_l = function (R) return dphi((R)) end
    dphi_r = function (R) return dphi((R-R_n)) end

    # perform integrations for each pair of diabatic states
    int_ll = function (R) return conj(phi_l(R))*dphi_l(R) end
    ll = riemann(int_ll,-10,10,1000)
    int_lr = function (R) return conj(phi_l(R))*dphi_r(R) end
    lr = riemann(int_lr,-10,10,1000)
    int_rl = function (R) return conj(phi_r(R))*dphi_l(R) end
    rl = riemann(int_rl,-10,10,1000)
    int_rr = function (R) return conj(phi_r(R))*dphi_r(R) end
    rr = riemann(int_rr,-10,10,1000)

    l_ge = dR_l*(conj(U_nk[1]*ll + U_nk[2]*rl)*U_nk[3])/dt
    r_ge = dR_r*(conj(U_nk[1]*lr + U_nk[2]*rr)*U_nk[4])/dt
    d_ge = l_ge+r_ge

    l_eg = dR_l*(conj(U_nk[3]*ll + U_nk[4]*rl)*U_nk[1])/dt
    r_eg = dR_r*(conj(U_nk[3]*lr + U_nk[4]*rr)*U_nk[2])/dt
    d_eg = l_eg+r_eg

    l_gg = dR_l*(conj(U_nk[1]*ll + U_nk[2]*rl)*U_nk[1])/dt
    r_gg = dR_r*(conj(U_nk[1]*lr + U_nk[2]*rr)*U_nk[2])/dt
    d_gg = l_gg+r_gg

    l_ee = dR_l*(conj(U_nk[3]*ll + U_nk[4]*rl)*U_nk[3])/dt
    r_ee = dR_r*(conj(U_nk[3]*lr + U_nk[4]*rr)*U_nk[4])/dt
    d_ee = l_ee+r_ee

    return d_ge,d_eg,d_gg,d_ee, [l_ge*dt/dR_l, r_ge*dt/dR_r, l_eg*dt/dR_l, r_eg*dt/dR_r]
end

"""
Function to determine adiabatic coefficients from diabatic coefficients
    and unitary transformation matrix
-------------
U_nk = Unitary Transformation Matrix
c1 = diabatic coefficient of left site
c2 = diabatic coefficient of right site
-------------
ag = adiabatic coefficient of ground state
ae = adiabatic coefficient of excited state
"""
function a_mn(U_nk::Array,cl::Complex{Float64},cr::Complex{Float64})
    ag = (cl*U_nk[4] - cr*U_nk[3])/(U_nk[1]*U_nk[4]-U_nk[2]*U_nk[3])
    ae = (cl*U_nk[2] - cr*U_nk[1])/(U_nk[2]*U_nk[3]-U_nk[1]*U_nk[4])
    return ag,ae
end

"""
Function to determine the probability of switching states
------------
a_m = adiabatic coefficient of state m
a_n = adiabatic coefficient of state n
d_mn = Non-adiabatic coupling vector for transition from state n to state m
dt = time step
------------
g_mn = probability of a transition from adiabatic surface m to adiabatic surface n
"""
function g_mn(a_m::Complex{Float64}, a_n::Complex{Float64}, d_mn::Float64, dt::Float64)
    g_mn =  -2*real(dt*(conj(a_n)*a_m*d_mn)/(conj(a_m)*a_m))
    return g_mn
end


text:
mg = Multigraph(3)
add_vertices!(mg, 3)
rem_vertices!(mg, [1, 3])
add_edge!(mg, 2, 5)
add_edge!(mg, 2, 4, 2)

@test outneighbors(mg, 2) == [4, 5]
eit = edges(mg)
iterate(eit)[1]
vertices(mg)
@test iterate(eit)[2] == (1, 2)
mes = [me for me in edges(mg)]
@test length(mes) == length(eit)


text:
using Plots, QuadGK, ColorSchemes, ForwardDiff

# Vector e.g. [1, 4, 5] or collected Range e.g. collect(1:5)
inputstates = [0, 2, 3]
sort!(inputstates)

function constant(n)
    1/sqrt(size(inputstates, 1))
    #2/(n*pi)*(1-cos(n*pi/2))
end

function calculate(states)
    dx = 0.25
    dt = 1/30
    elapse = 10
    m = 0.2
    w = 1

    h = 1.054571817
    A = sqrt(2h*(states[end]+0.5)/w/m)+4
    xran = (-A):dx:A
    yran = 0:dx:elapse
    clrs = ColorSchemes.GnBu_9
    grad1 = cgrad(ColorSchemes.rainbow)

    Cₙ = [constant(i) for i in states]

    function GenerateHermite(n)
        Hermite=Function[]

        push!(Hermite,x->1);
        push!(Hermite,x->2*x);

        for ni in 3:n
            push!(Hermite,x->2*x*Hermite[ni-1](x)-2*(ni-1)*Hermite[ni-2](x))
        end
        return Hermite
    end

    Hermite=GenerateHermite(states[end]+1)

    Ψₙ(x) = [(m*w/pi/h)^(1/4)*Hermite[i+1](sqrt(m*w/h)*x)*exp(-m*w*x^2/2/h)/sqrt(2^i*factorial(big(i))) for i in states]
    dΨₙ(x) = ForwardDiff.derivative(Ψₙ, x)

    φₙ(t) = [exp(-im*w*t*(i+0.5)) for i in states]

    radii = [[0,abs(i)] for i in Cₙ]

    E(n)=h*w*(n+0.5)

    Eₙ = [E(i) for i in states]

    anim = @animate for ti in 0:dt:elapse

            Ψ(x) = sum(Ψₙ(x).*φₙ(ti).*Cₙ)
            dΨ(x) = sum(dΨₙ(x).*φₙ(ti).*Cₙ)

            ref = real.(Ψ.(xran))
            imf = imag.(Ψ.(xran))
            arg = vcat(atan.(ref,imf), pi, -pi)
            ρ(x) = abs2(Ψ(x))
            V(x) = m*w^2*x^2/2
            J(x) = real(im*h/(2*m)*(Ψ(x)*conj(dΨ(x))-conj(Ψ(x))*dΨ(x)))
            xΨ(x) = x*ρ(x)
            x²Ψ(x) = x^2*ρ(x)
            #expcX = first(quadgk(xΨ, -Inf, Inf))
            #expcX² = first(quadgk(x²Ψ, -Inf, Inf))
            #σₓ = sqrt(expcX²-expcX^2)

            theta = [[0,-i] for i in angle.(φₙ(ti))]

            p1 = plot(xran, imf, ref, ylims = [-1,1], zlims = [-1,1], lc = grad1, labels = false, title = "Ψ(x, t)", zlabel = "Re(Ψ)", ylabel = "Im(Ψ)", line_z = arg, legend = false, camera = [20sin(ti)+55, 25])
            p2 = plot(xran, ρ, ylims = [0,1], color = clrs[6], labels = false, title = "∣Ψ(x, t)∣²", fill = true)
            p3 = plot(theta, radii, proj =:polar, lims = [0,1], showaxis = false, labels = false, palette = [clrs[4], clrs[5], clrs[6], clrs[7], clrs[8]], title = "φ(t)")
            #scatter!(p2, [expcX], [0.5], color = clrs[5], xerror = σₓ, markerstrokecolor=clrs[8], labels = false)
            plot!(p1, xran, zeros(size(xran,1),1), zeros(size(xran,1),1), label = false, alpha = 0.1, color = :black)
            plot!(p2, V, label = "V(x)", color = clrs[5], linestyle=:dash)

            p5 = plot(xran, J, ylims = [-A, A], labels = false, title = "J(x,t)", color = clrs[6])

            l1 = @layout[a b; c{0.4h} d]

            plot(p1, p2, p3, p5, fontfamily = :serif, layout = l1)

    end

    Eplot = vcat(E(states[begin]-2), E(states[begin]-1), Eₙ, E(states[end]+1), E(states[end]+2))
    Cplot = vcat(0, 0, abs2.(Cₙ), 0, 0)
    Psi(x, t) = abs2(sum(Ψₙ(x).*φₙ(t).*Cₙ))
    l2 = @layout[a{0.3w} [b; c]]

    p4 = plot(Cplot, Eplot, xlims = [0, 1], markershapes = :circle, ylabel = "E", xlabel = "∣C∣²", labels = false, line = clrs[4], markercolors = clrs[4])
    hline!(p4, [sum(abs2.(Cₙ).*Eₙ)], color = clrs[9], label = false)
    p8 = surface(xran, yran, Psi, camera = [15, 75], seriescolor = cgrad(clrs), ylabel = "t", zlabel = "∣Ψ(x, t)∣²")
    p6 = contourf(xran, yran, Psi, seriescolor = cgrad(clrs), ylabel = "t")

    stats = plot(p4, p8, p6, layout = l2, fontfamily = :serif)

    savefig(stats, "qho1_$inputstates _stats.png")

    gif(anim, "qho1_$inputstates.gif", fps = 1/dt)
end

@time calculate(inputstates)


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

@doc raw"""ReplicaSetStatus represents the current status of a ReplicaSet.

    IoK8sApiAppsV1beta2ReplicaSetStatus(;
        availableReplicas=nothing,
        conditions=nothing,
        fullyLabeledReplicas=nothing,
        observedGeneration=nothing,
        readyReplicas=nothing,
        replicas=nothing,
    )

    - availableReplicas::Int32 : The number of available replicas (ready for at least minReadySeconds) for this replica set.
    - conditions::Vector{IoK8sApiAppsV1beta2ReplicaSetCondition} : Represents the latest available observations of a replica set's current state.
    - fullyLabeledReplicas::Int32 : The number of pods that have labels matching the labels of the pod template of the replicaset.
    - observedGeneration::Int64 : ObservedGeneration reflects the generation of the most recently observed ReplicaSet.
    - readyReplicas::Int32 : The number of ready replicas for this replica set.
    - replicas::Int32 : Replicas is the most recently oberved number of replicas. More info: https://kubernetes.io/docs/concepts/workloads/controllers/replicationcontroller/#what-is-a-replicationcontroller
"""
mutable struct IoK8sApiAppsV1beta2ReplicaSetStatus <: SwaggerModel
    availableReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: availableReplicas
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiAppsV1beta2ReplicaSetCondition} } # spec name: conditions
    fullyLabeledReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: fullyLabeledReplicas
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration
    readyReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: readyReplicas
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas

    function IoK8sApiAppsV1beta2ReplicaSetStatus(;availableReplicas=nothing, conditions=nothing, fullyLabeledReplicas=nothing, observedGeneration=nothing, readyReplicas=nothing, replicas=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2ReplicaSetStatus, Symbol("availableReplicas"), availableReplicas)
        setfield!(o, Symbol("availableReplicas"), availableReplicas)
        validate_property(IoK8sApiAppsV1beta2ReplicaSetStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiAppsV1beta2ReplicaSetStatus, Symbol("fullyLabeledReplicas"), fullyLabeledReplicas)
        setfield!(o, Symbol("fullyLabeledReplicas"), fullyLabeledReplicas)
        validate_property(IoK8sApiAppsV1beta2ReplicaSetStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        validate_property(IoK8sApiAppsV1beta2ReplicaSetStatus, Symbol("readyReplicas"), readyReplicas)
        setfield!(o, Symbol("readyReplicas"), readyReplicas)
        validate_property(IoK8sApiAppsV1beta2ReplicaSetStatus, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        o
    end
end # type IoK8sApiAppsV1beta2ReplicaSetStatus

const _property_map_IoK8sApiAppsV1beta2ReplicaSetStatus = Dict{Symbol,Symbol}(Symbol("availableReplicas")=>Symbol("availableReplicas"), Symbol("conditions")=>Symbol("conditions"), Symbol("fullyLabeledReplicas")=>Symbol("fullyLabeledReplicas"), Symbol("observedGeneration")=>Symbol("observedGeneration"), Symbol("readyReplicas")=>Symbol("readyReplicas"), Symbol("replicas")=>Symbol("replicas"))
const _property_types_IoK8sApiAppsV1beta2ReplicaSetStatus = Dict{Symbol,String}(Symbol("availableReplicas")=>"Int32", Symbol("conditions")=>"Vector{IoK8sApiAppsV1beta2ReplicaSetCondition}", Symbol("fullyLabeledReplicas")=>"Int32", Symbol("observedGeneration")=>"Int64", Symbol("readyReplicas")=>"Int32", Symbol("replicas")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2ReplicaSetStatus }) = collect(keys(_property_map_IoK8sApiAppsV1beta2ReplicaSetStatus))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2ReplicaSetStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta2ReplicaSetStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2ReplicaSetStatus }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2ReplicaSetStatus[property_name]

function check_required(o::IoK8sApiAppsV1beta2ReplicaSetStatus)
    (getproperty(o, Symbol("replicas")) === nothing) && (return false)
    true
end

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


text:
#=
"""
    psfget(psf::BAMF.PSF)

psfget returns the contents of a PSF regardless of type
"""
function psfget(psf::BAMF.PSF_airy2D)
    return psf.ν
end

function psfget(psf::BAMF.PSF_gauss2D)
    return psf.σ
end
 =#
"""
    gendatastate(;seed::Integer=-1, psf::PSF=PSF_airy2D(.2*pi), n::Int32=6)

    Takes an optional Integer seed, set to random if not opted for
    An optional Int32 number of emitters n, defaulted to 6
    optional integer mean emitters per photon μ, defaulted to 1000
    optional Int32 pixel number sz, defaulted to 32
    optional Float32 intensity distribution shape parameter α, defaulted to 4f0
    optional Float32 background bg, defaulted to 1f-6
    Also takes an optional PSF, otherwise uses a PSF_airy2D with a v of .2*pi
    optional Float32 intensity standard deviation istd, defaulted to 10f0
    optional Float32 bndpixels, defaulted to -20f0

Outputs a StateFlatBg datastate, as well as all the information needed for an RJStruct sans the BAMFData object as a tuple, in the following order: datastate, psf, sz, xystd, istd, split_std, bndpixels, prior_photons
"""
function gendatastate(seed::Int32=Int32(-1); psf::PSF.PSF=PSF.Airy2D(1.2,0.6,0.1), n::Int32=Int32(6), μ::Int64=1000, sz::Int32=Int32(32), α::Float32=Float32(4), bg::Float32=1f-6, istd::Float32=10f0, bndpixels::Float32=-20f0)
    if seed == -1
        Random.seed!()
    else
        Random.seed!(seed)
    end

    σ=.42f0*pi/psf.ν
    # setup prior distribution on intensity
    θ=Float32(μ/α)
    g=Gamma(α,θ)
    len=Int32(1024)
    θ_start=Float32(1)
    θ_step=Float32(5.0)
    pdf_x=range(θ_start,step=θ_step,length=len)
    mypdf=pdf(g,pdf_x)
    prior_photons=BAMF.RJPrior(len,θ_start,θ_step,mypdf)

    # set emitter positions and intensity
    x=sz/2f0*ones(Float32,n)+2*σ*randn(Float32,n)
    y=sz/2f0*ones(Float32,n)+ 2*σ*randn(Float32,n)
    photons=Float32.(rand(g,n))
    datastate=BAMF.StateFlatBg(Int64(n),convert(Array{Float32},x),convert(Array{Float32},y),photons,bg)

    ## create a BAMF-type RJMCMC structure
    xystd=σ/10
    split_std=σ/2

    return datastate, psf, xystd, istd, split_std, bndpixels, prior_photons
end

"""
    genBAMFDD(T::Type{BAMF.BAMFData}, sz::Int32)

generates an empty BAMFData direct detection model of type T and size sz. The default size is 32.
"""
function genBAMFDD(T::Type{BAMF.ArrayDD}, sz::Int32=Int32(32))
    return BAMF.ArrayDD(sz)
end

#= function genBAMFDD(T::Type{BAMF.DataSLIVER}, sz::Int32=Int32(32))
    return BAMF.DataSLIVER(sz, [Int32(1)])
end

function genBAMFDD(T::Type{BAMF.AdaptData}, sz::Int32=Int32(32))
    return BAMF.AdaptData(sz, [(BAMF.DDMeasType, ())])
end =#

"""
    DDinfo(DD::BAMF.BAMFData)

returns a tuple of the information that ought to be contained in a direct detection BAMFData structure: size, an array of the appropriate dimensions, and the object itself.
"""
function DDinfo(DD::BAMF.BAMFData)
    dim= size(DD.data)
    sz= DD.sz
    return sz, dim, DD
end
#=
"""
    genBAMFSLIVER(T::Type{BAMF.BAMFData}, sz::Int32)

returns an empty BAMFData SLIVER model of type T and size sz. The default size is 32.
"""
function genBAMFSLIVER(T::Type{BAMF.DataSLIVER}, sz::Int32=Int32(32))
    return BAMF.DataSLIVER(sz, [Int32(2)])
end

function genBAMFSLIVER(T::Type{BAMF.AdaptData}, sz::Int32=Int32(32))
    return BAMF.AdaptData(sz, [(BAMF.SLIVERMeasType, (Float32(0f0), Float32(0f0)))])
end

"""
    SLIVERinfo(SLIVER::BAMF.BAMFData)

returns a tuple of the information that ought to be contained in a SLIVER BAMFData structure: size, an array of the appropriate dimensions, the object itself, and the x y inversion points alongside the integration time.
"""
function SLIVERinfo(SLIVER::BAMF.DataSLIVER)
    dim= size(SLIVER.data)
    sz= SLIVER.sz
    inttime, invx, invy= SLIVER.inttime[1], SLIVER.invx[1], SLIVER.invy[1]
    return sz, dim, inttime, invx, invy, SLIVER
end

function SLIVERinfo(SLIVER::BAMF.AdaptData)
    dim= size(SLIVER.data)
    sz= SLIVER.sz
    inttime, invx, invy= SLIVER.meastypes[1].inttime, SLIVER.meastypes[1].invx, SLIVER.meastypes[1].invy
    return sz, dim, inttime, invx, invy, SLIVER
end
 =#
"""
    genRJMCMC(burnin::Int32, iterations::Int32)

a function that produces an RJMCMCStruct with the usual DD accept and propose functions alongside some standard jump probabilities. The default burnin and iterations is 1000.
"""
function genRJMCMC(burnin::Int32=Int32(1000), iterations::Int32=Int32(1000))
    ## setup the RJMCMC.jl model
    # Jumptypes are: move, bg, add, remove, split, merge
    njumptypes=6
    jumpprobability=[1,0,.1,.1,.1,.1] # Model with no bg
    jumpprobability=jumpprobability/sum(jumpprobability)

    # create an RJMCMC structure with all model info
    acceptfuns=[BAMF.accept_move,BAMF.accept_bg,BAMF.accept_add,BAMF.accept_remove,BAMF.accept_split,BAMF.accept_merge] #array of functions
    propfuns=[BAMF.propose_move,BAMF.propose_bg,BAMF.propose_add,BAMF.propose_remove,BAMF.propose_split,BAMF.propose_merge] #array of functions
    myRJMCMC=RJMCMC.RJMCMCStruct(burnin,iterations,njumptypes,jumpprobability,propfuns,acceptfuns)
    return myRJMCMC
end


text:
"""
    struct DynamicDistributionModel{Tf<:Function} <: AbstractDistributionModel
        f::Tf
        ns::Int
        na::Int
    end

Using a general function `f` to store the transformations.
"""
struct DynamicDistributionModel{Tf<:Function} <: AbstractDistributionModel
    f::Tf
    ns::Int
    na::Int
end

function update!(::DynamicDistributionModel, ::AbstractTurnBuffer) end

getstates(m::DynamicDistributionModel) = m.ns
getactions(m::DynamicDistributionModel) = m.na
(m::DynamicDistributionModel)(s, a) = m.f(s, a)


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

@doc raw"""LimitedPriorityLevelConfiguration specifies how to handle requests that are subject to limits. It addresses two issues:  * How are requests for this priority level limited?  * What should be done with requests that exceed the limit?

    IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration(;
        assuredConcurrencyShares=nothing,
        limitResponse=nothing,
    )

    - assuredConcurrencyShares::Int32 : `assuredConcurrencyShares` (ACS) configures the execution limit, which is a limit on the number of requests of this priority level that may be exeucting at a given time.  ACS must be a positive number. The server's concurrency limit (SCL) is divided among the concurrency-controlled priority levels in proportion to their assured concurrency shares. This produces the assured concurrency value (ACV) --- the number of requests that may be executing at a time --- for each such priority level:              ACV(l) = ceil( SCL * ACS(l) / ( sum[priority levels k] ACS(k) ) )  bigger numbers of ACS mean more reserved concurrent requests (at the expense of every other PL). This field has a default value of 30.
    - limitResponse::IoK8sApiFlowcontrolV1alpha1LimitResponse : `limitResponse` indicates what to do with requests that can not be executed right now
"""
mutable struct IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration <: SwaggerModel
    assuredConcurrencyShares::Any # spec type: Union{ Nothing, Int32 } # spec name: assuredConcurrencyShares
    limitResponse::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1LimitResponse } # spec name: limitResponse

    function IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration(;assuredConcurrencyShares=nothing, limitResponse=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration, Symbol("assuredConcurrencyShares"), assuredConcurrencyShares)
        setfield!(o, Symbol("assuredConcurrencyShares"), assuredConcurrencyShares)
        validate_property(IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration, Symbol("limitResponse"), limitResponse)
        setfield!(o, Symbol("limitResponse"), limitResponse)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration

const _property_map_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration = Dict{Symbol,Symbol}(Symbol("assuredConcurrencyShares")=>Symbol("assuredConcurrencyShares"), Symbol("limitResponse")=>Symbol("limitResponse"))
const _property_types_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration = Dict{Symbol,String}(Symbol("assuredConcurrencyShares")=>"Int32", Symbol("limitResponse")=>"IoK8sApiFlowcontrolV1alpha1LimitResponse")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration)
    true
end

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


text:
module IssueReporter

using Pkg, Pkg.TOML, GitHub, URIParser, Documenter, DocStringExtensions

function generalregistrypath() :: String
  for i in DEPOT_PATH
    if isdir(joinpath(i, "registries", "General"))
      return joinpath(i, "registries", "General")
    end
  end

  ""
end

function generalregistry() :: Vector{Dict{String,Any}}
  if ! isempty(generalregistrypath())
    TOML.parsefile(joinpath(generalregistrypath(), "Registry.toml"))["packages"] |> values |> collect
  else
     Dict{String,Any}[]
   end
end

function searchregistry(pkgname::String) :: Dict{String,Any}
  for item in generalregistry()
    item["name"] == pkgname && return item
  end

  Dict{String,Any}()
end

"""
$(SIGNATURES)

Takes the name of a registered Julia package and returns the associated repo git URL.

Examples
```jldoctest
julia> IssueReporter.packageuri("IssueReporter")
"git://github.com/essenciary/IssueReporter.jl.git"
```
"""
function packageuri(pkgname::String) :: String
  pkg = searchregistry(pkgname)
  isempty(pkg) && return ""
  get!(TOML.parsefile(joinpath(generalregistrypath(), pkg["path"], "Package.toml")), "repo", "")
end

"""
$(SIGNATURES)

Checks if the required GitHub authentication token is defined.
"""
function tokenisdefined() :: Bool
  if ! haskey(ENV, "GITHUB_ACCESS_TOKEN")
    secrets_path = joinpath(@__DIR__, "secrets.jl")
    isfile(secrets_path) && include(secrets_path)
    haskey(ENV, "GITHUB_ACCESS_TOKEN") || return false
  end

  true
end

"""
$(SIGNATURES)

Returns the configured GitHub authentication token, if defined -- or throws an error otherwise.
"""
function token() :: String
  tokenisdefined() && return ENV["GITHUB_ACCESS_TOKEN"]
  error("""ENV["GITHUB_ACCESS_TOKEN"] is not set -- please make sure it's passed as a command line argument or defined in the `secrets.jl` file.""")
end

"""
$(SIGNATURES)

Performs GitHub authentication and returns the OAuth2 object, required by further GitHub API calls.
"""
function githubauth()
  token() |> GitHub.authenticate
end

"""
$(SIGNATURES)

Converts a registered Julia package name to the corresponding GitHub "username/repo_name" string.

Examples
```jldoctest
julia> IssueReporter.repoid("IssueReporter")
"essenciary/IssueReporter.jl"
```
"""
function repoid(package_name::String)
  pkg_url = packageuri(package_name) |> URIParser.parse_url
  repo_info = endswith(pkg_url.path, ".git") ? replace(pkg_url.path, r".git$"=>"") : pkg_url.path

  repo_info[2:end]
end

"""
$(SIGNATURES)

Creates a new GitHub issue with the title `title` and the content `body` onto the repo corresponding to the registered package called `pack-age_name`.
"""
function report(package_name::String, title::String, body::String)
  # GitHub.create_issue(repoid(package_name), auth = githubauth(), params = Dict(:title => title, :body => body))
  GitHub.create_issue("essenciary/julia-by-example-test-repo", auth = github_auth(), params = Dict(:title => title, :body => body))
end

end # module

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
# using Compat
using LinearOperators

"Abstract  for statistics returned by a solver"
abstract type KrylovStats end

"Type for statistics returned by non-Lanczos solvers"
mutable struct SimpleStats <: KrylovStats
  solved :: Bool
  inconsistent :: Bool
  residuals :: Array{Float64,1}
  Aresiduals :: Array{Float64,1}
  status :: String
end

# A standard implementation of the Conjugate Gradient method.
# The only non-standard point about it is that it does not check
# that the operator is definite.
# It is possible to check that the system is inconsistent by
# monitoring ‖p‖, which would cost an extra norm computation per
# iteration.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Salt Lake City, UT, March 2015.
#
# Try to adapt a stopping criterion based on regulα for the ARCq.
#
# First reformulate the TR case using the characteristic function to confirm it is
# equivalent to the usual implementation. A must be symmetric to define a quadratic objective
# q(x) = 0.5*x'*A*x - b'*x
#
#   JPD july 21 2016, Sherbrooke

include("krylov_aux.jl")

export cgARC
# Methods for various argument types.
#include("cg_methods.jl")
# cgARC{TA <: Number, Tb <: Number}(A :: Array{TA,2}, b :: Array{Tb,1};
#                                atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0, regulα :: Float64=1.0, verbose :: Bool=false)
cgARC(A::Array{TA, 2}, b::Array{Tb, 1}; atol::Float64 = 1e-08, rtol::Float64 = 1e-06, itmax::Int = 0, regulα::Float64 = 1, verbose::Bool = false) where {TA <: Number, Tb <: Number}  = cgARC(LinearOperator(A), b, atol=atol, rtol=rtol, itmax=itmax, regulα=regulα, verbose=verbose);

# cgARC{TA <: Number, Tb <: Number, IA <: Integer}(A :: SparseMatrixCSC{TA,IA}, b :: Array{Tb,1};
#                                               atol :: Float64=1.0e-8, rtol ::
#                                               Float64=1.0e-6, itmax :: Int=0, regulα :: Float64=1.0,  verbose :: Bool=false) =
cgARC(A::SparseMatrixCSC{TA, IA}, b::Array{Tb, 1}; atol::Float64 = 1e-08, rtol::Float64 = 1e-06, itmax::Int = 0, regulα::Float64 = 1, verbose::Bool = false) where {TA <: Number, Tb <: Number, IA <: Integer}  = cgARC(LinearOperator(A), b, atol=atol, rtol=rtol, itmax=itmax, regulα=regulα, verbose=verbose);

# function cgARC{T <: Real}(A :: LinearOperator, b :: Array{T,1};
#                        atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0,
#                        regulα :: Float64=1.0, verbose :: Bool=false)
"""The conjugate gradient method to solve the symmetric linear system Ax=b.

The method does _not_ abort if A is not definite.
"""
function cgARC(A::LinearOperator, b::Array{T, 1}; atol::Float64 = 1e-08, rtol::Float64 = 1e-06, itmax::Int = 0, regulα::Float64 = 1, verbose::Bool = false) where T <: Real
  n = size(b, 1);
  (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size");
  #isequal(triu(A)',tril(A)) || error("Must supply Hermitian matrix")
  regulα > 0.0  ||  error("Regularization must be strictly positive")

  verbose && @printf("CG: system of %d equations in %d variables\n", n, n);

  # Initial state.
  x = zeros(n)
  x̂ = copy(x)

  γ = dot(b, b);
  γ == 0 && return x;
  r = copy(b);
  p = copy(r);

  σ = 0.0

  iter = 0;
  itmax == 0 && (itmax = 2 * n);

  rNorm = sqrt(γ);
  rNorms = [rNorm;];
  ε = atol + rtol * rNorm;
  verbose && @printf("%5d  %8.1e ", iter, rNorm)

  solved = rNorm <= ε;
  tired = iter >= itmax;
  on_boundary = false;
  status = "unknown";

  q = s ->  0.5 * dot(s, copy(A * s)) - dot(b, s)

  m = s ->  q(s) + norm(s)^3/(3*regulα)
  #hO = α -> m(x+α*p)

  while ! (solved || tired)
    Ap = copy(A * p);  # Bug in LinearOperators? A side effect spoils the computation without the copy.
    pAp = BLAS.dot(n, p, 1, Ap, 1);

    α = γ / pAp;

    # Compute step size to the min of the regularized model.
    σ = to_minimum(A, b, x, p, Ap, pAp, α, regulα)

    verbose && @printf("%8.1e  %7.1e  %7.1e\n", pAp, α, σ);

    # Move along p from x to the min if either
    # the next step leads farther than the min of the regularized model or
    # we have nonpositive curvature.
    if ((pAp <= 0.0) | (α > (σ/2)))
      α = σ
      on_boundary = true
      verbose &&println("at minimum after ",iter," CG iterations.")
    end

    verbose && @printf("    %8.1e  %7.1e  %7.1e\n", pAp, α, σ);
    BLAS.axpy!(n,  α,  p, 1, x, 1);  # Faster than x = x + σ * p;
    BLAS.axpy!(n, -α, Ap, 1, r, 1);  # Faster than r = r - α * Ap;
    γ_next = BLAS.dot(n, r, 1, r, 1);
    rNorm = sqrt(γ);
    push!(rNorms, rNorm);

    solved = (rNorm <= ε) | on_boundary;
    tired = iter >= itmax;

    if !solved
      β = γ_next / γ;
      γ = γ_next;
      BLAS.scal!(n, β, p, 1)
      BLAS.axpy!(n, 1.0, r, 1, p, 1);  # Faster than p = r + β * p;
    end
    iter = iter + 1;
    verbose && @printf("%5d  %8.1e ", iter, rNorm);
  end
  verbose && @printf("\n");

  status = on_boundary ? "at min of the cubic regularization" : (tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol")
  stats = SimpleStats(solved, false, rNorms, [], status);
  return (x, stats);
end


text:
#=
    Model predictive control controller
=#

mutable struct OLFCOptions
    solver
    reducer::AbstractScenariosReducer
    generator::AbstractScenariosGenerator
    horizon::Int64
    nscenarios::Int64
    seasonal_targets::Union{Array{Float64,2}, Nothing}
    final_cost::Vector{Float64}

    OLFCOptions(; solver = CPLEX,
                  reducer = FeatureBasedReducer(),
                  generator = MarkovGenerator(),
                  horizon = 24,
                  nscenarios = 1,
                  seasonal_targets = nothing,
                  final_cost = [1e-3, 1e-3, 1e-2]) = new(solver, reducer, generator, horizon, nscenarios, seasonal_targets, final_cost)
end

mutable struct OLFC <: AbstractController
    options::OLFCOptions
    generations::Vector{Float64}
    storages::Vector{Float64}
    converters::Vector{Float64}
    decisions::NamedTuple
    model::JuMP.Model
    history::AbstractScenarios
    OLFC(; options = OLFCOptions(),
           generations = [0.],
           storages = [0.],
           converters = [0.]) =
           new(options, generations, storages, converters)
end

### Model
function build_model(mg::Microgrid, controller::OLFC, ω::Scenarios)
    # Sets
    nh, ns = controller.options.horizon, controller.options.nscenarios
    # Initialize
    m = Model(controller.options.solver.Optimizer)
    set_optimizer_attribute(m,"CPX_PARAM_SCRIND", 0)
    # Add investment variables
    add_investment_decisions!(m, mg.generations)
    add_investment_decisions!(m, mg.storages)
    add_investment_decisions!(m, mg.converters)
    # Fix their values
    fix_investment_decisions!(m, controller.generations, controller.storages, controller.converters)
    # Add operation decision variable with the non-anticipative structure
    add_operation_decisions!(m, mg.storages, nh, 1)
    add_operation_decisions!(m, mg.converters, nh, 1)
    # Add demand and generation variables to be fixed online, along with recourse variable (grid)
    add_operation_decisions!(m, mg.generations, nh, ns)
    add_operation_decisions!(m, mg.demands, nh, ns)
    add_operation_decisions!(m, mg.grids, nh, ns)
    # Add technical constraints
    add_technical_constraints!(m, mg.storages, mg.parameters.Δh, nh, 1)
    add_technical_constraints!(m, mg.converters, nh, 1)
    add_technical_constraints!(m, mg.grids, nh, ns)
    # Add power balance constraints
    add_power_balance!(m, mg, ω, Electricity, nh, ns, ispnet = true)
    add_power_balance!(m, mg, ω, Heat, nh, ns, ispnet = true)
    add_power_balance!(m, mg, ω, Hydrogen, nh, ns, ispnet = true)
    return m
end

### Offline
function initialize_controller!(mg::Microgrid, controller::OLFC, ω::Scenarios)
    # Preallocate
    preallocate!(mg, controller)
    # Build model
    println("Building the model...")
    controller.model = build_model(mg, controller, ω)
    # Scenario reduction
    println("Starting scenario reduction...")
    ω_reduced = reduce(controller.options.reducer, ω)[1]
    # Compute markov chain for scenario generation
    println("Initializing scenario generator...")
    controller.options.generator = initialize_generator!(controller.options.generator, [a for a in ω_reduced.generations]..., [a for a in ω_reduced.demands]...)
    # History
    controller.history = ω
    return controller
end

### Online
function compute_operation_decisions!(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::OLFC)
    # Forecast window
    window = h:min(mg.parameters.nh, h + controller.options.horizon - 1)
    # Compute forecast
    demands, generations, grids, probabilities = compute_forecast(h, y, s, mg, controller, window)
    # Fix model variables
    fix_operation_variables!(controller, demands, generations, [a.soc[h,y,s] * a.Erated[y,s] for a in mg.storages])
    # Add objective
    add_objective!(controller, mg, grids, probabilities, window)
    # Optimize
    optimize!(controller.model)
    # Assign controller values
    for k in 1:length(mg.storages)
        controller.decisions.storages[k][h,y,s] = value(controller.model[:p_dch][1,1,k] - controller.model[:p_ch][1,1,k])
    end
    for (k,a) in enumerate(mg.converters)
        if a isa Heater
            controller.decisions.converters[k][h,y,s] = - value(controller.model[:p_c][1,1,k])
        elseif a isa Electrolyzer
            controller.decisions.converters[k][h,y,s] = - value(controller.model[:p_c][1,1,k])
        elseif a isa FuelCell
            controller.decisions.converters[k][h,y,s] = value(controller.model[:p_c][1,1,k])
        end
    end
end
# OLFC forecast
function compute_forecast(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::OLFC, window::UnitRange{Int64})
    # To make the writting easier...
    nh, ns = controller.options.horizon, controller.options.nscenarios
    # Initial state
    s0 = vcat([a.carrier.power[h,y,s] / a.powerMax[y,s] for a in mg.generations], [a.carrier.power[h,y,s] for a in mg.demands])
    # Initial timestamp
    t0 = mg.demands[1].timestamp[h,y,s]
    # Demand and production forecast using the scenario generator
    forecast, probabilities = generate(controller.options.generator, s0, t0, nh, ny = ns)
    # Unpack the forecast vector
    generations, demands = [a.powerMax[y,s] * forecast[k] for (k,a) in enumerate(mg.generations)], forecast[length(mg.generations)+1:end]
    # Operating cost - unchanged along the horizon. We add zeros to guarantee a constant vector size
    grids = [(cost_in = vcat(a.cost_in[window,y,s], zeros(nh - length(window))), cost_out = vcat(a.cost_out[window,y,s], zeros(nh - length(window)))) for a in mg.grids]

    return demands, generations, grids, probabilities
end
# Fix OLFC variables
function fix_operation_variables!(controller::OLFC, demands::Vector{Array{Float64,3}}, generations::Vector{Array{Float64,3}}, soc_ini::Vector{Float64})
    # Fix demand values
    for (k,a) in enumerate(demands)
        fix.(controller.model[:p_d][:,:,k], a[:,:,1])
    end
    # Fix generation values
    for (k,a) in enumerate(generations)
        fix.(controller.model[:p_g][:,:,k], a[:,:,1])
    end
    # Fix initial SoC values
    for (k,a) in enumerate(soc_ini)
        set_normalized_coefficient.(controller.model[:soc_ini][:,k], controller.model[:r_sto][k], 0.)
        set_normalized_rhs.(controller.model[:soc_ini][:,k], a)
    end
end
# Add objective
function add_objective!(controller::OLFC, mg::Microgrid, tariff::Vector{NamedTuple{(:cost_in, :cost_out), Tuple{Vector{Float64},Vector{Float64}}}}, probabilities::Array{Float64,2}, window::UnitRange{Int64})
    # To make the writting easier...
    nh, ns = controller.options.horizon, controller.options.nscenarios
    # Initialize
    cost = AffExpr(0.)
    # Opex
    for (k,a) in enumerate(mg.grids)
        add_to_expression!(cost, sum(probabilities[s] * sum((controller.model[:p_in][h,s,k] * tariff[k].cost_in[h] - controller.model[:p_out][h,s,k] * tariff[k].cost_out[h]) * mg.parameters.Δh for h in 1:nh) for s in 1:ns))
    end
    # Add final cost
    for (k,a) in enumerate(mg.storages)
        if controller.options.seasonal_targets isa Nothing
            add_to_expression!(cost, - controller.options.final_cost[k] * controller.model[:soc][end,1,k])
        else
            if a isa H2Tank
                add_to_expression!(cost, - controller.options.final_cost[k] * (controller.model[:soc][end,1,k] - controller.options.seasonal_targets[window]))
            else
                add_to_expression!(cost, - controller.options.final_cost[k] * controller.model[:soc][end,1,k])
            end
        end
    end
    @objective(controller.model, Min, cost)
end


text:
using InfDecomp
using Read_PDF
using MathProgBase
# using Mosek
using Ipopt
using NLopt
using KNITRO
using ECOS
using SCS
#using PyCall
#unshift!(PyVector(pyimport("sys")["path"]), "")
#@pyimport cvxopt_solve
#-------------------------------
#READS P (CONTROL) PROBABILITY
#-------------------------------
function read_p(w::String)
	 pdf = Dict()
	 #pdf = read_pdf(w)
	 pdf = read_pdf("$w")
	 return pdf
	;
end

#-------------------------------------------------
#COMPUTES TOTAL VARIATION FROM TRUE DISTRIBUTION
#-------------------------------------------------
function tot_var{T1,T2,T3}(pdf::Dict{Tuple{T1,T2,T3},Float64}, true_pdf::Dict{Tuple{T1,T2,T3},Float64})
    a  = Float64[]
    b  = Float64[]
    u = 0
    w = Float64[]
    for (k,v) in pdf
        u = abs(v - get(true_pdf,k,0))
        push!(a, u)
    end
    for (k,v) in true_pdf
        if get(pdf,k,0) == 0
            push!(b, v)
        end
    end
    w = append!(a,b)
    return 0.5 * norm(w, 1)
end

#-------------------------------------------------
# Solutions & Statistics for Cvxopt
#-------------------------------------------------
# type Solution_and_Stats
#     var_num             :: Int64
#     x_sz                :: Int64
#     y_sz                :: Int64
#     z_sz                :: Int64
#     status              :: Symbol
#     obj_val             :: BigFloat
#     q_nonneg_viol       :: BigFloat
#     marginals_1         :: BigFloat
#     marginals_2         :: BigFloat
#     marginals_Inf       :: BigFloat
#     mu_nonneg_viol      :: BigFloat
#     complementarity_max :: BigFloat
#     complementarity_sum :: BigFloat
#     # MI_X_YZ           :: BigFloat
#     CI                  :: BigFloat
#     SI                  :: BigFloat
#     UI_Y                :: BigFloat
#     UI_Z                :: BigFloat
#     opt_time            :: String
#     # entropy_X           :: BigFloat
# end

Solution_and_Stats = InfDecomp.Solution_and_Stats

#-------
#RUN IT
#-------
function test(solver::Symbol,pdf_fn::String, w2::String, tmpFloatDatatype::DataType ; stats_file = "")
    p = read_p("../PDFs/Data/" * pdf_fn)
    true_p = read_p("../PDFs/Data/" * w2)
    time_l = 1.0
    if (solver != :Cvxopt)
        if (solver == :Mosek)
            InfDecomp.set_copy_sol_behaviour(true)
            println("Start optimization with Mosek...")
	    sd,myeval,model = InfDecomp.do_it(p, Mosek.solver(MSK_IPAR_INTPNT_MULTI_THREAD=0, MSK_IPAR_INTPNT_MAX_ITERATIONS=500, MSK_DPAR_OPTIMIZER_MAX_TIME=time_l),tmpFloatDatatype)
        elseif (solver == :Ipopt)
            InfDecomp.set_copy_sol_behaviour(false)
            println("Start optimization with Ipopt...")
	    sd,myeval,model = InfDecomp.do_it(p,Ipopt.IpoptSolver(max_cpu_time=time_l),tmpFloatDatatype)
            # elseif (solver == :NLopt_LD_SLSQP)
            #     println("Start optimization with NLopt/LD_SLSQP...")
            #     sd,myeval,model = InfDecomp.do_it(p,NLopt.NLoptSolver(algorithm=:LD_SLSQP))
        elseif (solver == :Knitro_Ip)
            InfDecomp.set_copy_sol_behaviour(false)
            println("Start optimization with Knitro/Interior_pt...")
	    sd,myeval,model = InfDecomp.do_it(p,KNITRO.KnitroSolver(algorithm=1,maxtime_cpu=time_l),tmpFloatDatatype)
        elseif (solver == :Knitro_IpCG)
            InfDecomp.set_copy_sol_behaviour(false)
            println("Start optimization with Knitro/Interior_pt_CG...")
	    sd,myeval,model = InfDecomp.do_it(p,KNITRO.KnitroSolver(algorithm=2,maxtime_cpu=time_l),tmpFloatDatatype)
        elseif (solver == :Knitro_AS)
            InfDecomp.set_copy_sol_behaviour(false)
            println("Start optimization with Knitro/Active_Set...")
	    sd,myeval,model = InfDecomp.do_it(p,KNITRO.KnitroSolver(algorithm=3,maxtime_cpu=time_l),tmpFloatDatatype)
        elseif (solver == :Knitro_SQP)
            InfDecomp.set_copy_sol_behaviour(false)
            println("Start optimization with Knitro/SQP...")
	    sd,myeval,model = InfDecomp.do_it(p,KNITRO.KnitroSolver(algorithm=4,maxtime_cpu=time_l),tmpFloatDatatype)

        elseif solver == :ECOS_L
            println("Starting conic optimization with ECOS (large).")
	    sd,myeval,model = InfDecomp.do_it(p,ECOS.ECOSSolver(),tmpFloatDatatype, model_type=solver)
        elseif solver == :ECOS_S
            println("Starting conic optimization with ECOS (small).")
	    sd,myeval,model = InfDecomp.do_it(p,ECOS.ECOSSolver(),tmpFloatDatatype, model_type=solver)
        elseif solver == :SCS_L
            println("Starting conic optimization with SCS (large).")
	    sd,myeval,model = InfDecomp.do_it(p,SCS.SCSSolver(),tmpFloatDatatype, model_type=solver)
        elseif solver == :SCS_S
            println("Starting conic optimization with SCS (small).")
	    sd,myeval,model = InfDecomp.do_it(p,SCS.SCSSolver(),tmpFloatDatatype, model_type=solver)
        elseif solver == :SCS_D
            println("Starting conic optimization with SCS (Dual).")
	    sd,myeval,model = InfDecomp.do_it(p,SCS.SCSSolver(),tmpFloatDatatype, model_type=solver)
        elseif solver == :My_GradDesc
            println("Starting My Gradient Descent.")
	    sd,myeval,model = InfDecomp.do_it(p,nothing,tmpFloatDatatype, model_type=solver)
        end

        feasstats = InfDecomp.check_feasibility(pdf_fn, model,myeval,solver)
        if stats_file != ""
            open(stats_file, "a") do ffile
                for i in fieldnames(feasstats)
                    print(ffile,"\t",getfield(feasstats,i),", ")
                end
                println(ffile,"")
            end
        end

        if solver ∈ [:ECOS_L,:SCS_L,:ECOS_S,:SCS_S,:SCS_D]
            # Conic Program:
            # do nothing
        else
            # IPM
            # do nothing
        end #^ if solver not ECOS/SCS

        return sd,myeval,model
    else # solver==:CVXOPT
        println("Start optimization with Cvxopt...")
        pdf = PyDict(p)
        fstat  = Solution_and_Stats( 0,0,0,0," ",  0,0,0,0,0,0,0,0,0,0,0,0," ")
        fstat.var_num, fstat.x_sz, fstat.y_sz, fstat.z_sz, fstat.status, fstat.obj_val, fstat.q_nonneg_viol, fstat.marginals_1, fstat.marginals_2, fstat.marginals_Inf, fstat.mu_nonneg_viol, fstat.complementarity_max, fstat.complementarity_sum, fstat.CI, fstat.SI, fstat.UI_Y, fstat.UI_Z = cvxopt_solve.solve_PDF(p,time_l*1000)
        fstat.opt_time = -1

        s = tot_var(p,true_p)
        open("feas_stats_64_gu_9.csv", "a") do ffile
            print(ffile,"# filename, Solver")
            for i in fieldnames(fstat)
                print(ffile,",\t",i)
            end
            print(ffile,", Tot_var\n",pdf_fn,",Cvxopt")
            for i in fieldnames(fstat)
                print(ffile,",\t",getfield(fstat,i))
            end
            print(ffile,", $s\n")
        end
        return 0
    end #^ if/else solver!=:CVXOPT
end

; # EOF

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
## Elementary Rand Variables
## =========================

abstract type ElementaryRandVar{T} <: SymbolicRandVar{T} end
dims(X::ElementaryRandVar) = union([Set(X.dim), map(dims, args(X))...]...)::Set{Int}
has_single_dim(::ElementaryRandVar) = true
num_params{T <: ElementaryRandVar}(X::Type{T}) = length(fieldnames(X)) - 1

abstract type ClosedFormQuantileRandVar{T} <: ElementaryRandVar{T} end

## Continuous RandVars
## ===================

"Arcsine"
immutable ArcsineRandVar{T <: Real, A <: SymbolicRandVar, B <: SymbolicRandVar} <: ElementaryRandVar{T}
  dim::Id
  a::A
  b::B
  ArcSineRandVar{T1<:Real}(dim::Id, a::SymbolicRandVar{T1}, b::SymbolicRandVar{T1}) = new{T1, T1, T1}(dim, a, b)
end

"Uniformly distributed RandVar"
immutable UniformRandVar{T <: Real, A <: Real, B <: Real} <: ClosedFormQuantileRandVar{T}
  dim::Id
  lb::SymbolicRandVar{A}
  ub::SymbolicRandVar{B}
end

quantile_expr(x::UniformRandVar) = (x.lb - x.ub) * omega_component(x.dim) + x.lb

args(X::UniformRandVar) = tuple(X.lb, X.ub)

# "Normally distributed RandVar"
# immutable NormalRandVar{T <: Real, A <: Real, B <: Real} <: ElementaryRandVar{T}
#   dim::Id
#   μ::SymbolicRandVar{A}
#   σ::SymbolicRandVar{B}
# end
#
# # param_types(X::Type{NormalRandVar}) = [Float64, Float64]
#
# args(X::NormalRandVar) = tuple(X.μ, X.σ)

# "Beta distributed RandVar"
# immutable BetaRandVar{T <: Real, A <: SymbolicRandVar, B <: SymbolicRandVar} <: ElementaryRandVar{T}
#   dim::Id
#   α::A
#   β::B
# end

## Discrete Distritbuions
## ======================
"Bernoulli Distribution"
immutable BernoulliRandVar{Bool, P <: Real} <: ClosedFormQuantileRandVar{Bool}
  dim::Id
  p::SymbolicRandVar{P}
end

args(X::BernoulliRandVar) = tuple(X.p)
quantile_expr(X::BernoulliRandVar) = X.p >= omega_component(X.dim)

# # FIXME: Closed form?
# "Binomial RandVar"
# immutable BinomialRandVar{T <: Integer, N <: Integer, P <: Real} <: ElementaryRandVar{T}
#   dim::Id
#   n::SymbolicRandVar{N}
#   p::SymbolicRandVar{P}
# end
#
# args(X::BinomialRandVar) = tuple(X.n, X.p)

# "Categorical Distribution"
# immutable CategoricalRandVar{T <: Integer, p::}
# end

"Discrete Uniform RandVar"
immutable DiscreteUniformRandVar{T <: Integer, A <: Integer, B <: Integer} <: ElementaryRandVar{T}
  dim::Id
  a::SymbolicRandVar{A}
  b::SymbolicRandVar{B}
end

args(X::DiscreteUniformRandVar) = tuple(X.a, X.b)

# "Uniformly distributed RandVar"
# immutable PoissonRandVar{T <: Real, A <: Real}
#   dim::Id
#   λ::SymbolicRandVar{A}
# end

# args(X::PoissonRandVar) = tuple(X.λ)