text:
using Lux, NNlib, Random, Test

include("../utils.jl")

rng = Random.default_rng()
Random.seed!(rng, 0)

@testset "Miscellaneous Layers" begin
    @testset "Reshape Layer" begin
        layer = ReshapeLayer((2, 3))
        println(layer)
        ps, st = Lux.setup(rng, layer)
        x = randn(rng, 6, 3)

        @test size(layer(x, ps, st)[1]) == (2, 3, 3)
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                      rtol=1.0f-3)
    end

    @testset "Flatten Layer" begin
        layer = FlattenLayer()
        println(layer)
        ps, st = Lux.setup(rng, layer)
        x = randn(rng, 6, 3, 2)

        @test size(layer(x, ps, st)[1]) == (18, 2)
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                      rtol=1.0f-3)
    end

    @testset "NoOpLayer" begin
        layer = NoOpLayer()
        println(layer)
        ps, st = Lux.setup(rng, layer)
        x = (x=2, b=5) # Something totally arbitrary

        @test layer(x, ps, st)[1] == x
        run_JET_tests(layer, x, ps, st)

        x = randn(rng, 6, 3)
        test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                      rtol=1.0f-3)
    end

    @testset "SelectDim Layer" begin
        layer = SelectDim(3, 1)
        println(layer)
        ps, st = Lux.setup(rng, layer)
        x = randn(rng, 6, 4, 3, 2)

        @test size(layer(x, ps, st)[1]) == (6, 4, 2)
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                      rtol=1.0f-3)
    end

    @testset "WrappedFunction" begin
        layer = WrappedFunction(x -> x .* x)
        println(layer)
        ps, st = Lux.setup(rng, layer)
        x = randn(rng, 6, 4, 3, 2)

        @test layer(x, ps, st)[1] == x .* x
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                      rtol=1.0f-3)
    end

    @testset "ActivationFunction" begin
        layer = ActivationFunction(tanh)
        println(layer)
        ps, st = Lux.setup(rng, layer)
        x = randn(rng, 6, 4, 3, 2)

        @test layer(x, ps, st)[1] == tanh.(x)
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                      rtol=1.0f-3)
    end
end

@testset "Containers" begin
    @testset "SkipConnection" begin
        @testset "zero sum" begin
            layer = SkipConnection(WrappedFunction(zero), (a, b) -> a .+ b)
            println(layer)
            ps, st = Lux.setup(rng, layer)
            x = randn(rng, 10, 10, 10, 10)

            @test layer(x, ps, st)[1] == x
            run_JET_tests(layer, x, ps, st)
            test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                          rtol=1.0f-3)
        end

        @testset "concat size" begin
            layer = SkipConnection(Dense(10, 10), (a, b) -> hcat(a, b))
            println(layer)
            ps, st = Lux.setup(rng, layer)
            x = randn(rng, 10, 2)

            @test size(layer(x, ps, st)[1]) == (10, 4)
            run_JET_tests(layer, x, ps, st)
            test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                          atol=1.0f-3, rtol=1.0f-3)
        end
    end

    @testset "Parallel" begin
        @testset "zero sum" begin
            layer = Parallel(+, WrappedFunction(zero), NoOpLayer())
            println(layer)
            ps, st = Lux.setup(rng, layer)
            x = randn(rng, 10, 10, 10, 10)

            @test layer(x, ps, st)[1] == x
            run_JET_tests(layer, x, ps, st)
            test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                          rtol=1.0f-3)
        end

        @testset "concat size" begin
            layer = Parallel((a, b) -> cat(a, b; dims=2), Dense(10, 10), NoOpLayer())
            println(layer)
            ps, st = Lux.setup(rng, layer)
            x = randn(rng, 10, 2)

            @test size(layer(x, ps, st)[1]) == (10, 4)
            run_JET_tests(layer, x, ps, st)
            test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                          rtol=1.0f-3)

            layer = Parallel(hcat, Dense(10, 10), NoOpLayer())
            println(layer)
            ps, st = Lux.setup(rng, layer)

            @test size(layer(x, ps, st)[1]) == (10, 4)
            run_JET_tests(layer, x, ps, st)
            test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                          rtol=1.0f-3)
        end

        @testset "vararg input" begin
            layer = Parallel(+, Dense(10, 2), Dense(5, 2), Dense(4, 2))
            println(layer)
            ps, st = Lux.setup(rng, layer)
            x = (randn(rng, 10, 1), randn(rng, 5, 1), randn(rng, 4, 1))

            @test size(layer(x, ps, st)[1]) == (2, 1)
            run_JET_tests(layer, x, ps, st)
            test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                          atol=1.0f-3, rtol=1.0f-3)
        end

        @testset "connection is called once" begin
            CNT = Ref(0)
            f_cnt = (x...) -> (CNT[] += 1; +(x...))
            layer = Parallel(f_cnt, WrappedFunction(sin), WrappedFunction(cos),
                             WrappedFunction(tan))
            ps, st = Lux.setup(rng, layer)
            Lux.apply(layer, 1, ps, st)
            @test CNT[] == 1
            run_JET_tests(layer, 1, ps, st)
            Lux.apply(layer, (1, 2, 3), ps, st)
            @test CNT[] == 2
            layer = Parallel(f_cnt, WrappedFunction(sin))
            Lux.apply(layer, 1, ps, st)
            @test CNT[] == 3
        end

        # Ref https://github.com/FluxML/Flux.jl/issues/1673
        @testset "Input domain" begin
            struct Input
                x::Any
            end

            struct L1 <: Lux.AbstractExplicitLayer end
            (::L1)(x, ps, st) = (ps.x * x, st)
            Lux.initialparameters(rng::AbstractRNG, ::L1) = (x=randn(rng, Float32, 3, 3),)
            Base.:*(a::AbstractArray, b::Input) = a * b.x

            par = Parallel(+, L1(), L1())
            ps, st = Lux.setup(rng, par)

            ip = Input(rand(Float32, 3, 3))
            ip2 = Input(rand(Float32, 3, 3))

            @test par(ip, ps, st)[1] ≈
                  par.layers[1](ip.x, ps.layer_1, st.layer_1)[1] +
                  par.layers[2](ip.x, ps.layer_2, st.layer_2)[1]
            @test par((ip, ip2), ps, st)[1] ≈
                  par.layers[1](ip.x, ps.layer_1, st.layer_1)[1] +
                  par.layers[2](ip2.x, ps.layer_2, st.layer_2)[1]
            gs = gradient((p, x...) -> sum(par(x, p, st)[1]), ps, ip, ip2)
            gs_reg = gradient(ps, ip, ip2) do p, x, y
                sum(par.layers[1](x.x, p.layer_1, st.layer_1)[1] +
                    par.layers[2](y.x, p.layer_2, st.layer_2)[1])
            end

            @test gs[1] ≈ gs_reg[1]
            @test gs[2].x ≈ gs_reg[2].x
            @test gs[3].x ≈ gs_reg[3].x
        end
    end
end

@testset "Dense" begin
    @testset "constructors" begin
        layer = Dense(10, 100)
        ps, st = Lux.setup(rng, layer)

        @test size(ps.weight) == (100, 10)
        @test size(ps.bias) == (100, 1)
        @test layer.activation == identity

        layer = Dense(10, 100, relu; bias=false)
        ps, st = Lux.setup(rng, layer)

        @test !haskey(ps, :bias)
        @test layer.activation == relu
    end

    @testset "dimensions" begin
        layer = Dense(10, 5)
        ps, st = Lux.setup(rng, layer)

        @test size(first(Lux.apply(layer, randn(10), ps, st))) == (5,)
        @test size(first(Lux.apply(layer, randn(10, 2), ps, st))) == (5, 2)
    end

    @testset "zeros" begin
        @test begin
            layer = Dense(10, 1, identity; init_weight=ones)
            first(Lux.apply(layer, ones(10, 1), Lux.setup(rng, layer)...))
        end == 10 * ones(1, 1)

        @test begin
            layer = Dense(10, 1, identity; init_weight=ones)
            first(Lux.apply(layer, ones(10, 2), Lux.setup(rng, layer)...))
        end == 10 * ones(1, 2)

        @test begin
            layer = Dense(10, 2, identity; init_weight=ones)
            first(Lux.apply(layer, ones(10, 1), Lux.setup(rng, layer)...))
        end == 10 * ones(2, 1)

        @test begin
            layer = Dense(10, 2, identity; init_weight=ones)
            first(Lux.apply(layer, [ones(10, 1) 2 * ones(10, 1)], Lux.setup(rng, layer)...))
        end == [10 20; 10 20]

        @test begin
            layer = Dense(10, 2, identity; init_weight=ones, bias=false)
            first(Lux.apply(layer, [ones(10, 1) 2 * ones(10, 1)], Lux.setup(rng, layer)...))
        end == [10 20; 10 20]
    end
end

@testset "Scale" begin
    @testset "constructors" begin
        layer = Scale(10, 100)
        ps, st = Lux.setup(rng, layer)

        @test size(ps.weight) == (10, 100)
        @test size(ps.bias) == (10, 100)
        @test layer.activation == identity

        layer = Scale(10, 100, relu; bias=false)
        ps, st = Lux.setup(rng, layer)

        @test !haskey(ps, :bias)
        @test layer.activation == relu
    end

    @testset "dimensions" begin
        layer = Scale(10, 5)
        ps, st = Lux.setup(rng, layer)

        @test size(first(Lux.apply(layer, randn(10), ps, st))) == (10, 5)
        @test size(first(Lux.apply(layer, randn(10, 5, 2), ps, st))) == (10, 5, 2)
    end

    @testset "zeros" begin
        @test begin
            layer = Scale(10, 1, identity; init_weight=ones)
            first(Lux.apply(layer, ones(10, 1), Lux.setup(rng, layer)...))
        end == ones(10, 1)

        @test begin
            layer = Scale(10, 1, identity; init_weight=ones)
            first(Lux.apply(layer, ones(10, 2), Lux.setup(rng, layer)...))
        end == ones(10, 2)

        @test begin
            layer = Scale(2, identity; init_weight=ones, init_bias=ones)
            first(Lux.apply(layer, [1 2; 3 4], Lux.setup(rng, layer)...))
        end == [2.0 3.0; 4.0 5.0]

        @test begin
            layer = Scale(2, tanh; bias=false, init_weight=zeros)
            first(Lux.apply(layer, [1 2; 3 4], Lux.setup(rng, layer)...))
        end == zeros(2, 2)
    end
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVRProjPO

using Test, RLinearAlgebra, LinearAlgebra, Random

@testset "LSVR Projection Partial Orthogonal -- Procedural" begin

    # Supertype
    @test supertype(LinSysVecRowProjPO) == LinSysVecRowProjection

    # Testing context
    Random.seed!(1010)

    A = rand(10,5)
    x = rand(5)
    b = A * x

    z = rand(5)

    rsub = LinSysVecRowProjPO()

    # Verify that Z (storage) is nothing
    @test isnothing(rsub.Z)

    # Verify that Z (storage) is appropriately initialized
    @test let
        RLinearAlgebra.rsubsolve!(rsub, z, (A[1,:], b[1]), 1)

        (length(rsub.Z) == rsub.m) & reduce(&, [length(z)==5 for z in rsub.Z])
    end

    # Verify that solution is orthogonalized against previous m search directions
    @test let
        for i = 1:6
            RLinearAlgebra.rsubsolve!(rsub, z, (A[i,:], b[i]), i)
        end

        norm( (A * z - b)[2:6] ) < 1e-15
    end

    @test let
        for i = 1:8
            RLinearAlgebra.rsubsolve!(rsub, z, (A[i,:], b[i]), i)
        end

        norm((A * z - b)[4:8]) < 1e-15
    end

end

end # End Module

ext: jl


text:
@testset "Packages list correctness" begin
    ds = RDatasets.datasets()
    dp = RDatasets.packages()

    @test dp[:Package] == sort(unique(ds[:Package]))

    # test caching works
    @test dp === RDatasets.packages()
end

max_stars_repo_name: blockchainhelppro/Wadren-Analysis

max_issues_repo_name: blockchainhelppro/Wadren-Analysis

max_forks_repo_name: blockchainhelppro/Wadren-Analysis


text:
module RLEVectors

using Base.Broadcast
using Requires
using Statistics
using StatsBase
using DataFrames
using Reexport
@reexport using Missings

# RLEVector type
export RLEVector,
    FloatRle,
    IntegerRle,
    BoolRle,
    StringRle,
    RLEVectorList,
    rfirst,
    rwidth,
    rlast,
    rvalue,
    nrun,
    similar,
    collect,
    similar,
    starts,
    widths,
    widths!,
    ends,
    values
import Base: show, length, size, first, last, step, convert, similar, collect, isequal, values, copy

# collections
import Base.==
import Base:
    eltype,
    vcat,
    pop!,
    push!,
    popfirst!,
    pushfirst!,
    insert!,
    deleteat!,
    splice!,
    resize!,
    empty!,
    lastindex,
    intersect,
    append!
export eltype,
    vcat,
    pop!,
    push!,
    popfirst!,
    pushfirst!,
    insert!,
    deleteat!,
    splice!,
    resize!,
    growat!,
    empty!,
    lastindex,
    intersect,
    append!
export deleterun!, decrement_run!

# indexing
import Base.Broadcast: BroadcastStyle, Broadcasted
import Base: getindex, setindex!
import Base: iterate
export getindex,
    setindex!, ind2run, setrun!, ind2runcontext, RLERangesIterator, eachrange, tapply, iterate

# describe
import Base: isempty
export endtype

# group_generics
import Base: broadcast, map
import Base: in, indexin, findmin, findmax
import Base: range, any, all, sum, eltype, unique, minimum, maximum, extrema, first, last
export range, any, all, eltype, sum, unique, minimum, maximum, extrema, first, last
export in, indexin, findmin, findmax
import Statistics: mean, median
export mean, median

# math
import StatsBase: mode, countmap
export mode, countmap

# ranges
export disjoin, disjoin_length, ree, ree!, inverse_ree, numruns, rangeMeans

# utils
export rep

# sorting
import Base.Order: Ordering
import Base.Sort: QuickSortAlg
import Base: sort, sort!, issorted, reverse, reverse!, sortperm, Algorithm
export sort, sort!, issorted, reverse, reverse!, sortperm, permute_runs

# DataFrames
import DataFrames: AbstractDataFrame, DataFrame, Index, index, nrow, ncol
export RLEDataFrame, nrow, ncol, columns, index, names
export rowSums, rowMeans, rowMedians, colSums, colMeans, colMedians

# Missings
import Missings: allowmissing, disallowmissing
export allowmissing, disallowmissing

### Includes
include("utils.jl")
include("runs.jl")
include("RLEVector-type.jl")
include("RLEDataFrame-type.jl")
include("ranges.jl")
include("describe.jl")
include("indexing.jl")
include("group_generics.jl")
include("collections_api.jl")
include("math.jl")
include("sorting.jl")
include("missing.jl")
include("precompile.jl")
_precompile_()

function __init__()
    @require RCall = "6f49c342-dc21-5d91-9882-a32aef131414" include("rcall.jl")
end

end # module RLEVectors


text:
import StaticArrays: SArray, SVector

##################################################################################
# Static simplex where vertices are represented by StaticArrays.SVector instances.
##################################################################################

struct SSimplex{dim, T} <: AbstractImmutableSimplex{dim, T}
    vertices::Vector{SVector{dim, T}}

    function SSimplex(pts::Vector{SVector{dim, T}}) where {dim, T}
        if !(length(pts) == dim + 1)
            err = """ The input cannot be converted to a simplex.
            Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
            """
            throw(DomainError(pts, err))
        end

        new{dim, T}(pts)
    end

    function SSimplex(pts::Vector{Vector{T}}) where {T}
        dim = length(pts[1])
        if !(length(pts) == dim + 1)
            err = """ The input cannot be converted to a simplex.
            Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
            """
            throw(DomainError(pts, err))
        end

        [SVector{dim, T}(pt) for pt in pts]
    end
end

export SSimplex


text:
using Adapt
using Test

# custom array type

struct CustomArray{T,N} <: AbstractArray{T,N}
    arr::AbstractArray
end

CustomArray(x::AbstractArray{T,N}) where {T,N} = CustomArray{T,N}(x)
Adapt.adapt_storage(::Type{<:CustomArray}, xs::AbstractArray) = CustomArray(xs)

Base.size(x::CustomArray, y...) = size(x.arr, y...)
Base.getindex(x::CustomArray, y...) = getindex(x.arr, y...)

const val = CustomArray{Float64,2}(rand(2,2))

# basic adaption
@test adapt(CustomArray, val.arr) == val
@test adapt(CustomArray, val.arr) isa CustomArray

# idempotency
@test adapt(CustomArray, val) == val
@test adapt(CustomArray, val) isa CustomArray

# custom wrapper
struct Wrapper{T}
    arr::T
end
Wrapper(x::T) where T = Wrapper{T}(x)
Adapt.adapt_structure(to, xs::Wrapper) = Wrapper(adapt(to, xs.arr))
@test adapt(CustomArray, Wrapper(val.arr)) == Wrapper(val)
@test adapt(CustomArray, Wrapper(val.arr)) isa Wrapper{<:CustomArray}

## base wrappers

@test adapt(CustomArray, (val.arr,)) == (val,)
@test @allocated(adapt(nothing, ())) == 0
@test @allocated(adapt(nothing, (1,))) == 0
@test @allocated(adapt(nothing, (1,2,3,4,5,6,7,8,9,10))) == 0

@test adapt(CustomArray, (a=val.arr,)) == (a=val,)

@test adapt(CustomArray, view(val.arr,:,:)) == view(val,:,:)
@test adapt(CustomArray, view(val.arr,:,:)) isa SubArray{<:Any,<:Any,<:CustomArray}

using LinearAlgebra

@test adapt(CustomArray, val.arr') == val'
@test adapt(CustomArray, val.arr') isa Adjoint{<:Any,<:CustomArray}


text:
using NLopt

# Firm value function if they do not change the markup
function V(para::MParameters,Vprime,μ_,ε)
    @unpack β,γ,W̄,μe,σe = para
    # Calculating this period's markup given this period's shock and last period's markup
    μ = μ_/exp(ε)
    v = W̄^(1-γ)*(μ - 1)*μ^(-γ) + β*expectation_normal(ε'->Vprime(μ,ε'),μe,σe,10)
    return v, μ
end

# Firm value function if they do change the markup
function J(para::MParameters)#,Πprime)
    @unpack β,γ,κ,W̄,μe,σe,mugrid = para
    # Value if the firm does change prices
    J(μo) = -(W̄^(1-γ)*(μo - 1)*μo^(-γ) - κ) #+ β*expectation_normal(ε->Πprime(μo,ε),μe,σe,10))
    Jresult = optimize(J,mugrid[1],mugrid[end])
    return return(μprime = Jresult.minimizer, Π = -Jresult.minimum)
end

J(para)

"""
    find_optimal_policy(αvec,Uf,hf)

Computes the optimal policy given policy fuctions hf and indirect utility Uf
"""
function find_optimal_policy(para::MParameters,Πprime, μ_)
    opt = Opt(:LN_COBYLA, 2)
    lower_bounds!(opt, [0., -Inf])
    upper_bounds!(opt, [20,Inf])
    ftol_rel!(opt,1e-8)

    # Setting up piecewise function
    min_objective!(opt, x)
    inequality_constraint!(opt, (μ,ε)->V(para, Πprime, μ_, ε) - x)
    inequality_constraint!(opt, (μ,ε)->J(para, Πprime) - x)

    min_objective!(opt, (x,g)->-government_welfare(x[1],x[2],αvec,Uf))
    equality_constraint!(opt, (x,g) -> -budget_residual(x[1],x[2],αvec,hf))

    minf,minx,ret = NLopt.optimize(opt, [0.3, 0.3])
    if ret == :FTOL_REACHED
        return minx
    end
end;

basis_x = SplineParams(LinRange(-1,1,5),0,3) #cubic splines along x
basis_y = ChebParams(3,-1,1)#Chebyshev polynomials along y

basis = Basis(basis_x,basis_y)

function optimalpolicy(para::MParameters,μ)
    @unpack β,γ,κ,W̄,μe,σe,mugrid = para
    mu_bounds = [mugrid[1],mugrid[end]]
    # Calculating this period's markup given this period's shock and last period's markup
    #μ = μ_/exp(ε)
    # Value if the firm does not change prices
    V = W̄^(1-γ)*(μ - 1)*μ^(-γ) #+ β*expectation_normal(εprime->Πprime([μ,εprime]),μe,σe,10)
    # Value if the firm does change prices
    J(μo) = - (W̄^(1-γ)*(μo - 1)*μo^(-γ) - κ) #+ β*expectation_normal(εprime->Πprime([μo,εprime]),μe,σe,10))
    Jresult = optimize(J,mu_bounds[1],mu_bounds[2])
    # Choosing whether to keep prices or change prices
    if V >= -Jresult.minimum
        return(μprime = μ, Π = V)
    else
        return(μprime = Jresult.minimizer, Π = -Jresult.minimum)
    end
end;

plot(
    layer(y->optimalpolicy(para,y).μprime,0,3,color=["FuncDtioDn"]),
    layer(y->y,0,3,color=["yellow"]),
            Guide.xlabel("μ"),Guide.ylabel("Optimal Policy"))


text:
# This file defines the Attribute type. This type defines an attribute of
#   personnel members (field in the database).

export Attribute

"""
The `Attribute` type defines an attribute of the personnel members in the simulation -- equivalent to a field in the personnel database.

The type contains the following fields:
* `name::String`: the name of the attribute.
* `possibleValues::Vector{String}`: the vector of possible values the attribute can take.
* `initValues::Vector{String}`: the values that the attribute can take upon initialisation.
* `initValueWeights::Vector{String}`: the weights of each initial value, in the same order as the `initValues` vector.

Two additional fields are used to speed up computations:
* `initValueDist::Categorical`: a categorical distribution object which can be sampled whenever an initial value is needed. This object is created in advance and stored to avoid the overhead of computing the probabilities and creating the distribution object every time an initial value must be sampled.
* `initRNG::MersenneTwister`: the random stream that generates the initial values for the attribute.

Constructor:
```
Attribute( name::String )
```
This constructor generates an `Attribute` object with name `name` and no defined possible values.
"""
mutable struct Attribute

    name::String
    # isOrdinal::Bool
    possibleValues::Vector{String}
    initValues::Vector{String}
    initValueWeights::Vector{Float64}

    initValueDist::Categorical
    initRNG::MersenneTwister

    # Basic constructor.
    function Attribute( name::String; # isOrdinal = false
        )::Attribute

        newAttr = new()
        newAttr.name = name
        # newAttr.isOrdinal = isOrdinal
        newAttr.possibleValues = Vector{String}()
        newAttr.initValues = Vector{String}()
        newAttr.initValueWeights = Vector{Float64}()
        newAttr.initValueDist = Categorical( [1.0] )
        newAttr.initRNG = MersenneTwister()
        return newAttr

    end  # Attribute( name )

end  # type Attribute

# ? Is the `isFixed` field still used anywhere? Not really.
# * `isOrdinal::Bool`: a flag that indicates if the attribute has ordinal values, that is, string values that have a hierarchical ordering.
# If the field `isOrdinal` is `true`, these values must be ordered from lowest to highest in the hierarchy.


text:
"""
    newpage()

Create a new package page.

Arguments:
* `path="page"`: the folder to place the source files of your site; these
                 are the markdown files and configuration files you will need
                 to adjust based on which your website will be generated.
* `overwrite=false`: whether to overwrite the folder if it already exists.
"""
function newpage(; path="page", overwrite=false)
    if isdir(path)
        if !overwrite
            error("Path '$path' already exists, use `overwrite=true` " *
                  "if you wish to remove that folder and start again.")
        else
            rm(path, recursive=true)
        end
    end
    mkdir(path)
    store = joinpath(dirname(dirname(pathof(PkgPage))), "page")
    for obj in readdir(store)
        obj == "DeployPage.yml" && continue
        obj == ".gitattributes" && continue
        src = joinpath(store, obj)
        dst = joinpath(path, obj)
        cp(src, dst)
    end
    # Pkg makes files read-only but we don't want that
    for (root, _, files) ∈ walkdir(path)
        for file in files
            chmod(joinpath(root, file), 0o644)
        end
    end

    # Try placing the `DeployPage.yml`
    name = "DeployPage.yml"
    gapath = mkpath(joinpath(".github", "workflows"))
    src = joinpath(store, name)
    dst = joinpath(gapath, name)
    cp(src, dst, force=overwrite)
    chmod(dst, 0o644)

    # Try placing the `.gitattributes`
    name = ".gitattributes"
    src = joinpath(store, name)
    dst = name
    cp(src, dst, force=overwrite)
    chmod(dst, 0o644)

    return nothing
end

function clever_cd(path)
    # try to be clever, in some circumstances, we might be in the "page"
    # folder where this expects us to be one level up; if that's the case
    # don't move.
    # cases covered:
    #   1. we're at the proper level, cd serve & cd back
    #   2. we're in the folder, don't cd and serve
    #   3. user is somewhere else --> error
    bkpath  = pwd()
    inside  = isfile("config.md") && isdir("_layout")
    outside = !inside && isdir(path)
    if !(inside || outside)
        error("Your current path does not seem to contain content for a " *
              "package page that matches `path=$path`; maybe you forgot to " *
              "call `instantiate`?")
    end
    outside && cd(path)
    return (bkpath, outside)
end

"""
    serve(path="page")

Launch the server to render the content of the `path` folder.
"""
function serve(path="page"; kw...)
    # Piracy to avoid code caching and reduce risks of silly
    # errors and hard-to-understand error messages.
    bkpath, outside = clever_cd(path)
    try
        F.serve(clear=true; kw...)
    catch e
        @show e
    finally
        outside && cd(bkpath)
    end
    return nothing
end

# Should only be called in the deploy (by github-action) or via `publish`
# assumes that `purgecss` is available to NodeJS.
function optimize(; input="page", output="", purge=true, kw...)
    occursin("/", output) &&
        error("No depth allowed in `output`: it should either be an " *
              "empty string or something like `web` or `page`.")
    bkpath, outside = clever_cd(input)

    # Optimization phase
    F.optimize(; minify=false, kw...)

    if purge
        # Purge CSS to decrease bootstrap size massively
        io = IOBuffer()
        run(pipeline(`$(NodeJS.npm_cmd()) root`, stdout=io))
        nodepath = String(take!(io))
        run(`$(strip(nodepath))/purgecss/bin/purgecss.js --css __site/css/bootstrap.min.css --content __site/index.html --output __site/css/bootstrap.min.css`)
    end

    # if required, copy the content of `__site` to a subfolder so that
    # that subfolder can be deployed. It's the user's responsibility to
    # ensure that there isn't already a folder $output.
    if !isempty(output)
        outp = mkpath(joinpath("__site", output))
        for obj in readdir("__site")
            obj == output && continue
            src = joinpath("__site", obj)
            dst = joinpath(outp, obj)
            mv(src, dst; force=true)
        end
    end
    outside && cd(bkpath)
    return nothing
end


text:
using CuYao
using Yao
using BenchmarkTools

sf(x, y) = abs(x-y)
a = randn(1024)
ca = a |> cu
b = randn(1024)
cb = b |> cu
@benchmark expect(StatFunctional{2}(sf), a, b) seconds=2
@benchmark expect(StatFunctional{2}(sf), a) seconds=2
@benchmark expect(StatFunctional{2}(sf), ca, cb) seconds=2
@benchmark expect(StatFunctional{2}(sf), ca) seconds=2


text:
function test_model_state_diff(model)
    n,m = RD.dims(model)
    e = RD.errstate_dim(model)
    x = Vector(rand(model)[1])
    x0 = Vector(rand(model)[1])
    dx = zeros(e)

    RD.state_diff!(model, dx, x, x0)
    @test dx ≈ RD.state_diff(model, x, x0)

    @test (@allocated RD.state_diff!(model, dx, x, x0)) == 0

    # Use static vectors for out-of-place
    x = rand(model)[1]
    x0 = rand(model)[1]
    @test (@allocated RD.state_diff(model, x, x0)) == 0
end

@testset "State Diff ($(typeof(model).name.name))" for model in (Cartpole(), Quadrotor())
    test_model_state_diff(model)
end


text:
@testset "convert_gt(vcfile)" begin
    vcffile = "test.08Jun17.d8b.vcf.gz"
    isfile(vcffile) || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        abspath(joinpath(dirname(pathof(VCFTools)), "..", "test/$vcffile")))
    #@code_warntype convert_gt(UInt8, vcffile; impute = false, center = false, scale = false)
    # @inferred convert_gt(UInt8, vcffile; impute = false, center = false, scale = false)
    # convert to a matrix of UInt8, no impute/center/scale (default)
    @time A = convert_gt(UInt8, vcffile)
    @test size(A) == (191, 1356)
    @test eltype(A) == Union{UInt8, Missing}
    # convert to a matrix of Float64, impute = center = scale = true
    @time B = convert_gt(Float64, vcffile; impute = true, center = true, scale = true)
    @test size(B) == (191, 1356)
    @test eltype(B) == Union{Float64, Missing}
    # convert to a matrix of Float32
    # @code_warntype convert_gt(Float32, vcffile)
    @time C = convert_gt(Float32, vcffile)
    @test size(C) == (191, 1356)
    @test eltype(C) == Union{Float32, Missing}
    # read the first record to a UInt8 vector, no impute/center/scale (default)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    A1 = Vector{Union{UInt8, Missing}}(undef, nsamples(vcffile))
    @time copy_gt!(A1, reader)
    @test all(A1[findall(!ismissing, A1)] .== A[findall(!ismissing, A[:, 1])])
    @test all(findall(ismissing, A1) .== findall(ismissing, A[:, 1]))
    # convert next 5 records to a Float64 matrix, impute = center = scale = true
    B5 = Matrix{Union{Float64, Missing}}(undef, nsamples(vcffile), 5)
    @time copy_gt!(B5, reader; impute = true, center = true, scale = true)
    B_copy = copy(B[:, 2:6])
    @test all(B5[findall(!ismissing, B5)] .== B_copy[findall(!ismissing, B[:, 2:6])])
    @test all(B5[findall(ismissing, B5)] .== B_copy[findall(ismissing, B_copy)])
    # convert next 5 records to a Float32 matrix
    C5 = Matrix{Union{Float64, Missing}}(undef, nsamples(vcffile), 5)
    @time copy_gt!(C5, reader)
    @test findall(x -> x === missing, C5) == []
    close(reader)
end

@testset "convert_ht(vcfile)" begin
    vcffile = "test.08Jun17.d8b.vcf.gz"
    # convert to a reference haplotype panel of Float32
    @time H1 = convert_ht(Float32, vcffile)
    @test size(H1) == (382, 1356)
    @test typeof(H1) == Matrix{Float32}
    # convert to a reference haplotype panel of Float64
    @time H2 = convert_ht(Float64, vcffile)
    @test size(H2) == (382, 1356)
    @test typeof(H2) == Matrix{Float64}
    # convert to bitarray
    Hb = convert_ht(Bool, vcffile)
    @test all(H1 .== convert(Matrix{Float32}, Hb))
    # convert first record into 2 haplotype vectors and check their sum is the genotype vector
    reader = VCF.Reader(openvcf(vcffile, "r"))
    h1h2 = Matrix{Float64}(undef, 2, nrecords(vcffile))
    copy_ht!(h1h2, reader)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    g1 = Matrix{Union{Float64, Missing}}(undef, 1, nrecords(vcffile))
    copy_gt!(g1, reader)
    @test all(sum(h1h2, dims=1) .== g1)
end

@testset "convert_ds" begin
    vcffile = "test.08Jun17.d8b.vcf.gz"
    D = convert_ds(Float64, vcffile, key = "DS", impute=false, center=false, scale=false)
    @test D[1, 5] == 1.0
    @test D[1, 1345] == 0.05
    @test D[2, 13] == 1.75
    @test eltype(D) == Union{Missing, Float64}
    @test size(D) == (191, 1356)
end

@testset "convert to transpose functions" begin
    vcffile = "test.08Jun17.d8b.vcf"

    # convert_gt_trans!
    @time A  = convert_gt(Float64, vcffile)
    @time At = convert_gt(Float64, vcffile, trans=true)
    @test all(At .== A')

    # test if eof(reader) is working
    out = Matrix{Union{Float64, Missing}}(undef, 191, 1400)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_gt!(out, reader)
    @test all(ismissing.(out[:, 1357:end]))

    out = Matrix{Union{Float64, Missing}}(undef, 1400, 191)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_gt_trans!(out, reader)
    @test all(ismissing.(out[1357:end, :]))

    # test impute, center, scale (these use mean and var functions, currently not imported)
    # A = Matrix{Union{Float64, Missing}}(undef, 191, 1400)
    # reader = VCF.Reader(openvcf(vcffile, "r"))
    # copy_gt!(A, reader, impute=true, center=true, scale=true)
    # @test isapprox(mean(A[:, 5]), 0, atol=10)
    # @test isapprox(var(A[:, 5]), 1, atol=10)
    # At = Matrix{Union{Float64, Missing}}(undef, 1400, 191)
    # reader = VCF.Reader(openvcf(vcffile, "r"))
    # copy_gt_trans!(At, reader, impute=true, center=true, scale=true)
    # @test isapprox(mean(At[5, :]), 0, atol=10)
    # @test isapprox(var(At[5, :]), 1, atol=10)

    # convert_ht_trans!
    @time H  = convert_ht(Float64, vcffile)
    @time Ht = convert_ht(Float64, vcffile, trans=true)
    @test all(Ht .== H')

    # test if eof(reader) is working
    out = Matrix{Float64}(undef, 382, 1400)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_ht!(out, reader)
    @test size(out) == (382, 1400)
    @test all(out[:, 1358:end] .== 0.0)

    out = Matrix{Float64}(undef, 1400, 382)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_ht_trans!(out, reader)
    @test size(out) == (1400, 382)
    @test all(out[1358:end, :] .== 0.0)

    # convert_ds_trans!
    @time D  = convert_ds(Float64, vcffile)
    @time Dt = convert_ds(Float64, vcffile, trans=true)
    @test all(Dt .== D')

    # test if eof(reader) is working
    out = Matrix{Union{Missing, Float64}}(undef, 191, 1400)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_ds!(out, reader)
    @test size(out) == (191, 1400)
    @test all(ismissing.(out[:, 1358:end]))

    out = Matrix{Union{Missing, Float64}}(undef, 1400, 191)
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_ds_trans!(out, reader)
    @test size(out) == (1400, 191)
    @test all(ismissing.(out[1358:end, :]))
end

@testset "record info" begin
    vcffile = "test.08Jun17.d8b.vcf"
    X, sampleID, chr, pos, SNPid, ref, alt = convert_gt(Float64, vcffile,
        trans=false, save_snp_info=true, msg = "Importing genotype file...")
    Xt, sampleIDt, chrt, post, SNPidt, reft, altt = convert_gt(Float64, vcffile,
        trans=true, save_snp_info=true, msg = "Importing genotype file...")

    @test all(X .== Xt')
    @test all(sampleID .== sampleIDt)
    @test all(chr .== chrt)
    @test all(pos .== post)
    @test all(SNPid .== SNPidt)
    @test all(ref .== reft)
    @test all(alt .== altt)

    D, sampleID, chr, pos, SNPid, ref, alt = convert_ds(Float64, vcffile,
        trans=false, save_snp_info=true, msg = "Importing genotype file...")
    Dt, sampleIDt, chrt, post, SNPidt, reft, altt = convert_ds(Float64, vcffile,
        trans=true, save_snp_info=true, msg = "Importing genotype file...")

    @test all(D .== Dt')
    @test all(sampleID .== sampleIDt)
    @test all(chr .== chrt)
    @test all(pos .== post)
    @test all(SNPid .== SNPidt)
    @test all(ref .== reft)
    @test all(alt .== altt)
end

@testset "write VCF" begin
    # write haplotypes
    H1 = bitrand(100, 200)
    H2 = bitrand(100, 200)
    write_vcf("test.write.vcf.gz", H1, H2)
    H = convert_ht(Bool, "test.write.vcf.gz")
    @test all(H1 .== view(H, 1:2:size(H, 1), :))
    @test all(H2 .== view(H, 2:2:size(H, 1), :))

    # write genotypes
    X = H1 + H2
    write_vcf("test.write.vcf.gz", X)
    G = convert_gt(Float64, "test.write.vcf.gz")
    @test all(G .== X)
end


text:
struct DSASignature{H,P} <: AbstractSignature
    hash::H # hash could have its own type
    r
    s
    pubkey::P #
end

import Base.mod
mod(G::CyclicGroup) = mod(value(G),order(G))
mod(n::Integer,G::CyclicGroup) = mod(n,order(G))

# https://github.com/zydeon/DSA/blob/master/DSA.py
function modinv_euclidean(z::Integer,a::Integer)
    if 0<z<a
	i = a
	j = z
	y1 = BigInt(1)
	y2 = BigInt(0)
	while j > 0
	    q = div(i,j)
	    r = i-j*q
	    y = y2 - y1*q
	    i, j = j, r
	    y2, y1 = y1, y
	    if i == 1
		return mod(y2,a)
            end
        end
    else
        error("Inverse Error")
    end
end

modinv_fermat(k::Integer,a::Integer) = powermod(k,a-2,a)

modinv(k::Integer,G::CyclicGroup) = modinv_euclidean(k,order(G))

function DSASignature(hash,signer::AbstractSigner)

    h = BigInt(hash)
    x = BigInt(signer.privkey)
    G = signer.G

    k = rand(1:order(G)) ### chooses a number from 0 to q
    r = mod(G^k)
    kinv = modinv(k,G)
    s = mod(kinv*(h + x*r),G)

    # This always works for modinv_fermat
    # s < order(G)/2 || s = order(G) - s

    if s==0
        return DSASignature(hash,signer,G)
    else
        return DSASignature(hash,r,s,signer.pubkey)
    end
end

function verify(sr::DSASignature,G::CyclicGroup)
    r,s = sr.r,sr.s
    h = BigInt(sr.hash)
    Y = typeof(G)(sr.pubkey,G)

    0<r<order(G) || return false
    0<s<order(G) || return false

    sinv = modinv(s,G)
    w = mod(sinv,G)
    u1 = mod(h*w,G)
    u2 = mod(r*w,G)
    v = mod(G^u1 * Y^u2)

    if v==r
        return true
    else
        return false
    end
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite  -  HS114
# This file is JuMP implementation of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# This problem has a quadratic objective with quadratic constraints as well
# as a couple constraints that have fractions of variables.
#############################################################################

using JuMP
using Base.Test

let

n = 10
a = 0.99
b = 0.9

lower = [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 85, 90, 3, 1.2, 145]
upper = [2000, 16000, 120, 5000, 2000, 93, 95, 12, 4, 162]
start = [1745, 12000, 110, 3048, 1974, 89.2, 92.8, 8, 3.6, 145]

m = Model()
@variable(m, lower[i] <= x[i=1:n] <= upper[i], start = start[i])

@NLobjective(m, Min, 5.04*x[1] + .035*x[2] + 10*x[3] + 3.36*x[5] - .063*x[4]*x[7])

@NLconstraint(m, 35.82 - .222*x[10] - b*x[9] >= 0)
@NLconstraint(m, -133 + 3*x[7] - a*x[10] >= 0)
@NLconstraint(m, -(35.82 - .222*x[10] - b*x[9]) + x[9]*(1/b - b) >= 0)
@NLconstraint(m, -(-133 + 3*x[7] - a*x[10]) + (1/a - a)*x[10] >= 0)
@NLconstraint(m, 1.12*x[1] + .13167*x[1]*x[8] - .00667*x[1]*x[8]^2 - a*x[4] >= 0)
@NLconstraint(m, 57.425 + 1.098*x[8] - .038*x[8]^2 + .325*x[6] - a*x[7] >= 0)
@NLconstraint(m, -(1.12*x[1] + .13167*x[1]*x[8] - .00667*x[1]*x[8]^2 - a*x[4]) + (1/a - a)*x[4] >= 0)
@NLconstraint(m, -(57.425 + 1.098*x[8] - .038*x[8]^2 + .325*x[6] - a*x[7]) + (1/a - a)*x[7] >= 0)
@NLconstraint(m, 1.22*x[4] - x[1] - x[5] == 0)
@NLconstraint(m, 98000*x[3]/(x[4]*x[9] + 1000*x[3]) - x[6] == 0)
@NLconstraint(m, (x[2] + x[5])/x[1] - x[8] == 0)

solve(m)

@test_approx_eq_eps getobjectivevalue(m) -1768.80696 1e-3

end


text:
include("Bus.jl")
include("Generator.jl")

export get_area
export get_base_kv
export get_bus
export get_cost_pg_model
export get_cost_pg_parameters
export get_cost_start_cold
export get_cost_start_hot
export get_cost_start_warm
export get_down_time
export get_fuel_type
export get_index
export get_max_run
export get_min_down
export get_min_run
export get_must_run
export get_name
export get_pg
export get_pg_0
export get_pg_lb
export get_pg_ub
export get_pgd_lb
export get_pgd_ub
export get_qg
export get_qg_lb
export get_qg_ub
export get_reserve_agc
export get_reserve_agc_pg_lb
export get_reserve_agc_pg_ub
export get_service_on
export get_status
export get_time_start_cold
export get_time_start_hot
export get_time_start_warm
export get_type
export get_up_time
export get_va
export get_vm
export get_vm_lb
export get_vm_ub
export get_zone

max_stars_repo_name: NREL-SIIP/CommonTransmissionModel

max_issues_repo_name: NREL-SIIP/CommonTransmissionModel

max_forks_repo_name: NREL-SIIP/CommonTransmissionModel


text:
function ilqr_solve!(solver::Solver)

    # printstyled("Iterative LQR\n",
	# 	color=:red, bold=true)

	# data
	policy = solver.policy
    problem = solver.problem
    reset!(problem.model)
    reset!(problem.objective)
	data = solver.data
    solver.options.reset_cache && reset!(data)

	cost!(data, problem,
        mode=:nominal)
    gradients!(problem,
        mode=:nominal)
    backward_pass!(policy, problem,
        mode=:nominal)

    obj_prev = data.objective[1]
    for i = 1:solver.options.max_iterations
        forward_pass!(policy, problem, data,
            min_step_size=solver.options.min_step_size,
            line_search=solver.options.line_search,
            verbose=solver.options.verbose)
        if solver.options.line_search != :none
            gradients!(problem,
                mode=:nominal)
            backward_pass!(policy, problem,
                mode=:nominal)
            lagrangian_gradient!(data, policy, problem)
        end

        # gradient norm
        gradient_norm = norm(data.gradient, Inf)

        # info
        data.iterations[1] += 1
        solver.options.verbose && println(
            "iter:                  $i
             cost:                  $(data.objective[1])
			 gradient_norm:         $(gradient_norm)
			 max_violation:         $(data.max_violation[1])
			 step_size:             $(data.step_size[1])")

        # check convergence
		gradient_norm < solver.options.lagrangian_gradient_tolerance && break
        abs(data.objective[1] - obj_prev) < solver.options.objective_tolerance ? break : (obj_prev = data.objective[1])
        !data.status[1] && break
    end

    return nothing
end

function ilqr_solve!(solver::Solver, states, actions; kwargs...)
    initialize_controls!(solver, actions)
    initialize_states!(solver, states)
    ilqr_solve!(solver; kwargs...)
end

"""
    gradient of Lagrangian
        https://web.stanford.edu/class/ee363/lectures/lqr-lagrange.pdf
"""
function lagrangian_gradient!(data::SolverData, policy::PolicyData, problem::ProblemData)
	p = policy.value.gradient
    Qx = policy.action_value.gradient_state
    Qu = policy.action_value.gradient_action
    H = length(problem.states)

    for t = 1:H-1
        Lx = @views data.gradient[data.indices_state[t]]
        Lx .= Qx[t]
        Lx .-= p[t]
        Lu = @views data.gradient[data.indices_action[t]]
        Lu .= Qu[t]
        # data.gradient[data.indices_state[t]] = Qx[t] - p[t] # should always be zero by construction
        # data.gradient[data.indices_action[t]] = Qu[t]
    end
    # NOTE: gradient wrt x1 is satisfied implicitly
end

"""
    augmented Lagrangian solve
"""
function constrained_ilqr_solve!(solver::Solver; augmented_lagrangian_callback!::Function=x->nothing)

	# verbose && printstyled("Iterative LQR\n",
	# 	color=:red, bold=true)

    # reset solver cache
    reset!(solver.data)

    # reset duals
    for (t, λ) in enumerate(solver.problem.objective.costs.constraint_dual)
        fill!(λ, 0.0)
	end

	# initialize penalty
	for (t, ρ) in enumerate(solver.problem.objective.costs.constraint_penalty)
        fill!(ρ, solver.options.initial_constraint_penalty)
	end

	for i = 1:solver.options.max_dual_updates
		solver.options.verbose && println("  al iter: $i")

		# primal minimization
		ilqr_solve!(solver)

		# update trajectories
		cost!(solver.data, solver.problem,
            mode=:nominal)

        # constraint violation
		solver.data.max_violation[1] <= solver.options.constraint_tolerance && break

        # dual ascent
		augmented_lagrangian_update!(solver.problem.objective.costs,
			scaling_penalty=solver.options.scaling_penalty,
            max_penalty=solver.options.max_penalty)

		# user-defined callback (continuation methods on the models etc.)
		augmented_lagrangian_callback!(solver)
	end

    return nothing
end

function constrained_ilqr_solve!(solver::Solver, states, actions; kwargs...)
    initialize_controls!(solver, actions)
    initialize_states!(solver, states)
    constrained_ilqr_solve!(solver; kwargs...)
end

function solve!(solver::Solver{T,N,M,NN,MM,MN,NNN,MNN,X,U,D,O}, args...; kwargs...) where {T,N,M,NN,MM,MN,NNN,MNN,X,U,D,O<:Objective{T}}
    ilqr_solve!(solver, args...; kwargs...)
end

function solve!(solver::Solver{T,N,M,NN,MM,MN,NNN,MNN,X,U,D,O}, args...; kwargs...) where {T,N,M,NN,MM,MN,NNN,MNN,X,U,D,O<:AugmentedLagrangianCosts{T}}
    constrained_ilqr_solve!(solver, args...; kwargs...)
end


text:
using Test, SpikeSynchrony

# These test are based on the analytic results of van Rossum (2001)
f = collect(1.0:10.0)

# Distance of f with f must be always 0
for τ in 10 .* rand(10)
    @test vanRossum(f, f, τ) == 0.0
    # @test vanRossum_fast(f, f, τ) == 0.0
end

# Distance of two identical vectors but the second has one extra spike,
# it is analytically known to be 1.0 (with our normalization)
# irrespectively of the length of the vectors, for any τ.: eq (7)
g = copy(f); push!(g, 11.0)
for τ in 10 .* rand(10)
    @test vanRossum(f, g, τ) ≈ 1.0
    # @test vanRossum_fast(f, g, 1.0) ≈ 1.0
end

# Shifting a spike from f to g produces a known result: eq. (8)
for δt in (0.25, 0.5)
    for τ in 10 .* rand(10)
        g = copy(f);
        g[end] = f[end] + δt
        # divide by 2 due to our normalization
        @test vanRossum(f, g, τ)^2 / 2 ≈ 1.0 - exp(-δt/τ)
        # @test vanRossum_fast(f, g, τ)^2 / 2 ≈ 1.0 - exp(-δt/τ)
    end
end


text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================

function transientTerm(phi_old::CellValue, dt::Real, alfa::Real)
transientTerm(phi_old, dt, alfa.*ones(Float64,tuple(phi_old.domain.dims...)))
end

function transientTerm(phi_old::CellValue, dt::Real)
transientTerm(phi_old, dt, 1.0)
end

function transientTerm(phi_old::CellValue, dt::Real, alfa::CellValue)
d = phi_old.domain.dimension
if d==1 || d==1.5
  transientTerm1D(phi_old, dt, alfa.value[2:end-1])
elseif d==2 || d==2.5 || d==2.8
  transientTerm2D(phi_old, dt, alfa.value[2:end-1,2:end-1])
elseif d==3 || d==3.2
  transientTerm3D(phi_old, dt, alfa[2:end-1,2:end-1,2:end-1])
end
end

function transientTerm(phi_old::CellValue, dt::Real, alfa::Array{T}) where T<:Real
d = phi_old.domain.dimension
if d==1 || d==1.5
  transientTerm1D(phi_old, dt, alfa)
elseif d==2 || d==2.5 || d==2.8
  transientTerm2D(phi_old, dt, alfa)
elseif d==3 || d==3.2
  transientTerm3D(phi_old, dt, alfa)
end

end

function transientTerm1D(phi_old::CellValue,
		    dt::Real, alfa::Array{T}) where T<:Real
# returns the matrix and RHS for a d(phi)/dt term

# extract data from the mesh structure
Nx = phi_old.domain.dims[1]
G = [1:Nx+2;]

# rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G[2:Nx+1],Nx) # main diagonal (only internal cells)
AP_diag = reshape(alfa/dt,Nx)
M = sparse(row_index, row_index, AP_diag, Nx+2, Nx+2)

# define the RHS Vector
RHS = zeros(Nx+2)

# assign the values of the RHS vector
RHS[row_index] = reshape(alfa.*phi_old.value[2:Nx+1]/dt,Nx)

(M, RHS)

end

function transientTerm2D(phi_old::CellValue,
		    dt::Real, alfa::Array{T}) where T<:Real
# returns the matrix and RHS for a d(phi)/dt term

# extract data from the mesh structure
Nx = phi_old.domain.dims[1]
Ny = phi_old.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)

# rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny) # main diagonal (only internal cells)
AP_diag = reshape(alfa/dt,Nx*Ny)
M = sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))

# define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2))

# assign the values of the RHS vector
RHS[row_index] = reshape(alfa.*phi_old.value[2:Nx+1,2:Ny+1]/dt,Nx*Ny)

(M, RHS)

end

function transientTerm3D(phi_old::CellValue,
		    dt::Real, alfa::Array{T}) where T<:Real
# returns the matrix and RHS for a d(phi)/dt term

# extract data from the mesh structure
Nx = phi_old.domain.dims[1]
Ny = phi_old.domain.dims[2]
Nz = phi_old.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)

# rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz) # main diagonal (only internal cells)
AP_diag = reshape(alfa/dt,Nx*Ny*Nz)
M = sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))

# define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2)*(Nz+2))

# assign the values of the RHS vector
RHS[row_index] = reshape(alfa.*phi_old.value[2:Nx+1,2:Ny+1,2:Nz+1]/dt,Nx*Ny*Nz)

(M, RHS)

end


text:
module Test_SignalAngle

using  Modia3D
import Modia3D.ModiaMath

@signal Signal(;A=5.0, B=2.0) begin
   y = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Signal, sim::ModiaMath.SimulationState)
    signal.y.value  = signal.A*sin(signal.B * sim.time)
end

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly PendulumDrivenKinematically(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx, m=1.0) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))

   # Pendulum
   body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))

   rev    = Modia3D.Revolute(world, frame1, axis = 1)

   sig    = Signal()
   signal = Modia3D.SignalToFlangeAngle(sig.y)
   Modia3D.connect(signal, rev)
end

pendulum = PendulumDrivenKinematically(Lx=1.6, m=0.5, sceneOptions = Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3))
model    = Modia3D.SimulationModel(pendulum, analysis=ModiaMath.KinematicAnalysis)
# ModiaMath.print_ModelVariables(model)
result   = ModiaMath.simulate!(model, stopTime=5.0)
ModiaMath.plot(result, ["sig.y", "rev.phi"] )

println("... success of Test_SignalAngle.jl!")
end


text:
@showbanner "Basic tria membrane test"
st=Structure()
lcset=LoadCaseSet()

add_uniaxial_metal!(st,"steel",2e11,0.3,7849.0474)
add_static_case!(lcset,"DL",0)

add_node!(st,1,0,0,0)
add_node!(st,2,6,0,0)
add_node!(st,3,12,0,0)
add_node!(st,4,18,0,0)
add_node!(st,5,0,6,0)
add_node!(st,6,6,6,0)
add_node!(st,7,12,6,0)
add_node!(st,8,18,6,0)

add_tria!(st,1,1,2,5,"steel",1e-3)
add_tria!(st,2,2,6,5,"steel",1e-3)
add_tria!(st,3,2,3,6,"steel",1e-3)
add_tria!(st,4,3,7,6,"steel",1e-3)
add_tria!(st,5,3,4,7,"steel",1e-3)
add_tria!(st,6,4,8,7,"steel",1e-3)

add_nodal_force!(lcset,"DL",8,0,-1e5,0,0,0,0)

set_nodal_restraint!(st,1,true,true,true,true,true,true)
set_nodal_restraint!(st,5,true,true,true,true,true,true)
for i in [2,3,4,6,7,8]
    set_nodal_restraint!(st,i,false,false,true,true,true,true)
end

assembly=assemble!(st,lcset,path=PATH)

solve(assembly)

r=result_nodal_displacement(assembly,"DL",8)

@test r≈[0.00345715, -0.0164571, 0.0, 0.0, 0.0, 0.0] atol=1e-4


text:
@testset "\e[34mEXAMPLES" begin
    println("\nTesting examples, stdout temporarily disabled")
    @suppress_out begin
        # @test_nothrow include("../examples/layout.jl")

        # @test_nothrow include("../examples/inspect.jl")
        include("../examples/inspect.jl")
        # @test_nothrow include("../examples/logging.jl")
        # include("../examples/logging.jl")
        # @test_nothrow include("../examples/progressbars.jl")
        include("../examples/progressbars.jl")
        # @test_nothrow include("../examples/text_box.jl")
        include("../examples/text_box.jl")
        # @test_nothrow include("../examples/text_panel.jl")
        include("../examples/text_panel.jl")
        # @test_nothrow include("../examples/text_style.jl")
        include("../examples/text_style.jl")
        # @test_nothrow include("../examples/tree.jl")
        include("../examples/tree.jl")
    end
end


text:
### Calculating average driving time: INPUT C, A_drive[N], transition_matrix and OUTPUT A_drive

function averagedrivingtime(C, A_drive, transition_matrix)
    sum_driving_time = 0;
    driving_time_entries = zeros(C,1);
    for t=1:T
        driving_time_entries = transition_matrix[:,t,3];
        sum_driving_time = sum_driving_time + sum(driving_time_entries);
    end
    A_drive = A_drive  + sum_driving_time/(C * T * 60 * 60);
    return A_drive
end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
# io/z-test.jl

using Test

@test caller_name(:test)
@test fld_read(:test)
@test fld_write(:test)
@test ir_dump(:test)
@test loadpfile(:test)
@test prompt(:test)

#@test shows(:test) # no - not exported
@shows ones(3,4)


text:
function solve!(solver::ProjectedNewtonSolver)
    # reset!(solver)
    update_constraints!(solver)
    copy_constraints!(solver)
    copy_multipliers!(solver)
    constraint_jacobian!(solver)
    copy_jacobians!(solver)
    TO.cost_expansion!(solver)
    update_active_set!(solver; tol=solver.opts.active_set_tolerance_pn)
    copy_active_set!(solver)

    if solver.opts.verbose_pn
        println("\nProjection:")
    end
    viol = projection_solve!(solver)
    # copyto!(solver.Z, solver.P)

    # Copy the multipliers back to the ALConSet
    copyback_multipliers!(solver.λ, solver)

    terminate!(solver)
    return solver
end

function projection_solve!(solver::ProjectedNewtonSolver)
    ϵ_feas = solver.opts.constraint_tolerance
    viol = norm(solver.d[solver.active_set], Inf)
    max_projection_iters = solver.opts.n_steps

    count = 0
    while count <= max_projection_iters && viol > ϵ_feas
        viol = _projection_solve!(solver)
        if solver.opts.multiplier_projection
            res = multiplier_projection!(solver)
        else
            res = Inf
        end
        count += 1
        record_iteration!(solver, viol, res)
    end
    return viol
end

function record_iteration!(solver::ProjectedNewtonSolver, viol, res)
    J = TO.cost(solver)
    J_prev = solver.stats.cost[solver.stats.iterations]
    record_iteration!(solver.stats, cost=TO.cost(solver), c_max=viol, is_pn=true,
        dJ=J_prev-J, gradient=res, penalty_max=NaN)
end

function _projection_solve!(solver::ProjectedNewtonSolver)
    Z = primals(solver)
    a = solver.active_set
    max_refinements = 10
    convergence_rate_threshold = solver.opts.r_threshold

    # regularization
    ρ_chol = solver.opts.ρ_chol
    ρ_primal = solver.opts.ρ_primal

    # Assume constant, diagonal cost Hessian (for now)
    H = solver.H

    # Update everything
    update_constraints!(solver)
    constraint_jacobian!(solver)
    update_active_set!(solver; tol=solver.opts.active_set_tolerance_pn)
    TO.cost_expansion!(solver)

    # Copy results from constraint sets to sparse arrays
    copyto!(solver.P, solver.Z)
    copy_constraints!(solver)
    copy_jacobians!(solver)
    copy_active_set!(solver)

    # Get active constraints
    D,d = active_constraints(solver)

    viol0 = norm(d,Inf)
    if solver.opts.verbose_pn
        println("feas0: $viol0")
    end

    if ρ_primal > 0.0
        dim_primal = size(solver.P.Z)[1]
        for i = 1:dim_primal
            H[i,i] += ρ_primal
        end
    end

    if isdiag(H)
        HinvD = Diagonal(H)\D'
    else
        HinvD = H \ Matrix(D')  # TODO: find a better way to do this
    end

    S = Symmetric(D*HinvD)
    Sreg = cholesky(S + ρ_chol*I) #TODO is this fast or slow? try above
    viol_prev = viol0
    count = 0
    while count < max_refinements
        viol = _projection_linesearch!(solver, (S,Sreg), HinvD)
        convergence_rate = log10(viol) / log10(viol_prev)
        viol_prev = viol
        count += 1

        if solver.opts.verbose_pn
            println("conv rate: $convergence_rate")
        end

        if convergence_rate < convergence_rate_threshold ||
                       viol < solver.opts.constraint_tolerance
            break
        end
    end
    copyto!(solver.Z, solver.P)
    return viol_prev
end

function _projection_linesearch!(solver::ProjectedNewtonSolver,
        S, HinvD)
    conSet = get_constraints(solver)
    a = solver.active_set
    d = solver.d[a]
    viol0 = norm(d,Inf)
    viol = Inf

    P = solver.P
    Z = solver.Z
    P̄ = solver.P̄
    Z̄ = solver.Z̄

    solve_tol = 1e-8
    refinement_iters = 25
    α = 1.0
    ϕ = 0.5
    count = 1
    while true
        δλ = reg_solve(S[1], d, S[2], solve_tol, refinement_iters)
        δZ = -HinvD*δλ
        P̄.Z .= P.Z + α*δZ

        copyto!(Z̄, P̄)
        update_constraints!(solver, Z̄)
        TO.max_violation!(conSet)
        viol_ = maximum(conSet.c_max)
        copy_constraints!(solver)
        d = solver.d[a]
        viol = norm(d,Inf)

        if solver.opts.verbose_pn
            println("feas: ", viol, " (α = ", α, ")")
        end
        if viol < viol0 || count > 10
            break
        else
            count += 1
            α *= ϕ
        end
    end
    copyto!(P.Z, P̄.Z)
    # copyto!(solver.Z, P.Z)
    return viol
end

reg_solve(A, b, reg::Real, tol=1e-10, max_iters=10) = reg_solve(A, b, A + reg*I, tol, max_iters)
function reg_solve(A, b, B, tol=1e-10, max_iters=10)
    x = B\b
    count = 0
    while count < max_iters
        r = b - A*x
        # println("r_norm = $(norm(r))")

        if norm(r) < tol
            break
        else
            x += B\r
            count += 1
        end
    end
    # println("iters = $count")

    return x
end

function active_constraints(solver::ProjectedNewtonSolver)
    return solver.D[solver.active_set, :], solver.d[solver.active_set]  # this allocates
end

function TO.cost_expansion!(solver::ProjectedNewtonSolver)
    Z = get_trajectory(solver)
    E = solver.E
    obj = get_objective(solver)
    init = !solver.opts.reuse_jacobians
    TO.cost_expansion!(E, obj, Z, init=init)

    xinds, uinds = solver.P.xinds, solver.P.uinds
    H = solver.H
    g = solver.g
    copy_expansion!(H, g, E, xinds, uinds)
    return nothing
end

function copy_expansion!(H, g, E, xinds, uinds)
    N = length(E)

    for k = 1:N-1
        H[xinds[k],xinds[k]] .= E[k].Q
        H[uinds[k],uinds[k]] .= E[k].R
        H[uinds[k],xinds[k]] .= E[k].H
        g[xinds[k]] .= E[k].q
        g[uinds[k]] .= E[k].r
    end
    H[xinds[N],xinds[N]] .= E[N].Q
    g[xinds[N]] .= E[N].q
    return nothing
end

function multiplier_projection!(solver::ProjectedNewtonSolver)
    # λ = view(solver.λ,solver.active_set)
    λ = solver.λ[solver.active_set]
    D,d = active_constraints(solver)
    g = solver.g
    res0 = g + D'λ
    A = D*D'
    Areg = A + I*solver.opts.ρ_primal
    b = D*res0
    δλ = -reg_solve(A, b, Areg)
    λ += δλ
    res = g + D'λ  # primal residual
    solver.λ[solver.active_set] = λ
    return norm(res)
end

function primal_residual(solver::ProjectedNewtonSolver, update::Bool=false)
    if update
        update_constraints!(solver)
        copy_constraints!(solver)
        update_active_set!(solver; tol=solver.opts.active_set_tolerance_pn)
        copy_active_set!(solver)
        constraint_jacobian!(solver)
        copy_jacobians!(solver)
        TO.cost_expansion!(solver)
    end
    λ = solver.λ[solver.active_set]
    D,d = active_constraints(solver)
    g = solver.g
    return norm(D'λ + g)
end

@inline copy_constraints!(solver::ProjectedNewtonSolver) = copy_constraints!(solver.d, solver)
@inline copy_multipliers!(solver::ProjectedNewtonSolver) = copy_multipliers!(solver.λ, solver)
@inline copy_jacobians!(solver::ProjectedNewtonSolver) = copy_jacobians!(solver.D, solver)
@inline copy_active_set!(solver::ProjectedNewtonSolver) = copy_active_set!(solver.active_set, solver)


text:
# ----------
# Geometry Module, providing Geometry convertion function
# ----------

# xyz Coordinetes -> rθz Coordinates
function xy_to_polar(ux,uy;Lx=2π,Ly=Lx,T=Float32)
  nx,ny,nz = size(ux);
  dev = CPU();
  grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)
  Ur,Uθ = xy_to_polar(ux,uy,grid;Lx=2π,Ly=Lx,T=Float32);
  return Ur,Uθ;
end

function xy_to_polar(ux::Array,uy::Array,grid;Lx=2π,Ly=Lx,T=Float32)
#=
  Function for converting x-y vector to r-θ vector, using linear transform
    [x']  =  [cos(θ) -rsin(θ)][r']
    [y']     [sin(θ)  rcos(θ)][θ']
    So e_r =  cosθ ̂i + sinθ ̂j
       e_θ = -sinθ ̂j + cosθ ̂j
=#
  nx,ny,nz = size(ux);
  Ur = zeros(T,nx,ny,nz);
  Uθ = zeros(T,nx,ny,nz);
  for j ∈ 1:ny, i ∈ 1:nx
    r = sqrt(grid.x[i]^2+grid.y[j]^2);
    θ = atan(grid.y[j],grid.x[i]) ;
    θ = isnan(θ) ? π/2 : θ;
    sinθ = sin(θ);
    cosθ = cos(θ);
    Uθ[i,j,:] .= @. -sinθ*ux[i,j,:] + cosθ*uy[i,j,:];
    Ur[i,j,:] .= @.  cosθ*ux[i,j,:] + sinθ*uy[i,j,:];
  end
  return Ur,Uθ;
end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
# utility for softmax
struct OneHotVector <: AbstractVector{Bool}
    ix::UInt32
    of::UInt32
end

Base.size(xs::OneHotVector) = (Int64(xs.of),)
Base.getindex(xs::OneHotVector, i::Integer) = i == xs.ix
Base.getindex(xs::OneHotVector, ::Colon) = OneHotVector(xs.ix, xs.of)

function onehot(l, labels)
    i = something(findfirst(isequal(l), labels), 0)
    i > 0 || error("Value $l is not in labels")
    OneHotVector(i, length(labels))
end

# linear
function update_grads!(::Linear, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L}
    @inbounds for i in eachindex(δ)
        δ[i] = SVector(2 * (pred[i][1] - target[i]) * 𝑤[i][1])
        δ²[i] = SVector(2 * 𝑤[i][1])
    end
end

# logistic - on linear predictor
function update_grads!(::Logistic, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L}
    ϵ = T(2e-7)
    @inbounds for i in eachindex(δ)
        p = clamp(sigmoid(pred[i][1]), ϵ, 1-ϵ)
        δ[i] = SVector((p - target[i]) * 𝑤[i][1])
        δ²[i] = SVector(max(ϵ, p * (1 - p)) * 𝑤[i][1])
    end
end

# Poisson
function update_grads!(::Poisson, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L}
    @inbounds for i in eachindex(δ)
        δ[i] = SVector((exp(pred[i][1]) .- target[i]) * 𝑤[i][1])
        δ²[i] = SVector(exp(pred[i][1]) * 𝑤[i][1])
    end
end

# L1
function update_grads!(::L1, α::T, pred::Vector{SVector{L,T}}, target::AbstractArray{T,1}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L}
    @inbounds for i in eachindex(δ)
        δ[i] =  SVector((α * max(target[i] - pred[i][1], 0) - (1 - α) * max(pred[i][1] - target[i], 0)) * 𝑤[i][1])
    end
end

# Softmax
function update_grads!(::Softmax, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{S}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L,S <: Integer}
    ϵ = T(2e-7)
    @inbounds for i in 1:size(pred, 1)
        # pred[i] = SVector{L,T}(pred[i] .- maximum(pred[i]))
        # p = exp.(pred[i])
        p = exp.(pred[i] .- maximum(pred[i]))
        p = clamp.(p / sum(p), ϵ, 1-ϵ)
        # δ[i] = SVector{L,T}(min.(1 - ϵ, max.(ϵ, p - (onehot(target[i], 1:L)))) .* 𝑤[i][1])
        # δ²[i] = SVector{L,T}(min.(1 - ϵ, max.(ϵ, 2 .* p .* (1 .- p))) .* 𝑤[i][1])
        δ[i] = SVector{L,T}(p - (onehot(target[i], 1:L)) .* 𝑤[i][1])
        δ²[i] = SVector{L,T}(max.(ϵ, p .* (1 .- p)) .* 𝑤[i][1])
        # δ²[i] = SVector{L,T}(p .* (1 .- p) .* 𝑤[i][1])

        # δ[i] = SVector{L,T}((exp.(pred[i]) / sums - (onehot(target[i], 1:L))) .* 𝑤[i][1])
        # δ²[i] = SVector{L,T}(1 / sums .* (1 .- exp.(pred[i]) / sums) .* 𝑤[i][1])
    end
end

# Quantile
function update_grads!(::Quantile, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L}
    @inbounds for i in eachindex(δ)
        δ[i] = target[i] > pred[i][1] ? SVector(α * 𝑤[i][1]) : SVector((α - 1) * 𝑤[i][1])
        δ²[i] = SVector(target[i] - pred[i][1]) # δ² serves to calculate the quantile value - hence no weighting on δ²
    end
end

# Gaussian - http://jrmeyer.github.io/machinelearning/2017/08/18/mle.html
# pred[i][1] = μ
# pred[i][2] = log(σ)
function update_grads!(::Gaussian, α, pred::Vector{SVector{L,T}}, target::AbstractArray{T,1}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat,L}
    ϵ = T(2e-7)
    @inbounds @threads for i in eachindex(δ)
        δ[i] = SVector((pred[i][1] - target[i]) / max(ϵ, exp(2 * pred[i][2])) * 𝑤[i][1], (1 - (pred[i][1] - target[i])^2 / max(ϵ, exp(2 * pred[i][2]))) * 𝑤[i][1])
        δ²[i] = SVector(𝑤[i][1] / max(ϵ, exp(2 * pred[i][2])),  2 * 𝑤[i][1] / max(ϵ, exp(2 * pred[i][2])) * (pred[i][1] - target[i])^2)
    end
end

# utility functions
function logit(x::AbstractArray{T,1}) where T <: AbstractFloat
    logit.(x)
end

function logit(x::T) where T <: AbstractFloat
    x = log(x / (1 - x))
    return x
end

function sigmoid(x::AbstractArray{T,1}) where T <: AbstractFloat
    sigmoid.(x)
end

function sigmoid(x::T) where T <: AbstractFloat
    x = 1 / (1 + exp(-x))
    return x
end

function softmax(x::AbstractVector{T}) where T <: AbstractFloat
    x .-= maximum(x)
    x = exp.(x) ./ sum(exp.(x))
    return x
end

##############################
# get the gain metric
##############################
# GradientRegression
function get_gain(::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: GradientRegression,T <: AbstractFloat,L}
    gain = sum((∑δ.^2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
return gain
end

# MultiClassRegression
function get_gain(::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: MultiClassRegression,T <: AbstractFloat,L}
    gain = sum((∑δ.^2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
    return gain
end

# L1 Regression
function get_gain(::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: L1Regression,T <: AbstractFloat,L}
    gain = sum(abs.(∑δ))
    return gain
end

# QuantileRegression
function get_gain(::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: QuantileRegression,T <: AbstractFloat,L}
    gain = sum(abs.(∑δ) ./ (1 .+ λ))
    return gain
end

# GaussianRegression
function get_gain(::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: GaussianRegression,T <: AbstractFloat,L}
    gain = sum((∑δ.^2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
    return gain
end


text:
using ADCME
using PyPlot
using Random
using MAT
using DelimitedFiles
matplotlib.use("agg")

include("fwi_data.jl")
reset_default_graph()
tf.set_random_seed(123)
tf.random.set_random_seed(123)

R_ = constant(R)

idx_var = pml+Int(round(0.1/Δx))

inv_sigmoid(y) = @. - log(1.0/y - 1.0)

case = "fwi"

## fwi
if case == "fwi"
if model == "step"
    var= Variable(zeros(N-1-idx_var)) + inv_sigmoid(1.3/2.0)
    var = -2.0*sigmoid(var)
elseif model == "slop"
    var = Variable(zeros(N-1-idx_var))
end
end

## FC
if case == "nnfwi"
x0 = collect( range(0, 1, length=N-1-idx_var) )
init_guess = fc_init([1, 30, 30, 30, 1])
θ = Variable(init_guess)
if model == "step"
    var= fc(x0, [30, 30, 30, 1], θ) * 0.1 + inv_sigmoid(1.3/2.0)
    var = -2.0*sigmoid(var)
elseif model == "slop"
    var= fc(x0, [30, 30, 30, 1], θ) * 0.1
end
end

## CNN
# x = constant(randn(1,4,1))
# x = tf.keras.layers.Conv1D(8, 3, padding="same")(x)
# x = tf.keras.layers.LeakyReLU(alpha=0.1)(x)
# x = tf.keras.layers.UpSampling1D(4)(x)
# x = tf.keras.layers.Conv1D(8, 3, padding="same")(x)
# x = tf.keras.layers.LeakyReLU(alpha=0.1)(x)
# x = tf.keras.layers.UpSampling1D(4)(x)
# x = tf.keras.layers.Conv1D(8, 3, padding="same")(x)
# x = tf.keras.layers.LeakyReLU(alpha=0.1)(x)
# x = tf.keras.layers.UpSampling1D(4)(x)
# x = tf.keras.layers.Conv1D(1, 3, padding="same")(x)
# x = squeeze(x)
# var = tf.slice(x, (length(x) - (N-1-idx_var),).÷2, (N-1-idx_var,)) * 0.1 - 1.0
# var = (3.5 + 0.5) * sigmoid(var) - 0.5

cH = scatter_add(constant(cH_init), idx_var+1:N-1, var)

cE = (cH[1:end-1]+cH[2:end])/2

reg_TV = sum(abs(cH[2:end]-cH[1:end-1]))
reg_smooth = sum((cH[1:end-2] + cH[3:end] - 2*cH[2:end-1])^2)

function condition(i, E_arr, H_arr)
    i<=NT+1
end

function body(i, E_arr, H_arr)
    E = read(E_arr, i-1)
    H = read(H_arr, i-1)

    ΔH = cH * (E[2:end]-E[1:end-1])/Δx - σH*H
    H += ΔH * Δt
    ΔH = 1/(24Δx) * cH[2:N-2] * (-E[4:end] + 3E[3:end-1] - 3E[2:end-2] + E[1:end-3])
    H = scatter_add(H, 2:N-2, ΔH * Δt)

    ΔE = cE * (H[2:end]-H[1:end-1])/Δx - σE[2:end-1]*E[2:end-1] + R_[i] * Z
    E = scatter_add(E, 2:N-1, ΔE * Δt)
    ΔE = 1/(24Δx) * cE[2:end-1] * (-H[4:end] + 3H[3:end-1] - 3H[2:end-2] + H[1:end-3])
    E = scatter_add(E, 3:N-2, ΔE*Δt)

    i+1, write(E_arr, i, E), write(H_arr, i, H)
end

E_arr = TensorArray(NT+1)
H_arr = TensorArray(NT+1)

E_arr = write(E_arr, 1, zeros(N))
H_arr = write(H_arr, 1, zeros(N-1))

i = constant(2, dtype = Int32)

_, E, H = while_loop(condition, body, [i, E_arr, H_arr])

E = stack(E); E = set_shape(E, (NT+1, N))
H = stack(H)

loss = sum((E[:, idx_rcv] - E_true[:, idx_rcv])^2 * [1.0, 1e-3])
loss1 = sum((E[:, idx_rcv[1]] - E_true[:, idx_rcv[1]])^2 * 1.0)
loss2 = sum((E[:, idx_rcv[2]] - E_true[:, idx_rcv[2]])^2 * 1e-3)

if case == "fwi"
    # loss += 1e3*reg_smooth + 0.0*reg_TV
    # case = string(case, "_reg_smooth")
    loss += 1e2*reg_smooth + 1e2*reg_TV
    case = string(case, "_reg1e2_smooth_TV")
end

if !ispath(joinpath("figures", model, case))
    mkpath(joinpath("figures", model, case))
end
if !ispath(joinpath("results", model, case))
    mkpath(joinpath("results", model, case))
end

function cb(vs, iter, loss)
    if (mod(iter, 20) == 1)
        E_est = vs[1]
        cH_est = vs[2]
        @info iter, vs[3], vs[4]
        close("all")
        figure(figsize = (10,6))
        subplot(221)
        plot(xH, cH_true, "C0", label="True")
        plot(xH, cH0, "k--", label="Initial")
        plot(xH, cH_est, "C3--", label="Estimated")
        legend()
        # xlim([0,1])
        autoscale(enable=true, axis="x", tight=true)
        xlabel("Location (km)")
        ylabel("Velocity (km/s)")
        subplot(222)
        plot((0:NT)*Δt, E_true[:,idx_rcv[1]], "C0", label="True #1",  alpha=0.8)
        plot((0:NT)*Δt, E_true[:,idx_rcv[2]], "C2", label="True #2", alpha=0.8)
        plot((0:NT)*Δt, E_est[:,idx_rcv[1]], "C3--",  label="Estimated #1", alpha=0.8)
        plot((0:NT)*Δt, E_est[:,idx_rcv[2]], "C1--",  label="Estimated #2", alpha=0.8)
        legend(loc="upper right")
        gca().ticklabel_format(axis="y", style="sci", scilimits=(0,0))
        autoscale(enable=true, axis="x", tight=true)
        xlabel("Time (s)")
        ylabel("Amplitude")
        subplot(223)
        pcolormesh(xE, (0:NT)*Δt, E_true, cmap="seismic", vmax=maximum(E_true)/3.0, vmin=-maximum(E_true)/3.0, rasterized=true)
        plot(xE[idx_src]*ones(NT+1), (0:NT)*Δt, "k", label="Source", alpha=0.8)
        plot(xE[idx_rcv[1]]*ones(NT+1), (0:NT)*Δt, "C0--", label="Receiver #1", alpha=0.8)
        plot(xE[idx_rcv[2]]*ones(NT+1), (0:NT)*Δt, "C2--", label="Receiver #2", alpha=0.8)
        # xlim([0,1])
        autoscale(enable=true, axis="x", tight=true)
        legend(loc="upper center")
        xlabel("Localtion (km)")
        ylabel("Time (s)")
        title("True")
        subplot(224)
        pcolormesh(xE, (0:NT)*Δt, E_est, cmap="seismic", vmax=maximum(E_true)/3.0, vmin=-maximum(E_true)/3.0, rasterized=true)
        plot(xE[idx_src]*ones(NT+1), (0:NT)*Δt, "k", label="Source", alpha=0.8)
        plot(xE[idx_rcv[1]]*ones(NT+1), (0:NT)*Δt, "C0--", label="Receiver #1", alpha=0.8)
        plot(xE[idx_rcv[2]]*ones(NT+1), (0:NT)*Δt, "C2--", label="Receiver #2", alpha=0.8)
        # xlim([0,1])
        autoscale(enable=true, axis="x", tight=true)
        legend(loc="upper center")
        xlabel("Localtion (km)")
        ylabel("Time (s)")
        title("Estimated")

        # subplot(325)
        # plot((0:NT)*Δt, E_true[:,idx_rcv[1]].-E_est[:,idx_rcv[1]], "C0", label="Diff #1",  alpha=0.8)
        # legend(loc="upper right")
        # gca().ticklabel_format(axis="y", style="sci", scilimits=(0,0))
        # autoscale(enable=true, axis="x", tight=true)
        # xlabel("Time (s)")
        # ylabel("Amplitude")
        # subplot(326)
        # plot((0:NT)*Δt, E_true[:,idx_rcv[2]].-E_est[:,idx_rcv[2]], "C3", label="Diff #2", alpha=0.8)
        # legend(loc="upper right")
        # gca().ticklabel_format(axis="y", style="sci", scilimits=(0,0))
        # autoscale(enable=true, axis="x", tight=true)
        # xlabel("Time (s)")
        # ylabel("Amplitude")

        tight_layout()
        savefig(joinpath("figures", model, case, "fwi$(lpad(iter,5,"0")).png"), bbox_inch="tight")
        if (mod(iter, 200) == 1)
        matwrite(joinpath("results", model, case, "fwi$(lpad(iter,5,"0")).mat"),
                 Dict("idx_rcv"=>idx_rcv, "idx_src"=>idx_src, "xE"=>xE,
                      "cH0"=>cH0, "cH_true"=>cH_true, "cH_est" => cH_est,
                      "E_true" => E_true, "E_est" => E_est))
        # writedlm(joinpath("results", model, case, "fwi$(lpad(iter,5,"0")).txt"), E_est)
        end
    end
end

lr = placeholder(Float64, shape=[])
opt = AdamOptimizer(learning_rate=lr).minimize(loss)
# opt = tf.train.AdagradOptimizer(learning_rate=0.01).minimize(loss)
# opt = tf.train.GradientDescentOptimizer(learning_rate=0.01).minimize(loss)
# opt = tf.train.RMSPropOptimizer(learning_rate=0.01).minimize(loss)

sess = Session(); init(sess)
cH0 = run(sess, cH)
cb(run(sess, [E, cH, loss1, loss2]), -1, run(sess, loss))
BFGS!(sess, loss, vars = [E, cH, loss1, loss2], callback = cb)

# for iter = 0:10000
#     _, l = run(sess, [opt, loss], feed_dict=Dict(lr=>(randexp()*0.001+0.00001)))
#     cb(run(sess, [E, cH]), iter, l)
#     if mod(iter, 20) == 1
#         @info iter, l
#     end
# end


text:
"""
Delegates loading functions & traits for the EXPRES spectrograph
Author: Eric Ford and collaborators
Created: August 2020
"""

"""
Module providing types and traits and customized functions for the EXPRES Spectrograph.
- [EXPRES Blog](http://exoplanets.astro.yale.edu/expresBlog/)
- ["EXPRES: a next generation RV spectrograph in the search for earth-like worlds" (Jurgenson et al. 2016)](https://ui.adsabs.harvard.edu/abs/2016SPIE.9908E..6TJ/abstract)
"""
module EXPRES
using RvSpectMLBase
import RvSpectMLBase: AbstractInstrument, AbstractInstrument1D, AbstractInstrument2D
import ..EchelleInstruments
import ..EchelleInstruments: read_manifest, read_metadata_from_fits, check_metadata_fields_expected_present, make_manifest
#import ..EchelleInstruments: read_header
import ..EchelleInstruments: default_Δv_to_avoid_tellurics

using DataFrames, Query, FITSIO
using Dates  # If need to use datetime2julian() to get jd.  Need to check about getting BJD.
using Polynomials  # For smoothing "blaze" that includes quantum efficiency variations
using Missings

""" Trait for 1D Extracted spectra from EXPRES """
struct EXPRES1D <: AbstractInstrument1D end

""" Trait for 2D Extracted spectra from EXPRES """
struct EXPRES2D <: AbstractInstrument2D end

"Trait to specify any 1D or 2D Extracted spectra from NEID "
const AnyEXPRES = Union{EXPRES1D,EXPRES2D}

export EXPRES1D, EXPRES2D, AnyEXPRES

# traits.jl imports from RvSpectMBase on its own
include("traits.jl")
export min_order, max_order, min_pixel_in_order, max_pixel_in_order
export orders_to_use_default, min_col_default, max_col_default
export metadata_symbols_default, metadata_strings_default
export default_ccf_mask_v_width
export get_inst_module

include("io.jl")
# not exported, since don't have parameter that specializes to instrument

#import RvSpectMLBase: continuum_normalize_spectrum, continuum_normalize_spectra
#import RvSpectMLBase: filter_line_list, find_worst_telluric_in_each_chunk
import RvSpectMLBase:  make_clean_line_list_from_tellurics, choose_obs_idx_for_init_guess
include("util.jl")
export continuum_normalize_spectrum!, continuum_normalize_spectra!
export filter_line_list, find_worst_telluric_in_each_chunk
export make_clean_line_list_from_tellurics
#export choose_obs_idx_for_init_guess

end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
# ===============================
# Written by AAE
# TU Delft, Winter 2014
# simulkade.com
# ===============================

function fluxLimiter(flName::ASCIIString)
# This function returns a function handle to a flux limiter of user's
# choice.
# available flux limiters are: 'CHARM', 'HCUS', 'HQUICK', 'VanLeer',
# 'VanAlbada1', 'VanAlbada2', 'MinMod', 'SUPERBEE', 'Sweby', 'Osher',
# 'Koren', 'smart', 'MUSCL', 'QUICK', 'MC', and 'UMIST'.
# Default limiter is 'SUPERBEE'. See:
# <http://en.wikipedia.org/wiki/Flux_limiter>

# find the flux limiter function

if flName=="CHARM"
  FL(r)=((r>0.0).*r.*(3.0*r+1.0)./(((r+1.0).^2.0)+eps()*(r==-1.0)))
elseif flName=="HCUS"
  FL(r)=(1.5*(r+abs(r))./(r+2.0))
elseif flName=="HQUICK"
  FL(r)=(2.0*(r+abs(r))./((r+3.0)+eps()*(r==-3.0)))
elseif flName=="ospre"
  FL(r)=((1.5*r.*(r+1.0))./(r.*(r+1.0)+1.0+eps()*((r.*(r+1.0)+1.0)==0.0)))
elseif flName=="VanLeer"
  FL(r)=((r+abs(r))./(1.0+abs(r)))
elseif flName=="VanAlbada1"
  FL(r)=((r+r.*r)./(1.0+r.*r))
elseif flName=="VanAlbada2"
  FL(r)=(2.0*r./(1+r.*r))
elseif flName=="MinMod"
  FL(r)=((r>0.0).*min(r,1.0))
elseif flName=="SUPERBEE"
  FL(r)=(max(0.0, max(min(2.0*r,1.0), min(r,2.0))))
elseif flName=="Osher"
  b=1.5
  FL(r)=(max(0.0, min(r,b)))
elseif flName=="Sweby"
  b=1.5
  FL(r)=(max(0.0, max(min(b*r,1.0), min(r,b))))
elseif flName=="smart"
  FL(r)=(max(0.0, min(4.0,min(0.25+0.75*r, 2.0*r))))
elseif flName=="Koren"
  FL(r)=(max(0.0, min(2.0*r, min((1.0+2.0*r)/3.0, 2.0))))
elseif flName=="MUSCL"
  FL(r)=(max(0.0, min(2.0*r, min(0.5*(1+r), 2.0))))
elseif flName=="QUICK"
  FL(r)=(max(0.0, min(2.0, min(2.0*r, (3.0+r)/4.0))))
elseif flName=="UMIST"
  FL(r)=(max(0.0, min(2.0, min(2.0*r, min((1.0+3.0*r)/4.0, (3.0+r)/4.0)))))
else
  println("The flux limiter of your choice is not available. The SUPERBEE flux limiter is used instead.")
  FL(r)=(max(0.0, max(min(2.0*r,1.0), min(r,2.0))))
end

end

function faceEval(f::Function, x::FaceValue)
FaceValue(x.domain,
    f(x.xvalue),
    f(x.yvalue),
    f(x.zvalue))
end

function faceEval(f::Function, x1::FaceValue, x2::FaceValue)
FaceValue(x1.domain,
    f(x1.xvalue, x2.xvalue),
    f(x1.yvalue, x2.yvalue),
    f(x1.zvalue, x2.zvalue))
end

function faceEval(f::Function, x1::FaceValue, x2::FaceValue, x3::FaceValue)
FaceValue(x1.domain,
    f(x1.xvalue, x2.xvalue, x3.xvalue),
    f(x1.yvalue, x2.yvalue, x3.yvalue),
    f(x1.zvalue, x2.zvalue, x3.zvalue))
end

function faceEval(f::Function, x1::FaceValue, x2::FaceValue, x3::FaceValue, x4::FaceValue)
FaceValue(x1.domain,
    f(x1.xvalue, x2.xvalue, x3.xvalue, x4.xvalue),
    f(x1.yvalue, x2.yvalue, x3.yvalue, x4.yvalue),
    f(x1.zvalue, x2.zvalue, x3.zvalue, x4.zvalue))
end

function faceEval(f::Function, x1::FaceValue, x2::FaceValue, x3::FaceValue, x4::FaceValue, x5::FaceValue)
FaceValue(x1.domain,
    f(x1.xvalue, x2.xvalue, x3.xvalue, x4.xvalue, x5.xvalue),
    f(x1.yvalue, x2.yvalue, x3.yvalue, x4.yvalue, x5.yvalue),
    f(x1.zvalue, x2.zvalue, x3.zvalue, x4.zvalue, x5.zvalue))
end

function faceEval(f::Function, x1::FaceValue, x2::FaceValue, x3::FaceValue, x4::FaceValue, x5::FaceValue, x6::FaceValue)
FaceValue(x1.domain,
    f(x1.xvalue, x2.xvalue, x3.xvalue, x4.xvalue, x5.xvalue, x6.xvalue),
    f(x1.yvalue, x2.yvalue, x3.yvalue, x4.yvalue, x5.yvalue, x6.yvalue),
    f(x1.zvalue, x2.zvalue, x3.zvalue, x4.zvalue, x5.zvalue, x6.zvalue))
end

function cellEval(f::Function, x::CellValue)
CellValue(x1.domain,
    f(x.value))
end

function cellEval(f::Function, x1::CellValue, x2::CellValue)
CellValue(x1.domain,
    f(x1.value, x2.value))
end

function cellEval(f::Function, x1::CellValue, x2::CellValue, x3::CellValue)
CellValue(x1.domain,
    f(x1.value, x2.value, x3.value))
end

function cellEval(f::Function, x1::CellValue, x2::CellValue, x3::CellValue, x4::CellValue)
CellValue(x1.domain,
    f(x1.value, x2.value, x3.value, x4.value))
end

function cellEval(f::Function, x1::CellValue, x2::CellValue, x3::CellValue, x4::CellValue, x5::CellValue)
CellValue(x1.domain,
    f(x1.value, x2.value, x3.value, x4.value, x5.value))
end

function cellEval(f::Function, x1::CellValue, x2::CellValue, x3::CellValue, x4::CellValue, x5::CellValue, x6::CellValue)
CellValue(x1.domain,
    f(x1.value, x2.value, x3.value, x4.value, x5.value, x6.value))
end

# ========================= Generate random perm field ======================
function permfieldlogrndg(Nx,k_avrg,V_dp,cl)
  # 1D random field generator:
  # hdf: Gaussian
  # acf: Gaussian
  Lx=1.0
  x = linspace(-Lx/2,Lx/2,Nx)
  s = -log(1-V_dp)
  mu = log(k_avrg)-s^2/2.0
  Z = s*randn(Nx)

  # Gaussian filter
  F = exp(-x.^2/(cl^2/2.0))

  # correlation of surface using convolution (faltung), inverse
  # Fourier transform and normalizing prefactors
  f = sqrt(2.0/sqrt(pi))*sqrt(Lx/Nx/cl)*ifft(fft(Z).*fft(F))
  perm = exp(mu+real(f))
end

function permfieldlogrnde(Nx,k_avrg,V_dp,cl)
  # 1D random field generator:
  # hdf: Gaussian
  # acf: Exponential
  Lx=1.0
  x = linspace(-Lx/2,Lx/2,Nx)
  s = -log(1-V_dp)
  mu = log(k_avrg)-s^2/2.0
  Z = s*randn(Nx)

  # Gaussian filter
  F = exp(-abs(x)/(cl/2.0))

  # correlation of surface using convolution (faltung), inverse
  # Fourier transform and normalizing prefactors
  f = sqrt(2.0)*sqrt(Lx/Nx/cl)*ifft(fft(Z).*fft(F))
  perm = exp(mu+real(f))
end

# ======================== 2D ==================================>
function permfieldlogrndg(Nx,Ny,k_avrg,V_dp,clx,cly)
# 2D random field generator:
# hdf: Gaussian
# acf: Gaussian
# The surface has a Gaussian height distribution function
# and Gaussian autocovariance functions
  Lx=1.0
  X = linspace(-Lx/2,Lx/2,Nx)
  Y = linspace(-Lx/2,Lx/2,Ny)'

  s = -log(1-V_dp)
  mu = log(k_avrg)-s^2/2
  Z = s*randn(Nx,Ny)
  # Gaussian filter
  F = exp(-(X.^2/(clx^2/2.0).+Y.^2/(cly^2/2.0)))
  # correlated surface generation including convolution (faltning) and inverse
  # Fourier transform and normalizing prefactors
  f = 2.0/sqrt(pi)*Lx/sqrt(Nx*Ny)/sqrt(clx)/sqrt(cly)*ifft(fft(Z).*fft(F))
  perm = exp(mu+real(f))
end

function permfieldlogrnde(Nx,Ny,k_avrg,V_dp,clx,cly)
# 2D random field generator:
# hdf: Gaussian
# acf: Exponential
# The surface has a Gaussian height distribution function
# and Exponential autocovariance functions
  Lx=1.0
  X = linspace(-Lx/2,Lx/2,Nx)
  Y = linspace(-Lx/2,Lx/2,Ny)'

  s = -log(1-V_dp)
  mu = log(k_avrg)-s^2/2
  Z = s*randn(Nx,Ny)
  # Gaussian filter
  F = exp(-(abs(X)/(clx/2).+abs(Y)/(cly/2)))
  # correlated surface generation including convolution (faltning) and inverse
  # Fourier transform and normalizing prefactors
  f = 2.0*Lx/sqrt(Nx*Ny)/sqrt(clx*cly)*ifft(fft(Z).*fft(F))
  perm = exp(mu+real(f))
end
# <========================= 2D ================================

# ======================== 3D ==================================>
function permfieldlogrndg(Nx,Ny,Nz,k_avrg,V_dp,clx,cly,clz)
# 2D random field generator:
# hdf: Gaussian
# acf: Gaussian
# The surface has a Gaussian height distribution function
# and Gaussian autocovariance functions
  Lx=1.0
  X = linspace(-Lx/2,Lx/2,Nx)
  Y = linspace(-Lx/2,Lx/2,Ny)'
  Z = zeros(1,1,Nz)
  Z[:]=linspace(-Lx/2,Lx/2,Nz)
  s = -log(1-V_dp)
  mu = log(k_avrg)-s^2/2
  z = s*randn(Nx,Ny,Nz)
  # Gaussian filter
  F = exp(-(X.^2/(clx^2/2.0).+Y.^2/(cly^2/2.0).+Z.^2/(clz^2/2.0)))
  # correlated surface generation including convolution (faltning) and inverse
  # Fourier transform and normalizing prefactors
  f = 2.0/sqrt(pi)*Lx/(Nx*Ny*Nz)^(1/3)/(clx*cly*clz)^(1/3)*ifft(fft(z).*fft(F))
  perm = exp(mu+real(f))
end

function permfieldlogrnde(Nx,Ny,Nz,k_avrg,V_dp,clx,cly,clz)
# 2D random field generator:
# hdf: Gaussian
# acf: Exponential
# The surface has a Gaussian height distribution function
# and Exponential autocovariance functions
  Lx=1.0
  X = linspace(-Lx/2,Lx/2,Nx)
  Y = linspace(-Lx/2,Lx/2,Ny)'
  Z = zeros(1,1,Nz)
  Z[:]=linspace(-Lx/2,Lx/2,Nz)
  s = -log(1-V_dp)
  mu = log(k_avrg)-s^2/2
  z = s*randn(Nx,Ny,Nz)
  # Gaussian filter
  F = exp(-(abs(X)/(clx/2.0).+abs(Y)/(cly/2.0).+abs(Z)/(clz/2.0)))
  # correlated surface generation including convolution (faltning) and inverse
  # Fourier transform and normalizing prefactors
  f = 2.0*Lx/(Nx*Ny*Nz)^(1/3)/(clx*cly*clz)^(1/3)*ifft(fft(z).*fft(F))
  perm = exp(mu+real(f))
end
# <============================== 3D ==============================

"""
function cellvol = cellVolume(m::MeshStructure)
returns the volume of each cell in the form of a cell variable
"""
function cellVolume(m::MeshStructure)
  dim = m.dimension
  BC = createBC(m)
  if dim==1
          c=m.cellsize.x[2:end-1]
  elseif dim==1.5
          c=2.0*pi()*m.cellsize.x[2:end-1].*m.cellcenters.x
  elseif dim==2
          c=m.cellsize.x[2:end-1]*m.cellsize.y[2:end-1]'
  elseif dim== 2.5 # cylindrical
          c=2.0*pi*m.cellcenters.x.*m.cellsize.x[2:end-1]*m.cellsize.y[2:end-1]'
  elseif dim==2.8 # radial
          c=m.cellcenters.x.*m.cellsize.x[2:end-1]*m.cellsize.y[2:end-1]'
  elseif dim==3
          error("Not available yet")
  elseif dim==3.2
          error("Not available yet")
  end
  cellvol= createCellVariable(m, c, BC)
end

"""
this function reshapes a vetorized cell variable to its domain shape
matrix based on the mesh structure data; it is assumed that the phi
includes the ghost cell data as well.
"""
function reshapeCell(m, phi)
  reshape(full(x), tuple(m.dims+2...))
end

"""
this function reshapes a vetorized cell variable to its domain shape
matrix based on the mesh structure data; it is assumed that the phi
does NOT include the ghost cell data.
"""
function reshapeInternalCell(m, phi)
  reshape(full(x), tuple(m.dims...))
end

"""
returns the internal cells of a cell variable as an array of the same shape
"""
function internalCells(phi::CellValue)
  d = phi.domain.dimension
  N = phi.domain.dims

  if (d==1) || (d==1.5)
  	cellvar= phi.value[2:N[1]+1]
  elseif (d==2) || (d==2.5) || (d==2.8)
  	cellvar= phi.value[2:N[1]+1, 2:N[2]+1]
  elseif (d==3) || (d==3.2)
      cellvar= phi.value[2:N[1]+1, 2:N[2]+1, 2:N[3]+1]
  end
  return cellvar
end

"""
Integrate variable phi over the domain it is defined
"""
function domainInt(phi::CellValue)
  return sum(internalCells(phi).*internalCells(cellVolume(phi.domain)))
end


text:
#__precompile__()
module S3Dicts

using JSON
using AWSCore
using AWSS3
using Retry
using Libz
using Memoize

const NEUROGLANCER_CONFIG_FILENAME = "info"
const AWS_CREDENTIAL = AWSCore.aws_config()
const IS_GZIP = true

const DEFAULT_METADATA = Dict{String, String}(
            "Content-Type"      => "binary/octet-stream",
            "Content-Encoding"  => "gzip")

# map datatype of python to Julia
const DATATYPE_MAP = Dict{String, String}(
    "uint8"     => "UInt8",
    "uint16"    => "UInt16",
    "uint32"    => "UInt32",
    "uint64"    => "UInt64",
    "float32"   => "Float32",
    "float64"   => "Float64"
)

export S3Dict

immutable S3Dict <: Associative
    dir         ::String
    configDict  ::Dict{Symbol,Any}
end

function get_config_dict( h::S3Dict )
    h.configDict
end

"""
split a s3 path to bucket name and key
"""
function splits3(path::AbstractString)
    path = replace(path, "s3://", "")
    bkt, key = split(path, "/", limit = 2)
    return String(bkt), String(key)
end

@memoize function get_config_dict( dir::String )
    configDict=Dict{Symbol,Any}()
    try
        finfo = joinpath( dirname(strip(dir, '/')), NEUROGLANCER_CONFIG_FILENAME )
        @show finfo
        bkt,key = splits3( finfo )
        data = s3_get( AWS_CREDENTIAL, bkt, key )
        if isa(data, Vector{UInt8})
            data = String(data)
        end
        @show data
        configDict = JSON.parse( data, dicttype=Dict{Symbol, Any} )
    catch e
        warn("this is not a neuroglancer formatted dict, did not find the info file: $e")
    end
    @show configDict

    if haskey(configDict, :data_type)
        # a neuroglancer format
        ## postprocessing for neuroglancer format
        configDict[:dataType] = DATATYPE_MAP[ configDict[:data_type] ]

        for d in configDict[:scales]
            if configDict[:num_channels] == 1
                configDict[:chunkSize] = d[:chunk_sizes][1]
            else
                configDict[:chunkSize] = [d[:chunk_sizes][1]..., configDict[:num_channels]]
            end
            configDict[:coding]     = d[:encoding]
            configDict[:size]  = d[:size]
            configDict[:offset]     = d[:voxel_offset]
        end
    end
    return configDict
end

"""
    S3Dict( dir::String )
construct S3Dict from a directory path of s3
"""
function S3Dict( dir::String,
                metadata::Dict{String,String}=DEFAULT_METADATA )
    configDict = get_config_dict(dir)
    configDict[:metadata] = metadata
    S3Dict(dir, configDict)
end

function Base.setindex!(h::S3Dict, v::Array, key::AbstractString)
    @assert ismatch(r"^s3://", h.dir)
    dstFileName = joinpath(h.dir, key)
    bkt, key = splits3(dstFileName)
    @repeat 4 try
        v = Libz.deflate(reinterpret(UInt8, v[:]))
        resp = s3_put(AWS_CREDENTIAL, bkt, key, v,
                      metadata=h.configDict[:metadata])
    catch e
        println("catch error while saving in BigArrays: $e")
        @show typeof(e)
        @delay_retry if true end
    end
end

function Base.getindex(h::S3Dict, key::AbstractString)
    @assert ismatch(r"^s3://", h.dir)
    bucket,key = splits3( joinpath(h.dir, key) )

	@repeat 2 try
        data = AWSS3.s3(AWS_CREDENTIAL, "GET", bucket; path = key, version="")
        return Libz.inflate(data)
    catch e
        @show e
        @show typeof(e)
        @delay_retry if e.code != "NoSuchKey" end
        if e.code == "NoSuchKey"
            throw( KeyError("no such key: $(key)") )
        end
    end
end

function Base.delete!( h::S3Dict, key::AbstractString)
    @assert ismatch(r"^s3://", h.dir)
    fileName = joinpath( h.dir, key )
    bkt, key = splits3( fileName )
    s3_delete(AWS_CREDENTIAL, bkt, key)
end

function Base.keys( h::S3Dict )
    bkt, key = splits3( h.dir )
    keyList = s3_list_objects(AWS_CREDENTIAL, bkt, key)
    return keyList
end

function Base.values(h::S3Dict)
    error("normally values are too large to get them all to RAM")
end

function Base.haskey(h::S3Dict, key::String)
    bkt, dir = splits3( h.dir )
    lst = s3_list_objects(AWS_CREDENTIAL, bkt, joinpath(dir, key))
    return !isempty(lst)
end

end # end of module S3Dicts


text:
function determine_patrol_location!(patrol_location::Matrix{Int64},
                                    patrol_area_time::Matrix{Int64},
                                    drivingtime::Array{Float64,2},
                                    ressource_flow::Array{Int64,3},
                                    locations::Vector{Int64},
                                    districts_grouped::GroupedDataFrame,
                                    mnt::Int64)
    patrol_area_time .+= 1
    for department = 1:size(patrol_location,2)
        for patrolcar = 1:ressource_flow[mnt,department,6]
            if patrol_location[patrolcar, department] == 0
                location = rand(1:nrow(districts_grouped[(locations[department],)]))
                patrol_location[patrolcar, department] = districts_grouped[(locations[department],)][location,:index]
                patrol_area_time[patrolcar, department] = 0
            elseif  patrol_area_time[patrolcar, department] > patrol_area
                location = rand(1:nrow(districts_grouped[(locations[department],)]))
                patrol_location[patrolcar, department] = districts_grouped[(locations[department],)][location,:index]
                patrol_area_time[patrolcar, department] = 0
            end
        end
        for patrolcar = ressource_flow[mnt,department,6]+1:size(patrol_location,1)
            patrol_location[patrolcar, department] = 0
            patrol_area_time[patrolcar, department] = 0
        end
    end
end

function determine_patrol_driving_time(district::Int64,
                                       incidents::DataFrame,
                                       traffic::Array{Float64,2},
                                       patrol_location::Matrix{Int64},
                                       patrol_drive::Vector{Float64},
                                       drivingtime::Array{Float64,2},
                                       current_incident::Int64,
                                       max_drive::Float64,
                                       mnt::Int64)
# calculate the driving time for each patrolling car in the district
    patrol_drive .= patrol_time
    new_time_calculated = false
    for patrol_car = 1:length(patrol_drive)
        if patrol_location[patrol_car,district] > 0
            patrol_drive[patrol_car] = dispatch_drivingtime(drivingtime[patrol_location[patrol_car,district],incidents[current_incident,:location]],traffic,mnt)
            new_time_calculated = true
        end
    end
# if a patrol car is faster than the maximal driving time, select the fastest car on patrol
    fastest_time = patrol_time
    if new_time_calculated == true
        fastest_time = findmin(patrol_drive)[1]
        patrol_car   = findmin(patrol_drive)[2]
        if fastest_time < patrol_time
            patrol_location[patrol_car,district] = 0
        end
    end
    fastest_time = Int(ceil(fastest_time))
    return fastest_time::Int64
end


text:
using SafeTestsets

@safetestset "Fitnessses Tests" begin include("fitnesses_tests.jl") end

@safetestset "Mutations Tests" begin include("mutations_tests.jl") end

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO


text:
# CUDA stream management

immutable CuStream
	handle::Ptr{Void}
	blocking::Bool
	priority::Int
end

function null_stream()
	CuStream(convert(Ptr{Void}, 0), true, 0)
end

function destroy(s::CuStream)
	@cucall(:cuStreamDestroy, (Ptr{Void},), s.handle)
end

function synchronize(s::CuStream)
	@cucall(:cuStreamSynchronize, (Ptr{Void},), s.handle)
end


text:
using BERL
using YAML
using Test
using Cambrian
using CartesianGeneticProgramming
using NeuroEvolution
using UUIDs
cd("..")

"Tests an environment with CGP"
function test_env(env_name::String; env_params...)
    if env_name == "gym"
        name = string("Gym_", env_params[:gym_env])
    elseif env_name == "atari_ram"
        name = string("Atari_", env_params[:atari_game])
    else
        name = env_name
    end
    @testset "$name" begin
        @test env_name in keys(environments)
        cfg_path = "./src/algorithms/CGP/$env_name.yaml"
        cfg = CartesianGeneticProgramming.get_config(cfg_path)
        cfg["n_population"]=10
        cfg["n_gen"]=10
        for (k, v) in pairs(env_params)
            cfg[string(k)] = v
        end

        # Create and setup the environment
        env::BERLenv = environments[env_name](cfg)
        @test env.name != ""

        ind = CGPInd(cfg)
        f = fitness(ind, env)

        if env_name in ["iris"]
            @test 0 <= f[1] <= 1
        end
        @test typeof(f)==Array{Float64, 1}

        new_gen!(env, cfg)
    end
end

"Tests an algo on XOR"
function test_algo(algo_name::String)
    @testset "$algo_name" begin
        cfg_path = "./src/algorithms/$algo_name/xor.yaml"

        if algo_name in ["NEAT", "HyperNEAT", "CPPN"]
            cfg = NeuroEvolution.get_config(cfg_path)
        else
            cfg = CartesianGeneticProgramming.get_config(cfg_path)
        end
        cfg["n_gen"]=10
        cfg["n_population"]=10

        # Create and setup the environment
        env::BERLenv = environments["xor"](cfg)

        # Fitness function
        fit::Function = i::Cambrian.Individual -> fitness(i, env)
        cfg["id"] = "TestAlgo_$(algo_name)_$(env.name)_$(string(UUIDs.uuid4()))"
        evo = algorithms[algo_name](cfg, fit)
        @test typeof(evo)==Cambrian.Evolution
        @test length(evo.population)>0
    end
end

function test_berl()
    # Algorithms
    algs = []
    for (a, v) in pairs(YAML.load_file("./run_config/algorithms.yaml"))
        if a != "runs"
            test_algo(a)
        end
    end

    # Environments
    envs=[]
    for (e, v) in pairs(YAML.load_file("./run_config/environments.yaml"))
        if e in gym_envs
            d = test_env("gym"; gym_env=e)
        elseif e in atari_envs
            d = test_env("atari_ram"; atari_game=e)
        else
            d = test_env(e)
        end
    end

    # Start BERL
    @testset "Start BERL" begin
        m = start_berl("NEAT", "iris")
        @test typeof(m)==Dict{Any,Any}
        @test m["Best fitness"]>0
        @test m["Run time"]>0
        @test m["Memory allocated"]>0
        print_metrics(m)
    end
end

println("Starting tests")

# test_berl()

run_berl()

@testset "plot" begin
    BERLplot("CGP_XOR_f7739035-ba96-4400-b9bf-0492f6ec1da0")
    BERLplot(["CGP_XOR_f7739035-ba96-4400-b9bf-0492f6ec1da0",
    "CGP_XOR_ea455af9-9423-452c-8946-f67bf39a7da4"]; name="test.svg")
end


text:
using StabTrussTopOpt
sto = StabTrussTopOpt
using Makie
using Test

# function example_relaxed_stab_opt()

file_dir = joinpath(pwd(), "test", "instances")
mesh_file = "tim.json"
load_file = "tim_load_support_case.json"

mesh_file = joinpath(file_dir, mesh_file)
load_file = joinpath(file_dir, load_file)

# parse ground mesh, supp, load from json
par = sto.parse_TO_problem(mesh_file, load_file)

# draw ground mesh
init_sc = Makie.Scene()
sto.draw_truss!(init_sc, par.X, par.T, par.Bound)
sto.draw_load!(init_sc, par.X, par.F)
init_sc

# define TO problem
# Steel-S235
E = 2.1e8 # kN/m^2 (210GPa)
τ = 1.0
σ_c = 3.5e5 # kN/m^2 (350MPa)
σ_t = 3.5e5
problem = sto.TO_define(par, E=E, τ=τ, σ_c=σ_c, σ_t=σ_t)

opt_a, opt_q = sto.TO_solve_relaxed_stab(problem)

# a_max = maximum(opt_a)
# for i=1:length(a)
#     if opt_a[i] < 0.001 * a_max
#         opt_a[i] = 0
#     end
# end

# draw solution
final_sc = Scene()
sto.draw_truss!(final_sc, par.X, par.T, par.Bound, area = opt_a, line_width=1e4, stress=opt_q)
# TODO: plot stress color
sto.draw_load!(final_sc, par.X, par.F, load_scale=1e-2)
display(final_sc)
# end

# example_relaxed_stab_opt()