text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

#Symmetric and Hermitian matrices
immutable Symmetric{T,S<:AbstractMatrix} <: AbstractMatrix{T}
    data::S
    uplo::Char
end
"""
    Symmetric(A, uplo=:U)

Construct a `Symmetric` matrix from the upper (if `uplo = :U`) or lower (if `uplo = :L`) triangle of `A`.

# Example

```jldoctest
julia> A = [1 0 2 0 3; 0 4 0 5 0; 6 0 7 0 8; 0 9 0 1 0; 2 0 3 0 4]
5×5 Array{Int64,2}:
 1  0  2  0  3
 0  4  0  5  0
 6  0  7  0  8
 0  9  0  1  0
 2  0  3  0  4

julia> Supper = Symmetric(A)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  2  0  3
 0  4  0  5  0
 2  0  7  0  8
 0  5  0  1  0
 3  0  8  0  4

julia> Slower = Symmetric(A, :L)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  6  0  2
 0  4  0  9  0
 6  0  7  0  3
 0  9  0  1  0
 2  0  3  0  4

julia> eigfact(Supper)
Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([-2.96684,-2.72015,0.440875,7.72015,14.526],[-0.302016 -2.22045e-16 … 1.11022e-16 0.248524; -6.67755e-16 0.596931 … -0.802293 1.93069e-17; … ; 8.88178e-16 -0.802293 … -0.596931 0.0; 0.772108 8.93933e-16 … 0.0 0.630015])
```

`eigfact` will use a method specialized for matrices known to be symmetric.
Note that `Supper` will not be equal to `Slower` unless `A` is itself symmetric (e.g. if `A == A.'`).
"""
Symmetric(A::AbstractMatrix, uplo::Symbol=:U) = (checksquare(A);Symmetric{eltype(A),typeof(A)}(A, char_uplo(uplo)))
immutable Hermitian{T,S<:AbstractMatrix} <: AbstractMatrix{T}
    data::S
    uplo::Char
end
"""
    Hermitian(A, uplo=:U)

Construct a `Hermitian` matrix from the upper (if `uplo = :U`) or lower (if `uplo = :L`) triangle of `A`.

# Example

```jldoctest
julia> A = [1 0 2+2im 0 3-3im; 0 4 0 5 0; 6-6im 0 7 0 8+8im; 0 9 0 1 0; 2+2im 0 3-3im 0 4];

julia> Hupper = Hermitian(A)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  2+2im  0+0im  3-3im
 0+0im  4+0im  0+0im  5+0im  0+0im
 2-2im  0+0im  7+0im  0+0im  8+8im
 0+0im  5+0im  0+0im  1+0im  0+0im
 3+3im  0+0im  8-8im  0+0im  4+0im

julia> Hlower = Hermitian(A, :L)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  6+6im  0+0im  2-2im
 0+0im  4+0im  0+0im  9+0im  0+0im
 6-6im  0+0im  7+0im  0+0im  3+3im
 0+0im  9+0im  0+0im  1+0im  0+0im
 2+2im  0+0im  3-3im  0+0im  4+0im

julia> eigfact(Hupper)
Base.LinAlg.Eigen{Complex{Float64},Float64,Array{Complex{Float64},2},Array{Float64,1}}([-8.32069,-2.72015,3.1496,7.72015,17.1711],Complex{Float64}[-0.231509+0.392692im -4.16334e-17+5.55112e-17im … -2.77556e-17-1.11022e-16im -0.129023-0.00628656im; 0.0+0.0im -0.532524+0.269711im … -0.521035+0.610079im 0.0+0.0im; … ; 0.0-1.38778e-17im 0.715729-0.362501im … -0.387666+0.453917im 0.0-6.93889e-18im; 0.67898-0.0im 0.0+0.0im … 0.0+0.0im -0.661651-0.0im])
```

`eigfact` will use a method specialized for matrices known to be Hermitian.
Note that `Hupper` will not be equal to `Hlower` unless `A` is itself Hermitian (e.g. if `A == A'`).
"""
function Hermitian(A::AbstractMatrix, uplo::Symbol=:U)
    n = checksquare(A)
    for i=1:n
        isreal(A[i, i]) || throw(ArgumentError(
            "Cannot construct Hermitian from matrix with nonreal diagonals"))
    end
    Hermitian{eltype(A),typeof(A)}(A, char_uplo(uplo))
end

typealias HermOrSym{T,S} Union{Hermitian{T,S}, Symmetric{T,S}}
typealias RealHermSymComplexHerm{T<:Real,S} Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}

size(A::HermOrSym, d) = size(A.data, d)
size(A::HermOrSym) = size(A.data)
@inline function getindex(A::Symmetric, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    @inbounds r = (A.uplo == 'U') == (i < j) ? A.data[i, j] : A.data[j, i]
    r
end
@inline function getindex(A::Hermitian, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    @inbounds r = (A.uplo == 'U') == (i < j) ? A.data[i, j] : conj(A.data[j, i])
    r
end

function setindex!(A::Symmetric, v, i::Integer, j::Integer)
    i == j || throw(ArgumentError("Cannot set a non-diagonal index in a symmetric matrix"))
    setindex!(A.data, v, i, j)
end

function setindex!(A::Hermitian, v, i::Integer, j::Integer)
    if i != j
        throw(ArgumentError("Cannot set a non-diagonal index in a Hermitian matrix"))
    elseif !isreal(v)
        throw(ArgumentError("Cannot set a diagonal entry in a Hermitian matrix to a nonreal value"))
    else
        setindex!(A.data, v, i, j)
    end
end

similar{T}(A::Symmetric, ::Type{T}) = Symmetric(similar(A.data, T))
# Hermitian version can be simplified when check for imaginary part of
# diagonal in Hermitian has been removed
function similar{T}(A::Hermitian, ::Type{T})
    B = similar(A.data, T)
    for i = 1:size(A,1)
        B[i,i] = 0
    end
    return Hermitian(B)
end

# Conversion
convert(::Type{Matrix}, A::Symmetric) = copytri!(convert(Matrix, copy(A.data)), A.uplo)
convert(::Type{Matrix}, A::Hermitian) = copytri!(convert(Matrix, copy(A.data)), A.uplo, true)
convert(::Type{Array}, A::Union{Symmetric,Hermitian}) = convert(Matrix, A)
full(A::Union{Symmetric,Hermitian}) = convert(Array, A)
parent(A::HermOrSym) = A.data
convert{T,S<:AbstractMatrix}(::Type{Symmetric{T,S}},A::Symmetric{T,S}) = A
convert{T,S<:AbstractMatrix}(::Type{Symmetric{T,S}},A::Symmetric) = Symmetric{T,S}(convert(S,A.data),A.uplo)
convert{T}(::Type{AbstractMatrix{T}}, A::Symmetric) = Symmetric(convert(AbstractMatrix{T}, A.data), Symbol(A.uplo))
convert{T,S<:AbstractMatrix}(::Type{Hermitian{T,S}},A::Hermitian{T,S}) = A
convert{T,S<:AbstractMatrix}(::Type{Hermitian{T,S}},A::Hermitian) = Hermitian{T,S}(convert(S,A.data),A.uplo)
convert{T}(::Type{AbstractMatrix{T}}, A::Hermitian) = Hermitian(convert(AbstractMatrix{T}, A.data), Symbol(A.uplo))

copy{T,S}(A::Symmetric{T,S}) = (B = copy(A.data); Symmetric{T,typeof(B)}(B,A.uplo))
copy{T,S}(A::Hermitian{T,S}) = (B = copy(A.data); Hermitian{T,typeof(B)}(B,A.uplo))

function copy!(dest::Symmetric, src::Symmetric)
    if src.uplo == dest.uplo
        copy!(dest.data, src.data)
    else
        transpose!(dest.data, src.data)
    end
    return dest
end

function copy!(dest::Hermitian, src::Hermitian)
    if src.uplo == dest.uplo
        copy!(dest.data, src.data)
    else
        ctranspose!(dest.data, src.data)
    end
    return dest
end

ishermitian(A::Hermitian) = true
ishermitian{T<:Real,S}(A::Symmetric{T,S}) = true
ishermitian{T<:Complex,S}(A::Symmetric{T,S}) = all(imag(A.data) .== 0)
issymmetric{T<:Real,S}(A::Hermitian{T,S}) = true
issymmetric{T<:Complex,S}(A::Hermitian{T,S}) = all(imag(A.data) .== 0)
issymmetric(A::Symmetric) = true
transpose(A::Symmetric) = A
ctranspose{T<:Real}(A::Symmetric{T}) = A
function ctranspose(A::Symmetric)
    AC = ctranspose(A.data)
    return Symmetric(AC, ifelse(A.uplo == 'U', :L, :U))
end
function transpose(A::Hermitian)
    AT = transpose(A.data)
    return Hermitian(AT, ifelse(A.uplo == 'U', :L, :U))
end
ctranspose(A::Hermitian) = A
trace(A::Hermitian) = real(trace(A.data))

Base.conj(A::HermOrSym) = typeof(A)(conj(A.data), A.uplo)
Base.conj!(A::HermOrSym) = typeof(A)(conj!(A.data), A.uplo)

#tril/triu
function tril(A::Hermitian, k::Integer=0)
    if A.uplo == 'U' && k <= 0
        return tril!(A.data',k)
    elseif A.uplo == 'U' && k > 0
        return tril!(A.data',-1) + tril!(triu(A.data),k)
    elseif A.uplo == 'L' && k <= 0
        return tril(A.data,k)
    else
        return tril(A.data,-1) + tril!(triu!(A.data'),k)
    end
end

function tril(A::Symmetric, k::Integer=0)
    if A.uplo == 'U' && k <= 0
        return tril!(A.data.',k)
    elseif A.uplo == 'U' && k > 0
        return tril!(A.data.',-1) + tril!(triu(A.data),k)
    elseif A.uplo == 'L' && k <= 0
        return tril(A.data,k)
    else
        return tril(A.data,-1) + tril!(triu!(A.data.'),k)
    end
end

function triu(A::Hermitian, k::Integer=0)
    if A.uplo == 'U' && k >= 0
        return triu(A.data,k)
    elseif A.uplo == 'U' && k < 0
        return triu(A.data,1) + triu!(tril!(A.data'),k)
    elseif A.uplo == 'L' && k >= 0
        return triu!(A.data',k)
    else
        return triu!(A.data',1) + triu!(tril(A.data),k)
    end
end

function triu(A::Symmetric, k::Integer=0)
    if A.uplo == 'U' && k >= 0
        return triu(A.data,k)
    elseif A.uplo == 'U' && k < 0
        return triu(A.data,1) + triu!(tril!(A.data.'),k)
    elseif A.uplo == 'L' && k >= 0
        return triu!(A.data.',k)
    else
        return triu!(A.data.',1) + triu!(tril(A.data),k)
    end
end

-{Tv,S<:AbstractMatrix}(A::Symmetric{Tv,S}) = Symmetric{Tv,S}(-A.data, A.uplo)

## Matvec
A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(y::StridedVector{T}, A::Symmetric{T,S}, x::StridedVector{T}) = BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y)
A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(y::StridedVector{T}, A::Hermitian{T,S}, x::StridedVector{T}) = BLAS.hemv!(A.uplo, one(T), A.data, x, zero(T), y)
##Matmat
A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(C::StridedMatrix{T}, A::Symmetric{T,S}, B::StridedMatrix{T}) = BLAS.symm!('L', A.uplo, one(T), A.data, B, zero(T), C)
A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Symmetric{T,S}) = BLAS.symm!('R', B.uplo, one(T), B.data, A, zero(T), C)
A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(C::StridedMatrix{T}, A::Hermitian{T,S}, B::StridedMatrix{T}) = BLAS.hemm!('L', A.uplo, one(T), A.data, B, zero(T), C)
A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Hermitian{T,S}) = BLAS.hemm!('R', B.uplo, one(T), B.data, A, zero(T), C)

*(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
*(A::StridedMatrix, B::HermOrSym) = A*full(B)

for T in (:Symmetric, :Hermitian), op in (:+, :-, :*, :/)
    # Deal with an ambiguous case
    @eval ($op)(A::$T, x::Bool) = ($T)(($op)(A.data, x), Symbol(A.uplo))
    S = T == :Hermitian ? :Real : :Number
    @eval ($op)(A::$T, x::$S) = ($T)(($op)(A.data, x), Symbol(A.uplo))
end

bkfact(A::HermOrSym) = bkfact(A.data, Symbol(A.uplo), issymmetric(A))
factorize(A::HermOrSym) = bkfact(A)

# Is just RealHermSymComplexHerm, but type alias seems to be broken
det{T<:Real,S}(A::Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}) = real(det(bkfact(A)))
det{T<:Real}(A::Symmetric{T}) = det(bkfact(A))
det(A::Symmetric) = det(bkfact(A))

\{T,S<:StridedMatrix}(A::HermOrSym{T,S}, B::StridedVecOrMat) = \(bkfact(A.data, Symbol(A.uplo), issymmetric(A)), B)

inv{T<:BlasFloat,S<:StridedMatrix}(A::Hermitian{T,S}) = Hermitian{T,S}(inv(bkfact(A.data, Symbol(A.uplo))), A.uplo)
inv{T<:BlasFloat,S<:StridedMatrix}(A::Symmetric{T,S}) = Symmetric{T,S}(inv(bkfact(A.data, Symbol(A.uplo), true)), A.uplo)

eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A)))

eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, irange::UnitRange) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange))

eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, vl::Real, vh::Real) = Eigen(LAPACK.syevr!('V', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)...)
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, vl::Real, vh::Real) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh))

eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A)))

eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, irange::UnitRange) = LAPACK.syevr!('N', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)[1]
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, irange::UnitRange) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange))

eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)[1]
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, vl::Real, vh::Real) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh))

eigmax{T<:Real,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = eigvals(A, size(A, 1):size(A, 1))[1]
eigmin{T<:Real,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = eigvals(A, 1:1)[1]

function eigfact!{T<:BlasReal,S<:StridedMatrix}(A::HermOrSym{T,S}, B::HermOrSym{T,S})
    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : B.data')
    GeneralizedEigen(vals, vecs)
end
function eigfact!{T<:BlasComplex,S<:StridedMatrix}(A::Hermitian{T,S}, B::Hermitian{T,S})
    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : B.data')
    GeneralizedEigen(vals, vecs)
end

eigvals!{T<:BlasReal,S<:StridedMatrix}(A::HermOrSym{T,S}, B::HermOrSym{T,S}) = LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : B.data')[1]
eigvals!{T<:BlasComplex,S<:StridedMatrix}(A::Hermitian{T,S}, B::Hermitian{T,S}) = LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : B.data')[1]

function svdvals!{T<:Real,S}(A::Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}) #  the union is the same as RealHermSymComplexHerm, but right now parametric typealiases are broken
    vals = eigvals!(A)
    for i = 1:length(vals)
        vals[i] = abs(vals[i])
    end
    return sort!(vals, rev = true)
end

#Matrix-valued functions
function expm(A::Symmetric)
    F = eigfact(A)
    return Symmetric((F.vectors * Diagonal(exp.(F.values))) * F.vectors')
end
function expm{T}(A::Hermitian{T})
    n = checksquare(A)
    F = eigfact(A)
    retmat = (F.vectors * Diagonal(exp.(F.values))) * F.vectors'
    if T <: Real
        return real(Hermitian(retmat))
    else
        for i = 1:n
            retmat[i,i] = real(retmat[i,i])
        end
        return Hermitian(retmat)
    end
end

for (funm, func) in ([:logm,:log], [:sqrtm,:sqrt])

    @eval begin

        function ($funm)(A::Symmetric)
            F = eigfact(A)
            if isposdef(F)
                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
            else
                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
            end
            return Symmetric(retmat)
        end

        function ($funm){T}(A::Hermitian{T})
            n = checksquare(A)
            F = eigfact(A)
            if isposdef(F)
                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
                if T <: Real
                    return Hermitian(retmat)
                else
                    for i = 1:n
                        retmat[i,i] = real(retmat[i,i])
                    end
                    return Hermitian(retmat)
                end
            else
                retmat = (F.vectors * Diagonal(($func).(complex(F.values)))) * F.vectors'
                return retmat
            end
        end

    end

end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# Autogenerated wrapper script for xxHash_jll for aarch64-linux-gnu
export libxxhash

JLLWrappers.@generate_wrapper_header("xxHash")
JLLWrappers.@declare_library_product(libxxhash, "libxxhash.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libxxhash,
        "lib/libxxhash.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
function prepcomp(;companame=COMP_A_NAME, compqname = COMP_Q_NAME, compname = COMP_NAME,
   datapath=DATA_PATH, refreshcomp=true,
   incsvstream::Function = IN_CSV_STREAM,
   inbinstream::Function = IN_BIN_STREAM,
   outbinstream::Function = OUT_BIN_STREAM,
   csvextension::String = CSV_EXTENSION,
   binextension::String = BIN_EXTENSION,
   usequarterlydates::Bool=USE_QUARTERLY_DATES,
   parallel::NBool=nothing,
   retainedcolscomp::Vector{Symbol} = RETAINED_COLS_COMP,
   archive::Bool = true)::DataFrame

   (!isnothing(parallel)) && (PARALLEL[] = parallel)

  local compa::DataFrame
  local comp::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcomp
    compa = incsvstream("$datapath\\$companame.$csvextension") |> CSV.File |> DataFrame
    usequarterlydates && (compq = incsvstream("$datapath\\$compqname.$csvextension") |> CSV.File |> DataFrame)

    comp = preprocesscomp(compa)
    comp = mergecompccm!(comp)
    comp = select!(comp,retainedcolscomp)
    outbinstream("$datapath\\$compname.$binextension", comp)

  else
    comp = inbinstream("$datapath\\$compname.$binextension")

  end

  @info "Compustat and CCM linking data loaded and saved into file $compname.$binextension"

  return comp
end

#helper function to clean up names in the quarterly file
function propername(s::Symbol)
  str::String = string(s)

  if str[end] == 'q' || (length(str) ≥ 3 && str[(end-2):end] == "qtr")
    return s
  end

  s = Symbol(s,:q)
  return s
end

#this processes both the q and the a comp frames, with as much commonly used code as possible
function preprocesscomp(comp::DataFrame;
  compdateformat::DateFormat = COMP_DATE_FORMAT,
  testoutput::Bool = TEST_OUTPUT)

  #do operations that are common to both the annual and the quarterly data
  comp = preprocesscomp!(comp) #methods common to compa and compq
  comp = preprocesscompa!(comp) #methods specific to compa

  println("rows: $(size(comp,1))")
  filter!(comp) do r::DataFrameRow
      r[:adate] ≥ eom(r[:fyrmonth])
  end

  println("rows adate≥fdate: $(size(comp,1))")
  validatefadates!(comp; Fgroup=:gvkey, Ffdate=:fdate, Fadate=:adate)
  sort!(comp, [:gvkey, :fdate])

  #form the initial interval candidates
  #these are accurate if the link doesn't changes (adjust for this case later)
  comp = formannouncementintervals!(comp)

  #println(describe(comp))
  if testoutput
    d = Date(2015,12,31)
    testoutput && CSV.write("output\\comp_$d.csv", comp[((bd,ed)->
      (d > bd) && (d ≤ ed)).(comp.begindate, comp.enddate), :])
  end

  return comp
end

#operations that are common to both the annual and the quarterly data
function preprocesscomp!(comp::DataFrame)
  #make all column names lower case
  rename!(comp, (s::String->lowercase(s)).(names(comp)))

  #format dates
  comp.ddate = (s->parsecomp(Date, s)).(comp.datadate)
  comp.dyear = (d::MDate-> ismissing(d) ? missing : Dates.year(d)).(comp.ddate)

  comp = comp[(!ismissing).(comp.dyear), :] #require the calendar date

  rename!(comp, :indfmt=>:indfmtold)
  comp.indfmt = (s::MString->parsecomp(Symbol, s)).(comp.indfmtold)
  select(comp, Not(:indfmtold))

  return comp
end

function preprocesscompa!(compa::DataFrame;
  usedcompafields::Vector{Symbol} = USED_COMPA_VALUE_FIELDS,
  yearrange::UnitRange{Int} = YEAR_RANGE[],
  minannouncementdelay::Int = MIN_ANNOUNCEMENT_DELAY)

  #println(describe(compa))
  #try to use the more efficient year-month format when possible
  compa = compa[(!ismissing).(compa.fyear), :]

  compa.fic = (s::MString->parsecomp(Symbol,s)).(compa.fic)
  compa = compa[compa.fic.==:USA, :]
  compa.fyrmonth = ((y::Integer,m::Integer)->
    YearMonth(y,m)).(compa.fyear, compa.fyr)
  compa.fdate = (eom).(compa.fyrmonth)

  filter!(r::DataFrameRow->year(r.fdate) ∈ yearrange, compa) #restrict to desired years

  #sets the fixed lead of each entry
  @inline function ayrmonth(ym::YearMonth)::YearMonth
    Δ::Int = minannouncementdelay + (ym.y < 2004)
    return ym + Δ#eom(ym + Δ)
  end

  sort!(compa, [:gvkey, :fdate])
  compa.ayrmonth = (ayrmonth).(compa.fyrmonth)
  compa.adate = (eom).(compa.ayrmonth)
  compa = dedupcomp(compa, [:gvkey, :fyear], :fyrmonth)
  #compa = dedupcomp(compa, [:gvkey, :fyear], :fyrmonth) #more of a validation check

  #make the cash field
  compa.cash = ((che::MFloat64, ch::MFloat64)-> coalesce(che, ch, 0.0)).(compa.che, compa.ch)
  compa.cash .= (floororx).(compa.cash)
  compa = compa[(!ismissing).(compa.cash), :] #this is an essential field #should this line go away?
  compa.cashat = (boundunit1).(compa.cash ./ compa.at)
  compa.bknegcash = (boundunit).(1.0 .- compa.cashat)

  compa.debt = (floororx).(compa.dltt .+ compa.dlc)
  compa.netdebt = compa.debt .- compa.cash
  compa.bkequity = (compa.ceq .+
    ((txditc::MFloat64, txdb::MFloat64)-> coalesce(txditc, txdb, 0.0)).(compa.txditc,compa.txdb) #.-
    #(pstk->coalesce(pstk,0.0)).(compa.pstk) #WARNING- commented out since excluded from ceq
    #this is different than how Ivo has it in the original
    )
  compa.bkequityold = (compa.ceq .+
    ((txditc::MFloat64, txdb::MFloat64)-> coalesce(txditc, txdb, 0.0)).(compa.txditc,compa.txdb) .-
    (pstk->coalesce(pstk,0.0)).(compa.pstk)
    )
  #compa.bkequity .= compa.bkequityold #WARNING- only enable for testing purposes
  @mpar for s ∈ [:netdebt, :bkequity, :bkequityold, :at, :lt] #can thread this a bit
    compa[!,s] .= (floororx).(compa[!,s])
  end

  #put in Ivo's minimum values
  compa = filter(r::DataFrameRow->
    (!ismissing(r.at)) &&
    (!ismissing(r.bkequity)) &&
    (r.at ≥ MIN_AT) &&
    (r.bkequity ≥ MIN_BKEQUITY), compa)

  compa.bkflev = (boundunit).(compa.netdebt ./ (compa.netdebt .+ compa.bkequity))

  #compute liab
  compa.netliab = compa.lt .- compa.cash
  compa.bkliab = (boundunit).(compa.netliab ./ compa.at)

  #gross profit is proxy for operating profits as in Novy Marx 2013
  #compa.op = (r::DataFrameRow->ismissing(r.gp) ? r.revt - r.cogs : r.gp).(eachrow(compa))
  compa.nop = (boundabsunit).(compa.pi ./ compa.bkequity)
  #compa.nop = compa.op ./ compa.bkequity #this is used for the operating profit ratio

  #make finite
  @mpar for s ∈ [:bknegcash, :bkflev, :bkliab, :nop]
    compa[!,s] = (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(compa[!,s])
  end

  #get the market value if available
  compa.mkequity = ((mkvalt::MFloat64, pso::MFloat64)->
    coalesce(mkvalt, pso)).(compa.mkvalt, compa.prcc_f .* compa.csho)

  #NOTE: we won't translate these into our measures until after we do the primary merge

  #sort the data
  sort!(compa, [:gvkey, :adate])
  issorted(compa, [:gvkey, :fdate]) || error("compa not sorted on fdate")
  lagwithin2!(compa, [:at], :gvkey, date=:fdate)

  compa.inv = compa.at ./ compa.Lat .- 1.

  #Now make the BRS metrics
  compa.Aat = (r::DataFrameRow->coalesce((r.at + r.Lat)/2, r.at, r.Lat)).(eachrow(compa))
  compa.brse = (compa.sstk .- compa.prstkc .- compa.dv) ./ compa.Aat
  compa.brsd = (r::DataFrameRow->(r.dltis - r.dltr - coalesce(r.dlcch,0.))/r.Aat).(eachrow(compa))
  compa.brsx = compa.brse .+ compa.brsd

  (f::Symbol->finiteormissing!(compa, f)).([:brse, :brsd, :brsx])

  #println(describe(compa))

  return compa
end

  #under no circumstances allow restatements past the subsequent quarter's announcement date
function validatefadates!(comp::DataFrame; Fgroup::NSymbol=nothing, Ffdate::NSymbol=nothing,
  Fadate::NSymbol=nothing, sort::Bool = true)::Nothing

  (isnothing(Fgroup)) && error("Fgroup missing from validatefadates")
  (isnothing(Ffdate)) && error("Ffdate missing from validatefadates")
  (isnothing(Fadate)) && error("Fadate missing from validatefadates")

  sort!(comp, [Fgroup, Fadate])
  T::Type = eltype(comp[!,Ffdate])
  comp.maxfdate_at_adate = Vector{Union{T,Missing}}(undef, size(comp,1))

  compgrp::GroupedDataFrame = groupby(comp, Fgroup)
  @mpar for i ∈ 1:length(compgrp)
    subdf::SubDataFrame = compgrp[i]
    maxfdate::T = subdf[1,Ffdate]
    for r ∈ eachrow(subdf)
      (r[Ffdate] > maxfdate) &&  (maxfdate = r[Ffdate])
      r.maxfdate_at_adate = maxfdate
    end
  end

  filter!(comp) do r::DataFrameRow
    r.maxfdate_at_adate ≤ r[Ffdate]
  end

  select!(comp, Not(:maxfdate_at_adate))

  return nothing
end

function dedupcomp(comp::DataFrame,
  groupbyfields::Vector{Symbol}, datefield::Symbol)

  comp = comp[comp.indfmt .== :INDL,:]

  issorted(comp, [groupbyfields; datefield]) || error("comp not sorted by $datefield")
  #keepflags
  #comp.keep = trues(size(comp,1))

  #this is to match ivo's code

  ncols::Int = size(comp,2)
  scomps::GroupedDataFrame = groupby(comp, groupbyfields)
  @mpar for i ∈ 1:length(scomps)
    scomp::SubDataFrame = scomps[i]
    (size(scomp,1) == 1) && continue #the most likely case

    #keep the most recent version
    maxdate::Union{Int, YearQuarter, Date, YearMonth} = maximum(scomp[!,datefield])
    mindate::Union{Int, YearQuarter, Date, YearMonth} = minimum(scomp[!,datefield])
    #println("got here")

    #Happens with change in fiscal year
    #=if maxdate ≠ mindate
      scomp.keep .= false
    end=#
    #println("got here")
    #scomp.keep .= (yd->isequal(yd,maxdate)).(scomp[!,datefield])

    #special case where we have two records of two different formats (indl and fs)
    #in this case, keep the format with the least missing values
    if (size(scomp,1) > 1) && (maxdate==mindate)
      error("duplicate rows after dropping fs records")
    end

    #if we still have a dup, delete all the records
    #=if sum(scomp.keep) > 1
      @warn "dedup failed group fields:
        $groupbyfields with values $((f->scomp[1,f]).(groupbyfields)) Records removed.\n
        $scomp"
      scomp.keep .= false
    end=#
  end

  #comp = comp[comp.keep, :] #bandaid until filter is fixed

  #select!(comp, Not([:keep, :valuesmissings]))

  return comp
end

#NOTE: messy function
#basically to merge in crsp, the crsp day or month must occur
#in the intersection of the link dates, the period six months following
#the announcement date, and the next announcement date
#this handles everything but the intersection with the link dates
function formannouncementintervals!(comp::DataFrame;
    lastdate::Date = LAST_DATE, months2stale::Int =  MONTHS_2_STALE_LONG)

  sort!(comp, [:gvkey, :adate])

  comp.Nadate = Vector{MDate}(undef, size(comp,1))

  scomps = groupby(comp, :gvkey)
  @mpar for i ∈ 1:length(scomps)
    scomp::SubDataFrame = scomps[i]
    Nsub::Int = size(scomp,1)

    for i ∈ 1:(Nsub-1)
      scomp[i,:Nadate] = scomp[i+1,:adate]

      #check some things that shouldn't happen
      (scomp[i,:adate] ≥ scomp[i,:Nadate]) && error(
        "announcement dates out of order.\n $(scomp[i:(i+1),:])")

      (scomp[i,:fyear] ≥ scomp[i+1,:fyear]) && error(
        "fiscal year dates out of order.\n $(scomp[i:(i+1),:])")
    end

    #final vliad interval extends to the present
    scomp[Nsub, :Nadate] = lastdate
  end

  comp.fdatemax = (eom).(comp.fyrmonth .+ months2stale)

  comp.begindate = comp.adate
  comp.enddate = (min).(comp.Nadate,comp.fdatemax)

  return comp
end

parsecomp(::Type{Date}, s::String;
    compdateformat::DateFormat = COMP_DATE_FORMAT) = Dates.Date(s, compdateformat)

#helper methods and special cases
parsecomp(::Type{<:Any}, v::Missing) = missing
parsecomp(::Type{Date}, i::Int) = parsecomp(Date, "$i")
parsecomp(::Type{Symbol}, s::String) = Symbol(s)

#parse a value where the type is known: Generic case
function parsecomp(::Type{T}, s::String;
    ##NOTE: uncomment when we have the missings parsedmissings = COMP_PARSED_MISSINGS
    ) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  #NOTE: uncomment when we have the missings
  #(!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

function parsecomp(::Type{YearQuarter}, s::String)
  local yq::MYearQuarter

  if length(s) ≠ 6
    yq = missing
  else
    y::MInt = parsecomp(Int, "$(s[1:4])")
    q::MInt = parsecomp(Int, "$(s[6])")
    (q==1) || (q==2) || (q==3) || (q==4) || (q=missing)

    yq = (!ismissing(y)) && (!ismissing(q)) ? YearQuarter(y,q) : missing
  end

  return yq
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Test
using GeoParams

@testset "CreepLaw" begin

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);

# Define a linear viscous creep law ---------------------------------
x1      =   LinearViscous(η=1e18Pa*s)
@test x1.η.val == 1e18Pa*s

x1_ND   =   LinearViscous(η=1e18Pa*s)
@test  isDimensional(x1_ND)==true
Nondimensionalize!(x1_ND,CharUnits_GEO)                 # check that we can nondimensionalize all entries within the struct
@test  isDimensional(x1_ND)==false
@test x1_ND.η*1.0==0.1
Dimensionalize!(x1_ND,CharUnits_GEO)                    # check that we can dimensionalize it again
@test x1_ND.η.val==1e18Pa*s
Nondimensionalize!(x1_ND,CharUnits_GEO)

# perform a computation with the viscous creep laws

# Given stress
@test ComputeCreepLaw_EpsII(1e6Pa, x1, CreepLawVariables())==5e-13/s                # dimensional input
@test ComputeCreepLaw_EpsII(1e0, x1_ND, CreepLawVariables())==5.0                   # non-dimensional
@test ComputeCreepLaw_EpsII([1e0; 2.0], x1_ND, CreepLawVariables())==[5.0; 10.0]    # vector input

# Given strainrate
@test ComputeCreepLaw_TauII(1e-13/s, x1, CreepLawVariables())==1e18*2*1e-13Pa       # dimensional input
@test ComputeCreepLaw_TauII(1e0, x1_ND, CreepLawVariables())==0.2                   # non-dimensional
@test ComputeCreepLaw_TauII([1e0; 2.0], x1_ND, CreepLawVariables())==[0.2; 0.4]     # vector input
# -------------------------------------------------------------------

# -------------------------------------------------------------------
# Define powerlaw viscous rheology
x2  =   PowerlawViscous()
Nondimensionalize!(x2,CharUnits_GEO)
@test x2.ε0.val ==0.001                                 # powerlaw

# perform computations

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

#include("test_DislocationCreep.jl")

end


text:
using DynamicGrids, DynamicGridsGtk, ColorSchemes, Colors

const DEAD = 1
const ALIVE = 2
const BURNING = 3

struct ForestFire{R,N,PC,PR} <: NeighborhoodRule{R}
    neighborhood::N
    prob_combustion::PC
    prob_regrowth::PR
end
ForestFire(; neighborhood=RadialNeighborhood{1}(), prob_combustion=0.0001, prob_regrowth=0.01) =
    ForestFire{DynamicGrids.radius(neighborhood),typeof.((neighborhood, prob_combustion, prob_regrowth))...
              }(neighborhood, prob_combustion, prob_regrowth)

@inline DynamicGrids.applyrule(rule::ForestFire, data, state::Integer, index, hood_buffer) =
    if state == ALIVE
        if BURNING in hood_buffer
            BURNING
        else
            rand() <= rule.prob_combustion ? BURNING : ALIVE
        end
    elseif state in BURNING
        DEAD
    else
        rand() <= rule.prob_regrowth ? ALIVE : DEAD
    end

init = fill(ALIVE, 400, 400)
ruleset = Ruleset(ForestFire(); init=init)

output = ArrayOutput(init, 50)
@time sim!(output, ruleset; tspan=(1, 50))


text:
################################################################################
# InteratomicPotentials API default generic implmentations
################################################################################

potential_energy(A::AbstractSystem, p::ArbitraryPotential) = energy_and_force(A, p).e
force(A::AbstractSystem, p::ArbitraryPotential) = energy_and_force(A, p).f
virial(A::AbstractSystem, p::ArbitraryPotential) = sum(virial_stress(A, p)[1:3])


text:
"""
    build_model()

Helper function to build simple MLP.

# Examples

```julia-repl
using BayesLaplace
nn = build_model()
```

"""
function build_model(;input_dim=2,n_hidden=32,output_dim=1,batch_norm=false,dropout=false,activation=Flux.relu)

    if batch_norm
        nn = Chain(
            Dense(input_dim, n_hidden),
            BatchNorm(n_hidden, activation),
            Dense(n_hidden, output_dim),
            BatchNorm(output_dim)
        )
    elseif dropout
        nn = Chain(
            Dense(input_dim, n_hidden, activation),
            Dropout(0.1),
            Dense(n_hidden, output_dim)
        )
    else
        nn = Chain(
            Dense(input_dim, n_hidden, activation),
            Dense(n_hidden, output_dim)
        )
    end

    return nn

end

using Flux
using Flux.Optimise: update!
"""
    forward_nn(nn, loss, data, opt; n_epochs=200, plotting=nothing)

Wrapper function to train neural network and generate an animation showing the training loss evolution.
"""
function forward_nn(nn, loss, data, opt; n_epochs=200, plotting=nothing)

    avg_l = []

    for epoch = 1:n_epochs
      for d in data
        gs = Flux.gradient(Flux.params(nn)) do
          l = loss(d...)
        end
        update!(opt, Flux.params(nn), gs)
      end
      if !isnothing(plotting)
        plt = plotting[1]
        anim = plotting[2]
        idx = plotting[3]
        avg_loss(data) = mean(map(d -> loss(d[1],d[2]), data))
        avg_l = vcat(avg_l,avg_loss(data))
        if epoch % plotting[4]==0
          plot!(plt, avg_l, color=idx)
          frame(anim, plt)
        end
      end
    end

end

"""
    build_ensemble(K::Int;kw=(input_dim=2,n_hidden=32,output_dim=1))

Helper function to build a simple ensemble composed of `K` MLPs.

# Examples

```julia-repl
using BayesLaplace
𝑬 = build_ensemble(5)
```

"""
function build_ensemble(K=5;kw=(input_dim=2,n_hidden=32,output_dim=1))
    ensemble = [build_model(;kw...) for i in 1:K]
    return ensemble
end

using Statistics
"""
    forward(ensemble, data, opt; loss_type=:logitbinarycrossentropy, plot_loss=true, n_epochs=200, plot_every=20)

Wrapper function to train deep ensemble and generate an animation showing the training loss evolution.
"""
function forward(ensemble, data, opt; loss_type=:logitbinarycrossentropy, plot_loss=true, n_epochs=10, plot_every=1)

    anim = nothing
    if plot_loss
        anim = Animation()
        plt = plot(ylim=(0,1), xlim=(0,n_epochs), legend=false, xlab="Epoch", title="Average (training) loss")
        for i in 1:length(ensemble)
            nn = ensemble[i]
            loss(x, y) = getfield(Flux.Losses,loss_type)(nn(x), y)
            forward_nn(nn, loss, data, opt, n_epochs=n_epochs, plotting=(plt, anim, i, plot_every))
        end
    else
        plt = nothing
        for nn in ensemble
            loss(x, y) = getfield(Flux.Losses,loss_type)(nn(x), y)
            forward_nn(nn, loss, data, opt, n_epochs=n_epochs, plotting=plt)
        end
    end

    return ensemble, anim
end;

using BSON: @save
"""
    save_ensemble(ensemble::AbstractArray; root="")

Saves all models in ensemble to disk.
"""
function save_ensemble(ensemble::AbstractArray; root="")
    if !isdir(root)
        mkdir(root)
    end
    for i in 1:length(ensemble)
        path = root * "/nn" * string(i) * ".bson"
        model = ensemble[i]
        @save path model
    end
end

using BSON: @load
"""
    load_ensemble(root="")

Loads all models in `root` folder and stores them in a list.
"""
function load_ensemble(;root="")
    all_files = Base.Filesystem.readdir(root)
    is_bson_file = map(file -> Base.Filesystem.splitext(file)[2][2:end], all_files) .== "bson"
    bson_files = all_files[is_bson_file]
    bson_files = map(file -> root * "/" * file, bson_files)
    ensemble = []
    for file in bson_files
        @load file model
        ensemble = vcat(ensemble, model)
    end
    return ensemble
end

# Plot data points:
using Plots
"""
    plot_data!(plt,X,y)

# Examples

```julia-repl
using CounterfactualExplanations
using CounterfactualExplanations.Data: toy_data_linear
using CounterfactualExplanations.Utils: plot_data!
X, y = toy_data_linear(100)
plt = plot()
plot_data!(plt, hcat(X...)', y)
```

"""
function plot_data!(plt,X,y)
    Plots.scatter!(plt, X[:,1],X[:,2],group=Int.(y),color=Int.(y))
end

# Plot contour of posterior predictive:
using Plots, CounterfactualExplanations.Models
"""
    plot_contour(X,y,M;colorbar=true,title="",length_out=50,type=:laplace,zoom=0,xlim=nothing,ylim=nothing)

Generates a contour plot for the posterior predictive surface.

# Examples

```julia-repl
using CounterfactualExplanations
using CounterfactualExplanations.Data: toy_data_linear
using CounterfactualExplanations.Utils: plot_contour
X, y = toy_data_linear(100)
X = hcat(X...)'
β = [1,1]
M =(β=β,)
predict(M, X) = σ.(M.β' * X)
plot_contour(X, y, M)
```

"""
function plot_contour(X,y,M;colorbar=true,title="",length_out=50,zoom=-1,xlim=nothing,ylim=nothing,linewidth=0.1)

    # Surface range:
    if isnothing(xlim)
        xlim = (minimum(X[:,1]),maximum(X[:,1])).+(zoom,-zoom)
    else
        xlim = xlim .+ (zoom,-zoom)
    end
    if isnothing(ylim)
        ylim = (minimum(X[:,2]),maximum(X[:,2])).+(zoom,-zoom)
    else
        ylim = ylim .+ (zoom,-zoom)
    end
    x_range = collect(range(xlim[1],stop=xlim[2],length=length_out))
    y_range = collect(range(ylim[1],stop=ylim[2],length=length_out))
    Z = [Models.probs(M,[x, y])[1] for x=x_range, y=y_range]

    # Plot:
    plt = contourf(
        x_range, y_range, Z';
        colorbar=colorbar, title=title, linewidth=linewidth,
        xlim=xlim,
        ylim=ylim
    )
    plot_data!(plt,X,y)

end

# Plot contour of posterior predictive:
using Plots, CounterfactualExplanations.Models
"""
    plot_contour_multi(X,y,M;colorbar=true,title="",length_out=50,type=:laplace,zoom=0,xlim=nothing,ylim=nothing)

Generates a contour plot for the posterior predictive surface.

# Examples

```julia-repl
using BayesLaplace, Plots
import BayesLaplace: predict
using Flux: σ
X, y = toy_data_linear(100)
X = hcat(X...)'
β = [1,1]
M =(β=β,)
predict(M, X) = σ.(M.β' * X)
plot_contour(X, y, M)
```

"""
function plot_contour_multi(X,y,M;title="",length_out=50,zoom=-1,xlim=nothing,ylim=nothing,linewidth=0.1)

    # Surface range:
    if isnothing(xlim)
        xlim = (minimum(X[:,1]),maximum(X[:,1])).+(zoom,-zoom)
    else
        xlim = xlim .+ (zoom,-zoom)
    end
    if isnothing(ylim)
        ylim = (minimum(X[:,2]),maximum(X[:,2])).+(zoom,-zoom)
    else
        ylim = ylim .+ (zoom,-zoom)
    end
    x_range = collect(range(xlim[1],stop=xlim[2],length=length_out))
    y_range = collect(range(ylim[1],stop=ylim[2],length=length_out))
    Z = reduce(hcat, [Models.probs(M,[x, y]) for x=x_range, y=y_range])

    # Plot:
    plt = plot()
    plot_data!(plt,X,y)
    out_dim = size(Z)[1]
    for d in 1:out_dim
        contour!(
            plt,
            x_range, y_range, Z[d,:];
            colorbar=false, title=title,
            xlim=xlim,
            ylim=ylim,
            colour=d
        )
    end

    return plot(plt)

end


text:
module TestFFIDatasets

using LightGBM
using Test
using Random

# we don't want the LightGBM vom
redirect_stderr()

verbosity = "verbose=-1"

@testset "LGBM_DatasetCreateFromMat -- floats/column major" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    @test created_dataset.handle != C_NULL
    @test LightGBM.LGBM_DatasetGetNumData(created_dataset) == 3
    @test LightGBM.LGBM_DatasetGetNumFeature(created_dataset) == 2

    @test LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b"]) == nothing
    # cause its not row major and has 2 cols
    @test_throws ErrorException LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b", "c"])

end

@testset "LGBM_DatasetCreateFromMat -- floats/row major" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity, true)

    @test created_dataset.handle != C_NULL
    @test LightGBM.LGBM_DatasetGetNumData(created_dataset) == 2
    @test LightGBM.LGBM_DatasetGetNumFeature(created_dataset) == 3

    # this kind of stuff is a bit circular but without a better way to interact with a dataset, its hard to test
    @test_throws ErrorException LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b"])
    # cause its row major with 2 initial cols and 3 rows
    @test LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b", "c"]) == nothing

end

@testset "LGBM_DatasetCreateFromMat -- ints/column major" begin

    mymat = [1 2; 3 4; 5 6]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    @test created_dataset.handle != C_NULL
    @test LightGBM.LGBM_DatasetGetNumData(created_dataset) == 3
    @test LightGBM.LGBM_DatasetGetNumFeature(created_dataset) == 2

    @test LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b"]) == nothing
    # cause its not row major and has 2 cols
    @test_throws ErrorException LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b", "c"])

end

@testset "LGBM_DatasetCreateFromMat -- ints/row major" begin

    mymat = [1 2; 3 4; 5 6]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity, true)

    @test created_dataset.handle != C_NULL
    @test LightGBM.LGBM_DatasetGetNumData(created_dataset) == 2
    @test LightGBM.LGBM_DatasetGetNumFeature(created_dataset) == 3

    @test_throws ErrorException LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b"])
    # cause its row major with 2 initial cols and 3 rows
    @test LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, ["a", "b", "c"]) == nothing

end

@testset "LGBM_DatasetGetSubset" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    get_idx = [1, 3] # check julia indices are correctly translated
    # gotta create the same so its not a reference, so we can check the mutation is correct (or at least it looks like it hasn't been mutated)
    get_idx_copy = [1, 3]
    faulty_idx = [1, 5]
    # this isn't a julia matrix -- its another dataset
    fetched = LightGBM.LGBM_DatasetGetSubset(created_dataset, Int32.(get_idx), verbosity)
    fetched2 = LightGBM.LGBM_DatasetGetSubset(created_dataset, get_idx, verbosity) # check int64s

    # check it isn't the same dataset...
    @test fetched.handle != created_dataset.handle
    @test get_idx == get_idx_copy

    @test_throws ErrorException LightGBM.LGBM_DatasetGetSubset(created_dataset, faulty_idx, verbosity)

end

@testset "LGBM_Dataset<Get|Set>FeatureNames" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    fieldnames = [randstring(22), randstring(50)]

    @test LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, fieldnames) == nothing
    @test LightGBM.LGBM_DatasetGetFeatureNames(created_dataset) == fieldnames

    # test string overlengths
    name_1 = randstring(3000)
    name_2 = randstring(3000)

    @test LightGBM.LGBM_DatasetSetFeatureNames(created_dataset, [name_1, name_2]) == nothing
    retrieved = LightGBM.LGBM_DatasetGetFeatureNames(created_dataset)

    @test retrieved[1] == name_1#[1:256]
    @test retrieved[2] == name_2#[1:256]

end

@testset "LGBM_DatasetFree" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    # These tests exposed a double-free
    @test LightGBM.LGBM_DatasetFree(created_dataset) == nothing
    @test created_dataset.handle == C_NULL

end

@testset "LGBM_DatasetSetField" begin
    # The current implementation for LGBM_DatasetSetField has questionable type promotions
    # but it touches too much code in order to tidy up right now
    @test_broken false

end

@testset "LGBM_DatasetGetField" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    created_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    # pre-emptively use the right types in anticipation of typing for SetField being fixed in the future
    weights = [0.25f0, 0.5f0, 0.75f0] # Float32
    labels = [0.2f0, 0.5f0, -1.2f0] # Float32
    init_scores = [2.2, 0.1, 0.7] # Float64
    groups = Int32.([2, 1]) # Int32

    LightGBM.LGBM_DatasetSetField(created_dataset, "weight", weights)
    LightGBM.LGBM_DatasetSetField(created_dataset, "label", labels)
    LightGBM.LGBM_DatasetSetField(created_dataset, "group", groups)
    LightGBM.LGBM_DatasetSetField(created_dataset, "init_score", init_scores)

    @test LightGBM.LGBM_DatasetGetField(created_dataset, "weight") == weights
    @test isa(LightGBM.LGBM_DatasetGetField(created_dataset, "weight"), Vector{Float32})

    @test LightGBM.LGBM_DatasetGetField(created_dataset, "label") == labels
    @test isa(LightGBM.LGBM_DatasetGetField(created_dataset, "label"), Vector{Float32})

    @test LightGBM.LGBM_DatasetGetField(created_dataset, "init_score") == init_scores
    @test isa(LightGBM.LGBM_DatasetGetField(created_dataset, "init_score"), Vector{Float64})

    # I dont understand the discrepancy between input and output
    # It looks like LightGBM returns the cumsum of the group query
    # counts with 0 prepended but that isn't mentioned in docs
    @test_broken LightGBM.LGBM_DatasetGetField(created_dataset, "group") == groups
    @test isa(LightGBM.LGBM_DatasetGetField(created_dataset, "group"), Vector{Int32})

end

end # module


text:
include(joinpath(@__DIR__, "pedagogy_spatial.jl"))
using Random

# Experimental setup refers to slide 9 of [E. Rubak](https://www-ljk.imag.fr/membres/Jean-Francois.Coeurjolly/documents/lecture4.pdf)
β, γ, r = 2.0, 0.2, 0.7  # γ = 0 ≡ Hard core

c, w = [0.0, 0.0], 10.0
win = PRS.SquareWindow(c, w)

strauss = PRS.StraussPointProcess(β, γ, r, win)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_dcftp(rng, strauss)

p = plot(strauss, sample; title="")
path = joinpath(@__DIR__, "output", "strauss", "strauss_dcftp.pdf")
Plots.savefig(p, path)


text:
using Pkg; Pkg.activate(@__DIR__)
using Comrade
using Distributions
using Pathfinder
using AdvancedHMC
using Plots
using Zygote
using Parameters
using GalacticOptim
using Optim
using NLopt
using BlackBoxOptim
using Dynesty
using AdaptiveMCMC

# load eht-imaging we use this to load eht data
load_ehtim()
# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about
# large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true).add_fractional_noise(0.01)
# extract log closure amplitudes and closure phases
dlcamp = extract_lcamp(obs)
dcphase = extract_cphase(obs)
lklhd = RadioLikelihood(dlcamp, dcphase)

# build the model here we fit a ring with a azimuthal
#brightness variation and a Gaussian

struct Model{C,F}
    cache::C
    fov::F
    npix::Int
end

function (model::Model)(θ)
    (;c, f, r, σ, τ, ξ, ma, mp, x, y) = θ
    #cmat = f*reshape(c, model.npix, model.npix)
    #img = IntensityMap(cmat, model.fov, model.fov, BSplinePulse{3}())
    #mimg = modelimage(img, model.cache)
    s,c = sincos(mp)
    α = ma*c
    β = ma*s
    #ring = (1-f)*rotated(smoothed(stretched(MRing((α,), (β,)), r, r*τ),σ), ξ)
    return MRing((α,), (β,))#stretched(MRing((α,), (β,)), r, r)#shifted(ring, x, y) #+ mimg
end

fovxy = μas2rad(90.0)
npix = 7
prior = (
          c = Dirichlet(npix*npix, 0.5),
          f = Uniform(0.0, 1.0),
          r = Uniform(μas2rad(10.0), μas2rad(30.0)),
          σ = truncated(Normal(0.0, μas2rad(0.5)), μas2rad(0.01), μas2rad(20.0)),
          τ = Uniform(0.5, 1.0),
          ξ = Uniform(-π/2, π/2),
          ma = Uniform(0.0, 0.5),
          mp = Uniform(-π, π),
          x = Uniform(-fovxy/2, fovxy/2),
          y = Uniform(-fovxy/2, fovxy/2)
        )
# Now form the posterior
cache = create_cache(Comrade.DFTAlg(dlcamp), IntensityMap(rand(npix,npix), fovxy, fovxy, BSplinePulse{3}()))
mms = Model(cache, fovxy, npix)

post = Posterior(lklhd, prior, mms)
tpost = asflat(post)
# We will use HMC to sample the posterior.
# First to reduce burn in we use pathfinder
ndim = dimension(tpost)
f, tr = GalacticOptim.OptimizationFunction(post, GalacticOptim.AutoForwardDiff{25}())
x0 = inverse(tr, rand(post.prior))
prob = GalacticOptim.OptimizationProblem(f, x0, nothing, lb=fill(-5.0, dimension(tr)), ub=fill(5.0, dimension(tr)))
sol, xopt = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(), tr; maxiters=1000_000)

prob = GalacticOptim.OptimizationProblem(f, 3*randn(ndim), nothing)
sol, xopt = solve(prob, LBFGS(), tr; show_trace=true, iterations=4000, g_tol=1e-2)

start = [transform(tr, 3*randn(ndim)) for _ in 1:50]
q, ϕ, compinds = multipathfinder(post, 50; init_params=start, importance=true, iterations=2000)

residual(mms(xopt), dlcamp)
residual(mms(xopt), dcphase)
img = intensitymap(mms(xopt), μas2rad(160.0), μas2rad(160.0), 512, 512)
plot(abs.(img) , xlims=(-60.0,60.0), ylims=(-60.0,60.0), colorbar_scale=:log10, clims=(maximum(img)/100, maximum(img)),)

# now we sample using hmc
metric = DenseEuclideanMetric(ndim)
chain, stats = sample(post, HMC(;metric, autodiff=Comrade.AD.ForwardDiffBackend(;chunksize=Val(25))), 10000; nadapts=7000, init_params=xopt)

serialize("m87hybrid_fits.jls", Dict(:chain=>chain, :stats=>stats, :fov=>fovxy, :npix=>npix, :prior=>priors))


text:
using BenchmarkTools
using DataFrames
using CSV
using Statistics
using StatsBase: sample, quantile
using EvoTrees
using EvoTrees: sigmoid, logit
using Plots
using GraphRecipes

# prepare a dataset
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# linear
params1 = EvoTreeRegressor(
    loss=:linear, metric=:mse,
    nrounds=100, nbins = 100,
    λ = 0.5, γ=0.1, η=0.05,
    max_depth = 4, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)

# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval, print_every_n = 25, metric=:mae)
@time pred_train_linear = predict(model, X_train)
# @time pred_eval_linear = predict(model, X_eval)
# mean(abs.(pred_train_linear .- Y_train))
# sqrt(mean((pred_train_linear .- Y_train) .^ 2))

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")

params1 = EvoTreeRegressor(
    loss=:linear, metric=:mse,
    nrounds=100, nbins = 100,
    λ = 0.0, γ=0.0, η=0.1,
    max_depth = 2, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

anim = @animate for i=1:20
    params1.nrounds = (i-1)*5+1
    model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = Inf)
    pred_train_linear = predict(model, X_train)
    x_perm = sortperm(X_train[:,1])
    plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
    plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")
end
gif(anim, "anim_fps1.gif", fps = 1)

# plot tree
function treemat(tree)
    mat = zeros(Int, length(tree), length(tree))
    for i in 1:length(tree)
        if tree[i].split
            mat[i,tree[i].left] = 1
            mat[tree[i].left, i] = 1
            mat[i,tree[i].right] = 1
            mat[tree[i].right, i] = 1
        end
        mat = sparse(mat)
    end
    return mat
end

# plot tree
function nodenames(tree)
    names = []
    for i in 1:length(tree)
        if tree[i].split
            push!(names, "feat: " * string(tree[i].feat) * "\n< " * string(round(tree[i].cond, sigdigits=3)))
        else
            push!(names, "pred: " * string(round(tree[i].pred[1], sigdigits=3)))
        end
    end
    return names
end

tree1 = model.trees[2].nodes
mat = treemat(tree1)
nodes = nodenames(tree1)
p = graphplot(mat, method=:tree, node_weights=ones(length(tree1)) .* 10, names = nodes, linecolor=:grey, nodeshape=:ellipse, nodecolor="#66ffcc")
savefig(p, "tree1.svg")


text:
using Pkg
Pkg.activate(".")

verbose = true
 make_plots = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments, EchelleInstruments.NEID
 using EchelleCCFs
 #=
 using RvSpectML
 using NeidSolarScripts
 using NeidSolarScripts.SolarRotation
 =#
 if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, StatsBase, Statistics, Dates
 using JLD2, FileIO
 using NaNMath

paths_to_search_for_param = [pwd(),"examples",joinpath(pkgdir(RvSpectMLBase),"..","RvSpectML","examples"), "/gpfs/scratch/jpn23/"]
  eval(read_data_paths(paths_to_search=paths_to_search_for_param))
  #target_subdir = "verygood_days"   # USER: Replace with directory of your choice
  target_subdir = "good_days"   # USER: Replace with directory of your choice
  fits_target_str = "Sun"
  output_dir = "output/continuum"
  #outputs = Dict{String,Any}()

if verbose println("# Finding what data files are avaliable.")  end
manifest_filename = joinpath(neid_data_path,target_subdir,"manifest.csv")
#if isfile("manifest.csv")
if isfile(manifest_filename)
    #df_files  = CSV.read("manifest.csv", DataFrame)
    df_files  = CSV.read(manifest_filename, DataFrame)
    @assert size(df_files,1) >= 1
    @assert hasproperty(df_files,:Filename)
    @assert hasproperty(df_files,:target)
    @assert hasproperty(df_files,:bjd)
    @assert hasproperty(df_files,:ssbz)
    @assert hasproperty(df_files,:exptime)
else
    @assert isdefined(Main,:neid_data_path)
    df_files = make_manifest(neid_data_path, target_subdir, NEID )
    #CSV.write("manifest.csv", df_files)
    CSV.write(manifest_filename, df_files)
end

if verbose println("# Reading in customized parameters from param.jl.")  end
   if !@isdefined(idx_day_to_use)
      if haskey(ENV,"PBS_ARRAYID")
         idx_day_to_use = parse(Int64,ENV["PBS_ARRAYID"])
         println("# Found PBS_ARRAYID.  Processing day ", idx_day_to_use, ".")
      elseif haskey(ENV,"PBS_JOBNAME")
         m = match(r"day=(\d+)", ENV["PBS_JOBNAME"])
         if length(m.captures) >= 1
            idx_day_to_use = parse(Int64,m.captures[1])
            println("# Found PBS_JOBNAME, processing day ", idx_day_to_use, ".")
         else
            idx_day_to_use = 1
            println("# No PBS_ARRAYID or day in PBS_JOBNAME provided, processing first day.")
         end
      else
         println("# No PBS_ARRAYID or PBS_JOBNAME provided, processing first day.")
         idx_day_to_use = 1
      end
   else
      println("# idx_day_to_use = ", idx_day_to_use)
   end
   eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))

   #=
   outputs_filename = joinpath(output_dir,"solar_" * date_str * "_new.jld2")
   if isfile(outputs_filename) && false
     times_already_processed = load(outputs_filename, "times")
     files_in_day_to_process = size(df_files_solar_by_day.data[idx_day_to_use],1)
      if files_in_day_to_process == length(times_already_processed)
         println("# Already processed all ", length(times_already_processed), " files for ", date_str)
         exit()
      end
   end
   =#

using Distributed
if haskey(ENV,"PBS_NUM_PPN")
  num_procs_to_add = parse(Int64,ENV["PBS_NUM_PPN"])
  addprocs(num_procs_to_add)
else
  addprocs()
end

@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using RvSpectMLBase
@everywhere using EchelleInstruments

@everywhere using RCall
@everywhere afs_src = joinpath(pwd(),"src","AFS.R")
@everywhere R"source($afs_src)"

@everywhere function calc_continuum_model(spectrum::AbstractSpectra2D; order_idx::Integer )
    possible_pix = get_pixel_range(get_inst(spectrum),order_idx)
    bad_pix = bad_col_ranges(get_inst(spectrum),order_idx)
    pix_rng = EchelleInstruments.calc_complement_index_ranges(possible_pix,bad_pix)
    pix = mapreduce(p->collect(p),vcat,pix_rng)
    afs_inputs = zeros(Float64,length(pix),2)
    afs_inputs[:,1] .= spectrum.λ[pix,order_idx]
    afs_inputs[:,2] .= spectrum.flux[pix,order_idx]
    @assert !any(isnan.(afs_inputs))
    #=
    wv = mapreduce(p->spec.λ[p,order_idx],vcat,pix_rng)
    @assert !any(isnan.(wv))
    inten = mapreduce(p->convert(Vector{Float64},spec.flux[p,order_idx]),vcat,pix_rng)
    @assert !any(isnan.(inten))
    afs_inputs = hcat(wv,inten)
    =#
    #df = DataFrame("wv"=>wv,"intes"=>inten)
    afs_output_R = R"AFS($afs_inputs,0.95,0.25)"
    afs_output = rcopy(afs_output_R)
    continuum = zeros(eltype(spectrum.flux),size(spectrum.flux[:,order_idx]))
    continuum = fill(NaN,size(spectrum.flux[:,order_idx]))
    continuum[pix] .= afs_output
    return continuum
end

@everywhere using EchelleInstruments.NEID
@everywhere function calc_continuum_model(spectrum::AbstractSpectra2D )
    vec_of_orders = pmap(ord->calc_continuum_model(spectrum,order_idx=ord), min_order(get_inst(spectrum)):max_order(get_inst(spectrum)) )
    output = fill(NaN, size(spectrum.flux))
    for (i,ord) in enumerate(min_order(get_inst(spectrum)):max_order(get_inst(spectrum)))
        output[:,ord] .= vec_of_orders[i]
    end
    return output
end

num_days_to_process = size(df_files_solar_by_day,1)
#for idx_day_to_use in 1:num_days_to_process

  df_files_use = df_files_solar_by_day[idx_day_to_use,:data] |> @orderby(_.bjd) |> DataFrame

  println("# *** Working on day ", idx_day_to_use, " with ", size(df_files_use,1), "files. ***" )
  @assert 1 <= idx_day_to_use <= size(df_files_solar_by_day,1)
  #@distributed
  for row in eachrow(df_files_use)
    spec = NEID.read_data(row)
    m = match(r"neidL1_(\d+)[T_](\d+)\.fits$",row.Filename)
    output_filename = "neidL1_" * m.captures[1] * "T" * m.captures[2] * "_continuum=afs.jld2"
    output_filename = joinpath(neid_data_path,target_subdir,output_dir,output_filename)
    if isfile(output_filename)
       println("# Skipping ", output_filename)
       continue
    end
    println("# Working on ", output_filename)
    continuum = calc_continuum_model(spec)
    @save output_filename continuum
  end

#end


text:
#special type which represents an LUMatrix with permutations
struct LUMatrix{T}
    m::Matrix{T}
    p::Vector{Int64}
end

#another testing function that makes life much easier:
luupper(lu::LUMatrix{T}) where {T} = lu.m |> UpperTriangular
lulower(lu::LUMatrix{T}) where {T} = (lu.m - luupper(lu) + eye(T, size(lu.m,1))) |> LowerTriangular

#just a testing function that tests the LU matrix.
function lureassemble(lu::LUMatrix{T}) where {T}
    lulower(lu) * luupper(lu)
end

# TODO: implement unit testing on this function.

@generated function solve_quire(A::LowerTriangular{T}, B::AbstractVector{T}, quire = Quire(T), ::Type{Val{onediag}} = Val{false}) where {T, onediag}
    first_factor, row_factor = if onediag
        :(B[1]), :(tmp)
    else
        :(B[1] / A[1,1]), :(tmp / A[row, row])
    end

    quote
        m,n = size(A); l = length(B)
        (l == m) || throw(DimensionMismatch("matrix is ($m x $n) and vector is $l"))
        res = Vector{T}(n)
        tmp = zero(T)

        #start with the first row, which does not require a quire.
        @inbounds res[1] = $first_factor
        #TODO: reimplement second row with a proper fused-multiply subtract instead of
        #quire.
        for row = 2:n
            set!(quire, B[row])
            for col = 1:row
                @inbounds tmp =fdp!(quire, res[col], -A[row, col])
            end
            @inbounds res[row] = $row_factor
        end
        res
    end
end

function solve_quire(A::UpperTriangular{T}, B::AbstractVector{T}, quire = Quire(T)) where {T}
    m,n = size(A); l = length(B)
    (l == m) || throw(DimensionMismatch("matrix is ($m x $n) and vector is $l"))
    res = Vector{T}(n)
    tmp = zero(T)

    #start with the last row, which does not require a quire.
    @inbounds res[m] = B[m] / A[m,m]
    #TODO: reimplement penultimate row with a proper fused-multiply subtract instead of
    #quire.
    for row = (m-1):-1:1
        @inbounds set!(quire, B[row])
        for col = (row+1):m
            @inbounds tmp =fdp!(quire, res[col], -A[row, col])
        end
        @inbounds res[row] = tmp / A[row, row]
    end
    res
end

function solve_quire(LM::LUMatrix{T}, v::AbstractVector{T}, quire = Quire(T)) where {T}
    M = LM.m
    p = LM.p
    #first, solve the lower matrix against the permuted values in the vector.

    m,n = size(M); l = length(v)
    (l == m) || throw(DimensionMismatch("matrix is ($m x $n) and vector is $l"))
    res = Vector{T}(m)
    tmp = zero(T)

    #next, permute the incoming vector
    for idx = 1:m
        res[idx] = v[p[idx]]
    end

    res = solve_quire(LowerTriangular(M), res, quire, Val{true})
    #CONTINUE WITH CODE FOR THE UPPER MATRIX.

    solve_quire(UpperTriangular(M), res, quire)
end

function lu_factors!(M::AbstractMatrix{T}, quire = Quire(T))::LUMatrix{T} where {T}
    # lufact_quire is only going to work on square matrices at this time.
    # this returns an LU factorized A along with a pivot permutation.
    m, n = size(M)
    info = 0
    #set up the pivot list.
    pivot_list = collect(1:m)
    @inbounds begin
        #here we use the doolittle algorithm for LU factorization because it
        #permits the use of the quire for numerical accuracy.  The doolittle
        #algorithm operates on one row of the upper matrix, followed by the
        #corresponding column of the lower matrix, so activity should be index
        #ed by diagonal.
        for diag = 1:m
            # here we're going to do pivots to optimize the process.
            # search for the maximum absolute value in this column.
            pidx = indmax(abs.(M[diag:m, diag])) + diag - 1
            if pidx != diag
                #swap these rows
                (pivot_list[diag], pivot_list[pidx]) = (pivot_list[pidx], pivot_list[diag])

                #do the row interchange.
                if M[pidx, diag] != zero(T)
                    for col = 1:m
                        (M[diag, col], M[pidx, col]) = (M[pidx, col], M[diag, col])
                    end
                else
                    pivot_list[diag] = 0
                end
            end

            # solve columns in the upper matrix
            for upper_col = diag:m
                #signal which dot product we're gonna do.
                #then set up a
                accum = M[diag, upper_col]
                set!(quire, accum)
                for lower_col = 1:(diag-1)
                    accum = fdp!(
                      quire,
                      M[diag, lower_col],
                      -M[lower_col, upper_col])
                end
                M[diag, upper_col] = accum
            end

            #solve rows in the lower matrix
            factor = inv(M[diag, diag])
            for lower_row = (diag + 1):n
                accum = M[lower_row, diag]
                set!(quire, accum)
                for upper_row = 1:(diag-1)
                    accum = fdp!(
                      quire,
                      M[lower_row, upper_row],
                      -M[upper_row, diag])
                end
                M[lower_row, diag] = accum * factor
            end
        end
    end
    LUMatrix(M, pivot_list)
end

function solve_quire(A::AbstractMatrix{T}, B::AbstractVector{T}, quire = Quire(T)) where {T}
    m, n = size(A)
    if m == n
        if istril(A)
            if istriu(A)
                return Diagonal(A) \ B
            else
                return solve_quire(LowerTriangular(A), B, quire)
            end
        end
        if istriu(A)
            return solve_quire(UpperTriangular(A), B, quire)
        end
        return solve_quire(lu_factors!(copy(A)), B, quire)
    end
    throw(ArgumentError("currently, nonsquare matrices are not supported"))
end

function solve_quire_refine(A::AbstractMatrix{T}, B::AbstractVector{T}) where {T}
    res = A \ B
    res = refine(res, A, B)
    res = refine(res, A, B)
end

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
# ==============================================================================
# Building model
# See https://easychair.org/publications/open/4cGr
# ==============================================================================

using MAT, Reachability, MathematicalSystems, SX
using SparseArrays, LinearAlgebra

# ==============================
# Load model
# ==============================
file = @relpath "SX/Building_more_decimals.xml"
H = readsxmodel(file, ST=ConstrainedLinearControlContinuousSystem)

# ===================
# Time-varying input
# ===================

n = size(H.modes[1].A, 1)-1 # the sx model has "time" as a state variable
@assert n == 48
A = H.modes[1].A[1:n, 1:n]
B = hcat(H.modes[1].B[1:n, 1])
X = nothing
U = Hyperrectangle(low=[0.8], high=[1.0])
S = ConstrainedLinearControlContinuousSystem(A, B, X, U)

# specify initial states
center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
X0 = Hyperrectangle(center_X0, radius_X0)
build_TV = InitialValueProblem(S, X0)

# specifications
pBDS01 = SafeStatesProperty(
    LinearConstraint(sparsevec([25], [1.0], n), 0.0051)  # x25 <= 0.0051
    )

time_horizon = 20.0

# ===================
# Constant input
# ===================

A = Reachability.add_dimension(A) # add an extra zero row and column
S = LinearContinuousSystem(A)
X0 = X0 × U
build_CONST = InitialValueProblem(S, X0)

# specifications
pBLDC01 = SafeStatesProperty(
    LinearConstraint(sparsevec([25], [1.0], n+1), 0.0051)  # x25 <= 0.0051
    )

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module MonotoneMaps

# Write your package code here.

import NLopt

#import IJulia # for notebook examples.

include("../src/endomorphisms/piece_wise_linear.jl")
include("../src/endomorphisms/composite_sigmoid.jl")

include("../src/fit/fit_to_linear.jl")

include("../src/utils.jl")

export getpiecewiselines,
    evalpiecewise2Dlinearfunc,
    evalinversepiecewise2Dlinearfunc,

    getendomorphismpiecewiselinear,
    getcompactsigmoidparameters,
    evalcompositelogisticprobit

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

# Test only on 64-bit, automated installation of HDF5 doesn't seem to work
# properly on 32-bit on CI system:
if Int == Int64
    import HDF5

    @testset "hdf5" begin
        mktempdir() do tmp_datadir
            results_filename = joinpath(tmp_datadir, "results.hdf5")
            samples = bat_sample(BAT.example_posterior(), MCMCSampling(nsteps = 1000, strict = false)).result
            bat_write(results_filename, samples)
            samples2 = bat_read(results_filename).result
            @test samples == samples2
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# TODO PR: Remove this line; this is a dummy change to trigger CI builds.

name = "Fontconfig"
version = v"2.13.1"

# Collection of sources required to build FriBidi
sources = [
    ArchiveSource("https://www.freedesktop.org/software/fontconfig/release/fontconfig-$(version).tar.bz2",
                  "f655dd2a986d7aa97e052261b36aa67b0a64989496361eca8d604e6414006741"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/fontconfig-*/

apk add gperf

# Ensure that `${prefix}/include` is..... included
export CPPFLAGS="-I${prefix}/include"

FLAGS=()
if [[ "${target}" == *-linux-* ]] || [[ "${target}" == *-freebsd* ]]; then
    FLAGS+=(--with-add-fonts="/usr/local/share/fonts")
    if [[ "${target}" == *-linux-* ]]; then
        FLAGS+=(--with-cache-dir="/var/cache/fontconfig")
    elif [[ "${target}" == *-freebsd* ]]; then
        FLAGS+=(--with-cache-dir="/var/db/fontconfig")
    fi
elif [[ "${target}" == *-apple-* ]]; then
    FLAGS+=(--with-add-fonts="/System/Library/Fonts,/Library/Fonts,~/Library/Fonts,/System/Library/Assets/com_apple_MobileAsset_Font4,/System/Library/Assets/com_apple_MobileAsset_Font5")
fi

atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0001-fix-config-linking.all.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0002-fix-mkdir.mingw.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0004-fix-mkdtemp.mingw.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0005-fix-setenv.mingw.patch"
autoreconf
./configure --prefix=$prefix --build=${MACHTYPE} --host=$target --disable-docs "${FLAGS[@]}"

# Disable tests
sed -i 's,all-am: Makefile $(PROGRAMS),all-am:,' test/Makefile

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libfontconfig", :libfontconfig),
    ExecutableProduct("fc-cache", :fc_cache),
    ExecutableProduct("fc-cat", :fc_cat),
    ExecutableProduct("fc-conflist", :fc_conflist),
    ExecutableProduct("fc-list", :fc_list),
    ExecutableProduct("fc-match", :fc_match),
    ExecutableProduct("fc-pattern", :fc_pattern),
    ExecutableProduct("fc-query", :fc_query),
    ExecutableProduct("fc-scan", :fc_scan),
    ExecutableProduct("fc-validate", :fc_validate),
    FileProduct("etc/fonts/fonts.conf", :fonts_conf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("FreeType2_jll"),
    Dependency("Bzip2_jll"),
    Dependency("Zlib_jll"),
    Dependency("Libuuid_jll"),
    Dependency("Expat_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               init_block = """
ENV["FONTCONFIG_FILE"] = get(ENV, "FONTCONFIG_FILE", fonts_conf)
    ENV["FONTCONFIG_PATH"] = get(ENV, "FONTCONFIG_PATH", dirname(ENV["FONTCONFIG_FILE"]))
""")


text:
@testset "SArray" begin
    @testset "Inner Constructors" begin
        @test SArray{(1,),Int,1,1}((1,)).data === (1,)
        @test SArray{(1,),Float64,1,1}((1,)).data === (1.0,)
        @test SArray{(2,2),Float64,2,4}((1, 1.0, 1, 1)).data === (1.0, 1.0, 1.0, 1.0)

        # Bad input
        @test_throws Exception SArray{(1,),Int,1,1}()
        @test_throws Exception SArray{(2,),Int,1,2}((1,))
        @test_throws Exception SArray{(1,),Int,1,1}(())

        # Bad parameters
        @test_throws Exception SArray{(1,),Int,1,2}((1,))
        @test_throws Exception SArray{(1,),Int,2,1}((1,))
        @test_throws Exception SArray{(1,),1,1,1}((1,))
        @test_throws Exception SArray{(2,),Int,1,1}((1,))
    end

    @testset "Outer constructors and macro" begin
        @test SArray{(1,),Int,1}((1,)).data === (1,)
        @test SArray{(1,),Int}((1,)).data === (1,)
        @test SArray{(1,)}((1,)).data === (1,)

        @test SArray{(2,2),Int,2}((1,2,3,4)).data === (1,2,3,4)
        @test SArray{(2,2),Int}((1,2,3,4)).data === (1,2,3,4)
        @test SArray{(2,2)}((1,2,3,4)).data === (1,2,3,4)

        @test ((@SArray [1])::SArray{(1,)}).data === (1,)
        @test ((@SArray [1,2])::SArray{(2,)}).data === (1,2)
        @test ((@SArray Float64[1,2,3])::SArray{(3,)}).data === (1.0, 2.0, 3.0)
        @test ((@SArray [1 2])::SArray{(1,2)}).data === (1, 2)
        @test ((@SArray Float64[1 2])::SArray{(1,2)}).data === (1.0, 2.0)
        @test ((@SArray [1 ; 2])::SArray{(2,1)}).data === (1, 2)
        @test ((@SArray Float64[1 ; 2])::SArray{(2,1)}).data === (1.0, 2.0)
        @test ((@SArray [1 2 ; 3 4])::SArray{(2,2)}).data === (1, 3, 2, 4)
        @test ((@SArray Float64[1 2 ; 3 4])::SArray{(2,2)}).data === (1.0, 3.0, 2.0, 4.0)

        @test ((@SArray [i*j*k for i = 1:2, j = 2:3, k = 3:4])::SArray{(2,2,2)}).data === (6, 12, 9, 18, 8, 16, 12, 24)
        @test ((@SArray Float64[i*j*k for i = 1:2, j = 2:3, k =3:4])::SArray{(2,2,2)}).data === (6.0, 12.0, 9.0, 18.0, 8.0, 16.0, 12.0, 24.0)

        @test (ex = macroexpand(:(@SArray [1 2; 3])); isa(ex, Expr) && ex.head == :error)

        @test ((@SArray zeros(2,2,1))::SArray{(2,2,1), Float64}).data === (0.0, 0.0, 0.0, 0.0)
        @test ((@SArray ones(2,2,1))::SArray{(2,2,1), Float64}).data === (1.0, 1.0, 1.0, 1.0)
        @test ((@SArray eye(2))::SArray{(2,2), Float64}).data === (1.0, 0.0, 0.0, 1.0)
        @test ((@SArray eye(2,2))::SArray{(2,2), Float64}).data === (1.0, 0.0, 0.0, 1.0)
        @test isa(@SArray(rand(2,2,1)), SArray{(2,2,1), Float64})
        @test isa(@SArray(randn(2,2,1)), SArray{(2,2,1), Float64})

        @test ((@SArray zeros(Float32, 2, 2, 1))::SArray{(2,2,1),Float32}).data === (0.0f0, 0.0f0, 0.0f0, 0.0f0)
        @test ((@SArray ones(Float32, 2, 2, 1))::SArray{(2,2,1),Float32}).data === (1.0f0, 1.0f0, 1.0f0, 1.0f0)
        @test ((@SArray eye(Float32, 2))::SArray{(2,2), Float32}).data === (1.0f0, 0.0f0, 0.0f0, 1.0f0)
        @test ((@SArray eye(Float32, 2, 2))::SArray{(2,2), Float32}).data === (1.0f0, 0.0f0, 0.0f0, 1.0f0)
        @test isa(@SArray(rand(Float32, 2, 2, 1)), SArray{(2,2,1), Float32})
        @test isa(@SArray(randn(Float32, 2, 2, 1)), SArray{(2,2,1), Float32})
    end

    @testset "Methods" begin
        m = @SArray [11 13; 12 14]

        @test isimmutable(m) == true

        @test m[1] === 11
        @test m[2] === 12
        @test m[3] === 13
        @test m[4] === 14

        @test Tuple(m) === (11, 12, 13, 14)

        @test size(m) === (2, 2)
        @test size(typeof(m)) === (2, 2)
        @test size(SArray{(2,2),Int,2}) === (2, 2)
        @test size(SArray{(2,2),Int}) === (2, 2)
        @test size(SArray{(2,2)}) === (2, 2)

        @test size(m, 1) === 2
        @test size(m, 2) === 2
        @test size(typeof(m), 1) === 2
        @test size(typeof(m), 2) === 2

        @test length(m) === 4

        @test_throws Exception m[1] = 1
    end
end


text:
# ---
# title: 1424. Diagonal Traverse II
# id: problem1424
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/diagonal-traverse-ii/description/>
# hidden: true
# ---
#
# Given a list of lists of integers, `nums`, return all elements of `nums` in
# diagonal order as shown in the below images.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/08/sample_1_1784.png)**
#
#
#
#     Input: nums = [[1,2,3],[4,5,6],[7,8,9]]
#     Output: [1,4,2,7,5,3,8,6,9]
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/08/sample_2_1784.png)**
#
#
#
#     Input: nums = [[1,2,3,4,5],[6,7],[8],[9,10,11],[12,13,14,15,16]]
#     Output: [1,6,2,8,7,3,9,4,12,10,5,13,11,14,15,16]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [[1,2,3],[4],[5,6,7],[8],[9,10,11]]
#     Output: [1,4,2,5,3,8,6,9,7,10,11]
#
#
# **Example 4:**
#
#
#
#     Input: nums = [[1,2,3,4,5,6]]
#     Output: [1,2,3,4,5,6]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `1 <= nums[i].length <= 10^5`
#   * `1 <= nums[i][j] <= 10^9`
#   * There at most `10^5` elements in `nums`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Mimi
using CSV
using DataFrames

# specify model settings
function set_globalbools()
    global use_variability = false

    # set random seed to have similar variability development in the base and the marginal model.
    # set variability seed.
    if use_variability
        global varseed = rand(1:1000000000000)
    end

    global use_linear = false
    global use_logburke = false
    global use_logpopulation = false
    global use_logwherepossible = true
end

# set global values for technical configuration options
set_globalbools()

# get main_model file
include("main_model_annual.jl")
include("mcs_annual.jl")
include("compute_scc_annual.jl")

if !@isdefined samplesize
    samplesize = 50000
end

for scenario in ["RCP1.9 & SSP1", "RCP2.6 & SSP1", "RCP4.5 & SSP2", "RCP8.5 & SSP5"]
    model = "PAGE-ANN"
    # define model, default settings: getpage(NDCs scenario, permafrost, no sea-ice, no page09damages)
    m = getpage(scenario, true, true)
    # run model
    run(m)

    # open up Explorer UI, for visual exploration of the variables
    # explore(m)

    # do general monte carlo simulation and save the output
    do_monte_carlo_runs(samplesize, scenario, joinpath(@__DIR__, "../../output", scenario, model))

    # get the social cost of carbon for the Monte Carlo simulations and save the output
    sccs = compute_scc_mcs(m, samplesize, year=2020)
    # store results in DataFrame
    df = DataFrame(Any[fill(model, samplesize), fill(scenario, samplesize), sccs], [:ModelName, :ScenarioName, :SCC])
    # write out to csv
    DIR = joinpath(@__DIR__, "..", "../output")
    CSV.write(joinpath(DIR, string(model, "_", scenario, "_", "scc.csv")), df)
end


text:
module OcalAPI

using CodecZlib
using HTTP
using Ipopt
using JSON
using JuMP
using Logging
using MLLabelUtils
using MLKernels
using OneClassActiveLearning
using Random
using Sockets
using SVDD

DEFAULT_SOLVER = with_optimizer(Ipopt.Optimizer, print_level=0)
SERVER_HOST = "OCAL_API_HOST" in keys(ENV) ? ENV["OCAL_API_HOST"] : ip"0.0.0.0"
SERVER_PORT = "OCAL_API_PORT" in keys(ENV) ? parse(Int, ENV["OCAL_API_PORT"]) : 8081

include("http_utils.jl")
include("parse_parameters.jl")
include("run_ocal.jl")
include("webserver.jl")

export
    SERVER_HOST, SERVER_PORT,
    run_ocal, start_webserver

end


text:
module MockServices

using SimpleTasks.Types
using SimpleTasksTests.Utils.TestTasks

import SimpleTasks.Services.Queue
import SimpleTasks.Services.Bucket
import SimpleTasks.Services.Cache
import SimpleTasks.Services.Datasource

export MockBucketService
type MockBucketService <: BucketService
    mockFiles::Dict{AbstractString, Any}
end
MockBucketService() = MockBucketService(Dict())

# local_file can also be a file location but we're not testing that for now
function Bucket.download(bucket::MockBucketService, key::ASCIIString,
        local_file::Union{ASCIIString, IO, Void})
    data = bucket.mockFiles[key]
    if local_file != nothing
        write(local_file, readbytes(data))
        return local_file
    end
    seekstart(data)

    return IOBuffer(data)
end
# local_file can also be a file location but we're not testing that for now
function Bucket.upload(bucket::MockBucketService,
    local_file::Union{ASCIIString, IO}, key::ASCIIString)
    seekstart(local_file)
    data = IOBuffer()
    write(data, readbytes(local_file))
    seekstart(data)
    bucket.mockFiles[key] = data
end

function Bucket.delete(bucket::MockBucketService, key::AbstractString)
    delete!(bucket.mockFiles, key)
end

export MockQueueService
type MockQueueService <: QueueService
end
function Queue.pop_message(mock_queue::MockQueueService)
end

export MockCacheService
type MockCacheService <: CacheService
    mockValues::Dict{AbstractString, IO}
end
MockCacheService() = MockCacheService(Dict())
function Cache.haskey(cache::MockCacheService, key::AbstractString)
    return haskey(cache.mockValues, key)
end
function Cache.put!(cache::MockCacheService, key::AbstractString,
        value_buffer::IO)
    data = IOBuffer()
    write(data, readbytes(value_buffer))
    seekstart(data)
    cache.mockValues[key] = data
end
function Cache.get(cache::MockCacheService, key::AbstractString)
    if haskey(cache.mockValues, key)
        data = cache.mockValues[key]
        seekstart(data)
        return data
    else
        return nothing
    end
end
function Cache.delete!(cache::MockCacheService, key::AbstractString)
    delete!(cache.mockValues, key)
end
function Cache.clear!(cache::MockCacheService)
    empty!(cache.mockValues)
end

#This mock datasource has no cache!
export MockDatasourceService
type MockDatasourceService <: DatasourceService
    mockSource::Dict{AbstractString, IO}
end
MockDatasourceService() = MockDatasourceService(Dict())
function Datasource.get(datasource::MockDatasourceService,
        key::AbstractString; override_cache::Bool=false)
    if haskey(datasource.mockSource, key)
        data = datasource.mockSource[key]
        seekstart(data)
        return data
    else
        return nothing
    end
end
function Datasource.put!(datasource::MockDatasourceService, key::AbstractString,
        new_value::Union{IO, Void}=nothing; only_cache::Bool=false)
    if new_value != nothing
        data = IOBuffer()
        write(data, readbytes(new_value))
        seekstart(data)
        datasource.mockSource[key] = data
    end
end
function Datasource.delete!(datasource::MockDatasourceService,
    key::AbstractString; only_cache::Bool=false)
    delete!(datasource.mockSource, key)
end

function Datasource.clear_cache(datasource::MockDatasourceService)
end

end # module MockServices

max_forks_repo_name: seung-lab/SimpleTasks


text: Dict{String, Any}("ints-and-floats" => Any[Dict{String, Any}("value" => "1", "type" => "integer"), Dict{String, Any}("value" => "1.1", "type" => "float")])


text:
# sample from a double well in 1D using ABOBA
using Printf
using Random
using LinearAlgebra
using ForwardDiff
using Statistics
using BasicMD

include("potentials.jl")

β = Float32(5.0);
γ = Float32(1.5);
M = Float32(2.1);
q₀ = Float32.([-1.0]);
p₀ = Float32.([0.0]);
Δt = Float32(1e-1);
n_iters = 10^4;
n_save_iters=10;

V = x->DoubleWell(x);
cfg = ForwardDiff.GradientConfig(V, q₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

# second moment and energy obsrevables
f₁ = x-> x[1][1]^2;
f₂ = x-> V(x[1][1]);
observables = (f₁, f₂);

sampler = ABOBA(gradV!, β, γ, M, Δt);

Random.seed!(100);
Xvals = sample_trajectory([q₀, p₀], sampler, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
f₁_estimate = mean(f₁.(Xvals));
f₂_estimate = mean(f₂.(Xvals));
@printf("f₁ estimate with %d samples: %g\n",length(Xvals), f₁_estimate);
@printf("f₂ estimate with %d samples: %g\n", length(Xvals), f₂_estimate);
Random.seed!(100);
observable_samples = sample_observables(Float32, [q₀, p₀], sampler, observables, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
@printf("f₁ estimate with %d samples: %g\n",length(observable_samples), mean(observable_samples[1,:]));
@printf("f₂ estimate with %d samples: %g\n", length(observable_samples), mean(observable_samples[2,:]));

max_stars_repo_name: gideonsimpson/JuBasicMD


text:
#BLAS.set_num_threads(6)

using TensorOperations
using JLD
using LinearMaps

include("src/PlotTools.jl")

include("src/EDutils.jl")
import EDutils
ED = EDutils

include("src/LatCFT.jl")
import LatCFT
LCFT = LatCFT

modelname = "ising";

λ=1.0; hx = 0.0; hz = 1.0
ham_comp = ED.ham_ising_comp
ham = ED.ham_ising
h = reshape(ham(λ, hz, hx), (2,2,2,2))
h_comp = ham_comp(λ, hz, hx)

nev = 41

eV_data = Dict{Int,Any}()
ed_data = Dict{Int,Any}()
Hs_data = Dict{Int,Any}()

datapath = "data/"

ED.ed!(eV_data, collect(4:1:18), h, nev)

save("$(datapath)$(modelname)_eV.jld", "eV", eV_data)

for (N, eVd) in eV_data
    ed_data[N] = (eVd[1:2]...)
end

Hs_data = Dict{Int,Any}()
LCFT.Hns_ed!(Hs_data, eV_data, h_comp, 2, 2)

LCFT.normalise_Hn!(Hs_data, ed_data)

save("$(datapath)$(modelname)_LCFTdata.jld", "ev", ed_data, "Hs", Hs_data)

max_stars_repo_name: amilsted/spinchain-virasoro-demo

max_issues_repo_name: amilsted/spinchain-virasoro-demo

max_forks_repo_name: amilsted/spinchain-virasoro-demo


text:
###
# UI Oracle
###

export UIOracle

type UIOracle <: AbstractOracle
    lbounds::Vector{Float64}
    ubounds::Vector{Float64}
    data_sort::Matrix{Float64}
    log_eps::Float64

    # Cutting plane algorithm
    qL::Vector{Float64}
    qR::Vector{Float64}
    cut_tol::Float64  ##defaults to 1e-6

    # Other options
    debug_printcut::Bool
end

#preferred constructor
function UIOracle(data, lbounds, ubounds, eps, alpha; cut_tol = 1e-6)
    N, d = size(data)
    data_sort = sort_data_cols(data)
    Gamma = KSGamma(alpha, N)
    qL, qR = gen_ql_qr(N, Gamma)
    UIOracle(vec(lbounds), vec(ubounds), data_sort, log(1/eps), qL, qR, 1e-6, false)
end

#returns bool, and ustar for degen case
is_degen(d, Gamma, log_eps) = d * log(1/Gamma) <= log_eps
degen_case(xs, lbounds::Vector{Float64}, ubounds::Vector{Float64}) = [xs[i] >= 0 ? ubounds[i] : lbounds[i] for i =1:length(xs)]

function gen_ql_qr(N::Int, Gamma)
    qL = zeros(Float64, N+2)
    qL[1] = Gamma
    qL[2:floor(Int, N * (1-Gamma)) + 1] = 1/N
    qL[floor(Int, N * (1-Gamma)) + 2] = 1-sum(qL)
    @assert (abs(sum(qL)-1) <= 1e-10) "QL not normalized $(sum(qL))"
    return qL, qL[N+2:-1:1]
end

function suppFcnUI(xs, data, lbounds, ubounds, log_eps, Gamma; cut_sense=:Max, lam_min=1e-8, lam_max=1e2, xtol=1e-8)
    data_sort = sort_data_cols(data)
    qL, qR = gen_ql_qr(size(data_sort, 1), Gamma)
    suppFcnUI(xs, data_sort, lbounds, ubounds, log_eps, qL, qR, cut_sense, lam_min, lam_max, xtol)
end

#returns zstar, ustar
function suppFcnUI(xs, data_sort, lbounds, ubounds,
                   log_eps, qL::Vector{Float64}, qR::Vector{Float64},
                   cut_sense, lam_min=1e-8, lam_max = 1e2, xtol=1e-8)
    sgn_flip = 1
    if cut_sense == :Min
        xs = -xs
        sgn_flip = -1
    end

    const N = size(data_sort, 1)
    const d = size(data_sort, 2)
    const Gamma = qL[1]

    if is_degen(d, Gamma, log_eps)
        ustar = degen_case(xs, lbounds, ubounds)
        return sgn_flip*dot(xs, ustar), ustar
    end

    #extend data with bounds, pre-condition for stability
    xdata = [lbounds'; data_sort; ubounds']
    cnst_term = 0.0
    for i = 1:d
        shift = xs[i] > 0.0 ? ubounds[i] : lbounds[i]
        cnst_term += xs[i]*shift
        for j = 1:N+2
            xdata[j, i] = xs[i]*(xdata[j, i] - shift)
        end
    end

    #objective for line-search
    #VG should you type annotate log_eps?
    function f(lam::Float64)
        val_out = lam*log_eps + cnst_term
        for i =1:d
            q = xs[i] > 0 ? qR : qL
            log_inner = 0.0
            for j = 1:N+2
                log_inner += q[j] * exp(xdata[j, i]/lam)
            end
            val_out += lam*log(log_inner)
        end
        val_out
    end

    res = Optim.optimize(f, lam_min, lam_max)
    !res.converged && error("Lambda linesearch did not converge")
    lamstar = res.minimum

    #reconstruct optimal sol.
    ustar = zeros(Float64, d)
    for i = 1:d
        if xs[i] >= 0
            qstar = qR .* exp(xdata[:, i]/lamstar)
        else
            qstar = qL .* exp(xdata[:, i]/lamstar)
        end
        qstar /= sum(qstar)

        #use the unconditioned data to construct sol
        for j= 2:N+1
            ustar[i] += qstar[j]*data_sort[j-1, i]
        end
        ustar[i] += qstar[1]*lbounds[i] + qstar[N+2]*ubounds[i]
    end

    #Not checked for performance, but should always be true
    # println(res.f_minimum, "\t", dot(ustar, xs))
    #@assert abs(res.f_minimum - dot(ustar, xs)) <= 1e-6

    sgn_flip*dot(ustar, xs), ustar
end

#preferred interface
suppFcn(xs, w::UIOracle, cut_sense) = suppFcnUI(xs, w.data_sort, w.lbounds, w.ubounds, w.log_eps, w.qL, w.qR, cut_sense)

function setup(w::UIOracle, rm::Model, prefs)
    # Extract preferences we care about
    w.debug_printcut = get(prefs, :debug_printcut, false)
    w.cut_tol        = get(prefs, :cut_tol, w.cut_tol)
    rd = JuMPeR.getRobust(rm)
    @assert (rd.numUncs == size(w.data_sort, 2)) "Num Uncertainties $(rd.numUncs) doesn't match columns in data $(size(w.data_sort, 2))"
    @assert (length(rd.uncertaintyset) == 0) "Auxiliary constraints on uncertainties not yet supported"
end


text:
# Autogenerated wrapper script for soapysdr_jll for armv6l-linux-musleabihf-cxx11
export libsoapysdr, soapysdrutil

JLLWrappers.@generate_wrapper_header("soapysdr")
JLLWrappers.@declare_library_product(libsoapysdr, "libSoapySDR.so.0.8")
JLLWrappers.@declare_executable_product(soapysdrutil)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsoapysdr,
        "lib/libSoapySDR.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        soapysdrutil,
        "bin/SoapySDRUtil",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("../common.jl")
build_julia(v"1.3.1")


text:
using Test
using Distributions
import Distributions: normpdf, normcdf, normlogpdf, normlogcdf

@testset "SkewNormal" begin
    @test_throws ArgumentError SkewNormal(0.0, 0.0, 0.0)
    d1 = SkewNormal(1, 2, 3)
    d2 = SkewNormal(1.0f0, 2, 3)
    @test partype(d1) == Float64
    @test partype(d2) == Float32
    @test params(d1) == (1.0, 2.0, 3.0)
    # Azzalini sn: sprintf("%.17f",dsn(3.3, xi=1, omega=2, alpha=3))
    @test pdf(d1, 3.3) ≈ 0.20587854616839998
    @test minimum(d1) ≈ -Inf
    @test maximum(d1) ≈  Inf
    @test logpdf(d1, 3.3) ≈ log(pdf(d1, 3.3))
    ## cdf and quantile: when we get Owen's T
    #@test cdf(d, 4.5) ≈ 1.0 #when we get Owen's T
    @test mean(d1) ≈ 2.513879513212096
    @test std(d1) ≈ 1.306969326142243
    #
    d0 = SkewNormal(0.0, 1.0, 0.0)
    @test SkewNormal() == d0
    d3 = SkewNormal(0.5, 2.2, 0.0)
    d4 = Normal(0.5, 2.2)
    #
    @test pdf(d3, 3.3) == Distributions.pdf(d4, 3.3)
    @test pdf.(d3, 1:3) == Distributions.pdf.(d4, 1:3)
    a = mean(d3), var(d3), std(d3)
    b = Distributions.mean(d4), Distributions.var(d4), Distributions.std(d4)
    @test a == b
    @test skewness(d3) == Distributions.skewness(d4)
    @test kurtosis(d3) == Distributions.kurtosis(d4)
    @test mgf(d3, 2.25) == Distributions.mgf(d4, 2.25)
    @test cf(d3, 2.25) == Distributions.cf(d4, 2.25)
end

# R code using Azzalini's package sn
#library("sn")
#sprintf("%.17f", dsn(3.3, xi=1, omega=2, alpha=3))
#f1 <- makeSECdistr(dp=c(1,2,3), family="SN", name="First-SN")
#sprintf("%.15f", mean(f1))
#sprintf("%.15f", sd(f1))


text:
function internalsolve(mdp::MDP, smcts::SerialMCTS)
  return smcts
end

function smctsaction!(mdp::MDP, smcts::SerialMCTS, actualstate::Vector)
  for iter in 1:smcts.niter
    simulate(mdp, smcts, smcts.maxdepth, actualstate)
  end
  statedim = length(mdp.statemap)
  actionargs = mdp.reward.argnames[1 + statedim:end]
  iaction = indmax(smcts.tree[actualstate].qval)
  return getvar(smcts.actiongrid, mdp.actionmap, actionargs, iaction)
end

function simulate(mdp::MDP, smcts::SerialMCTS, depth::Int64, actualstate::Vector)

  if depth == 0
    return 0.0
  end

  statedim = length(mdp.statemap)
  stateargs = mdp.reward.argnames[1:statedim]

  state = getgridstate(smcts.stategrid, mdp.statemap, stateargs, actualstate)

  # run rollout and add statistics to tree if state unexplored
  if !haskey(smcts.tree, actualstate)

    smcts.tree[deepcopy(actualstate)] = StateNode(length(smcts.actiongrid))
    return rollout(mdp, smcts, depth, state)

  else

    actionargs = mdp.reward.argnames[1 + statedim:end]
    snode = smcts.tree[actualstate]

    # pick action using UCT
    iaction = indmax(snode.qval + smcts.exex * real(sqrt(complex(log(sum(snode.n)) ./ snode.n))))
    action = ind2x(smcts.actiongrid, iaction)

    # transition to next state
    nextstate, reward = sampletransition(mdp, smcts, state, action)

    # update n and qval in state node
    actualstate = getvar(mdp.statemap, stateargs, state)
    actualaction = getvar(mdp.actionmap, actionargs, action)
    actualnextstate = getvar(mdp.statemap, stateargs, nextstate)

    q = reward + smcts.discount * simulate(mdp, smcts, depth - 1, actualnextstate)

    snode.n[iaction] += 1
    snode.qval[iaction] += (q - snode.qval[iaction]) / snode.n[iaction]

    return q

  end

end

function rollout(mdp::MDP, smcts::SerialMCTS, depth::Int64, state::Vector{Float64})
  if depth == 0
    return 0.0
  else
    iaction = rand(1:length(smcts.actiongrid))  # uniformly random rollout
    action = ind2x(smcts.actiongrid, iaction)
    nextstate, reward = sampletransition(mdp, smcts, state, action)
    return reward + smcts.discount * rollout(mdp, smcts, depth - 1, nextstate)
  end
end

function sampletransition(
    mdp::MDP,
    smcts::SerialMCTS,
    state::Vector{Float64},
    action::Vector{Float64})

  if length(mdp.transition.argnames) == length(mdp.statemap) + length(mdp.actionmap)
    return sampleset(mdp, smcts, state, action)
  elseif length(mdp.transition.argnames) > length(mdp.statemap) + length(mdp.actionmap)
    return sampleprob(mdp, smcts, state, action)
  else
    error(string(
      "unkown transition function of neither T(s,a,s') nor T(s,a) type\n",
      "number of transition arguments: ",
      length(mdp.transition.argnames), "\n",
      "number of state and action variables",
      length(mdp.statemap) + length(mdp.actionmap)))
  end
end

function sampleset(
    mdp::MDP,
    smcts::SerialMCTS,
    state::Vector{Float64},
    action::Vector{Float64})

  statedim = length(mdp.statemap)
  stateargs = mdp.reward.argnames[1:statedim]
  actionargs = mdp.reward.argnames[1 + statedim:end]

  actualstate = getvar(mdp.statemap, stateargs, state)
  actualaction = getvar(mdp.actionmap, actionargs, action)

  nextstateprobs = mdp.transition.fn(actualstate..., actualaction...)
  weightvec = WeightVec([nextstateprob[2]::Float64
    for nextstateprob in nextstateprobs])

  nextstateprob = sample(nextstateprobs, weightvec)
  nextstate = getgridstate(smcts.stategrid, mdp.statemap, stateargs, nextstateprob[1])
  reward = mdp.reward.fn(actualstate..., actualaction...)

  return nextstate, reward

end

function sampleprob(
    mdp::MDP,
    smcts::SerialMCTS,
    state::Vector{Float64},
    action::Vector{Float64})

  statedim = length(mdp.statemap)
  stateargs = mdp.reward.argnames[1:statedim]
  nstates = length(smcts.stategrid)

  actualstate = getvar(mdp.statemap, stateargs, state)
  actualaction = getvar(mdp.actionmap, actionargs, action)

  probs = zeros(nstates)
  for inextstate in 1:nstates
    actualnextstate = getvar(smcts.stategrid, mdp.statemap, stateargs, inextstate)
    probs[inextstate] = mdp.transition.fn(actualstate..., actualaction..., actualnextstate...)
  end
  weightvec = WeightVec(probs)

  nextstate = ind2x(smcts.stategrid, sample(1:nstates, weightvec))
  reward = mdp.reward.fn(actualstate..., actualaction...)

  return nextstate, reward

end

# Returns the actual variable from GridInterpolations indices
function getvar(
    grid::RectangleGrid,
    map::Dict{AbstractString, LazyVar},
    argnames::Vector{ASCIIString},
    index::Int64)

  raw = ind2x(grid, index)
  var = Array(Any, length(raw))

  for ivar in 1:length(raw)
    lazy = map[argnames[ivar]]
    if isa(lazy, RangeVar)
      var[ivar] = raw[ivar]
    elseif isa(lazy, ValuesVar)
      var[ivar] = lazy.values[Int64(raw[ivar])]
    else
      error(string(
        "unknown state/action variable definition type for ", argnames[ivar]))
    end
  end

  return var

end

# Returns the actual variable from GridInterpolations grid state
function getvar(
    map::Dict{AbstractString, LazyVar},
    argnames::Vector{ASCIIString},
    gridstate::Vector{Float64})

  var = Array(Any, length(gridstate))

  for ivar in 1:length(gridstate)
    lazy = map[argnames[ivar]]
    if isa(lazy, RangeVar)
      var[ivar] = gridstate[ivar]
    elseif isa(lazy, ValuesVar)
      var[ivar] = lazy.values[Int64(gridstate[ivar])]
    else
      error(string(
        "unknown state/action variable definition type for ", argnames[ivar]))
    end
  end

  return var

end

# Inverse of getvar; returns the GridInterpolations state for an actual state
function getgridstate(
    grid::RectangleGrid,
    map::Dict{AbstractString, LazyVar},
    argnames::Vector{ASCIIString},
    state::Vector)

  gridstate = Array(Float64, length(state))

  for i in 1:length(gridstate)
    lazy = map[argnames[i]]
    if isa(lazy, RangeVar)
      cutpoints = grid.cutPoints[i]
      gridstate[i] = cutpoints[indmin(abs(cutpoints - state[i]))]
    elseif isa(lazy, ValuesVar)
      gridstate[i] = findfirst(lazy.values, state[i])
    else
      error(string(
        "unknown state/action variable definition type for ", argnames[ivar]))
    end
  end

  return gridstate

end


text:
export WeightingOp

"""
  WeightingOp(weights::Vector{T}, rep::Int=1) where T

generates a `LinearOperator` which multiplies an input vector index-weise with `weights`

# Arguments
* `weights::Vector{T}` - weights vector
* `rep::Int=1`         - number of sub-arrays that need to be multiplied with `weights`
"""
function WeightingOp(weights::Vector{T}, rep::Int=1) where T
  weights_cat = repeat(weights,rep)
  return opDiagonal(weights_cat)
end


text:
@testset "testing `filtern!()`: " begin
	using DataFrames
	lenrow = 30;
	lencol = 8;
	maxeach = floor(lenrow*(lencol-1)/10) |> x -> convert(Int, x);
			# floor(lencol/4) |> x -> convert(Int, x);
	all_colid = collect(1:lencol);
	all_rowid = collect(1:lenrow)
	collist = ["col$i" for i in 1:lencol-1];

	pair_list = [];
	push!(pair_list, "absrow" => all_rowid);
	for colname in collist
		veci = Vector{Union{Float64, Missing, Nothing}}(undef, lenrow);
		veci[:] = randn(lenrow);
		push!(pair_list, colname => veci)
	end
	testDataFrame = DataFrame(pair_list);

	"""
	get random indices `(is, js)` for dataframe.
	"""
	function get_randind2(rrange, crange, numsamp)
	# numsamp = rand(0:5);
	idr = rand(rrange, numsamp);
	idc = rand(crange, numsamp);
	return [ (idr, idc[i]) for (i, idr) in enumerate(idr)]
	end

	function replace_by!(df, listidridc, sth)
		for ind2 in listidridc
			df[ind2...] = sth;
		end
	end

	function rand_replace_all!(df)
		rrange = 1:lenrow;
		crange = 2:lencol;
		ismissing_id = get_randind2(rrange, crange, rand(1:maxeach));
		isnothing_id = get_randind2(rrange, crange, rand(1:maxeach));
		isnan_id = get_randind2(rrange, crange, rand(1:maxeach));
		replace_by!(df, isnan_id, NaN)
		replace_by!(df, isnothing_id, nothing)
		replace_by!(df, ismissing_id, missing)

		fnx(id) = [v[1] for v in id];

		rowid_remain = collect(rrange);
		rowid_remain = setdiff!(rowid_remain, fnx(ismissing_id));
		rowid_remain = setdiff!(rowid_remain, fnx(isnothing_id));
		rowid_remain = setdiff!(rowid_remain, fnx(isnan_id));

		return rowid_remain
	end

# test begin:

	df = copy(testDataFrame);
	rowid_remain = rand_replace_all!(df);
	for colname in collist
		filtern!(colname, df)
	end
	@test df[!,"absrow"] == rowid_remain
end