text:
module JacobianTest

import Calculus

using Test
using ForwardDiff
using ForwardDiff: Dual, Tag, JacobianConfig
using StaticArrays
using DiffTests

include(joinpath(dirname(@__FILE__), "utils.jl"))

##################
# hardcoded test #
##################

f! = (y, x) -> begin
    y[1] = x[1] * x[2]
    y[1] *= sin(x[3]^2)
    y[2] = y[1] + x[3]
    y[3] = y[1] / y[2]
    y[4] = x[3]
    return nothing
end
f = x -> (y = fill(zero(promote_type(eltype(x), Float64)), 4); f!(y, x); return y)
x = [1, 2, 3]
v = f(x)
j = [0.8242369704835132  0.4121184852417566  -10.933563142616123
     0.8242369704835132  0.4121184852417566  -9.933563142616123
     0.169076696546684   0.084538348273342   -2.299173530851733
     0.0                 0.0                 1.0]

for c in (1, 2, 3), tags in ((nothing, nothing),
                             (Tag(f, eltype(x)), Tag(f!, eltype(x))))
    println("  ...running hardcoded test with chunk size = $c and tag = $(repr(tags))")
    cfg = JacobianConfig(f, x, ForwardDiff.Chunk{c}(), tags[1])
    ycfg = JacobianConfig(f!, fill(0.0, 4), x, ForwardDiff.Chunk{c}(), tags[2])

    @test eltype(cfg)  == Dual{typeof(tags[1]), eltype(x), c}
    @test eltype(ycfg) == Dual{typeof(tags[2]), eltype(x), c}

    # testing f(x)
    @test isapprox(j, ForwardDiff.jacobian(f, x, cfg))
    @test isapprox(j, ForwardDiff.jacobian(f, x))

    out = fill(0.0, 4, 3)
    ForwardDiff.jacobian!(out, f, x, cfg)
    @test isapprox(out, j)

    out = fill(0.0, 4, 3)
    ForwardDiff.jacobian!(out, f, x)
    @test isapprox(out, j)

    out = DiffResults.JacobianResult(fill(0.0, 4), fill(0.0, 3))
    ForwardDiff.jacobian!(out, f, x, cfg)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.jacobian(out), j)

    # testing f!(y, x)
    y = fill(0.0, 4)
    @test isapprox(j, ForwardDiff.jacobian(f!, y, x, ycfg))
    @test isapprox(v, y)

    y = fill(0.0, 4)
    @test isapprox(j, ForwardDiff.jacobian(f!, y, x))
    @test isapprox(v, y)

    out, y = fill(0.0, 4, 3), fill(0.0, 4)
    ForwardDiff.jacobian!(out, f!, y, x, ycfg)
    @test isapprox(out, j)
    @test isapprox(y, v)

    out, y = fill(0.0, 4, 3), fill(0.0, 4)
    ForwardDiff.jacobian!(out, f!, y, x)
    @test isapprox(out, j)
    @test isapprox(y, v)

    out = DiffResults.JacobianResult(fill(0.0, 4), fill(0.0, 3))
    y = fill(0.0, 4)
    ForwardDiff.jacobian!(out, f!, y, x, ycfg)
    @test DiffResults.value(out) == y
    @test isapprox(y, v)
    @test isapprox(DiffResults.jacobian(out), j)

    out = DiffResults.JacobianResult(fill(0.0, 4), fill(0.0, 3))
    y = fill(0.0, 4)
    ForwardDiff.jacobian!(out, f!, y, x)
    @test DiffResults.value(out) == y
    @test isapprox(y, v)
    @test isapprox(DiffResults.jacobian(out), j)
end

cfgx = ForwardDiff.JacobianConfig(sin, x)
@test_throws ForwardDiff.InvalidTagException ForwardDiff.jacobian(f, x, cfgx)
@test ForwardDiff.jacobian(f, x, cfgx, Val{false}()) == ForwardDiff.jacobian(f,x)

########################
# test vs. Calculus.jl #
########################

for f in DiffTests.ARRAY_TO_ARRAY_FUNCS
    v = f(X)
    j = ForwardDiff.jacobian(f, X)
    @test isapprox(j, Calculus.jacobian(x -> vec(f(x)), X, :forward), atol=1.3FINITEDIFF_ERROR)
    for c in CHUNK_SIZES, tag in (nothing, Tag)
        if tag == Tag
            tag = Tag(f, eltype(X))
        end
        println("  ...testing $f with chunk size = $c and tag = $(repr(tag))")
        cfg = JacobianConfig(f, X, ForwardDiff.Chunk{c}(), tag)

        out = ForwardDiff.jacobian(f, X, cfg)
        @test isapprox(out, j)

        out = similar(X, length(v), length(X))
        ForwardDiff.jacobian!(out, f, X, cfg)
        @test isapprox(out, j)

        out = DiffResults.DiffResult(similar(v, length(v)), similar(v, length(v), length(X)))
        ForwardDiff.jacobian!(out, f, X, cfg)
        @test isapprox(DiffResults.value(out), v)
        @test isapprox(DiffResults.jacobian(out), j)
    end
end

for f! in DiffTests.INPLACE_ARRAY_TO_ARRAY_FUNCS
    v = fill!(similar(Y), 0.0)
    f!(v, X)
    j = ForwardDiff.jacobian(f!, fill!(similar(Y), 0.0), X)
    @test isapprox(j, Calculus.jacobian(x -> (y = fill!(similar(Y), 0.0); f!(y, x); vec(y)), X, :forward), atol=FINITEDIFF_ERROR)
    for c in CHUNK_SIZES, tag in (nothing, Tag(f!, eltype(X)))
        println("  ...testing $(f!) with chunk size = $c and tag = $(repr(tag))")
        ycfg = JacobianConfig(f!, fill!(similar(Y), 0.0), X, ForwardDiff.Chunk{c}(), tag)

        y = fill!(similar(Y), 0.0)
        out = ForwardDiff.jacobian(f!, y, X, ycfg)
        @test isapprox(y, v)
        @test isapprox(out, j)

        y = fill!(similar(Y), 0.0)
        out = similar(Y, length(Y), length(X))
        ForwardDiff.jacobian!(out, f!, y, X)
        @test isapprox(y, v)
        @test isapprox(out, j)

        y = fill!(similar(Y), 0.0)
        out = DiffResults.JacobianResult(y, X)
        ForwardDiff.jacobian!(out, f!, y, X)
        @test DiffResults.value(out) == y
        @test isapprox(y, v)
        @test isapprox(DiffResults.jacobian(out), j)

        y = fill!(similar(Y), 0.0)
        out = DiffResults.JacobianResult(y, X)
        ForwardDiff.jacobian!(out, f!, y, X, ycfg)
        @test DiffResults.value(out) == y
        @test isapprox(y, v)
        @test isapprox(DiffResults.jacobian(out), j)
    end
end

##########################################
# test specialized StaticArray codepaths #
##########################################

println("  ...testing specialized StaticArray codepaths")

x = rand(3, 3)
for T in (StaticArrays.SArray, StaticArrays.MArray)
    sx = T{Tuple{3,3}}(x)

    cfg = ForwardDiff.JacobianConfig(nothing, x)
    scfg = ForwardDiff.JacobianConfig(nothing, sx)

    _diff(A) = diff(A; dims=1)

    actual = ForwardDiff.jacobian(_diff, x)
    @test ForwardDiff.jacobian(_diff, sx) == actual
    @test ForwardDiff.jacobian(_diff, sx, cfg) == actual
    @test ForwardDiff.jacobian(_diff, sx, scfg) == actual
    @test ForwardDiff.jacobian(_diff, sx, scfg) isa StaticArray
    @test ForwardDiff.jacobian(_diff, sx, scfg, Val{false}()) == actual
    @test ForwardDiff.jacobian(_diff, sx, scfg, Val{false}()) isa StaticArray

    out = similar(x, 6, 9)
    ForwardDiff.jacobian!(out, _diff, sx)
    @test out == actual

    out = similar(x, 6, 9)
    ForwardDiff.jacobian!(out, _diff, sx, cfg)
    @test out == actual

    out = similar(x, 6, 9)
    ForwardDiff.jacobian!(out, _diff, sx, scfg)
    @test out == actual

    result = DiffResults.JacobianResult(similar(x, 6), x)
    result = ForwardDiff.jacobian!(result, _diff, x)

    result1 = DiffResults.JacobianResult(similar(sx, 6), sx)
    result2 = DiffResults.JacobianResult(similar(sx, 6), sx)
    result3 = DiffResults.JacobianResult(similar(sx, 6), sx)
    result1 = ForwardDiff.jacobian!(result1, _diff, sx)
    result2 = ForwardDiff.jacobian!(result2, _diff, sx, cfg)
    result3 = ForwardDiff.jacobian!(result3, _diff, sx, scfg)
    @test DiffResults.value(result1) == DiffResults.value(result)
    @test DiffResults.value(result2) == DiffResults.value(result)
    @test DiffResults.value(result3) == DiffResults.value(result)
    @test DiffResults.jacobian(result1) == DiffResults.jacobian(result)
    @test DiffResults.jacobian(result2) == DiffResults.jacobian(result)
    @test DiffResults.jacobian(result3) == DiffResults.jacobian(result)

    sy = @SVector fill(zero(eltype(sx)), 6)
    sresult1 = DiffResults.JacobianResult(sy, sx)
    sresult2 = DiffResults.JacobianResult(sy, sx)
    sresult3 = DiffResults.JacobianResult(sy, sx)
    sresult1 = ForwardDiff.jacobian!(sresult1, _diff, sx)
    sresult2 = ForwardDiff.jacobian!(sresult2, _diff, sx, cfg)
    sresult3 = ForwardDiff.jacobian!(sresult3, _diff, sx, scfg)
    @test DiffResults.value(sresult1) == DiffResults.value(result)
    @test DiffResults.value(sresult2) == DiffResults.value(result)
    @test DiffResults.value(sresult3) == DiffResults.value(result)
    @test DiffResults.jacobian(sresult1) == DiffResults.jacobian(result)
    @test DiffResults.jacobian(sresult2) == DiffResults.jacobian(result)
    @test DiffResults.jacobian(sresult3) == DiffResults.jacobian(result)
end

@testset "dimension errors for jacobian" begin
    @test_throws DimensionMismatch ForwardDiff.jacobian(identity, 2pi) # input
    @test_throws DimensionMismatch ForwardDiff.jacobian(sum, fill(2pi, 2)) # vector_mode_jacobian
    @test_throws DimensionMismatch ForwardDiff.jacobian(sum, fill(2pi, 10^6)) # chunk_mode_jacobian
end

end # module


text:
@testset "Global" begin

function kernel(X)
    ptr = AMDGPU.get_global_pointer(Val(:myglobal), Float32)
    Base.unsafe_store!(ptr, 3f0)
    nothing
end

hk = AMDGPU.rocfunction(kernel, Tuple{Int32})
exe = hk.mod.exe
gbl = AMDGPU.get_global(exe.exe, :myglobal)
gbl_ptr = Base.unsafe_convert(Ptr{Float32}, gbl.ptr)
@test Base.unsafe_load(gbl_ptr) == 0f0
Base.unsafe_store!(gbl_ptr, 2f0)
@test Base.unsafe_load(gbl_ptr) == 2f0

dev = AMDGPU.default_device()
queue = AMDGPU.default_queue(dev)
kern = AMDGPU.create_kernel(dev, exe, hk.fun.entry, (Int32(3),))
signal = AMDGPU.create_event(hk.mod.exe)
AMDGPU.launch_kernel(queue, kern, signal; groupsize=(1,1,1), gridsize=(1,1,1))
wait(signal)
@test Base.unsafe_load(gbl_ptr) == 3f0

end

ext: jl


text:
# B-spline manifold

abstract type AbstractBSplineManifold end

function ⊗(X::AbstractArray{<:Real}, Y::AbstractArray{<:Real})
    # TODO: remove this function and use Tensor.jl
    m = size(X)
    n = size(Y)
    reshape(reshape(X, length(X)) * reshape(Y, length(Y))', m..., n...)
end

function tensorprod(X::AbstractArray)
    n = length(X)
    # X[1] ⊗ … ⊗ X[n]
    Y = X[1]
    for i in 2:n
        Y = Y ⊗ X[i]
    end
    return Y
end

"""
B-spline manifold for general polynomial degree
"""
struct BSplineManifold{T} <: AbstractBSplineManifold
    bsplinespaces::Vector{BSplineSpace}
    controlpoints::Array{T} where T<:Point
    function BSplineManifold(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{T}) where T
        Ps = BSplineSpace.(Ps)
        if collect(size(a)) ≠ dim.(Ps)
            throw(DimensionMismatch())
        else
            P = convert(Vector{BSplineSpace}, Ps)
            a′ = float(a)
            new{T}(P, a′)
        end
    end
end

"""
convert AbstractBSplineManifold to BSplineManifold
"""
function BSplineManifold(M::AbstractBSplineManifold)
    BSplineManifold(BSplineSpace.(M.bsplinespaces), M.controlpoints)
end

@doc raw"""
Multi-dimensional B-spline basis function.
```math
B_{i^1,\dots,i^d}(t^1,\dots,t^d)
=B_{(i^1,p^1,k^1)}(t^1)\cdots B_{(i^d,p^d,k^d)}(t^d)
```
"""
function bsplinebasis(Ps::AbstractVector{BSplineSpace}, t::AbstractVector{<:Real})
    d = length(t)
    Bs = [bsplinebasis(Ps[i], t[i]) for i in 1:d]
    return tensorprod(Bs)
end

@doc raw"""
Multi-dimensional B-spline basis function.
```math
B_{i^1,\dots,i^d}(t^1,\dots,t^d)
=B_{(i^1,p^1,k^1)}(t^1)\cdots B_{(i^d,p^d,k^d)}(t^d)
```
"""
function bsplinebasis(I::CartesianIndex, Ps::AbstractVector{BSplineSpace}, t::AbstractVector{<:Real})
    @warn "The method `bsplinebasis(I::CartesianIndex, Ps::AbstractVector{BSplineSpace}, t::AbstractVector{<:Real})` will be deprecated."
    d = length(Ps)
    Bs = prod(bsplinebasis(I[i], Ps[i], t[i]) for i in 1:d)
    return tensorprod(Bs)
end

function bsplinesupport(I::CartesianIndex, Ps::AbstractVector{BSplineSpace})
    @warn "The method `bsplinesupport(I::CartesianIndex, Ps::AbstractVector{BSplineSpace})` will be deprecated."
    d = length(Ps)
    return [bsplinesupport(I[i], Ps[i]) for i in 1:d]
end

@doc raw"""
Calculate the mapping of B-spline manifold for given parameter.
```math
\bm{p}(t^1,\dots,t^d)
=\sum_{i^1,\dots,i^d}B_{i^1,\dots,i^d}(t^1,\dots,t^d) \bm{a}_{i^1,\dots,i^d}
```
"""
function (M::BSplineManifold)(t::AbstractVector{<:Real})
    Ps = M.bsplinespaces
    a = M.controlpoints
    d = length(Ps)
    N = prod(dim.(Ps))

    B = bsplinebasis(Ps, t)
    B_flat = reshape(B,N)
    a_flat = reshape(a,N)
    return sum(B_flat .* a_flat[:])
end

@doc raw"""
Calculate the dimension of B-spline manifold.
"""
dim

dim(M::AbstractBSplineManifold) = length(M.bsplinespaces)

function bsplinespaces(M::BSplineManifold)
    return M.bsplinespaces
end

function controlpoints(M::BSplineManifold)
    return M.controlpoints
end


text:
module Tethys

using StaticArrays # Gotta go quick
using Printf

println("Loading Tethys...")

greet() = println("Hello Polarons! ♓")

include("BareExpansion.jl")

greet()

end # module


text:
function constantfunc(x)::Int
    return 1
end

function type2weightfunc(x::T, tol::T = eps(T))::T where T
    return cos(acos(x)/2)
end

function type3weightfunc(x::T, tol::T = eps(T))::T where T
    return sin(acos(x)) + 0.0001
end

function type4weightfunc(x::T, tol::T = eps(T))::T where T
    return sin(acos(x)/2) + 0.0001
end

function getweightfunc(type_number::Int)::Function
    @assert 1 <= type_number <= 4
    if type_number == 1
        return constantfunc
    elseif type_number == 2
        return type2weightfunc
    elseif type_number == 3
        return type3weightfunc
    end

    return type4weightfunc
end

# uses the Barycentric interpolation, Equation 3.17 of Filip's thesis.
function levellederror(wfunc::Function, w::Vector{T}, 𝓧, 𝓨::Vector{T})::T where T

    tmp_numerator = zero(T)
    tmp_denominator = zero(T)
    for k = 1:length(𝓧)

        tmp_numerator += 𝓨[k]*w[k]/wfunc(𝓧[k])

        tmp = w[k]/wfunc(𝓧[k])
        if iseven(k)
            tmp_denominator += tmp
        else
            tmp_denominator -= tmp
        end

    end

    @assert tmp_denominator != zero(T)
    return tmp_numerator/tmp_denominator
end

function levellederror(wfunc::Function, 𝓧, 𝓨::Vector{T})::T where T
    𝑤 = get𝑤scaled(𝓧)
    return levellederror(wfunc, w, 𝓧, 𝓨)
end

# implements || f-g ||_\infty
function discretizeL∞norm(f::Function, g::Function, xq::Vector{T})::T where T
    return maximum( abs(f(xq[i])-g(xq[i])) for i = 1:length(xq))
end

function discretizeL∞normparallel!(𝑒_xq::SharedArray{Float64},
                f::Function, g::Function, xq::Vector{T}) where T

    @assert length(𝑒_xq) == length(xq)

    @sync @distributed for i = 1:length(xq)
        𝑒_xq[i] = convert(Float64, f(xq[i])-g(xq[i]) )
    end

    k = argmax(abs.(𝑒_xq))
    return abs(𝑒_xq[k]), k, 𝑒_xq[k]
end

function discretizeL∞normparallel(f::Function, g::Function, xq::Vector{T}) where T

    𝑒_xq = SharedArray{Float64}(length(xq))

   return discretizeL∞normparallel!(𝑒_xq, f, g, xq)
end

function discretizeL∞normparallel(𝑒::Function, xq::Vector{T}) where T

    𝑒_xq = SharedArray{Float64}(length(xq))

    @sync @distributed for i = 1:length(xq)
        𝑒_xq[i] = convert(Float64, 𝑒(xq[i]) )
    end

    k = argmax(abs.(𝑒_xq))

    return abs(𝑒_xq[k]), k, 𝑒_xq[k]
end

### take 2.
function geth(  𝑤::Vector{T},
                wfunc::Function,
                dfunc::Function,
                𝓧)::T where T

    numerator = zero(T)
    denominator = zero(T)
    for i = 1:length(𝓧)
        k = i-1
        numerator += 𝑤[i]*dfunc(𝓧[i])
        denominator += (-1)^k*𝑤[i]/wfunc(𝓧[i])
    end
    out = numerator/denominator

    @assert isfinite(out)

    return out
end

function getp(x, wfunc::Function, dfunc::Function, 𝓧)
    𝑓 = get𝑓(wfunc,dfunc,𝓧)
    𝑤 = get𝑤scaled(𝓧)

    return wfunc(x)*Barycentric2nditp(x,𝑤,𝓧,𝑓)
end

function get𝑓(wfunc::Function, dfunc::Function, 𝓧::Vector{T}) where T

    𝑤 = get𝑤scaled(𝓧)
    h = geth(𝑤, wfunc, dfunc, 𝓧)

    return collect( dfunc(𝓧[i]) - (-1)^(i-1)*h/wfunc(𝓧[i]) for i = 1:length(𝓧))
end

# efficient version if 𝑤 and h were computed already.
function getp(x, 𝑤::Vector{T}, h::T, wfunc::Function, dfunc::Function, 𝓧::Vector{T}) where T

    𝑓 = collect( dfunc(𝓧[i]) - (-1)^(i-1)*h/wfunc(𝓧[i]) for i = 1:length(𝓧))

    return wfunc(x)*Barycentric2nditp(x,𝑤,𝓧,𝑓)
end

# even more efficient version, when 𝑓 is already computed.
function getp(x, 𝑓::Vector{T}, 𝑤::Vector{T}, h::T, wfunc::Function, dfunc::Function, 𝓧::Vector{T}) where T

    return wfunc(x)*Barycentric2nditp(x,𝑤,𝓧,𝑓)
end

function raisedcosinefreqrsp(x::T, β, 𝑇)::T where T
    return raisedcosinefreqrsp(x,convert(T,β),convert(T,𝑇))
end

function raisedcosinefreqrsp(x::T, β::T, 𝑇::T)::T where T
    two = one(T) + one(T)
    low = (one(T)-β)/(two*𝑇)
    high = (one(T)+β)/(two*𝑇)

    if abs(x) <= low
        return one(T)
    elseif low < abs(x) <= high
        return one(T)/two*( one(T) + cos(π*𝑇/β *(abs(x)-low)) )
    end

    return zero(T)
end

# generates L+2 positions in [-1,1]
function getinitialpositions(   L::Int,
                                candidate_multiple::Int)
    #
    N_candidates = candidate_multiple*L

    # generate canadiate positions. Take them as Chebyshev nodes.
    𝓧 = collect( Chebyshev2ndnode(k,N_candidates-1) for k = 0:N_candidates-1 )

    # select L positions.
    V = getunivariateChebyVmatrix(𝓧, L+1)

    B_unused, 𝓲 = maxvolviaqr(V)

    𝓧_chosen = 𝓧[𝓲][1:L+2]

    sort!(𝓧_chosen)

    return 𝓧_chosen, 𝓲
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
mutable struct FactorProportionStatistic <: MyXMLElement
    el::XMLOrNothing
    factor_model::IntegratedFactorsXMLElement
    trait_model::TraitLikelihoodXMLElement
    id::String

    function FactorProportionStatistic(factor_model::IntegratedFactorsXMLElement,
                                        trait_model::TraitLikelihoodXMLElement,
                                        id::String)
        return new(nothing, factor_model, trait_model, id)
    end
end

function name(::FactorProportionStatistic)
    return bn.FACTOR_PROPORTION_STATISTIC
end

function make_xml(fps::FactorProportionStatistic)
    el = new_element(name(fps))

    set_attribute(el, bn.ID, fps.id)
    add_ref_el(el, fps.factor_model)
    add_ref_el(el, fps.trait_model)

    fps.el = el
    return el
end

function get_loggables(fps::FactorProportionStatistic)
    return fps
end


text:
@testset "MultiMPIFile submodule" begin

fnMeasBruker = "./data/measurement"
fnMeasMultiV2 = "./data/mdf/measurement_multi_V2.mdf"

measBruker = MultiMPIFile([fnMeasBruker,fnMeasBruker,fnMeasBruker])

saveasMDF(fnMeasMultiV2, measBruker)

mdfv2 = MPIFile(fnMeasMultiV2)
@test typeof(mdfv2) <: MDFFileV2

for mdf in (measBruker,mdfv2)
  @info "Test $mdf"
  @test studyName(mdf) == "Wuerfelphantom"
  @test studyNumber(mdf) == 1
  @test studyDescription(mdf) == "n.a."
  @test studyTime(mdf) == DateTime( "2015-09-15T10:21:10.992" )

  @test experimentName(mdf) == "fuenf (E18)"
  @test experimentNumber(mdf) == 18
  @test experimentDescription(mdf) == "fuenf (E18)"
  @test experimentSubject(mdf) == "Wuerfelphantom"
  @test experimentIsSimulation(mdf) == false
  @test experimentIsCalibration(mdf) == false

  @test scannerFacility(mdf) == "Universitätsklinikum Hamburg Eppendorf"
  @test scannerOperator(mdf) == "nmrsu"
  @test scannerManufacturer(mdf) == "Bruker/Philips"
  @test scannerName(mdf) == "Preclinical MPI System"
  @test scannerTopology(mdf) == "FFP"

  @test tracerName(mdf) == ["Resovist"]
  @test tracerBatch(mdf) == ["0"]
  @test tracerVendor(mdf) == ["n.a."]
  @test tracerVolume(mdf) == [0.0]
  @test tracerConcentration(mdf) == [0.5]
  @test tracerInjectionTime(mdf) == [DateTime("2015-09-15T11:17:23.011")]

  @test acqStartTime(mdf) == DateTime("2015-09-15T11:17:23.011")
  @test acqGradient(mdf)[:,:,1,1] == [-1.25 0 0; 0 -1.25 0;0 0 2.5]
  @test acqFramePeriod(mdf) == 0.9792000000000001
  @test size(acqOffsetFieldShift(mdf)) == (3, 1, 1500)

  @test dfNumChannels(mdf) == 3
  @test dfWaveform(mdf) == "sine"
  @test dfStrength(mdf)[:,:,1] == [0.014 0.014 0.0]
  @test dfPhase(mdf)[:,:,1] == [1.5707963267948966 1.5707963267948966 1.5707963267948966]
  @test dfBaseFrequency(mdf) == 2500000.0
  @test dfDivider(mdf)[:,1] == [102; 96; 99]
  @test dfCycle(mdf) == 6.528E-4

  @test rxNumChannels(mdf) == 3
  @test rxBandwidth(mdf) == 1250000.0
  @test rxNumSamplingPoints(mdf) == 1632
  @test rxDataConversionFactor(mdf) == repeat([1.0, 0.0], outer=(1,rxNumChannels(mdf)))
  @test acqNumAverages(mdf) == 1

  @test acqNumFrames(mdf) == 1
  @test acqNumPeriodsPerFrame(mdf) == 1500
  @test acqNumPeriods(mdf) == 1500
  @test acqNumPatches(mdf) == 1
  @test acqNumPeriodsPerPatch(mdf) == 1500

  @test size( measData(mdf) ) == (1632,3,1500,1)
  @test size( measDataTDPeriods(mdf) ) == (1632,3,1500)
  @test size( measDataTDPeriods(mdf, 1001:1100) ) == (1632,3,100)

  N = acqNumFrames(mdf)

  @test size(getMeasurements(mdf, numAverages=1,
             spectralLeakageCorrection=false)) == (1632,3,1500,1)

  @test size(getMeasurements(mdf, numAverages=1, averagePeriodsPerPatch=true,
             spectralLeakageCorrection=false)) == (1632,3,1,1)

  #=@test size(getMeasurements(mdf, numAverages=10,
              spectralLeakageCorrection=false, fourierTransform=false)) == (1632,3,3,50)

  @test size(getMeasurements(mdf, numAverages=10, frames=1:100,
              fourierTransform=true)) == (817,3,3,10)

  @test size(getMeasurements(mdf, numAverages=10, frames=1:100,
              fourierTransform=true, loadasreal=true)) == (1634,3,3,10)

  @test size(getMeasurements(mdf,frequencies=1:10, numAverages=10)) == (10,3,50)=#

end

end


text:
using Plots
rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

linspace(a,b,N)= (a+ (b-a)*(i)/(N-1) for i=0:(N-1))

x = y = broadcast(exp,collect(linspace(log(0.5), log(1.5), 100)))

dpi = 600

l = Vector{Float64}()
for i in x
   for j in y
       v = [i j]
       append!(l, 1/rosenbrock(v))
   end
end

surface(x,y,l, dpi=dpi)
png("./inverse_rosenbrock_surface.png")
contour(x,y,l, dpi=dpi)
png("./inverse_rosenbrock_contour.png")


text:
module FETermsTests

using Test
using Gridap.Arrays
using Gridap.Algebra
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Integration
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData
using LinearAlgebra

domain =(0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 1
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]

u_sol(x) = x[1]+x[2]
f(x) = 0

dirichlet_tags = "boundary"
V = GradConformingFESpace(reffes,model,dirichlet_tags)
U = TrialFESpace(V,u_sol)

trian = get_triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

strian = SkeletonTriangulation(model)
degree = 2
squad = CellQuadrature(strian,degree)

assem = SparseMatrixAssembler(U,V)

uhd = zero(U)
v = get_cell_basis(V)
u = get_cell_basis(U)
uh = interpolate(u_sol,U)

a(u,v) = ∇(v)⊙∇(u)
l(v) = f*v

j(u,du,v) = a(du,v)
r(u,v) = a(u,v) - l(v)

z(u,v) = jump(v)*jump(u)
w(v) = jump(v)*f

t_linear = LinearFETerm(z,strian,squad)

t_source = FESource(w,strian,squad)

# AffineFETerm

t_affine = AffineFETerm(a,l,trian,quad)

matdata = collect_cell_matrix(u,v,[t_affine, t_linear])
vecdata = collect_cell_vector(uhd,v,[t_affine, t_linear, t_source])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ b
@test x ≈ get_free_values(uh)

data = collect_cell_matrix_and_vector(uhd,u,v,[t_affine, t_linear, t_source])
A, b = allocate_matrix_and_vector(assem,data)
assemble_matrix_and_vector!(A,b,assem,data)
x = A \ b
@test x ≈ get_free_values(uh)
A, b = assemble_matrix_and_vector(assem,data)
x = A \ b
@test x ≈ get_free_values(uh)

matdata = collect_cell_jacobian(uh,u,v,[t_affine])
vecdata = collect_cell_residual(uh,v,[t_affine])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ -b
@test (x.+1) ≈ ones(length(x))

# LinearFETerm and FESource

t_linear = LinearFETerm(a,trian,quad)
t_source = FESource(l,trian,quad)

matdata = collect_cell_matrix(u,v,[t_linear, t_source])
vecdata = collect_cell_vector(uhd,v,[t_linear, t_source])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ b
@test x ≈ get_free_values(uh)

matdata = collect_cell_jacobian(uh,u,v,[t_linear, t_source])
vecdata = collect_cell_residual(uh,v,[t_linear, t_source])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ -b
@test (x.+1) ≈ ones(length(x))

# NonlinearFETerm

t_nonlinear = FETerm(r,j,trian,quad)

matdata = collect_cell_jacobian(uh,u,v,[t_nonlinear])
vecdata = collect_cell_residual(uh,v,[t_nonlinear])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ -b
@test (x.+1) ≈ ones(length(x))

data = collect_cell_jacobian_and_residual(uh,u,v,[t_nonlinear])
A, b = allocate_matrix_and_vector(assem,data)
assemble_matrix_and_vector!(A,b,assem,data)
x = A \ -b
@test (x.+1) ≈ ones(length(x))

# AffineFETermFromCellMatVec

q = get_coordinates(quad)
w_q = get_weights(quad)
ϕ = get_cell_map(trian)
jac = ∇(ϕ)
jac_q = evaluate(jac,q)
x_q = evaluate(ϕ,q)

function matvecfun(u,v)
  cellmat = integrate(∇(u)⋅∇(v),trian,quad)
  cellvec = integrate(v*f,trian,quad)
  pair_arrays(cellmat,cellvec)
end

t_matvec_Ω = AffineFETermFromCellMatVec(matvecfun,trian)

matdata = collect_cell_matrix(u,v,[t_matvec_Ω,])
vecdata = collect_cell_vector(uhd,v,[t_matvec_Ω,])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ b
@test x ≈ get_free_values(uh)

data = collect_cell_matrix_and_vector(uhd,u,v,[t_matvec_Ω,])
A, b = allocate_matrix_and_vector(assem,data)
assemble_matrix_and_vector!(A,b,assem,data)
x = A \ b
@test x ≈ get_free_values(uh)

function jacresfun(uh,du,v)
  cellmat = integrate(∇(du)⋅∇(v),trian,quad)
  cellvec = integrate(∇(uh)⋅∇(v)-v*f,trian,quad)
  pair_arrays(cellmat,cellvec)
end

t_jacres_Ω = FETermFromCellJacRes(jacresfun,trian)

matdata = collect_cell_jacobian(uh,u,v,[t_jacres_Ω])
vecdata = collect_cell_residual(uh,v,[t_jacres_Ω])
A = assemble_matrix(assem,matdata)
b = assemble_vector(assem,vecdata)
x = A \ -b
@test (x.+1) ≈ ones(length(x))

data = collect_cell_jacobian_and_residual(uh,u,v,[t_jacres_Ω])
A, b = allocate_matrix_and_vector(assem,data)
assemble_matrix_and_vector!(A,b,assem,data)
x = A \ -b
@test (x.+1) ≈ ones(length(x))

end # module


text:
function dtduftau(t,
		  fft_u, fft_v,
		  fft_du, fft_dv,
                  A1, A2,
		  matr, conjmatr,
		  sigma, llambda, Ntaumm)

    sigma=1

    u = ifft(exp.(1im * t .* A1) .* fft_u .* matr)
    v = ifft(exp.(1im * t .* A1) .* fft_v .* matr)

    dtu = ifft(exp.(1im * t .* A1) .* (1im*A1) .* fft_u) .* matr
    dtv = ifft(exp.(1im * t .* A1) .* (1im*A1) .* fft_v) .* matr

    du = ifft(exp.(1im * t .* A1) .* fft_du .* matr)
    dv = ifft(exp.(1im * t .* A1) .* fft_dv .* matr)

    dtdu=ifft(exp.(1im * t .* A1) * (1im*A1) .* fft_du)*matr
    dtdv=ifft(exp.(1im * t .* A1) * (1im*A1) .* fft_dv)*matr

    z    = (u    + conj(v)    ) / 2
    dtz  = (dtu  + conj(dtv)  ) / 2
    dz   = (du   + conj(dv)   ) / 2
    dtdz = (dtdu + conj(dtdv) ) / 2

    fz1  = (  2 * dtz * conj(z) * dz + 2 * z * conj(dtz) * dz
            + 2 * abs(z)^2 * dtdz + 2 * z * dtz * conj(dz)
	    + z^2 * conj(dtdz))

    champu1 = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* fft(conjmatr .* fz1)
    champv1 = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* fft(conjmatr .* conj(fz1))

    fz1     = 2 * abs(z)^2 * dz + z^2 * conj(dz)
    champu2 = - 1im * llambda * A2 * exp.(-1im * t * A1) .* (-1im*A1) .* fft(conjmatr .* fz1)
    champv2 = - 1im * llambda * A2 * exp.(-1im * t * A1) .* (-1im*A1) .* fft(conjmatr .* conj(fz1))

    champu = champu1 .+ champu2
    champv = champv1 .+ champv2

    champu, champv

end


text:
# Implement experiments as described in the document "Research update 3_18_29.docx"
export avg_mutational_robustness, run_avg_mutational_robustness
using DataFrames
using CSV
using Statistics

#=
numinputs = 2
numoutputs = 2
nodearity = 2
numinteriors = 6
numlevelsback = 4
context = construct_contexts(numinputs)[numinputs]
p = Parameters(numinputs, numoutputs, nodearity, numinteriors, numlevelsback)
=#
#iterations = 20
#print_parameters(p)

function run_avg_mutational_robustness( iterations::Integer, numinputs::IntRange,
    numoutputs::IntRange, numinteriors::IntRange, numlevelsback::IntRange,
    funcs::Vector{Func}, csvfile::String; active_only::Bool=false )
  nodearity = 2
  df = DataFrame()
  df.numinputs=Int64[]
  df.numoutputs=Int64[]
  df.numinteriors=Int64[]
  df.numlevelsback=Int64[]
  df.mut_robustness=Float64[]
  df.redundancy=Float64[]
  df.complexity=Float64[]
  df.degeneracy=Float64[]
  df.fract_active=Float64[]
  for num_inputs = numinputs
    for num_outputs = numoutputs
      for num_interiors = numinteriors
        for num_levelsback = numlevelsback
          p = Parameters( num_inputs, num_outputs, nodearity, num_interiors, num_levelsback )
          #print_parameters( p )
          c = random_chromosome( p, funcs )
          (iterations,robustness,redundancy,complexity,degeneracy,fract_active) = avg_mutational_robustness( c, p, funcs,
              iterations, active_only=true )
          #println("robustness: ",robustness,"  fract_active: ",fract_active)
          new_row = (num_inputs,num_outputs,num_interiors,num_levelsback,robustness,redundancy,complexity,degeneracy,fract_active)
          Base.push!( df, new_row )
        end
      end
    end
  end
  df
  open( csvfile, "w" ) do f
    println(f,"funcs: ", funcs )
    println(f,"nodearity: ",nodearity)
    println(f,"iterations: ",iterations)
    CSV.write( f, df, append=true, writeheader=true )
  end
  println(df)
end

# Currently uses fract_active only in correlation
function avg_mutational_robustness( c::Chromosome, p::Parameters, funcs::Vector{Func}, iterations::Integer;
    active_only::Bool=false )
  #println("active_only: ",active_only," in:",p.numinputs," out:",p.numoutputs," ints:",p.numinteriors)
  base = 2.0
  context = construct_context( p.numinputs )
  sum_redundancy = 0.0
  sum_complexity = 0.0
  sum_degeneracy = 0.0
  sum_mutational_robustness = 0.0
  sum_fract_active = 0.0
  sum_num_mut_locations = 0
  # Lists are for the computation of correlations
  fract_active_list = fill(0.0,iterations)
  mutational_robustness_list = fill(0.0,iterations)
  for i = 1:iterations
    c = random_chromosome( p, funcs )
    redund = redundancy( c, base=base )
    complexity = complexity5( c, base=base )
    degen = degeneracy( c, base=base )
    fract_active = fraction_active(c)
    mut_robustness = mutational_robustness( c, funcs, active_only=active_only )
    sum_mutational_robustness += mut_robustness
    sum_redundancy += redund
    sum_complexity += complexity
    sum_degeneracy += degen
    sum_fract_active += fract_active
    #sum_num_mut_locations += num_mut_locs
    #println("i: ",i,"  average_fract_active: ",average_fract_active)
    fract_active_list[i] = fract_active
    #println("i: ",i,"  mutational_robustness: ",mutational_robustness)
    mutational_robustness_list[i] = mut_robustness
  end
  correlation = cor( mutational_robustness_list, fract_active_list )
  #println("correlation: ",correlation)
  #println("total mutational_robustness: ",sum_mutational_robustness/iterations)
  (iterations,sum_mutational_robustness/iterations,sum_redundancy/iterations,
      sum_complexity/iterations,sum_degeneracy/iterations,sum_fract_active/iterations)

end

#= Not used
# Computes the fraction of mutations that do not change the output
# If active_only==true, then only mutations that change active nodes are considered
# If active_only==false, then all mutations are considered
# Note that results are partially random since the result of mutations at a specific location can be different.
function mutational_robustness( c::Chromosome, funcs::Vector{Func}; active_only::Bool=false )
  context = construct_context( c.params.numinputs )
  prev_out = execute_chromosome( c, context )
  count_no_change = 0
  count_no_change = 0.0
  num_mut_locs = 0
  for j = 1:num_mutate_locations( c, funcs )
    cc = deepcopy(c)
    (cc,active) = mutate_chromosome!( cc, funcs, j )
    out_c = execute_chromosome( cc, context )
    #println("j: ",j,"  active: ",active)
    #print_chromosome( cc )
    if !active_only || active # consider all mutations if !active_only, only active mutations if active_only
      #out_c = execute_chromosome( cc, context )
      count_no_change += ((out_c == prev_out) ? 1 : 0)
      if active_only && active
        num_mut_locs += 1
      end
    end
    #println("j: ",j,"  active: ",active,"  no_change: ",(out_c == prev_out),"  count_no_change: ",count_no_change)
  end
  #println("num_mut_locs: ",num_mut_locs,"  num_mutate_locations(): ",num_mutate_locations(c,funcs),"  count_no_change: ",count_no_change)
  #@assert num_mut_locs == num_mutate_locations(c,funcs)
  if !active_only
    num_mut_locs = num_mutate_locations( c, funcs )
  end
  count_no_change/num_mut_locs
end
=#


text:
using DensityRatioEstimation
using Distributions
using LinearAlgebra
using Statistics
using Optim
using JuMP, Ipopt
using Convex, ECOS
using Plots, VisualRegressionTests
using Test, Pkg, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

# helper funcions
include("utils.jl")

# simple cases for testing
pair₁ = Normal(1,1), Normal(0,2)
pair₂ = MixtureModel([Normal(-2,1),Normal(2,2)], [0.2,0.8]), Normal(0,2)

# list of tests
testfiles = [
  "basic.jl",
  "kmm.jl",
  "kliep.jl",
  "lsif.jl"
]

@testset "DensityRatioEstimation.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
function test_square_loss_layer(backend::Backend, T, eps)
  println("-- Testing SquareLossLayer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  tensor_dim = abs(rand(Int)) % 4 + 2
  dims = tuple((abs(rand(Int,tensor_dim)) % 6 + 6)...)
  println("    > $dims")
  preds = rand(T, dims)
  labels = rand(T, dims)

  ############################################################
  # Setup
  ############################################################
  layer  = SquareLossLayer(; bottoms=[:predictions, :labels])
  pred_blob  = make_blob(backend, T, dims)
  label_blob = make_blob(backend, T, dims)
  diff_blob  = make_blob(backend, T, dims)
  copy!(pred_blob, preds)
  copy!(label_blob, labels)
  inputs = Blob[pred_blob, label_blob]
  diffs = Blob[diff_blob, NullBlob()]

  state  = setup(backend, layer, inputs, diffs)

  forward(backend, state, inputs)

  loss = 0.5*vecnorm(preds-labels)^2 / dims[end]
  @test -eps < loss-state.loss < eps

  backward(backend, state, inputs, diffs)
  grad = (preds - labels) / dims[end]
  diff = similar(grad)
  copy!(diff, diffs[1])

  @test all(-eps .< grad - diff .< eps)

  shutdown(backend, state)
end

function test_square_loss_layer(backend::Backend)
  test_square_loss_layer(backend, Float32, 1e-2)
  test_square_loss_layer(backend, Float64, 1e-8)
end

if test_cpu
  test_square_loss_layer(backend_cpu)
end
if test_gpu
  test_square_loss_layer(backend_gpu)
end


text:
using DrWatson
@quickactivate "OptimalEstimationProject"
using MATLAB
using IndirectTrajOpt
using StaticArrays
using LinearAlgebra
using DataFrames

## CORRESPONDS TO FIRST ME TRAJ GENERATED!
# Trajectory time span
#tspan = (0.0, 26.0)
# ME Initial CoStates
#λ0 = [-0.8778526555553886,-1.5568310993815173,-0.4120913077312581,0.005878876387420351,-0.002181056020767603,-0.00035144027756129617,0.0941110364795588]

## CORRESPONDS TO UPDATED TOF TRAJ with ϵ = 0.16662720666508826
tspan = (0.0, 33.1)
λ0 = [5.9258156810825735, 12.938611585295579, 1.3029560018476296, -0.04039806526574376, 0.018486336322689402, -0.00048478964517609035, 0.11840864999974321]

# Initial conditions
ics = [-0.019488511458668, -0.016033479812051, 0.0,
        8.918881923678198, -4.081793688818725, 0.0,
        1.0]

# Init parameters
ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
μ = ps.crp.μ
LU = ps.crp.LU
TU = ps.crp.TU

# Target final conditions
fcs = [(1.0 - μ) + 6060.483/LU, 19452.284/LU, -34982.968/LU, 0.082677*(TU/LU), 0.006820*(TU/LU), -0.368434*(TU/LU), 0.0]
xf = @SVector [fcs[1], fcs[2], fcs[3], fcs[4], fcs[5], fcs[6]]

# Generate homotopy parameter vector
N = 50
homotopyParamVec = [(j^2 - 1.)/(N^2 - 1.) for j in N:-1:1]

# Scale vector if not starting at ϵ = 1.0
homotopyParamVec *= 0.16662720666508826

# Solve with homotopy
prob = IndirectOptimizationProblem("Low Thrust 10 CR3BP", ics, fcs, tspan)
fss = FSSSolver(λ0, ics, fcs, prob.BVPFunc, prob.BVPWithSTMFunc;
    homotopy = true, homotopyParamVec = homotopyParamVec)
solve!(fss; factor = 1.0, ftol = 1e-10)

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
include("polyagamma_sampler.jl")

function _mcmc_iter_pg(
    latent::AbstractLatent,
    parameters::Vector{<:AbstractParameters},
    responses_val::Vector{Float64},
    W_val::Vector{Float64};
    sampling = true
    )
        latent.posterior = __posterior(latent, parameters, responses_val, W_val)
        _chain_append_and_set_val!(latent; sampling = sampling)
    return latent::AbstractLatent
end
function _mcmc_iter_pg(
    latent::AbstractLatent,
    parameters::Vector{<:AbstractParameters},
    responses_val::Vector{Float64},
    W_val::SharedArrays.SharedVector{Float64};
    sampling = true
    )
        latent.posterior = __posterior(latent, parameters, responses_val, W_val)
        _chain_append_and_set_val!(latent; sampling = sampling)
    return latent::AbstractLatent
end

function mcmc_iter_pg!(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem},
    responses::Vector{<:AbstractResponse},
    W::Vector{PolyaGammaSample};
    sampling = true,
    already_sorted = false,
)
    if !already_sorted
        sort!(items, by = i -> i.idx)
        sort!(responses, by = r -> r.examinee_idx)
    end
    update_posterior!(examinee, items, responses, W; already_sorted = true)
    set_val_from_posterior!(examinee; sampling = sampling)
end

function update_estimate!(item::AbstractItem; sampling = true)
    _update_estimate!(item.parameters; sampling = sampling)
end

function _assess_examinee_pg(
    latent::AbstractLatent,
    parameters::Vector{<:AbstractParameters},
    responses::Vector{Float64};
    mcmc_iter = 2_000,
    sampling = true
    )
    # same as, but without val update (faster)
    for iter = 1:mcmc_iter
        W = map( p -> _generate_w(latent, p), parameters)
        latent = _mcmc_iter_pg(latent, parameters, responses, W; sampling = sampling)
    end
    return latent::AbstractLatent
end

function assess_examinee_pg!(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem},
    responses::Vector{Response};
    mcmc_iter = 2_000,
    sampling = true,
    already_sorted = false,
    kwargs...
)
    if !already_sorted
        sort!(responses, by = r -> r.item_idx)
        sort!(items, by = i -> i.idx)
    end
    latent = examinee.latent
    _empty_chain!(latent)
    # chain = SharedArrays.SharedArray{Float64}(mcmc_iter)
    # # extract `mcmc_iter` samples from the polyagamma and from theta conditional posterior
    # Distributed.@sync Distributed.@distributed for iter in 1:mcmc_iter
    #     W = generate_w(examinee, items)
    #     chain[iter] = rand(posterior(examinee, items, responses, W))
    # end
    # examinee.latent.chain = copy(chain)

    # same as, but without val update (faster)
    latent = _assess_examinee_pg(
        latent,
        get_parameters(items),
        map( r -> r.val, responses),
        mcmc_iter = mcmc_iter,
        sampling = sampling,
        )
    update_estimate!(examinee; sampling = sampling)
    return nothing
end


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: AbstractRNG, gen_seed, split_uint, seed_type

"""
```julia
AbstractXoroshiro128 <: AbstractRNG{UInt64}
```

The base abstract type for `Xoroshiro128`, `Xoroshiro128Star`, `Xoroshiro128Plus` and `Xoroshiro128StarStar`.
"""
abstract type AbstractXoroshiro128 <: AbstractRNG{UInt64} end

for (star, plus, starstar) in (
        (false, false, false),
        (true, false, false),
        (false, true, false),
        (false, false, true),
    )
    rng_name = Symbol(string("Xoroshiro128", star ? "Star" : plus ? "Plus" : starstar ? "StarStar" : ""))
    @eval begin
        mutable struct $rng_name <: AbstractXoroshiro128
            x::UInt64
            y::UInt64
            function $rng_name(seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
                r = new(0, 0)
                seed!(r, seed)
                r
            end
        end

        $rng_name(seed::Integer) = $rng_name(split_uint(seed % UInt128))

        @inline function xorshift_next(r::$rng_name)
            $(plus ? :(p = r.x + r.y)
                : starstar ? :(p = xorshift_rotl(r.x * 5, 7) * 9) : nothing)
            s1 = r.y ⊻ r.x
            r.x = xorshift_rotl(r.x, 24) ⊻ s1 ⊻ (s1 << 16)
            r.y = xorshift_rotl(s1, 37)
            $(star ? :(r.y * 2685821657736338717) :
              (plus || starstar) ? :(p) : :(r.y))
        end
    end
end

@inline seed_type(::Type{T}) where T <: AbstractXoroshiro128 = NTuple{2, UInt64}

function copyto!(dest::T, src::T) where T <: AbstractXoroshiro128
    dest.x = src.x
    dest.y = src.y
    dest
end

copy(src::T) where T <: AbstractXoroshiro128 = copyto!(T(), src)

==(r1::T, r2::T) where T <: AbstractXoroshiro128 = r1.x == r2.x && r1.y == r2.y

seed!(r::AbstractXoroshiro128, seed::Integer) = seed!(r, split_uint(seed % UInt128))
function seed!(r::AbstractXoroshiro128, seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
    all(==(0), seed) && error("0 cannot be the seed")
    r.x = seed[1]
    r.y = seed[2]
    xorshift_next(r)
    xorshift_next(r)
    r
end

@inline rand(r::AbstractXoroshiro128, ::Type{UInt64}) = xorshift_next(r)


text:
module Manifolds

import Base:
    Array,
    +,
    -,
    *,
    \,
    /,
    ^,
    ==,
    angle,
    axes,
    convert,
    copy,
    copyto!,
    dataids,
    eltype,
    exp,
    getindex,
    identity,
    inv,
    isapprox,
    length,
    log,
    one,
    promote_rule,
    setindex!,
    show,
    similar,
    size,
    transpose,
    zero
import Distributions: _rand!, support
import LinearAlgebra: cross, det, Diagonal, dot, mul!, norm, I, UniformScaling
import ManifoldsBase:
    allocate,
    allocate_result,
    allocate_result_type,
    allocation_promotion_function,
    array_value,
    base_manifold,
    check_manifold_point,
    check_tangent_vector,
    combine_allocation_promotion_functions,
    decorated_manifold,
    decorator_transparent_dispatch,
    default_decorator_dispatch,
    distance,
    embed,
    embed!,
    exp,
    exp!,
    get_basis,
    get_coordinates,
    get_coordinates!,
    get_embedding,
    get_vector,
    get_vector!,
    get_vectors,
    geodesic,
    injectivity_radius,
    inner,
    isapprox,
    is_manifold_point,
    is_tangent_vector,
    inverse_retract,
    inverse_retract!,
    log,
    log!,
    manifold_dimension,
    norm,
    number_eltype,
    project,
    project!,
    representation_size,
    retract,
    retract!,
    shortest_geodesic,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    zero_tangent_vector,
    zero_tangent_vector!
import Random: rand
import Statistics: mean, mean!, median, median!, std, var
import StatsBase: kurtosis, mean_and_std, mean_and_var, moment, skewness

using Base.Iterators: repeated
using Distributions
using Einsum: @einsum
using FiniteDifferences
using HybridArrays
using LinearAlgebra
using LightGraphs
using LightGraphs: AbstractGraph
using ManifoldsBase: CoTVector, Manifold, MPoint, TVector, DefaultManifold
using ManifoldsBase
using ManifoldsBase:
    AbstractDecoratorManifold,
    AbstractNumbers,
    @decorator_transparent_fallback,
    @decorator_transparent_function,
    @decorator_transparent_signature,
    _euclidean_basis_vector,
    _extract_val,
    hat,
    hat!,
    @invoke_maker,
    is_decorator_transparent,
    is_default_decorator,
    manifold_function_not_implemented_message,
    vee,
    vee!
using ManifoldsBase:
    AbstractBasis,
    AbstractOrthogonalBasis,
    DefaultOrDiagonalizingBasis,
    DiagonalizingBasisData,
    VeeOrthogonalBasis,
    AbstractOrthonormalBasis,
    DefaultOrthonormalBasis
using ManifoldsBase:
    ArrayCoTVector, ArrayManifold, ArrayMPoint, ArrayTVector, ArrayCoTVector
using ManifoldsBase: AbstractRetractionMethod, ExponentialRetraction
using ManifoldsBase: QRRetraction, PolarRetraction, ProjectionRetraction
using ManifoldsBase: AbstractInverseRetractionMethod, LogarithmicInverseRetraction
using ManifoldsBase: QRInverseRetraction, PolarInverseRetraction, ProjectionInverseRetraction
using ManifoldsBase: AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport
using ManifoldsBase: ℝ, ℂ, ℍ
using ManifoldsBase:
    AbstractEmbeddingType,
    AbstractIsometricEmbeddingType,
    TransparentIsometricEmbedding,
    DefaultIsometricEmbeddingType,
    DefaultEmbeddingType
using ManifoldsBase:
    AbstractEmbeddedManifold,
    EmbeddedManifold

using Markdown: @doc_str
using Random: AbstractRNG
using Requires
using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, get_weight
using StaticArrays
using StatsBase: AbstractWeights, UnitWeights, values, varcorrection

include("utils.jl")
include("autodiff.jl")
include("SizedAbstractArray.jl")

include("statistics.jl")

include("manifolds/VectorBundle.jl")

include("distributions.jl")
include("projected_distribution.jl")
include("product_representations.jl")

# It's included early to ensure visibility of `Identity`
include("groups/group.jl")

include("manifolds/MetricManifold.jl")
include("manifolds/ProductManifold.jl")
include("manifolds/PowerManifold.jl")
include("manifolds/GraphManifold.jl")

include("manifolds/Euclidean.jl")
include("manifolds/Lorentz.jl")

include("manifolds/CholeskySpace.jl")
include("manifolds/Circle.jl")
include("manifolds/FixedRankMatrices.jl")
include("manifolds/GeneralizedGrassmann.jl")
include("manifolds/GeneralizedStiefel.jl")
include("manifolds/Grassmann.jl")
include("manifolds/Hyperbolic.jl")
include("manifolds/Rotations.jl")
include("manifolds/SkewSymmetric.jl")
include("manifolds/Stiefel.jl")
include("manifolds/Sphere.jl")
include("manifolds/Symmetric.jl")
include("manifolds/SymmetricPositiveDefinite.jl")
include("manifolds/SymmetricPositiveDefiniteLinearAffine.jl")
include("manifolds/SymmetricPositiveDefiniteLogCholesky.jl")
include("manifolds/SymmetricPositiveDefiniteLogEuclidean.jl")

include("manifolds/Torus.jl")
include("manifolds/Oblique.jl")

include("groups/metric.jl")
include("groups/group_action.jl")
include("groups/group_operation_action.jl")
include("groups/array_manifold.jl")
include("groups/product_group.jl")
include("groups/semidirect_product_group.jl")

include("groups/translation_group.jl")
include("groups/special_orthogonal.jl")
include("groups/circle_group.jl")

include("groups/translation_action.jl")
include("groups/rotation_action.jl")

include("groups/special_euclidean.jl")

function __init__()
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
        using .ForwardDiff
        include("forward_diff.jl")
    end

    @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
        using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve
        include("ode.jl")
    end
end
#
export CoTVector, Manifold, MPoint, TVector, Manifold
export Euclidean,
    CholeskySpace,
    Circle,
    FixedRankMatrices,
    GeneralizedGrassmann,
    GeneralizedStiefel,
    Grassmann,
    Hyperbolic,
    Lorentz,
    Oblique,
    Rotations,
    SkewSymmetricMatrices,
    Sphere,
    Stiefel,
    SymmetricMatrices,
    SymmetricPositiveDefinite,
    Torus
export SVDMPoint, UMVTVector, AbstractNumbers, ℝ, ℂ, ℍ
# decorator manifolds
export AbstractDecoratorManifold
export ArrayManifold, ArrayMPoint, ArrayTVector, ArrayCoTVector
export CotangentBundle,
    CotangentSpaceAtPoint, CotangentBundleFibers, CotangentSpace, FVector
export AbstractPowerManifold,
    AbstractPowerRepresentation,
    ArrayPowerRepresentation,
    NestedPowerRepresentation,
    PowerManifold
export ProductManifold
export GraphManifold, GraphManifoldType, VertexManifold, EdgeManifold
export ProjectedPointDistribution, ProductRepr, TangentBundle, TangentBundleFibers
export TangentSpace, TangentSpaceAtPoint, VectorSpaceAtPoint, VectorSpaceType, VectorBundle
export VectorBundleFibers
export AbstractVectorTransportMethod, ParallelTransport, ProjectedPointDistribution
export AbstractEmbeddedManifold
export Metric,
    RiemannianMetric,
    LorentzMetric,
    EmbeddedManifold,
    EuclideanMetric,
    LinearAffineMetric,
    LogCholeskyMetric,
    LogEuclideanMetric,
    MinkowskiMetric,
    PowerMetric,
    ProductMetric,
    MetricManifold
export AbstractEmbeddingType, AbstractIsometricEmbeddingType
export DefaultEmbeddingType, DefaultIsometricEmbeddingType, TransparentIsometricEmbedding
export AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport
export AbstractRetractionMethod,
    ExponentialRetraction,
    QRRetraction,
    PolarRetraction,
    ProjectionRetraction
export AbstractInverseRetractionMethod,
    LogarithmicInverseRetraction,
    QRInverseRetraction,
    PolarInverseRetraction,
    ProjectionInverseRetraction
export AbstractEstimationMethod,
    GradientDescentEstimation,
    CyclicProximalPointEstimation,
    GeodesicInterpolation,
    GeodesicInterpolationWithinRadius
export
    CachedBasis,
    DefaultBasis,
    DefaultOrthogonalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    ProjectedOrthonormalBasis
export ×,
    allocate,
    allocate_result,
    base_manifold,
    bundle_projection,
    check_manifold_point,
    check_tangent_vector,
    christoffel_symbols_first,
    christoffel_symbols_second,
    christoffel_symbols_second_jacobian,
    complex_dot,
    decorated_manifold,
    det_local_metric,
    distance,
    einstein_tensor,
    embed,
    embed!,
    exp,
    exp!,
    flat,
    flat!,
    gaussian_curvature,
    geodesic,
    get_embedding,
    hat,
    hat!,
    incident_log,
    injectivity_radius,
    inner,
    inverse_local_metric,
    inverse_retract,
    inverse_retract!,
    isapprox,
    is_group_decorator,
    is_decorator_transparent,
    is_default_metric,
    is_default_decorator,
    is_manifold_point,
    is_tangent_vector,
    isapprox,
    inner,
    kurtosis,
    local_metric,
    local_metric_jacobian,
    log,
    log!,
    log_local_metric_density,
    manifold_dimension,
    metric,
    mean,
    mean!,
    mean_and_var,
    mean_and_std,
    median,
    median!,
    minkowski_metric,
    moment,
    norm,
    normal_tvector_distribution,
    number_eltype,
    one,
    power_dimensions,
    project,
    project!,
    projected_distribution,
    real_dimension,
    ricci_curvature,
    ricci_tensor,
    representation_size,
    retract,
    retract!,
    riemann_tensor,
    sharp,
    sharp!,
    shortest_geodesic,
    skewness,
    std,
    sym_rem,
    submanifold,
    submanifold_component,
    submanifold_components,
    uniform_distribution,
    var,
    vector_space_dimension,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    vee,
    vee!,
    zero_vector,
    zero_vector!,
    zero_tangent_vector,
    zero_tangent_vector!
# Lie group types & functions
export AbstractGroupAction,
    AbstractGroupOperation,
    AbstractGroupManifold,
    ActionDirection,
    AdditionOperation,
    CircleGroup,
    GroupManifold,
    GroupOperationAction,
    Identity,
    InvariantMetric,
    LeftAction,
    LeftInvariantMetric,
    MultiplicationOperation,
    ProductGroup,
    ProductOperation,
    RightAction,
    RightInvariantMetric,
    RotationAction,
    SemidirectProductGroup,
    SpecialEuclidean,
    SpecialOrthogonal,
    TranslationGroup,
    TranslationAction
export affine_matrix,
    apply,
    apply!,
    apply_diff,
    apply_diff!,
    base_group,
    center_of_orbit,
    has_approx_invariant_metric,
    compose,
    compose!,
    direction,
    g_manifold,
    group_exp,
    group_exp!,
    group_log,
    group_log!,
    has_biinvariant_metric,
    has_invariant_metric,
    identity,
    identity!,
    inv,
    inv!,
    invariant_metric_dispatch,
    inverse_apply,
    inverse_apply!,
    inverse_apply_diff,
    inverse_apply_diff!,
    inverse_translate,
    inverse_translate!,
    inverse_translate_diff,
    inverse_translate_diff!,
    make_identity,
    optimal_alignment,
    optimal_alignment!,
    screw_matrix,
    switch_direction,
    translate,
    translate!,
    translate_diff,
    translate_diff!
# Orthonormal bases
export AbstractBasis,
    AbstractOrthonormalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    ProjectedOrthonormalBasis,
    CachedBasis,
    DiagonalizingBasisData,
    ProductBasisData,
    PowerBasisData
export get_basis,
    get_coordinates,
    get_coordinates!,
    get_vector,
    get_vector!,
    get_vectors,
    number_system
end # module


text:
function forward_stagewise_regression(X, y; epsilon = 0.01, numIterations = 100)
  n, p = size(X)
  ny = length(y)
  @assert n == ny
  forward_stagewise_regression_inner(X, y, n, p, epsilon, numIterations)
end

function forward_stagewise_regression_inner(X, y, n, p, epsilon, numIterations)
  betas = zeros(p, numIterations)
  for k in 2:numIterations
    maxabsinnerproduct = 0.0
    bestinnerproduct = 0.0
    residuals = y .- X * betas[:, k-1]
    bestcolumn = rand(1:p)
    for j in 1:p
      innerproduct = sum( X[:,j] .* residuals )
      if abs(innerproduct) > maxabsinnerproduct
        bestinnerproduct = innerproduct
        maxabsinnerproduct = abs(innerproduct)
        bestcolumn = j
      end
    end
    betas[:, k] = betas[:, k-1]
    betas[bestcolumn, k] += epsilon * sign(bestinnerproduct)
  end
  betas
end

N = 20
P = 5
X = randn(N, P)
y = 1.0 * X[:,1] + 2.0 * X[:,3]
betas = forward_stagewise_regression(X, y; epsilon = 0.10)


text:
const NEXT_UNIQUE_ID = Ref{Int}(0)
next_unique_id() = (NEXT_UNIQUE_ID[] += 1; NEXT_UNIQUE_ID[])

make_name(id::Int, prefix="") = Symbol("$(prefix)x$id")
make_name(op::AbstractOp, prefix="") = Symbol("$(prefix)x$(op.id)")
make_name(name::String, prefix="") = Symbol("$(prefix)$(name)")

arg2expr(v::Variable, prefix="") = make_name(v.id, prefix)
arg2expr(s::Symbol, prefix="") = QuoteNode(s)
arg2expr(c, prefix="") = c

function to_expr(op::Call, prefix="")
    call = Expr(:call, [arg2expr(v, prefix) for v in (op.fn, op.args...)]...)
    return Expr(:(=), make_name(op.id, prefix), call)
end

to_expr(op::Constant, prefix="") = :($(make_name(op.id, prefix)) = $(op.val))

function loop_exit_tuple_expr_at_point(op::Loop, id::Int, prefix::String, loop_prefix::String)
    exit_name = make_name(op.id, prefix)
    arg_vars = loop_exit_vars_at_point(op, id)
    arg_names = [make_name(v.id, loop_prefix) for v in arg_vars]
    return Expr(:(=), exit_name, Expr(:call, tuple, arg_names...))
end

function to_expr(op::Loop, prefix="")
    loop_prefix = "l$(next_unique_id())"
    exprs = []
    # map parent input ids to continue ids
    init_var_names = []
    for (inp, parent) in zip(inputs(op.subtape), op.parent_inputs)
        init_var_name = make_name(inp.id, loop_prefix)
        push!(init_var_names, init_var_name)
        ex = Expr(:(=), init_var_name, make_name(parent.id, prefix))
        push!(exprs, ex)
    end
    # add exit tuple which will be used in case of zero trip count
    init_exit_tuple_ex = loop_exit_tuple_expr_at_point(op, 0, prefix, loop_prefix)
    push!(exprs, init_exit_tuple_ex)
    loop_ex = :(while true end)
    body = loop_ex.args[2]
    for (id, subop) in enumerate(op.subtape)
        if !isa(subop, Input)
            subex = to_expr(subop, loop_prefix)
            if subex isa Vector
                push!(body.args, subex...)
            else
                push!(body.args, subex)
            end
            if id == op.condition.id
                exit_expr = :(if !$(make_name(op.condition.id, loop_prefix)) end)
                exit_body = exit_expr.args[2]
                # update exit tuple
                exit_tuple_ex = loop_exit_tuple_expr_at_point(op, id, prefix, loop_prefix)
                push!(exit_body.args, exit_tuple_ex)
                push!(exit_body.args, Expr(:break))
                push!(body.args, exit_expr)
            end
        end
    end
    # map continue vars to inputs
    for (inp, cont) in zip(inputs(op.subtape), op.cont_vars)
        ex = Expr(:(=), make_name(inp.id, loop_prefix), make_name(cont.id, loop_prefix))
        push!(body.args, ex)
    end
    push!(exprs, loop_ex)
end

"""
    to_expr(tape::Tape)

Generate a Julia expression corresponding to the tape.
"""
function to_expr(tape::Tape)
    fn_name = gensym("tape_$(tape[V(1)].val)")
    header = Expr(:call, fn_name)
    for v in inputs(tape)
        op = tape[v]
        push!(header.args, Expr(:(::), make_name(op), op.typ))
    end
    body = Expr(:block)
    for op in tape
        op isa Input && continue
        ex = to_expr(op)
        if ex isa Vector
            push!(body.args, ex...)
        else
            push!(body.args, ex)
        end
    end
    push!(body.args, Expr(:return, make_name(tape.result.id)))
    fn_ex = Expr(:function, header, body)
    return fn_ex
end

"""
    compile(tape::Tape)

Compile tape into a normal Julia function.
"""
compile(tape::Tape) = Base.eval(@__MODULE__, to_expr(tape))


text:
using Base.Test

source_dir = pwd()

include(joinpath(source_dir, "../src", "algorithms", "nn.jl"));
include(joinpath(source_dir, "../src", "utils", "parse.jl"));

# write your own tests here
function test_parser()

    ## assume for now that we have data in the "data" directory as given
    ## to us from the assignment

    pos_seq_fp = "../data/rap1-lieb-positives.txt";
    total_seq_fp = "../data/yeast-upstream-1k-negative.fa";

    tdata, labels = parse_training(pos_seq_fp, total_seq_fp; balance=1);

    @test size(tdata, 1) == size(labels, 1);

    @test size(tdata, 2) == 68;

end

function test_basic_autoencoder()

    results = encode_8x3x8(100000, 0.05, ""; write_output = false, verbose = true)

    # make sure results are the right dimension
    @test size(results[1], 1) == 8
    @test size(results[1], 2) == 8

    @test size(results[3], 1) == 8
    @test size(results[3], 2) == 3

    #make sure we were able to train correctly with these parameters
    @test results[2] == one_hot_8x3x8(8)

end

function test_feedforward()

    tdata = one_hot_8x3x8(8)

    tdata = hcat(ones(8, 1), tdata)

    labels = one_hot_8x3x8(8)

    # draw random weights for first layer to hidden layaer
    w1 = rand(9, 3)
    # draw random weights for hidden layer to output layer
    w2 = rand(4, 8)

    # make sure feedforward works correctly
    hidden, pred = feedforward(tdata[1,:], w1, w2)

    # number of columns should correspond to number of hidden nodes + bias term
    @test size(hidden, 2) == 4

    # remove the bias term dimension
    @test size(pred, 2) == size(tdata, 1)

    # make sure that we're outputing probabilities
    @test sum(pred .<= 1) == size(pred, 2)
    @test sum(hidden .<= 1) == size(hidden, 2)

end

function test_backprop()

    tdata = one_hot_8x3x8(8);

    tdata = hcat(ones(8, 1), tdata)

    labels = one_hot_8x3x8(8)

    # draw random weights for first layer to hidden layaer
    w1 = rand(9, 3)
    # draw random weights for hidden layer to output layer
    w2 = rand(4, 8)

    # make sure feedforward works correctly
    hidden, pred = feedforward(tdata[1,:], w1, w2)

    u1, u2 = backpropogate(pred, lab, hidden, w1, w2)
    u1 = u1[1:3]

    u1 = (u1.' * inp.').'
    u2 = (u2.' * hidden).'

    # make sure that the weight can be updated correctly
    @test size(u1) == w1;
    @test size(u2) == w2;

end

function test_threshold()

    tdata = one_hot_8x3x8(8);

    tdata = hcat(ones(8, 1), tdata)

    labels = one_hot_8x3x8(8)

    # draw random weights for first layer to hidden layaer
    w1 = rand(9, 3)
    # draw random weights for hidden layer to output layer
    w2 = rand(4, 8)

    # make sure feedforward works correctly
    hidden, pred = feedforward(tdata[1,:], w1, w2)

    pred = threshold(pred)

    # make sure we're outputing either zeros or ones here
    @test sum(pred .== 1 || pred .== 0) == size(pred, 2)

end

test_parser()
test_basic_autoencoder()
test_feedforward()


text:
using Pipe

const data = parse.(Int, reduce(hcat, collect.(readlines())))

neighbours(pos) = @pipe [
    CartesianIndex(Tuple(pos) .- (i, j))
    for i in (-1, 0, 1)
    for j in (-1, 0, 1)
    if (i, j) != (0, 0)
] |> filter(_) do I
    x, y = Tuple(I)
    x > 0 && y > 0 && x ≤ 10 && y ≤ 10
end

function step(grid)
    grid = copy(grid)
    grid .+= 1
    while any(grid .> 9)
        pos = argmax(grid)
        grid[pos] = typemin(Int64)

        nei = neighbours(pos)
        grid[nei] .+= 1
    end
    c = count(<(0), grid)
    grid[grid .< 0] .= 0

    grid, c
end

function step(grid, n)
    s = 0
    for i in 1:n
        grid, c = step(grid)
        s += c
    end
    grid, s
end

# Part One
step(data, 100)[2] |> println

# Part Two
let grid = data
    for i in 1:typemax(Int64)
        grid, c = step(grid)
        if c == 100
            println(i)
            break
        end
    end
end

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("OrdinaryDiffEq")
Pkg.add("DiffEqBase")
Pkg.add("PyCall")
using DifferentialEquations
using PyCall

max_stars_repo_name: jlapeyre/diffeqpy

max_issues_repo_name: jlapeyre/diffeqpy

max_forks_repo_name: jlapeyre/diffeqpy


text:
using EngEconomics

# Given
xInit = 80000
xMaintenance = 30000
xSalvage = 40000
yInit = 97000
yMaintenance = 27000
ySalvage = 50000
i = 0.15
N = 3 # years

# Find: The Future Worth Analysis
# Determine Alternatives
option1Init = xInit * 2
option2Init = yInit * 2
option3Init = xInit + yInit
option1Maintenance = xMaintenance * 2
option2Maintenance = yMaintenance * 2
option3Maintenance = xMaintenance + yMaintenance
option1Salvage = 2 * xSalvage
option2Salvage = 2 * ySalvage
option3Salvage = xSalvage + ySalvage

# Determine the FW of Option1
option1InitFW = -option1Init * compoundAmountFactor(i, N)
option1MaintenanceFW = -option1Maintenance * uniformSeriesCompoundAmoundFactor(i, N)
option1SalvageFW = option1Salvage
option1FW = option1InitFW + option1MaintenanceFW + option1SalvageFW

# Determine the FW of Option2
option2InitFW = -option2Init * compoundAmountFactor(i, N)
option2MaintenanceFW = -option2Maintenance * uniformSeriesCompoundAmoundFactor(i, N)
option2SalvageFW = option2Salvage
option2FW = option2InitFW + option2MaintenanceFW + option2SalvageFW

# Determine the FW of Option3
option3InitFW = -option3Init * compoundAmountFactor(i, N)
option3MaintenanceFW = -option3Maintenance * uniformSeriesCompoundAmoundFactor(i, N)
option3SalvageFW = option3Salvage
option3FW = option3InitFW + option3MaintenanceFW + option3SalvageFW

# Pick the Better one
maxOption = max(option1FW, option2FW, option3FW)

if maxOption == option1FW
	println("Pick XX")
elseif maxOption == option2FW
	println("Pick YY")
else
	println("Pick XY")
end


text:
type MktEvent
    d::Int32 # date
    t::Int64# time
    dt::DateTime # datetime
    isExec::Int32
    bid::Float32
    ask::Float32
    lastAmt::Int32
    BS::Int32
    OC::Int32
    CP::Int32
    oType::Int32
    AcctID::Int32
    function MktEvent()
        new(zero(Int32),   zero(Int64),   typemin(DateTime), zero(Int32),
            zero(Float32), zero(Float32), zero(Int32),       zero(Int32),
            zero(Int32),   zero(Int32),   zero(Int32),       zero(Int32))
    end
end

function datetimeParse(d::Int32, t::Int64)
    yr = floor(d / 10000)
    mn = floor((d - yr * 10000) / 100)
    dy = floor(d - (yr * 10000 + mn * 100))
    hr = floor(t/ 3600000)
    mi = floor((t - hr * 3600000) / 60000)
    sc = floor((t - (hr * 3600000 + mi * 60000)) /  1000)
    ms = floor(t - (hr * 3600000 + mi * 60000 + sc * 1000))
    DateTime(yr, mn, dy, hr, mi, sc, ms)
end

type IODataCache
    ints::Array{Int32,1}
    intsLrg::Array{Int64,1}
    floats::Array{Float32,1}
    IODataCache() = new(zeros(Int32,1),zeros(Int64,1),zeros(Float32,1))
end

function fillBinaryEvent!(io::IOStream, d::IODataCache, e::MktEvent)
    read!(io, d.ints); e.d = d.ints[1]
    read!(io, d.intsLrg); e.t = d.intsLrg[1]
    e.dt = datetimeParse(e.d,e.t)
    read!(io, d.ints); e.isExec = d.ints[1]
    read!(io, d.floats); e.bid = d.floats[1]
    read!(io, d.floats); e.ask = d.floats[1]
    read!(io, d.ints); e.lastAmt = d.ints[1]
    read!(io, d.ints); e.BS = d.ints[1] * -1
    read!(io, d.ints); e.OC = d.ints[1]
    read!(io, d.ints); e.CP = d.ints[1]
    read!(io, d.ints); e.oType = d.ints[1]
    read!(io, d.ints); e.AcctID = d.ints[1]
end

oTypes = ["UNDEFINED" => -1,"null" => 0, "Bulk Settlement"=> 1,"LIMIT"=> 2,"Loss Cut"=> 3,
          "Market"=> 4,"null"=> 5,"Quoted"=> 6,"STOP"=> 7]

max_stars_repo_name: Andy-P/PriceDist

max_issues_repo_name: Andy-P/PriceDist

max_forks_repo_name: Andy-P/PriceDist


text:
isdefined(Turf, :earth_radius) || const earth_radius = 6371008.8 # Earth Radius used with the Harvesine formula and approximated using a spherical (non-ellipsoid) Earth.

isdefined(Turf, :factors) || const factors = Dict(
    "centimeters" => earth_radius * 100,
    "centimetres" => earth_radius * 100,
    "degrees" => earth_radius / 111325,
    "feet" => earth_radius * 3.28084,
    "inches" => earth_radius * 39.370,
    "kilometers" => earth_radius / 1000,
    "kilometres" => earth_radius / 1000,
    "meters" => earth_radius,
    "metres" => earth_radius,
    "miles" => earth_radius / 1609.344,
    "millimeters" => earth_radius * 1000,
    "millimetres" => earth_radius * 1000,
    "nauticalmiles" => earth_radius / 1852,
    "radians" => 1,
    "yards" => earth_radius / 1.0936
)

isdefined(Turf, :units_factors) || const units_factors = Dict(
    "centimeters" => 100,
    "centimetres" => 100,
    "degrees" => 1 / 111325,
    "feet" => 3.28084,
    "inches" => 39.370,
    "kilometers" => 1 / 1000,
    "kilometres" => 1 / 1000,
    "meters" => 1,
    "metres" => 1,
    "miles" => 1 / 1609.344,
    "millimeters" => 1000,
    "millimetres" => 1000,
    "nauticalmiles" => 1 / 1852,
    "radians" => 1 / earth_radius,
    "yards" => 1 / 1.0936
)

isdefined(Turf, :area_factors) || const area_factors = Dict(
    "acres" => 0.000247105,
    "centimeters" => 10000,
    "centimetres" => 10000,
    "feet" => 10.763910417,
    "inches" => 1550.003100006,
    "kilometers" => 0.000001,
    "kilometres" => 0.000001,
    "meters" => 1,
    "metres" => 1,
    "miles" => 3.86e-7,
    "millimeters" => 1000000,
    "millimetres" => 1000000,
    "yards" => 1.195990046
)


text:
abstract type InterpolationGrid end

Base.minimum(g::InterpolationGrid) = g.min
Base.maximum(g::InterpolationGrid) = g.max
Base.min(g::InterpolationGrid) = g.min
Base.max(g::InterpolationGrid) = g.max
Base.range(g::InterpolationGrid) = g.min:g.step:g.max
Base.step(g::InterpolationGrid) = g.step

export InterpolationGrid


text:
using InfiniteLinearAlgebra, BandedMatrices, LazyArrays, FillArrays, BenchmarkTools, Test
import LazyArrays: MemoryLayout, arguments, resizedata!, partialqr!

@testset "Bessel timings" begin
    A = BandedMatrix(0 => -2*(0:∞)/1_000_000.0, 1 => Ones(∞), -1 => Ones(∞))
    n = 1_000_000
    V = view(A.data.args[2],:,1:n)
    b = similar(V)
    @test @belapsed(copyto!(b,V)) ≤ 0.01
    @test @belapsed(A.data[:,1:n]) ≤ 0.02
    @test @belapsed(A[1:n,1:n]) ≤ 0.02

    C = cache(A); @time resizedata!(C, n+1,n);
    s = 1000; C = cache(A); @time for m = s:s:1_000_000 resizedata!(C, m,m); end;
    s = 1000; F = qr(A); @time for m = s:s:1_000_000 partialqr!(F.factors.data,m); end;
    @time A \ Vcat([1.0],Zeros(∞));
end


text:
module SinkTest

using HPAT

@acc hpat function sink(file3)
    y = randn(10)
    DataSink(y,HDF5,"/y", file3)
end

@acc hpat function sink2(file3)
    y = randn(3,10)
    DataSink(y,HDF5,"/y2", file3)
end

@acc hpat function sink_txt(file3)
    y = randn(10)
    DataSink(y, TXT, file3)
end

@acc hpat function sink_txt2(file3)
    y = randn(3,10)
    DataSink(y, TXT, file3)
end

using MPI

function main()
    SinkTest.sink("y.hdf5")
    SinkTest.sink2("y.hdf5")
    SinkTest.sink_txt("y.txt")
    SinkTest.sink_txt2("y2.txt")
    if MPI.Comm_rank(MPI.COMM_WORLD)==0 rm("y.hdf5") end
end

end

SinkTest.main()


text:
## New ToeplitzOperator routines

ToeplitzOperator(L::UniformScaling) = ToeplitzOperator(eltype(L)[],[L.λ])

issymmetric(T::ToeplitzOperator) =
    length(T.negative)==length(T.nonnegative)-1&&T.negative==T.nonnegative[2:end]

for OP in (:+,:-)
    @eval begin
        $OP(A::ToeplitzOperator,c::UniformScaling) = $OP(A,ToeplitzOperator(c))
        $OP(c::UniformScaling,A::ToeplitzOperator) = $OP(ToeplitzOperator(c),A)
        function $OP(A::ToeplitzOperator,B::ToeplitzOperator)
            n=max(length(A.negative),length(B.negative))
            m=max(length(A.nonnegative),length(B.nonnegative))
            ToeplitzOperator($OP(pad(A.negative,n),pad(B.negative,n)),
                             $OP(pad(A.nonnegative,m),pad(B.nonnegative,m)))
        end
    end
end

## Symetric tridiagonal finite dimensional operators

struct SymTriOperator{T} <: TridiagonalOperator{T}
    dv::Vector{T}
    ev::Vector{T}
end
SymTriOperator(A::Vector,B::Vector) =
    SymTriOperator{promote_type(eltype(A),eltype(B))}(A,B)

for OP in (:domainspace,:rangespace)
    @eval $OP(::SymTriOperator) = ℓ⁰
end

function getindex(S::SymTriOperator,k::Integer,j::Integer)
    if k ≤ length(S.dv) && k == j
        S.dv[k]
    elseif k ≤ length(S.ev) && j==k+1
        S.ev[k]
    elseif 2 ≤ k ≤ length(S.ev)+1 && j==k-1
        S.ev[k-1]
    else
        zero(eltype(S))
    end
end

*(c::Number,A::SymTriOperator) = SymTriOperator(c*A.dv,c*A.ev)
*(A::SymTriOperator,c::Number) = SymTriOperator(c*A.dv,c*A.ev)

function SymTridiagonal(S::SymTriOperator,kr::UnitRange,jr::UnitRange)
    n=last(kr)
    @assert n==last(jr)
    SymTridiagonal(pad(S.dv,n),pad(S.ev,n-1))
end

# Represents a SymTriOperator + Symmetric ToeplitzOperator
struct SymTriPertToeplitz{T} <: TridiagonalOperator{T}
    dv::Vector{T}
    ev::Vector{T}
    a::T
    b::T

    SymTriPertToeplitz{T}(dv::Vector{T},ev::Vector{T},a::T,b::T) where T = new(dv,ev,a,b)
    SymTriPertToeplitz{T}(dv::Vector,ev::Vector,a,b) where T = new(Vector{T}(dv),Vector{T}(ev),T(a),T(b))
end

SymTriPertToeplitz(dv::Vector,ev::Vector,a,b) =
    SymTriPertToeplitz{promote_type(eltype(dv),eltype(dv),typeof(a),typeof(b))}(dv,ev,a,b)

function SymTriPertToeplitz(T::ToeplitzOperator,K::SymTriOperator)
    @assert bandwidths(T) == (1,1) && issym(T)
    SymTriPertToeplitz(K.dv+T.nonnegative[1],K.ev+T.nonnegative[2],T.nonnegative...)
end

function SymTriPertToeplitz(K::SymTriOperator{TT}) where TT
    SymTriPertToeplitz(K.dv,K.ev,zero(TT),zero(TT))
end

function SymTriPertToeplitz(T::ToeplitzOperator)
    @assert issym(T)

    if isdiag(T)
        SymTriPertToeplitz(eltype(T)[],eltype(T)[],T.nonnegative[1],zero(eltype(T)))
    elseif bandwidths(T) == (1,1)
        SymTriPertToeplitz(eltype(T)[],eltype(T)[],T.nonnegative...)
    else
        error("Not a tridiagonal operator")
    end
end

for OP in (:domainspace,:rangespace)
    @eval $OP(::SymTriPertToeplitz) = ℓ⁰
end

function getindex(S::SymTriPertToeplitz, kr::AbstractInfUnitRange, jr::AbstractInfUnitRange)
    k=first(kr)
    @assert k==first(jr)

    SymTriPertToeplitz(S.dv[k:end],S.ev[k:end],S.a,S.b)
end

function getindex(S::SymTriPertToeplitz,k::Integer,j::Integer)
        if 2 ≤ k && j ==k-1
            k ≤ length(S.ev)+1 ? S.ev[k-1] : S.b
        elseif j==k+1
            k ≤ length(S.ev) ? S.ev[k] : S.b
        elseif j==k
            k ≤ length(S.dv) ? S.dv[k] : S.a
        else
            zero(eltype(S))
        end
end

function SymTridiagonal(S::SymTriPertToeplitz,kr::UnitRange,jr::UnitRange)
    n=last(kr)
    @assert n==last(jr)
    dv= n>length(S.dv) ?
        [S.dv;S.a*ones(eltype(S),n-length(S.dv))] :
        S.dv[1:n]

    ev= n-1>length(S.ev) ?
        [S.ev;S.b*ones(eltype(S),n-1-length(S.ev))] :
        S.ev[1:n-1]

    SymTridiagonal(dv,ev)
end

## represents T + K where T is Toeplitz and K is finite-dimensional
struct PertToeplitz{S} <: Operator{S}
    T::ToeplitzOperator{S}
    K::FiniteOperator{BandedMatrix{S,Matrix{S},Base.OneTo{Int}},S}
end

for OP in (:domainspace,:rangespace)
    @eval $OP(::PertToeplitz) = ℓ⁰
end

bandwidths(P::PertToeplitz) =
    max(bandwidth(P.T,1),bandwidth(P.K,1)),max(bandwidth(P.T,2),bandwidth(P.K,2))

getindex(P::PertToeplitz,k::Integer,j::Integer) =
    P.T[k,j]+P.K[k,j]

getindex(P::PertToeplitz,k::AbstractInfUnitRange, j::AbstractInfUnitRange) =
    P.T[k,j]+P.K[k,j]

for OP in (:+,:-)
    @eval begin
        function $OP(A::SymTriPertToeplitz,B::SymTriPertToeplitz)
            n_dv = max(length(A.dv),length(B.dv))
            n_ev = max(length(A.ev),length(B.ev))
            SymTriPertToeplitz($OP([A.dv;fill(A.a,n_dv-length(A.dv))],
                               [B.dv;fill(B.a,n_dv-length(B.dv))]),
                               $OP([A.ev;fill(A.b,n_ev-length(A.ev))],
                                   [B.ev;fill(B.b,n_ev-length(B.ev))]),
                            $OP(A.a,B.a),$OP(A.b,B.b))
        end
        $OP(T::ToeplitzOperator,K::FiniteOperator) = PertToeplitz(T,$OP(K))
        $OP(T::ToeplitzOperator,K::SymTriOperator) = SymTriPertToeplitz(T,$OP(K))
        $OP(S::SymTriPertToeplitz,K::SymTriOperator) = $OP(S,SymTriPertToeplitz(K))
        $OP(K::Union{FiniteOperator,SymTriOperator},T::ToeplitzOperator) = $OP(T)+K

        $OP(A::SymTriPertToeplitz,c::UniformScaling) =
            SymTriPertToeplitz(broadcast($OP,A.dv,c.λ),A.ev,broadcast($OP,A.a,c.λ),A.b)
    end
end

-(A::SymTriPertToeplitz) = SymTriPertToeplitz(-A.dv,-A.ev,-A.a,-A.b)
+(c::UniformScaling,A::SymTriPertToeplitz) = SymTriPertToeplitz(c.λ+A.dv,A.ev,c.λ+A.a,A.b)
-(c::UniformScaling,A::SymTriPertToeplitz) = SymTriPertToeplitz(c.λ-A.dv,-A.ev,c.λ-A.a,-A.b)

*(c::Number,A::SymTriPertToeplitz) = SymTriPertToeplitz(c*A.dv,c*A.ev,c*A.a,c*A.b)
/(A::SymTriPertToeplitz,c::Number) = (1/c)*A

for OP in (:ql,:(eigvals),:(eig))
    @eval $OP(A::ToeplitzOperator)=$OP(SymTriPertToeplitz(A))
end

ql(A::SymTriPertToeplitz) = ql(A.dv,A.ev,A.a,A.b)

function spectralmeasure(J::SymTriPertToeplitz)
    μ = spectralmeasure(.5*(J.dv-J.a)/J.b,.5*J.ev/J.b)
    2*J.b*setdomain(μ,domain(μ) + J.a)
end

function principalresolvent(J::SymTriPertToeplitz)
    r = principalresolvent(.5*(J.dv-J.a)/J.b,.5*J.ev/J.b)
    λ -> (.5/J.b)*r(.5*(λ-J.a)/J.b)
end

eigvals(A::SymTriPertToeplitz) = domain(spectralmeasure(A))
spectrum(A::Operator) = eigvals(A)

function *(L::PertToeplitz,Q::HessenbergUnitary{'L'})
    n=max(size(L.K.matrix,1),length(Q.s)+3)

    if bandwidths(L) == (2,0)
         # We check if L*Q is tridiagin al
        tol=1E-14*(maximum(L.T)+maximum(L.K))
        istri=true
        for k=3:n
            if abs(L[k,k-2]*hc(Q,k-1)+L[k,k-1]*hs(Q,k-2)*hc(Q,k)+L[k,k]*hs(Q,k-2)*hs(Q,k-1)*hc(Q,k+1))>tol
                istri=false
                break
            end
        end
        if istri
            issym=true
            if !isapprox(-L[1,1]*hs(Q,1),L[2,1]*hc(Q,1)*hc(Q,2)+L[2,2]*hc(Q,1)*hc(Q,3)*hs(Q,1);atol=tol)
                issym=false
            end

            if issym
                for k=2:n+1  # kth row
                    if !isapprox(-L[k+1,k-1]*hs(Q,k-1)+L[k+1,k]*hc(Q,k)*hc(Q,k+1)+L[k+1,k+1]*hc(Q,k)*hc(Q,k+2)*hs(Q,k),
                    -L[k,k]*hs(Q,k);atol=tol)
                        issym=false
                        break
                    end
                end
            end

            if issym
               # result is SymTriToeplitxz

                ev=Array{Float64}(undef, max(min(size(L.K.matrix,1),size(L.K.matrix,2)),
                                     length(Q.s)))
                for k=1:length(ev)
                    ev[k]=-L[k,k]*hs(Q,k)
                end

                dv=Array{Float64}(undef, max(length(Q.s)+1,size(L.K.matrix,1)))
                dv[1]=hc(Q,1)*hc(Q,2)*L[1,1]
                for k=2:length(dv)
                    dv[k]=-hs(Q,k-1)*L[k,k-1]+hc(Q,k)*hc(Q,k+1)*L[k,k]
                end

                t1=-L.T[1,1]*Q.s∞
                t0=-Q.s∞*L.T[2,1]+Q.c∞^2*L.T[1,1]

                si=Q.sign ? 1 : -1
                return SymTriPertToeplitz(si*dv,si*ev,si*t0,si*t1)
            end
        end
    end

    # default constructor
    TimesOperator(L,Q)
end


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

# tests related to functional programming functions and styles

# map -- array.jl
@test isequal(map((x)->"$x"[end:end], 9:11), ["9", "0", "1"])
# TODO: @test map!() much more thoroughly
let a = [1.0, 2.0]
    map!(sin, a, a)
    @test isequal(a, sin.([1.0, 2.0]))
end
# map -- ranges.jl
@test isequal(map(sqrt, 1:5), [sqrt(i) for i in 1:5])
@test isequal(map(sqrt, 2:6), [sqrt(i) for i in 2:6])

# map on ranges should evaluate first value only once (#4453)
let io=IOBuffer(3)
    map(x->print(io,x), 1:2)
    @test String(take!(io))=="12"
end

# map over Bottom[] should return Bottom[]
# issue #6719
@test isequal(typeof(map(x -> x, Vector{Union{}}(uninitialized, 0))), Vector{Union{}})

# maps of tuples (formerly in test/core.jl) -- tuple.jl
@test map((x,y)->x+y,(1,2,3),(4,5,6)) == (5,7,9)
@test map((x,y)->x+y,
            (100001,100002,100003),
            (100004,100005,100006)) == (200005,200007,200009)

# maps of strings (character arrays) -- string.jl
@test map((c)->Char(c+1), "abcDEF") == "bcdEFG"

# issue #10633
@test isa(map(Integer, Any[1, 2]), Vector{Int})
@test isa(map(Integer, Any[]), Vector{Integer})

# filter -- array.jl
@test isequal(filter(x->(x>1), [0 1 2 3 2 1 0]), [2, 3, 2])
# TODO: @test_throws isequal(filter(x->x+1, [0 1 2 3 2 1 0]), [2, 3, 2])
@test isequal(filter(x->(x>10), [0 1 2 3 2 1 0]), [])
@test isequal(filter((ss)->length(ss)==3, ["abcd", "efg", "hij", "klmn", "opq"]), ["efg", "hij", "opq"])

# numbers
@test size(collect(1)) == size(1)

@test isa(collect(Any, [1,2]), Vector{Any})

# foreach
let a = []
    foreach(()->push!(a,0))
    @test a == [0]
    a = []
    foreach(x->push!(a,x), [1,5,10])
    @test a == [1,5,10]
    a = []
    foreach((args...)->push!(a,args), [2,4,6], [10,20,30])
    @test a == [(2,10),(4,20),(6,30)]
end

# generators (#4470, #14848)
@test sum(i/2 for i=1:2) == 1.5
@test collect(2i for i=2:5) == [4,6,8,10]
@test collect((i+10j for i=1:2,j=3:4)) == [31 41; 32 42]
@test collect((i+10j for i=1:2,j=3:4,k=1:1)) == reshape([31 41; 32 42], (2,2,1))

let A = collect(Base.Generator(x->2x, Real[1.5,2.5]))
    @test A == [3,5]
    @test isa(A,Vector{Float64})
end

let f(g) = (@test size(g.iter)==(2,3))
    f(i+j for i=1:2, j=3:5)
end

@test collect(Base.Generator(+, [1,2], [10,20])) == [11,22]

# generator ndims #16394
let gens_dims = [((i for i = 1:5),                    1),
                 ((i for i = 1:5, j = 1:5),           2),
                 ((i for i = 1:5, j = 1:5, k = 1:5),  3),
                 ((i for i = Array{Int,0}(uninitialized)),           0),
                 ((i for i = Vector{Int}(uninitialized, 1)),          1),
                 ((i for i = Matrix{Int}(uninitialized, 1, 2)),       2),
                 ((i for i = Array{Int}(uninitialized, 1, 2, 3)),    3)]
    for (gen, dim) in gens_dims
        @test ndims(gen) == ndims(collect(gen)) == dim
    end
end

# generator with destructuring
let d = Dict(:a=>1, :b=>2), a = Dict(3=>4, 5=>6)
    @test Dict( v=>(k,) for (k,v) in d) == Dict(2=>(:b,), 1=>(:a,))
    @test Dict( (x,b)=>(c,y) for (x,c) in d, (b,y) in a ) == Dict((:a,5)=>(1,6),(:b,5)=>(2,6),(:a,3)=>(1,4),(:b,3)=>(2,4))
end

let i = 1
    local g = (i+j for i=2:2, j=3:3)
    @test first(g) == 5
    @test i == 1
end

# generators and guards
let gen = (x for x in 1:10)
    @test gen.iter == 1:10
    @test gen.f(first(1:10)) == next(gen, start(gen))[1]
    for (a,b) in zip(1:10,gen)
        @test a == b
    end
end

let gen = (x * y for x in 1:10, y in 1:10)
    @test collect(gen) == collect(1:10) .* collect(1:10)'
    @test first(gen) == 1
    @test collect(gen)[1:10] == 1:10
end

let gen = Base.Generator(+, 1:10, 1:10, 1:10)
    @test first(gen) == 3
    @test collect(gen) == 3:3:30
end

let gen = (x for x in 1:10 if x % 2 == 0), gen2 = Iterators.filter(x->x % 2 == 0, x for x in 1:10)
    @test collect(gen) == collect(gen2)
    @test collect(gen) == 2:2:10
end

let gen = ((x,y) for x in 1:10, y in 1:10 if x % 2 == 0 && y % 2 == 0),
    gen2 = Iterators.filter(x->x[1] % 2 == 0 && x[2] % 2 == 0, (x,y) for x in 1:10, y in 1:10)
    @test collect(gen) == collect(gen2)
end

# inference on vararg generator of a type (see #22907 comments)
let f(x) = collect(Base.Generator(=>, x, x))
    @test @inferred(f((1,2))) == [1=>1, 2=>2]
end

# generators with nested loops (#4867)
@test [(i,j) for i=1:3 for j=1:i] == [(1,1), (2,1), (2,2), (3,1), (3,2), (3,3)]

@test [(i,j) for i=1:3 for j=1:i if j>1] == [(2,2), (3,2), (3,3)]

# issue #18707
@test [(q,d,n,p) for q = 0:25:100
                 for d = 0:10:100-q
                 for n = 0:5:100-q-d
                 for p = 100-q-d-n
                 if p < n < d < q] == [(50,30,15,5), (50,30,20,0), (50,40,10,0), (75,20,5,0)]

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using Gtk.ShortNames, GtkReactive, Graphics, Colors

win = Window("Drawing")
c = canvas(UserUnit)       # create a canvas with user-specified coordinates
push!(win, c)

const lines = Signal([])   # the list of lines that we'll draw
const newline = Signal([]) # the in-progress line (will be added to list above)

# Add mouse interactions
const drawing = Signal(false)  # this will be true if we're dragging a new line
sigstart = map(c.mouse.buttonpress) do btn
    if btn.button == 1 && btn.modifiers == 0
        push!(drawing, true)   # start extending the line
        push!(newline, [btn.position])
    end
end

const dummybutton = MouseButton{UserUnit}()
sigextend = map(filterwhen(drawing, dummybutton, c.mouse.motion)) do btn
    push!(newline, push!(value(newline), btn.position))
end

sigend = map(c.mouse.buttonrelease) do btn
    if btn.button == 1
        push!(drawing, false)  # stop extending the line
        push!(lines, push!(value(lines), value(newline)))
        push!(newline, [])
    end
end

# Draw on the canvas
redraw = draw(c, lines, newline) do cnvs, lns, newl
    fill!(cnvs, colorant"white")   # background is white
    set_coordinates(cnvs, BoundingBox(0, 1, 0, 1))  # set coords to 0..1 along each axis
    ctx = getgc(cnvs)
    for l in lns
        drawline(ctx, l, colorant"blue")
    end
    drawline(ctx, newl, colorant"red")
end

function drawline(ctx, l, color)
    isempty(l) && return
    p = first(l)
    move_to(ctx, p.x, p.y)
    set_source(ctx, color)
    for i = 2:length(l)
        p = l[i]
        line_to(ctx, p.x, p.y)
    end
    stroke(ctx)
end

Gtk.showall(win)


text:
using CLI
using Base.Test

include("Elements.jl")


text:
__precompile__()
module Holomorphic
    using Base, ApproxFun, SingularIntegralEquations

import ApproxFun: UnivariateSpace, domain, evaluate, ComplexBasis, spacescompatible,
                    Space, defaultFun, union_rule, ConstantSpace, points, checkpoints,
                    transform, plan_transform, conversion_rule, coefficients
import SingularIntegralEquations: stieltjes

export ℂ

# represents the whole Planess
immutable ComplexPlane <: Domain{Complex128,2} end

const ℂ=ComplexPlane()
Base.reverse(C::ComplexPlane)=C

Base.intersect(a::ComplexPlane,b::ComplexPlane)=a
Base.union(a::ComplexPlane,b::ComplexPlane)=a

# represents S\A
immutable Complement{SS,AA,T,d} <: Domain{T,d}
    full::SS
    del::AA
end

Base.intersect(a::Complement,b::Complement)=Complement(a.full ∩ b.full,a.del ∪ b.del)
Base.union(a::Complement,b::Complement)=Complement(a.full ∪ b.full,a.del ∩ b.del)

Complement(A::Domain,B::Domain)=Complement{typeof(A),typeof(B),eltype(A),ndims(A)}(A,B)
Base.setdiff(a::Domain,b::Domain)=Complement(a,b)
Base.reverse(a::Complement)=reverse(a.full)\reverse(a.del)

immutable StieltjesSpace{S,DD} <:  Space{ComplexBasis,Complement{ComplexPlane,DD},2}
    space::S
    function StieltjesSpace(sp::S)
        @assert isa(domain(sp),DD)
        new(sp)
    end
end

StieltjesSpace(sp::UnivariateSpace)=StieltjesSpace{typeof(sp),typeof(domain(sp))}(sp)
spacescompatible(a::StieltjesSpace,b::StieltjesSpace)=spacescompatible(a.space,b.space)

domain(sp::StieltjesSpace)=ℂ\domain(sp.space)
evaluate(v::AbstractVector,sp::StieltjesSpace,z)=stieltjes(sp.space,v,z)
stieltjes(f::Fun)=Fun(f.coefficients,StieltjesSpace(space(f)))

union_rule(A::StieltjesSpace,B::StieltjesSpace)=StieltjesSpace(A.space ∪ B.space)
coefficients(v::Vector,A::StieltjesSpace,B::StieltjesSpace)=coefficients(v,A.space,B.space)

# construct spaces for complement
# we assume boundedness

Space{DD<:Interval}(Γ::Complement{ComplexPlane,DD})=StieltjesSpace(JacobiWeight(0.5,0.5,Ultraspherical{1}(Γ.del)))
Space{DD<:Circle}(Γ::Complement{ComplexPlane,DD})=StieltjesSpace(Laurent(Γ.del))
Space{DD<:UnionDomain}(Γ::Complement{ComplexPlane,DD})=StieltjesSpace(PiecewiseSpace(map(d->Space(d).space,map(x->ComplexPlane()\x,Γ.del))))

## Interval

checkpoints{DD<:Interval}(C::Complement{ComplexPlane,DD})=[fromcanonical(C.del,1im),fromcanonical(C.del,2.)]

joukowsky(z)=(z+1./z)./2

points{DD<:Interval}(S::StieltjesSpace{JacobiWeight{Ultraspherical{1,DD},DD},DD},n)=
    fromcanonical(domain(S).del,joukowsky(points((1-1/n)*Circle(),n)))

plan_transform{DD<:Interval}(S::StieltjesSpace{JacobiWeight{Ultraspherical{1,DD},DD},DD},vals::Vector)=
        plan_transform(Taylor((1-1/length(vals))*Circle()),vals)

function transform{DD<:Interval}(S::StieltjesSpace{JacobiWeight{Ultraspherical{1,DD},DD},DD},vals,plan)
    @assert S.space.α==S.space.β==0.5
    n=length(vals)
    L=Laurent((1-1/n)*Circle())
    cfs1=transform(L,vals)
    for k=3:2:length(cfs1)
        cfs1[k]*=1/((1-1/n)^(div(k+1,2)-1)*π)
    end
    for k=2:2:length(cfs1)
        cfs1[k]*=(1-1/n)^(div(k,2))/(π)
    end

    L=Laurent(1/(1-1/n)*Circle())
    vals2=[vals[1];reverse!(vals[2:end])]
    cfs2=transform(L,vals2)

    for k=3:2:length(cfs2)
        cfs2[k]*=(1-1/n)^(div(k+1,2)-1)/(π)
    end
    for k=2:2:length(cfs2)
        cfs2[k]*=1/((1-1/n)^(div(k,2))*π)
    end
    (cfs1-cfs2)[3:2:end]
end

# function transform{DD<:Interval}(S::StieltjesSpace{JacobiWeight{Ultraspherical{1,DD},DD},DD},vals,plan)
#     @assert S.space.α==S.space.β==0.5
#     n=length(vals)
#     cfs=transform(Taylor((1-1/n)*Circle()),vals,plan)
#     for k=2:length(cfs)
#         cfs[k]*=1/((1-1/n)^(k-1)*π)
#     end
#     cfs
# end

## ℂ\Circle()

# represents strip between a < ℑ(z*exp(-im*θ)) < b
immutable Strip{T} <: Domain{Complex{T},2}
    a::T
    b::T
    θ::T
end

Strip(a,b,θ) = Strip(promote(a,b,θ)...)
Strip(a,b) = Strip(a,b,0)
Strip(b) = Strip(-b,b)

Base.show(io::IO,d::ComplexPlane)=print(io,"ℂ")
Base.show(io::IO,d::Complement)=print(io,"$(d.full) \\$(d.del)")

end #module


text:
transduce_ka(xf::Transducer, op, init, xs, device; groupsize, basesize) =
    transduce_ka(xf'(op), init, xs, device, groupsize, basesize)

function transduce_ka(op, init, xs, device, groupsize, basesize)
    groupsize = max(groupsize, 1)
    basesize = max(basesize, 2)
    xf0, coll = extract_transducer(xs)
    if coll isa Iterators.Zip
        arrays = coll.is
        xf = xf0
    else
        arrays = (coll,)
        xf = opcompose(Map(first), xf0)
    end
    rf = _reducingfunction(xf, op; init = init)
    acc = transduce_impl(device, groupsize, basesize, rf, init, arrays...)
    result = complete(rf, acc)
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
    end
    return result
end

function transduce_impl(device, groupsize, basesize, rf::F, init, arrays...) where {F}
    (ys, buf) = transduce!(nothing, device, groupsize, basesize, rf, init, arrays...)
    # @info "ys, = transduce!(nothing, rf, ...)" ys
    length(ys) == 1 && return @allowscalar ys[1]
    monoid = asmonoid(always_combine(rf))
    rf2 = Map(first)'(monoid)  # TODO: reduce wrapping
    dest = ys
    while true
        ys, = transduce!(buf, device, groupsize, basesize, rf2, InitialValue(monoid), ys)
        # @info "ys, = transduce!(buf, rf2, ...)" ys
        length(ys) == 1 && return @allowscalar ys[1]
        dest, buf = buf, dest
        # reusing buffer; is it useful?
    end
end

function fake_transduce(rf, xs, init)
    acc1 = next(rf, start(rf, init), first(xs))
    for x in xs
        acc1 = next(rf, acc1, x)
    end
    acc2 = acc1
    for x in xs
        acc2 = next(rf, acc2, x)
    end
    ys = [acc1, acc2]
    acc3 = acc2
    for y in ys
        acc3 = _combine(rf, acc3, y)
    end
    return acc3
end

valueof(::Val{x}) where {x} = x
valueof(x) = x

instantiate_kernel(f, device, ::Nothing) = f(device)
instantiate_kernel(f, device, groupsize) = f(device, valueof(groupsize))

Base.@propagate_inbounds getvalues(i) = ()
Base.@propagate_inbounds getvalues(i, a) = (a[i],)
Base.@propagate_inbounds getvalues(i, a, as...) = (a[i], getvalues(i, as...)...)

function transduce!(
    buf,
    device,
    groupsize,
    basesize,
    rf::F,
    init,
    arrays...,
) where {F}
    idx = eachindex(arrays...)
    n = Int(length(idx))  # e.g., `length(UInt64(0):UInt64(1))` is not an `Int`

    acctype = if buf === nothing
        # global _ARGS = (rf, zip(arrays...), init)
        # @show fake_transduce(rf, zip(arrays...), init)
        return_type(fake_transduce, Tuple{Typeof(rf),Typeof(zip(arrays...)),Typeof(init)})
        # Note: the result of `return_type` is not observable by the
        # caller of the API `transduce_impl`
    else
        eltype(buf)
    end
    # @show acctype

    function compute_sizes(n)
        local nbasecases = cld(n, basesize)
        local groupsize′ = min(nbasecases, nextpow(2, groupsize))
        local blocks = cld(nbasecases, groupsize′)
        return (blocks, groupsize′, nbasecases)
    end
    blocks, groupsize′, nbasecases = compute_sizes(n)
    if buf === nothing
        next_blocks, = compute_sizes(blocks)
        dest_buf = arrayfor(device){acctype}(undef, blocks + next_blocks)
        dest = view(dest_buf, 1:blocks)
        buf = view(dest_buf, blocks+1:length(dest_buf))
    else
        dest = view(buf, 1:blocks)
    end
    # @show nbasecases, blocks, basesize, groupsize, groupsize′

    ev = instantiate_kernel(transduce_kernel!, device, groupsize′)(
        dest,
        rf,
        init,
        basesize,
        nbasecases,
        idx,
        arrays...;
        ndrange = nbasecases,
    )
    wait(ev)
    # TODO: use `dependencies`

    return dest, buf
end

@kernel function transduce_kernel!(
    dest::AbstractArray{T},
    rf::F,
    init,
    basesize,
    nbasecases,
    idx,
    arrays...,
) where {F,T}
    # basecase
    n = length(idx)
    ill = @index(Local, Linear)
    igl = @index(Group, Linear)
    offset = ill - 1 + (igl - 1) * groupsize()[1]

    i1 = offset * basesize + 1
    if i1 <= n
        x1 = @inbounds getvalues(idx[i1], arrays...)
    else
        x1 = @inbounds getvalues(idx[1], arrays...) # random value (not used)
    end
    acc = next(rf, start(rf, init), x1)
    for i in offset*basesize+2:min((offset + 1) * basesize, n)
        x = @inbounds getvalues(idx[i], arrays...)
        acc = next(rf, acc, x)
    end

    if Base.issingletontype(T)
        # nothing to reduce (e.g., @floop w/o @reduce)
    else
        # combineblock
        offsetb = (igl - 1) * groupsize()[1]
        bound = max(0, nbasecases - offsetb)

        # shared mem for a complete reduction
        shared = @localmem(T, (2 * groupsize()[1]))
        @inbounds shared[ill] = acc

        m = ill - 1
        t = ill
        s = 1
        c = nextpow(2, groupsize()[1]) >> 1
        while c != 0
            @synchronize
            if t + s <= bound && iseven(m)
                @inbounds shared[t] = _combine(rf, shared[t], shared[t+s])
                m >>= 1
            end
            s <<= 1
            c >>= 1
        end

        if t == 1
            @inbounds dest[igl] = shared[1]
        end
    end
end

function always_combine(rf::F) where {F}
    @inline function op(a, b)
        _combine(rf, a, b)
    end
    return op
end

# Semantically correct but inefficient (eager) handling of `Reduced`.
@inline _combine(rf, a::Reduced, b::Reduced) = a
@inline _combine(rf, a::Reduced, b) = a
@inline _combine(rf::RF, a, b::Reduced) where {RF} = reduced(combine(rf, a, unreduced(b)))
@inline _combine(rf::RF, a, b) where {RF} = combine(rf, a, b)


text:
"""
    ForceDirectedLayout

The fields are, in order:

- `move`, a tuple to specify whether moves on the x and y axes are allowed
- `k`, a tuple (kₐ,kᵣ) giving the strength of attraction and repulsion
- `exponents`, a tuple (a,b,c,d) giving the exponents for the attraction and
  repulsion functions
- `gravity`, the strength of attraction towards the center, set to `0.0` as a
  default
- `δ`, a floating point constant regulating the attractive force of interaction
  strength -- when set to its default value of 0.0, all edges have the same
  attraction
- `degree`, a boolean to specificy whether the nodes repel one another according
  to their degree

The various coefficients are used to decide how strongly nodes will *attract* or
*repel* one another, as a function of their distance Δ. Specifically, the
default is that connected nodes will attract one another proportionally to
(Δᵃ)×(kₐᵇ), with a=2 and b=-1, and all nodes repel one another proportionally to
(Δᶜ)×(kᵣᵈ) with c=-1 and d=2.

The parameterization for the Fruchterman-Rheingold layout is the default one,
particularly if kₐ=kᵣ. The Force Atlas 2 parameters are kₐ=1 (or b=0), kᵣ set to
any value, a=1, c=-1, d=1. Note that in all cases, the gravity is a multiplying
constant of the resulting attraction force, so it will also be sensitive to
these choices. The `FruchtermanRheingold` and `ForceAtlas2` functions will
return a `ForceDirectedLayout` -- as this object is mutable, you can replace the
exponents at any time.

The δ parameter is particularly important for probabilistic networks, as these
tend to have *all* their interactions set to non-zero values. As such, setting a
value of δ=1 means that the interactions only attract as much as they are
probable.
"""
mutable struct ForceDirectedLayout
    move::Tuple{Bool,Bool}
    k::Tuple{Float64,Float64}
    exponents::Tuple{Float64,Float64,Float64,Float64}
    gravity::Float64
    δ::Float64
    degree::Bool
end

"""
    ForceDirectedLayout(ka::Float64, kr::Float64; gravity::Float64=0.75)

TODO
"""
ForceDirectedLayout(ka::Float64, kr::Float64; gravity::Float64=0.75) = ForceDirectedLayout((true,true), (ka,kr), (2.0, -1.0, -1.0, 2.0), gravity, 0.0, true)

"""
    FruchtermanRheingold(k::Float64; gravity::Float64=0.75)

The default `ForceDirectedLayout` uses the Fruchterman-Rheingold parameters -
this function is simply here to make the code more explicity, and to use a
"strict" version where kᵣ=kₐ.
"""
FruchtermanRheingold(k::Float64; gravity::Float64=0.75) = ForceDirectedLayout(k, k; gravity=gravity)

"""
    ForceAtlas2(k::Float64; gravity::Float64=0.75)

In the Force Atlas 2 layout, the attraction is proportional to the distance, and
the repulsion to the inverse of the distance. Note that kₐ in this layout is set
to 1, so kᵣ is the *relative* repulsion.
"""
ForceAtlas2(k::Float64; gravity::Float64=0.75) = ForceDirectedLayout((true, true), (1.0, k), (1.0, 0.0, -1.0, 1.0), gravity, 0.0, true)

"""
    SpringElectric(k::Float64; gravity::Float64=0.75)

In the spring electric layout, attraction is proportional to distance, and
repulsion to the inverse of the distance squared.
"""
SpringElectric(k::Float64; gravity::Float64=0.75) = ForceDirectedLayout((true,true), (k, k), (1.0, 1.0, -2.0, 1.0), gravity, 0.0, true)

"""
Stops the movement of a node position.
"""
function stop!(n::NodePosition)
    n.vx = 0.0
    n.vy = 0.0
end

"""
Repel two nodes
"""
function repel!(LA::T, n1::NodePosition, n2::NodePosition, fr) where {T <: ForceDirectedLayout}
    δx = n1.x - n2.x
    δy = n1.y - n2.y
    Δ = sqrt(δx^2.0+δy^2.0)
    Δ = Δ == 0.0 ? 0.0001 : Δ
    if LA.move[1]
        n1.vx += δx/Δ*fr(Δ)
        n2.vx -= δx/Δ*fr(Δ)
    end
    if LA.move[2]
        n1.vy += δy/Δ*fr(Δ)
        n2.vy -= δy/Δ*fr(Δ)
    end
end

"""
Attract two connected nodes
"""
function attract!(LA::T, n1::NodePosition, n2::NodePosition, fa) where {T <: ForceDirectedLayout}
    δx = n1.x - n2.x
    δy = n1.y - n2.y
    Δ = sqrt(δx^2.0+δy^2.0)
    if !iszero(Δ)
        if LA.move[1]
            n1.vx -= δx/Δ*fa(Δ)
            n2.vx += δx/Δ*fa(Δ)
        end
        if LA.move[2]
            n1.vy -= δy/Δ*fa(Δ)
            n2.vy += δy/Δ*fa(Δ)
        end
    end
end

"""
Update the position of a node
"""
function update!(n::NodePosition)
    Δ = sqrt(n.vx^2.0+n.vy^2.0)
    if !iszero(Δ)
        n.x += n.vx/Δ*min(Δ, 0.01)
        n.y += n.vy/Δ*min(Δ, 0.01)
    end
    stop!(n)
end

"""
    position!(LA::ForceDirectedLayout, L::Dict{K,NodePosition}, N::T) where {T <: EcologicalNetworks.AbstractEcologicalNetwork} where {K}

One iteration of the force-directed layout routine. Because these algorithms can
take some time to converge, it may be useful to stop every 500 iterations to
have a look at the results. Note that to avoid oscillations, the maximum
displacement at any given time is set to 0.01 units.

These layouts tend to have O(N³) complexity, where N is the number of nodes in
the network. This is because repulsion required to do (N×(N-1))/2 visits on
pairs of nodes, and an optimal layout usually requires s×N steps to converge.
With the maximal displacement set to 0.01, we have found that k ≈ 100 gives
acceptable results. This will depend on the complexity of the network, and its
connectance, as well as the degree and edge strengths distributions.
"""
function position!(LA::ForceDirectedLayout, L::Dict{K,NodePosition}, N::T) where {T <: EcologicalNetworks.AbstractEcologicalNetwork} where {K}

    degdistr = degree(N)

    # Exponents and forces - the attraction and repulsion functions are
    # (Δᵃ)×(kₐᵇ) and (Δᶜ)×(kᵣᵈ)
    a,b,c,d = LA.exponents
    ka, kr = LA.k
    fa(x) = (x^a)*(ka^b)
    fr(x) = (x^c)*(kr^d)

    plotcenter = NodePosition(0.0, 0.0, 0.0, 0.0)

    for (i, s1) in enumerate(species(N))
        attract!(LA, L[s1], plotcenter, (x) -> LA.gravity*fa(x))
        for (j, s2) in enumerate(species(N))
            if j > i
                if LA.degree
                    repel!(LA, L[s1], L[s2], (x) -> (degdistr[s1]+1)*(degdistr[s2]+1)*fr(x))
                else
                    repel!(LA, L[s1], L[s2], fr)
                end
            end
        end
    end

    for int in interactions(N)
        # We can do Bool^δ and it returns the Bool, so that's tight
        attract!(LA, L[int.from], L[int.to], (x) -> N[int.from, int.to]^LA.δ*fa(x))
    end

    for s in species(N)
        update!(L[s])
    end

end