text:
#   Generating #define constants
const CURL_TYPEOF_CURL_OFF_T             = Clong
const CURL_FORMAT_CURL_OFF_T             = "ld"
const CURL_FORMAT_CURL_OFF_TU            = "lu"
const CURL_SUFFIX_CURL_OFF_T             = Clong
const CURL_SUFFIX_CURL_OFF_TU            = Culong
const CURL_TYPEOF_CURL_SOCKLEN_T         = socklen_t
const CURL_PULL_SYS_TYPES_H              = 1
const CURL_PULL_SYS_SOCKET_H             = 1
const CURL_SOCKET_BAD                    = -1
const CURLSSLBACKEND_LIBRESSL            = CURLSSLBACKEND_OPENSSL
const CURLSSLBACKEND_BORINGSSL           = CURLSSLBACKEND_OPENSSL
const CURLSSLBACKEND_CYASSL              = CURLSSLBACKEND_WOLFSSL
const CURL_HTTPPOST_FILENAME             = (1<<0)
const CURL_HTTPPOST_READFILE             = (1<<1)
const CURL_HTTPPOST_PTRNAME              = (1<<2)
const CURL_HTTPPOST_PTRCONTENTS          = (1<<3)
const CURL_HTTPPOST_BUFFER               = (1<<4)
const CURL_HTTPPOST_PTRBUFFER            = (1<<5)
const CURL_HTTPPOST_CALLBACK             = (1<<6)
const CURL_HTTPPOST_LARGE                = (1<<7)
const CURL_MAX_READ_SIZE                 = 524288
const CURL_MAX_WRITE_SIZE                = 16384
const CURL_MAX_HTTP_HEADER               = (100*1024)
const CURL_WRITEFUNC_PAUSE               = 0x10000001
const CURLFINFOFLAG_KNOWN_FILENAME       = (1<<0)
const CURLFINFOFLAG_KNOWN_FILETYPE       = (1<<1)
const CURLFINFOFLAG_KNOWN_TIME           = (1<<2)
const CURLFINFOFLAG_KNOWN_PERM           = (1<<3)
const CURLFINFOFLAG_KNOWN_UID            = (1<<4)
const CURLFINFOFLAG_KNOWN_GID            = (1<<5)
const CURLFINFOFLAG_KNOWN_SIZE           = (1<<6)
const CURLFINFOFLAG_KNOWN_HLINKCOUNT     = (1<<7)
const CURL_CHUNK_BGN_FUNC_OK             = 0
const CURL_CHUNK_BGN_FUNC_FAIL           = 1
const CURL_CHUNK_BGN_FUNC_SKIP           = 2
const CURL_CHUNK_END_FUNC_OK             = 0
const CURL_CHUNK_END_FUNC_FAIL           = 1
const CURL_FNMATCHFUNC_MATCH             = 0
const CURL_FNMATCHFUNC_NOMATCH           = 1
const CURL_FNMATCHFUNC_FAIL              = 2
const CURL_SEEKFUNC_OK                   = 0
const CURL_SEEKFUNC_FAIL                 = 1
const CURL_SEEKFUNC_CANTSEEK             = 2
const CURL_READFUNC_ABORT                = 0x10000000
const CURL_READFUNC_PAUSE                = 0x10000001
const CURL_SOCKOPT_OK                    = 0
const CURL_SOCKOPT_ERROR                 = 1
const CURL_SOCKOPT_ALREADY_CONNECTED     = 2
const CURLE_OBSOLETE16                   = CURLE_HTTP2
const CURLE_OBSOLETE10                   = CURLE_FTP_ACCEPT_FAILED
const CURLE_OBSOLETE12                   = CURLE_FTP_ACCEPT_TIMEOUT
const CURLOPT_ENCODING                   = CURLOPT_ACCEPT_ENCODING
const CURLE_FTP_WEIRD_SERVER_REPLY       = CURLE_WEIRD_SERVER_REPLY
const CURLE_UNKNOWN_TELNET_OPTION        = CURLE_UNKNOWN_OPTION
const CURLE_SSL_PEER_CERTIFICATE         = CURLE_PEER_FAILED_VERIFICATION
const CURLE_OBSOLETE                     = CURLE_OBSOLETE50
const CURLE_BAD_PASSWORD_ENTERED         = CURLE_OBSOLETE46
const CURLE_BAD_CALLING_ORDER            = CURLE_OBSOLETE44
const CURLE_FTP_USER_PASSWORD_INCORRECT  = CURLE_OBSOLETE10
const CURLE_FTP_CANT_RECONNECT           = CURLE_OBSOLETE16
const CURLE_FTP_COULDNT_GET_SIZE         = CURLE_OBSOLETE32
const CURLE_FTP_COULDNT_SET_ASCII        = CURLE_OBSOLETE29
const CURLE_FTP_WEIRD_USER_REPLY         = CURLE_OBSOLETE12
const CURLE_FTP_WRITE_ERROR              = CURLE_OBSOLETE20
const CURLE_LIBRARY_NOT_FOUND            = CURLE_OBSOLETE40
const CURLE_MALFORMAT_USER               = CURLE_OBSOLETE24
const CURLE_SHARE_IN_USE                 = CURLE_OBSOLETE57
const CURLE_URL_MALFORMAT_USER           = CURLE_NOT_BUILT_IN
const CURLE_FTP_ACCESS_DENIED            = CURLE_REMOTE_ACCESS_DENIED
const CURLE_FTP_COULDNT_SET_BINARY       = CURLE_FTP_COULDNT_SET_TYPE
const CURLE_FTP_QUOTE_ERROR              = CURLE_QUOTE_ERROR
const CURLE_TFTP_DISKFULL                = CURLE_REMOTE_DISK_FULL
const CURLE_TFTP_EXISTS                  = CURLE_REMOTE_FILE_EXISTS
const CURLE_HTTP_RANGE_ERROR             = CURLE_RANGE_ERROR
const CURLE_FTP_SSL_FAILED               = CURLE_USE_SSL_FAILED
const CURLE_OPERATION_TIMEOUTED          = CURLE_OPERATION_TIMEDOUT
const CURLE_HTTP_NOT_FOUND               = CURLE_HTTP_RETURNED_ERROR
const CURLE_HTTP_PORT_FAILED             = CURLE_INTERFACE_FAILED
const CURLE_FTP_COULDNT_STOR_FILE        = CURLE_UPLOAD_FAILED
const CURLE_FTP_PARTIAL_FILE             = CURLE_PARTIAL_FILE
const CURLE_FTP_BAD_DOWNLOAD_RESUME      = CURLE_BAD_DOWNLOAD_RESUME
const CURLE_ALREADY_COMPLETE             = 99999
const CURLOPT_FILE                       = CURLOPT_WRITEDATA
const CURLOPT_INFILE                     = CURLOPT_READDATA
const CURLOPT_WRITEHEADER                = CURLOPT_HEADERDATA
const CURLOPT_WRITEINFO                  = CURLOPT_OBSOLETE40
const CURLOPT_CLOSEPOLICY                = CURLOPT_OBSOLETE72
const CURLAUTH_NONE                      = (0)
const CURLAUTH_BASIC                     = ((1)<<0)
const CURLAUTH_DIGEST                    = ((1)<<1)
const CURLAUTH_NEGOTIATE                 = ((1)<<2)
const CURLAUTH_GSSNEGOTIATE              = CURLAUTH_NEGOTIATE
const CURLAUTH_GSSAPI                    = CURLAUTH_NEGOTIATE
const CURLAUTH_NTLM                      = ((1)<<3)
const CURLAUTH_DIGEST_IE                 = ((1)<<4)
const CURLAUTH_NTLM_WB                   = ((1)<<5)
const CURLAUTH_BEARER                    = ((1)<<6)
const CURLAUTH_ONLY                      = ((1)<<31)
const CURLAUTH_ANY                       = (~CURLAUTH_DIGEST_IE)
const CURLAUTH_ANYSAFE                   = (~(CURLAUTH_BASIC|CURLAUTH_DIGEST_IE))
const CURLSSH_AUTH_ANY                   = ~0
const CURLSSH_AUTH_NONE                  = 0
const CURLSSH_AUTH_PUBLICKEY             = (1<<0)
const CURLSSH_AUTH_PASSWORD              = (1<<1)
const CURLSSH_AUTH_HOST                  = (1<<2)
const CURLSSH_AUTH_KEYBOARD              = (1<<3)
const CURLSSH_AUTH_AGENT                 = (1<<4)
const CURLSSH_AUTH_GSSAPI                = (1<<5)
const CURLSSH_AUTH_DEFAULT               = CURLSSH_AUTH_ANY
const CURLGSSAPI_DELEGATION_NONE         = 0
const CURLGSSAPI_DELEGATION_POLICY_FLAG  = (1<<0)
const CURLGSSAPI_DELEGATION_FLAG         = (1<<1)
const CURL_ERROR_SIZE                    = 256
const CURLSSLOPT_ALLOW_BEAST             = (1<<0)
const CURLSSLOPT_NO_REVOKE               = (1<<1)
const CURL_HET_DEFAULT                   = Clong(200)
const CURLFTPSSL_NONE                    = CURLUSESSL_NONE
const CURLFTPSSL_TRY                     = CURLUSESSL_TRY
const CURLFTPSSL_CONTROL                 = CURLUSESSL_CONTROL
const CURLFTPSSL_ALL                     = CURLUSESSL_ALL
const CURLFTPSSL_LAST                    = CURLUSESSL_LAST
const CURLHEADER_UNIFIED                 = 0
const CURLHEADER_SEPARATE                = (1<<0)
const CURLPROTO_HTTP                     = (1<<0)
const CURLPROTO_HTTPS                    = (1<<1)
const CURLPROTO_FTP                      = (1<<2)
const CURLPROTO_FTPS                     = (1<<3)
const CURLPROTO_SCP                      = (1<<4)
const CURLPROTO_SFTP                     = (1<<5)
const CURLPROTO_TELNET                   = (1<<6)
const CURLPROTO_LDAP                     = (1<<7)
const CURLPROTO_LDAPS                    = (1<<8)
const CURLPROTO_DICT                     = (1<<9)
const CURLPROTO_FILE                     = (1<<10)
const CURLPROTO_TFTP                     = (1<<11)
const CURLPROTO_IMAP                     = (1<<12)
const CURLPROTO_IMAPS                    = (1<<13)
const CURLPROTO_POP3                     = (1<<14)
const CURLPROTO_POP3S                    = (1<<15)
const CURLPROTO_SMTP                     = (1<<16)
const CURLPROTO_SMTPS                    = (1<<17)
const CURLPROTO_RTSP                     = (1<<18)
const CURLPROTO_RTMP                     = (1<<19)
const CURLPROTO_RTMPT                    = (1<<20)
const CURLPROTO_RTMPE                    = (1<<21)
const CURLPROTO_RTMPTE                   = (1<<22)
const CURLPROTO_RTMPS                    = (1<<23)
const CURLPROTO_RTMPTS                   = (1<<24)
const CURLPROTO_GOPHER                   = (1<<25)
const CURLPROTO_SMB                      = (1<<26)
const CURLPROTO_SMBS                     = (1<<27)
const CURLPROTO_ALL                      = (~0)
const CURLOPTTYPE_LONG                   = 0
const CURLOPTTYPE_OBJECTPOINT            = 10000
const CURLOPTTYPE_STRINGPOINT            = 10000
const CURLOPTTYPE_FUNCTIONPOINT          = 20000
const CURLOPTTYPE_OFF_T                  = 30000
const CURLOPT_XFERINFODATA               = CURLOPT_PROGRESSDATA
const CURLOPT_SERVER_RESPONSE_TIMEOUT    = CURLOPT_FTP_RESPONSE_TIMEOUT
const CURLOPT_POST301                    = CURLOPT_POSTREDIR
const CURLOPT_SSLKEYPASSWD               = CURLOPT_KEYPASSWD
const CURLOPT_FTPAPPEND                  = CURLOPT_APPEND
const CURLOPT_FTPLISTONLY                = CURLOPT_DIRLISTONLY
const CURLOPT_FTP_SSL                    = CURLOPT_USE_SSL
const CURLOPT_SSLCERTPASSWD              = CURLOPT_KEYPASSWD
const CURLOPT_KRB4LEVEL                  = CURLOPT_KRBLEVEL
const CURL_IPRESOLVE_WHATEVER            = 0
const CURL_IPRESOLVE_V4                  = 1
const CURL_IPRESOLVE_V6                  = 2
const CURLOPT_RTSPHEADER                 = CURLOPT_HTTPHEADER
const CURL_HTTP_VERSION_2                = CURL_HTTP_VERSION_2_0
const CURL_REDIR_GET_ALL                 = 0
const CURL_REDIR_POST_301                = 1
const CURL_REDIR_POST_302                = 2
const CURL_REDIR_POST_303                = 4
const CURL_REDIR_POST_ALL                = (CURL_REDIR_POST_301|CURL_REDIR_POST_302|CURL_REDIR_POST_303)
const CURL_ZERO_TERMINATED               = reinterpret(Csize_t, -1)
const CURLINFO_STRING                    = 0x100000
const CURLINFO_LONG                      = 0x200000
const CURLINFO_DOUBLE                    = 0x300000
const CURLINFO_SLIST                     = 0x400000
const CURLINFO_PTR                       = 0x400000
const CURLINFO_SOCKET                    = 0x500000
const CURLINFO_OFF_T                     = 0x600000
const CURLINFO_MASK                      = 0x0fffff
const CURLINFO_TYPEMASK                  = 0xf00000
const CURLINFO_HTTP_CODE                 = CURLINFO_RESPONSE_CODE
const CURL_GLOBAL_SSL                    = (1<<0)
const CURL_GLOBAL_WIN32                  = (1<<1)
const CURL_GLOBAL_ALL                    = (CURL_GLOBAL_SSL|CURL_GLOBAL_WIN32)
const CURL_GLOBAL_NOTHING                = 0
const CURL_GLOBAL_DEFAULT                = CURL_GLOBAL_ALL
const CURL_GLOBAL_ACK_EINTR              = (1<<2)
const CURLVERSION_NOW                    = CURLVERSION_FIFTH
const CURL_VERSION_IPV6                  = (1<<0)
const CURL_VERSION_KERBEROS4             = (1<<1)
const CURL_VERSION_SSL                   = (1<<2)
const CURL_VERSION_LIBZ                  = (1<<3)
const CURL_VERSION_NTLM                  = (1<<4)
const CURL_VERSION_GSSNEGOTIATE          = (1<<5)
const CURL_VERSION_DEBUG                 = (1<<6)
const CURL_VERSION_ASYNCHDNS             = (1<<7)
const CURL_VERSION_SPNEGO                = (1<<8)
const CURL_VERSION_LARGEFILE             = (1<<9)
const CURL_VERSION_IDN                   = (1<<10)
const CURL_VERSION_SSPI                  = (1<<11)
const CURL_VERSION_CONV                  = (1<<12)
const CURL_VERSION_CURLDEBUG             = (1<<13)
const CURL_VERSION_TLSAUTH_SRP           = (1<<14)
const CURL_VERSION_NTLM_WB               = (1<<15)
const CURL_VERSION_HTTP2                 = (1<<16)
const CURL_VERSION_GSSAPI                = (1<<17)
const CURL_VERSION_KERBEROS5             = (1<<18)
const CURL_VERSION_UNIX_SOCKETS          = (1<<19)
const CURL_VERSION_PSL                   = (1<<20)
const CURL_VERSION_HTTPS_PROXY           = (1<<21)
const CURL_VERSION_MULTI_SSL             = (1<<22)
const CURL_VERSION_BROTLI                = (1<<23)
const CURLPAUSE_RECV                     = (1<<0)
const CURLPAUSE_RECV_CONT                = (0)
const CURLPAUSE_SEND                     = (1<<2)
const CURLPAUSE_SEND_CONT                = (0)
const CURLPAUSE_ALL                      = (CURLPAUSE_RECV|CURLPAUSE_SEND)
const CURLPAUSE_CONT                     = (CURLPAUSE_RECV_CONT|CURLPAUSE_SEND_CONT)
const CURLM_CALL_MULTI_SOCKET            = CURLM_CALL_MULTI_PERFORM
const CURLPIPE_NOTHING                   = Clong(0)
const CURLPIPE_HTTP1                     = Clong(1)
const CURLPIPE_MULTIPLEX                 = Clong(2)
const CURL_WAIT_POLLIN                   = 0x0001
const CURL_WAIT_POLLPRI                  = 0x0002
const CURL_WAIT_POLLOUT                  = 0x0004
const CURL_POLL_NONE                     = 0
const CURL_POLL_IN                       = 1
const CURL_POLL_OUT                      = 2
const CURL_POLL_INOUT                    = 3
const CURL_POLL_REMOVE                   = 4
const CURL_SOCKET_TIMEOUT                = CURL_SOCKET_BAD
const CURL_CSELECT_IN                    = 0x01
const CURL_CSELECT_OUT                   = 0x02
const CURL_CSELECT_ERR                   = 0x04
const CURL_PUSH_OK                       = 0
const CURL_PUSH_DENY                     = 1


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

const standard_confidence_vals = [0.683, 0.955, 0.997]
const standard_colors = [:chartreuse2, :yellow, :red]

include("recipes_stats.jl")
include("recipes_histograms_1d.jl")
include("recipes_histograms_2d.jl")
include("recipes_samples_overview.jl")
include("recipes_prior_overview.jl")
include("split_histograms.jl")
include("localmodes.jl")
include("recipes_ahmi.jl")
include("recipes_prior.jl")
include("recipes_samples_1d.jl")
include("recipes_samples_2d.jl")
include("recipes_diagnostics.jl")

ext: jl


text:
"""

    Implements all three variants of the anomaly detectors based on k-nearest neighbors from
    From outliers to prototypes: Ordering data,
    Stefan Harmelinga and Guido Dornhegea and David Tax and Frank Meinecke and Klaus-Robert Muller, 2005

"""
mutable struct KNNAnomaly{V<:Val}
    t::NNTree
    X::Matrix
    v::V
    k::Int
    tt::Symbol
end

"""
    function KNNAnomaly(k::Int, v::Symbol, tree_type::Symbol = :BruteTree)

Create the k-nn anomaly detector with distance variant v::Symbol, where
v = :kappa - the radius of a ball containing all points
v = :gamma - the average distance to all k-nearest neighbors
v = :delta - the average distance to all k-nearest neighbors
"""
KNNAnomaly(k::Int, v::Symbol, tree_type::Symbol = :BruteTree) =
    KNNAnomaly(eval(tree_type)(Array{Float32,2}(undef, 1, 0)), Array{Float32,2}(undef, 1, 0), Val(v), k,
        tree_type)

"""
    fit!(model::KNNAnomaly, X::Array{T, 2})
"""
function StatsBase.fit!(model::KNNAnomaly, X::Array{T, 2}) where T<:Real
    model.X = X
    model.t = eval(model.tt)(X)
    return nothing
end

"""
    StatsBase.predict(model::KNNAnomaly, x, [k::Int, v::Symbol])

v = :kappa - the radius of a ball containing all points
v = :gamma - the average distance to all k-nearest neighbors
v = :delta - the average distance to all k-nearest neighbors
"""
function StatsBase.predict(model::KNNAnomaly, x, k, v::V) where {V<:Val{:kappa}}
    if size(model.X,2) == 0
        error("kNN model not fitted, call fit! before predict")
    end
    inds, dists = NearestNeighbors.knn(model.t, x, k, true)
    map(d -> d[end],dists)
end
function StatsBase.predict(model::KNNAnomaly, x, k, v::V) where {V<:Val{:gamma}}
    if size(model.X,2) == 0
        error("kNN model not fitted, call fit! before predict")
    end
    inds, dists = NearestNeighbors.knn(model.t, x, k)
    map(d -> Statistics.mean(d),dists)
end
function StatsBase.predict(model::KNNAnomaly, x, k, v::V) where {V<:Val{:delta}}
    if size(model.X,2) == 0
        error("kNN model not fitted, call fit! before predict")
    end
    inds, dists = NearestNeighbors.knn(model.t, x, k)
    map(i -> LinearAlgebra.norm(x[:,i[1]] - Statistics.mean(model.X[:,i[2]],dims=2)) ,enumerate(inds))
end
StatsBase.predict(model::KNNAnomaly, x, k, v::Symbol) = StatsBase.predict(model, x, k, Val(v))
StatsBase.predict(model::KNNAnomaly, x, k) = StatsBase.predict(model, x, k, model.v)
StatsBase.predict(model::KNNAnomaly, x, v::Symbol) = StatsBase.predict(model, x, model.k, v)
StatsBase.predict(model::KNNAnomaly, x) = StatsBase.predict(model, x, model.k, model.v)

knn_constructor(;k::Int=1,v::Symbol=:kappa) = KNNAnomaly(k, v)


text:
export eeInteraction, eeInteractions

@inline function twoBodyBFTensorCore(libcinFunc::Symbol,
                                     bf1::FloatingGTBasisFuncs, bf2::FloatingGTBasisFuncs,
                                     bf3::FloatingGTBasisFuncs, bf4::FloatingGTBasisFuncs;
                                     isGradient::Bool=false)
    env = Float64[]
    atm = Int32[]
    bas = Int32[]
    subSize = basisSize([bf1.subshell, bf2.subshell, bf3.subshell, bf4.subshell])

    id, uniqueBFs = markUnique([bf1, bf2, bf3, bf4])

    for bf in uniqueBFs
        addToDataChain!(env, atm, bas, bf)
    end

    cintFunc!(Val(libcinFunc), (subSize..., 1+isGradient*2)|>zeros,
              id .- 1, atm, 0, bas, length(uniqueBFs), env)
end

"""

    twoBodyBFTensor(libcinFunc::Symbol,
                    b1::AbstractGTBasisFuncs, b2::AbstractGTBasisFuncs,
                    b3::AbstractGTBasisFuncs, b4::AbstractGTBasisFuncs;
                    isGradient::Bool=false) ->
    Array{Float64, 5}

Core function for one-electron integrals.

`libcinFunc::Symbol` specifies the backend [libcint](https://github.com/sunqm/libcint)
function name, e.g. `"cint2e_cart"` should be converted to `:cint2e_cart` as the input
argument.
"""
@inline function twoBodyBFTensor(libcinFunc::Symbol,
                                 b1::AbstractGTBasisFuncs,
                                 b2::AbstractGTBasisFuncs,
                                 b3::AbstractGTBasisFuncs,
                                 b4::AbstractGTBasisFuncs;
                                 isGradient::Bool=false)
    f = @inline function (i,j,k,l)
        ints = twoBodyBFTensorCore(libcinFunc, i, j, k, l; isGradient)
        ints[ijkIndex(i), ijkIndex(j), ijkIndex(k), ijkIndex(l),:]
    end
    sum([f(i,j,k,l) for i in unpackBasisFuncs(b1), j in unpackBasisFuncs(b2),
                        k in unpackBasisFuncs(b3), l in unpackBasisFuncs(b4)])
end

@inline function twoBodyBSTensor(BasisSet::Vector{<:AbstractGTBasisFuncs},
                                 intFunc::F; outputUniqueIndices::Bool=false) where
                                {F<:Function}
    subSize = basisSize(BasisSet) |> collect
    accuSize = vcat(0, accumulate(+, subSize))
    totalSize = subSize |> sum
    nPage = (intFunc(BasisSet[1], BasisSet[1], BasisSet[1], BasisSet[1]) |> size)[5]
    buf = ones(totalSize, totalSize, totalSize, totalSize, nPage)
    for i = 1:length(BasisSet), j = 1:i, k = 1:i, l = 1:(k==i ? j : k)
        I = accuSize[i]+1 : accuSize[i+1]
        J = accuSize[j]+1 : accuSize[j+1]
        K = accuSize[k]+1 : accuSize[k+1]
        L = accuSize[l]+1 : accuSize[l+1]
        subBuf = intFunc(BasisSet[i], BasisSet[j], BasisSet[k], BasisSet[l])
        for page = 1:nPage
            buf[I,J,K,L,page:page] .= subBuf
            buf[J,I,K,L,page:page] .= PermutedDimsArray(subBuf, [2,1,3,4,5])
            buf[J,I,L,K,page:page] .= PermutedDimsArray(subBuf, [2,1,4,3,5])
            buf[I,J,L,K,page:page] .= PermutedDimsArray(subBuf, [1,2,4,3,5])
            buf[L,K,I,J,page:page] .= PermutedDimsArray(subBuf, [4,3,1,2,5])
            buf[K,L,I,J,page:page] .= PermutedDimsArray(subBuf, [3,4,1,2,5])
            buf[K,L,J,I,page:page] .= PermutedDimsArray(subBuf, [3,4,2,1,5])
            buf[L,K,J,I,page:page] .= PermutedDimsArray(subBuf, [4,3,2,1,5])
        end
    end
    if outputUniqueIndices
        s = sum(subSize)
        uniqueInts = fill(Int[0,0,0,0,0],
                          (3*binomial(s, 4)+6*binomial(s, 3)+4*binomial(s, 2)+s)*nPage)
        index = 1
        for i = 1:s, j = 1:i, k = 1:i, l = 1:(k==i ? j : k), p=1:nPage
            uniqueInts[index] = [i, j, k, l, p]
            index += 1
        end
    end
    outputUniqueIndices ? (return buf, uniqueInts) : (return buf)
end

@inline eeInteractionCore(bf1::AbstractGTBasisFuncs,
                          bf2::AbstractGTBasisFuncs,
                          bf3::AbstractGTBasisFuncs,
                          bf4::AbstractGTBasisFuncs) =
twoBodyBFTensor(:cint2e_cart, bf1, bf2, bf3, bf4)

"""

    eeInteraction(bf1::AbstractGTBasisFuncs,
                  bf2::AbstractGTBasisFuncs,
                  bf3::AbstractGTBasisFuncs,
                  bf4::AbstractGTBasisFuncs) ->
    Array{Float64, 4}

Return the electron-electron interaction tensor (an N×N×N×N Tensor where N is the number of
spatial orbitals) given 4 basis functions.
"""
eeInteraction(bf1::AbstractGTBasisFuncs,
              bf2::AbstractGTBasisFuncs,
              bf3::AbstractGTBasisFuncs,
              bf4::AbstractGTBasisFuncs) =
dropdims(eeInteractionCore(bf1, bf2, bf3, bf4), dims=5)

"""

    eeInteractionsCore(BSet::Array{<:AbstractGTBasisFuncs, 1};
                       outputUniqueIndices::Bool=false) ->
    Array{Float64, 5}, [Array{<:Array{Int, 1}, 1}]

Return the electron-electron interaction tensor (an N×N×N×N×1 Tensor where N is the number
of spatial orbitals) given a basis set in the form of an `Array`.

If `outputUniqueIndices=true`, additionally return the indices for all the unique integrals.
"""
@inline eeInteractionsCore(BSet::Vector{<:AbstractGTBasisFuncs};
                          outputUniqueIndices::Bool=false) =
        twoBodyBSTensor(BSet, eeInteractionCore; outputUniqueIndices)

"""

    eeInteractions(BSet::Array{<:AbstractGTBasisFuncs, 1}) -> Array{Float64, 4}

Return the electron-electron interaction tensor (an N×N×N×N Tensor where N is the number
of spatial orbitals) given a basis set in the form of an `Array`.
"""
eeInteractions(BSet::Vector{<:AbstractGTBasisFuncs}) =
dropdims(eeInteractionsCore(BSet), dims=5)


text:
using DifferentialEquations
using Plots
plotlyjs()

# include the case files
include("ex2-case-1.jl")
include("ex2-case-2.jl")
include("ex2-case-3.jl")

#initial conditions
r₀=6
v₀=-0.25
u₀=[17*π/18,0.1,0.1,0.1]

#calling the functions to run them
IC=[r₀,v₀,u₀]
IC2=ex2_case_1(IC)
IC3=ex2_case_2(IC2)
FC=ex2_case_3(IC3);

#overall trajectory
t1=IC2[4]
t2=IC3[4]
t3=FC[4]
t=vcat(t1,t2,t3)
display(plot(t[:,1],t[:,2],t[:,3],label="trajectory",
    xaxis="x (m)",yaxis="y (m)",zaxis="z (m)",grid=(:on,:black)))
savefig("./plots/trajectory.png")

max_stars_repo_path: Modeling the Dragon Module as a Pendulum/Scripts/ex2.jl


text:
module Coexist

using DataFrames
using CSVFiles
using LinearAlgebra
using Dates
import StatsFuns: logistic, gammapdf
using PyCall

const DATA_DIR = joinpath(dirname(@__FILE__), "..", "data")

include("diseaseProg.jl")
include("utils.jl")
include("generated_dydt.jl")
include("pyCoexist.jl")

end # module

max_stars_repo_name: vollmersj/coexist-julia

max_forks_repo_name: vollmersj/coexist-julia


text:
@testset "Resize" begin
    @test (Resize <: Augmentor.AffineOperation) == false
    @test typeof(@inferred(Resize())) <: Resize <: Augmentor.ImageOperation
    @testset "constructor" begin
        @test_throws MethodError Resize(())
        @test_throws MethodError Resize(1.,2.)
        @test_throws MethodError Resize(:a)
        @test_throws MethodError Resize([:a])
        @test_throws ArgumentError Resize(-1)
        @test_throws ArgumentError Resize(0,2)
        op = @inferred(Resize(10))
        @test str_show(op) == "Augmentor.Resize{1}((10,))"
        @test str_showconst(op) == "Resize(10)"
        @test str_showcompact(op) == "Resize to (10,)"
        op = @inferred(Resize(20,30))
        @test op === Resize(width=30, height=20)
        @test op.size == (20,30)
        @test str_show(op) == "Augmentor.Resize{2}((20, 30))"
        @test str_showconst(op) == "Resize(20, 30)"
        @test str_showcompact(op) == "Resize to 20×30"
        op = @inferred(Resize(20,30,40))
        @test op === @inferred(Resize((20,30,40)))
        @test op === @inferred(Resize{3}((20,30,40)))
        @test op.size == (20,30,40)
        @test str_show(op) == "Augmentor.Resize{3}((20, 30, 40))"
        @test str_showconst(op) == "Resize(20, 30, 40)"
        @test str_showcompact(op) == "Resize to (20, 30, 40)"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(Resize(10,10), nothing)
        @test Augmentor.supports_eager(Resize) === true
        ref = Gray{N0f8}[0.624 0.686 0.733 0.686 0.612; 0.667 0.055 0.71 0.675 0.596; 0.639 0.043 0.227 0.631 0.604; 0.569 0.031 0.518 0.553 0.529; 0.392 0.145 0.392 0.443 0.369]
        imgs = [
            (camera),
            (Augmentor.prepareaffine(camera)),
            (OffsetArray(camera, -10, -30)),
            (view(camera, IdentityRange(1:512), IdentityRange(1:512))),
        ]
        @testset "single image" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(Resize(5,5), img_in))
                @test res == ref
                @test typeof(res) == typeof(ref)
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs
                img_in2 = N0f8.(img_in1)
                img_out1 = ref
                img_out2 = N0f8.(ref)
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(Resize(5,5), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(Resize) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Resize) === true
        @test_throws MethodError Augmentor.applyaffineview(Resize(10,10), nothing)
        @test @inferred(Augmentor.toaffinemap(Resize(4,9), rect)) ≈ AffineMap([2. 0.; 0. 3], [-0.5,-1.0])
        for h in (1,2,3,4,5,9), w in (1,2,3,4,5,9)
            wv = @inferred Augmentor.applyaffineview(Resize(h,w), Augmentor.prepareaffine(square))
            @test eltype(wv) == eltype(square)
            @test typeof(wv) <: SubArray
            @test typeof(wv.indices) <: Tuple{Vararg{IdentityRange}}
            @test typeof(parent(wv)) <: InvWarpedView
            @test parent(parent(wv)).itp.coefs === square
            # round because `imresize` computes as float space,
            # while applyaffineview doesn't
            @test round.(Float64.(wv), digits=1) == round.(Float64.(imresize(square, h, w)), digits=1)
        end
        for h in (1,2,3,4,5,9), w in (1,2,3,4,5,9) # bigger show drift
            wv = @inferred Augmentor.applyaffineview(Resize(h,w), Augmentor.prepareaffine(checkers))
            @test eltype(wv) == eltype(checkers)
            @test typeof(wv) <: SubArray
            @test typeof(wv.indices) <: Tuple{Vararg{IdentityRange}}
            @test typeof(parent(wv)) <: InvWarpedView
            @test parent(parent(wv)).itp.coefs === checkers
            @test wv == imresize(checkers, h, w)
        end
        for h in (3,10,29,30,64), w in (3,10,29,30,64)
            (h, w) == (30, 3) && continue # weird tiny artifact in one pixel at []
            wv = @inferred Augmentor.applyaffineview(Resize(h,w), Augmentor.prepareaffine(camera))
            @test eltype(wv) == eltype(camera)
            @test typeof(wv) <: SubArray
            @test typeof(wv.indices) <: Tuple{Vararg{IdentityRange}}
            @test typeof(parent(wv)) <: InvWarpedView
            @test parent(parent(wv)).itp.coefs === camera
            @test wv == imresize(camera, h, w)
        end
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Resize) === true
        @testset "single image" begin
            wv = @inferred Augmentor.applylazy(Resize(2,3), square)
            @test eltype(wv) == eltype(square)
            @test typeof(wv) <: SubArray
            @test typeof(wv.indices) <: Tuple{Vararg{IdentityRange}}
            @test typeof(parent(wv)) <: InvWarpedView
            @test typeof(parent(parent(wv))) <: Interpolations.Extrapolation
            @test parent(parent(wv)).itp.coefs === square
            @test wv == imresize(square, 2, 3)
        end
        @testset "multiple images" begin
            wv1, wv2 = @inferred Augmentor.applylazy(Resize(2,3), (square, square2))
            @test typeof(wv1) <: SubArray
            @test typeof(wv1.indices) <: Tuple{Vararg{IdentityRange}}
            @test typeof(parent(wv1)) <: InvWarpedView
            @test typeof(parent(parent(wv1))) <: Interpolations.Extrapolation
            @test parent(parent(wv1)).itp.coefs === square
            @test wv1 == imresize(square, 2, 3)
            @test typeof(wv2) <: SubArray
            @test typeof(wv2.indices) <: Tuple{Vararg{IdentityRange}}
            @test typeof(parent(wv2)) <: InvWarpedView
            @test typeof(parent(parent(wv2))) <: Interpolations.Extrapolation
            @test parent(parent(wv2)).itp.coefs === square2
            @test wv2 == imresize(square2, 2, 3)
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(Resize) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Resize) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Resize) === false
    end
end


text:
using JSON3, HTTP

url = "https://raw.githubusercontent.com/spdx/license-list-data/v3.10/json/licenses.json"
r = HTTP.request("GET", url)
spdx = JSON3.read(String(r.body))

@assert spdx.licenseListVersion == "3.10"

osi_licenses = sort!([lic.licenseId
                      for lic in spdx.licenses
                      if lic.isOsiApproved])

open(joinpath(@__DIR__, "..", "src", "OSI_LICENSES.jl"); write=true) do io
    header = """
            # THIS FILE IS MACHINE-GENERATED
            # See `../generate_osi_list` for more information."""
    start = "const OSI_LICENSES = Set(["
    spaces = " "^textwidth(start)
    println(io, header)
    println(io, start, "\"", osi_licenses[1], "\",")
    for id in osi_licenses[2:(end - 1)]
        println(io, spaces, "\"", id, "\",")
    end
    return println(io, spaces, "\"", osi_licenses[end], "\"])")
end


text:
#############################################################
# This file was automatically generated on 2021-12-15.      #
#                                                           #
# JuliaPy Bindings Version 2.1.29                           #
#                                                           #
# If you have a bugfix for this file and want to commit it, #
# please fix the bug in the generator. You can find a link  #
# to the generators git repository on tinkerforge.com       #
#############################################################

export get_device_display_name
function get_device_display_name(device_identifier::Integer)
    device_display_names = Dict{Integer, String}(
        11 => "DC Brick",
        13 => "Master Brick",
        14 => "Servo Brick",
        15 => "Stepper Brick",
        16 => "IMU Brick",
        17 => "RED Brick",
        18 => "IMU Brick 2.0",
        19 => "Silent Stepper Brick",
        21 => "Ambient Light Bricklet",
        23 => "Current12 Bricklet",
        24 => "Current25 Bricklet",
        25 => "Distance IR Bricklet",
        26 => "Dual Relay Bricklet",
        27 => "Humidity Bricklet",
        28 => "IO-16 Bricklet",
        29 => "IO-4 Bricklet",
        111 => "HAT Brick",
        112 => "HAT Zero Brick",
        210 => "Joystick Bricklet",
        211 => "LCD 16x2 Bricklet",
        212 => "LCD 20x4 Bricklet",
        213 => "Linear Poti Bricklet",
        214 => "Piezo Buzzer Bricklet",
        215 => "Rotary Poti Bricklet",
        216 => "Temperature Bricklet",
        217 => "Temperature IR Bricklet",
        218 => "Voltage Bricklet",
        219 => "Analog In Bricklet",
        220 => "Analog Out Bricklet",
        221 => "Barometer Bricklet",
        222 => "GPS Bricklet",
        223 => "Industrial Digital In 4 Bricklet",
        224 => "Industrial Digital Out 4 Bricklet",
        225 => "Industrial Quad Relay Bricklet",
        226 => "PTC Bricklet",
        227 => "Voltage/Current Bricklet",
        228 => "Industrial Dual 0-20mA Bricklet",
        229 => "Distance US Bricklet",
        230 => "Dual Button Bricklet",
        231 => "LED Strip Bricklet",
        232 => "Moisture Bricklet",
        233 => "Motion Detector Bricklet",
        234 => "Multi Touch Bricklet",
        235 => "Remote Switch Bricklet",
        236 => "Rotary Encoder Bricklet",
        237 => "Segment Display 4x7 Bricklet",
        238 => "Sound Intensity Bricklet",
        239 => "Tilt Bricklet",
        240 => "Hall Effect Bricklet",
        241 => "Line Bricklet",
        242 => "Piezo Speaker Bricklet",
        243 => "Color Bricklet",
        244 => "Solid State Relay Bricklet",
        246 => "NFC/RFID Bricklet",
        249 => "Industrial Dual Analog In Bricklet",
        250 => "Accelerometer Bricklet",
        251 => "Analog In Bricklet 2.0",
        253 => "Load Cell Bricklet",
        254 => "RS232 Bricklet",
        255 => "Laser Range Finder Bricklet",
        256 => "Analog Out Bricklet 2.0",
        258 => "Industrial Analog Out Bricklet",
        259 => "Ambient Light Bricklet 2.0",
        260 => "Dust Detector Bricklet",
        262 => "CO2 Bricklet",
        263 => "OLED 128x64 Bricklet",
        264 => "OLED 64x48 Bricklet",
        265 => "UV Light Bricklet",
        266 => "Thermocouple Bricklet",
        267 => "Motorized Linear Poti Bricklet",
        268 => "Real-Time Clock Bricklet",
        270 => "CAN Bricklet",
        271 => "RGB LED Bricklet",
        272 => "RGB LED Matrix Bricklet",
        276 => "GPS Bricklet 2.0",
        277 => "RS485 Bricklet",
        278 => "Thermal Imaging Bricklet",
        279 => "XMC1400 Breakout Bricklet",
        282 => "RGB LED Button Bricklet",
        283 => "Humidity Bricklet 2.0",
        284 => "Industrial Dual Relay Bricklet",
        285 => "DMX Bricklet",
        286 => "NFC Bricklet",
        288 => "Outdoor Weather Bricklet",
        289 => "Remote Switch Bricklet 2.0",
        290 => "Sound Pressure Level Bricklet",
        291 => "Temperature IR Bricklet 2.0",
        292 => "Motion Detector Bricklet 2.0",
        293 => "Industrial Counter Bricklet",
        294 => "Rotary Encoder Bricklet 2.0",
        295 => "Analog In Bricklet 3.0",
        296 => "Solid State Relay Bricklet 2.0",
        297 => "Air Quality Bricklet",
        298 => "LCD 128x64 Bricklet",
        299 => "Distance US Bricklet 2.0",
        2100 => "Industrial Digital In 4 Bricklet 2.0",
        2101 => "PTC Bricklet 2.0",
        2102 => "Industrial Quad Relay Bricklet 2.0",
        2103 => "LED Strip Bricklet 2.0",
        2104 => "Load Cell Bricklet 2.0",
        2105 => "Voltage/Current Bricklet 2.0",
        2106 => "Real-Time Clock Bricklet 2.0",
        2107 => "CAN Bricklet 2.0",
        2108 => "RS232 Bricklet 2.0",
        2109 => "Thermocouple Bricklet 2.0",
        2110 => "Particulate Matter Bricklet",
        2111 => "IO-4 Bricklet 2.0",
        2112 => "OLED 128x64 Bricklet 2.0",
        2113 => "Temperature Bricklet 2.0",
        2114 => "IO-16 Bricklet 2.0",
        2115 => "Analog Out Bricklet 3.0",
        2116 => "Industrial Analog Out Bricklet 2.0",
        2117 => "Barometer Bricklet 2.0",
        2118 => "UV Light Bricklet 2.0",
        2119 => "Dual Button Bricklet 2.0",
        2120 => "Industrial Dual 0-20mA Bricklet 2.0",
        2121 => "Industrial Dual Analog In Bricklet 2.0",
        2122 => "Isolator Bricklet",
        2123 => "One Wire Bricklet",
        2124 => "Industrial Digital Out 4 Bricklet 2.0",
        2125 => "Distance IR Bricklet 2.0",
        2127 => "RGB LED Bricklet 2.0",
        2128 => "Color Bricklet 2.0",
        2129 => "Multi Touch Bricklet 2.0",
        2130 => "Accelerometer Bricklet 2.0",
        2131 => "Ambient Light Bricklet 3.0",
        2132 => "Hall Effect Bricklet 2.0",
        2137 => "Segment Display 4x7 Bricklet 2.0",
        2138 => "Joystick Bricklet 2.0",
        2139 => "Linear Poti Bricklet 2.0",
        2140 => "Rotary Poti Bricklet 2.0",
        2144 => "Laser Range Finder Bricklet 2.0",
        2145 => "Piezo Speaker Bricklet 2.0",
        2146 => "E-Paper 296x128 Bricklet",
        2147 => "CO2 Bricklet 2.0",
        2152 => "Energy Monitor Bricklet",
        2153 => "Compass Bricklet",
        2156 => "Performance DC Bricklet",
        2157 => "Servo Bricklet 2.0",
        2161 => "IMU Bricklet 3.0",
        2162 => "Industrial Dual AC Relay Bricklet",
        2164 => "Industrial PTC Bricklet",
        2165 => "DC Bricklet 2.0",
        2166 => "Silent Stepper Bricklet 2.0"
    )

    try
        device_display_name = dict["d"]
    catch e
        if e isa KeyError
            device_display_name = "Unknown Device [{device_identifier}]"
        end
    end

    return device_display_name
end


text:
module AoC2020

# include("util.jl")
include("console.jl")

for day in 1:10
    if day < 10
        include("day0$day.jl")
    else
        include("day$(day).jl")
    end
end

include("day12.jl")
include("day13.jl")
include("day14.jl")

# run all
function runalldays()
    display(Day1.day1())
    display(Day2.day2())
    display(Day3.day3())
    display(Day4.day4())
    display(Day5.day5())
    display(Day6.day6())
    display(Day7.day7())
    display(Day8.day8())
    display(Day9.day9())
    display(Day10.day10())
    # display(Day11.day11())
    display(Day12.day12())
    display(Day13.day13())
    display(Day13.day14())
end

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_SimpleStateEvents

Simulate SimpleStateEvents model.

Demonstrates/tests simple state events
"""
module Simulate_SimpleStateEvents

include(joinpath("models", "SimpleStateEvents.jl"))
import .SimpleStateEvents
import ModiaMath

model  = SimpleStateEvents.Model()
result = ModiaMath.simulate!(model, stopTime=10.0; log=true)
ModiaMath.plot(result, (:s, :v, :f), heading="sSimulate_SimpleStateEvents.jl")

end


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
using EntropyStableEuler
using EntropyStableEuler.Fluxes1D

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end
unorm(U) = sum(map((x->x.^2),U))
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function primitive_to_conservative_hardcode(rho,U,p)
    rhoU = rho.*U
    Unorm = unorm(U)
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rhoU,E)
end

function euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
    rhoR,uR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))

    unorm = (@. uL*uR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. f4aux*uavg)

    return (FxS1,FxS2,FxS3)
end

function euler_fluxes(rhoL,uL,betaL,rhoR,uR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                           rhoR,uR,betaR,rhologR,betalogR)
end

const TOL = 1e-16
"Approximation parameters"
N = 2 # The order of approximation
K = 50#800
T = 1.0#3.0

# Becker viscous shocktube
const γ = 1.4
const M_0 = 20.0
const mu = 0.001
const lambda = 2/3*mu
const Pr = 3/4
const cp = γ/(γ-1)
const cv = 1/(γ-1)
const kappa = mu*cp/Pr

const v_inf = 0.2
const rho_0 = 1.0
const v_0 = 1.0
const m_0 = rho_0*v_0
const v_1 = (γ-1+2/M_0^2)/(γ+1)
const v_01 = sqrt(v_0*v_1)

const uL = v_0+v_inf
const uR = v_1+v_inf
const rhoL = m_0/v_0
const rhoR = m_0/v_1
const eL = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_0^2)
const eR = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_1^2)
const pL = (γ-1)*rhoL*eL
const pR = (γ-1)*rhoR*eR
const EL = pL/(γ-1)+0.5*rhoL*uL^2
const ER = pR/(γ-1)+0.5*rhoR*uR^2

const Bl = -1.0
const Br = 1.5
gr(size=(300,300),ylims=(0,6.0),legend=false,markerstrokewidth=1,markersize=2)
plot()

"Mesh related variables"
VX = LinRange(Bl,Br,K+1)
EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K))

"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
Vq = vandermonde_1D(N,rq)/VDM
M = Vq'*diagm(wq)*Vq
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end
rf = [-1.0;1.0]
nrJ = [-1.0;1.0]
Vf = vandermonde_1D(N,rf)/VDM
LIFT = Mlump_inv*Vf'

# operators
Qr = M*Dr
B = zeros(N+1,N+1)
B[1,1] = -1
B[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Qr0 = S0+1/2*B

# Drop zeros
Qr = Matrix(droptol!(sparse(Qr),TOL))
Qr0 = Matrix(droptol!(sparse(Qr0),TOL))
B = Matrix(droptol!(sparse(B),TOL))
L = Matrix(droptol!(sparse(L),TOL))
S0 = Matrix(droptol!(sparse(S0),TOL))
S = Matrix(droptol!(sparse((Qr-Qr')/2),TOL))

"""High order mesh"""
x = V1*VX[transpose(EToV)]
xf = Vf*x
mapM = reshape(1:2*K,2,K)
mapP = copy(mapM)
mapP[1,2:end] .= mapM[2,1:end-1]
mapP[2,1:end-1] .= mapM[1,2:end]

# """Periodic"""
# mapP[1] = mapM[end]
# mapP[end] = mapP[1]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Initial condition"""
function bisection_solve_velocity(x,max_iter,tol)
    v_L = v_1
    v_R = v_0
    num_iter = 0

    L_k = kappa/m_0/cv
    f(v) = -x+2*L_k/(γ+1)*(v_0/(v_0-v_1)*log((v_0-v)/(v_0-v_01))-v_1/(v_0-v_1)*log((v-v_1)/(v_01-v_1)))

    v_new = (v_L+v_R)/2
    while num_iter < max_iter
        v_new = (v_L+v_R)/2

        if abs(f(v_new)) < tol
            return v_new
        elseif sign(f(v_L)) == sign(f(v_new))
            v_L = v_new
        else
            v_R = v_new
        end
        num_iter += 1
    end

    return v_new
end

const max_iter = 10000
const tol = 1e-14

function exact_sol_viscous_shocktube(x,t)
    u   = bisection_solve_velocity(x-v_inf*t,max_iter,tol)
    rho = m_0/u
    e   = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-u^2)
    return rho, rho*(v_inf+u), rho*(e+1/2*(v_inf+u)^2)
end

U = exact_sol_viscous_shocktube.(x,0.0)
U = ([x[1] for x in U], [x[2] for x in U], [x[3] for x in U])

function zhang_wavespd(rhoL,rhouL,EL,sigma2L,sigma3L,pL,rhoR,rhouR,ER,sigma2R,sigma3R,pR)
    uL   = rhouL/rhoL
    eL   = (EL - .5*rhoL*uL^2)/rhoL
    tauL = sigma2L
    qL   = uL*tauL-sigma3L
    uR   = rhouR/rhoR
    eR   = (ER - .5*rhoR*uR^2)/rhoR
    tauR = sigma2R
    qR   = uR*tauR-sigma3R
    wavespdL = abs(uL)+1/(2*rhoL^2*eL)*(sqrt(rhoL^2*qL^2+2*rhoL^2*eL*abs(tauL-pL)^2)+rhoL*abs(qL))
    wavespdR = abs(uR)+1/(2*rhoR^2*eR)*(sqrt(rhoR^2*qR^2+2*rhoR^2*eR*abs(tauR-pR)^2)+rhoR*abs(qR))
    return max(wavespdL,wavespdR)
end

function zhang_wavespd(rhoL,rhouL,EL,sigma2L,sigma3L,pL)
    uL   = rhouL/rhoL
    eL   = (EL - .5*rhoL*uL^2)/rhoL
    tauL = sigma2L
    qL   = uL*tauL-sigma3L
    wavespdL = abs(uL)+1/(2*rhoL^2*eL)*(sqrt(rhoL^2*qL^2+2*rhoL^2*eL*abs(tauL-pL)^2)+rhoL*abs(qL))
    return wavespdL
end

function limiting_param(U_low, P_ij)
    l = 1.0
    # Limit density
    if U_low[1] + P_ij[1] < -TOL
        l = min(abs(U_low[1])/(abs(P_ij[1])+1e-14), 1.0)
    end

    # limiting internal energy (via quadratic function)
    a = P_ij[1]*P_ij[3]-1.0/2.0*P_ij[2]^2
    b = U_low[3]*P_ij[1]+U_low[1]*P_ij[3]-U_low[2]*P_ij[2]
    c = U_low[3]*U_low[1]-1.0/2.0*U_low[2]^2

    l_eps_ij = 1.0
    if b^2-4*a*c >= 0
        r1 = (-b+sqrt(b^2-4*a*c))/(2*a)
        r2 = (-b-sqrt(b^2-4*a*c))/(2*a)
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

function flux_lowIDP(U_i,U_j,f_i,f_j,c_ij,wavespd)
    return c_ij*(f_i+f_j)-abs(c_ij)*wavespd*(U_j-U_i)
end

function flux_ES(rho_i,u_i,beta_i,rho_j,u_j,beta_j,c_ij)
    return 2*c_ij.*euler_fluxes(rho_i,u_i,beta_i,rho_j,u_j,beta_j)
end

function flux_viscous(sigma_i,sigma_j,c_ij)
    return c_ij*(sigma_i+sigma_j)
end

function rhs_viscous(U,K,N,Mlump_inv,S)
    J = (Br-Bl)/K/2 # assume uniform interval
    Nc = 3
    rhsU  = [zeros(N+1,K) for i = 1:Nc]
    theta = [zeros(N+1,K) for i = 1:Nc]
    sigma = [zeros(N+1,K) for i = 1:Nc]

    p = pfun_nd.(U[1],U[2],U[3])

    VU = v_ufun(U...)
    VUf = (x->Vf*x).(VU)
    VUP = (x->x[mapP]).(VUf)
    VL = v_ufun(rhoL,rhoL*uL,EL)
    VUP[1][1] = VL[1]
    VUP[2][1] = VL[2]
    VUP[3][1] = VL[3]
    VUP[1][end] = VUf[1][end]
    VUP[2][end] = VUf[2][end]
    VUP[3][end] = VUf[3][end]

    surfx(uP,uf) = LIFT*(@. .5*(uP-uf)*nxJ)
    # Strong form, dv/dx
    VUx = (x->Dr*x).(VU)
    VUx = VUx .+ surfx.(VUP,VUf)
    VUx = VUx./J

    # σ = K dv/dx
    for k = 1:K
        for i = 1:N+1
            Kx = zeros(3,3)
            v1 = VU[1][i,k]
            v2 = VU[2][i,k]
            v4 = VU[3][i,k]
            Kx[2,2] = -(2*mu-lambda)/v4
            Kx[2,3] = (2*mu-lambda)*v2/v4^2
            Kx[3,2] = Kx[2,3]
            Kx[3,3] = -(2*mu-lambda)*v2^2/v4^3+kappa/cv/v4^2

            sigma[2][i,k] += Kx[2,2]*VUx[2][i,k] + Kx[2,3]*VUx[3][i,k]
            sigma[3][i,k] += Kx[3,2]*VUx[2][i,k] + Kx[3,3]*VUx[3][i,k]
        end
    end

    sxf = (x->Vf*x).(sigma)
    sxP = (x->x[mapP]).(sxf)
    sxP[1][1] = sxf[1][1]
    sxP[2][1] = sxf[2][1]
    sxP[3][1] = sxf[3][1]
    sxP[1][end] = sxf[1][end]
    sxP[2][end] = sxf[2][end]
    sxP[3][end] = sxf[3][end]
    # strong form, dσ/dx
    penalization(uP,uf) = LIFT*(@. -.5*(uP-uf))
    sigmax = (x->Dr*x).(sigma)
    sigmax = sigmax .+ surfx.(sxP,sxf) #.+ penalization.(VUP,VUf)
    sigmax = sigmax./J

    return sigmax, sigma
end

function rhs_IDP(U,K,N,Mlump_inv,S0,S,wq)
    p = pfun_nd.(U[1],U[2],U[3])
    flux = zero.(U)
    @. flux[1] = U[2]
    @. flux[2] = U[2]^2/U[1]+p
    @. flux[3] = U[3]*U[2]/U[1]+p*U[2]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[1]/(2*p)

    J = (Br-Bl)/K/2 # assume uniform interval

    # Low order and high order algebraic fluxes
    F_low       = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_high      = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_low_P     = [zeros(2),zeros(2),zeros(2)] # 1: left boundary, 2: right boundary
    F_high_P    = [zeros(2),zeros(2),zeros(2)]
    L           = zeros(N+1,N+1) # Array of limiting params
    rhsU        = [zeros(N+1,K),zeros(N+1,K),zeros(N+1,K)]

    # TODO: redundant!
    _,sigma = rhs_viscous(U,K,N,Mlump_inv,S)
    wavespd_arr = zeros(N+1,K)
    beta_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            wavespd_arr[i,k] = wavespeed_1D(U[1][i,k],U[2][i,k],U[3][i,k])
            beta_arr[i,k] = zhang_wavespd(U[1][i,k],U[2][i,k],U[3][i,k],sigma[2][i,k],sigma[3][i,k],p[i,k])
        end
    end

    L_plot = zeros(N+1,K)

    #dt = Inf
    # TODO: hardcoded
    dt = 1e-4
    d_ii_arr = zeros(N+1,K)
    for k = 1:K
        L = zeros(N+1,N+1)
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k],beta_arr[i,k],beta_arr[j,k])
                    fluxS = flux_ES(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],S[i,j])
                    for c = 1:3
                        F_low[c][i,j]  = -flux_viscous(sigma[c][i,k],sigma[c][j,k],S[i,j])+flux_lowIDP(U[c][i,k],U[c][j,k],flux[c][i,k],flux[c][j,k],S0[i,j],wavespd)
                        F_high[c][i,j] = -flux_viscous(sigma[c][i,k],sigma[c][j,k],S[i,j])+fluxS[c]
                    end
                    d_ii_arr[i,k] += wavespd
                end
            end
        end

        EL = pL/(γ-1)+0.5*rhoL*uL^2
        ER = pR/(γ-1)+0.5*rhoR*uR^2
        U_left   = (k == 1) ? [rhoL; rhoL*uL; EL]                 : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
        Ub_left  = (k == 1) ? [rhoL; uL; rhoL/(2*pL)]             : [Ub[1][end,k-1]; Ub[2][end,k-1]; Ub[3][end,k-1]]
        f_left   = (k == 1) ? [rhoL*uL; rhoL*uL^2+pL; uL*(pL+EL)] : [flux[1][end,k-1]; flux[2][end,k-1]; flux[3][end,k-1]]
        U_right  = (k == K) ? [rhoR; rhoR*uR; ER]                 : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
        Ub_right = (k == K) ? [rhoR; uR; rhoR/(2*pR)]             : [Ub[1][1,k+1]; Ub[2][1,k+1]; Ub[3][1,k+1]]
        f_right  = (k == K) ? [rhoR*uR; rhoR*uR^2+pR; uR*(pR+ER)] : [flux[1][1,k+1]; flux[2][1,k+1]; flux[3][1,k+1]]
        p_left   = (k == 1) ? pL : p[end,k-1]
        p_right  = (k == K) ? pR : p[1,k+1]
        # Assume velocity doesn't change at boundary, isothermal
        TL = pL/rhoL/(γ-1)/cv
        TR = pR/rhoR/(γ-1)/cv
        sigma_left  = (k == 1) ? [sigma[1][1,1];sigma[2][1,1];sigma[3][1,1]]       : [sigma[1][end,k-1]; sigma[2][end,k-1]; sigma[3][end,k-1]]
        sigma_right = (k == K) ? [sigma[1][end,k];sigma[2][end,k];sigma[3][end,k]] : [sigma[1][1,k+1]; sigma[2][1,k+1]; sigma[3][1,k+1]]
        wavespd_l = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]),beta_arr[1,k],zhang_wavespd(U_left[1],U_left[2],U_left[3],sigma_left[2],sigma_left[3],p_left))
        wavespd_r = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]),beta_arr[end,k],zhang_wavespd(U_right[1],U_right[2],U_right[3],sigma_right[2],sigma_right[3],p_right))

        d_ii_arr[1,k] += wavespd_l
        d_ii_arr[end,k] += wavespd_r

        for c = 1:3
            F_low_P[c][1] = -flux_viscous(sigma[c][1,k],sigma_left[c],-1/2)+flux_lowIDP(U[c][1,k],U_left[c],flux[c][1,k],f_left[c],-0.5,wavespd_l)
            F_low_P[c][2] = -flux_viscous(sigma[c][end,k],sigma_right[c],1/2)+flux_lowIDP(U[c][end,k],U_right[c],flux[c][end,k],f_right[c],0.5,wavespd_r)

            if k == 1
                F_low_P[c][1] = -flux_viscous(sigma[c][1,k],sigma_left[c],-1/2)-0.5*(flux[c][1,k]+f_left[c])
            end

            if k == K
                F_low_P[c][2] = -flux_viscous(sigma[c][end,k],sigma_right[c],1/2)+0.5*(flux[c][end,k]+f_right[c])
            end

            F_high_P[c][1] = F_low_P[c][1]
            F_high_P[c][2] = F_low_P[c][2]
        end

        P_ij = zeros(3,1)
        # TODO: redundant
        U_low =  [zeros(N+1,1),zeros(N+1,1),zeros(N+1,1)]
        for c = 1:3
            U_low[c] .= sum(-F_low[c],dims=2)
            U_low[c][1] -= F_low_P[c][1]
            U_low[c][end] -= F_low_P[c][2]
            U_low[c] .= U[c][:,k]+dt/J*Mlump_inv*U_low[c]
        end

        for i = 1:N+1
            lambda_j = 1/N
            m_i = J*wq[i]
            for j = 1:N+1
                if i != j
                    for c = 1:3
                        P_ij[c] = dt/(m_i*lambda_j)*(F_low[c][i,j]-F_high[c][i,j])
                    end
                    L[i,j] = limiting_param([U_low[1][i]; U_low[2][i]; U_low[3][i]],P_ij)
                end
            end
        end

        # Symmetrize limiting parameters
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    l_ij = min(L[i,j],L[j,i])
                    L[i,j] = l_ij
                    L[j,i] = l_ij
                end
            end
        end

        # elementwise limiting
        l = 1.0
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    l = min(l,L[i,j])
                end
            end
        end
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    L[i,j] = l
                    L[j,i] = l
                end
            end
        end

        L = ones(size(L))

        for c = 1:3
            # With limiting
            rhsU[c][:,k] = sum((L.-1).*F_low[c] - L.*F_high[c],dims=2)

            rhsU[c][1,k] += -F_high_P[c][1]
            rhsU[c][N+1,k] += -F_high_P[c][2]

            rhsU[c][:,k] .= 1/J*Mlump_inv*rhsU[c][:,k]
        end

        for i = 1:N+1
            L_plot[i,k] = l
        end
    end

    dt = min(dt,minimum(J*wq./d_ii_arr./2.0))

    return rhsU,dt,L_plot
end

function rhs_IDP(U,K,N,Mlump_inv,S0,S,wq,dt)
    p = pfun_nd.(U[1],U[2],U[3])
    flux = zero.(U)
    @. flux[1] = U[2]
    @. flux[2] = U[2]^2/U[1]+p
    @. flux[3] = U[3]*U[2]/U[1]+p*U[2]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[1]/(2*p)

    J = (Br-Bl)/K/2 # assume uniform interval

    # Low order and high order algebraic fluxes
    F_low       = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_high      = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_low_P     = [zeros(2),zeros(2),zeros(2)] # 1: left boundary, 2: right boundary
    F_high_P    = [zeros(2),zeros(2),zeros(2)]
    L           = zeros(N+1,N+1) # Array of limiting params
    rhsU        = [zeros(N+1,K),zeros(N+1,K),zeros(N+1,K)]

    # TODO: redundant!
    _,sigma = rhs_viscous(U,K,N,Mlump_inv,S)
    wavespd_arr = zeros(N+1,K)
    beta_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            wavespd_arr[i,k] = wavespeed_1D(U[1][i,k],U[2][i,k],U[3][i,k])
            beta_arr[i,k] = zhang_wavespd(U[1][i,k],U[2][i,k],U[3][i,k],sigma[2][i,k],sigma[3][i,k],p[i,k])
        end
    end

    L_plot = zeros(N+1,K)

    for k = 1:K
        L = zeros(N+1,N+1)
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k],beta_arr[i,k],beta_arr[j,k])
                    fluxS = flux_ES(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],S[i,j])
                    for c = 1:3
                        F_low[c][i,j]  = -flux_viscous(sigma[c][i,k],sigma[c][j,k],S[i,j])+flux_lowIDP(U[c][i,k],U[c][j,k],flux[c][i,k],flux[c][j,k],S0[i,j],wavespd)
                        F_high[c][i,j] = -flux_viscous(sigma[c][i,k],sigma[c][j,k],S[i,j])+fluxS[c]
                    end
                end
            end
        end

        U_left   = (k == 1) ? [rhoL; rhoL*uL; EL]                 : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
        Ub_left  = (k == 1) ? [rhoL; uL; rhoL/(2*pL)]             : [Ub[1][end,k-1]; Ub[2][end,k-1]; Ub[3][end,k-1]]
        f_left   = (k == 1) ? [rhoL*uL; rhoL*uL^2+pL; uL*(pL+EL)] : [flux[1][end,k-1]; flux[2][end,k-1]; flux[3][end,k-1]]
        U_right  = (k == K) ? [rhoR; rhoR*uR; ER]                 : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
        Ub_right = (k == K) ? [rhoR; uR; rhoR/(2*pR)]             : [Ub[1][1,k+1]; Ub[2][1,k+1]; Ub[3][1,k+1]]
        f_right  = (k == K) ? [rhoR*uR; rhoR*uR^2+pR; uR*(pR+ER)] : [flux[1][1,k+1]; flux[2][1,k+1]; flux[3][1,k+1]]
        p_left   = (k == 1) ? pL : p[end,k-1]
        p_right  = (k == K) ? pR : p[1,k+1]
        # Assume velocity doesn't change at boundary, isothermal
        TL = pL/rhoL/(γ-1)/cv
        TR = pR/rhoR/(γ-1)/cv
        sigma_left  = (k == 1) ? [sigma[1][1,1];sigma[2][1,1];sigma[3][1,1]]       : [sigma[1][end,k-1]; sigma[2][end,k-1]; sigma[3][end,k-1]]
        sigma_right = (k == K) ? [sigma[1][end,k];sigma[2][end,k];sigma[3][end,k]] : [sigma[1][1,k+1]; sigma[2][1,k+1]; sigma[3][1,k+1]]
        wavespd_l = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]),beta_arr[1,k],zhang_wavespd(U_left[1],U_left[2],U_left[3],sigma_left[2],sigma_left[3],p_left))
        wavespd_r = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]),beta_arr[end,k],zhang_wavespd(U_right[1],U_right[2],U_right[3],sigma_right[2],sigma_right[3],p_right))

        for c = 1:3
            F_low_P[c][1] = -flux_viscous(sigma[c][1,k],sigma_left[c],-1/2)+flux_lowIDP(U[c][1,k],U_left[c],flux[c][1,k],f_left[c],-0.5,wavespd_l)
            F_low_P[c][2] = -flux_viscous(sigma[c][end,k],sigma_right[c],1/2)+flux_lowIDP(U[c][end,k],U_right[c],flux[c][end,k],f_right[c],0.5,wavespd_r)

            F_high_P[c][1] = F_low_P[c][1]
            F_high_P[c][2] = F_low_P[c][2]
        end

        P_ij = zeros(3,1)
        # TODO: redundant
        U_low =  [zeros(N+1,1),zeros(N+1,1),zeros(N+1,1)]
        for c = 1:3
            U_low[c] .= sum(-F_low[c],dims=2)
            U_low[c][1] -= F_low_P[c][1]
            U_low[c][end] -= F_low_P[c][2]
            U_low[c] .= U[c][:,k]+dt/J*Mlump_inv*U_low[c]
        end

        for i = 1:N+1
            lambda_j = 1/N
            m_i = J*wq[i]
            for j = 1:N+1
                if i != j
                    for c = 1:3
                        P_ij[c] = dt/(m_i*lambda_j)*(F_low[c][i,j]-F_high[c][i,j])
                    end
                    L[i,j] = limiting_param([U_low[1][i]; U_low[2][i]; U_low[3][i]],P_ij)
                end
            end
        end

        # Symmetrize limiting parameters
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    l_ij = min(L[i,j],L[j,i])
                    L[i,j] = l_ij
                    L[j,i] = l_ij
                end
            end
        end

        # elementwise limiting
        l = 1.0
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    l = min(l,L[i,j])
                end
            end
        end
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    L[i,j] = l
                    L[j,i] = l
                end
            end
        end

        # L = ones(size(L))

        for c = 1:3
            # With limiting
            rhsU[c][:,k] = sum((L.-1).*F_low[c] - L.*F_high[c],dims=2)

            rhsU[c][1,k] += -F_high_P[c][1]
            rhsU[c][N+1,k] += -F_high_P[c][2]

            rhsU[c][:,k] .= 1/J*Mlump_inv*rhsU[c][:,k]
        end

        for i = 1:N+1
            L_plot[i,k] = l
        end
    end

    return rhsU,L_plot
end

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resU = [zeros(size(x)),zeros(size(x)),zeros(size(x))]
resW = [zeros(size(x)),zeros(size(x)),zeros(size(x))]
resZ = [zeros(size(x)),zeros(size(x)),zeros(size(x))]

Vp = vandermonde_1D(N,LinRange(-1,1,10))/VDM

anim = Animation()
const ptL = Bl+(Br-Bl)/K/(N+1)/2
const ptR = Br-(Br-Bl)/K/(N+1)/2
const hplot = (Br-Bl)/K/(N+1)
i = 1

while t < T

    # SSPRK(3,3)
    # rhsU,dt,_ = rhs_IDP(U,K,N,Mlump_inv,S0,S,wq)
    # dt = min(dt,T-t)
    dt = min(1e-5,T)
    rhsU,_ = rhs_IDP(U,K,N,Mlump_inv,S0,S,wq,dt)
    @. resW = U + dt*rhsU
    rhsU,_ = rhs_IDP(resW,K,N,Mlump_inv,S0,S,wq,dt)
    @. resZ = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resZ
    rhsU,L_plot = rhs_IDP(resW,K,N,Mlump_inv,S0,S,wq,dt)
    @. resZ = resW+dt*rhsU
    @. U = 1/3*U+2/3*resZ

    global t = t + dt
    global i = i + 1
    println("Current time $t with time step size $dt, and final time $T")
    if mod(i,1000) == 1
        plot(x,U[1])
        # for k = 1:K
        #     plot!(ptL+(k-1)*hplot*(N+1):hplot:ptL+k*hplot*(N+1), 1 .-L_plot[:,k],st=:bar,alpha=0.2)
        # end
        frame(anim)
    end
end

plot(x,U[1])
gif(anim,"~/Desktop/tmp.gif",fps=15)

exact_U = @. exact_sol_viscous_shocktube.(x,T)
exact_rho = [x[1] for x in exact_U]
exact_u = [x[2] for x in exact_U]
exact_u = exact_u./exact_rho
exact_E = [x[3] for x in exact_U]

rho = U[1]
u = U[2]./U[1]
E = U[3]
p = pfun_nd.(U[1],U[2],U[3])
J = (Br-Bl)/K/2

Linferr = maximum(abs.(exact_rho-rho))/maximum(abs.(exact_rho)) +
          maximum(abs.(exact_u-u))/maximum(abs.(exact_u)) +
          maximum(abs.(exact_E-E))/maximum(abs.(exact_E))

L1err = sum(J*wq.*abs.(exact_rho-rho))/sum(J*wq.*abs.(rho)) +
        sum(J*wq.*abs.(exact_u-u))/sum(J*wq.*abs.(u)) +
        sum(J*wq.*abs.(exact_E-E))/sum(J*wq.*abs.(E))
println("N = $N, K = $K")
println("L1 error is $L1err")
println("Linf error is $Linferr")

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
@defcomp grosseconomy begin
    K       = Variable(index=[time])    #Capital stock (trillions 2005 US dollars)
    YGROSS  = Variable(index=[time])    #Gross world product GROSS of abatement and damages (trillions 2005 USD per year)

    al      = Parameter(index=[time])   #Level of total factor productivity
    I       = Parameter(index=[time])   #Investment (trillions 2005 USD per year)
    l       = Parameter(index=[time])   #Level of population and labor
    dk      = Parameter()               #Depreciation rate on capital (per year)
    gama    = Parameter()               #Capital elasticity in production function
    k0      = Parameter()               #Initial capital value (trill 2005 USD)

    function run_timestep(p, v, d, t)
        #Define function for K
        if is_first(t)
            v.K[t] = p.k0
        else
            v.K[t] = (1 - p.dk)^5 * v.K[t-1] + 5 * p.I[t-1]
        end

        #Define function for YGROSS
        v.YGROSS[t] = (p.al[t] * (p.l[t]/1000)^(1-p.gama)) * (v.K[t]^p.gama)
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Test

# =========================================================================================
# tess_pipeline_unlv(utf8 = true) => Filename
@testset "tess_pipeline_unlv(utf8 = true) => Filename" begin
    local inst       = TessInst("eng", datadir)
    local pipeline   = TessPipeline(inst)
    local output, io = mktemp(;cleanup=false); close(io)

    @test tess_pipeline_unlv(pipeline, output, true) == true

    @test (tess_run_pipeline(pipeline) do add
        @test add(pix_with([Line("This is image one. ®")]), 72) == true
        @test add(pix_with([Line("This is image two. ®")]), 72) == true
        @test add(pix_with([Line("This is image three. ®")]), 72) == true
        @test add(pix_with([Line("This is image four. ®")]), 72) == true
    end) == true

    @test isfile(output) == true
    @test filesize(output) > 0

    local text  = read(output, String)
    local lines = split(text, "\n"; keepempty=false)
    @test length(lines) == 4

    @test lines[1] == "This is image one. ®"
    @test lines[2] == "This is image two. ®"
    @test lines[3] == "This is image three. ®"
    @test lines[4] == "This is image four. ®"

    rm(output)
end

# =========================================================================================
# tess_pipeline_unlv(utf8 = true) => Filename
@testset "tess_pipeline_unlv(utf8 = false) => Filename" begin
    local inst       = TessInst("eng", datadir)
    local pipeline   = TessPipeline(inst)
    local output, io = mktemp(;cleanup=false); close(io)

    @test tess_pipeline_unlv(pipeline, output, false) == true

    @test (tess_run_pipeline(pipeline) do add
        @test add(pix_with([Line("This is image one. ®")]), 72) == true
        @test add(pix_with([Line("This is image two. ®")]), 72) == true
        @test add(pix_with([Line("This is image three. ®")]), 72) == true
        @test add(pix_with([Line("This is image four. ®")]), 72) == true
    end) == true

    @test isfile(output) == true
    @test filesize(output) > 0

    local text  = read(output, String)
    local lines = split(text, "\n"; keepempty=false)
    @test length(lines) == 4

    @test lines[1] == "This is image one. \xae"
    @test lines[2] == "This is image two. \xae"
    @test lines[3] == "This is image three. \xae"
    @test lines[4] == "This is image four. \xae"

    rm(output)
end

# =========================================================================================
# tess_pipeline_unlv() => String
@testset "tess_pipeline_unlv() => String" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local output   = tess_pipeline_unlv(pipeline)

    @test is_available(output) == false

    @test (tess_run_pipeline(pipeline) do add
        @test add(pix_with([Line("This is image one. ®")]), 72) == true
        @test add(pix_with([Line("This is image two. ®")]), 72) == true
        @test add(pix_with([Line("This is image three. ®")]), 72) == true
        @test add(pix_with([Line("This is image four. ®")]), 72) == true
    end) == true

    @test is_available(output) == true
    @test output[] != nothing

    local lines = split(output[], "\n"; keepempty=false)
    @test length(lines) == 4

    @test lines[1] == "This is image one. ®"
    @test lines[2] == "This is image two. ®"
    @test lines[3] == "This is image three. ®"
    @test lines[4] == "This is image four. ®"
end

# =========================================================================================
# tess_pipeline_text() => Dispatch
@testset "tess_pipeline_text() => Dispatch" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local lines    = Vector{String}()

    tess_pipeline_unlv(pipeline) do line
        local clean = strip(line)
        if isempty(clean) == false
            push!(lines, clean)
        end
    end

    @test (tess_run_pipeline(pipeline) do add
        @test add(pix_with([Line("This is image one. ®")]), 72) == true
        @test add(pix_with([Line("This is image two. ®")]), 72) == true
        @test add(pix_with([Line("This is image three. ®")]), 72) == true
        @test add(pix_with([Line("This is image four. ®")]), 72) == true
    end) == true

    @test length(lines) == 4

    @test lines[1] == "This is image one. ®"
    @test lines[2] == "This is image two. ®"
    @test lines[3] == "This is image three. ®"
    @test lines[4] == "This is image four. ®"
end

# =========================================================================================
# tess_pipeline_unlv_latin1() => Byte Array
@testset "tess_pipeline_unlvlatin1() => Byte Array" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local output   = tess_pipeline_unlv_latin1(pipeline)

    @test is_available(output) == false

    @test (tess_run_pipeline(pipeline) do add
        @test add(pix_with([Line("This is image one. ®")]), 72) == true
        @test add(pix_with([Line("This is image two. ®")]), 72) == true
        @test add(pix_with([Line("This is image three. ®")]), 72) == true
        @test add(pix_with([Line("This is image four. ®")]), 72) == true
    end) == true

    @test is_available(output) == true
    @test output[] != nothing

    local lines = split(String(output[]), "\n"; keepempty=false)
    @test length(lines) == 4

    @test lines[1] == "This is image one. \xae"
    @test lines[2] == "This is image two. \xae"
    @test lines[3] == "This is image three. \xae"
    @test lines[4] == "This is image four. \xae"
end


text:
export ParquetFile,
    write_parquet

import Parquet
import Tables: schema, Schema

struct ParquetFileIter{I} <: AbstractGorIter{I}
    path
    rows::I
end

#Base.eltype(::Type{Parquet.RecordCursor{T}}) where {T} = T
Tables.schema(p::Parquet.RecordCursor{T}) where {T} = Tables.Schema(T)

"""
    ParquetFile(path)

Open Parquet file at `path` as genome ordered stream.

This uses implementation in Parquet.jl, which is not very mature yet.
Iterator struct contains state, so best to use
`data = () -> ParquetFile("data.parquet")`.
"""
function ParquetFile(path)

    p = Parquet.File(path)
    rows = Parquet.RecordCursor(p)
    ParquetFileIter{typeof(rows)}(path, rows)

end

# TODO: change this so iterate creates RecordCursor

"""
    write_parquet(rows, path)
    rows |> write_parquet(path)

Write genome ordered stream `rows` as Parquet file to `path`.
"""
write_parquet(rows, path) = Parquet.write_parquet(path, rows)
write_parquet(path) = rows -> GOR.write_parquet(rows, path)


text:
using DFControl, Test
import DFControl: searchdir, QEFlagInfo, QEControlBlockInfo, QEInputInfo, QEDataBlockInfo, qe_flaginfo, QEInputInfos

@test eltype(qe_flaginfo(Exec("projwfc.x"), :calculation)) == Nothing
@test eltype(qe_flaginfo(Exec("pw.x"), :calculation)) == eltype(qe_flaginfo(:calculation))


text:
@sjdoc StringInterpolation """
 "string1  \$a string2" interpolates the value of `a` into the string.
"""

set_sysattributes(:StringInterpolation,:Protected)
function apprules(mx::Mxpr{:StringInterpolation})
    symstring(margs(mx)...)
end

### StringLength

@sjdoc StringLength """
    StringLength(s)

returns the length of the string `s`.
"""
apprules(mx::Mxpr{:StringLength}) = length(mx[1])

### ToString

@sjdoc ToString """
    ToString(expr)

return a string of the printed form or `expr`.
"""
apprules(mx::Mxpr{:ToString}) = symstring(mx[1])

### StringJoin

@sjdoc StringJoin """
    StringJoin(s1,s2,...)
    StringJoin([s1,s2,...])

return the concatenated strings.
"""
@mkapprule StringJoin
@doap StringJoin(s::ListT) = join(margs(s))
@doap StringJoin(args::String...) = join(args)

### Characters

@sjdoc Characters """
    Characters(s)

return a `List` of the characters in `String` `s`.
"""
@mkapprule Characters
@doap Characters(s::String) = mxpr(:List, [string(x) for x in s]...)

### StringReverse

@sjdoc StringReverse """
    StringReverse(s)

return string `s` with characters in reverse order.
"""

@mkapprule StringReverse

@doap StringReverse(s::String) = reverse(s)

### ToUpperCase

@mkapprule ToUpperCase :nargs => 1
@doap ToUpperCase(s::String) = uppercase(s)

### ToLowerCase

@mkapprule ToLowerCase :nargs => 1
@doap ToLowerCase(s::String) = lowercase(s)

### UpperCaseQ

@mkapprule UpperCaseQ :nargs => 1
@doap UpperCaseQ(s::String) = all(isuppercase, s)

### LowerCaseQ

@mkapprule LowerCaseQ :nargs => 1
@doap LowerCaseQ(s::String) = all(islowercase, s)

### ToCharacterCode

@mkapprule ToCharacterCode :nargs => 1
@doap ToCharacterCode(s::Char) = Int(s)
@doap ToCharacterCode(s::String) = MListA(MxprArgT[Int(c) for c in s])

### FromCharacterCode

@mkapprule FromCharacterCode
@doap FromCharacterCode(i::Int) = Char(i)

## This is where we need pattern matching for instance.
@doap function FromCharacterCode(s::ListT)
    listofpredq(s,integerq) ? _fromcharactercode_list(s) :
       listoflistsofpredq(s,integerq) ? MListA(MxprArgT[ _fromcharactercode_list(x) for x in s]) : mx
end

_fromcharactercode_list(s) = String([Char(i) for i in s])

### LetterNumber

@mkapprule LetterNumber :nargs => 1

@doap LetterNumber(x::String) = length(x) > 1 ? mx : letternumber(x[1])

function letternumber(c::Char)
    Int(uppercase(c)) - 64
end

### DigitQ

@mkapprule DigitQ :nargs =>1
@doap DigitQ(s::Union{Char,AbstractString}) = isdigit(s)

### StringSplit

@mkapprule StringSplit
@doap StringSplit(s::String) = MList(split(s)...)
@doap StringSplit(s::String, d::Union{String,Regex}) = MList(split(s,d)...)

@mkapprule MathematicaFullFormString  :nargs => 1

@doap MathematicaFullFormString(x) = SymataIO.symata_to_mma_fullform_string(x)


text:
module TestTree2DFDSBP

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_fdsbp")

@testset "Linear scalar advection" begin
  @trixi_testset "elixir_advection_extended.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [2.29497319e-06],
      linf = [8.48027655e-06],
      rtol = 5.0e-7) # These results change a little bit and depend on the CI system
  end
end

end # module


text:
function mean_markers_temperature(r::Robot)
    sum_temperatures, num_markers = get_sum_temperatures(r)
    return sum_temperatures/num_markers
end

function get_sum_temperatures(r::Robot)
    side=Ost
    sum_temperatures, num_markers = get_sum_temperatures(r,side)
    while isborder(r,Nord)==false
        move!(r,Nord)
        side=inverse(side)
        part_sum, part_num = get_sum_temperatures(r,side)
        sum_temperatures+=part_sum
        num_markers+=part_num
    end
    return sum_temperatures, num_markers
end

function get_sum_temperatures(r::Robot,side::HorizonSide)

    init_values(r) = ismarkers(r) ? (0,0) : (temperanure(r), 1)
    sum_temperatures, num_markers = init_values(r)

    while isborder(r,side)==false
        if ismarker(r)==true
            sum_temperatures += temperature(r)
            num_markers += 1
        end
    end
    return sum_temperatures, num_markers
end


text:
function generateweeklycomomentum(;
    minprice=PARAM[:plminprice],
    minmarketcap=PARAM[:plminmarketcap],
    minnysequantile=PARAM[:plminnysequantile],
    Fplmomentum::Symbol=PARAM[:plmomentum],
    momentum = PARAM[:momentumstrategies][Fplmomentum],
    lpfilterperiods::Vector{<:Union{Day, Nothing}} = PARAM[:analysislpfilterperiods],
    Fplret::Symbol=PARAM[:plret])

  #this check helps sync up the measure output frequency (weekly/monthly) with that of
  #the comomentum
  PARAM[:crspfrequencysuffix] == "w" || throw("Weekly comomentum should only be
    created when crspfrequencysuffix=\"w\"")

  #housekeeping

  selectsummarydata_pl(df; kwargs...) = selectsummarydata(df; minprice, minmarketcap, minnysequantile, kwargs...)

  crspw::DataFrame = makeweekly(selectsummarydata_pl,
    binprefix="wpl_",
    mainpanel=false,
    computeivol=false)
  validatedates(crspw.date |> unique |> sort!, frequency=:week)
  @assert issorted(crspw, [:permno, :date])
  conditionmomentumstrategies!(crspw)

  crspw = crspw[completecases(crspw, Fplret), :]

  #first do the weekly data
  comomw = generatecomomentum(deepcopy(crspw))
  simpledif(v) = [missing; v[2:end] .- v[1:(end-1)]]

  #comomw = lpfilter(comomw;
  #  focalcols = [:comom10, :comom1, :comom, :lGcomom10, :lGcomom1, :lGcomom], lpfilterperiods)
  comomw.Gcomom =[missing; comomw.comom[2:end] .- comomw.comom[1:(end-1)];]
  return comomw
end

function generatemonthlycomomentum(;
    minprice=PARAM[:plminprice],
    minmarketcap=PARAM[:plminmarketcap],
    minnysequantile=PARAM[:plminnysequantile],
    Fplmomentum::Symbol=PARAM[:plmomentum],
    momentum = PARAM[:momentumstrategies][Fplmomentum],
    lpfilterperiods::Vector{<:Union{Day, Nothing}} = PARAM[:analysislpfilterperiods],
    patchmonthly::Bool = PARAM[:plpatchmonthly],
    Fplret::Symbol=PARAM[:plret])

  PARAM[:crspfrequencysuffix] == "m" || throw("Monthly comomentum should only be
    created when crspfrequencysuffix=\"m\"")

  #housekeeping
  selectsummarydata_pl(df; kwargs...) = selectsummarydata(
    df; minprice, minmarketcap, minnysequantile, kwargs...)

  crspw::DataFrame = makeweekly(selectsummarydata_pl,
    binprefix="wpl_",
    mainpanel=false,
    computeivol=false)
  validatedates(crspw.date |> unique |> sort!, frequency=:week)
  @assert issorted(crspw, [:permno, :date])

  crspw = crspw[completecases(crspw, Fplret), :]

  if patchmonthly #this implies using the monthly momentum
    (PARAM[:refreshweekly] ≠ PARAM[:refreshmonthly]) && @warn(
      "PARAM[:refreshweekly]=$(PARAM[:refreshweekly]) while
      PARAM[:refreshmonthly]=$(PARAM[:refreshmonthly])!!! This could imply
      out of sync crsp files as used to compute comomentum.")
    crspm::DataFrame = makemonthly(selectsummarydata_pl,
      binprefix="mpl_",
      mainpanel=false,
      computeivol=false)
    validatedates(crspm.date |> unique |> sort!, frequency=:month)
    @assert issorted(crspm, [:permno, :date])

    crspm = crspm[completecases(crspm, Fplret), :]

    #for the monthly data, we first need to splice in the monthly deciles into the weekly data
    crspw.yearmonth = year.(crspw.date) .+ month.(crspw.date) ./ 100
    crspm.yearmonth = year.(crspm.date) .+ month.(crspm.date) ./ 100
    Nrowold = nrow(crspw)
    select!(crspw, Not(Fplmomentum))
    conditionmomentumstrategies!(crspm)
    crspw = leftjoin(crspw, crspm[!, [:permno, :yearmonth, Fplmomentum]], on=[:permno, :yearmonth],
      validate=(false, true))
    if !issorted(crspw, [:permno, :yearmonth])
      sort!(crspw, [:permno, :yearmonth])
    end
    @assert Nrowold == nrow(crspw)
  else
    conditionmomentumstrategies!(crspw)
  end

  comomm = generatecomomentum(crspw)

  #drop the extraneous rows
  @assert issorted(comomm, :date)
  @assert allunique(comomm.date)
  comomm.yearmonth = year.(comomm.date) .+ month.(comomm.date) ./ 100
  comomm = combine(last, groupby(comomm, :yearmonth), renamecols=false)
  select!(comomm, Not(:yearmonth))

  #run through an LP filter
  #comomm = lpfilter(comomm;
  #  focalcols = [:comom10, :comom1, :comom, :lGcomom10, :lGcomom1, :lGcomom], lpfilterperiods)
  comomm.Gcomom =[missing; comomm.comom[2:end] .- comomm.comom[1:(end-1)];]

  return comomm
end

function generatecomomentum(crsp::DataFrame;
  testdate::Date = PARAM[:testdate],
  minwindownonmissing = PARAM[:plminwindownonmissing],
  Fplret::Symbol=PARAM[:plret],
  Fplmomentum::Symbol=PARAM[:plmomentum]
  )

  crspprefix = Dict("m"=>PARAM[:crspmonthly], "w"=>PARAM[:crspweekly])[PARAM[:crspfrequencysuffix]]

  #merge the monthly momentum deciles with the weekly data

  #merge in ff
  crsp = mergeff(crsp, datefrequency=:week, Fgroup=:permno)
  if !issorted(crsp, [:permno, :date])
    sort!(crsp, [:permno, :date])
  end

  #need a symbol id for each column
  crsp.spermno = crsp.permno .|> (i)->Symbol(:p, i)

  #these are conditioned on when computing the correlations
  crsp.some1s = ones(nrow(crsp))
  Fcontrols = [:F_ff3m_mktrf, :F_ff3m_smb, :F_ff3m_hml]
  Fwindowcols = [:date; :spermno; Fplret; Fcontrols]
  crsp = crsp[completecases(crsp, Fwindowcols), :]

  #compute number missing in past year
  #this will be useful when we compute the cross-sectional correlations
  crsp.plretnmissing1y = missings(Float64, nrow(crsp))
  trailingmeasure!(length, crsp, Year(1), 1,
    Ftarget=Fplret,
    Fmeasure=:plretnmissing1y,)

  function deciles(v::AbstractVector{Float64})::Vector{Int} where T
    q = ecdf(v)
    #the extra min is due to the top value only
    d = (q.(v) .* 10 .+ 1) .|> x->min(floor(Int, x),10)
    @assert (1,10) ≡ extrema(d) "extrema(d): $(extrema(d))"
    return d
  end

  #create a version for iterating on deciles
  crspdeciles = crsp[completecases(crsp, [Fplmomentum]), [:date, :spermno, Fplmomentum, :plretnmissing1y, :mc]]
  crspdeciles.plmomentum = Vector{Float64}(crspdeciles[!, Fplmomentum])
  transform!(groupby(crspdeciles, :date),
    :plmomentum => deciles => :d_plmomentum,
    :plretnmissing1y => maximum => :maxnonmissing)

  #this is useful to examine the completeness
  completeness = combine(groupby(crspdeciles, :date),
    nrow => :totalpermno,
    [:plretnmissing1y,:maxnonmissing] =>
      ((plretnmissing1y,maxnonmissing) ->
        sum(plretnmissing1y .≥ minwindownonmissing * maxnonmissing)) =>
      :completepermno,
    :maxnonmissing => last => :maxnonmissing,
    [:plretnmissing1y,:maxnonmissing, :mc] =>
      ((plretnmissing1y,maxnonmissing, mc) ->
        sum((plretnmissing1y .≥ minwindownonmissing * maxnonmissing) .* mc)) =>
      :completemc,
    :mc => sum => :totalmc)
  completeness.perkept = completeness.completepermno ./ completeness.totalpermno
  completeness.permckept = completeness.completemc ./ completeness.totalmc
  #println("complete by date\n", completeness)

  @info "Average comomentum $crspprefix completeness: $(mean(completeness.perkept))
    Value-weighted: $(mean(completeness.permckept))"

  #drop unused deciles
  crspdeciles = crspdeciles[(crspdeciles.d_plmomentum .== 1) .| (crspdeciles.d_plmomentum .== 10), :]
  crspdeciles = crspdeciles[
    crspdeciles.plretnmissing1y .≥ minwindownonmissing .* crspdeciles.maxnonmissing, :]
  select!(crspdeciles, Not([:maxnonmissing, :mc]))

  #now compute the cross-correlations
  crspdeciles.xcor = missings(Float64, nrow(crspdeciles))

  #force controls to exclude missing values
  controls = sort!(unique(crsp[!, [:date; Fcontrols;]]), :date)
  for Fcontrol ∈ Fcontrols
    controls[!, Fcontrol] = Vector{Float64}(controls[!, Fcontrol])
  end

  #calls the partial correlation function after accounting for missing
  function mpartialcor(x::AbstractVector{<:MFloat64},
      y::AbstractVector{<:MFloat64},
      Z::AbstractMatrix{Float64})
    pcorinds = (x .!== missing) .& (y .!== missing)

    return partialcor(x[pcorinds], y[pcorinds], Z[pcorinds, :])
  end

  #Approach below is as follows:
  #1) For each date, select a lagging 1 year window
  #2) For each of the relevant deciles, identify the permnos with complete enough series over the window
  #3) Unstack these into a TxN matrix. Compute the EW portoflio returns after accounting for missings
  #4) Compute the partial correlation between each column and the EW portfolio returns

  gcrspdeciles = groupby(crspdeciles, [:date])
  @assert allunique(controls.date)
  Threads.@threads for i ∈ 1:length(gcrspdeciles)
    scrspdeciles = gcrspdeciles[i]
    dt = scrspdeciles[1,:date]
    dtM1y=dt - Year(1)

    #lagging one year window of returns
    crspwindow = view(crsp, (dtM1y .< crsp.date) .& (crsp.date .≤ dt),Fwindowcols)

    #compute the control matrix- only need to do this once per decile
    controlswindow = view(controls, (dtM1y .< controls.date) .& (controls.date .≤ dt),Fcontrols)

    controlmat = Matrix{Float64}(controlswindow) #control matrix for correlations
    for dec ∈ groupby(scrspdeciles, :d_plmomentum)
      #select only complete permno from the particular decile
      decwindow = innerjoin(crspwindow, dec[!, [:spermno]], on=:spermno)
      if !issorted(decwindow, [:spermno])
        sort!(decwindow, [:spermno])
      end
      @assert dec[1, :plretnmissing1y] .≥ length(unique(crspwindow.date)) .* minwindownonmissing

      #create the TxN matrix and compute the correlations
      crspbypermno = unstack(decwindow, :date, :spermno, Fplret)

      if (dt == testdate) && (dec[1, :d_plmomentum] == 10)
        crspwindow |> CSV.write(
          "$(PARAM[:testpath])\\$(crspprefix)_comom-crspwindowd-10.csv")
        controlswindow |> CSV.write(
          "$(PARAM[:testpath])\\$(crspprefix)_comom-ctrlwindowd-10.csv")
        crspbypermno  |> CSV.write(
          "$(PARAM[:testpath])\\$(crspprefix)_comom-crsppermnod-10.csv")
      end

      decilemat = Matrix{MFloat64}(crspbypermno[!, Not(:date)])

      #compute the total returns while taking care with missing values
      deciletotret = sum((x->coalesce(x,0.0)).(decilemat), dims=2) |> vec
      Ndecile = sum(decilemat .!== missing, dims=2) |> vec

      #NOTE- far better way to do the below, probably...
      colnames = setdiff(propertynames(crspbypermno), [:date])
      for (i, ri) ∈ enumerate(eachcol(decilemat))
        try
          dec[i, :xcor] = mpartialcor(ri, (deciletotret .- ri) ./ (Ndecile .- 1), controlmat)
        catch err
          dec[i, :xcor] = missing
          @warn "Failed to compute partial cor $err
            i=$i, colname: $(colnames[i]), dt: $dt"
        end
      end
    end
  end

  #print out some data for testing purposes
  testview = view(crspdeciles,
    ((s->parse(Int,string(s)[2:end])).(crspdeciles[!, :spermno])
      .% PARAM[:testpermnomult]) .== 0, :)
  testview |> CSV.write(
    "$(PARAM[:testpath])\\$(crspprefix)_comomdeciles.csv")

  #average the cross-sectional correlations, first within decile, then across top/bottom deciles
  crspdeciles = crspdeciles[crspdeciles.xcor .!== missing, :]
  crspdeciles = crspdeciles[isfinite.(crspdeciles.xcor), :]
  meanif(v, cond) = sum(v .* cond) ./ sum(cond)
  crspdeciles.d1 = crspdeciles.d_plmomentum .== 1
  crspdeciles.d10 = crspdeciles.d_plmomentum .== 10
  comom = combine(groupby(crspdeciles, :date),
    [:xcor, :d10] => meanif => :comom10,
    [:xcor, :d1] => meanif => :comom1)
  comom.comom = (comom.comom10 .+ comom.comom1) ./ 2
  sort!(comom, :date)

  #create log versions
  #comom.lGcomom = [missing; log.(comom.comom[2:end] ./ comom.comom[1:(end-1)])]
  #comom.lGcomom10 = [missing; log.(comom.comom10[2:end] ./ comom.comom10[1:(end-1)])]
  #comom.lGcomom1 = [missing; log.(comom.comom1[2:end] ./ comom.comom1[1:(end-1)])]

  return comom
end

function adjustindustry(crsp;
  industryfilename = PARAM[:plindustryfilename],
  industryfilepath = PARAM[:plindustryfilepath],
  sicfilename = PARAM[:plsicfilename],
  sicfilepath = PARAM[:plsicfilepath],
  missingunmatched=PARAM[:plindustrymissingunmatched],
  plretmin::Float64=PARAM[:plretmin],
  Fmc=:mc)

  lowercasesymbol(s) = s |> string |> lowercase |> Symbol

  #stack each sic code and with an ff code
  sic::DataFrame = CSV.File("$sicfilepath\\$(sicfilename).csv") |> DataFrame
  sort!(sic, [:ffcode, :beginsic])
  sic.ffdesc = sic.ffdesc .|> cleanpropertyname

  #create a lookup table of ffcd vs ffdesc
  fflookup = unique(sic[!, [:ffcode, :ffdesc]])
  @assert all(fflookup.ffcode .== collect(1:30))
  ffdesc = fflookup.ffdesc

  sic = reduce(append!!,
    (r->DataFrame(ffcode=r.ffcode, ffdesc=r.ffdesc,
      siccd=r.beginsic:r.endsic)).(eachrow(sic)))

  #validate
  @assert sum(nonunique(sic, :siccd)) == 0
  crsp.siccd = crsp.siccd .|> (i)->parsecrsp(Int, i)
  unknownsic = setdiff(unique(crsp.siccd), sic.siccd)

  #join the siccds
  crsp = leftjoin(crsp, sic, on=:siccd, matchmissing=:equal, validate=(false, true))
  if !issorted(crsp, [:permno, :date])
    sort!(crsp, [:permno, :date])
  end
  @assert sum(nonunique(crsp, [:permno, :date])) == 0

  #check what we are missing
  avgmissingsicmc = mean(combine(groupby(crsp, :date), [:ffcode, Fmc] =>
    ((ffcode, mc)->sum((ffcode .=== missing) .* mc)) => :missingsicmc).missingsicmc)
  avgtotalmc =  mean(combine(groupby(crsp, :date), [Fmc] => sum =>  :totmc).totmc)
  #@info "unmatched siccds: $unknownsic\n average mc missing: $avgmissingsicmc vs $avgtotalmc"

  #need to handle the mismatches- probably just link it to missing
  industry = loadff(industryfilename,
    ffpath=industryfilepath,
    fieldprefix = "",
    ffscale=true)

  Frets = ffdesc |> deepcopy
  if (!missingunmatched)
    crsp.ffcode[(ismissing).(crsp.ffcode)] .= 31 #new special code
    industry.unknown = zeros(nrow(industry))
    push!(Frets, :unknown)
  end

  if !issorted(crsp, [:permno, :date])
    sort!(crsp, [:permno, :date])
  end

  #now merge the returns
  crspindustry = mergereturns(crsp[!, [:permno, :date, :ret, :siccd]], industry,
    Fgroup=:permno,
    Frets=Frets,
    datefrequency=:day,
    testlabel="industry")

  industrymatrix = Matrix(crspindustry[!, Frets])

  #acquire the industry returns by using ffcode as the column#
  crspindustry.industryret = ((i,j)->
    j≡missing ? missing : industrymatrix[i,j]).(1:nrow(crsp), crsp.ffcode)

  #the below basically just drops the first row for which no industry returns are found
  nfound = combine(groupby(crspindustry, :date), [:industryret]=>
    ((iret)->(iret .!== missing)|>sum) => :numnonmissing)
  crspindustry = innerjoin(crspindustry, nfound, on=:date)
  sort!(crspindustry, :date)

  #optionally set fields with missing industry ret to 0.0
  if !missingunmatched
    crspindustry.industryret[ismissing.(crspindustry.industryret) .&
      (crspindustry.numnonmissing .> 0)] .= 0.0
  end

  crspindustry.plret = crsp.ret .- crspindustry.industryret

  crspindustry.plret .= crspindustry.plret .|> finiteormissing
  crspindustry.plret .=  ((x)->max(x,plretmin)).(crspindustry.plret)

  testpermno::Int = PARAM[:testpermnomult]
  testview = view(crspindustry, (crspindustry[!, :permno] .% testpermno) .== 0, :)
  testview |> CSV.write(
    "$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_industrypl.csv")

  return crspindustry.plret
end

#widget to create the cumulative return indices for plret
#=function computecumulativeindustryreturns(crsp)
  returns = crsp[!, [:permno, :date, :plret]] |> deepcopy

  #compute the cumulative adjusted returns
  returns.plcumret = missings(Float64, nrow(returns))
  returnsplretcomplete = view(returns, crsp.plret .!== missing, :)
  returnsplretcomplete.plcumret .= computecrspcumulativereturns(
    returnsplretcomplete|>DataFrame, Fret = :plret)

  return returns.plcumret
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# See https://rdrr.io/rforge/copula/man/Sibuya.html

struct Sibuya{T<:Real} <: Distributions.DiscreteUnivariateDistribution
    p::T
    function Sibuya(p::T) where {T <: Real}
        new{T}(p)
    end
end
function Distributions.rand(rng::Distributions.AbstractRNG, d::Sibuya{T}) where {T <: Real}
    u = rand(rng, T)
    if u <= d.p
        return T(1)
    end
    xMax = 1/eps(T)
    Ginv = ((1-u)*SpecialFunctions.gamma(1-d.p))^(-1/d.p)
    fGinv = floor(Ginv)
    if Ginv > xMax
        return fGinv
    end
    if 1-u < 1/(fGinv*SpecialFunctions.beta(fGinv,1-d.p))
        return ceil(Ginv)
    end
    return fGinv
end


text:
# ## `DatasetLoader` interface

"""
    abstract type DatasetLoader

A `DatasetLoader` defines how a dataset can made available and loaded.
See [`DataDepLoader`](#) as an example.

A `DatasetLoader` has to implement the following functions:

- [`loaddata`](#)
- [`makeavailable`](#)
- [`isavailable`](#)
"""
abstract type DatasetLoader end

function makeavailable end
function isavailable end
function loaddata end

# ## Loader for `DataDep`s

"""
    struct DataDepLoader(datadep) <: DatasetLoader

A dataset loader that uses DataDeps.jl to load datasets.
The DataDep has to be registered before creating the loader, and will
error otherwise.
"""
struct DataDepLoader <: DatasetLoader
    datadep::String
    function DataDepLoader(datadep)
        if !haskey(DataDeps.registry, datadep)
            throw(ArgumentError("DataDep \"$datadep\" does not exist."))
        end
        return new(datadep)
    end
end

isavailable(loader::DataDepLoader) = !isnothing(DataDeps.try_determine_load_path(loader.datadep, @__FILE__))

function makeavailable(loader::DataDepLoader)
    return DataDeps.resolve(loader.datadep, @__FILE__)
end

function loaddata(loader::DataDepLoader)
    isavailable(loader) || makeavailable(loader)
    return DataDeps.resolve(loader.datadep, @__FILE__)
end


text:
# Gondzio's multiple centrality correctors method

function update_rxs!(rxs_l, rxs_u, Hmin, Hmax, x_m_l_αΔp, u_m_x_αΔp, s_l_αΔp, s_u_αΔp, nlow, nupp)
  @inbounds @simd for i = 1:nlow
    rxs_l[i] = s_l_αΔp[i] * x_m_l_αΔp[i]
    if Hmin <= rxs_l[i] <= Hmax
      rxs_l[i] = 0
    elseif rxs_l[i] < Hmin
      rxs_l[i] -= Hmin
    else
      rxs_l[i] -= Hmax
    end
    if rxs_l[i] > Hmax
      rxs_l[i] = Hmax
    end
  end
  @inbounds @simd for i = 1:nupp
    rxs_u[i] = -s_u_αΔp[i] * u_m_x_αΔp[i]
    if Hmin <= -rxs_u[i] <= Hmax
      rxs_u[i] = 0
    elseif -rxs_u[i] < Hmin
      rxs_u[i] += Hmin
    else
      rxs_u[i] += Hmax
    end
    if rxs_u[i] < -Hmax
      rxs_u[i] = -Hmax
    end
  end
end

function multi_centrality_corr!(
  dda::DescentDirectionAllocsPC{T},
  pad::PreallocatedData{T},
  pt::Point{T},
  α_pri::T,
  α_dual::T,
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  cnts::Counters,
  res::AbstractResiduals{T},
  T0::DataType,
) where {T <: Real}
  iter_c = 0 # current number of correction iterations
  corr_flag = true #stop correction if false
  # for storage issues Δ_aff = Δp  and Δ_cc = Δm
  dda.Δxy_aff .= itd.Δxy
  dda.Δs_l_aff .= itd.Δs_l
  dda.Δs_u_aff .= itd.Δs_u
  @inbounds while iter_c < cnts.kc && corr_flag
    # Δp = Δ_aff + Δ_cc
    δα, γ, βmin, βmax = T(0.1), T(0.1), T(0.1), T(10)
    α_p2, α_d2 = min(α_pri + δα, one(T)), min(α_dual + δα, one(T))
    update_pt_aff!(
      dda.x_m_l_αΔ_aff,
      dda.u_m_x_αΔ_aff,
      dda.s_l_αΔ_aff,
      dda.s_u_αΔ_aff,
      dda.Δxy_aff,
      dda.Δs_l_aff,
      dda.Δs_u_aff,
      itd.x_m_lvar,
      itd.uvar_m_x,
      pt.s_l,
      pt.s_u,
      α_p2,
      α_d2,
      id.ilow,
      id.iupp,
    )
    μ_p = compute_μ(
      dda.x_m_l_αΔ_aff,
      dda.u_m_x_αΔ_aff,
      dda.s_l_αΔ_aff,
      dda.s_u_αΔ_aff,
      id.nlow,
      id.nupp,
    )

    σ = (μ_p / itd.μ)^3
    Hmin, Hmax = βmin * σ * itd.μ, βmax * σ * itd.μ

    # corrector-centering step
    update_rxs!(
      dda.rxs_l,
      dda.rxs_u,
      Hmin,
      Hmax,
      dda.x_m_l_αΔ_aff,
      dda.u_m_x_αΔ_aff,
      dda.s_l_αΔ_aff,
      dda.s_u_αΔ_aff,
      id.nlow,
      id.nupp,
    )
    itd.Δxy .= 0
    itd.Δxy[id.ilow] .+= dda.rxs_l ./ itd.x_m_lvar
    itd.Δxy[id.iupp] .+= dda.rxs_u ./ itd.uvar_m_x
    out = solver!(itd.Δxy, pad, dda, pt, itd, fd, id, res, cnts, T0, :cc)
    itd.Δs_l .= @views .-(dda.rxs_l .+ pt.s_l .* itd.Δxy[id.ilow]) ./ itd.x_m_lvar
    itd.Δs_u .= @views (dda.rxs_u .+ pt.s_u .* itd.Δxy[id.iupp]) ./ itd.uvar_m_x

    itd.Δxy .+= dda.Δxy_aff
    itd.Δs_l .+= dda.Δs_l_aff
    itd.Δs_u .+= dda.Δs_u_aff
    α_p2, α_d2 =
      compute_αs(pt.x, pt.s_l, pt.s_u, fd.lvar, fd.uvar, itd.Δxy, itd.Δs_l, itd.Δs_u, id.nvar)

    if α_p2 >= α_pri + γ * δα && α_d2 >= α_dual + γ * δα
      iter_c += 1
      dda.Δxy_aff .= itd.Δxy
      dda.Δs_l_aff .= itd.Δs_l
      dda.Δs_u_aff .= itd.Δs_u
      α_pri, α_dual = α_p2, α_d2
    else
      itd.Δxy .= dda.Δxy_aff
      itd.Δs_l .= dda.Δs_l_aff
      itd.Δs_u .= dda.Δs_u_aff
      corr_flag = false
    end
  end
  return α_pri, α_dual
end

# function to determine the number of centrality corrections
function nb_corrector_steps!(cnts::Counters, time_fact::T, time_solve::T) where {T <: Real}
  rfs = time_fact / time_solve
  if rfs <= T(10)
    cnts.kc = 0
  elseif T(10) < rfs <= T(30)
    cnts.kc = 1
  elseif T(30) < rfs <= T(50)
    cnts.kc = 2
  elseif rfs > T(50)
    cnts.kc = 3
  else
    p = Int(round(rfs / 50))
    cntd.kc = p + 2
    if cnts.kc > 10
      cnts.kc = 10
    end
  end
end


text:
module Hiccup

using Lazy, MacroTools, Compat

export Node, tag, @tags, @exporttags

# Void elements; not allowed to contain content
# See: http://www.w3.org/TR/html5/syntax.html#void-elements
const VOID_ELEMENTS = Set([:area, :base, :br, :col, :embed, :hr, :img, :input,
  :keygen, :link, :meta, :param, :source, :track, :wbr])

# Basic Types

type Node{tag}
  attrs::Dict{Any, Any}
  children::Vector{Any}
end

tag{T}(node::Node{T}) = T
attrs(node::Node) = node.attrs
children(node::Node) = node.children

isvoid(tag::Symbol) = tag in VOID_ELEMENTS
isvoid(node::Node) = isvoid(tag(node))

function Node(tag::Symbol, attrs::Dict = Dict(), content::AbstractVector = c(); kws...)
  if isvoid(tag) && !isempty(content)
    throw(ArgumentError("Void tag <$tag> cannot have content."))
  end
  Node{tag}(isempty(kws) ? attrs : merge(attrs, Dict([kws...])), content)
end

Node(tag::Symbol, attrs::Associative, content; kws...) =
  Node(tag, attrs, c(content); kws...)

Node(tag::Symbol, content; kws...) =
  Node(tag, Dict(), content; kws...)

Node(tag::Symbol, content::Node...) = Node(tag, collect(content))

# CSS selector parsing

function cssparse(s)
  trimfirst(s) = s[2:end]
  attrs = Dict()
  id = match(r"#-?[_a-zA-Z][_a-zA-Z0-9-]*", s)
  id == nothing || (attrs[:id] = trimfirst(id.match))
  classes = matchall(r"\.-?[_a-zA-Z][_a-zA-Z0-9-]*", s)
  isempty(classes) || (attrs[:class] = map(trimfirst, classes))
  return attrs
end

Node(tag::Symbol, selector::AbstractString, props::Dict, args...; kws...) =
  Node(tag, merge!(cssparse(selector), props), args...; kws...)

Node(tag::Symbol, selector::AbstractString, content, args...; kws...) =
  Node(tag, cssparse(selector), content, args...; kws...)

# Rendering

export htmlescape

htmlescape(s::AbstractString) =
  @> s replace("&", "&") replace("<", "<") replace(">", ">")

attrstring(xs::Vector) = join(xs, " ")
attrstring(x) =
  @> x string htmlescape replace("\"", """) replace("'", "'")
attrstring(d::Dict) = @as _ d ["$(t[1])=\"$(attrstring(t[2]))\"" for t in _] join(_, " ")

render(io::IO, s::AbstractString) = print(io, htmlescape(s))

function render(io::IO, node::Node)
  print(io, "<", tag(node))
  isempty(attrs(node)) || print(io, " ", attrstring(attrs(node)))
  if isvoid(node)
    print(io, " />")
  else
    print(io, ">")
    render(io, children(node))

    print(io, "</", tag(node), ">")
  end
end

function render(io::IO, xs::Vector)
  for x in xs
    render(io, x)
  end
end

render(io::IO, x) = @compat show(io, MIME"text/html"(), x)

@compat Base.show(io::IO, ::MIME"text/html", node::Node) = render(io, node)

Base.show(io::IO, node::Node) = render(io, node)

Node(tag::Symbol, io::IO, args...) = render(io, Node(tag, args...))

# Specific elements

tags(t) = :(($t)(args...; kws...) = Node($(Expr(:quote, t)), args...; kws...))

macro tags(ts)
  isexpr(ts, Symbol) && (ts = Expr(:tuple, ts))
  @assert isexpr(ts, :tuple)
  quote
    $([tags(t) for t in ts.args]...)
  end |> esc
end

macro exporttags(ts)
  quote
    @tags $(esc(ts))
    $(Expr(:export, (isexpr(ts, Symbol) ? [ts] : ts.args)...))
  end
end

@exporttags span, a,
            h1, h2, h3,
            html, head, body,
            pre, code,
            img, style,
            ol, ul, li, table, tr, td,
            strong

if VERSION < v"0.4"
  @exporttags div
else
  @tags div
end

end # module


text:
e
π, pi
sqrt(x)
log(x)
exp(x)
abs(x)
floor(x)
ceil(x)
x^y

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode