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

module BaseTestDeprecated

export @test, @test_throws, @test_approx_eq, @test_approx_eq_eps, @inferred

abstract Result
type Success <: Result
    expr
    resultexpr
    res
    Success(expr, resultexpr=nothing, res=nothing) = new(expr, resultexpr, res)
end
type Failure <: Result
    expr
    resultexpr
end
Failure(expr) = Failure(expr, nothing)
type Error <: Result
    expr
    err
    backtrace
end

default_handler(r::Success) = r.res
function default_handler(r::Failure)
    if r.resultexpr !== nothing
        error("test failed: $(r.resultexpr)\n in expression: $(r.expr)")
    else
        error("test failed in expression: $(r.expr)")
    end
end
default_handler(r::Error) = rethrow(r)

handler() = get(task_local_storage(), :TEST_HANDLER, default_handler)

with_handler(f::Function, handler) =
    task_local_storage(f, :TEST_HANDLER, handler)

import Base.showerror

showerror(io::IO, r::Error) = showerror(io, r, [])
function showerror(io::IO, r::Error, bt)
    println(io, "test error in expression: $(r.expr)")
    showerror(io, r.err, r.backtrace)
end

function do_test(body,qex)
    handler()(try
        rex, val = body()
        val ? Success(qex, rex) : Failure(qex,rex)
    catch err
        Error(qex,err,catch_backtrace())
    end)
end

function do_test_throws(body, qex, bt, extype)
    handler()(try
        body()
        Failure(qex, "$qex did not throw $(extype === nothing ? "anything" : extype)")
    catch err
        if extype === nothing
            Base.warn("""
            @test_throws without an exception type is deprecated;
            Use `@test_throws $(typeof(err)) $(qex)` instead.
            """, bt = bt)
            Success(qex, nothing, err)
        else
            if isa(err, extype)
                Success(qex, nothing, err)
            else
                if isa(err,Type)
                    Failure(qex, "the type $err was thrown instead of an instance of $extype")
                else
                    Failure(qex, "$err was thrown instead of $extype")
                end
            end
        end
    end)
end

macro test(ex)
    if typeof(ex) == Expr && ex.head == :comparison
        syms = [gensym() for i = 1:length(ex.args)]
        func_block = Expr(:block)
        # insert assignment into a block
        func_block.args = [:($(syms[i]) = $(esc(ex.args[i]))) for i = 1:length(ex.args)]
        # finish the block with a return
        push!(func_block.args, Expr(:return, :(Expr(:comparison, $(syms...)), $(Expr(:comparison, syms...)))))
        :(do_test(()->($func_block), $(Expr(:quote,ex))))
    else
        :(do_test(()->($(Expr(:quote,ex)), $(esc(ex))), $(Expr(:quote,ex))))
    end
end

macro test_throws(args...)
    ex = nothing
    extype = nothing
    # Users should pass (ExceptionType, Expression) but we give a warning to users that only pass (Expression)
    if length(args) == 1
        ex = args[1]
    elseif length(args) == 2
        ex = args[2]
        extype = args[1]
    end
    :(do_test_throws(()->($(esc(ex))),$(Expr(:quote,ex)),backtrace(),$(esc(extype))))
end

approx_full(x::AbstractArray) = x
approx_full(x::Number) = x
approx_full(x) = full(x)

function test_approx_eq(va, vb, Eps, astr, bstr)
    va = approx_full(va)
    vb = approx_full(vb)
    if length(va) != length(vb)
        error("lengths of ", astr, " and ", bstr, " do not match: ",
              "\n  ", astr, " (length $(length(va))) = ", va,
              "\n  ", bstr, " (length $(length(vb))) = ", vb)
    end
    diff = real(zero(eltype(va)))
    for i = 1:length(va)
        xa = va[i]; xb = vb[i]
        if isfinite(xa) && isfinite(xb)
            diff = max(diff, abs(xa-xb))
        elseif !isequal(xa,xb)
            error("mismatch of non-finite elements: ",
                  "\n  ", astr, " = ", va,
                  "\n  ", bstr, " = ", vb)
        end
    end

    if !isnan(Eps) && !(diff <= Eps)
        sdiff = string("|", astr, " - ", bstr, "| <= ", Eps)
        error("assertion failed: ", sdiff,
              "\n  ", astr, " = ", va,
              "\n  ", bstr, " = ", vb,
              "\n  difference = ", diff, " > ", Eps)
    end
end

array_eps{T}(a::AbstractArray{Complex{T}}) = eps(float(maximum(x->(isfinite(x) ? abs(x) : T(NaN)), a)))
array_eps(a) = eps(float(maximum(x->(isfinite(x) ? abs(x) : oftype(x,NaN)), a)))

test_approx_eq(va, vb, astr, bstr) =
    test_approx_eq(va, vb, 1E4*length(va)*max(array_eps(va), array_eps(vb)), astr, bstr)

macro test_approx_eq_eps(a, b, c)
    :(test_approx_eq($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b))))
end

macro test_approx_eq(a, b)
    :(test_approx_eq($(esc(a)), $(esc(b)), $(string(a)), $(string(b))))
end

macro inferred(ex)
    ex.head == :call || error("@inferred requires a call expression")
    quote
        vals = ($([esc(ex.args[i]) for i = 2:length(ex.args)]...),)
        inftypes = Base.return_types($(esc(ex.args[1])), Base.typesof(vals...))
        @assert length(inftypes) == 1
        result = $(esc(ex.args[1]))(vals...)
        rettype = isa(result, Type) ? Type{result} : typeof(result)
        rettype == inftypes[1] || error("return type $rettype does not match inferred return type $(inftypes[1])")
        result
    end
end

# Test approximate equality of vectors or columns of matrices modulo floating
# point roundoff and phase (sign) differences.
#
# This function is design to test for equality between vectors of floating point
# numbers when the vectors are defined only up to a global phase or sign, such as
# normalized eigenvectors or singular vectors. The global phase is usually
# defined consistently, but may occasionally change due to small differences in
# floating point rounding noise or rounding modes, or through the use of
# different conventions in different algorithms. As a result, most tests checking
# such vectors have to detect and discard such overall phase differences.
#
# Inputs:
#     a, b:: StridedVecOrMat to be compared
#     err :: Default: m^3*(eps(S)+eps(T)), where m is the number of rows
#
# Raises an error if any columnwise vector norm exceeds err. Otherwise, returns
# nothing.
function test_approx_eq_modphase{S<:Real,T<:Real}(
        a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, err=nothing)

    m, n = size(a)
    @test n==size(b, 2) && m==size(b, 1)
    err === nothing && (err=m^3*(eps(S)+eps(T)))
    for i=1:n
        v1, v2 = a[:, i], b[:, i]
        @test_approx_eq_eps min(abs(norm(v1-v2)), abs(norm(v1+v2))) 0.0 err
    end
end

end # module


text:
using SparseArrays
using LinearAlgebra
using QuadGK

mutable struct modelParameters
    N
    c1
    c2
    K
    xf
    τ
    A1
    solDir
    T
    tMax
    solAdj
    phaseSensNorm
end

function setParams(N,c1,c2,K,xf,τ,print=false)
	D  = zeros(N,N)
	Om = zeros(N,N)
	Id = zeros(N,N)
    for i in 1:N
        D[i,i]  = 2*damping(c1,c2,i)*i*π
        Om[i,i] = i*i*π*π
        Id[i,i] = 1
    end
    A1 = [zeros(N,N) Id;-Om -D]
    A1 = sparse(A1)
    params = modelParameters(N,c1,c2,K,xf,τ,A1,false,false,false,false,false)
    if(print)
	println("#####################################")
	println("## Rijke tube Galerkin modes model ##")
	println("#####################################")
	println("\n##  Parameters ##")
	println("N  = ",params.N)
	println("c1 = ",params.c1)
	println("c2 = ",params.c2)
	println("K  = ",params.K)
	println("xf = ",params.xf)
	println("τ  = ",params.τ)
	println("#################\n")
    end
    return params
end

function damping(c1,c2,i)
    return (1/(2π))*(c1*i+c2*sqrt(1/i))
end

function flameVel(xd,p)
    uf = 0
    for i in 1:p.N
        uf = uf + xd[i]*cos(i*π*p.xf)
    end
    return uf
end

function RHS(u,h,p,t)

    uf = flameVel(h(p,t-p.τ),p)

    y = p.A1*u

    for i in 1:p.N
        y[N+i] = y[N+i] - i*π*p.K*(sqrt(abs(1/3+uf))-sqrt(1/3))*sin(i*π*p.xf)
    end
    return y
end

function pressureVelocity(modes,x,params)
    # Returns the pressure and velocity from the Galerkin modes at x
    u = zeros(length(x),)
    p = zeros(length(x),)
    for i in 1:params.N
	    u += modes[i]*cos.(i*π*x)
        p -= modes[i+params.N]*sin.(i*π*x)*γ*Ma/(i*π)
    end

    return u,p
end

function energy(modes,params)
    # Calculates the non-dimensional acoustic energy per unit volume
    # from the Galerkin modes
    E = 0
    for i=1:params.N
        E += modes[i]^2+(modes[params.N+i]/(i*π))^2
    end
    E *= 0.5
    return E
end

## Linear and Adjoint routines
function linHeatRelease(t,p)
    # Creates the time-delayed system matrix A2 from the direct solution at t.
    Q = zeros(p.N,p.N)
    z = p.solDir(t-p.τ)
    ufb = flameVel(z,p)

    for i in 1:p.N
        for j in 1:p.N
            Q[i,j]=0.5*i*π*p.K*sin(i*π*p.xf)*sign(1/3+ufb)*cos(j*π*p.xf)/sqrt(abs(1/3+ufb))
        end
    end
    A2 = -[zeros(p.N,p.N) zeros(p.N,p.N);Q zeros(p.N,p.N)]
    return A2
end

function linHeatReleaseMat(base,p)
    # Creates the linearised matrix A2 from the state base
    Q = zeros(p.N,p.N)
    z = base
    ufb = flameVel(z,p)

    for i in 1:p.N
        for j in 1:p.N
            Q[i,j]=0.5*i*π*p.K*sin(i*π*p.xf)*sign(1/3+ufb)*cos(j*π*p.xf)/sqrt(abs(1/3+ufb))
        end
    end
    A2 = -[zeros(p.N,p.N) zeros(p.N,p.N);Q zeros(p.N,p.N)]
    return A2
end

function adjTime(tAdj)
    # Allows us to solve the adjoint forwards in time
    return params.tMax .- tAdj
end

function RHSadj(u,h,p,t)
    A2 = linHeatRelease(adjTime(t-p.τ),p)
    return p.A1'*u + A2'*h(p,t-p.τ)
end

# Misc functions

function bilinearForm(t,p)
    # Inner product part
    innerProduct = dot(p.solAdj(adjTime(t)),p.solDir(t,Val{1}))
    # Integral part
    integrand(ξ) = dot(p.solAdj(adjTime(t+ξ+p.τ)),linHeatRelease(t+ξ+p.τ,p)*p.solDir(t+ξ,Val{1}))
    integral,error = quadgk(integrand, -p.τ, 0., rtol=1e-4)

    norm = innerProduct + integral
    return norm, innerProduct, integral
end

function phaseCouplingFunction(ϕ,ωf,f,p)
    Tf = 2*π/ωf
    integrand(s) = dot(p.solAdj(adjTime(ϕ+params.T*s/Tf)),f(s))
    integral,error = quadgk(integrand, 0, Tf, rtol=1e-4)
    Γ = integral / Tf
    return Γ
end

function phaseCouplingFunctionmn(ϕ,ωf,f,m,n,p)
    Tf = 2*π/ωf
    integrand(s) = dot(p.solAdj(adjTime(ϕ+(n/m)*params.T*s/Tf)),f(s))
    integral,error = quadgk(integrand, 0, m*Tf, rtol=1e-4)
    Γ = integral / (m*Tf)
    return Γ
end

ext: jl


text:
@testset "Validating a Schema descriptor" begin

    @testset "Read in from JSON" begin
        s = Schema("../data/schema_valid_infer.json")
        @test length(s.fields) == 2
        @test s.fields[1].name == "id"
        @test !s.fields[2].constraints.required
        @test TableSchema.is_valid(s) == true
    end

    @testset "Created from scratch" begin
        f = Field("width")
        f.typed = "integer"
        f.constraints.required = true
        s = Schema()
        TableSchema.add_field(s, f)
        @test length(s.fields) == 1
        @test s.fields[1].constraints.required
        @test_throws SchemaError validate(s, true)
    end

    @testset "Check primary foreign keys" begin
        s = Schema("../data/schema_valid_full.json")
        @test length(s.primary_key) == 4
        @test length(s.foreign_keys) == 1
    end

    @testset "Handle schema errors" begin
        @test_throws SchemaError Schema("../data/schema_invalid_empty.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_wrong_type.json", true)
        s = Schema("../data/schema_invalid_multiple_errors.json")
        @test !(TableSchema.is_valid(s))
        @test length(s.errors) == 5
    end

    @testset "Invalid foreign keys in schema" begin
        @test_throws SchemaError Schema("../data/schema_invalid_fk_string.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_fk_array.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_fk_array_string.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_fk_array_string_ref.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_fk_array_wrong_number.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_fk_no_reference.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_fk_string_array_ref.json", true)
    end

    @testset "Invalid primary key in schema" begin
        @test_throws SchemaError Schema("../data/schema_invalid_pk_string.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_pk_array.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_pk_is_wrong_type.json", true)
        @test_throws SchemaError Schema("../data/schema_invalid_pk_no_fields.json", true)
    end

    @testset "Detailed validation errors" begin
        # TODO see https://github.com/frictionlessdata/tableschema-py/blob/c7b2125108d82029dab7bf6a31c894497837a6ac/tests/test_validate.py#L98
        s = Schema("../data/schema_invalid_pk_string.json")
        @test !(TableSchema.is_valid(s))
        @test length(s.errors) == 1
        @test s.errors[1].key == "identifier"
    end

end
@testset "Validating Table data with Schema" begin

    @testset "Check constraints" begin
        s = Schema("../data/schema_valid_missing.json")
        t = Table("../data/data_types.csv")
        tr = TableSchema.read(t, cast=false)
        @test s.fields[1].constraints.required
        @test TableSchema.checkrow(s.fields[1], tr[1,1])
        @test TableSchema.checkrow(s.fields[2], tr[2,2])
        @test TableSchema.checkrow(s.fields[3], tr[3,3])
        @test_throws ConstraintError TableSchema.checkrow(s.fields[1], "")
    end

    @testset "Schema constraint errors" begin
        s = Schema("../data/schema_valid_missing.json")
        t = Table("../data/data_constraints.csv", s)
        @test !(validate(t))
        @test length(t.errors) > 0
        @test t.errors[2].message == "unique"
        @test t.errors[2].field.name == "name"
        @test t.errors[1].message == "required"
        @test t.errors[1].field.name == "id"
    end

end


text:
"""
plot recipe for ::MDCSolution
kwargs: pnames are array of parameter names
idxs: are parameter indices to plot
what ∈ (:trajectory, :final_changes) determines the plot type
"""
@recipe function f(mdc::MDCSolution; pnames=nothing, idxs=nothing, what=:trajectory)
    if idxs === nothing
        num = min(5, mdc.N)
        idxs = biggest_movers(mdc, num)
    end
    # if !(names === nothing)
    #     labels --> names[idxs]
    # end
    # ["hi" "lo" "lo" "hi" "lo"]

    tfirst = mdc.sol.t[1]
    tend = mdc.sol.t[end]

    layout := (1, 1)
    bottom_margin := :match

    if what == :trajectory
        @series begin
            if !(pnames === nothing)
                label -->  reshape(pnames[idxs], 1, :)
            end
            title --> "change in parameters over minimally disruptive curve"
            xguide --> "distance"
            yguide --> "Δ parameters"
            distances(mdc), Δ(mdc)[idxs,:]'
        end
    end

    if what == :final_changes
        @series begin
            title --> "biggest changers"
        seriestype := :bar
            label --> "t=$tend"
            xticks --> (1:5, reshape(pnames[idxs], 1, :))
            xrotation --> 90
            Δ(mdc, tend)[idxs]
        end
        if tfirst < 0.
            @series begin
                label --> "t=$tfirst"
                seriestype := :bar
                xticks --> (1:5, reshape(pnames[idxs], 1, :))
                xrotation --> 90
                Δ(mdc, tfirst)[idxs]
            end
        end
    end
end

"""
    output_on_curve(f, mdc, t)

Useful when building an animation of f(p) as the parameters p vary along the curve.
"""
function output_on_curve(f, mdc, t)
    return f(mdc(t)[:states])
end


text:
include("./params.jl")

function fitgamma(isi)

    ################# WRITE SOMETHING #################

    return lambda, k

end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
function readinput(filename)
    map(group -> filter(!isempty, split(group, "\n")), split(read(filename, String), "\n\n"))
end

function parserules(value)
    matches = map(x -> match(r"^([\w\s]+): (\d+)-(\d+) or (\d+)-(\d+)$", x), value)
    return Dict(map(m -> (m.captures[1], (parse(Int, m.captures[2]):parse(Int, m.captures[3]),parse(Int, m.captures[4]):parse(Int, m.captures[5]))), matches))
end

function tointarray(value)
    map(x -> parse(Int,x), split(value, ","))
end

function inrange(rule, field)
    return field in rule[1] || field in rule[2]
end

function rulesapplies(rules, field)
    return any(rule -> inrange(rule, field), values(rules))
end

function findinvalid(rules, nearbytickets)
    invalidtickets = Set()
    invalidfields = Int[]
    for ticket in values(nearbytickets)
        for field in ticket
            if !rulesapplies(rules, field)
                push!(invalidtickets, ticket)
                push!(invalidfields, field)
            end
        end
    end
    return invalidtickets, invalidfields
end

function rulesbyposition(rules, myticket, validtickets)
    possibilities = Dict()
    for i in 1:length(myticket)
        for (label, rule) in rules
            if all(field -> inrange(rule, field), ticket[i] for ticket in validtickets)
                possibilities[i] = push!(get(possibilities, i, String[]), label)
            end
        end
    end
    return possibilities
end

function rulesorder(possiblerules)
    rulespositions = Dict()
    while !isempty(possiblerules)
        for (i, labels) in possiblerules
            if length(labels) == 1
                label = labels[1]
                rulespositions[label] = i
                pop!(possiblerules, i)
                for (i, labels) in possiblerules
                    setdiff!(labels, [label])
                end
                break
            end
        end
    end
    return rulespositions
end

function findfields(myticket, rulespositions, searchstring)
    return collect(myticket[i] for (label, i) in rulespositions if contains(label, searchstring))
end

input = readinput("d16.in")

rules = parserules(input[1])
myticket = tointarray(input[2][2])
nearbytickets = Dict((index, val) for (index, val) in enumerate(map(x -> tointarray(x), input[3][2:end])))

invalidtickets, invalidfields = findinvalid(rules, nearbytickets)
validtickets = setdiff(values(nearbytickets), invalidtickets)
possiblerules = rulesbyposition(rules, myticket, validtickets)
rulespositions = rulesorder(possiblerules)
departurefields = findfields(myticket, rulespositions, "departure")

println("P1: ", sum(invalidfields))
println("P2: ", prod(departurefields))

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
function ols(
        data::ModelSelection.ModelSelectionData;
        fixedvariables::Union{Nothing, Array}=FIXEDVARIABLES_DEFAULT,
        outsample::Union{Nothing, Int, Array}=OUTSAMPLE_DEFAULT,
        criteria::Array=CRITERIA_DEFAULT,
        ttest::Bool=TTEST_DEFAULT,
        modelavg::Bool=MODELAVG_DEFAULT,
        residualtest::Bool=RESIDUALTEST_DEFAULT,
        orderresults::Bool=ORDERRESULTS_DEFAULT
    )

    return ols!(
        ModelSelection.copy_data(data),
        fixedvariables=fixedvariables,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        modelavg=modelavg,
        residualtest=residualtest,
        orderresults=orderresults
    )
end

function ols!(
    data::ModelSelection.ModelSelectionData;
    fixedvariables::Union{Nothing, Array}=FIXEDVARIABLES_DEFAULT,
    outsample::Union{Nothing, Int, Array}=OUTSAMPLE_DEFAULT,
    criteria::Array=CRITERIA_DEFAULT,
    ttest::Bool=TTEST_DEFAULT,
    modelavg::Bool=MODELAVG_DEFAULT,
    residualtest::Bool=RESIDUALTEST_DEFAULT,
    orderresults::Bool=ORDERRESULTS_DEFAULT
)

    result = create_result(data, fixedvariables, outsample, criteria, ttest, modelavg, residualtest, orderresults)
    execute!(data, result)

    ModelSelection.addresult!(data, result)
    data = addextras(data, result)

    return data
end

function execute!(data::ModelSelection.ModelSelectionData, result::AllSubsetRegressionResult)
    if !data.removemissings
        data = ModelSelection.filter_data_by_empty_values(data)
    end

    expvars_num = size(data.expvars, 1)
    if data.intercept
        expvars_num = expvars_num - 1
    end
    num_operations = 2 ^ expvars_num - 1

    depvar_data = convert(SharedArray, data.depvar_data)
    expvars_data = convert(SharedArray, data.expvars_data)
    result_data = fill!(SharedArray{data.datatype}(num_operations, size(result.datanames, 1)), NaN)
    datanames_index = ModelSelection.create_datanames_index(result.datanames)
    if nprocs() == nworkers()
        for order = 1:num_operations
            execute_row!(order, data.depvar, data.expvars, datanames_index, depvar_data, expvars_data, result_data, data.intercept, data.time, data.datatype, result.outsample, result.criteria, result.ttest, result.residualtest, result.fixedvariables)
        end
    else
        ops_per_worker = div(num_operations, nworkers())
        if (nworkers() > num_operations)
            num_jobs = num_operations
        else
            num_jobs = nworkers()
        end
        jobs = []
        for num_job = 1:num_jobs
            push!(jobs, @spawnat num_job+1 execute_job!(num_job, num_jobs, ops_per_worker, data.depvar, data.expvars, datanames_index, depvar_data, expvars_data, result_data, data.intercept, data.time, data.datatype, result.outsample, result.criteria, result.ttest, result.residualtest, result.fixedvariables))
        end
        for job in jobs
            fetch(job)
        end

        remainder = num_operations - ops_per_worker * num_jobs
        if remainder > 0
            for j = 1:remainder
                order = j + ops_per_worker * num_jobs
                execute_row!(order, data.depvar, data.expvars, datanames_index, depvar_data, expvars_data, result_data, data.intercept, data.time, data.datatype, result.outsample, result.criteria, result.ttest, result.residualtest, result.fixedvariables)
            end
        end
    end

    result.data = Array(result_data)

    if :cp in result.criteria
        result.data[ :, datanames_index[:cp] ] = (result.data[ :, datanames_index[:nobs] ] .- maximum(result.data[ :, datanames_index[:ncoef] ]) .- 2) .* (result.data[ :,datanames_index[:rmse] ] ./ minimum(result.data[ :,datanames_index[:rmse] ])) .- (result.data[ :, datanames_index[:nobs] ] .- 2 .* result.data[ :, datanames_index[:ncoef] ])
    end

    len_criteria = length(result.criteria)
    for criteria in result.criteria
        result.data[ :,datanames_index[:order] ] += AVAILABLE_CRITERIA[criteria]["index"] * (1 / len_criteria) * ( (result.data[ :, datanames_index[criteria] ] .- mean(result.data[ :, datanames_index[criteria] ]) ) ./ std(result.data[ :, datanames_index[criteria] ]) )
    end

    if result.modelavg
        delta = maximum(result.data[:, datanames_index[:order]]) .- result.data[:, datanames_index[:order]]
        w1 = exp.(-delta/2)
        result.data[:, datanames_index[:weight]] = w1./sum(w1)
        result.modelavg_data = Vector{Float64}(undef, size(result.datanames))
        weight_pos = (result.ttest) ? 4 : 2
        for expvar in data.expvars
            obs = result.data[:, datanames_index[Symbol(string(expvar, "_b"))]]
            if result.ttest
                obs = hcat(obs, result.data[:, datanames_index[Symbol(string(expvar, "_bstd"))]])
                obs = hcat(obs, result.data[:, datanames_index[Symbol(string(expvar, "_t"))]])
            end
            obs = hcat(obs, result.data[:, datanames_index[:weight]])

            obs = obs[findall(x -> !isnan(obs[x,1]), 1:size(obs,1)),:]
            obs[:, weight_pos] /= sum(obs[:, weight_pos])

            result.modelavg_data[datanames_index[Symbol(string(expvar, "_b"))]] = sum(obs[:, 1] .* obs[:, weight_pos])
            if result.ttest
                result.modelavg_data[datanames_index[Symbol(string(expvar, "_bstd"))]] = sum(obs[:, 2] .* obs[:, weight_pos])
                result.modelavg_data[datanames_index[Symbol(string(expvar, "_t"))]] = sum(obs[:, 3] .* obs[:, weight_pos])
            end
        end

        for criteria in [:nobs, :r2adj, :F, :order]
            if criteria in keys(datanames_index)
                result.modelavg_data[datanames_index[criteria]] = sum(result.data[:, datanames_index[criteria]] .* result.data[:, datanames_index[:weight]])
            end
        end
    end

    if result.orderresults
        result.data = sortrows(result.data, [datanames_index[:order]]; rev=true)
        result.bestresult_data = result.data[1, :]
    else
        max_order = result.data[1, datanames_index[:order]]
        best_result_index = 1
        for i = 1:num_operations
            if result.data[i, datanames_index[:order]] > max_order
                max_order = result.data[i, datanames_index[:order]]
                best_result_index = i
            end
        end
        result.bestresult_data = result.data[best_result_index, :]
    end

    return result
end

function execute_job!(
    num_job,
    num_jobs,
    ops_per_worker,
    depvar,
    expvars,
    datanames_index,
    depvar_data,
    expvars_data,
    result_data,
    intercept,
    time,
    datatype,
    outsample,
    criteria,
    ttest,
    residualtest,
    fixedvariables
)
    for j = 1:ops_per_worker
        order = (j-1) * num_jobs + num_job
        execute_row!(
            order,
            depvar,
            expvars,
            datanames_index,
            depvar_data,
            expvars_data,
            result_data,
            intercept,
            time,
            datatype,
            outsample,
            criteria,
            ttest,
            residualtest,
            fixedvariables,
            num_jobs=num_jobs,
            num_job=num_job,
            iteration_num=j
        )
    end
end

function execute_row!(
    order,
    depvar,
    expvars,
    datanames_index,
    depvar_data,
    expvars_data,
    result_data,
    intercept,
    time,
    datatype,
    outsample,
    criteria,
    ttest,
    residualtest,
    fixedvariables;
    num_jobs=nothing,
    num_job=nothing,
    iteration_num=nothing
)

    selected_variables_index = ModelSelection.get_selected_variables(order, expvars, intercept, fixedvariables=fixedvariables, num_jobs=num_jobs, num_job=num_job, iteration_num=iteration_num)
    depvar_subset, expvars_subset = get_insample_subset(depvar_data, expvars_data, outsample, selected_variables_index)
    outsample_enabled = size(depvar_subset, 1) < size(depvar_data, 1)

    nobs = size(depvar_subset, 1)
    ncoef = size(expvars_subset, 2)
    qrf = qr(expvars_subset)
    b = qrf \ depvar_subset               # estimate
    ŷ = expvars_subset * b                # predicted values
    er = depvar_subset - ŷ                # in-sample residuals
    er2 = er .^ 2                         # squared errors
    sse = sum(er2)                        # residual sum of squares
    df_e = nobs - ncoef                   # degrees of freedom
    rmse = sqrt(sse / nobs)               # root mean squared error
    r2 = 1 - var(er) / var(depvar_subset) # model R-squared

    if ttest
        bstd = sqrt.( sum( (UpperTriangular(qrf.R) \ Matrix(1.0LinearAlgebra.I, ncoef, ncoef) ) .^ 2, dims=2) * (sse / df_e) )
    end

    if outsample_enabled > 0
        depvar_outsample_subset, expvars_outsample_subset = get_outsample_subset(depvar_data, expvars_data, outsample, selected_variables_index)
        erout = depvar_outsample_subset - expvars_outsample_subset * b  # out-of-sample residuals
        sseout = sum(erout .^ 2)                                        # residual sum of squares
        outsample_count = outsample
        if (isa(outsample, Array))
            outsample_count = size(outsample, 1)
        end
        rmseout = sqrt(sseout / outsample_count)                        # root mean squared error
        result_data[ order, datanames_index[:rmseout] ] = rmseout
    end

    result_data[ order, datanames_index[:index] ] = order
    for (index, selected_variable_index) in enumerate(selected_variables_index)
        result_data[ order, datanames_index[Symbol(string(expvars[selected_variable_index], "_b"))] ] = datatype(b[index])
        if ttest
            result_data[ order, datanames_index[Symbol(string(expvars[selected_variable_index], "_bstd"))] ] = datatype(bstd[index])
            result_data[ order, datanames_index[Symbol(string(expvars[selected_variable_index], "_t"))] ] = result_data[order, datanames_index[Symbol(string(expvars[selected_variable_index], "_b"))]] / result_data[order, datanames_index[Symbol(string(expvars[selected_variable_index], "_bstd"))]]
        end
    end

    result_data[ order, datanames_index[:nobs] ] = nobs
    result_data[ order, datanames_index[:ncoef] ] = ncoef
    result_data[ order, datanames_index[:sse] ] = datatype(sse)
    result_data[ order, datanames_index[:r2 ] ] = datatype(r2)
    result_data[ order, datanames_index[:rmse] ] = datatype(rmse)
    result_data[ order, datanames_index[:order] ] = 0

    if :aic in criteria || :aicc in criteria
        aic = 2 * result_data[ order, datanames_index[:ncoef] ] + result_data[ order, datanames_index[:nobs] ] * log(result_data[ order, datanames_index[:sse] ] / result_data[ order, datanames_index[:nobs] ])
    end

    if :aic in criteria
        result_data[ order, datanames_index[:aic] ] = aic
    end

    if :aicc in criteria
        result_data[ order, datanames_index[:aicc] ] = aic + (2(result_data[ order, datanames_index[:ncoef] ] + 1) * (result_data[ order, datanames_index[:ncoef] ] + 2)) / (result_data[ order, datanames_index[:nobs] ] - (result_data[ order, datanames_index[:ncoef] ] + 1 ) - 1)
    end

    if :bic in criteria
        result_data[ order, datanames_index[:bic] ] = result_data[ order, datanames_index[:nobs] ] * log.(result_data[ order, datanames_index[:rmse] ]) + ( result_data[ order, datanames_index[:ncoef] ] - 1 ) * log.(result_data[ order, datanames_index[:nobs] ]) + result_data[ order, datanames_index[:nobs] ] + result_data[ order, datanames_index[:nobs] ] * log(2π)
    end

    if :r2adj in criteria || :r2adj in keys(datanames_index)
        result_data[ order, datanames_index[:r2adj] ] = 1 - (1 - result_data[ order, datanames_index[:r2] ]) * ((result_data[ order, datanames_index[:nobs] ] - 1) / (result_data[ order, datanames_index[:nobs] ] - result_data[ order, datanames_index[:ncoef] ]))
    end

    result_data[ order, datanames_index[:F] ] = (result_data[ order, datanames_index[:r2] ] / (result_data[ order, datanames_index[:ncoef] ] - 1)) / ((1 - result_data[ order, datanames_index[:r2] ]) / (result_data[ order, datanames_index[:nobs] ] - result_data[ order, datanames_index[:ncoef] ]))

    if residualtest
        x = er
        n = length(x)
        m1 = sum(x) / n
        m2 = sum((x .- m1) .^ 2) / n
        m3 = sum((x .- m1) .^ 3) / n
        m4 = sum((x .- m1) .^ 4) / n
        b1 = (m3 / m2 ^ (3 / 2)) ^ 2
        b2 = (m4 / m2 ^ 2)
        statistic = n * b1 / 6 + n * (b2 - 3) ^ 2 / 24
        d = Chisq(2.)
        jbtest = 1 .- cdf(d, statistic)

        regmatw = hcat((ŷ .^ 2), ŷ , ones(size(ŷ, 1)))
        qrfw = qr(regmatw)
        regcoeffw = qrfw \ er2
        residw = er2 - regmatw * regcoeffw
        rsqw = 1 - dot(residw, residw) / dot(er2, er2) # uncentered R^2
        statisticw = n * rsqw
        wtest = ccdf(Chisq(2), statisticw)

        result_data[ order, datanames_index[:wtest] ] = wtest
        result_data[ order, datanames_index[:jbtest] ] = jbtest
        if time != nothing
            e = er
            lag = 1
            xmat = expvars_subset

            n = size(e,1)
            elag = zeros(Float64,n,lag)
            for ii = 1:lag
                elag[ii+1:end, ii] = e[1:end-ii]
            end

            offset = lag
            regmatbg = [xmat[offset+1:end,:] elag[offset+1:end,:]]
            qrfbg = qr(regmatbg)
            regcoeffbg = qrfbg \ e[offset+1:end]
            residbg = e[offset+1:end] .- regmatbg * regcoeffbg
            rsqbg = 1 - dot(residbg,residbg) / dot(e[offset+1:end], e[offset+1:end]) # uncentered R^2
            statisticbg = (n - offset) * rsqbg
            bgtest = ccdf(Chisq(lag), statisticbg)
            result_data[ order, datanames_index[:bgtest] ] = bgtest
        end
    end

end


text:
# Mask operator functions ##################################

struct opMaskFactored{T} <: AbstractLinearTransform{T} end
"""
    op = opMaskFactored(M::SparseMatrixCSC)

Create a masking linear operator that applies to factored matrices.

- forward mode: `op*(U,V)` returns a vector of the elements of the matrix `U*V'` indexed by the mask;

- adjoint mode: `op'*y` returns a sparse matrix with the elements of `y` scattered into the mask pattern.
"""
function opMaskFactored{T}(M::SparseMatrixCSC) where {T}

    # forward operation is a gather.
    function forward(nnzM, I, J, U, V)
        z = zeros(nnzM)
        for (u,v) in zip(eachcol(U), eachcol(V))
            z += [ u[i]*v[j] for (i,j) in zip(I, J) ]
        end
        return z
    end

    # Adjoint operation is scatter.
    adjoint(I, J, w) = sparse(I, J, w)

    dimsDomain = size(M)
    dimsRange = nnz(M)
    I, J = findnz(M)[1:2]
    prod = uv -> forward(nnz(M), I, J, uv[1], uv[2])
    tprod = w -> adjoint(I, J, w)
    LinearTransform{T}(dimsRange, dimsDomain, prod, tprod)
end
opMaskFactored(M::SparseMatrixCSC) = opMaskFactored{Float64}(M::SparseMatrixCSC)

struct opMask{T} <: AbstractLinearTransform{T} end
"""
    op = opMask(M::SparseMatrixCSC)

Create a masking linear operator.

- forward mode: `op*(X)` returns a vector of the elements of the matrix `X` indexed by the mask;

- adjoint mode: `op'*y` returns a sparse matrix with the elements of `y` scattered into the mask pattern.
"""
function opMask{T}(M::SparseMatrixCSC) where {T}
    dimsDomain = size(M)
    dimsRange = nnz(M)
    I, J = findnz(M)[1:2]
    prod = X -> [X[i,j] for (i,j) in zip(I,J)]
    tprod = w -> sparse(I, J, w, dimsDomain...)
    LinearTransform{T}(dimsRange, dimsDomain, prod, tprod)
end
opMask(M::SparseMatrixCSC) = opMask{Float64}(M)


text:
using RevisedAffineArithmetic
using Documenter

DocMeta.setdocmeta!(
    RevisedAffineArithmetic, :DocTestSetup, :(using RevisedAffineArithmetic); recursive=true
)

makedocs(;
    modules=[RevisedAffineArithmetic],
    authors="Luca Ferranti",
    repo="https://github.com/lucaferranti/RevisedAffineArithmetic.jl/blob/{commit}{path}#{line}",
    sitename="RevisedAffineArithmetic.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lucaferranti.github.io/RevisedAffineArithmetic.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/lucaferranti/RevisedAffineArithmetic.jl", devbranch="main")


text:
"""
    log_inverse_gamma(x; α=1., β=1., d=0)

Log of the InverseGamma PDF.
Equivalent to `using Distributions; logpdf(InverseGamma(α, β), x)`
https://en.wikipedia.org/wiki/Inverse-gamma_distribution

# Keyword Arguments
- `x::Real`: Input
- `α::Real=1.`: Shape parameter
- `β::Real=1.`: Scale parameter
- `d::Integer=0`: How many derivatives to take
"""
function log_inverse_gamma(x::Real; α::Real=1., β::Real=1., d::Integer=0)
    @assert 0 <= d <= 2 "Can only differentiate up to two times"
    @assert α > 0 "α needs to be > 0"
    @assert β > 0 "β needs to be > 0"
    if d == 0
        x > 0 ? val = -(β / x) - (1 + α) * log(x) + α * log(β) - loggamma(α) : val = -Inf
    elseif d == 1
        x > 0 ? val = (β / x - (1 + α)) / x : val = 0
    else
        x > 0 ? val = (-2 * β / x + (1 + α)) / (x * x) : val = 0
    end
    return val
end

"""
    log_gamma(x, parameters; d=0, passed_mode_std=false)

Log of the Gamma PDF.
Equivalent to `using Distributions; logpdf(Gamma(α, θ), x)`
https://en.wikipedia.org/wiki/Gamma_distribution

# Arguments
- `x::Real`: Input
- `parameters::Vector`: Either the shape and scale parameters (i.e. [α, θ]) or a mode and standard deviation
- `d::Integer=0`: How many derivatives to take
- `passed_mode_std::Bool=false`: Whether to the `parameters` need to be converted to shape and scale parameters
"""
function log_gamma(x::Real, parameters::Vector{<:Real}; d::Integer=0, passed_mode_std::Bool=false)
    @assert 0 <= d <= 2 "Can only differentiate up to two times"
    @assert length(parameters) == 2
    assert_positive(parameters)
    if passed_mode_std
        parameters = gamma_mode_std_to_α_θ(parameters[1], parameters[2])
    end
    α = parameters[1]
    θ = parameters[2]
    if d == 0
        x > 0 ? val = -(x / θ) + (α - 1) * log(x) - α * log(θ) - loggamma(α) : val = -Inf
    elseif d == 1
        x > 0 ? val = (α - 1) / x - 1 / θ : val = 0
    else
        x > 0 ? val = -(α - 1) / (x * x) : val = 0
    end
    return val
end

"""
    gamma_mode_std_to_α_θ(m, σ)

Convert a gamma distribution's mode and standard deviation into it's shape and
scale parameters

# Arguments
- `m::Real`: Desired mode of gamma distribution
- `σ::Real`: Desired standard deviation of gamma distribution
"""
function gamma_mode_std_to_α_θ(m::Real, σ::Real)
    θ = (sqrt(m ^ 2 + 4 * σ ^ 2) - m) / 2
    α = m / θ + 1
    return [α, θ]
end

"""
    gauss_cdf(x)

The CDF of a Gaussian (i.e. (1 + erf(x))/2)
"""
gauss_cdf(x::Real) = (1 + erf(x))/2

"""
    log_gaussian(x, parameters; d=0, min=-Inf, max=Inf)

Log of the Gaussian PDF.
Equivalent to `using Distributions; logpdf(Gaussian(μ, σ), x)`

# Arguments
- `x::Real`: Input
- `parameters::Vector`: The mean and standard deviation parameters (i.e. [μ, σ])
- `d::Integer=0`: How many derivatives to take
- `min::Real=-Inf`: Where to minimally truncate the Gaussian
- `max::Real=Inf`: Where to maximally truncate the Gaussian
"""
function log_gaussian(x::Real, parameters::Vector{<:Real}; d::Integer=0, min::Real=-Inf, max::Real=Inf)
    @assert 0 <= d <= 2 "Can only differentiate up to two times"
    @assert length(parameters) == 2
    μ = parameters[1]
    σ = parameters[2]
    assert_positive(σ)
    @assert min < max
    if d == 0
        (min < x < max) ? val = (-((x - μ) / σ) ^ 2 - log(2 * π * σ * σ)) / 2 - log(1 - gauss_cdf(min - μ) - gauss_cdf(μ - max)) : val = -Inf
    elseif d == 1
        (min < x < max) ? val = -(x - μ)/(σ * σ) : val = 0
    else
        (min < x < max) ? val = -1 / (σ * σ) : val = 0
    end
    return val
end

"""
    log_uniform(x; d=0, min=0, max=1)

Log of the Uniform PDF.

# Arguments
- `x::Real`: Input
- `min_max::Vector=[0,1]`: Where to truncate the Uniform
- `d::Integer=0`: How many derivatives to take
"""
function log_uniform(x::Real; min_max::Vector{<:Real}=[0,1], d::Integer=0)
    @assert 0 <= d <= 2 "Can only differentiate up to two times"
    @assert length(min_max) == 2
    min, max = min_max
    @assert min < max
    if d == 0
        min <= x <= max ? -log(max - min) : -Inf
    else
        return 0
    end
end

"""
    log_loguniform(x::Real, min_max::Vector{<:Real}; d::Integer=0, shift::Real=0)

Log of the log-Uniform PDF. Flattens out in log space starting at shift
Also known as a (modified in shifted case) Jeffrey's prior

# Arguments
- `x::Real`: Input
- `min_max::Vector`: Where to truncate the log-Uniform
- `d::Integer=0`: How many derivatives to take
- `shift::Real=0`: Where to shift the peak of the distribution
"""
function log_loguniform(x::Real, min_max::Vector{<:Real}; d::Integer=0, shift::Real=0)
    @assert 0 <= d <= 2 "Can only differentiate up to two times"
    @assert length(min_max) == 2
    min, max = min_max
    @assert 0 < min + shift < max + shift
    xpshift = x + shift
    if d == 0
        min <= x <= max ? val = -log(xpshift) - log(log((max + shift)/(min + shift))) : val = -Inf
    elseif d == 1
        min <= x <= max ? val = -1 / xpshift : val = 0
    elseif d == 2
        min <= x <= max ? val = 1 / (xpshift * xpshift) : val = 0
    end
    return val
end

"""
    log_Rayleigh(x, σ; d=0, cutoff=Inf)

Log of the Rayleigh PDF.
https://en.wikipedia.org/wiki/Rayleigh_distribution

# Arguments
- `x::Real`: Input
- `σ::Real`: The mode
- `d::Integer=0`: How many derivatives to take
- `cutoff::Real=Inf`: Where to cutoff the tail of the distribution
"""
function log_Rayleigh(x::Real, σ::Real; d::Integer=0, cutoff::Real=Inf)
    @assert 0 <= d <= 2
    @assert cutoff > 0
    cutoff == Inf ? normalization = 0 : normalization = -log(1 - exp(-cutoff * cutoff / 2 / σ / σ))
    if d == 0
        0 <= x <= cutoff ? val = normalization - (x * x / 2 / σ / σ) + log(x / σ / σ) : val = -Inf
    elseif d == 1
        0 <= x <= cutoff ? val = 1 / x - x / σ / σ : val = 0
    elseif d == 2
        0 <= x <= cutoff ? val = -1 / x / x - 1 / σ / σ : val = 0
    end
    return val
end

"""
    log_circle(xs, min_max_r; d=[0,0])

Log of the 2D circle PDF.

# Arguments
- `xs::Real`: Inputs
- `min_max_r::Vector{<:Real}`: How far the inner and outer edge of the circle extends from the origin
- `d::Vector{<:Integer}=[0,0]`: How many derivatives to take
"""
function log_circle(xs::Vector{<:Real}, min_max_r::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) == 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(xs) == length(min_max_r) == length(d) == 2
    min_r, max_r = min_max_r

    # @assert min_r < sqrt(dot(x, x)) < max_r
    # @assert min_θ < atan(x[1], x[2]) < max_θ
    if all(d .== 0)
        min_r < sqrt(dot(xs, xs)) < max_r ? val = -log(2 * π * (max_r ^ 2 - min_r ^ 2)) : val = -Inf
    else
        return 0
    end
    return val
end

"""
    log_cone(xs; d=[0,0])

Log of the 2D unit cone PDF.

# Arguments
- `xs::Real`: Inputs
- `d::Vector{<:Integer}=[0,0]`: How many derivatives to take
"""
function log_cone(xs::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(xs) == length(d) == 2

    r_sq = dot(xs, xs)  # x^2 + y^2
    r = sqrt(r_sq)
    if d == [0,0]
        0 <= r < 1 ? val = log(3 / π * (1 - r)) : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = xs[2] / (r_sq - r) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = (-xs[2] ^ 2 * r + xs[1] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = xs[1] / (r_sq - r) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = (xs[1] * xs[2] * (1 - 2 * r)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = (-xs[1] ^ 2 * r + xs[2] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    end
    return val
end

"""
    log_quad_cone(xs; d=[0,0])

Log of the 2D unit quadratic cone PDF.

# Arguments
- `xs::Real`: Inputs
- `d::Vector{<:Integer}=[0,0]`: How many derivatives to take
"""
function log_quad_cone(xs::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(xs) == length(d) == 2

    r_sq = dot(xs, xs)  # x^2 + y^2
    r = sqrt(r_sq)
    if d == [0,0]
        0 <= r < 1 ? val = log(6 / π * (1 - 2 * r + r_sq)) : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = 2 * xs[2] / (r_sq - r) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = 2 * (-xs[2] ^ 2 * r + xs[1] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = 2 * xs[1] / (r_sq - r) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = (2 * xs[1] * xs[2] * (1 - 2 * r)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = 2 * (-xs[1] ^ 2 * r + xs[2] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    end
    return val
end

"""
    log_cubic_cone(xs; d=[0,0])

Log of the 2D unit cubic cone PDF.

# Arguments
- `xs::Real`: Inputs
- `d::Vector{<:Integer}=[0,0]`: How many derivatives to take
"""
function log_cubic_cone(xs::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(xs) == length(d) == 2

    r_sq = dot(xs, xs)  # x^2 + y^2
    r = sqrt(r_sq)
    if d == [0,0]
        0 <= r < 1 ? val = log(10 / π * (1 - r)^3) : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = 3 * xs[2] / (r_sq - r) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = 3 * (-xs[2] ^ 2 * r + xs[1] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = 3 * xs[1] / (r_sq - r) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = (3 * xs[1] * xs[2] * (1 - 2 * r)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = 3 * (-xs[1] ^ 2 * r + xs[2] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    end
    return val
end

"""
    log_rot_Rayleigh(xs; d=[0,0], σ=1/5, cutoff=Inf)

Log of the 2D rotated Rayleigh PDF

# Arguments
- `xs::Real`: Inputs
- `d::Vector{<:Integer}=[0,0]`: How many derivatives to take
- `σ::Real=1/5`: the radial mode of the distribution
- `cutoff::Real=Inf`: Where to cutoff the tail of the distribution
"""
function log_rot_Rayleigh(xs::Vector{<:Real}, σ::Real; d::Vector{<:Integer}=[0,0], cutoff::Real=Inf)
    @assert all(0 .<= d .<= 2)
    @assert sum(d) <= 2
    @assert cutoff > 0

    @assert length(xs) == length(d) == 2
    r_sq = dot(xs, xs)  # x^2 + y^2
    r = sqrt(r_sq)
    σ_sq = σ ^ 2
    cutoff == Inf ? normalization = -log(2 * σ * σ * π * π * π) / 2 : normalization = -log(2 * π) - log(sqrt(π / 2) * σ * erf(cutoff / sqrt(2) / σ) - cutoff * exp(- cutoff * cutoff / 2 / σ / σ))
    if d == [0,0]
        0 <= r < cutoff ? val = normalization - r_sq / (2 * σ_sq) + log(r) - (2 * log(σ)) : val = -Inf
    elseif d == [0,1]
        0 <= r < cutoff ? val = -xs[2] * (r_sq - σ_sq) / (r_sq * σ_sq) : val = 0
    elseif d == [0,2]
        0 <= r < cutoff ? val = -(xs[1] ^ 4 + xs[1] ^ 2 * (2 * xs[2] ^ 2 - σ_sq) + xs[2] ^ 2 * (xs[2] ^ 2 + σ_sq)) / (r_sq ^ 2 * σ_sq) : val = 0
    elseif d == [1,0]
        0 <= r < cutoff ? val = -xs[1] * (r_sq - σ_sq) / (r_sq * σ_sq) : val = 0
    elseif d == [1,1]
        0 <= r < cutoff ? val = -2 * xs[1] * xs[2] / r_sq ^ 2 : val = 0
    elseif d == [2,0]
        0 <= r < cutoff ? val = -(xs[2] ^ 4 + xs[2] ^ 2 * (2 * xs[1] ^ 2 - σ_sq) + xs[1] ^ 2 * (xs[1] ^ 2 + σ_sq)) / (r_sq ^ 2 * σ_sq) : val = 0
    end
    return val
end

"""
    log_bvnormal(xs, Σ; μ=zeros(T, 2), d=[0,0], lows=zeros(T, 2) .- Inf)

Log of the bivariate normal PDF
NOTE THAT THAT WHEN USING lows!=[-∞,...], THIS IS NOT PROPERLY NORMALIZED

# Arguments
- `xs::Real`: Inputs
- `Σ::Colesky`: The covariance matrix of the distribution
- `μ::Vector=zeros(T, 2)`: The mean of the distribution
- `d::Vector{<:Integer}=[0,0]`: How many derivatives to take
- `lows::Vector=zeros(T, 2) .- Inf`: The lower cutoffs of the distribution
"""
function log_bvnormal(xs::Vector{T}, Σ::Cholesky; μ::Vector{T}=zeros(T, 2), d::Vector{<:Integer}=[0,0], lows::Vector{T}=zeros(T, 2) .- Inf) where {T<:Real}
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(xs) == length(d) == length(μ) == 2
    y = xs - μ

    if sum(d) == 0
        all(xs .>= lows) ? val = -nlogL(Σ, y) : val = -Inf
    elseif sum(d) == 1
        y1 = Float64.(d)
        all(xs .>= lows) ? val = -dnlogLdθ(y1, Σ \ y) : val = 0
    elseif d == [0,2]
        y1 = y2 = [0, 1.]
        all(xs .>= lows) ? val = -d2nlogLdθ(y2, [0,0.], Σ \ y, Σ \ y1) : val = 0
    elseif d == [1,1]
        y1 = [1, 0.]
        y2 = [0, 1.]
        all(xs .>= lows) ? val = -d2nlogLdθ(y2, [0,0.], Σ \ y, Σ \ y1) : val = 0
    elseif d == [2,0]
        y1 = y2 = [1, 0.]
        all(xs .>= lows) ? val = -d2nlogLdθ(y2, [0,0.], Σ \ y, Σ \ y1) : val = 0
    end
    return val
end

"""
    bvnormal_covariance(σ11, σ22, ρ)

Converting 2 standard deviations and their correlation into a bivariate
covariance matrix and Cholesky factorizing the result.
"""
function bvnormal_covariance(σ11::Real, σ22::Real, ρ::Real)
    assert_positive(σ11, σ22)
    @assert 0 <= ρ <= 1
    v12 = σ11*σ22*ρ
    return ridge_chol([σ11^2 v12;v12 σ22^2])
end

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
################################################################################
#
#  Infinite places for number fields
#
#
#  (C) 2017 Tommy Hofmann
#
################################################################################

export iscomplex, ispositive, istotally_positive, signs, sign, real_places,
       complex_places, infinite_places, infinite_place

################################################################################
#
#  Equality and hashing
#
################################################################################

function Base.:(==)(P::InfPlc, Q::InfPlc)
  P.K != Q.K && error("Places are of different number fields")
  return P.i == Q.i
end

function Base.hash(P::InfPlc, h::UInt)
  return xor(Base.hash(P.K, h), Base.hash(P.i, h))
end

################################################################################
#
#  Printing
#
################################################################################

function Base.show(io::IO, P::InfPlc)
  if P.isreal
    print(io, "Real ")
  else
    print(io, "Complex ")
  end

  print(io, "place of \n$(P.K)\n")
  print(io, "corresponding to the root $(P.r)")
end

################################################################################
#
#  Predicates
#
################################################################################

doc"""
***
    isreal(P::InfPlc) -> Bool

> Returns whether the embedding into $\mathbf{C}$ defined by P is real or not.
"""
function Base.isreal(P::InfPlc)
  return P.isreal
end

doc"""
***
    iscomplex(P::InfPlc) -> Bool

> Returns whether the embedding into $\mathbf{C}$ defined by P is complex or not.
"""
function iscomplex(P::InfPlc)
  return !isreal(P)
end

################################################################################
#
#  Construction
#
################################################################################

doc"""
***
    infinite_place(K::AnticNumberField, i::Int) -> InfPlc

> This function returns the infinite place of $K$ corresponding to the root
> `conjugates_arb(a)[i]`, where `a` is the primitive element of $K$.
"""
function infinite_place(K::AnticNumberField, i::Int)
  !(1 <= i <= degree(K)) && error("Index must be between 1 and $(degree(K))")
  return InfPlc(K, i)
end

doc"""
***
    infinite_places(K::AnticNumberField) -> Vector{InfPlc}

> This function returns all infinite places of $K$.
"""
function infinite_places(K::AnticNumberField)
  r1, r2 = signature(K)
  return [ InfPlc(K, i) for i in 1:(r1 + r2)]
end

doc"""
***
    real_places(K::AnticNumberField) -> Vector{InfPlc}

> This function returns all infinite real places of $K$.
"""
function real_places(K::AnticNumberField)
  r1, r2 = signature(K)
  return [ InfPlc(K, i) for i in 1:r1]
end

doc"""
***
    complex_places(K::AnticNumberField) -> Vector{InfPlc}

> This function returns all infinite complex places of $K$.
"""
function complex_places(K::AnticNumberField)
  r1, r2 = signature(K)
  return [ InfPlc(K, i) for i in (r1 + 1):r1 + r2]
end

################################################################################
#
#  Signs at places
#
################################################################################

doc"""
***
    signs(a::nf_elem)          -> Dict{InfPlc, Int}
    signs(a::FacElem{nf_elem}) -> Dict{InfPlc, Int}

> This function returns a dictionary of the signs of $a$ at all infinite places
> of the ambient number field. The keys are infinite places of the ambient
> number field. The value is $1$ if the sign is positive and $-1$ if the sign
> is negative.
"""
function signs(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}})
  K = _base_ring(a)
  r1, r2 = signature(K)
  D = Dict{InfPlc, Int}()
  s = _signs(a)
  for i in 1:r1
    P = InfPlc(K, i)
    D[P] = s[i]
  end
  return D
end

doc"""
***
    sign(a::nf_elem, P::InfPlc)          -> Int
    sign(a::FacElem{nf_elem}, P::InfPlc) -> Int

> This function returns the sign of $a$ at the place $P$. The value is $1$ if
> the sign is positive and $-1$ if the sign is negative.

"""
function sign(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, P::InfPlc)
  !isreal(P) && error("Place must be real")
  return signs(a, [P])[P]
end

doc"""
***
    signs(a::nf_elem, l::Vector{InfPlc})          -> Dict{InfPlc, Int}
    signs(a::FacElem{nf_elem}, l::Vector{InfPlc}) -> Dict{InfPlc, Int}

> This function returns a dictionary of the signs of $a$ at places in $l$. The
> keys are the elements of $l$. The value is $1$ if the sign is positive and
> $-1$ if the sign is negative. The result will contain as many signs as there
> are real places contained in $l$.
"""
function signs(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, l::Array{InfPlc, 1})
  K = _base_ring(a)
  r1, r2 = signature(K)
  D = Dict{InfPlc, Int}()
  s = _signs(a)
  for i in 1:r1
    P = InfPlc(K, i)
    if P in l
      D[P] = s[i]
    end
  end
  return D
end

# extend functionality to elements of orders

function signs(a::NfOrdElem, args...)
  return signs(a.elem_in_nf, args...)
end

function sign(a::NfOrdElem, args...)
  return sign(a.elem_in_nf, args...)
end

################################################################################
#
#  Positivity
#
################################################################################

doc"""
***
    ispositive(a::nf_elem, P::InfPlc)          -> Bool
    ispositive(a::FacElem{nf_elem}, P::InfPlc) -> Bool

> Returns whether the element $a$ is positive at the embedding corresponding to
> $P$. The place $P$ must be real.
"""
function ispositive(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, P::InfPlc)
  !isreal(P) && error("Place must be real")
  return sign(a, P) > 0
end

doc"""
***
    ispositive(a::nf_elem, l::Vector{InfPlc})          -> Bool
    ispositive(a::FacElem{nf_elem}, l::Vector{InfPlc}) -> Bool

> Returns whether the element $a$ is positive at the embeddings corresponding to
> the real places of $l$.
"""
function ispositive(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, l::Array{InfPlc, 1})
  return all(x -> ispositive(a, x), (y for y in l if isreal(y)))
end

doc"""
***
    istotally_positive(a::nf_elem)          -> Bool
    istotally_positive(a::FacElem{nf_elem}) -> Bool

> Returns whether the element $a$ is totally positive, that is, whether it is
> positive at all places of the ambient number field.
"""
function istotally_positive(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}})
  K = _base_ring(a)
  return ispositive(a, real_places(K))
end

# extend functionality to elements of orders

function ispositive(a::NfOrdElem, args...)
  return ispositive(a.elem_in_nf, args...)
end

function istotally_positive(a::NfOrdElem, args...)
  return istotally_positive(a.elem_in_nf, args...)
end


text:
export getData

function getData(m,pFor::FWIparam,doClear::Bool=false)

    # extract pointers
    M       	= pFor.Mesh
    omega   	= pFor.omega
	wavelet 	= pFor.WaveletCoef;
    gamma   	= pFor.gamma
    Q       	= pFor.Sources
    P       	= pFor.Receivers
	Ainv    	= pFor.ForwardSolver;
	batchSize 	= pFor.forwardSolveBatchSize;
	select  	= pFor.sourceSelection;

    nrec  		= size(P,2)
    nsrc  		= size(Q,2)

	# allocate space for data and fields
	n = prod(M.n+1);
	# ALL AT ONCE DIRECT CODE
	H = GetHelmholtzOperator(M,m,omega, gamma, true,useSommerfeldBC);

	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		Ainv.helmParam = HelmholtzParam(M,gamma,m,omega,true,useSommerfeldBC);
		H = H + GetHelmholtzShiftOP(m, omega,Ainv.shift[1]);
		H = H';
		# H is actually shifted laplacian now...
		Ainv.MG.relativeTol *= 1e-4;
	end

	if select==[]
		Qs = Q*wavelet;
	else
		Qs = Q[:,select]*wavelet;
	end

	nsrc 		= size(Qs,2);

	if batchSize > nsrc
		batchSize = nsrc;
	end

	Fields = [];

	if doClear==false
		if pFor.useFilesForFields
			tfilename = getFieldsFileName(omega);
			tfile     = matopen(tfilename, "w");
		else
			Fields    = zeros(FieldsType,n   ,nsrc);
		end
	end

	numBatches 	= ceil(Int64,nsrc/batchSize);
	D 			= zeros(Complex128,nrec,nsrc);
	U 			= zeros(Complex128,n,batchSize);

	Ainv.doClear = 1;
	for k_batch = 1:numBatches
		batchIdxs = (k_batch-1)*batchSize + 1 : min(k_batch*batchSize,nsrc);
		if length(length(batchIdxs))==batchSize
			U[:] = full(Qs[:,batchIdxs]);
		else
			U = convert(Array{Complex128},full(Qs[:,batchIdxs]));
		end
		U,Ainv = solveLinearSystem(H,U,Ainv,0)

		# pic = round(Int64,rand()*1000);
		# u = reshape(real(U[:,1]),M.n[1]+1,M.n[2]+1);
		# figure(pic)
		# imshow(u'); colorbar();

		Ainv.doClear = 0;
		D[:,batchIdxs]      = (P'*U);

		if doClear==false
			if pFor.useFilesForFields
				write(tfile,string("Ubatch_",k_batch),convert(Array{Complex64},U));
			else
				Fields[:,batchIdxs] = U;
			end
		end
	end

	# pic = round(Int64,rand()*1000);
	# figure(pic)
	# imshow(real(D')); colorbar();

	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		Ainv.MG.relativeTol *= 1e+4;
	end

	pFor.ForwardSolver = Ainv;

	if doClear==false
		if pFor.useFilesForFields
			close(tfile);
		else
			pFor.Fields = Fields;
		end
	end

	if doClear
		clear!(pFor);
	elseif isa(Ainv,ShiftedLaplacianMultigridSolver)
		clear!(Ainv.MG);
	end
    return D,pFor
end


text:
@inline identity(i, j, k, grid, c) = @inbounds c[i, j, k]
@inline identity(i, j, k, grid, a::Number) = a

"""Evaluate the function `F` with signature `F(i, j, k, grid, args...)` at index `i, j, k` without
interpolation."""
@inline identity(i, j, k, grid, F::TF, args...) where TF<:Function = F(i, j, k, grid, args...)

# Utilities for inferring the interpolation function needed to
# interpolate a field from one location to the next.
interpolation_code(::Type{Face}) = :ᶠ
interpolation_code(::Type{Cell}) = :ᶜ
interpolation_code(::Face) = :ᶠ
interpolation_code(::Cell) = :ᶜ
interpolation_code(from::L, to::L) where L = :ᵃ
interpolation_code(from, to) = interpolation_code(to)

for ξ in ("x", "y", "z")
    ▶sym = Symbol(:ℑ, ξ, :sym)
    @eval begin
        $▶sym(s::Symbol) = $▶sym(Val(s))
        $▶sym(::Union{Val{:ᶠ}, Val{:ᶜ}}) = $ξ
        $▶sym(::Val{:ᵃ}) = ""
    end
end

instantiate(L::NTuple{N, <:DataType}) where N = Tuple(X() for X in L)
instantiate(x, y...) = instantiate(tuple(x, y...))

"""
    interpolation_operator(from, to)

Returns the function to interpolate a field `from = (XA, YZ, ZA)`, `to = (XB, YB, ZB)`,
where the `XA`s and `XB`s are types `Face` or `Cell` (rather than type instances).
"""
interpolation_operator(from::NTuple{N, <:DataType}, to::NTuple{N, <:DataType}) where N =
    interpolation_operator(instantiate(from), instantiate(to))

"""
    interpolation_operator(from, to)

Returns the function to interpolate a field `from = (XA, YZ, ZA)`, `to = (XB, YB, ZB)`,
where the `XA`s and `XB`s are `Face()` or `Cell()` instances.
"""
function interpolation_operator(from, to)
    x, y, z = (interpolation_code(X, Y) for (X, Y) in zip(from, to))

    if all(ξ === :ᵃ for ξ in (x, y, z))
        return identity
    else
        return eval(Symbol(:ℑ, ℑxsym(x), ℑysym(y), ℑzsym(z), x, y, z))
    end
end

"""
    interpolation_operator(::Nothing, to)

Return the `identity` interpolator function. This is needed to obtain the interpolation
opertator for fields that have no instrinsic location, like numbers or functions.
"""
interpolation_operator(::Nothing, to) = identity

"""
    insert_location(ex::Expr, location)

Insert a symbolic representation of `location` into the arguments of an `expression`.

Used in the `@at` macro for specifying the location of an `AbstractOperation`.
"""
function insert_location!(ex::Expr, location)
    if ex.head === :call && ex.args[1] ∈ operators
        push!(ex.args, ex.args[end])
        ex.args[3:end-1] .= ex.args[2:end-2]
        ex.args[2] = location
    end

    for arg in ex.args
        insert_location!(arg, location)
    end

    return nothing
end

"Fallback for when `insert_location` is called on objects other than expressions."
insert_location!(anything, location) = nothing

"""
    @at location abstract_operation

Modify the `abstract_operation` so that it returns values at
`location`, where `location` is a 3-tuple of `Face`s and `Cell`s.
"""
macro at(location, abstract_operation)
    insert_location!(abstract_operation, location)
    return esc(abstract_operation)
end


text:
using LinearAlgebra

#Givens from Problem Statement
#Shear Modulus(psi),diameter(in),length(in),torque(lb*in)
G = 3.9*10^6
G2 = 4.0*10^6
d = 1.5
d2 = 1
l = 2.5*12
l2 = 1*12
l3 =2*12
τ = -200*12

# Polar Second Area Moment of a Cylinder (in^4)
J_cyl(d) = 0.5*pi*(d/2)^4
J = J_cyl(d)
J2 = J_cyl(d2)

# Torsional Stiffness Function
k(G,J,l) = (G*J)/l

# k_local(Elastic Modulus,thickness,height,length,number of nodes,
# row/column one,row/column two, ==1 left fixed/ ==2 right fixed)
function k_local(G,J,l,i,n1,n2,fixed=0)
    k_global = zeros(i,i)
    kin = k(G,J,l)
    k_global[n1,n1] = kin
    k_global[n1,n2] = -kin
    k_global[n2,n1] = -kin
    k_global[n2,n2] = kin
    if fixed == 1
      adj = zeros(1,i)
      adj[n1] = 1
      k_global[n1,:] = adj
      k_global[:,n1] = adj
    elseif fixed == 2
      adj = zeros(1,i)
      adj[n2] = 1
      k_global[n2,:] = adj
      k_global[:,n2] = adj
    end
    return k_global
end

# Global for angular displacement calculation
k_1 = k_local(G,J,l,4,1,2,1)
k_2 = k_local(G,J,l2,4,2,3)
k_3 = k_local(G2,J2,l3,4,3,4,2)
kG = k_1+k_2+k_3

#Global for Reaction calculation
k_11 = k_local(G,J,l,4,1,2)
k_22 = k_local(G,J,l2,4,2,3)
k_33 = k_local(G2,J2,l3,4,3,4)
kGG = k_11+k_22+k_33

# Torque Vector
t = Float64[0 τ 0 0]
T = transpose(t)

# Angular Displacement Vector
θ = kG\T
println("Nodal Angular Displacement Vector (radians)")
display(θ)

# Moment Reaction Vector
τ_app = kGG*θ
println("Moment Reaction Vector (lb*in):")
display(τ_app)


text:
# features specifically related to Flux.jl
@info "RoMEPlotting is adding Flux related functionality."

@eval PlotTypesPose2 = Union{Type{Pose2Pose2}, Type{Pose2Point2BearingRange}, Type{Pose2Point2Range}, Type{Pose2Point2Bearing}, Type{MixtureFluxPose2Pose2}}

@eval ExtendedPose2Pose2Types = Union{Pose2Pose2, MixtureFluxPose2Pose2}

# rerun since a few functions need to be reloaded with Flux
include("PlotFactorsReload.jl")

#


text:
module GlobalDiffEq

using Reexport
@reexport using DiffEqBase

import OrdinaryDiffEq, Richardson

abstract type GlobalDiffEqAlgorithm <: DiffEqBase.AbstractODEAlgorithm end

struct GlobalRichardson{A} <: GlobalDiffEqAlgorithm
    alg::A
end

function DiffEqBase.__solve(prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDAEProblem},
                            alg::GlobalRichardson, args...;
                            dt,kwargs...)
    opt = Dict(kwargs)
    otheropts = delete!(copy(opt), :dt)
    tstops = get(opt, :tstops, range(prob.tspan[1], stop=prob.tspan[2], step=dt))
    local sol
    val, err = Richardson.extrapolate(dt, rtol=get(opt, :reltol, 1e-3), atol=get(opt, :abstol, 1e-6), contract=0.5) do _dt
        sol = solve(prob, alg.alg, args...; dt=_dt, adaptive=false, otheropts...)
        sol.(tstops)
    end
    return sol
end

export GlobalRichardson

end


text:
using ..ChemistryFeaturization.Utils.ElementFeatureUtils
using DataFrames

# TODO: figure out what scheme would look like that is flexible to direct-value encoding (may just need a different feature type since it'll have to handle normalization, etc. too)
"""
    ElementFeatureDescriptor(feature_name, encode_f, decode_f, categorical, contextual, length, encodable_elements)

Construct a feature object that encodes features associated with individual atoms that depend only upon their elemental identity (if you want to encode a feature that depends upon an atom's environment, you shold use SpeciesFeatureDescriptor!)

## Arguments
- `name::String`: the name of the feature
- `categorical::Bool`: flag for whether the feature is categorical or continuous-valued
- `length::Int`: length of encoded vector
- `logspaced::Bool`: whether onehot-style bins should be logarithmically spaced or not
- `lookup_table::DataFrame`: table containing values of feature for every encodable element
"""
struct ElementFeatureDescriptor <: AbstractAtomFeatureDescriptor
    name::String
    length::Integer
    nbins::Integer
    logspaced::Bool
    categorical::Bool
    lookup_table::DataFrame
end

function ElementFeatureDescriptor(
    feature_name::String,
    lookup_table::DataFrame = atom_data_df;
    nbins::Integer = default_nbins,
    logspaced::Bool = default_log(feature_name, lookup_table),
    categorical::Bool = default_categorical(feature_name, lookup_table),
)
    colnames = names(lookup_table)
    @assert feature_name in colnames && "Symbol" in colnames "Your lookup table must have a column called :Symbol and one with the same name as your feature to be usable!"

    local vector_length
    if categorical
        vector_length = length(unique(skipmissing(lookup_table[:, Symbol(feature_name)])))
    else
        vector_length = nbins
    end

    lookup_table = lookup_table[:, ["Symbol", feature_name]]
    dropmissing!(lookup_table)

    ElementFeatureDescriptor(
        feature_name,
        vector_length,
        nbins,
        logspaced,
        categorical,
        lookup_table,
    )
end

# pretty printing, short version
Base.show(io::IO, af::ElementFeatureDescriptor) = print(io, "ElementFeature $(af.name)")

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", af::ElementFeatureDescriptor)
    st = "ElementFeature $(af.name):\n   categorical: $(af.categorical)\n   encoded length: $(af.length)"
    print(io, st)
end

encodable_elements(f::ElementFeatureDescriptor) = f.lookup_table[:, :Symbol]

function encodable_elements(feature_name::String, lookup_table::DataFrame = atom_data_df)
    info = lookup_table[:, [Symbol(feature_name), :Symbol]]
    return info[
        findall(x -> !ismissing(x), getproperty(info, Symbol(feature_name))),
        :Symbol,
    ]
end

function (f::ElementFeatureDescriptor)(a::AbstractAtoms)
    @assert all([el in encodable_elements(f) for el in a.elements]) "Feature $(f.name) cannot encode some element(s) in this structure!"
    reduce(
        hcat,
        map(
            e -> onehot_lookup_encoder(
                e,
                f.name,
                f.lookup_table;
                nbins = f.nbins,
                logspaced = f.logspaced,
                categorical = f.categorical,
            ),
            a.elements,
        ),
    )
end

decode(f::ElementFeatureDescriptor, encoded_feature) = onecold_decoder(
    encoded_feature,
    f.name,
    f.lookup_table;
    nbins = f.nbins,
    logspaced = f.logspaced,
    categorical = f.categorical,
)


text:
function dummy_node(up, nid)
    act = -1
    states = LightDark1DState(0, 0)
    ae = [MCVI.AlphaEdge([1], nid)]
    return MCVI.create_node(up, act, states, ae)
end

function dummy_node2(up, nid1, nid2)
    act = 1
    states = [LightDark1DState(0, 0), LightDark1DState(0, 3), LightDark1DState(0, 4.8), LightDark1DState(0, 5.2)]
    ae = [MCVI.AlphaEdge([-1, -1, 200, 200], nid1), MCVI.AlphaEdge([0, -1, -1, -1], nid2)]
    return MCVI.create_node(up, act, states, ae)
end

function build_dummy_graph(pomdp)
    up = MCVIUpdater(pomdp)
    n0 = MCVI.init_node(up, pomdp)
    MCVI.addnode!(up, n0)
    n1 = dummy_node(up, n0.id)
    MCVI.addnode!(up, n1)
    n2 = dummy_node(up, n1.id)
    MCVI.addnode!(up, n2)
    n3 = dummy_node(up, n2.id)
    MCVI.addnode!(up, n3)
    up.root = n3
    return up
end

function test_dummy_graph()
    rng = MersenneTwister(42)
    p = LightDark1D()
    up = build_dummy_graph(p)
    # s0 = initialstate(p)
    # sim = MCVISimulator(MersenneTwister(420), 1)
    sim = MCVISimulator()
    policy = MCVIPolicy(p, up)

    sumv = 0.0
    for i in 1:1000
        ss = initialstate(p, rng)
        sumv += MCVI.simulate(sim, p, policy, up, up.root, ss)
    end
    sumv /= 1000
    return abs(sumv - (-3.0)) <= 1
end

function test_dummy_graph2()
    rng = MersenneTwister(42)
    pomdp = LightDark1D()
    up = MCVIUpdater(pomdp)
    n0 = MCVI.init_node(up, pomdp)
    MCVI.addnode!(up, n0)
    n1 = dummy_node(up, n0.id)
    MCVI.addnode!(up, n1)
    n2 = dummy_node(up, n1.id)
    MCVI.addnode!(up, n2)
    n3 = dummy_node(up, n2.id)
    MCVI.addnode!(up, n3)
    n4 = dummy_node(up, n3.id)
    MCVI.addnode!(up, n4)
    n5 = dummy_node(up, n4.id)
    MCVI.addnode!(up, n5)
    # Move lefts 5 steps
    rn0 = dummy_node2(up, n5.id, n0.id)
    MCVI.addnode!(up, rn0)
    rn1 = dummy_node2(up, n5.id, rn0.id)
    MCVI.addnode!(up, rn1)
    rn2 = dummy_node2(up, n5.id, rn1.id)
    MCVI.addnode!(up, rn2)
    up.root = rn2

    sim = MCVISimulator()
    policy = MCVIPolicy(pomdp, up)
    dump_json(policy, joinpath(tempdir(), "test_policy.json"))
    sumv = 0.0
    for i in 1:1000
        ss = initialstate(pomdp, rng)
        sumv += MCVI.simulate(sim, pomdp, policy, up, up.root, ss)
    end
    sumv /= 1000
    println(sumv)
    return abs(sumv - (-1)) <= 1
end


text:
using .Turing

include("turing_plotcoef.jl")


text:
############################################################################################
"""
$(TYPEDEF)

Default arguments for PMCMC constructor.

# Fields
$(TYPEDFIELDS)
"""
struct ParticleGibbs{P<:ParticleFilter,M<:MCMC} <: PMCMCKernel
    "Particle Filter kernel to estimate latent trajectory."
    pf::P
    "MCMC kernel to sample continuous model parameter."
    mcmc::M
    function ParticleGibbs(pf::P, mcmc::M) where {P<:ParticleFilter,M<:MCMC}
        @argcheck !isa(pf.tune.referencing, Marginal) "Cannot use Marginal Filter in Conditional setting, use Ancestral or Conditional referencing"
        return new{P,M}(pf, mcmc)
    end
end
function ParticleGibbs(
    filter::F, mcmc::M; kwargs...
) where {F<:ParticleFilterConstructor,M<:MCMCConstructor}
    return PMCMC(ParticleGibbs, filter, mcmc; kwargs...)
end

############################################################################################
"""
$(SIGNATURES)
Propose new parameter with mcmc psampler. If update=true, objective function will be updated with input model and data.

# Examples
```julia
```

"""
function propose!(
    _rng::Random.AbstractRNG,
    pmcmc::ParticleGibbs,
    model::ModelWrapper,
    data::D,
    temperature::F = model.info.reconstruct.default.output(1.0),
    update::U=BaytesCore.UpdateTrue(),
) where {D,F<:AbstractFloat, U<:BaytesCore.UpdateBool}
    ## Get trajectory via PF - always update data.latent in model
    _, pf_diagnostics = propose!(_rng, pmcmc.pf, model, data, temperature, update)
    ## Propose new θₜ - if accepted, model is updated accordingly
    _, mcmc_diagnostics = propose!(_rng, pmcmc.mcmc, model, data, temperature, update)
    ## Assign base diagnostics - ℓobjective and predictions are taken from particle filter
    diagnostics = BaytesCore.BaseDiagnostics(
        pf_diagnostics.base.ℓobjective, pf_diagnostics.base.temperature, pf_diagnostics.base.prediction,
        mcmc_diagnostics.base.iter
    )
    ## Return pmcmc output
    return model.val, PMCMCDiagnostics(diagnostics, pf_diagnostics, mcmc_diagnostics)
end

############################################################################################
#export
export ParticleGibbs, propose!


text:
using LowRankModels

m,n,k = 10,20,3
Y = rand(3,20)
A = rand(10,20)

ry = Regularizer[fixed_latent_features(Y[:,i]) for i=1:n]
glrm = GLRM(A, QuadLoss(), SimplexConstraint(), ry, k)
X, Yp, ch = fit!(glrm)

@assert(Yp == Y)


text:
ui = plot_dict(simParams)
title = """
Simulation parameters
"""
content!(w, "div#content", title*ui)


text:
@testset "linesegments" begin
    @testset "intersection" begin
        p1 = Point2D(1.0, 3.0)
        p2 = Point2D(3.0, 1.0)
        p3 = Point2D(1.0, 1.0)
        p4 = Point2D(3.0, 3.0)

        s1 = LineSegment2D(p1, p2)
        s2 = LineSegment2D(p3, p4)

        @test intersection(s1,s2)
    end
end

max_stars_repo_name: helgith/GeometricStructures

max_issues_repo_name: helgith/GeometricStructures

max_forks_repo_name: helgith/GeometricStructures


text:
## Algorithms
## =======

# SMT Based
immutable SigmaSMT <: Algorithm
  capture::Vector{Symbol}
  randvartype::DataType
  preimage_sampler::Function
  ncores::Int
  split::Function
end

immutable SigmaRej <: Algorithm
  capture::Vector{Symbol}
end

immutable SigmaBFS <: Algorithm
  capture::Vector{Symbol}
end

immutable SigmaMH <: Algorithm
  capture::Vector{Symbol}
end

# eq_f{T}(x::T,y::T) = T.names == () ? (==) : deepequiv

# function deepequiv{T}(a::T,b::T)
#   for f in T.names
#     if Base.isdefined(a,f) && Base.isdefined(b,f)
#       same = eq_f(getfield(a,f), getfield(b,f))(getfield(a,f),getfield(b,f))
#       !same && return false
#     elseif Base.isdefined(a,f) $ Base.isdefined(b,f)
#       return false
#     end
#   end
#   return true
# end

# function deephash{T}(x::T, h = zero(Uint))
#   h += uint(0x7f53e68ceb575e76)
#   for t in T.names
#     if Base.isdefined(x,t)
#       h = hash(getfield(x,t),h)
#     end
#   end
#   return h
# end

# equiv{T}(a::T,b::T) =
#   all([getfield(a,f) == getfield(b,f) for f in T.names])

# # ==(a::SigmaAI,b::SigmaAI) = equiv(a,b)
# # hash(a::SigmaAI, h::Uint) = deephash(a,h)
# ==(a::SigmaSMT,b::SigmaSMT) = equiv(a,b)
# hash(a::SigmaSMT, h::Uint) = deephash(a,h)

# function isequal(a::SigmaSMT, b::SigmaSMT)
#


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize_nonconvex.jl
# Contains the optimize! routine and subroutines needed in the branch and
# bound routine called by EAGO.
#############################################################################

"""
$(SIGNATURES)

Selects node with the lowest lower bound in stack & loads `_lower_variable_bound`
and `_upper_variable_bound` fields.
"""
function node_selection!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    m._node_count -= 1
    m._current_node = popmin!(m._stack)
    i = 1
    for j in m._branch_variable
        m._lower_variable_bound[j] = lower_variable_bound(m._current_node, i)
        m._upper_variable_bound[j] = upper_variable_bound(m._current_node, i)
        i += 1
    end
    return nothing
end

"""
$(SIGNATURES)

Stores the current node to the stack after updating lower/upper bounds.
"""
function single_storage!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    n = m._current_node
    m._node_repetitions += 1
    m._node_count += 1
    lower_bound = max(n.lower_bound, m._lower_objective_value)
    upper_bound = min(n.upper_bound, m._upper_objective_value)
    push!(m._stack, NodeBB(ntuple(i -> _lower_bound(BranchVar, m, i), N),
                           ntuple(i -> _upper_bound(BranchVar, m, i), N),
                           lower_bound, upper_bound, n.depth, n.id,
                           n.last_branch, n.branch_direction))
    return nothing
end

"""
$(SIGNATURES)

Selects and deletes nodes from stack with lower bounds greater than global
upper bound.
"""
function fathom!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    upper = m._global_upper_bound
    continue_flag = !isempty(m._stack)
    while continue_flag
        max_node = maximum(m._stack)
        max_check = (max_node.lower_bound > upper)
        if max_check
            popmax!(m._stack)
            m._node_count -= 1
            if isempty(m._stack)
                continue_flag = false
            end
        else
            if !max_check
                continue_flag = false
            elseif isempty(m._stack)
                continue_flag = false
            end
        end
    end
    return nothing
end

"""
$(SIGNATURES)

Checks to see if current node should be reprocessed.
"""
repeat_check(t::ExtensionType, m::GlobalOptimizer) = false

relative_gap(L::Float64, U::Float64) = ((L > -Inf) && (U < Inf)) ?  abs(U - L)/(max(abs(L), abs(U))) : Inf
relative_tolerance(L::Float64, U::Float64, tol::Float64) = relative_gap(L, U)  > tol || ~(L > -Inf)

function _terminate_store!(m::GlobalOptimizer, t::MOI.TerminationStatusCode,
                           r::MOI.ResultStatusCode, s::String)
    m._termination_status_code = t
    m._result_status_code = r
    (_verbosity(m) >= 3) && println(s)
    return nothing
end

"""
$(SIGNATURES)

Checks for termination of algorithm due to satisfying absolute or relative
tolerance, infeasibility, or a specified limit, returns a boolean valued true
if algorithm should continue.
"""
function termination_check(t::ExtensionType, m::GlobalOptimizer)
    if length(m._stack) == 0
        if m._first_solution_node > 0
            _terminate_store!(m, MOI.OPTIMAL, MOI.FEASIBLE_POINT, "Empty Stack: Exhaustive Search Finished")
        else
            _terminate_store!(m, MOI.INFEASIBLE,MOI.INFEASIBILITY_CERTIFICATE, "Empty Stack: Infeasible")
        end
    elseif length(m._stack) >= m._node_limit
        _terminate_store!(m, MOI.NODE_LIMIT, MOI.UNKNOWN_RESULT_STATUS, "Node Limit Exceeded")
    elseif m._iteration_count >= m.iteration_limit
        _terminate_store!(m, MOI.ITERATION_LIMIT, MOI.UNKNOWN_RESULT_STATUS, "Maximum Iteration Exceeded")
    elseif _relative_gap(m) < m._relative_tolerance
        _terminate_store!(m, MOI.OPTIMAL, MOI.FEASIBLE_POINT, "Relative Tolerance Achieved")
    elseif _relative_gap(m) < m._absolute_tolerance
        _terminate_store!(m, MOI.OPTIMAL, MOI.FEASIBLE_POINT, "Absolute Tolerance Achieved")
    elseif m._run_time > m.time_limit
        _terminate_store!(m, MOI.TIME_LIMIT, MOI.UNKNOWN_RESULT_STATUS, "Time Limit Exceeded")
    else
        return false
    end
    return true
end

"""
$(SIGNATURES)

Checks for convergence of algorithm with respect to absolute and/or relative
tolerances.
"""
function convergence_check(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

  L = m._lower_objective_value
  U = m._global_upper_bound
  t = (U - L) <= _absolute_tolerance(m)
  if (U < Inf) && (L > Inf)
      t |= (abs(U - L)/(max(abs(L), abs(U))) <= _relative_tolerance(m))
  end

  if t && m._min_converged_value < Inf
      m._min_converged_value = min(m._min_converged_value, L)
  else
      m._min_converged_value = L
  end

  return t
end

"""
$(SIGNATURES)

Provides a hook for extensions to EAGO as opposed to standard global, local,
or linear solvers.
"""
optimize_hook!(t::ExtensionType, m::GlobalOptimizer) = nothing

function store_candidate_solution!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    if m._upper_feasibility && (m._upper_objective_value < m._global_upper_bound)
        m._feasible_solution_found = true
        m._first_solution_node = m._maximum_node_id
        m._global_upper_bound = m._upper_objective_value
        @__dot__ m._solution = m._upper_solution
    end
    return nothing
end

function set_global_lower_bound!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    if !isempty(m._stack)
        min_node = minimum(m._stack)
        l = min_node.lower_bound
        if m._global_lower_bound < l
            m._global_lower_bound = l
        end
    end
    return nothing
end

parse_global!(t::ExtensionType, m::GlobalOptimizer) = nothing

# wraps subroutine call to isolate ExtensionType
for f in (:parse_global!, :presolve_global!, :termination_check, :cut_condition,
          :convergence_check, :repeat_check, :node_selection!, :preprocess!,
          :lower_problem!, :add_cut!, :upper_problem!, :postprocess!,
          :single_storage, :branch_node!, :fathom!)
    @eval function ($f)(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
        ($f)(m._ext_type, m)
    end
end

"""
$(TYPEDSIGNATURES)

Solves the branch and bound problem with the input EAGO optimizer object.
"""
function global_solve!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    m._iteration_count = 1
    m._node_count = 1

    parse_global!(m)
    presolve_global!(m)

    logging_on = m.log_on
    verbosity = m.verbosity

    # terminates when max nodes or iteration is reach, or when node stack is empty
    while !termination_check(m)

        # Selects node, deletes it from stack, prints based on verbosity
        node_selection!(m)
        (verbosity >= 3) && print_node!(m)

        # Performs prepocessing and times
        logging_on && (start_time = time())
        preprocess!(m)
        if logging_on
            m._last_preprocess_time = time() - start_time
        end

        if m._preprocess_feasibility

            # solves & times lower bounding problem
            logging_on && (start_time = time())
            m._cut_iterations = 1
            lower_problem!(m)
            while cut_condition(m)
                add_cut!(m)
            end
            if logging_on
                m._last_lower_problem_time = time() - start_time
            end
            print_results!(m, true)
            print_results_post_cut!(m)

            # checks for infeasibility stores solution
            if m._lower_feasibility
                if !convergence_check(m)

                    logging_on && (start_time = time())
                    upper_problem!(m)
                    if logging_on
                        m._last_upper_problem_time = time() - start_time
                    end
                    print_results!(m, false)
                    store_candidate_solution!(m)
                    if m._input_problem._optimization_sense === MOI.FEASIBILITY_SENSE
                        if !m.feasible_local_continue || m.local_solve_only
                            break
                        end
                    end

                    # Performs and times post processing
                    logging_on && (start_time = time())
                    postprocess!(m)
                    if logging_on
                        m._last_postprocessing_time = time() - start_time
                    end

                    # Checks to see if the node
                    if m._postprocess_feasibility
                        if repeat_check(m)
                            single_storage!(m)
                        else
                            branch_node!(m)
                        end
                    end
                end
            end
            fathom!(m)
        else
            m._lower_objective_value = -Inf
            m._lower_feasibility = false
            m._upper_feasibility = false
        end
        set_global_lower_bound!(m)
        m._run_time = time() - m._start_time
        m._time_left = m.time_limit - m._run_time
        log_iteration!(m)
        print_iteration!(m)
        m._iteration_count += 1
    end

    m._objective_value = m._global_upper_bound

    # Prints the solution
    print_solution!(m)

    return nothing
end

optimize!(::Val{MINCVX}, m::Optimizer) = global_solve!(m._solver)


text:
# ---
# title: 363. Max Sum of Rectangle No Larger Than K
# id: problem363
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Dynamic Programming, Queue
# link: <https://leetcode.com/problems/max-sum-of-rectangle-no-larger-than-k/description/>
# hidden: true
# ---
#
# Given a non-empty 2D matrix _matrix_ and an integer _k_ , find the max sum of
# a rectangle in the _matrix_ such that its sum is no larger than _k_.
#
# **Example:**
#
#
#
#     Input: matrix = [[1,0,1],[0,-2,3]], k = 2
#     Output: 2
#     Explanation:  Because the sum of rectangle [[0, 1], [-2, 3]] is 2,
#                  and 2 is the max number no larger than k (k = 2).
#
# **Note:**
#
#   1. The rectangle inside the matrix must have an area > 0.
#   2. What if the number of rows is much larger than the number of columns?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end