text:
abstract type AbstractPath end

# Required methods
"""
	curves(P::AbstractPath)
Return an array of the curves that make up the path `P`.
"""
curves(p::AbstractPath) = @error "No curves() method defined for type $(typeof(p))"

# Default implementations
"""
	curve(P::AbstractPath,k::Integer)
Return the `k`th curve in the path `P`.
"""
curve(p::AbstractPath,k::Integer) = curves(p)[k]

"""
	vertex(P::AbstractPath,k::Integer)
Return the `k`th vertex of the path `P`.
"""
function vertex(P::AbstractPath,k::Integer)
	C = curves(P)
	n = length(C)
	if 1 ≤ k ≤ n
		point(C[k],0)
	elseif k==n+1
		point(C[n],1)
	else
		throw(BoundsError(P,k))
	end
end

"""
	vertices(P::AbstractPath)
Return an array of the vertices (endpoints of the curves) of the path `P`. The length is one greater than the number of curves in `P`.
"""
function vertices(P::AbstractPath)
	[ vertex(P,k) for k = 1:length(P)+1]
end

"""
	isfinite(P::AbstractPath)
Return `true` if the path is bounded in the complex plane (i.e., does not pass through infinity).
"""
isfinite(p::AbstractPath) = all(isfinite(s) for s in curves(p))

# iteration interface
eltype(::Type{AbstractPath}) = AbstractCurve
length(p::AbstractPath) = length(curves(p))
getindex(p::AbstractPath,k) = curve(p,k)
iterate(p::AbstractPath,state=1) = state > length(curves(p)) ? nothing : (p[state], state+1)

"""
	point(P::AbstractPath,t::Real)
	P(t)
Compute the point along path `P` at parameter value `t`. Values of `t` in [k,k+1] correspond to values in [0,1] along curve k of the path, for k = 1,2,...,length(P)-1.

	point(P::AbstractPath,t::AbstractVector)
Vectorize the `point` method for path `P`.
"""
point(p::AbstractPath,t::AbstractArray{T}) where T<:Real = [point(p,t) for t in t]

# This determines how to parse a parameter of a path. Overloaded later for the case of a closed path.
function sideargs(p::AbstractPath,t)
	n = length(p)
	if (t < 0) || (t > n)
		throw(BoundsError(p,t))
	end
	if t==n
		return n,1
	else
		return 1+floor(Int,t),t%1
	end
end

function point(p::AbstractPath,t::Real)
	k,s = sideargs(p,t)
	point(curve(p,k),s)
end

"""
	tangent(P::AbstractPath,t::Real)
Compute the complex-valued tangent along path `P` at parameter value `t`. Values of `t` in [k,k+1] correspond to values in [0,1] along curve k of the path, for k = 1,2,...,length(P)-1. The result is not well-defined at an integer value of `t`.
"""
function tangent(p::AbstractPath,t::Real)
	k,s = sideargs(p,t)
	tangent(curve(p,k),s)
end

"""
	unittangent(P::AbstractPath,t::Real)
Compute the complex-valued unit tangent along path `P` at parameter value `t`. Values of `t` in [k,k+1] correspond to values in [0,1] along curve k of the path, for k = 1,2,...,length(P)-1. The result is not well-defined at an integer value of `t`.
"""
function unittangent(p::AbstractPath,t::Real)
	k,s = sideargs(p,t)
	unittangent(curve(p,k),s)
end

"""
	normal(P::AbstractPath,t::Real)
Compute a complex-valued normal to path `P` at parameter value `t`. Values of `t` in [k,k+1] correspond to values in [0,1] along curve k of the path, for k = 1,2,...,length(P)-1. The result is not well-defined at an integer value of `t`.
"""
function normal(p::AbstractPath,t::Real)
	k,s = sideargs(p,t)
	normal(curve(p,k),s)
end

conj(p::AbstractPath) = typeof(p)(conj.(curves(p)))

reverse(p::AbstractPath) = typeof(p)(reverse(reverse.(curves(p))))

+(p::AbstractPath,z::Number) = typeof(p)([c+z for c in curves(p)])
+(z::Number,p::AbstractPath) = typeof(p)([z+c for c in curves(p)])

-(p::AbstractPath) = typeof(p)([-c for c in curves(p)])
-(p::AbstractPath,z::Number) = typeof(p)([c-z for c in curves(p)])
-(z::Number,p::AbstractPath) = typeof(p)([z-c for c in curves(p)])

*(p::AbstractPath,z::Number) = typeof(p)([c*z for c in curves(p)])
*(z::Number,p::AbstractPath) = typeof(p)([z*c for c in curves(p)])

/(p::AbstractPath,z::Number) = typeof(p)([c/z for c in curves(p)])
/(z::Number,p::AbstractPath) = z*inv(p)
inv(p::AbstractPath) = typeof(p)([inv(c) for c in curves(p)])

"""
	isapprox(P1::AbstractPath,R2::AbstractPath; tol=<default>)
	P1 ≈ P2       (type "\\approx" followed by tab key)
Determine whether `P1` and `P2` represent the same path, up to tolerance `tol`, irrespective of the parameterization of its curves.
"""
function isapprox(P1::AbstractPath,P2::AbstractPath;tol=DEFAULT[:tol])
	if length(P1) != length(P2)
		return false
	else
		c1,c2 = curves(P1),curves(P2)
		all( isapprox(c1[k],c2[k]) for k in eachindex(c1) )
	end
end
isapprox(::AbstractCurve,::AbstractPath;kw...) = false
isapprox(::AbstractPath,::AbstractCurve;kw...) = false

"""
	dist(z,P::AbstractPath)
Find the distance from the path `P` to the point `z`.
"""
dist(z::Number,P::AbstractPath) = minimum(dist(z,C) for C in P)

"""
	closest(z,P::AbstractPath)
Find the point on the path `P` that lies closest to `z`.
"""
function closest(z::Number,P::AbstractPath)
	k = argmin( [dist(z,s) for s in curves(P)] )
	closest(z,side(P,k))
end

intersect(P::AbstractPath,C::AbstractCurve) = ∪( [intersect(s,C) for s in curves(P)]...  )
intersect(P1::AbstractPath,P2::AbstractPath) = ∪( [intersect(P1,s) for s in curves(P2)]...  )

function show(io::IO,P::AbstractPath)
	print(IOContext(io,:compact=>true),typeof(P)," with ",length(P)," curves")
end
function show(io::IO,::MIME"text/plain",P::AbstractPath)
	print(io,typeof(P)," with ",length(P)," curves")
end

plotdata(P::AbstractPath) = vcat(plotdata.(P)...)

abstract type AbstractClosedPath <: AbstractPath end

"""
	curve(P::AbstractClosedPath,k::Integer)
Return the `k`th curve in the path `P`. The index is applied circularly; e.g, if the closed path has n curves, then ...,1-n,1,1+n,... all refer to the first curve.
"""
curve(p::AbstractClosedPath,k::Integer) = curves(p)[mod(k-1,length(p))+1]

"""
	vertex(P::AbstractPath,k::Integer)
Return the `k`th vertex of the path `P`. The index is applied circularly; e.g, if the closed path has n curves, then ...,1-n,1,1+n,... all refer to the first vertex.
"""
vertex(P::AbstractClosedPath,k::Integer) = point(curve(P,k),0)

"""
	vertices(P::AbstractClosedPath)
Return an array of the unique vertices (endpoints of the curves) of the closed path `P`. The length is equal the number of curves in `P`, i.e., the first/last vertex is not duplicated.
"""
vertices(P::AbstractClosedPath) = [ vertex(P,k) for k = 1:length(P) ]

function sideargs(p::AbstractClosedPath,t)
	n = length(p)
	return 1+mod(floor(Int,t),n), t%1
end

function winding(P::AbstractClosedPath,z::Number)
	# Integrate around the boundary
	w = 0
	for s in P
		f = t -> imag(tangent(s,t)/(point(s,t)-z))
		w += intadapt(f,0,1,1e-4)
	end
	return round(Int,w/(2π))
end

isinside(z::Number,P::AbstractClosedPath) = winding(P,z) != 0
isoutside(z::Number,P::AbstractClosedPath) = winding(P,z) == 0

#
# Concrete implementations
#

# Path
"""
	(type) Path
Generic implementation of an `AbstractPath`.
"""
struct Path <: AbstractPath
	curve
	function Path(c::AbstractVector{T};tol::Real=DEFAULT[:tol]) where T<:AbstractCurve
		n = length(c)
		for k = 1:n-1
			@assert isapprox(point(c[k],1.0),point(c[k+1],0.0),rtol=tol,atol=tol) "Curve endpoints do not match for pieces $(k) and $(k+1)"
		end
		new(c)
	end
end
"""
	Path(c::AbstractVector; tol=<default>)
Given a vector `c` of curves, construct a path. The path is checked for continuity (to tolerance `tol`) at the interior vertices.
"""
Path(c::AbstractCurve) = Path([c])

curves(p::Path) = p.curve
"""
	arclength(P::Path)
Compute the arclength of the path `P`.
"""
arclength(p::Path) = sum(arclength(c) for c in p)
(p::Path)(t) = point(p,t)

# ClosedPath
"""
	(type) ClosedPath
Generic implementation of an `AbstractClosedPath`.
"""
struct ClosedPath <: AbstractClosedPath
	curve
	function ClosedPath(p::AbstractVector{T};tol=DEFAULT[:tol]) where T<:AbstractCurve
		q = Path(p)
		zi,zf = point(q,0),point(q,length(q))
		@assert isapprox(zi,zf,rtol=tol,atol=tol) || (isinf(zi) && isinf(zf)) "Path endpoints do not match"
		new(p)
	end
end
"""
	ClosedPath(c::AbstractVector; tol=<default>)
	ClosedPath(P::Path; tol=<default>)
Given a vector `c` of curves, or an existing path, construct a closed path. The path is checked for continuity (to tolerance `tol`) at all of the vertices.
"""
ClosedPath(c::AbstractCurve) = ClosedPath([c])
ClosedPath(c::AbstractClosedPath) = c
ClosedPath(p::Path;kw...) = ClosedPath(p.curve;kw...)

curves(p::ClosedPath) = p.curve
arclength(p::ClosedPath) = sum(arclength(c) for c in p)
(p::ClosedPath)(t) = point(p,t)

# Find a circle that fully encloses all the finite vertices and some points of a path.
function enclosing_circle(p::AbstractPath,expansion=2)
	#z = [vertices(p);point(p,0:1/50:1)]
	z = vertices(p)
	return Circle(enclosing_circle(filter(isfinite,z),expansion)...)
end

include("polygons.jl")


text:
@testset "lambda pattern" begin
    @testset "example in docstring" begin
        xs = [(1, 2), (1, 3), (1, 4)]
        @test map((@λ (1, x) -> x), xs) == [2, 3, 4]

        @test (2, 3) |> @λ begin
            1 -> 2
            2 -> 7
            (a, b) -> (a + b) == 5
        end
    end

    @testset "gen_lambda case `:(\$a -> \$(b...))`" begin
        @test (2, 3) |> @λ (a, b) -> a == 2
    end
end

ext: jl


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Reformulate a R-DRO problem with Wasserstein-metric based ambiguity set into a standard form LP.
See Appendix C in the CLVR paper.
"""
function droreformuation_wmetric_hinge_standardformnormalized(
    yX_T::SparseMatrixCSC{Float64, Int},
    κ::Float64,
    ρ::Float64
)

    dim_dataset, num_dataset = size(yX_T)

    v1_n_T = transpose(sparse(ones(num_dataset)))
    v1_d_T = transpose(sparse(ones(dim_dataset)))
    v0_n_T = transpose(spzeros(num_dataset))

    I_nn = sparse(I, num_dataset, num_dataset)
    I_dd = sparse(I, dim_dataset, dim_dataset)
    O_nn = spzeros(num_dataset, num_dataset)
    O_dn = spzeros(dim_dataset, num_dataset)
    O_nd = spzeros(num_dataset, dim_dataset)
    O_dd = spzeros(dim_dataset, dim_dataset)

    A_T = [
        -I_nn       I_nn      I_nn      O_nd      O_nd  ;
        I_nn       O_nn      I_nn      O_nd      O_nd  ;
        O_nn      -I_nn     -I_nn      O_nd      O_nd  ;
        O_nn       O_nn     -I_nn      O_nd      O_nd  ;
        O_dn       yX_T      O_dn      I_dd      I_dd  ;
        O_dn      -yX_T      O_dn     -I_dd     -I_dd  ;
        O_dn       O_dn      O_dn     -I_dd      O_dd  ;
        O_dn       O_dn      O_dn      O_dd      I_dd  ;
        -2*κ*v1_n_T v0_n_T    v0_n_T    v1_d_T   -v1_d_T;
        2*κ*v1_n_T v0_n_T    v0_n_T   -v1_d_T    v1_d_T;
    ]
    b = sparsevec([
        spzeros(num_dataset);
        sparse(ones(num_dataset));
        2 * sparse(ones(num_dataset));
        spzeros(dim_dataset);
        spzeros(dim_dataset);
    ])

    # Normalize all row norms of A to 1
    divrownorm_A = map(x -> 1 / norm(x), eachcol(A_T))
    range_A = collect(1:A_T.n)
    diag_divrownorm_A = sparse(range_A, range_A, divrownorm_A)
    A_T *= diag_divrownorm_A
    b = diag_divrownorm_A * b

    # for i in 1:size(A_T)[2]
    #     A_T[:, i] = A_T[:, i] / rownorm_A[i]
    #     b[i] = b[i] / rownorm_A[i]
    # end

    c = [
        1 / num_dataset * ones(num_dataset);
        zeros(num_dataset);
        zeros(num_dataset);
        zeros(num_dataset);
        zeros(dim_dataset);
        zeros(dim_dataset);
        zeros(dim_dataset);
        zeros(dim_dataset);
        ρ;
        -ρ;
    ]
    A_T, b, c
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
"""
    circmean(x)

Compute the circular mean of `x` in radians.

[1] N. I. Fisher, Statistical Analysis of Circular Data. Cambridge University Press, 1993.
"""
function circmean(x)
    s = mean(sin, x)
    c = mean(cos, x)

    atan(s, c)
end

"""
    circmeand(x)

Compute the circular mean of `x` in degrees.

[1] N. I. Fisher, Statistical Analysis of Circular Data. Cambridge University Press, 1993.
"""
function circmeand(x)
    s = mean(sind, x)
    c = mean(cosd, x)

    atand(s, c)
end

"""
    circstd(x)

Compute the circular standard deviation of `x` in radians.

[1] N. I. Fisher, Statistical Analysis of Circular Data. Cambridge University Press, 1993.
"""
function circstd(x)
    s = mean(sin, x)
    c = mean(cos, x)

    sqrt(-2*log(hypot(c, s)))
end

"""
    circstdd(x)

Compute the circular standard deviation of `x` in degrees.

[1] N. I. Fisher, Statistical Analysis of Circular Data. Cambridge University Press, 1993.
"""
function circstdd(x)
    s = mean(sind, x)
    c = mean(cosd, x)

    rad2deg(sqrt(-2*log(hypot(c, s))))
end

"""
    mean_std_range(x, events) -> Tuple(avg, std)

Find the average range and range variability of `x` for the set of all intervals given by
`events`.
"""
function mean_std_range(x, events::AbstractVector{Int})
    mi, mx = intervalextrema(x, events)
    roms = mx - mi

    return mean(roms), std(roms)
end

"""
    avgextrema(x, events) -> Tuple(min, max)

Find the average minima and maxima for the set of all intervals given by `events`.
"""
function avgextrema(x, events::AbstractVector{Int})
    mi, mx = intervalextrema(x, events)

    return mean(mi), mean(mx)
end

function intervalextrema(x::AbstractVector{T}, events) where T
    rgs = intervals(events; endincluded=false)
    mi = Vector{T}(undef, length(rgs))
    mx = Vector{T}(undef, length(rgs))

    intervalextrema!(mi, mx, x, rgs)
end

function intervalextrema!(mi, mx, x, rgs)
    @inbounds for (i, rg) in enumerate(rgs)
        mi[i], mx[i] = extrema(view(x, rg))
    end

    return mi, mx
end


text:
"""

```
random_move!(x_ref::AbstractVector{T},
             irefatom::Int,
             sides::T,
             x_new::AbstractVector{T}, RNG) where T
```

Function that generates a new random position for a molecule.

The new position is returned in `x`, a previously allocated array.

`x_solvent_random` might be a view of the array that contains all the solvent molecules.

"""
function random_move!(
    x_ref::AbstractVector{T},
    irefatom::Int,
    sides::T,
    x_new::AbstractVector{T},
    RNG,
) where {T}

    # To avoid boundary problems, the center of coordinates are generated in a
    # much larger region, and wrapped aftwerwards
    scale = 100.0

    # Generate random coordiantes for the center of mass
    ncm1 = -scale * sides[1] / 2 + random(RNG, Float64) * scale * sides[1]
    ncm2 = -scale * sides[2] / 2 + random(RNG, Float64) * scale * sides[2]
    ncm3 = -scale * sides[3] / 2 + random(RNG, Float64) * scale * sides[3]
    newcm = T(ncm1, ncm2, ncm3)

    # Generate random rotation angles
    beta = (2 * π) * random(RNG, Float64)
    gamma = (2 * π) * random(RNG, Float64)
    theta = (2 * π) * random(RNG, Float64)

    # Copy the coordinates of the molecule chosen to the random-coordinates vector
    @. x_new = x_ref

    # Take care that this molecule is not split by periodic boundary conditions, by
    # wrapping its coordinates around its reference atom
    wrap!(x_new, sides, x_ref[irefatom])

    # Move molecule to new position
    move!(x_new, newcm, beta, gamma, theta)

    nothing
end


text:
"""
    RegEdit.RegistryCache(path, registries=Dict())

Represents a local cache of registry repositories rooted at `path`, where each registry is
stored in a subdirectory corresponding to the registry's UUID.

Maintains a dictionary `registries` which maps registry repository URLs to UUIDs.
"""
struct RegistryCache
    path::String
    registries::Dict{String, UUID}
end

RegistryCache(path) = RegistryCache(path, Dict())

const REGISTRY_CACHE = RegistryCache("registries")

path(cache::RegistryCache) = cache.path
path(cache::RegistryCache, reg_uuid::UUID) = joinpath(path(cache), string(reg_uuid))
function path(cache::RegistryCache, registry_url::AbstractString)
    path(cache, cache.registries[registry_url])
end

"""
    RegEdit.get_registry(registry_url)

Return a `GitRepo` object for an up-to-date copy of `registry`.
Update the existing copy if available.
"""
function get_registry(
    registry_url::String;
    gitconfig::Dict=Dict(),
    cache::RegistryCache=REGISTRY_CACHE,
)
    if haskey(cache.registries, registry_url)
        registry_path = path(cache, registry_url)

        if !ispath(registry_path)
            mkpath(path(cache))
            LibGit2.clone(registry_url, registry_path, branch="master")
        else
            # this is really annoying/impossible to do with LibGit2
            git = gitcmd(registry_path, gitconfig)
            run(`$git config remote.origin.url $registry_url`)
            run(`$git checkout -q -f master`)
            # uses config because git versions <2.17.0 did not have the -P option
            run(`$git -c fetch.pruneTags fetch -q origin master`)
            run(`$git reset -q --hard origin/master`)
        end
    else
        registry_temp = mktempdir(mkpath(path(cache)))
        try
            LibGit2.clone(registry_url, registry_temp)
            reg = parse_registry(joinpath(registry_temp, "Registry.toml"))
            registry_uuid = cache.registries[registry_url] = reg.uuid
            registry_path = path(cache, registry_uuid)
            rm(registry_path, recursive=true, force=true)
            mv(registry_temp, registry_path)
        finally
            rm(registry_temp, recursive=true, force=true)
        end
    end

    return LibGit2.GitRepo(registry_path)
end

@auto_hash_equals struct RegistryData
    name::String
    uuid::UUID
    repo::Union{String, Nothing}
    description::Union{String, Nothing}
    packages::Dict{String, Dict{String, Any}}
    extra::Dict{String, Any}
end

function RegistryData(
    name::AbstractString,
    uuid::Union{UUID, AbstractString};
    repo::Union{AbstractString, Nothing}=nothing,
    description::Union{AbstractString, Nothing}=nothing,
    packages::Dict=Dict(),
    extras::Dict=Dict(),
)
    RegistryData(name, UUID(uuid), repo, description, packages, extras)
end

function Base.copy(reg::RegistryData)
    RegistryData(
        reg.name,
        reg.uuid,
        reg.repo,
        reg.description,
        deepcopy(reg.packages),
        deepcopy(reg.extra),
    )
end

function parse_registry(io::IO)
    data = TOML.parse(io)

    name = pop!(data, "name")
    uuid = pop!(data, "uuid")
    repo = pop!(data, "repo", nothing)
    description = pop!(data, "description", nothing)
    packages = pop!(data, "packages", Dict())

    RegistryData(name, UUID(uuid), repo, description, packages, data)
end

parse_registry(str::AbstractString) = open(parse_registry, str)

function TOML.print(io::IO, reg::RegistryData)
    println(io, "name = ", repr(reg.name))
    println(io, "uuid = ", repr(string(reg.uuid)))

    if reg.repo !== nothing
        println(io, "repo = ", repr(reg.repo))
    end

    if reg.description !== nothing
        # print long-form string if there are multiple lines
        if '\n' in reg.description
            print(io, """\n
                description = \"\"\"
                $(reg.description)\"\"\"
                """
            )
        else
            println(io, "description = ", repr(reg.description))
        end
    end

    for (k, v) in pairs(reg.extra)
        TOML.print(io, Dict(k => v), sorted=true)
    end

    println(io, "\n[packages]")
    for (uuid, data) in sort!(collect(reg.packages), by=first)
        print(io, uuid, " = { ")
        print(io, "name = ", repr(data["name"]))
        print(io, ", path = ", repr(data["path"]))
        println(io, " }")
    end

    nothing
end

function package_relpath(reg::RegistryData, pkg::Pkg.Types.Project)
    joinpath("$(uppercase(pkg.name[1]))", pkg.name)
end

function Base.push!(reg::RegistryData, pkg::Pkg.Types.Project)
    reg.packages[string(pkg.uuid)] = Dict(
        "name" => pkg.name, "path" => package_relpath(reg, pkg)
    )
    reg
end


text:
module TensorNetworkStatesGates

    export OneQubitGateTensor, TwoQubitGateTensor, TNS_Gate, TNS_OneQubitGate, TNS_TwoQubitGate, GateTensor
    abstract type GateTensor end
    abstract type TNS_Gate end
    abstract type TNS_OneQubitGate <: TNS_Gate end
    abstract type TNS_TwoQubitGate <: TNS_Gate end

    mutable struct OneQubitGateTensor <: GateTensor
        item::Array{ComplexF64, 2}
        target::Int
    end

    mutable struct TwoQubitGateTensor <: GateTensor
        item::Array{ComplexF64, 2}
        control::Int
        target::Int
    end

    mutable struct X <: TNS_OneQubitGate
        _target::Int64
    end

    mutable struct Y <: TNS_OneQubitGate
        _target::Int64
    end

    mutable struct Z <: TNS_OneQubitGate
        _target::Int64
    end

    mutable struct H <: TNS_OneQubitGate
        _target::Int64
    end

    mutable struct T <: TNS_OneQubitGate
        _target::Int64
    end

    mutable struct S <: TNS_OneQubitGate
        _target::Int64
    end

    mutable struct RX <: TNS_OneQubitGate
        _target::Int64
        _theta::Float64
    end

    mutable struct RY <: TNS_OneQubitGate
        _target::Int64
        _theta::Float64
    end

    mutable struct RZ <: TNS_OneQubitGate
        _target::Int64
        _theta::Float64
    end

    mutable struct U3 <: TNS_OneQubitGate
        _target::Int64
        _theta::Float64
        _phi::Float64
        _lambd::Float64
    end

    mutable struct CZ <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    mutable struct CX <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    mutable struct CY <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    mutable struct SWAP <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    mutable struct CP <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    mutable struct CRX <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    mutable struct CRY <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    mutable struct CRZ <: TNS_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    # export all UG_Gate and TwoQubitGate
    for n in names(@__MODULE__; all=true)
        if Base.isidentifier(n) && n ∉ (Symbol(@__MODULE__), :eval, :include)
            if @eval typeof($n) <: DataType
                if @eval $n <: TNS_Gate
                    @eval export $n
                end
            end
        end
    end
end

max_stars_repo_name: ryuNagai/blueqat-julia

max_issues_repo_name: ryuNagai/blueqat-julia

max_forks_repo_name: ryuNagai/blueqat-julia


text:
#=
It was proposed by Christian Goldbach that every odd composite number can be written as the sum of a prime and twice a square.

9 = 7 + 2×1^2
15 = 7 + 2×2^2
21 = 3 + 2×3^2
25 = 7 + 2×3^2
27 = 19 + 2×2^2
33 = 31 + 2×1^2

It turns out that the conjecture was false.

What is the smallest odd composite that cannot be written as the sum of a prime and twice a square?
=#

using Primes

function nextoddcomposite(x)
  n = x
  while true
    n += 2
    if isprime(n) == false
      return n
    end
  end
end

function calc()
  composite = 7
  while true
    composite = nextoddcomposite(composite)
    primesless = reverse(primes(composite))
    canbecalculated = false
    for p in primesless, i in 1:54
      if p + 2*(i^2) == composite
        canbecalculated = true
      end
    end
    if canbecalculated == false
      return composite
    end
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using Random; Random.seed!(1234)

include("src/data.jl")
using .DataUtil
include("src/forward_stepwise/selective_inference.jl")
using .SIforSFS
include("src/forward_stepwise/data_splitting.jl")
import .DS
using LinearAlgebra:I
using Distributions
using ProgressMeter

const N = 100
const β₀ = 0.0
const 𝛃 = (Float64)[0.25, 0.25, 0, 0, 0]
const k = 3

# TPR experiment
function TPR_experiment(β₀, 𝛃, N, k; event=Active(), α=0.05, denominator_upper=100, iter_size=5)
    H_1 = filter(i -> !iszero(𝛃[i]), 1:length(𝛃)) # {j ∈ {1,...,d} | H_{1,j} is true}
    TPR = Array{Float64}(undef, length(H_1), iter_size)
    @showprogress "[$(N), $(event)]: " for l = 1:iter_size
        numerator = zeros(Int, length(H_1))
        denominator = zeros(Int, length(H_1))
        for _ = 1:denominator_upper
            𝐗, 𝐲_obs = centering(make_dataset(β₀, 𝛃, N)...)
            𝐴_obs, p_values = parametric_SFS_SI(𝐗, 𝐲_obs, k, I; selection_event=event)
            for (j, p_value) in zip(𝐴_obs, p_values)
                for (i, j′) in enumerate(H_1)
                    if j == j′ # H_{1,j} is true
                        numerator[i] += p_value < α ? 1 : 0
                        denominator[i] += 1
                        break
                    end
                end
            end
        end
        TPR[:, l] = numerator ./ denominator
    end
    return TPR
end

α = 0.05
nonzero_num = length(filter(x -> !iszero(x), 𝛃))
Ns = [50, 100, 150, 200]
iter_size = 5
TPR_A = fill(NaN, nonzero_num, iter_size, length(Ns))
TPR_Ao = fill(NaN, nonzero_num, iter_size, length(Ns))
TPR_As = fill(NaN, nonzero_num, iter_size, length(Ns))
TPR_Aso = fill(NaN, nonzero_num, iter_size, length(Ns))
TPR_DS = fill(NaN, nonzero_num, iter_size, length(Ns))
for (i, N) in enumerate(Ns)
    TPR_A[:,:,i] = Main.TPR_experiment(β₀, 𝛃, N, k; ℰ=Active(), α=α, iter_size=iter_size)
    TPR_Ao[:,:,i] = Main.TPR_experiment(β₀, 𝛃, N, k; ℰ=ActiveOrder(), α=α, iter_size=iter_size)
    TPR_As[:,:,i] = Main.TPR_experiment(β₀, 𝛃, N, k; ℰ=ActiveSign(), α=α, iter_size=iter_size)
    TPR_Aso[:,:,i] = Main.TPR_experiment(β₀, 𝛃, N, k; ℰ=ActiveSignOrder(), α=α, iter_size=iter_size)
    TPR_DS[:,:,i] = DS.TPR_experiment(β₀, 𝛃, N, k; α=α, iter_size=iter_size)
end

using Statistics
TPR_A_mean = reshape(mean(TPR_A, dims=2), nonzero_num, length(Ns))
TPR_A_SE = reshape(std(TPR_A, dims=2) ./ √iter_size, nonzero_num, length(Ns))
TPR_Ao_mean = reshape(mean(TPR_Ao, dims=2), nonzero_num, length(Ns))
TPR_Ao_SE = reshape(std(TPR_Ao, dims=2) ./ √iter_size, nonzero_num, length(Ns))
TPR_As_mean = reshape(mean(TPR_As, dims=2), nonzero_num, length(Ns))
TPR_As_SE = reshape(std(TPR_As, dims=2) ./ √iter_size, nonzero_num, length(Ns))
TPR_Aso_mean = reshape(mean(TPR_Aso, dims=2), nonzero_num, length(Ns))
TPR_Aso_SE = reshape(std(TPR_Aso, dims=2) ./ √iter_size, nonzero_num, length(Ns))
TPR_DS_mean = reshape(mean(TPR_DS, dims=2), nonzero_num, length(Ns))
TPR_DS_SE = reshape(std(TPR_DS, dims=2) ./ √iter_size, nonzero_num, length(Ns))

using Plots
default(:bglegend, plot_color(default(:bg), 0.5))
default(:fglegend, plot_color(ifelse(isdark(plot_color(default(:bg))), :white, :black), 0.6))
Cs = [:red, :orange, :green, :blue, :purple]

p = Array{typeof(plot())}(undef, nonzero_num)
for j in 1:nonzero_num
    plot()
    # Active
    plot!(1:length(Ns), (TPR_A_mean - TPR_A_SE)[j,:]; fillrange=(TPR_A_mean + TPR_A_SE)[j,:], fillalpha=0.2, α=0, label=nothing, color=Cs[1])
    plot!(1:length(Ns), TPR_A_mean[j,:], label="Homotopy", color=Cs[1], markershape=:circle)
    # ActiveOrder
    plot!(1:length(Ns), (TPR_Ao_mean - TPR_Ao_SE)[j,:]; fillrange=(TPR_Ao_mean + TPR_Ao_SE)[j,:], fillalpha=0.2, α=0, label=nothing, color=Cs[2])
    plot!(1:length(Ns), TPR_Ao_mean[j,:], label="Homotopy-H", color=Cs[2], markershape=:circle)
    # ActiveSign
    plot!(1:length(Ns), (TPR_As_mean - TPR_As_SE)[j,:]; fillrange=(TPR_As_mean + TPR_As_SE)[j,:], fillalpha=0.2, α=0, label=nothing, color=Cs[3])
    plot!(1:length(Ns), TPR_As_mean[j,:], label="Homotopy-S", color=Cs[3], markershape=:circle)
    # ActiveSignOrder
    plot!(1:length(Ns), (TPR_Aso_mean - TPR_Aso_SE)[j,:]; fillrange=(TPR_Aso_mean + TPR_Aso_SE)[j,:], fillalpha=0.2, α=0, label=nothing, color=Cs[4])
    plot!(1:length(Ns), TPR_Aso_mean[j,:], label="Polytope", color=Cs[4], markershape=:circle)
    # DS
    plot!(1:length(Ns), (TPR_DS_mean - TPR_DS_SE)[j,:]; fillrange=(TPR_DS_mean + TPR_DS_SE)[j,:], fillalpha=0.2, α=0, label=nothing, color=Cs[5])
    plot!(1:length(Ns), TPR_DS_mean[j,:], label="DS", color=Cs[5], markershape=:circle)
    # plot setting
    plot!(; ylims=(0, 1))
    plot!(; legend=:topleft)
    plot!(; xlabel="Sample size", ylabel="True Positive Rate (TPR)")
    plot!(; guidefontsize=16, legendfontsize=12, tickfontsize=10)
    p[j] = plot!(; xticks=(1:length(Ns), Ns))
    savefig(p[j], "img/TPR/TPR$(j).pdf")
end


text:
# Stream management

export
    CuStream, CuDefaultStream, synchronize

const CuStream_t = Ptr{Cvoid}

mutable struct CuStream
    handle::CuStream_t
    ctx::CuContext
end

Base.unsafe_convert(::Type{CuStream_t}, s::CuStream) = s.handle

Base.:(==)(a::CuStream, b::CuStream) = a.handle == b.handle
Base.hash(s::CuStream, h::UInt) = hash(s.handle, h)

@enum(CUstream_flags, STREAM_DEFAULT      = 0x00,
                      STREAM_NON_BLOCKING = 0x01)

# FIXME: EnumSet from JuliaLang/julia#19470
Base.:|(x::CUstream_flags, y::CUstream_flags) =
    reinterpret(CUstream_flags, Base.cconvert(Unsigned, x) | Base.cconvert(Unsigned, y))

"""
    CuStream(flags=STREAM_DEFAULT)

Create a CUDA stream.
"""
function CuStream(flags::CUstream_flags=STREAM_DEFAULT)
    handle_ref = Ref{CuStream_t}()
    @apicall(:cuStreamCreate, (Ptr{CuStream_t}, Cuint),
                              handle_ref, flags)

    ctx = CuCurrentContext()
    obj = CuStream(handle_ref[], ctx)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(s::CuStream)
    if isvalid(s.ctx)
        @apicall(:cuStreamDestroy, (CuStream_t,), s)
    end
end

"""
    CuDefaultStream()

Return the default stream.
"""
@inline CuDefaultStream() = CuStream(convert(CuStream_t, C_NULL), CuContext(C_NULL))

"""
    synchronize(s::CuStream)

Wait until a stream's tasks are completed.
"""
synchronize(s::CuStream) = @apicall(:cuStreamSynchronize, (CuStream_t,), s)


text:
module SparseVolterraExamples

    using ApproxFun, NLsolve, BandedMatrices, BlockBandedMatrices, LinearAlgebra, BlockArrays, MultivariateOrthogonalPolynomials, Test, Plots

    include("triIntegroDiff.jl")
    include("triQEgen.jl")
    include("triVolterraFullKernel.jl")
    include("triVolterraMonKernel.jl")

    export triNonLinearVolterraObjective,
        triNonLinearIntegroDiffVolterraObjective,
        pochhammer,
        alternI,
        EvalJPat1mat,
        DirectEvalLHSP10at0,
        DirectEvalLHSP10atPrime0,
        DirectEvalLHSP10atPPrime0,
        DirectEvalLHSP10atPPPrime0,
        DirectEvalLHSP10atPPPPrime0,
        triQgen,
        triEygenP01,
        triQEygenP01,
        WLoweringP01P00,
        reflectPabtoPba,
        WLoweringP01Pm11,
        ReflectorP01,
        triEygenP00,
        triEygenP11,
        triExgenP01,
        triVolterraEQ1FullKernelSolver,
        triVolterraEQ2FullKernelSolver,
        triVolterraFullKernelFunP01,
        triVolterraFullKernelOpP01,
        triVolterraFullKernelOpP01,
        triBandedSgen,
        triVolterraEQ1Solver,
        triVolterraEQ2Solver,
        triVolterraFunP01,
        triVolterraOpP01,
        triJQEJbasislexP01,
        triJQEJbasisgridgenP01,
        triNoKernelDirect

end # module


text:
macro catchall(expr)
    quote
        try
            $(expr |> esc)

        catch e
            println(Base.stderr, e)

            for (exc, bt) in Base.catch_stack()
                showerror(Base.stderr, exc, bt)
                println()
            end
        end
    end
end


text:
module Configs

using ..Utilities

export Config

abstract type ConfigValue <: ComponentValue end
abstract type ConfigType <: ComponentType end

struct Edits <: ConfigType
    annotationPosition::Bool
    annotationTail::Bool
    annotationText::Bool
    axisTitleText::Bool
    colorbarPosition::Bool
    legendPosition::Bool
    legendText::Bool
    shapePosition::Bool
    titleText::Bool

    function Edits(;
        annotationPosition::Bool = false,
        annotationTail::Bool = false,
        annotationText::Bool = false,
        axisTitleText::Bool = false,
        colorbarPosition::Bool = false,
        legendPosition::Bool = false,
        legendText::Bool = false,
        shapePosition::Bool = false,
        titleText::Bool = false)
        return new(
            annotationPosition,
            annotationTail,
            annotationText,
            axisTitleText,
            colorbarPosition,
            legendPosition,
            legendText,
            shapePosition,
            titleText)
    end
end

const FrameMargins = RangeValue{Float64,(0.,0.5),(true,true)}

const ScrollZoom = FlagList{(:cartesian,:gl3d,:geo,:mapbox),(true,false)}

const DoubleClick = Enumerated{:DoubleClick,
                        (false,:reset,:autosize,Symbol("reset+autosize"))}

const DisplayModeBar = Enumerated{:DisplayModeBar,(:hover,true,false)}

struct Config <: Component
    staticPlot::Bool
    plotlyServerURL::String
    editable::Bool
    edits::Edits
    autosizable::Bool
    responsive::Bool
    fillFrame::Bool
    frameMargins::FrameMargins
    scrollZoom::ScrollZoom
    doubleClick::DoubleClick
    showAxisDraghandles::Bool
    showAxisRangeEntryBoxes::Bool
    showTips::Bool
    showLink::Bool
    linkText::String
    sendData::Bool
    showSources::Any
    displayModeBar::DisplayModeBar
    showSendToCloud::Bool
    modeBarButtonsToRemove::Any
    modeBarButtonsToAdd::Any
    modeBarButtons::Any
    toImageButtonOptions::Any
    displaylogo::Bool
    watermark::Bool
    plotGlPixelRatio::RangeValue{Float64,(1,4),(true,true)}
    setBackground::Any
    topojsonURL::String
    mapboxAccessToken::String
    logging::Bool
    queueLength::Size{Int}
    globalTransforms::Any
    locale::String
    locales::Any

    function Config(;
        staticPlot::Bool = false,
        plotlyServerURL::String = "https://plot.ly",
        editable::Bool = false,
        edits::Edits = Edits(),
        autosizable::Bool = false,
        responsive::Bool = false,
        fillFrame::Bool = false,
        frameMargins = 0.,
        scrollZoom::Union{String,Bool} = "gl3d+geo+mapbox",
        doubleClick::String = "reset+autosize",
        showAxisDraghandles::Bool = true,
        showAxisRangeEntryBoxes::Bool = true,
        showTips::Bool = true,
        showLink::Bool = false,
        linkText::String = "Edit chart",
        sendData::Bool = true,
        showSources::Any = false,
        displayModeBar::Union{String,Bool} = "hover",
        showSendToCloud::Bool = false,
        modeBarButtonsToRemove::Any = [],
        modeBarButtonsToAdd::Any = [],
        modeBarButtons::Any = false,
        toImageButtonOptions::Any = Dict(),
        displaylogo::Bool = true,
        watermark::Bool = false,
        plotGlPixelRatio = 2,
        setBackground::Any = "transparent",
        topojsonURL::String = "https://cdn.plot.ly/",
        mapboxAccessToken::String = "",
        logging::Bool = true,
        queueLength::Int = 0,
        globalTransforms::Any = [],
        locale::String = "en-US",
        locales::Any = Dict())
        return new(
            staticPlot,
            plotlyServerURL,
            editable,
            edits,
            autosizable,
            responsive,
            fillFrame,
            FrameMargins(frameMargins),
            ScrollZoom(scrollZoom),
            DoubleClick(doubleClick),
            showAxisDraghandles,
            showAxisRangeEntryBoxes,
            showTips,
            showLink,
            linkText,
            sendData,
            showSources,
            DisplayModeBar(displayModeBar),
            showSendToCloud,
            modeBarButtonsToRemove,
            modeBarButtonsToAdd,
            modeBarButtons,
            toImageButtonOptions,
            displaylogo,
            watermark,
            RangeValue((1,4),Float64(plotGlPixelRatio)),
            setBackground,
            topojsonURL,
            mapboxAccessToken,
            logging,
            Size{Int}(queueLength),
            globalTransforms,
            locale,
            locales)
    end
end
default["config"] = Config()

end


text:
const tracetype = Dict{Symbol,Int}(:live=>1, :dead=>2, :aux=>3)

const stockdomain = Dict{Symbol,String}(
:ALACRITY       =>"alacrity",
:AMPLITUDE      =>"amplitude",
:COHERENCE      =>"coherence",
:DELTA          =>"delta",
:DENSITY        =>"density",
:DEPTH          =>"depth",
:DIP            =>"dip",
:ENVELOPE       =>"envelope",
:EPSILON        =>"epsilon",
:ETA            =>"eta",
:FLEX_BINNED    =>"flex_binned",
:FOLD           =>"fold",
:FREQUENCY      =>"frequency",
:IMPEDANCE      =>"impedence",
:INCIDENCE_ANGLE=>"incidence_angle",
:MODEL_TRANSFORM=>"model_transform",
:ROTATION_ANGLE =>"rotation_angle",
:SEMBLANCE      =>"semblence",
:SLOTH          =>"sloth",
:SLOWNESS       =>"slowness",
:SPACE          =>"space",
:TIME           =>"time",
:UNKNOWN        =>"unknown",
:VELOCITY       =>"velocity",
:VS             =>"vs",
:VSVP           =>"vsvp",
:WAVENUMBER     =>"wavenumber"
)

const stockunit = Dict{Symbol,String}(
:DEGREES     =>"degrees",
:FEET        =>"feet",
:FT          =>"ft",
:HERTZ       =>"hertz",
:HZ          =>"hz",
:M           =>"m",
:METERS      =>"meters",
:MICROSEC    =>"microseconds",
:MILLISECONDS=>"milliseconds",
:MS          =>"ms",
:MSEC        =>"msec",
:SECONDS     =>"seconds",
:S           =>"seconds",
:NULL        =>"null",
:UNKNOWN     =>"unknown"
)

const stockdatatype = Dict{Symbol,String}(
:CMP       =>"CMP",
:CUSTOM    =>"CUSTOM",
:OFFSET_BIN=>"OFFSET_BIN",
:RECEIVER  =>"RECEIVER",
:SOURCE    =>"SOURCE",
:STACK     =>"STACK",
:UNKNOWN   =>"UNKNOWN"
)

const stockprop = Dict{Symbol,TracePropertyDef}(
:AMP_NORM =>TracePropertyDef("AMP_NORM", "Amplitude normalization factor",    Float32, 1),
:AOFFSET  =>TracePropertyDef("AOFFSET",  "Absolute value of offset",          Float32, 1),
:CDP      =>TracePropertyDef("CDP",      "CDP bin number",                    Int32,   1),
:CDP_ELEV =>TracePropertyDef("CDP_ELEV", "Elevation of CDP",                  Float32, 1),
:CDP_NFLD =>TracePropertyDef("CDP_NFLD", "Number of traces in CDP bin",       Int32,   1),
:CDP_SLOC =>TracePropertyDef("CDP_SLOC", "External CDP number",               Int32,   1),
:CDP_X    =>TracePropertyDef("CDP_X",    "X coordinate of CDP (float)",       Float32, 1),
:CDP_XD   =>TracePropertyDef("CDP_XD",   "X coordinate of CDP (double)",      Float64, 1),
:CDP_Y    =>TracePropertyDef("CDP_Y",    "Y coordinate of CDP (float)",       Float32, 1),
:CDP_YD   =>TracePropertyDef("CDP_YD",   "Y coordinate of CDP (double)",      Float64, 1),
:CHAN     =>TracePropertyDef("CHAN",     "Recording channel number",          Int32,   1),
:CMP_X    =>TracePropertyDef("CMP_X",    "Average of shot and receiver x",    Float32, 1),
:CMP_Y    =>TracePropertyDef("CMP_Y",    "Average of shot and receiver y",    Float32, 1),
:CR_STAT  =>TracePropertyDef("CR_STAT",  "Corr. autostatics receiver static", Float32, 1),
:CS_STAT  =>TracePropertyDef("CS_STAT",  "Corr. autostatics source static",   Float32, 1),
:DEPTH    =>TracePropertyDef("DEPTH",    "Source depth",                      Float32, 1),
:DISKITER =>TracePropertyDef("DISKITER", "Disk data iteration*",              Int32,   1),
:DMOOFF   =>TracePropertyDef("DMOOFF",   "Offset bin for DMO",                Int32,   1),
:DS_SEQNO =>TracePropertyDef("DS_SEQNO", "Input dataset sequence number*",    Int32,   1),
:END_ENS  =>TracePropertyDef("END_ENS",  "End-of-ensemble flag*",             Int32,   1),
:END_VOL  =>TracePropertyDef("END_VOL",  "End-of-volume flag*",               Int32,   1),
:EOJ      =>TracePropertyDef("EOJ",      "End of job flag*",                  Int32,   1),
:FB_PICK  =>TracePropertyDef("FB_PICK",  "First break pick",                  Float32, 1),
:FFID     =>TracePropertyDef("FFID",     "Field file ID number",              Int32,   1),
:FILE_NO  =>TracePropertyDef("FILE_NO",  "Sequential file number",            Int32,   1),
:FK_WAVEL =>TracePropertyDef("FK_WAVEL", "Wavelength of F-K domain trace",    Float32, 1),
:FK_WAVEN =>TracePropertyDef("FK_WAVEN", "Wavenumber of F-K domain trace",    Float32, 1),
:FNL_STAT =>TracePropertyDef("FNL_STAT", "Static to move to final datum",     Float32, 1),
:FRAME    =>TracePropertyDef("FRAME",    "Frame index in framework",          Int32,   1),
:FT_FREQ  =>TracePropertyDef("FT_FREQ",  "Frequency of F-T domain trace",     Float32, 1),
:GEO_COMP =>TracePropertyDef("GEO_COMP", "Geophone component (x,y,z)",        Int32,   1),
:HYPRCUBE =>TracePropertyDef("HYPRCUBE", "Hypercube index in framework",      Int32,   1),
:IF_FLAG  =>TracePropertyDef("IF_FLAG",  "ProMax IF_FLAG",                    Int32,   1),
:ILINE_NO =>TracePropertyDef("ILINE_NO", "3D iline number",                   Int32,   1),
:LEN_SURG =>TracePropertyDef("LEN_SURG", "Length of surgical mute taper",     Float32, 1),
:LINE_NO  =>TracePropertyDef("LINE_NO",  "Line number (hased line name)*",    Int32,   1),
:LSEG_END =>TracePropertyDef("LSEG_END", "Line segment end*",                 Int32,   1),
:LSEG_SEQ =>TracePropertyDef("LSEG_SEQ", "Line segment sequence number*",     Int32,   1),
:NA_STAT  =>TracePropertyDef("NA_STAT",  "Portion of static not applied",     Float32, 1),
:NCHANS   =>TracePropertyDef("NCHANS",   "Number of channels of source",      Int32,   1),
:NDATUM   =>TracePropertyDef("NDATUM",   "Floating NMO Datum",                Float32, 1),
:NMO_STAT =>TracePropertyDef("NMO_STAT", "NMO datum static",                  Float32, 1),
:NMO_APLD =>TracePropertyDef("NMO_APLD", "NMO applied to traces",             Int32,   1),
:OFB_CNTR =>TracePropertyDef("OFB_CNTR", "Offset bin center",                 Float32, 1),
:OFB_NO   =>TracePropertyDef("OFB_NO",   "Offset bin number",                 Int32,   1),
:OFFSET   =>TracePropertyDef("OFFSET",   "Signed source-receiver offset",     Float32, 1),
:PAD_TRC  =>TracePropertyDef("PAD_TRC",  "Artifically padded trace",          Int32,   1),
:PR_STAT  =>TracePropertyDef("PR_STAT",  "Power autostatics receiver static", Float32, 1),
:PS_STAT  =>TracePropertyDef("PS_STAT",  "Power autostatics source static",   Float32, 1),
:R_LINE   =>TracePropertyDef("R_LINE",   "Receiver line number",              Int32,   1),
:REC_ELEV =>TracePropertyDef("REC_ELEV", "Receiver elevation",                Float32, 1),
:REC_H2OD =>TracePropertyDef("REC_H2OD", "Water depth at receiver",           Float32, 1),
:REC_NFLD =>TracePropertyDef("REC_NFLD", "Receiver fold",                     Int32,   1),
:REC_SLOC =>TracePropertyDef("REC_SLOC", "Receiver index number (internal)*", Int32,   1),
:REC_STAT =>TracePropertyDef("REC_STAT", "Total static for receiver",         Float32, 1),
:REC_X    =>TracePropertyDef("REC_X",    "Receiver X coordinate (float)",     Float32, 1),
:REC_XD   =>TracePropertyDef("REC_XD",   "Receiver X coordinate (double)",    Float64, 1),
:REC_Y    =>TracePropertyDef("REC_Y",    "Receiver Y coordinate (float)",     Float32, 1),
:REC_YD   =>TracePropertyDef("REC_YD",   "Receiver Y coordinate (double)",    Float64, 1),
:REPEAT   =>TracePropertyDef("REPEAT",   "Repeated data copy number",         Int32,   1),
:S_LINE   =>TracePropertyDef("S_LINE",   "Swath or sail line number",         Int32,   1),
:SAMPLE   =>TracePropertyDef("SAMPLE",   "Sample index in framework",         Int32,   1),
:SEQNO    =>TracePropertyDef("SEQNO",    "Sequence number in ensemble",       Int32,   1),
:SEQ_DISK =>TracePropertyDef("SEQ_DISK", "Trace sequence number from disk",   Int32,   1),
:SG_CDP   =>TracePropertyDef("SG_CDP",   "Super gather CDP number",           Int32,   1),
:SIN      =>TracePropertyDef("SIN",      "Source index number*",              Int32,   1),
:SLC_TIME =>TracePropertyDef("SLC_TIME", "Time slice input",                  Float32, 1),
:SMH_CDP  =>TracePropertyDef("SMH_CDP",  "Number of CDP's in supergather",    Int32,   1),
:SOU_COMP =>TracePropertyDef("SOU_COMP", "Source component (x,y,z)",          Int32,   1),
:SOU_ELEV =>TracePropertyDef("SOU_ELEV", "Source elevation",                  Float32, 1),
:SOU_H2OD =>TracePropertyDef("SOU_H2OD", "Water depth at source",             Float32, 1),
:SOU_SLOC =>TracePropertyDef("SOU_SLOC", "External source location number",   Int32,   1),
:SOU_STAT =>TracePropertyDef("SOU_STAT", "Total static for source",           Float32, 1),
:SOU_X    =>TracePropertyDef("SOU_X",    "Source X coordinate (float)",       Float32, 1),
:SOU_XD   =>TracePropertyDef("SOU_XD",   "Source X coordinate (double)",      Float64, 1),
:SOU_Y    =>TracePropertyDef("SOU_Y",    "Source Y coordinate (float)",       Float32, 1),
:SOU_YD   =>TracePropertyDef("SOU_YD",   "Source Y coordinate (double)",      Float64, 1),
:SOURCE   =>TracePropertyDef("SOURCE",   "Live source number (user-defined)", Int32,   1),
:SKEWSTAT =>TracePropertyDef("SKEWSTAT", "Multiplex skew static",             Float32, 1),
:SR_AZIM  =>TracePropertyDef("SR_AZIM",  "Source to receiver azimuth",        Float32, 1),
:SRF_SLOC =>TracePropertyDef("SRF_SLOC", "External receiver location number", Int32,   1),
:TFULL_E  =>TracePropertyDef("TFULL_E",  "End time of full samples",          Float32, 1),
:TFULL_S  =>TracePropertyDef("TFULL_S",  "Start time of full samples",        Float32, 1),
:TIME_IND =>TracePropertyDef("TIME_IND", "Time sample index",                 Int32,   1),
:TLIVE_E  =>TracePropertyDef("TLIVE_E",  "End time of live samples",          Float32, 1),
:TLIVE_S  =>TracePropertyDef("TLIVE_S",  "Start time of live samples",        Float32, 1),
:TOT_STAT =>TracePropertyDef("TOT_STAT", "Total static for this trace",       Float32, 1),
:TRACE    =>TracePropertyDef("TRACE",    "Trace index in framework",          Int32,   1),
:TRC_TYPE =>TracePropertyDef("TRC_TYPE", "Trace type (data, aux, etc.)",      Int32,   1),
:TR_FOLD  =>TracePropertyDef("TR_FOLD",  "Actual trace fold",                 Float32, 1),
:TRACENO  =>TracePropertyDef("TRACENO",  "Trace number in seismic line*",     Int32,   1),
:TRIMSTAT =>TracePropertyDef("TRIMSTAT", "Trim static",                       Float32, 1),
:UPHOLE   =>TracePropertyDef("UPHOLE",   "Source uphole time",                Float32, 1),
:VOLUME   =>TracePropertyDef("VOLUME",   "Volume index in framework",         Int32,   1),
:WB_TIME  =>TracePropertyDef("WB_TIME",  "Water bottom time",                 Float32, 1),
:XLINE_NO =>TracePropertyDef("XLINE_NO", "3D crossline number",               Int32,   1)
)


text:
#####
### gapanalysis_stay_change.jl
# this file produces graphs and the final cubic data needed for the
# player session cutoff threshold.
####
#####
### test similarity to see if there are major differences between days. This
### should not vary much.
#####

# using DataFrames
# df = readtable("C:/Users/Peter.Mancini/Documents/Datafiles/playersessiongap_post.csv")
### uncomment above if you are running this after having already run PSG_01.jl and are
### simply developing new code. Saves several minutes.

for i in 1:7
  println("on day ", i, " similarity is ", quantile(df[df[:day] .== i, :similarity], [0.05, 0.10, 0.15, 0.20, 0.80, 0.90, 0.95]))
end

#####
### bin the similarity data into 1% sized bins
### create other metrics that will be used. Assumption, 25%/75% is a good
### starting point for dividing the data in to players we are certain
### are the same person, and players we are certain are different.
#####
include( "./DataFrameBinning.jl" )
df[:idx] = find_bins(df[:similarity], range(0.0,0.01,101))
change_mid = sum(counts(sample(df[df[:similarity] .< 0.25,:idx], 10000), 0:20))/2
change_cumsum = cumsum(counts(sample(df[df[:similarity] .< 0.25,:idx], 10000), 0:20))
stay_mid = sum(counts(sample(df[df[:similarity] .> 0.75,:idx], 10000), 92:102))/2
stay_cumsum = cumsum(counts(sample(df[df[:similarity] .> 0.75,:idx], 10000), 92:102))

sim_stay_suggestion = Int(100)
sim_change_suggestion = Int(0)

for i in length(stay_cumsum):-1:1
  if stay_cumsum[i] >= stay_mid
    sim_stay_suggestion = 100 - i
  end
end

for i in 1:length(change_cumsum)
  if change_cumsum[i] >= change_mid
    sim_change_suggestion = i
  end
end

##### similarity thresholds
sim_stay = 85 ## or you could set it to sim_stay_suggestion
sim_stay = sim_stay_suggestion
sim_change = 25 ## or you could set it to sim_change_suggestion
sim_change = sim_change_suggestion

###
# Now let's look at the values to get a clearer understanding
###
StatsBase.summarystats(df[df[:similarity] .<= sim_change/100, :gap])
StatsBase.summarystats(df[df[:similarity] .>= sim_stay/100, :gap])
change_gap = by(df[df[:similarity] .< sim_change/100, :], :hour, x->median(x[:gap]))
Stay_gap = by(df[df[:similarity] .> sim_stay/100, :], :hour, x->median(x[:gap]))

###
# Start plotting the results
###
using Gadfly # this is basically ggplot2, very similar syntax
# hours go from 0 to 23, not 1 to 24.

l1 = layer(change_gap, x=0:23, y=:x1, Geom.point, Theme(default_color=colorant"orange"))

l2 = layer(Stay_gap,x=0:23, y=:x1, Geom.point, Theme(default_color=colorant"purple"))

# Plot 1, hour by hour gap analysis. Pretty basic.

P1 = plot(l1, l2, Guide.xlabel("Hour"),
  Guide.ylabel("seconds"),
  Guide.title("Gap Analysis"),
  Guide.manual_color_key("Legend", ["New Player", "Same Player"], ["orange", "purple"]))
draw(PNG("GapAnalysis_byhour.png", 6inch, 4inch), P1)
draw(PNG("GapAnalysis_byhour.png", 6inch, 4inch), P1)

stay_seconds = by(df[df[:similarity] .> sim_stay/100, :], :hour, x->quantile(x[:gap],[0.5]))
change_seconds = by(df[df[:similarity] .< sim_change/100, :], :hour, x->quantile(x[:gap],[0.5]))

stay_seconds_bydayhour = by(df[df[:similarity] .> sim_stay/100, :], [:day, :hour], x->quantile(x[:gap],[0.5]))
change_seconds_bydayhour = by(df[df[:similarity] .< sim_change/100, :], [:day, :hour], x->quantile(x[:gap],[0.5]))
# Plot 2, hour by hour and day by day. More sophisticated. Proves we need a
# 24x7 grid vs just an hour by hour array of thresholds.

P2 = plot(stay_seconds_bydayhour, ygroup=:day,x=:hour, y=:x1,  Geom.subplot_grid(Geom.bar), Guide.ylabel("seconds by day"))
draw(PNG("GapAnalysis_bydayhour.png", 6inch, 8inch), P2)
draw(PNG("GapAnalysis_bydayhour.png", 6inch, 8inch), P2)

#####
### The final analysis, we need similarity and gaps, hour by hour and day by day
### We also need data for our error bars. That data will be really useful in
### building our 24x7 grid of values.
### I chose 90%+ similarity for the players who stayed on the machine,
### and 25%- similarity for the change over in players. This produced a useful
### distinction between the two groups.

stay_gapmean = by(df[df[:similarity] .> sim_stay/100,:], [:day, :hour], x -> [quantile(x[:gap],[0.5]) std(x[:gap])])
Ns = by(df[df[:similarity] .> sim_stay/100,:], [:day, :hour], x-> size(x)[1])
stay_ymaxes = Array{Float64}(7,24)
stay_ymaxes2 = Array{Float64}(7,24)
stay_ymins = Array{Float64}(7,24)
for d in 1:7, h in 1:24
  i = (stay_gapmean[:day] .== d) .&  (stay_gapmean[:hour] .== (h-1))
  stay_ymaxes[d,h] = stay_gapmean[i,:x1][1] + (1.96 * stay_gapmean[i,:x2][1]) / sqrt(Ns[i,:x1][1])
  stay_ymaxes2[d,h] = stay_gapmean[i,:x1][1] + 1.5*(1.96 * stay_gapmean[i,:x2][1]) / sqrt(Ns[i,:x1][1])
  if stay_ymaxes2[d,h] > 300 stay_ymaxes2[d,h] = 300 end
  stay_ymins[d,h] = stay_gapmean[i,:x1][1] - (1.96 * stay_gapmean[i,:x2][1]) / sqrt(Ns[i,:x1][1])
end

change_gapmean = by(df[df[:similarity] .< sim_change/100,:], [:day, :hour], x -> [quantile(x[:gap],[0.5]) std(x[:gap])])
Ns = by(df[df[:similarity] .< sim_change/100,:], [:day, :hour], x-> size(x)[1])
change_ymaxes = Array{Float64}(7,24)
change_ymins = Array{Float64}(7,24)
for d in 1:7, h in 1:24
  i = (change_gapmean[:day] .== d) .&  (change_gapmean[:hour] .== (h-1))
  change_ymaxes[d,h] = change_gapmean[i,:x1][1] + (1.96 * change_gapmean[i,:x2][1]) / sqrt(Ns[i,:x1][1])
  change_ymins[d,h] = change_gapmean[i,:x1][1] - (1.96 * change_gapmean[i,:x2][1]) / sqrt(Ns[i,:x1][1])
end

### Plotting

stay_layer = layer(stay_gapmean[stay_gapmean[:day] .== 6, :],
  x=0:23, y=:x1, ymin=stay_ymins[6,:], ymax=stay_ymaxes[6,:],
  Geom.point, Geom.errorbar,
  Theme(default_color=colorant"orange"))
change_layer = layer(change_gapmean[change_gapmean[:day] .== 6, :],
  x=0:23, y=:x1, ymin=change_ymins[6,:], ymax=change_ymaxes[6,:],
  Geom.point, Geom.errorbar,
  Theme(default_color=colorant"purple"))
smooth_layer = layer(x=0:23, y=stay_ymaxes2[6,:], Geom.smooth(method=:loess, smoothing=0.85))
P3 = plot(stay_layer, change_layer, smooth_layer, Theme(background_color=colorant"old lace"),
  Guide.title("Player Session Gap Analysis"),
  Guide.XLabel("Hour of Day"),
  Guide.YLabel("Typical Gap Time (Seconds)"),
  Guide.manual_color_key("Legend", ["Same Player", "New Player"], ["orange", "purple"]))

draw(PNG("GapAnalysis_bydayhour_WinStar_final.png", 6inch, 8inch), P3)
draw(PNG("GapAnalysis_bydayhour_Winstar_final.png", 6inch, 8inch), P3)

#####
### Achievement Unlocked! <bing!>
#####
function estimate_session_gap()
  # if the stay top error index is below the change error index
  # then process normally
  return_array = Array{Float64}(7,24)
  for day in 1:7, hour in 1:24
    change_mean = (change_gapmean[(change_gapmean[:hour] .== hour-1) .& (change_gapmean[:day] .== day), :x1])[1]
    change_min = change_ymins[day,hour]
    stay_mean = (stay_gapmean[(stay_gapmean[:hour] .== hour-1) .& (stay_gapmean[:day] .== day), :x1])[1]
    stay_max = stay_ymaxes[day,hour]
    if stay_max <= change_min # no collistion in error
      return_array[day,hour] = 300.0 < (change_min - stay_max)/2 + stay_max ? 300.0 : (change_min - stay_max)/2 + stay_max
    else
      # colision in error, find point between means
      return_array[day,hour] = 300.0 < (change_mean - stay_mean)/2 + stay_mean ? 300.0 : (change_mean - stay_mean)/2 + stay_mean
    end
  end
  return return_array
end

Base.writedlm("Array_Winstar.csv", estimate_session_gap(), ",")


text:
using Test
using Cracker
using ChainRules: rrule, unthunk

@testset "basic propagate ($T)" for T in [Float64, ComplexF64]
    A, B, C = track(rand(T, 2, 2)), track(rand(T, 2, 2)), track(rand(T, 4))
    Z = A + B * reshape(C, 2, 2)

    ret = abs2(sum(Z))
    Cracker.backpropagate!(ret, 1.0)

    uA,uB,uC = untrack.((A, B, C))
    T1, pb1 = rrule(reshape, uC, 2, 2)
    T2, pb2 = rrule(*, uB, T1)
    T3, pb3 = rrule(+, uA, T2)
    T4, pb4 = rrule(sum, T3)
    T5, pb5 = rrule(abs2, T4)

    dT5 = 1.0
    _, dT4 = pb5(dT5)
    _, dT3 = pb4(dT4)
    _, duA, dT2 = pb3(dT3)
    _, duB, dT1 = pb2(dT2)
    _, duC, _, _ = pb1(dT1)

    @test unthunk(duA) ≈ A.record.grad
    @test unthunk(duB) ≈ B.record.grad
    @test unthunk(duC) ≈ C.record.grad
end


text:
# ======================= SHOWING ====================

Base.show(io::IO, x::PatVar) = print(io, x.name)

function Base.show(io::IO, x::PatLiteral)
    if x.val isa Symbol
        print(io, ":")
    end
    print(io, x.val)
end

Base.show(io::IO, x::PatTypeAssertion) = print(io, x.var, "::", x.type)

Base.show(io::IO, x::PatSplatVar) = print(io, x.name, "...")

Base.show(io::IO, x::PatEquiv) = print(io, x.left, "≡ₙ", x.right)

function Base.show(io::IO, x::PatTerm)
    print(io, Expr(x.head, x.args...))

    # show(io, Expr(x.head, x.args...))
    # n = length(x.args)
    # if x.head isa Symbol
    #     if Base.isbinaryoperator(x.head) && n == 2
    #         print(io, "(", x.args[1], x.head, x.args[2], ")")
    #         return
    #     elseif Base.isunaryoperator(x.head) && n == 1
    #         print(io, "(", x.head, x.args[1], ")")
    #         return
    #     end
    # end

    # print(io, x.head)
    # print(io, "(")
    # for i ∈ 1:n
    #     @inbounds print(io, x.args[i])
    #     if i < n
    #         print(io, ",")
    #     end
    # end
    # print(io, ")")
end

function Base.show(io::IO, x::PatAllTerm)
    n = length(x.args)

    # TODO change me ?
    print(io, "~", x.head)
    print(io, "(")
    for i ∈ 1:n
        @inbounds print(io, x.args[i])
        if i < n
            print(io, ",")
        end
    end
    print(io, ")")
end

# ======================= READING ====================

function Pattern(ex::Expr, mod=@__MODULE__)
    ex = preprocess(ex)
    head = gethead(ex)
    args = getargs(ex)
    meta = getmetadata(ex)

    n = length(args)
    patargs = Vector{Pattern}(undef, n)
    for i ∈ 1:n
        @inbounds patargs[i] = Pattern(args[i], mod)
    end

    if head == :call
        # println(:aaa)
        patargs[1] = PatLiteral(args[1])
    elseif head == :(::)
        v = patargs[1]
        t = patargs[2]
        ty = Union{}
        if t isa PatVar
            ty = getfield(mod, t.name)
        elseif t isa PatLiteral{<:Type}
            ty = t.val
        end
        return PatTypeAssertion(v, ty)
    elseif head == :(...)
        v = patargs[1]
        if v isa PatVar
            return PatSplatVar(v)
        end
    end

    PatTerm(head, patargs, meta)
end

function Pattern(x::Symbol, mod=@__MODULE__)
    PatVar(x)
end

function Pattern(x::QuoteNode, mod=@__MODULE__)
    if x.value isa Symbol
        PatLiteral(x.value)
    else
        PatLiteral(x)
    end
end

# Generic fallback
function Pattern(ex, mod=@__MODULE__)
    ex = preprocess(ex)
    if istree(ex)
        head = gethead(ex)
        args = getargs(ex)
        meta = getmetadata(ex)

        n = length(args)
        patargs = Vector{Pattern}(undef, n)
        for i ∈ 1:n
            @inbounds patargs[i] = Pattern(args[i], mod)
        end

        PatTerm(head, patargs, meta)
    end
    PatLiteral(ex)
end

function Pattern(p::Pattern, mod=@__MODULE__)
    p
end

macro pat(ex)
    Pattern(ex, __module__)
end


text:
module Kaleido

# Use README as the docstring of the module `Kaleido`:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest README")
end Kaleido

export
    @batchlens,
    FLens,
    FlatLens,
    IndexBatchLens,
    KeyBatchLens,
    MultiLens,
    PropertyBatchLens,
    ToField,
    batch,
    constraining,
    converting,
    getting,
    gettingasℝ₊,
    gettingasℝ₋,
    gettingas𝕀,
    nullsetter,
    prettylens,
    setting,
    settingasℝ₊,
    settingasℝ₋,
    settingas𝕀

using Setfield
using Setfield: ComposedLens, IdentityLens, PropertyLens
using Requires

include("base.jl")
include("lensutils.jl")
include("batchsetters.jl")
include("batchlenses.jl")
include("multilens.jl")
include("flatlens.jl")
include("batching.jl")
include("bijection.jl")
include("converterlens.jl")
include("constraininglens.jl")
include("functorlens.jl")
include("setters.jl")
include("prettylens.jl")

function __init__()
    @require(TransformVariables="84d833dd-6860-57f9-a1a7-6da5db126cff",
             include("transformvariables.jl"))
end

end # module


text:
using DataFrames: DataFrame
import CSV

using StatisticalRethinking: sr_datadir, PI
using StatisticalRethinkingCommon

@info "Reading foxes dataset"
@time foxes = CSV.read(sr_datadir("foxes.csv"), DataFrame)

"""
The first two problems are based on the same data. The data in data(foxes)
are 116 foxes from 30 different urban groups in England. These fox groups
are like street gangs. Group size (groupsize) varies from 2 to 8 individuals.
Each group maintains its own (almost exclusive) urban territory. Some ter-
ritories are larger than others. The area variable encodes this information.
Some territories also have more avgfood than others. And food influences
the weight of each fox. Assume this DAG:

A -> F -> G
     |   /
     |  /
     v v
      W

where F is avgfood, G is groupsize, A is area, and W is weight.
Use the backdoor criterion and estimate the total causal influence of A on
F. What effect would increasing the area of a territory have on the amount
of food inside it?
"""
function problem1()
    data = standardize_columns(foxes, (:area, :avgfood))

end

@gen function problem1_model(areas)

    a ~ normal(0, 1)
    b ~ normal(0, 1)

end


text:
"""
    abstract type GNNLayer end

An abstract type from which graph neural network layers are derived.

See also [`GNNChain`](@ref).
"""
abstract type GNNLayer end

# Forward pass with graph-only input.
# To be specialized by layers also needing edge features as input (e.g. NNConv).
(l::GNNLayer)(g::GNNGraph) = GNNGraph(g, ndata=l(g, node_features(g)))

"""
    WithGraph(model, g::GNNGraph; traingraph=false)

A type wrapping the `model` and tying it to the graph `g`.
In the forward pass, can only take feature arrays as inputs,
returning `model(g, x...; kws...)`.

If `traingraph=false`, the graph's parameters, won't be collected
when calling `Flux.params` on a `WithGraph` object.

# Examples

```julia
g = GNNGraph([1,2,3], [2,3,1])
x = rand(Float32, 2, 3)
model = SAGEConv(2 => 3)
wg = WithGraph(model, g)
# No need to feed the graph to `wg`
@assert wg(x) == model(g, x)

g2 = GNNGraph([1,1,2,3], [2,4,1,1])
x2 = rand(Float32, 2, 4)
# WithGraph will ignore the internal graph if fed with a new one.
@assert wg(g2, x2) == model(g2, x2)
```
"""
struct WithGraph{M, G<:GNNGraph}
    model::M
    g::G
    traingraph::Bool
end

WithGraph(model, g::GNNGraph; traingraph=false) = WithGraph(model, g, traingraph)

@functor WithGraph
Flux.trainable(l::WithGraph) = l.traingraph ? (l.model, l.g) : (l.model,)

# Work around
# https://github.com/FluxML/Flux.jl/issues/1733
# Revisit after
# https://github.com/FluxML/Flux.jl/pull/1742
function Flux.destructure(m::WithGraph)
    @assert m.traingraph == false # TODO
    p, re = Flux.destructure(m.model)
    function  re_withgraph(x)
        WithGraph(re(x), m.g, m.traingraph)
    end

    return p, re_withgraph
end

(l::WithGraph)(g::GNNGraph, x...; kws...) = l.model(g, x...; kws...)
(l::WithGraph)(x...; kws...) = l.model(l.g, x...; kws...)

"""
    GNNChain(layers...)
    GNNChain(name = layer, ...)

Collects multiple layers / functions to be called in sequence
on given input graph and input node features.

It allows to compose layers in a sequential fashion as `Flux.Chain`
does, propagating the output of each layer to the next one.
In addition, `GNNChain` handles the input graph as well, providing it
as a first argument only to layers subtyping the [`GNNLayer`](@ref) abstract type.

`GNNChain` supports indexing and slicing, `m[2]` or `m[1:end-1]`,
and if names are given, `m[:name] == m[1]` etc.

# Examples

```juliarepl
julia> m = GNNChain(GCNConv(2=>5), BatchNorm(5), x -> relu.(x), Dense(5, 4));

julia> x = randn(Float32, 2, 3);

julia> g = GNNGraph([1,1,2,3], [2,3,1,1]);

julia> m(g, x)
4×3 Matrix{Float32}:
  0.157941    0.15443     0.193471
  0.0819516   0.0503105   0.122523
  0.225933    0.267901    0.241878
 -0.0134364  -0.0120716  -0.0172505
```
"""
struct GNNChain{T} <: GNNLayer
    layers::T

    GNNChain(xs...) = new{typeof(xs)}(xs)

    function GNNChain(; kw...)
        :layers in Base.keys(kw) && throw(ArgumentError("a GNNChain cannot have a named layer called `layers`"))
        isempty(kw) && return new{Tuple{}}(())
        new{typeof(values(kw))}(values(kw))
    end
end

@forward GNNChain.layers Base.getindex, Base.length, Base.first, Base.last,
    Base.iterate, Base.lastindex, Base.keys

Flux.functor(::Type{<:GNNChain}, c) = c.layers, ls -> GNNChain(ls...)

# input from graph
applylayer(l, g::GNNGraph) = GNNGraph(g, ndata=l(node_features(g)))
applylayer(l::GNNLayer, g::GNNGraph) = l(g)

# explicit input
applylayer(l, g::GNNGraph, x) = l(x)
applylayer(l::GNNLayer, g::GNNGraph, x) = l(g, x)

# Handle Flux.Parallel
applylayer(l::Parallel, g::GNNGraph) = GNNGraph(g, ndata=applylayer(l, g, node_features(g)))
applylayer(l::Parallel, g::GNNGraph, x::AbstractArray) = mapreduce(f -> applylayer(f, g, x), l.connection, l.layers)

# input from graph
applychain(::Tuple{}, g::GNNGraph) = g
applychain(fs::Tuple, g::GNNGraph) = applychain(tail(fs), applylayer(first(fs), g))

# explicit input
applychain(::Tuple{}, g::GNNGraph, x) = x
applychain(fs::Tuple, g::GNNGraph, x) = applychain(tail(fs), g, applylayer(first(fs), g, x))

(c::GNNChain)(g::GNNGraph, x) = applychain(Tuple(c.layers), g, x)
(c::GNNChain)(g::GNNGraph) = applychain(Tuple(c.layers), g)

Base.getindex(c::GNNChain, i::AbstractArray) = GNNChain(c.layers[i]...)
Base.getindex(c::GNNChain{<:NamedTuple}, i::AbstractArray) =
    GNNChain(; NamedTuple{Base.keys(c)[i]}(Tuple(c.layers)[i])...)

function Base.show(io::IO, c::GNNChain)
    print(io, "GNNChain(")
    _show_layers(io, c.layers)
    print(io, ")")
end
_show_layers(io, layers::Tuple) = join(io, layers, ", ")
_show_layers(io, layers::NamedTuple) = join(io, ["$k = $v" for (k, v) in pairs(layers)], ", ")

"""
    DotDecoder()

A graph neural network layer that
for given input graph `g` and node features `x`,
returns the dot product `x_i ⋅ xj` on each edge.

# Examples

```juliarepl
julia> g = rand_graph(5, 6)
GNNGraph:
    num_nodes = 5
    num_edges = 6

julia> dotdec = DotDecoder()
DotDecoder()

julia> dotdec(g, rand(2, 5))
1×6 Matrix{Float64}:
 0.345098  0.458305  0.106353  0.345098  0.458305  0.106353
```
"""
struct DotDecoder <: GNNLayer end

function (::DotDecoder)(g, x)
    check_num_nodes(g, x)
    apply_edges(xi_dot_xj, g, xi=x, xj=x)
end


text:
"""
    Quick MinFlux preprocessing test for MinFlux constraints.
"""

include("testmodel_construction.jl")
m # model
reaction_dG # ΔG'° of reaction dictionary
minflux_bounds = MinFlux_preprocessing(m, reaction_dG)

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
module MapReduce

import Bench.Args

import LinearAlgebra: tr, det, logdet, inv, isposdef

import DataStructures: counter

import Random

import Base.Threads: @threads, nthreads

counter(s::AbstractString) = DataStructures.counter((s,))

const fs  = [ exp, sqrt, tr, det, logdet, inv, isposdef, counter ]

const ops = [ (+), (*), max, merge ]

const benchargs = Args.Bencharg.([
                      n -> map(x -> Random.randstring('a':'d', 2), 1:n),
                      n -> map(x -> Random.randstring('a':'d', 4), 1:n),
                      n -> map(x -> Random.randstring('a':'d', 6), 1:n),
                      n -> Random.randperm(n),
                      n -> convert(Vector{Float16}, Random.randperm(n)),
                      n -> convert(Vector{Float32}, Random.randperm(n)),
                      n -> convert(Vector{Float64}, Random.randperm(n)),
                      n -> convert(Vector{BigInt}, Random.randperm(n)),
                      n -> convert(Vector{BigFloat}, Random.randperm(n)),
                      n -> map(x -> Random.randn(3  , 3  ), 1:n),
                      n -> map(x -> convert(Matrix{BigFloat}, Random.randn(3 , 3 )), 1:n),
                      n -> map(x -> Random.randn(42 , 42 ), 1:n),
                      n -> map(x -> convert(Matrix{BigFloat}, Random.randn(42 , 42 )), 1:n),
                      n -> map(x -> Random.randn(200, 200), 1:n),
                  ])

function threaded_mapreduce(f, op, x)
    @assert length(x) % nthreads() == 0
    results = zeros(eltype(x), nthreads())
    @threads for tid in 1:nthreads()
        # split work
        start = 1 + ((tid - 1) * length(x)) ÷ nthreads()
        stop = (tid * length(x)) ÷ nthreads()
        domain = start:stop

        results[tid] = mapreduce(f, op, view(x, domain))
    end
    foldl(op, results)
end

function atomic_mapreduce(f,op,x)
    results = zeros(eltype(x), nthreads())

    total = Threads.Atomic{Float64}(0.0)
    @threads for tid in 1:nthreads()
        # split work
        start = 1 + ((tid - 1) * length(x)) ÷ nthreads()
        stop = (tid * length(x)) ÷ nthreads()
        domain = (start+1):stop

        acc = f(x[start])
        for j in domain
            @inbounds acc = op(acc, f(x[j]))
        end

        Threads.atomic_add!(total, acc)
    end

    total.value
end

function naive_simd(f, op, x::AbstractArray{T}) where {T}
   s = f(x[1])
   @simd for xᵢ ∈ view(x, 2:length(x))
       s = op(s, f(xᵢ))
   end
   s
end

const mrs = [ mapreduce, threaded_mapreduce, atomic_mapreduce, naive_simd ]

end # module


text:
module RandomBenchmarks

if VERSION >= v"0.7.0-DEV.3406"
    using Random
    using Random: RangeGenerator, GLOBAL_RNG
else
    using Base.Random: RangeGenerator, GLOBAL_RNG
end
using BenchmarkTools
using Compat

const SUITE = BenchmarkGroup()

const BITINTS = [Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128]
const INTS    = [BITINTS; BigInt; Bool]
const FLOATS  = [Float16, Float32, Float64]
const CFLOATS = [Complex{T} for T in FLOATS]
const NUMS    = [INTS; FLOATS; CFLOATS; [Complex{T} for T in BITINTS];]

const MT = MersenneTwister(0)
const RD = RandomDevice()
const vectint = Vector{Int}(undef, 1000)

function set_tolerance!(g, tol=0.25)
    for b in values(g)
        b.params.time_tolerance = tol
    end
end

##########
# Ranges #
##########

g = addgroup!(SUITE, "ranges", ["rand", "rand!", "RangeGenerator"])

g["rand",  "ImplicitRNG",     "Int", "1:1000"] = @benchmarkable rand($(1:1000))
g["rand",  "MersenneTwister", "Int", "1:1000"] = @benchmarkable rand($MT, $(1:1000))
g["rand",  "RandomDevice",    "Int", "1:1000"] = @benchmarkable rand($RD, $(1:1000))
g["rand!", "ImplicitRNG",     "Int", "1:1000"] = @benchmarkable rand!($vectint, $(1:1000))
g["rand!", "MersenneTwister", "Int", "1:1000"] = @benchmarkable rand!($MT, $vectint, $(1:1000))
g["rand!", "RandomDevice",    "Int", "1:1000"] = @benchmarkable rand!($RD, $vectint, $(1:1000))

const b2 = big(2)

for T in INTS
    tstr = string(T)
    onestr = string(one(T))
    for n in T[filter(x -> T === BigInt || x <= typemax(T),
                      [1, 2^32-1, b2^32+1, b2^64-1, b2^64, b2^127, b2^10000]);]
        nstr = n == b2^10000 ? "2^10000" : string(n)
        g["RangeGenerator", tstr, "$onestr:$nstr"] = @benchmarkable RangeGenerator($(T(1):n))
        g["rand", "MersenneTwister", tstr, "RangeGenerator($onestr:$nstr)"] =
            @benchmarkable rand($MT, $(RangeGenerator(T(1):n)))
    end
end

set_tolerance!(g)

#########
# Types #
#########

g = addgroup!(SUITE, "types", ["rand", "rand!", "randn", "randn!", "randexp", "randexp!"])

# these functions are used to interpolate values instead of Types in @benchmarkable,
# to measure the best possible performance, cf. #124
@inline _rand(        ::Val{T}) where {T} = rand(GLOBAL_RNG, T)
@inline _rand(   rng, ::Val{T}) where {T} = rand(rng, T)
@inline _randn(  rng, ::Val{T}) where {T} = randn(rng, T)
@inline _randexp(rng, ::Val{T}) where {T} = randexp(rng, T)

g["rand",    "ImplicitRNG",     "ImplicitFloat64"] = @benchmarkable rand()
g["rand",    "RandomDevice",    "ImplicitFloat64"] = @benchmarkable rand($RD)
g["rand",    "MersenneTwister", "ImplicitFloat64"] = @benchmarkable rand($MT)
g["randn",   "ImplicitRNG",     "ImplicitFloat64"] = @benchmarkable randn()
g["randn",   "RandomDevice",    "ImplicitFloat64"] = @benchmarkable randn($RD)
g["randn",   "MersenneTwister", "ImplicitFloat64"] = @benchmarkable randn($MT)
g["randexp", "ImplicitRNG",     "ImplicitFloat64"] = @benchmarkable randexp()
g["randexp", "RandomDevice",    "ImplicitFloat64"] = @benchmarkable randexp($RD)
g["randexp", "MersenneTwister", "ImplicitFloat64"] = @benchmarkable randexp($MT)

for T = [Int, Float64]
    tstr = string(T)
    dst = Vector{T}(undef, 1000)
    g["rand",     "ImplicitRNG",  tstr] = @benchmarkable _rand($(Val{T}()))
    g["rand",     "RandomDevice", tstr] = @benchmarkable _rand($RD, $(Val{T}()))
    g["rand!",    "ImplicitRNG",  tstr] = @benchmarkable rand!($dst)
    g["rand!",    "RandomDevice", tstr] = @benchmarkable rand!($RD, $dst)
    T === Float64 && VERSION >= v"0.5.0-pre+5657" || continue
    g["randn",    "ImplicitRNG",  tstr] = @benchmarkable randn(Float64)
    g["randn",    "RandomDevice", tstr] = @benchmarkable randn($RD, Float64)
    g["randn!",   "ImplicitRNG",  tstr] = @benchmarkable randn!($dst)
    g["randn!",   "RandomDevice", tstr] = @benchmarkable randn!($RD, $dst)
    g["randexp",  "ImplicitRNG",  tstr] = @benchmarkable randexp(Float64)
    g["randexp",  "RandomDevice", tstr] = @benchmarkable randexp($RD, Float64)
    g["randexp!", "ImplicitRNG",  tstr] = @benchmarkable randexp!($dst)
    g["randexp!", "RandomDevice", tstr] = @benchmarkable randexp!($RD, $dst)
end

for T in NUMS
    T === BigInt && continue
    tstr = string(T)
    dst = Vector{T}(undef, 1000)
    g["rand",     "MersenneTwister", tstr] = @benchmarkable _rand($MT, $(Val{T}()))
    g["rand!",    "MersenneTwister", tstr] = @benchmarkable rand!($MT, $dst)
    VERSION >= v"0.5.0-pre+5657" || continue
    T <: AbstractFloat || T in CFLOATS && VERSION >= v"0.7.0-DEV.973" || continue
    g["randn",    "MersenneTwister", tstr] = @benchmarkable _randn($MT, $(Val{T}()))
    g["randn!",   "MersenneTwister", tstr] = @benchmarkable randn!($MT, $dst)
    T <: AbstractFloat || continue
    g["randexp",  "MersenneTwister", tstr] = @benchmarkable _randexp($MT, $(Val{T}()))
    g["randexp!", "MersenneTwister", tstr] = @benchmarkable randexp!($MT, $dst)
end

set_tolerance!(g)

###############
# Collections #
###############

g = addgroup!(SUITE, "collections", ["rand", "rand!"])

collections = Pair[[1:3;]   => "small Vector",
                   [1:900;] => "large Vector",
                   'a':'z'  => "'a':'z'"]

if VERSION >= v"0.7.0-DEV.973"
    push!(collections,
          Dict(1=>2, 3=>4, 5=>6)  => "small Dict",
          Dict(zip(1:900, 1:900)) => "large Dict",
          Set(1:3)                => "small Set",
          Set(1:900)              => "large Set",
          BitSet(1:3)             => "small BitSet",
          BitSet(1:900)           => "large BitSet",
          "qwèrtï"                => "small String",
          randstring(900)         => "large String")
end

for (collection, collstr) in collections
    dst = Vector{eltype(collection)}(undef, 1000)
    g["rand",  "ImplicitRNG",     collstr] = @benchmarkable rand($collection)
    g["rand",  "MersenneTwister", collstr] = @benchmarkable rand($MT, $collection)
    g["rand",  "RandomDevice",    collstr] = @benchmarkable rand($RD, $collection)
    g["rand!", "ImplicitRNG",     collstr] = @benchmarkable rand!($dst, $collection)
    g["rand!", "MersenneTwister", collstr] = @benchmarkable rand!($MT, $dst, $collection)
    g["rand!", "RandomDevice",    collstr] = @benchmarkable rand!($RD, $dst, $collection)
end

set_tolerance!(g)

##############
# randstring #
##############

g = addgroup!(SUITE, "randstring")

qwerty, qwertystr = collect(UInt8, "qwerty"), "collect(UInt8, \"qwerty\""
g["randstring", "MersenneTwister"]                        = @benchmarkable randstring($MT)
g["randstring", "MersenneTwister", 100]                   = @benchmarkable randstring($MT, 100)
if VERSION >= v"0.7.0-DEV.973"
    g["randstring", "MersenneTwister", "\"qwèrtï\""]      = @benchmarkable randstring($MT, "qwèrtï")
    g["randstring", "MersenneTwister", "\"qwèrtï\"", 100] = @benchmarkable randstring($MT, "qwèrtï", 100)
    g["randstring", "MersenneTwister", qwertystr]         = @benchmarkable randstring($MT, $qwerty)
    g["randstring", "MersenneTwister", qwertystr, 100]    = @benchmarkable randstring($MT, $qwerty, 100)
end

set_tolerance!(g)

#############
# sequences #
#############

g = addgroup!(SUITE, "sequences", ["randsubseq!", "shuffle!", "randperm", "randcycle"])

# for randsubseq!, vectint is a big enough vector as output arg to avoid re-allocations
src = [1:1000;]
g["randsubseq!", "MersenneTwister", "0.2"]  = @benchmarkable randsubseq!($MT, $vectint, $src, 0.2)
g["randsubseq!", "MersenneTwister", "0.8"]  = @benchmarkable randsubseq!($MT, $vectint, $src, 0.8)
g["shuffle!",    "MersenneTwister"]         = @benchmarkable shuffle!($MT, $src)
g["randperm",    "MersenneTwister", "5"]    = @benchmarkable randperm($MT, 5)
g["randperm",    "MersenneTwister", "1000"] = @benchmarkable randperm($MT, 1000)
g["randcycle",   "MersenneTwister", "5"]    = @benchmarkable randcycle($MT, 5)
g["randcycle",   "MersenneTwister", "1000"] = @benchmarkable randcycle($MT, 1000)

set_tolerance!(g)

end # module


text:
const include_regex = Regex(raw"""^#include [<"](([A-Z][a-z0-9]+)+\.(hpp|inl))[">]$""", "m")

const header = """
/*

Juliapp
=======
Seamlessly embed Julia code in C++ with an API using modern language features.
https://github.com/kim366/Juliapp

=======

$(read("LICENSE", String))
*/

#ifndef JULIAPP_SINGLE_INCLUDE_HPP
#define JULIAPP_SINGLE_INCLUDE_HPP"""

const footer = """

#endif // JULIAPP_SINGLE_INCLUDE_HPP
"""

readheader(filename) = replace(read("src/$filename", String), "#pragma once" => "")
stripnewlines(content) = replace(content, Regex(raw"^\n\n+$", "m") => "\n")

let sourcecode = readheader("Julia.hpp")
  included_headers = []

  while true
    current_match = match(include_regex, sourcecode)
    current_match == nothing && break
    filename_to_include = current_match.captures[1]

    if filename_to_include ∉ included_headers
      content_to_substitute = readheader(filename_to_include)
      push!(included_headers, filename_to_include)
    else
      content_to_substitute = ""
    end

    sourcecode = sourcecode[1:current_match.offset - 1] *
      content_to_substitute *
      sourcecode[current_match.offset + length(current_match.match):end]
  end

  write("single_include/juliapp/julia.hpp", header * stripnewlines(sourcecode) * footer)
end

max_stars_repo_name: cdsousa/Juliapp

max_issues_repo_name: cdsousa/Juliapp

max_forks_repo_name: cdsousa/Juliapp


text:
import .Zygote

function ∂ℓπ∂θ_zygote(ℓπ, θ::AbstractVector)
    res, back = Zygote.pullback(ℓπ, θ)
    return res, first(back(Zygote.sensitivity(res)))
end

function ∂ℓπ∂θ_zygote(ℓπ, θ::AbstractMatrix)
    res, back = Zygote.pullback(ℓπ, θ)
    return res, first(back(ones(eltype(res), size(res))))
end

function ZygoteADHamiltonian(metric::AbstractMetric, ℓπ)
    ∂ℓπ∂θ(θ::AbstractVecOrMat) = ∂ℓπ∂θ_zygote(ℓπ, θ)
    return Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
end

ADAVAILABLE[Zygote] = ZygoteADHamiltonian


text:
function part1(p1, p2, dice=100, goal=1000)
    s1, s2 = 0, 0
    die_value, rolls = 0, 0
    die() = rolls += 1; mod1(die_value += 1, dice)
    roll() = die() + die() + die()
    while s1 < goal || s2 < goal
        s1 += (p1 = mod1(p1 + roll(), 10))
        s1 >= goal && break
        s2 += (p2 = mod1(p2 + roll(), 10))
    end
    return min(s1, s2) * rolls
end

function part2(P1, P2, goal=22)  # Add 1 to goal, since julia is 1 indexed, so dp lookup is different
    dp = zeros(Int64, (goal+11, goal+11, 10, 10, 2))
    dp[1, 1, P1, P2, 1] = 1
    dice = [sum(1 for a=1:3, b=1:3, c=1:3 if a+b+c==s) for s=3:9]
    for s1=1:goal-1, s2=1:goal-1
        for d1=1:10, d2=1:10, add=3:9
            d1a = mod1(d1+add, 10)
            dp[s1+d1a, s2, d1a, d2, 2] += dp[s1, s2, d1, d2, 1] * dice[add-2]
            d2a = mod1(d2+add, 10)
            dp[s1, s2+d2a, d1, d2a, 1] += dp[s1, s2, d1, d2, 2] * dice[add-2]
        end
    end
    dps = sum(dp, dims=(3, 4, 5))
    S1, S2 = 0, 0
    for s1=1:goal+11, s2=1:goal+11
        if max(s1, s2) >= goal
            s1 > s2 ? S1 += dps[s1, s2, 1, 1, 1] : S2 += dps[s1, s2, 1, 1, 1]
        end
    end
    return max(S1, S2)
end

function main()
    p1, p2 = parse.(Int, getindex.(split.(readlines(), ": "), 2))
    println(part1(p1, p2))
    println(part2(p1, p2))
end
main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode