text:
# In Base/intfuncs.jl, x^p returns zero(x) when p == 0
# Since one(PolyVar) and one(Monomial) do not return
# a PolyVar and a Monomial, this results in type instability
# Defining the specific methods solve this problem and also make
# them a lot faster
Base.:(^)(x::PolyVar{C}, i::Int) where {C} = Monomial{C}([x], [i])
Base.:(^)(x::Monomial{true}, i::Int) = Monomial{true}(x.vars, i*x.z)

myminivect(x::T, y::T) where {T} = [x, y]
function myminivect(x::S, y::T) where {S,T}
    U = promote_type(S, T)
    [U(x), U(y)]
end

Base.:(+)(x::DMonomialLike, y::DMonomialLike) = Term(x) + Term(y)
Base.:(-)(x::DMonomialLike, y::DMonomialLike) = Term(x) - Term(y)

_getindex(p::Polynomial, i) = p[i]
_getindex(t::Term, i) = t
function _plusorminus_to!(a::Vector{U}, Z::Vector{Vector{Int}}, op::Function, p::TermPoly{C}, q::TermPoly{C}, maps, nvars) where {C, U}
    i = j = 1
    while i <= nterms(p) || j <= nterms(q)
        z = zeros(Int, nvars)
        if j > nterms(q) || (i <= nterms(p) && _getindex(p, i).x > _getindex(q, j).x)
            t = _getindex(p, i)
            z[maps[1]] = t.x.z
            α = MA.scaling_convert(U, MA.copy_if_mutable(t.α))
            i += 1
        elseif i > nterms(p) || _getindex(q, j).x > _getindex(p, i).x
            t = _getindex(q, j)
            z[maps[2]] = t.x.z
            α = MA.scaling_convert(U, MA.operate(op, t.α))
            j += 1
        else
            t = _getindex(p, i)
            z[maps[1]] = t.x.z
            s = _getindex(q, j)
            α = op(t.α, s.α)
            i += 1
            j += 1
        end
        push!(a, α)
        push!(Z, z)
    end
end
function plusorminus(p::TermPoly{C, S}, q::TermPoly{C, T}, op::Function) where {C, S, T}
    varsvec = [_vars(p), _vars(q)]
    allvars, maps = mergevars(varsvec)
    U = MA.promote_operation(op, S, T)
    a = U[]
    Z = Vector{Int}[]
    _plusorminus_to!(a, Z, op, p, q, maps, length(allvars))
    Polynomial(a, MonomialVector{C}(allvars, Z))
end

function MA.mutable_operate_to!(output::Polynomial{C}, op::Union{typeof(+), typeof(-)},
                                p::TermPoly{C}, q::TermPoly{C}) where C
    if output === p || output === q
        # Otherwise, `_plusorminus_to!` never finishes
        error("Cannot call `mutable_operate_to!(output, $op, p, q)` with `output` equal to `p` or `q`, call `mutable_operate!` instead.")
    end
    varsvec = [_vars(p), _vars(q)]
    allvars, maps = mergevars(varsvec)
    Future.copy!(output.x.vars, allvars)
    empty!(output.a)
    empty!(output.x.Z)
    _plusorminus_to!(output.a, output.x.Z, op, p, q, maps, length(allvars))
    return output
end

function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial,
                             q::Union{PolyVar, Monomial, Term})
    return MA.mutable_operate!(op, p, polynomial(q))
end
const _NoVarTerm{T} = Tuple{T, Vector{Int}}
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{false},
                             q::Polynomial{false})
    return MA.mutable_operate_to!(p, op, MA.mutable_copy(p), q)
end
# TODO need to check that this also works for non-commutative
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{true},
                             q::Polynomial{true})
    if _vars(p) != _vars(q)
        varsvec = [_vars(p), _vars(q)]
        allvars, maps = mergevars(varsvec)
        if length(allvars) != length(_vars(p))
            _add_variables!(p.x, allvars, maps[1])
        end
        if length(allvars) == length(_vars(q))
            rhs = q
        else
            # We could avoid promoting `q` to the same variables
            # like in `plusorminus` to avoid extra allocation but it then
            # gives slower comparison. There is a tradeoff and the approach used here
            # should be better of `q` has less terms and then the same term is compared
            # many times.
            rhs = Polynomial(q.a, copy(q.x))
            _add_variables!(rhs.x, allvars, maps[2])
        end
        return MA.mutable_operate!(op, p, rhs)
    end
    get1(i) = (p.a[i], p.x.Z[i])
    get2(i) = (MA.scaling_convert(eltype(p.a), MA.operate(op, q.a[i])), copy(q.x.Z[i]))
    function set(i, t::_NoVarTerm)
        p.a[i] = t[1]
        p.x.Z[i] = t[2]
    end
    function push(t::_NoVarTerm)
        push!(p.a, t[1])
        push!(p.x.Z, t[2])
    end
    compare_monomials(t::_NoVarTerm, j::Int) = samevars_grlex(q.x.Z[j], t[2])
    compare_monomials(i::Int, j::Int) = compare_monomials(get1(i), j)
    combine(i::Int, j::Int) = p.a[i] = MA.operate!(op, p.a[i], q.a[j])
    combine(t::_NoVarTerm, j::Int) = (MA.operate!(op, t[1], q.a[j]), t[2])
    function resize(n)
        resize!(p.a, n)
        resize!(p.x.Z, n)
    end
    # We can modify the coefficient since it's the result of `combine`.
    keep(t::_NoVarTerm) = !MA.iszero!(t[1])
    keep(i::Int) = !MA.iszero!(p.a[i])
    MP.polynomial_merge!(
        nterms(p), nterms(q), get1, get2, set, push,
        compare_monomials, combine, keep, resize
    )
    return p
end

Base.:(+)(x::TermPoly{C}, y::TermPoly{C}) where C = plusorminus(x, y, +)
Base.:(-)(x::TermPoly{C}, y::TermPoly{C}) where C = plusorminus(x, y, -)
Base.:(+)(x::TermPoly{C}, y::Union{Monomial,PolyVar}) where C = x + Term{C}(y)
Base.:(+)(x::Union{Monomial,PolyVar}, y::TermPoly{C}) where C = Term{C}(x) + y

Base.:(-)(x::TermPoly{T}, y::DMonomialLike) where T = x - Term{T}(y)
Base.:(-)(x::DMonomialLike, y::TermPoly{T}) where T = Term{T}(x) - y

# `MA.operate(-, p)` redirects to `-p` as it assumes that `-p` can be modified
# through the MA API without modifying `p`. We should either copy `p.x` here
# or implement a `MA.operate(-, p)` that copies it. We choose the first option.
Base.:(-)(p::Polynomial) = Polynomial(-p.a, copy(p.x))
# TODO use that with MA v0.2
#Base.:(-)(p::Polynomial) = Polynomial(map(α -> MA.operate(-, α), p.a), copy(p.x))

include("mult.jl")


text:
@testset "MASH distances" begin
    a = minhash(dna"ATCGCCA-", 4, 3)
    b = minhash(dna"ATCGCCTA", 4, 3)
    @test mash(a, b) ≈ 0.2745 atol=0.001
    @test mash(a, a) == 0
    @test a.sketch == sort(a.sketch)
end

ext: jl


text:
using MPI
using Trilinos
using Test

function build_graph!(graph)
  rowmap = Tpetra.getRowMap(graph)
  n_my_elms = Tpetra.getNodeNumElements(rowmap)
  gmin = Tpetra.getMinGlobalIndex(rowmap)
  gmax = Tpetra.getMaxGlobalIndex(rowmap)

  row_indices = zeros(Int,5)
  row_indices_view = Teuchos.ArrayView(row_indices)

  for i in 0:n_my_elms-1
    global_row =  Tpetra.getGlobalElement(rowmap,i)
    for j in 1:5
      row_indices[j] = global_row - (3-j)
      if row_indices[j] < gmin || row_indices[j] > gmax
        row_indices[j] = global_row
      end
    end
      Tpetra.insertGlobalIndices(graph, global_row, row_indices)
  end
end

function fill_matrix!(A)
  rowmap = Tpetra.getRowMap(A)
  n_my_elms = Int(Tpetra.getNodeNumElements(rowmap))

  row_indices_arr = zeros(Int32,5)
  row_values_arr = zeros(5)

  row_indices = Teuchos.ArrayView(row_indices_arr)
  row_values = Teuchos.ArrayView(row_values_arr)
  n_elems = Ref(UInt(0))

  Tpetra.resumeFill(A)
  for row in 0:n_my_elms-1
    Tpetra.getLocalRowCopy(A, row, row_indices, row_values, n_elems)
    for j in 1:n_elems[]
      row_values_arr[j] = j
    end
    Tpetra.replaceLocalValues(A, row, row_indices, row_values)
  end
end

"""
Assemble the 2D Poisson problem on a structured grid
"""
function bench_assembly(comm, nnodes)
  # Construct map
  rowmap = Tpetra.Map(nnodes, 0, comm)

  # Construct graph (i.e. the sparsity pattern)
  matrix_graph = Tpetra.CrsGraph(rowmap, 5)
  println("-----------------------------------")
  println("Graph construction time:")
  @time build_graph!(matrix_graph)
  Tpetra.fillComplete(matrix_graph)

  # Matrix construction
  A = Tpetra.CrsMatrix(matrix_graph)
  Tpetra.resumeFill(A)
  println("Matrix fill time:")
  @time fill_matrix!(A)
  Tpetra.fillComplete(A)

  # This prints the matrix if uncommented
  #Tpetra.describe(A, Teuchos.VERB_EXTREME)
end

# MPI setup
if !MPI.Initialized()
  MPI.Init()
  MPI.finalize_atexit()
end
comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))

my_rank = Teuchos.getRank(comm)
if my_rank != 0
  redirect_stdout(open("/dev/null", "w"))
  redirect_stderr(open("/dev/null", "w"))
end

N = 1000000
bench_assembly(comm, N)
bench_assembly(comm, N)
bench_assembly(comm, N)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterface <: SwaggerModel
    iface_id::Union{ Nothing, String } # iface_id
    guest_mac::Union{ Nothing, String } # guest_mac
    host_dev_name::Union{ Nothing, String } # host_dev_name
    allow_mmds_requests::Union{ Nothing, Bool } # allow_mmds_requests
    rx_rate_limiter::Union{ Nothing, RateLimiter } # rx_rate_limiter
    tx_rate_limiter::Union{ Nothing, RateLimiter } # tx_rate_limiter

    function NetworkInterface(;iface_id=nothing, guest_mac=nothing, host_dev_name=nothing, allow_mmds_requests=nothing, rx_rate_limiter=nothing, tx_rate_limiter=nothing)
        o = new()
        set_field!(o, :iface_id, iface_id)
        set_field!(o, :guest_mac, guest_mac)
        set_field!(o, :host_dev_name, host_dev_name)
        set_field!(o, :allow_mmds_requests, allow_mmds_requests)
        set_field!(o, :rx_rate_limiter, rx_rate_limiter)
        set_field!(o, :tx_rate_limiter, tx_rate_limiter)
        o
    end
end # type NetworkInterface

const _name_map_NetworkInterface = Dict{String,Symbol}(["iface_id"=>:iface_id, "guest_mac"=>:guest_mac, "host_dev_name"=>:host_dev_name, "allow_mmds_requests"=>:allow_mmds_requests, "rx_rate_limiter"=>:rx_rate_limiter, "tx_rate_limiter"=>:tx_rate_limiter])
const _field_map_NetworkInterface = Dict{Symbol,String}([:iface_id=>"iface_id", :guest_mac=>"guest_mac", :host_dev_name=>"host_dev_name", :allow_mmds_requests=>"allow_mmds_requests", :rx_rate_limiter=>"rx_rate_limiter", :tx_rate_limiter=>"tx_rate_limiter"])
Swagger.name_map(::Type{ NetworkInterface }) = _name_map_NetworkInterface
Swagger.field_map(::Type{ NetworkInterface }) = _field_map_NetworkInterface

function check_required(o::NetworkInterface)
    (o.iface_id === nothing) && (return false)
    (o.host_dev_name === nothing) && (return false)
    true
end

function validate_field(o::NetworkInterface, name::Symbol, val)
end


text:
using MatchCore

# example assuming * operation is always binary

# ENV["JULIA_DEBUG"] = Metatheory

abstract type NumberFold <: AbstractAnalysis end

# This should be auto-generated by a macro
function EGraphs.make(an::Type{NumberFold}, g::EGraph, n::ENode)
    # @show n
    n.head isa Number && return n.head

    op_times = addexpr!(g, :*)
    op_plus = addexpr!(g, :+)

    if n.head == :call && arity(n) == 3
        op = geteclass(g, n.args[1])
        l = geteclass(g, n.args[2])
        r = geteclass(g, n.args[3])
        ldata = getdata(l, an, nothing)
        rdata = getdata(r, an, nothing)

        # @show ldata rdata

        if ldata isa Number && rdata isa Number
            if in_same_class(g, op, op_times)
                return ldata * rdata
            elseif in_same_class(g, op, op_plus)
                return ldata + rdata
            end
        end
    end

    return nothing
end

function EGraphs.join(an::Type{NumberFold}, from, to)
    # println("joining!")
    if from isa Number
        if to isa Number
            @assert from == to
        else return from
        end
    end
    return to
end

function EGraphs.modify!(an::Type{NumberFold}, g::EGraph, id::Int64)
    # !haskey(an, id) && return nothing
    eclass = g.classes[id]
    d = getdata(eclass, an, nothing)
    if d isa Number
        newclass = addexpr!(g, d)
        merge!(g, newclass.id, id)
    end
end

EGraphs.islazy(x::Type{NumberFold}) = false

comm_monoid = @theory begin
    a * b => b * a
    a * 1 => a
    a * (b * c) => (a * b) * c
end

G = EGraph(:(3 * 4))
analyze!(G, NumberFold)

# exit(0)

@testset "Basic Constant Folding Example - Commutative Monoid" begin
    @test (true == @areequalg G comm_monoid 3 * 4 12)

    @test (true == @areequalg G comm_monoid 3 * 4 12 4*3  6*2)
end

# Metatheory.EGraphs.PRINTIT[] = true
# Metatheory.options.verbose = true
# Metatheory.options.printiter = true

@testset "Basic Constant Folding Example 2 - Commutative Monoid" begin
    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    analyze!(G, NumberFold)
    addexpr!(G, :(12*a))
    println(saturate!(G, comm_monoid))
    display(G.classes); println()
    @test (true == @areequalg G comm_monoid (12*a)*b ((6*2)*b)*a)
    @test (true == @areequalg G comm_monoid (3 * a) * (4 * b) (12*a)*b ((6*2)*b)*a)
end

@testset "Basic Constant Folding Example - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    # addexpr!(G, 12)
    saturate!(G, comm_monoid)
    addexpr!(G, :(a * 2))
    analyze!(G, NumberFold)
    saturate!(G, comm_monoid)

    # display(G.classes); println()
    # println(G.root)
    # display(G.analyses[1].data); println()

    @test (true == areequal(G, comm_monoid, :(3 * 4), 12, :(4*3), :(6*2)))

    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    analyze!(G, NumberFold)
    params=SaturationParams(timeout=15)
    @test areequal(G, comm_monoid, :((3 * a) * (4 * b)), :((12*a)*b),
        :(((6*2)*b)*a); params=params)
end

@testset "Infinite Loops analysis" begin
    boson = @theory begin
        1 * x => x
    end

    G = EGraph(Util.cleanast( :(1 * x) ))
    params=SaturationParams(timeout=100)
    saturate!(G,boson, params)
    ex = extract!(G, ExtractionAnalysis{astsize})

    # println(ex)

    using Metatheory.EGraphs
    boson = @theory begin
        (:c * :cdag) => :cdag * :c + 1
        a * (b + c) => (a * b) + (a * c)
        (b + c) * a => (b * a) + (c * a)
        # 1 * x => x
        (a * b) * c => a * (b * c)
        a * (b * c) => (a * b) * c
    end

    G = EGraph(Util.cleanast( :(c * c * cdag * cdag) ))
    saturate!(G,boson)
    ex = extract!(G, astsize_inv)

    # println(ex)
end


text:
using ..DGmethods.NumericalFluxes:
    Rusanov, CentralNumericalFluxGradient, CentralNumericalFluxDiffusive

abstract type OceanBoundaryCondition end

"""
    Defining dummy structs to dispatch on for boundary conditions.
"""
struct CoastlineFreeSlip <: OceanBoundaryCondition end
struct CoastlineNoSlip <: OceanBoundaryCondition end
struct OceanFloorFreeSlip <: OceanBoundaryCondition end
struct OceanFloorNoSlip <: OceanBoundaryCondition end
struct OceanSurfaceNoStressNoForcing <: OceanBoundaryCondition end
struct OceanSurfaceStressNoForcing <: OceanBoundaryCondition end
struct OceanSurfaceNoStressForcing <: OceanBoundaryCondition end
struct OceanSurfaceStressForcing <: OceanBoundaryCondition end

"""
    CoastlineFreeSlip

applies boundary condition ν∇u = 0 and κ∇θ = 0
"""

"""
    ocean_boundary_state!(::HBModel, ::CoastlineFreeSlip, ::Rusanov)

apply free slip boundary conditions for velocity
apply no penetration boundary for temperature

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::Rusanov,
    ::CoastlineFreeSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    u⁻ = Q⁻.u
    n = @SVector [n⁻[1], n⁻[2]]

    Q⁺.u = u⁻ - 2 * (n ⋅ u⁻) * n

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::CoastlineFreeSlip, ::CentralNumericalFluxGradient)

apply free slip boundary conditions for velocity
apply no penetration boundary for temperature

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxGradient,
    ::CoastlineFreeSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    u⁻ = Q⁻.u
    n = @SVector [n⁻[1], n⁻[2]]

    Q⁺.u = u⁻ - (n ⋅ u⁻) * n

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::CoastlineFreeSlip, ::CentralNumericalFluxDiffusive)

apply free slip boundary conditions for velocity
apply no penetration boundary for temperature
sets ghost point to have no numerical flux on the boundary for ν∇u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::CoastlineFreeSlip,
    ::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    D⁺.ν∇u = -D⁻.ν∇u

    D⁺.κ∇θ = -D⁻.κ∇θ

    return nothing
end

"""
    CoastlineNoSlip

applies boundary condition u = 0 and κ∇θ = 0
"""

"""
    ocean_boundary_state!(::HBModel, ::CoastlineNoSlip, ::Rusanov)

apply no slip boundary condition for velocity
apply no penetration boundary for temperature
set sets ghost point to have no numerical flux on the boundary for u

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::Rusanov,
    ::CoastlineNoSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    Q⁺.u = -Q⁻.u

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::CoastlineNoSlip, ::CentralNumericalFluxGradient)

apply no slip boundary condition for velocity
apply no penetration boundary for temperature
set numerical flux to zero for u

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxGradient,
    ::CoastlineNoSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    FT = eltype(Q⁺)
    Q⁺.u = SVector(-zero(FT), -zero(FT))

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::CoastlineNoSlip, ::CentralNumericalFluxDiffusive)

apply no slip boundary condition for velocity
apply no penetration boundary for temperature
sets ghost point to have no numerical flux on the boundary for u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::CoastlineNoSlip,
    ::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = -Q⁻.u

    D⁺.κ∇θ = -D⁻.κ∇θ

    return nothing
end

"""
    OceanFloorFreeSlip

applies boundary condition ν∇u = 0 and κ∇θ = 0
"""

"""
    ocean_boundary_state!(::HBModel, ::OceanFloorFreeSlip, ::Rusanov)

apply free slip boundary conditions for velocity
apply no penetration boundary for temperature
set ghost point to have no numerical flux on the boundary for w

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::Rusanov,
    ::OceanFloorFreeSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    A⁺.w = -A⁻.w

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanFloorFreeSlip, ::CentralNumericalFluxGradient)

apply free slip boundary condition for velocity
apply no penetration boundary for temperature
set numerical flux to zero for w

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxGradient,
    ::OceanFloorFreeSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    FT = eltype(Q⁺)
    A⁺.w = -zero(FT)

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanFloorFreeSlip, ::CentralNumericalFluxDiffusive)

apply free slip boundary conditions for velocity
apply no penetration boundary for temperature
sets ghost point to have no numerical flux on the boundary for ν∇u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::OceanFloorFreeSlip,
    ::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    A⁺.w = -A⁻.w
    D⁺.ν∇u = -D⁻.ν∇u

    D⁺.κ∇θ = -D⁻.κ∇θ

    return nothing
end

"""
    OceanFloorNoSlip

applies boundary condition u = 0 and κ∇θ = 0
"""

"""
    ocean_boundary_state!(::HBModel, ::OceanFloorNoSlip, ::Rusanov)

apply no slip boundary condition for velocity
apply no penetration boundary for temperature
set sets ghost point to have no numerical flux on the boundary for u and w

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::Rusanov,
    ::OceanFloorNoSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    Q⁺.u = -Q⁻.u
    A⁺.w = -A⁻.w

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanFloorNoSlip, ::CentralNumericalFluxGradient)

apply no slip boundary condition for velocity
apply no penetration boundary for temperature
set numerical flux to zero for u and w

# Arguments
- `Q⁺`: state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxGradient,
    ::OceanFloorNoSlip,
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    FT = eltype(Q⁺)
    Q⁺.u = SVector(-zero(FT), -zero(FT))
    A⁺.w = -zero(FT)

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanFloorNoSlip, ::CentralNumericalFluxDiffusive)

apply no slip boundary condition for velocity
apply no penetration boundary for temperature
sets ghost point to have no numerical flux on the boundary for u,w and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::OceanFloorNoSlip,
    ::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)

    Q⁺.u = -Q⁻.u
    A⁺.w = -A⁻.w

    D⁺.κ∇θ = -D⁻.κ∇θ

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::Union{OceanSurface*}, ::Union{Rusanov, CentralNumericalFluxGradient})

applying neumann boundary conditions, so don't need to do anything for these numerical fluxes
"""
@inline function ocean_boundary_state!(
    ::Union{Rusanov, CentralNumericalFluxGradient},
    ::Union{
        OceanSurfaceNoStressNoForcing,
        OceanSurfaceStressNoForcing,
        OceanSurfaceNoStressForcing,
        OceanSurfaceStressForcing,
    },
    ::HBModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanSurfaceNoStressNoForcing, ::CentralNumericalFluxDiffusive)

apply no flux boundary condition for velocity
apply no flux boundary condition for temperature
set ghost point to have no numerical flux on the boundary for ν∇u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::OceanSurfaceNoStressNoForcing,
    ::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    D⁺.ν∇u = -D⁻.ν∇u
    D⁺.κ∇θ = -D⁻.κ∇θ

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanSurfaceStressNoForcing, ::CentralNumericalFluxDiffusive)

apply wind-stress boundary condition for velocity
apply no flux boundary condition for temperature
set ghost point for numerical flux on the boundary for ν∇u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::OceanSurfaceStressNoForcing,
    m::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    τᶻ = kinematic_stress(m.problem, A⁻.y, m.ρₒ)
    τ = @SMatrix [-0 -0; -0 -0; τᶻ -0]
    D⁺.ν∇u = -D⁻.ν∇u + 2 * τ

    D⁺.κ∇θ = -D⁻.κ∇θ

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanSurfaceNoStressForcing, ::CentralNumericalFluxDiffusive)

apply no flux boundary condition for velocity
apply forcing boundary condition for temperature
set ghost point for numerical flux on the boundary for ν∇u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::OceanSurfaceNoStressForcing,
    m::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    D⁺.ν∇u = -D⁻.ν∇u

    σᶻ = temperature_flux(m.problem, A⁻.y, Q⁻.θ)
    σ = @SVector [-0, -0, σᶻ]
    D⁺.κ∇θ = -D⁻.κ∇θ + 2 * σ

    return nothing
end

"""
    ocean_boundary_state!(::HBModel, ::OceanSurfaceStressForcing, ::CentralNumericalFluxDiffusive)

apply wind-stress boundary condition for velocity
apply forcing boundary condition for temperature
set ghost point for numerical flux on the boundary for ν∇u and κ∇θ

# Arguments
- `Q⁺`: state vector at ghost point
- `D⁺`: diffusive state vector at ghost point
- `A⁺`: auxiliary state vector at ghost point
- `n⁻`: normal vector, not used
- `Q⁻`: state vector at interior
- `D⁻`: diffusive state vector at interior point
- `A⁻`: auxiliary state vector at interior point
- `t`:  time, not used
"""
@inline function ocean_boundary_state!(
    ::CentralNumericalFluxDiffusive,
    ::OceanSurfaceStressForcing,
    m::HBModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    τᶻ = kinematic_stress(m.problem, A⁻.y, m.ρₒ)
    τ = @SMatrix [-0 -0; -0 -0; τᶻ -0]
    D⁺.ν∇u = -D⁻.ν∇u + 2 * τ

    σᶻ = temperature_flux(m.problem, A⁻.y, Q⁻.θ)
    σ = @SVector [-0, -0, σᶻ]
    D⁺.κ∇θ = -D⁻.κ∇θ + 2 * σ

    return nothing
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
struct CartPole{S, A, R} <: DiscEnv{S, A, R}
    current::Vector{S}
    observationspace::ObservationSpace{Vector{S}}
    actionspace::ActionSpace{A}
    function CartPole{T}() where T <: AbstractFloat
        observationspace = ObservationSpace{Vector{T}}(4)
        actionspace = ActionSpace([-one(T), one(T)])
        new{T, T, T}(rand(Uniform(-0.05, 0.05), 4), observationspace, actionspace)
    end
end

CartPole() = CartPole{Float32}()

struct CartPoleContinuous{S, A, R} <: ContEnv{S, A, R}
    current::Vector{S}
    observationspace::ObservationSpace{Vector{S}}
    actionspace::ActionSpace{Tuple{A, A}}
    function CartPoleContinuous{T}() where T <: AbstractFloat
        observationspace = ObservationSpace{Vector{T}}(4)
        actionspace = ActionSpace([(-one(T), one(T))])
        new{T, T, T}(rand(Uniform(-0.05, 0.05), 4), observationspace, actionspace)
    end
end

CartPoleContinuous() = CartPoleContinuous{Float32}()

step!(env::CartPole, action::Integer) = step!(env, env.actionspace.actions[action])
step!(env::CartPoleContinuous{S, A, R}, action::Vector{A}) where {S, A, R} = step!(env, action[1])
function step!(env::Union{CartPole{S, A, R}, CartPoleContinuous{S, A, R}}, action::A) where {S, A, R}
    x, ẋ, θ, θ̇ = env.current
    f = 10.0 * action
    θcos = cos(θ)
    θsin = sin(θ)
    tmp = (f + 0.05 * θ̇ * θ̇ * θsin) / 1.1
    θacc = (9.8 * θsin - θcos * tmp) / (0.5 * (4.0/3.0 - 0.1 * θcos * θcos / 1.1))
    xacc = tmp - 0.05 * θacc * θcos / 1.1

    x += 0.02 * ẋ
    ẋ += 0.02 * xacc
    θ += 0.02 * θ̇
    θ̇ += 0.02 * θacc

    env.current .= x, ẋ, θ, θ̇
    θthr = 2.0 * 12.0 * π / 360.0
    done = !(-2.4 <= x <= 2.4 && -θthr <= θ <= θthr)
    return env.current, R(!done), done
end

reset!(env::Union{CartPole, CartPoleContinuous}) = (env.current .= rand(Uniform(-0.05, 0.05), 4))


text:
data = openBF.loadSimulationFiles("test.yml")
# openBF.makeResultsFolder(data)
blood = openBF.buildBlood(data["blood"])
vessels, edges = openBF.buildArterialNetwork(data["network"], blood, 100)

Ccfl = data["solver"]["Ccfl"]
heart = vessels[1].heart

Pref = vessels[1].P[1]
Qref = vessels[1].Q[1]
dt = openBF.calculateDeltaT(vessels, Ccfl)

@test heart.inlet_type == "Q"
@test_nowarn openBF.setInletBC(0.0, dt, vessels[1])
@test vessels[1].Q[1] != Qref
@test vessels[1].P[1] == Pref

@test_nowarn openBF.inletCompatibility(dt, vessels[1], heart)
@test isapprox(vessels[1].A[1], 1.8e-4, atol=1e-5)
@test isapprox(vessels[1].P[1], 1e4, atol=1)

W1, W2 = openBF.riemannInvariants(1, vessels[1])
@test isapprox(W1, -W2, atol=1)

u, c = openBF.inverseRiemannInvariants(W1, W2)
@test isapprox(u, vessels[1].u[1], atol=1e-4)
@test c == vessels[1].c[1]

vessels[1].heart.inlet_type = "P"
openBF.setInletBC(0.0, dt, vessels[1])
@test vessels[1].Q[1] != Qref
@test vessels[1].P[1] != Pref

newA = openBF.areaFromPressure(vessels[1].P[1], vessels[1].A0[1], vessels[1].beta[1],
                                vessels[1].Pext)
@test newA == vessels[1].A[1]

@test isapprox(openBF.inputFromData(0.0, heart), -5.239e-07, atol=1e-7)

@test_nowarn openBF.setOutletBC(dt, vessels[3])
@test isapprox(vessels[3].A[end], 9.5e-5, atol=1e-5)
@test vessels[3].u[end] == 0.0

@test_nowarn openBF.threeElementWindkessel(dt, vessels[1])

M = vessels[1].M
@inferred openBF.updateGhostCells(vessels[1])
@test vessels[1].U00A == vessels[1].A[1]
@test vessels[1].U00Q == vessels[1].Q[1]
@test vessels[1].U01A == vessels[1].A[2]
@test vessels[1].U01Q == vessels[1].Q[2]
@test vessels[1].UM1A == vessels[1].A[M]
@test vessels[1].UM1Q == vessels[1].Q[M]
@test vessels[1].UM2A == vessels[1].A[M-1]
@test vessels[1].UM2Q == vessels[1].Q[M-1]

@test_nowarn openBF.updateGhostCells(vessels)

# cd("..")
# rm("test_results", recursive=true)


text:
using PlasmaDispersionFunctions
using ArbNumerics, QuadGK, Random, SpecialFunctions, Test

import SpecialFunctions.erfcx
erfcx(x::Complex{BigFloat}) = erfcx(ArbComplex(x))

@testset "Tests" begin
Random.seed!(0)

@testset "Plasma dispersion function" begin
  @test plasma_dispersion_function(0.0) ≈ im*sqrt(pi) rtol=1e-5
  @test plasma_dispersion_function(im) ≈ im*0.757872156141312 rtol=1e-5
  @test plasma_dispersion_function(ComplexF64(-1.52, 0.47), 0) ≈
    ComplexF64(0.6088888957234254, 0.33494583882874024) rtol=1e-5
  @test plasma_dispersion_function(big(0.0+0im)) ≈
    ArbComplex(im*sqrt(pi)) rtol=1e-5
end

@testset "using Z₋₁" begin
  Z0 = plasma_dispersion_function(0.0, 0)
  Z1 = plasma_dispersion_function(0.0, 1)
  @test Z1 == plasma_dispersion_function(0.0, 1, Z0)
  Z2 = plasma_dispersion_function(0.0, 2)
  @test Z2 == plasma_dispersion_function(0.0, 2, Z1)
end

@testset "vs quadrature" begin
  function test_quad(Z, pow)
    for z in Z
      principal(x) = x.^pow .* exp(-x.^2)/sqrt(pi)
      integrand(x) = principal(x) / (x-z)
      principalpart = if iszero(imag(z))
        folded(v) = (principal(v + real(z)) - principal(-v + real(z))) / v
        QuadGK.quadgk(folded, nextfloat(0.0), Inf, rtol=eps())[1]
      else
        QuadGK.quadgk(integrand, -Inf, Inf, rtol=eps())[1]
      end
      σ = imag(z) < 0 ? 2 : imag(z) == 0 ? 1 : 0
      cauchyresidue = im * (σ * π * principal(z))
      b = plasma_dispersion_function(z, pow)
      a = principalpart + cauchyresidue
      @test real(a) ≈ real(b) rtol=sqrt(eps()) atol=0.0
      @test imag(a) ≈ imag(b) rtol=sqrt(eps()) atol=0.0
    end
  end
  Zs = []
  push!(Zs,  1.0)
  push!(Zs, -1.0)
  push!(Zs,  1.0 + im)
  push!(Zs, -1.0 + im)
  push!(Zs,  1.0 - im)
  push!(Zs, -1.0 - im)
  for j in 1:10
    push!(Zs, 4 * (rand() - 0.5) + im * 4 * (rand() - 0.5))
  end
  for i in Unsigned.(0:10)
    @testset "QuadGK of maxwellian with $(i)th moment" begin
      test_quad(Zs, i)
    end
  end
end

@testset "check errors are caught" begin
  @test_throws ArgumentError plasma_dispersion_function(1.0, -1)
end

end


text:
using GGPlot
using RDatasets

diamonds = dataset("ggplot2", "diamonds")
g = ggplot(diamonds, aes(x=:Carat, y=:Price)) + geom_point()
display(g)
sleep(20)


text:
function initializeexasim()

version = "src";

# create pde model
pde = Preprocessing.initializepde(version);

# create a mesh structure
mesh = Preprocessing.initializemesh(version);

return pde, mesh

end


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using Pkg; Pkg.activate("JUDI")
using LinearAlgebra, Random
using JUDI, JUDI.TimeModeling, SeisIO, PyCall

pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(dirname(pathof(JUDI)), "Python"))
@pyimport devito as dv
@pyimport numpy as np
@pyimport PyModel as pm
@pyimport JAcoustic_codegen as cg

## Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v0 = ones(Float32,n) .+ 0.4f0
v[:,Int(round(end/2)):end] .= 3f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

## Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
timeR = 1000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell([500f0])
ysrc = convertToCell([0f0])
zsrc = convertToCell([20f0])

# source sampling and number of time steps
timeS = 1000f0  # ms
dtS = 4f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.01f0     # MHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

######################## WITH DENSITY ############################################

# Write shots as segy files to disk
opt = Options(save_data_to_disk=false, file_path=pwd(), file_name="observed_shot", optimal_checkpointing=false)

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*adjoint(Ps), q)

# Nonlinear modeling
dobs = Pr*F*adjoint(Ps)*q
qad = Ps*adjoint(F)*adjoint(Pr)*dobs

# Return wavefields
u = F*adjoint(Ps)*q
v = adjoint(F)*adjoint(Pr)*dobs

# Compute norm
println("forward wavefield 2-norm: ", norm(u))
println("adjoint wavefield 1-norm: ", norm(v, 1))

# Wavefields as source
dnew = Pr*F*v
qnew = Ps*adjoint(F)*u

# Wavefields as source + return wavefields
u2 = F*u
v2 = F*v

# Linearized modeling
dD = J*dm
rtm1 = adjoint(J)*dD


text:
module LazyViews

export AbstractLazyView, MapView, CView, LazyView

using FastAnonymous
typealias LVFunction Union{Function, FastAnonymous.AbstractClosure}

abstract AbstractLazyView;

# Lazy view of a function mapped to an array
immutable MapView{T, N} <: AbstractLazyView
  f::LVFunction;
  a::Array{T, N};

  MapView(f::LVFunction, a::Array{T, N}) = new(f, a);
end

Base.getindex(mv::MapView, idxs...) = mv.f(mv.a[idxs...]);

# Lazy view of a constant multiplied by each element of an array
type CView{T, N} <: AbstractLazyView
  c::Number;
  a::Array{T, N};

  CView(c::Number, a::Array{T, N}) = new(c, a);
end

Base.getindex(cv::CView, idxs...) = cv.c * cv.a[idxs...];

# "Weaker" typed, nestable lazy view
immutable LazyView{T, N} <: AbstractLazyView
  f::LVFunction;
  a::Union{Array{T, N}, LazyView{T, N}};

  LazyView(f::LVFunction, a::Union{Array{T, N}, LazyView{T, N}}) = new(f, a);
end

Base.getindex(lv::LazyView, idxs...) = lv.f(lv.a[idxs...]);
import Base.*
*{T, N}(lv::LazyView{T, N}, c::Number) = LazyView{T, N}(x -> c*x, lv);
*{T, N}(c::Number, lv::LazyView{T, N}) = lv * c;

# Forward size calls to base array
Base.size(amv::AbstractLazyView)          =   size(amv.a);
Base.size(amv::AbstractLazyView, ds...)   =   size(amv.a, ds...);
Base.eltype(amv::AbstractLazyView)        =   eltype(amv.a);
Base.elsize(amv::AbstractLazyView)        =   elsize(amv.a);
Base.ndims(amv::AbstractLazyView)         =   ndims(amv.a);
Base.length(amv::AbstractLazyView)        =   length(amv.a);
Base.endof(amv::AbstractLazyView)         =   endof(amv.a);
Base.first(amv::AbstractLazyView)         =   amv[1];
Base.last(amv::AbstractLazyView)          =   getindex(amv, length(amv));
Base.start(amv::AbstractLazyView)         =   1;
Base.next(amv::AbstractLazyView, i::Int)  =   (amv[i], i+1);
Base.done(amv::AbstractLazyView, i::Int)  =   i > length(amv);

import Base.==
function =={T, N}(mv::AbstractLazyView, a::Array{T, N})
  if size(mv) != size(a)
    return false;
  else
    for i in eachindex(a)
      if mv[i] != a[i]; return false; end
    end
    return true;
  end
end
=={T, N}(a::Array{T, N}, mv::AbstractLazyView) = mv == a;

function ==(mv1::AbstractLazyView, mv2::AbstractLazyView)
  if size(mv1) != size(mv2)
    return false;
  else
    for i in eachindex(mv1.a)
      if mv1[i] != mv2[i]; return false; end
    end
    return true;
  end
end

import Base.!=
!={T, N}(mv::AbstractLazyView, a::Array{T, N}) = !(mv == a);
!={T, N}(a::Array{T, N}, mv::AbstractLazyView) = !(mv == a);
!=(mv1::AbstractLazyView, mv2::AbstractLazyView) = !(mv1==mv2);

end # module


text:
function growth_bottleneck_snr(filepath; debug=false, quiet=false)
    if !debug
        parameters = Dict{Symbol, Vector}(
            :representation => [10, 100, 1000],
            :bottleneck_representation => [10, 100, 1000],
            :seq_depth => [10, 100, 1000],
            :num_bottlenecks => collect(1:20)
        )
        num_runs = 25
    else
        parameters = Dict{Symbol, Vector}(
            :representation => [100],
            :bottleneck_representation => [100],
            :seq_depth => [100],
            :num_bottlenecks => [1, 10]
        )
        num_runs = 1
    end

    runs = grouped_param_space(GrowthScreen(), parameters, [:num_bottlenecks], num_runs)

    max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
        :inactive => (0.85, Delta(0.0)),
        :negcontrol => (0.05, Delta(0.0)),
        :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
    )

    lib = Library(max_phenotype_dists, CRISPRi())
    before = time()
    results = pmap(args -> run_exp(args[1], lib, compute_snr; run_idx=args[2]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results = DataFrame(permutedims(hcat(results...), [2, 1]))
    results[:crisprtype] = "CRISPRi"
    lib = Library(max_phenotype_dists, CRISPRn())
    before = time()
    results2 = pmap(args -> run_exp(args[1], lib, compute_snr; run_idx=args[2]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results2 = DataFrame(permutedims(hcat(results2...), [2, 1]))
    results2[:crisprtype] = "CRISPRn"
    results = vcat(results, results2)

    hierarchy = reshape([:snr, :signal, :noise], (3, 1))
    new_names = [[:technique, :score]...; fieldnames(GrowthScreen)...; :run_idx; :crisprtype;]

    results = construct_hierarchical_label(hierarchy, results, new_names)
    CSV.write(filepath, results)
end


text:
#=
# function utils
=#
export GenFunction, CallbackFunction, stringToFunction, add_genfunction, makeFunction, makeFunctions
# export @stringToFunction, @makeFunction, @makeFunctions

# Stores everything we could need about a generated function.
struct GenFunction      # example:
    name::String        # "genfunction_7"
    args::String        # "x"
    str::String         # "sin(x)"
    expr                # expression: sin(x) NOTE: could be an Expr, Symbol, Number,
    func                # handle for: genfunction_7(x) = sin(x)
end

struct CallbackFunction
    name::String        # "myFunction"
    args::Array{String} # ["x", "u", "c"] list of strings for args should match defined symbols
    body::String        # String of the function body: "a=u+2; return a;"
    func                # The function
end

# Generates a function from a string
function stringToFunction(name, args, fun)
    return eval(Meta.parse(name*"("*args*")="*fun));
end

# Adds the GenFunction to a global array of generated functions
function add_genfunction(genfun)
    global genfunc_count += 1;
    global genfunctions = [genfunctions ; genfun];
    log_entry("Generated function: "*genfun.name);
end

# Makes a GenFunction and adds it to the
# args is a string like "x,y,z"
# fun is a string like "sin(x)*y + 3*z"
function makeFunction(args, fun)
    name = "genfunction_"*string(Finch.genfunc_count+1);
    if typeof(fun) == Expr
        ex = fun;
        strfun = string(fun);
        nf = GenFunction(name, args, strfun, ex, stringToFunction(name, args, strfun));
    else
        ex = Meta.parse(fun);
        nf = GenFunction(name, args, fun, ex, stringToFunction(name, args, fun));
    end
    add_genfunction(nf);
end

# Makes either: a constant number, a genfunction, or an array of genfunctions
function makeFunctions(ex; args="x=0,y=0,z=0,t=0,node_index=1,face_index=1")
    nfuns = 0;
    if typeof(ex) <: Array
        for i=1:length(ex)
            if typeof(ex[i]) == String
                tmp = makeFunction(args, ex[i]);
                nfuns += 1;
            # else # It could be a number or function handle
            end
        end
    else
        if typeof(ex) == String
            makeFunction(args, ex);
            nfuns = 1;
        # else # It could be a number or function handle
        end
    end
    return nfuns;
end

# Interpret symbols in the arguments to a function for boundary or initial conditions.
# Replace them with something meaningful.
# Everything must be in terms of x, y, z, t, node_index, face_index
# Symbols could be:
# numbers/math ops -> no change
# x,y,z,t -> no change (provided as input arguments)
# variable -> Finch.variables[1].values[node_index]
# coefficient -> evaluate_coefficient(Finch.coefficients[1], x,y,z,t)
# indexer -> Finch.indexers[1].value
# parameter -> TODO replaced with its expression
# normal -> Finch.grid_data.faceNormals[:,face_index]
# callbackFunction -> callback_functions[1].func
function replace_symbols_in_conditions(ex)
    if typeof(ex) <: Array
        for i=1:length(ex)
            ex[i] = replace_symbols_in_conditions(ex[i]);
        end
    elseif typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = replace_symbols_in_conditions(ex.args[i]);
        end

    elseif typeof(ex) == Symbol
        # Check against the possibilities and replace as needed
        if !(ex in [:x, :y, :z, :t, :node_index, :face_index]) # don't change these
            v_index = variable_index_from_symbol(ex);
            c_index = coefficient_index_from_symbol(ex);
            i_index = indexer_index_from_symbol(ex);
            cb_index = callback_index_from_symbol(ex);
            if v_index > 0
                if variables[v_index].type == SCALAR
                    newex = :(Finch.variables[$v_index].values[node_index]);
                else
                    newex = :(Finch.variables[$v_index].values[:,node_index]);
                end
                ex = newex;
            elseif c_index > 0
                if coefficients[c_index].type == SCALAR
                    newex = :(evaluate_coefficient(Finch.coefficients[$c_index], 1, [x,y,z], t, node_index, face_index));
                else
                    num_comps = length(Finch.coefficients[c_index].value);
                    newex = :([evaluate_coefficient(Finch.coefficients[$c_index], comp, [x,y,z], t, node_index, face_index) for comp in 1:$num_comps]);
                end
                ex = newex;
            elseif i_index > 0
                ex = :(Finch.indexers[$i_index].value);
            elseif cb_index > 0
                ex = :(callback_functions[$cb_index].func);
            elseif ex === :normal
                ex = :(Finch.grid_data.facenormals[:,face_index]);
            end

        end
    else
        # other types will not be changed
    end
    return ex;
end

function variable_index_from_symbol(s)
    for v in variables
        if s === v.symbol
            return v.index;
        end
    end
    return 0;
end

function coefficient_index_from_symbol(s)
    for c in coefficients
        if s === c.symbol
            return c.index;
        end
    end
    return 0;
end

function indexer_index_from_symbol(s)
    for i=1:length(indexers)
        if s === indexers[i].symbol
            return i;
        end
    end
    return 0;
end

function callback_index_from_symbol(s)
    for i=1:length(callback_functions)
        if string(s) == callback_functions[i].name
            return i;
        end
    end
    return 0;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module BayesianDataFusion

using Distributed: Future, @spawnat, remotecall_wait, remotecall_fetch, myid
using Printf: @printf, @sprintf
using SharedArrays: SharedVector, SharedArray
using Statistics: mean

include("ROC.jl")
include("normal_wishart.jl")
include("RelationData.jl")
include("parallel_matrix.jl")
include("parallel_cg.jl")
include("sampling.jl")
include("macau.jl")

end # module


text:
using Test, GoogleSheets
using GoogleSheets: CellRangeValues, UpdateSummary
using ColorTypes
using Colors

spreadsheet_id = "1pG4OyAdePAelCT2fSBTVJ9lVYo6M-ApuTyeEPz49DOM"

if haskey(ENV, "SPREADSHEET_ID")
    spreadsheet_id = ENV["SPREADSHEET_ID"]
end

client = sheets_client(AUTH_SCOPE_READWRITE)

spreadsheet = Spreadsheet(spreadsheet_id)
sheet_name = "TestSheet"

function init_test(;add_values::Bool=true)
    if sheet_name in sheet_names(client,spreadsheet)
        delete_sheet!(client, spreadsheet, sheet_name)
    end

    add_sheet!(client, spreadsheet, sheet_name)
    @test ["Sheet1", "TestSheet"] == sheet_names(client,spreadsheet)
    @test [Sheet(client,spreadsheet,"Sheet1"),Sheet(client,spreadsheet,"TestSheet")] == sheets(client,spreadsheet)

    if(add_values)
        # Add values to the sheet
        result = update!(client, CellRange(spreadsheet, sheet_name), fill(11, 5, 5))
        @test result == UpdateSummary(CellRange(spreadsheet, "$(sheet_name)!A1:E5"), 5, 5, 25)
        result = get(client, CellRange(spreadsheet, sheet_name))
        @test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), fill("11", 5, 5), "ROWS")
    end
end

# Define struct equality.  The default equality uses === for comparisions.
Base.:(==)(x::CellRangeValues, y::CellRangeValues) = x.range == y.range && x.values == y.values && x.major_dimension == y.major_dimension

################################################################################

values = ["A" "B" "C"; "1" "2" "3"; "4" "5" "6"]
crv = CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:"), values, "TEST_DIM")
df = DataFrame(crv)
@test df == DataFrame(A= ["1", "4"], B=["2","5"], C=["3","6"])

################################################################################

init_test(add_values=false)

# Get the empty sheet
result = get(client, CellRange(spreadsheet, sheet_name))
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), nothing, "ROWS")

# Add values to the sheet
df = DataFrame(A=[1,2], B=["X","Y"])
result = update!(client, CellRange(spreadsheet, sheet_name), df)
@test result == UpdateSummary(CellRange(spreadsheet, "$(sheet_name)!A1:B3"), 2, 3, 6)
result = get(client, CellRange(spreadsheet, sheet_name))
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), ["A" "B"; "1" "X"; "2" "Y"], "ROWS")

################################################################################

init_test(add_values=false)

# Get the empty sheet
result = get(client, CellRange(spreadsheet, sheet_name))
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), nothing, "ROWS")

# Add values to the sheet
result = update!(client, CellRange(spreadsheet, sheet_name), fill(11, 5, 5))
@test result == UpdateSummary(CellRange(spreadsheet, "$(sheet_name)!A1:E5"), 5, 5, 25)
result = get(client, CellRange(spreadsheet, sheet_name))
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), fill("11", 5, 5), "ROWS")

result = get(client, CellRanges(spreadsheet, ["$(sheet_name)!A1:B2", "$(sheet_name)!C2:D5"]))
@test result == [
    CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:B2"), fill("11",2,2), "ROWS"),
    CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!C2:D5"), fill("11",4,2), "ROWS")
]

################################################################################

init_test()

# Add rows and columns to the sheet
sheet = Sheet(client, spreadsheet, sheet_name)
append!(client, sheet, 1000, 3)
result = get(client, CellRange(spreadsheet, sheet_name))
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:AC2000"), fill("11", 5, 5), "ROWS")

################################################################################

init_test()

# Insert rows
sheet = Sheet(client, spreadsheet, sheet_name)
insert_rows!(client, CellIndexRange1D(sheet, 2, 3))
result = get(client, CellRange(spreadsheet, sheet_name))
values = fill("11", 6, 5)
values[3,:] .= ""
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1001"), values, "ROWS")

# Delete rows
sheet = Sheet(client, spreadsheet, sheet_name)
delete_rows!(client, CellIndexRange1D(sheet, 2, 3))
result = get(client, CellRange(spreadsheet, sheet_name))
values = fill("11", 5, 5)
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), values, "ROWS")

################################################################################

init_test()

# Insert columns
sheet = Sheet(client, spreadsheet, sheet_name)
insert_cols!(client, CellIndexRange1D(sheet, 2, 3))
result = get(client, CellRange(spreadsheet, sheet_name))
values = fill("11", 5, 6)
values[:,3] .= ""
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:AA1000"), values, "ROWS")

# Delete columns
sheet = Sheet(client, spreadsheet, sheet_name)
delete_cols!(client, CellIndexRange1D(sheet, 2, 3))
result = get(client, CellRange(spreadsheet, sheet_name))
values = fill("11", 5, 5)
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), values, "ROWS")

################################################################################

init_test()

result = clear!(client, CellRange(spreadsheet, "$(sheet_name)!B2:C3"))
@test result == UpdateSummary(CellRange(spreadsheet, "$(sheet_name)!B2:C3"), 2, 2, 4)
result = get(client, CellRange(spreadsheet, sheet_name))
values = fill("11", 5, 5)
values[2:3,2:3] .= ""
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), values, "ROWS")

result = clear!(client, CellRange(spreadsheet, sheet_name))
@test result == UpdateSummary(CellRange(spreadsheet, "$(sheet_name)!A1:E5"), 5, 5, 25)
result = get(client, CellRange(spreadsheet, sheet_name))
@test result == CellRangeValues(CellRange(spreadsheet, "$(sheet_name)!A1:Z1000"), nothing, "ROWS")

################################################################################

init_test()

sheet = Sheet(client, spreadsheet, sheet_name)
freeze!(client, sheet, 2, 3)

m = meta(client, spreadsheet)
@test spreadsheet_id == m["spreadsheetId"]
show(client, spreadsheet)

m = meta(client, spreadsheet, sheet_name)
@test sheet_name == m["title"]
@test 2 == m["gridProperties"]["frozenRowCount"]
@test 3 == m["gridProperties"]["frozenColumnCount"]
show(client, spreadsheet, sheet_name)

sheet_id = m["sheetId"]
m = meta(client, spreadsheet, sheet_id)
@test sheet_name == m["title"]
@test 2 == m["gridProperties"]["frozenRowCount"]
@test 3 == m["gridProperties"]["frozenColumnCount"]
show(client, spreadsheet, sheet_id)

m = meta(client, sheet)
@test sheet_name == m["title"]
@test 2 == m["gridProperties"]["frozenRowCount"]
@test 3 == m["gridProperties"]["frozenColumnCount"]
show(client, sheet)

################################################################################

init_test()

sheet = Sheet(client, spreadsheet, sheet_name)

cf = CellFormat(
    background_color=RGBA(0.5, 0.5, 0.5, 0.8), number_format_type=NUMBER_FORMAT_TYPE_NUMBER, number_format_pattern="0.0",
    text_italic=true, text_bold=true, text_strikethrough=true, text_color=RGB(0.5, 0.5, 0.5), text_font_size=14)

cf2 = CellFormat(background_color=RGBA(0.5, 0.5, 0.5, 0.8), text_italic=true, text_bold=true, text_strikethrough=true)

format!(client, CellIndexRange2D(sheet, 2, 3, 1, 2), CellFormat(number_format_type=NUMBER_FORMAT_TYPE_NUMBER, number_format_pattern="0.0"))
format!(client, CellIndexRange2D(sheet, 4, 5, 1, 2), CellFormat(number_format_type=NUMBER_FORMAT_TYPE_DATE, number_format_pattern="hh:mm:ss am/pm, ddd mmm dd yyyy"))
format!(client, CellIndexRange2D(sheet, 4, 5, 1, 2), CellFormat(background_color=RGBA(0.5, 0.5, 0.5, 0.8)))
format!(client, CellIndexRange2D(sheet, 4, 5, 1, 2), CellFormat(background_color=RGB(0.5, 0.5, 0.5)))
format!(client, CellIndexRange2D(sheet, 4, 5, 1, 2), CellFormat(background_color=Gray(0.5)))
format!(client, CellIndexRange2D(sheet, 2, 3, 1, 2), cf)

format_conditional!(client, CellIndexRange2D(sheet, 2, 3, 1, 2), cf2, CONDITION_TYPE_NUMBER_GREATER_THAN_EQ, 0.5)
format_conditional!(client, CellIndexRange2D(sheet, 2, 3, 1, 2), cf2, CONDITION_TYPE_NUMBER_BETWEEN, 0.5, 1.5)

format_color_gradient!(client, CellIndexRange2D(sheet, 4, 5, 1, 2))
format_color_gradient!(client, CellIndexRange2D(sheet, 4, 5, 1, 2); min_value_type=VALUE_TYPE_NUMBER, min_value=-3, max_value_type=VALUE_TYPE_NUMBER, max_value=3, mid_color=colorant"white", mid_value_type=VALUE_TYPE_NUMBER, mid_value=0)

################################################################################


text:
function nullsim(file::String, seed::Int, reps::Int, model::String, treemethod::String, ancmethod::String)
    simfile = "nullsim_$file"
    ancfile = "ancestor_$file"
    prepare_r(seed)
    R"""
    sequences <- as.phyDat(read.dna($file, format = "fasta"))
    tree <- initialTree(sequences, $model, $treemethod)
    mlFit <- refineML(tree, sequences, $model)
    anc <- reconstructAncestor(mlFit, tree, $ancmethod)
    sims <- generateAlignments(mlFit, anc, $reps)
    rownames(sims) <- paste(names(sequences), rep(1:$reps, each = length(sequences)))
    sendToFile(sims, $simfile)
    sendToFile(as.list(as.DNAbin(anc)), $ancfile)
    """
end

function prepare_r(seed::Int)
    R"""
    library(ape)
    library(phangorn)
    set.seed($seed)

    BASES <- c('a', 'c', 'g', 't')

    initialTree <- function(data, model, method){
        distance <- dist.ml(data, model)
        tree <- switch(tolower(method),
                       nj = NJ(distance),
                       upgma = upgma(distance),
                       stop("Invalid tree building method")
                )
        return(tree)
    }

    refineML <- function(phy, seq, mod){
        mlFit <- pml(phy, seq, model = mod)
        mlFit <- optim.pml(mlFit, model = mod, control = pml.control(trace = 0))
        return(mlFit)
    }

    reconstructAncestor <- function(fitmodel, tree, method){
        anc <- ancestral.pml(fitmodel, tolower(method))
        ancMat <- anc[[getRoot(tree)]]
        winners <- apply(ancMat, 1, function(x) sample(BASES, size = 1, prob = x))
        ancestor <- winners[attr(anc, "index")]
        return(ancestor)
    }

    generateAlignment <- function(mlfit, root){
        return(as.DNAbin(simSeq(mlfit, l = length(root), rootseq = root)))
    }

    generateAlignments <- function(mlfit, root, reps){
        return(do.call(rbind, lapply(1:reps, function(x) generateAlignment(mlfit, root))))
    }

    sendToFile <- function(seqs, filename){
        write.dna(seqs, file = filename, format = "fasta")
    }
    """
end


text:
@enum PType normal fighting flying poison ground rock bug ghost steel fire water grass electric psychic ice dragon dark fairy
@enum Habitat forest grassland rainforest
@enum Diet herbivore phototroph omnivore
@enum HeightClass small medium large
@enum EggGroup monster plant

mutable struct Pokemon
	number
	name
	baseStats
	megaStats
	typeOne
	typeTwo
	basicAbilities
	highAbilities
	megaAbilities
	height
	heightClass
	weight
	weightClass
	genderRatio
	eggGroupOne
	eggGroupTwo
	avgHatchRate
	diet
	habitat
	skillList
	currentMoves
	moves
	eggMoves
	tutorMoves
	captureRate
	expDrop
end

function setupPokemon()
	setupMoves()	# TODO: Write
	setupSkills()	# TODO: Write
	setupAbilities()

	global pokemonList
	pokemonList = {}

	pokemonList["Bulbasaur"] = Pokemon()
	pokemonList["Bulbasaur"].number = 1
	pokemonList["Bulbasaur"].name = "Bulbasaur"
	pokemonList["Bulbasaur"].baseStats = Stats(5, 5, 5, 7, 7, 5)
	pokemonList["Bulbasaur"].typeOne = grass
	pokemonList["Bulbasaur"].typeTwo = poison
	pokemonList["Bulbasaur"].basicAbilities = [abilities["Overgrow"], abilities["Venom"], abilities["Chlorophyll"]]
	pokemonList["Bulbasaur"].highAbilities = [abilities["Leaf Guard"]]
	pokemonList["Bulbasaur"].height = 0.7
	pokemonList["Bulbasaur"].heightClass = small
	pokemonList["Bulbasaur"].weight = 6.9
	pokemonList["Bulbasaur"].weightClass = 1
	pokemonList["Bulbasaur"].genderRatio = 0.875
	pokemonList["Bulbasaur"].eggGroupOne = monster
	pokemonList["Bulbasaur"].eggGroupTwo = plant
	pokemonList["Bulbasaur"].avgHatchRate = 10
	pokemonList["Bulbasaur"].diet = [herbivore, phototroph]
	pokemonList["Bulbasaur"].habitat = [forest, grassland, rainforest]
	pokemonList["Bulbasaur"].skillList = {"overland": 7, "surface": 2, "jump": 2, "power": 1, "intelligence": 3, "extra": [skills["Sprouter"]]}
	pokemonList["Bulbasaur"].moves = [{"level": 1, "move": moves["Tackle"]}, {"level": 3, "move": moves["Growl"]}, {"level": 7, "move": moves["Leech Seed"]}, {"level": 9, "move": moves["Vine Whip"]}, {"level": 13, "move": moves["Poison Powder"]}, {"level": 13, "move": moves["Sleep Powder"]}, {"level": 15, "move": moves["Take Down"]}, {"level": 19, "move": moves["Razor Leaf"]}, {"level": 21, "move": moves["Sweet Scent"]}, {"level": 25, "move": moves["Growth"]}, {"level": 27, "move": moves["Double-Edge"]}, {"level": 31, "move": moves["Worry Seed"]}, {"level": 33, "move": moves["Synthesis"]}, {"level": 37, "move": moves["Seed Bomb"]}]
	pokemonList["Bulbasaur"].eggMoves = [moves["Amnesia"], moves["Charm"], moves["Curse"], moves["Endure"], moves["Giga Drain"], moves["Grass Whistle"], moves["Grassy Terrain"], moves["Ingrain"], moves["Leaf Storm"], moves["Magical Leaf"], moves["Nature Power"], moves["Petal Dance"], moves["Power Whip"], moves["Skull Bash"], moves["Sludge"]]
	pokemonList["Bulbasaur"].tutorMoves = [moves["Bind"], moves["Frenzy Plant"], moves["Giga Drain"], moves["Grass Pledge"], moves["Knock Off"], moves["Seed Bomb"]]
	pokemonList["Bulbasaur"].captureRate = 20
	pokemonList["Bulbasaur"].expDrop = 15

	pokemonList["Ivysaur"] = Pokemon()
	pokemonList["Ivysaur"].number = 2
	pokemonList["Ivysaur"].name = "Ivysaur"
	pokemonList["Ivysaur"].baseStats = Stats(6, 6, 6, 8, 8, 6)
	pokemonList["Ivysaur"].typeOne = grass
	pokemonList["Ivysaur"].typeTwo = poison
	pokemonList["Ivysaur"].basicAbilities = [abilities["Overgrow"], abilities["Venom"], abilities["Chlorophyll"]]
	pokemonList["Ivysaur"].highAbilities = [abilities["Leaf Guard"], abilities["Flower Gift"]]
	pokemonList["Ivysaur"].height = 1
	pokemonList["Ivysaur"].heightClass = medium
	pokemonList["Ivysaur"].weight = 13
	pokemonList["Ivysaur"].weightClass = 2
	pokemonList["Ivysaur"].genderRatio = 0.875
	pokemonList["Ivysaur"].eggGroupOne = monster
	pokemonList["Ivysaur"].eggGroupTwo = plant
	pokemonList["Ivysaur"].diet = [omnivore, phototroph]
	pokemonList["Ivysaur"].habitat = [forest, grassland, rainforest]
	pokemonList["Ivysaur"].skillList = {"overland": 8, "surface": 2, "jump": 2, "power": 2, "intelligence": 3, "extra": [skills["Sprouter"]]}
	pokemonList["Ivysaur"].moves = [{"level": 3, "move": moves["Growl"]}, {"level": 7, "move": move["Leech Seed"]}, {"level": 9, "move": moves["Vine Whip"]}, {"level": 13, "move": moves["Poison Powder"]}, {"level": 13, "move": moves["Sleep Powder"]}, {"level": 15, "move": moves["Take Down"]}, {"level": 20, "move": moves["Razor Leaf"]}, {"level": 23, "move": moves["Sweet Scent"]}, {"level": 28, "move": moves["Growth"]}, {"level": 31, "move": moves["Double-Edge"]}, {"level": 36, "move": moves["Worry Seed"]}, {"level": 39, "move": moves["Synthesis"]}, {"level": 44, "move": moves["Solar Beam"]}]
	pokemonList["Ivysaur"].eggMoves = []
	pokemonList["Ivysaur"].tutorMoves = [moves["Bind"], moves["Frenzy Plant"], moves["Giga Drain"], moves["Grass Pledge"], moves["Knock Off"], moves["Seed Bomb"], moves["Snore"], moves["Synthesis"], moves["Worry Seed"]]
	pokemonList["Ivysaur"].captureRate = 10
	pokemonList["Ivysaur"].expDrop = 30

	pokemonList["Venusaur"] = Pokemon()
	pokemonList["Venusaur"].number = 3
	pokemonList["Venusaur"].name = "Venusaur"
	pokemonList["Venusaur"].baseStats = Stats(8, 8, 8, 10, 10, 8)
	pokemonList["Venusaur"].megaStats = Stats(0, 2, 4, 2, 2, 0)
	pokemonList["Venusaur"].typeOne = grass
	pokemonList["Venusaur"].typeTwo = poison
	pokemonList["Venusaur"].basicAbilities = [abilities["Overgrow"], abilities["Venom"], abilities["Chlorophyll"]]
	pokemonList["Venusaur"].highAbilities = [abilities["Leaf Guard"], abilities["Flower Gift"]]
	pokemonList["Venusaur"].megaAbilities = [abilities["Thick Fat"]]
	pokemonList["Venusaur"].height = 2
	pokemonList["Venusaur"].heightClass = large
	pokemonList["Venusaur"].weight = 100
	pokemonList["Venusaur"].weightClass = 5
	pokemonList["Venusaur"].genderRatio = 0.875
	pokemonList["Venusaur"].eggGroupOne = monster
	pokemonList["Venusaur"].eggGroupTwo = plant
	pokemonList["Venusaur"].diet = [omnivore, phototroph]
	pokemonList["Venusaur"].habitat = [forest, grassland, rainforest]
	pokemonList["Venusaur"].skillList = {"overland": 9, "surface": 2, "jump": 1, "power": 4, "intelligence": 4, "extra": [skills["Sprouter"]]}
	pokemonList["Venusaur"].moves = [{"level": "e", "move": moves["Petal Dance"]}, {"level": 3, "move": moves["Growl"]}, {"level": 7, "move": moves["Leech Seed"]}, {"level": 9, "move": moves["Vine Whip"]}, {"level": 13, "move": moves["Poison Powder"]}, {"level": 13, "move": moves["Sleep Powder"]}, {"level": 15, "move": moves["Take Down"]}, {"level": 20, "move": moves["Razor Leaf"]}, {"level": 23, "move": moves["Sweet Scent"]}, {"level": 28, "move": moves["Growth"]}, {"level": 31, "move": moves["Double-Edge"]}, {"level": 32, "move": moves["Petal Dance"]}, {"level": 39, "move": moves["Worry Seed"]}, {"level": 45, "move": moves["Synthesis"]}, {"level": 53, "move": moves["Solar Beam"]}]
	pokemonList["Venusaur"].eggMoves = []
	pokemonList["Venusaur"].tutorMoves = [moves["Bind"], moves["Block"], moves["Frenzy Plant"], moves["Giga Drain"], moves["Grass Pledge"], moves["Knock Off"], moves["Outrage"], moves["Petal Dance"], moves["Seed Bomb"], moves["Snore"], moves["Synthesis"], moves["Worry Seed"]]
	pokemonList["Venusaur"].captureRate = 5
	pokemonList["Venusaur"].expDrop = 50
end

max_stars_repo_name: skeetcha/pta-battle-sim

max_issues_repo_name: skeetcha/pta-battle-sim

max_forks_repo_name: skeetcha/pta-battle-sim


text:
function reverse_number(n::Integer)
    m = 0
    while (n != 0)
        m *= 10
        m += n % 10
        n ÷= 10
    end
    return m
end

function is_palindromic(n)
    return n == reverse_number(n)
end

function main()
    max = 0

    for a in 100:1000
        for b in a:1000
            p = a * b
            if (is_palindromic(p) && max < p)
                max = p
            end
        end
    end

    println(max)
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
Base.BroadcastStyle(::Type{<:ComponentArray{T, N, A, Axes}}) where {T, N, A, Axes} = Broadcast.BroadcastStyle(A)

# Need special case here for adjoint vectors in order to avoid type instability in axistype
Broadcast.combine_axes(a::ComponentArray, b::AdjOrTransComponentVector) = (axes(a)[1], axes(b)[2])
Broadcast.combine_axes(a::AdjOrTransComponentVector, b::ComponentArray) = (axes(b)[2], axes(a)[1])

Broadcast.axistype(a::CombinedAxis, b::AbstractUnitRange) = a
Broadcast.axistype(a::AbstractUnitRange, b::CombinedAxis) = b
Broadcast.axistype(a::CombinedAxis, b::CombinedAxis) = CombinedAxis(FlatAxis(), Base.Broadcast.axistype(_array_axis(a), _array_axis(b)))
Broadcast.axistype(a::T, b::T) where {T<:CombinedAxis} = a

Base.promote_shape(a::Tuple{Vararg{CombinedAxis}}, b::NTuple{N,AbstractUnitRange}) where {N} = Base.promote_shape(_array_axis.(a), b)
Base.promote_shape(a::NTuple{N,AbstractUnitRange}, b::Tuple{Vararg{CombinedAxis}}) where {N} = Base.promote_shape(a, _array_axis.(b))
Base.promote_shape(a::Tuple{Vararg{CombinedAxis}}, b::Tuple{Vararg{CombinedAxis}}) = Base.promote_shape(_array_axis.(a), _array_axis.(b))
Base.promote_shape(a::T, b::T) where {T<:Tuple{Vararg{CombinedAxis}}} = a

# From https://github.com/JuliaArrays/OffsetArrays.jl/blob/master/src/OffsetArrays.jl
Base.dataids(A::ComponentArray) = Base.dataids(parent(A))
Broadcast.broadcast_unalias(dest::ComponentArray, src) = getdata(dest) === getdata(src) ? src : Broadcast.unalias(dest, src)

max_forks_repo_name: jonniedie/ComponentArrays


text:
export rsvd_fnkz

immutable OuterProduct{T}
    X :: Matrix{T}
    Y :: Matrix{T}
end

X ⊗ Y = OuterProduct{promote_type(eltype(X), eltype(Y))}(X, Y)

@doc """
Compute the randomized SVD by iterative refinement from randomly selected
columns/rows.

TODO: Currently works only on tall and skinny A.

Arguments:

    A: Matrix whose SVD is desired
    k: Desired rank of approximation. Must be k ≤ min(m, n).
    l: Number of columns/rows to sample at each iteration.
       Must be 1 ≤ l ≤ k. Default: k
    N: Maximum number of iterations. Default: min(size(A))
    ϵ: Relative threshold for convergence, as measured by growth of the spectral norm.
       Default: prod(size(A))*eps(real(float(one(eltype(A)))))
    method: :eig - Solve eigenproblem (Default)
            :svd - Solve singular problem
    verbose: If true, prints convergence information at each iteration. Default: false

Returns:

    S: An SVD object of rank ≤ k.

Reference:

    @inproceedings{,
        author={Friedland, S. and Niknejad, A. and Kaveh, Mostafa and Zare, H.},
        booktitle={System of Systems Engineering, 2006 IEEE/SMC International Conference on},
        title={Fast Monte-Carlo low rank approximations for matrices},
        year={2006},
        month={April},
        pages={218--223},
        doi={10.1109/SYSOSE.2006.1652299}
    }
""" ->
function rsvd_fnkz(A, k::Int;
    l::Int=k, N::Int=minimum(size(A)), verbose::Bool=false,
    method::Symbol=:eig,
    ϵ::Real=prod(size(A))*eps(real(float(one(eltype(A))))))

    const m, n = size(A)
    const dosvd = method == :svd
    @assert 1 ≤ l ≤ k ≤ min(m, n)

    #Initialize k-rank approximation B₀ according to (III.9)
    const tallandskinny = m ≥ n
    B₀ = if tallandskinny
        #Select k columns of A
        π = randperm(n)[1:k]
        A[:, π]
    else
        warn("Not implemented")
        #Select k rows of A
        π = randperm(m)[1:k]
        A[π, :]
    end
    X, RRR = qr(B₀)
    X = X[:, abs(diag(RRR)) .> ϵ] #Remove linear dependent columns
    B = tallandskinny ? X ⊗ A'X : X ⊗ X'A

    oldnrmB = 0.0
    Λ = 0
    X̃ = 0
    V = 0
    for t=1:N
        π = randperm(k)[1:l]
        #Update B using Theorem 2.4
        X, RRR = qr([B.X A[:, π]]) #We are doing more work than needed here
        X = X[:, abs(diag(RRR)) .> ϵ] #Remove linearly dependent columns
        Y = A'X
        if dosvd
            S = svdfact!(Y)
            Λ = S[:S].^2
            B = (X * S[:V]) ⊗ (S[:U]*Diagonal(S[:S]))
        else
            E = eigfact!(Symmetric(Y'Y))
            π = sortperm(E[:values], rev=true)[1:min(length(E[:values]), k)]
            Λ = E[:values][π]
            O = E[:vectors][:, π] #Eq. 2.6
            X̃ = X * O
            B = X̃ ⊗ A'X̃
        end

        length(Λ)==0 && break

        nrmB = √Λ[1]
        verbose && println("iteration $t: Norm: $nrmB")
        oldnrmB > (1-ϵ)*nrmB && break
        oldnrmB = nrmB
    end
    for i=1:size(B.Y, 2)
        scale!(slice(B.Y, :, i), 1/√Λ[i])
    end
    Base.LinAlg.SVD(B.X, √Λ, B.Y')
end


text:
module TestDay12

import AOC2021.Day12.part1, AOC2021.Day12.part2
using Test

test_input_short = """
dc-end
HN-start
start-kj
dc-start
dc-HN
LN-dc
HN-end
kj-sa
kj-HN
kj-dc"""

test_input_long = """
fs-end
he-DX
fs-he
start-DX
pj-DX
end-zg
zg-sl
zg-pj
pj-he
RW-he
fs-DX
pj-RW
zg-RW
start-pj
he-WI
zg-he
pj-fs
start-RW"""

@testset "Day 12" begin
    @testset "part 1" begin
        @test part1(test_input_short) == 19
        @test part1(test_input_long) == 226
    end
    @testset "part 2" begin
        @test part2(test_input_short) == 103
        @test part2(test_input_long) == 3509
    end
end

end # module


text:
## Contains functions for animating Delta Robot solution
using Plots

"""
transforms the solution's state variables for the delta robot to cartesian
coordinates
"""
function transform_to_cartesian(t, phi1, phi2, phi3, p)
    # make sure lengths of all vectors are the same
    @assert(length(t) == length(phi1))
    @assert(length(t) == length(phi2))
    @assert(length(t) == length(phi3))

    # get parameters
    l1 = p.l1;
    l2 = p.l2;
    lc1 = l1 / 2.0;
    lc2 = l2 / 2.0;

    # end of link 1
    x1 = zeros(length(t))
    y1 = l1 * cos.(phi1) .+ p.rb
    z1 = l1 * sin.(phi1)

    # end of link 2 (end-effector)
    x2 = l2 .* sin.(phi2) .* sin.(phi3)
    y2 = y1 + l2 * cos.(phi2)
    z2 = z1 + l2 * sin.(phi2) .* cos.(phi3)

    return x1, y1, z1, x2, y2, z2
end

"""
generates an animation object from the solution to the delta robot arm IVP and
delta robot parameters

todo: turn into plot recipe and remove dependency on plots and gr()
"""
function animate_arm(sol, p)
    gr()
    t = sol.t
    x1, y1, z1, x2, y2, z2 = transform_to_cartesian(t, sol[1,:], sol[2,:], sol[3,:], p)

    anim = Animation()
    for i in 1:length(t)
        t_str = string("Time = ", round(t[i]), " sec");
        p1 = plot3d([0,-y1[i]], [0,-x1[i]], [0, -z1[i]], size=(200,200), xlim=(-1,1), ylim=(-1,1), zlim=(-1, 1), markersize = 5, markershape=:circle, camera = (45, 60))
        plot3d!(p1, [-y1[i],-y2[i]], [-x1[i],-x2[i]], [-z1[i], -z2[i]], markersize=5, markershape = :circle, title = t_str, title_location=:left)

        p2 = plot3d([0,-y1[i]], [0,-x1[i]], [0, -z1[i]], size=(200,200), xlim=(-1,1), ylim=(-1,1), zlim=(-1, 1), markersize = 5, markershape=:circle, camera = (0, 0))
        plot3d!(p2, [-y1[i],-y2[i]], [-x1[i],-x2[i]], [-z1[i], -z2[i]], markersize=5, markershape = :circle, title="Side View")

        p3 = plot3d([0,-y1[i]], [0,-x1[i]], [0, -z1[i]], size=(200,200), xlim=(-1,1), ylim=(-1,1), zlim=(-1, 1), markersize = 5, markershape=:circle, camera = (0, 90))
        plot3d!(p3, [-y1[i],-y2[i]], [-x1[i],-x2[i]], [-z1[i], -z2[i]], markersize=5, markershape = :circle, title="Bottom View")

        p = plot(p1, p2, p3, layout=(3,1), size=(400, 1200), legend=:none)
        frame(anim, p)
    end
    return anim
end

"""
generates an animation object from the solution to the delta robot IVP and
delta robot parameters

todo: turn into plot recipe and remove dependency on plots and gr()
"""
function animate_robot(sol, p)
    gr()
    t = sol.t;
    x11, y11, z11, x12, y12, z12 = transform_to_cartesian(sol.t, sol[1,:], sol[2,:], sol[3,:], p)
    x21, y21, z21, x22, y22, z22 = transform_to_cartesian(sol.t, sol[4,:], sol[5,:], sol[6,:], p)
    x31, y31, z31, x32, y32, z32 = transform_to_cartesian(sol.t, sol[7,:], sol[8,:], sol[9,:], p)

    Rz(θ) = [cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1];
    x21, y21, z21 = Rz(-2 * pi / 3) * [x21, y21, z21];
    x22, y22, z22 = Rz(-2 * pi / 3) * [x22, y22, z22];
    x31, y31, z31 = Rz(2 * pi / 3) * [x31, y31, z31];
    x32, y32, z32 = Rz(2 * pi / 3) * [x32, y32, z32];

    anim = Animation()
    for i in 1:length(sol.t)
        t_str = string("Time = ", round(sol.t[i]), " sec");
        p1 = plot3d([0,-x11[i]], [-p.rb,-y11[i]], [0.0, -z11[i]], size=(400,300), xlim=(-1,1), ylim=(-1,1), zlim=(-1, 1), markersize = 5, markershape=:circle, camera = (45, 60), legend=:none, linecolor=:blue, markercolor=:blue);
        plot3d!([-x11[i],-x12[i]], [-y11[i],-y12[i]], [-z11[i], -z12[i]], markersize=5, markershape = :circle, title = t_str, title_location=:left, aspect_ratio=:auto, linecolor=:blue, markercolor=:blue);
        plot3d!([-p.rb * -sin(-2 * pi / 3),-x21[i]], [-p.rb * cos(-2 * pi / 3),-y21[i]], [0, -z21[i]], markersize = 5, markershape=:circle, linecolor=:red, markercolor=:red);
        plot3d!([-x21[i],-x22[i]], [-y21[i],-y22[i]], [-z21[i], -z22[i]], markersize=5, markershape = :circle, linecolor=:red, markercolor=:red);
        plot3d!([-p.rb * -sin(2 * pi / 3),-x31[i]], [-p.rb * cos(2 * pi / 3),-y31[i]], [0, -z31[i]], markersize = 5, markershape=:circle, linecolor=:green, markercolor=:green);
        plot3d!([-x31[i],-x32[i]], [-y31[i],-y32[i]], [-z31[i], -z32[i]], markersize=5, markershape = :circle, linecolor=:green, markercolor=:green);

        p2 = plot3d([0,-x11[i]], [-p.rb,-y11[i]], [0.0, -z11[i]], size=(400,300), xlim=(-1,1), ylim=(-1,1), zlim=(-1, 1), markersize = 5, markershape=:circle, camera = (0, 0), legend=:none, linecolor=:blue, markercolor=:blue);
        plot3d!([-x11[i],-x12[i]], [-y11[i],-y12[i]], [-z11[i], -z12[i]], markersize=5, markershape = :circle, aspect_ratio=:auto, linecolor=:blue, markercolor=:blue);
        plot3d!([-p.rb * -sin(-2 * pi / 3),-x21[i]], [-p.rb * cos(-2 * pi / 3),-y21[i]], [0, -z21[i]], markersize = 5, markershape=:circle, linecolor=:red, markercolor=:red);
        plot3d!([-x21[i],-x22[i]], [-y21[i],-y22[i]], [-z21[i], -z22[i]], markersize=5, markershape = :circle, linecolor=:red, markercolor=:red);
        plot3d!([-p.rb * -sin(2 * pi / 3),-x31[i]], [-p.rb * cos(2 * pi / 3),-y31[i]], [0, -z31[i]], markersize = 5, markershape=:circle, linecolor=:green, markercolor=:green);
        plot3d!([-x31[i],-x32[i]], [-y31[i],-y32[i]], [-z31[i], -z32[i]], markersize=5, markershape = :circle, linecolor=:green, markercolor=:green, title="Side View");

        p3 = plot3d([0,-x11[i]], [-p.rb,-y11[i]], [0.0, -z11[i]], size=(400,300), xlim=(-1,1), ylim=(-1,1), zlim=(-1, 1), markersize = 5, markershape=:circle, camera = (0, 90), legend=:none, linecolor=:blue, markercolor=:blue);
        plot3d!([-x11[i],-x12[i]], [-y11[i],-y12[i]], [-z11[i], -z12[i]], markersize=5, markershape = :circle, aspect_ratio=:auto, linecolor=:blue, markercolor=:blue);
        plot3d!([-p.rb * -sin(-2 * pi / 3),-x21[i]], [-p.rb * cos(-2 * pi / 3),-y21[i]], [0, -z21[i]], markersize = 5, markershape=:circle, linecolor=:red, markercolor=:red);
        plot3d!([-x21[i],-x22[i]], [-y21[i],-y22[i]], [-z21[i], -z22[i]], markersize=5, markershape = :circle, linecolor=:red, markercolor=:red);
        plot3d!([-p.rb * -sin(2 * pi / 3),-x31[i]], [-p.rb * cos(2 * pi / 3),-y31[i]], [0, -z31[i]], markersize = 5, markershape=:circle, linecolor=:green, markercolor=:green);
        plot3d!([-x31[i],-x32[i]], [-y31[i],-y32[i]], [-z31[i], -z32[i]], markersize=5, markershape = :circle, linecolor=:green, markercolor=:green, title="Bottom View");

        plt = plot(p1, p2, p3, layout=(3,1), size=(400, 1200), legend=:none)
        frame(anim, plt)
    end

    return anim
end

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
push!(LOAD_PATH, "/Users/czhang/Desktop/Projects/dw_/julialib/")
import DimmWitted
DimmWitted.set_libpath("/Users/czhang/Desktop/Projects/dw_/libdw_julia")

######################################
# The following piece of code creates a
# synthetic data set:
#    - Data type is Cdouble
#    - Modle type is Array{Cdouble}
#
nexp = 100000
nfeat = 1024
examples = Array(Cdouble, nexp, nfeat+1)
for row = 1:nexp
	for col = 1:nfeat
		examples[row, col] = 1
	end
	if rand() > 0.8
		examples[row, nfeat+1] = 0
	else
		examples[row, nfeat+1] = 1
	end
end
model = Cdouble[0 for i = 1:nfeat]

######################################
# Define the loss function and gradient
# function for logistic regression
#
function loss(row::Array{Cdouble,1}, model::Array{Cdouble,1})
	const label = row[length(row)]
	const nfeat = length(model)
	d = 0.0
	for i = 1:nfeat
		d = d + row[i]*model[i]
	end
	return (-label * d + log(exp(d) + 1.0))
end

function grad(row::Array{Cdouble,1}, model::Array{Cdouble,1})
	const label = row[length(row)]
	const nfeat = length(model)
	d = 0.0
	for i = 1:nfeat
		d = d + row[i]*model[i]
	end
	d = exp(-d)
		Z = 0.00001 * (-label + 1.0/(1.0+d))
  	for i = 1:nfeat
  		model[i] = model[i] - row[i] * Z
  	end
	return 1.0
end

######################################
# Create a DimmWitted object using data
# and model. You do not need to specify
# the type, they are infer'ed by the
# open() function, which is parametric.
#
dw = DimmWitted.open(examples, model,
                DimmWitted.MR_PERMACHINE,
                DimmWitted.DR_SHARDING,
                DimmWitted.AC_ROW)

######################################
# Register functions.
#
handle_loss = DimmWitted.register_row(dw, loss)
handle_grad = DimmWitted.register_row(dw, grad)

######################################
# Run 10 epoches.
#
for iepoch = 1:10
	rs = DimmWitted.exec(dw, handle_loss)
	println("LOSS: ", rs/nexp)
	rs = DimmWitted.exec(dw, handle_grad)
end

max_stars_repo_name: GloriaGo/dw

max_issues_repo_name: GloriaGo/dw

max_forks_repo_name: GloriaGo/dw


text:
using NestedTuples
using Documenter

makedocs(;
    modules=[NestedTuples],
    authors="Chad Scherrer <chad.scherrer@gmail.com> and contributors",
    repo="https://github.com/cscherrer/NestedTuples.jl/blob/{commit}{path}#L{line}",
    sitename="NestedTuples.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://cscherrer.github.io/NestedTuples.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/cscherrer/NestedTuples.jl",
)


text:
export State, Trajectory
export primary_body, secondary_body
export collision, check_distance, crashed

#------------------#
# ORBITAL PROBLEMS #
#------------------#

struct State{M<:Abstract_DynamicalModel,F<:Abstract_ReferenceFrame,uType,tType,isinplace,O<:SciMLBase.AbstractODEProblem{uType,tType,isinplace}} <: SciMLBase.AbstractODEProblem{uType,tType,isinplace}
    model :: M
    frame :: F  # Reference frame that the problem's u0 is defined in
    prob :: O
end
State(model::Abstract_DynamicalModel, reference_frame::Abstract_ReferenceFrame, u0::AbstractArray, tspan) =
    State(model, reference_frame, MArray{Tuple{size(u0)...}}(u0), tspan)
State(model::Abstract_DynamicalModel, reference_frame::Abstract_ReferenceFrame, u0::StaticArray, tspan) =
    State(model, reference_frame, ODEProblem(model, u0, tspan, parameters(model)))
State(model::Abstract_DynamicalModel, u0::AbstractArray, tspan) = State(model, default_reference_frame(model), u0, tspan)

struct Trajectory{M<:Abstract_DynamicalModel,F<:Abstract_ReferenceFrame,T,N,A,O<:DiffEqBase.AbstractTimeseriesSolution{T,N,A},} <: DiffEqBase.AbstractTimeseriesSolution{T,N,A}
    model :: M
    frame :: F  # Reference frame that the solution is defined in
    sol :: O
end

State(traj::Trajectory) = State(traj.model, traj.frame, traj.sol.prob)

primary_body(state::State) = primary_body(state.model)
primary_body(traj::Trajectory) = primary_body(traj.model)
secondary_body(state::State) = secondary_body(state.model)
secondary_body(traj::Trajectory) = secondary_body(traj.model)

DiffEqBase.remake(state::State; kwargs...) = State(state.model, state.frame, remake(state.prob; kwargs...))
ModelingToolkit.parameters(state::State) = ModelingToolkit.parameters(state.model)

#---------------#
# INTERPOLATION #
#---------------#

# Interpolation
(traj::Trajectory)(t::Number) = State(traj.model, traj.frame, traj.sol(t), (t, t))
function (traj::Trajectory)(t::AbstractArray{<:Number})
    new_sol = DiffEqBase.build_solution(traj.prob, traj.alg, t, traj.sol(t); interp=traj.interp, retcode=traj.retcode)
    return Trajectory(traj.model, traj.frame, new_sol)
end

# Indexing
Base.getindex(state::State, idx...) = getindex(state.prob.u0, idx...)
Base.getindex(traj::Trajectory, idx...) = State(traj.model, traj.frame, getindex(traj.sol, idx...), (traj.sol.t[idx...], traj.sol.t[idx...]))
Base.getindex(traj::Trajectory, idx::Int) = State(traj.model, traj.frame, getindex(traj.sol, idx), (traj.sol.t[idx], traj.sol.t[idx]))
Base.getindex(traj::Trajectory, idx::AbstractArray{Int}) = State(traj.model, traj.frame, getindex(traj.sol, idx), (traj.sol.t[idx], traj.sol.t[idx]))
Base.axes(state::State, idx...) = axes(state.prob, idx...)
Base.axes(traj::Trajectory, idx...) = axes(traj.sol, idx...)
Base.firstindex(traj::Trajectory) = firstindex(traj.sol)
Base.firstindex(traj::Trajectory, idx) = firstindex(traj.sol, idx)
Base.lastindex(traj::Trajectory) = lastindex(traj.sol)
Base.lastindex(traj::Trajectory, idx) = lastindex(traj.sol, idx)
Base.size(traj::Trajectory) = size(traj.sol)

function Base.getproperty(x::T, b::Symbol) where {T<:Union{State,Trajectory}}
    if hasfield(T, b)
        return getfield(x, b)
    else
        return getproperty(isa(x, State) ? x.prob : x.sol, b)
    end
end

#---------#
# DISPLAY #
#---------#

Base.show(io::IO, A::State) = println(io, summary(A))
Base.summary(state::State) = string(
    SciMLBase.TYPE_COLOR, nameof(typeof(state)), SciMLBase.NO_COLOR, " in ",
    SciMLBase.TYPE_COLOR, state.model, SciMLBase.NO_COLOR, " in ",
    SciMLBase.TYPE_COLOR, state.frame, "\n",
    "    ", SciMLBase.NO_COLOR, "Underlying ", summary(state.prob), "\n",
    "    ", SciMLBase.NO_COLOR, "tspan = ", state.prob.tspan, "\n",
    "    u0    = ", state.prob.u0)
Base.show(io::IO, m::MIME"text/plain", A::Trajectory) = show(io, A) # XXX: Required because also defined in DiffEqBase
function Base.show(io::IO, A::Trajectory)
    println(io, string(
        SciMLBase.TYPE_COLOR, nameof(typeof(A)), SciMLBase.NO_COLOR, " in ",
        SciMLBase.TYPE_COLOR, A.model, SciMLBase.NO_COLOR, " in ",
        SciMLBase.TYPE_COLOR, A.frame, SciMLBase.NO_COLOR))
    println(io, string(
        "    retcode  = $(A.retcode) [$(length(A.t)) timesteps]\n",
        "    t        = ($(A.t[begin]), $(A.t[end]))\n",
        "    u[begin] = $(A.u[begin])\n",
        "    u[end]   = $(A.u[end])\n"
    ))
end

#-------#
# SOLVE #
#-------#

const DEFAULT_ALG = Vern7();

# XXX: Required to support solving a State problem.
DiffEqBase.solve(state::State, args...; reltol=1e-10, abstol=1e-10, kwargs...) =
    DiffEqBase.__solve(state, args...; reltol, abstol, kwargs...)

DiffEqBase.__solve(state::State; kwargs...) = DiffEqBase.__solve(state, DEFAULT_ALG; kwargs...)

# The __solve() method does the actual heavy lifting, including converting to a Trajectory.
function DiffEqBase.__solve(state::State, alg::OrdinaryDiffEqAlgorithm; userdata=Dict(), callback=nothing, kwargs...)
    default_frame = default_reference_frame(state.model)
    real_state = convert_to_frame(state, default_frame)

    # Pass the default state into the underlying solver
    # TODO: Remove the need for this in DiffCorrectAxisymmetric
    userdata = deepcopy(userdata)
    userdata[:real_state] = real_state

    # Copy the callbacks (for thread-safety)
    callback = deepcopy(callback)

    # Call the underlying solver
    raw_sol = solve(real_state.prob, alg; userdata, callback, kwargs...)
    return Trajectory(state.model, default_frame, raw_sol)
end


text:
function derive_pulse(average_power, rep_rate, pulse_time, pulse_type::Symbol=:gaussian, chirp=0)
    pulse_energy = sqrt(average_power / rep_rate / pulse_time)
    #u0(x) = (1+0im) * pulse_energy * exp(x^2/pulse_time^2)
    u0(x) = (1+0im) * pulse_energy * sech(x/pulse_time)
    return u0
end

function derive_pulse(peak_power, pulse_FWHM)
    sqrtpower = sqrt(peak_power)
    pulse_0 = pulse_FWHM/(2*log(1+sqrt(2)))
    u0(x) = (1+0im)*sqrtpower*sech(x/pulse_0)
end


text:
struct Fluct
    cons::Dictionary{Int, Float64}
    vars::Dictionary{Int, Float64}
    Fluct(cons, vars) = new(zeros(cons), zeros(vars))
end

function reset!(fluct)
    zeros! = d -> foreach(k -> set!(d, k, 0.0), keys(d))
    zeros!(fluct.cons)
    zeros!(fluct.vars)
end

function copy_to!(fluct, cons, vars)
    foreach(k -> set!(fluct.cons, k, cons[k]), keys(cons))
    foreach(k -> set!(fluct.vars, k, vars[k]), keys(vars))
end

function copy_from!(fluct, cons, vars)
    foreach(k -> set!(cons, k, fluct.cons[k]), keys(cons))
    foreach(k -> set!(vars, k, fluct.vars[k]), keys(vars))
end


text:
using DispatchedTuples
using BenchmarkTools

# (N, N_buf) are low for quick doc builds
# (N, N_buf) = (5, 3) is reasonable
N = 0
N_buf = 1

foo(i) = Symbol(:Foo, i)
expr(i) = :(struct $(foo(i)) end)
for i in 1:N+N_buf
    eval(expr(i))
end

fooinst(i::Int) = eval(:($(foo(i))))()

footup(n::Int) = ntuple(n) do i
    (fooinst(i), i)
end

function perf_dt(N, N_buf)
    println("********* DispatchedTuple")
    for n in N_buf:N+N_buf
        println("--- n = $n")
        tup = footup(n)
        @btime footup($n)

        dtup = DispatchedTuple(tup)
        @btime DispatchedTuple($tup)
        # @btime $dtup[Foo3()] # always results in 0.028 ns (0 allocations: 0 bytes)
    end
end

function perf_ds(N, N_buf)
    println("********* DispatchedSet")
    for n in N_buf:N+N_buf
        println("--- n = $n")
        tup = footup(n)
        @btime footup($n)

        dtup = DispatchedSet(tup)
        @btime DispatchedSet($tup)
        # @btime $dtup[Foo3()] # always results in 0.028 ns (0 allocations: 0 bytes)
    end
end

perf_dt(N, N_buf)
perf_ds(N, N_buf)


text:
function test_marginal_inference_complete_evidence()
	s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s)
    x = [true, false]
    @test marginalInference!(spn, x) ≈ log(0.522)
end

function test_marginal_inference_incomplete_evidence()
	s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s)
    x = [true, false]
    e = BitVector([true, false])
    @test marginalInference!(spn, x, e) ≈ log(0.6*0.5 + 0.6*0.2 + 0.9*0.3)

    # no evidence at all should lead to log(1)=0
    e = BitVector([false, false])
    p = marginalInference!(spn, x, e)
    @test p ≈ 0 atol=10e-6
end

function test_conditional_inference()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s)
    x = [true, false]

    # empty query and evidence should lead to log(1)=0
    q = e = falses(2)
    p = conditionalInference!(spn, x, q, e)
    @test p ≈ 0 atol=10e-6

    # if the evidence contains variables of the query an error is thrown

    # compute P(X_2 = false | X_1 = true):
    x = [true, false]
    e = BitVector([true,false])
    q = BitVector([false,true])
    # Expected result:
    # P(X_1=t,x_2=f) = 0.522
    # P(X_1=t) = 0.6*0.5 + 0.6*0.2 + 0.9*0.3 (see above)
    expected = log(0.522 / (0.6*0.5 + 0.6*0.2 + 0.9*0.3))
    @test conditionalInference!(spn, x, q, e) ≈ expected
end

test_marginal_inference_complete_evidence()
test_marginal_inference_incomplete_evidence()
test_conditional_inference()


text:
module WeightedEnsemble

using LinearAlgebra
using SparseArrays
using Printf
using NearestNeighbors
using Distributed
using SharedArrays
using StatsBase
using Distributions
using Arpack

# Load types
include("types.jl")
# Load data structures
include("ensemble.jl")
include("bins.jl")
# Load resampling algorithms
include("resampling.jl")
# Load allocation algorithms
include("allocation.jl")
export trivial_allocation!, minimal_bin_allocation!,
    targeted_bin_allocation!, optimal_bin_allocation!,
    uniform_bin_allocation!, within_bin_allocation!
# Load selection algorithms
include("selection.jl")
export trivial_selection!, uniform_selection!,
    optimal_selection!, targeted_selection!, repopulate!
# Load coarse model algorithms
include("coarse_serial.jl")
include("coarse_parallel.jl")
include("coarse_threads.jl")
# Load utility functions
include("utils.jl")
export setup_Voronoi_bins, Dirac_to_Ensemble
# Load sampler structures
include("sampler.jl")
# Load WE methods
include("we_serial.jl")
export run_we, run_we!, run_we_observables
include("we_parallel.jl")
export prun_we, prun_we!, prun_we_observables
include("we_threads.jl")
export trun_we, trun_we!, trun_we_observables

export Ensemble, Bins, WEsampler, DistributedWEsampler

end # end module

max_stars_repo_name: gideonsimpson/JuWeightedEnsemble

max_issues_repo_name: gideonsimpson/JuWeightedEnsemble

max_forks_repo_name: gideonsimpson/JuWeightedEnsemble


text:
using DataFrames

mutable struct problem
    data::Array{Float64}
    Sigma::Array{Float64}
end

################ Pitprops
thepitprops=[[1,0.954,0.364,0.342,-0.129,0.313,0.496,0.424,0.592,0.545,0.084,-0.019,0.134];
       [0.954,1,0.297,0.284,-0.118,0.291,0.503,0.419,0.648,0.569,0.076,-0.036,0.144];
       [0.364,0.297,1,0.882,-0.148,0.153,-0.029,-0.054,0.125,-0.081,0.162,0.22,0.126];
       [0.342,0.284,0.882,1,0.22,0.381,0.174,-0.059,0.137,-0.014,0.097,0.169,0.015];
       [-0.129,-0.118,-0.148,0.22,1,0.364,0.296,0.004,-0.039,0.037,-0.091,-0.145,-0.208];
       [0.313,0.291,0.153,0.381,0.364,1,0.813,0.09,0.211,0.274,-0.036,0.024,-0.329];
       [0.496,0.503,-0.029,0.174,0.296,0.813,1,0.372,0.465,0.679,-0.113,-0.232,-0.424];
       [0.424,0.419,-0.054,-0.059,0.004,0.09,0.372,1,0.482,0.557,0.061,-0.357,-0.202];
       [0.592,0.648,0.125,0.137,-0.039,0.211,0.465,0.482,1,0.526,0.085,-0.127,-0.076];
       [0.545,0.569,-0.081,-0.014,0.037,0.274,0.679,0.557,0.526,1,-0.319,-0.368,-0.291];
       [0.084,0.076,0.162,0.097,-0.091,-0.036,-0.113,0.061,0.085,-0.319,1,0.029,0.007];
       [-0.019,-0.036,0.22,0.169,-0.145,0.024,-0.232,-0.357,-0.127,-0.368,0.029,1,0.184];
       [0.134,0.144,0.126,0.015,-0.208,-0.329,-0.424,-0.202,-0.076,-0.291,0.007,0.184,1];]
thepitprops=reshape(thepitprops, (13,13));
B=sqrt(thepitprops);
pitprops=problem(B, thepitprops)

max_stars_repo_name: ryancorywright/Optimal-SPCA

max_issues_repo_name: ryancorywright/Optimal-SPCA

max_forks_repo_name: ryancorywright/Optimal-SPCA


text:
__precompile__()
module ConvUtils

using Images
using MosaicViews
using MappedArrays
using MLDataUtils

const DEFAULT_COLORNEG = RGB(0.230, 0.299, 0.754)
const DEFAULT_COLORZRO = RGB(0.95,  0.95,  0.95)
const DEFAULT_COLORPOS = RGB(0.706, 0.016, 0.150)

export

    falsecolor,
    upscale,
    immosaic

function falsecolor(A::Array{T,N};
                    symmetric = true,
                    colorneg  = DEFAULT_COLORNEG,
                    colorzero = DEFAULT_COLORZRO,
                    colorpos  = DEFAULT_COLORPOS) where {T<:Number,N}
    mi, ma = extrema(A)
    val = max(abs(mi), ma)
    ss = if symmetric
        scalesigned(-val, T(0), val)
    else
        scalesigned(mi, T(0), ma)
    end
    x = ss(one(T))
    cs = colorsigned(map(c->convert(RGB{typeof(x)},c), (colorneg, colorzero, colorpos))...)
    mappedarray(c->cs(ss(c)), A)
end

falsecolor(A; kw...) = falsecolor(Array(A); kw...)

function upscale(A::AbstractArray{T,N}, val::Int) where {T,N}
    repeat(A, inner=(val, val, ntuple(i->1, N-2)...))
end

function immosaic(A::AbstractArray, fill=RGB(1.,1.,1.); rep=1, npad=5, ncol=-1, nrow=-1)
    if ndims(A) == 3 && ncol == -1 && nrow == -1
        ncol = max(ceil(Int,sqrt(size(A,3))), 5)
    end
    mosaicview(upscale(A, rep), fill; npad=npad, nrow=nrow, ncol=ncol, rowmajor=true)
end

end # module


text:
@testset "Triangular" begin
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A, VA, L = randn.(Ref(rng), [N, N, N], [N, N, N])
            @test check_errs(LowerTriangular, LowerTriangular(L), A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A, VA, U = randn.(Ref(rng), [N, N, N], [N, N, N])
            @test check_errs(UpperTriangular, UpperTriangular(U), A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = LowerTriangular(randn(rng, N, N))
            VA = LowerTriangular(randn(rng, N, N))
            @test check_errs(det, 10.0, A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 3
        for _ in 1:10
            A = UpperTriangular(randn(rng, N, N))
            VA = UpperTriangular(randn(rng, N, N))
            @test check_errs(det, 10.0, A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = LowerTriangular(exp.(randn(rng, N, N)))
            VA = LowerTriangular(randn(rng, N, N))
            @test check_errs(logdet, 10.0, A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 3
        for _ in 1:10
            A = UpperTriangular(exp.(randn(rng, N, N)))
            VA = UpperTriangular(randn(rng, N, N))
            @test check_errs(logdet, 10.0, A, VA)
        end
    end

    # Check that the optimisations occur correctly and produce the required types when
    # everything is Diagonal.
    let rng = MersenneTwister(123456)
        A = UpperTriangular(exp.(randn(rng, 10, 10)))

        @test ∇(det)(A)[1] isa Diagonal
        @test ∇(A->det(A) + det(A))(A)[1] isa Diagonal
        @test ∇(logdet)(A)[1] isa Diagonal
        @test ∇(A->logdet(A) + det(A))(A)[1] isa Diagonal
    end
end