text:
#=
    EodData Web Service Calls
=#

# ===================================
# Make functions available externally
export country_list, data_client_latest_version, data_formats, exchange_get, exchange_list,
        exchange_months, fundamental_list, login, login_2, membership, quote_get, quote_list,
        quote_list_2, quote_list_by_date, quote_list_by_date_2, quote_list_by_date_period,
        quote_list_by_date_period_2, split_list_by_exchange, split_list_by_symbol,
        symbol_changes_by_exchange, symbol_chart, symbol_get, symbol_history, symbol_history_period,
        symbol_history_period_by_date_range, symbol_list, symbol_list_2, technical_list,
        top_10_gains, top_10_losses, validate_access

# =========
# Functions

# country_list()
# --------------
# Returns a list of available countries.
# INPUT: Token (Login Token)
# OUTPUT: Dict() of countries of type ::Dict{String, Country}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=CountryList
function country_list(token::String)
    call = "/CountryList"
    args = ["Token"=>"$token"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_countries(xml_tree)
end

# data_client_latest_version()
# ----------------------------
# Returns the latest version information of Data Client.
# INPUT: Token (Login Token)
# OUTPUT: Date Client Version of type ::String
# REFERENCE: http://ws.eoddata.com/data.asmx?op=DataClientLatestVersion
function data_client_latest_version(token::String)
    call = "/DataClientLatestVersion"
    args = ["Token"=>"$token"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return find(xml_tree, "/RESPONSE/VERSION[1]#string")
end

# data_formats()
# --------------
# Returns the list of data formats.
# INPUT: Token (Login Token)
# OUTPUT: Dict() of DataFormats of the type ::Dict{String, DataFormat}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=DataFormats
function data_formats(token::String)
    call = "/DataFormats"
    args = ["Token"=>"$token"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_data_formats(xml_tree)
end

# exchange_get()
# --------------
# Returns detailed information of a specific exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Exchange of type ::Exchange
# REFERENCE: http://ws.eoddata.com/data.asmx?op=ExchangeGet
function exchange_get(token::String, exchange_code::String)
    call = "/ExchangeGet"
    args = ["Token"=>"$token", "Exchange"=>"$exchange_code"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return Exchange(find(xml_tree, "/RESPONSE/EXCHANGE[1]"))
end

# exchange_list()
# ---------------
# Returns a list of available exchanges.
# INPUT: Token (Login Token)
# OUTPUT: Dict() of exchanges of the type ::Dict{String, Exchange}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=ExchangeList
function exchange_list(token::String)
    call = "/ExchangeList"
    args = ["Token"=>"$token"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_exchanges(xml_tree)

end

# exchange_months()
# -----------------
# Returns the number of Months history a user is allowed to download.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Number of Months as an ::Int
# REFERENCE: http://ws.eoddata.com/data.asmx?op=ExchangeMonths
function exchange_months(token::String, exchange_code::String)
    call = "/ExchangeMonths"
    args = ["Token"=>"$token", "Exchange"=>"$exchange_code"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return int(find(xml_tree, "/RESPONSE/MONTHS[1]#string"))
end

# fundamental_list()
# ------------------
# Returns a complete list of fundamental data for an entire exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Dict() of fundamentals of type ::Dict{String, Fundamental}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=FundamentalList
function fundamental_list(token::String, exchange_code::String)
    call = "/FundamentalList"
    args = ["Token"=>"$token", "Exchange"=>"$exchange_code"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_fundamentals(xml_tree)
end

# login()
# -------
# Login to EODData Financial Information Web Service. Used for Web Authentication.
# INPUT: Username, Password
# OUTPUT: Login Token, which is a field in the type ::LoginResponse
# REFERENCE: http://ws.eoddata.com/data.asmx?op=Login
function login(username::String, password::String)
    call = "/Login"
    args = ["Username"=>"$username", "Password"=>"$password"]
    xml_tree = get_response(call, args)

    validate_xml_login(xml_tree) && return LoginResponse(xml_tree)
end

# login_2()
# ---------
# Login to EODData Financial Information Web Service. Used for Application Authentication.
# INPUT: Username, Password, Version (Application Version)
# OUTPUT: Login Token, which is a field in the type ::LoginResponse
# REFERENCE: http://ws.eoddata.com/data.asmx?op=Login2
function login_2(username::String, password::String, version::String)
    call = "/Login"
    args = ["Username"=>"$username", "Password"=>"$password", "Version"=>"$version"]
    xml_tree = get_response(call, args)

    validate_xml_login(xml_tree) && return LoginResponse(xml_tree)
end

# membership()
# ------------
# Returns membership of user.
# INPUT: Token (Login Token)
# OUTPUT: Membership of type ::String
# REFERENCE: http://ws.eoddata.com/data.asmx?op=Membership
function membership(token::String)
    call = "/Membership"
    args = ["Token"=>"$token"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return find(xml_tree, "/RESPONSE/MEMBERSHIP[1]#string")
end

# quote_get()
# -----------
# Returns an end of day quote for a specific ticker.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), ticker (eg:"MSFT")
# OUTPUT: End of day quote of type ::Quote
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteGet
function quote_get(token::String, exchange::String, ticker::String)
    call = "/QuoteGet"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return Quote(find(xml_tree, "/RESPONSE/QUOTE[1]"))
end

# quote_list()
# ------------
# Returns a complete list of end of day quotes for an entire exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Dict() of end of day quotes of type ::Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteList
function quote_list(token::String, exchange::String)
    call = "/QuoteList"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# quote_list_2()
# --------------
# Returns end of day quotes for a list of tickers of a specific exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), Symbols (eg:"MSFT,INTC")
# OUTPUT: Dict() of end of day quotes of type ::Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteList2
function quote_list_2(token::String, exchange::String, tickers::String)
    call = "/QuoteList2"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbols"=>"$tickers"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# quote_list_by_date()
# --------------------
# Returns a complete list of end of day quotes for an entire exchange and a specific date.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), QuoteDate (format:yyyyMMdd eg:"20080225")
# OUTPUT: Dict() of end of day quotes of type ::Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteListByDate
function quote_list_by_date(token::String, exchange::String, quote_date::String)
    call = "/QuoteListByDate"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "QuoteDate"=>"$quote_date"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# quote_list_by_date_2()
# ----------------------
# Returns a complete list of end of day quotes for an entire exchange and a specific date.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), QuoteDate (format:yyyyMMdd eg:"20080225")
# OUTPUT: Dict() of end of day quotes of type ::Dict{String, Quote_2}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteListByDate2
function quote_list_by_date_2(token::String, exchange::String, quote_date::String)
    call = "/QuoteListByDate2"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "QuoteDate"=>"$quote_date"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes_2(xml_tree)
end

# quote_list_by_date_period()
# ---------------------------
# Returns a complete list of end of day quotes for an entire exchange, specific date, and specific period.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), QuoteDate (format:yyyyMMdd eg:"20080225"),
#         Period (eg: "1", "5", "10", "15", "30", "h", "d", "w", "m")
# OUTPUT: Dict() of end of period quotes of type ::Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteListByDatePeriod
function quote_list_by_date_period(token::String, exchange::String, quote_date::String, period::String)
    call = "/QuoteListByDatePeriod"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "QuoteDate"=>"$quote_date", "Period"=>"$period"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# quote_list_by_date_period_2()
# -----------------------------
# Returns a complete list of end of day quotes for an entire exchange, specific date, and specific period.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), QuoteDate (format:yyyyMMdd eg:"20080225"),
# Period ("1", "5", "10", "15", "30", "h", "d", "w", "m")
# OUTPUT: Dict() of end of period quotes of type ::Dict{String, Quote_2}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=QuoteListByDatePeriod2
function quote_list_by_date_period_2(token::String, exchange::String, quote_date::String, period::String)
    call = "/QuoteListByDatePeriod2"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "QuoteDate"=>"$quote_date", "Period"=>"$period"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes_2(xml_tree)
end

# split_list_by_exchange()
# ------------------------
# Returns a list of Splits of a specific exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Dict() of splits of type::Dict{String, Split}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SplitListByExchange
function split_list_by_exchange(token::String, exchange::String)
    call = "/SplitListByExchange"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_splits(xml_tree)
end

# split_list_by_symbol()
# ----------------------
# Returns a list of splits of a specific ticker.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), Ticker (eg:"MSFT")
# OUTPUT: Dict() of splits of type::Dict{String, Split}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SplitListBySymbol
function split_list_by_symbol(token::String, exchange::String, ticker::String)
    call = "/SplitListBySymbol"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_splits(xml_tree)
end

# symbol_changes_by_exchange()
# ----------------------------
# Returns a list of ticker changes of a given exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Dict() of ticker changes of type::Dict{String, TickerChange}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolChangesByExchange
function symbol_changes_by_exchange(token::String, exchange::String)
    call = "/SymbolChangesByExchange"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_ticker_changes(xml_tree)
end

# symbol_chart()
# --------------
# Returns a URL to a chart image of a specific ticker.
# INPUT: Token (Login Token), Exchange (eg: NASDAQ), Ticker (eg:MSFT)
# OUTPUT: Chart URL
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolChart
function symbol_chart(token::String, exchange::String, ticker::String)
    symbol_get(token, exchange, ticker)
    call = "/SymbolChart"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return find(xml_tree, "/RESPONSE/CHART[1]#string")
end

# symbol_get()
# ------------
# Returns detailed information of a specific ticker.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), Ticker (eg:"MSFT")
# OUTPUT: Ticker of type ::Ticker
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolGet
function symbol_get(token::String, exchange::String, ticker::String)
    call = "/SymbolGet"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return Ticker(find(xml_tree, "/RESPONSE/SYMBOL[1]"))
end

# symbol_history()
# ----------------
# Returns a list of historical end of day data of a specified symbol and specified start date up to today's date.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), Ticker (eg:"MSFT"), StartDate (format:yyyyMMdd eg:"20080225")
# OUTPUT: Dict() of quotes of type ::Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolHistory
function symbol_history(token::String, exchange::String, ticker::String, start_date::String)
    call = "/SymbolHistory"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker", "StartDate"=>"$start_date"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# symbol_history_period()
# -----------------------
# Returns a list of historical data of a specified symbol, specified date and specified period.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), Symbol (eg:"MSFT"), Date (format:yyyyMMdd eg:"20080225"),
#           Period ("1", "5", "10", "15", "30", "h", "d", "w", "m")
# OUTPUT: Dict() of quotes of type ::Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolHistoryPeriod
function symbol_history_period(token::String, exchange::String, ticker::String, quote_date::String, period::String)
    call = "/SymbolHistoryPeriod"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker", "Date"=>"$quote_date", "Period"=>"$period"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# symbol_history_period_by_date_range()
# -------------------------------------
# Returns a list of historical data of a specified symbol, specified date range and specified period.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), Symbol (eg:"MSFT"),
#           StartDate (format:yyyyMMdd eg:"20080225"), EndDate (format:yyyyMMdd eg:"20080225"),
#           Period ("1", "5", "10", "15", "30", "h", "d", "w", "m")
# OUTPUT: Dict() of historical quotes of type Dict{String, Quote}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolHistoryPeriodByDateRange
function symbol_history_period_by_date_range(token::String, exchange::String, ticker::String, start_date::String, end_date::String, period::String)
    call = "/SymbolHistoryPeriodByDateRange"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "Symbol"=>"$ticker",
            "StartDate"=>"$start_date", "EndDate"=>"$end_date","Period"=>"$period"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# symbol_list()
# -------------
# Returns a list of symbols of a specified exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Dict() of tickers of type ::Dict{String, Ticker}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolList
function symbol_list(token::String, exchange::String)
    call = "/SymbolList"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_tickers(xml_tree)
end

# symbol_list_2()
# ---------------
# Returns a list of symbols of a specified exchange.
# INPUT: Token (Login Token), Exchange (eg: NASDAQ)
# OUTPUT: Dict() of tickers of type ::Dict{String, Ticker_2}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=SymbolList2
function symbol_list_2(token::String, exchange::String)
    call = "/SymbolList2"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_tickers_2(xml_tree)
end

# technical_list()
# ----------------
# Returns a complete list of technical data for an entire exchange.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ")
# OUTPUT: Dict() of technical indicators for each ticker of type ::Dict{String, Technical}
# REFERENCE: http://ws.eoddata.com/data.asmx?op=TechnicalList
function technical_list(token::String, exchange::String)
    call = "/TechnicalList"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_technicals(xml_tree)
end

# top_10_gains()
# --------------
# Returns a list of the Top 10 Gains of a specified exchange.
# INPUT: Token (Login Token), Exchange (eg: NASDAQ)
# OUTPUT: List of quotes
# REFERENCE: http://ws.eoddata.com/data.asmx?op=Top10Gains
function top_10_gains(token::String, exchange::String)
    call = "/Top10Gains"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# top_10_losses()
# ---------------
# Returns a list of the Top 10 Losses of a specified exchange.
# INPUT: Token (Login Token), Exchange (eg: NASDAQ)
# OUTPUT: List of quotes
# REFERENCE: http://ws.eoddata.com/data.asmx?op=Top10Losses
function top_10_losses(token::String, exchange::String)
    call = "/Top10Losses"
    args = ["Token"=>"$token", "Exchange"=>"$exchange"]
    xml_tree = get_response(call, args)

    validate_xml(xml_tree) && return set_quotes(xml_tree)
end

# update_data_format()
# --------------------
# Update preferred Data Format
# INPUT: Token (Login Token), IncludeHeader, IncludeSuffix
# OUTPUT: List of DataFormats
# REFERENCE: http://ws.eoddata.com/data.asmx?op=UpdateDataFormat
# function update_data_format()
    # This function is not implemented
# end

# validate_access()
# -----------------
# Validate access for an entire exchange, specific date, and specific period.
# INPUT: Token (Login Token), Exchange (eg: "NASDAQ"), QuoteDate (format:yyyyMMdd eg:"20080225"),
# Period ("1", "5", "10", "15", "30", "h", "d", "w", "m")
# OUTPUT: Boolean
# REFERENCE: http://ws.eoddata.com/data.asmx?op=ValidateAccess
function validate_access(token::String, exchange::String, quote_date::String, period::String)
    call = "/ValidateAccess"
    args = ["Token"=>"$token", "Exchange"=>"$exchange", "QuoteDate"=>"$quote_date", "Period"=>"$period"]
    xml_tree = get_response(call, args)

    return lowercase(strip(find(xml_tree, "/RESPONSE[1]{Message}"))) == "success" ? true : false
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Topology

A data structure for constructing topological relations in a [`Mesh`](@ref).

## References

* Floriani, L. & Hui, A. 2007. [Shape representations
  based on simplicial and cell complexes]
  (https://diglib.eg.org/handle/10.2312/egst.20071055.063-087)
"""
abstract type Topology end

"""
    vertices(topology)

Return the vertices of the `topology`.
"""
vertices(t::Topology) = 1:nvertices(t)

"""
    nvertices(topology)

Return the number of vertices in the `topology`.
"""
function nvertices(::Topology) end

"""
    faces(topology, rank)

Return an iterator with `rank`-faces of the `topology`.

## Examples

Consider a mesh of tetrahedrons embedded in a 3D space. We can loop over
all 3-faces (a.k.a. elements) or over all 2-faces to handle the interfaces
(i.e. triangles) between adjacent elements:

```julia
tetrahedrons = faces(topology, 3)
triangles    = faces(topology, 2)
segments     = faces(topology, 1)
```
"""
function faces(::Topology, rank) end

"""
    elements(topology)

Return the top-faces (a.k.a. elements) of the `topology`.

## Examples

The elements of a volume embedded in 3D space can be tetrahedrons, hexahedrons,
or any 3-face. The elements of a surface embedded in 3D space can be triangles,
quadrangles or any 2-face.
"""
elements(t::Topology) = (element(t, i) for i in 1:nelements(t))

"""
    element(topology, ind)

Return the element of the `topology` at index `ind`.
"""
function element(::Topology, ind) end

"""
    nelements(topology)

Return the number of elements in the `topology`.
"""
function nelements(::Topology) end

"""
    facets(topology)

Return the (top-1)-faces (a.k.a. facets) of the `topology`.
"""
facets(t::Topology) = (facet(t, i) for i in 1:nfacets(t))

"""
    facet(topology, ind)

Return the facet of the `topology` at index `ind`.
"""
function facet(::Topology, ind) end

"""
    nfacets(topology)

Return the number of facets in the `topology`.
"""
function nfacets(::Topology) end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("topologies/full.jl")
include("topologies/grid.jl")
include("topologies/halfedge.jl")

ext: jl


text:
@testset "RCF" begin
  Qx, x = PolynomialRing(FlintQQ)
  k, a = NumberField(x - 1, "a")
  Z = maximal_order(k)

  function doit(u::UnitRange, p::Int = 3)
    cnt = 0
    for i in u
      I = ideal(Z, i)
      r, mr = ray_class_group(I, n_quo=p)
      for s in index_p_subgroups(r, fmpz(p), (A,x) -> quo(A, x)[2])
        a = ray_class_field(mr, s)
        if isconductor(a, I, check=false)
          K = number_field(a)
          cnt += 1
        end
      end
    end
    return cnt
  end

  @test doit(1:100) == 16
  @test doit(10^18:10^18+100) == 18
  @test doit(10^18:10^18+1000, 11) == 2
end


text:
cd(@__DIR__); include("setups/rgc100.jl");
pyplot(dpi = 200)

##
ϵ = 0.3
pair_inds, 𝛟1 = find_pairinds(W; ϵ = ϵ)
v_bar = (𝛟1[pair_inds[1, end-1:end]] - 𝛟1[pair_inds[2, end-1:end]]) ./ 2
band = ϵ * norm(𝛟1, Inf)
pos_ar = findall(0 .< 𝛟1 .< band)
neg_ar = findall(-band .< 𝛟1 .<= 0)

# soft partition 1D
plt = plot(layout = Plots.grid(2, 1), size = (520, 500))
    plot!(-0.055:0.11:0.055, zeros(2), c = :black, lw = 2, yticks = false,
    subplot = 1)
    plot!(zeros(2), [-0.1, 0.1], lw = 3, c = :grey, linestyle = :solid,
    xlab = latexstring("\\phi^{rw}_1(v)"), subplot = 1)
    scatter!(𝛟1[𝛟1 .> 0], 0:0, ylim = [-0.2, 0.2], grid = false,
    ms = 4, mswidth = 0,  legend = false, frame = :box, c = :red, subplot = 1)
    scatter!(𝛟1[𝛟1 .< 0], 0:0, m = (1, 4), mswidth = 0, c = :blue, subplot = 1)
    plot!(ones(2) * ϵ * norm(𝛟1, Inf), [-0.08, 0.08], lw = 2, c = :grey,
    linestyle = :dash, subplot = 1)
    plot!(-ones(2) * ϵ * norm(𝛟1, Inf), [-0.08, 0.08], lw = 2, c = :grey,
    linestyle = :dash, subplot = 1)
    annotate!(0.02, -0.12,
    text(latexstring("\\epsilon \\cdot \\Vert \\phi^{rw}_1 \\Vert_{\\infty}"), 11),
    subplot = 1)
    annotate!(-0.02, -0.12,
    text(latexstring("-\\epsilon \\cdot \\Vert \\phi^{rw}_1 \\Vert_{\\infty}"), 11),
    subplot = 1)
    annotate!(0, 0.15, text("action region", 10), subplot = 1)

    plot!(-0.016:0.032:0.016, zeros(2), c = :black, lw = 2, yticks = false,
    subplot = 2)
    plot!(zeros(2), [-0.1, 0.1], lw = 3, c = :grey, linestyle = :solid,
    xlab = latexstring("\\phi^{rw}_1(v)"), subplot = 2)
    scatter!(𝛟1[pos_ar], 0:0, ylim = [-0.2, 0.2], grid = false,
    ms = 4, mswidth = 0,  legend = false, frame = :box, c = :red, subplot = 2)
    scatter!(𝛟1[neg_ar], 0:0, m = (1, 4), mswidth = 0, c = :blue, subplot = 2)
    plot!(ones(2) * ϵ * norm(𝛟1, Inf), [-0.08, 0.08], lw = 2, c = :grey,
    linestyle = :dash, subplot = 2)
    plot!(-ones(2) * ϵ * norm(𝛟1, Inf), [-0.08, 0.08], lw = 2, c = :grey,
    linestyle = :dash, subplot = 2)
    annotate!(0.012, -0.12,
    text(latexstring("\\epsilon \\cdot \\Vert \\phi^{rw}_1 \\Vert_{\\infty}"), 11),
    subplot = 2)
    annotate!(-0.012, -0.12,
    text(latexstring("-\\epsilon \\cdot \\Vert \\phi^{rw}_1 \\Vert_{\\infty}"), 11),
    subplot = 2)
    annotate!(-0.0075, 0.15, text("negative action region", 10), subplot = 2)
    annotate!(0.0075, 0.15, text("positive action region", 10), subplot = 2)
savefig(plt, "../figs/RGC100_fiedler_embedding_soft_clustering.png")

# zoom-in: find relection points about the nodal point
plt = plot(layout = Plots.grid(2, 1), size = (520, 500))
    plot!(-0.005:0.01:0.005, zeros(2), c = :black, lw = 2, yticks = false,
    subplot = 1)
    scatter!(𝛟1[pair_inds[1, end-1:end]], 0:0, ylim = [-0.2, 0.2], grid = false,
    ms = 6, mswidth = 0,  legend = false, frame = :box, c = :red, subplot = 1)
    scatter!(𝛟1[pair_inds[2, end-1:end]], 0:0, ms = 6, mswidth = 0, c = :blue,
    subplot = 1)
    plot!(zeros(2), [-0.1, 0.1], lw = 3, c = :grey, linestyle = :solid,
    xlab = latexstring("\\phi^{rw}_1(v)"), subplot = 1)
    annotate!(𝛟1[pair_inds[1, end]], 0.05, text(latexstring("v^{+}_{1}"), 11),
    subplot = 1)
    annotate!(𝛟1[pair_inds[1, end-1]], 0.05, text(latexstring("v^{+}_{2}"), 11),
    subplot = 1)
    annotate!(𝛟1[pair_inds[2, end]], 0.05, text(latexstring("v^{-}_{1}"), 11),
    subplot = 1)
    annotate!(𝛟1[pair_inds[2, end-1]], 0.05, text(latexstring("v^{-}_{2}"), 11),
    subplot = 1)

    plot!(vec(-0.005:0.01:0.005), zeros(2), c = :black, lw = 2,
    yticks = false, subplot = 2)
    scatter!(vcat(v_bar, -v_bar), 0:0, ylim = [-0.2, 0.2], grid = false,
    ms = 6, mswidth = 0,  legend = false, frame = :box, c = :purple, subplot = 2)
    plot!(zeros(2), [-0.1, 0.1], lw = 3, c = :grey, linestyle = :solid,
    xlab = latexstring("r"), subplot = 2)
    annotate!(v_bar[end], -0.05, text(latexstring("r_{1}"), 11),
    subplot = 2)
    annotate!(v_bar[end-1], -0.05, text(latexstring("r_{2}"), 11),
    subplot = 2)
    annotate!(-v_bar[end], -0.05, text(latexstring("-r_{1}"), 11),
    subplot = 2)
    annotate!(-v_bar[end-1], -0.05, text(latexstring("-r_{2}"), 11),
    subplot = 2)
savefig(plt, "../figs/RGC100_fiedler_embedding_find_reflection_points.png")


text:
# Autogenerated wrapper script for Emoslib_jll for i686-linux-gnu-libgfortran3
export emos

using eccodes_jll
using FFTW_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Emoslib")
JLLWrappers.@declare_file_product(emos)
function __init__()
    JLLWrappers.@generate_init_header(eccodes_jll, FFTW_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_file_product(
        emos,
        "lib/libemos.a",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function finitehess{T <: Number}(
    f::Function,
    x::T,
    direction::Symbol = :central,
)
    ε = @hessianrule(x)
    return (f(x + ε) - 2 * f(x) + f(x - ε)) / (ε * ε)
end

function finitehess!{T <: Number}(
    H::Array{T},
    f::Function,
    x::Vector{T},
    direction::Symbol = :central,
)
    n = length(x)
    ε = nan(T)
    # TODO: Remove all these copies
    xpp, xpm, xmp, xmm = copy(x), copy(x), copy(x), copy(x)
    f_x = f(x)
    for i in 1:n
        x_i = x[i]
        ε = @hessianrule(x[i])
        xpp[i], xmm[i] = x_i + ε, x_i - ε
        H[i, i] = (f(xpp) - 2 * f_x + f(xmm)) / (ε * ε)
        ε_i = @centralrule(x[i])
        xp = x_i + ε_i
        xm = x_i - ε_i
        xpp[i], xpm[i], xmp[i], xmm[i] = xp, xp, xm, xm
        for j in (i + 1):n
            x_j = x[j]
            ε_j = @centralrule(x[j])
            xp = x_j + ε_j
            xm = x_j - ε_j
            xpp[j], xpm[j], xmp[j], xmm[j] = xp, xm, xp, xm
            H[i, j] = (f(xpp) - f(xpm) - f(xmp) + f(xmm)) / (4 * ε_i * ε_j)
            xpp[j], xpm[j], xmp[j], xmm[j] = x_j, x_j, x_j, x_j
        end
        xpp[i], xpm[i], xmp[i], xmm[i] = x_i, x_i, x_i, x_i
    end
    Base.LinAlg.copytri!(H, 'U')
    return H
end

function finitehess{T <: Number}(
    f::Function,
    x::Vector{T},
    direction::Symbol = :central
)
    n = length(x)
    H = Array(T, n, n)
    finitehess!(f, x, H, direction)
    return H
end

function taylor_finite_difference_hessian(
    f::Function,
    x::Real,
    h::Real,
)
    f_x = f(x)
    d = 4^6 * (2^4 * (f(x + h) + f(x - h) - 2 * f_x) - (f(x + 2 * h) + f(x - 2 * h) - 2 * f_x))
    d += - (2^4 * (f(x + 4 * h) + f(x - 4 * h) - 2 * f_x) - (f(x + 8 * h) + f(x - 8 * h) - 2 * f_x))
    return d / (3 * 2^6 * (2^8 - 1) * h^2)
end


text:
using Plots

include("models/Dicke.jl")
include("modules/ClassicalDynamics.jl")

dickeParameters = [2, 0.5, 1, 1]
energy = -1
numTrajectories = 20

savePath = "d:\\results\\Dicke"

pyplot(size = (1800, 1000))
fig = LyapunovExponents(energy, dickeParameters, numTrajectories; verbose=true)
savefig(fig, savePath * "_PS_E=$(energy)_$dickeParameters.png")

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
@testset "segwit" begin

    # Points to the very first segwit transaction, probably only works on my
    # personal copy of the blockchain:
    #
    # The block is at height 481825 (one based, if you ask bitcoin-cli, height is
    # 481824)
    link = Link(
        BTCParser.to_unsigned((0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
                               0x00, 0x1c, 0x80, 0x18, 0xd9, 0xcb, 0x3b, 0x74,
                               0x2e, 0xf2, 0x51, 0x14, 0xf2, 0x75, 0x63, 0xe3,
                               0xfc, 0x4a, 0x19, 0x02, 0x16, 0x7f, 0x98, 0x93)[end:-1:1]),
        UInt64(976),
        UInt64(0x577479)
    )

    Block(link)

end


text:
using Singular
using SignatureGB
R, (x1, x2, x3, x4, x5, x6, x7, x8) = Singular.PolynomialRing(Fp(101), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"])
I = spoly{n_Zp}[x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8, x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 + x6*x7 + x1*x8 + x7*x8, x1*x2*x3 + x2*x3*x4 + x3*x4*x5 + x4*x5*x6 + x5*x6*x7 + x1*x2*x8 + x1*x7*x8 + x6*x7*x8, x1*x2*x3*x4 + x2*x3*x4*x5 + x3*x4*x5*x6 + x4*x5*x6*x7 + x1*x2*x3*x8 + x1*x2*x7*x8 + x1*x6*x7*x8 + x5*x6*x7*x8, x1*x2*x3*x4*x5 + x2*x3*x4*x5*x6 + x3*x4*x5*x6*x7 + x1*x2*x3*x4*x8 + x1*x2*x3*x7*x8 + x1*x2*x6*x7*x8 + x1*x5*x6*x7*x8 + x4*x5*x6*x7*x8, x1*x2*x3*x4*x5*x6 + x2*x3*x4*x5*x6*x7 + x1*x2*x3*x4*x5*x8 + x1*x2*x3*x4*x7*x8 + x1*x2*x3*x6*x7*x8 + x1*x2*x5*x6*x7*x8 + x1*x4*x5*x6*x7*x8 + x3*x4*x5*x6*x7*x8, x1*x2*x3*x4*x5*x6*x7 + x1*x2*x3*x4*x5*x6*x8 + x1*x2*x3*x4*x5*x7*x8 + x1*x2*x3*x4*x6*x7*x8 + x1*x2*x3*x5*x6*x7*x8 + x1*x2*x4*x5*x6*x7*x8 + x1*x3*x4*x5*x6*x7*x8 + x2*x3*x4*x5*x6*x7*x8, x1*x2*x3*x4*x5*x6*x7*x8 - 1]
comp_ideal = gens(R)
SignatureGB.f5(comp_ideal)
gb = SignatureGB.f5(I, verbose = true)
println(SignatureGB.is_gb(gb))


text:
using ExtendableGrids
import Makie, PyPlot,Plots,VTKView

X=collect(0.0:0.05:1)

grid1d=simplexgrid(X)
grid2d=simplexgrid(X,X)

func1d=map(x->sin(5*x),grid1d)
func2d=map((x,y)->sin(5*x)*sin(5*y),grid2d)

function gridloop(Plotter)
    ctx=PlotterContext(Plotter)
    @time for i=1:20
        X=collect(0.0:0.05:1)
        Y=collect(0.0:0.05:1+0.01*i)
        grid=simplexgrid(X,Y)
        plot!(ctx,grid)
    end
end

function gridloops()
    gridloop(PyPlot)
    gridloop(VTKView)
#    gridloop(Plots)
    gridloop(Makie)
end

function funcloop(Plotter)
    ctx=PlotterContext(Plotter)
    X=collect(0.0:0.05:1)
    grid=simplexgrid(X,X)
    @time for i=1:20
        func2d=map((x,y)->sin(0.5*i*x)*sin(0.5*i*y),grid2d)
        plot!(ctx,grid,func2d)
    end
end

function funcloops()
    funcloop(PyPlot)
    funcloop(VTKView)
    funcloop(Plots)
    funcloop(Makie)
end

twofigloop(Plotter)=twofigloop(plottertype(Plotter), Plotter)

function twofigloop(::Type{MakieType},Makie)
    MakieLayout=Makie.AbstractPlotting.MakieLayout
    ## Layouting mit leeren scenen auch nicht
    ## tricontour mit linien ?
    scene, layout = MakieLayout.layoutscene(resolution = (600,300))
    ax1=layout[1, 1] = MakieLayout.LAxis(scene)
    ax2=layout[1, 2] = MakieLayout.LAxis(scene)
    Makie.lines!(ax1,[0.0,0.0001],[0.0,0.0001],color=:white)
    Makie.lines!(ax2,[0.0,0.0001],[0.0,0.0001],color=:white)
    ctx1=PlotterContext(Makie,scene=ax1)
    ctx2=PlotterContext(Makie,scene=ax2)
    Makie.display(scene)

    X=collect(0.0:0.05:1)
    grid=simplexgrid(X,X)

    for i=1:200
        func1=map((x,y)->sin(0.1*i*x)*sin(0.1*i*y),grid)
        func2=map((x,y)->exp(sin(0.1*i*x)*sin(0.1*i*y)),grid)
        plot!(ctx1,grid,func1,show=false)
        plot!(ctx2,grid,func2,show=false)
        Makie.update!(scene)
        Makie.sleep(1.0e-10)
    end
    scene
end

#
# Works ok this way
#
function twofigloop(::Type{PyPlotType},PyPlot)

    figure=PyPlot.figure(1,figsize=(6,3))
    ctx1=PlotterContext(PyPlot,figure=figure,clear=false,legend=false,subplot=121)
    ctx2=PlotterContext(PyPlot,figure=figure,clear=false,legend=false,subplot=122)

    for i=1:20
        X=collect(0.0:0.05:1-0.01*i)
        Y=collect(0.0:0.05:1+0.01*i)
        grid=simplexgrid(X,Y)
        PyPlot.clf()
        func2d=map((x,y)->sin(0.5*i*x)*sin(0.5*i*y),grid)
        plot!(ctx1,grid)
        plot!(ctx2,grid,func2d)
    end
end

function twowinloop()
    # in order to make this work we neet to switch form
    # PyPlot.plot to ax.plot.
    ctx1=PlotterContext(PyPlot,fignumber=1,legend=false)
    ctx2=PlotterContext(PyPlot,fignumber=2,legend=false)

    for i=1:20
        X=collect(0.0:0.05:1-0.01*i)
        Y=collect(0.0:0.05:1+0.01*i)
        grid=simplexgrid(X,Y)
        PyPlot.clf()
        func2d=map((x,y)->sin(0.5*i*x)*sin(0.5*i*y),grid)
        plot!(ctx1,grid)
        plot!(ctx2,grid,func2d)
    end
end

#
# Warnings from vtk when updating grid + functions
#
function twofigloop(::Type{VTKViewType},VTKView)

    frame=VTKView.StaticFrame()
    VTKView.clear!(frame)
    VTKView.size!(frame,600,300)
    VTKView.layout!(frame,2,1)
    dataset=VTKView.DataSet()
    ctx1=PlotterContext(VTKView,frame=frame,clear=false,dataset=dataset,framepos=1, label="func1")
    ctx2=PlotterContext(VTKView,frame=frame,clear=false,dataset=dataset,framepos=2, label="func2")

    X=collect(0.0:0.05:1)
    grid=simplexgrid(X,X)
    for i=1:200
        func1=map((x,y)->sin(0.1*i*x)*sin(0.1*i*y),grid)
        func2=map((x,y)->exp(sin(0.1*i*x)*sin(0.1*i*y)),grid)
        plot!(ctx1,grid,func1)
        plot!(ctx2,grid,func2)
    end
end

#
# Works well this way, missing tricontourf and trigrid plotting...
#
function twofigloop(::Type{PlotsType},Plots)

    ctx1=PlotterContext(Plots,plot=Plots.plot(),show=false,clear=false)
    ctx2=PlotterContext(Plots,plot=Plots.plot(),show=false,clear=false)
    p=Plots.plot(ctx1[:plot],ctx2[:plot])

    X=collect(0.0:0.05:1)
    Y=collect(0.0:0.05:1)
    grid=simplexgrid(X,Y)
    for i=1:20
        func1=map((x,y)->sin(0.5*i*x)*sin(0.5*i*y),grid)
        func2=map((x,y)->exp(sin(0.5*i*x)*sin(0.5*i*y)),grid)
        plot!(ctx1,grid,func1)
        plot!(ctx2,grid,func2)

        Plots.gui(p)
    end
end

function twofigloops()
    twofigloop(PyPlot)
    twofigloop(VTKView)
    twofigloop(Plots)
    twofigloop(Makie)
end

function alltests()
    gridloops()
    funcloops()
    twofigloops()
end
"""
Layout plan:

PlotterContext(layout=(2,3)) creates array of contexts
ctx[1,2]  contains subcontexts for ploting into

update!(ctx) updates at end of two figure loop

"""


text:
using Measures

include("helper_functions.jl")

begin
    λ1 = 1.4
    λ2 = 1
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    qs = 0:0.001:1
end

plot(qs, pstar.(qs, parameters...), label = "Best response (p*)", lw = 1.5)
plot!(xlims = (0, 1), ylims = (0, 1))
μ_boolvec, ν_boolvec = calculate_assumptions(pstar.(qs, parameters...), qs, parameters...)
plot!(shapes_from_intervals(shaded_regions(μ_boolvec, qs)), alpha = 0.3, fc = :red, label = "Assumptions violated")
plot!(xlabel = "Opponent probability (q)", ylabel = "Own probability (p)")
plot!(size=(600,300), margin = 2mm)

# savefig("best_response_plot.pdf")


text:
module RegularizedLeastSquares

import Base: length, iterate
using LinearAlgebra
import LinearAlgebra.BLAS: gemv, gemv!
import LinearAlgebra: BlasFloat, normalize!, norm, rmul!, lmul!
using SparseArrays
using IterativeSolvers
using Random
using SparsityOperators
using VectorizationBase
using VectorizationBase: shufflevector, zstridedpointer
@reexport using SparsityOperators
@reexport using ProgressMeter

export createLinearSolver, init, deinit, solve, linearSolverList,linearSolverListReal

abstract type AbstractLinearSolver end
const Trafo = Union{AbstractMatrix, AbstractLinearOperator, Nothing}
const FuncOrNothing = Union{Function, Nothing}

# Fallback function
setlambda(S::AbstractMatrix, λ::Real) = nothing

include("proximalMaps/ProxL1.jl")
include("proximalMaps/ProxL2.jl")
include("proximalMaps/ProxL21.jl")
include("proximalMaps/ProxLLR.jl")
# include("proximalMaps/ProxSLR.jl")
include("proximalMaps/ProxPositive.jl")
include("proximalMaps/ProxProj.jl")
include("proximalMaps/ProxTV.jl")
include("proximalMaps/ProxTVCondat.jl")
include("proximalMaps/ProxNuclear.jl")

include("Regularization.jl")
include("Utils.jl")
include("Kaczmarz.jl")
include("DAXKaczmarz.jl")
include("DAXConstrained.jl")
include("CGNR.jl")
include("Direct.jl")
include("FusedLasso.jl")
include("FISTA.jl")
include("ADMM.jl")
include("SplitBregman.jl")
include("PrimalDualSolver.jl")

"""
Return a list of all available linear solvers
"""
function linearSolverList()
  Any["kaczmarz","cgnr"] # These are those passing the tests
    #, "fusedlasso"]
end

function linearSolverListReal()
  Any["kaczmarz","cgnr","daxkaczmarz","daxconstrained","primaldualsolver"] # These are those passing the tests
    #, "fusedlasso"]
end

"""
    createLinearSolver(solver::AbstractString, A; log::Bool=false, kargs...)

This method creates a solver. The supported solvers are methods typically used in MPI/MRI.
All solvers return an approximate solution to STx = u.
Function returns choosen solver.

# solvers:
* `"kaczmarz"`        - kaczmarz method (the default)
* `"cgnr`             - CGNR
* `"direct"`          - A direct solver using the backslash operator
* `"daxkaczmarz"`     - Dax algorithm (with Kaczmarz) for unconstrained problems
* `"daxconstrained"`  - Dax algorithm for constrained problems
* `"pseudoinverse"`   - approximates a solution using the More-Penrose pseudo inverse
* `"fusedlasso"`      - solver for the Fused-Lasso problem
* `"fista"`           - Fast Iterative Shrinkage Thresholding Algorithm
* `"admm"`            - Alternating Direcion of Multipliers Method
* `"splitBregman"`    - Split Bregman method for constrained & regularized inverse problems
* `"primaldualsolver"`- First order primal dual method
"""
function createLinearSolver(solver::AbstractString, A, x=zeros(eltype(A),size(A,2));
                            log::Bool=false, kargs...)

  log ? solverInfo = SolverInfo(;kargs...) : solverInfo=nothing

  if solver == "kaczmarz"
    return Kaczmarz(A; kargs...)
  elseif solver == "cgnr"
    return CGNR(A, x; kargs...)
  elseif solver == "direct"
    return DirectSolver(A; kargs...)
  elseif solver == "daxkaczmarz"
    return DaxKaczmarz(A; kargs...)
  elseif solver == "daxconstrained"
    return DaxConstrained(A; kargs...)
  elseif solver == "pseudoinverse"
    return PseudoInverse(A; kargs...)
  elseif solver == "fusedlasso"
    return FusedLasso(A; kargs...)
  elseif solver == "fista"
    return FISTA(A, x; kargs...)
  elseif solver == "admm"
    return ADMM(A, x; kargs...)
  elseif solver == "splitBregman"
    return SplitBregman(A, x; kargs...)
  elseif solver == "primaldualsolver"
    return PrimalDualSolver(A; kargs...)
  else
    error("Solver $solver not found.")
    return Kaczmarz(A; kargs...)
  end
end

end


text:
function download_(filenameServer, filenameLocal)
  if !isfile(filenameLocal)
    @info "download $(filenameLocal)..."
    HTTP.open("GET", "http://media.tuhh.de/ibi/openMPIData/data/"*filenameServer) do http
      open(filenameLocal, "w") do file
        write(file, http)
      end
    end
  end
end

mkpath("./data/")

# Download data
download_("calibrations/16.mdf", "./data/SM.mdf")


text:
using HPAT

@acc hpat function movingavg_test(file_name)
    df = DataSource(DataTable{:id=Int64,:x=Float64,:y=Float64}, HDF5, file_name)
    avg = stencil(x->(x[-1]+x[0]+x[1])/3.0, df[:x])
    return sum(avg)
end

movingavg_test(ENV["HOME"]*"/tmp/cumsum.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
using PolaronPathIntegrals
using Test

@testset "PolaronPathIntegrals" begin
    include("test_coupling.jl")
    include("test_free_energy.jl")
    include("test_variation.jl")
end


text:
function load_netlist(netlist::AbstractArray{T}) where T <: AbstractString
    append!(netlist, ["C_NULL"])
    t = netlist |> ngSpice_Circ
end

function load_netlist(netlist::AbstractString)
    # The strings in Julia always use \n instead of \r \r\n
    # (even if the editor does).
    # But if the circuit-string has \r in it,
    # uncomment the following line
    # occursin("\r", netlist) && replace(netlist, "\r"=>"")
    load_netlist(split(netlist, "\n"))
end

source(netpath) = cmd("source $netpath")

display() = cmd("display")

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
#=**************************************************************************
This is a wrapper for a photochemical mechanism used to generate training
data for machine learning and restricted inverse algorithms.  It loads the
module "FormO3NOx" as mechanism.gas.  An adjustable parameter is number of days
to run the mechanism.

last edited June 30 2020
Obin Sturm posturm@ucdavis.edu
**************************************************************************=#

#mechanism = include("/Users/obin/Dropbox/ML Share/Photochemical Mechanisms/FormO3NOx/assign3_driver.jl
# For when in the same path
mechanism = include("assign3_driver.jl")
using DelimitedFiles

# The system ∆C = AS
A_reduced = [ 0  1 -1  0  0  0  0  0  0  0;  # O3
              1  0 -1  0  0  0 -1  0  0  0;  # NO
             -1  0  1  0  0  0  1 -1  0  0;  # NO2
              0  0  0 -1 -1 -1  0  0  0  0;  # HCHO
              0  0  0  2  0  1 -1  0  0  1;  # HO2.
              0  0  0  0  0  0  0  0 -1 -1];  # HO2H
#              0  0  0  0  0 -1  1 -1  2 -1;  # OH.
#              1 -1  0  0  0  0  0  0  0  0;  # O
#              0  0  0  0  0  0  0  1  0  0;  # HNO3
#              0  0  0  1  1  1  0  0  0  0;  # CO
#              0  0  0  0  1  0  0  0  0  0]; # H2

# Change this to adjust the number of full days the mechanism should run for
numdays = 365

# P
C,S,J = mechanism.gas(0.0,24.0*60,6.0,numdays)

# Create a data structure for day, time, daylight constant, and concentrations
Cday = repeat(1:numdays, inner = 241)
Ctime = repeat(collect(0.0:6.0:1440.0), outer = numdays)
Clabel = permutedims(["Day","Time [min]","Daylight Constant [0-1]","O3","NO","NO2","HCHO","HO2.","HO2H","OH.","O","HNO3","CO","H2"])
C_and_J = vcat(Clabel, hcat(Cday,Ctime,J,C))

# Create a data structure for S, related to adjacent C's as A*S = (C(t+∆t)-C(∆t))
Sday = repeat(1:numdays, inner = 240)
Stime = repeat(collect(0.0:6.0:1434.0), outer = numdays)
Slabel = permutedims(["Day","Time [min]","R1","R2","R3","R4","R5","R6","R7","R8","R9","R10"])
S = vcat(Slabel, hcat(Sday,Stime,S))

open("C_and_J.txt", "w") do io
                        writedlm(io, C,',')
                  end
open("S.txt", "w") do io
                        writedlm(io, S,',')
                  end


text:
struct enum_fdb_network_option
    LOCAL_ADDRESS::Cint
    CLUSTER_FILE::Cint
    TRACE_ENABLE::Cint
    TRACE_ROLL_SIZE::Cint
    TRACE_MAX_LOGS_SIZE::Cint
    TRACE_LOG_GROUP::Cint
    KNOB::Cint
    TLS_PLUGIN::Cint
    TLS_CERT_BYTES::Cint
    TLS_CERT_PATH::Cint
    TLS_KEY_BYTES::Cint
    TLS_KEY_PATH::Cint
    TLS_VERIFY_PEERS::Cint
    BUGGIFY_ENABLE::Cint
    BUGGIFY_DISABLE::Cint
    BUGGIFY_SECTION_ACTIVATED_PROBABILITY::Cint
    BUGGIFY_SECTION_FIRED_PROBABILITY::Cint
    DISABLE_MULTI_VERSION_CLIENT_API::Cint
    CALLBACKS_ON_EXTERNAL_THREADS::Cint
    EXTERNAL_CLIENT_LIBRARY::Cint
    EXTERNAL_CLIENT_DIRECTORY::Cint
    DISABLE_LOCAL_CLIENT::Cint
    DISABLE_CLIENT_STATISTICS_LOGGING::Cint
    ENABLE_SLOW_TASK_PROFILING::Cint
end

struct enum_fdb_cluster_option
    DUMMY_DO_NOT_USE::Cint
end

struct enum_fdb_database_option
    LOCATION_CACHE_SIZE::Cint
    MAX_WATCHES::Cint
    MACHINE_ID::Cint
    DATACENTER_ID::Cint
end

struct enum_fdb_transaction_option
    CAUSAL_WRITE_RISKY::Cint
    CAUSAL_READ_RISKY::Cint
    CAUSAL_READ_DISABLE::Cint
    NEXT_WRITE_NO_WRITE_CONFLICT_RANGE::Cint
    COMMIT_ON_FIRST_PROXY::Cint
    CHECK_WRITES_ENABLE::Cint
    READ_YOUR_WRITES_DISABLE::Cint
    READ_AHEAD_DISABLE::Cint
    DURABILITY_DATACENTER::Cint
    DURABILITY_RISKY::Cint
    DURABILITY_DEV_NULL_IS_WEB_SCALE::Cint
    PRIORITY_SYSTEM_IMMEDIATE::Cint
    PRIORITY_BATCH::Cint
    INITIALIZE_NEW_DATABASE::Cint
    ACCESS_SYSTEM_KEYS::Cint
    READ_SYSTEM_KEYS::Cint
    DEBUG_DUMP::Cint
    DEBUG_RETRY_LOGGING::Cint
    TRANSACTION_LOGGING_ENABLE::Cint
    TIMEOUT::Cint
    RETRY_LIMIT::Cint
    MAX_RETRY_DELAY::Cint
    SNAPSHOT_RYW_ENABLE::Cint
    SNAPSHOT_RYW_DISABLE::Cint
    LOCK_AWARE::Cint
    USED_DURING_COMMIT_PROTECTION_DISABLE::Cint
    READ_LOCK_AWARE::Cint
end

struct enum_fdb_streaming_mode
    WANT_ALL::Cint
    ITERATOR::Cint
    EXACT::Cint
    SMALL::Cint
    MEDIUM::Cint
    LARGE::Cint
    SERIAL::Cint
end

struct enum_fdb_mutation_type
    ADD::Cint
    AND::Cint
    BIT_AND::Cint
    OR::Cint
    BIT_OR::Cint
    XOR::Cint
    BIT_XOR::Cint
    MAX::Cint
    MIN::Cint
    SET_VERSIONSTAMPED_KEY::Cint
    SET_VERSIONSTAMPED_VALUE::Cint
    BYTE_MIN::Cint
    BYTE_MAX::Cint
end

struct enum_fdb_conflict_range_type_t
    READ::Cint
    WRITE::Cint
end

struct enum_fdb_error_predicate_t
    RETRYABLE::Cint
    MAYBE_COMMITTED::Cint
    RETRYABLE_NOT_COMMITTED::Cint
end

const FDBNetworkOption = enum_fdb_network_option(
    Cint(10), # FDB_NET_OPTION_LOCAL_ADDRESS
    Cint(20), # FDB_NET_OPTION_CLUSTER_FILE
    Cint(30), # FDB_NET_OPTION_TRACE_ENABLE
    Cint(31), # FDB_NET_OPTION_TRACE_ROLL_SIZE
    Cint(32), # FDB_NET_OPTION_TRACE_MAX_LOGS_SIZE
    Cint(33), # FDB_NET_OPTION_TRACE_LOG_GROUP
    Cint(40), # FDB_NET_OPTION_KNOB
    Cint(41), # FDB_NET_OPTION_TLS_PLUGIN
    Cint(42), # FDB_NET_OPTION_TLS_CERT_BYTES
    Cint(43), # FDB_NET_OPTION_TLS_CERT_PATH
    Cint(45), # FDB_NET_OPTION_TLS_KEY_BYTES
    Cint(46), # FDB_NET_OPTION_TLS_KEY_PATH
    Cint(47), # FDB_NET_OPTION_TLS_VERIFY_PEERS
    Cint(48), # FDB_NET_OPTION_BUGGIFY_ENABLE
    Cint(49), # FDB_NET_OPTION_BUGGIFY_DISABLE
    Cint(50), # FDB_NET_OPTION_BUGGIFY_SECTION_ACTIVATED_PROBABILITY
    Cint(51), # FDB_NET_OPTION_BUGGIFY_SECTION_FIRED_PROBABILITY
    Cint(60), # FDB_NET_OPTION_DISABLE_MULTI_VERSION_CLIENT_API
    Cint(61), # FDB_NET_OPTION_CALLBACKS_ON_EXTERNAL_THREADS
    Cint(62), # FDB_NET_OPTION_EXTERNAL_CLIENT_LIBRARY
    Cint(63), # FDB_NET_OPTION_EXTERNAL_CLIENT_DIRECTORY
    Cint(64), # FDB_NET_OPTION_DISABLE_LOCAL_CLIENT
    Cint(70), # FDB_NET_OPTION_DISABLE_CLIENT_STATISTICS_LOGGING
    Cint(71)  # FDB_NET_OPTION_ENABLE_SLOW_TASK_PROFILING
)

const FDBClusterOption = enum_fdb_cluster_option(
    Cint(-1) # FDB_CLUSTER_OPTION_DUMMY_DO_NOT_USE
)

const FDBDatabaseOption = enum_fdb_database_option(
    Cint(10), # FDB_DB_OPTION_LOCATION_CACHE_SIZE
    Cint(20), # FDB_DB_OPTION_MAX_WATCHES
    Cint(21), # FDB_DB_OPTION_MACHINE_ID
    Cint(22)  # FDB_DB_OPTION_DATACENTER_ID
)

const FDBTransactionOption = enum_fdb_transaction_option(
    Cint(10),  # FDB_TR_OPTION_CAUSAL_WRITE_RISKY
    Cint(20),  # FDB_TR_OPTION_CAUSAL_READ_RISKY
    Cint(21),  # FDB_TR_OPTION_CAUSAL_READ_DISABLE
    Cint(30),  # FDB_TR_OPTION_NEXT_WRITE_NO_WRITE_CONFLICT_RANGE
    Cint(40),  # FDB_TR_OPTION_COMMIT_ON_FIRST_PROXY
    Cint(50),  # FDB_TR_OPTION_CHECK_WRITES_ENABLE
    Cint(51),  # FDB_TR_OPTION_READ_YOUR_WRITES_DISABLE
    Cint(52),  # FDB_TR_OPTION_READ_AHEAD_DISABLE
    Cint(110), # FDB_TR_OPTION_DURABILITY_DATACENTER
    Cint(120), # FDB_TR_OPTION_DURABILITY_RISKY
    Cint(130), # FDB_TR_OPTION_DURABILITY_DEV_NULL_IS_WEB_SCALE
    Cint(200), # FDB_TR_OPTION_PRIORITY_SYSTEM_IMMEDIATE
    Cint(201), # FDB_TR_OPTION_PRIORITY_BATCH
    Cint(300), # FDB_TR_OPTION_INITIALIZE_NEW_DATABASE
    Cint(301), # FDB_TR_OPTION_ACCESS_SYSTEM_KEYS
    Cint(302), # FDB_TR_OPTION_READ_SYSTEM_KEYS
    Cint(400), # FDB_TR_OPTION_DEBUG_DUMP
    Cint(401), # FDB_TR_OPTION_DEBUG_RETRY_LOGGING
    Cint(402), # FDB_TR_OPTION_TRANSACTION_LOGGING_ENABLE
    Cint(500), # FDB_TR_OPTION_TIMEOUT
    Cint(501), # FDB_TR_OPTION_RETRY_LIMIT
    Cint(502), # FDB_TR_OPTION_MAX_RETRY_DELAY
    Cint(600), # FDB_TR_OPTION_SNAPSHOT_RYW_ENABLE
    Cint(601), # FDB_TR_OPTION_SNAPSHOT_RYW_DISABLE
    Cint(700), # FDB_TR_OPTION_LOCK_AWARE
    Cint(701), # FDB_TR_OPTION_USED_DURING_COMMIT_PROTECTION_DISABLE
    Cint(702)  # FDB_TR_OPTION_READ_LOCK_AWARE
)

const FDBStreamingMode = enum_fdb_streaming_mode(
    Cint(-2), # FDB_STREAMING_MODE_WANT_ALL
    Cint(-1), # FDB_STREAMING_MODE_ITERATOR
    Cint(0),  # FDB_STREAMING_MODE_EXACT
    Cint(1),  # FDB_STREAMING_MODE_SMALL
    Cint(2),  # FDB_STREAMING_MODE_MEDIUM
    Cint(3),  # FDB_STREAMING_MODE_LARGE
    Cint(4)   # FDB_STREAMING_MODE_SERIAL
)

const FDBMutationType = enum_fdb_mutation_type(
    Cint(2),  # FDB_MUTATION_TYPE_ADD
    Cint(6),  # FDB_MUTATION_TYPE_AND
    Cint(6),  # FDB_MUTATION_TYPE_BIT_AND
    Cint(7),  # FDB_MUTATION_TYPE_OR
    Cint(7),  # FDB_MUTATION_TYPE_BIT_OR
    Cint(8),  # FDB_MUTATION_TYPE_XOR
    Cint(8),  # FDB_MUTATION_TYPE_BIT_XOR
    Cint(12), # FDB_MUTATION_TYPE_MAX
    Cint(13), # FDB_MUTATION_TYPE_MIN
    Cint(14), # FDB_MUTATION_TYPE_SET_VERSIONSTAMPED_KEY
    Cint(15), # FDB_MUTATION_TYPE_SET_VERSIONSTAMPED_VALUE
    Cint(16), # FDB_MUTATION_TYPE_BYTE_MIN
    Cint(17)  # FDB_MUTATION_TYPE_BYTE_MAX
)

const FDBConflictRangeType = enum_fdb_conflict_range_type_t(
    Cint(0), # FDB_CONFLICT_RANGE_TYPE_READ
    Cint(1)  # FDB_CONFLICT_RANGE_TYPE_WRITE
)

const FDBErrorPredicate = enum_fdb_error_predicate_t(
    Cint(50000), # FDB_ERROR_PREDICATE_RETRYABLE
    Cint(50001), # FDB_ERROR_PREDICATE_MAYBE_COMMITTED
    Cint(50002)  # FDB_ERROR_PREDICATE_RETRYABLE_NOT_COMMITTED
)

const FDB_API_VERSION = Cint(510)
const fdb_error_t = Cint
const fdb_bool_t = Cint
const fdb_transaction_ptr_t = Ptr{Nothing}
const fdb_future_ptr_t = Ptr{Nothing}
const fdb_cluster_ptr_t = Ptr{Nothing}
const fdb_database_ptr_t = Ptr{Nothing}
const fdb_network_option_t = Cint
const fdb_cluster_option_t = Cint
const fdb_database_option_t = Cint
const fdb_transaction_option_t = Cint
const fdb_streaming_mode_t = Cint
const fdb_mutation_type_t = Cint
const fdb_conflict_range_type_t = Cint
const fdb_error_predicate_t = Cint

"""
A pointer to a function which takes FDBFuture* and void* and returns void.
"""
const FDBCallback = Ptr{Nothing}

"""
Represents a single key-value pair in the output of
fdb_future_get_keyvalue_array().

- key: A pointer to a key.
- key_length: The length of the key pointed to by key.
- value: A pointer to a value.
- value_length: The length of the value pointed to by value.
"""
struct fdb_kv_t
    key::Ptr{Nothing}
    key_length::Cint
    value::Ptr{Nothing}
    value_length::Cint
end


text:
# ---
# title: 869. Reordered Power of 2
# id: problem869
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/reordered-power-of-2/description/>
# hidden: true
# ---
#
# Starting with a positive integer `N`, we reorder the digits in any order
# (including the original order) such that the leading digit is not zero.
#
# Return `true` if and only if we can do this in a way such that the resulting
# number is a power of 2.
#
#
#
# **Example 1:**
#
#
#
#     Input: 1
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: 10
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: 16
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: 24
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: 46
#     Output: true
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.453305806680525 <= q <= 37.480198521751326)

                     add_NL_constraint(m, :(log(1 + exp(-0.5060529352779994 + 0.34605798298901425*log(1 + exp(-0.49627455608242954 + -0.46375151869975495*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + 0.4918746175115887*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.38157591143072933*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.5439511060803288*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.8504066842806357*log(1 + exp(-0.2321501407718678 + 0.5600091754803458*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.0034071169835434567*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.6765465160367254*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + 0.8220752699067027*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.462211279142271*log(1 + exp(-0.9807289687619161 + -0.5968595287514287*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.3897253830589995*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.5397807246416386*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.24558591797668194*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + -0.0428411186394575*log(1 + exp(0.7899732054652109 + -0.5363355140855166*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.6727460735415476*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.4129142282119118*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.4150508411017726*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + log(1 + exp(0.02992678336456356 + -0.8781954986928167*log(1 + exp(-0.49627455608242954 + -0.46375151869975495*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + 0.4918746175115887*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.38157591143072933*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.5439511060803288*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + -0.84297522225665*log(1 + exp(-0.2321501407718678 + 0.5600091754803458*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.0034071169835434567*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.6765465160367254*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + 0.8220752699067027*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.44937020244474857*log(1 + exp(-0.9807289687619161 + -0.5968595287514287*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.3897253830589995*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.5397807246416386*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.24558591797668194*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + -0.012040065223148932*log(1 + exp(0.7899732054652109 + -0.5363355140855166*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.6727460735415476*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.4129142282119118*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.4150508411017726*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + log(1 + exp(-0.25761034404046734 + -0.19911102915318102*log(1 + exp(-0.49627455608242954 + -0.46375151869975495*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + 0.4918746175115887*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.38157591143072933*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.5439511060803288*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.7457041629864394*log(1 + exp(-0.2321501407718678 + 0.5600091754803458*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.0034071169835434567*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.6765465160367254*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + 0.8220752699067027*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.7638178262546895*log(1 + exp(-0.9807289687619161 + -0.5968595287514287*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.3897253830589995*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.5397807246416386*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.24558591797668194*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.07889233901952375*log(1 + exp(0.7899732054652109 + -0.5363355140855166*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.6727460735415476*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.4129142282119118*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.4150508411017726*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + log(1 + exp(0.1153458350428731 + 0.9297971624965538*log(1 + exp(-0.49627455608242954 + -0.46375151869975495*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + 0.4918746175115887*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.38157591143072933*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.5439511060803288*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + -0.7017891371439648*log(1 + exp(-0.2321501407718678 + 0.5600091754803458*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.0034071169835434567*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.6765465160367254*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + 0.8220752699067027*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + -0.8093911002605614*log(1 + exp(-0.9807289687619161 + -0.5968595287514287*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.3897253830589995*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + 0.5397807246416386*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.24558591797668194*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))) + 0.43763294749909276*log(1 + exp(0.7899732054652109 + -0.5363355140855166*log(1 + exp(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5]))) + -0.6727460735415476*log(1 + exp(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5]))) + -0.4129142282119118*log(1 + exp(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5]))) + -0.4150508411017726*log(1 + exp(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module AcceleratedArrays

using SplitApplyCombine

using Base: @propagate_inbounds, Fix2, promote_op

export accelerate, accelerate!
export AcceleratedArray, AcceleratedVector, AcceleratedMatrix, MaybeVector, SingleVector
export AbstractIndex, AbstractUniqueIndex, HashIndex, UniqueHashIndex, SortIndex, UniqueSortIndex
export islessequal, isgreater, isgreaterequal

include("predicates.jl") # Add some predicates that are "missing" from Base
include("MaybeVector.jl")
include("SingleVector.jl")
include("Interval.jl")

include("AbstractIndex.jl")
include("AcceleratedArray.jl")
include("UniqueIndex.jl")
include("UniqueHashIndex.jl")
include("HashIndex.jl")
include("UniqueSortIndex.jl")
include("SortIndex.jl")

end # module

# TODO
#
# * Deal with <, ==, NaN, -0.0, missing, etc?
# * findfirst, findlast, findnext, findprev (requires stable sort?)
# * findmin, findmax (different behavior w.r.t NaN and missing?)
# * group, groupview, groupinds, groupreduce
# * leftgroupjoin, innerjoin, etc


text:
#=##############################################################################
# DESCRIPTION
    Windcraft vehicle kinematic maneuver definition.

# AUTHORSHIP
  * Author    : Judd Mehr and Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

function generate_maneuver_windcraft_kinematic(nrevs;
                                                 disp_plot=false,
                                                 includetail=true,
                                                 includewing=true,
                                                 includecontrols=true,
                                                 includerotors=true,
                                                 optargs...)

    omegamean = pi/180 * mean(omegastar.(range(0, 1, length=361); optargs...))

    ############################################################################
    # AIRCRAFT VELOCITY
    ############################################################################
    """
     Receives a nondimensional time between 0 and 1, and returns the
     vector of velocity of the vehicle at that instant.
    """
    function Vvehicle(t)

        theta = thetastar_periodic(t, nrevs; optargs...)*pi/180
        omega = omegastar_periodic(t, nrevs; optargs...)*pi/180
        scaling = omega/omegamean               # Scales velocity to a maximum value of 1

        Vcomp = [0.0, cos(theta), -sin(theta)]  # Counter-clockwise rotation
        return scaling*Vcomp
    end

    ############################################################################
    # AIRCRAFT ANGLES
    ############################################################################
    """
        Receives a nondimensional time between 0 and 1, and returns the angle
        (in degrees) of the vehicle.
        Returns: (angle_aircraft_x, angle_aircraft_y, angle_aircraft_z)
    """
    function anglevehicle(t)

        angle = [0.0, 0.0, thetastar_periodic(abs(t), nrevs; optargs...)]

        if t < 0.0
            return -angle
        else
            return angle
        end

    end

    ############################################################################
    # ANGLES OF CONTROL SURFACES
    ############################################################################

    """
        Receives a nondimensional time between 0 and 1, and returns the angles
        (in degrees) of the tilting system relative to the aircraft axis.
    """
    function angle_default(t)
        return [0, 0, 0]
    end

    if includetail == true
        angle = Tuple(angle_default for i in 1:6)  # Angle of each tilting system

    elseif includewing == true && includecontrols == true
        angle = Tuple(angle_default for i in 1:4)

    else
        angle = ()
    end

    ############################################################################
    # ROTORS RPM
    ############################################################################
    """
        Receives a nondimensional time between 0 and 1, and returns the RPM
        non-dimensionalized by the hover RPM.
    """
    function RPM_default(t)
        tsr = 4.0
        R = 1.0

        v = Vvehicle(t)
        vmag = sqrt(v[1]^2 + v[2]^2 + v[3]^2)
        omega = tsr*vmag/R
        rpm = omega*pi/30

        return rpm
    end

    if includerotors == true
        RPM = (RPM_default,)        # RPM of each rotor system
    else
        RPM = ()
    end

    ############################################################################
    # MANEUVER OBJECT
    ############################################################################
    maneuver = uns.KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)

    if disp_plot
        uns.plot(maneuver)
    end

    return maneuver

end

"""
angle theta output in degrees
"""
function thetastar(tstar; theta0=0.0, thetan=360.0, omegan=[4.0/9.0;20.0/27.0;4.0/9.0], tn=[0.0;0.5;1.0])
    #find index of max tn less than tstar and set to index variable m.
    if tstar==tn[end]
        m = length(tn)
    elseif tstar==tn[1]
        m = 1
    else

        m = maximum(find(x->x<tstar,tn))
    end
    # m = searchsortednearest(tn,tstar)
    #loop from 1 to m and sum areas.
    areas = 0.0
    den = 0.0
    for i = 2:length(tn)
        if i <= m
            areas += (tn[i]-tn[i-1])*(omegan[i]+omegan[i-1])/2.0
        end
        den += (tn[i]-tn[i-1])*(omegan[i]+omegan[i-1])/2.0
    end

    #add integral of final partial area up to tstar
    if m < length(tn) #if we haven't already integrated over everything
        integral = ( tstar - tn[m] ) * ( omegan[m] + (omegan[m+1]-omegan[m])/2.0 * (tstar-tn[m])/(tn[m+1]-tn[m]) )
    else
        integral = 0.0
    end

    #solve for scale factor, which
    sf = thetan / den
    #add theta0 (constant from integration)
    theta = sf * (areas + integral) + theta0

    return theta
end

function omegastar(tstar; h=1e-8, optargs...)
    return (thetastar(tstar+h; optargs...) - thetastar(tstar; optargs...))/h
end

"""
Transform time of thetastar and make it periodic to extend the angle output
outside the range [0, 360]
"""
function thetastar_periodic(t, nrevs; optargs...)
    tstar = (t*nrevs)%1    # Convert general time to time of one revolution
    return 360*floor(t*nrevs) + thetastar(tstar; optargs...)
end

function omegastar_periodic(tstar, nrevs; h=1e-8, optargs...)
    return (thetastar_periodic(tstar+h, nrevs; optargs...)
                - thetastar_periodic(tstar, nrevs; optargs...))/h
end


text:
"""
utility functions to randomly generate random graphical models with desired properties
"""

using Base.Iterators: enumerate
using Random: randperm

function rand_regular_graph(
    n::Integer,
    k::Integer,
    d::Integer;
    β::Real = 1.0,
    α::Real = 0.2,
    glassy::Bool = false,
)
    """
    random regular graph where the degree of factor node is fixed to be k, and the degree
    of variable node is `d` on average
    """

    # Singleton potential coefficients
    node_attr =
        [vcat(rand(Normal(0.0, 0.25), 1), rand(Normal(0.0, 1.0), k - 1)) for _ = 1:n]
    nfactor = Integer(ceil(n * d / k))
    all_indexk = collect(combinations(1:n, k))
    if nfactor > length(all_indexk)
        nfactor = length(all_indexk)
    end
    edge_index = all_indexk[randperm(length(all_indexk))[1:nfactor]]    # variables involved for each k-factor

    #  k-factor potential coefficients
    edge_attr = [rand(Uniform(0.02, 3.0), 1) * β for _ = 1:length(edge_index)]
    if glassy
        edge_attr .*= rand([-1.0, 1.0], length(edge_attr))
    end

    edge_type = [k for _ = 1:length(edge_index)]

    return DiscreteSpinGlass(edge_index, edge_attr .* β, edge_type, node_attr .* β, α)
end

## regular graphs

function randomfcspinglass(n::Integer;beta::Real=1.0)
    # generate a fully connected 2,3-spin glass distribution
    node_attr = [rand(Normal(0.0,1.0),1)*beta for _ in 1:n]
    edge_attr = [rand(Normal(0.0,1.0),1)*beta for _ in 1:(Int(n*(n-1)/2+n*(n-1)*(n-2)/6))]
    all_index2 = collect(combinations(1:n,2))
    all_index3 = collect(combinations(1:n,3))
    edge_index = vcat(all_index2,all_index3)
    edge_type = vcat([2 for _ in 1:length(all_index2)],[3 for _ in 1:length(all_index3)])

    return DiscreteSpinGlass(edge_index,edge_attr,edge_type,node_attr,0.0)
end

function randomfcferro(n::Integer;beta::Real=1.0)
    # generate a fully connected 2,3-spin glass distribution
    node_attr = [[0.0].*beta for _ in 1:n]
    edge_attr = [rand(Uniform(0.0,1.0),1)*beta for _ in 1:(Int(n*(n-1)/2+n*(n-1)*(n-2)/6))]
    all_index2 = collect(combinations(1:n,2))
    all_index3 = collect(combinations(1:n,3))
    edge_index = vcat(all_index2,all_index3)
    edge_type = vcat([2 for _ in 1:length(all_index2)],[3 for _ in 1:length(all_index3)])

    return DiscreteSpinGlass(edge_index,edge_attr,edge_type,node_attr,0.0)
end

function randomkfactorgraph(n::Integer,k::Integer,ρ::Real;β::Real=1.0)
    # generate a random k-factor graph without singleton potentials
    # ρ is the density of factors
    node_attr = [
        vcat(rand(Normal(0.0,0.25),1),rand(Normal(0.0,1.0),k-1)) for _ in 1:n
    ]
    all_indexk = collect(combinations(1:n,k))
    Cnk = length(all_indexk)
    Nfactor = Integer(ceil(ρ*length(all_indexk)))
    idx = randperm(length(all_indexk))[1:Nfactor]
    edge_index = all_indexk[idx]
    edge_attr = [rand(Normal(0.0,1.0),1)*β for _ in 1:length(edge_index)]
    edge_type = [k for _ in 1:length(edge_index)]

    return DiscreteSpinGlass(edge_index,edge_attr,edge_type,node_attr,0.0)
end

function randomregulargraph(n::Integer,k::Integer,d::Integer;β::Real = 1.0)
    # random regular graph where the degree of factor node is fixed to be k, and the degree
    # of variable node is also fixed to be d

    node_attr = [
        vcat(rand(Normal(0.0,0.25),1),rand(Normal(0.0,1.0),k-1)) for _ in 1:n
    ]
    nfactor = Integer(ceil(n*d/k))

    all_indexk = collect(combinations(1:n,k))
    edge_index = all_indexk[randperm(length(all_indexk))[1:nfactor]]

    edge_attr = [rand(Normal(0.0,1.0),1)*β for _ in 1:length(edge_index)]
    edge_type = [k for _ in 1:length(edge_index)]

    return DiscreteSpinGlass(edge_index,edge_attr,edge_type,node_attr,0.0)
end

function compute_count(channel::Integer, group::Integer)
    divide, remain = divrem(channel, group)
    out = zeros(Int,group)
    out[1:remain] .= divide +1
    out[remain+1:end] .= divide
    return out
end

function random_ws_graph(n::Integer,k::Real,p::Real)
    @assert k>=2 && k<=n
    edge_num = Int(round(k*n/2))
    count = compute_count(edge_num,n)
    idx = []
    # ring graph
    for i in 0:n-1
        source = repeat([i],count[i+1])
        target = collect(i+1:i+count[i+1]+1)
        target = [mod(t,n) for t in target]
        append!(idx, collect(zip(source,target)))
    end
    # rewire
    for i in 0:n-1
        u = i
        target = collect(i+1:i+count[i+1]+1)
        target = [mod(t,n) for t in target]
        for v in target
            if rand() < p
                w = rand(0:n-1)
                while w==u || (u,w) in idx
                    w = rand(0:n-1)
                    if length(findall(x->x[1]==u, idx)) >= n-1
                        break
                    end
                end
                deleteat!(idx, findall(x->x==1, (u,v)))
                push!(idx, (u,w))
            end
        end
    end
    # reshape into (nfactor, 2)
    nfactor = length(idx)
    factors = zeros(Int, nfactor,2)
    for i in 1:nfactor
        factors[i,:] = [idx[i][1], idx[i][2]]
    end
    return factors .+ 1
end

#### WS-style random graph and utilities
matchrow(a,B) = findfirst(i->all(j->a[j] == B[i,j],1:size(B,2)),1:size(B,1))

function random_ws_factor_structure(n::Integer, k::Real, p::Real)
    """
    k: in [2.0, floor((n-1)*(n-2)/6)]
    p: in [0,1]
    """
    f = floor(n*k/3) |> Int
    k̂ = floor(f/n) |> Int
    factors = zeros(Int, f,3)
    offsets = []
    for l in 1:n-1
        for i in 1:div(l,2)
            if l==2*i
                push!(offsets,[i,i])
            else
                if mod(n,2)==0 & l-2*i==1
                    push!(offsets,[i,l-i])
                else
                    push!(offsets,[i,l-i])
                    push!(offsets,[l-i,i])
                end
            end
        end
    end

    # 1. create "ring" factor graph
    counter = 0
    for i in 1:n, j in 1:k̂
        counter += 1
        off = offsets[j]
        factors[counter,:] = [mod(i-off[1], n), mod(i,n), mod(i+off[2],n)]
    end

    # 2. randomly connect remaining factors
    max_factor_num = length(combinations(1:n-1,2))
    for _ in 1:rem(f,n)
        counter += 1
        # find available center node
        found = false
        icenter = 0
        while !found
            i = rand(Categorical(n))
            if sum(factors[:,2]==i) < max_factor_num
                found = true
                icenter = i
            end
        end
        # find closest neibouring nodes to form a 3-factor
        for j in 1:length(offsets)
            off = offsets[j]
            candidate = [mod(icenter-off[1], n), mod(icenter,n), mod(icenter+off[2],n)]
            ioff = matchrow(candidate, factors)
            if ioff === nothing
                factors[counter,:] = candidate
                break
            end
        end
    end
    # keep each factor neighbor list sorted
    factors = mapslices(sort, factors, dims=2)

    # 3. randomly rewire each factor
    all_comb = collect(combinations(0:n-1,3))
    for i in 1:f
        if rand()< p
            all_remain = setdiff(all_comb, [factors[i,:] for i in 1:size(factors,1)])
            if length(all_remain) > 0
                candidate = rand(all_remain)
                factors[i,:] = candidate
            end
        end
    end

    return factors .+ 1
end

function random_ws_factor_graph(n::Integer, k2::Real, k3::Real, p::Real;γ::Real = 1.0,γ1=nothing, mode=:spinglass, degrees::Vector{<:Integer}=[3], singleton_dist=:uniform)
    if γ1 === nothing
        γ1 = γ
    end

    if singleton_dist == :uniform
        node_attr = [
                rand(Uniform(0.0, 1.0),3).*γ1  for _ in 1:n
            ]
    elseif singleton_dist == :randn
        node_attr = [
                rand(Normal(0.0,1.0),3).*γ1 for _ in 1:n
            ]
    elseif singleton_dist == :zero
        node_attr = [
                zeros(3).*γ1 for _ in 1:n
            ]
    end
    edge_index_tot = []
    edge_attr_tot = []
    edge_type_tot = []

    if 3 in degrees
        factors = random_ws_factor_structure(n,k3,p)
        edge_index = [factors[i,:] for i in 1:size(factors,1)]
        nfactor = length(edge_index)
        if mode==:spinglass
            edge_attr = [rand(Normal(0.0,1.0),1)*γ for _ in 1:length(edge_index)]
        elseif mode==:ferro
            edge_attr = [rand(Gamma(2.0),1)*γ for _ in 1:length(edge_index)]
        end
        edge_type = [3 for _ in 1:length(edge_index)]
        append!(edge_type_tot, edge_type)
        append!(edge_attr_tot, edge_attr)
        append!(edge_index_tot, edge_index)
    end
    if 2 in degrees
        factors = random_ws_graph(n,k2,p)
        edge_index = [factors[i,:] for i in 1:size(factors,1)]
        nfactor = length(edge_index)
        if mode==:spinglass
            edge_attr = [rand(Normal(0.0,1.0),1)*γ for _ in 1:length(edge_index)]
        elseif mode==:ferro
            edge_attr = [rand(Gamma(2.0),1)*γ for _ in 1:length(edge_index)]
        end
        edge_type = [2 for _ in 1:length(edge_index)]
        append!(edge_type_tot, edge_type)
        append!(edge_attr_tot, edge_attr)
        append!(edge_index_tot, edge_index)
    end

    return DiscreteSpinGlass(edge_index_tot,edge_attr_tot,edge_type_tot,node_attr,1.0)
end
function incident_matrix(d::Union{ThirdOrderExpFamily,DiscreteSpinGlass})
    A = zeros(Bool, d.nfactor, d.nnode)
    for i in 1:d.nfactor
        A[i, d.edge_index[i]] .= 1
    end
    return A
end


text:
module TestDoctest

using Documenter
using Test
using Julio

function test()
    doctest(Julio; manual = true)
end

end  # module


text:
struct ObsAdaptiveParticleFilter{P<:POMDP,S,R,RNG<:AbstractRNG} <: Updater
    pomdp::P
    resample::R
    max_frac_replaced::Float64
    rng::RNG
    _pm::Vector{S}
    _wm::Vector{Float64}
end

function ObsAdaptiveParticleFilter(p::POMDP, resample, max_frac_replaced, rng::AbstractRNG)
    S = statetype(p)
    return ObsAdaptiveParticleFilter(p, resample, max_frac_replaced, rng, S[], Float64[])
end

POMDPs.initialize_belief(up::ObsAdaptiveParticleFilter{S}, d::Any) where S = resample(up.resample, d, up.rng)
POMDPs.update(up::ObsAdaptiveParticleFilter, b, a, o) = update(up, resample(up.resample, b, up.rng), a, o)

function POMDPs.update(up::ObsAdaptiveParticleFilter, b::ParticleFilters.ParticleCollection, a, o)
    if n_particles(b) > 2*up.resample.n
        b = resample(up.resample, b, up.rng)
    end

    ps = particles(b)
    pm = up._pm
    wm = up._wm
    resize!(pm, 0)
    resize!(wm, 0)

    all_terminal = true
    for i in 1:n_particles(b)
        s = ps[i]
        if !isterminal(up.pomdp, s)
            all_terminal = false
            sp = generate_s(up.pomdp, s, a, up.rng)
            push!(pm, sp)
            od = observation(up.pomdp, s, a, sp)
            push!(wm, pdf(od, o))
        end
    end
    ws = sum(wm)
    if all_terminal || ws < eps(1.0/length(wm))
        # warn("All states in particle collection were terminal.")
        return initialize_belief(up, reset_distribution(up.pomdp, b, a, o))
    end

    pc = resample(up.resample, WeightedParticleBelief{statetype(up.pomdp)}(pm, wm, ws, nothing), up.rng)
    ps = particles(pc)

    mpw = max_possible_weight(up.pomdp, a, o)
    frac_replaced = up.max_frac_replaced*max(0.0, 1.0 - maximum(wm)/mpw)
    n_replaced = floor(Int, frac_replaced*length(ps))
    is = randperm(up.rng, length(ps))[1:n_replaced]
    for i in is
        ps[i] = new_particle(up.pomdp, b, a, o, up.rng)
    end
    return pc
end

ParticleFilters.weighted_posterior_belief(f::ObsAdaptiveParticleFilter) = WeightedParticleBelief(f._pm, f._wm)
ParticleFilters.predicted_belief(f::ObsAdaptiveParticleFilter) = ParticleCollection(f._pm)


text:
module TestTableTraits

using Test, DataFrames, IteratorInterfaceExtensions, TableTraits, DataValues

struct ColumnSource
end

TableTraits.supports_get_columns_copy_using_missing(::ColumnSource) = true
TableTraits.isiterabletable(::ColumnSource) = true

function TableTraits.get_columns_copy_using_missing(x::ColumnSource)
    return (a=[1,2,3], b=[4.,5.,6.], c=["A", "B", "C"])
end

@testset "TableTraits" begin
    df = DataFrame(a=[1,2,3], b=[1.,missing,3.])
    @test IteratorInterfaceExtensions.isiterable(df)
    @test TableTraits.isiterabletable(df)
    @test collect(IteratorInterfaceExtensions.getiterator(df)) ==
        [(a=1, b=DataValue(1.0)), (a=2, b=DataValue{Float64}()), (a=3, b=DataValue(3.0))]
    sdf = view(df, 1:2, :)
    @test IteratorInterfaceExtensions.isiterable(sdf)
    @test TableTraits.isiterabletable(sdf)
    @test collect(IteratorInterfaceExtensions.getiterator(sdf)) ==
        [(a=1, b=DataValue(1.0)), (a=2, b=DataValue{Float64}())]

    if VERSION ≥ v"1.0.0"
        df = DataFrame(ColumnSource())

        @test size(df)==(3,3)
        @test df.a==[1,2,3]
        @test df.b==[4.,5.,6.]
        @test df.c==["A", "B", "C"]
    end
end

end # module


text:
using GameZero
game_include("PointQuad.jl")

const WIDTH = 800
const HEIGHT = 800
const BACKGROUND = colorant"black"

mutable struct P
    x::Real
    y::Real
    angle::Real
end

const qt = PointQuad.QuadTree{P}(PointQuad.Sqr(WIDTH / 2, HEIGHT / 2, WIDTH), 4)

PointQuad.position(p::P) = (trunc(Int, p.x), trunc(Int, p.y))

function draw()
    agents = PointQuad.get_points(qt)
    PointQuad.clear!(qt)
    for p in agents
        if length(PointQuad.query(qt, PointQuad.Sqr(p.x, p.y, 12))) > 1
            draw(Circle(trunc(Int, p.x), trunc(Int, p.y), 0), colorant"green", fill = true)
        else
            draw(Circle(trunc(Int, p.x), trunc(Int, p.y), 0), colorant"white", fill = true)
        end
        new_x, new_y = p.x + cos(p.angle), p.y + sin(p.angle)
        if new_x < 5 || new_x >= WIDTH - 5 || new_y < 5 || new_y >= HEIGHT - 5
            new_x = min(WIDTH - 5, max(5, new_x))
            new_y = min(HEIGHT - 5, max(5, new_y))
            p.angle = rand() * 2 * π
        end
        p.x, p.y = new_x, new_y
        PointQuad.insert!(qt, p)
    end
    agents = nothing
    for sec in PointQuad.get_trees(qt)
        draw(Rect(trunc(Int, sec.boundary.x - sec.boundary.s / 2), trunc(Int, sec.boundary.y - sec.boundary.s / 2), trunc(Int, sec.boundary.s), trunc(Int, sec.boundary.s)), colorant"white", fill = false)
    end
end

function on_mouse_move(g::Game, pos)
    PointQuad.insert!(qt, P(pos[1], pos[2], rand() * 2 * π))
end