text:
export rand, profit!, evaluate!, project, initialize

## Constructors  ################################################

"Creation of a Gaussian copula with marginal distributions"
function GaussCopula(
  margs::Array{ContinuousUnivariateDistribution},
  Σ::Array{Real,2})
  if (size(Σ,1) != size(Σ,2)) | (size(Σ,1) != length(margs))
    error("dimensions are different")
  end
  GaussCopula(size(Σ,1), Σ, margs)
end

function PLInsurance(input::DataFrame,
                     i::Int, n_scen::Int, is_net::Bool)
  pl = PLInsurance(0, 0, 0, Array{Real}(undef, n_scen), 0, 0, 0)
  pl.ceded = (is_net ? input[i, :re_ceded] : 0)
  pl.premium = input[i, :premium] * (1 - pl.ceded)
  pl.costs =
    input[i, :cost_ratio] * pl.premium +
    input[i, :re_costs] * (input[i, :premium] - pl.premium)
  return pl
end

function PLTotal(n_scen)
  PLTotal(Array{Real}(undef, n_scen), 0, 0, 0)
end

function BuInvestments(id::Symbol,
                       name::AbstractString,
                       bu_ins::Array{BusinessUnit},
                       cost_ratio::Real,
                       invest_init::Real,
                       n_scen::Int)
  gross = PLInvestments(invest_init, 0,
                        zeros(Real, n_scen), 0, 0, 0)
  net = deepcopy(gross)
  for 𝑖𝑛𝑠 ∈ bu_ins
    gross.invest_bop += 𝑖𝑛𝑠.gross.premium
    net.invest_bop += 𝑖𝑛𝑠.net.premium
  end
  gross.costs = cost_ratio * gross.invest_bop
  net.costs = cost_ratio * net.invest_bop
  return BuInvestments(id, name, invest_init, gross, net)
end

## Interface ####################################################

"""
`rand(gc::GaussCopula, n::Int)`

`n` random samples from a Gaussian copula `gc`
"""
function rand(gc::GaussCopula, n::Int)
  u = zeros(Float64, gc.n, n)
  x = zeros(Float64, gc.n, n)
  z = rand(MvNormal(zeros(Float64, gc.n), gc.Σ), n)
  for 𝑖 ∈ 1:gc.n
    u[𝑖,:] = cdf.(Ref(Normal()), z[𝑖,:])
    x[𝑖,:] = quantile.(Ref(gc.marginals[𝑖]), u[𝑖,:])
  end
  return x'
end

"""
`profit!(pl::PLInsurance, r_distr::Vector{Float64}, s::Real)`

 Updates `pl.profit`, where `r_distr` is the loss distribution
 and `s` is the risk free interest rate
"""
function profit!(pl::PLInsurance,
                 r_distr::Vector{Float64},
                 s::Real)
  pl.profit =
    (1+s) * pl.premium .- (1 - pl.ceded) * r_distr .- pl.costs
  return pl
end

"""
`profit!(pl::PLInvestments, r_distr::Vector{Float64}, s::Real)`

 Updates `pl.profit`, where `r_distr` is the investment result
 and `s` is the risk free interest rate. Only the investment
 return above `s` counts as profit.
"""
function profit!(pl::PLInvestments,
                 r_distr::Vector{Float64},
                 s::Real)
  pl.profit = (r_distr .- s) * (pl.invest_bop) .- pl.costs
  return pl
end

"""
`profit!(pl::PLTotal, pl_bu::Array{ProfitLoss},
  cap_init::Real, costs_fixed::Real, s::Real)`

 Updates `pl.profit`, where `cap_init` is the initial capital,
 `pl_bu` are the profit loss accounts of the business units,
 `costs_fixed` are the fixed costs, and `s` is the risk free
 interest rate.
"""
function profit!(pl::PLTotal, pl_bu::Array{ProfitLoss},
                 cap_init::Real, costs_fixed::Real, s::Real)
  fill!(pl.profit, 0.0)
  for  𝑝𝑙𝑏𝑢 ∈ pl_bu
    pl.profit .+= 𝑝𝑙𝑏𝑢.profit
  end
  pl.profit .+= (s * cap_init - costs_fixed)
  return pl
end

"""
`evaluate!(pl::ProfitLoss, α::Real)`

 Calculate expected profit, economic capital (expected shortfall
   at safety level `α`), and RORAC for `pl`
"""
function evaluate!(pl::ProfitLoss, α::Real)
  pl.profit_mean = mean(pl.profit)
  pl.eco_cap = es(-pl.profit, α)
  pl.rorac =
    pl.eco_cap < eps() ? NaN : pl.profit_mean / pl.eco_cap
end

"""
`initialize(insurance_input::DataFrame, invest_input::DataFrame,
  tau_kendall::Matrix{Real}, n_scen::Int)`

 Set up an insurance company
"""
function initialize(insurance_input::DataFrame,
                    invest_input::DataFrame,
                    tau_kendall::Matrix{Real},
                    n_scen::Int)
  n_bu = nrow(insurance_input) + 1
  bu = Array{BusinessUnit}(undef, n_bu)
  distr = Array{ContinuousUnivariateDistribution}(undef, n_bu)

  for 𝑖 ∈ insurance_input[:ctr]
    bu[𝑖] =
      BuInsurance(insurance_input[𝑖, :id],
                  insurance_input[𝑖, :name],
                  PLInsurance(insurance_input, 𝑖, n_scen, false),
                  PLInsurance(insurance_input, 𝑖, n_scen, true))
    lognorm_sd = sqrt(log(1 + insurance_input[𝑖,:var_coeff]^2 ))
    lognorm_mean =
      log(insurance_input[𝑖, :loss_ratio] * bu[𝑖].gross.premium) -
      0.5lognorm_sd^2
    distr[𝑖] = LogNormal(lognorm_mean, lognorm_sd)
  end

  bu[invest_input[1, :ctr]] =
    BuInvestments(invest_input[1, :id],
                  invest_input[1, :name],
                  bu[1:(n_bu - 1)],
                  invest_input[1, :cost_ratio],
                  invest_input[1, :init],
                  n_scen)
  distr[invest_input[1, :ctr]] =
    Normal(invest_input[1, :mean], invest_input[1, :sd])
  gc = GaussCopula(distr, convert(Array{Real, 2},
                                  sin.(π/2 * tau_kendall)))
  return bu, gc
end

"""
`project(ins_input::DataFrame, inv_input::DataFrame,
  tau_kendall::Matrix{Real}, n_scen::Int, α::Real, s::Real
  costs_fixed::Real)`

 Set up an insurance company and project its results
"""
function project(ins_input::DataFrame,
                 inv_input::DataFrame,
                 tau_kendall::Matrix{Real},
                 n_scen::Int,
                 α::Real,
                 s::Real,
                 costs_fixed::Real)

  bu, gc = initialize(ins_input, inv_input, tau_kendall, n_scen)
  n_bu = length(bu)

  rand_distr = rand(gc, n_scen)
  for 𝑖 ∈ 1:n_bu
    for 𝑔𝑟𝑜𝑠𝑠_𝑛𝑒𝑡 ∈ [bu[𝑖].gross, bu[𝑖].net]
      profit!(𝑔𝑟𝑜𝑠𝑠_𝑛𝑒𝑡, rand_distr[:,𝑖], s)
      evaluate!(𝑔𝑟𝑜𝑠𝑠_𝑛𝑒𝑡,α)
    end
  end
  total = Total(PLTotal(n_scen), PLTotal(n_scen))
  ## get easier access for the following for loop
  bu_gross = ProfitLoss[bu[𝑖].gross for 𝑖 ∈ 1:n_bu]
  bu_net = ProfitLoss[bu[𝑖].net for 𝑖 ∈ 1:n_bu]
  for (𝑔𝑛, 𝑏𝑢𝑔𝑛) ∈ [(total.gross, bu_gross),
                     (total.net, bu_net)]
    profit!(𝑔𝑛, 𝑏𝑢𝑔𝑛, bu[end].init, costs_fixed, s)
    evaluate!(𝑔𝑛,α)
  end

  return bu, total
end


text:
deleteBlob(ctx, subscription_id::String, resource_group_name::String, diskuri::String; kwargs...) = deleteResource(ctx, subscription_id, resource_group_name, diskuri; kwargs...)

export deleteBlob

ext: jl


text:
module BrowseTables

export write_html_table, open_html_table, TableOptions, HTMLTable

using ArgCheck: @argcheck
import DefaultApplication
using DocStringExtensions: SIGNATURES, TYPEDEF
using Parameters: @unpack
import Tables

####
#### options and customizations
####

const HTMLHEADSTART = """
<!DOCTYPE html>
  <meta content="text/html; charset=utf-8" http-equiv="Content-Type">
  <meta content="width=device-width, initial-scale=1, shrink-to-fit=no" name="viewport">
"""

const DEFAULTCSSPATH = abspath(@__DIR__, "..", "assets", "BrowseTables.css")

Base.@kwdef struct TableOptions
    css_path::AbstractString = DEFAULTCSSPATH
    css_inline::Bool = true
end

####
#### low-level HTML construction
####

function write_escaped(io, str::AbstractString)
    for char in str
        if char == '&'
            write(io, "&")
        elseif char == '<'
            write(io, "<")
        elseif char == '>'
            write(io, ">")
        elseif char == '"'
            write(io, """)
        else
            write(io, char)
        end
    end
    nothing
end

escape_string(str::AbstractString) = (io = IOBuffer(); write_escaped(io, str); String(take!(io)))

"""
$(SIGNATURES)

Write tags `<tag>` and `</tag>` to `io` before and after a call to `f(io)`.

`attributes` are emitted in the opening tag. `bropen` and `brclose` determine if a newlines
are emitted after the opening and closing tags, respectively.
"""
function write_tags(f, io::IO, tag::AbstractString; bropen::Bool = false,
                    brclose::Bool = bropen, attributes::NamedTuple = NamedTuple())
    write(io, "<", tag)
    for (k, v) in pairs(attributes)
        write(io, " ", string(k), "=\"")
        write_escaped(io, v)
        write(io, "\"")
    end
    write(io, ">")
    bropen && println(io)
    f(io)
    write(io, "</", tag, ">")
    brclose && println(io)
    nothing
end

####
#### HTML writing API
####

"""
$(SIGNATURES)

Write the second argument to `io` using its HTML representation. This is used for formatting
string-like objectsinside cells, captions, etc, extend accordingly for custom types.
Defaults to whatever is emitted by `string`.

Extend `make_cell` for higher-level formatting using CSS.
"""
write_html(io::IO, object) = write_html(io, string(object))

write_html(io::IO, content::AbstractString) = write_escaped(io, content)

"""
$(TYPEDEF)

Wrapper for emitting contents which are already formatted as HTML.
"""
struct RawHTML{S <: String}
    html::S
end

write_html(io::IO, raw_html::RawHTML) = write(io, raw_html.html)

"""
$(SIGNATURES)

Content (written using `write_html`) and HTML attributes. HTML tag is specified by caller of
`write_html`.
"""
struct Cell{C, A <: NamedTuple}
    content::C
    attributes::A
end

Cell(content; kwargs...) = Cell(content, NamedTuple{keys(kwargs)}(values(kwargs)))

function write_html(io::IO, kind::Symbol, cell::Cell)
    @unpack content, attributes = cell
    write_tags(io -> write_html(io, content), io, string(kind); attributes = attributes)
end

"""
$(SIGNATURES)

Convert argument to `Cell`, which will be written using `write_html`.
"""
make_cell(::TableOptions, x) = Cell(x, NamedTuple())

make_cell(::TableOptions, x::Real) =
    Cell(x, isfinite(x) ? NamedTuple() : (class = "nonfinite", ))

make_cell(::TableOptions, str::AbstractString) = Cell(str; class = "alignleft")

make_cell(::TableOptions, x::Union{Missing,Nothing}) = Cell(repr(x); class = "likemissing")

####
#### table structure
####

function write_row(io, options::TableOptions, rowid, row)
    write_tags(io, "tr"; brclose = true) do io
        write_html(io, :td, Cell(string(rowid); class = "rowid"))
        for name in propertynames(row)
            write_html(io, :td, make_cell(options, getproperty(row, name)))
        end
    end
end

write_caption(io, ::Nothing) = nothing

write_caption(io, str) = write_tags(io -> write_html(io, str), io, "caption")

"""
$(SIGNATURES)

Write the schema returned by `Tables.schema` (for rows) to `io`, wrapped by `tag`.
"""
write_schema(io, tag, ::Nothing) = nothing

function write_schema(io, tag, sch::Tables.Schema)
    write_tags(io, tag; brclose = true) do io
        write_tags(io, "tr") do io
            write_html(io, :th, Cell("#"; class = "rowid")) # row ID
            for name in sch.names
                write_html(io, :th, Cell(string(name)))
            end
        end
    end
end

"""
$(SIGNATURES)

Write the stylesheet `inline` or via a `<link>`.
"""
function write_style(io::IO, path::AbstractString, inline::Bool)
    if inline
        write_tags(io -> write(io, read(path, String)), io, "style")
    else
        println(io, raw"<link href=\"", path, raw"\" rel=\"stylesheet\" type=\"text/css\">")
    end
    nothing
end

####
#### high-level API
####

require_table(table) =
    @argcheck Tables.istable(table) "The table should support the interface of Tables.jl."

"""
$(SIGNATURES)

Write a HTML representation of `table` to `io`.  A file path (string) may be given
instead of an `io`.

`title` and `caption` determine respective parts of the table. They will be escaped.

`options` can be used to specify table options, such as CSS.
"""
function write_html_table(io::IO, table;
                          title = "Table", caption = nothing,
                          options::TableOptions = TableOptions())
    require_table(table)
    write(io, HTMLHEADSTART) # contains an opening
    write_tags(io -> write_escaped(io, title), io, "title")
    write_style(io, options.css_path, options.css_inline)
    println(io, "")  # close manually, opened in HTMLHEADSTART
    write_tags(io, "body"; bropen = true) do io
        write_tags(io, "table"; bropen = true) do io
            write_caption(io, caption)
            rows = Tables.rows(table)
            sch = Tables.schema(rows)
            write_schema(io, "thead", sch)
            write_tags(io, "tbody"; bropen = true) do io
                for (id, row) in enumerate(rows)
                    write_row(io, options, id,  row)
                end
            end
            write_schema(io, "tfoot", sch)
        end
    end
    println(io, "")
    nothing
end

function write_html_table(filename::AbstractString, table; kwargs...)
    require_table(table)  # check it before creating a file
    open(io -> write_html_table(io, table; kwargs...), filename, "w")
    nothing
end

"""
$(SIGNATURES)

Write `table` to `filename` (temporary file is generated by default), then open it using the
default application (hopefully a browser). Keyword arguments are passed on to
[`write_html_table`](@ref).
"""
function open_html_table(table; filename = tempname() * ".html", kwargs...)
    write_html_table(filename, table; kwargs...)
    DefaultApplication.open(filename)
end

Base.@kwdef struct HTMLTable
    table
    title::AbstractString = "Table"
    caption = nothing
    options::TableOptions = TableOptions()
end

"""
    HTMLTable(table; title, caption, options)

Wrap `table` in an object that is showable as `text/html` in Julia's display system.
This is analogous to `Base.Docs.HTML`.

See [`write_html_table`](@ref) for how keyword arguments are interpreted.
"""
HTMLTable(table; kwargs...) = HTMLTable(; table = table, kwargs...)

Base.show(io::IO, ::MIME"text/html", html::HTMLTable) =
    write_html_table(io, html.table;
                     title = html.title,
                     caption = html.caption,
                     options = html.options)

end # module


text:
using Makie
using MixedModels
using MixedModelsMakie
using Random # we don't depend on exact PRNG vals, so no need for StableRNGs
using Test

using MixedModelsMakie: confint_table

@testset "There are no graphical tests" begin
    @test true
end

@testset "Simple linear regression" begin
    a, b = 1, 2
    x = collect(1:10)
    y = randn(MersenneTwister(42), 10) * 0.1
    @. y += a + b * x
    result = simplelinreg(x, y)
    @test result isa Tuple
    @test a ≈ result[1] atol=0.05
    @test b ≈ result[2] atol=0.05
end

m1 = fit(
    MixedModel,
    @formula(1000/reaction ~ 1 + days + (1+days|subj)),
    MixedModels.dataset(:sleepstudy),
)

@testset "confint_table" begin
    wald = confint_table(m1, 0.68)
    bsamp = parametricbootstrap(MersenneTwister(42), 1000, m1)
    boot = confint_table(bsamp, 0.68)

    @test wald.coefname == boot.coefname
    @test wald.estimate ≈ boot.estimate rtol=0.05
    @test wald.lower ≈ boot.lower rtol=0.05
    @test wald.upper ≈ boot.upper rtol=0.05
end

@testset "ranefinfo" begin
    reinfo = ranefinfo(m1)
    @test isone(length(reinfo))
    @test keys(reinfo) == (:subj, )
    re1 = only(reinfo)
    @test isa(re1, RanefInfo)
    @test re1.cnames == ["(Intercept)", "days"]
    @test first(re1.levels) == "S308"
    @test first(re1.ranef) ≈ -0.081830283935302 atol=1.0e-5
    @test first(re1.stddev) ≈ 0.14035486016644683 atol=1.0e-5
    f = Figure();
    @test f.content == Any[]
    caterpillar!(f, re1);
    @test length(f.content) == 2
    @test isa(first(f.content), Axis)
    @test isone(f.layout.nrows)
    @test f.layout.ncols == 2
    tbl = ranefinfotable(re1)
    @test keys(tbl) == (:name, :level, :cmode, :cstddev)
    @test length(tbl.cmode) == length(re1.cnames) * length(re1.levels)
end

@testset "shrinkageplot" begin
    f = shrinkageplot(m1)
    @test isone(length(f.content))
end

@testset "utilities" begin
    ppts = MixedModelsMakie.ppoints(64)
    @test length(ppts) == 64
    @test first(ppts) ≈ inv(128)
end


text:
module JungleHelperTorch

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Torch" Torch(x::Integer, y::Integer, flag::String="torch_flag", sprite::String="")

const placements = Ahorn.PlacementDict(
    "Torch (Jungle Helper)" => Ahorn.EntityPlacement(
        Torch
    )
)

function Ahorn.selection(entity::Torch)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 9, y - 22, 17, 22)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Torch, room::Maple.Room)
    Ahorn.drawSprite(ctx, "JungleHelper/TorchNight/TorchNightOff", 0, -11)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libim_/include/im.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const libim_ = "im") : (const libim_ = "libim")  # Name of im shared lib.

function imFileOpen(file_name::String, error::Ptr{Cint})
  ccall( (:imFileOpen, libim_), Ptr{imFile}, (Ptr{UInt8}, Ptr{Cint}), file_name, error)
end
function imFileOpenAs(file_name::String, format::String, error::Ptr{Cint})
  ccall( (:imFileOpenAs, libim_), Ptr{imFile}, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), file_name, format, error)
end
function imFileNew(file_name::String, format::String, error::Ptr{Cint})
  ccall( (:imFileNew, libim_), Ptr{imFile}, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), file_name, format, error)
end
function imFileClose(ifile::Ptr{imFile})
  ccall( (:imFileClose, libim_), Void, (Ptr{imFile},), ifile)
end
function imFileHandle(ifile::Ptr{imFile}, index::Integer)
  ccall( (:imFileHandle, libim_), Ptr{Void}, (Ptr{imFile}, Cint), ifile, index)
end
function imFileGetInfo(ifile::Ptr{imFile}, format::String, compression::String, image_count::Ptr{Cint})
  ccall( (:imFileGetInfo, libim_), Void, (Ptr{imFile}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), ifile, format, compression, image_count)
end
function imFileSetInfo(ifile::Ptr{imFile}, compression::String)
  ccall( (:imFileSetInfo, libim_), Void, (Ptr{imFile}, Ptr{UInt8}), ifile, compression)
end
function imFileSetAttribute(ifile::Ptr{imFile}, attrib::String, data_type::Integer, count::Integer, data::Ptr{Void})
  ccall( (:imFileSetAttribute, libim_), Void, (Ptr{imFile}, Ptr{UInt8}, Cint, Cint, Ptr{Void}), ifile, attrib, data_type, count, data)
end
function imFileGetAttribute(ifile::Ptr{imFile}, attrib::String, data_type::Ptr{Cint}, count::Ptr{Cint})
  ccall( (:imFileGetAttribute, libim_), Ptr{Void}, (Ptr{imFile}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}), ifile, attrib, data_type, count)
end
function imFileGetAttributeList(ifile::Ptr{imFile}, attrib::Ptr{Ptr{UInt8}}, attrib_count::Ptr{Cint})
  ccall( (:imFileGetAttributeList, libim_), Void, (Ptr{imFile}, Ptr{Ptr{UInt8}}, Ptr{Cint}), ifile, attrib, attrib_count)
end
function imFileGetPalette(ifile::Ptr{imFile}, palette::Ptr{Clong}, palette_count::Ptr{Cint})
  ccall( (:imFileGetPalette, libim_), Void, (Ptr{imFile}, Ptr{Clong}, Ptr{Cint}), ifile, palette, palette_count)
end
function imFileSetPalette(ifile::Ptr{imFile}, palette::Ptr{Clong}, palette_count::Integer)
  ccall( (:imFileSetPalette, libim_), Void, (Ptr{imFile}, Ptr{Clong}, Cint), ifile, palette, palette_count)
end
function imFileReadImageInfo(ifile::Ptr{imFile}, index::Integer, width::Ptr{Cint}, height::Ptr{Cint}, file_color_mode::Ptr{Cint}, file_data_type::Ptr{Cint})
  ccall( (:imFileReadImageInfo, libim_), Cint, (Ptr{imFile}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), ifile, index, width, height, file_color_mode, file_data_type)
end
function imFileWriteImageInfo(ifile::Ptr{imFile}, width::Integer, height::Integer, user_color_mode::Integer, user_data_type::Integer)
  ccall( (:imFileWriteImageInfo, libim_), Cint, (Ptr{imFile}, Cint, Cint, Cint, Cint), ifile, width, height, user_color_mode, user_data_type)
end
function imFileReadImageData(ifile::Ptr{imFile}, data::Ptr{Void}, convert2bitmap::Integer, color_mode_flags::Integer)
  ccall( (:imFileReadImageData, libim_), Cint, (Ptr{imFile}, Ptr{Void}, Cint, Cint), ifile, data, convert2bitmap, color_mode_flags)
end
function imFileWriteImageData(ifile::Ptr{imFile}, data::Ptr{Void})
  ccall( (:imFileWriteImageData, libim_), Cint, (Ptr{imFile}, Ptr{Void}), ifile, data)
end
function imFormatRegisterInternal()
  ccall( (:imFormatRegisterInternal, libim_), Void, (), )
end
function imFormatRemoveAll()
  ccall( (:imFormatRemoveAll, libim_), Void, (), )
end
function imFormatList(format_list::Ptr{Ptr{UInt8}}, format_count::Ptr{Cint})
  ccall( (:imFormatList, libim_), Void, (Ptr{Ptr{UInt8}}, Ptr{Cint}), format_list, format_count)
end
function imFormatInfo(format::String, desc::String, ext::String, can_sequence::Ptr{Cint})
  ccall( (:imFormatInfo, libim_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), format, desc, ext, can_sequence)
end
function imFormatInfoExtra(format::String, extra::String)
  ccall( (:imFormatInfoExtra, libim_), Cint, (Ptr{UInt8}, Ptr{UInt8}), format, extra)
end
function imFormatCompressions(format::String, comp::Ptr{Ptr{UInt8}}, comp_count::Ptr{Cint}, color_mode::Integer, data_type::Integer)
  ccall( (:imFormatCompressions, libim_), Cint, (Ptr{UInt8}, Ptr{Ptr{UInt8}}, Ptr{Cint}, Cint, Cint), format, comp, comp_count, color_mode, data_type)
end
function imFormatCanWriteImage(format::String, compression::String, color_mode::Integer, data_type::Integer)
  ccall( (:imFormatCanWriteImage, libim_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint, Cint), format, compression, color_mode, data_type)
end
function imEncodeColor(red::Cuchar, green::Cuchar, blue::Cuchar)
  ccall( (:imEncodeColor, libim_), Clong, (Cuchar, Cuchar, Cuchar), red, green, blue)
end
function imDecodeColor(red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, palette::Clong)
  ccall( (:imDecodeColor, libim_), Void, (Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Clong), red, green, blue, palette)
end
function imFileFormat(filename::String, format::Ptr{Cint})
  ccall( (:imFileFormat, libim_), Cint, (Ptr{UInt8}, Ptr{Cint}), filename, format)
end
function imImageInfo(filename::String, width::Ptr{Cint}, height::Ptr{Cint}, _type::Ptr{Cint}, palette_count::Ptr{Cint})
  ccall( (:imImageInfo, libim_), Cint, (Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), filename, width, height, _type, palette_count)
end
function imLoadRGB(filename::String, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar})
  ccall( (:imLoadRGB, libim_), Cint, (Ptr{UInt8}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}), filename, red, green, blue)
end
function imSaveRGB(width::Integer, height::Integer, format::Integer, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, filename::String)
  ccall( (:imSaveRGB, libim_), Cint, (Cint, Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{UInt8}), width, height, format, red, green, blue, filename)
end
function imLoadMap(filename::String, map::Ptr{Cuchar}, palette::Ptr{Clong})
  ccall( (:imLoadMap, libim_), Cint, (Ptr{UInt8}, Ptr{Cuchar}, Ptr{Clong}), filename, map, palette)
end
function imSaveMap(width::Integer, height::Integer, format::Integer, map::Ptr{Cuchar}, palette_count::Integer, palette::Ptr{Clong}, filename::String)
  ccall( (:imSaveMap, libim_), Cint, (Cint, Cint, Cint, Ptr{Cuchar}, Cint, Ptr{Clong}, Ptr{UInt8}), width, height, format, map, palette_count, palette, filename)
end
function imRGB2Map(width::Integer, height::Integer, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, map::Ptr{Cuchar}, palette_count::Integer, palette::Ptr{Clong})
  ccall( (:imRGB2Map, libim_), Void, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cint, Ptr{Clong}), width, height, red, green, blue, map, palette_count, palette)
end
function imMap2RGB(width::Integer, height::Integer, map::Ptr{Cuchar}, palette_count::Integer, colors::Ptr{Clong}, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar})
  ccall( (:imMap2RGB, libim_), Void, (Cint, Cint, Ptr{Cuchar}, Cint, Ptr{Clong}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}), width, height, map, palette_count, colors, red, green, blue)
end
function imRGB2Gray(width::Integer, height::Integer, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, map::Ptr{Cuchar}, grays::Ptr{Clong})
  ccall( (:imRGB2Gray, libim_), Void, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Clong}), width, height, red, green, blue, map, grays)
end
function imMap2Gray(width::Integer, height::Integer, map::Ptr{Cuchar}, palette_count::Integer, colors::Ptr{Clong}, grey_map::Ptr{Cuchar}, grays::Ptr{Clong})
  ccall( (:imMap2Gray, libim_), Void, (Cint, Cint, Ptr{Cuchar}, Cint, Ptr{Clong}, Ptr{Cuchar}, Ptr{Clong}), width, height, map, palette_count, colors, grey_map, grays)
end
function imResize(src_width::Integer, src_height::Integer, src_map::Ptr{Cuchar}, dst_width::Integer, dst_height::Integer, dst_map::Ptr{Cuchar})
  ccall( (:imResize, libim_), Void, (Cint, Cint, Ptr{Cuchar}, Cint, Cint, Ptr{Cuchar}), src_width, src_height, src_map, dst_width, dst_height, dst_map)
end
function imStretch(src_width::Integer, src_height::Integer, src_map::Ptr{Cuchar}, dst_width::Integer, dst_height::Integer, dst_map::Ptr{Cuchar})
  ccall( (:imStretch, libim_), Void, (Cint, Cint, Ptr{Cuchar}, Cint, Cint, Ptr{Cuchar}), src_width, src_height, src_map, dst_width, dst_height, dst_map)
end
function imRegisterCallback(cb::imCallback, cb_id::Integer, format::Integer)
  ccall( (:imRegisterCallback, libim_), Cint, (imCallback, Cint, Cint), cb, cb_id, format)
end


text:
VAR, REF, PAIR, CONST = 1, 2, 3, 4

println(VAR, REF, PAIR, CONST)

of_type = ["VAR", "REF", "PAIR", "CONST"]

mutable struct Var
    val::Any
    age::Int
    name::String
end

empty="#"
function ivar(age, name)
    Var(empty, age, name)
end

function type_of(x)
    if x isa Var
        if x.val == empty
            VAR
        else
            REF
        end
    elseif x isa Pair
        PAIR
    else
        CONST
    end
end

function bindVar(v,t,trail,ctr)
    v.val = t
    if v.age<ctr
        push!(trail,v)
    end
end

function deref(x)
    #println("DEREF<<<:",x,':',type_of(x))
    while true
        t = type_of(x)
        if t == REF
            x = x.val
        else
            #println("DEREF>>>:",x)
            @assert(x!=empty,"!!!!NOT EMPTY")
            return x, t
        end
    end
end

function unify(x,y,trail,ctr)
    ustack=[]
    push!(ustack,y)
    push!(ustack,x)
    while(!isempty(ustack))
        x1,t1=deref(pop!(ustack))
        x2,t2=deref(pop!(ustack))
        #println("### UNIFY ");pt(x1);pt(x2);println()
        if x1===x2
            continue
        elseif VAR==t1 && VAR==t2
          if x1.age>x2.age
              bindVar(x1,x2,trail,ctr)
          else
              bindVar(x2,x1,trail,ctr)
          end
        elseif VAR==t1
            bindVar(x1,x2,trail,ctr)
        elseif VAR==t2
            bindVar(x2,x1,trail,ctr)
        elseif t1!=t2
            return false
        elseif t1==CONST
            if x1!=x2
                return false
            end
        else
            @assert(PAIR==t1,"PAIR expected")
            a1,b1=x1
            a2,b2=x2
            push!(ustack,b2)
            push!(ustack,b1)
            push!(ustack,a2)
            push!(ustack,a1)
        end
    end
    return true
end

function unwind(trail,ttop)
  l=length(trail)
  #println("ENTER UNW ",l,',',ttop)
  for _ in ttop+1:l
      v=pop!(trail)
      #println("UNW",pt(v))
      v.val=nothing
  end
end

fuel=100

function step_(G,i,code,trail,varctr)
    ttop=length(trail)
    l=length(code)+1
    fail=(FAIL,nothing,ttop,i)
    global fuel;fuel-=1;if fuel<=0;return fail;end
    while i<l
        unwind(trail,ttop)
        clause,k=code[i]
        vs=[]
        i+=1
        #println("\n!!! ENTERING step: i=",i-1)
        for templ in clause
            (op,c1,c2,c3)=instr=activate(templ,k,vs)
            #println("ACTIVATE:",length(vs));ptrail(vs);println(vs);println("END\n")

            print("!!! step instr: i=",i-1,"--> ");pti(instr)

            if op=='u'
                old,oldt=deref(c3)
                if VAR==oldt
                    old.val=Pair(c1,c2)
                    #println("OLD---------");println(pt(old));println(">>>")
                    continue
                else
                    @assert(PAIR==oldt,pt(c3))

                    ok= unify(c1,old.first,trail,varctr.val)
                    if !ok
                        break
                    end
                    ok= unify(c2,old.second,trail,varctr.val)
                    if !ok
                        break
                    end
                end
                continue
            elseif op=='b'
                old,oldt=deref(c3)
                @assert(VAR==oldt,"VAR!!!")
                old.val=Pair(c1,c2)
                continue
            elseif op=='d'
                c1.val=G
                continue
            else
                @assert(op=='p',"EXPECTED: p")
                varctr.val+=k
                NewG,tg=deref(c1)
                println("p:----------:",pt(NewG),",TAG====",tg)
                if NewG=="true" && CONST == tg
                    return (DONE,G,ttop,i)
                else
                    #println("\nDO ERROR *** ",instr," NewG=",pt(NewG))
                    #ptrail(vs)
                    #println("@@@@@",length(vs),vs)
                    @assert(PAIR==tg,"!!!PAIR")
                    return (DO,(NewG,G),ttop,i)
                end
            end
        end
    end
    fail
end

FAIL,DONE,DO,UNDO=1,2,3,4

opn=["FAIL","DONE","DO","UNDO"]

function ensure_undo(G,todo,ttop,i,l)
    instr = (UNDO, G, ttop, i)
    push!(todo,instr)
end

function get_goal()
    answer=ivar(0,"ANSWER")
    Pair(answer,Pair("true","goal"))
end

function interp()
  global fuel;fuel-=1;if fuel<=0;return;end
  code=tload()
  l=length(code)+1
  todo=[];trail=[]
  goal=get_goal()
  answer=goal.first
  push!(todo,(DO,(goal,nothing),0,l+1))
  println("===========");ptodo(todo)
  varctr=ctr(1)
  while(!isempty(todo))
      action=pop!(todo)
      (op,G,ttop,i)=action
      #print("\n### interp: ");pta(action)
      #ptrail(trail)
      if DO==op
          NewG,OldG=G
          if i<l
              ensure_undo(OldG,todo,ttop,i,l)
          end
          r=step_(NewG,1,code,trail,varctr)

          push!(todo,r)
      elseif DONE==op
          print(">>>ANSWER:")
          pt(answer)
          println("")
          if i<l
              ensure_undo(G,todo,ttop,i,l)
          end
      elseif UNDO==op
          unwind(trail,ttop)
          if i<l
              push!(todo,step_(G,i,code,trail,varctr))
          end
      else # FAIL
          #println("FAILING: i=",i," l=",l," todo:",
            #length(todo)," trail=",length(trail), " END")
          #ptodo(todo)
      end
  end
end

# IO

function pt(z, size = 2000)
    if z==nothing;return z;end
    function pt0(a, k)
        if k <= 0
            return "<?>"
        else
            x, t = deref(a)
            if t == PAIR
                u, v = x
                k -= 1
                s1=pt0(u, k)
                s2=pt0(v, k)
                return string("(", s1, "->", s2, ")")
            elseif t == VAR
                #return string("_", string(objectid(x)))
                return string(x.name,"",x.age)
            else
                @assert(CONST==t,"should be CONST!")
                return string(x)
            end
        end
    end
    return pt0(z, size)
end

mutable struct ST
    d::Any
    a::Any
end

function new_st()
    ST(Dict(), [])
end

function add_cvar(w,st,ctr)
    if haskey(st.d, w)
        st.d[w]
    else
        ctr.val+=1
        push!(st.a, w)
        l = length(st.a)
        v = ivar(l,w)
        st.d[w] = ivar(l, w)
        v
    end
end

function vc(w,st,ctr)::Any
    if isuppercase(w[1])
        add_cvar(w,st,ctr)
    else
        w
    end
end

function to_instr(toks,st,ctr)
    op = toks[1]
    @assert(op!=nothing,"OP!!!")
    if (op == "d") || (op == "p")
        instr = (op[1], vc(toks[2],st,ctr), nothing, nothing)
    else
        instr = (op[1], vc(toks[2],st,ctr),
        vc(toks[3],st,ctr), vc(toks[4],st,ctr))
    end
end

mutable struct ctr
    val
end

function nv(t,k,vs)
    if t isa Var
        if(t.age<=length(vs))
            return vs[t.age]
        else
          v=ivar(t.age+k,t.name)
          push!(vs,v)
          return v
      end
    else
        t
    end
end

function activate(ts,k,vs)
   ts[1],nv(ts[2],k,vs),nv(ts[3],k,vs),nv(ts[4],k,vs)
end

function tload(fname = "../out/tnf_asm.txt")
    lines = readlines(fname)
    st=new_st()
    code=[]
    cls=[]
    count=ctr(0)
    for line in lines
        if length(line) < 2
            continue
        end
        #println("LINE:",line)
        toks = split(line)
        #for t in toks;println("TOK=",String(t));end;println()
        push!(cls,to_instr(toks,st,count))
        if toks[1]=="p"
            push!(code,(cls,count.val))
            st=new_st()
            cls=[]
            count=ctr(0)
        end
    end

    for i in 1:length(code)
        c,count=code[i]
        vs=[]
        for x in c
            #println('X',deref(x[1]))
            pti(x)
            a=activate(x,100,vs)

        end
        println("count: i=",i,"var count=",count) #;ptrail(vs)
    end

    code
end

function ptrail(xs)
    println("TRAIL:")
       for x in xs
           println(pt(x))
       end
       println("")
   end

   function ptodo(xs)
       println("TODO:")
          for x in xs
              pta(x)
          end
          println("")
      end

function pti(instr)
    #println("&&&_________",instr)
    s=" , "
    op,x1,x2,x3=instr
    #println("<",op,">::",pt(x1),s,pt(x2),s,pt(x3))
end

function pta(action)
    s=" , "
    act,G,ttop,i=action
    if act==DO
        println(opn[act],": ",pt(G[1]),s,pt(G[2]),s,ttop,i)
    else
        println(opn[act],": ",pt(G),s,s,ttop,i)
    end
end

########## tests ##########
function test()
    println("TEST")
    k=10
    trail=[]
    v=ivar(1,"V")
    u=ivar(2,"U")
    x=ivar(3,"X")
    y=ivar(4,"Y")
    bindVar(x,v,trail,k)
    bindVar(y,u,trail,k)
    println("DEREF: ",pt(x))
    p1=Pair(x,"a")
    p2=Pair("b",y)

    println("UNIF BEFORE: ",pt(p1)," ==== ",pt(p2))
    unify(p1,p2,trail,k)
    println("UNIF AFTER:  ",pt(p1)," ==== ",pt(p2))
    ptrail(trail)
    println(length(trail))
    unwind(trail,0)
    println("UNWINDED")
    ptrail(trail)

end

@assert(1==1,"we got  problem")

clearconsole()
#test()
interp()

max_stars_repo_name: ptarau/LogicTransformers

max_issues_repo_name: ptarau/LogicTransformers

max_forks_repo_name: ptarau/LogicTransformers


text:
function a_res(model::ogSAFTFamily,z,v,T)
    return a_seg(model,z,v,T) + a_chain(model,z,v,T) + a_assoc(model,z,v,T)
end

function a_seg(model::ogSAFTFamily,z,v,T)
    m = model.params.segment
    x = z/sum(z[i] for i in model.components)
    m̄ = sum(x[i]*m[i] for i in model.components)

    return m̄*(a_hs(model,z,v,T)+a_disp(model,z,v,T))
end

function a_chain(model::ogSAFTFamily,z,v,T)
    x = z/sum(z[i] for i in model.components)
    m = model.params.segment
    return sum(x[i]*(1-m[i])*log(g_hsij(model,z,v,T, i, i)) for i in model.components)
end

function d(model::ogSAFTFamily,z,v,T, component)
    ϵ = model.params.epsilon[union(component,component)]
    σ = model.params.sigma[union(component,component)]
    m = model.params.segment[component]
    fm = 0.0010477+0.025337*(m-1)/m
    f = (1+0.2977T/ϵ)/(1+0.33163T/ϵ+fm*(T/ϵ)^2)
    return σ * f
end

function dx(model::ogSAFTFamily,z,v,T)
    x = z/sum(z[i] for i in model.components)
    m = model.params.segment
    σ = model.params.sigma
    ϵ = model.params.epsilon

    mx = sum(x[i]*m[i] for i in model.components)
    σx = (sum(x[i]*x[j]*m[i]*m[j]*σ[union(i,j)]^3 for i in model.components for j in model.components)/mx^2)^(1/3)
    ϵx = (sum(x[i]*x[j]*m[i]*m[j]*σ[union(i,j)]^3*ϵ[union(i,j)] for i in model.components for j in model.components)/mx^2)/σx^3

    fm = 0.0010477+0.025337*(mx-1)/mx
    f = (1+0.2977T/ϵx)/(1+0.33163T/ϵx+fm*(T/ϵx)^2)
    return σx * f
end

function ζn(model::ogSAFTFamily,z,v,T, n)
    x = z/sum(z[i] for i in model.components)
    m = model.params.segment
    return N_A*sum(z[i] for i in model.components)*π/6/v * sum(x[i]*m[i]*d(model,z,v,T, i)^n for i in model.components)
end

function η(model::ogSAFTFamily,z,v,T)
    m = model.params.segment
    x = z/sum(z[i] for i in model.components)
    m̄ = sum(x[i]*m[i] for i in model.components)
    return N_A*sum(z[i] for i in model.components)*π/6/v*dx(model,z,v,T)^3*m̄
end

function g_hsij(model::ogSAFTFamily,z,v,T, i, j)
    di = d(model,z,v,T, i)
    dj = d(model,z,v,T, j)
    ζ2 = ζn(model,z,v,T, 2)
    ζ3 = ζn(model,z,v,T, 3)
    return 1/(1-ζ3) + di*dj/(di+dj)*3ζ2/(1-ζ3)^2 + (di*dj/(di+dj))^2*2ζ2^2/(1-ζ3)^3
end

function a_hs(model::ogSAFTFamily,z,v,T)
    ηx = η(model,z,v,T)
    return (4ηx-3ηx^2)/(1-ηx)^2
end

function a_disp(model::ogSAFTFamily,z,v,T)
    m = model.params.segment
    σ = model.params.sigma
    ϵ = model.params.epsilon
    x = z/sum(z[i] for i in model.components)
    ϵx = sum(x[i]*x[j]*m[i]*m[j]*σ[union(i,j)]^3*ϵ[union(i,j)] for i in model.components for j in model.components)/sum(x[i]*x[j]*m[i]*m[j]*σ[union(i,j)]^3 for i in model.components for j in model.components)
    ηx = η(model,z,v,T)
    ρR = (6/sqrt(2)/π)*ηx
    TR = T/ϵx
    a_seg1 = ρR*(-8.5959-4.5424ρR-2.1268ρR^2+10.285ρR^3)
    a_seg2 = ρR*(-1.9075+9.9724ρR-22.216ρR^2+15.904ρR^3)
    return 1/TR*(a_seg1+a_seg2/TR)
end

## This is an attempt to make Twu et al.'s segment term; does not work yet
# function a_seg(model::ogSAFTFamily,z,v,T)
#     Bo = [1.31024,-3.80636,-2.37238,-0.798872,0.198761,1.47014,-0.786367,2.19465,5.75429,6.7822,-9.94904,-15.6162,86.643,18.527,9.04755,8.68282]
#     Ba = [3.79621,-6.14518,-1.84061,-2.77584,-0.420751,-5.66128,19.2144,-3.33443,33.0305,-5.90766,9.55619,-197.883,-61.2535,77.1802,-6.57983,0.0]
#     ω = 0.199
#     A = []
#     for i in 1:16
#         append!(A,Bo[i]+ω*Ba[i])
#     end
#     m = model.params.segment
#     σ = model.params.sigma
#     ϵ = model.params.epsilon
#     x = z/sum(z[i] for i in model.components)
#     mx = sum(x[i]*m[i] for i in model.components)
#     σx = (sum(x[i]*x[j]*m[i]*m[j]*σ[union(i,j)]^3 for i in model.components for j in model.components)/mx^2)^(1/3)
#     ϵx = (sum(x[i]*x[j]*m[i]*m[j]*σ[union(i,j)]^3*ϵ[union(i,j)] for i in model.components for j in model.components)/mx^2)/σx^3
#
#     ρ  = N_A*sum(z[i] for i in model.components)/v
#     ρR = ρ*σx^3
#     TR = T/ϵx
#
#     u_res = (A[2]/TR+2A[3]/TR^2+3A[4]/TR^3+5A[5]/TR^5)*ρR+1/2*A[7]/TR*ρR^2+
#             1/(2*A[16])*(3A[9]/TR^3+4A[10]/TR^4+5A[11]/TR^5)*(1-exp(-A[16]*ρR^2))+
#             1/(2*A[16]^2)*(3A[12]/TR^3+4A[13]/TR^4+5A[14]/TR^5)*(1-(1+A[16]*ρR^2)*exp(-A[16]*ρR^2))+
#             1/5*A[15]/TR*ρR^5
#     s_res = -log(ρ*R*T)-(A[1]-A[3]/TR^2-2A[4]/TR^3-4A[5]/TR^5)*ρR-1/2*A[6]*ρR^2-1/3*A[8]*ρR^3+
#             1/(2*A[16])*(2A[9]/TR^3+3A[10]/TR^4+4A[11]/TR^5)*(1-exp(-A[16]*ρR^2))+
#             1/(2*A[16]^2)*(2A[12]/TR^3+3A[13]/TR^4+4A[14]/TR^5)*(1-(1+A[16]*ρR^2)*exp(-A[16]*ρR^2))
#     a_res = u_res-s_res
#     println(A)
#     println(u_res)
#     println(s_res)
#     return mx*(a_res)
# end

function a_assoc(model::ogSAFTFamily, z, v, T)
    x = z/sum(z[i] for i in model.components)
    X_iA = X_assoc(model,z,v,T)
    return sum(x[i]*sum(log(X_iA[i,a])-X_iA[i,a]/2 + model.params.n_sites[i][a]/2 for a in keys(model.params.n_sites[i])) for i in model.components)
end

function X_assoc(model::ogSAFTFamily, z, v, T)
    x = z/sum(z[i] for i in model.components)
    ρ = N_A*sum(z[i] for i in model.components)/v
    X_iA = Dict()
    X_iA_old = Dict()
    tol = 1.
    iter = 1
    while tol > 1e-12 && iter < 100
        for i in model.components
            for a in keys(model.params.n_sites[i])
                A = 0.
                for j in model.components
                    B = 0
                    for b in keys(model.params.n_sites[j])
                        if haskey(model.params.epsilon_assoc,Set([(i,a),(j,b)]))
                            if iter!=1
                                B+=X_iA_old[j,b]*Δ(model,z,v,T,i,j,a,b)
                            else
                                B+=Δ(model,z,v,T,i,j,a,b)
                            end
                        end
                    end
                    A += ρ*x[j]*B
                end
                if iter == 1
                    X_iA[i,a] =0.5+0.5*(1+A)^-1
                else
                    X_iA[i,a] =0.5*X_iA_old[i,a]+0.5*(1+A)^-1
                end
            end
        end
        if iter == 1
            tol = sqrt(sum(sum((1. -X_iA[i,a])^2 for a in keys(model.params.n_sites[i])) for i in model.components))
        else
            tol = sqrt(sum(sum((X_iA_old[i,a] -X_iA[i,a])^2 for a in keys(model.params.n_sites[i])) for i in model.components))
        end
        X_iA_old = deepcopy(X_iA)
        iter += 1
    end

    return X_iA
end

function Δ(model::ogSAFTFamily, z, v, T, i, j, a, b)
    ϵ_assoc = model.params.epsilon_assoc[Set([(i,a),(j,b)])]
    κ = model.params.bond_vol[Set([(i,a),(j,b)])]
    g = g_hsij(model,z,v,T,i,j)
    return (d(model, z, v, T, i)+d(model, z, v, T, j))^3/2^3*g*(exp(ϵ_assoc/T)-1)*κ
end


text:
import ..Properties: @property
import ..Vortex

@property begin
    signature = induce_acc(targ::Target, targvel::Target, src::Source, srcvel::Source)
    preallocator = allocate_acc
    stype = ComplexF64
end

function induce_acc(z::Complex{T}, ż::Complex{T},
                    p::Vortex.Point, point_vel::Complex{S}) where {T,S}
    -circulation(p)*Points.cauchy_kernel((z - p.z)^2)*conj(ż - point_vel)
end

function induce_acc(z::Complex{T}, ż::Complex{T},
                           b::Vortex.Blob, blob_vel::Complex{S}) where {T,S}
    induce_acc(z, ż,
               Vortex.Point(Elements.position(b), circulation(b)),
               blob_vel)
end

function surface_pressure_inst(p::Plate{T}, ṗ, ambient_sys, z_new, t, Δt, lesp, tesp, ss::AbstractArray{S}) where {T, S <: Real}
    @get p (L, C, α, dchebt!)
    @get ṗ (ċ, c̈, α̇ , α̈ )

    press = zeros(T, length(ss))

    # Get Ċ from movement of existing vortex blobs (without vortex shedding)
    enforce_no_flow_through!(p, ṗ, ambient_sys, t)

    srcvel = self_induce_velocity(ambient_sys, t)
    induce_velocity!(srcvel, ambient_sys, p, t)

    targvel = fill(ċ, length(p))
    Ċ = zeros(Complex{T},length(targvel))

    induce_acc!(Ċ, p.zs, targvel, ambient_sys, srcvel)

    z₊, z₋ = z_new
    point₊ = Vortex.Point(z₊, 1.0)
    point₋ = Vortex.Point(z₋, 1.0)

    Γ₊, Γ₋, ∂C₊, ∂C₋ = vorticity_flux!(p, point₊, point₋, t, lesp, tesp)

    n̂ = exp(-im*α)
    rmul!(Ċ, n̂)
    dchebt! * Ċ

    @. Ċ += (∂C₊ + ∂C₋)/Δt - im*α̇*C

    Ȧ = MappedVector(_chebyshev_coefficient, Ċ, 1)

    # Plate is subject to translation and rotation
    #B₀ = n⋅ẋc and B₁ = L/2α with n = i*exp(iα)
    Ḃ₀ = real(-α̇ *exp(-im*α)*ċ  + -im*exp(-im*α)*c̈)
    Ḃ₁ = L/2*α̈

    # Γ̇ = [Γ₋/Δt + _bound_circulation(Ȧ, Ḃ₀, Ḃ₁, L, -(Γ₊ + Γ₋)/Δt, s) for s in p.ss]
    map!(s-> Γ₋/Δt + _bound_circulation(Ȧ, Ḃ₀, Ḃ₁, L, -(Γ₊ + Γ₋)/Δt, s), press, ss)

    # strength(p) .* (Chebyshev.firstkind(real.(C), ss) .- tangent(ċ, α)) .+ Γ̇
    press .+=  strength(p, ss) .* (Chebyshev.firstkind(real.(C), ss) .- tangent(ċ, α))
    return press
end

surface_pressure_inst(p::Plate{T}, ṗ, ambient_sys, z_new, t, Δt, lesp, tesp) where {T} = surface_pressure_inst(p, ṗ, ambient_sys, z_new, t, Δt, lesp, tesp, p.ss)


text:
using DataFrames
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/TotalAdaptationCosts.jl")

    add_comp!(m, TotalAdaptationCosts)

    set_param!(m, :TotalAdaptationCosts, :pop_population, readpagedata(m, "test/validationdata/$valdir/pop_population.csv"))
    set_param!(m, :TotalAdaptationCosts, :ac_adaptationcosts_economic, readpagedata(m, "test/validationdata/ac_adaptationcosts_economic.csv"))
    set_param!(m, :TotalAdaptationCosts, :ac_adaptationcosts_noneconomic, readpagedata(m, "test/validationdata/ac_adaptationcosts_noneconomic.csv"))
    set_param!(m, :TotalAdaptationCosts, :ac_adaptationcosts_sealevelrise, readpagedata(m, "test/validationdata/$valdir/ac_adaptationcosts_sealevelrise.csv"))

    run(m)

    # Generated data
    adapt_cost = m[:TotalAdaptationCosts, :act_adaptationcosts_total]
    adapt_cost_per_cap = m[:TotalAdaptationCosts, :act_percap_adaptationcosts]

    # Recorded data
    cost_compare = readpagedata(m, "test/validationdata/$valdir/act_adaptationcosts_tot.csv")
    cost_cap_compare = readpagedata(m, "test/validationdata/$valdir/act_percap_adaptationcosts.csv")

    @test adapt_cost ≈ cost_compare rtol = 1e-3
    @test adapt_cost_per_cap ≈ cost_cap_compare rtol = 1e-6
end


text:
# It is possible to show that the square root of two can be expressed as an
# infinite continued fraction.
#
# √2=1+1/(2+1/(2+1/(2+...)))
#
# By expanding this for the first four iterations, we get:
#
# 1+1/2 = 3/2 = 1.5
# 1+1/(2+1/2) = 7/5 = 1.4
# 1+1/(2+1/(2+1/2)) = 17/12 =1.41666...
# 1+1/(2+1/(2+1/(2+1/2))) = 41/29 = 1.41379...
#
# The next three expansions are 99/70, 239/169, and 577/408, but the eighth
# expansion, 1393/985, is the first example where the number of digits in the
# numerator exceeds the number of digits in the denominator.
#
# In the first one-thousand expansions, how many fractions contain a numerator
# with more digits than the denominator?

using ProjectEulerSolutions

# Using Julia's built-in rational numbers only works up to n=47 or so, when we
# begin to run into over-runs.  Really this is a pattern matching game.  The
# numerator for iteration i (n_i) and the denominator for iteration i (d_i)
# follow a pattern: d_i = n_{i-1} + d_{i-1}, n_i = d_i + d_{i-1}
function p057solution(n_max::Integer=100)::Integer

    num = BigInt(3)
    denom = BigInt(2)
    count = 0
    for n in 2:n_max
        denom, num = num + denom, num + 2*denom
        if ndigits(num) > ndigits(denom)
            count += 1
        end
    end
    return count
end

p057 = Problems.Problem(p057solution)

Problems.benchmark(p057, 1_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# ================================================================
# Sparse reach set
# ================================================================

"""
    SparseReachSet{N, ST<:LazySet{N}, D} <: AbstractReachSet{N}

Type that wraps a reach-set using a `LazySet` as underlying representation,
together with a tuple of variables associated to this reach-set.

### Fields

- `X`        -- set
- `Δt`       -- time interval
- `vars` -- tuple of variable indices represented by the set `X`

### Notes

A `SparseReachSet` is a struct representing (an approximation of) the reachable
states for a given time interval. The type of the representation is `ST`, which
may be any subtype of `LazySet` (ideally, concrete). Moreover, this type also
stores information about the variables (also named coordinates, or by abuse of
notation, *dimensions*) corresponding to the set `X`.

For instance in the ambient space `n=5`, one may have a `SparseReachSet` whose
variables tuple is `vars = (4, 5, 6)`, i.e. representing a three-dimensional
projection of the full-dimensional reach-set. In consequence, the dimension of
`X` doesn't match the length of `vars`, in general

In this type, the parameter `N` represents the numerical type of the `LazySet`
(typically, `Float64`), the type `ST` represents the set representation used,
and `D` denotes the dimension of this sparse reach set. Note that, in contrast
to `ReachSet`, for `SparseReachSet` the number of dimensions is part of the type
information.
"""
struct SparseReachSet{N, ST<:LazySet{N}, D} <: AbstractLazyReachSet{N}
    X::ST
    Δt::IA.Interval{Float64}
    vars::NTuple{D, Int}
    # TODO: inner constructor that the dimension of vars matches that of X ?
end

# interface functions
set(R::SparseReachSet) = R.X
setrep(::SparseReachSet{N, ST}) where {N, ST<:LazySet{N}} = ST
setrep(::Type{<:SparseReachSet{N, ST}}) where {N, ST<:LazySet{N}} = ST
tstart(R::SparseReachSet) = inf(R.Δt)
tend(R::SparseReachSet) = sup(R.Δt)
tspan(R::SparseReachSet) = R.Δt
dim(R::SparseReachSet{N, ST, D}) where {N, ST<:LazySet{N}, D} = D
vars(R::SparseReachSet) = R.vars

# constructor from vector of dimensions
function SparseReachSet(X::ST, Δt::IA.Interval{Float64},
                        vars::AbstractVector) where {N, ST<:LazySet{N}}
    SparseReachSet(X, Δt, Tuple(vars))
end

function shift(R::SparseReachSet, t0::Number)
    return SparseReachSet(set(R), tspan(R) + t0, vars(R))
end

function reconstruct(R::SparseReachSet, Y::LazySet)
    return SparseReachSet(Y, tspan(R), vars(R))
end

# constructor with a time point
function SparseReachSet(X::ST, t::Real, vars::AbstractVector) where {N, ST<:LazySet{N}}
    return SparseReachSet(X, interval(t), vars)
end


text:
# =================================================
# DummySet: a set to test default methods
# =================================================
D = DummySet()
@test volume(D) == nothing
@test surface_area(D) == nothing

# =================================================
# Singleton: a singleton is a set with one element
# =================================================
s = Singleton([1.0, 2.0, 3.0])

# these are not-implemented functions
@test volume(s) == 0.
@test surface_area(s) == 0.

# =================================================
# Ball: a ball in n-dimensional space
# =================================================
B = Ball(zeros(3), 1.0)

@test volume(B) ≈ 4/3 * pi
@test surface_area(B) == 4 * pi


text:
"""
    plot_events(h)

Plot an event history with its marks if those are of a real subtype.
"""
function plot_events(h::TemporalHistory{M}) where {M<:Real}
    scatter(
        h.times,
        h.marks,
        xlim = (h.tmin, h.tmax),
        title = "Event history",
        xlabel = "Time",
        ylabel = "Marks",
        label = nothing,
    )
end

"""
    plot_intensity(pp, h; npoints)

Plot the conditional intensity function of a temporal point process along a given history.
"""
function plot_intensity(
    pp::TemporalPointProcess{M},
    h::TemporalHistory{M};
    npoints = 100,
) where {M}
    curve_times = h.tmin:(duration(h)/npoints):h.tmax
    curve_vals = [ground_intensity(pp, h, t) for t in curve_times]
    history_times = h.times
    history_vals = [ground_intensity(pp, h, t) for t in history_times]
    plot(curve_times, curve_vals, label = "intensity")
    scatter!(history_times, history_vals, label = "events")
    plot!(
        xlim = (h.tmin, h.tmax),
        xlabel = "Time",
        ylabel = "Ground intensity",
        title = "Intensity evolution",
    )
end

"""
    qqplot_interevent_times(h)

Compare the distribution of inter-event times in a given history to the exponential distribution.
"""
function qqplot_interevent_times(h::TemporalHistory)
    plot(
        qqplot(Dists.Exponential, diff(vcat(0., event_times(h)))),
        xlabel = "Quantiles of exponential distribution",
        ylabel = "Quantiles of inter-event times",
        title = "Residual analysis",
    )
end


text:
module MLJXGBoostInterface

export XGBoostRegressor, XGBoostClassifier, XGBoostCount

import MLJModelInterface
import MLJModelInterface: Table, Continuous, Count, Finite, OrderedFactor,
    Multiclass

const PKG="MLJXGBoostInterface"
const MMI = MLJModelInterface

import Tables: schema

import XGBoost

# helper for feature importances:
# XGBoost used "f
function named_importance(fi::XGBoost.FeatureImportance, features)

    new_fname = features[parse(Int, fi.fname[2:end]) + 1] |> string
    return XGBoost.FeatureImportance(new_fname, fi.gain, fi.cover, fi.freq)
end

# TODO: Why do we need this?
generate_seed() = mod(round(Int, time()*1e8), 10000)

# helper to preprocess hyper-parameters:
function kwargs(model, silent, seed, objective, eval_metric)

    kwargs = (booster = model.booster
              , silent = silent
              , disable_default_eval_metric = model.disable_default_eval_metric
              , eta = model.eta
              , gamma = model.gamma
              , max_depth = model.max_depth
              , min_child_weight = model.min_child_weight
              , max_delta_step = model.max_delta_step
              , subsample = model.subsample
              , colsample_bytree = model.colsample_bytree
              , colsample_bylevel = model.colsample_bylevel
              , lambda = model.lambda
              , alpha = model.alpha
              , tree_method = model.tree_method
              , sketch_eps = model.sketch_eps
              , scale_pos_weight = model.scale_pos_weight
              , refresh_leaf = model.refresh_leaf
              , process_type = model.process_type
              , grow_policy = model.grow_policy
              , max_leaves = model.max_leaves
              , max_bin = model.max_bin
              , predictor = model.predictor
              , sample_type = model.sample_type
              , normalize_type = model.normalize_type
              , rate_drop = model.rate_drop
              , one_drop = model.one_drop
              , skip_drop = model.skip_drop
              , feature_selector = model.feature_selector
              , top_k = model.top_k
              , tweedie_variance_power = model.tweedie_variance_power
              , objective = objective
              , base_score = model.base_score
              , eval_metric=eval_metric
              , seed = seed)

    if model.updater != "auto"
        return merge(kwargs, (updater=model.updater,))
    else
        return kwargs
    end
end

## REGRESSOR

mutable struct XGBoostRegressor <:MMI.Deterministic
    num_round::Int
    booster::String
    disable_default_eval_metric::Int
    eta::Float64
    gamma::Float64
    max_depth::Int
    min_child_weight::Float64
    max_delta_step::Float64
    subsample::Float64
    colsample_bytree::Float64
    colsample_bylevel::Float64
    lambda::Float64
    alpha::Float64
    tree_method::String
    sketch_eps::Float64
    scale_pos_weight::Float64
    updater::String
    refresh_leaf::Union{Int,Bool}
    process_type::String
    grow_policy::String
    max_leaves::Int
    max_bin::Int
    predictor::String
    sample_type::String
    normalize_type::String
    rate_drop::Float64
    one_drop
    skip_drop::Float64
    feature_selector::String
    top_k::Int
    tweedie_variance_power::Float64
    objective
    base_score::Float64
    eval_metric
    seed::Int
end

"""
    XGBoostRegressor(; objective="linear", seed=0, kwargs...)

The XGBoost model for univariate targets with `Continuous` element
scitype. Gives deterministic (point) predictions. For possible values
for `objective` and `kwargs`, see
[https://xgboost.readthedocs.io/en/latest/parameter.html](https://xgboost.readthedocs.io/en/latest/parameter.html).

For a time-dependent random seed, use `seed=-1`.

See also: XGBoostCount, XGBoostClassifier

"""
function XGBoostRegressor(
    ;num_round=100
    ,booster="gbtree"
    ,disable_default_eval_metric=0
    ,eta=0.3
    ,gamma=0
    ,max_depth=6
    ,min_child_weight=1
    ,max_delta_step=0
    ,subsample=1
    ,colsample_bytree=1
    ,colsample_bylevel=1
    ,lambda=1
    ,alpha=0
    ,tree_method="auto"
    ,sketch_eps=0.03
    ,scale_pos_weight=1
    ,updater="auto"
    ,refresh_leaf=1
    ,process_type="default"
    ,grow_policy="depthwise"
    ,max_leaves=0
    ,max_bin=256
    ,predictor="cpu_predictor" #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type="uniform"
    ,normalize_type="tree"
    ,rate_drop=0.0
    ,one_drop=0
    ,skip_drop=0.0
    ,feature_selector="cyclic"
    ,top_k=0
    ,tweedie_variance_power=1.5
    ,objective="reg:squarederror"
    ,base_score=0.5
    ,eval_metric="rmse"
    ,seed=0)

    model = XGBoostRegressor(
    num_round
    ,booster
    ,disable_default_eval_metric
    ,eta
    ,gamma
    ,max_depth
    ,min_child_weight
    ,max_delta_step
    ,subsample
    ,colsample_bytree
    ,colsample_bylevel
    ,lambda
    ,alpha
    ,tree_method
    ,sketch_eps
    ,scale_pos_weight
    ,updater
    ,refresh_leaf
    ,process_type
    ,grow_policy
    ,max_leaves
    ,max_bin
    ,predictor #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type
    ,normalize_type
    ,rate_drop
    ,one_drop
    ,skip_drop
    ,feature_selector
    ,top_k
    ,tweedie_variance_power
    ,objective
    ,base_score
    ,eval_metric
    ,seed)

     message = MMI.clean!(model)
     isempty(message) || @warn message

    return model
end

function MMI.clean!(model::XGBoostRegressor)
    warning = ""
    if model.objective == "count:poisson"
        warning *= "Your `objective` suggests prediction of a "*
        "`Count` variable.\n You may want to consider XGBoostCount instead. "
    elseif model.objective in ["reg:logistic", "binary:logistic",
                               "binary:logitraw", "binary:hinge",
                               "multi:softmax", "multi:softprob"]
        warning *="Your `objective` suggests prediction of a "*
        "`Finite` variable.\n You may want to consider XGBoostClassifier "*
        "instead. "
    end
    return warning
end

function MMI.fit(model::XGBoostRegressor
             , verbosity::Int
             , X
             , y)

             silent =
                 verbosity > 0 ?  false : true
    Xmatrix = MMI.matrix(X)
    dm = XGBoost.DMatrix(Xmatrix,label=y)

    objective =
        model.objective in ["linear", "gamma", "tweedie"] ?
            "reg:"*model.objective : model.objective

    seed =
        model.seed == -1 ? generate_seed() : model.seed

    fitresult = XGBoost.xgboost(dm, model.num_round;
                                kwargs(model, silent, seed, objective,model.eval_metric)...)

    features = schema(X).names
    importances = [named_importance(fi, features) for
                   fi in XGBoost.importance(fitresult)]
    cache = nothing

    report = (feature_importances=importances, )

    return fitresult, cache, report

end

function MMI.predict(model::XGBoostRegressor
        , fitresult
        , Xnew)
    Xmatrix = MMI.matrix(Xnew)
    return XGBoost.predict(fitresult, Xmatrix)
end

## COUNT REGRESSOR

mutable struct XGBoostCount <:MMI.Deterministic
    num_round::Int
    booster::String
    disable_default_eval_metric::Int
    eta::Float64
    gamma::Float64
    max_depth::Int
    min_child_weight::Float64
    max_delta_step::Float64
    subsample::Float64
    colsample_bytree::Float64
    colsample_bylevel::Float64
    lambda::Float64
    alpha::Float64
    tree_method::String
    sketch_eps::Float64
    scale_pos_weight::Float64
    updater::String
    refresh_leaf::Union{Int,Bool}
    process_type::String
    grow_policy::String
    max_leaves::Int
    max_bin::Int
    predictor::String
    sample_type::String
    normalize_type::String
    rate_drop::Float64
    one_drop
    skip_drop::Float64
    feature_selector::String
    top_k::Int
    tweedie_variance_power::Float64
    objective
    base_score::Float64
    eval_metric
    seed::Int
end

"""
    XGBoostCount(; seed=0, kwargs...)

The XGBoost model for targets with `Count` scitype. Gives
deterministic (point) predictions. For admissible `kwargs`, see
[https://xgboost.readthedocs.io/en/latest/parameter.html](https://xgboost.readthedocs.io/en/latest/parameter.html).

For a time-dependent random seed, use `seed=-1`.

See also: XGBoostRegressor, XGBoostClassifier

"""
function XGBoostCount(
    ;num_round=100
    ,booster="gbtree"
    ,disable_default_eval_metric=0
    ,eta=0.3
    ,gamma=0
    ,max_depth=6
    ,min_child_weight=1
    ,max_delta_step=0
    ,subsample=1
    ,colsample_bytree=1
    ,colsample_bylevel=1
    ,lambda=1
    ,alpha=0
    ,tree_method="auto"
    ,sketch_eps=0.03
    ,scale_pos_weight=1
    ,updater="auto"
    ,refresh_leaf=1
    ,process_type="default"
    ,grow_policy="depthwise"
    ,max_leaves=0
    ,max_bin=256
    ,predictor="cpu_predictor" #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type="uniform"
    ,normalize_type="tree"
    ,rate_drop=0.0
    ,one_drop=0
    ,skip_drop=0.0
    ,feature_selector="cyclic"
    ,top_k=0
    ,tweedie_variance_power=1.5
    ,objective="count:poisson"
    ,base_score=0.5
    ,eval_metric="rmse"
    ,seed=0)

    model = XGBoostCount(
    num_round
    ,booster
    ,disable_default_eval_metric
    ,eta
    ,gamma
    ,max_depth
    ,min_child_weight
    ,max_delta_step
    ,subsample
    ,colsample_bytree
    ,colsample_bylevel
    ,lambda
    ,alpha
    ,tree_method
    ,sketch_eps
    ,scale_pos_weight
    ,updater
    ,refresh_leaf
    ,process_type
    ,grow_policy
    ,max_leaves
    ,max_bin
    ,predictor #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type
    ,normalize_type
    ,rate_drop
    ,one_drop
    ,skip_drop
    ,feature_selector
    ,top_k
    ,tweedie_variance_power
    ,objective
    ,base_score
    ,eval_metric
    ,seed)

     message = MMI.clean!(model)
     isempty(message) || @warn message

    return model
end

function MMI.clean!(model::XGBoostCount)
    warning = ""
    if(!(model.objective in ["count:poisson"]))
        warning *= "Changing objective to \"poisson\", "*
                       "the only supported value. "
        model.objective="poisson"
    end
    return warning
end

function MMI.fit(model::XGBoostCount
             , verbosity::Int     #> must be here even if unsupported in pkg
             , X
             , y)

    silent = verbosity > 0 ?  false : true

    Xmatrix = MMI.matrix(X)
    dm = XGBoost.DMatrix(Xmatrix,label=y)

    seed =
        model.seed == -1 ? generate_seed() : model.seed

    fitresult = XGBoost.xgboost(dm, model.num_round;
                                kwargs(model, silent, seed, "count:poisson",model.eval_metric)...)
    features = schema(X).names
    importances = [named_importance(fi, features) for
                   fi in XGBoost.importance(fitresult)]

    cache = nothing
    report = (feature_importances=importances, )

    return fitresult, cache, report

end

function MMI.predict(model::XGBoostCount
        , fitresult
        , Xnew)
    Xmatrix = MMI.matrix(Xnew)
    return XGBoost.predict(fitresult, Xmatrix)
end

## CLASSIFIER

mutable struct XGBoostClassifier <:MMI.Probabilistic
    num_round::Int
    booster::String
    disable_default_eval_metric::Int
    eta::Float64
    gamma::Float64
    max_depth::Int
    min_child_weight::Float64
    max_delta_step::Float64
    subsample::Float64
    colsample_bytree::Float64
    colsample_bylevel::Float64
    lambda::Float64
    alpha::Float64
    tree_method::String
    sketch_eps::Float64
    scale_pos_weight::Float64
    updater::String
    refresh_leaf::Union{Int,Bool}
    process_type::String
    grow_policy::String
    max_leaves::Int
    max_bin::Int
    predictor::String
    sample_type::String
    normalize_type::String
    rate_drop::Float64
    one_drop
    skip_drop::Float64
    feature_selector::String
    top_k::Int
    tweedie_variance_power::Float64
    objective
    base_score::Float64
    eval_metric
    seed::Int
end

"""
    XGBoostClassifier(; seed=0, kwargs...)

The XGBoost model for targets with `Finite` scitype (which includes
`Binary=Finite{2}`). Gives probabilistic predictions. For admissible
`kwargs`, see
[https://xgboost.readthedocs.io/en/latest/parameter.html](https://xgboost.readthedocs.io/en/latest/parameter.html).

For a time-dependent random seed, use `seed=-1`.

See also: XGBoostCount, XGBoostRegressor

"""
function XGBoostClassifier(
    ;num_round=100
    ,booster="gbtree"
    ,disable_default_eval_metric=0
    ,eta=0.3
    ,gamma=0
    ,max_depth=6
    ,min_child_weight=1
    ,max_delta_step=0
    ,subsample=1
    ,colsample_bytree=1
    ,colsample_bylevel=1
    ,lambda=1
    ,alpha=0
    ,tree_method="auto"
    ,sketch_eps=0.03
    ,scale_pos_weight=1
    ,updater="auto"
    ,refresh_leaf=1
    ,process_type="default"
    ,grow_policy="depthwise"
    ,max_leaves=0
    ,max_bin=256
    ,predictor="cpu_predictor" #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type="uniform"
    ,normalize_type="tree"
    ,rate_drop=0.0
    ,one_drop=0
    ,skip_drop=0.0
    ,feature_selector="cyclic"
    ,top_k=0
    ,tweedie_variance_power=1.5
    ,objective="automatic"
    ,base_score=0.5
    ,eval_metric="mlogloss"
    ,seed=0)

    model = XGBoostClassifier(
    num_round
    ,booster
    ,disable_default_eval_metric
    ,eta
    ,gamma
    ,max_depth
    ,min_child_weight
    ,max_delta_step
    ,subsample
    ,colsample_bytree
    ,colsample_bylevel
    ,lambda
    ,alpha
    ,tree_method
    ,sketch_eps
    ,scale_pos_weight
    ,updater
    ,refresh_leaf
    ,process_type
    ,grow_policy
    ,max_leaves
    ,max_bin
    ,predictor #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type
    ,normalize_type
    ,rate_drop
    ,one_drop
    ,skip_drop
    ,feature_selector
    ,top_k
    ,tweedie_variance_power
    ,objective
    ,base_score
    ,eval_metric
    ,seed)

     message = MMI.clean!(model)           #> future proof by including these
     isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::XGBoostClassifier)
    warning = ""
    if(!(model.objective =="automatic"))
            warning *="Changing objective to \"automatic\", the only supported value. "
            model.objective="automatic"
    end
    return warning
end

function MMI.fit(model::XGBoostClassifier
                     , verbosity::Int     #> must be here even if unsupported in pkg
                     , X
                     , y)
    Xmatrix = MMI.matrix(X)

    a_target_element = y[1] # a CategoricalValue or CategoricalString
    num_class = length(MMI.classes(a_target_element))

    # The eval metric is different depending on the class size.
    eval_metric = model.eval_metric
    if num_class == 2 && eval_metric == "mlogloss"
        eval_metric = "logloss"
    end
    if num_class > 2 && eval_metric == "logloss"
        eval_metric = "mlogloss"
    end

    y_plain = MMI.int(y) .- 1 # integer relabeling should start at 0

    if(num_class==2)
        objective="binary:logistic"
        y_plain = convert(Array{Bool}, y_plain)
    else
        objective="multi:softprob"
    end

    silent =
        verbosity > 0 ?  false : true

    seed =
        model.seed == -1 ? generate_seed() : model.seed

    if num_class == 2
        # XGBoost API requires that num_class is not passed in when n=2.
        result = XGBoost.xgboost(Xmatrix, label=y_plain, model.num_round;
                             kwargs(model, silent, seed, objective,eval_metric)...)
    else
        result = XGBoost.xgboost(Xmatrix, label=y_plain, model.num_round;
                                 num_class=num_class,
                                 kwargs(model, silent, seed, objective,eval_metric)...)
    end
    fitresult = (result, a_target_element)

    features = schema(X).names

    importances = [named_importance(fi, features) for
                   fi in XGBoost.importance(result)]

    cache = nothing
    report = (feature_importances=importances, )

    return fitresult, cache, report

end

function MMI.predict(model::XGBoostClassifier
        , fitresult
        , Xnew)

    result, a_target_element = fitresult
    decode = MMI.decoder(a_target_element)
    classes = MMI.classes(a_target_element)

    Xmatrix = MMI.matrix(Xnew)
    XGBpredictions = XGBoost.predict(result, Xmatrix)

    nlevels = length(classes)
    npatterns = MMI.nrows(Xnew)

    if nlevels == 2
        true_class_probabilities = reshape(XGBpredictions, 1, npatterns)
        false_class_probabilities = 1 .- true_class_probabilities
        XGBpredictions = vcat(false_class_probabilities, true_class_probabilities)
    end

    prediction_probabilities = reshape(XGBpredictions, nlevels, npatterns)

    # note we use adjoint of above:
    predictions = MMI.UnivariateFinite(classes, prediction_probabilities')

    return predictions
end

## SERIALIZATION - REGRESSOR AND COUNT

const XGBoostInfinite = Union{XGBoostRegressor,XGBoostCount}

function MLJModelInterface.save(filename,
                                ::XGBoostInfinite,
                                fitresult; # `XGBoost.Booster` object
                                kwargs...)
    xgb_filename = string(filename, ".xgboost.model")
    XGBoost.save(fitresult, xgb_filename)
    serializable_fitresult = read(xgb_filename)
    @info "Additional XGBoost serialization file \"$xgb_filename\" generated. "
    return serializable_fitresult
end

function MLJModelInterface.restore(filename,
                                   ::XGBoostInfinite,
                                   serializable_fitresult)
    xgb_filename = string(filename, ".tmp")
    open(xgb_filename, "w") do file
        write(file, serializable_fitresult)
    end
    fitresult = XGBoost.Booster(model_file=xgb_filename)
    rm(xgb_filename)
    return fitresult
end

## SERIALIZATION - CLASSIFIER

function MLJModelInterface.save(filename,
                                ::XGBoostClassifier,
                                fitresult;
                                kwargs...)
    booster, a_target_element = fitresult

    xgb_filename = string(filename, ".xgboost.model")
    XGBoost.save(booster, xgb_filename)
    persistent_booster = read(xgb_filename)
    @info "Additional XGBoost serialization file \"$xgb_filename\" generated. "
    return (persistent_booster, a_target_element)
end

function MLJModelInterface.restore(filename,
                                   ::XGBoostClassifier,
                                   serializable_fitresult)
    persistent_booster, a_target_element = serializable_fitresult

    xgb_filename = string(filename, ".tmp")
    open(xgb_filename, "w") do file
        write(file, persistent_booster)
    end
    booster = XGBoost.Booster(model_file=xgb_filename)
    rm(xgb_filename)
    fitresult = (booster, a_target_element)
    return fitresult
end

## METADATA

XGTypes=Union{XGBoostRegressor,XGBoostCount,XGBoostClassifier}

MMI.package_name(::Type{<:XGTypes}) = "XGBoost"
MMI.package_uuid(::Type{<:XGTypes}) = "009559a3-9522-5dbb-924b-0b6ed2b22bb9"
MMI.package_url(::Type{<:XGTypes}) = "https://github.com/dmlc/XGBoost.jl"
MMI.is_pure_julia(::Type{<:XGTypes}) = false

MMI.load_path(::Type{<:XGBoostRegressor}) = "$PKG.XGBoostRegressor"
MMI.input_scitype(::Type{<:XGBoostRegressor}) = Table(Continuous)
MMI.target_scitype(::Type{<:XGBoostRegressor}) = AbstractVector{Continuous}
MMI.docstring(::Type{<:XGBoostRegressor}) =
    "The XGBoost gradient boosting method, for use with "*
    "`Continuous` univariate targets. "

MMI.load_path(::Type{<:XGBoostCount}) = "$PKG.XGBoostCount"
MMI.input_scitype(::Type{<:XGBoostCount}) = Table(Continuous)
MMI.target_scitype(::Type{<:XGBoostCount}) = AbstractVector{Count}
MMI.docstring(::Type{<:XGBoostCount}) =
    "The XGBoost gradient boosting method, for use with "*
    "`Count` univariate targets, using a Poisson objective function. "

MMI.load_path(::Type{<:XGBoostClassifier}) = "$PKG.XGBoostClassifier"
MMI.input_scitype(::Type{<:XGBoostClassifier}) = Table(Continuous)
MMI.target_scitype(::Type{<:XGBoostClassifier}) = AbstractVector{<:Finite}
MMI.docstring(::Type{<:XGBoostClassifier}) =
    "The XGBoost gradient boosting method, for use with "*
    "`Finite` univariate targets (`Multiclass`, "*
    "`OrderedFactor` and `Binary=Finite{2}`)."

end # module


text:
using Test
using HyperbolicDiffEq

@elapsed begin

  @time @testset "Riemann Problems" begin include("riemann_problems_test.jl") end
  @time @testset "Finite Volumes" begin include("finite_volumes_test.jl") end
  @time @testset "Monomial Fluxes" begin include("monomial_fluxes_test.jl") end
  @time @testset "Flux Difference" begin include("flux_difference_test.jl") end
  @time @testset "ENO Reconstruction" begin include("modified_eno_test.jl") end
  @time @testset "TeCNO" begin include("tecno_test.jl") end
  @time @testset "Quartic Nonconvex Flux" begin include("balance_laws/quartic_nonconvex_test.jl") end
  @time @testset "Keyfitz-Kranzer System" begin include("balance_laws/keyfitz_kranzer_test.jl") end

end


text:
function molzahn_heuristic(clique::AbstractVector, edge; merge_kwargs=merge_kwargs)
    size_intersect = length(intersect(clique[edge[1]], clique[edge[2]]))
    size_clique1 = length(clique[edge[1]])
    size_clique2 = length(clique[edge[2]])
    size_fusion = size_clique1 + size_clique2 - size_intersect
    delta = size_fusion * (2 * size_fusion + 1) -
            size_clique1 * (2 * size_clique1 + 1) -
            size_clique2 * (2 * size_clique2 + 1) -
            size_intersect * (2 * size_intersect + 1)
    return delta
end

function update_molzahn(;kwargs...)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


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

using Cxx

cxxinclude("TStyle.h")

export set_default_palette

const color_palette_ids = Dict(
    :deep_sea => 51,
    :grey_scale => 52,
    :dark_body_radiator => 53,
    :blue_yellow => 54,
    :rain_bow => 55,
    :inverted_dark_body_radiator => 56,
    :bird => 57,
    :cubehelix => 58,
    :green_red_violet => 59,
    :blue_red_yellow => 60,
    :ocean => 61,
    :color_printable_on_grey => 62,
    :alpine => 63,
    :aquamarine => 64,
    :army => 65,
    :atlantic => 66,
    :aurora => 67,
    :avocado => 68,
    :beach => 69,
    :black_body => 70,
    :blue_green_yellow => 71,
    :brown_cyan => 72,
    :cmyk => 73,
    :candy => 74,
    :cherry => 75,
    :coffee => 76,
    :dark_rain_bow => 77,
    :dark_terrain => 78,
    :fall => 79,
    :fruit_punch => 80,
    :fuchsia => 81,
    :grey_yellow => 82,
    :green_brown_terrain => 83,
    :green_pink => 84,
    :island => 85,
    :lake => 86,
    :light_temperature => 87,
    :light_terrain => 88,
    :mint => 89,
    :neon => 90,
    :pastel => 91,
    :pearl => 92,
    :pigeon => 93,
    :plum => 94,
    :red_blue => 95,
    :rose => 96,
    :rust => 97,
    :sandy_terrain => 98,
    :sienna => 99,
    :solar => 100,
    :south_west => 101,
    :starry_night => 102,
    :sunset => 103,
    :temperature_map => 104,
    :thermometer => 105,
    :valentine => 106,
    :visible_spectrum => 107,
    :water_melon => 108,
    :cool => 109,
    :copper => 110,
    :gist_earth => 111,
    :viridis => 112
)

set_default_palette(palette::Symbol, alpha::Real = 1.0) = begin
    id = color_palette_ids[palette]
    icxx"gStyle->SetPalette($id, 0, $alpha);"
    nothing
end


text:
function runGillMatsuno(
    S  :: Simulation{FT},
    G  :: Grid{FT},
    Qv :: Vector{QStructs{FT}},
    D  :: Domain{FT},
) where FT <: Real

	@info "$(now()) - Initializing the Gill-Matsuno Model ..."

    nx,ny = G.nx,G.ny

    Q  = createQ(Qv,G)
    rt = floor(Int32,S.tt/S.δt) + 1	# Number of running steps
    nt = floor(Int32,S.tt/S.ft) + 1	# Number of output steps
    st = 0.                 	# Simulation time
    ot = 0.                  	# Indicator as to whether to save output
    mt = 1                  	# Index in output matrix

	@info "$(now()) - The simulation will run for $(rt-1) timesteps, and output the fields for $(nt-1) timesteps"

	@info "$(now()) - Preallocating arrays for the wind and potential fields ..."

    ϕ  = zeros(FT,nx,ny);    u  = zeros(FT,nx,ny);    v  = zeros(FT,nx,ny+1);
    ϕn = zeros(FT,nx,ny);    un = zeros(FT,nx,ny);    vn = zeros(FT,nx,ny+1);
    ϕf = zeros(FT,nx,ny,nt); uf = zeros(FT,nx,ny,nt); vf = zeros(FT,nx,ny+1,nt);
	t  = zeros(nt)

    for it = 1 : rt

        ϕfield!(ϕn,ϕ,u,v,G,D,S,Q)
        ufield!(un,u,ϕ,v,G,D,S)
        vfield!(vn,v,ϕ,u,G,D,S)
        ϕ .= ϕn
		u .= un
		v .= vn

		if isnan(sum(ϕ)) || isnan(sum(u)) || isnan(sum(v))
			error("There are NaN values in the fields, indicating that the model is unstable.  Please reduce the timestep")
		end

        st += S.δt
        if st >= ot
            ot += S.ft
            if !isone(it)
            	mt += 1
				@info "$(now()) - Extracting output at $(@sprintf("%06.2f",st)) model seconds ..."
                ϕf[:,:,mt] .= ϕ
                uf[:,:,mt] .= u
                vf[:,:,mt] .= v
				t[mt] = st
            end
        end

    end

	savefields(uf,vf,ϕf,t,G,S)

end

function savefields(
	u :: Array{FT},
	v :: Array{FT},
	ϕ :: Array{FT},
	t :: Array{<:Real},
	G :: Grid{FT},
	S :: Simulation{FT}
) where FT <: Real

	if isfile(S.fnc)
        @info "$(now()) - Stale NetCDF file $(S.fnc) detected.  Overwriting ..."
        rm(S.fnc);
    end
	ds = NCDataset(S.fnc,"c",attrib = Dict(
        "Conventions" => "CF-1.6",
        "history"     => "Created as output by the GillMatsuno.jl package on $(Dates.now())"
    ));

	ds.dim["xC"] = G.nx
    ds.dim["yC"] = G.ny
    ds.dim["xF"] = G.nx
    ds.dim["yF"] = G.ny + 1
    ds.dim["time"] = length(t)

	ncxC = defVar(ds,"xC",FT,("xC",),attrib = Dict("long_name" => "X center",))
	ncyC = defVar(ds,"yC",FT,("yC",),attrib = Dict("long_name" => "Y center",))
	ncxF = defVar(ds,"xF",FT,("xF",),attrib = Dict("long_name" => "X face",))
	ncyF = defVar(ds,"yF",FT,("yF",),attrib = Dict("long_name" => "Y face",))
	nct  = defVar(ds,"time",Float64,("time",),attrib = Dict("long_name" => "time",))

	ncu = defVar(ds,"u",FT,("xF","yC","time"),attrib = Dict(
		"long_name" => "u_component_of_wind",
		"full_name" => "Nondimensional u-Wind",
		"units"		=> "L/T",
	))

	ncv = defVar(ds,"v",FT,("xC","yF","time"),attrib = Dict(
		"long_name" => "v_component_of_wind",
		"full_name" => "Nondimensional v-Wind",
		"units"		=> "L/T",
	))

	ncϕ = defVar(ds,"ϕ",FT,("xC","yC","time"),attrib = Dict(
		"full_name" => "geopotential",
		"full_name" => "Nondimensional Height",
		"units"		=> "H",
	))

	ncxC[:] = G.xc
	ncyC[:] = G.yc
	ncxF[:] = G.xf
	ncyF[:] = G.yf

	nct[:] = t
	ncu[:] = u
	ncv[:] = v
	ncϕ[:] = ϕ

	close(ds)

	@info "$(now()) - Output fields from the Gill-Matsuno run has been saved."

end

function ϕfield!(
    ϕn::Array{FT}, ϕ::Array{FT}, u::Array{FT}, v::Array{FT},
    G::Grid{FT}, D::Domain{FT}, S::Simulation{FT}, Q::Array{FT}
) where FT <: Real

    for jj = 1 : G.ny, ii = 1 : G.nx

		vy = (v[ii,jj+1] - v[ii,jj]) / G.δy

        if ii != 1
			ux = (u[ii,jj] - u[ii-1,jj]) / G.δx
		else
			ux = (u[1,jj] - u[G.nx,jj]) / G.δx
        end

        ϕn[ii,jj] = ϕ[ii,jj] - S.δt * (D.α * ϕ[ii,jj] + D.g * D.H * (ux+vy) + Q[ii,jj])

    end

    return

end

function ufield!(
    un::Array{FT}, u::Array{FT}, ϕ::Array{FT}, v::Array{FT},
    G::Grid{FT}, D::Domain{FT}, S::Simulation{FT}
) where FT <: Real

    for jj = 1 : G.ny, ii = 1 : G.nx

        if ii != G.nx
			ϕx = (ϕ[ii+1,jj] - ϕ[ii,jj]) / G.δx
			vi = 0.25 * (v[ii,jj]+v[ii+1,jj]+v[ii,jj+1]+v[ii+1,jj+1])
		else
			ϕx = (ϕ[1,jj] - ϕ[G.nx,jj]) / G.δx
			vi = 0.25 * (v[G.nx,jj]+v[1,jj]+v[G.nx,jj+1]+v[1,jj+1])
        end

        un[ii,jj] = u[ii,jj] - S.δt * (D.α * u[ii,jj] - D.β * G.yc[jj] * vi + ϕx)

    end

    return

end

function vfield!(
    vn::Array{FT}, v::Array{FT}, ϕ::Array{FT}, u::Array{FT},
    G::Grid{FT}, D::Domain{FT}, S::Simulation{FT}
) where FT <: Real

    for jj = 2 : G.ny, ii = 1 : G.nx

		ϕy = (ϕ[ii,jj] - ϕ[ii,jj-1]) / G.δy

        if ii != 1
			ui = 0.25 * (u[ii,jj]+u[ii-1,jj]+u[ii,jj-1]+u[ii-1,jj-1])
		else
			ui = 0.25 * (u[1,jj]+u[G.nx,jj]+u[1,jj-1]+u[G.nx,jj-1])
        end

        vn[ii,jj] = v[ii,jj] - S.δt * (D.α * v[ii,jj] + D.β * G.yf[jj] * ui + ϕy)

    end

    return

end


text:
"""
Can we avoid stepping trough R to leverage functionality form ergm R package? Most of it is written in C++.

Not doable because Cxx.jl does not work with juia 1.5

"""

using Cxx
import ScoreDrivenERGM:ErgmRcall

using ScoreDrivenERGM.Scalings
using ScoreDrivenERGM.Utilities
using ScoreDrivenERGM:StaticNets

using PyPlot
using RCall
using Statistics

ErgmRcall.clean_start_RCall()
ergmTermsString = "edges +  mutual"
R"""options(warn=-1) """

model = StaticNets.ErgmDirBin0Rec0()

##-------------------- Test and COmpare MLE and MPLE estimates
#Compare for a single value of the parameters
N=100
Ldgp = n_pox_dir_links(N)/4
Rdgp = n_pox_pairs(N)/4-100

θ_0, η_0 = Tuple(StaticNets.estimate(model, Ldgp, Rdgp, N))

nSample = 100
diadProb = StaticNets.diadProbFromPars(model, [θ_0, η_0])
A_vec = [StaticNets.samplSingMatCan(model, diadProb, N) for i=1:nSample]
statsVec = reduce(hcat,[StaticNets.stats_from_mat(model, A) for A in A_vec])

jnum = 10

cxx"""
        void printme(int x) {
            std::cout << x << std::endl;
        }
    """
@cxx printme(jnum)

function c_list(model, )


text:
#
#   Auxiliarfunctions.jl file:
#       ErrorEst:
#       MyNorm
#       Rdigits

function ErrorEst(U,F,dt,beta,abstol,reltol)

    uiType=eltype(U[1])

    (s,)=size(F)
	D=length(U[1])

	est=zero(uiType)

    @inbounds begin
	for k in eachindex(U[1])

		sum=zero(uiType)
		maxU=zero(uiType)
		for is in 1:s
        	sum+=beta[is]*F[is][k]
			maxU=max(maxU,abs(U[is][k]))
        end

		est+=(abs(dt*sum))^2/(abstol+maxU^2*reltol)

	end
    end

    return(est/D)

end

function MyNorm(u,abstol,reltol)

    uiType=eltype(u)
	norm=zero(uiType)

    @inbounds begin
	for k in eachindex(u)
	    aux=u[k]/(abstol+abs(u[k])*reltol)
		norm+=aux*aux
	end
    end

	norm=sqrt(norm/(length(u)))

    return(norm)

end

function Rdigits(x::Real,r::Integer)

	aux=copy(x)
	mx=r*aux
	mxx=mx+aux
	return (mxx-mx)

end


text:
using JLD
using PyPlot
function plot_sensitivity_cumsum()
X = load("../data/sensitivity_t2000_dt01_n500.jld")

dJds = X["dJdbeta"]
vsh = X["vsh"]
d, n, n_s = size(vsh)
t = 0.01*n
thres = 20
n_count = 0
mean_dJds = zeros(2)
var_dJds = zeros(2)
valid_samples = zeros(Int64, n_s)
for i = 1:n_s
	flag = 0
	for j = 1:n
		if norm(vsh[:,j,i]) > thres
			flag = 1
		end
	end
	if flag == 0
		valid_samples[i] = 1
		n_count += 1
		mean_dJds .+= dJds[:,i]
		var_dJds .+= dJds[:,i].*dJds[:,i]
	end
end

mean_dJds ./= n_count
var_dJds ./= n_count
var_dJds .-= mean_dJds.*mean_dJds

fig, ax = subplots(1,2)
ax[1].xaxis.set_tick_params(labelsize=20)
ax[1].yaxis.set_tick_params(labelsize=20)
ax[2].xaxis.set_tick_params(labelsize=20)
ax[2].yaxis.set_tick_params(labelsize=20)

dJds = dJds[:,(!iszero).(valid_samples)]
n_s = n_count

ax[1].plot(t*(1:n_s), cumsum(dJds[1,:])./(1:n_s), "bo", label=L"$\dfrac{d \langle J_{\rm ac}\rangle}{d \beta}$")
ax[2].plot(t*(1:n_s), cumsum(dJds[2,:])./(1:n_s), "ro", label=L"$\dfrac{d \langle J_{\rm ray}\rangle}{d \beta}$")

ax[1].set_xlabel("Time",fontsize=20)
ax[1].set_ylabel("Shadowing sensitivities",fontsize=20)

ax[2].set_xlabel("Time",fontsize=20)
#ax[2].set_ylabel("Shadowing sensitivities",fontsize=20)
fig.legend(fontsize=20)
ax[1].grid(true)
ax[2].grid(true)

fig1, ax1 = subplots(1,1)
ax1.xaxis.set_tick_params(labelsize=20)
ax1.yaxis.set_tick_params(labelsize=20)

ax1.errorbar(1:n_count, dJds[1,:], sqrt(var_dJds[1]), marker="o", mfc="blue", mec="blue", label=L"$\dfrac{d \langle J_{\rm ac}\rangle}{d \beta}$")
ax1.set_xlabel("Sample #",fontsize=20)
ax1.set_ylabel("Shadowing sensitivities",fontsize=20)

fig1.legend(fontsize=20)
ax1.grid(true)

fig2, ax2 = subplots(1,1)
ax2.xaxis.set_tick_params(labelsize=20)
ax2.yaxis.set_tick_params(labelsize=20)

ax2.errorbar((1:n_count), dJds[2,:], sqrt(var_dJds[2]), marker="o", mfc="red", mec="red", color="red", label=L"$\dfrac{d \langle J_{\rm ray}\rangle}{d \beta}$")

ax2.set_xlabel("Sample #",fontsize=20)
ax2.set_ylabel("Shadowing sensitivities",fontsize=20)

fig2.legend(fontsize=20)
ax2.grid(true)
return fig1, ax1, fig2, ax2, dJds
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
# Cell Index
# ==========

# Runtime Statistics
# ------------------

mutable struct RuntimeStats
    preproc_time::UInt64
    search_time::UInt64
    rank_time::UInt64
    timer::UInt64
end

RuntimeStats() = RuntimeStats(0, 0, 0, 0)

function reset!(s::RuntimeStats)
    s.preproc_time = 0
    s.search_time = 0
    s.rank_time = 0
    s.timer = 0
    return s
end

@inline function tic!(s::RuntimeStats)
    t = time_ns()
    ns = t - s.timer
    s.timer = t
    return ns
end

# Locality-Sensitive Hash
# -----------------------

struct LSHash{T<:BitVec}
    projection::Matrix{Float32}
    hammingindex::HammingIndexes.HammingIndex{T}
end

bitvectype(::LSHash{T}) where T = T

function findneighbors!(neighbors::Matrix{Int},
                        Z::Vector{T},
                        X::Matrix{Float32},
                        lshash::LSHash{T}) where T
    @assert length(Z) == size(neighbors, 2) == size(X, 2)
    sketch!(Z, X, lshash)
    # search for k-nearest neighbors
    k = size(neighbors, 1)
    nns = HammingIndexes.NearestNeighbors(k)
    for j in 1:length(Z)
        HammingIndexes.findknn!(nns, Z[j], lshash.hammingindex)
        for i in 1:k
            neighbors[i,j] = nns.indexes[i]
        end
    end
    return neighbors, Z
end

# Draw sketches of X.
function sketch!(Z::Vector{T}, X::Matrix{Float32}, lshash::LSHash{T}) where T
    return sketch!(Z, X, lshash.projection)
end

function sketch!(Z::Vector{T}, X::Matrix{Float32}, P::Matrix{Float32}) where {T<:BitVec}
    n = size(X, 2)
    @assert n == length(Z)
    PX = P * X
    @assert size(PX, 1) == bitsof(T)
    for j in 1:n
        z::T = 0
        i = (j - 1) * bitsof(T) + 1
        @inbounds for k in 0:bitsof(T)-1
            z |= T(PX[i+k] ≥ 0) << k
        end
        Z[j] = z
    end
    return Z
end

function generate_random_projections(K::Integer, D::Integer, superbit::Integer)
    @assert 1 ≤ superbit ≤ min(K, D)
    q, r = divrem(K, superbit)
    Ps = Matrix{Float32}[]
    for _ in 1:q
        # generate an orthogonal random matrix.
        M = randn(Float32, D, superbit)
        push!(Ps, Matrix(qr!(M).Q)')
    end
    if r > 0
        M = randn(Float32, D, r)
        push!(Ps, Matrix(qr!(M).Q)')
    end
    return vcat(Ps...)
end

# CellIndex
# ---------

struct CellIndex
    # preprocessor
    preproc::Preprocessor
    # locality-sensitive hashes
    lshashes::Vector{LSHash{T}} where T
    # any user-provided metadata (e.g. cellnames)
    metadata::Any
    # compressed count matrix (optional)
    counts::Union{CMatrix{Int32},Nothing}
    # run-time statistics
    rtstats::RuntimeStats
end

nbits(index::CellIndex) = bitsof(bitvectype(index.lshashes[1]))
ncells(index::CellIndex) = length(index.lshashes[1].hammingindex)

function Base.show(io::IO, index::CellIndex)
    print(io, summary(index), "(<#cells=$(ncells(index)), hash=$(nbits(index))×$(length(index.lshashes))>)")
end

"""
    CellIndex(counts, features; <keyword arguments>...)

Create a cell index from a count matrix `counts` along with `features`.

Dimensionality reduction and data hashing generate puseudo-random numbers, and
therefore it is recommended to fix the random number seed by calling
`Random.seed!` for reproducibility.

Arguments
---------

Required arguments (positional):
- `counts`: Transcriptome expression matrix (features x cells).
- `features`:
    Features used to compare expression profiles.
    It must be an object of `Features`.

Additional data:
- `metadata=nothing`: An arbitrary metadata (e.g., cell types).
- `index=true`: Create bit index(es) or not.
- `keep_counts=false`: Keep filtered counts in the index object.

Parameters for preprocessing:
- `normalize=true`: Normalize library sizes or not.
- `scalefactor=1.0e4`: The scale factor of cell-wise counts (or library sizes).
- `transformer=:log1p`: Variance-stabilizing transformer (`:log1p` or `:ftt`).
- `n_dims=50`: The number of dimensions of principal components.
               This value must be at least as large as the number of selected features.
- `randomize=true`: Use the randomized version of SVD or not.
- `standardize=true`: Standardize features or not.

Parameters for locality-sensitive hashing:
- `n_bits=128`: The number of bits (64, 128, 256, or 512).
- `n_lshashes=4`: The number of independent locality-sensitive hashes.
- `superbit=min(n_dims, n_bits)`: The depth of super-bits.
"""
function CellIndex(
        # required arguments
        counts::AbstractMatrix{<:Union{S,Missing}},
        features::Features;
        # additional data
        metadata=nothing,
        index::Bool=true,
        keep_counts::Bool=false,
        # parameters for preprocessing
        normalize::Bool=true,
        scalefactor::Real=1.0e4,
        transformer::Symbol=:log1p,
        n_dims::Integer=50,
        randomize::Bool=true,
        standardize::Bool=true,
        # parameters for LSH
        n_bits::Integer=128,
        n_lshashes::Integer=4,
        superbit::Integer=min(n_dims, n_bits),
       ) where S <: Real
    # check arguments
    m, n = size(counts)
    if nfeatures(features) != m
        throw(ArgumentError("mismatching features size"))
    end
    if nselected(features) == 0
        throw(ArgumentError("no selected features"))
    end
    if !(1 ≤ n_dims ≤ nselected(features))
        throw(ArgumentError("invalid n_dims (possibly too large)"))
    end
    if n_bits ∉ (64, 128, 256, 512)
        throw(ArgumentError("invalid n_bits"))
    end
    if !(1 ≤ superbit ≤ min(n_dims, n_bits))
        throw(ArgumentError("invalid superbit"))
    end
    if !(scalefactor > 0)
        throw(ArgumentError("invalid scalefactor"))
    end
    if !(1 ≤ superbit ≤ min(n_dims, n_bits))
        throw(ArgumentError("invalid superbit"))
    end
    if transformer ∉ (:log1p, :ftt)
        throw(ArgumentError("invalid transformer"))
    end
    # filter features
    featurenames = selectedfeatures(features)
    counts = Matrix{S}(counts[features.selected,:])
    Y = ExpressionMatrix(convert(Matrix{Float32}, counts), featurenames)
    # make cell sketches
    T = n_bits ==  64 ? BitVec64  :
        n_bits == 128 ? BitVec128 :
        n_bits == 256 ? BitVec256 :
        n_bits == 512 ? BitVec512 :
        @assert(false)
    if transformer == :log1p
        transformer = LogT(1.0f0)
    else
        @assert(transformer == :ftt)
        transformer = FTT()
    end
    # preprocess expression profiles
    preproc = Preprocessor(Y, transformer, PCA(n_dims, randomize=randomize), normalize, standardize, Float32(scalefactor))
    X = preprocess(preproc, Y, :database)
    # hash preprocessed data
    lshashes = LSHash{T}[]
    for _ in 1:n_lshashes
        P = generate_random_projections(n_bits, n_dims, superbit)
        Z = Vector{T}(undef, n)
        sketch!(Z, X, P)
        push!(lshashes, LSHash(P, HammingIndexes.HammingIndex(Z, index=index)))
    end
    # compress counts if required
    if keep_counts
        ccounts = CMatrix(counts, mmap=true)
    else
        ccounts = nothing
    end
    # create a cell index
    return CellIndex(preproc, lshashes, metadata, ccounts, RuntimeStats())
end

"""
    save(filename::AbstractString, index::CellIndex)

Save `index` to a file.

See also `load`.
"""
function save(filename::AbstractString, index::CellIndex)
    Base.open(filename, "w") do output
        serialize(output, index.preproc)
        serialize(output, index.lshashes)
        serialize(output, index.metadata)
        serialize(output, index.rtstats)
        if index.counts !== nothing
            counts = index.counts
            write(output, Int64(counts.size[1]))
            write(output, Int64(counts.size[2]))
            write(output, Int64(sizeof(counts.data)))
            write(output, counts.data)
        end
    end
    return nothing
end

"""
    load(filename::AbstractString; mmap::Bool=true)

Load an index from `filename`.

If `mmap` is `true`, the expression matrix stored in the file will be loaded
using memory-mapped file to save the memory space.

See also `save`.
"""
function load(filename::AbstractString; mmap::Bool=true)
    Base.open(filename, "r") do input
        preproc = deserialize(input)
        lshashes = deserialize(input)
        metadata = deserialize(input)
        rtstats = deserialize(input)
        if eof(input)
            counts = nothing
        else
            m = read(input, Int64)
            n = read(input, Int64)
            len = read(input, Int64)
            if mmap
                data = Mmap.mmap(filename, Vector{UInt8}, len, position(input), grow=false)
            else
                data = read(input, len)
            end
            counts = CMatrix{Int32}(data, (m, n))
        end
        return CellIndex(preproc, lshashes, metadata, counts, rtstats)
    end
end


text:
export Uniform

import Distributions: Uniform

vague(::Type{ <: Uniform }) = Uniform(-100.0, 100.0)


text:
include("../src/ccsa.jl")
using LinearAlgebra
using CatViews
using BenchmarkTools
using Profile
using ProfileView
using PProf
using SparseArrays
using BenchmarkTools
ncon = 0
#A=sprand(Float64, nvar, nvar, 0.3)
β=0.1
#x_true=sprand(Float64, nvar, 1, 0.3)
time=Array{Float64}(undef,10000)
pp=[]
iter=Array{Int}(undef,10000)
for n in 5000:100:5000
    A=Matrix(SymTridiagonal(2*ones(n),ones(n)))
    #A=SymTridiagonal(2*ones(n),ones(n))
    x_true=ones(n)*10
    y=A*x_true
    function f_and_∇f(x)
        f=[sum((A*x-y).^2)+β*norm(x,1)]
        ∇f=(A*x-y)'*A.+β.*signbit.(x)'
        return f,∇f
    end
    sparse_opt = CCSAState(n, ncon, f_and_∇f,xtol_rel=10e-4)
    time[n]=@belapsed begin
        sparse_opt = CCSAState($n, 0, $f_and_∇f)
        optimize(sparse_opt)
    end samples=5 evals=1
    #pp=optimize(sparse_opt)
    println(n)
    println(time[n])
end
plot(1:length(pp),pp,yscale=:log10,xlim=(-10,550),ylim=(0.5e3,1e7),label="CCSA",xlabel="iterations",ylabel="objective function value")
using Plots
plot(250:1000,time[250:1000],xlim=(150,1100),ylim=(0.00005,0.0002),label="time-dim graph")
plot(100:100:1000,time[100:100:1000],label="time-dim graph")
plot(1:1000,iter[1:1000])

time_cun
time_jiu
@belapsed begin
    1==1
end samples=100 evals=1
# Dense Jacobian should not be allocated if sparsity is provided

#
#=
a1=2;b1=0;a2=-1;b2=1;
function f(x)
    return [sqrt(x[2]),(a1*x[1]+b1)^3-x[2],(a2*x[1]+b2)^3-x[2]]
end
function ∇f(x) #(m+1)*n
    return [0 1/2/x[2]; 3*a1*(a1*x[1]+b1)^2 -1; 3*a2*(a2*x[1]+b2)^2 -1]
end
function f_and_∇f(x)
    return f(x),∇f(x)
end
n=2
m=2
σ=[30.0, 30.0] #n
ρ=[100.0, 100.0, 100.0]
x=[1.234, 5.678]
lb=[-Inf, 0.]
opt=CCSAState(n,m,f_and_∇f,ρ,σ,x)
opt.lb=lb
opt.xtol_rel = 1e-4
@time optimize(opt)
ProfileView.@profview optimize(opt)
=#


text:
function simp13(f,a,b,h)
    t = collect(Float64,a:h:b)
    oddsum, evensum = 0,0 #initialize
    for  j in 2:2:(length(t) -1)
        oddsum += f(t[j])
    end
    for j in 3:2:(length(t) -1)
        evensum += f(t[j])
    end
    I = (b-a)* ( f(t[1]) + 4*oddsum + 2*evensum + f(t[end]) )/(3*(length(t)-1))
    return I
end

function simp38(f,a,b,h)
    t = collect(Float64,a:h:b)
    return (b-a)/8*( f(t[1]) + 3*f(t[2]) + 3*f(t[3]) + f(t[4]) )
end

function testsimp2(f,a,b,h)
    t = collect(Float64,a:h:b)
    if mod(length(t)-1,2) == 0 #i.e. we have an even number of intervals
        return simp13(f,a,b,h)
    else #i.e. we have an odd number of intervals
        I1 = simp13(f,a,t[end-3],h)
        I2 = simp38(f,t[end-3],b,h)
        return I1 + I2
    end
end

#This seems to be the fastest
function intsimp(f,a,b,h)
    t = collect(Float64,a:h:b)

    if mod(length(t)-1,2) == 0 #i.e. we have an even number of intervals
        oddsum, evensum = 0,0 #initialize
        for  j in 2:2:(length(t) -1)
            oddsum += f(t[j])
        end
        for j in 3:2:(length(t) -1)
            evensum += f(t[j])
        end
        I = (b-a)* ( f(t[1]) + 4*oddsum + 2*evensum + f(t[end]) )/(3*(length(t)-1))
        return I
    else #i.e. we have an odd number of intervals
        # First we need to calculate the simpson 1/3 up until... end-3
        # since we have N-3 points for simpson 1/3 and N-4 intervals
        oddsum, evensum = 0,0 #initialize
        for  j in 2:2:(length(t) - 4)
            oddsum += f(t[j])
        end
        for j in 3:2:(length(t) - 4)
            evensum += f(t[j])
        end
        I1 = (t[end-3]-a)* ( f(t[1]) + 4*oddsum + 2*evensum + f(t[end-3]) )/(3*(length(t)-4))

        I2 = (b-t[end-3])/8*( f(t[end-3]) + 3*f(t[end-2]) + 3*f(t[end-1]) + f(t[end]) )

        return I1 + I2

    end
end

# Integration whenever a pair of vectors is given as a function, x and y
function intsimp(x::Array{Float64,1},y::Array{Float64,1})
    h = abs( x[2] - x[1] )

    if mod(length(x)-1,2) == 0 #i.e. we have an even number of intervals
        oddsum, evensum = 0,0 #initialize
        for  j in 2:2:(length(x) -1)
            oddsum += y[j]
        end
        for j in 3:2:(length(x) -1)
            evensum += y[j]
        end
        I = (x[end]-x[1])* ( y[1] + 4*oddsum + 2*evensum + y[end] )/(3*(length(x)-1))
        return I
    else #i.e. we have an odd number of intervals
        # First we need to calculate the simpson 1/3 up until... end-3
        # since we have N-3 points for simpson 1/3 and N-4 intervals
        oddsum, evensum = 0,0 #initialize
        for  j in 2:2:(length(x) - 4)
            oddsum += y[j]
        end
        for j in 3:2:(length(x) - 4)
            evensum += y[j]
        end
        I1 = (x[end-3]-x[1])* ( y[1] + 4*oddsum + 2*evensum + y[end-3] )/(3*(length(x)-4))

        I2 = (x[end]-x[end-3])/8*( y[end-3] + 3*y[end-2] + 3*y[end-1] + y[end] )

        return I1 + I2

    end
end


text:
struct FashionMNIST <: GrayImage{MLImage} end
getModule(::FashionMNIST) = MLDatasets.FashionMNIST
size(::FashionMNIST) = (60000, 0, 10000)
problem(::FashionMNIST) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
@plottest begin
    # LIF params
    τm = 100
    vreset = 0.0
    vth = 0.1
    R = 1.75

    # Input spike train params
    rate = 0.05
    T = 1000

    lif = Neuron(QueuedSynapse(Synapse.Delta()), LIF(τm, vreset, R), Threshold.Ideal(vth))
    input = ConstantRate(rate)
    spikes = excite!(lif, input, T)

    # callback to record voltages
    voltages = Float64[]
    record = let lif = lif
        () -> push!(voltages, getvoltage(lif))
    end

    # compile
    simulate!(lif, T; cb = record)
    reset!(lif)
    voltages = Float64[]
    excite!(lif, spikes)

    # simulate
    simulate!(lif, T; cb = record, dense = true)

    # plot dense voltage recording
    plot(collect(1:T), voltages, title = "LIF Membrane Potential Over Time", xlabel = "Time (sec)", ylabel = "Potential (V)", label = "Dense")

    # repeat with dense simulation
    reset!(lif)
    voltages = Float64[]
    excite!(lif, spikes)
    output = simulate!(lif, T; cb = record)

    # plot sparse voltage recording
    voltage_plot = plot!(spikes, voltages,
        title = "LIF Membrane Potential Over Time", xlabel = "Time (sec)", ylabel = "Potential (V)", label = "Sparse (default)")

    # plot raster plot
    raster_plot = rasterplot(spikes, output, label = ["Input", "Output"], title = "Raster Plot", xlabel = "Time (sec)")

    plot(raster_plot, voltage_plot, layout = grid(2, 1))
end joinpath(datadir, "Lif-Test.png") !isci


text:
isdefined(Base, :__precompile__) && __precompile__()

module JOS3
    export clipdb, fold, mps
    include("./filters/filters.jl")
end