text:
using ValueOrientedRiskManagementInsurance

using Distributions
using DataFrames

include("ECModel_Input.jl")

println("Start ECModel ...")
#################################################################
"Create DataFrame from information on business units and total"
function getdf(bu, total)
  df =  DataFrame(BU = Array(Symbol, 0),
                  GrossNet = Array(Symbol, 0),
                  CumProb = Array(Float64, 0),
                  Profit = Array(Real, 0))

  for i = 1:length(bu)
    append!(df, DataFrame(
      BU = fill!(Array(Symbol, n_scen), bu[i].id),
      GrossNet = fill!(Array(Symbol, n_scen), :gross),
      CumProb = real(collect(1:n_scen)) /  n_scen,
      Profit = sort(bu[i].gross.profit)))
    append!(df, DataFrame(
      BU = fill!(Array(Symbol, n_scen), bu[i].id),
      GrossNet = fill!(Array(Symbol, n_scen), :net),
      CumProb = real(collect(1:n_scen)) /  n_scen,
      Profit = sort(bu[i].net.profit)))
  end
  append!(df, DataFrame(
    BU = fill(:total, n_scen),
    GrossNet = fill(:gross, n_scen),
    CumProb = real(collect(1:n_scen)) /  n_scen,
    Profit = sort(total.gross.profit)))
  append!(df, DataFrame(
    BU = fill(:total, n_scen),
    GrossNet = fill(:net, n_scen),
    CumProb = real(collect(1:n_scen)) /  n_scen,
    Profit = sort(total.net.profit)))
  append!(df, DataFrame(
    BU = fill(:risk_cap_total, 4),
    GrossNet = [:gross, :gross, :net, :net],
    CumProb = real([0., 1., 0., 1.]),
    Profit = convert(Vector{Real},
                     [fill(-total.gross.eco_cap, 2);
                      fill(-total.net.eco_cap,2)])))
  return df
end

"Construct efficient frontier"
function effline(x₀, y₀, a, y_min, y_max)
  x_min = x₀ + (y_min-y₀) / a
  x_max = x₀ + (y_max-y₀) / a
  return hcat([x_min, x₀, x_max], [y_min, y₀, y_max])
end

eff(x₀, y₀, a, x) = y₀ + a * (x-x₀)

"Equalizing risk premium"
function riskprice(ins_input::DataFrame,
                   inv_inout::DataFrame,
                   tau_kendall::Matrix{Real},
                   n_scen::Int,
                   α::Real,
                   s::Real,
                   costs_fixed::Real,
                   Δp_min::Real,
                   Δp_max::Real,
                   Δp_init::Real)
  ins = deepcopy(ins_input)
  i_max = 50 ## max 50 iterations
  rorac_gross_i = zeros(Real, i_max)
  rorac_gross_f = zeros(Real, i_max)
  Δp = zeros(Real, i_max + 1)
  Δp_tol  = 1e-5      ## break criterion : minimum price change
  row_l = ins[ins[:name] .== "Liability", :ctr][1,1]
  row_f = ins[ins[:name] .== "Fire", :ctr][1,1]
  if (row_l != ins[row_l, :ctr]) | (row_f != ins[row_f, :ctr])
    error("row_l ($row_l) or row_f ($row_f) does not match.")
  end
  n_steps = 0
  Δp[1] = Δp_init
  for i = 1:i_max
    srand(seed)
    ins[row_l, :premium] =
      ins_input[row_l, :premium] * (1 + Δp[i])
    ins[row_f, :premium] =
      ins_input[row_f, :premium] -
      ins_input[row_l, :premium] * Δp[i]
    fac_i = ins_input[row_l, :premium] / ins[row_l, :premium]
    fac_f = ins_input[row_f, :premium] / ins[row_f, :premium]
    ins[row_l, :loss_ratio] =
      ins_input[row_l, :loss_ratio] * fac_i
    ins[row_f, :loss_ratio] =
      ins_input[row_f, :loss_ratio] * fac_f
    ins[row_l, :cost_ratio] =
      ins_input[row_l, :cost_ratio] * fac_i
    ins[row_f, :cost_ratio] =
      ins_input[row_f, :cost_ratio] * fac_f
    ins[row_l, :re_costs] = ins_input[row_l, :re_costs] * fac_i
    ins[row_f, :re_costs] = ins_input[row_f, :re_costs] * fac_f

    bu, total =
      project(ins, inv_input, tau_kendall, n_scen, α, s,
              costs_fixed)
    rorac_gross_i[i] = bu[row_l].gross.rorac
    rorac_gross_f[i] = bu[row_f].gross.rorac

    if Δp_max - Δp_min < Δp_tol
      n_steps = i
      break
    end
    if rorac_gross_i[i] < rorac_gross_f[i]
      Δp_min = Δp[i]
    else
      Δp_max = Δp[i]
    end
    Δp[i+1] = (Δp_min + Δp_max) / 2
  end
  return ins, Δp[1:n_steps],
         rorac_gross_i[1:n_steps], rorac_gross_f[1:n_steps]
end

"Optimizing Capitalization"
function optcap(ins_input::DataFrame,
                inv_input::DataFrame,
                tau_kendall::Matrix{Real},
                n_scen::Int,
                α::Real,
                s::Real,
                cost_fixed::Real)
  inv = deepcopy(inv_input)
  i_max = 8
  invest_init = inv[1,:init] - collect(0:(i_max-1)) * 100.
  profit_net = zeros(Real, i_max)
  ec_net = zeros(Real, i_max)
  rorac_net = zeros(Real, i_max)
  roc_net = zeros(Real, i_max)

  for i in 1:i_max
    srand(seed)
    inv[1,:init] = invest_init[i]
    bu, total =
      project(ins_input, inv, tau_kendall,
              n_scen, α, s, cost_fixed)
    profit_net[i] = total.net.profit_mean
    ec_net[i] = total.net.eco_cap
    rorac_net[i] = total.net.rorac
    roc_net[i] = total.net.profit_mean / invest_init[i]
  end
  inv_init_opt = 0
  ec_net_opt = 0
  profit_net_opt = 0
  i_opt = 0
  ## only works if invest_init[1] >  ec_cap_net[1]
  for i in 1:i_max
    if invest_init[i] < ec_net[i]
      inv_init_opt = invest_init[i-1]
      ec_net_opt = ec_net[i-1]
      profit_net_opt = profit_net[i-1]
      i_opt = i-1
      break
    end
  end
  inv[1, :init] = inv_init_opt
  return inv, i_opt, invest_init, ec_net, rorac_net, roc_net,
  inv_init_opt, ec_net_opt, profit_net_opt
end

"Optimizing Fire Reinsurance"
function optrefire(ins_input::DataFrame,
                   inv_input::DataFrame,
                   tau_kendall::Matrix{Real},
                   n_scen::Int,
                   α::Real,
                   s::Real,
                   cost_fixed::Real)
  ins = deepcopy(ins_input)
  n_points = 21
  points = zeros(Real, n_points+1, 4)

  i_opt = 0
  ceded_opt = 0
  profit_net_opt = 0
  ec_net_opt = 0
  rorac_net_opt = 0
  ceded = zeros(Real, n_points)
  profit_net = zeros(Real, n_points)
  ec_net = zeros(Real, n_points)
  rorac_net = zeros(Real, n_points)
  row_f = ins[ins[:name] .== "Fire", :ctr][1,1]
  if row_f != ins[row_f, :ctr]
    error("row_f ($row_f) does not match.")
  end

  for i in 1:n_points
    srand(seed)
    ceded[i] = (i-1) / (n_points-1)
    ins[row_f, :re_ceded] = ceded[i]
    bu, total =
      project(ins, inv_input, tau_kendall,
      n_scen, α, s, cost_fixed)
    ec_net[i] = total.net.eco_cap
    profit_net[i] = total.net.profit_mean
    rorac_net[i] = total.net.rorac
    if rorac_net[i] > rorac_net_opt
      ceded_opt = ceded[i]
      profit_net_opt = profit_net[i]
      ec_net_opt = ec_net[i]
      rorac_net_opt = rorac_net[i]
      i_opt = i
    end
  end
  ins[row_f, :re_ceded] = ceded_opt
  return ins, ceded, profit_net, ec_net,
  i_opt, ceded_opt, profit_net_opt, ec_net_opt, rorac_net_opt
end

"Optimizing Reinsurance RAROC no constraints"
function optreraroc(ins_input::DataFrame,
                    inv_input::DataFrame,
                    tau_kendall::Matrix{Real},
                    n_scen::Int,
                    α::Real,
                    s::Real,
                    cost_fixed::Real,
                    n_points::Int)
  ins = deepcopy(ins_input)
  unif = Uniform()
  ceded = Array(Real, nrow(ins), n_points)
  profit_net = Array(Real, n_points)
  ec_net = Array(Real, n_points)
  rorac_net = Array(Real, n_points)
  i_opt, profit_net_opt, ec_net_opt, rorac_net_opt =
    1, 0.0, 0.0, 0.0

  for b in 1:nrow(ins)
    ceded[b,:] = rand(unif, n_points)
  end
  for i = 1:n_points
    srand(seed)
    for b = 1:nrow(ins)
      ins[b, :re_ceded] = ceded[b, i]
    end
    bu, total =
      project(ins, inv_input, tau_kendall,
              n_scen, α, s, cost_fixed)
    profit_net[i] = total.net.profit_mean
    ec_net[i] = total.net.eco_cap
    rorac_net[i] = total.net.rorac
    if rorac_net[i] > rorac_net_opt
      i_opt = i
      profit_net_opt = profit_net[i]
      ec_net_opt = ec_net[i]
      rorac_net_opt = rorac_net[i]
    end
  end
  for b = 1:nrow(ins)
    ins[b, :re_ceded] = ceded[b, i_opt]
  end
  return ins, profit_net, ec_net, rorac_net,
  profit_net_opt, ec_net_opt, rorac_net_opt
end

"Optimizing Reinsurance EVA no constraints"
function optreeva(ins_input::DataFrame,
                  inv_input::DataFrame,
                  tau_kendall::Matrix{Real},
                  n_scen::Int,
                  α::Real,
                  s::Real,
                  cost_fixed::Real,
                  n_points::Int,
                  hurdle::Real)
  ins = deepcopy(ins_input)
  unif = Uniform()
  ceded = Array(Real, nrow(ins), n_points)
  avg_ceded = Array(Real, n_points)
  eva_net = Array(Real, n_points)
  i_opt = 0
  avg_ceded_opt  = 0.0
  eva_net_opt = 0.0

  for b in 1:nrow(ins)
    ceded[b,:] = rand(unif, n_points)
  end
  for i = 1:n_points
    srand(seed)
    for b = 1:nrow(ins)
      ins[b, :re_ceded] = ceded[b, i]
    end
    bu, total =
      project(ins, inv_input, tau_kendall,
              n_scen, α, s, cost_fixed)
    avg_ceded[i] =
      ins[:re_ceded] ⋅ ins[:premium] /sum(ins[:premium])
    eva_net[i] =
      total.net.profit_mean - hurdle * total.net.eco_cap
    if eva_net[i] > eva_net_opt
      i_opt = i
      avg_ceded_opt = avg_ceded[i]
      eva_net_opt = eva_net[i]
    end
  end
  for b = 1:nrow(ins)
    ins[b, :re_ceded] = ceded[b, i_opt]
  end
  return ins, avg_ceded, eva_net, avg_ceded_opt, eva_net_opt
end

"Optim. reins. RAROC prescribed average reinsurance quote"
function optre(ins_input::DataFrame,
               inv_input::DataFrame,
               tau_kendall::Matrix{Real},
               n_scen::Int,
               α::Real,
               s::Real,
               cost_fixed::Real,
               n_points::Int,
               avg_ceded::Real)
  ins = deepcopy(ins_input)
  unif = Uniform()
  ceded = Array(Real, nrow(ins), n_points)
  random =  Array(Real, nrow(ins), n_points)
  profit_net = Array(Real, n_points)
  ec_net = Array(Real, n_points)
  rorac_net = Array(Real, n_points)
  i_opt, profit_net_opt, ec_net_opt, rorac_net_opt =
    0, 0.0, 0.0, 0.0

  prem_gross = sum(ins_input[:, :premium])
  for i in 1:n_points
    rand_vector_admissable = false
    while !rand_vector_admissable
      random[:,i] = rand(unif, nrow(ins))
      prem_temp = random[:,i] ⋅ ins_input[:, :premium]
      for b in 1:nrow(ins)
        ceded[b, i] =
          avg_ceded * random[b, i] * prem_gross / prem_temp
      end
      # We cannot cede more than 100%
      rand_vector_admissable =
        all(ceded[:,i] .<= ones(nrow(ins)))
    end
  end
  for i in 1:n_points
    srand(seed)
    # for b = 1:nrow(ins)
    #   ceded[b, i] =
    #     avg_ceded * random[b, i] * prem_gross / prem_temp
    #   ins[b, :re_ceded] = ceded[b, i]
    # end
    ins[:, :re_ceded] = ceded[:, i]
    bu, total =
      project(ins, inv_input, tau_kendall,
              n_scen, α, s, cost_fixed)
    profit_net[i] = total.net.profit_mean
    ec_net[i] = total.net.eco_cap
    rorac_net[i] = total.net.rorac
    if rorac_net[i] > rorac_net_opt
      i_opt = i
      profit_net_opt = profit_net[i]
      ec_net_opt = ec_net[i]
      rorac_net_opt = rorac_net[i]
    end
  end
  ins[:, :re_ceded] = ceded[:, i_opt]
  return ins, ceded, profit_net, ec_net, rorac_net,
  profit_net_opt, ec_net_opt, rorac_net_opt
end

#################################################################

row_f =
  insurance_input[insurance_input[:name] .== "Fire", :ctr][1,1]
row_l =
  insurance_input[insurance_input[:name] .==
  "Liability", :ctr][1,1]
row_t =
  insurance_input[insurance_input[:name] .== "Theft", :ctr][1,1]

pl_fire_exact_gross =
  (1 + s -
     insurance_input[row_f, :loss_ratio] -
     insurance_input[row_f, :cost_ratio]) *
  insurance_input[row_f, :premium]

pl_cap_rf = s * invest_input[1, :init]
pl_diff = pl_cap_rf - costs_fixed

srand(seed) ## fix random seed for repeatable results
bu, total =
  project(insurance_input, invest_input, tau_kendall,
          n_scen, α, s, costs_fixed)

n_bu = length(bu)
abbrevs = AbstractString[string(bu[b].name[1]) for b = 1:n_bu]

diff_ec_gross =
  total.gross.eco_cap -
  sum([bu[i].gross.eco_cap  for i = 1:length(bu)])
diff_ec_net =
  total.net.eco_cap -
  sum([bu[i].net.eco_cap for i = 1:length(bu)])

## BU results
profit_gross = Real[bu[i].gross.profit_mean for i = 1:length(bu)]
profit_net = Real[bu[i].net.profit_mean for i = 1:length(bu)]

ec_gross = Real[bu[i].gross.eco_cap for i = 1:length(bu)]
ec_net = Real[bu[i].net.eco_cap for i = 1:length(bu)]

rorac_gross = Real[bu[i].gross.rorac for i = 1:length(bu)]
rorac_net = Real[bu[i].net.rorac for i = 1:length(bu)]

x = hcat(profit_gross, ec_gross, ec_net)

##  Portfolio Optimizations #####################################

ins_input = deepcopy(insurance_input)
inv_input = deepcopy(invest_input)

## Risk-Adjusted Price-Setting ==================================
println("Risk-Adjusted Price-Setting ...")
Δp_min = 0.0
Δp_max = 0.05
Δp_init = 0.0

ins_input_rp, Δp, rorac_price_gross_i, rorac_price_gross_f =
  riskprice(ins_input, inv_input, tau_kendall, n_scen, α, s,
            costs_fixed, Δp_min, Δp_max, Δp_init)

srand(seed)
bu_rp, total_rp =
  project(ins_input_rp, invest_input, tau_kendall,
          n_scen, α, s, costs_fixed)

## Optimizing Capitalization ====================================
println("Optimizing Capitalization ... ")
inv_input_oc, i_oc_opt, invest_init_oc, ec_net_oc, rorac_net_oc,
roc_net_oc, inv_init_oc_opt, ec_net_oc_opt, profit_net_oc_opt =
  optcap(ins_input_rp, inv_input, tau_kendall,
         n_scen, α, s, costs_fixed)

srand(seed)
bu_oc, total_oc =
  project(ins_input_rp, inv_input_oc, tau_kendall,
          n_scen, α, s, costs_fixed)

## Optimizing Fire Reinsurance ==================================
println("Optimizing Fire Reinsurance ...")

ins_input_ofr, ceded_ofr, profit_net_ofr, ec_net_ofr,
i_ofr, ceded_ofr_opt, profit_net_ofr_opt, ec_net_ofr_opt,
rorac_net_ofr_opt =
  optrefire(ins_input_rp, inv_input_oc, tau_kendall, n_scen, α,
            s, costs_fixed)

## Optimizing Reinsurance =======================================

## Optimizing Reinsurance RAROC no constraints ==================
println("Optimizing Reinsurance RAROC no constraints ... ")

ins_input_orraroc, profit_net_orraroc, ec_net_orraroc,
rorac_net_orraroc, profit_net_orraroc_opt, ec_net_orraroc_opt,
rorac_net_orraroc_opt =
  optreraroc(ins_input_ofr, inv_input_oc, tau_kendall,
             n_scen, α, s, costs_fixed, n_cloud)

ec_min = minimum(ec_net_orraroc)
p_min = findin(ec_net_orraroc, ec_min)

## Optimizing Reinsurance EVA no constraints ====================
println("Optimizing Reinsurance EVA no constraints ... ")

ins_input_oreva, avg_ceded_oreva, eva_net_oreva,
avg_ceded_opt_oreva, eva_net_opt_oreva =
  optreeva(ins_input_ofr, inv_input_oc, tau_kendall, n_scen,
           α,  s, costs_fixed, n_cloud, hurdle)

avg_ceded_ofr =
  ins_input_ofr[:re_ceded] ⋅ ins_input_ofr[:premium] /
  sum(ins_input_ofr[:premium])
eva_net_ofr = profit_net_ofr_opt - hurdle * ec_net_ofr_opt

## Optim. Reins. RAROC prescribed average reinsurance quote =====
println("Optim. Reins. RAROC prescr. average reins. quote ...")

ins_input_or, ceded_or, profit_net_or, ec_net_or, rorac_net_or,
profit_net_or_opt, ec_net_or_opt, rorac_net_or_opt =
  optre(ins_input_ofr, inv_input_oc, tau_kendall, n_scen,
        α,  s, costs_fixed, n_cloud, avg_ceded_ofr)

srand(seed)
bu_oc, total_oc =
  project(ins_input_or, inv_input_oc, tau_kendall,
          n_scen, α, s, costs_fixed)

## Other levels of confidence ===================================
println("Other levels of confidence ...")

srand(seed)
bu_final, total_final =
  project(ins_input_or, inv_input_oc, tau_kendall,
          n_scen, α, s, costs_fixed)

var_α = [0.5, 0.75, 0.9, 0.99]
var_eco_cap = [es(-total_final.net.profit, β) for β in var_α]

print("Mean total gross profit:      ")
println(round(total.gross.profit_mean, 2))
print("Mean total net profit:        ")
println(round(total.net.profit_mean, 2))
print("Total gross economic capital: ")
println(round(total.gross.eco_cap, 2))
print("Total net economic capital:   ")
println(round(total.net.eco_cap, 2))
print("Total gross RORAC:            ")
println(round(total.gross.rorac, 2))
print("Total net RORAC:              ")
println(round(total.net.rorac, 2))

print("Risk adj. pricing, fire loss ratio     : ")
println(round(ins_input_rp[1, :loss_ratio], 4))
print("Risk adj. pricing, liability loss ratio: ")
println(round(ins_input_rp[2, :loss_ratio], 4))
print("Risk adj. pricing, theft loss ratio:     ")
println(round(ins_input_rp[3, :loss_ratio], 4))

print("Optimal capital, net RORAC:  ")
println(round(ins_input_rp[3, :loss_ratio], 4))

print("Optimal product mix with quota constraint, net RORAC:  ")
println(round(rorac_net_or_opt, 4))

#################################################################
println("End ECModel \n")


text:
using ManifoldLearning
import ManifoldLearning: BruteForce, knn, swiss_roll
using Test
import Random

Random.seed!(3483743871)

@testset "Nearest Neighbors" begin

    # setup parameters
    k = 12
    X, L = swiss_roll(100)
    DD, EE = knn(X,k)
    @test_throws AssertionError knn(rand(3,10), k)

    NN = fit(BruteForce, X, k)
    D, E = knn(NN, X)
    @test size(X,2) == size(D,2) && size(D, 1) == k
    @test size(X,2) == size(E,2) && size(E, 1) == k
    @test E == EE
    @test D ≈ DD

    D, E = knn(NN, X[:,1:k+1])
    @test k+1 == size(D,2) && size(D, 1) == k
    @test k+1 == size(E,2) && size(E, 1) == k

    D, E = knn(NN, X[:,1:k+1], self=true)
    @test D[1,:] == zeros(k+1)
    @test E[1,:] == collect(1:k+1)

    @test_throws AssertionError knn(NN, X[:,1:k])
end

@testset "Manifold Learning" begin

    # setup parameters
    k = 12
    d = 2
    X, L = swiss_roll()

    # test algorithms
    @testset for algorithm in [Isomap, LEM, LLE, HLLE, LTSA, DiffMap]
        for (k, T) in zip([5, 12], [Float32, Float64])
            # construct KW parameters
            kwargs = [:maxoutdim=>d]
            if algorithm == DiffMap
                push!(kwargs, :t => k)
            else
                push!(kwargs, :k => k)
            end

            # call transformation
            Y = fit(algorithm, convert(Array{T}, X); kwargs...)

            # test results
            @test outdim(Y) == d
            @test size(transform(Y), 2) == size(X, 2)
            @test length(split(sprint(show, Y), '\n')) > 1
            @test length(eigvals(Y)) == d
            if algorithm !== DiffMap
                @test neighbors(Y) == k
                @test length(vertices(Y)) > 1
            end
        end
    end

end

@testset "OOS" begin
    k = 10
    d = 2
    X, _ = swiss_roll()
    M = fit(Isomap, X; k=k, maxoutdim=d)

    @test all(sum(abs2, transform(M) .- transform(M,X), dims=1) .< eps())
end

ext: jl


text:
using HTTP

"""
    register_client!(cm::ClientManager, req::HTTP.Request)

Register a new client in the `ClientManager`. The body of the request should
contain the client URL: `http://hostname:port`.

!!! info

    **API Endpoint**: this is an HTTP server endpoint.
"""
function register_client!(cm::ClientManager, req::HTTP.Request)::HTTP.Response
    client_url = String(req.body)
    register!(cm, client_url)
    @info "client registered: $client_url"
    return HTTP.Response(200, "Client registered successfully")
end


text:
using Test
using WavePropBase
using WavePropBase.Geometry

# Two points
clear_entities!()
p1 = ElementaryEntity(0, 1, ElementaryEntity[])
p2 = ElementaryEntity(0, 2, ElementaryEntity[])
@test p1 == p1
@test p1 != p2
points = [p1, p2]
# Three lines
l1 = ElementaryEntity(1, 1, ElementaryEntity[p1, p2])
l2 = ElementaryEntity(1, 2, ElementaryEntity[p1, p2])
l3 = ElementaryEntity(1, 3, ElementaryEntity[p1, p2])
lines = [l1, l2, l3]
# Two surfaces
s1 = ElementaryEntity(2, 1, ElementaryEntity[l1, l2])
s2 = ElementaryEntity(2, 2, ElementaryEntity[l2, l3])
surfaces = [s1, s2]

@testset "ElementaryEntity" begin
    @test geometric_dimension(p1) == 0
    @test tag(p1) == 1
    @test boundary(l1) == points
    for e in vcat(points, lines, surfaces)
        for b in boundary(e)
            @test geometric_dimension(e)-1 == geometric_dimension(b)
        end
    end
end

# Domains
Ω1 = Domain(s1)
Ω2 = Domain(s2)
Ω  = Domain(surfaces)

@testset "Domain" begin
    @test entities(Ω) == surfaces
    @test length(Ω) == 2
    @test !isempty(Ω)
    @test Ω2 == setdiff(Ω, Ω1)
    @test Ω  == union(Ω1, Ω2)
    @test s1 in Ω
    @test l1 in Ω
    @test Ω[1] == s1
    @test Ω[end] == s2
    for ω in Ω
        for Γ in boundary(ω)
            @test geometric_dimension(ω)-1 == geometric_dimension(Γ)
        end
    end
    @test Ω1 == intersect(Ω1, Ω)
    @test Domain(l2) == intersect(Ω1, Ω2)
    @test issubset(Ω1, Ω)
    @test Ω2 == setdiff(Ω, Ω1)
    @test union(Domain.(lines)...) == skeleton(Ω)
    @test Domain(l2) == internal_boundary(Ω)
    @test union(Domain(l1), Domain(l3)) == external_boundary(Ω)
    @test keys(Ω, 1) == [(1,1), (1,2), (1,3)]
    @test keys(Ω) == [(2,1), (2,2)]
    @test Domain([s1,s2]) == Domain([s2,s1])
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


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

mutable struct VirtualMachineScaleSetIPConfigurationProperties <: SwaggerModel
    subnet::Any # spec type: Union{ Nothing, ApiEntityReference } # spec name: subnet
    primary::Any # spec type: Union{ Nothing, Bool } # spec name: primary
    publicIPAddressConfiguration::Any # spec type: Union{ Nothing, VirtualMachineScaleSetPublicIPAddressConfiguration } # spec name: publicIPAddressConfiguration
    privateIPAddressVersion::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddressVersion
    applicationGatewayBackendAddressPools::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: applicationGatewayBackendAddressPools
    applicationSecurityGroups::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: applicationSecurityGroups
    loadBalancerBackendAddressPools::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: loadBalancerBackendAddressPools
    loadBalancerInboundNatPools::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: loadBalancerInboundNatPools

    function VirtualMachineScaleSetIPConfigurationProperties(;subnet=nothing, primary=nothing, publicIPAddressConfiguration=nothing, privateIPAddressVersion=nothing, applicationGatewayBackendAddressPools=nothing, applicationSecurityGroups=nothing, loadBalancerBackendAddressPools=nothing, loadBalancerInboundNatPools=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("primary"), primary)
        setfield!(o, Symbol("primary"), primary)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("publicIPAddressConfiguration"), publicIPAddressConfiguration)
        setfield!(o, Symbol("publicIPAddressConfiguration"), publicIPAddressConfiguration)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("privateIPAddressVersion"), privateIPAddressVersion)
        setfield!(o, Symbol("privateIPAddressVersion"), privateIPAddressVersion)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("applicationGatewayBackendAddressPools"), applicationGatewayBackendAddressPools)
        setfield!(o, Symbol("applicationGatewayBackendAddressPools"), applicationGatewayBackendAddressPools)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("applicationSecurityGroups"), applicationSecurityGroups)
        setfield!(o, Symbol("applicationSecurityGroups"), applicationSecurityGroups)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("loadBalancerBackendAddressPools"), loadBalancerBackendAddressPools)
        setfield!(o, Symbol("loadBalancerBackendAddressPools"), loadBalancerBackendAddressPools)
        validate_property(VirtualMachineScaleSetIPConfigurationProperties, Symbol("loadBalancerInboundNatPools"), loadBalancerInboundNatPools)
        setfield!(o, Symbol("loadBalancerInboundNatPools"), loadBalancerInboundNatPools)
        o
    end
end # type VirtualMachineScaleSetIPConfigurationProperties

const _property_map_VirtualMachineScaleSetIPConfigurationProperties = Dict{Symbol,Symbol}(Symbol("subnet")=>Symbol("subnet"), Symbol("primary")=>Symbol("primary"), Symbol("publicIPAddressConfiguration")=>Symbol("publicIPAddressConfiguration"), Symbol("privateIPAddressVersion")=>Symbol("privateIPAddressVersion"), Symbol("applicationGatewayBackendAddressPools")=>Symbol("applicationGatewayBackendAddressPools"), Symbol("applicationSecurityGroups")=>Symbol("applicationSecurityGroups"), Symbol("loadBalancerBackendAddressPools")=>Symbol("loadBalancerBackendAddressPools"), Symbol("loadBalancerInboundNatPools")=>Symbol("loadBalancerInboundNatPools"))
const _property_types_VirtualMachineScaleSetIPConfigurationProperties = Dict{Symbol,String}(Symbol("subnet")=>"ApiEntityReference", Symbol("primary")=>"Bool", Symbol("publicIPAddressConfiguration")=>"VirtualMachineScaleSetPublicIPAddressConfiguration", Symbol("privateIPAddressVersion")=>"String", Symbol("applicationGatewayBackendAddressPools")=>"Vector{SubResource}", Symbol("applicationSecurityGroups")=>"Vector{SubResource}", Symbol("loadBalancerBackendAddressPools")=>"Vector{SubResource}", Symbol("loadBalancerInboundNatPools")=>"Vector{SubResource}")
Base.propertynames(::Type{ VirtualMachineScaleSetIPConfigurationProperties }) = collect(keys(_property_map_VirtualMachineScaleSetIPConfigurationProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetIPConfigurationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetIPConfigurationProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetIPConfigurationProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetIPConfigurationProperties[property_name]

const _allowed_VirtualMachineScaleSetIPConfigurationProperties_privateIPAddressVersion = ["IPv4", "IPv6"]

function check_required(o::VirtualMachineScaleSetIPConfigurationProperties)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetIPConfigurationProperties }, name::Symbol, val)
    if name === Symbol("privateIPAddressVersion")
        Swagger.validate_param(name, "VirtualMachineScaleSetIPConfigurationProperties", :enum, val, _allowed_VirtualMachineScaleSetIPConfigurationProperties_privateIPAddressVersion)
    end
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, GrammianAngularFields

makedocs(sitename="Grammian Angular Fields Documentation")
deploydocs(repo = "github.com/DoktorMike/GrammianAngularFields.jl.git")


text:
module HWunconstrained

	# imports: which packages are we going to use in this module?
	using Distributions, Optim, Plots, DataFrames
	using Random
	using Statistics
	using LinearAlgebra

    export maximize_like_grad, makeData

	# methods/functions
	# -----------------

	# data creator
	# should/could return a dict with beta,numobs,X,y,norm)
	# true coeff vector, number of obs, data matrix X, response vector y, and a type of parametric distribution for G.
	function makeData(n=10_000)

 		# your turn!

		return Dict("beta"=>beta,"n"=>numobs,"X"=>X,"y"=>y,"dist"=>norm)
	end

	# log likelihood function at x
	function loglik(betas::Vector,d::Dict)

	end

	# gradient of the likelihood at x
	function grad!(storage::Vector,betas::Vector,d)

	end

	# hessian of the likelihood at x
	function hessian!(storage::Matrix,betas::Vector,d)

	end

	function info_mat(betas::Vector,d)

	end

	function inv_Info(betas::Vector,d)

	end

	"""
	inverse of observed information matrix
	"""
	function inv_observedInfo(betas::Vector,d)

	end

	"""
	standard errors
	"""
	function se(betas::Vector,d::Dict)

	end

	# function that maximizes the log likelihood without the gradient
	# with a call to `optimize` and returns the result
	function maximize_like(x0=[0.8,1.0,-0.1],meth=NelderMead())

	end
	function maximize_like_helpNM(x0=[ 1; 1.5; -0.5 ],meth=NelderMead())

	end

	# function that maximizes the log likelihood with the gradient
	# with a call to `optimize` and returns the result
	function maximize_like_grad(x0=[0.8,1.0,-0.1],meth=BFGS())

	end

	function maximize_like_grad_hess(x0=[0.8,1.0,-0.1],meth=Newton())

	end

	function maximize_like_grad_se(x0=[0.8,1.0,-0.1],meth=BFGS())

	end

	# visual diagnostics
	# ------------------

	# function that plots the likelihood
	# we are looking for a figure with 3 subplots, where each subplot
	# varies one of the parameters, holding the others fixed at the true value
	# we want to see whether there is a global minimum of the likelihood at the the true value
	function plotLike()

	end
	function plotGrad()

	end

	function saveplots()
		p1 = plotLike()
		savefig(p1,joinpath(dirname(@__FILE__),"..","likelihood.png"))
		p2 = plotGrad()
		savefig(p2,joinpath(dirname(@__FILE__),"..","gradient.png"))
		@info("saved both plots.")
	end

end


text:
using MesoscaleML
using Test

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


text:
# Composition
# ===========
#
# Sequence composition counter.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

"""
Sequence composition.

This is a subtype of `Associative{T,Int}`, and the `getindex` method returns the
number of occurrences of a symbol or a k-mer.
"""
struct Composition{T} <: AbstractDict{T,Int}
    counts::Dict{T,Int}
end

function Composition(seq::BioSequence{A}) where A <: NucAlphs
    counts = zeros(Int, 16)
    @inbounds for x in seq
        counts[reinterpret(UInt8, x) + 1] += 1
    end
    return Composition{eltype(A)}(count_array2dict(counts, alphabet(A)))
end

function Composition(seq::ReferenceSequence)
    counts = zeros(Int, 16)
    @inbounds for x in seq
        counts[reinterpret(UInt8, x) + 1] += 1
    end
    return Composition{DNA}(count_array2dict(counts, ACGTN))
end

function Composition(kmer::DNAKmer)
    counts = Dict{DNA,Int}()
    counts[DNA_A] = count_a(kmer)
    counts[DNA_C] = count_c(kmer)
    counts[DNA_G] = count_g(kmer)
    counts[DNA_T] = count_t(kmer)
    return Composition(counts)
end

function Composition(kmer::RNAKmer)
    counts = Dict{RNA,Int}()
    counts[RNA_A] = count_a(kmer)
    counts[RNA_C] = count_c(kmer)
    counts[RNA_G] = count_g(kmer)
    counts[RNA_U] = count_t(kmer)
    return Composition(counts)
end

function Composition(seq::AminoAcidSequence)
    counts = zeros(Int, length(alphabet(AminoAcid)))
    @inbounds for x in seq
        counts[reinterpret(UInt8, x) + 1] += 1
    end
    return Composition{AminoAcid}(count_array2dict(counts, alphabet(AminoAcid)))
end

function Composition(iter::EachKmerIterator{T}) where {T<:Kmer}
    counts = Dict{T,Int}()
    if kmersize(T) ≤ 8
        # This is faster for short k-mers.
        counts′ = zeros(Int, 4^kmersize(T))
        for (_, x) in iter
            @inbounds counts′[reinterpret(Int, x)+1] += 1
        end
        for x in 1:lastindex(counts′)
            @inbounds c = counts′[x]
            if c > 0
                counts[reinterpret(T, x-1)] = c
            end
        end
    else
        for (_, x) in iter
            get!(counts, x, 0)
            counts[x] += 1
        end
    end
    return Composition{T}(counts)
end

"""
    composition(seq | kmer_iter)

Calculate composition of biological symbols in `seq` or k-mers in `kmer_iter`.
"""
function composition(iter::Union{Sequence,EachKmerIterator})
    return Composition(iter)
end

"""
    composition(iter)

A generalised composition algorithm, which computes the number of unique items
produced by an iterable.

# Example

```jlcon

# Example, counting unique sequences.

julia> a = dna"AAAAAAAATTTTTT"
14nt DNA Sequence:
AAAAAAAATTTTTT

julia> b = dna"AAAAAAAATTTTTT"
14nt DNA Sequence:
AAAAAAAATTTTTT

julia> c = a[5:10]
6nt DNA Sequence:
AAAATT

julia> composition([a, b, c])
Vector{BioSequences.BioSequence{BioSequences.DNAAlphabet{4}}} Composition:
  AAAATT         => 1
  AAAAAAAATTTTTT => 2
```
"""
function composition(iter)
    counts = Dict{eltype(iter), Int}()
    @inbounds for item in iter
        counts[item] = get(counts, item, 0) + 1
    end
    return Composition(counts)
end

function Base.:(==)(x::Composition{T}, y::Composition{T}) where T
    return x.counts == y.counts
end

function Base.length(comp::Composition)
    return length(comp.counts)
end

Base.iterate(comp::Composition) = iterate(comp.counts)
Base.iterate(comp::Composition, state) = iterate(comp.counts, state)

function Base.getindex(comp::Composition{T}, x) where {T}
    return get(comp.counts, convert(T, x), 0)
end

function Base.copy(comp::Composition)
    return Composition(copy(comp.counts))
end

function Base.merge(comp::Composition{T}, other::Composition{T}) where T
    return merge!(copy(comp), other)
end

function Base.merge!(comp::Composition{T}, other::Composition{T}) where T
    for (x, c) in other
        comp.counts[x] = comp[x] + c
    end
    return comp
end

function Base.summary(::Composition{T}) where T
    if T == DNA
        return "DNA Composition"
    elseif T == RNA
        return "RNA Composition"
    elseif T == AminoAcid
        return "Amino Acid Composition"
    else
        return string(T, " Composition")
    end
end

function count_array2dict(counts, alphabet)
    counts′ = Dict{eltype(alphabet),Int}()
    sizehint!(counts′, count(x -> !iszero(x), counts))
    for x in alphabet
        @inbounds c = counts[convert(Int64,x)+1]
        if c > 0
            counts′[x] = c
        end
    end
    return counts′
end


text:
using PkgBenchmark
using Reactive

Reactive.stop()

@benchgroup "Signal creation" begin
    @bench "int" Signal(0)
    @bench "string" Signal($("x"))
    @bench "standalone update" (push!($(Signal(0)), 1); Reactive.run_till_now())
end

@benchgroup "map 1" begin
    x = Signal(0)
    y = map(-, x)
    @bench "2 node" (push!($x, 1); Reactive.run_till_now())

    z = map(+, x)

    @bench "3 nodes" (push!($x, 1); Reactive.run_till_now())

    a = map(+, x, y)
    @bench "4 nodes" (push!($x, 1); Reactive.run_till_now())
end

@benchgroup "map 2" begin
    a = Signal(0.0)

    b = map(/, a, Signal(23.0))
    c = map(/, a, Signal(8.0))
    f = foldp(+, 0.0, b)

    d = map(Vec3f0, b)
    e = map(Vec3f0, c)
    g = map(Vec3f0, f)

    m = map(translationmatrix, d)
    m2 = map(translationmatrix, e)

    m3 = map(*, m, m2)
    # I don't know why, but Mat*Vec is broken right now
    result = map(m3, g) do a, b
         r = a * Vec4f0(b, 1)
         Vec3f0(r[1], r[2], r[3])
    end
    @bench "10 nodes" (push!($x, 1); Reactive.run_till_now())
end


text:
# Henon gridding for dimension explanation
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

fig = figure()
ax = subplot(1,2,1)
ax2 = subplot(1,2,2)
he = Systems.henon()
N = 10000
tr = trajectory(he, N, Ttr = 100)
es = [0.2, 0.05, 0.01]
ns = zero(es)
ax.grid(false)
mini = minima(tr)

rectcols = ["C0", "C1", "C3"]

for (i, e) ∈ enumerate(es)
    if i ≤ 2
        for x in mini[1]:e:1.4
            ax.axvline(x; color = "k", lw = 0.5, alpha = 1/(3^i))
        end
        for y in mini[2]:e:0.4
            ax.axhline(y; color = "k", lw = 0.5, alpha = 1/(3^i))
        end
    end
    p, bins = binhist(tr, e)
    for b in bins
        r = matplotlib.patches.Rectangle(b, e, e; color = rectcols[i], ec = "k", lw = 1/i)
        ax.add_artist(r)
    end
    ns[i] = length(probabilities(tr, e))
end

ax.plot(tr[:, 1], tr[:, 2], ls = "None", marker = ".", color = COLORS[3], ms = 0.75, zorder = 99)
ax.set_yticks([])
ax.set_ylabel("\$y\$")
ax.set_xticks([])
ax.set_xlabel("\$x\$")

ax2.scatter(log.(1 ./ es), log.(ns); c = rectcols, s = 200, zorder = 99)
s = linreg(log.(1 ./ es), log.(ns))[2]
ax2.plot(log.(1 ./ es), log.(1 ./ es) .* s .+ 2, color = "C2")
ax2.text(3.5, 6, "\$\\Delta = $(rdspl(s))\$"; color = "C2")
ax2.set_xlabel("\$\\log ( 1/\\varepsilon)\$"; labelpad = -15)
ax2.set_ylabel("\$\\log ( M)\$")
ax2.set_xticks([2, 3, 4])
ax2.set_xticklabels([2, "", 4])
fig.tight_layout(pad = 0.3)
fig.subplots_adjust(wspace = 0.2)
wsave(plotsdir("5", "henon_gridding"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module ParallelRadixSort

export MSDRadixSort, ParallelMSDRadixSort, ParallelStableMSDRadixSort, StableMSDRadixSort

# macro inbounds(ex)
#     esc(ex)
# end

using ArgCheck: @check
using Base.Threads: @spawn
using Base.Order: Order, Ordering, ForwardOrdering, ReverseOrdering
using ConstructionBase: setproperties
using StaticArrays: MVector
using StaticNumbers: @stat, StaticInteger, maybe_static, static

abstract type RadixSortAlgorithm <: Base.Sort.Algorithm end

(alg::RadixSortAlgorithm)(; kw...) = setproperties(alg; kw...)

include("utils.jl")
include("pool.jl")
include("msd.jl")

end # module


text:
using DelimitedFiles

function givessum(list, sum)
    # Returns the first pair of numbers in the list that summed gives sum,
    # returns nothing is there is no one that satisfies

    for (i, n) in enumerate(list)
        for m in list[i + 1:end]
            if n + m == sum
                return n, m
            end
        end
    end
    nothing
end

function firstinvalid(list, pmblsz)
    # Returns first invalid number, if none returns nothing.

    for (i, n) ∈ enumerate(list[pmblsz + 1:end])
        if isnothing(givessum(list[i:i + pmblsz - 1], n))
            # XMAS starts by transmitting a preamble of pmblsz numbers.
            # After that, each number you receive should be the sum of any two of
            # the pmblsz immediately previous numbers.

            return n
        end
    end
    nothing
end

# First Half

puzzleinput = readdlm("inputs/day09.txt", Int)
firstinvalid(puzzleinput, 25)

# Second Half

function findweakness(list, pmblsz)
    # Returns the encryption weakness in your XMAS-encrypted list of numbers,
    # if nonte returns nothing

    invalnum = firstinvalid(list, pmblsz)
    if !isnothing(invalnum)
        for i ∈ 1:length(list) - 1
            for j ∈ i + 1:length(list)
                contigset = list[i:j]
                if sum(contigset) == invalnum
                    # Find a contiguous set of at least two numbers in your list
                    # which sum to invalnum (the invalid number found before)

                    return minimum(contigset) + maximum(contigset)
                    # To find the encryption weakness, add together the smallest
                    # and largest number in this contiguous set

                end
            end
        end
    end
    nothing
end

puzzleinput = readdlm("inputs/day09.txt", Int)
findweakness(puzzleinput, 25)


text:
################################################################################
#
#  Map/NfRel.jl : Types for maps with domains of type NfRel
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

#why only simple relative?

mutable struct NfRelToNf <: Map{NfRel{nf_elem}, AnticNumberField, HeckeMap, NfRelToNf}
  header::MapHeader{NfRel{nf_elem}, AnticNumberField}

  function NfRelToNf(K::NfRel{nf_elem}, L::AnticNumberField, a::nf_elem, b::nf_elem, c::NfRelElem{nf_elem})
    # let K/k, k absolute number field
    # k -> L, gen(k) -> a
    # K -> L, gen(K) -> b
    # L -> K, gen(L) -> c

    k = K.base_ring
    Ly, y = PolynomialRing(L, cached = false)
    R = parent(k.pol)
    S = parent(L.pol)

    function image(x::NfRelElem{nf_elem})
      # x is an element of K
      f = data(x)
      # First evaluate the coefficients of f at a to get a polynomial over L
      # Then evaluate at b
      r = [ evaluate(R(coeff(f, i)), a) for i in 0:degree(f)]
      return evaluate(Ly(r), b)
    end

    function preimage(x::nf_elem)
      # x is an element of L
      f = S(x)
      return evaluate(f, c)
    end

    z = new()
    z.header = MapHeader(K, L, image, preimage)
    return z
  end
end

function show(io::IO, h::NfRelToNf)
  println(io, "Isomorphism between ", domain(h), "\nand ", codomain(h))
end

mutable struct NfRelRelToNfRel{T} <: Map{NfRel{NfRelElem{T}}, NfRel{T}, HeckeMap, NfRelRelToNfRel}
  header::MapHeader{NfRel{NfRelElem{T}}, NfRel{T}}

  function NfRelRelToNfRel(K::NfRel{NfRelElem{T}}, L::NfRel{T}, a::NfRelElem{T}, b::NfRelElem{T}, c::NfRelElem{NfRelElem{T}}) where T
    # let K/k, k absolute number field
    # k -> L, gen(k) -> a
    # K -> L, gen(K) -> b
    # L -> K, gen(L) -> c

    k = K.base_ring
    Ly, y = PolynomialRing(L, cached = false)
    R = parent(k.pol)
    S = parent(L.pol)

    function image(x::NfRelElem{NfRelElem{T}}) where T
      # x is an element of K
      f = data(x)
      # First evaluate the coefficients of f at a to get a polynomial over L
      # Then evaluate at b
      r = [ R(coeff(f, i))( a) for i in 0:degree(f)]
      return Ly(r)(b)
    end

    function preimage(x::NfRelElem{T}) where T
      # x is an element of L
      f = S(x)
      return f(c)
    end

    z = new{T}()
    z.header = MapHeader(K, L, image, preimage)
    return z
  end
end

function show(io::IO, h::NfRelRelToNfRel)
  println(io, "Isomorphism between ", domain(h), "\nand ", codomain(h))
end

mutable struct NfRelToNfRelMor{T, S} <: Map{NfRel{T}, NfRel{S}, HeckeMap, NfRelToNfRelMor}
  header::MapHeader{NfRel{T}, NfRel{S}}
  prim_img ::NfRelElem{S}
  coeff_aut::NfToNfMor

  function NfRelToNfRelMor{T, S}() where {T, S}
    z = new{T, S}()
    return z
  end

  function NfRelToNfRelMor{T, S}(K::NfRel{T}, L::NfRel{S}, a::NfRelElem{S}) where {T, S}
    function image(x::NfRelElem{S})
      # x is an element of K
      f = data(x)
      # First evaluate the coefficients of f at a to get a polynomial over L
      # Then evaluate at b
      return f(a)
    end

    z = new{T, S}()
    z.prim_img = a
    z.header = MapHeader(K, L, image)
    return z
  end
end

  #so far, only for single relative.
function NfRelToNfRelMor(K::NfRel{nf_elem}, L::NfRel{nf_elem}, A::NfToNfMor, a::NfRelElem{nf_elem})
  function image(x::NfRelElem{nf_elem})
    # x is an element of K
    f = data(x)
    g = zero(f)
    for i=0:degree(f)
      setcoeff!(g, i, A(coeff(f, i)))
    end
    # First evaluate the coefficients of f at a to get a polynomial over L
    # Then evaluate at b
    return g(a)
  end

  z = NfRelToNfRelMor{nf_elem, nf_elem}()
  z.prim_img = a
  z.coeff_aut = A
  z.header = MapHeader(K, L, image)
  return z
end

function NfRelToNfRelMor(K::NfRel{nf_elem}, L::NfRel{nf_elem}, a::NfRelElem{nf_elem})
  function image(x::NfRelElem{nf_elem})
    # x is an element of K
    f = data(x)
    g = zero(f)
    for i=0:degree(f)
      setcoeff!(g, i, coeff(f, i))
    end
    return g(a)
  end

  z = NfRelToNfRelMor{nf_elem, nf_elem}()
  z.prim_img = a
  z.header = MapHeader(K, L, image)
  return z
end

function show(io::IO, h::NfRelToNfRelMor)
  if domain(h) == codomain(h)
    println(io, "Automorphism of ", domain(h))
  else
    println(io, "Injection of ", domain(h), " into ", codomain(h))
  end
  println(io, "defined by ", gen(domain(h)), " -> ", h.prim_img)
  if isdefined(h, :coeff_aut)
    println(io, "using coefficient map: ", h.coeff_aut)
  end
end

mutable struct NfRelToFqMor{T} <: Map{NfRel{T}, FqFiniteField, HeckeMap, NfRelToFqMor}
  header::MapHeader{NfRel{T}, FqFiniteField}

  function NfRelToFqMor{T}() where {T}
    z = new{T}()
    z.header = MapHeader{NfRel{T}, FqFiniteField}()
    return z
  end
end


text:
abstract type Tax end

net(x, T::Tax)=x-tax(x, T)
gross(x, T::Tax)=x+tax(x, T)

# * Gouveia and Strauss (1994)
struct GouveiaStrauss <: Tax
    maxmargtax :: Float64
    coefexp    :: Float64
    coefcons   :: Float64
end

function GouveiaStrauss(;maxmargtax=0.275,
                        coefexp=5.,
                        coefcons=0.3)
    GouveiaStrauss(maxmargtax, coefexp, coefcons)
end

tax(y, gs::GouveiaStrauss)=_gouveiastrauss(y, gs.maxmargtax, gs.coefexp, gs.coefcons)

@doc doc"
    _gouveiastrauss(y, φ₀, φ₁, φ₂)

Inspired by Gouveia and Strauss (1994)

Parametrizes income taxes with a differentiable function characterized by three
parameters (the phis). If
```math
  \tau_y(y) = \phi_0 \left[ y - \left( y^{-\phi_1} + \phi_2  \right)^{-\frac{1}{\phi_1}} \right]
```

The $\phi_0$ coefficient is the limit $\lim \tau_y^{\prime}(y)$ when
$y\to\infty$. The other coefficients govern the degree of progressiveness of the
function. (For example, when $\phi_2 \to 0$, no one is taxed.)
"
_gouveiastrauss(y, φ₀, φ₁, φ₂)=φ₀*(y-(y^(-φ₁)+φ₂)^(-1/φ₁))

# * Simple linear schedule

struct AdValorem <: Tax
    t::Float64
end

tax(y, lav::AdValorem)=y*lav.t

# * Null tax

struct NullTax <: Tax end
tax(y, ::NullTax)=zero(y)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_HloReducePrecisionOptions_Location <: ProtoEnum
    OP_INPUTS::Int32
    OP_OUTPUTS::Int32
    UNFUSED_OP_OUTPUTS::Int32
    FUSION_INPUTS_BY_CONTENT::Int32
    FUSION_OUTPUTS_BY_CONTENT::Int32
    __enum_HloReducePrecisionOptions_Location() = new(0,1,2,3,4)
end #struct __enum_HloReducePrecisionOptions_Location
const HloReducePrecisionOptions_Location = __enum_HloReducePrecisionOptions_Location()

mutable struct HloReducePrecisionOptions <: ProtoType
    location::Int32
    exponent_bits::UInt32
    mantissa_bits::UInt32
    opcodes_to_suffix::Base.Vector{UInt32}
    opname_substrings_to_suffix::Base.Vector{AbstractString}
    HloReducePrecisionOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloReducePrecisionOptions
const __pack_HloReducePrecisionOptions = Symbol[:opcodes_to_suffix]
meta(t::Type{HloReducePrecisionOptions}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_HloReducePrecisionOptions, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct DebugOptions_XlaBackendExtraOptionsEntry <: ProtoType
    key::AbstractString
    value::AbstractString
    DebugOptions_XlaBackendExtraOptionsEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DebugOptions_XlaBackendExtraOptionsEntry (mapentry)

mutable struct DebugOptions <: ProtoType
    xla_generate_hlo_graph::AbstractString
    xla_hlo_graph_addresses::Bool
    xla_hlo_graph_path::AbstractString
    xla_hlo_dump_as_graphdef::Bool
    xla_log_hlo_text::AbstractString
    xla_generate_hlo_text_to::AbstractString
    xla_dump_optimized_hlo_proto_to::AbstractString
    xla_hlo_profile::Bool
    xla_dump_computations_to::AbstractString
    xla_dump_executions_to::AbstractString
    xla_disable_hlo_passes::Base.Vector{AbstractString}
    xla_backend_optimization_level::Int32
    xla_embed_ir_in_executable::Bool
    xla_dump_ir_to::AbstractString
    xla_eliminate_hlo_implicit_broadcast::Bool
    xla_cpu_multi_thread_eigen::Bool
    xla_gpu_cuda_data_dir::AbstractString
    xla_gpu_ftz::Bool
    xla_gpu_disable_multi_streaming::Bool
    xla_llvm_enable_alias_scope_metadata::Bool
    xla_llvm_enable_noalias_metadata::Bool
    xla_llvm_enable_invariant_load_metadata::Bool
    xla_llvm_disable_expensive_passes::Bool
    hlo_reduce_precision_options::Base.Vector{HloReducePrecisionOptions}
    xla_test_all_output_layouts::Bool
    xla_test_all_input_layouts::Bool
    xla_hlo_graph_sharding_color::Bool
    xla_hlo_tfgraph_device_scopes::Bool
    xla_gpu_use_cudnn_batchnorm::Bool
    xla_dump_unoptimized_hlo_proto_to::AbstractString
    xla_dump_per_pass_hlo_proto_to::AbstractString
    xla_cpu_use_mkl_dnn::Bool
    xla_gpu_max_kernel_unroll_factor::Int32
    xla_cpu_enable_fast_math::Bool
    xla_gpu_enable_fast_min_max::Bool
    xla_gpu_crash_on_verification_failures::Bool
    xla_force_host_platform_device_count::Int32
    xla_backend_extra_options::Base.Dict{AbstractString,AbstractString} # map entry
    DebugOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DebugOptions
const __fnum_DebugOptions = Int[1,2,4,5,6,7,8,9,10,11,30,31,33,34,35,60,61,62,63,70,71,72,73,80,90,91,92,93,94,95,96,97,98,99,100,101,102,500]
meta(t::Type{DebugOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_DebugOptions, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ExecutionOptions <: ProtoType
    shape_with_output_layout::ShapeProto
    seed::UInt64
    debug_options::DebugOptions
    device_handles::Base.Vector{DeviceHandle}
    ExecutionOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecutionOptions
const __fnum_ExecutionOptions = Int[2,3,4,5]
meta(t::Type{ExecutionOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ExecutionOptions, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GetDeviceHandlesRequest <: ProtoType
    device_count::Int64
    GetDeviceHandlesRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetDeviceHandlesRequest

mutable struct GetDeviceHandlesResponse <: ProtoType
    device_handles::Base.Vector{DeviceHandle}
    GetDeviceHandlesResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetDeviceHandlesResponse

mutable struct TransferToClientRequest <: ProtoType
    data::GlobalDataHandle
    shape_with_layout::ShapeProto
    TransferToClientRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferToClientRequest

mutable struct TransferToClientResponse <: ProtoType
    literal::LiteralProto
    TransferToClientResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferToClientResponse

mutable struct TransferToServerRequest <: ProtoType
    literal::LiteralProto
    device_handle::DeviceHandle
    TransferToServerRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferToServerRequest

mutable struct TransferToServerResponse <: ProtoType
    data::GlobalDataHandle
    TransferToServerResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferToServerResponse

mutable struct TransferToInfeedRequest <: ProtoType
    literal::LiteralProto
    replica_id::Int64
    device_handle::DeviceHandle
    TransferToInfeedRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferToInfeedRequest

mutable struct TransferToInfeedResponse <: ProtoType
    TransferToInfeedResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferToInfeedResponse

mutable struct TransferFromOutfeedRequest <: ProtoType
    shape_with_layout::ShapeProto
    replica_id::Int64
    device_handle::DeviceHandle
    TransferFromOutfeedRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferFromOutfeedRequest

mutable struct TransferFromOutfeedResponse <: ProtoType
    literal::LiteralProto
    TransferFromOutfeedResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TransferFromOutfeedResponse

mutable struct ResetDeviceRequest <: ProtoType
    device_handle::DeviceHandle
    ResetDeviceRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ResetDeviceRequest

mutable struct ResetDeviceResponse <: ProtoType
    ResetDeviceResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ResetDeviceResponse

mutable struct ComputationGraphStatsRequest <: ProtoType
    computation::HloModuleProto
    debug_options::DebugOptions
    ComputationGraphStatsRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ComputationGraphStatsRequest

mutable struct ComputationStatsResponse <: ProtoType
    stats::ComputationStats
    ComputationStatsResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ComputationStatsResponse

mutable struct CreateChannelHandleRequest <: ProtoType
    channel_type::Int32
    CreateChannelHandleRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CreateChannelHandleRequest

mutable struct CreateChannelHandleResponse <: ProtoType
    channel::ChannelHandle
    CreateChannelHandleResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CreateChannelHandleResponse

mutable struct UnregisterRequest <: ProtoType
    data::Base.Vector{GlobalDataHandle}
    UnregisterRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UnregisterRequest

mutable struct UnregisterResponse <: ProtoType
    UnregisterResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UnregisterResponse

mutable struct CompileRequest <: ProtoType
    computation::HloModuleProto
    execution_options::ExecutionOptions
    input_shape_with_layout::Base.Vector{ShapeProto}
    CompileRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CompileRequest

mutable struct CompileResponse <: ProtoType
    handle::ExecutionHandle
    CompileResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CompileResponse

mutable struct ExecuteRequest <: ProtoType
    handle::ExecutionHandle
    arguments::Base.Vector{GlobalDataHandle}
    ExecuteRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecuteRequest

mutable struct ExecuteGraphRequest <: ProtoType
    computation::HloModuleProto
    arguments::Base.Vector{GlobalDataHandle}
    execution_options::ExecutionOptions
    ExecuteGraphRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecuteGraphRequest

mutable struct ExecuteGraphParallelRequest <: ProtoType
    requests::Base.Vector{ExecuteGraphRequest}
    ExecuteGraphParallelRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecuteGraphParallelRequest

mutable struct ExecuteResponse <: ProtoType
    output::GlobalDataHandle
    profile::ExecutionProfile
    ExecuteResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecuteResponse

mutable struct ExecuteParallelResponse <: ProtoType
    responses::Base.Vector{ExecuteResponse}
    ExecuteParallelResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecuteParallelResponse

mutable struct WaitForExecutionRequest <: ProtoType
    execution::ExecutionHandle
    WaitForExecutionRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct WaitForExecutionRequest

mutable struct WaitForExecutionResponse <: ProtoType
    output::GlobalDataHandle
    profile::ExecutionProfile
    WaitForExecutionResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct WaitForExecutionResponse

mutable struct ComputeConstantGraphRequest <: ProtoType
    computation::HloModuleProto
    output_layout::Layout
    ComputeConstantGraphRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ComputeConstantGraphRequest

mutable struct ComputeConstantResponse <: ProtoType
    literal::LiteralProto
    ComputeConstantResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ComputeConstantResponse

mutable struct DeconstructTupleRequest <: ProtoType
    tuple_handle::GlobalDataHandle
    DeconstructTupleRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DeconstructTupleRequest
const __fnum_DeconstructTupleRequest = Int[2]
meta(t::Type{DeconstructTupleRequest}) = meta(t, ProtoBuf.DEF_REQ, __fnum_DeconstructTupleRequest, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct DeconstructTupleResponse <: ProtoType
    element_handles::Base.Vector{GlobalDataHandle}
    DeconstructTupleResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DeconstructTupleResponse

mutable struct LoadDataRequest <: ProtoType
    columnio_tablet_path::AbstractString
    columnio_field::AbstractString
    element_shape::ShapeProto
    offset::Int64
    limit::Int64
    zip::Bool
    LoadDataRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LoadDataRequest

mutable struct LoadDataResponse <: ProtoType
    data::GlobalDataHandle
    data_shape::ShapeProto
    available_rows::Int64
    rows_loaded::Int64
    nanoseconds::Int64
    LoadDataResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LoadDataResponse

mutable struct GetShapeRequest <: ProtoType
    data::GlobalDataHandle
    GetShapeRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetShapeRequest

mutable struct GetShapeResponse <: ProtoType
    shape::ShapeProto
    GetShapeResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetShapeResponse

mutable struct UnpackRequest <: ProtoType
    data::GlobalDataHandle
    UnpackRequest(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UnpackRequest

mutable struct UnpackResponse <: ProtoType
    tied_data::Base.Vector{GlobalDataHandle}
    UnpackResponse(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UnpackResponse

export HloReducePrecisionOptions_Location, HloReducePrecisionOptions, DebugOptions_XlaBackendExtraOptionsEntry, DebugOptions, ExecutionOptions, GetDeviceHandlesRequest, GetDeviceHandlesResponse, TransferToClientRequest, TransferToClientResponse, TransferToServerRequest, TransferToServerResponse, TransferToInfeedRequest, TransferToInfeedResponse, TransferFromOutfeedRequest, TransferFromOutfeedResponse, ResetDeviceRequest, ResetDeviceResponse, ComputationGraphStatsRequest, ComputationStatsResponse, CreateChannelHandleRequest, CreateChannelHandleResponse, UnregisterRequest, UnregisterResponse, CompileRequest, CompileResponse, ExecuteRequest, ExecuteGraphRequest, ExecuteGraphParallelRequest, ExecuteResponse, ExecuteParallelResponse, WaitForExecutionRequest, WaitForExecutionResponse, ComputeConstantGraphRequest, ComputeConstantResponse, DeconstructTupleRequest, DeconstructTupleResponse, LoadDataRequest, LoadDataResponse, GetShapeRequest, GetShapeResponse, UnpackRequest, UnpackResponse
# mapentries: "DebugOptions_XlaBackendExtraOptionsEntry" => ("AbstractString", "AbstractString")


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
Base.@kwdef struct ElecStorage <: AbstractStorage
    min_kw::Float64 = 0.0
    max_kw::Float64 = 1.0e4
    min_kwh::Float64 = 0.0
    max_kwh::Float64 = 1.0e6
    internal_efficiency_pct::Float64 = 0.975
    inverter_efficiency_pct::Float64 = 0.96
    rectifier_efficiency_pct::Float64 = 0.96
    soc_min_pct::Float64 = 0.2
    soc_init_pct::Float64 = 0.5
    can_grid_charge::Bool = true
    installed_cost_per_kw::Float64 = 840.0
    cost_per_kwh::Float64 = 420.0
    replace_cost_per_kw::Float64 = 410.0
    replace_cost_per_kwh::Float64 = 200.0
    inverter_replacement_year::Int = 10
    battery_replacement_year::Int = 10
    macrs_option_years::Int = 7
    macrs_bonus_pct::Float64 = 1.0
    macrs_itc_reduction::Float64 = 0.5
    total_itc_pct::Float64 = 0.0
    total_rebate_per_kw::Float64 = 0.0
end

# TODO change all DenseAxisArray's into Dicts
struct Storage <: AbstractStorage
    types::Array{Symbol,1}
    min_kw::DenseAxisArray{Float64,1}
    max_kw::DenseAxisArray{Float64,1}
    min_kwh::DenseAxisArray{Float64,1}
    max_kwh::DenseAxisArray{Float64,1}
    charge_efficiency::DenseAxisArray{Float64,1}
    discharge_efficiency::DenseAxisArray{Float64,1}
    soc_min_pct::DenseAxisArray{Float64,1}
    soc_init_pct::DenseAxisArray{Float64,1}
    installed_cost_per_kw::DenseAxisArray{Float64,1}
    cost_per_kwh::DenseAxisArray{Float64,1}
    can_grid_charge::Array{Symbol,1}
    grid_charge_efficiency::Dict{Symbol, Float64}
end

"""
    function Storage(d::Dict{Symbol,Dict}, f::Financial)

Construct Storage struct from Dict with keys for each storage type (eg. :elec) and values with
    input dicts for each storage type. Note that the REopt inputs are indexed on the storage type.
"""
function Storage(d::Dict, f::Financial)  # nested dict
    types = Symbol[]
    can_grid_charge = Symbol[]
    raw_vals = Dict(zip(fieldnames(Storage), [Float64[] for _ in range(1, stop=fieldcount(Storage))]))
    # raw_vals has keys = fieldnames(Storage) and values = arrays of values for each type in types

    for (storage_type, input_dict) in d

        push!(types, storage_type)
        struct_name = string(titlecase(string(storage_type)) * "Storage")  # eg. ElecStorage
        storage_instance = eval(Meta.parse(struct_name * "(;$input_dict...)"))

        if storage_instance.can_grid_charge
            push!(can_grid_charge, storage_type)
        end
        fill_storage_vals!(raw_vals, storage_instance, storage_type, f)
    end

    storage_args = Dict(
        :types => types,
        :can_grid_charge => can_grid_charge,
    )
    d2 = Dict()  # Julia won't let me use storage_args: "unable to check bounds for indices of type Symbol"
    for k in keys(raw_vals)
        d2[k] = DenseAxisArray(raw_vals[k], types)
    end

    grid_charge_efficiency = Dict(:elec => convert(Float64, d2[:charge_efficiency][:elec]))

    return Storage(
        storage_args[:types],
        d2[:min_kw],
        d2[:max_kw],
        d2[:min_kwh],
        d2[:max_kwh],
        d2[:charge_efficiency],
        d2[:discharge_efficiency],
        d2[:soc_min_pct],
        d2[:soc_init_pct],
        d2[:installed_cost_per_kw],
        d2[:cost_per_kwh],
        storage_args[:can_grid_charge],
        grid_charge_efficiency
    )
    # TODO expand for smart thermostat
end

function fill_storage_vals!(d::Dict{Symbol, Array{Float64,1}}, s::AbstractStorage, t::Symbol, f::Financial)
    push!(d[:min_kw], s.min_kw)
    push!(d[:max_kw], s.max_kw)
    push!(d[:min_kwh], s.min_kwh)
    push!(d[:max_kwh], s.max_kwh)
    push!(d[:soc_min_pct], s.soc_min_pct)
    push!(d[:soc_init_pct], s.soc_init_pct)

    push!(d[:charge_efficiency], s.rectifier_efficiency_pct * s.internal_efficiency_pct^0.5)
    push!(d[:discharge_efficiency], s.inverter_efficiency_pct * s.internal_efficiency_pct^0.5)

    push!(d[:installed_cost_per_kw], effective_cost(;
        itc_basis=s.installed_cost_per_kw,
        replacement_cost=s.replace_cost_per_kw,
        replacement_year=s.inverter_replacement_year,
        discount_rate=f.owner_discount_pct,
        tax_rate=f.owner_tax_pct,
        itc=s.total_itc_pct,
        macrs_schedule = s.macrs_option_years == 7 ? f.macrs_seven_year : f.macrs_five_year,
        macrs_bonus_pct=s.macrs_bonus_pct,
        macrs_itc_reduction = s.macrs_itc_reduction,
        rebate_per_kw = s.total_rebate_per_kw
    ))
    push!(d[:cost_per_kwh], effective_cost(;
        itc_basis=s.cost_per_kwh,
        replacement_cost=s.replace_cost_per_kwh,
        replacement_year=s.inverter_replacement_year,
        discount_rate=f.owner_discount_pct,
        tax_rate=f.owner_tax_pct,
        itc=s.total_itc_pct,
        macrs_schedule=s.macrs_option_years == 7 ? f.macrs_seven_year : f.macrs_five_year,
        macrs_bonus_pct=s.macrs_bonus_pct,
        macrs_itc_reduction = s.macrs_itc_reduction
    ))
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


text:
function gti2df(hdu::TableHDU)
    gti_start=read(hdu,"START")
    gti_stop = read(hdu,"STOP")
    gti=DataFrame(START=gti_start,STOP=gti_stop)
    return gti
end

function gti_unique_or(gti_list::DataFrame)
# join all GTIs
    merged_gti=vcat(gti_list)
    sort!(merged_gti,[:START,:STOP])
    gti_start = merged_gti[!,:START]
    gti_stop = merged_gti[!,:STOP]
# Determine unique GTIs
    unique_gti_start=Float64[]
    unique_gti_stop=Float64[]
    push!(unique_gti_start,gti_start[1])
    for i=1:length(gti_start)-1
        if (gti_stop[i] <= gti_start[i+1])
            push!(unique_gti_stop, gti_stop[i])
            push!(unique_gti_start, gti_start[i+1])
        end
    end
    push!(unique_gti_stop, gti_stop[end])
    unique_gti_df=DataFrame(START=unique_gti_start,STOP=unique_gti_stop)
    return unique_gti_df
end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
module LogGoodies

using Dates
using DataFrames
using Logging
using Logging: BelowMinLevel, Debug, Error, Info, Warn
using LoggingExtras
import Base.CoreLogging: shouldlog, min_enabled_level, handle_message
using Printf

export FlushedLogger
export @debug, @error, @info, @warn
export Debug, Error, Info, Warn, BelowMinLevel
export TeeLogger, global_logger, set_logger, append_default_flushed_logger
export dataframe_to_md

"""
FlushedLogger(stream=stderr, min_level=Info)
FlushedLogger(path, min_level=Info, append=false)

Always flushed logger.

Based on `Logging.SimpleLogger`.
"""
struct FlushedLogger <: AbstractLogger
    stream::IO
    min_level::LogLevel
    message_limits::Dict{Any,Int}
    dateformat::String
    showmodule::Bool
    showfile::Symbol
end

function FlushedLogger(stream::IO = stderr, level = BelowMinLevel; dateformat="HH:MM:SS", showmodule=false, showfile=:short)
    @assert showfile in [:short, :full, :none]
    FlushedLogger(stream, level, Dict{Any,Int}(), dateformat, showmodule, showfile)
end

function FlushedLogger(path::AbstractString, level = BelowMinLevel, append = false; kwargs...)
    FlushedLogger(open(path, append ? "a" : "w"), level; kwargs...)
end

shouldlog(logger::FlushedLogger, level, _module, group, id) = get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::FlushedLogger) = logger.min_level

catch_exceptions(logger::FlushedLogger) = false

function handle_message(logger::FlushedLogger, level, message, _module, group, id,
                        filepath, line; maxlog = nothing, kwargs...)

    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end
    buf = IOBuffer()
    iob = IOContext(buf, logger.stream)
    levelstr = uppercase(string(level))
    datestr = logger.dateformat == "" ? "" : " $(Dates.format(now(), logger.dateformat))"
    modulestr = logger.showmodule ? " $(something(_module, "nothing"))" : ""
    if logger.showfile == :none
        filestr = ""
    elseif logger.showfile == :short
        filestr = @sprintf " %s:%s"  basename(something(filepath, "nothing")) something(line, "nothing")
    else
        filestr = @sprintf " %s:%s"  something(filepath, "nothing") something(line, "nothing")
    end
    msglines = split(chomp(string(message)), '\n')

    println(iob, levelstr, datestr, modulestr, filestr, " ", msglines[1])
    for i in 2:length(msglines)
        println(iob, " ", msglines[i])
    end
    for (key, val) in kwargs
        println(iob, "   ", key, " = ", val)
    end
    write(logger.stream, take!(buf))
    flush(logger.stream)
    nothing
end

function append_default_flushed_logger(logfile = "log.txt", level = Info)
    old_logger = global_logger()
    tee_logger = TeeLogger(
        old_logger,
        FlushedLogger(logfile, level, showmodule=false, showfile=:short)
    )
    global_logger(tee_logger)

    @info "host name: $(gethostname())"
    @info "SLURM_JOB_ID: $(get(ENV, "SLURM_JOB_ID", "N/A"))"
    @info "SLURM_JOB_NAME: $(get(ENV, "SLURM_JOB_NAME", "N/A"))"
    @info "SLURM_SUBMIT_HOST: $(get(ENV, "SLURM_SUBMIT_HOST", "N/A"))"
    @info "PBS_JOBID: $(get(ENV, "PBS_JOBID", "N/A"))"
    @info "PBS_JOBDIR: $(get(ENV, "PBS_JOBDIR", "N/A"))"
end

# ============= FORMATTERS =============

# TODO: move to a dedicated package?
function dataframe_to_md(df::AbstractDataFrame)
    s = string(df)
    s = foldl(replace, [(c => "|" for c in ["│", "┼", "├", "┤"])..., "─" => "-"]; init = s) # replace UTF symbols
    s = [l * "\n" for (i, l) in enumerate(split(s, "\n")) if i != 3] |> join
    println(s)
end

end # module

max_stars_repo_name: dhonza/LogGoodies

max_issues_repo_name: dhonza/LogGoodies

max_forks_repo_name: dhonza/LogGoodies


text:
@testset "c\"union ...\"" begin
	@eval c``

	@eval c"""
	union OpaqueUnion;
	"""
	@eval c"""
	union OpaqueUnion;
	"""
	@test_throws ErrorException sizeof(c"union OpaqueUnion")
	@test sizeof(c"union OpaqueUnion *") == sizeof(Ptr{Cvoid})

	@eval c"""
	union EmptyUnion {
	} __attribute__((packed));
	"""
	@test sizeof(c"union EmptyUnion") == 0

	@eval c"""
	union CcharUnion {
		char c;
	} __attribute__((packed));
	"""
	@test sizeof(c"union CcharUnion") == sizeof(Cchar)
	@test :c in fieldnames(c"union CcharUnion")

	@eval c"""
	union CintUnion {
		int i;
	} __attribute__((packed));
	"""
	@test sizeof(c"union CintUnion") == sizeof(Cint)
	@test :i in fieldnames(c"union CintUnion")

	@eval c"""
	union CunionCintUnion {
		union CintUnion i;
	} __attribute__((packed));
	"""
	@test sizeof(c"union CunionCintUnion") == sizeof(c"union CintUnion")
	@test :i in fieldnames(c"union CunionCintUnion")

	@eval c"""
	union CintCcharUnion {
		union CintUnion i;
		union CcharUnion c;
	} __attribute__((packed));
	"""
	@test sizeof(c"union CintCcharUnion") == max(sizeof(c"union CintUnion"), sizeof(c"union CcharUnion"))
	@test :i in fieldnames(c"union CintCcharUnion")
	@test :c in fieldnames(c"union CintCcharUnion")

	@eval c"""
	union CcharCintPackedUnion {
		union CcharUnion c;
		union CintUnion i;
	} __attribute__((packed));
	"""
	@test sizeof(c"union CcharCintPackedUnion") == max(sizeof(c"union CintUnion"), sizeof(c"union CcharUnion"))
	@test :c in fieldnames(c"union CcharCintPackedUnion")
	@test :i in fieldnames(c"union CcharCintPackedUnion")

	@eval c"""
	union CcharCintAlignedUnion {
		union CcharUnion c;
		_Alignas(sizeof(union CintUnion)*2) union CintUnion i;
	} __attribute__((packed));
	"""
	@test sizeof(c"union CcharCintAlignedUnion") == max(max(sizeof(c"union CintUnion"), sizeof(Cint))*2, max(sizeof(c"union CcharUnion"), sizeof(Cint)))
	@test :c in fieldnames(c"union CcharCintAlignedUnion")
	@test :i in fieldnames(c"union CcharCintAlignedUnion")

	@eval c"""
	union CunionUnion {
		union {
			unsigned char c[4];
			unsigned short s[2];
			unsigned int i;
		} __attribute__((packed));
	} __attribute__((packed));
	"""
	@test :c in fieldnames(c"union CunionUnion")
	@test :s in fieldnames(c"union CunionUnion")
	@test :i in fieldnames(c"union CunionUnion")

	cuu = c"union CunionUnion"()
	@test cuu.i == 0
	@test cuu.c[1] == 0
	@test cuu.c[4] == 0
	@test cuu.s[1] == 0
	@test cuu.s[2] == 0

	cuu = c"union CunionUnion"(cuu, i = 0xff0000ff)
	@test cuu.i == 0xff0000ff
	@test cuu.c[1] == 0xff
	@test cuu.c[4] == 0xff
	@test cuu.s[1] != 0
	@test cuu.s[2] != 0

	@eval c"""
	union PtrToUnionUnion;
	"""
	@eval c"""
	union PtrToPtrToUnionUnion {
		union PtrToUnionUnion *p;
	} __attribute__((packed));
	"""
	@eval c"""
	union PtrToUnionUnion {
		union PtrToUnionUnion *p;
	} __attribute__((packed));
	"""
	@test sizeof(c"union PtrToPtrToUnionUnion") == sizeof(Ptr)
	@test :p in fieldnames(c"union PtrToPtrToUnionUnion")
	@test sizeof(c"union PtrToUnionUnion") == sizeof(Ptr)
	@test :p in fieldnames(c"union PtrToUnionUnion")

	@eval c"""
	union Cint32BitfieldUnion {
		int i:32;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cint32BitfieldUnion") == sizeof(Cint)
	@test :i in fieldnames(c"union Cint32BitfieldUnion")

	@eval c"""
	union Cuint32BitfieldUnion {
		unsigned int u:32;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cuint32BitfieldUnion") == sizeof(Cuint)
	@test :u in fieldnames(c"union Cuint32BitfieldUnion")

	@eval c"""
	union Cint2BitfieldUnion {
		int i:2;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cint2BitfieldUnion") == 1
	@test :i in fieldnames(c"union Cint2BitfieldUnion")

	@eval c"""
	union Cuint2BitfieldUnion {
		unsigned int u:2;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cuint2BitfieldUnion") == 1
	@test :u in fieldnames(c"union Cuint2BitfieldUnion")

	@eval c"""
	union Cuint32Cint32BitfieldsUnion {
		unsigned int u:32;
		int i:32;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cuint32Cint32BitfieldsUnion") == sizeof(Cuint)
	@test :u in fieldnames(c"union Cuint32Cint32BitfieldsUnion")
	@test :i in fieldnames(c"union Cuint32Cint32BitfieldsUnion")

	@eval c"""
	union Cuint16Cint16BitfieldsUnion {
		unsigned int u:16;
		int i:16;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cuint16Cint16BitfieldsUnion") == 2
	@test :u in fieldnames(c"union Cuint16Cint16BitfieldsUnion")
	@test :i in fieldnames(c"union Cuint16Cint16BitfieldsUnion")

	bfUnion = c"union Cuint16Cint16BitfieldsUnion"()
	@test bfUnion.u == 0
	@test bfUnion.i == 0

	bfUnion = c"union Cuint16Cint16BitfieldsUnion"(bfUnion, i = Cint(-1))
	@test bfUnion.u == 0xffff
	@test bfUnion.i == -1

	bfUnion = c"union Cuint16Cint16BitfieldsUnion"(bfUnion, u = Cuint(0x0001))
	@test bfUnion.u == 0x0001
	@test bfUnion.i == 1

	@eval c"""
	union Cint16Cint16BitfieldsUnion {
		int i1:16;
		int i2:16;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cint16Cint16BitfieldsUnion") == 2

	@eval c"""
	union Cint16CintBitfieldsUnion {
		int bf:16;
		int i;
	} __attribute__((packed));
	"""
	@test sizeof(c"union Cint16CintBitfieldsUnion") == sizeof(Cint)
end


text:
using Compat
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using LuxurySparse

n = 16
Id = IMatrix(1 << n)
V = randn(ComplexF64, 1 << n)
Dv = Diagonal(V)

Pm = pmrand(ComplexF64, 1 << n)
Sp = SparseMatrixCSC(Pm)

Ds = randn(ComplexF64, 1 << 10, 1 << 10)

mats = [Id, Sp, Pm, Dv]
for A in mats
    bg = BenchmarkGroup()
    for B in mats
        if !(A === B === Ds)
            println("============ $(typeof(A))-$(typeof(B)) =============")
            display(@benchmark $A * $B)
            println()
        end
    end
    display(run(bg, verbose = true))
end


text:
# generate LLVM wrappers

using LLVM_full_jll

using Clang.Generators

cd(@__DIR__)

@add_def off_t

function main()
    options = load_options(joinpath(@__DIR__, "wrap.toml"))

    includedir = LLVM_full_jll.llvm_config() do config
        readchomp(`$config --includedir`)
    end
    cppflags = LLVM_full_jll.llvm_config() do config
        split(readchomp(`$config --cppflags`))
    end

    args = get_default_args("x86_64-linux-gnu")
    push!(args, "-I$includedir")
    append!(args, cppflags)

    header_files = detect_headers(joinpath(includedir, "llvm-c"), args)

    ctx = create_context(header_files, args, options)

    build!(ctx, BUILDSTAGE_NO_PRINTING)

    build!(ctx, BUILDSTAGE_PRINTING_ONLY)
end

isinteractive() || main()