text:
# OpenCL.Kernel
type Kernel <: CLObject
    id :: CL_kernel

    function Kernel(k::CL_kernel, retain=false)
        if retain
            @check api.clRetainKernel(k)
        end
        kernel = new(k)
        finalizer(kernel, _finalize)
        return kernel
    end
end

function _finalize(k::Kernel)
    if k.id != C_NULL
        @check api.clReleaseKernel(k.id)
        k.id = C_NULL
    end
end

Base.pointer(k::Kernel) = k.id

Base.show(io::IO, k::Kernel) = begin
    print(io, "OpenCL.Kernel(\"$(k[:name])\" nargs=$(k[:num_args]))")
end

Base.getindex(k::Kernel, kinfo::Symbol) = info(k, kinfo)

function Kernel(p::Program, kernel_name::String)
    for (dev, status) in info(p, :build_status)
        if status != CL_BUILD_SUCCESS
            msg = "OpenCL.Program has to be built before Kernel constructor invoked"
            throw(ArgumentError(msg))
        end
    end
    err_code = Ref{CL_int}()
    kernel_id = api.clCreateKernel(p.id, kernel_name, err_code)
    if err_code[] != CL_SUCCESS
        throw(CLError(err_code[]))
    end
    return Kernel(kernel_id)
end

immutable LocalMem{T}
    nbytes::Csize_t
end

function LocalMem{T}(::Type{T}, len::Integer)
    @assert len > 0
    nbytes = sizeof(T) * len
    return LocalMem{T}(convert(Csize_t, nbytes))
end

Base.ndims(l::LocalMem) = 1
Base.eltype{T}(l::LocalMem{T}) = T
Base.sizeof{T}(l::LocalMem{T}) = l.nbytes
Base.length{T}(l::LocalMem{T}) = Int(l.nbytes ÷ sizeof(T))

function set_arg!(k::Kernel, idx::Integer, arg::Void)
    @assert idx > 0
    @check api.clSetKernelArg(k.id, cl_uint(idx-1), sizeof(CL_mem), C_NULL)
    return k
end

function set_arg!(k::Kernel, idx::Integer, arg::Ptr{Void})
    if arg != C_NULL
        throw(AttributeError("set_arg! for void pointer $arg is undefined"))
    end
    set_arg!(k, idx, nothing)
end

function set_arg!(k::Kernel, idx::Integer, arg::CLMemObject)
    @assert idx > 0
    arg_boxed = Ref{typeof(arg.id)}(arg.id)
    @check api.clSetKernelArg(k.id, cl_uint(idx-1), sizeof(CL_mem), arg_boxed)
    return k
end

function set_arg!(k::Kernel, idx::Integer, arg::LocalMem)
    @assert idx > 0 "Kernel idx must be bigger 0"
    @check api.clSetKernelArg(k.id, cl_uint(idx-1), arg.nbytes, C_NULL)
    return k
end

function _contains_different_layout{T}(::Type{T})
    sizeof(T) == 0 && return true
    nfields(T) == 0 && return false
    for fname in fieldnames(T)
        contains_different_layout(fieldtype(T, fname)) && return true
    end
    return false
end

function contains_different_layout{T <: Union{Float32, Float64, Int8, Int32, Int64, UInt8, UInt32, UInt64}}(::Type{NTuple{3, T}})
    true
end

"""
    contains_different_layout(T)

Empty types and NTuple{3, CLNumber} have different layouts and need to be replaced
(Where `CLNumber <: Union{Float32, Float64, Int8, Int32, Int64, UInt8, UInt32, UInt64}`)
TODO: Float16 + Int16 should also be in CLNumbers
"""
@generated function contains_different_layout{T}(::Type{T})
    :($(_contains_different_layout(T)))
end

function struct2tuple{T}(x::T)
    ntuple(Val{nfields(T)}) do i
        getfield(x, i)
    end
end

"""
    replace_different_layout(x::T) where T

Replaces types with a layout different from OpenCL.
See [contains_different_layout(T)](@ref) for information what types those are!
"""
function replace_different_layout{T}(x::T)
    !contains_different_layout(T) && return x
    if nfields(x) == 0
        replace_different_layout((), (x,))
    elseif T <: Tuple
        replace_different_layout((), x)
    else
        replace_different_layout((), struct2tuple(x))
    end
end

replace_different_layout{N}(red::NTuple{N, Any}, rest::Tuple{}) = red
function replace_different_layout{N}(red::NTuple{N, Any}, rest)
    elem1 = first(rest)
    T = typeof(elem1)
    repl = if sizeof(T) == 0 && nfields(T) == 0
        Int32(0)
    elseif contains_different_layout(T)
        replace_different_layout(elem1)
    else
        elem1
    end
    replace_different_layout((red..., repl), Base.tail(rest))
end

# TODO UInt16/Float16?
# Handle different sizes of OpenCL Vec3, which doesn't agree with julia
function replace_different_layout{T <: Union{Float32, Float64, Int8, Int32, Int64, UInt8, UInt32, UInt64}}(arg::NTuple{3, T})
    pad = T(0)
    (arg..., pad)
end

function to_cl_ref{T}(arg::T)
    if !Base.datatype_pointerfree(T)
        error("Types should not contain pointers: $T")
    end
    if contains_different_layout(T)
        x = replace_different_layout(arg)
        return Base.RefValue(x), sizeof(x)
    end
    Base.RefValue(arg), sizeof(arg)
end

Base.@pure datatype_align(x::T) where {T} = datatype_align(T)
Base.@pure function datatype_align(::Type{T}) where {T}
    # typedef struct {
    #     uint32_t nfields;
    #     uint32_t alignment : 9;
    #     uint32_t haspadding : 1;
    #     uint32_t npointers : 20;
    #     uint32_t fielddesc_type : 2;
    # } jl_datatype_layout_t;
    field = T.layout + sizeof(UInt32)
    unsafe_load(convert(Ptr{UInt16}, field)) & convert(Int16, 2^9-1)
end

function set_arg!{T}(k::Kernel, idx::Integer, arg::T)
    @assert idx > 0 "Kernel idx must be bigger 0"
    ref, tsize = to_cl_ref(arg)
    err = api.clSetKernelArg(k.id, cl_uint(idx - 1), tsize, ref)
    if err == CL_INVALID_ARG_SIZE
        error("""
            Julia and OpenCL type don't match at kernel argument $idx: Found $T.
            Please make sure to define OpenCL structs correctly!
            You should be generally fine by using `__attribute__((packed))`, but sometimes the alignment of fields is different from Julia.
            Consider the following example:
                ```
                //packed
                // Tuple{NTuple{3, Float32}, Void, Float32}
                struct __attribute__((packed)) Test{
                    float3 f1;
                    int f2; // empty type gets replaced with Int32 (no empty types allowed in OpenCL)
                    // you might need to define the alignement of fields to match julia's layout
                    float f3; // for the types used here the alignement matches though!
                };
                // this is a case where Julia and OpenCL packed alignment would differ, so we need to specify it explicitely
                // Tuple{Int64, Int32}
                struct __attribute__((packed)) Test2{
                    long f1;
                    int __attribute__((aligned (8))) f2; // opencl would align this to 4 in packed layout, while Julia uses 8!
                };
                ```
            You can use `c.datatype_align(T)` to figure out the alignment of a Julia type!
        """)
    end
    @check err
    return k
end

function set_args!(k::Kernel, args...)
    for (i, a) in enumerate(args)
        set_arg!(k, i, a)
    end
end

function work_group_info(k::Kernel, winfo::CL_kernel_work_group_info, d::Device)
    if (winfo == CL_KERNEL_LOCAL_MEM_SIZE ||
        winfo == CL_KERNEL_PRIVATE_MEM_SIZE)
        result1 = Ref{CL_ulong}(0)
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo,
                                            sizeof(CL_ulong), result1, C_NULL)
        return Int(result1[])
    elseif winfo == CL_KERNEL_COMPILE_WORK_GROUP_SIZE
        # Intel driver has a bug so we can't query the required size.
        # As specified by [1] the return value in this case is size_t[3].
        # [1] https://www.khronos.org/registry/OpenCL/sdk/1.2/docs/man/xhtml/clGetKernelWorkGroupInfo.html
        @assert sizeof(Csize_t) == sizeof(Int)
        result2 = Vector{Int}(3)
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo, 3*sizeof(Int), result2, C_NULL)
        return result2
    else
        result = Ref{Csize_t}(0)
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo,
                                            sizeof(CL_ulong), result, C_NULL)
        return Int(result[])
    end
end

function work_group_info(k::Kernel, winfo::Symbol, d::Device)
    if winfo == :size
        work_group_info(k, CL_KERNEL_WORK_GROUP_SIZE, d)
    elseif winfo == :compile_size
        work_group_info(k, CL_KERNEL_COMPILE_WORK_GROUP_SIZE, d)
    elseif winfo == :local_mem_size
        work_group_info(k, CL_KERNEL_LOCAL_MEM_SIZE, d)
    elseif winfo == :private_mem_size
        work_group_info(k, CL_KERNEL_PRIVATE_MEM_SIZE, d)
    elseif winfo == :prefered_size_multiple
        work_group_info(k, CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE, d)
    else
        throw(ArgumentError(("Unknown work_group_info flag: :$winfo")))
    end
end

# produce a cl.call thunk with kernel queue, global/local sizes
Base.getindex(k::Kernel, args...) = begin
    if length(args) < 2 || length(args) > 3
        throw(ArgumentError("kernel must be called with a queue & global size as arguments"))
    end
    if !(isa(args[1], CmdQueue))
        throw(ArgumentError("kernel first argument must a a CmdQueue"))
    end
    if !(isa(args[2], Dims)) || length(args[2]) > 3
        throw(ArgumentError("kernel global size must be of Dims type (dim <= 3)"))
    end
    if length(args) == 3 && (!(isa(args[3], Dims)) || length(args[3]) > 3)
        throw(ArgumentError("kernel local size must be of Dims type (dim <= 3)"))
    end
    queue = args[1]
    global_size = args[2]
    local_size  = length(args) == 3 ? args[3] : nothing
    # TODO: we cannot pass keywords in anon functions yet, return kernel call thunk
    return (args...) -> queue(k, global_size, local_size, args...)
end

# blocking kernel call that finishes queue
function (q::CmdQueue)(k::Kernel, global_work_size, local_work_size,
                      args...; global_work_offset=nothing,
                      wait_on::Union{Void,Vector{Event}}=nothing)
    set_args!(k, args...)
    evt = enqueue_kernel(q, k,
                         global_work_size,
                         local_work_size,
                         global_work_offset=global_work_offset,
                         wait_on=wait_on)
    finish(q)
    return evt
end

function enqueue_kernel(q::CmdQueue, k::Kernel, global_work_size)
    enqueue_kernel(q, k, global_work_size, nothing)
end

function enqueue_kernel(q::CmdQueue,
                                k::Kernel,
                                global_work_size,
                                local_work_size;
                                global_work_offset=nothing,
                                wait_on::Union{Void,Vector{Event}}=nothing)
    device = q[:device]
    max_work_dim = device[:max_work_item_dims]
    work_dim     = length(global_work_size)
    if work_dim > max_work_dim
        throw(ArgumentError("global_work_size has max dim of $max_work_dim"))
    end
    gsize = Array{Csize_t}(work_dim)
    for (i, s) in enumerate(global_work_size)
        gsize[i] = s
    end

    goffset = C_NULL
    if global_work_offset !== nothing
        if length(global_work_offset) > max_work_dim
            throw(ArgumentError("global_work_offset has max dim of $max_work_dim"))
        end
        if length(global_work_offset) != work_dim
            throw(ArgumentError("global_work_size and global_work_offset have differing dims"))
        end
        goffset = Array{Csize_t}(work_dim)
        for (i, o) in enumerate(global_work_offset)
            goffset[i] = o
        end
    end

    lsize = C_NULL
    if local_work_size !== nothing
        if length(local_work_size) > max_work_dim
            throw(ArgumentError("local_work_offset has max dim of $max_work_dim"))
        end
        if length(local_work_size) != work_dim
            throw(ArgumentError("global_work_size and local_work_size have differing dims"))
        end
        lsize = Array{Csize_t}(work_dim)
        for (i, s) in enumerate(local_work_size)
            lsize[i] = s
        end
    end

    if wait_on !== nothing
        n_events = cl_uint(length(wait_on))
        wait_event_ids = [evt.id for evt in wait_on]
    else
        n_events = cl_uint(0)
        wait_event_ids = C_NULL
    end

    ret_event = Ref{CL_event}()
    @check api.clEnqueueNDRangeKernel(q.id, k.id, cl_uint(work_dim), goffset, gsize, lsize,
                                      n_events, wait_event_ids, ret_event)
    return Event(ret_event[], retain=false)
end

function enqueue_task(q::CmdQueue, k::Kernel; wait_for=nothing)
    n_evts  = 0
    evt_ids = C_NULL
    #TODO: this should be split out into its own function
    if wait_for !== nothing
        if isa(wait_for, Event)
            n_evts = 1
            evt_ids = [wait_for.id]
        else
            @assert all([isa(evt, Event) for evt in wait_for])
            n_evts = length(wait_for)
            evt_ids = [evt.id for evt in wait_for]
        end
    end
    ret_event = Ref{CL_event}()
    @check api.clEnqueueTask(q.id, k.id, n_evts, evt_ids, ret_event)
    return ret_event[]
end

let name(k::Kernel) = begin
        size = Ref{Csize_t}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_FUNCTION_NAME,
                                   0, C_NULL, size)
        result = Vector{Cchar}(size[])
        @check api.clGetKernelInfo(k.id, CL_KERNEL_FUNCTION_NAME,
                                   size[], result, size)
        return CLString(result)
    end

    num_args(k::Kernel) = begin
        ret = Ref{CL_uint}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_NUM_ARGS,
                                   sizeof(CL_uint), ret, C_NULL)
        return ret[]
    end

    reference_count(k::Kernel) = begin
        ret = Ref{CL_uint}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_REFERENCE_COUNT,
                                   sizeof(CL_uint), ret, C_NULL)
        return ret[]
    end

    program(k::Kernel) = begin
        ret = Ref{CL_program}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_PROGRAM,
                                   sizeof(CL_program), ret, C_NULL)
        return Program(ret[], retain=true)
    end

    attributes(k::Kernel) = begin
        size = Ref{Csize_t}()
        api.clGetKernelInfo(k.id, CL_KERNEL_ATTRIBUTES,
                            0, C_NULL, size)
        if size[] <= 1
            return ""
        end
        result = Vector{CL_char}(size[])
        @check api.clGetKernelInfo(k.id, CL_KERNEL_ATTRIBUTES,
                                   size[], result, size)
        return CLString(result)
    end

    info_map = Dict{Symbol, Function}(
        :name => name,
        :num_args => num_args,
        :reference_count => reference_count,
        :program => program,
        :attributes => attributes
    )

    function info(k::Kernel, kinfo::Symbol)
        try
            func = info_map[kinfo]
            func(k)
        catch err
            if isa(err, KeyError)
                error("OpenCL.Kernel has no info for: $kinfo")
            else
                throw(err)
            end
        end
    end
end

#TODO set_arg sampler...
# OpenCL 1.2 function
#TODO: get_arg_info(k::Kernel, idx, param)
#TODO: enqueue_async_kernel()


text:
module CachedValues

export CachedValue

mutable struct CachedValue{F}
  value::F
end

end # module

ext: jl


text:
let
    mktemp() do path, _

        a=1
        b=2.0
        c=[1.0,2.0]
        d=rand(ComplexF64,100,100)
        e=Dict(s => 2 for s in ["a", "b", "c"])
        f="A String"
        g=ntuple(i->i^2, 5)
        h=rand(2000)

        # jldsave
        jldsave(path; a, b, c, d, e, f, g, h)

        jldopen(path) do f
            for k in keys(f)
                f[k]
            end
        end

        nothing
    end
end


text:
module Sensors
export humidity, temperature, pressure

const I2C_DEVICE_PATH = "/dev/i2c-1"
const I2C_DEVICE = Ref{IOStream}()
const I2C_ADDR = Ref{UInt8}(0xff)

function __init__()
    try
        I2C_DEVICE[] = open("/dev/i2c-1","r+")
        HTS221_calibrate()
    catch
        @warn("I2C: Sense Hat not found")
    end
end

import ..ioctl
include("smbus.jl")

const LPS25H_ADDRESS = 0x5c # Pressure/temp
const HTS221_ADDRESS = 0x5f # Humidity/temp

function setaddr(addr::UInt8)
    if I2C_ADDR[] != addr
        ioctl(I2C_DEVICE[], I2C_SLAVE, addr)
        I2C_ADDR[] = addr
    end
end

include("HTS221.jl")
include("LPS25H.jl")

"""
    humidity()

The relative humidity (as a percentage between 0 and 100).
"""
humidity()    = HTS221_humidity()

"""
    temperature()

The temperature (in °C).
"""
temperature() = HTS221_temperature()

"""
    pressure()

The atmospheric pressure (in millibars).
"""
pressure()    = LPS25H_pressure()

end # module


text:
"""
    reweight(p::StaticParticle{T}, w::Real, logweight = false) where {T<:Any}

Reweight a particle - update the weight.

- ```p``` Particle to update
- ```w``` New weight of particle
- ```logweight``` Is ```w``` on log scale? Default is ```false```

"""
function reweight(p::StaticParticle{T}, w::Real, logweight = false) where {T<:Any}

    Particle(p.x, w, logweight)

end

# function reweight(p::DynamicParticle{T}, τ::Real) where {T<:Any}
#     if p.τ == τ
#         error("Current particle and  proposal particle have same temperatures τ")
#     end
#
#     newlogdensval = p.logdensity(p.x, τ)
#
#     logweight = newlogdensval -
#         p.cacheℓdens +
#         p.ℓw
#
#     DynamicParticle(p.x, logweight, p.ℓdens, τ, newlogdensval)
#
# end


text:
import LinearAlgebra
import ILog2

export flipbit, control
export swap, proj, ket, ketbra, zeroket, oneket
export ⊗

const ⊗ = kron

# FIXME: I think I copied this from somewhere. Probably get rid of it.
"""
    flipbit(x)

Return `1` if `x` is `0` and vice versa.
"""
flipbit(x) = x == 1 ? zero(x) : one(x)

"""
    swap(G::AbstractMatrix)

Swap the first and second qubits in the two-qubit operator `G`.
"""
swap(G) = SWAP * G * SWAP

"""
    ket(b::Integer)
    ket(bs::Integer...)
    ket(s::AbstractString)

Return a vector representation the single-qubit computational basis state indexed by `b`, either `0` or `1`.
This returns either `z0` or `z1`.

Return the `length(bs)`-qubit computational basis state corresponding to bits `bs`.
Each argument must be either `0` or `1`.

Construct a computational-basis ket from bitstring `s`.
"""
ket(b::Integer) = b == 0 ? z0 : b == 1 ? z1 : throw(DomainError(b))

function _ket(position::Integer, n_qubits::Integer)
    v = zeros(Int, 2^n_qubits)
    v[position] = 1
    return v
end

ket(b0::Integer, bs::Integer...) = _ket_from_int_args(b0, bs...)

function _ket_from_int_args(bs::Integer...)
    pos = 1
    n = length(bs)
    for (i, val) in enumerate(bs)
        pos += val << (n-i)
    end
    return _ket(pos, length(bs))
end

# Example: ket("01101")
function ket(s::AbstractString)
    pos = parse(Int, s, base=2) + 1
    n_qubits = length(s)
    return _ket(pos, n_qubits)
end

"""
    zeroket(n_qubits)

Return the `n_qubits`-qubit zero ket.
"""
zeroket(n_qubits::Integer) = _ket(1, n_qubits)

"""
    oneket(n_qubits)

Return the `n_qubits`-qubit one ket.
"""
oneket(n_qubits::Integer) = _ket(2^n_qubits, n_qubits)

"""
    ketbra(ket::AbstractVector, bra::LinearAlgebra.Adjoint)

Return `ket * bra`. That is the matrix element `|ket><bra|`.
"""
ketbra(ket::AbstractVector, bra::LinearAlgebra.Adjoint) = ket * bra

"""
    ketbra(ket1::AbstractVector, ket2::AbstractVector)

Return `ket1 * ket2'`
"""
ketbra(ket1::AbstractVector, ket2::AbstractVector) = ketbra(ket1, ket2')

ketbra(ket1::Tuple, ket2::Tuple) = ketbra(ket(ket1...), ket(ket2...))

"""
    ketbra(b0::Integer, b1::Integer)

Return `ketbra(ket(b0), ket(b1))`
"""
ketbra(b0::Integer, b1::Integer) = ketbra(ket(b0), ket(b1))

# Use a materialized covector. This is 8 or so times slower than using the adjoint wrapper (On Julia 1.2)
ketbra(ket1::AbstractVector, ket2::AbstractMatrix) = ket1 * ket2 # Assume dimensions are correct. i.e. ket2 is an adjoint

"""
    proj(ket::AbstractVector)

Return `ketbra(ket, ket)`.
"""
proj(ket::AbstractVector) = ketbra(ket, ket)

"""
    proj(ket::Integer)

Return `proj(ket(b))`.
"""
proj(b::Integer) = proj(ket(b))

"""
    control(A::AbstractMatrix, n_control_bits=1)

Return a controlled-`A` gate with `n_control_bits` control bits.
"""
function control(A::AbstractMatrix, n_control_bits=1)
    T = eltype(A)
    dA = checksquare(A)
    nqubits = ILog2.ilog2(dA) + n_control_bits
    d = 2^nqubits
    controlled_A = zeros(T, (d, d))
    @inbounds for i in 1:(d-dA)
        controlled_A[i, i] = one(T)
    end
    @inbounds controlled_A[d-dA+1:d, d-dA+1:d] .= A
    return controlled_A
end


text:
# 1- read in text file:
str = read("words1.txt", String)
# println(str) #> to be or not to be that is the question
# 2- replace non alphabet characters and digits from text with a space:
nonalpha = r"(\W\s?)"
str = replace(str, nonalpha => ' ')
digits = r"(\d+)"
str = replace(str, digits => ' ')
#> "to be or not to be that is the question "
# 3- split text in words:
word_list = split(str, ' ')
# println(word_list)
#> "to" "be" "or" "not" "to" "be" "that" "is" "the" "question" ""
# 4- make a dictionary with the words and count their frequencies:
word_freq = Dict{String, Int64}()
for word in word_list
    word = strip(word)
    if isempty(word) continue end
    haskey(word_freq, word) ?
    	word_freq[word] += 1 :
    	word_freq[word] = 1
end
# 5- sort the words (the keys) and print out the frequencies:
println("Word : frequency \n")
words = sort!(collect(keys(word_freq)))
for word in words
    println("$word : $(word_freq[word])")
end


text:
# using Revise
using Test
using MLJRegistry
using MLJBase

#@test !(isempty(MLJRegistry.metadata()))

d = Dict()
d[:test] = Tuple{Union{Continuous,Missing},Discrete}
d["junk"] = Dict{Any,Any}("H" => Missing, :cross => "lemon", :t => :w, "r" => "r")
d["a"] = "b"
d[:f] = true
d["j"] = :post

@test MLJRegistry.decode_dic(MLJRegistry.encode_dic(d)) == d


text:
using Test
using EquationEditor
const EE = EquationEditor
@test length(EquationEditor.MIMETYPES) == 985
reqstr = "svg\\Mickey%20mouse.svg/"
reqstr = reqstr |> EE._bck2fwdslash
reqstr = reqstr |> EE._unescapeuri
@test "hh%20%za" |> EE._unescapeuri == ""
reqstr = reqstr |> EE._nobadchars
reqstr = reqstr |> EE.lowercase
reqstr = reqstr |> EE._nojumpup
reqstr = reqstr |> EE._nokeywords
reqstr = reqstr |> EE._maybeindexhtml
reqstr = reqstr |> EE._doubleslash2slash
@test reqstr |> EE._startslash == "/svg/mickey mouse.svg/index.html"
@test reqstr |> EE._startslash == "/svg/mickey mouse.svg/index.html"
# Note: Query strings should be handled elsewhere. This file handler will just return a 404 file not found.
reqstr = "/fonts/fa.woff?ver=4.7.0;duck=bird&cat=feline"
@test  reqstr |> EE._noschemechars == "/fonts/fa.woff?ver4.7.0duckbirdcatfeline"
@test  reqstr |> EE._censoredrequest == "/fonts/fa.woffver4.7.0duckbirdcatfeline"
@test EE._censoredrequest("/svg/../others/") == "/svg/others/index.html"
@test EE._censoredrequest("../fonts/others/") == "/fonts/others/index.html"
@test EE._censoredrequest("/svg/Mickey%20%mouse.svg/") == "/index.html"
@test EE._censoredrequest("/html/touch-position/index.mjs") == "/html/touch-position/index.mjs"
@test EE._censoredrequest("mehekje.dot.no") == "/mehekje.dot.no"

show(EE._censoredrequest("/fonts/fa.woff?ver=4.7.0;duck=bird&cat=feline"))
show(EE._noschemechars("?A:;!"))
@test EE._censoredrequest("/svg/") == "/svg/index.html"

@test EE._longfilename("XX") == EE.FOLDERHOME * "XX"
@test EE._longfilename("/svg/Julia.svg") |> EE._pubfilename == "/svg/Julia.svg"
@test EE._longfilename("/fonts/fa.woff?ver=4.7.0") |> EE._pubfilename == "/fonts/fa.woff?ver=4.7.0"
@test EE._locfolder("/svg/special/Julia.svg") == "/svg/"
@test EE._locfolder("/stylesheet.css") == "/css/"
@test EE._locfolder("/stylesheet.woff") == "/html/"
@test EE._locfolder("/fonts/stylesheet.woff") == "/fonts/"
@test EE._locfolder("/jsmodule.mjs") == "/mjs/"
@test EE._locfolder("/gl-vec2/copy.mjs") == "/mjs/"
@test EE._locfolder("/html/touch-position/index.mjs") == "/mjs/"
@test EE._locfolder("/svg/index.html") == "/svg/"
@test EE._locfolder("/fonts/fa.woff?ver=4.7.0") == "/fonts/"

@test EE._topfolder_modreq("/Mickey.jpg") == ("/img/", "/img/mickey.jpg")
@test EE._topfolder_modreq("/img/Mickey.jpg") == ("/img/", "/img/mickey.jpg")
@test EE._topfolder_modreq("/html/touch-position/index.mjs") == ("/mjs/", "/mjs/touch-position/index.mjs")
@test EE._topfolder_modreq("/sub/touch-position/index.mjs") == ("/mjs/", "/mjs/sub/touch-position/index.mjs")
@test EE._topfolder_modreq("/sub/touch-position/index.mjs?ver=4.0.2") == ("/html/", "/html/sub/touch-position/index.mjsver4.0.2")
@test EE._topfolder_modreq("/sub/") == ("/html/", "/html/sub/index.html")
@test EE._topfolder_modreq("/html/index.html") == ("/html/", "/html/index.html")
@test EE._topfolder_modreq("/html/html/index.html") == ("/html/", "/html/html/index.html")
@test EE._topfolder_modreq("/svg/") == ("/svg/", "/svg/index.html")
@test EE._topfolder_modreq("/fonts/") == ("/fonts/", "/fonts/index.html")

pubstr = "/html/favicon.ico"
tf = EE.FOLDERHOME * pubstr
topfolder, modreq = pubstr |>  EE._topfolder_modreq
@test tf == EE._longfilename(modreq)
@test EE._pubfilename(tf) == pubstr

status, body, mime, tstamp = EE._rawresponse(pubstr)
skeleton = EE.WebSockets.Response()
skeleton.body = body
EE.WebSockets.setheader(skeleton, EE.WebSockets.Header("Content-Type", mime))
EE.WebSockets.setheader(skeleton, EE.WebSockets.Header("Content-Length",  string(length(body))))
EE.WebSockets.setheader(skeleton, EE.WebSockets.Header("Last-Modified", tstamp))
skeleton.status = status

pubindex = "/html/index.html"
@test EE._topfolder_modreq("") == ("/html/", pubindex)
@test EE._topfolder_modreq("/") == ("/html/", pubindex)
longindex = splitdir(EE.FOLDERHOME * pubindex)[1]
@test EE._pubfilename(longindex) == "/html"

@test typeof(EE._htmlindex(longindex)) == String

skeleton = EE.WebSockets.Response()
req = EE.WebSockets.Request("get", "/")
response = EE._handle_file(req, skeleton)
response = EE.handle_file(EE.WebSockets.Request("get", "/"))
req = EE.WebSockets.Request("get", "/favicon.ico")
response = EE.handle_file(req)
@test response.status == 200
response = EE.handle_file(EE.WebSockets.Request("get", "/favicon.ico"))
@test response.status == 200

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
# ==================================================================================================================
function solveStokes_penalty( U, P,  θ, r, Ucartesian,Upolar, g, ρ, η, 𝓒,
                      theta, rr, TT,
                      ∂Ωu, ufix, ifree,
                      PhaseID, EL2NOD, EL2NODP,
                      KKidx, GGidx, invMMidx,
                      to)
    @timeit to "Stokes" begin
        @timeit to "Assembly" KK, GG, invMM, Rhs =
            assembly_stokes_penalty(EL2NOD, EL2NODP, theta, rr, g, ρ, η, 𝓒, PhaseID, KKidx, GGidx, invMMidx)
        @timeit to "BCs" begin
            KK, GG, Rhs = _prepare_matrices_penalty(KK, GG, Rhs, TT)
            U, Rhs = _apply_bcs(U,KK,Rhs,∂Ωu,ufix)
        end
        @timeit to "Solve" U, P = PowellHesteness(U, P, KK, invMM, GG, Rhs, ifree)
        @timeit to "Remove net rotation" U, Ucart, Upolar, Ucartesian = updatevelocity(U,Ucartesian, Upolar,  θ, r,TT)
        # @timeit to "Remove nullspace"  U, Ucart, Upolar, Ucartesian = updatevelocity2(U, Ucartesian, Upolar, ρ, EL2NOD, TT, theta, rr, r)

    end
    return Ucartesian,Upolar,U,Ucart,P,to
end
# ==================================================================================================================

@inline function _apply_bcs_penalty(U, KK, Rhs, ∂Ω, vfix)
    U[∂Ω] = vfix # write prescribed temperature values
    Rhs -= KK[:,∂Ω] * vfix
    return U,Rhs
end

@inline function _prepare_matrices_penalty(KK, GG, Fb,TT)

    dropzeros!(GG)
    KK = SparseMatrixCSC(Symmetric(KK,:L))
    dropzeros!(KK)
    transTT = SparseMatrixCSC(TT')
    KK = transTT*KK*TT
    GG = transTT*GG
    Fb = transTT*Fb
    KK = SparseMatrixCSC(Symmetric(KK,:L))
    return KK,GG,Fb
end

function assembly_stokes_penalty(EL2NOD, EL2NODP, theta, r, g, ρ, η, 𝓒, PhaseID, KKidx, GGidx, invMMidx)
    penalty = 1e6
    # ============================================ MODEL AND BLOCKING PARAMETERS
    ndim = 2
    # nvert = 3
    # EL2NOD0 = deepcopy(EL2NOD)
    nnodel = size(EL2NOD,1)
    nel = size(EL2NOD,2)
    nelblk = min(nel, 1200)
    nblk = ceil(Int,nel/nelblk)
    il = one(nelblk); iu = nelblk
    # =========== PREPARE INTEGRATION POINTS & DERIVATIVES wrt LOCAL COORDINATES
    nip = 6
    # nUdof = ndim*maximum(EL2NOD)
    # nPdof = maximum(EL2NODP)
    nUdofel = ndim * nnodel
    nPdofel = size(EL2NODP,1)
    ni, nn, nnP, nn3 = Val(nip), Val(nnodel), Val(3), Val(3)
    N,dNds,_,w_ip = _get_SF(ni,nn)
    NP,_,_,_ = _get_SF(ni,nnP)
    N3,_,dN3ds,_ = _get_SF(ni,nn3)
    # ============================================================= ALLOCATIONS
    detJ_PL = Vector{Float64}(undef,nelblk)
    R_31 = similar(detJ_PL) # =  detJa_PL*dxi_dth
    R_21 = similar(detJ_PL) # = -detJa_PL*deta_dth
    Th_31 = similar(detJ_PL) # = -detJa_PL*dxi_dr
    Th_21 = similar(detJ_PL) # =  detJa_PL*deta_dr
    # NxN = Array{Float64,2}(undef,3,3)
    # ang_dist        = Array{Float64,2}(undef,3,nelblk)
    # VCOORD_th       = similar(ang_dist)
    # VCOORD_r        = similar(ang_dist)
    dNdx = Array{Float64,2}(undef,nelblk,size(dNds[1],2))
    dNdy = similar(dNdx)
    ω = Vector{Float64}(undef,Int(nelblk))
    invJx_double = Array{Float64,2}(undef,Int(nelblk),ndim)    # storage for x-components of Jacobi matrix
    invJz_double = similar(invJx_double)                     # storage for z-components of Jacobi matrix
    # ==================================== STORAGE FOR DATA OF ELEMENTS IN BLOCK
    K_blk = fill(0.0, nelblk, Int(nUdofel*(nUdofel+1)/2))
        # symmetric stiffness matrix, dim=vel dofs, but only upper triangle
    M_blk = fill(0.0, nelblk, Int(nPdofel*(nPdofel+1)/2))
        # symmetric mass matrix, dim=pressure nodes, but only upper triangle
    G_blk = fill(0.0, nelblk, nPdofel*nUdofel)
        # asymmetric gradient matrix, vel dofs x pressure nodes
    Fb_blk = fill(0.0, nelblk, nUdofel)
        # storage for entries in bouyancy force vector
    invM_blk = fill(0.0, nelblk, nPdofel*nPdofel)
    invMG_blk = fill(0.0, nelblk, nPdofel*nUdofel)

    # ====================================  STORAGE FOR DATA OF ALL ELEMENT MATRICES/VECTORS
    K_all = Array{Float64,2}(undef, Int(nUdofel*(nUdofel+1)/2),nel)
    M_all = Array{Float64,2}(undef, Int(nPdofel*(nPdofel+1)/2),nel)
    G_all = Array{Float64,2}(undef, nUdofel*nPdofel,nel)
    Fb_all = Array{Float64,2}(undef, nUdofel,nel)
    invM_all = Array{Float64,2}(undef, nPdofel*nPdofel,nel)

    #=========================================================================
    BLOCK LOOP - MATRIX COMPUTATION
    =========================================================================#
    for ib in 1:nblk
        #===========================================================================
        CALCULATE JACOBIAN, ITS DETERMINANT AND INVERSE
        ===========================================================================#
        """
        NOTE: For triangular elements with non-curved edges the Jacobian is
              the same for all integration points (i.e. calculated once
              before the integration loop). Further, linear 3-node shape are
              sufficient to calculate the Jacobian.
        """
        # idx       = @views EL2NOD[1:nvert,il:iu]
        # VCOORD_th = reshape(theta[idx],nvert,nelblk)
        # VCOORD_r  = reshape(r[idx],nvert,nelblk)
        VCOORD_th = view(theta,:,il:iu)
        VCOORD_r = view(r,:,il:iu)

        J_th    = (VCOORD_th'*dN3ds')
        J_r     = (VCOORD_r'*dN3ds')

        _fill_R_J!(J_th,J_r,VCOORD_th,VCOORD_r,
            R_31,R_21,Th_31,Th_21,detJ_PL)

        # ------------------------ NUMERICAL INTEGRATION LOOP (GAUSS QUADRATURE)
        _ip_loop_penalty!(K_blk, M_blk, G_blk, Fb_blk,
                  g, η, 𝓒, ρ, EL2NOD, PhaseID,il:iu,
                  dNds, w_ip, nip, nnodel, nelblk, nPdofel,nUdofel,
                  dNdx,dNdy,N,N3,NP,
                  ω,invJx_double, invJz_double, detJ_PL,
                  R_21,R_31,Th_21,Th_31,
                  VCOORD_th,VCOORD_r)

        _invM!(invM_blk, M_blk, detJ_PL)

        _fill_invMxG!(invMG_blk, invM_blk, G_blk, nPdofel, nUdofel)

        _fill_KplusGxinvMxGt(K_blk, invMG_blk, G_blk, nPdofel, nUdofel, penalty)

        # --------------------------------------- WRITE DATA INTO GLOBAL STORAGE
        @views begin
            K_all[:,il:iu] .= K_blk'
            G_all[:,il:iu] .= G_blk'
            M_all[:,il:iu] .= M_blk'
            Fb_all[:,il:iu] .= Fb_blk'
            invM_all[:,il:iu] .= invM_blk'
        end

        # -------------------------------------------------- RESET BLOCK MATRICES
        fill!(K_blk,0.0)
        fill!(G_blk,0.0)
        fill!(M_blk,0.0)
        fill!(Fb_blk,0.0)
        fill!(invM_blk,0.0)
        fill!(invMG_blk,0.0)

        # -------------------------------- READJUST START, END AND SIZE OF BLOCK
        il += nelblk;
        if ib == nblk-1
           # ------------ Account for different number of elements in last block
           nelblk = nel - il + 1  # new block size
           # --------------------------------- Reallocate at blocks at the edge
           detJ_PL = Vector{Float64}(undef,nelblk)
           R_31 = similar(detJ_PL)
           R_21 = similar(detJ_PL)
           Th_31 = similar(detJ_PL)
           Th_21 = similar(detJ_PL)
        #    ang_dist = Array{Float64,2}(undef,3,nelblk)
           dNdx = Array{Float64,2}(undef,nelblk,size(dNds[1],2))
           dNdy = similar(dNdx)
           K_blk = fill(0.0, nelblk, Int(nUdofel*(nUdofel+1)/2))
           M_blk = fill(0.0, nelblk, Int(nPdofel*(nPdofel+1)/2))
           G_blk = fill(0.0, nelblk, nPdofel*nUdofel)
           Fb_blk = fill(0.0, nelblk, nUdofel)
           ω = Vector{Float64}(undef,Int(nelblk))
           invJx_double = Array{Float64,2}(undef,Int(nelblk),ndim)
           invJz_double = similar(invJx_double)
           invM_blk = fill(0.0, nelblk, nPdofel*nPdofel)
           invMG_blk = fill(0.0, nelblk, nPdofel*nUdofel)
        end
        iu += nelblk

    end # end block loop

    # #===========================================================================
    # ASSEMBLY OF GLOBAL SPARSE MATRICES AND RHS-VECTOR
    # ===========================================================================#
    KK, GG, invMM, Fb =
        _create_stokes_matrices_penalty(EL2NOD, nUdofel, ndim,
                                        KKidx, GGidx, invMMidx,
                                        K_all, G_all, invM_all, Fb_all)

    return KK, GG, invMM, Fb

end # END OF ASSEMBLY FUNCTION

function _ip_loop_penalty!( K_blk, M_blk, G_blk, Fb_blk,
                            g, η, 𝓒, ρ, EL2NOD, PhaseID,els,
                            dNds, w_ip, nip, nnodel, nelblk, nPdofel,nUdofel,
                            dNdx,dNdy,N,N3,NP,
                            ω,invJx_double, invJz_double, detJ_PL,
                            R_21,R_31,Th_21,Th_31,
                            VCOORD_th,VCOORD_r)

    sin_ip  = similar(detJ_PL)
    cos_ip  = similar(detJ_PL)

    for ip=1:nip

        N_ip = N3[ip]
        # NP_blk = repeat(NP[ip],nelblk,1)
        NP_blk = NP[ip].*ones(nelblk)

        #=======================================================================
        PROPERTIES OF ELEMENTS AT ip-TH EVALUATION POINT
        =======================================================================#
        Dens_blk    = _element_density(ρ,EL2NOD,PhaseID,els,NP[ip])
        Visc_blk    = _element_viscosity(η,EL2NOD,PhaseID,els,NP[ip])
        #  Gravitational force at ip-th integration point
        Fg_blk      = g * Dens_blk

        #==========================================================================================================
        CALCULATE 2nd JACOBIAN (FROM CARTESIAN TO POLAR COORDINATES --> curved edges), ITS DETERMINANT AND INVERSE
        ===========================================================================================================
        NOTE: For triangular elements with curved edges the Jacobian needs to be computed at each integration
               point (inside the integration loop). =#
        th_ip   = VCOORD_th'*N_ip'
        r_ip    = VCOORD_r'*N_ip'

        _derivative_weights!(dNds[ip],ω,dNdx,dNdy,w_ip[ip],th_ip,r_ip,sin_ip,cos_ip,
            R_21,R_31,Th_21,Th_31,detJ_PL,invJx_double,invJz_double)

        _fill_Kblk!(K_blk, 𝓒, Visc_blk, Val(η), ω, dNdx, dNdy, nnodel, els, ip)

        _fill_Gblk!(G_blk, NP_blk, ω, dNdx, dNdy, nPdofel,nUdofel,nnodel)

        _fill_Mblk_penalty!(M_blk, NP_blk, ω, nPdofel)

        _fill_Fbblk!(Fb_blk,Fg_blk, N[ip],
            ω, sin_ip,cos_ip,nUdofel)

    end # end integration point loop
end # ENF OF IP_LOOP FUNCTION

@inline function _fill_Mblk_penalty!(M_blk, NP_blk, ω, nPdofel)
    # "weight" divided by viscosity at integration point in all elements
    indx = 1;
    for i in 1:nPdofel
        @inbounds @simd for j in i:nPdofel
            M_blk[:,indx] .+=  ω.*view(NP_blk,:,i).*view(NP_blk,:,j)
            indx += 1
        end
    end
end

function _invM!(invM_blk, M_blk, detJ)
    TMP   = 1.0./detJ
    M_blk = M_blk .* TMP
    # M_blk = M_blk .* TMP[:,ones[1,nPdofel*[nPdofel+1]/2]];

    # How to calculate the determinante of several symmetric 3x3 matrices:
    # det[M] =   M11*M22*M33 + M12*M23*M31 + M13*M21*M32
    #          - M13*M22*M31 - M12*M21*M33 - M11*M23*M32
    #   written in 1-index notation:
    #        =   M1 *M5 *M9  + M4 *M8 *M3  + M7 *M2 *M6
    #          - M7 *M5 *M3  - M4 *M2 *M9  - M1 *M8 *M6
    #   considering symmetry and using lower triangular part only
    #   [M4=M2, M7=M3, M8=M6]
    #        =   M1 *M5 *M9  + M2 *M6 *M3  + M3 *M2 *M6
    #          - M3 *M5 *M3  - M2 *M2 *M9  - M1 *M6 *M6
    #   and knowing where the 6 different values are stored in M_blk
    #   [1-->M1, 2-->M2, 3-->M3, 4-->M5, 5-->M6, 6-->M9]
    #        =   M_blk1*M_blk4*M_blk6 + M_blk2*M_blk5*M_blk3 + M_blk3*M_blk2*M_blk5
    #          - M_blk3*M_blk4*M_blk3 - M_blk2*M_blk2*M_blk6 - M_blk1*M_blk5*M_blk5
    #   re-arranging
    #        =   M_blk1 * [M_blk4*M_blk6 - M_blk5*M_blk5]
    #          + M_blk2 * [M_blk5*M_blk3 - M_blk2*M_blk6]
    #          + M_blk3 * [M_blk2*M_blk5 - M_blk4*M_blk3]
    detM_blk = @. M_blk[:,1] * (M_blk[:,4]*M_blk[:,6] - M_blk[:,5]*M_blk[:,5]) +
                  M_blk[:,2] * (M_blk[:,5]*M_blk[:,3] - M_blk[:,2]*M_blk[:,6]) +
                  M_blk[:,3] * (M_blk[:,2]*M_blk[:,5] - M_blk[:,4]*M_blk[:,3])
    detM_blk ./=TMP

    # The determinante is used to calculate the inverse of the symmetric
    # 3x3 element mass matrices. The same logic as above is used.
    invM_blk[:,1] = @views @. (M_blk[:,4]*M_blk[:,6] - M_blk[:,5]*M_blk[:,5])/detM_blk;
    invM_blk[:,2] = @views @. (M_blk[:,5]*M_blk[:,3] - M_blk[:,2]*M_blk[:,6])/detM_blk;
    invM_blk[:,3] = @views @. (M_blk[:,2]*M_blk[:,5] - M_blk[:,4]*M_blk[:,3])/detM_blk;
    invM_blk[:,4] = @views invM_blk[:,2];
    invM_blk[:,5] = @views @. (M_blk[:,1]*M_blk[:,6] - M_blk[:,3]*M_blk[:,3])/detM_blk;
    invM_blk[:,6] = @views @. (M_blk[:,2]*M_blk[:,3] - M_blk[:,1]*M_blk[:,5])/detM_blk;
    invM_blk[:,7] = @views invM_blk[:,3];
    invM_blk[:,8] = @views invM_blk[:,6];
    invM_blk[:,9] = @views @. (M_blk[:,1]*M_blk[:,4] - M_blk[:,5]*M_blk[:,5])/detM_blk;

end

function _fill_invMxG!(invMG_blk, invM_blk, G_blk, nPdofel, nUdofel)
    # --------------------------invM*G'----------------------------
    for i=1:nPdofel
        for j=1:nUdofel
            for k=1:nPdofel
                invMG_blk[:,(i-1)*nUdofel+j] = @views invMG_blk[:,(i-1)*nUdofel+j] +
                    invM_blk[:,(i-1)*nPdofel+k].*G_blk[:,(k-1)*nUdofel+j];
            end
        end
    end
end

function _fill_KplusGxinvMxGt(K_blk, invMG_blk, G_blk, nPdofel, nUdofel, penalty)
    # -----------------K = K + penalty*G*invM*G'-------------------
    indx = 1;
    for i=1:nUdofel
        for j=i:nUdofel
            for k=1:nPdofel
                K_blk[:,indx] = @views K_blk[:,indx] +
                     penalty*G_blk[:,(k-1)*nUdofel+i].*invMG_blk[:,(k-1)*nUdofel+j];
            end
            indx += 1;
        end
    end
end

@inline function _create_stokes_matrices_penalty(EL2NOD,nUdofel,ndim,
    KKidx, GGidx, invMMidx,
    K_all, G_all, invM_all, Fb_all)

    nel     = size(EL2NOD,2)

    #==========================================================================
    ASSEMBLY OF GLOBAL STIFFNESS MATRIX
    ==========================================================================#
    #  convert triplet data to sparse global stiffness matrix (assembly)
    KK = sparse(KKidx._i ,KKidx._j, vec(K_all))

    #==========================================================================
    ASSEMBLY OF GLOBAL GRADIENT MATRIX
    ==========================================================================#
    # convert triplet data to sparse global gradient matrix (assembly)
    GG = sparse(GGidx._i ,GGidx._j ,vec(G_all))

    #==========================================================================
    ASSEMBLY OF GLOBAL FORCE VECTOR
    ==========================================================================#
    EL2DOF = Array{Int64}(undef,nUdofel, nel)
    @views EL2DOF[1:ndim:nUdofel,:] .= @. ndim*(EL2NOD-1)+1
    @views EL2DOF[2:ndim:nUdofel,:] .= @. ndim*(EL2NOD-1)+2
    Fb = accumarray(vec(EL2DOF),vec(Fb_all));

    #==========================================================================
    ASSEMBLY OF GLOBAL (1/VISCOSITY)-SCALED MASS MATRIX
    ==========================================================================#
    # convert triplet data to sparse global matrix
    invMM = sparse(invMMidx._i ,invMMidx._j , vec(invM_all))

    return KK, GG, invMM, Fb

end

# POWELL-HESTENESS SOLVER ================================================================
@inline  function PowellHesteness(U, P, KK, M⁻¹, GG, Rhs, ifree)

    itmax_PH = 50
    itmin_PH = 2
    # itnum = 0
    norm∇U = Vector{Float64}(undef, itmax_PH)

    Rhs .+= GG*P

    LL = factorize(KK[ifree, ifree])
    # LL = cholesky(KK[ifree, ifree])
    GGᵀ = GG'
    # Uifree = U[ifree]
    # Rhsifree = Rhs[ifree]
    penalty = 1e6

    for itPH = 1:itmax_PH
        @inbounds U[ifree] .= LL\Rhs[ifree]
        ∇U = -GGᵀ*U
        # rm0space!(∇U)
        ΔP = penalty*M⁻¹*∇U
        Rhs .+= GG*ΔP
        P .+= ΔP
        # Check convergence -----------------------------------------------
        norm∇U[itPH] = mynorm(∇U)

        if (itPH > itmin_PH) && (norm∇U[itPH]/norm∇U[itPH-1] > 0.1)
            println("\n", itPH," Powell-Hesteness iterations\n")
            break
        end

    end

    # Ux = U[1:2:end]
    # Uz = U[2:2:end]
    # M = sum(mean(ρ[EL2NOD],dims=1)'.*area_el)
    # Ax =  sum(mean(ρ[EL2NOD].*Ux[EL2NOD], dims=1)'.*area_el)
    # Az =  sum(mean(ρ[EL2NOD].*Uz[EL2NOD], dims=1)'.*area_el)
    # A = @. (Ax+Az)/M
    # U .-=A

    return U, P

end # END POWELL-HESTENESS SOLVER

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
using QuantumSDPs
using Test
using SCS, Convex
using Random
using LinearAlgebra, SparseArrays

if VERSION < v"1.2pre"
    import LinearAlgebra.I
    (I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))
end
const ⊗ = kron
ket(i, d) = sparsevec([i], [true], d)
bra(i, d) = ket(i, d)'
ketbra(i,j, d) = ket(i,d) ⊗ bra(j, d)

@testset "QuantumSDPs.jl" begin

    @testset "permuteddims_matrix" begin
        for n in (2, 3, 4, 5)
            dims = ntuple( i -> rand(2:5), n)
            d = prod(dims)
            v = rand(d)
            p = randperm(n)
            out1 = vec(permutedims(reshape(v, dims), p))
            out2 = QuantumSDPs.permutedims_matrix(dims, p) * v
            @test out1 ≈ out2
        end
    end

    @testset "choi2so" begin
        for (dA, dB) in ((2,2), (2,3), (3,2), (3,3), (3, 7), (7, 3))
            J = Choi(dA, dB)
            T = ComplexVariable(dB^2, dA^2)
            p = minimize(real(tr( rand(dA^2, dB^2) * T)), [ choi2so(J, dA, dB) == T ])
            solve!( p, SCSSolver(verbose=false))

            # check conic_form agrees with the numeric implementation
            @test evaluate(T) ≈ choi2so(evaluate(J), dA, dB) atol=1e-4

            # check we can map back
            @test so2choi(evaluate(T), dA, dB) ≈ evaluate(J) atol=1e-4

            # check that `T` is unital
            @test evaluate(T)' * vec(I(dB)) ≈ vec(I(dA)) atol=1e-4

            ρ = randn(dA,dA) + im*randn(dA,dA)
            ρ = ρ/tr(ρ)

            for R in (evaluate(T), choi2so(evaluate(J), dA, dB))
                # define a function to implement of the action
                # of R on a matrix
                f_R = x ->  reshape(R*vec(x), (dB, dB))

                # check that the map is trace preserving
                @test tr(f_R(ρ)) ≈ tr(ρ) atol=1e-4

                # check that the method for using J as a function
                # agrees with `f_R`
                # @test f_R(ρ) ≈ evaluate(J(ρ)) atol=1e-4

                # check against the naive definition of the Choi matrix
                J_R = sum( f_R(ketbra(i,j, dA)) ⊗ ketbra(i,j, dA) for i = 1:dA, j = 1:dA )
                @test J_R ≈ evaluate(J) atol=1e-4

            end

        end
    end

end


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type UninterpretedOption_NamePart
    name_part::AbstractString
    is_extension::Bool
    UninterpretedOption_NamePart(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type UninterpretedOption_NamePart
const __req_UninterpretedOption_NamePart = Symbol[:name_part,:is_extension]
meta(t::Type{UninterpretedOption_NamePart}) = meta(t, __req_UninterpretedOption_NamePart, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::UninterpretedOption_NamePart) = ProtoBuf.protohash(v)
isequal(v1::UninterpretedOption_NamePart, v2::UninterpretedOption_NamePart) = ProtoBuf.protoisequal(v1, v2)
==(v1::UninterpretedOption_NamePart, v2::UninterpretedOption_NamePart) = ProtoBuf.protoeq(v1, v2)

type UninterpretedOption
    name::Array{UninterpretedOption_NamePart,1}
    identifier_value::AbstractString
    positive_int_value::UInt64
    negative_int_value::Int64
    double_value::Float64
    string_value::Array{UInt8,1}
    aggregate_value::AbstractString
    UninterpretedOption(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type UninterpretedOption
const __fnum_UninterpretedOption = Int[2,3,4,5,6,7,8]
meta(t::Type{UninterpretedOption}) = meta(t, ProtoBuf.DEF_REQ, __fnum_UninterpretedOption, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::UninterpretedOption) = ProtoBuf.protohash(v)
isequal(v1::UninterpretedOption, v2::UninterpretedOption) = ProtoBuf.protoisequal(v1, v2)
==(v1::UninterpretedOption, v2::UninterpretedOption) = ProtoBuf.protoeq(v1, v2)

type __enum_FieldOptions_CType <: ProtoEnum
    STRING::Int32
    CORD::Int32
    STRING_PIECE::Int32
    __enum_FieldOptions_CType() = new(0,1,2)
end #type __enum_FieldOptions_CType
const FieldOptions_CType = __enum_FieldOptions_CType()

type __enum_FieldOptions_JSType <: ProtoEnum
    JS_NORMAL::Int32
    JS_STRING::Int32
    JS_NUMBER::Int32
    __enum_FieldOptions_JSType() = new(0,1,2)
end #type __enum_FieldOptions_JSType
const FieldOptions_JSType = __enum_FieldOptions_JSType()

type FieldOptions
    ctype::Int32
    packed::Bool
    jstype::Int32
    lazy::Bool
    deprecated::Bool
    weak::Bool
    uninterpreted_option::Array{UninterpretedOption,1}
    FieldOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type FieldOptions
const __val_FieldOptions = @compat Dict(:ctype => FieldOptions_CType.STRING, :jstype => FieldOptions_JSType.JS_NORMAL, :lazy => false, :deprecated => false, :weak => false)
const __fnum_FieldOptions = Int[1,2,6,5,3,10,999]
meta(t::Type{FieldOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FieldOptions, __val_FieldOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::FieldOptions) = ProtoBuf.protohash(v)
isequal(v1::FieldOptions, v2::FieldOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::FieldOptions, v2::FieldOptions) = ProtoBuf.protoeq(v1, v2)

type MessageOptions
    message_set_wire_format::Bool
    no_standard_descriptor_accessor::Bool
    deprecated::Bool
    map_entry::Bool
    uninterpreted_option::Array{UninterpretedOption,1}
    MessageOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type MessageOptions
const __val_MessageOptions = @compat Dict(:message_set_wire_format => false, :no_standard_descriptor_accessor => false, :deprecated => false)
const __fnum_MessageOptions = Int[1,2,3,7,999]
meta(t::Type{MessageOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_MessageOptions, __val_MessageOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::MessageOptions) = ProtoBuf.protohash(v)
isequal(v1::MessageOptions, v2::MessageOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::MessageOptions, v2::MessageOptions) = ProtoBuf.protoeq(v1, v2)

type EnumOptions
    allow_alias::Bool
    deprecated::Bool
    uninterpreted_option::Array{UninterpretedOption,1}
    EnumOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type EnumOptions
const __val_EnumOptions = @compat Dict(:deprecated => false)
const __fnum_EnumOptions = Int[2,3,999]
meta(t::Type{EnumOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_EnumOptions, __val_EnumOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::EnumOptions) = ProtoBuf.protohash(v)
isequal(v1::EnumOptions, v2::EnumOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::EnumOptions, v2::EnumOptions) = ProtoBuf.protoeq(v1, v2)

type MethodOptions
    deprecated::Bool
    uninterpreted_option::Array{UninterpretedOption,1}
    MethodOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type MethodOptions
const __val_MethodOptions = @compat Dict(:deprecated => false)
const __fnum_MethodOptions = Int[33,999]
meta(t::Type{MethodOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_MethodOptions, __val_MethodOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::MethodOptions) = ProtoBuf.protohash(v)
isequal(v1::MethodOptions, v2::MethodOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::MethodOptions, v2::MethodOptions) = ProtoBuf.protoeq(v1, v2)

type __enum_FileOptions_OptimizeMode <: ProtoEnum
    SPEED::Int32
    CODE_SIZE::Int32
    LITE_RUNTIME::Int32
    __enum_FileOptions_OptimizeMode() = new(1,2,3)
end #type __enum_FileOptions_OptimizeMode
const FileOptions_OptimizeMode = __enum_FileOptions_OptimizeMode()

type FileOptions
    java_package::AbstractString
    java_outer_classname::AbstractString
    java_multiple_files::Bool
    java_generate_equals_and_hash::Bool
    java_string_check_utf8::Bool
    optimize_for::Int32
    go_package::AbstractString
    cc_generic_services::Bool
    java_generic_services::Bool
    py_generic_services::Bool
    deprecated::Bool
    cc_enable_arenas::Bool
    objc_class_prefix::AbstractString
    csharp_namespace::AbstractString
    uninterpreted_option::Array{UninterpretedOption,1}
    FileOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type FileOptions
const __val_FileOptions = @compat Dict(:java_multiple_files => false, :java_generate_equals_and_hash => false, :java_string_check_utf8 => false, :optimize_for => FileOptions_OptimizeMode.SPEED, :cc_generic_services => false, :java_generic_services => false, :py_generic_services => false, :deprecated => false, :cc_enable_arenas => false)
const __fnum_FileOptions = Int[1,8,10,20,27,9,11,16,17,18,23,31,36,37,999]
meta(t::Type{FileOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FileOptions, __val_FileOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::FileOptions) = ProtoBuf.protohash(v)
isequal(v1::FileOptions, v2::FileOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::FileOptions, v2::FileOptions) = ProtoBuf.protoeq(v1, v2)

type EnumValueOptions
    deprecated::Bool
    uninterpreted_option::Array{UninterpretedOption,1}
    EnumValueOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type EnumValueOptions
const __val_EnumValueOptions = @compat Dict(:deprecated => false)
const __fnum_EnumValueOptions = Int[1,999]
meta(t::Type{EnumValueOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_EnumValueOptions, __val_EnumValueOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::EnumValueOptions) = ProtoBuf.protohash(v)
isequal(v1::EnumValueOptions, v2::EnumValueOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::EnumValueOptions, v2::EnumValueOptions) = ProtoBuf.protoeq(v1, v2)

type OneofOptions
    uninterpreted_option::Array{UninterpretedOption,1}
    OneofOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type OneofOptions
const __fnum_OneofOptions = Int[999]
meta(t::Type{OneofOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_OneofOptions, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::OneofOptions) = ProtoBuf.protohash(v)
isequal(v1::OneofOptions, v2::OneofOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::OneofOptions, v2::OneofOptions) = ProtoBuf.protoeq(v1, v2)

type ServiceOptions
    deprecated::Bool
    uninterpreted_option::Array{UninterpretedOption,1}
    ServiceOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type ServiceOptions
const __val_ServiceOptions = @compat Dict(:deprecated => false)
const __fnum_ServiceOptions = Int[33,999]
meta(t::Type{ServiceOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ServiceOptions, __val_ServiceOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::ServiceOptions) = ProtoBuf.protohash(v)
isequal(v1::ServiceOptions, v2::ServiceOptions) = ProtoBuf.protoisequal(v1, v2)
==(v1::ServiceOptions, v2::ServiceOptions) = ProtoBuf.protoeq(v1, v2)

type __enum_FieldDescriptorProto_Type <: ProtoEnum
    TYPE_DOUBLE::Int32
    TYPE_FLOAT::Int32
    TYPE_INT64::Int32
    TYPE_UINT64::Int32
    TYPE_INT32::Int32
    TYPE_FIXED64::Int32
    TYPE_FIXED32::Int32
    TYPE_BOOL::Int32
    TYPE_STRING::Int32
    TYPE_GROUP::Int32
    TYPE_MESSAGE::Int32
    TYPE_BYTES::Int32
    TYPE_UINT32::Int32
    TYPE_ENUM::Int32
    TYPE_SFIXED32::Int32
    TYPE_SFIXED64::Int32
    TYPE_SINT32::Int32
    TYPE_SINT64::Int32
    __enum_FieldDescriptorProto_Type() = new(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)
end #type __enum_FieldDescriptorProto_Type
const FieldDescriptorProto_Type = __enum_FieldDescriptorProto_Type()

type __enum_FieldDescriptorProto_Label <: ProtoEnum
    LABEL_OPTIONAL::Int32
    LABEL_REQUIRED::Int32
    LABEL_REPEATED::Int32
    __enum_FieldDescriptorProto_Label() = new(1,2,3)
end #type __enum_FieldDescriptorProto_Label
const FieldDescriptorProto_Label = __enum_FieldDescriptorProto_Label()

type FieldDescriptorProto
    name::AbstractString
    number::Int32
    label::Int32
    _type::Int32
    type_name::AbstractString
    extendee::AbstractString
    default_value::AbstractString
    oneof_index::Int32
    json_name::AbstractString
    options::FieldOptions
    FieldDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type FieldDescriptorProto
const __fnum_FieldDescriptorProto = Int[1,3,4,5,6,2,7,9,10,8]
meta(t::Type{FieldDescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FieldDescriptorProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::FieldDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::FieldDescriptorProto, v2::FieldDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::FieldDescriptorProto, v2::FieldDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type MethodDescriptorProto
    name::AbstractString
    input_type::AbstractString
    output_type::AbstractString
    options::MethodOptions
    client_streaming::Bool
    server_streaming::Bool
    MethodDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type MethodDescriptorProto
const __val_MethodDescriptorProto = @compat Dict(:client_streaming => false, :server_streaming => false)
meta(t::Type{MethodDescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, __val_MethodDescriptorProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::MethodDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::MethodDescriptorProto, v2::MethodDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::MethodDescriptorProto, v2::MethodDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type EnumValueDescriptorProto
    name::AbstractString
    number::Int32
    options::EnumValueOptions
    EnumValueDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type EnumValueDescriptorProto
hash(v::EnumValueDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::EnumValueDescriptorProto, v2::EnumValueDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::EnumValueDescriptorProto, v2::EnumValueDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type EnumDescriptorProto
    name::AbstractString
    value::Array{EnumValueDescriptorProto,1}
    options::EnumOptions
    EnumDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type EnumDescriptorProto
hash(v::EnumDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::EnumDescriptorProto, v2::EnumDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::EnumDescriptorProto, v2::EnumDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type OneofDescriptorProto
    name::AbstractString
    options::OneofOptions
    OneofDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type OneofDescriptorProto
hash(v::OneofDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::OneofDescriptorProto, v2::OneofDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::OneofDescriptorProto, v2::OneofDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type DescriptorProto_ExtensionRange
    start::Int32
    _end::Int32
    DescriptorProto_ExtensionRange(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type DescriptorProto_ExtensionRange
hash(v::DescriptorProto_ExtensionRange) = ProtoBuf.protohash(v)
isequal(v1::DescriptorProto_ExtensionRange, v2::DescriptorProto_ExtensionRange) = ProtoBuf.protoisequal(v1, v2)
==(v1::DescriptorProto_ExtensionRange, v2::DescriptorProto_ExtensionRange) = ProtoBuf.protoeq(v1, v2)

type DescriptorProto_ReservedRange
    start::Int32
    _end::Int32
    DescriptorProto_ReservedRange(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type DescriptorProto_ReservedRange
hash(v::DescriptorProto_ReservedRange) = ProtoBuf.protohash(v)
isequal(v1::DescriptorProto_ReservedRange, v2::DescriptorProto_ReservedRange) = ProtoBuf.protoisequal(v1, v2)
==(v1::DescriptorProto_ReservedRange, v2::DescriptorProto_ReservedRange) = ProtoBuf.protoeq(v1, v2)

type DescriptorProto
    name::AbstractString
    field::Array{FieldDescriptorProto,1}
    extension::Array{FieldDescriptorProto,1}
    nested_type::Array{DescriptorProto,1}
    enum_type::Array{EnumDescriptorProto,1}
    extension_range::Array{DescriptorProto_ExtensionRange,1}
    oneof_decl::Array{OneofDescriptorProto,1}
    options::MessageOptions
    reserved_range::Array{DescriptorProto_ReservedRange,1}
    reserved_name::Array{AbstractString,1}
    DescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type DescriptorProto
const __fnum_DescriptorProto = Int[1,2,6,3,4,5,8,7,9,10]
meta(t::Type{DescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_DescriptorProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::DescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::DescriptorProto, v2::DescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::DescriptorProto, v2::DescriptorProto) = ProtoBuf.protoeq(v1, v2)

type ServiceDescriptorProto
    name::AbstractString
    method::Array{MethodDescriptorProto,1}
    options::ServiceOptions
    ServiceDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type ServiceDescriptorProto
hash(v::ServiceDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::ServiceDescriptorProto, v2::ServiceDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ServiceDescriptorProto, v2::ServiceDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type SourceCodeInfo_Location
    path::Array{Int32,1}
    span::Array{Int32,1}
    leading_comments::AbstractString
    trailing_comments::AbstractString
    leading_detached_comments::Array{AbstractString,1}
    SourceCodeInfo_Location(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type SourceCodeInfo_Location
const __fnum_SourceCodeInfo_Location = Int[1,2,3,4,6]
const __pack_SourceCodeInfo_Location = Symbol[:path,:span]
meta(t::Type{SourceCodeInfo_Location}) = meta(t, ProtoBuf.DEF_REQ, __fnum_SourceCodeInfo_Location, ProtoBuf.DEF_VAL, true, __pack_SourceCodeInfo_Location, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::SourceCodeInfo_Location) = ProtoBuf.protohash(v)
isequal(v1::SourceCodeInfo_Location, v2::SourceCodeInfo_Location) = ProtoBuf.protoisequal(v1, v2)
==(v1::SourceCodeInfo_Location, v2::SourceCodeInfo_Location) = ProtoBuf.protoeq(v1, v2)

type SourceCodeInfo
    location::Array{SourceCodeInfo_Location,1}
    SourceCodeInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type SourceCodeInfo
hash(v::SourceCodeInfo) = ProtoBuf.protohash(v)
isequal(v1::SourceCodeInfo, v2::SourceCodeInfo) = ProtoBuf.protoisequal(v1, v2)
==(v1::SourceCodeInfo, v2::SourceCodeInfo) = ProtoBuf.protoeq(v1, v2)

type FileDescriptorProto
    name::AbstractString
    package::AbstractString
    dependency::Array{AbstractString,1}
    public_dependency::Array{Int32,1}
    weak_dependency::Array{Int32,1}
    message_type::Array{DescriptorProto,1}
    enum_type::Array{EnumDescriptorProto,1}
    service::Array{ServiceDescriptorProto,1}
    extension::Array{FieldDescriptorProto,1}
    options::FileOptions
    source_code_info::SourceCodeInfo
    syntax::AbstractString
    FileDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type FileDescriptorProto
const __fnum_FileDescriptorProto = Int[1,2,3,10,11,4,5,6,7,8,9,12]
meta(t::Type{FileDescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FileDescriptorProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::FileDescriptorProto) = ProtoBuf.protohash(v)
isequal(v1::FileDescriptorProto, v2::FileDescriptorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::FileDescriptorProto, v2::FileDescriptorProto) = ProtoBuf.protoeq(v1, v2)

type FileDescriptorSet
    file::Array{FileDescriptorProto,1}
    FileDescriptorSet(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type FileDescriptorSet
hash(v::FileDescriptorSet) = ProtoBuf.protohash(v)
isequal(v1::FileDescriptorSet, v2::FileDescriptorSet) = ProtoBuf.protoisequal(v1, v2)
==(v1::FileDescriptorSet, v2::FileDescriptorSet) = ProtoBuf.protoeq(v1, v2)

type GeneratedCodeInfo_Annotation
    path::Array{Int32,1}
    source_file::AbstractString
    _begin::Int32
    _end::Int32
    GeneratedCodeInfo_Annotation(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type GeneratedCodeInfo_Annotation
const __pack_GeneratedCodeInfo_Annotation = Symbol[:path]
meta(t::Type{GeneratedCodeInfo_Annotation}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_GeneratedCodeInfo_Annotation, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
hash(v::GeneratedCodeInfo_Annotation) = ProtoBuf.protohash(v)
isequal(v1::GeneratedCodeInfo_Annotation, v2::GeneratedCodeInfo_Annotation) = ProtoBuf.protoisequal(v1, v2)
==(v1::GeneratedCodeInfo_Annotation, v2::GeneratedCodeInfo_Annotation) = ProtoBuf.protoeq(v1, v2)

type GeneratedCodeInfo
    annotation::Array{GeneratedCodeInfo_Annotation,1}
    GeneratedCodeInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type GeneratedCodeInfo
hash(v::GeneratedCodeInfo) = ProtoBuf.protohash(v)
isequal(v1::GeneratedCodeInfo, v2::GeneratedCodeInfo) = ProtoBuf.protoisequal(v1, v2)
==(v1::GeneratedCodeInfo, v2::GeneratedCodeInfo) = ProtoBuf.protoeq(v1, v2)

export FileDescriptorSet, FileDescriptorProto, DescriptorProto_ExtensionRange, DescriptorProto_ReservedRange, DescriptorProto, FieldDescriptorProto_Type, FieldDescriptorProto_Label, FieldDescriptorProto, OneofDescriptorProto, EnumDescriptorProto, EnumValueDescriptorProto, ServiceDescriptorProto, MethodDescriptorProto, FileOptions_OptimizeMode, FileOptions, MessageOptions, FieldOptions_CType, FieldOptions_JSType, FieldOptions, OneofOptions, EnumOptions, EnumValueOptions, ServiceOptions, MethodOptions, UninterpretedOption_NamePart, UninterpretedOption, SourceCodeInfo_Location, SourceCodeInfo, GeneratedCodeInfo_Annotation, GeneratedCodeInfo


text:
const QUAD_SOLVER = optimizer_with_attributes(OSQP.Optimizer, "verbose" => false, "max_iter" => 15000)
const LIN_SOLVER = optimizer_with_attributes(GLPK.Optimizer)
const TWO_PI_OVER_365 = 2 * π / 365

function calc_J(interest_yield::Number, maturity_interval::UnitRange{Int64})
    return sum((1 + interest_yield)^(-j) for j in maturity_interval)
end

function get_index_and_length_of_traded_contracts(contracts_by_day)
    pos_of_traded_contracts = Int[]
    for i = 1:length(contracts_by_day.contracts)
        if !isnan(contracts_by_day.contracts[i].price)
            push!(pos_of_traded_contracts, i)
        end
    end
    return pos_of_traded_contracts, length(pos_of_traded_contracts)
end

function construct_objective_function_from_adjustment(ω)
    objective = 0.0
    for t in 1:length(ω), adjust_var in ω[t]
        objective += adjust_var^2
    end
    return objective
end

function structural_expression(args_structural_expression, X)
    return sum(args_structural_expression[:beta] .* X)
end

function save_estimated_parameters!(args_structural_expression::Dict,
                                    coeficients::Vector{Float64})
    args_structural_expression[:beta] = coeficients
    return
end

function create_F_stack(input::StructuralModelInputs)
    vec_contracts_by_day = input.vec_contracts_by_day
    T                    = length(date_of_trades(vec_contracts_by_day))
    I                    = number_of_contracts(vec_contracts_by_day)
    F_matrix             = Matrix(undef, T, I)
    contracts_idx_matrix = Int.(hcat([ones(T) * i for i = 1:I]...))

    for i in 1:I, t in 1:T
        F_matrix[t, i] = vec_contracts_by_day[t].contracts[i].price
    end

    F_stack             = vcat(F_matrix...)
    contracts_idx_stack = vcat(contracts_idx_matrix...)
    remove_missing_idx  = isnan.(F_stack) .== false
    F_stack             = F_stack[remove_missing_idx]
    contracts_idx_stack = contracts_idx_stack[remove_missing_idx]

    return F_stack, remove_missing_idx, contracts_idx_stack
end

function from_X_t_j_to_X_t_i(input::StructuralModelInputs)
    interest_yield       = input.parameters["interest_yield"]
    vec_contracts_by_day = input.vec_contracts_by_day
    X                    = input.X
    I                    = number_of_contracts(vec_contracts_by_day)
    (T, J, L)            = size(X)
    X_t_i_matrix         = Array{Any,3}(undef, T, I, L)
    for (t, contract_by_day) in enumerate(vec_contracts_by_day)
        for (i, contract_info) in enumerate(contract_by_day.contracts)
            duration         = contract_info.duration
            time_to_maturity = contract_info.time_to_maturity
            maturity_range   = time_to_maturity:time_to_maturity + duration - 1
            for l = 1:L
                X_t_i_matrix[t, i, l] = sum(
                    X[t, j, l] * (1 + interest_yield)^(-j) for j in maturity_range
                ) / calc_J(interest_yield, maturity_range)
            end
        end
    end
    return X_t_i_matrix
end

function stack_X_matrix(X_t_i_matrix::Array)
    (T, I, L) = size(X_t_i_matrix)
    X_t_i_stack = Matrix{Float64}(undef, T * I, L)
    for i in 1:I, t in 1:T, l in 1:L
        X_t_i_stack[T * (i - 1) + t, l] = X_t_i_matrix[t, i, l]
    end
    return X_t_i_stack
end

function create_X_stack(input::StructuralModelInputs, remove_missing_idx::BitArray)
    X_t_i_matrix = from_X_t_j_to_X_t_i(input)
    X_t_i_stack  = stack_X_matrix(X_t_i_matrix)

    return  X_t_i_stack[remove_missing_idx, :]
end

function estimate_coeficients(F_stack::Vector, X_stack::Matrix)
    glm_model             = lm(X_stack, F_stack)
    estimated_coeficients = coef(glm_model)
    return estimated_coeficients, glm_model
end

"""
    Verify if X matrix has a column with all entries 'ones'.
        If true: return the matrix without this column and the idx
        If false: return the X matrix and a black vector
"""
function check_for_intercept(X::Matrix)
    has_intercept = false
    num_intercepts = 0
    for i = 1:size(X, 2)
        if all(X[:, i] .== 1)
            num_intercepts += 1
            has_intercept = true
        end
    end
    if num_intercepts > 1
        @error("Number of intercepts in structural model is greater than one. Revise the structural variables")
    end
    if has_intercept
        lin, col = size(X)
        X_wout_intercept = Matrix{Float64}(undef, lin, col - 1)
        (j, removed_idx) = (1, 0)
        for i = 1:size(X, 2)
            if !all(X[:, i] .== 1)
                X_wout_intercept[:,j] = X[:, i]
                j += 1
            else
                removed_idx = i
            end
        end
        return X_wout_intercept, removed_idx
    end
    return X, []
end

function structural_lasso_selection(input::StructuralModelInputs,
                                    X_stack::Matrix{Float64},
                                    F_stack::Vector{Float64})
    if input.parameters["run_lasso"]
        println("Selecting variables through AdaLasso.")
        return structural_lasso_selection(X_stack, F_stack)
    else
        println("No variable selection.")
        all_idx = collect(1:size(X_stack, 2))
        return "", all_idx
    end
end

function structural_lasso_selection(X_stack::Matrix{Float64}, y::Vector{Float64})
    X, rmv_idx = check_for_intercept(X_stack)
    T, L = size(X)
    cp_idx = collect(1:L)
    LASSO = try
        Lasso.fit(LassoModel, X, y; select=MinBIC())
    catch
        "Failed to converge Lasso algorithm in structural model. Running structural model only with intercept"
    end
    if typeof(LASSO) == String
        @warn(LASSO)
        return "", setdiff(cp_idx, rmv_idx)
    end
    LASSO_coef_without_intercept = coef(LASSO)[2:end]
    cp_idx_nonzeros = cp_idx[abs.(LASSO_coef_without_intercept) .> 0]
    pen = abs.(LASSO_coef_without_intercept[cp_idx_nonzeros]).^(-1)
    X_ada = X[:, cp_idx_nonzeros]

    # AdaLasso
    AdaLASSO = try
        Lasso.fit(LassoModel, X_ada, y, penalty_factor=pen, select=MinBIC())
    catch
        "Failed to converge AdaLasso algorithm in structural model. Running structural model only with intercept"
    end
    if typeof(AdaLASSO) == String
        @warn(AdaLASSO)
        return "", setdiff(cp_idx, rmv_idx)
    end
    AdaLASSO_coef_without_intercept = coef(AdaLASSO)[2:end]
    ada_idx_nonzeros = AdaLASSO_coef_without_intercept .!= 0
    AdaLASSO_coef_without_intercept = AdaLASSO_coef_without_intercept[ada_idx_nonzeros]
    variables_idx = cp_idx_nonzeros[ada_idx_nonzeros]
    variables_idx_final = union(setdiff(collect(1:size(X_stack, 2)), rmv_idx)[variables_idx], rmv_idx)
    return AdaLASSO, sort(variables_idx_final)
end

"""
"""
function first_level_estimate_parameters(input::StructuralModelInputs)
    println(stdout, "Estimating seasonality and explanatory coefficients.")
    args_structural_expression = Dict{Symbol,Any}()
    # Create F_stack
    F_stack, remove_missing_idx, contracts_idx = create_F_stack(input)
    # Create X_stack
    X_stack = create_X_stack(input, remove_missing_idx)
    struct_AdaLASSO, select_variables_idx = structural_lasso_selection(input, X_stack, F_stack)
    X_stack = X_stack[:, select_variables_idx]
    estimated_coeficients, glm_model = estimate_coeficients(F_stack, X_stack)

    save_estimated_parameters!(args_structural_expression, estimated_coeficients)
    return (struct_AdaLASSO, select_variables_idx, args_structural_expression, glm_model)
end

function maximum_maturity_from_vec_contracts_by_day(vec_contracts_by_day::Vector{ContractsByDay})
    maximum_maturity = 1
    for contracts_by_day in vec_contracts_by_day, contract_info in contracts_by_day.contracts
        if contract_info.time_to_maturity + contract_info.duration - 1 > maximum_maturity
            maximum_maturity = contract_info.time_to_maturity + contract_info.duration - 1
        end
    end
    return maximum_maturity
end

function maximum_maturity_from_contracts_by_day(contracts_by_day::ContractsByDay)
    maximum_maturity = 1
    for contract_info in contracts_by_day.contracts
        if contract_info.time_to_maturity + contract_info.duration - 1 > maximum_maturity
            maximum_maturity = contract_info.time_to_maturity + contract_info.duration - 1
        end
    end
    maximum_maturity
end

function create_second_level_sub_problem_variables!(model::JuMP.Model,
                                                    maximum_maturity::Int,
                                                    N::Int)
    @variable(model, ε[1:maximum_maturity])
    @variable(model, ζ[1:N])
    @variable(model, ζ_abs[1:N])
    return
end

function create_second_level_sub_problem_constraints!(model::JuMP.Model,
                                                      args_structural_expression::Dict{Symbol,Any},
                                                      X_t::Matrix,
                                                      contracts_by_day::ContractsByDay,
                                                      N::Int,
                                                      interest_yield::Number,
                                                      index_of_traded_contracts::Vector{Int})
    for i = 1:N
        contract_info              = contracts_by_day.contracts[index_of_traded_contracts[i]]
        price                      = contract_info.price
        contract_maturity_interval = contract_info.time_to_maturity:contract_info.time_to_maturity + contract_info.duration - 1

        @constraint(model, contract_info.price == (
            sum(
                (structural_expression(args_structural_expression, X_t[j,:]) + model[:ε][j]) *
                (1 + interest_yield)^(-j) for j in contract_maturity_interval
            ) /
            calc_J(interest_yield, contract_maturity_interval)
            ) +
                model[:ζ][i]
            )
        @constraint(model, model[:ζ_abs][i] >=  model[:ζ][i])
        @constraint(model, model[:ζ_abs][i] >= -model[:ζ][i])
    end
    return
end

function create_second_level_sub_problem_objective!(model)
    return @objective(model, Min, sum(model[:ζ_abs]))
end

function sub_problem_daily_minimum_arbitrage(args_structural_expression::Dict{Symbol,Any},
                                             X_t::Matrix,
                                             contracts_by_day::ContractsByDay,
                                             interest_yield::Number)

    maximum_maturity = maximum_maturity_from_contracts_by_day(contracts_by_day)
    index_of_traded_contracts, N = get_index_and_length_of_traded_contracts(contracts_by_day)

    if N != 0
        model = Model(LIN_SOLVER)
        create_second_level_sub_problem_variables!(model, maximum_maturity, N)
        create_second_level_sub_problem_constraints!(model, args_structural_expression,
                                                     X_t, contracts_by_day, N,
                                                     interest_yield,
                                                     index_of_traded_contracts)
        create_second_level_sub_problem_objective!(model)

        optimize!(model)

        if (termination_status(model) != MOI.OPTIMAL)
            error("O problema não foi resolvido até a otimalidade. $(termination_status(model))")
        end

        return objective_value(model)
    else
        return 0.0
    end
end

function second_level_daily_minimum_arbitrage(args_structural_expression::Dict{Symbol,Any},
                                              input::StructuralModelInputs)
    println(stdout, "Minimum arbitrage calculation.")
    vector_of_daily_minimum_arbitrage = Vector{Float64}(undef, length(input.vec_contracts_by_day))
    X = input.X
    for (t, contracts_by_day) in enumerate(input.vec_contracts_by_day)
        daily_minimum_arbitrage = sub_problem_daily_minimum_arbitrage(
                                                args_structural_expression,
                                                X[t,:,:],
                                                contracts_by_day,
                                                input.parameters["interest_yield"]
                                            )
        vector_of_daily_minimum_arbitrage[t] = daily_minimum_arbitrage
    end

    return vector_of_daily_minimum_arbitrage
end

function handle_f0_constraint!(model::JuMP.Model,
                               input::StructuralModelInputs,
                               args_structural_expression::Dict{Symbol,Any},
                               X_t_j::Vector,
                               t::Int)
    if haskey(input.parameters, "f0_spot")
        @constraint(
                model,
                structural_expression(args_structural_expression, X_t_j)
                +
                model[:ε][t,0]
                == input.parameters["f0_spot"].vals[t]
                )
    end
    return
end

maturity_0(X_t::Matrix) = X_t[end, :]

function third_level_common_constraints!(model::JuMP.Model,
                                         input::StructuralModelInputs,
                                         args_structural_expression::Dict{Symbol,Any})
    T              = length(input.vec_contracts_by_day)
    interest_yield = input.parameters["interest_yield"]
    ζ              = Vector{Vector{VariableRef}}(undef, T)
    ζ_abs          = Vector{Vector{VariableRef}}(undef, T)
    X              = input.X
    arbitrage_constraint = 0.0
    for (t, contracts_by_day) in enumerate(input.vec_contracts_by_day)
        index_of_traded_contracts, N = get_index_and_length_of_traded_contracts(contracts_by_day)
        ζ[t]     = @variable(model, [1:N])
        ζ_abs[t] = @variable(model, [1:N])
        handle_f0_constraint!(
            model,
            input,
            args_structural_expression,
            maturity_0(X[t,:,:]),
            t)

        for i = 1:N
            contract_info              = contracts_by_day.contracts[index_of_traded_contracts[i]]
            contract_maturity_interval = contract_info.time_to_maturity:contract_info.time_to_maturity +
                                         contract_info.duration - 1

            @constraint(model, contract_info.price == (
                sum(
                    (structural_expression(args_structural_expression, X[t, j, :]) + model[:ε][t,j]) *
                    (1 + interest_yield)^(-j) for j in contract_maturity_interval
                ) /
                calc_J(interest_yield, contract_maturity_interval)
                ) +
                    ζ[t][i]
                )
            @constraint(model, ζ_abs[t][i] >=  ζ[t][i])
            @constraint(model, ζ_abs[t][i] >= -ζ[t][i])
            arbitrage_constraint += ζ_abs[t][i]
        end
    end
    return arbitrage_constraint
end

function sub_problem_smooth_through_maturiy(args_structural_expression::Dict{Symbol,Any},
                                            input::StructuralModelInputs,
                                            vector_of_daily_minimum_arbitrage::Vector{Float64})
    println(stdout, "Getting the smoothing in maturities weight.")

    T = length(input.vec_contracts_by_day)
    maximum_maturity_contracts = maximum_maturity_from_vec_contracts_by_day(input.vec_contracts_by_day)
    maximum_maturity = max(maximum_maturity_contracts, input.parameters["minimum_maturity"])

    model = Model(QUAD_SOLVER)
    # structural error
    @variable(model, ε[0:T, 0:maximum_maturity])
    # smooth through maturity
    @variable(model, Γ[1:T, 1:maximum_maturity])
    arbitrage_constraint = third_level_common_constraints!(model, input, args_structural_expression)
    @constraint(model, [t = 1:T, j = 1:maximum_maturity - 1], Γ[t, j] == ε[t, j + 1] - 2 * ε[t, j] + ε[t, j - 1])
    @constraint(model, arbitrage_constraint <= sum(vector_of_daily_minimum_arbitrage))
    @constraint(model, [t = 1:T], ε[t, end] - ε[t, end - 1] == 0.0)
    @objective(model, Min, sum(Γ.^2))
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    return objective_value(model)
end

function sub_problem_smooth_through_time(args_structural_expression::Dict{Symbol,Any},
                                         input::StructuralModelInputs,
                                         vector_of_daily_minimum_arbitrage::Vector{Float64})
    println(stdout, "Getting the smoothing in time weight.")
    T = length(input.vec_contracts_by_day)
    maximum_maturity_contracts = maximum_maturity_from_vec_contracts_by_day(input.vec_contracts_by_day)
    maximum_maturity = max(maximum_maturity_contracts, input.parameters["minimum_maturity"])

    model = Model(QUAD_SOLVER)
    # structural error
    @variable(model, ε[0:T, 0:maximum_maturity])
    # smooth through time
    @variable(model, Λ[2:T, 1:maximum_maturity])
    arbitrage_constraint = third_level_common_constraints!(model, input, args_structural_expression)
    @constraint(model, [t = 2:T - 1, j = 1:maximum_maturity], Λ[t,j] == ε[t + 1, j] - 2 * ε[t, j] + ε[t - 1, j])
    @constraint(model, arbitrage_constraint <= sum(vector_of_daily_minimum_arbitrage))
    @constraint(model, [t = 1:T], ε[t, end] - ε[t, end - 1] == 0.0)
    @objective(model, Min, sum(Λ.^2))
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    return objective_value(model)
end

function handle_structural_error_value(value_ε)
    ## since we created a f0, the structural error must be treated
    lin, col = size(value_ε)
    structural_error = Array{Float64}(undef, lin - 1, col - 1)
    for l in 1:lin - 1, c in 1:col - 1
        structural_error[l, c] = value_ε[l, c]
    end
    return structural_error
end

function estimate_structural_error(args_structural_expression::Dict{Symbol,Any},
                                   input::StructuralModelInputs,
                                   vector_of_daily_minimum_arbitrage::Vector{Float64},
                                   time_weigth::Float64,
                                   maturity_weigth::Float64)
    println(stdout, "Estimating the structural model error.")
    T = length(input.vec_contracts_by_day)
    maximum_maturity_contracts = maximum_maturity_from_vec_contracts_by_day(input.vec_contracts_by_day)
    maximum_maturity = max(maximum_maturity_contracts, input.parameters["minimum_maturity"])

    model = Model(QUAD_SOLVER)
    # structural error
    @variable(model, ε[0:T, 0:maximum_maturity])
    # smooth through maturity
    @variable(model, Γ[1:T, 1:maximum_maturity])
    # smooth through time
    @variable(model, Λ[2:T, 1:maximum_maturity])
    arbitrage_constraint = third_level_common_constraints!(model, input, args_structural_expression)
    @constraint(model, [t = 1:T, j = 1:maximum_maturity - 1], Γ[t,j] == ε[t, j + 1] - 2 * ε[t, j] + ε[t, j - 1])
    @constraint(model, [t = 2:T - 1, j = 1:maximum_maturity], Λ[t,j] == ε[t + 1, j] - 2 * ε[t, j] + ε[t - 1, j])
    @constraint(model, [t = 1:T], ε[t, end] - ε[t, end - 1] == 0.0)
    @constraint(model, arbitrage_constraint <= sum(vector_of_daily_minimum_arbitrage))
    @objective(model, Min, (input.parameters["λ"] / time_weigth) * sum(Λ.^2) + ((1 - input.parameters["λ"]) / maturity_weigth) * sum(Γ.^2))
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    ε_val = handle_structural_error_value(JuMP.value.(ε))
    return objective_value(model), ε_val
end

function third_level_estimate_structural_error(args_structural_expression,
                                               input,
                                               vector_of_daily_minimum_arbitrage::Vector{Float64})
    maturity_smooth_time = @elapsed begin
        maturity_weigth = sub_problem_smooth_through_maturiy(args_structural_expression,
                                                         input,
                                                         vector_of_daily_minimum_arbitrage)
    end
    println(stdout, "Finished in $(round(maturity_smooth_time, digits=2)) s")

    time_smooth_time = @elapsed begin
        time_weigth = sub_problem_smooth_through_time(args_structural_expression,
                                                         input,
                                                         vector_of_daily_minimum_arbitrage)
    end
    println(stdout, "Finished in $(round(time_smooth_time, digits=2)) s")

    estimate_strutural_error_time = @elapsed begin
        (struct_err_obj_function,
            structural_error) = estimate_structural_error(args_structural_expression,
                                                        input,
                                                        vector_of_daily_minimum_arbitrage,
                                                        time_weigth, maturity_weigth)
    end
    println(stdout, "Finished in $(round(estimate_strutural_error_time, digits=2)) s")

    return (maturity_weigth, time_weigth, struct_err_obj_function, structural_error,
            ThirdLevelSolveTime(maturity_smooth_time,
                                    time_smooth_time,
                                    estimate_strutural_error_time))
end

function retrieve_elementary_forward_contracts(args_structural_expression, structural_error, input)
    T, J = size(structural_error)
    X = input.X
    elementary_forward_contracts = Array{Float64}(undef, T, J)
    for t in 1:T, j in 1:J
        elementary_forward_contracts[t,j] = (structural_expression(args_structural_expression,
                                                          X[t, j, :])
                                                            +
                                                          structural_error[t,j]
                                                          )
    end
    return elementary_forward_contracts
end

function reconstruct_forward_contracts(elementary_forward_contracts::Matrix{Float64},
                                       input::StructuralModelInputs)

    reconstructed_vec_contracts_by_day = deepcopy(input.vec_contracts_by_day)
    num_of_contracts = number_of_contracts(reconstructed_vec_contracts_by_day)
    for (t, contracts_by_day) in enumerate(reconstructed_vec_contracts_by_day)
        for i = 1:num_of_contracts
            contract_info = contracts_by_day.contracts[i]
            maturity_interval = contract_info.time_to_maturity:contract_info.time_to_maturity + contract_info.duration - 1
            reconstructed_price = sum(elementary_forward_contracts[t, j]
                            * (1 + input.parameters["interest_yield"])^(-j)
                            for j in maturity_interval) / calc_J(input.parameters["interest_yield"], maturity_interval)
            reconstructed_vec_contracts_by_day[t].contracts[i].price = reconstructed_price
        end
    end

    return reconstructed_vec_contracts_by_day
end

function create_ground_level_sub_problem_variables!(model::JuMP.Model,
                                                    maximum_maturity::Int,
                                                    N::Int)
    @variable(model, ε[1:maximum_maturity])
    @variable(model, ζ[1:N])
    @variable(model, s[1:N] == 0)
    @variable(model, Δ[1:N])
    @variable(model, θ)
    return
end

function create_ground_level_sub_problem_constraints!(model::JuMP.Model,
                                                      contracts_by_day::ContractsByDay,
                                                      N::Int,
                                                      interest_yield::Number,
                                                      index_of_traded_contracts::Vector{Int})
    for i = 1:N
        contract_info              = contracts_by_day.contracts[index_of_traded_contracts[i]]
        price                      = contract_info.price
        contract_maturity_interval = contract_info.time_to_maturity:contract_info.time_to_maturity + contract_info.duration - 1

        @constraint(model, contract_info.price == (
            sum(
                (model[:ε][j]) *
                (1 + interest_yield)^(-j) for j in contract_maturity_interval
            ) /
            calc_J(interest_yield, contract_maturity_interval)
            ) +
                model[:ζ][i]
            )
        @constraint(model, model[:Δ][i] >= model[:ζ][i] + model[:s][i])
        @constraint(model, model[:Δ][i] >= -model[:ζ][i] + model[:s][i])
    end
    @constraint(model, infinity_norm[i=1:N], model[:θ] >= model[:Δ][i])
    return
end

function create_ground_level_sub_problem_objective!(model::JuMP.Model)
    return @objective(model, Min, model[:θ])
end

function optimize_ground_level(model::JuMP.Model, N::Int)
    optimize!(model)
    if (termination_status(model) != MOI.OPTIMAL)
        error("O problema não foi resolvido até a otimalidade. $(termination_status(model))")
    end
    minimum_arbitrage_by_contract = Vector{Float64}(undef, N)
    while termination_status(model) == MOI.OPTIMAL
        _, i_max = findmax(JuMP.dual.(model[:infinity_norm]))
        minimum_arbitrage_by_contract[i_max] = value.(model[:ζ])[i_max]
        fix(model[:ζ][i_max], minimum_arbitrage_by_contract[i_max])
        unfix(model[:s][i_max])
        optimize!(model)
    end
    return minimum_arbitrage_by_contract
end

function sub_problem_contracts_minimum_arbitrage(contracts_by_day,
                                                  interest_yield)
    maximum_maturity = maximum_maturity_from_contracts_by_day(contracts_by_day)
    index_of_traded_contracts, N = get_index_and_length_of_traded_contracts(contracts_by_day)

    if N != 0
        model = Model(LIN_SOLVER)
        create_ground_level_sub_problem_variables!(model, maximum_maturity, N)
        create_ground_level_sub_problem_constraints!(model, contracts_by_day, N,
                                                     interest_yield,
                                                     index_of_traded_contracts)
        create_ground_level_sub_problem_objective!(model)

        minimum_arbitrage_by_contract = optimize_ground_level(model, N)

        return minimum_arbitrage_by_contract
    else
        return []
    end
end

function create_arbitrage_free_input(input::StructuralModelInputs,
                                     vector_of_minimum_arbitrage_by_contract::Vector{Vector{Float64}})

    arbitrage_free_input = deepcopy(input)
    vec_contracts_by_day = input.vec_contracts_by_day

    for (t, contracts_by_day) in enumerate(vec_contracts_by_day)
        j = 1
        for (i, contracts_info) in enumerate(contracts_by_day.contracts)
            if !isnan(contracts_info.price)
                arbitrage_free_price = (
                        contracts_info.price - vector_of_minimum_arbitrage_by_contract[t][j]
                )
                arbitrage_free_input.vec_contracts_by_day[t].contracts[i].price = arbitrage_free_price
                j += 1
            end
        end
    end

    return arbitrage_free_input
end

function ground_level_arbitrage_free_input(input::StructuralModelInputs)
    ### criar problema de otimização para recuperar a arbitragem de cada contratos
    println(stdout, "Creating arbitrage-free forward contracts.")
    vector_of_minimum_arbitrage_by_contract = Vector{Vector{Float64}}(undef, length(input.vec_contracts_by_day))
    for (t, contracts_by_day) in enumerate(input.vec_contracts_by_day)
        minimum_arbitrage_by_contract = sub_problem_contracts_minimum_arbitrage(
                                                contracts_by_day,
                                                input.parameters["interest_yield"]
                                            )
        vector_of_minimum_arbitrage_by_contract[t] = minimum_arbitrage_by_contract
    end

    return (create_arbitrage_free_input(input, vector_of_minimum_arbitrage_by_contract),
        vector_of_minimum_arbitrage_by_contract)
end

function forward_structural_model(input::StructuralModelInputs)
    ground_level_time = @elapsed begin
        (arbitrage_free_input,
            vector_minimum_arbitrage) = ground_level_arbitrage_free_input(input)
    end
    println(stdout, "Finished in $(round(ground_level_time, digits=2)) s")

    first_level_time = @elapsed begin
        (struct_AdaLASSO,
            select_variables_idx,
              args_structural_expression,
                glm_model) = first_level_estimate_parameters(arbitrage_free_input)

        arbitrage_free_input.X = arbitrage_free_input.X[:,:,select_variables_idx] # removing non selected variables by adaLasso
    end
    println(stdout, "Finished in $(round(first_level_time, digits=2)) s")

    second_level_time = @elapsed begin
        vector_of_daily_minimum_arbitrage = second_level_daily_minimum_arbitrage(
                                                args_structural_expression,
                                                arbitrage_free_input
                                            )
    end
    println(stdout, "Finished in $(round(second_level_time, digits=2)) s")

    (maturity_weigth, time_weigth,
        struct_err_obj_function, structural_error,
            third_level_time) = third_level_estimate_structural_error(
                                                args_structural_expression,
                                                arbitrage_free_input,
                                                vector_of_daily_minimum_arbitrage)

    elementary_forward_contracts = retrieve_elementary_forward_contracts(
                                                args_structural_expression,
                                                structural_error, arbitrage_free_input)

    structural_dinamic = elementary_forward_contracts .- structural_error

    reconstructed_forward_prices = reconstruct_forward_contracts(
                                                elementary_forward_contracts,
                                                arbitrage_free_input)

    return StructuralModelOutputs(
                glm_model,
                args_structural_expression,
                vector_of_daily_minimum_arbitrage,
                maturity_weigth,
                time_weigth,
                struct_err_obj_function,
                structural_error,
                structural_dinamic,
                reconstructed_forward_prices,
                StructuralSolveTime(
                        ground_level_time,
                        first_level_time,
                        second_level_time,
                        third_level_time
                    ),
                struct_AdaLASSO,
                select_variables_idx,
                vector_minimum_arbitrage
                )
end

max_stars_repo_name: dietzemarina/ForwardCurveSmoother

max_issues_repo_name: dietzemarina/ForwardCurveSmoother

max_forks_repo_name: dietzemarina/ForwardCurveSmoother


text:
#!/usr/bin/env julia --project=@. --startup-file=no

const sysimage_dir = "sysimage"
isdir(sysimage_dir) || mkdir(sysimage_dir)

import Pkg;
Pkg.add("PackageCompiler")
project_name = Pkg.project().name |> Symbol

import PackageCompiler
PackageCompiler.create_sysimage(
    project_name,
    cpu_target="generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)",
    sysimage_path = joinpath(sysimage_dir, "image.so"),
    # precompile_statements_file=joinpath(sysimage_dir, "precompile_trace.jl"),
    precompile_execution_file = "test/runtests.jl",
)


text:
using RandomizedLinAlg
using Test, LinearAlgebra, Random

@testset "rsvd" begin

    Random.seed!(1234321)

    @testset "Small wide rectangular" begin
        A = [1. 2 3 4; 5 6 7 8]
        S1 = svd(A)
        S2 = rsvd(A, 2, 0)

        @test norm(abs.(S1.U) - abs.(S2.U)) ≤ √(eps())
        @test norm(abs.(S1.Vt) - abs.(S2.Vt)) ≤ √(eps())
        @test norm(S1.S - S2.S) ≤ √(eps())
    end

    @testset "Small tall rectangular" begin
        A = [1. 2; 3 4; 5 6; 7 8]
        S1 = svd(A)
        S2 = rsvd(A, 2, 0)

        @test norm(abs.(S1.U) .- abs.(S2.U)) ≤ √(eps())
        @test norm(abs.(S1.Vt) .- abs.(S2.Vt)) ≤ √(eps())
        @test norm(S1.S - S2.S) ≤ √(eps())
    end

    @testset "Small square" begin
        A = [1. 2 3; 4 5 6; 7 8 9]
        S1 = svd(A)
        S2 = rsvd(A, 3, 0)

        @test norm(abs.(S1.U) - abs.(S2.U)) ≤ √(eps())
        @test norm(abs.(S1.Vt) - abs.(S2.Vt)) ≤ √(eps())
        @test norm(S1.S - S2.S) ≤ √(eps())
    end

    @testset "Low rank" begin #Issue #42
        n = 10
        r = 2
        A = randn(n, r) * randn(r, n)
        S = svdvals(A)
        for nvals = 1 : r
            S1 = RandomizedLinAlg.rsvdvals(A, nvals, r - nvals)
            for i = 1 : nvals
                @test abs(S[i] - S1[i]) ≤ n^2 * r * eps()
            end
        end
    end

    @testset "rrange_adaptive" begin
        A = [1. 2 3; 4 5 6; 7 8 9]
        @test size(RandomizedLinAlg.rrange_adaptive(A, 3, 1e-3)) == (3,2)
    end

    @testset "rrange" begin
        A = [1. 2 3; 4 5 6; 7 8 9]
        @test_throws ArgumentError RandomizedLinAlg.rrange(A, 20)
    end
end