text:
### L2 - obsolete
function L2_prod(basis::Union{Basis, TensorProductBasis})
    return mass_matrix(basis)
end

function L2_prod(basis1::Basis, basis2::Basis)
    return mass_matrix(basis1, basis2)
end

function ⊙(f::Function, b::Union{Basis, TensorProductBasis})
    return L2_prod(f, b)
end

function ⊙(b::Union{Basis, TensorProductBasis}, f::Function)
    return L2_prod(f, b)
end

function ⊙(a::Union{Basis, TensorProductBasis}, b::Union{Basis, TensorProductBasis})
    if a == b
        return L2_prod(a)
    else
        L2_prod(a, b)
    end
end

### derivative bases
# Maybe need to generalize this to d1, d2, d3,...
abstract type grad end
abstract type curl end
abstract type div end

struct d_basis{D, T <: Union{Basis, TensorProductBasis}}
    basis::T
end

function ∇(basis::Union{Basis, TensorProductBasis})
    return d_basis{grad, typeof(basis)}(basis)
end

function L2_prod(gbasis::d_basis{grad, T}) where T <: Basis
    return stiffness_matrix(derivative(gbasis.basis))
end

# Needed for poisson problem
function L2_prod(gbasis::d_basis{grad, BSplineTensorProductBasis{N, periodic}}) where N
    components = []
    basis = gbasis.basis
    for j in 1:N
        matrices = []
        for i in 1:N
            if i == j
                b = derivative(basis[i])
                push!(matrices, stiffness_matrix(b))
            else
                push!(matrices, mass_matrix(basis[i]))
            end
        end

        comp = ⊗(matrices...)
        push!(components, comp)
    end

    return sum(components)
end

function L2_prod(gbasis::d_basis{grad, BSplineTensorProductBasis{N, dirichlet}}) where N
    components = []
    basis = gbasis.basis
    ϵ = 1e-14
    for j in 1:N
        matrices = []
        for i in 1:N
            if i == j
                b = derivative(basis[i])
                sm = stiffness_matrix(b)
                sm[1,1] *= (1 +1/ϵ)
                sm[end, end] *= (1 +1/ϵ)
                push!(matrices, sm)
            else
                mm = mass_matrix(basis[i])
                mm[1,1] *= (1 +1/ϵ)
                mm[end, end] *= (1 +1/ϵ)
                push!(matrices, mm)
            end
        end

        comp = ⊗(matrices...)
        push!(components, comp)
    end

    return sum(components)
end

function ⊙(a::Union{d_basis{grad, T}, d_basis{grad, D}}, b::Union{d_basis{grad, T}, d_basis{grad, D}}) where {T <: Basis, D <: TensorProductBasis}
    if a == b
        return L2_prod(a)
    else
        error("todo")
    end
end

### discrete differentials
function der_mat(N::Int)
    c = zeros(N)
    c[1] = 1
    c[2] = -1

    return Circulant(c)
end

function der_mat(basis::BSplineBasis{periodic})

    return der_mat(basis.N)
end

function der_mat(basis::BSplineBasis{dirichlet})
    D = der_mat(length(basis.B))[2:end, :]

    return D
end

#This way its hard to generalize the notation I think..
#Would be better to give them names like d{t, N} corresponding to d^t in N dimensions
struct discrete_differential{T, N}
    D::SArray{Tuple{N}} #1d diff matrices
end

function disc_grad(basis::BSplineBasis{T}) where T
    return der_mat(basis)
end

function disc_grad(N::Int)
    return der_mat(N)
end

function disc_grad(N::Union{NTuple, BSplineTensorProductBasis{D, T}}) where {D, T}
    return disc_diff(grad, N)
end

function disc_curl(N::Union{NTuple, BSplineTensorProductBasis{D, T}}) where {D, T}
    return disc_diff(curl, N)
end

function disc_div(N::Union{NTuple, BSplineTensorProductBasis{D, T}}) where {D, T}
    return disc_diff(div, N)
end

#1D case
#function disc_diff(N)
#    D = der_mat(N)#
#
#    return D
#end

#2D case
#here we only have gradient and curl in the primal diagram, curl and div in the dual diagram, which we write similarly
#todo add dual operators properly
function disc_diff(T, N::NTuple{2})
    D1 = der_mat(N[1]) ⊗ circ_id(N[2])
    D2 = circ_id(N[1]) ⊗ der_mat(N[2])

    return discrete_differential{T, 2}(SVector(D1, D2))
end

function disc_diff(T, basis::BSplineTensorProductBasis{2, ty}) where ty
    if T == grad
        N1 = length(basis[1].B)
        N2 = length(basis[2].B)

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    elseif T == curl
        N1 = length(basis[1]'.B)
        N2 = length(basis[2]'.B)

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    end

    return discrete_differential{T, 2}(SVector(D1, D2))
end

function *(C::discrete_differential{grad, 2}, fh::SArray)
    return SVector(C.D[1] * fh, C.D[2] * fh)
end

function *(C::discrete_differential{curl, 2}, fh::SArray)
    return SVector{size(C.D[1])[1]}(C.D[1] * fh[2] - C.D[2] * fh[1])
end

function matrix(C::discrete_differential{grad, 2})
    return [C.D[1];
            C.D[2]]
end

function matrix(C::discrete_differential{curl, 2})
    return [-1 .* C.D[2] C.D[1]]
end

function adjoint(C::discrete_differential{grad, 2})
    return discrete_differential{curl, 2}(SVector(C.D[2]',  -1 .* C.D[1]'))
end

function adjoint(C::discrete_differential{curl, 2})
    return discrete_differential{grad, 2}(SVector( -1 .* C.D[2]', C.D[1]'))
end

#3D case
# brackets so that the kronecker product ordering coincides with the mass matrix
function disc_diff(T, N::NTuple{3})
    D1 = der_mat(N[1]) ⊗ (circ_id(N[2]) ⊗ circ_id(N[3]))
    D2 = circ_id(N[1]) ⊗ (der_mat(N[2]) ⊗ circ_id(N[3]))
    D3 = circ_id(N[1]) ⊗ (circ_id(N[2]) ⊗ der_mat(N[3]))

    return discrete_differential{T, 3}(SVector(D1, D2, D3))
end

##TODO FINISH FOR DIRICHLET 3D
function disc_diff(T, basis::BSplineTensorProductBasis{3, ty}) where ty
    if T == grad
        N1 = length(basis[1].B)
        N2 = length(basis[2].B)
        N3 = length(basis[3].B)

        D1 = der_mat(basis[1]) ⊗ (circ_id(N2) ⊗ circ_id(N3))
        D2 = circ_id(N1) ⊗ (der_mat(basis[2]) ⊗ circ_id(N3))
        D3 = circ_id(N1) ⊗ (circ_id(N2) ⊗ der_mat(basis[3]))
    elseif T == curl
        Nd1 = length(basis[1]'.B)
        Nd2 = length(basis[2]'.B)
        Nd3 = length(basis[3]'.B)

        D1 = circ_id(N2) ⊗ der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    elseif T == div
        N1 = length(basis[1]'.B)
        N2 = length(basis[2]'.B)

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    end

    return discrete_differential{T, 2}(SVector(D1, D2))
end

function *(C::discrete_differential{grad, 3}, fh::SArray)
    return SVector(C.D[1] * fh, C.D[2] * fh , C.D[3] * fh)
end

function *(C::discrete_differential{curl, 3}, fh::SArray{Tuple{3}})
    return SVector(C.D[3] * fh[2] - C.D[2] * fh[3], -C.D[3]* fh[1] + C.D[1] * fh[3], C.D[2] * fh[1] - C.D[1] * fh[2])
end

function *(C::discrete_differential{div, 3}, fh::SArray{Tuple{3}})
    return SVector{length(fh[1])}(C.D[1] * fh[1] + C.D[2] * fh[2] + C.D[3] * fh[3])
end

function matrix(C::discrete_differential{grad, 3})
    return [C.D[1];
            C.D[2];
            C.D[3]]
end

function matrix(C::discrete_differential{curl, 3})
    return [zero(C.D[1]) -C.D[3] C.D[2];
            C.D[3]  zero(C.D[2]) -C.D[1];
            -C.D[2] C.D[1] zero(C.D[3])]
end

function matrix(C::discrete_differential{div, 3})
    return [C.D[1] C.D[2] C.D[3]]
end

#since we don't have a matrix, we have to manually implement the transpose
function adjoint(C::discrete_differential{grad, 3})
    return discrete_differential{div, 3}(SVector(C.D[1]', C.D[2]', C.D[3]'))
end

function adjoint(C::discrete_differential{curl, 3})
    return discrete_differential{curl, 3}(-1 .* SVector(C.D[1]', C.D[2]', C.D[3]'))
end

function adjoint(C::discrete_differential{div, 3})
    return discrete_differential{grad, 3}(SVector(C.D[1]', C.D[2]', C.D[3]'))
end


text:
function OEMParser()
    return OEMParser(())
end

function OEMParser(arg0::DataContext)
    return OEMParser((DataContext,), arg0)
end

function parse(obj::OEMParser, arg0::BufferedReader, arg1::JString)
    return jcall(obj, "parse", OEMFile, (BufferedReader, JString), arg0, arg1)
end

function parse(obj::OEMParser, arg0::InputStream)
    return jcall(obj, "parse", OEMFile, (InputStream,), arg0)
end

function parse(obj::OEMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", OEMFile, (InputStream, JString), arg0, arg1)
end

function parse(obj::OEMParser, arg0::JString)
    return jcall(obj, "parse", OEMFile, (JString,), arg0)
end

function with_conventions(obj::OEMParser, arg0::IERSConventions)
    return jcall(obj, "withConventions", OEMParser, (IERSConventions,), arg0)
end

function with_data_context(obj::OEMParser, arg0::DataContext)
    return jcall(obj, "withDataContext", OEMParser, (DataContext,), arg0)
end

function with_international_designator(obj::OEMParser, arg0::jint, arg1::jint, arg2::JString)
    return jcall(obj, "withInternationalDesignator", OEMParser, (jint, jint, JString), arg0, arg1, arg2)
end

function with_mission_reference_date(obj::OEMParser, arg0::AbsoluteDate)
    return jcall(obj, "withMissionReferenceDate", OEMParser, (AbsoluteDate,), arg0)
end

function with_mu(obj::OEMParser, arg0::jdouble)
    return jcall(obj, "withMu", OEMParser, (jdouble,), arg0)
end

function with_simple_eop(obj::OEMParser, arg0::jboolean)
    return jcall(obj, "withSimpleEOP", OEMParser, (jboolean,), arg0)
end

ext: jl


text:
using .Pairs

"""
    HTTP.Client(;args...)

A type to facilitate connections to remote hosts, send HTTP requests, and manage state between requests.
Additional keyword arguments can be passed that will get transmitted with each HTTP request:

  * `chunksize::Int`: if a request body is larger than `chunksize`, the "chunked-transfer" http mechanism will be used and chunks will be sent no larger than `chunksize`; default = `nothing`
  * `connecttimeout::Float64`: sets a timeout on how long to wait when trying to connect to a remote host; default = Inf. Note that while setting a timeout will affect the actual program control flow, there are current lower-level limitations that mean underlying resources may not actually be freed until their own timeouts occur (i.e. libuv sockets only timeout after 75 seconds, with no option to configure)
  * `readtimeout::Float64`: sets a timeout on how long to wait when receiving a response from a remote host; default = Int
  * `tlsconfig::TLS.SSLConfig`: a valid `TLS.SSLConfig` which will be used to initialize every https connection; default = `nothing`
  * `maxredirects::Int`: the maximum number of redirects that will automatically be followed for an http request; default = 5
  * `allowredirects::Bool`: whether redirects should be allowed to be followed at all; default = `true`
  * `forwardheaders::Bool`: whether user-provided headers should be forwarded on redirects; default = `false`
  * `retries::Int`: # of times a request will be tried before throwing an error; default = 3
  * `managecookies::Bool`: whether the request client should automatically store and add cookies from/to requests (following appropriate host-specific & expiration rules); default = `true`
  * `statusraise::Bool`: whether an `HTTP.StatusError` should be raised on a non-2XX response status code; default = `true`
  * `insecure::Bool`: whether an "https" connection should allow insecure connections (no TLS verification); default = `false`
  * `canonicalizeheaders::Bool`: whether header field names should be canonicalized in responses, e.g. `content-type` is canonicalized to `Content-Type`; default = `true`
  * `logbody::Bool`: whether the request body should be logged when `verbose=true` is passed; default = `true`
"""
mutable struct Client
    # cookies are stored in-memory per host and automatically sent when appropriate
    cookies::Dict{String, Set{Cookie}}
    # global request settings
    options::(VERSION > v"0.7.0-DEV.2338" ? NamedTuple : Vector{Pair{Symbol,Any}})
end

if VERSION > v"0.7.0-DEV.2338"
Client(;options...) = Client(Dict{String, Set{Cookie}}(), merge(NamedTuple(), options))
else
Client(;options...) = Client(Dict{String, Set{Cookie}}(), options)
end
global const DEFAULT_CLIENT = Client()

# build Request
function request(client::Client, method, url::URI;
                 headers=Header[],
                 body="",
                 enablechunked::Bool=true,
                 stream::Bool=false,
                 verbose::Bool=false,
                 args...)

    # Add default values from client options to args...
    if VERSION > v"0.7.0-DEV.2338"
    args = merge(client.options, args)
    getarg = Base.get
    else
    for option in client.options
        defaultbyfirst(args, option)
    end
    getarg = getkv
    end
    newargs = Pair{Symbol,Any}[]

    if getarg(args, :chunksize, nothing) != nothing
        Base.depwarn(
        "The chunksize= option is deprecated and has no effect.\n" *
        "Use a HTTP.open and pass chunks of the desired size to `write`.",
        :chunksize)
    end

    if getarg(args, :connecttimeout, Inf) != Inf
        Base.depwarn(
        "The connecttimeout= is deprecated and has no effect.\n" *
        "See https://github.com/JuliaWeb/HTTP.jl/issues/114\n",
        :connecttimeout)
    end

    if getarg(args, :tlsconfig, nothing) != nothing
        Base.depwarn(
        "The tlsconfig= option is deprecated. Use sslconfig=::MbedTLS.SSLConfig",
        :tlsconfig)
        setkv(newargs, :sslconfig, getarg(args, :tlsconfig))
    end

    if getarg(args, :allowredirects, nothing) != nothing
        Base.depwarn(
        "The allowredirects= option is deprecated. Use redirect=::Bool",
        :allowredirects)
        setkv(newargs, :redirect, getarg(args, :allowredirects))
    end

    if getarg(args, :managecookies, nothing) != nothing
        Base.depwarn(
        "The managecookies= option is deprecated. Use cookies=::Bool",
        :managecookies)
        setkv(newargs, :cookies, getarg(args, :managecookies))
    end
    setkv(newargs, :cookiejar, client.cookies)

    if getarg(args, :statusraise, nothing) != nothing
        Base.depwarn(
        "The statusraise= options is deprecated. Use status_exception=::Bool",
        :statusraise)
        setkv(newargs, :status_exception, getarg(args, :statusraise))
    end

    if getarg(args, :insecure, nothing) != nothing
        Base.depwarn(
        "The insecure= option is deprecated. Use require_ssl_verification=::Bool",
        :insecure)
        setkv(newargs, :require_ssl_verification, !getarg(args, :insecure))
    end

    m = string(method)
    h = mkheaders(headers)
    if stream
        setkv(newargs, :response_stream, BufferStream())
    end

    if isa(body, Dict)
        body = HTTP.Form(body)
        setbyfirst(h, "Content-Type" =>
                            "multipart/form-data; boundary=$(body.boundary)")
        setkv(newargs, :bodylength, length(body))
    end

    if !enablechunked && isa(body, IO)
        body = read(body)
    end

    if VERSION > v"0.7.0-DEV.2338"
    args = merge(args, newargs)
    else
    for newarg in newargs
        defaultbyfirst(args, newarg)
    end
    end

    return request(m, url, h, body; args...)
end

for f in [:get, :post, :put, :delete, :head,
          :trace, :options, :patch, :connect]
    meth = f_str = uppercase(string(f))
    @eval begin
        ($f)(client::Client, url::URI; kw...) = request(client, $meth, url; kw...)
        ($f)(client::Client, url::AbstractString; kw...) = request(client, $meth, URI(url); kw...)
        ($f)(url::URI; kw...) = request(DEFAULT_CLIENT, $meth, url; kw...)
        ($f)(url::AbstractString; kw...) = request(DEFAULT_CLIENT, $meth, URI(url); kw...)
    end
end

function download(url::AbstractString, file; threshold::Int=50000000, verbose::Bool=false, query="", args...)
    res = request(GET, url; verbose=verbose, query=query, stream=true, args...)
    body = HTTP.body(res)
    file = Base.get(HTTP.headers(res), "Content-Encoding", "") == "gzip" ? string(file, ".gz") : file
    threshold_step = threshold
    nbytes = 0
    open(file, "w") do f
        while !eof(body)
            nbytes += write(f, readavailable(body))
            if verbose && nbytes > threshold
                println("[$(Dates.now())]: downloaded $nbytes bytes...")
                flush(STDOUT)
                threshold += threshold_step
            end
        end
        length(body) > 0 && write(f, readavailable(body))
    end
    return file
end


text:
# Factorization Type definitions

# Reference: I.S. Dhillon, B.N. Parlett. Orthogonal Eigenvectors and Relative Gaps.
# Siam J. Matrix Anal. Appl. Vol. 25, No. 3, pp. 858-899 (2003)

module mr3

export LDL, NDN, UDU, getvec

# LDL factorization type for tridiagonal matrices
type LDL{T}
    d::Vector{T} # diagonal entries of D
    l::Vector{T} # lower diagonal entries of L
    function LDL{T}(d::Vector{T}, l::Vector{T})
        length(d)==length(l)+1 ? new(d,l) : throw(DimensionMismatch())
    end
end

# UDU factorization type for tridiagonal matrices
type UDU{T}
    d::Vector{T} # diagonal entries of D
    u::Vector{T} # upper diagonal entries of U
    function UDU{T}(d::Vector{T}, u::Vector{T})
        length(d)==length(u)+1 ? new(d,u) : throw(DimensionMismatch())
    end
end

# Storage of all twisted factorizations
type NDN{T}
    dl::Vector{T} # diagonal entries of D+
    l::Vector{T} # sub-diagonal entries of L+
    du::Vector{T} # diagonal entries of D-
    u::Vector{T} # super-diagonal entries of U-
    γ::Vector{T} # one choice for each twisted factorization
    function NDN{T}(dl::Vector{T}, l::Vector{T}, du::Vector{T}, u::Vector{T}, γ::Vector{T})
        b = length(dl)==length(du) && length(l)==length(u) && length(dl)==(length(l)+1) && length(dl)==length(γ)
        b ? new(dl, l, du, u, γ) : throw(DimensionMismatch())
    end
end

# LDL constructors
LDL{T}(d::AbstractVector{T}, l::AbstractVector{T})=LDL{T}(copy(d), copy(l))
LDL(n::Integer)=LDL(zeros(n), zeros(n-1)) # blank constructor

# UDU constructors
UDU{T}(d::AbstractVector{T}, u::AbstractVector{T})=UDU{T}(copy(d), copy(u))
UDU(n::Integer)=UDU(zeros(n), zeros(n-1)) # blank constructor

#NDN constructors
function NDN{T}(dl::AbstractVector{T}, l::AbstractVector{T}, du::AbstractVector{T}, u::AbstractVector{T}, γ::AbstractVector{T})
    return NDN{T}(copy(dl), copy(l), copy(du), copy(u), copy(γ))
end
NDN(n::Integer)=NDN(zeros(n), zeros(n-1), zeros(n), zeros(n-1), zeros(n)) # blank constructor

# get LDL factorization of a SymTridiagonal factorization
function LDL{T}(M::SymTridiagonal{T})
    n = length(M.dv)
    Md = copy(M.dv)
    Me = copy(M.ev)
    # vectors for LDL factorization
    d = Vector{T}(n)
    l = Vector{T}(n-1)
    for i=1:n-1
        d[i] = Md[i];
        l[i] = Me[i]/d[i];
        Md[i+1] = Md[i+1] - l[i]*d[i]*l[i]
    end
    d[n] = Md[n]
    return LDL(d,l)
end

function fullLDL(F::LDL)
    n = length(F.d)
    L = diagm(F.l,-1)
    L = L + I
    D = diagm(F.d)
    return L, D
end

function dstqds(F::LDL, μ::Real)
    # Transform F = LDL to G = LDL-μI
    # where G is also a LDL factorization
    # ref Dhillon Parlett 2003 pp.869
    n = length(F.d)
    G = LDL(n)
    s = -μ
    for i = 1:n-1
        G.d[i] = s + F.d[i]
        G.l[i] = F.l[i]*F.d[i]/G.d[i]
        s = G.l[i]*F.l[i]*s - μ
    end
    G.d[n] = s + F.d[n]
    return G
end

function dqds(F::LDL, μ::Real)
    # transform F = LDL to G = LDL - μI
    # where G is a UDU factorization
    # ref Dhillon Parlett 2003 pp.870
    n = length(F.d)
    G = UDU(n)
    p = zeros(n)
    p[n] = F.d[n] - μ
    for i = n-1:-1:1
        G.d[i+1] = F.d[i]*F.l[i]^2 + p[i+1]
        t = F.d[i]/G.d[i+1]
        G.u[i] = F.l[i]*t
        p[i] = p[i+1]*t - μ
    end
    G.d[1] = p[1]
    return G
end

function fullUDU(F::UDU)
    n = length(F.d)
    U = diagm(F.u,1)
    U = U + I
    D = diagm(F.d)
    return D, U
end

function dtwqds(F::LDL, μ::Real)
    # transform F = LDL to G = LDL - μI
    # where G is a NDN factorization
    # ref Dhillon Parlett 2003 pp.871
    n = length(F.d)
    G = NDN(zeros(n), zeros(n-1), zeros(n), zeros(n-1), zeros(n))
    s = Vector{Real}(n)
    p = Vector{Real}(n)

    # L+D+L+ factorization
    s[1] = -μ
    for i = 1:n-1
        G.dl[i] = s[i] + F.d[i]
        G.l[i] = F.l[i]*F.d[i]/G.dl[i]
        s[i+1] = G.l[i]*F.l[i]*s[i] - μ
    end
    G.dl[n] = s[n] + F.d[n]

    # U-D-U- factorization
    p[n] = F.d[n] - μ
    for i = n-1:-1:1
        G.du[i+1] = F.d[i]*F.l[i]^2 + p[i+1]
        t = F.d[i]/G.du[i+1]
        G.u[i] = F.l[i]*t
        p[i] = p[i+1]*t - μ
    end
    G.du[1] = p[1]

    # compute γ
    for r  = 1:n-1
        G.γ[r] = s[r] + F.d[r]*p[r+1]/G.du[r+1]
    end
    G.γ[n] = s[n] + F.d[n]

    return G
end

function twisted_eigenvector(F::NDN, r::Integer, G::LDL)
    # transform find smallest eigenvector of F = NDN
    # r tells which NDN factorization to use
    # INPUTS:
    #   F - NDN factorizations of LDL-μI.  Generate using dtwqds
    #   r - which twisted factorization should be used
    #   G - unshifted LDL factorization.  Used if shifted factorization has 0 on diagonal
    # ref Dhillon Parlett 2003 pp.866

    n = length(F.dl)
    z = zeros(n)
    z[r] = 1
    for i = r-1:-1:1
        if F.dl[i]!=0
            z[i] = -F.l[i]*z[i+1]
        else
            z[i] = -G.d[i+1] * G.l[i+1] / (G.d[i] * G.l[i]) * z[i+2]
        end
    end
    for j = r:n-1
        if F.du[j+1]!=0
            z[j+1] = -F.u[j]*z[j]
        else
            z[j] = -G.d[j-1] * G.l[j-1] / (G.d[j] * G.l[j]) * z[j-1]
        end
    end

    # normalize
    z = z / norm(z)
    return z
end

function argmin(v::AbstractVector)
    # return index with smallest value

    curmin = Inf # current minimum
    amin = 0 # current argmin
    for i = 1:length(v)
        if v[i] <= curmin
            amin = i
            curmin = v[i]
        end
    end
    return amin
end

function getvec(F::LDL, λ::Real)
    # get eigenvector with approximate eigenvalue λ
    # use twisted factorizations
    # ref Dhillon Parlett 2003 pp.865-866

    G = dtwqds(F, λ)
    r = argmin(abs(G.γ))
    return twisted_eigenvector(G, r, F)
end

end # module


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    delete_duplicates!(db)
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
# Please run from the `basic_gppp` directory.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Random, Plots, Stheno
using Stheno: @model, EQ

###########################  Define our model  ###########################

# Define a distribution over f₁, f₂, and f₃, where f₃(x) = f₁(x) + f₂(x).
@model function model()
    f₁ = GP(randn(), EQ())
    f₂ = GP(EQ())
    f₃ = f₁ + f₂
    return f₁, f₂, f₃
end

# Construct the processes in the model.
f₁, f₂, f₃ = model()

# Randomly sample `N₁` and `N₃` locations at which to observe `f₁` and `f₃` respectively.
rng, N₁, N₃ = MersenneTwister(123546), 10, 11
x₁, x₃ = sort(rand(rng, N₁) * 10), sort(rand(rng, N₃) * 10)

# Generate some toy of `f₁` and `f₃`, `y₁` and `y₃` respectively.
y₁, y₃ = rand(rng, [f₁(x₁), f₃(x₃)])

# Compute the posterior processes.
(f₁_post, f₂_post, f₃_post) = (f₁, f₂, f₃) | (f₁(x₁)←y₁, f₃(x₃)←y₃)

# Define some plotting stuff.
Np, S = 500, 25;
xp = range(-2.5, stop=12.5, length=Np);

# Sample jointly from the posterior over each process.
f₁_post_xp, f₂_post_xp, f₃_post_xp =
    rand(rng, [f₁_post(xp, 1e-9), f₂_post(xp, 1e-9), f₃_post(xp, 1e-9)], S);

# Compute posterior marginals.
ms₁ = marginals(f₁_post(xp));
ms₂ = marginals(f₂_post(xp));
ms₃ = marginals(f₃_post(xp));

μf₁′, σf₁′ = mean.(ms₁), std.(ms₁);
μf₂′, σf₂′ = mean.(ms₂), std.(ms₂);
μf₃′, σf₃′ = mean.(ms₃), std.(ms₃);

###########################  Plot results  ###########################

gr();
posterior_plot = plot();

# Plot posterior marginal variances
plot!(posterior_plot, xp, [μf₁′ μf₁′];
    linewidth=0.0,
    fillrange=[μf₁′ .- 3 .* σf₁′, μf₁′ .+ 3 * σf₁′],
    fillalpha=0.3,
    fillcolor=:red,
    label="",
);
plot!(posterior_plot, xp, [μf₂′ μf₂′];
    linewidth=0.0,
    fillrange=[μf₂′ .- 3 .* σf₂′, μf₂′ .+ 3 * σf₂′],
    fillalpha=0.3,
    fillcolor=:green,
    label="",
);
plot!(posterior_plot, xp, [μf₃′ μf₃′];
    linewidth=0.0,
    fillrange=[μf₃′ .- 3 .* σf₃′, μf₃′ .+ 3 * σf₃′],
    fillalpha=0.3,
    fillcolor=:blue,
    label="",
);

# Plot joint posterior samples
plot!(posterior_plot, xp, f₁_post_xp,
    linecolor=:red,
    linealpha=0.2,
    label="",
);
plot!(posterior_plot, xp, f₂_post_xp,
    linecolor=:green,
    linealpha=0.2,
    label="",
);
plot!(posterior_plot, xp, f₃_post_xp,
    linecolor=:blue,
    linealpha=0.2,
    label="",
);

# Plot posterior means
plot!(posterior_plot, xp, μf₁′;
    linecolor=:red,
    linewidth=2.0,
    label="f1",
);
plot!(posterior_plot, xp, μf₂′;
    linecolor=:green,
    linewidth=2.0,
    label="f2",
);
plot!(posterior_plot, xp, μf₃′;
    linecolor=:blue,
    linewidth=2.0,
    label="f3",
);

# Plot observations
scatter!(posterior_plot, x₁, y₁;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);
scatter!(posterior_plot, x₃, y₃;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

display(posterior_plot);


text:
# Precompile linting some strings for performance

lintstr("""
x, y = 1, 2
for i in 1:10
    println(i)
end
if x == 0
    let z = 1
        println(x + z - y)
    end
elseif x == 1 && y == 2
    println("precompile")
end
""")


text:
function _string_to_bool(s::AbstractString)
    b = tryparse(Bool, s)
    if b isa Nothing
        return false
    else
        return b::Bool
    end
end

function _is_suppress_warning()
    s = get(ENV, "SUPPRESS_THREADINGUTILITIES_WARNING", "")
    b = _string_to_bool(s)::Bool
    return b
end

function _is_julia_exclusive()
    s = strip(get(ENV, "JULIA_EXCLUSIVE", "0"))
    # From https://github.com/JuliaLang/julia/blob/efc986003356be5daea6ccd9e38008de961c18aa/doc/src/manual/environment-variables.md:
    # If set to anything besides 0, then Julia's thread policy is consistent
    # with running on a dedicated machine: the master thread is on proc 0, and
    # threads are affinitized. Otherwise, Julia lets the operating system handle
    # thread policy.
    i = tryparse(Int, s)
    if isempty(s)
        return false
    elseif s == "0" || i == 0
        return false
    else
        return true
    end
end

function _print_exclusivity_warning()
    is_julia_exclusive = _is_julia_exclusive()::Bool
    if !is_julia_exclusive
        suppress_warning = _is_suppress_warning()::Bool
        if !suppress_warning
            msg = string(
                "The JULIA_EXCLUSIVE environment variable is not set to 1. ",
                "Therefore, the kernel is allowed to move threads to different ",
                "cores. We recommend that you set the JULIA_EXCLUSIVE ",
                "environment variable to 1. To suppress this warning, set the ",
                "SUPPRESS_THREADINGUTILITIES_WARNING environment ",
                "variable to 1.",
            )
            @debug(msg)
        end
    end
    return nothing
end


text:
"""
    struct TsptwReward <: AbstractReward end
TsptwReward is one of the already implemented rewards of SeaPearl.jl. A user can use it directly.
This reward is adapted to the tsptw problem and is inspired from the one used by Quentin Cappart in
his recent paper: Combining RL & CP for Combinatorial Optimization, https://arxiv.org/pdf/2006.01610.pdf.
"""
mutable struct TsptwReward <: AbstractReward
    value::Float32
    positiver::Float32
    normalizer::Float32
    max_dist::Float32
end

"""
    TsptwReward(model::CPModel)
Initialize a TsptwReward instance from a CPModel. The value is set to 0.
The positiver and normalizer are inspired from Quentin's paper.
Suggestions for the positiver:
    Mathematically exact and tighter than n * max(dist) :
n_biggest = partialsort(vec(dist), 1:n, rev = true)
n_minus_one_smallest = partialsort(vec(dist), 1:n-1)
positiver = sum(n_biggest) - sum(n_minus_one_smallest)
"""
function TsptwReward(model::CPModel)
    dist = nothing
    for constraint in model.constraints
        if isnothing(dist) && isa(constraint, Element2D) && size(constraint.matrix, 2) > 1
            dist = constraint.matrix
        end
    end
    n = size(dist, 1)
    max_dist = Float32(Base.maximum(dist))
    positiver = Float32(1 + (2^0.5) * n * max_dist)
    normalizer = positiver ^ (-1)
    TsptwReward(0, positiver, normalizer, max_dist)
end

"""
set_reward!(::StepPhase, lh::LearnedHeuristic{SR, R::TsptwReward, A}, model::CPModel, symbol::Union{Nothing, Symbol})
Change the "current_reward" attribute of the LearnedHeuristic at the StepPhase.
"""
function set_reward!(::Type{StepPhase}, lh::LearnedHeuristic{SR, TsptwReward, A}, model::CPModel, symbol::Union{Nothing, Symbol}) where {
    SR <: AbstractStateRepresentation,
    A <: ActionOutput
}
end

"""
    set_reward!(::DecisionPhase, lh::LearnedHeuristic{TsptwReward, O}, model::CPModel)
Change the current reward at the DecisionPhase. This is called right before making the next decision, so you know you have the very last state before the new decision
and every computation like fixPoints and backtracking has been done.
"""
function set_reward!(::Type{DecisionPhase}, lh::LearnedHeuristic{SR, TsptwReward, A}, model::CPModel) where {
    SR <: AbstractStateRepresentation,
    A <: ActionOutput
}
    n = size(lh.current_state.dist, 1)
    if !isnothing(model.statistics.lastVar)
        x = model.statistics.lastVar
        s = x.id
        current = parse(Int,split(x.id,'_')[2])
        if isbound(model.variables["a_"*string(current)])
            a_i = assignedValue(model.variables["a_"*string(current)])
            v_i = assignedValue(model.variables["v_"*string(current)])
            last_dist = lh.current_state.dist[v_i, a_i] * lh.reward.max_dist
            lh.reward.value += lh.reward.normalizer* (lh.reward.positiver - last_dist)
        end
    else lh.reward.value+=0
    end
end

"""
    set_reward!(::EndingPhase, lh::LearnedHeuristic{TsptwReward, A}, model::CPModel, symbol::Union{Nothing, Symbol})
Increment the current reward at the EndingPhase. Called when the search is finished by an optimality proof or by a limit in term of nodes or
in terms of number of solution. This is useful to add some general results to the reward like the number of ndoes visited during the episode for instance.
"""
function set_reward!(::Type{EndingPhase}, lh::LearnedHeuristic{SR, TsptwReward, A}, model::CPModel, symbol::Union{Nothing, Symbol}) where {
    SR <: AbstractStateRepresentation,
    A <: ActionOutput
}
    if symbol == :Feasible || symbol == :FoundSolution
        lh.reward.value += 20
    elseif symbol == :Infeasible
        lh.reward.value -= 10/(lh.search_metrics.total_steps+1)
    end
end


text:
export HomotopyWithCache

struct HomotopyWithCache{H<:AbstractHomotopy, C<:AbstractHomotopyCache} <: AbstractHomotopy
    homotopy::H
    cache::C
end

HomotopyWithCache(H::AbstractHomotopy, x, t) = HomotopyWithCache(H, cache(H, x, t))

cache(::HomotopyWithCache, x, t) = HomotopyNullCache()

(H::HomotopyWithCache)(x, t) = evaluate(H, x, t)

Base.size(H::HomotopyWithCache) = size(H.homotopy)
evaluate(H::HomotopyWithCache, x, t, ::HNC=HNC()) = evaluate(H.homotopy, x, t, H.cache)
evaluate!(u, H::HomotopyWithCache, x, t, ::HNC=HNC()) = evaluate!(u, H.homotopy, x, t, H.cache)
jacobian(H::HomotopyWithCache, x, t, ::HNC=HNC()) = jacobian(H.homotopy, x, t, H.cache)
jacobian!(U, H::HomotopyWithCache, x, t, ::HNC=HNC()) = jacobian!(U, H.homotopy, x, t, H.cache)
dt(H::HomotopyWithCache, x, t, ::HNC=HNC()) = dt(H.homotopy, x, t, H.cache)
dt!(u, H::HomotopyWithCache, x, t, ::HNC=HNC()) = dt!(u, H.homotopy, x, t, H.cache)
jacobian_and_dt(H::HomotopyWithCache, x, t, ::HNC=HNC()) = jacobian_and_dt(H.homotopy, x, t, H.cache)
jacobian_and_dt!(U, u, H::HomotopyWithCache, x, t, ::HNC=HNC()) = jacobian_and_dt!(U, u, H.homotopy, x, t, H.cache)
evaluate_and_jacobian(H::HomotopyWithCache, x, t, ::HNC=HNC()) = evaluate_and_jacobian(H.homotopy, x, t, H.cache)
evaluate_and_jacobian!(u, U, H::HomotopyWithCache, x, t, ::HNC=HNC()) = evaluate_and_jacobian!(u, U, H.homotopy, x, t, H.cache)
basehomotopy(H::HomotopyWithCache) = H.homotopy


text:
using SamplingHS3

# inverted egg crate function
f_ec = (x::Vector{Float64}) -> ( sin(x[1])^2 + sin(x[2])^2 ) / ( 0.1 + x[1]*x[1] + x[2]*x[2] )^(1/2.5)

# polytope
xlim = [-20.,20.]; ylim = [-20.,20.]
bounding_box_G =[eye(2);-eye(2)]
bounding_box_h =[-xlim[1];xlim[2];-ylim[1];ylim[2]]

x0 = [5.0;5.0]

# non-adaptive mcmc sampling using hr and mh:
conf_only_hr = SamplingHS3.DecorrelationConfig( 10000 , 1.0*eye(2) , 0.5 , 0.5 )
f_density    = f_ec
# this is quite slow.. most time is spent for the generation of the multivar. normal random numbers..
tic()
mcmc_result_01 = SamplingHS3.run_mcmc( x0 , f_density(x0) , f_density , (x) -> 0 , 1.0 , bounding_box_G  , bounding_box_h , conf_only_hr )
toc()

# run with adaptive mcmc
tic()
config_mcmc = SamplingHS3.MCMCConfig() # default config
mcmc_result_adaptive_01 = SamplingHS3.run_adaptive_mcmc( 10000 , x0 , f_density(x0) , f_density , (x) -> 0 , 1.0 , bounding_box_G  , bounding_box_h , config_mcmc )
toc()

# Plot results

using Plots, GR
Plots.gr()
Plots.plot(layout=(2,3))
z_sampled_1 = SamplingHS3.bin2d( mcmc_result_01.dd.xx , 40 , 40 )
z_sampled_2 = SamplingHS3.bin2d( mcmc_result_adaptive_01.dd.xx , 40 , 40 )
layout_fancy = @layout [a{0.33w} b{0.33w} c{0.33w}; d{0.5h}]
Plots.plot(layout=layout_fancy)
Plots.plot!( linspace(xlim[1],xlim[2],40), linspace(ylim[1],ylim[2],40) , f_likelihood_2p , subplot=1 , st = [:contourf])
Plots.plot!( z_sampled_1[2], [z_sampled_1[3]] , z_sampled_1[1][:]  , subplot=2 , st = [:contourf] , xlim=[-20;20] , ylim=[-20;20])
Plots.plot!( z_sampled_2[2], [z_sampled_2[3]] , z_sampled_2[1][:]  , subplot=3 , st = [:contourf] , xlim=[-20;20] , ylim=[-20;20])
# evolution of estimated covariance matrix:
ecov_steps = sort(collect(keys(mcmc_result_adaptive_01.ecov_regularized)))
get_ecov_entries = (i::Int) -> map( x -> mcmc_result_adaptive_01.ecov_regularized[x][i] , sort(ecov_steps) )
Plots.plot!( ecov_steps , get_ecov_entries(1) ,  subplot=4 , lab="a")
Plots.plot!( ecov_steps , get_ecov_entries(2) ,  subplot=4 , lab="b/c")
#Plots.plot!( ecov_steps , get_ecov_entries(3) ,  subplot=4 , lab="b/c")
Plots.plot!( ecov_steps , get_ecov_entries(4) ,  subplot=4 , lab="d")
Plots.title!( "Estimated Covariance Matrix" , subplot=4 )

#Plots.plot!( linspace(xlim[1],xlim[2],40), linspace(ylim[1],ylim[2],40) , z_sampled[:]    , subplot=2 , st = [:contourf])
abc = Plots.savefig("C:\\Temp\\non_adaptive_vs_adaptive_01.png" )


text:
using GeneticsMakie
using Documenter
# using other stuff

DocMeta.setdocmeta!(GeneticsMakie, :DocTestSetup, :(using GeneticsMakie); recursive=true)

makedocs(;
    modules=[GeneticsMakie],
    authors="Minsoo <mmkim1210@gmail.com> and contributors",
    repo="https://github.com/mmkim1210/GeneticsMakie.jl/blob/{commit}{path}#{line}",
    sitename="GeneticsMakie.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mmkim1210.github.io/GeneticsMakie.jl",
        assets=String[],
        sidebar_sitename = false,
    ),
    pages=[
        "Home" => "index.md",
        "Examples" =>
            ["Parsing GENCODE" => "examples/gtf.md",
            "Munging summary statistics" => "examples/summary.md",
            "Plotting genes" => "examples/genes.md",
            "Plotting isoforms" => "examples/isoforms.md",
            "Plotting LocusZoom" => "examples/locus.md",
            "Plotting GWAS" => "examples/gwas.md",
            "Plotting TWAS" => "examples/twas.md"],
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/mmkim1210/GeneticsMakie.jl",
)


text:
using Colors

"""
This function is a utility function that takes an agent and overwrites its color if it is of type :T (target) and gives it a grey color for display purposes.
"""
function PlotABM_ColorUtil(a::AbstractAgent)
    if a.type == :A || a.type == :O
        return a.color
    else
        return "#ffffff"
    end
end

"""
This function is a utility function for assigning agent display type based on agent type.
"""
function PlotABM_ShapeUtil(a::AbstractAgent)
    # potential options here:
    # https://gr-framework.org/julia-gr.html (search for "marker type")
    if a.type == :A
        return :circle
    elseif a.type == :O
        return :circle
    elseif a.type == :T
        return :circle
    else
        return :circle
    end
end

"""
This function is a utility function for setting agent plot size - note that this is for display purposes only and does not impact calculations involving agent radius.
"""
function PlotABM_RadiusUtil(a::AbstractAgent)
    # this is for display purposes only and does not impact FMP algorithm results
    # the scaling values are empirically selected
    # the object scale is based on the agent scaling

    # magic number appears to be scaling factor for plotting
    #   ex, an agent/object with radius=1 place in center of SS
    #   would take up the entire state space. This was empirically
    #   tested. The 1/a.SSdims is also tested and works well.

    SS_scale = 380*1/minimum(a.SSdims)  # technically a.SSdims[1] and [2] should be equal but just in case.

    if a.type == :O
        return a.radius*SS_scale
    else
        return a.radius*SS_scale
    end
end

"""
This function is a scheduler to determine draw order of agents. Draw order (left to right) is :T, :O, :A
"""
function PlotABM_Scheduler(model::ABM)

    # init blank lists
    agent_list = []
    object_list = []
    target_list = []
    for agent in values(model.agents)
        if agent.type == :A
            append!(agent_list, agent.id)
        elseif agent.type == :T
            append!(target_list, agent.id)
        elseif agent.type == :O
            append!(object_list, agent.id)
        end
    end

    # make composite list [targets, objects, agents]
    draw_order = []
    append!(draw_order, target_list)
    append!(draw_order, object_list)
    append!(draw_order, agent_list)

    return draw_order
end

"""
This function is a utility function for coloring agents.
"""
function AgentInitColor(i, num_agents)
    color_range = range(HSV(0,1,1), stop=HSV(-360,1,1), length=num_agents)
    agent_color = color_range[i]
    return string("#", hex(agent_color))

end

max_stars_repo_name: rmcsqrd/antfarm

max_issues_repo_name: rmcsqrd/antfarm

max_forks_repo_name: rmcsqrd/antfarm


text:
using Bokeh

autoopen(true)

x = linspace(0, 2pi)
y = sin(x)
plot(x, y, title="Hello World")


text:
#! /usr/local/julia/bin/julia
# Julia program to show a CSV file as a Julia dataframe as text.
using DataFrames
dfe_neutral(x) = 1.0
if length(ARGS) == 0
  simname = "../experiments/examples/in_example1"
else
  simname = ARGS[1]
end
include("$(simname).jl")
println("simname: ",simname)

#simname = "examples/in_example2"
df = readtable("$(simname).csv", makefactors=true, allowcomments=true)
println("df: ",df)

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using NGWP, LightGraphs, MTSG, Plots, LaTeXStrings

## Build Graph
N = 16; G = path_graph(N)
X = zeros(N,2); X[:, 1] = 1:N
L = Matrix(laplacian_matrix(G))
𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1,:])'
W = 1.0 * adjacency_matrix(G)

## preliminaries for PC-NGWP
Gstar_Sig = GraphSig(W)
G_Sig = GraphSig(W, xy = X)
GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
GP_primal = pairclustering(𝚽, GP_dual)

## Plots
function path_pc_plot!(W, X; ind = 1:size(X, 1), c = :teal)
    # plot(size = (500, 50), framestyle = :none, xlim = [1, N+1.5], ylim = [-0.2, 0.3])
    gplot!(W, X, width = 1, color = :blue, style = :solid)
    scatter_gplot!(X; ms = 6, c = :teal)
    annotate!(X[:, 1], X[:, 2] .- 0.2, [text(string(n), :top, 7) for n = 1:N])
    scatter_gplot!(X[ind, :]; ms = 6, c = c)
    plt = plot!(ratio = :auto)
    return plt
end


text:
#!/usr/bin/env julia

using Pkg
using MPI
using Documenter
using ICI_simulateur: Parallelisation, City, Population, Simulateur, Save

######################## Calcul du temps reel ##################################
deb = time() ;
################################################################################


text:
using LinSpaces2

function bench_fill(N, L)
    B = zeros(Float64, length(L))
    for i=1:N
        for j=1:length(L)
            B[j] = L[j]
        end
    end
end

function bench_collect(N, L)
    for i=1:N
        collect(L)
    end
end

function bench_sum(N, L)
    s = 0.0
    for i=1:N
        for j=1:length(L)
            s += L[j]
        end
        s /= length(L)
    end
    s
end

function bench_simd_sum(L)
    s = zero(eltype(L))
    @inbounds @simd for x in L
        s += x
    end
    s
end

L1 = linspace(0.0, 100000.0, 100000)
L2 = linspace2(0.0, 100000.0, 100000)
N = 100

println("bench_fill")
bench_fill(1, L1)
bench_fill(1, L2)
@time bench_fill(N, L1)
@time bench_fill(N, L2)

println("bench_sum")
bench_sum(1, L1)
bench_sum(1, L2)
@time bench_sum(N, L1)
@time bench_sum(N, L2)

println("bench_collect")
bench_collect(1, L1)
bench_collect(1, L2)
@time bench_collect(N, L1)
@time bench_collect(N, L2)

println("bench_simd_sum")
bench_simd_sum(linspace(1,5,10^7))
bench_simd_sum(linspace2(1,5,10^7))
@time bench_simd_sum(linspace(1,5,10^7))
@time bench_simd_sum(linspace2(1,5,10^7))

function bench_mult_sum(N, L)
    a = first(L)
    b = last(L)
    len = length(L)
    mult = Float64(1/(N-1))
    s = 0.0
    for i=1:N
        for j=1:length(L)
            l = (len-j)*mult*a + (j-1)*mult*b
            s += l
        end
        s /= length(L)
    end
    s
end

function bench_div_sum(N, L)
    a = first(L)
    b = last(L)
    len = length(L)
    divisor = Float64(N-1.0)
    s = 0.0
    for i=1:N
        for j=1:length(L)
            l = ((len-j)*a + (j-1)*b)/divisor
            s += l
        end
        s /= length(L)
    end
    s
end

println("bench_div_sum & bench_mult_sum")
bench_div_sum(1, L1)
bench_mult_sum(1, L1)
@time bench_div_sum(N, L1)
@time bench_mult_sum(N, L1)


text:
using RvSpectMLBase
using Test

@testset "Constructors and bare-bones functions" begin
    include("chunk_of_spectrum.jl")
    include("theoretical_instrument.jl")
    include("spectra.jl")
    include("chunk_list.jl")
    include("spectra_timeseries.jl")
    #include("spectra_timeseries_common_wavelengths.jl")  # TODO: Add tests
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
import brainflow

brainflow.enable_dev_board_logger()

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(256, board_shim)
brainflow.release_session(board_shim)

max_stars_repo_name: kusumikakd/Brainflow

max_issues_repo_name: kusumikakd/Brainflow

max_forks_repo_name: kusumikakd/Brainflow


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UnstructuredMesh(points, connec)

An unstructured mesh with `points` and connectivities
`connec`. The i-th element of the mesh is lazily built
based on the connectivity list `connec[i]`.
"""
struct UnstructuredMesh{Dim,T,Connectivity} <: Mesh{Dim,T}
  # input fields
  points::Vector{Point{Dim,T}}
  connec::Vector{Connectivity}

  # state fields
  ranks::Vector{Int}
  elms::Vector{Int}

  function UnstructuredMesh{Dim,T,C}(points, connec) where {Dim,T,C}
    ranks = [paramdim(polytopetype(c)) for c in connec]
    elms  = [i for i in eachindex(ranks) if ranks[i] == Dim]
    new(points, connec, ranks, elms)
  end
end

function UnstructuredMesh(points, connec)
  p = first(points)
  Dim = embeddim(p)
  T = coordtype(p)
  C = eltype(connec)
  UnstructuredMesh{Dim,T,C}(points, connec)
end

==(m1::UnstructuredMesh, m2::UnstructuredMesh) =
  m1.points == m2.points && m1.connec == m2.connec

vertices(m::UnstructuredMesh) = m.points

function faces(m::UnstructuredMesh{Dim}, r) where {Dim}
  @assert 0 ≤ r ≤ Dim "invalid rank for mesh"
  ps, cs, rs = m.points, m.connec, m.ranks
  r == 0 && return ps
  (materialize(cs[i], ps) for i in eachindex(cs) if rs[i] == r)
end

# -----------------
# DOMAIN INTERFACE
# -----------------

Base.getindex(m::UnstructuredMesh, ind::Int) =
  materialize(m.connec[m.elms[ind]], m.points)

nelements(m::UnstructuredMesh) = length(m.elms)

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, ::MIME"text/plain", m::UnstructuredMesh{Dim,T}) where {Dim,T}
  nvert = length(m.points)
  nface = length(m.connec)
  println(io, m)
  println(io, "  $nvert vertices")
  println(io, _lines(m.points, "    "))
  println(io, "  $nface faces")
  print(  io, _lines(m.connec, "    "))
end

function _lines(vec, tab="  ")
  N = length(vec)
  I, J = N > 10 ? (5, N-4) : (N, N+1)
  lines = [["$(tab)└─$(vec[i])" for i in 1:I]
           (N > 10 ? ["$(tab)⋮"] : [])
           ["$(tab)└─$(vec[i])" for i in J:N]]
  join(lines, "\n")
end


text:
"""
Compute the partial density at the indicated ``k``-Point and return it (in Fourier space).
"""
function compute_partial_density(basis, kpt, ψk, occupation)
    @assert length(occupation) == size(ψk, 2)

    # Build the partial density for this k-Point
    ρk_real = similar(ψk[:, 1], basis.fft_size)
    ρk_real .= 0
    for (ist, ψik) in enumerate(eachcol(ψk))
        ψik_real = G_to_r(basis, kpt, ψik)
        ρk_real .+= occupation[ist] .* abs2.(ψik_real)
    end

    # Check sanity of the density (real, positive and normalized)
    T = real(eltype(ρk_real))
    real_checked(ρk_real)
    if all(occupation .> 0)
        minimum(real(ρk_real)) < 0 && @warn("Negative ρ detected",
                                            min_ρ=minimum(real(ρk_real)))
    end
    n_electrons = sum(ρk_real) * basis.model.unit_cell_volume / prod(basis.fft_size)
    if abs(n_electrons - sum(occupation)) > sqrt(eps(T))
        @warn("Mismatch in number of electrons", sum_ρ=n_electrons,
              sum_occupation=sum(occupation))
    end

    # FFT and return
    r_to_G(basis, ρk_real)
end

"""
    compute_density(basis::PlaneWaveBasis, ψ::AbstractVector, occupation::AbstractVector)

Compute the density for a wave function `ψ` discretized on the plane-wave grid `basis`,
where the individual k-Points are occupied according to `occupation`. `ψ` should
be one coefficient matrix per k-Point.
"""
@timing function compute_density(basis::PlaneWaveBasis, ψ::AbstractVector,
                                 occupation::AbstractVector)
    n_k = length(basis.kpoints)

    # Sanity checks
    @assert n_k == length(ψ)
    @assert n_k == length(occupation)
    for ik in 1:n_k
        @assert length(G_vectors(basis.kpoints[ik])) == size(ψ[ik], 1)
        @assert length(occupation[ik]) == size(ψ[ik], 2)
    end
    @assert n_k > 0

    # Allocate an accumulator for ρ in each thread
    ρaccus = [similar(ψ[1][:, 1], basis.fft_size) for ithread in 1:Threads.nthreads()]

    # TODO Better load balancing ... the workload per kpoint depends also on
    #      the number of symmetry operations. We know heuristically that the Gamma
    #      point (first k-Point) has least symmetry operations, so we will put
    #      some extra workload there if things do not break even
    kpt_per_thread = [ifelse(i <= n_k, [i], Vector{Int}()) for i in 1:Threads.nthreads()]
    if n_k >= Threads.nthreads()
        kblock = floor(Int, length(basis.kpoints) / Threads.nthreads())
        kpt_per_thread = [collect(1:length(basis.kpoints) - (Threads.nthreads() - 1) * kblock)]
        for ithread in 2:Threads.nthreads()
            push!(kpt_per_thread, kpt_per_thread[end][end] .+ collect(1:kblock))
        end
        @assert kpt_per_thread[end][end] == length(basis.kpoints)
    end

    Threads.@threads for (ikpts, ρaccu) in collect(zip(kpt_per_thread, ρaccus))
        ρaccu .= 0
        for ik in ikpts
            ρ_k = compute_partial_density(basis, basis.kpoints[ik], ψ[ik], occupation[ik])
            # accumulates all the symops of ρ_k into ρaccu
            accumulate_over_symops!(ρaccu, ρ_k, basis, basis.ksymops[ik])
        end
    end

    count = sum(length(basis.ksymops[ik]) for ik in 1:length(basis.kpoints))
    from_fourier(basis, sum(ρaccus) / count)
end


text:
"""
    pairwise_sqd(x)

Pairwise squared distances between each colum vector in `x`.
"""
function pairwise_sqd(x)
    n = size(x, 2)
    xixj = transpose(x) * x
    xsq = sum(x .^ 2; dims=1)
    return transpose(xsq) .+ xsq - 2xixj
end

"""
    pairwise_sqd(x, y)

Pairwise squared distances between each colum vector in `x` and `y`.
"""
function pairwise_sqd(x, y)
    nx = size(x, 2)
    ny = size(y, 2)
    xiyj = transpose(x) * y
    xsq = sum(x .^ 2; dims=1)
    ysq = sum(y .^ 2; dims=1)
    return transpose(xsq) .+ ysq - 2xiyj
end

gaussian_gram_by_pairwise_sqd(pdot, σ) = exp.(-pdot ./ 2(σ ^ 2))
gaussian_gram(x, σ) = gaussian_gram_by_pairwise_sqd(pairwise_sqd(x), σ)
gaussian_gram(x, y, σ) = gaussian_gram_by_pairwise_sqd(pairwise_sqd(x, y), σ)

abstract type AbstractMMD end

function estimate_ratio(mmd::AbstractMMD, pdot_dede, pdot_denu, σ)
    Kdede = gaussian_gram_by_pairwise_sqd(pdot_dede, σ)
    Kdenu = gaussian_gram_by_pairwise_sqd(pdot_denu, σ)
    return _estimate_ratio(mmd, Kdede, Kdenu)
end

"""
    estimate_ratio(mmd::AbstractMMD, x_de, x_nu; σs=nothing)

Estimate `p_nu(x_de) / p_de(x_de)` using mutiple `σ` in `σs`.
If `σs` is not provided, the median of pairwise distances will be used.
"""
function estimate_ratio(mmd::AbstractMMD, x_de, x_nu; σs=nothing)
    pdot_dede = pairwise_sqd(x_de)
    pdot_denu = pairwise_sqd(x_de, x_nu)

    if isnothing(σs)
        σ = sqrt(median([pdot_dede..., pdot_denu...]))
        @info "Automatically choose σ using the median of pairwise distances: $σ."
        σs = [σ]
    end

    r_de = mapreduce(σ -> estimate_ratio(mmd, pdot_dede, pdot_denu, σ), +, σs)

    return inv(length(σs)) * r_de
end

struct MMDAnalytical{T<:AbstractFloat, S, M} <: AbstractMMD
    ϵ::T
    function MMDAnalytical(ϵ::T=1f-3; method::Symbol=:solve) where {T<:Number}
        @assert method in (:solve, :inv)
        S = iszero(ϵ) ? Val{:false} : Val{:true}
        if T == Int
            ϵ = float(ϵ)
        end
        return new{typeof(ϵ), S, Val{method}}(ϵ)
    end
end

function adddiag(mmd::MMDAnalytical{T, Val{:true}, M}, Kdede) where {T, M}
    return Kdede + diagm(0 => mmd.ϵ * fill(one(T), size(Kdede, 1)))
end
adddiag(mmd::MMDAnalytical{T, Val{:false}, M}, Kdede) where {T, M} = Kdede

function _estimate_ratio(mmd::MMDAnalytical{T, S, Val{:inv}}, Kdede, Kdenu) where {T, S}
    n_de, n_nu = size(Kdenu)
    return convert(T, n_de / n_nu) * inv(adddiag(mmd, Kdede)) * Kdenu * fill(one(T), n_nu)
end

function _estimate_ratio(mmd::MMDAnalytical{T, S, Val{:solve}}, Kdede, Kdenu) where {T, S}
    n_de, n_nu = size(Kdenu)
    return convert(T, n_de / n_nu) * (adddiag(mmd, Kdede) \ sum(Kdenu; dims=2)[:,1])
end


text:
# ---
# title: 162. Find Peak Element
# id: problem162
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/find-peak-element/description/>
# hidden: true
# ---
#
# A peak element is an element that is strictly greater than its neighbors.
#
# Given an integer array `nums`, find a peak element, and return its index. If
# the array contains multiple peaks, return the index to **any of the peaks**.
#
# You may imagine that `nums[-1] = nums[n] = -∞`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,1]
#     Output: 2
#     Explanation: 3 is a peak element and your function should return the index number 2.
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,1,3,5,6,4]
#     Output: 5
#     Explanation: Your function can return either index number 1 where the peak element is 2, or index number 5 where the peak element is 6.
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 1000`
#   * `-231 <= nums[i] <= 231 - 1`
#   * `nums[i] != nums[i + 1]` for all valid `i`.
#
#
#
# **Follow up:**  Could you implement a solution with logarithmic complexity?
#
#
## @lc code=start
using LeetCode

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


text:
# require

function find_in_path(name::String)
    isabspath(name) && return name
    isfile(name) && return abspath(name)
    base = name
    if endswith(name,".jl")
        base = name[1:end-3]
    else
        name = string(base,".jl")
        isfile(name) && return abspath(name)
    end
    for prefix in [Pkg.dir(), LOAD_PATH]
        path = joinpath(prefix, name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, base, "src", name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, name, "src", name)
        isfile(path) && return abspath(path)
    end
    return nothing
end

find_in_node1_path(name) = myid()==1 ?
    find_in_path(name) : remotecall_fetch(1, find_in_path, name)

function find_source_file(file)
    (isabspath(file) || isfile(file)) && return file
    file2 = find_in_path(file)
    file2 != nothing && return file2
    file2 = "$JULIA_HOME/../share/julia/base/$file"
    isfile(file2) ? file2 : nothing
end

# Store list of files and their load time
package_list = (ByteString=>Float64)[]
# to synchronize multiple tasks trying to require something
package_locks = (ByteString=>Any)[]
require(fname::String) = require(bytestring(fname))
require(f::String, fs::String...) = (require(f); for x in fs require(x); end)

# only broadcast top-level (not nested) requires and reloads
toplevel_load = true

function require(name::String)
    path = find_in_node1_path(name)
    path == nothing && error("$name not found")

    if myid() == 1 && toplevel_load
        refs = { @spawnat p _require(path) for p in filter(x->x!=1, procs()) }
        _require(path)
        for r in refs; wait(r); end
    else
        _require(path)
    end
    nothing
end

function _require(path)
    global toplevel_load
    if haskey(package_list,path)
        wait(package_locks[path])
    else
        last = toplevel_load
        toplevel_load = false
        try
            reload_path(path)
        finally
            toplevel_load = last
        end
    end
end

function reload(name::String)
    global toplevel_load
    path = find_in_node1_path(name)
    path == nothing && error("$name not found")
    refs = nothing
    if myid() == 1 && toplevel_load
        refs = { @spawnat p reload_path(path) for p in filter(x->x!=1, procs()) }
    end
    last = toplevel_load
    toplevel_load = false
    try
        reload_path(path)
    finally
        toplevel_load = last
    end
    if refs !== nothing
        for r in refs; wait(r); end
    end
    nothing
end

# remote/parallel load

include_string(txt::String, fname::String) =
    ccall(:jl_load_file_string, Any, (Ptr{Uint8},Ptr{Uint8}), txt, fname)

include_string(txt::String) = include_string(txt, "string")

function source_path(default::Union(String,Nothing)="")
    t = current_task()
    while true
        s = t.storage
        if !is(s, nothing) && haskey(s, :SOURCE_PATH)
            return s[:SOURCE_PATH]
        end
        if is(t, t.parent)
            return default
        end
        t = t.parent
    end
end

macro __FILE__() source_path() end

function include_from_node1(path::String)
    prev = source_path(nothing)
    path = (prev == nothing) ? abspath(path) : joinpath(dirname(prev),path)
    tls = task_local_storage()
    tls[:SOURCE_PATH] = path
    local result
    try
        if myid()==1
            # sleep a bit to process file requests from other nodes
            nprocs()>1 && sleep(0.005)
            result = Core.include(path)
            nprocs()>1 && sleep(0.005)
        else
            result = include_string(remotecall_fetch(1, readall, path), path)
        end
    finally
        if prev == nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
    result
end

function reload_path(path::String)
    had = haskey(package_list, path)
    if !had
        package_locks[path] = RemoteRef()
    end
    package_list[path] = time()
    tls = task_local_storage()
    prev = pop!(tls, :SOURCE_PATH, nothing)
    try
        eval(Main, :(Base.include_from_node1($path)))
    catch e
        had || delete!(package_list, path)
        rethrow(e)
    finally
        if prev != nothing
            tls[:SOURCE_PATH] = prev
        end
    end
    if !isready(package_locks[path])
        put!(package_locks[path],nothing)
    end
    nothing
end

function evalfile(path::String, args::Vector{UTF8String}=UTF8String[])
    return eval(Module(:__anon__),
                Expr(:toplevel,
                     :(const ARGS = $args),
                     :(eval(x) = Core.eval(__anon__,x)),
                     :(eval(m,x) = Core.eval(m,x)),
                     :(include($path))))
end
evalfile(path::String, args::Vector) = evalfile(path, UTF8String[args...])

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function bestimmePixelnummer(MittelpunkteX1,MittelpunkteX2,Pixelmitten)

PixelNummer = zeros(size(MittelpunkteX1,1),size(MittelpunkteX1,2));
AbstandPixelmitte = zeros(size(Pixelmitten,1),1);
Sortiert = zeros(size(Pixelmitten,1),1);

for Spalte = 1:size(MittelpunkteX1,2)
    for Zeile = 1:size(MittelpunkteX1,1)-1
        if (MittelpunkteX2[Zeile,Spalte] == 0) && (MittelpunkteX2[Zeile+1,Spalte] == 0)
            PixelNummer[Zeile,Spalte] = 0;
        else
            for Pixelzeile = 1:size(Pixelmitten,1)
                AbstandPixelmitte[Pixelzeile] = sqrt((Pixelmitten[Pixelzeile,2]-MittelpunkteX1[Zeile,Spalte])^2+(Pixelmitten[Pixelzeile,3]-MittelpunkteX2[Zeile,Spalte])^2);
            end
            Sortiert= sortieren(AbstandPixelmitte,1,0); # suche Minimum
            Minimum = Sortiert[1,1];
            Ort = finde(AbstandPixelmitte, Minimum);
            x = Int(Ort[1,1]);
            PixelNummer[Zeile,Spalte] = x;
        end
    end
    for Zeile = size(MittelpunkteX1,1)
        if (MittelpunkteX2[Zeile,Spalte] == 0) && (MittelpunkteX2[Zeile-1,Spalte] == 0)
            PixelNummer[Zeile,Spalte] = 0;
        else
            for Pixelzeile = 1:size(Pixelmitten,1)
                AbstandPixelmitte[Pixelzeile] = sqrt((Pixelmitten[Pixelzeile,2]-MittelpunkteX1[Zeile,Spalte])^2+(Pixelmitten[Pixelzeile,3]-MittelpunkteX2[Zeile,Spalte])^2);
            end
            Sortiert= sortieren(AbstandPixelmitte,1,0); # suche Minimum
            Minimum = Sortiert[1,1];
            Ort = finde(AbstandPixelmitte, Minimum);
            x = Int(Ort[1,1]);
            PixelNummer[Zeile,Spalte] = x;
        end
    end
end
print("Pixelnummer: ");
show(stdout,"text/plain",Int.(PixelNummer));

return PixelNummer;

end


text:
using CLHMM,HMMBase,Serialization,Distributions, BenchmarkTools, Test,Random,MS_HMMBase,Profile,ProfileView
import StatsFuns:logsumexp, logaddexp
include("ref_fns.jl")
Random.seed!(1)
π = fill((1/6),6,6)
D = [Categorical(ones(4)/4), Categorical([.7,.05,.15,.1]),Categorical([.15,.35,.4,.1]), Categorical([.6,.15,.15,.1]),Categorical([.1,.4,.4,.1]), Categorical([.2,.2,.3,.3])]
hmm = HMM(π, D)
log_π = log.(hmm.π)

#function modstep(hmm::HMM{Univariate,Float64}, observations::Matrix{Int64}, obs_lengths::Vector{Int64})
    # O,T = size(observations);
    # a = log.(hmm.π); π0 = transpose(log.(hmm.π0))
    # N = length(hmm.D); Γ = length(hmm.D[1].support);
    # mask=observations.!=0
    # #INITIALIZATION
    # βoi_T = zeros(O,N); βoi_t = zeros(O,N) #log betas at T initialised as zeros
    # Eoγim_T = fill(-Inf,O,Γ,N,N); Eoγim_t = fill(-Inf,O,Γ,N,N)
    # @inbounds for m in 1:N, i in 1:N, γ in 1:Γ, o in 1:O
    #     observations[o, obs_lengths[o]] == γ && m == i && (Eoγim_T[o, γ, i, m] = 0)
    # end
    # Tijm_T = fill(-Inf,O,N,N,N); Tijm_t = fill(-Inf,O,N,N,N) #Ti,j(T,m) = 0 for all m; in logspace

#     #RECURRENCE
#     modsweep!(hmm,a,N,Γ,βoi_T,βoi_t,Tijm_T,Tijm_t,Eoγim_T,Eoγim_t,observations,mask,obs_lengths)

#     #TERMINATION
#     lls = CLHMM.llhs(hmm,observations[:,1])
#     α1om = lls .+ π0 #first position forward msgs
#     Toij = [logsumexp([lps(view(Tijm_T,o,i,j,m), view(α1om,o,m)) for m in 1:N]) for o in 1:O, i in 1:N, j in 1:N] #terminate Tijs with forward messages
#     Eoiγ=[logsumexp([lps(view(Eoγim_T,o,γ,i,m), view(α1om,o,m)) for m in 1:N]) for o in 1:O, i in 1:N, γ in 1:Γ] #terminate Eids with forward messages

#     #INTEGRATE ACROSS OBSERVATIONS AND SOLVE FOR NEW HMM PARAMS
#     obs_penalty=log(O) #broadcast subtraction to normalise log prob vals by obs number
#     #INITIAL STATE DIST
#     π0_o=α1om.+βoi_T.-logsumexp.(eachrow(α1om.+βoi_T)) #estimate π0 for each o
#     new_π0=logsumexp.(eachcol(π0_o)).-obs_penalty #sum over obs and normalise by number of obs
#     #TRANSITION MATRIX
#     a_int = Toij.-logsumexp.([view(Toij,o,i,:) for o in 1:O, i in 1:N])
#     new_a = logsumexp.([a_int[:,i,j] for i in 1:N, j in 1:N]).-obs_penalty
#     #EMISSION MATRIX
#     e_int=Eoiγ.-logsumexp.([view(Eoiγ,o,j,:) for o in 1:O, j in 1:N])
#     new_b=logsumexp.([view(e_int,:,j,γ) for γ in 1:Γ, j in 1:N]).-obs_penalty
#     new_D::Vector{Categorical}=[Categorical(exp.(new_b[:,i])) for i in 1:N]

#     println(new_a)

#     return typeof(hmm)(exp.(new_π0), exp.(new_a), new_D), lps([logsumexp(lps.(α1om[o,:], βoi_T[o,:])) for o in 1:O])
# end

function modsweep!(hmm::HMM{Univariate,Float64}, a::Matrix{Float64},N::Int64,Γ::Int64,βoiT::Matrix{Float64},βoit::Matrix{Float64},TijmT::Array{Float64},Tijmt::Array{Float64},EoγimT::Array{Float64},Eoγimt::Array{Float64}, observations::Matrix{Int64}, mask::BitMatrix, obs_lengths::Vector{Int64})
    @inbounds for t in maximum(obs_lengths)-1:-1:1
        lls = CLHMM.llhs(hmm,observations[:,t+1])
        omask = findall(mask[:,t+1])
        for m in 1:N
            βoit[omask,m] .= logsumexp.(eachrow((view(βoiT,omask,:).+view(lls,omask,:)).+transpose(view(a,m,:))))
            for j in 1:N, i in 1:N
                Tijmt[omask, i, j, m] .= logsumexp.(eachrow(lps.(view(TijmT,omask,i,j,:), view(lls,omask,:), transpose(view(a,m,:)))))
                i==m && (Tijmt[omask, i, j, m] .= logaddexp.(Tijmt[omask, i, j, m], (βoiT[omask,j].+ a[m,j].+ lls[omask,j])))
            end
            for i in 1:N, γ in 1:Γ
                Eoγimt[omask, γ, i, m] .= logsumexp.(eachrow(lps.(view(EoγimT,omask,γ,i,:),view(lls,omask,:),transpose(view(a,m,:)))))
                if i==m
                    symmask = findall(observations[:,t].==γ)
                    Eoγimt[symmask, γ, i, m] .= logaddexp.(Eoγimt[symmask, γ, i, m], βoit[symmask,m])
                end
            end
        end
        βoiT=copy(βoit); TijmT=copy(Tijmt); EoγimT = copy(Eoγimt);
    end
    return βoiT, TijmT, EoγimT
end

no_obs=500
observations=zeros(Int64,no_obs,1001)
obs_lengths=Vector{Int64}()
for o in 1:size(observations)[1]
    obsl=rand(100:1000)
    push!(obs_lengths,obsl)
    observations[o,1:obsl]=rand(1:4,obsl)
end

O,T = size(observations);
a = log.(hmm.π); π0 = transpose(log.(hmm.π0))
N = length(hmm.D); Γ = length(hmm.D[1].support);
mask=observations.!=0
#INITIALIZATION
βoi_T = zeros(O,N); βoi_t = zeros(O,N) #log betas at T initialised as zeros
Eoγim_T = fill(-Inf,O,Γ,N,N); Eoγim_t = fill(-Inf,O,Γ,N,N)
@inbounds for m in 1:N, i in 1:N, γ in 1:Γ, o in 1:O
    observations[o, obs_lengths[o]] == γ && m == i && (Eoγim_T[o, γ, i, m] = 0)
end
Tijm_T = fill(-Inf,O,N,N,N); Tijm_t = fill(-Inf,O,N,N,N) #Ti,j(T,m) = 0 for all m; in logspace

# linear_step(hmm,observations,obs_lengths)

# Profile.clear()

# @profile linear_step(hmm,observations,obs_lengths)
# ProfileView.view()

@info "test sweeps"
ob, oa, od = CLHMM.backwards_sweep!(hmm,a,N,Γ,copy(βoi_T),copy(βoi_t),copy(Tijm_T),copy(Tijm_t),copy(Eoγim_T),copy(Eoγim_t),observations,mask,obs_lengths)
nb, na, nd = modsweep!(hmm,a,N,Γ,copy(βoi_T),copy(βoi_t),copy(Tijm_T),copy(Tijm_t),copy(Eoγim_T),copy(Eoγim_t),observations,mask,obs_lengths)

@test ob==nb
@test oa==na
@test od==nd

@info "bench sweeps"
modsweep = median(@benchmark (modsweep!(hmm,a,N,Γ,copy(βoi_T),copy(βoi_t),copy(Tijm_T),copy(Tijm_t),copy(Eoγim_T),copy(Eoγim_t),observations,mask,obs_lengths)))
sweep = median(@benchmark ( CLHMM.backwards_sweep!(hmm,a,N,Γ,copy(βoi_T),copy(βoi_t),copy(Tijm_T),copy(Tijm_t),copy(Eoγim_T),copy(Eoγim_t),observations,mask,obs_lengths)))

display(judge(modsweep,sweep))

# new_hmm, new_logpobs = linear_step(hmm, observations, obs_lengths)
# println(new_hmm.π)
# mod_hmm, mod_pobs = modstep(hmm, observations, obs_lengths)

# modsweep = median(@benchmark (modstep(hmm, observations, obs_lengths)))
# sweep = median(@benchmark (linear_step(hmm, observations, obs_lengths)))
# display(judge(modsweep,sweep))
# old_hmm, old_logpobs = old_linear(hmm, Array(transpose(observations)),obs_lengths)
# MS_hmm, MS_logpobs = MS_HMMBase.mle_step(hmm, Array(transpose(observations)),obs_lengths)

# for n in fieldnames(typeof(new_hmm))
#     if n == :D
#         for (d, dist) in enumerate(new_hmm.D)
#             @test new_hmm.D[d].support==mod_hmm.D[d].support #/==MS_hmm.D[d].support
#             #@test isapprox(new_hmm.D[d].support,old_hmm.D[d].support)=
#             @test isapprox(new_hmm.D[d].p,mod_hmm.D[d].p)
#             #@test isapprox(new_hmm.D[d].p,MS_hmm.D[d].p)
#         end
#     else
#         @test isapprox(getfield(new_hmm,n), getfield(mod_hmm,n))
#         #@test isapprox(getfield(new_hmm,n), getfield(MS_hmm,n))
#     end
# end

# @test isapprox(new_logpobs,mod_pobs)
# @test isapprox(new_logpobs,MS_logpobs)
# @info "Judging new linear vs old linear"

# @info "Judging new linear vs old B-W"
# new = median(@benchmark (modstep($hmm, $observations, $obs_lengths)))
# ms = median(@benchmark (MS_HMMBase.mle_step($hmm, $Array(transpose(observations)), $obs_lengths)))
# display(judge(new,ms))
# @info "Judging old linear vs old B-W"
# old = median(@benchmark (linear_step($hmm, $observations, $obs_lengths)))
# display(judge(old,ms))
# display(judge(new,old))

# obsl=1500
# sorted_obs = zeros(Int64, 1000,obsl+1)
# sorted_obsl = [i for i in obsl:-1:obsl-size(sorted_obs)[1]]
# unsorted_obs = zeros(Int64, 1000,obsl+1)
# unsorted_obsl = Vector{Int64}()
# for o in 1:size(sorted_obs)[1]
#     sorted_obs[o,1:obsl]=rand(1:4,obsl)
#     unsorted_obs[rand(findall(iszero,unsorted_obs[:,1])), 1:obsl]=rand(1:4,obsl)
#     push!(unsorted_obsl,obsl)
#     global obsl-=1
# end

# O,T = size(unsorted_obs);
# a = log.(hmm.π); π0 = transpose(log.(hmm.π0))
# N = length(hmm.D); D = length(hmm.D[1].support);
# umask=unsorted_obs.!=0
# smask=sorted_obs.!=0
# #INITIALIZATION
# βoi_T = zeros(O,N); βoi_t = zeros(O,N) #log betas at T initialised as zeros
# Eoγim_T = fill(-Inf,O,N,N,D); Eoγim_t = fill(-Inf,O,N,N,D)
# mEoimd_T = fill(-Inf,O,D,N,N); mEoimd_t = fill(-Inf,O,D,N,N)
# @inbounds for m in 1:N, i in 1:N, γ in 1:D, o in 1:O
#     unsorted_obs[o, unsorted_obsl[o]] == γ && m == i && (Eoγim_T[o, i, m, γ] = 0)
# end
# Tijm_T = fill(-Inf,O,N,N,N); Tijm_t = fill(-Inf,O,N,N,N) #Ti,j(T,m) = 0 for all m; in logspace

# @info "Judging modsweep vs sweep"
# modsweep = median(@benchmark (modsweep!(hmm,a,N,D,βoi_T,βoi_t,Tijm_T,Tijm_t,mEoimd_T,mEoimd_t,sorted_obs,smask,sorted_obsl)))
# sweep = median(@benchmark (CLHMM.backwards_sweep!(hmm,a,N,D,βoi_T,βoi_t,Tijm_T,Tijm_t,Eoγim_T,Eoγim_t,unsorted_obs,umask,unsorted_obsl)))
# display(judge(modsweep,sweep))


text:
# fitting noisy data to an exponential model
model(x, p) = p[1]*exp(-x.*p[2])

# some example data
srand(12345)
xdata = linspace(0,10,20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))

fit = curve_fit(model, xdata, ydata, [0.5, 0.5])
@assert norm(fit.param - [1.0, 2.0]) < 0.05

# can also get error estimates on the fit parameters
errors = estimate_errors(fit)
@assert norm(errors - [0.017, 0.075]) < 0.01


text:
using OOPStyle, PyCall, Random
using Test

@testset "OOPStyle.jl" begin
    Random.seed!(0)
    v = rand(10)
    @test reshape(v, 2, 5) == v.reshape(2,5)
    @test length(v) == v.length()

    np = pyimport("numpy")
    np.random.seed(0)
    @test sprint(print, np.random.rand(10)) == "[0.548814, 0.715189, 0.602763, 0.544883, 0.423655, 0.645894, 0.437587, 0.891773, 0.963663, 0.383442]"
end


text:
using Random: shuffle, seed!

seed!(1729)

@testset "Numeric" begin
    @syms a::Integer b c d x::Real y::Number
    @test simplify(x - y) == x + -1*y
    @test simplify(x - sin(y)) == x + -1*sin(y)
    @test simplify(-sin(x)) == -sin(x)
    @test simplify(1 * x * 2) == 2 * x
    @test simplify(1 + x + 2) == 3 + x
    @test simplify(b*b) == b^2 # tests merge_repeats
    @test simplify((a*b)^c) == a^c * b^c

    @test simplify(1x + 2x) == 3x
    @test simplify(3x + 2x) == 5x

    @test simplify(a + b + (x * y) + c + 2 * (x * y) + d)     == (3 * x * y) + a + b + c + d
    @test simplify(a + b + 2 * (x * y) + c + 2 * (x * y) + d) == (4 * x * y) + a + b + c + d

    @test simplify(a * x^y * b * x^d) == (a * b * (x ^ (d + y)))

    @test simplify(a + b + 0*c + d) == a + b + d
    @test simplify(a * b * c^0 * d) == a * b * d
    @test simplify(a * b * 1*c * d) == a * b * c * d
end

@testset "Pythagorean Identities" begin
    @syms a::Integer x::Real y::Number

    @test simplify(cos(x)^2 + 1 + sin(x)^2) == 2
    @test simplify(cos(y)^2 + 1 + sin(y)^2) == 2
    @test simplify(sin(y)^2 + cos(y)^2 + 1) == 2

    @test simplify(1 + y + tan(x)^2) == sec(x)^2 + y
    @test simplify(1 + y + cot(x)^2) == csc(x)^2 + y
end

#@testset "2π Identities" begin
#    @syms a::Integer x::Real y::Number
#
#    @test simplify(cos(x + 2π + a)) == cos(a + x)
#    @test simplify(tan(x + 2π * a)) == tan(x)
#end

@testset "Depth" begin
    @syms x
    R = RuleSet([@rule(sin(~x) => cos(~x)),
                 @rule( ~x + 1 => ~x - 1)])
    @test R(sin(sin(sin(x + 1)))) == cos(cos(cos(x - 1)))
    @test R(sin(sin(sin(x + 1))), depth=2) == cos(cos(sin(x + 1)))
end

@testset "RuleRewriteError" begin
    pred(x) = error("Fail")

    @syms a b

    rs = RuleSet([@rule ~x + ~y::pred => ~x])
    @test_throws SymbolicUtils.RuleRewriteError rs(a+b)
    err = try rs(a+b) catch err; err; end
    @test sprint(io->Base.showerror(io, err)) == "Failed to apply rule ~x + ~(y::pred) => ~x on expression a + b"
end

@testset "timerwrite" begin
    @syms a b c d
    expr1 = foldr((x,y)->rand([*, /])(x,y), rand([a,b,c,d], 100))
    SymbolicUtils.@timerewrite simplify(expr1)
end

@testset "Shuffle Rules" begin
    @syms a::Integer b c d x::Real y::Number
    R1 = RuleSet(SymbolicUtils.SIMPLIFY_RULES)
    R2 = RuleSet(shuffle(R1.rules))
    simplify_shuffle_tester(ex) = (R2 ∘ R1)(ex) == (R1 ∘ R2)(ex)

    for ex ∈ [foldr((x,y)->rand([*, +, -, ^, /])(x,y),
                        rand([a, b, c, d, x, y, 0, 1.0, 2], 5))
              for _ ∈ 1:20]
        @test isequal((R2 ∘ R1)(ex), (R1 ∘ R2)(ex))
    end
end


text:
using Gadfly
using ElectronDisplay
using Electron
using VegaDatasets
using Test

@testset "ElectronDisplay" begin

p1 = plot(y=[1,2,3])

f = display(p1)

@test f isa Electron.Window

ElectronDisplay.CONFIG.single_window = true

f = display(p1)

@test f === ElectronDisplay._getglobalwindow()

p2 = plot(y=[1,2,3])
f2 = display(p2)

@test f2 === f   # Window is reused

@testset "smoke test: single_window=$single_window focus=$focus " for
        single_window in [false, true],
        focus in [false, true]
    ElectronDisplay.CONFIG.single_window = single_window
    ElectronDisplay.CONFIG.focus = focus
    @test electrondisplay(dataset("cars")) isa Electron.Window
    @test electrondisplay(@doc reduce) isa Electron.Window
end

end


text:
using Test
using CliOptions

include("testutils.jl")

@testset "print_usage()" begin
    @testset "example 1" begin
        spec = CliOptionSpec(
            Option("-w", "--window-function", help = "window function to use"),
            OptionGroup(
                Option("-x", help = "explanatory variables"),
                Option("-y", help = "objective variable"),
                name = "Input data",
            ),
            Option("-t", help = "dtype (float or int)"),
            program = "",
            onerror = error,
        )
        buf = IOBuffer()
        print_usage(buf, spec; verbose = false)
        usage = String(take!(buf))
        @test usage == """Usage: PROGRAM -w WINDOW_FUNCTION -x X -y Y -t T
                       """

        buf = IOBuffer()
        print_usage(buf, spec; verbose = true)
        usage = String(take!(buf))
        # print(usage)  # DEBUG
        @test usage == """
                       Usage: PROGRAM -w WINDOW_FUNCTION -x X -y Y -t T

                       Options:
                           -w, --window-function WINDOW_FUNCTION
                                                  window function to use

                         Input data:
                           -x X                   explanatory variables

                           -y Y                   objective variable

                           -t T                   dtype (float or int)

                       """
    end
end


text:
include("Linear.jl")
include("Tanh.jl")
include("SoftMax.jl")
include("CrossEntropy.jl")

include("utils.jl")

println("Loading data...")
(Xtr, ttr), (Xva, tva), (Xte, tte) = loadMNIST(length(ARGS) == 1 ? ARGS[1] : "./data/mnist", Float32)

# And define a few training hyperparameters
const epochs = 100
const batch_size = 600
const λ = 0.13
const Nhid = 200

# And a few short-hands for later.
const Ntr, Ndim = size(Xtr)
const Nva = size(Xva, 1)
const Ncls = size(ttr, 2)

# Define the network (just a softmax-regression here)
const net = [
    Linear{Float32}(Ndim, Nhid, σ=1/sqrt(mean([Ndim, Nhid]))),
    Tanh{Float32}(),
    Linear{Float32}(Nhid, Ncls, σ=0),
    SoftMax{Float32}(),
]
const crit = CrossEntropy{Float32}()

# Train for a hundred epochs.
println("Starting training...")
for ep=1:epochs
    # Fit the network.
    cost = 0
    for i=1:batch_size:Ntr
        Xb = Xtr[i:min(i+batch_size-1,Ntr),:]
        tb = ttr[i:min(i+batch_size-1,Ntr),:]
        cost += train_batch(net, crit, Xb, tb, λ=λ)
    end

    # Evaluate on the validation set.
    nerr = 0
    for i=1:batch_size:Nva
        Xb = Xva[i:min(i+batch_size-1,Nva),:]
        tb = tva[i:min(i+batch_size-1,Nva),:]
        # Compute predicted probabilities.
        py = forward(net, Xb)
        # Find the most probable class label.
        _, y = findmax(py, 2)
        _, t = findmax(tb, 2)
        # Count the mistakes.
        nerr += sum(y .!= t)
    end

    @printf("Errors at epoch %d: %d (%.2f%% accuracy, %.2fNLL)\n", ep, nerr, 100*(1 - nerr/Nva), cost/length(1:batch_size:Ntr))
end


text:
mutable struct Contract
  conId::Int
  symbol::String
  secType::String
  lastTradeDateOrContractMonth::String
  strike::Float64
  right::String
  multiplier::String
  exchange::String
  primaryExchange::String
  currency::String
  localSymbol::String
  tradingClass::String
  includeExpired::Bool
  secIdType::String
  secId::String
  comboLegsDescrip::String
  comboLegs::Vector{ComboLeg}
  deltaNeutralContract::Union{DeltaNeutralContract,Nothing}
end
Contract(; conId=        0,
           symbol=      ns,
           secType=     ns,
           exchange=    ns,
           currency=    ns,
           localSymbol= ns) = Contract(conId, symbol, secType, ns, 0.0, ns, ns,
                                       exchange, ns, currency, localSymbol, ns,
                                       false, ns, ns, ns, [], nothing)

mutable struct Order
  orderId::Int
  clientId::Int
  permId::Int
  action::String
  totalQuantity::Float64
  orderType::String
  lmtPrice::Union{Float64,Nothing}
  auxPrice::Union{Float64,Nothing}
  tif::String
  activeStartTime::String
  activeStopTime::String
  ocaGroup::String
  ocaType::Int
  orderRef::String
  transmit::Bool
  parentId::Int
  blockOrder::Bool
  sweepToFill::Bool
  displaySize::Union{Int,Nothing}
  triggerMethod::Int
  outsideRth::Bool
  hidden::Bool
  goodAfterTime::String
  goodTillDate::String
  rule80A::String
  allOrNone::Bool
  minQty::Union{Int,Nothing}
  percentOffset::Union{Float64,Nothing}
  overridePercentageConstraints::Bool
  trailStopPrice::Union{Float64,Nothing}
  trailingPercent::Union{Float64,Nothing}
  faGroup::String
  faProfile::String
  faMethod::String
  faPercentage::String
  openClose::String
  origin::Origin
  shortSaleSlot::Int
  designatedLocation::String
  exemptCode::Int
  discretionaryAmt::Float64
  optOutSmartRouting::Bool
  auctionStrategy::AuctionStrategy
  startingPrice::Union{Float64,Nothing}
  stockRefPrice::Union{Float64,Nothing}
  delta::Union{Float64,Nothing}
  stockRangeLower::Union{Float64,Nothing}
  stockRangeUpper::Union{Float64,Nothing}
  randomizeSize::Bool
  randomizePrice::Bool
  volatility::Union{Float64,Nothing}
  volatilityType::Union{Int,Nothing}
  deltaNeutralOrderType::String
  deltaNeutralAuxPrice::Union{Float64,Nothing}
  deltaNeutralConId::Int
  deltaNeutralSettlingFirm::String
  deltaNeutralClearingAccount::String
  deltaNeutralClearingIntent::String
  deltaNeutralOpenClose::String
  deltaNeutralShortSale::Bool
  deltaNeutralShortSaleSlot::Int
  deltaNeutralDesignatedLocation::String
  continuousUpdate::Bool
  referencePriceType::Union{Int,Nothing}
  basisPoints::Union{Float64,Nothing}
  basisPointsType::Union{Int,Nothing}
  scaleInitLevelSize::Union{Int,Nothing}
  scaleSubsLevelSize::Union{Int,Nothing}
  scalePriceIncrement::Union{Int,Nothing}
  scalePriceAdjustValue::Union{Float64,Nothing}
  scalePriceAdjustInterval::Union{Int,Nothing}
  scaleProfitOffset::Union{Float64,Nothing}
  scaleAutoReset::Bool
  scaleInitPosition::Union{Int,Nothing}
  scaleInitFillQty::Union{Int,Nothing}
  scaleRandomPercent::Bool
  scaleTable::String
  hedgeType::String
  hedgeParam::String
  account::String
  settlingFirm::String
  clearingAccount::String
  clearingIntent::String
  algoStrategy::String
  algoParams::NamedTuple
  smartComboRoutingParams::NamedTuple
  algoId::String
  whatIf::Bool
  notHeld::Bool
  solicited::Bool
  modelCode::String
  orderComboLegs::Vector{Float64}
  orderMiscOptions::NamedTuple
  referenceContractId::Int
  peggedChangeAmount::Float64
  isPeggedChangeAmountDecrease::Bool
  referenceChangeAmount::Float64
  referenceExchangeId::String
  adjustedOrderType::String
  triggerPrice::Union{Float64,Nothing}
  adjustedStopPrice::Union{Float64,Nothing}
  adjustedStopLimitPrice::Union{Float64,Nothing}
  adjustedTrailingAmount::Union{Float64,Nothing}
  adjustableTrailingUnit::Int
  lmtPriceOffset::Union{Float64,Nothing}
  conditions::Vector{AbstractCondition}
  conditionsCancelOrder::Bool
  conditionsIgnoreRth::Bool
  extOperator::String
  softDollarTier::SoftDollarTier
  cashQty::Union{Float64,Nothing}
  mifid2DecisionMaker::String
  mifid2DecisionAlgo::String
  mifid2ExecutionTrader::String
  mifid2ExecutionAlgo::String
  dontUseAutoPriceForHedge::Bool
  isOmsContainer::Bool
  discretionaryUpToLimitPrice::Bool
  autoCancelDate::String
  filledQuantity::Union{Float64,Nothing}
  refFuturesConId::Union{Int,Nothing}
  autoCancelParent::Bool
  shareholder::String
  imbalanceOnly::Bool
  routeMarketableToBbo::Bool
  parentPermId::Union{Int,Nothing}
  usePriceMgmtAlgo::Union{Bool,Nothing}
  duration::Union{Int,Nothing}
  postToAts::Union{Int,Nothing}
end
Order() = Order(0, 0, 0, ns, 0, ns, nothing, nothing, ns, ns, ns, ns, 0, ns, true, 0,
                false, false, nothing, 0, false, false, ns, ns, ns, false, nothing, nothing,
                false, nothing, nothing, ns, ns, ns, ns, ns, CUSTOMER, 0, ns, -1, 0,
                false, UNSET, fill(nothing, 5)..., false, false,
                nothing, nothing, ns, nothing, 0, ns, ns, ns, ns, false, 0, ns, false,
                fill(nothing, 9)..., false, nothing, nothing, false, fill(ns, 8)...,
                (;), (;), ns, false, false, false, ns, [], (;), 0, 0, false, 0, ns, ns,
                fill(nothing, 4)..., 0, nothing, [], false, false, ns, SoftDollarTier(),
                nothing, ns, ns, ns, ns, false, false, false, ns, nothing, nothing,
                false, ns, false, false, nothing, nothing, nothing, nothing)

mutable struct ScannerSubscription
  numberOfRows::Int
  instrument::String
  locationCode::String
  scanCode::String
  abovePrice::Union{Float64,Nothing}
  belowPrice::Union{Float64,Nothing}
  aboveVolume::Union{Int,Nothing}
  marketCapAbove::Union{Float64,Nothing}
  marketCapBelow::Union{Float64,Nothing}
  moodyRatingAbove::String
  moodyRatingBelow::String
  spRatingAbove::String
  spRatingBelow::String
  maturityDateAbove::String
  maturityDateBelow::String
  couponRateAbove::Union{Float64,Nothing}
  couponRateBelow::Union{Float64,Nothing}
  excludeConvertible::Bool
  averageOptionVolumeAbove::Union{Int,Nothing}
  scannerSettingPairs::String
  stockTypeFilter::String
end
ScannerSubscription() = ScannerSubscription(-1, ns, ns, ns, fill(nothing, 5)..., fill(ns, 6)...,
                                            nothing, nothing, false, nothing, ns, ns)


text:
using ITensors
using Test

@testset "Scalar ITensors" begin
  A = ITensor(2.4)
  @test storage(A) isa ITensors.Dense{Float64}
  @test ndims(A) == 0
  @test order(A) == 0
  @test A[] == 2.4
  @test A[1] == 2.4
  @test scalar(A) == 2.4
  @test ITensors.symmetrystyle(A) == ITensors.NonQN()

  A[] = 3.4
  @test ndims(A) == 0
  @test order(A) == 0
  @test A[] == 3.4
  @test A[1] == 3.4
  @test scalar(A) == 3.4
  @test ITensors.symmetrystyle(A) == ITensors.NonQN()

  A[1] = 4.4
  @test ndims(A) == 0
  @test order(A) == 0
  @test A[] == 4.4
  @test A[1] == 4.4
  @test scalar(A) == 4.4
  @test ITensors.symmetrystyle(A) == ITensors.NonQN()

  A = ITensor()
  @test storage(A) isa ITensors.EmptyStorage{ITensors.EmptyNumber}
  @test ndims(A) == 0
  @test order(A) == 0
  @test A[] == 0.0
  @test A[1] == 0.0
  @test scalar(A) == 0.0
  @test ITensors.symmetrystyle(A) == ITensors.NonQN()

  A = ITensor()
  @test storage(A) isa ITensors.EmptyStorage{ITensors.EmptyNumber}
  A[] = 3.4
  @test storage(A) isa ITensors.Dense{Float64}
  @test ndims(A) == 0
  @test order(A) == 0
  @test A[] == 3.4
  @test A[1] == 3.4
  @test scalar(A) == 3.4
  @test ITensors.symmetrystyle(A) == ITensors.NonQN()

  A = ITensor()
  @test storage(A) isa ITensors.EmptyStorage{ITensors.EmptyNumber}
  A[1] = 4.4
  @test storage(A) isa ITensors.Dense{Float64}
  @test ndims(A) == 0
  @test order(A) == 0
  @test A[] == 4.4
  @test A[1] == 4.4
  @test scalar(A) == 4.4
  @test ITensors.symmetrystyle(A) == ITensors.NonQN()

  x = 2.3
  ITensor(fill(x, ())) == ITensor(x)
  ITensor(fill(x, (1))) == ITensor(x)
  ITensor(fill(x, (1, 1))) == ITensor(x)
  ITensor(fill(x, (1, 1, 1))) == ITensor(x)
  @test_throws ErrorException ITensor(fill(x, (2, 2)))
end

nothing


text:
# ---------------- #
# Fitting routines #
# ---------------- #

# Tensor representation, single function method; calls function that computes coefficients below below
function get_coefs(basis::Basis, bs::BasisMatrix{Tensor}, y::Vector{T}) where T
    _get_coefs_deep(basis, bs, y)[:,1]
end

# Tensor representation, multiple function method; calls function that computes coefficients below
function get_coefs(basis::Basis, bs::BasisMatrix{Tensor}, y::Matrix{T}) where T
    _get_coefs_deep(basis, bs, y)
end

function _get_coefs_deep(basis::Basis, bs::BasisMatrix{Tensor}, y)
    if any(bs.order[1, :] .!= 0)
        error("invalid basis structure - first elements must be order 0")
    end
    to_kron = bs.vals[1, :]  # 68
    ckronxi(to_kron, y, ndims(basis):-1:1)  # 66
end

# convert to expanded and call the method below
function get_coefs(basis::Basis, bs::BasisMatrix{Direct}, y)
    get_coefs(basis, convert(Expanded, bs), y)
end

get_coefs(basis::Basis, bs::BasisMatrix{Expanded}, y) = bs.vals[1] \ y

# common checks to be run at the top of each funfit
function check_funfit(basis::Basis, x, y)
    m = size(y, 1)
    length(basis) > m && error("Can't be more basis funcs than points in y")
    return m
end

# get_coefs(::Basis, ::BasisMatrix, ::Array) does almost all the work for
# these methods
function funfitxy(basis::Basis, bs::BasisMatrix, y)
    check_funfit(basis, bs, y)
    c = get_coefs(basis, bs, y)
    c, bs
end

# use tensor form
function funfitxy(basis::Basis, x::TensorX, y)
    m = check_funfit(basis, x, y)

    bs = BasisMatrix(basis, Tensor(), x, 0)
    c = get_coefs(basis, bs, y)
    c, bs
end

function funfitxy(basis::Basis, x, y)
    # check input sizes
    m = check_funfit(basis, x, y)

    # additional check
    size(x, 1) != m && error("x and y are incompatible")

    bs = BasisMatrix(basis, Direct(), x, 0)
    c = get_coefs(basis, bs, y)
    c, bs
end

function funfitf(basis::Basis, f::Function, args...)
    X, xn = nodes(basis)
    y = f(X, args...)
    funfitxy(basis, xn, y)[1]
end

function Base.:\(b::Basis, y::AbstractArray)
    x123 = nodes(b)[2]
    funfitxy(b, x123, y)[1]
end

Base.:\(b::Basis, f::Function) = funfitf(b, f)

# ---------- #
# Evaluation #
# ---------- #

function _funeval(c, bs::BasisMatrix{Tensor}, order::AbstractMatrix{Int})  # funeval1
    kk, d = size(order)  # 95

    # 98 reverse the order of evaluation: B(d) × B(d-1) × ⋯ × B(1)
    order = reverse(order .+ (size(bs.vals, 1).*(0:d-1)' - bs.order .+ 1), dims=2)

    # 99
    nx = prod([size(bs.vals[1, j], 1) for j=1:d])

    _T = promote_type(eltype(c), eltype(bs))
    f = Array{_T,3}(undef, nx, size(c, 2), kk)  # 100

    for i=1:kk
        f[:, :, i] = ckronx(bs.vals, c, order[i, :])  # 102
    end
    f
end

function _funeval(c, bs::BasisMatrix{Direct}, order::AbstractMatrix{Int})  # funeval2
    kk, d = size(order)  # 95
    # 114 reverse the order of evaluation: B(d)xB(d-1)x...xB(1)
    order = reverse(order .+ (size(bs.vals, 1).*(0:d-1)' - bs.order.+1), dims=2)

    _T = promote_type(eltype(c), eltype(bs))
    f = Array{_T,3}(undef, size(bs.vals[1], 1), size(c, 2), kk)  # 116

    for i in 1:kk
        f[:, :, i] = cdprodx(bs.vals, c, order[i, :])  # 118
    end
    f
end

function _funeval(c, bs::BasisMatrix{Expanded}, order::AbstractMatrix{Int})  # funeval3
    nx = size(bs.vals[1], 1)
    kk = size(order, 1)

    _T = promote_type(eltype(c), eltype(bs))
    f = Array{_T,3}(undef, nx, size(c, 2), kk)
    for i=1:kk
        this_order = order[i, :]
        ind = something(findfirst(x->bs.order[x, :] == this_order, 1:kk), 0)
        if ind == 0
            msg = string("Requested order $(this_order) not in BasisMatrix ",
                         "with order $(bs.order)")
            error(msg)
        end
        f[:, :, i] = bs.vals[ind]*c  # 154
    end

    f
end

# 1d basis + x::Number + c::Mat => 1 point, many func ==> out 1d
funeval(c::AbstractMatrix, basis::Basis{1}, x::Real, order=0) =
    vec(funeval(c, basis, fill(x, 1, 1), order))

# 1d basis + x::Number + c::Vec => 1 point, 1 func ==> out scalar
funeval(c::AbstractVector, basis::Basis{1}, x::Real, order=0) =
    funeval(c, basis, fill(x, 1, 1), order)[1]

# 1d basis + x::Vec + c::Mat => manypoints, many func ==> out 2d
funeval(c::AbstractMatrix, basis::Basis{1}, x::AbstractVector{T}, order=0) where {T<:Number} =
    funeval(c, basis, x[:, :], order)

# 1d basis + x::Vec + c::Vec => manypoints, 1 func ==> out 1d
funeval(c::AbstractVector, basis::Basis{1}, x::AbstractVector{T}, order=0) where {T<:Number} =
    vec(funeval(c, basis, reshape(x, length(x), 1), order))

# N(>1)d basis + x::Vec + c::Vec ==> 1 point, 1 func ==> out scalar
funeval(c::AbstractVector, basis::Basis{N}, x::AbstractVector{T}, order=0) where {N,T<:Number} =
    funeval(c, basis, reshape(x, 1, N), order)[1]

# N(>1)d basis + x::Vec + c::Mat ==> 1 point, many func ==> out vec
funeval(c::AbstractMatrix, basis::Basis{N}, x::AbstractVector{T}, order=0) where {N,T<:Number} =
    vec(funeval(c, basis, reshape(x, 1, N), order))

function funeval(c, basis::Basis{N}, x::TensorX, order::Int=0) where N
    # check inputs
    size(x, 1) == N ||  error("x must have d=$N elements")

    if order != 0
        msg = string("passing order as integer only allowed for $(order=0).",
                     " Try calling the version where `order` is a matrix")
        error(msg)
    end

    _order = fill(0, 1, N)
    bs = BasisMatrix(SplineSparse, basis, Tensor(), x, _order)  # 67
    funeval(c, bs, _order)
end

function funeval(c, basis::Basis{N}, x::TensorX, _order::AbstractMatrix) where N
    # check inputs
    size(x, 1) == N ||  error("x must have d=$N elements")
    order = _check_order(N, _order)

    # construct tensor form
    bs = BasisMatrix(SparseMatrixCSC, basis, Tensor(), x, order)  # 67

    # pass to specialized method below
    return funeval(c, bs, order)
end

function funeval(c, basis::Basis{N}, x::AbstractMatrix, order::Int=0) where N
    # check inputs
    @boundscheck size(x, 2) == N || error("x must have d=$(N) columns")

    if order != 0
        msg = string("passing order as integer only allowed for $(order=0).",
                     " Try calling the version where `order` is a matrix")
        error(msg)
    end

    _order = fill(0, 1, N)
    bs = BasisMatrix(SplineSparse, basis, Direct(), x, _order)  # 67
    _out = funeval(c, bs, _order)

    # we only had one order, so we want to collapse the third dimension of _out
    return _out[:, :, 1]
end

function funeval(c, basis::Basis{N}, x::AbstractMatrix, _order::AbstractMatrix) where N
    # check that inputs are conformable
    @boundscheck size(x, 2) == N || error("x must have d=$(N) columns")  # 62
    order = _check_order(N, _order)

    # construct BasisMatrix in Direct for
    bs = BasisMatrix(SplineSparse, basis, Direct(), x, order)  # 67

    # pass of to specialized method below
    funeval(c, bs, order)
end

function funeval(c::AbstractVector, bs::BasisMatrix, order::AbstractMatrix{Int})
    _funeval(c, bs, order)[:, 1, :]
end

function funeval(c::AbstractMatrix, bs::BasisMatrix, order::AbstractMatrix{Int})
    _funeval(c, bs, order)
end

# default method
function funeval(c::AbstractVector, bs::BasisMatrix,
                 order::Vector{Int}=fill(0, size(bs.order, 2)))
    _funeval(c, bs, reshape(order, 1, length(order)))[:, 1, 1]
end

function funeval(c::AbstractMatrix, bs::BasisMatrix,
                 order::Vector{Int}=fill(1, size(bs.order, 2)))
    _funeval(c, bs, reshape(order, 1, length(order)))[:, :, 1]
end

# ------------------------------ #
# Convenience `Interpoland` type #
# ------------------------------ #

mutable struct Interpoland{TB<:Basis,TC<:AbstractArray,TBM<:BasisMatrix{Tensor}}
    basis::TB  # the basis -- can't change
    coefs::TC  # coefficients -- might change
    bmat::TBM  # BasisMatrix at nodes of `b` -- can't change
end

function Interpoland(basis::Basis, bs::BasisMatrix{Tensor}, y::AbstractArray)
    c = get_coefs(basis, bs, y)
    Interpoland(basis, c, bs)
end

# compute Tensor form and hand off to method above
function Interpoland(basis::Basis, y::AbstractArray)
    bs = BasisMatrix(basis, Tensor())
    Interpoland(basis, bs, y)
end

"""
Construct an Interpoland from a function.

The function must have the signature `f(::AbstractMatrix)::AbstractArray`
where each column of the input matrix is a vector of values along a single
dimension
"""
function Interpoland(basis::Basis, f::Function)
    x, xd = nodes(basis)
    y = f(x)
    bs = BasisMatrix(basis, Tensor(), xd)
    Interpoland(basis, bs, y)
end

Interpoland(p::BasisParams, f::Function) = Interpoland(Basis(p), f)

# let funeval take care of order and such. This just exists to make it so the
# user doesn't have to keep track of the coefficient vector
(itp::Interpoland)(x, order=0) = funeval(itp.coefs, itp.basis, x, order)

# now, given a new vector of `y` data we construct a new coefficient vector
function update_coefs!(interp::Interpoland, y::AbstractArray)
    # leverage the BasisMatrix we kept around
    c = funfitxy(interp.basis, interp.bmat, y)[1]
    copy!(interp.coefs, c)  # update c inplace b/c Interpoland is immutable
end

# similar for a function -- just hand off to above
update_coefs!(interp::Interpoland, f::Function) =
    update_coefs!(interp, f(nodes(interp.basis)[1]))

# alias update_coefs! to fit!
fit!(interp::Interpoland, y::AbstractArray) = update_coefs!(interp, y)
fit!(interp::Interpoland, f::Function) = update_coefs!(interp, f)

Base.show(io::IO, ::Interpoland{T,N,BST}) where {T,N,BST<:ABSR} =
    print(io, "$N dimensional interpoland")


text:
##### auxillary.jl
# Created by Alexander Dowling

################
# Hidden Auxiliary Functions

function checkJacobianScaling(m::Model, dd::DegenData, f=stdout)

	# Determine location of non-zeros. This is important as fixed variables
	# are still in the Jacobian
	nz_ = findall(dd.J .!= 0)

	# Determine location of smallest and largest elements
	s_ = nz_[argmin(abs.(dd.J[nz_]))]
	l_ = argmax(abs.(dd.J))

	println(f," ")
	println(f,"Smallest non-zero element in Jacobian = ",dd.J[s_])
	println(f,"Variable: ")
	printVariable(m, dd.jC[s_], f)
	println(f,"Equation: ")
	printEquation(m, dd, dd.iR[s_], f)
	printVariablesInEquation(m, dd, dd.iR[s_], true, true, f)

	println(f," ")
	println(f,"Largest non-zero element in Jacobian = ",dd.J[l_])
	println(f,"Variable: ")
	printVariable(m, dd.jC[l_], f)
	println(f,"Equation: ")
	printEquation(m, dd, dd.iR[l_], f)
	printVariablesInEquation(m, dd, dd.iR[l_], true, true,  f)
	println(f," ")

end

function checkActiveJacobianRows(m::Model, dd::DegenData, f=stdout)
	for i = 1:dd.nLambda
		if(sum(abs.(dd.J_active[i,:])) < 1E-4)
			println(f,"Warning!!! (Near) singular Jacobian row for ")
			printRows(m, dd, i, f)
			println(f," ")
		end
	end

	return nothing
end

function checkActiveJacobianEntries(m::Model, dd::DegenData, f=stdout)
	for i = 1:size(dd.J_active,1)
		nan_ = sum(isnan.(dd.J_active[i,:])) > 0
		inf_ = sum(isinf.(dd.J_active[i,:])) > 0
		if(nan_ || inf_)
			if(nan_)
				print(f,"NaN")
			end
			if(nan_ && inf_)
				print(f," and ")
			end
			if(inf_)
				print(f,"Inf")
			end
			println(f," in Jacobian")
			printJacobianForEquation(m, dd, f)
			println(f," ")
		end
	end

	return nothing
end

function printEquation(m::Model, dd::DegenData, i::Int64, f=stdout)

	nLin = MathProgBase.numlinconstr(m)
	nQuad = MathProgBase.numquadconstr(m)

	if(i <= nLin)
		#println(f,ConstraintRef{LinearConstraint}(m, i))
		println(f, m.linconstr[i])
	elseif(i <= nLin + nQuad)
		#println(f,ConstraintRef{QuadConstraint}(m, i - nLin))
		println(f, m.quadconstr[i - nLin])
	else
		println(f,MathProgBase.constr_expr(dd.d,i))
	end

	return nothing

end

#=
function printRow(m::Model, dd::DegenData, j::Int64, epsiActive::Float64, f=stdout)

	if(j <= length(dd.gMap))
		i = dd.gMap[j]

		printEquation(m, dd, i, f)

	else
		i = dd.bMap[j - length(dd.gMap)]

		printBound(m, i, epsiActive, f)

	end

	return nothing

end
=#

function printVariable(m::Model, dd::DegenData, i::Int64, f=stdout)

	v = Variable(m, dd.vMap[i])
	println(f,"x[",string(i),"] = ",getname(v))

	return nothing

end

function printVariable(m::Model, i::Int64, f=stdout)
	v = Variable(m, i)
	println(f,"x[",string(i),"] = ",getname(v))
	return nothing
end

function printJacobianForEquation(m::Model, dd::DegenData, r::Int64, f = stdout)

	printEquation(m, dd, f)

	if(dd.Jactive === nothing)
		# Use sparse tuple representation of full Jacobian

		k = findall(r == dd.iR)
		for i=1:length(k)
			J_ = dd.J[k[i]]
			if(isnan(J_) || J_ != 0)
				c = jC[k]
				v = Variable(m,c)
				println(f,"J = ",J_," for x[",c,"] = ",getname(v)," = ",dd.x[c])
			end
		end

	else
		# Sparse matrix representation of active Jacobian
		for c = 1:size(dd.Jactive,2)
			if(isnan(dd.Jactive[r,c]) || dd.Jactive[r,c] != 0)
				i = vMap[c]
				v = Variable(m, i)
				println(f,"J = ",Jactive[r,c]," for x[",i,"] = ",getname(v)," = ",dd.x[i])
			end
		end

	end

	return nothing

end

function printVariablesInEquation(m::Model, r::Int64)

	dd = DegenData(m)
	printVariablesInEquation(m, dd, r, true)

end

function printVariablesInEquation(m::Model, dd::DegenData, r::Int64, printVarNumber::Bool, printVarBounds::Bool=true, f=STDOUT)
	printVariablesInEquation(m, dd, Array([r]), printVarNumber, printVarBounds, f)
end

function printVariablesInEquation(m::Model, dd::DegenData, r::Array{Int64,1}, printVarNumber::Bool, printVarBounds::Bool=true, f=STDOUT)

	if(length(r) > 1)
		k = Array{Int64}(0)

		for i = 1:length(r)
			t = findall(dd.iR .== r[i])
			for j = 1:length(t)
				push!(k, t[j])
			end
		end

	else
		k = findall(dd.iR .== r)
	end

	if(length(unique(dd.jC[k])) < 50)
		for j in unique(dd.jC[k])
			v = Variable(m,j)

			if(printVarNumber)
				print(f,"x[",string(j),"] = ")
			end

			print(f,getname(v)," = ",dd.x[j])

			if(printVarBounds)
				print(f,"  \t [")
				lo = getlowerbound(v)
				up = getupperbound(v)

				if(abs(dd.x[j] - lo) < 1E-4)
					printstyled(color=:red,f, string(lo))
				else
					print(f,string(lo))
				end

				print(f,",")

				if(abs(dd.x[j] - up) < 1E-4)
					printstyled(color=:blue,f, string(up))
				else
					print(f,string(up))
				end

				print(f,"]")
			end

			print(f,"\n")

		end
	else
		println(f,"Over 50 variables in this equation... not printing variable values")
	end

	return nothing

end

function printIDS(ids::IrreducibleDegenerateSet, m::Model, dd::DegenData, ds::DegenSettings, f=stdout)

	epsiLambda = ds.epsiLambda
	epsiActive = ds.epsiActive

	if(ids.nElements > 0)

		println(f,"IDS ",string(ids.c),"...")

		nLambda = length(ids.lambda)

		r = findall(abs.(ids.lambda) .> epsiLambda)

		for j = r

			s = @sprintf "%.4f" ids.y[j]
			print(f,"y = ", s,", ")

			s = @sprintf "%.4f" ids.lambda[j]
			print(f,"l = ", s ,"  \t ")

			printRows(m, dd, j, f, lite=true)
		end

		println(f," ")

		println(f,"Involved variables: ")

		r2 = intersect(r, 1:length(dd.gMap))
		printVariablesInEquation(m, dd, dd.gMap[r2], true, true, f)

		println(f," ")

	else

		#println("The following "," does not contribute significantly to a IDS")
	end

	return nothing

end