text:
"""
check_or_assume(value,units)

If `value` is not a Unitful quantity, assume it is in the given `units`.
"""
function check_or_assume(value,units)
    if isa(value,Unitful.Quantity)
        value
    else
        Unitful.Quantity(value,units)
    end
end

"""
    force_energy_units!(value)

If `value` is in temperature units, convert it to energy units.
"""
function force_energy_units!(value)
    if isa(value,Unitful.Temperature)
        value = kB*value
    end
    value
end


text:
#= -------------------------------------------------------------------------
# @file problema2_tiempo_60.jl
#
# @date 12/07/17 18:00:56
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
using DifferentialEquations
#=------------------------------------------------------------------------------
                              constants
------------------------------------------------------------------------------=#
const  A      =  2.0
const  R      =  2.390057361376673e-5
const  cₑ     =  4.184e6
const  V      =  0.5
const  θ_ref  =  333.0
const  θ₀     =  288.0
const  θ_a    =  293.0
G = (θ_ref - θ_a) / (R * A)
tspan = (0.0, 1000.0)
θ_dot(t, θ) = ((G * A) / (cₑ * V)) + ((θ_a / (R * cₑ * V))) - ((θ / (R * cₑ * V)))
prob = ODEProblem(θ_dot, θ₀, tspan)
condition(t,θ,integrator) = (θ - 333.0)
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
sol = solve(prob, callback=cb)

ext: jl


text:
module GoogleTrans
import HTTP
import JSON
export translate

function package_rpc(text::AbstractString, lang_src::String, lang_tgt::String)::String
    parameter = [[text, lang_src, lang_tgt, true], [1]]
    escaped_parameter = JSON.json(parameter)
    rpc = [[["MkEWBc", escaped_parameter, nothing, "generic"]]]
    escaped_rpc = JSON.json(rpc)
    freq = "f.req=" * HTTP.escapeuri(escaped_rpc) * "&"
    return freq
end

function translate(text::AbstractString, lang_tgt = "auto"::String, lang_src = "auto"::String)::String
    headers = Dict("Referer" => "http://translate.google.com", "User-Agent" => "Mozilla/5.0 (Windows NT 10.0; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/47.0.2526.106 Safari/537.36", "Content-type" => "application/x-www-form-urlencoded;charset=utf-8")
    freq = package_rpc(text, lang_src, lang_tgt)
    r = HTTP.request("POST", "https://translate.google.com/_/TranslateWebserverUi/data/batchexecute", headers, freq)
    JSON.parse(JSON.parse(split(String(r.body), '\n')[end])[1][3])[2][1][1][6][1][1]
end

translate(inputstrings::AbstractVector{S}, lang_tgt = "auto", lang_src = "auto") where {S<:AbstractString} = translate.(inputstrings, lang_tgt, lang_src)

end


text:
#=
  https://dtai.cs.kuleuven.be/problog/tutorial/basic/06_more_features.html
  """
  Suppose we are packing our bags to go on a trip. We have a set of items, each having a
  particular weight, and we pack each item with probability inversely proportional to its
  weight. We want to compute the probability that we will have excess baggage, i.e.,
  that the total weight of our baggage will exceed a given limit. We can model this with
  the following ProbLog program.

  """

  Note: This not the same approach as the ProbLog model, it's actually more like MiniZinc than Prolog.

  Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

   selected[1]    0.1588    0.3655     0.0037    0.0052   6215.6490    1.0011     1156.3998
   selected[2]    0.2455    0.4304     0.0043    0.0054   6377.7829    1.0001     1186.5643
   selected[3]    0.3322    0.4710     0.0047    0.0059   6058.6116    1.0000     1127.1836
   selected[4]    0.4870    0.4999     0.0050    0.0060   6763.9998    1.0000     1258.4186
    sumWeights    3.9054    3.3332     0.0333    0.0419   6509.3562    1.0002     1211.0430
  excess_limit    0.1237    0.3293     0.0033    0.0041   6781.6167    1.0001     1261.6961

  Distributions of variable sumWeights (num:0)
  0.00000 =>    2215  (0.221500)
  2.00000 =>    2112  (0.211200)
  3.00000 =>    1072  (0.107200)
  6.00000 =>     992  (0.099200)
  5.00000 =>     965  (0.096500)
  4.00000 =>     708  (0.070800)
  9.00000 =>     623  (0.062300)
  8.00000 =>     373  (0.037300)
  7.00000 =>     326  (0.032600)
  11.00000 =>     207  (0.020700)
  12.00000 =>     146  (0.014600)
  10.00000 =>     132  (0.013200)
  13.00000 =>      76  (0.007600)
  15.00000 =>      53  (0.005300)

  Distributions of variable (num:0)
  Any[0.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 0.0], [3.0, 3.0, 0.0], [4.0, 2.0, 0.0]]]	=>	2215 (0.2215)
  Any[2.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 0.0], [3.0, 3.0, 0.0], [4.0, 2.0, 1.0]]]	=>	2112 (0.2112)
  Any[3.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 0.0], [3.0, 3.0, 1.0], [4.0, 2.0, 0.0]]]	=>	1072 (0.1072)
  Any[5.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 0.0], [3.0, 3.0, 1.0], [4.0, 2.0, 1.0]]]	=>	965 (0.0965)
  Any[4.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 1.0], [3.0, 3.0, 0.0], [4.0, 2.0, 0.0]]]	=>	708 (0.0708)
  Any[6.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 1.0], [3.0, 3.0, 0.0], [4.0, 2.0, 1.0]]]	=>	606 (0.0606)
  Any[9.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 1.0], [3.0, 3.0, 1.0], [4.0, 2.0, 1.0]]]	=>	408 (0.0408)
  Any[6.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 0.0], [3.0, 3.0, 0.0], [4.0, 2.0, 0.0]]]	=>	386 (0.0386)
  Any[8.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 0.0], [3.0, 3.0, 0.0], [4.0, 2.0, 1.0]]]	=>	373 (0.0373)
  Any[7.0, [[1.0, 6.0, 0.0], [2.0, 4.0, 1.0], [3.0, 3.0, 1.0], [4.0, 2.0, 0.0]]]	=>	326 (0.0326)
  Any[9.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 0.0], [3.0, 3.0, 1.0], [4.0, 2.0, 0.0]]]	=>	215 (0.0215)
  Any[11.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 0.0], [3.0, 3.0, 1.0], [4.0, 2.0, 1.0]]]	=>	207 (0.0207)
  Any[12.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 1.0], [3.0, 3.0, 0.0], [4.0, 2.0, 1.0]]]	=>	146 (0.0146)
  Any[10.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 1.0], [3.0, 3.0, 0.0], [4.0, 2.0, 0.0]]]	=>	132 (0.0132)
  Any[13.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 1.0], [3.0, 3.0, 1.0], [4.0, 2.0, 0.0]]]	=>	76 (0.0076)
  Any[15.0, [[1.0, 6.0, 1.0], [2.0, 4.0, 1.0], [3.0, 3.0, 1.0], [4.0, 2.0, 1.0]]]	=>	53 (0.0053)

  Cf ~/blog/trip.blog
     ~/webppl/trip.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function trip()
    limit = 8
    # limit = 4.0
    # limit = 10.0
    # limit = 18.0

    skis = 1
    boots = 2
    helmet = 3
    gloves = 4
    items = [skis,boots,helmet,gloves]

    weights = Dict(skis=> 6,
                   boots=> 4,
                   helmet=> 3,
                   gloves=> 2
                  )
    n = 4
    selected = tzeros(n)
    for i in 1:n
        selected[i] ~ flip(1/weights[i])
    end

    sumWeights ~ Dirac(sum([weights[i]*selected[i] for i in 1:n]))

    # The probability the we excess the weight.
    function excess(v)
        return sumWeights > v
    end

    # Here we can enforce that we don't excess the weights
    # true ~ Dirac(sumWeights > 0 && excess(limit) == false)

    # [name,weight,selected]
    items_with_weights = [sumWeights, [[i,weights[i],selected[i]] for i in 1:n]]

    excess_limit ~ Dirac(excess(limit))

    return items_with_weights
end

model = trip()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:sumWeights)

chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chains_params)
show_var_dist_pct(genq)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using StatsBase

function random_subspaces(n_subspaces, min_dim = 1, max_dim = 2; min_ss_size = 2, max_ss_size=2)
    @assert min_dim < max_dim && min_dim >= 1
    @assert min_ss_size <= max_ss_size && min_ss_size >=2
    @assert max_ss_size <= max_dim - min_dim + 1

    max_possible_subspaces = sum([binomial(max_dim - min_dim + 1, k) for k in min_ss_size:max_ss_size])
    n_subspaces > max_possible_subspaces && error("Maximum number subspaces within the specification (min/max dim and ss_size) is $max_possible_subspaces,
        but method is called with n_subspaces = $n_subspaces.")
    subspaces = Set{Vector{Int}}()
    while length(subspaces) < n_subspaces
        subspace_size = rand(min_ss_size:max_ss_size)
        push!(subspaces, sort(sample(min_dim:max_dim, subspace_size; replace=false), rev=false))
    end
    return [s for s in subspaces]
end

function generate_randsub_from_data(n::Int; kwargs...)
    (data::Array{Float64, 2}) -> random_subspaces(n, 1, size(data, 1); kwargs...)
end

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
using BinDeps
using Base.Sys: isapple,isbsd,islinux,isunix,iswindows

@BinDeps.setup

libmjd_siggen = library_dependency("libmjd_siggen", aliases=[])

prefix = joinpath(BinDeps.depsdir(libmjd_siggen), "usr")
srcdir = joinpath(BinDeps.srcdir(libmjd_siggen), "mjd_siggen")
builddir = joinpath(BinDeps.builddir(libmjd_siggen), "mjd_siggen")

# BSD systems (other than macOS) use BSD Make rather than GNU Make by default
# We need GNU Make, and on such systems GNU make is invoked as `gmake`
make = isbsd() && !isapple() ? "gmake" : "make"

provides(SimpleBuild,
    (@build_steps begin
        if isdir(builddir)
            rm(builddir, recursive = true)
        end
        CreateDirectory(builddir)
        @build_steps begin
            ChangeDirectory(builddir)
            `cmake -DCMAKE_INSTALL_PREFIX="$prefix" $srcdir`
            `$make -j$(min(Sys.CPU_THREADS, 8)) install`
        end
    end), [libmjd_siggen], os = :Unix)

@BinDeps.install Dict(:libmjd_siggen => :libmjd_siggen)


text:
function AMFMdemod(ψ::Vector{ComplexF64}, t::Vector{Float64}; derivMethod="center11")::AMFMcompN
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  return AMFMcompN( abs.(ψ),
                    derivApprox(unwrap(angle.(ψ)),fs=1/(t[2]-t[1]), method=derivMethod),
                    real.(ψ),
                    imag.(ψ),
                    unwrap(angle.(ψ)),
                    t,
                    1/(t[2]-t[1]) )
end
function AMFMdemod(ψ::Vector{ComplexF64}, t::StepRangeLen; derivMethod="center11")::AMFMcompN
  return AMFMdemod(ψ::Vector{ComplexF64}, collect(t), derivMethod=derivMethod)
end

function derivApprox(x::Vector{Float64}; fs=1.0, method="center11")::Vector{Float64}
#References: http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences/
    if method == "forward"
        if length(x)<2; error("derivApprox:vector too short for selected method"); end
        x′ = append!([NaN],fs.*diff(x))
    end
    if method == "backward"
        if length(x)<2; error("derivApprox:vector too short for selected method"); end
        x′ = append!(fs.*diff(x),[NaN])
    end
    if method == "center3"
        if length(x)<3; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!([NaN], fs.*(x[3:end]-x[1:end-2])./2.0 ) ,[NaN])
    end
    if method == "center5"
        if length(x)<5; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,2), fs.*(x[1:end-4]-8.0.*x[2:end-3]+8.0.*x[4:end-1]-x[5:end])./12.0 ) ,NaN*ones(Float64,2))
    end
    if method == "center7"
        if length(x)<7; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,3), fs.*(-x[1:end-6]+ 9.0.*x[2:end-5]- 45.0.*x[3:end-4] +45.0.*x[5:end-2] -9.0.*x[6:end-1] + x[7:end]   )./60.0
 ) ,NaN*ones(Float64,3))
    end
    if method == "center9"
        if length(x)<9; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,4), fs.*(3.0.*x[1:end-8] - 32.0.*x[2:end-7] + 168.0.*x[3:end-6] -672.0.*x[4:end-5] + 672.0.*x[6:end-3] -168.0.*x[7:end-2] + 32.0.*x[8:end-1] -3.0.*x[9:end]  )./840.0
 ) ,NaN*ones(Float64,4))
    end
    if method == "center11"
        if length(x)<11; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,5), fs.*(-2.0.*x[1:end-10] +25.0.*x[2:end-9] -150.0.*x[3:end-8] +600.0.*x[4:end-7] -2100.0.*x[5:end-6] + 2100.0.*x[7:end-4]  -600.0.*x[8:end-3] +150.0.*x[9:end-2] -25.0.*x[10:end-1] +2.0.*x[11:end]  )./2520.0
 ) ,NaN*ones(Float64,5))
    end
    if method == "center13"
        if length(x)<13; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,6), fs.*(5.0.*x[1:end-12] - 72.0.*x[2:end-11] + 495.0.*x[3:end-10] -2200.0.*x[4:end-9]   +7425.0.*x[5:end-8]  -23760.0.*x[6:end-7]  +23760.0.*x[8:end-5]   + -7425.0.*x[9:end-4]   + 2200.0.*x[10:end-3] -495.0.*x[11:end-2] + 72.0.*x[12:end-1] -5.0.*x[13:end]  )./27720.0
 ) ,NaN*ones(Float64,6))
    end
    return x′
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
println("Welcome to rate3.jl, where rates for CIRIS are plotted")
using DataFrames
using Gadfly

path = "/scratch/cns7ae/$DIRNAME/"
outfile = "rates-plot.pdf"
basefile = "_reactions.wsv"
numspecies = 16
xmax = 1e15
xmin = 1e13

plottheme = Gadfly.Theme(
    highlight_width=0pt,
    default_point_size=2pt,
    line_width=2pt,
    minor_label_font_size=13pt,
    minor_label_color=colorant"black",
    major_label_font_size=15pt,
    major_label_color=colorant"black",
    key_label_color=colorant"black",
    key_title_color=colorant"black",
    key_title_font_size=8pt,
    key_label_font_size=6pt,
    grid_color=colorant"black")

p = Array(Gadfly.Plot,numspecies)
d = Array(Gadfly.Plot,numspecies)

for n=1:numspecies

    println("Plotting data for species $n")

    df = readtable(
        path*string(n)*basefile,
        header = true,
        names = [:Percent,:Time, :Fluence, :Label, :Type],
        eltypes = [Float64,Float64,Float64, String, String]
    )

    proddf = DataFrame()
    proddf = df[df[:Type].=="Production",:]

    destdf = DataFrame
    destdf = df[df[:Type].=="Destruction",:]

    p[n] = plot(
        proddf,
        x="Fluence",
        y="Percent",
        color = "Label",
        Geom.line,
        Scale.x_log10,
        Scale.y_log10,
        Guide.xlabel("Fluence"),
        Guide.ylabel("Fractional Importance")
    )

    push!(p[n],plottheme)

    d[n] = plot(
        destdf,
        x="Fluence",
        y="Percent",
        color = "Label",
        Geom.line,
        Scale.x_log10,
        Scale.y_log10,
        Guide.xlabel("Fluence"),
        Guide.ylabel("Fractional Importance")
    )

    push!(d[n],plottheme)
end

println("Ending rate3.jl!")


text:
module PP2PM
export load_pm_from_json, get_model, get_solver

import JSON
using PowerModels

import Cbc
import Ipopt
import Juniper
import JuMP

try
    import Gurobi
catch e
    if isa(e, LoadError)
        println("Cannot import Gurobi. That's fine if you do not plan to use it")
    end
end

try
    import KNITRO
catch e
    if isa(e, LoadError)
        println("Cannot import KNITRO. That's fine if you do not plan to use it")
    end
end

try
    import SCIP
catch e
    if isa(e, LoadError)
        println("Cannot import SCIP. That's fine if you do not plan to use it")
    end
end

function get_model(model_type)
    """
    gets the model function
    model_type (str) - examples: "ACPPowerModel", "DCPPowerModel", "SOCWRPowerModel"...
    see: https://lanl-ansi.github.io/PowerModels.jl/stable/formulation-details/
    """

    s = Symbol(model_type)
    return getfield(Main, s)
end

function get_solver(optimizer::String, nl::String="ipopt", mip::String="cbc",
    log_level::Int=0, time_limit::Float64=Inf, nl_time_limit::Float64=Inf,
    mip_time_limit::Float64=Inf)

    if optimizer == "gurobi"
            solver = JuMP.with_optimizer(Gurobi.Optimizer, TimeLimit=time_limit)
    end

    if optimizer == "ipopt"
                solver = JuMP.with_optimizer(Ipopt.Optimizer, print_level=log_level, max_cpu_time=time_limit)
    end

    if optimizer == "juniper" && nl == "ipopt" && mip == "cbc"
        mip_solver = JuMP.with_optimizer(Cbc.Optimizer, logLevel=log_level, seconds=mip_time_limit)
        nl_solver = JuMP.with_optimizer(Ipopt.Optimizer, print_level=log_level, max_cpu_time=nl_time_limit)
        solver = JuMP.with_optimizer(Juniper.Optimizer,
                     nl_solver = nl_solver,
                     mip_solver = mip_solver,
                     log_levels = [],
                     time_limit = time_limit)
    end

    if optimizer == "juniper" && nl == "gurobi" && mip == "cbc"
        mip_solver = JuMP.with_optimizer(Cbc.Optimizer, logLevel=log_level, seconds=mip_time_limit)
        nl_solver = JuMP.with_optimizer(Gurobi.Optimizer, TimeLimit=nl_time_limit)
        solver = JuMP.with_optimizer(Juniper.Optimizer,
                     nl_solver = nl_solver,
                     mip_solver = mip_solver,
                     log_levels = [],
                     time_limit = time_limit)
    end

    if optimizer == "knitro"
        solver = JuMP.with_optimizer(KNITRO.Optimizer)
    end

    if optimizer == "cbc"
        solver = JuMP.with_optimizer(Cbc.Optimizer, seconds=time_limit)
    end

    if optimizer == "scip"
        solver = JuMP.with_optimizer(SCIP.Optimizer)
    end

    return solver

end

function load_pm_from_json(json_path)
    pm = Dict()
    open(json_path, "r") do f
        pm = JSON.parse(f)  # parse and transform data
    end

    for (idx, gen) in pm["gen"]
        if gen["model"] == 1
            pm["gen"][idx]["cost"] = convert(Array{Float64,1}, gen["cost"])
        end
    end
    if pm["correct_pm_network_data"]
        correct_network_data!(pm)
    end
    return pm
end

end

max_stars_repo_name: andreabragantini/pandapower

max_issues_repo_name: andreabragantini/pandapower

max_forks_repo_name: venidera/pandapower


text:
""" Generate and parse computer programs representing morphisms.
"""
module Programs

using Reexport

include("GenerateJuliaPrograms.jl")
include("ParseJuliaPrograms.jl")
include("RelationalPrograms.jl")

@reexport using .GenerateJuliaPrograms
@reexport using .ParseJuliaPrograms
@reexport using .RelationalPrograms

end


text:
# RecipesPipeline API

## Warnings

function RecipesPipeline.warn_on_recipe_aliases!(
    plt::Plot,
    plotattributes::AKW,
    recipe_type,
    signature_string
)
    for k in keys(plotattributes)
        if !is_default_attribute(k)
            dk = get(_keyAliases, k, k)
            if k !== dk
                @warn "Attribute alias `$k` detected in the $recipe_type recipe defined for the signature $signature_string. To ensure expected behavior it is recommended to use the default attribute `$dk`."
            end
            plotattributes[dk] = RecipesPipeline.pop_kw!(plotattributes, k)
        end
    end
end

## Grouping

RecipesPipeline.splittable_attribute(plt::Plot, key, val::SeriesAnnotations, len) =
    RecipesPipeline.splittable_attribute(plt, key, val.strs, len)

function RecipesPipeline.split_attribute(plt::Plot, key, val::SeriesAnnotations, indices)
    split_strs = RecipesPipeline.split_attribute(plt, key, val.strs, indices)
    return SeriesAnnotations(split_strs, val.font, val.baseshape, val.scalefactor)
end

## Preprocessing attributes
function RecipesPipeline.preprocess_axis_args!(plt::Plot, plotattributes, letter)
    # Fix letter for seriestypes that are x only but data gets passed as y
    if treats_y_as_x(get(plotattributes, :seriestype, :path))
        if get(plotattributes, :orientation, :vertical) == :vertical
            letter = :x
        end
    end

    plotattributes[:letter] = letter
    RecipesPipeline.preprocess_axis_args!(plt, plotattributes)
end

RecipesPipeline.preprocess_attributes!(plt::Plot, plotattributes) =
    RecipesPipeline.preprocess_attributes!(plotattributes) # in src/args.jl

RecipesPipeline.is_axis_attribute(plt::Plot, attr) = is_axis_attr_noletter(attr) # in src/args.jl

RecipesPipeline.is_subplot_attribute(plt::Plot, attr) = is_subplot_attr(attr) # in src/args.jl

## User recipes

function RecipesPipeline.process_userrecipe!(plt::Plot, kw_list, kw)
    _preprocess_userrecipe(kw)
    warn_on_unsupported_scales(plt.backend, kw)
    # add the plot index
    plt.n += 1
    kw[:series_plotindex] = plt.n

    push!(kw_list, kw)
    _add_errorbar_kw(kw_list, kw)
    _add_smooth_kw(kw_list, kw)
    return
end

function _preprocess_userrecipe(kw::AKW)
    _add_markershape(kw)

    # map marker_z if it's a Function
    if isa(get(kw, :marker_z, nothing), Function)
        # TODO: should this take y and/or z as arguments?
        kw[:marker_z] = isa(kw[:z], Nothing) ? map(kw[:marker_z], kw[:x], kw[:y]) :
            map(kw[:marker_z], kw[:x], kw[:y], kw[:z])
    end

    # map line_z if it's a Function
    if isa(get(kw, :line_z, nothing), Function)
        kw[:line_z] = isa(kw[:z], Nothing) ? map(kw[:line_z], kw[:x], kw[:y]) :
            map(kw[:line_z], kw[:x], kw[:y], kw[:z])
    end

    # convert a ribbon into a fillrange
    if get(kw, :ribbon, nothing) !== nothing
        make_fillrange_from_ribbon(kw)
    end
    return
end

function _add_errorbar_kw(kw_list::Vector{KW}, kw::AKW)
    # handle error bars by creating new recipedata data... these will have
    # the same recipedata index as the recipedata they are copied from
    for esym in (:xerror, :yerror, :zerror)
        if get(kw, esym, nothing) !== nothing
            # we make a copy of the KW and apply an errorbar recipe
            errkw = copy(kw)
            errkw[:seriestype] = esym
            errkw[:label] = ""
            errkw[:primary] = false
            push!(kw_list, errkw)
        end
    end
end

function _add_smooth_kw(kw_list::Vector{KW}, kw::AKW)
    # handle smoothing by adding a new series
    if get(kw, :smooth, false)
        x, y = kw[:x], kw[:y]
        β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)
        sx = [ignorenan_minimum(x), ignorenan_maximum(x)]
        sy = β .* sx .+ α
        push!(
            kw_list,
            merge(
                copy(kw),
                KW(
                    :seriestype => :path,
                    :x => sx,
                    :y => sy,
                    :fillrange => nothing,
                    :label => "",
                    :primary => false,
                ),
            ),
        )
    end
end

RecipesPipeline.get_axis_limits(plt::Plot, f, letter) = axis_limits(plt[1], letter)

## Plot recipes

RecipesPipeline.type_alias(plt::Plot) = get(_typeAliases, st, st)

## Plot setup

function RecipesPipeline.plot_setup!(plt::Plot, plotattributes, kw_list)
    _plot_setup(plt, plotattributes, kw_list)
    _subplot_setup(plt, plotattributes, kw_list)
    return nothing
end

function RecipesPipeline.process_sliced_series_attributes!(plt::Plots.Plot, kw_list)
    # swap errors
    err_inds = findall(kw -> get(kw, :seriestype, :path) in (:xerror, :yerror, :zerror), kw_list)
    for ind in err_inds
        if get(kw_list[ind-1],:seriestype,:path) == :scatter
            tmp = copy(kw_list[ind])
            kw_list[ind] = copy(kw_list[ind-1])
            kw_list[ind-1] = tmp
        end
    end
    return nothing
end

# TODO: Should some of this logic be moved to RecipesPipeline?
function _plot_setup(plt::Plot, plotattributes::AKW, kw_list::Vector{KW})
    # merge in anything meant for the Plot
    for kw in kw_list, (k, v) in kw
        haskey(_plot_defaults, k) && (plotattributes[k] = pop!(kw, k))
    end

    # TODO: init subplots here
    _update_plot_args(plt, plotattributes)
    if !plt.init
        plt.o = Base.invokelatest(_create_backend_figure, plt)

        # create the layout and subplots from the inputs
        plt.layout, plt.subplots, plt.spmap = build_layout(plt.attr)
        for (idx, sp) in enumerate(plt.subplots)
            sp.plt = plt
            sp.attr[:subplot_index] = idx
        end

        plt.init = true
    end

    # handle inset subplots
    insets = plt[:inset_subplots]
    if insets !== nothing
        if !(typeof(insets) <: AVec)
            insets = [insets]
        end
        for inset in insets
            parent, bb = is_2tuple(inset) ? inset : (nothing, inset)
            P = typeof(parent)
            if P <: Integer
                parent = plt.subplots[parent]
            elseif P == Symbol
                parent = plt.spmap[parent]
            else
                parent = plt.layout
            end
            sp = Subplot(backend(), parent = parent)
            sp.plt = plt
            push!(plt.subplots, sp)
            push!(plt.inset_subplots, sp)
            sp.attr[:relative_bbox] = bb
            sp.attr[:subplot_index] = length(plt.subplots)
        end
    end
    plt[:inset_subplots] = nothing
end

function _subplot_setup(plt::Plot, plotattributes::AKW, kw_list::Vector{KW})
    # we'll keep a map of subplot to an attribute override dict.
    # Subplot/Axis attributes set by a user/series recipe apply only to the
    # Subplot object which they belong to.
    # TODO: allow matrices to still apply to all subplots
    sp_attrs = Dict{Subplot, Any}()
    for kw in kw_list
        # get the Subplot object to which the series belongs.
        sps = get(kw, :subplot, :auto)
        sp = get_subplot(
            plt,
            _cycle(
                sps == :auto ? plt.subplots : plt.subplots[sps],
                series_idx(kw_list, kw),
            ),
        )
        kw[:subplot] = sp

        # extract subplot/axis attributes from kw and add to sp_attr
        attr = KW()
        for (k, v) in collect(kw)
            if is_subplot_attr(k) || is_axis_attr(k)
                attr[k] = pop!(kw, k)
            end
            if is_axis_attr_noletter(k)
                v = pop!(kw, k)
                for letter in (:x, :y, :z)
                    attr[Symbol(letter, k)] = v
                end
            end
            for k in (:scale,), letter in (:x, :y, :z)
                # Series recipes may need access to this information
                lk = Symbol(letter, k)
                if haskey(attr, lk)
                    kw[lk] = attr[lk]
                end
            end
        end
        sp_attrs[sp] = attr
    end

    # override subplot/axis args.  `sp_attrs` take precendence
    for (idx, sp) in enumerate(plt.subplots)
        attr = if !haskey(plotattributes, :subplot) || plotattributes[:subplot] == idx
            merge(plotattributes, get(sp_attrs, sp, KW()))
        else
            get(sp_attrs, sp, KW())
        end
        _update_subplot_args(plt, sp, attr, idx, false)
    end

    # do we need to link any axes together?
    link_axes!(plt.layout, plt[:link])
end

function series_idx(kw_list::AVec{KW}, kw::AKW)
    Int(kw[:series_plotindex]) - Int(kw_list[1][:series_plotindex]) + 1
end

## Series recipes

function RecipesPipeline.slice_series_attributes!(plt::Plot, kw_list, kw)
    sp::Subplot = kw[:subplot]
    # in series attributes given as vector with one element per series,
    # select the value for current series
    _slice_series_args!(kw, plt, sp, series_idx(kw_list, kw))
    return nothing
end

RecipesPipeline.series_defaults(plt::Plot) = _series_defaults # in args.jl

RecipesPipeline.is_seriestype_supported(plt::Plot, st) = is_seriestype_supported(st)

function RecipesPipeline.add_series!(plt::Plot, plotattributes)
    sp = _prepare_subplot(plt, plotattributes)
    _expand_subplot_extrema(sp, plotattributes, plotattributes[:seriestype])
    _update_series_attributes!(plotattributes, plt, sp)
    _add_the_series(plt, sp, plotattributes)
end

# getting ready to add the series... last update to subplot from anything
# that might have been added during series recipes
function _prepare_subplot(plt::Plot{T}, plotattributes::AKW) where {T}
    st::Symbol = plotattributes[:seriestype]
    sp::Subplot{T} = plotattributes[:subplot]
    sp_idx = get_subplot_index(plt, sp)
    _update_subplot_args(plt, sp, plotattributes, sp_idx, true)

    st = _override_seriestype_check(plotattributes, st)

    # change to a 3d projection for this subplot?
    if RecipesPipeline.needs_3d_axes(st)
        sp.attr[:projection] = "3d"
    end

    # initialize now that we know the first series type
    if !haskey(sp.attr, :init)
        _initialize_subplot(plt, sp)
        sp.attr[:init] = true
    end
    sp
end

function _override_seriestype_check(plotattributes::AKW, st::Symbol)
    # do we want to override the series type?
    if !RecipesPipeline.is3d(st) && !(st in (:contour, :contour3d))
        z = plotattributes[:z]
        if !isa(z, Nothing) &&
           (size(plotattributes[:x]) == size(plotattributes[:y]) == size(z))
            st = (st == :scatter ? :scatter3d : :path3d)
            plotattributes[:seriestype] = st
        end
    end
    st
end

function _expand_subplot_extrema(sp::Subplot, plotattributes::AKW, st::Symbol)
    # adjust extrema and discrete info
    if st == :image
        xmin, xmax = ignorenan_extrema(plotattributes[:x])
        ymin, ymax = ignorenan_extrema(plotattributes[:y])
        expand_extrema!(sp[:xaxis], (xmin, xmax))
        expand_extrema!(sp[:yaxis], (ymin, ymax))
    elseif !(st in (:histogram, :bins2d, :histogram2d))
        expand_extrema!(sp, plotattributes)
    end
    # expand for zerolines (axes through origin)
    if sp[:framestyle] in (:origin, :zerolines)
        expand_extrema!(sp[:xaxis], 0.0)
        expand_extrema!(sp[:yaxis], 0.0)
    end
end

function _add_the_series(plt, sp, plotattributes)
    extra_kwargs = warn_on_unsupported_args(plt.backend, plotattributes)
    if (kw = plt[:extra_kwargs]) isa AbstractDict
        plt[:extra_plot_kwargs] = get(kw,:plot,KW())
        sp[:extra_kwargs] = get(kw,:subplot, KW())
        plotattributes[:extra_kwargs] = get(kw,:series,KW())
    elseif plt[:extra_kwargs] == :plot
        plt[:extra_plot_kwargs] = extra_kwargs
    elseif plt[:extra_kwargs] == :subplot
        sp[:extra_kwargs] = extra_kwargs
    elseif plt[:extra_kwargs] == :series
        plotattributes[:extra_kwargs] = extra_kwargs
    else
        ArgumentError("Unsupported type for extra keyword arguments")
    end
    warn_on_unsupported(plt.backend, plotattributes)
    series = Series(plotattributes)
    push!(plt.series_list, series)
    push!(sp.series_list, series)
    _series_added(plt, series)
end


text:
module SaveREPL

export saveREPL
export printREPL

immutable REPLEntry
    time::String
    mode::String
    command::String
end

"""
    saveREPL(filename::String, n::Int=10)

Save `n` last executed Julia commands to script file `filename`.

Order: oldest command first. Ignores help and shell mode commands as well as
`printREPL` and `saveREPL` calls.
"""
function saveREPL(filename::String, n::Int=10)
    open(filename, "w") do f
        # write(f, script(julia_history, lines))
        entries = history(n+1)
        write(f, join((x->x.command).(entries[1:end-1]), "\n"))
    end
    nothing
end

"""
    printREPL(n::Int=10)

Prints the `n` last executed Julia commands.

Order: latest command first. Ignores help and shell mode commands as well as
`printREPL` and `saveREPL` calls.
"""
function printREPL(n::Int=10)
    entries = history(n+1)
    for e in reverse(entries[1:end-1])
        println(e.command)
    end
    nothing
end

"""
    history(n::Int)

Load last `n` executed Julia commands from `.julia_history` file.
"""
function history(n::Int)
    h = reverse(readlines(Base.REPL.find_hist_file()))
    entries = REPLEntry[]
    i = 1
    N = length(h)
    c = 0
    while i<=N && c<n
        line = h[i]
        !contains(line, "saveREPL(") || (begin i+=1; continue; end)
        !contains(line, "printREPL(") || (begin i+=1; continue; end)

        cmdlines = String[]
        while startswith(line, "\t")
            push!(cmdlines, replace(line, "\t", "", 1))
            i+=1
            line = h[i]
        end
        command = join(reverse(cmdlines), "\n")

        contains(line, "# mode:") || warn("wrong order: expected mode")
        mode = replace(chomp(line), "# mode: ", "")

        i+=1
        line = h[i]
        contains(line, "# time: ") || warn("wrong order: expected time")
        time = replace(chomp(line), "# time: ", "")
        i+=1

        contains(mode, "julia") || continue
        push!(entries, REPLEntry(time, mode, command))
        c+=1
    end
    return reverse(entries)
end

end # module


text:
# computes the pdf of the 2SLS estimator under multivariate t errors
using QuadraticFormsMGHyp
using Random
using ToeplitzMatrices
using LinearAlgebra

Random.seed!(0)
order = 2
n = 25
beta = 0.
k = 1
Z = randn(n, k)
In = diagm(0=>ones(n))
Z = In[:, 1:k] * (Z'*Z)^(.5) # distribution depends on Z only through Z'Z if R=I
alpha = 0.
if abs(alpha) >= 1.
     bb = 0.
else
     bb = 1 / sqrt(1-alpha^2)
end
R = Toeplitz(alpha .^(0:n-1), vcat(1., zeros(n-1))) * I
R[:, 1] = bb * R[:, 1]
R = R ./ sqrt(tr(R * R') / n)
mu = .5
Pi = sqrt(mu / k) * (Z' * Z)^(-.5) * ones(k, 1)
s2u = 1.
s2v = 1.
rho = 1.
suv = sqrt(s2u * s2v) * rho
nuvec = [3, 9]

Zp = Z * Pi
Pz = Z * inv(Z' * Z) * Z'
xvec = beta-3:.01:beta+3
#E = eigen(.5*Pz+.5*Pz')
#[pp,ll]=eig(.5*Pz+.5*Pz')
nx = length(xvec)
cdf = zeros(length(nuvec), nx)
pm = similar(cdf)
spacdf = similar(cdf)
spapm =  similar(cdf)

for nuloop=1:length(nuvec)
    nu = nuvec[nuloop]
    for loop = 1 : nx
        global cdf, spacdf, pm, spapm
        Si = ((nu-2) / nu * [s2u suv; suv s2v])^.5 # rescale to make Si^2 the covariance matrix
        S = kron(Si, R*R')
        local b = xvec[loop] - beta
        local a0 = (-b * (Zp' * Zp))[1]
        local a = [Zp; -2*b*Zp]
        local A = .5*[Pz*0 Pz; Pz -2*b*Pz]
        SAS = S * A * S
        SAS = .5 * (SAS + SAS')
        E = eigen(SAS)
        P = E.vectors
        omega = E.values
        local d = a' * S * P
        ccdf, p = qfmgh(0., a0, a[:], A, S, zeros(2n), zeros(2n), -nu/2, nu, 0., do_spa=false)
        cdf[nuloop, loop] = 1.0 - ccdf
        pm[nuloop, loop] = p
        ccdf, p = qfmgh(0., a0, a[:], A, S, zeros(2n), zeros(2n), -nu/2, nu, 0., do_spa=true, order=order)
        spacdf[nuloop, loop] = 1.0 - ccdf
        spapm[nuloop, loop] = p
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation

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 Intel Corporation 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 OWNER 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.
=#

using Latte

net = Net(64)
data, label  = HDF5DataLayer(net, "data/train.txt", "data/val.txt")

conv1_1 = ConvolutionLayer( :conv1_1, net, data,    64,  3, 1, 1)
relu1_1 = ReLULayer(        :relu1_1, net, conv1_1)
conv1_2 = ConvolutionLayer( :conv1_2, net, relu1_1, 64,  3, 1, 1)
relu1_2 = ReLULayer(        :relu1_2, net, conv1_2)
pool1   = MaxPoolingLayer(  :pool1,   net, relu1_2, 2,   2, 0)

conv2_1 = ConvolutionLayer( :conv2_1, net, pool1,   128, 3, 1, 1)
relu2_1 = ReLULayer(        :relu2_1, net, conv2_1)
conv2_2 = ConvolutionLayer( :conv2_2, net, relu2_1, 128, 3, 1, 1)
relu2_2 = ReLULayer(        :relu2_2, net, conv2_2)
pool2   = MaxPoolingLayer(  :pool2,   net, relu2_2, 2,   2, 0)

conv3_1 = ConvolutionLayer( :conv3_1, net, pool2,   256, 3, 1, 1)
relu3_1 = ReLULayer(        :relu3_1, net, conv3_1)
conv3_2 = ConvolutionLayer( :conv3_2, net, relu3_1, 256, 3, 1, 1)
relu3_2 = ReLULayer(        :relu3_2, net, conv3_2)
conv3_3 = ConvolutionLayer( :conv3_3, net, relu3_2, 256, 3, 1, 1)
relu3_3 = ReLULayer(        :relu3_3, net, conv3_3)
conv3_4 = ConvolutionLayer( :conv3_4, net, relu3_3, 256, 3, 1, 1)
relu3_4 = ReLULayer(        :relu3_4, net, conv3_4)
pool3   = MaxPoolingLayer(  :pool3,   net, relu3_4, 2,   2, 0)

fc4     = InnerProductLayer(:fc4,     net, pool3,   1024)
relu4   = ReLULayer(        :relu4,   net, fc4)
# drop4   = DropoutLayer(     :drop4,   net, relu4, .5f0)
# fc5     = InnerProductLayer(:fc5,     net, drop4,   1024)
fc5     = InnerProductLayer(:fc5,     net, relu4,   1024)
relu5   = ReLULayer(        :relu5,   net, fc5)
# drop5   = DropoutLayer(     :drop5,   net, relu5, .5f0)
# fc6     = InnerProductLayer(:fc6,     net, drop5,   10)
fc6     = InnerProductLayer(:fc6,     net, relu5,   200)

loss     = SoftmaxLossLayer(:loss, net, fc6, label)
accuracy = AccuracyLayer(:accuracy, net, fc6, label)

params = SolverParameters(
    LRPolicy.Step(.01, .1, 100000),
    MomPolicy.Fixed(0.9),
    100000,
    .0005,
    1000)
sgd = SGD(params)
solve(sgd, net)


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

struct MetricsV1beta1Api <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_listMetricsV1beta1NodeMetrics(_api::MetricsV1beta1Api; allowWatchBookmarks=nothing, __continue__=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, resourceVersion=nothing, timeoutSeconds=nothing, watch=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApiMetricsV1beta1NodeMetricsList, "/apis/metrics.k8s.io/v1beta1/nodes", ["BearerToken"])
    Swagger.set_param(_ctx.query, "allowWatchBookmarks", allowWatchBookmarks)  # type Bool
    Swagger.set_param(_ctx.query, "continue", __continue__)  # type String
    Swagger.set_param(_ctx.query, "fieldSelector", fieldSelector)  # type String
    Swagger.set_param(_ctx.query, "labelSelector", labelSelector)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "resourceVersion", resourceVersion)  # type String
    Swagger.set_param(_ctx.query, "timeoutSeconds", timeoutSeconds)  # type Int32
    Swagger.set_param(_ctx.query, "watch", watch)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf", "application/json;stream=watch", "application/vnd.kubernetes.protobuf;stream=watch"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    return _ctx
end

"""

list or watch objects of kind NodeMetrics
Param: allowWatchBookmarks::Bool
Param: __continue__::String
Param: fieldSelector::String
Param: labelSelector::String
Param: limit::Int32
Param: resourceVersion::String
Param: timeoutSeconds::Int32
Param: watch::Bool
Return: IoK8sApiMetricsV1beta1NodeMetricsList
"""
function listMetricsV1beta1NodeMetrics(_api::MetricsV1beta1Api; allowWatchBookmarks=nothing, __continue__=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, resourceVersion=nothing, timeoutSeconds=nothing, watch=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_listMetricsV1beta1NodeMetrics(_api; allowWatchBookmarks=allowWatchBookmarks, __continue__=__continue__, fieldSelector=fieldSelector, labelSelector=labelSelector, limit=limit, resourceVersion=resourceVersion, timeoutSeconds=timeoutSeconds, watch=watch, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function listMetricsV1beta1NodeMetrics(_api::MetricsV1beta1Api, response_stream::Channel; allowWatchBookmarks=nothing, __continue__=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, resourceVersion=nothing, timeoutSeconds=nothing, watch=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_listMetricsV1beta1NodeMetrics(_api; allowWatchBookmarks=allowWatchBookmarks, __continue__=__continue__, fieldSelector=fieldSelector, labelSelector=labelSelector, limit=limit, resourceVersion=resourceVersion, timeoutSeconds=timeoutSeconds, watch=watch, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_listMetricsV1beta1PodMetrics(_api::MetricsV1beta1Api; allowWatchBookmarks=nothing, __continue__=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, resourceVersion=nothing, timeoutSeconds=nothing, watch=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApiMetricsV1beta1PodMetricsList, "/apis/metrics.k8s.io/v1beta1/pods", ["BearerToken"])
    Swagger.set_param(_ctx.query, "allowWatchBookmarks", allowWatchBookmarks)  # type Bool
    Swagger.set_param(_ctx.query, "continue", __continue__)  # type String
    Swagger.set_param(_ctx.query, "fieldSelector", fieldSelector)  # type String
    Swagger.set_param(_ctx.query, "labelSelector", labelSelector)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "resourceVersion", resourceVersion)  # type String
    Swagger.set_param(_ctx.query, "timeoutSeconds", timeoutSeconds)  # type Int32
    Swagger.set_param(_ctx.query, "watch", watch)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf", "application/json;stream=watch", "application/vnd.kubernetes.protobuf;stream=watch"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    return _ctx
end

"""

list or watch objects of kind PodMetrics
Param: allowWatchBookmarks::Bool
Param: __continue__::String
Param: fieldSelector::String
Param: labelSelector::String
Param: limit::Int32
Param: resourceVersion::String
Param: timeoutSeconds::Int32
Param: watch::Bool
Return: IoK8sApiMetricsV1beta1PodMetricsList
"""
function listMetricsV1beta1PodMetrics(_api::MetricsV1beta1Api; allowWatchBookmarks=nothing, __continue__=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, resourceVersion=nothing, timeoutSeconds=nothing, watch=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_listMetricsV1beta1PodMetrics(_api; allowWatchBookmarks=allowWatchBookmarks, __continue__=__continue__, fieldSelector=fieldSelector, labelSelector=labelSelector, limit=limit, resourceVersion=resourceVersion, timeoutSeconds=timeoutSeconds, watch=watch, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function listMetricsV1beta1PodMetrics(_api::MetricsV1beta1Api, response_stream::Channel; allowWatchBookmarks=nothing, __continue__=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, resourceVersion=nothing, timeoutSeconds=nothing, watch=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_listMetricsV1beta1PodMetrics(_api; allowWatchBookmarks=allowWatchBookmarks, __continue__=__continue__, fieldSelector=fieldSelector, labelSelector=labelSelector, limit=limit, resourceVersion=resourceVersion, timeoutSeconds=timeoutSeconds, watch=watch, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_readMetricsV1beta1NamespacedPodMetrics(_api::MetricsV1beta1Api, name::String, namespace::String; pretty=nothing, exact=nothing, __export__=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApiMetricsV1beta1PodMetrics, "/apis/metrics.k8s.io/v1beta1/namespaces/{namespace}/pods/{name}", ["BearerToken"])
    Swagger.set_param(_ctx.path, "name", name)  # type String
    Swagger.set_param(_ctx.path, "namespace", namespace)  # type String
    Swagger.set_param(_ctx.query, "pretty", pretty)  # type String
    Swagger.set_param(_ctx.query, "exact", exact)  # type Bool
    Swagger.set_param(_ctx.query, "export", __export__)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf", "application/json;stream=watch", "application/vnd.kubernetes.protobuf;stream=watch"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    return _ctx
end

"""

read the specified PodMetrics
Param: name::String (required)
Param: namespace::String (required)
Param: pretty::String
Param: exact::Bool
Param: __export__::Bool
Return: IoK8sApiMetricsV1beta1PodMetrics
"""
function readMetricsV1beta1NamespacedPodMetrics(_api::MetricsV1beta1Api, name::String, namespace::String; pretty=nothing, exact=nothing, __export__=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_readMetricsV1beta1NamespacedPodMetrics(_api, name, namespace; pretty=pretty, exact=exact, __export__=__export__, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function readMetricsV1beta1NamespacedPodMetrics(_api::MetricsV1beta1Api, response_stream::Channel, name::String, namespace::String; pretty=nothing, exact=nothing, __export__=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_readMetricsV1beta1NamespacedPodMetrics(_api, name, namespace; pretty=pretty, exact=exact, __export__=__export__, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_readMetricsV1beta1NodeMetrics(_api::MetricsV1beta1Api, name::String; pretty=nothing, exact=nothing, __export__=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApiMetricsV1beta1NodeMetrics, "/apis/metrics.k8s.io/v1beta1/nodes/{name}", ["BearerToken"])
    Swagger.set_param(_ctx.path, "name", name)  # type String
    Swagger.set_param(_ctx.query, "pretty", pretty)  # type String
    Swagger.set_param(_ctx.query, "exact", exact)  # type Bool
    Swagger.set_param(_ctx.query, "export", __export__)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    return _ctx
end

"""

read the specified NodeMetrics
Param: name::String (required)
Param: pretty::String
Param: exact::Bool
Param: __export__::Bool
Return: IoK8sApiMetricsV1beta1NodeMetrics
"""
function readMetricsV1beta1NodeMetrics(_api::MetricsV1beta1Api, name::String; pretty=nothing, exact=nothing, __export__=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_readMetricsV1beta1NodeMetrics(_api, name; pretty=pretty, exact=exact, __export__=__export__, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function readMetricsV1beta1NodeMetrics(_api::MetricsV1beta1Api, response_stream::Channel, name::String; pretty=nothing, exact=nothing, __export__=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_readMetricsV1beta1NodeMetrics(_api, name; pretty=pretty, exact=exact, __export__=__export__, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export listMetricsV1beta1NodeMetrics, listMetricsV1beta1PodMetrics, readMetricsV1beta1NamespacedPodMetrics, readMetricsV1beta1NodeMetrics


text:
using PowerModelsDistribution
using PowerModels
using Ipopt
using JuMP
using DataFrames
using CSV
# using Memento

include("support_functions.jl")

_PMs = PowerModels

function constraint_tp_load_current_delta_mod(pm::_PMs.GenericPowerModel{T}, nw::Int, load_id::Int, load_bus_id::Int, cp::_PMs.MultiConductorVector, cq::_PMs.MultiConductorVector) where T <: _PMs.AbstractACPForm
    cp_ab, cp_bc, cp_ca = cp
    cq_ab, cq_bc, cq_ca = cq
    vm_a, vm_b, vm_c = [_PMs.var(pm, nw, c, :vm, load_bus_id) for c in 1:3]
    va_a, va_b, va_c = [_PMs.var(pm, nw, c, :va, load_bus_id) for c in 1:3]
    # v_xy = v_x - v_y
    vre_xy(vm_x, va_x, vm_y, va_y) = JuMP.@NLexpression(pm.model, vm_x*cos(va_x)-vm_y*cos(va_y))
    vim_xy(vm_x, va_x, vm_y, va_y) = JuMP.@NLexpression(pm.model, vm_x*sin(va_x)-vm_y*sin(va_y))
    vre_ab = vre_xy(vm_a, va_a, vm_b, va_b)
    vim_ab = vim_xy(vm_a, va_a, vm_b, va_b)
    vre_bc = vre_xy(vm_b, va_b, vm_c, va_c)
    vim_bc = vim_xy(vm_b, va_b, vm_c, va_c)
    vre_ca = vre_xy(vm_c, va_c, vm_a, va_a)
    vim_ca = vim_xy(vm_c, va_c, vm_a, va_a)
    # i_xy = conj(s_xy/v_xy)
    ire_xy(cp_xy, cq_xy, vre_xy, vim_xy) = JuMP.@NLexpression(pm.model, 1/sqrt(vre_xy^2+vim_xy^2)*(cp_xy*vre_xy+cq_xy*vim_xy))
    iim_xy(cp_xy, cq_xy, vre_xy, vim_xy) = JuMP.@NLexpression(pm.model, 1/sqrt(vre_xy^2+vim_xy^2)*(cp_xy*vim_xy-cq_xy*vre_xy))
    ire_ab = ire_xy(cp_ab, cq_ab, vre_ab, vim_ab)
    iim_ab = iim_xy(cp_ab, cq_ab, vre_ab, vim_ab)
    ire_bc = ire_xy(cp_bc, cq_bc, vre_bc, vim_bc)
    iim_bc = iim_xy(cp_bc, cq_bc, vre_bc, vim_bc)
    ire_ca = ire_xy(cp_ca, cq_ca, vre_ca, vim_ca)
    iim_ca = iim_xy(cp_ca, cq_ca, vre_ca, vim_ca)
    # s_x = v_x*conj(i_xy-i_zx)
    # p_x = vm_x*cos(va_x)*(ire_xy-ire_zx) + vm_x*sin(va_x)*(iim_xy-iim_zx)
    # q_x = vm_x*sin(va_x)*(ire_xy-ire_zx) - vm_x*cos(va_x)*(iim_xy-iim_zx)
    p_x(vm_x, va_x, ire_xy, iim_xy, ire_zx, iim_zx) = JuMP.@NLexpression(pm.model, vm_x*cos(va_x)*(ire_xy-ire_zx) + vm_x*sin(va_x)*(iim_xy-iim_zx))
    q_x(vm_x, va_x, ire_xy, iim_xy, ire_zx, iim_zx) = JuMP.@NLexpression(pm.model, vm_x*sin(va_x)*(ire_xy-ire_zx) - vm_x*cos(va_x)*(iim_xy-iim_zx))
    # s_x = s_x,ref
    _PMs.var(pm, nw, 1, :pd)[load_id] = p_x(vm_a, va_a, ire_ab, iim_ab, ire_ca, iim_ca)
    _PMs.var(pm, nw, 2, :pd)[load_id] = p_x(vm_b, va_b, ire_bc, iim_bc, ire_ab, iim_ab)
    _PMs.var(pm, nw, 3, :pd)[load_id] = p_x(vm_c, va_c, ire_ca, iim_ca, ire_bc, iim_bc)
    _PMs.var(pm, nw, 1, :qd)[load_id] = q_x(vm_a, va_a, ire_ab, iim_ab, ire_ca, iim_ca)
    _PMs.var(pm, nw, 2, :qd)[load_id] = q_x(vm_b, va_b, ire_bc, iim_bc, ire_ab, iim_ab)
    _PMs.var(pm, nw, 3, :qd)[load_id] = q_x(vm_c, va_c, ire_ca, iim_ca, ire_bc, iim_bc)
end

function constraint_tp_load_impedance_delta_mod(pm::_PMs.GenericPowerModel{T}, nw::Int, load_id::Int, load_bus_id::Int, cp::_PMs.MultiConductorVector, cq::_PMs.MultiConductorVector) where T <: _PMs.AbstractACPForm
    cp_ab, cp_bc, cp_ca = cp
    cq_ab, cq_bc, cq_ca = cq
    vm_a, vm_b, vm_c = [_PMs.var(pm, nw, c, :vm, load_bus_id) for c in 1:3]
    va_a, va_b, va_c = [_PMs.var(pm, nw, c, :va, load_bus_id) for c in 1:3]
    # v_xy = v_x - v_y
    vre_xy(vm_x, va_x, vm_y, va_y) = JuMP.@NLexpression(pm.model, vm_x*cos(va_x)-vm_y*cos(va_y))
    vim_xy(vm_x, va_x, vm_y, va_y) = JuMP.@NLexpression(pm.model, vm_x*sin(va_x)-vm_y*sin(va_y))
    vre_ab = vre_xy(vm_a, va_a, vm_b, va_b)
    vim_ab = vim_xy(vm_a, va_a, vm_b, va_b)
    vre_bc = vre_xy(vm_b, va_b, vm_c, va_c)
    vim_bc = vim_xy(vm_b, va_b, vm_c, va_c)
    vre_ca = vre_xy(vm_c, va_c, vm_a, va_a)
    vim_ca = vim_xy(vm_c, va_c, vm_a, va_a)
    # i_xy = conj(s_xy/v_xy)
    ire_xy(cp_xy, cq_xy, vre_xy, vim_xy) = JuMP.@NLexpression(pm.model, cp_xy*vre_xy+cq_xy*vim_xy)
    iim_xy(cp_xy, cq_xy, vre_xy, vim_xy) = JuMP.@NLexpression(pm.model, cp_xy*vim_xy-cq_xy*vre_xy)
    ire_ab = ire_xy(cp_ab, cq_ab, vre_ab, vim_ab)
    iim_ab = iim_xy(cp_ab, cq_ab, vre_ab, vim_ab)
    ire_bc = ire_xy(cp_bc, cq_bc, vre_bc, vim_bc)
    iim_bc = iim_xy(cp_bc, cq_bc, vre_bc, vim_bc)
    ire_ca = ire_xy(cp_ca, cq_ca, vre_ca, vim_ca)
    iim_ca = iim_xy(cp_ca, cq_ca, vre_ca, vim_ca)
    # s_x = v_x*conj(i_xy-i_zx)
    # p_x = vm_x*cos(va_x)*(ire_xy-ire_zx) + vm_x*sin(va_x)*(iim_xy-iim_zx)
    # q_x = vm_x*sin(va_x)*(ire_xy-ire_zx) - vm_x*cos(va_x)*(iim_xy-iim_zx)
    p_x(vm_x, va_x, ire_xy, iim_xy, ire_zx, iim_zx) = JuMP.@NLexpression(pm.model, vm_x*cos(va_x)*(ire_xy-ire_zx) + vm_x*sin(va_x)*(iim_xy-iim_zx))
    q_x(vm_x, va_x, ire_xy, iim_xy, ire_zx, iim_zx) = JuMP.@NLexpression(pm.model, vm_x*sin(va_x)*(ire_xy-ire_zx) - vm_x*cos(va_x)*(iim_xy-iim_zx))
    # s_x = s_x,ref
    _PMs.var(pm, nw, 1, :pd)[load_id] = p_x(vm_a, va_a, ire_ab, iim_ab, ire_ca, iim_ca)
    _PMs.var(pm, nw, 2, :pd)[load_id] = p_x(vm_b, va_b, ire_bc, iim_bc, ire_ab, iim_ab)
    _PMs.var(pm, nw, 3, :pd)[load_id] = p_x(vm_c, va_c, ire_ca, iim_ca, ire_bc, iim_bc)
    _PMs.var(pm, nw, 1, :qd)[load_id] = q_x(vm_a, va_a, ire_ab, iim_ab, ire_ca, iim_ca)
    _PMs.var(pm, nw, 2, :qd)[load_id] = q_x(vm_b, va_b, ire_bc, iim_bc, ire_ab, iim_ab)
    _PMs.var(pm, nw, 3, :qd)[load_id] = q_x(vm_c, va_c, ire_ca, iim_ca, ire_bc, iim_bc)
end

function constraint_load_current_wye_mod(pm::_PMs.GenericPowerModel{T}, nw::Int, cnd::Int, load_id::Int, load_bus_id::Int, scale_p::Real, scale_q::Real) where T <: _PMs.AbstractACPForm
    vm = _PMs.var(pm, nw, cnd, :vm, load_bus_id)
    # this has to be a NLexpression and not an expression;
    # might be mixed in KCL with NLexpression, so has to also be NLexpression
    # original expression
    # _PMs.var(pm, nw, cnd, :pd)[load_id] = JuMP.@NLexpression(pm.model, scale_p*vm)
    # _PMs.var(pm, nw, cnd, :qd)[load_id] = JuMP.@NLexpression(pm.model, scale_q*vm)
    # end of original expression
    # _PMs.var(pm, nw, cnd, :pd)[load_id] = JuMP.@NLexpression(pm.model, scale_p*vm*1.02)
    # _PMs.var(pm, nw, cnd, :qd)[load_id] = JuMP.@NLexpression(pm.model, scale_q*vm*1.02)
    _PMs.var(pm, nw, cnd, :pd)[load_id] = JuMP.@variable(pm.model)
    _PMs.var(pm, nw, cnd, :qd)[load_id] = JuMP.@variable(pm.model)
    JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :pd)[load_id] == scale_p*vm)
    JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :qd)[load_id] == scale_q*vm)

    # JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :pd)[load_id] <=  scale_p*vm*(1.02))
    # JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :qd)[load_id] <=  scale_q*vm*(1.02))
    # JuMP.@constraint(pm.model, scale_p*vm*(0.98) <= _PMs.var(pm, nw, cnd, :pd)[load_id])
    # JuMP.@constraint(pm.model, scale_q*vm*(0.98) <= _PMs.var(pm, nw, cnd, :qd)[load_id])
    # JuMP.@NLexpression(pm.model, scale_q*vm*(0.98)) <=
    # JuMP.@NLexpression(pm.model, scale_p*vm*(0.98)) <=
end

function constraint_load_impedance_wye_mod(pm::_PMs.GenericPowerModel{T}, nw::Int, cnd::Int, load_id::Int, load_bus_id::Int, cp::Real, cq::Real) where T <: _PMs.AbstractACPForm
    vm = _PMs.var(pm, nw, cnd, :vm, load_bus_id)
    _PMs.var(pm, nw, cnd, :pd)[load_id] = JuMP.@NLexpression(pm.model, cp*vm^2)
    _PMs.var(pm, nw, cnd, :qd)[load_id] = JuMP.@NLexpression(pm.model, cq*vm^2)
end

# this is the only function really required
function constraint_load_power_wye_mod(pm::_PMs.GenericPowerModel{T}, nw::Int, cnd::Int, load_id::Int, pd::Real, qd::Real, adc_flex) where T <: _PMs.AbstractACPForm
    # _PMs.var(pm, nw, cnd, :pd)[load_id] = pd
    # _PMs.var(pm, nw, cnd, :qd)[load_id] = qd

    # println("we are in constatn load power y model")
    # @show load_id
    _PMs.var(pm, nw, cnd, :pd)[load_id] = JuMP.@variable(pm.model, start = pd)
    _PMs.var(pm, nw, cnd, :qd)[load_id] = JuMP.@variable(pm.model, start = qd)

    load_bus = pm.data["load"][string(load_id)]["load_bus"]
    if string(load_bus) in keys(adc_flex)
        flex_region = adc_flex[string(load_bus)]*1e-5
        # @show flex_region

#        println("Splitting flex region into phases ...")

        pd_original = pm.data["load"][string(load_id)]["pd"]
        qd_original = pm.data["load"][string(load_id)]["qd"]

        pmin = flex_region[1,2]
        pmax = flex_region[1,1]
        qmin = flex_region[2,2]
        qmax = flex_region[2,1]
        #
        pmin = pmin*pd_original[cnd]/sum(pd_original)
        pmax = pmax*pd_original[cnd]/sum(pd_original)
        qmin = qmin*qd_original[cnd]/sum(qd_original)
        qmax = qmax*qd_original[cnd]/sum(qd_original)
        #
        # pmin = pd*0.8
        # pmax = pd*1.2
        # qmin = qd*0.8
        # qmax = qd*1.2
        #
        # pmin = pd
        # pmax = pd
        # qmin = qd
        # qmax = qd

        # @show pd, qd
        # if pmin > pd || pmax < pd || qmin > qd || qmax < qd
        #     println("Infeasible constraints!!!")
        # end

        # # println("Adding flex constraints for ADC ", string(load_bus))
        JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :pd)[load_id] >=  pmin)
        JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :pd)[load_id] <=  pmax)
        JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :qd)[load_id] >=  qmin)
        JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :qd)[load_id] <=  qmax)
    else
        JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :pd)[load_id] ==  pd)
        JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :qd)[load_id] ==  qd)
    end

    # JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :pd)[load_id] <=  pd*1.02)
    # JuMP.@constraint(pm.model, _PMs.var(pm, nw, cnd, :qd)[load_id] <=  qd*1.02)
    # JuMP.@constraint(pm.model, pd*0.98 <= _PMs.var(pm, nw, cnd, :pd)[load_id])
    # JuMP.@constraint(pm.model, qd*0.98 <= _PMs.var(pm, nw, cnd, :qd)[load_id])
end

function constraint_tp_load_power_delta_mod(pm::_PMs.GenericPowerModel{T}, nw::Int, load_id::Int, load_bus_id::Int, pd::_PMs.MultiConductorVector, qd::_PMs.MultiConductorVector) where T <: _PMs.AbstractACPForm
    p_ab, p_bc, p_ca = pd
    q_ab, q_bc, q_ca = qd
    vm_a, vm_b, vm_c = [_PMs.var(pm, nw, c, :vm, load_bus_id) for c in 1:3]
    va_a, va_b, va_c = [_PMs.var(pm, nw, c, :va, load_bus_id) for c in 1:3]
    # v_xy = v_x - v_y
    vre_xy(vm_x, va_x, vm_y, va_y) = JuMP.@NLexpression(pm.model, vm_x*cos(va_x)-vm_y*cos(va_y))
    vim_xy(vm_x, va_x, vm_y, va_y) = JuMP.@NLexpression(pm.model, vm_x*sin(va_x)-vm_y*sin(va_y))
    vre_ab = vre_xy(vm_a, va_a, vm_b, va_b)
    vim_ab = vim_xy(vm_a, va_a, vm_b, va_b)
    vre_bc = vre_xy(vm_b, va_b, vm_c, va_c)
    vim_bc = vim_xy(vm_b, va_b, vm_c, va_c)
    vre_ca = vre_xy(vm_c, va_c, vm_a, va_a)
    vim_ca = vim_xy(vm_c, va_c, vm_a, va_a)
    # i_xy = conj(s_xy/v_xy)
    ire_xy(p_xy, q_xy, vre_xy, vim_xy) = JuMP.@NLexpression(pm.model, (p_xy*vre_xy+q_xy*vim_xy)/(vre_xy^2+vim_xy^2))
    iim_xy(p_xy, q_xy, vre_xy, vim_xy) = JuMP.@NLexpression(pm.model, (p_xy*vim_xy-q_xy*vre_xy)/(vre_xy^2+vim_xy^2))
    ire_ab = ire_xy(p_ab, q_ab, vre_ab, vim_ab)
    iim_ab = iim_xy(p_ab, q_ab, vre_ab, vim_ab)
    ire_bc = ire_xy(p_bc, q_bc, vre_bc, vim_bc)
    iim_bc = iim_xy(p_bc, q_bc, vre_bc, vim_bc)
    ire_ca = ire_xy(p_ca, q_ca, vre_ca, vim_ca)
    iim_ca = iim_xy(p_ca, q_ca, vre_ca, vim_ca)
    # s_x = v_x*conj(i_xy-i_zx)
    # p_x = vm_x*cos(va_x)*(ire_xy-ire_zx) + vm_x*sin(va_x)*(iim_xy-iim_zx)
    # q_x = vm_x*sin(va_x)*(ire_xy-ire_zx) - vm_x*cos(va_x)*(iim_xy-iim_zx)
    p_x(vm_x, va_x, ire_xy, iim_xy, ire_zx, iim_zx) = JuMP.@NLexpression(pm.model, vm_x*cos(va_x)*(ire_xy-ire_zx) + vm_x*sin(va_x)*(iim_xy-iim_zx))
    q_x(vm_x, va_x, ire_xy, iim_xy, ire_zx, iim_zx) = JuMP.@NLexpression(pm.model, vm_x*sin(va_x)*(ire_xy-ire_zx) - vm_x*cos(va_x)*(iim_xy-iim_zx))
    # s_x = s_x,ref
    _PMs.var(pm, nw, 1, :pd)[load_id] = p_x(vm_a, va_a, ire_ab, iim_ab, ire_ca, iim_ca)
    _PMs.var(pm, nw, 2, :pd)[load_id] = p_x(vm_b, va_b, ire_bc, iim_bc, ire_ab, iim_ab)
    _PMs.var(pm, nw, 3, :pd)[load_id] = p_x(vm_c, va_c, ire_ca, iim_ca, ire_bc, iim_bc)
    _PMs.var(pm, nw, 1, :qd)[load_id] = q_x(vm_a, va_a, ire_ab, iim_ab, ire_ca, iim_ca)
    _PMs.var(pm, nw, 2, :qd)[load_id] = q_x(vm_b, va_b, ire_bc, iim_bc, ire_ab, iim_ab)
    _PMs.var(pm, nw, 3, :qd)[load_id] = q_x(vm_c, va_c, ire_ca, iim_ca, ire_bc, iim_bc)
end

function constraint_tp_load_mod(pm::_PMs.GenericPowerModel, id::Int, adc_flex; nw=pm.cnw)
    load = _PMs.ref(pm, nw, :load, id)
    model = load["model"]
    conn = _PMs.ref(pm, nw, :load, id, "conn")
    @assert(conn in ["delta", "wye"])

    if model=="constant_power"
        pd = load["pd"]
        qd = load["qd"]

        if conn=="wye"
            for c in _PMs.conductor_ids(pm)
                constraint_load_power_wye_mod(pm, nw, c, id, pd[c], qd[c], adc_flex)
            end
        elseif conn=="delta"
            @assert(_PMs.ref(pm, 0, :conductors)==3)
            constraint_tp_load_power_delta_mod(pm, nw, id, load["load_bus"], pd, qd)
        end

    elseif model=="constant_current"
        vnom_kv = load["vnom_kv"]
        vbase_kv_LL = _PMs.ref(pm, nw, :bus, load["load_bus"])["base_kv"]
        vbase_kv_LN = vbase_kv_LL/sqrt(3)

        pd = load["pd"]
        qd = load["qd"]
        cp = pd/(vnom_kv/vbase_kv_LN)
        cq = qd/(vnom_kv/vbase_kv_LN)

        if conn=="wye"
            for c in _PMs.conductor_ids(pm)
                constraint_load_current_wye_mod(pm, nw, c, id, load["load_bus"], cp[c], cq[c])
            end
        elseif conn=="delta"
            @assert(_PMs.ref(pm, 0, :conductors)==3)
            constraint_tp_load_current_delta_mod(pm, nw, id, load["load_bus"], cp, cq)
        end

    elseif model=="constant_impedance"
        vnom_kv = load["vnom_kv"]
        vbase_kv_LL = _PMs.ref(pm, nw, :bus, load["load_bus"])["base_kv"]
        vbase_kv_LN = vbase_kv_LL/sqrt(3)

        pd = load["pd"]
        qd = load["qd"]
        cp = pd/(vnom_kv/vbase_kv_LN)^2
        cq = qd/(vnom_kv/vbase_kv_LN)^2

        if conn=="wye"
            for c in _PMs.conductor_ids(pm)
                constraint_load_impedance_wye_mod(pm, nw, c, id, load["load_bus"], cp[c], cq[c])
            end
        elseif conn=="delta"
            @assert(_PMs.ref(pm, 0, :conductors)==3)
            constraint_tp_load_impedance_delta_mod(pm, nw, id, load["load_bus"], cp, cq)
        end

    else
        # Memento.@error(LOGGER, "Unknown model $model for load $id.")
    end
end

function post_tp_opf_lm_mod(pm::_PMs.GenericPowerModel)
    variable_tp_voltage(pm)
    variable_tp_branch_flow(pm)

    for c in _PMs.conductor_ids(pm)
        _PMs.variable_generation(pm, cnd=c)
        variable_load(pm, cnd=c)
        _PMs.variable_dcline_flow(pm, cnd=c)
    end
    variable_tp_trans_flow(pm)

    constraint_tp_model_voltage(pm)

    for i in _PMs.ids(pm, :ref_buses)
        constraint_tp_theta_ref(pm, i)
    end

    # loads should be constrained before KCL, or Pd/Qd undefined
    adc_flex = build_adc_flex_model(filename_flex, pm.data)
    for id in _PMs.ids(pm, :load)
        constraint_tp_load_mod(pm, id, adc_flex)
    end

    for i in _PMs.ids(pm, :bus), c in _PMs.conductor_ids(pm)
        constraint_tp_power_balance_shunt_trans_load(pm, i, cnd=c)
    end

    for i in _PMs.ids(pm, :branch)
        for c in _PMs.conductor_ids(pm)
            constraint_tp_ohms_yt_from(pm, i, cnd=c)
            constraint_tp_ohms_yt_to(pm, i, cnd=c)

            _PMs.constraint_voltage_angle_difference(pm, i, cnd=c)

            _PMs.constraint_thermal_limit_from(pm, i, cnd=c)
            _PMs.constraint_thermal_limit_to(pm, i, cnd=c)
        end
    end

    for i in _PMs.ids(pm, :dcline), c in _PMs.conductor_ids(pm)
        _PMs.constraint_dcline(pm, i, cnd=c)
    end

    for i in _PMs.ids(pm, :trans)
        constraint_tp_trans(pm, i)
    end

    # _PMs.objective_min_fuel_cost(pm)
    JuMP.@NLobjective(pm.model,Min, 100*(_PMs.var(pm, 0, 1, :pg)[1] + _PMs.var(pm, 0, 2, :pg)[1] + _PMs.var(pm, 0, 3, :pg)[1])
        + 1e5*sum((_PMs.var(pm, 0, i, :pd)[load_id] - pm.data["load"][string(load_id)]["pd"][i])^2  for load_id in _PMs.ids(pm, :load), i=1:3)
        + 1e5*sum((_PMs.var(pm, 0, i, :qd)[load_id] - pm.data["load"][string(load_id)]["qd"][i])^2  for load_id in _PMs.ids(pm, :load), i=1:3)
        )
end

# result_new = _PMs.run_model(data_new, ACPPowerModel, with_optimizer(Ipopt.Optimizer), post_tp_opf_lm_mod; multiconductor=true, ref_extensions=[ref_add_arcs_trans!])

# Changing load bounds

max_stars_repo_name: areiman/PFO-ADC-DER-Testbed

max_issues_repo_name: areiman/PFO-ADC-DER-Testbed

max_forks_repo_name: areiman/PFO-ADC-DER-Testbed


text:
function jacobian(fun, x, h)
    values = fun.(x)
    N = length(values)
    M = zeros(N, N)

    for i=1:length(values)
        M[i, i] = imag(fun(x[i] + im*h))/h
    end
    return M
end


text:
# bounds.jl

"""
   bound_remainder(::Type{TaylorModel1}f::Function, polf::Taylor1, polfI::Taylor1, x0::Interval, I::Interval)

Bound the absolute remainder of the polynomial approximation of `f` given
by the Taylor polynomial `polf` around `x0` on the interval `I`. It requires
the interval extension `polfI` of the polynomial that approximates `f` for
the whole interval `I`, in order to compute the Lagrange remainder.

If `polfI[end]` has a definite sign, then it is monotonic in the intervals
[I.lo, x0] and [x0.hi, I.hi], which is exploited; otherwise, it is used
to compute the Lagrange remainder. This corresponds to Prop 2.2.1 in Mioara
Joldes PhD thesis (pp 52).

"""
function bound_remainder(::Type{TaylorModel1}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)

    _order = get_order(polf) + 1
    fTIend = polfI[_order]
    if (sup(fTIend) < 0 || inf(fTIend) > 0)
        # Absolute remainder is monotonic
        a = interval(I.lo)
        b = interval(I.hi)
        Δlo = f(a) - polf(a-x0)
        # Δlo = f(a) - bound_taylor1(polf, a-x0)
        Δhi = f(b) - polf(b-x0)
        # Δhi = f(b) - bound_taylor1(polf, b-x0)
        Δx0 = f(x0) - polf[0]
        Δ = hull(hull(Δlo, Δx0), Δhi)
    else
        # Lagrange bound
        Δ = fTIend * (I-x0)^_order
    end
    return Δ
end

"""
   bound_remainder(::Type{RTaylorModel1}, f::Function, polf::Taylor1, polfI::Taylor1, x0::Interval, I::Interval)

Bound the relative remainder of the polynomial approximation of `f` given
by the Taylor polynomial `polf` around `x0` on the interval `I`. It requires
an the interval extension `polfI` of a polynomial that approximates `f` for
the whole interval `I`, in order to compute the Lagrange remainder.

If `polfI[end]` has a definite sign, then it is monotonic in the interval `I`,
which is exploited; otherwise, the last coefficients bounds the relative
remainder. This corresponds to Prop 2.3.7 in Mioara Joldes' PhD thesis (pp 67).

"""
function bound_remainder(::Type{RTaylorModel1}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)

    _order = get_order(polf) + 1
    fTIend = polfI[_order+1]
    a = interval(I.lo)
    b = interval(I.hi)
    if (sup(fTIend) < 0 || inf(fTIend) > 0) && isempty(a ∩ x0) && isempty(b ∩ x0)
        # Error is monotonic
        denom_lo = (a-x0)^_order
        Δlo = f(a) - polf(a-x0)
        # Δlo = f(a) - bound_taylor1(polf, a-x0)
        Δlo = Δlo / denom_lo
        denom_hi = (b-x0)^_order
        Δhi = f(b) - polf(b-x0)
        # Δhi = f(b) - bound_taylor1(polf, b-x0)
        Δhi = Δhi / denom_hi
        Δ = hull(Δlo, Δhi)
    else
        # Lagrange coefficient
        Δ = polfI[_order]
    end
    return Δ
end

"""
    bound_taylor1(fT::Taylor1, I::Interval)

Compute a *tight* polynomial bound for the Taylor polynomial `fT`
in the interval `I`.

Note: Algorithm 2.1.1 corresponds to `evaluate(fT, I)` or simply `fT(I).
This function uses the roots of the derivative of `ft` to obtain a
tighter bound.

"""
function bound_taylor1(fT::Taylor1, I::Interval)

    # Check if the fT is monotonous (the derivative has a given sign)
    fTd  = TaylorSeries.derivative(fT)
    range_deriv = fTd(I)
    (sup(range_deriv) ≤ 0 || inf(range_deriv) ≥ 0) && return bound_taylor1(fT, fTd, I)

    # Compute roots of the derivative using the second derivative
    # Fix some sort of relative tolerance for Newton root search
    fTd2 = TaylorSeries.derivative(fTd)
    rootsder = roots(x->fTd(x), x->fTd2(x), I, Newton, 1.0e-5*mag(I))

    # Bound the range of fT using the roots and end points
    num_roots = length(rootsder)
    num_roots == 0 && return fT(I)
    rangepoly = hull( fT(interval(I.lo)), fT(interval(I.hi)) )
    for ind in 1:num_roots
        rangepoly = hull(rangepoly, fT(rootsder[ind].interval))
    end

    return rangepoly
end

"""
    bound_taylor1(fT::Taylor1, fTd::Taylor1, I::Interval)

Compute a *tight* polynomial bound for the Taylor polynomial `fT`
in the interval `I`, considering whether its derivative `ftd` has
a definite sign.

"""
function bound_taylor1(fT::Taylor1{T}, fTd::Taylor1{T}, I::Interval{T}) where {T}
    #
    if inf(fTd(I)) ≥ 0
        return Interval(fT(I.lo), fT(I.hi))
    elseif sup(fTd(I)) ≤ 0
        return Interval(fT(I.hi), fT(I.lo))
    end
    return fT(I)
end
function bound_taylor1(fT::Taylor1{Interval{T}}, fTd::Taylor1{Interval{T}},
        I::Interval{S}) where {T,S}
    #
    if inf(fTd(I)) ≥ 0
        return hull(fT(I.lo), fT(I.hi))
    elseif sup(fTd(I)) ≤ 0
        return hull(fT(I.hi), fT(I.lo))
    end
    return fT(I)
end

"""
    bound_taylor1(fT::TaylorModel1, I=domain(fT)::Interval)

Compute a *tight* polynomial bound for the Taylor model `fT`
in the interval `I`, considering whether its derivative `ftd` has
a definite sign.

"""
bound_taylor1(fT::TaylorModel1, I=domain(fT)::Interval) = bound_taylor1(polynomial(fT), I)

"""
    linear_dominated_bounder(fT::TaylorModel1, ϵ=1e-3::Float64, max_iter=5::Int)

Compute a tighter polynomial bound for the Taylor model `fT` by the linear
dominated bounder algorithm. The linear dominated algorithm is applied until
the bound of `fT` gets tighter than `ϵ` or the number of steps reachs `max_iter`.
The returned bound corresponds to the improved polynomial bound with the remainder
of the `TaylorModel1` included.
"""
function linear_dominated_bounder(fT::TaylorModel1{T, S}; ϵ=1e-3, max_iter=5) where {T, S}
    d = 1.
    Dn = fT.dom
    Dm = fT.x0
    Pm = Taylor1(copy(fT.pol.coeffs))
    bound = interval(0.)
    n_iter = 0
    while d > ϵ && n_iter < max_iter
        x0 = mid(Dn - Dm)
        c = mid(Dn)
        update!(Pm, x0)
        linear = Taylor1(Pm.coeffs[1:2], Pm.order)
        non_linear = Pm - linear
        if T <: Interval
            Li = mid(linear[1])
        else
            Li = linear[1]
        end
        I1 = linear(Dn - c)
        Ih = non_linear(Dn - c)
        bound = I1.lo + Ih
        d = diam(bound)
        n_iter += 1
        if Li == 0
            break
        elseif Li > 0
            new_hi = min(Dn.lo + (d / abs(Li)), Dn.hi)
            Dm = Dn
            Dn = interval(Dn.lo, new_hi)
        else
            new_lo = max(Dn.hi - (d / abs(Li)), Dn.lo)
            Dm = Dn
            Dn = interval(new_lo, Dn.hi)
        end
    end
    hi = fT.pol(fT.dom - fT.x0).hi
    return interval(bound.lo, hi) + fT.rem
end

"""
    linear_dominated_bounder(fT::TaylorModelN, ϵ=1e-3::Float64, max_iter=5::Int)

Compute a tighter polynomial bound for the Taylor model `fT` by the linear
dominated bounder algorithm. The linear dominated algorithm is applied until
the bound of `fT` gets tighter than `ϵ` or the number of steps reachs `max_iter`.
The returned bound corresponds to the improved polynomial bound with the remainder
of the `TaylorModelN` included.
"""
function linear_dominated_bounder(fT::TaylorModelN{N, T, S}; ϵ=1e-5, max_iter=5) where {N, T, S}
    d = 1.
    Dn = fT.dom
    Dm = fT.x0
    Pm = TaylorN(deepcopy(fT.pol.coeffs))
    bound = interval(0.)
    n_iter = 0
    x0 = Array{S}(undef, N)
    linear_coeffs = Array{Float64}(undef, N)
    while d > ϵ && n_iter < max_iter
        x0 .= mid(Dn - Dm)
        c = mid(Dn)
        update!(Pm, x0)
        linear_part = Pm[1]
        if T <: Interval
            linear_coeffs .= mid.(linear_part.coeffs)
        else
            linear_coeffs .= linear_part.coeffs
        end
        linear = TaylorN(deepcopy(Pm.coeffs[1:2]), Pm.order)
        non_linear = Pm - linear
        centered_domain = Dn .- c
        I1 = linear(centered_domain)
        Ih = non_linear(centered_domain)
        bound = I1.lo + Ih
        d = diam(bound)
        n_iter += 1
        new_boxes = Interval[]
        for (idx, box) in enumerate(Dn)
            Li = linear_coeffs[idx]
            if Li == 0
                Dni = box
            elseif Li < 0
                lo = max(box.hi - (d / abs(Li)), box.lo)
                Dni = interval(lo, box.hi)
            else
                hi = min(box.lo + (d / abs(Li)), box.hi)
                Dni = interval(box.lo, hi)
            end
            push!(new_boxes, Dni)
        end
        Dm = Dn
        Dn = IntervalBox(new_boxes...)
    end
    hi = fT.pol(fT.dom - fT.x0).hi
    return interval(bound.lo, hi) + fT.rem
end

"""
    quadratic_fast_bounder(fT::TaylorModel1)

Compute a *tighter* polynomial bound by the quadratic fast bounder.
The returned bound corresponds to the "improved" polynomial bound
with the remainder of the `TaylorModel1` included. This "improved" bound
can be one of the following two:
    1) An improved bound: if the domain of `fT` has a local minimizer,
       then an improved bound is returned.
    2) Original TaylorModel bound: if the local minimizer criteria is not
       satisfied, then the original bound of `fT` is returned.

This algorithm is a slight modification to the Makino & Berz algorithm.
For this algorithm the linear part is bounded by solving a simple
set of linear equations (compared to the iterative procedure by Makino & Berz).
"""
function quadratic_fast_bounder(fT::TaylorModel1)
    P = fT.pol
    bound_tm = fT(fT.dom - fT.x0)
    if signbit(P[2])
        return bound_tm
    else
        x0 = -P[1] / (2 * P[2])
        x = Taylor1(fT.pol.order)
        Qx0 = (x - x0) * P[2] * (x - x0)
        bound_qfb = (P - Qx0)(fT.dom - fT.x0)
        hi = P(fT.dom - fT.x0).hi
        bound_qfb = interval(bound_qfb.lo, hi) + fT.rem
        bound = bound_qfb ∩ bound_tm
        return bound
    end
end

"""
    quadratic_fast_bounder(fT::TaylorModelN)

Compute a *tighter* polynomial bound by the quadratic fast bounder.
The returned bound corresponds to the "improved" polynomial bound
with the remainder of the `TaylorModelN` included. This "improved" bound
can be one of the following two:
    1) An improved bound: if the domain of `fT` has a local minimizer,
       then an improved bound is returned.
    2) Original TaylorModel bound: if the local minimizer criteria is not
       satisfied, then the original bound of `fT` is returned.

This algorithm is a slight modification to the Makino & Berz algorithm.
For this algorithm the linear part is bounded by solving a simple
set of linear equations (compared to the iterative procedure by Makino & Berz).
"""
function quadratic_fast_bounder(fT::TaylorModelN)
    P = fT.pol
    H = Matrix(TaylorSeries.hessian(P))
    bound_tm = fT(fT.dom - fT.x0)
    if isposdef(H)
        P1 = -P[1].coeffs
        xn = H \ P1
        x = set_variables("x", numvars=length(xn))
        Qxn = 0.5 * (x - xn)' * H * (x - xn)
        bound_qfb = (P - Qxn)(fT.dom - fT.x0)
        hi = P(fT.dom - fT.x0).hi
        bound_qfb = interval(bound_qfb.lo, hi) + fT.rem
        bound = bound_qfb ∩ bound_tm
        return bound
    else
        return bound_tm
    end
end


text:
## NOTE - Uses code from src/a_star_implicit_constrained.jl

struct AStarEpsilonMCSPHEntry{D <: Number}
    a_star_entry::AStarMCSPHEntry{D}
    focal_heuristic::D
end

# Base.isless(e1::AStarEpsilonMCSPHEntry, e2::AStarEpsilonMCSPHEntry) = (e1.a_star_entry.fvalue, e1.a_star_entry.gvalue) < (e2.a_star_entry.fvalue, e2.a_star_entry.gvalue)

compare(comp::CompareFocalHeap, e1::AStarEpsilonMCSPHEntry, e2::AStarEpsilonMCSPHEntry) = (e1.focal_heuristic, e1.a_star_entry.fvalue, -e1.a_star_entry.gvalue) < (e2.focal_heuristic, e2.a_star_entry.fvalue, -e2.a_star_entry.gvalue)

@with_kw mutable struct AStarEpsilonMCSPStates{D <: Number}
    a_star_states::AStarMCSPStates{D}                                           = AStarMCSPStates{D}()
    # Subset of heap
    focal_heap::MutableBinaryHeap{AStarEpsilonMCSPHEntry{D},CompareFocalHeap}   = MutableBinaryHeap{AStarEpsilonMCSPHEntry{D},CompareFocalHeap}()
    ol_handles_in_focal::Set{Int64}                                             = Set{Int64}()
    focal_to_ol_handle::Dict{Int64,Int64}                                       = Dict{Int64,Int64}()
    best_fvalue::D = zero(D)
    focal_heur_value::Dict{Int64,D}                                         = Dict{Int64,D}()
end

function process_neighbors_implicit!(
    state::AStarEpsilonMCSPStates{D},
    graph::AbstractGraph{V},
    edge_wt_fn::F1,
    neighbors::Vector{Int64},
    parent_entry::AStarEpsilonMCSPHEntry{D},
    entry_cl_handle::Int64,    # For closed list
    visitor::AbstractDijkstraVisitor,
    eps_weight::Float64,
    admissible_heuristic::F2,
    focal_state_heuristic::F3,
    focal_transition_heuristic::F4,
    weight_functions::Vector{F5},
    weight_heuristics::Vector{F6},
    weight_constraints::Vector{D}
    ) where {V, D <: Number, F1 <: Function, F2 <: Function, F3 <: Function, F4 <: Function, F5 <: Function, F6 <: Function}

    for nbr_idx in neighbors

        is_elig, new_weight_vector =  is_eligible(graph, parent_entry.a_star_entry, nbr_idx, weight_functions, weight_heuristics, weight_constraints)

        if is_elig
            # Create new entry for open list
            # Updated gvalues, fvalues, weights AND focal stuff
            new_gvalue = parent_entry.a_star_entry.gvalue + edge_wt_fn(graph.vertices[parent_entry.a_star_entry.v_idx], graph.vertices[nbr_idx])
            new_fvalue = new_gvalue + admissible_heuristic(graph.vertices[nbr_idx])
            new_focal_heur = parent_entry.focal_heuristic + focal_state_heuristic(graph.vertices[nbr_idx]) +
                             focal_transition_heuristic(graph.vertices[parent_entry.a_star_entry.v_idx], graph.vertices[nbr_idx])

            new_entry = AStarEpsilonMCSPHEntry(AStarMCSPHEntry(nbr_idx, parent_entry.a_star_entry.v_idx, new_gvalue, new_fvalue, new_weight_vector, entry_cl_handle),
                                                new_focal_heur)

            # Insert into open list if not dominated
            # And accordingly into focal list if valid
            if ~(is_dominated(state.a_star_states, new_entry.a_star_entry))
                new_entry_handle = push!(state.a_star_states.open_list, new_entry.a_star_entry)
                state.focal_heur_value[new_entry_handle] = new_focal_heur

                if ~(haskey(state.a_star_states.open_list_hmap, nbr_idx))
                    Graphs.discover_vertex!(visitor, graph.vertices[parent_entry.a_star_entry.v_idx], graph.vertices[nbr_idx], new_gvalue)
                    state.a_star_states.open_list_hmap[nbr_idx] = Set{Int64}(new_entry_handle)
                else
                    push!(state.a_star_states.open_list_hmap[nbr_idx], new_entry_handle)
                end

                # Now enter to focal list if valid according to weight
                if new_fvalue <= eps_weight * state.best_fvalue
                    new_focal_handle = push!(state.focal_heap, new_entry)
                    push!(state.ol_handles_in_focal, new_entry_handle)
                    state.focal_to_ol_handle[new_focal_handle] = new_entry_handle
                end # end if
            end # end if
        end
    end
end

"""
Runs the Focal-MCSP (Focal search with Multiple Constraints bounded-suboptimal Shortest Path) algorithm on the given list-of-indices graph.
As with all A* versions (except a_star_spath), the graph is implicit and the graph visitor must generate the neighbors of the expanded
vertex just-in-time.

Arguments:
    - `graph::AbstractGraph{V} where {V}`
    - `edge_wt_fn::F1 where {F1 <: Function}` Maps (u,v) to the edge weight of the u -> v edge
    - `source::Int64`
    - `visitor::AbstractDijkstraVisitor` The graph visitor that implements include_vertex!, which is called when
                                         a vertex is expanded
    - `eps_weight::Float64` w sub-optimality factor for focal search
    - `admissible_heuristic::F2 where {F2 <: Function}` Maps u to h(u); an admissible heuristic for the cost-to-go
    - `focal_state_heuristic::F3 where {F3 <: Function}`
    - `focal_transition_heuristic::F4 where {F4 <: Function}`
    - `weight_functions::Vector{F5} where {F5 <: Function}` A list of constraint weight functions for edges
    - `weight_heuristics::Vector{F6} where {F6 <: Function}` A list of admissible heuristics for the weight constraint for edges
    - `weight_constraints::Vector{D} where {D <: Number}` A list of correspoding weight constraint values
    - `::Type{D} = Float64` The type of the edge weight value

Returns:
    - `::AStarMCSPEpsilonStates` The result of the A* search
    - `::AStarMCSP` The MCSP heap entry for the target node (or the source if no path found)
"""
function a_star_implicit_constrained_epsilon_path!(
    graph::AbstractGraph{V},                # the graph
    edge_wt_fn::F1, # distances associated with edges
    source::Int64,             # the source index
    visitor::AbstractDijkstraVisitor,# visitor object\
    eps_weight::Float64,
    admissible_heuristic::F2,      # Heuristic function for vertices
    focal_state_heuristic::F3,
    focal_transition_heuristic::F4,
    weight_functions::Vector{F5},
    weight_heuristics::Vector{F6},
    weight_constraints::Vector{D},
    ::Type{D} = Float64
    ) where {V, D <: Number, F1 <: Function, F2 <: Function, F3 <: Function, F4 <: Function, F5 <: Function, F6 <: Function}

    state = AStarEpsilonMCSPStates{D}()

    n_weights = length(weight_constraints)

    source_heur = admissible_heuristic(graph.vertices[source])
    state.best_fvalue = source_heur

    source_entry = AStarEpsilonMCSPHEntry(AStarMCSPHEntry(source, source, zero(D), source_heur, zeros(D, n_weights), 0),
                                         focal_state_heuristic(source))

    source_handle = push!(state.a_star_states.open_list, source_entry.a_star_entry)
    state.a_star_states.open_list_hmap[source] = Set{Int64}(source_handle)
    state.focal_heur_value[source_handle] = focal_state_heuristic(source)

    focal_handle = push!(state.focal_heap, source_entry)
    push!(state.ol_handles_in_focal, source_handle)
    state.focal_to_ol_handle[focal_handle] = source_handle

    while ~(isempty(state.a_star_states.open_list))

        old_best_fvalue = state.best_fvalue
        state.best_fvalue = top(state.a_star_states.open_list).fvalue

        # Update focal list
        if state.best_fvalue > old_best_fvalue
            # Iterate over open set  in increasing order of fvalue and insert in focal list if valid
            for node in sort(state.a_star_states.open_list.nodes, by = x->x.value.fvalue)
                fvalue = node.value.fvalue

                if fvalue > eps_weight * old_best_fvalue && fvalue <= eps_weight * state.best_fvalue &&
                    ~(node.handle in state.ol_handles_in_focal)
                    # Need to insert open list entry to focal list
                    # And assign node handle to focal handle
                    new_focal_entry = AStarEpsilonMCSPHEntry(node.value, state.focal_heur_value[node.handle])
                    focal_handle = push!(state.focal_heap, new_focal_entry)
                    push!(state.ol_handles_in_focal, node.handle)
                    state.focal_to_ol_handle[focal_handle] = node.handle
                end

                if fvalue > eps_weight * state.best_fvalue
                    break
                end
            end
        end

        focal_entry, focal_handle = top_with_handle(state.focal_heap)

        # Get corresponding open list handle and entry
        ol_handle = state.focal_to_ol_handle[focal_handle]

        # @show state.focal_heap

        # Remove from open list and focal list
        pop!(state.focal_heap)
        delete!(state.focal_to_ol_handle, focal_handle) # Don't really need this?
        delete!(state.a_star_states.open_list, ol_handle)
        delete!(state.a_star_states.open_list_hmap[focal_entry.a_star_entry.v_idx], ol_handle)

        nbrs = Vector{Int64}(undef, 0)

        if Graphs.include_vertex!(visitor, graph.vertices[focal_entry.a_star_entry.parent_idx],
                                  graph.vertices[focal_entry.a_star_entry.v_idx], focal_entry.a_star_entry.gvalue, nbrs) == false
            return state, focal_entry.a_star_entry
        end

        # Insert into closed list to get new handle
        cl_entry_handle = push!(state.a_star_states.closed_list, focal_entry.a_star_entry)
        if ~(haskey(state.a_star_states.closed_list_hmap, focal_entry.a_star_entry.v_idx))
            state.a_star_states.closed_list_hmap[focal_entry.a_star_entry.v_idx] = Set{Int64}(cl_entry_handle)
        else
            push!(state.a_star_states.closed_list_hmap[focal_entry.a_star_entry.v_idx], cl_entry_handle)
        end

        # Process nbrs
        process_neighbors_implicit!(state, graph, edge_wt_fn, nbrs,
                                    focal_entry, cl_entry_handle, visitor, eps_weight,
                                    admissible_heuristic, focal_state_heuristic, focal_transition_heuristic,
                                    weight_functions, weight_heuristics, weight_constraints)
        Graphs.close_vertex!(visitor, graph.vertices[focal_entry.a_star_entry.v_idx])
    end

    return state, source_entry.a_star_entry
end


text:
# Power electronics Problem

# 30
function SOPWM_for_3_level_invereters(x)
    D = length(x)
    ## SOPWM for 3-level Invereters
    m = 0.32;
    s = (-ones(25)) .^ (2:26);
    k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];
    ## objective function
    f = zeros(2);
    su = 0;
    for j = 1:31
        su2 = 0;
        for l = 1:D
            su2 = su2 + s[l] .* cos(k[j] .* x[l]*pi/180);
        end
        su = su + su2 .^ 2 ./ k[j] .^ 4;
    end
    f[1] = (su) .^ 0.5 ./ (sum(1 ./ k .^ 4)) .^ 0.5;
    f[2] = (sum(s .* cos.(x*pi/180))-m) .^ 2;
    ## constraints
    g = zeros(D-1);
    for i = 1:D-1
        g[i] = x[i]-x[i+1]+1e-6;
    end
    h = zeros(1);

    return f, g, h
end

# 31
function SOPWM_for_5_level_Inverters(x)
    D = length(x)
    ## SOPWM for 5-level Inverters
    m = 0.32;
    s = [1,-1,1,1,-1,1,-1,1,-1,-1,1,-1,1,1,-1,1,-1,1,-1,-1,1,-1,1,1,-1];
    k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];
    ## objective function
    f = zeros(2);
    su = 0;
    for j = 1:31
        su2 = 0;
        for l = 1:D
            su2 = su2 + s[l] .* cos(k[j] .* x[l]*pi/180);
        end
        su = su + su2 .^ 2 ./ k[j] .^ 4;
    end
    f[1] = (su) .^ 0.5 ./ (sum(1 ./ k .^ 4)) .^ 0.5;
    f[2] = (sum(s .* cos.(x*pi/180))-m) .^ 2;
    ## constraints
    g = zeros(D-1);
    for i = 1:D-1
        g[i] = x[i]-x[i+1]+1e-6;
    end
    h = zeros(1);
    return f, g, h
end

# 32
function SOPWM_for_7_level_inverters(x)
    D = length(x)
    ## SOPWM for 7-level Inverters
    m = 0.36;
    s = [1,-1,1,1,1,-1,-1,-1,1,1,-1,-1,1,1,1,-1,-1,-1,1,1,-1,-1,1,1,1];
    k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];
    ## objective function
    f = zeros(2);
    su = 0;
    for j = 1:31
        su2 = 0;
        for l = 1:D
            su2 = su2 + s[l] .* cos(k[j] .* x[l]*pi/180);
        end
        su = su + su2 .^ 2 ./ k[j] .^ 4;
    end
    f[1] = (su) .^ 0.5 ./ (sum(1 ./ k .^ 4)) .^ 0.5;
    f[2] = (sum(s .* cos.(x*pi/180))-m) .^ 2;
    ## constraints
    g = zeros(D-1);
    for i = 1:D-1
        g[i] = x[i]-x[i+1]+1e-6;
    end
    h = zeros(1);
    return f, g, h
end

# 33
function SOPWM_for_9_level_inverters(x)
    D = length(x)
    ## SOPWM for 9-level Inverters
    m = 0.32;
    s = [1,1,1,1,-1,1,-1,-1,-1,1,-1,-1,1,1,1,1,-1,1,-1,-1,-1,1,-1,-1,1,1,1,1,-1,1];
    k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];
    ## objective function
    f = zeros(2);
    su = 0;
    for j = 1:31
        su2 = 0;
        for l = 1:D
            su2 = su2 + s[l] .* cos(k[j] .* x[l]*pi/180);
        end
        su = su + su2 .^ 2 ./ k[j] .^ 4;
    end
    f[1] = (su) .^ 0.5 ./ (sum(1 ./ k .^ 4)) .^ 0.5;
    f[2] = (sum(s .* cos.(x*pi/180))-m) .^ 2;
    ## constraints
    g = zeros(D-1);
    for i = 1:D-1
        g[i] = x[i]-x[i+1]+1e-6;
    end
    h = zeros(1);

    return f, g, h
end

# 34
function SOPWM_for_11_level_inverters(x)
    D = length(x)
    ## SOPWM for 11-level Inverters
    m = 0.3333;
    s = [1,-1,1,1,1,-1,-1,-1,1,1,1,1,-1,-1,1,-1,-1,-1,1,1,1,1,-1,1,1,-1,-1,1,-1,-1];
    k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];
    ## objective function
    f = zeros(2);
    su = 0;
    for j = 1:31
        su2 = 0;
        for l = 1:D
            su2 = su2 + s[l] .* cos(k[j] .* x[l]*pi/180);
        end
        su = su + su2 .^ 2 ./ k[j] .^ 4;
    end
    f[1] = (su) .^ 0.5 ./ (sum(1 ./ k .^ 4)) .^ 0.5;
    f[2] = (sum(s .* cos.(x*pi/180))-m) .^ 2;
    ## constraints
    g = zeros(D-1);
    for i = 1:D-1
        g[i] = x[i]-x[i+1]+1e-6;
    end

    h = zeros(1);
    return f, g, h
end

# 35
function SOPWM_for_13_level_inverters(x)
    D = length(x)
    ## SOPWM for 13-level Inverters
    m = 0.32;
    s = [1,1,1,-1,1,-1,1,-1,1,1,1,1,-1,-1,-1,-1,1,-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,1];
    k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];

    ## objective function
    f = zeros(2);
    su = 0;
    for j = 1:31
        su2 = 0;
        for l = 1:D
            su2 = su2 + s[l] .* cos(k[j] .* x[l]*pi/180);
        end
        su = su + su2 .^ 2 ./ k[j] .^ 4;
    end
    f[1] = (su) .^ 0.5 ./ (sum(1 ./ k .^ 4)) .^ 0.5;
    f[2] = (sum(s .* cos.(x*pi/180))-m) .^ 2;

    ## constraints
    g = zeros(D-1);
    for i = 1:D-1
        g[i] = x[i]-x[i+1]+1e-6;
    end
    h = zeros(1);

    return f, g, h
end


text:
"""
The `mode` type parameter must be either `:inflate` or `:deflate`.
"""
type Source{mode,T<:BufferedInputStream}
    input::T
    zstream::ZStream
    state::State
    reset_on_end::Bool
end

# inflate source constructors
# ---------------------------

function InflateSource{T<:BufferedInputStream}(input::T, gzip::Bool,
                                               reset_on_end::Bool)
    return Source{:inflate,T}(
        input,
        init_inflate_zstream(gzip),
        initialized,
        reset_on_end)
end

function InflateSource(input::BufferedInputStream, bufsize::Integer, gzip::Bool,
                       reset_on_end::Bool)
    return InflateSource(input, gzip, reset_on_end)
end

function InflateSource(input::IO, bufsize::Integer, gzip::Bool, reset_on_end::Bool)
    input_stream = BufferedInputStream(input, bufsize)
    return InflateSource(input_stream, gzip, reset_on_end)
end

function InflateSource(input::Vector{UInt8}, bufsize::Integer, gzip::Bool,
                       reset_on_end::Bool)
    return InflateSource(BufferedInputStream(input), gzip, reset_on_end)
end

"""
    ZlibInflateInputStream(input[; <keyword arguments>])

Construct a zlib inflate input stream to decompress gzip/zlib data.

# Arguments
* `input`: a byte vector, IO object, or BufferedInputStream containing compressed data to inflate.
* `bufsize::Integer=8192`: input and output buffer size.
* `gzip::Bool=true`: if true, data is gzip compressed; if false, zlib compressed.
* `reset_on_end::Bool=true`: on stream end, try to find the start of another stream.
"""
function ZlibInflateInputStream(input; bufsize::Integer=8192, gzip::Bool=true,
                                reset_on_end::Bool=true)
    return BufferedInputStream(
        InflateSource(input, bufsize, gzip, reset_on_end),
        bufsize)
end

# deflate source constructors
# ---------------------------

function DeflateSource{T<:BufferedInputStream}(
        input::T, gzip::Bool, level::Integer, mem_level::Integer, strategy)
    return Source{:deflate,T}(
        input,
        init_deflate_zstream(gzip, level, mem_level, strategy),
        initialized,
        false)
end

function DeflateSource(input::BufferedInputStream, bufsize::Integer, gzip::Bool,
                       level::Integer, mem_level::Integer, strategy)
    return DeflateSource(input, gzip, level, mem_level, strategy)
end

function DeflateSource(input::IO, bufsize::Integer, gzip::Bool, level::Integer,
                       mem_level::Integer, strategy)
    input_stream = BufferedInputStream(input, bufsize)
    return DeflateSource(input_stream, gzip, level, mem_level, strategy)
end

function DeflateSource(input::Vector{UInt8}, bufsize::Integer, gzip::Bool,
                       level::Integer, mem_level::Integer, strategy)
    return DeflateSource(BufferedInputStream(input), gzip, level, mem_level, strategy)
end

"""
    ZlibDeflateInputStream(input[; <keyword arguments>])

Construct a zlib deflate input stream to compress gzip/zlib data.

# Arguments
* `input`: a byte vector, IO object, or BufferedInputStream containing data to compress.
* `bufsize::Integer=8192`: input and output buffer size.
* `gzip::Bool=true`: if true, data is gzip compressed; if false, zlib compressed.
* `level::Integer=6`: compression level in 1-9.
* `mem_level::Integer=8`: memory to use for compression in 1-9.
* `strategy=Z_DEFAULT_STRATEGY`: compression strategy; see zlib documentation.
"""
function ZlibDeflateInputStream(input;
                                bufsize::Integer=8192,
                                gzip::Bool=true,
                                level::Integer=6,
                                mem_level::Integer=8,
                                strategy=Z_DEFAULT_STRATEGY)
    return BufferedInputStream(
        DeflateSource(input, bufsize, gzip, level, mem_level, strategy), bufsize)
end

"""
    readbytes!(source, buffer, from, to)

Read bytes from the zlib stream to a buffer. Satisfies the BufferedStreams
source interface.
"""
function BufferedStreams.readbytes!{mode}(
        source::Source{mode},
        buffer::Vector{UInt8},
        from::Int, to::Int)
    if source.state == finalized
        return 0
    elseif source.state == finished && source.reset_on_end
        reset!(source)
    end

    @trans source (
        initialized => inprogress,
        inprogress  => inprogress
    )

    fillbuffer!(source.input)
    source.zstream.next_out = pointer(buffer, from)
    source.zstream.avail_out = to - from + 1
    _, n_out = process(
        source,
        mode == :deflate && eof(source.input) ? Z_FINISH : Z_NO_FLUSH)
    return n_out
end

function process{mode}(source::Source{mode}, flush)
    @assert source.state == inprogress
    # counter of processed input/output bytes
    n_in = n_out = 0
    input = source.input
    zstream = source.zstream

    #println("--- Source{", mode, "} ---")
    @label process
    zstream.next_in = pointer(input)
    zstream.avail_in = BufferedStreams.available_bytes(input)
    old_avail_in = zstream.avail_in
    old_avail_out = zstream.avail_out
    if mode == :inflate
        ret = inflate!(zstream, flush)
    else
        ret = deflate!(zstream, flush)
    end
    n_in += old_avail_in - zstream.avail_in
    n_out += old_avail_out - zstream.avail_out
    input.position += old_avail_in - zstream.avail_in

    if ret == Z_OK
        if zstream.avail_in == 0
            if BufferedStreams.fillbuffer!(input) == 0
                flush = Z_FINISH
            end
            @goto process
        end
    elseif ret == Z_STREAM_END
        @trans source inprogress => finished
    elseif ret == Z_BUF_ERROR
        # could not consume more input or produce more output
    elseif ret < 0
        zerror(zstream, ret)
    else
        @assert false
    end

    return n_in, n_out
end

@inline function Base.eof{mode}(source::Source{mode})
    if source.state == initialized ||
        (mode == :inflate && source.state == finished && source.reset_on_end)
        return eof(source.input)
    end
    return source.state == finished || source.state == finalized
end

function Base.close{mode}(source::Source{mode})
    if source.state == finalized
        isopen(source.input) && close(source.input)
        return
    end
    if mode == :inflate
        @zcheck end_inflate!(source.zstream)
    else
        @zcheck end_deflate!(source.zstream)
    end
    @trans source (
        initialized => finalized,
        inprogress  => finalized,
        finished    => finalized
    )
    close(source.input)
    return
end

function reset!{mode}(source::Source{mode})
    if mode == :inflate
        @zcheck reset_inflate!(source.zstream)
    else
        @zcheck reset_deflate!(source.zstream)
    end
    @trans source (
        initialized => initialized,
        inprogress  => initialized,
        finished    => initialized,
        finalized   => initialized
    )
    return source
end


text: load_path_function() = 1

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
function ntd()

end


text:
"Update the header bounding box and count based on point data"
function update!(h::LasHeader, pvec::Vector{T}) where T <: LasPoint
    x_min, y_min, z_min = Inf, Inf, Inf
    x_max, y_max, z_max = -Inf, -Inf, -Inf
    for p in pvec
        x, y, z = xcoord(p, h), ycoord(p, h), zcoord(p, h)
        if x < x_min
            x_min = x
        end
        if y < y_min
            y_min = y
        end
        if z < z_min
            z_min = z
        end
        if x > x_max
            x_max = x
        end
        if y > y_max
            y_max = y
        end
        if z > z_max
            z_max = z
        end
    end
    h.x_min = x_min
    h.y_min = y_min
    h.z_min = z_min
    h.x_max = x_max
    h.y_max = y_max
    h.z_max = z_max
    h.records_count = length(pvec)
    nothing
end


text:
"""

Contains embedding of data under diffusion maps.
"""
struct DiffMap
    Y  # m x n
    L
    Λ
    Φ
end

"""

"""
function DiffMap(X; m::Integer=2, σ = 2)
    D, n = size(X)

    # # Step 1: Compute similarity matrix
    𝒦 = gaussian_similarity(X, σ)

    # Step 2: compute laplacian
    L = rlap(𝒦)

    # Step 3: Compute eigenvectors of L
    # Arpack returns complex-valued eigenvalues.
    Λ, Φ = eigs(L, nev=m+1, which=:SR)
    Λ = real(Λ)
    Φ = real(Φ)

    # Remember that (1,...,1) is an eigenvector with eigenvalue 0. This is uninformative,
    # hence we discard it as usual.
    Φ ./= Φ[:, 1]
    Φ = Φ[:, 2:m+1]
    Λ = Λ[2:m+1]
    Y = permutedims(Φ .* reshape(Λ, (1, :)))
    return DiffMap(Y, L, Λ, Φ)
end

getm(DM::DiffMap) = size(DM.Y, 1)
getn(DM::DiffMap) = size(DM.Y, 2)
size(DM::DiffMap) = size(DM.Y)
size(DM::DiffMap, d::Integer) = size(DM.Y, d)

"""
    Return measure of density for laplacian. (#nonzeros / #total)
"""
sparsity(DM::DiffMap) = sparsity(DM.L)
sparsity(A) = sum(iszero.(A)) / length(A)

function Base.iterate(DM::DiffMap, state = 1)
    if state > 4
        return nothing
    else
        getfield(DM, state), state + 1
    end
end