text:
# This file was generated, do not modify it. # hide
train, test = partition(collect(eachindex(y)), 0.70, shuffle=true, rng=52)
XGBC = @load XGBoostClassifier
xgb_model = XGBC()

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PrivateEndpointConnectionProperties <: SwaggerModel
    privateEndpoint::Any # spec type: Union{ Nothing, PrivateEndpoint } # spec name: privateEndpoint
    privateLinkServiceConnectionState::Any # spec type: Union{ Nothing, PrivateLinkServiceConnectionState } # spec name: privateLinkServiceConnectionState
    provisioningState::Any # spec type: Union{ Nothing, PrivateEndpointConnectionProvisioningState } # spec name: provisioningState

    function PrivateEndpointConnectionProperties(;privateEndpoint=nothing, privateLinkServiceConnectionState=nothing, provisioningState=nothing)
        o = new()
        validate_property(PrivateEndpointConnectionProperties, Symbol("privateEndpoint"), privateEndpoint)
        setfield!(o, Symbol("privateEndpoint"), privateEndpoint)
        validate_property(PrivateEndpointConnectionProperties, Symbol("privateLinkServiceConnectionState"), privateLinkServiceConnectionState)
        setfield!(o, Symbol("privateLinkServiceConnectionState"), privateLinkServiceConnectionState)
        validate_property(PrivateEndpointConnectionProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type PrivateEndpointConnectionProperties

const _property_map_PrivateEndpointConnectionProperties = Dict{Symbol,Symbol}(Symbol("privateEndpoint")=>Symbol("privateEndpoint"), Symbol("privateLinkServiceConnectionState")=>Symbol("privateLinkServiceConnectionState"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_PrivateEndpointConnectionProperties = Dict{Symbol,String}(Symbol("privateEndpoint")=>"PrivateEndpoint", Symbol("privateLinkServiceConnectionState")=>"PrivateLinkServiceConnectionState", Symbol("provisioningState")=>"PrivateEndpointConnectionProvisioningState")
Base.propertynames(::Type{ PrivateEndpointConnectionProperties }) = collect(keys(_property_map_PrivateEndpointConnectionProperties))
Swagger.property_type(::Type{ PrivateEndpointConnectionProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PrivateEndpointConnectionProperties[name]))}
Swagger.field_name(::Type{ PrivateEndpointConnectionProperties }, property_name::Symbol) =  _property_map_PrivateEndpointConnectionProperties[property_name]

function check_required(o::PrivateEndpointConnectionProperties)
    (getproperty(o, Symbol("privateLinkServiceConnectionState")) === nothing) && (return false)
    true
end

function validate_property(::Type{ PrivateEndpointConnectionProperties }, name::Symbol, val)
end

ext: jl


text:
# May need to use build.jl to rebuild PyCall then
# restart Julia to make this work
using PyCall

sys = pyimport("sys")
qws = pyimport("PySide2.QtWidgets")

app = qws.QApplication([])
label = qws.QLabel("Qt for Python!")
label.show()
#sys.exit(app.exec_())


text:
# This version is focused in Julia 1.0 version
module Liga
	# Dependencies
	import Base.show
	import Base.reverse
	import Base.convert
	#   import OhMyREPL
	using SparseArrays,LinearAlgebra
	using OhMyREPL
	# New definitions
	export  BasisBlade,MultiVector,
	    GAAbstractType,Blade,
	    layout,layout_info,
	    tree,buildtree,
	    grade,grademinus,
	    gradeplus,gradeprojection,
	    bbgeometric,bbinner,bbouter,
	    geometric,scalar,inner,
	    outer,mvsum,mvdiff,
	    reverse,conjugate,⋅,
	    dual,inverse,OperationTable,
	    embedding,copy,isequal,==,
	    multivector,≈,magnitude,blade

	# Main files
	# Objects.jl produces two objects: BasisBlade and MultiVectors
	include("Objects.jl")
	# Layout.jl produces the enviroment to work. In current version
	# three enviroments are suported: GA, Conformal and Projetive
	include("Layout.jl")
	# OperationalTables.jl produces arrays of BasisBlade geometric,
	# inner and outer products. This is an important function because
	# this is the optimization of engine of Liga
	include("OperationalTables.jl")
	# LigaFunctions.jl contains all implemented functions like MultiVector
	# geometric product.
	include("Functions.jl")
	# LigaShow.jl contains some modifications to print BasisBlade and
	# MultiVectors.
	include("LigaShow.jl")
end


text:
# -----------------------------------------
# update [I - TV] operator based on current prices
# -----------------------------------------

function update_IminusTVp!(tmpv::DCDPTmpVars{T,M}, ddm::DynamicDrillModel{S,PF,M}, q0::AbstractVector) where {T,S,M<:Matrix,PF}
    IminusTVp = IminusTEVp(tmpv)
    ztrans = ztransition(ddm)
    β = discount(ddm)

    n = checksquare(IminusTVp)
    n == length(q0) == checksquare(ztrans) || throw(DimensionMismatch())

    @inbounds for j in OneTo(n)
        @simd for i in OneTo(n)
            x = -ztrans[i,j] * β * q0[j]
            IminusTVp[i,j] = i==j  ?  1+x  : x
        end
    end
end

function update_IminusTVp!(tmpv::DCDPTmpVars{T,M}, ddm::DynamicDrillModel{S,PF,M}, q0::AbstractVector) where {T,S,M<:SparseMatrixCSC,PF}
    IminusTVp = IminusTEVp(tmpv)
    ztrans = ztransition(ddm)

    n = checksquare(IminusTVp)
    n == length(q0) == checksquare(ztrans) || throw(DimensionMismatch())

    ztrans_rows = rowvals(ztrans)
    ztrans_vals = nonzeros(ztrans)
    IminusTVp_rows = rowvals( IminusTVp)
    IminusTVp_vals = nonzeros(IminusTVp)

    # consider eliminmating this check?
    length(ztrans_rows) == length(IminusTVp_rows) || throw(DimensionMismatch(
        "length(ztrans_rows) = $(length(ztrans_rows)) but length(IminusTVp_rows) = $(length(IminusTVp_rows))"
    ))

    # IminusTVp_rows == ztrans_rows || throw(error("IminusTVp rows not the same as ztrans"))

    fill!(IminusTVp, 0)
    @inbounds for j in OneTo(n)
        @simd for nzi in nzrange(IminusTVp, j)
            x = -ztrans_vals[nzi] * discount(ddm) * q0[j]
            i = IminusTVp_rows[nzi]
            IminusTVp_vals[nzi] = i==j  ?  1+x  :  x
        end
    end

end


text:
"""
    update!(∂, x, optim)

Computes parameter update step according to `optim` algorithm (e.g. ADAM), and gradient `∂`
of parameters `x`.
Overwrites `∂` with update step and returns it.
Note that this does not update parameters.
"""
function update!(∂::AbstractArray, x::AbstractArray, optim)
    ∂ .= Flux.Optimise.apply!(optim, x, ∂)
    return ∂
end

function update!(∂::NamedTuple, rbm::RBM, optim)
    update!(∂.visible, rbm.visible, optim)
    update!(∂.hidden, rbm.hidden, optim)
    update!(∂.w, rbm.w, optim)
    return ∂
end

function update!(∂::NamedTuple, l::Union{Binary,Spin,Potts}, optim)
    update!(∂.θ, l.θ, optim)
    return ∂
end

function update!(∂::NamedTuple, l::Union{Gaussian,ReLU}, optim)
    update!(∂.θ, l.θ, optim)
    update!(∂.γ, l.γ, optim)
    return ∂
end

function update!(∂::NamedTuple, l::dReLU, optim)
    update!(∂.θp, l.θp, optim)
    update!(∂.θn, l.θn, optim)
    update!(∂.γp, l.γp, optim)
    update!(∂.γn, l.γn, optim)
    return ∂
end

function update!(∂::NamedTuple, l::pReLU, optim)
    update!(∂.θ, l.θ, optim)
    update!(∂.γ, l.γ, optim)
    update!(∂.Δ, l.Δ, optim)
    update!(∂.η, l.η, optim)
    return ∂
end

function update!(∂::NamedTuple, l::xReLU, optim)
    update!(∂.θ, l.θ, optim)
    update!(∂.γ, l.γ, optim)
    update!(∂.Δ, l.Δ, optim)
    update!(∂.ξ, l.ξ, optim)
    return ∂
end

"""
    update!(x, ∂)

Updates parameters according to steps `∂`.
"""
function update!(x::AbstractArray, ∂::AbstractArray)
    x .-= ∂
    return x
end

function update!(rbm::RBM, ∂::NamedTuple)
    update!(rbm.visible, ∂.visible)
    update!(rbm.hidden, ∂.hidden)
    update!(rbm.w, ∂.w)
    return rbm
end

function update!(l::Union{Binary,Spin,Potts}, ∂::NamedTuple)
    update!(l.θ, ∂.θ)
    return l
end

function update!(l::Union{Gaussian,ReLU}, ∂::NamedTuple)
    update!(l.θ, ∂.θ)
    update!(l.γ, ∂.γ)
    return l
end

function update!(l::dReLU, ∂::NamedTuple)
    update!(l.θp, ∂.θp)
    update!(l.θn, ∂.θn)
    update!(l.γp, ∂.γp)
    update!(l.γn, ∂.γn)
    return l
end

function update!(l::pReLU, ∂::NamedTuple)
    update!(l.θ, ∂.θ)
    update!(l.γ, ∂.γ)
    update!(l.Δ, ∂.Δ)
    update!(l.η, ∂.η)
    return l
end

function update!(l::xReLU, ∂::NamedTuple)
    update!(l.θ, ∂.θ)
    update!(l.γ, ∂.γ)
    update!(l.Δ, ∂.Δ)
    update!(l.ξ, ∂.ξ)
    return l
end

update!(x::Any, ::Nothing) = x

"""
    gradnorms(∂)

Computes gradient norms.
"""
gradnorms(∂::AbstractArray) = norm(∂)
gradnorms(∂::NamedTuple) = map(gradnorms, ∂)
gradnorms(::Nothing) = false

"""
    default_optimizer(nsamples, batchsize, epochs; decay_after, decay_final, clip, optim)

Sane defaults for optimizer. All keyword arguments are optional and have default settings.
Trains with constant learning rate for first `decay_after` period, then decays
learning rate exponentially every epoch, starting after
`decay_after` of training time, until reaching `decay_final` at the last epoch.
Clips gradients by `clip`.
"""
function default_optimizer(
    nsamples::Int, batchsize::Int, epochs::Int;
    decay_final::Real = 0.01, decay_after::Real = 0.5, clip = 1,
    optim = ADAM(5e-3, (0, 0.99), 1e-3)
)
    # Defaults from https://github.com/jertubiana/PGM
    # See also 10.1016/j.cels.2020.11.005 (search RMSprop)
    steps_per_epoch = minibatch_count(nsamples; batchsize)
    startepoch = round(Int, epochs * decay_after)
    start = startepoch * steps_per_epoch - 1
    decay = decay_final^(1/(epochs - startepoch))
    return Flux.Optimise.Optimiser(
        Flux.ClipValue(clip),
        optim,
        Flux.ExpDecay(1, decay, steps_per_epoch, decay_final, start)
    )
end

subtract_gradients(∂1::NamedTuple, ∂2::NamedTuple) = map(subtract_gradients, ∂1, ∂2)
subtract_gradients(∂1::AbstractArray, ∂2::AbstractArray) = ∂1 - ∂2

add_gradients(∂1::NamedTuple, ∂2::NamedTuple) = map(add_gradients, ∂1, ∂2)
add_gradients(∂1::AbstractArray, ∂2::AbstractArray) = ∂1 + ∂2
add_gradients(∂1::Union{NamedTuple,AbstractArray}, ::Nothing) = ∂1
add_gradients(::Nothing, ∂2::Union{NamedTuple,AbstractArray}) = ∂2

accum_gradients!(∂::NamedTuple, ∂other::NamedTuple) = map(accum_gradients!, ∂, ∂other)
accum_gradients!(∂::AbstractArray, ∂other::AbstractArray) = ∂ .+= ∂other
accum_gradients!(∂::Union{AbstractArray,Nothing}, ::Nothing) = ∂

function combine_gradients(op, ∂1::NamedTuple, ∂2::NamedTuple)
    _op(∂1::NamedTuple, ∂2::NamedTuple) = map(_op, ∂1, ∂2)
    _op(∂1::AbstractArray, ∂2::AbstractArray) = op(∂1, ∂2)
    return map(_op, ∂1, ∂2)
end

"""
    gradmult(∂, λ)

Multiplies gradients by a scalar `λ`.
"""
gradmult(∂::AbstractArray, λ::Real) = oftype(∂, λ * ∂)
gradmult(∂::Nothing, ::Real) = ∂
gradmult(∂::NamedTuple, λ::Real) = map(∂) do x
    gradmult(x, λ)
end


text:
module LibSnappy

using snappy_jll
export snappy_jll

using CEnum

@cenum snappy_status::UInt32 begin
    SNAPPY_OK = 0
    SNAPPY_INVALID_INPUT = 1
    SNAPPY_BUFFER_TOO_SMALL = 2
end

function snappy_compress(input, input_length, compressed, compressed_length)
    ccall((:snappy_compress, libsnappy), snappy_status, (Ptr{Cchar}, Csize_t, Ptr{Cchar}, Ptr{Csize_t}), input, input_length, compressed, compressed_length)
end

function snappy_uncompress(compressed, compressed_length, uncompressed, uncompressed_length)
    ccall((:snappy_uncompress, libsnappy), snappy_status, (Ptr{Cchar}, Csize_t, Ptr{Cchar}, Ptr{Csize_t}), compressed, compressed_length, uncompressed, uncompressed_length)
end

function snappy_max_compressed_length(source_length)
    ccall((:snappy_max_compressed_length, libsnappy), Csize_t, (Csize_t,), source_length)
end

function snappy_uncompressed_length(compressed, compressed_length, result)
    ccall((:snappy_uncompressed_length, libsnappy), snappy_status, (Ptr{Cchar}, Csize_t, Ptr{Csize_t}), compressed, compressed_length, result)
end

function snappy_validate_compressed_buffer(compressed, compressed_length)
    ccall((:snappy_validate_compressed_buffer, libsnappy), snappy_status, (Ptr{Cchar}, Csize_t), compressed, compressed_length)
end

const SNAPPY_MAJOR = 1

const SNAPPY_MINOR = 1

const SNAPPY_PATCHLEVEL = 9

const SNAPPY_VERSION = (SNAPPY_MAJOR << 16 | SNAPPY_MINOR << 8) | SNAPPY_PATCHLEVEL

end # module


text:
module HybridizedSBPUtils

using LinearAlgebra
using SparseArrays
using UnPack

using StartUpDG

# constructing reference hybridized SBP operators
export build_hSBP_ops

# for constructing global DG matrices
export build_rhs_matrix
export global_SBP_matrices_2D

"function build_hSBP_ops(rd::RefElemData)
    Build reference hybridized SBP operators"
function build_hSBP_ops(rd::RefElemData)
        @unpack M,Dr,Ds,Vq,Pq,Vf,wf,nrJ,nsJ = rd

        Qr = Pq'*M*Dr*Pq
        Qs = Pq'*M*Ds*Pq
        Ef = Vf*Pq
        Br = diagm(wf.*nrJ)
        Bs = diagm(wf.*nsJ)
        Qrh = .5*[Qr-Qr' Ef'*Br;
                -Br*Ef  Br]
        Qsh = .5*[Qs-Qs' Ef'*Bs;
                -Bs*Ef  Bs]

        Vh = [Vq;Vf]
        Ph = M\transpose(Vh)
        VhP = Vh*Pq

        # make skew symmetric versions of the operators"
        Qrhskew = .5*(Qrh-transpose(Qrh))
        Qshskew = .5*(Qsh-transpose(Qsh))

        return Qrhskew,Qshskew,Vh,Ph,VhP
end

"""
function build_rhs_matrix(applyRHS,Np,K,vargs...)

flexible but slow function to construct global matrices based on rhs evals
applyRHS = function to evaluate rhs f(u(t)) given a solution vector u(t)
Np,K = number dofs and elements
vargs = other args for applyRHS(...)

# Examples
```jldoctest
"""
function build_rhs_matrix(applyRHS,Np,K,vargs...)
    u = zeros(Np,K)
    A = spzeros(Np*K,Np*K)
    for i in eachindex(u)
        u[i] = one(eltype(u))
        r_i = applyRHS(u,vargs...)
        A[:,i] = droptol!(sparse(r_i[:]),1e-12)
        u[i] = zero(eltype(u))
    end
    return A
end

"""
function global_SBP_matrices_2D(rd::RefElemData, md::MeshData,
                                         Qrhskew, Qshskew, dtol=1e-12)

rd,md = ref elem and mesh data
Ax,Ay,Bx,By = sparse global operators for vol, surface terms in DG.
B = global operator for surface term with only sJ

note that md::MeshData needs FToF to also be periodic

# Examples
```jldoctest
"""

function global_SBP_matrices_2D(rd::RefElemData, md::MeshData,
                                 Qrhskew, Qshskew, dtol=1e-12)

    @unpack rxJ,sxJ,ryJ,syJ,nxJ,nyJ,sJ,mapP,FToF = md
    @unpack Nfaces,wf,wq = rd
    Nh = size(Qrhskew,1)
    NfqNfaces,K = size(nxJ)
    Nfq = convert(Int,NfqNfaces/Nfaces)

    @assert Nh==size(rxJ,1) "size of geofacs = $(size(rxJ)), while SBP ops are size $Nh."

    EToE = @. (FToF.-1) ÷ Nfaces + 1
    mapPerm = ((mapP.-1) .% NfqNfaces) .+ 1 # mod out face offset
    fids = length(wq)+1:Nh # last indices correspond to face nodes

    ids(e) = (1:Nh) .+ (e-1)*Nh # block offsets
    Block(e1,e2) = CartesianIndices((ids(e1),ids(e2))) # emulating BlockArrays, but faster
    face_ids = (x,f)->reshape(x,Nfq,Nfaces)[:,f]

    Ax,Ay,Bx,By,B = ntuple(x->spzeros(Nh*K,Nh*K),5)
    for e = 1:K # loop over elements

        # self-contributions - assume affine for now
        Ax_local = rxJ[1,e]*Qrhskew + sxJ[1,e]*Qshskew
        Ay_local = ryJ[1,e]*Qrhskew + syJ[1,e]*Qshskew
        Ax[Block(e,e)] .= sparse(Ax_local)
        Ay[Block(e,e)] .= sparse(Ay_local)

        # neighboring contributions (ignore self-neighbors)
        for (f,enbr) in enumerate(EToE[:,e])
            if enbr!=e
                fperm = face_ids(reshape(mapPerm[:,e],Nfq,Nfaces),f)
                Axnbr = spdiagm(0 => face_ids(.5*wf.*nxJ[:,e],f))
                Aynbr = spdiagm(0 => face_ids(.5*wf.*nyJ[:,e],f))
                Bx[Block(e,enbr)[face_ids(fids,f),fids[fperm]]] .= Axnbr
                By[Block(e,enbr)[face_ids(fids,f),fids[fperm]]] .= Aynbr
                B[ Block(e,enbr)[face_ids(fids,f),fids[fperm]]] .= spdiagm(0 => face_ids(.5*wf.*sJ[:,e],f))
            end
        end
    end
    return droptol!.((Ax+Bx,Ay+By,Bx,By,B),dtol)
end

end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE

"""
    calculate_normals(elements::Dict{Int, Vector{Int}},
                      element_types::Dict{Int, Symbol},
                      X::Dict{Int, Vector{Float64})

Given elements, element types and node locations, calculate nodal normals by
first calculating normal directions for each element and then averaging them
in nodes. As a result we get unique normal direction defined to each node.

# Notes
Only linear elements supported.

# Example

```jldoctest
X = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 0.0, 0.0], 3 => [0.0, 1.0, 0.0])
elements = Dict(1 => [1, 2, 3])
element_types = Dict(1 => :Tri3)
normals = calculate_normals(elements, element_types, X)

# output

Dict{Int64,Array{Float64,1}} with 3 entries:
  2 => [0.0, 0.0, 1.0]
  3 => [0.0, 0.0, 1.0]
  1 => [0.0, 0.0, 1.0]

```
"""
function calculate_normals(elements, element_types, X)
    normals = similar(X)
    for (elid, elcon) in elements
        @assert element_types[elid] == :Tri3
        p1 = X[elcon[1]]
        p2 = X[elcon[2]]
        p3 = X[elcon[3]]
        n = cross(p2-p1, p3-p1)
        n /= norm(n)
        for nid in elcon
            haskey(normals, nid) || (normals[nid] = zeros(3))
            normals[nid] += n
        end
    end
    for (nid, n) in normals
        normals[nid] /= norm(normals[nid])
    end
    return normals
end


text:
#=============================#
#== put needed modules here ==#
#=============================#
import JuliaChem

#================================#
#== JuliaChem execution script ==#
#================================#
function full_rhf(input_file)
  println("--------------------------------------------------------------------------------")
  println("                       ========================================                 ")
  println("                               WELCOME TO JULIACHEM.JL!                         ")
  println("                       ========================================                 ")
  println("                                                                                ")
  println("              JuliaChem.jl is an electronic structure theory package            ")
  println("         developed by the Gordon group at Iowa State University. It is          ")
  println("        designed to apply the strengths of the Julia programming language       ")
  println("          (high-level semantics and low-level performance) to a quantum         ")
  println("                                chemistry package.                              ")
  println("                                                                                ")
  println("         For more information, refer to the following paper: Poole, D.          ")
  println("        Galvaz Vallejo, J. L.; Gordon, M. S. \"A New Kid on the Block:          ")
  println("         Application of Julia to Hartree-Fock Calculations.\" J. Chem.          ")
  println("                     Theory Compute. 2020, 16, 8, 5006-5013.                    ")
  println("                                                                                ")
  println("       For questions on usage, email David Poole at davpoole@iastate.edu.       ")
  println("--------------------------------------------------------------------------------")

  try
    #== initialize JuliaChem ==#
    JuliaChem.initialize()

    #== read in input file ==#
    molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
      output=2)

    #== generate basis set ==#
    mol, basis = JuliaChem.JCBasis.run(molecule, model;
      output=2)

    #== molecule info ==#
    JuliaChem.JCMolecule.run(mol)

    #== calculation driver ==#
    if driver == "energy"
      if model["method"] == "RHF"
        #== perform scf calculation ==#
        if haskey(keywords, "scf")
          rhf_energy = JuliaChem.JCRHF.Energy.run(mol, basis, keywords["scf"];
            output=2)
        else
          rhf_energy = JuliaChem.JCRHF.Energy.run(mol, basis;
            output=2)
        end

        #== compute molecular properties such as dipole moment ==#
        properties = JuliaChem.JCRHF.Properties.run(mol, basis, rhf_energy,
          keywords["prop"]; output=2)
      else
        throw("Exception: Methods other than RHF are not supported yet!")
      end
    else
      throw("Exception: Only energy calculations are currently supported!")
    end
  catch e
    bt = catch_backtrace()
    msg = sprint(showerror, e, bt)
    println(msg)

    JuliaChem.finalize()
    exit()
  end

  println("--------------------------------------------------------------------------------")
  println("                      Your calculation has run to completion!                   ")
  println("                                                                                ")
  println("                       ========================================                 ")
  println("                                   HAVE A NICE DAY!                             ")
  println("                       ========================================                 ")
  println("--------------------------------------------------------------------------------")
end

full_rhf(ARGS[1])


text:
using MixedModels,Base.Test
const m1 = LowerTriangular(eye(3))

@test size(m1) == (3,3)
@test size(m1,1) == 3
@test size(m1,2) == 3
@test size(m1,3) == 1
const m1c = copy(m1)
@test MixedModels.nlower(m1) == 6
@test m1c == m1
@test getθ(m1) == [1., 0, 0, 1, 0, 1]
@test lowerbd(m1) == [0., -Inf, -Inf, 0., -Inf, 0.]

setθ!(m1, collect(1.:6))                  # assignment of parameters
@test getθ(m1) == collect(1.:6)           # check assignment and extraction
@test m1c != m1
@test copy!(m1c,m1) == m1

@test full(m1) == reshape([1.,2,3,0,4,5,0,0,6],(3,3))


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
    _add_net_injection_eqs!(model)
    _add_reserve_eqs!(model)
    return
end

function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
    T = model[:instance].time
    net_injection = _init(model, :net_injection)
    eq_net_injection = _init(model, :eq_net_injection)
    eq_power_balance = _init(model, :eq_power_balance)
    for t in 1:T, b in model[:instance].buses
        n = net_injection[b.name, t] = @variable(model)
        eq_net_injection[b.name, t] =
            @constraint(model, -n + model[:expr_net_injection][b.name, t] == 0)
    end
    for t in 1:T
        eq_power_balance[t] = @constraint(
            model,
            sum(net_injection[b.name, t] for b in model[:instance].buses) == 0
        )
    end
    return
end

function _add_reserve_eqs!(model::JuMP.Model)::Nothing
    eq_min_reserve = _init(model, :eq_min_reserve)
    instance = model[:instance]
    for t in 1:instance.time
        # Equation (68) in Kneuven et al. (2020)
        # As in Morales-España et al. (2013a)
        # Akin to the alternative formulation with max_power_avail
        # from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
        shortfall_penalty = instance.shortfall_penalty[t]
        eq_min_reserve[t] = @constraint(
            model,
            sum(model[:reserve][g.name, t] for g in instance.units) +
            (shortfall_penalty >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
            instance.reserves.spinning[t]
        )

        # Account for shortfall contribution to objective
        if shortfall_penalty >= 0
            add_to_expression!(
                model[:obj],
                shortfall_penalty,
                model[:reserve_shortfall][t],
            )
        end
    end
    return
end


text:
using PiCraft, Test

@testset "NBT/schematic File reading" begin
    @test parseNBT("../example/hello_world.nbt", devnull) == 0
    @test parseNBT("../example/bigtest.nbt", devnull) == 0
    m = Array{Block, 3}(undef, 3, 2, 4)
    m[:, :, 1] = Block.([41 0; 41 0; 41 0])
    m[:, :, 4] = m[:, :, 3] = m[:, :, 2] = m[:, :, 1]
    m[3, 2, 3] = m[3, 2, 4] = Block(57)
    @test importSchematic("../example/sample.schematic") == m
end


text:
using MSM
using Test

@testset "MSM.jl" begin
    # Write your tests here.
end


text:
@testset "DataMatrix" begin
    a = @data [1.0, 2.0, 3.0]
    v_a = [1.0, 2.0, 3.0]

    b = @data eye(3, 3)
    m_b = eye(3, 3)

    #
    # Transposes
    #

    @test all(a' .== v_a')
    @test all(a'' .== v_a'')
    @test all(b' .== m_b')
    @test all(b'' .== m_b'')

    #
    # DataVector * DataMatrix
    #

    # TODO: Get indexing for b[1, :] to work
    # @test all(a * b[1, :] .== v_a * m_b[1, :])

    #
    # DataMatrix * DataVector
    #

    @test all(b * a .== m_b * v_a)
    @test all(convert(Array, b * a) .== m_b * v_a)

    #
    # DataMatrix * DataMatrix
    #

    @test all(b * b .== m_b * m_b)

    #
    # DataVector * DataMatrix w/ missings
    #

    b[1, 1] = missing
    res = a * b[1:1, :]
    @test all(ismissing.(res[:, 1]))
    @test all(.!(ismissing.(res[:, 2])))
    @test all(.!(ismissing.(res[:, 3])))
    res = a * b[2:2, :]
    @test all(.!(ismissing.(res)))

    #
    # DataMatrix w missings * DataVector
    #

    res = b * a
    @test ismissing.(res[1])
    @test .!(ismissing.(res[2]))
    @test .!(ismissing.(res[3]))

    #
    # DataMatrix * DataMatrix
    #

    res = b * b
    # 3x3 Float64 DataMatrix:
    #  missing   missing   missing
    #  missing  1.0  0.0
    #  missing  0.0  1.0
    @test ismissing.(res[1, 1])
    @test ismissing.(res[1, 2])
    @test ismissing.(res[1, 3])
    @test ismissing.(res[2, 1])
    @test .!(ismissing.(res[2, 2]))
    @test .!(ismissing.(res[2, 3]))
    @test ismissing.(res[3, 1])
    @test .!(ismissing.(res[3, 2]))
    @test .!(ismissing.(res[3, 3]))

    res = b * @data eye(3)
    # 3x3 Float64 DataMatrix:
    #   missing   missing   missing
    #  0.0  1.0  0.0
    #  0.0  0.0  1.0
    @test ismissing.(res[1, 1])
    @test ismissing.(res[1, 2])
    @test ismissing.(res[1, 3])
    @test .!(ismissing.(res[2, 1]))
    @test .!(ismissing.(res[2, 2]))
    @test .!(ismissing.(res[2, 3]))
    @test .!(ismissing.(res[3, 1]))
    @test .!(ismissing.(res[3, 2]))
    @test .!(ismissing.(res[3, 3]))

    res = (@data eye(3)) * b
    # julia> dataeye(3) * b
    # 3x3 Float64 DataMatrix:
    #  missing  0.0  0.0
    #  missing  1.0  0.0
    #  missing  0.0  1.0
    @test ismissing.(res[1, 1])
    @test .!(ismissing.(res[1, 2]))
    @test .!(ismissing.(res[1, 3]))
    @test ismissing.(res[2, 1])
    @test .!(ismissing.(res[2, 2]))
    @test .!(ismissing.(res[2, 3]))
    @test ismissing.(res[3, 1])
    @test .!(ismissing.(res[3, 2]))
    @test .!(ismissing.(res[3, 3]))

    # Test row operations
    dm = @data eye(6, 2)
    # rowmeans(dm)

    # Test column operations
    dm = @data eye(6, 2)
    # colmeans(dm)

    # Test linear algebra
    du, dd, dv = svd((@data eye(3, 3)))
    u, d, v = svd(eye(3, 3))
    @test all(du .== u)
    @test all(dd .== d)
    @test all(dv .== v)

    # Test elementary functions
    dm = -(@data eye(5, 5))
    @test all(abs(dm) .== eye(5, 5))
end


text:
#=
BFSParents:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

#function to insert inside the result array
function insert(p, result)
    for i=1:length(p)
        if(p[i]!=nothing && result[i]==nothing)
            result[i] = p[i]
        end
    end

    return result

end

#function to update visited nodes
function updateVisited(f, visited)
    for i=1:length(f)
        if(f[i]!=nothing)
            visited[i] = true
        end
    end

    return visited

end

#A is the input matrix, s is the source node, and n is the number of nodes, Output: parent matrix
function bfs_parents(A, s, n)

    index = GBVector{Int64}(n)
        for i = 1:n
            index[i] = i
        end

    #wavefront vector -> w[s]=source node insrted in wavefront
    f = GBVector{Int64}(n)
    f[s] = 1

    #parent vector -> p[s]=1 because the source is already visited
    p = GBVector{Int64}(n)

    #result vector -> result[s]=0 because the parent of the source is considered node 0
    result = GBVector{Int64}(n)
    result[s] = 0

    #visited vector -> v[s]=1 because the source is already visited
    visited = GBVector{Bool}(n)
    visited[s] = true

    print("STARTING ...\n\n")
    for i = 1:n-1
            f = mul(f, A, Semirings.MIN_FIRST, mask=p, desc=Descriptors.SC)
            p = f[:, mask=f, desc=Descriptors.S]
            f = index[:, mask=f, desc=Descriptors.S]
            result = insert(p, result)
            visited = updateVisited(f, visited)
            unvisited = filter(x -> x!=true, visited)
            if(isempty(unvisited))
                break
            end
            print("\nITERATION FINISHED\n\n")
            print(result)

    end

    print("COMPLETED\n\n")
    print("RESULT\n\n")
    print(result')

    return result

end


text:
# Make sure the examples in the README continue to function

@testset "$(rpad("examples",80))" begin
    using DomainSets, StaticArrays
    @test repr(UnitInterval()) == "0.0..1.0 (Unit)"
    @test repr(ChebyshevInterval()) == "-1.0..1.0 (Chebyshev)"
    @test repr(HalfLine()) == "0.0..Inf (closed–open) (HalfLine)"

    using DomainSets: ×
    @test repr((-1..1) × (0..3) × (4.0..5.0)) == "(-1.0..1.0) × (0.0..3.0) × (4.0..5.0)"
    @test SVector(1,2) in (-1..1) × (0..3)

    @test SVector(0,0,1.0) in UnitSphere(Val(3))
    @test [0.0,1.0,0.0,0.0] in UnitSphere(4)
    @test SVector(1,0) in UnitCircle()

    @test SVector(0.1,0.2,0.3) in UnitBall(Val(3))
    @test [0.1,0.2,0.3,-0.1] in UnitBall(4)
    @test SVector(0.1,0.2) in UnitDisk()

    @test 1:5 in ProductDomain([0..i for i in 1:5])
    @test ("a", 0.4) ∈ ProductDomain(["a","b"], 0..1)

    d = UnitCircle() ∪ 2UnitCircle()
    @test in.([SVector(1,0),SVector(0,2), SVector(1.5,1.5)], Ref(d)) == [1,1,0]
    d = UnitCircle() ∩ (2UnitCircle() .+ SVector(1.0,0.0))
    @test !(SVector(1,0) in d)
    @test SVector(-1,0) in d

    d = LevelSet{SVector{2,Float64}}(prod, 1.0)
    @test [0.5,2] ∈ d

    d = IndicatorFunction{Float64}( t ->  cos(t) > 0)
    @test (0.5 ∈ d, 3.1 ∈ d) == (true, false)

    d = Domain(x>0 for x in -1..1)
    @test (0.5 ∈ d, -0.5 ∈ d) == (true, false)

    d = Domain( x*y > 0 for (x,y) in UnitDisk())
    @test ([0.2, 0.3] ∈ d, [0.2, -0.3] ∈ d) == (true, false)

    d = Domain( x+y+z > 0 for (x,y,z) in ProductDomain(UnitDisk(), 0..1))
    @test [0.3,0.2,0.5] ∈ d
end


text:
module DiffOperatorsTests

using Gridap
using Test
using Gridap.ODEs
using Gridap.ODEs.ODETools: ∂t, ∂tt

using ForwardDiff

f(x,t) = 5*x[1]*x[2]+x[2]^2*t^3
∂tf(x,t) = x[2]^2*3*t^2

tv = rand(Float64)
xv = Point(rand(Float64,2)...)
@test ∂tf(xv,tv) ≈ ∂t(f)(xv,tv)

F(x,t) = VectorValue([5*x[1]*x[2],x[2]^2*t^3])
∂tF(x,t) = VectorValue([0.0,x[2]^2*3*t^2])

tv = rand(Float64)
xv = Point(rand(Float64,2)...)
@test ∂tF(xv,tv) ≈ ∂t(F)(xv,tv)

# Time derivatives

f(x,t) = t^2
dtf = (x,t) -> ForwardDiff.derivative(t->f(x,t),t)
@test dtf(xv,tv) ≈ ∂t(f)(xv,tv) ≈ ∂t(f)(xv)(tv) ≈ ∂t(f)(tv)(xv)

f2(x,t) = x[1]^2
dtf2 = (x,t) -> ForwardDiff.derivative(t->f2(x,t),t)
@test dtf2(xv,tv) ≈ ∂t(f2)(xv,tv) ≈ ∂t(f2)(xv)(tv) ≈ ∂t(f2)(tv)(xv)

f2(x,t) = x[1]^t^2
dtf2 = (x,t) -> ForwardDiff.derivative(t->f2(x,t),t)
@test dtf2(xv,tv) ≈ ∂t(f2)(xv,tv) ≈ ∂t(f2)(xv)(tv) ≈ ∂t(f2)(tv)(xv)

f2(x,t) = VectorValue(x[1]^2,0.0)
dtf2 = (x,t) -> VectorValue(ForwardDiff.derivative(t -> get_array(f2(x,t)),t))
@test dtf2(xv,tv) ≈ ∂t(f2)(xv,tv) ≈ ∂t(f2)(xv)(tv) ≈ ∂t(f2)(tv)(xv)

f2(x,t) = VectorValue(x[1]^2,t)
dtf2 = (x,t) -> VectorValue(ForwardDiff.derivative(t -> get_array(f2(x,t)),t))
@test dtf2(xv,tv) ≈ ∂t(f2)(xv,tv) ≈ ∂t(f2)(xv)(tv) ≈ ∂t(f2)(tv)(xv)

f6(x,t) = TensorValue(x[1]*t,x[2]*t,x[1]*x[2],x[1]*t^2)
dtf6 = (x,t) -> TensorValue(ForwardDiff.derivative(t->get_array(f6(x,t)),t))
@test dtf6(xv,tv) ≈ ∂t(f6)(xv,tv) ≈ ∂t(f6)(xv)(tv) ≈ ∂t(f6)(tv)(xv)

# Spatial derivatives
f2(x,t) = VectorValue(x[1]^2,t)
∇f2(x,t) = ∇(y->f2(y,t))(x)
∇f2(xv,tv)
# @santiagobadia : Is there any way to make this transparent to the user
# I guess not unless we create a type for these analytical (space-only or
# space-time via a trait) functions
# Probably a try-catch?

# 2nd time derivative
f(x,t) = t^2
dtf = (x,t) -> ForwardDiff.derivative(t->f(x,t),t)
dttf = (x,t) -> ForwardDiff.derivative(t->dtf(x,t),t)
@test dttf(xv,tv) ≈ ∂tt(f)(xv,tv) ≈ ∂tt(f)(xv)(tv) ≈ ∂tt(f)(tv)(xv)
@test ∂tt(f)(xv,tv) ≈ 2.0

f2(x,t) = x[1]*t^2
dtf2 = (x,t) -> ForwardDiff.derivative(t->f2(x,t),t)
dttf2 = (x,t) -> ForwardDiff.derivative(t->dtf2(x,t),t)
@test dttf2(xv,tv) ≈ ∂tt(f2)(xv,tv) ≈ ∂tt(f2)(xv)(tv) ≈ ∂tt(f2)(tv)(xv)
@test ∂tt(f2)(xv,tv) ≈ 2.0*xv[1]

end #module


text:
module Rebugger

using UUIDs
using REPL
import REPL.LineEdit, REPL.Terminals
using REPL.LineEdit: buffer, bufend, content, edit_splice!
using REPL.LineEdit: transition, terminal, mode, state

using Revise
using Revise: ExLike, RelocatableExpr, get_signature, funcdef_body, get_def, striplines!
using Revise: printf_maxsize
using HeaderREPLs

const msgs = []  # for debugging. The REPL-magic can sometimes overprint error messages

include("debug.jl")
include("ui.jl")
include("deepcopy.jl")

# Set up keys that enter rebug mode from the regular Julia REPL
# This should be called from your ~/.julia/config/startup.jl file
function repl_init(repl)
    repl.interface = REPL.setup_interface(repl; extra_repl_keymap = rebugger_modeswitch)
end

function __init__()
    # Set up the Rebugger REPL mode with all of its key bindings
    repl_inited = isdefined(Base, :active_repl)
    @async begin
        while !isdefined(Base, :active_repl)
            sleep(0.05)
        end
        sleep(0.1) # for extra safety
        # Set up the custom "rebug" REPL
        main_repl = Base.active_repl
        repl = HeaderREPL(main_repl, RebugHeader())
        interface = REPL.setup_interface(repl; extra_repl_keymap=[rebugger_modeswitch, rebugger_keys])
        rebug_prompt_ref[] = interface.modes[end]
        # Add F5 to the history prompt
        history_prompt = find_prompt(main_repl.interface, LineEdit.PrefixHistoryPrompt)
        add_key_stacktrace!(history_prompt.keymap_dict)
        # If the REPL was already initialized, add the keys to the julia> prompt now
        # (This will already be done if the user turned on the key bindings in her startup.jl file)
        if repl_inited
            julia_prompt = find_prompt(main_repl.interface, "julia")
            add_key_stacktrace!(julia_prompt.keymap_dict)
            add_key_stepin!(julia_prompt.keymap_dict)
        end
    end
end

end # module


text:
module HuntTest

using Gridap
using Test

using GridapMHD: hunt
using GridapMHD: analytical_hunt_u
using GridapMHD: analytical_hunt_j
using GridapMHD: compute_u_j_errors

# Problem setting
ρ = 1.0
ν = 1.0
σ = 1.0

U0 = 10.0
B0 = 10.0
L = 1.0

Re = U0 * L / ν
Ha = B0 * L * sqrt(σ/(ρ*ν))
K = Ha / (1-0.95598*Ha^(-1/2)-Ha^(-1))
∂p∂z = -Re * K / L^3

f_u = VectorValue(0.0,0.0, -∂p∂z) * L/U0^2

xh, trian, dΩ = hunt(nx=5,ny=5,Re=Re,Ha=Ha,f_u=f_u)

# Postprocess
uh, ph, jh, φh = xh
divj = (∇⋅jh)

uh = uh * U0
jh = jh * σ * B0 * U0

# Analyical solutions
side_wall_semilength = 1.0
hartmann_wall_semilength = 1.0
number_fourier_sumands = 10
K = Ha / (1-0.95598*Ha^(-1/2)-Ha^(-1))
∂p∂z = -Re * K / L^3
u0(x) = analytical_hunt_u(side_wall_semilength, hartmann_wall_semilength,
  ρ*ν, ∂p∂z, Ha, number_fourier_sumands, x)

j0(x) = analytical_hunt_j(side_wall_semilength, hartmann_wall_semilength,
  σ, ρ*ν, ∂p∂z, Ha, number_fourier_sumands, x)

eu_l2, ej_l2 = compute_u_j_errors(uh, jh, u0, j0, dΩ)
e_divj = sqrt(sum(∫(divj*divj)* dΩ))

@test eu_l2 < 0.06
@test ej_l2 < 0.8
@test e_divj < 1e-9
end #module


text:
module BijectiveHilbert

include("bitops.jl")
include("gray_code.jl")
include("hilbert_algorithm.jl")
export index_type
export encode_hilbert
export decode_hilbert!

include("global_gray.jl")
export GlobalGray
export axis_type
export encode_hilbert_zero
export decode_hilbert_zero!

include("hamilton.jl")
export SpaceGray
export Compact

# include("bijective.jl")
include("simple2d.jl")
import .Bijective: Simple2D

@doc """
Simple2D(::Type{T})

The type is a data type to hold the Hilbert index. It should have
enough bits to hold all bits of integer axes to encode.

The variant algorithm used differs from the usual Hilbert code because it
doesn't need to know the size of the whole grid before computing the code.
It looks like a slightly-rotated version of the Hilbert curve, but it
has the benefit that it is 1-1 between `(x, y)` and `z`, so you can translate
back and forth.

If the size of the axis values or the size of the Hilbert index are too large
to be stored in the data types of the axis vector or index, then you'll
see an error from a failure of `trunc`, which tries to copy values into the
smaller datatype. Solve this by using a larger datatype, so UInt64 instead
of UInt32.

It comes from a paper:
N. Chen, N. Wang, B. Shi, A new algorithm for encoding and decoding
the Hilbert order. Software—Practice and Experience 2007; 37(8): 897–908.
"""
Simple2D

export Simple2D

include("facecontinuous.jl")
export FaceContinuous
end


text:
#=##############################################################################
# DESCRIPTION
    Methods for generation and evaluation of multigrids.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : May 2018
  * License   : MIT License
=###############################################################################

################################################################################
# MULTIGRID
################################################################################

mutable struct MultiGrid <: AbstractGrid

  # User inputs
  dims::Int64                         # Number of dimensions

  # Properties
  grids::Array{T,1} where {T<:AbstractGrid}     # Grids
  grid_names::Array{String, 1}        # Grid names
  nnodes::Int64                       # Number of nodes
  ncells::Int64                       # Number of cells
  ngrids::Int64                       # Number of grids
  # bbox::Array{Int64, 1}               # Bounding box (cells in each dimension)

  # Internal data
  _grid_origins::Array{Array{Int64,1},1}  # Origin of each grid

  MultiGrid(          dims,
                      grids=AbstractGrid[],
                        grid_names=String[],
                        nnodes=0,
                        ncells=0,
                        ngrids=0,
                        # bbox=zeros(Int64, dims),
                      _grid_origins=Array{Int64,1}[]
            ) = new(
                      dims,
                      grids,
                        grid_names,
                        nnodes,
                        ncells,
                        ngrids,
                        # bbox,
                      _grid_origins
            )
end

function addgrid(self::MultiGrid, grid_name::String, grid::T) where {T<:AbstractGrid}

  if grid_name in self.grid_names
    error("Grid $grid_name already exists!")
  end

  self.nnodes += grid.nnodes
  self.ncells += grid.ncells
  self.ngrids += 1

  push!(self.grids, grid)
  push!(self.grid_names, grid_name)
end

"Returns the requested grid"
function get_grid(self::MultiGrid, grid_name::String)

  grids_found = find(x->x==grid_name, self.grid_names)

  if length(grids_found)==0
    error("$grid_name not found. Grids available are $(self.grid_names).")
  elseif length(grids_found)!=1
    error("Logic error!")
  end

  return get_grid(self, grids_found[1])
end

"Returns the requested grid"
function get_grid(self::MultiGrid, grid_index::Int64)
  if grid_index>self.ngrids
    error("Requested invalid grid index $grid_index; max is $(self.ngrids).")
  end
  return self.grids[grid_index]
end

function get_node(self::MultiGrid, i::Int64)
  if i>self.nnodes
    error("Requested invalid node index $i; max is $(self.nnodes).")
  end

  ind = i
  for grid in self.grids
    if ind<=grid.nnodes
      return get_node(grid, ind)
    end
    ind -= grid.nnodes
  end

  error("Logic error!")
end
function get_node(self::MultiGrid, coor::Array{Int64,1})
  error("Not Implemented yet!")
end

function get_cell(self::MultiGrid, i::Int64)
  if i>self.ncells
    error("Requested invalid cell index $i; max is $(self.ncells).")
  end

  ind = i
  for grid in self.grids
    if ind<=grid.ncells
      return get_ncell(grid, ind)
    end
    ind -= grid.ncells
  end

  error("Logic error!")

end
function get_cell(self::MultiGrid, coor::Array{Int64,1})
  error("Not Implemented yet!")
end

function get_fieldval(self::MultiGrid, field_name::String, i::Int64)
  if i>self.nnodes
    error("Requested invalid node index $i; max is $(self.nnodes).")
  end

  ind = i
  for grid in self.grids
    if ind<=grid.nnodes
      return get_fieldval(grid, field_name, ind)
    end
    ind -= grid.nnodes
  end

  error("Logic error!")
end
function get_fieldval(self::MultiGrid, field_name::String, coor::Array{Int64,1})
  error("Not Implemented yet!")
end

function add_field(self::MultiGrid, field_name::String, field_type::String,
                                                                        field_data)
  ind = 0
  for grid in self.grids
    add_field(grid, field_name, field_type, field_data[ind+1:ind+grid.nnodes])
    ind += grid.nnodes
  end
end
function calculate_field(self::MultiGrid, args...)
  for grid in self.grids
    calculate_fields(grid, args...)
  end
end
function lintransform!(self::MultiGrid, args...; optargs...)
  for grid in self.grids
    lintransform!(grid, args...; optargs...)
  end
end
function transform!(self::MultiGrid, args...; optargs...)
end
function save(self::MultiGrid, file_name::String, args...; optargs...)
  str = ""
  for (i, grid) in enumerate(self.grids)
    str *= save(grid, file_name*"_"*self.grid_names[i], args...; optargs...)
  end
  return str
end
function plot(self::MultiGrid; optargs...)
  for grid in self.grids
    save(grid, args...; optargs...)
  end
end
##### INTERNAL FUNCTIONS  ######################################################

##### END OF MULTIGRID #########################################################


text:
function FGrad3Vec!(F,cCG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;

D1cCG = Param.CacheC1
D2cCG = Param.CacheC2
D3cCG = Param.CacheC3
D3cCGE =Param.CacheC4

mul!(reshape(D1cCG,OP,OP*NF*nz),CG.DS,reshape(cCG,OP,OP*nz*NF))
mul!(reshape(PermutedDimsArray(D2cCG,(2,1,3,4)),OP,OP*NF*nz),CG.DS,reshape(PermutedDimsArray(cCG,(2,1,3,4)),OP,OP*nz*NF))
@views D3cCG[:,:,:,1:nz-1] .= 0.5*(cCG[:,:,:,2:nz] .- cCG[:,:,:,1:nz-1])

@views F[:,:,:,1:nz-1,Param.wPos] .= F[:,:,:,1:nz-1,Param.wPos] .- Param.dXdxIF33.*D3cCG[:,:,:,1:nz-1];

if nz>1
@views  D3cCGE[:,:,:,1] .= D3cCG[:,:,:,1];
@views  D3cCGE[:,:,:,2:end-1] .= 0.5.*(D3cCG[:,:,:,1:end-2] .+ D3cCG[:,:,:,2:end-1]);
@views  D3cCGE[:,:,:,end] .= D3cCG[:,:,:,end-1];
else
    @views  D3cCGE[:,:,:,1] .= 0
end

@views F[:,:,:,:,Param.uPos] .= F[:,:,:,:,Param.uPos] .- (Param.dXdxIC11.*D1cCG .+
  Param.dXdxIC21.*D2cCG .+
  Param.dXdxIC31.*D3cCGE)

@views F[:,:,:,:,Param.vPos] .= F[:,:,:,:,Param.vPos] .- (Param.dXdxIC12.*D1cCG .+
  Param.dXdxIC22.*D2cCG .+ Param.dXdxIC32.*D3cCGE)

end

function FGrad3RhoVec!(F,cCG,RhoCG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;

D1cCG = Param.CacheC1
D2cCG = Param.CacheC2
D3cCG = Param.CacheC3
D3cCGE = Param.CacheC4

mul!(reshape(D1cCG,OP,OP*NF*nz),CG.DS,reshape(cCG,OP,OP*nz*NF))
mul!(reshape(PermutedDimsArray(D2cCG,(2,1,3,4)),OP,OP*NF*nz),CG.DS,reshape(PermutedDimsArray(cCG,(2,1,3,4)),OP,OP*nz*NF))
@views D3cCG[:,:,:,1:end-1] .= 0.5.*(cCG[:,:,:,2:end] .- cCG[:,:,:,1:end-1])

@views F[:,:,:,1:nz-1,Param.wPos] .= F[:,:,:,1:nz-1,Param.wPos] .-
    Param.dXdxIF33.*D3cCG[:,:,:,1:nz-1] ./ (0.5.*(RhoCG[:,:,:,1:nz-1]+RhoCG[:,:,:,2:nz]))
if nz>1
@views  D3cCGE[:,:,:,1] .= D3cCG[:,:,:,1];
@views  D3cCGE[:,:,:,2:end-1] .= 0.5.*(D3cCG[:,:,:,1:nz-2] .+ D3cCG[:,:,:,2:nz-1]);
@views  D3cCGE[:,:,:,end] .= D3cCG[:,:,:,nz-1];
else
    @views  D3cCGE[:,:,:,1] .= 0
end

@views F[:,:,:,:,Param.uPos] .= F[:,:,:,:,Param.uPos] .- (Param.dXdxIC11.*D1cCG .+
  Param.dXdxIC21.*D2cCG .+
  Param.dXdxIC31.*D3cCGE) ./ RhoCG

@views F[:,:,:,:,Param.vPos] .= F[:,:,:,:,Param.vPos] .- (Param.dXdxIC12.*D1cCG .+
  Param.dXdxIC22.*D2cCG .+
  Param.dXdxIC32.*D3cCGE) ./ RhoCG

end

function FGrad3Vec(cCG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
D1cCG=reshape(
  CG.DS*reshape(cCG,OP,OP*NF*nz)
  ,OP,OP,NF,nz);
D2cCG= permute(
  reshape(
  CG.DS *reshape(permute(cCG
  ,[2 1 3 4])
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz)
  ,[2 1 3 4]);

D3cCG=0.5*(cCG[:,:,:,2:nz]-cCG[:,:,:,1:nz-1]);

gradCG=zeros(OP,OP,NF,nz,3);
gradCG[:,:,:,1:nz-1,3]=Param.dXdxIF[:,:,:,2:nz,3,3].*D3cCG;
D3cCGE=zeros(OP,OP,NF,nz);
if nz>1
  D3cCGE[:,:,:,1]=D3cCG[:,:,:,1];
  D3cCGE[:,:,:,2:nz-1]=0.5*(D3cCG[:,:,:,1:end-1]+D3cCG[:,:,:,2:end]);
  D3cCGE[:,:,:,nz]=D3cCG[:,:,:,nz-1];
end

gradCG[:,:,:,:,1]=Param.dXdxIC[:,:,:,:,1,1].*D1cCG+
  Param.dXdxIC[:,:,:,:,2,1].*D2cCG+
  Param.dXdxIC[:,:,:,:,3,1].*D3cCGE;
gradCG[:,:,:,:,2]=Param.dXdxIC[:,:,:,:,1,2].*D1cCG+
  Param.dXdxIC[:,:,:,:,2,2].*D2cCG+
  Param.dXdxIC[:,:,:,:,3,2].*D3cCGE;

return gradCG
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

"""
```julia
struct FinitePolynomial <: FiniteField
```

Has parameters `p`, which is an abstract polynomial, and `n` which is the modulus of the field under which the molynomial is defined.

---

```julia
FinitePolynomial(p::AbstractPolynomial, n::Int)
```

A constructor method for `FinitePolynomial`.  Takes in a polynomial `p` and a number `n`, and constructs a polynomial under modulo n.
"""
struct FinitePolynomial <: FiniteField
	p::AbstractPolynomial
	n::Int

	function FinitePolynomial(p::AbstractPolynomial, n::Int)
		p = Polynomial(mod.(p.coeffs, n))
		new(p, n)
	end
end

"""
```julia
mod(p::Polynomial, n::Int) -> Polynomial
```

Uses the `FinitePolynomial` constructor to return a polynomial `p` under modulus `n`.

Parameters:
  - `p::Polynomial`: The input polynomial.
  - `n::Int`: The modulus of the field.

Returns
  - Polynomial: A polynomial modulo n.
"""
Base.mod(p::Polynomial, n::Int) = FinitePolynomial(p, n).p

"""
```julia
Polynomial(A::Union{NTuple{N, T}, Vector{T}}, n::Int) -> Polynomial
```

Constructs a polynomial under modulo `n`.

Parameters:
  - `A::Union{Tuple, AbstractArray}`: The polynomial coefficients.
  - `n::Int`: The modulus of the field.

Returns
  - `Polynomial`: A polynomial modulo n.
"""
Polynomial(A::Union{NTuple{N, T}, Vector{T}}, n::Int) where {N, T} = mod(Polynomial(A), n)

"""
```julia
list_polys(n::Int, m::Int) -> Array
```

Lists all polynomials of degree less than to `n` under modulo `m`.

Parameters:
  - `n::Int`: Highest degree of polynomial.
  - `m::Int`: The modulus of the field.

Returns:
  - `Array`: An array of polynomials of degree less than n, under modulo m.
"""
function list_polys(n::Int, m::Int)
	return collect(Polynomial(collect(t)) for t in Iterators.product([0:(m-1) for i in 1:n]...))
end

"""
```julia
multiplication_table(degree::Int, modulo::Int) -> Matrix
```

Returns a table (matrix) of the multiplication of all combinations of polynomials for degree less than `degree`, under modulo `modulo`.

Parameters:
  - `degree::Int`: Highest degree of polynomial.
  - `modulo::Int`: The modulus of the field.

Returns:
  - `Matrix`: A multiplication table of all polynomials with degree less than n, under modulus.

---

### Examples

```julia
julia> julia> multiplication_table(2, 3) # multiplication table of all polynomials of degree less than 3 modulo 2
9×9 Array{Polynomial,2}:
 Polynomial(0)  Polynomial(0)        Polynomial(0)        …  Polynomial(0)                Polynomial(0)
 Polynomial(0)  Polynomial(1)        Polynomial(2)           Polynomial(1 + 2*x)          Polynomial(2 + 2*x)
 Polynomial(0)  Polynomial(2)        Polynomial(1)           Polynomial(2 + x)            Polynomial(1 + x)
 Polynomial(0)  Polynomial(x)        Polynomial(2*x)         Polynomial(x + 2*x^2)        Polynomial(2*x + 2*x^2)
 Polynomial(0)  Polynomial(1 + x)    Polynomial(2 + 2*x)     Polynomial(1 + 2*x^2)        Polynomial(2 + x + 2*x^2)
 Polynomial(0)  Polynomial(2 + x)    Polynomial(1 + 2*x)  …  Polynomial(2 + 2*x + 2*x^2)  Polynomial(1 + 2*x^2)
 Polynomial(0)  Polynomial(2*x)      Polynomial(x)           Polynomial(2*x + x^2)        Polynomial(x + x^2)
 Polynomial(0)  Polynomial(1 + 2*x)  Polynomial(2 + x)       Polynomial(1 + x + x^2)      Polynomial(2 + x^2)
 Polynomial(0)  Polynomial(2 + 2*x)  Polynomial(1 + x)       Polynomial(2 + x^2)          Polynomial(1 + 2*x + x^2)
```
"""
function multiplication_table(degree::Int, modulo::Int)
	polys = list_polys(degree, modulo)
	number_of_polys = length(polys)
	poly_matrix = Matrix{Polynomial}(undef, number_of_polys, number_of_polys)

	for i in 1:number_of_polys, j in 1:number_of_polys
		poly_matrix[i,j] = mod(polys[i] * polys[j], modulo)
	end

	return poly_matrix
end

function __list_span_inner(modulo::Int, u̲::Vector{T}...) where T <: Number
    n_vec = length(u̲)
    span = Vector{T}[zero(u̲[1])]

    for n in 1:n_vec
        new_span = copy(span)
        for v in span, λ in 1:(modulo - 1)
            w̲ = mod.(v + λ * u̲[n], modulo)
            if w̲ ∉ new_span
                push!(new_span, w̲)
            end
        end
        span = new_span
    end

	return span
end

"""
```julia
list_span(u̲::Vector, v̲::Vector,... modulo::Int) -> Array
```

Given any number of vectors `u̲`, `v̲`,... , prints all linear combinations of those vectors, modulo `modulo`.

Parameters:
  - `u̲::Vector`: One vector.
  - `v̲::Vector`: Another vector.
  - `...`: Other vectors.
  - `modulo::Int`: The modulus of the field.

Returns:
  - `Array`: All vectors in the span of u̲ and v̲, under modulo.

---

### Examples

```julia
julia> list_span([2, 1, 1], [1, 1, 1], 3) # list the span of two vectors modulo 3
9-element Array{Array{T,1} where T,1}:
 [0, 0, 0]
 [1, 1, 1]
 [2, 2, 2]
 [2, 1, 1]
 [0, 2, 2]
 [1, 0, 0]
 [1, 2, 2]
 [2, 0, 0]
 [0, 1, 1]
```
"""
list_span(args...) = __list_span_inner(last(args), args[1:end-1]...)

"""
```julia
islinear(C::Vector, modulo::Int; verbose::Bool = false) -> Bool
```

Determines whether a code `C` is a linear code (i.e., if it is closed under addition, scalar multiplication, and has the zero vector in it).

Parameters:
  - `C::Vector`: A code, typically consisting of multiple vectors or strings.
  - `modulo::Int`: The modulus of the field under which you are working.
  - `verbose::Bool` (kwarg): print the point at which C fails, if it does.

Returns:
  - `Bool`: Whether or not the code `C` is linear (true or false).

---

### Examples

```julia

julia> islinear([[0,0,0],[1,1,1],[1,0,1],[1,1,0]], 2) # checks whether a vector of vectors is linear/a subspace (modulo 2)
false
```
"""
function islinear(C::Vector, modulo::Int; verbose::Bool = false)
	allequal_length(C) || return false # not all codes are of the same length
	block_length = length(C[1])
	𝟎 = fill(0, block_length)

	if 𝟎 ∉ C
		verbose && println("The zero vector 0̲ is not in C.\n")
		return false # the zero vector is not in the code
	end

	for c̲ ∈ C
		for λ in 0:modulo-1
			if mod.(λ*c̲, modulo) ∉ C
				verbose && println(λ, " ⋅ ", c̲, " = ", mod.(λ*c̲, modulo), " ∉ C\n")
				return false # this code isn't closed under scalar multiplication
			end
		end

		for c̲ ∈ C, c̲′ ∈ C
			if c̲ ≠ c̲′
				if mod.(c̲ + c̲′, modulo) ∉ C
					verbose && println(c̲, " + ", c̲′, " = ", mod.(c̲ + c̲′, modulo), " ∉ C\n")
					return false # this code isn't closed under addition
				end
			end
		end
	end

	verbose && println()

	return true
end

"""
```julia
isirreducible(f::AbstractPolynomial, modulo::Int) -> Bool
```

Checks if a polynomial is irreducible.

Parameters:
  - `f::Polynomial`: The polynomial you need to check.
  - `modulo::Int`: The modulus under which you are working.

Returns:
  - `Bool`: Whether or not the polynomial is irreducible (true or false).

---

### Examples

```julia
julia> isirreducible(Polynomial([1, 1, 0, 0, 1]), 2) # is 1 + x + x^4 mod 2 irreducible?
true
```
"""
function isirreducible(f::Polynomial, modulo::Int)
	deg = length(f.coeffs) - 1
	f = mod(f, deg)
	polys = list_polys(deg, modulo)

	for a in polys, b in polys
		isequal(f, mod(a*b, modulo)) && return false
	end

	return true
end

"""
```julia
normal_form!(M::AbstractArray{Int}, n::Int) -> Matrix{Int}
```

Convert a matrix M into normal form under modulo n via Gauss-Jordan elimination.  *This directly changes the matrix M.  Use `normal_form` for a non-mutating version of this function.*

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.

Returns:
  - `Matrix{Int}`: A matrix in normal form from Gauss-Jordan elimination.

---

### Examples

```julia

julia> normal_form!([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) # computes rref colswap = false
4×7 Array{Int64,2}:
1  0  0  0  2  2  2
0  1  0  0  2  0  1
0  0  1  0  1  0  2
0  0  0  1  2  2  1
```
"""
function normal_form!(M::AbstractArray{Int}, n::Int)
	return rref!(M, n, colswap = false)
end

"""
```julia
normal_form(M::AbstractArray{Int}, n::Int) -> Matrix{Int}
```

Convert a matrix M into normal form under modulo n via Gauss-Jordan elimination.

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.

Returns:
  - `Matrix{Int}`: A matrix in normal form from Gauss-Jordan elimination.

---

### Examples

```julia

julia> normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) # computes rref colswap = false
4×7 Array{Int64,2}:
1  0  0  0  2  2  2
0  1  0  0  2  0  1
0  0  1  0  1  0  2
0  0  0  1  2  2  1
```
"""
function normal_form(M::AbstractArray{Int}, n::Int)
	normal_form!(copy(M), n)
end

"""
```julia
equivalent_code!(M::AbstractArray{Int}, n::Int) -> Matrix{Int}
```

Peforms Gauss-Jordan elimination on a matrix M, but allows for column swapping.  This constructs an "equivalent" matrix.  *This directly changes the matrix M.  Use `equivalent_code` for a non-mutating version of this function.*

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.

Returns:
  - `Matrix{Int}`: A which represents an "equivalent" code to that of the matrix M.

---

### Examples

```julia
julia> equivalent_code!([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) # computes rref colswap = true
4×7 Array{Int64,2}:
1  0  0  0  2  2  2
0  1  0  0  2  0  1
0  0  1  0  1  0  2
0  0  0  1  2  2  1
```
"""
function equivalent_code!(M::AbstractArray{Int}, n::Int)
	return rref!(M, n, colswap=true)
end

"""
```julia
equivalent_code(M::AbstractArray{Int}, n::Int) -> Matrix{Int}
```

Peforms Gauss-Jordan elimination on a matrix M, but allows for column swapping.

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.

Returns:
  - `Matrix{Int}`: A which represents an "equivalent" code to that of the matrix M.

---

### Examples

```julia
julia> equivalent_code([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) # computes rref colswap = true
4×7 Array{Int64,2}:
1  0  0  0  2  2  2
0  1  0  0  2  0  1
0  0  1  0  1  0  2
0  0  0  1  2  2  1
```
"""
function equivalent_code(M::AbstractArray{Int}, n::Int)
	equivalent_code!(copy(M), n)
end

"""
```julia
generator!(M::AbstractArray{Int}, n::Int; colswap::Bool = true) -> Matrix{Int}
```

Constructs a generator matrix of the code, depending on if you allow for column swapping or not.   This function uses `normal_form!` or `equivalent_code!`.  *This directly changes the matrix M.  Use `generator` for a non-mutating version of this function.*

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.
  - `colswap::Bool` (kwarg): A boolean flag indicating whether or not you allow for swapping of columns when constructing the generating matrix.

Returns:
  - `Matrix{Int}`: A generating matrix.
"""
function generator!(M::AbstractArray{Int}, n::Int; colswap::Bool = false)
	return ifelse(colswap, equivalent_code!(M, n), normal_form!(M, n))
end

"""
```julia
generator(M::AbstractArray{Int}, n::Int; colswap::Bool = true) -> Matrix{Int}
```

Constructs a generator matrix of the code, depending on if you allow for column swapping or not.  This function uses `normal_form` or `equivalent_code`.

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.
  - `colswap::Bool` (kwarg): A boolean flag indicating whether or not you allow for swapping of columns when constructing the generating matrix.

Returns:
  - `Matrix{Int}`: A generating matrix.
"""
function generator(M::AbstractArray{Int}, n::Int; colswap::Bool = true)
	generator!(copy(M), n, colswap=colswap)
end

"""
```julia
parity_check(M::AbstractArray{Int}, n::Int) -> Matrix{Int}
```

Constructs a parity check matrix.  This is calculated from taking the non-identity part of a matrix in normal form (or equivalent — see `generator`), transposing it, multiplying it by negative one, and appending to it an appropriate sized identity matrix.

Parameters:
  - `M::AbstractArray{Int}`: A matrix of Ints.
  - `n::Int`: The modulus of the finite field.

Returns:
  - `Matrix{Int}`: A parity check matrix.
"""
function parity_check(M::AbstractArray{Int}, n::Int)
	has_identity_on_left(M) || throw(error("This matrix is not in normal form.  Use normal_form or equivalent_code."))

	minus_Dᵀ = mod.(-one(Int) .* transpose(M[:, (size(M, 1) + 1):end]), n)
	return H = Int[minus_Dᵀ I(size(minus_Dᵀ, 1))]
end

"""
```julia
syndrome(v̲::Vector, Hᵀ::AbstractArray{Int}, n::Int) -> Matrix{Int}
```

Calculates the syndrome of a given word v̲ and a parity check matrix, transposed (Hᵀ), under modulo n.

Parameters:
  - `v̲::Vector`: A word in the code.
  - `Hᵀ::AbstractArray{Int}`: The transpose of a parity check matrix.

Returns:
  - `Vector`: The syndrome of a word in the code.

---

### Examples

```julia
julia> syndrome([0, 2, 1, 2, 0, 1, 0], [1 1 1; 1 0 2; 2 0 1; 1 1 2; 1 0 0; 0 1 0; 0 0 1], 3)
1×3 Array{Int64,2}:
0  0  0

julia> syndrome([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3)
1×3 Array{Int64,2}:
 0  0  0
```
"""
function syndrome(v̲::Vector, Hᵀ::AbstractArray{Int}, n::Int)
	return mod.(v̲' * Hᵀ, n)
end

"""
```julia
isincode(v̲::Vector, Hᵀ::AbstractArray{Int}, n::Int) -> Bool
```

If the syndrome of a code is the zero vector, then the word used to calculate the syndrome is in the code.

Parameters:
  - `v̲::Vector`: A word.
  - `Hᵀ::AbstractArray{Int}`: The transpose of a parity check matrix.

Returns:
  - `Bool`: If the word is in the code or not (true or false).

---

### Examples

```julia
julia> isincode([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) # tests if the syndrome is equal to the zero vector, and is thus in the code
true
```
"""
function isincode(v̲::Vector, Hᵀ::AbstractArray{Int}, n::Int)
	return iszero(syndrome(v̲, Hᵀ, n))
end


text:
function OpenIntToFieldHashMap(arg0::Field)
    return OpenIntToFieldHashMap((Field,), arg0)
end

function OpenIntToFieldHashMap(arg0::Field, arg1::FieldElement)
    return OpenIntToFieldHashMap((Field, FieldElement), arg0, arg1)
end

function OpenIntToFieldHashMap(arg0::Field, arg1::jint)
    return OpenIntToFieldHashMap((Field, jint), arg0, arg1)
end

function OpenIntToFieldHashMap(arg0::Field, arg1::jint, arg2::FieldElement)
    return OpenIntToFieldHashMap((Field, jint, FieldElement), arg0, arg1, arg2)
end

function OpenIntToFieldHashMap(arg0::OpenIntToFieldHashMap)
    return OpenIntToFieldHashMap((OpenIntToFieldHashMap,), arg0)
end

function contains_key(obj::OpenIntToFieldHashMap, arg0::jint)
    return jcall(obj, "containsKey", jboolean, (jint,), arg0)
end

function get(obj::OpenIntToFieldHashMap, arg0::jint)
    return jcall(obj, "get", FieldElement, (jint,), arg0)
end

function iterator(obj::OpenIntToFieldHashMap)
    return jcall(obj, "iterator", OpenIntToFieldHashMap_JIterator, ())
end

function put(obj::OpenIntToFieldHashMap, arg0::jint, arg1::FieldElement)
    return jcall(obj, "put", FieldElement, (jint, FieldElement), arg0, arg1)
end

function remove(obj::OpenIntToFieldHashMap, arg0::jint)
    return jcall(obj, "remove", FieldElement, (jint,), arg0)
end

function size(obj::OpenIntToFieldHashMap)
    return jcall(obj, "size", jint, ())
end


text:
initialize{uType}(integrator,cache::OrdinaryDiffEqCache,::Type{uType}) =
                error("This algorithm does not have an initialization function")

@inline function loopheader!(integrator)
  # Apply right after iterators / callbacks

  # Accept or reject the step
  if integrator.iter > 0
    if (integrator.opts.adaptive && integrator.accept_step) || !integrator.opts.adaptive
      apply_step!(integrator)
    elseif integrator.opts.adaptive && !integrator.accept_step
      if integrator.isout
        integrator.dt = integrator.dt*integrator.opts.qmin
      else
        integrator.dt = integrator.dt/min(inv(integrator.opts.qmin),integrator.q11/integrator.opts.gamma)
      end
    end
  end

  integrator.iter += 1
  fix_dt_at_bounds!(integrator)
  modify_dt_for_tstops!(integrator)
  choose_algorithm!(integrator,integrator.cache)

end

@def ode_exit_conditions begin
  if integrator.iter > integrator.opts.maxiters
    if integrator.opts.verbose
      warn("Interrupted. Larger maxiters is needed.")
    end
    postamble!(integrator)
    integrator.sol.retcode = :MaxIters
    return integrator.sol
  end
  if !integrator.opts.force_dtmin && integrator.opts.adaptive && abs(integrator.dt) <= abs(integrator.opts.dtmin)
    if integrator.opts.verbose
      warn("dt <= dtmin. Aborting. If you would like to force continuation with dt=dtmin, set force_dtmin=true")
    end
    postamble!(integrator)
    integrator.sol.retcode = :DtLessThanMin
    return integrator.sol
  end
  if integrator.opts.unstable_check(integrator.dt,integrator.t,integrator.u)
    if integrator.opts.verbose
      warn("Instability detected. Aborting")
    end
    postamble!(integrator)
    integrator.sol.retcode = :Unstable
    return integrator.sol
  end
end

@inline function modify_dt_for_tstops!(integrator)
  tstops = integrator.opts.tstops
  if !isempty(tstops)
    if integrator.opts.adaptive
      if integrator.tdir > 0
        integrator.dt = min(abs(integrator.dt),abs(top(tstops)-integrator.t)) # step! to the end
      else
        integrator.dt = -min(abs(integrator.dt),abs(top(tstops)-integrator.t))
      end
    elseif integrator.dtcache == zero(integrator.t) && integrator.dtchangeable # Use integrator.opts.tstops
      integrator.dt = integrator.tdir*abs(top(tstops)-integrator.t)
    elseif integrator.dtchangeable # always try to step! with dtcache, but lower if a tstops
      integrator.dt = integrator.tdir*min(abs(integrator.dtcache),abs(top(tstops)-integrator.t)) # step! to the end
    end
  end
end

@inline function savevalues!(integrator::ODEIntegrator)
  while !isempty(integrator.opts.saveat) && integrator.tdir*top(integrator.opts.saveat) <= integrator.tdir*integrator.t # Perform saveat
    integrator.saveiter += 1
    curt = pop!(integrator.opts.saveat)
    if curt!=integrator.t # If <t, interpolate
      ode_addsteps!(integrator)
      Θ = (curt - integrator.tprev)/integrator.dt
      val = ode_interpolant(Θ,integrator,indices(integrator.uprev),Val{0}) # out of place, but no force copy later
      copyat_or_push!(integrator.sol.t,integrator.saveiter,curt)
      if eltype(integrator.sol.u) <: DEDataArray
        save_val = copy_non_array_fields(integrator.uprev,val)
      else
        save_val = val
      end
      copyat_or_push!(integrator.sol.u,integrator.saveiter,save_val,Val{false})
      if typeof(integrator.alg) <: OrdinaryDiffEqCompositeAlgorithm
        copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
      end
    else # ==t, just save
      copyat_or_push!(integrator.sol.t,integrator.saveiter,integrator.t)
      copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
      if typeof(integrator.alg) <: Discrete || integrator.opts.dense
        integrator.saveiter_dense +=1
        copyat_or_push!(integrator.notsaveat_idxs,integrator.saveiter_dense,integrator.saveiter)
        if integrator.opts.dense
          copyat_or_push!(integrator.sol.k,integrator.saveiter_dense,integrator.k)
        end
      end
      if typeof(integrator.alg) <: OrdinaryDiffEqCompositeAlgorithm
        copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
      end
    end
  end
  if integrator.opts.save_timeseries && integrator.iter%integrator.opts.timeseries_steps==0
    integrator.saveiter += 1
    copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
    copyat_or_push!(integrator.sol.t,integrator.saveiter,integrator.t)
    if typeof(integrator.alg) <: Discrete || integrator.opts.dense
      integrator.saveiter_dense +=1
      copyat_or_push!(integrator.notsaveat_idxs,integrator.saveiter_dense,integrator.saveiter)
      if integrator.opts.dense
        copyat_or_push!(integrator.sol.k,integrator.saveiter_dense,integrator.k)
      end
    end
    if typeof(integrator.alg) <: OrdinaryDiffEqCompositeAlgorithm
      copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
    end
  end
  resize!(integrator.k,integrator.kshortsize)
end

@inline function postamble!(integrator)
  solution_endpoint_match_cur_integrator!(integrator)
  resize!(integrator.sol.t,integrator.saveiter)
  resize!(integrator.sol.u,integrator.saveiter)
  resize!(integrator.sol.k,integrator.saveiter_dense)
  !(typeof(integrator.prog)<:Void) && Juno.done(integrator.prog)
end

@inline function solution_endpoint_match_cur_integrator!(integrator)
  if integrator.sol.t[integrator.saveiter] !=  integrator.t
    integrator.saveiter += 1
    copyat_or_push!(integrator.sol.t,integrator.saveiter,integrator.t)
    copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
    if typeof(integrator.alg) <: Discrete || integrator.opts.dense
      integrator.saveiter_dense +=1
      copyat_or_push!(integrator.notsaveat_idxs,integrator.saveiter_dense,integrator.saveiter)
      if integrator.opts.dense
        copyat_or_push!(integrator.sol.k,integrator.saveiter_dense,integrator.k)
      end
    end
  end
end

@inline function loopfooter!(integrator)
  if integrator.opts.adaptive
    integrator.q11 = integrator.EEst^integrator.opts.beta1
    q = integrator.q11/(integrator.qold^integrator.opts.beta2)
    q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),q/integrator.opts.gamma))
    dtnew = integrator.dt/q
    ttmp = integrator.t + integrator.dt
    integrator.isout = integrator.opts.isoutofdomain(ttmp,integrator.u)
    integrator.accept_step = (!integrator.isout && integrator.EEst <= 1.0) || (integrator.opts.force_dtmin && abs(integrator.dt) <= abs(integrator.opts.dtmin))
    if integrator.accept_step # Accept
      integrator.tprev = integrator.t
      integrator.t = ttmp
      integrator.qold = max(integrator.EEst,integrator.opts.qoldinit)
      calc_dt_propose!(integrator,dtnew)
      handle_callbacks!(integrator)
    end
  else #Not adaptive
    integrator.tprev = integrator.t
    integrator.t += integrator.dt
    integrator.accept_step = true
    integrator.dtpropose = integrator.dt
    handle_callbacks!(integrator)
  end
  if !(typeof(integrator.prog)<:Void) && integrator.opts.progress && integrator.iter%integrator.opts.progress_steps==0
    Juno.msg(integrator.prog,integrator.opts.progress_message(integrator.dt,integrator.t,integrator.u))
    Juno.progress(integrator.prog,integrator.t/integrator.sol.prob.tspan[2])
  end
end

@inline function handle_callbacks!(integrator)
  discrete_callbacks = integrator.opts.callback.discrete_callbacks
  continuous_callbacks = integrator.opts.callback.continuous_callbacks
  atleast_one_callback = false

  continuous_modified = false
  discrete_modified = false
  if !(typeof(continuous_callbacks)<:Tuple{})
    time,upcrossing,idx,counter = find_first_continuous_callback(integrator,continuous_callbacks...)
    if time != zero(typeof(integrator.t)) && upcrossing != 0 # if not, then no events
      atleast_one_callback = true
      continuous_modified = apply_callback!(integrator,continuous_callbacks[idx],time,upcrossing)
    end
  end
  if !(typeof(discrete_callbacks)<:Tuple{})
    atleast_one_callback = true
    discrete_modified = apply_discrete_callback!(integrator,discrete_callbacks...)
  end
  if !atleast_one_callback
    savevalues!(integrator)
  end

  integrator.u_modified = continuous_modified || discrete_modified
  if integrator.u_modified
    handle_callback_modifiers!(integrator)
  end
end

@inline function handle_callback_modifiers!(integrator::ODEIntegrator)
  integrator.reeval_fsal = true
end

@inline function apply_step!(integrator)

  integrator.accept_step = false # yay we got here, don't need this no more

  #Update uprev
  if alg_extrapolates(integrator.alg)
    if typeof(integrator.u) <: AbstractArray
      recursivecopy!(integrator.uprev2,integrator.uprev)
    else
      integrator.uprev2 = integrator.uprev
    end
  end
  if typeof(integrator.u) <: AbstractArray
    recursivecopy!(integrator.uprev,integrator.u)
  else
    integrator.uprev = integrator.u
  end

  #Update dt if adaptive or if fixed and the dt is allowed to change
  if integrator.opts.adaptive || integrator.dtchangeable
    integrator.dt = integrator.dtpropose
  elseif integrator.dt != integrator.dtpropose && !integrator.dtchangeable
    error("The current setup does not allow for changing dt.")
  end

  # Update fsal if needed
  if isfsal(integrator.alg)
    if !isempty(integrator.opts.d_discontinuities) && top(integrator.opts.d_discontinuities) == integrator.t
      pop!(integrator.opts.d_discontinuities)
      reset_fsal!(integrator)
    elseif integrator.reeval_fsal || (typeof(integrator.alg)<:DP8 && !integrator.opts.calck) || (typeof(integrator.alg)<:Union{Rosenbrock23,Rosenbrock32} && !integrator.opts.adaptive)
      reset_fsal!(integrator)
    else # Do not reeval_fsal, instead copy! over
      if typeof(integrator.fsalfirst) <: AbstractArray
        recursivecopy!(integrator.fsalfirst,integrator.fsallast)
      else
        integrator.fsalfirst = integrator.fsallast
      end
    end
  end
end

@inline function calc_dt_propose!(integrator,dtnew)
  integrator.dtpropose = integrator.tdir*min(abs(integrator.opts.dtmax),abs(dtnew))
  integrator.dtpropose = integrator.tdir*max(abs(integrator.dtpropose),abs(integrator.opts.dtmin))
end

@inline function fix_dt_at_bounds!(integrator)
  if integrator.tdir > 0
    integrator.dt = min(integrator.opts.dtmax,integrator.dt)
  else
    integrator.dt = max(integrator.opts.dtmax,integrator.dt)
  end
  if integrator.tdir > 0
    integrator.dt = max(integrator.dt,integrator.opts.dtmin) #abs to fix complex sqrt issue at end
  else
    integrator.dt = min(integrator.dt,integrator.opts.dtmin) #abs to fix complex sqrt issue at end
  end
end

@inline function handle_tstop!(integrator)
  tstops = integrator.opts.tstops
  if !isempty(tstops)
    t = integrator.t
    ts_top = top(tstops)
    if t == ts_top
      pop!(tstops)
      integrator.just_hit_tstop = true
    elseif integrator.tdir*t > integrator.tdir*ts_top
      if !integrator.dtchangeable
        change_t_via_interpolation!(integrator, pop!(tstops), Val{true})
        integrator.just_hit_tstop = true
      else
        error("Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error.")
      end
    end
  end
end

@inline function reset_fsal!(integrator)
  # Under these condtions, these algorithms are not FSAL anymore
  if typeof(integrator.cache) <: OrdinaryDiffEqMutableCache
    integrator.f(integrator.t,integrator.u,integrator.fsalfirst)
  else
    integrator.fsalfirst = integrator.f(integrator.t,integrator.u)
  end
  integrator.reeval_fsal = false
end

function (integrator::ODEIntegrator)(t,deriv::Type=Val{0};idxs=size(integrator.uprev))
  current_interpolant(t,integrator,idxs,deriv)
end

(integrator::ODEIntegrator)(val::AbstractArray,t::Union{Number,AbstractArray},deriv::Type=Val{0};idxs=eachindex(integrator.uprev)) = current_interpolant!(val,t,integrator,idxs,deriv)


text:
using SpecialFunctions

n = 2000

probAgivenB(k) = 1- 1/(k+1)
probB(k) = 6/(pi*(k+1))^2

numerical= sum([probAgivenB(k)*probB(k) for k in 0:n])
analytic = 1 - 6*zeta(3)/pi^2

println("Analytic: ", analytic, "\tNumerical: ", numerical)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
export QMatrix, subtract!, getcolumn

"""

    QMatrix{T} <: AbstractMatrix{T}

Efficient implementation of a matrix with mixed binary and non-binary
entries. Binary entries are stored in a dense BitMatrix whereas
non-binary columns are stored as separate vectors.

"""
mutable struct QMatrix{T} <: AbstractMatrix{T}
    m::Int # number of rows
    n::Int # number of columns
    binary::BitMatrix # binary data
    qary::Dict{Int,Vector{T}} # qary data
    function QMatrix{T}(m::Int, n::Int) where T
        if m % 64 != 0 # allows for working over the BitMatrix chunks directly
            throw(ArgumentError("the number of QMatrix rows must be a multiple of 64"))
        end
        new(m, n, falses(m, n), Dict{Int,Vector{T}}())
    end
end

## for finding chunk indices. copied from the BitArray module ##
@inline _div64(l) = l >> 6
@inline _mod64(l) = l & 63
@inline get_chunks_id(i::Integer) = _div64(Int(i)-1)+1, _mod64(Int(i)-1)

function Base.size(M::QMatrix)
    return M.m, M.n
end

"""
    getindex{T}(M::QMatrix{T}, r::Int, c::Int)

Implements M[i,j]

"""
function Base.getindex(M::QMatrix{T}, r::Int, c::Int) where T
    @boundscheck checkbounds(M.binary, r, c)
    if haskey(M.qary, c)
        return M.qary[c][r]
    end
    return T(M.binary[r,c])
end

"""
    setindex!{T}(M::QMatrix{T}, d::T, r::Int, c::Int)

Implements M[i,j] = d

"""
function Base.setindex!(M::QMatrix{T}, d::T, r::Int, c::Int) where T
    @boundscheck checkbounds(M.binary, r, c)
    if iszero(d) || d == one(T)
        if haskey(M.qary, c)
            M.qary[c][r] = d
        else
            M.binary[r,c] = d
        end
    else
        if !haskey(M.qary, c)
            M.qary[c] = Vector{T}(M.binary[:,c])
        end
        M.qary[c][r] = d
    end
    return d
end

"""
    setindex!{T}(M::QMatrix{T}, d::T, r::Int, c::Int)

Implements M[i,j] = d

"""
function Base.setindex!(M::QMatrix{T}, d::T, ::Colon, c::Int) where T
    @boundscheck checkbounds(M.binary, 1, c)
    if iszero(d) || d == one(T)
        delete!(M.qary, c)
        M.binary[:,c] .= d
        return d
    else
        if !haskey(M.qary, c)
            M.qary[c] = Vector{T}(M.binary[:,c])
        end
        M.qary[c][:] .= d
    end
    return d
end

"""
    resize!(M::QMatrix{T}, m, n)

Resize the matrix to have m rows and n columns. The given m and n must
be larger than the original values.

"""
function Base.resize!(M::QMatrix{T}, m, n) where T
    if m < M.m
        error("the given m must not be smaller than its original value")
    end
    if n < M.n
        error("the given n must not be smaller than its original value")
    end
    if m % 64 != 0 # allows for working over the BitMatrix chunks directly
        throw(ArgumentError("the number of QMatrix rows must be a multiple of 64"))
    end
    for v in values(M.qary)
        resize!(v, m)
        for i in M.m+1:m # zero out the new values
            v[i] = zero(T)
        end
    end
    binary = falses(m, n)
    cm, cmo = _div64(m), _div64(M.m)
    for i in 1:length(M.binary.chunks)
        j = cm*(div(i-1, cmo)) + i % cmo + 1
        binary.chunks[j] = M.binary.chunks[i]
    end
    M.binary = binary
    M.m, M.n = m, n
    return M
end

"""
    getcolumn{T}(M::QMatrix{T}, c::Int)

Return the c-th column of the matrix.

"""
function getcolumn(M::QMatrix{T}, c::Int) where T
    @boundscheck checkbounds(M.binary, 1, c)
    if haskey(M.qary, c)
        return M.qary[c]
    end
    return Vector{T}(M.binary[:,c])
end

"""
    getcolumn!{T}(v::Vector{T}, M::QMatrix{T}, c::Int)

Return the c-th column of the matrix in-place.

"""
function getcolumn!(v::Vector{T}, M::QMatrix{T}, c::Int) where T
    @boundscheck checkbounds(M.binary, 1, c)
    if haskey(M.qary, c)
        v[:] .= M.qary[c]
        return v
    end
    v[:] .= M.binary[:,c]
    return v
end

"""
    countnz(M::QMatrix, c::Int)

Return the number of non-zero entries in the c-th column of M.

TODO: no longer part of the standard library

"""
function countnz(M::QMatrix, c::Int)
    if haskey(M.qary, c)
        n = 0
        for d in M.qary[c]
            if !iszero(d)
                n += 1
            end
        end
        return n
    else
        n = sum(M.binary[:,c])
        return n
    end
end

"""
    isbinary(M::QMatrix)

Return true if the c-th column is all-binary.

"""
function isbinary(M::QMatrix, c::Int)
    @boundscheck checkbounds(M.binary, 1, c)
    return !haskey(M.qary, c)
end

"""
    subtract!(M::QMatrix{T}, c1::Int, c2::Int)

Subtract column c2 from column c1, i.e., M[:,c1] = M[:,c1] - M[:,c2].

"""
function subtract!(M::QMatrix{T}, c1::Int, c2::Int) where T
    @boundscheck checkbounds(M.binary, 1, c1)
    @boundscheck checkbounds(M.binary, 1, c2)
    h1, h2 = haskey(M.qary, c1), haskey(M.qary, c2)
    if h1 && !h2
        @views M.qary[c1] .-= M.binary[:, c2]
    elseif !h1 && h2
        @views M.qary[c1] = M.binary[:,c1] .- M.qary[c2]
    elseif h1 && h2
        @views q1, q2 = M.qary[c1], M.qary[c2]
        q1 .-= q2
    else
        s2i = LinearIndices(M)
        kd0, ld0 = get_chunks_id(s2i[1, c1])
        kd1, ld1 = get_chunks_id(s2i[M.m, c1])
        offset = (c2-c1)*div(M.m,64)
        @simd for i in kd0:kd1
            @inbounds M.binary.chunks[i] = xor(
                M.binary.chunks[i],
                M.binary.chunks[i+offset],
            )
        end
    end
    return
end

"""
    subtract!(M::QMatrix{T}, d::T, c1::Int, c2::Int)

Subtract column d*c2 from column c1, i.e., M[:,c1] .-= d*M[:,c2].

"""
function subtract!(M::QMatrix{T}, d::T, c1::Int, c2::Int) where T
    if iszero(d)
        return
    end
    if d == one(d)
        return subtract!(M, c1, c2)
    end
    h1, h2 = haskey(M.qary, c1), haskey(M.qary, c2)
    if h1 && !h2
        @views M.qary[c1] .-= d .* M.binary[:, c2]
    elseif !h1 && h2
        @views M.qary[c1] = M.binary[:,c1] .- d.*M.qary[c2]
    elseif h1 && h2
        @views q1, q2 = M.qary[c1], M.qary[c2]
        q1 .-= d.*q2
    else
        @views M.qary[c1] = M.binary[:,c1] .- d.*M.binary[:,c2]
    end
    return
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
(+)(td::TimeDate) = td
(+)(tdz::TimeDateZone) = tdz

function (-)(td1::TimeDate, td2::TimeDate)
    dtm1 = DateTime(td1)
    dtm2 = DateTime(td2)
    dtm0 = canonical(Microsecond(td1) - Microsecond(td2) + Nanosecond(td1) - Nanosecond(td2))
    dtm = dtm1 - dtm2
    dtm += dtm0
    return sum(fldmod(dtm))
end

(-)(td::TimeDate, dt::DateTime) = td - TimeDate(dt)
(-)(dt::DateTime, td::TimeDate) = TimeDate(dt) - td
(-)(td::TimeDate, dt::Date) = td - TimeDate(dt)
(-)(dt::Date, td::TimeDate) = TimeDate(dt) - td

function (-)(td::TimeDate, tm::Time)
    tm_td, dt_td = at_time(td), on_date(td)
    if tm_td < tm
        dt_td -= Day(1)
        tm_td = tm_td - tm + Hour(24)
    else
        tm_td = tm_td - tm
    end
    tm_td = canonical(tm_td)
    tm = Time(tm_td)
    return TimeDate(tm, dt_td)
end

for P in (:Year, :Month, :Day)
    @eval begin
        function (+)(td::TimeDate, p::$P)
            tm, dt = at_time(td), on_date(td)
            dt = dt + p
            return TimeDate(tm, dt)
        end
        function (-)(td::TimeDate, p::$P)
            p = -p
            return td + p
        end
    end
end

for P in (:Hour, :Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)
    @eval begin
        function (+)(td::TimeDate, p::$P)
            tm, dt = at_time(td), on_date(td)
            cptm = canonical(CompoundPeriod(tm) + p)
            extradays = Day(cptm)
            cptm -= extradays
            dt += extradays
            tim = Time(cptm)
            return TimeDate(tim, dt)
        end
        function (-)(td::TimeDate, p::$P)
            p = -p
            return td + p
        end
    end
end

function(+)(tdz::TimeDateZone, p1::Period, p2::Period)
    p = p1 + p2
    return tdz + p
end

function(-)(tdz::TimeDateZone, p1::Period, p2::Period)
    p = p1 + p2
    return tdz - p
end

function(+)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period)
    p = p1 + p2 + p3
    return tdz + p
end

function(-)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period)
    p = p1 + p2 + p3
    return tdz - p
end

function(+)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period, p4::Period)
    p = p1 + p2 + p3 + p4
    return tdz + p
end

function(-)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period, p4::Period)
    p = p1 + p2 + p3 + p4
    return tdz - p
end

function(+)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period, p4::Period, p5::Period)
    p = p1 + p2 + p3 + p4 + p5
    return tdz + p
end

function(-)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period, p4::Period, p5::Period)
    p = p1 + p2 + p3 + p4 + p5
    return tdz - p
end

function(+)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period, p4::Period, p5::Period, p6::Period)
    p = p1 + p2 + p3 + p4 + p5 + p6
    return tdz + p
end

function(-)(tdz::TimeDateZone, p1::Period, p2::Period, p3::Period, p4::Period, p5::Period, p6::Period)
    p = p1 + p2 + p3 + p4 + p5 + p6
    return tdz - p
end


text:
@testset "struct ClustData" begin
    @test 1==1
    # constructors are all tested by run_clust
    # TODO: add unit tests for constructors anyways.
end

@testset "SimpleExtremeValueDescr" begin
    @testset "normal" begin
        ev1 = SimpleExtremeValueDescr("wind-dena42","max","absolute")
        @test ev1.data_type =="wind-dena42"
        @test ev1.extremum =="max"
        @test ev1.peak_def =="absolute"
        @test ev1.consecutive_periods ==1
    end
    @testset "edge cases" begin
        @test_throws ErrorException ev1 = SimpleExtremeValueDescr("wind-dena42","maximum","absolute")
        @test_throws ErrorException ev1 = SimpleExtremeValueDescr("wind-dena42","max","abs")
    end
end


text:
module MCCGA

import Optim
import Statistics

export Functions

export bits
export floats
export initialprobs
export sample
export mccga

include("stats.jl")
include("bitworks.jl")
include("thealgorithm.jl")
include("funcset.jl")

end # module

max_stars_repo_name: jbytecode/MCCGA

max_issues_repo_name: jbytecode/MCCGA

max_forks_repo_name: jbytecode/MCCGA


text:
module MYASMrun1

using Modia
@usingModiaPlot
using StaticArrays

include("../src/Utils.jl")
include("../src/SIM8.jl")

# logData = true

ZoneDataMap = ReadDataSet("ZONEdata.yaml")

instModel = Main.MYASM1.instModel

include("ASMinit.start.jl")

function mySimulation(;reuse=false)
    printstyled("\nNEW simulation\n",color=:light_blue)

    simulate!(instModel,
#         IDA(),
        CVODE_BDF(),
        startTime=-0.1,
        stopTime = 1.5,
        interval=2e-4,
        tolerance=1e-6,
        log=false,
        logStates=true,
        logEvents=false,
        logTiming = false,
        logParameters    = false,
        logEvaluatedParameters   = false,

        merge=Map(
            transformer = ZoneDataMap["setGEN"],
            meterGEN    = ZoneDataMap["setGEN"],
            varvolts    = Map(Tstart=1e-3, Vmin=0.3),     # Event muss T>0.0

        ),
    )

     plot(instModel, [( "asmUnit.current"),("asmUnit.tau")], figure=1)
     plot(instModel, [ ("J.w")], figure=2)
#      plot(instModel, [("synUnit.SM.i"),("GridVolts.volts")], figure=2)
#       plot(instModel, [ ("GridVolts.amps")], figure=2)

end

# open("LOG.txt", "a") do io
#     redirect_stdout(io) do
        mySimulation()
#     end
# end

Ergebnis(instModel, [
    "asmUnit.Psi1",
    "asmUnit.Psi2",
    "asmUnit.voltage",
    "meterGEN.power",
    "meterGRID.power",
    "J.w"])

printstyled(length(instModel.result_x[1,:]),"\n", color=:red)

# Speichern(instModel, ["time","synUnit.exciter.y", "synUnit.SM.ifd", "synUnit.SM.tau", "meterGEN.voltage", "meterGEN.current", "meterGRID.voltage", "synUnit.SM.phiKR"])

end

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
using InstantiateFromURL

# activate the QuantEcon environment
activate_github("QuantEcon/QuantEconLecturePackages", tag = "v0.9.5");

max_stars_repo_name: fratrik/lecture-source-jl

max_issues_repo_name: fratrik/lecture-source-jl

max_forks_repo_name: fratrik/lecture-source-jl


text:
# This file is a part of BinaryCommitteeMachineFBP.jl. License is MIT: http://github.com/carlobaldassi/BinaryCommitteeMachineFBP.jl/LICENCE.md

module Util

export exclusive, checkdims, chgeltype,
       @readmagvec, @dumpmagvecs,
       IVec, Vec, Vec2, MagVec, MagVec2, MagVec3

using ..Magnetizations: Mag64, showinner, parseinner

if success(`which lockfile`)
    function exclusive(f::Function, fn::AbstractString = "lock.tmp")
        run(`lockfile -1 $fn`)
        try
            f()
        finally
            run(`rm -f $fn`)
        end
    end
else
    @warn("could not find the `lockfile` program. Try installing it (it's included with `procmail` – www.procmail.org).")
    exclusive(f::Function, fn::AbstractString) = nothing
end

annotate_new!(T, ex) = ex
function annotate_new!(T, ex::Expr)
    if ex.head == :call && ex.args[1] == :new
        ex.args[1] = Expr(:curly, :new, T)
    end
    map!(x->annotate_new!(T, x), ex.args, ex.args)
    return ex
end

const IVec = Vector{Int}
const Vec = Vector{Float64}
const Vec2 = Vector{Vec}
const MagVec{M<:Mag64} = Vector{M}
const MagVec2{M<:Mag64} = Vector{MagVec{M}}
const MagVec3{M<:Mag64} = Vector{MagVec2{M}}

function checkdims(x::AbstractArray{<:AbstractArray}, N::Integer, r::Integer...)
    @assert length(x) == N
    @assert all(v->checkdims(v, r...), x)
    return true
end
function checkdims(x::AbstractArray, N::Integer)
    @assert length(x) == N
    return true
end
checkdims(x::AbstractArray{<:AbstractArray}, N::Integer) = error("too few dimensions")
checkdims(x::AbstractArray, N::Integer, r::Integer...) = error("not enough dimensions")

changedeltype(::Type{Vector{V}}, ::Type{T}) where {V<:Vector,T} = Vector{changedeltype(V, T)}
changedeltype(::Type{Vector{X}}, ::Type{T}) where {X,T} = Vector{T}

function chgeltype(x::Vector{<:Vector}, ::Type{T}) where {T}
    eltype(changedeltype(typeof(x), T))[chgeltype(y, T) for y in x]
end
chgeltype(x::Vector, ::Type{T}) where {T} = convert(Vector{T}, x)

macro readmagvec(l, fmt, vs...)
    ex = :()
    for v in vs
        vn = Regex(string(v, "(\\[|\\s)"))
        ex = quote
            $ex
            if occursin($vn, $(esc(l)))
                _readmagvec($(esc(l)), $(esc(fmt)), $(esc(v)))
                @goto found
            end
        end
    end
    ex = quote
        $ex
        error("unrecognized line: ", $(esc(l)))
        @label found
    end
    ex
end

function _readmagvec(l::AbstractString, fmt::Type{Val{F}}, v::Array) where {F}
    sl = split(l)
    pre = sl[1]
    occursin(r"^[^[]+(?:\[\d+\])*$", pre) || error("invalid messages file")
    inds = map(x->parse(Int,x), split(replace(pre, r"[][]" => " "))[2:end])
    for i in inds
        v = v[i]
    end
    length(v) == length(sl) - 1 || error("invalid messages file")
    for i = 1:length(v)
        v[i] = parseinner(fmt, sl[i+1])
    end
end

macro dumpmagvecs(io, vs...)
    ex = :()
    for v in vs
        vn = string(v)
        ex = quote
            $ex
            _dumpmagvec($(esc(io)), $(esc(v)), $vn)
        end
    end
    ex
end

function _dumpmagvec(io::IO, a::Array{<:Array}, s::AbstractString)
    for (i,v) in enumerate(a)
        _dumpmagvec(io, v, "$s[$i]")
    end
end
function _dumpmagvec(io::IO, a::Array, s::AbstractString)
    print(io, s)
    for x in a
        print(io, ' ')
        showinner(io, x)
    end
    println(io)
end

end


text:
"""DNN Feature decoding - feature prediction script"""

using Statistics

using BrainDecoder
using FastL2LiR
using HDF5
using Printf
using Glob
using Debugger

# Settings ###################################################################

# Brain data
brain_dir = "/home/kiss/data/fmri_shared/datasets/Deeprecon/fmriprep"
subjects = Dict(
    "AM" => "AM_ImageNetTest_volume_native.h5",
    # "ES" => "ES_ImageNetTest_volume_native.h5",
    # "FC" => "FC_ImageNetTest_volume_native.h5",
    # "JK" => "JK_ImageNetTest_volume_native.h5",
    # "JP" => "JP_ImageNetTest_volume_native.h5",
    # "KS" => "KS_ImageNetTest_volume_native.h5",
    # "TH" => "TH_ImageNetTest_volume_native.h5",
)

label_key = "stimulus_name"

rois = Dict(
    "VC"  => "ROI_VC",
    # "LVC" => "ROI_LVC",
    # "HVC" => "ROI_HVC",
    # "V1"  => "ROI_V1",
    # "V2"  => "ROI_V2",
    # "V3"  => "ROI_V3",
    # "V4"  => "ROI_hV4",
    # "LOC" => "ROI_LOC",
    # "FFA" => "ROI_FFA",
    # "PPA" => "ROI_PPA",
)

# DNN features
features_dir = "/home/nu/data/contents_shared/ImageNetTest/derivatives/features/default"
network = "caffe/bvlc_alexnet"
layers = [
    "conv1", #"relu1", "pool1",
    "conv2", #"relu2", "pool2",
    "conv3", #"relu3",
    "conv4", #"relu4",
    "conv5", #"relu5", "pool5",
    "fc6", #"relu6",
    "fc7", #"relu7",
    "fc8",
]

chunk_axis = 0
# Note that Julia's indexing is one-based.

# Feature decoders
models_dir_root = "./results/feature_decoders/deeprecon_fmriprep_rep5_500voxel_allunits_fastl2lir_alpha100"

# Results directory
#results_dir_root = "/home/nu/data/contents_shared/ImageNetTraining/derivatives/feature_decoders/deeprecon_fmriprep_rep5_500voxel_allunits_fastl2lir_alpha100"
results_dir_root = "./results/decoded_features/deeprecon_fmriprep_rep5_500voxel_allunits_fastl2lir_alpha100"

# Functions ##################################################################

function predict_features(model_dir, x; x_mean=[], x_norm=[], y_mean=[], y_norm=[], chunk_axis=0)

    # Normalize X
    x = (x .- x_mean) ./ x_norm

    model = load_model(model_dir, chunk_axis=chunk_axis)

    @show size(model.W)
    @show size(model.b)

    y = predict(model, x)

    # Denormalize Y
    y = y .* y_norm .+ y_mean

    return y
end

# Main #######################################################################

# Load data ------------------------------------------------------------------
features = Features(joinpath(features_dir, network))

# Analysis loop --------------------------------------------------------------
println("----------------------------------------")
println("Analysis loop")

for lay in layers

    # Load DNN features
    println("Loading features: " * lay)
    y = features.get_features(lay)  # TODO: speed-up
    y_labels = features.labels

    y = Float32.(y)

    @show size(y)
    @show size(y_labels)

    for sub in keys(subjects)

        # Load fMRI data
        brain_data_file = joinpath(brain_dir, subjects[sub])
        println("Loading " * brain_data_file)
        brain_data = BData(brain_data_file)

        for roi in keys(rois)
            println("Layer:   " * lay)
            println("Subject: " * sub)
            println("ROI:     " * roi)

            # Feature prediction
            save_dir = joinpath(results_dir_root, "decoded_features", network, lay, sub, roi)
            if !isdir(save_dir)
                mkpath(save_dir)
            end
            # TODO: check existing files

            model_dir = joinpath(models_dir_root, network, lay, sub, roi, "model")

            x = brain_data.select(rois[roi])
            x_labels = brain_data.get_labels(label_key)

            x = Float32.(x)

            @show size(x)
            @show size(x_labels)

            # Average fMRI samples
            x_labels_unique = unique(x_labels)
            x_ave = cat([mean(x[findall(x -> x == xl, x_labels), :], dims=1)
                         for xl in unique(x_labels)]..., dims=1)

            @show size(x_ave)
            @show size(x_labels_unique)

            # Preprocessing parameters
            x_mean = load_array(joinpath(model_dir, "x_mean.mat"), "x_mean")
            x_norm = load_array(joinpath(model_dir, "x_norm.mat"), "x_norm")
            y_mean = load_array(joinpath(model_dir, "y_mean.mat"), "y_mean")
            y_norm = load_array(joinpath(model_dir, "y_norm.mat"), "y_norm")

            # Get sample index of Y
            y_index = [findall(x -> x == xl, y_labels)[1] for xl in x_labels_unique]

            # Call training routine
            @time y_pred = predict_features(
                model_dir, x_ave,
                x_mean=x_mean, x_norm=x_norm,
                y_mean=y_mean, y_norm=y_norm,
                chunk_axis=chunk_axis
            )

            @show size(y_pred)

            # Save decoded features
            y_idx = axes(y_pred)
            for (i, xl) in enumerate(x_labels_unique)
                save_file = joinpath(save_dir, xl * ".mat")
                save_array(save_file, "feat", y_pred[[i], y_idx[2:end]...])
            end
        end
    end
end

println("All done")


text:
#!/usr/bin/julia
using AssociationsGUI, Gtk
f = open_dialog("")
folder, _ = splitdir(f)
#win = Gtk.Window("")
#folder = open_dialog("Select videos-folder", action=Gtk.GtkFileChooserAction.SELECT_FOLDER)
main(folder)


text:
struct Factor{T<:Unsigned} <: AbstractFactor{T}
    name::String
    levels::AbstractVector{<:AbstractString}
    data::AbstractVector{T}
end

Base.length(factor::Factor{T}) where {T<:Unsigned} = length(factor.data)

function Factor{T}(name::String, length::Integer, levels::Vector{String}, datpath::String) where {T<:Unsigned}
    data = Vector{T}(length)
    open(datpath) do f
        read!(f, data)
    end
    Factor{T}(name, levels, data)
end

function Factor(name::String, data::AbstractVector{<:AbstractString})
    maxlevelcount = length(unique(data))
    if maxlevelcount <= typemax(UInt8) + 1
        levelmap = Dict{String, UInt8}()
        levels = Vector{String}()
        resdata = Vector{UInt8}()
        for v in data
            if v == "" || v == MISSINGLEVEL
                push!(resdata, 0)
            else
                levelcount = length(levelmap)
                levelindex = get(levelmap, v, levelcount + 1)
                if levelindex > levelcount
                    push!(levels, v)
                end
                push!(resdata, levelindex)
            end
        end
        Factor{UInt8}(name, levels, resdata)
    elseif maxlevelcount <= typemax(UInt16) + 1
        levelmap = Dict{String, UInt16}()
        levels = Vector{String}()
        resdata = Vector{UInt16}()
        for v in data
            if v == "" || v == MISSINGLEVEL
                push!(resdata, 0)
            else
                levelcount = length(levelmap)
                levelindex = get(levelmap, v, levelcount + 1)
                if levelindex > levelcount
                    push!(levels, v)
                end
                push!(resdata, levelindex)
            end
        end
        Factor{UInt16}(name, levels, resdata)
    else
        levelmap = Dict{String, UInt32}()
        levels = Vector{String}()
        resdata = Vector{UInt32}()
        for v in data
            if v == "" || v == MISSINGLEVEL
                push!(resdata, 0)
            else
                levelcount = length(levelmap)
                levelindex = get(levelmap, v, levelcount + 1)
                if levelindex > levelcount
                    push!(levels, v)
                end
                push!(resdata, levelindex)
            end
        end
        Factor{UInt32}(name, levels, resdata)
    end
end

function slice(factor::Factor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    slice(factor.data, fromobs, toobs, slicelength)
end

function Base.map(factor::Factor, dataframe::AbstractDataFrame)
    dataframe[getname(factor)]
    #PermuteFactor(factor, dataframe[getname(factor)])
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using TraitWrappers
import Base: IteratorEltype, HasEltype

fn_holy(itr) = fn_holy(IteratorEltype(itr), itr)

fn_holy(::HasEltype, itr) = begin
	println("this itr has `eltype()` defined")
	println(itr)
end

fn_holy(1)
fn_tw(1)