text:
using Test
using AdaptiveSparseGrids
using ForwardDiff
using StaticArrays
using QuadGK

@testset "Adaptive Sparse Grid Tests" begin

@testset "Basis Functions" begin
    import AdaptiveSparseGrids: m, Y, h, ϕ, Dϕ, leftchild, rightchild

    @testset "m(i) : Dimension of basis at level i" begin
        @test m(1) == 1
        @test m(2) == 2
        @test m(3) == 4
        @test m(4) == 8
        @test_throws ArgumentError  m(0)
        @test_throws MethodError    m(1.0)
        @test_throws MethodError    m(:test)
    end

    @testset "Y(i,j) : Get the coordinate of the jth node in the ith layer" begin
        @test Y(1, 1) == 0.5
        @test Y(2, 0) == 0.0
        @test Y(2, 1) == 0.5
        @test Y(2, 2) == 1.0

        @testset "Check raw numbers" begin
            for i = 3:10, j = 0:m(i)
                @test Y(i, j) == j/2^(i-1)
            end
        end

        @testset "Check Nesting" begin
            for i = 3:10, j = 0:m(i)
                !iseven(j) && continue
                @test Y(i,j) == Y(i-1, Int(j/2))
            end
        end

        @test_throws ArgumentError Y(0, 3)
        @test_throws ArgumentError Y(1, 2)
        @test_throws ArgumentError Y(1, -1)
        @test_throws MethodError   Y(1.0, 1)
        @test_throws MethodError   Y(1, 1.0)
    end

    @testset "h(i): Width of basis elements in ith level" begin
        @test h(1) == 1.0
        @test h(2) == 0.5
        @test h(3) == 0.25
        for i in 4:10
            @test h(i) == 2.0^(-i + 1)
        end

        @test_throws MethodError h(1.0)
        @test_throws ArgumentError h(-1)
    end

    @testset "ϕ(x): Hat function" begin
        @test ϕ(-2.0) == 0.0
        @test ϕ(-1.0) == 0.0
        @test ϕ(0.0)  == 1.0
        @test ϕ(1.0)  == 0.0
        @test ϕ(2.0)  == 0.0

        @test ϕ(0.5)  == 0.5
        @test ϕ(-0.5) == 0.5
        @test ϕ(0.25) == 0.75
        @test ϕ(0.75) == 0.25

        for x in LinRange(0, 1, 5)
            # Test Symmetry
            @test ϕ(x) == ϕ(-x)

            # Test Values
            if x <= 1
                @test ϕ(x) == 1 - x
            else
                @test ϕ(x) == 0.0
            end
        end

        # Technically it should work with integers
        @test typeof(ϕ(1.0)) == Float64
        @test typeof(ϕ(1))   == Int
        @test ϕ(-1)== ϕ(-1.0)
        @test ϕ(0) == ϕ(0.0)
        @test ϕ(1) == ϕ(1.0)
        @test ϕ(0) == ϕ(0.0)
        @test ϕ(2) == ϕ(2.0)
    end

    @testset "Dϕ: Derivatives" begin
        for x in LinRange(-1, 1, 5)
            @test Dϕ(x) == ForwardDiff.derivative(ϕ, x)
        end
    end

    @testset "ϕ(i,j,x): Hierarchical Hat Functions" begin
        @test ϕ(1, 1, 0.0) == 1.0
        @test ϕ(1, 1, 0.5) == 1.0
        @test ϕ(1, 1, 1.0) == 1.0

        @test ϕ(2, 0, Y(2, 0)) == 1.0
        @test ϕ(2, 0, Y(2, 1)) == 0.0
        @test ϕ(2, 2, Y(2, 1)) == 0.0
        @test ϕ(2, 2, Y(2, 2)) == 1.0

        for i = 3:10, j = 1:m(i) - 1
            if isodd(j)
                # Boundaries
                @test ϕ(i, j, Y(i, j-1))    == 0.0
                @test ϕ(i, j, Y(i, j+1))    == 0.0

                # Interior
                @test ϕ(i, j, Y(i, j))      == 1.0
                @test ϕ(i, j, Y(i+1, 2j-1)) == 0.5
                @test ϕ(i, j, Y(i+1, 2j+1)) == 0.5
                @test ϕ(i, j, Y(i+2, 4j-3)) == 0.25
                @test ϕ(i, j, Y(i+2, 4j-1)) == 0.75
                @test ϕ(i, j, Y(i+2, 4j+1)) == 0.75
                @test ϕ(i, j, Y(i+2, 4j+3)) == 0.25
            end
        end
    end

    @testset "lefchild and rightchild" begin

        @testset "Base Values" begin
            @test leftchild(1, 1)  == (2, 0)
            @test rightchild(1, 1) == (2, 2)
            @test rightchild(2, 0) == (3, 1)
            @test leftchild(2, 2)  == (3, 3)
            @test leftchild(2, 0)  == (3, -1)
            @test rightchild(2, 2) == (3, -1)
        end

        @testset "Higher Level Hierarchy" begin
            for i in 3:10, j in 1:2:m(i)
                @test ϕ(i,j, Y(leftchild(i,j)...))  == 0.5
                @test ϕ(i,j, Y(rightchild(i,j)...)) == 0.5

                @test ϕ(leftchild(i,j)..., Y(i,j))  == 0.0
                @test ϕ(rightchild(i,j)..., Y(i,j)) == 0.0

                # They do overlap in their support with their children
                for c in (leftchild(i,j), rightchild(i,j))
                    @test ϕ(i,j,  (Y(i,j) + Y(c...))/2) == 0.75
                    @test ϕ(c..., (Y(i,j) + Y(c...))/2) == 0.5
                end

            end

        end

        @testset "Base Level Errors" begin
            @test_throws ArgumentError leftchild(2, 1)
            @test_throws ArgumentError rightchild(2, 1)
            @test_throws ArgumentError leftchild(0, 1)
            @test_throws ArgumentError rightchild(0, 1)
        end

        @testset "No Floats" begin
            @test_throws MethodError leftchild(1.0, 1)
            @test_throws MethodError leftchild(1,  1.0)
            @test_throws MethodError rightchild(1.0, 1)
            @test_throws MethodError rightchild(1,  1.0)
        end

    end
end

@testset "Nodes" begin

    @testset "KTuples" begin
        import AdaptiveSparseGrids: KTuple
        @test (1,2) isa KTuple
        @test (1,2) isa KTuple{2}
        @test (1,2) isa KTuple{2, Int}
        @test (1.0,2.0) isa KTuple
        @test (1.0,2.0) isa KTuple{2}
        @test (1.0,2.0) isa KTuple{2, Float64}

        @test (1,2,3) isa KTuple
        @test (1,2,3) isa KTuple{3}
        @test (1,2,3) isa KTuple{3, Int}
        @test (1.0,2.0,3.0) isa KTuple
        @test (1.0,2.0,3.0) isa KTuple{3, Float64}

        @test ! ((1, 2.0) isa KTuple)

        @test (a = 1.0, b = 2.0) isa KTuple
        @test (a = 1.0, b = 2.0) isa KTuple{2}
        @test (a = 1.0, b = 2.0) isa KTuple{2, Float64}

        @test !( (a = 1.0, b = 2) isa KTuple )
        @test !( (a = 1.0, b = 2) isa KTuple{2} )
        @test !( (a = 1.0, b = 2) isa KTuple{2, Float64} )
    end

    @testset "getzero" begin
        import AdaptiveSparseGrids: getzero
        @test getzero((a=1, b=2))               == (a = 0,   b = 0)
        @test getzero((a=1.0, b=2.0))           == (a = 0.0, b = 0.0)
        @test getzero(typeof((a=1, b=2)))       == (a = 0,   b = 0)
        @test getzero(typeof((a=1.0, b=2.0)))   == (a = 0.0, b = 0.0)
        @test getzero((1, 2, 3))                == (0,0,0)
        @test getzero((1.0, 2.0, 3.0))          == (0.0,0.0,0.0)
        @test getzero([1, 2, 3])                == (0,0,0)
        @test getzero([1, 2.0, 3])              == (0.0,0.0,0.0)
        @test getzero(1.0)                      == (0.0,)
        @test getzero(0)                        == (0,)
    end

    @testset "Node 1D" begin
        import AdaptiveSparseGrids: Node, getx, ϕ, makeleftchild, makerightchild
        n = Node(1,1)
        @test ϕ(n, (0.0,), 1) == 1.0
        @test ϕ(n, (0.5,), 1) == 1.0
        @test ϕ(n, (1.0,), 1) == 1.0

        # Check the children
        lc = makeleftchild(n, 1)
        @test getx(lc) == @SVector [0.0]
        @test ϕ(lc, 0.0, 1) == 1.0
        @test ϕ(lc, 0, 1)   == 1.0
        @test ϕ(lc, 0.5, 1) == 0.0
        @test ϕ(lc, 0.25, 1) == 0.5

        rc = makerightchild(n, 1)
        @test getx(rc) == @SVector [1.0]
        @test ϕ(rc, 1.0, 1) == 1.0
        @test ϕ(rc, 1, 1)   == 1.0
        @test ϕ(rc, 0.5, 1) == 0.0
        @test ϕ(rc, 0.75, 1) == 0.5

        for l in 3:10, j in 1:2:m(l)
            n = Node(l, j)
            @test getx(n) == @SVector [Y(l, j)]
            @test ϕ(n, getx(n), 1) == 1.0
            @test ϕ(n, Y(leftchild(l,j)...), 1)  == 0.5
            @test ϕ(n, Y(rightchild(l,j)...), 1) == 0.5
            j > 1       && @test ϕ(n, Y(l, j-1), 1) == 0.0
            j < m(l)- 1 && @test ϕ(n, Y(l, j+1), 1) == 0.0
            j > 2       && @test ϕ(n, Y(l, j-2), 1) == 0.0
            j < m(l)-1  && @test ϕ(n, Y(l, j+2), 1) == 0.0
        end
    end

    @testset "Node 2D" begin
        import AdaptiveSparseGrids: Node, getx, ϕ, m, Y
        for l1 in 3:5, l2 in 3:5, i1 in 1:2:m(l1), i2 in 1:2:m(l2)
            n = Node((l1,l2), (i1, i2))

            # Check that we're linearly interpolating to the corners
            for d1 in -1:1, d2 in -1:1
                c = @SVector [Y(l1, i1 + d1), Y(l2, i2+d2)]
                for θ in LinRange(0,1, 5)
                    x = θ * c + (1-θ) * getx(n)
                    @test ϕ(n, x) == (1-θ)^(abs(d1) + abs(d2))
                end
            end
        end
    end
end

@testset "Scaling" begin
    import AdaptiveSparseGrids: scale, rescale

    f = AdaptiveSparseGrid([0.0,1.0], [2.0, 3.0]) do (x,)
        x^2
    end

    @test scale(f, [0.0, 1.0]) == [0.0, 0.0]
    @test scale(f, (0.0, 1.0)) == [0.0, 0.0]
    @test scale(f, @SVector [0.0, 1.0]) == [0.0, 0.0]

    @test scale(f, [1.0, 1.0]) == [0.5, 0.0]
    @test scale(f, [1.0, 2.0]) == [0.5, 0.5]
    @test scale(f, [1.0, 3.0]) == [0.5, 1.0]

    @test scale(f, [2.0, 1.0]) == [1.0, 0.0]
    @test scale(f, [2.0, 2.0]) == [1.0, 0.5]
    @test scale(f, [2.0, 3.0]) == [1.0, 1.0]

    @test rescale(f, [0.5, 0.0]) == [1.0, 1.0]
    @test rescale(f, [0.5, 0.5]) == [1.0, 2.0]
    @test rescale(f, [0.5, 1.0]) == [1.0, 3.0]

    @test rescale(f, [1.0, 0.0]) == [2.0, 1.0]
    @test rescale(f, [1.0, 0.5]) == [2.0, 2.0]
    @test rescale(f, [1.0, 1.0]) == [2.0, 3.0]
end

@testset "Interpolation Tests (1D)" begin

    for g in [x -> x^2, sin, cos, exp]
        f = AdaptiveSparseGrid([0.0], [10.0], tol=1e-8, max_depth = 20) do x
            g.(x)
        end

        for x in LinRange(0, 10, 1_000)
            @test abs(f(x) - g(x))/max(g(x),1) < 1e-8
        end
    end

    # Test the named tuple features
    f = AdaptiveSparseGrid([0.0], [10.0], tol=1e-8, max_depth = 20) do (x,)
        (sin = sin(x), cos = cos(x), x2 = x^2)
    end

    @test f(pi/2, :sin) ≈ 1.0
    @test f(pi/2, :cos) ≈ 0.0 atol=1e-8
    @test f(pi  , :sin) ≈ 0.0 atol=1e-8
    @test f(pi  , :cos) ≈ -1.0

    @test f(pi/2, :sin) == f(pi/2, 1)
    @test f(pi/2, :cos) == f(pi/2, 2)

    # Test that everything still works when you construct with Tuple Boundaries
    for g in [x -> x^2, sin, cos, exp]
        f = AdaptiveSparseGrid((0.0,), (10.0,), tol=1e-8, max_depth = 20) do x
            g.(x)
        end

        for x in LinRange(0, 10, 1_000)
            @test abs(f(x) - g(x))/max(g(x),1) < 1e-8
        end
    end

end

@testset "Integration Tests (1D)" begin

    for g in [x -> x^2, sin, cos, exp]

        f = AdaptiveIntegral([0.0], [10.0], 1, tol=1e-8, max_depth = 20) do x
            g.(x)
        end

        v = quadgk(g, 0, 10)[1]
        @test abs(f()[1] - v)/max(v,1) < 1e-8
    end
end

@testset "Interpolation Tests (2D)" begin

    h((x,y)) = sin(x) * cos(x)
    h(x...)  = h(x)

    f = AdaptiveSparseGrid(h, [0.0, 0.0], [2 * pi, 2 * pi],
                           tol=1e-6, max_depth = 20)

    for x in LinRange(0,2*pi,100), y in LinRange(0,2*pi,100)
        @test f((x,y)) ≈ h(x,y) atol = 1e-6
        @test f((x,y)) == f(x,y)
    end

    # Check that using integers works fine
    for x in LinRange(0, 2*pi, 100)
        @test f((x,1)) == f(x, 1.0)
        @test f((1,x)) == f(1.0, x)
        @test f((x,2)) == f(x, 2.0)
        @test f((2,x)) == f(2.0, x)
    end

    # Check that the max depth parameter works
    @test mapreduce(x -> x.depth, max, values(nodes(f))) <= AdaptiveSparseGrids.max_depth(f)

    # Check that Tuple Constructors Work
    f = AdaptiveSparseGrid(h, (0.0, 0.0), (2 * pi, 2 * pi),
                           tol=1e-6, max_depth = 20)

    for x in LinRange(0,2*pi,100), y in LinRange(0,2*pi,100)
        @test f((x,y)) ≈ h(x,y) atol = 1e-6
        @test f((x,y)) == f(x,y)
    end

    # Now a version from R^2 →  R^2
    h((x,y)) = (a = sin(x) * cos(y), b = x^2 + y^2)

    f = AdaptiveSparseGrid(h, [0.0, 0.0], [2 * pi, 2 * pi],
                           tol=1e-6, max_depth = 30)

    for x in LinRange(0,2*pi,100), y in LinRange(0,2*pi,100)
        @test f((x,y),:a) ≈ h(x,y).a atol = 1e-5
    end

    # Check that using integers works fine
    for x in LinRange(0, 2*pi, 3)
        @test f((x,1)) == f((x, 1.0))
    end

    for x in LinRange(0, 2*pi, 3)
        @test f((1,x)) == f((1.0, x))
    end

end

end


text:
@reexport using Distributions
import Distributions: mean, var, logpdf, pdf, cdf, maximum, minimum, insupport, quantile
export ContinuousUnivariateCopula

"""
    ContinuousUnivariateCopula(d, c0, c1, c2)

The distribution with density `c * f(x) * (c0 + c1 * x + c2 * x^2)`, where `f`
is the density of the base distribution `d` and `c` is the normalizing constant.
"""
struct ContinuousUnivariateCopula{
    DistT <: ContinuousUnivariateDistribution,
    T     <: Real
    } <: ContinuousUnivariateDistribution
    d  :: DistT
    c0 :: T
    c1 :: T
    c2 :: T
    c  :: T # normalizing constant
end

function ContinuousUnivariateCopula(
    d  :: DistT,
    c0 :: T,
    c1 :: T,
    c2 :: T) where {DistT <: ContinuousUnivariateDistribution, T <: Real}
    c  = inv(c0 + c1 * mean(d) + c2 * (var(d) + abs2(mean(d))))
    Tc = typeof(c)
    ContinuousUnivariateCopula(d, Tc(c0), Tc(c1), Tc(c2), c)
end

minimum(d::ContinuousUnivariateCopula) = minimum(d.d)
maximum(d::ContinuousUnivariateCopula) = maximum(d.d)
insupport(d::ContinuousUnivariateCopula) = insupport(d.d)

"""
    mvsk_to_absm(μ, σ², sk, kt)

Convert mean `μ`, variance `σ²`, skewness `sk`, and kurtosis `kt` to first four
moments about zero. See formula at <https://en.wikipedia.org/wiki/Central_moment#Relation_to_moments_about_the_origin>.
"""
function mvsk_to_absm(μ, σ², sk, kt)
    σ = sqrt(σ²)
    m1 = μ
    m2 = σ² + abs2(μ)
    m3 = sk * σ^3 + 3μ * m2 - 2μ^3
    m4 = kt * abs2(σ²) + 4μ * m3 - 6 * abs2(μ) * m2 + 3μ^4
    m1, m2, m3, m4
end

function mean(d::ContinuousUnivariateCopula)
    μ, σ², sk, kt = mean(d.d), var(d.d), skewness(d.d), kurtosis(d.d, false)
    m1, m2, m3, _ = mvsk_to_absm(μ, σ², sk, kt)
    d.c * (d.c0 * m1 + d.c1 * m2 + d.c2 * m3)
end

function var(d::ContinuousUnivariateCopula)
    μ, σ², sk, kt = mean(d.d), var(d.d), skewness(d.d), kurtosis(d.d, false)
    _, m2, m3, m4 = mvsk_to_absm(μ, σ², sk, kt)
    d.c * (d.c0 * m2 + d.c1 * m3 + d.c2 * m4) - abs2(mean(d))
end

function logpdf(d::ContinuousUnivariateCopula, x::Real)
    log(d.c * (d.c0 + d.c1 * x + d.c2 * abs2(x))) + logpdf(d.d, x)
end

function pdf(d::ContinuousUnivariateCopula, x::Real)
    d.c * pdf(d.d, x) * (d.c0 + d.c1 * x + d.c2 * abs2(x))
end

# this function specialized to Normal base distribution
function cdf(d::ContinuousUnivariateCopula{Normal{T},T}, x::Real) where T <: Real
    μ, σ    = d.d.μ, d.d.σ
    z       = (x - μ) / σ
    result  = (d.c0 + d.c1 * μ + d.c2 * abs2(μ)) * cdf(Normal(), z)
    result += (d.c1 + 2d.c2 * μ) * σ / sqrt(2π) *
        ifelse(z < 0, - ccdf(Chi(2), -z), cdf(Chi(2), z) - 1)
    result += d.c2 * abs2(σ) / 2 * ifelse(z < 0, ccdf(Chi(3), -z), cdf(Chi(3), z) + 1)
    result *= d.c
end

quantile(d::ContinuousUnivariateCopula, p::Real) = Distributions.quantile_newton(d, p, mean(d))

ext: jl


text:
#=
Model: CD-Player (120 variables, 2 inputs)
=#
using Reachability, MAT, SparseArrays

cdplayer(o::Pair{Symbol, <:Any}...) = cdplayer(Options(Dict{Symbol, Any}(o)))

function cdplayer(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    file = matopen(@relpath "cdplayer.mat")
    A = read(file, "A")

    # initial set
    X0 = BallInf(zeros(120), 1.0)

    # input set
    B = read(file, "B")
    U = BallInf([0.0, 0.0], 1.0)

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # prpoperty: 2*x1 -3*x2 < 450.8
    property =
        SafeStatesProperty(HalfSpace(sparsevec([1, 2], [2., -3.], 120), 450.8))

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:60] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [1],
                                  :partition => partition,
                                  :plot_vars => [0, 1],
                                  :assume_sparse => true)
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [1, 2],
                                  :partition => partition,
                                  :property => property,
                                  :assume_sparse => true)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
export Eye

"""
`Eye([domainType=Float64::Type,] dim_in::Tuple)`

`Eye([domainType=Float64::Type,] dims...)`

Create the identity operator.

```julia
julia> op = Eye(Float64,(4,))
I  ℝ^4 -> ℝ^4

julia> op = Eye(2,3,4)
I  ℝ^(2, 3, 4) -> ℝ^(2, 3, 4)

julia> op*ones(2,3,4) == ones(2,3,4)
true

```

"""
struct Eye{T, N} <: LinearOperator
	dim::NTuple{N, Integer}
end

# Constructors
###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int})
Eye(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} = Eye{DomainType,N}(DomainDim)
###

Eye(t::Type, dims::Vararg{Integer}) = Eye(t,dims)
Eye(dims::NTuple{N, Integer}) where {N} = Eye(Float64,dims)
Eye(dims::Vararg{Integer}) = Eye(Float64,dims)
Eye(x::A) where {A <: AbstractArray} = Eye(eltype(x), size(x))

# Mappings

mul!(y::AbstractArray{T, N}, L::Eye{T, N}, b::AbstractArray{T, N}) where {T, N} = y .= b
mul!(y::AbstractArray{T, N}, L::AdjointOperator{Eye{T, N}}, b::AbstractArray{T, N}) where {T, N} = mul!(y, L.A, b)

# Properties
diag(L::Eye) = 1.
diag_AcA(L::Eye) = 1.
diag_AAc(L::Eye) = 1.

domainType(L::Eye{T, N}) where {T, N} = T
codomainType(L::Eye{T, N}) where {T, N} = T

size(L::Eye) = (L.dim, L.dim)

fun_name(L::Eye) = "I"

is_eye(L::Eye) = true
is_diagonal(L::Eye) = true
is_AcA_diagonal(L::Eye) = true
is_AAc_diagonal(L::Eye) = true
is_orthogonal(L::Eye) = true
is_invertible(L::Eye) = true
is_full_row_rank(L::Eye) = true
is_full_column_rank(L::Eye) = true


text:
using MPISimulations
using Test

include("MagneticMoments.jl")
include("MagneticFields.jl")
include("Magnetization.jl")
include("TOML.jl")
include("TracerDistributions.jl")
include("Receivers.jl")
include("ExcitationFields.jl")
include("Quadrature.jl")


text:
module Day01
export part1, part2

fuel1(m) = m÷3 - 2

function fuel2(m)
    acc = fuel1(m)
    f = acc
    while true
        f = fuel1(f)
        f <= 0 && return acc

        acc += f
    end
end

function main(input, fuel)
    open(input) do f
        acc = 0
        for line in readlines(f)
            m = parse(Int, line)
            acc += fuel(m)
        end
        acc
    end
end

part1(input) = main(input, fuel1)
part2(input) = main(input, fuel2)

end # module


text:
#CircuitAnalysis: function tools
#-------------------------------------------------------------------------------

#==Basic tools
===============================================================================#

#Generate error that a function using dispatchable symbols is not supported:
function nosupport(fn::Function, symlist::Vector{Symbol}, args...; kwargs...)
	argstr = [typeof(a) for a in args]
		arglist = join(argstr, ", ")
	kvstr = ["$k=$(typeof(v))" for (k,v) in kwargs]
		kwarglist = join(kvstr, ", ")
	symstr = [string(sym) for sym in symlist]
		symlist = join(symstr, ", ")
	return ArgumentError("Unsupported: $fn{$symlist}($arglist; $kwarglist)")
end

#Last line


text:
# Modelo de flash provis򱨯rio, pois o PP tem cૣlculo de flash mas
# esta fun袯⭠o n⭠o est� disponibilizada no plugin (esta fun袯⭠o seria mais
# eficiente nos cૣlculos, n⭠o teria cૣlculos repetitivos)
type Flash
	Flash()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Fracao(Dict{Symbol,Any}(
				:Upper=>2
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}",
				:PosX=>1,
				:PosY=>0.5
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{outL}",
				:PosX=>0.4,
				:PosY=>1
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{outV}",
				:PosX=>0.4,
				:PosY=>0
			)),
			[
				:(Fin.F = FoutL.F + FoutV.F),
				:(FoutV.F = TIT * Fin.F),
				:(FoutL.T = FoutV.T),
				:(FoutL.P = FoutV.P),
				:([FoutL.S,FoutL.H,FoutV.S,FoutV.H,FoutL.T,TIT] = PP2.FlashPH(FoutL.P,Fin.H)),
			],
			[
				"","","","","",
			],
			[:PP2,],
			[:TIT,:Fin,:FoutL,:FoutV,]
		)
	end
	PP2::DanaPlugin
	TIT::Fracao
	Fin::Corrente
	FoutL::Corrente
	FoutV::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Flash
function setEquationFlow(in::Flash)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	#	[S_ad,FoutL.T] = PP2.propPH(FoutL.P,Fin.H);
	#	[FoutV.S,FoutV.H] = PP2.propPTv(FoutV.P,FoutV.T+0.1*'K'); # perturbado para evitar ir para liq.
	#	[FoutL.S,FoutL.H] = PP2.propPTl(FoutL.P,FoutL.T-0.1*'K'); # perturbado para evitar ir para vap.
	addEquation(5)
	#	TIT * (FoutV.H - FoutL.H) = Fin.H - FoutL.H;

end
function atributes(in::Flash,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/flash"
	drive!(fields,_)
	return fields
end
Flash(_::Dict{Symbol,Any})=begin
	newModel=Flash()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module BitArrayTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)

@acc function bitarrtest(A::Array{Float64,1})
    pcond = A .> 2.0
    A[pcond] = 0.1
    return A
end

@acc function bitarrtest2(A::Array{Float64,1})
    A[A.>2.0] = 0.1
    return A
end

@acc function bitarrtest3(A::Array{Float64,1},B::Array{Float64,1},C::Array{Float64,1})
    pcond = A .> C
    A[pcond] = B[pcond]
    return A
end

@acc function bitarrtest4(A::Array{Float64,1}, C::Float64)
    pcond = A .> C
    return sum(A[pcond])
end

function test1()
    return bitarrtest([1.1; 2.2; 3.2; 1.9])
end

function test2()
    return bitarrtest2([1.1; 2.2; 3.2; 1.9])
end

function test3()
    return bitarrtest3([1.1; 2.2; 3.2; 1.9], [9.0; 9.0; 10.0; 11.0], [2.0; 2.0; 2.0; 2.0])
end

function test4()
    return bitarrtest4([1.1; 2.2; 3.2; 1.9], 2.0)
end

end

using Base.Test

println("Testing BitArrays...")
@test BitArrayTest.test1() == [1.1; 0.1; 0.1; 1.9]
@test BitArrayTest.test2() == [1.1; 0.1; 0.1; 1.9]
@test BitArrayTest.test3() == [1.1; 9.0; 10.0; 1.9]
@test BitArrayTest.test4() == 2.2 + 3.2
println("Done testing BitArrays.")


text:
const default_prepad=5*48000
# gain levels set from running scaling.jl
# TODO: by default everything should output full-scale, but for
# MLS, Golay, and RPMS include an option `peaknorm` to normalize by the true
# peak. Alternatively we could just add a `peaknorm` function separately
const mls_gain = 0.38911592741900447
const golay_gain = 0.4143208008205667
const rpms_gain = 0.2316778907504581

const optiondoc_gain =
"""
- `gain::AbstractFloat`: Gain applied to the signal to reduce the chance of
  distortion. The default provides enough headroom so that the 99.99% of the
  analog stimulus is below -1dbFS. This means that if you set the gain of your
  playback system such that a full-scale sinusoid does not distort, this signal
  should not distort either.
"""

const optiondoc_prepad =
"""
- `prepad::Int`: Defaults to `$default_prepad`. Puts a period of silence at the
  beginning of the measurement that can be used to estimate the noise floor.
"""


text:
using Plots

# https://github.com/JuliaLang/julia-logo-graphics
blue, green, purple, red = :royalblue, :forestgreen, :mediumorchid3, :brown3

function archimedes!(p, n, xy=(0,0), radius=1; color=blue)

    x₀,y₀=xy
    ts = range(0, 2pi, length=100)

    plot!(p, x₀ .+ sin.(ts), y₀ .+ cos.(ts), linewidth=2)

    α = ((2π)/n)/2
    αs = (-pi/2 + α):2α:(3pi/2 + α)
    r = radius/cos(α)

    xs = x₀ .+ r*cos.(αs)
    ys = y₀ .+ r*sin.(αs)

    plot!(p, xs, ys,
          fill=true,
          fillcolor=color,
          alpha=0.4)

    r = radius
    xs = x₀ .+ r*cos.(αs)
    ys = y₀ .+ r*sin.(αs)

    plot!(p, xs, ys,
          fill=true,
          fillcolor=color,
          alpha=0.8)

    p
end

gr()
Δ = 2.75
p = plot(;xlims=(-Δ,Δ), ylims=(-Δ,Δ),
         axis=nothing,
         xaxis=false,
         yaxis=false,
         legend=false,
         padding = (0.0, 0.0),
         background_color = :transparent,
         foreground_color = :black,
         aspect_ratio=:equal)
archimedes!(p, 5, (-1.5, -1); color=red )
archimedes!(p, 8, (0,   1); color=green )
archimedes!(p, 13, (1.5, -1); color=purple )

savefig(p, "logo.png")


text:
# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# Contributed by David Campbell

# FIXME(davekong) Is there support in Julia for doing more efficient IO and
# handling of byte arrays?

const revcompdata = Dict(
   'A'=> 'T', 'a'=> 'T',
   'C'=> 'G', 'c'=> 'G',
   'G'=> 'C', 'g'=> 'C',
   'T'=> 'A', 't'=> 'A',
   'U'=> 'A', 'u'=> 'A',
   'M'=> 'K', 'm'=> 'K',
   'R'=> 'Y', 'r'=> 'Y',
   'W'=> 'W', 'w'=> 'W',
   'S'=> 'S', 's'=> 'S',
   'Y'=> 'R', 'y'=> 'R',
   'K'=> 'M', 'k'=> 'M',
   'V'=> 'B', 'v'=> 'B',
   'H'=> 'D', 'h'=> 'D',
   'D'=> 'H', 'd'=> 'H',
   'B'=> 'V', 'b'=> 'V',
   'N'=> 'N', 'n'=> 'N',
)

function print_buff(b)
    br = reverse(b)
    l = length(br)
    for i = 1:60:l
        if i+59 > l
            write(br[i:end]); println()
        else
            write(br[i:i+59]); println()
        end
    end
end

function perf_revcomp()
    infile= joinpath(SHOOTOUT_DATA_PATH, "revcomp-input.txt")
    open(infile, "r") do input
        buff = UInt8[]
        while true
            line = Vector{UInt8}(readline(input))
            if isempty(line)
    #            print_buff(buff)
                return
            elseif line[1] == UInt8('>')
    #            print_buff(buff)
                buff = UInt8[]
    #            write(line)
            else
                l = length(line)-1
                append!(buff, [UInt8(revcompdata[Char(line[i])]) for i=1:l])
            end
        end
    end
end


text:
# ## Loading modules
# using Distributed
# using DataFrames
# using CSV
#
# addprocs(6)
#
# ## Load modules in parallel
# @everywhere include("RP_Cournot/modules/MyModules.jl")
#
# @everywhere begin
#     using .MyModules
#     using Random
#     using SharedArrays
#     ## Setting seed
#     Random.seed!(66636)
# end
#
# ## Settin up directories
#
# dir = @__DIR__ ;
# datadir=dir*"/data/oildata.csv" ;
# resultsdir=dir*"/results/"
#
# ## Loading data
#
# pt, qt = load_data_select(datadir)

## Initializing array for results

mchain=SharedArray{Float64,3}(zeros(Integer(size(pt,1)/6),size(qt,2),Param.nsim))
decade=SharedArray{Float64,2}(zeros(4,2))
acc_r_t=SharedArray{Float64,1}(zeros(5))
TS, pvalue = nothing, nothing
## Testing

# Markov Chain
@time MyModules.shared_fill_chain!(mchain, pt, qt, fill_chain_sem_alt!)
# Joint test
@time TS, pvalue, acc_r_t[1] = MyModules.RP_Cournot_test_par_alt(mchain)
# Testing by decade
@time begin
    dec = [0, 16, 36, 56, 72]
    for i=1:4
        m = mchain[dec[i]+1:dec[i+1],:,:]
        decade[i,1], decade[i,2], acc_r_t[i+1] = MyModules.RP_Cournot_test_par_alt(m)
    end
end

## Saving results to csv

joint_sem_alt=DataFrame(TS=TS,
                 pvalue=pvalue,
                 stars=mystars.(pvalue))
CSV.write(resultsdir*"joint_sem_alt.csv", joint_sem_alt)

decade_sem_alt=DataFrame(TS=decade[:,1],
                 pvalue=decade[:,2],
                 stars=mystars.(decade[:,2]))
CSV.write(resultsdir*"decade_sem_alt.csv", decade_sem_alt)

acc_rates_alt=DataFrame(Ar=acc_r_t)
CSV.write(resultsdir*"accp_r_alt.csv",acc_rates_alt)

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
@testset "Helper" begin
    @testset "VarString" begin
        @testset "Serialize" begin
            want = hex2bytes("0f2f5361746f7368693a302e372e322f")
            vstr = Bitcoin.VarString("/Satoshi:0.7.2/")
            @test Bitcoin.serialize(vstr) == want
        end
    end
end


text: printstyled("SeisIO.Quake submodule\n", color=:light_green)


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################
"""
    get_cmp_id(pm, nw, i)
Retrieves the id of component i. This is a tuple if the component is a branch. Otherwise, it is a singleton.
"""
function get_cmp_id(pm, nw, i)
    if  _PMD.ref(pm, nw, :meas, i, "cmp") == :branch
        branch_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
        cmp_id = (branch_id, _PMD.ref(pm,nw,:branch, branch_id)["f_bus"], _PMD.ref(pm,nw,:branch,branch_id)["t_bus"])
    else
        cmp_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
    end
    return cmp_id
end
"""
    update_voltage_bounds!(data; v_min, v_max)
Function that allows to automatically set upper (v_max) and lower (v_min) voltage bounds for all buses.
It assumes that all bus terminals have the same bounds and that there are at most four active terminals.
"""
function update_voltage_bounds!(data::Dict; v_min::Float64=0.0, v_max::Float64=Inf)
    for (_,bus) in data["bus"]
        bus["vmin"] = [v_min, v_min, v_min, v_min]
        bus["vmax"] = [v_max, v_max, v_max, v_max]
    end
end
"""
    update_generator_bounds!(data; p_min, p_max, q_min, q_max)
Function that allows to automatically set upper (p_max, q_max) and lower (p_min, q_min)
active and reactive power bounds for all generators. It assumes that there are at most
four active phases and all have the same bounds.
"""
function update_generator_bounds!(data::Dict; p_min::Float64=0.0, p_max::Float64=Inf, q_min::Float64=-Inf, q_max::Float64=Inf)
    for (_,gen) in data["gen"]
        gen["pmin"] = [p_min, p_min, p_min]
        gen["pmax"] = [p_max, p_max, p_max]
        gen["qmin"] = [q_min, q_min, q_min]
        gen["qmax"] = [q_max, q_max, q_max]
    end
end
"""
    update_load_bounds!(data::Dict; p_min::Float64=0.0, p_max::Float64=Inf, q_min::Float64=-Inf, q_max::Float64=Inf)
Function that allows to automatically set upper (p_max, q_max) and lower (p_min, q_min) active and reactive power bounds for all loads.
It assumes that there are at most four active phases and all have the same bounds.
"""
function update_load_bounds!(data::Dict; p_min::Float64=0.0, p_max::Float64=Inf, q_min::Float64=-Inf, q_max::Float64=Inf)
    for (_,load) in data["load"]
        load["pmin"] = [p_min, p_min, p_min]
        load["pmax"] = [p_max, p_max, p_max]
        load["qmin"] = [q_min, q_min, q_min]
        load["qmax"] = [q_max, q_max, q_max]
    end
end
"""
    update_all_bounds!(data; v_min, v_max, pg_min, pg_max, qg_min, qg_max, pd_min, pd_max, qd_min, qd_max)
Function that combines `update_voltage_bounds!`, `update_generator_bounds!` and `update_load_bounds!` and assigns bounds to all bus voltages and load and generator power.
"""
function update_all_bounds!(data::Dict; v_min::Float64=0.0, v_max::Float64=Inf, pg_min::Float64=0.0, pg_max::Float64=Inf, qg_min::Float64=-Inf, qg_max::Float64=Inf, pd_min::Float64=0.0, pd_max::Float64=Inf, qd_min::Float64=-Inf, qd_max::Float64=Inf)
    update_voltage_bounds!(data; v_min=v_min, v_max=v_max)
    update_generator_bounds!(data; p_min=pg_min, p_max=pg_max, q_min=qg_min, q_max=qg_max)
    update_load_bounds!(data; p_min=pd_min, p_max=pd_max, q_min=qd_min, q_max=qd_max)
end
"""
    assign_basic_individual_criteria!(data::Dict; chosen_criterion::String="rwlav")
Basic function that assigns individual criteria to measurements in data["meas"].
For each measurement, if the distribution type of at least one of the phases is normal, the criterion defaults to the chosen_criterion.
Otherwise, it is assigned the 'mle' criterion.
The function takes as input either a single measurement dictionary, e.g., data["meas"]["1"]
or the full MATHEMATICAL data model.
"""
function assign_basic_individual_criteria!(data::Dict; chosen_criterion::String="rwlav")
    if haskey(data, "meas")
        for (_, meas) in data["meas"]
            dst_type = [typeof(i) for i in meas["dst"]]
            if any(x->x==Distributions.Normal{Float64}, dst_type)
                meas["crit"] = chosen_criterion
            else
                meas["crit"] = "mle"
            end
        end
    elseif haskey(data, "dst")
        dst_type = [typeof(i) for i in data["dst"]]
        if any(x->x==Distributions.Normal{Float64}, dst_type)
            data["crit"] = chosen_criterion
        else
            data["crit"] = "mle"
        end
    else
        error("Unrecognized data input.")
    end
end
"""
    assign_residual_ub!(data::Dict; chosen_upper_bound::Float64=100, rescale::Bool=false)
Basic function that assigns upper bounds to the residual variables, by adding a `"res_max"` entry to the measurement dictionary.
The function takes as input either a single measurement dictionary, e.g., `data["meas"]["1"]` or the full measurement dictionary, `data["meas]`.
    # Arguments
    - data: `ENGINEERING` data model of the feeder, or dictionary corresponding to a single measurement
    - chosen_upper_bound: finite upper bound, defaults to 100 if no indication provided
    - rescale: `false` by default. If `true`, the `chosen_upper_bound` is multiplied by the value of `data["se_settings"]["rescaler"]`
               and their product is used as the upper bound. Otherwise,
"""
function assign_residual_ub!(data::Dict; chosen_upper_bound::Float64=100.0, rescale::Bool=false)
    upp_bound = rescale ? chosen_upper_bound*data["se_settings"]["rescaler"] : chosen_upper_bound
    if haskey(data, "meas")
        for (_, meas) in data["meas"]
            meas["res_max"] = upp_bound
        end
    elseif haskey(data, "dst")
        data["res_max"] = upp_bound
    else
        error("Unrecognized data input.")
    end
end
"""
    vm_to_w_conversion!(data::Dict)

This function should be called after measurements are added to the data dictionary. It converts voltage magnitude
    measurements into their square, so :vm is transformed into :w. It is useful when using the LinDist3Flow or SDP formulation.
    The conversion is exact if applied to a Normal distribution, while does not necessarily apply to other distributions.
    In the SDP case, :vm is currently not supported as input measurement, so this is necessary.
    In the LinDist3Flow it allows to ignore the square vm conversion constraint.
"""
function vm_to_w_conversion!(data::Dict)
    for (m,meas) in data["meas"]
        if meas["var"] == :vm
            for c in 1:length(meas["dst"])
                if meas["dst"][c] != 0.0
                    @assert (isa(meas["dst"][c], Distributions.Normal{Float64})) "vm_to_w conversion only available for the Normal distribution"
                    current_μ = _DST.mean(meas["dst"][c])
                    current_σ = _DST.std(meas["dst"][c])
                    data["meas"][m]["dst"][c] = _DST.Normal(current_μ^2, current_σ)
                end
            end
            data["meas"][m]["var"] = :w
        end
    end
end
"""
    assign_load_pseudo_measurement_info!(data::Dict, pseudo_load_list::Array, cluster_list::Array; time_step::Int64=1, day::Int64=1)

This function is a helper function to associate pseudo measurement information to a list of loads.
    #Arguments:
    - data: `ENGINEERING` data model of the feeder, or dictionary corresponding to a single measurement
    - pseudo_load_list: list of loads that are described by pseudo measurements
    - cluster_list: list of clusters to be associated one-on-one with the pseudo_load_list
    - time_step: time step to extract the probability distribuion function for, if applicable
    - day: day to extract the probability distribuion function for, if applicable
"""
function assign_load_pseudo_measurement_info!(data::Dict, pseudo_load_list::Array, cluster_list::Array; time_step::Int64=1, day::Int64=1)
    for idx in 1:length(pseudo_load_list)
        data["load"]["$(pseudo_load_list[idx])"]["pseudo"] = Dict{String, Any}("day" => day,
                                                                               "time_step" => time_step,
                                                                               "cluster" => cluster_list[idx]
                                                                               )
    end
end
"""
    assign_unique_individual_criterion!(data::Dict)
    - data: `MATHEMATICAL` data model of the network
    Assigns the criterion in data["se_settings"]["criterion"] to all individual measurements.
"""
function assign_unique_individual_criterion!(data::Dict)
    for (m, meas) in data["meas"]
        meas["crit"] = data["se_settings"]["criterion"]
    end
end
"""
    reduce_single_phase_loadbuses!(data::Dict; exclude = [])
Reduces the dimensions of voltage variables for load buses in which the connected load(s) only have one active phase.

    - data: `MATHEMATICAL` data model of the network
    - exclude: buses to which this transformation is not applied
"""
function reduce_single_phase_loadbuses!(data::Dict; exclude=[])
    load_info = get_load_info(data, exclude)
    already_checked = []
    for (bus_id,load_id,nr_conn_loads) in load_info
        load_connections = data["load"][load_id]["connections"]
        if nr_conn_loads == 1 && length(load_connections) == 1
            perform_dimension_reduction!(data, bus_id, load_id, load_connections)
        elseif nr_conn_loads > 1 && bus_id ∉ already_checked
            push!(already_checked, bus_id)
            loads_at_bus_id = [load_id for (x, load_id, z) in load_info if x == bus_id]
            used_connections = unique([data["load"][lid]["connections"] for lid in loads_at_bus_id])
            if length(used_connections) == 1 perform_dimension_reduction!(data, bus_id, load_id, load_connections) end
        end
    end
end
"reduces the dimension of bus terminals and branches f_ and t_connections to match those of the connected load(s)"
function perform_dimension_reduction!(data::Dict, bus_id::Int64, load_id, load_connections)

    conn_branches = find_branch_t_bus(data["branch"], bus_id)

    idx = load_connections[1]

    if !isempty(conn_branches)
        data["bus"]["$bus_id"]["terminals"] = load_connections
        data["bus"]["$bus_id"]["grounded"] = data["bus"]["$bus_id"]["grounded"][load_connections]
        for br_id in conn_branches
            branch = data["branch"][br_id]
            for (key,value) in branch
                if isa(value, Matrix)
                    branch[key] = reshape([value[idx,idx]], 1, 1)
                elseif isa(value, Vector)
                    branch[key] = load_connections
                end
            end
        end
    else # let's assume for simplicity of calculations of the d.o.f. that if a user is connected to a "f_bus", i.e., three-phase, the user is three phase too
        for (key, value) in data["load"][load_id]
            if isa(value, Vector) data["load"][load_id][key] = fill(value[1]/3, (3,)) end
        end
        data["load"][load_id]["connections"] = [1,2,3]
    end
end
"returns a tuple with all loads' information. Every load is assigned a tuple with the following content: (bus the load is connected to, the load index, total number of loads connected to the same bus) "
function get_load_info(data::Dict, exclude=[])
    load_buses = []
    load_idx = []
    for (l, load) in data["load"]
        if load["load_bus"] ∉ exclude
            push!(load_buses, load["load_bus"])
            push!(load_idx, l)
        end
    end
    loads_per_bus = [count(x->(x == i), load_buses) for i in load_buses]
    return zip(load_buses, load_idx, loads_per_bus)
end
"find the branches that have the load bus as t_bus. These will be the \"extremities\" of connection cables, so their number of
connections equals that of the load. If the load is directly connected to the main three-phase cable, this is not the case."
function find_branch_t_bus(branches, bus_id)
    conn_branches = []
    for (b, branch) in branches
        if branch["t_bus"] == bus_id push!(conn_branches, b) end
    end
    return conn_branches
end
"""
    function get_active_connections(pm::_PMD.AbstractUnbalancedPowerModel, nw::Int, cmp_type::Symbol, cmp_id::Int)
Returns the list of terminals, connections or t_ and f_connections, depending on the type of the component.
"""
function get_active_connections(pm::_PMD.AbstractUnbalancedPowerModel, nw::Int, cmp_type::Symbol, cmp_id::Int)
    if cmp_type == :bus
       active_conn = _PMD.ref(pm, nw, :bus, cmp_id)["terminals"]
   elseif cmp_type ∈ [:gen, :load]
       active_conn = _PMD.ref(pm, nw, cmp_type, cmp_id)["connections"]
   elseif cmp_type == :branch
       active_conn = intersect(_PMD.ref(pm, nw, :branch, cmp_id)["f_connections"], _PMD.ref(pm, nw, :branch, cmp_id)["t_connections"])
   else
       error("Measurements for component of type $cmp_type are not supported")
   end
   return active_conn
end
"""
    function add_zib_virtual_meas!(data::Dict, σ::Float64; exclude::Array=[])
Adds zero-injection buses measurements to the "meas" dictionary if they are not present.
Any bus_id given within `exclude` will not be assigned a virtual measurement.
The virtual measurements added are :pd = :qd = 0, with high accuracy σ, for each phase
of the bus these virtual loads are connected to.
"""
function add_zib_virtual_meas!(data::Dict, σ::Float64; exclude::Array=[])
    non_zibs = unique(vcat([load["load_bus"] for (_, load) in data["load"]], [gen["gen_bus"] for (_, gen) in data["gen"]]))
    zibs = [b for (b,_) in data["bus"] if parse(Int64, b) ∉ non_zibs]
    measured_buses = [meas["cmp_id"] for (_, meas) in data["meas"] if (meas["cmp"] == :bus)] # to check if the zib virtual measurement is already there
    highest_meas_idx = maximum(parse.(Int64, collect(keys(data["meas"]))))
    highest_load_idx = maximum(parse.(Int64, collect(keys(data["load"]))))
    for zib in zibs
        if parse(Int64, zib) ∉ measured_buses && parse(Int64, zib) ∉ exclude
            phases = length(data["bus"][zib]["terminals"])
            data["load"]["$(highest_load_idx+1)"] = Dict("load_bus" => parse(Int64, zib), "note"=>"virtual load")
            data["meas"]["$(highest_meas_idx+1)"] = Dict( "cmp_id" => highest_load_idx+1,
                                                          "cmp" => :load,
                                                          "var" => :pd,
                                                          "dst" => fill(_DST.Normal(0.0, σ), (phases,))
                                                        )
            data["meas"]["$(highest_meas_idx+2)"] = Dict( "cmp_id" => highest_load_idx+1,
                                                          "cmp" => :load,
                                                          "var" => :qd,
                                                          "dst" => fill(_DST.Normal(0.0, σ), (phases,))
                                                        )
            highest_meas_idx+=2
            highest_load_idx+=1
        end
    end
end
"""
finds the set of all the buses in `data` which are adjacent to `bus_idx`.
"""
function adjacent_branch_and_buses(bus_idx::Int64, data::Dict)
    conn_branches = [b for (b, branch) in data["branch"] if (branch["f_bus"] == bus_idx || branch["t_bus"] == bus_idx)]
    adj_buses = []
    for br in conn_branches
        push!(adj_buses, data["branch"][br]["f_bus"])
        push!(adj_buses, data["branch"][br]["t_bus"])
    end
    return conn_branches, (filter(x->x!=bus_idx, adj_buses) |> unique)
end
"""
Helper function to add a measurement. Currently can assign Gaussian measurements only.
"""
function add_measurement!(data::Dict, var::Symbol, cmp::Symbol, cmp_id::Int64, measured_value::Array, σ::Array; id::String="", crit=missing)
    m = isempty(id) ? "$(maximum([parse(Int64, m) for (m,meas) in data["meas"]])+1)" : id

    data["meas"][m] = Dict( "var" => var,
                            "crit" => crit,
                            "cmp" => cmp,
                            "dst" => [_DST.Normal(mu, sigma) for(mu, sigma) in zip(measured_value, σ)],
                            "cmp_id"=> cmp_id
                            )
end


text:
"""
    ConvDual

ConvDual uses convex relaxation to compute over-approximated bounds for a network

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hypercube
3. Output: halfspace

# Return
`BasicResult`

# Method
Convex relaxation with duality.

# Property
Sound but not complete.

# Reference
[E. Wong and J. Z. Kolter, "Provable Defenses against Adversarial Examples via the Convex Outer Adversarial Polytope,"
*ArXiv Preprint ArXiv:1711.00851*, 2017.](https://arxiv.org/abs/1711.00851)

[https://github.com/locuslab/convex_adversarial](https://github.com/locuslab/convex_adversarial)
"""
struct ConvDual end

function solve(solver::ConvDual, problem::Problem)
    o = dual_value(solver, problem.network, problem.input, problem.output)
    # Check if the lower bound satisfies the constraint
    if o >= 0.0
        return BasicResult(:SAT)
    end
    return BasicResult(:Unknown)
end

# compute lower bound of the dual problem.
function dual_value(solver::ConvDual, network::Network, input::Hyperrectangle{N}, output::HPolytope{N}) where N

    @assert all(iszero.(input.radius .- input.radius[1])) "input.radius must be uniform. Got $(input.radius)"

    layers = network.layers
    L, U  = get_bounds(network, input.center, input.radius[1])
    v0, d = tosimplehrep(output)

    v = vec(v0)
    o = d[1]

    for i in reverse(1:length(layers))
        o -= v'*layers[i].bias
        v = layers[i].weights'*v
        if i>1
            o += backprop!(v, U[i-1], L[i-1])
        end
    end
    o -= input.center' * v + input.radius[1] * sum(abs.(v))
    return o
end

#=
modifies v and returns o
=#
function backprop!(v::Vector{Float64}, u::Vector{Float64}, l::Vector{Float64})
    o = 0.0
    for j in 1:length(v)
        val = relaxed_ReLU(l[j], u[j])
        if val < 1.0 # if val is 1, it means ReLU result is identity so do not update (NOTE is that the right reasoning?)
            v[j] = v[j] * val
            o += v[j] * l[j]
        end
    end
    return o
end

# Forward_network and forward_layer:
# This step is similar to reachability method
function get_bounds(nnet::Network, input::Vector{Float64}, ϵ::Float64)
    layers  = nnet.layers
    n_layer = length(layers)

    l = Vector{Vector{Float64}}()
    u = Vector{Vector{Float64}}()
    γ = Vector{Vector{Float64}}()
    μ = Vector{Vector{Vector{Float64}}}()

    v1 = layers[1].weights'
    push!(γ, layers[1].bias)
    # Bounds for the first layer
    l1, u1 = input_layer_bounds(layers[1], input, ϵ)
    push!(l, l1)
    push!(u, u1)

    for i in 2:n_layer
        n_input  = length(layers[i-1].bias)
        n_output = length(layers[i].bias)

        input_ReLU = relaxed_ReLU.(last(l), last(u))
        D = Diagonal(input_ReLU)   # a matrix whose diagonal values are the relaxed_ReLU values (maybe should be sparse?)

        # Propagate existing terms
        WD = layers[i].weights*D
        v1 = v1 * WD' # TODO CHECK
        map!(g -> WD*g,   γ, γ)
        for M in μ
            map!(m -> WD*m,   M, M)
        end
        # New terms
        push!(γ, layers[i].bias)
        push!(μ, new_μ(n_input, n_output, input_ReLU, WD))

        # Compute bounds
        ψ = v1' * input + sum(γ)
        eps_v1_sum = ϵ * vec(sum(abs, v1, dims = 1))
        neg, pos = all_neg_pos_sums(input_ReLU, l, μ, n_output)
        push!(l,  ψ - eps_v1_sum + neg )
        push!(u,  ψ + eps_v1_sum - pos )
    end

    return l, u
end

# TODO rename function and inputs
function all_neg_pos_sums(slopes, l, μ, n_output)
    # n_output = length(last(l))
    neg = zeros(n_output)
    pos = zeros(n_output)
    # Need to debug
    for (i, ℓ) in enumerate(l)                # ℓ::Vector{Float64}
        for (j, M) in enumerate(μ[i])         # M::Vector{Float64}
            if 0 < slopes[j] < 1              # if in the triangle region of relaxed ReLU
                #posind = M .> 0
                neg .+= ℓ[j] * min.(M, 0) #-M .* !posind  # multiply by boolean to set the undesired values to 0.0
                pos .+= ℓ[j] * max.(M, 0) #M .* posind
            end
        end
    end
    return neg, pos
end

function input_layer_bounds(input_layer, input, ϵ)
    W, b = input_layer.weights, input_layer.bias

    out1 = vec(W * input + b)
    Δ    = ϵ * vec(sum(abs, W, dims = 2))

    l = out1 - Δ
    u = out1 + Δ
    return l, u
end

function new_μ(n_input, n_output, input_ReLU, WD)
    sub_μ = Vector{Vector{Float64}}(undef, n_input)
    for j in 1:n_input
        if 0 < input_ReLU[j] < 1 # negative region  ## TODO CONFIRM. Previously input_ReLU[j] == 0
            sub_μ[j] = WD[:, j] # TODO CONFIRM
        else
            sub_μ[j] = zeros(n_output)
        end
    end
    return sub_μ
end

function relaxed_ReLU(l::Float64, u::Float64)
    u <= 0.0 && return 0.0
    l >= 0.0 && return 1.0
    return u / (u - l)
end


text:
include("../functions/analysisPULSE.jl")

####################################################################
step_max_perturb_Xsensor=step_max_perturb÷(Ts_Xsensor/Ts_MTiXsens)

# Energysignal_Xsensor_back,
# rms_Xsensor_back,
# vdv_Xsensor_back,
# crest_factor_Xsensor_back=analysisPULSE_All(
#                         Xsensor_reduce_sync_perturb,
#                         filenameXsensor,
#                         Ts_Xsensor,
#                         step_max_perturb_Xsensor,
#                         10)
#
# plotAnalysisPULSE_All(
#             Energysignal_Xsensor_back,
#             rms_Xsensor_back,
#             vdv_Xsensor_back,
#             crest_factor_Xsensor_back,
#             filenameXsensor,
#             "Xsensor_back",
#             ["Planche 1" "Planche 2"],
#             [1 1 1 2 2 2 3 3 3]',
#             [:blue :red],
#             [:circle :utriangle])

Energysignal_Xsensor_seat,
rms_Xsensor_seat,
vdv_Xsensor_seat,
crest_factor_Xsensor_seat,
null_,
null_,
null_=analysisPULSE_All(
                    Xsensor_reduce_sync_perturb,
                    filenameXsensor,
                    Ts_Xsensor,
                    step_max_perturb_Xsensor,
                    4)

null_,
null_,
null_,
null_,
std_Xsensor_seat_x,
min_Xsensor_seat_x,
max_Xsensor_seat_x,
range_Xsensor_seat_x=analysisPULSE_All(
                    Xsensor_reduce_sync_perturb,
                    filenameXsensor,
                    Ts_Xsensor,
                    step_max_perturb_Xsensor,
                    5)

null_,
null_,
null_,
null_,
std_Xsensor_seat_y,
min_Xsensor_seat_y,
max_Xsensor_seat_y,
range_Xsensor_seat_y=analysisPULSE_All(
                    Xsensor_reduce_sync_perturb,
                    filenameXsensor,
                    Ts_Xsensor,
                    step_max_perturb_Xsensor,
                    6)

min_Xsensor_seat=Dict()
max_Xsensor_seat=Dict()
range_Xsensor_seat=Dict()
std_Xsensor_seat=Dict()
for k in collect(keys(min_Xsensor_seat_x))
    std_Xsensor_seat[k]=[std_Xsensor_seat_x[k] std_Xsensor_seat_y[k]]
    min_Xsensor_seat[k]=[min_Xsensor_seat_x[k] min_Xsensor_seat_y[k]]
    max_Xsensor_seat[k]=[max_Xsensor_seat_x[k] max_Xsensor_seat_y[k]]
    range_Xsensor_seat[k]=[range_Xsensor_seat_x[k] range_Xsensor_seat_y[k]]
end

# plotAnalysisPULSE_All(
#         Energysignal_Xsensor_seat,
#         rms_Xsensor_seat,
#         vdv_Xsensor_seat,
#         crest_factor_Xsensor_seat,
#         filenameXsensor[1:end-1,:],
#         "Xsensor_seat",
#         ["Planche 1"],
#         [ones(1,10) ones(1,10).*3]',
#         [:blue],
#         [:circle])

# using JLD
#
# jldopen("results/PULSE_analysis_gel.jld", "r+") do file
#     write(file,
#     "filenameXsensor",filenameXsensor)
#     write(file,
#     "step_max_perturb_Xsensor",step_max_perturb_Xsensor)
#     write(file,
#     "Energysignal_Xsensor_seat",Energysignal_Xsensor_seat)
#     write(file,
#     "rms_Xsensor_seat",rms_Xsensor_seat)
#     write(file,
#     "vdv_Xsensor_seat",vdv_Xsensor_seat)
#     write(file,
#     "crest_factor_Xsensor_seat",crest_factor_Xsensor_seat )
#     write(file,
#     "std_Xsensor_seat",std_Xsensor_seat)
# end

## #
toto=load("results/PULSE_analysis_gel.jld")

toto["filenameXsensor"]=filenameXsensor
toto["step_max_perturb_Xsensor"]=step_max_perturb_Xsensor
toto["Energysignal_Xsensor_seat"]=Energysignal_Xsensor_seat
toto["rms_Xsensor_seat"]=rms_Xsensor_seat
toto["vdv_Xsensor_seat"]=vdv_Xsensor_seat
toto["crest_factor_Xsensor_seat"]=crest_factor_Xsensor_seat
toto["std_Xsensor_seat"]=std_Xsensor_seat
toto["min_Xsensor_seat"]=min_Xsensor_seat
toto["max_Xsensor_seat"]=max_Xsensor_seat
toto["range_Xsensor_seat"]=range_Xsensor_seat

jldopen("results/PULSE_analysis_gel.jld", "w") do file
    toto_key=collect(keys(toto))
    for k in toto_key
        write(file,k,toto[k])
    end
end

toto=load("results/PULSE_analysis_gel.jld")

# ## #
# for k in collect(keys(titi))
#     print(k)
#     print(" : ")
#     println(haskey(toto,k))
# end


text:
#FileIO2 MIME support.
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
#TODO: Not sure if mime type strings are correct.
const MIMEASSOCIATIONS = Tuple{String,Symbol}[
	("image/bmp", :BMPFmt),
	("image/png", :PNGFmt),
	("image/gif", :GIFFmt),
	("image/jpeg", :JPEGFmt),
	("image/tiff", :TIFFFmt),

	("image/svg+xml", :SVGFmt),
	("image/cgm", :CGMFmt),
	("image/eps", :EPSFmt),
	("image/emf", :EMFFmt),
]

for (mtstr, fmt) in MIMEASSOCIATIONS; @eval begin #CODEGEN-----------------------

Base.show(stream, ::MIME{Symbol($mtstr)}, f::File{$fmt}) = write(stream, read(f.path, String))

end; end #CODEGEN---------------------------------------------------------------

#Last line


text:
## Functions to serve metadata about the PkgServer itself

function get_pkgserver_version()
    # Get PkgServer.jl's version and git sha
    version = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
    try
        repo = LibGit2.GitRepo(dirname(@__DIR__))
        gitsha = string(LibGit2.GitHash(LibGit2.GitCommit(repo, "HEAD")))
        return "$(version)-$(gitsha)"
    catch
        return "$(version)"
    end
end
# We store this once, since it's not exactly going to change
const pkgserver_version = get_pkgserver_version()

function get_num_hashnamed_files(dir)
    # If this directory doesn't exist, then we haven't cached anything!
    if !isdir(dir)
        return 0
    end

    num_files = 0
    for f in readdir(dir)
        # only consider things with a name length that is exactly a tree-hash (filters out in-progress downloads, etc...)
        # TODO: if we change hashes (e.g. SHA256), add that hash length here as well
        if length(f) != 40
            continue
        end

        # Only consider actual files
        f_path = joinpath(dir, f)
        if !isfile(f_path)
            continue
        end

        num_files += 1
    end
    return num_files
end

function get_num_pkgs_cached(cache_dir)
    # If we don't even have a `package` directory, then don't worry about it
    if !isdir(joinpath(cache_dir, "package"))
        return 0
    end

    num_pkgs_cached = 0
    for d in readdir(joinpath(cache_dir, "package"))
        # Only consider things with a name that is exactly a UUID (filters out in-progress downloads, etc...)
        if length(d) != 36
            continue
        end

        # Only consider directories
        pkg_dir = joinpath(cache_dir, "package", d)
        if !isdir(pkg_dir)
            continue
        end

        num_pkgs_cached += get_num_hashnamed_files(pkg_dir)
    end
    return num_pkgs_cached
end

function get_num_artifacts_cached(cache_dir)
    num_artifacts_cached = get_num_hashnamed_files(joinpath(cache_dir, "artifact"))
end

function collect_meta(config::Config)
    # Get a count of how many things we've cached
    cache_dir = config.cache_dir
    num_pkgs = get_num_pkgs_cached(cache_dir)
    num_artifacts = get_num_artifacts_cached(cache_dir)

    metadata = Dict(
        "pkgserver_version" => pkgserver_version,
        "julia_version" => string(VERSION),
        "packages_cached" => num_pkgs,
        "artifacts_cached" => num_artifacts,
    )
    return metadata
end

function serve_meta(config, http::HTTP.Stream)
    metadata = collect_meta(config)

    # We serve a JSON representation of some metadata about this PkgServer.
    # We generate an extremely simple subset of JSON here.
    metadata_json = string(
        "{",
        join(["$(repr(k)): $(repr(v))" for (k, v) in metadata], ","),
        "}",
    )

    HTTP.setheader(http, "Content-Length" => string(length(metadata_json)))
    HTTP.setheader(http, "Content-Type" => "application/json")
    startwrite(http)
    write(http, metadata_json)
end


text:
## Flow algorithms: maxflow, multirouteflow, connectivity

# Classical max-flow algorithm
function edmondsKarp(g::Network)
  flows = zeros(length(g.links))
  flow = 0.0

  while true
    u = g.source
    p, x = bfs(g,:Path,flows)
    if x == 0
      break
    end
    for v in p[2:end]
      i = findArc(g,u,v)
      rev = (g.links[i]).sym
      flows[i] = flows[i] + x
      flows[rev] = flows[rev] - x
      u = v
    end
    flow = flow + x
  end
  return (flow,flows)
end

#Connectivity
function connectivity(g::Network)
  function unweight(a::Arc)
    if a.cap == 0.
      return a
    else
      return Arc(a.sym,a.head,1.)
    end
  end
  arcs = map(unweight,g.links)
  gaux = Network(g.source,g.target,g.tails,arcs)
  λ, flows = edmondsKarp(gaux)
  return λ
end

# Multirouteflow
function kishimoto(g0::Network,k::Int)
  #Init
  i = 1
  flow,flows = edmondsKarp(g0)
  p = flow / k
  function restrict(a::Arc)
    Arc(a.sym,a.head,min(p,a.cap))
  end
  arcs = map(restrict,g0.links)
  g = Network(g0.source,g0.target,g0.tails,arcs)
  flow,flows = edmondsKarp(g)

  #Loop
  while flow < k * p
    p = (flow - i * p)/(k - i)
    i = i + 1
    g.links = map(restrict,g0.links)
    flow,flows = edmondsKarp(g)
  end
  #End
  return flow,flows
end


text:
# test various model functions

# word vectors
model = wordvectors("bigvecs.txt", Float32)
println(model)
rm("bigvecs.txt")
rm("vocab.txt")

modelbin = wordvectors("bigvecs.bin", kind = :binary)
rm("bigvecs.bin")

#@test_throws ArgumentError wordvectors("bigvecs.txt", kind = :unkown)

len_vecs, num_words = size(model)
wordvecs = model.vectors

len_vecs_bin, num_words_bin = size(modelbin)
@test size(wordvecs) == (len_vecs, num_words)
@test len_vecs == len_vecs_bin

words = vocabulary(model)
word1 = words[rand(1:end)]
word2 = words[rand(1:end)]
word3 = words[rand(1:end)]

@test in_vocabulary(model, word1)

n = rand(1:100)

indxs, mes = cosine(model, word1, n)
@test words[indxs] == cosine_similar_words(model, word1, n)
w4_indx = indxs[rand(1:end)]
loc = findall(in(w4_indx), indxs)
word4 = words[w4_indx]
@test index(model, word4) == w4_indx

s = similarity(model, word1, word4)
@test mes[loc[1]] ≈ s

inx, mes = analogy(model, [word1, word2], [word3], n)
@test words[inx] == analogy_words(model, [word1, word2], [word3], n)

# word clusters
model = wordclusters("bigclus.txt")
println(model)
rm("bigclus.txt")
words = vocabulary(model)
word1 = words[rand(1:end)]
@test in_vocabulary(model, word1)
idx = index(model, word1)
@test words[idx] == word1

clus_num = get_cluster(model, word1)
clus_words = get_words(model, clus_num)
@test word1 in clus_words

clus = clusters(model)
clu = clus[rand(1:end)]
get_words(model, clu)

println("model passed test...")


text:
@testset "combinedims" begin
    m = [1 5 9; 2 6 10; 3 7 11; 4 8 12]
    vv = [[1,2,3,4],[5,6,7,8],[9,10,11,12]]

    @testset for combine in (combinedims, combinedimsview)
        @test @inferred(size(combine(vv))) == size(m)
        @test @inferred(axes(combine(vv))) == axes(m)
        @test @inferred(combine(vv)) == m
        @test @inferred(combine(vv)[2, 1]) == m[2, 1]

        @test @inferred(size(combine(vv, 2))) == size(m)
        @test @inferred(axes(combine(vv, 2))) == axes(m)
        @test @inferred(combine(vv, 2)) == m
        @test @inferred(combine(vv, 2)[2, 1]) == m[2, 1]
        @test @inferred(combine(vv, 1)) == permutedims(m)

        A = rand(2, 3, 5, 7)
        @test @inferred(combine(splitdims(A))) == A
        @testset for d in [1:ndims(A); (1, 2); (1, 3); (2, 3); (2, 4); (1, 2, 3); (1, 2, 4); (2, 3, 4); (1, 2, 3, 4)]
            @test @inferred(combine(splitdims(A, d), d)) == A
            @test @inferred(combine(splitdims(A, d), d)[1, 2, 3, 4]) == A[1, 2, 3, 4]
        end
    end

    c_copy = combinedims(vv)
    c_view = combinedimsview(vv)
    vv[1][3] = 100
    @test c_copy[3, 1] == 3
    @test c_view[3, 1] == 100

    stat = [SVector(1, 2), SVector(3, 4), SVector(5, 6)]
    @test @inferred(combinedims(stat)) == [1 3 5; 2 4 6]
    # infer to a union:
    @test combinedims(stat, 1) == [1 2; 3 4; 5 6]
    @test combinedims(stat, 2) == [1 3 5; 2 4 6]
    @testset for ds in [(), (1,), (2,)]
        cv = @inferred(combinedimsview(stat, ds...))
        @inferred(size(cv))
        @inferred(cv[1, 2])
    end
end


text:
#=
    generate_state.jl

This file contains several methods to generate the state associated to a
Neural Network and an hilbert space.
While the basic logic is simple, handling batches of states and diagonal
states of the density matrix complicates it.
=#

state(op::AbsLinearOperator, net, args...) =
    state(basis(op), net, args...)
state(hilb::AbstractHilbert, net::NeuralNetwork, args...) =
    state(input_type(net), hilb, net, args...)

state(T::Type{<:Number}, op::AbsLinearOperator, net, args...) =
    state(T, basis(op), net, args...)
state(T::Type{<:Number}, hilb::AbstractHilbert, net::CachedNet, args...) =
    state(T, hilb, net.net, net.cache, args...)

# If it's a cached net enforce the batch_size
state(T::Type{<:Number}, hilb::AbstractHilbert, net::NeuralNetwork, cache::NNCache, dims::Vararg{Int,n}) where n =
    state(T, hilb, net, dims...)
state(T::Type{<:Number}, hilb::AbstractHilbert, net::NeuralNetwork, cache::NNBatchedCache, dims::Vararg{Int,n}) where n =
    state(T, hilb, net, batch_size(cache), dims...)

# Those are the inner-most dispatches
state(T::Type{<:Number}, hilb::AbstractHilbert,
      net::Union{MatrixNeuralNetwork, KetNeuralNetwork},
      dims::Vararg{Int,n}) where n =
    state(trainable_first(net), T, hilb, dims...)


text:
module MyExample

# Write your package code here.
great() = print("Hello world")
include("extra_file.jl")
export my_f
end


text:
"""
    liftingline_section_model()

Construct a lifting line section model with state variables ``v_x, v_y, v_z, \\omega_x,
\\omega_y, \\omega_z``, inputs ``f_x, f_y, f_z, m_x, m_y, m_z``, and parameters
``\\rho, c``.  Two-dimensional aerodynamic models may be extended to three-dimensional
models by coupling with this model.  Note that this model has no rate equations of its own
since its state variables are defined as functions of the 3D structural model's state
variables.  Also note that additional parameters may not be used when coupling with this
model, this is a limitation of the [`LiftingLine`](@ref) model.
"""
function liftingline_section_model()

    # residual function (no residual function)
    f = nothing

    # number of state, input, and parameters (use Val(N) to use inferrable dimensions)
    nx = Val(6)
    ny = Val(6)
    np = Val(2)

    # convenience functions for setting states, inputs, and parameters
    setstate = liftingline_section_setstate!
    setinput = liftingline_section_setinput!
    setparam = liftingline_section_setparam!

    # convenience functions for separating states, inputs, and parameters
    sepstate = liftingline_section_sepstate
    sepinput = liftingline_section_sepinput
    sepparam = liftingline_section_sepparam

    return Model{false}(f, nx, ny, np;
        setstate = setstate,
        setinput = setinput,
        setparam = setparam,
        sepstate = sepstate,
        sepinput = sepinput,
        sepparam = sepparam)
end

# --- Internal Methods for this Model --- #

# convenience function for defining this model's state vector
function liftingline_section_setstate!(x; v, omega)
    x[1] = v[1]
    x[2] = v[2]
    x[3] = v[3]
    x[4] = omega[1]
    x[5] = omega[2]
    x[6] = omega[3]
    return x
end

# convenience function for defining this model's input vector
function liftingline_section_setinput!(y; f, m)
    y[1] = f[1]
    y[2] = f[2]
    y[3] = f[3]
    y[4] = m[1]
    y[5] = m[2]
    y[6] = m[3]
    return y
end

# convenience function for defining this model's parameter vector
function liftingline_section_setparam!(p; rho, c)
    p[1] = rho
    p[2] = c
    return p
end

# convenience function for separating this model's state vector
liftingline_section_sepstate(x) = (v=view(x, 1:3), omega=view(x, 4:6))

# convenience function for separating this model's input vector
liftingline_section_sepinput(y) = (f=view(y, 1:3), m=view(y, 4:6))

# convenience function for separating this model's parameter vector
liftingline_section_sepparam(p) = (rho=p[1], c=p[2])

# --- Internal Methods for Couplings with this Model --- #

# local freestream linear/angular velocities
function liftingline_section_velocities(vx, vz, ωy)
    u = vx
    v = vz
    ω = ωy
    return u, v, ω
end

# local freestream linear/angular accelerations
function liftingline_section_accelerations(dvx, dvz, dωy)
    udot = dvx
    vdot = dvz
    ωdot = dωy
    return udot, vdot, ωdot
end


text:
#: Abstract functional types

export MinimizableFunction, DifferentiableFunction, ProximableFunction, ProjectionableSet
export OptimOptions, minimize, minimize!, proxy, proxy!, project, project!, grad, grad!

# Abstract type declarations

"""Expected behavior for MinimizableFunction:
- minimize!(f::MinimizableFunction{T,N}, x0::AbstractArray{T,N}) where {T,N}
It approximates the solution of the optimization problem: min_x f(x)
"""
abstract type MinimizableFunction{T,N} end

abstract type OptimOptions end

minimize(fun::MinimizableFunction{T,N}, x0::AbstractArray{T,N}, opt::OptimOptions) where {T,N} = minimize!(fun, x0, opt, similar(x0))
minimize_debug(fun::MinimizableFunction{T,N}, x0::AbstractArray{T,N}, opt::OptimOptions) where {T,N} = minimize_debug!(fun, x0, opt, similar(x0))

"""Expected behavior for DifferentiableFunction:
- fval::T = grad!(f::DifferentiableFunction{DT}, x::DT, g::DT) where {T,N,DT<:AbstractArray{T,N}}
"""
abstract type DifferentiableFunction{T,N} end

function grad(f::DifferentiableFunction{T,N}, x::AbstractArray{T,N}) where {T,N}
    g = similar(x)
    fval = grad!(f, x, g)
    return fval, g
end

(f::DifferentiableFunction{T,N})(x::AbstractArray{T,N}) where {T,N} = grad!(f, x, nothing)

"""Expected behavior for ProximableFunction:
- proxy!(y::AbstractArray{T,N}, λ::T, g::ProximableFunction{T,N}, x::AbstractArray{T,N}) where {T,N}
- proxy!(y::AbstractArray{T,N}, λ::T, g::ProximableFunction{T,N}, x::AbstractArray{T,N}, opt::OptimOptions) where {T,N}
It approximates the solution of the optimization problem: min_x 0.5*norm(x-y)^2+λ*g(x)
- project!(y::AbstractArray{T,N}, ε::T, g::ProximableFunction{T,N}, x::AbstractArray{T,N}) where {T,N}
- project!(y::AbstractArray{T,N}, ε::T, g::ProximableFunction{T,N}, x::AbstractArray{T,N}, opt::OptimOptions) where {T,N}
It approximates the solution of the optimization problem: min_{g(x)<=ε} 0.5*norm(x-y)^2
"""
abstract type ProximableFunction{T,N} end

proxy(y::AbstractArray{CT,N}, λ::T, g::ProximableFunction{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = proxy!(y, λ, g, similar(y))
project(y::AbstractArray{CT,N}, ε::T, g::ProximableFunction{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = project!(y, ε, g, similar(y))
proxy(y::AbstractArray{CT,N}, λ::T, g::ProximableFunction{CT,N}, opt::OptimOptions) where {T<:Real,N,CT<:RealOrComplex{T}} = proxy!(y, λ, g, similar(y), opt::OptimOptions)
project(y::AbstractArray{CT,N}, ε::T, g::ProximableFunction{CT,N}, opt::OptimOptions) where {T<:Real,N,CT<:RealOrComplex{T}} = project!(y, ε, g, similar(y), opt::OptimOptions)

"""Projectional sets
Expected behavior for convex sets: y = project!(x, C, y), y = Π_C(x)
"""
abstract type ProjectionableSet{T,N} end

project(x::AbstractArray{T,N}, C::ProjectionableSet{T,N}) where {T,N} = project!(x, C, similar(x))
project(x::AbstractArray{T,N}, C::ProjectionableSet{T,N}, opt::OptimOptions) where {T,N} = project!(x, C, similar(x), opt)

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
@enum L_Incomes begin
    WAGES = 1
    SELF_EMPLOYMENT_INCOME = 2
    ODD_JOBS = 3
    PRIVATE_PENSIONS = 4
    NATIONAL_SAVINGS = 5
    BANK_INTEREST = 6
    STOCKS_SHARES = 7
    INDIVIDUAL_SAVINGS_ACCOUNT = 8
    PROPERTY = 9
    ROYALTIES = 10
    BONDS_AND_GILTS = 11
    OTHER_INVESTMENT_INCOME = 12
    OTHER_INCOME = 13
    ALIMONY_AND_CHILD_SUPPORT_RECEIVED = 14
    PRIVATE_SICKNESS_SCHEME_BENEFITS = 15
    ACCIDENT_INSURANCE_SCHEME_BENEFITS = 16
    HOSPITAL_SAVINGS_SCHEME_BENEFITS = 17
    UNEMPLOYMENT_OR_REDUNDANCY_INSURANCE = 18
    PERMANENT_HEALTH_INSURANCE = 19
    ANY_OTHER_SICKNESS_INSURANCE = 20
    CRITICAL_ILLNESS_COVER = 21
    TRADE_UNION_SICK_OR_STRIKE_PAY = 22
    SPARE_INC_1 = 23
    SPARE_INC_2 = 24
    SPARE_INC_3 = 25
    SPARE_INC_4 = 26
    SPARE_INC_5 = 27
    HEALTH_INSURANCE = 28
    ALIMONY_AND_CHILD_SUPPORT_PAID = 29
    TRADE_UNIONS_ETC = 30
    FRIENDLY_SOCIETIES = 31
    WORK_EXPENSES = 32
    AVCS = 33
    OTHER_DEDUCTIONS = 34
    LOAN_REPAYMENTS = 35
    PENSION_CONTRIBUTIONS_EMPLOYEE = 36
    PENSION_CONTRIBUTIONS_EMPLOYER = 37
    SPARE_DEDUCT_1 = 38
    SPARE_DEDUCT_2 = 39
    SPARE_DEDUCT_3 = 40
    SPARE_DEDUCT_4 = 41
    SPARE_DEDUCT_5 = 42

    INCOME_TAX = 43
    NATIONAL_INSURANCE = 44
    LOCAL_TAXES = 45
    SOCIAL_FUND_LOAN_REPAYMENT = 46
    STUDENT_LOAN_REPAYMENTS = 47
    CARE_INSURANCE = 48
    SPARE_TAX_1 = 49
    SPARE_TAX_2 = 50
    SPARE_TAX_3 = 51
    SPARE_TAX_4 = 52
    SPARE_TAX_5 = 53

    CHILD_BENEFIT = 54
    STATE_PENSION = 55
    BEREAVEMENT_ALLOWANCE = 56
    ARMED_FORCES_COMPENSATION_SCHEME = 57
    WAR_WIDOWS_PENSION = 58
    SEVERE_DISABILITY_ALLOWANCE = 59
    ATTENDANCE_ALLOWANCE = 60
    CARERS_ALLOWANCE = 61
    INDUSTRIAL_INJURY_BENEFIT = 62
    INCAPACITY_BENEFIT = 63
    PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING = 64
    PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY = 65
    DLA_SELF_CARE = 66
    DLA_MOBILITY = 67
    EDUCATION_ALLOWANCES = 68
    FOSTER_CARE_PAYMENTS = 69
    MATERNITY_ALLOWANCE = 70
    MATERNITY_GRANT = 71
    FUNERAL_GRANT = 72
    ANY_OTHER_NI_OR_STATE_BENEFIT = 73
    FRIENDLY_SOCIETY_BENEFITS = 74
    GOVERNMENT_TRAINING_ALLOWANCES = 75
    CONTRIB_JOBSEEKERS_ALLOWANCE = 76
    GUARDIANS_ALLOWANCE = 77
    WIDOWS_PAYMENT = 78
    WINTER_FUEL_PAYMENTS = 79
   # legacy mt benefits
    WORKING_TAX_CREDIT = 80
    CHILD_TAX_CREDIT = 81
    NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE = 82
    INCOME_SUPPORT = 83
    PENSION_CREDIT = 84
    SAVINGS_CREDIT = 85
    NON_CONTRIB_JOBSEEKERS_ALLOWANCE = 86 # FIXME JOBSEEKERS->JOB_SEEKERS everwhere
    HOUSING_BENEFIT = 87

    FREE_SCHOOL_MEALS = 88
    UNIVERSAL_CREDIT = 89
    OTHER_BENEFITS = 90
    STUDENT_GRANTS = 91
    STUDENT_LOANS = 92
    COUNCIL_TAX_BENEFIT = 93
    CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE = 94
    SCOTTISH_CARERS_SUPPLEMENT = 95
    SPARE_BEN_1 = 96
    SPARE_BEN_2 = 97
    SPARE_BEN_3 = 98
    SPARE_BEN_4 = 99
    SPARE_BEN_5 = 100
    SPARE_BEN_6 = 101
    SPARE_BEN_7 = 102
    SPARE_BEN_8 = 103
    SPARE_BEN_9 = 104
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function test_bracket()
    f(x) = cos(x)

    a = 0
    b = 0.1

    bracketing_int = opti.bracket_min(f, a, b, )
    show(bracketing_int)

end

function test_golden()
    # fcn(x) = cos(x)
    fcn(x) = x^2

    a = 0
    b = 0.1
    bracketing_int = opti.bracket_min(fcn, a, b)

    stat = opti.golden_section_search(fcn, bracketing_int.x_lo, bracketing_int.x_hi)

    println("The found min: $(stat.x_sol) with value: $(stat.y_sol) with iterations: $(stat.iter)")

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Prevent Zygote from differentiating push!
# See https://github.com/TuringLang/Turing.jl/issues/1199
ZygoteRules.@adjoint function push!(
    vi::VarInfo,
    vn::VarName,
    r,
    dist::Distribution,
    gidset::Set{Selector}
)
    return push!(vi, vn, r, dist, gidset), _ -> nothing
end


text:
abstract type HilbertAlgorithm{T} end

index_type(::HilbertAlgorithm{T}) where {T} = T

"""
    encode_hilbert(ha::HilbertAlgorithm{T}, X::Vector{A})

A 1-based Hilbert encoding. Both the Hilbert index and the axes start counting
at 1 instead of 0.
"""
function encode_hilbert(gg::HilbertAlgorithm{T}, X::Vector{A}) where {A, T}
    encode_hilbert_zero(gg, X .- one(A)) + one(T)
end

"""
    decode_hilbert!(ha::HilbertAlgorithm{T}, X::Vector{A})

A 1-based Hilbert decode, from [`decode_hilbert_zero!`](@ref). Both the Hilbert
index and the axes start counting at 1 instead of 0.
"""
function decode_hilbert!(gg::HilbertAlgorithm{T}, X::Vector{A}, h::T) where {A,T}
    decode_hilbert_zero!(gg, X, h - one(T))
    X .+= one(A)
end


text:
export Elf64_Ehdr,
    Elf64_Shdr,
    Elf64_Phdr,
    Elf64_Rel,
    Elf64_Rela,
    Elf64_Sym,
    ELF64_ST_TYPE,
    ELF64_R_SYM,
    Elf64_Chdr,
    Elf64_Dyn

const Elf64_Word = UInt32
const Elf64_Sword = Int32
const Elf64_Xword = UInt64
const Elf64_Sxword = Int64
const Elf64_Half = UInt16
const Elf64_Addr = UInt64
const Elf64_Off = UInt64
const Elf64_Section = UInt16
const Elf64_Versym = Elf64_Half

"""
    Elf64_Ehdr

A type represents 64bit ELF header.
"""
struct Elf64_Ehdr
    e_ident::SVector{16,UInt8}
    e_type::Elf64_Half
    e_machine::Elf64_Half
    e_version::Elf64_Word
    e_entry::Elf64_Addr
    e_phoff::Elf64_Off
    e_shoff::Elf64_Off
    e_flags::Elf64_Word
    e_ehsize::Elf64_Half
    e_phentsize::Elf64_Half
    e_phnum::Elf64_Half
    e_shentsize::Elf64_Half
    e_shnum::Elf64_Half
    e_shstrndx::Elf64_Half
end

Elf64_Ehdr(bin::Vector{UInt8}) = pointer(bin) |> Ptr{Elf64_Ehdr} |> unsafe_load
Elf64_Ehdr(bin::Vector{UInt8}, off::UInt) = pointer(bin, off) |> Ptr{Elf64_Ehdr} |> unsafe_load

"""
    Elf64_Shdr

A type represents 64bit ELF section header.
"""
struct Elf64_Shdr
    sh_name::Elf64_Word
    sh_type::Elf64_Word
    sh_flags::Elf64_Xword
    sh_addr::Elf64_Addr
    sh_offset::Elf64_Off
    sh_size::Elf64_Xword
    sh_link::Elf64_Word
    sh_info::Elf64_Word
    sh_addralign::Elf64_Xword
    sh_entsize::Elf64_Xword
end

Elf64_Shdr(bin::Vector{UInt8}, off::UInt) =
    pointer(bin, off + 1) |> Ptr{Elf64_Shdr} |> unsafe_load

"""
    Elf64_Phdr

A type represents 64bit ELF program header.
"""
struct Elf64_Phdr
    p_type::Elf64_Word
    p_flags::Elf64_Word
    p_offset::Elf64_Off
    p_vaddr::Elf64_Addr
    p_paddr::Elf64_Addr
    p_filesz::Elf64_Xword
    p_memsz::Elf64_Xword
    p_align::Elf64_Xword
end

Elf64_Phdr(bin::Vector{UInt8}, off::UInt) =
    pointer(bin, off + 1) |> Ptr{Elf64_Phdr} |> unsafe_load

"""
    Elf64_Sym

A type represents 64bit ELF Symbol Table entry.
"""
struct Elf64_Sym
    st_name::Elf64_Word
    st_info::UInt8
    st_other::UInt8
    st_shndx::Elf64_Section
    st_value::Elf64_Addr
    st_size::Elf64_Xword
end

Elf64_Sym(bin::Vector{UInt8}, off::UInt) =
    pointer(bin, off + 1) |> Ptr{Elf64_Sym} |> unsafe_load

st_info_bind(st_info::UInt8) = st_info >>> 4
st_info_type(st_info::UInt8) = st_info & 0xf
st_info_info(bind::UInt8, type::UInt8) = ((bind << 4) + (type & 0xf))

st_visibility(st_other::UInt8) = st_other & 0x3

struct Elf64_SymInfo
    si_boundto::Elf64_Half
    si_flags::Elf64_Half
end

Elf64_SymInfo(bin::Vector{UInt8}, off::UInt8) =
    pointer(bin, off + 1) |> Ptr{Elf64_SymInfo} |> unsafe_load

"""
    Elf64_Rel

A type represents 64bit ELF relocasion table entry (SHT_REL).
"""
struct Elf64_Rel
    r_offset::Elf64_Addr
    r_info::Elf64_Xword
end

Elf64_Rel(bin::Vector{UInt8}, off::UInt) =
    pointer(bin, off + 1) |> Ptr{Elf64_Rel} |> unsafe_load

"""
    Elf64_Rela

A type represents 64bit ELF relocasion table entry (SHT_RELA).
"""
struct Elf64_Rela
    r_offset::Elf64_Addr
    r_info::Elf64_Xword
    r_addend::Elf64_Sxword
end

Elf64_Rela(bin::Vector{UInt8}, off::UInt) =
    pointer(bin, off + 1) |> Ptr{Elf64_Rela} |> unsafe_load

# struct SymbolTable64
    # syms::Vector{Union{Elf64_Rel, Elf64_Rela}}

    # function SymbolTable64(relsec::Elf64_Shdr)
        # for i=
# end

struct Elf64_Chdr
    ch_type::Elf64_Word # Compression format.
    ch_size::Elf64_Word # Uncompressed size.
    ch_addralign::Elf64_Xword # Uncompressed data alignment.
end

Elf64_Chdr(bin::Vector{UInt8}, off::UInt) = pointer(bin, off + 1) |> Ptr{Elf64_Chdr} |> unsafe_load

struct Elf64_Dyn
    d_tag::Elf64_Xword
    d_val_or_ptr::UInt64
end

Elf64_Dyn(bin::Vector{UInt8}, off::UInt) = pointer(bin, off + 1) |> Ptr{Elf64_Dyn} |> unsafe_load


text:
using Test, CommonTests

include("example.jl")
using Main.Example

@testset "MyTypeA test" begin
    global TERGET = MyTypeA
    include("common.jl")

    @testset "greet test" begin
        @runcommon

        @test startswith(greet(MyTypeA), "Hello, ")
    end

end

@testset "MyTypeB test" begin
    global TERGET = MyTypeB
    include("common.jl")

    @testset "greet test" begin
        @runcommon

        @test startswith(greet(MyTypeB), "Hi, ")
    end

end


text:
# ============================================================
# Thermodynamics
# ============================================================

"""
    heat_capacity_ratio(K, D)
    heat_capacity_ratio(K, Nr, D)

Calculate heat capacity ratio

"""
function heat_capacity_ratio(K, D::T) where {T<:Integer}
    γ = begin
        if D == 1
            (K + 3.0) / (K + 1.0)
        elseif D == 2
            (K + 4.0) / (K + 2.0)
        elseif D == 3
            (K + 5.0) / (K + 3.0)
        end
    end

    return γ
end

#--- diatomic ---#
function heat_capacity_ratio(K, Nr, D::T) where {T<:Integer}
    γ = begin
        if D == 1
            (K + 3.0 + Nr) / (K + 1.0 + Nr)
        elseif D == 2
            (K + 4.0 + Nr) / (K + 2.0 + Nr)
        elseif D == 3
            (K + 5.0 + Nr) / (K + 3.0 + Nr)
        end
    end

    return γ
end

"""
    sound_speed(λ, γ)
    sound_speed(prim, γ)

Calculate speed of sound

"""
sound_speed(λ::Real, γ::Real) = (0.5 * γ / λ)^0.5

sound_speed(prim::AbstractVector{T}, γ) where {T<:Real} = sound_speed(prim[end], γ)

#--- mixture ---#
function sound_speed(prim::AbstractMatrix{T}, γ) where {T<:Real}
    c = similar(prim, axes(prim, 2))
    for j in eachindex(c)
        c[j] = sound_speed(prim[end, j], γ)
    end

    return maximum(c)
end


text:
"""
	struct ExtractArray{T}
		item::T
	end

Convert array of values to a `Mill.BagNode` with items converted
by `item`. The entire array is assumed to be a single bag.

# Examples

```jldoctest
julia> ec = ExtractArray(ExtractCategorical(2:4));

julia> ec([2,3,1,4]).data
4×4 Mill.ArrayNode{Array{Bool,2},Nothing}:
 1  0  0  0
 0  1  0  0
 0  0  0  1
 0  0  1  0

julia> es = ExtractArray(ExtractScalar());

julia> es([2,3,4])
BagNode with 1 obs
  └── ArrayNode(1×3 Array with Float32 elements) with 3 obs

julia> es([2,3,4]).data
1×3 Mill.ArrayNode{Array{Float32,2},Nothing}:
 2.0  3.0  4.0
```
"""
struct ExtractArray{T} <: AbstractExtractor
	item::T
end

function (s::ExtractArray)(v::V) where {V<:Union{Missing, Nothing}}
	Mill._emptyismissing[] && return(BagNode(missing, [0:-1]))
	BagNode(s.item(extractempty), [0:-1])
end

(s::ExtractArray)(v::V) where {V<:Vector} = isempty(v) ? s(missing) : BagNode(reduce(catobs, map(s.item, v)),[1:length(v)])
(s::ExtractArray)(v) = s(missing)
(s::ExtractArray)(v::ExtractEmpty) = BagNode(s.item(extractempty), Mill.AlignedBags(Array{UnitRange{Int64},1}()))

Base.hash(e::ExtractArray, h::UInt) = hash(e.item, h)
Base.:(==)(e1::ExtractArray, e2::ExtractArray) = e1.item == e2.item


text:
function printcircle(lo::Integer, hi::Integer, ndots::Integer; pad::Integer = 2)
    canvas = falses(2hi + 1, 2hi + 1)
    i = 0
    while i < ndots
        x, y = rand(-hi:hi, 2)
        if lo ^ 2 - 1 < x ^ 2 + y ^ 2 < hi ^ 2 + 1
            canvas[x + hi + 1, y + hi + 1] = true
            i += 1
        end
    end
    # print
    for i in 1:(2hi + 1)
        row = map(j -> j ? "\u25cf " : "  ", canvas[i, :])
        println(" " ^ pad, join(row))
    end
    return canvas
end

printcircle(10, 15, 100)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Pkg
ipks = Pkg.installed()
pkgs = ["BenchmarkTools","UnicodePlots","PyPlot"]
for p in pkgs
    !haskey(ipks,p) && Pkg.add(p)
end
#
# We can define these here as as required
using BenchmarkTools, UnicodePlots, PyPlot
using Printf
#
# Define a function to compute the sum of the squares of (x,y)
sumsq(x,y) = x*x + y*y;
#
# Set the number of trials and initialise the counter
# Using zero(Integer) will ensure the counter is an integer

N = 10^8;
K = zero(Integer);
#
# Sum the number of pairs which lie within the inscribed circle
for i = 1:N
  if sumsq(rand(), rand()) < 1.0
    global K += 1
  end
end
#
# ... and output the value for PI
@printf "Estimate of PI for %d trials is %8.5f\n" N 4.0*(K / N);

# We can time this using the @time macro
# Note that his takes some time, will be different in the REPL
@time for i = 1:N
  if sumsq(rand(), rand()) < 1.0
    global K += 1
  end
end

#=
We can rewrite as a generator expressions
This does not require the variable 'K'
Comprehensions can also be written without the enclosing square brackets,
producing an object known as a generator.
=#
((rand()^2 + rand()^2) < 1.0 for i = 1:N)

# The generator can be "evaluated" as the argument to an aggregate function
count((rand()^2 + rand()^2) < 1.0 for i = 1:N) * 4 / N

# And this runs much more quickly with a much better memory allocation
@time count((rand()^2 + rand()^2) < 1.0 for i = 1:N) * 4 / N

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

# Rather than just using the @time macro we can get a more accuracy timing
# ... using the BenchmarkTools package

using BenchmarkTools

# Define the function to sum the series
# It converges very slowly

function basel(N::Integer)
  @assert N > 0
  s = 0.0
  for i = 1:N
    s += 1.0/float(i)^2
  end
  return s
end

basel(10^8)   # Evaluate it over 10^8 terms
π*π/6.0    # Compare this to the analytic solution

# Use the BenchmarkTools macro @benchmark to get more unbiased timings
@benchmark basel(10^8) samples=10

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

using UnicodePlots

# Generate an array of the numbers from 1 to 100
# The ordinate value is create using a list comprehensive

x = collect(1:100);
y = [x[i]*sin(0.3*x[i])*exp(-0.03*x[i]) for i = 1:length(x)];

myPlot = lineplot(x, y, title = "My Plot", name = "chp-01")

# Alternatively this can be done using a map() construction

t = collect(0.0:0.1:10.0);
y =  map(x -> x*sin(3.0*x)*exp(-0.3*x), t);
myPlot = lineplot(t, y, title = "Second Plot", name = "chp-01")

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

using PyPlot
plt = PyPlot

S0  = 100;      # Spot price
K   = 102;      # Strike price
r   = 0.05;     # Risk free rate
q   = 0.0;      # Dividend yield
v   = 0.2;      # Volatility
tma = 0.25;     # Time to maturity
T   = 90;       # Number of time steps
dt  = tma/T;

N = T + 1;
x = collect(0:T);

plt.title("Asian Option trajectories");
plt.xlabel("Time");
plt.ylabel("Stock price");

for k = 1:5
  S = zeros(Float64,N)
  S[1] = S0;
  dW = randn(N)*sqrt(dt);
  [ S[t] = S[t-1] * (1 + (r - q - 0.5*v*v)*dt + v*dW[t] + 0.5*v*v*dW[t]*dW[t]) for t=2:N ]
  plt.plot(x,S)
end


text:
using Random
using LinearAlgebra
using GalacticOptim,Optim,NLopt#,BlackBoxOptim
using DifferentialEquations#,DiffEqFlux
#using ForwardDiff
using Zygote
using Plots
#using ReverseDiff
include("../src/GLOQ.jl")
pyplot()
# BlackBoxOptim somehow downgrade some packages and as a result breaks the auto-differentiation with Zygote
# we should avoid it.
# NLopt seems to be really fast. Their LBFGS with box constrains is faster than Optim's FMinbox(LBFGS)

N_states = 4;
freqs = [4.10817777; 3.88303940; 3.61937188]
omegas = 2*pi.*freqs
gamma1   = [2.59451982e-05; 4.54296537e-05; 0.0]
gamma2   = [2.58005604e-05; 9.00000000e-05; 0.0]
omr_ramsey = 2.0*pi*(3.8830355 - 2.0e-3)
omr_echo = 2.0*pi*3.8830355
omr_T1 = 2.0*pi*3.8830355

width = 17.00
TC = 2.5*width
Ω = 0.5*pi/(TC*sqrt(2.0))
θ = 0.0

fromState = 1
rho_u0 = [0.0;0.0;0.0;0.0]
rho_v0 = [0.0;0.0;0.0;0.0]
rho_u0[fromState+1] = 1.0

T_Ramsey = 1.0*GLOQ.GLOQ_MICRO_SEC#10.0 * 1000
N_Ramsey = 101
t_ramsey = zeros(Float64,N_Ramsey)
dt_ramsey = T_Ramsey/(N_Ramsey-1)
for i = 1:N_Ramsey
	dark_time = (i-1)*dt_ramsey
	t_ramsey[i] = dark_time
end

# Ramsey Forward Solve
@time rho_ramsey_u_exp,rho_ramsey_v_exp = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
			     omegas,omr_ramsey,
				 gamma1,gamma2,
				 fromState,
				 TC,t_ramsey,N_states;
				 initial_type="states")
@time rho_ramsey_u,rho_ramsey_v = GLOQ.RamseyForwardSolve(rho_u0,rho_v0,
			     omegas,omr_ramsey,
				 gamma1,gamma2,
				 fromState,
				 TC,t_ramsey,N_states;
				 method = Trapezoid(),
				 abstol = 1e-10,
				 restol = 1e-10
				 )
population_Ramsey = GLOQ.get_population(rho_ramsey_u)
population_Ramsey_exp = GLOQ.get_population(rho_ramsey_u_exp)

println("Ramsey max error:",maximum(abs.(population_Ramsey-population_Ramsey_exp)) )

fig=plot(t_ramsey./GLOQ.GLOQ_MICRO_SEC,population_Ramsey)
plot!(fig,t_ramsey./GLOQ.GLOQ_MICRO_SEC,population_Ramsey_exp)
display(fig)

fig2=plot(t_ramsey./GLOQ.GLOQ_MICRO_SEC,population_Ramsey-population_Ramsey_exp)
display(fig2)

# Echo Forward Solve
T_Echo = 10.0*GLOQ.GLOQ_MICRO_SEC#10.0 * 1000
N_Echo = 1001
t_echo = zeros(Float64,N_Echo)
dt_echo = T_Echo/(N_Echo-1)
for i = 1:N_Echo
	dark_time = (i-1)*dt_echo
	t_echo[i] = dark_time
end
@time rho_echo_u_exp,rho_echo_v_exp = GLOQ.EchoForwardSolve(rho_u0,rho_v0,
			     omegas,omr_echo,
				 gamma1,gamma2,
				 fromState,
				 TC,t_echo,N_states;
				 initial_type = "states")
@time rho_echo_u,rho_echo_v = GLOQ.EchoForwardSolve(rho_u0,rho_v0,
			     omegas,omr_echo,
				 gamma1,gamma2,
				 fromState,
				 TC,t_echo,N_states;
				 initial_type = "states",
				 method = Trapezoid(),
				 abstol = 1e-10,
				 restol = 1e-10)
population_Echo = GLOQ.get_population(rho_echo_u)
population_Echo_exp = GLOQ.get_population(rho_echo_u_exp)

println( maximum(abs.(population_Echo-population_Echo_exp)) )

fig_echo=plot(t_echo./GLOQ.GLOQ_MICRO_SEC,population_Echo)
plot!(fig_echo,t_echo./GLOQ.GLOQ_MICRO_SEC,population_Echo_exp,line=(:dash))
display(fig_echo)

fig_echo2=plot(t_echo./GLOQ.GLOQ_MICRO_SEC,population_Echo-population_Echo_exp)
display(fig_echo2)

# T1 Forward Solve
T_T1 = 10.0*1000
N_T1 = 1001
t_T1 = zeros(Float64,N_T1)
dt_T1 = T_T1/(N_T1-1)
for i = 1:N_T1
	dark_time = (i-1)*dt_T1
	t_T1[i] = dark_time
end
@time rho_t1_u_exp,rho_t1_v_exp = GLOQ.T1ForwardSolve(rho_u0,rho_v0,
			     omegas,omr_T1,
				 gamma1,gamma2,
				 fromState,
				 TC,t_T1,N_states;
				 initial_type = "states")
population_T1_exp = GLOQ.get_population(rho_t1_u_exp)
@time rho_t1_u,rho_t1_v = GLOQ.T1ForwardSolve(rho_u0,rho_v0,
			     omegas,omr_T1,
				 gamma1,gamma2,
				 fromState,
				 TC,t_T1,N_states;
				 initial_type = "states",
				 method = Trapezoid(),
				 abstol = 1e-10,
				 restol = 1e-10)
population_T1 = GLOQ.get_population(rho_t1_u)

println( maximum(abs.(population_T1-population_T1_exp)) )

fig_t1=plot(t_T1./GLOQ.GLOQ_MICRO_SEC,population_T1)
plot!(fig_t1,t_T1./GLOQ.GLOQ_MICRO_SEC,population_T1_exp,line=(:dash))
display(fig_t1)

fig_t1_error=plot(t_T1./GLOQ.GLOQ_MICRO_SEC,population_T1-population_T1_exp)
display(fig_t1_error)


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 5b2ee40e-a2b8-11ea-0fef-c35fe6918860
md"""
# The tower of Hanoi Harry want to understand this

The tower of hanoi is a famous puzzle.

![setup of the tower of a hanoi](https://upload.wikimedia.org/wikipedia/commons/0/07/Tower_of_Hanoi.jpeg)

The game consists of three rods with disks stacked on top of them. The puzzle will start with all disks in a stack on one of the rods (like in the picture). The goal is to move all the discs to a single stack on  different rod.

To move the disks, you have to follow the following rules:

* You can move only one disk at a time.
* For each move, you have to take the upper disk from one of the stacks, and place it on top of another stack or empty rod.
* You cannot place a larger disk on top of a smaller disk.

This notebook will define a Julia implementation of the puzzle. It's up to you to write an algorithm that solves it.
"""

# ╔═╡ 95fbd0d2-a2b9-11ea-0682-fdf783251797
md"""
## Setting up the game pieces

What does a Julia implementation look like? We're not really interested in writing code that will manipulate physical disks. Our final goal is a function that will give us a _recipe_ to solve the tower of hanoi, which is just a list of moves to make. Because of that, we can use a lot of abstraction in our implementation, and keep the data structures as simple as possible.

To start, we have to define some representation of the disks and the stacks. The disks have one important property, which is that they are ordered. We can use integers to represent them.
"""

# ╔═╡ 620d6834-a2ba-11ea-150a-2132bb54e4b3
num_disks = 8

# ╔═╡ 35ada214-a32c-11ea-0da3-d5d494b28467
md"""(Side note: the number of disks is arbitrary. When testing your function, you may want to set it to 1 or 2 to start.)"""

# ╔═╡ 7243cc8e-a2ba-11ea-3f29-99356f0cdcf4
all_disks = 1:num_disks

# ╔═╡ 7e1ba2ac-a2ba-11ea-0134-2f61ed75be18
md"""
A single stack can be represented as an array with all the disks in it. We will list them from top to bottom.
"""

# ╔═╡ 43781a52-a339-11ea-3803-e56e2d08aa83
first_stack = collect(all_disks)

# ╔═╡ b648ab70-a2ba-11ea-2dcc-55b630e44325
md"""
Now we have to make three of those.
"""

# ╔═╡ 32f26f80-a2bb-11ea-0f2a-3fc631ada63d
starting_stacks = [first_stack, [], []]

# ╔═╡ e347f1de-a2bb-11ea-06e7-87cca6f2a240
md"""
## Defining the rules

Now that we have our "game board", we can implement the rules.

To start, we make two functions for states. A state of the game is just an array of stacks.

We will define a function that checks if a state is okay according to the rules. To be legal, all the stacks should be in the correct order, so no larger disks on top of smaller disks.

Another good thing to check: no disks should have appeared or disappeared since we started!
"""

# ╔═╡ 512fa6d2-a2bd-11ea-3dbe-b935b536967b
function islegal(stacks)
	order_correct = all(issorted, stacks)

	# check if we use the same disk set that we started with

	disks_in_state = sort([disk for stack in stacks for disk in stack])
	disks_complete = disks_in_state == all_disks

	order_correct && disks_complete
end

# ╔═╡ c56a5858-a2bd-11ea-1d96-77eaf5e74925
md"""
Another function for states: check if we are done! We can assume that we already checked if the state was legal. So we know that all the disks are there and they are ordered correctly.  To check if we are finished, we just need to check if the last stack contains all the disks.
"""

# ╔═╡ d5cc41e8-a2bd-11ea-39c7-b7df8de6ae3e
function iscomplete(stacks)
	last(stacks) == all_disks
end

# ╔═╡ 53374f0e-a2c0-11ea-0c91-97474780721e
md"""
Now the only rules left to implement are the rules for moving disks.

We could implement this as another check on states, but it's easier to just write a legal `move` function. Your solution will specify moves for the `move` function, so this will be the only way that the stacks are actually manipulated. That way, we are sure that nothing fishy is happening.

We will make our `move` function so that its input consists of a state of the game, and instructions for what to do. Its output will be the new state of the game.

So what should those instructions look like? It may seem intuitive to give a _disk_ that should be moved, but that's more than we need. After all, we are only allowed to take the top disk from one stack, and move it to the top of another. So we only have to say which _stacks_ we are moving between.

(Note that the `move` function is okay with moving a larger disk on top of a smaller disk. We already implemented that restriction in `islegal`.)
"""

# ╔═╡ e915394e-a2c0-11ea-0cd9-1df6fd3c7adf
function move(stacks, source::Int, target::Int)
	# check if the from stack if not empty
	if isempty(stacks[source])
		error("Error: attempted to move disk from empty stack")
	end

	new_stacks = deepcopy(stacks)

	disk = popfirst!(new_stacks[source]) # take disk
	pushfirst!(new_stacks[target], disk) # put on new stack

	return new_stacks
end

# ╔═╡ 87b2d164-a2c4-11ea-3095-916628943879
md"""
## Solving the problem

We have implemented the game pieces and the rules, so you can start working on your solution.

To do this, you can fill in the `solve(stacks)` function. This function should give a solution for the given `stacks`, by moving all the disks from stack 1 to stack 3.

As output, `solve` should give a recipe, that tells us what to do. This recipe should be an array of moves. Each moves is a `(source, target)` tuple, specifying from which stack to which stack you should move.

For example, it might look like this:
"""

# ╔═╡ 29b410cc-a329-11ea-202a-795b31ce5ad5
function wrong_solution(stacks)::Array{Tuple{Int,Int}}
	return [(1, 2), (2, 3), (2, 1), (1, 3)]
end

# ╔═╡ ea24e778-a32e-11ea-3f11-dbe9d36b1011
md"""
Now you can work on building an actual solution. Some tips:
* `solve(stacks)` can keep track of the board if you want, but it doesn't have to.
* The section below will actually run your moves, which is very useful for checking them.
* If you want, you can change `num_disks` to 1 or 2. That can be a good starting point.
"""

# ╔═╡ 010dbdbc-a2c5-11ea-34c3-837eae17416f
function solve(stacks)::Array{Tuple{Int,Int}}

	# what to do?

	return []
end

# ╔═╡ 3eb3c4c0-a2c5-11ea-0bcc-c9b52094f660
md"""
## Checking solutions

This is where we can check a solution. We start with a function that takes our recipe and runs it.
"""

# ╔═╡ 4709db36-a327-11ea-13a3-bbfb18da84ce
function run_solution(solver::Function, start = starting_stacks)
	moves = solver(deepcopy(start)) # apply the solver

	all_states = Array{Any,1}(undef, length(moves) + 1)
	all_states[1] = starting_stacks

	for (i, m) in enumerate(moves)
		try
			all_states[i + 1] = move(all_states[i], m[1], m[2])
		catch
			all_states[i + 1] = missing
		end
	end

	return all_states
end

# ╔═╡ 372824b4-a330-11ea-2f26-7b9a1ad018f1
md"""
You can use this function to see what your solution does.

If `run_solution` tries to make an impossible move, it will give `missing` from that point onwards. Look at what happens in the `wrong_solution` version and compare it to the moves in `wrong_solution`.
"""

# ╔═╡ d2227b40-a329-11ea-105c-b585d5fcf970
run_solution(wrong_solution)

# ╔═╡ 9173b174-a327-11ea-3a69-9f7525f2e7b4
run_solution(solve)

# ╔═╡ bb5088ec-a330-11ea-2c41-6b8b92724b3b
md"""
Now that we have way to run a recipe, we can check if its output is correct. We will check if all the intermediate states are legal and the final state is the finished puzzle.
"""

# ╔═╡ 10fb1c56-a2c5-11ea-2a06-0d8c36bfa138
function check_solution(solver::Function, start = starting_stacks)
	try
		# run the solution
		all_states = run_solution(solver, start)

		# check if each state is legal
		all_legal = all(islegal, all_states)

		# check if the final state is is the completed puzzle
		complete = (iscomplete ∘ last)(all_states)

		all_legal && complete
	catch
		# return false if we encountered an error
		return false
	end
end

# ╔═╡ 8ea7f944-a329-11ea-22cc-4dbd11ec0610
check_solution(solve)

# ╔═╡ e54add0a-a330-11ea-2eeb-1d42f552ba38
if check_solution(solve)
	if num_disks >= 8
		md"""
		Congratulations, your solution works!
		"""
	else
		md"""
		Your solution works for $(num_disks) disks. Change `num_disks` to see if it works for 8 or more.
		"""
	end
else
	md"""
	The `solve` function doesn't work yet. Keep working on it!
	"""
end

# ╔═╡ Cell order:
# ╟─5b2ee40e-a2b8-11ea-0fef-c35fe6918860
# ╟─95fbd0d2-a2b9-11ea-0682-fdf783251797
# ╠═620d6834-a2ba-11ea-150a-2132bb54e4b3
# ╟─35ada214-a32c-11ea-0da3-d5d494b28467
# ╠═7243cc8e-a2ba-11ea-3f29-99356f0cdcf4
# ╟─7e1ba2ac-a2ba-11ea-0134-2f61ed75be18
# ╠═43781a52-a339-11ea-3803-e56e2d08aa83
# ╟─b648ab70-a2ba-11ea-2dcc-55b630e44325
# ╠═32f26f80-a2bb-11ea-0f2a-3fc631ada63d
# ╟─e347f1de-a2bb-11ea-06e7-87cca6f2a240
# ╠═512fa6d2-a2bd-11ea-3dbe-b935b536967b
# ╟─c56a5858-a2bd-11ea-1d96-77eaf5e74925
# ╠═d5cc41e8-a2bd-11ea-39c7-b7df8de6ae3e
# ╟─53374f0e-a2c0-11ea-0c91-97474780721e
# ╠═e915394e-a2c0-11ea-0cd9-1df6fd3c7adf
# ╟─87b2d164-a2c4-11ea-3095-916628943879
# ╠═29b410cc-a329-11ea-202a-795b31ce5ad5
# ╟─ea24e778-a32e-11ea-3f11-dbe9d36b1011
# ╠═010dbdbc-a2c5-11ea-34c3-837eae17416f
# ╟─3eb3c4c0-a2c5-11ea-0bcc-c9b52094f660
# ╠═4709db36-a327-11ea-13a3-bbfb18da84ce
# ╟─372824b4-a330-11ea-2f26-7b9a1ad018f1
# ╠═d2227b40-a329-11ea-105c-b585d5fcf970
# ╠═9173b174-a327-11ea-3a69-9f7525f2e7b4
# ╟─bb5088ec-a330-11ea-2c41-6b8b92724b3b
# ╠═10fb1c56-a2c5-11ea-2a06-0d8c36bfa138
# ╠═8ea7f944-a329-11ea-22cc-4dbd11ec0610
# ╟─e54add0a-a330-11ea-2eeb-1d42f552ba38

max_stars_repo_path: My cool notebook.jl

max_issues_repo_path: My cool notebook.jl

max_forks_repo_path: My cool notebook.jl