text:
using Simulacrum
using HyperGraphs, Catalyst

# setting up some Catalyst Reaction and ChemicalHyperEdge variables
# r1 and r2 are equivalent to che1 and che2, respectively
@variables t X(t) Y(t) Z(t)
r1 = Reaction(2, [X], [Y])					# Catalyst reaction unimolecular
r2 = Reaction(2, [X, Y], [Z])				# Catalyst reaction multimolecular
che1 = ChemicalHyperEdge([X], [Y], 2)		# chemical hyperedge unimolecular
che2 = ChemicalHyperEdge([X, Y], [Z], 2)	# chemical hyperedge multimolecular
@named rsys = ReactionSystem([r1, r2], t, [X, Y, Z], [])
chx = ChemicalHyperGraph([che1, che2])

# testing equivalence
@test incidence_matrix(ChemicalHyperGraph(r1)) == incidence_matrix(ChemicalHyperGraph(che1))
@test incidence_matrix(chx) == incidence_matrix(ChemicalHyperGraph(rsys))

# conversions
@test convert(ChemicalHyperEdge, r1) == che1 && convert(ChemicalHyperEdge, r2) == che2
rx = [Reaction(1, nothing, [X]), Reaction(1, [Y], nothing)]
@test typeof(convert(ChemicalHyperEdge, rx)) == Vector{ChemicalHyperEdge{eltype(rx[1].products)}}
@test typeof(convert(Reaction, che1).substrates) == Vector{eltype(r1.substrates)}
@test eltype(convert(ChemicalHyperGraph, rsys)) == eltype(equations(rsys)[1].substrates)
catalyst_reactions(chx)
convert(ReactionSystem, chx)
ChemicalHyperGraph(r1)
ChemicalHyperGraph([r1, r2])
ChemicalHyperGraph(rsys)
convert(ODESystem, convert(ReactionSystem, chx))

# these convert reactions to chemical hyperedges under the hood
@test isunimolecular(r1) && ismultimolecular(r2)
@test !isunimolecular(r2) && !ismultimolecular(r1)


text:
module FieldOperationsTests

using BlockArrays
using Gridap.Arrays
using Gridap.Fields
using Gridap.Fields: IntKernel
using Gridap.Fields: FieldOpKernel
using Gridap.Fields: trialize_basis_value
using Gridap.TensorValues
using Gridap.Fields: MockField, MockBasis

using Test
using FillArrays
using LinearAlgebra

s = field_operation_metasize((),(:,),(1,:))
@test s == (:,:)

s = field_operation_metasize((),(),(:,))
@test s == (:,)

s = field_operation_metasize((1,:),(1,:))
@test s == (1,:)

ndofs1 = 5
ndofs2 = 4
np = 7
ax1 = (Base.OneTo(ndofs1),)
ax2 = (1,Base.OneTo(ndofs2),)
ax3 = ()
ax = field_operation_axes(ax1,ax2,ax3)
@test ax == (Base.OneTo(ndofs1),Base.OneTo(ndofs2))

ax1 = (Base.OneTo(np),Base.OneTo(ndofs1),)
ax2 = (Base.OneTo(np),1,Base.OneTo(ndofs2),)
ax3 = (Base.OneTo(np),)
ax = field_operation_axes(ax3,ax2,ax1)
@test ax == (Base.OneTo(np),Base.OneTo(ndofs1),Base.OneTo(ndofs2))

ax1 = (blockedrange([ndofs1,ndofs2]),)
ax2 = (blockedrange([1]),blockedrange([ndofs2,ndofs1]))
ax3 = ()
ax = field_operation_axes(ax3,ax2,ax1)
@test ax == (blockedrange([ndofs1,ndofs2]),blockedrange([ndofs2,ndofs1]))

l = 10
np = 3
ndofs = 4

# Operations between fields and array of fields

p = Point(1,2)
x = fill(p,np)
z = 2.0

v = VectorValue(3.0,1.5)
w = VectorValue(3.4,3.5)
test_basis = MockBasis{2}(v,ndofs)
test_basis_2 = MockBasis{2}(w,ndofs)

t1x = evaluate(test_basis,x)
t2x = evaluate(test_basis_2,x)
∇t1x = evaluate(∇(test_basis),x)
∇t2x = evaluate(∇(test_basis_2),x)

b = operate_fields(+,test_basis,test_basis_2)
r = broadcast(+,t1x,t2x)
∇r = broadcast(+,∇t1x,∇t2x)
test_field(b,x,r,grad=∇r)

b = operate_fields(-,test_basis,test_basis_2)
r = broadcast(-,t1x,t2x)
∇r = broadcast(-,∇t1x,∇t2x)
test_field(b,x,r,grad=∇r)

b = operate_fields(⋅,test_basis,test_basis_2)
r = broadcast(⋅,t1x,t2x)
∇r = broadcast(⋅,∇t1x,t2x) + broadcast(⋅,t1x,∇t2x)
test_field(b,x,r,grad=∇r)

b = operate_fields(*,test_basis,z)
r = broadcast(*,t1x,z)
∇r = broadcast(*,∇t1x,z)
test_field(b,x,r,grad=∇r)

b = operate_fields(+,test_basis,z)
test_field(b,x,fill(v+z,np,ndofs))

trial_basis = trialize_basis(test_basis)
r = reshape(evaluate(test_basis,x),(np,1,ndofs))
∇r = reshape(evaluate(∇(test_basis),x),(np,1,ndofs))
test_field(trial_basis,x,r,grad=∇r)

xl = Fill(x,l)
fl = [ z for  i in 1:l]

test_basis_array = Fill(test_basis,l)
test_basis_2_array = fill(test_basis_2,l)

bl = operate_arrays_of_fields(Nothing,*,test_basis_array,fl)
test_array_of_fields(bl,xl,fill(fill(z*v,np,ndofs),l))

bl = operate_arrays_of_fields(*,test_basis_array,fl)
test_array_of_fields(bl,xl,fill(fill(z*v,np,ndofs),l))

bl = operate_arrays_of_fields(+,test_basis_array,test_basis_2_array)
r = fill(broadcast(+,t1x,t2x),l)
∇r = fill(broadcast(+,∇t1x,∇t2x),l)
test_array_of_fields(bl,xl,r,grad=∇r)

bl = operate_arrays_of_fields(⋅,test_basis_array,test_basis_2_array)
r = fill(broadcast(⋅,t1x,t2x),l)
∇r = fill(broadcast(⋅,∇t1x,t2x) + broadcast(⋅,t1x,∇t2x),l)
test_array_of_fields(bl,xl,r,grad=∇r)

trial_basis_array = trialize_array_of_bases(operate_arrays_of_fields(*,test_basis_array,fl))

trial_basis_array_x = evaluate(trial_basis_array,xl)
@test trial_basis_array_x.g.value === trialize_basis_value

bl = operate_arrays_of_fields(⋅,trial_basis_array,test_basis_array)
r = fill(fill(z*v⋅v,np,ndofs,ndofs),l)
bl_x = evaluate(bl,xl)
test_array_of_fields(bl,xl,r)
@test bl_x.g.value == FieldOpKernel(⋅)

# Operations between values

al = [rand(np,ndofs) for k in 1:l]
bl = [rand(np) for k in 1:l]
cl = [rand(np,ndofs) for k in 1:l]

f(a,b) = 2*a-b*a
dl = apply(FieldOpKernel(f),bl,bl)
test_array(dl,map((a,b)->f.(a,b),bl,bl))

f(a,b) = 2*a-b*a
dl = apply(FieldOpKernel(f),al,bl)
test_array(dl,map((a,b)->f.(a,b),al,bl))

f(a,b) = 2*a-b
dl = apply(FieldOpKernel(f),al,cl)
test_array(dl,map((a,b)->f.(a,b),al,cl))

f(a,b,c) = b*(2*a-c)
dl = apply(FieldOpKernel(f),al,bl,cl)
test_array(dl,map((a,b,c)->f.(a,b,c),al,bl,cl))

dl = apply(trialize_basis_value,al)
test_array(dl,map(a->reshape(a,(size(a,1),1,size(a,2))),al))

atl = apply(trialize_basis_value,al)
ctl = apply(trialize_basis_value,cl)

f(a,b) = 2*a*b
dl = apply(FieldOpKernel(f),al,atl)
test_array(dl,map((a,b)->f.(a,b),al,atl))

f(a,c,at,ct) = 2*(a+c)*(2*at-ct)
dl = apply(FieldOpKernel(f),al,cl,atl,ctl)
test_array(dl,map((a,c,at,ct)->f.(a,c,at,ct),al,cl,atl,ctl))

@test size(dl[1]) == (np,ndofs,ndofs)

# Blocks

blocks = (al,)
blockids = [(1,1)]
axs_i = (blockedrange([np]),blockedrange([ndofs,ndofs]))
axs = Fill(axs_i,l)
aBl = VectorOfBlockArrayCoo(blocks,blockids,axs)

blocks = (cl,)
blockids = [(1,2)]
cBl = VectorOfBlockArrayCoo(blocks,blockids,axs)

atBl = apply(trialize_basis_value,aBl)
ctBl = apply(trialize_basis_value,cBl)
@test isa(atBl,VectorOfBlockArrayCoo)

f(a) = 2*a
dl = apply(FieldOpKernel(f),aBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a)->f.(a),aBl))

f(a,b) = 2*a + a*b
dl = apply(FieldOpKernel(f),aBl,bl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),aBl,bl))

f(b,a) = 2*a + a*b
dl = apply(FieldOpKernel(f),bl,atBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),bl,atBl))

f(b,a) = 2*a + b
dl = apply(FieldOpKernel(f),atBl,ctBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),atBl,ctBl))

dl = apply(FieldOpKernel(+),atBl,ctBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->a+b,atBl,ctBl))

dl = apply(FieldOpKernel(-),atBl,ctBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->a-b,atBl,ctBl))

f(b,a) = 2*a*b
dl = apply(FieldOpKernel(f),aBl,ctBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),aBl,ctBl))

f(b,a) = 2*a*b
dl = apply(FieldOpKernel(f),atBl,cBl)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),atBl,cBl))

# BlockWise integration

cj = [ fill(TensorValue(1,0,0,2),np) for cell in 1:l ]
cw = [ rand(np) for cell in 1:l ]

f(a,b) = 2*a + a*b
vl = apply(FieldOpKernel(f),aBl,bl)
dl = apply(IntKernel(),vl,cw,cj)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map( (v,w,j) -> reshape(sum( broadcast(*,v,w,det.(j)), dims=1),(2*ndofs,)), vl,cw,cj ))

vl = apply(FieldOpKernel(*),aBl,ctBl)
dl = apply(IntKernel(),vl,cw,cj)
@test isa(dl,VectorOfBlockArrayCoo)
test_array(dl,map( (v,w,j) -> reshape(sum( broadcast(*,v,w,det.(j)), dims=1),(2*ndofs,2*ndofs)), vl,cw,cj ))

# Blocks of Blocks

a0Bl = zeros_like(aBl)
c0Bl = zeros_like(cBl)
@test isa(a0Bl,VectorOfBlockArrayCoo)
@test blocksize(a0Bl) == blocksize(aBl)
@test length(a0Bl.blocks) == 0

blockids = [(1,1)]
ax1 = blockedrange([np])
ax2 = blockedrange([ndofs,ndofs])
axs_i  = (blockedrange([ax1]),blockedrange([ax2,ax2]))
axs = Fill(axs_i,l)
aLl = VectorOfBlockArrayCoo((aBl,),blockids,axs,(a0Bl,))
cLl = VectorOfBlockArrayCoo((cBl,),blockids,axs,(c0Bl,))

@test isa(aLl,VectorOfBlockArrayCoo)
@test aBl === aLl[Block(1,1)]
@test al === aLl[Block(1,1)][Block(1,1)]
@test isa(aLl[Block(1,1)],VectorOfBlockArrayCoo)
@test isa(aLl[Block(1,2)],VectorOfBlockArrayCoo)

blockids = [(1,2)]
aRl = VectorOfBlockArrayCoo((aBl,),blockids,axs,(a0Bl,))
cRl = VectorOfBlockArrayCoo((cBl,),blockids,axs,(c0Bl,))

a0tBl = zeros_like(atBl)
c0tBl = zeros_like(ctBl)

blockids = [(1,1,1)]
axs_i  = (blockedrange([ax1]),blockedrange([blockedrange([1])]),blockedrange([ax2,ax2]))
axs = Fill(axs_i,l)
atLl = VectorOfBlockArrayCoo((atBl,),blockids,axs,(a0tBl,))
ctLl = VectorOfBlockArrayCoo((ctBl,),blockids,axs,(c0tBl,))

@test isa(atLl,VectorOfBlockArrayCoo)
@test atBl === atLl[Block(1,1,1)]
@test atl === atLl[Block(1,1,1)][Block(1,1,1)]
@test isa(atLl[Block(1,1,1)],VectorOfBlockArrayCoo)
@test isa(atLl[Block(1,1,2)],VectorOfBlockArrayCoo)

blockids = [(1,1,2)]
atRl = VectorOfBlockArrayCoo((atBl,),blockids,axs,(a0tBl,))
ctRl = VectorOfBlockArrayCoo((ctBl,),blockids,axs,(c0tBl,))

f(a) = 2*a
dl = apply(FieldOpKernel(f),aLl)
@test isa(dl,VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,2)],BlockArrayCoo)
test_array(dl,map((a)->f.(a),aLl))

f(a,b) = 2*a + a*b
dl = apply(FieldOpKernel(f),aRl,bl)
@test isa(dl,VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,2)],BlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),aRl,bl))

f(b,a) = 2*a + b
dl = apply(FieldOpKernel(f),atRl,ctLl)
@test isa(dl,VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,1,2)],BlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),atRl,ctLl))

f(b,a) = 2*a + b
dl = apply(FieldOpKernel(f),aRl,cLl)
@test isa(dl,VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,2)],BlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),aRl,cLl))

f(b,a) = 2*a + b
dl = apply(FieldOpKernel(f),aRl,aRl)
@test isa(dl,VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,2)],BlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),aRl,aRl))

f(b,a) = a * b
dl = apply(FieldOpKernel(f),aLl,ctRl)
@test isa(dl,VectorOfBlockArrayCoo)
@test is_zero_block(dl,Block(1,1,1))
@test is_nonzero_block(dl,Block(1,1,2))
@test is_zero_block(dl,Block(1,2,1))
@test is_zero_block(dl,Block(1,2,2))
@test isa(dl[Block(1,1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1,2)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,1,2)],BlockArrayCoo)
@test isa(dl[1][Block(1,2,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,2,2)],BlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),aLl,ctRl))

f(b,a) = a * b
v = apply(FieldOpKernel(-),ctLl,ctRl)
dl = apply(FieldOpKernel(f),v,aLl)
@test isa(dl,VectorOfBlockArrayCoo)
@test is_nonzero_block(dl,Block(1,1,1))
@test is_nonzero_block(dl,Block(1,1,2))
@test is_zero_block(dl,Block(1,2,1))
@test is_zero_block(dl,Block(1,2,2))
@test isa(dl[Block(1,1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,1,2)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2,2)],VectorOfBlockArrayCoo)
@test isa(dl[1][Block(1,1,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,1,2)],BlockArrayCoo)
@test isa(dl[1][Block(1,2,1)],BlockArrayCoo)
@test isa(dl[1][Block(1,2,2)],BlockArrayCoo)
test_array(dl,map((a,b)->f.(a,b),v,aLl))

# Integration of Blocks of Blocks

f(a,b) = 2*a + a*b
vl = apply(FieldOpKernel(f),aRl,bl)
dl = apply(IntKernel(),vl,cw,cj)
test_array(dl,map( (v,w,j) -> reshape(sum( broadcast(*,v,w,det.(j)), dims=1),(4*ndofs,)), vl,cw,cj ))
@test isa(dl,VectorOfBlockArrayCoo)
@test is_zero_block(dl,Block(1))
@test is_nonzero_block(dl,Block(2))
@test isa(dl[Block(1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(2)],VectorOfBlockArrayCoo)

f(b,a) = a * b
v = apply(FieldOpKernel(-),ctLl,ctRl)
vl = apply(FieldOpKernel(f),v,aLl)
dl = apply(IntKernel(),vl,cw,cj)
test_array(dl,map( (v,w,j) -> reshape(sum( broadcast(*,v,w,det.(j)), dims=1),(4*ndofs,4*ndofs)), vl,cw,cj ))
@test isa(dl,VectorOfBlockArrayCoo)
@test is_nonzero_block(dl,Block(1,1))
@test is_nonzero_block(dl,Block(1,2))
@test is_zero_block(dl,Block(2,1))
@test is_zero_block(dl,Block(2,2))
@test isa(dl[Block(1,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(1,2)],VectorOfBlockArrayCoo)
@test isa(dl[Block(2,1)],VectorOfBlockArrayCoo)
@test isa(dl[Block(2,2)],VectorOfBlockArrayCoo)

#using BenchmarkTools
#cache = array_cache(dl)
#@btime getindex!($cache,$dl,3)

end # module

ext: jl


text:
#####
##### Divergence operators
#####

"""
    hdivᶜᶜᵃ(i, j, k, grid, u, v)

Calculates the horizontal divergence ∇ₕ·(u, v) of a 2D velocity field (u, v) via

    1/V * [δxᶜᵃᵃ(Ax * u) + δyᵃᶜᵃ(Ay * v)]

which will end up at the location `cca`.
"""
@inline function hdivᶜᶜᵃ(i, j, k, grid, u, v)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v))
end

"""
    divᶜᶜᶜ(i, j, k, grid, u, v, w)

Calculates the divergence ∇·U of a vector field U = (u, v, w),

    1/V * [δxᶜᵃᵃ(Ax * u) + δxᵃᶜᵃ(Ay * v) + δzᵃᵃᶜ(Az * w)],

which will end up at the cell centers `ccc`.
"""
@inline function divᶜᶜᶜ(i, j, k, grid, u, v, w)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, w))
end


text:
struct AddNewCard
    value::Observable{String}
    list::List
end

function AddNewCard(keys::Observable{Vector{String}}, value=Observable(""), list::List=List(keys, value))
    return AddNewCard(value, list)
end

function jsrender(session::Session, add::AddNewCard)
    list = JSServe.jsrender(session, add.list)
    hidden, keydown = add.list.hidden, add.list.keydown
    onfocus = js"JSServe.update_obs($(hidden), false);"
    onblur=js"""
        const tgt = event.relatedTarget;
        tgt && $(list).contains(tgt) || JSServe.update_obs($(hidden), true);
    """
    onkeydown = js"event.key == 'Escape' ? this.blur() : JSServe.update_obs($(keydown), event.key)"
    box = DOM.button(
        "+";
        class="w-full p-8 cursor-pointer text-left text-blue-800 text-2xl hover:bg-gray-200 hover:text-blue-900",
        onfocus,
        onblur,
        onkeydown
    )
    onjs(session, add.value, js"function (value) {JSServe.update_obs($(hidden), true);}")
    ui = DOM.div(
        box,
        list,
        dataType="add",
        dataId=string(objectid(add)),
        dataSelected="false",
    )
    return jsrender(session, ui)
end

struct MoveSelected end

struct EditableList
    options::Observable{SimpleList}
    steps::Observable{SimpleList}
    selected::Observable{Vector{String}}
    list::Observable{SimpleList}
end

function get_step_indices(steps::Vector, selected::Vector{String})
    ids = @. string(objectid(getproperty(steps, :card)))
    return filter(!isnothing, indexin(selected, ids))
end

updatecards(el::EditableList, idx, card) = insert_item(el.steps[], idx, card)

function updatecards(el::EditableList, idx, ::MoveSelected)
    _steps, _selected = el.steps[], el.selected[]
    selected_idx = @maybereturn findfirst(_steps) do step
        return string(objectid(step.card)) in _selected
    end
    return move_item(_steps, selected_idx => idx - (idx > selected_idx))
end

function AddNewCard(keys::Observable{Vector{String}}, el::EditableList)
    add = AddNewCard(keys)
    on(add.value) do key
        thunk = @maybereturn getatkey(el.options[], key)
        idx = @maybereturn indexoftype(AddNewCard, el.list[], add)
        steps = @maybereturn updatecards(el, idx, thunk())
        el.steps[] = steps
    end
    return add
end

function EditableList(options::Observable, steps::Observable)
    selected = Observable(String[])
    keys = @lift getkey.($options)
    list = Observable{SimpleList}()
    el = EditableList(options, steps, selected, list)
    map!(list, steps) do steps
        elements = Any[]
        push!(elements, AddNewCard(keys, el))
        for step in steps
            push!(elements, step)
            push!(elements, AddNewCard(keys, el))
        end
        push!(elements, DOM.div(class="h-64"))
        return elements
    end
    return el
end

function jsrender(session::Session, el::EditableList)
    ui = map(DOM.div, session, el.list)
    return jsrender(session, ui)
end


text:
"Initialize LIKWIDs affinity domain module."
function init_affinity()
    LibLikwid.affinity_init()
    _affinity[] = unsafe_load(LibLikwid.get_affinityDomains())
    _build_jl_affinity()
    affinity_initialized[] = true
    return true
end

"Close and finalize LIKWIDs affinity domain module."
function finalize_affinity()
    LibLikwid.affinity_finalize()
    affinity_initialized[] = false
    _affinity[] = nothing
    return nothing
end

function _build_jl_affinity()
    af = _affinity[]

    ndomains = af.numberOfAffinityDomains
    _domains = unsafe_wrap(Array, af.domains, ndomains)
    domains = Vector{AffinityDomain}(undef, ndomains)
    for (i, d) in enumerate(_domains)
        # bstring: http://mike.steinert.ca/bstring/doc/
        bstr = unsafe_load(d.tag)
        domains[i] = AffinityDomain(
            unsafe_string(bstr.data),
            d.numberOfProcessors,
            d.numberOfCores,
            unsafe_wrap(Array, d.processorList, d.numberOfProcessors),
        )
    end

    affinity[] = AffinityDomains(
        af.numberOfSocketDomains,
        af.numberOfNumaDomains,
        af.numberOfProcessorsPerSocket,
        af.numberOfCacheDomains,
        af.numberOfCoresPerCache,
        af.numberOfProcessorsPerCache,
        ndomains,
        domains
    )
    return nothing
end

"Query affinity domain information"
function get_affinity()
    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    if !numa_initialized[]
        init_numa() || error("Couldn't init numa.")
    end
    if !affinity_initialized[]
        init_affinity() || error("Couldn't init affinity.")
    end
    return affinity[]
end

"Transform a valid cpu string in LIKWID syntax into a list of CPU IDs"
function cpustr_to_cpulist(cpustr::AbstractString)
    if !config_initialized[]
        init_configuration() || error("Couldn't init configuration.")
    end
    config = get_configuration()
    cpulist = zeros(Int32, config.maxNumThreads)
    LibLikwid.cpustr_to_cpulist(cpustr, cpulist, length(cpulist))
    return cpulist
end

"""
Returns the ID of the currently executing CPU.
"""
get_processor_id() = Int(LibLikwid.likwid_getProcessorId())

"""
Returns the ID of the currently executing CPU via `glibc`s `sched_getcpu` function.
"""
get_processor_id_glibc() = Int(@ccall sched_getcpu()::Cint)

"""
Pins the current process to the CPU given as `cpuid`.
"""
function pinprocess(cpuid::Integer)
    ret = LibLikwid.likwid_pinProcess(cpuid)
    return ret == 0
end

"""
Pins the current thread to the CPU given as `cpuid`.
"""
function pinthread(cpuid::Integer)
    ret = LibLikwid.likwid_pinThread(cpuid)
    return ret == 0
end

"""
Get the CPU core IDs of the Julia threads.
"""
function get_processor_ids()
    N = Threads.nthreads()
    coreids = zeros(Int, N)
    Threads.@threads :static for i in 1:N
        coreids[i] = LIKWID.get_processor_id()
    end
    return coreids
end

"""
Pin all Julia threads to the CPU cores `coreids`.
Note that `length(coreids) == Threads.nthreads()` must hold!
"""
function pinthreads(coreids::AbstractVector{<:Integer})
    N = Threads.nthreads()
    @assert length(coreids) == N
    @assert minimum(coreids) ≥ 0
    Threads.@threads :static for i in 1:N
        LIKWID.pinthread(coreids[i])
    end
    return nothing
end


text:
#!/usr/local/bin/julia
# -*- coding: utf-8 -*-

# $ julia --precompiled=yes test_tk.jl
# julia> include("test/test_tk.jl")

using Tk

widgets = Dict{Symbol, Union{Void, Tk.Widget}}()

function univPath(p)
  return replace(p, "\\", "/")
end

function test_tk(p::AbstractString)
  widgets[:w] = w = Toplevel(p, 640, 480)
  tcl("pack", "propagate", w, false) # or pack_stop_propagate(w)
  bind(w, "<Destroy>", path -> widgets[:w] = nothing)

  cbk_file = function (path)
    println(path)
    d = ChooseDirectory()
    length(d) == 0 && (println("skipped: cd"); return 0)
    println(d)
    cd(d)
    set_value(widgets[:lbl], d)
    # println(get_items(widgets[:rb].buttons[2]))
    # set_items(widgets[:rb].buttons[2], d)
    o = GetOpenFile() # can not set directory (needs cd() ?)
    length(o) == 0 && (println("canceled: open"); return 0)
    println(o)
    s = GetSaveFile() # can not set directory (needs cd() ?)
    length(s) == 0 && (println("canceled: save"); return 0)
    println(s)
  end

  mbar = Menu(w)
  fmenu = menu_add(mbar, "(&F)ile")
  omenu = menu_add(mbar, "(&O)ptions")
  hmenu = menu_add(mbar, "(&H)elp")
  menu_add(fmenu, "開く(&O)...", path -> cbk_file(path))
  menu_add(fmenu, Separator(w))
  menu_add(fmenu, "終了(&X)", path -> destroy(w))
  cb = Checkbutton(w, "チェックボックス")
  set_value(cb, true)
  menu_add(omenu, cb)
  menu_add(omenu, Separator(w))
  # menu_add(omenu, Label(w, "ラベル"))
  widgets[:rb] = rb = Radio(w, ["option 1", "オプション 2"])
  set_value(rb, "option 1")
  menu_add(omenu, rb)
  menu_add(hmenu, "(&A)bout", path -> cbk_dsp(path))

  nb = Notebook(w)
  pack(nb, expand=true, fill="both")

  img = Image(univPath(Pkg.dir("Tk", "examples", "logo.gif")))
  f1 = Frame(nb, padding=[16,16,8,8], relief="groove")
  page_add(f1, "Tab 1")
  widgets[:lbl] = lbl = Label(f1, "label ラベル")
  grid(lbl, 1, 1, sticky="sw")
  b = Button(f1, "表示 selected options", img)
  grid(b, 2, 1, sticky="news")
  cbk_dsp(path) = println(@sprintf "%s: %s" path map(get_value, (cb, rb)))
  # callback_add(b, cbk_dsp)
  bind(b, "command", cbk_dsp)
  bind(b, "<Return>", cbk_dsp)
  sc = Slider(f1, 1:20)
  v = Label(f1, "value 選択ディレクトリ")
  v[:textvariable] = sc[:variable]
  grid(v, 3, 1, sticky="ne")
  grid(sc, 4, 1, sticky="ew")
  grid_columnconfigure(f1, 1, weight=1)
  grid_rowconfigure(f1, 1, weight=3)
  grid_rowconfigure(f1, 3, weight=1)

  sash_dir = "vertical" # "horizontal"
  f2 = Frame(nb, padding=[8,16,8,16], relief="groove")
  page_add(f2, "椨 2")
  p21 = Panedwindow(f2, sash_dir)
  # pack(p21, expand=true, fill="both")
  grid(p21, 1, 1, sticky="news")
  grid_columnconfigure(f2, 1, weight=1)
  grid_rowconfigure(f2, 1, weight=1)
    pp1 = Frame(p21, padding=[4,4,4,4], relief="sunken")
    page_add(pp1, 2)
      b2 = Button(pp1, "鹽🐢 2")
      pack(b2, expand=true, fill="both")
      # callback_add(b2, cbk_dsp)
      bind(b2, "command", cbk_dsp)
      bind(b2, "<Return>", cbk_dsp)
    pp2 = Frame(p21, padding=[4,4,4,4], relief="sunken")
    page_add(pp2, 1)
      l2 = Label(pp2, "paned2ラベル")
      pack(l2, expand=true, fill="both")
    pp3 = Frame(p21, padding=[4,4,4,4], relief="sunken")
    page_add(pp3, 3)
      ff = Frame(pp3)
      formlayout(Entry(ff), "Name:")
      formlayout(Entry(ff), "Rank:")
      formlayout(Entry(ff), "Serial Number:")
      pack(ff, expand=true, fill="both")
    pp4 = Frame(p21, padding=[4,4,4,4], relief="sunken")
    page_add(pp4, 2)
      l3 = Label(pp4, "paned3ラベル")
      pack(l3, expand=true, fill="both")
  set_value(p21, 100)
  tcl(p21, "sashpos", 1, 200)
  tcl(p21, "sashpos", 2, 300)

  # load gif sash image http://wiki.tcl.tk/20058
  img_sash = (
    ("horizontal", "xsash", "ew", """
R0lGODlhBQB5AKIAAP///9LS0s7OzjU1NTExMf4BAgAAAAAAACH5BAUUAAUALAAAAAAFAHkAAANZ
KLoztDDKucJQ9ommc71UGAHBR5HmhGKgtArWWbJy6s5wO+KxyuuQV+9m2xWDvxrNd1wIgc4kc0mk
Gq3IpuKpzHWHV+9UXCWHwVlsVCuS9DwbB0tjbts5kgQAOw==
"""),
    ("vertical", "ysash", "ns", """
R0lGODlheQAFAKIAAP///9LS0s7OzjU1NTExMf4BAgAAAAAAACH5BAUUAAUALAAAAAB5AAUAAANO
KLrc/jDKSau9agSnsf9gKHXNkG3MoKJKAAgu/Mb0bMt4nd96z/+7oC61KikILAFhSSAJOtBN9Cmt
Uq/TrFWL3Xq7YKeSKSaIzui0+pEAADs=
"""),
  )
  for (d, hv, news, img) in img_sash
    d != sash_dir && continue
    try
      tcl_eval("image create photo img:$hv -data {$img}") # => "img:?sash"
      tcl_eval("""ttk::style element create Sash.$hv image [list img:$hv] \\
        -border {1 1} -sticky $news -padding {1 1}""") # => "Sash.?sash"
    catch e
      if ! isa(e, Tk.TclError) throw(e) end
    finaly
      tcl_eval("ttk::style layout TPanedwindow { Sash.$hv }") # => ""
    end
  end

  f3 = Frame(nb, padding=[16,16,8,8], relief="groove")
  page_add(f3, "椨 3")

  set_value(nb, 2)

  # Messagebox(title=p, message="Hello, 世界!") # to keep

  exists(w) && cbk_dsp(w) # skip after destroy(w)
  # widgets[:w] != nothing && cbk_dsp(w) # skip after destroy(w)

  return w
end

if ! isinteractive()
  w = test_tk("TEST日本語漢字表示申能 A")
  Messagebox(title=get_value(w), message="Hello, 世界! A") # to keep
end
if current_module() == Main
  w = test_tk("TEST日本語漢字表示申能 B")
  Messagebox(title=get_value(w), message="Hello, 世界! B") # to keep
end


text:
using OrdinaryDiffEq, ParameterizedFunctions, ODE, ODEInterface,
      ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools

f = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

p = [0.7,0.8,1/12.5,0.5]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10);

sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()

plot(sol)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) #fails
          Dict(:alg=>dopri5())
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) # fails
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=10000,error_estimate=:L2,dense_errors=true)
plot(wp)

setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78()) # fails
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>DP8())
          Dict(:alg=>Vern7())
          Dict(:alg=>CVODE_Adams())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
          Dict(:alg=>lsoda())
          Dict(:alg=>odex())
          Dict(:alg=>ddeabm())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78()) # fails
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=1000,error_estimate=:L2,dense_errors=true)
plot(wp)

setups = [Dict(:alg=>Tsit5())
          Dict(:alg=>Vern9())
          Dict(:alg=>VCABM())
          Dict(:alg=>AitkenNeville(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))]
solnames = ["Tsit5","Vern9","VCABM","AitkenNeville","Midpoint Deuflhard","Midpoint Hairer Wanner"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=9, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :romberg, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :bulirsch, threading=true))]
solnames = ["Deuflhard","No threads","standard","Romberg","Bulirsch"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=15, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=7, init_order=6, threading=true))]
solnames = ["1","2","3","4","5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
"""
    initialize_hbv_model(config::Config)

Initial part of the SBM model concept. Reads the input settings and data as defined in the
Config object. Will return a Model that is ready to run.
"""
function initialize_hbv_model(config::Config)
    # unpack the paths to the NetCDF files
    tomldir = dirname(config)
    static_path = joinpath(tomldir, config.input.path_static)
    cyclic_path = joinpath(tomldir, config.input.path_static)
    dynamic_path = joinpath(tomldir, config.input.path_forcing)
    instate_path = joinpath(tomldir, config.state.path_input)
    output_path = joinpath(tomldir, config.output.path)

    Δt = Second(config.timestepsecs)

    sizeinmetres = get(config.model, "sizeinmetres", false)
    reinit = get(config.model, "reinit", true)
    do_reservoirs = get(config.model, "reservoirs", false)
    do_lakes = get(config.model, "lakes", false)
    do_pits = get(config.model, "pits", false)
    set_kquickflow = get(config.model, "set_kquickflow", false)

    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
    kw_land_tstep = get(config.model, "kw_land_tstep", 0)

    nc = NCDataset(static_path)
    dims = dimnames(nc[param(config, "input.subcatchment")])

    # There is no need to permute the dimensions of the data, since the active indices are
    # correctly calculated in both ways.
    # The dimension order only needs to be known for interpreting the LDD directions
    # and creating the coordinate maps.
    dims_xy = dims[2] in ("y", "lat")
    subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
    # indices based on catchment
    inds, rev_inds = active_indices(subcatch_2d, missing)
    n = length(inds)

    cfmax = ncread(
        nc,
        param(config, "input.vertical.cfmax", nothing);
        sel = inds,
        defaults = 3.75653,
        type = Float64,
    ) .* (Δt / basetimestep)
    tt = ncread(
        nc,
        param(config, "input.vertical.tt", nothing);
        sel = inds,
        defaults = -1.41934,
        type = Float64,
    )
    tti = ncread(
        nc,
        param(config, "input.vertical.tti", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )
    ttm = ncread(
        nc,
        param(config, "input.vertical.ttm", nothing);
        sel = inds,
        defaults = -1.41934,
        type = Float64,
    )
    whc = ncread(
        nc,
        param(config, "input.vertical.whc", nothing);
        sel = inds,
        defaults = 0.1,
        type = Float64,
    )
    # glacier parameters
    g_tt = ncread(
        nc,
        param(config, "input.vertical.g_tt", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float64,
        fill = 0.0,
    )
    g_cfmax = ncread(
        nc,
        param(config, "input.vertical.g_cfmax", nothing);
        sel = inds,
        defaults = 3.0,
        type = Float64,
        fill = 0.0,
    ).* (Δt / basetimestep)
    g_sifrac = ncread(
        nc,
        param(config, "input.vertical.g_sifrac", nothing);
        sel = inds,
        defaults = 0.001,
        type = Float64,
        fill = 0.0,
    )
    glacierfrac = ncread(
        nc,
        param(config, "input.vertical.glacierfrac", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float64,
        fill = 0.0,
    )
    glacierstore = ncread(
        nc,
        param(config, "input.vertical.glacierstore", nothing);
        sel = inds,
        defaults = 5500.0,
        type = Float64,
        fill = 0.0,
    )
    fc = ncread(
        nc,
        param(config, "input.vertical.fc", nothing);
        sel = inds,
        defaults = 260.0,
        type = Float64,
    )
    betaseepage = ncread(
        nc,
        param(config, "input.vertical.betaseepage", nothing);
        sel = inds,
        defaults = 1.8,
        type = Float64,
    )
    lp = ncread(
        nc,
        param(config, "input.vertical.lp", nothing);
        sel = inds,
        defaults = 0.53,
        type = Float64,
    )
    k4 = ncread(
        nc,
        param(config, "input.vertical.k4", nothing);
        sel = inds,
        defaults = 0.02307,
        type = Float64,
    ) .* (Δt / basetimestep)
    kquickflow = ncread(
        nc,
        param(config, "input.vertical.kquickflow", nothing);
        sel = inds,
        defaults = 0.09880,
        type = Float64,
    ) .* (Δt / basetimestep)
    suz = ncread(
        nc,
        param(config, "input.vertical.suz", nothing);
        sel = inds,
        defaults = 100.0,
        type = Float64,
    )
    k0 = ncread(
        nc,
        param(config, "input.vertical.k0", nothing);
        sel = inds,
        defaults = 0.30,
        type = Float64,
    ) .* (Δt / basetimestep)
    khq = ncread(
        nc,
        param(config, "input.vertical.khq", nothing);
        sel = inds,
        defaults = 0.09880,
        type = Float64,
    ) .* (Δt / basetimestep)
    hq = ncread(
        nc,
        param(config, "input.vertical.hq", nothing);
        sel = inds,
        defaults = 3.27,
        type = Float64,
    ) .* (Δt / basetimestep)
    alphanl = ncread(
        nc,
        param(config, "input.vertical.alphanl", nothing);
        sel = inds,
        defaults = 1.1,
        type = Float64,
    )
    perc = ncread(
        nc,
        param(config, "input.vertical.perc", nothing);
        sel = inds,
        defaults = 0.4,
        type = Float64,
    ) .* (Δt / basetimestep)
    cfr = ncread(
        nc,
        param(config, "input.vertical.cfr", nothing);
        sel = inds,
        defaults = 0.05,
        type = Float64,
    )
    pcorr = ncread(
        nc,
        param(config, "input.vertical.pcorr", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )
    rfcf = ncread(
        nc,
        param(config, "input.vertical.rfcf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )
    sfcf = ncread(
        nc,
        param(config, "input.vertical.sfcf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )
    cflux = ncread(
        nc,
        param(config, "input.vertical.cflux", nothing);
        sel = inds,
        defaults = 2.0,
        type = Float64,
    ) .* (Δt / basetimestep)
    icf = ncread(
        nc,
        param(config, "input.vertical.icf", nothing);
        sel = inds,
        defaults = 2.0,
        type = Float64,
    )
    cevpf = ncread(
        nc,
        param(config, "input.vertical.cevpf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )
    epf = ncread(
        nc,
        param(config, "input.vertical.epf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )
    ecorr = ncread(
        nc,
        param(config, "input.vertical.ecorr", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float64,
    )

    # read x, y coordinates and calculate cell length [m]
    y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
    x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
    if dims_xy
        y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
    else
        y = repeat(y_nc, outer = (1, length(x_nc)))[inds]
    end
    cellength = abs(mean(diff(x_nc)))

    xl = fill(mv, n)
    yl = fill(mv, n)
    for i=1:n
        xl[i] = sizeinmetres ? cellength : lattometres(y[i])[1] * cellength
        yl[i] = sizeinmetres ? cellength : lattometres(y[i])[2] * cellength
    end

    threshold = fc .* lp

    altitude =
    ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float64)

    hbv = HBV{Float64}(
        n = n,
        yl = yl,
        xl = xl,
        altitude = altitude,
        fc = fc,
        betaseepage = betaseepage,
        lp = lp,
        threshold = threshold,
        k4 = k4,
        kquickflow = set_kquickflow ? kquickflow : pow.(khq, 1.0 .+ alphanl) .* pow.(hq, -alphanl),
        suz = suz,
        k0 = k0,
        khq = khq,
        hq = hq,
        alphanl = alphanl,
        perc = perc,
        cfr = cfr,
        pcorr = pcorr,
        rfcf = rfcf,
        sfcf = sfcf,
        cflux = cflux,
        icf = icf,
        cevpf = cevpf,
        epf = epf,
        ecorr = ecorr,
        tti = tti,
        tt = tt,
        ttm = ttm,
        cfmax = cfmax,
        whc = whc,
        # glacier parameters
        g_tt = g_tt,
        g_sifrac = g_sifrac,
        g_cfmax = g_cfmax,
        glacierstore = glacierstore,
        glacierfrac = glacierfrac,
        # default (cold) states:
        interceptionstorage = fill(0.0, n),
        snow = fill(0.0, n),
        snowwater = fill(0.0, n),
        soilmoisture = copy(fc),
        upperzonestorage = 0.2 .* fc,
        lowerzonestorage = 1.0 ./ (3.0 .* k4),
    )

    modelsize_2d = size(subcatch_2d)
    river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
    river = river_2d[inds]
    riverwidth_2d =
        ncread(nc, param(config, "input.lateral.river.width"); type = Float64, fill = 0)
    riverwidth = riverwidth_2d[inds]
    riverlength_2d =
        ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
    riverlength = riverlength_2d[inds]

    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
    nriv = length(inds_riv)

     # reservoirs
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoirs, resindex, reservoir, pits = initialize_simple_reservoir(config, nc, inds_riv, nriv, pits)
     else
         reservoir = ()
     end

     # lakes
     if do_lakes
         lakes, lakeindex, lake, pits = initialize_natural_lake(config, nc, inds_riv, nriv, pits)
     else
         lake = ()
     end

     ldd_2d = ncread(nc, param(config, "input.ldd"); allow_missing = true)
     ldd = ldd_2d[inds]
     if do_pits
        pits_2d = ncread(nc, param(config, "input.pits"); type = Bool, fill = false)
        ldd = set_pit_ldd(pits_2d, ldd, inds)
    end

    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float64)
    clamp!(βₗ, 0.00001, Inf)
    dl = fill(mv, n)
    dw = fill(mv, n)

    for i = 1:n
        dl[i] = detdrainlength(ldd[i], xl[i], yl[i])
        dw[i] = detdrainwidth(ldd[i], xl[i], yl[i])
    end

    n_land = ncread(
        nc,
        param(config, "input.lateral.land.n", nothing);
        sel = inds,
        defaults = 0.072,
        type = Float64,
    )

    olf = SurfaceFlow(
        sl = βₗ,
        n = n_land,
        dl = dl,
        Δt = tosecond(Δt),
        its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
        width = dw,
        wb_pit = pits[inds],
    )

    pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
    graph = flowgraph(ldd, inds, pcr_dir)

    riverslope = ncread(
        nc,
        param(config, "input.lateral.river.slope");
        sel = inds_riv,
        type = Float64,
    )
    clamp!(riverslope, 0.00001, Inf)
    riverlength = riverlength_2d[inds_riv]
    riverwidth = riverwidth_2d[inds_riv]
    n_river = ncread(
        nc,
        param(config, "input.lateral.river.n", nothing);
        sel = inds_riv,
        defaults = 0.036,
        type = Float64,
    )
    ldd_riv = ldd_2d[inds_riv]
    if do_pits
        ldd_riv = set_pit_ldd(pits_2d, ldd_riv, inds_riv)
    end
    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)

    # the indices of the river cells in the land(+river) cell vector
    index_river = filter(i -> !isequal(river[i], 0), 1:n)
    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)

    rf = SurfaceFlow(
        sl = riverslope,
        n = n_river,
        dl = riverlength,
        Δt = tosecond(Δt),
        its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) : kw_river_tstep,
        width = riverwidth,
        reservoir_index = do_reservoirs ? resindex : fill(0, nriv),
        lake_index = do_lakes ? lakeindex : fill(0, nriv),
        reservoir = do_reservoirs ? reservoirs : nothing,
        lake = do_lakes ? lakes : nothing,
        rivercells = river,
    )

    state_ncnames = ncnames(config.state)

    reader = prepare_reader(dynamic_path, cyclic_path, config)

    modelmap = (vertical = hbv, lateral = (land = olf, river = rf))
    indices_reverse = (
        land = rev_inds,
        river = rev_inds_riv,
        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
        lake = isempty(lake) ? nothing : lake.reverse_indices,
    )
    writer = prepare_writer(
        config,
        reader,
        output_path,
        modelmap,
        state_ncnames,
        indices_reverse,
        x_nc,
        y_nc,
        dims_xy,
        nc,
    )

    # for each domain save the directed acyclic graph, the traversion order,
    # and the indices that map it back to the two dimensional grid
    land = (
        graph = graph,
        order = topological_sort_by_dfs(graph),
        indices = inds,
        reverse_indices = rev_inds,
    )
    river =
        (graph = graph_riv, order = topological_sort_by_dfs(graph_riv), indices = inds_riv, reverse_indices = rev_inds_riv)

    model = Model(
        config,
        (; land, river, reservoir, lake, index_river, frac_toriver),
        (land = olf, river = rf),
        hbv,
        Clock(config.starttime, 1, Δt),
        reader,
        writer,
    )

    # read and set states in model object if reinit=true
    if reinit == false
        set_states(
            instate_path,
            model,
            state_ncnames;
            type = Float64,
        )
    end

    # make sure the forcing is already loaded
    # it's fine to run twice, and may help catching errors earlier
    update_forcing!(model)
    if haskey(config.input, "cyclic")
        update_cyclic!(model)
    end
    return model
end

function update(model::Model{N,L,V,R,W}) where {N,L,V<:HBV,R,W}
    @unpack lateral, vertical, network, clock, config = model

    kinwave_it = get(config.model, "kin_wave_iteration", false)

    update_forcing!(model)
    if haskey(config.input, "cyclic")
        update_cyclic!(model)
    end

    update_until_snow(vertical, config)

    if get(config.model, "masswasting", false)
        snowflux_frac =
            min.(0.5, lateral.land.sl ./ 5.67) .* min.(1.0, vertical.snow ./ 10000.0)
        maxflux = snowflux_frac .* vertical.snow
        vertical.snow .= accucapacityflux(network.land, vertical.snow, maxflux)
        vertical.snowwater .= accucapacityflux(network.land, vertical.snowwater, vertical.snowwater .* snowflux_frac)
    end

    update_after_snow(vertical, config)

    lateral.land.qlat .=
        (vertical.runoff .* vertical.xl .* vertical.yl .* 0.001) ./ lateral.land.Δt ./
        lateral.land.dl

    update(
        lateral.land,
        network.land,
        frac_toriver = network.frac_toriver,
        river = lateral.river.rivercells,
        do_iter = kinwave_it,
    )

    lateral.river.qlat .=
        (lateral.land.to_river[network.index_river]) ./ lateral.river.dl

    update(lateral.river, network.river, do_iter = kinwave_it, doy = dayofyear(clock.time))

    write_output(model)

    # update the clock
    advance!(clock)

    return model
end


text:
using ArenaSim
using Statistics

include("utils.jl")

files = readdir("decks/top8")

decks = ArenaSim.deckreader_mtgo_format.( joinpath.("decks/top8", files),
	strip_sideboard = false)

decks = deckinfo.(decks)

pars = SimParameters(
	icrs_per_pack = ArenaSim.icrgen_qc(0.5,0,0),
	welcome_bundle = true,
	bonus_packs = Dict{Symbol,Int}(:M19 => 5, :GRN => 4)
	)

sets = ArenaSim.sets[1:5]

perm = [7, 2, 5, 8, 6, 3, 4, 1]
decks = decks[perm]
files = files[perm]

function f(i)
    res = simulate(10000, decks[i:i]; parameters =pars, sets = sets)
    res
end

zz = map(i -> f(i), 1:length(decks))

using PyCall
pygui()
using PyPlot

using Colors

pal = distinguishable_colors(9)
pal = "#" .* lowercase.(hex.(pal))

tot = ArenaSim.reduce_deck(reduce(vcat, decks); additive = false)
res2 = simulate(10000, [tot]; parameters = pars, sets = sets)

PyPlot.plt[:close]("all")
xmax = extrema( extrema(x.total_earned_packs) for x in zz)[2][2]
fig = figure(figsize = (11,8.5))
for i in 1:length(zz)
    subplot(3,3,i)
    PyPlot.plt[:hist]( zz[i].total_earned_packs, bins = -0.5:5:175.5, density = true, color = pal[i])
    title(files[i] * "\n" * string(rm_count(decks[i])), fontsize = 10)
    xticks(0:25:175)
end
subplot(3,3,9)
xr = extrema(res2.total_earned_packs)
PyPlot.plt[:hist](res2.total_earned_packs, color = pal[9], bins = 220:5:400, density = true)
title("All 8 lists crafted in bulk" * "\n" * string(rm_count(tot)), fontsize = 10)
suptitle("Pack Cost Estimates for Top 8 GP Lille Decklists in MTGA (assuming fresh account, welcome bundle purchased)")
figtext(0.5, 0.02, "Packs needed after welcome bundle and starting packs", horizontalalignment = "center", fontsize = "large")
figtext(0.015,  0.5, "Density",
rotation = "vertical", horizontalalignment = "right", verticalalignment="center",
fontsize = "large")
tight_layout(rect=[0.03, 0.03, 1, 0.95])
savefig("earned_pack_estimates.png")

max_stars_repo_name: asieev/GP-Lille

max_issues_repo_name: asieev/GP-Lille

max_forks_repo_name: asieev/GP-Lille


text:
using Revise

using DataFrames, Dates, Random

#=function mwe(;N)
  alldates = collect(Date(1950,1,1):Day(1):Date(2020,12,31))
  df = reduce(vcat, [DataFrame(;grp = i, date=alldates[1:(rand(1:length(alldates)))]) for i ∈ 1:N])
  df.yearmonth = (year).(df.date) .+ (month).(df.date)
  df.val = rand([1,2,3,missing],nrow(df))
  sort!(df, [:grp, :yearmonth])
  dfeom = combine(groupby(df, [:grp, :yearmonth]), propertynames(df) .=> last, nrow => :ndays)

  dfeom.tokeep = trues(nrow(dfeom))
  dfbyeom = groupby(dfeom, :yearmonth)
  Threads.@threads for i ∈ length(dfbyeom)
    sdf = dfbyeom[i]
    mindays = floor(median(sdf.ndays)) |> Int
    sdf.tokeep .= (sdf.ndays .≥ mindays) .& (mindays ≥ 2)
  end

  dfeom = dfeom[dfeom.tokeep, :]

  #=print("size(df): $(size(df))")
  for s ∈ propertynames(df)
    (!(eltype(df[!,s]) <: Union{Real, Missing})) && continue
    print(" $s $(mean(skipmissing(df[!, s]))),")
  end
  print("\n")=#

  return dfeom
end

function testmwe(;loops, N)
  Random.seed!(1111)
  base = mwe(;N)
  for i ∈ 1:loops
    Random.seed!(1111)
    iter = mwe(;N)
    for s ∈ propertynames(base)
      if !all( iter[!,s] .=== base[!,s])
        throw("Column mismatch! i=$i, col=$s")
      end
    end
  end
end

@time testmwe(loops=10, N=1000)=#
#=function mweold()
  df = DataFrame(G=[:A,:B,:C,:A,:B,:C], X = rand(6))
  local combinedlast

  println(df)

  try
    combineddf = combine(last, groupby(df, :G))
  catch err
    @warn "(:) .=> last failed. Error: $err"
  end

  #workaround
  combineddf = combine(groupby(df, :G), setdiff!(propertynames(df), [:G]) .=> last)

  println(combineddf)
end

function mwe()
  df = DataFrame(X = rand(4), Y=rand(4))

  local transformed

  try
    transformed = transform(cumsum, df)
  catch err
    @warn "(:) .=> last failed. Error: $err"
  end

  #workaround
  transformed = transform(df, propertynames(df) .=> cumsum)

  println(transformed)
end

mweold()=#

using DataFrames, BenchmarkTools, Dates, StatsBase
function mwedates()
  #build the sample
  dts = reduce(vcat, [[Date(2011,11,11) + Day(i) for j in 1:10^4] for i in 1:100])
  mdts = dts |> Vector{Union{Date, Missing}}
  id = reduce(vcat, [[j for j in 1:10^4] for i in 1:100])
  df = DataFrame(date=dts, mdate = mdts, id=id)

  #shuffle
  df = df[randperm(10^6), :]

  print("sort date and id: ")
  @btime sort($df, [:date, :id])
  print("sort date(with missings) and id: ")
  @btime sort($df, [:mdate, :id])
  print("work around performance: ")
  @btime begin
    $df.mdateconverted = $df.mdate |> Vector{Date}
    sort($df, [:mdateconverted, :id])
  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
Minimization wrapper for unsupervised learning.
"""
function minimize_unsupervised(L::LossUnsupervised, R::RegularizerUnsupervised, X, k,
                               beta=0.8, alpha=0.5, init=nothing, t_init=1.0;
                               max_iters=5000, verbose=false, tol=1e-8)

    theta_X, theta_Y, losses = optimize_unsupervised(L, R, X, k,
                                                     alpha, beta, init, t_init,
                                                     t_min=1e-15,
                                                     max_iters=max_iters, verbose=verbose, tol=tol)
    return theta_X, theta_Y, losses
end

"""
Alternating minimization for unsupervised problems.
"""
function optimize_unsupervised(L::LossUnsupervised, R::RegularizerUnsupervised,
                               C, k, alpha=.5, beta=.8, init=nothing, t_init=1.0;
                               t_min=1e-15, verbose=true, max_iters=5000, tol=1e-8)

    n, d = size(C)
    decay = true

    if verbose
        info("Solving problem. Shape is $n by $d.")
    end

    # convenience functions
    LOSS(u, v) = eval(L, C, u, v)
    GRAD_X(u, v) = deriv(L, C, u, v, "X")
    GRAD_Y(u, v) = deriv(L, C, u, v, "Y")
    PROX(u) = prox(R, u)

    theta_X, theta_Y, losses = [], [], []

    if init == nothing
        init_X = rand(n, k)
        init_Y = rand(d, k)
    else
        init_X = init["X"]
        init_Y = init["Y"]
    end

    if verbose
        info("The initial loss is $(LOSS(init_X, init_Y))")
        info("With gradient X $(GRAD_X(init_X, init_Y))")
        info("With gradient Y $(GRAD_Y(init_X, init_Y))")
    end

    assert(size(init_X)[1] == n)
    assert(size(init_Y)[1] == d)
    assert(size(init_X)[2] == k && size(init_Y)[2] == k)

    push!(theta_X, init_X)
    push!(theta_Y, init_Y)

    tx, ty = t_init, t_init

    progress_x, progress_y = true, true

    for k = 1:max_iters
        # Step for X
        curr_grad_x = GRAD_X(theta_X[end], theta_Y[end])
        grad_step_X = PROX(theta_X[end] - tx*curr_grad_x)

        progress_x = true
        prev_loss = LOSS(theta_X[end], theta_Y[end])

        while (LOSS(grad_step_X, theta_Y[end]) >= prev_loss)
            tx *= beta
            grad_step_X = PROX(theta_X[end] - tx*curr_grad_x)
            if tx < t_min
                progress_x = false
                break
            end
        end
        tx *= 2

        push!(theta_X, grad_step_X)

        # Step for Y
        curr_grad_y = GRAD_Y(theta_X[end], theta_Y[end])
        grad_step_Y = PROX(theta_Y[end] - ty*curr_grad_y)

        progress_y = true
        prev_loss = LOSS(theta_X[end], theta_Y[end])

        while (LOSS(theta_X[end], grad_step_Y) >= prev_loss)
            ty *= beta
            grad_step_Y = PROX(theta_Y[end] - ty*curr_grad_y)
            if ty < t_min
                progress_y = false
                break
            end
        end
        ty *= 2

        push!(theta_Y, grad_step_Y)
        push!(losses, LOSS(theta_X[end], theta_Y[end]))

        if verbose
            info("Iteration: $k,  Loss: $(losses[end])")
        end

        if (!progress_x && !progress_y)
            warn("Algorithm is not making any progress... breaking.")
            return theta_X, theta_Y, losses
        end

        if k > 4 && maximum(abs.(losses[end-4:end-1] - losses[end-3:end])) < tol
            info("Done.")
            return theta_X, theta_Y, losses
        end
    end
    print("Did not converge. Loss: $(losses[end])")
    return -1
end


text:
Base.@kwdef mutable struct Stop{TI, TE, TT, TS, TF, TX, TFT}
    it::Int = 0
    t0::Float64 = time()
    maxiter::TI = nothing
    maxfevals::TE = nothing
    maxtime::TT = nothing
    stagnation::TS = 100
    ftarget::TF = nothing
    xtol::TX = nothing
    ftol::TFT = 1e-11
    reason::Symbol = :none
end
Stop(n, λ; kwargs...) = Stop(; stagnation = floor(Int, 100 + 100 * n^1.5/λ), kwargs...)
function Base.show(io::IO, ::MIME"text/plain", s::Stop)
    println(io, "Stop")
    for f in fieldnames(Stop)
        v = getproperty(s, f)
        v !== nothing && println(io, "  $f = $v")
    end
end
start!(s) = s.t0 = time()
maxiter(::Nothing, ::Any, ::Any) = false
maxiter(m, it, s) = it >= m && (s.reason = :maxiter; true)
maxtime(::Nothing, ::Any, ::Any) = false
maxtime(m, t0, s) = time() - t0 > m && (s.reason = :maxtime; true)
maxfevals(::Nothing, ::Any, ::Any) = false
maxfevals(m, it, s) = it > m && (s.reason = :maxfevals; true)
ftarget(::Nothing, ::Any, ::Any) = false
ftarget(f, fbest, s) = fbest < f && (s.reason = :ftarget; true)
xtol(::Nothing, ::Any, ::Any, ::Any) = false
xtol(tol, p, sigma, s) = (&)((sigma * p .< tol)...) && (s.reason = :xtol; true)
ftol(::Nothing, ::Any, ::Any) = false
ftol(tol, logger, s) = length(logger.frange) > 10 && (&)((logger.frange[end-2:end] .< tol)...) && (s.reason = :ftol; true)
stagnation(::Nothing, ::Any, ::Any) = false
function stagnation(tol, logger, o)
    o.stop.it > o.p.n * (5 + 100/o.p.λ) &&
    length(logger.fmedian) >= 2*tol &&
    median(logger.fmedian[end-tol+1:end]) >= median(logger.fmedian[end-2tol+1:end-tol]) &&
    (o.stop.reason = :stagnation; true)
end
function terminate!(o::Optimizer)
    s = o.stop
    s.it += 1
    maxiter(s.maxiter, s.it, s) ||
    maxtime(s.maxtime, s.t0, s) ||
    maxfevals(s.maxfevals, s.it * o.p.λ + noisefevals(o.p.noise_handling), s) ||
    ftarget(s.ftarget, best(o.logger), s) ||
    xtol(s.xtol, o.p.cov.p, sigma(o.p), s) ||
    ftol(s.ftol, o.logger, s) ||
    stagnation(s.stagnation, o.logger, o)
end

max_stars_repo_name: jbrea/CMAEvolutionStrategy

max_issues_repo_name: jbrea/CMAEvolutionStrategy


text:
module TestResources

using AMGX, Defer, Test
using AMGX: Config, Resources

@scope @testset "Resources" begin
    cfg = @! Config("")
    resources = @! Resources(cfg)
end

end # module


text:
function _info(model::Optimizer, key::MOI.VariableIndex)
    if haskey(model.variable_info, key)
        return model.variable_info[key]
    end
    throw(MOI.InvalidIndex(key))
    return
end

function _make_var(model::Optimizer, variable::Variable)
    # Initialize `VariableInfo` with a dummy `VariableIndex` and a column,
    # because we need `add_item` to tell us what the `VariableIndex` is.
    index = CleverDicts.add_item(
        model.variable_info,
        VariableInfo(MOI.VariableIndex(0), variable),
    )
    _info(model, index).index = index
    return index
end

function _make_var(
    model::Optimizer,
    variable::Variable,
    set::MOI.AbstractScalarSet,
)
    index = _make_var(model, variable)
    return index,
    MOI.ConstraintIndex{MOI.VariableIndex, typeof(set)}(index.value)
end

function _make_vars(model::Optimizer, variables::Vector{<:Variable})
    # Barely used, because add_constrained_variables may have variable sets (except for AbstractVectorSet).
    # Only implemented in the unconstrained case.
    indices = Vector{MOI.VariableIndex}(undef, length(variables))
    for i in 1:length(variables)
        indices[i] = CleverDicts.add_item(
            model.variable_info,
            VariableInfo(MOI.VariableIndex(0), variables[i]),
        )
        _info(model, indices[i]).index = indices[i]
    end
    return indices
end

function _make_numvar(
    model::Optimizer,
    set::MOI.AbstractScalarSet;
    lb::Float64=-IloInfinity,
    ub::Float64=IloInfinity,
)
    return _make_var(model, cpo_java_numvar(model.inner, lb, ub), set)
end

function _make_intvar(
    model::Optimizer,
    set::MOI.AbstractScalarSet;
    lb::Int32=-IloMinInt,
    ub::Int32=IloMaxInt,
)
    vindex, cindex = _make_var(
        model,
        cpo_java_intvar(model.inner, Int32(lb), Int32(ub)),
        set,
    )
    _info(model, vindex).type = INTEGER
    return vindex, cindex
end

function _make_boolvar(model::Optimizer, set::MOI.AbstractScalarSet)
    vindex, cindex = _make_var(model, cpo_java_boolvar(model.inner), set)
    _info(model, vindex).type = BINARY
    return vindex, cindex
end

function supports_add_constrained_variables(
    ::Optimizer,
    ::Type{F},
) where {
    F <: Union{
        MOI.EqualTo{Float64},
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.Interval{Float64},
        MOI.EqualTo{Int},
        MOI.LessThan{Int},
        MOI.GreaterThan{Int},
        MOI.Interval{Int},
        MOI.ZeroOne,
        MOI.Integer,
    },
}
    return true
end

function MOI.add_variable(model::Optimizer)
    return _make_var(
        model,
        cpo_java_numvar(model.inner, -IloInfinity, IloInfinity),
    )
end

function MOI.add_variables(model::Optimizer, N::Int)
    return _make_vars(
        model,
        cpo_java_numvararray(model.inner, Int32(N), -IloInfinity, IloInfinity),
    )
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.GreaterThan{T},
) where {T <: Real}
    return _make_numvar(model, set, lb=set.lower)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.LessThan{T},
) where {T <: Real}
    return _make_numvar(model, set, lb=set.upper)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.EqualTo{T},
) where {T <: Real}
    return _make_numvar(model, set, lb=set.value, ub=set.value)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.Interval{T},
) where {T <: Real}
    return _make_numvar(model, set, lb=set.lower, ub=set.upper)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.GreaterThan{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.lower)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.LessThan{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.upper)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.EqualTo{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.value, ub=set.value)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.Interval{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.lower, ub=set.upper)
end

function MOI.add_constrained_variable(model::Optimizer, set::MOI.ZeroOne)
    return _make_boolvar(model, set)
end

function MOI.add_constrained_variable(model::Optimizer, set::MOI.Integer)
    return _make_intvar(model, set)
end

function MOI.is_valid(model::Optimizer, v::MOI.VariableIndex)
    return haskey(model.variable_info, v)
end

function MOI.delete(model::Optimizer, v::MOI.VariableIndex)
    # Actual deletion of the variable.
    info = _info(model, v)
    if info.n_socs > 0
        throw(MOI.DeleteNotAllowed(v))
    end
    cpo_java_remove(model.inner, info.variable)
    delete!(model.variable_info, v)

    # We throw away name_to_constraint so we will rebuild VariableIndex
    # constraint names (without v) on the next access to it.
    model.name_to_variable = nothing
    model.name_to_constraint = nothing

    return
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    if model.name_to_variable === nothing
        _rebuild_name_to_variable(model)
    end
    return get(model.name_to_variable, name, nothing)
end

function _rebuild_name_to_variable(model::Optimizer)
    model.name_to_variable = Dict{String, MOI.VariableIndex}()
    for (index, info) in model.variable_info
        if info.name == ""
            continue
        end

        if haskey(model.name_to_variable, info.name)
            error("Duplicate variable name detected: $(info.name)")
        end
        model.name_to_variable[info.name] = index
    end
    return
end

function MOI.get(model::Optimizer, ::MOI.VariableName, v::MOI.VariableIndex)
    return _info(model, v).name
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariableName,
    v::MOI.VariableIndex,
    name::String,
)
    info = _info(model, v)
    info.name = name
    cpo_java_addable_setname(model.inner, info.variable, name)

    model.name_to_variable = nothing
    model.name_to_constraint = nothing

    return
end


text:
lines = readlines(ARGS[1])
vals = [parse(Int, x) for x in split(lines[1], ",")]

function solution1(vals)
    min_val = min(vals...)
    max_val = max(vals...)

    target = nothing
    for candidate = min_val:max_val
        total_fuel = sum([abs(x - candidate) for x in vals])
        if isnothing(target)
            target = total_fuel
        elseif total_fuel < target
            target = total_fuel
        end
    end
    return target
end

function fuel_burnt(val1, val2)
    return sum([x for x = 1:abs(val1 - val2)])
end

function solution2(vals)
    min_val = min(vals...)
    max_val = max(vals...)

    target = nothing
    for candidate = min_val:max_val
        total_fuel = sum([fuel_burnt(x, candidate) for x in vals])
        if isnothing(target)
            target = total_fuel
        elseif total_fuel < target
            target = total_fuel
        end
    end
    return target
end

println(solution1(vals))
println(solution2(vals))


text:
using Test
using DataFrames, CSV, CategoricalArrays, HypothesisTests

path     = dirname(@__FILE__)
io       = IOBuffer();

#dat.li2007
#dat.hine1989
#dat.graves2010
#dat.bourassa1996

@testset "  Proportion Confidence Intarvals                          " begin
    ct = MetidaFreq.contab([38, 62])
    # Wald
    ci = MetidaFreq.propci(ct; level = 0.95, method = :wald)
    @test collect(ci)  ≈ [0.28486600512143206, 0.47513399487856794] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :wald)
    @test collect(ci)  ≈ [0.28486600512143206, 0.47513399487856794] atol=1E-6
    # Wald cc
    ci = MetidaFreq.propci(ct; level = 0.95, method = :waldcc)
    @test collect(ci)  ≈ [0.27986600512143206, 0.48013399487856795] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :waldcc)
    @test collect(ci)  ≈ [0.27986600512143206, 0.48013399487856795] atol=1E-6
    # Wilson
    ci = MetidaFreq.propci(ct; level = 0.95, method = :wilson)
    @test collect(ci)  ≈ [0.2909759925247873, 0.47790244704488943] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :wilson)
    @test collect(ci)  ≈ [0.2909759925247873, 0.47790244704488943] atol=1E-6
    # Wilson cc - check
    ci = MetidaFreq.propci(ct; level = 0.95, method = :wilsoncc)
    @test collect(ci)  ≈ [0.28639471634406627, 0.48294114679105093] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :wilsoncc)
    @test collect(ci)  ≈ [0.28639471634406627, 0.48294114679105093] atol=1E-6
    # Clopper-Pirson
    ci = MetidaFreq.propci(ct; level = 0.95, method = :cp)
    @test collect(ci)  ≈ [0.2847674761414794, 0.48253930575080645] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :cp)
    @test collect(ci)  ≈ [0.2847674761414794, 0.48253930575080645] atol=1E-6
    # SOC
    ci = MetidaFreq.propci(ct; level = 0.95, method = :soc)
    @test collect(ci)  ≈ [0.2891917018839231, 0.4775592393403469] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :soc)
    @test collect(ci)  ≈ [0.2891917018839231, 0.4775592393403469] atol=1E-6
    # Blaker - check!
    ci = MetidaFreq.propci(ct; level = 0.95, method = :blaker)
    @test collect(ci)  ≈ [0.2881924139040972, 0.4798218835425426] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :blaker)
    @test collect(ci)  ≈ [0.2881924139040972, 0.4798218835425426] atol=1E-6
    # arcsine
    ci = MetidaFreq.propci(ct; level = 0.95, method = :arc)
    @test collect(ci)  ≈ [0.2877714314998773, 0.47682358116201534] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :arc)
    @test collect(ci)  ≈ [0.2877714314998773, 0.47682358116201534] atol=1E-6
    # Jeffrey
    ci = MetidaFreq.propci(ct; level = 0.95, method = :jeffrey)
    @test collect(ci)  ≈ [0.2893837310094326, 0.4774506923359312] atol=1E-6
    ci = MetidaFreq.propci(38, 100; level = 0.95, method = :jeffrey)
    @test collect(ci)  ≈ [0.2893837310094326, 0.4774506923359312] atol=1E-6

    ############################################################################
    # default = wilson, level = 0.95
    ci = MetidaFreq.propci(MetidaFreq.Proportion(38, 100))
    @test collect(ci)  ≈ [0.2909759925247873, 0.47790244704488943] atol=1E-6
end

@testset "  Proportion difference Confidence Intarvals               " begin
    ct = MetidaFreq.contab([30 70; 40 50])
    # Miettinen & Nurminen
    ci = MetidaFreq.diffci(ct; level = 0.95, method = :mn)
    @test collect(ci)  ≈ [-0.2781290897168457, -0.006708341755865329] atol=1E-6
    # Wald
    ci = MetidaFreq.diffci(ct; level = 0.95, method = :wald)
    @test collect(ci)  ≈ [-0.28084842238, -0.00804046650] atol=1E-6
    # Wald CC
    ci = MetidaFreq.diffci(ct; level = 0.95, method = :waldcc)
    @test collect(ci)  ≈ [-0.29140397794, 0.00251508905] atol=1E-6
    # Newcombe Hybrid Score
    ci = MetidaFreq.diffci(ct; level = 0.95, method = :nhs)
    @test collect(ci)  ≈ [-0.275381800, -0.007158419] atol=1E-6
    # Newcombe Hybrid Score CC

    # Agresti-Caffo interval
    ci = MetidaFreq.diffci(ct; level = 0.95, method = :ac)
    @test collect(ci)  ≈ [-0.276944506, -0.006516705] atol=1E-6
    # ha

    # mnmee

    # mover

end

@testset "  Odd ratio Confidence Intarvals                           " begin
    ct = MetidaFreq.contab([30 70; 40 50])
    # Miettinen & Nurminen
    ci = MetidaFreq.orci(ct; level = 0.95, method = :mn)
    @test collect(ci)  ≈ [0.29537414, 0.97166697] atol=1E-6

    # woolf

    # awoolf

    #mover

end

@testset "  Risk ratio Confidence Intarvals                          " begin
    ct = MetidaFreq.contab([30 70; 40 50])
    # Miettinen & Nurminen
    ci = MetidaFreq.rrci(ct; level = 0.95, method = :mn)
    @test collect(ci)  ≈ [0.46636099123297575, 0.9799258384796817] atol=1E-6

    # cli

    # li

    #mover

end

@testset "  Contab from tabular data                                 " begin
    freqdat  = CSV.File(path*"/csv/freqdat.csv") |> DataFrame
    ct = MetidaFreq.contab(freqdat, :row, :col)
    @test ct.tab[1,1] == 21
    @test ct.tab[1,2] == 5
    @test ct.tab[2,1] == 8
    @test ct.tab[2,2] == 9

    Base.show(io, ct)

    HypothesisTests.ChisqTest(ct)

    HypothesisTests.MultinomialLRTest(ct)

    HypothesisTests.FisherExactTest(ct)

    transform!(freqdat, :row => categorical, renamecols=false)
    transform!(freqdat, :col => categorical, renamecols=false)

    ct = MetidaFreq.contab(freqdat, :row, :col)
    @test ct.tab[1,1] == 8
    @test ct.tab[1,2] == 9
    @test ct.tab[2,1] == 21
    @test ct.tab[2,2] == 5

    freqdat2  = CSV.File(path*"/csv/ft.csv") |> DataFrame

    ct = MetidaFreq.contab(freqdat2, :row, :col)
    @test ct.tab[1,1] == 20
    @test ct.tab[1,2] == 24
    @test ct.tab[2,1] == 83
    @test ct.tab[2,2] == 44

    ct = MetidaFreq.contab(freqdat2, :row, :col; sort = :s1)

    @test ct[1].tab[1,1] == 0
    @test ct[1].tab[1,2] == 2
    @test ct[1].tab[2,1] == 65
    @test ct[1].tab[2,2] == 34

    @test ct[2].tab[1,1] == 20
    @test ct[2].tab[1,2] == 22
    @test ct[2].tab[2,1] == 18
    @test ct[2].tab[2,2] == 10

    ct = MetidaFreq.contab(freqdat2, :row, :col; sort = [:s1, :s2])
    @test  length(ct) == 6

    @test ct[6].id[:s1] == "e"
    @test ct[6].id[:s2] == "i"

    @test ct[6].tab[1,1] == 6
    @test ct[6].tab[1,2] == 12
    @test ct[6].tab[2,1] == 18
    @test ct[6].tab[2,2] == 10

    transform!(freqdat2, :row => categorical, renamecols=false)
    transform!(freqdat2, :col => categorical, renamecols=false)
    transform!(freqdat2, :s1 => categorical, renamecols=false)
    transform!(freqdat2, :s2 => categorical, renamecols=false)

    ct = MetidaFreq.contab(freqdat2, :row, :col; sort = [:s1, :s2])

    @test ct[5].tab[1,1] == 12
    @test ct[5].tab[1,2] == 6
    @test ct[5].tab[2,1] == 10
    @test ct[5].tab[2,2] == 18

end

@testset "  Meta proportions                                         " begin
    # Make DataSet
    ############################################################################
    pf1 = MetidaFreq.contab([15 8; 5 14])
    pf2 = MetidaFreq.contab([45 72; 23 95])
    mds = MetidaFreq.DataSet([pf1, pf2])
    ############################################################################

    # Check all metrics
    mp = MetidaFreq.metaprop(mds, :rr)
    mp = MetidaFreq.metaprop(mds, :or)
    mp = MetidaFreq.metaprop(mds, :diff)

    # Fixed effect MH (diff)
    mp  = MetidaFreq.metaprop(mds, :diff)
    mpf = MetidaFreq.metapropfixed(mp; weights = :mh)
    @test mpf.est ≈ 0.2196889005445779 atol=1E-6
    @test mpf.var ≈ 0.0028728509817330943 atol=1E-6
    ci = confint(mpf; level = 0.95)
    @test collect(ci)  ≈ [0.1146368241999458, 0.32474097688921] atol=1E-6

    # Fixed effect IV (diff)
    mp  = MetidaFreq.metaprop(mds, :diff)
    mpf = MetidaFreq.metapropfixed(mp; weights = :iv)

    # Check random effect tau (diff)
    mp  = MetidaFreq.metaprop(mds, :diff)
    mpf = MetidaFreq.metaproprandom(mp; tau = :dl)
    mpf = MetidaFreq.metaproprandom(mp; tau = :ho)
    mpf = MetidaFreq.metaproprandom(mp; tau = :hm)
    mpf = MetidaFreq.metaproprandom(mp; tau = :sj)

end

@testset "  Goodman CI                                               " begin

    ci = MetidaFreq.ci_prop_goodman([44,55,43,32,67,78], 0.05)
    @test ci[1][1] ≈ 0.09468368035184335  atol=1E-6
    @test ci[1][2] ≈ 0.1966412812730604  atol=1E-6
    @test ci[6][1] ≈ 0.18692731624998904  atol=1E-6
    @test ci[6][2] ≈ 0.3130119423857655  atol=1E-6

    ci = MetidaFreq.ci_prop_goodman([91,49,37,43], 0.05)
    @test collect(ci[1]) ≈ [0.3342024783268433, 0.4978332073846343]  atol=1E-6
    @test collect(ci[3]) ≈ [0.11455134046157951, 0.24011208358778185]  atol=1E-6

    ct = MetidaFreq.contab(permutedims([91,49,37,43]))
    ci = MetidaFreq.mpropci(ct; level = 0.95, method = :default)
    @test collect(ci[1]) ≈ [0.3342024783268433, 0.4978332073846343]  atol=1E-6
    @test collect(ci[3]) ≈ [0.11455134046157951, 0.24011208358778185]  atol=1E-6
end


text:
@testset "MultipleBranches" begin

# ((freq THz)) ((IR Activity / e^2 amu^-1))
# These data from MAPbI3-Cubic_PeakTable.csv
# https://github.com/WMD-group/Phonons/tree/master/2015_MAPbI3/SimulatedSpectra
# Data published in Brivio2015 (PRB)
# https://doi.org/10.1103/PhysRevB.92.144308
MAPI= [
96.20813558773261 0.4996300522819191
93.13630357703363 1.7139631746083817
92.87834578121567 0.60108592692181
92.4847918585963 0.0058228799414729
92.26701437594754 0.100590086574602
89.43972834606603 0.006278895133832249
46.89209141511332 0.2460894564364346
46.420949316788 0.14174282581124137
44.0380222871706 0.1987196948553428
42.89702947649343 0.011159939465770681
42.67180170168193 0.02557751102757614
41.46971205834201 0.012555230726601503
37.08982543385215 0.00107488277468418
36.53555265689563 0.02126940080871224
30.20608114002676 0.009019481779712388
27.374810898415028 0.03994453721421388
26.363055017011728 0.05011922682554448
9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05   # Acoustic modes!
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]

vol=(6.29E-10)^3
ϵ_o=4.5
meff=0.12

ϵ_i=IRtoDielectric(MAPI,vol)
#ϵ_ired=ϵ_i/ϵ_0
#ϵ_s=ϵ_ired + ϵ_o

ϵ_s=sum(ϵ_i)+ϵ_o # total (static) dielectric = sum of ionic, and optical

println("Sum of ionic dielectric: $(ϵ_s)")

IRtoalpha(MAPI, volume=vol, ϵ_o=ϵ_o, ϵ_s=ϵ_s, meff=meff)

println()
splat=DielectricFromIRmode.(eachrow(MAPI), volume=vol)
println(splat)
println("Sum of dieletric: ", sum(splat))

f_dielectric=hcat( MAPI[:,1], splat)
println("f_dielectric: ",f_dielectric)

alphas=frohlichPartial.(eachrow(f_dielectric), ϵ_o = ϵ_o, ϵ_s = ϵ_o+sum(splat), meff=meff)
println("alphas: ",alphas)
println("sum alphas: ", sum(alphas))

println("Feynman v,w for alphas: ", feynmanvw.(alphas))

mobilityproblem=hcat(alphas, feynmanvw.(alphas), MAPI[:,1])
println("Specify mobility problem: ",mobilityproblem)

inverse_μ=Hellwarth1999mobilityRHS.(eachrow(mobilityproblem), meff, 300)

μ=sum(inverse_μ)^-1
@printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
println()

ħ=1.05457162825e-34
kB=1.3806504e-23

MAPI=[
      5     20.0
      0.5   0.67
     ] # fictious two-frequency MAPI-a-like

splat=DielectricFromIRmode.(eachrow(MAPI), volume=vol)
println(splat)
println("Sum of dieletric: ", sum(splat))

f_dielectric=hcat( MAPI[:,1], splat)

for T in 10:1:500
    alphas=frohlichPartial.(eachrow(f_dielectric), ϵ_o = ϵ_o, ϵ_s = ϵ_o+sum(splat), meff=meff)
    βred=ħ*MAPI[:,1]*1E12*2π/(kB*T)

    mobilityproblem=hcat(alphas, feynmanvw.(alphas, βred), MAPI[:,1])
    inverse_μ=Hellwarth1999mobilityRHS.(eachrow(mobilityproblem), meff, T)
    μ=sum(inverse_μ)^-1

    @printf("\nT= %d μ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",T,μ,μ*100^2)
end

println()
# A bit cyclical, but this is extrapolating back to if MAPI had a single LO
# mode, at the point where the Hellwarth1999 averaging approximation puts it.
MAPI_singlemode = [
2.25 1.67502750447212449732 ]

# By comparison to 2017 PRB, this should be 24.1-4.5 = 19.6
#  Nb: however, that dielectric is from my 2014 NanoLetters, whereas the
#  2015PRB above is a different calculation, with different convergence.
ϵ_i=IRtoDielectric(MAPI_singlemode,vol)
ϵ_s=sum(ϵ_i)+ϵ_o
println("Sum of ionic dielectric: $(ϵ_s)")
_,αmode_MAPIe=IRtoalpha(MAPI_singlemode,volume=vol, ϵ_o=ϵ_o, ϵ_s=ϵ_s, meff=meff)

α_MAPIe=frohlichalpha(4.5, 24.1, 2.25E12, meff)

@test αmode_MAPIe ≈ α_MAPIe atol=0.01

f,α = IRtoalpha(MAPI, volume=vol, ϵ_o=ϵ_o, ϵ_s=ϵ_s, meff=meff)

end # @testset


text:
module Embedding
export
    umap,
    tsne

    include("umap.jl")
    include("tsne.jl")

end


text:
# @doc doc"""function lrtest(Ns::Int,N::Array{Int,1},L::Array{Int,1},
#                dataType::Array{DataType,1},distType)
#   Monte-Carlo evaluation of lattice reduction techniques. NxN matrices of
#   dataType are decomposed with each LR technique, with execution time,
#   memory, orthogonality defect, and Hermite factor shown.  Accepted
#   datatypes include Integer types such as Int8; FloatingPoint types such
#   as Float32; and Complex versions of these types. Accepted distTypes are
#   "randn" (Gaussian distributed with variance L), "rand" (uniformly
#   distributed over -L:L), "colUnif" (first column is uniform over -L:L).
#
# Examples
#   # As function of Integer width
#   lrtest(50,2.^[4;],[100;],[Int32,Int64,Int128,Float64,BigInt,BigFloat],"rand")
#   # As function of N
#   lrtest(50,2.^[0:6;],[100;],[Float64],"rand")
# """
function lrtest(Ns::Int,N::Array{Int,1},L::Array{Int,1},
                dataType::Array{DataType,1},distType)
#                dataType,distType)

# Packages that need to be loaded for lrtest to work include PyPlot, and
# CPUTime.

#lrAlgs = [lll, lllrecursive,seysen]
lrAlgs = [lll, seysen]

@printf("      Ns      N      L   dataType")
for ax = 1:min(length(lrAlgs),6)
    @printf(" %13s",lrAlgs[ax])
end
@printf("\n")

xtickStrs = [];
out = cell(0)
if length(N)>1
    for s = 1:length(N)
        push!(out,lrsim(Ns,N[s],L[1],dataType[1], distType, lrAlgs));
    end
    plotfun = PyPlot.loglog;
    xval = N;
    xlab = "Matrix Size";
    tstr = @sprintf("Ns=%d,L=%d,Type=%s,dist=%s",
                    Ns,L[1],string(dataType[1]),distType);
elseif length(L)>1
    for s = 1:length(L)
        push!(out,lrsim(Ns,N[1],L[s],dataType[1], distType, lrAlgs));
    end
    plotfun = semilogy;
    xval = L;
    xlab = "L";
    tstr = @sprintf("Ns=%d,N=%d,Type=%s,dist=%s",
                    Ns,N[1],string(dataType[1]),distType);
elseif length(dataType)>1
    for s = 1:length(dataType)
        push!(out,lrsim(Ns,N[1],L[1],dataType[s], distType, lrAlgs));
    end
    plotfun = semilogy;
    xval = 1:length(dataType)
    xtickStrs = map(string,dataType)
    xlab = "dataType";
    tstr = @sprintf("Ns=%d,N=%d,L=%d,dist=%s",Ns,N[1],L[1],distType);
end

pColor = {"r-","b.-","k-","g-","c-","m-",
          "r--","b.-.","k--","g.-.","c--","m*-.",
          "rs--","gp-.","bv-","kh--","c+-.","m.-",};
pIdx   = 1;

Nout = size(out,1);

f=figure(num=1,figsize=(6.5,4.5))
clf()
plt.style[:use]("ggplot")

times = zeros(Nout);
for a=1:length(out[1][2])
    for k=1:Nout
        times[k] = out[k][1][a];
    end
    plotfun(xval,times,pColor[pIdx],label=out[1][2][a]);
    hold(true);
    pIdx = pIdx==length(pColor)? 1: pIdx + 1;
end
hold(false);

xlabel(xlab);
ylabel("execution time (sec)");
legend(loc=2);
#title(tstr)
if ~isempty(xtickStrs)
    xticks(xval,xtickStrs)
end
grid(figure=f,which="both",axis="y")
grid(figure=f,which="major",axis="x")
plt.tight_layout()

return
end

#######################################################################
function lrsim(Ns,N,L,dataType,distType,lrAlgs)
# sortDict = Dict("Insertion" => d->sort(d,alg=InsertionSort),

distDict =     {"rand"=> d -> rand!(-L:L,d),
                "randmax"=> d -> rand(dataType,N,N),
                "randn"=> d -> L*randn!(d::AbstractArray{dataType,2}),
                 };
# distDict =     {"rand"=> d -> rand!(-L:L,d),
#                 "randmax"=> d -> rand!(dataType,d),
#                 "randn"=> d -> L*randn!(d::AbstractArray{dataType,2}),
#                  };
randf! = distDict[distType];

Nalgs = length(lrAlgs);
times = ones(Nalgs,Ns)*Inf;
hermitef = ones(Nalgs,Ns)*Inf;
orthf = ones(Nalgs,Ns)*Inf;
algNames = cell(Nalgs)

data = Array(dataType,N,N);

for ix = 1:Ns
    cnum = Inf;
    while cnum>1e5  # With Ints need to check for singular matrices
        data = randf!(data);
        if !(dataType<:BigInt || dataType<:BigFloat)
            cnum = cond(data)
        else
            cnum = 0
        end
    end

    #
    # Lattice-reduction algorithms
    #
    for ax = 1:length(lrAlgs)
        CPUtic();
        (B,T) = lrAlgs[ax](data);
        times[ax,ix] = CPUtoq();
        # detB = abs(det(B))
        # # Hermite factor
        # hermitef[ax,ix] = norm(B[:,1])/detB^(1/N)
        # # Orthogonality defect
        # orthf[ax,ix] = prodBi(B,N)/detB

        # if abs(abs(det(T))-1.0)>1e-6
        #     println("For $(lrAlgs[ax]), det(T) is $(abs(det(T)))")
        #     times[ax,ix] = Inf
        # elseif abs(sum(B-data*T))>1e-6
        #     println("For $(lrAlgs[ax]), sum(B-H*T) "*
        #             "is $(abs(sum(B-data*T))), ")
        #     times[ax,ix] = Inf
        # end
    end
end

#algNames[Nalgs] = "randf";
for ax = 1:length(lrAlgs)
    algNames[ax] = string(lrAlgs[ax]);
end
@printf("%8d %6d %6d %10s", Ns,  N,  L, string(dataType))

#outtimes = mean(times,2);
outtimes = zeros(Nalgs,1);
outHF = zeros(Nalgs,1);
outOF = zeros(Nalgs,1);
for ax = 1:Nalgs
#    stimes = sort(vec(times[ax,:]));
#    outtimes[ax] = mean(stimes[1:floor(Ns*2/3)]);
    outtimes[ax] = mean(times[ax,:]);
    outHF[ax] = mean(hermitef[ax,:]);
    outOF[ax] = mean(orthf[ax,:]);
end

# BOLD= "\x1B[1;30m"
# RESET="\x1B[0;0m"
# (mn,mx)=findmin(outtimes)
for ax = 1:min(Nalgs,7)
    # if ax==mx
    #     @printf("    %s%10.4f%s",BOLD,outtimes[ax]*1000,RESET)
    # else
        @printf("    %10.4f",outtimes[ax]*1000)
#        @printf("    %10.4f",outHF[ax])
    # end
end
@printf("\n")

return outtimes, algNames
end

#######################################################################
function prodBi(B,N)
    prod = 1.0
    for n=1:N
        pn = 0.0
        for m=1:N
            pn+=conj(B[m,n])*B[m,n]
        end
        println("pn=$(pn), typeof(pn)=$(typeof(pn))")
        prod*=sqrt(pn)
    end
    return prod
end


text:
a = BigInt("123456789012345678901234567890")
b = BigInt("123456789012345678901234567891")

@test typeof(a+1) == BigInt
@test a+1 == b
@test isequal(a+1, b)
@test b == a+1
@test !(b == a)
@test b > a
@test b >= a
@test !(b < a)
@test !(b <= a)

c = BigInt("246913578024691357802469135780")
@test typeof(a * 2) == BigInt
@test a*2 == c
@test c-a == a
@test c == a + a
@test c+1 == a+b

d = BigInt("-246913578024691357802469135780")
@test typeof(d) == BigInt
@test d == -c

ee = typemax(Int64)
@test typeof(BigInt(ee)) == BigInt
@test BigInt(ee)+1 == BigInt("9223372036854775808")

#Multiple calls for sanity check, since we're doing direct memory manipulation
@test string(a) == "123456789012345678901234567890"
@test string(b) == "123456789012345678901234567891"
@test string(c) == "246913578024691357802469135780"
@test string(d) == "-246913578024691357802469135780"
@test string(a) == "123456789012345678901234567890"

for i = -10:10, j = [-10:-1,1:10]
    @test div(BigInt(i), BigInt(j)) == div(i,j)
    @test fld(BigInt(i), BigInt(j)) == fld(i,j)
    @test mod(BigInt(i), BigInt(j)) == mod(i,j)
    @test rem(BigInt(i), BigInt(j)) == rem(i,j)
end

@test typeof(BigInt(typemax(Int8))) == BigInt
@test typeof(BigInt(typemax(Int16))) == BigInt
@test typeof(BigInt(typemax(Int32))) == BigInt
@test typeof(BigInt(typemax(Int64))) == BigInt
@test typeof(BigInt(typemax(Int128))) == BigInt

@test typeof(BigInt(true)) == BigInt
@test typeof(BigInt(typemax(Uint8))) == BigInt
@test typeof(BigInt(typemax(Uint16))) == BigInt
@test typeof(BigInt(typemax(Uint32))) == BigInt
@test typeof(BigInt(typemax(Uint64))) == BigInt
@test typeof(BigInt(typemax(Uint128))) == BigInt

@test typeof(BigInt(BigInt(1))) == BigInt

# Signed addition
@test a+int8(1) == b
@test a+int16(1) == b
@test a+int32(1) == b
@test a+int64(1) == b
@test a+int128(1) == b
@test int8(1)+ a == b
@test int16(1)+a == b
@test int32(1)+a == b
@test int64(1)+a == b
@test b+int8(-1) == a
@test b+int16(-1) == a
@test b+int32(-1) == a
@test b+int64(-1) == a
@test int8(-1)+ b == a
@test int16(-1)+b == a
@test int32(-1)+b == a
@test int64(-1)+b == a

# Unsigned addition
@test a+true == b
@test a+uint8(1) == b
@test a+uint16(1) == b
@test a+uint32(1) == b
@test a+uint64(1) == b
@test a+uint128(1) == b
@test true+a == b
@test uint8(1)+ a == b
@test uint16(1)+a == b
@test uint32(1)+a == b
@test uint64(1)+a == b

# Signed subtraction
@test b-int8(1) == a
@test b-int16(1) == a
@test b-int32(1) == a
@test b-int64(1) == a
@test int8(1)- b == -a
@test int16(1)-b == -a
@test int32(1)-b == -a
@test int64(1)-b == -a
@test a-int8(-1) == b
@test a-int16(-1) == b
@test a-int32(-1) == b
@test a-int64(-1) == b
@test int8(-1)- a == -b
@test int16(-1)-a == -b
@test int32(-1)-a == -b
@test int64(-1)-a == -b

# Unsigned subtraction
@test b-true == a
@test b-uint8(1) == a
@test b-uint16(1) == a
@test b-uint32(1) == a
@test b-uint64(1) == a
@test true-b == -a
@test uint8(1)- b == -a
@test uint16(1)-b == -a
@test uint32(1)-b == -a
@test uint64(1)-b == -a

# Signed multiplication
@test a*int8(1) == a
@test a*int16(1) == a
@test a*int32(1) == a
@test a*int64(1) == a
@test int8(1)* a == a
@test int16(1)*a == a
@test int32(1)*a == a
@test int64(1)*a == a
@test a*int8(-1) == -a
@test a*int16(-1) == -a
@test a*int32(-1) == -a
@test a*int64(-1) == -a
@test int8(-1)* a == -a
@test int16(-1)*a == -a
@test int32(-1)*a == -a
@test int64(-1)*a == -a

# Unsigned multiplication
@test a*true == a
@test a*uint8(1) == a
@test a*uint16(1) == a
@test a*uint32(1) == a
@test a*uint64(1) == a
@test true*a == a
@test uint8(1)* a == a
@test uint16(1)*a == a
@test uint32(1)*a == a
@test uint64(1)*a == a

@test a+BigInt(1) == b

@test BigInt(5) << 3 == 40
@test BigInt(5) >> 1 == 2
@test BigInt(-5) << 3 == -40
@test BigInt(-5) >> 1 == -3

@test ~BigInt(123) == -124
@test BigInt(123) & BigInt(234) == 106
@test BigInt(123) | BigInt(234) == 251
@test BigInt(123) $ BigInt(234) == 145

@test gcd(BigInt(48), BigInt(180)) == 12
@test lcm(BigInt(48), BigInt(180)) == 720

@test factorial(BigInt(40)) == BigInt("815915283247897734345611269596115894272000000000")
@test binomial(BigInt(-53), 42) == BigInt("959509335087854414441273718")
@test binomial(BigInt(113), BigInt(42)) == BigInt("18672199984318438125634054194360")

a = rand(1:100, 10000)
b = map(BigInt, a)
@test sum(a) == sum(b)

# Iterated arithmetic
a = BigInt("315135")
b = BigInt("12412")
c = BigInt("3426495623485904783478347")
d = BigInt("-1398984130")
f = BigInt("2413804710837418037418307081437315263635345357386985747464")
g = BigInt("-1")

@test +(a, b) == BigInt("327547")
@test +(a, b, c) == BigInt("3426495623485904783805894")
@test +(a, b, c, d) == BigInt("3426495623485903384821764")
@test +(a, b, c, d, f) == BigInt("2413804710837418037418307081437318690130968843290370569228")
@test +(a, b, c, d, f, g) == BigInt("2413804710837418037418307081437318690130968843290370569227")

@test *(a, b) == BigInt("3911455620")
@test *(a, b, c) == BigInt("13402585563389346256121263521460140")
@test *(a, b, c, d) == BigInt("-18750004504148804423388563022070650287578200")
@test *(a, b, c, d, f) == BigInt("-45258849200337190631492857400003938881995610529251881450243326128168934937055005474972396281351684800")
@test *(a, b, c, d, f, g) == BigInt("45258849200337190631492857400003938881995610529251881450243326128168934937055005474972396281351684800")

@test ($)(a, b) == BigInt("327299")
@test ($)(a, b, c) == BigInt("3426495623485904783798472")
@test ($)(a, b, c, d) == BigInt("-3426495623485906178489610")
@test ($)(a, b, c, d, f) == BigInt("-2413804710837418037418307081437316711364709261074607933698")
@test ($)(a, b, c, d, f, g) == BigInt("2413804710837418037418307081437316711364709261074607933697")

@test (&)(a, b) == BigInt("124")
@test (&)(a, b, c) == BigInt("72")
@test (&)(a, b, c, d) == BigInt("8")
@test (&)(a, b, c, d, f) == BigInt("8")
@test (&)(a, b, c, d, f, g) == BigInt("8")

@test (|)(a, b) == BigInt("327423")
@test (|)(a, b, c) == BigInt("3426495623485904783802111")
@test (|)(a, b, c, d) == BigInt("-1396834561")
@test (|)(a, b, c, d, f) == BigInt("-1358954753")
@test (|)(a, b, c, d, f, g) == BigInt("-1")

@test isprime(BigInt(1000000007))
@test isprime(BigInt(1000000007), 1)
@test isprime(BigInt(10000000019))
@test isprime(BigInt("359334085968622831041960188598043661065388726959079837"))
@test !isprime(BigInt(1))
@test !isprime(BigInt(10000000020))

@test trailing_ones(a) == 8
@test trailing_zeros(b) == 2
@test count_ones(a) == 14

# Large Fibonacci to exercise BigInt
# from Bill Hart, https://groups.google.com/group/julia-dev/browse_frm/thread/798e2d1322daf633
function mul(a::Vector{BigInt}, b::Vector{BigInt})
   x = a[2]*b[2]
   c = Array(BigInt,3)
   c[1] = a[1]*b[1] + x
   c[2] = a[1]*b[2] + a[2]*b[3]
   c[3] = x + a[3]*b[3]
   return c
end

function bigfib(n)
    n == 0 && return BigInt(0)
    n -= 1
    r = [BigInt(1), BigInt(1), BigInt(0)]
    s = [BigInt(1), BigInt(0), BigInt(1)]
    while true
       (n & 1) == 1 && (s = mul(s,r))
       (n >>= 1) == 0 && return s[1]
       r = mul(r,r)
    end
end
@test [bigfib(n) for n=0:10] == [0,1,1,2,3,5,8,13,21,34,55]

n = bigfib(1000001)
@test ndigits(n) == 208988
@test mod(n,big(10)^15) == 359244926937501
@test div(n,big(10)^208973) == 316047687386689

s = string(n)
@test length(s) == 208988
@test endswith(s, "359244926937501")
@test beginswith(s, "316047687386689")

# serialization (#5133)
let
    n = BigInt("359334085968622831041960188598043661065388726959079837")
    b = IOBuffer()
    serialize(b,n)
    seek(b,0)
    @test deserialize(b) == n
end

# issue #5873
@test ndigits(big(90)) == 2
@test ndigits(big(99)) == 2
ndigits_mismatch(n) = ndigits(n) != ndigits(BigInt(n))
@test !any(ndigits_mismatch, 8:9)
@test !any(ndigits_mismatch, 64:99)
@test !any(ndigits_mismatch, 512:999)
@test !any(ndigits_mismatch, 8192:9999)

# conversion from float
@test BigInt(2.0) == BigInt(2.0f0) == BigInt(big(2.0)) == 2
@test_throws convert(BigInt, 2.1)
@test_throws convert(BigInt, big(2.1))

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using BeaData
using Test

@testset "metadata methods" begin
    ds = bea_datasets()
    @test "NIPA" ∈ ds[!, :DatasetName]

    pl = bea_parameterlist("FixedAssets")
    @test "TableName" ∈ pl[!, :ParameterName]

    pv = bea_parametervalues("NIUnderlyingDetail", "TableName")
    @test "U001A" ∈ pv[!, :Value]
end

@testset "data retrieval" begin
    nipa111 = bea_table("NIPA", "T10101", "A", 0 , 0)
    @test typeof(nipa111) <: BeaTable

    ud1a = bea_table("NIUnderlyingDetail", "U001A", "Q", 1996, 1996)
    @test typeof(ud1a) <: BeaTable

    faat101 = bea_table("FixedAssets", "FAAt101", "A", 0, 0)
    @test typeof(faat101) <: BeaTable
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdprogress

"""
    'feffery'_antdprogress(;kwargs...)

An AntdProgress component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `format` (optional): . format has the following type: lists containing elements 'prefix', 'suffix', 'content'.
Those elements have the following types:
  - `prefix` (String; optional)
  - `suffix` (String; optional)
  - `content` (String; optional)
- `gapDegree` (Real; optional)
- `gapPosition` (a value equal to: 'top', 'bottom', 'left', 'right'; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `percent` (Real; optional)
- `showInfo` (Bool; optional)
- `size` (a value equal to: 'default', 'small'; optional)
- `status` (a value equal to: 'success', 'exception', 'normal', 'active'; optional)
- `steps` (Real; optional)
- `strokeColor` (optional): . strokeColor has the following type: String | lists containing elements 'from', 'to'.
Those elements have the following types:
  - `from` (String; optional)
  - `to` (String; optional)
- `strokeLinecap` (a value equal to: 'round', 'square'; optional)
- `strokeWidth` (Real; optional)
- `style` (Dict; optional)
- `trailColor` (String; optional)
- `type` (a value equal to: 'line', 'circle', 'dashboard'; optional)
- `width` (Real; optional)
"""
function 'feffery'_antdprogress(; kwargs...)
        available_props = Symbol[:id, :className, :format, :gapDegree, :gapPosition, :loading_state, :percent, :showInfo, :size, :status, :steps, :strokeColor, :strokeLinecap, :strokeWidth, :style, :trailColor, :type, :width]
        wild_props = Symbol[]
        return Component("'feffery'_antdprogress", "AntdProgress", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
@symbol_func function G_PB(cur_reactor::AbstractReactor)
  cur_numerator = ( 1 - K_CD(cur_reactor) * cur_reactor.sigma_v )
  cur_denominator = ( 1 - K_CD(cur_reactor) * cur_reactor.sigma_v )

  cur_numerator /= sqrt(cur_reactor.T_bar) ^ 3
  cur_denominator /= sqrt(cur_reactor.T_bar) ^ 3

  cur_numerator *= K_omega(cur_reactor)
  cur_denominator *= K_omega(cur_reactor)

  cur_numerator += K_P(cur_reactor) * cur_reactor.sigma_v
  cur_denominator += K_P(cur_reactor) * cur_reactor.sigma_v

  # cur_numerator -= K_BR(cur_reactor) * sqrt(cur_reactor.T_bar)
  cur_denominator -= K_BR(cur_reactor) * sqrt(cur_reactor.T_bar)

  ( cur_numerator > 0 ) || return NaN

  cur_numerator ^= cur_reactor.mode_scaling[:P]

  cur_G = cur_reactor.T_bar

  cur_G *= cur_numerator

  cur_G /= cur_denominator

  cur_G
end


text:
default_namestrat(f) = name_runningnr()
default_namestrat(g::CompGraph) = default_namestrat(vertices(g))
function default_namestrat(vs::AbstractVector{<:AbstractVertex})
    # Even if all vertices have unique names, we can't be certain that no vertex produces more than one node
    # Therefore, we must take a pass through name_runningnr for each op even after we have mapped the vertex to a nextname. This is the reason for the v -> f -> namegen(v, "") wierdness
    namegen = name_runningnr()
    all(isnamed, vs) && length(unique(name.(vs))) == length(name.(vs)) && return v -> f -> namegen(v, "")
    ng(v::AbstractVertex) = namegen
    ng(f) = namegen(f)
    return ng
end

isnamed(v::AbstractVertex) = isnamed(base(v))
isnamed(v::CompVertex) = false
isnamed(v::InputVertex) = true
isnamed(v::SourceVertex) = true
isnamed(v::MutationVertex) = isnamed(trait(v))
isnamed(t::DecoratingTrait) = isnamed(base(t))
isnamed(t) = false
isnamed(::NamedTrait) = true

function name_runningnr(namefun = genname)
    exists = Set{String}()

    return function(f, init="_0")
        bname = namefun(f)
        candname = bname * init
        next = -1
        while candname in exists
            next += 1
            candname = bname * "_" * string(next)
        end
        push!(exists, candname)
        return candname
    end
end

genname(v::AbstractVertex) = name(v)
genname(f::F) where F = lowercase(string(nameof(F)))
genname(s::AbstractString) = s
genname(f::Function) = lowercase(string(f))

recursename(f, namestrat) = recursename(f, namestrat(f))
recursename(f, fname::String) = fname


text:
"""
    TranslationUnit(idx, source, args)
    TranslationUnit(idx, source, args, unsavedFiles, options)
Parse the given source file and the translation unit corresponding to that file.
"""
mutable struct TranslationUnit
    ptr::CXTranslationUnit
    idx::Index
    function TranslationUnit(idx::Index, source_filename, command_line_args, num_command_line_args, unsaved_files, num_unsaved_files, options)
        ptr = clang_parseTranslationUnit(idx, source_filename, command_line_args, num_command_line_args, unsaved_files, num_unsaved_files, options)
        @assert ptr != C_NULL "failed to parse file: $source_filename"
        obj = new(ptr, idx)
        finalizer(obj) do x
            if x.ptr != C_NULL
                clang_disposeTranslationUnit(x)
                x.ptr = C_NULL
            end
        end
        return obj
    end
end
TranslationUnit(idx, source, args, unsavedFiles, options) = TranslationUnit(idx, source, args, length(args), unsavedFiles, length(unsavedFiles), options)
TranslationUnit(idx, source, args, options) = TranslationUnit(idx, source, args, length(args), C_NULL, 0, options)
TranslationUnit(idx, source, args) = TranslationUnit(idx, source, args, length(args), C_NULL, 0, CXTranslationUnit_None)

Base.unsafe_convert(::Type{CXTranslationUnit}, x::TranslationUnit) = x.ptr

"""
    spelling(tu::TranslationUnit) -> String
Return the original translation unit source file name.
"""
function spelling(tu::TranslationUnit)
    cxstr = clang_getTranslationUnitSpelling(tu)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    getcursor(tu::TranslationUnit) -> CLCursor
    getcursor(tu::CXTranslationUnit) -> CXCursor
Return the cursor that represents the given translation unit.
"""
getcursor(tu::CXTranslationUnit) = clang_getTranslationUnitCursor(tu)
getcursor(tu::TranslationUnit)::CLCursor = clang_getTranslationUnitCursor(tu)

## TODO:
# clang_parseTranslationUnit2
# clang_createTranslationUnit2
# clang_parseTranslationUnit2FullArgv
# clang_saveTranslationUnit
# clang_suspendTranslationUnit
# clang_defaultSaveOptions
# clang_defaultEditingTranslationUnitOptions
# clang_defaultReparseOptions
# clang_reparseTranslationUnit

# helper
"""
    parse_header(header::AbstractString; index::Index=Index(), args::Vector{String}=String[],
                 includes::Vector{String}=String[], flags=CXTranslationUnit_None) -> TranslationUnit
Return the TranslationUnit for a given header. This is the main entry point for parsing.
See also [`parse_headers`](@ref).

# Arguments
- `header::AbstractString`: the header file to parse.
- `index::Index`: CXIndex pointer (pass to avoid re-allocation).
- `args::Vector{String}`: compiler switches as string array, eg: ["-x", "c++", "-fno-elide-type"].
- `includes::Vector{String}`: vector of extra include directories to search.
- `flags`: bitwise OR of CXTranslationUnit_Flags.
"""
function parse_header(header::AbstractString; index::Index=Index(), args::Vector{String}=String[],
                      includes::Vector{String}=String[], flags=CXTranslationUnit_None)
    # TODO: support parsing in-memory with CXUnsavedFile arg
    #       to _parseTranslationUnit.jj
    #unsaved_file = CXUnsavedFile()
    !isfile(header) && throw(ArgumentError("$header not found."))

    !isempty(includes) && foreach(includes) do x
        push!(args, "-I", x)
    end

    return TranslationUnit(index, header, args, flags)
end

"""
    parse_headers(headers::Vector{String}; index::Index=Index(), args::Vector{String}=String[], includes::Vector{String}=String[],
        flags = CXTranslationUnit_DetailedPreprocessingRecord | CXTranslationUnit_SkipFunctionBodies) -> Vector{TranslationUnit}
Return a [`TranslationUnit`](@ref) vector for the given headers. See also [`parse_header`](@ref).
"""
parse_headers(headers::Vector{String}; index::Index=Index(), args::Vector{String}=String[],
    includes::Vector{String}=String[], flags=CXTranslationUnit_DetailedPreprocessingRecord |
    CXTranslationUnit_SkipFunctionBodies) = [parse_header(header; index=index, args=args,
                                                          includes=includes, flags=flags) for header in unique(headers)]