text:
# byproducts.jl

"""
    StartWorkers(nwrkrs::Int)

Start workers if needed.
"""
function StartWorkers(nwrkrs::Int)
   set_workers = nwrkrs
   nworkers() < set_workers ? addprocs(set_workers) : nothing
   nworkers()
end

"""
    TaskDriver(indx,fn)

Broacast / distribute task (fn; e.g. loop_task1) over indices (indx; e.g. file indices)

Examples:

```
using CbiomesProcessing, Distributed, SparseArrays
TaskDriver(1,CbiomesProcessing.loop_task1)

StartWorkers(4)
@everywhere using CbiomesProcessing, SparseArrays
TaskDriver(1:4,CbiomesProcessing.loop_task1)
```

Visualize results:

```
using FortranFiles, Plots
k=1
recl=720*360*4
fil="diags_interp/ETAN/ETAN.0000000732.data"
f =  FortranFile(fil,"r",access="direct",recl=recl,convert="big-endian")
tmp=read(f,rec=k,(Float32,(720,360))); close(f)
heatmap(tmp)
```
"""
function TaskDriver(indx::Union{UnitRange{Int},Array{Int,1},Int},fn::Function)
    i=collect(indx)
    length(i)>1 ? i=distribute(i) : nothing
    isa(i,DArray) ? println(i.indices) : nothing
    fn.(i)
end


text:
#!/usr/bin/env julia

#A few utilities to estimate the percentage of detached
# dendritic spines

"""

    find_trunk_mapping( om, trunk_list )

Given an overlap matrix and list of segment ids, find the segments
which overlap the most in the compared segmentation. If the segment
arguments are trunks, the results are likely to correspond to the
same trunks in the comparison segmentation.
"""
function find_trunk_mapping( om, trunk_list )

    trunk_om = om[trunk_list,:];

    maxs, inds = findmax(trunk_om,2)
    matches = ind2sub( size(trunk_om), inds[:] )[2]

    Dict( trunk_list[i] => matches[i] for i in eachindex(matches) )
end

"""

    synapse_recovery( uncorr_es, corr_es, trunk_mapping )

Given synapse tables (as read by load_edges) and a trunk mapping,
determine the ratio between the number of synapses innervating the
mapped segments in each segmentation.
"""
function synapse_recovery( uncorr_es, corr_es, trunk_mapping )

    corr_segs = corr_es[:,2]; uncorr_segs = uncorr_es[:,2];

    corr_trunks = Set(keys(trunk_mapping))
    uncorr_trunks = Set(values(trunk_mapping))

    num_corr_trunk_syns = count( x -> x[2] in corr_trunks, corr_segs )
    num_uncorr_trunk_syns = count( x -> x[2] in uncorr_trunks, uncorr_segs )

    num_uncorr_trunk_syns / num_corr_trunk_syns
end

ext: jl


text:
# Introspection for kernel benchmarks to see native code.
function introspect(A = nothing)
    ### Kernel
    options = [
        "Read Only"    => ReadOnly,
        "Write Only"   => WriteOnly,
        "Read + Write" => ReadWrite,
    ]

    menu = RadioMenu(first.(options); pagesize = 10)
    choice = request("Choose a kernel to inspect: ", menu)
    checkexit(choice) && return nothing

    # Pull out the function
    kernel = last(options[choice])

    ### Iterator
    options = [
        "Sequential" => Sequential,
        "Pseudo Random" => PseudoRandom,
    ]

    menu = RadioMenu(first.(options); pagesize = 10)
    choice = request("Choose an Iterator: ", menu)
    checkexit(choice) && return nothing

    # Pull out the function
    iterator = last(options[choice])

    ### Element Type
    # Now, ask which element type they's like to use and create
    if isnothing(A)
        options = [
            "Float32" => Float32,
            "Float64" => Float64,
        ]

        menu = RadioMenu(first.(options); pagesize = 10)
        choice = request("Choose a Element Type: ", menu)
        checkexit(choice) && return nothing
        T = last(options[choice])
        A = Vector{T}(undef, 1024)
    end

    ### Choose a vector size
    options = [
        "4" => 4,
        "8" => 8,
    ]

    # If we have a 4-byte datatype, we can have a vector width of 16 as well.
    if sizeof(eltype(A)) < 8
        push!(options, "16" => 16)
    end

    menu = RadioMenu(first.(options); pagesize = 10)
    choice = request("Choose a vector size ", menu)
    checkexit(choice) && return nothing
    vectorsize = last(options[choice])

    ### Create an array view?
    options = [
        "Yes" => true,
        "No" => false,
    ]

    menu = RadioMenu(first.(options); pagesize = 10)
    choice = request("Create Array View? ", menu)
    checkexit(choice) && return nothing
    makeview = last(options[choice])
    if makeview
        # Create a SubArray of the array.
        # What's really important is that we get the correct `SubArray` type.
        A = view(A, 1:div(length(A),2))
    end

    ### Now, we ask some options for non-temporal and aligned
    options = [
        "Yes" => Nontemporal,
        "No" => Standard,
    ]

    menu = RadioMenu(first.(options); pagesize = 10)
    choice = request("Non-temporal Instructions? ", menu)
    checkexit(choice) && return nothing
    nontemporal = last(options[choice])

    # If this function supports unrolling, make this an option.
    options = [
        "1" => 1,
        "2" => 2,
        "4" => 4,
        "8" => 8,
    ]

    menu = RadioMenu(first.(options); pagesize = 10)
    choice = request("Unroll Kernel? ", menu)
    checkexit(choice) && return nothing
    unroll = last(options[choice])

    # Generate the "KernelParam" type
    params = KernelParam(
        kernel = kernel,
        loadtype = nontemporal,
        storetype = nontemporal,
        iterator = iterator,
        eltype = eltype(A),
        vectorsize = vectorsize,
        unroll = unroll
    )

    # Now that we have everything ready to go, we create the appropriate calls to the
    # innermost kernels.
    generate_code_native_call(A, params)
end

checkexit(x) = (x == -1)

function generate_code_native_call(A, K::KernelParam)
    # No reason to specialize on these arguments
    @nospecialize

    V = reinterpret(K, A)
    itr = makeitr(V, K)

    return native_call(_execute!, (V, itr, K))
end

function native_call(f, args)
    @nospecialize f args
    types = typeof.(args)
    @show types
    call = """
        code_native(KernelBenchmarks.$f, $types; syntax=:intel, debuginfo=:none)
    """
    code_native(f, types; syntax=:intel, debuginfo=:none)
    printstyled(stdout, "\n\nRun this to recreate\n"; color = :green)
    println(call)
    return nothing
end


text:
#==============================================================================#
# LexRuntime.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/runtime.lex-2016-11-28.normal.json
#==============================================================================#

__precompile__()

module LexRuntime

using AWSCore

"""
    using AWSSDK.LexRuntime.post_content
    post_content([::AWSConfig], arguments::Dict)
    post_content([::AWSConfig]; botName=, botAlias=, userId=, *header:* Content-Type=, inputStream=, <keyword arguments>)

    using AWSCore.Services.runtime_lex
    runtime_lex([::AWSConfig], "POST", "/bot/{botName}/alias/{botAlias}/user/{userId}/content", arguments::Dict)
    runtime_lex([::AWSConfig], "POST", "/bot/{botName}/alias/{botAlias}/user/{userId}/content", botName=, botAlias=, userId=, *header:* Content-Type=, inputStream=, <keyword arguments>)

# PostContent Operation

Sends user input (text or speech) to Amazon Lex. Clients use this API to send text and audio requests to Amazon Lex at runtime. Amazon Lex interprets the user input using the machine learning model that it built for the bot.

The `PostContent` operation supports audio input at 8kHz and 16kHz. You can use 8kHz audio to achieve higher speech recognition accuracy in telephone audio applications.

In response, Amazon Lex returns the next message to convey to the user. Consider the following example messages:

*   For a user input "I would like a pizza," Amazon Lex might return a response with a message eliciting slot data (for example, `PizzaSize`): "What size pizza would you like?".

*   After the user provides all of the pizza order information, Amazon Lex might return a response with a message to get user confirmation: "Order the pizza?".

*   After the user replies "Yes" to the confirmation prompt, Amazon Lex might return a conclusion statement: "Thank you, your cheese pizza has been ordered.".

Not all Amazon Lex messages require a response from the user. For example, conclusion statements do not require a response. Some messages require only a yes or no response. In addition to the `message`, Amazon Lex provides additional context about the message in the response that you can use to enhance client behavior, such as displaying the appropriate client user interface. Consider the following examples:

*   If the message is to elicit slot data, Amazon Lex returns the following context information:

    *   `x-amz-lex-dialog-state` header set to `ElicitSlot`

    *   `x-amz-lex-intent-name` header set to the intent name in the current context

    *   `x-amz-lex-slot-to-elicit` header set to the slot name for which the `message` is eliciting information

    *   `x-amz-lex-slots` header set to a map of slots configured for the intent with their current values

*   If the message is a confirmation prompt, the `x-amz-lex-dialog-state` header is set to `Confirmation` and the `x-amz-lex-slot-to-elicit` header is omitted.

*   If the message is a clarification prompt configured for the intent, indicating that the user intent is not understood, the `x-amz-dialog-state` header is set to `ElicitIntent` and the `x-amz-slot-to-elicit` header is omitted.

In addition, Amazon Lex also returns your application-specific `sessionAttributes`. For more information, see [Managing Conversation Context](http://docs.aws.amazon.com/lex/latest/dg/context-mgmt.html).

# Arguments

## `botName = ::String` -- *Required*
Name of the Amazon Lex bot.

## `botAlias = ::String` -- *Required*
Alias of the Amazon Lex bot.

## `userId = ::String` -- *Required*
The ID of the client application user. Amazon Lex uses this to identify a user's conversation with your bot. At runtime, each request must contain the `userID` field.

To decide the user ID to use for your application, consider the following factors.

*   The `userID` field must not contain any personally identifiable information of the user, for example, name, personal identification numbers, or other end user personal information.

*   If you want a user to start a conversation on one device and continue on another device, use a user-specific identifier.

*   If you want the same user to be able to have two independent conversations on two different devices, choose a device-specific identifier.

*   A user can't have two independent conversations with two different versions of the same bot. For example, a user can't have a conversation with the PROD and BETA versions of the same bot. If you anticipate that a user will need to have conversation with two different versions, for example, while testing, include the bot alias in the user ID to separate the two conversations.

## `*header:* x-amz-lex-session-attributes = ::String`
You pass this value as the `x-amz-lex-session-attributes` HTTP header.

Application-specific information passed between Amazon Lex and a client application. The value must be a JSON serialized and base64 encoded map with string keys and values. The total size of the `sessionAttributes` and `requestAttributes` headers is limited to 12 KB.

For more information, see [Setting Session Attributes](http://docs.aws.amazon.com/lex/latest/dg/context-mgmt.html#context-mgmt-session-attribs).

## `*header:* x-amz-lex-request-attributes = ::String`
You pass this value as the `x-amz-lex-request-attributes` HTTP header.

Request-specific information passed between Amazon Lex and a client application. The value must be a JSON serialized and base64 encoded map with string keys and values. The total size of the `requestAttributes` and `sessionAttributes` headers is limited to 12 KB.

The namespace `x-amz-lex:` is reserved for special attributes. Don't create any request attributes with the prefix `x-amz-lex:`.

For more information, see [Setting Request Attributes](http://docs.aws.amazon.com/lex/latest/dg/context-mgmt.html#context-mgmt-request-attribs).

## `*header:* Content-Type = ::String` -- *Required*
You pass this value as the `Content-Type` HTTP header.

Indicates the audio format or text. The header value must start with one of the following prefixes:

*   PCM format, audio data must be in little-endian byte order.

    *   audio/l16; rate=16000; channels=1

    *   audio/x-l16; sample-rate=16000; channel-count=1

    *   audio/lpcm; sample-rate=8000; sample-size-bits=16; channel-count=1; is-big-endian=false

*   Opus format

    *   audio/x-cbr-opus-with-preamble; preamble-size=0; bit-rate=256000; frame-size-milliseconds=4

*   Text format

    *   text/plain; charset=utf-8

## `*header:* Accept = ::String`
You pass this value as the `Accept` HTTP header.

The message Amazon Lex returns in the response can be either text or speech based on the `Accept` HTTP header value in the request.

*   If the value is `text/plain; charset=utf-8`, Amazon Lex returns text in the response.

*   If the value begins with `audio/`, Amazon Lex returns speech in the response. Amazon Lex uses Amazon Polly to generate the speech (using the configuration you specified in the `Accept` header). For example, if you specify `audio/mpeg` as the value, Amazon Lex returns speech in the MPEG format.

    The following are the accepted values:

    *   audio/mpeg

    *   audio/ogg

    *   audio/pcm

    *   text/plain; charset=utf-8

    *   audio/* (defaults to mpeg)

## `inputStream = blob` -- *Required*
User input in PCM or Opus audio format or text format as described in the `Content-Type` HTTP header.

You can stream audio data to Amazon Lex or you can create a local buffer that captures all of the audio data before sending. In general, you get better performance if you stream audio data rather than buffering the data locally.

# Returns

`PostContentResponse`

# Exceptions

`NotFoundException`, `BadRequestException`, `LimitExceededException`, `InternalFailureException`, `ConflictException`, `UnsupportedMediaTypeException`, `NotAcceptableException`, `RequestTimeoutException`, `DependencyFailedException`, `BadGatewayException` or `LoopDetectedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/runtime.lex-2016-11-28/PostContent)
"""
@inline post_content(aws::AWSConfig=default_aws_config(); args...) = post_content(aws, args)

@inline post_content(aws::AWSConfig, args) = AWSCore.Services.runtime_lex(aws, "POST", "/bot/{botName}/alias/{botAlias}/user/{userId}/content", args)

@inline post_content(args) = post_content(default_aws_config(), args)

"""
    using AWSSDK.LexRuntime.post_text
    post_text([::AWSConfig], arguments::Dict)
    post_text([::AWSConfig]; botName=, botAlias=, userId=, inputText=, <keyword arguments>)

    using AWSCore.Services.runtime_lex
    runtime_lex([::AWSConfig], "POST", "/bot/{botName}/alias/{botAlias}/user/{userId}/text", arguments::Dict)
    runtime_lex([::AWSConfig], "POST", "/bot/{botName}/alias/{botAlias}/user/{userId}/text", botName=, botAlias=, userId=, inputText=, <keyword arguments>)

# PostText Operation

Sends user input (text-only) to Amazon Lex. Client applications can use this API to send requests to Amazon Lex at runtime. Amazon Lex then interprets the user input using the machine learning model it built for the bot.

In response, Amazon Lex returns the next `message` to convey to the user an optional `responseCard` to display. Consider the following example messages:

*   For a user input "I would like a pizza", Amazon Lex might return a response with a message eliciting slot data (for example, PizzaSize): "What size pizza would you like?"

*   After the user provides all of the pizza order information, Amazon Lex might return a response with a message to obtain user confirmation "Proceed with the pizza order?".

*   After the user replies to a confirmation prompt with a "yes", Amazon Lex might return a conclusion statement: "Thank you, your cheese pizza has been ordered.".

Not all Amazon Lex messages require a user response. For example, a conclusion statement does not require a response. Some messages require only a "yes" or "no" user response. In addition to the `message`, Amazon Lex provides additional context about the message in the response that you might use to enhance client behavior, for example, to display the appropriate client user interface. These are the `slotToElicit`, `dialogState`, `intentName`, and `slots` fields in the response. Consider the following examples:

*   If the message is to elicit slot data, Amazon Lex returns the following context information:

    *   `dialogState` set to ElicitSlot

    *   `intentName` set to the intent name in the current context

    *   `slotToElicit` set to the slot name for which the `message` is eliciting information

    *   `slots` set to a map of slots, configured for the intent, with currently known values

*   If the message is a confirmation prompt, the `dialogState` is set to ConfirmIntent and `SlotToElicit` is set to null.

*   If the message is a clarification prompt (configured for the intent) that indicates that user intent is not understood, the `dialogState` is set to ElicitIntent and `slotToElicit` is set to null.

In addition, Amazon Lex also returns your application-specific `sessionAttributes`. For more information, see [Managing Conversation Context](http://docs.aws.amazon.com/lex/latest/dg/context-mgmt.html).

# Arguments

## `botName = ::String` -- *Required*
The name of the Amazon Lex bot.

## `botAlias = ::String` -- *Required*
The alias of the Amazon Lex bot.

## `userId = ::String` -- *Required*
The ID of the client application user. Amazon Lex uses this to identify a user's conversation with your bot. At runtime, each request must contain the `userID` field.

To decide the user ID to use for your application, consider the following factors.

*   The `userID` field must not contain any personally identifiable information of the user, for example, name, personal identification numbers, or other end user personal information.

*   If you want a user to start a conversation on one device and continue on another device, use a user-specific identifier.

*   If you want the same user to be able to have two independent conversations on two different devices, choose a device-specific identifier.

*   A user can't have two independent conversations with two different versions of the same bot. For example, a user can't have a conversation with the PROD and BETA versions of the same bot. If you anticipate that a user will need to have conversation with two different versions, for example, while testing, include the bot alias in the user ID to separate the two conversations.

## `sessionAttributes = ::Dict{String,String}`
Application-specific information passed between Amazon Lex and a client application.

For more information, see [Setting Session Attributes](http://docs.aws.amazon.com/lex/latest/dg/context-mgmt.html#context-mgmt-session-attribs).

## `requestAttributes = ::Dict{String,String}`
Request-specific information passed between Amazon Lex and a client application.

The namespace `x-amz-lex:` is reserved for special attributes. Don't create any request attributes with the prefix `x-amz-lex:`.

For more information, see [Setting Request Attributes](http://docs.aws.amazon.com/lex/latest/dg/context-mgmt.html#context-mgmt-request-attribs).

## `inputText = ::String` -- *Required*
The text that the user entered (Amazon Lex interprets this text).

# Returns

`PostTextResponse`

# Exceptions

`NotFoundException`, `BadRequestException`, `LimitExceededException`, `InternalFailureException`, `ConflictException`, `DependencyFailedException`, `BadGatewayException` or `LoopDetectedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/runtime.lex-2016-11-28/PostText)
"""
@inline post_text(aws::AWSConfig=default_aws_config(); args...) = post_text(aws, args)

@inline post_text(aws::AWSConfig, args) = AWSCore.Services.runtime_lex(aws, "POST", "/bot/{botName}/alias/{botAlias}/user/{userId}/text", args)

@inline post_text(args) = post_text(default_aws_config(), args)

end # module LexRuntime

#==============================================================================#
# End of file
#==============================================================================#


text:
# min X[1,1] + X[2,2]    max y
#     X[2,1] = 1         [0   y/2     [ 1  0
#                         y/2 0    <=   0  1]
#     X >= 0              y free
# Optimal solution:
#
#     ⎛ 1   1 ⎞
# X = ⎜       ⎟           y = 2
#     ⎝ 1   1 ⎠
@testset "Differentiating simple PSD back" begin

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    X = MOI.add_variables(model, 3)
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{Float64}(vov),
        MOI.PositiveSemidefiniteConeTriangle(2)
    )

    c  = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, X[2]))],
            [-1.0]
        ),
        MOI.Zeros(1)
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    x = MOI.get(model, MOI.VariablePrimal(), X)

    cone_types = unique([S for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())])
    conic_form = DiffOpt.GeometricConicForm{Float64}()
    cones = conic_form.constraints.sets
    DiffOpt.set_set_types(cones, cone_types)
    index_map = MOI.copy_to(conic_form, model)

    @test x ≈ ones(3) atol=ATOL rtol=RTOL

    MOI.set(model, DiffOpt.BackwardInVariablePrimal(), X[1], 1.0)

    DiffOpt.backward(model)

    db = MOI.constant(MOI.get(model, DiffOpt.BackwardOutConstraint(), c))

    @test db ≈ [-1.0]  atol=ATOL rtol=RTOL
end


text:
struct Edge
    id::Union{Integer, Nothing}
    relation::String
    src_node::Node
    dest_node::Node
    properties::Dict
end

Edge(relation::String, src_node::Union{Node, Integer}, dest_node::Union{Node, Integer}, properties::Dict) = Edge(nothing, relation, src_node, dest_node, properties)
Edge(relation::String, src_node::Union{Node, Integer}, dest_node::Union{Node, Integer}) = Edge(nothing, relation, src_node, dest_node, Dict())

function _get_node_str(node::Node)
    return string(node)
end

function _get_node_str(node::Integer)
    return "()"
end

function string(edge::Edge)
    src_node_str = _get_node_str(edge.src_node)
    dest_node_str = _get_node_str(edge.dest_node)

    relation = ""
    props_str = ""

    if edge.relation != ""
        relation = ":" * edge.relation
    end

    if length(edge.properties) != 0
        props = ["$prop_name: " * prop_value_to_string(prop_value) for (prop_name, prop_value) in pairs(edge.properties)]
        props_str = "{" * join(props, ",") * "}"
    end
    "$src_node_str-[$relation $props_str]->$dest_node_str"
end

function isequal(x::Edge, y::Edge)
    if x.id != nothing && y.id != nothing && x.id == y.id
        return true
    else
        if x.relation == y.relation && x.src_node == y.src_node && x.dest_node == y.dest_node && x.properties == y.properties
            return true
        end
    end
    return false
end

Base.:(==)(x::Edge, y::Edge) = isequal(x, y)


text:
export opposite_algebra

mutable struct AbsAlgAssMor{R, S, T} <: Map{R, S, HeckeMap, AbsAlgAssMor}
  header::MapHeader{R, S}

  mat::T
  imat::T
  c_t::T
  d_t::T

  function AbsAlgAssMor{R, S, T}(A::R, B::S) where {R, S, T}
    z = new{R, S, T}()
    z.header = MapHeader(A, B)
    return z
  end

  function AbsAlgAssMor{R, S, T}(A::R, B::S, M::T) where {R, S, T}
    z = new{R, S, T}()
    z.c_t = similar(M, 1, dim(A))
    z.d_t = similar(M, 1, dim(B))
    z.mat = M

    function image(a)
      for i in 1:dim(A)
        z.c_t[1, i] = coeffs(a, copy = false)[i]
      end
      s = Vector{elem_type(base_ring(B))}(undef, dim(B))
      #mul!(z.d_t, z.c_t, M) # there is no mul! for Generic.Mat
      z.d_t = z.c_t*M
      for i in 1:dim(B)
        s[i] = z.d_t[1, i]
      end

      return B(s)
    end

    z.header = MapHeader(A, B, image)
    return z
  end

  function AbsAlgAssMor{R, S, T}(A::R, B::S, M::T, N::T) where {R, S, T}
    z = new{R, S, T}()
    z.c_t = similar(M, 1, dim(A))
    z.d_t = similar(M, 1, dim(B))

    z.mat = M
    z.imat = N

    function image(a)
      for i in 1:dim(A)
        z.c_t[1, i] = coeffs(a, copy = false)[i]
      end
      s = Vector{elem_type(base_ring(B))}(undef, dim(B))
      #mul!(z.d_t, z.c_t, M) # there is no mul! for Generic.Mat
      z.d_t = z.c_t * M
      for i in 1:dim(B)
        s[i] = z.d_t[1, i]
      end

      return B(s)
    end

    function preimage(a)
      for i in 1:dim(B)
        z.d_t[1, i] = coeffs(a, copy = false)[i]
      end
      s = Vector{elem_type(base_ring(A))}(undef, dim(A))
      z.c_t = z.d_t * N
      for i in 1:dim(A)
        s[i] = z.c_t[1, i]
      end
      return A(s)
    end

    z.header = MapHeader(A, B, image, preimage)
    return z
  end
end

inv(a::AbsAlgAssMor) = pseudo_inv(a)

#mutable struct AlgAssMor{R, S, T} <: Map{AlgAss{R}, AlgAss{S}, HeckeMap, AlgAssMor}
#  header::MapHeader{AlgAss{R}, AlgAss{S}}
#
#  mat::T
#  imat::T
#  c_t::T
#  d_t::T
#
#  function AlgAssMor(A::AlgAss{R}, B::AlgAss{S}, M::T) where {R, S, T}
#    z = new{R, S, T}()
#    z.c_t = similar(M, 1, dim(A))
#    z.d_t = similar(M, 1, dim(B))
#    z.mat = M
#
#    function image(a::AlgAssElem)
#      for i in 1:dim(A)
#        z.c_t[1, i] = a.coeffs[i]
#      end
#      s = Vector{S}(undef, dim(B))
#      #mul!(z.d_t, z.c_t, M) # there is no mul! for Generic.Mat
#      z.d_t = z.c_t*M
#      for i in 1:dim(B)
#        s[i] = z.d_t[1, i]
#      end
#
#      return AlgAssElem{S}(B, s)
#    end
#
#    z.header = MapHeader(A, B, image)
#    return z
#  end
#
#  function AlgAssMor(A::AlgAss{R}, B::AlgAss{S}, M::T, N::T) where {R, S, T}
#    z = new{R, S, T}()
#    z.c_t = similar(M, 1, dim(A))
#    z.d_t = similar(M, 1, dim(B))
#
#    z.mat = M
#    z.imat = N
#
#    function image(a::AlgAssElem)
#      for i in 1:dim(A)
#        z.c_t[1, i] = a.coeffs[i]
#      end
#      s = Vector{S}(undef, dim(B))
#      #mul!(z.d_t, z.c_t, M) # there is no mul! for Generic.Mat
#      z.d_t = z.c_t * M
#      for i in 1:dim(B)
#        s[i] = z.d_t[1, i]
#      end
#
#      return AlgAssElem{S}(B, s)
#    end
#
#    function preimage(a::AlgAssElem)
#      for i in 1:dim(B)
#        z.d_t[1, i] = a.coeffs[i]
#      end
#      s = Vector{R}(undef, dim(A))
#      z.c_t = z.d_t * N
#      for i in 1:dim(A)
#        s[i] = z.c_t[1, i]
#      end
#      return AlgAssElem{R}(A, s)
#    end
#
#    z.header = MapHeader(A, B, image, preimage)
#    return z
#  end
#end

function compose_and_squash(f::AbsAlgAssMor{R, U, T}, g::AbsAlgAssMor{S, R, T}) where {R, T, S, U}
  if isdefined(f, :imat) && isdefined(g, :imat)
    return hom(domain(g), codomain(f), g.mat * f.mat, f.imat * g.imat)
  else
    return hom(domain(g), codomain(f), g.mat * f.mat)
  end
end

function hom(A::R, B::S, M::T) where {R <: AbsAlgAss, S <: AbsAlgAss, T}
  return AbsAlgAssMor{R, S, T}(A, B, M)
end

function hom(A::R, B::S, M::T, N::T) where {R <: AbsAlgAss, S <: AbsAlgAss, T}
  return AbsAlgAssMor{R, S, T}(A, B, M, N)
end

#function hom(A::AlgAss{R}, B::AlgAss{S}, M::T) where {R <: AlgAss, S <: AlgAss, T}
#  return AlgAssMor{R, S, T}(A, B, M)
#end
#
#function hom(A::AlgAss{R}, B::AlgAss{S}, M::T, N::T) where {R <: AlgAss, S <: AlgAss, T}
#  return AlgAssMor{R, S, T}(A, B, M, N)
#end

function haspreimage(m::AbsAlgAssMor, a::AbsAlgAssElem)
  if isdefined(m, :imat)
    return true, preimage(m, a)
  end

  A = parent(a)
  t = matrix(base_ring(A), 1, dim(A), coeffs(a))
  b, p = can_solve_with_solution(m.mat, t, side = :left)
  if b
    return true, domain(m)([ p[1, i] for i = 1:nrows(m.mat) ])
  else
    return false, zero(domain(m))
  end
end

################################################################################
#
#  Morphisms between algebras and number fields
#
################################################################################

# S is the type of the algebra, T the element type of the algebra
mutable struct AbsAlgAssToNfAbsMor{S, T, U, V} <: Map{S, U, HeckeMap, AbsAlgAssToNfAbsMor}
  header::MapHeader{S, U}
  mat::V
  imat::V
  t::V # dummy vector used in image and preimage
  tt::V # another dummy vector

  function AbsAlgAssToNfAbsMor{S, T, U, V}(A::S, K::U, M::V, N::V) where {S, T, U, V}

    z = new{S, T, U, V}()
    z.mat = M
    z.imat = N
    z.t = zero_matrix(base_field(K), 1, dim(A))
    z.tt = zero_matrix(base_field(K), 1, degree(K))

    function _image(x::T)
      for i = 1:dim(A)
        z.t[1, i] = coeffs(x, copy = false)[i]
      end
      s = z.t*M
      return K(parent(K.pol)([ s[1, i] for i = 1:degree(K) ]))
    end

    function _preimage(x)
      for i = 1:degree(K)
        z.tt[1, i] = coeff(x, i - 1)
      end
      s = z.tt*N
      return A([ s[1, i] for i = 1:dim(A) ])
    end

    z.header = MapHeader{S, U}(A, K, _image, _preimage)
    return z
  end
end

function AbsAlgAssToNfAbsMor(A::S, K::U, M::V, N::V) where {S, U, V}
  return AbsAlgAssToNfAbsMor{S, elem_type(A), U, V}(A, K, M, N)
end

function _abs_alg_ass_to_nf_abs_mor_type(A::AbsAlgAss{T}) where {T}
  return AbsAlgAssToNfAbsMor{typeof(A), elem_type(A), _ext_type(T), dense_matrix_type(T)}
end

################################################################################
#
#  Morphisms between algebras and finite fields
#
################################################################################

# Morphism between an algebra A and a finite field Fq.
# base_ring(A) can be a GaloisField, a Generic.ResField{fmpz} or a Fq(Nmod)FiniteField, Fq can be a
# Fq(Nmod)FiniteField.
# MatType is the type of matrices over base_ring(A), PolyRingType the type of a
# polynomial ring over base_ring(A)
mutable struct AbsAlgAssToFqMor{S, T, MatType, PolyRingType} <: Map{S, T, HeckeMap, AbsAlgAssToFqMor}
  header::MapHeader{S, T}
  mat::MatType
  imat::MatType
  t::MatType # dummy vector used in image and preimage
  tt::MatType # another dummy vector
  R::PolyRingType
  RtoFq::FqPolyRingToFqMor # only used if S == AbsAlgAss{fq} or S == AbsAlgAss{fq_nmod}

  # First case: base_ring(A) == F_p
  function AbsAlgAssToFqMor{S, T, MatType, PolyRingType}(A::S, Fq::T, M::MatType, N::MatType, R::PolyRingType) where {
           S, T, MatType, PolyRingType
           #S <: AbsAlgAss{S1} where { S1 <: Union{ gfp_elem, Generic.ResF{fmpz} } },
           #T <: Union{ FqNmodFiniteField, FqFiniteField },
           #MatType <: Union{ gfp_mat, Generic.MatSpaceElem{Generic.ResF{fmpz}} },
           #PolyRingType <: Union{ GFPPolyRing, GFPFmpzPolyRing }
    }

    z = new{S, T, MatType, PolyRingType}()
    z.mat = M
    z.imat = N
    z.t = zero_matrix(base_ring(A), 1, dim(A))
    z.tt = zero_matrix(base_ring(A), 1, dim(A))
    z.R = R

    function _image(x::AbsAlgAssElem)
      @assert typeof(x) == elem_type(A)
      for i = 1:dim(A)
        z.t[1, i] = coeffs(x, copy = false)[i]
      end
      s = z.t*M
      sR = z.R([ s[1, i] for i = 1:dim(A) ])
      return Fq(sR)
    end

    function _preimage(x::Union{ fq_nmod, fq })
      @assert typeof(x) == elem_type(T)
      for i = 1:dim(A)
        z.tt[1, i] = base_ring(A)(coeff(x, i - 1))
      end
      s = z.tt*N
      return A([ s[1, i] for i = 1:dim(A) ])
    end

    z.header = MapHeader{S, T}(A, Fq, _image, _preimage)
    return z
  end

  # Second case: base_ring(A) == F_p^n
  function AbsAlgAssToFqMor{S, T, MatType, PolyRingType}(A::S, Fq::T, M::MatType, N::MatType, R::PolyRingType, RtoFq::FqPolyRingToFqMor) where {
                                                                                S, T, MatType, PolyRingType
           #S <: AbsAlgAss{S1} where { S1 <: Union{ fq, fq_nmod } },
           #T <: Union{ FqNmodFiniteField, FqFiniteField },
           #MatType <: Union{ fq_nmod_mat, fq_mat },
           #PolyRingType <: Union{ FqNmodPolyRing, FqPolyRing }
    }

    z = new{S, T, MatType, PolyRingType}()
    z.mat = M
    z.imat = N
    z.t = zero_matrix(base_ring(A), 1, dim(A))
    z.tt = zero_matrix(base_ring(A), 1, dim(A))
    z.R = R
    z.RtoFq = RtoFq

    function _image(x::AbsAlgAssElem)
      @assert typeof(x) == elem_type(A)
      for i = 1:dim(A)
        z.t[1, i] = coeffs(x, copy = false)[i]
      end
      s = z.t*M
      sR = z.R([ s[1, i] for i = 1:dim(A) ])
      return Fq(z.RtoFq(sR))
    end

    function _preimage(x::Union{ fq_nmod, fq })
      @assert typeof(x) == elem_type(T)
      x = z.RtoFq\x
      for i = 1:dim(A)
        z.tt[1, i] = base_ring(A)(coeff(x, i - 1))
      end
      s = z.tt*N
      return A([ s[1, i] for i = 1:dim(A) ])
    end

    z.header = MapHeader{S, T}(A, Fq, _image, _preimage)
    return z
  end

end

function AbsAlgAssToFqMor(A::AbsAlgAss{gfp_elem}, Fq::FqNmodFiniteField, M::gfp_mat, N::gfp_mat, R::GFPPolyRing)
  return AbsAlgAssToFqMor{typeof(A), FqNmodFiniteField, gfp_mat, GFPPolyRing}(A, Fq, M, N, R)
end

function AbsAlgAssToFqMor(A::AbsAlgAss{fq_nmod}, Fq::FqNmodFiniteField, M::fq_nmod_mat, N::fq_nmod_mat, R::FqNmodPolyRing, RtoFq::FqPolyRingToFqMor)
  return AbsAlgAssToFqMor{typeof(A), FqNmodFiniteField, fq_nmod_mat, FqNmodPolyRing}(A, Fq, M, N, R, RtoFq)
end

#function AbsAlgAssToFqMor(A::AbsAlgAss{Generic.ResF{fmpz}}, Fq::FqFiniteField, M::Generic.MatSpaceElem{Generic.ResF{fmpz}}, N::Generic.MatSpaceElem{Generic.ResF{fmpz}}, R::GFPFmpzPolyRing)
function AbsAlgAssToFqMor(A, Fq::FqFiniteField, M, N, R::GFPFmpzPolyRing)
  return AbsAlgAssToFqMor{typeof(A), FqFiniteField, typeof(M), GFPFmpzPolyRing}(A, Fq, M, N, R)
end

function AbsAlgAssToFqMor(A::AbsAlgAss{fq}, Fq::FqFiniteField, M::fq_mat, N::fq_mat, R::FqPolyRing, RtoFq::FqPolyRingToFqMor)
  return AbsAlgAssToFqMor{typeof(A), FqFiniteField, fq_mat, FqPolyRing}(A, Fq, M, N, R, RtoFq)
end

################################################################################
#
#  Maps to orders
#
################################################################################

function (f::AbsAlgAssMor)(O::AlgAssAbsOrd)
  domain(f) != algebra(O) && throw(error("Order not an order of the domain"))
  B = codomain(f)
  C = Order(B, elem_type(B)[f(b) for b in O.basis_alg])
  return C
end

function (f::AbsAlgAssMor)(I::AlgAssAbsOrdIdl)
  domain(f) != algebra(I) && throw(error("Order not an order of the domain"))
  B = codomain(f)
  J = ideal_from_lattice_gens(B, elem_type(B)[f(b) for b in basis(I)])
  return J
end


text:
"""
    BoolDistribution

A distribution that provides the probabilities of true or false.
Can construct with `BoolDistribution(p_true)`.
"""
struct BoolDistribution
    p::Float64 # probability of true
end

pdf(d::BoolDistribution, s::Bool) = s ? d.p : 1.0-d.p

rand(rng::AbstractRNG, d::BoolDistribution) = rand(rng) <= d.p

iterator(d::BoolDistribution) = [true, false]

==(d1::BoolDistribution, d2::BoolDistribution) = d1.p == d2.p

Base.hash(d::BoolDistribution) = hash(d.p)

Base.length(d::BoolDistribution) = 2


text:
line_list = []
append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 3, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 3, Y = -1im / 0.02)])

node_list = []
append!(node_list, [SlackAlgebraic(U = 1, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = -0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = -0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
#!/usr/bin/env julia

# Draw some text to a PNG image.

using Compose

img = PNG("text.png", 400px, 400px)
c = compose(compose(context(), text(150px, 200px, "hello & goodbye")), fill("tomato"))
draw(img, c)


text:
module BetaTest

using Test
using ReactiveMP
using Distributions
using Random

import ReactiveMP: mirrorlog

@testset "Beta" begin

    # Beta comes from Distributions.jl and most of the things should be covered there
    # Here we test some extra ReactiveMP.jl specific functionality

    @testset "vague" begin
        d = vague(Beta)

        @test typeof(d) <: Beta
        @test mean(d) === 0.5
        @test params(d) === (1.0, 1.0)
    end

    @testset "prod" begin
        @test prod(ProdAnalytical(), Beta(3.0, 2.0), Beta(2.0, 1.0)) ≈ Beta(4.0, 2.0)
        @test prod(ProdAnalytical(), Beta(7.0, 1.0), Beta(0.1, 4.5)) ≈ Beta(6.1, 4.5)
        @test prod(ProdAnalytical(), Beta(1.0, 3.0), Beta(0.2, 0.4)) ≈ Beta(0.19999999999999996, 2.4)
    end

    @testset "mean(::typeof(log))" begin
        @test mean(log, Beta(1.0, 3.0)) ≈ -1.8333333333333335
        @test mean(log, Beta(0.1, 0.3)) ≈ -7.862370395825961
        @test mean(log, Beta(4.5, 0.3)) ≈ -0.07197681436958758
    end

    @testset "mean(::typeof(mirrorlog))" begin
        @test mean(mirrorlog, Beta(1.0, 3.0)) ≈ -0.33333333333333337
        @test mean(mirrorlog, Beta(0.1, 0.3)) ≈ -0.9411396776150167
        @test mean(mirrorlog, Beta(4.5, 0.3)) ≈ -4.963371962929249
    end

end

end


text:
# broadcasting vector y(w) to a N-dim array x(w,k,m,l...)
# y is expected to be a vector to broadcast over x
# this kernel tries to take advantage of data reuse of broadcast operation
# y vector devided by BLOCK_SIZE number of elements between each thread block,
# each thread block iterates over second dimension of matrix,
#so from (1,1) to (BLOCK_SIZE,k) belongs to first thread block
# if number of thread block

#performance limitations
#if size of vector y smaller than BLOCK_SIZE_x,
#(BLOCK_SIZE_y*(BLOCK_SIZE_x-size(y))) number of thread from each block will be idle

#                         for high performance we need > 42 thread block (14 warp*3)
# thread_block_count = size(y)/BLOCK_SIZE_y > 42
# size(y) > BLOCK_SIZE_y*42
# with BLOCK_SIZE=32, broadcast_dim_size should be > 1344
# so if we have less than 1344 elements in vector y, performance will sour
# we can call lower limit 672, for better optimisation,(TODO-enis) I will make performance tests
# worst 448
# (sliding factor is introduced for this problem)

# this kernel can handle all vector and matrix dimensions

# explanation of kernel code is not added to prevent increase size of cuda14.cu
# (TODO-enis) provide a link to explanation of kernel index calculations for development

# shufle is slow due to index Access pattern, so removed

fp = open("cuda14.cu","w")
#using Knet: binary_ops

function cuda14src(f, j=f, ex="$f(xi,yi)")
  sprint() do s
    print(s,"#define BLOCK_SIZE_x 32\n#define BLOCK_SIZE_y 32\n")
    for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
        print(s,

"""

__global__ void _$(F)_14_x_y($T *x, $T *y,$T *z, int firstdimsize, int x_N,int sf)
{
    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.y;
    int ty = threadIdx.x;

    //row and col of matrix that thread will work on
    int row = by * blockDim.y * (sf+1) + ty;
    int col = bx * blockDim.y + tx;

    __shared__ $T Ys[BLOCK_SIZE_y];

    while((col)<firstdimsize)
    {
      if( ty==0 )
      {
          Ys[tx]=y[col];
      }
      __syncthreads();
      int Start = (row * firstdimsize) + col;
      int Step = firstdimsize * BLOCK_SIZE_y;

      for (int k= 0; k<(sf+1); k++)
      {
          if (Start<x_N)
          {
            $T xi = x[Start];
            $T yi = Ys[tx];
            z[Start]=$ex;
            Start+=Step;
          }
      }
      col += BLOCK_SIZE_x*gridDim.x;
    }
}

extern "C" {
  $DLLEXPORT void $(F)_14_x_y($T *x,$T *y,$T *z, int firstdimsize, int x_N,int flat_dims) {
    //int flat_dims=x_N/firstdimsize;
    int sf;
    if(flat_dims<129)
      {sf=4;}
    else if(flat_dims<513)
      {sf=8;}
    else
      {sf=16;}
    int thread_block_number4cover_x = (firstdimsize+BLOCK_SIZE_x-1)/BLOCK_SIZE_x;
    int thread_block_number4cover_y = (flat_dims+BLOCK_SIZE_y-1)/BLOCK_SIZE_y;
    int thread_blockwsliding = (thread_block_number4cover_y+(sf)-1)/(sf);
    dim3 dimGrid(thread_block_number4cover_x, thread_blockwsliding);
    dim3 dimBlock(BLOCK_SIZE_x, BLOCK_SIZE_y);
    //x_N size of the x
    //firstdimsize is size of y
    _$(F)_14_x_y<<<dimGrid,dimBlock>>>(x,y,z,firstdimsize,x_N,sf);
  }

  $DLLEXPORT void $(F)_14_x_y_stream($T *x,$T *y,$T *z, int firstdimsize, int x_N,int flat_dims, cudaStream_t STR) {
    //int flat_dims=x_N/firstdimsize;
    int sf;
    if(flat_dims<129)
      {sf=4;}
    else if(flat_dims<513)
      {sf=8;}
    else
      {sf=16;}
    int thread_block_number4cover_x = (firstdimsize+BLOCK_SIZE_x-1)/BLOCK_SIZE_x;
    int thread_block_number4cover_y = (flat_dims+BLOCK_SIZE_y-1)/BLOCK_SIZE_y;
    int thread_blockwsliding = (thread_block_number4cover_y+(sf)-1)/(sf);
    dim3 dimGrid(thread_block_number4cover_x, thread_blockwsliding);
    dim3 dimBlock(BLOCK_SIZE_x, BLOCK_SIZE_y);
    //x_N size of the x
    //firstdimsize is size of y
    _$(F)_14_x_y<<<dimGrid,dimBlock,0,STR>>>(x,y,z,firstdimsize,x_N,sf);
  }
}
""")
    end

    for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
        print(s,

"""

__global__ void _$(F)_14_y_x($T *x, $T *y,$T *z, int firstdimsize, int x_N,int sf)
{
    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.y;
    int ty = threadIdx.x;

    //row and col of matrix that thread will work on
    int row = by * blockDim.y * (sf+1) + ty;
    int col = bx * blockDim.y + tx;

    __shared__ $T Ys[BLOCK_SIZE_y];

    while((col)<firstdimsize)
    {
      if( ty==0 )
      {
          Ys[tx]=y[col];
      }
      __syncthreads();
      int Start = (row * firstdimsize) + col;
      int Step = firstdimsize * BLOCK_SIZE_y;

      for (int k= 0; k<(sf+1); k++)
      {
          if (Start<x_N)
          {
            $T yi = x[Start];
            $T xi = Ys[tx];
            z[Start]=$ex;
            Start+=Step;
          }
      }
      col += BLOCK_SIZE_x*gridDim.x;
    }
}

extern "C" {
  $DLLEXPORT void $(F)_14_y_x($T *x,$T *y,$T *z, int firstdimsize, int x_N,int flat_dims) {
    //int flat_dims=x_N/firstdimsize;
    int sf;
    if(flat_dims<129)
      {sf=4;}
    else if(flat_dims<513)
      {sf=8;}
    else
      {sf=16;}
    int thread_block_number4cover_x = (firstdimsize+BLOCK_SIZE_x-1)/BLOCK_SIZE_x;
    int thread_block_number4cover_y = (flat_dims+BLOCK_SIZE_y-1)/BLOCK_SIZE_y;
    int thread_blockwsliding = (thread_block_number4cover_y+(sf)-1)/(sf);
    dim3 dimGrid(thread_block_number4cover_x, thread_blockwsliding);
    dim3 dimBlock(BLOCK_SIZE_x, BLOCK_SIZE_y);
    //x_N size of the x
    //firstdimsize is size of y
    _$(F)_14_y_x<<<dimGrid,dimBlock>>>(x,y,z,firstdimsize,x_N,sf);
  }

  $DLLEXPORT void $(F)_14_y_x_stream($T *x,$T *y,$T *z, int firstdimsize, int x_N,int flat_dims, cudaStream_t STR) {
    //int flat_dims=x_N/firstdimsize;
    int sf;
    if(flat_dims<129)
      {sf=4;}
    else if(flat_dims<513)
      {sf=8;}
    else
      {sf=16;}
    int thread_block_number4cover_x = (firstdimsize+BLOCK_SIZE_x-1)/BLOCK_SIZE_x;
    int thread_block_number4cover_y = (flat_dims+BLOCK_SIZE_y-1)/BLOCK_SIZE_y;
    int thread_blockwsliding = (thread_block_number4cover_y+(sf)-1)/(sf);
    dim3 dimGrid(thread_block_number4cover_x, thread_blockwsliding);
    dim3 dimBlock(BLOCK_SIZE_x, BLOCK_SIZE_y);
    //x_N size of the x
    //firstdimsize is size of y
    _$(F)_14_y_x<<<dimGrid,dimBlock,0,STR>>>(x,y,z,firstdimsize,x_N,sf);
  }
}
""")
    end
  end
end

for a in binary_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda14src(a...))
end

close(fp)


text:
function house!(h, j, X, k)
    T = eltype(h)
    n = size(X,1)
    σ = zero(T)
    @inbounds @simd for i in j+1:n
        σ += X[i,k]^2
    end
    @inbounds begin
        h[j] = one(T)
        @simd for i in j+1:n
            h[i] = X[i,k]
        end
    end
    if σ == zero(T)
        β = zero(T)
    else
        @inbounds begin
            μ = sqrt(X[j,k]^2 + σ)
            if X[j,k] <= zero(T)
                h[j] = X[j,k] - μ
            else
                h[j] = -σ/(X[j,k] + μ)
            end
            β = 2*(h[j]^2) / (σ + h[j]^2)
            d = h[j]
            @simd for i in j:n
                h[i] /= d
            end
        end
    end
    β
end

"""
accumulate!(Q, A, β)

Q: Matrix to be over-written with Q, gets returned
A: Has the Householder vectors on the lower triangle
β: Has the Householder coefficients
"""
function accumulate!(Q, A, β)
    T = eltype(Q)
    m, n = size(Q)
    @inbounds begin
        Q .= zero(T)
        @simd for i in 1:n
            Q[i,i] = one(T)
        end
    end
    @inbounds for j in n:-1:1
        @simd for t in j:n
            dp = Q[j,t]
            for i in j+1:m
                dp += A[i,j]*Q[i,t]
            end
            Q[j,t] -= β[j]*dp
            for i in j+1:m
                Q[i,t] -= β[j]*dp*A[i,j]
            end
        end
    end
    Q
end

"""
getR!(R, A)

Copies the right triangular part of A into the right triangular part of R.
"""
function getR!(R, A)
    T = eltype(R)
    m, n = size(A)
    @inbounds @simd for j in 1:n
        for i in 1:j
            R[i,j] = A[i,j]
        end
    end
    R
end

"""
qrfact!(A, β, h=zeros(eltype(A), size(A,1)))

A: Matrix to be QR factorized, returned by the function
β: Empty buffer to store the Householder coefficients, returned by the function
h: Empty buffer for in-place operations
"""
function qrfact!(A, β, h=zeros(eltype(A), size(A,1)))
    m, n = size(A)
    @inbounds for j in 1:n
        β[j] = house!(h, j, A, j)
        @simd for k in j:n
            dp = A[j,k]
            for i in j+1:m
                dp += h[i]*A[i,k]
            end
            for l in j:m
                A[l,k] -= (β[j] * dp) * h[l]
            end
        end
        if j < m
            @simd for k in 1:(m-j)
                A[k+j,j] = h[k+j]
            end
        end
    end
    A, β
end


text:
using SphericalCodesNMR
using Test

@testset "SphericalCodesNMR.jl" begin
    # Write your tests here.
end

max_stars_repo_name: itscoe/SphericalCodesNMR

max_issues_repo_name: itscoe/SphericalCodesNMR

max_forks_repo_name: itscoe/SphericalCodesNMR


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule rocminfo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("rocminfo")
JLLWrappers.@generate_main_file("rocminfo", UUID("5a766526-3cf8-5128-8c31-4f7b7ad60f0e"))
end  # module rocminfo_jll


text:
using Distributions
using OBC
import OBC: propose, energy, reject
using Base.Test

########### Test Sampler ########
type TestParams
    mu :: Float64
    sigma :: Float64
end
TestParams() = TestParams(0.0,1.0)

function copy!(to::TestParams, from::TestParams)
    to.mu = from.mu
    to.sigma = from.sigma
end

type TestSampler <: Sampler
    curr :: TestParams
    old :: TestParams
    data :: Vector{Float64}
end
TestSampler() = TestSampler(TestParams(), TestParams(), [0.1, 0.1, 0.1]*3)

function propose(s::TestSampler)
    copy!(s.old, s.curr)
    s.curr.mu += randn()*0.2
    s.curr.sigma += randn()*0.2
    s.curr.sigma = clamp(s.curr.sigma, 0.1, Inf)
    0
end

function energy(s::TestSampler)
    -sum(logpdf(Normal(s.curr.mu, s.curr.sigma), s.data)) #Energies are -logp
end

function reject(s::TestSampler)
    copy!(s.curr, s.old)
end
########### /Test Sampler ########

s = TestSampler()
#x = set_energy_limits(s)
mymh = MHRecord(s, verbose=true)
sample(mymh, int(2e4))

@test length(mymh.db) == 20000
@test eltype(mymh.db) == Any
@test typeof(mymh.db[1]) == TestParams


text:
include("current_search.jl")
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using Evolutionary, Test, Random
using Distributed
using SharedArrays
using SignalAnalysis

function get_trains(p)
    fire = p.records[:fire]
    x, y = Float32[], Float32[]
    for time in eachindex(fire)
        for neuron_id in findall(fire[time])
            push!(x, time)
            push!(y, neuron_id)
        end
    end
    cellsa = Array{Union{Missing,Any}}(undef, 1, Int(findmax(y)[1]))
    nac = Int(findmax(y)[1])
    for (inx, cell_id) in enumerate(1:nac)
        cellsa[inx] = []
    end
    @inbounds for cell_id in unique(y)
        @inbounds for (time, cell) in collect(zip(x, y))
            if Int(cell_id) == cell
                append!(cellsa[Int(cell_id)], time)

            end

        end
    end

    cellsa

end

function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::GA, options)
    idx = sortperm(state.fitpop)
    record["fitpop"] = state.fitpop[:]#idx[1:last(idx)]]
    record["pop"] = population[:]
end

###
# Code appropriated from:
# https://github.com/JuliaML/MLPlots.jl/blob/master/src/optional/onlineai.jl
###
function get_ranges(ranges)

    lower = Float32[]
    upper = Float32[]
    for (k,v) in ranges
        append!(lower,v[1])
        append!(upper,v[2])
    end
    lower,upper
end

function init_b(lower,upper)
    gene = []#Float32
    for (i,(l,u)) in enumerate(zip(lower,upper))
        p1 = rand(l:u, 1)
        append!(gene,p1)
    end
    gene
end

function initf(n)
    genesb = []#Float32[]
    for i in 1:n
        genes = init_b(lower,upper)
        append!(genesb,[genes])
    end
    genesb
end

function custom_raster(P1::Array,P2::Array)
    y0 = Int32[0]
    X = Float32[];
	Y = Float32[]
	N=1
	y=0
	for p in [P1,P2]
		y+=1
		append!(X, P1)
		append!(Y, 1 .+ sum(y0))
    	push!(y0, 1)
		append!(X, P2)
		append!(Y, 1 .+ sum(y0))
	    push!(y0, 1)
		N+=1
    end
    plt = scatter(X, Y,w=50, m = (1, :black), leg = :none, marker=:vline,
                  xaxis=("t", (0, Inf)), yaxis = ("neuron",))
    y0 = y0[2:end-1]
    !isempty(y0) && hline!(plt, cumsum(y0), linecolor = :red)
    return plt
end

function get_vm(p,simulation_duration)
    vm = p.records[:v]
	vm = [i[1] for i in vm]
	vm = signal(vm, length(vm)/simulation_duration)
	vm
end

function checkmodel(param,cell_type,ngt_spikes)
    if cell_type=="IZHI"
        pp = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(;N = 1, param = pp)
    end
    if cell_type=="ADEXP"
		adparam = SNN.ADEXParameter(;a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10])
			E = SNN.AD(;N = 1, param=adparam)
    end

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
	current = current_search(cell_type,param,ngt_spikes)
    E.I = [current*nA]#_search(param,ngt_spikes)*nA]

    SNN.monitor(E, [:v])
	SNN.monitor(E, [:fire])

	###
	# TODO buid t into neuron models
	# SNN.monitor(E, [:t])
    ###
	simulation_duration= ALLEN_DURATION+ALLEN_DELAY+343ms
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = simulation_duration)
	spikes = get_spikes(E)
    spikes = [s/1000.0 for s in spikes]
	vm = get_vm(E,simulation_duration)
	#vm = E.records[:v]
	#vm = [i[1] for i in vm]

	return (vm,spikes)
end
function Evolutionary.value!(::Val{:multiproc}, fitness, objfun, population::AbstractVector{IT}) where {IT}
	fitness = SharedArrays.SharedArray{Float32}(fitness)
	@time @sync @distributed for i in 1:length(population)
		fitness[i] = value(objfun, population[i])
	end
	#@show(fitness)
	fitness
end
function Evolutionary.value!(::Val{:serial}, fitness, objfun, population::AbstractVector{IT}) where {IT}
	Threads.@threads for i in 1:length(population)
		fitness[i] = value(objfun, population[i])
	end
end
function get_data()
	#print(pwd())
	file="../JLD/ground_truth.jld"
    if false#isfile(file)
        vmgtv = load(file,"vmgtv")
        ngt_spikes = load(file,"ngt_spikes")
        gt_spikes = load(file,"gt_spikes")

        ground_spikes = gt_spikes
        ngt_spikes = size(gt_spikes)[1]
        vmgtt = load(file,"vmgtt")

        plot(plot(vmgtv,vmgtt,w=1))

    else
        py"""
        from neo import AnalogSignal
        from neuronunit.allenapi import make_allen_tests_from_id

        """

        py"""
        specimen_id = (
            325479788,
            324257146,
            476053392,
            623893177,
            623960880,
            482493761,
            471819401
        )
        specimen_id = specimen_id[1]
        target_num_spikes=7
        sweep_numbers, data_set, sweeps = make_allen_tests_from_id.allen_id_to_sweeps(specimen_id)
        (vmm,stimulus,sn,spike_times) = make_allen_tests_from_id.get_model_parts_sweep_from_spk_cnt(
            target_num_spikes, data_set, sweep_numbers, specimen_id
        )
        """
        gt_spikes = py"spike_times"
        ground_spikes = gt_spikes

        ngt_spikes = size(gt_spikes)[1]
        vmgtv = py"vmm.magnitude"
        vmgtt = py"vmm.times"

		#s_a = signal(opt_vec, length(opt_vec)/3.5)
		julia_version_vm = signal(vmgtv, length(vmgtt)/last(vmgtt))

		plot(plot(vmgtv,vmgtt,w=1))

        save("ground_truth.jld", "vmgtv", vmgtv,"vmgtt",vmgtt, "ngt_spikes", ngt_spikes,"gt_spikes",gt_spikes,"julia_version_vm",julia_version_vm)
        #filename = string("ground_truth: ", py"target_num_spikes")#,py"specimen_id)
        #filename = string(filename,py"specimen_id")
        #filename = string(filename,".jld")
        #save(filename, "vmgtv", vmgtv,"vmgtt",vmgtt, "ngt_spikes", ngt_spikes,"gt_spikes",gt_spikes)

    end

        vmgtv = load("ground_truth.jld","vmgtv")
        ngt_spikes = load("ground_truth.jld","ngt_spikes")
        ngt_spikes = size(gt_spikes)[1]
        ground_spikes = load("ground_truth.jld","gt_spikes")

        vmgtt = load("ground_truth.jld","vmgtt")
		julia_version_vm = load("ground_truth.jld","julia_version_vm")
    return (vmgtv,vmgtt,ngt_spikes,ground_spikes,julia_version_vm)
end
function get_izhi_ranges()
    ranges_izhi = DataStructures.OrderedDict{Char,Float32}()
    ranges_izhi = ("a"=>(0.002,0.3),"b"=>(0.02,0.36),"c"=>(-75,-35),"d"=>(0.005,16))#,"I"=>[100,9000])
    lower,upper = get_ranges(ranges_izhi)
    return lower,upper
end
#=
function initd()
    population = initf(50)
    garray = zeros((length(population)[1], length(population[1])))
    for (i,p) in enumerate(population)
        garray[i,:] = p
    end
    garray[1,:]
end
=#
function get_adexp_ranges()
    ranges_adexp = DataStructures.OrderedDict{String,Tuple{Float32,Float32}}()
    ranges_adexp[:"a"] = (2.0, 10)
    ranges_adexp[:"b"] = (5.0, 10)
    ranges_adexp[:"cm"] = (700.0, 983.5)
    ranges_adexp[:"v0"] = (-70, -55)
    ranges_adexp[:"τ_m"] = (10.0, 42.78345)
    ranges_adexp[:"τ_w"] = (300.0, 454.0)  # Tau_w 0, means very low adaption
    ranges_adexp[:"θ"] = (-45.0,-10)
    ranges_adexp[:"delta_T"] = (1.0, 5.0)
    ranges_adexp[:"v_reset"] = (-70.0, -15.0)
    ranges_adexp[:"spike_delta"] = (1.25, 20.0)
    lower,upper = get_ranges(ranges_adexp)
    return lower,upper
end

function vecplot(p, sym)
    v = SNN.getrecord(p, sym)
    y = hcat(v...)'
    x = 1:length(v)
    plot(x, y, leg = :none,
    xaxis=("t", extrema(x)),
    yaxis=(string(sym), extrema(y)))
end

function vecplot(P::Array, sym)
    plts = [vecplot(p, sym) for p in P]
    N = length(plts)
    plot(plts..., size = (600, 400N), layout = (N, 1))
end

function spike_train_difference(spkd_ground, spkd_found)
	#@show(spkd_ground)
	#@show(spkd_found)
	if length(spkd_found)==0
		return sum(ones(Int(length(spkd_ground))))
	end
	if length(spkd_found)>0
	    maxi0 = findmax(spkd_ground[:,:])[1]
		#@show(maxi0)
	    maxi1 = findmax(spkd_found[:,:])[1]
		#@show([maxi0, maxi1])
		maxi1 = maxi1[1]
		maxi0 = maxi0[1]
		#@show(maxi1)
		#@show(maxi0)

		mini = findmin([maxi0, maxi1])[1]
		#@show(mini)
	    #spkd = ones(Int(mini)#SharedArrays.SharedArray{Float32}(mini)
	    #maxi = findmax([maxi0, maxi1])[1]
		l0 = length(spkd_ground)
		l1 = length(spkd_found)
		#@show(spkd_ground)
		#s0 = sum(spkd_ground[1][:,:])
		#s1 = sum(spkd_found[1][:,:])
		#@show(s0)

		maxi = findmin([l0,l1])[1]
		spkd = zeros(Int(maxi))
		#@show(spkd)
		#@show(maxi)
		#=
		if maxi > 0
			if maxi0 != maxi1
				println("gets here a")
				return sum(ones(Int(mini)))
			end
			if isempty(spkd_found[1, :])
				println("gets here b")

				return sum(ones(Int(mini)))
			end
		end
		=#
		#println(length(spkd))
		#println(length(spkd_ground))
		#println(length(spkd_found))
	    @inbounds for (_, i) in zip(spkd, eachindex(spkd))
	        if !isempty(spkd_ground[i]) && !isempty(spkd_found[i])

	            maxt1 = findmax(spkd_ground[i])[1]
	            maxt2 = findmax(spkd_found[i])[1]
	            maxt = findmax([maxt1, maxt2])[1]

	            if maxt1 > 0.0 && maxt2 > 0.0
	                t, S = SpikeSynchrony.SPIKE_distance_profile(
	                    unique(sort(spkd_ground[i])),
	                    unique(sort(spkd_found[i]));
	                    t0 = 0.0,
	                    tf = maxt,
	                )
	                #b = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1]) # == SPIKE_distance(y1, y2)
	                spkd[i] = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1]) # == SPIKE_distance(y1, y2)

	            end
	        #end
			else
				spkd[i]=0.0
				#spkd = ones(Int(mini))
				#println("gets here c")

			end
	    end
	    scatter([i for i in 1:mini],spkd)|>display
		#e1=abs(l0-l1)
		#e2=abs(s0-s1)
		#@show(e1)
		#@show(e2)
		#@show(spkd)

		error = sum(spkd)#+abs(l0-l1)+abs(s0-s1)#+ sum(spkd)
	end
end

function initd()
    population = initd(10)
    garray = zeros((length(population)[1], length(population[1])))
    for (i, p) in enumerate(population)
        garray[i, :] = p
    end
    garray[1, :]
end

function initd(n)
    genesb = []
    lower = Float32[0.0 0.0 0.0 0.0]
    upper = Float32[1.0 1.0 1.0 1.0]
    lower = vec(lower)
    upper = vec(upper)

    for i = 1:n
        genes = initd(lower, upper)
        append!(genesb, [genes])
    end
    genesb
end

function initd(lower, upper)
    gene = []

    for (i, (l, u)) in enumerate(zip(lower, upper))
        p1 = rand(l:u, 1)
        append!(gene, p1)
    end
    gene
end


text:
using RecursiveArrayTools

@testset "NBodyProblem Test" begin
    G = 2.95912208286e-4
    M = [1.00000597682, 0.000954786104043, 0.000285583733151, 0.0000437273164546, 0.0000517759138449, 1 / 1.3e8]
    invM = inv.(M)
    planets = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"]

    pos_x = [0.0,-3.5023653,9.0755314,8.3101420,11.4707666,-15.5387357]
    pos_y = [0.0,-3.8169847,-3.0458353,-16.2901086,-25.7294829,-25.2225594]
    pos_z = [0.0,-1.5507963,-1.6483708,-7.2521278,-10.8169456,-3.1902382]
    pos = ArrayPartition(pos_x, pos_y, pos_z)

    vel_x = [0.0,0.00565429,0.00168318,0.00354178,0.00288930,0.00276725]
    vel_y = [0.0,-0.00412490,0.00483525,0.00137102,0.00114527,-0.00170702]
    vel_z = [0.0,-0.00190589,0.00192462,0.00055029,0.00039677,-0.00136504]
    vel = ArrayPartition(vel_x, vel_y, vel_z)

    tspan = (0., 200_000)

    ∑ = sum
    N = 6
    potential(p, t, x, y, z, M) = -G * ∑(i -> ∑(j -> (M[i] * M[j]) / sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2), 1:i - 1), 2:N)
    nprob = NBodyProblem(potential, M, vel, pos, tspan)
    sol = solve(nprob, Yoshida6(), dt=100)

    # Make sure the distances from the sun stay small enough
    f = (x, y, z, j) -> sqrt((x[1] - x[j])^2 + (y[1] - y[j])^2 + (z[1] - z[j])^2)
    for i in 1:length(sol)
        x, y, z = sol[i].x[1].x
        @test f(x, y, z, 2) < 6
        @test f(x, y, z, 3) < 10.5
        @test f(x, y, z, 4) < 21
        @test f(x, y, z, 5) < 32
        @test f(x, y, z, 6) < 50
    end
end


text:
struct Cluster
    name::String
    status::Symbol
    status_explanation::String
    num_jobs_running::Int32
    s3_bucket_arn::String
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
# API demo
# This is for the layout here: http://postimage.org/image/s6pskvknf/

require("LayoutLP")
using LayoutLP

# Create a fake window
win = Window(0, 0, 800, 600)

# Create the overall layout inside window "win". This is specified in terms of named tabstops.
xtabs = [:scatterL, :scatterR, :margeEL, :margeER, :imageL, :imageR]
ytabs = [:margeNT, :margeNB, :scatterT, :scatterB, :buttonsT]
l = LayoutLP(win, xtabs, ytabs)

# Now we'd probably do something like this:
cimage = canvas(l, [:imageL, :imageR], [:scatterT, :scatterB])
himage = image(cimage, "peppers.png")
# ... and then rely on height(himage) and width(himage) to give us the image dimensions.
# But because the canvas & image functionality is not yet worked out,
# let's do it this way:
imsz = [384,512]
height = () -> imsz[1]
width = () -> imsz[2]
# These are only nominally hard-coding the image size; you could change
# imsz[1] or imsz[2] and it would update to the new dimensions.

# For plot axes, the "outer" canvas include all labels, the "inner" canvas is the graphing area
outerscatter, innerscatter = plotaxis(l, [:scatterL, :scatterR], [:scatterT, :scatterB])
outermargE, innermargE = plotaxis(l, [:margeEL, :margeER], [:scatterT, :scatterB])
outermargN, innermargN = plotaxis(l, [:scatterL, :scatterR], [:margeNT, :margeNB])

# Create a sub-layout for the buttons
# The first two refer to parent coordinates, the second introduce tab stops in the child
lbuttons = LayoutLP(l, [:W, :E], [:buttonsT, :S], [:cancelL, :cancelR, :doneL], [])

# Create the image object inside the image canvas, so we have something to refer to in getting its height and width
himg = image(cimg)   # note: no image data yet, all we need now is the handle

# Now specify the geometry. This is done by adding constraints and penalties. Note that a "subplot" function would take care of many of these for you, but here we're going low-level because this GUI has a lot of custom requirements.

# First we require that the tabstops are in increasing order. Can also use "decreasing". One can pick out a subset of the tab stops and do something different for the rest (partial ordering)
increasing(l, l.xsyms)
l = increasing(l, l.ysyms)
# Now add particular constraints
l = addconstraints(l,
        # Ensure the proper aspect ratio for the image
        :(height($himg)*(imageR-imageL) = width($himg)*(scatterB-scatterT)),
        # Leave enough room for scatterplot y labels
        :(left($outerscatter) > W),
        :(left($outermargN) > W),
        # Scatterplot bottom label
        :(bottom($outerscatter) < buttonsT),
        # Ensure a small gap between east marginal axis and the image axis
        :(imageL > margeER + max(10px, 5mm))
    )
# Add penalties: these encourage desirable behavior, but are not strictly enforced
l = addpenalties(l,
        # we'd like the image to occupy approximately 40% of the horizontal space
        :(5*abs((imageR-imageL)-0.4*w)),
        # small gaps for the marginal axes
        :(10*abs(margeEL-scatterR-5px)),
        :(10*abs(margeEL-scatterR-5px)),
    )


text:
module AdventOfCode

include("TestUtils.jl")
include("Computer.jl")

include("Day01.jl")
include("Day02.jl")
include("Day03.jl")
include("Day04.jl")
include("Day05.jl")
include("Day06.jl")
include("Day07.jl")
include("Day08.jl")
include("Day09.jl")

end # module


text:
__precompile__()

module WaveletsExt

include("mod/Utils.jl")
include("mod/DWT.jl")
include("mod/TIWT.jl")
include("mod/SIWPD.jl")
include("mod/ACWT.jl")
include("mod/BestBasis.jl")
include("mod/SWT.jl")
include("mod/Denoising.jl")
include("mod/LDB.jl")
include("mod/Visualizations.jl")

using Reexport
@reexport using .DWT,
                .BestBasis,
                .Denoising,
                .Utils,
                .LDB,
                .TIWT,
                .SWT,
                .SIWPD,
                .ACWT,
                .Visualizations

end


text:
# Autogenerated wrapper script for opusfile_jll for x86_64-linux-musl
export libopusfile, libopusurl

using Ogg_jll
using Opus_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("opusfile")
JLLWrappers.@declare_library_product(libopusfile, "libopusfile.so.0")
JLLWrappers.@declare_library_product(libopusurl, "libopusurl.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Ogg_jll, Opus_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libopusfile,
        "lib/libopusfile.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopusurl,
        "lib/libopusurl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function reorder(t,s, strat::CommonVertex)

    # Find the permutation P of t and s that make
    # Pt = [P, A1, A2]
    # Ps = [P, B1, B2]
    I = zeros(Int,1)
    J = zeros(Int,1)
    e = 1
    for i in 1:3
        v = t[i]
        for j in 1:3
            w = s[j]
            if norm(w-v) < eps(eltype(v)) * 1.0e3
                I[e] = i
                J[e] = j
                e += 1
                break
            end
        end
        e == 2 && break
    end

    append!(I, setdiff([1,2,3], I))
    append!(J, setdiff([1,2,3], J))

    # # inverse permutations
    # K = indexin([1,2,3], I)
    # L = indexin([1,2,3], J)

    K = zeros(Int,3)
    for i in 1:3
        for j in 1:3
            if I[j] == i
                K[i] = j
                break
            end
        end
    end

    L = zeros(Int,3)
    for i in 1:3
        for j in 1:3
            if J[j] == i
                L[i] = j
                break
            end
        end
    end

    return I, J, K, L
end

function reorder(t, s, strat::CommonEdge)

    I = zeros(Int,3)
    J = zeros(Int,3)
    e = 1
    for i in 1:3
        v = t[i]
        for j in 1:3
            w = s[j]
            if norm(w-v) < eps(eltype(v)) * 1.0e3
                I[e] = i
                J[e] = j
                e += 1
                break
            end
        end
    end
    I[3] = setdiff([1,2,3], I[1:2])[1]
    J[3] = setdiff([1,2,3], J[1:2])[1]

    I = circshift(I,-1)
    J = circshift(J,-1)

    # # inverse permutations
    # K = indexin([1,2,3], I)
    # L = indexin([1,2,3], J)

    K = zeros(Int,3)
    for i in 1:3
        for j in 1:3
            if I[j] == i
                K[i] = j
                break
            end
        end
    end

    L = zeros(Int,3)
    for i in 1:3
        for j in 1:3
            if J[j] == i
                L[i] = j
                break
            end
        end
    end

    return I, J, K, L
end

function reorder(t, s, strat::CommonFace)

    I = [1,2,3]
    J = [1,2,3]
    K = [1,2,3]
    L = [1,2,3]

    return I, J, K, L
end


text:
using nfl_squares
using Test

@testset "nfl_squares.jl" begin
    # Write your tests here.
end


text:
using BoundaryValueDiffEq
using DiffEqBase, DiffEqDevTools
using Test

# First order test
function func_1!(du,u,p,t)
    du[1] = u[2]
    du[2] = 0
end

# Not able to change the initial condition.
# Hard coded solution.
func_1 = ODEFunction(func_1!, analytic=(u0, p, t) -> [5-t,-1])

function boundary!(residual, u, p, t)
    residual[1] = u[1][1]-5
    residual[2] = u[end][1]
end

function boundary_two_point!(residual, u, p, t)
    ua = u[1]
    ub = u[end]
    residual[1] = ua[1]-5
    residual[2] = ub[1]
end

# Second order linear test
function func_2!(du, u, p, t)
    du[1] = u[2]
    du[2] = -u[1]
end

# Not able to change the initial condition.
# Hard coded solution.
func_2!(::Type{Val{:analytic}}, u0, p, t) = [5*(cos(t) - cot(5)*sin(t)),
                                             5*(-cos(t)*cot(5) - sin(t))]
func_2 = ODEFunction(func_2!, analytic=(u0, p, t) -> [5*(cos(t) - cot(5)*sin(t)),
                                                      5*(-cos(t)*cot(5) - sin(t))])
tspan = (0.,5.)
u0 = [5.,-3.5]
probArr = [BVProblem(func_1, boundary!, u0, tspan),
           BVProblem(func_2, boundary!, u0, tspan),
           TwoPointBVProblem(func_1, boundary_two_point!, u0, tspan),
           TwoPointBVProblem(func_2, boundary_two_point!, u0, tspan)]

testTol = 0.2
affineTol = 1e-2
dts = (1/2) .^ (5:-1:1)
order = 4

println("Collocation method (GeneralMIRK)")
println("Affineness Test")
prob = probArr[1]
@time sol = solve(prob, GeneralMIRK4(), dt=0.2)
@test norm(diff(map(x->x[1], sol.u)) .+ 0.2, Inf) + abs(sol[1][1]-5) < affineTol

println("Convergence Test on Linear")
prob = probArr[2]
@time sim = test_convergence(dts,prob,GeneralMIRK4())
@test abs(sim.𝒪est[:final]-order) < testTol

println("Collocation method (MIRK)")
println("Affineness Test")
prob = probArr[3]
@time sol = solve(prob, MIRK4(), dt=0.2)
@test norm(diff(map(x->x[1], sol.u)) .+ 0.2, Inf) .+ abs(sol[1][1]-5) < affineTol

println("Convergence Test on Linear")
prob = probArr[4]
@time sim = test_convergence(dts,prob,MIRK4())
@test abs(sim.𝒪est[:final]-order) < testTol


text:
module ErgmRcall

"""
Helper functions to call useful functions from R ergm package.
For further information on ergm R package, please refer to
Statnet Development Team (Pavel N. Krivitsky, Mark S. Handcock, David R. Hunter, Carter T. Butts, Chad Klumb, Steven M. Goodreau, and Martina Morris) (2003-2020). statnet: Software tools for the Statistical Modeling of Network Data. URL http://statnet.org
"""
ErgmRcall

using DataFrames
using RCall
using SparseArrays
using Logging

function clean_start_RCall()
    # Clean R enviroment
    R"rm(list = ls())
    sessionInfo()"
   R"set.seed(0)"
end
export clean_start_RCall

function install_and_load_R_package(pkgName::String)
    R"""
        if( !( " $pkgName " %in% rownames(installed.packages())) ){
            install.packages("$pkgName")
        }

        library("$pkgName")
    """

end

function init_R_settings()
    # R"""options(warn=-1) """

    R"""local({r <- getOption("repos")
        r["CRAN"] <- "http://cran.r-project.org"
        options(repos=r)})
        """
end

function __init__()

    init_R_settings()
    # install_and_load_R_package("sna")
    # install_and_load_R_package("coda")
    # install_and_load_R_package("network")
    # install_and_load_R_package("ergm")
    R"""library("network")"""
    R"""library("ergm")"""

    clean_start_RCall()
end

"""
Sample from a sequence of ergms with possibly different parameters using R package ergm
"""
function sample_ergm_RCall_sequence(ergmTermsString, N, parDgpT_in, Nsample, mcmcBurnIn=100000)

    @debug "[sample_ergm_RCall_sequence][init][$ergmTermsString, N=$N, Nsample=$Nsample, size(parDgpT_in) = $(size(parDgpT_in)), parDgpT_in[:,1] = $(parDgpT_in[:,1]) ]"

    reval("formula_ergm = net ~ "*ergmTermsString)

    T = size(parDgpT_in)[2]
    # For each t, sample the ERGM with parameters corresponding to the DGP at time t
    parDgpT = Float64.(parDgpT_in)
    @rput T
    @rput parDgpT
    @rput N
    @rput Nsample
    @rput mcmcBurnIn
    #create an empty network, the formula defining ergm, sample the ensemble and store the sufficient statistics and change statistics in R
    samplingTime = @elapsed begin
    R"
        net <- network.initialize(N)
        sampledMat_T_R =    array(0, dim=c(N,N,T,Nsample))
        for(t in 1:T){
            changeStats_t_R = list()
            stats_t_R = list()
            # print(t)
            for(n in 1:Nsample){
                net <- simulate(formula_ergm, nsim = 1, seed = sample(1:100000000,1), coef = parDgpT[,t], control = control.simulate.formula(MCMC.burnin = mcmcBurnIn))
                sampledMat_T_R[,,t,n] <- as.matrix.network( net)
                }
            }"
    end

    # import sampled networks in julia
    sampledMat_T = Int8.(@rget(sampledMat_T_R))

    @debug "[sample_ergm_RCall_sequence][end][samplingTime=$samplingTime]"

    return sampledMat_T
    end
export sample_ergm_RCall_sequence

function get_edge_list(A::Matrix)
    @debug "[get_edge_list][begin]"

    if size(A)[2] == 2
        return A
    end

    size(A)[1] == size(A)[2] ? () : error("Adjacency matrix must be squared")
    spA = sparse(A)
    spA_lists = findnz(spA)
    @debug "[get_edge_list][end]"
    return hcat(spA_lists[1], spA_lists[2])
end

"""
Function that estimates a sequence of ergm defined by ergmTermsString (according to the notation of R package ergm)
"""
function get_one_mle(A::Matrix{T} where T<:Integer, ergmTermsString::String)
    @debug "[estimate_mle_RCall][begin]"

    edge_list = get_edge_list(A)
    @rput edge_list
    reval("formula_ergm = net ~ "*ergmTermsString)
    try
        R"
            net <- as.network.matrix(edge_list)

            tmp <- ergm(formula_ergm)#,estimate = 'MPLE')#)#

            estPar_R <- tmp[[1]]
        "
        estPar = @rget(estPar_R)
    catch
        @warn "[estimate_mle_RCall][mle failed. returning mple]"
        estPar = get_one_mple(A, ergmTermsString)
    end

    @debug "[estimate_mle_RCall][end]"
    return estPar
end
export get_one_mle

function get_change_stats(A::Matrix{T} where T<:Integer, ergmTermsString::String)
    # given a matrix returns the change statistics wrt to a given formula
    edge_list = get_edge_list(A)
    @rput edge_list

    reval("formula_ergm = net ~ "* ergmTermsString)
    # store the sufficient statistics and change statistics in R
    R"""
        net <- network(edge_list)
        chStat_t <- ergmMPLE(formula_ergm)
        changeStats_t_R <- cbind(chStat_t$response, chStat_t$predictor,     chStat_t$weights)
            """

    changeStats = @rget changeStats_t_R;# tmp = Array{Array{Float64,2}}(T); for
    return changeStats
end
export get_change_stats

""" Get a single pmle from adjacency matrix"""
function get_one_mple(A::Matrix{T} where T<:Integer, ergmTermsString::String)
    edge_list = get_edge_list(A)
    @rput edge_list
    reval("formula_ergm = net ~ "* ergmTermsString)

    with_logger(NullLogger()) do

        R"""
            net <- network(edge_list)
            mple_t_R <- ergmMPLE(formula_ergm, output="fit")$coef
            # print(mple_t_R)
            """
    end

    mple_t = @rget mple_t_R;# tmp = Array{Array{Float64,2}}(T); for
    return mple_t
end
export get_one_mple

""" Get a single pmle from change statistics for a single network"""
function get_one_mple(changeStats_T::Array{T,2} where T<:Real, ergmTermsString::String)
    # and run a single snapeshot estimate
    @rput changeStats_T
    reval("formula_ergm = net ~ "* ergmTermsString)
    with_logger(NullLogger()) do
        R"""
            mple_static_R <- glm(changeStats_T[,1] ~ . -1, data = data.frame(changeStats_T[,c(2:(ncol(changeStats_T)-1))]), weights =changeStats_T[,ncol(changeStats_T)],family="binomial")$coefficients
            """
    end
    mple_static = @rget mple_static_R;# tmp = Array{Array{Float64,2}}(T); for

    return mple_static
end

""" Get a single pmle from change statistics for a sequence of networks"""
function get_one_mple(changeStats_T::Array{Array{T,2},1} where T<:Real, ergmTermsString::String)
    return get_one_mple(reduce(vcat,changeStats_T), ergmTermsString)
end
export get_one_mple

function decomposeMPLEmatrix(ergmMPLEmatrix)
    changeStat = ergmMPLEmatrix[:,2:end-1]
    response = ergmMPLEmatrix[:,1]
    weights = ergmMPLEmatrix[:,4]
    return changeStat, response, weights
end
export decomposeMPLEmatrix

end


text:
using Test
using ResourceBundles
using ResourceBundles.CLocales

function localized_isless(loc::LocaleId)
    function lt(a::AbstractString, b::AbstractString)
        strcollate(a, b, loc) < 0
    end
end


text:
module TarIterators

export TarIterator

using Tar
using BoundedStreams

struct TarIterator{T,F<:Function}
    stream::T
    filter::F
    closestream::Bool
end

"""
    TarIterator(io[, predicate])

Iterator over an input stream open for reading. `io` may also be pathname of an existing
tar file. The data in the stream or file must obey the `tar` file format understood by
package `Tar`.

If `predicate` is given, only tar elements with matching header data are processed.
* `String` or `Regex`: `predicate` is applied to the element names.
* `Symbol`: it must coincide with field `type` of `Tar.Header`.
* boolean function: it is applied to the `Tar.Header` of the elements.
* `Tuple`: AND of all predicates contained in tuple
* `Vector`: OR of all predicates contained in vector

Examples:
```
    io = open("/tmp/abc.tar")
    TarIterator(io)              # all entries
    TarIterator(io, :file)       # all entries of file type
    TarIterator(io, "y")         # only entry with name "y"
    TarIterator(io, r".*[.]txt") # all entries with name ending in ".txt"
    TarIterator(io, h -> h.size > 100) # only entries with data size > 100
```
"""
function TarIterator(source::T, a::F=nothing; close_stream::Bool=false) where {T,F}
    TarIterator(source, selector(a), close_stream)
end
function TarIterator(file::AbstractString, f=nothing; close_stream::Bool=false)
    TarIterator(open(file), f, close_stream=close_stream)
end

function Base.iterate(ti::TarIterator, status=nothing)
    stream = ti.stream
    status !== nothing && close(status)
    h = Tar.read_header(stream)
    while h !== nothing
        s = align(h.size)
        if ti.filter(h)
            closeop = ti.closestream ? BoundedStreams.CLOSE : s
            io = BoundedInputStream(stream, h.size, close=closeop)
            return (h, io), io
        end
        skip(stream, s)
        h = Tar.read_header(stream)
    end
    nothing
end

# align to TAR block size
align(pos::Integer) = mod(-pos, 512) + pos

# predicates processing
selector(f::Function) = f
selector(a::Tuple) = function AND(h::Tar.Header)
    for f in a
        if !selector(f)(h)
            return false
        end
    end
    true
end
selector(a::AbstractVector) = function OR(h::Tar.Header)
    for f in a
        if selector(f)(h)
            return true
        end
    end
    false
end
selector(a::Union{AbstractString,Nothing,Symbol,Regex}) = h::Tar.Header -> predicate(h, a)
predicate(h::Tar.Header, ::Nothing) = true
predicate(h::Tar.Header, a::AbstractString) = h.path == a
predicate(h::Tar.Header, s::Symbol) = h.type == s
function predicate(h::Tar.Header, r::Regex)
    fn = h.path
    m = match(r, fn)
    m !== nothing && m.match == fn
end

"""
    open(ti::TarIterator)::BoundedInputStream

Skip to the first entry in tar stream according to predicates of `ti` and return
an open input stream, which allows to read data part of this entry.
"""
function Base.open(ti::TarIterator)
    s = iterate(ti)
    s === nothing && throw(EOFError())
    s[1][2]
end

"""
    open(ti::TarIterator) do h, io; ... end

Process all tar entries selected by `ti` and provide `h::Tar.Header` and
`io::BoundedInputStream` to the called function.
"""
function Base.open(f::Function, ti::TarIterator)
    for (h, io) in ti
        f(h, io)
    end
end

Base.close(ti::TarIterator) = close(ti.stream)
Base.seekstart(ti::TarIterator) = seekstart(ti.stream)

end # module


text:
struct FullVar end
struct BranchVar end

Base.@kwdef mutable struct GlobalOptimizer{N,T<:Real,S<:ExtensionType} <: MOI.AbstractOptimizer

        _ext_type::S = S()
        _stack::BinaryMinMaxHeap{NodeBB{N,T}} = BinaryMinMaxHeap{NodeBB{N,T}}()
        _current_node::NodeBB{N,T} = NodeBB{N,T}()

        _input_to_local_map::MOIU.IndexMap              = MOIU.IndexMap()
        _constraint_primal::Dict{CI,Union{T,Vector{T}}} = Dict{CI,T}()
        _constraint_offset::Vector{Int}                 = Int[]
        _constraint_index_num::Int = 0
        _constraint_row_num::Int = 0

        _working_problem::ParsedProblem{N,T} = ParsedProblem{N,T}()
        _input_problem::InputModel{T} = InputModel{T}()

        _first_relax_point_set::Bool = false
        _current_xref::Vector{T} = T[]
        _candidate_xref::Vector{T} = T[]

        _use_prior_objective_xref::Bool = false
        _current_objective_xref::Vector{T} = T[]
        _prior_objective_xref::Vector{T} = T[]

        _user_branch_variables::Bool = false
        _fixed_variable::Vector{Bool} = Bool[]
        _branch_variable_num::Int = 0
        _branch_to_variable_map::Vector{Int} = Int[]
        _variable_to_branch_map::Vector{Int} = Int[]

        _lower_variable_bound::Vector{T} = T[]
        _upper_variable_bound::Vector{T} = T[]

        # all subproblem immutable subproblem status are set in global_solve in corresponding routines
        # in optimize_nonconvex.jl
        _preprocess_feasibility::Bool = true
        _preprocess_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
        _preprocess_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED

        _lower_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
        _lower_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
        _lower_feasibility::Bool = true
        _lower_objective_value::T = -Inf

        # set in TODO
        _lower_solution::Vector{T} = T[]
        _lower_lvd::Vector{T}      = T[]
        _lower_uvd::Vector{T}      = T[]

        _cut_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
        _cut_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
        _cut_objective_value::T = -Inf
        _cut_feasibility::Bool = true

        # set in TODO
        _cut_solution::Vector{T} = T[]

        _upper_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
        _upper_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
        _upper_feasibility::Bool = true
        _upper_objective_value::T = Inf

        # array is initialized to correct size in TODO, reset in single_nlp_solve! in optimize_convex.jl
        _upper_variables::Vector{VI} =  VI[]

        # set in TODO
        _upper_solution::Vector{T} = T[]

        _postprocess_feasibility::Bool = true

        _time_left::T = 1000*one(T)
        _start_time::T = zero(T)
        _run_time::T = zero(T)
        _parse_time::T = zero(T)
        _presolve_time::T = zero(T)
        _last_preprocess_time::T = zero(T)
        _last_lower_problem_time::T = zero(T)
        _last_upper_problem_time::T = zero(T)
        _last_postprocessing_time::T = zero(T)

        # reset in initial_parse! in parse.jl
        _min_converged_value::T = typemax(T)
        _global_lower_bound::T = typemin(T)
        _global_upper_bound::T = typemax(T)
        _maximum_node_id::Int = 0
        _iteration_count::Int = 0
        _node_count::Int = 0

        # Storage for output, reset in initial_parse! in parse.jl
        _solution_value::T = zero(T)
        _feasible_solution_found::Bool = false
        _first_solution_node::Int = -1
        _best_upper_value::T = Inf

        # Optimality-Based Bound Tightening (OBBT) Options
        # set in TODO
        _obbt_working_lower_index::Vector{Bool} = Bool[]
        _obbt_working_upper_index::Vector{Bool} = Bool[]
        _lower_indx_diff::Vector{Bool} = Bool[]
        _upper_indx_diff::Vector{Bool} = Bool[]
        _old_low_index::Vector{Bool} = Bool[]
        _old_upp_index::Vector{Bool} = Bool[]
        _new_low_index::Vector{Bool} = Bool[]
        _new_upp_index::Vector{Bool} = Bool[]
        _obbt_variables::Vector{VI} = VI[]
        _obbt_variable_num::Int = 0
        _obbt_performed_flag::Bool = false

        # Feasibility-Based Bound Tightening Options
        # set in set_constraint_propagation_fbbt in domain_reduction.jl
        _cp_improvement::T = zero(T)
        _cp_evaluation_reverse::Bool = false

        _cut_iterations::Int = 0
        _cut_add_flag::Bool = false

        # Options for Repetition (If DBBT Performed Well)
        # set in within preprocess in optimize_nonconvex.jl
        _node_repetitions::Int = 0
        _initial_volume::T = zero(T)
        _final_volume::T = zero(T)

        # Log
        _log::Log = Log()

        # set in TODO
        _buffered_quadratic_ineq_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]
        _buffered_quadratic_eq_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]
        _buffered_nonlinear_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]

        # need to retreive primal _relaxed_variable_index
        # set in TODO
        #"Number of variables actively branched on in B&B routine (excludes linear and fixed)"
        _relaxed_variable_number::Int = 0
        _relaxed_variable_index::Vector{VI} = VI[]
        _relaxed_variable_eq::Vector{Tuple{CI{SV, ET}, Int}} = Tuple{CI{SV, ET}, Int}[]
        _relaxed_variable_lt::Vector{Tuple{CI{SV, LT}, Int}} = Tuple{CI{SV, LT}, Int}[]
        _relaxed_variable_gt::Vector{Tuple{CI{SV, GT}, Int}} = Tuple{CI{SV, GT}, Int}[]

        # set as user-input
        _branch_variables::Vector{Bool} = Bool[]

        _new_eval_constraint::Bool = false
        _new_eval_objective::Bool = false

        _node_to_sv_leq_ci::Vector{CI{SV,LT}} = CI{SV,LT}[]
        _node_to_sv_geq_ci::Vector{CI{SV,GT}} = CI{SV,GT}[]

        #"Set to true if a nonlinear evaluator was created (NLconstraint or NLobjective specified)"
        _nonlinear_evaluator_created::Bool = false

        #_relaxed_evaluator::Evaluator = Evaluator{1,NS}()
        #_relaxed_constraint_bounds::Vector{MOI.NLPBoundsPair} = Vector{MOI.NLPBoundsPair}[]
end
GlobalOptimizer(::Type{T}) where T<:AbstractFloat = GlobalOptimizer{1,T,DefaultExt}()

function MOI.is_empty(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    empty_opt = GlobalOptimizer{N,T,S}()
    is_empty_flag = true
    for f in fieldnames(GlobalOptimizer)
        if f == :_stack
            is_empty_flag &= isempty(getfield(m, f))
        elseif f == :_current_node
            is_empty_flag &= getfield(m, f) == NodeBB{N,T}()
        elseif f == :_log
                # TODO
        elseif f == :_working_problem
                # TODO
        elseif f == :_input_problem
        else
            if getfield(empty_opt, f) != getfield(m, f)
                is_empty_flag = false
                break
            end
        end
    end
    return is_empty_flag
end

#=
TODO: StructArray access is faster here for most...
=#

Base.@propagate_inbounds function _is_unfixed_integer(::FullVar, m::GlobalOptimizer, i)
end
Base.@propagate_inbounds function _is_unfixed_integer(::BranchVar, m::GlobalOptimizer, i)
end

for (f,arr) in ((:_lower_bound, :_lower_variable_bound),
                (:_upper_bound, :_upper_variable_bound),
                (:_lower_solution, :_lower_solution),)
    @eval function ($f)(::FullVar, m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S<:ExtensionType}
        m.$arr
    end
    @eval Base.@propagate_inbounds function ($f)(::FullVar, m::GlobalOptimizer{N,T,S}, i) where {N,T<:Real,S<:ExtensionType}
        m.$arr[i]
    end
    @eval Base.@propagate_inbounds function ($f)(::BranchVar, m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S<:ExtensionType}
        m.$arr[m._branch_variable]
    end
    @eval Base.@propagate_inbounds function ($f)(::BranchVar, m::GlobalOptimizer{N,T,S}, i) where {N,T<:Real,S<:ExtensionType}
        m.$arr[m._branch_variable[i]]
    end
end

Base.@propagate_inbounds function _set_lower_bound!(::BranchVar, m::GlobalOptimizer{N,T,S}, i, v) where {N,T<:Real,S<:ExtensionType}
    m._lower_variable_bound[m._branch_variable[i]] = v
end
Base.@propagate_inbounds function _set_lower_bound!(::FullVar, m::GlobalOptimizer{N,T,S}, i, v) where {N,T<:Real,S<:ExtensionType}
    m._lower_variable_bound[i] = v
end
Base.@propagate_inbounds function _set_upper_bound!(::BranchVar, m::GlobalOptimizer{N,T,S}, i, v) where {N,T<:Real,S<:ExtensionType}
    m._upper_variable_bound[m._branch_variable[i]] = v
end
Base.@propagate_inbounds function _set_upper_bound!(::FullVar, m::GlobalOptimizer{N,T,S}, i, v) where {N,T<:Real,S<:ExtensionType}
    m._upper_variable_bound[i] = v
end

function _mid(::FullVar, m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S<:ExtensionType}
    0.5*(_upper_bound(FullVar, m) - _lower_bound(FullVar, m))
end
Base.@propagate_inbounds function _mid(::FullVar, m::GlobalOptimizer{N,T,S}, i) where {N,T<:Real,S<:ExtensionType}
    0.5*(_upper_bound(FullVar, m, i) - _lower_bound(FullVar, m, i))
end
function _mid(::BranchVar, m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S<:ExtensionType}
end
Base.@propagate_inbounds function _mid(::BranchVar, m::GlobalOptimizer{N,T,S}, i) where {N,T<:Real,S<:ExtensionType}
end

_input_problem(m::GlobalOptimizer) = m._input_problem
_working_problem(m::GlobalOptimizer) = m._working_problem