text:
#holds the time-indexed parameters and a pre-allocation for the row-indexed versions
struct VolumePartsMLE{
    TΠ<:AbstractVector,
    TA<:AbstractMatrix,
    TV₀<:AbstractVector,
    Tτ<:AbstractVector,
    Tts<:AbstractVector,
    Texpand<:NamedTuple,
    elTΠ<:Real} <: AbstractVolumeParts{TΠ, TA, Tts, Texpand, elTΠ}

  Π::TΠ

  A::TA
  V₀::TV₀
  τ::Tτ

  #not to be accessed directly- indices for the expand function
  ts::Tts
  expand::Texpand

  #inner constructor needed to identify elTΠ
  function VolumePartsMLE(Π::TΠ, A::TA, V₀::TV₀, τ::Tτ, ts::Tts, expand::Texpand) where {
      TΠ, TA, TV₀, Tτ, Tts, Texpand}

    @assert (Π == reduce(vcat, [vec(A), vec(V₀), vec(τ)]))

    return new{TΠ, TA, TV₀, Tτ, Tts, Texpand, eltype(TΠ)}(Π,A,V₀,τ, ts,expand)
  end
end

#pre-allocates for the vpexpand
function VolumePartsMLE(A::TA, V₀::TV₀, τ::Tτ, ts::Tts, expand::Texpand) where {
    TA<:AbstractMatrix,
    TV₀<:AbstractVector,
    Tτ<:AbstractVector,
    Tts<:AbstractVector,
    Texpand<:NamedTuple}

  #create vector versions for the parameter vector
  vecA = vec(A)
  vecV₀ = vec(V₀)
  vecτ = vec(τ)

  #this will hold the parameter data
  Π = reduce(vcat, [vecA, vecV₀, vecτ])
  start=1
  linkedA = less_unsafe_wrap(TA, source=Π, start=start, dims=size(A))
  start+=length(vecA)
  linkedV₀ = less_unsafe_wrap(TV₀, source=Π, start=start, dims=size(V₀))
  start+=length(vecV₀)
  linkedτ = less_unsafe_wrap(Tτ, source=Π, start=start, dims=size(τ))
  start+=length(vecτ)
  @assert start - 1 == length(Π)

  #final checks on the linked dataframes
  @assert sum(linkedA .== A) + sum(linkedV₀ .== V₀) + sum(linkedτ .== τ)== length(Π)

  Θ = VolumePartsMLE(Π, linkedA, linkedV₀, linkedτ, ts, expand)

  return Θ
end

#basic constructor from dimension arguments
function VolumePartsMLE(T::Int, K::Int, ts::Vector{Int}, ::Type{Tvp}=PARAM[:fluxtype];
    expander::Function=volumepartsexpanders[PARAM[:fluxmodel]]) where {Tvp<:Real}

  #note- includes a work-around until Flux.gpu works
  A = abs.(randn(Tvp, K,T))
  V₀ = abs.(randn(Tvp, T-1))
  τ = abs.(randn(Tvp, T-1))

  return VolumePartsMLE(A, V₀, τ, ts, expander(A, V₀, ts))
end

#creates new VolumeParts with copies of the trainable portion and pass-by-reference
# of the non-trainable indices
function copytrainable(Θ::VolumePartsMLE)
  #copy trainable portions by value
  A  = deepcopy(Θ.A)
  V₀ = deepcopy(Θ.V₀)
  τ = deepcopy(Θ.τ)

  #copy the non-trainable portions by reference
  return VolumePartsMLE(A, V₀, τ, Θ.ts, Θ.expand)
end

#needed due to the segemented array structure
function Base.deepcopy(Θ::VolumePartsMLE)
  A  = deepcopy(Θ.A)
  V₀  = deepcopy(Θ.V₀)
  τ  = deepcopy(Θ.τ)
  ts = deepcopy(Θ.ts)
  expand = deepcopy(Θ.expand)

  return VolumePartsMLE(A,V₀,τ,ts,expand)
end

Flux.@functor VolumePartsMLE
Flux.params(Θ::VolumePartsMLE) = Flux.params(Θ.A, Θ.V₀,Θ.τ,)

function testvolumepartsindicesmle()
  Θ = VolumePartsMLE(200,5, rand(2:200, 10^5))
  Θcopy = deepcopy(Θ)

  #first test getidx
  (Θ.A == reshape(Θ[1:1000], :, 200)) || error("inconsistent getidx dimensions Θ.A")
  (Θ.V₀ == Θ[1001:1199]) || error("inconsistent getidx dimensions Θ.V₀")
  (Θ.τ == Θ[1200:1398]) || error("inconsistent getidx dimensions Θ.τ")

  #now test setidx
  testvals = rand(PARAM[:fluxtype], 1398)
  testidx = Random.randperm(1398)

  @assert allunique(testidx)
  Θ[testidx] .= testvals

  testcombined = hcat(testidx, testvals)
  testcombined = testcombined[sortperm(testcombined[:,1]),:]
  orderedtestvals = vec(testcombined[:,2])
  display(zip(Θ.Π, orderedtestvals))
  @assert sum(orderedtestvals .== Θ.Π) == length(Θ.Π)

  vec(Θcopy.A) .= orderedtestvals[1:1000]
  vec(Θcopy.V₀) .= orderedtestvals[1001:1199]
  vec(Θcopy.τ) .= orderedtestvals[1200:1398]

  (Θcopy.A == Θ.A) || error("inconsistent setidx! dimensions Θ.A")
  (Θcopy.V₀ == Θ.V₀) || error("inconsistent setidx! dimensions Θ.V₀")
  (Θcopy.τ == Θ.τ) || error("inconsistent setidx! dimensions Θ.τ")

  @info "test vol parts indices completed successfully"
end

#takes a volume part object and writes the results to a df

function formresults(panel::AbstractDataFrame, ms::MeasureSpec, Θ::VolumePartsMLE,
    grossexposures::AbstractMatrix)
  Fcoefs::Vector{Symbol} = Fcoefficients(ms)

  tidx::Dict = gettidx(panel, ms)
  tidxinv::Vector{Date} = sort(collect(keys(tidx)))

  if (length(Fcoefs) ≠ size(Θ.A,1) + 1) || Fcoefs[end] ≠ ms.Fintercept
    error("Coeffcients not in expected form. Fcoefs: $Fcoefs")
  end

  results::DataFrame = DataFrame(date=tidxinv)

  #delog the appropriate series
  V₀ = (exp).(Θ.V₀)
  τ = (exp).(Θ.τ)

  results = hcat(results, DataFrame(hcat(grossexposures, [missing; V₀...;], [missing; τ...;]), [Fcoefs; :tau]))

  return results

end

###########∇VolumeParts
#this struct/function computes the Jacobian
vec∇(Θ::VolumePartsMLE, ∇Θ) = reduce(vcat, [vec(∇Θ[Θ.A]),vec(∇Θ[Θ.V₀]),vec(∇Θ[Θ.τ])])

function test∇VolumeParts!(loss, X::NamedTuple, y, Θ::VolumePartsMLE)
  @info "test ∇VolumeParts!"
  Π = Flux.params(Θ) #formatted parameters
  @info "length: $(length(Π))"
  gs = gradient(()->loss(X,y),Π)

  rms(v) = sqrt(mean(v.^2))
  vrms(m) = rms(vec(m))

  @info "test gradrms: Θ.A-$(vrms(gs[Θ.A])); Θ.V₀-$(vrms(gs[Θ.V₀])); Θ.τ-$(vrms(gs[Θ.τ]));"

  ∇Θ! = ∇VolumeParts!(loss, X, y, Θ)

  #out .= reduce(vcat, [∇Θ[Θ.A], vec(∇Θ[Θ.V₀])])#vec∇(Θ, ∇Θ)
  out = similar(Θ.Π)
  ∇Θ!(out, Θ)
  NvecA::Int = length(vec(Θ.A))
  NvecV₀::Int = length(vec(Θ.V₀))
  Nvecτ::Int = length(vec(Θ.τ))
  @info "test ∇Θ!: Θ.A-$(rms(out[1:NvecA]));" *
    " Θ.V₀-$(vrms(out[(NvecA+1):(NvecA + NvecV₀)]));" *
    " Θ.τ-$(vrms(out[(NvecA+NvecV₀+1):(NvecA+NvecV₀+Nvecτ)]));"
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module openQC

#Package Settings
pkg_dir = Pkg.dir("openQC")

#consts
include("utils/consts.jl")
include("utils/math.jl")
include("utils/qusat.jl")

####################
# type definitions #
####################
# quantum computation abstracts
include("base.jl")
# state types
include("states/states.jl")

##############
# Simulators #
##############

# AQC model
include("qc/adiabatic/adiabatic.jl")

end # module

ext: jl

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using Test

@testset "Longest paths" begin
  # Create a directed path graph.
  g = path_digraph(3)
  costs = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0)
  i = ElementaryPathInstance(g, costs, 1, 3)
  d = lp_dp(i)

  @test d.instance == i
  @test d.path == [Edge(1, 2), Edge(2, 3)]

  add_edge!(g, 1, 3)
  costs[Edge(1, 3)] = 3.0
  i = ElementaryPathInstance(g, costs, 1, 3)
  d = lp_dp(i)

  @test d.instance == i
  @test d.path == [Edge(1, 3)]

  # Positive-cost cycle.
  add_edge!(g, 2, 1)
  costs[Edge(2, 1)] = 25.0
  i = ElementaryPathInstance(g, costs, 1, 3)
  @test_logs (:warn, "The graph contains a positive-cost cycle around edge 2 -> 1.") lp_dp(i)
end

@testset "Budgeted longest paths" begin
  @testset "Basic" begin
    # Create a directed path graph.
    g = path_digraph(3)
    rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0)
    weights = Dict(Edge(1, 2) => 1, Edge(2, 3) => 1)
    i = BudgetedLongestPathInstance(g, rewards, weights, 1, 3, budget=4, max_weight=1)
    d = budgeted_lp_dp(i)

    @test d.instance == i
    @test d.path == [] # No path with a total weight of at least 4.

    for β in [0, 1, 2]
      @test d.solutions[1, β] == []
      @test d.solutions[2, β] == ((β == 2) ? [] : [Edge(1, 2)])
      @test d.solutions[3, β] == [Edge(1, 2), Edge(2, 3)]
    end
    for β in [3, 4]
      for v in [1, 2, 3]
        @test d.solutions[v, β] == []
      end
    end

    add_edge!(g, 1, 3)
    rewards[Edge(1, 3)] = 1.0
    weights[Edge(1, 3)] = 4
    i = BudgetedLongestPathInstance(g, rewards, weights, 1, 3, budget=4, max_weight=1)
    d = budgeted_lp_dp(i)

    @test d.instance == i
    @test d.path == [Edge(1, 3)]
  end

  @testset "Conformity" begin
    # More advanced tests to ensure the algorithm works as expected.
    g = complete_digraph(3)
    rewards = Dict(Edge(1, 2) => 1.0, Edge(3, 1) => -1.0, Edge(3, 2) => -1.0, Edge(2, 3) => 1.0, Edge(2, 1) => -1.0, Edge(1, 3) => 0.0)
    weights = Dict(Edge(1, 2) => 0, Edge(3, 1) => 0, Edge(3, 2) => 0, Edge(2, 3) => 0, Edge(2, 1) => 0, Edge(1, 3) => 2)
    i = BudgetedLongestPathInstance(g, rewards, weights, 1, 3, budget=2, max_weight=1)
    d = budgeted_lp_dp(i)

    @test d.path == [Edge(1, 3)]

    @test d.solutions[1, 0] == []
    @test d.solutions[2, 0] == [Edge(1, 2)]
    @test d.solutions[3, 0] == [Edge(1, 2), Edge(2, 3)]
    @test d.states[1, 0] == 0.0
    @test d.states[2, 0] == 1.0
    @test d.states[3, 0] == 2.0

    @test d.solutions[1, 1] == []
    @test d.solutions[2, 1] == [Edge(1, 3), Edge(3, 2)]
    @test d.solutions[3, 1] == [Edge(1, 3)]
    @test d.states[1, 1] == 0.0
    @test d.states[2, 1] == -1.0
    @test d.states[3, 1] == 0.0

    @test d.solutions[1, 2] == []
    @test d.solutions[2, 2] == [Edge(1, 3), Edge(3, 2)]
    @test d.solutions[3, 2] == [Edge(1, 3)]
    @test d.states[1, 2] == 0.0
    @test d.states[2, 2] == -1.0
    @test d.states[3, 2] == 0.0
  end
end


text:
using SpecialFunctions
using QuadGK

function ourquadgk(f, a, b)
    (value, _) = quadgk(f, a, b, )
    return value
end

"""
K(ρ, β, x, r) computes the integral original function, gives a more convinent way to calculste the cumulative sum of the inner function.
"""

"""
wright(ρ, β, z) computes the Wright function of z.

While in the paper, q₁ is a fixed small number which 0<q₁<1 and q₂ is a fixed big number,
Here we set q₁ = 0.5 and q₂ = 100
"""
q₁ = 2
q₂ = 100

function wright(ρ, β, z)
    ϵ = 0.2
    x=abs(z)
    if abs(z) ≤ q₁ #HACK: This part is working fine
        if ρ > 0
            k₀ = max((1-β)/ρ, log(ϵ)/log(abs(z))-1)
        elseif -1 < ρ < 0
            med = (1-β)/(1+ρ)
            k₀ = max(-β/ρ, med, log(π*(1+ρ)*ϵ*abs(z)^(-med))/(log(abs(z))-1-ρ)+med)
        end
        s = zero(z)
        k₀ = ceil(k₀)
        for k=0:k₀
            s += (z^k)/(factorial(k)*gamma(β + ρ*k))#+μz #NOTE: factorial might led to overflow when k is bigger than 21
        end
        return s

    elseif q₁ < abs(z) < q₂ #The most important in applications

        #TODO: What is x stand for?? The real part of z or just a defination??
        med1 = 1-z*cos(ρ*π)
        med2 = 1/ρ+1-β/ρ
        r₀=0
        if (-1 < ρ < -1/2) && (1+ρ ≤ β)
            r₀ = max(med1^(1/(1+ρ)), (-log(-π*ρ*ϵ))^(-1/ρ))
        elseif (-1 < ρ < -1/2) && (β<1+ρ)
            B = -1/(ρ*π)*(abs(med2)+2)*(-2*med2)^(-med2)
            r₀ = max((med1), (abs(med2)+1)^(-1/ρ), (-2*log(ϵ/B))^(-1/ρ))
        elseif (-1/2 <= ρ < 0) & (β ≥ 0)
            r₀ = -log(ϵ*π)
        elseif (-1 < ρ < -1/2) & (β<0)
            r₀ = max(abs(β)+1, -2*log(ϵ*π/(4(abs(β)+2)*β^2)))
        elseif (ρ > 0) & (β ≥ 0)
            r₀ = max(abs(x*cos(ρ*π))^(1/ρ), -log(ϵ*π)+1)
        elseif (ρ > 0) & (β < 0)
            r₀ = max(abs(x*cos(ϵ*π))^(1/ρ), abs(β)+1, -2*log(ϵ*π/(4*(abs(β)+2)*β^2))+2)
        end
        function K(r)
            return exp(-r-x*r^(-ρ)*cos(ρ*π))*r^(-β)*sin(-x*r^(-ρ)*sin(ρ*π)+π*β)
        end
        return ourquadgk(K, 1, r₀) # The improper integral require the lower limit ∈ {0, 1} and r₀≥1, we set 0 as the lower limit.

    elseif abs(z)>q₂
    end

    #FIXME: The third part abs(z) > q₂ is not implement yet

end


text:
"""
    RandomWalks

A module that supports random walks on both integer lattices (of arbitrary finite dimension)
and on graphs.
"""
module RandomWalks
using LinearAlgebra: Diagonal
using LightGraphs: AbstractGraph, adjacency_matrix, degree, grid, nv, vertices
using Distributions: rand, Sampleable
using StochasticProcesses
using StochasticProcesses.MarkovChains
using StochasticProcesses.Samplers

export

# concrete types
MultiRandomWalk, RandomWalk,

# methods
nwalkers, rand, randomwalkmc

# include scripts
include("randomwalks.jl")
include("multirandomwalks.jl")
include("simulaterandomwalks.jl")
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
using Dates

@testset "Ur.jl" begin
    filename = "data\\test.ur"
    dc = DDR2import.Ur.read(filename)
    @test Dates.day(dc["EB"].start_date) == 1
    @test Dates.month(dc["ED"].end_date) == 2
    @test dc["EE"].unitrate == 31.51
    @test dc["EG"].valutaconversion == 0.848504
    @test dc["EG"].valuta == "GBP"
    @test dc["EH"].country == "Netherlands"
    #LAST
    @test dc["UG"].unitrate == 28.30
end


text:
using PartialSvdStoch, LinearAlgebra, Test

# reduceToRank comparison with exact svd
# we check that the reconstructed matrix has nearly the same singular values
function testifastsvd_S(mat, rank)
    @info "doing exact svd ..."
    @time Uexact,Sexact,Vtexact = svd(mat)
    #
    #
    svdpb = LowRankApproxMc(mat)
    @time reduceToRank(svdpb, rank, nbiter = 2)
    U = PartialSvdStoch.getindex(svdpb,:U)
    S = PartialSvdStoch.getindex(svdpb,:S)
    Vt = PartialSvdStoch.getindex(svdpb,:Vt)
    rankReached = svdpb.currentRank
    matApprox = U.value * diagm(S.value) * Vt.value
    @info "relative residual error LowRankApproxMc " norm(matApprox - mat)/norm(mat)
    # doing exact svd of matApprox
    @time Urecons,Srecons,Vtrecons = svd(matApprox)
    deltaRecons  = norm(Srecons[1:rankReached]-Sexact[1:rankReached])
    relerr = deltaRecons/norm(Sexact[1:rankReached])
    @info "relative error on eigen values after reconstruction : " relerr
    relerr < 0.1 ? true : false
end

# reduceToEpsil comparison with exact svd
function testifastsvd_epsil(mat)
    @debug "doing exact svd ..."
    @time Uexact,Sexact,Vtexact = svd(mat)
    #
    @info "reducing to epsil " 0.1
    svdpb = LowRankApproxMc(mat)
    @time reduceToEpsil(svdpb, 50, 0.01)
    U = PartialSvdStoch.getindex(svdpb,:U)
    S = PartialSvdStoch.getindex(svdpb,:S)
    Vt = PartialSvdStoch.getindex(svdpb,:Vt)
    rankReached = svdpb.currentRank
    deltaS  = norm(S.value[1:rankReached]-Sexact[1:rankReached])/norm(Sexact[1:rankReached])
    matApprox = U.value * diagm(S.value) * Vt.value
    relerrMat = norm(matApprox - mat)/norm(mat)
    @info "delta mat after LowRankApproxMc reduce" deltaS  relerrMat
    relerrMat < 0.05 ? true : false
end

# ifsvd and epsil iterations
function testifastsvd_epsil(mat,  samplingSize, epsil)
    #
    @debug "doing LowRankApproxMc reduceToRank"
    svdpb = LowRankApproxMc(mat)
    @time reduceToEpsil(svdpb, samplingSize , epsil)
    U = PartialSvdStoch.getindex(svdpb,:U)
    S = PartialSvdStoch.getindex(svdpb,:S)
    Vt = PartialSvdStoch.getindex(svdpb,:Vt)
    rankReached = svdpb.currentRank
    matApprox = U.value * BLAS.gemm('T','N', U.value, mat)
    relerrMat = norm(matApprox - mat)/norm(mat)
    @info "delta mat after LowRankApproxMc reduce"  norm(matApprox - mat) norm(mat)
    @info "relative residual error" norm(matApprox - mat)/norm(mat)
    relerrMat < 0.05 ? true : false
end

# ifsvd and rank control
function testifastsvd_rank(mat,  k)
    #
    @info "\n doing LowRankApproxMc reduceToRank"
    svdpb = LowRankApproxMc(mat)
    @time reduceToRank(svdpb, k)
    U = PartialSvdStoch.getindex(svdpb,:U)
    S = PartialSvdStoch.getindex(svdpb,:S)
    Vt = PartialSvdStoch.getindex(svdpb,:Vt)
    rankReached = svdpb.currentRank
    matApprox = U.value * BLAS.gemm('T','N', U.value, mat)
    relerrMat = norm(matApprox - mat)/norm(mat)
    @info "delta mat after LowRankApproxMc reduce"  norm(matApprox - mat) norm(mat)
    @info "relative residual error" norm(matApprox - mat)/norm(mat)
    relerrMat < 0.05 ? true : false
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using SimilaritySearch

import SimilaritySearch:
    search, fit, push!

export Kvp, near_and_far, fit, search, push!

mutable struct Kvp{T} <: Index
    db::Vector{T}
    refs::Vector{T}
    sparsetable::Vector{Vector{Item}}
    k::Int
end

function fit(::Type{Kvp}, dist::Function, db::Vector{T}, k::Int, refs::Vector{T}) where T
    @info "Kvp, refs=$(typeof(db)), k=$(k), numrefs=$(length(refs)), dist=$(dist)"
    sparsetable = Vector{Item}[]

    for i in 1:length(db)
        if (i % 10000) == 0
            println(stderr, "advance $(i)/$(length(db))")
        end

        row = near_and_far(dist, db[i], refs, k)
        push!(sparsetable, row)
    end

    Kvp(db, refs, sparsetable, k)
end

function fit(::Type{Kvp}, dist::Function, db::Vector{T}, k::Int, numrefs::Int) where T
    refList = rand(1:length(db), numrefs)
    refs = [db[x] for x in refList]
    fit(Kvp, dist, db, k, refs)
end

function near_and_far(dist::Function, obj::T, refs::Vector{T}, k::Int) where T
    near = KnnResult(k)
    far = KnnResult(k)
    for (refID, ref) in enumerate(refs)
        d = dist(obj, ref)
        push!(near, refID, d)
        push!(far, refID, -d)
    end

    row = Vector{Item}(undef, k+k)
    for (j, item) in enumerate(near)
        row[j] = item
    end

    for (j, item) in enumerate(far)
        row[k+j] = Item(item.objID, -item.dist)
    end

    row
end

function search(index::Kvp{T}, dist::Function, q::T, res::KnnResult) where T
    d::Float64 = 0.0
    qI = [dist(q, piv) for piv in index.refs]

    for i in 1:length(index.db)
        obj::T = index.db[i]
        objSparseRow = index.sparsetable[i]

        discarded::Bool = false
        @inbounds for item in objSparseRow
            pivID = item.objID
            dop = item.dist
            if abs(dop - qI[pivID]) > covrad(res)
                discarded = true
                break
            end
        end

        if discarded
            continue
        end
        d = dist(q, obj)
        push!(res, i, d)
    end

    res
end

function push!(index::Kvp{T}, dist::Function, obj::T) where T
    push!(index.db, obj)
    row = near_and_far(dist, obj, index.refs, index.k)
    push!(index.sparsetable, row)
    length(index.db)
end


text:
"""
    MPI.install_mpiexecjl(; command::String = "mpiexecjl",
                          destdir::String = joinpath(DEPOT_PATH[1], "bin"),
                          force::Bool = false, verbose::Bool = true)

Install the `mpiexec` wrapper to `destdir` directory, with filename `command`.
Set `force` to `true` to overwrite an existing destination file with the same
path.  If `verbose` is `true`, the installation prints information about the
progress of the process.
"""
function install_mpiexecjl(; command::String = "mpiexecjl",
                           destdir::String = joinpath(DEPOT_PATH[1], "bin"),
                           force::Bool = false, verbose::Bool = true)
    # Adapted from https://github.com/fredrikekre/jlpkg.
    destdir = abspath(expanduser(destdir))
    exec = joinpath(destdir, command)
    if ispath(exec) && !force
        error("file `$(exec)` already exists; use `MPI.install_mpiexecjl(force=true)` to overwrite.")
    end
    mkpath(destdir)
    verbose && @info "Installing `$(command)` to `$(destdir)`..."
    cp(joinpath(@__DIR__, "..", "bin", "mpiexecjl"), exec; force = force)
    verbose && @info "Done!"
end


text:
istransparent(::Any) = false
istransparent(::MLJType) = true

"""
    params(m::MLJType)

Recursively convert any transparent object `m` into a named tuple,
keyed on the fields of `m`. An object is *transparent* if
`MLJBase.istransparent(m) == true`. The named tuple is possibly nested
because `params` is recursively applied to the field values, which
themselves might be transparent.

Most objects of type `MLJType` are transparent.

    julia> params(EnsembleModel(atom=ConstantClassifier()))
    (atom = (target_type = Bool,),
     weights = Float64[],
     bagging_fraction = 0.8,
     rng_seed = 0,
     n = 100,
     parallel = true,)

"""
params(m) = params(m, Val(istransparent(m)))
params(m, ::Val{false}) = m
function params(m, ::Val{true})
    fields = fieldnames(typeof(m))
    NamedTuple{fields}(Tuple([params(getfield(m, field)) for field in fields]))
end


text:
# "tests" are a sanity check that the code is runnable in each notebook

path = joinpath(@__DIR__(), "..", "notebooks")

for file in filter(x -> endswith(x, ".ipynb"), readdir(path))
    @info "including notebook $file"
    TempMod = Main.eval(:(module __Temp__; using NBInclude; end))
    ex = :(@nbinclude(joinpath($path, $file)))
    @eval(TempMod, $ex)
end


text:
load("../../crossover",                                 dirname(@__FILE__))
load("../../../representation/permutation/permutation", dirname(@__FILE__))

type OrderedCrossover <: Crossover
  rep::PermutationRepresentation
  rate::Float
  length::Int
  atom::Type

  OrderedCrossover(rp::PermutationRepresentation, rt::Float) =
    new(rp, rt, rp.size, atom(rp))
end

num_inputs(o::OrderedCrossover) = 2
num_outputs(o::OrderedCrossover) = 2

function operate!{T}(o::OrderedCrossover,
  inputs::Vector{IndividualStage{Vector{T}}})

  # Enforce the crossover rate.
  rand() > o.rate && return inputs

  # Buffer the parent chromosomes.
  p1 = get(inputs[1])
  p2 = get(inputs[2])

  # Create the child chromosomes.
  c1 = Array(o.atom, o.length)
  c2 = Array(o.atom, o.length)

  # Randomly select two cut points from the parent chromosomes.
  cut_from = rand(1:o.length)
  cut_to = rand(cut_from:o.length)

  # Calculate contents of each strip.
  s1 = p1[cut_from:cut_to]
  s2 = p2[cut_from:cut_to]

  # Fill in genes for both children, in sequence.
  for i in 1:o.length
    if i >= cut_from && i <= cut_to
      c1[i] = p1[i]
      c2[i] = p2[i]
    else
      c1[i] = in(p2[i], s1) ? p1[i] : p2[i]
      c2[i] = in(p1[i], s2) ? p2[i] : p1[i]
    end
  end

  # Swap the parents with their children.
  set(inputs[1], c1)
  set(inputs[2], c2)
  return inputs
end

register(joinpath(dirname(@__FILE__), "manifest.yml"))


text:
module ProfileAligner

  export AlignmentMatrix, Profile,
         score, align, getstrings, scoreprofiles,
         measurequality, setScoringMatrix

  import DataReader.FastaRecord, DataReader.ScoreMatrix

  scoringMatrix = ScoreMatrix()

  setScoringMatrix(sm :: ScoreMatrix) = scoringMatrix = sm

  getScoringMatrixValue(matrix :: ScoreMatrix, i1 ::Int64, i2::Int64) = i1 > i2 ? matrix.hsh[i1][i2] : matrix.hsh[i2][i1]

  function getScoringMatrixValue(aa1 :: Char, aa2 :: Char)
    haskey(scoringMatrix.keys, aa1) && haskey(scoringMatrix.keys, aa2) &&
      return getScoringMatrixValue(scoringMatrix, scoringMatrix.keys[aa1], scoringMatrix.keys[aa2])
    aa1 == aa2 ? 1 : -1
  end

  debugprint(str :: Any) = 0
  #debugprintln(str :: Any) = println(str)

  immutable Profile{T}
    rawdata :: Array{Char, 2}
    data :: Array{Dict{Char, T}, 1}
    stringsize :: Int
    numberofstrings :: Int
    descriptions :: Array{ASCIIString, 1}

    #Profile(raw :: Array{Char, 2}, desc :: Array{ASCIIString, 1}) = getprofile{T}(raw, desc)
    #Profile(raw :: Array{Char, 2}, desc :: Array{String, 1}) = getprofile(raw, desc)
    Profile(str :: ASCIIString, desc :: ASCIIString = "") = Profile{T}( reshape([ letter for letter in str ], length(str), 1), [desc] )
    Profile(record :: FastaRecord) = Profile{T}(record.sequence, record.description)

    Profile(raw :: Array{Char, 2},
            d :: Array{Dict{Char, T}, 1},
            s1 :: Int,
            s2 :: Int,
            desc :: Array{ASCIIString, 1}) = new{T}(raw, d, s1, s2, desc)
    #function getprofile(raw :: Array{Char, 2}, descriptions :: Array{ASCIIString, 1})

    #function call{T}(::Type{SummedArray}, a::Vector{T})
    function  call{T}(::Type{Profile{T}}, raw :: Array{Char, 2}, descriptions :: Array{ASCIIString, 1} = [])
      rawdata = copy(raw)
      size_1 = size(raw, 1)
      size_2 = size(raw, 2)
      data = Array(Dict{Char, T}, size_1)
      for row in 1 : size_1
        data[row] = Dict{Char, T}()
        s = zero(T)
        for col in 1 : size_2
          letter = raw[row, col]
          data[row][letter] = get!(data[row], letter, zero(T)) + one(T)
          s += one(T)
        end
        if s > 0
          for (key, value) in data[row]
            data[row][key] /= s
          end
        end
      end
      Profile{T}(rawdata, data, size_1, size_2, descriptions)
    end
  end

  type AlignmentMatrix{T}
    matrix :: Array{T, 2}
    path :: Array{Char, 2}

    AlignmentMatrix(m :: Array{T, 2}, p :: Array{Char, 2}) = new(m, p)

    function call{T}(::Type{AlignmentMatrix{T}}, n :: Int64, m :: Int64)
      matrix = zero(Array(T, n, m))
      path = Array(Char, n, m)
      AlignmentMatrix{T}(matrix, path)
    end
  end

  #method returns pair of previous coordinates
  function getprev(direction :: Char, i :: Int64, j :: Int64)
    i < 2 && direction != 'R' && error("1st index is less than 2")
    j < 2 && direction != 'D' && error("2nd index is less than 2")
    direction == 'D' && return (i - 1, j)
    direction == 'R' && return (i, j - 1)
    direction == 'M' && return (i - 1, j - 1)
    error("unknown direction")
  end

  get{T}(profile::Profile{T}, i::Int64, aa :: Char) = Base.get(profile.data[i], aa, zero(T))

  aminoacids = "ARNDCQEGHILKMFPSTWYVBZX"

  get_ri{T}(profile::Profile{T}, i::Int64) = length(keys(profile.data[i]))

  get_nij{T}(profile::Profile{T}, i::Int64, aa::Char) = Base.get(profile.data[i], aa, 0.0)

  # 1. define weighting schemes
  function sequenceWeight1{T}(profile::Profile{T}, sequence::ASCIIString)
    sum([
      1.0/(get_ri(profile, i) * get_nij(profile, i, sequence[i]))
      for i in 1:length(sequence)
    ])
  end

  # 2. define scores
  function score2{T}(P::Profile{T}, Q::Profile{T}, i, j)
    sum([ get(P, i, aa) * get(Q, j, aa) for aa in aminoacids ] )
  end

  function score{T}(P::Profile{T}, Q::Profile{T}, i, j)
    sum([
          get(P, i, aa1) * get(Q, j, aa2) * getScoringMatrixValue(aa1, aa2)
            for aa1 in aminoacids, aa2 in aminoacids
        ] )
  end

  GAP_COST = -1.0

  function prepareAlignmentMatrix{T}(
                P :: Profile{T},
                Q :: Profile{T},
                matrixdata :: AlignmentMatrix{T}
                )
    n = P.stringsize
    m = Q.stringsize
    matrixdata.matrix[1, 1] = zero(T)
    matrixdata.path[1, 1] = 'U'
    for i = 1 : n
      matrixdata.matrix[i + 1, 1] = matrixdata.matrix[i, 1] + GAP_COST
      matrixdata.path[i + 1, 1] = 'D'
    end
    for j = 1 : m
      matrixdata.matrix[1, j + 1] = matrixdata.matrix[1, j] + GAP_COST
      matrixdata.path[1, j + 1] = 'R'
    end
  end

  function buildAlignmentMatrix{T}(
                  P::Profile{T},
                  Q :: Profile{T},
                  matrixdata :: AlignmentMatrix{T},
                  scoreFunc
                  )
    n = P.stringsize
    m = Q.stringsize
    for j = 1 : m
      for i = 1 : n
        # 1. Match
        matrixdata.matrix[i + 1, j + 1] = matrixdata.matrix[i, j] + scoreFunc(P, Q, i, j)
        matrixdata.path[i + 1, j + 1] = 'M'
        # 2. Down
        temp_value = matrixdata.matrix[i, j + 1] + GAP_COST
        if (matrixdata.matrix[i + 1, j + 1] < temp_value)
          matrixdata.matrix[i + 1, j + 1] = temp_value
          matrixdata.path[i + 1, j + 1] = 'D'
        end
        # 3. Right
        temp_value = matrixdata.matrix[i + 1, j] + GAP_COST
        if (matrixdata.matrix[i + 1, j + 1] < temp_value)
          matrixdata.matrix[i + 1, j + 1] = temp_value
          matrixdata.path[i + 1, j + 1] = 'R'
        end

      end
    end
  end

  #returns score from precomputed matrix and last column/row coordinates
  function getscore{T}(P::Profile{T}, Q :: Profile{T}, matrixdata :: AlignmentMatrix{T})
    n = P.stringsize
    m = Q.stringsize
    return (n + 1, m + 1, matrixdata.matrix[n + 1, m + 1])
  end

  function backtracing{T}(P :: Profile{T}, Q :: Profile{T},
                          matrixdata :: AlignmentMatrix{T},
                          lastrow :: Int64, lastcol :: Int64)
    debugprint(matrixdata.path)
    path = Tuple{Char, Int64, Int64}[]
    (direction, row, col) = (matrixdata.path[lastrow, lastcol], lastrow, lastcol)
    while (row > 1 || col > 1)
      push!(path, (direction, row - 1, col - 1))
      debugprint(matrixdata.path[row, col])
      (row, col) = getprev(direction, row, col)
      direction = matrixdata.path[row, col]
    end
    path
  end

  construct{T}(a1 :: Array{T}, a2 :: Array{T}) = hcat(a1, a2)

  function mixprofilecolumn{T}(P :: Profile{T}, Q :: Profile{T},
                                 direction :: Char, i :: Int64, j :: Int64)
    direction == 'M' && return construct(P.rawdata[i, 1:end], Q.rawdata[j, 1:end])
    direction == 'R' && return construct(reshape(['-' for k in 1 : P.numberofstrings], 1, P.numberofstrings),
              Q.rawdata[j, 1 : end])
    direction == 'D' && return construct(
              P.rawdata[i, 1 : end],
              reshape(['-' for k in 1 : Q.numberofstrings], 1, Q.numberofstrings)
              )
    error("unknown direction in mix profile column")
  end

  function mixprofiles{T}(P :: Profile{T}, Q :: Profile{T}, indices :: Vector{Tuple{Char, Int64, Int64}})

    newProfileSize = length(indices)
    tempMatrix = Array(Char, newProfileSize, P.numberofstrings + Q.numberofstrings)
    for index in  1 : newProfileSize
      newProfileColumn = mixprofilecolumn(P, Q,
          indices[index][1], indices[index][2], indices[index][3])
      for j in 1 : length(newProfileColumn)
        tempMatrix[newProfileSize - index + 1, j] = newProfileColumn[j]
      end
    end

    Profile{T}(tempMatrix, append!(P.descriptions, Q.descriptions))
  end

  #
  # main alignment function for profiles.
  # returns consensus profile + alignment score
  #
  function align{T}(P :: Profile{T}, Q :: Profile{T})
    debugprint(P)
    debugprint(Q)
    n = P.stringsize
    m = Q.stringsize
    matrix = AlignmentMatrix{T}(n + 1, m + 1)
    prepareAlignmentMatrix(P, Q, matrix)
    buildAlignmentMatrix(P, Q, matrix, score)
    (lastrow, lastcol, sc) = getscore(P, Q, matrix)
    debugprint(sc)

    pairsofindices = backtracing(P, Q, matrix, lastrow, lastcol) # this should return new Profile{T} with resulting alignment and score
    debugprint(pairsofindices)
    res = mixprofiles(P, Q, pairsofindices)
    return res
  end

  function scoreprofiles{T}(P :: Profile{T}, Q :: Profile{T})
    debugprint(P)
    debugprint(Q)
    n = P.stringsize
    m = Q.stringsize
    matrix = AlignmentMatrix{T}(n + 1, m + 1)
    prepareAlignmentMatrix(P, Q, matrix)
    buildAlignmentMatrix(P, Q, matrix, score)
    (lastrow, lastcol, sc) = getscore(P, Q, matrix)
    debugprint(sc)
    return sc
  end

  #
  function measurequality{T}(P :: Profile{T})
    sum([
      sum([- Base.get(P.data[i], aa, zero(T))*log2(Base.get(P.data[i], aa, one(T)))
        for aa in setdiff(keys(P.data[i]), '-')])
      for i in 1 : P.stringsize])
  end

  #
  # method should return FastaRecord array for given profile
  #
  getstrings{T}(P :: Profile{T}) =  [ FastaRecord(P.descriptions[i], ascii(P.rawdata[1:end, i])) for i in 1 : P.numberofstrings]

end


text:
pyisnone(x) = pyis(x, pybuiltins.None)

pyconvert_rule_none(::Type{Nothing}, x::Py) = pyconvert_return(nothing)
pyconvert_rule_none(::Type{Missing}, x::Py) = pyconvert_return(missing)

pyconvert_rule_fast(::Type{Nothing}, x::Py) =
    if pyisnone(x)
        pyconvert_return(nothing)
    else
        pyconvert_unconverted()
    end

pyconvert_rule_fast(::Type{Missing}, x::Py) =
    if pyisnone(x)
        pyconvert_return(missing)
    else
        pyconvert_unconverted()
    end


text:
using Documenter, GLM

makedocs()


text:
# small scrip to compute the solution of Lipman Schinwer equation
# We test the types introduced in FastConvolution.jl
# we test that the application is fast and that the construction
# is performed fast.

using PyPlot
using IterativeSolvers
using SpecialFunctions

include("../src/FastConvolutionSlow.jl")
include("../src/FastConvolution.jl")
include("../src/Preconditioner.jl")

FFTW.set_num_threads(4);
BLAS.set_num_threads(4);

#Defining Omega
h = 0.001
k = 1/h

# size of box
a  = 1
x = -a/2:h:a/2
y = -a/2:h:a/2
(n,m) = length(x), length(y)
N = n*m
X = repmat(x, 1, m)[:]
Y = repmat(y', n,1)[:]
# we solve \triangle u + k^2(1 + nu(x))u = 0

# We use the modified quadrature in Ruan and Rohklin
(ppw,D) = referenceValsTrapRule();
D0 = D[1];

spline(y) = (y.<0) +  (y.>=0).*(y.<1).*(2*y.^3 - 3*y.^2 + 1)
filter_(a1,b1,b2,a2,y) = (y.>=a1).*( (y.<b1).*spline(1/(abs.(b1-a1)).*(-y+b1)) + (y.>=b1).*(y.<b2) + (y.>=b2).*spline(1/(abs.(b2-a2)).*(y-b2)))

fcut(y) = exp.(-1./y).*(y.>0)
fcut(y,alpha) = fcut(y)./(fcut(y) + fcut(alpha - y))
filter_(x,y) =  0*x + (y.<=(ylim[2]+h)).*(y.>=(ylim[1]-h)).*(1- fcut(abs.(y - width) - witdhy-h, width))  ; #add the rest of the filter_ in here
width = 0;
window(y,alpha) = 1*(abs.(y).<=alpha/2) + (abs.(y).>alpha/2).*(abs.(y).<alpha).*exp.(2*exp.(-0.5*alpha./(abs.(y)-alpha/2))./ ((abs.(y)-alpha/2)./(0.5*alpha)-1) )

# Defining the smooth perturbation of the slowness
nu(x,y) = 0.05*sin.(2*pi*x).*filter_(-0.48,-0.4,0.4,0.48,y).*filter_(-0.48,-0.05,0.05,0.48,x);

figure(1); clf();
imshow(real(reshape( (1-nu(X,Y)),n,m)), extent=[x[1], x[end], y[1], y[end]]);cb =  colorbar();

# Sampling the Green's function for the Toeplitz form
Ge = buildGConv(x,y,h,n,m,D0,k);
GFFT = fft(Ge);

fastconv     = FastM(GFFT,nu(X,Y),2*n-1,2*m-1,n, m, k);
fastconvslow = FastMslow(GFFT,nu(X,Y),X,Y,2*n-1,2*m-1,n, m, k, [1, 0]);

# assembling the sparsifiying preconditioner
As = buildSparseA(k,X,Y,D0, n ,m);

# assembling As*( I + k^2G*nu)
Mapproxsp = As + k^2*(buildSparseAG(k,X,Y,D0, n ,m)*spdiagm(nu(X,Y)));

# defining the preconditioner
precond = SparsifyingPreconditioner(Mapproxsp, As)

# building the RHS from the incident field
u_plane = exp.(k*im*sqrt(1-nu(0*Y,X)).*X);

u_inc = exp.(k*im*X);

rhs = -(fastconv*u_inc - u_inc);

# allocating the solution
u = zeros(Complex128,N);

# solving the system using GMRES
@time info =  gmres!(u, fastconv, rhs, Pl = precond)
# println(info[2].residuals[:])

# plotting the solution
figure(2)
clf()
imshow(real(reshape(u+u_inc,n,m)))

figure(3); clf()
imshow(real(reshape((u+u_inc)./u_inc, n,m))); colorbar()
title("\$u(x)/e^{ikx}\$ ")

intx = cumsum(sqrt(1-nu(x,0*x)))*h;
u_gem = exp.(k*im*intx);

UUgem = repmat(u_gem, 1,n)

figure(4); clf()
imshow(real(reshape((u+u_inc)./UUgem[:], n,m))); colorbar()
title("\$u(x)/e^{ik \\phi (x)}\$ where \$ \\phi \$ comes from \$ \\int_{0}^{x} \\sqrt{1 + n(y)} dy \$ ")

refFact = 10;
xRef = -a/2:h/refFact:a/2

intxRef = cumsum(sqrt(1-nu(xRef,0*xRef)))*h/refFact;
u_gem = exp.(k*im*intxRef)[1:refFact:end];

UUgem = repmat(u_gem, 1,n)

figure(14); clf()
imshow(real(reshape(-(u+u_inc)./UUgem[:], n,m))); colorbar()
title("\$u(x)/e^{ik \\phi (x)}\$ where \$ \\phi \$ comes from \$ \\int_{0}^{x} \\sqrt{1 + n(y)} dy \$ ")

UU = reshape(u+u_inc,n,m)
ucomp = UU[:,round(Integer, end/2)]

# there is a two that had to disappear... I need to find the bug... :/
phi_1 = -1im*rhs./u_inc

figure(5); clf()
imshow(real(reshape(phi_1, n,m))); colorbar()
title("phase function \$phi_1(x)\$ ")

u_rytov = u_inc.*exp.(1im*phi_1)

figure(6); clf()
imshow(real(reshape(u_rytov, n,m))); colorbar()
title("\$u_{rytov}\$" )

figure(7); clf()
imshow(real(reshape((u+u_inc)./ u_rytov, n,m))); colorbar()
title("\$u(x)/u_{rytov}(x)\$ ")

#Defining the smooth perturbation of the slowness
#nu1(x,y) = nu(x,y).*filter_(-0.28,-0.10,0.10,0.28,y);
nu1(x,y) = nu(x,y).*window(y, 0.4);

figure(24); clf();
imshow(real(reshape( (1-nu1(X,Y)),n,m)), extent=[x[1], x[end], y[1], y[end]]);cb =  colorbar();

fastconvWindow = FastM(GFFT,nu1(X,Y),2*n-1,2*m-1,n, m, k);

rhsWindow = -(fastconvWindow*u_inc - u_inc);

# there is a two that had to disappear... I need to find the bug... :/
phi_Window = -1im*rhsWindow./u_inc

figure(15); clf()
imshow(real(reshape(phi_1, n,m))); colorbar()
title("phase function phi_1(x) ")

u_rytovWindow = u_inc.*exp.(1im*phi_Window)

figure(16); clf()
imshow(real(reshape(u_rytovWindow, n,m))); colorbar()
title("u_rytov ")

figure(17); clf()
imshow(real(reshape((u+u_inc)./ u_rytovWindow, n,m))); colorbar()
title("u/u_rytov window ")

UU = reshape((u+u_inc)./u_rytov,n,m)
ucomp = UU[:,round(Integer, end/2)]

figure(20); clf();
plot(real(ucomp[:]))

rhsslow = -(fastconv*u_inc - u_inc)./u_inc;

# allocating the solution
uslow = zeros(Complex128,N);

@time info =  gmres!(uslow, fastconvslow, rhsslow);

UU = reshape((u)./(uslow.*exp.(1im*k*X)),n,m);
ucomp = UU[:,round(Integer, end/2)]

figure(21); clf();
plot(real(ucomp[:]))
title("\$u/(u_{slow}\\cdot e^{ikx})\$")

fastconvslowWindow = FastMslow(GFFT,nu1(X,Y),X,Y,3*n-2,3*m-2,n, m, k, [1 ,0]);

rhsslowWindow = -(fastconvWindow*u_inc - u_inc)./u_inc;

uslowWindow = zeros(Complex128,N);

@time info =  gmres!(uslowWindow , fastconvslowWindow, rhsslowWindow)

figure(30); clf();
imshow(real(reshape(uslowWindow, n,m))); colorbar()
title("u_slow with a windowed perturbation")

figure(31); clf();
imshow(real(reshape(uslow, n,m))); colorbar()
title("u_slow ")

UU = reshape(filter_(-0.48,-0.40,0.40,0.48,X).*(uslowWindow - uslow),n,m);
ucomp = UU[:,round(Integer, end/2)]

figure(32); clf();
plot(real(ucomp[:]))
title("uslow - u_slowwindow")

UU = reshape((u)./(uslowWindow.*exp.(1im*k*X)),n,m);
ucomp = UU[:,round(Integer, end/2)]

figure(23); clf();
plot(real(ucomp[:]))
title("u/u_slowwindow")

UU = reshape((u+u_inc)./(u_rytovWindow),n,m);
ucomp = UU[:,round(Integer, end/2)]

figure(29); clf();
plot(real(ucomp[:]))
title("u /u_rytovwindow")

convu =  -(fastconv*u - u);
convuslow = -(fastconvWindow*(uslowWindow.*exp.(1im*k*X)) - uslowWindow.*exp.(1im*k*X));

figure(26); clf();
 imshow(real(reshape(convu -convuslow, n,m))); colorbar()

UU = reshape(convu - convuslow,n,m);
ucomp = UU[:,round(Integer, end/2)]
figure(25); clf();
plot(real(ucomp[:]))
title("T_k (u- u_slowwindow)")

# u_inc = exp.(k*im*X);

# rhs1 = -(fastconv1*u_inc - u_inc);

# # allocating the solution
# u1 = zeros(Complex128,N);

# # solving the system using GMRES
# @time info =  gmres!(u1, fastconv1, rhs1, precond1)
# println(info[2].residuals[:])

# # plotting the solution
# figure(5)
# clf()
# imshow(real(reshape(u1+u_inc,n,m)))

# figure(6); clf()
# imshow(real(reshape((u1+u_inc)./u_inc, n,m))); colorbar()

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
"""
The base type used by all ephemeral random constants.
"""
abstract ERC

"""
Used to represent numerical ephemeral random constants.
"""
type NumericERC{N} <: ERC
  """
  The minimum value that the constant may take.
  """
  min::N

  """
  The maximum value that the constant may take.
  """
  max::N
end

"""
TODO ==========================================================================
This isn't right...
ERCs should be compiled and not called.
"""
call{N}(erc::NumericERC{N}) =
  rand(erc.min, erc.max)


text:
for test_file in readdir(joinpath(dirname(@__FILE__), "..", "examples/oldinterface"))
    # Do not test JuMP examples to avoid conflict between JuMP 0.18 and 0.19.
    if test_file[end-2:end] == ".jl" && !startswith(test_file, "jump")
        println("Running $test_file")
        include(joinpath(dirname(@__FILE__), "..", "examples/oldinterface", test_file))
    end
end


text:
import Obj
import Obj.Objective

type SvrgState <: State
    w_est::Array{Float64,1}     # w estimate
    gw_est::Array{Float64,1}    # μ (= full gradient of w estimate)
end

"""
    svrg([estiter, strategy])

    The estiter parameter describes the number of iterations before a new w estimate is
    obtained. The strategy to obtain a new w estimate is set by the strategy parameter.
    The following stategies exist:
    - :last     Snapshat of very last iteration (default)
    - :rand     Random snapshot of last iterations
    - :avg      Average of last iterations
"""
svrg(estiter::Int=10, strategy::Symbol=:last) = GdOptimizer(

    # name
    "svrg-"*string(estiter)*"-"*string(strategy),

    # initial state
    SvrgState(Float64[], Float64[]),    # correct values will be set in first iteration

    # update
    function (obj::Objective, k::Int, w::Array{Float64,1}, b::Array{Int,1}, inf, state::State)
        return _svrg_update(obj, k, w, b, inf, state, estiter, strategy)
    end

)

@fastmath function _svrg_update(obj::Objective, k::Int, w::Array{Float64,1}, b::Array{Int,1}, inf, state::State,
                                estiter::Int, strategy::Symbol)

    gw = Obj.g(obj, w, b)

    # svrg estimation update
    if (k-1) % estiter == 0
        if strategy == :last
            state.w_est  = w
        elseif strategy == :rand
            try
                state.w_est = inf[:w][ end-rand(1:estiter)+1 ]
            catch e
                state.w_est = w
            end
        elseif strategy == :avg
            try
                state.w_est = mean(inf[:w][ (end-estiter+1):end ], 2)
            catch e
                state.w_est = w
           end
        end
        state.gw_est = Obj.g(obj, state.w_est)
    end
    gw_est_b = Obj.g(obj, state.w_est, b)   # stochastic gradient of w estimate

   s = -gw + gw_est_b - state.gw_est # -g(w) + g(w_est) - μ

    return Obj.f(obj, w, b), gw, s, state

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
import YAML
import LinearAlgebra
import Base.length
import IterTools: groupby
using Printf
using Formatting

struct Structure
    lattice::Array{Float64,2}
    species::Array{String}
    fcoords::Array{Float64,2}
    replicate::Array{Int16}
    composition::Dict

    function Structure(lattice, species, fcoords, repl, comp)
        size(fcoords, 1) != length(species) && throw(ArgumentError("Length of species does not match the number of atomic positions"))
        new(lattice, species, fcoords, repl, comp)
    end
end

volume(s::Structure) = LinearAlgebra.det(s.lattice)
fractional(s::Structure) = s.fcoords
cartesian(s::Structure) = s.fcoords * s.lattice
length(s::Structure) = size(s.fcoords, 1)

slice_structure(s::Structure, indices::Array{Int, 1}) = Structure(s.lattice, s.species[indices], s.fcoords[indices, :], s.replicate, s.composition)
order(s::Structure, by=identity) = slice_structure(s, sortperm(s.species, by=by))

function sublattice(s::Structure, element::String)
    !(element in s.species) && throw(ArgumentError("Atomic species is $(element) part of the structure"))
    slice_structure(s, findall(el -> el == element , s.species))
end

exclude(s::Structure, lattice::String) = slice_structure(s, findall(el -> el != lattice , s.species))

function supercell(s::Structure, sort=false)
    scale = s.replicate
    inverse = collect(map(s -> 1.0/s, scale))
    num_sites = size(s.fcoords, 1)
    scaled_positions = hcat((inverse .* s.fcoords[row, :] for row in 1:num_sites)...)
    scale_matrix = [scale[1] 0 0; 0 scale[2] 0; 0 0 scale[3]]
    fcoords = vcat((
        (scaled_positions +
            hcat(repeat([inverse .* [i, j, k]], num_sites)...)
        )'
        for i=0:scale[1]-1
        for j=0:scale[2]-1
        for k=0:scale[3]-1
    )...)
    species = vcat((s.species for _ in 1:prod(scale))...)
    lattice = scale_matrix * s.lattice
    final = Structure(lattice, species, fcoords, s.replicate, s.composition)
    sort ? order(final) : final
end

function to_poscar(io::IO, s::Structure, direct=true, comment="automatically generated", format=".7f")
    @printf(io, "# %s\n", comment)
    @printf(io, "1.0\n")
    row_format(N) = join(("{$i:$format}" for i in 1:N), " ")
    print_row(row) = printfmtln(io, row_format(length(row)), row...)
    foreach(print_row, eachrow(s.lattice))
    species, amount =  T((g[1], length(g)) for g in groupby(identity, s.species))
    println(io, join(species, " "))
    println(io, join(amount, " "))
    println(io, direct ? "direct" : "cartesian")
    foreach(print_row, eachrow( direct ? fractional(s) : cartesian(s) ))
end

"""
The Struct 'Backpack' is a struct, which shall contain all variables that
have to be calculated just once. It serves to get rid of potentially disturbing
global variables.
"""
mutable struct Backpack
    supercell::Structure
    element_dict::Dict ## Dictionary with the element as a key and a tuple of
                       ## an assigned index and the occurrence of the element.
    mole_fractions::Dict ## Dictionary with all mole fractions (normed to sublattice)
    shell_matrix::Array ## Array with all relationshsips between two atoms
    triplet_shell_matrix::Array ## An array containing all three shell distances
                                ## between atoms i, j and k
    alpha::Array ## Pair SRO array
    beta::Array ## Triplet SRO array
    pair_prefactors::Array ##Array containg pair prefactors
    trip_prefactors::Dict ## Array containing the constant parts of the triplet
                            ##SRO parameter
end

max_stars_repo_name: ulipoto/WCSRO

max_issues_repo_name: ulipoto/WCSRO

max_forks_repo_name: ulipoto/WCSRO


text:
module KernelDensitySJ

using Statistics
using DataStructures
using Roots

export bwsj, smooth, density, GaussianKernelSmoother

include("sums.jl")
include("phi.jl")
include("bwsj.jl")
include("smooth.jl")

end


text:
#=
FastAI.jl:

Author: Peter Wolf (opus111@gmail.com)
=#

module FastAI

using Random
using StatsBase
using Statistics
using Flux
using Flux: update!
using Flux.Data
using Zygote
using Infiltrator
using Base: length, getindex
using Random: randperm
using DocStringExtensions

export AbstractLearner
export AbstractCallback
export IterableDataset
export MapDataset

export DataBunch
export train
export valid

export DummyCallback
export ProgressCallback

export Learner
export model
export data_bunch
export loss
export loss!
export opt
export opt!
export fit!
export add_cb!

export Recorder

@template (FUNCTIONS, METHODS) =
    """
    $(TYPEDSIGNATURES)
    $(DOCSTRING)
    """

@template (TYPES) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include("dataset.jl")
include("databunch.jl")
include("learner.jl")
include("callback.jl")
include("recorder.jl")

end


text:
module KhepriAutoCAD
using KhepriBase
using Sockets
using EzXML

# functions that need specialization
include(khepribase_interface_file())
include("AutoCAD.jl")

function __init__()
  add_current_backend(autocad)
end

end


text:
@testset "Measure" begin
    Mod.@polyvar x y
    @test_throws ArgumentError measure([1, 2], [x, x*y, y])
    @test_throws ArgumentError measure([1, 2, 3, 4], [x, x*y, y])
    μ = measure([1, 0, 2, 3], [x^2*y^2, y*x^2, x*y*x^2, x*y^2])
    @test monomials(μ) == [x^3*y, x^2*y^2, x^2*y, x*y^2]
    @test monomial.(moments(μ)) == [x^3*y, x^2*y^2, x^2*y, x*y^2]
    @test MultivariateMoments.moment_value.(moments(μ)) == [2, 1, 0, 3]
    @test all(nvariables.(moments(μ)) .== 2)
    @test degree.(moments(μ)) == [4, 4, 3, 3]
    @test μ.a == [2, 1, 0, 3]
    @test (-μ).a == [-2, -1, 0, -3]
    @test (2 * μ).a == [4, 2, 0, 6]
    @test (μ * 3).a == [6, 3, 0, 9]
    μ = measure([1, 1], [x, x])
    @test monomials(μ) == [x]
    @test MultivariateMoments.moment_value.(moments(μ)) == [1]
    err = ErrorException("The monomial `x` is occurs twice with different values: `1` and `2`")
    @test_throws err measure([1, 2], [x, x])
    #@test_throws ArgumentError measure([1], [x]) + measure([1], [y])
end


text:
using Vulkan
using Rhyolite
using ColorTypes
using FileIO
using Dictionaries
using Accessors
using Test

resource(filename) = joinpath(@__DIR__, "resources", filename)
render(filename) = joinpath(@__DIR__, "renders", filename)

function render_headless(output_png, points, colors; width = 1000, height = 1000)
    instance, device = init(; instance_extensions = ["VK_KHR_surface"], device_extensions = ["VK_KHR_synchronization2", "VK_KHR_swapchain"])
    config = dictionary([QUEUE_GRAPHICS_BIT | QUEUE_COMPUTE_BIT => 1])
    disp = QueueDispatch(device, config)
    pools = ThreadedCommandPool(device, disp, dictionary(1:1 .=> collect(keys(config))))
    shader_cache = ShaderCache(device)
    vert_shader = find_shader!(shader_cache, ShaderSpecification(resource("headless.vert"), GLSL))
    frag_shader = find_shader!(shader_cache, ShaderSpecification(resource("headless.frag"), GLSL))

    # compared to a display-based rendering, we have a lot more freedom over the format that we use.
    format = FORMAT_R32G32B32A32_SFLOAT
    format_props = get_physical_device_format_properties(handle(device).physical_device, format)

    if !(FORMAT_FEATURE_TRANSFER_SRC_BIT | FORMAT_FEATURE_TRANSFER_DST_BIT | FORMAT_FEATURE_COLOR_ATTACHMENT_BIT in format_props.optimal_tiling_features)
        error("Physical device $(get_physical_device_properties(handle(device).physical_device)) not supported")
    end

    target_attachment = AttachmentDescription(
        format,
        SAMPLE_COUNT_1_BIT,
        ATTACHMENT_LOAD_OP_CLEAR,
        ATTACHMENT_STORE_OP_STORE,
        ATTACHMENT_LOAD_OP_DONT_CARE,
        ATTACHMENT_STORE_OP_DONT_CARE,
        IMAGE_LAYOUT_UNDEFINED,
        IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL,
    )

    render_pass = RenderPass(
        device,
        [target_attachment],
        [
            SubpassDescription(
                PIPELINE_BIND_POINT_GRAPHICS,
                [],
                [AttachmentReference(0, IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL)],
                [],
            ),
        ],
        [
            SubpassDependency(
                SUBPASS_EXTERNAL,
                0;
                src_stage_mask = PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
                dst_stage_mask = PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
                dst_access_mask = ACCESS_COLOR_ATTACHMENT_WRITE_BIT,
            ),
        ],
    )

    # create image and framebuffer
    fb_image_info = ImageCreateInfo(
        IMAGE_TYPE_2D,
        format,
        Extent3D(width, height, 1),
        1,
        1,
        SAMPLE_COUNT_1_BIT,
        IMAGE_TILING_OPTIMAL,
        IMAGE_USAGE_COLOR_ATTACHMENT_BIT | IMAGE_USAGE_TRANSFER_SRC_BIT,
        SHARING_MODE_EXCLUSIVE,
        queue_family_indices(disp),
        IMAGE_LAYOUT_UNDEFINED,
    )
    fb_image = unwrap(create_image(device, fb_image_info))

    # some GPUs don't offer host-coherent host-visible memory for color attachments
    fb_image_memory = DeviceMemory(fb_image, MEMORY_PROPERTY_DEVICE_LOCAL_BIT)
    fb_image_view = ImageView(
        device,
        fb_image,
        IMAGE_VIEW_TYPE_2D,
        format,
        ComponentMapping(fill(COMPONENT_SWIZZLE_IDENTITY, 4)...),
        ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
    )
    framebuffer = Framebuffer(device, render_pass, [fb_image_view], width, height, 1)

    # prepare vertex buffer
    vdata = collect(zip(points, colors))
    vbuffer = Buffer(device, buffer_size(vdata), BUFFER_USAGE_VERTEX_BUFFER_BIT, SHARING_MODE_EXCLUSIVE, queue_family_indices(disp))
    vmemory = DeviceMemory(vbuffer, vdata)

    # build graphics pipeline
    shader_stage_cis = PipelineShaderStageCreateInfo.([vert_shader, frag_shader])
    vertex_input_state = PipelineVertexInputStateCreateInfo(
        [VertexInputBindingDescription(eltype(vdata), 0)],
        vertex_input_attribute_descriptions(eltype(vdata), 0, [FORMAT_R32G32_SFLOAT, FORMAT_R32G32B32_SFLOAT]),
    )
    input_assembly_state = PipelineInputAssemblyStateCreateInfo(PRIMITIVE_TOPOLOGY_TRIANGLE_STRIP, false)
    viewport_state = PipelineViewportStateCreateInfo(
        viewports = [Viewport(0, 0, width, height, 0, 1)],
        scissors = [Rect2D(Offset2D(0, 0), Extent2D(width, height))],
    )
    rasterizer = PipelineRasterizationStateCreateInfo(
        false,
        false,
        POLYGON_MODE_FILL,
        FRONT_FACE_CLOCKWISE,
        false,
        0.0,
        0.0,
        0.0,
        1.0,
        cull_mode = CULL_MODE_BACK_BIT,
    )
    multisample_state = PipelineMultisampleStateCreateInfo(SAMPLE_COUNT_1_BIT, false, 1.0, false, false)
    color_blend_attachment = PipelineColorBlendAttachmentState(
        true,
        BLEND_FACTOR_SRC_ALPHA,
        BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        BLEND_OP_ADD,
        BLEND_FACTOR_SRC_ALPHA,
        BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        BLEND_OP_ADD;
        color_write_mask = COLOR_COMPONENT_R_BIT | COLOR_COMPONENT_G_BIT | COLOR_COMPONENT_B_BIT,
    )
    color_blend_state = PipelineColorBlendStateCreateInfo(
        false,
        LOGIC_OP_CLEAR,
        [color_blend_attachment],
        Float32.((0.0, 0.0, 0.0, 0.0)),
    )
    pipeline_layout = PipelineLayout(device, [], [])
    (pipeline, _...), _ = unwrap(
        create_graphics_pipelines(
            device,
            [
                GraphicsPipelineCreateInfo(
                    shader_stage_cis,
                    rasterizer,
                    pipeline_layout,
                    render_pass,
                    0,
                    0;
                    vertex_input_state,
                    multisample_state,
                    color_blend_state,
                    input_assembly_state,
                    viewport_state,
                ),
            ],
        ),
    )

    dst_image_info = setproperties(fb_image_info, tiling = IMAGE_TILING_LINEAR, usage = IMAGE_USAGE_TRANSFER_DST_BIT)
    dst_image = unwrap(create_image(device, dst_image_info))
    dst_image_memory = DeviceMemory(dst_image, MEMORY_PROPERTY_HOST_COHERENT_BIT | MEMORY_PROPERTY_HOST_VISIBLE_BIT)

    # record commands
    command_buffer = first(unwrap(allocate_command_buffers(device, CommandBufferAllocateInfo(CommandPool(pools), COMMAND_BUFFER_LEVEL_PRIMARY, 1))))
    @record command_buffer begin
        cmd_bind_vertex_buffers([vbuffer], [0])
        cmd_bind_pipeline(PIPELINE_BIND_POINT_GRAPHICS, pipeline)
        cmd_pipeline_barrier(
            PIPELINE_STAGE_TOP_OF_PIPE_BIT,
            PIPELINE_STAGE_TRANSFER_BIT,
            [],
            [],
            [
                ImageMemoryBarrier(
                    AccessFlag(0),
                    ACCESS_MEMORY_READ_BIT,
                    IMAGE_LAYOUT_UNDEFINED,
                    IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
                    QUEUE_FAMILY_IGNORED,
                    QUEUE_FAMILY_IGNORED,
                    dst_image,
                    ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
                ),
            ],
        )
        cmd_begin_render_pass(
            RenderPassBeginInfo(
                render_pass,
                framebuffer,
                Rect2D(Offset2D(0, 0), Extent2D(width, height)),
                [ClearValue(ClearColorValue((0.1f0, 0.1f0, 0.15f0, 1.0f0)))],
            ),
            SUBPASS_CONTENTS_INLINE,
        )
        cmd_draw(length(vdata), 1, 0, 0)
        cmd_end_render_pass()
        cmd_pipeline_barrier(
            PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
            PIPELINE_STAGE_TRANSFER_BIT,
            [],
            [],
            [
                ImageMemoryBarrier(
                    ACCESS_COLOR_ATTACHMENT_WRITE_BIT,
                    ACCESS_MEMORY_READ_BIT,
                    IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL,
                    IMAGE_LAYOUT_TRANSFER_SRC_OPTIMAL,
                    QUEUE_FAMILY_IGNORED,
                    QUEUE_FAMILY_IGNORED,
                    fb_image,
                    ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
                ),
            ],
        )
        cmd_copy_image(
            fb_image,
            IMAGE_LAYOUT_TRANSFER_SRC_OPTIMAL,
            dst_image,
            IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
            [
                ImageCopy(
                    ImageSubresourceLayers(IMAGE_ASPECT_COLOR_BIT, 0, 0, 1),
                    Offset3D(0, 0, 0),
                    ImageSubresourceLayers(IMAGE_ASPECT_COLOR_BIT, 0, 0, 1),
                    Offset3D(0, 0, 0),
                    Extent3D(width, height, 1),
                ),
            ],
        )
    end

    # submit commands
    q = submit(disp, QUEUE_GRAPHICS_BIT, [SubmitInfo2KHR([], [CommandBufferSubmitInfoKHR(command_buffer, 0)], [])])
    GC.@preserve fb_image fb_image_view fb_image_memory framebuffer dst_image dst_image_memory vbuffer vmemory command_buffer pools begin
        unwrap(queue_wait_idle(q))
        data = download_data(Array{RGBA{Float32}}, dst_image_memory, (width, height))
    end

    save(output_png, data)
end

output_png = render("render_headless.png")

!ispath(output_png) || rm(output_png)

render_headless(output_png,
    NTuple{2,Float32}[
        (-1, -1),
        (1, -1),
        (-1, 1),
        (1, 1),
    ],
    RGB{Float32}[
        RGB(1,0,0),
        RGB(0,1,0),
        RGB(0,0,1),
        RGB(1,0,1),
    ]
)

@test stat(output_png).size > 200000


text:
module ThreadTest

	mutable struct Model
		a::Vector{Float64}
		r::Vector{Float64}
		id::Vector{Int}
	end

	mutable struct LModel
		b::Array
		function LModel(i::Int)
			new(ones(10,10))
		end
	end

	function doit0(n=10)
		m = Model(zeros(n),zeros(n),zeros(Int,n))
		Threads.@threads for i in 1:n
			z=ThreadTest.dothread!(m,i)
		end
		return m
	end

	function doit1(n=10)
		m = Model(zeros(n),rand(n),zeros(Int,n))
		x = 1.3
		Threads.@threads for i in 1:n
			Lm = LModel(i)  # create a local model
			lx = x + i
			id = Threads.threadid()  # note which id did which slot
			m.id[i] = id  # note which id did which slot
			m.a[i] = mean(Lm.b .* id)
		end
		return m,hcat(m.a,m.id)
	end
	function doit2()
		a = zeros(Int,10)
		y = zeros(10)
		x = 1.3
		Threads.@threads for i = 1:10
			lx = x + i
            a[i] = Threads.threadid()
            y[i] = lx
   		end
   		return (a,y)
	end

	function dothread!(m::Model,i::Int)
		Lm = LModel(i)  # create a local model
		# do some work on the local model and store the results in
		# the global model array
		# tmp = [Threads.threadid() for i in 1:14]
		m.id[i] = Threads.threadid()  # note which id did which slot
		m.a[i] = mean(Lm.b .* m.id[i] )
		m.r[i] = rand()
	end

end # module


text:
function EBFMI(Es::AbstractVector{<:AbstractScalarOrVec{<:AbstractFloat}})
    return mean(i -> (Es[i+1] - Es[i]) .^ 2, 1:length(Es)-1) ./ var(Es)
end


text:
"""
    `getindex(tbl::Table, fld::Symbol)`

Extract the column with field `fld` from `tbl`. Equivalent to `tbl[fld]`
"""
Base.getindex(tbl::Table, fld::Symbol) = columns(tbl)[index(tbl)[fld]]

"""
    `setindex(tbl::Table, col, fld::Symbol)`

Set `col` as the column respective to `fld` in `tbl`. Equivalent to
`tbl[fld] = col`

Notes: `col` will not be coerced. `length(col)` must equal `nrow(tbl)`.
"""
function Base.setindex!(tbl::Table, col::AbstractArray, field::Symbol)
    field = Symbol(replace_dots(string(field)))
    nrows, ncols = nrow(tbl), ncol(tbl)
    if (ncols > 0) & (length(col) != nrows)
        msg = "All columns in a Table must be the same length"
        throw(ArgumentError(msg))
    end
    j = get!(()->ncols+1, index(tbl), field)
    cols = columns(tbl)
    flds = fields(tbl)
    if j <= ncols
        cols[j] = convert(NullableArray, col)
        # tbl.allowsnulls[j] = _allowsnulls(col)
        tbl.hasnulls[j] = _hasnulls(col)
    else
        push!(cols, convert(NullableArray, col))
        push!(flds, field)
        # push!(tbl.allowsnulls, _allowsnulls(col))
        push!(tbl.hasnulls, _hasnulls(col))
    end
    return col
end


text:
using Test
using SymSpellChecker


text:
using LinuxPerf
using Test

using LinuxPerf: make_bench, enable!, disable!, reset!, reasonable_defaults, counters, EventType, EventTypeExt, parse_groups

@testset "LinuxPerf" begin

@testset "simple benchmark" begin
    @test begin
        bench = make_bench(reasonable_defaults);
        @noinline function g(a)
            enable!(bench)
            c = 0
            for x in a
                if x > 0
                    c += 1
                end
            end
            disable!(bench)
            c
        end
        g(zeros(10000))

        counters(bench)

        true  # Succeeded without any exceptions...
    end
end

@testset "@measure" begin
    c1 = @measure 2 + 2
    v, b1 = @measured 2 + 2

    @test v === 4
    @test typeof(c1) == typeof(counters(b1))
end

@testset "Parser" begin
    cycles = EventType(:hw, :cycles)
    insns = EventType(:hw, :instructions)

    @test parse_groups("") == []
    @test parse_groups("cpu-cycles") == [[EventTypeExt(cycles, false, 0)]]
    @test parse_groups("(cpu-cycles)") == parse_groups("cpu-cycles")
    @test parse_groups("cpu-cycles,instructions") == [[EventTypeExt(cycles, false, 0)], [EventTypeExt(insns, false, 0)]]
    @test parse_groups("(cpu-cycles,instructions)") == [[EventTypeExt(cycles, false, 0), EventTypeExt(insns, false, 0)]]
    @test parse_groups("  cpu-cycles,  instructions  ") == parse_groups("cpu-cycles,instructions")
    @test parse_groups("  (  cpu-cycles,  instructions  )  ") == parse_groups("(cpu-cycles,instructions)")

    # exclude flags
    u = LinuxPerf.exclude_flags(true, false, false)
    k = LinuxPerf.exclude_flags(false, true, false)
    h = LinuxPerf.exclude_flags(false, false, true)
    uk = LinuxPerf.exclude_flags(true, true, false)
    ukh = LinuxPerf.exclude_flags(true, true, true)

    # event-level modifiers
    @test parse_groups("cpu-cycles:u") == [[EventTypeExt(cycles, true, u)]]
    @test parse_groups("cpu-cycles:k") == [[EventTypeExt(cycles, true, k)]]
    @test parse_groups("cpu-cycles:h") == [[EventTypeExt(cycles, true, h)]]
    @test parse_groups("cpu-cycles:uk") == [[EventTypeExt(cycles, true, uk)]]
    @test parse_groups("cpu-cycles:ukh") == [[EventTypeExt(cycles, true, ukh)]]
    @test parse_groups("cpu-cycles:ku") == parse_groups("cpu-cycles:uk")
    @test parse_groups("cpu-cycles:uu") == parse_groups("cpu-cycles:u")
    @test parse_groups("cpu-cycles  :  u  ") == parse_groups("cpu-cycles:u")

    # group-level modifiers
    @test parse_groups("(cpu-cycles,instructions):u") == parse_groups("(cpu-cycles:u,instructions:u)")
    @test parse_groups("(cpu-cycles,instructions):k") == parse_groups("(cpu-cycles:k,instructions:k)")
    @test parse_groups("(cpu-cycles,instructions):h") == parse_groups("(cpu-cycles:h,instructions:h)")
    @test parse_groups("(cpu-cycles,instructions):uk") == parse_groups("(cpu-cycles:uk,instructions:uk)")
    @test parse_groups("(cpu-cycles,instructions):ukh") == parse_groups("(cpu-cycles:ukh,instructions:ukh)")
    @test parse_groups("(cpu-cycles,instructions):ku") == parse_groups("(cpu-cycles,instructions):uk")
    @test parse_groups("(cpu-cycles,instructions):uu") == parse_groups("(cpu-cycles,instructions):u")
    @test parse_groups("(cpu-cycles:k,instructions):u") == parse_groups("(cpu-cycles:k,instructions:u)")
    @test parse_groups("(cpu-cycles,instructions)  :  u  ") == parse_groups("(cpu-cycles,instructions):u")
end

@testset "@pstats" begin
    n = 10^3
    a = randn(n)
    b = randn(n)
    c = randn(n)
    dest = zero(a)

    function foo!(dest, a, b, c)
        @. dest = a + b * c
        sum(dest)
    end

    # Simple smoke tests
    @test_nowarn LinuxPerf.list()
    @pstats foo!(dest, a, b, c)

    @pstats "cpu-cycles,(instructions,branch-instructions,branch-misses),(task-clock,context-switches,cpu-migrations,page-faults),(L1-dcache-load-misses,L1-dcache-loads,L1-icache-load-misses),(dTLB-load-misses,dTLB-loads)" foo!(dest, a, b, c)
end

end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/19/2014

module CommonInterfaces

export
    addObserver,
    clearObservers!,

    initialize,
    step

function addObserver()
end

function clearObservers!()
end

function initialize()
end

function update()
end

end #module


text:
using TestPkgTemplates
using Documenter

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

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

deploydocs(;
    repo="github.com/rryi/TestPkgTemplates.jl",
)


text:
using SafeTestsets

include("sampler_tests.jl")
include("volume_tests.jl")

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
#=
# Linear elasticity
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("elasticitynew");

# Optionally generate a log
useLog("elasticitymatlablog", level=3)
generateFor(MATLAB)

n = [30,6,6]; # number of elements in x,y,z
interval = [0,1,0,0.2,0,0.2]; # domain bounds

# Set up the configuration
domain(3)
functionSpace(order=2)

# Specify the problem
mesh(HEXMESH, elsperdim=n, bids=4, interval=interval)

u = variable("u", type=VECTOR)          # make a vector variable with symbol u
testSymbol("v", type=VECTOR)            # sets the symbol for a test function

boundary(u, 1, DIRICHLET, [0,0,0]) # x=0
boundary(u, 2, NEUMANN, [0,0,0])   # elsewhere
boundary(u, 3, NEUMANN, [0,0,0])
boundary(u, 4, NEUMANN, [0,0,0])

# Write the weak form
# coefficient("mu", "x>0.5 ? 0.2 : 10") # discontinuous mu
coefficient("mu", 1) # constant mu
coefficient("lambda", 1.25)
coefficient("f", [0,0,-10], type=VECTOR)
weakForm(u, "inner( (lambda * div(u) .* [1 0 0; 0 1 0; 0 0 1] + mu .* (grad(u) + transpose(grad(u)))), grad(v)) - dot(f,v)")

exportCode("elasticitycode")

solve(u);

finalize_finch() # Finish writing and close any files

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# Tag model and embed data
#
export Model2Fit_TaggedModel
export get_tagged_data, get_tagged_model
export get_tagged_data_type, get_tagged_model_type

@doc raw"""
```julia
Model2Fit_TaggedModel(model,data)
```

Tag model and embed data

# Extra method

- [`get_tagged_data`](@ref)
- [`get_tagged_model`](@ref)
- [`get_tagged_data_type`](@ref)
"""
struct Model2Fit_TaggedModel{MODEL<:Abstract_Model2Fit,DATA} <: Abstract_Model2Fit
    _model::MODEL
    _data::DATA
end

# Extra methods  ================
#
@doc raw"""
```julia
get_tagged_data(m::Model2Fit_TaggedModel{MODEL,DATA)::DATA
```

Return embedded data
"""
get_tagged_data(m::Model2Fit_TaggedModel) = m._data
@doc raw"""
```julia
get_tagged_model(m::Model2Fit_TaggedModel{MODEL,DATA)::DATA
```

Return tagged model
"""
get_tagged_model(m::Model2Fit_TaggedModel) = m._model

@doc raw"""
```julia
get_tagged_data_type(m::Abstract_Model2Fit)::DataType
```

Return embedded data type, or `Nothing` if m is not a
[`Model2Fit_TaggedModel`](@ref).
"""
get_tagged_data_type(m::Abstract_Model2Fit) = Nothing
get_tagged_data_type(m::Model2Fit_TaggedModel{MODEL,DATA}) where {MODEL<:Abstract_Model2Fit,DATA} = DATA

@doc raw"""
```julia
get_tagged_model_type(m::Abstract_Model2Fit)::DataType
```

Return wrapped model type, or `Nothing` if m is not a
[`Model2Fit_TaggedModel`](@ref).
"""
get_tagged_model_type(m::Abstract_Model2Fit) = Nothing
get_tagged_model_type(m::Model2Fit_TaggedModel{MODEL,DATA}) where {MODEL<:Abstract_Model2Fit,DATA} = MODEL

# Visit  ================
#
visit_submodel_size(model::Model2Fit_TaggedModel) = 1
visit_get_submodel(model::Model2Fit_TaggedModel,submodel_idx::Int) = get_tagged_model(model)
visit_get_Y(model::Model2Fit_TaggedModel,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = Y
visit_get_X(model::Model2Fit_TaggedModel,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = X
visit_get_θ(model::Model2Fit_TaggedModel,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = θ

# Interface  ================
#
parameter_size(m::Model2Fit_TaggedModel) = parameter_size(visit_get_submodel(m,1))

function accumulate_y!(m::Model2Fit_TaggedModel,Y::AbstractVector,X::AbstractVector,θ::AbstractVector)
    @assert length(θ) == parameter_size(m)
    @assert length(X) == length(Y)

    accumulate_y!(get_tagged_model(m),Y,X,θ)
end