text:
module TestMPCRQ

using Base.Experimental: @sync
using ConcurrentCollections
using ConcurrentCollections.Implementations:
    DATA,
    IndirectMultiPolarityConcurrentRingQueueNode,
    MPCRQSlot,
    MPCRQ_CLOSED,
    MPCRQ_ENQUEUED,
    Waiter,
    denqueue!,
    tryput!
using ProgressLogging: @logprogress, @withprogress
using Test
using ..Utils: ⊏

function test_print_mpcrqslot()
    slot = MPCRQSlot(; index = 111, safe = false, polarity = DATA, storage = UInt32(0xaaa))
    str = sprint(show, "text/plain", slot)
    @test "MPCRQSlot" ⊏ str
    @test r"index *= *111" ⊏ str
    @test r"safe *= *false" ⊏ str
    @test r"polarity *=.*DATA" ⊏ str
    @test r"storage *= *0x0+aaa" ⊏ str
end

function test_close()
    crq = IndirectMultiPolarityConcurrentRingQueueNode{Int}(3)
    @testset for i in 1:8
        @test denqueue!(crq, i) === MPCRQ_ENQUEUED
    end
    @test denqueue!(crq, 9) === MPCRQ_CLOSED
    local ys, ylast
    @sync begin
        @async begin
            ys = [denqueue!(crq, Waiter{Int}())::Some{Int} for _ in 1:8]
            ylast = denqueue!(crq, Waiter{Int}())
        end
    end
    @test something.(ys) == 1:8
    @test ylast === MPCRQ_CLOSED
end

function unfair_sleep(seconds::Real)
    t0 = time_ns()
    ns = seconds * 1e9
    while time_ns() - t0 < ns
        GC.safepoint()
    end
end

function concurrent_denqueue!(
    crq::IndirectMultiPolarityConcurrentRingQueueNode,
    nitems::Integer,
    nsend::Integer,
    nrecv::Integer,
)
    @assert nsend > 0
    @assert nrecv > 0
    use_yield = nsend + nrecv > Threads.nthreads()
    # TODO: If `!use_yield`, make sure the tasks are spawned in different
    # threads and they are sticky.

    received = Vector{Int}[]
    senders = Task[]
    receivers = Task[]
    global SENDERS = senders
    global RECEIVERS = receivers
    ref = Threads.Atomic{Int}(0)
    @sync begin
        for offset in 1:nsend
            t = Threads.@spawn begin
                local y = nothing
                local i = 0
                for outer i in offset:nsend:nitems
                    local s = Threads.atomic_add!(ref, 1)
                    if s > nsend
                        while ref[] > 1
                            if use_yield
                                yield()
                            else
                                GC.safepoint()
                            end
                        end
                    end
                    local x = eltype(crq)(i)
                    y = denqueue!(crq, x)
                    y === MPCRQ_CLOSED && break
                    y === MPCRQ_ENQUEUED && continue
                    local ok = tryput!(y::Waiter, x)
                    @assert ok "tryput!(y::Waiter, x)"
                end
                return y, i
            end
            push!(senders, t)
        end
        for _ in 1:nrecv
            ys = Int[]
            push!(received, ys)
            t = Threads.@spawn begin
                local ys_nb = Int[]
                local ys_b = Int[]
                while true
                    local s = Threads.atomic_sub!(ref, 1)
                    if s < -nrecv
                        while ref[] < -1
                            if use_yield
                                yield()
                            else
                                GC.safepoint()
                            end
                        end
                    end
                    local w = Waiter{eltype(crq)}()
                    local y = denqueue!(crq, w)
                    if y === MPCRQ_CLOSED
                        return (; y, ys_nb, ys_b)
                    end
                    local i
                    if y isa Some
                        i = something(y)
                        push!(ys, i)
                        push!(ys_nb, i)
                    else
                        @assert y === MPCRQ_ENQUEUED
                        i = fetch(w)
                        push!(ys, i)
                        push!(ys_b, i)
                    end

                    if i > nitems - nrecv
                        return (; y, ys_nb, ys_b)
                    end
                end
            end
            push!(receivers, t)
        end
    end
    return received, senders, receivers
end

function check_consecutive(xs)
    notfound = Int[]
    dups = Int[]
    pre = xs[begin] - 1
    for x in xs
        e = pre + 1
        append!(notfound, e:x-1)
        append!(dups, x:e-1)
        pre = x
    end
    return (; notfound, dups)
end

function test_concurrent_push_pop(ntrials = 100)
    @withprogress name = "concurrent push-pop" begin
        @testset for trial in 1:ntrials
            @logprogress (trial - 1) / ntrials
            check_concurrent_push_pop()
        end
    end
end

function check_concurrent_push_pop()
    global received, notfound, dups, allreceived
    nsend = cld(Threads.nthreads(), 2)
    nrecv = max(1, Threads.nthreads() - nsend)
    crq = IndirectMultiPolarityConcurrentRingQueueNode{Int}(7)
    global CRQ = crq
    nitems = 2^20
    received, senders, receivers = concurrent_denqueue!(crq, nitems, nsend, nrecv)
    allreceived = reduce(vcat, received)

    @test length(allreceived) == nitems
    sort!(allreceived)
    (; notfound, dups) = check_consecutive(allreceived)
    @test length(allreceived) == nitems
    @test notfound == []
    @test dups == []
    @test allreceived == 1:nitems
end

function test_print()
    crq = IndirectMultiPolarityConcurrentRingQueueNode{Int}(3)
    str = sprint(show, "text/plain", crq)
    @test "MP-CRQ: " ⊏ str
    @test "empty" ⊏ str

    denqueue!(crq, 333)
    str = sprint(show, "text/plain", crq)
    @test "1 data item" ⊏ str

    denqueue!(crq, Waiter{Int}())  # pop
    w = Waiter{Int}()
    denqueue!(crq, w)
    str = sprint(show, "text/plain", crq)
    @test "1 waiter" ⊏ str
end

end  # module


text:
module CleverDicts

# The following two functions are overloaded for `MOI.VariableIndex` here
# it is the original use-case for `CleverDict`, and it would be type-piracy for
# solvers using `CleverDicts` to implement it themselves.

import MathOptInterface

function index_to_key(::Type{MathOptInterface.VariableIndex}, index::Int)
    return MathOptInterface.VariableIndex(index)
end

key_to_index(key::MathOptInterface.VariableIndex) = key.value

# Now, on with `CleverDicts`.

import OrderedCollections

"""
    CleverDict{K, V}

A smart storage type for managing sequential objects with non-decreasing integer
indices.

Provided no keys are deleted, the backing storage is a `Vector{V}`. Once a key
has been deleted, the backing storage switches to an `OrderedDict{K, V}`.

The i'th ordered element can be obtained with `c[LinearIndex(i)]`.

Note that querying a `LinearIndex` immediately after deleting a key via
`delete!` is very slow. (It requires a rebuild of an ordered list of variables.)

Store an item `val` using `add_item(c::CleverDict, val)`. `add_item` returns a
key corresponding to the stored item.

Overload the functions `index_to_key` and `key_to_index` to enable mappings
between the integer index of the vector and the dictionary key.

## Example

```julia
struct MyKey
    x::Int
end
index_to_key(::Type{MyKey}, i::Int) = MyKey(i)
key_to_index(key::MyKey) = key.x
```
"""
mutable struct CleverDict{K, V}
    last_index::Int
    vector::Union{Nothing, Vector{V}}
    dict::Union{Nothing, OrderedCollections.OrderedDict{K, V}}
    CleverDict{K, V}() where {K, V} = new{K, V}(0, V[], nothing)
end

"""
    index_to_key(::Type{K}, index::Int)

Create a new key associated with the integer value `index`.
"""
function index_to_key end

"""
    key_to_index(key::K)

Map `key` to an integer valued index, assuming that there have been no
deletions.
"""
function key_to_index end

"""
    add_item(c::CleverDict{K, V}, val::Val)::K where {K, V}

Set `val` in the next available key, and return that key.
"""
function add_item(c::CleverDict{K, V}, val::V)::K where {K, V}
    c.last_index += 1
    key = index_to_key(K, c.last_index)
    if c.dict === nothing
        push!(c.vector, val)
    else
        c.dict[key] = val
        # If there is a vector (e.g., because it has been rebuild for
        # `LinearIndex`), clear it.
        c.vector = nothing
    end
    return key
end

function Base.empty!(c::CleverDict{K, V})::Nothing where {K, V}
    c.vector = V[]
    c.last_index = 0
    c.dict = nothing
    return
end

function Base.getindex(c::CleverDict{K, V}, key::K)::V where {K, V}
    # Perform this `haskey` check up front to detect getting with keys that are
    # invalid (i.e., have previously been deleted).
    if !haskey(c, key)
        throw(KeyError(key))
    end
    # Case  I) no call to `Base.delete!`, so return the element:
    # Case II) `Base.delete!` must have been called, so return the element
    #          from the dictionary.
    return c.dict === nothing ? c.vector[key_to_index(key)] : c.dict[key]
end

function Base.setindex!(c::CleverDict{K, V}, val::V, key::K)::V where {K, V}
    # Perform this `haskey` check up front to detect setting with keys that are
    # invalid (i.e., have already been deleted). You can only call setindex!
    # with a key obtained from `new_key` that hasn't been deleted.
    if !haskey(c, key)
        throw(KeyError(key))
    elseif c.dict === nothing
        @assert c.vector !== nothing
        c.vector[key_to_index(key)] = val
    else
        c.dict[key] = val
    end
    return val
end

struct LinearIndex
    i::Int
end

function Base.getindex(c::CleverDict{K, V}, index::LinearIndex)::V where {K, V}
    if !(1 <= index.i <= length(c))
        throw(KeyError(index))
    end
    # Get the `index` linear element. If `c.vector` is currently `nothing`
    # (i.e., there has been a deletion), rebuild `c.vector`. This is a
    # trade-off: We could ensure `c.vector` is always updated, but this requires
    # a `splice!` in `delete!`, making deletions costly. However, it makes this
    # `getindex` operation trival because we would never have to rebuild the
    # vector.
    # The current implemented approach offers quick deletions, but an expensive
    # rebuild the first time you query a `LinearIndex` after a deletion or a new
    # key is added. Once the rebuild is done, there are quick queries until the
    # next deletion or addition. Thus, the worst-case is a user repeatedly
    # deleting a key and then querying a LinearIndex (e.g., getting the MOI
    # objective function).
    if c.vector === nothing
        c.vector = Vector{V}(undef, length(c))
        for (i, val) in enumerate(values(c.dict))
            c.vector[i] = val
        end
    end
    return c.vector[index.i]::V
end

function Base.delete!(c::CleverDict{K, V}, key::K)::Nothing where {K, V}
    if c.dict === nothing
        c.dict = OrderedCollections.OrderedDict{K, Union{Nothing, V}}()
        for (i, info) in enumerate(c.vector)
            c.dict[index_to_key(K, i)] = info
        end
    end
    delete!(c.dict, key)
    c.vector = nothing
    return
end

function Base.length(c::CleverDict)::Int
    return c.dict === nothing ? length(c.vector) : length(c.dict)
end

function Base.isempty(c::CleverDict)
    return c.dict === nothing ? isempty(c.vector) : isempty(c.dict)
end

Base.haskey(::CleverDict, key) = false
function Base.haskey(c::CleverDict{K, V}, key::K)::Bool where {K, V}
    if c.dict === nothing
        return 1 <= key_to_index(key) <= length(c.vector)
    else
        return haskey(c.dict, key)
    end
end

# Here, we implement the iterate functions for our `CleverDict`. If the backing
# datastructure is an `OrderedDict`, we just forward `iterate` to the dict. If
# it's the vector, we create a key-value pair so that `iterate` returns the same
# type regardless of the backing datastructure. To help inference, we annotate
# the return type.
#
# Also note that iterating an `OrderedDict` returns an `Int` state variable.
# This is identical to the type of the state variable that we return when
# iterating the vector, so we can add a type restriction on
# `iterate(c, s::Int)`.

function Base.iterate(
    c::CleverDict{K, V}
)::Union{Nothing, Tuple{Pair{K, V}, Int}} where {K, V}
    if c.dict === nothing
        @assert c.vector !== nothing
        if isempty(c.vector)
            return nothing
        end
        key = index_to_key(K, 1)
        return key => c.vector[1], 2
    else
        return iterate(c.dict)
    end
end

function Base.iterate(
    c::CleverDict{K, V}, s::Int
)::Union{Nothing, Tuple{Pair{K, V}, Int}} where {K, V}
    if c.dict === nothing
        @assert c.vector !== nothing
        if s > length(c.vector)
            return nothing
        end
        key = index_to_key(K, s)
        return key => c.vector[s], s + 1
    else
        return iterate(c.dict, s)
    end
end

function Base.values(c::CleverDict{K, V}) where {K, V}
    return c.dict === nothing ? c.vector : values(c.dict)
end

function Base.keys(c::CleverDict{K, V}) where {K, V}
    return c.dict === nothing ? index_to_key.(K, 1:length(c)) : keys(c.dict)
end

end

ext: jl


text:
#addprocs([("iriadmin@10.63.36.22", 1), ("iriadmin@10.63.36.23", 1)])
#addprocs(2)
#using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib
#@everywhere using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib

#root="/mnt/resource/analytics/models/Jenny-o"
#root="/mnt/resource/analytics/models/NatChoice"
#root="/mnt/resource/analytics/models/rev"
#root="/mnt/resource/analytics/models/ALL#10"
#root="/mnt/resource/analytics/models/ALL#11"
#root="/mnt/resource/analytics/models/CDW#6"
#root="/mnt/resource/analytics/models/HormelChili#8"
#root="/mnt/resource/analytics/models/Rev#8"
#root="/mnt/resource/analytics/models/Jenny-o"
using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib
root="/mnt/resource/analytics/models/ALL#10"
jmod_fname = root*"/dfd_model.json"
mod_fname = root*"/dfd_model.csv"
modelsDict = readModels(jmod_fname)
dfd = readtable(mod_fname,header=true);
poolit!(dfd,modelsDict[:factors])
for m in [modelsDict[:iocc], modelsDict[:idolocc], modelsDict[:ipen]]
    for r in m[:raneff]
        dfd[r] = relevel(dfd[r], "none")
    end
end
m=modelsDict[:iocc]

#GLM
f=genF(m[:y_var],m[:finalvars])
g= m[:Buyer_Pos_P1_is1] ? glm(f, dfd[(dfd[:buyer_pos_p1].==1),:] , m[:dist], m[:lnk]) : glm(f, dfd , m[:dist], m[:lnk])
#GLMM
f=genF(m[:y_var],setdiff(m[:finalvars],[:group]),[   r=m[:raneff][2]   ])
gmm = m[:Buyer_Pos_P1_is1] ? fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),:] , m[:dist], m[:lnk]),false,0)   : fit!(glmm(f, dfd, m[:dist], m[:lnk]),false,0)
raneffect(gmm)

ex_re = isBreak ? "& (df_data[:$ranef] .== \"$level\")" : "  "
cdict[:M] = cdict[:Mt] + cdict[:Mc]   # Recalculate totals as Test + the defaulted Control
cdict[:N] = cdict[:Nt] + cdict[:Nc]

#
cond="(dfd[:buyer_pos_p1].==1),"
r=string(coef(g)[1])*reduce(*,map((c,a)-> "+"*string(c)*"*dfd[$cond:"string(a)*"]" , coef(g)[2:end], setdiff(f.rhs.args,[:+,1])))
y=mean(exp(eval(parse(r))))                            # 1.238803537583254
raw = mean(dfd[(dfd[:buyer_pos_p1].==1),:trps_pos_p1]) # 1.2388035126234906
rawEQpred = round(mean(exp(eval(parse(r)))),6) == round(mean(dfd[(dfd[:buyer_pos_p1].==1),:trps_pos_p1]),6)

#MEAN SCORES
cond="(dfd[:buyer_pos_p1].==1)&(dfd[:group].==%%),"
r=string(coef(g)[1])*reduce(*,map((c,a)-> "+"*string(c)*"*dfd[$cond:"string(a)*"]" , coef(g)[2:end], setdiff(f.rhs.args,[:+,1])))
o_mean_score0=mean(exp(eval(parse( replace(r,"%%","0") ))))
o_mean_score1=mean(exp(eval(parse( replace(r,"%%","1") ))))

ex_re = "   "
isBreak=false
ex_re_single_col = isBreak ? ex_re[2:end]*"," : ""

mt=eval(parse("length(dfd[ (dfd[:group] .== 1) & (dfd[:buyer_pos_p1] .== 1 ) $ex_re ,1])" ))
n=eval(parse("length(dfd[ $ex_re_single_col 1])"))
nt=eval(parse("length(dfd[ (dfd[:group] .== 1) $ex_re ,1])"))
mc = eval(parse("length(dfd[ (dfd[:group] .== 0) & (dfd[:buyer_pos_p1] .== 1 ) ,1])"))
nc = eval(parse("length(dfd[ (dfd[:group] .== 0) ,1])"))

ctrl=o_mean_score0*(nc/n)
tst=o_mean_score1*(nt/n)

tst+ctrl == y  # it should be close

# ------------------------------------------------------------

NOTE: for fixed you're already got mean conditional on group
    SOOOO:  exp ::: mean(y estimate for only group X) *    B+-Z*stderr.....

b4*C8 = control proportions
B5*C7 = test proportions

  LB =   (B5*C7+    B4*C8*EXP( (B1-0.85*B2))  )   -
         (B4*C8+    B5*C7*EXP(-(B1-0.85*B2))  )
  UB =
         (B5*C7+    B4*C8*EXP( (B1+0.85*B2))  )   -
         (B4*C8+    B5*C7*EXP(-(B1+0.85*B2))  )

  LB =   (test+    ctrl*EXP( (B1-0.85*B2))  )   -
         (ctrl+    test*EXP(-(B1-0.85*B2))  )
  UB =
         (test+    ctrl*EXP( (B1+0.85*B2))  )   -
         (ctrl+    test*EXP(-(B1+0.85*B2))  )

LB

                Lb_pre = (   (mean_score1*(Mt/M))    -   (mean_score1*exp(-(B-(SE*zscore)))*(Mt/M))   )  +  ## ------------ Lower Bound ---------------
                         (   (mean_score0*exp((B-(SE*zscore)))*(Mc/M))    -   (mean_score0*(Mc/M))    )
                row[Symbol(zscore_key*"_lb")] = ( Lb_pre/mean_score0 ) * 100

UB             Ub_pre =  (     ( mean_score1*(Mt/M) )   -   ( mean_score1*exp(-(B+(SE*zscore)))*(Mt/M))   )  +  ## ------------ Upper Bound ---------------
                          (     ( mean_score0*exp((B+(SE*zscore)))*(Mc/M))  - (mean_score0*(Mc/M) )   )
                row[Symbol(zscore_key*"_ub")] = ( Ub_pre/mean_score0 ) * 100

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
#=
	This file implements FSA operations on other datatypes
	which are not defined in the respective package.
	This way, the benchmark functions stay simple
=#
import Base: +, -, .*, ./, *, /
ops = (:(+), :(-), :(.*), :(./), :(*), :(/))
for N=1:10, op in ops
	@eval begin
		@inline function $(op){T}(a::NTuple{$N, T}, b::NTuple{$N, T})
			@inbounds return tuple($([:($(op)(a[$i], b[$i])) for i=1:N]...))
		end
		@inline function $(op){T}(a::NTuple{$N, T}, b)
			@inbounds return tuple($([:($(op)(a[$i], b)) for i=1:N]...))
		end
		@inline function $(op){T}(a, b::NTuple{$N, T})
			@inbounds return tuple($([:($(op)(a, b[$i])) for i=1:N]...))
		end

		@inline function $(op){T<:VecElement}(a::NTuple{$N, T}, b::NTuple{$N, T})
			@inbounds return tuple($([:(VecElement($(op)(a[$i].value, b[$i].value))) for i=1:N]...))
		end
		@inline function $(op){T<:VecElement}(a::NTuple{$N, T}, b)
			@inbounds return tuple($([:(VecElement($(op)(a[$i].value, b))) for i=1:N]...))
		end
		@inline function $(op){T<:VecElement}(a, b::NTuple{$N, T})
			@inbounds return tuple($([:(VecElement($(op)(a, b[$i].value))) for i=1:N]...))
		end
	end
end
unary_ops = (:(+), :(-))
for N=1:10, op in ops
	@eval begin
		@inline function $(op){T}(a::NTuple{$N, T})
			@inbounds return tuple($([:($(op)(a[$i])) for i=1:N]...))::NTuple{$N, T}
		end
		@inline function $(op){T<:VecElement}(a::NTuple{$N, T})
			@inbounds return tuple($([:(VecElement($(op)(a[$i].value))) for i=1:N]...))::NTuple{$N, T}
		end
	end
end
Base.zero{T}(::VecElement{T}) = VecElement(zero(T))
function +(a::VecElement, b::VecElement)
	VecElement(a.value + b.value)
end
function *(a::VecElement, b::VecElement)
	VecElement(a.value * b.value)
end
function Base.sum(a::Tuple)
	reduce(+, a)
end

@inline function (*)(a::Vector, b::Vector)
	a .* b
end
@inline function (*)(a::StaticArrays.SVector, b::StaticArrays.SVector)
	a .* b
end
@inline function (*)(a::FixedSizeArrays.FixedVector, b::FixedSizeArrays.FixedVector)
	a .* b
end
@inline function (./)(a::SIMD.Vec, b::SIMD.Vec)
	a / b
end
@inline function (./){T}(a::SIMD.Vec, b::T)
	a / b
end

@inline function (/){N, T<:Integer, T2<:Integer}(a::SIMD.Vec{N, T}, b::SIMD.Vec{N, T2})
	TT = Base.promote_op(/, T, T2)
	SIMD.Vec{N, TT}(a) / SIMD.Vec{N, TT}(b)
end

@inline function (/){N,T1<:Integer,T2<:Integer}(a::SIMD.Vec{N, T1}, b::T2)
	TT =  Base.promote_op(/, T1, T2)
	SIMD.Vec{N, TT}(a) / TT(b)
end


text:
using Test

include("../src/simulator.jl")
using .simulator.fresnelltools

@testset "test_fresnelltools" begin

    @testset "FresnellBoundraryNoangle" begin
      # A randomized test verifying the angle-independent fresnell boundrary
      # is equivalent to the angle dependent at 90 degrees.
      for _ in 1:100
        n1 = rand(ComplexF64)
        n2 = rand(ComplexF64)

        independent = FresnellBoundrary(n1, n2)
        te_matrix, tm_matrix, θ = FresnellBoundrary(n1, n2, 0)

        @test independent == te_matrix == tm_matrix
        @test θ == 0
      end
    end

    @testset "FresnellBoundraryNonexisting" begin
      # Testing a variety of angles and refractive indices to check that nonexisting borders do no reflect.
      for _ in 1:100
        θ1 = rand() * π / 2
        n = rand(ComplexF64) * 2
        te, tm, θ2 = FresnellBoundrary(n, n, θ1)

        @test θ1 ≈ θ2
        @test round.(te, digits=5) ≈ round.(tm, digits=5) ≈ [1.0 0.0; 0.0 1.0]
      end
    end

    @testset "FresnellSlabTest" begin
      for _ in 1:100
        length = rand() * 10
        slab = FresnellSlab(1, 1, length, 0)

        @test slab ≈ [ℯ^(- length * 1im) 0; 0 ℯ^(- length * 1im)]
      end

      @test FresnellSlab(1, 2, 0, 0) == [1 0; 0 1]
    end

end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    StageDecompositionStructure

Vertical memory structure. Decomposes stochastic program by stages.

"""
struct StageDecompositionStructure{N, M, SP <: NTuple{M, AbstractScenarioProblems}} <: AbstractDecompositionStructure{N}
    decisions::Decisions{N}
    first_stage::JuMP.Model
    scenarioproblems::SP
    proxy::NTuple{N,JuMP.Model}

    function StageDecompositionStructure(decisions::Decisions{N}, scenarioproblems::NTuple{M,AbstractScenarioProblems}) where {N, M}
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        SP = typeof(scenarioproblems)
        proxy = ntuple(Val{N}()) do _
            Model()
        end
        return new{N,M,SP}(decisions, Model(), scenarioproblems, proxy)
    end
end

function StochasticStructure(decisions::Decisions{N}, scenario_types::ScenarioTypes{M}, instantiation::Union{StageDecomposition, DistributedStageDecomposition}) where {N, M}
    scenarioproblems = ntuple(Val(M)) do i
        ScenarioProblems(scenario_types[i], instantiation)
    end
    return StageDecompositionStructure(decisions, scenarioproblems)
end

function StochasticStructure(decisions::Decisions{N}, scenarios::NTuple{M, Vector{<:AbstractScenario}}, instantiation::Union{StageDecomposition, DistributedStageDecomposition}) where {N, M}
    scenarioproblems = ntuple(Val(M)) do i
        ScenarioProblems(scenarios[i], instantiation)
    end
    return StageDecompositionStructure(decisions, scenarioproblems)
end

# Base overloads #
# ========================== #
function Base.print(io::IO, structure::StageDecompositionStructure{N}) where N
    print(io, "Stage 1\n")
    print(io, "============== \n")
    print(io, structure.first_stage)
    for stage in 2:N
        print(io, "\nStage $stage\n")
        print(io, "============== \n")
        for (id, subproblem) in enumerate(subproblems(structure, stage))
            @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(structure, stage, id))
            print(io, subproblem)
            print(io, "\n")
        end
    end
end
function Base.print(io::IO, structure::StageDecompositionStructure{2})
    print(io, "First-stage \n")
    print(io, "============== \n")
    print(io, structure.first_stage)
    print(io, "\nSecond-stage \n")
    print(io, "============== \n")
    for (id, subproblem) in enumerate(subproblems(structure, 2))
        @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(structure, 2, id))
        print(io, subproblem)
        print(io, "\n")
    end
end
# ========================== #

# MOI #
# ========================== #
function MOI.get(structure::StageDecompositionStructure, attr::MOI.AbstractModelAttribute)
    return MOI.get(backend(structure.first_stage), attr)
end
function MOI.get(structure::StageDecompositionStructure, attr::MOI.AbstractVariableAttribute, index::MOI.VariableIndex)
    return MOI.get(backend(structure.first_stage), attr, index)
end
function MOI.get(structure::StageDecompositionStructure, attr::Type{MOI.VariableIndex}, name::String)
    return MOI.get(backend(structure.first_stage), attr, name)
end
function MOI.get(structure::StageDecompositionStructure, attr::MOI.AbstractConstraintAttribute, ci::MOI.ConstraintIndex)
    if _function_type(ci) <: SingleDecision
        # Need to map constraint
        con_ref = ConstraintRef(structure.first_stage, ci)
        return MOI.get(structure.first_stage, attr, con_ref)
    end
    return MOI.get(backend(structure.first_stage), attr, ci)
end
function MOI.get(structure::StageDecompositionStructure, attr::ScenarioDependentModelAttribute)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    return MOI.get(scenarioproblems(structure, attr.stage), attr)
end
function MOI.get(structure::StageDecompositionStructure, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    return MOI.get(scenarioproblems(structure, attr.stage), attr, index)
end
function MOI.get(structure::StageDecompositionStructure, attr::ScenarioDependentConstraintAttribute, ci::MOI.ConstraintIndex)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    return MOI.get(scenarioproblems(structure, attr.stage), attr, ci)
end

function MOI.set(structure::StageDecompositionStructure, attr::MOI.AbstractModelAttribute, value)
    MOI.set(backend(structure.first_stage), attr, value)
    return nothing
end
function MOI.set(structure::StageDecompositionStructure, attr::MOI.Silent, flag)
    # Silence master
    MOI.set(backend(structure.first_stage), attr, flag)
    # Silence subproblems
    MOI.set(scenarioproblems(structure), attr, flag)
    return nothing
end
function MOI.set(structure::StageDecompositionStructure, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    MOI.set(backend(structure.first_stage), attr, index, value)
    return nothing
end
function MOI.set(structure::StageDecompositionStructure, attr::MOI.AbstractConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    con_ref = ConstraintRef(structure.first_stage, ci)
    MOI.set(structure.first_stage, attr, con_ref, value)
    return nothing
end
function MOI.set(structure::StageDecompositionStructure, attr::ScenarioDependentModelAttribute, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    MOI.set(scenarioproblems(structure, attr.stage), attr, value)
    return nothing
end
function MOI.set(structure::StageDecompositionStructure, attr::ScenarioDependentVariableAttribute,
                 index::MOI.VariableIndex, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    MOI.set(scenarioproblems(structure, attr.stage), attr, index, value)
    return nothing
end
function MOI.set(structure::StageDecompositionStructure, attr::ScenarioDependentConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    MOI.set(scenarioproblems(structure, attr.stage), attr, ci, value)
    return nothing
end

function MOI.is_valid(structure::StageDecompositionStructure, index::MOI.VariableIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    return MOI.is_valid(backend(structure.first_stage), index)
end
function MOI.is_valid(structure::StageDecompositionStructure, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer)
    stage == 1 && error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return MOI.is_valid(scenarioproblems(structure, stage), index, scenario_index)
end
function MOI.is_valid(structure::StageDecompositionStructure, ci::MOI.ConstraintIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    return MOI.is_valid(backend(structure.first_stage), ci)
end
function MOI.is_valid(structure::StageDecompositionStructure{N}, ci::MOI.ConstraintIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return MOI.is_valid(scenarioproblems(structure, stage), ci, scenario_index)
end

function MOI.delete(structure::StageDecompositionStructure{N}, indices::Vector{MOI.VariableIndex}, stage::Integer) where N
    stage == 1 || error("No scenario index specified.")
    JuMP.delete(structure.first_stage, DecisionRef.(structure.first_stage, indices))
    # Remove known decisions
    for s in 2:N
        MOI.delete(scenarioproblems(structure, s), indices)
    end
    return nothing
end
function MOI.delete(structure::StageDecompositionStructure{N}, indices::Vector{MOI.VariableIndex}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), indices, scenario_index)
    return nothing
end
function MOI.delete(structure::StageDecompositionStructure, ci::MOI.ConstraintIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    MOI.delete(backend(structure.first_stage), ci)
    return nothing
end
function MOI.delete(structure::StageDecompositionStructure, cis::Vector{<:MOI.ConstraintIndex}, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    MOI.delete(backend(structure.first_stage), cis)
    return nothing
end
function MOI.delete(structure::StageDecompositionStructure{N}, ci::MOI.ConstraintIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), ci, scenario_index)
    return nothing
end
function MOI.delete(structure::StageDecompositionStructure{N}, cis::Vector{<:MOI.ConstraintIndex}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), cis, scenario_index)
    return nothing
end

# JuMP #
# ========================== #
function decision_dispatch(decision_function::Function,
                           structure::StageDecompositionStructure{N},
                           index::MOI.VariableIndex,
                           stage::Integer,
                           args...) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 || error("No scenario index specified.")
    dref = DecisionRef(structure.first_stage, index)
    return decision_function(dref, args...)
end
function decision_dispatch!(decision_function!::Function,
                            structure::StageDecompositionStructure{N},
                            index::MOI.VariableIndex,
                            stage::Integer,
                            args...) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 || error("No scenario index specified.")
    dref = DecisionRef(structure.first_stage, index)
    decision_function!(dref, args...)
    return nothing
end
function scenario_decision_dispatch(decision_function::Function,
                                    structure::StageDecompositionStructure{N},
                                    index::MOI.VariableIndex,
                                    stage::Integer,
                                    scenario_index::Integer,
                                    args...) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return scenario_decision_dispatch(decision_function,
                                      scenarioproblems(structure, stage),
                                      index,
                                      scenario_index,
                                      args...)
end
function scenario_decision_dispatch!(decision_function!::Function,
                                     structure::StageDecompositionStructure{N},
                                     index::MOI.VariableIndex,
                                     stage::Integer,
                                     scenario_index::Integer,
                                     args...) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    scenario_decision_dispatch!(decision_function!,
                                scenarioproblems(structure, stage),
                                index,
                                scenario_index,
                                args...)
    return nothing
end
function JuMP.fix(structure::StageDecompositionStructure{N}, index::MOI.VariableIndex, stage::Integer, val::Number) where N
    dref = DecisionRef(structure.first_stage, index)
    fix(dref, val)
    update_decision_state!(scenarioproblems(structure), index, Known)
    return nothing
end
function JuMP.fix(structure::StageDecompositionStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer, val::Number) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    fix(scenarioproblems(structure, stage), index, scenario_index, val)
    return nothing
end
function JuMP.unfix(structure::StageDecompositionStructure{N}, index::MOI.VariableIndex, stage::Integer) where N
    dref = DecisionRef(structure.first_stage, index)
    unfix(dref)
    return nothing
end
function JuMP.unfix(structure::StageDecompositionStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    unfix(scenarioproblems(structure, stage), index, scenario_index)
    return nothing
end
function JuMP.set_objective_sense(structure::StageDecompositionStructure, stage::Integer, sense::MOI.OptimizationSense)
    if stage == 1
        MOI.set(structure, MOI.ObjectiveSense(), sense)
    else
        # Every sub-objective in the given stage should be changed
        MOI.set(scenarioproblems(structure), MOI.ObjectiveSense(), sense)
    end
    return nothing
end
function JuMP.objective_function_type(structure::StageDecompositionStructure)
    return jump_function_type(structure.first_stage,
                              MOI.get(structure, MOI.ObjectiveFunctionType()))
end

function JuMP.objective_function(structure::StageDecompositionStructure, FunType::Type{<:AbstractJuMPScalar})
    MOIFunType = moi_function_type(FunType)
    func = MOI.get(structure.first_stage,
                   MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
    return JuMP.jump_function(structure, 1, func)
end

function JuMP.objective_function(structure::StageDecompositionStructure, stage::Integer, FunType::Type{<:AbstractJuMPScalar})
    if stage == 1
        return objective_function(structure, FunType)
    else
        return objective_function(structure.proxy, FunType)
    end
end

function JuMP._moi_optimizer_index(structure::StageDecompositionStructure, index::VI)
    return decision_index(backend(structure.first_stage), index)
end
function JuMP._moi_optimizer_index(structure::StageDecompositionStructure, index::VI, scenario_index::Integer)
    return JuMP._moi_optimizer_index(scenarioproblems(structure), index, scenario_index)
end
function JuMP._moi_optimizer_index(structure::StageDecompositionStructure, ci::CI)
    return decision_index(backend(structure.first_stage), ci)
end

function JuMP.set_objective_coefficient(structure::StageDecompositionStructure{N}, index::VI, var_stage::Integer, stage::Integer, coeff::Real) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    var_stage <= stage || error("Can only modify coefficient in current stage of decision or subsequent stages from where decision is taken.")
    if var_stage == 1
        if stage == 1
            dref = DecisionRef(structure.first_stage, index)
            set_objective_coefficient(structure.first_stage, dref, coeff)
        else
            set_objective_coefficient(scenarioproblems(structure, stage), index, coeff)
        end
    else
        for scenario_index in 1:num_scenarios(structure, stage)
            set_objective_coefficient(scenarioproblems(structure, stage), index, scenario_index, coeff)
        end
    end
    return nothing
end
function JuMP.set_objective_coefficient(structure::StageDecompositionStructure{N}, index::VI, var_stage::Integer, stage::Integer, scenario_index::Integer, coeff::Real) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    if var_stage == 1
        set_objective_coefficient(scenarioproblems(structure, stage), index, scenario_index, coeff)
    else
        set_objective_coefficient(scenarioproblems(structure, stage), index, scenario_index, coeff)
    end
    return nothing
end

function JuMP.set_normalized_coefficient(structure::StageDecompositionStructure,
                                         ci::CI{F,S},
                                         index::VI,
                                         value) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    MOI.modify(backend(structure.first_stage), ci,
               DecisionCoefficientChange(index, convert(T, value)))
    return nothing
end
function JuMP.set_normalized_coefficient(structure::StageDecompositionStructure{N},
                                         ci::CI{F,S},
                                         index::VI,
                                         var_stage::Integer,
                                         stage::Integer,
                                         scenario_index::Integer,
                                         value) where {N, T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    set_normalized_coefficient(scenarioproblems(structure, stage), ci, index, scenario_index, value)
    return nothing
end

function JuMP.normalized_coefficient(structure::StageDecompositionStructure{N},
                                     ci::CI{F,S},
                                     index::VI,
                                     var_stage::Integer,
                                     stage::Integer,
                                     scenario_index::Integer) where {N, T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return normalized_coefficient(scenarioproblems(structure, stage), ci, index, scenario_index)
end

function JuMP.set_normalized_rhs(structure::StageDecompositionStructure,
                                 ci::CI{F,S},
                                 value) where {T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    MOI.set(backend(structure.first_stage), MOI.ConstraintSet(), ci,
            S(convert(T, value)))
    return nothing
end
function JuMP.set_normalized_rhs(structure::StageDecompositionStructure{N},
                                 ci::CI{F,S},
                                 stage::Integer,
                                 scenario_index::Integer,
                                 value) where {N,
                                               T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    set_normalized_rhs(scenarioproblems(structure, stage), ci, scenario_index, value)
    return nothing
end

function DecisionRef(structure::StageDecompositionStructure, index::VI)
    return DecisionRef(structure.first_stage, index)
end

function JuMP.jump_function(structure::StageDecompositionStructure{N},
                            stage::Integer,
                            f::MOI.AbstractFunction) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    if stage == 1
        return JuMP.jump_function(structure.first_stage, f)
    else
        return JuMP.jump_function(structure.proxy[stage], f)
    end
end
function JuMP.jump_function(structure::StageDecompositionStructure{N},
                            stage::Integer,
                            scenario_index::Integer,
                            f::MOI.AbstractFunction) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return JuMP.jump_function(structure.proxy[stage], f)
end

# Getters #
# ========================== #
function structure_name(structure::StageDecompositionStructure)
    return "Stage-decomposition"
end

deferred_first_stage(structure::StageDecompositionStructure, ::Val{1}) = num_variables(first_stage(structure)) == 0

function decision(structure::StageDecompositionStructure{N}, index::MOI.VariableIndex, stage::Integer) where N
    stage == 1 || error("No scenario index specified.")
    return decision(structure.decisions, stage, index)
end
function decision(structure::StageDecompositionStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return decision(scenarioproblems(structure, stage), index, scenario_index)
end
# ========================== #

# Setters
# ========================== #
function update_known_decisions!(structure::StageDecompositionStructure)
    update_known_decisions!(structure.first_stage)
    return nothing
end


text:
# the current julia logo

"""
Draw the logo:

    julialogo(colored=true, masking=false)

Just create paths for the logo but don't apply :

    julialogo(false, true)

Typical usage:

```
using Luxor, Colors
# find this file, eg:
include("~/.julia/v0.5/Luxor/test/julia-logo.jl"))
Drawing(1000, 1000, "/tmp/heart-julia.pdf")
origin()
background("white")
julialogo()
finish()
preview()
```

"""

function julialogo(colored=true, masking=false)
    # set colors if colored is true
    # do the fills if masking is false
    colored && setcolor("black")
    move(72.872, 177.311)
    curve(72.872, 184.847, 72.024, 190.932, 70.329, 195.567)
    curve( 68.633, 200.202, 66.222, 203.8, 63.094, 206.362)
    curve( 59.967, 208.925, 56.217, 210.639, 51.846, 211.506)
    curve( 47.476, 212.373, 42.615, 212.806, 37.264, 212.806)
    curve( 30.029, 212.806, 24.49, 211.675, 20.647, 209.414)
    curve( 16.803, 207.154, 14.881, 204.441, 14.881, 201.275)
    curve( 14.881, 198.638, 15.955, 196.415, 18.103, 194.606)
    curve( 20.251, 192.797, 23.134, 191.893, 26.751, 191.893)
    curve( 29.464, 191.893, 31.631, 192.628, 33.251, 194.097)
    curve( 34.872, 195.567, 36.209, 197.018, 37.264, 198.449)
    curve( 38.47, 200.032, 39.487, 201.087, 40.316, 201.615)
    curve( 41.145, 202.142, 41.899, 202.406, 42.577, 202.406)
    curve( 44.009, 202.406, 45.102, 201.558, 45.855, 199.862)
    curve( 46.609, 198.167, 46.985, 194.87, 46.985, 189.971)
    line(46.985, 97.051)
    line(72.872, 89.929)
    line(72.872, 177.311)
    masking ? newsubpath() : fill()

    move(109.739, 92.416)
    line(109.739, 152.215)
    curve( 109.739, 153.874, 110.06, 155.437, 110.7, 156.907)
    curve( 111.341, 158.376, 112.226, 159.639, 113.357, 160.694)
    curve( 114.487, 161.749, 115.806, 162.596, 117.313, 163.237)
    curve( 118.821, 163.878, 120.441, 164.198, 122.174, 164.198)
    curve( 124.134, 164.198, 126.36, 163.101, 129.069, 161.202)
    curve( 133.363, 158.194, 135.965, 156.127, 135.965, 153.685)
    curve( 135.965, 153.099, 135.965, 92.416, 135.965, 92.416)
    line(161.738, 92.416)
    line(161.738, 177.311)
    line(135.965, 177.311)
    line(135.965, 169.398)
    curve( 132.574, 172.262, 128.956, 174.56, 125.113, 176.293)
    curve( 121.27, 178.027, 117.539, 178.893, 113.922, 178.893)
    curve( 109.702, 178.893, 105.783, 178.196, 102.165, 176.802)
    curve( 98.548, 175.408, 95.383, 173.505, 92.67, 171.093)
    curve( 89.957, 168.682, 87.828, 165.856, 86.283, 162.615)
    curve( 84.739, 159.375, 83.966, 155.908, 83.966, 152.215)
    line(83.966, 92.416)
    line(109.739, 92.416)
    masking ? newsubpath() : fill()

    move(197.881, 177.311)
    line(172.221, 177.311)
    line(172.221, 58.278)
    line(197.881, 51.156)
    line(197.881, 177.311)
    masking ? newsubpath() : fill()

    move(208.603, 97.051)
    line(234.376, 89.929)
    line(234.376, 177.311)
    line(208.603, 177.311)
    line(208.603, 97.051)
    masking ? newsubpath() : fill()

    move(288.225, 133.451)
    curve(285.738, 134.506, 283.232, 135.731, 280.707, 137.124)
    curve(278.183, 138.519, 275.884, 140.045, 273.812, 141.703)
    curve(271.74, 143.361, 270.062, 145.132, 268.782, 147.015)
    curve(267.501, 148.9, 266.86, 150.859, 266.86, 152.894)
    curve(266.86, 154.476, 267.067, 156.002, 267.481, 157.472)
    curve(267.896, 158.941, 268.48, 160.204, 269.234, 161.259)
    curve(269.988, 162.314, 270.816, 163.162, 271.721, 163.802)
    curve(272.625, 164.443, 273.605, 164.763, 274.66, 164.763)
    curve(276.77, 164.763, 278.899, 164.123, 281.047, 162.841)
    curve(283.194, 161.56, 285.587, 159.94, 288.225, 157.981)
    line(288.225, 133.451)
    newsubpath()
    move(314.111, 177.311)
    line(288.225, 177.311)
    line(288.225, 170.528)
    curve(286.792, 171.734, 285.399, 172.846, 284.042, 173.863)
    curve(282.686, 174.88, 281.16, 175.766, 279.464, 176.519)
    curve(277.768, 177.273, 275.866, 177.857, 273.755, 178.272)
    curve(271.645, 178.686, 269.158, 178.893, 266.295, 178.893)
    curve(262.375, 178.893, 258.853, 178.328, 255.725, 177.198)
    curve(252.597, 176.067, 249.941, 174.522, 247.756, 172.563)
    curve(245.571, 170.604, 243.893, 168.287, 242.725, 165.611)
    curve(241.558, 162.936, 240.973, 160.015, 240.973, 156.85)
    curve(240.973, 153.61, 241.595, 150.671, 242.838, 148.033)
    curve(244.082, 145.396, 245.778, 143.022, 247.925, 140.911)
    curve(250.073, 138.801, 252.579, 136.917, 255.443, 135.259)
    curve(258.306, 133.601, 261.377, 132.075, 264.655, 130.681)
    curve(267.934, 129.287, 271.344, 128.006, 274.886, 126.838)
    curve(278.427, 125.67, 281.932, 124.558, 285.399, 123.503)
    line(288.225, 122.825)
    line(288.225, 114.459)
    curve(288.225, 109.034, 287.188, 105.19, 285.116, 102.929)
    curve(283.044, 100.668, 280.274, 99.538, 276.807, 99.538)
    curve(272.738, 99.538, 269.912, 100.518, 268.329, 102.477)
    curve(266.747, 104.436, 265.955, 106.81, 265.955, 109.599)
    curve(265.955, 111.181, 265.786, 112.726, 265.447, 114.233)
    curve(265.108, 115.741, 264.523, 117.059, 263.695, 118.19)
    curve(262.866, 119.32, 261.679, 120.225, 260.134, 120.903)
    curve(258.589, 121.581, 256.649, 121.92, 254.312, 121.92)
    curve(250.695, 121.92, 247.756, 120.884, 245.495, 118.812)
    curve(243.234, 116.739, 242.104, 114.12, 242.104, 110.955)
    curve(242.104, 108.016, 243.102, 105.284, 245.099, 102.76)
    curve(247.097, 100.235, 249.79, 98.068, 253.182, 96.26)
    curve(256.573, 94.451, 260.492, 93.019, 264.938, 91.964)
    curve(269.384, 90.909, 274.094, 90.382, 279.068, 90.382)
    curve(285.173, 90.382, 290.429, 90.928, 294.838, 92.021)
    curve(299.246, 93.114, 302.883, 94.677, 305.746, 96.712)
    curve(308.609, 98.747, 310.72, 101.196, 312.076, 104.06)
    curve(313.433, 106.923, 314.111, 110.127, 314.111, 113.668)
    line(314.111, 177.311)
    masking ? newsubpath() : fill()

    colored && setcolor(0.835, 0.388, 0.361)
    move(222.772, 84.091)
    curve(213.936, 84.091, 206.772, 76.928, 206.772, 68.091)
    curve(206.772, 59.255, 213.936, 52.091, 222.772, 52.091)
    curve(231.609, 52.091, 238.772, 59.255, 238.772, 68.091)
    curve(238.772, 76.928, 231.609, 84.091, 222.772, 84.091)
    masking ? newsubpath() : fill()

    colored && setcolor(0.796, 0.235, 0.2)
    move(222.772, 50.591)
    curve(232.437, 50.591, 240.272, 58.426, 240.272, 68.091)
    curve(240.272, 77.756, 232.437, 85.591, 222.772, 85.591)
    curve(213.107, 85.591, 205.272, 77.756, 205.272, 68.091)
    curve(205.272, 58.426, 213.107, 50.591, 222.772, 50.591)
    newsubpath()
    move(222.772, 53.591)
    curve(214.764, 53.591, 208.272, 60.083, 208.272, 68.091)
    curve(208.272, 76.099, 214.764, 82.591, 222.772, 82.591)
    curve(230.781, 82.591, 237.272, 76.099, 237.272, 68.091)
    curve(237.272, 60.083, 230.781, 53.591, 222.772, 53.591)
    masking ? newsubpath() : fill()

    colored && setcolor(0.4, 0.51, 0.878)
    move(60.454, 84.091)
    curve(51.618, 84.091, 44.454, 76.928, 44.454, 68.091)
    curve(44.454, 59.255, 51.618, 52.091, 60.454, 52.091)
    curve(69.291, 52.091, 76.454, 59.255, 76.454, 68.091)
    curve(76.454, 76.928, 69.291, 84.091, 60.454, 84.091)

    masking ? newsubpath() : fill()
    colored && setcolor(0.251, 0.388, 0.847)
    move(60.454, 50.591)
    curve(70.119, 50.591, 77.954, 58.426, 77.954, 68.091)
    curve(77.954, 77.756, 70.119, 85.591, 60.454, 85.591)
    curve(50.789, 85.591, 42.954, 77.756, 42.954, 68.091)
    curve(42.954, 58.426, 50.789, 50.591, 60.454, 50.591)
    newsubpath()
    move(60.454, 53.591)
    curve(52.446, 53.591, 45.954, 60.083, 45.954, 68.091)
    curve(45.954, 76.099, 52.446, 82.591, 60.454, 82.591)
    curve(68.462, 82.591, 74.954, 76.099, 74.954, 68.091)
    curve(74.954, 60.083, 68.462, 53.591, 60.454, 53.591)
    masking ? newsubpath() : fill()

    colored && setcolor(0.667, 0.475, 0.757)
    move(264.821, 84.091)
    curve(255.984, 84.091, 248.821, 76.928, 248.821, 68.091)
    curve(248.821, 59.255, 255.984, 52.091, 264.821, 52.091)
    curve(273.658, 52.091, 280.821, 59.255, 280.821, 68.091)
    curve(280.821, 76.928, 273.658, 84.091, 264.821, 84.091)
    masking ? newsubpath() : fill()

    colored && setcolor(0.584, 0.345, 0.698)
    move(264.821, 50.591)
    curve(274.486, 50.591, 282.321, 58.426, 282.321, 68.091)
    curve(282.321, 77.756, 274.486, 85.591, 264.821, 85.591)
    curve(255.156, 85.591, 247.321, 77.756, 247.321, 68.091)
    curve(247.321, 58.426, 255.156, 50.591, 264.821, 50.591)
    newsubpath()
    move(264.821, 53.591)
    curve(256.813, 53.591, 250.321, 60.083, 250.321, 68.091)
    curve(250.321, 76.099, 256.813, 82.591, 264.821, 82.591)
    curve(272.829, 82.591, 279.321, 76.099, 279.321, 68.091)
    curve(279.321, 60.083, 272.829, 53.591, 264.821, 53.591)
    masking ? newsubpath() : fill()

    colored && setcolor(0.376, 0.678, 0.318)
    move(243.799, 47.672)
    curve(234.962, 47.672, 227.799, 40.509, 227.799, 31.672)
    curve(227.799, 22.836, 234.962, 15.672, 243.799, 15.672)
    curve(252.636, 15.672, 259.799, 22.836, 259.799, 31.672)
    curve(259.799, 40.509, 252.636, 47.672, 243.799, 47.672)
    masking ? newsubpath() : fill()

    colored && setcolor(0.22, 0.596, 0.149)
    move(243.799, 14.172)
    curve(253.464, 14.172, 261.299, 22.007, 261.299, 31.672)
    curve(261.299, 41.337, 253.464, 49.172, 243.799, 49.172)
    curve(234.134, 49.172, 226.299, 41.337, 226.299, 31.672)
    curve(226.299, 22.007, 234.134, 14.172, 243.799, 14.172)
    newsubpath()
    move(243.799, 17.172)
    curve(235.791, 17.172, 229.299, 23.664, 229.299, 31.672)
    curve(229.299, 39.681, 235.791, 46.172, 243.799, 46.172)
    curve(251.807, 46.172, 258.299, 39.681, 258.299, 31.672)
    curve(258.299, 23.664, 251.807, 17.172, 243.799, 17.172)
    if ! masking
      fill()
    end
end


text:
module GridapTypesTests

using Test
using Gridap.Helpers

#import Gridap.Helpers: operate

struct Foo{A} <: GridapType
  bar::A
end

object = Foo(nothing)

@test sprint(show,"text/plain",object) == "$(nameof(typeof(object)))()"

#function operate(op,a::Foo)
#  Foo(op(a.bar))
#end
#
#function operate(op,a::Foo,b::Foo)
#  Foo(op(a.bar,b.bar))
#end
#
#function operate(op,a::Number,b::Foo)
#  operate(op,Foo(a),b)
#end
#
#function operate(op,a::Foo,b::Number)
#  operate(op,a,Foo(b))
#end
#
#a = Foo(1)
#
#b = -a
#@test b.bar == -1
#
#c = a + b
#
#@test c.bar == 0
#
#c = a - 2
#@test c.bar == -1
#
#c = 2 - a
#@test c.bar == 1
#
#using Gridap.TensorValues
#
#c = inner(2,a)
#@test c.bar == 2
#
#c = a / b
#@test c.bar == -1

end # module


text:
using Test
using LinearAlgebra
using UnPack
using MOONS.Grids
using MOONS.Fields

"""
    convergence_rate(err, Δh)

Estimate convergence rate given
vectors `err` and `Δh`
```
err = C Δh^p+ H.O.T
err_k ≈ C Δh_k^p
err_k/err_m ≈ Δh_k^p/Δh_m^p
log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p
```
"""
convergence_rate(err, Δh) =
    [log(err[i]/err[i-1])/log(Δh[i]/Δh[i-1]) for i in 2:length(Δh)]

@testset "Face -> Center" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    err, Δh = zeros(length(nr)), zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n; warpfun=Roberts_both, args=(FT(1.01),))
        Δh[k] = grid.c[1].n.Δh[1]
        face = CellFace(grid, 1)
        cent_exact = CellCenter(grid)
        cent = CellCenter(grid)
        @inbounds for local_grid in GridField(grid,face)
            @unpack ijk,x,y,z = local_grid
            face[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,cent_exact)
            @unpack ijk,x,y,z = local_grid
            cent_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        interp!(cent, face, grid)
        err[k] = norm(cent .- cent_exact) / length(cent)
    end
    conv = convergence_rate(err, Δh)
    @test conv[1] ≈ 0.9049443448120211
    @test conv[2] ≈ 2.55205758617018
    @test conv[3] ≈ 2.992429087432938
end

@testset "Center -> Face" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    err, Δh = zeros(length(nr)), zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n; warpfun=Roberts_both, args=(FT(1.01),))
        Δh[k] = grid.c[1].n.Δh[1]
        cent = CellCenter(grid)
        face_exact = CellFace(grid, 1)
        face = CellFace(grid, 1)
        @inbounds for local_grid in GridField(grid,cent)
            @unpack ijk,x,y,z = local_grid
            cent[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,face_exact)
            @unpack ijk,x,y,z = local_grid
            face_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        face_exact[1,:,:] .= 0 # interpolation does not populate this
        face_exact[end,:,:] .= 0 # interpolation does not populate this
        interp!(face, cent, grid)
        err[k] = norm(face .- face_exact) / length(face)
    end
    conv = convergence_rate(err, Δh)
    @test conv[1] ≈ 1.8334675036878896
    @test conv[2] ≈ 2.521189589590259
    @test conv[3] ≈ 2.9686936734418015
end

@testset "Corner -> Edge" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    err, Δh = zeros(length(nr)), zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n; warpfun=Roberts_both, args=(FT(1.01),))
        Δh[k] = grid.c[1].n.Δh[1]
        corn = CellCorner(grid)
        edge_exact = CellEdge(grid, 1)
        edge = CellEdge(grid, 1)
        @inbounds for local_grid in GridField(grid,corn)
            @unpack ijk,x,y,z = local_grid
            corn[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,edge_exact)
            @unpack ijk,x,y,z = local_grid
            edge_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        interp!(edge, corn, grid)
        err[k] = norm(edge .- edge_exact) / length(edge)
    end
    conv = convergence_rate(err, Δh)
    @test conv[1] ≈ 1.2685430731433063
    @test conv[2] ≈ 2.482553166205038
    @test conv[3] ≈ 2.9528060057161647
end

@testset "Edge -> Corner" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    err, Δh = zeros(length(nr)), zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n; warpfun=Roberts_both, args=(FT(1.01),))
        Δh[k] = grid.c[1].n.Δh[1]
        edge = CellEdge(grid, 1)
        corn_exact = CellCorner(grid)
        corn = CellCorner(grid)
        @inbounds for local_grid in GridField(grid,edge)
            @unpack ijk,x,y,z = local_grid
            edge[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,corn_exact)
            @unpack ijk,x,y,z = local_grid
            corn_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        corn_exact[1,:,:] .= 0 # interpolation does not populate this
        corn_exact[end,:,:] .= 0 # interpolation does not populate this
        interp!(corn, edge, grid)
        err[k] = norm(corn .- corn_exact) / length(corn)
    end
    conv = convergence_rate(err, Δh)
    @test conv[1] ≈ 2.1970662320191745
    @test conv[2] ≈ 2.4516851696251174
    @test conv[3] ≈ 2.9290705917250293
end

@testset "Edge -> Face" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    err, Δh = zeros(length(nr)), zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n; warpfun=Roberts_both, args=(FT(1.01),))
        Δh[k] = grid.c[1].n.Δh[1]
        edge = CellEdge(grid, 2)
        face_exact = CellFace(grid, 1)
        face = CellFace(grid, 1)
        @inbounds for local_grid in GridField(grid,edge)
            @unpack ijk,x,y,z = local_grid
            edge[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,face_exact)
            @unpack ijk,x,y,z = local_grid
            face_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        interp!(face, edge, grid)
        err[k] = norm(face .- face_exact) / length(face)
    end
    conv = convergence_rate(err, Δh)
    @test conv[1] ≈ 1.0867437089776637
    @test conv[2] ≈ 2.5173053761876085
    @test conv[3] ≈ 2.972617546574552
end

@testset "Face -> Edge" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    err, Δh = zeros(length(nr)), zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n; warpfun=Roberts_both, args=(FT(1.01),))
        Δh[k] = grid.c[1].n.Δh[1]
        face = CellFace(grid, 1)
        edge_exact = CellEdge(grid, 2)
        edge = CellEdge(grid, 2)
        @inbounds for local_grid in GridField(grid,face)
            @unpack ijk,x,y,z = local_grid
            face[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,edge_exact)
            @unpack ijk,x,y,z = local_grid
            edge_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        interp!(edge, face, grid)
        edge_exact[:,:,1] .= 0 # interpolation does not populate this
        edge_exact[:,:,end] .= 0 # interpolation does not populate this
        err[k] = norm(edge .- edge_exact) / length(edge)
    end
    conv = convergence_rate(err, Δh)
    @test conv[1] ≈ 2.0152668678535317
    @test conv[2] ≈ 2.4864373796076884
    @test conv[3] ≈ 2.9488821325834156
end


text:
export BrickletThermocoupleTemperatureCallbackThreshold
struct BrickletThermocoupleTemperatureCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletThermocoupleConfiguration
struct BrickletThermocoupleConfiguration
    averaging::Integer
    thermocouple_type::Integer
    filter::Integer
end

export BrickletThermocoupleErrorState
struct BrickletThermocoupleErrorState
    over_under::Bool
    open_circuit::Bool
end

export BrickletThermocoupleIdentity
struct BrickletThermocoupleIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletThermocouple
"""
Measures temperature with thermocouples
"""
mutable struct BrickletThermocouple <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletThermocouple(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_thermocouple")
        deviceInternal = package.BrickletThermocouple(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of the thermocouple.

If you want to get the temperature periodically, it is recommended
to use the :cb:`Temperature` callback and set the period with
:func:`Set Temperature Callback Period`.
"""
function get_temperature(device::BrickletThermocouple)
    return device.deviceInternal.get_temperature()
end

export set_temperature_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Temperature` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Temperature` callback is only triggered if the temperature has changed
since the last triggering.
"""
function set_temperature_callback_period(device::BrickletThermocouple, period)
    device.deviceInternal.set_temperature_callback_period(period)
end

export get_temperature_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Temperature Callback Period`.
"""
function get_temperature_callback_period(device::BrickletThermocouple)
    return device.deviceInternal.get_temperature_callback_period()
end

export set_temperature_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Temperature Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the temperature is *outside* the min and max values"
 "'i'",    "Callback is triggered when the temperature is *inside* the min and max values"
 "'<'",    "Callback is triggered when the temperature is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the temperature is greater than the min value (max is ignored)"
"""
function set_temperature_callback_threshold(device::BrickletThermocouple, option, min, max)
    device.deviceInternal.set_temperature_callback_threshold(option, min, max)
end

export get_temperature_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Temperature Callback Threshold`.
"""
function get_temperature_callback_threshold(device::BrickletThermocouple)
    return device.deviceInternal.get_temperature_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Temperature Reached`

is triggered, if the threshold

* :func:`Set Temperature Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletThermocouple, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletThermocouple)
    return device.deviceInternal.get_debounce_period()
end

export set_configuration
"""
    $(SIGNATURES)

You can configure averaging size, thermocouple type and frequency
filtering.

Available averaging sizes are 1, 2, 4, 8 and 16 samples.

As thermocouple type you can use B, E, J, K, N, R, S and T. If you have a
different thermocouple or a custom thermocouple you can also use
G8 and G32. With these types the returned value will not be in °C/100,
it will be calculated by the following formulas:

* G8: ``value = 8 * 1.6 * 2^17 * Vin``
* G32: ``value = 32 * 1.6 * 2^17 * Vin``

where Vin is the thermocouple input voltage.

The frequency filter can be either configured to 50Hz or to 60Hz. You should
configure it according to your utility frequency.

The conversion time depends on the averaging and filter configuration, it can
be calculated as follows:

* 60Hz: ``time = 82 + (samples - 1) * 16.67``
* 50Hz: ``time = 98 + (samples - 1) * 20``
"""
function set_configuration(device::BrickletThermocouple, averaging, thermocouple_type, filter)
    device.deviceInternal.set_configuration(averaging, thermocouple_type, filter)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletThermocouple)
    return device.deviceInternal.get_configuration()
end

export get_error_state
"""
    $(SIGNATURES)

Returns the current error state. There are two possible errors:

* Over/Under Voltage and
* Open Circuit.

Over/Under Voltage happens for voltages below 0V or above 3.3V. In this case
it is very likely that your thermocouple is defective. An Open Circuit error
indicates that there is no thermocouple connected.

You can use the :cb:`Error State` callback to automatically get triggered
when the error state changes.
"""
function get_error_state(device::BrickletThermocouple)
    return device.deviceInternal.get_error_state()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletThermocouple)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletThermocouple, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
module UniqueVectors

include("delegate.jl")

import Base: copy, in, getindex, findfirst, findlast, length, size, isempty, iterate, empty!, push!, pop!, setindex!, getindex, indexin, findnext, findprev, findall, count, allunique, unique, unique!, permute!, invpermute!

EqualTo = Base.Fix2{typeof(isequal)}

abstract type AbstractUniqueVector{T} <: AbstractVector{T} end

struct UniqueVector{T} <: AbstractUniqueVector{T}
    items::Vector{T}
    lookup::Dict{T,Int}

    UniqueVector{T}() where {T} = new(T[], Dict{T,Int}())

    function UniqueVector{T}(items) where {T}
        uv = new(Vector{T}(items), Dict{T,Int}())
        sizehint!(uv.lookup, length(uv.items))
        for (i, item) in enumerate(uv.items)
            if item in uv
                throw(ArgumentError("cannot construct UniqueVector with duplicate items"))
            end
            uv.lookup[item] = i
        end
        @assert length(uv.items) == length(uv.lookup)
        return uv
    end
end

UniqueVector(items::AbstractVector{T}) where {T} = UniqueVector{T}(items)

@deprecate UniqueVector(items) UniqueVector(collect(items))

copy(uv::UniqueVector) = UniqueVector(uv.items)

@delegate UniqueVector.items [ length, size, isempty, getindex, iterate ]

function empty!(uv::UniqueVector)
    # NOTE: does not provide any exception safety guarantee
    empty!(uv.items)
    empty!(uv.lookup)
    return uv
end

in(item::T, uv::UniqueVector{T}) where {T} = haskey(uv.lookup, item)
in(item, uv::UniqueVector{T}) where {T} = in(convert(T, item), uv)

findfirst(p::EqualTo{<:T}, uv::UniqueVector{T}) where {T} =
    get(uv.lookup, p.x, nothing)

function findfirst!(p::EqualTo{<:T}, uv::UniqueVector{T}) where {T}
    rv = get!(uv.lookup, p.x) do
        # NOTE: does not provide any exception safety guarantee
        push!(uv.items, p.x)
        return length(uv.items)
    end
    @assert length(uv.items) == length(uv.lookup)
    return rv
end

findfirst(p::EqualTo, uv::UniqueVector{T}) where {T} =
    findfirst(isequal(convert(T, p.x)), uv)

function findfirst(p::EqualTo, uv::SubArray{<:Any,1,<:AbstractUniqueVector})
    ip = findfirst(p, uv.parent)
    ip === nothing ? nothing : findfirst(isequal(ip), first(uv.indices))
end

findfirst!(p::EqualTo, uv::UniqueVector{T}) where {T} =
    findfirst!(isequal(convert(T, p.x)), uv)

findlast(p::EqualTo, uv::AbstractUniqueVector) =
    findfirst(p, uv)

indexin(a::AbstractArray, b::AbstractUniqueVector) =
    [findlast(isequal(elt), b) for elt in a]

function findall(p::Base.Fix2{typeof(in),<:AbstractUniqueVector},
                 a::AbstractArray)
    # The version in Base creates a Set that is unnecessary in our case, hence
    # the override here.
    [i for (i, ai) in (@static if VERSION >= v"1.1.0-DEV.832" pairs else enumerate end)(a) if p(ai)]
end

function findall(p::Base.Fix2{typeof(in),<:AbstractUniqueVector},
                 a::Tuple)
    # The version in Base creates a Set that is unnecessary in our case, hence
    # the override here.
    [i for (i, ai) in (@static if VERSION >= v"1.1.0-DEV.832" pairs else enumerate end)(a) if p(ai)]
end

function findnext(p::EqualTo, A::AbstractUniqueVector, i::Integer)
    idx = findfirst(p, A)
    idx >= i ? idx : nothing
end

function findprev(p::EqualTo, A::AbstractUniqueVector, i::Integer)
    idx = findfirst(p, A)
    idx <= i ? idx : nothing
end

function findall(p::EqualTo, uv::AbstractUniqueVector)
    idx = findfirst(p, uv)
    (idx == nothing) ? Int[] : Int[idx]
end

count(p::EqualTo, uv::AbstractUniqueVector) =
    Int(p.x ∈ uv)

allunique(::AbstractUniqueVector) = true

unique!(uv::AbstractUniqueVector)::AbstractUniqueVector = uv
unique(uv::UniqueVector)::Vector = copy(uv.items)

function push!(uv::UniqueVector{T}, item::T) where {T}
    if item in uv
        throw(ArgumentError("cannot add duplicate item to UniqueVector"))
    end
    # NOTE: does not provide any exception safety guarantee
    push!(uv.items, item)
    uv.lookup[item] = length(uv)
    @assert length(uv.items) == length(uv.lookup)
    return uv
end

push!(uv::UniqueVector{T}, item) where {T} =
    push!(uv, convert(T, item))

function pop!(uv::UniqueVector)
    if isempty(uv.items)
        throw(ArgumentError("array must be non-empty"))
    end
    # NOTE: does not provide any exception safety guarantee
    delete!(uv.lookup, uv.items[end])
    rv = pop!(uv.items)
    @assert length(uv.items) == length(uv.lookup)
    return rv
end

function setindex!(uv::UniqueVector{T}, item::T, idx::Integer) where {T}
    @boundscheck checkbounds(uv, idx)
    uv[idx] == item && return uv # nothing to do
    item ∉ uv || throw(ArgumentError("cannot set an element that exists elsewhere in UniqueVector"))
    # NOTE: does not provide any exception safety guarantee
    delete!(uv.lookup, uv.items[idx])
    uv.items[idx] = item
    uv.lookup[item] = idx
    @assert length(uv.items) == length(uv.lookup)
    return uv
end

setindex!(uv::UniqueVector{T}, item, idx::Integer) where {T} =
    setindex!(uv, convert(T, item), idx)

"`swap!(uv::UniqueVector, to::Int, from::Int)` interchange/swap the values on the indices `to` and `from` in the `UniqueVector`"
function swap!(uv::UniqueVector, to::Int, from::Int)
    if to == from
        @boundscheck checkbounds(uv,to)
        return uv
    end

    # NOTE: does not provide any exception safety guarantee
    previous_id  = uv[to]
    future_id    = uv[from]

    uv.items[to]   = future_id
    uv.items[from] = previous_id

    uv.lookup[previous_id] = from
    uv.lookup[future_id]   = to

    return uv
end

function permute!(uv::UniqueVector, perm::AbstractVector)
    ip = invperm(perm)
    map!(i -> ip[i], values(uv.lookup))
    permute!(uv.items, perm)
    return uv
end

function invpermute!(uv::UniqueVector, perm::AbstractVector)
    map!(i -> perm[i], values(uv.lookup))
    invpermute!(uv.items, perm)
    return uv
end

export AbstractUniqueVector, UniqueVector, findfirst!, swap!

end # module


text:
using CircularArrays
using OffsetArrays
using Test

@testset "index style" begin
    @test IndexStyle(CircularArray) == IndexCartesian()
    @test IndexStyle(CircularVector) == IndexLinear()
end

@testset "construction" begin
    @testset "construction ($T)" for T = (Float64, Int)
        data = rand(T,10)
        arrays = [CircularVector(data), CircularVector{T}(data),
                CircularArray(data), CircularArray{T}(data), CircularArray{T,1}(data)]
        @test all(a == first(arrays) for a in arrays)
        @test all(a isa CircularVector{T,Vector{T}} for a in arrays)
    end
end

@testset "type stability" begin
    @testset "type stability $(n)d" for n in 1:10
        a = CircularArray(fill(1, ntuple(_->1, n)))

        @test @inferred(a[1]) isa Int64
        @test @inferred(a[[1]]) isa CircularVector{Int64}
        @test @inferred(a[[1]']) isa CircularArray{Int64,2}
        @test @inferred(axes(a)) isa Tuple{Vararg{AbstractUnitRange}}
        @test @inferred(similar(a)) isa typeof(a)
        @test @inferred(similar(typeof(a), axes(a))) isa typeof(a)
        @test @inferred(a[a]) isa typeof(a)
    end
end

@testset "display" begin
    @testset "display $(n)d" for n in 1:3
        data = rand(Int64, ntuple(_->3, n))
        v1 = CircularArray(data)
        io = IOBuffer()
        io_compare = IOBuffer()

        print(io, v1)
        print(io_compare, data)
        @test String(take!(io)) == String(take!(io_compare))

        print(io, summary(v1))
        print(io_compare, summary(data))

        text = String(take!(io_compare))
        text = replace(text, " Vector" => " CircularVector")
        text = replace(text, " Matrix" => " CircularArray")
        text = replace(text, " Array" => (n == 1 ? " CircularVector" : " CircularArray"))
        text = replace(text, r"{.+}" => "(::$(string(typeof(data))))")
        @test String(take!(io)) == text
    end
end

@testset "vector" begin
    data = rand(Int64, 5)
    v1 = CircularVector(data)

    @test size(v1, 1) == 5
    @test parent(v1) == data
    @test typeof(v1) == CircularVector{Int64,Vector{Int64}}
    @test isa(v1, CircularVector)
    @test isa(v1, AbstractVector{Int})
    @test !isa(v1, AbstractVector{String})
    @test v1[2] == v1[2 + length(v1)]

    @test IndexStyle(v1) == IndexStyle(typeof(v1)) == IndexLinear()
    @test v1[0] == data[end]
    @test v1[-4:10] == [data; data; data]
    @test v1[-3:1][-1] == data[end]
    @test v1[[true,false,true,false,true]] == v1[[1,3,0]]
    @test all(e in data for e in v1)
    @test all(e in v1 for e in data)

    v1copy = copy(v1)
    v1_2 = v1[2]
    v1[2] = 0
    v1[3] = 0
    @test v1[2] == v1[3] == 0
    @test v1copy[2] == v1_2
    @test v1copy[7] == v1_2
    @test_throws MethodError v1[2] = "Hello"

    v2 = CircularVector("abcde", 5)

    @test prod(v2) == "abcde"^5

    @test_throws MethodError push!(v1, 15)

    @testset "deleteat!" begin
        @test deleteat!(CircularVector([1, 2, 3]), 5) == CircularVector([1, 3])
        @test deleteat!(CircularVector([1, 2, 3, 4]), 1:5:10) == CircularVector([3, 4])
        @test deleteat!(CircularVector([1, 2, 3, 4]), [1, 5]) == CircularVector([2, 3, 4])
        @test deleteat!(CircularVector([1, 2, 3, 4]), (1, 5)) == CircularVector([2, 3, 4])
        @test deleteat!(CircularVector([1, 2, 3, 4]), (1, 6)) == CircularVector([3, 4])
        @test deleteat!(CircularVector([1, 2, 3, 4]), (1, 3, 5)) == CircularVector([2, 4])
        @test deleteat!(CircularVector([1, 2, 3, 4]), (1, 5, 7)) == CircularVector([2, 4])
    end

    @testset "doubly circular" begin
        a = CircularVector([1, 2, 3, 4, 5])
        b = CircularVector(a)

        @test all(a[i] == b[i] for i in -50:50)
    end

    @testset "type stability" begin
        v3 = @inferred(map(x -> x+1, CircularArray([1, 2, 3, 4])))
        @test v3 isa CircularVector{Int64}
        @test v3 == CircularArray([2, 3, 4, 5])
        @test similar(v3, Base.OneTo(4)) isa typeof(v3)
        @test similar(typeof(v3), Base.OneTo(4)) isa typeof(v3)
        @test similar(typeof(v3), 4) isa typeof(v3)
        @test similar(typeof(v3), 4) isa typeof(v3)

        v4 = @inferred(CircularArray([1, 2, 3, 4]) .+ 1)
        @test v4 isa CircularVector{Int64}
        @test v4 == CircularArray([2, 3, 4, 5])

        v5 = v4 .> 3
        @test v5 isa CircularVector{Bool, BitVector}
        @test v5 == CircularArray([0, 0, 1, 1])
    end
end

@testset "matrix" begin
    b_arr = [2 4 6 8; 10 12 14 16; 18 20 22 24]
    a1 = CircularArray(b_arr)
    @test size(a1) == (3, 4)
    @test parent(a1) == b_arr

    @test a1[2, 3] == 14
    @test a1[2, Int32(3)] == 14
    a1[2, 3] = 17
    @test a1[2, 3] == 17
    @test a1[-1, 7] == 17
    @test a1[CartesianIndex(-1, 7)] == 17
    @test a1[-1:5, 4:10][1, 4] == 17
    @test a1[:, -1:-1][2, 1] == 17
    a1[CartesianIndex(-2, 7)] = 99
    @test a1[1, 3] == 99

    a1[18] = 9
    @test a1[18] == a1[-6] == a1[6] == a1[3,2] == a1[0,6] == b_arr[3,2] == b_arr[6] == 9

    @test IndexStyle(a1) == IndexStyle(typeof(a1)) == IndexCartesian()
    @test a1[3] == a1[3,1]
    @test a1[Int32(4)] == a1[1,2]
    @test a1[-1] == a1[length(a1)-1]

    @test a1[2, 3, 1] == 17 # trailing index
    @test a1[2, 3, 99] == 17
    @test a1[2, 3, :] == [17]

    @test !isa(a1, CircularVector)
    @test !isa(a1, AbstractVector)
    @test isa(a1, AbstractArray)

    @test size(reshape(a1, (2, 2, 3))) == (2, 2, 3)

    a2 = CircularArray(4, (2, 3))
    @test isa(a2, CircularArray{Int, 2})

    a3 = @inferred(a2 .+ 1)
    @test a3 isa CircularArray{Int64, 2}
    @test a3 == CircularArray(5, (2, 3))

    @testset "doubly circular" begin
        a = CircularArray(b_arr)
        da = CircularArray(a)

        @test all(a[i, j] == da[i, j] for i in -8:8, j in -8:8)
        @test all(a[i] == da[i] for i in -50:50)
    end
end

@testset "3-array" begin
    t3 = collect(reshape('a':'x', 2,3,4))
    c3 = CircularArray(t3)

    @test parent(c3) == t3

    @test c3[1,3,3] == c3[3,3,3] == c3[3,3,7] == c3[3,3,7,1]

    c3[3,3,7] = 'Z'
    @test t3[1,3,3] == 'Z'

    @test c3[3, CartesianIndex(3,7)] == 'Z'
    c3[Int32(3), CartesianIndex(3,7)] = 'ζ'
    @test t3[1,3,3] == 'ζ'

    c3[34] = 'J'
    @test c3[34] == c3[-38] == c3[10] == c3[2,2,2] == c3[4,5,6] == t3[2,2,2] == t3[10] == 'J'

    @test vec(c3[:, [CartesianIndex()], 1, 5]) == vec(t3[:, 1, 1])

    @test IndexStyle(c3) == IndexStyle(typeof(c3)) == IndexCartesian()
    @test c3[-1] == t3[length(t3)-1]

    @test_throws BoundsError c3[2,3] # too few indices
    @test_throws BoundsError c3[CartesianIndex(2,3)]

    @testset "doubly circular" begin
        c = CircularArray(t3)
        dc = CircularArray(c)

        @test all(c[i, j, k] == dc[i, j, k] for i in -5:5, j in -5:5, k in -5:5)
        @test all(c[i] == dc[i] for i in -50:50)
    end
end

@testset "offset indices" begin
    i = OffsetArray(1:5,-3)
    a = CircularArray(i)
    @test axes(a) == axes(i)
    @test a[1] == 4
    @test a[10] == a[-10] == a[0] == 3
    @test a[-2:7] == [1:5; 1:5]
    @test a[0:9] == [3:5; 1:5; 1:2]
    @test a[1:10][-10] == 3
    @test a[i] == OffsetArray([4,5,1,2,3],-3)

    @testset "type stability" begin
        @test @inferred(similar(a)) isa CircularVector
        @test @inferred(similar(typeof(a), axes(a))) isa CircularVector

        b = CircularVector([1, 2, 3, 4, 5])
        @test @inferred(similar(b, 3:5)) isa CircularVector
        @test @inferred(similar(typeof(b), 3:5)) isa CircularVector
    end

    circ_a = circshift(a,3)
    @test axes(circ_a) == axes(a)
    @test circ_a[1:5] == 1:5

    j = OffsetArray([true,false,true],1)
    @test a[j] == [5,2]

    data = reshape(1:9,3,3)
    a = CircularArray(OffsetArray(data,-1,-1))
    @test collect(a) == data
    @test all(a[x,y] == data[mod1(x+1,3),mod1(y+1,3)] for x=-10:10, y=-10:10)
    @test a[i,1] == CircularArray(OffsetArray([5,6,4,5,6],-2:2))
    @test a[CartesianIndex.(i,i)] == CircularArray(OffsetArray([5,9,1,5,9],-2:2))
    @test a[a .> 4] == 5:9
end


text:
module OptiMo

import LinearAlgebra.norm

export AbstractOptiModelMeta, OptiModelMeta, AbstractOptiModel

export obj, grad!, cons!, jprod!, jtprod!, proj!, prox!, objprox!

export grad, objgrad, objgrad!, cons, objcons, objcons!,
       jprod, jtprod, prox, objprox, proj, dist

const Maybe{T} = Union{T,Nothing}

include("optimo_utils.jl")
include("optimo_types.jl")

##########################################################
# Methods to be necessarily overridden in other packages.
##########################################################
"""
    fx = obj(prob, x)
"""
function obj end

"""
    grad!(prob, x, dfx)
"""
function grad! end

"""
    cons!(prob, x, cx)
"""
function cons! end

"""
    jprod!(prob, x, v, Jv)
"""
function jprod! end

"""
    jtprod!(prob, x, v, Jtv)
"""
function jtprod! end

"""
    proj!(prob, v, p)
"""
function proj! end

"""
    prox!(prob, x, a, z)
"""
function prox! end

"""
    gz = objprox!(prob, x, a, z)
"""
function objprox! end

##########################################################
# Methods with basic default implementation
##########################################################
"""
    dfx = grad(prob, x)
"""
function grad(prob::AbstractOptiModel, x::AbstractArray)
  @lencheck prob.meta.nvar x
  dfx = similar(x)
  grad!(prob, x, dfx)
  return dfx
end

"""
    fx, dfx = objgrad(prob, x)
"""
function objgrad(prob::AbstractOptiModel, x::AbstractArray)
  @lencheck prob.meta.nvar x
  fx = obj(prob, x)
  dfx = grad(prob, x)
  return fx, dfx
end

"""
    fx = objgrad!(prob, x, dfx)
"""
function objgrad!(prob::AbstractOptiModel, x::AbstractArray, dfx::AbstractArray)
  @lencheck prob.meta.nvar x dfx
  fx = obj(prob, x)
  grad!(prob, x, dfx)
  return fx
end

"""
    cx = cons(prob, x)
"""
function cons(prob::AbstractOptiModel, x::AbstractArray)
  @lencheck prob.meta.nvar x
  cx = similar(x, prob.meta.ncon)
  cons!(prob, x, cx)
  return cx
end

"""
    fx, cx = objcons(prob, x)
"""
function objcons(prob::AbstractOptiModel, x::AbstractArray)
  @lencheck prob.meta.nvar x
  fx = obj(prob, x)
  cx = prob.meta.ncon > 0 ? cons(prob, x) : eltype(x)[]
  return fx, cx
end

"""
    fx = objcons!(prob, x, cx)
"""
function objcons!(prob::AbstractOptiModel, x::AbstractArray, cx::AbstractArray)
  @lencheck prob.meta.nvar x
  @lencheck prob.meta.ncon cx
  fx = obj(prob, x)
  prob.meta.ncon > 0 && cons!(prob, x, cx)
  return fx
end

"""
    Jv = jprod(prob, x, v)
"""
function jprod(prob::AbstractOptiModel, x::AbstractArray, v::AbstractArray)
  @lencheck prob.meta.nvar x v
  Jv = similar(v, prob.meta.ncon)
  jprod!(prob, x, v, Jv)
  return Jv
end

"""
    Jtv = jtprod(prob, x, v)
"""
function jtprod(prob::AbstractOptiModel, x::AbstractArray, v::AbstractArray)
  @lencheck prob.meta.nvar x
  @lencheck prob.meta.ncon v
  Jtv = similar(x)
  jtprod!(prob, x, v, Jtv)
  return Jtv
end

"""
    z = prox(prob, x, a)
"""
function prox(prob::AbstractOptiModel, x::AbstractArray, a::Real)
  @lencheck prob.meta.nvar x
  z = similar(x)
  prox!(prob, x, a, z)
  return z
end

"""
    gz, z = objprox(prob, x, a)
"""
function objprox(prob::AbstractOptiModel, x::AbstractArray, a::Real)
  @lencheck prob.meta.nvar x
  z = similar(x)
  gz = objprox!(prob, x, a, z)
  return gz, z
end

"""
    p = proj(prob, v)
"""
function proj(prob::AbstractOptiModel, v::AbstractArray)
  @lencheck prob.meta.ncon v
  p = similar(v)
  proj!(prob, v, p)
  return p
end

"""
    d = dist(prob, v)
"""
function dist(prob::AbstractOptiModel, v::AbstractArray)
  @lencheck prob.meta.ncon v
  p = proj(prob, v)
  return norm(p - v)
end

##########################################################
# Model subtypes
##########################################################
include("optimo_autodiff.jl")
include("optimo_slack.jl")
include("optimo_auglag.jl")
include("optimo_prox.jl")
include("optimo_feas.jl")
include("optimo_nlp.jl")

##########################################################
# Output types
##########################################################
include("optimo_output.jl")

end # module


text:
using Test

entry_dimesion = 2
number_of_hidden_units = 3
output_dimension = 2
OLNN = RandomWeightsNN(
    entry_dimesion,
    number_of_hidden_units,
    output_dimension
)

@testset "Dimension of one layer networks random initialization" begin
    # Weights have correct dimensions
    # Notemos que OLNN ha sido creada con la iniciacilización aleatoria,
    # Las única hipótesis que debe de cumplir es que:
    # 1. Inicialización con las dimensiones correctas
    @test size(OLNN.W1)==(number_of_hidden_units, 1+entry_dimesion)
    @test size(OLNN.W2)==(output_dimension, number_of_hidden_units)
    # 2. Por la aleatoriedad generada no todas las entradas debieran de ser iguales
    @test OLNN.W1[1,:] != OLNN.W1[2,:]
    @test OLNN.W2[1,:] != OLNN.W2[2,:]
end

@testset "One layer created from matrix" begin
    S = [1,2] #vector
    A = [3 4; 4 6] # matrix
    B = reshape([ 1 ; 1 ],1,2) # matrix 2 x 1
    h = FromMatrixNN(S, A, B)
    # Comprobación de tipo correcto
    @test typeof(h) <: AbstractOneLayerNeuralNetwork
    # Comprobación de tamaños correctos
    ### Para la matriz W_1
    (n_rows1, n_columns1) = size(h.W1)
    (n_rows2, n_columns2) = size(h.W2)
    (r_a, c_a) = size(A)
    @test n_rows1 == r_a
    @test n_columns1 == c_a+1
     ### Para la matriz W_1
     (n_rows2, n_columns2) = size(h.W2)
     (r_b, c_b) = size(B)
     @test n_rows2 == r_b
     @test n_columns2 == r_a
     @test n_columns2 == c_b
     println("Revisión ocular:")
     println(h)
end


text:
@snn_kw struct RateSynapseParameter{FT=Float32}
    lr::FT = 1e-3
end

@snn_kw mutable struct RateSynapse{VIT=Vector{Int32},VFT=Vector{Float32}}
    param::RateSynapseParameter = RateSynapseParameter()
    colptr::VIT # column pointer of sparse W
    I::VIT      # postsynaptic index of W
    W::VFT  # synaptic weight
    rI::VFT # postsynaptic rate
    rJ::VFT # presynaptic rate
    g::VFT  # postsynaptic conductance
    records::Dict = Dict()
end

"""
[Rate Synapse](https://brian2.readthedocs.io/en/2.0b4/resources/tutorials/2-intro-to-brian-synapses.html)
"""
RateSynapse

function RateSynapse(pre, post; σ = 0.0, p = 0.0, kwargs...)
    w = σ / √(p * pre.N) * sprandn(post.N, pre.N, p)
    rowptr, colptr, I, J, index, W = dsparse(w)
    rI, rJ, g = post.r, pre.r, post.g
    RateSynapse(;@symdict(colptr, I, W, rI, rJ, g)..., kwargs...)
end

function forward!(c::RateSynapse, param::RateSynapseParameter)
    @unpack colptr, I, W, rI, rJ, g = c
    @unpack lr = param
    fill!(g, zero(eltype(g)))
    @inbounds for j in 1:(length(colptr) - 1)
        rJj = rJ[j]
        for s = colptr[j]:(colptr[j+1] - 1)
            g[I[s]] += W[s] * rJj
        end
    end
end

function plasticity!(c::RateSynapse, param::RateSynapseParameter, dt::Float32, t::Float32)
    @unpack colptr, I, W, rI, rJ, g = c
    @unpack lr = param
    @inbounds for j in 1:(length(colptr) - 1)
        s_row = colptr[j]:(colptr[j+1] - 1)
        rIW = zero(Float32)
        for s in s_row
            rIW += rI[I[s]] * W[s]
        end
        Δ = lr * (rJ[j] - rIW)
        for s in s_row
            W[s] += rI[I[s]] * Δ
        end
    end
end


text:
using DataFrames
using CSV

#Exponents
a = 9^2

A = [10, 20, 30]

#Julia is 1-indexed
A[1]

#MATRICES

#Hacer una matriz 3x3
M = [1 2 3;4 5 6;7 8 9]

#In here, I access the position [2,2] of the matrix

M[2,2] = 28
M

#Here I traspose a matrix
M'

#And here I calculate the inverse of the matrix
inv(M)

#Dictionaries
D = Dict("name" => "Sofia Alvarez", "birthyear" => 2000, "university" => "Los Andes")

D["name"]

D.count

#Strings ARE immutable
nombre = "Sofia Alvarez"
println(nombre[1])
nombre[1:5]

#Loops
var = 0
stop = 3
for i=1:stop
    var = i + pi
    print(i, " ")
end
println("The value of the sum is: ", var)

variable = 28
    if variable == 28
        println("Tau day!")
    else
        println("Ups!")
end

personas = CSV.read("MOCK_DATA.csv")

size(personas)

names(personas)

#Summary stats of the column age
describe(personas[:age])

max_stars_repo_path: Learning Julia.jl

max_stars_repo_name: sofiaalvarezlopez/Learning-Julia

max_issues_repo_path: Learning Julia.jl

max_issues_repo_name: sofiaalvarezlopez/Learning-Julia

max_forks_repo_path: Learning Julia.jl

max_forks_repo_name: sofiaalvarezlopez/Learning-Julia


text:
# This file was generated, do not modify it. # hide
lr1 = machine(LinMod, Xhp, y) # poly of degree 1 (line)
fit!(lr1, rows=train)

LinMod.feature_selector.features = [:hp1, :hp2] # poly of degree 2
lr2 = machine(LinMod, Xhp, y)
fit!(lr2, rows=train)

LinMod.feature_selector.features = [:hp1, :hp2, :hp3] # poly of degree 3
lr3 = machine(LinMod, Xhp, y)
fit!(lr3, rows=train)


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

@testset "Scalarize" begin
    bridged_mock = MOIB.Constraint.Scalarize{Float64}(mock)

    MOIT.basic_constraint_tests(
        bridged_mock,
        config,
        include = [
            (F, S) for F in [
                MOI.VectorOfVariables,
                MOI.VectorAffineFunction{Float64},
                MOI.VectorQuadraticFunction{Float64},
            ] for S in [MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros]
        ],
    )

    # VectorOfVariables-in-Nonnegatives
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-3, -1],
        )
    MOIT.lin1vtest(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ),
    )
    test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.SingleVariable, MOI.GreaterThan{Float64}, 0),
        ),
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        ),
    )
    func = MOI.get(bridged_mock, MOI.ConstraintFunction(), ci)
    MOI.delete(bridged_mock, func.variables[2])
    new_func = MOI.VectorOfVariables(func.variables[[1, 3]])
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) == new_func
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.Nonnegatives(2)

    @testset "$attr" for attr in [
        MOI.ConstraintPrimalStart(),
        MOI.ConstraintDualStart(),
    ]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        values = [1.0, 2.0]
        MOI.set(bridged_mock, attr, ci, values)
        @test MOI.get(bridged_mock, attr, ci) == values
    end

    test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.SingleVariable, MOI.GreaterThan{Float64}, 0),
        ),
    )

    # VectorAffineFunction-in-Nonnegatives
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0, 2, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-3, -1],
        )
    MOIT.lin1ftest(bridged_mock, config)

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ),
    )

    @testset "$attr" for attr in [
        MOI.ConstraintPrimalStart(),
        MOI.ConstraintDualStart(),
    ]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        values = [1.0, -2.0]
        MOI.set(bridged_mock, attr, ci, values)
        @test MOI.get(bridged_mock, attr, ci) == values
    end

    test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        ),
    )

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonnegatives,
            }(),
        ),
    )

    @testset "$attr" for attr in [
        MOI.ConstraintPrimalStart(),
        MOI.ConstraintDualStart(),
    ]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        values = [1.0, -2.0, 3.0]
        MOI.set(bridged_mock, attr, ci, values)
        @test MOI.get(bridged_mock, attr, ci) == values
    end

    test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        ),
    )

    # Test setting VectorAffineFunction with nonzero constants
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.0, 0.0]),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.75]),
    )
    MOIT.solve_func_vectoraffine_nonneg(bridged_mock, config)

    # VectorOfVariables-in-Nonnegatives
    # VectorOfVariables-in-Nonpositives
    # VectorOfVariables-in-Zeros
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [-4, -3, 16, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [7, 2, -4],
        )
    MOIT.lin2vtest(bridged_mock, config)

    # VectorAffineFunction-in-Nonnegatives
    # VectorAffineFunction-in-Nonpositives
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [-4, -3, 16, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [7, 2, -4, 7],
        )
    MOIT.lin2ftest(bridged_mock, config)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Vec{Dim,T}

A vector in `Dim`-dimensional space with coordinates of type `T`.

A vector can be obtained by subtracting two [`Point`](@ref) objects.

## Examples

```julia
A = Point(0.0, 0.0)
B = Point(1.0, 0.0)
v = B - A
```

### Notes

- A `Vec` is a `SVector` from StaticArrays.jl
- Type aliases are `Vec1`, `Vec2`, `Vec3`, `Vec1f`, `Vec2f`, `Vec3f`
"""
const Vec = SVector

# type aliases for convenience
const Vec1  = Vec{1,Float64}
const Vec2  = Vec{2,Float64}
const Vec3  = Vec{3,Float64}
const Vec1f = Vec{1,Float32}
const Vec2f = Vec{2,Float32}
const Vec3f = Vec{3,Float32}