text:
# deals with CCP's crest API
module CREST
using Requests
using DataFrames
using Currencies
using ProgressMeter
using TimeSeries

using eve.SDE
using eve.MarketDatas
@usingcustomcurrency isk "Interstellar Kredit" 2
export getLatestSellOrders
export getLatestBuyOrders
export getLatestOrders
export MarketData
global bucket = 400
global last_update = 0.0
global connections = 20
global connectionChanged = Condition()

function addTokens()
  global bucket
  global last_update
  t = time()
  e = t - last_update
  tok = e * 150.0
  bucket += tok
  bucket = min(bucket, 400)
  last_update = t
end

macro ratelimit(f)
  quote
    global bucket
    global connections
    addTokens()
    while bucket < 1.0 || connections < 1
      if bucket < 1.0
        sleep(rand((1/150):1.0))
      elseif connections < 1
        wait(connectionChanged)
      end
      addTokens()
    end
    bucket -= 1.0
    connections -= 1.0
    r = $(esc(f))
    connections += 1.0
    notify(connectionChanged, all=false)
    r
  end
end

const publicEndpoint = "https://public-crest.eveonline.com/"
const authEndpoint = "https://crest-tq.eveonline.com/"

let
  global endpoints
  crestresp = ""
  function endpoints()
    if crestresp != ""
      return crestresp
    end
    r = @ratelimit get(publicEndpoint)
    crestresp = Requests.json(r)
    crestresp
  end
end
function refreshEndpoints()
  r = @ratelimit get(publicEndpoint)
  crestresp = Requests.json(r)
end
function regionsEndpoint()
  endpoints()["regions"]["href"]
end
function regionEndpoint(regionID :: Number)
  regionsEndpoint() * "$regionID/"
end
function regionEndpoint(regionName :: AbstractString)
  regionEndpoint(regionID(regionName))
end
function typeEndpoint(itemID)
  endpoints()["itemTypes"]["href"] * "$itemID/"
end
function marketBuyEndpoint(region :: Number)
  publicEndpoint * "market/$region/orders/buy/"
end
function marketSellEndpoint(region :: Number)
  publicEndpoint * "market/$region/orders/sell/"
end
function marketOrdersEndpoint(region :: Number)
  publicEndpoint * "market/$region/orders/"
end
function marketHistoryEndpoint(region :: Number, itemID :: Number)
  publicEndpoint * "market/$region/types/$itemID/history/"
end

function processMarketData(jsonData)
  count = jsonData["totalCount"]
  #result = DataFrame([typeof(isk), Int, Int, Int, Dates.Day, DateTime],
  #                   [:Price, :Location, :Quantity, :MinQuantity, :Duration, :Issued], count)
  result = DataFrame()
  df = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS")
  result[:Price] = map(x -> x["price"], jsonData["items"])
  result[:Location] = map(x -> x["location"]["id"], jsonData["items"])
  result[:Quantity] = map(x -> x["volume"], jsonData["items"])
  result[:MinQuantity] = map( x -> x["minVolume"], jsonData["items"])
  result[:Duration] = map(x -> Dates.Day(x["duration"]), jsonData["items"])
  result[:Issued] = map(x -> DateTime(x["issued"], df), jsonData["items"])
  result[:OrderID] = map(x -> x["id"], jsonData["items"])
  result[:Bid] = map(x -> x["buy"], jsonData["items"])
  result
end
function convertMarketHistory(df :: DataFrame)
  timestamps = map(DateTime, df[:Date])
  nrows, ncols = size(df)
  print(nrows)
  local values = Array(Float64, nrows, 7)
  values[:, 1] = df[:LowPrice]
  values[:, 2] = df[:HighPrice]
  values[:, 3] = df[:AvgPrice]
  values[:, 4] = df[:Volume]
  values[:, 5] = df[:OrderCount]
  values[:, 6] = df[:BuyVolume]
  values[:, 7] = df[:SellVolume]
  TimeArray(timestamps.data, values,
    ["LowPrice", "HighPrice", "AvgPrice", "Volume", "OrderCount", "BuyVolume", "SellVolume"])
end
function addAverageOrderVols(df :: DataFrame)
  function computevols(buy, sell, avg, vol)
    avg == buy && return [vol 0]
    avg == sell && return [0, vol]
    A = [buy sell; 1 1]
    B = [vol * avg; vol]
    X = A  \ B
    return [X[1] X[2]]
  end
  vols = map(computevols, df[:LowPrice].data,
                          df[:HighPrice].data,
                          df[:AvgPrice].data,
                          df[:Volume].data)
  vols = vcat(vols...)
  df[:BuyVolume] = vols[:, 1]
  df[:SellVolume] = vols[:, 2]
end
function processMarketHistoryData(jsonData)
  df = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS")
  result = DataFrame()
  result[:Date] = map(x -> DateTime(x["date"], df), jsonData["items"])
  #timestamps = map(x -> DateTime(x["date"], df), jsonData["items"])
  result[:LowPrice] = map(x -> x["lowPrice"], jsonData["items"])
  result[:HighPrice] = map(x -> x["highPrice"], jsonData["items"])
  result[:AvgPrice] = map(x -> x["avgPrice"], jsonData["items"])
  result[:Volume] = map(x -> x["volume"], jsonData["items"])
  result[:OrderCount] = map(x -> x["orderCount"], jsonData["items"])
  addAverageOrderVols(result)
  #print(values)
  convertMarketHistory(result)
end
function getLatestSellOrders(regionID, itemID)
  uri = marketSellEndpoint(regionID)
  typ = typeEndpoint(itemID)
  resp = @ratelimit get(uri, query = Dict("type" => typ))
  if resp.status != 200
    print(resp.status)
  end
  result = processMarketData(Requests.json(resp))
  sort!(result, cols=[:Price], rev=false)
  result
end
function getLatestBuyOrders(regionID, itemID)
  uri = marketBuyEndpoint(regionID)
  typ = typeEndpoint(itemID)
  resp = @ratelimit get(uri, query = Dict("type" => typ))
  if resp.status != 200
    print(resp.status)
  end
  result = processMarketData(Requests.json(resp))
  sort!(result, cols=[:Price], rev=true)
  result
end
function getLatestOrders(regionID, item :: Int)
  result = MarketData(item, now(Dates.UTC), DataFrame(), DataFrame())
  uri =
  resp = @ratelimit get()
  result.sells = getLatestSellOrders(regionID, item)
  result.buys = getLatestBuyOrders(regionID, item)
  result
end
function getLatestOrders{T<:Int}(regionID, items :: Array{T})
  result = Dict{Int64, MarketData}()
  p = Progress(length(items), 1, "Fetching Orders... ", 50)
  @sync for i in items
    @async begin
      result[i] = getLatestOrders(regionID, i)
      next!(p)
    end
  end
  result
end
function getMarketHistory(regionID, itemID)
  uri = marketHistoryEndpoint(regionID, itemID)
  resp = @ratelimit get(uri)
  if resp.status != 200
    print(resp.status)
  end
  processMarketHistoryData(Requests.json(resp))
end
function getMarketHistory{T<:Int}(regionID, items :: Array{T})
  result = Dict{Int64, TimeArray}()
  p = Progress(length(items), 1, "Fetching History... ", 50)
  @sync for i in items
    @async begin
      result[i] = getMarketHistory(regionID, i)
      next!(p)
    end
  end
  result
end

end


text:
module CBF

import ..FileFormats
import MathOptInterface

const MOI = MathOptInterface

MOI.Utilities.@model(Model,
    (MOI.Integer,),
    (),
    (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
        MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone, MOI.DualExponentialCone),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,)
)
function MOI.supports_constraint(
    ::Model{T}, ::Type{MOI.SingleVariable},
    ::Type{<:MOI.Utilities.SUPPORTED_VARIABLE_SCALAR_SETS{T}}) where T

    return false
end
function MOI.supports_constraint(
    ::Model, ::Type{MOI.SingleVariable}, ::Type{MOI.Integer})
    return true
end

struct Options end

get_options(m::Model) = get(m.ext, :CBF_OPTIONS, Options())

"""
    Model()

Create an empty instance of `FileFormats.CBF.Model`.
"""
function Model()
    model = Model{Float64}()
    model.ext[:CBF_OPTIONS] = Options()
    return model
end

Base.show(io::IO, ::Model) = print(io, "A Conic Benchmark Format (CBF) model")

# ==============================================================================
#
#   Base.write
#
# ==============================================================================

"""
    Base.write(io::IO, model::FileFormats.CBF.Model)

Write `model` to `io` in the Conic Benchmark Format.
"""
function Base.write(io::IO, model::Model)
    options = get_options(model)
    # Helper functions for MOI constraints.
    model_cons(con_func, con_set) = MOI.get(model,
        MOI.ListOfConstraintIndices{con_func, con_set}())
    con_function(con_idx) = MOI.get(model, MOI.ConstraintFunction(), con_idx)
    con_set(con_idx) = MOI.get(model, MOI.ConstraintSet(), con_idx)

    # CBF version number.
    println(io, "VER")
    println(io, 3)
    println(io)

    # Power cone parameters.
    powcone_cons = vcat(model_cons(MOI.VectorOfVariables,
        MOI.PowerCone{Float64}), model_cons(
        MOI.VectorAffineFunction{Float64}, MOI.PowerCone{Float64}))
    dpowcone_cons = vcat(model_cons(MOI.VectorOfVariables,
        MOI.DualPowerCone{Float64}), model_cons(
        MOI.VectorAffineFunction{Float64}, MOI.DualPowerCone{Float64}))
    for (cons, powcone_str) in
        ((powcone_cons, "POWCONES"), (dpowcone_cons, "POW*CONES"))
        if !isempty(cons)
            println(io, powcone_str)
            num_powcones = length(cons)
            println(io, num_powcones, " ", 2 * num_powcones)
            for powcone in cons
                println(io, 2)
                powcone_param = con_set(powcone).exponent
                println(io, powcone_param)
                println(io, 1.0 - powcone_param)
            end
            println(io)
        end
    end

    # Objective sense.
    println(io, "OBJSENSE")
    obj_sense = MOI.get(model, MOI.ObjectiveSense())
    if obj_sense == MOI.MAX_SENSE
        println(io, "MAX")
    else
        println(io, "MIN") # Includes the case of feasibility sense.
    end
    println(io)

    # Variables.
    num_var = MOI.get(model, MOI.NumberOfVariables())
    println(io, "VAR")
    println(io, num_var, " 1")
    println(io, "F ", num_var)
    println(io)

    # Integrality constraints.
    integer_cons = model_cons(MOI.SingleVariable, MOI.Integer)
    if length(integer_cons) > 0
        println(io, "INT")
        println(io, length(integer_cons))
        for con_idx in integer_cons
            println(io, con_function(con_idx).variable.value - 1)
        end
        println(io)
    end

    # Objective function terms.
    obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_function = MOI.get(model, MOI.ObjectiveFunction{obj_type}())
    if obj_type == MOI.SingleVariable
        println(io, "OBJACOORD")
        println(io, 1)
        println(io, obj_function.variable.value - 1, " ", 1.0)
        println(io)
    elseif obj_type == MOI.ScalarAffineFunction{Float64}
        if !isempty(obj_function.terms)
            println(io, "OBJACOORD")
            println(io, length(obj_function.terms))
            for t in obj_function.terms
                println(io, t.variable_index.value - 1, " ", t.coefficient)
            end
            println(io)
        end
        if !iszero(obj_function.constant)
            println(io, "OBJBCOORD")
            println(io, obj_function.constant)
            println(io)
        end
    else
        error("Objective function type $obj_type is not supported.")
    end

    # Non-PSD constraints.
    num_rows = 0 # Number of constraint rows.
    con_cones = Tuple{String, Int}[] # List of cone types/dimensions.
    acoord = Tuple{Int, Int, Float64}[] # Affine terms.
    bcoord = Tuple{Int, Float64}[] # Constant terms.
    powcone_idx = 0
    dpowcone_idx = 0

    for (set_type, cone_str) in (
        (MOI.Zeros, "L="),
        (MOI.Reals, "F"),
        (MOI.Nonnegatives, "L+"),
        (MOI.Nonpositives, "L-"),
        (MOI.SecondOrderCone, "Q"),
        (MOI.RotatedSecondOrderCone, "QR"),
        (MOI.ExponentialCone, "EXP"),
        (MOI.DualExponentialCone, "EXP*"),
        (MOI.PowerCone{Float64}, ""),
        (MOI.DualPowerCone{Float64}, ""),
        )
        for con_idx in model_cons(MOI.VectorOfVariables, set_type)
            vars = con_function(con_idx).variables
            if set_type in (MOI.ExponentialCone, MOI.DualExponentialCone)
                # Reverse order of indices.
                for vj in reverse(vars)
                    num_rows += 1
                    push!(acoord, (num_rows, vj.value, 1.0))
                end
            else
                for vj in vars
                    num_rows += 1
                    push!(acoord, (num_rows, vj.value, 1.0))
                end
            end
            if set_type == MOI.PowerCone{Float64}
                cone_str = "@$(powcone_idx):POW"
                powcone_idx += 1
            elseif set_type == MOI.DualPowerCone{Float64}
                cone_str = "@$(dpowcone_idx):POW*"
                dpowcone_idx += 1
            end
            push!(con_cones, (cone_str, MOI.dimension(con_set(con_idx))))
        end

        for con_idx in model_cons(MOI.VectorAffineFunction{Float64}, set_type)
            con_func = con_function(con_idx)
            con_dim = MOI.dimension(con_set(con_idx))
            if set_type in (MOI.ExponentialCone, MOI.DualExponentialCone)
                # Reverse order of indices.
                @assert con_dim == 3
                for t in con_func.terms
                    push!(acoord, (num_rows + 4 - t.output_index,
                        t.scalar_term.variable_index.value,
                        t.scalar_term.coefficient))
                end
                for row_constant in reverse(con_func.constants)
                    num_rows += 1
                    push!(bcoord, (num_rows, row_constant))
                end
            else
                for t in con_func.terms
                    push!(acoord, (num_rows + t.output_index,
                        t.scalar_term.variable_index.value,
                        t.scalar_term.coefficient))
                end
                for row_constant in con_func.constants
                    num_rows += 1
                    push!(bcoord, (num_rows, row_constant))
                end
            end
            if set_type == MOI.PowerCone{Float64}
                cone_str = "@$(powcone_idx):POW"
                powcone_idx += 1
            elseif set_type == MOI.DualPowerCone{Float64}
                cone_str = "@$(dpowcone_idx):POW*"
                dpowcone_idx += 1
            end
            push!(con_cones, (cone_str, con_dim))
        end
    end

    if !isempty(con_cones)
        println(io, "CON")
        println(io, num_rows, " ", length(con_cones))
        for (cone_str, con_dim) in con_cones
            println(io, cone_str, " ", con_dim)
        end
        println(io)
    end

    if !isempty(acoord)
        println(io, "ACOORD")
        println(io, length(acoord))
        for (row, var, coef) in acoord
            println(io, row - 1, " ", var - 1, " ", coef)
        end
        println(io)
    end

    if !isempty(bcoord)
        println(io, "BCOORD")
        println(io, length(bcoord))
        for (row, constant) in bcoord
            println(io, row - 1, " ", constant)
        end
        println(io)
    end

    # PSD constraints.
    psd_side_dims = Int[] # List of side dimensions.
    hcoord = Tuple{Int, Int, Int, Int, Float64}[] # Affine terms.
    dcoord = Tuple{Int, Int, Int, Float64}[] # Constant terms.

    for con_idx in model_cons(MOI.VectorOfVariables,
        MOI.PositiveSemidefiniteConeTriangle) # MOI variable constraints.
        side_dim = con_set(con_idx).side_dimension
        push!(psd_side_dims, side_dim)
        vars = con_function(con_idx).variables
        k = 0
        for i in 1:side_dim, j in 1:i
            k += 1
            push!(hcoord, (length(psd_side_dims), vars[k].value, i, j, 1.0))
        end
        @assert k == length(vars)
    end

    for con_idx in model_cons(MOI.VectorAffineFunction{Float64},
        MOI.PositiveSemidefiniteConeTriangle) # MOI affine constraints.
        side_dim = con_set(con_idx).side_dimension
        push!(psd_side_dims, side_dim)
        con_func = con_function(con_idx)
        for t in con_func.terms
            # Get (i,j) index in symmetric matrix from lower triangle index.
            i = div(1 + isqrt(8 * t.output_index - 7), 2)
            j = t.output_index - div(i * (i - 1), 2)
            push!(hcoord, (length(psd_side_dims),
                t.scalar_term.variable_index.value,
                i, j, t.scalar_term.coefficient))
        end
        k = 0
        for i in 1:side_dim, j in 1:i
            k += 1
            row_constant = con_func.constants[k]
            if !iszero(row_constant)
                push!(dcoord, (length(psd_side_dims), i, j, row_constant))
            end
        end
        @assert k == MOI.output_dimension(con_func)
    end

    if !isempty(psd_side_dims)
        println(io, "PSDCON")
        println(io, length(psd_side_dims))
        for side_dim in psd_side_dims
            println(io, side_dim)
        end
        println(io)
    end

    if !isempty(hcoord)
        println(io, "HCOORD")
        println(io, length(hcoord))
        for (psd_idx, var, i, j, coef) in hcoord
            println(io, psd_idx - 1, " ", var - 1, " ",
                i - 1, " ", j - 1, " ", coef)
        end
        println(io)
    end

    if !isempty(dcoord)
        println(io, "DCOORD")
        println(io, length(dcoord))
        for (psd_idx, i, j, constant) in dcoord
            println(io, psd_idx - 1, " ", i - 1, " ",
                j - 1, " ", constant)
        end
        println(io)
    end

    return
end

# ==============================================================================
#
#   MOI.read_from_file
#
# The CBF file format (version 3) is described at
# http://cblib.zib.de/doc/format3.pdf.
# Note CBF indices start at 0.
#
# ==============================================================================

# Convert a pair of row and column indices of a symmetric matrix into a vector
# index for the row-wise lower triangle
function mat_to_vec_idx(i::Int, j::Int)
    if i < j
        return div((j - 1) * j, 2) + i
    else
        return div((i - 1) * i, 2) + j
    end
end

function cbf_to_moi_cone(cone_str::AbstractString, cone_dim::Int)
    if cone_str == "F" # Free cones (redundant but add anyway).
        return MOI.Reals(cone_dim)
    elseif cone_str == "L=" # Zero cones.
        return MOI.Zeros(cone_dim)
    elseif cone_str == "L-" # Nonpositive cones.
        return MOI.Nonpositives(cone_dim)
    elseif cone_str == "L+" # Nonnegative cones.
        return MOI.Nonnegatives(cone_dim)
    elseif cone_str == "Q" # Second-order cones.
        @assert cone_dim >= 2
        return MOI.SecondOrderCone(cone_dim)
    elseif cone_str == "QR" # Rotated second-order cones.
        @assert cone_dim >= 3
        return MOI.RotatedSecondOrderCone(cone_dim)
    else
        error("CBF cone name $cone_str is not recognized or supported.")
    end
end

function powcone_to_moi_cone(cone_str::AbstractString,
    powcone_alphas::Vector{Vector{Float64}},
    dpowcone_alphas::Vector{Vector{Float64}})
    raw_powcone_info = split(cone_str[2:end], ":")
    powcone_idx = parse(Int, raw_powcone_info[1]) + 1
    powcone_type = raw_powcone_info[2]
    if powcone_type == "POW"
        alpha = powcone_alphas[powcone_idx]
        return MOI.PowerCone{Float64}(first(alpha)/sum(alpha))
    elseif powcone_type == "POW*"
        alpha = dpowcone_alphas[powcone_idx]
        return MOI.DualPowerCone{Float64}(first(alpha)/sum(alpha))
    else
        error("Failed to parse parametric cone $powcone_type information.")
    end
end

"""
    Base.read!(io::IO, model::FileFormats.CBF.Model)

Read `io` in the Conic Benchmark Format and store the result in `model`.
"""
function Base.read!(io::IO, model::Model)
    if !MOI.is_empty(model)
        error("Cannot read in file because model is not empty.")
    end

    scalar_vars = MOI.VariableIndex[]
    psd_vars = Vector{MOI.VariableIndex}[]
    obj_terms = MOI.ScalarAffineTerm{Float64}[]
    obj_constant = 0.0
    powcone_alphas = Vector{Float64}[]
    dpowcone_alphas = Vector{Float64}[]
    con_cones = Tuple{String, Int}[]
    row_terms = Vector{MOI.ScalarAffineTerm{Float64}}[]
    row_constants = Float64[]
    psd_side_dims = Int[]
    psd_row_starts = Int[]
    psd_row_terms = Vector{MOI.ScalarAffineTerm{Float64}}[]
    psd_row_constants = Float64[]

    while !eof(io)
        line = strip(readline(io))

        # Skip blank lines and comments (CBF comments start with #).
        if isempty(line) || startswith(line, "#")
            continue
        end

        # CBF version number.
        if line == "VER"
            ver = parse(Int, split(strip(readline(io)))[1])
            if ver < 1 || ver > 3
                error("CBF version number $ver is not yet supported.")
            end
            continue
        end

        # Power cone parameters.
        if line == "POWCONES" || line == "POW*CONES"
            alpha = (line == "POWCONES") ? powcone_alphas : dpowcone_alphas
            raw_powcone_info = split(strip(readline(io)))
            @assert length(raw_powcone_info) == 2
            num_powcone = parse(Int, raw_powcone_info[1])
            num_lines = parse(Int, raw_powcone_info[2])
            alpha_idx = 0
            for j in 1:num_powcone
                num_alphas = parse(Int, strip(readline(io)))
                if num_alphas != 2
                    error("Only 3-dimensional power cones are supported.")
                end
                push!(alpha, [parse(Float64, strip(readline(io)))
                    for k in 1:num_alphas])
                alpha_idx += num_alphas
            end
            @assert num_lines == alpha_idx
            continue
        end

        # Objective sense.
        if line == "OBJSENSE"
            obj_sense = strip(readline(io))
            if obj_sense == "MIN"
                MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
            elseif obj_sense == "MAX"
                MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
            else
                error("Objective sense $obj_sense is not supported.")
            end
            continue
        end

        # Non-PSD variable constraints.
        if line == "VAR"
            raw_var_info = split(strip(readline(io)))
            @assert length(raw_var_info) == 2
            num_var = parse(Int, raw_var_info[1])
            num_lines = parse(Int, raw_var_info[2])
            append!(scalar_vars, MOI.add_variables(model, num_var))
            var_idx = 0
            for k in 1:num_lines
                raw_cone_info = split(strip(readline(io)))
                @assert length(raw_cone_info) == 2
                cone_str = raw_cone_info[1]
                cone_dim = parse(Int, raw_cone_info[2])
                if cone_str == "F" # Free cones (no constraint).
                    var_idx += cone_dim
                    continue
                end
                if cone_str in ("EXP", "EXP*") # Exponential cones.
                    # Reverse order of indices.
                    @assert cone_dim == 3
                    con_func = MOI.VectorOfVariables(scalar_vars[
                        [var_idx + 3, var_idx + 2, var_idx + 1]])
                    con_set = (cone_str == "EXP") ? MOI.ExponentialCone() :
                        MOI.DualExponentialCone()
                else
                    con_func = MOI.VectorOfVariables(scalar_vars[
                        (var_idx + 1):(var_idx + cone_dim)])
                    if startswith(cone_str, "@") # Power cones (parametric).
                        @assert cone_dim == 3
                        con_set = powcone_to_moi_cone(cone_str,
                            powcone_alphas, dpowcone_alphas)
                    else
                        con_set = cbf_to_moi_cone(cone_str, cone_dim)
                    end
                end
                MOI.add_constraint(model, con_func, con_set)
                var_idx += cone_dim
            end
            @assert var_idx == num_var
            continue
        end

        # Integrality constraints.
        if line == "INT"
            for k in 1:parse(Int, strip(readline(io)))
                var_idx = parse(Int, strip(readline(io))) + 1
                MOI.add_constraint(model, MOI.SingleVariable(
                    scalar_vars[var_idx]), MOI.Integer())
            end
            continue
        end

        # PSD variable constraints.
        if line == "PSDVAR"
            for k in 1:parse(Int, strip(readline(io)))
                side_dim = parse(Int, strip(readline(io)))
                cone_dim = div(side_dim * (side_dim + 1), 2)
                psd_vars_k = MOI.add_variables(model, cone_dim)
                push!(psd_vars, psd_vars_k)
                MOI.add_constraint(model, MOI.VectorOfVariables(psd_vars_k),
                    MOI.PositiveSemidefiniteConeTriangle(side_dim))
            end
            continue
        end

        # Objective function terms.
        if line == "OBJFCOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 4
                (psd_var_idx, i, j) = (parse(Int, raw_coord[i]) + 1
                    for i in 1:3)
                coef = parse(Float64, raw_coord[end])
                if i != j
                    coef += coef # scale off-diagonals
                end
                push!(obj_terms, MOI.ScalarAffineTerm{Float64}(coef,
                    psd_vars[psd_var_idx][mat_to_vec_idx(i, j)]))
            end
            continue
        end

        if line == "OBJACOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 2
                var_idx = parse(Int, raw_coord[1]) + 1
                coef = parse(Float64, raw_coord[end])
                push!(obj_terms, MOI.ScalarAffineTerm{Float64}(coef,
                    scalar_vars[var_idx]))
            end
            continue
        end

        if line == "OBJBCOORD"
            obj_constant += parse(Float64, strip(readline(io)))
            continue
        end

        # Non-PSD constraints.
        if line == "CON"
            raw_con_info = split(strip(readline(io)))
            @assert length(raw_con_info) == 2
            num_rows = parse(Int, raw_con_info[1])
            num_lines = parse(Int, raw_con_info[2])
            row_idx = 0
            for k in 1:num_lines
                raw_cone_info = split(strip(readline(io)))
                @assert length(raw_cone_info) == 2
                cone_str = raw_cone_info[1]
                cone_dim = parse(Int, raw_cone_info[2])
                push!(con_cones, (cone_str, cone_dim))
                row_idx += cone_dim
            end
            @assert row_idx == num_rows
            append!(row_terms, Vector{MOI.ScalarAffineTerm{Float64}}()
                for k in 1:num_rows)
            append!(row_constants, zeros(num_rows))
            continue
        end

        if line == "FCOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 5
                (row_idx, psd_var_idx, i, j) = (parse(Int, raw_coord[i]) + 1
                    for i in 1:4)
                coef = parse(Float64, raw_coord[end])
                if i != j
                    coef += coef # scale off-diagonals
                end
                push!(row_terms[row_idx], MOI.ScalarAffineTerm{Float64}(
                    coef, psd_vars[psd_var_idx][mat_to_vec_idx(i, j)]))
            end
            continue
        end

        if line == "ACOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 3
                (row_idx, var_idx) = (parse(Int, raw_coord[i]) + 1
                    for i in 1:2)
                coef = parse(Float64, raw_coord[end])
                push!(row_terms[row_idx], MOI.ScalarAffineTerm{Float64}(
                    coef, scalar_vars[var_idx]))
            end
            continue
        end

        if line == "BCOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 2
                row_idx = parse(Int, raw_coord[1]) + 1
                row_constants[row_idx] = parse(Float64, raw_coord[end])
            end
            continue
        end

        # PSD constraints.
        if line == "PSDCON"
            idx = 0
            for k in 1:parse(Int, strip(readline(io)))
                side_dim = parse(Int, strip(readline(io)))
                push!(psd_side_dims, side_dim)
                push!(psd_row_starts, idx)
                idx += div(side_dim * (side_dim + 1), 2)
            end
            append!(psd_row_terms, Vector{MOI.ScalarAffineTerm{Float64}}()
                for i in 1:idx)
            append!(psd_row_constants, zeros(idx))
            continue
        end

        if line == "HCOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 5
                (psd_idx, var_idx, i, j) = (parse(Int, raw_coord[i]) + 1
                    for i in 1:4)
                coef = parse(Float64, raw_coord[end])
                row_idx = psd_row_starts[psd_idx] + mat_to_vec_idx(i, j)
                push!(psd_row_terms[row_idx], MOI.ScalarAffineTerm{Float64}(
                    coef, scalar_vars[var_idx]))
            end
            continue
        end

        if line == "DCOORD"
            for k in 1:parse(Int, strip(readline(io)))
                raw_coord = split(strip(readline(io)))
                @assert length(raw_coord) == 4
                (psd_idx, i, j) = (parse(Int, raw_coord[i]) + 1 for i in 1:3)
                row_idx = psd_row_starts[psd_idx] + mat_to_vec_idx(i, j)
                psd_row_constants[row_idx] += parse(Float64, raw_coord[end])
            end
            continue
        end

        error("Failed to parse CBF file due to corrupted line: $line")
    end

    # Objective function.
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(obj_terms, obj_constant))

    # Non-PSD constraints.
    row_idx = 0
    for (cone_str, cone_dim) in con_cones
        if cone_str in ("EXP", "EXP*") # Exponential cones.
            # Reverse order of indices.
            @assert cone_dim == 3
            con_func = MOI.VectorAffineFunction([MOI.VectorAffineTerm{Float64}(
                4 - l, t) for l in 1:cone_dim for t in row_terms[row_idx + l]],
                row_constants[[row_idx + 3, row_idx + 2, row_idx + 1]])
            con_set = (cone_str == "EXP") ? MOI.ExponentialCone() :
                MOI.DualExponentialCone()
        else
            con_func = MOI.VectorAffineFunction([MOI.VectorAffineTerm{Float64}(
                l, t) for l in 1:cone_dim for t in row_terms[row_idx + l]],
                row_constants[(row_idx + 1):(row_idx + cone_dim)])
            if startswith(cone_str, "@") # Power cones (parametric).
                @assert cone_dim == 3
                con_set = powcone_to_moi_cone(cone_str, powcone_alphas,
                    dpowcone_alphas)
            else
                con_set = cbf_to_moi_cone(cone_str, cone_dim)
            end
        end
        MOI.add_constraint(model, con_func, con_set)
        row_idx += cone_dim
    end

    # PSD constraints.
    for psd_idx in eachindex(psd_side_dims)
        row_start = psd_row_starts[psd_idx]
        side_dim = psd_side_dims[psd_idx]
        cone_dim = div(side_dim * (side_dim + 1), 2)
        con_func = MOI.VectorAffineFunction([MOI.VectorAffineTerm{Float64}(
            l, t) for l in 1:cone_dim for t in psd_row_terms[row_start + l]],
            psd_row_constants[(row_start + 1):(row_start + cone_dim)])
        MOI.add_constraint(model, con_func,
            MOI.PositiveSemidefiniteConeTriangle(side_dim))
    end

    return
end

end

ext: jl


text:
export getFaceAverageMatrix, getEdgeAverageMatrix, getNodalAverageMatrix
export ndgrid, meshgrid
export getDivergenceMatrix, getNodalGradientMatrix, getCurlMatrix, getNodalLaplacianMatrix,getNodalDivSigGradMatrix
export getEdgeMassMatrix, getdEdgeMassMatrix, getFaceMassMatrix,
       getdFaceMassMatrix,getNodalMassMatrix,getdNodalMassMatrix
export avN2C, avN2C_Nearest

## The matrix below properly averages Nodes to cells using 0.5,0.5
## However, we usually use this matrix the other way around: cells to nodes.
## Then, we need to decide what to do with the boundary. This matrix averages with 0.
## This matrix is used as default for all averging operators
function avN2C(n)
	# A = avN2C(n), 1D average operator
	I,J,V = SparseArrays.spdiagm_internal(0 => fill(.5,n), 1 => fill(.5,n))
	return sparse(I, J, V, n, n+1)
end

## The matrix below, uses a nearest neighbor interpolation on the boundaries, so that avN2C_Nearest'*ones = ones.
## This is important for geometric multigrid to work.
function avN2C_Nearest(n)
	# A = av(n), 1D average operator
	I,J,V = SparseArrays.spdiagm_internal(0 => fill(.5,n), 1 => fill(.5,n))
	av = sparse(I, J, V, n, n+1);
	av[1,1] = 1.0;
	av[end,end] = 1.0;
	return av;
end

function ddx(n)
# D = ddx(n), 1D derivative operator
	I,J,V = SparseArrays.spdiagm_internal(0 => fill(-1.0,n), 1 => fill(1.0,n))
	return sparse(I, J, V, n, n+1)
end

# --- linear operator
"""
 	function jInv.Mesh.getFaceAverageMatrix

	Returns Face-to-CellCenter average matrix from Mesh.Af.
	Matrix is constructed if Mesh.Af is empty.

	For 2D Mesh: Af = [A1 A2]
	For 3D Mesh: Af = [A1 A2 A3]

	Input:
		Mesh::Abstract Mesh

"""
function getFaceAverageMatrix(Mesh::AbstractTensorMesh; saveMat::Bool = true, avN2C::Function = avN2C)
	if isempty(Mesh.Af)
		if Mesh.dim==2
			A1 = kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),avN2C(Mesh.n[1]))
			A2 = kron(avN2C(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1]))
			Af = [A1 A2]
		elseif Mesh.dim==3
			I1 = sparse(1.0I, Mesh.n[1], Mesh.n[1]);
			I2 = sparse(1.0I, Mesh.n[2], Mesh.n[2]);
			I3 = sparse(1.0I, Mesh.n[3], Mesh.n[3]);
			A1 = kron(I3,kron(I2,avN2C(Mesh.n[1])))
			A2 = kron(I3,kron(avN2C(Mesh.n[2]),I1))
			A3 = kron(avN2C(Mesh.n[3]),kron(I2,I1))
			Af = [A1 A2 A3]
		end
		if saveMat
			Mesh.Af = Af;
		end
		return Af;
	else
		return Mesh.Af
	end
end

"""
 	function jInv.Mesh.getEdgeAverageMatrix

	Returns Edge-to-CellCenter average matrix from Mesh.Ae.
	Matrix is constructed if Mesh.Ae is empty.

	For 3D Mesh: Ae = [A1 A2 A3]

	Input:
		Mesh::Abstract Mesh

"""
function getEdgeAverageMatrix(Mesh::AbstractTensorMesh; saveMat::Bool = true, avN2C::Function = avN2C)
	if isempty(Mesh.Ae)
		if Mesh.dim==3
			Av1 = avN2C(Mesh.n[1]);
			Av2 = avN2C(Mesh.n[2]);
			Av3 = avN2C(Mesh.n[3]);
			A1 = kron(Av3,kron(Av2,speye(Mesh.n[1]) ))
			A2 = kron(Av3,kron(speye(Mesh.n[2]),Av1))
			A3 = kron(speye(Mesh.n[3]),kron(Av2,Av1))
			Ae = [A1 A2 A3]
		elseif Mesh.dim==2
			Av1 = avN2C(Mesh.n[1]);
			Av2 = avN2C(Mesh.n[2]);
			A1 = kron(Av2,sparse(1.0I, Mesh.n[1], Mesh.n[1]))
			A2 = kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),Av1)
			Ae = [A1 A2]
		else
			error("getEdgeAverageMatrix not implemented fot $(Mesh.dim)D Meshes")
		end
		if saveMat
			Mesh.Ae = Ae;
		end
		return Ae;
	else
		return Mesh.Ae
	end
end

"""
 	function jInv.Mesh.getNodalAverageMatrix

	Returns Nodal-to-CellCenter average matrix from Mesh.An.
	Matrix is constructed if Mesh.An is empty.

	Input:
		Mesh::Abstract Mesh

"""
function getNodalAverageMatrix(Mesh::AbstractTensorMesh; saveMat::Bool = true, avN2C::Function = avN2C)
# Mesh.An = getNodalAverageMatrix(Mesh::TensorMesh3D) builds nodal-to-cc average operator
	if isempty(Mesh.An)
		if Mesh.dim==2
			An = kron(avN2C(Mesh.n[2]),avN2C(Mesh.n[1]))
		elseif Mesh.dim==3
			An = kron(avN2C(Mesh.n[3]),kron(avN2C(Mesh.n[2]),avN2C(Mesh.n[1])))
		else
			error("getNodalAverageMatrix: Dimension must be 2 or 3")
		end
		if saveMat
			Mesh.An = An;
		end
		return An;
	else
		return Mesh.An
	end
end

"""
	function jInv.Mesh.getEdgeMassMatrix(mesh,sigma)

	Returns mass matrix on cell edges, weighted by vector sigma.
	Matrix is always constructed. Uses pre-constructed edge averaging
	and cell volume matrices if available.

	Input:
		Mesh::Abstract Mesh
	       sigma::Vector

	Output:
		SparseMatrixCSC{Float64,Int64}
"""
function getEdgeMassMatrix(M::AbstractMesh,sigma::Vector)

	Ae   = getEdgeAverageMatrix(M)
	V    = getVolume(M)
	Masse    = sdiag(Ae'*(V*vec(sigma)))
	return Masse
end

"""
	function jInv.Mesh.getdEdgeMassMatrix(M,v)

	Returns directional derivative of edge mass matrix w.r.t. sigma, i.e.,

			d_sigma (M(sigma)*v)

	Matrix is always constructed. Uses pre-constructed edge averaging
	and cell volume matrices if available.

	Input:
		Mesh::Abstract - Mesh
		   v::Vector   - edge vector defining derivative

	Output:
		SparseMatrixCSC{Float64,Int64}
"""
function getdEdgeMassMatrix(M::AbstractMesh,v::Vector)

	Ae   = getEdgeAverageMatrix(M)
	V    = getVolume(M)
	return sdiag(v)*Ae'*V
end

# for compatibility with applications supporting isotropic/anisotropic PDE coefficients
getdEdgeMassMatrix(M::AbstractMesh,sigma::Vector,v::Vector) = getdEdgeMassMatrix(M, v)

"""
	function jInv.Mesh.getFaceMassMatrix(M,sigma)

	Returns face mass matrix, weighted by vector sigma.
	Matrix is always constructed. Uses pre-constructed face averaging
	and cell volume matrices if available.

	Input:
		M::Abstract Mesh
	       sigma::Vector

	Output:
		SparseMatrixCSC{Float64,Int64}
"""
function getFaceMassMatrix(M::AbstractMesh,sigma::Vector; saveMat::Bool = true, avN2C::Function = avN2C)
	Af    = getFaceAverageMatrix(M;saveMat=saveMat,avN2C=avN2C)
	V     = getVolume(M;saveMat=saveMat)
	Massf = sdiag(Af'*(V*sigma));
	return Massf
end

"""
	function jInv.Mesh.getdFaceMassMatrix(M,v)

	Returns directional derivative of face mass matrix w.r.t sigma, i.e.

		d_sigma (M(sigma)*v)

	Matrix is always constructed. Uses pre-constructed face averaging
	and cell volume matrices if available.

	Input:
		M::Abstract  - Mesh
		   v::Vector -  face vector defining derivative

	Output:
		SparseMatrixCSC{Float64,Int64}
"""
function getdFaceMassMatrix(M::AbstractMesh,v::Vector)
	Af    = getFaceAverageMatrix(M)
	V     = getVolume(M)
	Massf = sdiag(v)*Af'*V
	return Massf
end

# for compatibility with applications supporting isotropic/anisotropic PDE coefficients
getdFaceMassMatrix(M::AbstractMesh,sigma::Vector,v::Vector) = getdFaceMassMatrix(M, v)

"""
	function jInv.Mesh.getNodalMassMatrix(M,sigma)

	Returns nodal mass matrix, weighted by vector sigma.
	Matrix is always constructed. Uses pre-constructed nodal averaging
	and cell volume matrices if available.

	Input:
		M::Abstract Mesh
	       sigma::Vector

	Output:
		SparseMatrixCSC{Float64,Int64}
"""
function getNodalMassMatrix(M::AbstractMesh,sigma::Vector)
    An = getNodalAverageMatrix(M)
    V  = getVolume(M)
    Mn = Diagonal(An'*(V*sigma))
  return Mn
end

"""
	function jInv.Mesh.getdNodalMassMatrix(M,v)

	Returns directional derivative of nodal mass matrix w.r.t. sigma, i.e.,

		d_sigma (M(sigma)*v)

	Matrix is always constructed. Uses pre-constructed nodal averaging
	and cell volume matrices if available.

	Input:
		M::Abstract  - Mesh
		   v::Vector - nodal vector defining derivative

	Output:
		SparseMatrixCSC{Float64,Int64}
"""
function getdNodalMassMatrix(M::AbstractMesh,v::Vector)
	An  = getNodalAverageMatrix(M)
	V   = getVolume(M)
	dMn = Diagonal(v)*An'*V
end

# for compatibility with applications supporting isotropic/anisotropic PDE coefficients
getdNodalMassMatrix(M::AbstractMesh,sigma::Vector,v::Vector) = getdNodalMassMatrix(M, v)

# --- ndgrid
ndgrid(v::AbstractVector) = copy(v)
# function ndgrid{T}(v1::AbstractVector{T}, v2::AbstractVector{T})
# 	m, n = length(v1), length(v2)
# 	v1 = reshape(v1, m, 1)
# 	v2 = reshape(v2, 1, n)
# 	(repeat(v1, 1, n), repeat(v2, m, 1))
# end

function ndgrid_fill(a, v, s, snext)
	for j = 1:length(a)
		a[j] = v[div(rem(j-1, snext), s)+1]
	end
end

function ndgrid(vs::AbstractVector{T}...) where {T}
	n = length(vs)
	sz = map(length, vs)
	out = ntuple(i->Array{T}(undef,sz), n)
	s = 1
	for i=1:n
		a = out[i]::Array
		v = vs[i]
		snext = s*size(a,i)
		ndgrid_fill(a, v, s, snext)
		s = snext
	end
	out
end

# --- meshgrid
meshgrid(v::AbstractVector) = meshgrid(v, v)
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
	m, n = length(vy), length(vx)
	vx = reshape(vx, 1, n)
	vy = reshape(vy, m, 1)
	(repeat(vx, m, 1), repeat(vy, 1, n))
end

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T},
	                 vz::AbstractVector{T}) where {T}
	m, n, o = length(vy), length(vx), length(vz)
	vx = reshape(vx, 1, n, 1)
	vy = reshape(vy, m, 1, 1)
	vz = reshape(vz, 1, 1, o)
	om = ones(Int, m)
	on = ones(Int, n)
	oo = ones(Int, o)
	(vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end

# --- differential operators
function getDivergenceMatrix(Mesh::AbstractTensorMesh;saveMat = true)
# Mesh.Div = getDivergenceMatrix(Mesh::AbstractTensorMesh) builds face-to-cc divergence operator
	if isempty(Mesh.Div)
		if Mesh.dim==2
			D1 = kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),ddx(Mesh.n[1]))
			D2 = kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1]))
			Div = [D1 D2]
		elseif Mesh.dim==3
			D1 = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),ddx(Mesh.n[1])))
			D2 = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
			D3 = kron(ddx(Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
			Div = [D1 D2 D3]
		end
		Vi = getVolumeInv(Mesh)
		F  = getFaceArea(Mesh)
		Div = Vi*(Div*F);
		if saveMat
			Mesh.Div = Div;
		end
		return Div;
	else
		return Mesh.Div;
	end
end

function getNodalGradientMatrix(Mesh::AbstractTensorMesh; saveMat::Bool = true)
# Mesh.Grad = getNodalGradientMatrix(Mesh::AbstractTensorMesh) builds nodal-to-edge gradient operator
	if isempty(Mesh.Grad)
		if Mesh.dim==2
			G1 = kron(sparse(1.0I, Mesh.n[2]+1, Mesh.n[2]+1),ddx(Mesh.n[1]))
			G2 = kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1]+1, Mesh.n[1]+1))
			Grad =[G1; G2]
		elseif Mesh.dim==3
			G1 = kron(sparse(1.0I, Mesh.n[3]+1, Mesh.n[3]+1),kron(sparse(1.0I, Mesh.n[2]+1, Mesh.n[2]+1),ddx(Mesh.n[1])))
			G2 = kron(sparse(1.0I, Mesh.n[3]+1, Mesh.n[3]+1),kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1]+1, Mesh.n[1]+1)))
			G3 = kron(ddx(Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2]+1, Mesh.n[2]+1),sparse(1.0I, Mesh.n[1]+1, Mesh.n[1]+1)))
			Grad =[G1; G2; G3]
		end
		Li  = getLengthInv(Mesh)
		Grad = Li*Grad;
		if saveMat
			Mesh.Grad = Grad
		end
		return Grad;
	else
		return Mesh.Grad
	end

end

function getCurlMatrix(Mesh::AbstractTensorMesh)
# Mesh.Curl = getCurlMatrix(Mesh::AbstractTensorMesh) builds edge-to-face curl operator
	if isempty(Mesh.Curl)
		if Mesh.dim==3
			# The Curl from edges to faces
			Dyz = kron(ddx(Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),sparse(1.0I, Mesh.n[1]+1, Mesh.n[1]+1)))
			Dzy = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1]+1, Mesh.n[1]+1)))

			Dxz = kron(ddx(Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2]+1, Mesh.n[2]+1),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
			Dzx = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2]+1, Mesh.n[2]+1),ddx(Mesh.n[1])))

			Dxy = kron(sparse(1.0I, Mesh.n[3]+1, Mesh.n[3]+1),kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
			Dyx = kron(sparse(1.0I, Mesh.n[3]+1, Mesh.n[3]+1),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),ddx(Mesh.n[1])))

			# curl on the edges
			Curl = [
				 spzeros(Mesh.nf[1],Mesh.ne[1])  -Dyz   Dzy
				 Dxz   spzeros(Mesh.nf[2],Mesh.ne[2])  -Dzx
				-Dxy   Dyx   spzeros(Mesh.nf[3],Mesh.ne[3])]

				Fi = getFaceAreaInv(Mesh)
				L  = getLength(Mesh)
				Mesh.Curl = Fi*(Curl*L)
		else
			error("CURL is only implemented for 3D meshes.")
		end

	end
	return Mesh.Curl
end

function getNodalLaplacianMatrix(Mesh::AbstractTensorMesh; saveMat::Bool = true)

	if isempty(Mesh.nLap)
		G = getNodalGradientMatrix(Mesh)
		L = G'*G;
		if saveMat
			Mesh.nLap = L;
		end
		return L;
	else
		return Mesh.nLap;
	end

end

function getNodalDivSigGradMatrix(M::AbstractTensorMesh,sig::Vector; saveMat::Bool = true,avN2C::Function = avN2C)
 	G       = getNodalGradientMatrix(M, saveMat = saveMat)
	Ae      = getEdgeAverageMatrix(M, saveMat = saveMat,avN2C = avN2C);
    L       = (G'*sparse(Diagonal(Ae'*(sig))))*G
	return L
end


text:
function add(g::GS, o, s::Symbol)
    v = getfield(o, s)
    val2str(s, v) |> g
    return
end

function add(g::GS, o, ss...)
    for s in ss
        add(g, o, s)
    end
end

val2str(s, x)         = "$s $(val2str(x))"
val2str(s, ::Nothing) = ""
val2str(s, b::Bool)   = ifelse(b, "$s", "")

val2str(x::Color)  = "\"$(col2str(x))\""
val2str(x::F64)    = fl2str(x)
val2str(x::Int)    = string(x)
val2str(x::String) = x


text:
export hadamard, sgate, sdggate, tgate, tdggate, phaseshift, singleunitary,
       controlsgate, cnot, toffoli, swap

"""
---
       sigmax()

```math
\\sigma_x = \\begin{bmatrix}
       0 & 1 \\\\
       1 & 0
\\end{bmatrix}
```
"""
sigmax() = sigmax(SpinBasis(1//2))

"""
---
       sigmay()

```math
\\sigma_y = \\begin{bmatrix}
       0 & -i \\\\
       i & 0
\\end{bmatrix}
```
"""
sigmay() = sigmay(SpinBasis(1//2))

"""
---
       sigmaz()

```math
\\sigma_z = \\begin{bmatrix}
       1 & 0 \\\\
       0 & -1
\\end{bmatrix}
```
"""
sigmaz() = sigmaz(SpinBasis(1//2))

"""
---
       sigmam()

```math
\\sigma_- = | 1 \\rangle \\langle 0 | = \\begin{bmatrix}
       0 & 0 \\\\
       1 & 0
\\end{bmatrix}
```
"""
sigmam() = sigmam(SpinBasis(1//2))

"""
---
       sigmap()

```math
\\sigma_+ = | 0 \\rangle \\langle 1 | = \\begin{bmatrix}
       0 & 1 \\\\
       0 & 0
\\end{bmatrix}
```
"""
sigmap() = sigmap(SpinBasis(1//2))

"""
       hadamard(n::Int=1)

n tensor product Hadamard gate.
```math
H^{\\otimes n} = \\frac{1}{\\sqrt{2^n}} \\begin{bmatrix}
       1 & 1 \\\\
       1 & -1
\\end{bmatrix}^{\\otimes n}
```
"""
hadamard(n::Int=1) = tensor([DenseOperator(SpinBasis(1//2), 1/sqrt(2)*[1 1; 1 -1]) for i in 1:n]...)

"""
       sgate()

```math
S = \\sqrt{\\sigma_z} = \\begin{bmatrix}
       1 & 0 \\\\
       0 & i
\\end{bmatrix}
```
"""
function sgate()
       x = spzeros(Complex{Float64}, 2, 2)
       x[1,1] = 1.
       x[2,2] = im
       return SparseOperator(SpinBasis(1//2), x)
end

"""
       sdggate()

```math
S^\\dagger = (\\sqrt{\\sigma_z})^\\dagger = \\begin{bmatrix}
       1 & 0 \\\\
       0 & -i
\\end{bmatrix}
```
"""
function sdggate()
       x = spzeros(Complex{Float64}, 2, 2)
       x[1,1] = 1.
       x[2,2] = -im
       return SparseOperator(SpinBasis(1//2), x)
end

"""
       tgate()

```math
T = \\begin{bmatrix}
       1 & 0 \\\\
       0 & \\exp(i \\pi/4)
\\end{bmatrix}
```
"""
function tgate()
       x = spzeros(Complex{Float64}, 2, 2)
       x[1,1] = 1.
       x[2,2] = (1 + im) / sqrt(2)
       return SparseOperator(SpinBasis(1//2), x)
end

"""
       tdggate()

```math
T^\\dagger = \\begin{bmatrix}
       1 & 0 \\\\
       0 & \\exp(- i \\pi/4)
\\end{bmatrix}
```
"""
function tdggate()
       x = spzeros(Complex{Float64}, 2, 2)
       x[1,1] = 1.
       x[2,2] = (1 - im) / sqrt(2)
       return SparseOperator(SpinBasis(1//2), x)
end

"""
       phaseshift(λ)

Phase shift gate. `λ` must be real number.
This is equivalent to `singleunitary(0, 0, λ)`.
```math
U(\\lambda) = \\begin{bmatrix}
       1 & 0 \\\\
       0 & e^{i \\lambda}
\\end{bmatrix}
```
"""
function phaseshift(λ::Real)
       x = spzeros(Complex{Float64}, 2, 2)
       x[1,1] = 1.
       x[2,2] = exp(im * λ)
       return SparseOperator(SpinBasis(1//2), x)
end

"""
       singleunitary(θ, ϕ, λ)

Single qubit unitary gate.
```math
U(\\theta, \\phi, \\lambda) = \\begin{bmatrix}
       \\cos(\\theta/2) & -e^{i \\lambda}\\sin(\\theta/2) \\\\
       e^{i \\phi}\\sin(\\theta/2) & e^{i (\\lambda + \\phi)} \\cos(\\theta/2)
\\end{bmatrix}
```
"""
function singleunitary(θ, ϕ, λ)
       x = spzeros(Complex{Float64}, 2, 2)
       x[1,1] = cos(θ/2)
       x[1,2] = - exp(im * λ) * sin(θ/2)
       x[2,1] = exp(im * ϕ) * sin(θ/2)
       x[2,2] = exp(im * (λ+ϕ)) * cos(θ/2)
       return SparseOperator(SpinBasis(1//2), x)
end

"""
       controlsgate(N, controls, target, op)

Multiple controlled gate.
![controlled gate](https://raw.githubusercontent.com/goropikari/QuantumInformation.jl/master/docs/pictures/controls.png)
# Arguments
- `N::Int`: the number of qubits.
- `controls::Vector{Int}`: list of controlled qubit positions.
- `target::Int`: position of target qubit.
- `op::Operator`: Operator that acts on target qubit.

# Example
```jldoctest
julia> nqubits, controls, target, op = 3, [1, 2], 3, sigmax();

julia> controlsgate(nqubits, controls, target, op) # equivalent to toffoli gate
SparseOperator(dim=8x8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
  [1, 1]  =  1.0+0.0im
  [2, 2]  =  1.0+0.0im
  [3, 3]  =  1.0+0.0im
  [4, 4]  =  1.0+0.0im
  [5, 5]  =  1.0+0.0im
  [6, 6]  =  1.0+0.0im
  [8, 7]  =  1.0+0.0im
  [7, 8]  =  1.0+0.0im
```
"""
function controlsgate(N::Int, controls::Vector{Int}, target::Int, op::QuantumOptics.operators.Operator)
       nc = length(controls) # number of control qubits

       N < nc + 1 && throw("the number of qubits must be larger or equal to $(maximum(vcat(controls, target)))")
       isempty(controls) && error("you don\'t specify control qubits.")
       target ∈ controls && error("target qubit overlaps a control qubit.")
       (target < 1 || target > N) && error("target qubit is out of bound.")

       if !issorted(controls)
              sort!(controls)
       end
       basis = SpinBasis(1//2)
       ctrlstate = sparse_spinalldown_dm(nc) # |11...1><11...1|
       if N == nc + 1
              # (I - |1...1><1...1|) ⊗ I + |1...1><1...1| ⊗ op
              basicstate = (identityoperator(basis^nc) - ctrlstate) ⊗ identityoperator(basis) + ctrlstate ⊗ op
              perm = collect(1:N)
              perm[end], perm[target] = perm[target], perm[end]
       else
              # I ⊗ ⋯ ⊗ I ⊗ ( (I - |1...1><1...1|) ⊗ I + |1...1><1...1| ⊗ op )
              basicstate = tensor([identityoperator(basis) for i in 1:N-(nc+1)]...) ⊗ ( (identityoperator(basis^nc) - ctrlstate) ⊗ identityoperator(basis) + ctrlstate ⊗ op )
              perm = collect(1:N)
              operator_position = N
              for (idx, cidx) in enumerate(controls)
                    perm[cidx], perm[N - nc - 1 + idx] = perm[N - nc - 1 + idx], perm[cidx]
                    if cidx == operator_position
                           operator_position = N - nc - 1 + idx
                    end
              end
              perm[target], perm[operator_position] = perm[operator_position], perm[target]
       end

       return permutesystems(basicstate, perm)
end

"""
       cnot(N::Int=2, control::Int=1, target::Int=2)

# Arguments
- `N::Int`: the number of qubits.
- `control::Int`: position of control qubit.
- `target::Int`: position of target qubit.

# Example
```julia
julia> cnot()
SparseOperator(dim=4x4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
  [1, 1]  =  1.0+0.0im
  [2, 2]  =  1.0+0.0im
  [4, 3]  =  1.0+0.0im
  [3, 4]  =  1.0+0.0im

julia> nqubits, control, target = 3, 1, 3;

julia> cnot(nqubits, control, target)
SparseOperator(dim=8x8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
  [1, 1]  =  1.0+0.0im
  [2, 2]  =  1.0+0.0im
  [3, 3]  =  1.0+0.0im
  [4, 4]  =  1.0+0.0im
  [6, 5]  =  1.0+0.0im
  [5, 6]  =  1.0+0.0im
  [8, 7]  =  1.0+0.0im
  [7, 8]  =  1.0+0.0im
```
"""
cnot(N::Int=2, control::Int=1, target::Int=2) = controlsgate(N, [control], target, sigmax())

"""
       toffoli(N::Int=3, controls::Vector{Int}=[1, 2], target::Int=3)

Toffoli gate.
![toffoli gate](https://raw.githubusercontent.com/goropikari/QuantumInformation.jl/master/docs/pictures/toffoli.png)

# Example
```julia
julia> toffoli() # equivalent to  toffoli(3, [1,2], 3)
SparseOperator(dim=8x8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
  [1, 1]  =  1.0+0.0im
  [2, 2]  =  1.0+0.0im
  [3, 3]  =  1.0+0.0im
  [4, 4]  =  1.0+0.0im
  [5, 5]  =  1.0+0.0im
  [6, 6]  =  1.0+0.0im
  [8, 7]  =  1.0+0.0im
  [7, 8]  =  1.0+0.0im

julia> nqubits, controls, target = 4, [1,2], 3;

julia> toffoli(nqubits, controls, target)
SparseOperator(dim=16x16)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
  [1 ,  1]  =  1.0+0.0im
  [2 ,  2]  =  1.0+0.0im
  [3 ,  3]  =  1.0+0.0im
  [4 ,  4]  =  1.0+0.0im
  [5 ,  5]  =  1.0+0.0im
  [6 ,  6]  =  1.0+0.0im
  [7 ,  7]  =  1.0+0.0im
  [8 ,  8]  =  1.0+0.0im
  [9 ,  9]  =  1.0+0.0im
  [10, 10]  =  1.0+0.0im
  [11, 11]  =  1.0+0.0im
  [12, 12]  =  1.0+0.0im
  [15, 13]  =  1.0+0.0im
  [16, 14]  =  1.0+0.0im
  [13, 15]  =  1.0+0.0im
  [14, 16]  =  1.0+0.0im
```
"""
function toffoli(N::Int=3, controls::Vector{Int}=[1, 2], target::Int=3)
       if N < 3
              error("the number of qubits must be larger or equal to 3.")
       end
       return controlsgate(N, controls, target, sigmax())
end

"""
       swap(N::Int=2, targets::Vector{Int}=[1,2])

# Exapmle
```julia
julia> swap()
SparseOperator(dim=4x4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
  [1, 1]  =  1.0+0.0im
  [3, 2]  =  1.0+0.0im
  [2, 3]  =  1.0+0.0im
  [4, 4]  =  1.0+0.0im

julia> swap(3, [1,2])
SparseOperator(dim=8x8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
  [1, 1]  =  1.0+0.0im
  [2, 2]  =  1.0+0.0im
  [5, 3]  =  1.0+0.0im
  [6, 4]  =  1.0+0.0im
  [3, 5]  =  1.0+0.0im
  [4, 6]  =  1.0+0.0im
  [7, 7]  =  1.0+0.0im
  [8, 8]  =  1.0+0.0im
```
"""
function swap(N::Int=2, targets::Vector{Int}=[1,2])
       if N < 2
              error("the number of qubits must be larger or equal to 2.")
       end
       if targets[1] == targets[2]
              error("target qubits mustn't be same.")
       end
       if length(targets) != 2
              error("the number of target qubits must be equal to 2")
       end
       x = spzeros(Complex{Float64}, 4, 4)
       x[1,1] = 1.
       x[3,2] = 1.
       x[2,3] = 1.
       x[4,4] = 1.
       basicswap = SparseOperator(SpinBasis(1//2)^2, x)

       if N == 2
              return basicswap
       else
              perm = collect(1:N)
              state = identityoperator(SpinBasis(1//2)^(N-2)) ⊗ basicswap
              target1, target2 = ifelse(targets[1] < targets[2], (targets[1], targets[2]), (targets[2], targets[1]))
              perm[target1], perm[end-1] = perm[end-1], perm[target1]
              perm[target2], perm[end] = perm[end], perm[target2]

              return permutesystems(state, perm)
       end
end


text:
module Game

import Util
import Moves
import Search

function pickMoveUCI(board,colour,history)
    Search.search(3,colour,history,board)
    return Search.best_move
end

function play(board,colour,history)
    states = []
    while(true)
        opp_colour = (colour === 'W') ? 'B' : 'W'
        board,states = @time takeTurn(board,history,colour,states)
        colour = opp_colour
        if(checkEndGame(board,states,history,colour)) break end
    end

end

function takeTurn(board,history,colour,states)
    move = pickMoveUCI(board,colour,history)
    push!(history,move)
    println("move picked: ", move)
    board = Moves.makeMoveUCI(move,board,colour)
    state = reduce((x,y) -> x | y, values(board))
    push!(states,state)
    Util.prettyPrintBoard(board)
    return board,states
end

# //=============================  END GAME LOGIC  =============================//
# //============================================================================//

function checkEndGame(board,states,history,colour)
    opp_colour = (colour ==='W') ? 'B' : 'W'
    if (checkForMate(board,colour,history) === 1 || checkForMate(board,colour,history) === 2 ||
        sameBoardStateFiveTimes(states) || fiftyMoves(history) || ~(hasEnoughMaterial(board,colour) || hasEnoughMaterial(board,opp_colour)))
        println("GAME OVER")
        print(board)
        return true
    end
    return false
end

function checkForMate(board,colour,history)
    check_mate = Moves.isCheckMate(board,colour,history)
    if(check_mate == true)
        println("checkmate")
    end
    return check_mate
end

function sameBoardStateFiveTimes(states)
    groups = Dict()
    for i in states
        if i in keys(groups)
            groups[i] = groups[i] + 1
            if(groups[i] > 5)
                println("same state five times")
                return true
            end
        else
            groups[i] = 1
        end
    end
    return false
end

# The fifty-move rule in chess states that a player can claim a draw if no capture has been made and no pawn has been moved in the last fifty moves
# (for this purpose a "move" consists of a player completing their turn followed by the opponent completing their turn)
function fiftyMoves(history)
    if(length(history) < 50)
        return false
    end
    fifty = history[length(history)-49:length(history)]
    pawn_moves = [] #how to get pawn moves from UCI notated?
    if(length(pawn_moves) > 0)
        return false
    end
    captures = ['l'] #how to get caps
    if(length(captures) > 0)
        return false
    end
    print("fifty moves")
    return true
end

function hasEnoughMaterial(board,colour)
    if(board[colour * 'P'] > 0 || (board[colour * 'N'] > 0 && board[colour * 'B'] > 0) || board[colour * 'R'] > 0
        || board[colour * 'Q'] > 0  || count_ones(board[colour * 'B']) === 2)
        return true
    end
    print("not enough material")
    return false
end

end

max_stars_repo_name: grahamclyne/chessEngine

max_issues_repo_name: grahamclyne/chessEngine

max_forks_repo_name: grahamclyne/chessEngine


text:
using Checkpointing
using Enzyme
using ReverseDiff
using Test

function f(x)
    return [(x[1]*(x[2]-2.0))^2, (x[1]-1.0)^2*x[2]^2]
end

function f2(x,res)
    y = f(x)
    copyto!(res,y)
    return nothing
end

x = [2.0,6.0]

J = ReverseDiff.jacobian(f, x)
@show J

dx = copy(x)
dx = [0.0,0.0]
y = [0.0,0.0]
dy = [0.0,1.0]
autodiff(f2, Duplicated(x,dx), Duplicated(y, dy))

J2 = similar(J)
fill!(J2, 0)
for i in 1:2
    x = [2.0,6.0]
    fill!(dx, 0)
    fill!(y, 0)
    dy[i] = 1.0
    autodiff(f2, Duplicated(x,dx), Duplicated(y, dy))
    J2[i,:] = dx[:]
end

# Is correct, but gives the memmove warning
@test J2 ≈ J

include("../examples/optcontrol.jl")

@testset "Testing Enzyme in Revolve..." begin
    global steps = 100
    global snaps = 3
    global info = 1

    function store(F_H, F_C,t, i)
        F_C[1,i] = F_H[1]
        F_C[2,i] = F_H[2]
        F_C[3,i] = t
        return
    end

    function restore(F_C, i)
        F_H = [F_C[1,i], F_C[2,i]]
        t = F_C[3,i]
        return F_H, t
    end
    revolve = Revolve(steps, snaps, store, restore; verbose=info)
    F_opt, F_final, L_opt, L = optcontrol(revolve, steps, ReverseDiffADTool())
    # Check whether ReverseDiff works
    @test isapprox(F_opt, F_final, rtol=1e-4)
    @test isapprox(L_opt, L, rtol=1e-4)
    revolve = Revolve(steps, snaps, store, restore; verbose=info)
    F_opt_enzyme, F_final_enzyme, L_opt_enzyme, L_enzyme = optcontrol(revolve, steps, EnzymeADTool())
    @test isapprox(F_final, F_final_enzyme)
    # Returns wrong adjoints
    @test_broken isapprox(L, L_enzyme, rtol=1e-4)
end


text:
include(joinpath(DATA_DIR,"data_5bus_pu.jl"))
include(joinpath(DATA_DIR,"data_14bus_pu.jl"))

checksys = false

@testset "Test System constructors from .jl files" begin
    tPowerSystem = System(nothing)

    for i in nodes5
          nodes5[i].angle = deg2rad(nodes5[i].angle)
    end

    sys5 = System(nodes5, thermal_generators5, loads5, nothing, nothing,  100.0, nothing,
                  nothing, nothing; runchecks = checksys)

    sys5b = System(nodes5, thermal_generators5, loads5, nothing, battery5,  100.0, nothing,
                   nothing, nothing; runchecks = checksys)

    # GitHub issue #234 - fix forecasts5 in data file, use new format
    #_sys5b = PowerSystems._System(nodes5, thermal_generators5, loads5, nothing, battery5,
    #                              100.0, forecasts5, nothing, nothing)
    #sys5b = System(_sys5b)

    sys5bh = System(nodes5, vcat(thermal_generators5, hydro_generators5), loads5, branches5,
                    battery5,  100.0, nothing, nothing, nothing; runchecks = checksys)

    sys5bh = System(; buses=nodes5,
                    generators=vcat(thermal_generators5, hydro_generators5),
                    loads=loads5,
                    branches=branches5,
                    storage=battery5,
                    basepower=100.0,
                    forecasts=nothing,
                    services=nothing,
                    annex=nothing,
                    runchecks = checksys)

    # Test Data for 14 Bus

    # GitHub issue #234 - fix forecasts5 in data file, use new format
    #_sys14 = PowerSystems._System(nodes14, thermal_generators14, loads14, nothing, nothing,
    #                            100.0, Dict{Symbol,Vector{<:Forecast}}(),nothing,nothing)
    #sys14 = System(_sys14)

    for i in nodes14
          nodes14[i].angle = deg2rad(nodes14[i].angle)
    end

    sys14b = PowerSystems.System(nodes14, thermal_generators14, loads14, nothing,
                                 battery14, 100.0, nothing, nothing, nothing; runchecks = checksys)
    sys14b = PowerSystems.System(nodes14, thermal_generators14, loads14, branches14,
                                 battery14, 100.0, nothing, nothing, nothing; runchecks = checksys)
end

@testset "Test System constructor from Matpower" begin
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "case5_re.m"))
end


text:
using ShapML
using Test
using DataFrames
using Random
using Statistics
using Distributed

@testset "shap" begin
  include("shap.jl")
end

max_stars_repo_name: nredell/shapML

max_issues_repo_name: nredell/shapML

max_forks_repo_name: nredell/shapML


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

mutable struct VirtualNetworkGatewayIPConfigurationPropertiesFormat <: SwaggerModel
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, IPAllocationMethod } # spec name: privateIPAllocationMethod
    subnet::Any # spec type: Union{ Nothing, SubResource } # spec name: subnet
    publicIPAddress::Any # spec type: Union{ Nothing, SubResource } # spec name: publicIPAddress
    privateIPAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddress
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function VirtualNetworkGatewayIPConfigurationPropertiesFormat(;privateIPAllocationMethod=nothing, subnet=nothing, publicIPAddress=nothing, privateIPAddress=nothing, provisioningState=nothing)
        o = new()
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("publicIPAddress"), publicIPAddress)
        setfield!(o, Symbol("publicIPAddress"), publicIPAddress)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("privateIPAddress"), privateIPAddress)
        setfield!(o, Symbol("privateIPAddress"), privateIPAddress)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type VirtualNetworkGatewayIPConfigurationPropertiesFormat

const _property_map_VirtualNetworkGatewayIPConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("subnet")=>Symbol("subnet"), Symbol("publicIPAddress")=>Symbol("publicIPAddress"), Symbol("privateIPAddress")=>Symbol("privateIPAddress"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_VirtualNetworkGatewayIPConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("privateIPAllocationMethod")=>"IPAllocationMethod", Symbol("subnet")=>"SubResource", Symbol("publicIPAddress")=>"SubResource", Symbol("privateIPAddress")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkGatewayIPConfigurationPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkGatewayIPConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayIPConfigurationPropertiesFormat[property_name]

function check_required(o::VirtualNetworkGatewayIPConfigurationPropertiesFormat)
    true
end

function validate_property(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }, name::Symbol, val)
end


text:
function Base.convert(::Type{HyperEdge}, e::ChemicalHyperEdge{T}; skipcatalysts=false) where {T}
    v = skipcatalysts ? unique(vertices(e)) : vertices(e)
    HyperEdge{T}(v)
end

function Base.convert(::Type{HyperGraph}, x::ChemicalHyperGraph; skipcatalysts=false)
    HyperGraph(convert.(HyperEdge, hyperedges(x), skipcatalysts=skipcatalysts))
end


text:
using Statistics
using Distributions
using OnlineStats
import Streamfall: run_node!

using Infiltrator

"""
Update base parameter values with a partial list of new values by index.

# Arguments
-----------
- base_params : the base parameter set to update
- target_idx : the index position of parameters to update
- updated_params : a flattened array of new parameter values in same order as given indices
"""
function update_partial(base_params, target_idx, updated_params)
    # Update base parameter values with given updated values
    n_selected_params = length(target_idx)
    tmp_idxs = collect(Iterators.partition(1:length(updated_params), n_selected_params))
    mod_params = []
    for idxs in tmp_idxs
        bin_params = updated_params[idxs]
        new_params = copy(base_params)
        for (idx, new_value) in zip(target_idx, bin_params)
            new_params[idx] = new_value
        end
        append!(mod_params, new_params)
    end

    return mod_params
end

"""Runs node online, determining state thresholds around quantiles.
"""
function run_online_node!(sn, v_id, climate, params, quantiles; releases=nothing)
    node = sn[v_id]

    # Get node parameters
    _, x0, __ = param_info(node; with_level=false)
    num_params = length(x0)

    timesteps = sim_length(climate)
    active_param_set = zeros(Int, timesteps)

    o = Quantile(quantiles, b=2000)
    param_idxs = nothing
    for ts in (1:timesteps)
        cmd = node.storage[ts]

        # Update param set based on CMD state
        # Quantiles
        fit!(o, cmd)
        thresholds = value(o)

        cmd_param_set_id, node_params, param_idxs = find_state_vars(cmd, thresholds, params, num_params, length(thresholds))
        if (ts == 1) || (cmd_param_set_id != active_param_set[ts-1])
            update_params!(node, node_params...)
        end

        # record timestep in which this param was active
        active_param_set[ts] = cmd_param_set_id

        Streamfall.run_node!(sn, v_id, climate, ts; extraction=releases)
    end

    return active_param_set, param_idxs
end

"""Generic run_node method that handles any single state variable.

Log-transforms the state value.
"""
function online_state_node!(sn, v_id, climate, state, params, quantiles; releases=nothing)
    node = sn[v_id]

    # Get node parameters
    _, x0, __ = param_info(node; with_level=false)
    num_params = length(x0)

    timesteps = sim_length(climate)
    active_param_set = zeros(Int, timesteps)

    if state != :rainfall
        state_var = getfield(node, state)
    else
        tmp_climate = subcatchment_data(node, climate)
        state_var = tmp_climate[:, "410730_P"]
    end

    o = Quantile(quantiles, b=2000)
    param_idxs = nothing
    for ts in (1:timesteps)
        # Update param set based on state
        tmp = state_var[ts]

        # Add constant to avoid log(0)
        state_val = log(tmp+1)

        # Update param set based on state
        fit!(o, state_val)
        thresholds = value(o)
        param_set_id, node_params, param_idxs = find_state_vars(state_val, thresholds, params, num_params, length(thresholds))
        if (ts == 1) || (param_set_id != active_param_set[ts-1])
            update_params!(node, node_params...)
        end

        # record timestep in which this param was active
        active_param_set[ts] = param_set_id

        Streamfall.run_node!(sn, v_id, climate, ts; extraction=releases)
    end

    return active_param_set, param_idxs
end

function state_obj_func(x, base_params, climate, sn, v_id, state, calib_data, metric, target_idx, offsets)
    param_idxs = nothing
    active_param_set = nothing
    extraction = nothing

    # Update base parameter values with given updated values
    mod_params = update_partial(base_params, target_idx, x)

    try
        active_param_set, param_idxs = online_state_node!(sn, v_id, climate, state, mod_params, offsets; releases=extraction)
    catch err
        if err isa AssertionError
            return 9999.0
        end

        throw(err)
    end

    node = sn[v_id]
    node_data = node.outflow
    h_data = calib_data[node.name]

    # Calculate score
    score = metric(active_param_set, param_idxs, h_data, node_data)

    # reset to clear stored values
    Streamfall.reset!(node)

    return score
end

"""Calibration with online statistics with thresholds based around specified standard deviation offsets.

Generalized to handle any one state variable.
"""
function state_based_calibrate(sn, v_id, state_name::Symbol, climate, calib_data, metric,
                                       target_idx, thresholds; kwargs...)

    # Set defaults as necessary
    defaults = (;
        MaxTime=CALIB_TIME,
        TraceInterval=300.0,
        PopulationSize=125
    )
    kwargs = merge(defaults, kwargs)

    node = sn[v_id]

    # Get node parameters
    _, x0, param_bounds = param_info(node; with_level=false)
    param_bounds = [param_bounds[i] for i in target_idx]  # subset bounds array

    # Create new optimization function
    opt_func = x -> state_obj_func(x, x0, climate, sn, v_id, state_name,
                                   calib_data, metric, target_idx, thresholds)

    # Set up parameters for each state threshold
    n_states = length(thresholds)  # lower, center, upper
    param_bounds = repeat(param_bounds, n_states)
    opt = bbsetup(opt_func; SearchRange=param_bounds,
                  kwargs...)

    res = bboptimize(opt)

    bs = best_candidate(res)
    @info "Calibrated $(v_id) ($(node.name)), with score: $(best_fitness(res))"
    @info "Best Params:" collect(bs)

    return res, opt
end

"""Generic run_node method that handles any single state variable.

Develops an "online" running statistic of states from the first N days,
which remain static for the rest of the analysis period.
"""
function online_burn_state_node!(sn, v_id, climate, state, params, quantiles; burn_in=1826, releases=nothing, log_transform=false)
    node = sn[v_id]

    # Get node parameters
    _, x0, __ = param_info(node; with_level=false)
    num_params = length(x0)

    timesteps = sim_length(climate)
    active_param_set = zeros(Int, timesteps)

    if state != :rainfall
        state_var = getfield(node, state)
    else
        # hard coded, I know :(
        state_var = Streamfall.subcatchment_data(node, climate)[:, "410730_P"]
        # TODO: Change over when ready
        # state_var = Streamfall.rainfall_data(node, climate)
    end

    o = Quantile(quantiles, b=2000)
    param_idxs = nothing
    thresholds = nothing
    for ts in (1:timesteps)
        # Update param set based on state
        state_val = state_var[ts]
        if log_transform
            # add constant to avoid log(0)
            tmp = log(state_val .+ 1e-3)
        else
            tmp = state_val
        end

        if ts < burn_in
            # Update thresholds based on state value so far
            fit!(o, tmp)
            thresholds = value(o)
        end

        param_set_id, node_params, param_idxs = find_state_vars(tmp, thresholds, params, num_params, length(thresholds))
        if (ts == 1) || (param_set_id != active_param_set[ts-1])
            update_params!(node, node_params...)
        end

        # record timestep in which this param was active
        active_param_set[ts] = param_set_id

        Streamfall.run_node!(sn, v_id, climate, ts; extraction=releases)
    end

    return active_param_set, param_idxs
end

function burn_state_obj_func(x, base_params, climate, sn, v_id, state, calib_data, metric, target_idx, offsets; burn_in=1826, log_transform=false)
    param_idxs = nothing
    active_param_set = nothing
    extraction = nothing

    # Update base parameter values with given updated values
    mod_params = update_partial(base_params, target_idx, x)

    try
        active_param_set, param_idxs = online_burn_state_node!(sn, v_id, climate, state, mod_params, offsets;
                                                               burn_in=burn_in, releases=extraction, log_transform=log_transform)
    catch err
        if err isa AssertionError
            return 9999.0
        end

        throw(err)
    end

    node = sn[v_id]
    node_data = node.outflow
    h_data = calib_data[node.name]

    # Calculate score
    score = metric(active_param_set, param_idxs, h_data, node_data)

    # reset to clear stored values
    Streamfall.reset!(sn)

    return score
end

"""Calibration with online statistics with thresholds based around specified standard deviation offsets.

Generalized to handle any one state variable.
"""
function burn_state_based_calibrate(sn, v_id, state_name::Symbol, climate, calib_data, metric,
                                       target_idx, thresholds; burn_in=1826, log_transform=false, kwargs...)

    # Set defaults as necessary
    defaults = (;
        MaxTime=CALIB_TIME,
        TraceInterval=300.0,
        PopulationSize=125
    )
    kwargs = merge(defaults, kwargs)

    node = sn[v_id]

    # Get node parameters
    _, x0, param_bounds = param_info(node; with_level=false)
    param_bounds = [param_bounds[i] for i in target_idx]  # subset bounds array

    # Create new optimization function
    opt_func = x -> burn_state_obj_func(x, x0, climate, sn, v_id, state_name,
                                        calib_data, metric, target_idx, thresholds; burn_in=burn_in, log_transform=log_transform)

    # Set up parameters for each state threshold
    n_states = length(thresholds)  # lower, center, upper
    param_bounds = repeat(param_bounds, n_states)
    opt = bbsetup(opt_func; SearchRange=param_bounds,
                  kwargs...)

    res = bboptimize(opt)

    bs = best_candidate(res)
    @info "Calibrated $(v_id) ($(node.name)), with score: $(best_fitness(res))"
    @info "Best Params:" collect(bs)

    return res, opt
end

"""
Meta-metric that applies the given metric (`metric`) to each threshold and returns
the combination of these with the metric applied to the entire time series
(using the provided `agg_func`).
"""
function bin_metric(active_param_set, param_idxs, obs, sim, metric, agg_func)
    n_len = length(sim)
    @assert length(obs) == n_len
    sub_scores = []
    for param_set in 1:length(param_idxs)
        out_set = findall(active_param_set .== param_set)
        obs_set = obs[out_set]
        sim_set = sim[out_set]
        if length(sim_set) > 0 && (length(obs_set) == length(sim_set))
            tmp = metric(obs_set, sim_set)
            if isnan(tmp)
                push!(sub_scores, 9999.9)
            else
                push!(sub_scores, tmp)
            end
        else
            push!(sub_scores, 9999.8)
        end
    end

    if isnothing(agg_func)
        return Tuple(sub_scores)
    end

    push!(sub_scores, metric(h_data, n_data))

    return agg_func(sub_scores)
end

"""
Meta-metric that applies the given metric (`metric`) to each threshold and returns
the combination of these with the metric applied to the entire time series
(using the provided `agg_func`).
"""
function weighted_bin_metric(active_param_set, param_idxs, h_data, n_data, metric, agg_func)
    n_len = length(n_data)
    @assert length(h_data) == n_len
    sub_scores = []
    for param_set in 1:length(param_idxs)
        out_set = findall(active_param_set .== param_set)
        obs = h_data[out_set]
        sim = n_data[out_set]
        n_sim = length(sim)
        push!(sub_scores, metric(obs, sim) * (n_sim / n_len))
    end

    push!(sub_scores, metric(h_data, n_data))

    return agg_func(sub_scores)
end

function targeted_threshold_metric(active_param_set, param_idxs, obs, sim, metric, target_threshold)
    out_set = findall(active_param_set .== target_threshold)
    obs_set = obs[out_set]
    sim_set = sim[out_set]

    if length(sim_set) == 0
        return 9999.9
    end

    tmp = metric(obs_set, sim_set)

    if isnan(tmp)
        @infiltrate
    end

    return tmp
end


text:
const DESCRIPTORS = joinpath(RTS_GMLC_DIR, "user_descriptors.yaml")
const MULTISTART_MAPPING = joinpath(RTS_GMLC_DIR, "generator_mapping_multi_start.yaml")

mutable struct TestDevice <: Device
    name::String
end

mutable struct TestRenDevice <: RenewableGen
    name::String
end

function create_rts_system(time_series_resolution = Dates.Hour(1))
    data = PowerSystemTableData(RTS_GMLC_DIR, 100.0, DESCRIPTORS)
    return System(data; time_series_resolution = time_series_resolution)
end

function create_rts_multistart_system(time_series_resolution = Dates.Hour(1))
    data = PowerSystemTableData(
        RTS_GMLC_DIR,
        100.0,
        DESCRIPTORS;
        generator_mapping_file = MULTISTART_MAPPING,
    )
    return System(data; time_series_resolution = time_series_resolution)
end

function create_rts_system_with_hybrid_system(; add_forecasts = true)
    sys = PSB.build_system(
        PSB.PSITestSystems,
        "test_RTS_GMLC_sys",
        add_forecasts = add_forecasts,
    )
    thermal_unit = first(get_components(ThermalStandard, sys))
    bus = get_bus(thermal_unit)
    electric_load = first(get_components(PowerLoad, sys))
    storage = first(get_components(GenericBattery, sys))
    renewable_unit = first(get_components(RenewableDispatch, sys))

    name = "Test H"
    h_sys = HybridSystem(
        name = name,
        available = true,
        status = true,
        bus = bus,
        active_power = 1.0,
        reactive_power = 1.0,
        thermal_unit = thermal_unit,
        electric_load = electric_load,
        storage = storage,
        renewable_unit = renewable_unit,
        base_power = 100.0,
        operation_cost = TwoPartCost(nothing),
    )
    add_component!(sys, h_sys)
    return sys
end

"""Allows comparison of structs that were created from different parsers which causes them
to have different UUIDs."""
function compare_values_without_uuids(x::T, y::T) where {T <: IS.InfrastructureSystemsType}
    match = true

    for (fieldname, fieldtype) in zip(fieldnames(T), fieldtypes(T))
        if fieldname === :internal
            continue
        end

        val1 = getfield(x, fieldname)
        val2 = getfield(y, fieldname)

        # Recurse if this is an InfrastructureSystemsType
        if val1 isa IS.InfrastructureSystemsType
            if !compare_values_without_uuids(val1, val2)
                match = false
            end
            continue
        end

        if val1 != val2
            @error "values do not match" fieldname repr(val1) repr(val2)
            match = false
        end
    end

    return match
end

"""Return the first component of type component_type that matches the name of other."""
function get_component_by_name(sys::System, component_type, other::Component)
    for component in get_components(component_type, sys)
        if get_name(component) == get_name(other)
            return component
        end
    end

    error("Did not find component $component")
end

"""Return the Branch in the system that matches another by case-insensitive arc
names."""
function get_branch(sys::System, other::Branch)
    for branch in get_components(Branch, sys)
        if lowercase(other.arc.from.name) == lowercase(branch.arc.from.name) &&
           lowercase(other.arc.to.name) == lowercase(branch.arc.to.name)
            return branch
        end
    end

    error("Did not find branch with buses $(other.arc.from.name) ", "$(other.arc.to.name)")
end

function create_system_with_dynamic_inverter()
    nodes_OMIB = [
        Bus(
            1, #number
            "Bus 1", #Name
            "REF", #BusType (REF, PV, PQ)
            0, #Angle in radians
            1.06, #Voltage in pu
            (min = 0.94, max = 1.06), #Voltage limits in pu
            69,
            nothing,
            nothing,
        ), #Base voltage in kV
        Bus(2, "Bus 2", "PV", 0, 1.045, (min = 0.94, max = 1.06), 69, nothing, nothing),
    ]

    battery = GenericBattery(
        name = "Battery",
        prime_mover = PrimeMovers.BA,
        available = true,
        bus = nodes_OMIB[2],
        initial_energy = 5.0,
        state_of_charge_limits = (min = 5.0, max = 100.0),
        rating = 0.0275, #Value in per_unit of the system
        active_power = 0.01375,
        input_active_power_limits = (min = 0.0, max = 50.0),
        output_active_power_limits = (min = 0.0, max = 50.0),
        reactive_power = 0.0,
        reactive_power_limits = (min = -50.0, max = 50.0),
        efficiency = (in = 0.80, out = 0.90),
        base_power = 100.0,
    )
    converter = AverageConverter(
        138.0, #Rated Voltage
        100.0,
    ) #Rated MVA

    branch_OMIB = [
        Line(
            "Line1", #name
            true, #available
            0.0, #active power flow initial condition (from-to)
            0.0, #reactive power flow initial condition (from-to)
            Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
            0.01, #resistance in pu
            0.05, #reactance in pu
            (from = 0.0, to = 0.0), #susceptance in pu
            18.046, #rate in MW
            1.04,
        ),
    ]  #angle limits (-min and max)

    dc_source = FixedDCSource(1500.0) #Not in the original data, guessed.

    filt = LCLFilter(
        0.08, #Series inductance lf in pu
        0.003, #Series resitance rf in pu
        0.074, #Shunt capacitance cf in pu
        0.2, #Series reactance rg to grid connection (#Step up transformer or similar)
        0.01,
    ) #Series resistance lg to grid connection (#Step up transformer or similar)

    pll = KauraPLL(
        500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
        0.084, #k_p: PLL proportional gain
        4.69,
    ) #k_i: PLL integral gain

    virtual_H = VirtualInertia(
        2.0, #Ta:: VSM inertia constant
        400.0, #kd:: VSM damping coefficient
        20.0, #kω:: Frequency droop gain in pu
        2 * pi * 50.0,
    ) #ωb:: Rated angular frequency

    Q_control = ReactivePowerDroop(
        0.2, #kq:: Reactive power droop gain in pu
        1000.0,
    ) #ωf:: Reactive power cut-off low pass filter frequency

    outer_control = OuterControl(virtual_H, Q_control)

    vsc = CurrentControl(
        0.59, #kpv:: Voltage controller proportional gain
        736.0, #kiv:: Voltage controller integral gain
        0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
        0.0, #rv:: Virtual resistance in pu
        0.2, #lv: Virtual inductance in pu
        1.27, #kpc:: Current controller proportional gain
        14.3, #kiv:: Current controller integral gain
        0.0, #kffi:: Binary variable enabling the current feed-forward in output of current controllers
        50.0, #ωad:: Active damping low pass filter cut-off frequency
        0.2,
    ) #kad:: Active damping gain

    sys = System(100)
    for bus in nodes_OMIB
        add_component!(sys, bus)
    end
    for lines in branch_OMIB
        add_component!(sys, lines)
    end
    add_component!(sys, battery)

    test_inverter = DynamicInverter(
        get_name(battery),
        1.0, #ω_ref
        converter, #Converter
        outer_control, #OuterControl
        vsc, #Voltage Source Controller
        dc_source, #DC Source
        pll, #Frequency Estimator
        filt,
    ) #Output Filter

    add_component!(sys, test_inverter, battery)

    return sys
end

function test_accessors(component)
    ps_type = typeof(component)
    @debug "test_accessors" ps_type

    for (field_name, field_type) in zip(fieldnames(ps_type), fieldtypes(ps_type))
        if field_name === :name || field_name === :time_series_container
            func = getfield(InfrastructureSystems, Symbol("get_" * string(field_name)))
            _func! =
                getfield(InfrastructureSystems, Symbol("set_" * string(field_name) * "!"))
        else
            getter_name = Symbol("get_" * string(field_name))
            if !hasproperty(PowerSystems, getter_name)
                continue
            end
            func = getfield(PowerSystems, getter_name)
            if !hasmethod(func, (ps_type,))
                continue
            end
            setter_name = Symbol("set_" * string(field_name) * "!")
            # In some cases there is a getter but no setter.
            if hasproperty(PowerSystems, setter_name)
                _func! = getfield(PowerSystems, setter_name)
            else
                _func! = nothing
            end
        end

        val = func(component)
        @test val isa field_type
        try
            if typeof(val) == Float64 || typeof(val) == Int
                if !isnan(val)
                    aux = val + 1
                    if _func! !== nothing
                        _func!(component, aux)
                        @test func(component) == aux
                    end
                end
            elseif typeof(val) == String
                aux = val * "1"
                if _func! !== nothing
                    _func!(component, aux)
                    @test func(component) == aux
                end
            elseif typeof(val) == Bool
                aux = !val
                if _func! !== nothing
                    _func!(component, aux)
                    @test func(component) == aux
                end
            else
                _func! !== nothing && _func!(component, val)
            end
        catch MethodError
            continue
        end
    end
end

function validate_serialization(
    sys::System;
    time_series_read_only = false,
    runchecks = false,
)
    test_dir = mktempdir()
    orig_dir = pwd()
    cd(test_dir)

    try
        path = joinpath(test_dir, "test_system_serialization.json")
        @info "Serializing to $path"
        sys_ext = get_ext(sys)
        sys_ext["data"] = 5
        ext_test_bus_name = ""
        IS.prepare_for_serialization!(sys.data, path; force = true)
        bus = collect(get_components(PSY.Bus, sys))[1]
        ext_test_bus_name = PSY.get_name(bus)
        ext = PSY.get_ext(bus)
        ext["test_field"] = 1
        to_json(sys, path)

        data = open(path, "r") do io
            JSON3.read(io)
        end
        @test data["data_format_version"] == PSY.DATA_FORMAT_VERSION

        sys2 = System(
            path;
            time_series_read_only = time_series_read_only,
            runchecks = runchecks,
        )
        isempty(get_bus_numbers(sys2)) && return false
        sys_ext2 = get_ext(sys2)
        sys_ext2["data"] != 5 && return false
        bus = PSY.get_component(PSY.Bus, sys2, ext_test_bus_name)
        ext = PSY.get_ext(bus)
        ext["test_field"] != 1 && return false
        return sys2, IS.compare_values(sys, sys2)
    finally
        cd(orig_dir)
    end
end


text:
## Julia program for Adaptive Rejection Sampling
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-20
##
## refer to Inferentialist's R code:
##   https://blog.inferentialist.com/2016/09/26/adaptive-sampling.html

module corears
# example 3.22 in Davison(2008)
r = 2; m = 10; mu = 0; sig2 = 1;
# range of y
yl = -2; yu = 2;
# function of h
function h(y)
    return(r * y - m * log(1 + exp(y)) - (y - mu)^2 / (2 * sig2))
end

function h(y::Array)
    return(r * y - m * log.(1 .+ exp.(y)) .- (y .- mu) .^2 ./ (2 * sig2))
end

# derivative of h
function dh(y)
    return(r - m * exp(y) / (1 + exp(y)) - (y - mu) / sig2)
end

function dh(y::Array)
    return(r .- m * exp.(y) ./ (1 .+ exp.(y)) .- (y .- mu)./sig2)
end

# intersection point
function zfix(yfixed::Array)
    yf0 = yfixed[1:end-1]
    yf1 = yfixed[2:end]
    zfixed = yf0 .+ (h(yf0) .- h(yf1) .+ (yf1 .- yf0) .* dh(yf1)) ./ (dh(yf1) .- dh(yf0))
    return(zfixed)
end

# evaluate log-density
function hplus(y, yfixed::Array)
    zfixed = zfix(yfixed)
    n = size(zfixed, 1)
    for i = 1:n
        if i == 1 && y < zfixed[i]
            return(h(yfixed[i]) + (y - yfixed[i]) * dh(yfixed[i]))
        elseif i < n && y >= zfixed[i] && y < zfixed[i+1]
            return(h(yfixed[i+1]) + (y - yfixed[i+1]) * dh(yfixed[i+1]))
        elseif i == n && y >= zfixed[n]
            return(h(yfixed[i+1]) + (y - yfixed[i+1]) * dh(yfixed[i+1]))
        end
    end
end

# calculate G_+(z_i)
function gplus_cdf(yfixed::Array, zfixed::Array)
    n = size(zfixed, 1)
    s = zeros(n+1)
    pr = zeros(n+1)
    for i = 1:(n+1)
        ## integral from -infty to zi
        if i == 1
        #    s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - )
            s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - exp((yl-yfixed[i]) * dh(yfixed[i])))
        elseif i == n+1
            s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((yu-yfixed[i]) * dh(yfixed[i])) - exp((zfixed[n]-yfixed[i]) * dh(yfixed[i])))
        else
            s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])))
        end
    end
    pr = s / sum(s)
    return cumsum(pr), sum(s)
end

# sample from gplus density
function gplus_sample(yfixed::Array)
    zfixed = zfix(yfixed)
    gp = gplus_cdf(yfixed, zfixed)
    zpr = gp[1]
    norm_const = gp[2]
    n = size(zfixed, 1)
    u = rand()
    # Invert the gplus pdf
    for i = 1:n
        if i == 1 && u < zpr[i]
            ey = u * dh(yfixed[i]) * norm_const / exp(h(yfixed[i])) + exp((yl-yfixed[i])*dh(yfixed[i]))
            return(yfixed[i] + log(ey)/dh(yfixed[i]))
        elseif i == n && u >= zpr[i]
            ey = (u - zpr[i]) * dh(yfixed[i+1]) * norm_const / exp(h(yfixed[i+1])) + exp((zfixed[i]-yfixed[i+1])*dh(yfixed[i+1]))
            return(yfixed[i+1] + log(ey)/dh(yfixed[i+1]))
        elseif i < n && u >= zpr[i] && u < zpr[i+1]
            ey = (u - zpr[i]) * dh(yfixed[i+1]) * norm_const / exp(h(yfixed[i+1])) + exp((zfixed[i]-yfixed[i+1])*dh(yfixed[i+1]))
            return(yfixed[i+1] + log(ey)/dh(yfixed[i+1]))
        end
    end
end

## adaptive rejection sampling
function ars(yfixed::Array)
    x = gplus_sample(yfixed)
    u = rand()
    if u <= exp(h(x)-hplus(x, yfixed))
        return(x)
    else
        return(ars(append!(yfixed, x)))
    end
end

## example
N = 100
res = ones(N);
for i = 1:N
    res[i] = ars([-1.8,-1.1,-0.5,-0.2])
end

export zfix, hplus, gplus_cdf, gplus_sample
end

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

 ## Basic functions ##

isreal(x::AbstractQuasiArray) = all(isreal,x)
iszero(x::AbstractQuasiArray) = all(iszero,x)
isreal(x::AbstractQuasiArray{<:Real}) = true
all(::typeof(isinteger), ::AbstractQuasiArray{<:Integer}) = true

## Constructors ##

dropdims(A; dims) = _dropdims(A, dims)
function _dropdims(A::AbstractQuasiArray, dims::Dims)
    for i in 1:length(dims)
        1 <= dims[i] <= ndims(A) || throw(ArgumentError("dropped dims must be in range 1:ndims(A)"))
        length(axes(A, dims[i])) == 1 || throw(ArgumentError("dropped dims must all be size 1"))
        for j = 1:i-1
            dims[j] == dims[i] && throw(ArgumentError("dropped dims must be unique"))
        end
    end
    d = ()
    for i = 1:ndims(A)
        if !in(i, dims)
            d = tuple(d..., axes(A, i))
        end
    end
    reshape(A, d::typeof(_sub(axes(A), dims)))
end
_dropdims(A::AbstractQuasiArray, dim::Integer) = _dropdims(A, (Int(dim),))

## Unary operators ##

conj(x::AbstractQuasiArray{<:Real}) = x
conj!(x::AbstractQuasiArray{<:Real}) = x

real(x::AbstractQuasiArray{<:Real}) = x
imag(x::AbstractQuasiArray{<:Real}) = zero(x)

+(x::AbstractQuasiArray{<:Number}) = x
*(x::AbstractQuasiArray{<:Number,2}) = x

# index A[:,:,...,i,:,:,...] where "i" is in dimension "d"

@inline selectdim(A::AbstractQuasiArray, d::Integer, i) = _selectdim(A, d, i, setindex(map(Slice, axes(A)), i, d))
@noinline function _selectdim(A, d, i, idxs)
    d >= 1 || throw(ArgumentError("dimension must be ≥ 1"))
    nd = ndims(A)
    d > nd && (i == 1 || throw(BoundsError(A, (ntuple(k->Colon(),d-1)..., i))))
    return view(A, idxs...)
end


text:
# The purpose of this script is to compare the filtering capability of the single snapshots
# estimates and the gas filter on a given

#function that Generates the dynamical fitnesses and estimates single snapshot

using Distributed
addprocs(4)

## Simulate estimate and save
@everywhere push!(LOAD_PATH,"../../DynNets/src/")
@everywhere push!(LOAD_PATH,"../../../DynNets/src/")
@everywhere push!(LOAD_PATH,"../DynNets/src/")
@everywhere push!(LOAD_PATH,"./DynNets/src/")

@everywhere using SharedArrays, Statistics, Utilities, DynNets, StaticNets, JLD

@everywhere begin
    function simWrap1(N,T, NTV;dynType = "SIN", degIOUncMeans =0, degb = [10, 40],
                      dynFitDgp = zeros(10,50), indsTVnodes=falses(10))

       if sum(dynFitDgp)==0
           dynFitDgp,indsTVnodes =  StaticNets.dgpDynamic(StaticNets.fooErgmDirBin1,dynType,N,T;
                                               NTV = NTV,degIOUncMeans = degIOUncMeans,degb = degb  )
       end
       dynExpDegsDgp = zero(dynFitDgp)
       for t =1:T  dynExpDegsDgp[:,t] = StaticNets.expValStats(StaticNets.fooErgmDirBin1,dynFitDgp[:,t]) end
       # sanmple dynamical networks from dgp fitnesses
       dynDegsSam = zeros(size(dynFitDgp))
       for t=1:T
           matSamp = StaticNets.sampl(StaticNets.ErgmDirBin1(dynFitDgp[:,t]),1;  parGroupsIO = dynFitDgp[:,t])
           dynDegsSam[:,t] = [sumSq(matSamp,2);sumSq(matSamp,1)]
       end
       #println(dynDegsSam[:,4])
       estFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(dynDegsSam); identPost = true,identIter= false,targetErr = 0.001 )
       for t=1:T
           dynFitDgp[:,t] =  StaticNets.identify(StaticNets.ErgmDirBin1(dynDegsSam[:,t]),
                                              dynFitDgp[:,t]; idType = "equalIOsums" )
           estFitSS[:,t] =  StaticNets.identify(StaticNets.ErgmDirBin1(dynDegsSam[:,t]),
                                              estFitSS[:,t]; idType =  "equalIOsums")
       end
       return dynFitDgp,dynDegsSam,estFitSS,indsTVnodes
       end

    N = 10# [500,1000,1500] #[50,60,80,90,100,125,150,200,300,400]# 200##40 #
  minDeg,maxDeg = (3,8)
  degb = [round(Int,minDeg), round(Int,maxDeg)]#
  unifDeg = round(Int,(maxDeg-  minDeg)/2)
  dgpType = "STEPS"
  Nsample = 100
  targetStatic = true
  T=250
    dynDegsSam = SharedArray{Float64,3}(2N,T,Nsample)
    estFitSS = SharedArray{Float64,3}(2N,T,Nsample)
    #dynFitDgp = Array{Array{Float64,2},2}(Nsample,2)
    filFitGas = SharedArray{Float64,3}(2N,T,Nsample)
    rmseSSandGas =  SharedArray{Float64,3}(2N,2,Nsample)
    storeGasPar = SharedArray{Float64,3}(2N,3,Nsample)#  SharedArray{Array{Array{Float64,1},1},2}(Nsample,2)
    # convGasFlag = SharedArrayBitArray{2}(Nsample,2)

     NTV = round.(Int,N/2)

     tmprmseSSandGas = SharedArray{Float64,3}(2N,2,Nsample) #zeros(2N,2,Nsample)
    prog=0
     dynFitDgp,~,~,indsTVnodes =
           simWrap1(N,T,NTV;dynType =dgpType,degIOUncMeans =unifDeg*ones(2N),degb = degb)
end

@sync @distributed for s=1:Nsample
        println((N,s))
        #sample the dgp
        @time ~,degsIO_T,estFitSS_T,~ =
               simWrap1(N,T,NTV;dynType =dgpType,degIOUncMeans =unifDeg*ones(2N),degb = degb, dynFitDgp = dynFitDgp,indsTVnodes = indsTVnodes)
       fooPar =[zeros(2N), zeros(2N) ,zeros(2N)]# [zeros(2N),zeros(1),zeros(1)]#
       indsGroups = [Int.(1:2N),Int.(1:2N)] #[Int.(1:2N),ones(2N)] #
       modGasDirBin1 = DynNets.SdErgmDirBin1(degsIO_T, fooPar, indsGroups,"FISHER-DIAG")

       rmseSSandGas[:,1,s] =sqrt.(meanSq((dynFitDgp - estFitSS_T).^2,2))
       estFitSS[:,:,s] = estFitSS_T
       if targetStatic
          estFitStatic = DynNets.estimateSnapSeq(modGasDirBin1, degsIO_T = mean(degsIO_T,dims=2) )
          estFitSS_T = repeat(estFitStatic, 1, 5)
        end
       @time        parGas,~  =  DynNets.estimateTarg(modGasDirBin1;SSest =estFitSS_T)

        #estimate single snapshot sequence
        #estimate score_driven_filter_or_dgp
        ## e ora stima il gasPar
        dynDegsSam[:,:,s] = degsIO_T
          storeGasPar[:,1,s] = parGas[1]; storeGasPar[:,2,s] = parGas[2]; storeGasPar[:,3,s] = parGas[3]
          display(parGas[2:3])
          gasFiltFit,~ = DynNets.score_driven_filter_or_dgp(modGasDirBin1,[parGas[1];parGas[2];parGas[3]])
          rmseSSandGas[:,2,s] =sqrt.(meanSq((dynFitDgp - gasFiltFit).^2,2))
          filFitGas[:,:,s] = gasFiltFit

      end

      save_fold = "./data/sdergm_estimates/"
      @save(save_fold*"FilterTestDynamicFixedSize_N_$(N)_Nsample_$(Nsample)_scaling_DGP_" * dgpType  *
        "targetStatic_$(targetStatic)"* ".jld",
              dynFitDgp,dynDegsSam,estFitSS,indsTVnodes,storeGasPar,rmseSSandGas,filFitGas,degb,unifDeg)

##
if false
  using PyCall; pygui(:qt); using PyPlot; pygui(true)

indsPlot = [ 6]
simsToDrop = dropdims(sum(isnan.(abs.(filFitGas)), dims = (1,2)).>0, dims = (1,2))
  Nplot =indsPlot[1]
    plot(1:T,dynFitDgp[Nplot,:], "k")
    plot(1:T,estFitSS[Nplot,:,.!simsToDrop],".b", alpha=0.5)
    plot(1:T, filFitGas[Nplot,:,.!simsToDrop],".r", alpha = 0.3)

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM