text:
struct WedderburnDecomposition{B,iV,DS<:DirectSummand,Hom}
    basis::B
    invariants::iV
    Uπs::Vector{DS}
    hom::Hom
end

function WedderburnDecomposition(
    T::Type,
    G::Group,
    action::Action,
    basis_full,
    basis_half,
    S = Rational{Int};
    semisimple = false,
)
    basis = StarAlgebras.Basis{UInt32}(basis_full)
    invariants = invariant_vectors(G, action, basis)

    tbl = CharacterTable(S, G)
    ehom = CachedExtensionHomomorphism(G, action, basis_half, precompute = true)

    Uπs = let
        sa_basis = symmetry_adapted_basis(T, tbl, ehom; semisimple = semisimple)
    end

    return WedderburnDecomposition(basis, invariants, Uπs, ehom)
end

function Base.show(io::IO, wbdec::SymbolicWedderburn.WedderburnDecomposition)
    ds = direct_summands(wbdec)
    simple = all(issimple.(ds))
    dims = size.(ds, 1)
    norbs = length(invariant_vectors(wbdec))

    print(io, "Wedderburn Decomposition into $norbs orbits and $(length(ds))")
    all(simple) && print(io, " simple")
    println(io, " summands of sizes")
    return print(io, dims)
end

invariant_vectors(wbdec::WedderburnDecomposition) = wbdec.invariants
StarAlgebras.basis(wbdec::WedderburnDecomposition) = wbdec.basis
direct_summands(wbdec::WedderburnDecomposition) = wbdec.Uπs
Base.eltype(wbdec::WedderburnDecomposition) =
    eltype(eltype(direct_summands(wbdec)))

_tmps(wbdec::WedderburnDecomposition) =
    zeros.(eltype(wbdec), size.(direct_summands(wbdec)))

_fillfactor(M::AbstractMatrix) = count(!iszero, M) / length(M)
_fillfactor(M::AbstractSparseMatrix) = nnz(M) / length(M)

function diagonalize(
    M::AbstractMatrix,
    wbdec::WedderburnDecomposition,
    tmps = _tmps(wbdec),
)
    # return [degree(Uπ).*(Uπ*(M*Uπ')) for Uπ in summands(wbdec)]

    T = eltype(eltype(direct_summands(wbdec)))
    Mπs = [(d = size(Uπ, 1); zeros(T, d, d)) for Uπ in direct_summands(wbdec)]
    return diagonalize!(Mπs, M, direct_summands(wbdec), tmps)
end

function diagonalize!(
    Mπs,
    M::AbstractMatrix,
    wbdec::WedderburnDecomposition,
    tmps = _tmps(wbdec),
)
    return diagonalize!(Mπs, M, direct_summands(wbdec), tmps)
end

function diagonalize!(
    Mπs::AbstractVector{<:AbstractMatrix},
    M::AbstractMatrix,
    Uπs::AbstractVector{<:DirectSummand},
    tmps,
)
    @assert length(Mπs) == length(Uπs)

    for (π, Uπ) in enumerate(Uπs)
        imUπ = image_basis(Uπ) # Base.Matrix to allow BLAS paths below
        LinearAlgebra.mul!(tmps[π], imUπ, M)
        LinearAlgebra.mul!(Mπs[π], tmps[π], imUπ')
        zerotol!(Mπs[π], atol = eps(eltype(imUπ)) * max(size(imUπ)...))
        Mπs[π] .*= degree(Uπ)
    end
    return Mπs
end

function invariant_vectors(
    G::Group,
    act::ByPermutations,
    basis::StarAlgebras.Basis{T,I},
) where {T,I}
    tovisit = trues(length(basis))
    invariant_vs = Vector{SparseVector{Rational{Int}}}()

    ordG = order(Int, G)
    elts = collect(G)
    orbit = zeros(I, ordG)

    for i in eachindex(basis)
        if tovisit[i]
            bi = basis[i]
            Threads.@threads for j in eachindex(elts)
                orbit[j] = basis[action(act, elts[j], bi)]
            end
            tovisit[orbit] .= false
            push!(
                invariant_vs,
                sparsevec(orbit, fill(1 // ordG, ordG), length(basis)),
            )
        end
    end
    return invariant_vs
end

function _orth_proj(A::AbstractMatrix, v, QR = qr(A))
    # return A * (QR \ v)
    y = QR.Q' * v
    y[size(A, 2)+1:end] .= 0 # project to y = Q̂'x
    return QR.Q * y
end

function invariant_vectors(
    G::Group,
    act::ByLinearTransformation,
    basis,
    atol = 1e-12,
)
    hom = CachedExtensionHomomorphism(G, act, basis, precompute = true)

    dim = 0
    T = coeff_type(action(hom))
    @assert T == Float64
    subspaces = SparseMatrixCSC{T,Int64}[]
    qrs = LinearAlgebra.QRCompactWY{Float64,Matrix{Float64}}[]

    inv_vectors = SparseVector{T,Int64}[]
    inv_v = Vector{T}(undef, length(basis))

    for b in basis
        b_orth = sparsevec(decompose(b, hom)..., length(basis))
        if dim > 0
            b_orth -= sum(
                _orth_proj(V, b_orth, QR) for (V, QR) in zip(subspaces, qrs)
            )
            @debug "orthogonalizing $b:" norm(b_orth)
            norm(b_orth, 2) < atol && continue
        end

        # we found a new subspace!
        b_orth ./= norm(b_orth, 2)
        orbit_subspace = Matrix{T}(undef, length(basis), 1)
        orbit_subspace[:, 1] .= b_orth
        inv_v .= zero(T)

        # for the orbit of b_orth find
        # * `inv_v`, the invariant (average) vector, and
        # * `tmp_subspace`, the subspace spanned by the orbit

        for g in G
            k = induce(hom, g) * b_orth
            @assert isapprox(norm(k, 2), 1.0, atol = atol)
            inv_v .+= k

            residual = k - _orth_proj(orbit_subspace, k)
            (rnorm = norm(residual, 2)) < atol && continue
            orbit_subspace = hcat(orbit_subspace, Vector(residual ./ rnorm))
        end

        push!(inv_vectors, inv_v / order(Int, G)) # a copy if inv_v
        push!(subspaces, orbit_subspace) # orbit_subspace is freshly allocated every nonzero b_orth
        dim += size(orbit_subspace, 2)
        dim == length(basis) && break
        dim > length(basis) && throw("This should never occur")
        # moving qr here to save one qr decomposition
        push!(qrs, qr(orbit_subspace))
    end
    @assert dim == length(basis) "Encountered incomplete decomposition into orbits!"
    return inv_vectors
end


text:
defaultdataregistry() = Datasets.FASTAI_DATA_REGISTRY
defaultmethodregistry() = FASTAI_METHOD_REGISTRY

Datasets.listdatasources() = listdatasources(defaultdataregistry())
Datasets.finddatasets(; kwargs...) = finddatasets(defaultdataregistry(); kwargs...)
Datasets.loaddataset(args...; kwargs...) = loaddataset(defaultdataregistry(), args...; kwargs...)

findlearningmethods(args...; kwargs...) = findlearningmethods(defaultmethodregistry(), args...; kwargs...)

ext: jl


text:
# import Base: show

"""
The structure used for robot requests.
"""
mutable struct RobotRequest
  id::String
  name::String
  description::String
  status::String
end

"""
The structure returned when any robot requests are made.
"""
mutable struct RobotResponse
  robotId::String
  userId::Union{Nothing, String}
  status::String
  id::String
  name::String
  description::String
  createdTimestamp::String
  lastUpdatedTimestamp::String
  links::Dict{String, Any}
end

function show(io::IO, r::RobotResponse)
    println(io, "GraffSDK Robot:")
    println(io, " - ID: $(r.id)")
    println(io, "   - Name: $(r.name)")
    println(io, "   - Description: $(r.description)")
    println(io, "   - Status: $(r.status)")
    println(io, "   - Created: $(r.createdTimestamp)")
    println(io, "   - Last Updated: $(r.lastUpdatedTimestamp)")
end


text:
import Base:
    +,
    -,
    *,
    /,
    ^,
    %,
    !,
    >,
    <,
    >=,
    <=,
    ==,
    !=,
    ===,
    !==,
    all,
    any,
    cumprod,
    cumsum,
    falses,
    findall,
    findfirst,
    findlast,
    isnan,
    log,
    maximum,
    minimum,
    ones,
    prod,
    rand,
    randn,
    round,
    sign,
    sum,
    trues,
    zeros

import Base.Broadcast:
    BroadcastStyle,
    DefaultArrayStyle,
    Broadcasted,
    broadcasted,
    broadcastable,
    broadcast,
    result_style

import Statistics:
    mean

# enable use of functions using dot operator
struct TemporalBroadcastStyle <: BroadcastStyle end
Base.BroadcastStyle(::Type{TS{V,T}}) where {V<:VALTYPE,T<:IDXTYPE} = TemporalBroadcastStyle
BroadcastStyle(::TemporalBroadcastStyle) = TemporalBroadcastStyle()
BroadcastStyle(::TemporalBroadcastStyle, ::TemporalBroadcastStyle) = TemporalBroadcastStyle()
BroadcastStyle(::DefaultArrayStyle, ::TemporalBroadcastStyle) = TemporalBroadcastStyle()
broadcastable(x::TS) = x
broadcast(f::Function, X::TS) = TS(f.(X.values), X.index, X.fields)
broadcast(f::Function, X::TS, Y::TS) = opjoined(X, Y, f)
broadcasted(::TemporalBroadcastStyle, f, X::TS) = TS(f.(X.values), X.index, X.fields)
broadcasted(::TemporalBroadcastStyle, f, X::TS, Y::TS) = opjoined(X, Y, f)
result_style(::Type{TemporalBroadcastStyle}) = TemporalBroadcastStyle()

# passthrough functions
all(x::TS) = all(x.values)
any(x::TS) = any(x.values)
findall(x::TS) = findall(x.values)
findfirst(x::TS) = findfirst(x.values)
findlast(x::TS) = findlast(x.values)

# equivalently-shaped convenience functions
ones(x::TS) = TS(ones(size(x)), x.index, x.fields)
ones(::Type{TS}, n::Int) = TS(ones(n))
ones(::Type{TS}, dims::Tuple{Int,Int}) = TS(ones(dims))
ones(::Type{TS}, r::Int, c::Int) = TS(ones(r, c))
zeros(x::TS) = TS(zeros(size(x.values)), x.index, x.fields)
zeros(::Type{TS}, n::Int) = TS(zeros(n))
zeros(::Type{TS}, r::Int, c::Int) = TS(zeros(r, c))
zeros(::Type{TS}, dims::Tuple{Int,Int}) = TS(zeros(dims))
trues(x::TS) = TS(trues(size(x)), x.index, x.fields)
falses(x::TS) = TS(falses(size(x)), x.index, x.fields)

# random numbers
rand(::Type{TS}, n::Int=1) = TS(rand(Float64, n))
rand(::Type{TS}, r::Int, c::Int) = TS(rand(Float64, r, c))
rand(::Type{TS}, dims::Tuple{Int,Int}) = TS(rand(Float64, dims))
randn(::Type{TS}, n::Int=1) = TS(randn(Float64, n))
randn(::Type{TS}, r::Int, c::Int) = TS(randn(Float64, r, c))
randn(::Type{TS}, dims::Tuple{Int,Int}) = TS(randn(Float64, dims))

# number functions
round(x::TS; digits::Int=0) = TS(round.(x.values, digits=digits), x.index, x.fields)
round(R::Type, x::TS) = TS(round.(R, x.values), x.index, x.fields)
sum(x::TS) = sum(x.values)
sum(x::TS, dim::Int) = sum(x.values, dim)
sum(f::Function, x::TS) = sum(f, x.values)
prod(x::TS) = prod(x.values)
prod(x::TS, dim::Int) = prod(x.values, dim)
maximum(x::TS) = maximum(x.values)
maximum(x::TS, dim::Int) = maximum(x.values, dim)
minimum(x::TS) = minimum(x.values)
minimum(x::TS, dim::Int) = minimum(x.values, dim)
cumsum(x::TS; dims::Int=1) = TS(cumsum(x.values, dims=dims), x.index, x.fields)
cumprod(x::TS; dims::Int=1) = TS(cumprod(x.values, dims=dims), x.index, x.fields)
mean(x::TS) = mean(x.values)

# artithmetic operators
function opjoined(x::TS, y::TS, f::Function)
    z = [x y]
    xcols = 1:size(x,2)
    ycols = size(x,2)+1:size(x,2)+size(y,2)
    a = z[:,xcols].values
    b = z[:,ycols].values
    return TS(f.(a, b), z.index)
end

# negatiion
-(x::TS) = TS(-x.values, x.index, x.fields)
!(x::TS) = TS(.!x.values, x.index, x.fields)

# ts + ts arithmetic
+(x::TS, y::TS) = opjoined(x, y, +)
-(x::TS, y::TS) = opjoined(x, y, -)
*(x::TS, y::TS) = opjoined(x, y, *)
/(x::TS, y::TS) = opjoined(x, y, /)
^(x::TS, y::TS) = opjoined(x, y, ^)
%(x::TS, y::TS) = opjoined(x, y, %)
# ts + ts logical
==(x::TS, y::TS) = x.values == y.values && x.index == y.index && x.fields == y.fields
!=(x::TS, y::TS) = x.values != y.values || x.index != y.index || x.fields != y.fields
>(x::TS, y::TS) = opjoined(x, y, >)
<(x::TS, y::TS) = opjoined(x, y, <)
>=(x::TS, y::TS) = opjoined(x, y, >=)
<=(x::TS, y::TS) = opjoined(x, y, <=)

# ts + array
+(x::TS, y::Y) where {Y<:VALARR} = x + TS(y, x.index, x.fields)
-(x::TS, y::Y) where {Y<:VALARR} = x - TS(y, x.index, x.fields)
*(x::TS, y::Y) where {Y<:VALARR} = x * TS(y, x.index, x.fields)
/(x::TS, y::Y) where {Y<:VALARR} = x / TS(y, x.index, x.fields)
%(x::TS, y::Y) where {Y<:VALARR} = x % TS(y, x.index, x.fields)
^(x::TS, y::Y) where {Y<:VALARR} = x ^ TS(y, x.index, x.fields)

+(y::Y, x::TS) where {Y<:VALARR} = TS(y, x.index, x.fields) + x
-(y::Y, x::TS) where {Y<:VALARR} = TS(y, x.index, x.fields) - x
*(y::Y, x::TS) where {Y<:VALARR} = TS(y, x.index, x.fields) * x
/(y::Y, x::TS) where {Y<:VALARR} = TS(y, x.index, x.fields) / x
%(y::Y, x::TS) where {Y<:VALARR} = TS(y, x.index, x.fields) % x
^(y::Y, x::TS) where {Y<:VALARR} = TS(y, x.index, x.fields) ^ x

+(x::TS, y::Y) where {Y<:VALTYPE} = TS(x.values .+ y, x.index, x.fields)
-(x::TS, y::Y) where {Y<:VALTYPE} = TS(x.values .- y, x.index, x.fields)
*(x::TS, y::Y) where {Y<:VALTYPE} = TS(x.values .* y, x.index, x.fields)
/(x::TS, y::Y) where {Y<:VALTYPE} = TS(x.values ./ y, x.index, x.fields)
%(x::TS, y::Y) where {Y<:VALTYPE} = TS(x.values .% y, x.index, x.fields)
^(x::TS, y::Y) where {Y<:VALTYPE} = TS(x.values .^ y, x.index, x.fields)

+(y::Y, x::TS) where {Y<:VALTYPE} = x + y
-(y::Y, x::TS) where {Y<:VALTYPE} = TS(y .- x.values, x.index, x.fields)
*(y::Y, x::TS) where {Y<:VALTYPE} = x * y
/(y::Y, x::TS) where {Y<:VALTYPE} = TS(y ./ x.values, x.index, x.fields)
%(y::Y, x::TS) where {Y<:VALTYPE} = TS(y .% x.values, x.index, x.fields)
^(y::Y, x::TS) where {Y<:VALTYPE} = TS(y .^ x.values, x.index, x.fields)


text:
module ITensorChemistry

using Fermi
using ITensors
using Suppressor

include(joinpath("molecules", "molecule.jl"))
for file in readdir(joinpath(@__DIR__, "molecules", "molecules"); join=true)
  if endswith(file, ".jl")
    include(file)
  end
end
include("molecular_orbital_hamiltonian.jl")

export molecular_orbital_hamiltonian

end


text:
ia> [RpcApi.SimpleOrganism(0x000000000005b0cf,:(function (o::Creature.Organism,)
     function func_52(var_51::Int8=-36)
         local var_74::String = "ubD5GK5v"
         if var_51 < var_51
         end
         Creature.stepUp(o)
         Creature.eatRight(o,Int(var_51))
         Creature.stepDown(o)
         return var_51
     end
     local var_1::Int8 = 121
     Creature.eatRight(o,Int(var_1))
     for i = 1:var_1
         local var_73::Int8 = 40
         Creature.eatRight(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         local var_76::Int64 = -3451009689120369765
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.stepRight(o)
         o.mem[var_8] = var_8
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepLeft(o)
         var_6 = func_52(var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         var_8 = var_8 / var_8
         o.mem[var_8] = var_8
         var_15 = var_15 | var_15
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = var_8 + var_8
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         o.mem[var_8] = var_8
         var_6 = func_52(var_6)
         Creature.stepUp(o)
         var_15 = var_15 | var_15
         Creature.eatLeft(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
     end
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.eatLeft(o,Int(var_1))
     Creature.stepUp(o)
     var_1 = var_1 * var_1
     for i = 1:var_1
         local var_75::Int16 = -10198
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         var_75 = var_75 + var_75
         Creature.stepDown(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,2332,1,3306751,0x014fb830,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[146,148],56),RpcApi.SimpleOrganism(0x000000000005b072,:(function (o::Creature.Organism,)
     function func_71(var_69::Bool=false,var_70::Bool=false)
         return var_69
     end
     function func_62(var_60::Bool=true,var_61::String="alf6iOMU")
         Creature.stepUp(o)
         return var_60
     end
     local var_59::Bool = true
     function func_52(var_51::Int8=-36)
         Creature.eatLeft(o,Int(var_51))
         for i = 1:var_51
             local var_72::String = "oTexJKZV"
             Creature.stepUp(o)
         end
         Creature.eatLeft(o,Int(var_51))
         return var_51
     end
     local var_1::Int8 = 121
     Creature.stepDown(o)
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepLeft(o)
         var_6 = func_52(var_6)
         Creature.stepUp(o)
         Creature.stepUp(o)
         Creature.eatUp(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         var_8 = var_8 + var_8
         var_8 = var_8 / var_8
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_6))
         var_15 = var_15 | var_15
         Creature.eatLeft(o,Int(var_6))
         o.mem[var_8] = var_8
         var_6 = func_52(var_6)
         Creature.eatUp(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         var_8 = var_8 * var_8
         var_15 = Bool(abs(var_15 - var_15))
     end
     var_1 = var_1 * var_1
     for i = 1:var_1
         local var_63::Bool = false
         var_63 = var_63 & var_63
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
         var_63 = var_63 | var_63
         Creature.stepDown(o)
         Creature.stepDown(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,2332,1,2953147,0x014b3278,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[194,291],53),RpcApi.SimpleOrganism(0x000000000005b284,:(function (o::Creature.Organism,)
     function func_71(var_69::Bool=false,var_70::Bool=false)
         local var_74::Int8 = -83
         return var_69
     end
     function func_62(var_60::Bool=true,var_61::String="alf6iOMU")
         return var_60
     end
     local var_59::Bool = true
     function func_52(var_51::Int8=-36)
         Creature.eatRight(o,Int(var_51))
         Creature.eatLeft(o,Int(var_51))
         return var_51
     end
     local var_1::Int8 = 121
     Creature.stepDown(o)
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepLeft(o)
         var_6 = func_52(var_6)
         Creature.stepUp(o)
         Creature.stepUp(o)
         Creature.eatUp(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         var_8 = var_8 + var_8
         var_8 = var_8 / var_8
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_6))
         var_15 = var_15 | var_15
         Creature.eatLeft(o,Int(var_6))
         o.mem[var_8] = var_8
         var_6 = func_52(var_6)
         Creature.eatUp(o,Int(var_6))
         o.mem[var_8] = var_8
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         var_8 = var_8 * var_8
         var_15 = Bool(abs(var_15 - var_15))
     end
     var_1 = var_1 * var_1
     for i = 1:var_1
         local var_63::Bool = false
         var_63 = var_63 & var_63
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
         var_63 = var_63 | var_63
         Creature.stepDown(o)
         Creature.stepDown(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,2332,1,2262614,0x014bb588,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[97,146],51),RpcApi.SimpleOrganism(0x000000000005b2a5,:(function (o::Creature.Organism,)
     function func_78(var_77::Int64=6096170113460993338)
         return var_77
     end
     local var_75::Int8 = -52
     function func_52(var_51::Int8=-36)
         local var_76::Bool = false
         Creature.eatUp(o,Int(var_51))
         for i = 1:var_51
         end
         Creature.eatRight(o,Int(var_51))
         Creature.eatRight(o,Int(var_51))
         var_76 = var_76 | var_76
         Creature.stepDown(o)
         Creature.eatLeft(o,Int(var_51))
         Creature.eatDown(o,Int(var_51))
         return var_51
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         local var_73::Int8 = 40
         var_73 = var_73 * var_73
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_73))
         Creature.eatUp(o,Int(var_73))
         Creature.eatRight(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.stepLeft(o)
         Creature.eatDown(o,Int(var_73))
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_68 = var_68 * var_68
         Creature.eatUp(o,Int(var_6))
         Creature.stepLeft(o)
         var_6 = func_52(var_6)
         Creature.eatUp(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         o.mem[var_8] = var_8
         Creature.stepLeft(o)
         Creature.eatUp(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = var_8 + var_8
         var_8 = var_8 / var_8
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepDown(o)
         o.mem[var_8] = var_8
         Creature.eatLeft(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
         Creature.eatDown(o,Int(var_6))
     end
     var_1 = func_52(var_1)
     Creature.stepDown(o)
     var_1 = var_75 - var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         var_63 = Bool(abs(var_63 - var_63))
         Creature.stepUp(o)
         Creature.stepDown(o)
     end
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],7,2332,1,2253280,0x01525d58,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[359,335],65),RpcApi.SimpleOrganism(0x000000000005ae6a,:(function (o::Creature.Organism,)
     local var_75::Int8 = -52
     function func_52(var_51::Int8=-36)
         local var_76::Bool = false
         Creature.eatRight(o,Int(var_51))
         var_76 = var_76 | var_76
         Creature.stepLeft(o)
         Creature.stepUp(o)
         return var_51
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         local var_73::Int8 = 40
         var_73 = var_73 * var_73
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.stepLeft(o)
         Creature.eatDown(o,Int(var_73))
     end
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatUp(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         Creature.stepLeft(o)
         var_6 = func_52(var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_6))
         o.mem[var_8] = var_8
         Creature.eatUp(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = var_8 + var_8
         var_8 = var_8 / var_8
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         o.mem[var_8] = var_8
         Creature.eatLeft(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
         Creature.eatDown(o,Int(var_6))
     end
     var_1 = func_52(var_1)
     Creature.stepDown(o)
     var_1 = var_1 * var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepUp(o)
         Creature.stepDown(o)
     end
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],7,2332,1,1944625,0x01503b40,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[98,264],56),RpcApi.SimpleOrganism(0x000000000005b29d,:(function (o::Creature.Organism,)
     function func_78(var_77::Int64=6096170113460993338)
         return var_77
     end
     local var_75::Int8 = -52
     function func_52(var_51::Int8=-36)
         local var_76::Bool = false
         Creature.eatUp(o,Int(var_51))
         for i = 1:var_51
         end
         Creature.eatRight(o,Int(var_51))
         Creature.eatRight(o,Int(var_51))
         var_76 = var_76 | var_76
         Creature.stepDown(o)
         Creature.eatLeft(o,Int(var_51))
         Creature.eatDown(o,Int(var_51))
         return var_51
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         local var_73::Int8 = 40
         var_73 = var_73 * var_73
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_73))
         Creature.eatUp(o,Int(var_73))
         Creature.eatRight(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.eatLeft(o,Int(var_73))
         Creature.stepLeft(o)
         Creature.eatDown(o,Int(var_73))
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_68 = var_68 * var_68
         Creature.eatUp(o,Int(var_6))
         Creature.stepLeft(o)
         var_6 = func_52(var_6)
         Creature.eatUp(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         o.mem[var_8] = var_8
         Creature.eatUp(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = var_8 + var_8
         var_8 = var_8 / var_8
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepLeft(o)
         o.mem[var_8] = var_8
         Creature.eatLeft(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
         Creature.eatDown(o,Int(var_6))
     end
     var_1 = func_52(var_1)
     Creature.stepDown(o)
     var_1 = var_75 - var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         var_63 = Bool(abs(var_63 - var_63))
         Creature.stepUp(o)
         Creature.stepDown(o)
     end
     Creature.eatLeft(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],7,2332,1,1895848,0x015231a8,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[273,335],64),RpcApi.SimpleOrganism(0x000000000005b214,:(function (o::Creature.Organism,)
     function func_66(var_65::Bool=true)
         local var_67::Int16 = -27495
         Creature.stepDown(o)
         Creature.stepRight(o)
         return var_65
     end
     function func_50(var_48::Int64=8246173392054715175,var_49::Bool=false)
         return var_48
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         var_18 = var_18 / var_18
         Creature.stepUp(o)
         return var_17
     end
     local var_1::Int8 = 121
     if var_1 <= var_1
     end
     for i = 1:var_1
         local var_61::Int16 = -15736
     end
     for i = 1:var_1
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepRight(o)
     end
     Creature.stepRight(o)
     for i = 1:var_1
         local var_52::String = "3Zf0ESw4"
         local var_37::Int8 = 66
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.eatLeft(o,Int(var_6))
         var_52 = var_52 * var_52
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatLeft(o,Int(var_37))
         Creature.eatDown(o,Int(var_37))
         Creature.eatLeft(o,Int(var_6))
         var_9 = var_9 | var_9
         var_8 = var_8 / var_8
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         var_52 = var_52[1:if length(var_52) > length(var_52) > 0
                         0
                     else
                         div(length(var_52),length(var_52))
                     end]
         Creature.eatDown(o,Int(var_37))
         Creature.eatUp(o,Int(var_6))
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
         Creature.stepUp(o)
         Creature.stepRight(o)
     end
     for i = 1:var_1
         local var_62::Bool = false
     end
     Creature.stepRight(o)
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],13,3200,1,1784773,0x01516d10,Dict(23832=>23832),[398,213],57)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SQLite_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SQLite")
JLLWrappers.@generate_main_file("SQLite", UUID("76ed43ae-9a5d-5a62-8c75-30186b810ce8"))
end  # module SQLite_jll


text:
#=

Preparation of the data

Ignacio Quintero Mächler

t(-_-t)

Created 16 03 2020

=#

"""
    prepare_data(cov_mod    ::NTuple{M,String},
                 tv         ::Dict{Int64,Array{Float64,1}},
                 x          ::Array{Float64,1},
                 y          ::Array{Float64,N},
                 ed         ::Array{Int64,2},
                 el         ::Array{Float64,1},
                 ρ          ::Array{Float64,1},
                 h          ::Int64,
                 constraints::NTuple{O,String},
                 mvpars     ::NTuple{P,String})

Prepare data for **EGeoHiSSE** likelihood calculations.
"""
function prepare_data(cov_mod    ::NTuple{M,String},
                      tv         ::Dict{Int64,Array{Float64,1}},
                      x          ::Array{Float64,1},
                      y          ::Array{Float64,N},
                      ed         ::Array{Int64,2},
                      el         ::Array{Float64,1},
                      ρ          ::Array{Float64,1},
                      h          ::Int64,
                      constraints::NTuple{O,String},
                      mvpars     ::NTuple{P,String}) where {N,M,O,P}

  # k areas
  k = length(tv[1])::Int64

  # number of covariates
  ny = size(y,2)

 # number of states
  ns = (2^k - 1)*h

  # number of tips
  ntip = length(tv)

  # add root of length 0
  ed = cat(ed, [2*ntip ntip + 1], dims = 1)
  push!(el, 0.0)

  # number of edges
  ned = size(ed,1)

  # make z(t) approximation from discrete data `af!()`
  af! = make_af(x, y, Val(ny))

  # define model
  model = define_mod(cov_mod, k, h, ny)

  # make dictionary with relevant parameters
  pardic = build_par_names(k, h, ny, model)

  # get number of parameters
  npars = length(pardic)

  # find hidden factors for hidden states
  phid = Int64[]
  if h > 1
    re   = Regex(".*_[1-"*string(h-1)*"]\$")
    for (k,v) in pardic
      occursin(re, k) && push!(phid, v)
    end
    sort!(phid)
  end

  # create factor parameter vector
  fp = zeros(npars)

  # generate initial parameter values
  p  = fill(1e-1, npars)
  βs = h*(k^2 + 2k + h) + 1
  δ  = Float64(ntip-1)/sum(el)
  p[βs:end]             .= 0.0                  # set βs
  p[1:(k+1)*h]          .= δ + rand()*δ         # set λs
  p[(k+1)*h+1:h*(2k+1)] .= p[1] - δ             # set μs

  # parameter update
  pupd = 1:npars

  # parameters constraints and fixed to 0
  dcp, dcfp, zp, zfp =
    set_constraints(constraints, pardic, k, h, ny, model)

  # set multivariate sampled parameters
  mvps0 = set_multivariate(mvpars, pardic)

  # remove hidden factors from being updated from `p`
  pupd = setdiff(pupd, phid)

  # force pars in zerp to 0
  for i in zp
    p[i] = 0.0
  end

  # remove constraints and fixed to zero parameters from being updated
  setdiff!(pupd, values(dcp), zp)
  setdiff!(phid, values(dcfp), zfp)

  # check if there are hidden factors forced to 0 that also have a forced equality
  for (k,v) in dcfp
    if in(k, zfp) || in(v, zfp)
      filter!(x -> x ≠ k, phid)
      filter!(x -> x ≠ v, phid)
    end
  end

  # divide between non-negative and negative values
  nnps = filter(x -> βs >  x, pupd)
  nps  = filter(x -> βs <= x, pupd)

  # divide multivariate updates into hidden and non-hidden updates
  mvps  = Array{Int64,1}[]
  nngps = Array{Bool,1}[]
  mvhfs = Array{Int64,1}[]
  for (p1, p2) in mvps0

    p1nn = in(p1, nnps)
    p1n  = in(p1, nps)
    p2nn = in(p2, nnps)
    p2n  = in(p2, nps)

    if p1nn && p2nn
      setdiff!(nnps, p1, p2)
      push!(mvps, [p1,p2])
      push!(nngps, [p1nn,p2nn])
    end

    if p1n && p2nn
      setdiff!(nps,  p1)
      setdiff!(nnps, p2)
      push!(mvps, [p1,p2])
      push!(nngps, [p1nn, p2nn])
    end

    if p1nn && p2n
      setdiff!(nnps, p1)
      setdiff!(nps,  p2)
      push!(mvps, [p1,p2])
      push!(nngps, [p1nn, p2nn])
    end

    if p1n && p2n
      setdiff!(nps, p1, p2)
      push!(mvps, [p1,p2])
      push!(nngps, [p1nn, p2nn])
    end

    if in(p1, phid) && in(p2, phid)
      setdiff!(phid, p1, p2)
      push!(mvhfs, [p1,p2])
    end
  end

  # make hidden factors assigning
  assign_hidfacs! = make_assign_hidfacs(Val(k), Val(h), Val(ny), Val(model))

  # force same parameter values for constraints
  for wp in keys(dcp)
    while haskey(dcp, wp)
      tp = dcp[wp]
      p[tp] = p[wp]
      wp = tp
    end
  end

  for wp in keys(dcfp)
    while haskey(dcfp, wp)
      tp = dcfp[wp]
      fp[tp] = fp[wp]
      wp = tp
    end
  end

  # assign hidden factors
  assign_hidfacs!(p, fp)

  # extinction at time 0 with sampling fraction `ρ_i`
  if isone(lastindex(ρ))
    E0 = fill(1.0 - ρ[1], ns)
  elseif lastindex(ρ) == ns
    E0 = ones(Float64, ns)
    for i in Base.OneTo(ns)
      E0[i] -= ρ[i]
    end
  else
    @error "Length of sampling fraction vector $(lastindex(ρ)) does not match length of the number of states $ns"
  end

  # create geographic & hidden states
  S = create_states(k, h)

  # get absolute times of branches as related to z(t)
  abts = abs_time_branches(el, ed, ntip)

  # get branching times
  bts = map(x -> round(x; digits = 9), abts[:,1])
  bts = unique!(bts)

  sort!(bts)

  # make trios
  trios = maketriads(ed, rev = true)

  # preallocate tip likelihoods
  X = [zeros(ns) for i in Base.OneTo(ned)]

  child = ed[:,2]
  wtp   = findall(child .<= ntip)

  # assign states to terminal branches
  for wi in wtp
    wig = Set(findall(map(x -> isone(x), tv[child[wi]])))
    X[wi][findall(map(x -> isequal(x.g, wig), S))] .= 1.0
  end

  return X, p, fp, trios, ns, ned, pupd, phid, nnps, nps, mvps, nngps, mvhfs,
    dcp, dcfp, pardic, k, h, ny, model, af!, assign_hidfacs!, abts, bts, E0

end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
## File contains the fitting funcs

using JuMP
using Gurobi
using Distributions

include("BertrandNash.jl")  #file generates data

function setUpFitting()
	m = Model(solver=GurobiSolver(OutputFlag=false))
	@defVar(m, vthetas1[1:3])
	@defVar(m, vthetas2[1:3])
	# @addConstraint(m, vthetas1[2])
	# @addConstraint(m, vthetas2[1])
	@defVar(m, vintercept1 >= 0)
	@defVar(m, vintercept2 >= 0)
	#normalization choice
	@addConstraint(m, vthetas1[1] <= -1.)
	@addConstraint(m, vthetas2[2] <= -1.)
	return m, vthetas1, vthetas2, vintercept1, vintercept2
end

#adds the constraints to the modl for the residual for this observation
function addResids(m, thetas1, thetas2, intercept1, intercept2, prices_t, GDP_t, pmax)
	@defVar(m, y[1:2] >= 0)
	@defVar(m, resid >= 0)
	@defVar(m, d[1:2] >=0)  #doesn't seem to work if we insist these are positive...
	@addConstraint(m, d[1] == prices_t[1] * thetas1[1] + prices_t[2] * thetas1[2] + GDP_t * thetas1[3] + intercept1 )
	@addConstraint(m, d[2] == prices_t[1] * thetas2[1] + prices_t[2] * thetas2[2] + GDP_t * thetas2[3] + intercept2 )
	@addConstraint(m, y[1] >= d[1] + prices_t[1] * thetas1[1])
	@addConstraint(m, y[2] >= d[2] + prices_t[2] * thetas2[2])
	@addConstraint(m, pmax * (y[1] + y[2]) - prices_t[1] * d[1] - prices_t[2] * d[2] -
					prices_t[1]^2 * thetas1[1] - prices_t[2]^2 * thetas2[2] <= resid )
	return resid
end

#calc an out of sample residual.  may be a better way to do this.
function calcResid(fitTheta1, fitTheta2, fitIntercept1, fitIntercept2, prices, GDP, pmax)
	m, vthetas1, vthetas2, vintercept1, vintercept2 = setUpFitting()
	for i = 1:3
		@addConstraint(m, vthetas1[i] == fitTheta1[i])
		@addConstraint(m, vthetas2[i] == fitTheta2[i])
	end
	@addConstraint(m, vintercept1 == fitIntercept1)
	@addConstraint(m, vintercept2 == fitIntercept2)
	resid = addResids(m, vthetas1, vthetas2, vintercept1, vintercept2, prices, GDP, pmax)
	@setObjective(m, Min, resid)
	status = solve(m);
	if status != :Optimal
		println(prices, status)
		error(status)
	end

	return getObjectiveValue(m)
end

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs


text:
# ---
# title: 553. Optimal Division
# id: problem553
# author: Indigo
# date: 2021-06-23
# difficulty: Medium
# categories: Math, String
# link: <https://leetcode.com/problems/optimal-division/description/>
# hidden: true
# ---
#
# Given a list of **positive integers** , the adjacent integers will perform the
# float division. For example, [2,3,4] -> 2 / 3 / 4.
#
# However, you can add any number of parenthesis at any position to change the
# priority of operations. You should find out how to add parenthesis to get the
# **maximum** result, and return the corresponding expression in string format.
# **Your expression should NOT contain redundant parenthesis.**
#
# **Example:**
#
#
#
#     Input: [1000,100,10,2]
#     Output: "1000/(100/10/2)"
#     Explanation:
#     1000/(100/10/2) = 1000/((100/10)/2) = 200
#     However, the bold parenthesis in "1000/( **(** 100/10 **)** /2)" are redundant,
#     since they don't influence the operation priority. So you should return "1000/(100/10/2)".
#
#     Other cases:
#     1000/(100/10)/2 = 50
#     1000/(100/(10/2)) = 50
#     1000/100/10/2 = 0.5
#     1000/100/(10/2) = 2
#
#
# **Note:**
#
#   1. The length of the input array is [1, 10].
#   2. Elements in the given array will be in range [2, 1000].
#   3. There is only one optimal division for each test case.
#
#
## @lc code=start
using LeetCode

function optimal_division(nums::Vector{Int})
    (length(nums) == 1) && return string(nums[1])
    (length(nums) == 2) && return string(nums[1]) * "/" * string(nums[1])
    res = "$(nums[1])/($(nums[2])"
    for i in 3:length(nums)
        res *= "/$(nums[i])"
    end
    return res * ")"
end
## @lc code=end


text:
broadcast(::typeof(identity), t::Table) = copy(t) # output shouldn't alias input
map(::typeof(identity), t::Table) = copy(t)

broadcast(::typeof(merge), t::Table) = copy(t) # output shouldn't alias input
map(::typeof(merge), t::Table) = copy(t)

function broadcast(::typeof(merge), t1::Table, t2::Table)
    map(merge, t1, t2)
end

function map(::typeof(merge), t1::Table, t2::Table)
    # TODO merge can overwrite a column. Need to invalidate the appropriate acceleration
    # indexes in this case
    return Table(merge(columns(t1), columns(t2)), (accelerators(t1)..., accelerators(t2)...))
end

# TODO broadcast/map Project (and stop doing autovectorization?)
# TODO broadcast/map Rename (and stop doing autovectorization?)

# TODO should we intercept all `map`s and `broadcast`s that map rows to rows (NamedTuples)?
# (we would need to rely on inference...)


text:
using CImPlot_jll

using CImGui

import CImGui:
    # Vector primitives:
    ImVec2, ImVec4,
    # Enums
    ImGuiMouseButton, ImGuiKeyModFlags, ImGuiCond, ImGuiDragDropFlags,
    # Primitive type aliases; uncomment after CImGui update
    #=ImS8,=# ImU8, ImS16, ImU16, ImS32, ImU32, ImS64, ImU64,
    ImTextureID,
    ImDrawList,
    ImGuiContext

#Temporary patch; CImGui.jl v1.79.0 aliases ImS8 incorrectly
const ImS8 = Int8

const IMPLOT_AUTO = Cint(-1)
const IMPLOT_AUTO_COL = ImVec4(0,0,0,-1)
export IMPLOT_AUTO, IMPLOT_AUTO_COL


text:
export GaussianNoise

"""
    GaussianNoise(σ = 1.0)

Returns a Gaussian noise with standard deviation `σ`.
"""
struct GaussianNoise <: Sequence
  σ
end

GaussianNoise() = GaussianNoise(1.0)

sample(seq::GaussianNoise, t::AbstractVector) = seq.σ * randn(length(t))


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

using Markdown
using InteractiveUtils

# ╔═╡ 9e4e93ea-8793-11eb-100e-4d3853237dde
begin
  using Pkg; Pkg.activate("MLJ_env", shared=true)
  using PlutoUI

    using CategoricalArrays
    import MLJLinearModels
    import DataFrames
    import CSV
    import DecisionTree
    using MLJ
    import MLJLinearModels
    import MultivariateStats
    import Distributions
    using Plots
    using Test
end

# ╔═╡ 1282be54-8793-11eb-3575-af6b5b5cc6eb
md"""
## A workshop introducing the machine learning toolbox MLJ - Part4 - Tuning Hyper-parameters

  - source:[MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
  - based on [Machine Learning in Julia using MLJ, JuliaCon2020](https://github.com/ablaom/MachineLearningInJulia2020)
"""

# ╔═╡ 3b3dc082-8793-11eb-28d4-a9ec7ae95230
html"""

"""

# ╔═╡ 3b1f0f8c-8793-11eb-01ae-43b0e6864086
md"""
##### TOC
  - [Naive tuning of a single Parameter](#naive-tuning)
  - [Self-Tuning Models](#self-tuning)
  - [Unbounded ranges and sampling](#unbounded-ranges)
  - [Resources or this part](#resources)
  - [Exercises for Part 4](#exercises)
"""

# ╔═╡ 2235d2c8-8793-11eb-18eb-ef8d66e5bb8e
html"""

  back to TOC

"""

# ╔═╡ dc7748d4-8792-11eb-28cf-635c2c924dc3
md"""
### Naive tuning of a single parameter

The most naive way to tune a single hyper-parameter is to use `learning_curve`, which we already saw in Part 2. Let's see this in the Horse Colic classification problem, in a case where the parameter to be tuned is *nested* (because the model is a pipeline):
"""

# ╔═╡ 7fb9f9ae-8793-11eb-161b-4131e1621b0a
begin
  const DIR = "."
    hfile = CSV.File(joinpath(DIR, "data", "horse.csv"));
    horse = DataFrames.DataFrame(hfile); # convert to data frame without copying columns
end

# ╔═╡ c3500cfa-8793-11eb-1f7b-6d23b07dd9a7
begin
  coerce!(horse, autotype(horse));
  coerce!(horse, Count => Continuous);
  coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
  schema(horse)
end

# ╔═╡ c331a116-8793-11eb-2de4-3130572c91cc
begin
  yₕ, Xₕ = unpack(horse, ==(:outcome), name -> true);
  scitype(yₕ)
end

# ╔═╡ c3147f0a-8793-11eb-34d1-35ad468cc0fa
@load LogisticClassifier pkg=MLJLinearModels

# ╔═╡ c2fd2da0-8793-11eb-16ce-4bf069a83cf0
begin
  modelₕ = @pipeline Standardizer ContinuousEncoder        MLJLinearModels.LogisticClassifier

  machₕ = machine(modelₕ, Xₕ, yₕ);
end

# ╔═╡ 2034a066-8794-11eb-21c1-af64ee1bdac8
rₕ = range(modelₕ, :(logistic_classifier.lambda), lower = 1e-2, upper=100,
  scale=:log10)

# ╔═╡ 34e486fc-8794-11eb-05a2-d19cd6695a73
md"""
If you're curious, you can see what `lambda` values this range will generate for a given resolution:
"""

# ╔═╡ 20190554-8794-11eb-21ad-5bb2a9c0016f
iterator(rₕ, 5)

# ╔═╡ 20015968-8794-11eb-37de-7f7e9074c251
begin
  _, _, lambdas, losses = learning_curve(machₕ, range=rₕ,
    resampling=CV(nfolds=6),
        resolution=30, # default
        measure=cross_entropy)

  plot(lambdas, losses, xscale=:log10,
    xlabel="lambda",
    ylabel="cross entropy using 6-fold CV")
end

# ╔═╡ 1fe99152-8794-11eb-35b7-ffc47f8024ea
best_lambda = lambdas[argmin(losses)]

# ╔═╡ e7e6beee-87c3-11eb-3ad0-41c7e2c612c0
html"""

  back to TOC

"""

# ╔═╡ 1fd19e44-8794-11eb-0e65-c935ea2d9867
md"""
### Self tuning models

A more sophisticated way to view hyper-parameter tuning (inspired by
MLR) is as a model *wrapper*. The wrapped model is a new model in
its own right and when you fit it, it tunes specified
hyper-parameters of the model being wrapped, before training on all
supplied data. Calling `predict` on the wrapped model is like
calling `predict` on the original model, but with the
hyper-parameters already optimized.

In other words, we can think of the wrapped model as a "self-tuning"
version of the original.

We now create a self-tuning version of the pipeline above, adding a
parameter from the `ContinuousEncoder` to the parameters we want
optimized.

First, let's choose a tuning strategy (from [these
options](https://github.com/alan-turing-institute/MLJTuning.jl#what-is-provided-here)). MLJ
supports ordinary `Grid` search (query `?Grid` for
details). However, as the utility of `Grid` search is limited to a
small number of parameters, and as `Grid` searches are demonstrated
elsewhere (see the [resources below](#resources-for-part-4)) we'll
demonstrate `RandomSearch` here:
"""

# ╔═╡ 1fb75a0a-8794-11eb-1157-f38619aaffac
tuning = RandomSearch(rng=123)

# ╔═╡ 1f99d98c-8794-11eb-3365-09608fee79cd
md"""
In this strategy each parameter is sampled according to a pre-specified prior distribution that is fit to the one-dimensional range object constructed using `range` as before. While one has a lot of control over the specification of the priors (run `?RandomSearch` for details) we'll let the algorithm generate these priors automatically.
"""

# ╔═╡ 9608ba16-8794-11eb-12ee-a905643b6de9
html"""

  back to TOC

"""

# ╔═╡ 843525f4-8794-11eb-262a-bd9839059599
md"""
#### Unbounded ranges and sampling

In MLJ a range does not have to be bounded. Furthermore, in `RandomSearch` a A positive unbounded range is specified will be sampled using a `Gamma` distribution, by default:
"""

# ╔═╡ 83fe4d2c-8794-11eb-1395-c5908bda104c
r₂ = range(modelₕ, :(logistic_classifier.lambda), lower=0, origin=6, unit=5)

# ╔═╡ 83e467cc-8794-11eb-0ed6-6900e640b3be
md"""
Let's see what sampling using a Gamma distribution is going to mean for this range:
"""

# ╔═╡ 19325af0-8795-11eb-182b-251ef4508684
begin
  # import Distributions
  sampler_r = sampler(r₂, Distributions.Gamma)
  histogram(rand(sampler_r, 10000), nbins=50)
end

# ╔═╡ 19125886-8795-11eb-12f2-e1f6258d0aeb
md"""
The second parameter we'll add to this is *nominal* (finite) and, by default, will be sampled uniformly. Since it is nominal, we specify `values` instead of `upper` and `lower` bounds:
"""

# ╔═╡ 18f8de56-8795-11eb-0cfd-d73c64d968bd
r₃ = range(modelₕ, :(continuous_encoder.one_hot_ordered_factors),
  values = [true, false])

# ╔═╡ 18e01c22-8795-11eb-26ab-05ecf1c7ff90
md"""
##### The tuning wrapper

Now for the wrapper, which is an instance of `TunedModel`:
"""

# ╔═╡ 18c7df90-8795-11eb-21f6-67e53e00a527
tuned_model = TunedModel(model=modelₕ, ranges=[r₂, r₃],
  resampling=CV(nfolds=6), measures=cross_entropy,
  tuning=tuning, n=15)

# ╔═╡ 18aed248-8795-11eb-054b-1558b74bc6ea
md"""
We can apply the `fit!/predict` work-flow to `tuned_model` just as for any other model:
"""

# ╔═╡ 83c9e1b8-8794-11eb-1133-77f9475b62b8
begin
  tuned_mach = machine(tuned_model, Xₕ, yₕ);
  fit!(tuned_mach);
  predict(tuned_mach, rows=1:3)
end

# ╔═╡ f047f48a-8795-11eb-1f0f-b5e2fe9d7a2a
md"""
The outcomes of the tuning can be inspected from a detailed report. For example, we have:
"""

# ╔═╡ f04773ac-8795-11eb-3d99-358e3fb52fcc
begin
  rep_tuned_mach = report(tuned_mach);
  rep_tuned_mach.best_model
end

# ╔═╡ 2cad2636-8796-11eb-26e5-3dea19ddc99f
md"""
In the special case of two-parameters, you can also plot the results:
"""

# ╔═╡ 2c8a3108-8796-11eb-2e24-1f46d8653318
plot(tuned_mach)

# ╔═╡ 2c6fd7ea-8796-11eb-3fdf-8f29bb83cc37
md"""
Finally, let's compare cross-validation estimate of the performance of the self-tuning model with that of the original model (an example  of [*nested resampling*](https://mlr3book.mlr-org.com/nested-resampling.html) here):
"""

# ╔═╡ 2c549336-8796-11eb-0ac9-db88b339e91a
errₕ = evaluate!(machₕ, resampling=CV(nfolds=3), measure=cross_entropy)

# ╔═╡ 2c3af49e-8796-11eb-3fe8-43ba51426411
tuned_err = evaluate!(tuned_mach, resampling=CV(nfolds=3), measure=cross_entropy)

# ╔═╡ 57a0fe58-8796-11eb-0fb7-332dfa29bcf5
html"""

  back to TOC

"""

# ╔═╡ 5516fd22-8796-11eb-20ba-ef554955d071
md"""
#### Resources for Part 4

- From the MLJ manual:
   - [Learning Curves](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_curves/)
   - [Tuning Models](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/)
- The [MLJTuning repo](https://github.com/alan-turing-institute/MLJTuning.jl#who-is-this-repo-for) - mostly for developers

- From Data Science Tutorials:
    - [Tuning a model](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/model-tuning/)
    - [Crabs with XGBoost](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/crabs-xgb/) `Grid` tuning in stages for a tree-boosting model with many parameters
    - [Boston with LightGBM](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/boston-lgbm/) -  `Grid` tuning for another popular tree-booster
    - [Boston with Flux](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/boston-flux/) - optimizing batch size in a simple neural network regressor
- [UCI Horse Colic Data Set](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)

"""

# ╔═╡ 54fc4464-8796-11eb-324b-0195abc8c4bc
html"""

  back to TOC

"""

# ╔═╡ 54e17b3e-8796-11eb-3b20-45a682c2065b
md"""
#### Exercises for Part 4

##### Exercise 1

This exercise continues our analysis of the King County House price prediction problem:
"""

# ╔═╡ 54c55cce-8796-11eb-3774-058c344a886a
begin
  hofile = CSV.File(joinpath(DIR, "data", "house.csv"));
  house = DataFrames.DataFrame(hofile);

  coerce!(house, autotype(hofile));
  coerce!(house, Count => Continuous, :zipcode => Multiclass);

  yₓ, Xₓ = unpack(house, ==(:price), name -> true, rng=123);
end

# ╔═╡ f04752ac-8795-11eb-03a3-a56cd548eca8
md"""
Your task will be to tune the following pipeline regression model:
"""

# ╔═╡ 1ae849ac-8797-11eb-2376-4189345388ff
@load(EvoTreeRegressor)

# ╔═╡ 1accb14c-8797-11eb-1c71-81ad55c19b1d
tree_booster = EvoTrees.EvoTreeRegressor(nrounds = 70);

# ╔═╡ 1ab6cef4-8797-11eb-2f81-6140be1f3b0d
modelₓ = @pipeline ContinuousEncoder tree_booster

# ╔═╡ 1a9b44e0-8797-11eb-3e92-2f8dbd5282ae
md"""
(a) Construct a bounded range `rₓ` for the `evo_tree_booster` parameter `max_depth`, varying between 1 and 12.

 (b) For the `nbins` parameter of the `EvoTreeRegressor`, define the range
"""

# ╔═╡ 1a7e108c-8797-11eb-11f8-f7be28d15ab9
# (a)
rₓ = range(modelₓ, :(evo_tree_regressor.max_depth),
  lower = 1.0, upper= 12.0)

# ╔═╡ 6c41590e-87c6-11eb-17e2-bd5698ae552f
# (b)
rₒ = range(modelₓ,
  :(evo_tree_regressor.nbins),
  lower = 2.5,
  upper= 7.5, scale=x->2^round(Int, x))

# ╔═╡ 43e21928-8797-11eb-0e87-cd5790ab59b1
md"""
(c) Optimize `model` over these parameter ranges `rₓ` and `r2`
using a random search with uniform priors (the default). Use
`Holdout()` resampling, and implement your search by first
constructing a "self-tuning" wrap of `model`, as described
above. Make `mae` (mean absolute error) the loss function that you
optimize, and search over a total of 40 combinations of
hyper-parameters.  If you have time, plot the results of your
search. Feel free to use all available data.

(d) Evaluate the best model found in the search using 3-fold
cross-validation and compare with that of the self-tuning model
(which is different!). Setting data hygiene concerns aside, feel
free to use all available data.
"""

# ╔═╡ 43cc0f16-8797-11eb-1cca-874f98f40302
tuned_modelₓ = TunedModel(model=modelₓ,
                         ranges=[rₓ, rₒ],
                         resampling=Holdout(),
                         measures=mae,
                         tuning=RandomSearch(rng=123),
                         n=40)

# ╔═╡ d60b7a5a-87c5-11eb-0685-5100caf67827
tuned_machₓ = machine(tuned_modelₓ, Xₓ, yₓ) |> fit!

# ╔═╡ d5ea675c-87c5-11eb-1bf4-d1ea71c95392
plot(tuned_machₓ)

# ╔═╡ ab758052-8797-11eb-1193-b1f04d8f328c
# (d)
begin
  best_modelₓ = report(tuned_machₓ).best_model;
  best_machₓ = machine(best_modelₓ, Xₓ, yₓ);
  best_errₓ = evaluate!(best_machₓ, resampling=CV(nfolds=3), measure=mae)
end

# ╔═╡ 0388a272-87c7-11eb-21e9-9190a2f07a1c
tuned_errₓ = evaluate!(tuned_machₓ, resampling=CV(nfolds=3), measure=mae)

# ╔═╡ Cell order:
# ╟─1282be54-8793-11eb-3575-af6b5b5cc6eb
# ╠═9e4e93ea-8793-11eb-100e-4d3853237dde
# ╟─3b3dc082-8793-11eb-28d4-a9ec7ae95230
# ╟─3b1f0f8c-8793-11eb-01ae-43b0e6864086
# ╟─2235d2c8-8793-11eb-18eb-ef8d66e5bb8e
# ╟─dc7748d4-8792-11eb-28cf-635c2c924dc3
# ╠═7fb9f9ae-8793-11eb-161b-4131e1621b0a
# ╠═c3500cfa-8793-11eb-1f7b-6d23b07dd9a7
# ╠═c331a116-8793-11eb-2de4-3130572c91cc
# ╠═c3147f0a-8793-11eb-34d1-35ad468cc0fa
# ╠═c2fd2da0-8793-11eb-16ce-4bf069a83cf0
# ╠═2034a066-8794-11eb-21c1-af64ee1bdac8
# ╟─34e486fc-8794-11eb-05a2-d19cd6695a73
# ╠═20190554-8794-11eb-21ad-5bb2a9c0016f
# ╠═20015968-8794-11eb-37de-7f7e9074c251
# ╠═1fe99152-8794-11eb-35b7-ffc47f8024ea
# ╟─e7e6beee-87c3-11eb-3ad0-41c7e2c612c0
# ╠═1fd19e44-8794-11eb-0e65-c935ea2d9867
# ╠═1fb75a0a-8794-11eb-1157-f38619aaffac
# ╟─1f99d98c-8794-11eb-3365-09608fee79cd
# ╟─9608ba16-8794-11eb-12ee-a905643b6de9
# ╟─843525f4-8794-11eb-262a-bd9839059599
# ╠═83fe4d2c-8794-11eb-1395-c5908bda104c
# ╟─83e467cc-8794-11eb-0ed6-6900e640b3be
# ╠═19325af0-8795-11eb-182b-251ef4508684
# ╟─19125886-8795-11eb-12f2-e1f6258d0aeb
# ╠═18f8de56-8795-11eb-0cfd-d73c64d968bd
# ╟─18e01c22-8795-11eb-26ab-05ecf1c7ff90
# ╠═18c7df90-8795-11eb-21f6-67e53e00a527
# ╟─18aed248-8795-11eb-054b-1558b74bc6ea
# ╠═83c9e1b8-8794-11eb-1133-77f9475b62b8
# ╟─f047f48a-8795-11eb-1f0f-b5e2fe9d7a2a
# ╠═f04773ac-8795-11eb-3d99-358e3fb52fcc
# ╟─2cad2636-8796-11eb-26e5-3dea19ddc99f
# ╠═2c8a3108-8796-11eb-2e24-1f46d8653318
# ╟─2c6fd7ea-8796-11eb-3fdf-8f29bb83cc37
# ╠═2c549336-8796-11eb-0ac9-db88b339e91a
# ╠═2c3af49e-8796-11eb-3fe8-43ba51426411
# ╟─57a0fe58-8796-11eb-0fb7-332dfa29bcf5
# ╟─5516fd22-8796-11eb-20ba-ef554955d071
# ╟─54fc4464-8796-11eb-324b-0195abc8c4bc
# ╟─54e17b3e-8796-11eb-3b20-45a682c2065b
# ╠═54c55cce-8796-11eb-3774-058c344a886a
# ╟─f04752ac-8795-11eb-03a3-a56cd548eca8
# ╠═1ae849ac-8797-11eb-2376-4189345388ff
# ╠═1accb14c-8797-11eb-1c71-81ad55c19b1d
# ╠═1ab6cef4-8797-11eb-2f81-6140be1f3b0d
# ╟─1a9b44e0-8797-11eb-3e92-2f8dbd5282ae
# ╠═1a7e108c-8797-11eb-11f8-f7be28d15ab9
# ╠═6c41590e-87c6-11eb-17e2-bd5698ae552f
# ╟─43e21928-8797-11eb-0e87-cd5790ab59b1
# ╠═43cc0f16-8797-11eb-1cca-874f98f40302
# ╠═d60b7a5a-87c5-11eb-0685-5100caf67827
# ╠═d5ea675c-87c5-11eb-1bf4-d1ea71c95392
# ╠═ab758052-8797-11eb-1193-b1f04d8f328c
# ╠═0388a272-87c7-11eb-21e9-9190a2f07a1c

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using Test
using UnROOT
using DataFrames
using StaticArrays
using MD5

const SAMPLES_DIR = joinpath(@__DIR__, "samples")

# @stack
struct A
    n::Int32
    m::Int64
end
struct B
    o::Float16
    p
    q::Bool
end
struct C
    r
    s
    t
    u
end
@UnROOT.stack D A B C
expected_fieldnames = (:n, :m, :o, :p, :q, :r, :s, :t, :u)
expected_fieldtypes = [Int32, Int64, Float16, Any, Bool, Any, Any, Any, Any]

@testset "utils" begin
    @test fieldnames(D) == expected_fieldnames
    @test [fieldtype(D, f) for f in fieldnames(D)] == expected_fieldtypes
end

# io.jl
UnROOT.@io struct Foo
    a::Int32
    b::Int64
    c::Float32
    d::SVector{5, UInt8}
end

@testset "@io" begin

    d = SA{UInt8}[1, 2, 3, 4, 5]

    foo = Foo(1, 2, 3, d)

    @test foo.a == 1
    @test foo.b == 2
    @test foo.c ≈ 3
    @test d == foo.d

    @test_skip 21 == sizeof(Foo)

    buf = IOBuffer(Vector{UInt8}(1:sizeof(Foo)))
    foo = UnROOT.unpack(buf, Foo)

    @test foo.a == 16909060
    @test foo.b == 361984551142689548
    @test foo.c ≈ 4.377526f-31
    @test foo.d == UInt8[0x11, 0x12, 0x13, 0x14, 0x15]
end

struct Bar
    x::Int8
    y::UInt16
end

@testset "io functions" begin
    @test 21 == UnROOT.packedsizeof(Foo)
    @test 3 == UnROOT.packedsizeof(Bar)
end

@testset "Header and Preamble" begin
    fobj = open(joinpath(SAMPLES_DIR, "km3net_online.root"))
    file_preamble = UnROOT.unpack(fobj, UnROOT.FilePreamble)
    @test "root" == String(file_preamble.identifier)

    file_header = UnROOT.unpack(fobj, UnROOT.FileHeader32)
    @test 100 == file_header.fBEGIN
end

@testset "ROOTFile" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_histos.root"))
    show(rootfile)
    @test 100 == rootfile.header.fBEGIN
    @test 1 == length(rootfile.directory.keys)
    @test "t1" ∈ keys(rootfile)
    for key in keys(rootfile)
        rootfile[key]
    end

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_custom_struct.root"))
    @test 100 == rootfile.header.fBEGIN
    @test 1 == length(rootfile.directory.keys)
    @test "T" ∈ keys(rootfile)
    for key in keys(rootfile)
        rootfile[key]
    end

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "histograms.root"))
    for branch in ["one", "two", "three"]
        @test branch in keys(rootfile)
    end

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    @test 100 == rootfile.header.fBEGIN
    @test 10 == length(rootfile.directory.keys)
    @test "E" ∈ keys(rootfile)
    @test "META" ∈ keys(rootfile)
    @test "JTRIGGER::JTriggerParameters" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_L0" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_L1" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_L2" ∈ keys(rootfile)
    @test "KM3NET_TIMESLICE_SN" ∈ keys(rootfile)
    @test "KM3NET_EVENT" ∈ keys(rootfile)
    @test "KM3NET_SUMMARYSLICE" ∈ keys(rootfile)
    # for key in keys(rootfile)
    #     @test !ismissing(rootfile[key])
    # end
end

@testset "Compressions" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array_lzma.root"))
    @test rootfile isa ROOTFile
    arr = array(rootfile, "t1/float_array")
    @test 100000 == length(arr)
    @test [0.0, 1.0588236, 2.1176472, 3.1764705, 4.2352943] ≈ arr[1:5] atol=1e-7
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array_lz4.root"))
    arr = array(rootfile, "t1/float_array")
    @test 100000 == length(arr)
    @test [0.0, 1.0588236, 2.1176472, 3.1764705, 4.2352943] ≈ arr[1:5] atol=1e-7
end

@testset "ROOTDirectoryHeader" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_histos.root"))
    header = rootfile.directory.header
    @test 5 == header.fVersion
    @test 1697049339 == header.fDatimeC
    @test 1697049339  == header.fDatimeM
    @test 111 == header.fNbytesKeys
    @test 78 == header.fNbytesName
    @test 100 == header.fSeekDir
    @test 0 == header.fSeekParent
    @test 1398 == header.fSeekKeys

    # rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    # header = rootfile.directory.header
    # @test 5 == header.fVersion
    # @test 1658540644 == header.fDatimeC
    # @test 1658540645 == header.fDatimeM
    # @test 629 == header.fNbytesKeys
    # @test 68 == header.fNbytesName
    # @test 100 == header.fSeekDir
    # @test 0 == header.fSeekParent
    # @test 1619244 == header.fSeekKeys
end

@testset "array()" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_histos.root"))
    arr = array(rootfile, "t1/mynum")
    @test 25 == length(arr)
    @test [0, 1, 2, 3, 4] ≈ arr[1:5] atol=0.1
    @test [10, 10, 10, 10, 10, 10] ≈ arr[20:end] atol=0.1

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array.root"))
    arr = array(rootfile, "t1/int32_array")
    @test 100000 == length(arr)
    @test [0, 1, 2, 3, 4] ≈ arr[1:5] atol=0.1
    arr = array(rootfile, "t1/float_array")
    @test 100000 == length(arr)
    @test [0.0, 1.0588236, 2.1176472, 3.1764705, 4.2352943] ≈ arr[1:5] atol=1e-7

    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    arr = array(rootfile, "KM3NET_EVENT/KM3NET_EVENT/snapshotHits"; raw=true)
end

@testset "DataFrame()" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "tree_with_large_array.root"))
    df = DataFrame(rootfile, "t1")
    @test 100000 == length(df.int32_array)
    @test [0, 1, 2, 3, 4] ≈ df.int32_array[1:5] atol=0.1
    @test [0.0, 1.0588236, 2.1176472, 3.1764705, 4.2352943] ≈ df.float_array[1:5] atol=1e-7
end

@testset "readbasketsraw()" begin
    array_md5 = [0xb4, 0xe9, 0x32, 0xe8, 0xfb, 0xff, 0xcf, 0xa0, 0xda, 0x75, 0xe0, 0x25, 0x34, 0x9b, 0xcd, 0xdf]
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "km3net_online.root"))
    data, offsets = array(rootfile, "KM3NET_EVENT/KM3NET_EVENT/snapshotHits"; raw=true)
    @test array_md5 == md5(data)
end

# Issues

@testset "issues" begin
    rootfile = ROOTFile(joinpath(SAMPLES_DIR, "issue7.root"))
    @test 2 == length(keys(rootfile))
    @test [1.0, 2.0, 3.0] == array(rootfile, "TreeD/nums")
    @test [1.0, 2.0, 3.0] == array(rootfile, "TreeF/nums")
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module ReadWriteSet

using CompilerTools
using CompilerTools.LambdaHandling
using CompilerTools.Helper

import ..DebugMsg
DebugMsg.init()

import Base.show

#LHSVar = Union{Symbol, GenSym}

"""
Holds which scalars and which array are accessed and for array which index expressions are used.
"""
type AccessSet
    scalars :: Set{LHSVar}
    arrays  :: Dict{LHSVar, Array{Array{Any,1},1}}
    AccessSet() = new(Set{LHSVar}(), Dict{LHSVar, Array{Array{Any,1},1}}())
end

"""
Stores which scalars and arrays are read or written in some code region.
"""
type ReadWriteSetType
    readSet  :: AccessSet
    writeSet :: AccessSet
    ReadWriteSetType() = new(AccessSet(),AccessSet())
end

function getArraysAccessed(rws :: ReadWriteSetType)
    return union(collect(keys(rws.readSet.arrays)), collect(keys(rws.writeSet.arrays)))
end

export from_exprs, ReadWriteSetType, AccessSet, set_debug_level, isRead, isWritten

function merge!(outer :: AccessSet, inner :: AccessSet)
    union!(outer.scalars, inner.scalars)
    for i in inner.arrays
        key   = i[1]
        value = i[2]
        if haskey(outer.arrays, key)
            append!(outer.arrays[key], value)
        else
            outer.arrays[key] = value
        end
    end
end

function merge!(outer :: ReadWriteSetType, inner :: ReadWriteSetType)
    merge!(outer.readSet,  inner.readSet)
    merge!(outer.writeSet, inner.writeSet)
end

"""
Return true if some symbol in "sym" is read either as a scalar or array within the computed ReadWriteSetType.
"""
function isRead(sym :: LHSVar, rws :: ReadWriteSetType)
    if in(sym, rws.readSet.scalars)
        return true
    elseif haskey(rws.readSet.arrays, sym)
        return true
    else
        return false
    end
end

"""
Return true if some symbol in "sym" is written either as a scalar or array within the computed ReadWriteSetType.
"""
function isWritten(sym :: LHSVar, rws :: ReadWriteSetType)
    if in(sym, rws.writeSet.scalars)
        return true
    elseif haskey(rws.writeSet.arrays, sym)
        return true
    else
        return false
    end
end

const CallbackType = Union{Function, Void}

"""
Walk through a lambda expression.
We just need to recurse through the lambda body and can ignore the rest.
"""
function from_lambda(ast :: Expr, depth, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY)
  # :lambda expression
  (linfo, body) = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast)
  from_expr(body, depth, rws, callback, cbdata, linfo)
  return ast
end

"""
Walk through an array of expressions.
Just recursively call from_expr for each expression in the array.
"""
function from_exprs(ast::Array, linfo :: LambdaVarInfo)
  from_exprs(ast, 1, ReadWriteSetType(), nothing, nothing, linfo)
end

"""
Walk through an array of expressions.
Just recursively call from_expr for each expression in the array.
Takes a callback and an opaque object so that non-standard Julia AST nodes can be processed via the callback.
"""
function from_exprs(ast :: Array, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
  from_exprs(ast, 1, ReadWriteSetType(), callback, cbdata, linfo)
end

"""
Walk through one AST node.
Calls the main internal walking function after initializing an empty ReadWriteSetType.
"""
function from_expr(ast :: ANY, linfo :: LambdaVarInfo)
  from_expr(ast, 1, ReadWriteSetType(), nothing, nothing, linfo)
end

"""
Walk through one AST node.
Calls the main internal walking function after initializing an empty ReadWriteSetType.
Takes a callback and an opaque object so that non-standard Julia AST nodes can be processed via the callback.
"""
function from_expr(ast :: ANY, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
  from_expr(ast, 1, ReadWriteSetType(), callback, cbdata, linfo)
end

"""
Walk through an array of expressions.
Just recursively call from_expr for each expression in the array.
Takes a callback and an opaque object so that non-standard Julia AST nodes can be processed via the callback.
Takes a ReadWriteSetType in "rws" into which information will be stored.
"""
function from_exprs(ast :: Array, depth :: Integer, rws :: ReadWriteSetType, callback :: Union{Function,Void}, cbdata :: ANY, linfo :: LambdaVarInfo)
  # sequence of expressions
  # ast = [ expr, ... ]
  local len  = length(ast)
  @dprintln(3, "RWS starting with ", ast)
  for i = 1:len
    @dprintln(2,"RWS Processing ast #",i," depth=", depth, " ", ast[i])
    from_expr(ast[i], depth, rws, callback, cbdata, linfo)
  end
  rws
end

"""
Walk through a tuple.
Just recursively call from_exprs on the internal tuple array to process each part of the tuple.
"""
function from_tuple(ast :: Array, depth :: Integer, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
  from_exprs(ast, depth+1, rws, callback, cbdata, linfo)
end

"""
Process a :(::) AST node.
Just process the symbol part of the :(::) node in ast[1] (which is args of the node passed in).
"""
function from_coloncolon(ast :: Array, depth :: Integer, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
  assert(length(ast) == 2)
  from_expr(ast[1], depth+1, rws, callback, cbdata, linfo)
end

"""
Process an assignment AST node.
The left-hand side of the assignment is added to the writeSet.
"""
function from_assignment(ast :: Array{Any,1}, depth :: Integer, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
  assert(length(ast) == 2)
  local lhs = ast[1]
  local rhs = ast[2]
  # lhs_type = typeof(lhs)
  push!(rws.writeSet.scalars, toLHSVar(lhs))
  from_expr(rhs, depth, rws, callback, cbdata, linfo)
end

"""
Takes a dictionary of symbol to an array of index expression.
Takes an array in "array_name" being accessed with expression "index_expr".
Makes sure there is an entry in the dictionary for this array and adds the index expression to this array.
"""
function addIndexExpr!(this_dict, array_name, index_expr)
  @dprintln(2,"addIndexExpr! ", typeof(array_name), " index_expr = ", index_expr, " typeof(index_expr) = ", typeof(index_expr))
  key = toLHSVar(array_name)
  if(!haskey(this_dict, key))
    this_dict[key] = Array{Any,1}[]
  end
  push!(this_dict[key], index_expr)
end

"""
Process :call Expr nodes to find arrayref and arrayset calls and adding the corresponding index expressions to the read and write sets respectively.
"""
function from_call(ast :: Expr, depth :: Integer, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
  local fun  = getCallFunction(ast)
  local args = getCallArguments(ast)
  @dprintln(2,fun)
  for i = 1:length(args)
    @dprintln(2,"RWS from_call first arg[",i,"] = ",args[i], " type = ", typeof(args[i]))
  end
  if(isBaseFunc(fun, :arrayref) || isBaseFunc(fun, :unsafe_arrayref))
    @dprintln(2,"Handling arrayref in from_call")
    # args[1]  = array
    # args[2+] = index expressions
    assert(length(args) >= 2)
    indices = args[2:end]
    @dprintln(3, "indices = ", indices, " typeof(indices) = ", typeof(indices))
    addIndexExpr!(rws.readSet.arrays, args[1], indices)
    for j = 1:length(indices)
      from_expr(indices[j], depth, rws, callback, cbdata, linfo)
    end
  elseif (isBaseFunc(fun, :arrayset) || isBaseFunc(fun, :unsafe_arrayset))
    @dprintln(2,"Handling arrayset in from_call, length(args) = ",length(args))
    # args[1]  = array
    # args[2]  = value
    # args[3+] = index expression
    assert(length(args) >= 3)
    indices = args[3:end]
    @dprintln(3, "indices = ", indices, " typeof(indices) = ", typeof(indices))
    addIndexExpr!(rws.writeSet.arrays, args[1], indices)
    from_expr(args[2], depth, rws, callback, cbdata, linfo)
    for j = 1:length(indices)
      from_expr(indices[j], depth, rws, callback, cbdata, linfo)
    end
  elseif isBaseFunc(fun, :setfield!)
    @dprintln(2,"Handling setfield! in from_call")
    obj = toLHSVar(args[1])
    if isa(obj,LHSVar)
      push!(rws.writeSet.scalars, obj)
    end
    from_exprs(args[2:end], depth+1, rws, callback, cbdata, linfo)
  else
    @dprintln(2,"Unhandled function ", fun, " in from_call")
    from_exprs(args, depth+1, rws, callback, cbdata, linfo)
  end
end

"""
If an AST node is not recognized then we try the passing the node to the callback to see if
it was able to process it.  If so, then we process the regular Julia statement returned by
the callback.
"""
function tryCallback(ast :: ANY, callback :: CallbackType, cbdata :: ANY, depth :: Integer, rws :: ReadWriteSetType, linfo :: LambdaVarInfo)
  if callback != nothing
    res = callback(ast, cbdata)
  else
    res = nothing
  end
  if res != nothing
    from_exprs(res, depth+1, rws, callback, cbdata, linfo)
    return false
  end
  return true
end

if VERSION >= v"0.6.0-pre"
import ..LambdaHandling.LambdaInfo
end

function from_expr(ast :: LambdaInfo, depth :: Integer, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY)
    (linfo, body) = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast)
    from_expr(body, depth, rws, callback, cbdata, linfo)
    return rws
end

function from_expr(ast :: Expr, depth :: Integer, rws :: ReadWriteSetType, callback :: CallbackType, cbdata :: ANY, linfo :: LambdaVarInfo)
    @dprintln(2,"RWS from_expr depth=", depth," ")

    @dprint(2,"RWS Expr ")
    local head = ast.head
    local args = ast.args
    local typ  = ast.typ
    @dprintln(2,head, " ", args)

    if !tryCallback(ast, callback, cbdata, depth, rws, linfo)
        return rws
    end

    if head == :lambda
        from_lambda(ast, depth, rws, callback, cbdata, linfo)
    elseif head == :body
        from_exprs(args, depth+1, rws, callback, cbdata, linfo)
    elseif head == :block || head == :(.)
        from_exprs(args, depth+1, rws, callback, cbdata, linfo)
    elseif head == :(=)
        from_assignment(args, depth, rws, callback, cbdata, linfo)
    elseif head == :return
        from_exprs(args, depth, rws, callback, cbdata, linfo)
    elseif head == :invoke || head == :call || head == :call1 || head == :foreigncall
        from_call(ast, depth, rws, callback, cbdata, linfo)
        # TODO?: tuple
    elseif head == Symbol("'")
        from_exprs(args, depth, rws, callback, cbdata, linfo)
    elseif head == :line
        # skip
    elseif head == :copy
        @dprintln(3,"RWS copyast type")
        # skip
    elseif head == :copyast
        @dprintln(3,"RWS copyast type")
        # skip
    elseif head == :arraysize
        @dprintln(3,"RWS arraysize")
        # skip
    elseif head == :assertEqShape
        @dprintln(3,"RWS assertEqShape")
        # skip
    elseif head == :alloc
        from_expr(args[2], depth, rws, callback, cbdata, linfo)
    elseif head == :tuple
        @dprintln(2,"RWS tuple")
        from_tuple(args, depth, rws, callback, cbdata, linfo)
        # skip
    elseif head == :(::)
        @dprintln(2,"RWS ::")
        from_coloncolon(args, depth, rws, callback, cbdata, linfo)
    elseif head == :getindex || head == :new || head == :type_goto || head == :ccall
        from_exprs(args, depth+1, rws, callback, cbdata, linfo)
    elseif head == :gotoifnot
        from_expr(args[1], depth, rws, callback, cbdata, linfo)
    elseif head == :meta || head == :inbounds || head == :static_parameter ||
           head == :enter || head == :leave || head == :curly || head == :the_exception ||
           head == :& || head == :static_typeof || head == :(->) || head == :(&&) ||
           head == :(||) || head == :break || head == :continue || head == Symbol("'") ||
           head == :(...) || head == :parameters || head == :kw || head == :macrocall ||
           head == :simdloop || head == :quote || head == :local || head == :let ||
           head == :while || head == :comparison || head == :if ||
           head in Set([:(+=), :(/=), :(*=), :(-=)]) || head == :for || head == :const ||
           head == :llvmcall
        # skip
    else
        #println("from_expr: unknown Expr head :", head)
#        if tryCallback(ast, callback, cbdata, depth, rws)
            throw(string("from_expr: unknown Expr head :", head))
#        end
    end

    return rws
end

function from_expr(ast::Union{LabelNode,GotoNode,LineNumberNode,DataType,Module,NewvarNode},
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    return rws
end

function from_expr(ast::RHSVar,
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    vtyp = CompilerTools.LambdaHandling.getType(ast, linfo)
    @dprintln(3,"RWS ", ast, " type = ", typeof(ast), " vtyp = ", vtyp)
    if isArrayType(vtyp)
        # If a RHSVar is an array type and accessed in this way then it could potentially be read in its entirety.
        # We put zeros in the indices to indicate this and this should trigger all code to realize that this array
        # is accessed outside any loop indices.
        @dprintln(3,"Adding array to readSet with all zero indices.")
        addIndexExpr!(rws.readSet.arrays, toLHSVar(ast), zeros(Int, ndims(vtyp)))
    else
        push!(rws.readSet.scalars, toLHSVar(ast))
    end

    return rws
end

function from_expr(ast::Union{TopNode,AbstractString},
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    # skip
    @dprintln(3,"RWS ", typeof(ast), " type")

    return rws
end

function from_expr(ast::GlobalRef,
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    local mod  = ast.mod
    local name = ast.name
    @dprintln(3,"RWS GlobalRef type ",typeof(mod))
    # warn(string("from_expr: GetfieldNode typeof(mod)=", typeof(mod)))

    return rws
end

function from_expr(ast::QuoteNode,
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    local value = ast.value
    # TODO: fields: value
    @dprintln(3,"RWS QuoteNode type ",typeof(value))
    # warn(string("from_expr: QuoteNode typeof(value)=", typeof(value)))

    return rws
end

function from_expr(ast::SimpleVector,
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    @dprintln(3,"RWS SimpleVector ast = ", ast)
    for i = 1:length(ast)
        @dprintln(2,"RWS Processing ast #",i," depth=", depth, " ", ast[i])
        from_expr(ast[i], depth, rws, callback, cbdata, linfo)
    end

    return rws
end

function from_expr(ast::Any,
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    asttyp = typeof(ast)

    # For Julia 0.5, Symbol won't be in RHSVar and so will come here.
    # All such uses can currently be ignored.
    if isbits(asttyp) || asttyp == Symbol
        #skip
    else
        if tryCallback(ast, callback, cbdata, depth, rws, linfo)
            throw(string("from_expr: unknown ast :", asttyp))
        end
        #@dprintln(2,"RWS from_expr: unknown AST (", typeof(ast), ",", ast, ")")
        #warn(string("from_expr: unknown AST (", typeof(ast), ",", ast, ")"))
    end
    return rws
end

function from_expr(ast::ReadWriteSetType,
                   depth :: Integer,
                   rws :: ReadWriteSetType,
                   callback :: CallbackType,
                   cbdata::ANY,
                   linfo :: LambdaVarInfo)
    merge!(rws, ast)
    return rws
end

function show(io::IO, pnode::ReadWriteSetType)
    println(io,"read set: ")
    print(io, pnode.readSet)
    println(io,"write set: ")
    print(io, pnode.writeSet)
end

function show(io::IO, pnode::AccessSet)
    println(io, "  scalars: ", pnode.scalars)
    print(io, "  arrays: ")
    for (k,v) in pnode.arrays
        print(io, k," => ", v,", ")
    end
    println(io)
end

end


text:
using PGA2D
using Test

@testset "PGA2D.jl" begin

    @testset "Coordinates" begin
        @test point_coordinates(point(1.0,2.0)) == (1.0,2.0)
        @test direction_coordinates(direction(1.0,2.0)) == (1.0,2.0)
        @test line_coordinates(line(1.0,2.0,3.0)) == (1.0,2.0,3.0)

        @test_throws DomainError point_coordinates(direction(1.0,2.0))
        @test_throws DomainError point_coordinates(line(1.0,2.0,3.0))
        @test_throws DomainError direction_coordinates(point(1.0,2.0))
        @test_throws DomainError direction_coordinates(line(1.0,2.0,3.0))
        @test_throws DomainError line_coordinates(point(1.0,2.0))
        @test_throws DomainError line_coordinates(direction(1.0,2.0))

        @test is_point(point(1,2))
        @test is_line(line(1,2,3))
        @test is_direction(direction(1,2,))

        @test !is_point(direction(1,2))
        @test !is_point(line(1,2,3))
        @test !is_direction(point(1,2))
        @test !is_direction(line(1,2,3))
        @test !is_line(point(1,2))
        @test !is_line(direction(1,2))
    end

   @testset "Geometric operations" begin
        P1 = point(0,0)
        P2 = point(2,0)
        P3 = point(3,1)
        P4 = point(1,2)

        l1 = join_pp(P1,P2)
        l2 = join_pp(P2,P3)
        l3 = join_pp(P3,P4)
        l4 = join_pp(P4,P1)

        @test point_coordinates(meet_ll(l1,l2)) == point_coordinates(P2)
        @test point_coordinates(meet_ll(l2,l3)) == point_coordinates(P3)
        @test point_coordinates(meet_ll(l3,l4)) == point_coordinates(P4)
        @test point_coordinates(meet_ll(l4,l1)) == point_coordinates(P1)

        @test point_coordinates(meet_ll(join_pp(P1,P2),join_pp(P3,P4))) == (5.0,0.0)

   end

end


text:
using CSV
using DataFrames
using RCall

@rimport readr as readr

# load functions to fit data via least squares
println(snakemake)
include(string(snakemake.scriptdir, "/optim.jl"))

println("Read in data...")
println(snakemake.input["oesophagusmetadata"])

DFdonor = rcopy(readr.read_csv(snakemake.input["oesophagusmetadata"]))

println("Reading in per gene data...")
DF = rcopy(readr.read_csv(snakemake.input["all"]))
DF[:A] = 2 .* DF[:cutoff];
DFnon = rcopy(readr.read_csv(snakemake.input["nonsense"], guess_max = 10^5))
DFnon[:A] = 2 .* DFnon[:cutoff];
DFmis = rcopy(readr.read_csv(snakemake.input["missense"], guess_max = 10^5))
DFmis[:A] = 2 .* DFmis[:cutoff];
println("Data read in")

myDFall = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :gene, :nmutations], 0)
myDFmiss = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :gene, :nmutations], 0)
myDFnon = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, String, String, Float64, String, Float64],
    [:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda,
    :rsq, :patient, :Age, :Age2, :gene, :nmutations], 0)

println("Analyse per patient per gene data for all mutations...")
for gene in unique(DF[:gene])
    println("Analysing gene: $gene")
    DFgene = filter(row -> row[:gene] == gene, DF);
    println("Analysing all mutations")
    for p in unique(DF[:ID])
        DFpatient = filter(row -> row[:ID] == p, DFgene);
        age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]

        nmuts = DFpatient[:muts][end]
        if nmuts < 5
            println("Patient $p, only $(nmuts) mutations, moving to next patient")
            continue
        end
        println("Analysing patient $p, $(nmuts) mutations")

        x = LLoptimizationresults(DFpatient[:dnds], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        x.DF[:patient] = p
        x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
        x.DF[:Age2] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = nmuts
        append!(myDFall, x.DF)
    end
end

println("Number of fits: $(length(unique(myDFall[:gene])))")

println("Analyse per patient per gene data for nonsense mutations...")
for gene in unique(DFnon[:gene])
    println("Analysing gene: $gene")
    DFgene = filter(row -> row[:gene] == gene, DFnon);
    println("Analysing nonsense mutations")
    println(DFgene)
    for p in unique(DFnon[:ID])
        DFpatient = filter(row -> row[:ID] == p, DFgene);
        age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]

        if length(DFpatient[:muts]) == 0
            println("Missing fit: $p, $gene")
            continue
        end

        nmuts = DFpatient[:muts][end]
        println(nmuts)
        if nmuts < 5
            println("Patient $p, only $(nmuts) mutations, moving to next patient")
            continue
        end
        println("Analysing patient $p, $(nmuts) mutations")

        x = LLoptimizationresults(DFpatient[:dnds], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        x.DF[:patient] = p
        x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
        x.DF[:Age2] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = nmuts
        append!(myDFnon, x.DF)
    end
end

println("Number of fits: $(length(unique(myDFnon[:gene])))")

println("Analyse per patient per gene data for missense mutations...")
for gene in unique(DF[:gene])
    println("Analysing gene: $gene")
    DFgene = filter(row -> row[:gene] == gene, DFmis);
    println("Analysing missense mutations")
    for p in unique(DFmis[:ID])
        DFpatient = filter(row -> row[:ID] == p, DFgene);
        age = filter(row -> row[:patient] == p, DFdonor)[:Age2][1]

        nmuts = DFpatient[:muts][end]
        if nmuts < 5
            println("Patient $p, only $(nmuts) mutations, moving to next patient")
            continue
        end
        println("Analysing patient $p, $(nmuts) mutations")

        x = LLoptimizationresults(DFpatient[:dnds], DFpatient[:A]; t = age, Amin = 2 * 0.01, ρ = 5000.0)
        x.DF[:patient] = p
        x.DF[:Age] = filter(row -> row[:patient] == p, DFdonor)[:Age][1]
        x.DF[:Age2] = age
        x.DF[:gene] = gene
        x.DF[:nmutations] = nmuts
        append!(myDFmiss, x.DF)
    end
end

println("Number of fits: $(length(unique(myDFmiss[:gene])))")

myDFall[:mutation_class] = "All"
myDFmiss[:mutation_class] = "Missense"
myDFnon[:mutation_class] = "Nonsense"

append!(myDFall, myDFmiss)
append!(myDFall, myDFnon)

@rput myDFall
R"""
library(readr)
write_csv(myDFall, $(snakemake.output["oesophagusfit"]))
""";


text:
#linearalgebra.jl
"""
a' * B * a
"""
function mulαtβα(a::AbstractVector, B::AbstractMatrix{T}) where T
    if length(a) != size(B, 2)::Int  || size(B, 2)::Int  != size(B, 1)::Int  error("Dimention error") end
    c = zero(T)
    for i ∈ axes(B, 1)::Base.OneTo{Int}
        ct = zero(T)
        for i2 ∈ axes(B, 2)::Base.OneTo{Int}
            @inbounds  ct += B[i, i2]*a[i2]
        end
        @inbounds c += ct*a[i]
    end
    c
end
"""
θ + A * B * A'

Change θ (only upper triangle).
"""
function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix) where T
    #1 alloc
    axb  = axes(B, 1)
    axa  = axes(A, 1)
    sa   = size(A, 1)
    c  = Vector{T}(undef, size(B, 1))
    for i ∈ axa
        fill!(c, zero(T))
        @inbounds  for n ∈ axb, m ∈ axb
            c[n] += B[m, n] * A[i, m]
        end
        #upper triangular n = i:p
        @inbounds  for n ∈ i:sa, m ∈ axb
            θ[i, n] += A[n, m] * c[m]
        end
    end
    θ
end
"""
    mulθ₃(y, X, β, V::AbstractMatrix{T})::T where T

(y - X * β)' * (-V) * (y - X * β)

use only upper triangle of V
"""
function mulθ₃(y, X, β, V::AbstractArray{T}) where T
    q = size(V, 1)
    p = size(X, 2)
    θ = zero(T)

    if q == 1
        c = zero(T)
        #=@turbo=# @inbounds  for m in 1:p
            c += X[1, m] * β[m]
        end
        return -V[1, 1] * (y[1] - c)^2
    end

    c = zeros(T, q)
    #=@turbo=# @inbounds  for n = 1:q, m = 1:p
        c[n] += X[n, m] * β[m]
    end

    @simd for n = 1:q
        @simd for m = n+1:q
            @inbounds θ -= V[n, m] * (y[n] - c[n]) * (y[m] - c[m]) * 2
        end
    end
    #=@turbo=# @inbounds  for m = 1:q
        θ -= V[m, m] * (y[m] - c[m]) ^ 2
    end

    return θ
end
"""
θ + A' * b

Change θ.
"""
function mulαtβinc!(θ::AbstractVector, A::AbstractMatrix, b::AbstractVector)
    q = size(A, 1)
    if q != length(b) throw(DimensionMismatch("size(A, 1) should be equal length(b)")) end
    p = size(A, 2)
    #=@turbo=# @inbounds  for n in 1:p, m in 1:q
        θ[n] += b[m] * A[m, n]
    end
    θ
end
################################################################################

#=
If length θ >= size(rz, 2)
!!! without checking
vec = rz * θ
=#
@inline function tmul_unsafe(rz, θ::AbstractVector{T}) where T
    vec = zeros(T, size(rz, 1))
    #=@turbo=# @inbounds  for r ∈ axes(rz, 1)
        for i ∈ axes(rz, 2)
            vec[r] += rz[r, i] * θ[i]
        end
    end
    vec
end


text:
#=

  Euler #31 in Julia.

  Problem 31
  """
  In England the currency is made up of pound, £, and pence, p, and
  there are eight coins in general circulation:

     1p, 2p, 5p, 10p, 20p, 50p, £1 (100p) and £2 (200p).

  It is possible to make £2 in the following way:

     1×£1 + 1×50p + 2×20p + 1×5p + 1×2p + 3×1p

  How many different ways can £2 be made using any number of coins?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# DP approach
function coins(c,money,m)
    len = length(c)
    if m == len
        return 1
    end

    sum1 = 0
    for i in m:len
        if money - c[i] == 0
            sum1 += 1
        end
        if  money - c[i] > 0
            sum1 += coins(c, money-c[i], i)
        end
    end
    return sum1;
end

# 0.00058022s
function euler31a()
    c = [200,100,50,20,10,5,2,1]
    return coins(c, 200, 1)
end

run_euler(euler31a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank