text:
using NiLang, Test

@testset "partial" begin
    x = Partial{:im}(3+2im)
    println(x)
    @test x === Partial{:im,Complex{Int64},Int64}(3+2im)
    @test value(x) == 2
    @test chfield(x, value, 4) == Partial{:im}(3+4im)
    @test zero(x) == Partial{:im}(0.0+0.0im)
    @test (~Partial{:im})(x) == 3+2im
end

@testset "value" begin
    x = 1.0
    @test value(x) === 1.0
    @assign (x |> value) 0.2
    @test x == 0.2
end

struct NiTypeTest{T} <: IWrapper{T}
    x::T
    g::T
end
NiTypeTest(x) = NiTypeTest(x, zero(x))
@fieldview NiLang.value(invtype::NiTypeTest) = invtype.x
@fieldview gg(invtype::NiTypeTest) = invtype.g

@testset "inv type" begin
    it = NiTypeTest(0.5)
    @test eps(typeof(it)) === eps(Float64)
    @test value(it) == 0.5
    @test it ≈ NiTypeTest(0.5)
    @test it > 0.4
    @test it < NiTypeTest(0.6)
    @test it < 7
    @test 0.4 < it
    @test 7 > it
    @test chfield(it, value, 0.3) == NiTypeTest(0.3)
    it = chfield(it, Val(:g), 0.2)
    @test almost_same(NiTypeTest(0.5+1e-15), NiTypeTest(0.5))
    @test !almost_same(NiTypeTest(1.0), NiTypeTest(1))
    it = NiTypeTest(0.5)
    @test chfield(it, gg, 0.3) == NiTypeTest(0.5, 0.3)
end


text:
using Printf
using SparseArrays

import PyPlot
const plt = PyPlot

include("../LF_common/m_LF1d.jl")
include("../LF_common/m_LF3d.jl")

function test_main( NN::Array{Int64} )
  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]
  LF = init_LF3d_c( NN, AA, BB )

  Npoints = prod(NN)
  println("Step 1: ")

  NNZ = ccall( (:get_nnz_laplacian3d_, "./Laplace3d_sparse.so"), Int64,
         (Ptr{Int64}, Ref{Int64}, Ptr{Int64},
          Ptr{Float64}, Ptr{Float64}, Ptr{Float64} ),
         NN, Npoints, LF.lin2xyz, LF.LFx.D2jl, LF.LFy.D2jl, LF.LFz.D2jl )

  #
  row    = zeros(Int64,NNZ)
  column = zeros(Int64,NNZ)
  values = zeros(Float64,NNZ)
  #
  println("Step 2: ")

  ccall( (:init_laplacian3d_, "./Laplace3d_sparse.so"), Nothing,
         ( Ref{Int64}, Ptr{Int64}, Ref{Int64}, Ptr{Int64},
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Int64}, Ptr{Int64}, Ptr{Float64} ),
         NNZ, NN, Npoints, LF.lin2xyz, LF.LFx.D2jl, LF.LFy.D2jl, LF.LFz.D2jl,
         row, column, values )

  #
  #for ii = 1:NNZ
  #  @printf("%5d %5d %18.10f\n", row[ii], column[ii], values[ii])
  #end
  println("Step 3: ")

  Laplacian = sparse( row, column, values )

  plt.clf()
  plt.spy(Matrix(Laplacian))
  filename = "Laplacian_"*string(NN[1])*"_"*string(NN[2])*"_"*string(NN[3])*".pdf"
  plt.savefig(filename)

end

test_main([3,5,5])
#test_main([20,20,20])
#test_main([30,30,30])
#test_main([40,40,40])
#test_main([50,50,50])
#for ii = 10:10:60
#  @time test_main( ii*ones(Int64,3) )
#end

ext: jl


text:
# GaussianVectors.jl

# Define and simulate normal vectors.
struct GaussianVector{R<:AbstractArray{Float64, 1}, S<:AbstractArray{Float64},
		T<:AbstractArray{Float64}}
	mean::R
	Σ::S
	var::T
	length::Int
	numberofGaussians::Int

	function GaussianVector(mean::R, Σ::S) where {R<:AbstractArray{Float64, 1}, S<:AbstractArray{Float64}}
		Base.depwarn("GaussianVector is deprecated. Use MvNormal or MvNormalCanon from the Distributions package instead.", :GaussianVector, force=true)
		length = Base.length(mean)
		size(Σ,1) == length || throw(DimensionMismatch("Dimensions of Σ and μ don't match."))
		numberofGaussians = size(Σ,2)
		var = Σ * Σ'
		new{R, S, typeof(var)}(mean, Σ, var, length, numberofGaussians)
	end
end

GaussianVector(Σ::T) where T<:AbstractArray{Float64} = GaussianVector(spzeros(size(Σ,1)), Σ)

rand(d::GaussianVector) = d.mean + d.Σ * randn((d.numberofGaussians, 1))

mean(d::GaussianVector) = d.mean

var(d::GaussianVector) = d.var

cov(d::GaussianVector) = Matrix(d.var)

invcov(d::GaussianVector) = inv(cholesky(cov(d)))

length(d::GaussianVector) = d.length

scaleGV(d::T, scale::Float64) where {T<:GaussianVector} = GaussianVector(scale * d.mean, scale * d.Σ)


text:
# function generate_windturbine(Rtip::Float64, h::Float64, blade_name::String,
#     hub_name::String, tower_name::String;
#     nblades::Int64=3, data_path::String=def_data_path,
#     save_path="/Users/dprocell/WF/files", file_name="windturbine",
#     paraview=true,
#     rot=nothing, random_rot::Bool=true, pitch=70)

H = 150.0
D = 200.0
nblades = 3

blade = "NREL5MW"
hub = "hub"
tower = "tower1"

module_path = splitdir(@__FILE__)[1]
data_path = joinpath(module_path, "../datav07/")

save_path = "test"
gt.create_path(save_path, true)

generate_windturbine(D/2, H, blade, hub, tower;
                    nblades=nblades, data_path=data_path,
                    save_path=save_path, file_name="windturbine",
                    paraview=false)

#end


text:
using Distributions, Plots; pyplot()

mu, sigma = 10, 4
n, N = 10, 10^6

sMeans = Array{Float64}(undef, N)
sVars  = Array{Float64}(undef, N)
tStats = Array{Float64}(undef, N)

for i in 1:N
    data       = rand(Normal(mu,sigma),n)
    sampleMean = mean(data)
    sampleVars = var(data)
    sMeans[i]  = sampleMean
    sVars[i]   = sampleVars
    tStats[i]  = (sampleMean - mu)/(sqrt(sampleVars/n))
end

xRangeMean = 5:0.1:15
xRangeVar = 0:0.1:60
xRangeTStat = -5:0.1:5

p1 = stephist(sMeans, bins=50, c=:blue, normed=true, legend=false)
p1 = plot!(xRangeMean, pdf.(Normal(mu,sigma/sqrt(n)), xRangeMean),
    c=:red, xlims=(5,15), ylims=(0,0.35), xlabel="Sample mean",ylabel="Density")

p2 = stephist(sVars, bins=50, c=:blue, normed=true, label="Simulated")
p2 = plot!(xRangeVar, (n-1)/sigma^2*pdf.(Chisq(n-1), xRangeVar*(n-1)/sigma^2),
    c=:red, label="Analytic", xlims=(0,60), ylims=(0,0.06),
    xlabel="Sample Variance",ylabel="Density")

p3 = stephist(tStats, bins=100, c=:blue, normed=true, legend=false)
p3 = plot!(xRangeTStat, pdf.(TDist(n-1), xRangeTStat),
    c=:red, xlims=(-5,5), ylims=(0,0.4), xlabel="t-statistic",ylabel="Density")

plot(p1, p2, p3, layout = (1,3), size=(1200, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
__precompile__()

"""A collection of physical constants used in the SoilTracers model."""
module Constants

# do not export anything; preserve the namespace

"""Zero Celsius in Kelvin."""
const T_0 = 273.15

"""Molar gas constant [J K^-1 mol^-1]."""
const R = 8.314_4621

"""Standard acceleration of gravity on earth [m s^-2]."""
const g = 9.806_65

"""Standard atmospheric pressure [Pa]."""
const atm = 1.01325e5

"""Molar mass of water [kg mol^-1]."""
const M_w = 18.015_28e-3

end


text:
export MtlComputeCommandEncoder
export set_function!, set_buffer!, set_bytes!, dispatchThreads!, endEncoding!
export set_buffers!, append_current_function!

const MTLComputeCommandEncoder = Ptr{MtComputeCommandEncoder}

mutable struct MtlComputeCommandEncoder <: MtlCommandEncoder
    handle::MTLComputeCommandEncoder
    cmdbuf::MtlCommandBuffer
end

Base.unsafe_convert(::Type{MTLComputeCommandEncoder}, q::MtlComputeCommandEncoder) = q.handle

function MtlComputeCommandEncoder(cmdbuf::MtlCommandBuffer; dispatch_type::Union{Nothing,MtDispatchType} = nothing)
    if isnothing(dispatch_type)
        handle = mtNewComputeCommandEncoder(cmdbuf)
    else
        handle = mtNewComputeCommandEncoderWithDispatchtype(cmdbuf, dispatchtype)
    end
    obj = MtlComputeCommandEncoder(handle, cmdbuf)
    #finalizer(unsafe_destroy!, obj)
    return obj
end

device(cce::MtlComputeCommandEncoder) = cce.cmdbuf.device

set_function!(cce::MtlComputeCommandEncoder, pip::MtlComputePipelineState) =
    mtComputeCommandEncoderSetComputePipelineState(cce, pip)

set_buffer!(cce::MtlComputeCommandEncoder, buf::MtlBuffer, offset::Integer, index::Integer) =
    mtComputeCommandEncoderSetBufferOffsetAtIndex(cce, buf, offset, index - 1)
#set_bufferoffset!(cce::MtlComputeCommandEncoder, offset::Integer, index::Integer) =
#    mtComputeCommandEncoderBufferSetOffsetAtIndex(cce, offset, index)
set_buffers!(cce::MtlComputeCommandEncoder, bufs::Vector{T},
             offsets::Vector{Int}, indices::UnitRange{Int}) where {T<:MtlBuffer} =
    mtComputeCommandEncoderSetBuffersOffsetsWithRange(cce, handle_array(bufs), offsets, indices .- 1)
#=set_buffers!(cce::MtlComputeCommandEncoder, bufs::Vector{MtlPtr{T}},
             offsets::Vector{Int}, indices::UnitRange{Int}) where {T} =
    mtComputeCommandEncoderSetBuffersOffsetsWithRange(cce, bufs, offsets, indices .- 1)=#

set_bytes!(cce::MtlComputeCommandEncoder, ptr, length::Integer, index::Integer) =
    mtComputeCommandEncoderSetBytesLengthAtIndex(cce, ptr, length, index - 1)

dispatchThreads!(cce::MtlComputeCommandEncoder, gridSize::MtSize, threadGroupSize::MtSize) =
    mtComputeCommandEncoderDispatchThreadgroups_threadsPerThreadgroup(cce, gridSize, threadGroupSize)

#####
# encode in the Command Encoder
function MtlComputeCommandEncoder(f::Base.Callable, cmdbuf::MtlCommandBuffer; kwargs...)
    encoder = MtlComputeCommandEncoder(cmdbuf; kwargs...)
    f(encoder)
    close(encoder)
    return encoder
end

append_current_function!(cce::MtlComputeCommandEncoder, gridSize::MtSize, threadGroupSize::MtSize) =
    dispatchThreads!(cce, gridSize, threadGroupSize)

#### use
use!(cce::MtlComputeCommandEncoder, buf::MtlBuffer, mode::MtResourceUsage=ReadWriteUsage) =
    mtComputeCommandEncoderUseResourceUsage(cce, buf, mode)

use!(cce::MtlComputeCommandEncoder, buf::Vector{MtlBuffer}, mode::MtResourceUsage=ReadWriteUsage) =
    mtComputeCommandEncoderUseResourceCountUsage(cce, handle_array(buf), length(buf), mode)


text:
using Documenter
using Kalman

makedocs(
    modules = [Kalman],
    format = :html,
    sitename = "Kalman.jl",
    authors = "Moritz Schauer and contributors",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
#        "Manual" => "manual.md",
        "Library" => "library.md",
        ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/mschauer/Kalman.jl.git",
    julia  = "0.6",
    target = "build",
    deps = nothing,
    make = nothing,
)


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

"""
Unit tests for the `EnhancedTestSet` type.

This set of unit tests checks the behavior of `EnhancedTestSet` when the test set type is
`EnhancedTestSet{FallbackTestSet}` (i.e., fail-fast).
"""

# --- Imports

# Standard library
using Test
using Test: DefaultTestSet, FallbackTestSet, FallbackTestSetException

# External packages
using Suppressor

# Local modules
using TestTools.jltest

# --- Tests

@testset "EnhancedTestSet{FallbackTestSet} Tests" begin
    test_set_type = EnhancedTestSet{FallbackTestSet}

    # Single-level test set with a failing test
    error_type = Nothing
    error_message = ""
    output = strip(@capture_out begin
        try
            @testset test_set_type "top-level tests" begin
                @test 1 == 2
                @test 1 == 1
            end
        catch error
            error_type = typeof(error)
            error_message = error.msg
            println(error)
        end
    end)

    @test error_type == EnhancedTestSetException
    @test error_message == "FallbackTestSetException occurred"

    expected_prefix = strip("""
                            =====================================================
                            Test Failed at $(@__FILE__):45
                              Expression: 1 == 2
                               Evaluated: 1 == 2
                            """)
    @test startswith(output, expected_prefix)

    # Nested test sets with a failing test
    error_type = Nothing
    error_message = ""
    output = strip(@capture_out begin
        try
            @testset test_set_type "top-level tests" begin
                @testset "Test set with failing test" begin
                    @test 1 == 2
                    @test 1 == 1
                end
                @testset "Test set with no failing tests" begin
                    @test 2 == 2
                    @test 3 == 3
                end
            end
        catch error
            error_type = typeof(error)
            error_message = error.msg
        end
    end)

    @test error_type == EnhancedTestSetException
    @test error_message == "FallbackTestSetException occurred"

    expected_prefix = strip("""
                            =====================================================
                            Test Failed at $(@__FILE__):73
                              Expression: 1 == 2
                               Evaluated: 1 == 2
                            """)
    @test startswith(output, expected_prefix)

    # --- Nested DefaultTestSet tests
    #
    # * Tests Test.record(EnhancedTestSet{FallbackTestSet}, DefaultTestSet) needed for
    #   backward compatibility with Julia<=1.3.

    default_test_set = DefaultTestSet

    # ------ DefaultTest nested under single EnhancedTestSet{FallbackTestSet} test set

    # With failing tests
    error_type = Nothing
    error_message = ""
    output = (@capture_out begin
        try
            @testset test_set_type "top-level tests" begin
                @testset default_test_set "Failing test" begin
                    @test 1 == 2
                    @test 1 == 1
                end
                @testset default_test_set "No failing tests" begin
                    @test 2 == 2
                    @test 3 == 3
                end
            end
        catch error
            error_type = typeof(error)
            error_message = error.msg
        end
    end)

    @test error_type == FallbackTestSetException

    expected_prefix = strip("""
                            Failing test: Test Failed at $(@__FILE__):114
                              Expression: 1 == 2
                               Evaluated: 1 == 2
                            """)
    @test startswith(output, expected_prefix)

    # With no failing tests
    error = nothing
    try
        @testset test_set_type "top-level tests" begin
            @testset default_test_set "No failing tests" begin
                @test 2 == 2
                @test 3 == 3
            end
        end
    catch error
    end

    @test isnothing(error)

    # ------ DefaultTest nested under multiple EnhancedTestSet{FallbackTestSet} test sets

    # With failing tests
    error_type = Nothing
    error_message = ""
    output = @capture_out begin
        try
            @testset test_set_type "top-level tests" begin
                @testset test_set_type "2nd-level tests" begin
                    @testset default_test_set "Failing test" begin
                        @test 1 == 2
                        @test 1 == 1
                    end
                    @testset default_test_set "No failing tests" begin
                        @test 2 == 2
                        @test 3 == 3
                    end
                end
            end
        catch error
            error_type = typeof(error)
            error_message = error.msg
        end
    end

    @test error_type == EnhancedTestSetException
    @test error_message == "FallbackTestSetException occurred"

    # With no failing tests
    error = nothing
    try
        @testset test_set_type "top-level tests" begin
            @testset test_set_type "2nd-level tests" begin
                @testset default_test_set "No failing tests" begin
                    @test 2 == 2
                    @test 3 == 3
                end
            end
        end
    catch error
    end

    @test isnothing(error)
end


text:
"""
parse_tuple(ps, ret)

`ret` is followed by a comma so tries to parse the rest of the
tuple.
"""
function parse_tuple end

@static if VERSION > v"1.1-"
    function parse_tuple(ps::ParseState, ret::EXPR)
        op = mPUNCTUATION(next(ps))
        if typof(ret) == TupleH
            if (isassignment(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
                push!(ret, op)
            elseif closer(ps)
                ps.errored = true
                push!(ret, mErrorToken(op, Unknown))
            else
                nextarg = @closer ps :tuple parse_expression(ps)
                if !(is_lparen(first(ret.args)))
                    push!(ret, op)
                    push!(ret, nextarg)
                else
                    ret = EXPR(TupleH, EXPR[ret, op, nextarg])
                end
            end
        else
            if (isassignment(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
                ret = EXPR(TupleH, EXPR[ret, op])
            elseif closer(ps)
                ps.errored = true
                ret = mErrorToken(EXPR(TupleH, EXPR[ret, op]), Unknown)
            else
                nextarg = @closer ps :tuple parse_expression(ps)
                ret = EXPR(TupleH, EXPR[ret, op, nextarg])
            end
        end
        return ret
    end
else
    function parse_tuple(ps::ParseState, ret::EXPR)
        op = mPUNCTUATION(next(ps))
        if typof(ret) === TupleH
            if closer(ps) || (isassignment(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
                push!(ret, op)
            else
                nextarg = @closer ps :tuple parse_expression(ps)
                if !(is_lparen(first(ret.args)))
                    push!(ret, op)
                    push!(ret, nextarg)
                else
                    ret = EXPR(TupleH, EXPR[ret, op, nextarg])
                end
            end
        else
            if closer(ps) || (isassignment(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
                ret = EXPR(TupleH, EXPR[ret, op])
            else
                nextarg = @closer ps :tuple parse_expression(ps)
                ret = EXPR(TupleH, EXPR[ret, op, nextarg])
            end
        end
        return ret
    end

end

"""
    parse_array(ps)
Having hit '[' return either:
+ A vect
+ A vcat
+ A comprehension
+ An array (vcat of hcats)
"""
function parse_array(ps::ParseState, isref = false)
    args = EXPR[mPUNCTUATION(ps)]

    if kindof(ps.nt) == Tokens.RSQUARE
        accept_rsquare(ps, args)
        ret = EXPR(Vect, args)
    else
        first_arg = @nocloser ps :newline @closesquare ps  @closer ps :insquare @closer ps :ws @closer ps :wsop @closer ps :comma parse_expression(ps)
        if isref && _do_kw_convert(ps, first_arg)
            first_arg = _kw_convert(first_arg)
        end

        if kindof(ps.nt) == Tokens.RSQUARE
            if typof(first_arg) === Generator || typof(first_arg) === Flatten
                accept_rsquare(ps, args)

                if typof(first_arg.args[1]) === BinaryOpCall && is_pairarrow(first_arg.args[1].args[2])
                    return EXPR(DictComprehension, EXPR[args[1], first_arg, INSTANCE(ps)])
                else
                    return EXPR(Comprehension, EXPR[args[1], first_arg, INSTANCE(ps)])
                end
            elseif kindof(ps.ws) == SemiColonWS
                push!(args, first_arg)
                accept_rsquare(ps, args)
                return EXPR(Vcat, args)
            else
                push!(args, first_arg)
                accept_rsquare(ps, args)
                ret = EXPR(Vect, args)
            end
        elseif kindof(ps.nt) == Tokens.COMMA
            etype = Vect
            push!(args, first_arg)
            accept_comma(ps, args)
            @closesquare ps parse_comma_sep(ps, args, isref)
            accept_rsquare(ps, args)
            return EXPR(etype, args)
        elseif kindof(ps.ws) == NewLineWS
            ps.closer.inref = false
            ret = EXPR(Vcat, args)
            push!(ret, first_arg)
            while kindof(ps.nt) != Tokens.RSQUARE
                if kindof(ps.nt) == Tokens.ENDMARKER
                    break
                end
                a = @closesquare ps  parse_expression(ps)
                push!(ret, a)
            end
            accept_rsquare(ps, ret)
            update_span!(ret)
            return ret
        elseif kindof(ps.ws) == WS || kindof(ps.ws) == SemiColonWS
            ps.closer.inref = false
            first_row = EXPR(Hcat, EXPR[first_arg])
            while kindof(ps.nt) != Tokens.RSQUARE && kindof(ps.ws) != NewLineWS && kindof(ps.ws) != SemiColonWS
                if kindof(ps.nt) == Tokens.ENDMARKER
                    break
                end
                a = @closesquare ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                push!(first_row, a)
            end
            if kindof(ps.nt) == Tokens.RSQUARE && kindof(ps.ws) != SemiColonWS
                if length(first_row.args) == 1
                    first_row = EXPR(Vcat, first_row.args)
                end
                push!(first_row, INSTANCE(next(ps)))
                pushfirst!(first_row, args[1])
                update_span!(first_row)
                return first_row
            else
                if length(first_row.args) == 1
                    first_row = first_row.args[1]
                else
                    first_row = EXPR(Row, first_row.args)
                end
                ret = EXPR(Vcat, EXPR[args[1], first_row])
                while kindof(ps.nt) != Tokens.RSQUARE
                    if kindof(ps.nt) == Tokens.ENDMARKER
                        break
                    end
                    first_arg = @closesquare ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                    push!(ret, EXPR(Row, EXPR[first_arg]))
                    while kindof(ps.nt) != Tokens.RSQUARE && kindof(ps.ws) != NewLineWS && kindof(ps.ws) != SemiColonWS
                        if kindof(ps.nt) == Tokens.ENDMARKER
                            break
                        end
                        a = @closesquare ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                        push!(last(ret.args), a)
                    end
                    # if only one entry dont use :row
                    if length(last(ret.args).args) == 1
                        ret.args[end] = setparent!(ret.args[end].args[1], ret)
                    end
                    update_span!(ret)
                end
                accept_rsquare(ps, ret)
                update_span!(ret)
                return ret
            end
        else
            ret = EXPR(Vect, args)
            push!(ret, first_arg)
            push!(ret, accept_rsquare(ps))
        end
    end
    return ret
end

"""
    parse_ref(ps, ret)

Handles cases where an expression - `ret` - is followed by
`[`. Parses the following bracketed expression and modifies it's
`.head` appropriately.
"""
function parse_ref(ps::ParseState, ret::EXPR)
    next(ps)
    ref = @closer ps :inref @nocloser ps :inwhere parse_array(ps, true)
    if typof(ref) === Vect
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(Ref, args)
    elseif typof(ref) === Hcat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(TypedHcat, args)
    elseif typof(ref) === Vcat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(TypedVcat, args)
    else
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(TypedComprehension, args)
    end
end

"""
parse_curly(ps, ret)

Parses the juxtaposition of `ret` with an opening brace. Parses a comma
seperated list.
"""
function parse_curly(ps::ParseState, ret::EXPR)
    args = EXPR[ret, mPUNCTUATION(next(ps))]
    parse_comma_sep(ps, args, true)
    accept_rbrace(ps, args)
    return EXPR(Curly, args)
end

function parse_braces(ps::ParseState)
    return @default ps @nocloser ps :inwhere parse_barray(ps)
end

function parse_barray(ps::ParseState)
    args = EXPR[mPUNCTUATION(ps)]

    if kindof(ps.nt) == Tokens.RBRACE
        accept_rbrace(ps, args)
        ret = EXPR(Braces, args)
    else
        first_arg = @nocloser ps :newline @closebrace ps  @closer ps :ws @closer ps :wsop @closer ps :comma parse_expression(ps)
        if kindof(ps.nt) == Tokens.RBRACE
            push!(args, first_arg)
            if kindof(ps.ws) == SemiColonWS
                push!(args, EXPR(Parameters, EXPR[]))
            end
            accept_rbrace(ps, args)
            ret = EXPR(Braces, args)
        elseif kindof(ps.nt) == Tokens.COMMA
            push!(args, first_arg)
            accept_comma(ps, args)
            @closebrace ps parse_comma_sep(ps, args, true)
            accept_rbrace(ps, args)
            return EXPR(Braces, args)
        elseif kindof(ps.ws) == NewLineWS
            ret = EXPR(BracesCat, args)
            push!(ret, first_arg)
            while kindof(ps.nt) != Tokens.RBRACE
                if kindof(ps.nt) == Tokens.ENDMARKER
                    break
                end
                a = @closebrace ps  parse_expression(ps)
                push!(ret, a)
            end
            accept_rsquare(ps, ret)
            update_span!(ret)
            return ret
        elseif kindof(ps.ws) == WS || kindof(ps.ws) == SemiColonWS
            first_row = EXPR(Row, EXPR[first_arg])
            while kindof(ps.nt) != Tokens.RBRACE && kindof(ps.ws) != NewLineWS && kindof(ps.ws) != SemiColonWS
                if kindof(ps.nt) == Tokens.ENDMARKER
                    break
                end
                a = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                push!(first_row, a)
            end
            if kindof(ps.nt) == Tokens.RBRACE && kindof(ps.ws) != SemiColonWS
                if length(first_row.args) == 1
                    first_row = EXPR(BracesCat, first_row.args)
                end
                push!(args, first_row)
                push!(args, INSTANCE(next(ps)))
                return EXPR(BracesCat, args)
            else
                if length(first_row.args) == 1
                    first_row = first_row.args[1]
                else
                    first_row = EXPR(Row, first_row.args)
                end
                ret = EXPR(BracesCat, EXPR[args[1], first_row])
                while kindof(ps.nt) != Tokens.RBRACE
                    if kindof(ps.nt) == Tokens.ENDMARKER
                        break
                    end
                    first_arg = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                    push!(ret, EXPR(Row, EXPR[first_arg]))
                    while kindof(ps.nt) != Tokens.RBRACE && kindof(ps.ws) != NewLineWS && kindof(ps.ws) != SemiColonWS
                        if kindof(ps.nt) == Tokens.ENDMARKER
                            break
                        end
                        a = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                        push!(last(ret.args), a)
                    end
                    # if only one entry dont use :row
                    if length(last(ret.args).args) == 1
                        ret.args[end] = setparent!(ret.args[end].args[1], ret)
                    end
                    update_span!(ret)
                end
                accept_rbrace(ps, ret)
                update_span!(ret)
                return ret
            end
        else
            ret = EXPR(Braces, args)
            push!(ret, first_arg)
            push!(ret, accept_rbrace(ps))
        end
    end
    return ret
end


text:
using CompScienceMeshes
using Test
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)

Sourcechart = Testchart = simplex(pI, pII, pIII)

Accuracy = 12
cf = CommonFace(SauterSchwabQuadrature._legendre(Accuracy,0.0,1.0))

function integrand(x,y)
			return(((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(û,v̂)
	n1 = neighborhood(Testchart, û)
	n2 = neighborhood(Sourcechart, v̂)
	x = cartesian(n1)
	y = cartesian(n2)
	output = integrand(x,y)*jacobian(n1)*jacobian(n2)

return(output)
end

result = sauterschwab_parameterized(INTEGRAND, cf)-
		   verifintegral1(Sourcechart, Testchart, integrand, Accuracy)

@test norm(result) < 1.e-3

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

# Test the use of SauterSchwabQuadrature with the kernel generator utility
# using BEAST
using CompScienceMeshes
using StaticArrays

t1 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.0, -0.980785, -0.19509],
    @SVector[0.0, -0.92388, -0.382683])

rt = RTRefSpace{Float64}()
kernel_cf(x,y) = 1/norm(cartesian(x)-cartesian(y))
igd = generate_integrand_uv(kernel_cf, rt, rt, t1, t1)

t1 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.0, -0.980785, -0.19509],
    @SVector[0.0, -0.92388, -0.382683])

i5 = sauterschwab_parameterized(igd, CommonFace(SauterSchwabQuadrature._legendre(5,0.0,1.0)))
i10 = sauterschwab_parameterized(igd, CommonFace(SauterSchwabQuadrature._legendre(10,0.0,1.0)))

iref = numquad_cf(kernel_cf, rt, rt, t1, t1, zero(i5))

# # BEAST will arbitrate
# tqd = BEAST.quadpoints(rt, [t1], (12,))
# bqd = BEAST.quadpoints(rt, [t1], (13,))
#
# SE_strategy = BEAST.WiltonSEStrategy(
#   tqd[1,1],
#   BEAST.DoubleQuadStrategy(
# 	tqd[1,1],
# 	bqd[1,1],
#   ),
# )
#
# op = BEAST.MWSingleLayer3D(0.0, 4.0π, 0.0)
# z2 = zeros(3,3)
# BEAST.momintegrals!(op, rt, rt, t1, t1, z2, SE_strategy)
#
# @test i5  ≈ iref atol=1e-3
# @test i10 ≈ iref atol=1e-3
# @test i10 ≈ i5   atol=1e-6


text:
__precompile__()

module UnitfulNuclear

  import Unitful
  using Unitful: @unit, @u_str, uconvert

  export eps_0, mu_0
  export @u_str, uconvert

  eps_0 = Unitful.ϵ0
  mu_0 = uconvert(u"T*m/MA", Unitful.μ0)

  @unit n20 "n20" PlasmaDensity 1e20/Unitful.m^3 false
  @unit b "b" Barn 1e-28*Unitful.m^2 true
  @unit u "u" AMU 1.660539040e-27*Unitful.kg false

  const localunits = Unitful.basefactors
  function __init__()
      merge!(Unitful.basefactors, localunits)
      Unitful.register(UnitfulNuclear)
  end

end


text:
module TestAccessors
import PerformanceTestTools
import Accessors

include("test_core.jl")
include("test_optics.jl")
include("test_delete.jl")
include("test_insert.jl")
include("test_examples.jl")
include("test_staticarrays.jl")
include("test_quicktypes.jl")
include("test_setmacro.jl")
include("test_setindex.jl")
include("test_functionlenses.jl")

using Documenter: doctest
if Base.thisminor(VERSION) == v"1.6"
    # ⨟ needs to be defined
    doctest(Accessors)
else
    @info "Skipping doctests, on VERSION = $VERSION"
end
PerformanceTestTools.@include("perf.jl")

end  # module


text:
module LiveSlice

using Interact
using Blink
using CMDimData
using CMDimData.EasyPlot
using MDDatasets

struct SliceParam
	sweep::PSweep
	choosevalue::Observable
	sweepindex::Observable
end

#==Helper functions
===============================================================================#
_throt(x) = throttle(0.1, x) #Default throttling of inputs for this tool.

#Sweep value string
_swvalstr(x::Int) = string(x)
_swvalstr(x::Float64) = EasyPlot.SI(x)

#==Main autoslice() function
===============================================================================#
"""
    autoslice(refreshfn, refparamlist::Vector{PSweep}, slicewfrmlist::Vector{Waveform})

Automatically updates provided waveforms, and calls provided function to refresh plots

#Arguments
 - refreshfn: Function to be called when caller should refresh plots
 - refparamlist: List of parameters to use to slice input data
 - slicewfrmlist: List of DataHR-waveforms to be sliced along different parameter dimensions
"""
function autoslice(refreshfn::Function, refparamlist::Vector{PSweep}, slicewfrmlist::Vector{EasyPlot.Waveform})
	origdatalist = DataHR[wfrm.data for wfrm in slicewfrmlist] #Only supports DataHR
	paramlist = SliceParam[]
	obslist = Observable[] #Index into each parameter sweep value
	widgetlist = [] #To be displayed in Blink

	#Build observables & sliders for each parameter provided by caller:
	for p in refparamlist
		choosevalue = Observable(true)
		sweepindex = Observable(1)
		push!(paramlist, SliceParam(p, choosevalue, sweepindex))
		push!(obslist, choosevalue)
		push!(obslist, _throt(sweepindex)) #Throttle observables to reduce load
		cb = checkbox(choosevalue)
		sld = slider(1:length(p.v), value=sweepindex, label=p.id)
		push!(widgetlist, hbox(cb, sld))
	end

	#Build up an expression string to see which parameter values we are selecting:
	_paramlist_str = map(obslist...) do varlist...
		isfirst = true
		paramstr = ""

		for p in paramlist
			if !isfirst
				paramstr *= " / "
			end
			valstr = "ALL"
			if p.choosevalue[]
				i = p.sweepindex[]
				valstr = _swvalstr(p.sweep.v[i])
			end
			paramstr *= p.sweep.id * "=" * valstr
			isfirst=false
		end
		return HTML(paramstr)
	end
	push!(widgetlist, _paramlist_str)

	#Swap out waveform data with a sliced version (dep. on slider values):
	_updater = map(obslist...) do varlist...
		for (od, wfrm) in zip(origdatalist, slicewfrmlist)
			idxlen = length(size(od.elem))
			Δidxlen = idxlen - length(paramlist)
			idxlist = fill!(Array{Any}(undef, idxlen), :)
			swlist = PSweep[]
			for (i, p) in enumerate(paramlist)
				if p.choosevalue[] #Want a specific index, not all swept values
					idxlist[i] = p.sweepindex[]
				else
					push!(swlist, p.sweep)
				end
			end
			append!(swlist, od.sweeps[(end+1-Δidxlen):end])
			newelem = od.elem[idxlist...]
			wfrm.data = DataHR(swlist, newelem)
		end
		refreshfn() #Ask user to refresh plots, etc
	end

	#Display main window:
	wnd = Window()
	body!(wnd, vbox()) #first call to body! does not always work.
	body!(wnd, hbox(vbox(widgetlist...), nothing))
	return
end

end #module


text:
### Generic MCBaseSample type used by most samplers

type MCBaseSample <: MCSample{NullOrder}
  sample::Vector{Float64}
  logtarget::Float64
end

MCBaseSample(s::Vector{Float64}) = MCBaseSample(s, NaN)
MCBaseSample() = MCBaseSample(Float64[], NaN)

MCBaseSample(l::Int) = MCBaseSample(fill(NaN, l), NaN)

logtarget!(s::MCSample, f::Function) = (s.logtarget = f(s.sample))

### MCGradSample is used by samplers that compute (up to) the gradient of the log-target (for ex HMC and MALA)

type MCGradSample <: MCSample{FirstOrder}
  sample::Vector{Float64}
  logtarget::Float64
  gradlogtarget::Vector{Float64}
end

MCGradSample(s::Vector{Float64}) = MCGradSample(s, NaN, Float64[])
MCGradSample() = MCGradSample(Float64[], NaN, Float64[])

MCGradSample(l::Int) = MCGradSample(fill(NaN, l), NaN, fill(NaN, l))

gradlogtargetall!(s::MCSample, f::Function) = ((s.logtarget, s.gradlogtarget) = f(s.sample))

### MCMCTensorSample is used by samplers that compute (up to) the tensor of the log-target (for ex SMMALA)

type MCTensorSample <: MCSample{SecondOrder}
  sample::Vector{Float64}
  logtarget::Float64
  gradlogtarget::Vector{Float64}
  tensorlogtarget::Matrix{Float64}
end

MCTensorSample(s::Vector{Float64}) = MCTensorSample(s, NaN, Float64[], Array(Float64, 0, 0))
MCTensorSample() = MCTensorSample(Float64[], NaN, Float64[], Array(Float64, 0, 0))

MCTensorSample(l::Int) = MCTensorSample(fill(NaN, l), NaN, fill(NaN, l), fill(NaN, l, l))

tensorlogtargetall!(s::MCSample, f::Function) = ((s.logtarget, s.gradlogtarget, s.tensorlogtarget) = f(s.sample))

### MCMCDTensorSample is used by samplers that compute (up to) the derivative of the tensor of the log-target (for ex
### RMHMC and PMALA)

type MCDTensorSample <: MCSample{ThirdOrder}
  sample::Vector{Float64}
  logtarget::Float64
  gradlogtarget::Vector{Float64}
  tensorlogtarget::Matrix{Float64}
  dtensorlogtarget::Array{Float64, 3}
end

MCDTensorSample(s::Vector{Float64}) = MCDTensorSample(s, NaN, Float64[], Array(Float64, 0, 0), Array(Float64, 0, 0, 0))
MCDTensorSample() = MCDTensorSample(Float64[], NaN, Float64[], Array(Float64, 0, 0), Array(Float64, 0, 0, 0))

MCDTensorSample(l::Int) = MCDTensorSample(fill(NaN, l), NaN, fill(NaN, l), fill(NaN, l, l), fill(NaN, l, l, l))

dtensorlogtargetall!(s::MCSample, f::Function) =
  ((s.logtarget, s.gradlogtarget, s.tensorlogtarget, s.dtensorlogtarget) = f(s.sample))


text:
using Test
using Printf: @printf
using LinearAlgebra: tr
using StaticArrays: SVector
using QuantumLattices.Essentials.Frameworks
using QuantumLattices.Essentials: update!, reset!, register!
using QuantumLattices.Essentials.Spatials: Point, PID, Bond, Bonds, Lattice, acrossbonds, zerothbonds
using QuantumLattices.Essentials.DegreesOfFreedom: SimpleIID, SimpleInternal, IIDSpace, Coupling, Subscript, Subscripts, SubscriptsID
using QuantumLattices.Essentials.DegreesOfFreedom: Term, Hilbert, Index, Table, OID, OIDToTuple, Operator, Operators, plain, @couplings
using QuantumLattices.Essentials.QuantumOperators: ID, id, idtype, Identity
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Interfaces:  expand!, expand, add!
using QuantumLattices.Prerequisites.Traits: contentnames
using QuantumLattices.Prerequisites.CompositeStructures: NamedContainer

import QuantumLattices.Essentials.Frameworks: Parameters
import QuantumLattices.Essentials: prepare!, run!, update!
import QuantumLattices.Essentials.DegreesOfFreedom: ishermitian, couplingcenters
import QuantumLattices.Prerequisites.VectorSpaces: shape, ndimshape

struct FID{N<:Union{Int, Symbol}} <: SimpleIID
    nambu::N
end
@inline Base.adjoint(sl::FID{Int}) = FID(3-sl.nambu)

struct FFock <: SimpleInternal{FID{Int}}
    nnambu::Int
end
@inline shape(f::FFock) = (1:f.nnambu,)
@inline ndimshape(::Type{FFock}) = 1
@inline FID(i::CartesianIndex, vs::FFock) = FID(i.I...)
@inline CartesianIndex(did::FID{Int}, vs::FFock) = CartesianIndex(did.nambu)
@inline shape(iidspace::IIDSpace{FID{Symbol}, FFock}) = (1:iidspace.internal.nnambu,)
@inline shape(iidspace::IIDSpace{FID{Int}, FFock}) = (iidspace.iid.nambu:iidspace.iid.nambu,)
@inline shape(iidspace::IIDSpace{FID{Symbol}, FFock}) = (1:iidspace.internal.nnambu,)
@inline shape(iidspace::IIDSpace{FID{Int}, FFock}) = (iidspace.iid.nambu:iidspace.iid.nambu,)

const FCoupling{V, I<:ID{FID}, C<:Subscripts, CI<:SubscriptsID} = Coupling{V, I, C, CI}
@inline couplingcenters(::(Coupling{V, <:ID{FID}} where V), ::Bond, ::Val) = (1, 2)
@inline FCoupling(value, nambus::Tuple{Vararg{Int}}) = Coupling(value, ID(FID, nambus), Subscripts((nambu=Subscript(nambus),)))
@inline FCoupling(value, nambus::Subscript) = Coupling(value, ID(FID, convert(Tuple, nambus)), Subscripts((nambu=nambus,)))

@inline ishermitian(::Type{<:Term{:Mu}}) = true
@inline ishermitian(::Type{<:Term{:Hp}}) = false

@testset "Parameters" begin
    ps1 = Parameters{(:t1, :t2, :U)}(1.0im, 1.0, 2.0)
    ps2 = Parameters{(:t1, :U)}(1.0im, 2.0)
    ps3 = Parameters{(:t1, :U)}(1.0im, 2.1)
    @test match(ps1, ps2) == true
    @test match(ps1, ps3) == false
end

@testset "Generator" begin
    @test contentnames(AbstractGenerator) == (:operators, :table, :boundary)
    @test contentnames(Generator) == (:operators, :terms, :bonds, :hilbert, :half, :table, :boundary)
    @test contentnames(SimplifiedGenerator) == (:parameters, :operators, :table, :boundary)

    lattice = Lattice("Tuanzi", [Point(PID(1), (0.0, 0.0), (0.0, 0.0)), Point(PID(2), (0.5, 0.0), (0.0, 0.0))], vectors=[[1.0, 0.0]], neighbors=1)
    bonds = Bonds(lattice)
    hilbert = Hilbert{FFock}(pid->FFock(2), lattice.pids)
    table = Table(hilbert, OIDToTuple(:scope, :site))
    t = Term{:Hp}(:t, 2.0, 1, couplings=@couplings(FCoupling(1.0, (2, 1))))
    μ = Term{:Mu}(:μ, 1.0, 0, couplings=@couplings(FCoupling(1.0, (2, 1))), modulate=true)
    tops₁ = expand(t, filter(acrossbonds, bonds, Val(:exclude)), hilbert, half=true, table=table)
    tops₂ = expand(one(t), filter(acrossbonds, bonds, Val(:include)), hilbert, half=true, table=table)
    μops = expand(one(μ), filter(zerothbonds, bonds, Val(:include)), hilbert, half=true, table=table)
    i = Identity()

    optp = Operator{Float, ID{OID{Index{PID, FID{Int}}, SVector{2, Float}}, 2}}
    entry = Entry(tops₁, NamedContainer{(:μ,)}((μops,)), NamedContainer{(:t, :μ)}((tops₂, Operators{optp|>idtype, optp}())))
    @test entry == deepcopy(entry) && isequal(entry, deepcopy(entry))
    @test entry == Entry((t, μ), bonds, hilbert, half=true, table=table)
    @test entry|>eltype == entry|>typeof|>eltype == Operators{idtype(optp), optp}
    @test expand!(Operators{idtype(optp), optp}(), entry, plain, t=2.0, μ=1.5) == tops₁+tops₂*2.0+μops*1.5
    @test empty(entry) == empty!(deepcopy(entry))
    @test empty(entry) == Entry(empty(μops), NamedContainer{(:μ,)}((empty(μops),)), NamedContainer{(:t, :μ)}((empty(μops), empty(μops))))
    @test merge!(empty(entry), entry) == merge(entry, entry) == entry
    @test reset!(deepcopy(entry), (t, μ), bonds, hilbert, half=true, table=table) == entry
    @test i(entry) == entry

    cgen = Generator((t, μ), bonds, hilbert; half=true, table=table, boundary=plain)
    @test cgen == deepcopy(cgen) && isequal(cgen, deepcopy(cgen))
    @test Parameters(cgen) == Parameters{(:t, :μ)}(2.0, 1.0)
    @test expand!(Operators{idtype(optp), optp}(), cgen) == expand(cgen) == tops₁+tops₂*2.0+μops
    @test expand(cgen, :t) == tops₁+tops₂*2.0
    @test expand(cgen, :μ) == μops
    @test expand(cgen, 1)+expand(cgen, 2)+expand(cgen, 3)+expand(cgen, 4) == expand(cgen)
    @test expand(cgen, :μ, 1)+expand(cgen, :μ, 2) == μops
    @test expand(cgen, :t, 3) == tops₁
    @test expand(cgen, :t, 4) == tops₂*2.0
    @test empty!(deepcopy(cgen)) == Generator((t, μ), empty(bonds), empty(hilbert), half=true, table=empty(table), boundary=plain) == empty(cgen)
    @test reset!(empty(cgen), lattice) == cgen
    @test update!(cgen, μ=1.5)|>expand == tops₁+tops₂*2.0+μops*1.5

    params = Parameters{(:t, :μ)}(2.0, 1.0)
    sgen = SimplifiedGenerator(params, entry, table=table, boundary=plain)
    @test Parameters(sgen) == params
    @test expand!(Operators{idtype(optp), optp}(), sgen) == expand(sgen) == tops₁+tops₂*2.0+μops
    @test empty!(deepcopy(sgen)) == SimplifiedGenerator(params, empty(entry), table=empty(table), boundary=plain) == empty(sgen)
    @test reset!(empty(sgen), entry, table=table) == sgen
    @test update!(sgen, μ=1.5)|>expand == tops₁+tops₂*2.0+μops*1.5
    @test i(cgen) == sgen
end

mutable struct VCA <: Engine
    t::Float
    U::Float
    dim::Int
end
function update!(vca::VCA; kwargs...)
    vca.t = get(kwargs, :t, vca.t)
    vca.U = get(kwargs, :U, vca.U)
    return vca
end
@inline Parameters(vca::VCA) = Parameters{(:t, :U)}(vca.t, vca.U)

mutable struct GF <: Action
    count::Int
end
@inline prepare!(gf::GF, vca::VCA) = Matrix{Float}(undef, vca.dim, vca.dim)
@inline run!(alg::Algorithm{VCA}, assign::Assignment{GF}) = (assign.action.count += 1; assign.data[:, :] .= alg.engine.t+alg.engine.U)

mutable struct DOS <: Action
    mu::Float
end
@inline update!(eb::DOS; kwargs...) = (eb.mu = get(kwargs, :mu, eb.mu); eb)
@inline prepare!(dos::DOS, vca::VCA) = 0.0
@inline function run!(alg::Algorithm{VCA}, assign::Assignment{DOS})
    rundependences!(alg, assign)
    gf = alg.assignments[first(assign.dependences)]
    assign.data = tr(gf.data)
end
@inline dosmap(params::Parameters) = Parameters{(:t, :U, :mu)}(params.t, params.U, -params.U/2)

@testset "Framework" begin
    vca = VCA(1.0, 8.0, 4)
    @test vca == deepcopy(vca)
    @test isequal(vca, deepcopy(vca))
    @test repr(vca) == "VCA"
    @test string(vca) == "VCA"

    gf = GF(0)
    @test gf == deepcopy(gf)
    @test isequal(gf, deepcopy(gf))
    @test update!(gf) == gf

    vca = Algorithm("test", vca)
    @test vca == deepcopy(vca)
    @test isequal(vca, deepcopy(vca))
    @test repr(vca, ≠(:U)) == "test(VCA)_1.0"
    @test string(vca) == repr(vca) == "test(VCA)_1.0_8.0"

    add!(vca, :GF, GF(0))
    rex = r"Action DOS\(DOS\)\: time consumed [0-9]*\.[0-9]*(e[+-][0-9]*)*s."
    @test_logs (:info, rex) register!(vca, :DOS, DOS(-3.5), parameters=(U=7.0,), map=dosmap, dependences=(:GF,))

    dos = vca.assignments[:DOS]
    @test dos == deepcopy(dos)
    @test isequal(dos, deepcopy(dos))
    @test valtype(dos) == valtype(typeof(dos)) == Float
    @test dos.data == 32.0
    @test dos.action.mu == -3.5
    @test nameof(vca, dos) == "test(VCA)_1.0_7.0_DOS"

    update!(dos, U=6.0)
    run!(vca, :DOS, false)
    @test dos.parameters == (t=1.0, U=6.0)
    @test dos.data == 28.0
    @test dos.action.mu == -3.0
end


text:
# Kernel Principal Component Analysis

using Arpack: eigs

"""Center a kernel matrix"""
struct KernelCenter{T<:Real}
    means::AbstractVector{T}
    total::T
end

"""Fit `KernelCenter` object"""
function fit(::Type{KernelCenter}, K::AbstractMatrix{<:Real})
    n = size(K, 1)
    means = vec(mean(K, dims=2))
    KernelCenter(means, sum(means) / n)
end

"""Center kernel matrix."""
function transform!(C::KernelCenter, K::AbstractMatrix{<:Real})
    r, c = size(K)
    tot = C.total
    means = mean(K, dims=1)
    @simd for i in 1:r
        for j in 1:c
            @inbounds K[i, j] -= C.means[i] + means[j] - tot
        end
    end
    return K
end

"""
This type contains kernel PCA model parameters.
"""
struct KernelPCA{T<:Real} <: NonlinearDimensionalityReduction
    X::AbstractMatrix{T}           # fitted data or precomputed kernel
    ker::Union{Nothing, Function}  # kernel function
    center::KernelCenter           # kernel center
    λ::AbstractVector{T}           # eigenvalues  in feature space
    α::AbstractMatrix{T}           # eigenvectors in feature space
    inv::AbstractMatrix{T}         # inverse transform coefficients
end

## properties
"""
    size(M::KernelPCA)

Returns a tuple with the input dimension ``d``, *i.e* the dimension of the observation space, and the output dimension ``p``, *i.e* the dimension of the principal subspace.
"""
size(M::KernelPCA) = (size(M.X, 1), length(M.λ))

"""
    eigvals(M::KernelPCA)

Return eigenvalues of the kernel matrix of the model `M`.
"""
eigvals(M::KernelPCA) = M.λ

"""
    eigvecs(M::KernelPCA)

Return eigenvectors of the kernel matrix of the model `M`.
"""
eigvecs(M::KernelPCA) = M.α

"""
    projection(M::KernelPCA)

Return the projection matrix (of size ``n \\times p``).
Each column of the projection matrix corresponds to an eigenvector, and ``n`` is a number of observations.

The principal components are arranged in descending order of the corresponding eigenvalues.
"""
projection(M::KernelPCA) = eigvecs(M) ./ sqrt.(eigvals(M)')

## use

"""
    predict(M::KernelPCA, x)

Transform out-of-sample transformation of `x` into a kernel space of the model `M`.

Here, `x` can be either a vector of length `d` or a matrix where each column is an observation.
"""
function predict(M::KernelPCA, x::AbstractVecOrMat{T}) where {T<:Real}
    k = pairwise(M.ker, eachcol(M.X), eachcol(x))
    transform!(M.center, k)
    return projection(M)'*k
end

"""
    predict(M::KernelPCA)

Transform the data fitted to the model `M` to a kernel space of the model.
"""
predict(M::KernelPCA) = sqrt.(eigvals(M)) .* eigvecs(M)'

"""
    reconstruct(M::KernelPCA, y)

Approximately reconstruct observations, given in `y`, to the original space using the kernel PCA model `M`.

Here, `y` can be either a vector of length `p` or a matrix where each column gives the principal components for an observation.
"""
function reconstruct(M::KernelPCA, y::AbstractVecOrMat{T}) where {T<:Real}
    if size(M.inv, 1) == 0
        throw(ArgumentError("Inverse transformation coefficients are not available, set `inverse` parameter when fitting data"))
    end
    Pᵗ = predict(M)
    k = pairwise(M.ker, eachcol(Pᵗ), eachcol(y))
    return M.inv*k
end

## show

function show(io::IO, M::KernelPCA)
    indim, outdim = size(M)
    print(io, "Kernel PCA(indim = $indim, outdim = $outdim)")
end

## interface functions
"""
    fit(KernelPCA, X; ...)

Perform kernel PCA over the data given in a matrix `X`. Each column of `X` is an observation.

This method returns an instance of [`KernelPCA`](@ref).

**Keyword arguments:**

Let `(d, n) = size(X)` be respectively the input dimension and the number of observations:

- `kernel`: The kernel function. This functions accepts two vector arguments `x` and `y`,
and returns a scalar value (*default:* `(x,y)->x'y`)
- `solver`: The choice of solver:
    - `:eig`: uses `LinearAlgebra.eigen` (*default*)
    - `:eigs`: uses `Arpack.eigs` (always used for sparse data)
- `maxoutdim`:  Maximum output dimension (*default* `min(d, n)`)
- `inverse`: Whether to perform calculation for inverse transform for non-precomputed kernels (*default* `false`)
- `β`: Hyperparameter of the ridge regression that learns the inverse transform (*default* `1` when `inverse` is `true`).
- `tol`: Convergence tolerance for `eigs` solver (*default* `0.0`)
- `maxiter`: Maximum number of iterations for `eigs` solver (*default* `300`)
"""
function fit(::Type{KernelPCA}, X::AbstractMatrix{T};
             kernel::Union{Nothing, Function} = (x,y)->x'y,
             maxoutdim::Int = min(size(X)...),
             remove_zero_eig::Bool = false, atol::Real = 1e-10,
             solver::Symbol = :eig,
             inverse::Bool = false,  β::Real = convert(T, 1.0),
             tol::Real = 0.0, maxiter::Real = 300) where {T<:Real}

    d, n = size(X)
    maxoutdim = min(min(d, n), maxoutdim)

    # set kernel function if available
    K = if isa(kernel, Function)
        pairwise(kernel, eachcol(X), symmetric=true)
    elseif kernel === nothing
        @assert issymmetric(X) "Precomputed kernel matrix must be symmetric."
        inverse = false
        X
    else
        throw(ArgumentError("Incorrect kernel type. Use a function or a precomputed kernel."))
    end

    # center kernel
    center = fit(KernelCenter, K)
    transform!(center, K)

    # perform eigenvalue decomposition
    evl, evc = if solver == :eigs || SparseArrays.issparse(K)
        evl, evc = eigs(K, nev=maxoutdim, which=:LR, v0=2.0*rand(n) .- 1.0, tol=tol, maxiter=maxiter)
        real.(evl), real.(evc)
    else
        Eg = eigen(Hermitian(K))
        Eg.values, Eg.vectors
    end

    # sort eigenvalues in descending order
    ord = sortperm(evl; rev=true)
    ord = ord[1:min(length(ord), maxoutdim)]

    # remove zero eigenvalues
    λ, α = if remove_zero_eig
        ez = map(!, isapprox.(evl[ord], zero(T), atol=atol))
        evl[ord[ez]], evc[:, ord[ez]]
    else
        evl[ord], evc[:, ord]
    end

    # calculate inverse transform coefficients
    Q = zeros(T, 0, 0)
    if inverse
        Pᵗ = α' .* sqrt.(λ)
        KT = pairwise(kernel, eachcol(Pᵗ), symmetric=true)
        Q = (KT + diagm(0 => fill(β, size(KT,1)))) \ X'
    end

    KernelPCA(X, kernel, center, λ, α, Q')
end


text:
using QMSimFiles
using Test

snpmap1 = "lm_mrk_001.txt"
qtlmap1 = "lm_qtl_001.txt"
qtleff1 = "effect_qtl_001.txt"
snpdata1 = "p1_mrk_001.txt"
snpdata2 = "p1_mrk_002.txt"   # saved as /snp_code
qtldata1 = "p1_qtl_001.txt"

@testset "get_number_of_chromosomes" begin
   nchr = QMSimFiles.get_number_of_chromosomes(snpmap1)
   @test nchr == 3
   nchr = QMSimFiles.get_number_of_chromosomes(qtlmap1)
   @test nchr == 3
end

@testset "get_number_of_loci" begin
   nchr = QMSimFiles.get_number_of_chromosomes(snpmap1)
   nloci = QMSimFiles.get_number_of_loci(snpmap1,nchr)
   @test all( nloci == [4,3,4] )
   nchr = QMSimFiles.get_number_of_chromosomes(qtlmap1)
   nloci = QMSimFiles.get_number_of_loci(qtlmap1,nchr)
   @test all( nloci == [3,4,2] )
end

@testset "get_number_of_QTL_allele" begin
   maxna,maxAllele = QMSimFiles.get_number_of_QTL_allele(qtleff1,3)
   @test maxna == 4
   @test all( maxAllele .== [4,4,3] )
end

@testset "load_SNP_QTL_maps!" begin
   maxna,maxAllele = QMSimFiles.get_number_of_QTL_allele(qtleff1,3)
   maps = Vector{QMSimFiles.QMSimChromosomeMap}(undef,3)
   maps[1] = QMSimFiles.QMSimChromosomeMap(7,4,3,zeros(Int,7),zeros(Float64,7),maxAllele[1],zeros(Int,maxAllele[1]),zeros(Float64,maxAllele[1],3))
   maps[2] = QMSimFiles.QMSimChromosomeMap(7,3,4,zeros(Int,7),zeros(Float64,7),maxAllele[2],zeros(Int,maxAllele[2]),zeros(Float64,maxAllele[2],4))
   maps[3] = QMSimFiles.QMSimChromosomeMap(6,4,2,zeros(Int,6),zeros(Float64,6),maxAllele[3],zeros(Int,maxAllele[3]),zeros(Float64,maxAllele[3],2))
   QMSimFiles.load_SNP_QTL_maps!(snpmap1,qtlmap1,maps)
   @test all( maps[1].seqQTL .== [0,0,0,0,1,2,3] )
   @test all( maps[2].seqQTL .== [0,0,0,4,5,6,7] )
   @test all( maps[3].seqQTL .== [0,0,0,0,8,9] )
   @test all( maps[1].pos .≈ [5.99360,46.98708,71.46866,78.57678, 23.76397,28.15399,89.81130] )
   @test all( maps[2].pos .≈ [8.92501,13.56385,43.67896, 29.88287,38.45590,66.06151,82.55315] )
   @test all( maps[3].pos .≈ [25.60614,86.19237,91.47393,98.69095, 14.34842,24.60904] )
   omaps = deepcopy(maps)

   QMSimFiles.sort_by_position!(maps[1].pos, maps[1].seqQTL)
   QMSimFiles.sort_by_position!(maps[2].pos, maps[2].seqQTL)
   QMSimFiles.sort_by_position!(maps[3].pos, maps[3].seqQTL)
   perm1 = [1,5,6,2,3,4,7]
   perm2 = [1,2,4,5,3,6,7]
   perm3 = [5,6,1,2,3,4]
   @test all( maps[1].seqQTL .== omaps[1].seqQTL[perm1] )
   @test all( maps[2].seqQTL .== omaps[2].seqQTL[perm2] )
   @test all( maps[3].seqQTL .== omaps[3].seqQTL[perm3] )
   @test all( maps[1].pos .≈ omaps[1].pos[perm1] )
   @test all( maps[2].pos .≈ omaps[2].pos[perm2] )
   @test all( maps[3].pos .≈ omaps[3].pos[perm3] )
end

@testset "load_QTL_effect! and read_maps" begin
   maxna,maxAllele = QMSimFiles.get_number_of_QTL_allele(qtleff1,3)
   maps = Vector{QMSimFiles.QMSimChromosomeMap}(undef,3)
   maps[1] = QMSimFiles.QMSimChromosomeMap(7,4,3,zeros(Int,7),zeros(Float64,7),maxAllele[1],zeros(Int,3),zeros(Float64,maxAllele[1],3))
   maps[2] = QMSimFiles.QMSimChromosomeMap(7,3,4,zeros(Int,7),zeros(Float64,7),maxAllele[2],zeros(Int,4),zeros(Float64,maxAllele[2],4))
   maps[3] = QMSimFiles.QMSimChromosomeMap(6,4,2,zeros(Int,6),zeros(Float64,6),maxAllele[3],zeros(Int,2),zeros(Float64,maxAllele[3],2))
   QMSimFiles.load_SNP_QTL_maps!(snpmap1,qtlmap1,maps)
   QMSimFiles.sort_by_position!(maps[1].pos, maps[1].seqQTL)
   QMSimFiles.sort_by_position!(maps[2].pos, maps[2].seqQTL)
   QMSimFiles.sort_by_position!(maps[3].pos, maps[3].seqQTL)

   QMSimFiles.load_QTL_effect!(qtleff1,maps)
   @test all( maps[1].naQTL .== [3,4,4] )
   @test all( maps[2].naQTL .== [3,4,2,2] )
   @test all( maps[3].naQTL .== [3,3] )
   @test all( maps[1].effQTL[1:3,1] .≈ [-0.177111,-0.223061,0.071867] )
   @test all( maps[1].effQTL[1:4,2] .≈ [0.034120,-0.053426,-0.032978,0.022941] )
   @test all( maps[1].effQTL[1:4,3] .≈ [-0.010957,-0.001999,0.033815,-0.010451] )
   @test all( maps[2].effQTL[1:3,1] .≈ [0.054960,-0.038202,-0.063630] )
   @test all( maps[2].effQTL[1:4,2] .≈ [0.130199,0.137865,-0.455851,0.296964] )
   @test all( maps[2].effQTL[1:2,3] .≈ [0.001240,-0.000023] )
   @test all( maps[2].effQTL[1:2,4] .≈ [-0.071812,0.063439] )
   @test all( maps[3].effQTL[1:3,1] .≈ [-0.020500,0.004283,-0.075465] )
   @test all( maps[3].effQTL[1:3,2] .≈ [0.004419,-0.010919,0.021111] )

   refmap = read_maps(snpmap1,qtlmap1,qtleff1)
   @test refmap.nchr == 3
   @test refmap.totalSNP == 11
   @test refmap.totalQTL == 9
   for i in 1:3
      @test refmap.chr[i].nLoci == maps[i].nLoci
      @test refmap.chr[i].nSNP == maps[i].nSNP
      @test refmap.chr[i].nQTL == maps[i].nQTL
      @test all( refmap.chr[i].seqQTL .== maps[i].seqQTL )
      @test all( refmap.chr[i].pos .≈ maps[i].pos )
      @test refmap.chr[i].maxAllele .== maps[i].maxAllele
      @test all( refmap.chr[i].naQTL .== maps[i].naQTL )
      @test all( refmap.chr[i].effQTL .≈ maps[i].effQTL )
   end
end

@testset "generate_chromosome_set" begin
   gmap = read_maps(snpmap1,qtlmap1,qtleff1)
   chromosome_set = QMSimFiles.generate_chromosome_set(gmap)
   @test length(chromosome_set)==gmap.nchr
end

@testset "parse SNP/QTL data" begin
   gmap = read_maps(snpmap1,qtlmap1,qtleff1)
   @test QMSimFiles.is_snpcode_file(snpdata1,gmap.totalSNP)==false
   @test QMSimFiles.is_snpcode_file(snpdata2,gmap.totalSNP)==true

   # for SNP (standard)
   stext = "2       2 1 2 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1"
   snp = zeros(Int8,gmap.totalSNP*2)
   QMSimFiles.text_to_code!(stext,gmap.totalSNP,snp,7)
   @test all( snp .== [2, 1, 2, 2, 1, 1, 1, 1,   1, 1, 1, 1, 1, 1,   1, 2, 1, 1, 1, 2, 1, 1] )

   chromosome_set1 = QMSimFiles.generate_chromosome_set(gmap)
   QMSimFiles.convert_snpdata_to_haplotype!(snp,gmap,chromosome_set1)
   @test all( chromosome_set1[1].gp[findall(x->gmap.chr[1].seqQTL[x]==0,1:gmap.chr[1].nLoci)] .== [1,1,0,0] )
   @test all( chromosome_set1[1].gm[findall(x->gmap.chr[1].seqQTL[x]==0,1:gmap.chr[1].nLoci)] .== [0,1,0,0] )
   @test all( chromosome_set1[2].gp[findall(x->gmap.chr[2].seqQTL[x]==0,1:gmap.chr[2].nLoci)] .== [0,0,0] )
   @test all( chromosome_set1[2].gm[findall(x->gmap.chr[2].seqQTL[x]==0,1:gmap.chr[2].nLoci)] .== [0,0,0] )
   @test all( chromosome_set1[3].gp[findall(x->gmap.chr[3].seqQTL[x]==0,1:gmap.chr[3].nLoci)] .== [0,0,0,0] )
   @test all( chromosome_set1[3].gm[findall(x->gmap.chr[3].seqQTL[x]==0,1:gmap.chr[3].nLoci)] .== [1,0,1,0] )

   # for SNP (snp_code)
   ctext = "2      42000003030"
   snp_code = zeros(Int8,gmap.totalSNP)
   QMSimFiles.text_to_snpcode!(ctext,gmap.totalSNP,snp_code,7)
   @test all( snp_code .== [4,2,0,0, 0,0,0, 3,0,3,0] )

   chromosome_set2 = QMSimFiles.generate_chromosome_set(gmap)
   QMSimFiles.convert_snpcode_to_haplotype!(snp_code,gmap,chromosome_set2)
   @test all( chromosome_set2[1].gp .== chromosome_set1[1].gp )
   @test all( chromosome_set2[1].gm .== chromosome_set1[1].gm )
   @test all( chromosome_set2[2].gp .== chromosome_set1[2].gp )
   @test all( chromosome_set2[2].gm .== chromosome_set1[2].gm )
   @test all( chromosome_set2[3].gp .== chromosome_set1[3].gp )
   @test all( chromosome_set2[3].gm .== chromosome_set1[3].gm )

   # for QTL
   qtext = "2       2 2 2 4 4 1 1 1 2 1 2 2 2 2 2 1 1 2"
   qtl = zeros(Int8,gmap.totalQTL*2)
   QMSimFiles.text_to_code!(qtext,gmap.totalQTL,qtl,7)
   @test all( qtl .== [2, 2, 2, 4, 4, 1,   1, 1, 2, 1, 2, 2, 2, 2,   2, 1, 1, 2] )

   QMSimFiles.convert_qtldata_to_haplotype!(qtl,gmap,chromosome_set2)
   @test all( chromosome_set2[1].gp[findall(x->gmap.chr[1].seqQTL[x]>0,1:gmap.chr[1].nLoci)] .== [2,2,4] )
   @test all( chromosome_set2[1].gm[findall(x->gmap.chr[1].seqQTL[x]>0,1:gmap.chr[1].nLoci)] .== [2,4,1] )
   @test all( chromosome_set2[2].gp[findall(x->gmap.chr[2].seqQTL[x]>0,1:gmap.chr[2].nLoci)] .== [1,2,2,2] )
   @test all( chromosome_set2[2].gm[findall(x->gmap.chr[2].seqQTL[x]>0,1:gmap.chr[2].nLoci)] .== [1,1,2,2] )
   @test all( chromosome_set2[3].gp[findall(x->gmap.chr[3].seqQTL[x]>0,1:gmap.chr[3].nLoci)] .== [2,1] )
   @test all( chromosome_set2[3].gm[findall(x->gmap.chr[3].seqQTL[x]>0,1:gmap.chr[3].nLoci)] .== [1,2] )

   # general function
   g1 = read_genotypes(snpdata1,qtldata1,gmap)
   g2 = read_genotypes(snpdata2,qtldata1,gmap)
   for i=1:3
      @test all( g1[2].chr[i].gp .== chromosome_set2[i].gp )
      @test all( g1[2].chr[i].gm .== chromosome_set2[i].gm )
      @test all( g2[2].chr[i].gp .== chromosome_set2[i].gp )
      @test all( g2[2].chr[i].gm .== chromosome_set2[i].gm )
   end
   @test isapprox(g1[1].tbv, -0.412531, atol=1e-5)
   @test isapprox(g1[2].tbv, -0.015914, atol=1e-5)
   @test isapprox(g1[3].tbv, -0.137409, atol=1e-5)
   @test isapprox(g1[4].tbv, -0.242572, atol=1e-5)
   for i=1:4
      @test g1[i].tbv ≈ g2[i].tbv
   end
end

@testset "number_of_crosses" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   chr = g.map.chr[1]
   refpos = [5.9936, 23.76397, 28.15399, 46.98708, 71.46866, 78.57678, 89.8113]
   n = QMSimFiles.number_of_crosses(chr.pos[end])
   @test (0 <= n) && (n <= 100)
end

@testset "location_of_crossover" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   chr = g.map.chr[1]
   ntests = 10
   for i=1:ntests
      chrlen = chr.pos[end]
      ncross = QMSimFiles.number_of_crosses(chrlen) + 1
      (ncross,gameteid,loccross) = QMSimFiles.location_of_crossover(chrlen,chr.pos,ncross)
      if ncross>0
         @test issorted(loccross) && (loccross[1]>=0) && (loccross[end]<=chrlen) && length(loccross)==ncross
         @test all( map(x -> ((gameteid[x][1] in [1,2]) && gameteid[x][2] in [3,4]), 1:ncross ) ) && length(gameteid)==ncross
      end
   end
   # adjustment to remove double cross etc
   gameteid = [(1,3),(1,3),(1,4),(1,4),(2,3),(2,3),(2,4),(2,4)]
   loccross = [1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0]
   (ncross2,gameteid2,loccross2) = QMSimFiles.adjusted_location_of_crossover(gameteid,loccross)
   @test ncross2==0

   gameteid = [(1,3),(1,3),(1,4),(1,4),(1,4),(1,4),(2,4),(2,4)]
   loccross = [1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0]
   (ncross2,gameteid2,loccross2) = QMSimFiles.adjusted_location_of_crossover(gameteid,loccross)
   @test ncross2==0

   gameteid = [(2,3),(1,3),(1,4),(1,3),(2,3),(1,3),(2,4),(1,4),(2,4)]
   loccross = [1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0]
   (ncross2,gameteid2,loccross2) = QMSimFiles.adjusted_location_of_crossover(gameteid,loccross)
   @test ncross2==1

   gameteid = [(2,3),(1,3),(1,4),(1,3),(2,3),(1,3),(2,4),(1,4),(2,4),(2,3),(1,4)]
   loccross = [1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0,  10.0, 11.0]
   (ncross2,gameteid2,loccross2) = QMSimFiles.adjusted_location_of_crossover(gameteid,loccross)
   @test ncross2==3
end

@testset "generate_gamete_with_crossover" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   chr = g.map.chr[1]
   chrlen = chr.pos[end]
   refpos = [5.9936, 23.76397, 28.15399, 46.98708, 71.46866, 78.57678, 89.8113]

   # refgp = [1, 3, 3, 0, 0, 0, 4]
   # refgm = [1, 3, 4, 0, 0, 1, 4]
   gp = g.individual[4].chr[1].gp
   gm = g.individual[4].chr[1].gm
   ncross = 1
   gameteid = [(1,3)]
   loccross = [3]
   ga = QMSimFiles.generate_gamete_with_crossover(gameteid,loccross,1,gp,gm)
   @test all( ga .== [1,3,4,0,0,1,4] )

   ncross = 2
   gameteid = [(1,3),(1,4)]
   loccross = [3,6]
   ga = QMSimFiles.generate_gamete_with_crossover(gameteid,loccross,1,gp,gm)
   @test all( ga .== [1,3,4,0,0,1,4] )
end

@testset "general mating" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   sumtbv = 0.0
   nanim = 1000
   for i=1:nanim
      animal = mating(g.map, g.individual[1], g.individual[2])
      sumtbv = sumtbv + animal.tbv
   end
   pa = (g.individual[1].tbv + g.individual[2].tbv)/2
   @test isapprox(sumtbv/nanim,pa, atol=1e-1)
end

@testset "equal sign" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   @test g.map==g.map
   @test g.individual[1]==g.individual[1]
   @test !(g.individual[1]==g.individual[2])
end

@testset "pack/unpack genotypes" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   for k=1:length(g.individual)
      # pack
      ref_packed = zeros(Int8,0)
      for i=1:g.map.nchr
         ref_packed = [ref_packed; g.individual[k].chr[i].gp]
         ref_packed = [ref_packed; g.individual[k].chr[i].gm]
      end
      ind_packed = QMSimFiles.pack_genome(g.map,g.individual[k])
      @test all( ind_packed .== ref_packed )
      packed = QMSimFiles.pack_genome(g.map,g.individual[k].chr)
      @test all( packed .== ref_packed )

      # unpack
      chromosome_set = QMSimFiles.unpack_genome(g.map,packed)
      @test all( chromosome_set == g.individual[k].chr )
   end
end

@testset "save/read genotypes" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)

   #hdf5file = "runtests.h5"
   hdf5file = tempname()
   @info "temporary HDF5 file: $(hdf5file)"
   save_qmsim_data_hdf5(g,hdf5file)

   # map
   saved_map = read_qmsim_map_hdf5(hdf5file)
   @test saved_map == g.map

   # particular individual
   for i=1:length(g.individual)
      individual = QMSimFiles.read_qmsim_individual_hdf5(g.map,hdf5file,i)
      @test individual == g.individual[i]
   end

   # all genotypes
   genotypes = QMSimFiles.read_qmsim_all_genotypes_hdf5(g.map,hdf5file)
   @test g == QMSimPopulationGenome(g.map,genotypes)

   # whole data
   @test g == read_qmsim_data_hdf5(hdf5file)

   # added individual
   add_qmsim_individual_hdf5(g.map, hdf5file, g.individual[1])
   added = QMSimFiles.read_qmsim_individual_hdf5(g.map,hdf5file,5)
   @test added == g.individual[1]

   # size of saved genotypes
   @test get_qmsim_genotype_size(hdf5file) == (5, g.map.totalSNP + g.map.totalQTL)

   # empty hdf5 only with map
   hdf5file2 = tempname()
   @info "temporary HDF5 file: $(hdf5file2)"
   create_qmsim_map_hdf5(g.map, hdf5file2)

   saved_map = read_qmsim_map_hdf5(hdf5file2)
   @test saved_map == g.map

   add_qmsim_individual_hdf5(g.map, hdf5file2, g.individual[1])
   added = QMSimFiles.read_qmsim_individual_hdf5(g.map,hdf5file2,1)
   @test added == g.individual[1]
end

@testset "export to blupf90" begin
   g = read_qmsim_data(snpmap1,qtlmap1,qtleff1,snpdata2,qtldata1)
   @test String( QMSimFiles.pack_markrs(g.map,g.individual[2]) .+ UInt8(48)) == "12000001010"
   @test String( QMSimFiles.pack_markrs(g.map,g.individual[3]) .+ UInt8(48)) == "22220001120"

   snpfile = tempname()
   @info "temporary text file: $(snpfile)"
   if isfile(snpfile); rm(snpfile,force=true); end
   export_qmsim_individual_blupf90(g.map, 2, g.individual[2], snpfile)
   export_qmsim_individual_blupf90(g.map, 3, g.individual[3], snpfile)

   @test countlines(snpfile)==2
   open(snpfile,"r") do io
      # individual 2
      line = readline(io)
      items = split(line)
      @test items[1]=="2" && items[2]=="12000001010"
      # individual 3
      line = readline(io)
      items = split(line)
      @test items[1]=="3" && items[2]=="22220001120"
   end
end


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")

function mreplace(x)
    return replace(x, missing=>NaN)
end

function calΨ!(
    ilev  :: Array{Float64, 1},
    lat   :: Array{Float64, 1},
    v     :: Array{Float64, 2},
    Ψ :: Array{Float64, 2},
)

    Ny, Nz = size(v)
    Nz_bnd = length(ilev)
    dp = ilev[1:end-1] - ilev[2:end]

    Ψ[:, Nz_bnd] .= 0.0
    for j=1:Ny
        for k=Nz_bnd-1:-1:1
            Ψ[j, k] = Ψ[j, k+1] + dp[k] * v[j, k]
        end

        Ψ[j, :] .*= 2π * Re * cos(lat[j] * π/180.0) / g
    end

end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin
        "--input-data-file-prefix"
            help = "Input data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--output-data-file-prefix"
            help = "Output data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--ilev-varname"
            help = "Variable name of pressure level. Assume in unit of hPa."
            arg_type = String
            default = "ilev"

        "--V-varname"
            help = "Variable name of meridional velocity"
            arg_type = String
            default = "V"

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

        "--overwrite"
            action = :store_true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global lon = replace(ds["xc"][:, 1], missing=>NaN)
    global lat = replace(ds["yc"][1, :], missing=>NaN)
end

Ψ   = nothing
ilev = nothing
for y=parsed["beg-year"]:parsed["end-year"], m=1:12

    i_filename = format("{:s}{:04d}-{:02d}.nc", parsed["input-data-file-prefix"], y, m)
    o_filename = format("{:s}{:04d}-{:02d}.nc", parsed["output-data-file-prefix"], y, m)

    if !isfile(o_filename) || parsed["overwrite"]

        Dataset(i_filename, "r") do ds

            v = convert(Array{Float64}, nomissing(ds[parsed["V-varname"]][:, :, 1]))

            if ilev == nothing
                global Ny, Nz = size(v)
                global ilev = ( ds[parsed["ilev-varname"]][:] |> mreplace ) * 100.0

                global Ψ = zeros(Float64, Ny, length(ilev))
            end

            calΨ!(ilev, lat, v, Ψ)

            Dataset(o_filename, "c") do ds
                defDim(ds, "Ny", Ny)
                defDim(ds, "ilev", length(ilev))
                defDim(ds, "time", Inf)
                for (varname, vardata, vardim, attrib) in [
                    ("psi", Ψ, ("Ny", "ilev"), Dict()),
                    ("ilev", ilev, ("ilev",), Dict()),
                ]

                    var = defVar(ds, varname, Float64, vardim)
                    var.attrib["_FillValue"] = 1e20

                    var = ds[varname]

                    for (k, v) in attrib
                        var.attrib[k] = v
                    end

                    rng = []
                    for i in 1:length(vardim)-1
                        push!(rng, Colon())
                    end
                    push!(rng, 1:size(vardata)[end])
                    var[rng...] = vardata

                end
            end
            println("Output: ", o_filename)
        end

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
text = "Hello world"
println(text)

name = "Roger"
println("Hello $name")

println("1 + 1 = $(1 + 1)")

😄= "hi"
println(😄)

ϵ = 2.2

println( π)

info = (name="Shooter", age="18", )

println(info.name)

max_stars_repo_name: liushooter/OneDayOneCommit

max_issues_repo_name: liushooter/OneDayOneCommit

max_forks_repo_name: liushooter/OneDayOneCommit


text:
using CITIfile: rdciti
using Test

s,f = rdciti("Wire_mesh_knit1_scan_th045deg_ph345deg.cit")
@testset "MA-style file" begin
    @test size(s) == (4,4,12)
    @test size(f) == (12,)
    @test abs(s[1,1,1]) ≈ 9.993692243725852E-01
    @test rad2deg(angle(s[1,1,1])) ≈  1.776698273015162E+02
    @test abs(s[1,1,end]) ≈ 9.869994561177106E-01
    @test rad2deg(angle(s[1,1,end])) ≈ 1.517921770116529E+02
    @test abs(s[1,3,3]) ≈ 3.261798828943794E-02
    @test rad2deg(angle(s[1,3,3])) ≈ 8.669814811044907E+01
    @test abs(s[3,1,3]) ≈ 3.263875888744300E-02
    @test rad2deg(angle(s[3,1,3])) ≈ 8.703675893037459E+01
    @test abs(s[end,end,end]) ≈ 8.073910702105189E-01
    @test rad2deg(angle(s[end,end,end])) ≈ 1.474292363421648E+02
end


text:
module Biogeochemistry

export generate_nutrients
export nutrients_init, default_nut_init
export nut_update!

using KernelAbstractions
using CUDA

using PlanktonIndividuals.Architectures: device, array_type, Architecture
using PlanktonIndividuals.Grids
using PlanktonIndividuals.Fields

include("Advection/Advection.jl")

using .Advection

include("nutrient_fields.jl")
include("nutrient_forcings.jl")
include("nutrient_update.jl")

end


text:
function ConsistentElDims(X::AbstractVector)
    elDim = length(X[1])
    all(x -> length(x) == elDim, X) ? elDim : throw("Inconsistent element lengths for given Vector.")
end
ConsistentElDims(X::AbstractVector{<:Number}) = 1
ConsistentElDims(M::AbstractMatrix{<:Number}) = size(M,2)
ConsistentElDims(T::Tuple) = ConsistentElDims(collect(T))

function HealthyCovariance(M::AbstractMatrix{<:Number})
    M = size(M,1) > size(M,2) ? Unwind(M) : M
    M = if isdiag(M)
        Diagonal(floatify(M))
    elseif sum(abs, M - Diagonal(floatify(M))) < 1e-20
        Diagonal(floatify(M))
    else
        floatify(M)
    end
    if !isposdef(M)
        println("Given Matrix not perfectly positive-definite. Using only upper half and symmetrizing.")
        M = Symmetric(M)
        !isposdef(M) && throw("Matrix still not positive-definite after symmetrization.")
        M = convert(Matrix, M)
    end
    return M
end
HealthyCovariance(D::Diagonal) = all(x->x>0, D.diag) ? D : throw("Given covariance Matrix has non-positive values on diagonal: $(D.diag)")
# Interpret vector as uncertainties, therefore square before converting to Matrix
HealthyCovariance(X::AbstractVector{<:Number}) = all(x->x>0, X) ? Diagonal(floatify(X).^2) : throw("Not all given uncertainties positive: $(X)")
HealthyCovariance(X::AbstractVector{<:AbstractVector{<:Number}}) = Unwind(X)

# Need to implement for each DataSet:   xdata, ydata, sigma, xsigma, ysigma, InvCov, Npoints, xdim, ydim, dims
#                                       WoundX (already generic), logdetInvCov (already generic), length (already generic)
# Need to implement for each DataModel: the above and: Data, model (Predictor), dmodel (dPredictor),
#                                       pdim, MLE, LogLikeMLE,
#                                       EmbeddingMap (already (kind of) generic), EmbeddingMatrix (already (kind of) generic)
#                                       Score (already (kind of) generic), FisherMetric (already (kind of) generic)

# Generic Methods for AbstractDataSets      -----       May be superceded by more specialized functions!
import Base.length
length(DS::AbstractDataSet) = Npoints(DS)
WoundX(DS::AbstractDataSet) = Windup(xdata(DS),xdim(DS))
# Becomes dangerous now that there is a distinction between yInvCov and xInvCov
# logdetInvCov(DS::AbstractDataSet) = logdet(InvCov(DS))
DataspaceDim(DS::AbstractDataSet) = Npoints(DS) * ydim(DS)
#sigma(DS::AbstractDataSet) = ysigma(DS)
@deprecate sigma(x) ysigma(x) true

xdist(DS::AbstractDataSet) = xDataDist(DS)
ydist(DS::AbstractDataSet) = yDataDist(DS)

dist(DS::AbstractDataSet) = GeneralProduct([xdist(DS), ydist(DS)])

Npoints(DS::AbstractDataSet) = Npoints(dims(DS))
xdim(DS::AbstractDataSet) = xdim(dims(DS))
ydim(DS::AbstractDataSet) = ydim(dims(DS))
Npoints(dims::Tuple{Int,Int,Int}) = dims[1]
xdim(dims::Tuple{Int,Int,Int}) = dims[2]
ydim(dims::Tuple{Int,Int,Int}) = dims[3]

_WoundX(DS::AbstractDataSet, WoundX::Nothing) = xdata(DS)
_WoundX(DS::AbstractDataSet, WoundX::AbstractVector) = WoundX

# Can eliminate specialized methods
function InformNames(DS::T, xnames::AbstractVector{String}, ynames::AbstractVector{String}) where T <: AbstractDataSet
    @assert length(xnames) == xdim(DS) && length(ynames) == ydim(DS)
    remake(DS; xnames=xnames, ynames=ynames)
end

# Generic Methods for AbstractDataModels      -----       May be superceded by more specialized functions!
pdim(DM::AbstractDataModel) = pdim(Data(DM), Predictor(DM))
MLE(DM::AbstractDataModel) = FindMLE(DM)
LogLikeMLE(DM::AbstractDataModel) = loglikelihood(DM, MLE(DM))
LogPrior(DM::AbstractDataModel) = x->0.0

# Generic passthrough of queries from AbstractDataModel to AbstractDataSet for following functions:
for F in [  :xdata, :ydata, :xsigma, :ysigma, :xInvCov, :yInvCov,
            :dims, :length, :Npoints, :xdim, :ydim, :DataspaceDim,
            :logdetInvCov, :WoundX, :xnames, :ynames, :xdist, :ydist, :dist]
    @eval $F(DM::AbstractDataModel) = $F(Data(DM))
end

# Generic Methods which are not simply passed through
pnames(DM::AbstractDataModel) = pnames(DM, Predictor(DM))
pnames(DM::AbstractDataModel, M::ModelMap) = pnames(M)
pnames(DM::AbstractDataModel, F::Function) = CreateSymbolNames(pdim(DM),"θ")

Domain(DM::AbstractDataModel) = Predictor(DM) isa ModelMap ? Domain(Predictor(DM)) : FullDomain(pdim(DM))

function AutoDiffDmodel(DS::AbstractDataSet, model::Function; custom::Bool=false, ADmode::Union{Symbol,Val}=Val(:ForwardDiff), Kwargs...)
    Grad, Jac = _GetGrad(ADmode; Kwargs...), _GetJac(ADmode; Kwargs...)
    GradPass, JacPass = _GetGradPass, _GetJacPass
    ## Allow for symbolic passthrough here
    Autodmodel(x::Number,θ::AbstractVector{<:Number}; kwargs...) = transpose(Grad(z->model(x,z; kwargs...),θ))
    Autodmodel(x::Number,θ::AbstractVector{<:Num}; kwargs...) = transpose(GradPass(z->model(x,z; kwargs...),θ))
    NAutodmodel(x::AbstractVector{<:Number},θ::AbstractVector{<:Number}; kwargs...) = transpose(Grad(z->model(x,z; kwargs...),θ))
    NAutodmodel(x::AbstractVector{<:Number},θ::AbstractVector{<:Num}; kwargs...) = transpose(GradPass(z->model(x,z; kwargs...),θ))
    AutodmodelN(x::Number,θ::AbstractVector{<:Number}; kwargs...) = Jac(p->model(x,p; kwargs...), θ)
    AutodmodelN(x::Number,θ::AbstractVector{<:Num}; kwargs...) = JacPass(p->model(x,p; kwargs...), θ)
    NAutodmodelN(x::AbstractVector{<:Number},θ::AbstractVector{<:Number}; kwargs...) = Jac(p->model(x,p; kwargs...),θ)
    NAutodmodelN(x::AbstractVector{<:Number},θ::AbstractVector{<:Num}; kwargs...) = JacPass(p->model(x,p; kwargs...),θ)
    # Getting extract_gradient! error from ForwardDiff when using gradient method with observables
    # CustomAutodmodel(x::Union{Number,AbstractVector{<:Number}},θ::AbstractVector{<:Number}) = transpose(Grad(p->model(x,p),θ))
    CustomAutodmodelN(x::Union{Number,AbstractVector{<:Number}},θ::AbstractVector{<:Number}; kwargs...) = Jac(p->model(x,p; kwargs...),θ)
    CustomAutodmodelN(x::Union{Number,AbstractVector{<:Number}},θ::AbstractVector{<:Num}; kwargs...) = JacPass(p->model(x,p; kwargs...),θ)
    if ydim(DS) == 1
        custom && return CustomAutodmodelN
        return xdim(DS) == 1 ? Autodmodel : NAutodmodel
    else
        custom && return CustomAutodmodelN
        return xdim(DS) == 1 ? AutodmodelN : NAutodmodelN
    end
end

"""
    DetermineDmodel(DS::AbstractDataSet, model::Function; ADmode::Union{Symbol,Val}=:ForwardDiff)::Function
Returns appropriate function which constitutes the automatic derivative of the `model(x,θ)` with respect to the parameters `θ` depending on the format of the x-values and y-values of the DataSet.
"""
function DetermineDmodel(DS::AbstractDataSet, model::Function; custom::Bool=false, ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...)
    # For the symbolically generated jacobians to work with MArrays, it requires ≥ v0.11.3 of SymbolicUtils.jl:  https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/286
    if ADmode === :Symbolic || ADmode isa Val{:Symbolic}
        Symbolic_dmodel = Optimize(DS, model; inplace=false)[2]
        !isnothing(Symbolic_dmodel) && return Symbolic_dmodel
        # Fall back to ForwarDiff if Symbolic differentiation did not work
        @info "Falling back to ForwardDiff for model jacobian."
        ADmode = Val(:ForwardDiff)
    end
    AutoDiffDmodel(DS, model; custom=custom, ADmode=ADmode, kwargs...)
end
function DetermineDmodel(DS::AbstractDataSet, M::ModelMap; custom::Bool=iscustom(M), kwargs...)
    ModelMap(DetermineDmodel(DS, M.Map; custom=custom, kwargs...), M)
end

MeasureAutoDiffPerformance(DM::AbstractDataModel; kwargs...) = MeasureAutoDiffPerformance(Data(DM), Predictor(DM), MLE(DM); kwargs...)
function MeasureAutoDiffPerformance(DS::AbstractDataSet, model::ModelOrFunction, mle::AbstractVector; kwargs...)
    modes = [:ForwardDiff, :Zygote, :ReverseDiff, :Symbolic, :FiniteDiff]
    perfs = Vector{Float64}(undef, length(modes))
    for (i,mode) in enumerate(modes)
        perfs[i] = try
            dmodel = DetermineDmodel(DS, model; ADmode=mode, kwargs...)
            @assert EmbeddingMatrix(DS, dmodel, mle) isa AbstractMatrix
            @belapsed EmbeddingMatrix($DS, $dmodel, $mle)
        catch;
            Inf
        end
    end
    M = sortslices([round.(perfs;sigdigits=4) modes]; dims=1)
    Res = [M[:,2] M[:,1]]

    println("$(Res)")
    println("$(Res[1,1]) performed best at $(Res[1,2])s. Second best was $(Res[2,1]) at $(Res[2,2])s.")
    Res
end

function CheckModelHealth(DS::AbstractDataSet, model::ModelOrFunction)
    P = ones(pdim(DS,model));   X = xdim(DS) < 2 ? xdata(DS)[1] : xdata(DS)[1:xdim(DS)]
    try  model(X,P)   catch Err
        throw("Got xdim=$(xdim(DS)) but model appears to not accept x-values of this size.")
    end
    out = model(X,P)
    !(size(out,1) == ydim(DS)) && println("Got ydim=$(ydim(DS)) but output of model does not have this size.")
    !(out isa SVector || out isa MVector) && (1 < ydim(DS) < 90) && @info "It may be beneficial for the overall performance to define the model function such that it outputs static vectors, i.e. SVectors."
    return
end

DataDist(Y::AbstractVector, Sig::AbstractVector, Dist=Normal) = product_distribution([Dist(Y[i],Sig[i]) for i in eachindex(Y)])
DataDist(Y::AbstractVector, Sig::AbstractMatrix, Dist=MvNormal) = Dist(Y, HealthyCovariance(Sig))
yDataDist(DS::AbstractDataSet) = DataDist(ydata(DS), ysigma(DS))
xDataDist(DS::AbstractDataSet) = xsigma(DS) == zeros(Npoints(DS)*xdim(DS)) ? InformationGeometry.Dirac(xdata(DS)) : DataDist(xdata(DS), HealthyCovariance(xsigma(DS)))
yDataDist(DM::AbstractDataModel) = yDataDist(Data(DM))
xDataDist(DM::AbstractDataModel) = xDataDist(Data(DM))

"""
    pdim(DS::AbstractDataSet, model::ModelOrFunction) -> Int
Infers the (minimal) number of components that the given function `F` accepts as input by successively testing it on vectors of increasing length.
"""
pdim(DS::AbstractDataSet, model::ModelOrFunction) = xdim(DS) < 2 ? GetArgLength(p->model(xdata(DS)[1],p)) : GetArgLength(p->model(xdata(DS)[1:xdim(DS)],p))

# DataFrame(DM::DataModel) = DataFrame(Data(DM))
DataFrame(DS::AbstractDataSet; kwargs...) = SaveDataSet(DS; kwargs...)

import Base.join
function join(DS1::T, DS2::T) where T <: AbstractDataSet
    !(xdim(DS1) == xdim(DS2) && ydim(DS1) == ydim(DS2)) && throw("DataSets incompatible.")
    NewΣ = if typeof(ysigma(DS1)) <: AbstractVector && typeof(ysigma(DS2)) <: AbstractVector
        vcat(ysigma(DS1), ysigma(DS2))
    else
        BlockMatrix(ysigma(DS1), ysigma(DS2))
    end
    DataSet(vcat(xdata(DS1), xdata(DS2)), vcat(ydata(DS1), ydata(DS2)), NewΣ, (Npoints(DS1)+Npoints(DS2), xdim(DS1), ydim(DS1)))
end
join(DM1::AbstractDataModel, DM2::AbstractDataModel) = DataModel(join(Data(DM1),Data(DM2)), Predictor(DM1), dPredictor(DM1))
join(DS1::T, DS2::T, args...) where T <: Union{AbstractDataSet,AbstractDataModel} = join(join(DS1,DS2), args...)
join(DSVec::AbstractVector{T}) where T <: Union{AbstractDataSet,AbstractDataModel} = join(DSVec...)

SortDataSet(DS::AbstractDataSet) = DS |> DataFrame |> sort |> DataSet
SortDataModel(DM::AbstractDataModel) = DataModel(SortDataSet(Data(DM)), Predictor(DM), dPredictor(DM), MLE(DM))
function SubDataSet(DS::AbstractDataSet, range::Union{AbstractVector{<:Int},BoolVector})
    @assert DS isa DataSet || xdist(DS) isa InformationGeometry.Dirac
    Npoints(DS) < length(range) && throw("Length of given range unsuitable for DataSet.")
    X = WoundX(DS)[range] |> Unwind
    Y = Windup(ydata(DS),ydim(DS))[range] |> Unwind
    Σ = ysigma(DS)
    if typeof(Σ) <: AbstractVector
        Σ = Windup(Σ,ydim(DS))[range] |> Unwind
    elseif ydim(DS) == 1
        Σ = Σ[range,range]
    else
        throw("Under construction.")
    end
    InformNames(DataSet(X,Y,Σ,(Int(length(X)/xdim(DS)),xdim(DS),ydim(DS))), xnames(DS), ynames(DS))
end
SubDataModel(DM::AbstractDataModel, range::Union{AbstractVector{<:Int},BoolVector}) = DataModel(SubDataSet(Data(DM),range), Predictor(DM), dPredictor(DM), MLE(DM))

import Base: getindex, lastindex
getindex(DS::AbstractDataSet, x) = SubDataSet(DS, x)
lastindex(DS::AbstractDataSet) = Npoints(DS)

Sparsify(DS::AbstractDataSet) = SubDataSet(DS, rand(Bool,Npoints(DS)))
Sparsify(DM::AbstractDataModel) = SubDataSet(DS, rand(Bool,Npoints(DS)))


text:
mutable struct DuplexInterval{T}
   first::Interval{T}  # fwd helical region
   last::Interval{T}   # rev helical region
   duplex::AbstractDuplex   # rna duplex formed
end

distance( int::DuplexInterval ) = (int.last.first - int.first.last) - 1

const ZERO_DUPLEX_INTERVAL = DuplexInterval(Interval("", 0, 0), Interval("", 0, 0), RNADuplex())

# a DuplexCollection contains DuplexIntervals indexed by first interval
mutable struct DuplexCollection{T}
   names::Dict{String,Dict{UInt64,Vector{DuplexInterval{T}}}}
   binsize::Int64
   length::Int64

   DuplexCollection{T}() where T               = new(Dict{String,Dict{UInt64,Vector{DuplexInterval{T}}}}(), 50, 0 )
   DuplexCollection{T}( binsize::Int ) where T = new(Dict{String,Dict{UInt64,Vector{DuplexInterval{T}}}}(), binsize, 0)
end

function Base.minimum{T}( col::DuplexCollection{T}; default::Float64=0.0 )
   minval = default
   minentry = Nullable{DuplexInterval{T}}()
   for n in keys(col.names)
      for b in keys(col.names[n])
         idx    = indmin(map(x->energy(x.duplex), col.names[n][b]))
         idxdg  = energy(col.names[n][b][idx].duplex)
         if idxdg < minval
            minval   = idxdg
            minentry = Nullable(col.names[n][b][idx])
         end
      end
   end
   minentry
end

function Base.collect{T}( col::DuplexCollection{T}; minlength::Int=1, minenergy::Float64=0.0, range::UnitRange=1:typemax(Int) )
   res = Vector{DuplexInterval{T}}()
   for n in keys(col.names)
      for b in sort(collect(keys(col.names[n])))
         for d in col.names[n][b]
            if npairs( d.duplex ) >= minlength &&
               energy( d.duplex ) <= minenergy &&
               d.last.first - d.first.last in range
                push!( res, d )
            end
         end
      end
   end
   res
end

Base.length{T}( col::DuplexCollection{T}; keep::Bool=false ) = col.length > 0 && keep ? col.length : calculate_length!( col )

function calculate_length!{T}( col::DuplexCollection{T} )
   for n in keys(col.names)
      for b in keys(col.names[n])
         col.length += length(col.names[n][b])
      end
   end
   col.length
end

# hash by identity for speed
Base.hash(x::UInt64) = x

@inline bin{I <: Integer}( val::I, size::Int ) = (v = div(val, size)*size; return UInt64(v+1))

function isoverlapping{S,T}(x::Interval{S}, y::Interval{T}, identity::Float64)
   if isoverlapping( x, y )
      overlap  = min( x.last, y.last ) - max( x.first, y.first ) + 1
      smallest = min( x.last - x.first, y.last - y.first ) + 1
      if overlap / smallest >= identity
         return true
      end
   end
   return false
end

@inline isoverlapping(x::DuplexInterval, y::DuplexInterval) = isoverlapping(x.first, y.first) &&
                                                              isoverlapping(x.last, y.last) ? true : false
@inline isoverlapping(x::DuplexInterval, y::DuplexInterval, ident::Float64) = isoverlapping(x.first, y.first, ident) &&
                                                                              isoverlapping(x.last, y.last, ident) ? true : false

@inline precedes(x::DuplexInterval, y::DuplexInterval) = GenomicFeatures.precedes(x.first, y.first)
@inline Base.isless(x::DuplexInterval, y::DuplexInterval) = (x.first.first == y.first.first && x.first.last == y.first.last) ?
                                                             GenomicFeatures.isless(x.last, y.last) : GenomicFeatures.isless(x.first, y.first)

Base.:<(x::DuplexInterval, y::DuplexInterval) = isless( x, y )

# return boolean of whether the duplex collection is
# properly ordered or not.
function isordered{T}( col::DuplexCollection{T} )
   for n in keys(col.names)
      const dict = col.names[n]
      prev = ZERO_DUPLEX_INTERVAL
      for b in sort(collect(keys(dict)))
         const vect = dict[b]
         for i in 1:length(vect)
            if !isless( prev, vect[i] )
               return false
            end
            prev = vect[i]
         end
      end
   end
   return true
end

function Base.delete!{T}(col::DuplexCollection{T}, int::DuplexInterval{T})
   const key = bin( int.first.first, col.binsize )
   if haskey( col.names, int.first.seqname ) && length(col.names[int.first.seqname][key]) > 0
      const vect = col.names[int.first.seqname][key]
      delete!( vect, int )
      if length(vect) == 0
         delete!( col.names[int.first.seqname], key )
      end
   end
end

function Base.delete!{T}( vect::Vector{DuplexInterval{T}}, int::DuplexInterval{T} )
   j = searchsortedfirst( vect, int )
   while j <= length(vect)
      if vect[j] == int
         deleteat!( vect, j )
      else
         j += 1
      end
   end
end

@inline function Base.push!{T}(col::DuplexCollection{T}, int::DuplexInterval{T})
   if !haskey(col.names, int.first.seqname)
      col.names[int.first.seqname] = Dict{UInt64,Vector{DuplexInterval{T}}}()
   end
   key = bin( int.first.first, col.binsize )
   # check bins before and after for better matches
   if !findbetter(col.names[int.first.seqname], key-0x01, int, true) &&
      !findbetter(col.names[int.first.seqname], key+0x01, int, false)
        pushinterval!(col.names[int.first.seqname], key, int)
   end
end

@inline function findbetter{T, K}(dict::Dict{K,Vector{DuplexInterval{T}}}, key::K, int::DuplexInterval{T}, left::Bool=true, rate::Float64=0.99)
   !haskey( dict, key ) && (return false)
   const vect = dict[key]
   i = left ? searchsortedfirst( vect, int, lt=precedes ) : 1
   while i <= length(vect)
      if left && precedes( int, vect[j] )
         return false
      elseif isoverlapping( vect[i], int, rate )
         if energy(vect[i].duplex) < energy(int.duplex)
            return true
         else
            deleteat!( vect, i )
            continue
         end
      end
      i += 1
   end
   if length(vect) == 0
      delete!( dict, key )
   end
   return false
end

function pushinterval!{T}(vect::Vector{DuplexInterval{T}}, int::DuplexInterval{T}; rate::Float64=0.99)
   j = searchsortedfirst( vect, int, lt=precedes )
   while j <= length(vect)
      if isoverlapping( int, vect[j], rate )
         if energy(int.duplex) < energy(vect[j].duplex)
            splice!( vect, j )
            j -= 1
         else
            return
         end
      end
      j += 1
   end
   insert!( vect, searchsortedfirst( vect, int ), int )
end

@inline function pushinterval!{T, K}(dict::Dict{K,Vector{DuplexInterval{T}}}, key::K, int::DuplexInterval{T}; rate::Float64=0.99)
   if !haskey( dict, key )
      dict[key] = Vector{DuplexInterval{T}}()
   end
   const vect = dict[key]
   pushinterval!( vect, int, rate=rate )
end

join_duplex!( left::RNADuplex, right::RNADuplex, npairs, npairs_first, npairs_last ) = push!( left, path(right)[npairs+1:end] )

# This function 'stitches' two overlapping and compatible
# duplexes together into one.
@inline function _stitch{T}( a::DuplexInterval{T}, b::DuplexInterval{T}, max_bulge::Int, max_mis::Int )
   ret = Nullable{DuplexInterval{T}}()
   if isoverlapping( a, b )
      npairs_first = a.first.last - b.first.first + 1
      npairs_last  = b.last.last  - a.last.first  + 1
      afirst,alast = strings(path(a.duplex))
      bfirst,blast = strings(path(b.duplex))
      npairs = index_npairs(path(b.duplex), npairs_first, npairs_last)
      if abs(npairs_first - npairs_last) <= max_bulge &&
         length(bfirst) - npairs_first > 0 &&
         length(blast)  - npairs_last  > 0 &&
         npairs_first <= length(afirst) && npairs_first <= length(bfirst) &&
         npairs_last  <= length(alast)  && npairs_last  <= length(blast) &&
         afirst[end-npairs_first+1:end] == bfirst[1:npairs_first] &&
         alast[end-npairs_last+1:end]   == blast[1:npairs_last]

         #spliced = deepcopy(a)
         spliced_first = Interval{T}( a.first.seqname, a.first.first, b.first.last, a.first.strand, a.first.metadata )
         spliced_last  = Interval{T}( a.last.seqname,  b.last.first,  a.last.last,  a.last.strand,  a.last.metadata )
         spliced = DuplexInterval{T}( spliced_first, spliced_last, deepcopy(a.duplex) )
         join_duplex!( spliced.duplex, b.duplex, npairs, npairs_first, npairs_last )
         #spliced.first.last = b.first.last
         #spliced.last.first = b.last.first
         aspliced, bspliced = strings(spliced.duplex)
              #println("$a\n$b\n")
              #println("$afirst : $alast")
              #println("$bfirst : $blast")
              #println("$npairs_first : $npairs_last : $npairs")
              #println("$(b.duplex.path[npairs+1:end])")
              #println("RESULT: $spliced")

         if spliced.first.last - spliced.first.first + 1 == length(aspliced) &&
            spliced.last.last  - spliced.last.first  + 1 == length(bspliced) &&
            energy(spliced.duplex) < energy(a.duplex) &&
            energy(spliced.duplex) < energy(b.duplex) &&
            nbulges(spliced.duplex) <= max_bulge &&
            nmismatches(spliced.duplex) <= max_mis
             ret = Nullable(spliced)
         end
      end
   end
   ret
end

stitch{T}( a::DuplexInterval{T}, b::DuplexInterval{T}, max_bulge::Int, max_mis::Int ) = a.first < b.first ?
                                                                                _stitch( a, b, max_bulge, max_mis ) :
                                                                                _stitch( b, a, max_bulge, max_mis )
#perform stitching on duplex collection
function stitch{T}(col::DuplexCollection{T})
   colret = deepcopy(col)
   for n in keys(col.names)
      res  = Vector{DuplexInterval{T}}()
      bins = sort(collect(keys(col.names[n])))
      for i in 1:length(bins)
         if i < length(bins) && bins[i+1] == bins[i]+1
            rec_stitch!( res, [ col.names[n][bins[i]]; col.names[n][bins[i+1]] ] )
         else
            rec_stitch!( res, col.names[n][bins[i]] )
         end
      end
      for r in res
         push!( colret, r )
      end
   end
   colret
end

function stitch!{T}( res::Vector{DuplexInterval{T}}, vect::Vector{DuplexInterval{T}} )
   for i in 1:length(vect)
      for j in i:searchsortedlast( vect, vect[i], lt=precedes )
         j == i && continue
         sval = stitch( vect[i], vect[j], 3, 3 )
         if !isnull(sval)
            pushinterval!( res, sval.value )
         end
      end
   end
end

function rec_stitch!{T}( res::Vector{DuplexInterval{T}}, vect::Vector{DuplexInterval{T}} )
   newvec = Vector{DuplexInterval{T}}()
   stitch!( newvec, vect )
   if length(newvec) > 0
      rec_stitch!( res, newvec )
      for r in newvec
         pushinterval!( res, r )
      end
   end
   res
end


text:
using Graphics: @mustimplement

export SurveillanceUnit, getSurveillanceUnits, getSurveillanceUnit, enableACPower, disableACPower, getTemperatures, getACStatus, resetDAQ

abstract type SurveillanceUnit <: Device end

include("DummySurveillanceUnit.jl")
include("ArduinoSurveillanceUnit.jl")
include("ArduinoSurveillanceUnitExternalTemp.jl")
include("ArduinoSurveillanceUnitInternalTemp.jl")
#include("MPSSurveillanceUnit.jl") TODO

Base.close(su::SurveillanceUnit) = nothing

@mustimplement getTemperatures(su::SurveillanceUnit)
@mustimplement getACStatus(su::SurveillanceUnit)
@mustimplement enableACPower(su::SurveillanceUnit)
@mustimplement disableACPower(su::SurveillanceUnit)
@mustimplement resetDAQ(su::SurveillanceUnit)
@mustimplement hasResetDAQ(su::SurveillanceUnit) # TODO is this has as in "was reset successfull" or if it has the ability to do so?

getSurveillanceUnits(scanner::MPIScanner) = getDevices(scanner, SurveillanceUnit)
getSurveillanceUnit(scanner::MPIScanner) = getDevice(scanner, SurveillanceUnit)


text:
function _convert_to_svector(vector::AbstractVector{T}) where {T <: Real}

    n_elements = length(vector)
    type_elements = eltype(vector)

    return SVector{n_elements, type_elements}(vector)
end

function _convert_to_svector(vectors)

    vector_first = vectors[1]

    n_elements = length(vector_first)
    type_elements = eltype(vector_first)

    vectors_static = []

    for vector in vectors
        push!(vectors_static, SVector{n_elements, type_elements}(vector))
    end

    return vectors_static
end


text:
#reload("../../GaussianProcess/src/GaussianProcess.jl")
module t
using PyPlot
import GaussianProcess
#rf = SKLearnRF.RandomForest(100)
# fitc = GPyFITC.FITCModel(50)

k = GaussianProcess.MaternIso32(0, 0)
lik = GaussianProcess.GaussianLikelihood(0)
fullgp = GaussianProcess.FullModel(k, lik)
c = GaussianProcess.OneD(13)

X = linspace(-10, 10, 7)''
X = vcat(X, linspace(0, 3, 10)'')
y = exp(-(X[:].^2) / 5)  - exp(-((X[:] - 4).^2) / 3)
Xp = linspace(-10, 10, 500)''

close(1);figure(1, figsize=(6, 9.6));clf();
function plot_it(X, y, Xp, mu, sigma2)
	plot(X, y, "kx")
	plot(Xp, yp)
	fill_between(Xp[:], (yp + 2sqrt(sigma2))[:], (yp - 2sqrt(sigma2))[:], alpha=0.5)
	locs, labels = yticks()
	yticks([locs[1]*1.5, 0, locs[end]*1.5])
	xticks([])
	ylim([-2.5, 2.5])
	xlim([-10, 10])
end

inducing = 50
subplot(4, 1, 1)
title("Full GP")
GaussianProcess.observe!(fullgp, X, y)
GaussianProcess.infer!(GaussianProcess.MLE2([-10, -10, -10],
			[10, 10, 10]), fullgp)
yp, sigma2 = GaussianProcess.predict(fullgp, Xp)
plot_it(X, y, Xp, yp, sigma2)

subplot(4, 1, 2)
title("Nyström ($inducing basis functions)")
nystrom = GaussianProcess.NystromModel(k, lik, [10.], [-10.], neig_samples=inducing)
GaussianProcess.observe!(nystrom, X, y)
GaussianProcess.infer!(GaussianProcess.MLE2([-10, -10, -10],
			[10, 10, 10]), nystrom)
yp, nsigma2 = GaussianProcess.predict(nystrom, Xp)
plot_it(X, y, Xp, yp, nsigma2)

subplot(4, 1, 3)
title("Laplacian ($inducing basis functions)")
laplacian = GaussianProcess.LaplacianModel(k, lik, c, inducing)
GaussianProcess.observe!(laplacian, X, y)
GaussianProcess.infer!(GaussianProcess.MLE2([-10, -10, -10],
			[10, 10, 10]), laplacian)
yp, lsigma2 = GaussianProcess.predict(laplacian, Xp)
plot_it(X, y, Xp, yp, lsigma2)

subplot(4, 1, 4)
u = linspace(-10, 10, inducing)''
title("FITC ($inducing inducing points)")
fitc = GaussianProcess.FITCModel(k, lik, u)
GaussianProcess.observe!(fitc, X, y)
GaussianProcess.infer!(GaussianProcess.MLE2(
			[-10, -10, -10, -10*ones(size(u, 1))],
			[10, 10, 10, 10*ones(size(u, 1))]), fitc)
yp, sigma2 = GaussianProcess.predict(fitc, Xp)
plot_it(X, y, Xp, yp, sigma2)

tight_layout()

end


text:
using Test

import ShuffleProofs: prove, verify, Simulator, gen_shuffle, Verifier, PoSChallenge, Shuffle, shuffle, VShuffleProof, PoSProof, ProtocolSpec, gen_roprg

import CryptoGroups: ElGamal, PGroup, Enc, Dec

SPEC = "$(@__DIR__)/../ref/demo/protInfo.xml"
verifier = ProtocolSpec(SPEC)

(; g) = verifier
sk = 123
pk = g^sk

enc = Enc(pk, g)

𝐦 = [g^4, g^2, g^3]
𝐞 = enc(𝐦, [2, 3, 4])

N = length(𝐞)

𝐡 = [g^i for i in 2:N+1]

𝐫′ = [4, 2, 10]

proposition, secret = gen_shuffle(enc, 𝐞, 𝐫′) # In practice total of random factors can't match as it reveals

@test verify(proposition, secret)

simulator = prove(proposition, secret, verifier)
(; proof) = simulator

@test verify(simulator)

### Testing proof translation and verification with Verificatum notation written verifier
vproof = VShuffleProof(proof)
@test verify(proposition, vproof, verifier)

### Higher order API

roprg = gen_roprg(UInt8[2])

simulator2 = shuffle(𝐞, g, pk, verifier; roprg)
@test verify(simulator2)