text:
using Test
using OrdinaryDiffEq
using LinearAlgebra, LinearSolve

import OrdinaryDiffEq.dolinsolve

n = 8
dt = 1/16
u0 = ones(n)
tspan = (0.0,1.0)

M1 = 2ones(n) |> Diagonal #|> Array
M2 = 2ones(n) |> Diagonal #|> Array

f1 = M1 |> DiffEqArrayOperator
f2 = M2 |> DiffEqArrayOperator
prob = SplitODEProblem(f1,f2,u0,tspan)

for algname in (
                :SBDF2,
                :SBDF3,
                :KenCarp47,
               )
    @testset "$algname" begin

    alg0 = @eval $algname()
    alg1 = @eval $algname(linsolve=GenericFactorization())

    kwargs = (dt=dt,)

    # expected error message
    msg = "Split ODE problem do not work with factorization linear solvers. Bug detailed in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Defaulting to linsolve=KrylovJL()"
    @test_logs (:warn, msg) solve(prob, alg0; kwargs...)
    @test DiffEqBase.__solve(prob, alg0; kwargs...).retcode == :Success
    @test_broken DiffEqBase.__solve(prob, alg1; kwargs...).retcode == :Success

    end
end

#####
# deep dive
#####

alg0 = KenCarp47()                                # passing case
alg1 = KenCarp47(linsolve=GenericFactorization()) # failing case

## objects
ig0 = SciMLBase.init(prob, alg0; dt=dt)
ig1 = SciMLBase.init(prob, alg1; dt=dt)

nl0 = ig0.cache.nlsolver
nl1 = ig1.cache.nlsolver

lc0 = nl0.cache.linsolve
lc1 = nl1.cache.linsolve

W0 = lc0.A
W1 = lc1.A

# perform first step
OrdinaryDiffEq.loopheader!(ig0)
OrdinaryDiffEq.loopheader!(ig1)

OrdinaryDiffEq.perform_step!(ig0,ig0.cache)
OrdinaryDiffEq.perform_step!(ig1,ig1.cache)

@test_broken OrdinaryDiffEq.nlsolvefail(nl0) # fails
@test        OrdinaryDiffEq.nlsolvefail(nl1) # passes

# check operators
@test        W0._concrete_form == W1._concrete_form
@test_broken W0._func_cache    == W1._func_cache

# check operator application
b = ones(n)
@test W0 * b == W1 * b
@test mul!(rand(n), W0, b) == mul!(rand(n), W1, b)
@test W0 \ b == W1 \ b

# check linear solve
lc0.b .= 1.0
lc1.b .= 1.0

solve(lc0)
solve(lc1)

@test_broken lc0.u == lc1.u

# solve contried problem using OrdinaryDiffEq machinery
linres0 = dolinsolve(ig0, lc0; A = W0, b = b, linu = ones(n), reltol = 1e-8)
linres1 = dolinsolve(ig1, lc1; A = W1, b = b, linu = ones(n), reltol = 1e-8)

@test_broken linres0 == linres1

###
# custom linsolve function
###

function linsolve(A,b,u,p,newA,Pl,Pr,solverdata;kwargs...)
    # avoid preconditioner monkeybusiness
    prob = LinearProblem(A,b;u0=u)
    solve(prob, nothing)
    return u
end

alg = KenCarp47(linsolve=LinearSolveFunction(linsolve))

@test solve(prob, alg).retcode == :Success

nothing


text:
@testset "RandomSelector" begin
    st = Date(2019, 1, 1)
    ed = Date(2019, 2, 1)
    date_range = st:Day(1):ed

    @testset "construction" begin
        selector = RandomSelector(2, 1//3)
        @test selector.seed == 2
        @test selector.holdout_fraction == 1//3
        @test selector.block_size == Day(1)
        @test selector.offset == Day(0)

        selector = RandomSelector(2, 1//3, Day(15))
        @test selector.seed == 2
        @test selector.holdout_fraction == 1//3
        @test selector.block_size == Day(15)
        @test selector.offset == Day(0)

        selector = RandomSelector(2, 1//3, Day(15), Day(5))
        @test selector.seed == 2
        @test selector.holdout_fraction == 1//3
        @test selector.block_size == Day(15)
        @test selector.offset == Day(5)

        @test (
            RandomSelector(99, 1//2, Day(1), Day(0)) ==
            RandomSelector(99, 1//2, Day(1)) ==
            RandomSelector(99, 1//2) ==
            RandomSelector(99)
        )
    end

    @testset "basic" begin
        selector = RandomSelector(42)
        @test selector isa DateSelector

        result = partition(date_range, selector)
        @test result isa NamedTuple{(:validation, :holdout)}
        @test sort(union(result...)) == date_range
        @test isempty(intersect(result...))

        # Running partition twice should return the same result
        result2 = partition(date_range, selector)
        @test result.validation == result2.validation
        @test result.holdout == result2.holdout

        @testset "holdout fraction" begin
            # Setting holdout_fraction 1 all days leaves the validation set empty
            validation, holdout = partition(date_range, RandomSelector(42, 1))
            @test isempty(validation)
            @test collect(holdout) == date_range

            # Setting holdout_fraction 0 all days leaves the holdout set empty
            validation, holdout = partition(date_range, RandomSelector(42, 0))
            @test isempty(holdout)
            @test collect(validation) == date_range
        end

        @testset "block size" begin
            selector = RandomSelector(42, 1//2, Day(3))
            result = partition(date_range, selector)
            @test sort(union(result...)) == date_range
            @test isempty(intersect(result...))

            for subset in result
                # at vary least the first 3 items must be from same block
                @test subset[2] == subset[1] + Day(1)
                @test subset[3] == subset[1] + Day(2)

                # no gaps of 2 are possible
                @test !any(isequal(2), diff(subset))
            end
        end
    end

    @testset "Right holdout fraction RandomSelector($seed, $holdout_fraction, $block_size)" for
        seed in (1, 42),
        holdout_fraction in (1//2, 0.1, 0.9, 0.7),
        block_size in (Day(1), Day(2), Week(1), Week(2))

        selector = RandomSelector(seed, holdout_fraction, block_size)
        range = Date(2000, 1, 1):Day(1):Date(2010, 1, 31)  # 10 year total range

        r = partition(range, selector)
        @test length(r.holdout)/length(range) ≈ holdout_fraction atol=0.05
    end

    @testset "Set Seed (protects against julia RNG changing)" begin
        r1 = partition(date_range, RandomSelector(42, 1//5))
        r2 = partition(date_range, RandomSelector(42, 1//5))
        r3 = partition(date_range, RandomSelector(42, 1//5))

        # cannot directly equate NamedTuples of iterators
        @test collect(r1.validation) == collect(r2.validation) == collect(r3.validation)
        @test collect(r1.holdout) == collect(r2.holdout) == collect(r3.holdout)
        @test isequal(r1.holdout, [Date(2019, 1, 1), Date(2019, 1, 5), Date(2019, 1, 30)])
    end

    @testset "Different date inputs" begin
        exp = partition(date_range, RandomSelector(42))

        @testset "$(typeof(d))" for d in (
            # Interval
            st..ed,
            # AnchoredInterval should include both start and dates
            AnchoredInterval{Day(31), Date, Closed, Closed}(st)
        )
            result = partition(d, RandomSelector(42))
            @test result.validation == exp.validation
            @test result.holdout == exp.holdout
        end
    end

    @testset "Invarient to date_range RandomSelector($seed, $holdout_fraction, $block_size, $offset))" for
        seed in (1, 42),
        holdout_fraction in (1//2, 1/20),
        block_size in (Day(1), Day(2), Week(1), Week(2)),
        offset in Day.(-3:3)

        selector = RandomSelector(seed, holdout_fraction, block_size, offset)  #only varient to the selector
        initial_range = Date(2019, 1, 1):Day(1):Date(2020, 1, 31)  # 1 year total range

        initial_sets = map(Set, partition(initial_range, selector))

        later_range = initial_range + Day(20)  # 20 days later.
        later_sets = map(Set, partition(later_range, selector))

        @test initial_sets.holdout ∩ later_sets.validation == Set()
        @test initial_sets.validation ∩ later_sets.holdout == Set()
    end
end

ext: jl


text:
@testset "examples" begin
    @test include("../examples/test_function_optimisation.jl")
    @test include("../examples/test_function_optimisation_noise.jl")
    @test include("../examples/test_function_optimisation_categorical.jl")
    @test include("../examples/test_function_optimisation_discontinuous.jl"); sleep(3)
end


text:
using LaserTypes
using LinearAlgebra
using Test

ω = 0.057
T₀ = 2π/ω
c = 137.036
λ = c*T₀
w₀ = 75 * λ
a₀ = 2.
s = setup_laser(LaguerreGaussLaser, :atomic, m = 1, p = 1, λ = λ, a₀ = a₀, w₀ = w₀,
    profile=ConstantProfile())

fieldE(x,y,z,t) = E([x,y,z],t,s)
fieldB(x,y,z,t) = B([x,y,z],t,s)

wenergy(x,y,t) = fieldE(x,y,0,t)⋅fieldE(x,y,0,t) + c*c * fieldB(x,y,0,t)⋅fieldB(x,y,0,t)
wenergy₀(x,y) = wenergy(x,y,0)

domainXY = -4*w₀:λ:4*w₀
wXY = map(r -> wenergy₀(r[1],r[2]), Iterators.product(domainXY,domainXY))
wMax = maximum(wXY)
@test wMax < 1e6


text:
function solve()
    n = parse(Int, readline())
    as = [parse(Int, x) for x in split(readline())]

    bs = zeros(Int, n)
    bs[1] = as[1]
    for i in 2:n
        bs[i] = as[i] + bs[i-1]
    end
    ret = (n-1) * sum(x->x^2, as)
    for i in 2:n
        ret -= 2*as[i]*bs[i-1]
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


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

mutable struct IoK8sApiStorageV1VolumeAttachmentSpec <: SwaggerModel
    attacher::Any # spec type: Union{ Nothing, String } # spec name: attacher
    nodeName::Any # spec type: Union{ Nothing, String } # spec name: nodeName
    source::Any # spec type: Union{ Nothing, IoK8sApiStorageV1VolumeAttachmentSource } # spec name: source

    function IoK8sApiStorageV1VolumeAttachmentSpec(;attacher=nothing, nodeName=nothing, source=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1VolumeAttachmentSpec, Symbol("attacher"), attacher)
        setfield!(o, Symbol("attacher"), attacher)
        validate_property(IoK8sApiStorageV1VolumeAttachmentSpec, Symbol("nodeName"), nodeName)
        setfield!(o, Symbol("nodeName"), nodeName)
        validate_property(IoK8sApiStorageV1VolumeAttachmentSpec, Symbol("source"), source)
        setfield!(o, Symbol("source"), source)
        o
    end
end # type IoK8sApiStorageV1VolumeAttachmentSpec

const _property_map_IoK8sApiStorageV1VolumeAttachmentSpec = Dict{Symbol,Symbol}(Symbol("attacher")=>Symbol("attacher"), Symbol("nodeName")=>Symbol("nodeName"), Symbol("source")=>Symbol("source"))
const _property_types_IoK8sApiStorageV1VolumeAttachmentSpec = Dict{Symbol,String}(Symbol("attacher")=>"String", Symbol("nodeName")=>"String", Symbol("source")=>"IoK8sApiStorageV1VolumeAttachmentSource")
Base.propertynames(::Type{ IoK8sApiStorageV1VolumeAttachmentSpec }) = collect(keys(_property_map_IoK8sApiStorageV1VolumeAttachmentSpec))
Swagger.property_type(::Type{ IoK8sApiStorageV1VolumeAttachmentSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1VolumeAttachmentSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1VolumeAttachmentSpec }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1VolumeAttachmentSpec[property_name]

function check_required(o::IoK8sApiStorageV1VolumeAttachmentSpec)
    (getproperty(o, Symbol("attacher")) === nothing) && (return false)
    (getproperty(o, Symbol("nodeName")) === nothing) && (return false)
    (getproperty(o, Symbol("source")) === nothing) && (return false)
    true
end

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


text:
# take a given molecular dynamics trial and save every tenth PE, KE, T, virial, and
# speed distribution to a file whose name is based on ρ and Tdesired.
function extract(folder::String, r::MolecularDynamicsTrial, N::Integer)
    rho = r.ρ
    Td  = r.Td

    speeds = √sum(r.v[:,:,N:N:end] .* r.v[:,:,N:N:end], 1)
    speeds = reshape(speeds, size(speeds)[2], size(speeds)[3])
    writedlm("$folder/speeds-rho-$rho-temp-$Td.dat", speeds)

    temps = r.T[N:N:end]
    writedlm("$folder/temps-rho-$rho-temp-$Td.dat", temps)

    ke = r.ket[N:N:end]
    writedlm("$folder/ke-rho-$rho-temp-$Td.dat", ke)

    pe = r.pet[N:N:end]
    writedlm("$folder/pe-rho-$rho-temp-$Td.dat", pe)

    vir = r.vir[N:N:end]
    writedlm("$folder/vir-rho-$rho-temp-$Td.dat", vir)
end

# save to current folder
extract(r::MolecularDynamicsTrial, N::Integer) = extract("", r, N)


text:
module FastCap

# ghostcript
const ps2pdf = joinpath(dirname(dirname(readstring(`where gswin64c`))),"lib","ps2pdf.bat")
# vbscript for FastFieldSolvers programs on Windows
const fastcap2_vbs  = joinpath(dirname(dirname(@__FILE__)),"deps","fastcap2.vbs")
const fastercap_vbs = joinpath(dirname(dirname(@__FILE__)),"deps","fastercap.vbs")
# Line endings
const NL = is_windows()?"\r\n":"\n"

export Cube, Sphere,
       Model, add_rectangle!, add_polygon!, add_physicalsurf!,
       extrude!,
       Conductor, Dielectric, ConductorGroup, DielectricGroup, move,
       fastcap2, fastercap, fastcapwr, linedrawing

include("shapes.jl")
include("geometries.jl")
include("commands.jl")
include("utilities.jl")
include("gmsh/model.jl")
include("gmsh/msh.jl")
include("gmsh/helpers.jl")

type FastCapResult
    capmat::Matrix{Float64}
    condnames::Vector{String}
    memory::Int64
    time::Float64
    log::String
end

end # module

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MPItrampoline"
version = v"1.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/eschnett/MPItrampoline/archive/d31bba9c3de5b37049bf3dc5343122dca264841a.tar.gz",
                  "ba1646970ae6923d35b5b3aed8018c4e129386955a4a9bea41a1846d9b23e98b"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd MPItrampoline-*
mkdir build
cd build
cmake \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DBUILD_SHARED_LIBS=ON \
    ..
cmake --build . --config RelWithDebInfo --parallel $nproc
cmake --build . --config RelWithDebInfo --parallel $nproc --target install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# Apple: Dynamically loaded libraries (`dlopen`) are not handled correctly.
# Windows: Does not have `dlopen`.
# musl: Does not define `RTLD_DEEPBIND` for `dlopen`.
platforms = filter(p -> !(Sys.isapple(p) || Sys.iswindows(p) || libc(p) == "musl"), platforms)
platforms = expand_gfortran_versions(platforms)
# libgfortran3 does not support `!GCC$ ATTRIBUTES NO_ARG_CHECK`. (We
# could in principle build without Fortran support there.)
platforms = filter(p -> libgfortran_version(p) ≠ v"3", platforms)

# TODO: on 32-bit systems, Fortran `MPI_Status` has the wrong size

# The products that we will ensure are always built
products = [
    ExecutableProduct("mpicc", :mpicc),
    ExecutableProduct("mpicxx", :mpicxx),
    ExecutableProduct("mpifc", :mpifc),
    ExecutableProduct("mpifort", :mpifort),
    ExecutableProduct("mpiexec", :mpiexec),

    LibraryProduct("libmpitrampoline", :libmpitrampoline),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
function linsolve_cr( LF::LF3dGrid, b::Array{Float64,1};
                      x0 = nothing,
                      NiterMax = 1000, TOL=5.e-10,
                      convmsg=false, showprogress=false )
  #
  Npoints = size(b)[1]
  #
  if x0 == nothing
    x = zeros(Float64,Npoints)
  else
    x = copy(x0)
  end
  #
  L_x = apply_Laplacian( LF, x )
  #
  r = b[:] - L_x[:]
  p = r[:]

  r_old = r[:]
  p_old = p[:]
  #
  L_p = apply_Laplacian( LF, p )
  x_old = x[:]
  #
  for iter = 1 : NiterMax
    #
    L_r = apply_Laplacian( LF, r )
    #
    rLr = dot( r, L_r )
    α = rLr / dot( L_p, L_p )  # ???
    #
    x_old = x[:]
    x[:] = x[:] + α * p[:]  # FIXME use x[:] to force x to be copied, not referenced
    #
    r_old = r[:]
    rLr_old = rLr
    r[:] = r[:] - α * L_p[:]
    #
    norm_res = sqrt( dot( r, r ) )
    norm_p   = sqrt( dot( p, p ) )
    Δx = abs(x-x_old)
    norm_Δx = sqrt( dot( Δx, Δx ) )
    if showprogress
      @printf("%8d %18.10e %18.10e %18.10e\n", iter, norm_res, norm_p, norm_Δx)
    end
    #
    if norm_res < TOL || norm_p < TOL || norm_Δx < TOL
      if convmsg
        @printf("#Convergence achieved in linsolve_cr: %8d iterations.\n", iter )
      end
      break
    end
    #
    L_r = apply_Laplacian( LF, r ) # calculate new L_r
    β = dot( r, L_r ) / rLr_old
    #
    p_old = p[:]
    p = r[:] + β * p[:]
    #
    L_p = L_r[:] + β * L_p[:]
  end
  #
  return x
  #
end # of function


text:
# ---
# title: Supervised DDVFA Example
# id: ddvfa_supervised
# cover: ../assets/ddvfa.png
# date: 2021-11-30
# author: "[Sasha Petrenko](https://github.com/AP6YC)"
# julia: 1.6
# description: This demo shows how to use DDVFA for simple supervised learning by clustering Iris samples and mapping the modules internal categories to the true labels.
# ---

# DDVFA is an unsupervised clustering algorithm by definition, but it can be adaptived for supervised learning by mapping the module's internal categories to the true labels.
# ART modules such as DDVFA can also be used in simple supervised mode where provided labels are used in place of internal incremental labels for the clusters, providing a method of assessing the clustering performance when labels are available.

# We begin with importing AdaptiveResonance for the ART modules and MLDatasets for some data utilities.
using AdaptiveResonance # ART
using MLDatasets        # Iris dataset
using MLDataUtils       # Shuffling and splitting
using Printf            # Formatted number printing

# We will download the Iris dataset for its small size and benchmark use for clustering algorithms.
Iris.download(i_accept_the_terms_of_use=true)
features, labels = Iris.features(), Iris.labels()

# Because the MLDatasets package gives us Iris labels as strings, we will use the `MLDataUtils.convertlabel` method with the `MLLabelUtils.LabelEnc.Indices` type to get a list of integers representing each class:
labels = convertlabel(LabelEnc.Indices{Int}, labels)
unique(labels)

# Next, we will create a train/test split with the `MLDataUtils.stratifiedobs` utility:
(X_train, y_train), (X_test, y_test) = stratifiedobs((features, labels))

# Now, we can create our DDVFA module.
# We'll do so with the default contstructor, though the module itself has many options that you can alter during instantiation.
art = DDVFA()

# We can train the model in batch mode upon the data in a simple supervised mode.
# We do so by passing the integer vector of labels to the training method with the simple keyword `y`.
# Just as in unsupervised training, we can extract the module's prescribed labels from the training method, which should match up to the training labels as we will see later.

## Train in simple supervised mode by passing the labels as a keyword argument.
y_hat_train = train!(art, X_train, y=y_train)
println("Training labels: ",  size(y_hat_train), " ", typeof(y_hat_train))

# We can classify the testing data to see how we generalize.
# At the same time, we can see the effect of getting the best-matching unit in the case of complete mismatch (see the docs on [Mismatch vs. BMU](@ref mismatch-bmu))

## Classify both ways
y_hat = AdaptiveResonance.classify(art, X_test)
y_hat_bmu = AdaptiveResonance.classify(art, X_test, get_bmu=true)

## Check the shape and type of the output labels
println("Testing labels: ",  size(y_hat), " ", typeof(y_hat))
println("Testing labels with bmu: ",  size(y_hat_bmu), " ", typeof(y_hat_bmu))

# Finally, we can calculate the performances (number correct over total) of the model upon all three regimes:
# 1. Training data
# 2. Testing data
# 2. Testing data with `get_bmu=true`

## Calculate performance on training data, testing data, and with get_bmu
perf_train = performance(y_hat_train, y_train)
perf_test = performance(y_hat, y_test)
perf_test_bmu = performance(y_hat_bmu, y_test)

## Format each performance number for comparison
@printf "Training performance: %.4f\n" perf_train
@printf "Testing performance: %.4f\n" perf_test
@printf "Best-matching unit testing performance: %.4f\n" perf_test_bmu


text:
insts = OrderedDict()
insts["minimal"] = [
    ((false, 2, 3, 4, 0, 0, 0, 0, 0),),
    ((false, 2, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 5, 3, 4, 0, 0, 0, 0, 0),),
    ((true, 5, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ]
insts["fast"] = [
    ((false, 5, 3, 4, 0, 0, 0, 0, 0),),
    ((false, 5, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 5, 3, 4, 0, 0.1, 0.1, 0, 0),),
    ((false, 3, 4, 5, 0, 0, 0, 0, 0),),
    ((false, 3, 4, 5, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 3, 4, 5, 0, 0.1, 0.1, 0, 0),),
    ((true, 5, 3, 4, 0, 0, 0, 0, 0),),
    ((true, 5, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 5, 3, 4, 0, 0.1, 0.1, 0, 0),),
    ((true, 3, 4, 5, 0, 0, 0, 0, 0),),
    ((true, 3, 4, 5, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 3, 4, 5, 0, 0.1, 0.1, 0, 0),),
    ((false, 15, 10, 20, 0, 0, 0, 0, 0),),
    ((false, 15, 10, 20, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 15, 10, 20, 0, 0.1, 0.1, 0, 0),),
    ((true, 15, 10, 20, 0, 0, 0, 0, 0),),
    ((true, 15, 10, 20, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 15, 10, 20, 0, 0.1, 0.1, 0, 0),),
    ((false, 100, 8, 12, 0, 0, 0, 0, 0),),
    ((false, 100, 8, 12, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 100, 8, 12, 0, 0.1, 0.1, 0, 0),),
    ((true, 100, 8, 12, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 100, 8, 12, 0, 0.1, 0.1, 0, 0),),
    ]
return (MatrixRegressionNative, insts)


text:
export ghost

# GHOST's node similarity
function lapevals(A::DenseMatrix)
    D = Diagonal(vec(sum(A,1)))
    S = Diagonal(sqrt.(1./D.diag))
    L = S * (D - A) * S
    LAPACK.syevr!('N', 'V', 'U', L, -0.01, 2.01, 0, 0, 1e-3)[1]
end

# Calculates topological signature of a node
# based on its local neiborhood's spectral signature
# It induces a local neighbood graph and calculates
# its spectra. It does so for local neighborhoods graphs created
# by inducing on nodes near v from hops r=1..k
function multiscale_spectra(G::SparseMatrixCSC,k::Integer,v,debug=false)
    evals = Vector{Vector{Float64}}(k)
    dists = bfs(G,v)[1]
    pf = find(x -> (x <= k) & (x >= 0), dists)
    G = G[pf,pf]
    dists = dists[pf]
    for r = 1:k
        p = find(x -> x <= r, dists) # calculate neighbors
        if length(p) > 1
            evals[r] = lapevals(full(G[p,p]))
        else
            evals[r] = Float64[]
        end
    end
    evals
end

function ghost_signature(G::SparseMatrixCSC,k::Integer,u,n=50,mode=:bins)
    evs = multiscale_spectra(G,k,u)
    bins = -1/n:(1/n):(2+1/n)
    if mode == :struct
        s = 0.01
        p = pdf.(Normal(0,s),linspace(-7.5s,7.5s,round(7.5s*n)))
    end
    map(evs) do x
        hx = Float64.(fit(Histogram, x, bins, closed=:right).weights .> 0)

        if mode == :struct
            resize!(hx, length(hx)+length(p)-1)
            hx[length(bins):end] = 0.0
            filt!(hx, p, 1, hx) #conv(p,cx)
        end

        nx = norm(hx,1)
        nx==0 ? hx : hx./nx
    end
end

function ghost_signature(G::SparseMatrixCSC,k::Integer)
    n = size(G,1)
    evs = Vector{Vector{Vector{Float64}}}(n)
    for v = 1:n
        print("\rG: $v/$n")
        evs[u] = ghost_signature(G,k,v)
    end
    println()
    evs
end


text:
function make_lists(d::Dict{Int,Vector{Int}})
    n = maximum(keys(d))
    allow = [d[k] for k = 1:n]
    return allow
end

"""
    iter_perms

This unexported method is the engine beneth `GenPerms`. The difference is this
returns a list of vectors, where each vector is a permutation of `1:n` based
on the table in `allow`.
"""
function iter_perms(allow::Vector{Vector{Int}}, p::Vector{Int} = Int[])
    n = length(allow)
    k = length(p)

    if n == k
        return [p]
    end

    valid_nexts = setdiff(allow[k+1], p)
    nested_its = Iterators.map((x) -> iter_perms(allow, [p; x]), valid_nexts)
    return Iterators.flatten(nested_its)
end

iter_perms(d::Dict{Int,Vector{Int}}) = iter_perms(make_lists(d))

export PermGen

"""
    PermGen

Create an iterator for permutations.
* `PermGen(n::Int)` creates an iterator for all permutations of length `n`.
* `PermGen(d::Dict{Int,Vector{Int}})` creates an iterator for all permutations in which
   element `k` may be a value in `d[k]`. Alternatively, `d` may be a `Vector` of `Vectors`.

Examples
--------
```
julia> for p in PermGen(3)
       println(p)
       end
(1)(2)(3)
(1)(2,3)
(1,2)(3)
(1,2,3)
(1,3,2)
(1,3)(2)

julia> d = [[2,3,4], [1,3,4], [1,2,4], [1,2,3]];

julia> for p in PermGen(d)
       println(p)
       end
(1,4)(2,3)
(1,3,2,4)
(1,2,3,4)
(1,4,2,3)
(1,3)(2,4)
(1,3,4,2)
(1,2,4,3)
(1,4,3,2)
(1,2)(3,4)
```
"""
function PermGen(allow::Vector{Vector{Int}})
    X = iter_perms(allow)
    return Iterators.map(Permutation, X)
end
PermGen(d::Dict{Int,Vector{Int}}) = PermGen(make_lists(d))

PermGen(n::Int) = Base.Generator(Permutation, permutations(1:n))

function deranged_allow(n::Int)
    d = Dict{Int,Vector{Int}}()
    base = collect(1:n)
    for k = 1:n
        d[k] = setdiff(base, k)
    end
    return d
end

function all_allow(n::Int)
    d = Dict{Int,Vector{Int}}()
    base = collect(1:n)
    for k = 1:n
        d[k] = base
    end
    return d
end

export DerangeGen

"""
    DerangeGen(n::Int)

Create an iterator for all derangements of `{1,2,...,n}`. These are
all permutations without any fixed points.
"""
DerangeGen(n::Int) = PermGen(deranged_allow(n))


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

function GitRepo(path::AbstractString)
    repo_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    err = ccall((:git_repository_open, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Cstring), repo_ptr_ptr, path)
    if err != Int(Error.GIT_OK)
        if repo_ptr_ptr[] != C_NULL
            finalize(GitRepo(repo_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return GitRepo(repo_ptr_ptr[])
end

function GitRepoExt(path::AbstractString, flags::Cuint = Cuint(Consts.REPOSITORY_OPEN_DEFAULT))
    separator = @static is_windows() ? ";" : ":"
    repo_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    err = ccall((:git_repository_open_ext, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Cstring, Cuint, Cstring),
                 repo_ptr_ptr, path, flags, separator)
    if err != Int(Error.GIT_OK)
        if repo_ptr_ptr[] != C_NULL
            finalize(GitRepo(repo_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return GitRepo(repo_ptr_ptr[])
end

function cleanup(r::GitRepo)
    if r.ptr != C_NULL
        ccall((:git_repository__cleanup, :libgit2), Void, (Ptr{Void},), r.ptr)
    end
end

function init(path::AbstractString, bare::Bool=false)
    repo_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_repository_init, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Cstring, Cuint), repo_ptr_ptr, path, bare)
    return GitRepo(repo_ptr_ptr[])
end

function head_oid(repo::GitRepo)
    head_ref = head(repo)
    try
        return Oid(head_ref)
    finally
        finalize(head_ref)
    end
end

function headname(repo::GitRepo)
    with(head(repo)) do href
        if isattached(repo)
            shortname(href)
        else
            "(detached from $(string(Oid(href))[1:7]))"
        end
    end
end

function isbare(repo::GitRepo)
    return ccall((:git_repository_is_bare, :libgit2), Cint, (Ptr{Void},), repo.ptr) == 1
end

function isattached(repo::GitRepo)
    ccall((:git_repository_head_detached, :libgit2), Cint, (Ptr{Void},), repo.ptr) != 1
end

""" Returns a found object """
function revparse(repo::GitRepo, objname::AbstractString)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    err = ccall((:git_revparse_single, :libgit2), Cint,
            (Ptr{Ptr{Void}}, Ptr{Void}, Cstring), obj_ptr_ptr, repo.ptr, objname)
    err != 0 && return nothing
    return GitAnyObject(obj_ptr_ptr[])
end

""" Returns id of a found object """
function revparseid(repo::GitRepo, objname::AbstractString)
    obj = revparse(repo, objname)
    obj === nothing && return Oid()
    oid = Oid(obj.ptr)
    finalize(obj)
    return oid
end

function get{T <: GitObject}(::Type{T}, r::GitRepo, oid::Oid, oid_size::Int=OID_HEXSZ)
    id_ptr  = Ref(oid)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    git_otype = getobjecttype(T)

    err = if oid_size != OID_HEXSZ
        ccall((:git_object_lookup_prefix, :libgit2), Cint,
              (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Oid}, Csize_t, Cint),
              obj_ptr_ptr, r.ptr, id_ptr, Csize_t(oid_size), git_otype)
    else
        ccall((:git_object_lookup, :libgit2), Cint,
              (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Oid}, Cint),
              obj_ptr_ptr, r.ptr, id_ptr, git_otype)
    end
    if err == Int(Error.ENOTFOUND)
        return nothing
    elseif err != Int(Error.GIT_OK)
        if obj_ptr_ptr[] != C_NULL
            finalize(GitAnyObject(obj_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return T(obj_ptr_ptr[])
end

function get{T <: GitObject}(::Type{T}, r::GitRepo, oid::AbstractString)
    return get(T, r, Oid(oid), length(oid))
end

function gitdir(repo::GitRepo)
    return unsafe_string(ccall((:git_repository_path, :libgit2), Cstring,
                        (Ptr{Void},), repo.ptr))
end

function path(repo::GitRepo)
    rpath = gitdir(repo)
    return isbare(repo) ? rpath : splitdir(rpath[1:end-1])[1]*"/" # remove '.git' part
end

function peel(obj::GitObject, obj_type::Cint)
    peeled_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    git_otype = getobjecttype(obj_type)
    err = ccall((:git_object_peel, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cint), peeled_ptr_ptr, obj.ptr, obj_type)
    if err == Int(Error.ENOTFOUND)
        return Oid()
    elseif err != Int(Error.GIT_OK)
        if peeled_ptr_ptr[] != C_NULL
            finalize(GitAnyObject(peeled_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return git_otype(peeled_ptr_ptr[])
end

peel{T <: GitObject}(::Type{T}, obj::GitObject) = peel(obj, getobjecttype(T))

function checkout_tree(repo::GitRepo, obj::GitObject;
                       options::CheckoutOptions = CheckoutOptions())
    @check ccall((:git_checkout_tree, :libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}, Ptr{CheckoutOptions}),
                 repo.ptr, obj.ptr, Ref(options))
end

function checkout_index(repo::GitRepo, idx::Nullable{GitIndex} = Nullable{GitIndex}();
                        options::CheckoutOptions = CheckoutOptions())
    @check ccall((:git_checkout_index, :libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}, Ptr{CheckoutOptions}),
                 repo.ptr,
                 isnull(idx) ? C_NULL : Base.get(idx).ptr,
                 Ref(options))
end

function checkout_head(repo::GitRepo; options::CheckoutOptions = CheckoutOptions())
    @check ccall((:git_checkout_head, :libgit2), Cint,
                 (Ptr{Void}, Ptr{CheckoutOptions}),
                 repo.ptr, Ref(options))
end

"""Updates some entries, determined by the `pathspecs`, in the index from the target commit tree."""
function reset!{T<:AbstractString, S<:GitObject}(repo::GitRepo, obj::Nullable{S}, pathspecs::T...)
    with(StrArrayStruct(pathspecs...)) do sa
        @check ccall((:git_reset_default, :libgit2), Cint,
                (Ptr{Void}, Ptr{Void}, Ptr{StrArrayStruct}),
                repo.ptr,
                isnull(obj) ? C_NULL: Base.get(obj).ptr,
                Ref(sa))
    end
end

"""Sets the current head to the specified commit oid and optionally resets the index and working tree to match."""
function reset!(repo::GitRepo, obj::GitObject, mode::Cint;
               checkout_opts::CheckoutOptions = CheckoutOptions())
    @check ccall((:git_reset, :libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}, Cint, Ptr{CheckoutOptions}),
                  repo.ptr, obj.ptr, mode, Ref(checkout_opts))
end

function clone(repo_url::AbstractString, repo_path::AbstractString,
               clone_opts::CloneOptions)
    clone_opts_ref = Ref(clone_opts)
    repo_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_clone, :libgit2), Cint,
            (Ptr{Ptr{Void}}, Cstring, Cstring, Ref{CloneOptions}),
            repo_ptr_ptr, repo_url, repo_path, clone_opts_ref)
    return GitRepo(repo_ptr_ptr[])
end

function fetchheads(repo::GitRepo)
    fhr = Ref{Vector{FetchHead}}(FetchHead[])
    ffcb = fetchhead_foreach_cb()
    @check ccall((:git_repository_fetchhead_foreach, :libgit2), Cint,
                  (Ptr{Void}, Ptr{Void}, Ptr{Void}),
                   repo.ptr, ffcb, fhr)
    return fhr[]
end

function remotes(repo::GitRepo)
    out = Ref(StrArrayStruct())
    @check ccall((:git_remote_list, :libgit2), Cint,
                  (Ptr{Void}, Ptr{Void}), out, repo.ptr)
    return convert(Vector{AbstractString}, out[])
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# Sample postprocessing script for customized figure.
#
# This script creates contours from a series of input files with a fixed color range.
# Combined with ffmpeg, we can easily make animations from data.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Glob, PyPlot, Printf

files = glob("bulk*.vlsv")
nfile = length(files)
# Set output directory
outdir = "out/"

meta = load(files[1])
vardict = Dict("v"=>"proton/vg_v", "rho"=>"proton/vg_rho", "b"=>"vg_b_vol", "b2"=>"fg_b")
varname = "rho"

fig, ax = plt.subplots()

comp = :z # vector component for plotting (if applicable)
axisunit = EARTH
colorscale = Log
addcolorbar = true
# Choose colormap
cmap = matplotlib.cm.turbo
# Set data plotting range
vmin = 7.0e4
vmax = 2.5e6

pArgs = Vlasiator.set_args(meta, vardict[varname], axisunit; normal=:none)

x1, x2 = Vlasiator.get_axis(pArgs)

norm, ticks = Vlasiator.set_colorbar(colorscale, vmin, vmax)

fakedata = zeros(Float32, length(x2), length(x1))

c = ax.pcolormesh(x1, x2, fakedata; norm, cmap)

Vlasiator.set_plot(c, ax, pArgs, ticks, addcolorbar)

for (i, file) in enumerate(files)
   @info "$i out of $nfile"
   local meta = load(file)

   var = meta[vardict[varname]]
   t = readparameter(meta, "time")

   data = Vlasiator.prep2d(meta, vardict[varname], comp)'
   c.set_array(data)

   str_title = @sprintf "t= %4.1fs" t
   ax.set_title(str_title)

   savefig(outdir*lpad(i, 4, '0')*".png", bbox_inches="tight", dpi=100)
end

close()


text:
zero(::Type{ScalarExpr}) = scalar(0)

qsparse(arr::AbstractArray) = sparse(arr)
qsparse(arr::Array{InnerProduct}) = sparse(convert(Array{ScalarExpr}, arr))
qsparse(arr::Array{Union(InnerProduct, Int)}) = sparse(convert(Array{ScalarExpr}, arr))
qsparse(arr::Array{Union(Int, InnerProduct)}) = sparse(convert(Array{ScalarExpr}, arr))

showsp(x) = showsp(STDOUT::IO, x)

orient_error(a,b) = error("Multiplication $(typeof(a))*$(typeof(b)) is undefined. Perhaps you meant to use kron($(typeof(a)), $(typeof(b)))?")

samebasis(a::Dirac, b::Dirac)= bsym(a)==bsym(b)
samebasis(a, b::Dirac) = a==bsym(b)
samebasis(a::Dirac, b) = samebasis(b, a)

for t=(:DiracMatrix, :DiracVector, :State, :OuterProduct)
	@eval +(n::Number, d::($t)) = error("cannot add number and $(string($t))")
	@eval +(d::($t), n::Number) = +(n,d)
	@eval -(n::Number, d::($t)) = +(n,d)
	@eval -(d::($t), n::Number) = +(n,d)
end

*(a::Dirac, b::Dirac) = orient_error(a,b)
*(a, b::Dirac) = orient_error(a,b)
*(a::Dirac, b) = orient_error(a,b)

kron(a::AbstractScalar, b::AbstractScalar) = a*b
kron(a::Dirac, b::DiracCoeff) = a*b
kron(a::DiracCoeff, b::Dirac) = b*a

convert{A<:Array}(::Type{A}, d::DiracVector) = convert(A, vec(full(d.coeffs)))
convert{A<:Array}(::Type{A}, d::DiracMatrix) = convert(A, full(d.coeffs))
convert{A<:AbstractArray}(::Type{A}, d::Dirac) = convert(A, d.coeffs)

promote_rule{A<:AbstractArray, D<:Dirac}(::Type{A}, ::Type{D}) = A
#arithmetic functrions for generic arrays
for op in (:+, :-, :*, :.+, :.-, :.*, :./, :^, :.^, :kron)
	@eval begin
	($op)(a::AbstractArray, d::Dirac) = ($op)(a, d.coeffs)
	($op)(d::Dirac, a::AbstractArray) = ($op)(d.coeffs, a)
	end
end


text:
module OmegaGrad

using Reexport

include("gradinterface.jl")

include("reversediff.jl")
@reexport using .OmegaReverseDiff

include("zygote.jl")
@reexport using .OmegaZygote

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
#!/usr/bin/env julia

module HealPy

export hp

using PyCall

const hp = PyNULL()

function __init__()
    if Sys.iswindows()
        @warn "HealPy is not supported on Windows."
    else
        copy!(hp, pyimport_conda("healpy", "healpy", "conda-forge"))
    end
end

end

# vim: set sw=4 et sts=4 :


text:
# Default parameters

default(k::Symbol) = get(_defaults, k, nothing)
default() = copy(_defaults)
default!(k::Symbol, v) = (_defaults[k] = v)
default!(d::Dict{Symbol,Any}) = (_defaults = copy(d); nothing)
default!(v::Superfluid) = default!(:superfluid, v)
default!(v::Discretisation) = default!(:discretisation, v)

# TODO extend Documenter to format docstrings for default symbols, and
# warn about symbols in the default table that are not documented.

# TODO ask why optional positional arguments must come at the end.
# (Kind of similar to why Julia doesn't add multiple dispatch when it already handles ambiguity.)

# Abstract syntax types

struct Xpn{S} <: AbstractVector{Any}
    args
end

const MethodDefinition = Union{Xpn{:(=)}, Xpn{:function}}

# Generate default arguments
#
# If the expression m is a call, define the function and methods that allow
# initial arguments of types Superfluid and Discretisation to be elided.
#
# If m is a method definition, keyword arguments with no explicit
# default are looked up in defaults

macro defaults(E) defaults(xpn(E)) |> esc end

defaults(E::Xpn{:call}) =
    quote
        function $(called_function(E)) end
        $(default_methods(E)...)
    end

defaults(E::MethodDefinition) = default_kwargs(E)

# The whole generated code is escaped from hygene, so default is
# inserted as a function, not as a symbol to be evaluated in the
# module where @defaults is expanded.

default_methods(E) = default_methods(E, method_types(E))
default_methods(E, ts) = [elided(E, e) for e in power_set(elisions(ts...)...)]

elisions(::Type{Superfluid}, _...) = [1=>:superfluid]
elisions(::Type{Discretisation}, _...) = [1=>:discretisation]
elisions(::Type{Superfluid}, ::Type{Discretisation}, _...) =
    [1=>:superfluid, 2=>:discretisation]
elisions(_, t::Union{Type{Superfluid}, Type{Discretisation}}, ts...) =
    [(j+1)=>s for (j,s) in elisions(t, ts...)]

function elided(E, ds)
    # TODO convert xs::Vararg{T} to xs...
    f(E::Symbol) = E
    f(E::Xpn{:(::)}) =
        if E[2] isa Xpn{:curly} && E[2][1] == :Vararg
            :( $(E[1])... ) |> xpn
        else
            E[1]
        end
    f(E::Xpn{:(...)}) = E
    f(E::Xpn{:kw}) = f(E[1])
    as = args(E)
    bs = Any[f(a) for a in as]
    for d in ds
        bs[d[1]] = :( $default($(QuoteNode(d[2]))) ) |> xpn
    end
    as = as[setdiff(eachindex(as), first.(ds))]
    Xpn{:function}([method_call(E, as), method_call(E, bs)]) |> expr
end

power_set() = []
power_set(a) = [[a]]
power_set(a,b) = [[a], [b], [a,b]]

"Add defaults for keyword arguments"
function default_kwargs(E)
    as, ps = aps(E)
    Xpn{:function}([method_call(E, as, add_defaults(ps)), body(E)]) |> expr
end

# This doesn't handle typed keyword parameters
add_defaults(ps) = [
    (p isa Symbol && haskey(_defaults, p)) ?
    Xpn{:kw}([p, xpn(:( $default($(QuoteNode(p))) ))]) :
    p
    for p in ps]

# Abstract syntax

xpn(E) = E
xpn(E::Expr) = Xpn(E)
Xpn(E::Expr) = Xpn{E.head}(xpn.(E.args))
Base.size(E::Xpn) = size(E.args)
Base.getindex(E::Xpn, i) = E.args[i]
expr(E) = E
expr(E::Xpn{S}) where S = Expr(S, expr.(E)...)

called_function(E::Xpn{:call}) = E[1]
called_function(E::MethodDefinition) = called_function(E[1])
body(E::MethodDefinition) = E[2]

args(E::Xpn) = aps(E)[1]
kwargs(E::Xpn) = aps(E)[2]
aps(E::MethodDefinition) = aps(E[1])
function aps(E::Xpn{:call})
    arg_list = E[2:end]
    if arg_list[1] isa Xpn{:parameters}
        arg_list[2:end], collect(arg_list[1])
    else
        arg_list, Xpn[]
    end
end

method_call(E::MethodDefinition, xs...) = method_call(E[1], xs...)
method_call(E::Xpn{:call}, args, kwargs=[Xpn(:(kwargs...))]) =
    :( $(called_function(E))(
        $(expr.(args)...);
        $(expr.(kwargs)...))) |> Xpn

function method_types(m::Xpn)
    t(x::Symbol) = :Any
    t(E::Xpn{:(::)}) = E[2]
    t(E::Xpn{:(...)}) = t(only(E))
    t(E::Xpn{:kw}) = t(E[1])
    # In principle, type identifiers should be evaluated in the module
    # where @defaults is expanded.  In practice that is always Superfluids.
    eval.(t.(args(m)))
end

_defaults = Dict{Symbol,Any}(
    :hbm => 1,
    :g_tol => 1e-6,
    :iterations => 1000,
    :dt => 1e-4,
    :relaxer => Optim.ConjugateGradient,
    :integrator => DifferentialEquations.RK4,
    :diagonizer => nothing,# eigenproblem solver
    :xlims => nothing,
    :ylims => nothing,
)

const OPT_ARGS = (:g_tol, :iterations, :relaxer)
const DIFEQ_ARGS = (:dt, :formula)


text:
"""
    Parameter(type::String, value)

SQL query parameter, where `type` is one of the Druid SQL data types and `value`
is the parameter value.
"""
struct Parameter
    type::String
    value
    function Parameter(type, value)
        type ∈ [
            "CHAR", "VARCHAR",
            "DECIMAL", "FLOAT", "REAL", "DOUBLE",
            "BOOLEAN", "TINYINT", "SMALLINT", "INTEGER", "BIGINT",
            "TIMESTAMP", "DATE", "OTHER"
        ] || error("Invalid 'type'")
        new(type, value)
    end
end
JSON.lower(p::Parameter) = non_nothing_dict(p)

"""
    Sql(query::String; <keyword arguments>)

A method with all arguments as keyword arguments is also provided.

# Arguments
- parameters::Vector{Parameter} = nothing
- resultFormat::String = nothing
- header::Bool = true
- context::Dict = nothing
"""
mutable struct Sql <: Query
    query::String
    parameters
    resultFormat
    header
    context
    function Sql(query; parameters=nothing, resultFormat=nothing, header=nothing, context=nothing)
        nothing_or_type(parameters, Vector{Parameter})
        resultFormat === nothing || resultFormat ∈ ["object", "array", "objectLines", "arrayLines", "csv"] || error("Invalid resultFormat")
        nothing_or_type(header, Bool)
        nothing_or_type(context, Dict)
        new(query, parameters, resultFormat, header, context)
    end
end
Sql(; query, parameters=nothing, resultFormat=nothing, header=nothing, context=nothing) =
    Sql(query; parameters, resultFormat, header, context)

query_type(::Sql) = "SQL"

struct SqlResult{ResultFormat} <: QueryResult
    names::Vector{Symbol}
    inner_data::ResultFormat
end

Tables.istable(::SqlResult) = false

struct SqlRow{T} <: Tables.AbstractRow
    row::Int
    source::SqlResult{T}
end

# resultFormat = "object", "objectLines"
Tables.istable(::SqlResult{Vector{Dict{String, T}}}) where {T} = true
names(sr::SqlResult{Vector{Dict{String, T}}}) where {T} = getfield(sr, :names)

Base.summary(io::IO, sr::SqlResult{Vector{Dict{String, T}}}) where {T} = print(io, string(length(sr)) * "-element " * string(typeof(sr)))
Base.getindex(sr::SqlResult{Vector{Dict{String, T}}}, i::Int) where {T} = (i <= length(sr) || throw(BoundsError(sr, i))) && SqlRow(i, sr)

Tables.rowaccess(::SqlResult{Vector{Dict{String, T}}}) where {T} = true
Tables.rows(sr::SqlResult{Vector{Dict{String, T}}}) where {T} = sr

Base.eltype(::SqlResult{Vector{Dict{String, T}}}) where {T} = SqlRow{Vector{Dict{String, T}}}
Base.length(sr::SqlResult{Vector{Dict{String, T}}}) where {T} = length(getfield(sr, :inner_data))
Base.iterate(sr::SqlResult{Vector{Dict{String, T}}}, state=1) where {T} = state > length(sr) ? nothing : (SqlRow(state, sr), state + 1)

function Tables.schema(sr::SqlResult{Vector{Dict{String, T}}}) where {T}
    if length(sr) != 0
        types = [typeof(cl) for cl in sr[1]]
    else
        types = DataType[]
    end
    Tables.Schema(names(sr), types)
end

Tables.getcolumn(sr::SqlRow{Vector{Dict{String, T}}}, name::Symbol) where {T} = getfield(getfield(sr, :source), :inner_data)[getfield(sr, :row)][string(name)]
Tables.getcolumn(sr::SqlRow{Vector{Dict{String, T}}}, i::Int) where {T} = Tables.getcolumn(sr, names(getfield(sr, :source))[i])
Tables.columnnames(sr::SqlRow{Vector{Dict{String, T}}}) where {T} = names(getfield(sr, :source))

function transform_result(res, header, ::Val{:object}; names=:auto)
    local interm
    try
        interm = JSON.parse(res)
    catch e
        println("Corrupted response received from Druid, unable to parse")
        rethrow(e)
    end
    d = convert(Vector{Dict{String, Any}}, interm)
    if names == :auto
        if header == true
            names = Symbol.(keys(popfirst!(d)))
        elseif length(d) != 0
            names = Symbol.(keys(d[1]))
        else
            names = Symbol[]
        end
    end
    SqlResult(names, d)
end

function transform_result(res, header, ::Val{:objectLines}; names=:auto)
    length(res) >= 2 && res[end - 1:end] == "\n\n" || error("Corrupted response received from Druid, unable to parse")
    interm = strip(res)
    if interm == ""
        d = Dict{String, Any}[]
    else
        d = JSON.parse.(split(interm, '\n'))
    end
    if names == :auto
        if header == true
            names = Symbol.(keys(popfirst!(d)))
        elseif length(d) != 0
            names = Symbol.(keys(d[1]))
        else
            names = Symbol[]
        end
    end
    SqlResult(names, d)
end

# resultFormat = "array", "arrayLines"
Tables.istable(::SqlResult{Matrix{T}}) where {T} = true
names(sr::SqlResult{Matrix{T}}) where {T} = getfield(sr, :names)

Base.summary(io::IO, sr::SqlResult{Matrix{T}}) where {T} = print(io, string(length(sr)) * "-element " * string(typeof(sr)))
Base.getindex(sr::SqlResult{Matrix{T}}, i::Int) where {T} = (i <= length(sr) || throw(BoundsError(sr, i))) && SqlRow(i, sr)

Tables.rowaccess(::SqlResult{Matrix{T}}) where {T} = true
Tables.rows(sr::SqlResult{Matrix{T}}) where {T} = sr

Base.eltype(::SqlResult{T}) where {T<:Matrix} = SqlRow{T}
Base.length(sr::SqlResult{Matrix{T}}) where {T} = size(getfield(sr, :inner_data), 1)
Base.iterate(sr::SqlResult{Matrix{T}}, state=1) where {T} = state > length(sr) ? nothing : (SqlRow(state, sr), state + 1)

function Tables.schema(sr::SqlResult{Matrix{T}}) where {T}
    if length(sr) != 0
        types = [typeof(cl) for cl in sr[1]]
    else
        types = DataType[]
    end
    Tables.Schema(names(sr), types)
end

Tables.getcolumn(sr::SqlRow{Matrix{T}}, i::Int) where {T} = getfield(getfield(sr, :source), :inner_data)[getfield(sr, :row), i]
Tables.getcolumn(sr::SqlRow{Matrix{T}}, name::Symbol) where {T} = Tables.getcolumn(sr, findfirst(name .== names(getfield(sr, :source))))
Tables.columnnames(sr::SqlRow{Matrix{T}}) where {T} = names(getfield(sr, :source))

function transform_result(res, header, ::Val{:array}; names=:auto)
    local interm
    try
        interm = JSON.parse(res)
    catch e
        println("Corrupted response received from Druid, unable to parse")
        rethrow(e)
    end
    if names == :auto
        if header == true
            names = Symbol.(popfirst!(interm))
        elseif length(interm) != 0
            names = [Symbol("Column" * string(i)) for i ∈ 1:length(interm[1])]
        else
            names = Symbol[]
        end
    end
    if interm == []
        d = Matrix{Any}(undef, 0, 0)
    else
        d = permutedims(reduce(hcat, interm))
    end
    SqlResult(names, d)
end

function transform_result(res, header, ::Val{:arrayLines}; names=:auto)
    length(res) >= 2 && res[end - 1:end] == "\n\n" || error("Corrupted response received from Druid, unable to parse")
    stripped = strip(res)
    if stripped == ""
        interm = Vector{Any}[]
    else
        interm = JSON.parse.(split(stripped, '\n'))
    end
    if names == :auto
        if header == true
            names = Symbol.(popfirst!(interm))
        elseif length(interm) != 0
            names = [Symbol("Column" * string(i)) for i ∈ 1:length(interm[1])]
        else
            names = Symbol[]
        end
    end
    if interm == []
        d = Matrix{Any}(undef, 0, 0)
    else
        d = permutedims(reduce(hcat, interm))
    end
    SqlResult(names, d)
end

# resultFormat = "csv"
function transform_result(res, header, ::Val{:csv}; names=:auto)
    (length(res) >= 2 && res[end - 1:end] == "\n\n") ||
        (length(res) == 1 && res[end] == '\n') || error("Corrupted response received from Druid, unable to parse")
    if names == :auto
        if header == true
            return CSV.File(IOBuffer(res); header)
        else
            return CSV.File(IOBuffer(res); header=false)
        end
    else
        return CSV.File(IOBuffer(res); header=names)
    end
end

"""
    execute(client, query::Sql; names=:auto)

Executes the Druid SQL `query` on the `client::Client`. Returns the query
results as a String. Throws exception if query execution fails.

A Vector{Symbol} can be provided for the column names through `names`.

# Examples
```julia-repl
julia> query = Sql("SELECT * FROM some_datasource LIMIT 10")
julia> execute(client, query)
"[...\\n]"

julia> query = Sql("SELECT * FROM non_existent_datasource")
julia> execute(client, query)
ERROR: Druid error during query execution
Dict{String,Any}("errorClass" => "org.apache.calcite.tools.ValidationException"...
```
"""
function execute(client::Client, query::Sql, names=:auto)
    names == :auto || typeassert(names, Vector{Symbol})
    url = joinpath(client.url, "druid/v2/sql")
    header = query.header == true
    if query.resultFormat === nothing
        resultFormat = "object"
    else
        resultFormat = query.resultFormat
    end
    transform_result(do_query(url, query), header, Val(Symbol(resultFormat)); names=names)
end


text:
### input and output functions for jFlow.jl

using DelimitedFiles

function ReadMaterials()::Array{Material, 1}
    # read various model parameters from file
    material = Material[]
    data = readdlm("Materials.txt", '\t', header=true)
    for i = 1:size(data[1], 1)
        index = Int64(data[1][i, 1])
        Kh = Float64(data[1][i, 2])                 # horizontal and vertical hydraulic conductivities
        Kz = Float64(data[1][i, 3])
        Ss = Float64(data[1][i, 4])                 # specific storage
        phi = Float64(data[1][i, 5])                # porosity
        alpha = Float64(data[1][i, 6])              # Van Genuchten unsaturated hydraulic parameters
        N = Float64(data[1][i, 7])
        m = 1.0 - 1.0/N
        push!(material, Material(index, Kh, Kz, Ss, phi, alpha, N, m))
    end
    println("Read material properties.")
    return material
end

function ReadParams()::Params
    # read numerical model parameters
    data = readdlm("Params.txt", '\t', header=false)
    endTime = Float64(data[1, 2])                   # end of simulation
    dt0 = Float64(data[2, 2])                       # starting time step
    dtMax = Float64(data[3, 2])                     # maximum time step
    dhMax = Float64(data[4, 2])                     # maximum head change (anywhere) per time step
    f = Float64(data[5, 2])                         # time step adjustment factor
    c = Float64(data[6, 2])                         # explicit time-step weighting (0=explicit, 0.5=Crank-Nicholson, 1=fully-implicit)
    mesh = Bool(data[7, 2])                         # externally generated mesh (if == true)
    params = Params(endTime, dt0, dtMax, dhMax, f, c, mesh)
    println("Read model run parameters.")
    return params
end

function ReadCells1(material::Array{Material, 1})::Array{Cell, 1}
    # read grid properties from file generated by external grid generator
    cell = Cell[]
    data = readdlm("nodes.csv", ',', header=true)
    for iRow = 1:size(data[1], 1)
        cellIndex = Int64(data[1][iRow, 1])
        x = Float64(data[1][iRow, 2])               # location of volume element/cell node point
        y = Float64(data[1][iRow, 3])
        z = Float64(data[1][iRow, 4])
        vol = Float64(data[1][iRow, 5])             # cell volume
        b = Float64(data[1][iRow, 6])               # cell vertical thickness
        h = Float64(data[1][iRow, 7])               # hydraulic head (an initial condition, at start up)
        matIndex = Int64(data[1][iRow, 8])          # index number for material
        mat = material[matIndex]
        connects = Int64[]
        neighbors = Int64[]
        push!(cell, Cell(cellIndex, x, y, z, vol, b, h-z, 0.0, mat, connects, neighbors))
        cell[end].S = Sat(cell[end])                # calculate initial saturation
    end
    println("Read and process externally-generated node assembly.")
    return cell
end

function ReadCells2(material::Array{Material, 1})::Array{Cell, 1}
    # read cell setup file and process cell objects from raw constraints
    cell = Cell[]
    data = readdlm("Cells.txt", '\t', header=true)
    for iRow = 1:size(data[1], 1)
        cellIndex = Int64(data[1][iRow, 1])
        x0 = Float64(data[1][iRow, 2])              # location of volume element/cell node point
        y0 = Float64(data[1][iRow, 3])
        z0 = Float64(data[1][iRow, 4])
        vol = Float64(data[1][iRow, 5])             # cell volume; note: set equal to very large number for fixed-head condition
        b = Float64(data[1][iRow, 6])               # cell vertical thickness
        h = Float64(data[1][iRow, 7])               # hydraulic head (an initial condition, at start up)
        matIndex = Int64(data[1][iRow, 8])          # index number for material
        mat = material[matIndex]
        xNum = Int64(data[1][iRow, 9])              # number of identical cells along different axes
        yNum = Int64(data[1][iRow, 10])
        zNum = Int64(data[1][iRow, 11])
        xInc = Int64(data[1][iRow, 12])             # increment number for cellIndex while populating along different axes
        yInc = Int64(data[1][iRow, 13])
        zInc = Int64(data[1][iRow, 14])
        xStep = Float64(data[1][iRow, 15])          # change in node point coordinate while populating along different axes
        yStep = Float64(data[1][iRow, 16])
        zStep = Float64(data[1][iRow, 17])
        for k = 1:zNum, j = 1:yNum, i = 1:xNum
            index = cellIndex + (i-1)*xInc + (j-1)*yInc + (k-1)*zInc
            x = x0 + (i-1)*xStep
            y = y0 + (j-1)*yStep
            z = z0 + (k-1)*zStep
            connects = Int64[]
            neighbors = Int64[]
            push!(cell, Cell(index, x, y, z, vol, b, h-z, 0.0, mat, connects, neighbors))
            cell[end].S = Sat(cell[end])        # calculate initial saturation, now that cell has been defined
        end
    end
    println("Read and assigned cell properties.")
    return cell
end

function ReadCellMods(cell::Array{Cell, 1}, material::Array{Material, 1})::Array{Cell, 1}
    # modify select cell properties, if specified
    data = readdlm("CellMods.txt", '\t', header=true)
    for i = 1:size(data[1], 1)
        cellIndex = Int64(data[1][i, 1])
        vol = Float64(data[1][i, 2])
        h0 = Float64(data[1][i, 3])
        matIndex = Int64(data[1][i, 4])
        cell[cellIndex].vol = vol
        cell[cellIndex].P = h0 - cell[cellIndex].z
        cell[cellIndex].mat = material[matIndex]
    end
    println("Modified subset of cells, as specified.")
    return cell
end

function ReadGridConnects1()::Array{Connection, 1}
    # read grid cell connection parameters from file generated by external grid generator
    connect = Connection[]
    data = readdlm("connections.csv", ',', header=true)
    for iRow = 1:size(data[1], 1)
        cell1 = Int64(data[1][iRow, 1])                         # connected cells index numbers
        cell2 = Int64(data[1][iRow, 2])
        del1 = Float64(data[1][iRow, 3])                        # distances from cell centerpoints to interface
        del2 = Float64(data[1][iRow, 4])
        A = Float64(data[1][iRow, 5])                           # connection interface area
        zConnect = Bool(data[1][iRow, 6])                       # is this a vertical connection
        push!(connect, Connection(cell1, cell2, del1, del2, A, zConnect, 0.0))
    end
    println("Read and process externally-generated node assembly.")
    return connect
end

function ReadGridConnects2()::Array{Connection, 1}
    # read and assign connections, formatted for internal grid cells
    connect = Connection[]
    data = readdlm("gridConnections.txt", '\t', header=true)
    for iRow = 1:size(data[1], 1)
        cellStart = Int64(data[1][iRow, 1])                 # starting cell index for connection sequence
        numCellStart = Int64(data[1][iRow, 2])              # number of start cells (along row, col, etc.)
        stepCellStart = Int64(data[1][iRow, 3])             # increment between start cells
        numPopulate = Int64(data[1][iRow, 4])               # number of connections in direction of population
        stepPopulate = Int64(data[1][iRow, 5])              # increment between populated cells
        del1 = Float64(data[1][iRow, 6])                    # node point-to-interface distances and interface area
        del2 = Float64(data[1][iRow, 7])
        A = Float64(data[1][iRow, 8])
        zConnect = Bool(data[1][iRow, 9])                   # boolean flag indicating if connect is oriented vertically
        for j = 1:numCellStart, i = 1:numPopulate
            cell1 = cellStart + (j-1)*stepCellStart + (i-1)*stepPopulate
            cell2 = cellStart + (j-1)*stepCellStart + i*stepPopulate
            push!(connect, Connection(cell1, cell2, del1, del2, A, zConnect, 0.0))
        end
    end
    println("Read and assigned internal grid cell-to-cell connections.")
    return connect
end

function ReadSpecConnects(connect::Array{Connection, 1})::Array{Connection, 1}
    # read and assign special cell connections for boundaries, unique cells, etc.
    data = readdlm("specConnections.txt", '\t', header=true)
    for iRow = 1:size(data[1], 1)
        cellStart = Int64(data[1][iRow, 1])                 # starting cell index for connection sequence
        numCellStart = Int64(data[1][iRow, 2])              # number of start cells (along row, col, etc.)
        stepCellStart = Int64(data[1][iRow, 3])             # increment between start cells
        boundaryCell = Int64(data[1][iRow, 4])              # index number of cell completing the connection
        del1 = Float64(data[1][iRow, 5])
        del2 = Float64(data[1][iRow, 6])
        A = Float64(data[1][iRow, 7])
        zConnect = Bool(data[1][iRow, 8])
        for j = 1:numCellStart
            cell1 = cellStart + (j-1)*stepCellStart
            cell2 = boundaryCell
            push!(connect, Connection(cell1, cell2, del1, del2, A, zConnect, 0.0))
        end
    end
    println("Read and assigned special cell-to-cell connections.")
    return connect
end

function ReadMonitors()::Array{Int64, 1}
    # return a list of monitoring point (cell) indices
    data = readdlm("Monitors.txt", '\t', header=true)
    monitor = Int.(data[1])
    monitor = vec(monitor)
    return monitor
end

function ReadSources()
    # read source/sinks terms and assign to cells & stress periods
    data = readdlm("Sources.txt", '\t', header=true)
    stressPeriods = sort(collect(Set(float(data[1][:, 2]))))
    sourceCells = sort(collect(Set(data[1][:, 1])))
    sourceCells = floor.(Int, sourceCells)
    source = zeros(Float64, length(sourceCells), length(stressPeriods))
    for i = 1:size(data[1], 1)
        cellIndex = Int64(data[1][i, 1])
        t0 = Float64(data[1][i, 2])             # starting time for source
        Q = Float64(data[1][i, 3])              # volumetric flux (L^3/T)
        iCol = FindStressPeriod(t0, stressPeriods)
        iRow = findall(sourceCells.==cellIndex)[1]
        source[iRow, iCol:end].= Q
    end
    println("Processed source terms.")
    return stressPeriods, sourceCells, source
end

function ReadDrains()
    # read drain-designated cells and associated threshold pressures
    data = readdlm("Drains.txt", '\t', header=true)
	drainCell = Int64.(data[1][:, 1])
	drainP = Float64.(data[1][:, 2])
    println("Read designated drain cells.")
    return drainCell, drainP
end

function WriteCells(cell::Array{Cell, 1}, fileName::String, boundFlag::Int64, drainCell::Array{Int64, 1}, drainP::Array{Float64, 1})
    # summarize cell properties to file
    csvfile = open(fileName,"w")
    line_out = "cellIndex" * "," * "x" * "," * "y" * "," * "z" * "," *
        "vol" * "," * "P" * "," * "S" * "," * "material" * "," * "drain status"
    println(csvfile, line_out)
	# drain status
	drainStatus = zeros(Int64, length(cell)).-1 	# default value is no drain condition assigned
	for (i, indx) in enumerate(drainCell)
		if cell[indx].P == drainP[i]
			drainStatus[indx] = 1
		else
			drainStatus[indx] = 0
		end
	end
    for (i, ce) in enumerate(cell)
        if (ce.vol < 1e+10) || boundFlag == 1       # if applicable, exclude boundary elements from summary
            ce.S = Sat(ce)      # update saturation estimate
            line_out = string(ce.index) * "," * string(ce.x) * "," * string(ce.y) * "," * string(ce.z) * "," *
                string(ce.vol) * "," * string(ce.P) * "," * string(ce.S) * "," * string(ce.mat.index) * "," *
				string(drainStatus[i])
			println(csvfile, line_out)
        end
    end
    close(csvfile)
    println("Wrote " * fileName * ".")
end

function WriteConnections(connect::Array{Connection, 1})
    # summarize cell properties to file
    csvfile = open("ConnectionSummary.csv","w")
    line_out = "connection" * "," * "cell1" * "," * "cell2" * "," * "del1" * "," *
        "del2" * "," * "area" * "," * "zConnect"
    println(csvfile, line_out)
    for (i, cn) in enumerate(connect)
        line_out = string(i) * "," * string(cn.cell1) * "," * string(cn.cell2) * "," * string(cn.del1) * "," *
            string(cn.del2) * "," * string(cn.A) * "," * string(cn.zConnect)
        println(csvfile, line_out)
    end
    close(csvfile)
    println("Wrote connections file.")
end

function WriteSources(stressPeriod::Array{Float64, 1}, sourceCells::Array{Int64, 1}, source::Array{Float64, 2})
    # summarize source term(s) to file
    csvfile = open("SourceSummary.csv","w")
    line_out = "cell"
    for t in stressPeriod
        line_out *= "," * string(t)
    end
    println(csvfile, line_out)
    for cellIndex in sourceCells
        line_out = string(cellIndex)
        iRow = findall(sourceCells.==cellIndex)[1]
        for iCol = 1:length(stressPeriod)
            line_out *= "," * string(source[iRow, iCol])
        end
        println(csvfile, line_out)
    end
    close(csvfile)
    println("Wrote source terms file.")
end

function WriteTimeSeries(monitor::Array{Int64, 1}, timeMon::Array{Float64, 1}, timeSeries)
    # write monitoring point time series to file
    csvfile = open("TimeSeries.csv","w")
    line_out = "time"
    for m in monitor
        line_out *= "," * "cell" * "_" * string(m)
    end
    println(csvfile, line_out)
    for (i, t) in enumerate(timeMon)
        line_out = string(t)
        for j = 1:length(monitor)
            line_out *= "," * string(timeSeries[j][i])
        end
        println(csvfile, line_out)
    end
    close(csvfile)
    println("Wrote time series file.")
end

function WriteVectors(cell::Array{Cell, 1}, vx::Array{Float64, 1}, vy::Array{Float64, 1}, vz::Array{Float64, 1})
    # write velocity vectors at cell node points
    csvfile = open("Velocity.csv","w")
    line_out = "x" * "," * "y" * "," * "z" * "," *
        "vx" * "," * "vy" * "," * "vz"
    println(csvfile, line_out)
    for (i, ce) in enumerate(cell)
        if (ce.vol < 1e+10)         # exclude boundary elements
            line_out = string(ce.x) * "," * string(ce.y) * "," * string(ce.z) * "," *
                string(vx[i]) * "," * string(vy[i]) * "," * string(vz[i])
            println(csvfile, line_out)
        end
    end
    close(csvfile)
    println("Wrote velocity vectors file.")
end


text:
## This file contains functions and types used for GBS data geneartion
"function yp sample QTL positions (randomly across the genome)"
function sampleQTLPosition(totalQTL::Int64)
    println(
        "INFO: A total of $totalQTL QTLs are sampled randomly across $numChr chromosome(s).",
    )
    qtlChr = sample([1:numChr...], totalQTL)
    global numQTL = [sum(qtlChr .== c) for c = 1:numChr]
    qtlPos = [sort(Int64.(sample([1:(chrLen[c]*1e6)...], numQTL[c]))) for c = 1:numChr]
    qtlPos
end

"function to sample SNP positions (either randomly or via a two-stage appraoch)"
function sampleSNPPosition(totalSNP::Int64, winSize::Int64, mu::Float64, sigmasq::Float64)
    if totalSNP != 0 ## set totalSNP = 0 to trigger option 1
        println(
            "INFO: A total of $totalSNP SNPs are sampled randomly across $numChr chromosome(s).",
        )
        snpChr = sample([1:numChr...], totalSNP) # same method as sampling QTL positions
        global numSNP = [sum(snpChr .== c) for c = 1:numChr]
        ## option 1: sample SNP positions randomly across the genome
        snpPos = [sort(Int64.(sample([1:(chrLen[c]*1e6)...], numSNP[c]))) for c = 1:numChr]
    else
        ## option 2: sample SNP positions using densities
        global numSNP = Array{Int64}(undef, 0)
        snpPos = Array{Int64}(undef, 0)
        for c = 1:numChr
            numWin = Int(round(chrLen[c] * 1e6 / winSize)) # number of windows to be split in chromosome c (fixed window size, variable window number)
            # theta = sigmasq/mu # calcualte Gamma paramter: scale
            # alpha = mu^2/sigmasq # calcualte Gamma paramter: shape
            density = exp.(rand(Laplace(mu, sigmasq), numWin)) # sample SNP density for each window
            win = fill(winSize, numWin) # generate a vector (with length = numWin) contating windown size for each window
            winPos = [0; cumsum(win)[1:end-1]] # starting postion of each window
            sam = [rand(Bernoulli(1.4 * density[i]), winSize) for i = 1:numWin] # sampling the occurence of SNP at each bp using local SNP density (sampled from Gamma model for each win)
            sam2 = map(x -> findall(x .> 0), sam) # record SNP positions within each window
            sam3 = reduce(vcat, [sam2[i] .+ winPos[i] for i = 1:numWin]) # calcualte SNP position within chromosome
            numSampled = length(sam3) # number of sampled SNP postions
            numSNP = [numSNP; numSampled] # record the number of sampled SNP postions for each chromosome
            snpPos = [snpPos; [sort(Int64.(sam3[randperm(numSampled)[1:numSampled]]))]] # sort SNP positions
            println(
                "CHROMOSOME $c: $numSampled SNPs sampled with average SNP density = $(round(mean(density),digits=3)) (window size = $winSize bp).",
            )
        end
    end
    snpPos
end

"function to sample variants (incl. SNP and QTL) allele frequency"
function sampleAlleleFrequency(numLoci::Array{Int64}, mu::Float64, sigmasq::Float64)
    beta = ((mu - 1) * (mu^2 - mu + sigmasq)) / (sigmasq) # (mu * (1 - mu) / sigmasq -1) * (1 - mu)
    alpha = (-beta * mu) / (mu - 1) # (mu * (1 - mu) / sigmasq -1) * mu # (-beta * mu) / (mu - 1)
    af = Array{Float64}(undef, 0)
    for c = 1:size(numLoci, 1)
        tmp = rand(Beta(alpha, beta), numLoci[c] * 3)
        af = [af; [tmp[(tmp.>0.05).&(tmp.<0.95)][1:numLoci[c]]]] # sample allele frequency from a Beta distribution
    end
    af
end

"function to generate foudner SNPs"
function makeFounderSNPs(founders, snpAF::Array{Any,1})
    for c = 1:size(snpAF, 1)
        numHaps = 2 * size(founders, 1)
        haps =
            [rand(Bernoulli(snpAF[c][i]), 1)[1] for i = 1:size(snpAF[c], 1), j = 1:numHaps]
        for i = 1:size(founders, 1)
            founders[i].MatChrs[c].SNPs = haps[:, (2*i)-1]
            founders[i].PatChrs[c].SNPs = haps[:, 2*i]
        end
    end
    founders
end

"function to generate foudner QTLs"
function makeFounderQTL(founders, qtlAF::Array{Any,1})
    for c = 1:size(qtlAF, 1)
        numHaps = 2 * size(founders, 1)
        haps =
            [rand(Bernoulli(qtlAF[c][i]), 1)[1] for i = 1:size(qtlAF[c], 1), j = 1:numHaps]
        for i = 1:size(founders, 1)
            founders[i].MatChrs[c].QTL = haps[:, (2*i)-1]
            founders[i].PatChrs[c].QTL = haps[:, 2*i]
        end
    end
    founders
end

"function to fill haplotypes at individual level"
function fillHaplotypes(
    samples,
    founders,
    numChr::Int64,
    useChr::Array{Int64},
    snpPos,
    qtlPos::Array{Array{Int64,1},1},
)
    if size(founders[1].MatChrs[1].QTL, 1) == 0
        error("Haplotypes not made - need to run makeFounderSNPs(founders,snpAF)")
    else
        for c = 1:numChr
            println("CHROMOSOME $(useChr[c]): Filling haplotypes!")
            founderHaps = Array{Int64}(undef, size(snpPos[c], 1), 2 * size(founders, 1))
            for i = 1:size(founders, 1)
                founderHaps[:, 2*i-1] = founders[i].MatChrs[c].SNPs
                founderHaps[:, 2*i] = founders[i].PatChrs[c].SNPs
            end
            for an = 1:size(samples, 1)
                matOrig = [
                    samples[an].MatChrs[c].Origin[maximum(
                        findall(samples[an].MatChrs[c].Position .< snpPos[c][i]),
                    )] for i = 1:size(snpPos[c], 1)
                ]
                patOrig = [
                    samples[an].PatChrs[c].Origin[maximum(
                        findall(samples[an].PatChrs[c].Position .< snpPos[c][i]),
                    )] for i = 1:size(snpPos[c], 1)
                ]
                samples[an].MatChrs[c].SNPs =
                    [founderHaps[i, matOrig[i]] for i = 1:size(snpPos[c], 1)]
                samples[an].PatChrs[c].SNPs =
                    [founderHaps[i, patOrig[i]] for i = 1:size(snpPos[c], 1)]
            end
            if size(qtlPos, 1) > 0
                founderHaps = Array{Int64}(undef, size(qtlPos[c], 1), 2 * size(founders, 1))
                for i = 1:size(founders, 1)
                    founderHaps[:, 2*i-1] = founders[i].MatChrs[c].QTL
                    founderHaps[:, 2*i] = founders[i].PatChrs[c].QTL
                end
                for an = 1:size(samples, 1)
                    matOrig = [
                        samples[an].MatChrs[c].Origin[maximum(
                            findall(samples[an].MatChrs[c].Position .< qtlPos[c][i]),
                        )] for i = 1:size(qtlPos[c], 1)
                    ]
                    patOrig = [
                        samples[an].PatChrs[c].Origin[maximum(
                            findall(samples[an].PatChrs[c].Position .< qtlPos[c][i]),
                        )] for i = 1:size(qtlPos[c], 1)
                    ]
                    samples[an].MatChrs[c].QTL =
                        [founderHaps[i, matOrig[i]] for i = 1:size(qtlPos[c], 1)]
                    samples[an].PatChrs[c].QTL =
                        [founderHaps[i, patOrig[i]] for i = 1:size(qtlPos[c], 1)]
                end
            end
            println("CHROMOSOME $(useChr[c]): DONE!")
        end
    end
end

"function to extract haplotypes from each (diploid x2) individual"
function getHaplotypes(samples = ind)
    numInd = size(samples, 1)
    numChr = size(samples[1].MatChrs, 1)
    numSNP = [size(samples[1].MatChrs[i].SNPs, 1) for i = 1:numChr]
    totalSNP = cumsum(numSNP)
    starts = [2; totalSNP[1:(end-1)] .+ 2]
    ends = totalSNP .+ 1
    haplotypes = Array{Int64}(undef, 2 * size(samples, 1), sum(numSNP) + 1)
    for i = 1:numInd
        haplotypes[2*i-1:2*i, 1] .= samples[i].ID
        for c = 1:numChr
            mat = samples[i].MatChrs[c].SNPs
            pat = samples[i].PatChrs[c].SNPs
            haplotypes[collect(1:2:2*numInd)[i], starts[c]:ends[c]] = mat
            haplotypes[collect(2:2:2*numInd)[i], starts[c]:ends[c]] = pat
        end
    end
    haplotypes
end

"function to geneate SNP genotypes"
function getSNPGenotypes(samples = ind)
    numChr = size(samples[1].MatChrs, 1)
    numSNP = [size(samples[1].MatChrs[i].SNPs, 1) for i = 1:numChr]
    totalSNP = cumsum(numSNP) # count total number of SNPs across the genome
    starts = [2; totalSNP[1:(end-1)] .+ 2]
    ends = totalSNP .+ 1
    genotypes = Array{Int64}(undef, size(samples, 1), sum(numSNP) + 1) # define empty array to store SNP genotypes
    for i = 1:size(samples, 1)
        genotypes[i, 1] = samples[i].ID
        for c = 1:numChr
            genos = samples[i].MatChrs[c].SNPs + samples[i].PatChrs[c].SNPs
            genotypes[i, starts[c]:ends[c]] = genos
        end
    end
    genotypes
end

"function to geneate QTL genotypes"
function getQTLGenotypes(samples = ind)
    numChr = size(samples[1].MatChrs, 1)
    numQTL = [size(samples[1].MatChrs[c].QTL, 1) for c = 1:numChr]
    totalQTL = cumsum(numQTL)
    starts = ifelse(size(totalQTL, 1) == 1, 2, [2; totalQTL[1:(end-1)] .+ 2])
    ends = totalQTL .+ 1
    genotypes = Array{Int64}(undef, size(samples, 1), sum(numQTL) + 1)
    for i = 1:size(samples, 1)
        genotypes[i, 1] = samples[i].ID
        for c = 1:numChr
            genos = samples[i].MatChrs[c].QTL + samples[i].PatChrs[c].QTL
            genotypes[i, starts[c]:ends[c]] = genos
        end
    end
    genotypes
end

"function to replicate values over array - Milan Bouchet-Valat (https://github.com/JuliaLang/julia/issues/16443)"
function rep(x, lengths)
    if length(x) != length(lengths)
        throw(DimensionMismatch("vector lengths must match"))
    end
    res = similar(x, sum(lengths))
    i = 1
    for idx = 1:length(x)
        tmp = x[idx]
        for kdx = 1:lengths[idx]
            res[i] = tmp
            i += 1
        end
    end
    res
end

"function to sample read depth"
function sampleReadDepth(numLoci::Int64, numInd::Int64, meanDepth::Float64)
    # function sampleReadDepth(numLoci::Int64, numInd::Int64, meanDepth::Float64, sigmasqReadDepth::Float64, sigmasqSampleDepth::Float64, meanCallRate::Float64, sigmasqReadCallRate::Float64, sigmasqSampleCallRate::Float64)
    readDepth = rand(Gamma(meanDepth / 11, 11), numLoci)
    sampleDepth = reshape(rand(Gamma(meanDepth * 3.1, 1 / 3), numInd), numInd, 1)
    depthRaw = (sampleDepth * readDepth') ./ meanDepth
    depthMat = [rand(NegativeBinomial(depthRaw[i, j], 0.5)) for i = 1:numInd, j = 1:numLoci]
    # Caution: changes are needed here!
    dp = log10.(Base.dropdims(mean(depthMat, dims = 1); dims = 1))
    dp[isinf.(dp)] .= -1
    # sampleCallRate = reshape(rand(Laplace(meanCallRate, sqrt(sigmasqSampleCallRate / 2)), numInd), numInd, 1)
    # readCallRate = rand(Beta((meanCallRate * (1- meanCallRate) / sigmasqReadCallRate - 1) * meanCallRate, (meanCallRate * (1- meanCallRate) / sigmasqReadCallRate - 1) * (1 - meanCallRate)), numLoci)
    cr = Base.dropdims(mean(depthMat .!= 0, dims = 1); dims = 1)
    @. model(x, p) = p[1] / (1 + exp(-(x - p[2]) * p[3]))
    p0 = [1, 0.5, 5]
    fit_cr = curve_fit(model, dp, cr, p0)
    add_var1 = rand(Normal(0, 0.01), numLoci)
    add_var3 = rand(Beta(0.5, 1), numLoci) # rand(Uniform(-1,1), numLoci)
    for i = 1:numLoci
        par1 = fit_cr.param[1]
        par2 = fit_cr.param[2]
        par3 = fit_cr.param[3]
        output = (1.1 - add_var1[i]) / (1 + exp(-(dp[i] - par2 - add_var3[i]) * par3))
        if output > 1
            output = 1
        end
        #if output > cr[i]
        #    output = cr[i]
        #end
        zeros = findall(x -> x == 0, depthMat[:, i]) #depthMat[i, sample([1:numLoci...],Int(round(numLoci*rand(Laplace(1-callRate,0.sigmasqCallRate),1))), replace=false)] .= 0
        newzeros = Int64(floor((1 - output) * numInd))
        potentials = deleteat!([1:numInd...], zeros)
        diff = newzeros - length(zeros)
        if diff > 0
            setzero = sample(potentials, diff, replace = false)
            s = sum(depthMat[setzero, i])
            allzeros = sort([setzero; zeros...])
            k = deleteat!([1:numInd...], allzeros)
            l = sample(k, s, replace = true)
            depthMat[setzero, i] .= 0
            for j in l
                depthMat[j, i] = depthMat[j, i] + 1
            end
        end
    end
    depthMat
end

"function to generate key file for simulated GBS data"
function getKeyFile(barcodeFile, numInd, flowcell, lane, numRow = 24, numCol = 16)
    barcode = readdlm(barcodeFile, '\t', header = false)
    useBarcode = barcode[1:numInd]
    row = reduce(vcat, rep(string.(['A':'P'...]), fill(numRow, numCol)))
    col = reduce(vcat, fill([1:numRow...], numCol))
    key =
        [fill("$flowcell", numRow * numCol)[1:numInd] fill(lane, numRow * numCol)[1:numInd] useBarcode [
            "Ind_" * string(i) for i = 1:(numInd)
        ] fill("Plate1", numInd) row[1:numInd] col[1:numInd]]
    writedlm("keyFile_$(flowcell)_$(lane).txt", key)
    useBarcode
end

"""
	GBS(totalQTL, totalSNP, muDensity, sigmasqDensity, winSize, muAlleleFreq, sigmasqAlleleFreq, re, meanDepth, barcodeFile, useChr, plotOutput, writeOutput, onlyOutputGBS)

Simulate Genotyping-by-Sequencing (GBS) data.

This function generates GBS reads by inserting genomic variants into _in silico_ digested genomic fragments, ligates the polymorphic sequence with barcodes and replicates based on sequencing depth.

# Arguments
* `totalQTL`: total number of QTL to be simulated
* `totalSNP`: total number of SNPs to be simulated (set to "0" if sampling SNP positions based on density)
* `muDensity`: location parameter of log-Laplace distribution (for sampling SNP density)
* `sigmasqDensity`: scale parameter of log-Laplace distribution (for sampling SNP density)
* `winSize`: Size of window and bin for sampling SNP positions
* `muAlleleFreq`: mean of sampled allele frequency
* `sigmasqAlleleFreq`: variance of sampled allele frequency
* `re`: restriction enzyme(s) to be used
* `barcodeFile`: file containing GBS barcodes
* `useChr`: either the number of chromosomes or a set of chromosome(s) to be simulated
* `plotOutput`: set to true if graphical outputs are required
* `writeOutput`: set to true if text outputs are required
* `onlyOutputGBS`: set to true if only GBS data is kept

# Examples
```julia
julia> GBS(100, 0, -2.0, 0.001, 1000000, 0.5, 0.001, [SimGBS.ApeKI], 20.0, "GBS_Barcodes.txt", [1], false, true, true)
```
"""
function GBS(
    totalQTL::Int64,
    totalSNP::Int64,
    muSNPdensity::Float64,
    sigmasqSNPdensity::Float64,
    winSize::Int64,
    muAlleleFreq::Float64,
    sigmasqAlleleFreq::Float64,
    re,
    meanDepth::Float64,
    barcodeFile::String,
    useChr::Array{Int64},
    plotOutput::Bool,
    writeOutput::Bool,
    onlyOutputGBS::Bool,
)
    # function GBS(totalQTL::Int64, totalSNP::Int64, muSNPdensity::Float64, sigmasqSNPdensity::Float64, winSize::Int64, muAlleleFreq::Float64,sigmasqAlleleFreq::Float64, re, meanDepth::Float64, sigmasqReadDepth::Float64, sigmasqSampleDepth::Float64, meanCallRate::Float64, sigmasqReadCallRate::Float64, sigmasqSampleCallRate::Float64, barcodeFile::String, plotOutput::Bool, writeOutput::Bool,outputOnlyGBS::Bool)
    ## 1. sample variants
    ### 1.1 QTL positions, number and allele frequencies
    qtlPos = sampleQTLPosition(totalQTL)
    numQTL = length.(qtlPos)
    qtlAF = sampleAlleleFrequency(numQTL, muAlleleFreq, sigmasqAlleleFreq)
    ###  1.2 SNP positions
    snpPos = sampleSNPPosition(totalSNP, winSize, muSNPdensity, sigmasqSNPdensity)
    numSNP = length.(snpPos)
    numFrag = [length(findall(x -> x == useChr[c], GBSFrag.chr)) for c = 1:numChr] # number of GBS fragments found within each chromosome
    ### 1.3 SNP captured by GBS fragments
    snpGBS = Array{Int64}(undef, 0) # empty array to store only SNP found on GBS fragments
    hapSize = Array{Int64}(undef, 0) # empty array to store haplotype size (number of SNPs within each hap)
    fragSNP = Array{Int64}(undef, 0) # empty array to store GBS fragments (conataining SNPs) index
    for c = 1:numChr
        #### 1.3.1 extract GBS fragments from chromosome c
        chr = findall(x -> x == useChr[c], GBSFrag.chr) # index fragments from chromosome c
        numFragChr = numFrag[c] # total number of selected GBS fragments
        starts = GBSFrag.pos[chr] # starting position of selected GBS fragments
        len = GBSFrag.len[chr] # length of selected GBS fragments
        ends = starts + len .- 1 # ending position of selected GBS fragments
        #### 1.3.2 select only SNPs captured by GBS fragments and non-empty GBS fragments (contain SNPs)
        snpSelected = [
            findall(x -> x == 1, (snpPos[c] .< ends[t]) .& (snpPos[c] .> starts[t])) for
            t = 1:numFragChr
        ]
        fragSelected = findall(x -> x != [], snpSelected) # findall(x -> x != [], [snpPos[c][snpSelected[t]] for t = 1:numFragChr])
        hapSizeChr = length.(snpSelected)[fragSelected] # size of shortHaps (i.e., number of SNPs within each)
        numSNPChr = sum(hapSizeChr) # total number of SNPs found on GBS fragments
        numHapLoci = length(fragSelected) # total number of GBS fragments contating SNPs (or equivelently, number of shortHaps loci)
        #### 1.3.3 store SNP index, shortHap size and fragment index
        snpGBS = [snpGBS; [reduce(vcat, snpSelected[fragSelected])]] # store index of SNPs found on GBS fragments
        # [snpGBS; [(reduce(vcat,snpSelected[fragSelected]) .+ [0;length.(snpPos)][c])]] # store index of SNPs found on GBS fragments
        hapSize = [hapSize; hapSizeChr] # store shortHap size (number of SNPs with each shortHap)
        fragSNP = [fragSNP; fragSelected .+ [0; numFrag][c]] # store only fragments contain SNPs
        println(
            "CHROMOSOME $(useChr[c]): Found $numSNPChr SNPs on $numHapLoci GBS fragments, with an average of $(round(mean(hapSizeChr),digits=2)) SNPs per GBS fragment.",
        )
    end
    #### 1.3.4 (Optional) keep only SNPs captured by GBS fragments for the subsequent steps
    if onlyOutputGBS == true
        println("[INFO]: a total of $(sum(numSNP)) SNPs were sampled.")
        numSNP = length.(snpGBS) # update number of SNPs
        println("[INFO]: $(sum(numSNP)) SNPs captured by selected GBS SNPs.")
        snpPos = [snpPos[c][snpGBS[c]] for c = 1:numChr]
    end
    ### 1.4 SNP allele frequencies
    snpAF = sampleAlleleFrequency(numSNP, muAlleleFreq, sigmasqAlleleFreq)
    ### 1.5 founder QTL and SNPs
    makeFounderQTL(founders, qtlAF)
    makeFounderSNPs(founders, snpAF)
    ### 1.6 fill haplotypes
    @time fillHaplotypes(ind, founders, numChr, useChr, snpPos, qtlPos)
    ### 1.7 extract haplotypes
    haps = getHaplotypes(ind)
    ## 2. generate GBS reads
    ### 2.1 sample read depth
    totalInd = size(ind, 1) # number of samples to be simulated
    totalFrag = length(fragSNP) # number of loci
    totalDepth = Array{Int64}(undef, totalInd, totalFrag)
    ### 2.2 inster variants
    #### 2.2.1 subset haplotype alleel by chromosome
    hapStart = [1; numSNP[1:end-1] .+ 1] .+ 1
    hapEnd = cumsum(numSNP) .+ 1
    #### 2.2.2 genearte FASTQ file(s)
    flowcell = "ABC12AAXX" # flowcell name
    multiplex = 192 # 384 number of samples per lane
    numLane = Int(ceil(totalInd / multiplex)) # numebr of lanes in total
    for lane = 1:numLane
        if totalInd < multiplex * lane
            numInd = totalInd - multiplex * (lane - 1) # number of samples to be 'sequenced'
        else
            numInd = multiplex # number of samples to be 'sequenced'
        end
        readDepth = sampleReadDepth(totalFrag, numInd, meanDepth) # depth
        useInd = [i + (lane - 1) * multiplex for i = 1:numInd] # which individual gets sequenced (ascending order)
        totalDepth[useInd, :] = readDepth
        println(
            "[INFO] On Lane $lane of Flowcell $flowcell: Average GBS read depth equals to $(round(mean(readDepth),digits=2)).",
        )
        matHapDepth = [rand(Binomial(readDepth[i, j])) for i = 1:numInd, j = 1:totalFrag]
        patHapDepth = readDepth - matHapDepth # depth of parternal haplotypes
        matDepthFwd = [rand(Binomial(matHapDepth[i, j])) for i = 1:numInd, j = 1:totalFrag]
        matDepthRevComp = matHapDepth - matDepthFwd
        patDepthFwd = [rand(Binomial(patHapDepth[i, j])) for i = 1:numInd, j = 1:totalFrag]
        patDepthRevComp = patHapDepth - patDepthFwd
        println(
            "[INFO] On Lane $lane of Flowcell $flowcell: Generating GBS data for $numInd samples.",
        )
        barcodes = getKeyFile(barcodeFile, numInd, flowcell, lane) # extract barcodes
        io = GZip.open("$(flowcell)_$(lane)_fastq.txt.gz", "w") # output file
        numReadsTotal = 0 # reads counter
        for c = 1:numChr
            ##### 2.2.2.1 extract GBS fragments and SNPs
            chr = findall(x -> x == useChr[c], GBSFrag.chr) # index fragments from chromosome c
            numFragChr = numFrag[c] # total number of selected GBS fragments
            starts = GBSFrag.pos[chr] # starting position of selected GBS fragments
            len = GBSFrag.len[chr] # length of selected GBS fragments
            ends = starts + len .- 1 # ending position of selected GBS fragments
            frags = GBSFrag.frag[chr] # select only fragments from chromosome c
            ##### 2.2.2.2 select only SNPs that captured by each GBS fragment and non-empty (containing SNP(s)) GBS fragments
            snpSelected = [
                findall(x -> x == 1, (snpPos[c] .< ends[t]) .& (snpPos[c] .> starts[t])) for t = 1:numFragChr
            ]
            fragSelected = findall(x -> x != [], snpSelected) # findall(x -> x != [], [snpPos[c][snpSelected[t]] for t = 1:numFragChr])
            hapSizeChr = length.(snpSelected[fragSelected]) # number of SNPs within each non-empty GBS fragment
            numSNPChr = sum(hapSizeChr) # total number of SNPs found on GBS fragments
            numHapLoci = length(fragSelected) # total number of GBS fragments contating SNPs (or equivelently, number of shortHaps loci)
            ##### 2.2.2.3 genearte GBS reads with variants
            hapChr = haps[:, hapStart[c]:hapEnd[c]]
            for f = 1:length(fragSelected)
                k = fragSelected[f]
                frag = frags[k] # extract kth fragment in chromsome c
                sites = snpPos[c][snpSelected[k]] .- starts[k] # (starts[k].-1) # SNP sites on this GBS fragment
                startPos = [1; sites .+ 1] # starting position of segments of the specified GBS fragment
                endPos = [sites .- 1; len[k]] # ending position of segments of the specified GBS fragment
                seg = [SubString(frag, startPos[j], endPos[j]) for j = 1:length(startPos)] # poistional info of small segments of the specified GBS fragment
                bases = [split(map(x -> nucl[x], frag[sites]), "") split(
                    map(x -> comp[nucl[x]], frag[sites]),
                    "",
                )] # define SNPs (ref. allele = ref. allele found on the genome; alt. allele = complement of the ref. allele. E.g., A/T,C/G only)
                hap = hapChr[:, snpSelected[k]] .+ 1 # extract shortHaps (i.e. short haplotypes defined by GBS fragmentation), adding 1 here to covert
                matHap = hap[1:2:(totalInd*2), :]
                patHap = hap[2:2:(totalInd*2), :]
                if length(frag) < 101
                    stops = length(frag) # read the entire fragment if its less than 101 bp
                else
                    stops = 101 # if the length of fragment is longer than 101, stop at 101
                end
                matHapChrTemp = [
                    join(
                        [
                            [seg[j] * bases[j, matHap[i, j]] for j = 1:length(sites)]
                            seg[length(sites)+1]
                            re[1].overhang[1]
                        ],
                    ) for i in useInd
                ]
                patHapChrTemp = [
                    join(
                        [
                            [seg[j] * bases[j, patHap[i, j]] for j = 1:length(sites)]
                            seg[length(sites)+1]
                            re[1].overhang[1]
                        ],
                    ) for i in useInd
                ]
                matHapChr = [matHapChrTemp[i][1:stops] for i = 1:numInd]
                patHapChr = [patHapChrTemp[i][1:stops] for i = 1:numInd]
                matHapRevCompChr = [
                    reverse(map(x -> comp[nucl[x]], matHapChrTemp[i]))[1:stops] for
                    i = 1:numInd
                ] # reverse complement GBS reads of maternal haplotypes
                patHapRevCompChr = [
                    reverse(map(x -> comp[nucl[x]], patHapChrTemp[i]))[1:stops] for
                    i = 1:numInd
                ] # reverse complement GBS reads of paternal haplotypes
                readGBS = collect(
                    Iterators.flatten([
                        [
                            rep([barcodes[i] .* matHapChrTemp[i]], matDepthFwd[i, f])
                            rep(
                                [barcodes[i] .* matHapRevCompChr[i]],
                                matDepthRevComp[i, f],
                            )
                            rep([barcodes[i] .* patHapChr[i]], patDepthFwd[i, f])
                            rep(
                                [barcodes[i] .* patHapRevCompChr[i]],
                                patDepthRevComp[i, f],
                            )
                        ] for i = 1:numInd
                    ]),
                )
                totalReads = shuffle(readGBS)
                numReads = length(totalReads)
                numReadsTotal = numReadsTotal + numReads
                writedlm(
                    io,
                    [
                        "@SIM001:001:ABC12AAXX:$lane:0000:0000:$r 1:N:0:0\n" *
                        totalReads[r] *
                        "\n+\n" *
                        repeat("I", length(totalReads[r])) for r = 1:length(totalReads)
                    ],
                    quotes = false,
                )
            end
        end
        println("INFO: A total of $numReadsTotal GBS reads generated.")
        close(io)
    end
    ## 3. GBS data
    #### 3.1.1 SNP and QTL genotypes
    snpGeno = getSNPGenotypes(ind)
    qtlGeno = getQTLGenotypes(ind)
    ### 3.2 SNP and QTL datasets [ID chromosome position allele_frequency]
    qtlData =
        [["QTL_" * string(i) for i = 1:totalQTL] reduce(vcat, [repeat([useChr[i]], inner = numQTL[i]) for i = 1:numChr]) reduce(vcat, qtlPos) reduce(vcat, qtlAF)]
    snpData =
        [["SNP_" * string(i) for i = 1:sum(numSNP)] reduce(vcat, [repeat([useChr[i]], inner = numSNP[i]) for i = 1:numChr]) reduce(vcat, snpPos) reduce(vcat, snpAF)]
    ### 3.3 GBS fragments conatin SNPs
    snpFragGBS = [GBSFrag.id[collect(Iterators.flatten(fragSNP))] GBSFrag.chr[collect(
        Iterators.flatten(fragSNP),
    )] GBSFrag.pos[collect(Iterators.flatten(fragSNP))] GBSFrag.len[collect(
        Iterators.flatten(fragSNP),
    )] GBSFrag.frag[collect(Iterators.flatten(fragSNP))]]
    ### 3.4 , shortHaps
    hapIndex = [0; rep(fragSNP, hapSize)] # indexing shortHaps
    hapGBS = vcat(hapIndex', haps)
    ### 3.5 SNP depth, call rates
    snpDepth = transpose(
        reshape(
            vcat(
                [
                    rep(totalDepth[i, :], collect(Iterators.flatten(hapSize))) for
                    i = 1:totalInd
                ]...,
            ),
            sum(hapSize),
            totalInd,
        ),
    )
    if writeOutput == true
        writedlm("snpGeno.txt", snpGeno)
        writedlm("qtlGeno.txt", qtlGeno)
        writedlm("snpInfo.txt", snpData)
        writedlm("qtlInfo.txt", qtlData)
        writedlm("snpFragGBS.txt", snpFragGBS)
        writedlm("shortHap.txt", hapGBS)
        writedlm("readDepth.txt", totalDepth)
        writedlm("snpDepth.txt", snpDepth)
        # writedlm("allele.txt",allele[:,2:end])
    end
    if plotOutput == true
        histogram(
            snpAF,
            normalize = :probability,
            title = s"SNP Allele Frequency",
            bins = 1000,
            label = "",
            fmt = :png,
        )
        savefig("snpAF")
        histogram(
            qtlAF,
            normalize = :probability,
            title = "QTL Allele Frequency",
            bins = 1000,
            label = "",
            fmt = :png,
        )
        savefig("qtlAF")
        histogram(
            hapSize,
            normalize = :probability,
            title = "Number of SNP per GBS Fragment",
            bins = 1000,
            label = "",
            fmt = :png,
        )
        savefig("snpPerTag")
    end
end


text:
using Test;

include(joinpath("IGA","nurbs_toolbox","runtests.jl"));


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

"""
    GridTopology(dims)

A data structure for grid topologies with `dims` elements.
"""
struct GridTopology{D} <: Topology
  dims::Dims{D}
end

GridTopology(dims::Vararg{Int,D}) where {D} = GridTopology{D}(dims)

paramdim(::GridTopology{D}) where {D} = D

Base.size(t::GridTopology) = t.dims

"""
    elem2cart(t, e)

Return the Cartesian indices of the `e`-th element of
the grid topology `t`.
"""
elem2cart(t::GridTopology, e) = CartesianIndices(t.dims)[e].I

"""
    cart2elem(t, i, j, k, ...)

Return the linear index of the element of the
grid topology `t` with Cartesian indices
`i`, `j`, `k`, ...
"""
cart2elem(t::GridTopology, ijk...) = LinearIndices(t.dims)[ijk...]

"""
    corner2cart(t, v)

Return the Cartesian indices of the element of the
grid topology `t` with top left vertex `v`.
"""
corner2cart(t, v) = CartesianIndices(t.dims .+ 1)[v].I

"""
    cart2corner(t, i, j, k, ...)

Return the linear index of the top left vertex of the
element of the grid topology `t` with Cartesian indices
`i`, `j`, `k`, ...
"""
cart2corner(t::GridTopology, ijk...) = LinearIndices(t.dims .+ 1)[ijk...]

"""
    elem2corner(t, e)

Return the linear index of the top left vertex of the
`e`-th element of the grid topology `t`.
"""
elem2corner(t::GridTopology, e) = cart2corner(t, elem2cart(t, e)...)

"""
    corner2elem(t, v)

Return the linear index of the element of the grid
topology `t` with top left vertex `v`
"""
corner2elem(t::GridTopology, v) = cart2elem(t, corner2cart(t, v)...)

"""
    rank2type(t, rank)

Return the polytope type of given `rank` for the grid
topology `t`.
"""
function rank2type(::GridTopology{D}, rank::Integer) where {D}
  @assert rank ≤ D "invalid rank for grid topology"
  rank == 1 && return Segment
  rank == 2 && return Quadrangle
  rank == 3 && return Hexahedron
  throw(ErrorException("not implemented"))
end

# ---------------------
# HIGH-LEVEL INTERFACE
# ---------------------

nvertices(t::GridTopology) = prod(t.dims .+ 1)

function faces(t::GridTopology{D}, rank) where {D}
  if rank == D
    elements(t)
  elseif rank == D - 1
    facets(t)
  else
    throw(ErrorException("not implemented"))
  end
end

function element(t::GridTopology{D}, ind) where {D}
  ∂ = Boundary{D,0}(t)
  T = rank2type(t, D)
  connect(Tuple(∂(ind)), T)
end

nelements(t::GridTopology) = prod(t.dims)

function facet(t::GridTopology{D}, ind) where {D}
  ∂ = Boundary{D-1,0}(t)
  T = rank2type(t, D-1)
  connect(Tuple(∂(ind)), T)
end

nfacets(t::GridTopology{1}) =
  t.dims[1] + 1

nfacets(t::GridTopology{2}) =
  2prod(t.dims) +
  t.dims[1] +
  t.dims[2]

nfacets(t::GridTopology{3}) =
  3prod(t.dims) +
  prod(t.dims[[1,2]]) +
  prod(t.dims[[1,3]]) +
  prod(t.dims[[2,3]])