text:
uk_0_eval(x)= @evalpoly x 1.0

uk_1_eval(x)= @evalpoly x 0.0  0.125  0.0  -0.20833333333333334

uk_2_eval(x)= @evalpoly x 0.0  0.0  0.0703125  0.0  -0.4010416666666667  0.0  0.3342013888888889

uk_3_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0732421875  0.0  -0.8912109375  0.0  1.8464626736111112  0.0  -1.0258125964506173

uk_4_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0  0.112152099609375  0.0  -2.3640869140625  0.0  8.78912353515625  0.0  -11.207002616222995  0.0  4.669584423426247

uk_5_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0  0.0  0.22710800170898438  0.0  -7.368794359479631  0.0  42.53499874538846  0.0  -91.81824154324003  0.0  84.63621767460074  0.0  -28.212072558200244

uk_6_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0  0.0  0.0  0.5725014209747314  0.0  -26.491430486951554  0.0  218.1905117442116  0.0  -699.5796273761327  0.0  1059.9904525279999  0.0  -765.2524681411816  0.0  212.5701300392171

uk_7_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.7277275025844574  0.0  -108.09091978839464  0.0  1200.9029132163525  0.0  -5305.646978613405  0.0  11655.393336864536  0.0  -13586.550006434136  0.0  8061.722181737308  0.0  -1919.4576623184068

uk_8_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  6.074042001273483  0.0  -493.915304773088  0.0  7109.514302489364  0.0  -41192.65496889756  0.0  122200.46498301747  0.0  -203400.17728041555  0.0  192547.0012325315  0.0  -96980.5983886375  0.0  20204.29133096615

uk_9_eval(x)= @evalpoly x 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  24.380529699556064  0.0  -2499.830481811209  0.0  45218.76898136274  0.0  -331645.1724845636  0.0  1.2683652733216248e6  0.0  -2.813563226586534e6  0.0  3.763271297656404e6  0.0  -2.998015918538106e6  0.0  1.311763614662977e6  0.0  -242919.18790055133

const uk_10_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 110.01714026924674, 0.0, -13886.089753717039, 0.0, 308186.40461266245, 0.0, -2.785618128086455e6, 0.0, 1.328876716642182e7, 0.0, -3.756717666076335e7, 0.0, 6.634451227472903e7, 0.0, -7.410514821153264e7, 0.0, 5.095260249266463e7, 0.0, -1.970681911843222e7, 0.0, 3.2844698530720375e6, ]
uk_10_eval(x) = evalpoly(x, uk_10_coef)

const uk_11_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 551.3358961220206, 0.0, -84005.43360302408, 0.0, 2.24376817792245e6, 0.0, -2.4474062725738734e7, 0.0, 1.420629077975331e8, 0.0, -4.958897842750303e8, 0.0, 1.1068428168230145e9, 0.0, -1.621080552108337e9, 0.0, 1.5535968995705795e9, 0.0, -9.39462359681578e8, 0.0, 3.255730741857656e8, 0.0, -4.932925366450995e7, ]
uk_11_eval(x) = evalpoly(x, uk_11_coef)

const uk_12_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3038.0905109223845, 0.0, -549842.3275722886, 0.0, 1.739510755397817e7, 0.0, -2.2510566188941535e8, 0.0, 1.5592798648792577e9, 0.0, -6.563293792619284e9, 0.0, 1.79542137311556e10, 0.0, -3.302659974980072e10, 0.0, 4.128018557975397e10, 0.0, -3.463204338815877e10, 0.0, 1.868820750929582e10, 0.0, -5.866481492051846e9, 0.0, 8.14789096118312e8, ]
uk_12_eval(x) = evalpoly(x, uk_12_coef)

const uk_13_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 18257.75547429318, 0.0, -3.871833442572612e6, 0.0, 1.4315787671888906e8, 0.0, -2.167164983223796e9, 0.0, 1.763473060683497e10, 0.0, -8.786707217802325e10, 0.0, 2.879006499061506e11, 0.0, -6.453648692453765e11, 0.0, 1.008158106865382e12, 0.0, -1.098375156081223e12, 0.0, 8.19218669548577e11, 0.0, -3.990961752244664e11, 0.0, 1.1449823773202577e11, 0.0, -1.4679261247695614e10, ]
uk_13_eval(x) = evalpoly(x, uk_13_coef)

const uk_14_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 118838.42625678328, 0.0, -2.918838812222081e7, 0.0, 1.247009293512711e9, 0.0, -2.1822927757529232e10, 0.0, 2.0591450323241003e11, 0.0, -1.1965528801961816e12, 0.0, 4.612725780849132e12, 0.0, -1.2320491305598287e13, 0.0, 2.334836404458184e13, 0.0, -3.1667088584785152e13, 0.0, 3.056512551993531e13, 0.0, -2.051689941093443e13, 0.0, 9.109341185239896e12, 0.0, -2.406297900028503e12, 0.0, 2.8646403571767896e11, ]
uk_14_eval(x) = evalpoly(x, uk_14_coef)

const uk_15_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 832859.3040162894, 0.0, -2.3455796352225152e8, 0.0, 1.1465754899448242e10, 0.0, -2.2961937296824658e11, 0.0, 2.4850009280340854e12, 0.0, -1.663482472489248e13, 0.0, 7.437312290867914e13, 0.0, -2.3260483118893994e14, 0.0, 5.230548825784446e14, 0.0, -8.574610329828949e14, 0.0, 1.0269551960827622e15, 0.0, -8.894969398810261e14, 0.0, 5.427396649876595e14, 0.0, -2.2134963870252512e14, 0.0, 5.417751075510603e13, 0.0, -6.019723417234003e12, ]
uk_15_eval(x) = evalpoly(x, uk_15_coef)

const uk_16_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.252951493434797e6, 0.0, -2.0016469281917765e9, 0.0, 1.1099740513917906e11, 0.0, -2.5215584749128555e12, 0.0, 3.1007436472896465e13, 0.0, -2.3665253045164925e14, 0.0, 1.2126758042503475e15, 0.0, -4.3793258383640155e15, 0.0, 1.1486706978449752e16, 0.0, -2.226822513391114e16, 0.0, 3.213827526858623e16, 0.0, -3.4447226006485136e16, 0.0, 2.705471130619707e16, 0.0, -1.5129826322457674e16, 0.0, 5.705782159023669e15, 0.0, -1.301012723549699e15, 0.0, 1.3552215870309362e14, ]
uk_16_eval(x) = evalpoly(x, uk_16_coef)

const uk_17_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0069589531988926e7, 0.0, -1.807822038465807e10, 0.0, 1.1287091454108745e12, 0.0, -2.886383763141477e13, 0.0, 4.000444570430363e14, 0.0, -3.450385511846272e15, 0.0, 2.0064271476309532e16, 0.0, -8.270945651585064e16, 0.0, 2.4960365126160426e17, 0.0, -5.62631788074636e17, 0.0, 9.575335098169137e17, 0.0, -1.233611693196069e18, 0.0, 1.1961991142756303e18, 0.0, -8.592577980317544e17, 0.0, 4.4347954614171885e17, 0.0, -1.5552983504313898e17, 0.0, 3.3192764720355212e16, 0.0, -3.2541926196426675e15, ]
uk_17_eval(x) = evalpoly(x, uk_17_coef)

const uk_18_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.2593921650476694e8, 0.0, -1.7228323871735056e11, 0.0, 1.2030115826419195e13, 0.0, -3.4396530474307606e14, 0.0, 5.33510697870884e15, 0.0, -5.160509319348522e16, 0.0, 3.376676249790609e17, 0.0, -1.5736434765189596e18, 0.0, 5.402894876715981e18, 0.0, -1.3970803516443374e19, 0.0, 2.7572829816505184e19, 0.0, -4.1788614446568374e19, 0.0, 4.859942729324835e19, 0.0, -4.301555703831442e19, 0.0, 2.8465212251676553e19, 0.0, -1.3639420410571586e19, 0.0, 4.4702009640123085e18, 0.0, -8.966114215270461e17, 0.0, 8.301957606731907e16, ]
uk_18_eval(x) = evalpoly(x, uk_18_coef)

const uk_19_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.836255180230434e9, 0.0, -1.7277040123530002e12, 0.0, 1.3412416915180642e14, 0.0, -4.2619355104269e15, 0.0, 7.351663610930973e16, 0.0, -7.921651119323832e17, 0.0, 5.789887667664653e18, 0.0, -3.0255665989903716e19, 0.0, 1.1707490535797255e20, 0.0, -3.434621399768417e20, 0.0, 7.756704953461136e20, 0.0, -1.3602037772849937e21, 0.0, 1.8571089321463448e21, 0.0, -1.9677247077053117e21, 0.0, 1.601689857369359e21, 0.0, -9.824438427689853e20, 0.0, 4.39279220088871e20, 0.0, -1.3512175034359957e20, 0.0, 2.556380296052923e19, 0.0, -2.242438856186774e18, ]
uk_19_eval(x) = evalpoly(x, uk_19_coef)

const uk_20_coef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.646840080706557e10, 0.0, -1.8187262038511043e13, 0.0, 1.5613123930484675e15, 0.0, -5.484033603883292e16, 0.0, 1.0461721131134348e18, 0.0, -1.2483700995047234e19, 0.0, 1.0126774169536592e20, 0.0, -5.891794135069496e20, 0.0, 2.548961114664971e21, 0.0, -8.40591581710835e21, 0.0, 2.1487414815055883e22, 0.0, -4.302534303482378e22, 0.0, 6.7836616429518815e22, 0.0, -8.423222750084318e22, 0.0, 8.194331005435126e22, 0.0, -6.173206302884411e22, 0.0, 3.5284358439034075e22, 0.0, -1.478774352843361e22, 0.0, 4.285296082829493e21, 0.0, -7.671943936729004e20, 0.0, 6.393286613940834e19, ]
uk_20_eval(x) = evalpoly(x, uk_20_coef)

const UK_POLYS = [uk_0_eval, uk_1_eval, uk_2_eval, uk_3_eval, uk_4_eval,
    uk_5_eval, uk_6_eval, uk_7_eval, uk_8_eval, uk_9_eval, uk_10_eval,
    uk_11_eval, uk_12_eval, uk_13_eval, uk_14_eval, uk_15_eval, uk_16_eval,
    uk_17_eval, uk_18_eval, uk_19_eval, uk_20_eval]


text:
#! /bin/bash julia --project generator.jl
using Pkg
using Pkg.Artifacts
using Clang.Generators
using Clang.Generators.JLLEnvs
using SuiteSparse_jll

cd(@__DIR__)

# headers
SuiteSparse_toml = joinpath(dirname(pathof(SuiteSparse_jll)), "..", "StdlibArtifacts.toml")
SuiteSparse_dir = Pkg.Artifacts.ensure_artifact_installed("SuiteSparse", SuiteSparse_toml)

include_dir = joinpath(SuiteSparse_dir, "include") |> normpath

klu_h = joinpath(include_dir, "klu.h")
@assert isfile(klu_h)

# load common option
options = load_options(joinpath(@__DIR__, "generator.toml"))

# run generator for all platforms
for target in JLLEnvs.JLL_ENV_TRIPLES
    @info "processing $target"

    options["general"]["output_file_path"] = joinpath(@__DIR__, "..", "lib", "$target.jl")

    args = get_default_args(target)
    push!(args, "-I$include_dir")
    if startswith(target, "x86_64") || startswith(target, "powerpc64le") || startswith(target, "aarch64")
        push!(args, "-DSUN64 -DLONGBLAS='long long'")
    end

    header_files = [klu_h]

    ctx = create_context(header_files, args, options)

    build!(ctx)
end

ext: jl


text:
export KalmanSmoother, kalman_smoother

struct KalmanSmoother
    x
    P
    Po
    J
    S11
    S00
    S10
    Syy
    z
end

function kalman_smoother(m::LinearStateSpaceModel,z,u=fill([0.0],length(z)))
    kalman_smoother(kalman_filter(m,z,u),m)
end

function kalman_smoother(m::LinearStateSpaceModel,θ,kf::KalmanFilter)
    T = length(kf.xp)

    xf = kf.xf
    xp = kf.xp
    Pf = kf.Pf
    Pp = kf.Pp
    K  = kf.K

    F = state_transition_matrix(m,θ)
    H = observation_matrix(m,θ)

    x = fill(xf[end],T+1)
    P = fill(Pf[end],T+1)

    # Po[T] = P^T_{T,T-1} = (I-K[T])*F*P^{T-1}_{T-1}
    Po = fill((I-K[T]*H)*F*Pf[T],T)

    J = fill(Pf[end],T)
    S = fill(x[end]*x[end]' + P[end],T+1)

    S10 = zeros(size(Pf[end]))
    Syy = zeros(size(H*Pf[end]*H'))
    for t in T:-1:1

        if t == T
            J[t],x[t],P[t],_,S[t],_,syy = smoother_step(x[t+1],P[t+1],xf[t],Pf[t],xp[t],Pp[t],F,Pf[t],J[t],Po[t],H,kf.z[t])
            S10 += x[t+1]*x[t]' + Po[t]
            Syy += syy
        else
            J[t],x[t],P[t],Po[t],S[t],s10,syy = smoother_step(x[t+1],P[t+1],xf[t],Pf[t],xp[t],Pp[t],F,Pf[t+1],J[t+1],Po[t+1],H,kf.z[t])
            S10 += s10
            Syy += syy
        end

        #=
        # J[t] = J_{t-1}
        # x[t] = x^n_{t-1}
        # P[t] = P^n_{t-1}
        # xf[t] = x^{t-1}_{t-1}
        # Pf[t] = P^{t-1}_{t-1}
        # xp[t] = x^{t-1}_t
        # Pp[t] = P^{t-1}_t

        # J_{t-1} = P^{t-1}_{t-1}*F'*(P^{t-1}_t)^(-1)
        J[t] = Pf[t]*F'/Pp[t]

        # x^n_{t-1} = x^{t-1}_{t-1} + J_{t-1}*(x^n_t - x^{t-1}_t)
        x[t] = xf[t] + J[t]*(x[t+1]-xp[t])

        # P^n_{t-1} = P^{t-1}_{t-1} + J_{t-1}*(P^n_t - P^{t-1}_t)*J_{t-1}
        P[t] = Pf[t] + J[t]*(P[t+1]-Pp[t])*J[t]'

        # Po[t] = P^n_{t,t-1}
        # = P^{t}_{t} + J_{t}*(P^n_{t+1,t}-F*P^{t}_{t})*J_{t-1}
        if t == T
        else
                Po[t] = Pf[t+1]*J[t]' + J[t+1]*(Po[t+1]-F*Pf[t+1])*J[t]'
        end
        S[t] = x[t]*x[t]' + P[t]
        S10 += x[t+1]*x[t]' + Po[t]
        Syy += (kf.z[t].-H*x[t+1])*(kf.z[t].-H*x[t+1])' .+ H*P[t+1]*H'
        =#
    end

    S11 = sum(S[2:end])
    S00 = sum(S[1:end-1])

    KalmanSmoother(x,P,Po,J,S11,S00,S10,Syy,kf.z)
end

function kalman_smoother(kf::KalmanFilter,m::TimeVariantLinearStateSpaceModel)
    T = length(kf.xp)

    xf = kf.xf
    xp = kf.xp
    Pf = kf.Pf
    Pp = kf.Pp
    K  = kf.K

    F = m.F
    H = m.H

    x = fill(xf[end],T+1)
    P = fill(Pf[end],T+1)

    # Po[T] = P^T_{T,T-1} = (I-K[T])*F*P^{T-1}_{T-1}
    Po = fill((I-K[T]*H[end])*F[end]*Pf[T],T)

    J = fill(Pf[end],T)
    S = fill(x[end]*x[end]' + P[end],T+1)

    S10 = zeros(size(Pf[end]))
    Syy = zeros(size(H[end]*Pf[end]*H[end]'))
    for t in T:-1:1

        if t == T
            J[t],x[t],P[t],_,S[t],_,syy = smoother_step(x[t+1],P[t+1],xf[t],Pf[t],xp[t],Pp[t],F[t],Pf[t],J[t],Po[t],H[t],kf.z[t])
            S10 += x[t+1]*x[t]' + Po[t]
            Syy += syy
        else
            J[t],x[t],P[t],Po[t],S[t],s10,syy = smoother_step(x[t+1],P[t+1],xf[t],Pf[t],xp[t],Pp[t],F[t],Pf[t+1],J[t+1],Po[t+1],H[t],kf.z[t])
            S10 += s10
            Syy += syy
        end

    end

    S11 = sum(S[2:end])
    S00 = sum(S[1:end-1])

    KalmanSmoother(x,P,Po,J,S11,S00,S10,Syy,kf.z)
end

function smoother_step(x0,P0,xf,Pf,xp,Pp,F,Pf2,J2,Po2,H,z)
    # J[t] = J_{t-1}
    # x[t] = x^n_{t-1}
    # P[t] = P^n_{t-1}
    # xf[t] = x^{t-1}_{t-1}
    # Pf[t] = P^{t-1}_{t-1}
    # xp[t] = x^{t-1}_t
    # Pp[t] = P^{t-1}_t

    # J_{t-1} = P^{t-1}_{t-1}*F'*(P^{t-1}_t)^(-1)
    J = Pf*F'/Pp

    # x^n_{t-1} = x^{t-1}_{t-1} + J_{t-1}*(x^n_t - x^{t-1}_t)
    x = xf + J*(x0-xp)

    # P^n_{t-1} = P^{t-1}_{t-1} + J_{t-1}*(P^n_t - P^{t-1}_t)*J_{t-1}
    P = Pf + J*(P0-Pp)*J'

    # Po[t] = P^n_{t,t-1}
    # = P^{t}_{t} + J_{t}*(P^n_{t+1,t}-F*P^{t}_{t})*J_{t-1}
    Po = Pf2*J' + J2*(Po2-F*Pf2)*J'
    S = x*x' + P
    S10 = x0*x' + Po
    Syy = (z.-H*x0)*(z.-H*x0)' .+ H*P0*H'
    J,x,P,Po,S,S10,Syy
end


text:
"""
    apply(f,a::AbstractArray...) -> AbstractArray

Applies the kernel `f` to the entries of the arrays in `a` (see the definition of [`Kernel`](@ref)).

The resulting array `r` is such that `r[i]` equals to `apply_kernel(f,ai...)` where `ai`
is the tuple containing the `i`-th entry of the arrays in `a` (see function
[`apply_kernel`](@ref) for more details).
In other words, the resulting array is numerically equivalent to:

    map( (x...)->apply_kernel(f,x...), a...)

See the [`apply_kernel`](@ref) function for details.

# Examples

Using a function as kernel

```jldoctest
using Gridap.Arrays

a = collect(0:5)
b = collect(10:15)

c = apply(+,a,b)

println(c)

# output
[10, 12, 14, 16, 18, 20]
```

Using a user-defined kernel

```jldoctest
using Gridap.Arrays
import Gridap.Arrays: apply_kernel!

a = collect(0:5)
b = collect(10:15)

struct MySum <: Kernel end

apply_kernel!(cache,::MySum,x,y) = x + y

k = MySum()

c = apply(k,a,b)

println(c)

# output
[10, 12, 14, 16, 18, 20]
```

"""
function apply(f,a::AbstractArray...)
  s = common_size(a...)
  apply(Fill(f,s...),a...)
end

"""
    apply(::Type{T},f,a::AbstractArray...) where T

Like [`apply(f,a::AbstractArray...)`](@ref), but the user provides the element type
of the resulting array in order to circumvent type inference.
"""
function apply(::Type{T},f,a::AbstractArray...) where T
  s = common_size(a...)
  apply(T,Fill(f,s...),a...)
end

"""
    apply(f::AbstractArray,a::AbstractArray...) -> AbstractArray
Applies the kernels in the array of kernels `f` to the entries in the arrays in `a`.

The resulting array has the same entries as the one obtained with:

    map( apply_kernel, f, a...)

See the [`apply_kernel`](@ref) function for details.

# Example

"Evaluating" an array of functions

```jldoctest
using Gridap.Arrays

f = [+,-,max,min]
a = [1,2,3,4]
b = [4,3,2,1]

c = apply(f,a,b)

println(c)

# output
[5, -1, 3, 1]
```
"""
function apply(f::AbstractArray,a::AbstractArray...)
  AppliedArray(f,a...)
end

"""
    apply(::Type{T},f::AbstractArray,a::AbstractArray...) where T

Like [`apply(f::AbstractArray,a::AbstractArray...)`](@ref), but the user provides the element type
of the resulting array in order to circumvent type inference.
"""
function apply(::Type{T},f::AbstractArray,a::AbstractArray...) where T
  AppliedArray(T,f,a...)
end

"""
    apply_all(f::Tuple,a::AbstractArray...) -> Tuple

Numerically equivalent to

    tuple( ( apply(fi, a...) for fi in f)... )

# Examples

```jldoctest
using Gridap.Arrays

a = [1,2,3,4]
b = [4,3,2,1]

c = apply_all( (+,-), a, b)

# Equivalent to
# c = ( apply(+,a,b), apply(-,a,b) )

println(c)

# output
([5, 5, 5, 5], [-3, -1, 1, 3])

```
"""
function apply_all(f::Tuple,a::AbstractArray...)
  _apply_several(a,f...)
end

function _apply_several(a,f,g...)
  fa = apply(f,a...)
  ga = _apply_several(a,g...)
  (fa,ga...)
end

function _apply_several(a,f)
  fa = apply(f,a...)
  (fa,)
end

# Helpers

struct AppliedArray{T,N,F,G} <: AbstractArray{T,N}
  g::G
  f::F
  function AppliedArray(::Type{T},g::AbstractArray,f::AbstractArray...) where T
    G = typeof(g)
    F = typeof(f)
    f1, = f
    new{T,ndims(f1),F,G}(g,f)
  end
end

function AppliedArray(g::AbstractArray,f::AbstractArray...)
  gi = testitem(g) #Assumes that all kernels return the same type
  fi = testitems(f...)
  T = typeof(kernel_testitem(gi,fi...))
  AppliedArray(T,g,f...)
end

function uses_hash(::Type{<:AppliedArray})
  Val(true)
end

function array_cache(hash::Dict,a::AppliedArray)
    id = objectid(a)
    _cache = _array_cache(hash,a)
    if haskey(hash,id)
      cache = _get_stored_cache(_cache,hash,id)
    else
      cache = _cache
      hash[id] = cache
    end
    cache
end

#TODO
@static if VERSION >= v"1.1"
  @noinline function _get_stored_cache(cache::T,hash,id) where T
    c::T = hash[id]
    c
  end
else
  @noinline function _get_stored_cache(cache::T,hash,id) where T
    hash[id]
  end
end

function _array_cache(hash,a::AppliedArray)
  cg = array_cache(hash,a.g)
  gi = testitem(a.g)
  fi = testitems(a.f...)
  @notimplementedif ! all( map(isconcretetype,map(eltype,a.f)) )
  if ! (eltype(a.g)<:Function)
    @notimplementedif ! isconcretetype(eltype(a.g))
  end
  cf = array_caches(hash,a.f...)
  cgi = kernel_cache(gi,fi...)
  ai = kernel_testitem!(cgi,gi,fi...)
  i = -testitem(eachindex(a))
  e = Evaluation((i,),ai)
  c = (cg, cgi, cf)
  (c,e)
end

function testitem(a::AppliedArray)
  cg = array_cache(a.g)
  gi = testitem(a.g)
  fi = testitems(a.f...)
  kernel_testitem(gi,fi...)
end

function getindex!(cache,a::AppliedArray,i::Integer...)
  li = LinearIndices(a)
  getindex!(cache,a,li[i...])
end

function getindex!(cache,a::AppliedArray,i::Integer)
  _cached_getindex!(cache,a,(i,))
end

function getindex!(cache,a::AppliedArray,i::CartesianIndex)
  _cached_getindex!(cache,a,Tuple(i))
end

function _cached_getindex!(cache,a::AppliedArray,i::Tuple)
  c, e = cache
  v = e.fx
  if e.x != i
    v = _getindex!(c,a,i...)
    e.x = i
    e.fx = v
  end
   v
end

function _getindex!(cache,a::AppliedArray,i...)
  cg, cgi, cf = cache
  gi = getindex!(cg,a.g,i...)
  fi = getitems!(cf,a.f,i...)
  vi = apply_kernel!(cgi,gi,fi...)
  vi
end

function Base.getindex(a::AppliedArray,i...)
  ca = array_cache(a)
  getindex!(ca,a,i...)
end

function Base.IndexStyle(
  ::Type{AppliedArray{T,N,F,G}}) where {T,N,F,G}
  common_index_style(F)
end

function Base.size(a::AppliedArray)
  f, = a.f
  size(f)
end

function common_size(a::AbstractArray...)
  a1, = a
  c = all([length(a1) == length(ai) for ai in a])
  if !c
    error("Array sizes $(map(size,a)) are not compatible.")
  end
  l = length(a1)
  (l,)
end

#TODO not sure what to do with shape and index-style
function common_index_style(::Type{<:Tuple})
  IndexLinear()
end

mutable struct Evaluation{X,F} <: GridapType
  x::X
  fx::F
  function Evaluation(x::X,fx::F) where {X,F}
    new{X,F}(x,fx)
  end
end

# Particular implementations for Fill

function apply(f::Fill,a::Fill...)
  ai = getvalues(a...)
  r = apply_kernel(f.value,ai...)
  s = common_size(f,a...)
  Fill(r,s)
end

function apply(::Type{T},f::Fill,a::Fill...) where T
  apply(f,a...)
end

function getvalues(a::Fill,b::Fill...)
  ai = a.value
  bi = getvalues(b...)
  (ai,bi...)
end

function getvalues(a::Fill)
  ai = a.value
  (ai,)
end

function Base.sum(a::AppliedArray)
  cache = array_cache(a)
  _sum_array_cache(cache,a)
end

function _sum_array_cache(cache,a)
  r = zero(eltype(a))
  for i in eachindex(a)
    ai = getindex!(cache,a,i)
    r += ai
  end
  r
end

# TODO Think about iteration and sub-iteration

# Helper type to test caching of intermediate results

struct ArrayWithCounter{T,N,A,C} <: AbstractArray{T,N}
  array::A
  counter::C
  function ArrayWithCounter(a::AbstractArray{T,N}) where {T,N}
    c = zeros(Int,size(a))
    c[:] .= 0
    new{T,N,typeof(a),typeof(c)}(a,c)
  end
end

Base.size(a::ArrayWithCounter) = size(a.array)

function Base.getindex(a::ArrayWithCounter,i::Integer...)
  a.counter[i...] += 1
  a.array[i...]
end

Base.IndexStyle(::Type{<:ArrayWithCounter{T,N,A}}) where {T,A,N} = IndexStyle(A)

function reset_counter!(a::ArrayWithCounter)
  a.counter[:] .= 0
end


text:
using BinaryBuilder

name = "glog"
version = v"0.4.0"

sources = [
    ArchiveSource("https://github.com/google/glog/archive/v$(version).tar.gz",
                  "f28359aeba12f30d73d9e4711ef356dc842886968112162bc73002645139c39c")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glog-*

mkdir cmake_build
cd cmake_build/
FLAGS=()
if [[ "${target}" == i686-w64-mingw32 ]]; then
    # For some reasons, CMake can't find the symbol `UnDecorateSymbolName` in libdbghelp
    FLAGS+=(-DHAVE_DBGHELP=TRUE)
elif [[ "${target}" == *-freebsd* ]]; then
    # It looks like the build is broken with Clang for FreeBSD:
    #
    #   In file included from /workspace/srcdir/glog-0.4.0/src/utilities.cc:375:
    #   /workspace/srcdir/glog-0.4.0/src/stacktrace_x86_64-inl.h:64:6: error: use of undeclared identifier '_Unwind_Backtrace'
    #        _Unwind_Backtrace(nop_backtrace, NULL);
    #        ^
    #   /workspace/srcdir/glog-0.4.0/src/stacktrace_x86_64-inl.h:100:3: error: use of undeclared identifier '_Unwind_Backtrace'
    #     _Unwind_Backtrace(GetOneFrame, &targ);
    #     ^
    #   2 errors generated.
    #   make[2]: *** [CMakeFiles/glog.dir/build.make:115: CMakeFiles/glog.dir/src/utilities.cc.o] Error 1
    #
    CMAKE_TARGET_TOOLCHAIN=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake
fi

cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
      -DBUILD_SHARED_LIBS=ON \
      -DBUILD_TESTING=OFF \
      "${FLAGS[@]}" \
      ..

make -j${nproc}
make install

# Upstream PR from 30 Aug 2017 "Produce pkg-config file under cmake"
# See https://github.com/google/glog/pull/239 and #483, not released
mkdir ${prefix}/lib/pkgconfig
cat >${prefix}/lib/libglog.pc <<EOS
  prefix=${prefix}
  exec_prefix=${prefix}
  libdir=${libdir}
  includedir=${prefix}/include

  Name: libglog
  Description: Google Log (glog) C++ logging framework
  Version: 0.4.0
  Libs: -L${libdir} -lglog
  Cflags: -I${prefix}/include
EOS
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
# No products: Eigen is a pure header library
products = Product[
    LibraryProduct("libglog", :libglog),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using PCRP
using Test

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


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 374fe5f2-5731-11eb-1607-918d7e1432d5
TableOfContents(depth=10)

# ╔═╡ cea3f932-f935-11ea-0eb3-4f47e65da5fa
md"""
## $(@bind run_file_stats CheckBox()) File stats

Processes raw data in `DATA_DIR` and saves summary of header info to file
"""

# ╔═╡ 30a59602-3e46-11eb-083b-01049a15e927
if run_file_stats
	const DATA_DIR = "./Projects/HATP26b/data/ut190313"
end

# ╔═╡ 18536d9c-f922-11ea-02ed-5340b350c0c0
if run_file_stats
	let
		# Load
		df = fitscollection(
			DATA_DIR,
			abspath=false,
			exclude=fn"ift*[!c1].fits",
			exclude_key=("COMMENT"),
		)

		# Save
		header_keys = [
			"FILENAME",
			"UT-DATE",
			"UT-TIME",
			"UT-END",
			"EXPTIME",
			"OBJECT",
			"EXPTYPE",
			"RA",
			"DEC",
			"SLITMASK",
			"FILTER",
			"DISPERSR",
			"AIRMASS",
			"G-SEEING",
			"BINNING",
			"SPEED",

		]
		csv_path = dirname(DATA_DIR)
		date = basename(DATA_DIR)
		CSV.write(
			"$csv_path/logs/nightlog_$date.csv",
			df[!, header_keys]
		)

		# Preview
		df
	end
end

# ╔═╡ 5b08e224-3e4c-11eb-358e-8532b35c3151
md"""
## $(@bind run_raw_data_plot CheckBox()) Raw data plot

View summary data stored in `LOG_PATH`
"""

# ╔═╡ 4bf3a536-3e4e-11eb-394e-6326142d4746
if run_raw_data_plot
	const LOG_PATH = "./Projects/HATP26b/data/logs/nightlog_ut190313.csv"
end

# ╔═╡ d9d441d8-3e48-11eb-3dda-f334426e5b80
if run_raw_data_plot
	let
		key = :AIRMASS
		df = CSV.read(LOG_PATH, DataFrame)
		#df = @where(df, occursin.("sci", :OBJECT))
		plot(
			df[!, key],
			label = "some data",
			linecolor = :cyan,
			#markershape = :square,
			markercolor = :cyan,
			legend = :none, #:bottomright,
			xlabel = "Index",
			ylabel = key,
			#title = data_path,
		)
	end
end

# ╔═╡ 0110f78a-f945-11ea-004c-291cbfd9365c
md"""
## $(@bind run_tepspec CheckBox()) `tepspec` exploration

Plot divided LCs and chip traces from `tepspec` pickle stored under `RED_DIR`
"""

# ╔═╡ ac13cca4-3e35-11eb-083b-cba8e8b46c07
begin
	const TEP_PKL = "./Projects/HATP26b/data_reductions/ut190313_a15_25_noflat_LBR/LCs_hp26_bins_r.pkl"
	const DATA_RED = dirname(TEP_PKL)
end

# ╔═╡ 65530184-f968-11ea-3f55-7d8a470384ed
if run_tepspec
	data = load_pickle(TEP_PKL)
	flux_target = data["oLC"]
	flux_comps = data["cLC"]
	cNames = data["cNames"]
	with_terminal() do
		println("Keys: ", keys(data))
		println("\ncomp stars: ", cNames)
	end
end

# ╔═╡ 6f8188e8-f98f-11ea-0bd4-bd30f24b976e
if run_tepspec
	md"""
	**comp star index:**
	$(@bind idx_comp Slider(1:length(cNames), show_value=true))
	"""
end

# ╔═╡ 51666c34-39bd-11eb-0e2a-2371e207dce8
if run_tepspec
	scatter(
			flux_target ./ flux_comps[:, idx_comp],
			label = "target / $(cNames[idx_comp])",
			linecolor = :cyan,
			markershape = :circle,
			markercolor = :cyan,
			legend = :topright,
			xlabel = "Index",
			ylabel = "Flux",
			title = "Divided WLC",
	)
end

# ╔═╡ d6251ef2-fb4c-11ea-1ce3-97a485de4b95
if run_tepspec
	XX = load_pickle("$DATA_RED/XX.pkl")
    YY = load_pickle("$DATA_RED/YY.pkl")

	md"""
	**Number of traces to show for:**
	$(@bind trace_key Select(collect([Pair(key,key) for key in keys(XX)])))
	$(@bind num_traces Slider(2:7, show_value=true))
	"""
end

# ╔═╡ 6d2b41fc-fb96-11ea-0ca4-1d62d94d25d7
if run_tepspec
	trace_idxs = round.(Int, range(1, length(XX[trace_key]), length=num_traces))
	plot(
		XX[trace_key][trace_idxs],
		YY[trace_key][trace_idxs],
		label = reshape(trace_idxs, 1, :) .- 1,
		title = trace_key,
		legend = :left,
		legend_title = "trace 0-index",
		palette = palette(:diverging_isoluminant_cjo_70_c25_n256, num_traces),
	)
end

# ╔═╡ 93b51aba-fba4-11ea-2181-a9f87f9348e4
md"""
### $(@bind init_extracted_spectra CheckBox()) Intitial Extracted spectra

From spec.fits files
"""

# ╔═╡ 6b0c5d92-3bf7-11eb-03eb-871c9ffbd2e8
@bind extr_spec_key Select([
	"Wavelength",
	"Simple extracted object spectrum",
	"Simple extracted flat spectrum",
	"Pixel sensitivity (obtained by the flat)",
	"Simple extracted object spectrum/pixel sensitivity",
	"Sky flag (0 = note uneven sky, 1 = probably uneven profile,
		2 = Sky_Base failed)",
	"Optimally extracted object spectrum",
	"Optimally extracted object spectrum/pixel sensitivity",
		])

# ╔═╡ 25af8e98-fba8-11ea-2089-db810bb8778a
function plot_spectra!(p, data; key="Wavelength", label=nothing)
	i = Dict(
		"Wavelength" => 0,
		"Simple extracted object spectrum" => 1,
		"Simple extracted flat spectrum" => 2,
		"Pixel sensitivity (obtained by the flat)" => 3,
		"Simple extracted object spectrum/pixel sensitivity" => 4,
		"Sky flag (0 = note uneven sky, 1 = probably uneven profile,
		2 = Sky_Base failed)" => 5,
		"Optimally extracted object spectrum" => 6,
		"Optimally extracted object spectrum/pixel sensitivity" => 7,
	)[key]
	d = data[:, i+1, :]
	extracted_spectra_med = median(d, dims=2)
	σ = std(d, dims=2)
	plot!(
		p,
		extracted_spectra_med,
		ribbon = σ,
		label = label,
		legend = :topright,
	)
end

# ╔═╡ 9edfdb1e-fba4-11ea-2223-a5266f7ec2bb
if init_extracted_spectra
	let
		p = plot(xguide="Index", yguide="Value")
		for fpath in glob("$DATA_RED/*spec.fits")
			FITS(fpath) do f
				data = read(f[1])
				plot_spectra!(p, data, label=basename(fpath), key=extr_spec_key)
			end
		end
		p
	end
end

# ╔═╡ a8e8d9e2-3c1a-11eb-2bc3-bd1f05686d01
md"""
### $(@bind final_extracted_spectra CheckBox()) Final Extracted spectra

Loaded from tepspec pickle
"""

# ╔═╡ 07920dee-3c1e-11eb-2885-33fcf07a818f
if final_extracted_spectra
	let
		data = load_pickle(TEP_PKL)
		final_specs = data["optimal spectra"]
		wav_final_spec = final_specs["wavelengths"]
		p = plot()
		for (k, v) in final_specs
			if k != "wavelengths"
				flux_final_spec = vec(mean(v, dims=1))
				σ = vec(std(v, dims=1))
				if !any(isnan, flux_final_spec)
					plot!(
						p,
						wav_final_spec,
						flux_final_spec,
						ribbon = σ,
						label = k,
					)
				end
			end
		end
		p
	end
end

# ╔═╡ 767b67be-3714-11eb-2a8b-e13213320300
md"## Detrending"

# ╔═╡ bcda6c82-3e53-11eb-34f6-81e71a34a293
const GPT_DET_WLC = "./Projects/HATP26b/data_detrending/out_r/HATP26/hp26_190313_r/white-light/PCA_2/detrended_lc.dat"

# ╔═╡ 586f547e-3714-11eb-0fa1-cf1a1ebd5f50
md"""### $(@bind run_GPT_WLC CheckBox()) GPT WLC

Visualize GPT data stored in `GPT_DET_WLC`
"""

# ╔═╡ 72c7e8f4-3cc3-11eb-2f6a-07aff33cd8aa
if run_GPT_WLC
	let
		data = CSV.read(
			GPT_DET_WLC,
			DataFrame;
			comment = "#",
			header = ["Time", "DetFlux", "DetFluxErr", "Model"],
		)

		model = data.Model
		flux = data.DetFlux
		resid = (flux - model) * 1e6

		p_data = scatter(flux; label="data", yguide="Normalized flux")
		p_model_and_data = plot!(p_data, model; label="model", lw=4)
		p_resid = scatter(resid; label="residual", yguide="ppm")

		l = @layout [top; bottom{0.2h}]
		plot(p_model_and_data, p_resid; layout=l, link=:x, legend=nothing)
	end
end

# ╔═╡ e2387b52-572d-11eb-1460-0dc9ab7a1866
md"""### $(@bind run_GPT_WLC_corner CheckBox()) GPT WLC Corner Plot

Visualize GPT data stored in
"""

# ╔═╡ 03b4f076-3e56-11eb-133c-393802c43eb4
gpts_pkl_to_dict = load_pickle(
	"Projects/HATP23b/data_detrending/out_c/HATP23b/hp23b_160621_custom/white-light/BMA_posteriors.pkl"
);

# ╔═╡ 8b1bf7d4-56ee-11eb-2d79-e983c78df3ec
with_terminal() do
	println(keys(gpts_pkl_to_dict))
end

# ╔═╡ 6fa1fd52-573c-11eb-2597-d9ac3957e144
# The filter version below would be great if the dict keys were already keys, but they're not cause python =\
function sub_dict(dict, keys)
    Dict(
        Symbol(index)=>value for (index, value) in pairs(dict) if index ∈ keys
    )
end

# d_gpts = filter(
# 	p -> p.first ∈ ["p", "t0", "P", "rho", "inc", "b", "aRs", "q1"],
# 	d_pkl_gpts
# )

# ╔═╡ b1259430-56ee-11eb-3bf4-37e96bcc9e44
begin
	const gpts_params = ["p", "t0", "P", "rho", "inc", "b", "aRs", "q1"]
	gpts_params_dict = sub_dict(
		gpts_pkl_to_dict, gpts_params
	)
end

# ╔═╡ 6e4cf0ba-5728-11eb-289d-ebd928a3e9f1
gpts_params_dict[:t0] .-= 2.45485e6; gpts_params_dict

# ╔═╡ 7a7d7af0-5947-11eb-210b-cf4c787ed654
gpts_params_dict[:P] .-= 1.21288; gpts_params_dict

# ╔═╡ bd5fde4c-5705-11eb-3bd3-45befc4431ae
df_gpts = DataFrame(gpts_params_dict)#[1:50, :]

# ╔═╡ d829de14-56ef-11eb-16c2-c5ef0982e63f
function corner(
		df,
		vars=nothing;
		quantiles=[0.25,0.50,0.75],
		bandwidthx=100,
		bandwidthy=100,
		kwargs...,
	)

	# valid variables
	validvars = propertynames(df)

	# plot all variables by default
	isnothing(vars) && (vars = validvars)
	@assert vars ⊆ validvars "invalid variable name"

	plts = []
	n = length(vars)
	for i in 1:n, j in 1:n
		# Only label axes of first column and last row or plots
		xaxis = i == n
		yaxis = j == 1

		xticks = i == n
		xguide = i == n ? vars[j] : ""
		yticks = i > 1 && j == 1
		yguide = i > 1 && j == 1 ? vars[i] : ""

		# 1D histograms
		if i == j
			p = plot(xticks=xticks, yticks=yticks)
			# p = distplot1D(
			# 	df,
			# 	vars[i];
			# 	quantiles = quantiles,
			# 	xaxis = xaxis,
			# 	yaxis = yaxis,
			# 	xrotation = 45,
			# 	yrotation = 45,
			# 	title = vars[i],
			# 	# xticks=xticks,
			# 	# yticks=yticks,
			# 	# xguide=xguide,
			# 	# yguide=yguide,
			#)
		# 2D histograms
		elseif i > j
			p = plot(xticks=xticks, yticks=yticks)
			# p = distplot2D(
			# 	df,
			# 	vars[j],
			# 	vars[i],
			# 	xaxis = xaxis,
			# 	yaxis = yaxis,
			# 	xrotation = 45,
			# 	yrotation = 45,
			# 	# quantiles=quantiles,
			# 	# bandwidthx=bandwidthx,
			# 	# bandwidthy=bandwidthy,
			# 	# xticks=xticks,
			# 	# yticks=yticks,
			# 	# xguide=xguide,
			# 	# yguide=yguide,
			# )
		else
	  		p = plot(framestyle=:none)
		end

		push!(plts, p)
	end

	return plot(plts...; layout=(n, n), leg=:none, colorbar=:none, kwargs...)
end

# ╔═╡ c0960acc-5701-11eb-0183-7f3134eb4f1d
p = corner(df_gpts, Symbol.(gpts_params), size=(15_000, 30_000)); png("/home/mango/Desktop/test")

# ╔═╡ dc52d0fc-56f1-11eb-1f03-2d0c940b554c
function distplot1D(df, key; quantiles=[0.25, 0.50, 0.75], kwargs...)
	p = @df df stephist(cols(key), c=:blue, fill=true; kwargs...) # Histogram
	qs = @df df quantile(cols(key), quantiles) # Quantiles
	vline!(p, qs)
	return p
end

# ╔═╡ 7421e850-5738-11eb-0170-d92e031c6836
# distplot2D(df_gpts, "b", "aRs", fill=(true, :reds), levels=5)

# ╔═╡ c838cdb0-5700-11eb-1e17-b13f17cc64af
function distplot2D(df, k1, k2; kwargs...)
	dens = @df df_gpts kde((cols(Symbol(k1)), cols(Symbol(k2))))
	plot(dens; fill=(true, :Blues), levels=5, kwargs...)
end

# ╔═╡ 8fbbec74-f970-11ea-34d8-174a293a4107
md"### Photmetric monitoring"

# ╔═╡ 4494f844-f98b-11ea-330d-e3c18965972d
md"#### $(@bind run_phot_mon CheckBox()) ASAS-SN data"

# ╔═╡ b38c37b6-f970-11ea-033e-6d215c8b87df
if run_phot_mon
	df_phot = CSV.read(
		"projects/HATP23b/asas-sn_hatp23.csv",
	)
	HTMLTable(describe(df_phot))
end

# ╔═╡ 03cd4060-fab6-11ea-1298-ebbeccac9d09
function to_dt(s)
	# Converts y-m-d.x to y-m-dTh:m:s.x
	# Courtesy of @Yakir Luc Gagnon =]
	m = match(r"^(\d{4})-(\d\d)-(\d\d).(\d+)$", s)
	xs = parse.(Int, m.captures)
	x = pop!(xs)
	x /= exp10(ndigits(x))
	DateTime(xs...) + Millisecond(round(86_400_000*x))
end

# ╔═╡ 8041e34c-f971-11ea-0a9e-dd6415b3f33c
if run_phot_mon
	@df df_phot scatter(
		to_dt.(cols(Symbol("UT Date"))),
		cols(Symbol("flux(mJy)")),
		group = :Filter,
		legend = :bottomleft,
		xlabel = "UT Date",
		ylabel = "flux(mJy)",
	)

	transit_epochs = DateTime.(
		[
			"2016-06-22T08:18:00",
			"2017-06-10T07:05:00",
			"2018-06-04T07:24:00",
			"2018-06-21T06:56",
			"2018-08-22T03:30:00",
		]
	)
	vline!(transit_epochs, label="transits")
end

# ╔═╡ 7b2ef6f4-f9e3-11ea-22ed-0367935634be
# # Call Julia's copy of Python to load in pickles
# begin
# 	py"""
# 	import numpy as np
# 	import pickle

# 	def load_pickle(fpath):
# 		with open(fpath, "rb") as f:
# 			data = pickle.load(f)
# 		return data
# 	"""
# 	load_pickle = py"load_pickle"
# end

# ╔═╡ 6763632d-7a58-4d22-aeb7-b2242f5e920b
md"""
## Plot configs
"""

# ╔═╡ b34df520-e94e-4a29-ba42-1e265848819b
const pal = palette(["#f5d300", "#d55e00", "#08f7fe", "#0173b2", "#029e73"])

# ╔═╡ 6319a232-f8a3-11ea-0a58-7bf7af34ff4d
begin
	theme(:dark)
	default(
		titlefont = "Lato",
		guidefont = "Lato",
		markerstrokewidth = 0,
		palette=pal,
	)
end

# ╔═╡ 1497392c-f96d-11ea-3e29-15d8b60f4559
plotly()

# ╔═╡ 6b12e47b-74f0-414b-aa47-583e11c3ab10
md"""
## Packages
"""

# ╔═╡ d02112cc-5739-11eb-117f-1db80a9baafc
using PlotlyBase

# ╔═╡ de8e285c-f93c-11ea-1571-a5a5dbc48e3c
using Revise, PlutoUI, StatsPlots, CSV, CCDReduction, FITSIO, DataFrames, DataFramesMeta, Glob, AstroImages, HDF5, Dates, Printf, PyCall, RecursiveArrayTools, BenchmarkTools, Statistics, PaddedViews, PlotlyBase, KernelDensity

# ╔═╡ Cell order:
# ╟─374fe5f2-5731-11eb-1607-918d7e1432d5
# ╟─cea3f932-f935-11ea-0eb3-4f47e65da5fa
# ╠═30a59602-3e46-11eb-083b-01049a15e927
# ╟─18536d9c-f922-11ea-02ed-5340b350c0c0
# ╟─5b08e224-3e4c-11eb-358e-8532b35c3151
# ╠═4bf3a536-3e4e-11eb-394e-6326142d4746
# ╟─d9d441d8-3e48-11eb-3dda-f334426e5b80
# ╟─0110f78a-f945-11ea-004c-291cbfd9365c
# ╠═ac13cca4-3e35-11eb-083b-cba8e8b46c07
# ╟─65530184-f968-11ea-3f55-7d8a470384ed
# ╟─6f8188e8-f98f-11ea-0bd4-bd30f24b976e
# ╟─51666c34-39bd-11eb-0e2a-2371e207dce8
# ╟─d6251ef2-fb4c-11ea-1ce3-97a485de4b95
# ╠═6d2b41fc-fb96-11ea-0ca4-1d62d94d25d7
# ╟─93b51aba-fba4-11ea-2181-a9f87f9348e4
# ╟─6b0c5d92-3bf7-11eb-03eb-871c9ffbd2e8
# ╟─9edfdb1e-fba4-11ea-2223-a5266f7ec2bb
# ╠═25af8e98-fba8-11ea-2089-db810bb8778a
# ╟─a8e8d9e2-3c1a-11eb-2bc3-bd1f05686d01
# ╟─07920dee-3c1e-11eb-2885-33fcf07a818f
# ╟─767b67be-3714-11eb-2a8b-e13213320300
# ╠═bcda6c82-3e53-11eb-34f6-81e71a34a293
# ╟─586f547e-3714-11eb-0fa1-cf1a1ebd5f50
# ╟─72c7e8f4-3cc3-11eb-2f6a-07aff33cd8aa
# ╠═e2387b52-572d-11eb-1460-0dc9ab7a1866
# ╠═03b4f076-3e56-11eb-133c-393802c43eb4
# ╠═8b1bf7d4-56ee-11eb-2d79-e983c78df3ec
# ╠═b1259430-56ee-11eb-3bf4-37e96bcc9e44
# ╠═6e4cf0ba-5728-11eb-289d-ebd928a3e9f1
# ╠═7a7d7af0-5947-11eb-210b-cf4c787ed654
# ╠═6fa1fd52-573c-11eb-2597-d9ac3957e144
# ╠═bd5fde4c-5705-11eb-3bd3-45befc4431ae
# ╠═c0960acc-5701-11eb-0183-7f3134eb4f1d
# ╠═d829de14-56ef-11eb-16c2-c5ef0982e63f
# ╠═dc52d0fc-56f1-11eb-1f03-2d0c940b554c
# ╠═7421e850-5738-11eb-0170-d92e031c6836
# ╠═d02112cc-5739-11eb-117f-1db80a9baafc
# ╠═c838cdb0-5700-11eb-1e17-b13f17cc64af
# ╟─8fbbec74-f970-11ea-34d8-174a293a4107
# ╟─4494f844-f98b-11ea-330d-e3c18965972d
# ╠═b38c37b6-f970-11ea-033e-6d215c8b87df
# ╠═8041e34c-f971-11ea-0a9e-dd6415b3f33c
# ╟─03cd4060-fab6-11ea-1298-ebbeccac9d09
# ╠═7b2ef6f4-f9e3-11ea-22ed-0367935634be
# ╟─6763632d-7a58-4d22-aeb7-b2242f5e920b
# ╠═b34df520-e94e-4a29-ba42-1e265848819b
# ╠═6319a232-f8a3-11ea-0a58-7bf7af34ff4d
# ╠═1497392c-f96d-11ea-3e29-15d8b60f4559
# ╟─6b12e47b-74f0-414b-aa47-583e11c3ab10
# ╠═de8e285c-f93c-11ea-1571-a5a5dbc48e3c


text:
export sass_make_file_context, sass_delete_file_context, sass_compile_file_context
export sass_file_context_get_context, sass_file_context_get_options, sass_file_context_set_options
export sass_make_data_context, sass_delete_data_context, sass_compile_data_context
export sass_data_context_get_context, sass_data_context_get_options, sass_data_context_set_options
export sass_context_get_options, sass_context_set_options
export sass_context_get_output_string
export sass_context_get_source_map_string
export sass_context_get_error_status
export sass_context_get_error_json
export sass_context_get_error_text
export sass_context_get_error_message

# File

function sass_make_file_context(filename)
    ccall((:sass_make_file_context, libsass_so),
        Ptr{Cvoid}, (Cstring,), filename)
end

function sass_delete_file_context(context)
    ccall((:sass_delete_file_context, libsass_so),
        Cvoid, (Ptr{Cvoid},), context)
end

function sass_compile_file_context(context)
    ccall((:sass_compile_file_context, libsass_so),
        Cint, (Ptr{Cvoid},), context)
end

function sass_file_context_get_context(context)
    ccall((:sass_file_context_get_context, libsass_so),
        Ptr{Cvoid}, (Ptr{Cvoid},), context)
end

function sass_file_context_get_options(context)
    ccall((:sass_file_context_get_options, libsass_so),
        Ptr{Cvoid}, (Ptr{Cvoid},), context)
end

function sass_file_context_set_options(context, options)
    ccall((:sass_file_context_get_options, libsass_so),
        Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), context, options)
end

# Data

function sass_make_data_context(data)
    ccall((:sass_make_data_context, libsass_so),
        Ptr{Cvoid}, (Ptr{UInt8},), data)
end

function sass_delete_data_context(context)
    ccall((:sass_delete_data_context, libsass_so),
        Cvoid, (Ptr{Cvoid},), context)
end

function sass_compile_data_context(context)
    ccall((:sass_compile_data_context, libsass_so),
        Cint, (Ptr{Cvoid},), context)
end

function sass_data_context_get_context(context)
    ccall((:sass_data_context_get_context, libsass_so),
        Ptr{Cvoid}, (Ptr{Cvoid},), context)
end

function sass_data_context_get_options(context)
    ccall((:sass_data_context_get_options, libsass_so),
        Ptr{Cvoid}, (Ptr{Cvoid},), context)
end

function sass_data_context_set_options(context, options)
    ccall((:sass_data_context_get_options, libsass_so),
        Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), context, options)
end

# Context

function sass_context_get_options(context)
    ccall((:sass_context_get_options, libsass_so),
        Ptr{Cvoid}, (Ptr{Cvoid},), context)
end

function sass_context_set_options(context, options)
    ccall((:sass_context_get_options, libsass_so),
        Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), context, options)
end

function sass_context_get_output_string(context)
    s = ccall((:sass_context_get_output_string, libsass_so),
        Cstring, (Ptr{Cvoid},), context)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_context_get_source_map_string(context)
    s = ccall((:sass_context_get_source_map_string, libsass_so),
        Cstring, (Ptr{Cvoid},), context)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_context_get_error_status(context)
    ccall((:sass_context_get_error_status, libsass_so),
        Cint, (Ptr{Cvoid},), context)
end

function sass_context_get_error_json(context)
    s = ccall((:sass_context_get_error_json, libsass_so),
        Cstring, (Ptr{Cvoid},), context)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_context_get_error_text(context)
    s = ccall((:sass_context_get_error_text, libsass_so),
        Cstring, (Ptr{Cvoid},), context)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_context_get_error_message(context)
    s = ccall((:sass_context_get_error_message, libsass_so),
        Cstring, (Ptr{Cvoid},), context)
    s == C_NULL ? nothing : unsafe_string(s)
end


text:
components(sym::Symbol, n) =
    map(Base.OneTo(n)) do i
        Symbol(sym, subscripts[i])
    end

scalar(sym::Symbol, n::TupleN{Int}) =
    @variables($sym[Base.OneTo(prod(n))]) |> only |> collect

vector(sym::Symbol, n::TupleN{Int}) =
    scalar.(components(sym, length(n)), Ref(n))

vector(sym::Symbol, arr::TupleN{AbstractArray}) =
    scalar.(components(sym, length(arr)), arr)


text:
# create IOftuart

export FT_Location, FT_DeviceIndex, purge, UARTConfiguration, setlatencytimer

type IOftuart <: IO
  ft_handle  :: Culong
  readbuffer :: Array{UInt8,1}
  writebuffer :: Array{UInt8,1}
  rxbytesreturnedbuffer :: Array{Cuint,1}
  txbytesreturnedbuffer :: Array{Cuint,1}
  IOftuart(x) = new(x,Array(UInt8,1),Array(UInt8,1), Array(Cuint,1), Array(Cuint,1))
end

function Base.show(io::IO, x::IOftuart)
  print(io,"FTD2XX.IOftuart.ft_handle(")
  if handleis32
    @printf(io,"0x%08x",x.ft_handle)
  else
    @printf(io,"0x%16x",x.ft_handle)
  end
  println(io,")")
  return nothing
end

gethandle(io::IOftuart) = io.ft_handle

immutable FT_Location
  location :: Unsigned
end

immutable FT_DeviceIndex
  deviceindex :: Unsigned
end

immutable UARTConfiguration
  baud :: Cuint
  bits :: UInt8
  stop :: UInt8
  parity :: UInt8
  readtimeout :: Cuint
  writetimeout :: Cuint
  flowcontrol  :: Cushort
  xon :: UInt8
  xoff :: UInt8
  function UARTConfiguration(baud::Integer = 9600,
                             bits_user::Integer = 8,
                             stop_user::Integer = 1,
                             parity_user = 'n';
                             readtimeout::Integer = 100,  # in milliseconds
                             writetimeout::Integer = 100, # in milliseconds
                             flowcontrol = "none",
                             xon::UInt8 = 0x11,
                             xoff::UInt8 = 0x13)
    if bits_user == 8
      bits = FT_BITS_8
    elseif bits_user == 7
      bits = FT_BITS_7
    else
      error("ERROR: UART config bits must be 7 ot 8")
    end
    if stop_user == 1
      stop = FT_STOP_BITS_1
    elseif stop_user == 2
      stop = FT_STOP_BITS_2
    else
      error("ERROR: UART config stop bits must be 1 or 2")
    end
    upper_parity_user = uppercase(parity_user)
    if in(upper_parity_user,Set(['N',"N","NONE","FT_PARITY_NONE"]))
      parity = FT_PARITY_NONE
    elseif in(upper_parity_user,Set(['O',"O","ODD","FT_PARITY_ODD"]))
      parity = FT_PARITY_ODD
    elseif in(upper_parity_user,Set(['E',"E","EVEN","FT_PARITY_EVEN"]))
      parity = FT_PARITY_EVEN
    elseif in(upper_parity_user,Set(['M',"M","MARK","FT_PARITY_MARK"]))
      parity = FT_PARITY_MARK
    elseif in(upper_parity_user,Set(['S',"S","SPACE","FT_PARITY_SPACE"]))
      parity = FT_PARITY_SPACE
    else
      error("ERROR: UART config parity must be n, o, e, m or s")
    end
    upper_flowcontrol = uppercase(flowcontrol)
    if in(upper_flowcontrol,Set(['N',"N","NONE","FT_FLOW_NONE"]))
      flowcontrol = FT_FLOW_NONE
    elseif in(upper_flowcontrol,Set(['R',"R","RTS_CTS","FT_FLOW_RTS_CTS"]))
      flowcontrol = FT_FLOW_RTS_CTS
    elseif in(upper_flowcontrol,Set(['D',"D","DTR_DSR","FT_FLOW_DTR_DSR"]))
      flowcontrol = FT_FLOW_DTR_DSR
    elseif in(upper_flowcontrol,Set(['X',"X","XON_XOFF","FT_FLOW_XON_XOFF"]))
      flowcontrol = FT_FLOW_XON_XOFF
    else
      error("ERROR: UART config flow control must be n, r, d or x")
    end
    new(baud, bits, stop, parity, readtimeout, writetimeout, flowcontrol, xon, xoff)
  end
end

function configureuart(io::IOftuart, config :: UARTConfiguration)
  ft_setbitmode(io.ft_handle, 0x00, FT_BITMODE_RESET)  #  make sure configured for UART
  ft_setbaudrate(io.ft_handle, config.baud)
  ft_setdatacharacteristics(io.ft_handle, config.bits, config.stop, config.stop)
  ft_setflowcontrol(io.ft_handle, config.flowcontrol, config.xon, config.xoff)
  ft_settimeouts(io.ft_handle, config.readtimeout, config.writetimeout)
  return nothing
end

function Base.open(device::FT_DeviceIndex, config::UARTConfiguration)
  io = IOftuart(ft_open(device.deviceindex))
  configureuart(io, config)
  finalizer(io,close)
  return io
end

function Base.open(device::FT_Location, config::UARTConfiguration)
  io = IOftuart(ft_openex(device.location))
  configureuart(io, config)
  finalizer(io,close)
  return io
end

function Base.open(device::Union{FT_SerialNumber,FT_Description}, config::UARTConfiguration)
  io = IOftuart(ft_openex(device))
  configureuart(io, config)
  finalizer(io,close)
  return io
end

function Base.open(device::Union{FT_DeviceIndex, FT_Location, FT_SerialNumber, FT_Description},
                   baud::Integer = 9600,
                   bits_user::Integer = 8,
                   stop_user::Integer = 1,
                   parity_user = 'n')
  open(device, UARTConfiguration(baud, bits_user, stop_user, parity_user))
end

function Base.close(io::IOftuart)
  if io.ft_handle != 0
    ft_close(io.ft_handle)
    io.ft_handle = 0
  end
  return nothing
end

Base.isopen(io::IOftuart) = io.ft_handle!=0
Base.isreadable(io::IOftuart) = isopen(io)
Base.iswritable(io::IOftuart) = isopen(io)

function ft_read!(io::IOftuart,
                 buffer::Array{UInt8,1},
                 bytestoread::Integer = length(buffer))
  ft_status = ccall((:FT_Read, d2xx),
                      Cuint,
                      (Culong, Ptr{UInt8}, Cuint, Ptr{Cuint}),
                      io.ft_handle, buffer, bytestoread, io.rxbytesreturnedbuffer)
  checkstatus(ft_status)
  return (io.rxbytesreturnedbuffer[1])
end

function Base.read(s::IOftuart, ::Type{UInt8})
  ft_read!(s, s.readbuffer, 1)
  return s.readbuffer[1]
end

function Base.read!(s::IOftuart, a::Vector{UInt8})
  bytesread = ft_read!(s, a)
  return a  # IO retruns the data not bytes read?
end

function ft_write(io::IOftuart,
                  buffer::Array{UInt8,1},
                  bytestowrite::Integer = length(buffer))
  ft_status = ccall((:FT_Write, d2xx),
                     Cuint,
                     (Culong, Ptr{UInt8}, Cuint, Ptr{Cuint}),
                     io.ft_handle, buffer, bytestowrite, io.txbytesreturnedbuffer)
  checkstatus(ft_status)
  return io.txbytesreturnedbuffer[1]
end

function Base.write(s::IOftuart, x::UInt8)
  s.writebuffer[1] = x
  byteswritten = ft_write(s.ft_handle, s.writebuffer, 1)
end

function Base.write(s::IOftuart, x::Vector{UInt8})
  byteswritten = ft_write(s, x, length(x))
end

function purge(io::IOftuart)
  ft_purge(io.ft_handle, FT_PURGE_RX | FT_PURGE_TX)
end

function setlatencytimer(io::IOftuart, latencyms::Integer)
  ft_setlatencytimer(io.ft_handle, latencyms)
end

Base.eof(io::IOftuart) = false


text:
function layerMean( x :: Array{Float64, 3}, region = (1,2))
   return squeeze( mean(x, region), region)
end

export layerMean

function layermixing( x :: Array{Float64, 3}, cutoff :: Float64)
  h = zeros(Int64, size(x,1,2))
  for ix = 1:size(h,1)
    for iy = 1:size(h,2)
      h[ix,iy] = findlast( x[ix,iy,:]  .< cutoff)
    end
  end
  return h
end
export layermixing


text:
using GeomDAE
using Base.Test

# explicit midpoint
a = Rational{Int64}[[0    0]
                    [1//2 0]]
b = Rational{Int64}[0, 1   ]
c = Rational{Int64}[0, 1//2]
o = 2

@test TableauERK(:explicit_midpoint, o, a, b, c) == getTableauExplicitMidpoint()

tab_explicit_midpoint1 = TableauERK(:explicit_midpoint, o, a, b, c)
tab_explicit_midpoint2 = TableauDIRK(:explicit_midpoint, o, a, b, c)
tab_explicit_midpoint3 = TableauFIRK(:explicit_midpoint, o, a, b, c)

@test tab_explicit_midpoint1 == tab_explicit_midpoint2
@test tab_explicit_midpoint2 == tab_explicit_midpoint3
@test tab_explicit_midpoint3 == tab_explicit_midpoint1

@test !isequal(tab_explicit_midpoint1, tab_explicit_midpoint2)
@test !isequal(tab_explicit_midpoint2, tab_explicit_midpoint3)
@test !isequal(tab_explicit_midpoint3, tab_explicit_midpoint1)

tab_explicit_midpoint1 = TableauERK(:explicit_midpoint, o, a, b, c)
tab_explicit_midpoint2 = TableauERK(:explicit_midpoint, o, a, b, c)
@test tab_explicit_midpoint1 == tab_explicit_midpoint2
@test isequal(tab_explicit_midpoint1, tab_explicit_midpoint2)

tmp = mktempdir()
tab_explicit_midpoint1 = TableauERK(:explicit_midpoint, o, a, b, c)
writeTableauToFile(tmp, tab_explicit_midpoint1)
tab_explicit_midpoint2 = readTableauERKFromFile(tmp, "explicit_midpoint")
rm(tmp, recursive=true)
@test tab_explicit_midpoint1 == tab_explicit_midpoint2
@test !isequal(tab_explicit_midpoint1, tab_explicit_midpoint2) # fails because data arrays in tab_explicit_midpoint1 are of type Rational{Int64} while data arrays in tab_explicit_midpoint2 are of type Float

# Heun
a = [[0.0 0.0]
     [1.0 0.0]]
b = [0.5, 0.5]
c = [0.0, 1.0]
o = 2

tab_explicit_heun = TableauERK(:heun, o, a, b, c)
tab_explicit_heun = TableauDIRK(:heun, o, a, b, c)
tab_explicit_heun = TableauFIRK(:heun, o, a, b, c)

@test TableauERK(:heun, o, a, b, c) == getTableauHeun()

tmp = mktempdir()
tab_explicit_heun1 = TableauERK(:heun, o, a, b, c)
writeTableauToFile(tmp, tab_explicit_heun1)
tab_explicit_heun2 = readTableauERKFromFile(tmp, "heun")
rm(tmp, recursive=true)
@test tab_explicit_heun1 == tab_explicit_heun2
@test isequal(tab_explicit_heun1, tab_explicit_heun2)

# Crouzeix
fac = 0.5/√3
a = [[ 0.5+fac 0.0    ]
     [-2.0*fac 0.5+fac]]
b = [0.5,     0.5    ]
c = [0.5+fac, 0.5-fac]
o = 2

@test TableauDIRK(:crouzeix, o, a, b, c) == getTableauCrouzeix()

@test_throws AssertionError tab_explicit_crouzeix = TableauERK(:crouzeix, o, a, b, c)
tab_explicit_crouzeix = TableauDIRK(:crouzeix, o, a, b, c)
tab_explicit_crouzeix = TableauFIRK(:crouzeix, o, a, b, c)

# implicit midpoint
a = ones(Rational{Int64}, 1, 1) // 2
b = Rational{Int64}[1   ]
c = Rational{Int64}[1//2]
o = 2

@test TableauFIRK(:implicit_midpoint, o, a, b, c) == getTableauImplicitMidpoint()
@test TableauFIRK(:implicit_midpoint, o, a, b, c) == getTableauGLRK1()

@test_throws AssertionError tab_implicit_midpoint = TableauERK(:implicit_midpoint, o, a, b, c)
@test_throws AssertionError tab_implicit_midpoint = TableauDIRK(:implicit_midpoint, o, a, b, c)
tab_implicit_midpoint = TableauFIRK(:implicit_midpoint, o, a, b, c)

# Gauss-Legendre
a = [[0.25       0.25-√3/6]
     [0.25+√3/6  0.25     ]]
b = [0.5,      0.5     ]
c = [0.5-√3/6, 0.5+√3/6]
o = 4

@test TableauFIRK(:glrk2, o, a, b, c) == getTableauGLRK2()

@test_throws AssertionError tab_explicit_glrk2 = TableauERK(:glrk2, o, a, b, c)
@test_throws AssertionError tab_explicit_glrk2 = TableauDIRK(:glrk2, o, a, b, c)
tab_explicit_glrk2 = TableauFIRK(:glrk2, o, a, b, c)

# instantiate all explicit tableaus
@test typeof(getTableauExplicitEuler()) <: TableauERK
@test typeof(getTableauExplicitMidpoint()) <: TableauERK
@test typeof(getTableauHeun()) <: TableauERK
@test typeof(getTableauKutta()) <: TableauERK
@test typeof(getTableauERK4()) <: TableauERK
@test typeof(getTableauERK438()) <: TableauERK

# instantiate all diagonally implicit tableaus
@test typeof(getTableauCrouzeix()) <: TableauDIRK

# instantiate all fully implicit tableaus
@test typeof(getTableauImplicitEuler()) <: TableauFIRK
@test typeof(getTableauImplicitMidpoint()) <: TableauFIRK
@test typeof(getTableauGLRK1()) <: TableauFIRK
@test typeof(getTableauGLRK2()) <: TableauFIRK
@test typeof(getTableauGLRK3()) <: TableauFIRK

# instatiate all partitioned tableaus
@test typeof(getTableauSymplecticEulerA()) <: TableauEPRK
@test typeof(getTableauSymplecticEulerB()) <: TableauEPRK

# test instatiation of partioned tableau by composition of two RK tableaus
@test typeof(TableauEPRK(:PERK4, 4, getTableauERK4(), getTableauERK4())) <: TableauEPRK

# TODO Add tests for TableauIPRK, TableauSARK, TableauSPARK and TableauGLM.


text:
module LovoAlign

  using Random
  include("./PDB.jl")

  #export Pair
  #struct Pair
  #  p1 :: String
  #  p2 :: String
  #  TMScore :: Float64
  #end
  #Pair(;p1=p1,p2=p2) = Pair(p1,p2)

  #
  # Performs a alignment from coordinates, not pdb names
  #

  function align(x1,x2;
                 method=2,
                 lovoalign_exec="lovoalign",debug=false)

    rnd_chars=Random.randstring(8)

    # Write temporary PDB file to lovoalign

    n1 = size(x1)[1]
    file1_name = "LovoAlign_1_$(rnd_chars).pdb"
    file = open(file1_name,"w")
    for i in 1:n1
      line = PDB.printatom(i,"CA","GLY","A",i,x1[i,1],x1[i,2],x1[i,3])
      write(file,line,"\n")
    end
    close(file)

    n2 = size(x2)[1]
    file2_name = "LovoAlign_2_$(rnd_chars).pdb"
    file = open(file2_name,"w")
    for i in 1:n2
      line = PDB.printatom(i,"CA","GLY","A",i,x2[i,1],x2[i,2],x2[i,3])
      write(file,line,"\n")
    end
    close(file)

    if n1 < 15 || n2 < 15
      method = 3
    end

    output="LovoAlign_3_$(rnd_chars).pdb"
    lovoalign_output=read(
     `$lovoalign_exec -p1 $file1_name -p2 $file2_name -m $method -o $output`,
     String)

    if debug
      println(lovoalign_output)
    end

    y = Array{Float32}(undef,n1,3)
    file = open(output,"r")
    i = 0
    for line in eachline(file)
      if line[1:4] == "ATOM"
        i = i + 1
        y[i,1] = parse(Float64,line[31:38])
        y[i,2] = parse(Float64,line[39:46])
        y[i,3] = parse(Float64,line[47:54])
      end
    end
    close(file)

   run(`\rm -f $file1_name $file2_name $output`)

   return y

  end

  #
  # Perform a general, non-sequential alignment, using the Non-Bijective method
  #

  function nonseq(pdb1,pdb2,sel1,sel2,output;
                  vmd_exec="vmd",
                  lovoalign_exec="lovoalign",
                  debug=false)

    checkpdbs(pdb1,pdb2,output)

    # Define first selection using VMD
    sel1_pdb = set_selection_pdb(pdb1,sel1,vmd_exec=vmd_exec)

    # Define second selection using VMD
    sel2_pdb = set_selection_pdb(pdb2,sel2,vmd_exec=vmd_exec,debug=debug)

    # Run LovoAlign for general non-sequential alignment
    lovoalign_output=read(
     `$lovoalign_exec -p1 $sel1_pdb -p2 $sel2_pdb -g 0. -beta1 -beta2 -m 3 -all -nglobal 50 -dtri 10. -o $output`,
     String)

    # Remove temporary files
    run(`\rm -f $sel1_pdb`)
    run(`\rm -f $sel2_pdb`)

    if debug
      println(lovoalign_output)
      return
    else
      println("Wrote $output file with $pdb1 aligned to $pdb2")
      return
    end

  end

  #
  # Perform a maximization of the TM-score
  #

  function TMScore(pdb1,pdb2,sel1,sel2,output;
                   vmd_exec="vmd",
                   lovoalign_exec="lovoalign",
                   debug=false)

    checkpdbs(pdb1,pdb2,output)

    # Define first selection using VMD
    sel1_pdb = set_selection_pdb(pdb1,sel1,vmd_exec=vmd_exec)

    # Define second selection using VMD
    sel2_pdb = set_selection_pdb(pdb2,sel2,vmd_exec=vmd_exec,debug=debug)

    # Run LovoAlign using default (TMscore) computation
    lovoalign_output=read(
     `$lovoalign_exec -p1 $sel1_pdb -p2 $sel2_pdb -o $output`,
     String)

    # Remove temporary files
    run(`\rm -f $sel1_pdb`)
    run(`\rm -f $sel2_pdb`)

    if debug
      println(lovoalign_output)
      return
    else
      println("Wrote $output file with $pdb1 aligned to $pdb2")
      return
    end

  end

  #
  # Function that checks some properteis of the PDB files
  #

  function checkpdbs(pdb1,pdb2,output)

    if pdb1 == output || pdb2 == output
      error("ERROR: Output file has the same name of one of input files.")
    end

  end

  #
  # Function that calls VMD to define the selected atoms
  #

  function set_selection_pdb(pdb,sel;vmd_exec="vmd",debug=false)

    rnd_chars=Random.randstring(8)
    filename = "./LovoAlign__$(rnd_chars)_$pdb"
    vmd_input_name = "./LovoAlign_$(rnd_chars).vmd"
    vmd_input = open(vmd_input_name,"w")
    write(vmd_input,"""
                    mol new \"$pdb\"
                    set all [ atomselect top all ]
                    \$all set beta 0
                    set sel [ atomselect top \"$sel\" ]
                    \$sel set beta 1
                    \$all writepdb $filename
                    \$all delete
                    \$sel delete
                    exit
                    """)
    close(vmd_input)
    vmd_output = read(`$vmd_exec -dispdev text -e $vmd_input_name`,String)
    if debug
      println(vmd_output)
    end
    run(`\rm -f $vmd_input_name`)
    return filename

  end

end


text:
module Advection

export
    div_Uu, div_Uv, div_Uw, div_Uc,
    U_grad_u, U_grad_v, U_grad_w, U_grad_c,

    momentum_flux_uu,
    momentum_flux_uv,
    momentum_flux_uw,
    momentum_flux_vu,
    momentum_flux_vv,
    momentum_flux_vw,
    momentum_flux_wu,
    momentum_flux_wv,
    momentum_flux_ww,
    advective_tracer_flux_x,
    advective_tracer_flux_y,
    advective_tracer_flux_z,

    CenteredSecondOrder,
    UpwindBiasedThirdOrder,
    UpwindBiasedFifthOrder,
    CenteredFourthOrder,
    WENO5

using Oceananigans.Grids
using Oceananigans.Operators

import Oceananigans.Grids: required_halo_size

abstract type AbstractAdvectionScheme end
abstract type AbstractCenteredAdvectionScheme <: AbstractAdvectionScheme end
abstract type AbstractUpwindBiasedAdvectionScheme <: AbstractAdvectionScheme end

required_halo_size(scheme::AbstractAdvectionScheme) = boundary_buffer(scheme) + 1

include("topologically_conditional_interpolation.jl")

include("centered_advective_fluxes.jl")
include("upwind_biased_advective_fluxes.jl")
include("flat_advective_fluxes.jl")

include("centered_second_order.jl")
include("upwind_biased_third_order.jl")
include("centered_fourth_order.jl")
include("upwind_biased_fifth_order.jl")
include("weno_fifth_order.jl")

include("momentum_advection_operators.jl")
include("tracer_advection_operators.jl")

end # module


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using LinearAlgebra
using SparseArrays
using DifferentialEquations
using ProgressMeter

include("utils\\randutils.jl")
include("utils\\categorical_tree.jl")

"""
    contact_process_gillespie(g, Xi, β; kwargs)

Simulate SI spreading on network `g` using the Gillespie algorithm.

`g` is an `AbstractSimpleGraph{<:Integer}`, `β` is the infection rate and `Xi` is a bit vector whose true entries denote initially infected nodes. Additional arguments are
  * `nmax=1000`: the maximum number of iterations of the algorithm
  * `tmax=100.0`: the maximum time allowed
  * `sampling_method=:tree`: The algorithm for sampling which node gets infected at a given time step. Choose from `:array`, `:sparse` or `:tree`. Default is `:tree` which scales better with the number of nodes (log(N) vs. N^2).
  * `record=false`: Whether or not to record which node gets infected and by whom. If `true`, return a vector of pairs `(k,j)` describing which node `k` got infected and `j` its neighbor that infected it.
"""
function contact_process_gillespie(g::AbstractSimpleGraph{<:Integer}, Xi::BitVector, β::Real; nmax=length(Xi), tmax=100, sampling_method=:tree, record=false)
    N = nv(g)
    X = copy(Xi)
    #nreactions::Int64 = 0
    if sampling_method == :tree
        a = zeros(N)
    elseif sampling_method == :sparse
        a = spzeros(N)
    else
        a = zeros(N)
    end
    for k in 1:N
        if !X[k]
            l = length(filter(i->X[i], inneighbors(g,k)))
            a[k] = β*l# probability of k recieving infection is β time its number of infected incoming neighbors
            #nreactions += l
        end
    end
    if sampling_method == :tree # use a specialized binary tree
        a = CategoricalTree(a)
    end
    a0 = sum(a)
    #a0 = nreactions*β
    if record
        Es = Vector{Tuple{Int64,Int64}}(undef, nmax-1) # which node gets infected, and by whom
    end
    ts = zeros(nmax)
    t = 0.0
    n = 1
    while n < nmax && t < tmax && a0 > 0
        τ = log(1/rand())/a0
        k = rand_categorical(a, a0) # which node gets infected
        X[k] = true
        #nreactions -= Int(floor(a[k]/β))
        a[k] = 0.0
        for i in Iterators.filter(j->!X[j],outneighbors(g,k))
            a[i] += β
            #nreactions += 1
        end
        if sampling_method == :sparse
            SparseArrays.dropstored!(a,k) # O(1) cost?
        end
        #a0 = β*nreactions
        a0 = sum(a)
        t += τ
        ts[n+1] = t
        if record
            j = rand(filter(l->X[l],inneighbors(g,k)))
            Es[n] = (k,j)
        end
        n += 1
    end
    if record
        return ts[1:n], sum(Xi) .+ collect(0:n-1), Es[1:(n-1)]
    else
        return ts[1:n], sum(Xi) .+ collect(0:n-1)
    end
end

"""
    contact_process_ode(g, Xi, β; kwargs)

Integrate the mean field approximation for SI spreading on network `g`.

`Xi` is the initial condition, a vector whose elements denote infected nodes by 1 and susceptible by 0. `β` is the infection rate.

Keyword arguments:
  * `tmax = 100.0`: maximum time of the simulation.
"""
function contact_process_ode(g::AbstractSimpleGraph{<:Integer}, Xi, β::Real; tmax=100.0, saveat=[])
    A = adjacency_matrix(g)
    u0 = float(Xi)
    f! = function(du,u,p,t)
        du .= β*(1.0 .- u).*(A*u)
    end
    prob = ODEProblem(f!, u0, (0.0,tmax))
    sol = solve(prob, Tsit5(), saveat=saveat)
    return  sol.t, hcat(sol.u...)', sol
end

"""
    contact_process_montecarlo(g, Xi, β; kwargs)

Compute the average and variance of trajectories given by the gillespie algorithm over `nsims` simulations.

Keyword arguments:
  * `nmax=1000`: maximum number of iterations of each simulation
  * `tmax=100.0`: maximum time of the simulations
  * `nsims=100`: number of simulations to perform
  * `nbins`: the number of time steps at which the average is computed
"""
function contact_process_montecarlo(g::AbstractSimpleGraph{<:Integer}, Xi::BitVector, β::Real; nmax = length(Xi), tmax = 100.0, nsims=1000, nbins = 100)
    ts = LinRange(0.0, tmax, nbins)
    X_sum = zeros(nbins)
    M_sum = zeros(nbins) # running SSE
    p = Progress(nsims, dt=1.0)
    for n in 1:nsims
        t, X = contact_process_gillespie(g, Xi, β, nmax=nmax, tmax=tmax)
        l = 1
        for k in 1:nbins
            while l < length(t) && t[l+1] < ts[k]
                l += 1
            end
            Xn = X[l]
            Δn = Xn - X_sum[k]/(n>1 ? n-1 : 1)
            X_sum[k] += Xn
            M_sum[k] += Δn*(Xn-X_sum[k]/n)
        end
        ProgressMeter.next!(p; showvalues= [(:n,n)])
    end
    return ts, X_sum/nsims, sqrt.(M_sum/(nsims-1))
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
using LinearAlgebraicRepresentation, ViewerGL
using BenchmarkTools
using Distributed
using TimerOutputs

include("../../src/AlphaShapes.jl")

Lar = LinearAlgebraicRepresentation
GL = ViewerGL

filename = "./examples/examples3D/OBJ/teapot.obj";
W, EVs, FVs = Lar.obj2lar(filename);
WW = [[i] for i = 1:size(W, 2)];
V, VV = Lar.apply(Lar.r(pi / 2, 0, 0), (W, WW)); #object rotated

points = convert(Lar.Points, V')
GL.VIEW([
    GL.GLPoints(points)
    GL.GLAxis(GL.Point3d(-1, -1, -1), GL.Point3d(1, 1, 1))
]);

V, VV = Lar.apply(Lar.r(pi / 2, 0, 0), (W, WW)); #object rotated
@btime AlphaShapes.alphaFilter(V);
filtration = AlphaShapes.alphaFilter(V);
VV, EV, FV, TV = AlphaShapes.alphaSimplex(V, filtration, 3.7)

GL.VIEW([
    GL.GLGrid(V, EV, GL.COLORS[1], 0.6) # White
    GL.GLGrid(V, FV, GL.COLORS[2], 0.3) # Red
    GL.GLGrid(V, TV, GL.COLORS[3], 0.3) # Green
]);

filter_key = sort(unique(values(filtration)))

granular = 10

reduced_filter =
    filter_key[sort(abs.(rand(Int, granular) .% length(filter_key)))]
reduced_filter = [reduced_filter; max(filter_key...)]
α=0.0

for α in reduced_filter
    @show α
    VVV, EEV, FFV, TTV = AlphaShapes.alphaSimplex(V, filtration, α)

    # GL.VIEW(GL.GLExplode(
    #     V,
    #     [[[t] for t in TTV]; [[f] for f in FFV]; [[e] for e in EEV]],
    #     1.0,
    #     1.0,
    #     1.0,# Explode Ratio
    #     99,
    #     1,# Colors
    # ))

end

print_timer(AlphaShapes.to);


text:
module SegmentIntersections

include("point.jl")
include("segment.jl")
include("bentley_ottmann/queue.jl")
include("bentley_ottmann/event.jl")
include("bentley_ottmann/status.jl")
include("bentley_ottmann/base.jl")
include("brute.jl")

end


text:
export iauPv2p
"""
Discard velocity component of a pv-vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   pv      double[2][3]     pv-vector

Returned:
   p       double[3]        p-vector

Called:
   iauCp        copy p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPv2p(double pv[2][3], double p[3])

function iauPv2p(pv::AbstractMatrix{<:Real})
   # Preallocate return values
   p = zeros(Float64, 3)

   ccall((:iauPv2p, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
          convert(Array{Float64, 2}, pv'), p)

   return SVector{3}(p)
end


text:
@testset "conditional_probability_rules" begin
    @syms a b c d

    @eqtest M.identify(P(a + b + c, d)) == P(a, d) + P(b, d) + P(c, d)
    @eqtest M.identify(P(a + b, (c + d))) == P(a, (c + d)) + P(b, (c + d))
end


text:
using EnergyExpressions
using AtomicLevels
using Test

include("conjugate_orbitals.jl")
include("bit_configurations.jl")
include("minors.jl")


text:
module MultiHZ

using JuLIP
using JuLIP.Potentials: pairs
import JuLIP: energy, forces, cutoff

export Multi_HZ

@pot struct MHZ{T} <: AbstractCalculator
   Z2idx::Vector{Int}
   V::Matrix{T}
   rcut::Float64
end

cutoff(V::MHZ) = V.rcut
# k is the energy scale, r0 is the sum of the two contacted grain
Hertz(k, α, r0) = @analytic r -> k * (r - r0)^α

function Multi_HZ(Z, k, α, r0; rcutfact = 2.5)
    @assert length(Z) == length(unique(Z))
    # create a mapping from atomic numbers to indices
    Z2idx = zeros(Int, maximum(Z))
    Z2idx[Z] = 1:length(Z)
    # generate the potential
    rcut = rcutfact*maximum(r0) # rcut is 5*R_{max}, used when searching possible pairs
    V = [ Hertz(k[a,b], α[a,b], r0[a,b]) * HS(r0[a,b])
        for a = 1:length(Z), b = 1:length(Z)]
    return MHZ(Z2idx, V, rcut)
end

function energy(V::MHZ, at::Atoms)
   E = 0.0
   for (i, j, r, R) in pairs(at, cutoff(V))
      a = V.Z2idx[at.Z[i]]
      b = V.Z2idx[at.Z[j]]
      E += 0.5 * V.V[a, b](r)
   end
   return E
end

function forces(V::MHZ, at::Atoms)
   F = zeros(JVecF, length(at))
   for (i, j, r, R) in pairs(at, cutoff(V))
      a = V.Z2idx[at.Z[i]]
      b = V.Z2idx[at.Z[j]]
      f = 0.5 * grad(V.V[a,b], r, R)
      F[i] += f
      F[j] -= f
   end
   return F
end

end


text:
module TestExpressionsLift

import Test: @testset, @test
import Volcanito: lift_function_calls

# TODO: Test the proper behavior when it's implemented.
test_cases = (
    (
        input=:(f(x)),
        output=:(f(x)),
    ),
)

@testset "lift" begin
    for case in test_cases
        @test lift_function_calls(case.input) == case.output
    end
end

end


text:
using MPI
using P4est_wrapper
using Test

# Initialize MPI if not initialized yet
if !MPI.Initialized()
    MPI.Init()
end

#############################################################################
# Definition of data structures and function callbacks
#############################################################################
mpicomm = P4est_wrapper.P4EST_ENABLE_MPI ? MPI.COMM_WORLD : Cint(0)

# Dummy callback
dummy_callback( ::Ptr{p8est_t}, which_tree::p4est_topidx_t, quadrant::Ptr{p8est_quadrant_t}) = Cint(0)

# C-callable dummy callback
const dummy_callback_c = @cfunction(dummy_callback, Cint, (Ptr{p8est_t}, p4est_topidx_t, Ptr{p8est_quadrant_t}))

#############################################################################
# Main program
#############################################################################

# Create a connectivity structure for the unit square.
unitcube_connectivity = p8est_connectivity_new_unitcube()
@test unitcube_connectivity != C_NULL
@test p8est_connectivity_is_valid(unitcube_connectivity) == 1
@test p8est_connectivity_complete(unitcube_connectivity) == nothing
@test p8est_connectivity_reduce(unitcube_connectivity) == nothing

# Test new connectivities
for test_conn in [p8est_connectivity_new_periodic(),
                p8est_connectivity_new_twocubes(),
                p8est_connectivity_new_twowrap(),
                p8est_connectivity_new_rotcubes(),
                p8est_connectivity_new_shell(),
                p8est_connectivity_new_sphere(),
                p8est_connectivity_new_brick(3, 2, 8, 1, 0, 1)]
    @test test_conn != C_NULL
    @test p8est_connectivity_is_valid(test_conn) == 1
    @test p8est_connectivity_is_equivalent(unitcube_connectivity, test_conn) == 0
    @test p8est_connectivity_destroy(test_conn) == nothing
end

# Create a new forest
unitcube_forest = p8est_new(mpicomm, unitcube_connectivity, 0, C_NULL, C_NULL)
@test unitcube_forest != C_NULL

# Create/Destroy geometry
unitcube_geometry = p8est_geometry_new_connectivity(unitcube_forest)
@test unitcube_geometry != C_NULL
@test p8est_geometry_destroy(unitcube_geometry) == nothing

# Copy/Destroy forest
unitcube_forest_copy =  p8est_copy(unitcube_forest, 1)
@test p8est_destroy(unitcube_forest_copy) == nothing
unitcube_forest_copy =  p8est_copy_ext(unitcube_forest, 1, 0)
@test p8est_destroy(unitcube_forest_copy) == nothing

# Operations on forests
@test p8est_partition(unitcube_forest, 0, C_NULL) == nothing
@test p8est_partition_ext(unitcube_forest, 0, C_NULL) == 0
@test p8est_balance(unitcube_forest, P4est_wrapper.P8EST_CONNECT_FULL, C_NULL) == nothing
@test p8est_balance_ext(unitcube_forest, P4est_wrapper.P8EST_CONNECT_FULL, C_NULL, C_NULL) == nothing
@test p8est_refine(unitcube_forest, 1, dummy_callback_c, C_NULL) == nothing
@test p8est_refine_ext(unitcube_forest, 1, -1, dummy_callback_c, C_NULL, C_NULL) == nothing
@test p8est_coarsen(unitcube_forest, 1, dummy_callback_c, C_NULL) == nothing
@test p8est_coarsen_ext(unitcube_forest, 1, 0, dummy_callback_c, C_NULL, C_NULL) == nothing
@test p8est_revision(unitcube_forest) == 0
@test typeof(p8est_checksum(unitcube_forest)) == UInt32
@test p8est_connect_type_int(P4est_wrapper.P8EST_CONNECT_FACE) == 1
@test unsafe_string(p8est_connect_type_string(P4est_wrapper.P8EST_CONNECT_FACE)) == "FACE"

# Destroy the forest
@test p8est_destroy(unitcube_forest) == nothing
# Destroy the connectivity
@test p8est_connectivity_destroy(unitcube_connectivity) == nothing

# Finalize MPI if initialized and session is not interactive
if (MPI.Initialized() && !isinteractive())
    MPI.Finalize()
end


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

amat_data = reshape([
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
  0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
  0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
  0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0
], (16,16));

vars = [Symbol("n$i") for i in 1:size(amat_data, 1)]
a = NamedArray(Int.(amat_data), (vars, vars), ("Rows", "Cols"));

a |> display

dag = DAG("ag_example", a)

fn = joinpath(mktempdir(), "ag.dot")
to_graphviz(dag, fn)
#Sys.isapple() && run(`open -a GraphViz.app $(fn)`)

m = [:n3, :n5, :n6, :n15, :n16];
c = [:n4, :n7];

ag = ancestral_graph(dag; m=m, c=c)
indx = filter(x -> !isnothing(x), indexin(vars, names(ag, 1)))
ag = ag[indx, indx]

rg = ribbon_graph(dag; m=m, c=c)
indx = filter(x -> !isnothing(x), indexin(vars, names(rg, 1)))
rg = rg[indx, indx]

maximize(ag) |> display
maximize(rg) |> display


text:
function g_create_or_get(node, name)
    exists(node, name) ? node[name] : g_create(node, name)
end

function a_overwrite(node, name, value)
    # check if the hdf5 attribute already exists.
    # if it does we need to delete it in order to be able
    # to overwrite it.
    if HDF5.h5a_exists(node.id, name)
        a_delete(node, name)
    end
    attrs(node)[name] = value
end

function d_create_and_push(node, name, value)
    if exists(node, name)
        data = node[name]
        # simply expand the dimension by 1
        len = Int(HDF5.get_dims(data)[1][1])
        set_dims!(data, (len+1,))
        # write new entry at the new end of the arrays
        data[len+1] = [ value ]
    else
        # create the required dataset such that it
        # can be expanded indefinitely
        data = d_create(node, name, typeof(value), ((1,),(-1,)), "chunk", (50,))
        # write new entry as the first element of the arrays
        data[1] = [ value ]
    end
end


text:
c = 0
for n = 3:1000
  global c
  v = Rational{BigInt}(1, 2)
  for i = 1:n-1
    v += 2
    v = 1/v
  end
  v += 1
  if length(string(numerator(v))) > length(string(denominator(v)))
    c += 1
  end
end
println(c)

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
using Pkg

# Add julia packages to install like comments below when building the container image
#Pkg.add("ArgParse")
#Pkg.add(PackageSpec(url="https://github.com/staticfloat/InfluxDB.jl", rev="master"))

# Run `julia packages.jl` to install packages above

max_stars_repo_name: hassiweb/docker-julia

max_issues_repo_name: hassiweb/docker-julia

max_forks_repo_name: hassiweb/docker-julia


text:
include("NamedViewVector.jl")

using DifferentialEquations
using Plots

theme(:juno)
plotly()

# -- Constants -----------------------------------------------------------------------------
# Gravitational constant
const G0 = 9.80665

# -- Component equations -------------------------------------------------------------------
function pendulum!(dθ, θ, p, t)
    m, L, ΣM = p.mass, p.L, p.ΣM

    dθ.pos = θ.vel
    dθ.vel = (ΣM - m*G0*sin(θ.pos)) / (m*L)
    return nothing
end

function cart!(dx, x, p, t)
    m, ΣF = p.mass, p.ΣF

    dx.pos = x.vel
    dx.vel = ΣF/m
    return nothing
end

function composed!(du, u, p, t)
    m, L = p.pend.mass, p.pend.L
    M = p.cart.mass

    pendulum!(du.θ, u.θ, (mass=m, L=L, ΣM=), t)
    cart!(du.x, u.x, (mass=m+M, ΣF=), t)
end

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
using Test, MPI

using ClimateMachine
using ClimateMachine.MPIStateArrays
using ClimateMachine.MPIStateArrays: getstateview
using ClimateMachine.VariableTemplates: @vars, varsindex, varsindices
using StaticArrays

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD

const V = @vars begin
    a::Float32
    b::SVector{3, Float32}
    c::SMatrix{3, 8, Float32}
    d::Float32
    e::@vars begin
        a::Float32
        b::SVector{3, Float32}
        d::Float32
    end
end

@testset "MPIStateArray varsindex" begin
    # check with invalid vars size
    @test_throws ErrorException MPIStateArray{Float32, V}(
        mpicomm,
        ArrayType,
        4,
        1,
        8,
    )

    Q = MPIStateArray{Float32, V}(mpicomm, ArrayType, 4, 34, 8)
    @test Q.a === view(MPIStateArrays.realview(Q), :, 1:1, :)
    @test Q.b === view(MPIStateArrays.realview(Q), :, 2:4, :)
    @test Q.c === view(MPIStateArrays.realview(Q), :, 5:28, :)
    @test Q.d === view(MPIStateArrays.realview(Q), :, 29:29, :)
    @test Q.e === view(MPIStateArrays.realview(Q), :, 30:34, :)

    @test getstateview(Q, "a") === view(MPIStateArrays.realview(Q), :, 1:1, :)
    @test getstateview(Q, "b") === view(MPIStateArrays.realview(Q), :, 2:4, :)
    @test getstateview(Q, "c") === view(MPIStateArrays.realview(Q), :, 5:28, :)
    @test getstateview(Q, "d") === view(MPIStateArrays.realview(Q), :, 29:29, :)
    @test getstateview(Q, "e") === view(MPIStateArrays.realview(Q), :, 30:34, :)
    @test getstateview(Q, "e.a") ===
          view(MPIStateArrays.realview(Q), :, 30:30, :)
    @test getstateview(Q, "e.b") ===
          view(MPIStateArrays.realview(Q), :, 31:33, :)
    @test getstateview(Q, "e.d") ===
          view(MPIStateArrays.realview(Q), :, 34:34, :)

    @test getstateview(Q, :(a)) === view(MPIStateArrays.realview(Q), :, 1:1, :)
    @test getstateview(Q, :(b)) === view(MPIStateArrays.realview(Q), :, 2:4, :)
    @test getstateview(Q, :(c)) === view(MPIStateArrays.realview(Q), :, 5:28, :)
    @test getstateview(Q, :(d)) ===
          view(MPIStateArrays.realview(Q), :, 29:29, :)
    @test getstateview(Q, :(e)) ===
          view(MPIStateArrays.realview(Q), :, 30:34, :)
    @test getstateview(Q, :(e.a)) ===
          view(MPIStateArrays.realview(Q), :, 30:30, :)
    @test getstateview(Q, :(e.b)) ===
          view(MPIStateArrays.realview(Q), :, 31:33, :)
    @test getstateview(Q, :(e.d)) ===
          view(MPIStateArrays.realview(Q), :, 34:34, :)

    @test_throws ErrorException Q.aa
    @test_throws ErrorException getstateview(Q, "aa")

    P = similar(Q)
    @test P.a === view(MPIStateArrays.realview(P), :, 1:1, :)
    @test P.b === view(MPIStateArrays.realview(P), :, 2:4, :)
    @test P.c === view(MPIStateArrays.realview(P), :, 5:28, :)
    @test P.d === view(MPIStateArrays.realview(P), :, 29:29, :)
    @test P.e === view(MPIStateArrays.realview(P), :, 30:34, :)

    @test getstateview(P, "a") === view(MPIStateArrays.realview(P), :, 1:1, :)
    @test getstateview(P, "b") === view(MPIStateArrays.realview(P), :, 2:4, :)
    @test getstateview(P, "c") === view(MPIStateArrays.realview(P), :, 5:28, :)
    @test getstateview(P, "d") === view(MPIStateArrays.realview(P), :, 29:29, :)
    @test getstateview(P, "e") === view(MPIStateArrays.realview(P), :, 30:34, :)
    @test getstateview(P, "e.a") ===
          view(MPIStateArrays.realview(P), :, 30:30, :)
    @test getstateview(P, "e.b") ===
          view(MPIStateArrays.realview(P), :, 31:33, :)
    @test getstateview(P, "e.d") ===
          view(MPIStateArrays.realview(P), :, 34:34, :)

    @test getstateview(P, :(a)) === view(MPIStateArrays.realview(P), :, 1:1, :)
    @test getstateview(P, :(b)) === view(MPIStateArrays.realview(P), :, 2:4, :)
    @test getstateview(P, :(c)) === view(MPIStateArrays.realview(P), :, 5:28, :)
    @test getstateview(P, :(d)) ===
          view(MPIStateArrays.realview(P), :, 29:29, :)
    @test getstateview(P, :(e)) ===
          view(MPIStateArrays.realview(P), :, 30:34, :)
    @test getstateview(P, :(e.a)) ===
          view(MPIStateArrays.realview(P), :, 30:30, :)
    @test getstateview(P, :(e.b)) ===
          view(MPIStateArrays.realview(P), :, 31:33, :)
    @test getstateview(P, :(e.d)) ===
          view(MPIStateArrays.realview(P), :, 34:34, :)

    A = MPIStateArray{Float32}(mpicomm, ArrayType, 4, 29, 8)
    @test_throws ErrorException A.a
    @test_throws ErrorException getstateview(A, "a")
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# Check if the record data is in the four column BED format.
function isLikeRecord(line::String) :: Bool
    return  occursin(r"^\s*\S*(?=[A-Za-z0-9])\S*\s+(\d+)\s+(\d+)\s+(\S*\d)\s*$", line) # Note: is like a record.
end

function isBrowser(line::String) :: Bool
    return  occursin(r"^browser", lowercase(line))
end

function isComment(line::String) :: Bool
    return occursin(r"^\s*(?:#|$)", line)
end

function seekNextRecord(io::IO) :: Nothing

    pos = position(io)
    initial = pos == 0 ? -1 : pos # Note: Allows for the fist line of headerless bedGraph file to be read.
    line = ""

    while !eof(io) && (!isLikeRecord(line) || pos == initial)
        pos = position(io)
        line = readline(io)
    end

    seek(io, pos)

    return nothing

end

# Note: all options are placed in a single line separated by spaces.
function readParameters(io::IO) :: String
    seekstart(io)

    pos = position(io)

    while !eof(io) && !isLikeRecord(line) # Note: regex is used to limit the search by exiting the loop when a line matches the bedGraph record format.
        line = readline(io)

        if contains(line, "type=bedGraph") # Note: the track type is REQUIRED, and must be bedGraph.
            return line
        end

    end
end

function readRecord(io::IO) :: Union{Nothing, Record}

    line = readline(io)

    if isLikeRecord(line)
        return Record(line)
    end

    return nothing
end

function readRecords(io::IO) :: Vector{Record}
    seekstart(io)
    seekNextRecord(io)

    records = Vector{Record}()

    while !eof(io)
        record = readRecord(io)
        if record != nothing
            push!(records, record)
        end
    end

    return records

end

function Base.read(io::IO, ::Type{Vector{Record}}) :: Vector{Record}
    return readRecords(io)
end


text:
using Documenter, StructuralEquationModels

makedocs(
    sitename="StructuralEquationModels.jl",
    pages = [
        "index.md",
        "Tutorials" => [
            "A first model" => "tutorials/first_model.md",
            "Our Concept of a Structural Equation Model" => "tutorials/concept.md",
            "Model Specification" => [
                "tutorials/specification/specification.md",
                "tutorials/specification/graph_interface.md",
                "tutorials/specification/ram_matrices.md",
                "tutorials/specification/parameter_table.md"
            ],
            "Model Construction" => [
                "tutorials/construction/construction.md",
                "tutorials/construction/outer_constructor.md",
                "tutorials/construction/build_by_parts.md"
            ],
            "Optimization Backends" => [
                "tutorials/backends/optim.md",
                "tutorials/backends/nlopt.md"
            ],
            "Model Fitting" => "tutorials/fitting/fitting.md",
            "Model Inspection" => "tutorials/inspection/inspection.md",
            "Mean Structures" => "tutorials/meanstructure.md",
            "Collections" => [
                "tutorials/collection/collection.md",
                "tutorials/collection/multigroup.md"
            ],
            "Constraints" => "tutorials/constraints/constraints.md"
        ],
        "Developer documentation" => [
            "Extending the package" => "developer/extending.md",
            "Custom loss functions" => "developer/loss.md",
            "Custom imply types" => "developer/imply.md",
            "Custom optimizer types" => "developer/optimizer.md",
            "Custom observed types" => "developer/observed.md",
            "Custom model types" => "developer/sem.md"
        ],
        "Performance tips" => [
            "Model sorting" => "performance/sorting.md",
            "MKL" => "performance/mkl.md",
            "Simulation studies" => "performance/simulation.md",
            "Symbolic precomputation" => "performance/symbolic.md",
            "Starting values" => "performance/starting_values.md",
            "Parametric Types" => "performance/parametric.md"
        ],
        "Internals and design" => [
            "Internals and design" => "internals/internals.md",
            "files" => "internals/files.md",
            "types" => "internals/types.md"
        ],
        "Complementary material" => [
            "Mathematical appendix" => "complementary/maths.md"
        ]
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/favicon_zeta.ico"],
        ansicolor = true,
        collapselevel = 1
    ),
    doctest = false
)

# doctest(StructuralEquationModels, fix=true)

deploydocs(
    repo = "github.com/StructuralEquationModels/StructuralEquationModels.jl",
    devbranch = "devel"
)


text:
"""
$(SIGNATURES)
"""
function plot_as_quiver(
    mesh::FaceMesh;
    n_steps = 20,
    scale = 30,
    max_length = N_Pixel_Width / n_steps,
    fliprc = false,
    reverser = false,
    reversec = false,
)

    height, width = size(mesh)
    stride = Int64(ceil(max(height, width) / n_steps))

    row_s = Float64[]
    col_s = Float64[]
    v_rows = Float64[]
    v_cols = Float64[]

    # The gradient vectors reflect change of luminosity:outward flow = bigger cells = more luminosity
    # Plots are cosmetically look better when showing the light focusing to a specific location.
    # Those gradients are the opposite direction.
    ϕ = -mesh.corners.ϕ

    v_length = sqrt.(mesh.corners.vr .^ 2 + mesh.corners.vr .^ 2)
    v_max = maximum(v_length)

    mat_vr = (reverser ? -1 : 1) .* mesh.corners.vr * scale / v_max
    mat_vc = (reversec ? -1 : 1) .* mesh.corners.vc * scale / v_max

    mat_vr = clamp.(mat_vr, -max_length, max_length)
    mat_vc = clamp.(mat_vc, -max_length, max_length)

    for row = 1:stride:height, col = 1:stride:width
        reverser ? push!(row_s, row) : push!(row_s, -row)
        reversec ? push!(col_s, -col) : push!(col_s, col)

        push!(v_rows, mat_vr[row, col])
        push!(v_cols, mat_vc[row, col])
    end

    display(
        fliprc ? quiver(row_s, col_s, quiver = (v_rows, v_cols)) :
        quiver(col_s, row_s, quiver = (v_cols, v_rows)),
    )
    return nothing
end

"""
$(SIGNATURES)
"""
function plot_as_quiver(
    ϕ;
    n_steps = 20,
    scale = 300,
    max_length = N_Pixel_Width / n_steps,
    fliprc = false,
    reverser = false,
    reversec = false,
)

    height, width = size(ϕ)
    stride = Int64(ceil(max(height, width) / n_steps))

    row_s = Float64[]
    col_s = Float64[]
    Δrow = Float64[]
    Δcol = Float64[]

    for row = 1:stride:height, col = 1:stride:width
        dr = (-(ϕ[row, col+1] - ϕ[row, col]) * scale)
        dc = (ϕ[row+1, col] - ϕ[row, col]) * scale

        reverser ? push!(row_s, col) : push!(row_s, width + 1 - col)
        reversec ? push!(col_s, height + 1 - row) : push!(col_s, row)
        reverser && (dr = -dr)
        reversec && (dc = -dc)

        push!(Δrow, clamp(dr, -maxlength, max_length))
        push!(Δcol, clamp(dc, -maxlength, max_length))
    end

    # By default, we flip to match the original image.
    display(
        fliprc ? quiver(row_s, col_s, quiver = (Δrow, Δcol)) :
        quiver(col_s, row_s, quiver = (Δcol, Δrow)),
    )

    return nothing
end

"""
$(SIGNATURES)
"""
function plot_velocities_as_quiver(
    vr,
    vc;
    n_steps = 20,
    scale = 300,
    max_length = N_Pixel_Width / n_steps,
)

    height, width = size(vr)
    stride = Int64(ceil(max(height, width) / n_steps))

    row_s = Float64[]
    col_s = Float64[]
    Δrow = Float64[]
    Δcol = Float64[]

    dr = max.(vr, 0.001)
    dc = max.(vc, 0.001)

    for row = 1:stride:height, col = 1:stride:width
        push!(row_s, height - row)
        push!(col_s, col)

        push!(Δrow, dr[row, col])
        push!(Δcol, dc[row, col])
    end

    # readline()
    display(quiver(row_s, col_s, quiver = (Δrow, Δcol), aspect_ratio = :equal))
end

"""
$(SIGNATURES)
"""
function save_plot_scalar_field!(scalar_field, filename)
    blue = zeros(Float64, size(scalar_field))
    red = zeros(Float64, size(scalar_field))
    green = zeros(Float64, size(scalar_field))

    height, width = size(scalar_field)

    for r ∈ 1:height, c ∈ 1:width
        field = scalar_field[r, c]
        if field < -1.0
            # Green if luminosity error < -1.0
            red[r, c] = 0.0
            blue[r, c] = 0.0
            green[r, c] = 1.0

        elseif -1.0 <= field < 0.0
            # Red shade if negative
            red[r, c] = -field
            blue[r, c] = 0.0
            green[r, c] = 0.0

        elseif 0.0 <= field < 1.0
            # Blue shade if positive
            red[r, c] = 0.0
            blue[r, c] = field
            green[r, c] = 0.0

        elseif 1.0 <= field
            # White  if > 1
            red[r, c] = 1.0
            blue[r, c] = 1.0
            green[r, c] = 1.0
        end
    end

    rgbImg = RGB.(red, green, blue)
    save("./examples/$(filename).png", rgbImg)
end

"""
$(SIGNATURES)
"""
function save_plot_scalar_field(scalar_field, prefix, img)
    normalised_D_max = scalar_field
    normalised_D_min = scalar_field
    # normalised_D_max = scalar_field ./ maximum(scalar_field)
    # normalised_D_min = scalar_field ./ minimum(scalar_field)

    blue = zeros(size(scalar_field))
    red = zeros(size(scalar_field))
    green = zeros(size(scalar_field))

    blue[scalar_field.>0] = normalised_D_max[scalar_field.>0]
    red[scalar_field.<0] = -normalised_D_min[scalar_field.<0]

    rgbImg = map(clamp01nan, RGB.(red, green, blue))
    save("./examples/$(prefix)_loss.png", map(clamp01nan, rgbImg))

    println("Saving output image:")
    println(typeof(img))
    E = Gray.(scalar_field)
    println(typeof(E))
    outputImg = img - E
    save("./examples/$(prefix)_actual.png", outputImg)
end


text:
# ---
# title: 629. K Inverse Pairs Array
# id: problem629
# author: Indigo
# date: 2021-06-22
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/k-inverse-pairs-array/description/>
# hidden: true
# ---
#
# Given two integers `n` and `k`, find how many different arrays consist of
# numbers from `1` to `n` such that there are exactly `k` inverse pairs.
#
# We define an inverse pair as following: For `ith` and `jth` element in the
# array, if `i` < `j` and `a[i]` > `a[j]` then it's an inverse pair; Otherwise,
# it's not.
#
# Since the answer may be very large, the answer should be modulo 109 \+ 7.
#
# **Example 1:**
#
#
#
#     Input: n = 3, k = 0
#     Output: 1
#     Explanation:
#     Only the array [1,2,3] which consists of numbers from 1 to 3 has exactly 0 inverse pair.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: n = 3, k = 1
#     Output: 2
#     Explanation:
#     The array [1,3,2] and [2,1,3] have exactly 1 inverse pair.
#
#
#
#
# **Note:**
#
#   1. The integer `n` is in the range [1, 1000] and `k` is in the range [0, 1000].
#
#
#
#
## @lc code=start
using LeetCode

function k_inverse_pairs(n::Int, k::Int)
    dp = OffsetArray(fill(0, n + 1, k + 1), -1, -1)
    dp[:, 0] .= 1
    m = 1000000007
    @inbounds for i in 1:n, j in 1:min(k, i * (i - 1) ÷ 2)
        val = dp[i - 1, j] - ((j - i) >= 0 ? dp[i - 1, j - i] : 0)
        dp[i, j] = (dp[i, j - 1] + val) % m
    end
    return (dp[n, k] + m) % m
end
## @lc code=end


text:
module FieldMetadata

export @metadata, @chain

"""
Generate a macro that constructs methods of the same name.
These methods return the metadata information provided for each
field of the struct.

```julia
@metadata def_range (0, 0)
@def_range struct Model
    a::Int | (1, 4)
    b::Int | (4, 9)
end

model = Model(3, 5)
def_range(model, Val{:a})
(1, 4)

def_range(model)
((1, 4), (4, 9))
```
"""
macro metadata(name, default)
    symname = QuoteNode(name)
    default = esc(default)
    rename = esc(Meta.parse("re$name"))
    name = esc(name)
    return quote
        macro $name(ex)
            name = $symname
            return add_field_funcs(ex, name)
        end

        macro $rename(ex)
            name = $symname
            return add_field_funcs(ex, name; update=true)
        end

        # Single field methods
        $name(x, key) = $default
        $name(x::Type, key::Type) = $default
        $name(::X, key::Symbol) where X = $name(X, Val{key})
        $name(::X, key::Type) where X = $name(X, key)
        $name(::Type{X}, key::Symbol) where X = $name(X, Val{key})

        # All field methods
        $name(::X) where X = $name(X)
        $name(x::Type{X}) where X = $name(X, fieldname_vals(X))
        $name(::Type{X}, keys::Tuple) where X =
            ($name(X, keys[1]), $name(X, Base.tail(keys))...)
        $name(::Type{X}, keys::Tuple{}) where X = tuple()
    end
end

"""
Chain together any macros. Useful for combining @metadata macros.

### Example
```julia
@chain columns @label @units @default_kw

@columns struct Foo
  bar::Int | 7 | u"g" | "grams of bar"
end
```
"""
macro chain(name, ex)
    macros = chained_macros(ex)
    return quote
        macro $(esc(name))(ex)
            macros = $macros
            for mac in reverse(macros)
                ex = Expr(:macrocall, mac, LineNumberNode(75, "FieldMetadata.jl"), ex)
            end
            esc(ex)
        end
    end
end

Base.@pure fieldname_vals(::Type{X}) where X = ([Val{fn} for fn in fieldnames(X)]...,)

function add_field_funcs(ex, name; update=false)
    macros = chained_macros(ex)

    typ = firsthead(x -> namify(x.args[2]), ex, :struct)
    func_exps = Expr[]

    # Parse the block of lines inside the struct.
    # Function expressions are built for each field, and metadata removed.
    firsthead(ex, :block) do block
        for (i, line) in enumerate(block.args)
            :head in fieldnames(typeof(line)) || continue
            if line.head == :(=)
                call = line.args[2]
                call.args[1] == :(|) || continue
                # The fieldname is the first arg
                fn = line.args[1]
                addmethod_unlessdefault!(func_exps, name, typ, fn, call)
                # Replace the rest of the line after the = call
                line.args[2] = call.args[2]
            elseif line.head == :call || line.args[1] == :(|)
                # The fieldname is the second arg
                fn = line.args[2]
                if :head in fieldnames(typeof(fn))
                    addmethod_unlessdefault!(func_exps, name, typ, fn, line)
                    # Replace the line with the field
                    line.head = fn.head
                    line.args = fn.args
                else
                    addmethod_unlessdefault!(func_exps, name, typ, fn, line)
                    # Replace the line in the parent block
                    block.args[i] = line.args[2]
                end
            end
        end
    end
    if update && length(macros) == 0
        Expr(:block, func_exps...)
    else
        Expr(:block, :(Base.@__doc__ $(esc(ex))), func_exps...)
    end
end

function addmethod_unlessdefault!(func_exps, name, typ, fn, ex)
    # Get just the fieldname symbol from the full fieldname
    val = ex.args[3]
    key = getkey(fn)
    # Add a method expression unless this contains a _ default
    val == :_ || addmethod!(func_exps, name, typ, key, val)
end

function addmethod!(func_exps, name, typ, key, value)
    func = esc(:(function $name(::Type{<:$typ}, ::Type{Val{$(QuoteNode(key))}}) $value end))
    push!(func_exps, func)
end

getkey(ex::Expr) = begin
    key = firsthead(y -> y.args[1], ex, :(::))
    if key == nothing
        key = firsthead(y -> y.args[1], ex, :(=))
    end
    if key == nothing
        key = firsthead(y -> y.args[2], ex, :call)
    end
    key
end
getkey(ex::Symbol) = ex

chained_macros(ex) = chained_macros!(Symbol[], ex)

chained_macros!(macros, ex::Expr) = begin
    if ex.head == :macrocall
        push!(macros, ex.args[1])
        length(ex.args) > 2 && chained_macros!(macros, ex.args[3])
    end
    macros
end

firsthead(f, ex::Expr, sym) =
    if ex.head == sym
        out = f(ex)
        return out
    else
        for arg in ex.args
            x = firsthead(f, arg, sym)
            x == nothing || return x
        end
        return nothing
    end
firsthead(f, ex, sym) = nothing

namify(x) = x
namify(x::Expr) = namify(x.args[1])

# FieldMetadata api
@metadata default nothing
@metadata units 1
@metadata prior nothing
@metadata description ""
@metadata limits (1e-7, 1.0) # just above zero so log transform is possible
@metadata bounds (1e-7, 1.0) # just above zero so log transform is possible
@metadata label ""
@metadata logscaled false
@metadata flattenable true
@metadata plottable true
@metadata selectable Nothing

# Set the default label to be the field name
label(x::Type, ::Type{Val{F}}) where F = F

end # module


text:
mutable struct AbstractTree
    node_counter::Int
    nodes::Vector{AbstractNode}
    processed::Vector{AbstractNode}
    best_bound::Real     # e.g., lower bound for minimization
    best_incumbent::Real # e.g., upper bound for minimization

    function AbstractTree(node_counter::Int = 0)
        return new(node_counter, [], [], -Inf, Inf)
    end
end

function run(tree::AbstractTree)
    while !termination(tree)
        node = next_node(tree)
        update_best_bound!(tree)

        bound!(node)
        heuristics!(node)
        processed!(tree, node)
        branch!(tree, node)
    end
end

# basic termination function
termination(tree::AbstractTree) = isempty(tree)
isempty(tree::AbstractTree) = Base.isempty(tree.nodes)

# add node to tree
function push!(tree::AbstractTree, node::AbstractNode)
    tree.node_counter += 1
    node.id = tree.node_counter
    Base.push!(tree.nodes, node)
    sort!(tree)
end

# add a set of nodes to tree
function push!(tree::AbstractTree, nodes::Vector{T}) where T<:AbstractNode
    for node in nodes
        push!(tree, node)
    end
end

# mark node as processed
processed!(tree::AbstractTree, node::AbstractNode) = Base.push!(tree.processed, node)

function update_best_bound!(tree::AbstractTree)
    if !isempty(tree)
        tree.best_bound = Base.minimum([node.bound for node in tree.nodes])
    end
end

function branch!(tree::AbstractTree, node::AbstractNode)
    children = branch(node)
    push!(tree, children)
end

# return the next search node
function next_node(tree::AbstractTree)
    sort!(tree)
    node = Base.pop!(tree.nodes)
    apply_changes!(node)
    return node
end

# best bound
sort!(tree::AbstractTree) = Base.sort!(tree.nodes, by=x->x.bound, rev=true)