text:
using Test
using TimeseriesSurrogates
ENV["GKSwstype"] = "100"

N = 500
ts = cumsum(randn(N))
ts_nan = cumsum(randn(N))
ts_nan[1] = NaN
x = cos.(range(0, 20π, length = N)) .+ randn(N)*0.05

@testset "LombScargle" begin
    t = sort((0:N-1) + rand(N))
    tol = 10
    ls = IrregularLombScargle(t, tol = 10, n_total = 20000, n_acc = 5000)

    s = surrogate(x, ls)

    @test all(sort(s) .== sort(x))
end

@testset "WLS" begin
    wts_norescale = WLS(AAFT(), rescale = false)
    wts_rescale = WLS(AAFT())

    s_norescale = surrogate(x, wts_norescale)
    s_rescale = surrogate(x, wts_rescale)
    @test length(s_norescale) == length(x)
    @test length(s_rescale) == length(x)

    # If rescaling, the surrogate will have the same values as the original
    @test sort(x) ≈ sort(s_rescale)

end

@testset "PartialRandomization" begin
    pr = PartialRandomization(0.2)
    s = surrogate(x, pr)

    @test length(s) == length(x)

    @test_throws AssertionError PartialRandomization(-0.01)
    @test_throws AssertionError PartialRandomization(1.01)
end

@testset "PartialRandomizationAAFT" begin
    praaft = PartialRandomizationAAFT(0.5)
    s = surrogate(x, praaft)

    @test length(s) == length(x)
    @test sort(x) ≈ sort(s)

    @test_throws AssertionError PartialRandomizationAAFT(-0.01)
    @test_throws AssertionError PartialRandomizationAAFT(1.01)
end

@testset "RandomCascade" begin
    randomcascade = RandomCascade()
    s = surrogate(x, randomcascade)

    @test length(s) == length(x)
end

@testset "Periodic" begin
	pp = PseudoPeriodic(3, 25, 0.05)
	s = surrogate(x, pp)
	@test length(s) == length(ts)
	@test all(s[i] ∈ x for i in 1:N)
	# Perhaps a more advanced test, e.g. that both components have Fourier peak at
	# the same frequency, should be considered.

	#TODO: Test for noiseradius
end

@testset "PeriodicTwin" begin
    # A better test based on recurrence plots should be considered.
    d, τ = 2, 6
    δ = 0.2
    ρ = noiseradius(x, d, τ, 0.02:0.02:0.5)
    method = PseudoPeriodicTwin(d, τ, δ, ρ)

    sg = surrogenerator(x, method)
    s = sg()[:, 1]
    @test length(s) == length(ts)
    @test all(s[i] ∈ x for i in 1:N)
end

@testset "BlockShuffle" begin
    bs1 = BlockShuffle()
    bs2 = BlockShuffle(4)
    s1 = surrogate(x, bs1)
    s2 = surrogate(x, bs2)

    @test length(s1) == length(x)
    @test length(s2) == length(x)
    @test all([s1[i] ∈ x for i = 1:N])
    @test all([s2[i] ∈ x for i = 1:N])
end

@testset "RandomShuffle" begin
    rs = RandomShuffle()
    s = surrogate(x, rs)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "AutoRegressive" begin
	y = TimeseriesSurrogates.AR1(2000, 0.1, 0.5)
    sg = surrogenerator(y, AutoRegressive(1))
	@test 0.4 ≤ abs(sg.init.φ[1]) ≤ 0.6
	s = sg()
    @test length(s) == length(y)
end

@testset "AAFT" begin
    aaft = AAFT()
    s = surrogate(x, aaft)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "IAAFT" begin
    iaaft = IAAFT()
    s = surrogate(x, iaaft)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "TFTS" begin
    method_preserve_lofreq = TFTS(0.05)
    method_preserve_hifreq = TFTS(-0.05)

    s = surrogate(x, method_preserve_lofreq)
    @test length(s) == length(x)

    s = surrogate(x, method_preserve_hifreq)
    @test length(s) == length(x)

    @test_throws ArgumentError TFTS(0)
end

@testset "TFTDAAFT" begin
    tftdaaft = TFTDAAFT(0.03)
    s = surrogate(x, tftdaaft)
    @test length(s) == length(x)
    @test sort(x) ≈ sort(s)
    @test_throws ArgumentError TFTD(0)
    @test_throws ArgumentError TFTD(1.2)
end

@testset "TFTDIAAFT" begin
    tftdiaaft = TFTDAAFT(0.05)
    s = surrogate(x, tftdiaaft)
    @test length(s) == length(x)
    @test sort(x) ≈ sort(s)
    @test_throws ArgumentError TFTD(0)
    @test_throws ArgumentError TFTD(1.2)
end

@testset "TAAFT" begin
    method_preserve_lofreq = TAAFT(0.05)
    method_preserve_hifreq = TAAFT(-0.05)

    s = surrogate(x, method_preserve_lofreq)
    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])

    s = surrogate(x, method_preserve_hifreq)
    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])

    @test_throws ArgumentError TAAFT(0)
end

@testset "RandomFourier" begin
    @testset "random phases" begin
        phases = true
        rf = RandomFourier(phases)
        s = surrogate(x, rf)

        @test length(s) == length(x)
    end

    @testset "random amplitudes" begin
        phases = false
        rf = RandomFourier(phases)
        s = surrogate(x, rf)

        @test length(s) == length(x)
    end
end

@testset "TFTDRandomFourier" begin
    @testset "random phases" begin
        phases = true
        rf = TFTDRandomFourier(phases)
        s = surrogate(x, rf)

        @test length(s) == length(x)
    end

    @testset "random amplitudes" begin
        phases = false
        rf = TFTDRandomFourier(phases)
        s = surrogate(x, rf)

        @test length(s) == length(x)
    end
end

@testset "Circ/Cycle shuffle" begin
	x = random_cycles()
	s = surrogate(x, CycleShuffle())
	for a in s
		@test a ∈ x
	end
	s = surrogate(x, CircShift(1:length(x)))
	for a in s
		@test a ∈ x
	end
end

using DelayEmbeddings
@testset "ShufleDims" begin
	X = Dataset(rand(100, 3))
	Y = surrogate(X, ShuffleDimensions())
	for i in 1:100
		@test sort(X[i]) == sort(Y[i])
	end
end


text:
@testset "Testing the reparametrisation functions" begin
    δ, γ, σ = 1.1, 1.6, 2.3
    for pidx in eachindex([δ, γ, σ])
        @test [δ, γ, σ][pidx] ≈ DualOptimalFiltering.inverse_reparam_CIR(DualOptimalFiltering.reparam_CIR(δ, γ, σ)...)[pidx]
    end
    for pidx in eachindex([δ, γ, σ])
        @test [δ, γ, σ][pidx] ≈ DualOptimalFiltering.reparam_CIR(DualOptimalFiltering.inverse_reparam_CIR(δ, γ, σ)...)[pidx]
    end
end

ext: jl

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

#####################
# OptimizationState #
#####################

mutable struct OptimizationState
    linfo::MethodInstance
    result_vargs::Vector{Any}
    calledges::Vector{Any}
    src::CodeInfo
    mod::Module
    nargs::Int
    min_valid::UInt
    max_valid::UInt
    params::Params
    sp::SimpleVector # static parameters
    const_api::Bool
    function OptimizationState(frame::InferenceState)
        s_edges = frame.stmt_edges[1]
        if s_edges === ()
            s_edges = []
            frame.stmt_edges[1] = s_edges
        end
        src = frame.src
        return new(frame.linfo, frame.result.vargs,
                   s_edges::Vector{Any},
                   src, frame.mod, frame.nargs,
                   frame.min_valid, frame.max_valid,
                   frame.params, frame.sp, false)
    end
    function OptimizationState(linfo::MethodInstance, src::CodeInfo,
                               params::Params)
        # prepare src for running optimization passes
        # if it isn't already
        nssavalues = src.ssavaluetypes
        if nssavalues isa Int
            src.ssavaluetypes = Any[ Any for i = 1:nssavalues ]
        end
        if src.slottypes === nothing
            nslots = length(src.slotnames)
            src.slottypes = Any[ Any for i = 1:nslots ]
        end
        s_edges = []
        # cache some useful state computations
        toplevel = !isa(linfo.def, Method)
        if !toplevel
            meth = linfo.def
            inmodule = meth.module
            nargs = meth.nargs
        else
            inmodule = linfo.def::Module
            nargs = 0
        end
        result_vargs = Any[] # if you want something more accurate, set it yourself :P
        return new(linfo, result_vargs,
                   s_edges::Vector{Any},
                   src, inmodule, nargs,
                   min_world(linfo), max_world(linfo),
                   params, spvals_from_meth_instance(linfo), false)
        end
end

function OptimizationState(linfo::MethodInstance, params::Params)
    src = retrieve_code_info(linfo)
    src === nothing && return nothing
    return OptimizationState(linfo, src, params)
end

#############
# constants #
#############

# The slot has uses that are not statically dominated by any assignment
# This is implied by `SLOT_USEDUNDEF`.
# If this is not set, all the uses are (statically) dominated by the defs.
# In particular, if a slot has `AssignedOnce && !StaticUndef`, it is an SSA.
const SLOT_STATICUNDEF  = 1
const SLOT_ASSIGNEDONCE = 16 # slot is assigned to only once
const SLOT_USEDUNDEF    = 32 # slot has uses that might raise UndefVarError
# const SLOT_CALLED      = 64

const IR_FLAG_INBOUNDS = 0x01

# known affect-free calls (also effect-free)
const _PURE_BUILTINS = Any[tuple, svec, fieldtype, apply_type, ===, isa, typeof, UnionAll, nfields]

# known effect-free calls (might not be affect-free)
const _PURE_BUILTINS_VOLATILE = Any[getfield, arrayref, isdefined, Core.sizeof]

const TOP_TUPLE = GlobalRef(Core, :tuple)

#########
# logic #
#########

_topmod(sv::OptimizationState) = _topmod(sv.mod)

function update_valid_age!(min_valid::UInt, max_valid::UInt, sv::OptimizationState)
    sv.min_valid = max(sv.min_valid, min_valid)
    sv.max_valid = min(sv.max_valid, max_valid)
    @assert(!isa(sv.linfo.def, Method) ||
            (sv.min_valid == typemax(UInt) && sv.max_valid == typemin(UInt)) ||
            sv.min_valid <= sv.params.world <= sv.max_valid,
            "invalid age range update")
    nothing
end

update_valid_age!(li::MethodInstance, sv::OptimizationState) = update_valid_age!(min_world(li), max_world(li), sv)

function add_backedge!(li::MethodInstance, caller::OptimizationState)
    isa(caller.linfo.def, Method) || return # don't add backedges to toplevel exprs
    push!(caller.calledges, li)
    update_valid_age!(li, caller)
    nothing
end

function isinlineable(m::Method, me::OptimizationState, bonus::Int=0)
    # compute the cost (size) of inlining this code
    inlineable = false
    cost_threshold = me.params.inline_cost_threshold
    if m.module === _topmod(m.module)
        # a few functions get special treatment
        name = m.name
        sig = m.sig
        if ((name === :+ || name === :* || name === :min || name === :max) &&
            isa(sig,DataType) &&
            sig == Tuple{sig.parameters[1],Any,Any,Any,Vararg{Any}})
            inlineable = true
        elseif (name === :iterate || name === :unsafe_convert ||
                name === :cconvert)
            cost_threshold *= 4
        end
    end
    if !inlineable
        inlineable = inline_worthy(me.src.code, me.src, me.sp, me.params, cost_threshold + bonus)
    end
    return inlineable
end

# run the optimization work
function optimize(opt::OptimizationState, @nospecialize(result))
    def = opt.linfo.def
    nargs = Int(opt.nargs) - 1
    @timeit "optimizer" ir = run_passes(opt.src, nargs, opt)
    force_noinline = any(x -> isexpr(x, :meta) && x.args[1] == :noinline, ir.meta)
    replace_code_newstyle!(opt.src, ir, nargs)

    # compute inlining and other related optimizations
    if (isa(result, Const) || isconstType(result))
        proven_pure = false
        # must be proven pure to use const_api; otherwise we might skip throwing errors
        # (issue #20704)
        # TODO: Improve this analysis; if a function is marked @pure we should really
        # only care about certain errors (e.g. method errors and type errors).
        if length(opt.src.code) < 10
            proven_pure = true
            for i in 1:length(opt.src.code)
                if !statement_effect_free(opt.src.code[i], opt, opt.src.ssavaluetypes[i])
                    proven_pure = false
                    break
                end
            end
            if proven_pure
                for fl in opt.src.slotflags
                    if (fl & SLOT_USEDUNDEF) != 0
                        proven_pure = false
                        break
                    end
                end
            end
        end
        if proven_pure
            opt.src.pure = true
        end

        if proven_pure && !coverage_enabled()
            # use constant calling convention
            # Do not emit `jl_fptr_const_return` if coverage is enabled
            # so that we don't need to add coverage support
            # to the `jl_call_method_internal` fast path
            # Still set pure flag to make sure `inference` tests pass
            # and to possibly enable more optimization in the future
            if !(isa(result, Const) && !is_inlineable_constant(result.val))
                opt.const_api = true
            end
            force_noinline || (opt.src.inlineable = true)
        end
    end

    # determine and cache inlineability
    if !force_noinline
        # don't keep ASTs for functions specialized on a Union argument
        # TODO: this helps avoid a type-system bug mis-computing sparams during intersection
        sig = unwrap_unionall(opt.linfo.specTypes)
        if isa(sig, DataType) && sig.name === Tuple.name
            for P in sig.parameters
                P = unwrap_unionall(P)
                if isa(P, Union)
                    force_noinline = true
                    break
                end
            end
        else
            force_noinline = true
        end
    end
    if force_noinline
        opt.src.inlineable = false
    elseif !opt.src.inlineable && isa(def, Method)
        bonus = 0
        if result ⊑ Tuple && !isbitstype(widenconst(result))
            bonus = opt.params.inline_tupleret_bonus
        end
        opt.src.inlineable = isinlineable(def, opt, bonus)
    end
    nothing
end

# whether `f` is pure for inference
function is_pure_intrinsic_infer(f::IntrinsicFunction)
    return !(f === Intrinsics.pointerref || # this one is volatile
             f === Intrinsics.pointerset || # this one is never effect-free
             f === Intrinsics.llvmcall ||   # this one is never effect-free
             f === Intrinsics.arraylen ||   # this one is volatile
             f === Intrinsics.sqrt_llvm ||  # this one may differ at runtime (by a few ulps)
             f === Intrinsics.cglobal)  # cglobal lookup answer changes at runtime
end

# whether `f` is pure for optimizations
function is_pure_intrinsic_optim(f::IntrinsicFunction)
    return !(f === Intrinsics.pointerref || # this one is volatile
             f === Intrinsics.pointerset || # this one is never effect-free
             f === Intrinsics.llvmcall ||   # this one is never effect-free
             f === Intrinsics.arraylen ||   # this one is volatile
             f === Intrinsics.checked_sdiv_int ||  # these may throw errors
             f === Intrinsics.checked_udiv_int ||
             f === Intrinsics.checked_srem_int ||
             f === Intrinsics.checked_urem_int ||
             f === Intrinsics.cglobal)  # cglobal throws an error for symbol-not-found
end

function is_pure_builtin(@nospecialize(f))
    if isa(f, IntrinsicFunction)
        return is_pure_intrinsic_optim(f)
    elseif isa(f, Builtin)
        return (contains_is(_PURE_BUILTINS, f) ||
                contains_is(_PURE_BUILTINS_VOLATILE, f))
    else
        return f === return_type
    end
end

function statement_effect_free(@nospecialize(e), me::OptimizationState, @nospecialize(etype))
    if isa(e, Expr)
        if e.head === :(=)
            return !isa(e.args[1], GlobalRef) && effect_free(e.args[2], me, false, etype)
        elseif e.head === :gotoifnot
            return effect_free(e.args[1], me, false, etype)
        end
    elseif isa(e, GotoNode)
        return true
    end
    return effect_free(e, me, false, etype)
end

effect_free(@nospecialize(e), s::InferenceState, allow_volatile::Bool, @nospecialize(etype)) =
    effect_free(e, s.src, s.sp, allow_volatile, etype)

effect_free(@nospecialize(e), s::OptimizationState, allow_volatile::Bool, @nospecialize(etype)) =
    effect_free(e, s.src, s.sp, allow_volatile, etype)

# detect some important side-effect-free calls (allow_volatile=true)
# and some affect-free calls (allow_volatile=false) -- affect_free means the call
# cannot be affected by previous calls, except assignment nodes
function effect_free(@nospecialize(e), src, spvals::SimpleVector, allow_volatile::Bool, @nospecialize(etype))
    if isa(e, GlobalRef)
        return (isdefined(e.mod, e.name) && (allow_volatile || isconst(e.mod, e.name)))
    elseif isa(e, Slot)
        return src.slotflags[slot_id(e)] & SLOT_USEDUNDEF == 0
    elseif isa(e, Expr)
        e = e::Expr
        head = e.head
        if is_meta_expr_head(head)
            return true
        end
        if head === :static_parameter
            # if we aren't certain enough about the type, it might be an UndefVarError at runtime
            return isa(etype, Const) || issingletontype(widenconst(etype))
        end
        if etype === Bottom
            return false
        end
        ea = e.args
        if head === :call
            if is_known_call_p(e, is_pure_builtin, src, spvals)
                if !allow_volatile
                    if is_known_call(e, arrayref, src, spvals) || is_known_call(e, arraylen, src, spvals)
                        return false
                    elseif is_known_call(e, getfield, src, spvals)
                        nargs = length(ea)
                        (3 <= nargs <= 4) || return false
                        # TODO: check ninitialized
                        if !isa(etype, Const) && !isconstType(etype)
                            # first argument must be immutable to ensure e is affect_free
                            a = ea[2]
                            typ = unwrap_unionall(widenconst(argextype(a, src, spvals)))
                            if isType(typ)
                                # all fields of subtypes of Type are effect-free
                                # (including the non-inferrable uid field)
                            elseif !isa(typ, DataType) || typ.abstract || (typ.mutable && length(typ.types) > 0)
                                return false
                            end
                        end
                    end
                end
                # fall-through
            elseif is_known_call(e, _apply, src, spvals) && length(ea) > 1
                ft = argextype(ea[2], src, spvals)
                if !isa(ft, Const) || (!contains_is(_PURE_BUILTINS, ft.val) &&
                                       ft.val !== Core.sizeof)
                    return false
                end
                # fall-through
            else
                return false
            end
        elseif head === :new
            a = ea[1]
            typ = argextype(a, src, spvals)
            # `Expr(:new)` of unknown type could raise arbitrary TypeError.
            typ, isexact = instanceof_tfunc(typ)
            isexact || return false
            isconcretedispatch(typ) || return false
            typ = typ::DataType
            if !allow_volatile && typ.mutable
                return false
            end
            fieldcount(typ) >= length(ea) - 1 || return false
            for fld_idx in 1:(length(ea) - 1)
                eT = argextype(ea[fld_idx + 1], src, spvals)
                fT = fieldtype(typ, fld_idx)
                eT ⊑ fT || return false
            end
            # fall-through
        elseif head === :return
            # fall-through
        elseif head === :isdefined
            return allow_volatile
        elseif head === :the_exception
            return allow_volatile
        elseif head === :copyast
            return true
        else
            return false
        end
        for a in ea
            if !effect_free(a, src, spvals, allow_volatile, argextype(a, src, spvals))
                return false
            end
        end
    elseif isa(e, GotoNode)
        return false
    end
    return true
end

## Computing the cost of a function body

# saturating sum (inputs are nonnegative), prevents overflow with typemax(Int) below
plus_saturate(x, y) = max(x, y, x+y)

# known return type
isknowntype(@nospecialize T) = (T == Union{}) || isconcretetype(T)

function statement_cost(ex::Expr, line::Int, src::CodeInfo, spvals::SimpleVector, params::Params)
    head = ex.head
    if is_meta_expr_head(head) || head == :copyast # not sure if copyast is right
        return 0
    end
    argcost = 0
    for a in ex.args
        if a isa Expr
            argcost = plus_saturate(argcost, statement_cost(a, -1, src, spvals, params))
        end
    end
    if head == :return || head == :(=)
        return argcost
    end
    extyp = line == -1 ? Any : src.ssavaluetypes[line]
    if head == :call
        ftyp = argextype(ex.args[1], src, spvals)
        if isa(ftyp, Type)
            return argcost
        end
        if isa(ftyp, Const)
            f = (ftyp::Const).val
            if isa(f, IntrinsicFunction)
                iidx = Int(reinterpret(Int32, f::IntrinsicFunction)) + 1
                if !isassigned(T_IFUNC_COST, iidx)
                    # unknown/unhandled intrinsic
                    return plus_saturate(argcost, params.inline_nonleaf_penalty)
                end
                return plus_saturate(argcost, T_IFUNC_COST[iidx])
            end
            if isa(f, Builtin)
                # The efficiency of operations like a[i] and s.b
                # depend strongly on whether the result can be
                # inferred, so check the type of ex
                if f == Main.Core.getfield || f == Main.Core.tuple
                    # we might like to penalize non-inferrability, but
                    # tuple iteration/destructuring makes that
                    # impossible
                    # return plus_saturate(argcost, isknowntype(extyp) ? 1 : params.inline_nonleaf_penalty)
                    return argcost
                elseif f == Main.Core.arrayref && length(ex.args) >= 3
                    atyp = argextype(ex.args[3], src, spvals)
                    return plus_saturate(argcost, isknowntype(atyp) ? 4 : params.inline_nonleaf_penalty)
                end
                fidx = findfirst(x->x===f, T_FFUNC_KEY)
                if fidx === nothing
                    # unknown/unhandled builtin or anonymous function
                    # Use the generic cost of a direct function call
                    return plus_saturate(argcost, 20)
                end
                return plus_saturate(argcost, T_FFUNC_COST[fidx])
            end
        end
        return plus_saturate(argcost, params.inline_nonleaf_penalty)
    elseif head == :foreigncall || head == :invoke
        # Calls whose "return type" is Union{} do not actually return:
        # they are errors. Since these are not part of the typical
        # run-time of the function, we omit them from
        # consideration. This way, non-inlined error branches do not
        # prevent inlining.
        return extyp == Union{} ? 0 : plus_saturate(20, argcost)
    elseif head == :llvmcall
        return plus_saturate(10, argcost) # a wild guess at typical cost
    elseif head == :enter
        # try/catch is a couple function calls,
        # but don't inline functions with try/catch
        # since these aren't usually performance-sensitive functions,
        # and llvm is more likely to miscompile them when these functions get large
        return typemax(Int)
    elseif head == :gotoifnot
        target = ex.args[2]::Int
        # loops are generally always expensive
        # but assume that forward jumps are already counted for from
        # summing the cost of the not-taken branch
        return target < line ? plus_saturate(40, argcost) : argcost
    end
    return argcost
end

function inline_worthy(body::Array{Any,1}, src::CodeInfo, spvals::SimpleVector, params::Params,
                       cost_threshold::Integer=params.inline_cost_threshold)
    bodycost = 0
    for line = 1:length(body)
        stmt = body[line]
        if stmt isa Expr
            thiscost = statement_cost(stmt, line, src, spvals, params)::Int
        elseif stmt isa GotoNode
            # loops are generally always expensive
            # but assume that forward jumps are already counted for from
            # summing the cost of the not-taken branch
            thiscost = stmt.label < line ? 40 : 0
        else
            continue
        end
        bodycost = plus_saturate(bodycost, thiscost)
        bodycost == typemax(Int) && return false
    end
    return bodycost <= cost_threshold
end

function is_known_call(e::Expr, @nospecialize(func), src, spvals)
    if e.head !== :call
        return false
    end
    f = argextype(e.args[1], src, spvals)
    return isa(f, Const) && f.val === func
end

function is_known_call_p(e::Expr, @nospecialize(pred), src, spvals)
    if e.head !== :call
        return false
    end
    f = argextype(e.args[1], src, spvals)
    return (isa(f, Const) && pred(f.val)) || (isType(f) && pred(f.parameters[1]))
end

function renumber_stuff!(body::Vector{Any}, changemap::Vector{Int})
    for i = 1:length(body)
        el = body[i]
        if isa(el, GotoNode)
            body[i] = GotoNode(el.label + changemap[el.label])
        elseif isa(el, SSAValue)
            body[i] = SSAValue(el.id + changemap[el.id])
        elseif isa(el, Expr)
            if el.head === :gotoifnot
                cond = el.args[1]
                if isa(cond, SSAValue)
                    el.args[1] = SSAValue(cond.id + changemap[cond.id])
                end
                tgt = el.args[2]::Int
                el.args[2] = tgt + changemap[tgt]
            elseif el.head === :enter
                tgt = el.args[1]::Int
                el.args[1] = tgt + changemap[tgt]
            elseif !is_meta_expr_head(el.head)
                renumber_stuff!(el.args, changemap)
            end
        end
    end
end

include("compiler/ssair/driver.jl")

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
import Distributions.Bernoulli
using JLD
using ArgParse
using Random

@show q = 6
@show ρ = 2.0

let
  Random.seed!(123)

  include("BEMUtils.jl")

  λ = 1.0

  @show Threads.nthreads()

  hwv = haarCube(q)

  n_test = 100
  test_points = points([1.0, 1.0, 1.0] .- [1.0, 0.1, 0.1] .* rand(3, n_test))

  n_charges = 1000
  charges_weights = [sign.(rand(Bernoulli(0.8), n_charges) .- 0.5);
                     sign.(rand(Bernoulli(0.3), n_charges) .- 0.5)]
  charges_points = points([(rand(3, n_charges) .* [1.0, 1.0, 0.1] .+ [0.0, 0.0, 1.0])';
                           (rand(3, n_charges) .* [1.0, 0.1, 1.0] .+ [0.0, 1.0, 0.0])']')

  # Loads the ordering
  println("Loading ordering from file")
  ld = load("./out/jld/BEM_sort_q_$(q).jld")
  P = ld["P"]
  revP = ld["revP"]
  hwv = hwv[P]

  println("loading boundary data from file")
  ld = load("./out/jld/BEM_boundaryData_q_$(q).jld")
  boundaryData = ld["boundaryData"]

  colptr, rowval, P, revP, distances = sortSparseRev(hwv, ρ)

  # Change distance to ensure that no grouping into skeletons takes place.

  distances .+= 1e-12 * (1 : length(distances))
  length_rowval = length(rowval)
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  hwv = hwv[P]

  u_true = assembleMatrix(test_points, charges_points) * charges_weights

  # >>>>>>>>>>>>>>>>> Loading dense matrix and computing ref sol >>>>>>>>
  println("Loading full matrix from file")
  ld = load("./out/jld/BEM_exact_q_$(q).jld")
  Θ = ld["Θ"]
  time_assemble_dense = ld["timing_assemble_dense"]
  time_compute_dense = @elapsed u_BEM_exact = assembleMatrix(test_points, hwv) * (Θ \ boundaryData)
  time_dense = time_compute_dense + time_assemble_dense
  # <<<<<<<<<<<<<<<<< Loading dense matrix and computing the  <<<<<<<<

  # >>>>>>>>>>>>>>>>> computing prediction >>>>>>>>
  println("computing prediction from scratch")
  # time_predict = @elapsed u_predict = predict(skeletons, boundaryData, hwv, test_points, Θ)[1]
  time_predict = @elapsed u_predict = predict(skeletons, boundaryData, hwv, test_points, Θ)[1]
  # <<<<<<<<<<<<<<<<< computing prediction  <<<<<<<<

  # >>>>>>>>>>>>>>>>> computing prediction by approximating L >>>>>>>>
  println("computing prediction from scratch, by approximating L")
  @time L= assembleL(skeletons, hwv, Θ)
  u_approx_L = assembleMatrix(test_points, hwv) * (L * (L' * boundaryData))
  # <<<<<<<<<<<<<<<<< computing prediction by approximating L <<<<<<<<

  @show mean((u_predict - u_BEM_exact).^2)
  @show mean((u_approx_L - u_BEM_exact).^2)
  @show time_dense

  @show mean((u_true - u_BEM_exact).^2)
  @show mean((u_true - u_predict).^2)
  @show mean((u_true - u_approx_L).^2)
  @show time_predict

  @show 2 * length_rowval / length(hwv) / length(hwv)

  @show maximum([length(s.parents) for s in skeletons])
end


text:
export
    create,
    annihilate,
    displace,
    squeeze

############
# a† and a #
############
"""
    QuantumStateBase.Creation(T; dim::Integer)

Construct a creation operator in matrix form.

``\\hat{a}^{\\dagger}``
"""
Creation(T; dim::Integer) = .√(diagm(-1 => T.(Base.OneTo(dim-1))))

"""
    create(s::AbstractArray{T})

Apply creation operator on a quantum state.

## Arguments

* `s`: `s` can be a quantum state vector or a density matrix of a quantum state.

## Examples

```jldoctest
julia> state = VacuumState(dim=5)
5-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0

julia> create(state)
5-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0
```

```jldoctest
julia> state = VacuumState(Matrix, dim=5)
5×5 Matrix{Float64}:
 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

julia> create(state)
5×5 Matrix{Float64}:
 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  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
function create(v::AbstractVector{T}) where {T}
    dim = length(v)

    return Creation(T, dim=dim) * v
end

function create(ρ::AbstractMatrix{T}) where {T}
    dim = size(ρ, 1)
    c = Creation(T, dim=dim)

    return c * ρ * c'
end

"""
    QuantumStateBase.Annihilation(T; dim::Integer)

Construct an annihilation operator in matrix form.

``\\hat{a}``
"""
Annihilation(T; dim::Integer) = .√(diagm(1 => T.(Base.OneTo(dim-1))))

"""
    annihilate(s::AbstractArray{T})

Apply annihilation operator on a quantum state.

## Arguments

* `s`: `s` can be a quantum state vector or a density matrix of a quantum state.

## Examples

```jldoctest
julia> state = SinglePhotonState(dim=5)
5-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0

julia> annihilate(state)
5-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0
```

```jldoctest
julia> state = SinglePhotonState(Matrix, dim=5)
5×5 Matrix{Float64}:
 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  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> annihilate(state)
5×5 Matrix{Float64}:
 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
```
"""
function annihilate(v::AbstractVector{T}) where {T}
    dim = length(v)

    return Annihilation(T, dim=dim) * v
end

function annihilate(ρ::AbstractMatrix{T}) where {T}
    dim = size(ρ, 1)
    a = Annihilation(T, dim=dim)

    return a * ρ * a'
end

###########
# α and ξ #
###########

struct ComplexVec{T<:Real}
    r::T
    θ::T
end

"""
    QuantumStateBase.ComplexVec(r::T, θ::T)

Construct a vector in polar coordinate for complex plane.

``v = re^{-i\\theta}``

## Arguments

* `r`: Radius of the complex number.
* `θ`: Argument of the complex number.
"""
ComplexVec(r::T, θ::T) where {T} = ComplexVec{T}(r, θ)

function Base.show(io::IO, complexvec::ComplexVec{T}) where {T}
    print(io, "ComplexVec{$T}($(complexvec.r)exp(-$(complexvec.θ)im))")
end

"""
    QuantumStateBase.z(complexvec::ComplexVec)

Evaluate the complex vector.
"""
z(complexvec::ComplexVec) = complexvec.r * exp(-im * complexvec.θ)

################
# displacement #
################
"""
    QuantumStateBase.Displacement(T::Type{<:Complex}, r::Real, θ::Real; dim::Integer)

Construct a displacement operator in matrix form.

``\\hat{D}(\\alpha) = exp(\\alpha \\hat{a}^{\\dagger} - \\alpha^{*} \\hat{a})``
"""
function Displacement(T::Type{<:Complex}, r::Real, θ::Real; dim)
    U = real(T)
    α = ComplexVec(U(r), U(θ))

    return exp(
        z(α) * Creation(U, dim=dim) - z(α)' * Annihilation(U, dim=dim)
    )
end

"""
    displace(s::AbstractArray{T}, r, θ)

Apply displacement operator on a quantum state.

## Arguments

* `s`: `s` can be a quantum state vector or a density matrix of a quantum state.
* `r`: Radius of displacement.
* `θ`: Phase of displacement.

## Examples

```julia-repl
julia> state = VacuumState(dim=5)
5-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0

julia> displace(state, 2, π/4)
5-element Vector{ComplexF64}:
   0.14864196494651452 - 5.551115123125783e-17im
   0.15309939720606983 - 0.15309939720606955im
 8.326672684688674e-17 - 0.5201876236174409im
   -0.1359447277638306 - 0.13594472776383068im
   -0.7896009204837534 - 3.0531133177191805e-16im
```

```julia-repl
julia> state = VacuumState(Matrix, dim=5)
5×5 Matrix{Float64}:
 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

julia> displace(state, 2, π/4)
5×5 Matrix{ComplexF64}:
   0.0220944+0.0im              0.022757+0.022757im   4.12531e-17+0.0773217im    -0.0202071+0.0202071im    -0.117368+8.92137e-17im
    0.022757-0.022757im        0.0468789+0.0im          0.0796404+0.0796404im  -2.42861e-17+0.0416261im    -0.120887+0.120887im
 4.12531e-17-0.0773217im       0.0796404-0.0796404im     0.270595+0.0im           0.0707168+0.0707168im  9.30717e-17+0.410741im
  -0.0202071-0.0202071im    -2.42861e-17-0.0416261im    0.0707168-0.0707168im     0.0369619+0.0im           0.107342+0.107342im
   -0.117368-8.92137e-17im     -0.120887-0.120887im   9.30717e-17-0.410741im       0.107342-0.107342im       0.62347+0.0im
```
"""
function displace(v::AbstractVector{T}, r, θ) where {T}
    dim = length(v)

    return Displacement(complex(T), r, θ, dim=dim) * v
end

function displace(ρ::AbstractMatrix{T}, r, θ) where {T}
    dim = size(ρ, 1)
    d = Displacement(complex(T), r, θ, dim=dim)

    return d * ρ * d'
end

#############
# squeezing #
#############

"""
    QuantumStateBase.Squeezing(T::Type{<:Complex}, r::Real, θ::Real; dim::Integer)

Construct a squeezing operator in matrix form.

## Arguments

* `r`: Squeezing level.
* `θ`: Squeezing phase.

``\\hat{S}(\\xi) = exp(\\frac{1}{2}(\\xi^{*} \\hat{a}^{2} - \\xi \\hat{a}^{\\dagger 2}))``
"""
function Squeezing(T::Type{<:Complex}, r::Real, θ::Real; dim)
    U = real(T)
    ξ = ComplexVec(U(r), U(θ))

    return exp(
        (z(ξ)' * Annihilation(U, dim=dim)^2)/2 - (z(ξ) * Creation(U, dim=dim)^2)/2
    )
end

"""
    squeeze(s::AbstractArray{T})

Apply squeezing operator on a quantum state.

## Arguments

* `s`: `s` can be a quantum state vector or a density matrix of a quantum state.
* `r`: Squeezing level.
* `θ`: Squeezing phase.

## Examples

```julia-repl
julia> state = VacuumState(dim=5)
5-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0

julia> squeeze(state, 0.5, π/4)
5-element Vector{ComplexF64}:
     0.9419264274916698 + 0.0im
                    0.0 + 0.0im
   -0.21510388914370943 + 0.21510388914370937im
                    0.0 + 0.0im
 3.0228557826221164e-17 - 0.1422506201859301im
```

```julia-repl
julia> state = VacuumState(Matrix, dim=5)
5×5 Matrix{Float64}:
 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

julia> squeeze(state, 0.5, π/4)
5×5 Matrix{ComplexF64}:
    0.887225+0.0im       0.0+0.0im   -0.202612-0.202612im   0.0+0.0im  2.84731e-17+0.13399im
         0.0+0.0im       0.0+0.0im         0.0+0.0im        0.0+0.0im          0.0+0.0im
   -0.202612+0.202612im  0.0+0.0im   0.0925394+0.0im        0.0+0.0im   -0.0305987-0.0305987im
         0.0+0.0im       0.0+0.0im         0.0+0.0im        0.0+0.0im          0.0+0.0im
 2.84731e-17-0.13399im   0.0+0.0im  -0.0305987+0.0305987im  0.0+0.0im    0.0202352+0.0im
```
"""
function squeeze(v::AbstractVector{T}, r, θ) where {T}
    dim = length(v)

    return Squeezing(complex(T), r, θ, dim=dim) * v
end

function squeeze(ρ::AbstractMatrix{T}, r, θ) where {T}
    dim = size(ρ, 1)
    s = Squeezing(complex(T), r, θ, dim=dim)

    return s * ρ * s'
end


text:
const _GW_TASK_RUNME_FILE_NAME = "runme.jl"
_runme_file(taskdir::String) = joinpath(taskdir, _GW_TASK_RUNME_FILE_NAME)
_runme_file(gwt::GWTask) = _runme_file(task_dir(gwt))

## ------------------------------------------------------
function _write_runme(gwt::GWTask)
    rfile = _runme_file(gwt)
    write(rfile,
        string(
            "# It is the responsibility of the caller to use an appropriate environment", "\n",
            "import GitWorkers", "\n",
            "GitWorkers._runme(@__DIR__)",
        )
    )
end

## ------------------------------------------------------
# Try to find a GitWorker config file
function _resolve_gwtask(taskdir)
    gwtfile = _task_file(taskdir)
    tid = get(_read_toml(gwtfile), _GW_TASK_TID_KEY, "")
    return isempty(tid) ?
        GWTask(basename(taskdir), taskdir) :
        GWTask(tid, taskdir)
end

## ------------------------------------------------------
# runme

function _runme(taskdir::String)

    @eval begin

        # TODO: activate a temp env
        # A copy of a pre-generated one

        ## ------------------------------------------------------
        const __GWM = GitWorkers

        ## ------------------------------------------------------
        # GWTASK
        const __GW_TASK_DIR = $(taskdir)
        const __GWT = __GWM._resolve_gwtask(__GW_TASK_DIR)

        ## ------------------------------------------------------
        # HANDLE FlAGS
        __GWM._parse_args!(__GWT, ARGS)

        ## ------------------------------------------------------
        # MOVE TO HOMEDIR
        cd(homedir())

        ## ------------------------------------------------------
        # WELCOME
        __GWM._print_welcome(__GWT)

        ## ------------------------------------------------------
        # TASK STATUS
        @async __GWM._run_task_os(__GWT)
        # __GWM._run_task_os(__GWT)

        ## ------------------------------------------------------
        # RUN EXPR
        __GWM._fatal_err(__GWT) do
            __GWM._read_task_dat!(__GWT)
            __GW_EXPR = _task_expr(__GWT)
            isnothing(__GW_EXPR) && error("FATAL ERROR: expr object missing!")
            eval(__GW_EXPR) # run
            __GWM._flush()
        end

        ## ------------------------------------------------------
        # LABEL AS DONE
        _up_task_status!(__GWT, _GW_TASK_DONE_STATUS)

        ## ------------------------------------------------------
        # SAY GOOD BY
        __GWM._print_eotask(__GWT)

        ## ------------------------------------------------------
        # flush like a crazy
        for _ in 1:5
            __GWM._flush()
            sleep(0.05)
        end

        ## ------------------------------------------------------
        # EXIT
        exit()
    end # @eval begin
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
module Stheno

    using Distributions, Distances, BlockArrays, Statistics, Random,
        LinearAlgebra, Zygote
    import Base: length, map
    import Base.Broadcast: broadcasted, materialize, broadcast_shape
    import Statistics: mean, cov
    using LinearAlgebra: AbstractTriangular
    using ZygoteRules: @adjoint
    using Zygote: @nograd
    using BlockArrays: _BlockArray
    import LinearAlgebra: cholesky, cross
    import Distances: pairwise, colwise

    const AV{T} = AbstractVector{T}
    const AM{T} = AbstractMatrix{T}
    const AVM{T} = AbstractVecOrMat{T}

    const BlockLowerTriangular{T} = LowerTriangular{T, <:BlockMatrix{T}}
    const BlockUpperTriangular{T} = UpperTriangular{T, <:BlockMatrix{T}}
    const BlockTriangular{T} = Union{BlockLowerTriangular{T}, BlockUpperTriangular{T}}

    function elementwise end

    const pw = pairwise
    const ew = elementwise

    # Various bits of utility that aren't inherently GP-related. Often very type-piratic.
    include(joinpath("util", "zygote_rules.jl"))
    include(joinpath("util", "covariance_matrices.jl"))
    include(joinpath("util", "block_arrays", "dense.jl"))
    include(joinpath("util", "block_arrays", "diagonal.jl"))
    include(joinpath("util", "block_arrays", "triangular.jl"))
    include(joinpath("util", "abstract_data_set.jl"))
    include(joinpath("util", "distances.jl"))
    include(joinpath("util", "proper_type_piracy.jl"))

    # Supertype for GPs.
    include("abstract_gp.jl")

    # Atomic GP objects.
    include(joinpath("gp", "mean.jl"))
    include(joinpath("gp", "kernel.jl"))
    include(joinpath("gp", "gp.jl"))

    # Composite GPs, constructed via affine transformation of CompositeGPs and GPs.
    include(joinpath("composite", "composite_gp.jl"))
    include(joinpath("composite", "indexing.jl"))
    include(joinpath("composite", "cross.jl"))
    include(joinpath("composite", "conditioning.jl"))
    include(joinpath("composite", "approximate_conditioning.jl"))
    include(joinpath("composite", "product.jl"))
    include(joinpath("composite", "addition.jl"))
    include(joinpath("composite", "compose.jl"))
    # include(joinpath("composite", "gradient.jl"))
    # include(joinpath("composite", "integrate.jl"))

    # Various stuff for convenience.
    include(joinpath("util", "model.jl"))
    include(joinpath("util", "plotting.jl"))
end # module


text:
module NAssets
using Agents
using Plots
# using MetaGraphs
using LightGraphs, MetaGraphs, SimpleWeightedGraphs
using LightGraphs: YenState, complete_graph, add_edge!
using GraphPlot, GraphRecipes, NetworkLayout
# using AgentsPlots
using Tables, DataFrames
using CSV, JSON, Serialization, DelimitedFiles
using BritishNationalGrid
using ZipFile, Shapefile
using Random
using Match
using LinearAlgebra
using StatsBase
using StatsPlots
using SparseArrays
using DataStructures
using RollingFunctions
using BenchmarkTools
using Statistics
using Logging,LoggingExtras, LoggingFacilities
using Dates,TimeZones
using Distributions: Exponential, Uniform, Normal
using PyCall

#Core structures
include("core/core_structs.jl")
#graph-related
include("core/graph_functions.jl")
#events
include("eve/artificial_events.jl")
#various util functions
include("utils/util_functions.jl")
#logging
include("utils/logging_functions.jl")
#Plotting functions
include("utils/plotting_functions.jl")
# time-to-event functions
include("utils/tte_functions.jl")
#running functions
include("utils/running_functions.jl")

#maintenance model
include("phy/maintenance_model.jl")

include("ntw/of_switch.jl")
include("ntw/of_control.jl")
include("ctl/agent_control.jl")

include("phy/physical_model.jl")

include("ntw/network_model.jl")
include("phy/geo_model.jl")

#Main Functions
include("model/netManFunctions.jl")

#Agents.jl function implementation for this model
include("model/netManModel.jl")

#Basic queries
#include("model/queries_basic.jl")
#Multiple queries
include("ctl/queries_multiple.jl")

export load_run_configs, single_run_with_logging

end # module


text:
dir_name = "sample1"
l1fq1 = joinpath(dir_name, "l1fq1.fq.gz")
l1fq2 = joinpath(dir_name, "l1fq2.fq.gz")
l2fq1 = joinpath(dir_name, "l2fq1.fq.gz")
l2fq2 = joinpath(dir_name, "l2fq2.fq.gz")
files = [l1fq1, l1fq2, l2fq1, l2fq2]
try
    mkdir(dir_name)
    map(touch, reverse(files))
    d = @task dataprovider(dir_name)
    i = 1
    while true
        fq = consume(d)
        if fq == "END"
            break
        end
        fq1, fq2 = fq
        fq1 = files[i]
        fq2 = files[i+1]
        i += 2
    end
finally
    rm(dir_name, recursive=true)
end


text:
using CausalInference
using LightGraphs
using Test

include(joinpath("..", "docs", "make.jl"))

include("fci.jl")
include("klentropy.jl")
include("skeleton.jl")
include("dsep.jl")
include("pc.jl")
include("cpdag.jl")


text:
module AdventOfCode2020
using Printf: @sprintf

@info "Advent of Code 2020, Dan Burns"
@info ""

for day in 25:-1:1
    filename = @sprintf("day_%02d.jl", day)
    path = joinpath(@__DIR__, filename)
    isfile(path) && println("Loading file: ", filename)
    isfile(path) && include(path)
end

end


text:
using SafeTestsets

@safetestset "MMA approximation" begin include("approximation.jl") end
@safetestset "MMA algorithms" begin include("mma.jl") end


text:
module Netpbm

using FileIO, ImageCore

# Note: there is no endian standard, but netpbm is big-endian
const is_little_endian = ENDIAN_BOM == 0x04030201
const ufixedtype = Dict(10=>N6f10, 12=>N4f12, 14=>N2f14, 16=>N0f16)

function load(f::Union{File{format"PBMBinary"},File{format"PGMBinary"},File{format"PPMBinary"}})
    open(f) do s
        skipmagic(s)
        load(s)
    end
end

function load(f::Union{File{format"PBMText"},File{format"PGMText"},File{format"PPMText"}})
    open(f) do s
        skipmagic(s)
        load(s)
    end
end

function load(s::Stream{format"PBMBinary"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    dat = BitArray(undef,w, h)
    nbytes_per_row = ceil(Int, w/8)
    for irow = 1:h, j = 1:nbytes_per_row
        tmp = read(io, UInt8)
        offset = (j-1)*8
        for k = 1:min(8, w-offset)
            dat[offset+k, irow] = (tmp>>>(8-k))&0x01
        end
    end
    PermutedDimsArray(dat, (2,1))
end

function load(s::Stream{format"PGMBinary"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    maxval = parse_netpbm_maxval(io)
    if maxval <= 255
        dat8 = Array{UInt8}(undef, h, w)
        readio!(io, dat8)
        return reinterpret(Gray{N0f8}, dat8)
    elseif maxval <= typemax(UInt16)
        datraw = Array{UInt16}(undef, h, w)
        readio!(io, datraw)
        # Determine the appropriate Normed type
        T = ufixedtype[ceil(Int, log2(maxval)/2)<<1]
        return reinterpret(Gray{T}, datraw)
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
end

function load(s::Stream{format"PPMBinary"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    maxval = parse_netpbm_maxval(io)
    local dat
    if maxval <= 255
        dat8 = Array{UInt8}(undef, 3, h, w)
        readio!(io, dat8)
        return reshape(reinterpret(RGB{N0f8}, dat8), (h, w))
    elseif maxval <= typemax(UInt16)
        datraw = Array{UInt16}(undef, 3, h, w)
        readio!(io, datraw)
        # Determine the appropriate Normed type
        T = ufixedtype[ceil(Int, log2(maxval)/2)<<1]
        return reshape(reinterpret(RGB{T}, datraw), (h, w))
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
end

function load(s::Stream{format"PBMText"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    dat = BitArray(undef, h, w)
    nbytes_per_row = ceil(Int, w/8)
    readtextio!(io, dat)
    return(dat)
end

function load(s::Stream{format"PGMText"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    maxval = parse_netpbm_maxval(io)
    if maxval <= 255
        dat8 = Array{UInt8}(undef, h, w)
        readtextio!(io, dat8)
        return reinterpret(Gray{N0f8}, dat8)
    elseif maxval <= typemax(UInt16)
        datraw = Array{UInt16}(undef, h, w)
        readtextio!(io, datraw)
        # Determine the appropriate Normed type
        T = ufixedtype[ceil(Int, log2(maxval)/2)<<1]
        return reinterpret(Gray{T}, datraw)
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
end

function load(s::Stream{format"PPMText"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    maxval = parse_netpbm_maxval(io)
    local dat
    if maxval <= 255
        dat8 = Array{UInt8}(undef, 3, h, w)
        readtextio!(io, dat8)
        return reshape(reinterpret(RGB{N0f8}, dat8), (h, w))
    elseif maxval <= typemax(UInt16)
        datraw = Array{UInt16}(undef, 3, h, w)
        readtextio!(io, datraw)
        # Determine the appropriate Normed type
        T = ufixedtype[ceil(Int, log2(maxval)/2)<<1]
        return reshape(reinterpret(RGB{T}, datraw), (h, w))
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
end

@noinline function readio!(io, dat::AbstractMatrix{T}) where {T}
    axh, axw = axes(dat)
    for i = axh, j = axw  # io is stored in row-major format
        dat[i,j] = default_swap(read(io, T))
    end
    dat
end

@noinline function readio!(io, dat::AbstractArray{T,3}) where {T}
    size(dat, 1) == 3 || throw(DimensionMismatch("must be of size 3 in first dimension, got $(size(dat, 1))"))
    axh, axw = axes(dat, 2), axes(dat, 3)
    for i = axh, j = axw, k = 1:3  # io is stored row-major, color-first
        dat[k,i,j] = default_swap(read(io, T))
    end
    dat
end

@noinline function readtextio!(io, dat::AbstractMatrix{T}) where {T}
    axh, axw = axes(dat)
    for i = axh, j = axw  # io is stored in row-major format
        dat[i,j] = parsenextint(io)
    end
    dat
end

@noinline function readtextio!(io, dat::AbstractArray{T,3}) where {T}
    size(dat, 1) == 3 || throw(DimensionMismatch("must be of size 3 in first dimension, got $(size(dat, 1))"))
    axh, axw = axes(dat, 2), axes(dat, 3)
    for i = axh, j = axw, k = 1:3  # io is stored row-major, color-first
        dat[k,i,j] = parsenextint(io)
    end
    dat
end

function save(filename::File{format"PBMBinary"}, img; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P4\n")
        write(io, "# pbm file written by Julia\n")
        save(s, img, mapf=mapf)
    end
end

function save(filename::File{format"PGMBinary"}, img; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P5\n")
        write(io, "# pgm file written by Julia\n")
        save(s, img, mapf=mapf)
    end
end

function save(filename::File{format"PPMBinary"}, img; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P6\n")
        write(io, "# ppm file written by Julia\n")
        save(s, img, mapf=mapf)
    end
end

function save(filename::File{format"PBMText"}, img; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P1\n")
        write(io, "# pbm file written by Julia\n")
        save(s, img, mapf=mapf)
    end
end

function save(filename::File{format"PGMText"}, img; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P2\n")
        write(io, "# pgm file written by Julia\n")
        save(s, img, mapf=mapf)
    end
end

function save(filename::File{format"PPMText"}, img; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P3\n")
        write(io, "# ppm file written by Julia\n")
        save(s, img, mapf=mapf)
    end
end

function save(s::Stream, img::AbstractMatrix; mapf=identity, mapi=nothing)
    mapf = kwrename(:mapf, mapf, :mapi, mapi, :save)
    save(s, img, mapf)
end

@noinline function save(s::Stream{format"PBMBinary"}, img::AbstractMatrix{T}, mapf) where {T<:Number}
    axh, axw = axes(img)
    w = length(axw)
    write(s, "$(w) $(length(axh))\n")
    dat = PermutedDimsArray(img, (2,1))
    nbytes_per_row = ceil(Int, w/8)
    for irow = axh, j = 1:nbytes_per_row
        tmp = UInt8(0)
        offset = (j-1)*8
        for k = 1:min(8, w-offset)
            tmp |= round(UInt8, mapf(dat[offset+k, irow]) << (8-k))
        end
        write(s, UInt8(tmp))
    end
    nothing
end

@noinline function save(s::Stream{format"PGMBinary"}, img::AbstractMatrix{T}, mapf) where {T<:Union{Gray,Number}}
    axh, axw = axes(img)
    Tout, mx = pnmmax(img)
    if sizeof(Tout) > 2
        error("element type $Tout (from $T) not supported")
    end
    write(s, "$(length(axw)) $(length(axh))\n$mx\n")
    for i = axh, j = axw  # s is stored in row-major format
        write(s, default_swap(round(Tout, mx*gray(mapf(img[i,j])))))
    end
    nothing
end

@noinline function save(s::Stream{format"PPMBinary"}, img::AbstractMatrix{T}, mapf) where {T<:Color}
    axh, axw = axes(img)
    Tout, mx = pnmmax(img)
    if sizeof(Tout) > 2
        error("element type $Tout (from $T) not supported")
    end
    write(s, "$(length(axw)) $(length(axh))\n$mx\n")
    for i = axh, j = axw  # io is stored row-major, color-first
        c = RGB(mapf(img[i,j]))
        write(s, default_swap(round(Tout, mx*red(c))))
        write(s, default_swap(round(Tout, mx*green(c))))
        write(s, default_swap(round(Tout, mx*blue(c))))
    end
    nothing
end

@noinline function save(s::Stream{format"PBMText"}, img::AbstractMatrix{T}, mapf) where {T<:Number}
    axh, axw = axes(img)
    w = length(axw)
    write(s, "$(w) $(length(axh))\n")
    linelen = offset = 0
    maxdigits = 1
    for i = axh, j = axw  # s is stored in row-major format
        if linelen + maxdigits + 1 > 70
            write(s, "\n")
            linelen = 0
        elseif linelen > 0
            write(s, " ")
            linelen += 1
        end
        write(s, Bool(mapf(round(img[i,j]))) ? "1" : "0")
        linelen += maxdigits
        offset += 1
        if offset == w
            write(s, "\n")
            linelen = 0
        end
    end
    nothing
end

@noinline function save(s::Stream{format"PGMText"}, img::AbstractMatrix{T}, mapf) where {T<:Union{Gray,Number}}
    axh, axw = axes(img)
    Tout, mx = pnmmax(img)
    if sizeof(Tout) > 2
        error("element type $Tout (from $T) not supported")
    end
    w = length(axw)
    write(s, "$(w) $(length(axh))\n$mx\n")
    linelen = offset = 0
    maxdigits = round(Int, log10(mx + 1), RoundUp)
    for i = axh, j = axw  # s is stored in row-major format
        value = round(Tout, mx*gray(mapf(img[i,j])))
        if linelen + maxdigits + 1 > 70
            write(s, "\n")
            linelen = 0
        elseif linelen > 0
            write(s, " ")
            linelen += 1
        end
        write(s, lpad("$(Int(value))", maxdigits))
        linelen += maxdigits
        offset += 1
        if offset == w
            write(s, "\n")
            linelen = 0
        end
    end
    nothing
end

@noinline function save(s::Stream{format"PPMText"}, img::AbstractMatrix{T}, mapf) where {T<:Color}
    axh, axw = axes(img)
    Tout, mx = pnmmax(img)
    if sizeof(Tout) > 2
        error("element type $Tout (from $T) not supported")
    end
    w = length(axw)
    write(s, "$(w) $(length(axh))\n$mx\n")
    linelen = offset = 0
    maxdigits = round(Int, log10(mx + 1), RoundUp)
    for i = axh, j = axw  # io is stored row-major, color-first
        c = RGB(mapf(img[i,j]))
        r,g,b = round.(Tout, mx.*(red(c),green(c),blue(c)))
        if linelen + 3*maxdigits + 1 > 70
            write(s, "\n")
            linelen = 0
        elseif linelen > 0
            write(s, "  ")
            linelen += 1
        end
        write(s, string(lpad("$(Int(r))", maxdigits), ' ', lpad("$(Int(g))", maxdigits), ' ', lpad("$(Int(b))", maxdigits)))
        linelen += 2+3*maxdigits
        offset += 1
        if offset == w
            write(s, "\n")
            linelen = 0
        end
    end
    nothing
end

function parse_netpbm_size(stream::IO)
    (parsenextint(stream), parsenextint(stream))
end

function parse_netpbm_maxval(stream::IO)
    parsenextint(stream)
end

function parsenextint(stream::IO)
    # ikirill: ugly, but I can't figure out a better way
    skipchars(isspace, stream, linecomment='#')
    from = position(stream)
    mark(stream)
    skipchars(isdigit, stream)
    to = position(stream)
    reset(stream)
    parse(Int, String(read(stream, to-from+1)))
end

function pnmmax(img::AbstractArray{T}) where {T}
    if isconcretetype(T)
        return pnmmax(eltype(T))
    end
    # Determine the concrete type that can hold all the elements
    S = typeof(first(img))
    for val in img
        S = promote_type(S, typeof(val))
    end
    pnmmax(eltype(S))
end

pnmmax(::Type{T}) where {T<:AbstractFloat} = UInt8, 255
function pnmmax(::Type{U}) where {U<:Normed}
    FixedPointNumbers.rawtype(U), reinterpret(one(U))
end
pnmmax(::Type{T}) where {T<:Unsigned} = T, typemax(T)

mybswap(i::Integer)  = bswap(i)
mybswap(i::Normed)   = bswap(i)
mybswap(c::Colorant) = mapc(bswap, c)
mybswap(c::RGB24) = c

const default_swap = is_little_endian ? mybswap : identity

function kwrename(newname, newval, oldname, oldval, caller::Symbol)
    if oldval !== nothing
        Base.depwarn("keyword $oldname has been renamed $newname", caller)
        return oldval
    end
    newval
end

end # module


text:
using Test

include("../../src/utils/merge.jl")

@testset "utils/merge" begin
    @testset "merge_cliques" begin
        cliques = [[1, 2, 3],
                   [2, 3, 4],
                   [3, 4, 5],
                   [3, 4, 6],
                   [3, 6, 5],
                   [4, 6, 5],
                   [1, 2, 7],
                   [7, 8, 9],
                   [9, 10, 11]]
        cliquetree = [[1, 2],
                      [2, 3],
                      [3, 5],
                      [2, 4],
                      [4, 6],
                      [1, 7],
                      [7, 8],
                      [8, 9]]
        nb_cliques = 8
        new_cliques, new_cliquetree = merge_cliques(cliques, cliquetree, [1, 2])
        new_cliques = map(sort, new_cliques)
        new_cliquetree = map(sort, new_cliquetree)
        cliques_cliquetree = [[new_cliques[i], new_cliques[j]] for (i, j) in new_cliquetree]
        cliques_cliquetree = map(x -> map(sort, x), cliques_cliquetree)
        cliques_cliquetree = map(sort, cliques_cliquetree)
        sort!(cliques_cliquetree)
        @test cliques_cliquetree == [[[1, 2, 3, 4], [1, 2, 7]],
                                     [[1, 2, 3, 4], [3, 4, 5]],
                                     [[1, 2, 3, 4], [3, 4, 6]],
                                     [[1, 2, 7], [7, 8, 9]],
                                     [[3, 4, 5], [3, 5, 6]],
                                     [[3, 4, 6], [4, 5, 6]],
                                     [[7, 8, 9], [9, 10, 11]]]
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using PostgresQoL, Test

@test 1 == 1


text:
using ModelingToolkitStandardLibrary.Magnetic, ModelingToolkit, OrdinaryDiffEq, Test

@parameters t
@named ground = Ground()

@info "Testing basic magnetic components..."


text:
"""
    zero_one_pack!(capacity::N, weights::V, values::V, dp::V) where {N <: Number,V <: AbstractVector}
This does 0-1 (each item can be chosen only once) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values
dp array is what the function works on
It returns the ans (dp[capacity])
```
julia> dp=zeros(Int,30)
julia> zero_one_pack!(20,[1,3,11],[2,5,30],dp)
37
```
"""
function zero_one_pack!(
    capacity::N,
    weights::V,
    values::V,
    dp::V,
) where {N<:Number,V<:AbstractVector}
    for i in 1:length(weights)
        j = capacity
        while j > weights[i] # reversed loop
            dp[j] = max(dp[j], dp[j-weights[i]] + values[i])
            j -= 1
        end
        dp[weights[i]] = max(dp[weights[i]], values[i]) # dp[j]=max(dp[j],dp[0]+values[i])
    end
    return dp[capacity]
end

"""
This does complete/infinite (each item can be chosen infinite times) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values
dp array is what the function works on
It returns the ans (dp[capacity])
```
julia> dp=zeros(Int,30)
julia> complete_pack!(20,[1,2,9],[1,3,20],dp)
43
```
"""
function complete_pack!(
    capacity::N,
    weights::V,
    values::V,
    dp::V,
) where {N<:Number,V<:AbstractVector}
    for i in 1:length(weights)
        dp[weights[i]] = max(dp[weights[i]], values[i]) # dp[j]=max(dp[j],dp[0]+values[i])
        for j in (weights[i]+1):capacity
            dp[j] = max(dp[j], dp[j-weights[i]] + values[i])
        end
    end
    return dp[capacity]
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
    struct TemplateParameterList <: Any
Hold a pointer to a `clang::TemplateParameterList` object.
"""
struct TemplateParameterList
    ptr::CXTemplateParameterList
end

Base.unsafe_convert(::Type{CXTemplateParameterList}, x::TemplateParameterList) = x.ptr
Base.cconvert(::Type{CXTemplateParameterList}, x::TemplateParameterList) = x

"""
    struct TemplateArgumentList <: Any
Hold a pointer to a `clang::TemplateArgumentList` object.
"""
struct TemplateArgumentList
    ptr::CXTemplateArgumentList
end

Base.unsafe_convert(::Type{CXTemplateArgumentList}, x::TemplateArgumentList) = x.ptr
Base.cconvert(::Type{CXTemplateArgumentList}, x::TemplateArgumentList) = x

"""
    abstract type AbstractTemplateDecl <: AbstractNamedDecl
Supertype for `TemplateDecl`s.
"""
abstract type AbstractTemplateDecl <: AbstractNamedDecl end

"""
    struct TemplateDecl <: AbstractTemplateDecl
Hold a pointer to a `clang::TemplateDecl` object.
"""
struct TemplateDecl <: AbstractTemplateDecl
    ptr::CXTemplateDecl
end

Base.unsafe_convert(::Type{CXTemplateDecl}, x::TemplateDecl) = x.ptr
Base.cconvert(::Type{CXTemplateDecl}, x::TemplateDecl) = x

"""
    abstract type AbstractRedeclarableTemplateDecl <: AbstractTemplateDecl
Supertype for `RedeclarableTemplateDecl`s.
"""
abstract type AbstractRedeclarableTemplateDecl <: AbstractTemplateDecl end

"""
    struct RedeclarableTemplateDecl <: AbstractRedeclarableTemplateDecl
Hold a pointer to a `clang::RedeclarableTemplateDecl` object.
"""
struct RedeclarableTemplateDecl <: AbstractRedeclarableTemplateDecl
    ptr::CXRedeclarableTemplateDecl
end

Base.unsafe_convert(::Type{CXRedeclarableTemplateDecl}, x::RedeclarableTemplateDecl) = x.ptr
Base.cconvert(::Type{CXRedeclarableTemplateDecl}, x::RedeclarableTemplateDecl) = x

"""
    abstract type AbstractClassTemplateDecl <: AbstractRedeclarableTemplateDecl
Supertype for `ClassTemplateDecl`s.
"""
abstract type AbstractClassTemplateDecl <: AbstractRedeclarableTemplateDecl end

"""
    struct ClassTemplateDecl <: AbstractClassTemplateDecl
Hold a pointer to a `clang::ClassTemplateDecl` object.
"""
struct ClassTemplateDecl <: AbstractClassTemplateDecl
    ptr::CXClassTemplateDecl
end

Base.unsafe_convert(::Type{CXClassTemplateDecl}, x::ClassTemplateDecl) = x.ptr
Base.cconvert(::Type{CXClassTemplateDecl}, x::ClassTemplateDecl) = x

"""
    struct ClassTemplateSpecializationDecl <: AbstractCXXRecordDecl
Hold a pointer to a `clang::ClassTemplateSpecializationDecl` object.
"""
struct ClassTemplateSpecializationDecl <: AbstractCXXRecordDecl
    ptr::CXClassTemplateSpecializationDecl
end

Base.unsafe_convert(::Type{CXClassTemplateSpecializationDecl}, x::ClassTemplateSpecializationDecl) = x.ptr
Base.cconvert(::Type{CXClassTemplateSpecializationDecl}, x::ClassTemplateSpecializationDecl) = x


text:
# SPDX-License-Identifier: MIT

abstract type AbstractMixture{C} <: AbstractLatentVariableModel end

"""
    struct Mixture{C,M} <: AbstractMixture
        π
        components
    end

Mixture model. `C` is the number of  components and `M` is the type of
the components.
"""
mutable struct Mixture{C,M<:AbstractModel} <: AbstractMixture{C}
    π::P where P <: AbstractParameter
    components::ModelList{C,M}
end

basemeasure(m::Mixture, x::AbstractVector) = basemeasure(m.components[1], x)

function Mixture(T ;components, pstrength = 1)
    C = length(components)
    πprior = EFD.Dirichlet(T(pstrength) .* ones(T, C) ./ C)
    πposterior = EFD.Dirichlet(T(pstrength) .* ones(T, C) ./ C)
    π = BayesianParameter(πprior, πposterior)
    Mixture{C,eltype(components)}(π, tuple(components...))
end
Mixture(;kwargs...) = Mixture(Float64; kwargs...)

function vectorize(m::Mixture{C,M}, cache = Dict()) where {C,M}
    @cache cache lnπ = statistics(m.π)
    @cache cache component_caches = [Dict() for i in 1:C]
    vecs = [vectorize(comp, comp_cache)
            for (comp, comp_cache) in zip(m.components, component_caches)]
    vcat(hcat(vecs...), lnπ')
end

function statistics(m::Mixture, X::AbstractMatrix)
    one_const = similar(X, eltype(X), (1,size(X,2)))
    fill!(one_const, 1)
    vcat(statistics(m.components[1], X), one_const)
end

function loglikelihood(m::Mixture, X::AbstractMatrix, cache = Dict())
    TH = vectorize(m, cache)
    @cache cache TX = statistics(m, X)
    r = TH' * TX
    lnγ = r .- logsumexp(r, dims = 1)
    @cache cache γ = exp.(lnγ)
    sum(γ .* r, dims = 1)' .- sum(γ .* lnγ, dims = 1)'
end

function ∇sum_loglikelihood(m::Mixture, cache)
    TX = cache[:TX]
    γ = cache[:γ]
    acc_Tx = TX*γ'
    ∂Tπ = sum(eachcol(γ))
    grads = Dict{Any,Any}(m.π => ∂Tπ)
    for (comp, Tx, c_cache) in zip(m.components, eachcol(acc_Tx), cache[:component_caches])
        merge!(grads, ∇sum_loglikelihood(comp, Tx, c_cache))
    end
    grads
end

function posterior(m::Mixture, X::AbstractMatrix)
    TH = vectorize(m)
    TX = statistics(m, X)
    r = TH' * TX
    exp.(r .- logsumexp(r, dims = 1))
end

function predict(m::Mixture, X::AbstractMatrix)
    TH = vectorize(m)
    TX = statistics(m, X)
    r = TH' * TX
    _, maxinds = findmax(r, dims = 1)
    dropdims(getindex.(maxinds, 1), dims = 1)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
RO_data = Load_RO_data()
weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized")
weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control")

data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

treepar = TreePar()
cons = Constants()

gₛ_alt = 0.04 #mol CO₂ m⁻² s⁻¹
J_max_opt_season = 184.0e-6

A_c_F = Float64[]
Cᵢ_F = Float64[]
GPP_F = Float64[]
C_comp = Float64[]

for i = 1:length(weatherts_F.date)

    photorates = PhotoKineticRates()
    env = EnvironmentStruct(weatherts_F,i)
    photopar = PhotoPar(photorates,env.Tₐ)
    xₜ = weatherts_F.acclimation_fac[i]

    C_crit = (photopar.Γ+photopar.b_r)/(1-photopar.b_r)

    ind = Find_data_ind(i)
    LAI = data_F.Wf[ind]*data_F.N[ind]/treepar.LMA
    N = data_F.N[ind]

    daylight = weatherts_F.daylight[i]/7

    Iᵢ = CCPH.Calc_Iᵢ(env.I₀,treepar)
    gₜ = CCPH.Calc_gₜ(data_F.gₛ[i],treepar)
    Jₘₐₓ = xₜ*photopar.b_Jmax*J_max_opt_season
    photopar.α = CCPH.Calc_α(Jₘₐₓ,treepar.r_α)

    A_c,Cᵢ = CCPH.Farquhar(gₜ,Iᵢ,Jₘₐₓ,photopar,env)

    GPP = cons.M_C*(1-exp(-treepar.k*LAI))/(treepar.k)*A_c*daylight*10^3

    push!(A_c_F,A_c)
    push!(GPP_F,GPP)
    push!(Cᵢ_F,Cᵢ)
    push!(C_comp,C_crit)
end

plot(weatherts_F.date,data_F.GPP,ylabel="GPP")
plot!(weatherts_F.date,GPP_F)

@show cor(data_F.GPP,GPP_F)

plot(weatherts_F.date,A_c_F,ylabel="A_c")

plot(weatherts_F.date,data_F.gₛ,ylabel="gₛ")

plot(weatherts_F.date,Cᵢ_F,ylabel="Cᵢ")
plot!(weatherts_F.date,C_comp)
plot!([weatherts_F.date[1],weatherts_F.date[end]],[40,40])


text:
include(joinpath("..", "common", "utils.jl"))
include("cg_sem.jl")
import GaussQuadrature
using PGFPlotsX: @pgf, Plot, Axis, pgfsave, Table, LegendEntry, GroupPlot

Base.@kwdef struct Params{T}
    T_virt_surf::T = 280
    T_min_ref::T = 230
    MSLP::T = 10^5
    grav::T = 98 // 100
    R_d::T = 287
    γ::T = 7 // 5
end

# Taken from Daniel's code
function decaying_temperature_profile(params::Params, z)
    # Scale height for surface temperature
    H_sfc = params.R_d * params.T_virt_surf / params.grav
    H_t = H_sfc

    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv = params.T_virt_surf - params.T_min_ref
    Tv = params.T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / params.T_virt_surf
    p = -H_t * (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
    p /= H_sfc * (1 - ΔTv′^2)
    p = params.MSLP * exp(p)
    ρ = p / (params.R_d * Tv)
    return (p = p, ρ = ρ)
end

function plot_comparison(output_file, ρ, ρ_ref, ρe, ρe_ref, z)
    # Plot ρ_ref and ρ
    @pgf ρ_plot = Axis()
    @pgf push!(ρ_plot, Plot({color = "red", no_marks}, Table(ρ_ref[:], z[:])))
    @pgf push!(ρ_plot, LegendEntry(raw"$\rho_{ref}$"))
    @pgf push!(ρ_plot, Plot({color = "blue", no_marks}, Table(ρ[:], z[:])))
    @pgf push!(ρ_plot, LegendEntry(raw"$\rho_{hb}$"))

    # Plot ρ_ref - ρ
    @pgf Δρ_plot = Axis()
    @pgf push!(
        Δρ_plot,
        Plot({color = "red", no_marks}, Table(ρ_ref[:] - ρ[:], z[:])),
    )
    @pgf push!(Δρ_plot, LegendEntry(raw"$\rho_{ref} - \rho_{hb}$"))

    # Plot ρe_ref and ρe
    @pgf ρe_plot = Axis()
    @pgf push!(ρe_plot, Plot({color = "red", no_marks}, Table(ρe_ref[:], z[:])))
    @pgf push!(ρe_plot, LegendEntry(raw"$\rho e_{ref}$"))
    @pgf push!(ρe_plot, Plot({color = "blue", no_marks}, Table(ρe[:], z[:])))
    @pgf push!(ρe_plot, LegendEntry(raw"$\rho e_{hb}$"))

    # Plot ρe_ref - ρe
    @pgf ΔE_plot = Axis()
    @pgf push!(
        ΔE_plot,
        Plot(
            {color = "red", no_marks},
            Table((ρe_ref[:] - ρe[:]) ./ (ρe_ref[:]), z[:]),
        ),
    )
    @pgf push!(
        ΔE_plot,
        LegendEntry(raw"$(\rho e_{ref} - \rho e_{hb})/\rho e_{ref}$"),
    )

    # Make subplots
    @pgf hb_compare = GroupPlot(
        {
            group_style = {
                group_size = "2 by 2",
                horizontal_sep = "50pt",
                vertical_sep = "50pt",
            },
            legend_pos = "north east",
            ymin = minimum(z[:]),
            ymax = maximum(z[:]),
            ylabel = "height",
            # "scaled y ticks = false",
            # "scaled x ticks = false",
        },
        ρ_plot,
        ρe_plot,
        Δρ_plot,
        ΔE_plot,
    )

    pgfsave(output_file, hb_compare)
end

#=
let
    T = Float64
    N = 5
    Nq = N + 1
    K = 10
    #
    z0, z1 = 0, 30e3
    elem = element_operators(N)
    mesh = create_mesh(K, elem, z0, z1)
    #
    p = rand(T, Nq, K)
    #
    dg2cg_scatter(p, mesh, elem)
    #
    params = Params{T}()
    #
    ρ0 = -elem.D * p / (mesh.J * params.grav)
    #
    ρ = balanced_rho(p, params.grav, mesh, elem)
    #
    Pr = DG_cutoff_filter_matrix(elem.ξ, N)
    @assert all(Pr * ρ .≈ ρ0)
    ρ1 = dg2cg_scatter(copy(ρ), mesh, elem)
    @assert all(ρ .≈ ρ1)
end

let
    T = Float64
    N = 4
    Nq = N + 1
    K = 250
    #
    z0, z1 = 0, 30e4
    elem = element_operators(N)
    mesh = create_mesh(K, elem, z0, z1)
    #
    params = Params{T}()
    #
    p, ρref = similar(mesh.zdg), similar(mesh.zdg)
    for (i, z) in enumerate(mesh.zdg)
        (p[i], ρref[i]) = decaying_temperature_profile(params, z)
    end
    #
    ρ_dg = -elem.D * p / (mesh.J * params.grav)
    #
    ρ_cg = balanced_rho(p, params.grav, mesh, elem)
    Pr = DG_cutoff_filter_matrix(elem.ξ, N)
    @assert all(Pr * ρ_cg .≈ ρ_dg)
    #
    nothing
end
=#

let
    T = Float64
    N = 4
    Nq = N + 1
    K = 15
    #
    z0, z1 = 0, 30e4
    elem = element_operators(N)
    mesh = create_mesh(K, elem, z0, z1)
    #
    params = Params{T}()
    #
    p, ρ_ref = similar(mesh.zdg), similar(mesh.zdg)
    for (i, z) in enumerate(mesh.zdg)
        (p[i], ρ_ref[i]) = decaying_temperature_profile(params, z)
    end
    #
    ρ = balanced_rho(p, params.grav, mesh, elem)
    #
    w = fill!(similar(ρ), 0)
    #
    ρe_ref =
        ρ_ref .* w .^ 2 / 2 +
        p / (params.γ - 1) +
        params.grav * ρ_ref .* mesh.zdg
    ρe = ρ .* w .^ 2 / 2 + p / (params.γ - 1) + params.grav * ρ .* mesh.zdg
    #
    #=
    plot_comparison(
        "balance_compare_initial.pdf",
        ρ,
        ρ_ref,
        ρe,
        ρe_ref,
        mesh.zdg,
    )
    #
    ∂q = (
        ρ = fill!(similar(ρ), 0),
        w = fill!(similar(w), 0),
        ρe = fill!(similar(ρe), 0),
    )
    =#
    pres(ρ, w, ρe, z) =
        (params.γ - 1) * (ρe - ρ * w^2 / 2 - ρ * params.grav * z)
    #
    #
    # q = (ρ = ρ_ref, w = w, ρe = ρe_ref)
    q = (ρ = ρ, w = w, ρe = ρe)
    q0 = (ρ = copy(q.ρ), w = copy(q.w), ρe = copy(q.ρe))
    Pr = DG_cutoff_filter_matrix(elem.ξ, N)
    #=
    q_ref = (ρ = ρ_ref, w = w, ρe = ρe_ref)
    #
    element_tendency!(∂q, q_ref, elem, mesh.J, mesh.zdg, params.grav, pres)
    println("Using reference values without projection")
    @show extrema(∂q.w)
    println()
    #
    Pr = DG_cutoff_filter_matrix(elem.ξ, N)
    element_tendency!(∂q, q_ref, elem, mesh.J, mesh.zdg, params.grav, pres, Pr)
    println("Using reference values with projection")
    @show extrema(∂q.w)
    println()
    #
    #
    element_tendency!(∂q, q, elem, mesh.J, mesh.zdg, params.grav, pres)
    println("Using correct values without projection")
    @show extrema(∂q.w)
    println()
    #
    element_tendency!(∂q, q, elem, mesh.J, mesh.zdg, params.grav, pres, Pr)
    println("Using correct values with projection")
    @show extrema(∂q.w)
    println()
    =#
    #
    c = maximum(wave_speed.(ρ, pres.(q.ρ, q.w, q.ρe, mesh.zdg), params.γ))
    #
    dt = mesh.Δz / (N^2 * c)
    #
    @show extrema(q.ρ)
    rhs!(∂q, q, t) = tendency!(∂q, q, t, mesh, elem, params.grav, pres, Pr)
    ndays = 1
    tspan = (0, 86400 * ndays)
    tspan = (0, 100 * dt)
    timestep!(q, rhs!, dt, tspan)
    #
    plot_comparison("balanced.pdf", q0.ρ, q.ρ, q0.ρe, q.ρe, mesh.zdg)
    #
    @show extrema(q0.ρ)
    @show extrema(q.ρ)
    @show extrema(q.ρ - q0.ρ)
    @show extrema(q.w)
    @show sum(elem.ω .* q0.ρ) - sum(elem.ω .* q.ρ)
end
nothing


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra
using Statistics
using FiniteDiff

# ## visualizer
vis = Visualizer()
open(vis)

# ## system
include(joinpath(pathof(Dojo), "../../environments/atlas/methods/template.jl"))

function get_damper_impulses(mechanism::Mechanism{T}) where T
	joints = mechanism.joints
	# set the controls in the equality constraints
	off = 0
	nu = input_dimension(mechanism)
	u = zeros(nu)
	for joint in joints
		pbody = get_body(mechanism, joint.parent_id)
		if typeof(pbody) <: Body
			F = damper_impulses(mechanism, joint, pbody)
			oF = 0
			for joint in [joint.translational, joint.rotational]
				nf, nF = size(nullspace_mask(joint))
				u[off .+ (1:nf)] .= nullspace_mask(joint) * F[oF .+ (1:nF)]
				off += nf
				oF += nF
			end
		else
			for joint in [joint.translational, joint.rotational]
				nf, nF = size(nullspace_mask(joint))
				off += nf
			end
		end
	end
	return u
end

gravity = -9.81
timestep = 0.05
friction_coefficient = 0.8
damper = 50.0
spring = 0.0
model_type = :armless
env = get_environment(:atlas,
    representation=:minimal,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring,
	model_type=model_type)

# ## dimensions
n = env.num_states
m = env.num_inputs
d = 0

# ## reference trajectory
N = 1
initialize!(env.mechanism, :atlas)
xref = atlas_trajectory(env.mechanism;
	timestep=timestep,
	r=0.0,
	x=0.03,
	z=0.85,
	N=12,
	Ncycles=N)
zref = [minimal_to_maximal(env.mechanism, x) for x in xref]

# ## visualize reference
open(env.vis)
visualize(env, xref)

# ## gravity compensation TODO: solve optimization instead
mech = get_mechanism(:atlas,
	timestep=timestep,
	gravity=gravity,
	friction_coefficient=friction_coefficient,
	damper=1000.0,
	spring=spring,
	model_type=model_type)

initialize!(mech, :atlas)
storage = simulate!(mech, 0.1,
	record=true,
	verbose=false)

visualize(mech, storage,
	vis=env.vis)

F_damper = get_damper_impulses(env.mechanism)
u_damper = F_damper * env.mechanism.timestep
u_control = u_damper[6 .+ (1:15)]

mech = get_mechanism(:atlas,
	timestep=timestep,
	gravity=gravity,
	friction_coefficient=friction_coefficient,
	damper=0.0,
	spring=spring,
	model_type=model_type)

function controller!(mechanism, t)
    set_input!(mechanism, u_damper)
    return
end

initialize!(mech, :atlas,
	body_position=[0,0,0.0])
storage = simulate!(mech, 0.50, controller!,
	record=true,
	verbose=false)

visualize(mech, storage,
	vis=env.vis)

# ## horizon
T = N * (25 - 1) + 1

# ## model
dyn = IterativeLQR.Dynamics(
    (y, x, u, w) -> dynamics(y, env, x, u, w),
    (dx, x, u, w) -> dynamics_jacobian_state(dx, env, x, u, w),
    (du, x, u, w) -> dynamics_jacobian_input(du, env, x, u, w),
    n, n, m, d)

model = [dyn for t = 1:T-1]

# ## rollout
x1 = xref[1]
ū = [u_control for t = 1:T-1]
w = [zeros(d) for t = 1:T-1]
x̄ = IterativeLQR.rollout(model, x1, ū, w)
visualize(env, x̄)

# ## objective
qt = [0.3; 0.05; 0.05; 0.01 * ones(3); 0.01 * ones(3); 0.01 * ones(3); fill(0.002, 30)...]
ots = [(x, u, w) -> transpose(x - xref[t]) * Diagonal(timestep * qt) * (x - xref[t]) + transpose(u) * Diagonal(timestep * 0.002 * ones(m)) * u for t = 1:T-1]
oT = (x, u, w) -> transpose(x - xref[end]) * Diagonal(timestep * qt) * (x - xref[end])

cts = IterativeLQR.Cost.(ots, n, m; nw=d)
cT = IterativeLQR.Cost(oT, n, 0; nw=0)
obj = [cts..., cT]

# ## constraints
function goal(x, u, w)
    Δ = x - xref[end]
	Δ[3] -= 0.40
    return Δ[collect(1:3)]
end

cont = IterativeLQR.Constraint()
conT = IterativeLQR.Constraint(goal, n, 0)
cons = [[cont for t = 1:T-1]..., conT]

# ## solver
s = IterativeLQR.solver(model, obj, cons,
	opts=IterativeLQR.Options(
		verbose=true,
		linesearch=:armijo,
		α_min=1.0e-5,
		obj_tol=1.0e-3,
		grad_tol=1.0e-3,
		max_iter=100,
		max_al_iter=5,
		ρ_init=1.0,
		ρ_scale=10.0))
IterativeLQR.initialize_controls!(s, ū)
IterativeLQR.initialize_states!(s, x̄)

# ## solve
IterativeLQR.solve!(s)

vis= Visualizer()
open(env.vis)

# ## solution
x_sol, u_sol = IterativeLQR.get_trajectory(s)
@show IterativeLQR.eval_obj(s.m_data.obj.costs, s.m_data.x, s.m_data.u, s.m_data.w)
@show s.s_data.iter[1]
@show norm(goal(s.m_data.x[T], zeros(0), zeros(0)), Inf)

# ## visualize
x_view = [[x_sol[1] for t = 1:15]..., x_sol..., [x_sol[end] for t = 1:15]...]
visualize(env, x_view)


text:
const SemiSets{T} = Union{MOI.Semicontinuous{T}, MOI.Semiinteger{T}}

"""
    SemiToBinaryBridge{T, S <: MOI.AbstractScalarSet}

The `SemiToBinaryBridge` replaces an Semicontinuous constraint:
``x \\in \\mathsf{Semicontinuous}(l, u)``
is replaced by:
``z \\in \\{0, 1\\}``,
``x \\leq z \\cdot u ``,
``x \\geq z \\cdot l ``.

The `SemiToBinaryBridge` replaces an Semiinteger constraint:
``x \\in Semiinteger(l, u)``
is replaced by:
``z \\in \\{0, 1\\}``,
``x \\in \\Integer``,
``x \\leq z \\cdot u ``,
``x \\geq z \\cdot l ``.
"""
mutable struct SemiToBinaryBridge{T, S <: SemiSets{T}} <: AbstractBridge
    semi_set::S
    variable_index::MOI.VariableIndex
    binary_variable_index::MOI.VariableIndex
    binary_constraint_index::MOI.ConstraintIndex{MOI.SingleVariable, MOI.ZeroOne}
    lower_bound_index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}
    upper_bound_index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
    integer_index::Union{Nothing, MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer}}
end

function bridge_constraint(::Type{SemiToBinaryBridge{T,S}}, model::MOI.ModelLike, f::MOI.SingleVariable, s::S) where {T <: Real, S<:SemiSets{T}}
    binary, binary_con = MOI.add_constrained_variable(model, MOI.ZeroOne())

    # var - LB * bin >= 0
    lb = MOIU.operate(*, T, -s.lower, MOI.SingleVariable(binary))
    lb = MOIU.operate!(+, T, lb, f)
    lb_ci = MOI.add_constraint(model, lb, MOI.GreaterThan{T}(zero(T)))

    # var - UB * bin <= 0
    ub = MOIU.operate(*, T, -s.upper, MOI.SingleVariable(binary))
    ub = MOIU.operate!(+, T, ub, f)
    ub_ci = MOI.add_constraint(model, ub, MOI.LessThan{T}(zero(T)))

    if s isa MOI.Semiinteger{T}
        int_ci = MOI.add_constraint(model, f, MOI.Integer())
    else
        int_ci = nothing
    end

    return SemiToBinaryBridge{T,S}(s, f.variable, binary, binary_con, lb_ci, ub_ci, int_ci)
end

function MOIB.added_constrained_variable_types(::Type{<:SemiToBinaryBridge{T, S}}) where {T, S}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{<:SemiToBinaryBridge{T, S}}) where {T, S<:MOI.Semicontinuous{T}}
    return [
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
            ]
end

function MOIB.added_constraint_types(::Type{<:SemiToBinaryBridge{T, S}}) where {T, S <: MOI.Semiinteger{T}}
    return [
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
            (MOI.SingleVariable, MOI.Integer),
            ]
end

function concrete_bridge_type(::Type{<:SemiToBinaryBridge{T}},
                              ::Type{MOI.SingleVariable},
                              ::Type{S}) where {T, S<:SemiSets}
    return SemiToBinaryBridge{T, S}
end

function MOI.supports_constraint(::Type{<:SemiToBinaryBridge},
                                 ::Type{MOI.SingleVariable},
                                 ::Type{<:SemiSets})
    return true
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 b::SemiToBinaryBridge)
    return b.semi_set
end

function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintSet,
    bridge::SemiToBinaryBridge{T, S}, set::S) where {T, S}
    bridge.semi_set = set
    MOI.modify(model, bridge.upper_bound_index,
        MOI.ScalarCoefficientChange(bridge.binary_variable_index, -set.upper))
    MOI.modify(model, bridge.lower_bound_index,
        MOI.ScalarCoefficientChange(bridge.binary_variable_index, -set.lower))
    return
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 b::SemiToBinaryBridge{T}) where {T}
    return MOI.SingleVariable(b.variable_index)
end

function MOI.delete(model::MOI.ModelLike, bridge::SemiToBinaryBridge)
    if bridge.integer_index !== nothing
        MOI.delete(model, bridge.integer_index)
    end
    MOI.delete(model, bridge.upper_bound_index)
    MOI.delete(model, bridge.lower_bound_index)
    MOI.delete(model, bridge.binary_constraint_index)
    MOI.delete(model, bridge.binary_variable_index)
    return
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
    bridge::SemiToBinaryBridge)
    MOI.get(model, MOI.VariablePrimal(attr.N), bridge.variable_index)
end

function MOI.supports(
    ::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    ::Type{<:SemiToBinaryBridge})
    return true
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimalStart, bridge::SemiToBinaryBridge)
    return MOI.get(model, MOI.VariablePrimalStart(), bridge.variable_index)
end

function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintPrimalStart,
                 bridge::SemiToBinaryBridge{T}, value) where {T}
    MOI.set(model, MOI.VariablePrimalStart(), bridge.variable_index, value)
    bin_value = ifelse(iszero(value), 0.0, 1.0)
    MOI.set(model, MOI.VariablePrimalStart(), bridge.binary_variable_index, bin_value)
    MOI.set(model, MOI.ConstraintPrimalStart(),
        bridge.upper_bound_index, value - bridge.semi_set.upper * bin_value)
    MOI.set(model, MOI.ConstraintPrimalStart(),
        bridge.lower_bound_index, value - bridge.semi_set.lower * bin_value)
    return
end

# Attributes, Bridge acting as a model

function MOI.get(::SemiToBinaryBridge, ::MOI.NumberOfVariables)
    return 1
end

function MOI.get(b::SemiToBinaryBridge, ::MOI.ListOfVariableIndices)
    return [b.binary_variable_index]
end

function MOI.get(::SemiToBinaryBridge{T, S},
    ::MOI.NumberOfConstraints{MOI.SingleVariable, MOI.ZeroOne}) where {T, S}
    return 1
end

function MOI.get(::SemiToBinaryBridge{T, S},
    ::MOI.NumberOfConstraints{MOI.SingleVariable, MOI.Integer}) where {T, S<:MOI.Semiinteger}
    return 1
end

function MOI.get(::SemiToBinaryBridge{T, S},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}) where {T, S}
    return 1
end

function MOI.get(::SemiToBinaryBridge{T, S},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}) where {T, S}
    return 1
end

function MOI.get(b::SemiToBinaryBridge{T, S},
    ::MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.ZeroOne}) where {T, S}
    return [b.binary_constraint_index]
end

function MOI.get(b::SemiToBinaryBridge{T, S},
    ::MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.Integer}) where {T, S<:MOI.Semiinteger}
    return [b.integer_index]
end

function MOI.get(b::SemiToBinaryBridge{T, S},
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}) where {T, S}
    return [b.upper_bound_index]
end

function MOI.get(b::SemiToBinaryBridge{T, S},
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}) where {T, S}
    return [b.lower_bound_index]
end


text:
mutable struct BinaryTree
    depth::Int
    nnodes::Int
    edges::Vector{Tuple{Int,Int}}
    pathedges::Dict{Tuple{Int,Int},Vector{Tuple{Int,Int}}}
    codes::Vector{Vector{Int}}
end

function BinaryTree(depth::Int)
    nnodes = 2^(depth+1)-1
    edges = Tuple{Int,Int}[]
    allcodes = Any[[-1],[0],[1]]
    codes = 0:1
    for d in 1:(depth-1)
        codes = collect(IterTools.product(codes,0:1))
        push!(allcodes, codes)
    end
    allcodes = vcat([vec(c) for c in allcodes]...)
    function flatten(arr)
        rst = Any[]
        grep(v) =   for x in v
                    if isa(x, Tuple)
                    grep(x)
                    else push!(rst, x) end
                    end
        grep(arr)
        rst
    end
    allcodes = [reverse(flatten(code)) for code in allcodes]
    bt = BinaryTree(depth, nnodes,
        edges,
        Dict{Tuple{Int,Int},Vector{Tuple{Int,Int}}}(),
        allcodes)

    g = LightGraphs.Graph(nnodes)
    for d=0:(depth-1), n in getnodes(bt, d), c in getchildren(bt,n)
        push!(edges, (n,c))
        LightGraphs.add_edge!(g, (n,c))
    end
    fw = LightGraphs.floyd_warshall_shortest_paths(g)
    paths =  LightGraphs.enumerate_paths(fw) # [u][v]
    # (u,v) => list of edges on shortest path
    pathedges = Dict((u,v) => [getedge(bt, paths[u][v][i-1],paths[u][v][i])
                               for i in 2:length(paths[u][v])]
                     for u in 1:nnodes for v in 1:nnodes)
    bt.pathedges = pathedges
    bt
end

validnode(bt::BinaryTree, node::Int) = ((node >= 1) && (node <= bt.nnodes))
validlayer(bt::BinaryTree, layer::Int) = ((layer >= 0) && (layer <= bt.depth))

function getlayer(bt::BinaryTree, node::Int)
    @assert validnode(bt, node)
    floor(Int, log(2,node))
end

function getnodes(bt::BinaryTree, layer::Int)
    #layer 0    1
    #layer 1    2:3
    @assert validlayer(bt, layer)
    (2^layer):(2^(layer+1)-1)
end

getleaves(bt::BinaryTree) = getnodes(bt, bt.depth)

function getparent(bt::BinaryTree, node::Int)
    @assert validnode(bt, node)
    floor(Int, node/2)
end

function getchildren(bt::BinaryTree, node::Int)
    @assert validnode(bt, node)
    thislayer = getlayer(bt, node)
    thisid = findfirst(getnodes(bt, thislayer) .== node)
    nextnodes = getnodes(bt, thislayer+1)
    [nextnodes[2*thisid-1], nextnodes[2*thisid]]
end

function getdescendants(bt::BinaryTree, node::Int)
    @assert validnode(bt, node)
    thislayer = getlayer(bt, node)
    descendants = Int[]
    if thislayer < bt.depth
        descendants = getchildren(bt, node)
        for d in descendants
            dlayer = getlayer(bt, d)
            descendants = vcat(descendants, getdescendants(bt, d))
        end
    end
    return descendants
end

function getsubtreeleaves(bt::BinaryTree, node::Int)
    return intersect(vcat(node,getdescendants(bt, node)), getleaves(bt))
end

function getedge(bt::BinaryTree, u::Int, v::Int)
    @assert validnode(bt, u)
    @assert validnode(bt, v)
    return (min(u,v),max(u,v))
end

function Base.show(io::IO, bt::BinaryTree; offset::String="")
    println(io, offset, string(typeof(bt)))
    println(io, offset, "    * depth: $(bt.depth)  (nnodes: $(bt.nnodes))")
end


text:
using SignedDistanceFields
using Base.Test

T = true
F = false

@test_throws ArgumentError edf([])
@test_throws ArgumentError edf([T])

@test edf([T T]) == [0.0 0.0]
@test edf([T T F F F F]) == [0.0 0.0 1.0 2.0 3.0 4.0]
@test edf([F F T T F F]) == [2.0 1.0 0.0 0.0 1.0 2.0]
@test edf([F F F F T T]) == [4.0 3.0 2.0 1.0 0.0 0.0]
@test edf([F T F F T F]) == [1.0 0.0 1.0 1.0 0.0 1.0]

@test sdf([T T F F F F]) == [-2.0 -1.0   1.0  2.0  3.0  4.0]
@test sdf([F F T T F F]) == [ 2.0  1.0  -1.0 -1.0  1.0  2.0]
@test sdf([F F F F T T]) == [ 4.0  3.0   2.0  1.0 -1.0 -2.0]
@test sdf([F T F F T F]) == [ 1.0 -1.0   1.0  1.0 -1.0  1.0]

@test sdf([
	T F F;
	F F F;
	F F F
]) == [
	-1.0  1.0          2.0;
	 1.0  sqrt(2)      sqrt(1 + 4);
	 2.0  sqrt(1 + 4)  sqrt(4 + 4)
]

if VERSION < v"0.4.0-dev"
    bitrand(x...) = randbool(x...)
end

a = sdf(bitrand(20,20), 10, 10)
@test size(a) == (10,10)
a = edf(bitrand(20,20), 10, 10)
@test size(a) == (10,10)

if VERSION < v"0.4.0-dev"
	@test_throws ErrorException sdf(bitrand(20,20), 12, 10)
	@test_throws ErrorException edf(bitrand(20,20), 10, 12)
else
	@test_throws AssertionError sdf(bitrand(20,20), 12, 10)
	@test_throws AssertionError edf(bitrand(20,20), 10, 12)
end


text:
export Ellipsoid

mutable struct Ellipsoid{N,N²} <: Shape{N,N²}
    c::SVector{N,Float64}  # center of ellipsoid
    ri2::SVector{N,Float64}  # inverse squares of "radii" (semi-axes) in axis directions
    p::SMatrix{N,N,Float64,N²}  # projection matrix to Ellipsoid coordinates; must be orthonormal (see surfpt_nearby)
    Ellipsoid{N,N²}(c,ri2,p) where {N,N²} = new(c,ri2,p)  # suppress default outer constructor
end

Ellipsoid(c::SVector{N,<:Real},
          r::SVector{N,<:Real},
          axes::SMatrix{N,N,<:Real}=SMatrix{N,N,Float64}(I)
          ) where {N} =
    Ellipsoid{N,N*N}(c, float.(r).^-2, inv(axes ./ sqrt.(sum(abs2,axes,dims=Val(1)))))

Ellipsoid(c::AbstractVector{<:Real},  # center of ellipsoid
          r::AbstractVector{<:Real},  # ""
          axes::AbstractMatrix{<:Real}=Matrix{Float64}(I,length(c),length(c))) =  # columns are axes vector; assumed orthogonal
    (N = length(c); Ellipsoid(SVector{N}(c), SVector{N}(r), SMatrix{N,N}(axes)))

Ellipsoid(b::Cuboid{N,N²}) where {N,N²} = Ellipsoid{N,N²}(b.c, (b.r).^-2, b.p)

Base.:(==)(b1::Ellipsoid, b2::Ellipsoid) = b1.c==b2.c && b1.ri2==b2.ri2 && b1.p==b2.p
Base.isapprox(b1::Ellipsoid, b2::Ellipsoid) = b1.c≈b2.c && b1.ri2≈b2.ri2 && b1.p≈b2.p
Base.hash(b::Ellipsoid, h::UInt) = hash(b.c, hash(b.ri2, hash(b.p, hash(:Ellipsoid, h))))

level(x::SVector{N,<:Real}, b::Ellipsoid{N}) where {N} = √dot((b.p * (x - b.c)).^2, b.ri2) - 1.0

function surfpt_nearby(x::SVector{N,<:Real}, b::Ellipsoid{N}) where {N}
    if x == b.c
        _m, i = findmax(b.ri2)
        nout = b.p[i,:]  # assume b.p is orthogonal
        return b.c + nout/√b.ri2[i], nout
    end

    # For a given point x and equation of ellipsoid f(x) = 1, find t such that x₀ = x + t*∇f(x)
    # is on the ellipsoid.  Eventually this reduces to a quadratic equation for t.  The
    # following is evaluation of the quadratic formula in a numerically stable way.
    px = b.p * (x - b.c)  # in ellipsoid's coordinates
    px2 = px.^2

    px²r⁻² = px2 ⋅ b.ri2
    px²r⁻⁴ = px2 ⋅ (b.ri2.^2)
    px²r⁻⁶ = px2 ⋅ (b.ri2.^3)

    q24 = (px²r⁻² - 1) / px²r⁻⁴
    q64 = px²r⁻⁶ / px²r⁻⁴

    t = -q24 / (1 + √(1 - q24 * q64))

    # From t, recover x₀ = x + t*∇f(x).
    px₀ = (t*b.ri2 .+ 1) .* px  # surface point in ellipsoid coordinates

    # Transform back to the original coordinates.
    x₀ = b.p' * px₀ + b.c
    nout = normalize(b.p' * (px .* b.ri2))

    return x₀, nout
end

function boundpts(b::Ellipsoid{N}) where {N}
    # Return the points tangential to the bounding box.
    # For N = 3, it returns three points at which the direction normals are +x, +y, +z
    # directions, respectively.

    r2 = 1 ./ b.ri2
    ndir = b.p  # Cartesian directions in ellipsoid coordinates: b.p * I

    # In the ellipsoid coordinates, the point on the ellipsoid where the direction normal is
    # n is (n .* r2) / sqrt(r2' * n.^2).  Below is the broadcasted version of this over a
    # matrix n, whose each column is a direction normal.  Once calculated, we need to
    # change the coordinates back to the original coordinates.
    M = b.p' * ((ndir .* r2) ./ sqrt.(r2' * ndir.^2))

    return M
end

function bounds(b::Ellipsoid{N}) where {N}
    M = boundpts(b)

    # Note that when M does not include NaN, we can simply set m = diag(M), because the
    # first (second, third) column of M is the point on the ellipsoid at which the normal
    # vector is +x (+y, +z).  Therefore, the x (y, z) point of the first (second, third)
    # column has the largest x (y, z) coordinate.
    #
    # However, if one of a, b, c is zero, the shape is a disk.  Then one column of M can be
    # completely filled with NaN.  This column must not be counted as a bounding point, so
    # we apply NaN-ignoring maximum by StaticArrays.reducedim along the row direction.
    #
    # For the efficient implementation of NaN-ignoring maximum, see
    # https://discourse.julialang.org/t/inconsistency-between-findmax-and-maximum-with-respect-to-nan/4214/8
    m = reduce((x,y) -> x<y ? y : x, M, init=-Inf, dims=Val(2))[:,1]

    return (b.c-m,b.c+m)
end


text:
"""
This relatively light array wrapper allows us to access arrays with the quantum number
as the index, which means we can just type in the equations from the papers.
"""

using Base: @propagate_inbounds
using CustomUnitRanges: filename_for_urange
include(filename_for_urange)

struct WignerSymbolVector{T, Ti, AA<:AbstractArray{T,1}} <: AbstractArray{T,1}
    symbols::AA
    offset::Ti
end

Base.parent(A::WignerSymbolVector) = A.symbols
offset(axparent::AbstractUnitRange, ax::AbstractUnitRange) = first(ax) - first(axparent)

function WignerSymbolVector(A::AbstractArray{T,1}, inds::AbstractUnitRange) where {T}
    axparent = first(axes(A))
    WignerSymbolVector(A, offset(axparent, inds))
end

Base.eachindex(::IndexLinear, A::WignerSymbolVector)   = axes(A, 1)
Base.IndexStyle(::Type{<:WignerSymbolVector}) = IndexLinear()
@inline Base.size(A::WignerSymbolVector) = size(parent(A))

Base.@propagate_inbounds function Base.getindex(
        A::WignerSymbolVector{T,Ti,AA}, hi::Ti) where {T, Ti<:HalfInteger, AA}
    return parent(A)[div(hi.twice - A.offset.twice, 2)]
end
Base.@propagate_inbounds function Base.getindex(
        A::WignerSymbolVector{T,Ti,AA}, hi::Rational) where {T, Ti<:HalfInteger, AA}
    return parent(A)[div(2 * hi - A.offset.twice, 2)]
end
@propagate_inbounds Base.getindex(
    A::WignerSymbolVector{T,Ti,AA}, i::Ti) where {T, Ti<:Int, AA}  = parent(A)[i - A.offset]

Base.@propagate_inbounds function Base.setindex!(
        A::WignerSymbolVector{T,Int,AA}, val, i::Int) where {T, AA}
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[i - A.offset] = val
    val
end
Base.@propagate_inbounds function Base.setindex!(
        A::WignerSymbolVector{T,HalfInt,AA}, val, hi::HalfInt) where {T, AA}
    @boundscheck checkbounds(A, hi)
    @inbounds parent(A)[div(hi.twice - A.offset.twice, 2)] = val
    val
end

Base.axes(A::WignerSymbolVector{T,Ti,AA}) where {T, Ti, AA} =
    map(n->URange{Ti}(1+A.offset,n+A.offset), size(A))
Base.axes1(A::WignerSymbolVector{T,Ti,AA}) where {T, Ti, AA} = URange{Ti}(1+A.offset,size(A,1)+A.offset)

Base.dataids(A::WignerSymbolVector) = Base.dataids(parent(A))

function Base.show(io::IO, a::WignerSymbolVector)
    print(io, "$(length(a))-element WignerSymbolVector (")
    println(io, eltype(a), ") over ", Base.axes1(a))
    Base.print_array(io, parent(a))
end