text:
const directions_2d = [
    :x, :y, :xy, :yx
]

const directions_3d = [
    :x, :y, :z,
    :xy, :yx,
    :xz, :zx,
    :yz, :zy,
    :xyz, :yxz, :xzy, :zyx
]

const trackers = [
    S2Tracker(false),
    L2Tracker(false),
    L2Tracker(true),
    SSTracker(false),
    SVTracker(false),
    CCTracker(false, true),
    CCTracker(true, false)
]

create_tracker(array, trackers, periodic, directions) =
    CorrelationTracker(array;
                       tracking   = trackers,
                       periodic   = periodic,
                       directions = directions)

function test_tracker!(array    :: AbstractArray{T},
                       trackers :: Vector{<:AbstractTracker{T}},
                       periodic :: Bool,
                       directions) where T
    indices = CartesianIndices(array)
    tracker = create_tracker(array, trackers, periodic, directions)
    for n in 1:100
        idx = rand(indices)
        array[idx]   = 1 - array[idx]
        tracker[idx] = 1 - tracker[idx]
    end

    for data in tracked_data(tracker)
        expected = data(array;
                        periodic   = periodic,
                        directions = directions)
        got = data(tracker)
        for (direction, data) in expected
            @test data ≈ got[direction]
        end
    end
    @test tracker == array
end

function test_rollback!(array    :: AbstractArray{T},
                        trackers :: Vector{<:AbstractTracker{T}},
                        periodic :: Bool,
                        directions) where T
    indices = CartesianIndices(array)
    tracker = create_tracker(array, trackers, periodic, directions)

    for n in 1:100
        idx = rand(indices)
        token = update_corrfns!(tracker, 1 - tracker[idx], idx)
        if rand() < 0.5
            rollback!(tracker, token)
        else
            array[idx] = 1 - array[idx]
        end
    end

    for data in tracked_data(tracker)
        expected = data(array;
                        periodic   = periodic,
                        directions = directions)
        got = data(tracker)
        for (direction, data) in expected
            @test data ≈ got[direction]
        end
    end
    @test tracker == array
end

@testset "2D system" begin
    test_tracker!(rand(Bool, (100, 50)),  trackers, false, directions_2d)
    test_tracker!(rand(Bool, (100, 100)), trackers, true,  directions_2d)
end

@testset "3D system" begin
    test_tracker!(rand(Bool, (100, 50,  200)), trackers, false, directions_3d)
    test_tracker!(rand(Bool, (100, 100, 100)), trackers, true,  directions_3d)
end

@testset "2D system (rollback)" begin
    test_rollback!(rand(Bool, (100, 50)),  trackers, false, directions_2d)
    test_rollback!(rand(Bool, (100, 100)), trackers, true,  directions_2d)
end

@testset "3D system (rollback)" begin
    test_rollback!(rand(Bool, (100, 50,  200)), trackers, false, directions_3d)
    test_rollback!(rand(Bool, (100, 100, 100)), trackers, true,  directions_3d)
end

@testset "Extrapolations" begin
    # TODO: test diagonals
    array = rand(Bool, (200, 200))
    tracker = CorrelationTracker(array; periodic = true)
    extra = ExtrapolatedData(tracker, 3, [:x, :y, :z])

    s21 = Directional.s2(tracker, false)
    s22 = Directional.s2(extra, false)

    @test s22[:x] ≈ s21[:x]
    @test s22[:y] ≈ s21[:y]
    @test s22[:z] ≈ (s21[:x] + s21[:y]) / 2
end


text:
#=
triangle.jl
2D triangle shape
=#

#using ImagePhantoms #: Object, Object2d

export Triangle
export phantom, radon, spectrum

const sqrt3 = sqrt(3)

"""
    Triangle{T} <: AbstractShape2
By default an equilateral triangle pointing upward with the "center"
in the middle of the base, for the default parameter = `0.5`.

For parameter `p` the object is a triangle whose base is along the x-axis
going from (-p,0) to (1-p,0)` and with height=sqrt(3)/2.

Most methods currently support only the case `p=0.5`.
"""
struct Triangle <: AbstractShape2 end

#=
    Triangle{T} <: AbstractShape2
The (width and height can be scaled when constructing an object).
struct Triangle{T} <: AbstractShape2
    param::T # fraction in interval (0,1)
    function Triangle{T}(param::T = 0.5) where {T <: Real}
        0 < param < 1 || throw(ArgumentError("param=$param"))
        new{T}(param)
    end
end
=#

# constructors

"""
    Triangle(cx, cy, wx=1, wy=wx, ϕ=0, value::Number=1, param::Real=0.5)
    Triangle(center::NTuple{2,RealU}, width::NTuple{2,RealU}=(1,1), ϕ::RealU=0, v=1, param=0.5)
    Triangle([6-vector] or [7-vector])
Construct `Triangle` object from parameters.
In the typical case where `param=0.5` and `width[1] == width[2]`,
this is an equilateral triangle with base `width[1]` centered along the x axis.
"""
function Triangle(
    cx::RealU,
    cy::RealU,
    wx::RealU = oneunit(cx),
    wy::RealU = wx,
    ϕ::RealU = 0,
    value::Number = 1,
    param::Real = 0.5,
)
    (cx, cy, wx, wy) = promote(cx, cy, wx, wy)
    Object(Triangle(), (cx,cy), (wx,wy), ϕ, value, param)
end

function Triangle(
    center::NTuple{2,RealU},
    width::NTuple{2,RealU} = (1,1) .* oneunit(center[1]),
    ϕ::RealU = 0,
    value::Number = 1,
    param::Real = 0.5,
)
    Triangle(center..., width..., ϕ, value, param)
end

function Triangle(v::AbstractVector{<:Number})
    length(v) ∈ (6,7) || throw(ArgumentError("$v wrong length"))
    Triangle(v...)
end

#=
"""
    shape = Triangle(param=0.5)
"""
Triangle(param::T=0.5) where {T <: Real} = Triangle{T}(param)
=#

# helpers

function _trifun(x, y, param)
    param == 1/2 || throw("todo")
    return (0 ≤ y ≤ sqrt3/2) && y ≤ sqrt3 * (0.5 - abs(x))
end

"""
    (lower, upper) = _interval(a, b)
Determine the interval `[lower, upper]` corresponding to the set
`{x : b ≤ a x}` where `a` is unitless but `b` and `x` have same units.
"""
function  _interval(a::Real, b::RealU)
    Infu = Inf * oneunit(b)
    return a > 0 ? (b/a, Infu) :
        a < 0 ? (-Infu, b/a) :
        b ≤ zero(b) ? (-Infu, Infu) : # a == 0
        (T = typeof(oneunit(b)/one(a)); (zero(T),zero(T)))
end

"""
    radon_tri(r, sinϕ, cosϕ)

For a line integral at radial position `r` and angle `ϕ`,
the locus of points along the line is
`{(r cos(ϕ), r sin(ϕ)) + τ (-sin(ϕ), cos(ϕ)) : τ ∈ ℝ}`.
This function treats the equilateral triangle
with base [-1/2,1/2], pointing upwards,
as the intersection of three half planes:
* H0 = `{(x,y) : y ≥ 0}`
* H1 = `y ≤ √3 (1/2 - x)`
* H2 = `y ≤ √3 (1/2 + x)`.
Find the `τ` values where the line locus lies in each half planes,
then take the length of the intersection of those three intervals.

For example, for H1 we have
`r sin(ϕ) + τ cos(ϕ) ≤ √3 (1/2 - (r cos(ϕ) - τ sin(ϕ)))`
or equivalently
`r (sin(ϕ) + √3 cos(ϕ)) - √3/2 ≤ τ (√3 sin(ϕ) - cos(ϕ))`
which is a set the form `b1 ≤ a1 τ`, corresponding to some interval `(l1,u1)`.
Similarly for H0 and H2.

This approach might not be the most efficient, but it is simple.

See Peter Aundal Toft, "The Radon transform - theory and implementation", 1996
https://orbit.dtu.dk/en/publications/the-radon-transform-theory-and-implementation
for a different approach to finding the Radon transform of a triangle.
"""
function radon_tri(r, sinϕ, cosϕ)
    (l1, u1) = _interval( sqrt3*sinϕ - cosϕ, r*(sinϕ + sqrt3*cosϕ) - sqrt3/2)
    (l2, u2) = _interval(-sqrt3*sinϕ - cosϕ, r*(sinϕ - sqrt3*cosϕ) - sqrt3/2)
    (l0, u0) = _interval(cosϕ, -r*sinϕ)
    tmp = min(u0, u1, u2) - max(l0, l1, l2)
    return tmp > 0 ? tmp : zero(tmp)
end

# methods

"""
    phantom(ob::Object2d{Triangle})
Returns function of `(x,y)` for making image.
"""
phantom(ob::Object2d{Triangle}) = (x,y) ->
    ob.value * _trifun(coords(ob, x, y)..., ob.param)

"""
    radon_tri(r, ϕ, cx, cy, wx, wy, θ, p)
Radon transform at `(r,ϕ)` of triangle.
"""
function radon_tri(r, ϕ, cx::C, cy::C, wx::C, wy::C, θ, p) where {C <: RealU}
    T = promote_type(eltype(r), C)
    p == 1/2 || throw("todo")
    (sinϕ, cosϕ) = sincos(ϕ)
    r -= cx * cosϕ + cy * sinϕ # Radon translation property
    ϕ -= θ # Radon rotation property
    (sinϕ, cosϕ) = sincos(ϕ)
    # Radon affine scaling property
    denom = sqrt(abs2(cosϕ/wy) + abs2(sinϕ/wx))
    (sinϕ, cosϕ) = sincos(ϕ)
    r /= wx * wy * denom # unitless
    ϕ = atan(sinϕ/wx, cosϕ/wy)
    return abs(r) ≥ sqrt(3)/2 ? zero(T) : T(radon_tri(r, sincos(ϕ)...) / denom)
end

"""
    radon(ob::Object2d{Triangle})
Returns function of `(r,ϕ)` for making a sinogram.
"""
radon(ob::Object2d{Triangle}) = (r,ϕ) -> ob.value *
    radon_tri(r, ϕ, ob.center..., ob.width..., ob.angle[1], ob.param)

expimp(x::Real) = exp(1im*π * x)

function spectrum_tri(u, v)
    (u == 0 && v == 0) ? sqrt3 / 4 :
    (u == 0) ? 1im / (2π * v) * (expimp(-v*sqrt3) - 1) + 1/(2*sqrt3*abs2(π*v)) *
        (1 - expimp(-v * sqrt3) * (1im * π * v * sqrt3 + 1)) :
    im*sqrt3/(4π*u) * expimp(-sqrt3/2*v) * (
        expimp(-u/2) * sinc((v*sqrt3 - u)/2) -
        expimp( u/2) * sinc((v*sqrt3 + u)/2) )
end

function spectrum_tri(fx, fy, cx, cy, wx, wy, θ, param)
    param == 1/2 || throw("todo")
    (kx, ky) = rotate2d(fx, fy, θ) # rotate first, then translate
    return wx * exp(-2im*π*fx*cx) *
           wy * exp(-2im*π*fy*cy) * spectrum_tri(kx*wx, ky*wy)
end

"""
    spectrum(ob::Object2d{Triangle})
Returns function of ``(f_x,f_y)`` for the spectrum (2D Fourier transform).
"""
spectrum(ob::Object2d{Triangle}) = (fx,fy) -> ob.value *
    spectrum_tri(fx, fy, ob.center..., ob.width..., ob.angle[1], ob.param)

ext: jl


text:
# test pd matrix types
using PDMats
using Base.Test

call_test_pdmat(p::AbstractPDMat,m::Matrix) = test_pdmat(p,m,cmat_eq=true,verbose=1)

for T in [Float64,Float32]
  #test that all external constructors are accessible
  m = eye(T,2) ; @test PDMat(m,cholfact(m)).mat == PDMat(Symmetric(m)).mat == PDMat(m).mat == PDMat(cholfact(m)).mat
  d = ones(T,2) ; @test PDiagMat(d,d).inv_diag == PDiagMat(d).inv_diag
  x = one(T) ; @test ScalMat(2,x,x).inv_value == ScalMat(2,x).inv_value
  #for Julia v0.3, do not test sparse matrices with Float32 (see issue #14076)
  (T == Float64 || VERSION >= v"0.4.2") && (s = speye(T,2,2) ; @test PDSparseMat(s,cholfact(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholfact(s)).mat)

  #test the functionality
  M = convert(Array{T,2},[4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
  V = convert(Array{T,1},[1.5, 2.5, 2.0])
  X = convert(T,2.0)

  call_test_pdmat(PDMat(M),M) #tests of PDMat
  call_test_pdmat(PDiagMat(V),diagm(V)) #tests of PDiagMat
  call_test_pdmat(ScalMat(3,x),x*eye(T,3)) #tests of ScalMat
  #for Julia v0.3, do not test sparse matrices with Float32 (see issue #14076)
  (T == Float64 || VERSION >= v"0.4.2") && call_test_pdmat(PDSparseMat(sparse(M)),M)
end

m = eye(Float32,2)
@test convert(PDMat{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
m = ones(Float32,2)
@test convert(PDiagMat{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
x = one(Float32); d = 4
@test convert(ScalMat{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
if VERSION >= v"0.4.2"
    s = speye(Float32, 2, 2)
    @test convert(PDSparseMat{Float64}, PDSparseMat(s)).mat == PDSparseMat(convert(SparseMatrixCSC{Float64}, s)).mat
end


text:
@testset "Extensive2Matrix" begin
    game = CFR.IIEMatrixGame()
    mat = Matrix(game)

    @test all(
        all(mat[i] .≈ game.R[i]) for i in eachindex(mat, game.R)
    )

    game = CFR.IIEMatrixGame([
        (rand(),rand()) ((rand(),rand()))
        (rand(),rand()) ((rand(),rand()))
    ])
    mat = Matrix(game)

    @test all(
        all(mat[i] .≈ game.R[i]) for i in eachindex(mat, game.R)
    )
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SmoothingMethod

A method for smoothing meshes.
"""
abstract type SmoothingMethod end

"""
    smooth(mesh, method)

Smooth `mesh` with given method.
"""
function smooth end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("smoothing/taubin.jl")


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

function tohtml(io::IO, m::MIME"text/html", x)
    show(io, m, x)
end

function tohtml(io::IO, m::MIME"text/plain", x)
    htmlesc(io, sprint(show, m, x))
end

function tohtml(io::IO, m::MIME"image/png", img)
    print(io, """<img src="data:image/png;base64,""")
    print(io, stringmime(m, img))
    print(io, "\" />")
end

function tohtml(m::MIME"image/svg+xml", img)
    show(io, m, img)
end

# Display infrastructure

function bestmime(val)
    for mime in ("text/html", "image/svg+xml", "image/png", "text/plain")
        mimewritable(mime, val) && return MIME(Symbol(mime))
    end
    error("Cannot render $val to Markdown.")
end

tohtml(io::IO, x) = tohtml(io, bestmime(x), x)


text:
using Reexport

include("Meta.jl")
include("GAT.jl")
include("Syntax.jl")
include("Rewrite.jl")
include("Present.jl")

@reexport using .GAT
@reexport using .Syntax
@reexport using .Rewrite
@reexport using .Present


text:
using ThreeBodyDecay
using Test
using Parameters

@testset "BasicParityRecoupling" begin
    reaction = jp"1/2+"=>(jp"1/2+",jp"0-")

    H_pc = ParityRecoupling(1,0, jp"1/2+"=>(jp"1/2+",jp"0-"))
    @test H_pc == ParityRecoupling(1,0,'-')

    H_ls = RecouplingLS(possible_ls(reaction)[1] |> x2, reaction)

    factor = -1 / sqrt(2)
    @test amplitude(H_pc, 1, 0) * factor ≈ amplitude(H_ls, 1, 0)
    @test amplitude(H_pc,-1, 0) * factor ≈ amplitude(H_ls,-1, 0)
end

@testset "DecayChainCouplings" begin
    mΛb = 5.61960
    mπ = 0.14
    # intermediate resonances
    mΣb = 5.81065; ΓΣb = 0.005
    mΣb_x = 5.83032; ΓΣb_x = 0.0094

    mΛb2S = 6.3; # just a peak of the plot

    tbs = ThreeBodySystem(mπ,mΛb,mπ; m0=mΛb2S,
        two_js = ThreeBodySpins(0, 1, 0; two_h0=1))  # (0- 1/2+ 0- | 1/2+)

    dc_pc = DecayChain(
        k=1,
        Xlineshape=σ->BW(σ,mΣb,ΓΣb),
        tbs=tbs, two_s = 1,
        Hij=ParityRecoupling(1,0,'-'),
        HRk=ParityRecoupling(1,0,'-'))
    #
    σs = randomPoint(tbs.ms)
    #
    @test amplitude(σs, [0, 1,0, 1], dc_pc) != 0
    @test amplitude(σs, [0, 1,0,-1], dc_pc) != 0
    @test amplitude(σs, [0,-1,0, 1], dc_pc) != 0
    @test amplitude(σs, [0,-1,0,-1], dc_pc) != 0
    #
    dc_pv = DecayChain(
        k=1,
        Xlineshape=σ->BW(σ,mΣb,ΓΣb),
        tbs=tbs, two_s = 1,
        Hij=ParityRecoupling(1,0,'-'),
        HRk=NoRecoupling(1,0))
    #
    @test amplitude(σs, [0, 1,0, 1], dc_pv) != 0
    @test amplitude(σs, [0,-1,0, 1], dc_pv) != 0
    @test amplitude(σs, [0,-1,0,-1], dc_pv) == 0im
    @test amplitude(σs, [0,-1,0,-1], dc_pv) == 0im
    #
    dc_pv = DecayChain(
        k=1,
        Xlineshape=σ->BW(σ,mΣb,ΓΣb),
        tbs=tbs, two_s = 1,
        Hij=ParityRecoupling(1,0,'-'),
        HRk=NoRecoupling(1,0))
    #
    @test amplitude(σs, [0, 1,0, 1], dc_pv) != 0
    @test amplitude(σs, [0,-1,0, 1], dc_pv) != 0
    @test amplitude(σs, [0,-1,0,-1], dc_pv) == 0im
    @test amplitude(σs, [0,-1,0,-1], dc_pv) == 0im
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
## 01.jl
# 文字列の切り出し

str = "パタトクカシーー"
println(join(split(str,"")[1:2:end]))


text:
# Preliminary file to think through David's suggestion of splitting defs between "registered"
# readonly "templates"  and user-constructed models so it's clear which functions operate on
# templates vs defs within a model.

@class ComponentDef <: NamedObj begin
    comp_id::Union{Nothing, ComponentId}
    variables::OrderedDict{Symbol, VariableDef}
    parameters::OrderedDict{Symbol, ParameterDef}
    dim_names::Set{Symbol}
end

@class CompositeComponentDef <: ComponentDef begin
    comps_dict::OrderedDict{Symbol, AbstractComponentDef}
    exports::ExportsDict
    internal_param_conns::Vector{InternalParameterConnection}
    external_params::Dict{Symbol, ModelParameter}
end

# Define these for building out a ModelDef, which reference the
# central definitions using the classes above.

@class mutable ModelComponentDef <: NamedObj begin
    comp_id::ComponentId    # references registered component def
    comp_path::Union{Nothing, ComponentPath}
    dim_dict::OrderedDict{Symbol, Union{Nothing, Dimension}}
    first::Union{Nothing, Int}
    last::Union{Nothing, Int}
    is_uniform::Bool
end

@class mutable ModelCompositeComponentDef <: ModelComponentDef begin
    comps_dict::OrderedDict{Symbol, AbstractModelComponentDef}
    bindings::Vector{Binding}
    exports::ExportsDict
    external_param_conns::Vector{ExternalParameterConnection}
    external_params::Dict{Symbol, ModelParameter}

    # Names of external params that the ConnectorComps will use as their :input2 parameters.
    backups::Vector{Symbol}

    sorted_comps::Union{Nothing, Vector{Symbol}}
end


text:
using OrdinaryDiffEq, StochasticDiffEq, Flux, DiffEqFlux, Test

x = Float32[2.; 0.]
tspan = (0.0f0,25.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2))

neural_ode(dudt,x,tspan,Tsit5(),save_everystep=false,save_start=false)
neural_ode(dudt,x,tspan,Tsit5(),saveat=0.1)
neural_ode_rd(dudt,x,tspan,Tsit5(),saveat=0.1)

Flux.back!(sum(neural_ode(dudt,x,tspan,Tsit5(),save_everystep=false,save_start=false)))

# Adjoint
@testset "adjoint mode" begin
    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode(dudt,x,tspan,Tsit5(),save_everystep=false,save_start=false)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode(dudt,x,tspan,Tsit5(),saveat=0.0:0.1:10.0)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    @test_broken Flux.back!(sum(neural_ode(dudt,x,tspan,Tsit5(),saveat=0.1)))
    # @test ! iszero(Tracker.grad(dudt[1].W))
end;

# RD
@testset "reverse mode" begin
    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode_rd(dudt,x,tspan,Tsit5(),save_everystep=false,save_start=false)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode_rd(dudt,x,tspan,Tsit5(),saveat=0.0:0.1:10.0)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode_rd(dudt,x,tspan,Tsit5(),saveat=0.1)))
    @test ! iszero(Tracker.grad(dudt[1].W))
end;

mp = Float32[0.1,0.1]
Tracker.zero_grad!(dudt[1].W.grad)
neural_dmsde(dudt,x,mp,(0.0f0,2.0f0),SOSRI(),saveat=0.1)
Flux.back!(sum(neural_dmsde(dudt,x,mp,(0.0f0,2.0f0),SOSRI(),saveat=0.0:0.1:2.0)))
@test ! iszero(Tracker.grad(dudt[1].W))

# Batch
xs = Float32.(hcat([0.; 0.], [1.; 0.], [2.; 0.]))
tspan = (0.0f0,25.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2))

neural_ode(dudt,xs,tspan,Tsit5(),save_everystep=false,save_start=false)
neural_ode(dudt,xs,tspan,Tsit5(),saveat=0.1)
neural_ode_rd(dudt,xs,tspan,Tsit5(),saveat=0.1)

# Adjoint
@testset "adjoint mode batches" begin
    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode(dudt,xs,tspan,Tsit5(),save_everystep=false,save_start=false)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode(dudt,xs,tspan,Tsit5(),saveat=0.0:0.1:10.0)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    @test_broken Flux.back!(sum(neural_ode(dudt,xs,tspan,Tsit5(),saveat=0.1)))
    #@test ! iszero(Tracker.grad(dudt[1].W))
end;

# RD
@testset "reverse mode batches" begin
    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode_rd(dudt,xs,tspan,Tsit5(),save_everystep=false,save_start=false)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode_rd(dudt,xs,tspan,Tsit5(),saveat=0.0:0.1:10.0)))
    @test ! iszero(Tracker.grad(dudt[1].W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Flux.back!(sum(neural_ode_rd(dudt,xs,tspan,Tsit5(),saveat=0.1)))
    @test ! iszero(Tracker.grad(dudt[1].W))
end;

# Grads w.r.t. u0
tspan = (0.0f0,25.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2))
downsample = Dense(3,2)
x0 = Float32.(hcat([0.; 0; 0.], [1.,1.,1.],[2.,2.,2.]))
u0 = downsample(x0)

neural_ode(dudt,u0,tspan,Tsit5(),save_everystep=false,save_start=false)
neural_ode(dudt,u0,tspan,Tsit5(),saveat=0.1)
neural_ode_rd(dudt,u0,tspan,Tsit5(),saveat=0.1)

# Adjoint
@testset "adjoint mode trackedu0" begin
    Tracker.zero_grad!(dudt[1].W.grad)
    Tracker.zero_grad!(downsample.W.grad)
    m1 = Chain(downsample, u0->neural_ode(dudt,u0,tspan,Tsit5(),save_everystep=false,save_start=false))
    Flux.back!(sum(m1(x0)))
    @test ! iszero(Tracker.grad(dudt[1].W))
    @test ! iszero(Tracker.grad(downsample.W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Tracker.zero_grad!(downsample.W.grad)
    m2 = Chain(downsample, u0->neural_ode(dudt,u0,tspan,Tsit5(),saveat=0.0:0.1:10.0))
    Flux.back!(sum(m2(x0)))
    @test ! iszero(Tracker.grad(dudt[1].W))
    @test ! iszero(Tracker.grad(downsample.W))

    Tracker.zero_grad!(dudt[1].W.grad)
    Tracker.zero_grad!(downsample.W.grad)
    m3 = Chain(downsample, u0->neural_ode(dudt,u0,tspan,Tsit5(),saveat=0.1))
    @test_broken Flux.back!(sum(m3(x0)))
    #@test ! iszero(Tracker.grad(dudt[1].W))
    #@test ! iszero(Tracker.grad(downsample.W))
end;

# RD
@testset "reverse mode trackedu0" begin
    Tracker.zero_grad!(dudt[1].W.grad)
    u0 = downsample(x0)
    Flux.back!(sum(neural_ode_rd(dudt,u0,tspan,Tsit5(),save_everystep=false,save_start=false)))
    @test ! iszero(Tracker.grad(dudt[1].W))
    @test ! iszero(Tracker.grad(u0))

    Tracker.zero_grad!(dudt[1].W.grad)
    u0 = downsample(x0)
    Flux.back!(sum(neural_ode_rd(dudt,u0,tspan,Tsit5(),saveat=0.0:0.1:10.0)))
    @test ! iszero(Tracker.grad(dudt[1].W))
    @test ! iszero(Tracker.grad(u0))

    Tracker.zero_grad!(dudt[1].W.grad)
    u0 = downsample(x0)
    @test_broken Flux.back!(sum(neural_ode_rd(dudt,u0,tspan,Tsit5(),saveat=0.1)))
    #@test ! iszero(Tracker.grad(dudt[1].W))
    #@test ! iszero(Tracker.grad(downsample.W))
end;

@testset "struct NeuralODE" begin
    using BenchmarkTools
    no = DiffEqFlux.NeuralODE(dudt,tspan)
    struct_sol = no(u0,Tsit5(),save_everystep=false,save_start=false)
    func_sol = neural_ode(dudt,u0,tspan,Tsit5(),save_everystep=false,save_start=false)
    @test struct_sol == func_sol
end;


text:
@testset "Limited Length Cycles" begin
    completedg = complete_digraph(4)
    pathdg = path_digraph(5)
    cycledg = cycle_digraph(5)

    @testset "complete digraph" for g in testgraphs(completedg)
        @test length(simplecycles_limited_length(g, 0)) == 0
        @test length(simplecycles_limited_length(g, 1)) == 0
        @test length(simplecycles_limited_length(g, 2)) == 6
        @test length(simplecycles_limited_length(g, 3)) == 14
        @test length(simplecycles_limited_length(g, 4)) == 20
        @test length(simplecycles_limited_length(g, 4, 10)) == 10
        @test length(simplecycles_limited_length(g, 4, typemax(Int))) == 20
    end

    @testset "path digraph" for g in testgraphs(pathdg)
        @test length(simplecycles_limited_length(g, 1)) == 0
        @test length(simplecycles_limited_length(g, 2)) == 0
        @test length(simplecycles_limited_length(g, 3)) == 0
        @test length(simplecycles_limited_length(g, 4)) == 0
        @test length(simplecycles_limited_length(g, 5)) == 0
    end

    @testset "cycle digraph" for g in testgraphs(cycledg)
        @test length(simplecycles_limited_length(g, 1)) == 0
        @test length(simplecycles_limited_length(g, 2)) == 0
        @test length(simplecycles_limited_length(g, 3)) == 0
        @test length(simplecycles_limited_length(g, 4)) == 0
        @test length(simplecycles_limited_length(g, 5)) == 1
    end

    @testset "self loops" begin
        selfloopg = DiGraph([
            0 1 0 0;
            0 0 1 0;
            1 0 1 0;
            0 0 0 1;
        ])
        cycles = simplecycles_limited_length(selfloopg, nv(selfloopg))
        @test [3] in cycles
        @test [4] in cycles
        @test [1, 2, 3] in cycles || [2, 3, 1] in cycles || [3, 1, 2] in cycles
        @test length(cycles) == 3
    end

    @testset "octahedral graph" begin
        octag = smallgraph(:octahedral)
        octadg = DiGraph(octag)
        octalengths, _ = simplecycleslength(octadg)
        for k = 1:6
            @test sum(octalengths[1:k]) == length(simplecycles_limited_length(octag, k))
            @test sum(octalengths[1:k]) == length(simplecycles_limited_length(octadg, k))
        end
    end
end


text:
module PkgDev

using Compat

export Entry, Generate, GitHub

include("utils.jl")
include("github.jl")
include("entry.jl")
include("license.jl")
include("generate.jl")

const cd = Pkg.Dir.cd

"""
    dir(pkg, [paths...])

Returns package `pkg` directory location through search. Additional `paths` are appended.
"""
function dir(pkg::AbstractString)
    pkgsrc = Base.find_in_path(bytestring(pkg), Pkg.dir())
    pkgsrc === nothing && return ""
    abspath(dirname(pkgsrc), "..") |> realpath
end
dir(pkg::AbstractString, args...) = normpath(dir(pkg),args...)

"""
    register(pkg, [url])

Register `pkg` at the git URL `url`, defaulting to the configured origin URL of the git repo `Pkg.dir(pkg)`.
"""
register(pkg::AbstractString) = Entry.register(pkg)
register(pkg::AbstractString, url::AbstractString) = Entry.register(pkg,url)

"""
    tag(pkg, [ver, [commit]])

Tag `commit` as version `ver` of package `pkg` and create a version entry in `METADATA`. If not provided, `commit` defaults to the current commit of the `pkg` repo. If `ver` is one of the symbols `:patch`, `:minor`, `:major` the next patch, minor or major version is used. If `ver` is not provided, it defaults to `:patch`.
"""
tag(pkg::AbstractString, sym::Symbol=:patch) = cd(Entry.tag,pkg,sym)
tag(pkg::AbstractString, sym::Symbol, commit::AbstractString) = cd(Entry.tag,pkg,sym,false,commit)

tag(pkg::AbstractString, ver::VersionNumber; force::Bool=false) = cd(Entry.tag,pkg,ver,force)
tag(pkg::AbstractString, ver::VersionNumber, commit::AbstractString; force::Bool=false) =
    cd(Entry.tag,pkg,ver,force,commit)

submit(pkg::AbstractString) = cd(Entry.submit, pkg)
submit(pkg::AbstractString, commit::AbstractString) = cd(Entry.submit,pkg,commit)

"""
    publish()

For each new package version tagged in `METADATA` not already published, make sure that the tagged package commits have been pushed to the repo at the registered URL for the package and if they all have, open a pull request to `METADATA`.

A branch in user's remote fork of `METADATA` repository will be created that can be used to create pull request to `METADATA` package registry.
Optionally, function accepts a name for a pull request branch. If it isn't provided name will be automatically generated.
"""
publish(prbranch::AbstractString="") = Entry.publish(Pkg.Dir.getmetabranch(), prbranch)

doc"""
    generate(pkg,license)

Generate a new package named `pkg` with one of these license keys: `"MIT"`, `"BSD"`, `"ASL"` or `"MPL"`. If you want to make a package with a different license, you can edit it afterwards. Generate creates a git repo at `Pkg.dir(pkg)` for the package and inside it `LICENSE.md`, `README.md`, `REQUIRE`, and the julia entrypoint `$pkg/src/$pkg.jl`. Travis, AppVeyor CI configuration files `.travis.yml` and `appveyor.yml` with code coverage statistics using Coveralls or Codecov are created by default, but each can be disabled individually by setting `travis`, `appveyor` or `coverage` to `false`.
"""
generate(pkg::AbstractString, license::AbstractString;
         force::Bool=false, authors::Union{AbstractString,Array} = [],
         config::Dict=Dict(), path::AbstractString = Pkg.Dir.path(),
         travis::Bool = true, appveyor::Bool = true, coverage::Bool = true) =
    Generate.package(pkg, license, force=force, authors=authors, config=config, path=path,
                     travis=travis, appveyor=appveyor, coverage=coverage)

"""
    config()
Interactive configuration of the development environment.

PDK operations require `git` minimum configuration that keeps user signature (user.name & user.email).
"""
function config(force::Bool=false)
    # setup global git configuration
    cfg = LibGit2.GitConfig(LibGit2.Consts.CONFIG_LEVEL_GLOBAL)
    try
        println("Julia PDK configuration:")

        username = LibGit2.get(cfg, "user.name", "")
        if isempty(username) || force
            username = LibGit2.prompt("Enter user name", default=username)
            LibGit2.set!(cfg, "user.name", username)
        else
            println("User name: $username")
        end

        useremail = LibGit2.get(cfg, "user.email", "")
        if isempty(useremail) || force
            useremail = LibGit2.prompt("Enter user email", default=useremail)
            LibGit2.set!(cfg, "user.email", useremail)
        else
            println("User email: $useremail")
        end

        # setup github account
        ghuser = LibGit2.get(cfg, "github.user", "")
        if isempty(ghuser) || force
            ghuser = LibGit2.prompt("Enter GitHub user", default=(isempty(ghuser) ? username : ghuser))
            LibGit2.set!(cfg, "github.user", ghuser)
        else
            println("GitHub user: $ghuser")
        end
    finally
        finalize(cfg)
    end
    lowercase(LibGit2.prompt("Do you want to change this configuration?", default="N")) == "y" && config(true)
    return
end

function __init__()
    # Check if git configuration exists
    cfg = LibGit2.GitConfig(LibGit2.Consts.CONFIG_LEVEL_GLOBAL)
    try
        username = LibGit2.get(cfg, "user.name", "")
        if isempty(username)
            warn("Julia PDK is not configured. Please, run `PkgDev.config()` before performing any operations.")
        end
    finally
        finalize(cfg)
    end
end

end # module


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

mutable struct IoK8sApiCoreV1ExecAction <: SwaggerModel
    command::Any # spec type: Union{ Nothing, Vector{String} } # spec name: command

    function IoK8sApiCoreV1ExecAction(;command=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ExecAction, Symbol("command"), command)
        setfield!(o, Symbol("command"), command)
        o
    end
end # type IoK8sApiCoreV1ExecAction

const _property_map_IoK8sApiCoreV1ExecAction = Dict{Symbol,Symbol}(Symbol("command")=>Symbol("command"))
const _property_types_IoK8sApiCoreV1ExecAction = Dict{Symbol,String}(Symbol("command")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1ExecAction }) = collect(keys(_property_map_IoK8sApiCoreV1ExecAction))
Swagger.property_type(::Type{ IoK8sApiCoreV1ExecAction }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ExecAction[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ExecAction }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ExecAction[property_name]

function check_required(o::IoK8sApiCoreV1ExecAction)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ExecAction }, name::Symbol, val)
end


text:
fejer_plan1(μ) = FFTW.plan_r2r!(μ, FFTW.REDFT01)
fejer_plan2(μ) = FFTW.plan_r2r!(μ, FFTW.RODFT00)

"""
Compute nodes and weights of Fejer's first quadrature rule with a Jacobi weight.
"""
fejer1{T<:AbstractFloat}(N::Int,α::T,β::T) = fejer1(N,α,β,fejer_plan1(zeros(T,N)))

"""
Compute nodes and weights of Fejer's second quadrature rule with a Jacobi weight.
"""
fejer2{T<:AbstractFloat}(N::Int,α::T,β::T) = fejer2(N,α,β,fejer_plan2(zeros(T,N)))

fejer1{T<:AbstractFloat}(N::Int,α::T,β::T,plan) = T[sinpi((N-2k-one(T))/2N) for k=0:N-1],fejerweights1(N,α,β,plan)
fejer2{T<:AbstractFloat}(N::Int,α::T,β::T,plan) = T[cospi((k+one(T))/(N+one(T))) for k=0:N-1],fejerweights2(N,α,β,plan)

"""
Compute weights of Fejer's first quadrature rule with a Jacobi weight.
"""
fejerweights1{T<:AbstractFloat}(N::Int,α::T,β::T) = fejerweights1(N,α,β,fejer_plan1(zeros(T,N)))

"""
Compute weights of Fejer's second quadrature rule with a Jacobi weight.
"""
fejerweights2{T<:AbstractFloat}(N::Int,α::T,β::T) = fejerweights2(N,α,β,fejer_plan2(zeros(T,N)))

function fejerweights1{T<:AbstractFloat}(N::Int,α::T,β::T,plan)
    μ = chebyshevjacobimoments1(N,α,β)
    scale!(μ,inv(T(N)))
    return plan*μ
end

function fejerweights2{T<:AbstractFloat}(N::Int,α::T,β::T,plan)
    μ = chebyshevjacobimoments2(N,α,β)
    Np1 = N+one(T)
    scale!(μ,inv(Np1))
    plan*μ
    @inbounds for i=1:N μ[i] = sinpi(i/Np1)*μ[i] end
    return μ
end


text:
"""
    OlsSureFn(Y,X,NWQ=false,m=0)

LS of Y on X; where Y is Txn, and X is the same for all regressions

# Usage
(b,res,Yhat,Covb,R2) = OlsSureFn(Y,X,NWQ,m)

# Input
- `Y::Matrix`:     Txn, the n dependent variables
- `X::Matrix`:     Txk matrix of regressors (including deterministic ones)
- `NWQ:Bool`:      if true, then Newey-West's covariance matrix is used, otherwise Gauss-Markov
- `m::Int`:        scalar, bandwidth in Newey-West

# Output
- `b::Matrix`:     n*kx1, regression coefficients
- `u::Matrix`:     Txn, residuals Y - Yhat
- `Yhat::Matrix`:  Txn, fitted values X*b
- `V::Matrix`:     covariance matrix of vec(b)
- `R2::Vector`:    n-vector, R2 values

"""
function OlsSureFn(Y,X,NWQ=false,m=0)

    (T,n) = (size(Y,1),size(Y,2))
    k     = size(X,2)

    b     = X\Y
    Yhat  = X*b
    u     = Y - Yhat

    g     = zeros(T,n*k)
    for i = 1:n
      vv      = (1+(i-1)*k):(i*k)   #1:k,(1+k):2k,...
      g[:,vv] = X.*u[:,i]           #moment conditions for Y[:,i] regression
    end

    Sxx = X'X
    if NWQ
        S     = CovNWFn(g,m)            #Newey-West covariance matrix
        Sxx_1 = kron(I(n),inv(Sxx))
        V     = Sxx_1 * S * Sxx_1
    else
        V = kron(cov(u),inv(Sxx))      #traditional covariance matrix, Gauss-Markov
    end

    R2   = 1 .- var(u,dims=1)./var(Y,dims=1)

    return b, u, Yhat, V, R2

end
#------------------------------------------------------------------------------

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
"""
This file is used to run the full simulation pipeline, including running the generative
model, merging and saving files and creating a
"""

include("simulate_states.jl")

# parameters to pass to the main function
const file_date = ARGS[1]
const state = ARGS[2]

const run_simulation = true
if length(ARGS) > 2
    run_simulation = ARGS[3] == "false" ? false : true
end

# const nsims = parse(Int, ARGS[2])

const nsims = Dict{String, Int}(
    "NSW" => 1000,
    "QLD" => 1000,
    "SA" => 1000,
    "VIC" => 1000,
    "WA" => 1000,
    "ACT" => 1000,
    "NT" => 1000,
    "TAS" => 1000,
)

# run main
# run main simulation script
simulate_single_state(
    file_date,
    state,
    nsims[state],
    run_simulation,
    adjust_TP = true,
)


text:
using StatsBase
export PutBlock, put, Swap, swap, PSwap, pswap

"""
    PutBlock{D,C,GT<:AbstractBlock} <: AbstractContainer{GT,D}

Type for putting a block at given locations.
"""
struct PutBlock{D,C,GT<:AbstractBlock} <: AbstractContainer{GT,D}
    n::Int
    content::GT
    locs::NTuple{C,Int}

    function PutBlock(n::Int,block::GT, locs::NTuple{C,Int}) where {D,C,GT<:AbstractBlock{D}}
        @assert_locs_safe n locs
        @assert nqudits(block) == C "number of locations doesn't match the size of block"
        return new{D,C,GT}(n, block, locs)
    end
end
nqudits(pb::PutBlock) = pb.n

"""
    put(total::Int, pair)

Create a [`PutBlock`](@ref) with total number of active qubits, and a pair of
location and block to put on.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> put(4, 1=>X)
nqudits: 4
put on (1)
└─ X
```

If you want to put a multi-qubit gate on specific locations, you need to write down all possible locations.

```jldoctest; setup=:(using YaoBlocks)
julia> put(4, (1, 3)=>kron(X, Y))
nqudits: 4
put on (1, 3)
└─ kron
   ├─ 1=>X
   └─ 2=>Y
```

The outter locations creates a scope which make it seems to be a contiguous two qubits for the block inside `PutBlock`.

!!! tips
    It is better to use [`subroutine`](@ref) instead of `put` for large blocks, since put will use the matrix of its contents
    directly instead of making use of what's in it. `put` is more efficient for small blocks.
"""
put(total::Int, pa::Pair{NTuple{M,Int},<:AbstractBlock}) where {M} =
    PutBlock(total, pa.second, pa.first)
put(total::Int, pa::Pair{Int,<:AbstractBlock}) = PutBlock(total, pa.second, (pa.first,))
put(total::Int, pa::Pair{<:Any,<:AbstractBlock}) =
    PutBlock(total, pa.second, Tuple(pa.first))

"""
    put(pair) -> f(n)

Lazy curried version of [`put`](@ref).

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> put(1=>X)
(n -> put(n, 1 => X))
```
"""
put(pa::Pair) = @λ(n -> put(n, pa))

occupied_locs(x::PutBlock) = map(i -> x.locs[i], x.content |> occupied_locs)
chsubblocks(x::PutBlock, b::AbstractBlock) = PutBlock(x.n, b, x.locs)
PropertyTrait(::PutBlock) = PreserveAll()
cache_key(pb::PutBlock) = cache_key(pb.content)

mat(::Type{T}, pb::PutBlock{2,1}) where {T} = u1mat(pb.n, mat(T, pb.content), pb.locs...)
mat(::Type{T}, pb::PutBlock{2,C}) where {T,C} = unmat(pb.n, mat(T, pb.content), pb.locs)

function _apply!(r::AbstractRegister, pb::PutBlock{2})
    instruct!(r, mat_matchreg(r, pb.content), pb.locs)
    return r
end

# NOTE: Roger: these specialization should be removed after the new interpret
# mechanism is implemented
# specialization
for G in [:X, :Y, :Z, :T, :S, :Sdag, :Tdag, :H]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))
    @eval function _apply!(r::AbstractRegister, pb::PutBlock{2,C,<:$GT}) where {C}
        instruct!(r, Val($(QuoteNode(G))), pb.locs)
        return r
    end
end

Base.adjoint(x::PutBlock) = PutBlock(nqudits(x), adjoint(content(x)), x.locs)
Base.copy(x::PutBlock) = PutBlock(x.n, x.content, x.locs)
function Base.:(==)(lhs::PutBlock{D,C,GT}, rhs::PutBlock{D,C,GT}) where {D,C,GT}
    return (lhs.n == rhs.n) && (lhs.content == rhs.content) && (lhs.locs == rhs.locs)
end

function YaoBase.iscommute(x::PutBlock{D}, y::PutBlock{D}) where {D}
    _check_block_sizes(x, y)
    if x.locs == y.locs
        return iscommute(x.content, y.content)
    else
        return iscommute_fallback(x, y)
    end
end

const Swap = PutBlock{2,2,G} where {G<:ConstGate.SWAPGate}
const PSwap{T} = PutBlock{2,2,RotationGate{2,T,G}} where {G<:ConstGate.SWAPGate}
Swap(n::Int, locs::Tuple{Int,Int}) = PutBlock(n, ConstGate.SWAPGate(), locs)
PSwap(n::Int, locs::Tuple{Int,Int}, θ::Real) =
    PutBlock(n, rot(ConstGate.SWAPGate(), θ), locs)

"""
    swap(n, loc1, loc2)

Create a `n`-qubit [`Swap`](@ref) gate which swap `loc1` and `loc2`.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> swap(4, 1, 2)
nqudits: 4
put on (1, 2)
└─ SWAP
```
"""
swap(n::Int, loc1::Int, loc2::Int) = Swap(n, (loc1, loc2))

"""
    swap(loc1, loc2) -> f(n)

Create a lambda that takes the total number of active qubits as input. Lazy curried
version of `swap(n, loc1, loc2)`. See also [`Swap`](@ref).

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> swap(1, 2)
(n -> swap(n, 1, 2))
```
"""
swap(loc1::Int, loc2::Int) = @λ(n -> swap(n, loc1, loc2))

function mat(::Type{T}, g::Swap) where {T}
    mask = bmask(g.locs[1], g.locs[2])
    orders = map(b -> swapbits(b, mask) + 1, basis(g.n))
    return PermMatrix(orders, ones(T, nlevel(g)^g.n))
end

_apply!(r::AbstractRegister, g::Swap) = (instruct!(r, Val(:SWAP), g.locs); r)
occupied_locs(g::Swap) = g.locs

"""
    pswap(n::Int, i::Int, j::Int, α::Real)
    pswap(i::Int, j::Int, α::Real) -> f(n)

parametrized swap gate.
"""
pswap(n::Int, i::Int, j::Int, α::Real) = PSwap(n, (i, j), α)
pswap(i::Int, j::Int, α::Real) = n -> pswap(n, i, j, α)

for (G, GT) in [
    (:Rx, :(PutBlock{2,1,RotationGate{T,XGate}} where {T})),
    (:Ry, :(PutBlock{2,1,RotationGate{T,YGate}} where {T})),
    (:Rz, :(PutBlock{2,1,RotationGate{T,ZGate}} where {T})),
    (:PSWAP, :(PSwap)),
]
    @eval function _apply!(reg::AbstractRegister, g::$GT)
        instruct!(reg, Val($(QuoteNode(G))), g.locs, g.content.theta)
        return reg
    end
end


text:
using PosteriorSummary
using Test

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


text:
module Bravais

using Formatting,
RecipesBase

export BravaisLattice,
bravais_coordinates_unit_cell,
bravais_coordinates_unit_cell!,
bravais_coordinates

"""
    lattice = BravaisLattice(; a=Inf, α=0, b=Inf, β=π/2, x0=0, y0=0)

bravais lattice with:

- `a` is the length of the first vector (by default aligned along x-axis).

- `b` is length of second vector.

- `α` is the angle of first primitive vector.

- `β` is the angle of second primitive vector

- `x0`, `y0` define the origin of the lattice.
"""
struct BravaisLattice
    a::Float64
    b::Float64
    α::Float64
    β::Float64
    x0::Float64
    y0::Float64

    sin²θ::Float64
    cosθ::Float64
    R::Array{Float64,2}#SArray{Tuple{2,2},Float64,2,4}#
    Rᵀ::Array{Float64,2}#SArray{Tuple{2,2},Float64,2,4}#
    v1::Array{Float64,1}#SArray{Tuple{2},Float64,1,2}#
    v2::Array{Float64,1}#SArray{Tuple{2},Float64,1,2}#

    function BravaisLattice(;a::Number=Inf, α::Number=0, b::Number=Inf, β::Number=π/2, x0::Number=0, y0::Number=0)
        θ = β-α
        if iszero(θ)
            throw(ArgumentError("lattice angle θ=$(θ) cannot be zero"))
        end
        sinθ, cosθ = sincos(θ)
        # R = @SMatrix [ cos(α) -sin(α); sin(α)  cos(α)]
        R = [ cos(α) -sin(α); sin(α)  cos(α)]
        Rᵀ = transpose(R)
        # v1 = R*SVector(1,0)
        v1 = R*[1,0]
        # v2 = R*SVector(cosθ,sinθ)
        v2 = R*[cosθ,sinθ]
        new(float(a), float(b), float(α), float(β), float(x0), float(y0),
            sinθ^2, cosθ, R, Rᵀ, v1, v2)
    end

    function Base.show(io::IO, bvl::BravaisLattice)
        if !get(io, :sub, false)
            print(io, typeof(bvl), ": \n",
            "\tprimitive vector 1: ", fmt("3.2f",bvl.a), ", ∠", fmt("3.2f",(mod2pi(bvl.α))*180/pi), "°\n",
            "\tprimitive vector 2: ", fmt("3.2f",bvl.b), ", ∠", fmt("3.2f",(mod2pi(bvl.β))*180/pi), "°\n",
            "\torigin: (", fmt("3.2f",bvl.x0), ", ", fmt("3.2f",bvl.y0), ")")
        elseif !get(io, :sub1, false)
            print(io,
            "\n\tprimitive vector 1: ", fmt("3.2f",bvl.a), ", ∠", fmt("3.2f",(mod2pi(bvl.α))*180/pi), "°\n",
            "\tprimitive vector 2: ", fmt("3.2f",bvl.b), ", ∠", fmt("3.2f",(mod2pi(bvl.β))*180/pi), "°\n",
            "\torigin: (", fmt("3.2f",bvl.x0), ", ", fmt("3.2f",bvl.y0), ")")
        else
            print(io,
            "\n\t\tprimitive vector 1: ", fmt("3.2f",bvl.a), ", ∠", fmt("3.2f",(mod2pi(bvl.α))*180/pi), "°\n",
            "\t\tprimitive vector 2: ", fmt("3.2f",bvl.b), ", ∠", fmt("3.2f",(mod2pi(bvl.β))*180/pi), "°\n",
            "\t\torigin: (", fmt("3.2f",bvl.x0), ", ", fmt("3.2f",bvl.y0), ")")
        end
    end
end # end of struct BravaisLatice

"""
    lattice = BravaisLattice(lattice; :key1 => value1, :key2 => value2, ...)

new lattice from old lattice with modified fields.
"""
BravaisLattice(lattice::BravaisLattice; a=lattice.a, b=lattice.b, α=lattice.α, β=lattice.β, x0=lattice.x0, y0=lattice.y0) = BravaisLattice(;a=a, b=b, α=α, β=β, x0=x0, y0=y0)

################################################################################
# COORDINATE TRANSFORMATIONS
################################################################################
"""
    xb, yb = bravais_coordinates_unit_cell(x, y, lattice)

maps cartesian (`x`,`y`) into cartesian (`xb`,`yb`) unit cell specified by
lattice
"""
function bravais_coordinates_unit_cell(x::Real, y::Real, lattice::BravaisLattice)
    a, b, α, β, v1, v2, x0, y0 = lattice.a, lattice.b, lattice.α, lattice.β, lattice.v1, lattice.v2, lattice.x0, lattice.y0
    if (isinf(a) && iszero(β-π/2)) || (isinf(b) && iszero(α)) || (isinf(a) && isinf(b))
        return float(x), float(y)
    end

    p1, p2 = bravais_coordinates(x, y, lattice)
    p1 = mod(p1, a)
    p2 = mod(p2, b)

    xb = v1[1]*p1 + v2[1]*p2 + x0
    yb = v1[2]*p1 + v2[2]*p2 + y0

    return xb, yb
end

function bravais_coordinates_unit_cell(x::AbstractArray, y::AbstractArray, lattice::Union{Tlat,AbstractArray{Tlat}}) where Tlat<:BravaisLattice
    xb = Array{Float64}(undef, size(x)...)
    yb = Array{Float64}(undef, size(y)...)
    bravais_coordinates_unit_cell!(xb,yb,x,y,lattice)
    return xb, yb
end

function bravais_coordinates_unit_cell!(xb::AbstractArray, yb::AbstractArray, x::AbstractArray, y::AbstractArray, lattice::BravaisLattice)
    for k ∈ eachindex(x)
        xb[k], yb[k] = bravais_coordinates_unit_cell(x[k],y[k],lattice)
    end
    return nothing
end
function bravais_coordinates_unit_cell!(xb::AbstractArray, yb::AbstractArray, x::AbstractArray, y::AbstractArray, lattices::AbstractArray{Tlat}) where Tlat<:BravaisLattice
    for k ∈ eachindex(x)
        xb[k], yb[k] = bravais_coordinates_unit_cell(x[k],y[k],lattices[k])
    end
    return nothing
end

"""
    p1, p2 = bravais_coordinates(x,y,lattice)

coordinates in bravais frame (i.e. (x,y) = p1*v1 + p2*v2)
"""
function bravais_coordinates(x::Real, y::Real, lattice::BravaisLattice)
    a, b, α, β, v1, v2, x0, y0 = lattice.a, lattice.b, lattice.α, lattice.β, lattice.v1, lattice.v2, lattice.x0, lattice.y0
    if (isinf(a) && iszero(β-π/2)) || (isinf(b) && iszero(α)) || (isinf(a) && isinf(b))
        return float(x-x0), float(y-y0)
    end

    x += -x0
    y += -y0

    rv1 = v1[1]*x + v1[2]*y
    rv2 = v2[1]*x + v2[2]*y

    p1 = (rv1 - rv2*lattice.cosθ)/lattice.sin²θ
    p2 = (rv2 - rv1*lattice.cosθ)/lattice.sin²θ
    return p1, p2
end
function bravais_coordinates(x::AbstractArray,y::AbstractArray,lattice::BravaisLattice)
    p1 = Array{Float64,2}(undef,length(x),length(y))
    p2 = Array{Float64,2}(undef,length(x),length(y))
    bravais_coordinates!(p1,p2,x,y,lattice)
    return p1, p2
end
function bravais_coordinates!(p1::AbstractArray,p2::AbstractArray,x::AbstractArray,y::AbstractArray,lattice::BravaisLattice)
    if size(x)!==size(y)
        for i ∈ eachindex(x), j ∈ eachindex(y)
            p1[i,j],p2[i,j] = bravais_coordinates(x[i],y[j],lattice)
        end
    else
        for i ∈ eachindex(x)
            p1[i],p2[i] = bravais_coordinates(x[i],y[i],lattice)
        end
    end
    return nothing
end
bravais_coordinates(x::Real, y::AbstractArray, args...) = bravais_coordinates(hcat(x),y,args...)
bravais_coordinates(x::AbstractArray, y::Real, args...) = bravais_coordinates(x,vcat(y),args...)

########################################################################
### PLOTTING
########################################################################
"""
    plot(lattice::BravaisLattice, N=[4,4])

scatter plot of lattice vectors from -`N[i]`:`+N[i]` for each dimension i.

If `N` is a scalar, uses same span for both dimensions.
"""
@recipe function f(lattice::BravaisLattice, N::Array{Int,1}=[4,4])
    seriestype --> :scatter
    legend --> false
    aspect_ratio --> 1
    xlims --> [-N[1]*(lattice.a*lattice.v1[1]+lattice.b*lattice.v2[1]), +N[1]*(lattice.a*lattice.v1[1]+lattice.b*lattice.v2[1])]/2
    ylims --> [-N[2]*(lattice.a*lattice.v1[2]+lattice.b*lattice.v2[2]), +N[2]*(lattice.a*lattice.v1[2]+lattice.b*lattice.v2[2])]/2
    @series begin
        x, y = Float64[], Float64[]
        for n ∈ -N[1]:N[1], m ∈ -N[2]:N[2]
            temp = n*lattice.a*lattice.v1 + m*lattice.b*lattice.v2
            push!(x,temp[1]); push!(y,temp[2])
        end
        (x, y)
    end
    @series begin
        x, y = Float64[], Float64[]
        for i ∈ 1:5
            n = [0,1,1,0,0][i]
            m = [0,0,1,1,0][i]
            temp = n*lattice.a*lattice.v1 + m*lattice.b*lattice.v2
            push!(x,temp[1]); push!(y,temp[2])
        end
        seriestype --> :path
        fill --> (0,.15, :red)
        (x, y)
    end
end
@recipe f(lattice::BravaisLattice, N::Int) = (lattice, [N,N])

end # module


text:
# Practical Machine learning
# Decision Trees based learning - Random Forests example
# Chapter 5

# Criteria
abstract Criterion
abstract ClassificationCriterion <: Criterion
abstract RegressionCriterion <: Criterion

immutable GiniCriterion <: ClassificationCriterion; end  # Gini index
immutable CrossEntropyCriterion <: ClassificationCriterion; end  # cross entropy
immutable MSECriterion <: RegressionCriterion; end  # Mean Sequared Error
const Gini = GiniCriterion()
const CrossEntropy = CrossEntropyCriterion()
const MSE = MSECriterion()

# Splitters
abstract Splitter

type ClassificationSplitter{T<:ClassificationCriterion} <: Splitter
    samples::Vector{Int}
    feature::AbstractVector
    range::UnitRange{Int}
    example::Example

    n_left_samples::Float64
    n_right_samples::Float64
    left_label_counts::Vector{Float64}
    right_label_counts::Vector{Float64}

    function ClassificationSplitter(samples, feature, range, example)
        left_label_counts = zeros(example.n_labels)
        right_label_counts = zeros(example.n_labels)
        n_right_samples = 0.
        for i in 1:endof(range)
            s = samples[range[i]]
            w = example.sample_weight[s]
            l = example.y[s]
            n_right_samples += w
            right_label_counts[l] += w
        end
        new(samples, feature, range, example, 0., n_right_samples, left_label_counts, right_label_counts)
    end
end

type RegressionSplitter{T<:RegressionCriterion} <: Splitter
    samples::Vector{Int}
    feature::AbstractVector
    range::UnitRange{Int}
    example::Example

    n_left_samples::Float64
    n_right_samples::Float64
    sum_left::Float64
    sum_right::Float64
    sqsum_left::Float64
    sqsum_right::Float64

    function RegressionSplitter(samples, feature, range, example)
        n_right_samples = 0.
        sum_right = 0.
        sqsum_right = 0.
        for i in 1:endof(range)
            s = samples[range[i]]
            w = example.sample_weight[s]
            v = example.y[s]
            n_right_samples += w
            mass = v * w
            sum_right += mass
            sqsum_right += v * mass
        end
        new(samples, feature, range, example, 0., n_right_samples, 0., sum_right, 0., sqsum_right)
    end
end

immutable Split{T}
    threshold::T
    boundary::Int
    left_range::Range{Int}
    right_range::Range{Int}
    n_left_samples::Float64
    n_right_samples::Float64
    left_impurity::Float64
    right_impurity::Float64
end

Base.start(sp::Splitter) = 1

function Base.done(sp::Splitter, state::Int)
    # constant feature
    sp.feature[state] == sp.feature[end]
end

function Base.next(sp::ClassificationSplitter, state::Int)
    # seek for the next boundary
    local threshold, boundary = 0
    @inbounds for i in state:endof(sp.range)-1
        # transfer sample `s` from right to left
        s = sp.samples[sp.range[i]]
        w = sp.example.sample_weight[s]
        # as for classification the type of a label must be Int
        # this type annotation is useful for improving performance
        l = sp.example.y[s]::Int
        sp.n_left_samples += w
        sp.n_right_samples -= w
        sp.left_label_counts[l] += w
        sp.right_label_counts[l] -= w

        if sp.feature[i] != sp.feature[i+1]
            boundary = i
            threshold = sp.feature[i]
            break
        end
    end

    @assert boundary > 0

    left_range = sp.range[1:boundary]
    right_range = sp.range[boundary+1:end]
    left_impurity, right_impurity = impurity(sp)
    Split(threshold, boundary, left_range, right_range, sp.n_left_samples, sp.n_right_samples, left_impurity, right_impurity), boundary + 1
end

function Base.next(sp::RegressionSplitter, state::Int)
    # seek for the next boundary
    local threshold, boundary = 0
    @inbounds for i in state:endof(sp.range)-1
        # transfer sample `s` from right to left
        s = sp.samples[sp.range[i]]
        w = sp.example.sample_weight[s]
        v = sp.example.y[s]
        sp.n_left_samples += w
        sp.n_right_samples -= w
        mass = v * w
        sp.sum_left += mass
        sp.sum_right -= mass
        mass *= v  # mass = (v * v) * w
        sp.sqsum_left += mass
        sp.sqsum_right -= mass

        if sp.feature[i] != sp.feature[i+1]
            boundary = i
            threshold = sp.feature[i]
            break
        end
    end

    @assert boundary > 0

    left_range = sp.range[1:boundary]
    right_range = sp.range[boundary+1:end]
    left_impurity, right_impurity = impurity(sp)
    Split(threshold, boundary, left_range, right_range, sp.n_left_samples, sp.n_right_samples, left_impurity, right_impurity), boundary + 1
end

function count_labels(samples::Vector{Int}, example::Example)
    counts = zeros(Float64, example.n_labels)
    n_samples = 0.
    for s in samples
        n_samples += example.sample_weight[s]
        label = example.y[s]
        counts[label] += example.sample_weight[s]
    end

    counts, n_samples
end

function impurity(sp::ClassificationSplitter{GiniCriterion})
    left_impurity = 0.
    right_impurity = 0.

    for i in 1:endof(sp.left_label_counts)
        r = sp.left_label_counts[i] / sp.n_left_samples
        left_impurity += r * r

        r = sp.right_label_counts[i] / sp.n_right_samples
        right_impurity += r * r
    end

    1. - left_impurity, 1. - right_impurity
end

function impurity(samples::Vector{Int}, example::Example, ::GiniCriterion)
    counts, n_samples = count_labels(samples, example)
    gini_index = 0.
    for c in counts
        r = c / n_samples
        gini_index += r * r
    end
    1. - gini_index
end

function impurity(sp::ClassificationSplitter{CrossEntropyCriterion})
    left_impurity = 0.
    right_impurity = 0.

    for i in 1:endof(sp.left_label_counts)
        p = sp.left_label_counts[i] / sp.n_left_samples
        left_impurity -= p == 0. ? 0. : p * log(p)

        p = sp.right_label_counts[i] / sp.n_right_samples
        right_impurity -= p == 0. ? 0. : p * log(p)
    end

    left_impurity, right_impurity
end

function impurity(samples::Vector{Int}, example::Example, ::CrossEntropyCriterion)
    counts, n_samples = count_labels(samples, example)
    cross_entropy = 0.
    for c in counts
        p = c / n_samples
        if p != 0.
            cross_entropy -= p * log(p)
        end
    end
    cross_entropy
end

function impurity(sp::RegressionSplitter{MSECriterion})
    left_mean = sp.sum_left / sp.n_left_samples
    right_mean = sp.sum_right / sp.n_right_samples

    # sample variance
    (sp.sqsum_left / sp.n_left_samples - left_mean * left_mean,
     sp.sqsum_right / sp.n_right_samples - right_mean * right_mean)
end

function impurity(samples::Vector{Int}, example::Example, ::MSECriterion)
    n_samples = 0.
    sum = 0.
    sqsum = 0.

    for s in samples
        n_samples += example.sample_weight[s]
        v = example.y[s]
        sum += v
        sqsum += v * v
    end

    mean = sum / n_samples
    sqsum / n_samples - mean * mean
end


text:
@testset "Esat_from_Tair" begin
  Tair = 15
  Esat_formula=Sonntag1990()
  constants=BigleafConstants()
  eSat = Esat_from_Tair(Tair; Esat_formula, constants)
  # regression test
  @test isapprox(eSat, 1.702, atol=1e-3)
  eSat2 = Esat_from_Tair(Tair; Esat_formula = Alduchov1996(), constants)
  eSat3 = Esat_from_Tair(Tair; Esat_formula = Allen1998(), constants)
  # different results with different methods
  @test !(eSat2 ≈ eSat)
  @test !(eSat3 ≈ eSat)
  @test !(eSat3 ≈ eSat2)
  # but very similar
  @test ≈(eSat2, eSat, rtol = 1e-2)
  @test ≈(eSat3, eSat, rtol = 1e-2)
end

@testset "Esat_from_Tair_slope" begin
    step = 0.1
    Tair = range(20.0, 22.0, step = step)
    Esat_formula=Sonntag1990()
    constants=BigleafConstants()
    Esat = Esat_from_Tair.(Tair)
    delta = Esat_from_Tair_deriv.(Tair)
    delta2 = diff(Esat)/step
    @test all(isapprox.(delta[2:end] - delta2, 0, atol=1e-3))
    Esat3, delta3 = Esat_slope(Tair[1])
    @test Esat3 ≈ Esat[1]
    @test delta3 ≈ delta[1]
end

@testset "LE_to_ET" begin
  ET0 = 200.0
  Tair = 25.0
  res = LE_to_ET(ET0, Tair)
  # regression test
  @test ≈(res, 8.19e-5, atol =1e-7)
  ET = ET_to_LE(res, Tair)
  @test ET ≈ ET0
end

@testset "LE_to_ET" begin
  G_ms,Tair,pressure = 0.005,25.0,100.0
  rmol = ms_to_mol(G_ms,Tair,pressure)
  # regression test
  @test ≈(rmol, 0.2017, atol =1e-4)
  rms = mol_to_ms(rmol, Tair, pressure)
  @test rms ≈ G_ms
end

@testset "e_to_rH" begin
  Tair = 25.0
  eSat = Esat_from_Tair(Tair)
  e = range(0,eSat, length = 5)
  vpd = e_to_VPD.(e, Tair)
  @test VPD_to_e.(vpd, Tair) ≈ e
  rH0 = VPD_to_rH.(vpd,Tair)
  @test rH_to_VPD.(rH0, Tair) ≈ vpd
  rH = e_to_rH.(e, Tair)
  @test  rH ≈ rH0
  #
  rHOversat = @test_logs (:warn,r"higher than saturation") e_to_rH(eSat + 1e-3, Tair)
  @test rHOversat == 1.0
  #
  vpd2 = @test_logs (:warn,r"between 0 and 1") rH_to_VPD(1.2, Tair)
  @test vpd2 < 0 # garbage in garbage out
end

@testset "e_to_q" begin
  pressure = 0.1 #kPa
  Tair = 25.0
  eSat = Esat_from_Tair(Tair)
  e = range(0,eSat, length = 5)
  q = e_to_q.(e, pressure)
  @test q_to_e.(q, pressure) ≈ e
  vpd = q_to_VPD.(q, Tair, pressure)
  @test VPD_to_q.(vpd, Tair, pressure) ≈ q
end

@testset "Rg_to_PPFD" begin
  Rg = 500.0
  PPFD = Rg_to_PPFD(Rg)
  @test PPFD_to_Rg(PPFD) ≈ Rg
end

@testset "kg_to_mol" begin
  mass = 10.0
  molarMass =BigleafConstants().H2Omol
  amountOfSubstance = kg_to_mol(mass, molarMass)
  @test amountOfSubstance == mass/molarMass
end

@testset "umolCO2_to_gC" begin
  CO2_flux = 50.0
  C_flux = umolCO2_to_gC(CO2_flux)
  @test gC_to_umolCO2(C_flux) ≈ CO2_flux
end


text:
module TopologyFromGeometry

#using dataset_pb2

include("topology_from_geom.jl")

export MatchingParameters
export bond_topologies_from_geom

end  # module TopologyFromGeom

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_gbifautosuggest

"""
    ''_gbifautosuggest(;kwargs...)

A GbifAutosuggest component.

Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `label` (String; required): A label that will be printed when this component is rendered.
- `value` (Bool | Real | String | Dict | Array; optional): The value displayed in the input.
"""
function ''_gbifautosuggest(; kwargs...)
        available_props = Symbol[:id, :label, :value]
        wild_props = Symbol[]
        return Component("''_gbifautosuggest", "GbifAutosuggest", "gbif_autosuggest", available_props, wild_props; kwargs...)
end

max_stars_repo_name: nikitcha/ceebios-biowser

max_issues_repo_name: nikitcha/ceebios-biowser

max_forks_repo_name: nikitcha/ceebios-biowser


text:
"""
A mutable type that bundles multi-component records at receiver array.
```julia
records=Records(tgrid, ageom, [:p, :vx])

```
Here, we initialized records, of `:p` and `:vx` fields, in time domain for receivers and supersources in `ageom`.

## Indexing
* `records[i]` : the records due to ith supersource of all fields
* `records[i][:p]` : extract field `:p`
* `records[i][:nr]` : number of receivers (same as in `ageom`)
* `records.grid` : returns `tgrid`
As mutable objects in Julia are like containers that might hold different values over time, `records` can be modified.
"""
Records = Array{NamedD{Recs},1}

function Array{NamedD{Recs},1}(grid::StepRangeLen, ageom::AGeom, fields::Vector{Symbol})
    return [NamedD(grid, Recs(ageom[i].nr), fields) for i = 1:length(ageom)]
end


text:
# This file was generated, do not modify it. # hide
knnc.K = 3
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
@show mean(ŷ .== y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
import JUDI.Geometry, JUDI.get_nsrc, JUDI.n_samples, JUDI.subsample, JUDI.compareGeometry, JUDI.isequal

export H5GeometryOOC

mutable struct H5GeometryOOC{T} <: Geometry{T}
	h5geo::PyCall.PyObject
	container::H5SeisCon
  # JUDI GeometryOOC members start
  dt::Array{T,1}
  nt::Array{<:Integer,1}
  t::Array{T,1}
  nrec::Array{<:Integer,1}
  # JUDI GeometryOOC members end
	key::String
	xkey::String
	ykey::String
	zkey::String
  do_coord_transform::Bool
  model_origin_x::Number
  model_origin_y::Number
	model_orientation::Number
end

# mutable struct MyStruct{T}
#   a::Array{<:Integer,1}
#   b::Array{<:Integer,1}
#   c::Array{T,1}
# end

# a = ones(Int64, 2)
# b = ones(Integer, 3)
# c = ones(Float32, 3)

# MyStruct(a, b, c)

function H5GeometryOOC(;
  h5geo::PyCall.PyObject,
  container::H5SeisCon,
  key::String,
  xkey::String,
  ykey::String,
  zkey::String,
  do_coord_transform::Bool,
  model_origin_x::Number,
  model_origin_y::Number,
  model_orientation::Number)
  if isnothing(container)
    @error "H5SeisCon is Nothing\n"
    return
  end

  if isnothing(key)
    @error "key is Nothing\n"
    return
  end

  if key != "source" && key != "receiver"
    @error "key must be either 'source' or 'receiver'\n"
    return
  end

  if isnothing(xkey) || isnothing(ykey) || isnothing(zkey)
    @error "xkey/ykey/zkey is Nothing\n"
    return
  end

  nsrc = length(container)

  dt = zeros(Float32, nsrc) .+ Float32(abs(container.seis.getSampRate("ms")))
  nt = zeros(Integer, nsrc) .+ Integer(container.seis.getNSamp())
  t = zeros(Float32, nsrc) .+ Float32.((nt.-1).*dt)

  nrec = zeros(Integer, nsrc)
  i = 1
  for val in container.pkeyvals
    nrec[i] = Integer(container.seis.getPKeyTraceSize(container.pkey, val, val))
    i+=1
  end

  return H5GeometryOOC(h5geo, container, dt, nt, t, nrec,
                      key, xkey, ykey, zkey, do_coord_transform,
                      model_origin_x, model_origin_y, model_orientation)
end

######################## shapes easy access ################################
get_nsrc(g::H5GeometryOOC) = length(g.container.pkeyvals)

function n_samples(g::H5GeometryOOC, info::Info)
  return n_samples(g, info.nsrc)
end

n_samples(g::H5GeometryOOC, nsrc::Integer) = sum([g.nrec[j]*g.nt[j] for j=1:nsrc])

# # Subsample out-of-core geometry structure
# function subsample(geometry::H5GeometryOOC, srcnum::Number)
#   con = deepcopy(geometry.container)
#   con.pkeyvals = [con.pkeyvals[srcnum]]
#   return H5GeometryOOC(h5geo=geometry.h5geo, container=con, key=geometry.key,
#                       xkey=geometry.xkey, ykey=geometry.ykey, zkey=geometry.zkey,
#                       do_coord_transform=geometry.do_coord_transform,
#                       model_origin_x=geometry.model_origin_x,
#                       model_origin_y=geometry.model_origin_y,
#                       model_orientation=geometry.model_orientation)
# end

# # srcnum maybe vector
# function subsample(geometry::H5GeometryOOC, srcnum::AbstractArray)
#   con = deepcopy(geometry.container)
#   con.pkeyvals = con.pkeyvals[srcnum]
#   return H5GeometryOOC(h5geo=geometry.h5geo, container=con, key=geometry.key,
#                       xkey=geometry.xkey, ykey=geometry.ykey, zkey=geometry.zkey,
#                       do_coord_transform=geometry.do_coord_transform,
#                       model_origin_x=geometry.model_origin_x,
#                       model_origin_y=geometry.model_origin_y,
#                       model_orientation=geometry.model_orientation)
# end

# getindex out-of-core geometry structure
function getindex(geometry::H5GeometryOOC{T}, srcnum::JUDI.RangeOrVec) where T
  con = deepcopy(geometry.container)
  con.pkeyvals = con.pkeyvals[srcnum]
  return H5GeometryOOC(h5geo=geometry.h5geo, container=con, key=geometry.key,
                      xkey=geometry.xkey, ykey=geometry.ykey, zkey=geometry.zkey,
                      do_coord_transform=geometry.do_coord_transform,
                      model_origin_x=geometry.model_origin_x,
                      model_origin_y=geometry.model_origin_y,
                      model_orientation=geometry.model_orientation)
end

getindex(geometry::H5GeometryOOC{T}, srcnum::Integer) where T = getindex(geometry, srcnum:srcnum)

# Compare geometries
function compareGeometry(geometry_A::H5GeometryOOC, geometry_B::H5GeometryOOC)
  if geometry_A.key != geometry_B.key
    return false
  end

  if geometry_A.xkey != geometry_B.xkey ||
    geometry_A.ykey != geometry_B.ykey ||
    geometry_A.zkey != geometry_B.zkey
    return false
  end

  if !geometry_A.container.seis.isEqual(geometry_B.container.seis)
    return false
  end

  if geometry_A.container.pkey != geometry_B.container.pkey
    return false
  end

  if geometry_A.dt != geometry_B.dt
    return false
  end

  if geometry_A.nt != geometry_B.nt
    return false
  end

  if geometry_A.t != geometry_B.t
    return false
  end

  if geometry_A.nrec != geometry_B.nrec
    return false
  end

  return geometry_A.container.pkeyvals == geometry_B.container.pkeyvals
end

isequal(geometry_A::H5GeometryOOC, geometry_B::H5GeometryOOC) = compareGeometry(geometry_A, geometry_B)

compareGeometry(geometry_A::H5GeometryOOC, geometry_B::Geometry) = true
compareGeometry(geometry_A::Geometry, geometry_B::H5GeometryOOC) = true

######################## shapes easy access ################################
# Load geometry from out-of-core Geometry container
function Geometry(geometry::H5GeometryOOC)
  if isnothing(geometry.container.seis) || isnothing(geometry.container.seis) ||
    (geometry.key != "source" && geometry.key != "receiver")
    @error "Geometry key is Nothing or it is neither 'source' nor 'receiver'\n"
    return
  end

  if isnothing(geometry.container.seis)
    @error "Seis is Nothing\n"
    return
  end

  if isnothing(geometry.xkey) ||
    isnothing(geometry.ykey) ||
    isnothing(geometry.zkey)
    @error "xkey/ykey/zkey is Nothing\n"
    return
  end

  if isnothing(geometry.container.pkeyvals)
    @error "pkeyvals is Nothing\n"
    return
  end

  nsrc = length(geometry.container)
  xCell = Vector{Vector{Float32}}(undef, nsrc)
  yCell = Vector{Vector{Float32}}(undef, nsrc)
  zCell = Vector{Vector{Float32}}(undef, nsrc)
  dtCell = Vector{Float32}(undef, nsrc)
  ntCell = Vector{Integer}(undef, nsrc)
  tCell = Vector{Float32}(undef, nsrc)
  keylist = [geometry.container.pkey]
  for block in 1:nsrc
    minlist = [geometry.container.pkeyvals[block]]
    maxlist = [geometry.container.pkeyvals[block]]
    _, _, ind = geometry.container.seis.getSortedData(keylist, minlist, maxlist, 0, 0)
    if length(ind) < 1
      continue
    end
    if geometry.key == "source"
      ind = [ind[1]]
    end
    xy = Float32.(geometry.container.seis.getXYTraceHeaders(
      [geometry.xkey, geometry.ykey], ind, "m", geometry.do_coord_transform))

    xy[:,1] = xy[:,1] .- geometry.model_origin_x
    xy[:,2] = xy[:,2] .- geometry.model_origin_y
    xy = rotate_Nx2_array(xy, -geometry.model_orientation)
    xy[:,1] = xy[:,1] .+ geometry.model_origin_x
    xy[:,2] = xy[:,2] .+ geometry.model_origin_y

    xCell[block] = xy[:,1]
    yCell[block] = xy[:,2]
    zCell[block] = Float32.(geometry.container.seis.getTraceHeader(
      [geometry.zkey], ind, [geometry.container.seis.getLengthUnits()], ["m"])[:])
    dtCell[block] = geometry.dt[block]
    ntCell[block] = geometry.nt[block]
    tCell[block] = geometry.t[block]
  end

  return GeometryIC{Float32}(xCell, yCell, zCell, dtCell, ntCell, tCell)
end

function Geometry(;
  container::H5SeisCon,
  xkey::String,
  ykey::String,
  zkey::String,
  do_coord_transform::Bool,
  model_origin_x::Number,
  model_origin_y::Number,
  model_orientation::Number)

  if isnothing(container.seis)
    @error "Seis is Nothing\n"
    return
  end

  if isnothing(container.pkey) || isnothing(xkey) || isnothing(ykey) || isnothing(zkey)
    @error "pkey/xkey/ykey/zkey is Nothing\n"
    return
  end

  if isnothing(container.pkeyvals)
    @error "pkeyvals is Nothing\n"
    return
  end

  dt = Float32(abs(container.seis.getSampRate("ms")))
  nt = Integer(container.seis.getNSamp())
  t = Float32((nt-1)*dt)

  nsrc = length(container.pkeyvals)
  xCell = Vector{Vector{Float32}}(undef, nsrc)
  yCell = Vector{Vector{Float32}}(undef, nsrc)
  zCell = Vector{Vector{Float32}}(undef, nsrc)
  dtCell = Vector{Float32}(undef, nsrc)
  ntCell = Vector{Integer}(undef, nsrc)
  tCell = Vector{Float32}(undef, nsrc)
  keylist = [container.pkey]
  for block in 1:nsrc
    minlist = [container.pkeyvals[block]]
    maxlist = [container.pkeyvals[block]]
    _, _, ind = container.seis.getSortedData(keylist, minlist, maxlist, 0, 0)
    xy = Float32.(container.seis.getXYTraceHeaders([xkey, ykey], ind, "m", do_coord_transform))

    xy[:,1] = xy[:,1] .- model_origin_x
    xy[:,2] = xy[:,2] .- model_origin_y
    xy = rotate_Nx2_array(xy, -model_orientation)
    xy[:,1] = xy[:,1] .+ model_origin_x
    xy[:,2] = xy[:,2] .+ model_origin_y

    xCell[block] = xy[:,1]
    yCell[block] = xy[:,2]
    zCell[block] = Float32.(container.seis.getTraceHeader([zkey], ind, [container.seis.getLengthUnits()], ["m"])[:])
    dtCell[block] = dt
    ntCell[block] = nt
    tCell[block] = t
  end

  return GeometryIC{Float32}(xCell, yCell, zCell, dtCell, ntCell, tCell)
end


text:
include_string("""
    x = sum([1, 2, 3])
    @show x
    """)

@show typeof(x) # Int64

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
__precompile__(true)

module Transformations

# using Reexport
using RecipesBase

import CatViews: CatView, splitview
import Base: rand
using LearnBase
import LearnBase: transform, transform!,
                  grad, grad!, addgrad!,
                  value, params, update!,
                  learn!
# import StatsBase: logistic, logit
import OnlineStats
import OnlineStats: weight, Weight, BoundedEqualWeight

export
    input_node,
    output_node,
    input_length,
    output_length,
    input_value,
    output_value,
    input_grad,
    output_grad,
    params_length,
    consolidate_params,
    # params,
    totalcost,

    # Learnable,
    # Node,
    link_nodes!,
    Linear,
    Affine,
    LayerNorm,
    InputNorm,
    Activation,
    Concat,
    Gate,
    PreprocessStep,
    NoPreprocessing,
    Whiten,
    Chain,
    nnet,
    ResidualLayer,
    resnet,
    ConvFilter,
    ConvLayer,
    MaxPooling,
    Differentiable,
    tfunc,  # shorthand constructor for Differentiable
    NonDifferentiable,
    OnceDifferentiable,
    TwiceDifferentiable,
    MvNormalTransformation
    # Params

function zero!{T,N}(v::AbstractArray{T,N})
    for i in eachindex(v)
        v[i] = zero(T)
    end
end

input_node(t::Transformation) = t.input
output_node(t::Transformation) = t.output
input_length(t::Transformation) = length(input_node(t))
output_length(t::Transformation) = length(output_node(t))
input_value(t::Transformation) = value(input_node(t))
output_value(t::Transformation) = value(output_node(t))
input_grad(t::Transformation) = grad(input_node(t))
output_grad(t::Transformation) = grad(output_node(t))

# abstract Learnable <: Transformation

"The length of the parameter vector"
params_length(t::Transformation) = 0
params_length(t::Learnable) = length(params(t))

"Learnable parameter values"
params(t::Transformation) = zeros(0)
params(t::Learnable) = value(t.params)

"Gradients of the learnable parameters"
grad(t::Transformation) = zeros(0)
grad(t::Learnable) = grad(t.params)

# notes:
#   Transformations will be updated in a forward (transform) and backward (grad) pass.
#   The input nodes of a larger comp graph will be called with `transform!(t,input)` and all other
#   nodes will be called with `transform!(t)`, assuming they are properly "linked" beforehand.
#   This is because `output` is computed, which shares a reference to the arrays in the following
#   Transformation's `input`, so it's ready to compute with `transform!`.

#   The same happens in reverse.  An `input` node's `∇` is the same array as the child node's `output.∇`,
#   so the gradients can flow backwards with one call to `grad!` in the proper (reverse) order.
#   In this case, the

# TODO:
#   - DAGs of Transformations
#   - Then handle cycles

# ----------------------------------------------------------------

# Most Transformations can share these methods as long as they have the required fields.

# Copy input values into the input node, then transform
function transform!(t::Transformation, x::AbstractArray)
    # copy!(input_value(t), input)
    transform!(t.input, x)
    transform!(t)
end

# # # allow a transformation to be "called" to transform:  output = t(input)
# function (t::Transformation)(input::AbstractArray)
#     transform!(t, input)
# end
# function (t::Transformation)()
#     transform!(t)
# end

# Copy the gradient into the output node, and propagate it back.
function grad!(t::Transformation, ∇out::AbstractArray)
    # copy!(output_grad(t), ∇out)
    grad!(t.output, ∇out)
    grad!(t)
end

# do a forward and backward pass for an observation
function update!{T,S}(t::Minimizable, obs::Tuple{T,S})
    transform!(t, obs[2], obs[1])
    grad!(t)
end

# when obs == nothing, don't copy input data in transform!
function update!(t::Learnable, ::Void)
    transform!(t)
    grad!(t)
end

# ----------------------------------------------------------------

function initialize_weights!{T}(w::AbstractArray{T})
    nout, nin = size(w)
    # scalar = sqrt(T(2.0 / (nin + nout)))
    scalar = sqrt(T(1.0 / nin))
    for i in eachindex(w)
        w[i] = scalar * randn(T)
    end
end

function initialize_bias!{T}(b::AbstractArray{T})
    fill!(b, zero(T))
end

function smooth!{T,N}(a::AbstractArray{T,N}, b::AbstractArray{T,N}, λ::T)
    @assert length(a) == length(b)
    oml = one(λ) - λ
    @inbounds for i=1:length(a)
        a[i] = λ * b[i] + oml * a[i]
    end
end

# ----------------------------------------------------------------

include("params.jl")
include("nodes.jl")
include("linear.jl")
include("affine.jl")
include("misc.jl")
include("layernorm.jl")
# include("inputnorm.jl")
include("activations.jl")
include("whiten.jl")
include("chain.jl")
include("nnet.jl")
include("convolutions.jl")
include("pooling.jl")
include("functions.jl")
include("distributions.jl")

# ----------------------------------------------------------------

module Testing

using LearnBase
using ..Transformations

function check_gradient(t::Transformations.Transformation,
                        x = randn(input_length(t));
                        ϵ::Number = 1e-5)
    # first get: y = f(x)
    y = transform!(t, x)

    #= note: we assume the loss is:
        L = sum(y)
        ∂L/∂yᵢ = 1
    =#
    L = sum(y)
    grad!(t, ones(output_length(t)))

    # @show L x y t

    # now add a slight perterbation and check that it approximately
    # matches the analytic gradient

    if typeof(t) <: Learnable
        # compute the gradient
        Θ = copy(params(t))
        ∇ = grad(t)
        # @show Θ ∇

        Θ̃ = params(t)
        perr = zeros(length(Θ))
        for i=1:length(Θ)
            Θ̃[i] = Θ[i] + ϵ
            y = transform!(t, x)
            L1 = sum(y)
            Θ̃[i] = Θ[i] - ϵ
            y = transform!(t, x)
            L2 = sum(y)

            true_grad = (L1 - L2) / (2ϵ)
            perr[i] = true_grad - ∇[i]
            # @show i, L1, L2, true_grad, ∇[i], perr[i]
            if abs(true_grad) > 1e-8
                perr[i] /= true_grad
            end
            # @show perr[i]
        end
    else
        perr = zeros(0)
    end

    x = copy(x)
    x̃ = copy(x)
    ∇x = input_grad(t)
    xerr = zeros(length(x))
    for i=1:length(x)
        x̃[i] = x[i] + ϵ
        y = transform!(t, x̃)
        L1 = sum(y)
        x̃[i] = x[i] - ϵ
        y = transform!(t, x̃)
        L2 = sum(y)

        true_grad = (L1 - L2) / (2ϵ)
        xerr[i] = true_grad - ∇x[i]
        # @show L L1 L2 true_grad ∇x[i] xerr[i]
        if abs(true_grad) > 1e-8
            xerr[i] /= true_grad
        end
        # @show xerr[i]
    end
    perr, xerr
end

end # Testing

# ----------------------------------------------------------------

end # module


text:
struct MyAwesomeClass
    a::Int;
    b::Int;

    MyAwesomeClass(a::Int,b::Int) = new(a,b);

    function sum()::Int
        return (a + b);
    end;

    function subs()::Int
        return (a - b);
    end;
end;

myAwesomeClass = MyAwesomeClass(1,2);
println(myAwesomeClass.a," ",myAwesomeClass.b);
println(myAwesomeClass.subs());
println(myAwesomeClass.sum());

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
# Bogumił Kamiński, 2021

# Codes for appendix A

# As opposed to all other files in this repository in the presented codes I
# do not assume that you have activated the project from the repository folder.
# Instead, please follow the instructions contained in the appendix A.

# starting Julia session

julia

# exiting Julia session

exit()

# getting help

?&&

# setting up the project environment for working with examples in the book
# D:\\JuliaForDataAnalysis is an example folder

cd("D:\\JuliaForDataAnalysis")
isfile("Project.toml")
isfile("Manifest.toml")
# switch to package manager mode by pressing ]
activate .
instantiate
# exit package manager mode by pressing backspace

# activating project environment in a given folder
activate D:\\JuliaForDataAnalysis\\

# setting up a new project folder
# make sure your working directory is an empty folder

pwd()
# press ]
activate .
status
add BenchmarkTools
status
update
remove BenchmarkTools
status

# setting up integration with Python

# make sure you are in package manager mode
add PyCall
# press backspace
using PyCall
PyCall.python

# setting up integration with R

# make sure you are in package manager mode
add RCall
# press backspace
using RCall
RCall.Rhome

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis