text:
# MathOptInterface constraints

##################################################
# Generic constraint definition
#--------------------------------------------------
# Bound constraint on variables.

# TODO(Xhub) is this sufficient?
function chk_var_cons(model::Optimizer, vi::MOI.VariableIndex)
    check_inbounds(model, vi)
end

function chk_var_cons(model::Optimizer, vi::MOI.VariableIndex, rel::SS)
    check_inbounds(model, vi)

    if isa(rel, MOI.LessThan{Float64})
        if isnan(rel.upper)
           error("Invalid upper bound value $(rel.upper).")
       end
        if has_upper_bound(model, vi)
            error("Upper bound on variable $vi already exists.")
        end
        if is_fixed(model, vi)
            error("Variable $vi is fixed. Cannot also set upper bound.")
        end
    elseif isa(rel, MOI.GreaterThan{Float64})
        if isnan(rel.lower)
            error("Invalid lower bound value $(rel.lower).")
        end
        if has_lower_bound(model, vi)
            error("Lower bound on variable $vi already exists.")
        end
        if is_fixed(model, vi)
            error("Variable $vi is fixed. Cannot also set lower bound.")
        end
    elseif isa(rel, MOI.EqualTo{Float64})
        if isnan(rel.value)
            error("Invalid fixed value $(eq.value).")
        end
        if has_lower_bound(model, vi)
            error("Variable $vi has a lower bound. Cannot be fixed.")
        end
        if has_upper_bound(model, vi)
            error("Variable $vi has an upper bound. Cannot be fixed.")
        end
        if is_fixed(model, vi)
            error("Variable $vi is already fixed.")
        end
    elseif isa(rel, MOI.Interval{Float64})
        if isnan(rel.lower)
            error("Invalid lower value $(eq.lower).")
        end
        if isnan(rel.upper)
            error("Invalid upper value $(eq.upper).")
        end

        if has_lower_bound(model, vi)
            error("Variable $vi already has a lower bound.")
        end
        if has_upper_bound(model, vi)
            error("Variable $vi already has an upper bound.")
        end
        if is_fixed(model, vi)
            error("Variable $vi is already fixed.")
        end
    end

end

function reshop_set_bound(ctx::Ptr{context}, vidx::Integer, rel::MOI.LessThan{Float64})
    ctx_setvarub(ctx, vidx, rel.upper)
end

function reshop_set_bound(ctx::Ptr{context}, vidx::Integer, rel::MOI.GreaterThan{Float64})
    ctx_setvarlb(ctx, vidx, rel.lower)
end

function reshop_set_bound(ctx::Ptr{context}, vidx::Integer, rel::MOI.EqualTo{Float64})
    ctx_setvarfx(ctx, vidx, rel.value)
end

function reshop_set_bound(ctx::Ptr{context}, vidx::Integer, rel::MOI.Interval{Float64})
    ctx_setvarub(ctx, vidx, rel.upper)
    ctx_setvarlb(ctx, vidx, rel.lower)
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, rel::SS)
    vi = v.variable
    chk_var_cons(model, vi, rel)
    reshop_set_bound(model.ctx, vi.value-1, rel)
    return MOI.ConstraintIndex{MOI.SingleVariable, typeof(rel)}(vi.value)
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.ScalarAffineFunction{Float64}, set::LS)
    check_inbounds(model, func)
    eidx = rhp_add_equ(model.ctx)
    # Parse structure of constraint.
    vidx, coefs = canonical_linear_reduction(func)
    avar = _ensure_avar(model)
    rhp_avar_set(avar, vidx)
    rhp_equ_add_linear(model.ctx, eidx, avar, coefs)
    _set_rhs(model, eidx, func.constant, set)
    # Add constraint to index.
    return MOI.ConstraintIndex{typeof(func), typeof(set)}(eidx+1)
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.ScalarQuadraticFunction{Float64}, set::LS)
    check_inbounds(model, func)
    eidx = rhp_add_equ(model.ctx)

    # Parse linear structure of constraint.
    # Parse quadratic term.
    qvar1, qvar2, qcoefs = canonical_quadratic_reduction(func)

    rhp_equ_add_quadratic(model.ctx, eidx, qvar1, qvar2, qcoefs)
    vidx, coefs = canonical_linear_reduction(func)

    avar = _ensure_avar(model)
    rhp_avar_set(avar, vidx)
    # TODO(Xhub) change this when quadratic constrains are supported
    rhp_equ_add_lin_tree(model.ctx, eidx, coefs, avar, 1.)
    _set_rhs(model, eidx, func.constant, set)
    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(eidx+1)
    model.quadfn_mapping[ci] = eidx
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.VectorAffineFunction, set::MOI.AbstractVectorSet)
    # TODO: add check inbounds for VectorAffineFunction.
    current_m = ctx_numequ(model.ctx)
    num_cons = length(func.constants)
    # Add constraints inside ReSHOP.
    rhp_add_equs(model.ctx, num_cons)

    # Parse vector affine expression.
    # The indices are already 0-based
    cidx, vidx, coeffs = canonical_vector_affine_reduction(func)
    # Get true indices
    cidx .+= current_m

    for (idx, eidx) in enumerate(cidx)
        ctx_add_lin_var(model.ctx, eidx, vidx[idx], coeffs[idx])
    end

    for (idx, cst) in enumerate(func.constants)
        _set_rhs(model, current_m-1+idx, cst, set)
    end

    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(cidx[1]+1)
    model.vaf_mapping[ci] = unique(cidx)
    return ci
end

# Add second order cone constraint.
#function MOI.add_constraint(model::Optimizer,
#                            func::MOI.VectorAffineFunction, set::MOI.SecondOrderCone)
#    # Add constraints inside ReSHOP.
#    rhp_add_equ(model.ctx)
#    eidx = ctx_numequ(model.ctx)
#
#    # Parse vector affine expression.
#    indexcoords, vidx, coefs = canonical_vector_affine_reduction(func)
#    constants = func.constants
#    # Distinct two parts of secondordercone.
#    # First row corresponds to linear part of SOC.
#    indlinear = indexcoords .== 0
#    indcone = indexcoords .!= 0
#    ncoords = length(constants) - 1
#    @assert ncoords == set.dimension - 1
#
#    # Load Second Order Conic constraint.
#    ## i) linear part
#    KN_set_con_upbnd(model.ctx, eidx, constants[1])
#    rhp_equ_add_linear_struct(model.ctx, eidx,
#                             vidx[indlinear], -coefs[indlinear])
#
#    ## ii) soc part
#    index_var_cone = vidx[indcone]
#    nnz = length(index_var_cone)
#    index_coord_cone = convert.(Cint, indexcoords[indcone] .- 1)
#    coefs_cone = coefs[indcone]
#    const_cone = constants[2:end]
#
#    rhp_equ_add_L2norm(model.ctx,
#                      eidx, ncoords, nnz,
#                      index_coord_cone,
#                      index_var_cone,
#                      coefs_cone,
#                      const_cone)
#
#    # Add constraints to index.
#    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(eidx)
#    model.constraint_mapping[ci] = vidx
#    return ci
#end

function MOI.add_constraint(model::Optimizer, vov::VOV, set::T) where T <: VLS
    vidx = convert.(Cint, [v.value for v in vov.variables] .- 1)
    chk_var_cons.(model, vov.variables)
    bnd = zeros(Float64, length(vidx))

    if isa(set, MOI.Zeros)
        ctx_setvarfx.(model.ctx, vidx, bnd)
    elseif isa(set, MOI.Nonnegatives)
        ctx_setvarlb.(model.ctx, vidx, bnd)
    elseif isa(set, MOI.Nonpositives)
        ctx_setvarub.(model.ctx, vidx, bnd)
    end

    # TODO is this necessary?
#    ncons = MOI.get(model, MOI.NumberOfConstraints{VOV, T}())
    ci = MOI.ConstraintIndex{VOV, T}(vidx[1]+1)
    model.vov_mapping[ci] = vidx
    return ci
end

function MOI.add_constraint(model::Optimizer, vov::VOV, set::MOI.SecondOrderCone)

    # the variables have the following structure
    # (t, x) where t ≥ ||x||_2
    #

    vidx = [v.value - 1 for v in func.variables]

    soc_data = rhp_cone_soc(vidx[2:end])
    reshop_set_vartype(model.ctx, vidx[1], RHP_VARTYPE_SOC, soc_data)

    # Add constraints to index.
    # The index of the t variable is used, since the multiplier
    # would only be defined for it.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(vidx[1]+1)
    model.vov_mapping[ci] = vidx
    return ci
end

function MOI.add_constraint(model::Optimizer, vov::VOV,
                            set::Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}})
    @assert length(vov.variables) == length(set.weights)
    vidx = convert.(Cint, [v.value for v in vov.variables] .- 1)
    avar = _ensure_avar(model)
    rhp_avar_set(avar, vidx)
    reshop_set_vartype(model.ctx, avar, set, set.weights)
    ci = MOI.ConstraintIndex{VOV, typeof(set)}(vidx[1]+1)
    model.vov_mapping[ci] = vidx
    model.sos_sets[ci] = set
    return ci
end

##################################################
## Change some constraints
function MOI.set(model::Optimizer, ::MOI.ConstraintSet, ci::MOI.ConstraintIndex{MOI.SingleVariable, S}, new_set::S) where S <: SS
    reshop_set_bound(model.ctx, ci.value-1, new_set);
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S <: SS
    vidx = ci.value-1
    ctx_setvarlb(model.ctx, vidx, -Inf64)
    ctx_setvarub(model.ctx, vidx, Inf64)
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, <:Union{MOI.ZeroOne, MOI.Integer}})
    reshop_set_vartype(model.ctx, ci.value - 1, RHP_VARTYPE_X)
end

##################################################
## Binary & Integer constraints.

# Define integer and boolean constraints.
# Those don't need to be stored in the constraint index
function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, S::Union{MOI.ZeroOne, MOI.Integer})
    vi = v.variable
    check_inbounds(model, vi)
    reshop_set_vartype(model.ctx, vi.value - 1, set_to_reshop[typeof(S)])
    return MOI.ConstraintIndex{MOI.SingleVariable, typeof(S)}(vi.value)
end

##################################################
## Constraint DualStart
function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable, S}, value::Union{Real, Nothing}) where S <: SS

    val = isa(value, Real) ? Cdouble(value) : 0.
    ctx_setvarmult(model.ctx, ci.value - 1, val)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex, value::Union{Real, Nothing})
    check_inbounds(model, ci)
    val = isa(value, Real) ? Cdouble(value) : 0.
    if ci in model.vaf_mapping
        ctx_setequmult(model.ctx, ci.value - 1, val)
    else
        error("Unknown contraint $ci")
    end
    return
end

function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{VOV, T}, value::Union{Real, Nothing}) where T
    check_inbounds(model, ci)
    val = isa(value, Real) ? Cdouble(value) : 0.

    if ci in model.vov_mapping
        for vidx in model.vov_mapping[ci]
            ctx_setvarmult(model.ctx, vidx, val)
        end
    else
        error("Unknown contraint $ci")
    end
    return
end

function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{VOV, T}, value::AbstractVector{Real}) where T
    check_inbounds(model, ci)

    if ci in model.vov_mapping
        @assert length(value) == length(model.vov_mapping[ci])
        for (idx, vidx) in enumerate(model.vov_mapping[ci])
            ctx_setvarmult(model.ctx, vidx, Cdouble(value[idx]))
        end
    else
        error("Unknown contraint $ci")
    end
    return
end

function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{VAF, T}, value::Union{Real, Nothing}) where T
    check_inbounds(model, ci)
    val = isa(value, Real) ? Cdouble(value) : 0.

    if ci in model.vaf_mapping
        for eidx in model.vaf_mapping[ci]
            ctx_setequmult(model.ctx, eidx, val)
        end
    else
        error("Unknown contraint $ci")
    end
    return
end

function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{VAF, T}, value::AbstractVector{Real}) where T
    check_inbounds(model, ci)

    if ci in model.vaf_mapping
        @assert length(value) == length(model.vaf_mapping[ci])
        for (idx, eidx) in enumerate(model.vaf_mapping[ci])
            ctx_setequmult(model.ctx, eidx, Cdouble(value[idx]))
        end
    else
        error("Unknown contraint $ci")
    end
    return
end

##################################################
# getters
#

function MOI.get(model::Optimizer, loc::MOI.ListOfConstraints)
    @warn "MOI.get(model::Optimizer, loc::MOI.ListOfConstraints) is not yet implemented for ReSHOP"
    list = Vector{Tuple{DataType,DataType}}()
    return list
    for S in (
        MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.Interval{Float64},
        MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64}, MOI.Integer, MOI.ZeroOne
       )
        nb = MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable, S}())
        if nb > 0
            push!(list, (MOI.SingleVariable, S))
        end
    end
    for S in (MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros, MOI.SecondOrderCone)
        nb = MOI.get(model, MOI.NumberOfConstraints{VAF, S}())
        if nb > 0
            push!(list, (VAF, S))
        end
    end
    for S in (MOI.Nonnegatives,  MOI.Nonpositives, MOI.Zeros, MOI.SOS1, MOI.SOS2)
         nb = MOI.get(model, MOI.NumberOfConstraints{VOV, S}())
         if nb > 0
             push!(list, (VOV, S))
         end
    end
    for F in (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction)
        for S in (MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64})
            nb = MOI.get(model, MOI.NumberOfConstraints{F, S}())
            if nb > 0
                push!(list, (F, S))
            end
        end
    end

    return list
#    for S in (MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, MOI.ExponentialCone,
#              MOI.DualExponentialCone, MOI.PowerCone, MOI.DualPowerCone)
#        nb = MOI.get(model, MOI.ScalarAffineFunction
end

function _get_var_set(model::Optimizer, S::Type{<:Union{MOI.EqualTo, MOI.GreaterThan}}, lb::Float64, ub::Float64)
    return S(lb)
end

function _get_var_set(model::Optimizer, S::Type{<:MOI.LessThan}, lb::Float64, ub::Float64)
    return S(ub)
end

# TODO semiinteger should have a different codepath
function _get_var_set(model::Optimizer, S::Type{<:Union{MOI.Interval, MOI.Semicontinuous, MOI.Semiinteger}}, lb::Float64, ub::Float64)
    return S(lb, ub)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintSet, ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S
    lb, ub = ctx_getvarbounds(model.ctx, ci.value-1)
    return _get_var_set(model, S, lb, ub)
end

##################################################
#Constraint validity check
#

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex)
    error("MOI.is_valid not implemented for type $(ci)")
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}})
    return has_upper_bound(model, MOI.VariableIndex(ci.value)) && !is_fixed(model, MOI.VariableIndex(ci.value))
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}})
    return has_lower_bound(model, MOI.VariableIndex(ci.value)) && !is_fixed(model, MOI.VariableIndex(ci.value))
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}})
    return is_fixed(model, MOI.VariableIndex(ci.value))
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}})
    return has_box_bounds(model, MOI.VariableIndex(ci.value))
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer})
    return reshop_getvartype(model.ctx, ci.value-1) == RHP_VARTYPE_I
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.ZeroOne})
    return reshop_getvartype(model.ctx, ci.value-1) == RHP_VARTYPE_B
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Semicontinuous{Float64}})
    return reshop_getvartype(model.ctx, ci.value-1) == RHP_VARTYPE_SC
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Semiinteger{Float64}})
    return reshop_getvartype(model.ctx, ci.value-1) == RHP_VARTYPE_SI
end

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{VOV, <:VLS})
    allvars = get(model.vov_mapping, ci, Cuint[])
    if length(allvars) == 0
        return false
    end

    return !any([MOI.VariableIndex(1+vi) in model.delvars for vi in allvars])
end

##################################################
## Constraint naming
# TODO
function MOI.set(model::Optimizer, ::MOI.ConstraintName, ci::MOI.ConstraintIndex{<:SF,<:LS}, name::String)
    println("saving constraint named $name at index $(ci.value)")
    ctx_setequname(model.ctx, ci.value-1, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, ci::MOI.ConstraintIndex{MOI.SingleVariable,<:SS}, name::String)
    ctx_setequname(model.ctx, ci.value-1, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, ci::MOI.ConstraintIndex, name::String)
    println("unsupported constraint named $name at index $(ci.value) of type $ci")
    model.fake_cons_name[ci] = name
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ci::MOI.ConstraintIndex{<:SF, <:LS})
    return ctx_getequname(model.ctx, ci.value-1)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ci::MOI.ConstraintIndex)
    return get(model.fake_cons_name, ci, "")
end

function MOI.get(model::Optimizer, ::MOI.ConstraintFunction, ci::MOI.ConstraintIndex{MOI.SingleVariable,<: SS})
    MOI.throw_if_not_valid(model, ci)
    return MOI.SingleVariable(MOI.VariableIndex(ci.value))
end


text:
function arg_decl(x::Spec)
    T = idiomatic_julia_type(x)
    T in struct_name.(spec_handles, true) && return wrap_identifier(x)
    :($(wrap_identifier(x))::$(signature_type(T)))
end

function kwarg_decl(x::Spec)
    val = default(x)
    Expr(:kw, wrap_identifier(x), val)
end

drop_arg(x::Spec) = is_length(x) && !is_length_exception(x) && is_inferable_length(x) || is_pointer_start(x) || x.type == :(Ptr{Ptr{Cvoid}})

"""
Function pointer arguments for a handle.
Includes one `fptr_create` for the constructor (if applicable),
and one `fptr_destroy` for the destructor (if applicable).
"""
function func_ptr_args(spec::SpecHandle)
    args = Expr[]
    spec ∈ spec_create_funcs.handle && push!(args, :(fptr_create::FunctionPtr))
    destructor(spec) ≠ :identity && push!(args, :(fptr_destroy::FunctionPtr))
    args
end

"""
Function pointer arguments for a function.
Takes the function pointers arguments of the underlying handle if it is a Vulkan constructor,
or a unique `fptr` if that's just a normal Vulkan function.
"""
function func_ptr_args(spec::SpecFunc)
    if spec.type ∈ [FTYPE_CREATE, FTYPE_ALLOCATE]
        func_ptr_args(create_func(spec).handle)
    else
        [:(fptr::FunctionPtr)]
    end
end

"""
Corresponding pointer argument for a Vulkan function.
"""
func_ptrs(spec::Spec) = name.(func_ptr_args(spec))

function add_func_args!(p::Dict, spec, params; with_func_ptr = false)
    params = filter(!drop_arg, params)
    arg_filter = if spec.type ∈ [FTYPE_DESTROY, FTYPE_FREE]
        destroyed_type = destroy_func(spec).handle.name
        x -> !is_optional(x) || x.type == destroyed_type
    else
        !is_optional
    end

    p[:args] = convert(Vector{ExprLike}, map(arg_decl, filter(arg_filter, params)))
    p[:kwargs] = map(kwarg_decl, filter(!arg_filter, params))

    with_func_ptr && append!(p[:args], func_ptr_args(spec))
end

ext: jl


text:
using NSGAIII
using Base.Test

const C1 = [2 5 4 7 ; 3 3 5 7 ; 3 8 4 2 ; 6 5 2 5]
const C2 = [3 3 6 2 ; 5 3 7 3 ; 5 2 7 4 ; 4 6 3 5]
const C3 = [4 2 5 3 ; 5 3 4 3 ; 4 3 5 2 ; 6 4 7 3]

z(x, C) = sum(inds->C[inds...], enumerate(x))
z(x::Vector{Int}) = z(x, C1), z(x, C2), z(x, C3)

res = unique(nsga(500, 100, ()->randperm(4), z, 15))
display(res)
@test length(res) == 7

const d = RealCoding(6, [-10], [10])
z1(x) = x^2
z2(x) = (x-2)^2
z(bits) = begin
    x = decode(bits, d)[1]
    z1(x), z2(x)
end
seed = encode.([-10 + rand()*20 for _ =1:100] ,d)
res = nsga(500, 200, ()->rand(Bool, d.nbbitstotal), z, 5, seed = seed)

@test maximum(x -> x.y[1], res) >= 3.99
@test minimum(x -> x.y[1], res) <= 1e-4
@test maximum(x -> x.y[2], res) >= 3.99
@test minimum(x -> x.y[2], res) <= 1e-4


text:
module TurnTable
using Printf
using PyCall
using Libaudio

const PySerial = PyCall.PyNULL()
const PySerialListPorts = PyCall.PyNULL()
const SerialString = String

struct SerialPort <: IO
    port::SerialString
    baudrate::Int
    bytesize::Int
    parity::SerialString
    stopbits::Int
    timeout
    xonxoff::Bool
    rtscts::Bool
    dsrdtr::Bool
    python_ptr::PyObject
end

function __init__()
    copy!(PySerial, pyimport_conda("serial", "pyserial"))
    copy!(PySerialListPorts, pyimport("serial.tools.list_ports"))
end

function serialport(port, baudrate)
    py_ptr = PySerial[:Serial](port, baudrate)
    SerialPort(port,
               baudrate,
               py_ptr[:bytesize],
               py_ptr[:parity],
               py_ptr[:stopbits],
               py_ptr[:timeout],
               py_ptr[:xonxoff],
               py_ptr[:rtscts],
               py_ptr[:dsrdtr], py_ptr)
end

@static if VERSION >= v"0.5"
   function (::Type{SerialPort})(port, baudrate)
       serialport(port, baudrate)
   end
end
# serialport.python_ptr[:open]()

"""
List available serialports on the system.
"""
function listserialports()
    @static if Sys.isunix()
        ports = readdir("/dev/")
        f = is_apple() ? _valid_darwin_port : _valid_linux_port
        filter!(f, ports)
        return [string("/dev/", port) for port in ports]
    end
    @static if Sys.iswindows()
        [i[1] for i in collect(PySerialListPorts[:comports]())]
    end
end

function setorigin(id; baudrate = 19200)
    s = SerialPort(id, baudrate)
    s.python_ptr[:write](b"Get BaudRate\r")
    y = s.python_ptr[:read](6)

    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(root, :blink, Libaudio.nows() * " | TurnTable.setorigin: baudrate = $(y) bps")

    # disable analog input to prevent noise input (must)
    s.python_ptr[:write](b"Set AnalogInput OFF \r")
    y = s.python_ptr[:read](3)
    s.python_ptr[:write](b"Set PulseInput OFF \r")
    y = s.python_ptr[:read](3)
    s.python_ptr[:write](b"Set Torque 70.0 \r")
    y = s.python_ptr[:read](3)
    s.python_ptr[:write](b"Set SmartTorque ON \r")
    y = s.python_ptr[:read](3)
    s.python_ptr[:write](b"Set Velocity 2.00 \r")
    y = s.python_ptr[:read](3)
    s.python_ptr[:write](b"Set Origin \r")
    y = s.python_ptr[:read](3)
    s.python_ptr[:close]()
    Libaudio.printl(root, :blink, Libaudio.nows() * " | TurnTable.setorigin: set current position as origin")
    return y
end

function rotate(id, degree; direction="CCW", baudrate=19200)
    d = @sprintf("%3.1f", degree)
    for i = 1:5-length(d)
        d = "0" * d
    end
    s = SerialPort(id, baudrate)
    s.python_ptr[:write](Vector{UInt8}("GoTo " * direction * " +" * d * " \r"))
    y = s.python_ptr[:read](3)
    s.python_ptr[:close]()

    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(root, :blink, Libaudio.nows() * " | TurnTable.rotate: moved to degree $(degree) $(direction)")
    return y
end

end # module


text:
module TestStationary

using Test
using LinearAlgebra
# using LinearAlgebraExtensions: LowRank

using CovarianceFunctions
using CovarianceFunctions: MercerKernel, StationaryKernel, isstationary, isisotropic,
                    Constant, EQ, RQ, Exp, γExp, Delta, Cosine, MaternP, Matern#, SM
using CovarianceFunctions: iscov, enorm, Normed
using ForwardDiff: derivative

# isotropic kernels
const iso_k_strings = ["Exponentiated Quadratic", "Exponential", "δ",
            "Constant", "Rational Quadratic",
            "γ-Exponential", "Matern"]
# all stationary kernels
const k_strings = vcat(iso_k_strings, "Cosine")
                #, "Spectral Mixture"]

r = 2*rand()
const iso_k_arr = [EQ(), Exp(), Delta(),
        Constant(r), RQ(r),
        γExp(r), Matern(r)]
const k_arr = vcat(iso_k_arr, Cosine(r))
        #, SM]

const tol = 1e-12

# TODO:
# test higher input dimensions
# test type stability

@testset "basic properties" begin
    T = Float64
    n = 16
    x = randn(T, n) # IDEA: include higher-dimensional test cases
    Σ = zeros(T, (n, n))
    r² = sum(abs2, x[1] - x[2])
    for k in iso_k_arr
        @test k(x[1], x[2]) ≈ k(r²)
    end
    for (k, k_str) in zip(k_arr, k_strings)
        @testset "$k_str" begin
            Σ .= k.(x, permutedims(x))
            @test iscov(Σ, tol)
            @test isstationary(k, x)
        end
    end

    @testset "Matern" begin
        # testing different inputs
        for p = 0:8
            k = MaternP(p)
            Σ .= k.(x, permutedims(x))
            @test iscov(Σ, tol)
            @test isstationary(k, x)
        end

        # testing against control implementation
        # r² = x.^2
        r² = 10.0.^(1:16) * eps()
        p = 0
        k = MaternP(p)
        @test k.(r²) ≈ MaternP.(r², p)
        for p = 2:3
            k = MaternP(p)
            @test k.(0) ≈ MaternP.(0, p) # evaluation at zero
            @test k.(r²) ≈ MaternP.(r², p) # for small numbers up to O(1)
            # testing differentiability
            # k_naive = z->MaternP(z, p)
            k_naive = z->Matern(z, p + 1/2)
            dk_taylor = derivative.((k,), r²)
            dk_naive = derivative.((k_naive,), r²)
            #
            @test isapprox(dk_taylor, dk_naive, atol = 1e-6)

            # second order problem: naïve implementation is not as accurate anymore
            dk_taylor = derivative.((z->derivative(k, z),), r²)
            dk_naive = derivative.((z->derivative(k_naive, z),), r²)
            @test isapprox(dk_taylor, dk_naive, atol = 1e-5)
        end

        # testing constructor for invalid inputs
        @test_throws DomainError MaternP(-1)

        # can be constructed via Matern constructor with appropriate ν
        # we might not want to do that if we want to optimize ν
        for p = 0:4
            @test MaternP(Matern(p+1/2)) isa MaternP
        end

        # k = Matern(1.5)
        # @test derivative(k, 0) ≈ ... # this needs to be fixed with taylor expansion
    end

end

@testset "multidimensional inputs" begin
    n = 16
    Σ = zeros(n, n)
    for (k, k_str) in zip(k_arr, k_strings)
        if typeof(k) <: CovarianceFunctions.IsotropicKernel
            @testset "$k_str" begin
                for d = 2:3
                    x = [randn(d) for _ in 1:n]
                    K = k.(x, permutedims(x))
                    Σ = k.(x, permutedims(x))
                    @test iscov(Σ, tol)
                    @test isstationary(k, x)
                end
            end
        end
    end
end

@testset "kernel modifications" begin

    using CovarianceFunctions: Lengthscale
    l = exp(randn())
    for d = 1:3
        r = randn(d)
        for (k, k_str) in zip(iso_k_arr, iso_k_strings)
            kl = Lengthscale(k, l)
            @test kl(r) ≈ k(sum(abs2, r)/l^2)
            @test typeof(kl) <: CovarianceFunctions.IsotropicKernel
            @test isstationary(kl)
        end
    end

    ########### Testing ARD
    d = 2
    k = CovarianceFunctions.EQ()
    l = ones(d)
    kl = CovarianceFunctions.ARD(k, l)
    x = randn(d)
    y = randn(d)

    @test kl(x, y) ≈ k(x, y)

    @test sum(abs2, x-y) ≈ kl.n²(x-y) # this tests Metrics
    c = 2
    l = c^2*l
    kl = CovarianceFunctions.ARD(k, l)
    @test 1/c^2*sum(abs2, x-y) ≈ kl.n²(x-y)

    l = exp.(randn(d))
    w = @. sqrt(1/l)
    kl = CovarianceFunctions.ARD(k, l)
    @test sum(abs2, (x-y).*w) ≈ kl.n²(x-y)
    xl = x .* w
    yl = y .* w
    @test kl(x, y) ≈ k(xl, yl)

    d = 3
    U = randn(1, d)
    n²(x) = CovarianceFunctions.enorm2(U'U, x)
    kl = CovarianceFunctions.Normed(k, n²)
    x = randn(d)
    y = randn(d)
    @test kl(x, y) ≈ k(U*x, U*y)
    @test kl(x, y) ≈ k((x-y)'*(U'U)*(x-y))

    S = U'U # LowRank(U')
    kSLR = CovarianceFunctions.Energetic(k, S)
    @test kSLR(x, y) ≈ kl(x, y)

    # periodic CovarianceFunctions
    k = CovarianceFunctions.EQ()
    p = CovarianceFunctions.Periodic(k)
    x, y = randn(2)
    @test p(x, y) isa Real
    @test p(1+x) ≈ p(x) # 1 periodic
    @test p(x, 1+y) ≈ p(x, y)  # 1 periodic
    @test p(x+1, y-1) ≈ p(x, y)  # 1 periodic
    @test p(x+7, y-3) ≈ p(x, y)  # 1 periodic

    # TODO: write benchmarks
    # println(size(U*x))
    # println(size(x))
    # L = randn(d, 16)
    # @time k(L'x, L'y)
    # @time k(U*x, U*y)
    # @time kl(x, y)
    # using MyFactorizations: SymmetricLowRank
    # S = SymmetricLowRank(U)
    # kSLR = CovarianceFunctions.Subspace(k, S)

end # CovarianceFunctions modifications

end # TestStationary

# @testset "Rational Quadratic" begin
#     α = exp(randn())
#     k = RQ(α)
#     Σ .= k.(x, permutedims(x))
#     @test iscov(Σ)
#     @test isstationary(k, x)
#     # test constructor for negative input
# end
#
# @testset "Exponential" begin
#     k = Exp()
#     Σ .= k.(x, permutedims(x))
#     @test iscov(Σ)
#     @test isstationary(k, x)
# end
#
# @testset "γ-Exponential CovarianceFunctions" begin
#     γ = 2rand()
#     k = ΓExp(γ)
#     Σ .= k.(x, permutedims(x))
#     @test iscov(Σ)
#     @test isstationary(k, x)
#     # test constructor for γ outside of [0,2]
# end
#
# @testset "δ CovarianceFunctions" begin
#     k = Delta()
#     Σ .= k.(x, permutedims(x))
#     @test iscov(Σ)
#     @test isstationary(k, x)
# end


text:
function test_jump_lp(solver)
    @testset "Linear Programming" begin
        @testset "LP1" begin
            # simple 2 variable, 1 constraint problem
            # min -x
            # st   x + y <= 1   (x + y - 1 ∈ Nonpositives)
            #       x, y >= 0   (x, y ∈ Nonnegatives)

            m = Model(solver=solver)
            @variable(m, x >= 0.0)
            @variable(m, y >= 0.0)
            @objective(m, Min, -x)

            @constraint(m, x + y <= 1)

            JuMP.solve(m)

            @test JuMP.isattached(m)
            @test JuMP.hasvariableresult(m)

            @test JuMP.terminationstatus(m) == MOI.Success
            @test JuMP.primalstatus(m) == MOI.FeasiblePoint

            @test JuMP.resultvalue(x) ≈ 1.0 atol=1e-6
            @test JuMP.resultvalue(y) ≈ 0.0 atol=1e-6
            @test JuMP.objectivevalue(m) ≈ -1.0 atol=1e-6

        end
    end
end


text:
# gradient 计算函数f(x)微分 一元微分直接在gradient函数后取第一个元组数据
function df1()
    f(x) = 3x^2 + 2x + 1
    df(x) = gradient(f, x)[1]
    df(2)
end

# 如果f(x)的参数是矩阵 gradient调用的函数需要包含在sum()函数里面 f(x)函数运算符前要加.作为广播运算符 计算结果为一个矩阵
function df_arg_matrix()
    f(x) = 3 .* x .^2 + 2 .*x .+ 1
    sf(x) = sum(f(x))
    dsf(x) = gradient(sf, x)[1]
    dsf([2,1])
end

# gradient 计算多元微分 分别计算各元的偏微分
function df_m()
    g(x, y) = (x - y)^2
    dg(x, y) = gradient(g, x, y)
    dg(2,1)
end

# 简单模型
function simple_models()
    W = rand(2, 5)
    b = rand(2)

    predict(x) = W*x .+ b

    function loss(x, y)
        ŷ = predict(x)
        sum((y .- ŷ).^2)
    end

    x, y = rand(5), rand(2)
    bf = loss(x, y)
    gs = gradient(() -> loss(x, y), params(W, b))
    W̄ = gs[W]
    W .-= 0.1 .* W̄
    af = loss(x, y)
    bf, af
end

# 建立分成模型
function build_layers()
    W1 = rand(3, 5)
    b1 = rand(3)
    layer1(x) = W1 * x .+ b1

    W2 = rand(2, 3)
    b2 = rand(2)
    layer2(x) = W2 * x .+ b2

    model(x) = layer2(σ.(layer1(x)))

    model(rand(5))
end

# 建立分成模型 使用linear函数优化
function build_layers_o1()
    function linear(in, out)
      W = randn(out, in)
      b = randn(out)
      x -> W * x .+ b
    end

    linear1 = linear(5, 3) # we can access linear1.W etc
    linear2 = linear(3, 2)

    model(x) = linear2(σ.(linear1(x)))

    model(rand(5))
end

# 建立分成模型 使用Dense
function build_layers_dense()
    model = Chain(
      Dense(10, 5, σ),
      Dense(5, 2),
      softmax)

    model(rand(10))
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module SK

using Random, LinearAlgebra
using ExtractMacro
using ..Interface
using ..Common
using ..QT

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphSK, GraphSKNormal

import ..Interface: energy, delta_energy, update_cache!, neighbors

function gen_J(N::Integer)
    J = BitVector[bitrand(N) for i = 1:N]
    for i = 1:N
        J[i][i] = 0
        for j = (i+1):N
            J[j][i] = J[i][j]
        end
    end
    return J
end

struct GraphSK <: SimpleGraph{Float64}
    N::Int
    sN::Float64
    J::Vector{BitVector}
    #tmps::BitVector
    cache::LocalFields{Int}
    function GraphSK(J::Vector{BitVector}; check::Bool = true)
        N = length(J)
        if check
            all(Jx->length(Jx) == N, J) || throw(ArgumentError("invalid J inner length, expected $N, given: $(unique(map(length,J)))"))
            for i = 1:N
                J[i][i] == 0 || throw(ArgumentError("diagonal entries of J must be 0, found: J[$i][$i] = $(J[i][i])"))
                for j = (i+1):N
                    J[i][j] == J[j][i] || throw(ArgumentError("J must be symmetric, found: J[$i][$j] = $(J[i][j]), J[$j][$i] = $(J[j][i])"))
                end
            end
        end
        #tmps = BitArray(undef, N)
        cache = LocalFields{Int}(N)
        return new(N, √N, J, cache)
    end
end

"""
    GraphSK(N::Integer) <: SimpleGraph{Float64}

A `SimpleGraph` implementing a Sherrington-Kirkpatrick fully-connected Ising model
with `N` spins and random binary interactions (\$J ∈ \\{-1/√N,1/√N\\}\$) and no
external fields.

Same as [`GraphSKNormal`](@ref), but with binary interactions.
"""
GraphSK(N::Integer) = GraphSK(gen_J(N), check=false)

function energy(X::GraphSK, C::Config)
    @assert X.N == C.N
    @extract C : s
    @extract X : N sN J cache
    @extract cache : lfields lfields_last
    tmps = BitArray(undef, N)
    n = -2 * sum(s)
    for i = 1:N
        Ji = J[i]
        sc = sum(map!(⊻, tmps, Ji, s))
        si = s[i]
        lf = -(2si-1) * (N-1 - 2sc)
        lfields[i] = 2 * (-lf + 2si)
        n += lf
    end
    @assert n % 2 == 0
    n ÷= 2
    cache.move_last = 0
    fill!(lfields_last, 0)

    # altn = 0
    # for i = 1:N
    #     Ji = J[i]
    #     for j = 1:N
    #         j == i && continue
    #         altn -= (2Ji[j] - 1) * (2s[j] - 1) * (2s[i] - 1)
    #     end
    # end
    # altn ÷= 2
    # @assert n == altn

    return n / sN
end

function update_cache!(X::GraphSK, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : J cache

    @extract cache : lfields lfields_last move_last

    if move_last == move
        cache.lfields, cache.lfields_last = cache.lfields_last, cache.lfields
        return
    end

    @inbounds begin
        Ji = J[move]
        si = s[move]
        lfm = lfields[move]
        @simd for j = 1:N
            # note: we don't check move ≠ j to avoid branching
            Jσij = si ⊻ s[j] ⊻ Ji[j]
            lfj = lfields[j]
            lfields_last[j] = lfj
            lfields[j] = lfj + 8 * Jσij - 4
        end
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache.move_last = move

    # lfields_bk = copy(lfields)
    # lfields_last_bk = copy(lfields_last)
    # energy(X, C)
    # @assert lfields_bk == lfields
    # copy!(lfields_last, lfields_last_bk)
    # cache.move_last = move

    return
end

function delta_energy(X::GraphSK, C::Config, move::Int)
    @extract X : sN cache
    @extract cache : lfields

    @inbounds Δ = lfields[move] / sN
    return Δ

    # @extract X : N J tmps
    # @extract C : s
    # @assert N == length(s)
    # @assert 1 ≤ move ≤ N

    # Ji = J[move]
    # si = s[move]
    # sc = sum(map!(⊻, tmps, Ji, s)) - si
    # @assert Δ == 2 * (2si-1) * (N-1 - 2sc)
    # return Δ
end

# function check_delta(X::GraphSK, C::Config, move::Int)
#     @extract C : s
#     delta = delta_energy(X, s, move)
#     e0 = energy(X, s)
#     s[move] ⊻= 1
#     e1 = energy(X, s)
#     s[move] ⊻= 1
#
#     (e1-e0) == delta || (@show e1,e0,delta,e1-e0; error())
# end

neighbors(X::GraphSK, i::Int) = AllButOne(X.N, i)

#####

function gen_J_gauss(N::Integer)
    J = Vec[rmul!(randn(N), 1/√N) for i = 1:N]
    for i = 1:N
        J[i][i] = 0
        for j = (i+1):N
            J[j][i] = J[i][j]
        end
    end
    return J
end

struct GraphSKNormal <: SimpleGraph{Float64}
    N::Int
    J::Vec2
    cache::LocalFields{Float64}
    function GraphSKNormal(J::Vec2; check::Bool = true)
        N = length(J)
        if check
            all(Jx->length(Jx) == N, J) || throw(ArgumentError("invalid J inner length, expected $N, given: $(unique(map(length,J)))"))
            for i = 1:N
                J[i][i] == 0 || throw(ArgumentError("diagonal entries of J must be 0, found: J[$i][$i] = $(J[i][i])"))
                for j = (i+1):N
                    J[i][j] == J[j][i] || throw(ArgumentError("J must be symmetric, found: J[$i][$j] = $(J[i][j]), J[$j][$i] = $(J[j][i])"))
                end
            end
        end
        cache = LocalFields{Float64}(N)
        return new(N, J, cache)
    end
end

"""
    GraphSK(N::Integer) <: SimpleGraph{Float64}

A `SimpleGraph` implementing a Sherrington-Kirkpatrick fully-connected Ising model
with `N` spins and random interactions extracted from a normal distribution with
zero mean and \$1/N\$ variance, and no external fields.

Same as [`GraphSK`](@ref), but with Gaussian interactions.
"""
GraphSKNormal(N::Integer) = GraphSKNormal(gen_J_gauss(N), check=false)

function energy(X::GraphSKNormal, C::Config)
    @assert X.N == C.N
    @extract C : s
    @extract X : N J cache
    @extract cache : lfields lfields_last

    n = 0.0
    for i = 1:N
        Ji = J[i]
        si = s[i]
        lf = 0.0
        for j = 1:N
            # j == i && continue # avoid branching
            #n -= (2s[j] - 1) * (2s[i] - 1) * Ji[j]
            lf += (1 - 2 * (si ⊻ s[j])) * Ji[j]
        end
        lfields[i] = 2lf
        n -= lf
    end
    n /= 2

    cache.move_last = 0
    fill!(lfields_last, 0)

    return n
end

function update_cache!(X::GraphSKNormal, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : J cache

    @extract cache : lfields lfields_last move_last

    if move_last == move
        cache.lfields, cache.lfields_last = cache.lfields_last, cache.lfields
        return
    end

    @inbounds begin
        Ji = J[move]
        si = s[move]
        lfm = lfields[move]
        @simd for j = 1:N
            # note: we don't check move ≠ j to avoid branching
            Jσij = (1 - 2 * (si ⊻ s[j])) * Ji[j]
            lfj = lfields[j]
            lfields_last[j] = lfj
            lfields[j] = lfj + 4 * Jσij
        end
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache.move_last = move

    # lfields_bk = copy(lfields)
    # lfields_last_bk = copy(lfields_last)
    # energy(X, C)
    # @assert maximum(abs.(lfields_bk - lfields)) < 1e-10
    # copy!(lfields_last, lfields_last_bk)
    # cache.move_last = move

    return
end

function delta_energy(X::GraphSKNormal, C::Config, move::Int)
    @extract X : cache
    @extract cache : lfields

    @inbounds Δ = lfields[move]
    return Δ
end

function check_delta(X::GraphSKNormal, C::Config, move::Int)
    @extract C : s
    delta = delta_energy(X, C, move)
    s[move] = s[move] ⊻ 1
    e1 = energy(X, C)
    s[move] = s[move] ⊻ 1
    e0 = energy(X, C)

    abs((e1-e0) - delta) < 1e-10 || (@show e1,e0,delta,e1-e0; error())
end

neighbors(X::GraphSKNormal, i::Int) = AllButOne(X.N, i)

end


text:
# # Arbitrary floating-point types
#
# Since DFTK is completely generic in the floating-point type
# in its routines, there is no reason to perform the computation
# using double-precision arithmetic (i.e.`Float64`).
# Other floating-point types such as `Float32` (single precision)
# are readily supported as well.
# On top of that we already reported[^HLC2020] calculations
# in DFTK using elevated precision
# from [DoubleFloats.jl](https://github.com/JuliaMath/DoubleFloats.jl)
# or interval arithmetic
# using [IntervalArithmetic.jl](https://github.com/JuliaIntervals/IntervalArithmetic.jl).
# In this example, however, we will concentrate on single-precision
# computations with `Float32`.
#
# The setup of such a reduced-precision calculation is basically identical
# to the regular case, since Julia automatically compiles all routines
# of DFTK at the precision, which is used for the lattice vectors.
# Apart from setting up the model with an explicit cast of the lattice
# vectors to `Float32`, there is thus no change in user code required:
#
# [^HLC2020]:
#     M. F. Herbst, A. Levitt, E. Cancès.
#     *A posteriori error estimation for the non-self-consistent Kohn-Sham equations*
#     [ArXiv 2004.13549](https://arxiv.org/abs/2004.13549)

using DFTK

## Setup silicon lattice
a = 10.263141334305942  # lattice constant in Bohr
lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]

## Cast to Float32, setup model and basis
model = model_DFT(Array{Float32}(lattice), atoms, positions, [:lda_x, :lda_c_vwn])
basis = PlaneWaveBasis(model, Ecut=7, kgrid=[4, 4, 4])

## Run the SCF
scfres = self_consistent_field(basis, tol=1e-4);

# To check the calculation has really run in Float32,
# we check the energies and density are expressed in this floating-point type:

scfres.energies
#-
eltype(scfres.energies.total)
#-
eltype(scfres.ρ)

#
# !!! note "Generic linear algebra routines"
#     For more unusual floating-point types (like IntervalArithmetic or DoubleFloats),
#     which are not directly supported in the standard `LinearAlgebra` library of Julia
#     one additional step is required: One needs to explicitly enable the generic versions
#     of standard linear-algebra operations like `cholesky` or `qr`, which are needed
#     inside DFTK by loading the `GenericLinearAlgebra` package in the user script
#     (i.e. just add ad `using GenericLinearAlgebra` next to your `using DFTK` call).
#


text:
"""
Performance statistics of greedy_multi_path for 1 random end-user pair
over a range of graph sizes and for different numbers of max_path (1,2,3) (ie. the
maximum number of paths that can be purified by a given end-user pair.)
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

datafile = "data/maxpaths"
generate_new_data = false
if generate_new_data == true

    # Params
    # 1000
    num_trials = 10000::Int64
    min_size = 5::Int64
    max_size = 30::Int64

    @assert min_size < max_size

    perf_data1 = []
    perf_data2 = []
    perf_data3 = []
    perf_data4 = []

    path_data4 = []
    path_data4_err = []

    size_list = collect(min_size:1:max_size)
    for (j, data_array) in enumerate([perf_data1, perf_data2, perf_data3, perf_data4])
        println("Collecting for max_paths: $j")
        for i in size_list
            println("Collecting for gridsize: $i")
            # Generate ixi graph:
            net = GridNetwork(i, i)

            # Collect performance statistics
            performance, dummy, pat, pat_e = net_performance(net, num_trials, 1, max_paths=j)
            push!(data_array, performance)

            if j == 4
                push!(path_data4, pat)
                push!(path_data4_err, pat_e)
            end
        end
    end

    # Save data
    d = Dict{Symbol, Any}()
    @pack! d = num_trials, min_size, max_size, perf_data1, perf_data2, perf_data3, perf_data4, path_data4, path_data4_err
    save("$datafile.jld", "data", d)

else
    # Load data
    if !isfile("$datafile.jld")
        error("$datafile.jld not found in working directory")
    end
    d = load("$datafile.jld")["data"]
    @unpack num_trials, min_size, max_size, perf_data1, perf_data2, perf_data3, perf_data4, path_data4, path_data4_err = d
end

# Get values for x axis
x = collect(min_size:1:max_size)

# Extract data from performance data
loss_arr1 = collect(map(x->x["loss"], perf_data1))
z_arr1 = collect(map(x->x["Z"], perf_data1))
loss_arr2 = collect(map(x->x["loss"], perf_data2))
z_arr2 = collect(map(x->x["Z"], perf_data2))
loss_arr3 = collect(map(x->x["loss"], perf_data3))
z_arr3 = collect(map(x->x["Z"], perf_data3))
loss_arr4 = collect(map(x->x["loss"], perf_data4))
z_arr4 = collect(map(x->x["Z"], perf_data4))

##### Plot cost_maxpaths #####
plot1 = plot(x, loss_arr1, seriestype = :scatter, label=L"$\eta_1$", xlims=(0, max_size+11),
markershape = :utriangle, legend=:best, guidefontsize=14, tickfontsize=12, legendfontsize=8,
fontfamily="computer modern", markersize=5, color_palette = palette(:Spectral_4, 4))
plot!(x, loss_arr2, seriestype = :scatter, label=L"$\eta_2$", markershape=:utriangle,
markersize=5)

# [3:end]
plot!(x, loss_arr3, seriestype = :scatter, label=L"$\eta_3$", markershape=:utriangle,
markersize=5)
plot!(x, loss_arr4, seriestype = :scatter, label=L"$\eta_3$", markershape=:utriangle,
markersize=5)
plot!(x, z_arr1, seriestype = :scatter, label=L"$F_1$", markersize = 5)
plot!(x, z_arr2, seriestype = :scatter, label=L"$F_2$", markersize = 5)
plot!(x, z_arr3, seriestype = :scatter, label=L"$F_3$", markersize = 5)
plot!(x, z_arr4, seriestype = :scatter, label=L"$F_4$", markersize = 5)
yaxis!("Costs")

# Plot analytic function for average cost:
function ave_pathlength(n)
    return 2/3 * n
end

function ave_same(n)
    return (n+1)/3
end

function ave_dif(n)
    return 2*(1+n)/3
end

ave_e = []
ave_f = []
for size in min_size:max_size
    e = dB_to_P(ave_pathlength(size))
    f = dB_to_Z(ave_pathlength(size))
    push!(ave_e, e)
    push!(ave_f, f)
end

# plot!(x, ave_e, linewidth=2, label=L"$\textrm{ave } \eta$")
# plot!(x, ave_f, linewidth=2, label=L"$\textrm{ave } F$", linestyle=:dash)

# # Get data for naive 2 path purification
# naive_e2 = []
# naive_f2 = []
# for size in min_size:max_size
#     e = dB_to_P(ave_pathlength(size))
#     f = dB_to_Z(ave_pathlength(size))
#     p = 2 / (size + 1)
#     epur, fpur = QuNet.purify(e, e, f, f)
#     push!(naive_e2, epur)
#     push!(naive_f2, fpur)
# end
# plot!(x, naive_e2, linewidth=2, label=L"$\textrm{Naive E2}$", linestyle=:dot, color =:blue)
# plot!(x, naive_f2, linewidth=2, label=L"$\textrm{Naive F2}$", linestyle=:dot, color =:blue)

# Get analytic data for 2 path purification
# ave_e2 = []
# ave_f2 = []

# 2 path-purification
ave_e2_same = []
ave_f2_same = []
ave_e2_dif = []
ave_f2_dif = []
# 3-path purification
ave_e3_same = []
ave_f3_same = []
ave_e3_dif = []
ave_f3_dif = []
# 4-path purification
ave_e4_same = []
ave_f4_same = []
ave_e4_dif = []
ave_f4_dif = []

for size in min_size:max_size
    # Method 1
    # p = 2 / (size + 1)
    # e = dB_to_P(ave_pathlength(size))
    # eplus = dB_to_P(ave_pathlength(size) + 2)
    # f = dB_to_Z(ave_pathlength(size))
    # fplus = dB_to_Z(ave_pathlength(size) + 2)
    # e2, f2 = QuNet.purify(e, eplus, f, fplus) .* p .+ QuNet.purify(e, e, f, f) .* (1-p)
    # # end of method 1

    # Method 2:
    # e = dB_to_P(ave_pathlength(size))
    # f = dB_to_Z(ave_pathlength(size))
    # eplus = dB_to_P(ave_pathlength(size) + 2)
    # fplus = dB_to_Z(ave_pathlength(size) + 2)
    # p = 2 / (size + 1)
    # e2, f2 = QuNet.purify(e, eplus*p + e*(1-p), f, fplus*p + f*(1-p))
    # # end of method 2

    # # Method 3
    # # Probability that pair lies on the same row or column
    # p = 2 / (size + 1)
    # # Average costs when users lie on different rows/columns
    # e_dif = dB_to_P(ave_dif(size))
    # f_dif = dB_to_Z(ave_dif(size))
    # # Average costs when users lie on same row/column
    # e_same = dB_to_P(ave_same(size))
    # f_same = dB_to_Z(ave_same(size))
    # # Average cost when users lie on same row/column plus two average edges
    # e_same2 = dB_to_P(ave_same(size) + 2)
    # f_same2 = dB_to_Z(ave_same(size) + 2)
    # # Average purification cost when users lie on different rows/columns
    # e_dif, f_dif = QuNet.purify(e_dif, e_dif, f_dif, f_dif)
    # # Average purification cost when same row/column
    # e_same, f_same = QuNet.purify(e_same, e_same2, f_same, f_same2)
    # e2 = (1-p)*e_dif + p*e_same
    # f2 = (1-p)*f_dif + p*f_same

    # Method 4:
    # Average costs when users lie on different row/columns
    e_dif = dB_to_P(ave_dif(size))
    f_dif = dB_to_Z(ave_dif(size))
    # Average cost when users lie on same row/column plus four average edges
    e_dif4 = dB_to_P(ave_dif(size) + 4)
    f_dif4 = dB_to_Z(ave_dif(size) + 4)
    # Average costs when users lie on same row/column
    e_same = dB_to_P(ave_same(size))
    f_same = dB_to_Z(ave_same(size))
    # Average cost when users lie on same row/column plus two average edges
    e_same2 = dB_to_P(ave_same(size) + 2)
    f_same2 = dB_to_Z(ave_same(size) + 2)
    # "" + 8 average edges
    e_same8 = dB_to_P(ave_same(size) + 8)
    f_same8 = dB_to_Z(ave_same(size) + 8)

    # 2 Purification
    # Average purification cost when users lie on different rows/columns
    e2_dif, f2_dif = QuNet.purify(e_dif, e_dif, f_dif, f_dif)
    # Average purification cost when same row/column
    e2_same, f2_same = QuNet.purify(e_same, e_same2, f_same, f_same2)

    push!(ave_e2_same, e2_same)
    push!(ave_f2_same, f2_same)
    push!(ave_e2_dif, e2_dif)
    push!(ave_f2_dif, f2_dif)

    # 3 Purification
    e3_dif, f3_dif = QuNet.purify(e_dif4, e2_dif, f_dif4, f2_dif)
    e3_same, f3_same = QuNet.purify(e2_same, e_same2, f2_same, f_same2)

    push!(ave_e3_dif, e3_dif)
    push!(ave_f3_dif, f3_dif)
    push!(ave_e3_same, e3_same)
    push!(ave_f3_same, f3_same)

    # 4 Purification
    e4_dif, f4_dif = QuNet.purify(e_dif4, e3_dif, f_dif4, f3_dif)
    e4_same, f4_same = QuNet.purify(e3_same, e_same8, f3_same, f_same8)

    push!(ave_e4_dif, e4_dif)
    push!(ave_f4_dif, f4_dif)
    push!(ave_e4_same, e4_same)
    push!(ave_f4_same, f4_same)
end

# plot!(x, ave_e2_best, linewidth=2, label=L"$\textrm{Best case } E_2$", linestyle=:dot, color =:blue)
# plot!(x, ave_f2_best, linewidth=2, label=L"$\textrm{Best case } F_2$", linestyle=:dot, color =:blue)
# plot!(x, ave_e2_worst, linewidth=2, label=L"$\textrm{Worst case } E_2$", color =:blue)
# plot!(x, ave_f2_worst, linewidth=2, label=L"$\textrm{Worst case } F_2$", linestyle=:dot, color =:blue)
# # Plot 3 path purification
# plot!(x, ave_e3_worst, linewidth=2, label=L"$\textrm{Worst case } E_3$", color =:green)
# plot!(x, ave_f3_worst, linewidth=2, label=L"$\textrm{Worst case } F_3$", linestyle=:dot, color =:green)
# # Plot 4 path purification
# plot!(x, ave_e4_worst, linewidth=2, label=L"$\textrm{Worst case } E_4$", color =:purple)
# plot!(x, ave_f4_worst, linewidth=2, label=L"$\textrm{Worst case } F_4$", linestyle=:dot, color =:purple)

# plot!(x, ave_e2, linewidth=2, label=L"$\textrm{ave } F_2$", linestyle=:dot, color =:black)
# plot!(x, ave_f2, linewidth=2, label=L"$\textrm{ave } F_2$", linestyle=:dot, color =:black)

# xaxis!(L"$\textrm{Grid Size}$")
savefig("plots/cost_maxpaths.png")
savefig("plots/cost_maxpaths.pdf")

##### plot elog_cost_maxpaths #####
plot2 = plot(x, loss_arr1, seriestype = :scatter, label=L"$\eta_1$", xlims=(0, max_size+11), color_palette = palette(:Spectral_4, 4),
markershape = :utriangle, legend=:best, guidefontsize=14, tickfontsize=12, legendfontsize=8,
fontfamily="computer modern", markersize=5, yaxis=:log)
plot!(x, loss_arr2, seriestype = :scatter, label=L"$\eta_2$", markershape=:utriangle,
markersize=5)
plot!(x, loss_arr3, seriestype = :scatter, label=L"$\eta_3$", markershape=:utriangle,
markersize=5)
plot!(x, loss_arr4, seriestype = :scatter, label=L"$\eta_3$", markershape=:utriangle,
markersize=5)
yaxis!("Costs")
plot!(x, ave_e, linewidth=2, label=L"$\textrm{Analytic } \eta_1$")
plot!(x, ave_e2_dif, fillrange = ave_e2_same, alpha=0.3, linewidth=2, label=L"$\eta_2 \textrm{ Est.}$")
plot!(x, ave_e3_dif, fillrange = ave_e3_same, alpha=0.3, linewidth=2, label=L"$\eta_3 \textrm{ Est.}$")
plot!(x, ave_e4_dif, fillrange = ave_e4_same, alpha=0.3, linewidth=2, label=L"$\eta_4 \textrm{ Est.}$")

savefig("plots/elog_cost_maxpaths.png")
savefig("plots/elog_cost_maxpaths.pdf")

##### plot flog_cost_maxpaths #####
plot3 = plot(x, z_arr1 .- .5, seriestype = :scatter, label=L"$F_1$", xlims=(0, max_size+11), color_palette = palette(:Spectral_4, 4),
legend=:best, guidefontsize=14, tickfontsize=12, legendfontsize=8,
fontfamily="computer modern", markersize=5, yaxis=:log)
plot!(x, z_arr2 .- .5, seriestype = :scatter, label=L"$F_2$", markersize=5)
plot!(x, z_arr3 .- .5, seriestype = :scatter, label=L"$F_3$", markersize=5)
plot!(x, z_arr4 .- .5, seriestype = :scatter, label=L"$F_4$", markersize=5)
plot!(x, ave_f .- .5, linewidth=2, label=L"$\textrm{Analytic } F$")
plot!(x, ave_f2_dif .- .5, fillrange = ave_f2_same .- .5, alpha=0.3, linewidth=2, label=L"$F_2 \textrm{ Est.} $")
plot!(x, ave_f3_dif .- .5, fillrange = ave_f3_same .- .5, alpha=0.3, linewidth=2, label=L"$F_3 \textrm{ Est.}$")
plot!(x, ave_f4_dif .- .5, fillrange = ave_f4_same .- .5, alpha=0.3, linewidth=2, label=L"$F_4 \textrm{ Est.}$")

yaxis!("Costs")

savefig("plots/flog_cost_maxpaths.png")
savefig("plots/flog_cost_maxpaths.pdf")

# Average of purifications not purification of averages

# Method outline
# Given d, the gridsize, we consider a vertex on an infinite grid lattice
# and the set of all other vertices that lie within d steps.
# This set can be subdivided into the set of all points of distance one, two, etc from the central vertex
# To consider the average purification, we get the costs associated with pathlength 1
# And perform 2, 3, and 4 purification. We assume the worst case path lengths of
# L & L & 4+L & 4+L
# With these costs we calculate the mean, taking into account weights which result from
# The number of paths present (4*d paths)
# ave_e2 = []
# ave_f2 = []
# ave_e3 = []
# ave_f3 = []
# ave_e4 = []
# ave_f4 = []
#
# for max in min_size:max_size
#     e2_arr = []
#     f2_arr = []
#     e3_arr = []
#     f3_arr = []
#     e4_arr = []
#     f4_arr = []
#
#     for size in 1:max
#         num_paths = 4 * (size-1)
#         e = dB_to_P(Float64(size))
#         f = dB_to_Z(Float64(size))
#         ep4 = dB_to_P(Float64(size+4))
#         fp4 = dB_to_Z(Float64(size+4))
#
#         e2, f2 = QuNet.purify(e, e, f, f)
#         e3, f3 = QuNet.purify(e2, ep4, f2, fp4)
#         e4, f4 = QuNet.purify(e3, ep4, f3, fp4)
#
#         # Push once for each possible path
#         # This will balance out the averages later
#         for np in 1:num_paths
#             push!(e2_arr, e2)
#             push!(e3_arr, e3)
#             push!(e4_arr, e4)
#             push!(f2_arr, f2)
#             push!(f3_arr, f3)
#             push!(f4_arr, f4)
#         end
#     end
#     # Take averages of each
#     push!(ave_e2, mean(e2_arr))
#     push!(ave_e3, mean(e3_arr))
#     push!(ave_e4, mean(e4_arr))
#     push!(ave_f2, mean(f2_arr))
#     push!(ave_f3, mean(f3_arr))
#     push!(ave_f4, mean(f4_arr))
# end
#
# # Plot average of purifications
# plot(x, loss_arr1, seriestype = :scatter, label=L"$\eta_1$", xlims=(0, max_size), color = :red,
# legend=:best, legendfontsize = 7, xguidefontsize = 10, tickfontsize = 10)
# plot!(x, z_arr1, seriestype = :scatter, label=L"$F_1$", color =:red)
# plot!(x, loss_arr2, seriestype = :scatter, label=L"$\eta_2$", color =:blue)
# plot!(x, z_arr2, seriestype = :scatter, label=L"$F_2$", color =:blue)
# plot!(x, loss_arr3, seriestype = :scatter, label=L"$\eta_3$", color =:green)
# plot!(x, z_arr3, seriestype = :scatter, label=L"$F_3$", linestyle=:dash, color =:green)
# plot!(x, loss_arr4, seriestype = :scatter, label=L"$\eta_3$", color =:purple)
# plot!(x, z_arr4, seriestype = :scatter, label=L"$F_4$", linestyle=:dash, color =:purple)
#
# plot!(x, ave_e2, linewidth=2, label=L"$\textrm{Average } E_2$", linestyle=:dot, color =:blue)
# plot!(x, ave_f2, linewidth=2, label=L"$\textrm{Average } F_2$", linestyle=:dot, color =:blue)
# plot!(x, ave_e3, linewidth=2, label=L"$\textrm{Average } E_3$", linestyle=:dot, color =:black)
# plot!(x, ave_f3, linewidth=2, label=L"$\textrm{Average } F_3$", linestyle=:dot, color =:black)
# plot!(x, ave_e4, linewidth=2, label=L"$\textrm{Average } E_4$", linestyle=:dot, color =:purple)
# plot!(x, ave_f4, linewidth=2, label=L"$\textrm{Average } F_4$", linestyle=:dot, color =:purple)

# Average 2: Nathan's Method:

dif_e2 = []
dif_f2 = []
dif_e3 = []
dif_f3 = []
dif_e4 = []
dif_f4 = []

same_e2 = []
same_f2 = []
same_e3 = []
same_f3 = []
same_e4 = []
same_f4 = []

for size in min_size:max_size
    L = 2/3 * size
    e = dB_to_P(Float64(L))
    f = dB_to_Z(Float64(L))
    ep2 = dB_to_P(Float64(L+2))
    fp2 = dB_to_Z(Float64(L+2))
    ep4 = dB_to_P(Float64(L+4))
    fp4 = dB_to_Z(Float64(L+4))
    ep8 = dB_to_P(Float64(L+8))
    fp8 = dB_to_Z(Float64(L+8))

    e2, f2 = QuNet.purify(e, e, f, f)
    e3, f3 = QuNet.purify(e2, ep4, f2, fp4)
    e4, f4 = QuNet.purify(e3, ep4, f3, fp4)

    se2, sf2 = QuNet.purify(e, ep2, f, fp2)
    se3, sf3 = QuNet.purify(se2, ep2, sf2, fp2)
    se4, sf4 = QuNet.purify(se3, ep8, sf3, fp8)

    push!(dif_e2, e2)
    push!(dif_e3, e3)
    push!(dif_e4, e4)
    push!(dif_f2, f2)
    push!(dif_f3, f3)
    push!(dif_f4, f4)

    push!(same_e2, se2)
    push!(same_e3, se3)
    push!(same_e4, se4)
    push!(same_f2, sf2)
    push!(same_f3, sf3)
    push!(same_f4, sf4)
end

# Plot "testcost_maxpaths"
plot(x, loss_arr1, seriestype = :scatter, label=L"$\eta_1$", xlims=(0, max_size+10), color = :red,
legend=:best, legendfontsize = 7, xguidefontsize = 10, tickfontsize = 10, yaxis=:log)
plot!(x, z_arr1 .- .5, seriestype = :scatter, label=L"$F_1$", color =:red)
plot!(x, loss_arr2, seriestype = :scatter, label=L"$\eta_2$", color =:blue)
plot!(x, z_arr2 .- .5, seriestype = :scatter, label=L"$F_2$", color =:blue)
plot!(x, loss_arr3, seriestype = :scatter, label=L"$\eta_3$", color =:green)
plot!(x, z_arr3 .- .5, seriestype = :scatter, label=L"$F_3$", linestyle=:dash, color =:green)
plot!(x, loss_arr4, seriestype = :scatter, label=L"$\eta_3$", color =:purple)
plot!(x, z_arr4 .- .5, seriestype = :scatter, label=L"$F_4$", linestyle=:dash, color =:purple)

plot!(x, dif_e2, linewidth=2, label=L"$\textrm{Diff case } E_2$", linestyle=:dot, color =:blue)
plot!(x, dif_f2 .- .5, linewidth=2, label=L"$\textrm{Diff case } F_2$", linestyle=:dot, color =:blue)
plot!(x, dif_e3, linewidth=5, label=L"$\textrm{Diff case } E_3$", linestyle=:dot, color =:green)
plot!(x, dif_f3 .- .5, linewidth=2, label=L"$\textrm{Diff case } F_3$", linestyle=:dot, color =:green)
plot!(x, dif_e4, linewidth=2, label=L"$\textrm{Diff case } E_4$", linestyle=:dot, color =:purple)
plot!(x, dif_f4 .- .5, linewidth=2, label=L"$\textrm{Diff case } F_4$", linestyle=:dot, color =:purple)

plot!(x, same_e2, linewidth=2, label=L"$\textrm{Same case } E_2$", color =:blue)
plot!(x, same_f2 .- .5, linewidth=2, label=L"$\textrm{Same case } F_2$", color =:blue)
plot!(x, same_e3, linewidth=2, label=L"$\textrm{Same case } E_3$", color =:green)
plot!(x, same_f3 .- .5, linewidth=2, label=L"$\textrm{Same case } F_3$", color =:green)
plot!(x, same_e4, linewidth=2, label=L"$\textrm{Same case } E_4$", color =:purple)
plot!(x, same_f4 .- .5, linewidth=2, label=L"$\textrm{Same case } F_4$", color =:purple)

xaxis!(L"$\textrm{Grid Size}$")
savefig("plots/testcost_maxpaths.png")
savefig("plots/testcost_maxpaths.pdf")

# Plot path data:
P0 = [path_data4[i][1] for i in 1:length(path_data4)]
P1 = [path_data4[i][2] for i in 1:length(path_data4)]
P2 = [path_data4[i][3] for i in 1:length(path_data4)]
P3 = [path_data4[i][4] for i in 1:length(path_data4)]
P4 = [path_data4[i][5] for i in 1:length(path_data4)]

P0e = [path_data4_err[i][1] for i in 1:length(path_data4_err)]
P1e = [path_data4_err[i][2] for i in 1:length(path_data4_err)]
P2e = [path_data4_err[i][3] for i in 1:length(path_data4_err)]
P3e = [path_data4_err[i][4] for i in 1:length(path_data4_err)]
P4e = [path_data4_err[i][5] for i in 1:length(path_data4_err)]

# Plot analytical path probabilities
probs = []
probs_corner = []
# probs_line = []
for n in x
    N_m = n^2 - 4*(n-1)
    prob = N_m / n^2 * (N_m-1) / (n^2 - 1)
    push!(probs, prob)

    # Prob at least one node is a corner node
    prob_corner = 4/n^2 + (n^2 - 4)/(n^2) * (4/(n^2-1))
    push!(probs_corner, prob_corner)

    # Prob that users end on same row or column
    # prob_line = 2*(n-1)/(n^2 - 1)
    # push!(probs_line, prob_line)
end

##### Plot paths_maxpaths #####
plot4 = plot(x, zeros(length(x)), linewidth=2, xlims=(0, max_size+10),
color_palette = palette(:plasma, 4), legend=:topright,
guidefontsize=14, tickfontsize=12, legendfontsize=8, fontfamily="computer modern",
markersize=5, label=false)
plot!(x, probs_corner, linewidth=2, label=false)
plot!(x, 1 .- probs .- probs_corner, linewidth=2, label=false)
plot!(x, probs, linewidth=2, label=false )

# colorarray = (collect(0:length(P1)-1)./(length(P1)-1))
# palette([:purple, :green], 7)
# color_palette = cgrad(:rainbow, 4, categorical = true)
# cgrad(cs, [z], alpha = nothing, rev = false, scale = nothing, categorical = nothing)
# cgrad(:plasma, [0.1, 0.9])
plot!(x, P1, yerr = P1e, seriestype = :scatter, label=L"P_1", xlims=(0, max_size+10),
color_palette = palette(:plasma, 4), legend=:topright,
guidefontsize=14, tickfontsize=12, legendfontsize=8, fontfamily="computer modern",
markersize=5)
# cgrad(:matter, 4, categorical = true)
#:jet
# Specral_10
# :seaborn_bright
# plot!(x, P1, yerr = P1e, seriestype = :scatter, label=L"P_1", markersize=5)
plot!(x, P2, yerr = P2e, seriestype = :scatter, label=L"P_2", markersize=5)
plot!(x, P3, yerr = P3e, seriestype = :scatter, label=L"P_3", markersize=5)
plot!(x, P4, yerr = P4e, seriestype = :scatter, label=L"P_4", markersize=5)

xlabel!("Gridsize")
ylabel!("Path probability")
savefig("plots/path_maxpaths.png")
savefig("plots/path_maxpaths.pdf")

# # Plot efficiencies
# plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), yaxis=:log, color = :red, legend=:left)
# plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", yaxis=:log, color =:blue)
# plot!(x[3:end], loss_arr3[3:end], linewidth=2, label=L"$\eta_3$", yaxis=:log, color =:green)
# plot!(x[3:end], loss_arr4[3:end], linewidth=2, label=L"$\eta_4$", yaxis=:log, color =:purple)
#
# # # Plot fidelities
# z_arr1 = z_arr1 .- 0.5
# z_arr2 = z_arr2 .- 0.5
# z_arr3 = z_arr3 .- 0.5
# z_arr4 = z_arr4 .- 0.5
#
# plot!(x[1:end], z_arr1[1:end], linewidth=2, label=L"$F_1 - 0.5$", linestyle=:dash, color =:red, yaxis=:log)
# plot!(x[1:end], z_arr2[1:end], linewidth=2, label=L"$F_2 - 0.5$", linestyle=:dash, color =:blue, yaxis=:log)
# plot!(x[3:end], z_arr3[3:end], linewidth=2, label=L"$F_3 - 0.5$", linestyle=:dash, color =:green, yaxis=:log)
# plot!(x[3:end], z_arr4[3:end], linewidth=2, label=L"$F_4 - 0.5$", linestyle=:dash, color =:purple, yaxis=:log)
# xaxis!(L"$\textrm{Grid Size}$")
#
# savefig("plots/cost_maxpathslog.png")
# savefig("plots/cost_maxpathslog.pdf")

##### Plot bigfig_maxpaths #####
scale = 1.75
plot(plot1, plot2, plot4, plot3, layout = 4, size = (scale * 600, scale * 400))

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# Test the example scripts

module TestSIMPExample
    println("PointLoadCantilever Example")
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../../docs/src/literate/simp.jl"))
        end
    end
end

module TestBESOExample
    println("BESO Example")
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../../docs/src/literate/beso.jl"))
        end
    end
end

module TestGESOExample
    println("GESO Example")
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../../docs/src/literate/geso.jl"))
        end
    end
end


text:
module term

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Tesseract"
version = v"5.1.0"

# Collection of sources required to build Tesseract
sources = [
    ArchiveSource("https://github.com/tesseract-ocr/tesseract/archive/$(version).tar.gz",
                  "fdec8528d5a0ecc28ab5fff985e0b8ced60726f6ef33f54126f2868e323d4bd2"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/tesseract-*/
atomic_patch -p1 "$WORKSPACE/srcdir/patches/disable_fast_math.patch"
./autogen.sh
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
install_license ./LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libtesseract", :libtesseract),
    ExecutableProduct("tesseract", :tesseract),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Giflib_jll"),
    Dependency("JpegTurbo_jll"),
    Dependency("libpng_jll"),
    Dependency("Libtiff_jll"; compat="4.3.0"),
    Dependency("Zlib_jll"),
    Dependency("Leptonica_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    # Optional dependencies
    # Dependency("ICU_jll"),
    Dependency("Cairo_jll"; compat="1.16.1"),
    Dependency("Pango_jll"; compat="1.47.0"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"7",
    julia_compat="1.6"
)

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
import Base: sum, max
export argmax

function sum{T,N}(x::CuArray{T,N}, dim::Int)
    f = reducedim_kernel("v0 + v1", "0")
    f(x, dim)
end

function argmax{T,N}(x::CuArray{T,N}, dim::Int)
    f = reducedim_kernel("v0 < v1 ? v0 : v1", "0")
    f(x, dim)
end

function reshape3d(x::CuArray, dim::Int)
    dim1, dim2, dim3 = 1, size(x,dim), 1
    for i = 1:dim-1
        dim1 *= size(x,i)
    end
    for i = dim+1:ndims(x)
        dim3 *= size(x,i)
    end
    reshape(x, dim1, dim2, dim3)
end

function reducedim_kernel(op::String, v0::String)
    function kernel{T,N}(x::CuArray{T,N}, dim::Int)
        t = ctype(T)
        f = @nvrtc """
        $array_h
        template<typename T>
        __inline__ __device__ T warpReduce(T v0) {
            for (int offset = warpSize/2; offset > 0; offset /= 2) {
                T v1 = __shfl_down(v0, offset);
                v0 = $op;
            }
            return v0;
        }

        template<typename T>
        __inline__ __device__ T blockReduce(T val) {
            static __shared__ T shared[32];
            int lane = threadIdx.y % warpSize;
            int wid = threadIdx.y / warpSize;

            val = warpReduce<T>(val);
            if (lane == 0) shared[wid] = val;
            __syncthreads();

            val = (lane < 4) ? shared[lane] : 0;
            //if (wid == 0) val = warpReduce<T>(val);
            val = warpReduce<T>(val);
            return val;
        }

        __global__ void f(Array<$t,3> x, Array<$t,3> y) {
            int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
            int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
            int idx_z = blockIdx.z * blockDim.z + threadIdx.z;
            if (idx_y >= x.dims[1]) return;

            $t val = x(idx_x, idx_y, idx_z);
            val = blockReduce<$t>(val);
            if (threadIdx.y == 0) y(blockIdx.x, 0, blockIdx.z) = val;
        } """

        y = similar(x, ntuple(i -> i==dim ? 1 : size(x,i), N))
        x = reshape3d(x, dim)
        f(x, reshape3d(y,dim), dx=size(x,1), dy=size(x,2), dz=size(x,3), bx=1, by=128, sharedmem=sizeof(T)*32)
        y
    end
    kernel
end


text:
@testset "Peptide" begin
    n_steps = 100
    temp = 298.0u"K"
    s = System(
        joinpath(data_dir, "5XER", "gmx_coords.gro"),
        joinpath(data_dir, "5XER", "gmx_top_ff.top");
        loggers=Dict(
            "temp"   => TemperatureLogger(10),
            "coords" => CoordinateLogger(10),
            "energy" => TotalEnergyLogger(10),
            "writer" => StructureWriter(10, temp_fp_pdb),
        ),
    )
    simulator = VelocityVerlet(dt=0.0002u"ps", coupling=AndersenThermostat(temp, 10.0u"ps"))

    true_n_atoms = 5191
    @test length(s.atoms) == true_n_atoms
    @test length(s.coords) == true_n_atoms
    @test size(s.neighbor_finder.nb_matrix) == (true_n_atoms, true_n_atoms)
    @test size(s.neighbor_finder.matrix_14) == (true_n_atoms, true_n_atoms)
    @test length(s.pairwise_inters) == 2
    @test length(s.specific_inter_lists) == 3
    @test s.box_size == SVector(3.7146, 3.7146, 3.7146)u"nm"
    show(devnull, first(s.atoms))

    s.velocities = [velocity(a.mass, temp) .* 0.01 for a in s.atoms]
    @time simulate!(s, simulator, n_steps; parallel=false)

    traj = read(temp_fp_pdb, BioStructures.PDB)
    rm(temp_fp_pdb)
    @test BioStructures.countmodels(traj) == 10
    @test BioStructures.countatoms(first(traj)) == 5191
end

@testset "Peptide Float32" begin
    n_steps = 100
    temp = 298.0f0u"K"
    s = System(
        Float32,
        joinpath(data_dir, "5XER", "gmx_coords.gro"),
        joinpath(data_dir, "5XER", "gmx_top_ff.top");
        loggers=Dict(
            "temp"   => TemperatureLogger(typeof(1.0f0u"K"), 10),
            "coords" => CoordinateLogger(typeof(1.0f0u"nm"), 10),
            "energy" => TotalEnergyLogger(typeof(1.0f0u"kJ * mol^-1"), 10),
        ),
    )
    simulator = VelocityVerlet(dt=0.0002f0u"ps", coupling=AndersenThermostat(temp, 10.0f0u"ps"))

    s.velocities = [velocity(a.mass, Float32(temp)) .* 0.01f0 for a in s.atoms]
    @time simulate!(s, simulator, n_steps; parallel=false)
end

@testset "OpenMM protein comparison" begin
    openmm_dir = joinpath(data_dir, "openmm_6mrr")
    ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...)
    sys = System(joinpath(data_dir, "6mrr_equil.pdb"), ff; centre_coords=false)
    neighbors = find_neighbors(sys, sys.neighbor_finder)

    for inter in ("bond", "angle", "proptor", "improptor", "lj", "coul", "all")
        if inter == "all"
            pin = sys.pairwise_inters
        elseif inter == "lj"
            pin = sys.pairwise_inters[1:1]
        elseif inter == "coul"
            pin = sys.pairwise_inters[2:2]
        else
            pin = ()
        end

        if inter == "all"
            sils = sys.specific_inter_lists
        elseif inter == "bond"
            sils = sys.specific_inter_lists[1:1]
        elseif inter == "angle"
            sils = sys.specific_inter_lists[2:2]
        elseif inter == "proptor"
            sils = sys.specific_inter_lists[3:3]
        elseif inter == "improptor"
            sils = sys.specific_inter_lists[4:4]
        else
            sils = ()
        end

        sys_part = System(
            atoms=sys.atoms,
            pairwise_inters=pin,
            specific_inter_lists=sils,
            coords=sys.coords,
            box_size=sys.box_size,
            neighbor_finder=sys.neighbor_finder,
        )

        forces_molly = forces(sys_part, neighbors; parallel=false)
        forces_openmm = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "forces_$(inter)_only.txt"))))u"kJ * mol^-1 * nm^-1"
        # All force terms on all atoms must match at some threshold
        @test !any(d -> any(abs.(d) .> 1e-6u"kJ * mol^-1 * nm^-1"), forces_molly .- forces_openmm)

        E_molly = potential_energy(sys_part, neighbors)
        E_openmm = readdlm(joinpath(openmm_dir, "energy_$(inter)_only.txt"))[1] * u"kJ * mol^-1"
        # Energy must match at some threshold
        @test E_molly - E_openmm < 1e-5u"kJ * mol^-1"
    end

    # Run a short simulation with all interactions
    n_steps = 100
    simulator = VelocityVerlet(dt=0.0005u"ps")
    velocities_start = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_300K.txt"))))u"nm * ps^-1"
    sys.velocities = deepcopy(velocities_start)
    @test kinetic_energy(sys) ≈ 65521.87288132431u"kJ * mol^-1"
    @test temperature(sys) ≈ 329.3202932884933u"K"

    simulate!(sys, simulator, n_steps; parallel=true)

    coords_openmm = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "coordinates_$(n_steps)steps.txt"))))u"nm"
    vels_openmm   = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_$(n_steps)steps.txt" ))))u"nm * ps^-1"

    coords_diff = sys.coords .- wrap_coords_vec.(coords_openmm, (sys.box_size,))
    vels_diff = sys.velocities .- vels_openmm
    # Coordinates and velocities at end must match at some threshold
    @test maximum(maximum(abs.(v)) for v in coords_diff) < 1e-9u"nm"
    @test maximum(maximum(abs.(v)) for v in vels_diff  ) < 1e-6u"nm * ps^-1"

    # Test with no units
    ff_nounits = OpenMMForceField(
        joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...;
        units=false,
    )
    sys_nounits = System(
        joinpath(data_dir, "6mrr_equil.pdb"),
        ff_nounits;
        velocities=deepcopy(ustrip_vec.(velocities_start)),
        units=false,
        centre_coords=false,
    )
    simulator_nounits = VelocityVerlet(dt=0.0005)
    @test kinetic_energy(sys_nounits)u"kJ * mol^-1" ≈ 65521.87288132431u"kJ * mol^-1"
    @test temperature(sys_nounits)u"K" ≈ 329.3202932884933u"K"

    simulate!(sys_nounits, simulator_nounits, n_steps; parallel=true)

    coords_diff = sys_nounits.coords * u"nm" .- wrap_coords_vec.(coords_openmm, (sys.box_size,))
    vels_diff = sys_nounits.velocities * u"nm * ps^-1" .- vels_openmm
    @test maximum(maximum(abs.(v)) for v in coords_diff) < 1e-9u"nm"
    @test maximum(maximum(abs.(v)) for v in vels_diff  ) < 1e-6u"nm * ps^-1"

    params_dic = extract_parameters(sys_nounits, ff_nounits)
    @test length(params_dic) == 639
    atoms_grad, pis_grad, sis_grad = inject_gradients(sys_nounits, params_dic)
    @test atoms_grad == sys_nounits.atoms
    @test pis_grad == sys_nounits.pairwise_inters

    # Test the same simulation on the GPU
    if run_gpu_tests
        sys = System(
            joinpath(data_dir, "6mrr_equil.pdb"),
            ff;
            velocities=cu(deepcopy(velocities_start)),
            gpu=true,
            centre_coords=false,
        )
        simulate!(sys, simulator, n_steps)

        coords_diff = Array(sys.coords) .- wrap_coords_vec.(coords_openmm, (sys.box_size,))
        vels_diff = Array(sys.velocities) .- vels_openmm
        @test maximum(maximum(abs.(v)) for v in coords_diff) < 1e-9u"nm"
        @test maximum(maximum(abs.(v)) for v in vels_diff  ) < 1e-6u"nm * ps^-1"

        sys_nounits = System(
            joinpath(data_dir, "6mrr_equil.pdb"),
            ff_nounits;
            velocities=cu(deepcopy(ustrip_vec.(velocities_start))),
            units=false,
            gpu=true,
            centre_coords=false,
        )
        simulate!(sys_nounits, simulator_nounits, n_steps)

        coords_diff = Array(sys_nounits.coords * u"nm") .- wrap_coords_vec.(coords_openmm, (sys.box_size,))
        vels_diff = Array(sys_nounits.velocities * u"nm * ps^-1") .- vels_openmm
        @test maximum(maximum(abs.(v)) for v in coords_diff) < 1e-9u"nm"
        @test maximum(maximum(abs.(v)) for v in vels_diff  ) < 1e-6u"nm * ps^-1"

        params_dic_gpu = extract_parameters(sys_nounits, ff_nounits)
        @test params_dic == params_dic_gpu
        atoms_grad, pis_grad, sis_grad = inject_gradients(sys_nounits, params_dic_gpu)
        @test atoms_grad == sys_nounits.atoms
        @test pis_grad == sys_nounits.pairwise_inters
    end
end

@testset "Implicit solvent" begin
    openmm_dir = joinpath(data_dir, "openmm_6mrr")
    ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "his.xml"])...)
    sys = System(
        joinpath(data_dir, "6mrr_nowater.pdb"),
        ff;
        box_size=SVector(100.0, 100.0, 100.0)u"nm",
        implicit_solvent="obc2",
        dist_cutoff=5.0u"nm",
        nl_dist=5.0u"nm",
    )
    neighbors = find_neighbors(sys)

    forces_molly = forces(sys, neighbors)
    forces_openmm = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "forces_obc2.txt"))))u"kJ * mol^-1 * nm^-1"
    @test !any(d -> any(abs.(d) .> 1e-3u"kJ * mol^-1 * nm^-1"), forces_molly .- forces_openmm)

    E_molly = potential_energy(sys, neighbors)
    E_openmm = readdlm(joinpath(openmm_dir, "energy_obc2.txt"))[1] * u"kJ * mol^-1"
    @test E_molly - E_openmm < 1e-2u"kJ * mol^-1"
end


text:
using JLD, DataFrames

star_cat = load("q1q17_dr25_gaia_fgk.jld", "stellar_catalog")

star_catF = star_cat[(6000.<=star_cat[:teff].<7300),:]
star_catG = star_cat[(5300.<=star_cat[:teff].<6000),:]
star_catK = star_cat[(3900.<=star_cat[:teff].<5300),:]

save("q1q17_dr25_gaia_f.jld", "stellar_catalog", star_catF)
save("q1q17_dr25_gaia_g.jld", "stellar_catalog", star_catG)
save("q1q17_dr25_gaia_k.jld", "stellar_catalog", star_catK)


text:
module SpikingNN

using Distributions
using RecipesBase
using StructArrays
using CUDA
using Adapt
using TensorCast
using LoopVectorization
using Random

using Reexport

export  AbstractSynapse,
        AbstractThreshold,
        AbstractCell,
        AbstractInput,
        AbstractLearner

export  excite!, simulate!, evaluate!, reset!,
        getvoltage,
        Soma, Neuron,
        LIF, SRM0,
        George, STDP,
        Population,
        neurons, synapses,
        prespike!, postspike!, record!, update!,
        ConstantRate, StepCurrent, PoissonInput,
        InputPopulation,
        Network, connect!,
        cpu, gpu

include("utils.jl")

# prototypes
function excite! end
function spike! end
function evaluate! end
function reset! end
function isactive end

include("synapse.jl")
include("threshold.jl")
include("neuron.jl")
include("models/lif.jl")
include("models/srm0.jl")
include("inputs.jl")
include("learning.jl")
include("population.jl")
include("network.jl")
include("gpu.jl")

end


text:
using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Molenaar", "Figure 3", "Plots")
resultspath = joinpath("Molenaar", "Results")

x = [res["mu"] for res in results]
y = [res["used_ribo"] for res in results]

a = [] # transporter rate
b = [] # matabolic rate
c = [] # ribosomale rate
d = [] # lipid biosynthesis rate

for i = 1:31
    push!(a, y[i][1])
    push!(b, y[i][2])
    push!(c, y[i][3])
    push!(d, y[i][4])
end

a = Float64.(a)
b = Float64.(b)
c = Float64.(c)
d = Float64.(d)

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a) # blue
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b) # yellow
line2 = lines!(ax, x, b)
scatter3 = scatter!(ax, x, c) # green
line3 = lines!(ax, x, c)
scatter4 = scatter!(ax, x, d) # pink
line4 = lines!(ax, x, d)

ax.xlabel = "Growth rate (1/h)"
ax.ylabel = "Protein fraction"
fig[1, 1] = Legend(
    fig,
    [[scatter1, line1], [scatter2, line2], [scatter3, line3], [scatter4, line4]],
    ["Transporter", "Metabolism", "Ribosomes", "Lipids"],
    tellheight = false,
    tellwidth = false,
    halign = :right,
    valign = :top,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Figure_3C.pdf"), fig)


text:
@defcomp emissions begin
    regions = Index()

    mitigationcost = Variable(index=[time,regions])
    ch4cost = Variable(index=[time,regions])
    ch4costindollar = Variable(index=[time,regions])
    n2ocost = Variable(index=[time,regions])
    ryg = Variable(index=[time,regions])
    energint = Variable(index=[time,regions])
    emissint = Variable(index=[time,regions])
    emission = Variable(index=[time,regions])
    emissionwithforestry = Variable(index=[time,regions])
    sf6 = Variable(index=[time,regions])
    reei = Variable(index=[time,regions])
    rcei = Variable(index=[time,regions])
    energuse = Variable(index=[time,regions])
    seei = Variable(index=[time,regions])
    scei = Variable(index=[time,regions])
    ch4 = Variable(index=[time,regions])
    n2o = Variable(index=[time,regions])
    n2ored = Variable(index=[time,regions])
    taxpar = Variable(index=[time,regions])
    co2red = Variable(index=[time,regions])
    know = Variable(index=[time,regions])
    perm = Variable(index=[time,regions])
    cumaeei = Variable(index=[time,regions])
    ch4red = Variable(index=[time,regions])

    minint = Variable(index=[time])
    globknow = Variable(index=[time])

    mco2 = Variable(index=[time])
    globch4 = Variable(index=[time])
    globn2o = Variable(index=[time])
    globsf6 = Variable(index=[time])

    cumglobco2 = Variable(index=[time])
    cumglobch4 = Variable(index=[time])
    cumglobn2o = Variable(index=[time])
    cumglobsf6 = Variable(index=[time])

    taxmp = Parameter(index=[regions])
    sf60 = Parameter(index=[regions])
    gdp90 = Parameter(index=[regions])
    pop90 = Parameter(index=[regions])
    ch4par1 = Parameter(index=[regions])
    ch4par2 = Parameter(index=[regions])
    n2opar1 = Parameter(index=[regions])
    n2opar2 = Parameter(index=[regions])
    emissint0 = Parameter(index=[regions])

    forestemm = Parameter(index=[time,regions])
    aeei = Parameter(index=[time,regions])
    acei = Parameter(index=[time,regions])
    ch4em = Parameter(index=[time,regions])
    n2oem = Parameter(index=[time,regions])
    currtax = Parameter(index=[time,regions])
    currtaxch4 = Parameter(index=[time,regions])
    currtaxn2o = Parameter(index=[time,regions])
#    pgrowth = Parameter(index=[time,regions])
    ypcgrowth = Parameter(index=[time,regions])
    income = Parameter(index=[time,regions])
    population = Parameter(index=[time,regions])

    sf6gdp = Parameter(default = 0.00022628870292887)
    sf6ypc = Parameter(default = -2.46698744769874e-6)
    knowpar = Parameter(default = 0.9)
    knowgpar = Parameter(default = 0.1)
    gwpch4 = Parameter(default = 0.077)
    gwpn2o = Parameter(default = 0.361666666666667)

    taxconstant = Parameter(default = 0.784)
    taxemint = Parameter(default = 0.084)
    taxthreshold = Parameter(default = 100)
    taxdepreciation = Parameter(default = 0.1)
    maxcostfall = Parameter(default = 10)

    ch4add = Parameter(default = 0)
    n2oadd = Parameter(default = 0)
    sf6add = Parameter(default = 0)

    function run_timestep(p, v, d, t)

        if is_first(t)
            for r in d.regions
                v.energint[t, r] = 1
                v.energuse[t, r] = p.income[t,r]
                v.emissint[t, r] = p.emissint0[r]
                v.emission[t, r] = v.emissint[t, r] / v.energuse[t, r]
                v.ch4cost[t, r] = 0
                v.n2ocost[t, r] = 0
                v.ryg[t, r] = 0
                v.reei[t, r] = 0
                v.rcei[t, r] = 0
                v.seei[t, r] = 0
                v.scei[t, r] = 0
                v.co2red[t, r] = 0
                v.know[t, r] = 1
                v.ch4red[t, r] = 0
                v.n2ored[t, r] = 0
                v.mitigationcost[t, r] = 0
            end

            v.globknow[t] = 1
            v.cumglobco2[t] = 0.0
            v.cumglobch4[t] = 0.0
            v.cumglobn2o[t] = 0.0
            v.cumglobsf6[t] = 0.0

            # SocioEconomicState.minint[t]=Inf
            minint = Inf
            for r in d.regions
                if v.emission[t, r] / p.income[t, r] < minint
                    minint = v.emission[t, r] / p.income[t, r]
                end
            end
            # v.minint[t] = minint
            v.minint[t] = 0
        else
            # Calculate emission and carbon intensity
            for r in d.regions
                v.energint[t, r] = (1.0 - 0.01 * p.aeei[t, r] - v.reei[t - 1, r]) * v.energint[t - 1, r]
                v.emissint[t, r] = (1.0 - 0.01 * p.acei[t, r] - v.rcei[t - 1, r]) * v.emissint[t - 1, r]
            end

            # Calculate sf6 emissions
            for r in d.regions
                v.sf6[t, r] = (p.sf60[r] + p.sf6gdp * (p.income[t, r] - p.gdp90[r]) + p.sf6ypc * (p.income[t - 1, r] / p.population[t - 1, r] - p.gdp90[r] / p.pop90[r])) * (gettime(t) <= 2010 ? 1 + (gettime(t) - 1990) / 40.0 : 1.0 + (60.0 - 40.0) / 40.0) * (gettime(t) > 2010 ? 0.99^(gettime(t) - 2010) : 1.0)
            end

            # Check for unrealistic values
            for r in d.regions
                if v.sf6[t, r] < 0.0
                    v.sf6[t, r] = 0
                end
            end

            # Calculate energy use
            for r in d.regions
                v.energuse[t, r] = (1 - v.seei[t - 1, r]) * v.energint[t, r] * p.income[t, r]
            end

            # Calculate co2 emissions
            for r in d.regions
                v.emission[t, r] = (1 - v.scei[t - 1, r]) * v.emissint[t, r] * v.energuse[t, r]
                v.emissionwithforestry[t, r] = v.emission[t, r] + p.forestemm[t, r]
            end

            # Calculate ch4 emissions
            for r in d.regions
                v.ch4[t, r] = p.ch4em[t, r] * (1 - v.ch4red[t - 1, r])
            end

            # Calculate n2o emissions
            for r in d.regions
                v.n2o[t, r] = p.n2oem[t, r] * (1 - v.n2ored[t - 1, r])
            end

            # TODO RT check
            for r in d.regions
                if (v.emission[t, r] / p.income[t, r] - v.minint[t - 1] <= 0)
                    v.taxpar[t, r] = p.taxconstant
                else
                    v.taxpar[t, r] = p.taxconstant - p.taxemint * sqrt(v.emission[t, r] / p.income[t, r] - v.minint[t - 1])
                end
            end

            for r in d.regions
                v.co2red[t, r] = p.currtax[t, r] * v.emission[t, r] * v.know[t - 1, r] * v.globknow[t - 1] / 2 / v.taxpar[t, r] / p.income[t, r] / 1000

                if (v.co2red[t, r] < 0)
                    v.co2red[t, r] = 0
                elseif (v.co2red[t, r] > 0.99)
                    v.co2red[t, r] = 0.99
                end
            end

            for r in d.regions
                v.ryg[t, r] = v.taxpar[t, r] * v.co2red[t, r]^2 / v.know[t - 1, r] / v.globknow[t - 1]
            end

            # TODO RT check
            for r in d.regions
                v.perm[t, r] = 1.0 - 1.0 / p.taxthreshold * p.currtax[t, r] / (1 + 1.0 / p.taxthreshold * p.currtax[t, r])
            end

            for r in d.regions
                v.reei[t, r] = v.perm[t, r] * 0.5 * v.co2red[t, r]
            end

            # TODO RT check
            for r in d.regions
                if (p.currtax[t, r] < p.taxthreshold)
                    v.rcei[t, r] = v.perm[t, r] * 0.5 * v.co2red[t, r]^2
                else
                    v.rcei[t, r] = v.perm[t, r] * 0.5 * v.co2red[t, r]
                end
            end

            # TODO RT check
            # TODO RT what is the 1.7?
            for r in d.regions
                v.seei[t, r] = (1.0 - p.taxdepreciation) * v.seei[t - 1, r] + (1.0 - v.perm[t, r]) * 0.5 * v.co2red[t, r] * 1.7
            end

            for r in d.regions
                if (p.currtax[t, r] < 100)
                    v.scei[t, r] = 0.9 * v.scei[t - 1, r] + (1 - v.perm[t, r]) * 0.5 * v.co2red[t, r]^2
                else
                    v.scei[t, r] = 0.9 * v.scei[t - 1, r] + (1 - v.perm[t, r]) * 0.5 * v.co2red[t, r] * 1.7
                end
            end

            # TODO RT check
            for r in d.regions
                v.know[t, r] = v.know[t - 1, r] * sqrt(1 + p.knowpar * v.co2red[t, r])

                if (v.know[t, r] > sqrt(p.maxcostfall))
                    v.know[t, r] = sqrt(p.maxcostfall)
                end
            end

            v.globknow[t] = v.globknow[t - 1]
            for r in d.regions
                v.globknow[t] = v.globknow[t] * sqrt(1 + p.knowgpar * v.co2red[t, r])
            end

            if (v.globknow[t] > 3.16)
                v.globknow[t] = 3.16
            end

            for r in d.regions
                v.ch4red[t, r] = p.currtaxch4[t, r] * p.ch4em[t, r] / 2 / p.ch4par1[r] / p.ch4par2[r] / p.ch4par2[r] / p.income[t, r] / 1000

                if (v.ch4red[t, r] < 0)
                    v.ch4red[t, r] = 0
                elseif (v.ch4red[t, r] > 0.99)
                    v.ch4red[t, r] = 0.99
                end
            end

            for r in d.regions
                v.n2ored[t, r] = p.gwpn2o * p.currtaxn2o[t, r] * p.n2oem[t, r] / 2 / p.n2opar1[r] / p.n2opar2[r] / p.n2opar2[r] / p.income[t, r] / 1000

                if (v.n2ored[t, r] < 0)
                    v.n2ored[t, r] = 0
                elseif (v.n2ored[t, r] > 0.99)
                    v.n2ored[t, r] = 0.99
                end
            end

            for r in d.regions
                v.ch4cost[t, r] = p.ch4par1[r] * p.ch4par2[r]^2 * v.ch4red[t, r]^2
                v.ch4costindollar[t, r] = v.ch4cost[t, r] * p.income[t, r]
            end

            for r in d.regions
                v.n2ocost[t, r] = p.n2opar1[r] * p.n2opar2[r]^2 * v.n2ored[t, r]^2
            end

            minint = Inf
            for r in d.regions
                if (v.emission[t, r] / p.income[t, r] < minint)
                    minint = v.emission[t, r] / p.income[t, r]
                end
            end
            v.minint[t] = minint

            for r in d.regions
                if gettime(t) > 2000
                    v.cumaeei[t, r] = v.cumaeei[t - 1, r] * (1.0 - 0.01 * p.aeei[t, r] - v.reei[t, r] + v.seei[t - 1, r] - v.seei[t, r])
                else
                    v.cumaeei[t, r] = 1.0
                end
            end

            for r in d.regions
                #v.mitigationcost[t, r] = (p.taxmp[r] * v.ryg[t, r] /*+ v.ch4cost[t, r]*/ + v.n2ocost[t, r]) * p.income[t, r]
                v.mitigationcost[t, r] = (p.taxmp[r] * v.ryg[t, r] + v.n2ocost[t, r]) * p.income[t, r]
            end

            globco2 = 0.
            globch4 = 0.
            globn2o = 0.
            globsf6 = 0.

            for r in d.regions
                globco2 = globco2 + v.emissionwithforestry[t, r]
                globch4 = globch4 + v.ch4[t, r]
                globn2o = globn2o + v.n2o[t, r]
                globsf6 = globsf6 + v.sf6[t, r]
            end

            v.mco2[t] = globco2
            v.globch4[t] = max(0.0, globch4 + (gettime(t) > 2000 ? p.ch4add * (gettime(t) - 2000) : 0.0))
            v.globn2o[t] = max(0.0, globn2o + (gettime(t) > 2000 ? p.n2oadd * (gettime(t) - 2000) : 0.0))
            v.globsf6[t] = max(0.0, globsf6 + (gettime(t) > 2000 ? p.sf6add * (gettime(t) - 2000) : 0.0))

            v.cumglobco2[t] = v.cumglobco2[t - 1] + v.mco2[t]
            v.cumglobch4[t] = v.cumglobch4[t - 1] + v.globch4[t]
            v.cumglobn2o[t] = v.cumglobn2o[t - 1] + v.globn2o[t]
            v.cumglobsf6[t] = v.cumglobsf6[t - 1] + v.globsf6[t]
        end
    end
end


text:
"""
Modia module with rotational component models (inspired from Modelica Standard Library.

* Developer: Hilding Elmqvist, Mogram AB
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)

"""
module Rotational

export Flange, Inertia, Spring, SpringDamper, EMF, IdealGear, Torque, CurrentSensor
using ..Instantiation
using ..Electric
using ..Blocks

using Modia

@model Flange begin
  phi=Var(T=Float64, size=())
  tau=Float(flow=true, size=())
end

@model Inertia begin
    # 1D-rotational component with inertia"
  J=Parameter(0, min=0)   # Moment of inertia

  flange_a=Flange()  # Left flange of shaft
  flange_b=Flange()  # Right flange of shaft

  phi=Float(start=0.0)
  w=Float(start=0.0)
  a=Float()
@equations begin
  phi = flange_a.phi
  phi = flange_b.phi
  w = der(phi)
  a = der(w)
  J*a = flange_a.tau + flange_b.tau
  end
end

@model PartialCompliant begin
    # Partial model for the compliant connection of two rotational 1-dim. shaft flanges
  phi_rel=Float()
  tau=Float(size=())
  flange_a=Flange()
  flange_b=Flange()
@equations begin
  phi_rel = flange_b.phi - flange_a.phi
  flange_b.tau = tau
  flange_a.tau = -tau
  end
end

@model Spring begin
    # Linear 1D rotational spring
  @extends PartialCompliant()
  @inherits tau, phi_rel
  c=Parameter(min=0, start=1.0e5)       # Spring constant
  phi_rel0=0   # Unstretched spring angle
@equations begin
  tau = c*(phi_rel - phi_rel0)
  end
end

@model SpringDamper begin
    # Linear 1D rotational spring
  @extends PartialCompliant(phi_rel=Float(size=(), state=false))
  @inherits tau, phi_rel
  c=Parameter(min=0, start=1.0e5)       # Spring constant
  d=Parameter()
  phi_rel0=0   # Unstretched spring angle
@equations begin
  tau = c*(phi_rel - phi_rel0) + d*der(phi_rel)
  end
end

@model EMF begin
    # Electromotoric force (electric/mechanic transformer)
  k=1

  p=Pin()
  n=Pin()
  flange=Flange()

  v=Float(size=())
  i=Float(size=())
  phi=Float(state=false, size=())
  w=Float(size=())
@equations begin
  v = p.v - n.v
  0 = p.i + n.i
  i = p.i

  phi = flange.phi
  w = der(phi)
  k*w = v
  flange.tau = -k*i
  end
end

@model IdealGear begin
    # Ideal gear without inertia
#  @extends Gear()

#  @extends PartialElementaryTwoFlangesAndSupport2()
  flange_a=Flange()  # Left flange of shaft
  flange_b=Flange()  # Right flange of shaft
  ratio=1     # Transmission ratio (flange_a.phi/flange_b.phi)
  phi_a=Float(size=())
  phi_b=Float(size=())
@equations begin
  phi_a = flange_a.phi
  phi_b = flange_b.phi
  phi_a = ratio*phi_b
  0 = ratio*flange_a.tau + flange_b.tau
  end
end

@model Torque begin # Input signal acting as external torque on a flange
  flange=Flange()  # Flange of shaft
  tau=Float()
@equations begin
  flange.tau = -tau
  end
end

@model RotationalSensor begin

end

@model CurrentSensor begin
    # Sensor to measure the current in a branch
  @extends RotationalSensor()
  p=Pin()
  n=Pin()
  i=Float(size=())
@equations begin
  p.v = n.v
  p.i = i
  n.i = -i
  end
end

end