text:
# How to create an extra prompt and attach it to the current REPL was learned
# by reading the source code of @Keno's CXX.jl

module Prompt

import REPL
import REPL.LineEdit
import REPL.Terminals

import REPL: respond, return_callback
import REPL.LineEdit: buffer, cmove_col, cmove_up, InputAreaState, transition,
                      terminal, buffer, on_enter, move_input_end, add_history, state, mode, edit_insert

import Tokenize.Lexers

import REPL.Terminals: raw!, width, height, cmove, getX, TerminalBuffer,
                  getY, clear_line, beep, disable_bracketed_paste, enable_bracketed_paste

using OhMyREPL
import OhMyREPL: untokenize_with_ANSI, apply_passes!, PASS_HANDLER

@nospecialize # use only declared type signatures

function rewrite_with_ANSI(s, cursormove::Bool = false)
        if isa(s, LineEdit.SearchState)
            return
        end
        # Clear input area
        p = position(buffer(s))
        if isa(s, LineEdit.PrefixSearchState)
            s = s.mi
        end
        if isa(s, LineEdit.MIState)
            mode = s.mode_state[s.current_mode]
        else
            mode = s
        end

        io = IOBuffer()
        outbuf = IOContext(io, stdout)
        termbuf = Terminals.TerminalBuffer(outbuf)
        # Hide the cursor
        LineEdit.write(outbuf, "\e[?25l")
        LineEdit.clear_input_area(termbuf, mode)
        # Extract the cursor index in character count
        cursoridx = length(String(buffer(s).data[1:p]))

        l = textwidth(get_prompt(s))
        if !isa(s, LineEdit.SearchState)
            # xref: https://github.com/JuliaLang/julia/pull/36689
            @static if VERSION ≥ v"1.6.0-DEV.517"
                LineEdit.write_prompt(termbuf, mode, LineEdit.hascolor(terminal(s)))
            else
                LineEdit.write_prompt(termbuf, mode)
            end
            LineEdit.write(termbuf, "\e[0m") # Reset any formatting from Julia so that we start with a clean slate
        end

        # Insert colorized text from running the passes
        seekstart(buffer(s))
        tokens = collect(Lexers.Lexer(buffer(s)))
        apply_passes!(PASS_HANDLER, tokens, cursoridx, cursormove)
        untokenize_with_ANSI(outbuf, PASS_HANDLER , tokens, l)

        # Reset the buffer since the Lexer messed with it (maybe the Lexer should reset it on done)
        seek(buffer(s), p)

        # Our cursor now seems to be out of place, we run the already existing refresh_multi_line code to put it where it belongs.
        # Maybe it is possible to save the cursor and just restore it but that is probably Terminal dependent...
        mode.ias = refresh_multi_line(termbuf, terminal(s), buffer(s), mode.ias, l)
        LineEdit.write(outbuf, "\e[?25h")  # Show the cursor
        write(terminal(s), take!(io))
        flush(terminal(s))
end

function create_keybindings()

    D = Dict{Any, Any}()
    D['\b']   = (s, data, c) -> if LineEdit.edit_backspace(buffer(s))
        rewrite_with_ANSI(s)
    else
        beep(terminal(s))
    end
    D["*"]    = (s, data, c) ->  (LineEdit.edit_insert(buffer(s), c); rewrite_with_ANSI(s))
    D["^B"]   = (s, data, c) -> (LineEdit.edit_move_left(buffer(s)) ;rewrite_with_ANSI(s))
    D["^F"]   = (s, data, c) -> (LineEdit.edit_move_right(buffer(s)) ;rewrite_with_ANSI(s))
    # Meta B
    D["\eb"]  = (s, data, c) -> (LineEdit.edit_move_word_left(s) ; rewrite_with_ANSI(s))
    # Meta F
     D["\ef"]  = (s, data, c) -> (LineEdit.edit_move_word_right(s); rewrite_with_ANSI(s))
    # Meta Enter
    D["\e\r"] = (s, data, c) -> (LineEdit.edit_insert(buffer(s), '\n'); rewrite_with_ANSI(s))
    D["^A"]   = (s, data, c) -> (LineEdit.move_line_start(s); rewrite_with_ANSI(s))
    D["^E"]   = (s, data, c) -> (LineEdit.move_line_end(s); rewrite_with_ANSI(s))
    D["\e[H"] = (s, data, c) -> (LineEdit.move_input_start(s); rewrite_with_ANSI(s))
    D["\e[F"] = (s, data, c) -> (LineEdit.move_input_end(buffer(s)); rewrite_with_ANSI(s))
    D["^L"]   = (s, data, c) -> (Terminals.clear(terminal(s)); rewrite_with_ANSI(s))
    D["^W"]   = (s, data, c) -> LineEdit.edit_werase(s)
    # Right Arrow
    D["\e[C"] = (s, data, c)->(LineEdit.edit_move_right(buffer(s)); rewrite_with_ANSI(s))
    # Left Arrow
    D["\e[D"] = (s, data, c)->(LineEdit.edit_move_left(buffer(s)); rewrite_with_ANSI(s))
    # Up Arrow
    # Delete
    D["\e[3~"] = (s, data, c)->(LineEdit.edit_delete(buffer(s)); rewrite_with_ANSI(s))
    D["^T"] = (s, data, c)->(LineEdit.edit_transpose_chars(buffer(s)); rewrite_with_ANSI(s))
    D["\et"] = (s, data, c)->(LineEdit.edit_transpose_words(buffer(s)); rewrite_with_ANSI(s))
    D["\ed"] = (s, data, c)->(LineEdit.edit_delete_next_word(buffer(s)); rewrite_with_ANSI(s))
    D["\e\b"] = (s, data, c)->(LineEdit.edit_delete_prev_word(buffer(s)); rewrite_with_ANSI(s))
    D["^D"] = (s, data, c)->begin
        if buffer(s).size > 0
            LineEdit.edit_delete(buffer(s)); rewrite_with_ANSI(s)
        else
            println(terminal(s))
            return :abort
        end
    end

    # Hack around a bit to make enter not remove syntax highlighting above
    D["\r"] = (s, data, c) -> begin
        if on_enter(s) || (eof(buffer(s)) && s.key_repeats > 1)
            # Disable bracket highlighting before entering
            brackidx = OhMyREPL._find_pass(OhMyREPL.PASS_HANDLER, "BracketHighlighter")
            brackstatus = false
            if brackidx != -1
                brackstatus = OhMyREPL.PASS_HANDLER.passes[brackidx][2].enabled
                OhMyREPL.enable_pass!(PASS_HANDLER, "BracketHighlighter", false)
            end
            _commit_line(s, data, c)
            if brackidx != -1 && brackstatus == true
                OhMyREPL.enable_pass!(PASS_HANDLER, "BracketHighlighter", true)
            end
            return :done
        else
            edit_insert(buffer(s), '\n')
            rewrite_with_ANSI(s)
        end
    end

    # Hack around a bit to make Ctrl + C work
    D["^C"] = (s, data, c) -> begin
        try # raise the debugger if present
            ccall(:jl_raise_debugger, Int, ())
        catch
        end
        move_input_end(s)
        rewrite_with_ANSI(s)
        print(terminal(s), "^C\n\n")
        transition(s, :reset)
        rewrite_with_ANSI(s)
    end

    # Fixup bracket paste a bit
    D["\e[200~"] = (s, data, c) ->begin
        input = LineEdit.bracketed_paste(s) # read directly from s until reaching the end-bracketed-paste marker
        sbuffer = LineEdit.buffer(s)
        curspos = position(sbuffer)
        seek(sbuffer, 0)
        shouldeval = (bytesavailable(sbuffer) == curspos && !occursin(UInt8('\n'), sbuffer))
        seek(sbuffer, curspos)
        if curspos == 0
            # if pasting at the beginning, strip leading whitespace
            input = lstrip(input)
        end
        if !shouldeval
            # when pasting in the middle of input, just paste in place
            # don't try to execute all the WIP, since that's rather confusing
            # and is often ill-defined how it should behave
            edit_insert(s, input)
            rewrite_with_ANSI(s)
            return
        end
        LineEdit.push_undo(s)
        edit_insert(sbuffer, input)
        input = String(take!(sbuffer))
        oldpos = firstindex(input)
        firstline = true
        isprompt_paste = false
        prompt = get_prompt(s)
        default_prompt = "julia> "
        jl_prompt_len = textwidth(prompt)
        jl_default_len = textwidth(default_prompt)
        while oldpos <= lastindex(input) # loop until all lines have been executed
            # Check if the next statement starts with "julia> ", in that case
            # skip it. But first skip whitespace
            while input[oldpos] in ('\n', ' ', '\t')
                oldpos = nextind(input, oldpos)
                oldpos >= sizeof(input) && return
            end
            # Check if input line starts with "julia> ", remove it if we are in prompt paste mode
            sstr = SubString(input, oldpos)
            match_default = startswith(sstr, default_prompt)
            match_prompt = startswith(sstr, prompt)
            if (firstline || isprompt_paste) && (match_default || match_prompt)
                isprompt_paste = true
                match_prompt ? (oldpos += jl_prompt_len) : (oldpos += jl_default_len)
            # If we are prompt pasting and current statement does not begin with julia> , skip to next line
            elseif isprompt_paste
                while input[oldpos] != '\n'
                    oldpos = nextind(input, oldpos)
                    oldpos >= sizeof(input) && return
                end
                continue
            end
            ast, pos = Meta.parse(input, oldpos, raise=false, depwarn=false)
            if (isa(ast, Expr) && (ast.head == :error || ast.head == :incomplete)) ||
                    (pos > ncodeunits(input) && !endswith(input, '\n'))
                # remaining text is incomplete (an error, or parser ran to the end but didn't stop with a newline):
                # Insert all the remaining text as one line (might be empty)
                tail = input[oldpos:end]
                if !firstline
                    # strip leading whitespace, but only if it was the result of executing something
                    # (avoids modifying the user's current leading wip line)
                    tail = lstrip(tail)
                end
                if isprompt_paste # remove indentation spaces corresponding to the prompt
                    tail = replace(tail, r"^"m * ' '^jl_prompt_len => "")
                end
                LineEdit.replace_line(s, tail, true)
                rewrite_with_ANSI(s)
                break
            end
            # get the line and strip leading and trailing whitespace
            line = strip(input[oldpos:prevind(input, pos)])
            if !isempty(line)
                if isprompt_paste # remove indentation spaces corresponding to the prompt
                    line = replace(line, r"^"m * ' '^jl_prompt_len => "")
                end
                # put the line on the screen and history
                LineEdit.replace_line(s, line)
                _commit_line(s, data, c)
                # execute the statement
                terminal = LineEdit.terminal(s) # This is slightly ugly but ok for now
                raw!(terminal, false) && disable_bracketed_paste(terminal)
                LineEdit.mode(s).on_done(s, LineEdit.buffer(s), true)
                raw!(terminal, true) && enable_bracketed_paste(terminal)
                LineEdit.push_undo(s) # when the last line is incomplete
            end
            oldpos = pos
            firstline = false
        end
    end

    # Tab
    D['\t'] = (s, data, c) -> begin
        buf = buffer(s)
        # Yes, we are ignoring the possiblity
        # the we could be in the middle of a multi-byte
        # sequence, here but that's ok, since any
        # whitespace we're interested in is only one byte
        i = position(buf)
        if i != 0
            c = buf.data[i]
            if c == UInt8('\n') || c == UInt8('\t') ||
               # hack to allow path completion in cmds
               # after a space, e.g., `cd <tab>`, while still
               # allowing multiple indent levels
               (c == UInt8(' ') && i > 3 && buf.data[i-1] == UInt8(' '))
                edit_insert(s, " "^4)
                return
            end
        end
        LineEdit.complete_line(s)
        rewrite_with_ANSI(s)
    end

    return D
end
NEW_KEYBINDINGS = create_keybindings()

import Pkg
function insert_keybindings(repl = Base.active_repl)
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    main_mode = mirepl.interface.modes[1]
    p = mirepl.interface.modes[5]

    NEW_KEYBINDINGS["\e[A"] = (s,o...)-> begin
        LineEdit.edit_move_up(buffer(s)) || LineEdit.enter_prefix_search(s, p, true)
        Prompt.rewrite_with_ANSI(s)
    end
    # Down Arrow
    NEW_KEYBINDINGS["\e[B"] = (s,o...)-> begin
        LineEdit.edit_move_down(buffer(s)) || LineEdit.enter_prefix_search(s, p, false)
        Prompt.rewrite_with_ANSI(s)
    end

    main_mode.keymap_dict = LineEdit.keymap([NEW_KEYBINDINGS, main_mode.keymap_dict])
end

function _commit_line(s, data, c)
    move_input_end(s)
    rewrite_with_ANSI(s)
    println(terminal(s))
    add_history(s)
    state(s, mode(s)).ias = InputAreaState(0, 0)
end

function get_prompt(s)
    if isa(s, LineEdit.PromptState)
        prompt = s.p.prompt
    elseif isa(s, LineEdit.MIState)
        mode = s.current_mode
        if isa(mode, LineEdit.PrefixHistoryPrompt)
            prompt = mode.parent_prompt.prompt
        else
            prompt = mode.prompt
        end
    else
        error("Bug: $(typeof(s)) not accounted for")
    end
    isa(prompt, String) ? (return prompt) : (return prompt())
end

# Pasted from LineEdit.jl but the writes to the Terminal have been removed.
function refresh_multi_line(termbuf, terminal, buf, state, promptlength)
    cols = width(terminal)
    curs_row = -1 # relative to prompt (1-based)
    curs_pos = -1 # 1-based column position of the cursor
    cur_row = 0   # count of the number of rows
    buf_pos = position(buf)
    line_pos = buf_pos

    # Count the '\n' at the end of the line if the terminal emulator does (specific to DOS cmd prompt)
    miscountnl = @static Sys.iswindows() ? (isa(Terminals.pipe_reader(terminal), Base.TTY) && !Base.ispty(Terminals.pipe_reader(terminal))) : false
    lindent = promptlength
    indent = promptlength # TODO this gets the cursor right but not the text
    # Now go through the buffer line by line
    seek(buf, 0)
    moreinput = true # add a blank line if there is a trailing newline on the last line
    while moreinput
        l = readline(buf, keep=true)
        moreinput = endswith(l, "\n")
        # We need to deal with on-screen characters, so use textwidth to compute occupied columns
        llength = textwidth(l)
        slength = sizeof(l)
        cur_row += 1
        cmove_col(termbuf, lindent + 1)
        # We expect to be line after the last valid output line (due to
        # the '\n' at the end of the previous line)
        if curs_row == -1
            # in this case, we haven't yet written the cursor position
            line_pos -= slength # '\n' gets an extra pos
            if line_pos < 0 || !moreinput
                num_chars = (line_pos >= 0 ? llength : textwidth(l[1:prevind(l, line_pos + slength + 1)]))
                curs_row, curs_pos = divrem(lindent + num_chars - 1, cols)
                curs_row += cur_row
                curs_pos += 1
                # There's an issue if the cursor is after the very right end of the screen. In that case we need to
                # move the cursor to the next line, and emit a newline if needed
                if curs_pos == cols
                    # only emit the newline if the cursor is at the end of the line we're writing
                    if line_pos == 0
                        write(termbuf, '\n')
                        cur_row += 1
                    end
                    curs_row += 1
                    curs_pos = 0
                    cmove_col(termbuf, 1)
                end
            end
        end
        cur_row += div(max(lindent + llength + miscountnl - 1, 0), cols)
        lindent = indent
    end
    seek(buf, buf_pos)

    # Let's move the cursor to the right position
    # The line first
    n = cur_row - curs_row
    if n > 0
        cmove_up(termbuf, n)
    end

    #columns are 1 based
    cmove_col(termbuf, curs_pos + 1)

    # Updated cur_row,curs_row
    return InputAreaState(cur_row, curs_row)
end

end # module


text:
@testset "spgl1" begin

    @testset "Explicit A" begin

        DTs = (Float32, Float64, Complex{Float32}, Complex{Float64})

        for DTA in DTs, DTx in DTs
            @testset "SPGL1 - $DTA $DTx" begin
                A = rand(DTA, 100, 10)
                x = rand(DTx, 10)
                b = A*x
                x_out = spgl1(A,b, options = spgOptions(verbosity = 0))[1]
                @test snr(x, x_out) > 80
            end
        end
    end

end

ext: jl


text:
immutable LinearSpline <: AbstractInterpolation
    LinearSpline() = LinearSpline
end

immutable LinearSplineInterpolation{T,N}
    x::Vector{T}
    y::Vector{T}
    m::Vector{T}
end

function default_settings{D,T}(::Type{LinearSpline},axes::NTuple{D,T})
    N = length(axes[1])
    return Irregular{D,N}, Nearest
end

function LinearSpline(x::AbstractVector,y::AbstractVector)
  issorted(x) || throw(ArgumentError("x points must be in ascending order"))
  nx = length(x)
  x = convert(Array,x)
  y = convert(Array,y)
  m = zeros(nx)

  for i=1:nx-1
      m[i] = (y[i+1]-y[i])/(x[i+1] - x[i])
  end
  x, y, m = promote(x,y,m)
  return LinearSplineInterpolation{eltype(x),nx}(x,y,m)
end

@inline function preprocess(::Type{LinearSpline}, axes, f)
    return LinearSpline(axes...,f)
end

function apply_bc(::Type{LinearSpline},BC::Type{Nearest}, I::Symbol, x::Symbol)
    apply_nearest(I,1,x,1)
end

function apply_bc(::Type{LinearSpline},BC::Type{Reflect}, I::Symbol, x::Symbol)
    apply_reflect(I,1,x,1)
end

function apply_bc{V}(::Type{LinearSpline},BC::Type{Value{V}}, I::Symbol, x::Symbol)
    apply_value(BC, I, 1, x, 1)
end

function apply_bc(::Type{LinearSpline},BC::Type{None}, I::Symbol, x::Symbol)
    apply_none()
end

function apply_bc(::Type{LinearSpline},BC::Type{Error}, I::Symbol, x::Symbol)
    apply_error(I, 1, x, 1)
end

@inline function get_indices{S}(::Type{LinearSpline}, G::Type{Uniform{1,S}}, I, x)
    dx = I.x[2]-I.x[1]
    i1 = clamp(Int(floor((x[1] - I.x[1])/dx)) + 1, 1, S.parameters[1][1]-1)
    return i1
end

@inline function get_indices{N}(::Type{LinearSpline}, G::Type{Irregular{1,N}}, I, x)
    xr = searchsortedlast(I.x,x[1],Base.Order.Forward)
    i1 = clamp(xr, 1, N-1)
    return i1
end

@inline function get_value(::Type{LinearSpline}, I, x, ix)
    return I.y[ix] + I.m[ix]*(x[1] - I.x[ix])
end


text:
# Forward BLAS calls to generic implementation
#
using LinearAlgebra
import LinearAlgebra.BLAS

# generic implementations

@reroute LinearAlgebra.BLAS.dot(x,y) LinearAlgebra.dot(Any, Any)
@reroute LinearAlgebra.BLAS.axpy!(a, x, y) LinearAlgebra.axpy!(Any,
                                                      AbstractArray,
                                                      AbstractArray)

gengemv!(tA, α, A, x, β, y) = LinearAlgebra.generic_matvecmul!(y, tA, A, x, LinearAlgebra.MulAddMul(α, β))

@reroute LinearAlgebra.BLAS.gemv!(tA, α, A, x, β, y) gengemv!(Any, Any, Any, Any, Any, Any)


text:
include("../src/hypergraph-affinity-functions.jl")
using MAT
using Plots
M = matread("Total_Ns.mat")
Ng = M["Ng"]
Nf = M["Nf"]
Nw = M["Nw"]

N = Nf
r = 10

## Bootstrap
for k = 2:4
B = 100
for j = 1:3
if  j == 1
    M = Nf
elseif j == 2
    M = Nw
else
    M = Ng
end

H1, H2 = affinities_from_N(M)
alpha = 1/2
B1 = baselines(r,alpha)
B2 = baselines(r,1-alpha)
tR1 = H1./B1
tR2 = H2./B2

# Set of types of hyperedges
Nk = weights(M[k,1:k+1])

TR1 = tR1[k,1:k]
TR0 = tR2[k,1:k]

R0, R1, MR1, MR0, SR1, SR0, A0, A1, MA1, MA0, SA1, SA0 = Bootstrap_From_N(Nk,B,alpha)

# Plot
sc = 1
lw = 2
ms = 9
s1 = 375
s2 = 300
gfs = 18
tfs = 18
titlesize = 18
color1 = :green
color2 = :blue
title = ""
leg = false
if j == 3
    xlab = "Affinity type t"
else
    xlab = ""
end

if k == 2
    yl = [.2,1.7]
elseif k ==3
    yl = [.2,2]
elseif k == 4
    yl = [.1,2.5]
end

if k == 2
    ylab = "Affinity / Baseline"
else
    ylab = ""
end
p = Plots.plot([1; k],[1; 1],yscale = :identity,fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = 1:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs,ylim = yl)
plot!(p,[TR0 TR0], fillrange=[MR0-sc*SR0 MR0+sc*SR0], fillalpha=0.3,color = color2, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color2, label = class1)
plot!(p,1:k,[TR1 TR1], fillrange=[MR1-sc*SR1 MR1+sc*SR1],fillalpha=0.3,color = color1, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color1, label = class2)

savefig("Figures/Ratios_$(k)_$(j)_sc_$sc.pdf")
end
end

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
using QuantEcon
using PyPlot

include("jv.jl")

wp = JvWorker(grid_size=25)
v_init = wp.x_grid .* 0.5

f(x) = bellman_operator(wp, x)
V = compute_fixed_point(f, v_init, max_iter=40)

s_policy, phi_policy = bellman_operator(wp, V, return_policies=true)

# === plot policies === #
fig, ax = subplots()
ax[:set_xlim](0, maximum(wp.x_grid))
ax[:set_ylim](-0.1, 1.1)
ax[:plot](wp.x_grid, phi_policy, "b-", label="phi")
ax[:plot](wp.x_grid, s_policy, "g-", label="s")
ax[:set_xlabel]("x")
ax[:legend]()

max_stars_repo_name: gsulis/structural

max_issues_repo_name: gsulis/structural

max_forks_repo_name: gsulis/structural


text:
# ---
# title: 1163. Last Substring in Lexicographical Order
# id: problem1163
# author: Indigo
# date: 2021-03-06
# difficulty: Hard
# categories: String, Suffix Array
# link: <https://leetcode.com/problems/last-substring-in-lexicographical-order/description/>
# hidden: true
# ---
#
# Given a string `s`, return the last substring of `s` in lexicographical order.
#
#
#
# **Example 1:**
#
#
#
#     Input: "abab"
#     Output: "bab"
#     Explanation: The substrings are ["a", "ab", "aba", "abab", "b", "ba", "bab"]. The lexicographically maximum substring is "bab".
#
#
# **Example 2:**
#
#
#
#     Input: "leetcode"
#     Output: "tcode"
#
#
#
#
# **Note:**
#
#   1. `1 <= s.length <= 4 * 10^5`
#   2. s contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

function last_substring(s::String)
    m_ch = maximum(s)
    mc_pos = findall(==(m_ch), s)
    pos, window_len, new_window_len = first(mc_pos), 1, 1
    for idx in 2:length(mc_pos)
        if mc_pos[idx] == mc_pos[idx - 1] + 1
            new_window_len += 1
        else
            if new_window_len > window_len
                window_len = new_window_len
                pos = mc_pos[idx - 1] - new_window_len + 1
            end
            new_window_len = 1
        end
    end
    return s[pos:end]
end
## @lc code=end


text:
#        field symbol   name        is array
fields = [ (:X,       "positions",    true),
           (:P,       "momenta",      true),
           (:M,       "masses",       true),
           (:Z,       "numbers",      true),
           (:cell,    "cell",         false),
           (:calc,    "calculator",   false),
           (:cons,    "constraint",   false),
           (:pbc,     "pbc",          false),
           ]

for (S, name, isarray) in fields
   set_name = Meta.parse("set_$(name)!")
   get_name = Meta.parse("$name")
   if isarray
      @eval begin
         function $(get_name)(at::Atoms)
            return copy(at.$S)
         end
         function $(set_name)(at::Atoms, Q)
            if length(at.$S) != length(Q)
               at.$S = copy(Q)
            else
               at.$S .= Q
            end
            return at
         end
      end
   else
      function $(get_name)(at::Atoms)
         return at.$S
      end
      function $(set_name)(at::Atoms, Q)
         at.$S = Q
         return at
      end
   end
end

# # an FF preconditioner for pair potentials
# function precon!(tmp, V::PairPotential, r::T, R::JVec{T}, innerstab=T(0.1)
#                  ) where {T <: Number}
#    r = norm(R)
#    dV = evaluate_d!(tmp, V, r)
#    ddV = evaluate_dd!(tmp, V, r)
#    R̂ = R/r
#    return (1-innerstab) * (abs(ddV) * R̂ * R̂' + abs(dV / r) * (I - R̂ * R̂')) +
#              innerstab  * (abs(ddV) + abs(dV / r)) * I
# end

# struct MPairPotential <: MSitePotential
#    Z2V::Dict{Tuple{Int16,Int16}, WrappedPairPotential}
# end
#
# MPairPotential() =
#       MPairPotential(Dict{(Int16, Int16), WrappedPairPotential}())
#
# function MPairPotential(args...)
#    Vm = MPairPotential()
#    for p in args
#       @assert p isa Pair
#       z, V = _convert_multi_pp(sym::Symbol, p)
#       Vm.Z2V[z] = V
#    end
#    return Vm
# end
#
# F64fun(V::F64fun) = V
# _convert_multi_pp(sym::Tuple{Symbol, Symbol}, V) =
#       _convert_multi_pp(atomic_number.(sym), V)
# _convert_multi_pp(z::Tuple{<:Integer,<:Integer}, V) =
#       Int16.(z), V
#
# cutoff(V::MPairPotential) =
#    maximum( cutoff(V) for (key, V) in V.Z2V )


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    SymbolicNode

SymbolicNode is the basic building block of the symbolic graph in MXNet.jl.
It's a callable object and supports following calls:

    (s::SymbolicNode)(args::SymbolicNode...)
    (s::SymbolicNode)(; kwargs...)

Make a new node by composing `s` with `args`. Or the arguments
can be specified using keyword arguments.
"""
mutable struct SymbolicNode
  handle::MX_SymbolHandle
end

const SymbolicNodeOrReal = Union{SymbolicNode,Real}

Base.unsafe_convert(::Type{MX_handle}, s::SymbolicNode) =
  Base.unsafe_convert(MX_handle, s.handle)
Base.convert(T::Type{MX_handle}, s::SymbolicNode) = Base.unsafe_convert(T, s)
Base.cconvert(T::Type{MX_handle}, s::SymbolicNode) = Base.unsafe_convert(T, s)

"""
    deepcopy(s::SymbolicNode)

Make a deep copy of a SymbolicNode.
"""
function Base.deepcopy(s::SymbolicNode)
  ref_hdr = Ref{MX_handle}(C_NULL)
  @mxcall(:MXSymbolCopy, (MX_handle, Ref{MX_handle}), s, ref_hdr)
  SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    copy(s::SymbolicNode)

Make a copy of a SymbolicNode. The same as making a deep copy.
"""
Base.copy(s::SymbolicNode) = Base.deepcopy(s)

function (s::SymbolicNode)(args::SymbolicNode...)
  s = deepcopy(s)
  _compose!(s, args...)
end

function (s::SymbolicNode)(; kwargs...)
  s = deepcopy(s)
  _compose!(s; kwargs...)
end

"""
    Variable(name::Union{Symbol,AbstractString}; attrs)

Create a symbolic variable with the given name. This is typically used as a placeholder.
For example, the data node, acting as the starting point of a network architecture.

## Arguments

* `attrs::Dict{Symbol,<:AbstractString}`: The attributes associated with this `Variable`.
"""
function Variable(name::Union{Symbol,AbstractString}; attrs = Dict())
  attrs = convert(Dict{Symbol, AbstractString}, attrs)
  hdr_ref = Ref{MX_handle}(C_NULL)
  @mxcall(:MXSymbolCreateVariable, (char_p, Ref{MX_handle}), name, hdr_ref)
  node = SymbolicNode(MX_SymbolHandle(hdr_ref[]))
  for (k, v) in attrs
    set_attr(node, k, v)
  end
  node
end

"""
    @var <symbols>...

A handy macro for creating `mx.Variable`.

```julia
julia> x = @mx.var x
MXNet.mx.SymbolicNode x

julia> x, y, z = @mx.var x y z
(MXNet.mx.SymbolicNode x, MXNet.mx.SymbolicNode y, MXNet.mx.SymbolicNode z)
```
"""
macro var(n::Symbol)
  Expr(:call, :Variable, QuoteNode(n))
end

macro var(names::Symbol...)
  Expr(:tuple, map(n -> Expr(:call, :Variable, QuoteNode(n)), names)...)
end

"""
    Group(nodes::SymbolicNode...)

Create a `SymbolicNode` by grouping nodes together.
"""
function Group(nodes::SymbolicNode...)
  handles = MX_handle[nodes...]
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateGroup, (MX_uint, Ptr{MX_handle}, Ref{MX_handle}),
          length(handles), handles, ref_hdr)
  SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end


text:
abstract type ActivationFunction end

"""
    ReLU <: ActivationFunction

    (ReLU())(x) -> max.(x, 0)
"""
struct ReLU <: ActivationFunction end

"""
    Max <: ActivationFunction

    (Max())(x) -> max(maximum(x), 0)
"""
struct Max <: ActivationFunction end

"""
    Id <: ActivationFunction
Identity operator

    (Id())(x) -> x
"""
struct Id <: ActivationFunction end

"""
    GeneralAct <: ActivationFunction
Wrapper type for a general activation function.

### Usage
```julia
act = GeneralAct(tanh)

act(0) == tanh(0)           # true
act(10.0) == tanh(10.0)     # true
```
```julia
act = GeneralAct(x->tanh.(x))

julia> act(-2:2)
5-element Array{Float64,1}:
 -0.9640275800758169
 -0.7615941559557649
  0.0
  0.7615941559557649
  0.9640275800758169
```
"""
struct GeneralAct <: ActivationFunction
    f::Function
end

#=
TODO: consider writing our own interpolation scheme to avoid a dependency for this one thing.
Should only require a handful of functions.
Also NOTE: inherently not type stable unless parameterized.
=#
"""
    PiecewiseLinear <: ActivationFunction
Activation function that uses linear interpolation between supplied `knots`.
An extrapolation condition can be set for values outside the set of knots. Default is `Linear`.

    PiecewiseLinear(knots_x, knots_y, [extrapolation = Line()])

### Usage
```julia
kx = [0.0, 1.2, 1.7, 3.1]
ky = [0.0, 0.5, 1.0, 1.5]
act = PiecewiseLinear(kx, ky)

act(first(kx)) == first(ky) == 0.0
act(last(kx))  == last(ky)  == 1.5

act(1.0)    # 0.4166666666666667
act(-102)   # -42.5
```
```julia
act = PiecewiseLinear(kx, ky, Flat())

act(-102)   # 0.0
act(Inf)    # 1.5
```

### Extrapolations
- Flat()
- Line()
- constant (supply a number as the argument)
- Throw() (throws bounds error)

`PiecewiseLinear` uses [Interpolations.jl](http://juliamath.github.io/Interpolations.jl/latest/).
"""
struct PiecewiseLinear <: ActivationFunction
    f::Interpolations.Extrapolation
end

# default extrapolation is Line(). Can also do Flat() or supply a constant,
function PiecewiseLinear(knots_x::AbstractVector,
                         knots_y::AbstractVector,
                         extrapolation = Interpolations.Line())
    PiecewiseLinear(LinearInterpolation(knots_x, knots_y, extrapolation_bc = extrapolation))
end

# the type stable definitions probably don't need to go in the paper as-is
(f::ReLU)(x) = max.(x, zero(eltype(x)))
(f::Max)(x) = max(maximum(x), zero(eltype(x)))
(f::Id)(x) = x
(G::GeneralAct)(x) = G.f(x)
(PL::PiecewiseLinear)(x) = PL.f(x)


text:
module RocketSubjectInstanceTest

using Test
using Rocket

@testset "Subject" begin

    println("Testing: Subject")

    @testset begin
        subject1 = Subject(Int)
        @test eltype(subject1) === Int

        subject2 = Subject(Float64)
        @test eltype(subject2) === Float64
    end

    @testset begin
        subject = Subject(Int)

        actor1 = keep(Int)
        actor2 = keep(Int)

        subscription1 = subscribe!(subject, actor1)

        next!(subject, 0)
        next!(subject, 1)

        subscription2 = subscribe!(subject, actor2)

        next!(subject, 3)
        next!(subject, 4)

        unsubscribe!(subscription1)

        next!(subject, 5)
        next!(subject, 6)

        complete!(subject)

        unsubscribe!(subscription2)

        @test actor1.values == [ 0, 1, 3, 4 ]
        @test actor2.values == [ 3, 4, 5, 6 ]
    end

    @testset begin
        subject = Subject(Int)

        actor1 = keep(Int)
        actor2 = keep(Int)

        values = Int[]
        source = from(1:5) |> tap(d -> push!(values, d))

        subscription1 = subscribe!(subject, actor1)
        subscription2 = subscribe!(subject, actor2)

        subscribe!(source, subject)

        complete!(subject)

        @test values        == [ 1, 2, 3, 4, 5 ]
        @test actor1.values == [ 1, 2, 3, 4, 5 ]
        @test actor2.values == [ 1, 2, 3, 4, 5 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
    end

    @testset begin
        subject = Subject(Int)

        values      = []
        errors      = []
        completions = []

        actor = lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(errors, e),
            on_complete = ()  -> push!(completions, 0)
        )

        subscribe!(subject, actor)

        @test values      == [ ]
        @test errors      == [ ]
        @test completions == [ ]

        error!(subject, "err")

        @test values      == [ ]
        @test errors      == [ "err" ]
        @test completions == [ ]

        subscribe!(subject, actor)

        @test values      == [ ]
        @test errors      == [ "err", "err" ]
        @test completions == [ ]

    end

    @testset begin
        subject_factory = SubjectFactory(AsapScheduler())
        subject = create_subject(Int, subject_factory)

        actor1 = keep(Int)
        actor2 = keep(Int)

        values = Int[]
        source = from(1:5) |> tap(d -> push!(values, d))

        subscription1 = subscribe!(subject, actor1)
        subscription2 = subscribe!(subject, actor2)

        subscribe!(source, subject)

        complete!(subject)

        @test values        == [ 1, 2, 3, 4, 5 ]
        @test actor1.values == [ 1, 2, 3, 4, 5 ]
        @test actor2.values == [ 1, 2, 3, 4, 5 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
    end

    @testset begin
        subject1 = Subject(Int)
        subject2 = similar(subject1)

        @test subject1 !== subject2
        @test typeof(subject2) <: Subject
        @test eltype(subject2) === Int

        actor1 = keep(Int)
        actor2 = keep(Int)

        subscription1 = subscribe!(subject1, actor1)
        subscription2 = subscribe!(subject2, actor2)

        @test subscription1 !== subscription2

        next!(subject1, 1)

        @test actor1.values == [ 1 ]
        @test actor2.values == [  ]

        next!(subject2, 2)

        @test actor1.values == [ 1 ]
        @test actor2.values == [ 2 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
    end

    @testset begin
        subject = Subject(Int)

        limited = subject |> take(1)

        actor1 = keep(Int)
        actor2 = keep(Int)

        subscription1 = subscribe!(limited, actor1)
        subscription2 = subscribe!(subject, actor2)

        next!(subject, 1)

        @test actor1.values == [ 1 ]
        @test actor2.values == [ 1 ]

        next!(subject, 2)

        @test actor1.values == [ 1 ]
        @test actor2.values == [ 1, 2 ]
    end

    @testset begin

        main = Subject(Int)

        listener1 = Subject(Int)
        listener2 = Subject(Int)
        listener3 = Subject(Int)

        keep1 = keep(Int)
        keep2 = keep(Int)
        keep3 = keep(Int)

        k1 = subscribe!(listener1, keep1)
        k2 = subscribe!(listener2, keep2)
        k3 = subscribe!(listener3, keep3)

        sub1 = subscribe!(main, listener1)
        sub2 = subscribe!(main, listener2)
        sub3 = subscribe!(main, listener3)

        next!(main, 1)

        @test keep1.values == [ 1 ]
        @test keep2.values == [ 1 ]
        @test keep3.values == [ 1 ]

        unsubscribe!(sub2)

        next!(main, 2)

        @test keep1.values == [ 1, 2 ]
        @test keep2.values == [ 1 ]
        @test keep3.values == [ 1, 2 ]

        unsubscribe!(sub3)

        next!(main, 3)

        @test keep1.values == [ 1, 2, 3 ]
        @test keep2.values == [ 1 ]
        @test keep3.values == [ 1, 2 ]

        unsubscribe!(sub2) # sub2 here is intentional

        next!(main, 4)

        @test keep1.values == [ 1, 2, 3, 4 ]
        @test keep2.values == [ 1 ]
        @test keep3.values == [ 1, 2 ]
    end

end

end


text:
#
# Computes the norm of the projected gradient
# One method is defined for each possible scenario
# of bounds (all, only lower, only upper, none). This is
# needed such that the lower and upper bounds vectors can
# be of type Nothing if those bounds are not defined
#
# Both bounds
function pr_gradnorm(g, x, lower, upper)
    gnorm = zero(eltype(x))
    for i in eachindex(x)
        z = max(lower[i], min(upper[i], x[i] - g[i])) - x[i]
        gnorm = max(gnorm, abs(z))
    end
    return gnorm
end
# Only lower
function pr_gradnorm(g, x, lower, upper::Nothing)
    gnorm = zero(eltype(x))
    for i in eachindex(x)
        z = max(lower[i], x[i] - g[i]) - x[i]
        gnorm = max(gnorm, abs(z))
    end
    return gnorm
end
# Only upper
function pr_gradnorm(g, x, lower::Nothing, upper)
    gnorm = zero(eltype(x))
    for i in eachindex(x)
        z = min(upper[i], x[i] - g[i]) - x[i]
        gnorm = max(gnorm, abs(z))
    end
    return gnorm
end
# No bounds
pr_gradnorm(g, x, lower::Nothing, upper::Nothing) = maximum(abs, g)

#
# Compute a trial point
#
function compute_xn!(xn,x,t,g,lower,upper)
    for i in eachindex(x)
        xn[i] = x[i] - t * g[i]
        xn[i] = max(xn[i], lower[i])
        xn[i] = min(xn[i], upper[i])
    end
end
function compute_xn!(xn,x,t,g,lower,upper::Nothing)
    for i in eachindex(x)
        xn[i] = x[i] - t * g[i]
        xn[i] = max(xn[i], lower[i])
    end
end
function compute_xn!(xn,x,t,g,lower::Nothing,upper)
    for i in eachindex(x)
        xn[i] = x[i] - t * g[i]
        xn[i] = min(xn[i], upper[i])
    end
end
compute_xn!(xn,x,t,g,lower::Nothing,upper::Nothing) = @. xn = x - t * g


text:
function RecursiveArrayTools.recursive_unitless_bottom_eltype(
    ::Type{T},
) where {T<:ProductRepr}
    return recursive_unitless_bottom_eltype(T.parameters[1])
end
function RecursiveArrayTools.recursive_bottom_eltype(a::ProductRepr)
    return recursive_bottom_eltype(a.parts[1])
end
RecursiveArrayTools.recursive_unitless_eltype(::Type{T}) where {T<:ProductRepr} = T
RecursiveArrayTools.recursive_unitless_eltype(a::ProductRepr) = typeof(a)

function RecursiveArrayTools.recursivecopy!(x, y)
    map(copyto!, submanifold_components(x), submanifold_components(y))
    return x
end

function DiffEqBase.ODE_DEFAULT_NORM(u::ProductRepr, t)
    return sqrt(
        real(sum(p -> DiffEqBase.ODE_DEFAULT_NORM(p, t)^2, u.parts)) / length(u.parts),
    )
end


text:
mutable struct GFSimulateState
    trace::DynamicDSLTrace
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFSimulateState(gen_fn::GenerativeFunction, args::Tuple, params)
    trace = DynamicDSLTrace(gen_fn, args)
    GFSimulateState(trace, AddressVisitor(), params)
end

function traceat(state::GFSimulateState, dist::Distribution{T},
              args, key) where {T}
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    retval = random(dist, args...)

    # compute logpdf
    score = logpdf(dist, retval, args...)

    # add to the trace
    add_choice!(state.trace, key, retval, score)

    retval
end

function traceat(state::GFSimulateState, gen_fn::GenerativeFunction{T,U},
              args, key) where {T,U}
    local subtrace::U
    local retval::T

    # check key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # get subtrace
    subtrace = simulate(gen_fn, args)

    # add to the trace
    add_call!(state.trace, key, subtrace)

    # get return value
    retval = get_retval(subtrace)

    retval
end

function splice(state::GFSimulateState, gen_fn::DynamicDSLFunction,
                args::Tuple)
    prev_params = state.params
    state.params = gen_fn.params
    retval = exec(gen_fn, state, args)
    state.params = prev_params
    retval
end

function simulate(gen_fn::DynamicDSLFunction, args::Tuple)
    state = GFSimulateState(gen_fn, args, gen_fn.params)
    retval = exec(gen_fn, state, args)
    set_retval!(state.trace, retval)
    state.trace
end


text:
#=
Base:
- Julia version: 1.4.2
- Author: Bill Winnett
- Author email: bwinnett12@gmail.com
- Date: 2020-07-31
=#

#=
This algorithm is designed to pull out the names and sequences from Clustal files (.aln)
The good thing about these files is that they are already aligned but follows a different scheme

Scheme of .aln (Repeats until end of sequences):
Aligned Sequences   [1 -> n] - Line for match/mismatch identification - blank line
=#

# Parses the file for the name of each sequence
# Scheme runs several times but picks up scheme after first iteration
function parseSequenceNames(raw::Array)
    # Sets up an array to add names to
    SeqNameArr = String[]

    for i in eachindex(raw)
        # Finds a line of match signatures
        if startswith(raw[i], " ")
            # Try block is to catch the error if it went to the end of the document
            # (Which means its not an aln)
            try
                # Finds the empty line
                if raw[i+1] == ""

                    # The next couple of lines of the file should be the actual data
                    # The next conditionals filter out the non sequences
                    if raw[i+2] == ""
                        continue
                    # These two are probably redundant
                    elseif startswith(raw[i+2], " ")
                        continue
                    elseif occursin("*", raw[i+2])
                        continue

                    # Iterates through each line until its the match signature line
                    else
                        # Sets starting index at two away from the current index (i)
                        n = i + 2
                        # Blank character should be the end of the actual data
                        while !startswith(raw[n], " ")
                            # Adds each name to array
                            push!(SeqNameArr, split(raw[n], " ")[1])
                            n += 1
                        end

                        # Only one loop through is needed since scheme should run whole file
                        break
                    end
                end
            # This will be hit if it doesn't match the scheme. Thus not a clustal
            catch BoundsError
                @warn("Out of bounds. Probably not a clustal", BoundsError)
                return "Not an aln"
            end
        end
    end

    return SeqNameArr
end

# Parses the file for the sequences and puts them into an array
function parseSequences(raw::Array, names::Array)

SequencesArray = []
    # Iterates through each name
    for id in names
        SequenceById::String = ""

        # Runs loop through line of the clustal file once per name
        for i in eachindex(raw)

            # Appends each bit of the sequence to SequenceById
            if occursin(id, raw[i])
                SequenceById = string(SequenceById, split(raw[i], " ")[2])

            end
        end
        # Adds the sequence to an array. The index will be shared with the name
        push!(SequencesArray, SequenceById)
    end

    return SequencesArray
end

"""
    parseClustal(file_location)

- Call this function to convert the .aln or Clustal format to a tuple of (SeqID, Sequence)
- Julia version: 1.4.2

# Arguments

- `file_location::String`:

# Returns
- `Tuple(SequenceID, Sequence)`:
# Examples

```jldoctest
julia> parseClustal("aligned_ATP6.aln")
("staratp6.1", "--atatgaggaacttgga...")
("staratp6.2", "--atatgaggatctt_ga...")
("staratp6.3", "--atatgagcatcttgga...")
```
"""
function parseClustal(file_location::String)
    file_read = readlines(open(file_location))

    ArrayNames = parseSequenceNames(file_read)
    ArraySequences = parseSequences(file_read, ArrayNames)

    # Interweaves the name and sequence into a tuple
    interweavedArray = []
    for loc in eachindex(ArraySequences)
        push!(interweavedArray, (ArrayNames[loc], ArraySequences[loc]))
    end

    return interweavedArray
end


text:
############################################################################################
"""
$(SIGNATURES)
Recursively collect constraints of 'Param' struct. Not exported.

# Examples
```julia
```

"""
function _get_constraint(param::Param)
    return param.constraint
end
function _get_constraint(val::NamedTuple{names}) where {names}
    return NamedTuple{names}(Tuple(map(_get_constraint, val)))
end
"""
$(SIGNATURES)
Recursively collect values of 'Param' struct. Not exported.

# Examples
```julia
```

"""
function _get_val(param::Param)
    return param.val
end
function _get_val(val::NamedTuple{names}) where {names}
    return NamedTuple{names}(Tuple(map(_get_val, val)))
end

############################################################################################
"""
$(SIGNATURES)
Returns NamedTuple of true/false given parameter is not fixed. Not exported.

# Examples
```julia
```

"""
function _allparam(val::Bool)
    return val
end
function _allparam(val::A) where {A<:NamedTuple}
    return all(map(_allparam, val))
end
"""
$(SIGNATURES)
Returns NamedTuple of true/false given parameter is not fixed. Not exported.

# Examples
```julia
```

"""
_anyparam(val::Bool) = _allparam(val)
function _anyparam(val::A) where {A<:NamedTuple}
    return any(map(_anyparam, val))
end

############################################################################################
"""
$(SIGNATURES)
Return parameter names as a string in increasing order. Not exported.

# Examples
```julia
```

"""
function _paramnames(sym::Symbol, len::Integer)
    return [len == 1 ? string(sym) : string(sym, i) for i in Base.OneTo(len)]
end
function _paramnames(sym::NTuple{N,Symbol}, len::NTuple{N,Integer}) where {N}
    ArgCheck.@argcheck length(sym) == length(len)
    return [
        len[n] == 1 ? string(sym[n]) : string(sym[n], i)
        for n in Base.OneTo(length(sym)) for i in Base.OneTo(len[n])
    ]
end
"""
$(SIGNATURES)
Return all parameter names in increasing order. Not exported.

# Examples
```julia
```

"""
function paramnames(
    sym::Symbol,
    types::F,
    constraint,
    val
) where {F<:FlattenDefault}
    flatten, _ = construct_flatten(types, constraint, val)
    val_flattened = flatten(val)
    return _paramnames(sym, length(val_flattened))
end
function paramnames(
    sym::NTuple{N,Symbol},
    types::F,
    constraint,
    val,
) where {F<:FlattenDefault,N}
    flatten, unflatten = construct_flatten(types, constraint, val)
    return _paramnames(sym, unflatten._unflatten.lngth)
end
function paramnames(
    types::F,
    constraint::NamedTuple,
    val::NamedTuple{names},
) where {F<:FlattenDefault,N,names}
    return reduce(
        vcat, map(fld -> paramnames(fld, types, constraint[fld], val[fld]), names)
    )
end

############################################################################################
"""
$(SIGNATURES)
Count length of nested parameter tuple. Not exported.

# Examples
```julia
```

"""
function paramcount(sym::Symbol, types::F, val) where {F<:FlattenDefault}
    _flatten, _ = construct_flatten(types, val)
    val_flattened = _flatten(val)
    return length(val_flattened)
end
function paramcount(sym::NTuple{N,Symbol}, types::F, val) where {F<:FlattenDefault,N}
    _flatten, _unflatten = construct_flatten(types, val)
    return _unflatten._unflatten.lngth
end
function paramcount(types::F, val::NamedTuple{names}) where {F<:FlattenDefault,names}
    return map(fld -> paramcount(fld, types, val[fld]), names)
end

############################################################################################
# Export


text:
"""
Grammar

Represents a grammar and its production rules.
Use the @grammar macro to create a Grammar object.
"""
struct ContextFreeGrammar <: Grammar
	rules::Vector{Any}    # list of RHS of rules (subexpressions)
	types::Vector{Symbol} # list of LHS of rules (types, all symbols)
	isterminal::BitVector # whether rule i is terminal
	iseval::BitVector     # whether rule i is an eval rule
	bytype::Dict{Symbol,Vector{Int}}   # maps type to all rules of said type
	childtypes::Vector{Vector{Symbol}} # list of types of the children for each rule. Empty if terminal
end
"""
@grammar

Define a grammar and return it as a Grammar. For example:
```julia-repl
grammar = @grammar begin
R = x
R = 1 | 2
R = R + R
end
```
"""
macro cfgrammar(ex)
	rules = Any[]
	types = Symbol[]
	bytype = Dict{Symbol,Vector{Int}}()
	for e in ex.args
		if isa(e, Expr)
			if e.head == :(=)
				s = e.args[1] # name of return type
				rule = e.args[2] # expression?
				rvec = Any[]
				_parse_rule!(rvec, rule)
				for r in rvec
					push!(rules, r)
					push!(types, s)
					bytype[s] = push!(get(bytype, s, Int[]), length(rules))
				end
			end
		end
	end
	alltypes = collect(keys(bytype))
	is_terminal = [isterminal(rule, alltypes) for rule in rules]
	is_eval = [iseval(rule) for rule in rules]
	childtypes = [get_childtypes(rule, alltypes) for rule in rules]
	return ContextFreeGrammar(rules, types, is_terminal, is_eval, bytype, childtypes)
end

_parse_rule!(v::Vector{Any}, r) = push!(v, r)

function _parse_rule!(v::Vector{Any}, ex::Expr)
	if ex.head == :call && ex.args[1] == :|
		terms = length(ex.args) == 2 ?
		collect(interpret(ex.args[2])) :    #|(a:c) case
		ex.args[2:end]                      #a|b|c case
		for t in terms
			_parse_rule!(v, t)
		end
	else
		push!(v, ex)
	end
end

function Base.display(rulenode::RuleNode, grammar::ContextFreeGrammar)
	root = get_executable(rulenode, grammar)
	if isa(root, Expr)
	    walk_tree(root)
	else
	    root
	end
end

function _next_state!(node::RuleNode, grammar::Grammar, max_depth::Int)

	if max_depth < 1
	    return (node, false) # did not work
	elseif isterminal(grammar, node)
	    # do nothing
	    if iseval(grammar, node.ind) && (node._val === nothing)  # evaluate the rule
		node._val = eval(grammar.rules[node.ind].args[2])
	    end
	    return (node, false) # cannot change leaves
	else # !isterminal
	    # if node is not terminal and doesn't have children, expand every child
	    if isempty(node.children)
		if max_depth ≤ 1
		    return (node,false) # cannot expand
		end

		# build out the node
		for c in child_types(grammar, node)

		    worked = false
		    i = 0
		    child = RuleNode(0)
		    child_rules = grammar[c]
		    while !worked && i < length(child_rules)
			i += 1
			child = RuleNode(child_rules[i])

			if iseval(grammar, child.ind) # if rule needs to be evaluated (_())
			    child._val = eval(grammar.rules[child.ind].args[2])
			end
			worked = true
			if !isterminal(grammar, child)
			    child, worked = _next_state!(child, grammar, max_depth-1)
			end
		    end
		    if !worked
			return (node, false) # did not work
		    end
		    push!(node.children, child)
		end

		return (node, true)
	    else # not empty
		# make one change, starting with rightmost child
		worked = false
		child_index = length(node.children) + 1
		while !worked && child_index > 1
		    child_index -= 1
		    child = node.children[child_index]

		    # this modifies the node if succesfull
		    child, child_worked = _next_state!(child, grammar, max_depth-1)
		    while !child_worked
			child_type = return_type(grammar, child)
			child_rules = grammar[child_type]
			i = something(findfirst(isequal(child.ind), child_rules), 0)
			if i < length(child_rules)
			    child_worked = true
			    child = RuleNode(child_rules[i+1])

			    # node needs to be evaluated
			    if iseval(grammar, child.ind)
				child._val = eval(grammar.rules[child.ind].args[2])
			    end

			    if !isterminal(grammar, child)
				child, child_worked = _next_state!(child, grammar, max_depth-1)
			    end
			    node.children[child_index] = child
			else
			    break
			end
		    end

		    if child_worked
			worked = true

			# reset remaining children
			for child_index2 in child_index+1 : length(node.children)
			    c = child_types(grammar, node)[child_index2]
			    worked = false
			    i = 0
			    child = RuleNode(0)
			    child_rules = grammar[c]
			    while !worked && i < length(child_rules)
				i += 1
				child = RuleNode(child_rules[i])

				if iseval(grammar, child.ind)
				    child._val = eval(grammar.rules[child.ind].args[2])
				end

				worked = true
				if !isterminal(grammar, child)
				    child, worked = _next_state!(child, grammar, max_depth-1)
				end
			    end
			    if !worked
				break
			    end
			    node.children[child_index2] = child
			end
		    end
		end

		return (node, worked)
	    end
	end
    end


text:
"""
mutable struct InverseDistanceSurrogate{X,Y,P,L,U} <: AbstractSurrogate

The square polynomial model can be expressed by 𝐲 = 𝐗β + ϵ, with β = 𝐗ᵗ𝐗⁻¹𝐗ᵗ𝐲
"""
mutable struct SecondOrderPolynomialSurrogate{X,Y,B,L,U} <: AbstractSurrogate
    x::X
    y::Y
    β::B
    lb::L
    ub::U
end

function SecondOrderPolynomialSurrogate(x,y,lb,ub)
    X = _construct_2nd_order_interp_matrix(x, first(x))
    Y = _construct_y_matrix(y, first(y))
    β = X\Y
    return SecondOrderPolynomialSurrogate(x, y, β, lb, ub)
end

function _construct_2nd_order_interp_matrix(x, x_el)
    n = length(x)
    d = length(x_el)
    D = 1 + 2*d + d*(d-1)÷2
    X = ones(eltype(x_el), n, D)
    for i = 1:n, j = 1:d
        X[i, j+1] = x[i][j]
    end
    idx = d + 1
    for j = 1:d, k = j+1:d
        idx += 1
        for i = 1:n
            X[i, idx] = x[i][j]*x[i][k]
        end
    end
    for i = 1:n, j = 1:d
        X[i, j+1+d+d*(d-1)÷2] = x[i][j]^2
    end
    return X
end

_construct_y_matrix(y, y_el::Number) = y
_construct_y_matrix(y, y_el) = [y[i][j] for i=1:length(y), j=1:length(y_el)]

function (my_second_ord::SecondOrderPolynomialSurrogate)(val)
    #just create the val vector as X and multiply
    d = length(val)

    y = my_second_ord.β[1, :]
    for j = 1:d
        y += val[j]*my_second_ord.β[j+1, :]
    end
    idx = d + 1
    for j = 1:d, k = j+1:d
        idx += 1
        y += val[j] * val[k] * my_second_ord.β[idx, :]
    end
    for j = 1:d
        y += val[j]^2 * my_second_ord.β[j+1+d+d*(d-1)÷2, :]
    end
    return _match_container(y, first(my_second_ord.y))
end

function add_point!(my_second::SecondOrderPolynomialSurrogate, x_new, y_new)
    if eltype(x_new) == eltype(my_second.x)
        append!(my_second.x, x_new)
        append!(my_second.y, y_new)
    else
        push!(my_second.x, x_new)
        push!(my_second.y, y_new)
    end
    X = _construct_2nd_order_interp_matrix(my_second.x, first(my_second.x))
    Y = _construct_y_matrix(my_second.y, first(my_second.y))
    β = X\Y
    my_second.β = β
    nothing
end


text:
@with_kw struct LQR
    γ::Float64 = 1.0
end

discount(mdp::LQR) = mdp.γ

generate_start_state(mdp::LQR) = rand(Normal(0.3,0.1))

function transition(mdp::LQR, s::Float64, a::Float64)
    # NOTE: Truncated to prevent going off to infinity with poor policies
    return truncated(Normal(s + a, 0.1), -10.0, 10.0)
end
reward(mdp::LQR, s::Float64, a::Float64) = -s^2

function MDP(mdp::LQR; γ::Float64=mdp.γ)
    return MDP(
            γ,
            nothing, # continuous state
            nothing, # continuous action space
            (s,a, s′=nothing) -> begin
                S′ = transition(mdp, s, a)
                if s′ == nothing
                    return S′
                end
                return pdf(S′, s′)
            end,
            (s,a) -> reward(mdp, s, a),
            (s, a)->begin
                s′ = rand(transition(mdp,s,a))
                r = reward(mdp, s, a)
                return (s′, r)
            end
        )
end

const LINEAR_1D_GAUSSIAN_POLICY = (θ::Vector{Float64}, s::Float64) -> begin
            λ = θ[1]
            σ = abs(θ[2]) + 0.00001
            return Normal(λ*s, σ)
        end
const LINEAR_1D_GAUSSIAN_POLICY_∇ = (θ::Vector{Float64}, a::Float64, s::Float64) -> begin
            λ = θ[1]
            σ = θ[2] # note: standard dev would be abs(θ[2])
            ∇₁ = s*(a - λ*s)*exp(-(a-λ*s)^2/(2*σ^2)) / (sqrt(2π)*abs(σ)^3)
            ∇₂ = sign(σ)*exp(-(a-λ*s)^2/(2*σ^2)) * ((a-λ*s)^2 - σ^2)/(sqrt(2π)*σ^4)
            if abs(σ) < 1e-6
                ∇₁ = 0.0
                ∇₂ = 0.0
            end
            return [∇₁, ∇₂]
        end
const LINEAR_1D_GAUSSIAN_POLICY_∇logπ = (θ::Vector{Float64}, a::Float64, s::Float64) -> begin
            λ = θ[1]
            σ = θ[2]
            ∇₁ = (a*s - λ*s^2) / (σ^2)
            ∇₂ = (a^2-2a*λ*s - σ^2 + λ^2*s^2)/(σ^3)
            if abs(σ) < 1e-6
                ∇₁ = 0.0
                ∇₂ = 0.0
            end
            return [∇₁, ∇₂]
        end
const LINEAR_1D_GAUSSIAN_POLICY_KLDIV = (θ::Vector{Float64}, θ′::Vector{Float64}, s::Float64) -> begin
            μ₁ = θ[1]*s
            μ₂ = θ′[1]*s
            σ₁ = abs(θ[2])
            σ₂ = abs(θ′[2])
            return log(σ₂/σ₁) + (σ₁^2 + (μ₁ - μ₂)^2)/(2*σ₂^2) - 1/2
        end
const LINEAR_1D_GAUSSIAN_POLICY_H_KL = (θ::Vector{Float64}, θ′::Vector{Float64}, s::Float64) -> begin
            F = zeros(Float64,2,2)
            λ₁ = θ[1]
            λ₂ = θ′[1]
            σ₁ = θ[2]
            σ₂ = θ′[2]
            F[1,1] = s^2/(σ₂^2)
            F[2,1] = F[1,2] = 2*s^2*(λ₁ - λ₂)/(σ₂^3)
            F[2,2] = 3*((λ₁ - λ₂)^2*s^2+σ₁^2)/(σ₂^4) - 1/(σ₂)^2
            return F
        end


text:
import Pkg; Pkg.activate(".")

using Mangal
using DataFrames
using CSV

total_networks = count(MangalNetwork)
count_per_page = 100
number_of_pages = convert(Int, ceil(total_networks/count_per_page))

LS = DataFrame(fill(Int64, 5), [:id, :S, :L, :P, :H], total_networks)

global cursor = 1
@progress "Paging networks" for page in 1:number_of_pages
    global cursor
    networks_in_page = Mangal.networks("count" => count_per_page, "page" => page-1)
    @progress "Counting items" for current_network in networks_in_page
        S = count(MangalNode, current_network)
        L = count(MangalInteraction, current_network)
        P = count(MangalInteraction, current_network, "type" => "predation")
        H = count(MangalInteraction, current_network, "type" => "herbivory")
        LS[cursor,:] .= (current_network.id, S, L, P, H)
        cursor = cursor + 1
    end
end

CSV.write("ls.csv", LS)


text:
using AnalogDiscovery

function go(fs::Float64)
  cDev = enumDevices()
  @show ("Available",cDev)
  myDev = Int32(0)
  if cDev > 0
    if enumDeviceIsOpened(myDev)
      error("Error opening device, already in use.")
    end
    hdwf = deviceOpen(myDev)
  else
    error("No devices available.")
  end

  analogInFrequencySet(hdwf,fs)
  bmin,bmax = analogInBufferSizeInfo(hdwf)
  bsize = min(Int32(round(fs)),bmax)
  analogInBufferSizeSet(hdwf,bsize)
  analogInAcquisitionModeSet(hdwf,ACQMODESINGLE)

  chVoltMin,chVoltMax,chSteps = analogInChannelRangeInfo(hdwf)
  chCount = analogInChannelCount(hdwf)
  for ch in Int32(0):Int32(chCount-1)
    analogInChannelEnableSet(hdwf,ch,true)
    analogInChannelFilterSet(hdwf,ch,FILTERAVERAGE)
    analogInChannelRangeSet(hdwf,ch,chVoltMax)
  end

  # Collect 10-sec worth of data from all channels
  blks = Int(ceil(10*fs/bsize))
  chData = [Vector{Float64}(blks*bsize) for i in 1:chCount]

  for b in 1:blks
    analogInConfigure(hdwf,false,true) # start a single acquistion
    while analogInStatus(hdwf,true) != DWFSTATEDONE
      sleep(0.001)
    end

    chDataBuf = Vector{Float64}(bsize)
    for ch in Int32(0):Int32(chCount-1)
      analogInStatusData!(chDataBuf,hdwf,ch)
      chData[ch+1][(b-1)*bsize+(1:bsize)] = chDataBuf[:]
      #showall(chDataBuf)
    end
  end

  deviceClose(hdwf)
  return chData
end

const fs = 128.0
chData = go(fs)


text:
using FHIRClient
using Test

import Dates
import HTTP
import JSON3
import Test

@testset "FHIRClient.jl" begin
    include("unit.jl")
    include("integration.jl")
end


text:
# Arvind Prasadan
# 2019 January
# DMF Package Code
# Optshrink Function

"""
	optshrink(Y, r)

Perform a rank-r denoising of a data matrix Y via the OptShrink algorithm, described in http://doi.org/10.1109/TIT.2014.2311661 (R R Nadakuditi, 2014)

# Arguments
- `Y`: 2D Array (Data Matrix), where Y = X + Noise, and we want to find X
- `r`: An estimate of the rank of X (Integral)

# Output
- `Xh`: An estimate of X
- `w`: The shrunken singular values; if r = rank(Y), w is all zero
"""
function optshrink(Y, r)

    r = round(abs(r[1])) # Sanity check rank
    r = (1 <= r) ? r : 1

    (m, n) = size(Y)
    r = minimum([r, m, n]) # ensure r <= min(m, n)
    if r == m || r == n
	return Y, zeros(r)
    end

    # U is [m, min(m, n)], s is min(m, n), V is [n, min(n, m)]
    (U, s, V) = svd(Y, full = false)

    sv = s[(r + 1):end] # tail singular values for noise estimation

    w = zeros(r)
    for k = 1:1:r
        (D, Dder) = D_transform_from_vector(s[k], sv, max(m, n) - r, min(m, n) - r)
        w[k] = -2.0 * D / Dder
    end

    w[isnan.(w)] .= 0.0

    Xh = U[:, 1:r] * Diagonal(w) * V[:, 1:r]'

    return Xh, w
end

"""
	 D_transform_from_vector(z, sn, m, n)

Find the D Transform of the singular value spectrum sn at a point z, for an m x n matrix.

# Arguments
- `z`: Point to evaluate transform
- `sn`: Vector of singular values
- `m`: Matrix dimension
- `n`: Matrix dimension

# Outputs
- `D`: D Transform
- `D_der`: D Transform Derivative
"""
function D_transform_from_vector(z, sn, m, n)
    # sn is of length n <= m

    sm = [sn; zeros(m - n)] # m x 1

    inv_n = 1.0 ./ (z^2.0 .- sn.^2.0) # vector corresponding to diagonal
    inv_m = 1.0 ./ (z^2.0 .- sm.^2.0)

    D1 = (1.0 / n) * sum(z * inv_n)
    D2 = (1.0 / m) * sum(z * inv_m)

    D = D1 * D2 # eq (16a) in paper

    # derivative of D transform
    D1_der = (1.0 / n) * sum(-2.0 * z^2.0 .* inv_n.^2.0 + inv_n)
    D2_der = (1.0 / m) * sum(-2.0 * z^2.0 .* inv_m.^2.0 + inv_m)

    D_der = D1 * D2_der + D2 * D1_der # eq (16b) in paper

    return (D, D_der)
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Random
using ArraysOfArrays
import TypedTables

@testset "scalar_shape" begin
    @test @inferred(ValueShapes._valshapeoftype(Int)) == ScalarShape{Int}()
    @test @inferred(ValueShapes._valshapeoftype(Complex{Float64})) == ScalarShape{Complex{Float64}}()

    @test @inferred(valshape(3)) == ScalarShape{Int}()

    @test @inferred(size(ScalarShape{Real}())) == ()
    @test @inferred(size(ScalarShape{Complex}())) == ()

    @test @inferred(size(ScalarShape{Real}())) == ()
    @test @inferred(size(ScalarShape{Complex}())) == ()

    @test @inferred(eltype(ScalarShape{Real}())) == Real
    @test @inferred(ValueShapes.nonabstract_eltype(ScalarShape{Real}())) == Float64

    @test @inferred(totalndof(ScalarShape{Int}())) == 1
    @test @inferred(totalndof(ScalarShape{Complex{Float64}}())) == 2

    @test @inferred(ScalarShape{Real}()(undef)) === zero(Float64)
    @test @inferred(ScalarShape{Complex}()(undef)) === zero(Complex{Float64})

    @test @inferred(ScalarShape{Real}()([42])) == 42

    @test begin
        shape = ScalarShape{Real}()
        valshape(shape.(push!(@inferred(VectorOfSimilarVectors(shape)), @inferred(Vector(undef, shape))))[1]) == valshape(shape(undef))
    end
end


text:
# include("julia/yinsGraph/graphAlgs.jl")
# include("julia/yinsGraph/graphGenerators.jl")
# include("julia/yinsGraph/treeAlgs.jl")

using yinsGraph

#partition into clusters of radii O(x(n)lnn) where 1/x(n) edges border the cluster
  # means radius <= ln(n)/#edges on border
#construct shortest path spanning tree in each cluster
  #what vertex is the base for the shortest spanning tree??
#add edges of spanning tree to output tree T
#recurse after collapsing each cluster into one vertex

# QUESTIONS:

#store cluster as array..?
#yes

# store new combined graph as a new graph though.
#better to just copy onto a new sparseMatrix..? Then we have to keep track of all of those matrices..?

#show proposal

#what's the start for the shortest path tree??

#AFTER TALKING WITH SPIELMAN:
# imagine unweighted graph!!!!!
  # start with some node
  # breathfirst search on neighbors of nodes (n, n2, n3 each is group of neighbors)
  # go until
  # number of nodes in cluster / number of neighbors of cluster is small (<= ln(n)/r).
      # this ratio is a function of R somehow
      #their x(n) is
      #grow ball until nodes in cluster / number of neighbors <= ln(n)/r

      # |
      # V
  #start with some arbitrary node
  #volume = sum of all the degrees of nodes in cluster
  #boundary = # edges leaving cluster
  #grow balls (using basically breadth-first search -- gives us shortest path) until "boundary" of the ball/volume <= 1/x(n).
    # use their function of x(n)? make it parameter of code to be fixed later..?
  #don't really need radius function

  # in unweighted, two big issued:
    # 1. what is X(n)?
    # 2. how to we collapse..?
    # don't do data structures and stuff yet..?

    #to compute

#returns volume = sum of all the degrees of nodes in cluster
function getVolume(mat::SparseMatrixCSC, cluster)
  count = 0
  for v in 1:length(cluster)
    if cluster[v] == 1
      for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
        count += 1
      end #for
    end #if
  end #for
  return count
end #getVolume

#returns boundary = # edges leaving cluster
function getBoundary(mat::SparseMatrixCSC, cluster)
  count = 0
  for v in 1:length(cluster)
    if cluster[v] == 1
      for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
        if cluster[mat.rowval[ind]] == 0
          count += 1
        end #if
      end #for
    end #if
  end #for
  return count
end #getBoundary

# takes a matrix from which partition will be constructed (referred to as the parent matrix)
# a partition is a double array of nodes. each sub-array c (cluster) has n vertices (the n of the parent mat)
# and either a 1 or 0 indicating whether or not v is in cluster c
# this function also returns an array of edge indices (of the parent mat) with the edges in the tree added
function partitionMatrix(mat::SparseMatrixCSC)
  nVertices = mat.n
  partition = Array{Float64,1}[]
  nClusters = 0
  visited = zeros(Bool, nVertices)
  nAddedToClusters = 0
  treeInds = zeros(Int64, nVertices)

  # build all the clusters
  while (nAddedToClusters != nVertices)
    nClusters += 1
    push!(partition, zeros(Bool, nVertices))

    #select arbitrary unvisited node
    start = 1
    for node in 1:nVertices
      if !visited[node]
        start = node
        break
      end #if
    end #for
    visited[start] = true

    partition[nClusters][start] = true

    nNodesInCluster = 1
    nAddedToClusters += 1

    lastnNodesInCluster = 0

    #building a single cluster
    while !(nAddedToClusters == nVertices)
      volume = getVolume(mat, partition[nClusters])
      boundary = getBoundary(mat, partition[nClusters])

      if nNodesInCluster > 2
        # x = 0.7199*exp(sqrt(log(nNodesInCluster) * log(log(nNodesInCluster))))
        x = .75*exp(sqrt(log(nNodesInCluster) * log(log(nNodesInCluster))))
        # x = log(nNodesInCluster)
      else
        x = 10
      end

      if (convert(Float64, volume) / convert(Float64, boundary) >= x) || lastnNodesInCluster == nNodesInCluster
        break
      end
      # else
      lastnNodesInCluster = nNodesInCluster

      newVertices = zeros(Bool, nVertices)
      for v in 1:length(partition[nClusters])
        if partition[nClusters][v] == true
          for eInd in mat.colptr[v]:(mat.colptr[v+1]-1) #eInd is the edge
            otherV = mat.rowval[eInd] #this gives us the other vertex

            #doing this step before for boundary... maybe can combine?
            if !visited[otherV]
              newVertices[otherV] = true

              treeInds[mat.rowval[eInd]] = eInd #..?
            end #if
          end #for
        end #if
      end #for

      for v in 1:length(newVertices)
        if newVertices[v] == 1
          partition[nClusters][v] = true
          nAddedToClusters += 1
          nNodesInCluster += 1
          visited[v] = true
        end #if
      end #for
      # end #if/else

      # lastnNodesInCluster = nNodesInCluster
    end #while
  end #while

  return partition, treeInds

end #partitionMatrix

function collapsePartition(mat::SparseMatrixCSC, partition, map)
  nVertices = mat.n
  nClusters = length(partition)
  iVector = Int64[]
  jVector = Int64[]
  vVector = Int64[]
  for c in 1:nClusters
    for v in 1:nVertices #all vertices in cluster i
      if partition[c][v] == 1
        for eInd in mat.colptr[v]:(mat.colptr[v+1]-1) #all edges connecting to vertex j
          otherV = mat.rowval[eInd]

          #if the two ends of the edge are in the same cluster don't add edge
          if partition[c][otherV] == 1
            continue

          #otherwise... add edge to matrix
          else
            push!(iVector, c)
            push!(jVector, map[otherV])
            push!(vVector, 1)

            # push!(iVector, map[otherV])
            # push!(jVector, c)
            # push!(vVector, 1)
          end #if
        end
      end #if
    end #for
  end #for

  newMat = sparse(iVector, jVector, vVector, nClusters, nClusters, min)

  return newMat
end #collapsePartition

#start with some arbitrary node
  #volume = sum of all the degrees of nodes in cluster
  #boundary = # edges leaving cluster
  #grow balls (using basically breadth-first search -- gives us shortest path) until "boundary" of the ball/volume <= 1/x(n).
    # use their function of x(n)? make it parameter of code to be fixed later..?
    # i'm converting it to volume/boundary >= x(n)
  #don't really need radius function

  # in unweighted, two big issued:
    # 1. what is X(n)?
    # 2. how to we collapse..?
    # don't do data structures and stuff yet..?

  # To collapse, will build a while new matrix, and keep a mapping to the original matrix
    # mapping: n-degree array. each value is the cluster index (v's new compressed value)

function akpw_unweighted(mat::SparseMatrixCSC)
  nVertices = mat.n
  nEdges = nnz(mat)
  treeInds = zeros(Int64, nEdges)

  newMat = mat
  bigMapC = zeros(Int64, nVertices) #assuming map does C->D, bigMapC does A->C
  bigMapD = zeros(Int64, nVertices) #assuming map does C->D, bigMapD does A->D. don't think I need this
  map = zeros(Int64, nVertices)

  first = true

  while (true)
    partition, partitionTreeInds = partitionMatrix(newMat)
    println(length(partition))

    if !first
      for i = 1:length(bigMapD)
        bigMapC[i] = bigMapD[i]
      end #for
    end #if

    #making the new map
    for c = 1:length(partition)
      for v = 1:length(partition[c])
        # println("checking: ", c, ", ", v)
        if partition[c][v] == 1
          # println("printing: ", v)
          map[v] = c
        end #if
      end #for
    end #for

    if first
      for i = 1:length(bigMapD)
        bigMapD[i] = map[i]
        bigMapC[i] = map[i]
      end #for

    else
      for i = 1:length(bigMapD)
        bigMapD[i] = map[bigMapD[i]]
      end #for
    end #if

    #CAN WE SPEED THIS UP BY USING partition INSTEAD OF iterating through the maps..?
    # append every edge in partitionTree to treeInds
    for e in 1:length(partitionTreeInds)
      if partitionTreeInds[e] == 0
        continue
      end

      if first
        treeInds[partitionTreeInds[e]] = 1
      else

        a = findn(newMat)
        newr = a[1][partitionTreeInds[e]]
        newc = a[2][partitionTreeInds[e]]

        for i = 1:length(bigMapC)
          found = false

          #this finds an edge in big mat corresponding to e in the partitionTree
          if bigMapC[i] == newr
            for eInd in mat.colptr[i]:(mat.colptr[i+1]-1)
              edgeV = mat.rowval[eInd]
              if bigMapC[edgeV] == newc
                found = true
                # println("treeInd: ", eInd)
                treeInds[eInd] = 1
                break
              end #if
            end #for
          end #if

          if found
            break
          end #if
        end #for
      end #if/else

    end #for

    if (length(partition) == 1)
      break
    end #if

    newMat = collapsePartition(newMat, partition, map)

    first = false

  end #while

  t2 = zeros(Int64, nVertices - 1)
  t2Counter = 1

  for i in 1:length(treeInds)
    if treeInds[i] == 1
      t2[t2Counter] = i
      t2Counter += 1
    end
  end

  (ai,aj,av) = findnz(mat);
  tr2 = sparse(ai[t2],aj[t2],av[t2],nVertices,nVertices)
  tr2 = tr2 + tr2';

  return tr2

end #akpw

# k = 30
# ani = productGraph(pathGraph(k),pathGraph(k))
# akpw_unweighted(ani)

# a = grid2(10)
# n = size(a)[1]
# (ai,aj,av) = findnz(triu(a))

# @time stretchTree = akpw(a)
# primTree = prim(a)
# spy(primTree)
# println("done")
# println(stretchTree)

# println(compStretches(stretchTree, a))
# println(compStretches(stretchTree, a))

# BUGS:

# why do clusters have weird edge weights


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
    using ViewerGL
    GL = ViewerGL
end

function viewExplode(model; complete = false)
    V,CVs,FVs,EVs = Lar.pols2tria(model.G, model.T[1], model.T[2], model.T[3])

    GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1));
    GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
    GL.VIEW(GL.GLExplode(V,CVs[2 - complete:end],5,5,5,99,0.5));
end

function displayModel(model, exp = 1.5)
    V = model.G
    EV = Lar.cop2lar(model.T[1])
    FE = Lar.cop2lar(model.T[2])
    CF = Lar.cop2lar(model.T[3])
    FV = Lar.cop2lar(map(x -> Int8(x/2), abs.(model.T[2]) * abs.(model.T[1])))
    CV = []
    triangulated_faces = Lar.triangulate(convert(Lar.Points, V'), [model.T[1], model.T[2]])
    FVs = convert(Array{Lar.Cells}, triangulated_faces)
    EVs = Lar.FV2EVs(model.T[1], model.T[2])

    GL.VIEW([
        GL.GLPol(V,EV, GL.COLORS[1])
        GL.GLFrame
    ]);
    GL.VIEW(GL.GLExplode(V,FVs,exp,exp,exp,99));
    GL.VIEW(GL.GLExplode(V,EVs,exp,exp,exp,99,1));
end

function viewNum(model, numsize = 0.25)
	V = model.G
	EV = Lar.cop2lar(model.T[1])
	FV = Lar.cop2lar(map(x -> floor(Int8, x / 2), abs.(model.T[2]) * abs.(model.T[1])))
	VV = [[k] for k=1:size(V,2)]
	GL.VIEW([
        GL.numbering(numsize)((V,[VV, EV, FV]))
        GL.GLFrame
    ])
end

function viewNumEdge(model::CAGD.Model, τ::Int, numsize=0.25)
    V  = model.G
    EV = Lar.cop2lar(model.T[1])
    VV = [[k] for k=1:size(V,2)]
    faces = model.T[2][:, τ].nzind
    EVidx = unique(vcat([model.T[2][f, :].nzind for f in faces]...))
    EVs = Lar.cop2lar(model.T[1][EVidx, :])
    FVs = Lar.cop2lar(map(x -> floor(Int8, x / 2), abs.(model.T[2][faces, :]) * abs.(model.T[1])))
    VV = [[k] for k in unique(vcat(EVs...))]
    GL.VIEW([
        GL.GLPol(V, EV, GL.COLORS[1])
        GL.numbering(numsize)((V, [VV, EVs, FVs]))
        GL.GLFrame
    ])
end

function viewNumFace(model::CAGD.Model, τ::Int, numsize=0.25)
    V  = model.G
    EV = Lar.cop2lar(model.T[1])
    VV = [[k] for k=1:size(V,2)]
    EVidx = model.T[2][τ, :].nzind
    EVs = Lar.cop2lar(model.T[1][EVidx, :])
    FVs = Lar.cop2lar(map(x -> floor(Int8, x / 2), abs.(model.T[2][τ:τ, :]) * abs.(model.T[1])))
    VV = [[k] for k in unique(vcat(EVs...))]
    GL.VIEW([
        GL.GLPol(V, EV, GL.COLORS[1])
        GL.numbering(numsize)((V, [VV, EVs, FVs]))
        GL.GLFrame
    ])
end


text:
using MathOptInterface, COSMO, Test, LinearAlgebra
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges
const MOIU = MOI.Utilities
MOIU.@model(COSMOModelData,
        (),
        (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
        (MOI.Zeros, MOI.Nonnegatives, MOI.SecondOrderCone,
         MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone, MOI.DualExponentialCone),
        (MOI.PowerCone, MOI.DualPowerCone),
        (),
        (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
        (MOI.VectorOfVariables,),
        (MOI.VectorAffineFunction,),);

struct UnsupportedModelAttribute  <: MOI.AbstractModelAttribute end

get_inner_optimizer(bridged) = bridged.model.optimizer

@testset "MOI Wrapper" begin

    @testset "Basic properties" begin
        optimizer =  COSMO.Optimizer(check_termination = 1, verbose = false);
        @test sprint(show, optimizer) == string("Empty COSMO - Optimizer")
        @test MOI.supports(optimizer, MOI.ObjectiveFunction{MOI.SingleVariable}())
        @test MOI.supports(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test MOI.supports(optimizer, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
        @test MOI.supports(optimizer, MOI.ObjectiveSense())
        @test MOI.supports(optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)
        @test MOI.supports(optimizer, MOI.ConstraintPrimalStart(), MOI.ConstraintIndex)
        @test MOI.supports(optimizer, MOI.ConstraintDualStart(), MOI.ConstraintIndex)

    end
    # Solve the following problem:
    # min c'*x
    # vec(A1)' * x == b1
    # vec(A2)' * x == b2
    # mat(x) is posdeftriangle
    # with:
    A1_t = [1.0; 0; 3; 2; 14; 5];
    A2_t = [0.0; 4; 6; 16; 0; 4];
    C_t = [1.; 4; 9; 6; 0; 7];

    A1 = [1.0 0 1; 0 3 7; 1 7 5];
    A2 = [0.0 2 8; 2 6 0; 8 0 4];
    C = [1.0 2 3; 2 9 0; 3 0 7];
    b1 = 11.0;
    b2 = 19.0;

    cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
    optimizer =  COSMO.Optimizer(check_termination = 1, verbose = false);
    cached = MOIU.CachingOptimizer(cache, optimizer)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)

    x = MOI.add_variables(bridged, 6);
    # define objective function:
    objectiveFunction = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(vec(C_t),x[1:6]),0.0);
    MOI.set(bridged, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),objectiveFunction);
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE);

    # eq constraints
    con1 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A1_t,x[1:6]),0.0), MOI.EqualTo(b1));
    con2 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A2_t,x[1:6]),0.0), MOI.EqualTo(b2));

    # SDP constraint
    con3 = MOI.add_constraint(bridged, MOI.VectorOfVariables([x[1], x[2], x[3], x[4], x[5], x[6]]), MOI.PositiveSemidefiniteConeTriangle(3));

    # copy model into optimizer and solve
    @test MOI.get(bridged, MOI.ListOfVariableIndices()) == MOI.VariableIndex.(1:6)
    MOI.optimize!(bridged);
    @test sprint(show, optimizer) != nothing

    t_cold = MOI.get(bridged, MOI.SolveTime())
    iter_cold = MOI.get(bridged, COSMO.ADMMIterations())
    x_sol = MOI.get(bridged, MOI.VariablePrimal(), x)
    y_c1 = MOI.get(bridged, MOI.ConstraintDual(), con1)
    y_c2 = MOI.get(bridged, MOI.ConstraintDual(), con2)
    y_c3 = MOI.get(bridged, MOI.ConstraintDual(), con3)
    # Constraint Primal: For a constraint F-S, constraint primal is the value of the function evaluated at the primal solution
    s_c1 = MOI.get(bridged, MOI.ConstraintPrimal(), con1)
    s_c2 = MOI.get(bridged, MOI.ConstraintPrimal(), con2)
    s_c3 = MOI.get(bridged, MOI.ConstraintPrimal(), con3)

    @testset "SolverAttributes" begin
        @test MOI.get(bridged, MOI.SolverName()) == "COSMO";
        @test MOI.get(bridged, MOI.SolveTime()) >= 0.;
        # @test typeof(MOI.get(bridged, MOI.RawSolver())) <: COSMO.Workspace
        @test MOI.get(bridged, MOI.ResultCount()) == 1
        @test MOI.get(bridged, MOI.NumberOfVariables()) == 6
        @test isapprox(s_c1, b1, atol = 1e-3)
        @test isapprox(s_c2, b2, atol = 1e-3)
    end

    @testset "Warm starting" begin

        # solve the same problem again
        A1_t = [1.0; 0; 3; 2; 14; 5];
        A2_t = [0.0; 4; 6; 16; 0; 4];
        C_t = [1.; 4; 9; 6; 0; 7];

        cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
        optimizer =  COSMO.Optimizer(check_termination = 1, verbose = false);
        cached = MOIU.CachingOptimizer(cache, optimizer)
        bridged = MOIB.full_bridge_optimizer(cached, Float64)

        x = MOI.add_variables(bridged, 6);
        objectiveFunction = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(C_t, x[1:6]), 0.0);
        MOI.set(bridged, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objectiveFunction);
        MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE);
        con1 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A1_t, x[1:6]), 0.0), MOI.EqualTo(b1));
        con2 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A2_t, x[1:6]), 0.0), MOI.EqualTo(b2));
        con3 = MOI.add_constraint(bridged, MOI.VectorOfVariables([x[1], x[2], x[3], x[4], x[5], x[6]]), MOI.PositiveSemidefiniteConeTriangle(3));
        MOI.optimize!(bridged);

        x_sol_tri = MOI.get(bridged, MOI.VariablePrimal(), x)
        X_sol_tri = zeros(3, 3)
        optimizer = get_inner_optimizer(bridged)
        internal_scaled_s = copy(optimizer.inner.vars.s.data)
        internal_scaled_μ = copy(optimizer.inner.vars.μ)
        y_c3 = MOI.get(bridged, MOI.ConstraintDual(), con3)
        s_c3 = MOI.get(bridged, MOI.ConstraintPrimal(), con3)

        # warm start the PSD triangle related constraint primal and dual variables
        MOI.set(bridged, MOI.ConstraintDualStart(), con3, y_c3)
        MOI.set(bridged, MOI.ConstraintPrimalStart(), con3, s_c3)
        MOI.optimize!(bridged);

        optimizer = get_inner_optimizer(bridged)
        # @test isapprox(optimizer.inner.vars.μ[3:end], internal_scaled_μ[3:end], atol = 1e-3)
        # @test isapprox(optimizer.inner.vars.s.data[3:end], internal_scaled_s[3:end], atol = 1e-4)

        # Warm start triggered by copy_to function (warm starting values from the model)
        # Different QP problem:
        # min 1/2 x'Px + q'x
        # with  P = [4. 1;1 2]; q = [1; 1.]
        # s.t. Ax <= u   <=> -Ax + u in Nonnegatives
        #      Ax >= l   <=> Ax - l in Nonnegatives
        # with A = [1. 1;1 0; 0 1];
        l = [1.; 0; 0];
        u = [1.; 0.7; 0.7]

        cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
        optimizer =  COSMO.Optimizer();
        cached = MOIU.CachingOptimizer(cache, optimizer)
        bridged = MOIB.full_bridge_optimizer(cached, Float64)
        MOI.set(optimizer, MOI.Silent(), true)

        x = MOI.add_variables(bridged, 2);
        objectiveFunction = MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(1.0, x[1]); MOI.ScalarAffineTerm(1.0, x[2])], [MOI.ScalarQuadraticTerm(4.0, x[1], x[1]); MOI.ScalarQuadraticTerm(1.0, x[1], x[2]); MOI.ScalarQuadraticTerm(2.0, x[2], x[2])] , 0);
        MOI.set(bridged, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objectiveFunction);
        MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE);
        A = [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x[1])),MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x[2])),MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x[1])),MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, x[2]))];
        con1 = MOI.add_constraint(bridged, MOI.VectorAffineFunction(A, -u), MOI.Nonpositives(3));
        con2 = MOI.add_constraint(bridged, MOI.VectorAffineFunction(A, -l),MOI.Nonnegatives(3));

        MOI.optimize!(bridged);

        # solve once to get optimal solution
        x_sol = MOI.get(bridged, MOI.VariablePrimal(), x);
        y_c1 = MOI.get(bridged, MOI.ConstraintDual(), con1);
        y_c2 = MOI.get(bridged, MOI.ConstraintDual(), con2);
        optimizer = get_inner_optimizer(bridged)
        y_c1_rows = optimizer.rowranges[1];
        y_c2_rows = optimizer.rowranges[2];

        # provide warm start values to the model
        MOI.set.(bridged, MOI.VariablePrimalStart(), x, x_sol)
        MOI.set.(bridged, MOI.ConstraintDualStart(), [con1, con2], [y_c1, y_c2])
        MOI.empty!(optimizer);

        MOIU.attach_optimizer(bridged.model)
;

        # check that internal variables are set correctly
        @test optimizer.inner.vars.x == x_sol
        @test optimizer.inner.vars.μ[y_c1_rows] == y_c1
        @test optimizer.inner.vars.μ[y_c2_rows] == -y_c2

    end

    @testset "Small edge cases" begin
        cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
        optimizer =  COSMO.Optimizer(max_iter = 2, verbose = false);
        cached = MOIU.CachingOptimizer(cache, optimizer)
        bridged = MOIB.full_bridge_optimizer(cached, Float64)

        x = MOI.add_variables(bridged, 1);
        objf = MOI.ScalarAffineFunction{Float64}([MOI.ScalarAffineTerm(-1.0, x[1])] , 0);
        MOI.set(bridged, MOI.ObjectiveFunction{typeof(objf)}(), objf);
        MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], x[1:1]), 0.0), MOI.GreaterThan{Float64}(10.));

        MOI.optimize!(bridged);
        @test MOI.get(bridged, MOI.TerminationStatus()) == MOI.ITERATION_LIMIT
        @test MOI.get(bridged, MOI.PrimalStatus()) == MOI.NEARLY_FEASIBLE_POINT
        @test MOI.get(bridged, MOI.DualStatus()) == MOI.NEARLY_FEASIBLE_POINT

    end
    @testset "Set merging" begin
        # based on HS21 from Maros Meszaros QP Set but with a range of redundant constraints to test set merging
        # the original problem is:
        # min 1/2 x' P x + r
        # s.t. l <= Ax <= u with the following data
        P =[0.02 0; 0 2.0]
        r = -100.;
        A = [-10.0 1; -1 0; 0 -1]
        l = [ 10.0; 2; -50];
        u = [1e20; 50.0; 50.0]

        # the known solution is
        obj_true = -99.96
        x_true = [-2.; 0]

        cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
        optimizer =  COSMO.Optimizer(verbose = false);
        cached = MOIU.CachingOptimizer(cache, optimizer)
        bridged = MOIB.full_bridge_optimizer(cached, Float64)

        x = MOI.add_variables(bridged, 2);
        objectiveFunction = MOI.ScalarQuadraticFunction{Float64}([], [MOI.ScalarQuadraticTerm(0.02, x[1], x[1]); MOI.ScalarQuadraticTerm(2.0, x[2], x[2])], r);
        MOI.set(bridged, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objectiveFunction);
        MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE);

        # Add the constraints in an unsorted order to make merging more "challenging"
        # A[1, :] * x >= 10
        nn1 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[1, :], x[1:2]), 0.0), MOI.GreaterThan(l[1]));
        # Redundant SOC constraint: -x[1] >= ||x[2]||_2
        soc = MOI.add_constraint(bridged, MOI.VectorAffineFunction(MOI.VectorAffineTerm.( [1; 2], MOI.ScalarAffineTerm.([-1. ;1. ], x[1:2])), zeros(2)), MOI.SecondOrderCone(2));
        # Interval: 2 <= A[2, :] * x <= 50
        box1 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[2, :], x[1:2]), 0.0), MOI.Interval(l[2], u[2]));
        # Equal To: x[1] == - 2
        zero1 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], x[1:1]), 0.0), MOI.EqualTo(x_true[1]));
        # Nonpositives: x[1; 2] - [10; 10] <= [0; 0]
        nn2 = MOI.add_constraint(bridged, MOI.VectorAffineFunction(MOI.VectorAffineTerm.( [1; 2], MOI.ScalarAffineTerm.([1. ;1. ], x[1:2])), -10. .* ones(2)), MOI.Nonpositives(2));
        # Interval: -50 <= A[3, :] * x <= 50
        box2 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[3, :], x[1:2]), 0.0), MOI.Interval(l[3], u[3]));
        # ZeroSet: x[1; 2] - x_true == [0; 0]
        zeroset = MOI.add_constraint(bridged, MOI.VectorAffineFunction(MOI.VectorAffineTerm.( [1; 2], MOI.ScalarAffineTerm.([1. ;1. ], x[1:2])), -x_true), MOI.Zeros(2));

        # copy model into optimizer
        MOI.optimize!(bridged);

        # check that the sets have been merged
        optimizer = get_inner_optimizer(bridged)
        C = optimizer.inner.p.C
        @test length(C.sets) == 4
        @test typeof(C.sets[1]) <: COSMO.ZeroSet
        @test typeof(C.sets[2]) <: COSMO.Nonnegatives
        @test typeof(C.sets[3]) <: COSMO.Box
        @test typeof(C.sets[4]) <: COSMO.SecondOrderCone
        x_sol = MOI.get(optimizer, MOI.VariablePrimal(), x)
        @test isapprox( MOI.get(optimizer, MOI.ObjectiveValue()), obj_true, atol = 1e-3)
        @test norm(x_sol - x_true, Inf) < 1e-3

     end

     @testset "Unsupported Attributes" begin

        cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
        optimizer =  COSMO.Optimizer{Float64}();
        cached = MOIU.CachingOptimizer(cache, optimizer)
        bridged = MOIB.full_bridge_optimizer(cached, Float64)

         x = MOI.add_variable(bridged)
         y = MOI.add_variable(bridged)

         cf = MOI.ScalarAffineFunction{Float64}(MOI.ScalarAffineTerm{Float64}.([1., 1.], [x, y]), 0.)
         c = MOI.add_constraint(bridged, cf, MOI.LessThan(0.))

         MOI.set(bridged, UnsupportedModelAttribute(), 0.)
         @test_throws MOI.UnsupportedAttribute{UnsupportedModelAttribute} MOI.optimize!(bridged);

     end

     @testset "BigFloat Test" begin
         T = BigFloat

         # min 1/2 x'Px + q'x
         # with  P = [4. 1;1 2]; q = [1; 1.]
         # s.t. Ax <= u   <=> -Ax + u in Nonnegatives
         #      Ax >= l   <=> Ax - l in Nonnegatives
         # with A = [1. 1;1 0; 0 1];
         l = T[1.; 0; 0];
         u = T[1.; 0.7; 0.7]

        cache = MOIU.UniversalFallback(MOIU.Model{T}());
        optimizer =  COSMO.Optimizer{T}(verbose = false);
        cached = MOIU.CachingOptimizer(cache, optimizer)
        bridged = MOIB.full_bridge_optimizer(cached, T)

        MOI.set(bridged, MOI.Silent(), true)

         x = MOI.add_variables(bridged, 2);
         objectiveFunction = MOI.ScalarQuadraticFunction{T}([MOI.ScalarAffineTerm(T(1.0), x[1]); MOI.ScalarAffineTerm(T(1.0), x[2])], [MOI.ScalarQuadraticTerm(T(4.0), x[1], x[1]); MOI.ScalarQuadraticTerm(T(1.0), x[1], x[2]); MOI.ScalarQuadraticTerm(T(2.0), x[2], x[2])] , zero(T));
         MOI.set(bridged, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}(), objectiveFunction);
         MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE);
         A = [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1.0), x[1])),MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1.0), x[2])),MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(1.0), x[1])),MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(T(1.0), x[2]))];
         con1 = MOI.add_constraint(bridged, MOI.VectorAffineFunction(A, -u), MOI.Nonpositives(3));
         con2 = MOI.add_constraint(bridged, MOI.VectorAffineFunction(A, -l),MOI.Nonnegatives(3));

         MOI.optimize!(bridged);
         optimizer = get_inner_optimizer(bridged)
         @test optimizer.results.status == :Solved
         @test typeof(optimizer.results) == COSMO.Result{BigFloat}
         @test isapprox(optimizer.results.obj_val, 1.88, atol = 1e-3)
     end

    # # ------------------
    # # MOI - Test sets
    # --------------------

    optimizer = COSMO.Optimizer(eps_abs = 1e-5, eps_rel = 1e-5, check_termination = 25, check_infeasibility = 50)
    MOI.set(optimizer, MOI.Silent(), true)

    cache = MOIU.UniversalFallback(MOIU.Model{Float64}())
    cached = MOIU.CachingOptimizer(cache, optimizer)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)
    config = MOIT.TestConfig(atol = 1e-3, rtol = 1e-3, duals = false)

    @testset "Unit" begin
        MOIT.unittest(bridged, config,
                        [# Quadratic constraints are not supported
                        "solve_qcp_edge_cases",
                       # ArgumentError: The number of constraints in SCSModel must be greater than 0
                        "solve_unbounded_model",
                       # Integer and ZeroOne sets are not supported
                       "solve_integer_edge_cases", "solve_objbound_edge_cases",
                        "solve_zero_one_with_bounds_1",
                        "solve_zero_one_with_bounds_2",
                        "solve_zero_one_with_bounds_3",
                       # Can not set the # of threads
                        "number_threads",
                        "solve_farkas_interval_upper"
                         ])
    end

    @testset "Continuous Linear" begin
        MOIT.contlineartest(bridged, config, [
            "linear8a", # It expects `ResultCount` to be 0 as we disable `duals`.
            ])
    end

    @testset "Continuous Quadratic" begin

        MOIT.contquadratictest(bridged, config, [ "socp",
                                                 "ncqcp"])
    end

    exclude_conic_test_sets = ["rootdet", "logdet"]
    @testset "Continuous Conic" begin
        MOIT.contconictest(bridged, config, exclude_conic_test_sets)
    end

 end