text:
using SimilaritySearch
using SimilaritySearch.Graph
using Test
#
# This file contains a set of tests for SearchGraph over databases of vectors (of Float32)
#

function test_index_search_with_hint(dist::Function, ksearch::Int, search_algo, neighborhood_algo)
    @testset "indexing with different algorithms" begin
        index = fit(SearchGraph, dist, Vector{Float32}[], search_algo=search_algo, neighborhood_algo=neighborhood_algo)

        n = 10_000
        dim = 3
        @info "inserting items to the index"
        for i in 1:n
            vec = rand(Float32, dim)
            push!(index, dist, vec)
        end

        @info "done; now testing with hint"

        m = 1

        start = time()
        for i in 1:m
            q = index.db[i]
            #search(index, dist, q, KnnResult(ksearch), hints=Int.(index.links[i]))
            search(index, dist, q, KnnResult(ksearch), hints=Int.(index.links[i]))
        end

        @info "Hints=true, noise=false: Query time $((time() - start) / m)"

        start = time()
        for i in 1:m
            q = index.db[i] + rand(Float32, dim) .* Float32(0.0001)
            search(index, dist, q, KnnResult(ksearch), hints=Int.(index.links[i]))
        end

        @info "Hints=true, noise=true: Query time $((time() - start) / m)"

        start = time()
        for i in 1:m
            q = index.db[i]
            search(index, dist, q, KnnResult(ksearch))
        end
        @info "Hints=false, noise=false: Query time $((time() - start) / m)"

        start = time()
        for i in 1:m
            q = index.db[i] + rand(Float32, dim) .* Float32(0.0001)
            search(index, dist, q, KnnResult(ksearch))
        end
        @info "Hints=false, noise=true: Query time $((time() - start) / m)"

    end
end

function test_index(dist::Function, ksearch::Int, search_algo, neighborhood_algo)
    @testset "indexing with different algorithms" begin
        index = fit(SearchGraph, dist, Vector{Float32}[], recall=0.9, search_algo=search_algo, neighborhood_algo=neighborhood_algo)
        n = 10_000
        dim = 3

        @info "inserting items to the index"
        for i in 1:n
            vec = rand(Float32, dim)
            push!(index, dist, vec)
        end

        @info "done; now testing"
        @test length(index.db) == n
        res = search(index, dist, rand(Float32, dim), KnnResult(ksearch))
        @show res
        return index, length(res)
    end
end

@testset "some vector indexing" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are requiered
    ksearch = 10
    acc = 0
    expected_acc = 0

    for search_algo in [TIHCSearch(), IHCSearch(), BeamSearch()]
        for neighborhood_algo in [SatNeighborhood()]
        #for neighborhood_algo in [EssencialNeighborhood(), FixedNeighborhood(8), GallopingNeighborhood(), GallopingSatNeighborhood(), LogNeighborhood(), LogSatNeighborhood(), SatNeighborhood(), VorNeighborhood()]
            # for dist in Any[l2_distance, L2Distance(), L1Distance(), LInfDistance(), LpDistance(0.5)]
            dist = l2_distance
            index, numres = test_index(dist, ksearch, search_algo, neighborhood_algo)
            acc += numres
            expected_acc += ksearch
        end
    end

    # this is not really an error, but we test it anyway, it is more about the quality of the results
    # @test acc / expected_acc > 0.9

    test_index_search_with_hint(l2_distance, ksearch, BeamSearch(), SatNeighborhood())
    # @show "Showing AKNN ($k)"
    ## n = length(index.db)
    ## k = 3
    ## aknn = compute_aknn(index, l2_distance, k)
    ## @test n == length(aknn)
    ## for p in aknn
    ##     @test length(p) > 0
    ## end
end


text:
#***********************************************************************
#  Copyright 2006 John A. Trangenstein
#
#  This software is made available for research and instructional use
#  only.
#  You may copy and use this software without charge for these
#  non-commercial purposes, provided that the copyright notice and
#  associated text is reproduced on all copies.
#  For all other uses (including distribution of modified versions),
#  please contact the author at
#    John A. Trangenstein
#    Department of Mathematics
#    Duke University
#    Durham, NC 27708-0320
#    USA
#  or
#    johnt@math.duke.edu
#
#  This software is made available "as is" without any assurance that it
#  is completely correct, or that it will work for your purposes.
#  Use the software at your own risk.
#***********************************************************************

using OffsetArrays

function fluxderv(fc::Int,lc::Int,ifirst::Int,ilast::Int,
                  conservd::OffsetArray{Float64,1}, dfdu::OffsetArray{Float64,1})
#   ******************************************************************
#   compute derivative of flux with respect to conserved variable
#   ******************************************************************
    @inbounds for ic=ifirst:ilast
        dfdu[ic]=velocity
    end
end

function riemann(fs::Int,ls::Int,ifirst::Int,ilast::Int,
                 left::OffsetArray{Float64,1},right::OffsetArray{Float64,1},
                 flux::OffsetArray{Float64,1})
#   ******************************************************************
#   the following riemann solver is hard-wired for linear advection
#   with velocity > 0
#   ******************************************************************
    @inbounds for ie=ifirst:ilast+1
        flux[ie]=velocity*left[ie]
    end
end

ext: jl


text:
Base.@kwdef mutable struct PortUsed <:Leaf
    id::Int         = 0
    proto::String   = ""
    state::String   = ""
end
StructTypes.names(::Type{PortUsed}) = ((:id, :portid),)

Base.@kwdef mutable struct OSFingerprint <:Leaf
    fingerprint::String = ""
end
fingerprint(osf::OSFingerprint) = osf.fingerprint

Base.@kwdef mutable struct OSClass <:Marsh
    accuracy::Int           = 0
    family::String          = ""
    gen::String             = ""
    type::String            = ""
    vendor::String          = ""
    cpes::Vector{String}    = Vector{String}()
end
StructTypes.names(::Type{OSClass}) = ((:cpes, :cpe), (:family, :osfamily), (:gen, :osgen))

accuracy(class::OSClass)    = class.accuracy
family(class::OSClass)      = class.family
gen(class::OSClass)         = class.gen
vendor(class::OSClass)      = class.vendor
type(class::OSClass)        = class.type

Base.@kwdef mutable struct OSMatch <:Marsh
    accuracy::Int               = 0
    line::Int                   = 0
    name::String                = ""
    classes::Vector{OSClass}    = Vector{OSClass}()
end
StructTypes.names(::Type{OSMatch}) = ((:classes, :osclass),)

accuracy(match::OSMatch)    = match.accuracy
line(match::OSMatch)        = match.line
name(match::OSMatch)        = match.name
classes(match::OSMatch)     = match.classes
Base.string(match::OSMatch) = name(match)

Base.@kwdef mutable struct OS <:Marsh
    ports::Vector{PortUsed}              = Vector{PortUsed}()
    matches::Vector{OSMatch}             = Vector{OSMatch}()
    fingerprints::Vector{OSFingerprint}  = Vector{OSFingerprint}()
end
StructTypes.names(::Type{OS}) = (
    (:ports,        :portused),
    (:matches,      :osmatch),
    (:fingerprints, :osfingerprint))

ports(os::OS)       = os.ports
matches(os::OS)     = os.matches
fingerprint(os::OS) = os.fingerprints


text:
include("../src/NeuralNetwork.jl")
include("../src/MNIST/src/MNIST.jl")

using .NeuralNetwork
using .MNIST

trd = traindata()
tsd = testdata()

net = read(joinpath(@__DIR__, "net1"), Chain)

test(net, tsd)


text:
@testset "gridgen" begin
rt = R_TARGET_DEFAULT  # 1.9
rmax = R_MAX_DEFAULT  # 2.0

@testset "gen_sublprim1d stability" begin
    domain = OpenInterval((-100,100),10)
    intvprim = [ClosedInterval((-100,-20)), ClosedInterval((20,100)), ClosedInterval((-30,30),1)]
    intvdual = ClosedInterval[]

    lprim0 = Float64[-75, 75]
    ldual0 = Float64[-10,0,10]

    sublprim = MaxwellFDM.gen_sublprim1d(domain, DUAL, intvprim, intvdual, lprim0, ldual0)
    @test @inferred(MaxwellFDM.gen_sublprim1d(domain, DUAL, intvprim, intvdual, lprim0, ldual0)) == sublprim
end  # @testset "gen_sublprim1d stability"

@testset "gen_lprim1d for fixed intervals to fail" begin
    domain = OpenInterval((-3.,1.), Inf)
    lprim0 = Float64[]
    ldual0 = Float64[]
    intvdual = ClosedInterval[]

    intvprim = [ClosedInterval((-0.2,0.5), 0.05), ClosedInterval((-0.3,1.2), 0.15)]
    sublprim = MaxwellFDM.gen_sublprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0)
    # @info "sublprim = $sublprim"
    @test_throws ArgumentError MaxwellFDM.gen_lprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0, rt, rmax)

    intvprim = [ClosedInterval((-3.,-0.1), 1.0), ClosedInterval((0.1,1.), 1.0)]
    sublprim = MaxwellFDM.gen_sublprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0)
    @test_throws ErrorException MaxwellFDM.gen_lprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0, rt, rmax)
end  # @testset "gen_lprim1d for fixed intervals to fail"

@testset "gen_lprim1d for fixed intervals" begin
    domain = OpenInterval((-100,100),10)
    intvprim = [ClosedInterval((-100,-20)), ClosedInterval((20,100)), ClosedInterval((-30,30),1)]
    intvdual = ClosedInterval[]

    lprim0 = Float64[-75, 75]
    ldual0 = Float64[-10,0,10]

    sublprim = MaxwellFDM.gen_sublprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0)
    # @info "sublprim = $sublprim"
    @test typeof(sublprim) == Vector{Vector{Float64}}
    @test isodd(length(sublprim))
    subgrid = sublprim[1:2:end]
    ∆lt = sublprim[2:2:end]
    @test length.(∆lt) == ones(Int, length(∆lt))
    # @info "domain = $domain"
    # @info "intvprim = $intvprim"
    # @info "lprim0 = $lprim0"
    # @info "ldual0 = $ldual0"
    # @info "rt = $rt"
    # @info "rmax = $rmax"
    lprim = MaxwellFDM.gen_lprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0, rt, rmax)
    # @info "lprim = $lprim"
    @test all(lprim[[1,end]] .== bounds(domain))
    @test MaxwellFDM.issmooth(diff(lprim), rmax)
end  # @testset "gen_lprim1d for fixed intervals"

@testset "gen_lprim1d for fixed intervals, dual domain" begin
    domain = OpenInterval((-100,100),10)
    intvprim = [ClosedInterval((-100,-20)), ClosedInterval((20,100)), ClosedInterval((-30,30),1)]
    intvdual = ClosedInterval[]

    lprim0 = Float64[-75, 75]
    ldual0 = Float64[-10,0,10]

    sublprim = MaxwellFDM.gen_sublprim1d(domain, DUAL, intvprim, intvdual, lprim0, ldual0)
    # @info "sublprim = $sublprim"
    @test typeof(sublprim) == Vector{Vector{Float64}}
    @test isodd(length(sublprim))
    subgrid = sublprim[1:2:end]
    ∆lt = sublprim[2:2:end]
    @test length.(∆lt) == ones(Int, length(∆lt))
    # @info "domain = $domain"
    # @info "intvprim = $intvprim"
    # @info "lprim0 = $lprim0"
    # @info "ldual0 = $ldual0"
    # @info "rt = $rt"
    # @info "rmax = $rmax"
    lprim = MaxwellFDM.gen_lprim1d(domain, DUAL, intvprim, intvdual, lprim0, ldual0, rt, rmax)
    # @info "lprim = $lprim"
    @test mean(lprim[1:2]) == bounds(domain)[1]
    @test mean(lprim[end-1:end]) == bounds(domain)[2]
    @test MaxwellFDM.issmooth(diff(lprim), rmax)
end  # @testset "gen_lprim1d for fixed intervals, dual domain"

@testset "gen_lprim1d for randomly generated intervals" begin  # case generated by @testset "gen_lprim1d for random intervals"
    domain = OpenInterval((-2.4275370034248813,1.5578784568403026),0.0038696573133163836)
    intvprim = [ClosedInterval((-2.42754,-0.105743),0.00429962),
                ClosedInterval((-0.433823,1.49839),0.00357817),
                ClosedInterval((0.751153,1.55788),0.00149394),
                ClosedInterval((-0.409778,-0.21345),0.00036357)]
    intvdual = ClosedInterval[]

    lprim0 = [-0.314801,-0.314801,-0.770512,-0.685944,-0.664589]
    ldual0 = [-0.395004,-0.225236,0.081886,-0.477425]

    sublprim = MaxwellFDM.gen_sublprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0)
    # @info "sublprim = $sublprim"
    @test typeof(sublprim) == Vector{Vector{Float64}}
    @test isodd(length(sublprim))
    subgrid = sublprim[1:2:end]
    ∆lt = sublprim[2:2:end]
    @test length.(∆lt) == ones(Int, length(∆lt))
    lprim = MaxwellFDM.gen_lprim1d(domain, PRIM, intvprim, intvdual, lprim0, ldual0, rt, rmax)
    # @info "lprim = $lprim"
    @test all(lprim[[1,end]] .== bounds(domain))
    @test MaxwellFDM.issmooth(diff(lprim), rmax)
end  # @testset "gen_lprim1d for randomly generated intervals"

@testset "gen_lprim1d for randomly generated intervals, dual domain" begin  # case generated by @testset "gen_lprim1d for random intervals, dual domain"
    domain = OpenInterval((-2.4275370034248813,1.5578784568403026),0.0038696573133163836)
    intvprim = [ClosedInterval((-2.42754,-0.105743),0.00429962),
                ClosedInterval((-0.433823,1.49839),0.00357817),
                ClosedInterval((0.751153,1.55788),0.00149394),
                ClosedInterval((-0.409778,-0.21345),0.00036357)]
    intvdual = ClosedInterval[]

    lprim0 = [-0.314801,-0.314801,-0.770512,-0.685944,-0.664589]
    ldual0 = [-0.395004,-0.225236,0.081886,-0.477425]

    sublprim = MaxwellFDM.gen_sublprim1d(domain, DUAL, intvprim, intvdual, lprim0, ldual0)
    # @info "sublprim = $sublprim"
    @test typeof(sublprim) == Vector{Vector{Float64}}
    @test isodd(length(sublprim))
    subgrid = sublprim[1:2:end]
    ∆lt = sublprim[2:2:end]
    @test length.(∆lt) == ones(Int, length(∆lt))
    lprim = MaxwellFDM.gen_lprim1d(domain, DUAL, intvprim, intvdual, lprim0, ldual0, rt, rmax)
    # @info "lprim = $lprim"

    # Entries of diff(wprim) can exceed max∆l(domain)[w], but not much.  See gridgen.jl:fill_constant
    @test all(diff(lprim) .< max∆l(domain)*rmax)
    @test mean(lprim[1:2]) == bounds(domain)[1]
    @test mean(lprim[end-1:end]) == bounds(domain)[2]
    @test MaxwellFDM.issmooth(diff(lprim), rmax)
end  # @testset "gen_lprim1d for randomly generated intervals, dual domain"

@testset "comp_lprim1d" begin
    sublprim = [[0., 0.5], [1.], [10., 11., 12], [1.5], [20.5, 21.], [2.], [30., 31.]];
    lprim = MaxwellFDM.comp_lprim1d(sublprim, rt, rmax)
    @test MaxwellFDM.issmooth(diff(lprim), rmax)
    lprim_result = [0.0,0.5,1.20711,2.20711,3.32038,4.43365,5.54692,6.66019,7.77346,8.88673,10.0,11.0,12.0,13.6982,15.3964,17.0947,18.7929,19.7929,20.5,21.0,21.7071,22.7071,23.7071,24.909,26.3536,27.7981,29.0,30.0,31.0]
    @test isapprox(lprim, lprim_result, rtol = 1e-6)
end  # @testset "comp_lprim1d"

# Test with shapes
@testset "gen_lprim" begin
    ∆L = 5
    ∆l = 1

    vac = Material("vacuum")
    domain = Object(Box([0,0,0], [200,200,200]), vac, 10)

    h = 50
    r = 20
    sprim = [Object(Cylinder([0,0,0],r,h,[0,0,1]), vac, [∆l,∆l,∆L]), Object(Sphere([0,0,h/2], r), vac, ∆l), Object(Sphere([0,0,-h/2], r), vac, ∆l)]
    # sprim = [Object(Sphere([0,0,0],r), vac, [∆l,∆l,∆L]), Object(Sphere([0,0,h/2], r), vac, ∆l), Object(Sphere([0,0,-h/2], r), vac, ∆l)]

    (xprim, yprim, zprim) = gen_lprim(domain, (PRIM,PRIM,PRIM), sprim)
    # @info "xprim = $xprim"
    # @info "yprim = $yprim"
    # @info "zprim = $zprim"
    @test xprim == yprim

    # Entries of diff(wprim) can exceed max∆l(domain)[w], but not much.  See gridgen.jl:fill_constant
    @test all(diff(xprim) .< max∆l(domain)[1]*rmax)
    @test all(diff(yprim) .< max∆l(domain)[2]*rmax)
    @test all(diff(zprim) .< max∆l(domain)[3]*rmax)

    @test MaxwellFDM.issmooth(diff(xprim), rmax)
    @test MaxwellFDM.issmooth(diff(yprim), rmax)
    @test MaxwellFDM.issmooth(diff(zprim), rmax)
end

end  # @testset "gridgen"


text:
#Correr este archivo para correr todos los tests
#Desde el directorio principal correr "julia --project test/runtests.jl"
using Test

using Reporstat, DataFrames, CSV

#agregar aquí el nombre del archivo donde hagan los test de su módulo
tests = [
"utils.jl",
"select.jl"
]

for test in tests
  @testset "$test" begin
    include(test)
  end
end

max_stars_repo_name: yollotltamayo/Covid

max_issues_repo_name: yollotltamayo/Covid

max_forks_repo_name: yollotltamayo/Covid


text:
module Jacobi

# package code goes here
export jacobi, djacobi, jacobi_zeros, jacobi_zeros!
export zgj, zglj, zgrjm, zgrjp, wgj, wglj, wgrjm, wgrjp
export dgj, dglj, dgrjm, dgrjp
export lagrange, lagrange!, interp_mat
export chebyshev, chebyshev2, dchebyshev, dchebyshev2
export chebyshev_zeros!, chebyshev_zeros, chebyshev2_zeros!, chebyshev2_zeros
export legendre, dlegendre, legendre_zeros
export chebyshev!, chebyshev2!, dchebyshev!, dchebyshev2!
export legendre!, jacobi!, djacobi!
export poly_chebyshev, poly_chebyshev2, poly_dchebyshev, poly_dchebyshev2
export poly_legendre, poly_dlegendre, poly_jacobi, poly_djacobi

export QUADRATURE_TYPE, GJ, GLJ, GRJM, GRJP
export qzeros, qweights, qdiff, Quadrature
export qalpha, qbeta

using SpecialFunctions: gamma

include("jac_poly.jl")
include("chebyshev.jl")
include("legendre.jl")
include("gauss_quad.jl")
include("polynomials.jl")

end # module


text:
function constraint_frequency_stab_OPF(pm::_PM.AbstractPowerModel) #gurobi
    if pm.setting["Permanentloss"] == true
        constraint_frequency_stab_OPF_MP_PL(pm)
    elseif pm.setting["FSprotection"] == true || pm.setting["NSprotection"] == true
        constraint_frequency_stab_OPF_MP_FSNS(pm)
    end
end

function constraint_frequency_stab_OPF_MP_PL(pm::_PM.AbstractPowerModel)
    base_nws = pm.setting["base_list"]
    cont_nws = pm.setting["Cont_list"]
    ev_syncarea = pm.setting["syncarea"]

    for (base, nw, br) in cont_nws
        load = _PM.ref(pm, nw, :load, 1) # get index later on
        syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
        Zb1 = _PM.var(pm, nw)[:zb1]

        # for i in syncarea # when there is more than one synchornous areas
        i = ev_syncarea
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
        JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
        JuMP.@constraint(pm.model,  Phvdccaux[i]  == 0)
        JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] )
        JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])

        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        H = reserves[bi_bp[i]]["H"]
        Tcl = reserves[bi_bp[i]]["Tcl"]
        M = 1 # largest possible Pl=0.15 for now

        z1 = _PM.var(pm, nw, :z1, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))

        ############## frequency minimum in interval 1 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

         k11 = _PM.var(pm, nw, :k11, i)
         k12 = _PM.var(pm, nw, :k12, i)
         k11_dup = _PM.var(pm, nw, :k11_dup, i)
         k12_dup = _PM.var(pm, nw, :k12_dup, i)
         variable_on_off_switch(pm, k11, k11_dup, M, z1)
         variable_on_off_switch(pm, k12, k12_dup, M, z1)

         JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
         JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
         JuMP.@constraint(pm.model,k12_dup^2 <=  k11_dup)

         ############## frequency minimum in interval 2 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  M*(1-z41))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg   - M*(z41))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tf)   <=   Phvdcoaux[i]  + M*(1-z42))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tf)   >= e +  Phvdcoaux[i]  - M*(z42))
         JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
         JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

         k41 = _PM.var(pm, nw, :k41, i)
         k42 = _PM.var(pm, nw, :k42, i)
         k43 = _PM.var(pm, nw, :k43, i)
         k41_dup = _PM.var(pm, nw, :k41_dup, i)
         k42_dup = _PM.var(pm, nw, :k42_dup, i)
         k43_dup = _PM.var(pm, nw, :k43_dup, i)
         variable_on_off_switch(pm, k41, k41_dup, M, z4)
         variable_on_off_switch(pm, k42, k42_dup, M, z4)
         variable_on_off_switch(pm, k43, k43_dup, M, z4)

         JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4) )
         JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
         JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf )/2 )
         JuMP.@constraint(pm.model,[k41_dup/sqrt(2), k42_dup/sqrt(2), k43_dup, 0] in JuMP.RotatedSecondOrderCone() )

         JuMP.@constraint(pm.model, z1 + z4 == 1)
     end
end

function constraint_frequency_stab_OPF_MP_FSNS(pm::_PM.AbstractPowerModel)
    base_nws = pm.setting["base_list"]
    cont_nws = pm.setting["Cont_list"]
    ev_syncarea = pm.setting["syncarea"]

    for (base, nw, br) in cont_nws
        load = _PM.ref(pm, nw, :load, 1)
        syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
        Zb1 = _PM.var(pm, nw)[:zb1]

        # for i in syncarea
        i = ev_syncarea
        conv_conn = _PM.ref(pm, nw, :bus_arcs_conv, i)

        if pm.setting["FSprotection"] == true
            # display("FSprotection")
            Mmax = 72/ load["pd"]
            # for k in conv_conn
            #     display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= - (Pconv1[k]) /  load["pd"]) )
            #     display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  <= - (Pconv1[k]) / load["pd"] + Mmax*(1-Zb1[k]) ) )
            #     display(JuMP.@constraint(pm.model,  Phvdccaux[i]  >= - Pconv2[k] /  load["pd"] -  Mmax*(1-Zb1[k]) ) )
            #     display(JuMP.@constraint(pm.model,  Phvdccaux[i]  <= - Pconv2[k] / load["pd"] +   Mmax*(1-Zb1[k]) ) )
            # end
            # JuMP.@constraint(pm.model,  sum(Zb1[k] for k in conv_conn)==1 )

            display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_convs_dirct, i)) / load["pd"]) )
            display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_convs_dirct, i)) /load["pd"]) )
            # display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"] + sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"] ) )
            JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] )
            JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])

         elseif pm.setting["NSprotection"] == true
            # display("NSprotection")
            # display(load["pd"])
            display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"]))
            display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"]))
            JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4])
            JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])

        end

        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])

        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]

        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        Tcl = reserves[bi_bp[i]]["Tcl"]
        H = reserves[bi_bp[i]]["H"]

        z1 = _PM.var(pm, nw, :z1, i)
        z2 = _PM.var(pm, nw, :z2, i)
        z3 = _PM.var(pm, nw, :z3, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z21 = _PM.var(pm, nw, :z21, i)
        z31 = _PM.var(pm, nw, :z31, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z22 = _PM.var(pm, nw, :z22, i)
        z32 = _PM.var(pm, nw, :z32, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))
        M = 1
        # M = ub(Pf + Pg +  Phvdccaux[i])
       ############## frequency minimum in interval 1 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

         k11 = _PM.var(pm, nw, :k11, i)
         k12 = _PM.var(pm, nw, :k12, i)
         k11_dup = _PM.var(pm, nw, :k11_dup, i)
         k12_dup = _PM.var(pm, nw, :k12_dup, i)
         variable_on_off_switch(pm, k11, k11_dup, M, z1)
         variable_on_off_switch(pm, k12, k12_dup, M, z1)

         JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
         JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
         JuMP.@constraint(pm.model,k12_dup^2 <=  k11_dup)
         ############## frequency minimum in interval 2 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e  <=  Pf + (Pg/Tg)*Tcl + M*(1-z21))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tcl - M*(z21))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tf <=   Phvdcoaux[i]  + M*(1-z22))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tf  >= e +  Phvdcoaux[i]  - M*(z22))
         JuMP.@constraint(pm.model,  0 <= z21 + z22 - 2*z2)
         JuMP.@constraint(pm.model,  z21 + z22 - 2*z2 <= 1)

         k21 = _PM.var(pm, nw, :k21, i)
         k22 = _PM.var(pm, nw, :k22, i)
         k23 = _PM.var(pm, nw, :k23, i)

         k21_dup = _PM.var(pm, nw, :k21_dup, i)
         k22_dup = _PM.var(pm, nw, :k22_dup, i)
         k23_dup = _PM.var(pm, nw, :k23_dup, i)
         variable_on_off_switch(pm, k21, k21_dup, M, z2)
         variable_on_off_switch(pm, k22, k22_dup, M, z2)
         variable_on_off_switch(pm, k23, k23_dup, M, z2)

         JuMP.@constraint(pm.model,k21 ==  (H/50 - Pf*Tf/4))
         JuMP.@constraint(pm.model,k22 ==  Pg/Tg)
         JuMP.@constraint(pm.model,k23 ==  (Phvdcoaux[i] - Pf)/2)

         JuMP.@constraint(pm.model,[k21_dup/sqrt(2), k22_dup/sqrt(2), k23_dup, 0] in JuMP.RotatedSecondOrderCone() )

         ############## frequency minimum in interval 3 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e  <=   Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  + M*(1-z31))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  - M*(z31))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tcl <=   Phvdcoaux[i]  + M*(1-z32))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tcl  >= e+  Phvdcoaux[i]  - M*(z32))
         JuMP.@constraint(pm.model,  0 <= z31 + z32 - 2*z3)
         JuMP.@constraint(pm.model,  z31 + z32 - 2*z3 <= 1)

         k31 = _PM.var(pm, nw, :k31, i)
         k32 = _PM.var(pm, nw, :k32, i)
         k33 = _PM.var(pm, nw, :k33, i)
         k31_dup = _PM.var(pm, nw, :k31_dup, i)
         k32_dup = _PM.var(pm, nw, :k32_dup, i)
         k33_dup = _PM.var(pm, nw, :k33_dup, i)
         variable_on_off_switch(pm, k31, k31_dup, M, z3)
         variable_on_off_switch(pm, k32, k32_dup, M, z3)
         variable_on_off_switch(pm, k33, k33_dup, M, z3)

         JuMP.@constraint(pm.model, k31 ==  (H/50 - Pf*Tf/4 + Phvdccaux[i]*Tcl^2/(Td*4) ) )
         JuMP.@constraint(pm.model, k32 ==  (Pg/Tg + Phvdccaux[i]/Td) )
         JuMP.@constraint(pm.model, k33 ==  (Phvdcoaux[i] - Pf + Phvdccaux[i]*Tcl/Td)/2 )

         display(JuMP.@constraint(pm.model,[k31_dup/sqrt(2), k32_dup/sqrt(2), k33_dup, 0] in JuMP.RotatedSecondOrderCone() ))

        ############## frequency minimum in interval 4 ##########################################
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  Phvdccaux[i]  + M*(1-z41))
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg +  Phvdccaux[i]  - M*(z41))
        JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  <=   Phvdcoaux[i]  + M*(1-z42))
        JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]   >= e +  Phvdcoaux[i]  - M*(z42))
        JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
        JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

        k41 = _PM.var(pm, nw, :k41, i)
        k42 = _PM.var(pm, nw, :k42, i)
        k43 = _PM.var(pm, nw, :k43, i)
        k41_dup = _PM.var(pm, nw, :k41_dup, i)
        k42_dup = _PM.var(pm, nw, :k42_dup, i)
        k43_dup = _PM.var(pm, nw, :k43_dup, i)
        variable_on_off_switch(pm, k41, k41_dup, M, z4)
        variable_on_off_switch(pm, k42, k42_dup, M, z4)
        variable_on_off_switch(pm, k43, k43_dup, M, z4)

        JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4 - Phvdccaux[i]*(Td + 2*Tcl)/4) )
        JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
        JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf - Phvdccaux[i])/2 )
        JuMP.@constraint(pm.model,[k41_dup/sqrt(2), k42_dup/sqrt(2), k43_dup, 0] in JuMP.RotatedSecondOrderCone() )

         JuMP.@constraint(pm.model, z1 +z2 + z3+ z4 == 1)
     end
end

function constraint_frequency_stab_TNEP(pm::_PM.AbstractPowerModel)
    if pm.setting["Permanentloss"] == true
        constraint_frequency_stab_TNEP_MP_PL(pm)
    elseif pm.setting["FSprotection"] == true || pm.setting["NSprotection"] == true
        constraint_frequency_stab_TNEP_MP_FSNS(pm)
    end
end

function constraint_frequency_stab_TNEP_MP_PL(pm::_PM.AbstractPowerModel)
        base_nws = pm.setting["base_list"]
        cont_nws = pm.setting["Cont_list"]
		ev_syncarea = pm.setting["syncarea"]

        for (base, n, br) in cont_nws
        nw = n
        conv = _PM.ref(pm, :convdc)
        load = _PM.ref(pm, nw, :load, 1)
		syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn_ne)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Phvdcoaux_dup = _PM.var(pm, nw)[:Phvdcoaux_dup]
        Phvdccaux_dup = _PM.var(pm, nw)[:Phvdccaux_dup]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr_ne)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr_ne)
        Zb1 = _PM.var(pm, nw)[:zb1]

        i = ev_syncarea
        # for i in syncarea
        # JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv_ne, i)) / load["pd"])
        # JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv_ne, i)) / load["pd"])
        # JuMP.@constraint(pm.model,  Phvdccaux[i]  == 0)

        display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_incident_ne, br)) / load["pd"]) ) #sync. area is not implemented
        display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_incident_ne, br)) /  load["pd"]) ) #sync. area is not implemented
        for a in _PM.ref(pm, nw, :bus_arcs_incident_ne, br)
            display(JuMP.@constraint(pm.model,  Pconv2[a] == 0 ) ) #chagne in OPF but shouldnt make too much difference. The problem here is that it causes other converters not to build
        end
        # JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4])
        # JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])

         #
        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
        zbr = _PM.var(pm, nw, :branchdc_ne, br)

        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >=  Phvdcoaux[i] - JuMP.upper_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  <=  Phvdcoaux[i] - JuMP.lower_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >=  JuMP.lower_bound(Phvdcoaux[i])*(zbr) )
        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  <=  JuMP.upper_bound(Phvdcoaux[i])*(zbr) )

        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  >=  Phvdccaux[i] - JuMP.upper_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  <=  Phvdccaux[i] - JuMP.lower_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  >=  JuMP.lower_bound(Phvdccaux[i])*(zbr) )
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  <=  JuMP.upper_bound(Phvdccaux[i])*(zbr) )

        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        H = reserves[bi_bp[i]]["H"]
        Tcl = reserves[bi_bp[i]]["Tcl"]
        M = 1 # largest possible Pl=0.15 for now
        z1 = _PM.var(pm, nw, :z1, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))
        M = 1

        ####################################################################################
         JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux_dup[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux_dup[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

         k11 = _PM.var(pm, nw, :k11, i)
         k12 = _PM.var(pm, nw, :k12, i)
         k11_dup = _PM.var(pm, nw, :k11_dup, i)
         k12_dup = _PM.var(pm, nw, :k12_dup, i)
         variable_on_off_switch(pm, k11, k11_dup, M, z1)
         variable_on_off_switch(pm, k12, k12_dup, M, z1)

         JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
         JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux_dup[i])/2)
         JuMP.@constraint(pm.model,k12_dup^2 <=  k11_dup)

       ####################################################################################
       JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  + e <=   Pf + Pg +  M*(1-z41))
       JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >= Pf + Pg   - M*(z41))
       JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tf)   <=   Phvdcoaux_dup[i]  + M*(1-z42))
       JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tf)   >= e +  Phvdcoaux_dup[i]  - M*(z42))
       JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
       JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

       k41 = _PM.var(pm, nw, :k41, i)
       k42 = _PM.var(pm, nw, :k42, i)
       k43 = _PM.var(pm, nw, :k43, i)
       k41_dup = _PM.var(pm, nw, :k41_dup, i)
       k42_dup = _PM.var(pm, nw, :k42_dup, i)
       k43_dup = _PM.var(pm, nw, :k43_dup, i)
       variable_on_off_switch(pm, k41, k41_dup, M, z4)
       variable_on_off_switch(pm, k42, k42_dup, M, z4)
       variable_on_off_switch(pm, k43, k43_dup, M, z4)

       JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4) )
       JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
       JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux_dup[i] - Pf )/2 )
       JuMP.@constraint(pm.model,[k41_dup/sqrt(2), k42_dup/sqrt(2), k43_dup, 0] in JuMP.RotatedSecondOrderCone() )

       JuMP.@constraint(pm.model, z1 + z4 == 1)
        end
    end

function constraint_frequency_stab_TNEP_MP_FSNS(pm::_PM.AbstractPowerModel)
    base_nws = pm.setting["base_list"]
    cont_nws = pm.setting["Cont_list"]
    ev_syncarea = pm.setting["syncarea"]
    for (base, nw, br) in cont_nws
        load = _PM.ref(pm, nw, :load, 1) # get index later on
        syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn_ne)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Phvdcoaux_dup = _PM.var(pm, nw)[:Phvdcoaux_dup]
        Phvdccaux_dup = _PM.var(pm, nw)[:Phvdccaux_dup]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr_ne)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr_ne)
        Zb1 = _PM.var(pm, nw)[:zb1]

        i = ev_syncarea
        conv_conn = _PM.ref(pm, nw, :bus_arcs_conv_ne, i)
        Mmax = 72/ load["pd"] # 72MW = 24MW*3 highest onshroe rating for all test cases
        if pm.setting["FSprotection"] == true
            display(_PM.ref(pm, nw, :bus_convs_dirct_ne) )
            display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_convs_dirct_ne, i)) / load["pd"]) )
            display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_convs_dirct_ne, i)) /load["pd"]) )
            # for k in conv_conn
            #     display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= - (Pconv1[k]) /  load["pd"]) )
            #     display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  <= - (Pconv1[k]) / load["pd"] + Mmax*(1-Zb1[k]) ) )
            #     display(JuMP.@constraint(pm.model,  Phvdccaux[i]  >= - Pconv2[k] /  load["pd"] -  Mmax*(1-Zb1[k]) ) )
            #     display(JuMP.@constraint(pm.model,  Phvdccaux[i]  <= - Pconv2[k] / load["pd"] +   Mmax*(1-Zb1[k]) ) )
            # end
            JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4])
            JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])
            # display(JuMP.@constraint(pm.model,  sum(Zb1[k] for k in conv_conn)==1 ) )
            # display(JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] ) ) # what if a conv is not built, find a way around to know what is built
        elseif pm.setting["NSprotection"] == true
            #old
            # display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv_ne, i)) / load["pd"]) )
            # display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv_ne, i)) /  load["pd"]) )
            #new
            display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_incident_ne, br)) / load["pd"]) ) #sync. area is not implemented
            display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_incident_ne, br)) /  load["pd"]) ) #sync. area is not implemented
            # display(JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] ) ) # what if a conv not built, find a way around to know what is built
            JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4])
            JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])
        end

        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
        zbr = _PM.var(pm, nw, :branchdc_ne, br)

        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >=  Phvdcoaux[i] - JuMP.upper_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  <=  Phvdcoaux[i] - JuMP.lower_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >=  JuMP.lower_bound(Phvdcoaux[i])*(zbr) )
        JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  <=  JuMP.upper_bound(Phvdcoaux[i])*(zbr) )
        #
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  >=  Phvdccaux[i] - JuMP.upper_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  <=  Phvdccaux[i] - JuMP.lower_bound(Phvdcoaux[i])*(1-zbr) )
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  >=  JuMP.lower_bound(Phvdccaux[i])*(zbr) )
        JuMP.@constraint(pm.model,  Phvdccaux_dup[i]  <=  JuMP.upper_bound(Phvdccaux[i])*(zbr) )

        #check if pf and pg need to be disabled when a branch is not built since the FFR_list still accounts the pf,pg for nonbuilt branches..
        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        H = reserves[bi_bp[i]]["H"]
        Tcl = reserves[bi_bp[i]]["Tcl"]
        M = 1 # largest possible Pl=0.15 for now
        z1 = _PM.var(pm, nw, :z1, i)
        z2 = _PM.var(pm, nw, :z2, i)
        z3 = _PM.var(pm, nw, :z3, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z21 = _PM.var(pm, nw, :z21, i)
        z31 = _PM.var(pm, nw, :z31, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z22 = _PM.var(pm, nw, :z22, i)
        z32 = _PM.var(pm, nw, :z32, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))
        M = 1

        ############## frequency minimum in interval 1 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux_dup[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux_dup[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

         k11 = _PM.var(pm, nw, :k11, i)
         k12 = _PM.var(pm, nw, :k12, i)
         k11_dup = _PM.var(pm, nw, :k11_dup, i)
         k12_dup = _PM.var(pm, nw, :k12_dup, i)
         variable_on_off_switch(pm, k11, k11_dup, M, z1)
         variable_on_off_switch(pm, k12, k12_dup, M, z1)

         JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
         JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux_dup[i])/2)
         JuMP.@constraint(pm.model,k12_dup^2 <=  k11_dup)

         ############## frequency minimum in interval 2 ##########################################
           JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  +e  <=  Pf + (Pg/Tg)*Tcl + M*(1-z21))
           JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >= Pf + (Pg/Tg)*Tcl - M*(z21))
           JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tf <=   Phvdcoaux_dup[i]  + M*(1-z22))
           JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tf  >= e +  Phvdcoaux_dup[i]  - M*(z22))
           JuMP.@constraint(pm.model,  0 <= z21 + z22 - 2*z2)
           JuMP.@constraint(pm.model,  z21 + z22 - 2*z2 <= 1)

           k21 = _PM.var(pm, nw, :k21, i)
           k22 = _PM.var(pm, nw, :k22, i)
           k23 = _PM.var(pm, nw, :k23, i)
           k21_dup = _PM.var(pm, nw, :k21_dup, i)
           k22_dup = _PM.var(pm, nw, :k22_dup, i)
           k23_dup = _PM.var(pm, nw, :k23_dup, i)
           variable_on_off_switch(pm, k21, k21_dup, M, z2)
           variable_on_off_switch(pm, k22, k22_dup, M, z2)
           variable_on_off_switch(pm, k23, k23_dup, M, z2)

           JuMP.@constraint(pm.model,k21 ==  (H/50 - Pf*Tf/4))
           JuMP.@constraint(pm.model,k22 ==  Pg/Tg)
           JuMP.@constraint(pm.model,k23 ==  (Phvdcoaux_dup[i] - Pf)/2)

           JuMP.@constraint(pm.model,[k21_dup/sqrt(2), k22_dup/sqrt(2), k23_dup, 0] in JuMP.RotatedSecondOrderCone() )

           ############## frequency minimum in interval 3 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  +e  <=   Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux_dup[i]  + M*(1-z31))
         JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >= Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux_dup[i]  - M*(z31))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tcl <=   Phvdcoaux_dup[i]  + M*(1-z32))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tcl  >= e+  Phvdcoaux_dup[i]  - M*(z32))
         JuMP.@constraint(pm.model,  0 <= z31 + z32 - 2*z3)
         JuMP.@constraint(pm.model,  z31 + z32 - 2*z3 <= 1)

         k31 = _PM.var(pm, nw, :k31, i)
         k32 = _PM.var(pm, nw, :k32, i)
         k33 = _PM.var(pm, nw, :k33, i)
         k31_dup = _PM.var(pm, nw, :k31_dup, i)
         k32_dup = _PM.var(pm, nw, :k32_dup, i)
         k33_dup = _PM.var(pm, nw, :k33_dup, i)
         variable_on_off_switch(pm, k31, k31_dup, M, z3)
         variable_on_off_switch(pm, k32, k32_dup, M, z3)
         variable_on_off_switch(pm, k33, k33_dup, M, z3)

         JuMP.@constraint(pm.model, k31 ==  (H/50 - Pf*Tf/4 + Phvdccaux_dup[i]*Tcl^2/(Td*4) ) )
         JuMP.@constraint(pm.model, k32 ==  (Pg/Tg + Phvdccaux_dup[i]/Td) )
         JuMP.@constraint(pm.model, k33 ==  (Phvdcoaux_dup[i] - Pf + Phvdccaux_dup[i]*Tcl/Td)/2 )

         JuMP.@constraint(pm.model,[k31_dup/sqrt(2), k32_dup/sqrt(2), k33_dup, 0] in JuMP.RotatedSecondOrderCone() )

         ############## frequency minimum in interval 4 ##########################################
           JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  + e <=   Pf + Pg +  Phvdccaux_dup[i]  + M*(1-z41))
           JuMP.@constraint(pm.model,  Phvdcoaux_dup[i]  >= Pf + Pg +  Phvdccaux_dup[i]  - M*(z41))
           JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux_dup[i]  <=   Phvdcoaux_dup[i]  + M*(1-z42))
           JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux_dup[i]   >= e +  Phvdcoaux_dup[i]  - M*(z42))
           JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
           JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

           k41 = _PM.var(pm, nw, :k41, i)
           k42 = _PM.var(pm, nw, :k42, i)
           k43 = _PM.var(pm, nw, :k43, i)
           k41_dup = _PM.var(pm, nw, :k41_dup, i)
           k42_dup = _PM.var(pm, nw, :k42_dup, i)
           k43_dup = _PM.var(pm, nw, :k43_dup, i)
           variable_on_off_switch(pm, k41, k41_dup, M, z4)
           variable_on_off_switch(pm, k42, k42_dup, M, z4)
           variable_on_off_switch(pm, k43, k43_dup, M, z4)

           JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4 - Phvdccaux_dup[i]*(Td + 2*Tcl)/4) )
           JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
           JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux_dup[i] - Pf - Phvdccaux_dup[i])/2 )
           JuMP.@constraint(pm.model,[k41_dup/sqrt(2), k42_dup/sqrt(2), k43_dup, 0] in JuMP.RotatedSecondOrderCone() )

        JuMP.@constraint(pm.model, z1 +z2 + z3+ z4 == 1)

    end
end

function constraint_frequency_stab_OPF_MP_PLdim(pm::_PM.AbstractPowerModel)
    base_nws = pm.setting["base_list"]
    cont_nws = pm.setting["Cont_list"]
    ev_syncarea = pm.setting["syncarea"]

    for (base, nw, br) in cont_nws
        load = _PM.ref(pm, nw, :load, 1) # get index later on
        syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
        Zb1 = _PM.var(pm, nw)[:zb1]

        # for i in syncarea # when there is more than one synchornous areas
        i = ev_syncarea
        display("connstraint: post_acdcscopf_nocl_PLdim")
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
        JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
        JuMP.@constraint(pm.model,  Phvdccaux[i]  == 0)
        JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] )
        JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])

        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        H = reserves[bi_bp[i]]["H"]

        M = 1 # largest possible Pl=0.15 for now

        z1 = _PM.var(pm, nw, :z1, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))

        ############## frequency minimum in interval 1 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

         k11 = _PM.var(pm, nw, :k11, i)
         k12 = _PM.var(pm, nw, :k12, i)

         JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
         JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
         JuMP.@constraint(pm.model,k12^2 <=  k11)

         ############## frequency minimum in interval 2 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  M*(1-z41))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg   - M*(z41))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tf)   <=   Phvdcoaux[i]  + M*(1-z42))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tf)   >= e +  Phvdcoaux[i]  - M*(z42))
         JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
         JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

         k41 = _PM.var(pm, nw, :k41, i)
         k42 = _PM.var(pm, nw, :k42, i)
         k43 = _PM.var(pm, nw, :k43, i)

         JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4) )
         JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
         JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf )/2 )
         JuMP.@constraint(pm.model,[k41/sqrt(2), k42/sqrt(2), k43, 0] in JuMP.RotatedSecondOrderCone() )

         JuMP.@constraint(pm.model, z1 + z4 == 1)
     end
end


text:
using Augmentor, ISICArchive, Plots
pyplot(reuse=true);
mkpath(fname) = joinpath(dirname(@__FILE__()), fname)

srand(123)

img = get(ImageThumbnailRequest(id = "5592ac599fc3c13155a57a85"))

pipeline =
    Either(1=>FlipX(), 1=>FlipY(), 2=>NoOp()) |>
    Rotate(0:360) |>
    ShearX(-5:5) * ShearY(-5:5) |>
    CropSize(165, 165) |>
    Zoom(1:0.05:1.2) |>
    Resize(64, 64)

# Create image that shows the input
plot(img, size=(256,169), xlim=(1,255), ylim=(1,168), grid=false, ticks=true)
Plots.png(mkpath("readme/isic_in.png"))

# create animate gif that shows 10 outputs
anim = @animate for i=1:10
    plot(augment(img, pipeline), size=(169,169), xlim=(1,63), ylim=(1,63), grid=false, ticks=true);
end
Plots.gif(anim, mkpath("readme/isic_out.gif"), fps = 2)

max_stars_repo_name: JuliaTagBot/FileStorage

max_issues_repo_name: JuliaTagBot/FileStorage

max_forks_repo_name: JuliaTagBot/FileStorage


text:
@testset "Fetch Function Tests" begin

  @test isdefined(SetAndForget, :fetch) == true

end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using GLMakie
GLMakie.activate!() #HIDE
let
    t = 0:0.1:15
    u = -1:0.1:1
    x = [u*sin(t) for t in t, u in u]
    y = [u*cos(t) for t in t, u in u]
    z = [t/4      for t in t, u in u]
    fig, _ = surface(x, y, z, colorrange = (-2,-1.1), highclip = :orangered,
    lightposition = Vec3f0(0, 0, 0), ambient = Vec3f0(0.65, 0.65, 0.65),
    backlight = 5f0) # the colorrange must be outside the range of z
    wireframe!(x,y,z, overdraw = false, linewidth = 0.1) # try overdraw = true
    fig
    save(joinpath(@__DIR__, "output", "surface1.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Dates: Year, Date

toDate(dstr::String)::Date = Date(dstr, "mm/dd/yy") + Year(2000)

# Same dates for all test cases
# Test dates and values from https://www.isda.org/2008/12/22/30-360-day-count-conventions/ , https://www.isda.org/a/mIJEE/30-360-2006ISDADefs.xls
startDates = map(toDate, ["01/15/07", "01/15/07", "01/15/07", "09/30/07", "09/30/07", "09/30/07", "01/15/07", "01/31/07", "02/28/07", "08/31/06", "02/28/07", "02/14/07", "02/26/07", "02/29/08", "02/29/08", "02/29/08", "02/28/07", "10/31/07", "08/31/07", "02/29/08", "08/31/08", "02/28/09"])
endDates   = map(toDate, ["01/30/07", "02/15/07", "07/15/07", "03/31/08", "10/31/07", "09/30/08", "01/31/07", "02/28/07", "03/31/07", "02/28/07", "08/31/07", "02/28/07", "02/29/08", "02/28/09", "03/30/08", "03/31/08", "03/05/07", "11/28/07", "02/29/08", "08/31/08", "02/28/09", "08/31/09"])
terminationDate = toDate("02/28/09")

function test_helper(daysResults, daysFracResults, ::Type{T}) where T <: Learn.AbstractDayCountConvention
    for (sDate, eDate, daysR, yearFracR) in zip(startDates, endDates, daysResults, daysFracResults)
        @test days(sDate, eDate, T) == daysR
        @test yearfraction(sDate, eDate, T) ≈ yearFracR atol=0.0001
    end
end

function test_helper(daysResults, daysFracResults, tDate, ::Type{T}) where T <: Learn.AbstractDayCountConvention
    for (sDate, eDate, daysR, yearFracR) in zip(startDates, endDates, daysResults, daysFracResults)
        @test days(sDate, eDate, tDate, T) == daysR
        @test yearfraction(sDate, eDate, tDate, T) ≈ yearFracR atol=0.0001
    end
end

@testset "DayCounts" begin

    @testset "30/360 Bond Basis, 30A/360" begin
        # Test cases from https://www.isda.org/a/mIJEE/30-360-2006ISDADefs.xls
        testType = DayCount30A_360
        daysResults = [15,30,180,180,30,360,16,28,33,178,183,14,363,359,31,32,7,28,179,182,178,183]
        daysFracResults = [0.0417,0.0833,0.5000,0.5000,0.0833,1.0000,0.0444,0.0778,0.0917,0.4944,0.5083,0.0389,1.0083,0.9972,0.0861,0.0889,0.0194,0.0778,0.4972,0.5056,0.4944,0.5083]
        test_helper(daysResults, daysFracResults, testType)
    end

    @testset "30E/360, 30/360 ICMA, 30/360 ISMA, 30S/360, Eurobond basis (ISDA 2006), Special German" begin
        # Test cases from https://www.isda.org/a/mIJEE/30-360-2006ISDADefs.xls
        testType = DayCount30E_360
        daysResults = [15, 30, 180, 180, 30, 360, 15, 28, 32, 178, 182, 14, 363, 359, 31, 31, 7, 28, 179, 181, 178, 182]
        daysFracResults = [0.0417, 0.0833, 0.5000, 0.5000, 0.0833, 1.0000, 0.0417, 0.0778, 0.0889, 0.4944, 0.5056, 0.0389, 1.0083, 0.9972, 0.0861, 0.0861, 0.0194, 0.0778, 0.4972, 0.5028, 0.4944, 0.5056]
        test_helper(daysResults, daysFracResults, testType)
    end

    @testset "30E/360ISDA, Eurobond basis (ISDA 2000), German" begin
        # Test cases from https://www.isda.org/a/mIJEE/30-360-2006ISDADefs.xls
        testType = DayCount30E_360ISDA
        daysResults = [15, 30, 180, 180, 30, 360, 15, 30, 30, 180, 180, 16, 364, 358, 30, 30, 5, 28, 180, 180, 178, 180]
        daysFracResults = [0.0417, 0.0833, 0.5000, 0.5000, 0.0833, 1.0000, 0.0417, 0.0833, 0.0833, 0.5000, 0.5000, 0.0444, 1.0111, 0.9944, 0.0833, 0.0833, 0.0139, 0.0778, 0.5000, 0.5000, 0.4944, 0.5000]
        test_helper(daysResults, daysFracResults, terminationDate, testType)
    end

    @testset "Act/360, French" begin
        # Test cases from https://www.isda.org/a/mIJEE/30-360-2006ISDADefs.xls
        testType = DayCountAct_360
        daysResults = [15, 31, 181, 183, 31, 366, 16, 28, 31, 181, 184, 14, 368, 365, 30, 31, 5, 28, 182, 184, 181, 184]
        daysFracResults = [0.0417, 0.0861, 0.5028, 0.5083, 0.0861, 1.0167, 0.0444, 0.0778, 0.0861, 0.5028, 0.5111, 0.0389, 1.0222, 1.0139, 0.0833, 0.0861, 0.0139, 0.0778, 0.5056, 0.5111, 0.5028, 0.5111]
        test_helper(daysResults, daysFracResults, testType)
    end

    @testset "Act/364" begin
        testType = DayCountAct_364
        daysResults = [15, 31, 181, 183, 31, 366, 16, 28, 31, 181, 184, 14, 368, 365, 30, 31, 5, 28, 182, 184, 181, 184]
        daysFracResults = map(x-> x*360.0/364.0, [0.0417, 0.0861, 0.5028, 0.5083, 0.0861, 1.0167, 0.0444, 0.0778, 0.0861, 0.5028, 0.5111, 0.0389, 1.0222, 1.0139, 0.0833, 0.0861, 0.0139, 0.0778, 0.5056, 0.5111, 0.5028, 0.5111])
        test_helper(daysResults, daysFracResults, testType)
    end

    @testset "Act/365Fixed" begin
        testType = DayCountAct_365Fixed
        daysResults = [15, 31, 181, 183, 31, 366, 16, 28, 31, 181, 184, 14, 368, 365, 30, 31, 5, 28, 182, 184, 181, 184]
        daysFracResults = map(x-> x*360.0/365.0, [0.0417, 0.0861, 0.5028, 0.5083, 0.0861, 1.0167, 0.0444, 0.0778, 0.0861, 0.5028, 0.5111, 0.0389, 1.0222, 1.0139, 0.0833, 0.0861, 0.0139, 0.0778, 0.5056, 0.5111, 0.5028, 0.5111])
        test_helper(daysResults, daysFracResults, testType)
    end

    @testset "Act/365Act" begin
        testType = DayCountAct_365Actual
        daysResults = [15, 31, 181, 183, 31, 366, 16, 28, 31, 181, 184, 14, 368, 365, 30, 31, 5, 28, 182, 184, 181, 184]
        daysFracResults = map(x-> x*360.0/365.0, [0.0417, 0.0861, 0.5028, 0.5083, 0.0861, 1.0167, 0.0444, 0.0778, 0.0861, 0.5028, 0.5111, 0.0389, 1.0222, 1.0139, 0.0833, 0.0861, 0.0139, 0.0778, 0.5056, 0.5111, 0.5028, 0.5111])

        for i in [4, 6, 13, 19]
            daysFracResults[i] = daysFracResults[i]*365.0/366.0
        end
        test_helper(daysResults, daysFracResults, testType)
    end

    # @testset "Act/ActISDA" begin
    #     testType = DayCountAct_Act_ISDA
    #     daysResults = [15, 31, 181, 183, 31, 366, 16, 28, 31, 181, 184, 14, 368, 365, 30, 31, 5, 28, 182, 184, 181, 184]
    #     daysFracResults = map(x-> x*360.0/365.0, [0.0417, 0.0861, 0.5028, 0.5083, 0.0861, 1.0167, 0.0444, 0.0778, 0.0861, 0.5028, 0.5111, 0.0389, 1.0222, 1.0139, 0.0833, 0.0861, 0.0139, 0.0778, 0.5056, 0.5111, 0.5028, 0.5111])

    #     test_helper(daysResults, daysFracResults, testType)
    # end

end

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
using Microscopy
using Test

include("microscope_modes.jl")

@testset "Microscopy.jl" begin
end

max_stars_repo_name: JuliaMicroscopy/Microscopy

max_forks_repo_name: JuliaMicroscopy/Microscopy


text:
Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg) = _check_index_any(axis, arg)
_check_index_any(axis, arg) = checkindex(Bool, parent(axis), arg)
_check_index_any(axis::Axis, arg) = in(arg, keys(axis))

Base.checkindex(::Type{Bool}, axis::AbstractAxis, ::Interval) = true
Base.checkindex(::Type{Bool}, axis::AbstractAxis, ::Colon) = true
Base.checkindex(::Type{Bool}, axis::AbstractAxis, ::Slice) = true
function Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::AbstractArray{T}) where {T}
    return _check_index_array(axis, arg)
end
_check_index_array(axis, arg) = checkindex(Bool, eachindex(axis), arg)
function _check_index_array(axis::Axis, arg::AbstractArray{T}) where {T}
    if T <: Integer
        return checkindex(Bool, eachindex(axis), arg)
    else
        return all(in(keys(axis)), arg)
    end
end

Base.checkindex(::Type{Bool}, axis::AbstractAxis, ::AbstractArray{Bool}) = false
function Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::AbstractVector{Bool})
    return checkindex(Bool, eachindex(axis), arg)
end
Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::Real) = _check_index_real(axis, arg)
_check_index_real(axis, arg) = checkindex(Bool, parent(axis), arg)
function _check_index_real(axis::AbstractOffsetAxis, arg)
    return checkindex(Bool, parent(axis), _sub_offset(axis, arg))
end
function _check_index_real(axis::Axis, arg)
    if arg isa Integer
        return checkindex(Bool, parent(axis), arg)
    else
        return _check_index_any(axis, arg)
    end
end
function Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::Union{<:Equal,<:Approx})
    return !(find_first(arg, keys(axis)) === nothing)
end
function Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::AbstractRange{T}) where {T}
    return _check_index_range(axis, arg)
end
_check_index_range(axis, arg) = checkindex(Bool, parent(axis), arg)
function _check_index_range(axis::AbstractOffsetAxis, arg)
    return checkindex(Bool, parent(axis), _sub_offset(axis, arg))
end
function _check_index_range(axis::Axis, arg::AbstractRange{T}) where {T}
    if T <: Integer
        return checkindex(Bool, parent(axis), arg)
    else
        return length(find_all_in(arg, keys(axis))) == length(axis)
    end
end

Base.checkindex(::Type{Bool}, axis::AbstractAxis, ::Fix2) = true
function Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::LogicalIndex)
    return (axis,) == axes(arg.mask)
end


text:
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Swifter
import Swifter: chains

@query vc.view

text = "Hello world"
@query vc.label.text = text

@test isdefined(Main, :text)

@query intext = "inside"
@query vc.label.text = intext

@test Any[:text] == chains(:text)

pairs = chains(:(UIFont(name: "Courier", size: 90)))
@test Any[(:call, (:UIFont,Any[Any[:name,"Courier"],Any[:size,90]]))] == pairs

outrect = "{{33, 50}, {531, 200}}"
@query vc.label.frame = outrect

@query inrect = "{{33, 50}, {131, 200}}"
@test "{{33, 50}, {131, 200}}" == Swifter.env[:inrect]

@query inrect = outrect
@test "{{33, 50}, {531, 200}}" == Swifter.env[:inrect]

@query rect = CGRectMake(1,2,5,6)
@test (:call,(:CGRectMake,Any[1,2,5,6])) == Swifter.env[:rect]

@query vc.view.subviews[0].backgroundColor = UIColor.redColor()

@test Any[(:call,:tap)] == chains(:(tap()))

@test Any[:view, (:call,:tap)] == chains(:(view.tap()))

@query vc.view.alpha = 0.3

if VERSION >= v"0.5-" @eval begin

@query intext = text
@test "Hello world" == @query intext
@query vc.label.text = intext

end end


text:
# simple example of a harmonic power flow execution
using HarmonicPowerFlow
using StatsPlots
using LinearAlgebra

harmonics = [h for h in 1:2:51]
# settings
coupled = false
settings1 = init_settings(coupled, harmonics)

# import nodes and lines from csv files
net3 = init_power_grid(
    import_nodes_from_csv("examples\\multiple_pv\\net3_buses.csv"), import_lines_from_csv("examples\\multiple_pv\\net3_lines.csv"),
    settings1)

# run harmonic power flow
u, err_h_final, n_iter_h = hpf(net3, settings1)

THD = HarmonicPowerFlow.THD(u)

# show barplot of harmonics at bus 4
HarmonicPowerFlow.barplot_THD(u, 4, h_max= 41)

# Analyze NE parameters
NE = HarmonicPowerFlow.import_Norton_Equivalents(net3.nodes, settings1)
l = @layout [a{0.66w}  b{1h}]
if coupled
    p1 = heatmap(harmonics, harmonics, abs.(NE["SMPS"][2]), ticks = harmonics, title = "Y_Norton Matrix (coupled)", yflip=true)
    p2 = heatmap(transpose(abs.(NE["SMPS"][1])), aspect_ratio=1, yflip=true, title = "I_Norton", yticks = nothing, border = :none)
    plot(p1, p2, layout = l)
else
    Y_N = heatmap(harmonics, harmonics, diagm(vec(abs.(NE["SMPS"][2]))), ticks = harmonics, title = "Y_Norton Vector (uncoupled)", yflip=true)
end


text:
#==============================================================================================#
#                            functions.jl -- TypeTree functionality                            #
#==============================================================================================#

using InteractiveUtils: subtypes

"""
`function tt(TY::Type, pref=("",); uni=true, concrete=true)`\n
Returns a `Vector{String}` of type entries formatted as a tree for human (programmer)
inspection.

- The `TY` argument is the type for which to produce the tree;
- The `pref` optional argument is a text prefix;
- The `uni` optional keyword argument controls whether Unicode characters will be used;
- The `concrete` optional keyword argument controls whether to filter off concrete subtypes.
"""
function tt(TY::Type, pref=("",); uni=true, concrete=true)
    ret = Array{String,1}()
    ELL, FRK, BAR = uni ? (" └─ ", " ├─ ", " │  ") : (" \\-- ", " +-- ", " |  ")
    SPC = "    "
    LEN = length(pref) - 1
    PREF = LEN > 0 ? ( i == ELL ? SPC : i == FRK ? BAR : i for i in pref[1:LEN]) : ("",)
    append!(ret, ["$(join((PREF..., pref[end])))$(TY)\n"])
    ST = [i for i in subtypes(TY) if isabstracttype(i) || concrete]
    LE = length(ST)
    for i in 1:LE
        append!(
            ret, tt(ST[i], (pref..., i < LE ? FRK : ELL), uni=uni, concrete=concrete)
        )
    end
    ret
end


text:
export translate, translate_optim

function translationclosure(mm::AbstractArray{T,N}, movement) where {T,N}
    function translation_loss(matrix)
        dims = N - 1
        c = reshape(matrix, dims, dims) * movement
        loss = zero(eltype(mm))
        for (i, j) in enumerate(axes(mm, N))
            idx = c[:, i]
            checkbounds(Bool, mm, idx..., j) || return typemax(eltype(loss))
            loss += mm(idx..., j)
        end
        return loss
    end
    return translation_loss
end

function translate(fixed, moving, maxshift, thresh = nothing)
    mm = RegisterMismatch.mismatch(fixed, moving, maxshift)
    _, denom = RegisterMismatch.separate(mm)
    if thresh == nothing
        thresh = 0.25maximum(denom)
    end
    return RegisterCore.ratio.(mm, thresh, Inf)
end

function translate_optim(
    mm::AbstractArray{T,N},
    initial_matrix,
    movement;
    kwargs...,
) where {T,N}
    # The final dimension of mm is assumed to represent "image number" and is not a continuous variable.
    # Hence we should not interpolate over it.
    qfc = BSpline(Quadratic(Free(OnCell())))
    itp_mm = interpolate(mm, (ntuple(i->qfc, N-1)..., NoInterp()))
    f = translationclosure(itp_mm, movement)
    result = optimize(f, initial_matrix, LBFGS(), Optim.Options(; kwargs...))
    Optim.converged(result) || @warn "Optimization failed to converge"
    return result
end


text:
struct Noaccel end

function update!(_::Noaccel, _, _) end

import Base: *

function (*)(L::Noaccel, v)
    w = similar(v)
    mul!(w, L, v)
end

import LinearAlgebra: mul!

mul!(d::T, _::Noaccel, v::T) where {T} = copyto!(d, v)


text:
# SPDX-FileCopyrightText: 2021 Johannes Hampp <johannes.hampp@zeu.jlug.de>
#
# SPDX-License-Identifier: MIT

# Activate julia environment
using Pkg
Pkg.activate("./workflow/envs")
using GlobalEnergyGIS

println("Downloading aux. datasets.")
download_datasets()

println("Rasterising datasets.")
rasterize_datasets(cleanup=:all)

max_stars_repo_name: euronion/synde

max_issues_repo_name: euronion/synde

max_forks_repo_name: euronion/synde


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using StaticArrays
using Distributions # to define priors
using Random        # to seed the random number generator
# Let's generate the data
# -----------------------
using Bridge
include(joinpath(SRC_DIR, "auxiliary", "data_simulation_fns.jl"))
Random.seed!(4)
#θ₀ = [10.0, 28.0, 8.0/3.0, 3.0, 3.0, 3.0]
θˣ = [10.0, 28.0, 8.0/3.0, 3.0]
Pˣ = LorenzCV(θˣ...)

x0, dt, T = ℝ{3}(1.5, -1.5, 25.0), 1/5000, 10.0
tt = 0.0:dt:T
XX, _ = simulate_segment(ℝ{3}(0.0, 0.0, 0.0), x0, Pˣ, tt)

θ₀ = [5.0, 15.0, 6.0, 8.0]
Pˣ = LorenzCV(θ₀...)

skip = 1000

Σdiagel = 10^0
Σ = SMatrix{2,2}(1.0I)*Σdiagel
L = @SMatrix[1.0 0.0 0.0;
             0.0 1.0 0.0]

obsTime, obsVals = XX.tt[1:skip:end], [rand(Gaussian(L*x, Σ)) for x in XX.yy[1:skip:end]]

fptOrPartObs = PartObs()
fpt = [NaN for _ in obsTime[2:end]]

auxFlag = Val{(true,true,false)}()
P̃ = [LorenzCVAux(θ₀..., t₀, u, T, v, auxFlag, x0[3]) for (t₀, T, u, v)
     in zip(obsTime[1:end-1], obsTime[2:end], obsVals[1:end-1], obsVals[2:end])]

Ls = [L for _ in P̃]
Σs = [Σ for _ in P̃]
τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
numSteps=1*10^4
saveIter=1*10^3

tKernel = RandomWalk([2.0, 1.0, 0.64, 0.3],
                     [false, false, false, true])

priors = Priors((ImproperPrior(), ImproperPrior(), ImproperPrior(),
                 ImproperPrior()))

𝔅 = NoBlocking()
blockingParams = ([], 0.1, NoChangePt())
changePt = NoChangePt()
#x0Pr = KnownStartingPt(x0)
x0Pr = GsnStartingPt(x0, x0, @SMatrix [20.0 0.0 0.0; 0.0 20.0 0.0; 0.0 0.0 400.0])
warmUp = 100

Random.seed!(4)
start = time()
(chain, accRateImp, accRateUpdt,
    paths, time_) = mcmc(eltype(x0), fptOrPartObs, obsVals, obsTime, x0Pr,
                         ℝ{3}(0.0, 0.0, 0.0), Pˣ, P̃, Ls, Σs, numSteps,
                         tKernel, priors, τ;
                         fpt=fpt,
                         ρ=0.995,
                         dt=1/2000,
                         saveIter=saveIter,
                         verbIter=10^2,
                         updtCoord=(Val((true, false, false, false)),
                                    Val((false, true, false, false)),
                                    Val((false, false, true, false)),
                                    Val((false, false, false, true))
                                    ),
                         paramUpdt=true,
                         updtType=(MetropolisHastingsUpdt(),
                                   MetropolisHastingsUpdt(),
                                   MetropolisHastingsUpdt(),
                                   MetropolisHastingsUpdt()
                                   ),
                         skipForSave=10^0,
                         blocking=𝔅,
                         blockingParams=blockingParams,
                         solver=Vern7(),
                         changePt=changePt,
                         warmUp=warmUp,
                         adaptiveProp=NoAdaptation())
elapsed = time() - start
print("time elapsed: ", elapsed, "\n")

print("imputation acceptance rate: ", accRateImp,
      ", parameter update acceptance rate: ", accRateUpdt)

using Plots
pTp = [[[x[i] for x in path] for path in paths] for i in 1:3]

function plotPaths(j, obsIdxS, obsIdxE, show_obs=true)
    idxS = div((obsIdxS-1)*skip,5)+1
    idxE = div((obsIdxE-1)*skip,5)+1
    p = plot()
    for i in 1:length(paths)
        plot!(time_[idxS:idxE], pTp[j][i][idxS:idxE], label="", color="steelblue", alpha=0.2, linewidth=0.2)
    end
    if show_obs
        scatter!(obsTime[obsIdxS:obsIdxE], [x[j] for x in obsVals][obsIdxS:obsIdxE],
                 color="orange", label="")
    end
    p
end

plotPaths(1, 1, 10)
plotPaths(2, 1, 10)
plotPaths(3, 1, 10, false)

plot([θ[1] for θ in chain])
plot([θ[2] for θ in chain])
plot([θ[3] for θ in chain])
plot([θ[4] for θ in chain])


text:
module DiagTree
using AbstractTrees
using Printf, PyCall, DataFrames

@enum TwoBodyChannel Alli = 1 PHr PHEr PPr AnyChan
@enum Permutation Di = 1 Ex DiEx

export TwoBodyChannel, Alli, PHr, PHEr, PPr, AnyChan
export Permutation, Di, Ex, DiEx

Base.length(r::TwoBodyChannel) = 1
Base.iterate(r::TwoBodyChannel) = (r, nothing)
function Base.iterate(r::TwoBodyChannel, ::Nothing) end

Base.length(r::Permutation) = 1
Base.iterate(r::Permutation) = (r, nothing)
function Base.iterate(r::Permutation, ::Permutation) end

include("common.jl")
include("traits.jl")
include("tree.jl")
include("operation.jl")
include("io.jl")
include("eval.jl")
include("optimize.jl")

const INL, OUTL, INR, OUTR = 1, 2, 3, 4

export Diagram
export addSubDiagram!
export toDataFrame
export evalDiagTree!
export Operator, Sum, Prod
export DiagramId, GenericId, Ver4Id, Ver3Id, GreenId, SigmaId, PolarId
export PropagatorId, BareGreenId, BareInteractionId
export BareGreenNId, BareHoppingId, GreenNId, ConnectedGreenNId
export uidreset, toDataFrame, mergeby, plot_tree

end


text:
module RigidBodyMotions

export RigidBodyMotion, Kinematics, d_dt, assign_velocity!

using DocStringExtensions
import ForwardDiff
import Base: +, *, -, >>, <<, show

"""
An abstract type for types that takes in time and returns `(c, ċ, c̈, α, α̇, α̈)`.
"""
abstract type Kinematics end

"""
An abstract type for real-valued functions of time.
"""
abstract type Profile end

"""
    RigidBodyMotion

A type to store the body's current kinematics

# Fields

- `c`: current centroid position (relative to initial position)
- `ċ`: current centroid velocity
- `c̈`: current centroid acceleration
- `α`: current angle (relative to initial angle)
- `α̇`: current angular velocity
- `α̈`: current angular acceleration
- `kin`: a [`Kinematics`](@ref) structure

The first six fields are meant as a cache of the current kinematics
while the `kin` field can be used to find the plate kinematics at any time.
"""
mutable struct RigidBodyMotion
    c::ComplexF64
    ċ::ComplexF64
    c̈::ComplexF64
    α::Float64
    α̇::Float64
    α̈::Float64

    kin::Kinematics
end

RigidBodyMotion(ċ, α̇) = RigidBodyMotion(0.0im, complex(ċ), 0.0im, 0.0, float(α̇),
                                          0.0, Constant(ċ, α̇))
RigidBodyMotion(kin::Kinematics) = RigidBodyMotion(kin(0)..., kin)
(m::RigidBodyMotion)(t) = m.kin(t)

function (m::RigidBodyMotion)(t,x̃::Tuple{Real,Real})
  # This expects coordinates in body's own coordinate system
  #
  z̃ = ComplexF64(x̃[1],x̃[2])
  m.c, m.ċ, m.c̈, m.α, m.α̇, m.α̈ = m.kin(t)
  z = exp(im*m.α)*z̃
  return m.c + z, m.ċ + im*m.α̇*z, m.c̈ + (im*m.α̈-m.α̇^2)*z
end

"""
    assign_velocity!(u::AbstractVector{Float64},v::AbstractVector{Float64},
                     x::AbstractVector{Float64},y::AbstractVector{Float64},
                     xc::Real,yc::Real,α::Real,
                     motion,t::Real)

Assign the components of rigid body velocity `u` and `v` (in inertial coordinate system)
at positions described by coordinates `x`, `y` (also in inertial coordinate system) at time `t`,
based on supplied motion `motion` for the body.
"""
function assign_velocity!(u::AbstractVector{Float64},v::AbstractVector{Float64},
                          x::AbstractVector{Float64},y::AbstractVector{Float64},
                          xc::Real,yc::Real,α::Real,m::RigidBodyMotion,t::Real)
   _,ċ,_,_,α̇,_ = m(t)
  for i = 1:length(x)
      Δz = (x[i]-xc)+im*(y[i]-yc)
      ċi = ċ + im*α̇*Δz
      u[i] = real(ċi)
      v[i] = imag(ċi)
  end
  nothing
end

function show(io::IO, m::RigidBodyMotion)
    println(io, "Rigid Body Motion:")
    println(io, "  ċ = $(round(m.ċ, digits=2))")
    println(io, "  c̈ = $(round(m.c̈, digits=2))")
    println(io, "  α̇ = $(round(m.α̇, digits=2))")
    println(io, "  α̈ = $(round(m.α̈, digits=2))")
    print(io, "  $(m.kin)")
end

#=
Kinematics
=#

struct Constant{C <: Complex, A <: Real} <: Kinematics
    ċ::C
    α̇::A
end
Constant(ċ, α̇) = Constant(complex(ċ), α̇)
(c::Constant{C})(t) where C = zero(C), c.ċ, zero(C), 0.0, c.α̇, 0.0
show(io::IO, c::Constant) = print(io, "Constant (ċ = $(c.ċ), α̇ = $(c.α̇))")

"""
    Pitchup <: Kinematics

Kinematics describing a pitchup motion (horizontal translation with rotation)

# Constructors
# Fields
$(FIELDS)
"""
struct Pitchup <: Kinematics
    "Freestream velocity"
    U₀::Float64
    "Axis of rotation, relative to the plate centroid"
    a::Float64

    "Non-dimensional pitch rate ``K = \\dot{\\alpha}_0\\frac{c}{2U_0}``"
    K::Float64

    "Initial angle of attack"
    α₀::Float64
    "Nominal start of pitch up"
    t₀::Float64

    "Total pitching angle"
    Δα::Float64

    α::Profile
    α̇::Profile
    α̈::Profile
end

function Pitchup(U₀, a, K, α₀, t₀, Δα, ramp)
    Δt = 0.5Δα/K
    p = ConstantProfile(α₀) + 2K*((ramp >> t₀) - (ramp >> (t₀ + Δt)))
    ṗ = d_dt(p)
    p̈ = d_dt(ṗ)
    Pitchup(U₀, a, K, α₀, t₀, Δα, c, ċ, c̈, p, ṗ, p̈)
end

function (p::Pitchup)(t)
    α = p.α(t)
    α̇ = p.α̇(t)
    α̈ = p.α̈(t)

    c = p.U₀*t - p.a*exp(im*α)
    ċ = p.U₀ - p.a*im*α̇*exp(im*α)
    if (t - p.t₀) > p.Δα/p.K
        c̈ = 0.0im
    else
        c̈ = p.a*exp(im*α)*(α̇^2 - im*α̈)
    end

    return c, ċ, c̈, α, α̇, α̈
end

function show(io::IO, p::Pitchup)
    print(io, "Pitch-up kinematics with rate K = $(p.K)")
end

"""
    PitchHeave <: Kinematics

Kinematics describing an oscillatory pitching and heaving (i.e. plunging) motion

# Constructors
# Fields
$(FIELDS)
"""
struct PitchHeave <: Kinematics
    "Freestream velocity"
    U₀::Float64

    "Axis of pitch rotation, relative to the plate centroid"
    a::Float64

    "Reduced frequency ``K = \\frac{\\Omega c}{2U_0}``"
    K::Float64

    "Phase of pitch (in radians)"
    ϕp::Float64

    "Phase of heave (in radians)"
    ϕh::Float64

    "Mean angle of attack"
    α₀::Float64

    "Amplitude of pitching"
    Δα::Float64

    "Amplitude of translational heaving"
    A::Float64

    Y::Profile
    Ẏ::Profile
    Ÿ::Profile

    α::Profile
    α̇::Profile
    α̈::Profile
end

function PitchHeave(U₀, a, K, ϕp, α₀, Δα, A, ϕh)
    p = A*(Sinusoid(2K) >> (ϕp/(2K)))
    ṗ = d_dt(p)
    p̈ = d_dt(ṗ)
    α = ConstantProfile(α₀) + Δα*(Sinusoid(2K) >> (ϕh/(2K)))
    α̇ = d_dt(α)
    α̈ = d_dt(α̇)
    PitchHeave(U₀, a, K, ϕp, ϕh, α₀, Δα, A, p, ṗ, p̈, α, α̇, α̈)
end

function (p::PitchHeave)(t)
    α = p.α(t)
    α̇ = p.α̇(t)
    α̈ = p.α̈(t)

    c = p.U₀*t + im*p.Y(t) - p.a*exp(im*α)
    ċ = p.U₀ + im*p.Ẏ(t) - p.a*im*α̇*exp(im*α)
    c̈ = im*p.Ÿ(t) + p.a*exp(im*α)*(α̇^2 - im*α̈)

    return c, ċ, c̈, α, α̇, α̈
end

function show(io::IO, p::PitchHeave)
    println(io, "Oscillatory pitch-heave kinematics with")
    println(io, "     Reduced frequency K = $(p.K)")
    println(io, "     Heaving amplitude A = $(p.A)")
    println(io, "     Pitching amplitude Δα = $(p.Δα)")
    println(io, "     Pitch lag ϕp = $(p.ϕp)")
    println(io, "     Heave lag ϕh = $(p.ϕh)")
end

struct Oscillation <: Kinematics
    "Angular frequency"
    Ω :: Float64

    "Amplitude x direction"
    Ax:: Float64

    "Phase in x direction."
    ϕx :: Float64

    "Amplitude y direction"
    Ay:: Float64

    "Phase in y direction."
    ϕy :: Float64

    cx::Profile
    ċx::Profile
    c̈x::Profile

    cy::Profile
    ċy::Profile
    c̈y::Profile

end

function Oscillation(Ω,Ax,ϕx,Ay,ϕy)
    Δtx = ϕx/Ω
    px = Ax*(Sinusoid(Ω) << Δtx)
    ṗx = d_dt(px)
    p̈x = d_dt(ṗx)

    Δty = ϕy/Ω
    py = Ay*(Sinusoid(Ω) << Δty)
    ṗy = d_dt(py)
    p̈y = d_dt(ṗy)
    Oscillation(Ω, Ax, ϕx, Ay, ϕy, px, ṗx, p̈x, py, ṗy, p̈y)
end

function (p::Oscillation)(t)
    α = 0.0
    α̇ = 0.0
    α̈ = 0.0

    c = ComplexF64(p.cx(t)) + im*ComplexF64(p.cy(t))
    ċ = ComplexF64(p.ċx(t)) + im*ComplexF64(p.ċy(t))
    c̈ = ComplexF64(p.c̈x(t)) + im*ComplexF64(p.c̈y(t))
    return c, ċ, c̈, α, α̇, α̈

    #return [p.ċ(t),0.0], [p.c̈(t),0.0], α̇
end

struct OscilX <: Kinematics
    "Angular frequency"
    Ω :: Float64

    "Mean velocity"
    Umean :: Float64

    "Velocity amplitude"
    Ux:: Float64

    "Velocity phase"
    ϕx :: Float64

    cx::Profile
    ċx::Profile
    c̈x::Profile

end

function OscilX(Ω,Umean,Ux,ϕx)
    Δtx = ϕx/Ω
    px = ConstantProfile(0.0)
    ṗx = ConstantProfile(Umean) + Ux*(Sinusoid(Ω) << Δtx)
    p̈x = d_dt(ṗx)

    OscilX(Ω, Umean, Ux, ϕx, px, ṗx, p̈x)
end

function (p::OscilX)(t)
    α = 0.0
    α̇ = 0.0
    α̈ = 0.0

    c = ComplexF64(p.cx(t))
    ċ = ComplexF64(p.ċx(t))
    c̈ = ComplexF64(p.c̈x(t))
    return c, ċ, c̈, α, α̇, α̈

    #return [p.ċ(t),0.0], [p.c̈(t),0.0], α̇
end

abstract type Switch end
abstract type SwitchOn <: Switch end
abstract type SwitchOff <: Switch end

"""
    SwitchedKinematics <: Kinematics

Modulates a given set of kinematics between simple on/off states. The velocity
specified by the given kinematics is toggled on/off.

# Fields
$(FIELDS)
"""
struct SwitchedKinematics{S <: Switch} <: Kinematics

    "time at which the kinematics should be turned on"
    t_on :: Float64

    "time at which the kinematics should be turned off"
    t_off :: Float64

    "kinematics to be followed in the on state"
    kin :: Kinematics

    off :: Kinematics

    SwitchedKinematics(t_on,t_off,kin) = t_on > t_off ?
            new{SwitchOn}(t_on,t_off,kin,RigidBodyMotions.Constant(0,0)) :
            new{SwitchOff}(t_on,t_off,kin,RigidBodyMotions.Constant(0,0))
end

# note that these do not introduce impulsive changes into the derivatives
(p::SwitchedKinematics{SwitchOn})(t) = t <= p.t_on ? p.off(t) : p.kin(t-p.t_on)

(p::SwitchedKinematics{SwitchOff})(t) = t <= p.t_off ? p.kin(t-p.t_on) : p.off(t)

#=
Profiles
=#

"""
    ConstantProfile(c::Number)

Create a profile consisting of a constant `c`.

# Example

```jldoctest
julia> p = RigidBodyMotions.ConstantProfile(1.0)
Constant (2.3)
```
"""
struct ConstantProfile <: Profile
    c::Number
end

function show(io::IO, p::ConstantProfile)
    print(io, "Constant ($(p.c))")
end

(p::ConstantProfile)(t) = p.c

struct DerivativeProfile{P} <: Profile
    p::P
end

function show(io::IO, ṗ::DerivativeProfile)
    print(io, "d/dt ($(ṗ.p))")
end

(ṗ::DerivativeProfile)(t) = ForwardDiff.derivative(ṗ.p, t)

"""
    d_dt(p::Profile)

Take the time derivative of `p` and return it as a new profile.

# Example

```jldoctest
julia> s = Plates.RigidBodyMotions.Sinusoid(π)
Sinusoid (ω = 3.14)

julia> s.([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 0.0
 1.0
 0.707107

julia> c = Plates.RigidBodyMotions.d_dt(s)
d/dt (Sinusoid (ω = 3.14))

julia> c.([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
  3.14159
  1.92367e-16
 -2.22144
```
"""
d_dt(p::Profile) = DerivativeProfile(p)

struct ScaledProfile{N <: Real, P <: Profile} <: Profile
    s::N
    p::P
end
function show(io::IO, p::ScaledProfile)
    print(io, "$(p.s) × ($(p.p))")
end

"""
    s::Number * p::Profile

Returns a scaled profile with `(s*p)(t) = s*p(t)`

# Example

```jldoctest
julia> s = Plates.RigidBodyMotions.Sinusoid(π)
Sinusoid (ω = 3.14)

julia> 2s
2 × (Sinusoid (ω = 3.14))

julia> (2s).([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 0.0
 2.0
 1.41421
```
"""
s::Number * p::Profile = ScaledProfile(s, p)

"""
    -(p₁::Profile, p₂::Profile)

```jldoctest
julia> s = Plates.RigidBodyMotions.Sinusoid(π)
Sinusoid (ω = 3.14)

julia> 2s
2 × (Sinusoid (ω = 3.14))

julia> (2s).([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 0.0
 2.0
 1.41421

julia> s = Plates.RigidBodyMotions.Sinusoid(π);

julia> s.([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 0.0
 1.0
 0.707107

julia> (-s).([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 -0.0
 -1.0
 -0.707107

julia> (s - s).([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 0.0
 0.0
 0.0
```
"""
-(p::Profile) = ScaledProfile(-1, p)
(p::ScaledProfile)(t) = p.s*p.p(t)

struct ShiftedProfile{N <: Real, P <: Profile} <: Profile
    Δt::N
    p::P
end
function show(io::IO, p::ShiftedProfile)
    print(io, "$(p.p) >> $(p.Δt)")
end

(p::ShiftedProfile)(t) = p.p(t - p.Δt)

"""
    p::Profile >> Δt::Number

Shift the profile in time so that `(p >> Δt)(t) = p(t - Δt)`

# Example

```jldoctest
julia> s = Plates.RigidBodyMotions.Sinusoid(π);

julia> s >> 0.5
Sinusoid (ω = 3.14) >> 0.5

julia> (s >> 0.5).([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
 -1.0
  0.0
  0.707107

julia> (s << 0.5).([0.0, 0.5, 0.75])
3-element Array{Float64,1}:
  1.0
  1.22465e-16
 -0.707107
```
"""
p::Profile >> Δt::Number = ShiftedProfile(Δt, p)
p::Profile << Δt::Number = ShiftedProfile(-Δt, p)

struct AddedProfiles{T <: Tuple} <: Profile
    ps::T
end
function show(io::IO, Σp::AddedProfiles)
    println(io, "AddedProfiles:")
    for p in Σp.ps
        println(io, "  $p")
    end
end

"""
    p₁::Profile + p₂::Profile

Add the profiles so that `(p₁ + p₂)(t) = p₁(t) + p₂(t)`.

# Examples

```jldoctest
julia> ramp₁ = Plates.RigidBodyMotions.EldredgeRamp(5)
logcosh ramp (aₛ = 5.0)

julia> ramp₂ = Plates.RigidBodyMotions.ColoniusRamp(5)
power series ramp (n = 5.0)

julia> ramp₁ + ramp₂
AddedProfiles:
  logcosh ramp (aₛ = 5.0)
  power series ramp (n = 5.0)

julia> ramp₁ + (ramp₂ + ramp₁) == ramp₁ + ramp₂ + ramp₁
true

```
"""
+(p::Profile, Σp::AddedProfiles) = AddedProfiles((p, Σp.ps...))
+(Σp::AddedProfiles, p::Profile) = AddedProfiles((Σp.ps..., p))
function +(Σp₁::AddedProfiles, Σp₂::AddedProfiles)
    AddedProfiles((Σp₁..., Σp₂...))
end

-(p₁::Profile, p₂::Profile) = p₁ + (-p₂)

+(p::Profile...) = AddedProfiles(p)

function (Σp::AddedProfiles)(t)
    f = 0.0
    for p in Σp.ps
        f += p(t)
    end
    f
end

struct MultipliedProfiles{T <: Tuple} <: Profile
    ps::T
end
function show(io::IO, Πp::MultipliedProfiles)
    println(io, "MultipliedProfiles:")
    for p in Πp.ps
        println(io, "  $p")
    end
end
*(p::Profile, Πp::MultipliedProfiles) = MultipliedProfiles((p, Πp.ps...))
*(Πp::MultipliedProfiles, p::Profile) = MultipliedProfiles((Πp.ps..., p))
function *(Πp₁::MultipliedProfiles, Πp₂::MultipliedProfiles)
    MultipliedProfiles((Πp₁..., Πp₂...))
end

function (Πp::MultipliedProfiles)(t)
    f = 1.0
    for p in Πp.ps
        f *= p(t)
    end
    f
end

struct Sinusoid <: Profile
    ω::Float64
end
(s::Sinusoid)(t) = sin(s.ω*t)
show(io::IO, s::Sinusoid) = print(io, "Sinusoid (ω = $(round(s.ω, digits=2)))")

struct EldredgeRamp <: Profile
    aₛ::Float64
end
(r::EldredgeRamp)(t) = 0.5(log(2cosh(r.aₛ*t)) + r.aₛ*t)/r.aₛ
show(io::IO, r::EldredgeRamp) = print(io, "logcosh ramp (aₛ = $(round(r.aₛ, digits=2)))")

struct ColoniusRamp <: Profile
    n::Int
end
function (r::ColoniusRamp)(t)
    Δt = t + 0.5
    if Δt ≤ 0
        0.0
    elseif Δt ≥ 1
        Δt - 0.5
    else
        f = 0.0
        for j = 0:r.n
            f += binomial(r.n + j, j)*(r.n - j + 1)*(1 - Δt)^j
        end
        f*Δt^(r.n + 2)/(2r.n + 2)
    end
end
show(io::IO, r::ColoniusRamp) = print(io, "power series ramp (n = $(round(r.n, digits=2)))")

end


text:
@testset "AbstractTensor" begin
    struct Point{dim, T} <: AbstractVec{dim, T}
        x::NTuple{dim, T}
    end
    Base.Tuple(p::Point) = p.x
    Base.getindex(p::Point, i::Int) = p.x[i]
    for T in (Float32, Float64)
        x = Vec{2, T}(1, 2)
        p = Point(Tuple(x))
        @test (@inferred p + p)::Vec{2, T} |> Tuple == (x + x).data
        @test (@inferred p - p)::Vec{2, T} |> Tuple == (x - x).data
        @test (@inferred p ⋅ p)::T == x ⋅ x
        @test (@inferred p ⊗ p)::Mat{2, 2, T} == x ⊗ x
    end

    # size/axes
    TT = typeof(rand(Mat{3, 2}))
    @test (@inferred size(TT)) == (3, 2)
    @test (@inferred size(TT, 1)) == 3
    @test (@inferred size(TT, 2)) == 2
    @test (@inferred size(TT, 3)) == 1
    @test (@inferred size(TT, 4)) == 1
    @test (@inferred axes(TT)) == (Base.OneTo(3), Base.OneTo(2))
    @test (@inferred axes(TT, 1)) == Base.OneTo(3)
    @test (@inferred axes(TT, 2)) == Base.OneTo(2)
    @test (@inferred axes(TT, 3)) == Base.OneTo(1)
    @test (@inferred axes(TT, 4)) == Base.OneTo(1)
end

@testset "Operations for array interface" begin
    @testset "getindex for slicing" begin
        v = Vec(5,4,3,2,1)
        @test @inferred(v[:]) === Vec(5,4,3,2,1)
        @test @inferred(v[SOneTo(3)]) === Vec(5,4,3)
        @test @inferred(v[SUnitRange(2,4)]) === Vec(4,3,2)
        @test @inferred(v[SVector(3,4,1)]) === Vec(3,2,5)
        A = Mat{2,3}(1,2,3,4,5,6)
        @test @inferred(A[:]) === Vec(1,2,3,4,5,6)
        @test @inferred(A[1,:]) === Vec(1,3,5)
        @test A[end,:] === Vec(2,4,6)
        @test @inferred(A[:,1]) === Vec(1,2)
        @test A[:,end] === Vec(5,6)
        @test @inferred(A[1,SVector(1,3)]) === Vec(1,5)
        # by @Tensor macro
        @test @Tensor(A[:,[1,3]]) === @Mat [1 5; 2 6]
        @test @Tensor(A[1,:]) === Vec(1,3,5)
        @test @Tensor(A[end,:]) === Vec(2,4,6)
        @test @Tensor(A[1,2:3]) === Vec(3,5)
        @test_throws Exception @Tensor(Array(A)[1,2:3])
    end
    @testset "vcat/hcat" begin
        @test @inferred(vcat(Vec(1,2,3))) === Vec(1,2,3)
        @test @inferred(vcat(Vec(1,2,3)')) === Vec(1,2,3)'
        @test @inferred(vcat(Mat{1,3}(1,2,3))) === Mat{1,3}(1,2,3)
        @test @inferred(hcat(Vec(1,2,3))) === Mat{3,1}(1,2,3)
        @test @inferred(hcat(Vec(1,2,3)')) === Vec(1,2,3)'
        @test @inferred(hcat(Mat{3,1}(1,2,3))) === Mat{3,1}(1,2,3)

        @test @inferred(vcat(Vec(1,2,3), Vec(4,5,6))) === Vec(1,2,3,4,5,6)
        @test @inferred(vcat(Vec(1,2,3), Mat{3,1}(4,5,6))) === Mat{6,1}(1,2,3,4,5,6)
        @test @inferred(vcat(Mat{1,3}(1,2,3), Mat{1,3}(4,5,6))) === @Mat [1 2 3; 4 5 6]
        @test @inferred(vcat(Mat{1,3}(1,2,3), Vec(4,5,6)')) === @Mat [1 2 3; 4 5 6]
        @test @inferred(vcat(symmetric(Mat{2,2}(1,2,2,3), :U), Vec(4,5)')) === @Mat [1 2; 2 3; 4 5]
        @test @inferred(hcat(Vec(1,2,3), Vec(4,5,6))) === Mat{3,2}(1,2,3,4,5,6)
        @test @inferred(hcat(Vec(1,2,3)', Vec(4,5,6)')) === Mat{1,6}(1,2,3,4,5,6)
        @test @inferred(hcat(Vec(1,2,3), Mat{3,2}(4,5,6,7,8,9))) === Mat{3,3}(1,2,3,4,5,6,7,8,9)
        @test @inferred(hcat(Mat{3,1}(1,2,3), Mat{3,2}(4,5,6,7,8,9))) === Mat{3,3}(1,2,3,4,5,6,7,8,9)
        @test @inferred(hcat(symmetric(Mat{2,2}(1,2,2,3), :U), Vec(4,5))) === @Mat [1 2 4; 2 3 5]

        @test @inferred(vcat(Vec(1), Vec(2), Vec(3), Mat{2,1}(4, 5))) === Mat{5,1}(1,2,3,4,5)
        @test @inferred(vcat(Vec(1), Vec(2), Vec(3), Vec(4), Vec(5))) === Vec(1,2,3,4,5)
        @test @inferred(hcat(Vec(1), Vec(2), Vec(3), Vec(4), Vec(5))) === Mat{1,5}(1,2,3,4,5)

        # special cases for vcat with vector and real
        ## start from Vec
        @test @inferred(vcat(Vec(1), 2)) === Vec(1,2)
        @test @inferred(vcat(Vec(1), 2, Vec(3))) === Vec(1,2,3)
        ## start from Real
        n = 10
        for I in 1:n
            @test @inferred(vcat(1:I..., Vec(I+1))) === Vec((1:I+1)...)
            @test @inferred(vcat(1:I..., Vec(I+1), I+2)) === Vec((1:I+2)...)
            @test @inferred(vcat(1:I..., Vec(I+1), Vec(I+2))) === Vec((1:I+2)...)
        end
        n += 1
        @test @inferred(vcat(1:n..., Vec(n+1)))::Vector == 1:n+1
        @test @inferred(vcat(1:n..., Vec(n+1), n+2))::Vector == 1:n+2
        @test @inferred(vcat(1:n..., Vec(n+1), Vec(n+2)))::Vector == 1:n+2

        @test vcat(Vec(1.0f0), Vec(1.0)) === Vec(1.0, 1.0)
        @test hcat(Vec(1.0f0), Vec(1.0)) === Mat{1,2}(1.0, 1.0)
    end
    @testset "reverse" begin
        @test @inferred(reverse(Vec(1, 2, 3))) === Vec(3, 2, 1)
        A = rand(Mat{3,4})
        if VERSION ≥ v"1.6"
            @test @inferred(reverse(A))::typeof(A) == reverse(reverse(collect(A), dims = 1), dims = 2)
        end
    end
end


text:
@testset "ReLU Last Layer" begin

    small_nnet_file = "$(@__DIR__)/../examples/networks/small_nnet.nnet"
    # small_nnet encodes the simple function 24*max(x + 1.5, 0) + 18.5
    small_nnet = read_nnet(small_nnet_file, last_layer_activation = ReLU())

    # The input set is [-0.9:0.9]
    in_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])
    in_hpoly  = convert(HPolytope, in_hyper)

    # Output region is entirely contained in this interval:
    out_superset    = Hyperrectangle(low = [30.0], high = [80.0])    # 20.0 ≤ y ≤ 90.0
    # Includes some points in the output region but not all:
    out_overlapping = Hyperrectangle(low = [-1.0], high = [50.0])    # -1.0 ≤ y ≤ 50.0

    @testset "Group 1" begin
        problem_holds    = Problem(small_nnet, in_hpoly, convert(HPolytope, out_superset))
        problem_violated = Problem(small_nnet, in_hpoly, convert(HPolytope, out_overlapping))

        for solver in [MaxSens(resolution = 0.6), ExactReach(), Ai2(), Ai2h(), Box()]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end

    end

    @testset "Group 2, 3, 4, 6" begin
        # halfspace only constraints:
        problem_holds    = Problem(small_nnet, in_hyper, HPolytope([HalfSpace([1.], 100.)]))     # y < 100.0
        problem_violated = Problem(small_nnet, in_hyper, HPolytope([HalfSpace([1.], 10.)]))      # y < 10.0

        group2 = [NSVerify(), MIPVerify(), ILP()]
        group3 = [Duality()] # hypothetically also ConvDual
        group4 = [FastLin(), FastLip()]
        group6 = [Reluplex(), Planet()]

        for solver in [group2; group3; group4; group6]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end

        # ConvDual can not handle ReLU networks at present.
        # We should ignore the result even if this particular network is trivial
        holds    = solve(ConvDual(), problem_holds)
        violated = solve(ConvDual(), problem_violated)
        @test_skip holds.status    ∈ (:holds, :unknown)
        @test_skip violated.status ∈ (:violated, :unknown)
    end

    @testset "Group 5" begin
        problem_holds    = Problem(small_nnet, in_hyper, out_superset)
        problem_violated = Problem(small_nnet, in_hyper, out_overlapping)

        for solver in [ReluVal(max_iter = 10), DLV(), Sherlock(ϵ = 0.5), BaB(), Neurify()]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end
    end

    @testset "Certify" begin
        ### Certify - only works for single hidden layer
        tiny_nnet = read_nnet("$(@__DIR__)/../examples/networks/tiny_nnet.nnet")
        solver_certify = Certify()
        in_set  = Hyperrectangle([2.0], [.5])
        out_set = HPolytope(ones(1,1), [2.5])
        problem_certify = Problem(tiny_nnet, in_set, out_set)
        @test @no_error solve(solver_certify, problem_certify)
    end
end