text:
using LibPQ
# using Query
using JSON

function getproperties_asdict(object::Any,include_missing::Bool)
    properties_symbols = propertynames(object) # Get the properties as a vector of Symbol
    props::Dict{Symbol,Any} =
        Dict(zip(properties_symbols, # The properties symbols
        getfield.([object],properties_symbols)) # the properties values
                    )
    # Remove the missing properties if needed
    if (!include_missing)
        props_as_vectoroftuples = dict2vectoroftuples(props)
        props_as_vectoroftuples_womissing = filter(t -> !ismissing(t[2]),
                                                   props_as_vectoroftuples)
        props = vectoroftuples2dict(props_as_vectoroftuples_womissing)
    end
    return props

end

function getpropertiesvalues(obj::Any,props::Vector{Symbol})
    values = []
    for p in props
        push!(values,getproperty(obj,p))
    end
    return values
end

function getdictvalues(dict::Dict,props::Vector{Symbol})
    # NOTE: Cannot use the following line because we want the values in the same
    #           order as the props argument
    # collect(values(filter(x -> x.first in props, dict)))

    values = []
    for p in props
        push!(values,get(dict,p,missing))

    end
    return values
end

function getdictvalues(dict::Dict,props::Vector{String})
    # NOTE: Cannot use the following line because we want the values in the same
    #           order as the props argument
    # collect(values(filter(x -> x.first in props, dict)))

    values = []
    for p in props
        push!(values,get(dict,p,missing))

    end
    return values
end

function setpropertiesvalues!(obj::Any,props::Vector{Symbol},values::Vector)
    for (p,v) in zip(props,values)
        setproperty!(obj,p,v)
    end
end

function dataframe2vector_of_namedtuples(df::DataFrame)
    result = []
    for row in eachrow(df)
        push!(result, dataframerow2namedtuple(row))
    end
    return result
end

function dataframerow2namedtuple(row::DataFrameRow)
    vals = values(row)
    cols = Symbol.(names(row))
    dict2namedtuple(Dict(zip(cols,vals)))
end

function vectoroftuples2dict(t::Vector{Tuple{T,K}}) where {T <: Any, K <: Any}
    Dict(zip(
            first.(t), # keys
            last.(t) # values
            )
        )
end

function dict2vectoroftuples(d::Dict)
    collect(zip(collect(keys(d)),collect(values(d))))
end

# https://discourse.julialang.org/t/how-to-make-a-named-tuple-from-a-dictionary/10899/10?u=tencnivel
function dict2namedtuple(d::Dict)
    (; d...)
end

dictstringkeys2symbol_helper(x) = x
dictstringkeys2symbol_helper(d::Array) = dictstringkeys2symbol_helper.(d)
dictstringkeys2symbol_helper(d::Dict) =
    Dict(Symbol(k) => dictstringkeys2symbol_helper(v) for (k, v) in d)
function dictstringkeys2symbol(d::Dict)
    dictstringkeys2symbol_helper(d)
end

dictnothingvalues2missing_helper(x) = x
dictnothingvalues2missing_helper(x::Nothing) = missing
dictnothingvalues2missing_helper(d::Array) = dictnothingvalues2missing_helper.(d)
dictnothingvalues2missing_helper(d::Dict) =
    Dict(k => dictnothingvalues2missing_helper(v) for (k, v) in d)
function dictnothingvalues2missing(d::Dict)
    dictnothingvalues2missing_helper(d)
end

# https://discourse.julialang.org/t/how-to-make-a-named-tuple-from-a-dictionary/10899/11?u=tencnivel
function namedtuple2dict(nt::NamedTuple)
    Dict{Symbol,Any}(pairs(nt))
end

function string2enum(enumType::DataType,str::SubString{String})
    string2enum(enumType,string(str))
end

function string2enum(enumType::DataType,str::String)
    if ismissing(str)
        return missing
    end

    result = try
        enum_ints = Int.(instances(enumType))
        enum_names = string.(instances(enumType))

        index_oi = findfirst(x -> x == str,enum_names)

        enum_int = enum_ints[index_oi]
        enumType(enum_int)
    catch e
        throw(DomainError("Unable to convert[$str] to Enum[$enumType]"))
    end

    return result

end

"""
"""
function string2vector_of_enums(vectorOfEnumsTypes::Type{Vector{T}},
                                str::Union{String,Missing}) where T <: Base.Enums.Enum
    enumType = eltype(vectorOfEnumsTypes)

    if ismissing(str)
        return missing
    end

    # If no element in the array return an empty array
    if str == "{}"
        return T[]
    end

    chop(str,head = 1,tail = 1) |> n -> split(n,",") |>
        n -> string.(n) |>
        n -> string2enum.(enumType,n)

end

function string2zoneddatetime(str)
    # eg. "2019-09-03T11:00:00.000Z"
    date_match_GMT =
        match(r"^([0-9]{4}-[0-9]{2}-[0-9]{2}T[0-9]{2}:[0-9]{2}:[0-9]{2}.[0-9]{3}Z)$", str)

    if !isnothing(date_match_GMT)
        # "2019-07-24T00:41:49.732Z" becomes "2019-07-24T00:41:49.732"
        date_match_remove_endingZ =
          match(r"^([0-9]{4}-[0-9]{2}-[0-9]{2}T[0-9]{2}:[0-9]{2}:[0-9]{2}.[0-9]{3})",
                str)
        utc_zdt =
            ZonedDateTime(DateTime(date_match_remove_endingZ.match),
                          TimeZone("UTC"))
        return utc_zdt
    else
        return nothing
    end

end

function int2enum(enumType::DataType,enum_int::Missing)
    return missing
end

function int2enum(enumType::DataType,enum_int::Integer)
    enumType(enum_int)
end

function enum2int(enumobject::T) where T <: Enum
    Int(enumobject)
end

function tovector(x;elementstype = missing)

    # Create vector
    if isa(x,Vector)
        result = x
    elseif isa(x,Tuple)
        result = collect(x)
    elseif isa(x,Number) || isa(x,String) || isa(x,Bool) || isa(x,Symbol)
        result = [x]
    else
        throw(ArgumentError("Unsupported type[$(typeof(x))]"))
    end

    # Convert vector
    if (!ismissing(elementstype))
        if elementstype == Symbol
            result = Symbol.(result)
        elseif elementstype == String
            result = string.(result)
        else
            throw(ArgumentError("Unsupported elementstype[$elementstype]"))
        end
    end

    return result

end

"""
    postgresql_string_array_2_string_vector(str::String)

Transform the following string "{"Deleted,. 3, Scenes a2e","Behind the Scenes",""}"
(that's what we get from postrgesql) into a string vector fo size 3
"""
function postgresql_string_array_2_string_vector(str::String)

    #TODO handle the case "{\"bla bla\",bb,cc}"

    if occursin("\"",str)

        my_matches =
            #              1            2          3          4           5
            eachmatch(r"(?:\"([^\"]*)\"|,([^\"]*),|{([^\"]*),|{([^\"]*)}|,([^\"]*)})",
            # Where 1: "Deleted,. 3, Scenes a2e" => Deleted,. 3, Scenes a2e
            #       2: ,dwed, => dwed
            #       3: {dwed, => dwed
            #       4: {dwed} => dwed
            #       5: ,dwed} => dwed
                       str
                      )
        return map(m ->filter(x -> !isnothing(x), m.captures) |> n -> (string ∘ first)(n),
            collect(my_matches))
      else
          return  str[2:length(str)-1] |>
                   n -> split(n,",") |>
                   n -> string.(n)
      end

end

# This method exists so that we can call it even without having to test if the
#   argument is a Union
function get_nonmissing_typeof_uniontype(arg::Any)
    return arg
end

# NOTE: We cannot say that this method returns a DataType (::DataType) because
#         arrays are not datatypes
function get_nonmissing_typeof_uniontype(arg::Union)
    if arg.a != Missing && arg.a != Nothing
        return get_nonmissing_typeof_uniontype(arg.a)
    else
        return get_nonmissing_typeof_uniontype(arg.b)
    end
end

function diff_dict(old_dict::Dict, new_dict::Dict)

    old_dict_keys = collect(keys(old_dict))
    new_dict_keys = collect(keys(new_dict))

    diff_result = Dict()

    # Loop over the old dictionary keys
    for (k,v) in old_dict
        # Values that have changed without being nulled
        if k in new_dict_keys

            # Dirty quickfix to handle byte arrays
            if isa(v,Vector{UInt8})
                diff_result[k] = (old = "old bytea", new = "new bytea")
                continue
            end

            # For vectors of things that are not IEntity it is possible that
            #   we have a missing value although the property has been loaded
            if ((isa(v, Vector) && !isa(v,Vector{<:PostgresORM.IEntity}))
             || (isa(new_dict[k], Vector) && !isa(new_dict[k],Vector{<:PostgresORM.IEntity})))
                if old_dict[k] !== new_dict[k]
                    diff_result[k] = (old = old_dict[k], new = new_dict[k])
                end
                continue
            end

            # Sanity check, we do not allow comparing a vector with missing
            #   because it's a symptom of the user not comparing objects that
            #   can be compared. Recall that when we load a vector property from
            #   the database, it will be at least an empty array never missing.
            if (isa(v, Vector) && ismissing(new_dict[k]))
                throw(DomainError(
                    "You are trying to compare Missing with a Vector for "
                    *"property[$k]. Missing property is in the new version."))
            end
            if (isa(new_dict[k], Vector) && ismissing(v))
                throw(DomainError(
                    "You are trying to compare Missing with a Vector for "
                    *"property[$k]. Missing property is in the old version."))
            end

            if (ismissing(old_dict[k]) && ismissing(new_dict[k]))
                # do nothing
            elseif (ismissing(old_dict[k]) && !ismissing(new_dict[k])) ||
               (!ismissing(old_dict[k]) && ismissing(new_dict[k])) ||
               old_dict[k] != new_dict[k]
               diff_result[k] = (old = old_dict[k], new = new_dict[k])
            end
        # Values that have been nulled
        elseif  !(k in new_dict_keys)
            diff_result[k] = (old = old_dict[k], new = missing)
        end
    end

    # Loop over the new dictionary for new keys
    for k in new_dict_keys
        if !(k in old_dict_keys)
            diff_result[k] = (old = missing, new = new_dict[k])
        end
    end

    diff_result

end

function remove_spaces_and_split(str::String)
    str = replace(str, " " => "")
    return split(str,',')
end


text:
julia> n = 20
20

julia> [(x, y, z) for x = 1:n for y = x:n for z = y:n if x^2 + y^2 == z^2]
6-element Array{Tuple{Int64,Int64,Int64},1}:
 (3,4,5)
 (5,12,13)
 (6,8,10)
 (8,15,17)
 (9,12,15)
 (12,16,20)

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import .PyCall: PyObject, NpyArray, npy_type, NPY_STRING, @npyinitialize, @pycheck
import .PyCall: npy_api, NPY_ARRAY_ALIGNED, NPY_ARRAY_WRITEABLE, array2py, PyPtr

npy_type(::Type{T}) where T <: MaxLenString = NPY_STRING

function NpyArray(a::StridedArray{T}, revdims::Bool) where T <: MaxLenString
    @npyinitialize
    size_a = revdims ? reverse(size(a)) : size(a)
    strides_a = revdims ? reverse(strides(a)) : strides(a)
    p = @pycheck ccall(npy_api[:PyArray_New], PyPtr,
          (PyPtr, Cint, Ptr{Int}, Cint, Ptr{Int}, Ptr{T}, Cint, Cint, PyPtr),
          npy_api[:PyArray_Type],
          ndims(a), Int[size_a...], npy_type(T),
          Int[strides_a...] * sizeof(eltype(a)), a, sizeof(eltype(a)),
          NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE,
          C_NULL)
    return PyObject(p, a)
end

function PyObject(a::StridedArray{T}) where T <: MaxLenString
    try
        return NpyArray(a, false)
    catch
        return array2py(a) # fallback to non-NumPy version
    end
end


text:
#=
The four adjacent digits in the 1000-digit number that have the greatest product are 9 × 9 × 8 × 9 = 5832.

73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450
Find the thirteen adjacent digits in the 1000-digit number that have the greatest product. What is the value of this product?
=#

function greatest_adjacent_digit_product(adjacent_digits)
  i = 7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450
  i_digits = digits(i)
  largest_product = 1
  for x in 1:length(i_digits)-adjacent_digits
    product = 1
    for y in 1:adjacent_digits
      product *= i_digits[x+y]
    end
    largest_product = product > largest_product ? product : largest_product
  end
  largest_product
end
@time println(greatest_adjacent_digit_product(13))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module PiMath
import Showoff
import Gadfly
import Base.+, Base.-, Base.*, Base./, Base.<, Base.<=

immutable MulOfPi{T} <: Real
    f::T
end
immutable PiRange{S} <: Range{S}
    r::LinSpace{S} #This should be a Range object
end

function Showoff.showoff{T}(xs::AbstractArray{MulOfPi{T}},style=:none)
    return [string(x, "π") for x in Showoff.showoff([x.f for x in xs],style)]
end
function Showoff.showoff(xs::AbstractArray{MulOfPi},style=:none)
    return [string(x, "π") for x in Showoff.showoff([x.f for x in xs],style)]
end

Base.getindex(x::PiRange,i::Integer)=MulOfPi(getindex(x.r,i))
Base.getindex(x::PiRange,i::Colon)=PiRange(getindex(x.r,i))
Base.getindex(x::PiRange,i::Range)=PiRange(getindex(x.r,i))
Base.length(x::PiRange)=length(x.r)
Base.step(x::PiRange)=MulOfPi(step(x.r))
Base.isempty(x::PiRange)=isempty(x.r)
Base.first(x::PiRange) = MulOfPi(first(x.r))
Base.last(x::PiRange) = MulOfPi(last(x.r))
Base.start(x::PiRange) = start(x.r)
Base.done(x::PiRange, i::Int) = done(x.r,i)
Base.next(x::PiRange, i::Int) = begin a,b=next(x.r,i); (MulOfPi(a),b) end

Base.show{T<:Rational}(io::IO,x::MulOfPi{T})=print(io,x.f.num==1 ? "" : x.f.num,"π", x.f.den==1 ? "" : string("/",x.f.den))
Base.show{T<:Real}(io::IO,x::MulOfPi{T})=print(io,x.f,"π")
Base.isfinite(x::MulOfPi)=isfinite(x.f)
Base.convert{T}(::Type{MulOfPi{T}}, x::Real)=MulOfPi(convert(T,x/convert(Float64,pi)))
Base.convert(::Type{MulOfPi}, x::Real)=error("")

+(a::MulOfPi,b::MulOfPi)=MulOfPi(a.f+b.f)
-(a::MulOfPi)=MulOfPi(-a.f)
-(a::MulOfPi,b::MulOfPi)=MulOfPi(a.f-b.f)
*(a::Bool,b::Irrational{:π})=a ? zero(Float64) : convert(Float64,Irrational{:π})
*{T}(a::Bool,b::MulOfPi{T})=a ? MulOfPi(zero(T)) : MulOfPi(one(T))
*(a::MulOfPi,b::Union{Integer,AbstractFloat})=MulOfPi(a.f*b)
*(b::Union{Integer,AbstractFloat,Rational},a::MulOfPi)=MulOfPi(a.f*b)
*(a::Union{Integer,AbstractFloat,Rational},b::Irrational{:π})=MulOfPi(a)
*(b::Irrational{:π},a::Union{Integer,AbstractFloat})=MulOfPi(a)
*(a::MulOfPi,b::MulOfPi)=convert(Float64,a)*convert(Float64,b)
/(a::MulOfPi,b::Union{Integer,AbstractFloat})=MulOfPi(a.f/b)
/(a::MulOfPi,b::MulOfPi)=a.f/b.f
/(b::Irrational{:π},a::Union{Integer,AbstractFloat})=MulOfPi(inv(a))
<(a::MulOfPi,b::MulOfPi)=a.f<b.f
<=(a::MulOfPi,b::MulOfPi)=a.f<=b.f

Base.convert{S}(::Type{MulOfPi{S}},x::Irrational{:π})=MulOfPi(1)
Base.convert(::Type{Bool}, x::MulOfPi)=convert(Bool,x.f)
Base.convert{T<:AbstractFloat,S}(::Type{T},x::MulOfPi{S})=begin a,b=promote(x.f,pi);convert(T,a*b) end
Base.convert{T<:Integer,S}(::Type{T},x::MulOfPi{S})=begin a,b=promote(x.f,pi);convert(T,a*b) end
Base.convert{T<:Rational,S}(::Type{T},x::MulOfPi{S})=begin a,b=promote(x.f,pi);convert(T,a*b) end
Base.convert(::Type{MulOfPi},x::Real)=MulOfPi(x/pi)
Base.promote_rule{S}(::Type{Irrational{:π}},::Type{MulOfPi{S}})=MulOfPi{S}
Base.promote_rule{T<:AbstractFloat,S<:Any}(::Type{MulOfPi{S}},::Type{T})=T
Base.promote_rule{T<:Union{Integer,Rational},S<:Any}(::Type{MulOfPi{S}},::Type{T})=Float64

function Gadfly.optimize_ticks(x_min::MulOfPi, x_max::MulOfPi; kwargs...)
    ticks, viewmin, viewmax = Gadfly.optimize_ticks(x_min.f, x_max.f;kwargs...)
    return [MulOfPi(t) for t in ticks],MulOfPi(viewmin),MulOfPi(viewmax)
end

Base.linspace(start::MulOfPi,stop::MulOfPi,len::Real)=PiRange(linspace(start.f,stop.f,len))
Base.linspace{T<:AbstractFloat}(start::T,stop::MulOfPi,len::Real)=PiRange(linspace(start,stop.f,len))
Base.linspace{T<:Integer}(start::T,stop::MulOfPi,len::Real)=PiRange(linspace(start,stop.f,len))
Base.linspace{T<:AbstractFloat}(start::MulOfPi,stop::T,len::Real)=PiRange(linspace(start.f,stop,len))
Base.linspace{T<:Integer}(start::MulOfPi,stop::T,len::Real)=PiRange(linspace(start.f,stop,len))

#Some sin and cosine methods
Base.sin{T<:Integer}(x::MulOfPi{T})=zero(Float64)
Base.cos{T<:Integer}(x::MulOfPi{T})=iseven(x.f) ? one(Float64) : -one(Float64)

end # module


text:
using FileIO: load, save, loadstreaming, savestreaming
using Unitful

export streamFormat, LE_BOM, BE_BOM

struct streamFormat{eltype <: Number}
    nchannels::Int
    samplerate::Unitful.Frequency
    endianBOM::UInt32
end

const LE_BOM = UInt32(0x04030201)
const BE_BOM = UInt32(0x01020304)

function endianBswap(endianBOM::UInt32, x)
    if endianBOM != ENDIAN_BOM
        x = bswap(x)
    end
    return x
end

function endianBswap(endianBOM::UInt32, array::AbstractArray)
    if endianBOM != ENDIAN_BOM
        array = bswap.(array)
    end
    return array
end

function endianBswap!(endianBOM::UInt32, array::AbstractArray)
    if endianBOM != ENDIAN_BOM
        array .= bswap.(array)
    end
    return
end

function endianRead(IO, endianBOM::UInt32, ::Type{T}) where {T}
    res = read(IO, T)
    return endianBswap(endianBOM, res)
end

function endianRead!(IO, endianBOM::UInt32, array::AbstractArray)
    read!(IO, array)
    endianBswap!(endianBOM, array)
    return
end

function endianWrite(IO, endianBOM::UInt32, x)
    write(IO, endianBswap(endianBOM, x))
end

function endianWrite!(IO, endianBOM::UInt32, array::AbstractArray)
    endianBswap!(endianBOM, array)
    write(IO, array)
end

function channelRead(IO, endianBOM::UInt32, nchannels::Int, ::Type{T}, n::Int) where {T}
    data = Array{T, 2}(undef, nchannels, n)
    endianRead!(IO, endianBOM, data)
    return vec(view(data, 1:1, :))
end

function normalizedChannelRead(IO, endianBOM::UInt32, nchannels::Int, ::Type{T}, ::Type{R}, n::Int) where {T <: Number, R <: Real}
    data = Array{T, 2}(undef, nchannels, n)
    endianRead!(IO, endianBOM, data)
    if T <: Int
        return vec(view(data, 1:1, :)) .* 1.0 ./ typemax(T)
    end
    return vec(view(data, 1:1, :))
end


text:
# ____________________________________________________________________
# SpecLineLorentz
#
mutable struct SpecLineLorentz <: AbstractComponent
    norm::Parameter
    center::Parameter
    fwhm::Parameter
    voff::Parameter
    index::Vector{Int}  # optimization
    norm_integrated::Bool

    function SpecLineLorentz(center::Number)
        out = new(Parameter(1),
                  Parameter(center),
                  Parameter(3000),
                  Parameter(0),
                  Vector{Int}(), true)
        @assert center > 0
        out.norm.low = 0
        out.center.low = 0
        out.fwhm.low = 0
        out.voff.low = 0
        out.center.fixed = true
        return out
    end
end

function prepare!(comp::SpecLineLorentz, domain::Domain{1})
    comp.index = collect(1:length(domain))
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::SpecLineLorentz, x::Domain{1},
                   norm, center, fwhm, voff)
    buffer[comp.index] .= 0.
    empty!(comp.index)

    x0 = center - (voff / 3.e5) * center
    hwhm = fwhm / 3.e5 * center / 2  # Note: this is in `center` units

    X = (x .- x0) ./ hwhm
    i = findall(abs.(X) .< 20) # optimization
    append!(comp.index, i)
    buffer[i] .= norm ./ (1 .+ X[i].^2.)
    if comp.norm_integrated
        buffer[i] ./= pi * hwhm
    end
end

#=
    x = Domain(500:1:1500.)
    comp = QSFit.SpecLineLorentz(1000.)
    comp.fwhm.val = 3e4
    ceval = GFit.CompEval(comp, x)
    GFit.evaluate_cached(ceval)
    @gp x ceval.buffer ./ maximum(ceval.buffer) "w l"
=#


text:
# This file is not readable yet, is used for unit testing purposes needs to be organizedo
include("preprocess.jl")

function test_datatype(datafile)
    serve_type = "bitarray"
    batchsize = 57
    data = Data(datafile; serve_type=serve_type, batchsize=batchsize)
    return data
end

""" test the sequence field of datatype"""
function test_sequence_field(testfile)
    # Based on the property of the testdata.txt we only have sentences of length 3-4-5, corresponding
    # lengths of 5-6-7 that is being hold on to the test_sentence_sequence
    # preprocess the data
    tdata = test_datatype(testfile)
    seq5 = tdata.sequences[5]
    seq6 = tdata.sequences[6]
    seq7 = tdata.sequences[7]

    sets = icreate_sets(testfile)
    set5 = sets[1]
    set6 = sets[2]
    set7 = sets[3]

    # Not fancy but works for now
    for i=1:length(seq5)
        assert(map(x->tdata.index_to_word[x], seq5[i])[2:end-1] in set5) # ignore <s> and </s> for now
    end
    for i=1:length(seq6)
        assert(map(x->tdata.index_to_word[x], seq6[i])[2:end-1] in set6) # ignore <s> and </s> for now
    end
    for i=1:length(seq7)
        assert(map(x->tdata.index_to_word[x], seq7[i])[2:end-1] in set7) # ignore <s> and </s> for now
    end
end

""" reserves all the batches of a given length-sequence for test purposes """
function icreate_batches(sentences::Array{Any,1}, batchsize::Int, vocabsize::Int)
    batches = Array{Any, 1}()
    (d, from) = sentenbatch(sentences, 1, batchsize, vocabsize) #what we expect from sentencebatch is to give us a sequence, from is 1 in the first call
    push!(batches, d)
    while from != 1
        (d, from) = sentenbatch(sentences, from, batchsize, vocabsize)
        push!(batches, d)
    end
    if nothing in batches
        info("Removing nothing entry")
        pop!(batches)
    end
    return batches
end

""" Returns the sets that contains all the same length sentences """
function icreate_sets(testfile::AbstractString)
    # Be careful in our testfile there exists only the length of 5-6-7 sentences!
    # Here in the length of 6 there are two sentences of the same : "part of a series"
    set5 = Set()
    set6 = Set()
    set7 = Set()

    open(testfile) do file
        for line in eachline(file)
            len = length(split(line))
            (len == 3) && push!(set5, split(line)) #push!(lines5, split(line))
            (len == 4) && push!(set6, split(line))
            (len == 5) && push!(set7, split(line))
        end
    end
    return (set5, set6, set7)
end

""" Builds the kth sentence from a given sequence """
function ibuild_sentence(tdata::Data, sequence::Array{BitArray{2},1}, kth::Int)
    sentence = Array{Any, 1}()
    for i=1:length(sequence)
        z = find(x->x==true, sequence[i][kth,:])
        append!(sentence, z)
    end
    ret = map(x->tdata.index_to_word[x], sentence)
    return ret[2:end-1] # because first character and the last character are padded.
end

"""Compares the whether batched sentence in the given test or not, we expect it to return empty array"""
function ihelp_sentence_batch(tdata, batches, set; test_surplus=true)
    ret = Any[]
    for sequence in batches
        l = size(sequence[1])[1]
        sentences = map(x->ibuild_sentence(tdata, sequence, x), collect(1:l))
        (test_surplus) && (l < 3) && push!(ret, sentences)
        for item in sentences
            assert(item in set)
        end
    end
    return ret
end

function test_sentence_batch(testfile)
    # get ready for the sentence batch
    tdata = test_datatype(testfile)
    batchsize = 3
    test_cleanseqdict!(tdata.sequences, batchsize)
    vocabsize = length(tdata.word_to_index)
    seq5 = tdata.sequences[5]
    seq6 = tdata.sequences[6]
    seq7 = tdata.sequences[7]

    sets = icreate_sets(testfile)
    set5 = sets[1]
    set6 = sets[2]
    set7 = sets[3]

    batches5 = icreate_batches(seq5, batchsize, vocabsize)
    batches6 = icreate_batches(seq6, batchsize, vocabsize)
    batches7 = icreate_batches(seq7, batchsize, vocabsize)

    x1 = ihelp_sentence_batch(tdata, batches5, set5)
    x2 = ihelp_sentence_batch(tdata, batches6, set6)
    x3 = ihelp_sentence_batch(tdata, batches7, set7)
    return (x1, x2, x3)
end

function test_cleanseqdict!(seqdict, batchsize)
    clean_seqdict!(seqdict, batchsize)
end

function test_iteration(s::Data)
    sequences = Any[]
    vocabsize = length(s.word_to_index)
    for sequence in s
        push!(sequences, sequence)
    end

    # Those part of the code is very costly and use it just to make sure batching done correctly
    # for seq in sequences
    #     for item in seq
    #         @assert (size(item)==(s.batchsize, vocabsize)) "$(size(item))"
    #     end
    # end
    return sequences
end

function main()
    datafile = "ptb/ptb.train.txt"
    testfile = "testdata.txt"

    # tests the vocabulary creation correctedness
    # vocabfile = "test_vocab"
    # global V = vocab_from_file(vocabfile)
    # global x = length(V)
    # words = split(readstring(testfile))
    # global vocab = Dict()
    # for word in words; get!(vocab, word, 1+length(vocab));end
    # for k in keys(vocab); assert(k in V);end

    # test of Data type
    # global tdata = test_datatype(testfile)

    # tests the sequence property that is being hold by the Data type

    # test of sentenbatch
    # global surplused = test_sentence_batch(testfile)

    # test of surplus cleaning
    # test_cleanseqdict!(tdata.sequences, 3)

    # test of iterable implementation
    global tdata = test_datatype(testfile)

    # test the iteration
    global data = test_datatype(datafile)
    global sequences = test_iteration(data)

    # TODO: check that it <unk>s the unknown words if an external vocabulary is given
end
# main()

# Rough documentation:
# --
# Read the all sentence and put them into unique sets by icreate_sets
# After iterating over the Data type ibuild_sentence can build ith sentence by going over the whole sequence of sequence in Datatype


text:
export CopositiveInner

"""
    struct CopositiveInner{S} <: PolyJuMP.PolynomialSet
        # Inner approximation of the PSD cone, i.e. typically either
        # `SOSCone`, `DSOSCone` or `SDSOSCone`,
        psd_inner::S
    end

A symmetric matrix ``Q`` is copositive if ``x^{\\top} Q x \\ge 0`` for all
vector ``x`` in the nonnegative orthant. Checking copositivity is a
co-NP-complete problem [MK87] and this cone is only the inner approximation of
the cone of copositive symmetric matrices given by Minknowski sum of `psd_inner`
and the cone of symmetric matrices with nonnegative entries (the diagonal
entries can be chosen to be zero) [Lemma 3.164, BPT12].

The matrix with nonnegative entries can be interpreted as lagrangian
multipliers. For instance,
```julia
@polyvar x y
@constraint(model, x^2 - 2x*y + y^2 in CopositiveInner(SOSCone()))
```
is equivalent to
```julia
# Matrix that we require to be copositive
Q = [ 1 -1
     -1  1]
λ = @variable(model, lower_bound=0)
# Symmetric matrix of nonnegative entries
Λ = [0 λ
     λ 0]
using LinearAlgebra # For `Symmetric`
@constraint(model, Symmetric(Q - Λ) in PSDCone())
```
which is equivalent to
```julia
@polyvar x y
λ = @variable(model, lower_bound=0)
@constraint(model, x^2 - 2x*y + y^2 - 2*λ * x*y in SOSCone())
```
which is the same as, using the `domain` keyword,
```julia
@polyvar x y
@constraint(model, x^2 - 2x*y + y^2 in SOSCone(), domain = @set x*y ≥ 0)
```

For consistency with its equivalent forms, the [`GramMatrixAttribute`](@ref) for
this constraint is given by the gram matrix in the `psd_inner` cone, i.e. which
should be equal to `Q - Λ`.

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
*Semidefinite Optimization and Convex Algebraic Geometry*.
Society for Industrial and Applied Mathematics, **2012**.

[MK87] K. G. Murty and S. N. Kabadi.
*Some NP-complete problems in quadratic and nonlinear programming*.
Mathematical programming, 39:117–129, **1987**.
"""
struct CopositiveInner{S} <: SOSLikeCone
    psd_inner::S
end

function matrix_cone_type(::Type{CopositiveInner{S}}) where S
    return CopositiveInner{matrix_cone_type(S)}
end
function matrix_cone(set::Type{CopositiveInner{S}}, side_dimension) where S
    return CopositiveInner(matrix_cone(S, side_dimension))
end
side_dimension(set::CopositiveInner) = side_dimension(set.psd_inner)


text:
"""
    Perplexification{T <: Number} <: Plexification{T}

A perplexification as a pair.
"""
struct Perplexification{T <: Number} <: Plexification{T}
    l::T
    r::T
end

function (==)(x::Perplexification, y::Perplexification)
    x.l == y.l && x.r == y.r
end

function real(z::Perplexification)
    real(z.l)
end

function zero(::Type{Perplexification{T}}) where T <: Number
    Perplexification{T}(zero(T), zero(T))
end

function zero(::Perplexification{T}) where T <: Number
    zero(Perplexification{T})
end

function one(::Type{Perplexification{T}}) where T <: Number
    Perplexification{T}(one(T), zero(T))
end

function one(::Perplexification{T}) where T <: Number
    one(Perplexification{T})
end

function (-)(z::Perplexification)
    Perplexification(-(z.l), -(z.r))
end

function (+)(x::Perplexification{T}, y::Perplexification{T}) where T <: Number
    Perplexification{T}(x.l + y.l, x.r + y.r)
end

function (+)(a::T, z::Perplexification{T}) where T <: Number
    Perplexification{T}(a + z.l, z.r)
end

function (+)(z::Perplexification{T}, a::T) where T <: Number
    Perplexification{T}(z.l + a, z.r)
end

function (+)(a::Real, z::Perplexification)
    Perplexification(a + z.l, z.r)
end

function (+)(z::Perplexification, a::Real)
    Perplexification(z.l + a, z.r)
end

function (-)(x::Perplexification{T}, y::Perplexification{T}) where T <: Number
    Perplexification{T}(x.l - y.l, x.r - y.r)
end

function (-)(a::T, z::Perplexification{T}) where T <: Number
    Perplexification{T}(a - z.l, -(z.r))
end

function (-)(z::Perplexification{T}, a::T) where T <: Number
    Perplexification{T}(z.l - a, z.r)
end

function (-)(a::Real, z::Perplexification)
    Perplexification(a - z.l, -(z.r))
end

function (-)(z::Perplexification, a::Real)
    Perplexification(z.l - a, z.r)
end

function (*)(x::Perplexification{T}, y::Perplexification{T}) where T <: Number
    Perplexification{T}(
        (x.l * y.l) + (x.r * y.r),
        (x.l * y.r) + (x.r * y.l)
    )
end

function (*)(a::T, z::Perplexification{T}) where T <: Number
    Perplexification{T}(a * z.l, a * z.r)
end

function (*)(z::Perplexification{T}, a::T) where T <: Number
    Perplexification{T}(z.l * a , z.r * a)
end

function (*)(a::Real, z::Perplexification)
    Perplexification(a * z.l, a * z.r)
end

function (*)(z::Perplexification, a::Real)
    Perplexification(z.l * a, z.r * a)
end

function (/)(z::Perplexification, a::Real)
    if iszero(a)
        error(ZeroDenominator)
    end

    Perplexification(z.l / a, z.r / a)
end

function (\)(a::Real, z::Perplexification)
    if iszero(a)
        error(ZeroDenominator)
    end

    Perplexification(a \ z.l, a \ z.r)
end


text:
export triangulation

function get_angle_from_sides(a::Float64, b::Float64, c::Float64)
    return acosd((a ^ 2 + b ^ 2 - c ^ 2) / (2 * a * b))
end

function checkangle(a::Float64, b::Float64, c::Float64)
    angle = []
    try
        push!(angle, get_angle_from_sides(a, b, c))
        push!(angle, get_angle_from_sides(c, a, b))
        push!(angle, get_angle_from_sides(b, c, a))
        return sort(angle)
    catch
        return [Nan, Nan, Nan]
    end
end

function triangulation(InputDict::Dict)

    finame = InputDict["finame"]
    fogmtname = InputDict["fogmtname"]
    fodvvname = InputDict["fodvvname"]
    maxdistance = InputDict["maxdistance"]
    mindistance = InputDict["mindistance"]
    maxangle = InputDict["maxangle"]
    minangle = InputDict["minangle"]
    studyarea = InputDict["studyarea"]
    plotfig = InputDict["plotfig"]

    # for convert from LLA coordinate to ENU coordinate
    origin_lla = LLA((studyarea[1]+studyarea[3])/2, (studyarea[2]+studyarea[4])/2, 0.0)
    trans = ENUfromLLA(origin_lla, wgs84)

    t = jldopen(finame)["stationlist"]

    dist = []
    stationnum=length(t)

    # compute Delaunay triangulation of stations
    sta_lonlat = Array{Float64, 2}(undef, stationnum, 2)
    for i = 1:stationnum
        sta_lonlat[i, :] = Float64[t[i]["coord_lon"] t[i]["coord_lat"]]
    end

    println("compute Delaunay.")
    tri_temp = Triangle.basic_triangulation_vertices(sta_lonlat)

    traceall = GenericTrace{Dict{Symbol,Any}}[]

    tri = []
    dist = []

    for i = 1:length(tri_temp)
        # compute distances of all pairs
        p1 = [tri_temp[i][1 ,1], tri_temp[i][1,2]] #[lon, lat]
        p2 = [tri_temp[i][2 ,1], tri_temp[i][2,2]]
        p3 = [tri_temp[i][3 ,1], tri_temp[i][3,2]]

        pp1 = LLA(p1[2], p1[1]) #LLA(lat, lon, alt = 0.0)
        pp2 = LLA(p2[2], p2[1])
        pp3 = LLA(p3[2], p3[1])

        # pe1 = trans(pp1)
        # pe2 = trans(pp2)
        # pe3 = trans(pp3)
        #
        # v1 = diff(hcat([pe1[1], pe2[1]], [pe1[2], pe2[2]]), dims=1)
        # v2 = diff(hcat([pe2[1], pe3[1]], [pe2[2], pe3[2]]), dims=1)
        # v3 = diff(hcat([pe3[1], pe1[1]], [pe3[2], pe1[2]]), dims=1)
        #
        # dist1 = norm(v1)
        # dist2 = norm(v2)
        # dist3 = norm(v3)

        #check
        dist1= distance(pp1, pp2);
        dist2= distance(pp2, pp3);
        dist3= distance(pp3, pp1);
        #
        # if dist11 == dist1 && dist22 == dist2 && dist33 == dist3
        #     println("Pass Distance test")
        # else
        #     println([dist11, dist1, dist22, dist2, dist33, dist3])
        #     @warn "distance wrong"
        # end

        angle = checkangle(dist1, dist2, dist3)

        # println([dist1, dist2, dist3])
        # println(angle)
        # println( any(x -> x>maxdistance, [dist1, dist2, dist3]))
        # println( any(x -> x<mindistance, [dist1, dist2, dist3]))
        # println( any(x -> x<minangle, angle))
        # println( any(x -> x>maxangle, angle))

        if any(x -> x>maxdistance, [dist1, dist2, dist3]) ||
            any(x -> x<mindistance, [dist1, dist2, dist3]) ||
            any(x -> x<minangle, angle) ||
            any(x -> x>maxangle, angle)
            continue
        else

            push!(tri, tri_temp[i])
            push!(dist, [dist1, dist2, dist3])

            trace1 = scatter(;x=tri_temp[i][[1,2],1], y=tri_temp[i][[1,2],2], mode="lines+markers")
            trace2 = scatter(;x=tri_temp[i][[2,3],1], y=tri_temp[i][[2,3],2], mode="lines+markers")
            trace3 = scatter(;x=tri_temp[i][[3,1],1], y=tri_temp[i][[3,1],2], mode="lines+markers")

            push!(traceall, trace1)
            push!(traceall, trace2)
            push!(traceall, trace3)
        end
    end

    if plotfig
        fodir = "./fig"
        if !ispath(fodir) mkpath(fodir); end

        mul = 400
        height= abs(diff([studyarea[1], studyarea[3]])[1]) * mul
        width = abs(diff([studyarea[2], studyarea[4]])[1]) * mul
        layout = Layout(;width=width, height=height, title="Triangulation of station network",
                       xaxis_range=[studyarea[4], studyarea[2]],
                       yaxis_range=[studyarea[3], studyarea[1]],
                       xaxis_title="lon",
                       yaxis_title="lat",
                       showlegend=false)

        p = PlotlyJS.plot(traceall, layout)
        PlotlyJS.savefig(p, fodir*"/"*InputDict["figname"])
        #display(p)
        #println("type return...")
        #readline()
    end

    # compute max, min, average distance
    maxedge= maximum(maximum(dist))
    minedge= minimum(minimum(dist))
    meanedge = mean(mean(dist))
    # paring with station id
    triall = []
    for i = 1:length(tri)
        edgeinfo = []
        for (j, k) in [(1,2), (2,3), (3,1)]
            edge = [tri[i][j, 1] tri[i][j, 2]; tri[i][k, 1] tri[i][k, 2]]
            sta1 = findfirst(x -> x == edge[1], sta_lonlat)
            sta2 = findfirst(x -> x == edge[2], sta_lonlat)
            sta1id = t[sta1]["stationname"]
            sta2id = t[sta2]["stationname"]
            staloc1 =  edge[1, 1:2]
            staloc2 =  edge[2, 1:2]
            midpoint = mean(edge, dims=1)
            push!(edgeinfo, (sta1id, sta2id, staloc1, staloc2, midpoint))
        end
        centerxy= mean(tri[i], dims=1)
        push!(triall, Dict("edgeinfo" => edgeinfo, "centerxy" => centerxy))
    end

    #---output for dv/v---#
    jldopen(fodvvname, "w") do file
        file["tri"] = triall
        file["maxedgesize"] = maxedge
        file["minedgesize"] = minedge
        file["meanedgesize"] = meanedge
    end

    #---output for gmt---#
    jldopen(fogmtname, "w") do file
        edges = []
        for i = 1:length(tri)
            for (j, k) in [(1,2), (2,3), (3,1)]
                edge = [tri[i][j, 1] tri[i][j, 2]; tri[i][k, 1] tri[i][k, 2]]
                if !any(x -> x == edge, edges)
                    #to avoid overlap of line
                    push!(edges, edge)
                end
            end
        end
        file["edges"] = edges
    end

    println("analysis has been successfully done.")
    return nothing

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

# Measurements

"""
    angle = planarRotationAngle(e, v1, v2)
          = planarRotationAngle(frame1, frame2)

Return angle of a planar rotation, given the rotation axis e (a unit vector)
and the representations of a vector in frame 1 (v1) and frame 2 (v2).

Under the assumption that the z-axes of frame1 and frame2 coincide, return the
angle between the x-axis of frame1 and the position vector from frame1 to frame2.
"""
planarRotationAngle(e::AbstractVector, v1::AbstractVector, v2::AbstractVector) =
     atan2( dot(-cross(e,v1), v2), dot(v1,v2) - dot(e,v1)*dot(e,v2) )

function planarRotationAngle(frame1::Object3D, frame2::Object3D)
   # Extract needed vectors
   e1z = frame1.R_abs[3,:]
   e2z = frame2.R_abs[3,:]
   e1x = frame1.R_abs[1,:]
   r1  = frame1.r_abs
   r2  = frame2.r_abs

   # Check that z-axes of both frames coincide
   angle_z_axes = acos( dot(e1z,e2z) )
   if abs(angle_z_axes) > Basics.neps
      error("\nError from Modia3D.planarRotationAngle(", ModiaMath.fullName(frame1), ", ", ModiaMath.fullName(frame2), "):\n",
                 "The z-axes of the two frames do not coincide.\n",
                 "( angle(zaxis(frame1), zaxis(frame2)) = ", string(angle_z_axes*Basics.radToDeg), " > ", Basics.neps, ")")
   end

   # Angle between frame2.r_abs - frame1.r_abs and x_axis of frame1
   r12 = r2-r1
   if norm(r12) < 100*eps()
      error("\nError from Modia3D.planarRotationAngle(", ModiaMath.fullName(frame1), ", ", ModiaMath.fullName(frame2), "):\n",
              "The origins of the two frames coincide (this is not allowed for this measurement.")
   end
   r12n = r12/norm(r12)

   return atan(dot(cross(e1x,r12n), e1z), dot(e1x,r12n) )
end

"""
    d = distance(frame1, frame2)

Return the distance between the origin of frame1 and the origin of frame2
"""
distance(frame1::Object3D, frame2::Object3D) = norm(frame2.r_abs - frame1.r_abs)

function distanceAndAngles(frame1::Composition.Object3D, frame2::Composition.Object3D)
   # Extract needed vectors
   e1z = frame1.R_abs[3,:]
   e2z = frame2.R_abs[3,:]
   r1  = frame1.r_abs
   r2  = frame2.r_abs

   # Check that z-axes of both frames are pointing in the same direction
   @assert(abs(dot(e1z,e2z) - 1.0) <= 1e-10)

   # Angle between frame1 and frame2 along z-axis
   rdiff0 = r2-r1
   ey = normalize(rdiff0)

   ey1  = frame1.R_abs*ey
   ey2  = frame2.R_abs*ey
   ex1b = cross(ey1, e1z)
   ex2b = cross(ey2, e2z)
   ex1a = frame1.R_abs[:,1]
   ex2a = frame2.R_abs[:,1]

   phi1 = planarRotationAngle(e1z, ex1a, ex1b)
   phi2 = planarRotationAngle(e2z, ex2a, ex2b)
   d    = norm(rdiff0)
   return (d,phi1,phi2)
end


text:
using LightGraphs, GraphMakie, NetworkLayout, GLMakie

function to_graph(is::InventorySystem)
    g = SimpleDiGraph(0)
    dict = IdDict()
    dict2 = Dict()
    for (i,v) in enumerate(is.bom)
        add_vertex!(g)
        dict[v] = i
        dict2[i] = v
    end
    for v in is.bom
        for c in InventoryModels.children(v)
            add_edge!(g, dict[c], dict[v])
        end
    end
    return g,dict2
end

struct ISLayout
end

function (l::ISLayout)(g)
    layout(l, g)
end

function layout(::ISLayout, g)
    rows = Vector{Int}[]
    laidnodes = Set{Int}()
    while length(laidnodes) < maximum(vertices(g))
        toprow = Int[]
        push!(rows, toprow)
        for v in vertices(g)
            if (isempty(g.fadjlist[v]) || (all(f-> (f in laidnodes), g.fadjlist[v]) && all(f -> !(f in toprow), g.fadjlist[v]))) && !(v in laidnodes)
                push!(toprow, v)
                push!(laidnodes, v)
            end
        end
    end
    nrows = length(rows)
    pos = Vector{Tuple{Float64,Float64}}(undef, nv(g))
    largest = maximum(length.(rows))
    center = iseven(largest) ? largest/2 : ceil(largest/2)
    for (i,row) in enumerate(rows)
        for (j,v) in enumerate(row)
            x = length(row) == 1 ? center : largest/length(row) + j - 1
            pos[v] = (x, nrows - i)
        end
    end
    pos
end

function draw_graph(is::InventorySystem)
    g,dict = to_graph(is)
    layout = ISLayout()
    f,a,p = graphplot(g, layout = layout,nlabels=[dict[i].name for i in 1:nv(g)], nodesize = 15, arrow_shift = 0.9, )
    hidespines!(a); hidedecorations!(a)
    display(f)
    f,a,p
end


text:
using Revise,LinearAlgebra, Random, StatsBase, Statistics, Distributions

Random.seed!(1111)

printm(m) = show(IOContext(stdout, :limit=>false), MIME"text/plain"(), m)
function printmln(m)
  show(IOContext(stdout, :limit=>false), MIME"text/plain"(), m)
  println()
end

function testex2(;
  Ξ = ([4,8], [4,-8], [8,4], [-4, -8], [-4, 8], [-8, -4]),
  dims=(N=6, T=10^5, K=2))

  #set up the factor structure and dividend realizations
  Θ0 = reduce(hcat, [rand(Normal(), dims.T) for k in 1:dims.K])
  @assert size(Θ0) ≡ (dims.T, dims.K)
  D̄ = 1.0
  ϵ = [rand(Normal(), dims.T) for n ∈ 1:dims.N]
  Dvec = [D̄ .+ Θ0*Ξ[n] .+ ϵ[n] for n ∈ 1:dims.N]
  D = reduce(hcat, Dvec)
  Σ = cov(D)

  #select a baseline and a bunch of comp portfolios
  w1 = [1/6 for n ∈ 1:6]
  walt = [rand(6) for i ∈ 1:10^5]
  walt ./= sum.(walt)

  varbase = w1'*Σ*w1
  @info "baseline: w1 $w1: $(w1'*Σ*w1)"

  vars = (w-> w'*Σ*w).(walt)
  for (i,w) ∈ enumerate(walt)
    if vars[i] > varbase
      #@info "w$i $w: $(vars[i])"
      nothing
    else
      @warn "w$i $w: $(vars[i])"
    end
  end

  printmln(Σ)
end

testex2()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Documenter, InitialValues

makedocs(;
    modules=[InitialValues],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaFolds/InitialValues.jl/blob/{commit}{path}#L{line}",
    sitename="InitialValues.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    strict = true,
)

deploydocs(;
    repo="github.com/JuliaFolds/InitialValues.jl",
    push_preview = true,
)


text:
########################################################################################################
## Details
# This file will create the structure for the UQ inference problems.
# This package is designed to be used with a family of julia packages,
#           - Potential.jl (*)
#               - Creates potential types, has abstract 'Potential' type (*), let p <: 'Potential'
#               - Assume mutable struct Potential{T<:AbstractFloat}
#                            x       # Potential trainable parameters
#                            p       # Potential nontrainable parameters
#                        end
#               - Assume mutable struct SNAP{T<:AbstractFloat}
#                            A
#                            b
#                            β
#                        end
#           - PotentialLearning.jl
#               - Fits data-driven potential
#               - We assume that potential learning updates parameters p.x
#           - PotentialUQ.jl (this package)
#               - Defines distribution over potential parameters.
#               - Initially must provide distribution over p.x (assume user provides distribution, in form of log likelihood)
#               - In case of SNAP, potential is predetermined.
#               - Allow users to modify distribution with priors by adding to log-likelihood.
#
#########################################################################################################
module PotentialUQ

using Base: NamedTuple
using LinearAlgebra: tr, logdet, transpose, Symmetric, diag, I, pinv
using Statistics: cov, std, mean
using SciMLBase, GalacticOptim, Optim
using TransformVariables
import Turing
using Distributions

include("distribution_types.jl")
include("distributions.jl")
include("utils.jl")
include("map.jl")
include("sampling.jl")
end


text:
module BlinkPOMDPSimulator

using POMDPs
using POMDPModelTools
using POMDPSimulators
using POMDPPolicies
using Parameters
using Blink
using Random

export BlinkSimulator

struct TestWindow end

include("testing.jl")

const WindowLike = Union{Window, TestWindow}

@with_kw mutable struct BlinkSimulator <: Simulator
    max_steps::Int         = typemax(Int)
    max_fps::Float64       = 2.0
    extra_initial::Bool    = false
    extra_final::Bool      = false
    window::WindowLike     = Window()
    render_kwargs          = NamedTuple()
    rng::AbstractRNG       = Random.GLOBAL_RNG
end

function POMDPs.simulate(sim::BlinkSimulator,
                         m::MDP,
                         p::Policy=RandomPolicy(m, rng=sim.rng),
                         is=initialstate(m, sim.rng)
                        )

    dt = 1/sim.max_fps
    tm = time()

    if sim.extra_initial
        step = (t=0, sp=is)
        gfx = render(m, step; sim.render_kwargs...)
        blink!(sim.window, gfx)
        sleep_until(tm += dt)
    end

    nsteps = 0
    last_sp = missing
    for step in stepthrough(m, p, is; max_steps=sim.max_steps, rng=sim.rng)
        gfx = render(m, step; sim.render_kwargs...)
        blink!(sim.window, gfx)
        sleep_until(tm + dt)
        tm = max(time(), tm + dt)
        nsteps += 1
        last_sp = get(step, :sp, missing)
    end

    if sim.extra_final
        step = (t=nsteps+1, s=last_sp, done=true)
        gfx = render(m, step; sim.render_kwargs...)
        blink!(sim.window, gfx)
    end
end

function POMDPs.simulate(sim::BlinkSimulator,
                         m::POMDP,
                         p::Policy=RandomPolicy(m, rng=sim.rng),
                         u::Updater=updater(p),
                         isd=initialstate_distribution(m),
                         is=initialstate(m, sim.rng)
                        )

    b0 = initialize_belief(u, isd)

    dt = 1/sim.max_fps
    tm = time()

    if sim.extra_initial
        step = (t=0, sp=is, bp=b0)
        gfx = render(m, step; sim.render_kwargs...)
        blink!(sim.window, gfx)
        sleep_until(tm += dt)
    end

    nsteps = 0
    last = (s=missing, b=missing)
    for step in stepthrough(m, p, u, b0, is; max_steps=sim.max_steps, rng=sim.rng)
        gfx = render(m, step; sim.render_kwargs...)
        blink!(sim.window, gfx)
        sleep_until(tm + dt)
        tm = max(time(), tm + dt)
        nsteps += 1
        last = (s=get(step, :sp, missing), b=get(step, :bp, missing))
    end

    if sim.extra_final
        step = merge((t=nsteps+1, done=true), last)
        gfx = render(m, step; sim.render_kwargs...)
        blink!(sim.window, gfx)
    end
end

sleep_until(t) = sleep(max(t-time(), 0.0))

blink!(win::Window, rendering) = body!(win, rendering, fade=false)

end # module


text:
import Polyhedra.polyhedron
# Fullin also provides a certificate so it is interesting so it is exported
export getextremerays, fullin

function unlift(h::EntropyConeLift)
    poly = getpoly(h)
    EntropyCone(eliminate(poly, [(ntodim(h.n[1])+1):Polyhedra.fulldim(h)]))
end

function fullin(h::AbstractPrimalEntropy, H::AbstractEntropyCone)
    Ray(h.h) in H.poly
    #all(H.A*h.h .>= 0)
end
function partialin(h::AbstractPrimalEntropy{S}, H::AbstractEntropyCone{T}) where {S, T}
    hps = HyperPlane{T, SparseVector{T, Int}}[]
    offseth = 0
    offsetsH = [0; cumsum(map(ntodim, H.n))]
    for i in eachindex(collect(h.n)) # use of collect in case h.n is scalar
        for j in indset(h, i)
            col = offsetsH[h.liftid[i]]+j
            push!(hps, HyperPlane(sparsevec([col], [one(T)], Polyhedra.fulldim(H)), T(h.h[offseth+j])))
        end
        offseth += ntodim(h.n[i])
    end
    !isempty(H.poly ∩ hrep(hps))
end

function Base.in(h::PrimalEntropy, H::EntropyCone)
    if Polyhedra.fulldim(h) > Polyhedra.fulldim(H)
        error("The vector has a higher dimension than the cone")
    elseif Polyhedra.fulldim(h) == Polyhedra.fulldim(H)
        fullin(h, H)[1]
    else
        partialin(h, H)
    end
end

function Base.in(h::PrimalEntropyLift, H::EntropyConeLift)
    if length(h.n) > length(H.n) || any(h.n .> H.n)
        error("The vector has a higher dimension than the cone")
    elseif h.n == H.n
        fullin(h, H)[1]
    else
        partialin(h, H)
    end
end

function Base.in(h::PrimalEntropy, H::EntropyConeLift)
    if h.liftid < 1 || h.liftid > length(H.n) || h.n > H.n[h.liftid]
        error("The vector has a higher dimension than the cone")
    elseif h.n == H.n
        fullin(h, H)[1]
    else
        partialin(h, H)
    end
end

#function redundant(h::AbstractDualEntropy{L, S}, H::AbstractEntropyCone{T}) where {L, S, T}
#    (isin, certificate, vertex) = ishredundant(H.poly, HRepElement(h))
#    (isin, certificate, vertex)
#end

Base.in(h::DualEntropy, H::EntropyCone) = H.poly ⊆ HRepElement(h)

Base.in(h::DualEntropyLift, H::EntropyConeLift) = H.poly ⊆ HRepElement(h)

function Base.in(h::DualEntropy, H::EntropyConeLift)
    Base.in(DualEntropyLift(h, length(H.n)), H)
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Functions for map coordinate transformations ###

###############################################
### Conversion from LLA to ECEF coordinates ###
###############################################

# For dictionary of nodes
function Geodesy.ECEF(nodes::Dict{Int,LLA}, datum::Ellipsoid = WGS84)
    r = Dict{Int,ECEF}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = ECEF(node, datum)
    end

    return r
end

###############################################
### Conversion from ECEF to LLA coordinates ###
###############################################

# For dictionary of nodes
function Geodesy.LLA(nodes::Dict{Int,ECEF}, datum::Ellipsoid = WGS84)
    r = Dict{Int,LLA}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = LLA(node, datum)
    end

    return r
end

#####################################
### Conversion to ENU coordinates ###
#####################################

# Given a reference point
function Geodesy.ENU{T<:@compat Union{LLA,ECEF}}(nodes::Dict{Int,T},
                                         lla_ref::LLA,
                                         datum::Ellipsoid = WGS84)
    r = Dict{Int,ENU}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = ENU(node, lla_ref, datum)
    end

    return r
end

# Given Bounds
function Geodesy.ENU(nodes::Dict, bounds::Bounds{LLA}, datum::Ellipsoid = WGS84)
    ENU(nodes, center(bounds), datum)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c7952ee6-45b5-11eb-1158-5bb2ff274ce9
using DrWatson

# ╔═╡ 0c145632-3927-11eb-19b9-877e05c1bcdc
begin
	import Pkg
	#Pkg.activate("Mage")

	using GameOne
	using JSON
	using HTTP
	using Plots
	using Images: load, ARGB
	using PlutoUI
	using Serialization
	using PlaymatSimulator
	using ImageTransformations: imresize

	AC = PlaymatSimulator.Actors

	plotly()

	md"""
	*game settings loaded!*

	## MtG EDH Deck: Vannifar's Circus (UG Creature Ramp Combo) by Dustin Irwin
	For simple games, a "deck" in PlaymatSimulator is simply a collection of images, one for each card in the deck. Let's build an example EDH deck `Vannifar's Circus`.

	To get started, define a `deck` Dict object below of type Dict{String,Int} where the key is the official card name and the value is the quantity of that card in the deck.
	"""
end

# ╔═╡ 621b08a4-384e-11eb-0109-61e9b9ecf125
deck = Dict{Symbol,Any}(
	:name => split(@__DIR__, "/")[end],
	:commander_names => [
	    "Kynaios and Tiro of Meletis",
	    ],
	:card_names => [
	    "Altar of the Pantheon",
	    "Angel of Sanctions",
	    "Armillary Sphere",
	    "Ash Barrens",
	    "Azorius Chancery",
	    "Back from the Brink",
	    "Boros Guildgate",
	    "Brudiclad, Telchor Engineer",
	    "Caller of the Pack",
	    "Commander's Sphere",
	    "Command Tower",
	    "Cultivate",
	    "Desolation Twin",
	    "Doomed Artisan",
	    "Draconic Disciple",
	    "Dragonmaster Outcast",
	    "Druid's Deliverance",
	    "Emmara Tandris",
	    "Ephara, God of the Polis",
	    "Evolving Wilds",
	    "Exotic Orchard",
	    "Feldon of the Third Path",
	    "Forbidden Orchard",
	    "Forest",
	    "Forest",
	    "Forest",
	    "Forest",
	    "Forest",
	    "Forest",
	    "Full Flowering",
	    "Gargoyle Castle",
	    "Garruk, Primal Hunter",
	    "Ghired, Conclave Exile",
	    "Ghired's Belligerence",
	    "Giant Adephage",
	    "God-Pharaoh's Gift",
	    "Golden Guardian // Gold-Forge Garrison",
	    "Growing Ranks",
	    "Gruul Turf",
	    "Heart-Piercer Manticore",
	    "Hellion Crucible",
	    "Helm of the Host",
	    "Hour of Reckoning",
	    "Idol of Oblivion",
	    "Intangible Virtue",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Island",
	    "Izzet Boilerworks",
	    "Kazandu Tuskcaller",
	    "Kiora, the Crashing Wave",
	    "Metallurgic Summonings",
	    "Mimic Vat",
	    "Mirror Match",
	    "Mist-Syndicate Naga",
	    "Moonsilver Spear",
	    "Mountain",
	    "Mountain",
	    "Mountain",
	    "Mountain",
	    "Mountain",
	    "Mountain",
	    "Myriad Landscape",
	    "Overwhelming Stampede",
	    "Parhelion II",
	    "Phyrexian Rebirth",
	    "Rogue's Passage",
	    "Rampaging Baloths",
	    "Rootborn Defenses",
	    "Saheeli Rai",
	    "Saheeli's Artistry",
	    "Second Harvest",
	    "Selesnya Eulogist",
	    "Selesnya Sanctuary",
	    "Simic Growth Chamber",
	    "Song of the Worldsoul",
	    "Soul Foundry",
	    "Soul of Eternity",
	    "Spawning Grounds",
	    "Spectral Searchlight",
	    "Spitting Image",
	    "Stolen Identity",
	    "Sundering Growth",
	    "Tempt with Discovery",
	    "Terramorphic Expanse",
	    "Titan Forge",
	    "Trostani, Selesnya's Voice",
	    "Trostani's Judgment",
	    "Vitu-Ghazi Guildmage",
	    "Wayfaring Temple",
	    "Wingmate Roc",
    ]
)

# ╔═╡ fb61d01c-458d-11eb-2c2a-f711dc7ab7f4
(length(deck[:card_names]) + length(deck[:commander_names]))

# ╔═╡ c61fa79e-4583-11eb-3b71-2d334aca843d
begin
	mtg_cards = JSON.parsefile("$(projectdir())/games/MtG/MtG.jl/json/oracle-cards-20201224220555.json")

	md"""Look OK? Keep in mind that images that do not require in-game scaling will suffer less distortion.

	Alright, lets load up a JSON file with the URIs we need to grab the card images. For MtG, we can use the .json file available here: TODO

	Save the json file to the /json directory in the MtG project directory and modify the following cell to point at the json file.

	##### MtG database loaded! Found $(length(mtg_cards)) unique cards (by name).

	mtg_cards is of type Array{Any}. The dicts contained within are of type Dict{String,Any}.

	Alrighty, let's collect the data we need to download the card images!
	"""
end

# ╔═╡ 7420cf10-45cc-11eb-2780-4f320bd8a2cf
begin  # note: this func only downloads the first card with a matching name and then moves to the next.
	deck_cards = []
	commander_cards = []

	for n in sort(deck[:card_names])

		for c in mtg_cards

			if n == c["name"]
				push!(deck_cards, c)
			end
		end
	end

	for n in deck[:commander_names]

		for c in mtg_cards

			if n == c["name"]
				push!(commander_cards, c)
			end
		end
	end

	all_cards = vcat(commander_cards, deck_cards)

	md"""
	Found $(length(deck_cards) + length(commander_cards)) matching cards in mtg_cards!
	"""
end

# ╔═╡ c31dd202-50c6-11eb-0631-13c70535635e
missing_cards = filter!(x->!(x in [ c["name"] for c in all_cards ]), vcat(deck[:card_names], deck[:commander_names]))

# ╔═╡ 2775088a-4648-11eb-2218-af69e0e95f1f
@bind i Slider(1:length(all_cards), show_value=true)

# ╔═╡ c5374766-4ef1-11eb-2555-c159dba953f0
md"""
##### *Adjust this slider to shrink / grow the cards while preserving the aspect ratio*

$(@bind card_ratio Slider(0.1:0.05:1.25, default=0.5, show_value=true))
"""

# ╔═╡ 5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
begin
	CARD_BACK_PATH = "$(projectdir())/games/MtG/MtG.jl/ui/cards/card_back.png"
	CARD_BACK_IMG = imresize(load(CARD_BACK_PATH), ratio=card_ratio)
end

# ╔═╡ 614764b8-4648-11eb-0493-732a00df7bca
md"""
#### Look good? These images will be displayed in-game!
TODO: write support for previews of double-sided cards
"""

# ╔═╡ dfc9b56e-50ce-11eb-0e7f-83ec6e831901
begin
	#=
	CARD_FRONT_IMGS = []
	COMMANDER_FRONT_IMGS = []

	for c in deck_cards
		push!(CARD_FRONT_IMGS, imresize(get_mtg_card_front_img(c), ratio=card_ratio))
		sleep(0.1)
	end

	for c in commander_cards
		push!(COMMANDER_FRONT_IMGS, imresize(get_mtg_card_front_img(c), ratio=card_ratio+0.05))
		sleep(0.1)
	end

	COMMANDER_FRONT_IMGS, CARD_FRONT_IMGS
	=#
end

# ╔═╡ ce216c54-468a-11eb-13b8-7f3dac7af44a
function get_card_img(img_uri::String)
	img_resp = HTTP.get(img_uri)
	card_img = img_resp.body |> IOBuffer |> load
end

# ╔═╡ 2ab53d00-50cd-11eb-1cd4-5bf94ce53692
function get_mtg_card_img(c)
	if haskey(all_cards[i], "card_faces") && haskey(all_cards[i]["card_faces"][1], "image_uris")
		hcat([
			get_card_img(f["image_uris"]["border_crop"]) for f in all_cards[i]["card_faces"]
			]...)
	else
		get_card_img(all_cards[i]["image_uris"]["border_crop"])
	end
end

# ╔═╡ 73d1cd18-4647-11eb-3994-7d4eb92eddca
if (@isdefined all_cards) && length(all_cards) > 0
	deck_card_preview = imresize(get_mtg_card_img(all_cards[i]), ratio=card_ratio)
else
	nothing
end

# ╔═╡ 97bc3768-50ce-11eb-3f74-95d4fefe3792
function get_mtg_card_front_img(c)
	if haskey(c, "card_faces") && haskey(c["card_faces"][1], "image_uris")
		get_card_img(c["card_faces"][1]["image_uris"]["border_crop"])
	else
		get_card_img(c["image_uris"]["border_crop"])
	end
end

# ╔═╡ d654ad1e-468a-11eb-2348-695621b7b9b0
function search_mtg_cards_by_keyword(q::String, mtg_cards::Array)
	[ n for n in [ c["name"] for c in mtg_cards ] if occursin(q, n) ]
end

# ╔═╡ cec924ac-50c7-11eb-3795-85b3c183a8eb
search_mtg_cards_by_keyword("Meletis", mtg_cards)

# ╔═╡ 5e6f7046-4da5-11eb-0122-bd82397aab4f
begin
	#=
	deck[:CARD_WIDTH] = Int32(size(CARD_BACK_IMG)[1])
	deck[:CARD_HEIGHT] = Int32(size(CARD_BACK_IMG)[2])
	deck[:CARD_FRONT_IMGS] = CARD_FRONT_IMGS
	deck[:CARD_BACK_IMG] = CARD_BACK_IMG
	deck[:COMMANDER_FRONT_IMGS] = COMMANDER_FRONT_IMGS

	fn = "$(projectdir())/games/MtG/EDH/decks/$(deck[:name])/$(deck[:name]).jls"
	serialize(fn, deck)
	=#
end

# ╔═╡ Cell order:
# ╟─c7952ee6-45b5-11eb-1158-5bb2ff274ce9
# ╟─0c145632-3927-11eb-19b9-877e05c1bcdc
# ╟─621b08a4-384e-11eb-0109-61e9b9ecf125
# ╠═fb61d01c-458d-11eb-2c2a-f711dc7ab7f4
# ╟─5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
# ╟─c61fa79e-4583-11eb-3b71-2d334aca843d
# ╟─7420cf10-45cc-11eb-2780-4f320bd8a2cf
# ╟─c31dd202-50c6-11eb-0631-13c70535635e
# ╠═cec924ac-50c7-11eb-3795-85b3c183a8eb
# ╟─73d1cd18-4647-11eb-3994-7d4eb92eddca
# ╟─2775088a-4648-11eb-2218-af69e0e95f1f
# ╠═c5374766-4ef1-11eb-2555-c159dba953f0
# ╟─614764b8-4648-11eb-0493-732a00df7bca
# ╠═dfc9b56e-50ce-11eb-0e7f-83ec6e831901
# ╟─ce216c54-468a-11eb-13b8-7f3dac7af44a
# ╟─2ab53d00-50cd-11eb-1cd4-5bf94ce53692
# ╟─97bc3768-50ce-11eb-3f74-95d4fefe3792
# ╟─d654ad1e-468a-11eb-2348-695621b7b9b0
# ╠═5e6f7046-4da5-11eb-0122-bd82397aab4f

max_stars_repo_path: EDH/decks/Kynaios and Tiro of Meletis/Kynaios and Tiro of Meletis.jl


text:
# TODO: KEEP
# function compareAll(Al::T1, Bl::T2; show::Bool=true, skip::Vector{Symbol}=Symbol[])::Bool where {T1 <: Union{SingleThreaded, MultiThreaded}, T2 <: Union{SingleThreaded, MultiThreaded}}
#   return T1 == T2
# end
import DistributedFactorGraphs: compare, compareAllSpecial

# These are now moved to DistributedFactorGraphs, with the exceptions of
# the functions with IIF-specific parameters.
# To extend these, import the relevant DFG compareX function and overload it.

function compare(p1::BallTreeDensity, p2::BallTreeDensity)::Bool
  return compareAll(p1.bt,p2.bt, skip=[:calcStatsHandle; :data]) &&
         compareAll(p1,p2, skip=[:calcStatsHandle; :bt])
end

function compareAllSpecial(A::T1, B::T2;
                    skip=Symbol[], show::Bool=true) where {T1 <: CommonConvWrapper, T2 <: CommonConvWrapper}
  #
  T1 != T2 && return false
  return compareAll(A, B, skip=skip, show=show)
end


text:
# this script runs the Hooghoudt model benchmark
using Athesis

using TimerOutputs

include("../benchmark/benchmark_utils.jl")

to = TimerOutput()

disabledTimer = TimerOutput()
disable_timer!(disabledTimer)

archs = ["CPU"]
@withCUDA archs = ["CPU", "GPU"]
myFloat = Float32

dims = [(16, 16, 16), (32, 32, 16), (64, 64, 16), (128, 128, 16), (256, 256, 16)] #, (512, 512, 16)]

for arch in archs
    @synctimeit to arch begin
        for dim in dims
            input = getDefaultInput(myFloat, "../Athesis.toml")
            input.nx = dim[1]
            input.ny = dim[2]
            input.nz = dim[3]

            # enforce precompilation
            @info "starting precompilation..."
            preSim = initSimulation(input, arch=="GPU", myFloat, disabledTimer)
            doExplicitIter!(1, preSim)

            # run timing
            label = string(dim)
            @info "run benchmark: " * label
            @synctimeit to label begin
                sim = initSimulation(input, arch == "GPU", myFloat, disabledTimer)

                nmax = 200
                for i in 1:nmax
                    hasCvg, = doExplicitIter!(i, sim, disabledTimer)
                    if hasCvg
                        break
                    end
                end

                @test true
            end
        end
    end
end

print_timer(to)

modelName = "Hooghoudt"
plotTimerOutput(archs, dims, to, modelName)


text:
__precompile__(true)
module PandasLite

using PyCall
using Compat
using TableTraits
using Dates
using Requires

export values, DataFrame, Series, np, pd
export @pytype, @pyasvec

const np = PyNULL()
const pd = PyNULL()

function __init__()
    copy!(np, pyimport_conda("numpy", "numpy"))
    copy!(pd, pyimport_conda("pandas", "pandas"))
    empty!(type_map)  # for behaving nicely in system image
    for (pandas_expr, julia_type) in pre_type_map
        type_map[pandas_expr()] = julia_type
    end
    for (pytype, jltype) in type_map
        PyCall.pytype_mapping(pytype, jltype)
    end
    PyCall.pytype_mapping(np.float32, Float32)
    PyCall.pytype_mapping(np.float64, Float64)
    PyCall.pytype_mapping(np.int32, Int32)
    PyCall.pytype_mapping(np.int64, Int64)
    if get(ENV, "PD_NO_CONSOLID", "0") == "1"
        noconsolidation()
    end
    @require Atom="c52e3926-4ff0-5f6e-af25-54175e0327b1" include("juno.jl")
end

version() = VersionNumber(pd.__version__)

const pre_type_map = []

# Maps a python object corresponding to a PandasLite class to a Julia type which
# wraps that class.
const type_map = Dict{PyObject, Type}()

abstract type PandasWrapped end

PyCall.PyObject(x::PandasWrapped) = x.pyo

macro pytype(name, class)
    quote
        mutable struct $(name) <: PandasWrapped
            pyo::PyObject
            $(esc(name))(pyo::PyObject) = new(pyo)
            $(esc(name))(args...; kwargs...) = ($class)()(args...; kwargs...)
        end

        # This won't work until PyCall is updated to support
        # the Julia 1.0 iteration protocol.
        function Base.iterate(x::$name, state...)
            res = Base.iterate(x.pyo, state...)
            if res === nothing
                return nothing
            else
                value, state = res
                return value, state
            end
        end

        Base.convert(::Type{$name}, o::PyObject) = $name(o)

        push!(pre_type_map, ($class, $name))
    end
end

function Base.Array(x::PyObject)
    if x.dtype.kind == "M"
        map(z -> unix2datetime(z / 1e9), x.view("int64"))
    elseif x.dtype.kind == "m"
        map(z -> Millisecond(z / 1e6), x.view("int64"))
    elseif x.dtype.kind == "O" && get(x, 0) isa String
        convert(Array{String}, x)
    elseif x.dtype.kind == "O"
        map(PyAny, PyAny(x))
    else
        convert(PyAny, x)
    end
end

Base.Array(x::PandasWrapped) = Array(x."values")

function Base.values(x::PandasWrapped)
    # Zero-copy conversion to a Julia native type is possible
    x_kind = x.dtype.kind
    if x_kind in ["i", "u", "f", "b"]
        pyarray = convert(PyArray, x."values")
        if pyarray.f_contig
            dims, T, ptr = size(pyarray), eltype(pyarray), pyarray.data
            if Int(ptr) % Base.datatype_alignment(T) == 0
                unsafe_wrap(Array, ptr, dims)
            else
                Aflat = unsafe_wrap(Array, Ptr{UInt8}(ptr), prod(dims) * sizeof(T))
                A = reshape(reinterpret(T, Aflat), dims)
            end
        else
            Array(pyarray)
        end
    else
        Array(x)
    end
end

struct StringRange{T <: AbstractString}
    start::T
    stop::T
end
(::Colon)(start::T, stop::T) where T <: AbstractString = StringRange{T}(start, stop)

Base.show(io::IO, r::StringRange) =
    print(io, '"', r.start, '"', ':', '"', r.stop, '"')

PyCall.PyObject(x::Colon) = pybuiltin("slice")(nothing, nothing, nothing)
PyCall.PyObject(x::StringRange) = pybuiltin("slice")(x.start, x.stop)

fix_arg(x) = x
fix_arg(x::StepRange) = pybuiltin("slice")(x.start, x.start + length(x) * x.step, x.step)
fix_arg(x::UnitRange) = fix_arg(StepRange(x.start, 1, x.stop))

fix_arg(x, offset) = fix_arg(x)
fix_arg(x::Union{Integer, AbstractArray{<:Integer}}, offset) = offset ? fix_arg(x .- 1) : fix_arg(x)

macro pyasvec(class)
    index_expr = quote
        function Base.getindex(pyt::$class, args...)
            offset = should_offset(pyt, args...)
            new_args = tuple([fix_arg(arg, offset) for arg in args]...)
            new_args = (length(new_args) == 1 ? new_args[1] : new_args)
            pyt.__getitem__(new_args)
        end

        Base.view(pyt::$class, is...) = Base.getindex(x, is...)

        function Base.setindex!(pyt::$class, value, idxs...)
            offset = should_offset(pyt, idxs...)
            new_idx = [fix_arg(idx, offset) for idx in idxs]
            if length(new_idx) > 1
                pyt.__setitem__(tuple(new_idx...), value)
            else
                pyt.__setitem__(new_idx[1], value)
            end
        end
    end

    length_expr = if class in [:Iloc, :Loc]
        :(Base.length(x::$class) = x.obj.__len__())
    else
        :(Base.length(x::$class) = x.__len__())
    end

    quote
        $index_expr
        $length_expr
        Base.lastindex(x::$class) = length(x)
        Base.lastindex(x::$class, i) = size(x, i)
    end
end

@pytype DataFrame () -> pd.core.frame.DataFrame
@pytype Series () -> pd.core.series.Series
@pytype Iloc () -> pd.core.indexing._iLocIndexer
@pytype Loc () -> pd.core.indexing._LocIndexer
@pytype Index () -> version() < VersionNumber(1) ? pd.core.index.Index : pd.core.indexes.multi.Index
@pytype MultiIndex () -> version() < VersionNumber(1) ? pd.core.index.MultiIndex : pd.core.indexes.multi.MultiIndex

Base.size(x::Union{Loc, Iloc}) = x.obj.shape
Base.size(df::PandasWrapped, i::Integer) = size(df)[i]
Base.size(df::PandasWrapped) = df.shape
Base.ndims(df::PandasWrapped) = length(size(df))

Base.isempty(df::PandasWrapped) = df.empty
Base.empty!(df::PandasWrapped) = df.drop(df.index, inplace = true)

should_offset(::Any, args...) = false
should_offset(::Union{Iloc, Index, MultiIndex}, args...) = true

@pyasvec Series
@pyasvec Loc
@pyasvec Iloc
@pyasvec DataFrame
@pyasvec Index
@pyasvec MultiIndex

Base.show(io::IO, df::PandasWrapped) = println(io, df.__str__())

for (jl_op, py_op, py_opᵒ) in [(:+, :__add__, :__add__), (:*, :__mul__, :__mul__),
                               (:/, :__truediv__, :__rtruediv__), (:-, :__sub__, :__rsub__),
                               (:(==), :__eq__, :__eq__), (:!=, :__ne__, :__ne__),
                               (:>, :__gt__, :__lt__), (:<, :__lt__, :__gt__),
                               (:>=, :__ge__, :__le__), (:<=, :__le__, :__ge__),
                               (:&, :__and__, :__and__), (:|, :__or__, :__or__)]
    @eval begin
        Base.$jl_op(x::PandasWrapped, y) = x.$py_op(y)
        Base.$jl_op(y, x::PandasWrapped) = x.$py_opᵒ(y)
        Base.$jl_op(x::PandasWrapped, y::PandasWrapped) = invoke($jl_op, Tuple{PandasWrapped, Any}, x, y)
    end
end

for (jl_op, py_op) in [(:-, :__neg__), (:!, :__neg__)]
    @eval Base.$jl_op(x::PandasWrapped) = x.$py_op()
end

include("operators_v6.jl")

DataFrame(pairs::Pair...) = DataFrame(Dict(pairs...))

function Base.eltype(s::Series)
    dtype_map = Dict(
        np.dtype("int8") => Int8,
        np.dtype("uint8") => UInt8,
        np.dtype("int16") => Int16,
        np.dtype("uint16") => UInt16,
        np.dtype("int32") => Int32,
        np.dtype("uint32") => UInt32,
        np.dtype("int64") => Int64,
        np.dtype("uint64") => UInt64,
        np.dtype("float16") => Float16,
        np.dtype("float32") => Float32,
        np.dtype("float64") => Float64,
        np.dtype("object") => String,
    )
    get(dtype_map, s.dtype, Any)
end

Base.eltype(df::DataFrame) = Tuple{[eltype(df[c]) for c in df.columns]...}

function Base.map(f::Function, s::Series)
    if eltype(s) <: Real
        Series([f(_) for _ in values(s)])
    else
        Series([f(_) for _ in s])
    end
end

Base.get(df::PandasWrapped, key, default) = df.get(key, default)

Base.getindex(s::Series, c::CartesianIndex{1}) = s[c[1]]

Base.copy(df::PandasWrapped) = df.copy()

include("tabletraits.jl")

function DataFrame(obj)
    y = _construct_pandas_from_iterabletable(obj)
    if y === nothing
        return invoke(DataFrame, Tuple{Vararg{Any}}, obj)
    else
        return y
    end
end

Base.getproperty(x::PandasWrapped, s::Symbol) = s == :pyo ? getfield(x, s) : getproperty(x.pyo, s)

Base.getproperty(x::PandasWrapped, s::String) = getproperty(x.pyo, s)

function Base.setproperty!(x::PandasWrapped, s::Symbol, v)
    if s == :pyo
        return setfield!(x, s, v)
    else
        setproperty!(x.pyo, s, v)
    end
end

include("miscellaneous.jl")

end


text:
module egrss
using LinearAlgebra

export generators, symv, gemv, potrf, ldl, trmv, trsv, trtri, trtri2, trnrms, full, tril

"""
    generators(t, p)

Computes generator representation of p'th order spline kernel matrix generated
by a strictly monotonic vector t of length n.

Ut,Vt = generators(t,p) returns two matrices Ut and Vt of size p-by-n
(with p > 0) such that K = tril(Ut'*Vt) + triu(Vt'*Ut,1) is the kernel
matrix with elements

    K[i,j] = sum_{k=0}^{p-1} (-1)^k/(factorial(p-1-k)*factorial(p+k)
                *(t[i]*t[j])^(p-1-k)*min(t[i],t[j])^(2*k+1)

where t is a nonnegative vector.

# Example

```
t = range(0,1,101)
Ut, Vt = generators(t, 2)
```
"""
function generators(t, p)
    Ut = repeat(t',p,1).^collect(p-1:-1:0)./map(factorial,p-1:-1:0)
    Vt = ((-1).^(0:p-1)).*(repeat(t',p,1).^collect(p:2*p-1)./map(factorial,p:2*p-1))
    return Ut,Vt
end

"""
    full(Ut, Vt[, d])

Forms dense symmetric matrix from generator representation.

K = full(Ut,Vt) forms the symmetric matrix

    K = tril(Ut'*Vt) + triu(Vt'*Ut,1).

K = full(Ut,Vt,d) forms the symmetric matrix

    K = tril(Ut'*Vt) + triu(Vt'*Ut,1) + diag(d).
"""
function full(Ut, Vt, d=Nothing)
    K = tril(Ut'*Vt,-1)
    K = K + K' + Diagonal(sum(Ut.*Vt,dims=1)[:])
    if d !== Nothing
        if isa(d,Number)
            K = K + Diagonal(d*ones(eltype(d),size(Ut,2)))
        else
            K = K + Diagonal(d[:])
        end
    end
    return Symmetric(K)
end

"""
    full_tril(Ut, Wt[, d])

Forms dense lower triangular matrix from generator representation.

L = full_tril(Ut,Wt) forms the lower triangular matrix

    L = tril(Ut'*Wt)

L = full_tril(Ut,Wt,d) forms the lower triangular matrix

    L = tril(Ut'*Wt,-1) + diag(d).
"""
function full_tril(Ut, Wt, c=Nothing)
    if c === Nothing
        return LowerTriangular(tril(Ut'*Wt))
    else
        if isa(c,Number)
            c = c*ones(eltype(c),size(Ut,2))
        end
        return LowerTriangular(tril(Ut'*Wt,-1) + Diagonal(c[:]))
    end
end

"""
    gemv(Ut, Vt, Pt, Qt, x)

Computes matrix-vector product A*x where A is an extended generator
representable semiseparable matrix given by

   A = tril(Ut'*Vt) + triu(Pt'*Qt,1).

# Example

```
b = gemv(Ut,Vt,Pt,Qt,x)
```
"""
function gemv(Ut,Vt,Pt,Qt,x)
    p, m = size(Ut)
    q, n = size(Qt)
    min_mn = min(m,n)

    b = [x[1:min_mn]; zeros(eltype(Ut),max(0,m-n))]
    v = zeros(eltype(Ut),size(Ut,1))
    q = Qt[:,1:min_mn]*x[1:min_mn]

    for k = 1:min_mn
        v = v + Vt[:,k]*b[k]
        q = q - Qt[:,k]*b[k]
        b[k] = Ut[:,k]'*v + Pt[:,k]'*q
    end

    if n > m
        b = b + Pt'*(Qt[:,min_mn+1:end]*x[min_mn+1:end])
    elseif m > n
        b[min_mn+1:end] = Ut[:,min_mn+1:end]'*(Vt*x)
    end

    return b
end

"""

Computes matrix-vector product A*x where A is a symmetric and extended generator
representable semiseparable matrix given by

    A = tril(Ut'*Vt) + triu(Vt'*Ut,1)

# Example

```
b = symv(Ut,Vt,x)
```
"""
function symv(Ut,Vt,x)
    p, n = size(Ut)
    b = copy(x)
    z = zeros(eltype(Ut),p)
    y = Ut*x
    for k = 1:n
        z = z + Vt[:,k]*b[k]
        y = y - Ut[:,k]*b[k]
        b[k] = Ut[:,k]'*z + Vt[:,k]'*y
    end
    return b
end

"""
    potrf(Ut, Vt[, d])

Computes implicit Cholesky factorization of the sum of a symmetric extended
generator representable semiseparable matrix and a diagonal matrix.

Wt = potrf(Ut,Vt) computes a matrix Wt such that L = tril(Ut'*Wt) is
the Cholesky factor of the symmetric matrix A = tril(Ut'*Vt)+triu(Vt'*Ut,1),
i.e., A = L*L'.

Wt,c = potrf(Ut,Vt,d) computes a matrix Wt and a vector c such that
L = tril(Ut'*Wt,-1) + diag(c) is the Cholesky factor of the symmetric matrix
A = tril(Ut'*Vt) + triu(Vt'*Ut,1) + diag(d), i.e., A = L*L'. The input d must
either be a vector of length size(Ut,2) or a scalar; in the latter case diag(d)
is interpreted as the identity matrix scaled by d.
"""
function potrf(Ut,Vt,d=Nothing)
    p, n = size(Ut)
    Wt = copy(Vt)
    P = zeros(eltype(Ut),p,p)
    if d === Nothing
        for k = 1:n
            Wt[:,k] = Wt[:,k] - P*Ut[:,k]
            Wt[:,k] = Wt[:,k]/sqrt(Ut[:,k]'*Wt[:,k])
            P = P + Wt[:,k]*Wt[:,k]'
        end
        return Wt
    else
        if isa(d,Number)
            c = d*ones(eltype(d),n)
        else
            c = copy(d)
        end
        for k = 1:n
            Wt[:,k] = Wt[:,k] - P*Ut[:,k]
            c[k] = sqrt(Ut[:,k]'*Wt[:,k] + c[k])
            Wt[:,k] = Wt[:,k]/c[k]
            P = P + Wt[:,k]*Wt[:,k]'
        end
        return Wt, c
    end
end

"""
    ldl(Ut, Vt[, d])

Computes implicit LDL factorization of the sum of a symmetric extended
generator representable semiseparable matrix and a diagonal matrix.

Wt,c = ldl(Ut,Vt) computes a matrix Wt and a vector c such that the LDL
factorization of the symmetric matrix A = tril(Ut'*Vt)+triu(Vt'*Ut,1) + diag(d)
is given by L = tril(Ut'*Wt,-1)+I and D = diag(c), i.e., A = L*D*L'.

Wt,c = ldl(Ut,Vt,d) computes a matrix Wt and a vector c such that the LDL
factorization of the symmetric matrix A = tril(Ut'*Vt)+triu(Vt'*Ut,1) is given
by L = tril(Ut'*Wt,-1)+I and D = diag(c), i.e., A = L*D*L'.
"""
function ldl(Ut,Vt,d=Nothing)
    p, n = size(Ut)
    Wt = copy(Vt)
    P = zeros(eltype(Ut),p,p)
    if d === Nothing
        c = zeros(eltype(Ut),n);
        for k = 1:n
            Wt[:,k] = Wt[:,k] - P*Ut[:,k]
            c[k] = Ut[:,k]'*Wt[:,k]
            P = P + (Wt[:,k]*Wt[:,k]')/c[k]
            Wt[:,k] = Wt[:,k]/c[k]
        end
    else
        if isa(d,Number)
            c = d*ones(eltype(d),n)
        else
            c = copy(d)
        end
        for k = 1:n
            Wt[:,k] = Wt[:,k] - P*Ut[:,k]
            c[k] = Ut[:,k]'*Wt[:,k] + c[k]
            P = P + (Wt[:,k]*Wt[:,k]')/c[k]
            Wt[:,k] = Wt[:,k]/c[k]
        end
    end
    return Wt, c
end

"""
Computes matrix-vector product L*x or L'*x where L is a lower triangular
extended generator representable semiseparable or quasi-separable matrix
defined in terms of matrices Ut and Wt, which are of size p-by-n (with p > 0),
and in the quasi-separable case, a vector c of length n.

b = trmv(Ut,Wt,x) computes b = L*x with L = tril(Ut'*Wt).

b = trmv(Ut,Wt,x,'N') is the same as b = trmv(Ut,Wt,x).

b = trmv(Ut,Wt,x,'T') computes b = L'*x with L = tril(Ut'*Wt).

b = trmv(Ut,Wt,c,x) computes b = L*x with L = tril(Ut'*Wt,-1) + diag(c).

b = trmv(Ut,Wt,c,x,'N') is the same as b = trmv(Ut,Wt,c,x).

b = trmv(Ut,Wt,c,x,'T') computes b = L'*x with L = tril(Ut'*Wt,-1) + diag(c).
"""
function trmv(Ut,Wt,args...)
    p, n = size(Ut)

    if length(args) == 3
        c = args[1]
        x = args[2]
        trans = args[3]
    elseif length(args) == 2
        if args[2] isa Char
            c = Nothing
            x = args[1]
            trans = args[2]
        else
            c = args[1]
            x = args[2]
            trans = 'N'
        end
    elseif length(args) == 1
        c = Nothing
        x = args[1]
        trans = 'N'
    else
        error("Invalid number of input arguments.")
    end

    b = copy(x)
    z = zeros(eltype(Ut),p)
    if c === Nothing
        if trans == 'N'
            for k = 1:n
                z = z + Wt[:,k]*x[k]
                b[k] = Ut[:,k]'*z
            end
        elseif trans == 'T'
            for k = n:-1:1
                z = z + Ut[:,k]*x[k]
                b[k] = Wt[:,k]'*z
            end
        else
            error("Expected character 'N or 'T'.")
        end
    else
        if trans == 'N'
            for k = 1:n
                b[k] = c[k]*b[k] + Ut[:,k]'*z
                z = z + Wt[:,k]*x[k]
            end
        elseif trans == 'T'
            for k = n:-1:1
                b[k] = c[k]*b[k] + Wt[:,k]'*z
                z = z + Ut[:,k]*x[k]
            end
        else
            error("Expected character 'N or 'T'.")
        end
    end

    return b
end

"""
Solves the equation L*x = b or L'*x = b where L is a lower triangular extended
generator representable semiseparable or quasi-separable matrix defined in terms
of matrices Ut and Wt, which are of size p-by-n (with p > 0), and in the
quasi-separable case, a vector c of length n.

x = trsv(Ut,Wt,b) solves L*x = b with L = tril(Ut'*Wt).

x = trsv(Ut,Wt,b,'N') is the same as x = trsv(Ut,Wt,b).

x = trsv(Ut,Wt,b,'T') solves L'*x = b with L = tril(Ut'*Wt).

x = trsv(Ut,Wt,c,b) solves L*x = b with L = tril(Ut'*Wt,-1) + diag(c).

x = trsv(Ut,Wt,c,b,'N') is the same as x = trsv(Ut,Wt,c,b).

x = trsv(Ut,Wt,c,b,'T') solves L'*x = b with L = tril(Ut'*Wt,-1) + diag(c).
"""
function trsv(Ut,Wt,args...)
    p, n = size(Ut)

    if length(args) == 3
        c = args[1]
        b = args[2]
        trans = args[3]
    elseif length(args) == 2
        if args[2] isa Char
            c = Nothing
            b = args[1]
            trans = args[2]
        else
            c = args[1]
            b = args[2]
            trans = 'N'
        end
    elseif length(args) == 1
        c = Nothing
        b = args[1]
        trans = 'N'
    else
        error("Invalid number of input arguments.")
    end

    x = copy(b)
    z = zeros(eltype(Ut),p);
    if c === Nothing
        if trans == 'N'
            for k = 1:n
                x[k] = (x[k]-Ut[:,k]'*z)/(Wt[:,k]'*Ut[:,k])
                z = z + Wt[:,k]*x[k]
            end
        elseif trans == 'T'
            for k = n:-1:1
                x[k] = (x[k]-Wt[:,k]'*z)/(Wt[:,k]'*Ut[:,k])
                z = z + Ut[:,k]*x[k]
            end
        else
            error("Expected character 'N or 'T'.")
        end
    else
        if trans == 'N'
            for k = 1:n
                x[k] = (x[k]-Ut[:,k]'*z)/c[k]
                z = z + Wt[:,k]*x[k]
            end
        elseif trans == 'T'
            for k = n:-1:1
                x[k] = (x[k]-Wt[:,k]'*z)/c[k]
                z = z + Ut[:,k]*x[k]
            end
        else
            error("Expected character 'N or 'T'.")
        end
    end

    return x
end

"""
    trtri(Ut, Wt, c)

Computes an implicit inverse of a lower triangular matrix of the form

    L = tril(Ut'*Wt,-1) + diag(c).

Yt,Zt = trtri(Ut,Wt,c) computes Yt and Zt such that

    L^{-1} = tril(Yt'*Zt,-1) + diag(1./c).
"""
function trtri(Ut,Wt,c)

    p, n = size(Ut)
    Yt = copy(Ut)
    Zt = copy(Wt)

    for k = 1:p
        Yt[k,:] = trsv(Ut,Wt,c,Yt[k,:]')';
    end

    Zt[:,end] .*= 0.0;
    for k = 1:p
        Zt[k,:] = trsv(Ut,Wt,c,Zt[k,:]','T')';
    end
    Zt = (Zt[:,1:end-1]*Ut[:,1:end-1]' - I)\Zt;

    return Yt,Zt
end

"""
    trtri2(Ut, Wt[, c])

Computes the explicit inverse of a lower triangular matrix of either the form
L = tril(Ut'*Wt) or L = tril(Ut'*Wt,-1) + diag(c).

"""
function trtri2(Ut,Wt,c=Nothing)
    p, n = size(Ut)
    Linv = zeros(eltype(Ut),n,n)
    e = zeros(eltype(Ut),n)
    e[1] = 1
    if c === Nothing
        for k = 1:n
            Linv[k:end,:] = trsv(Ut,Wt,e[1:end+1-k])
        end
    else
        for k = 1:n
            Linv[k:end,k] = trsv(Ut[:,k:end],Wt[:,k:end],c[k:end],e[1:end+1-k])
        end
    end
    return LowerTriangular(Linv)
end

"""
    trnrms(Ut, Wt, c)

Computes the squared column norms of a lower triangular matrix of the form
L = tril(Ut'*Wt,-1) + diag(c) where Ut and Wt are of size p-by-n (with p > 0),
and c is a positive vector of length n.

b = trnrms(Ut,Wt,c) computes a vector b and with entries
b(i) = norm(L(:,i))^2 where L = tril(Ut'*Wt,-1) + diag(c).
"""
function trnrms(Ut,Wt,c)
    p, n = size(Ut)
    colnrms = zeros(eltype(Ut),n)
    P = zeros(eltype(Ut),p,p);
    for k = n:-1:1
        colnrms[k] = c[k]^2 + Wt[:,k]'*P*Wt[:,k];
        P = P + Ut[:,k]*Ut[:,k]';
    end
    return colnrms
end

end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
eye(n) = LinearAlgebra.I + zeros(n, n)

mutable struct PropCov
  value::Matrix{Float64}
  buffer::Matrix{Float64}
  iter::Int
  batchsize::Int
  scale::Float64
  momentum::Float64
  eps::Float64
end

function PropCov(value::Matrix{Float64}; iter::Int=0,
                 batchsize::Int=200, momentum::Float64=0.1,
                 eps::Float64=1e-6)
  # Get dimensions of prop cov
  K, K2 = size(value)
  @assert K == K2

  # Momentum needs to be between 0 and 1
  @assert 0 <= momentum <= 1

  # Create buffer
  buffer = zeros(Float64, batchsize, K)

  # Scale for proposal covariance
  scale = 2.38 * 2.38 / K

  return PropCov(value, buffer, iter, batchsize, scale, momentum, eps)
end

function PropCov(K::Int; sigma::Float64=1.0, iter::Int=0, batchsize::Int=200,
                 momentum::Float64=0.1, eps::Float64=1e-6)
  return PropCov(eye(K) * sigma, iter=iter, batchsize=batchsize,
                 momentum=momentum, eps=eps)
end

nparam(pc::PropCov) = size(pc.buffer, 2)

function update!(pc::PropCov, x::Vector{Float64})
  # Increment iteration
  pc.iter += 1

  # Compute current buffer index
  buffer_index = let
    idx = mod(pc.iter, pc.batchsize)
    idx == 0 ? pc.batchsize : idx
  end

  # Cache current value in buffer
  pc.buffer[buffer_index, :] .= deepcopy(x)

  if buffer_index == pc.batchsize
    # Update prop cov
    pc.value .= let
      S = cov(pc.buffer) * pc. scale * pc.momentum + pc.value * (1 - pc.momentum)
      S + eye(nparam(pc)) * pc.eps
    end
  end
end

# K = 30
# pc = Metropolis.PropCov(K, sigma=1e-2)
# for i in 1:10000 Metropolis.update!(pc, randn(K)) end
# pc.value


text:
mutable struct SmallLLLRelationsCtx{T}
  A::NfOrdIdl
  b::Array{T, 1}
  bd::Int
  cnt::Int
  elt::T
  function SmallLLLRelationsCtx(a::T) where {T}
    n = new{T}()
    n.bd = 1
    n.cnt = 0
    return n
  end
end

mutable struct NormCtx_split <: NormCtx
  O::NfAbsOrd{AnticNumberField, nf_elem}
  lp::Array{Int, 1}  #the primes
  lR::Array{GaloisField, 1} #the local (finite field) splitting field
  nb::Int #bound on the number of bits the norm is allowed to have
  lC::Array{gfp_mat, 1} # for each p in lp, the conjugates of the basis of O
  mp::Array{gfp_mat, 1} # for each p in lp, the conjugates of the basis of O
  np::Array{gfp_mat, 1} # for each p in lp, the conjugates of the basis of O
  e::crt_env
  function NormCtx_split(O::NfAbsOrd, nb::Int)
    p = p_start
    NC = new()
    NC.O = O
    NC.nb = nb
    NC.lp = Array{Int, 1}()
    NC.lR = Array{GaloisField, 1}()
    NC.lC = Array{gfp_mat, 1}()
    NC.mp = Array{gfp_mat, 1}()
    NC.np = Array{gfp_mat, 1}()
    B = basis(O)
    n = degree(O)
    while (true)
      p = next_prime(p)
      lP = prime_decomposition(O, p)
      if length(lP) < degree(O)
        continue
      end
      push!(NC.lp, p)
      R = GF(p, cached = false)
      push!(NC.lR, R)
      push!(NC.mp, zero_matrix(R, 1, degree(O)))
      push!(NC.np, zero_matrix(R, 1, degree(O)))
      M = Array{Nemo.elem_type(R), 1}()
      for P = lP
        F, mF = ResidueFieldSmall(O, P[1])
        for x = B
          xp = mF(x)
          push!(M, R(coeff(xp, 0)))
        end
      end
      push!(NC.lC, matrix(R, n, n, M)')
      nb -= nbits(p)
      if nb <= 0
        NC.e = crt_env(fmpz[p for p = NC.lp])
        return NC
      end
    end
  end
end

function NormCtx(O::NfAbsOrd, nb::Int, normal::Bool = false)
  if normal
    return NormCtx_split(O, nb)
  else
    return NormCtx_gen(O, nb)
  end
end

function norm(m::fmpz_mat, NC::NormCtx_split, div::fmpz = fmpz(1))
  l = Array{fmpz, 1}()
  for i = 1:length(NC.lp)
    ccall((:fmpz_mat_get_nmod_mat, libflint), Cvoid, (Ref{gfp_mat}, Ref{fmpz_mat}), NC.mp[i], m)
    mul!(NC.np[i], NC.mp[i], NC.lC[i])
    n = NC.np[i]
    p = n[1,1]
    for j=2:ncols(m)
      p *= n[1,j]
    end
    push!(l, lift(inv(NC.lR[i](div))*p))
  end
  no = crt_signed(l, NC.e)
  if nbits(no) > NC.nb - 10
    return nothing
  else
    return no
  end
end

#not used.
function norms(a::fmpz_mat, NC::NormCtx_split, div::fmpz = fmpz(1))
  no = Array{fmpz, 1}()
  nr = Array{gfp_mat, 1}()
  for i=1:length(NC.lp)
    n = matrix(NC.lR[i], a)*NC.lC[i]
    m = zero_matrix(NC.lR[i], nrows(a), 1)
    for k=1:nrows(n)
      for j=2:ncols(n)
        n[k, 1] *= n[k, j]
      end
      m[k, 1] = n[k, 1]*inv(NC.lR[i](div))
    end
    push!(nr, m)
  end
  mm = induce_crt(nr, NC.e, true)
end

mutable struct NormCtx_gen <: NormCtx
  nb::Int
  O::NfAbsOrd
  function NormCtx_gen(O::NfAbsOrd, nb::Int)
    NC = new()
    NC.nb = nb
    NC.O = O
    return NC
  end
end

function norm(a::fmpz_mat, NC::NormCtx_gen, div::fmpz = fmpz(1))
  O = NC.O
  k = nf(O)
  no = numerator(norm_div(O(fmpz[a[1, i] for i = 1:degree(k)]).elem_in_nf, div, NC.nb))
  if nbits(no) > NC.nb-10
    return nothing
  else
    return no
  end
end

function class_group_small_lll_elements_relation_start(clg::ClassGrpCtx{T},
                O::NfOrd; prec::Int = 200, val::Int = 0,
                limit::Int = 0) where {T}
  return class_group_small_lll_elements_relation_start(clg, ideal(O, parent(basis_matrix(O).num)(1)), prec = prec)
end

function class_group_small_lll_elements_relation_start(clg::ClassGrpCtx{T},
                A::NfOrdIdl; prec::Int = 200, val::Int = 0,
                limit::Int = 0) where {T}
  global _start = A
  O = order(A)
  n = degree(O)
  L, Tr = lll(A, prec = prec)
  I = SmallLLLRelationsCtx(zero_matrix(FlintZZ, 1, 1))
  S = Tr*basis_matrix(A, copy = false)
  bd = abs(discriminant(O))*norm(A)^2
  bd = root(bd, n)
  bd *= L.den
  f = Int[i for i = 1:n if compare_index(L.num, i, i, bd) < 0]
  m = div(n, 4)
  if m < 2
    m = n
  end
  while length(f) < m
    f = Int[i for i = 1:n if compare_index(L.num, i, i, bd) < 0]
    bd *= 2
  end
  I.b = Vector{fmpz_mat}(undef, length(f))
  for j = 1:length(f)
    I.b[j] = view(S, f[j]:f[j], 1:ncols(S))
  end
  #println([Float64(numerator(L)[i,i]//denominator(L)*1.0) for i=1:degree(K)])
  #now select a subset that can yield "small" relations, where
  #small means of effective norm <= sqrt(disc)
  I.A = A
  I.elt = zero_matrix(FlintZZ, 1, n)
  return I
end

function class_group_small_lll_elements_relation_next(I::SmallLLLRelationsCtx)
  #the 1st n relations are the basis elements
  if I.cnt < length(I.b)
    I.cnt += 1
    I.elt = I.b[I.cnt]
    return deepcopy(I.b[I.cnt])
  end
  #the next 2*n*(n-1)/2 relations are the ones of weight 2
  #1st b[i] + b[j], all combinations, then b[i]-b[j]
  #it may (or may not) be helpful
  if I.cnt <= length(I.b)^2
    n = length(I.b)
    c = I.cnt - n # we already did n relations in the 1st if
    if c > n*(n-1)/2
      c -= div(n*(n-1), 2)
      s = -1
    else
      s = 1
    end
    i = 1
    while c+i-1 >= n
      c -= (n-i)
      i += 1
    end
    I.cnt += 1
    I.elt = I.b[i] + s*I.b[c+i]
    return I.elt
  end

  if I.cnt > (2*I.bd+1)^div(length(I.b), 2)
    I.bd += 1
  end

  I.cnt += 1
  while true
    # Before: rand!(I.elt, I.b, -I.bd:I.bd, min(length(I.b), 5))
    # TODO: make this non-allocating
    I.elt = rand(-I.bd:I.bd) * rand(I.b)
    for i in 2:min(length(I.b), 5)
      I.elt = I.elt + rand(-I.bd:I.bd) * rand(I.b)
    end
    if !iszero(I.elt)
      return deepcopy(I.elt)
    end
  end
end


text:
#!/bin/julia

using DataFrames
using CSV
using Dates

include("moveSim_config")

### Load data
println("Loading data...")

#commuter network
cnet = DataFrame(CSV.File(CNET_FILE))

#population estimates
pop = DataFrame(CSV.File(POP_FILE))

### Construct population table columns
## Required: {PersonID, Origin, Age} + {PersonID, Origin, Age, Worker, Destination}

p = sum(pop.All) #total population

#Initialise DataFrame
ss = repeat([""],p) #Empty string column
bs = repeat([false],p) #Empty Bool column
df = DataFrame(PersonID=1:p, Origin=ss, Age=ss, Worker=bs, Destination=ss)

#Origins
function add_origins!(df,pop)
    areas = nrow(pop)
    i = 1
    for area in 1:areas
        n = pop[area,:All]
        df[i:i+n-1,:Origin] .= pop.RefArea[area]
        i += n
    end
end
println("Origins...")
@time add_origins!(df,pop)

#Age
function add_ages!(df,pop)
    n_areas = nrow(pop)
    i = 1
    for area in 1:n_areas
        for age in [:Young,:Adult,:Elderly] #Assumes age cats
            n = pop[area,age]
            df[i:i+n-1,:Age] .= string(age)
            i += n
        end
    end
end
println("Ages...")
@time add_ages!(df,pop)

#Worker
#Assuming all adults work
println("Workers...")
@time df[df.Age.=="Adult",:Worker] .= true

#Destination
function add_destinations!(df,pop,cnet)
    #Get weights
    #If no weights in file then calculate
    if !(:weightedFromOrigin in names(cnet))
        wfo = by(cnet,:origin, weightedFromOrigin = :numberFromOrigin => x -> x ./ sum(x))
        cnet.weightedFromOrigin = wfo.weightedFromOrigin
    end
    #make weight lookup
    weight = Dict((cnet[i,:origin],cnet[i,:destination])
                  => cnet[i,:weightedFromOrigin] for i = 1:nrow(cnet))
    #get number of workers
    n_workers = sum(pop.Working)
    #initialise worker destiantion array
    worker_dests = Array{String}(undef,n_workers) .= ""
    #assign workers
    n_areas = nrow(pop)
    i = 1
    for orig in 1:n_areas
        i_origin = i
        n_from_origin = pop[orig,:Working]
        for dest in 1:n_areas
            o_d = (pop[orig,:RefArea],pop[dest,:RefArea])
            if (orig != dest) && (o_d in keys(weight))
                n::Int64 = floor(n_from_origin * weight[o_d])
                worker_dests[i:i+n-1] .= pop[dest,:RefArea]
                i += n
            end
        end
        i = i_origin + n_from_origin
        #print(round((i/n_workers)*100,digits=2),"%","\r")
        print(i," / ",n_workers," workers, ",round((i/n_workers)*100,digits=2),"%","\r")
    end
    df[df.Worker,:Destination] = worker_dests
end

println("Destinations...")
@time add_destinations!(df,pop,cnet)

println("Writing population files...")
CSV.write(OUT_DIR*"population_data.csv",df)
CSV.write(OUT_DIR*"population.csv",df[:,[:PersonID,:Origin,:Age]])

max_stars_repo_name: Kao-Group/SCoVMod

max_issues_repo_name: Kao-Group/SCoVMod

max_forks_repo_name: Kao-Group/SCoVMod


text:
import LinearAlgebra

"""
LowPrecArray

An n-dimensional array that behaves essentially like a Float32 array for all
scalar operations, but matmuls are performed in BFloat16 with Float32
accumulation. This is intended to match the behavior of TPUs.
"""
struct LowPrecArray{N, A<:AbstractArray{Float32}} <: AbstractArray{Float32, N}
    storage::A
end
LowPrecArray(A::AbstractArray{Float32, N}) where {N} =
    LowPrecArray{N, typeof(A)}(A)

Base.similar(a::LowPrecArray, b::Type{Float32}, dims::Dims) =
    LowPrecArray(similar(a.storage, b, dims))
Base.size(a::LowPrecArray) = size(a.storage)
Base.getindex(A::LowPrecArray, inds...) = getindex(A.storage, inds...)
Base.setindex!(A::LowPrecArray, args...) = setindex!(A.storage, args...)

# A BFloat16 whose multiplication is a widemul
struct ExpandingBFloat16
    a::BFloat16
end
ExpandingBFloat16(a::Float32) = ExpandingBFloat16(BFloat16(a))
*(a::ExpandingBFloat16, b::ExpandingBFloat16) = widemul(a.a, b.a)

struct MatMulView{A} <: AbstractArray{ExpandingBFloat16, 2}
    a::A
end
Base.size(a::MatMulView) = size(a.a)
Base.getindex(A::MatMulView, inds...) = ExpandingBFloat16(getindex(A.a, inds...))
# setindex! deliberately not defined

import LinearAlgebra: generic_matmatmul!
function generic_matmatmul!(C::AbstractMatrix{Float32}, ta, tb,
    A::Union{LowPrecArray, AbstractMatrix{BFloat16}},
    B::Union{LowPrecArray, AbstractMatrix{BFloat16}})
    generic_matmatmul!(C, ta, tb,
        MatMulView(A), MatMulView(B))
end

# For now do gemv! in Float32 precision - we'll have to evaluate if this makes
# sense on the real hardware (which would have to do this on the vector units)
function LinearAlgebra.BLAS.gemv!(trans::AbstractChar, alpha::Float32,
                                  A::AbstractVecOrMat{Float32}, X::AbstractVector{Float32},
                                  beta::Float32, Y::LowPrecArray)
    LinearAlgebra.BLAS.gemv!(trans, alpha, A, X, beta, Y.storage)
end