text:
__precompile__()

module BinDeps

using Compat
using Compat.Libdl
using Compat.Pkg

export @build_steps, find_library, download_cmd, unpack_cmd,
    Choice, Choices, CCompile, FileDownloader, FileRule,
    ChangeDirectory, FileUnpacker, prepare_src,
    autotools_install, CreateDirectory, MakeTargets,
    MAKE_CMD, glibc_version

function find_library(pkg,libname,files)
    Base.warn_once("BinDeps.find_library is deprecated; use Base.find_library instead.")
    dl = C_NULL
    for filename in files
        dl = Libdl.dlopen_e(joinpath(Pkg.dir(),pkg,"deps","usr","lib",filename))
        if dl != C_NULL
            ccall(:add_library_mapping,Cint,(Ptr{Cchar},Ptr{Cvoid}),libname,dl)
            return true
        end

        dl = Libdl.dlopen_e(filename)
        if dl != C_NULL
            ccall(:add_library_mapping,Cint,(Ptr{Cchar},Ptr{Cvoid}),libname,dl)
            return true
        end
    end

    dl = Libdl.dlopen_e(libname)
    dl != C_NULL ? true : false
end

macro make_rule(condition,command)
    quote
        if(!$(esc(condition)))
            $(esc(command))
            @assert $(esc(condition))
        end
    end
end

abstract type BuildStep end

downloadcmd = nothing
function download_cmd(url::AbstractString, filename::AbstractString)
    global downloadcmd
    if downloadcmd === nothing
        for download_engine in (Compat.Sys.iswindows() ? ("C:\\Windows\\System32\\WindowsPowerShell\\v1.0\\powershell",
                :powershell, :curl, :wget, :fetch) : (:curl, :wget, :fetch))
            if endswith(string(download_engine), "powershell")
                checkcmd = `$download_engine -NoProfile -Command ""`
            else
                checkcmd = `$download_engine --help`
            end
            try
                if success(checkcmd)
                    downloadcmd = download_engine
                    break
                end
            catch
                continue # don't bail if one of these fails
            end
        end
    end
    if downloadcmd == :wget
        return `$downloadcmd -O $filename $url`
    elseif downloadcmd == :curl
        return `$downloadcmd -f -o $filename -L $url`
    elseif downloadcmd == :fetch
        return `$downloadcmd -f $filename $url`
    elseif endswith(string(downloadcmd), "powershell")
        tls_cmd = "[Net.ServicePointManager]::SecurityProtocol = [Net.SecurityProtocolType]::Tls12"
        download_cmd = "(new-object net.webclient).DownloadFile(\"$(url)\", \"$(filename)\")"
        return `$downloadcmd -NoProfile -Command "$(tls_cmd); $(download_cmd)"`
    else
        extraerr = Compat.Sys.iswindows() ? "check if powershell is on your path or " : ""
        error("No download agent available; $(extraerr)install curl, wget, or fetch.")
    end
end

if Compat.Sys.isunix() && Sys.KERNEL != :FreeBSD
    function unpack_cmd(file,directory,extension,secondary_extension)
        if ((extension == ".gz" || extension == ".Z") && secondary_extension == ".tar") || extension == ".tgz"
            return (`tar xzf $file --directory=$directory`)
        elseif (extension == ".bz2" && secondary_extension == ".tar") || extension == ".tbz"
            return (`tar xjf $file --directory=$directory`)
        elseif extension == ".xz" && secondary_extension == ".tar"
            return pipeline(`unxz -c $file `, `tar xv --directory=$directory`)
        elseif extension == ".tar"
            return (`tar xf $file --directory=$directory`)
        elseif extension == ".zip"
            return (`unzip -x $file -d $directory`)
        elseif extension == ".gz"
            return pipeline(`mkdir $directory`, `cp $file $directory`, `gzip -d $directory/$file`)
        end
        error("I don't know how to unpack $file")
    end
end

if Sys.KERNEL == :FreeBSD
    # The `tar` on FreeBSD can auto-detect the archive format via libarchive.
    # The supported formats can be found in libarchive-formats(5).
    # For NetBSD and OpenBSD, libarchive is not available.
    # For macOS, it is. But the previous unpack function works fine already.
    function unpack_cmd(file, dir, ext, secondary_ext)
        tar_args = ["--no-same-owner", "--no-same-permissions"]
        return pipeline(
            `mkdir -p $dir`,
            `tar -xf $file -C $dir $tar_args`)
    end
end

if Compat.Sys.iswindows()
    const exe7z = joinpath(Compat.Sys.BINDIR, "7z.exe")

    function unpack_cmd(file,directory,extension,secondary_extension)
        if ((extension == ".Z" || extension == ".gz" || extension == ".xz" || extension == ".bz2") &&
                secondary_extension == ".tar") || extension == ".tgz" || extension == ".tbz"
            return pipeline(`$exe7z x $file -y -so`, `$exe7z x -si -y -ttar -o$directory`)
        elseif (extension == ".zip" || extension == ".7z" || extension == ".tar" ||
                (extension == ".exe" && secondary_extension == ".7z"))
            return (`$exe7z x $file -y -o$directory`)
        end
        error("I don't know how to unpack $file")
    end
end

mutable struct SynchronousStepCollection
    steps::Vector{Any}
    cwd::AbstractString
    oldcwd::AbstractString
    SynchronousStepCollection(cwd) = new(Any[],cwd,cwd)
    SynchronousStepCollection() = new(Any[],"","")
end

import Base: push!, run, |
push!(a::SynchronousStepCollection,args...) = push!(a.steps,args...)

mutable struct ChangeDirectory <: BuildStep
    dir::AbstractString
end

mutable struct CreateDirectory <: BuildStep
    dest::AbstractString
    mayexist::Bool
    CreateDirectory(dest, me) = new(dest,me)
    CreateDirectory(dest) = new(dest,true)
end

struct RemoveDirectory <: BuildStep
    dest::AbstractString
end

mutable struct FileDownloader <: BuildStep
    src::AbstractString     # url
    dest::AbstractString    # local_file
end

mutable struct ChecksumValidator <: BuildStep
    sha::AbstractString
    path::AbstractString
end

mutable struct FileUnpacker <: BuildStep
    src::AbstractString     # archive file
    dest::AbstractString    # directory to unpack into
    target::AbstractString  # file or directory inside the archive to test
                            # for existence (or blank to check for a.tgz => a/)
end

mutable struct MakeTargets <: BuildStep
    dir::AbstractString
    targets::Vector{String}
    env::Dict
    MakeTargets(dir,target;env = Dict{AbstractString,AbstractString}()) = new(dir,target,env)
    MakeTargets(target::Vector{<:AbstractString};env = Dict{AbstractString,AbstractString}()) = new("",target,env)
    MakeTargets(target::String;env = Dict{AbstractString,AbstractString}()) = new("",[target],env)
    MakeTargets(;env = Dict{AbstractString,AbstractString}()) = new("",String[],env)
end

mutable struct AutotoolsDependency <: BuildStep
    src::AbstractString     #src direcory
    prefix::AbstractString
    builddir::AbstractString
    configure_options::Vector{AbstractString}
    libtarget::Vector{AbstractString}
    include_dirs::Vector{AbstractString}
    lib_dirs::Vector{AbstractString}
    rpath_dirs::Vector{AbstractString}
    installed_libpath::Vector{String} # The library is considered installed if any of these paths exist
    config_status_dir::AbstractString
    force_rebuild::Bool
    env
    AutotoolsDependency(;srcdir::AbstractString = "", prefix = "", builddir = "", configure_options=AbstractString[], libtarget = AbstractString[], include_dirs=AbstractString[], lib_dirs=AbstractString[], rpath_dirs=AbstractString[], installed_libpath = String[], force_rebuild=false, config_status_dir = "", env = Dict{String,String}()) =
        new(srcdir,prefix,builddir,configure_options,isa(libtarget,Vector) ? libtarget : AbstractString[libtarget],include_dirs,lib_dirs,rpath_dirs,installed_libpath,config_status_dir,force_rebuild,env)
end

### Choices

mutable struct Choice
    name::Symbol
    description::AbstractString
    step::SynchronousStepCollection
    Choice(name,description,step) = (s=SynchronousStepCollection();lower(step,s);new(name,description,s))
end

mutable struct Choices <: BuildStep
    choices::Vector{Choice}
    Choices() = new(Choice[])
    Choices(choices::Vector{Choice}) = new(choices)
end

push!(c::Choices, args...) = push!(c.choices, args...)

function run(c::Choices)
    println()
    Compat.@info("There are multiple options available for installing this dependency:")
    while true
        for x in c.choices
            println("- "*string(x.name)*": "*x.description)
        end
        while true
            print("Plese select the desired method: ")
            method = Symbol(chomp(readline(STDIN)))
            for x in c.choices
                if(method == x.name)
                    return run(x.step)
                end
            end
            Compat.@warn("Invalid method")
        end
    end
end

mutable struct CCompile <: BuildStep
    srcFile::AbstractString
    destFile::AbstractString
    options::Vector{String}
    libs::Vector{String}
end

lower(cc::CCompile,c) = lower(FileRule(cc.destFile,`gcc $(cc.options) $(cc.srcFile) $(cc.libs) -o $(cc.destFile)`),c)
##

mutable struct DirectoryRule <: BuildStep
    dir::AbstractString
    step
end

mutable struct PathRule <: BuildStep
    path::AbstractString
    step
end

function meta_lower(a::Expr,blk::Expr,collection)
    if(a.head == :block || a.head == :tuple)
        for x in a.args
            if(isa(x,Expr))
                if(x.head == :block)
                    new_collection = gensym()
                    push!(blk.args,quote
                        $(esc(new_collection)) = SynchronousStepCollection($(esc(collection)).cwd)
                        push!($(esc(collection)),$(esc(new_collection)))
                    end)
                    meta_lower(x,blk,new_collection)
                 elseif(x.head != :line)
                     push!(blk.args,quote
                         lower($(esc(x)), $(esc(collection)))
                     end)
                 end
            elseif(!isa(x,LineNumberNode))
                meta_lower(x,blk,collection)
            end
        end
    else
        push!(blk.args,quote
            $(esc(collection)),lower($(esc(a)), $(esc(collection)))
        end)
    end
end

function meta_lower(a::Tuple,blk::Expr,collection)
    for x in a
        meta_lower(a,blk,collection)
    end
end

function meta_lower(a,blk::Expr,collection)
    push!(blk.args,quote
        $(esc(collection)), lower($(esc(a)), $(esc(collection)))
    end)
end

macro dependent_steps(steps)
    blk = Expr(:block)
    meta_lower(steps,blk,:collection)
    blk
end

macro build_steps(steps)
    collection = gensym()
    blk = Expr(:block)
    push!(blk.args,quote
        $(esc(collection)) = SynchronousStepCollection()
    end)
    meta_lower(steps,blk,collection)
    push!(blk.args, quote; $(esc(collection)); end)
    blk
end

src(b::BuildStep) = b.src
dest(b::BuildStep) = b.dest

(|)(a::BuildStep,b::BuildStep) = SynchronousStepCollection()
function (|)(a::SynchronousStepCollection,b::SynchronousStepCollection)
    if a.cwd == b.cwd
        append!(a.steps,b.steps)
    else
        push!(a.steps,b)
    end
    a
end
(|)(a::SynchronousStepCollection,b::Function) = (lower(b,a);a)
(|)(a::SynchronousStepCollection,b) = (lower(b,a);a)

(|)(b::Function,a::SynchronousStepCollection) = (c=SynchronousStepCollection(); ((c|b)|a))
(|)(b,a::SynchronousStepCollection) = (c=SynchronousStepCollection(); ((c|b)|a))

# Create any of these files
mutable struct FileRule <: BuildStep
    file::Array{AbstractString}
    step
    FileRule(file::AbstractString,step) = FileRule(AbstractString[file],step)
    function FileRule(files::Vector{AbstractString},step)
        new(files,@build_steps (step,) )
    end
end
FileRule(files::Vector{T},step) where {T <: AbstractString} = FileRule(AbstractString[f for f in files],step)

function lower(s::ChangeDirectory,collection)
    if !isempty(collection.steps)
        error("Change of directory must be the first instruction")
    end
    collection.cwd = s.dir
end
lower(s::Nothing,collection) = nothing
lower(s::Function,collection) = push!(collection,s)
lower(s::CreateDirectory,collection) = @dependent_steps ( DirectoryRule(s.dest,()->(mkpath(s.dest))), )
lower(s::RemoveDirectory,collection) = @dependent_steps ( `rm -rf $(s.dest)` )
lower(s::BuildStep,collection) = push!(collection,s)
lower(s::Base.AbstractCmd,collection) = push!(collection,s)
lower(s::FileDownloader,collection) = @dependent_steps (
    CreateDirectory(dirname(s.dest), true),
    ()->Compat.@info("Downloading file $(s.src)"),
    FileRule(s.dest, download_cmd(s.src, s.dest)),
    ()->Compat.@info("Done downloading file $(s.src)")
)
lower(s::ChecksumValidator,collection) = isempty(s.sha) || @dependent_steps ()->sha_check(s.path, s.sha)
function splittarpath(path)
    path,extension = splitext(path)
    base_filename,secondary_extension = splitext(path)
    if extension == ".tgz" || extension == ".tbz" || extension == ".zip" && !isempty(secondary_extension)
        base_filename *= secondary_extension
        secondary_extension = ""
    end
    (base_filename,extension,secondary_extension)
end
function lower(s::FileUnpacker,collection)
    base_filename,extension,secondary_extension = splittarpath(s.src)
    target = if isempty(s.target)
        basename(base_filename)
    elseif s.target == "."
        ""
    else
        s.target
    end
    @dependent_steps begin
        CreateDirectory(dirname(s.dest),true)
        PathRule(joinpath(s.dest,target),unpack_cmd(s.src,s.dest,extension,secondary_extension))
    end
end

adjust_env(env) = merge(ENV,env)  # s.env overrides ENV

if Compat.Sys.isunix()
    function lower(a::MakeTargets,collection)
        cmd = `make -j8`

        if Sys.KERNEL == :FreeBSD
            jobs = readchomp(`make -V MAKE_JOBS_NUMBER`)
            if isempty(jobs)
                jobs = readchomp(`sysctl -n hw.ncpu`)
            end
            # Tons of project have written their Makefile in GNU Make only syntax,
            # but the implementation of `make` on FreeBSD system base is `bmake`
            cmd = `gmake -j$jobs`
        end

        if !isempty(a.dir)
            cmd = `$cmd -C $(a.dir)`
        end
        if !isempty(a.targets)
            cmd = `$cmd $(a.targets)`
        end
        @dependent_steps ( setenv(cmd, adjust_env(a.env)), )
    end
end
Compat.Sys.iswindows() && (lower(a::MakeTargets,collection) = @dependent_steps ( setenv(`make $(!isempty(a.dir) ? "-C "*a.dir : "") $(a.targets)`, adjust_env(a.env)), ))
lower(s::SynchronousStepCollection,collection) = (collection|=s)

lower(s) = (c=SynchronousStepCollection();lower(s,c);c)

#run(s::MakeTargets) = run(@make_steps (s,))

function lower(s::AutotoolsDependency,collection)
    prefix = s.prefix
    if Compat.Sys.iswindows()
        prefix = replace(replace(s.prefix, "\\" => "/"), "C:/" => "/c/")
    end
    cmdstring = "pwd && ./configure --prefix=$(prefix) "*join(s.configure_options," ")

    env = adjust_env(s.env)

    for path in s.include_dirs
        if !haskey(env,"CPPFLAGS")
            env["CPPFLAGS"] = ""
        end
        env["CPPFLAGS"]*=" -I$path"
    end

    for path in s.lib_dirs
        if !haskey(env,"LDFLAGS")
            env["LDFLAGS"] = ""
        end
        env["LDFLAGS"]*=" -L$path"
    end

    for path in s.rpath_dirs
        if !haskey(env,"LDFLAGS")
            env["LDFLAGS"] = ""
        end
        env["LDFLAGS"]*=" -Wl,-rpath -Wl,$path"
    end

    if s.force_rebuild
        @dependent_steps begin
            RemoveDirectory(s.builddir)
        end
    end

    @static if Compat.Sys.isunix()
        @dependent_steps begin
            CreateDirectory(s.builddir)
            begin
                ChangeDirectory(s.builddir)
                FileRule(isempty(s.config_status_dir) ? "config.status" : joinpath(s.config_status_dir,"config.status"), setenv(`$(s.src)/configure $(s.configure_options) --prefix=$(prefix)`,env))
                FileRule(s.libtarget,MakeTargets(;env=s.env))
                MakeTargets("install";env=env)
            end
        end
    end

    @static if Compat.Sys.iswindows()
        @dependent_steps begin
            ChangeDirectory(s.src)
            FileRule(isempty(s.config_status_dir) ? "config.status" : joinpath(s.config_status_dir,"config.status"),setenv(`sh -c $cmdstring`,env))
            FileRule(s.libtarget,MakeTargets())
            MakeTargets("install")
        end
    end
end

function run(f::Function)
    f()
end

function run(s::FileRule)
    if !any(map(isfile,s.file))
        run(s.step)
        if !any(map(isfile,s.file))
            error("File $(s.file) was not created successfully (Tried to run $(s.step) )")
        end
    end
end
function run(s::DirectoryRule)
    Compat.@info("Attempting to create directory $(s.dir)")
    if !isdir(s.dir)
        run(s.step)
        if !isdir(s.dir)
            error("Directory $(s.dir) was not created successfully (Tried to run $(s.step) )")
        end
    else
        Compat.@info("Directory $(s.dir) already exists")
    end
end

function run(s::PathRule)
    if !ispath(s.path)
        run(s.step)
        if !ispath(s.path)
            error("Path $(s.path) was not created successfully (Tried to run $(s.step) )")
        end
    else
        Compat.@info("Path $(s.path) already exists")
    end
end

function run(s::BuildStep)
    error("Unimplemented BuildStep: $(typeof(s))")
end
function run(s::SynchronousStepCollection)
    for x in s.steps
        if !isempty(s.cwd)
            Compat.@info("Changing directory to $(s.cwd)")
            cd(s.cwd)
        end
        run(x)
        if !isempty(s.oldcwd)
            Compat.@info("Changing directory to $(s.oldcwd)")
            cd(s.oldcwd)
        end
    end
end

const MAKE_CMD = Compat.Sys.isbsd() && !Compat.Sys.isapple() ? `gmake` : `make`

function prepare_src(depsdir,url, downloaded_file, directory_name)
    local_file = joinpath(joinpath(depsdir,"downloads"),downloaded_file)
    @build_steps begin
        FileDownloader(url,local_file)
        FileUnpacker(local_file,joinpath(depsdir,"src"),directory_name)
    end
end

function autotools_install(depsdir,url, downloaded_file, configure_opts, directory_name, directory, libname, installed_libname, confstatusdir)
    prefix = joinpath(depsdir,"usr")
    libdir = joinpath(prefix,"lib")
    srcdir = joinpath(depsdir,"src",directory)
    dir = joinpath(joinpath(depsdir,"builds"),directory)
    prepare_src(depsdir,url, downloaded_file,directory_name) |
    @build_steps begin
        AutotoolsDependency(srcdir=srcdir,prefix=prefix,builddir=dir,configure_options=configure_opts,libtarget=libname,installed_libpath=[joinpath(libdir,installed_libname)],config_status_dir=confstatusdir)
    end
end
autotools_install(depsdir,url, downloaded_file, configure_opts, directory_name, directory, libname, installed_libname) = autotools_install(depsdir,url, downloaded_file, configure_opts, directory_name, directory, libname, installed_libname, "")
autotools_install(depsdir,url, downloaded_file, configure_opts, directory, libname)=autotools_install(depsdir,url,downloaded_file,configure_opts,directory,directory,libname,libname)

autotools_install(args...) = error("autotools_install has been removed")

function eval_anon_module(context, file)
    m = Module(:__anon__)
    if isdefined(Base, Symbol("@__MODULE__"))
        Core.eval(m, :(ARGS=[$context]))
        Base.include(m, file)
    else
        body = Expr(:toplevel, :(ARGS=[$context]), :(include($file)))
        Core.eval(m, body)
    end
    return
end

"""
    glibc_version()

For Linux-based systems, return the version of glibc in use. For non-glibc Linux and
other platforms, returns `nothing`.
"""
function glibc_version()
    Compat.Sys.islinux() || return
    libc = ccall(:jl_dlopen, Ptr{Cvoid}, (Ptr{Cvoid}, UInt32), C_NULL, 0)
    ptr = Libdl.dlsym_e(libc, :gnu_get_libc_version)
    ptr == C_NULL && return # non-glibc
    v = unsafe_string(ccall(ptr, Ptr{UInt8}, ()))
    occursin(Base.VERSION_REGEX, v) ? VersionNumber(v) : nothing
end

include("dependencies.jl")
include("debug.jl")
include("show.jl")

# deprecations

@Base.deprecate_binding shlib_ext Libdl.dlext

const has_sudo = Ref{Bool}(false)
function __init__()
    if lowercase(get(Base.ENV, "NOSUDO", "false")) == "false"
        has_sudo[] = try success(`sudo -V`) catch err false end
    end
end

end


text:
import Markdown

convert(::Type{Tile}, md::Markdown.MD) = blocktile(md)

blocktile(md::Markdown.MD) = vbox(map(blocktile, md.content))
inlinetile(x::String) = x

blocktile{n}(md::Markdown.Header{n}) = heading(n, map(inlinetile, md.text))

inlinetile(md::Markdown.Code) = code(md.code, language=md.language)
blocktile(md::Markdown.Code)  = codeblock(md.code, language=md.language)

blocktile(md::Markdown.BlockQuote) = blockquote(map(blocktile, md.content))
blocktile(md::Markdown.List) = list(map(item -> map(inlinetile, item), md.items), ordered=md.ordered)
blocktile(md::Markdown.Paragraph) = paragraph(map(inlinetile, md.content))

inlinetile(md::Markdown.Italic) = emph(map(inlinetile, md.text))
inlinetile(md::Markdown.Bold) = font(bold, map(inlinetile, md.text))
inlinetile(md::Markdown.Link) = link(md.url, map(inlinetile, md.text))

inlinetile(md::Markdown.Image) = img(md.url, alt=md.alt)
blocktile(md::Markdown.Image) = img(md.url, alt=md.alt)

inlinetile(md::Markdown.LaTeX) = latex(md.formula, block=false)
blocktile(md::Markdown.LaTeX) = latex(md.formula, block=true)

ext: jl


text:
module ShaleDrillingLikelihood_Production_Test

using ShaleDrillingLikelihood
using Test
using StatsFuns
# using Calculus
using FiniteDiff: finite_difference_gradient
using Optim
using Random
using LinearAlgebra
using BenchmarkTools
using ForwardDiff

using ShaleDrillingLikelihood: _num_x,
    idx_produce_ψ, idx_produce_β, idx_produce_σ2η, idx_produce_σ2u,
    theta_produce, theta_produce_ψ, theta_produce_β, theta_produce_σ2η, theta_produce_σ2u,
    simloglik_produce!, grad_simloglik_produce!, loglik_produce_scalars,
    DataProduce, ObservationGroupProduce, ObservationProduce,
    _x, _y, _xsum, _nu, _i, update_nu!, update_xpnu!, _qm, _psi2,
    SimulationDraws, psi2_wtd_sum_and_sumsq, Observation,
    _nparm, update!

fdgradient(f, x0) = finite_difference_gradient(f, x0, Val{:central})

println("Starting production likelihood tests")

@testset "Production simulated likelihood" begin

    Random.seed!(1234)

    k = 3
    beta = rand(k)
    sigmas = (0.5, 0.3, 0.4,)
    theta = vcat(sigmas[1], beta, sigmas[2:3]...)
    num_i = 10_000
    M = 500

    data = DataProduce(num_i, 10, 10:20, theta)
    allsim = SimulationDraws(M,num_i)

    @test length(theta) == _nparm(data)
    @test _num_x(data) == k

    @test     theta_produce_ψ(  data,theta)  == theta[1]
    @test all(theta_produce_β(  data,theta) .== theta[1 .+ (1:k)])
    @test     theta_produce_σ2η(data,theta)  == theta[end-1]
    @test     theta_produce_σ2u(data,theta)  == theta[end]

    ff(x)          = grad_simloglik_produce!(zeros(length(x)), data, x, allsim, false)
    ffgg!(grad, x) = grad_simloglik_produce!(grad,             data, x, allsim, true)

    tmpgrad = similar(theta)
    ff(theta)

    fd = fdgradient(ff, theta)
    # ad = ForwardDiff.gradient(ff, theta)
    fill!(tmpgrad, 0)
    ffgg!(tmpgrad, theta)
    # @test ad ≈ tmpgrad
    @test norm(fd .- tmpgrad, Inf) < 2e-3
    # @test norm(ad .- tmpgrad, Inf) < 2e-5
    @show fd .- tmpgrad

    od = OnceDifferentiable(ff, ffgg!, ffgg!, theta)
    res = optimize(od, theta, BFGS(), Optim.Options(time_limit = 40.0))
    @test norm(res.minimizer .- theta, Inf) < 0.02
end

println("Beginning speed tests")

M = 500
num_i = 2000
sigmas = (0.5, 0.3, 0.4,)
beta = randn(3)
theta = vcat(sigmas[1], beta, sigmas[2:3]...)
data = DataProduce(num_i, 10, 10:20, theta)
sim = SimulationDraws(M, num_i)
simi = view(sim, 1)
_qm(sim) .= softmax(randn(M))

# @show @benchmark psi2_wtd_sum_and_sumsq(simi)
psi2_wtd_sum_and_sumsq(simi)

obs = Observation(data,1)
simloglik_produce!(obs, theta, simi)

end # module


text:
#
#  CreateRunDirectory.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 02/05/2020.
#
#

module CreateRunDirectory

# Import Julia packages
using Dates

# Import local program modules

function createRunDirectory(nMonomers,L,nDomains,μ,kT,ϵLJ,σ,k,rₑ,Ebend,D,tMax,outputInterval,nGrid,boxSize,intrctnThrshld)

    # Create directory for run data labelled with current time.
    foldername = Dates.format(Dates.now(),"yyyy-mm-dd-HH-MM-SS")
    mkdir("output/$(foldername)")
    # Store system parameters.
    open("output/$(foldername)/conditions.txt","w") do conditionsfile
        println(conditionsfile,"nMonomers      $nMonomers      ")
        println(conditionsfile,"L              $L              ")
        println(conditionsfile,"σ              $σ              ")
        println(conditionsfile,"nDomains       $nDomains       ")
        println(conditionsfile,"nGrid          $nGrid          ")
        println(conditionsfile,"intrctnThrshld $intrctnThrshld ")
        println(conditionsfile,"μ              $μ              ")
        println(conditionsfile,"kT             $kT             ")
        println(conditionsfile,"ϵLJ            $ϵLJ            ")
        println(conditionsfile,"k              $k              ")
        println(conditionsfile,"rₑ             $rₑ             ")
        println(conditionsfile,"Ebend          $Ebend          ")
        println(conditionsfile,"D              $D              ")
        println(conditionsfile,"tMax           $tMax           ")
        println(conditionsfile,"outputInterval $outputInterval ")
        println(conditionsfile,"boxSize        $boxSize        ")
    end

    return foldername
end

export createRunDirectory

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
# using LazySets

function simplify_overapprox(P::LazySet)
    return overapproximate(P, Zonotope, OctDirections(dim(P)))
end

function simplify_underapprox(P::LazySet; ignore::Bool=false)
    if ignore || P isa EmptySet
        return P
    end
    return underapproximate(P, OctDirections(dim(P)))
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
using BenchmarkTools, BioSimulator, StatsBase

function summarize(results, benchmark_set)

  for m in benchmark_set
    fname = Threads.nthreads() > 1 ? "biosimulator-parallel/$(m).txt" : "biosimulator-serial/$(m).txt"

    file = open(fname, "a")

    for a in ALGOS
      for o in OPTIONS
        result = results[m][a][o]

        println(file, "model: $(m) / algorithm: $(a) / output: $(o)")
        println(file)
        println(file, "n = $(length(result.times))")
        println(file, "mean:               $(mean(result.times) * 1e-6) ms")
        println(file, "std:                $(std(result.times) * 1e-6) ms")
        println(file, "median:             $(median(result.times) * 1e-6) ms")
        println(file, "mad (normal=true):  $(mad(result.times, normalize=true) * 1e-6) ms")
        println(file, "mad (normal=false): $(mad(result.times, normalize=false) * 1e-6) ms")
        println(file)
      end
    end
    close(file)
  end

  return nothing
end

# initialize test suite
suite = BenchmarkGroup()

# load the following variables:
# - parameters
# - MODELS
# - SEED
include("benchmark_parameters.jl")

# build handy collections
const algtypemap = Dict(
                  "dm"  => Direct(),
                  "frm" => FirstReaction(),
                  "nrm" => NextReaction(),
                  "odm" => OptimizedDirect(),
                  "otl" => TauLeaping(),
                  "sal" => StepAnticipation()
                )

const outtypemap = Dict("fixed" => Val(:fixed), "full" => Val(:full))

const ALGOS   = collect(keys(algtypemap))
const OPTIONS = collect(keys(outtypemap))
const benchmark_set = intersect(ARGS, MODELS)

# build up the benchmark suite
for m in benchmark_set
  # build the set for the model m
  suite[m] = BenchmarkGroup()

  for a in ALGOS
    # build the set for the algorithm a
    suite[m][a] = BenchmarkGroup()

    # select algorithm
    algtype = algtypemap[a]

    for o in OPTIONS
      # select output option
      outtype = outtypemap[o]

      # unpack parameters
      t_final, n_saves, n_trial, n_sample, t_limit = parameters[m]

      # load model
      include(joinpath(m, "biosimulator.jl"))

      suite[m][a][o] =
        @benchmarkable simulate($model, $algtype, $outtype,
          time   = $t_final,
          epochs = $n_saves,
          trials = $n_trial) samples=n_sample seconds=t_limit setup=(srand($SEED))
    end
  end
end

results = run(suite, verbose = true)

# summarize(results, benchmark_set)

fname = Threads.nthreads() > 1 ? "biosimulator-parallel.json" : "biosimulator-serial.json"

BenchmarkTools.save(fname, results)


text:
using CondaPkg
using Test

# Only run the gc tests on CI (because it's annoying to do it locally)
const testgc = get(ENV, "CI", "") == "true"

# output more than usual when testing
ENV["JULIA_CONDAPKG_VERBOSITY"] = "0"

status() = sprint(io -> CondaPkg.status(io=io))

const backend = get(ENV, "JULIA_CONDAPKG_BACKEND", "MicroMamba")

const isnull = backend == "Null"

@testset "CondaPkg" begin

    @testset "internals" begin

        @testset "PkgSpec" begin
            @test_throws Exception CondaPkg.PkgSpec("")
            @test_throws Exception CondaPkg.PkgSpec("foo!")
            @test_throws Exception CondaPkg.PkgSpec("foo", version="foo")
            spec = CondaPkg.PkgSpec("  F...OO_-0  ", version="  =1.2.3  ", channel="  SOME_chaNNEL  ")
            @test spec.name == "f...oo_-0"
            @test spec.version == "=1.2.3"
            @test spec.channel == "SOME_chaNNEL"
        end

        @testset "ChannelSpec" begin
            @test_throws Exception CondaPkg.ChannelSpec("")
            spec = CondaPkg.ChannelSpec("  SOME_chaNNEL  ")
            @test spec.name == "SOME_chaNNEL"
        end

        @testset "PipPkgSpec" begin
            @test_throws Exception CondaPkg.PipPkgSpec("")
            @test_throws Exception CondaPkg.PipPkgSpec("foo!")
            @test_throws Exception CondaPkg.PipPkgSpec("foo", version="1.2")
            spec = CondaPkg.PipPkgSpec("  F...OO_-0  ", version="  @./SOME/Path  ")
            @test spec.name == "f-oo-0"
            @test spec.version == "@./SOME/Path"
        end

        @testset "meta IO" begin
            specs = Any[
                CondaPkg.PkgSpec("foo", version="=1.2.3", channel="bar"),
                CondaPkg.ChannelSpec("fooo"),
                CondaPkg.PipPkgSpec("foooo", version="==2.3.4")
            ]
            for spec in specs
                io = IOBuffer()
                CondaPkg.write_meta(io, spec)
                seekstart(io)
                spec2 = CondaPkg.read_meta(io, typeof(spec))
                @test spec == spec2
                for k in propertynames(spec)
                    @test getproperty(spec, k) == getproperty(spec2, k)
                end
            end
        end

        @testset "abspathurl" begin
            @test startswith(CondaPkg.abspathurl("foo"), "file://")
            @test endswith(CondaPkg.abspathurl("foo"), "/foo")
            if Sys.iswindows()
                @test CondaPkg.abspathurl("D:\\Foo\\bar.TXT") == "file:///D:/Foo/bar.TXT"
            else
                @test CondaPkg.abspathurl("/Foo/bar.TXT") == "file:///Foo/bar.TXT"
            end
        end

    end

    isnull && @testset "Null backend" begin
        @test CondaPkg.backend() == :Null
        @test CondaPkg.activate!(copy(ENV)) == ENV
        @test occursin("Null", status())
        @test_throws ErrorException CondaPkg.envdir()
    end

    @testset "add/remove channel" begin
        @test !occursin("conda-forge", status())
        CondaPkg.add_channel("conda-forge")
        @test occursin("conda-forge", status())
        CondaPkg.rm_channel("conda-forge")
        @test !occursin("conda-forge", status())
    end

    @testset "install python" begin
        @test !occursin("python", status())
        CondaPkg.add("python", version="==3.10.2")
        isnull || CondaPkg.withenv() do
            pythonpath = joinpath(CondaPkg.envdir(), Sys.iswindows() ? "python.exe" : "bin/python")
            @test isfile(pythonpath)
        end
        @test occursin("python", status())
    end

    @testset "install/remove python package" begin
        # verify package isn't already installed
        @test !occursin("six", status())
        CondaPkg.withenv() do
            isnull || @test_throws Exception run(`python -c "import six"`)
        end

        # install package
        CondaPkg.add("six", version="==1.16.0")
        @test occursin("six", status())
        @test occursin("(==1.16.0)", status())
        CondaPkg.withenv() do
            isnull || run(`python -c "import six"`)
        end
        @test occursin("v1.16.0", status()) == !isnull

        # remove package
        CondaPkg.rm("six")
        @test !occursin("six", status())
        CondaPkg.withenv() do
            isnull || @test_throws Exception run(`python -c "import six"`)
        end
    end

    @testset "pip install/remove python package" begin
        # verify package isn't already installed
        @test !occursin("six", status())
        CondaPkg.withenv() do
            isnull || @test_throws Exception run(`python -c "import six"`)
        end

        # install package
        CondaPkg.add_pip("six", version="==1.16.0")
        @test occursin("six", status())
        @test occursin("(==1.16.0)", status())
        CondaPkg.withenv() do
            isnull || run(`python -c "import six"`)
        end
        @test occursin("v1.16.0", status()) == !isnull

        # remove package
        CondaPkg.rm_pip("six")
        @test !occursin("six", status())
        CondaPkg.withenv() do
            isnull || @test_throws Exception run(`python -c "import six"`)
        end
    end

    isnull || @testset "install/remove executable package" begin
        # install package, verify that executable exists
        CondaPkg.add("curl")
        curl_path = CondaPkg.which("curl")
        @test isfile(curl_path)

        # test uninstall
        CondaPkg.rm("curl")
        CondaPkg.resolve(force=true)
        @test !isfile(curl_path)
    end

    if testgc
        @testset "gc()" begin
            # verify that micromamba clean runs without errors
            CondaPkg.gc()
        end
    end

    @testset "validation" begin
        @test_throws Exception CondaPkg.add("!invalid!package!")
        @test_throws Exception CondaPkg.add_pip("!invalid!package!")
        @test_throws Exception CondaPkg.add("valid-package", version="*invalid*version*")
        @test_throws Exception CondaPkg.add_pip("valid-package", version="*invalid*version*")
        @test_throws Exception CondaPkg.add_channel("")
        @test !occursin("valid", status())
    end

    @testset "PkgREPL" begin
        # The functions in CondaPkg.PkgREPL are not API
        # but calling them is the closest thing to calling
        # the Pkg REPL.

        # add
        CondaPkg.PkgREPL.add(["  six ==1.16.0 "])
        @test occursin("six", status())
        @test occursin("(==1.16.0)", status())

        # rm
        CondaPkg.PkgREPL.rm([" six "])
        @test !occursin("six", status())

        # add channel
        CondaPkg.PkgREPL.channel_add([" numba "])
        @test occursin("numba", status())

        # rm channel
        CondaPkg.PkgREPL.channel_rm([" numba "])
        @test !occursin("numba", status())

        # add pip
        CondaPkg.PkgREPL.pip_add([" six ==1.16.0 "])
        @test occursin("six", status())
        @test occursin("(==1.16.0)", status())

        # rm pip
        CondaPkg.PkgREPL.pip_rm([" six "])
        @test !occursin("six", status())

        # status
        # TODO: capture the output and check it equals status()
        CondaPkg.PkgREPL.status()

        # resolve
        CondaPkg.PkgREPL.resolve()

        # gc
        testgc && CondaPkg.PkgREPL.gc()

        # run
        # TODO: capture the output and check it contains "Python 3.10.2"
        isnull || CondaPkg.PkgREPL.run(["python", "--version"])
    end

end


text:
using jInv.Utils
using Base.Test
using jInv.ForwardShare
using jInv.Mesh
using OED

n = 100
nrhs = 2
# A = abs(sprandn(N,n,.2))
A = sdiag(1./linspace(1,n,n).^4)
A = A[1:2:end,:]
N = size(A,1)
x = rand(n,nrhs); x[4:7,:] = 1;
D = A*x  + 1e-1* randn(N,nrhs)/1e3
OEDparamA = OEDTypeAParam(A,D,zeros(0,0))

# setup OED type B param
r = 20;
l = 8;
Mimg  = getRegularMesh([-1 1 -1 1.],[10,10])
Mgrid = getRegularMesh([-1.2 2.4 -1.6 1.4],[20,20])
OEDparamB = getOEDTypeBParam(r,l,x,Mimg,Mgrid,noiseLevel=0.)

OEDs = (OEDparamA, OEDparamB)

println("--- solve inequality constrained problem ---")
bounds = [.3 .9]
L = .1*speye(n)
Ce = ones(1,n)
ce = Ce*x
Ci = [speye(n);-speye(n)]
ci = [bounds[1]*ones(n,nrhs); -bounds[2]*ones(n,nrhs)]
pForBox    = getIneqconstrainedOEDParam(OEDparamA, L,Ce,ce,Ci,ci,maxIterIPQP=40,outIPQP=-1)
pForNonNeg = getIneqconstrainedOEDParam(OEDparamA, L,[],[],speye(n),0*ones(n,nrhs),maxIterIPQP=40,outIPQP=-1)

pFors = (pForBox, pForNonNeg)
for j=1:length(OEDs)
    for k=1:length(pFors)
		pFors[k].Design = OEDs[j]
		display(pFors[k])

		np  =  getNumberOfDesignParameters(OEDs[j])
		p0      = rand(np)
		xr, = getData(p0,pFors[k])
		@test norm(pFors[k].Ce*xr-pFors[k].ce) < 1e-12
		@test all(pFors[k].Ci*xr.>= pFors[k].ci)
		println("rel.error: $(norm(xr-x)/norm(x))")

		println("--- test derivative of LS solution w.r.t. design ---")
		checkDerivative(rand(np),pFors[k])

		println("--- adjoint test ---")
		(nw,nv) = getSensMatSize(pFors[k])
		v = randn(nv)
		w = randn(nw)
		t1 = dot(w,getSensMatVec(v,p0,pFors[k]))
		t2 = dot(v,getSensTMatVec(w,p0,pFors[k]))
		@test abs(t1-t2)/abs(t1) < 1e-10
	end
end


text:
using DataFrames
using GLM
using RegressionDiscontinuity
using Test

Zs = randn(1000)

ZsR = RunningVariable(Zs, cutoff=0.0)
Zsgeg0 = RunningVariable(Zs, cutoff=1.0)
Zsless0 = RunningVariable(Zs, cutoff=1.0, treated=:<)
@test all(xor.(Zsgeg0.Ws + Zsless0.Ws) .== 1)

tmp_df = DataFrame(Ws=ZsR.Ws, Zs=ZsR.Zs, ZsC=ZsR.ZsC)
@test Tables.schema(tmp_df) == Tables.schema(ZsR)

Ys = rand(1000)
rdd_data = RegressionDiscontinuity.RDData(Ys, ZsR)
myfit = fit(LinearModel, @formula(Ys ~ Zs + Ws), rdd_data)


text:
import Base.Threads.@threads
using DataFrames
"""
    Nucleotide_Essentials.FastqRecord

# Components

* ID: The unique sequence identifier associated with that entry
* sequence: The nucleotide sequence of that entry
* quality: The quality scores of that entry
* filename: The original file name
"""
struct FastqRecord
    ID::Vector{SubString{String}}
    sequence::Vector{SubString{String}}
    quality::Vector{SubString{String}}
    dfquality::DataFrame
    #filename::SubString{String}
end

"""
    Nucleotide_Essentials.readFastq

    readFastq(Path::String)
.fastq file => readFastq(Path) => FastqRecord(ID, sequence, quality, filename)

Supported keyword arguments include:

* 'Path::String': The full or relative path to a .fastq file

# Example:

```julia
# Supply the path to a .Fastq file that you would like to import
myfastq = readFastq("myfastq.fastq")
```
"""
function readFastq(Path::String)
    size = filesize(Path)
    #file = rsplit(Path, "/", limit = 2)[2];
    input = read(Path, String);
    input2 = replace(input, "\n@" => "Z");
    input3 = replace(input2, "@" => "Z", count = 1);
    input4 = replace(input3, "Z" => "\n");
    split_fastq = split(input4, "\n");
    og_quals = split_fastq[5:4:end];
    string_quals = og_quals
    max_length = maximum(length.(og_quals));
    if occursin("^", join(og_quals)) == true && occursin("a", join(og_quals)) == true && occursin("]", join(og_quals)) == true && occursin("f", join(og_quals)) == true
        encoding = 64;
        println("Detected Phred+64 quality encoding")
    end
    if occursin("^", join(og_quals)) == false && occursin("a", join(og_quals)) == false && occursin("]", join(og_quals)) == false && occursin("f", join(og_quals)) == false
        println("Detected Phred+33 quality encoding")
        encoding = 33;
    end

    if size > 50000000
        merged_df = DataFrame()
        println("Large file detected, beginning staged analysis")
        splitamount = Int64(round(length(og_quals)/50))
        og_quals_sub1 = og_quals[1:splitamount]
        og_quals_sub2 = og_quals[(splitamount+1):(2*splitamount)]
        og_quals_sub3 = og_quals[(2*splitamount+1):(3*splitamount)]
        og_quals_sub4 = og_quals[(3*splitamount+1):(4*splitamount)]
        og_quals_sub5 = og_quals[(4*splitamount+1):(5*splitamount)]
        og_quals_sub6 = og_quals[(5*splitamount+1):(6*splitamount)]
        og_quals_sub7 = og_quals[(6*splitamount+1):(7*splitamount)]
        og_quals_sub8 = og_quals[(7*splitamount+1):(8*splitamount)]
        og_quals_sub9 = og_quals[(8*splitamount+1):(9*splitamount)]
        og_quals_sub10 = og_quals[(9*splitamount+1):(10*splitamount)]
        og_quals_sub11 = og_quals[(10*splitamount+1):(11*splitamount)]
        og_quals_sub12 = og_quals[(11*splitamount+1):(12*splitamount)]
        og_quals_sub13 = og_quals[(12*splitamount+1):(13*splitamount)]
        og_quals_sub14 = og_quals[(13*splitamount+1):(14*splitamount)]
        og_quals_sub15 = og_quals[(14*splitamount+1):(15*splitamount)]
        og_quals_sub16 = og_quals[(15*splitamount+1):(16*splitamount)]
        og_quals_sub17 = og_quals[(16*splitamount+1):(17*splitamount)]
        og_quals_sub18 = og_quals[(17*splitamount+1):(18*splitamount)]
        og_quals_sub19 = og_quals[(18*splitamount+1):(19*splitamount)]
        og_quals_sub20 = og_quals[(19*splitamount+1):(20*splitamount)]
        og_quals_sub21 = og_quals[(20*splitamount+1):(21*splitamount)]
        og_quals_sub22 = og_quals[(21*splitamount+1):(22*splitamount)]
        og_quals_sub23 = og_quals[(22*splitamount+1):(23*splitamount)]
        og_quals_sub24 = og_quals[(23*splitamount+1):(24*splitamount)]
        og_quals_sub25 = og_quals[(24*splitamount+1):(25*splitamount)]
        og_quals_sub26 = og_quals[(25*splitamount+1):(26*splitamount)]
        og_quals_sub27 = og_quals[(26*splitamount+1):(27*splitamount)]
        og_quals_sub28 = og_quals[(27*splitamount+1):(28*splitamount)]
        og_quals_sub29 = og_quals[(28*splitamount+1):(29*splitamount)]
        og_quals_sub30 = og_quals[(29*splitamount+1):(30*splitamount)]
        og_quals_sub31 = og_quals[(30*splitamount+1):(31*splitamount)]
        og_quals_sub32 = og_quals[(31*splitamount+1):(32*splitamount)]
        og_quals_sub33 = og_quals[(32*splitamount+1):(33*splitamount)]
        og_quals_sub34 = og_quals[(33*splitamount+1):(34*splitamount)]
        og_quals_sub35 = og_quals[(34*splitamount+1):(35*splitamount)]
        og_quals_sub36 = og_quals[(35*splitamount+1):(36*splitamount)]
        og_quals_sub37 = og_quals[(36*splitamount+1):(37*splitamount)]
        og_quals_sub38 = og_quals[(37*splitamount+1):(38*splitamount)]
        og_quals_sub39 = og_quals[(38*splitamount+1):(39*splitamount)]
        og_quals_sub40 = og_quals[(39*splitamount+1):(40*splitamount)]
        og_quals_sub41 = og_quals[(40*splitamount+1):(41*splitamount)]
        og_quals_sub42 = og_quals[(41*splitamount+1):(42*splitamount)]
        og_quals_sub43 = og_quals[(42*splitamount+1):(43*splitamount)]
        og_quals_sub44 = og_quals[(43*splitamount+1):(44*splitamount)]
        og_quals_sub45 = og_quals[(44*splitamount+1):(45*splitamount)]
        og_quals_sub46 = og_quals[(45*splitamount+1):(46*splitamount)]
        og_quals_sub47 = og_quals[(46*splitamount+1):(47*splitamount)]
        og_quals_sub48 = og_quals[(47*splitamount+1):(48*splitamount)]
        og_quals_sub49 = og_quals[(48*splitamount+1):(49*splitamount)]
        og_quals_sub50 = og_quals[(49*splitamount+1):end]

        subgroups = [og_quals_sub1,
        og_quals_sub2,
        og_quals_sub3,
        og_quals_sub4,
        og_quals_sub5,
        og_quals_sub6,
        og_quals_sub7,
        og_quals_sub8,
        og_quals_sub9,
        og_quals_sub10,
        og_quals_sub11,
        og_quals_sub12,
        og_quals_sub13,
        og_quals_sub14,
        og_quals_sub15,
        og_quals_sub16,
        og_quals_sub17,
        og_quals_sub18,
        og_quals_sub19,
        og_quals_sub20,
        og_quals_sub21,
        og_quals_sub22,
        og_quals_sub23,
        og_quals_sub24,
        og_quals_sub25,
        og_quals_sub26,
        og_quals_sub27,
        og_quals_sub28,
        og_quals_sub29,
        og_quals_sub30,
        og_quals_sub31,
        og_quals_sub32,
        og_quals_sub33,
        og_quals_sub34,
        og_quals_sub35,
        og_quals_sub36,
        og_quals_sub37,
        og_quals_sub38,
        og_quals_sub39,
        og_quals_sub40,
        og_quals_sub41,
        og_quals_sub42,
        og_quals_sub43,
        og_quals_sub44,
        og_quals_sub45,
        og_quals_sub46,
        og_quals_sub47,
        og_quals_sub48,
        og_quals_sub49,
        og_quals_sub50]

        for s in 1:length(subgroups)
            println(string("Beginning batched reading stage ", s, " of ", length(subgroups)))
            if encoding == 33
                og_quals = replace.(subgroups[s],
                            "!" => "00?",
                            "\"" => "01?",
                            "#" => "02?",
                            "\$" => "03?",
                            "%" => "04?",
                            "&" => "05?",
                            "'" => "06?",
                            "(" => "07?",
                            ")" => "08?",
                            "*" => "09?",
                            "+" => "10?",
                            "," => "11?",
                            "-" => "12?",
                            "." => "13?",
                            "/" => "14?",
                            "0" => "15?",
                            "1" => "16?",
                            "2" => "17?",
                            "3" => "18?",
                            "4" => "19?",
                            "5" => "20?",
                            "6" => "21?",
                            "7" => "22?",
                            "8" => "23?",
                            "9" => "24?",
                            ":" => "25?",
                            ";" => "26?",
                            "<" => "27?",
                            "=" => "28?",
                            ">" => "29?",
                            "?" => "30?",
                            "@" => "31?",
                            "A" => "32?",
                            "B" => "33?",
                            "C" => "34?",
                            "D" => "35?",
                            "E" => "36?",
                            "F" => "37?",
                            "G" => "38?",
                            "H" => "39?",
                            "I" => "40?",
                            "J" => "41?",
                            "K" => "42?")
            end
            if encoding == 64
                og_quals = replace.(subgroups[s],
                            "@" => "00?",
                            "A" => "01?",
                            "B" => "02?",
                            "C" => "03?",
                            "D" => "04?",
                            "E" => "05?",
                            "F" => "06?",
                            "G" => "07?",
                            "H" => "08?",
                            "I" => "09?",
                            "J" => "10?",
                            "K" => "11?",
                            "L" => "12?",
                            "M" => "13?",
                            "N" => "14?",
                            "O" => "15?",
                            "P" => "16?",
                            "Q" => "17?",
                            "R" => "18?",
                            "S" => "19?",
                            "T" => "20?",
                            "U" => "21?",
                            "V" => "22?",
                            "W" => "23?",
                            "X" => "24?",
                            "Y" => "25?",
                            "Z" => "26?",
                            "[" => "27?",
                            "\\" => "28?",
                            "]" => "29?",
                            "^" => "30?",
                            "_" => "31?",
                            "`" => "32?",
                            "a" => "33?",
                            "b" => "34?",
                            "c" => "35?",
                            "d" => "36?",
                            "e" => "37?",
                            "f" => "38?",
                            "g" => "39?",
                            "h" => "40?",
                            "i" => "41?",
                            "j" => "42?")
            end
            og_quals = Array(split.(og_quals, "?"))
            pop!.(og_quals)
            df = DataFrame([String[] for i in 1:max_length], :auto)
            for i in 1:length(og_quals)
                #println(string("starting loop #", i, " of ", length(og_quals)))
                temp = collect(og_quals[i])
                if length(temp) < max_length
                    difference = max_length - length(temp)
                    for r in 1:difference
                        append!(temp, ["NaN"])
                    end
                end
                parse.(Float64, temp)
                push!(df, temp)
            end
            df = passmissing(parse).(Float64, df);
            append!(merged_df,df)
        end
        for col in eachcol(merged_df)
            replace!(col,NaN => Float64(100))
        end
        out = FastqRecord(split_fastq[2:4:end], split_fastq[3:4:end], string_quals, merged_df)#, file)
    else
        if encoding == 33
            og_quals = replace.(og_quals,
                        "!" => "00?",
                        "\"" => "01?",
                        "#" => "02?",
                        "\$" => "03?",
                        "%" => "04?",
                        "&" => "05?",
                        "'" => "06?",
                        "(" => "07?",
                        ")" => "08?",
                        "*" => "09?",
                        "+" => "10?",
                        "," => "11?",
                        "-" => "12?",
                        "." => "13?",
                        "/" => "14?",
                        "0" => "15?",
                        "1" => "16?",
                        "2" => "17?",
                        "3" => "18?",
                        "4" => "19?",
                        "5" => "20?",
                        "6" => "21?",
                        "7" => "22?",
                        "8" => "23?",
                        "9" => "24?",
                        ":" => "25?",
                        ";" => "26?",
                        "<" => "27?",
                        "=" => "28?",
                        ">" => "29?",
                        "?" => "30?",
                        "@" => "31?",
                        "A" => "32?",
                        "B" => "33?",
                        "C" => "34?",
                        "D" => "35?",
                        "E" => "36?",
                        "F" => "37?",
                        "G" => "38?",
                        "H" => "39?",
                        "I" => "40?",
                        "J" => "41?",
                        "K" => "42?")
        end
        if encoding == 64
            og_quals = replace.(og_quals,
                        "@" => "00?",
                        "A" => "01?",
                        "B" => "02?",
                        "C" => "03?",
                        "D" => "04?",
                        "E" => "05?",
                        "F" => "06?",
                        "G" => "07?",
                        "H" => "08?",
                        "I" => "09?",
                        "J" => "10?",
                        "K" => "11?",
                        "L" => "12?",
                        "M" => "13?",
                        "N" => "14?",
                        "O" => "15?",
                        "P" => "16?",
                        "Q" => "17?",
                        "R" => "18?",
                        "S" => "19?",
                        "T" => "20?",
                        "U" => "21?",
                        "V" => "22?",
                        "W" => "23?",
                        "X" => "24?",
                        "Y" => "25?",
                        "Z" => "26?",
                        "[" => "27?",
                        "\\" => "28?",
                        "]" => "29?",
                        "^" => "30?",
                        "_" => "31?",
                        "`" => "32?",
                        "a" => "33?",
                        "b" => "34?",
                        "c" => "35?",
                        "d" => "36?",
                        "e" => "37?",
                        "f" => "38?",
                        "g" => "39?",
                        "h" => "40?",
                        "i" => "41?",
                        "j" => "42?")
        end
        og_quals = Array(split.(og_quals, "?"))
        pop!.(og_quals)
        df = DataFrame([String[] for i in 1:max_length], :auto)
        for i in 1:length(og_quals)
            println(string("starting loop #", i, " of ", length(og_quals)))
            temp = collect(og_quals[i])
            if length(temp) < max_length
                difference = max_length - length(temp)
                for r in 1:difference
                    append!(temp, ["NaN"])
                end
            end
            parse.(Float64, temp)
            push!(df, temp)
        end
        df = passmissing(parse).(Float64, df);
        for col in eachcol(df)
            replace!(col,NaN => Float64(100))
        end
        out = FastqRecord(split_fastq[2:4:end], split_fastq[3:4:end], string_quals, df)#, file)
    end
    return out
end # end function readFastq

# Some benchmarking timing stats based on .fastq size
    # 2 MB - 4.049943 seconds (54.18 M allocations: 3.418 GiB, 17.70% gc time, 0.11% compilation time)
    # 34.8 MB - 61.989768 seconds (900.48 M allocations: 56.452 GiB, 10.68% gc time)
    # 33.2 MB - 53.795583 seconds (825.51 M allocations: 51.120 GiB, 8.06% gc time)
    # 1.41 GB - 2646.207447 seconds (32.56 G allocations: 1.949 TiB, 21.35% gc time)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractAggregation end
abstract type AbstractAggregator end
"""
    RawAggregationParameter

An optimizer attribute used for raw parameters of the aggregator. Defers to `RawParameter`.
"""
struct RawAggregationParameter <: AggregationParameter
    name::Any
end

function MOI.get(aggregator::AbstractAggregator, param::RawAggregationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(aggregator)))
        error("Unrecognized parameter name: $(name) for aggregator $(typeof(aggregator)).")
    end
    return getfield(aggregator, name)
end

function MOI.set(aggregator::AbstractAggregator, param::RawAggregationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(aggregator)))
        error("Unrecognized parameter name: $(name) for aggregator $(typeof(aggregator)).")
    end
    setfield!(aggregator, name, value)
    return nothing
end

# ------------------------------------------------------------
include("cut_collection.jl")
include("distance_measures.jl")
include("selection_rules.jl")
include("cluster_rules.jl")
include("lock_rules.jl")
include("no_aggregation.jl")
include("partial_aggregation.jl")
include("dynamic_aggregation.jl")
include("cluster_aggregation.jl")
include("granulated_aggregation.jl")
include("hybrid_aggregation.jl")


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Physical_Properties_Heatex
	Physical_Properties_Heatex()=begin
		new(
			molweight(Dict{Symbol,Any}(
				:Brief=>"Average Mol Weight",
				:Default=>75,
				:Lower=>1,
				:Upper=>1e8
			)),
			cp_mol(Dict{Symbol,Any}(
				:Brief=>"Average Molar Heat Capacity",
				:Upper=>1e10
			)),
			Properties_In_Out(Dict{Symbol,Any}(
				:Brief=>"Properties at Inlet Stream",
				:Symbol=>"_{in}"
			)),
			Properties_In_Out(Dict{Symbol,Any}(
				:Brief=>"Properties at Outlet Stream",
				:Symbol=>"_{out}"
			)),
			[:Mw,:Cp,:Inlet,:Outlet,]
		)
	end
	Mw::molweight
	Cp::cp_mol
	Inlet::Properties_In_Out
	Outlet::Properties_In_Out
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Physical_Properties_Heatex
function atributes(in::Physical_Properties_Heatex,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Physical_Properties_Heatex(_::Dict{Symbol,Any})=begin
	newModel=Physical_Properties_Heatex()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module Thrift

using Distributed
using Sockets

import Sockets: TCPServer, listen, accept
import Base: open, close, isopen, read, read!, write, flush, skip, show, copy!

export open, close, isopen, read, read!, write, flush, skip, listen, accept, show, copy!

# from base.jl
export TSTOP, TVOID, TBOOL, TBYTE, TI08, TDOUBLE, TI16, TI32, TI64, TSTRING, TUTF7, TSTRUCT, TMAP, TSET, TLIST, TUTF8, TUTF16
export TType, TProcessor, TTransport, TServerTransport, TServer, TProtocol
export writeMessageBegin, writeMessageEnd, writeStructBegin, writeStructEnd, writeFieldBegin, writeFieldEnd, writeFieldStop, writeMapBegin, writeMapEnd, writeListBegin, writeListEnd, writeSetBegin, writeSetEnd, writeBool, writeByte, writeI16, writeI32, writeI64, writeDouble, writeString
export readMessageBegin, readMessageEnd, readStructBegin, readStructEnd, readFieldBegin, readFieldEnd, readMapBegin, readMapEnd, readListBegin, readListEnd, readSetBegin, readSetEnd, readBool, readByte, readI16, readI32, readI64, readDouble, readString
export ApplicationExceptionType, MessageType, TException, TApplicationException
export ThriftMetaAttribs, ThriftMeta, meta
export isinitialized, set_field, set_field!, get_field, clear, has_field, fillunset, fillset, filled, isfilled, thriftbuild, enumstr

# from transports.jl
export TFramedTransport, TSASLClientTransport, TSocket, TServerSocket, TSocketBase, TMemoryTransport, TFileTransport
export TransportExceptionTypes, TTransportException

# from sasl.jl
export SASL_MECH_PLAIN, SASL_MECH_KERB, SASL_MECH_LDAP, SASLException

# from protocols.jl
export TBinaryProtocol, TCompactProtocol

# from processor.jl
export ThriftProcessor, ThriftHandler, process, handle, extend, distribute

# from server.jl
export TSimpleServer, TTaskServer, TProcessPoolServer, serve

if !isdefined(Base, :fieldtypes)
    fieldtypes(T::Type) = ntuple(i -> fieldtype(T, i), fieldcount(T))
end

include("base.jl")
include("codec.jl")
include("sasl.jl")
include("transports.jl")
include("protocols.jl")
include("processor.jl")
include("server.jl")

end # module


text:
#versioninfo()

#import Pkg; Pkg.add("Arpack")
#import Pkg; Pkg.add("ArgParse")

using LinearAlgebra
using Arpack
using SparseArrays
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--Lx"
            help = "set Lx (should be >=2)"
            arg_type = Int
            default = 4
        "--Ly"
            help = "set Ly (=2)"
            arg_type = Int
            default = 2
        "--Jleg"
            help = "set Jleg"
            arg_type = Float64
            default = 1.0
        "--Jrung"
            help = "set Jrung"
            arg_type = Float64
            default = 1.0
        "--Jising"
            help = "set Jising"
            arg_type = Float64
            default = 1.0
        "--J4"
            help = "set J4"
            arg_type = Float64
            default = 0.0
        "--twoSz"
            help = "set twoSz"
            arg_type = Int
            default = 0
        "--momkx"
            help = "set momkx"
            arg_type = Int
            default = 0
        "--momky"
            help = "set momky"
            arg_type = Int
            default = 0
    end
    return parse_args(s)
end

function get_next_same_nup_state(state)
    next = 0
    if (state>0)
        smallest = state & -(state)
        ripple = state + smallest
        ones = xor(state,ripple)
        ones = div((ones >> 2),smallest)
        next = ripple | ones
    end
    return next
end

function shift_child_BA2AB(s,a,b)
    return ((s<<b)&(((1<<a)-1)<<b))|((s>>a)&((1<<b)-1))
end

function shift_child_CBA2CAB(s,a,b,c)
    left = s&(((1<<c)-1)<<(a+b))
    right = s&((1<<(a+b))-1)
    right = shift_child_BA2AB(right,a,b)
    return left|right
end

function shift_child_DCBA2DBCA(s,a,b,c,d)
    left = s&(((1<<d)-1)<<(a+b+c))
    right = s&((1<<a)-1)
    mid = (s>>a)&((1<<(b+c))-1)
    mid = shift_child_BA2AB(mid,b,c)
    mid = mid<<a
    return left|(mid|right)
end

function shift_y_1spin(state,Lx,Ly,Ns)
    return shift_child_BA2AB(state,Ns-Lx,Lx)
end

function shift_x_1spin(state,Lx,Ly,Ns)
    s = shift_child_CBA2CAB(state,Lx-1,1,Ns-Lx)
    for n in 1:Ly-1
        s = shift_child_DCBA2DBCA(s,n*Lx,Lx-1,1,Ns-(n+1)*Lx)
    end
    return s
end

function check_state(state,Lx,Ly,Ns)
    list_allstate = Int[]
    r = state
    t = state
    push!(list_allstate,t)
    for i in 0:Lx-2
        t = shift_x_1spin(t,Lx,Ly,Ns)
        push!(list_allstate,t)
        if (t < r)
            r = t
        end
    end
    for j in 0:Ly-2
        t = shift_y_1spin(t,Lx,Ly,Ns)
        for i in 0:Lx-1
            t = shift_x_1spin(t,Lx,Ly,Ns)
            push!(list_allstate,t)
            if (t < r)
                r = t
            end
        end
    end
    D = length(unique(list_allstate)) ## remove duplication and count length of a list
#    println(D," ",list_allstate)
    if (r == state)
        return +D
    else
        return -D
    end
end

function calc_exp(Lx,Ly,momkx,momky)
#    return [[exp(-1im*2.0*pi*(expx*momkx/Lx+expy*momky/Ly)) for expy in 0:Ly-1] for expx in 0:Lx-1]
#    return [exp(-1im*2.0*pi*(expx*momkx/Lx+expy*momky/Ly)) for expy in 0:Ly-1, expx in 0:Lx-1]
    @inbounds return [exp(-1im*2.0*pi*(expx*momkx/Lx+expy*momky/Ly)) for expx in 0:Lx-1, expy in 0:Ly-1]
end

function calc_sum_exp(r,Lx,Ly,Ns,expk)
    t = r
#    F2 = expk[0,0]
    F2 = expk[1,1]
    @inbounds for i in 0:Lx-2
        t = shift_x_1spin(t,Lx,Ly,Ns)
        if (t == r)
#            F2 += expk[i+1,0]
            F2 += expk[i+2,1]
        end
    end
    @inbounds for j in 0:Ly-2
        t = shift_y_1spin(t,Lx,Ly,Ns)
        for i in 0:Lx-1
            t = shift_x_1spin(t,Lx,Ly,Ns)
            if (t == r)
#                F2 += expk[i,j+1]
                F2 += expk[i+1,j+2]
            end
        end
    end
    return F2
end

function make_basis(Lx,Ly,Ns,nup,Nbinom,momkx,momky,expk)
#function make_basis(Lx,Ly,Ns,nup,Nbinom,momkx,momky)
    list_state = Int[]
    list_R = Int[]
    list_F2 = Complex{Float64}[]
    first = (1<<(Ns-nup))-1
    last = ((1<<(Ns-nup))-1)<<(nup)
    Nrep = 0
    state = first
    for i in 0:Nbinom-1
        R = check_state(state,Lx,Ly,Ns)
        if (R>=0)
            push!(list_state,state)
            push!(list_R,R)
            Nrep += 1
            F2 = calc_sum_exp(state,Lx,Ly,Ns,expk)
            push!(list_F2,F2)
        end
        state = get_next_same_nup_state(state)
    end
    return list_state, list_R, list_F2, Nrep
#    return list_state, list_R, Nrep
end

function find_representative(state,Lx,Ly,Ns)
    r = state
    t = state
    expx = 0
    expy = 0
    for i in 0:Lx-2
        t = shift_x_1spin(t,Lx,Ly,Ns)
        if (t < r)
            r = t
            expx = i+1
        end
    end
    for j in 0:Ly-2
        t = shift_y_1spin(t,Lx,Ly,Ns)
        for i in 0:Lx-1
            t = shift_x_1spin(t,Lx,Ly,Ns)
            if (t < r)
                r = t
                expx = i
                expy = j+1
            end
        end
    end
    return r, expx, expy
end

function get_spin(state,site)
    return (state>>site)&1
end

function flip_2spins(state,i1,i2)
    return xor(state,((1<<i1)+(1<<i2)))
end

function find_state_2(state,list_1,maxind)
    imin = 0
    imax = maxind-1
    i = div(imin+imax,2) ## make "i" in the loop global
    while true
        i = div(imin+imax,2)
#        if (state < list_1[i])
        if (state < list_1[i+1])
            imax = i-1
#        elseif (state > list_1[i])
        elseif (state > list_1[i+1])
            imin = i+1
        else
            break
        end
        if (imin > imax)
            return -1
        end
    end
    return i
end

function make_hamiltonian_child(Jxx,Jzz,list_site1,list_site2,Nbond1,
    J4xx,J4zz,J4xz,list_site3,list_site4,list_site5,list_site6,Nbond2,
    Nrep,list_state,list_sqrtnorm,Lx,Ly,Ns,momkx,momky,expk)
    Nbond = Nbond1 + 3*Nbond2
#
    listki = ones(Int,(Nbond+1)*Nrep)
    loc = ones(Int,(Nbond+1)*Nrep)
    elemnt = zeros(Complex{Float64},(Nbond+1)*Nrep)
    @inbounds for k in 0:Nbond
        for i in 0:Nrep-1
            listki[k*Nrep+i+1] = i+1
        end
    end
##
## index of loc[]
##   Nbond*Nrep+a: diag
##   i*Nrep+a: offdiag such as [s2,s1] = 01, 10
##   (Nbond1+i)*Nrep+a: offdiag such as [s6,s5,s4,s3] = 0100, 1000, 0111, 1011
##   (Nbond1+Nbond2+i)*Nrep+a: offdiag such as [s6,s5,s4,s3] = 0001, 0010, 1101, 1110
##   (Nbond1+2*Nbond2+i)*Nrep+a: offdiag such as [s6,s5,s4,s3] = 0101, 0110, 1001, 1010
##
    @inbounds for a in 0:Nrep-1
        sa = list_state[a+1]
        loc[Nbond*Nrep+a+1] = a+1
        for i in 0:Nbond1-1
            i1 = list_site1[i+1]
            i2 = list_site2[i+1]
            wght = 2.0*Jxx[i+1]
            diag = Jzz[i+1]
#            loc[Nbond*Nrep+a+1] = a+1
            if get_spin(sa,i1) == get_spin(sa,i2)
                elemnt[Nbond*Nrep+a+1] += diag
            else
                elemnt[Nbond*Nrep+a+1] -= diag
                bb = flip_2spins(sa,i1,i2)
                sb, expx, expy = find_representative(bb,Lx,Ly,Ns)
                b = find_state_2(sb,list_state,Nrep)
                if (b>=0)
                    elemnt[i*Nrep+a+1] += wght*list_sqrtnorm[b+1]/list_sqrtnorm[a+1]*expk[expx+1,expy+1]
                    loc[i*Nrep+a+1] = b+1
                end
            end
        end
        for i in 0:Nbond2-1
            i3 = list_site3[i+1]
            i4 = list_site4[i+1]
            i5 = list_site5[i+1]
            i6 = list_site6[i+1]
            czz = J4zz[i+1]
            cxz = 2.0*J4xz[i+1]
            cxx = 4.0*J4xx[i+1]
#            loc[Nbond*Nrep+a+1] = a+1
            if get_spin(sa,i3) == get_spin(sa,i4) && get_spin(sa,i5) == get_spin(sa,i6)
                elemnt[Nbond*Nrep+a+1] += czz
            elseif get_spin(sa,i3) == get_spin(sa,i4) && get_spin(sa,i5) != get_spin(sa,i6)
                elemnt[Nbond*Nrep+a+1] -= czz
#
                bb = flip_2spins(sa,i5,i6)
                sb, expx, expy = find_representative(bb,Lx,Ly,Ns)
                b = find_state_2(sb,list_state,Nrep)
                if (b>=0)
                    elemnt[(Nbond1+i)*Nrep+a+1] += cxz*list_sqrtnorm[b+1]/list_sqrtnorm[a+1]*expk[expx+1,expy+1]
                    loc[(Nbond1+i)*Nrep+a+1] = b+1
                end
            elseif get_spin(sa,i3) != get_spin(sa,i4) && get_spin(sa,i5) == get_spin(sa,i6)
                elemnt[Nbond*Nrep+a+1] -= czz
#
                bb = flip_2spins(sa,i3,i4)
                sb, expx, expy = find_representative(bb,Lx,Ly,Ns)
                b = find_state_2(sb,list_state,Nrep)
                if (b>=0)
                    elemnt[(Nbond1+Nbond2+i)*Nrep+a+1] += cxz*list_sqrtnorm[b+1]/list_sqrtnorm[a+1]*expk[expx+1,expy+1]
                    loc[(Nbond1+Nbond2+i)*Nrep+a+1] = b+1
                end
            elseif get_spin(sa,i3) != get_spin(sa,i4) && get_spin(sa,i5) != get_spin(sa,i6)
                elemnt[Nbond*Nrep+a+1] += czz
#
                bb = flip_2spins(sa,i3,i4)
                sb, expx, expy = find_representative(bb,Lx,Ly,Ns)
                b = find_state_2(sb,list_state,Nrep)
                if (b>=0)
                    elemnt[(Nbond1+Nbond2+i)*Nrep+a+1] += -cxz*list_sqrtnorm[b+1]/list_sqrtnorm[a+1]*expk[expx+1,expy+1]
                    loc[(Nbond1+Nbond2+i)*Nrep+a+1] = b+1
                end
#
                bb = flip_2spins(sa,i5,i6)
                sb, expx, expy = find_representative(bb,Lx,Ly,Ns)
                b = find_state_2(sb,list_state,Nrep)
                if (b>=0)
                    elemnt[(Nbond1+i)*Nrep+a+1] += -cxz*list_sqrtnorm[b+1]/list_sqrtnorm[a+1]*expk[expx+1,expy+1]
                    loc[(Nbond1+i)*Nrep+a+1] = b+1
                end
#
                bbtmp = flip_2spins(sa,i3,i4)
                bb = flip_2spins(bbtmp,i5,i6)
                sb, expx, expy = find_representative(bb,Lx,Ly,Ns)
                b = find_state_2(sb,list_state,Nrep)
                if (b>=0)
                    elemnt[(Nbond1+2*Nbond2+i)*Nrep+a+1] += cxx*list_sqrtnorm[b+1]/list_sqrtnorm[a+1]*expk[expx+1,expy+1]
                    loc[(Nbond1+2*Nbond2+i)*Nrep+a+1] = b+1
                end
            end
        end
    end
    return elemnt, listki, loc
end

##
## ...- L+1 - L+2 - L+3 -...- 2L -...
##       |     |     |         |
## ...-  0  -  1  -  2  -...-  L -...
##
## - 0 - 2 -     - ix+L - ix+L+1 -
##   |   |   -->     |       |
## - 1 - 3 -     -  ix  -  ix+1  -
##
function make_lattice(Lx,Ly,Jleg,Jrung,Jising,J4)
    Jxx = Float64[]
    Jzz = Float64[]
    J4xx = Float64[]
    J4zz = Float64[]
    J4xz = Float64[]
    list_isite1 = Int[]
    list_isite2 = Int[]
    list_isite3 = Int[]
    list_isite4 = Int[]
    list_isite5 = Int[]
    list_isite6 = Int[]
    Nint = 0
    Nint2 = 0
    @inbounds for ix in 0:Lx-1
        site0 = ix+Lx
        site1 = ix
        site2 = (ix+1)%Lx+Lx
        site3 = (ix+1)%Lx
#
## 2-body terms
        push!(list_isite1,site0)
        push!(list_isite2,site2)
        push!(Jxx,Jleg)
        push!(Jzz,Jleg*Jising)
        Nint += 1
#
        push!(list_isite1,site1)
        push!(list_isite2,site3)
        push!(Jxx,Jleg)
        push!(Jzz,Jleg*Jising)
        Nint += 1
#
        push!(list_isite1,site0)
        push!(list_isite2,site1)
        push!(Jxx,Jrung)
        push!(Jzz,Jrung*Jising)
        Nint += 1
#
## 4-body terms
        push!(list_isite3,site0)
        push!(list_isite4,site2)
        push!(list_isite5,site1)
        push!(list_isite6,site3)
        push!(J4xx,J4)
        push!(J4zz,J4)
        push!(J4xz,J4)
        Nint2 += 1
#
    end
    return Jxx, Jzz, list_isite1, list_isite2, Nint,
        J4xx, J4zz, J4xz, list_isite3, list_isite4, list_isite5, list_isite6, Nint2
end

#function make_hamiltonian(Nrep,elemnt,listki,loc)
function make_hamiltonian(elemnt,listki,loc)
    return dropzeros(sparse(listki,loc,elemnt))
end

function main()
    parsed_args = parse_commandline()
    Lx = parsed_args["Lx"]
    Ly = parsed_args["Ly"]
    twoSz = parsed_args["twoSz"]
    Jleg = parsed_args["Jleg"]
    Jrung = parsed_args["Jrung"]
    Jising = parsed_args["Jising"]
    J4 = parsed_args["J4"]
    momkx = parsed_args["momkx"]
    momky = parsed_args["momky"]
    Ns = Lx*Ly
    nup = div(Ns+twoSz,2)
    Nbinom = binomial(Ns,nup)

    println("# make list")
    println("Lx=",Lx)
    println("Ly=",Ly)
    println("Jleg=",Jleg)
    println("Jrung=",Jrung)
    println("Jising=",Jising)
    println("J4=",J4)
    println("Ns=",Ns)
    println("twoSz=",twoSz)
    println("nup=",nup)
    println("Nbinom=",Nbinom)
    println("momkx=",momkx)
    println("momky=",momky)
    @time expk = calc_exp(Lx,Ly,momkx,momky)
#    println(expk)
#    println(expk[0+1,0+1])
#    println(expk[0+1,1+1])
#    println(expk[1+1,0+1])
#    println(expk[1+1,1+1])
#    println(expk[2+1,0+1])
#    println(expk[2+1,1+1])
#    println(expk[3+1,0+1])
#    println(expk[3+1,1+1])
    @time list_state, list_R, list_F2, Nrep = make_basis(Lx,Ly,Ns,nup,Nbinom,momkx,momky,expk)
#    list_state, list_R, Nrep = make_basis(Lx,Ly,Ns,nup,Nbinom,momkx,momky)
#    println(list_state)
#    println(list_R)
#    println(list_F2)
#    println(abs.(list_F2))
#    println(abs.(list_F2).>1e-12)
    mask = abs.(list_F2).>1e-12
#    println(mask)
    println("Nrep=",Nrep)
    println("# num of of excluded elements from |F|^2:",Nrep-sum(mask))
    Nrep = sum(mask)
    println("Nrep=",Nrep)
    list_state = list_state[mask]
    list_R = list_R[mask]
    list_F2 = list_F2[mask]
    list_sqrtnorm = sqrt.(list_R.*list_F2)
#    println(list_state)
#    println(list_R)
#    println(list_F2)
#    println(list_sqrtnorm)
    println()

    println("# make interactions")
    @time Jxx, Jzz, list_site1, list_site2, Nbond1,
        J4xx, J4zz, J4xz, list_site3, list_site4, list_site5, list_site6, Nbond2 =
        make_lattice(Lx,Ly,Jleg,Jrung,Jising,J4)
##
## make 2-body interactions for 1/4*(sigma^z.sigma^z + 2(S^+.S^- + S^-.S^+))
    Jxx = Jxx/4.0
    Jzz = Jzz/4.0
##
## make 4-body interactions for 1/16*(sigma^z.sigma^z + 2(S^+.S^- + S^-.S^+))^2
    J4xx = J4xx/16.0
    J4zz = J4zz/16.0
    J4xz = J4xz/16.0
##
    println("Jxx=",Jxx)
    println("Jzz=",Jzz)
    println("list_site1=",list_site1)
    println("list_site2=",list_site2)
    println("Nbond1=",Nbond1)
    println("J4xx=",J4xx)
    println("J4zz=",J4zz)
    println("J4xz=",J4xz)
    println("list_site3=",list_site3)
    println("list_site4=",list_site4)
    println("list_site5=",list_site5)
    println("list_site6=",list_site6)
    println("Nbond2=",Nbond2)
    println()

    println("# make Hamiltonian")
    @time elemnt, listki, loc =
        make_hamiltonian_child(Jxx,Jzz,list_site1,list_site2,Nbond1,
        J4xx,J4zz,J4xz,list_site3,list_site4,list_site5,list_site6,Nbond2,
        Nrep,list_state,list_sqrtnorm,Lx,Ly,Ns,momkx,momky,expk)
#    println(elemnt)
#    println(listki)
#    println(loc)
#    @time Ham = make_hamiltonian(Nrep,elemnt,listki,loc)
    @time Ham = make_hamiltonian(elemnt,listki,loc)
#    println(Ham)
    println()

    println("# diag Hamiltonian")
    Neig = 5
#    @time ene,vec = eigs(Ham,nev=min(Neig,Nrep-1),which=:SR)
    @time ene,vec = eigs(Ham,nev=min(Neig,Nrep-1),which=:SR,ritzvec=false)
    ene = sort(real(ene))/Ns
    println("# energy:",ene[1]," ",ene[2]," ",ene[3]," ",ene[4]," ",ene[5])
    println()

end

main()


text:
###########
# Imports #
###########

# Julia Parallel Packages
using Distributed

# Julia Packages
using Dates
using Random
using Profile
using Printf
using Logging
using Statistics
using Distributions

# Register Number of processes
# Julia must be started with `julia -p auto -e "using IJulia; jupyterlab();"`
# On Computer with 6 cores add 6-1 = 5 processes
# NOTE: This must be called _before_ @everywhere is used those functions get propery distributed
if length(procs()) == 1
    # Remove existing process IDs so reloading works
    addprocs(length(Sys.cpu_info()));
    println("Started $(nprocs()) total processes - $(nworkers()) worker processes.")
else
    println("Found $(nprocs()) total processes - $(nworkers()) worker processes.")
end

# Satellite Dynamics Packages
@everywhere using SatelliteDynamics

# Load SatelliteTasking
@everywhere using SatelliteTasking

# Add Map Plots
using PGFMapPlots

#########################
# Simulation Definition #
#########################

# Set simulation time window
t_start = Epoch(2019, 1, 1, 0, 0, 0, tsys="UTC"); # Start of time span
t_end = Epoch(2019, 1, 2, 0, 0, 0, tsys="UTC"); # End of simulation time span

# Initialize Planning Problem
problem = SatPlanningProblem(t_start=t_start, t_end=t_end);

# Add constraints to problem
push!(problem.constraints, constraint_agility_single_axis)

# Create Spacecraft
line1 = "1     1U          19  1.00000000  .00000000  00000-0  00000-0 0    05"
line2 = "2     1  90.0000   0.0000 0010000   0.0000   0.0000 15.24162312    00"
sc1 = Spacecraft(id=1, tle=TLE(line1, line2), slew_rate=1.0);
push!(problem.spacecraft, sc1);

# Compute Orital Period
T = orbit_period(state(problem.spacecraft[1].tle, problem.spacecraft[1].tle.epoch)[1]);

println("Loadded $(length(problem.spacecraft)) spacecraft")

# Set Solve Parameters
problem.solve_discount = 0.99 # Typical values 0.99 - 0.9999
problem.solve_depth = 3
problem.solve_breadth = 3
problem.solve_horizon = T

println("Default Solve Parameters:\nsolve_discount: $(problem.solve_discount)\nsolve_depth: $(problem.solve_depth)\nsolve_breadth: $(problem.solve_breadth)\nsolve_horizon: $(problem.solve_horizon)")

######################
# Load Location Data #
######################

# Load Groundstation Data
stations = load_stations("./data/groundstations.json");
num_stations = length(stations);

# Load Image Data - [100, 200, 500, 1000, 2000, 5000, 10000, 12000, 15000]
num_requests = 1000
requests = load_requests("./data/landsat_test_$num_requests.json", id_offset=num_stations);
@assert num_requests == length(requests);

# Set Request Parameters
collect_duration = 30.0
set_collect_duration(requests, collect_duration)

# Add Locations To Problem
# add_locations(problem, stations)
add_locations(problem, requests)

#########################
# Compute Opportunities #
#########################

# Compute problem accesses
compute_access(problem, orbit_fraction=0.75);

println("")
println("Found $(length(problem.contacts)) passes in planning horizon")
println("Found $(length(problem.collects)) collects in planning horizon")
num_accessible_requests = sum([1 for req in problem.requests if length(problem.lt_loc_opps[req.id]) > 0])
pct_accessible_requests = floor(num_accessible_requests/length(problem.requests) * 10000)/100
println("Found $num_accessible_requests out of $(length(problem.requests)) requests have accesses. $pct_accessible_requests%")

# MDP Planning Preparation - (Pre-compute Actions)
println("Precomputing action space")
precompute_action_space(problem, enable_resources=false);

########################
# Define Work Function #
########################

@everywhere function run_mcts_sim(problem::SatPlanningProblem, n::Int, params::Tuple{Real, Real, Real, Real})
    println("Reporting in!")
    # Create reward storage
    hyper_time = Dict{Tuple{Real, Real, Real, Real}, AbstractVector{Real}}()
    hyper_reward = Dict{Tuple{Real, Real, Real, Real}, AbstractVector{Real}}()

    # Extract search point from params
    ec, ni, sd, ds = params

    # Copy problem
    # problem = copy(problem)
    problem.solve_discount = ds # Typical values 0.99 - 0.9999
    problem.solve_depth = sd
    problem.mcts_n_iterations = ni
    problem.mcts_exploration_constant = ec

    hyper_time[(ec, ni, sd, ds)]   = Real[]
    hyper_reward[(ec, ni, sd, ds)] = Real[]

    # Run specified number of iterations
    for idx in 1:n
        start_time = time()

        # Solve problem
        println("Iteration $idx - Solving Problem with parameters:($(problem.solve_discount), $(problem.solve_depth), $(problem.mcts_n_iterations), $(problem.mcts_exploration_constant))")
        mdp_mcts_states, mdp_mcts_plan, mdp_mcts_reward = satellite_plan_mdp_mcts(problem, parallel=false)
        end_time = time()

        # Extract reward
        mdp_mcts_feasible, mdp_mcts_realized_reward, mdp_mcts_nur, mdp_mcts_ndr, mdp_mcts_nr, mdp_mcts_nc = analyze_plan(problem, mdp_mcts_plan);

        # Update serach values
        push!(hyper_time[(ec, ni, sd, ds)], end_time-start_time)
        push!(hyper_reward[(ec, ni, sd, ds)], mdp_mcts_realized_reward)
    end

    return hyper_time, hyper_reward
end

######################################
# Define Hyperparameter Search Space #
######################################

# Define parameter space
ec_space = [0.01, 0.1, 1.0, 1.5, 2.0, 3.0, 5.0, 10.0]
ni_space = [10, 20, 30, 50, 100, 200, 500, 1000]
sd_space = [2, 3, 5, 10, 20, 30, 50] # solve depth
ds_space = [0.99, 0.995, 0.999] # discount
bin_iterations = 10

ec_space = [1.5, 2.0, 3.0]
ni_space = [10, 20]
sd_space = [10, 20] # solve depth
ds_space = [0.99, 0.995] # discount
bin_iterations = 3

# ec_space = [0.01, 0.1, 0.2]
# ni_space = [5, 10]
# sd_space = [2, 3, 5]
# ds_space = [0.99]
# bin_iterations = 10

# Create reward storage
hyper_time = Dict{Tuple{Real, Real, Real, Real}, AbstractVector{Real}}()
hyper_reward = Dict{Tuple{Real, Real, Real, Real}, AbstractVector{Real}}()

println("Total number of required simulations: $(length(ec_space)*length(ni_space)*length(sd_space)*length(ds_space)*bin_iterations)")
num_jobs = length(ec_space)*length(ni_space)*length(sd_space)*length(ds_space)
println("Total number of jobs: $num_jobs")

#######################
# Run Parallel Search #
#######################

fn(p, n) = x -> run_mcts_sim(p, n, x)

# Create work assignments
nw = nworkers()        # Number of workers
lw = num_jobs # Length of work

wa = floor(Int, lw/nw) # Average work per worker
wr = lw - nw*wa        # Remaining worker

println("Have $lw items of work and $nw workers")
println("Average work: $wa, Remaining work: $wr")

# Create search points
work = []
for ec in ec_space
    for ni in ni_space
        for sd in sd_space
            for ds in ds_space
                push!(work, (ec, ni, sd, ds))
            end
        end
    end
end

println("Work: $work")
# Construct assignments (function inputs) for Workers
# assignments = []
# println("Assignments: $assignments")

# push!(assignments, work[1:(wa+wr)])
# for i in 2:min(nw, length(work))
#     # @debug "Worker $i convering $(1+wr+(i-1)*wa):$(i*wa+wr)")
#     push!(assignments, work[(1+wr+(i-1)*wa):(i*wa+wr)])
# end

# Execute Work in parallel

println("Starting hyper parameter search")
hp_start = time()
results = pmap(fn(problem, bin_iterations), work)
# results = pmap(fn(problem, bin_iterations), assignments)
hp_end = time()
println("Parameter search finished after $(hp_end-hp_start) seconds.")

# Merge results
hyper_time = Dict{Tuple{Real, Real, Real, Real}, AbstractVector{Real}}()
hyper_reward = Dict{Tuple{Real, Real, Real, Real}, AbstractVector{Real}}()

for r in results
    ht, hr = r
    for params, time_values in ht
        hyper_time[params] = time_values
        hyper_reward[params] = hr[params]
    end
end

#################################
# Extract Best Hyper Parameters #
#################################

# Extract Maximum reward
max_max = -Inf
max_mean = -Inf
max_min = - Inf
max_hp = (-Inf, -Inf, -Inf, -Inf)

for (hparams, rvec) in hyper_reward
    minval, maxval = extrema(rvec)
    if length(rvec) > 0 && maxval > max_max
        max_max = maxval
        max_mean = mean(rvec)
        max_min = minval
        max_hp = hparams
    end
end

# Print Best Point
println("Best Hyperparameter Cominbation -- Min: $(max_min) - Mean: $(max_mean) - Max: $(max_max)")
println("Exploration Constant: $(max_hp[1])")
println("Number of Iterations: $(max_hp[2])")
println("Solve Depth: $(max_hp[3])")
println("Discount Facgtor: $(max_hp[4])")

# Save Rest to File
open("hyperparam_search_rewards_$(Dates.format(now(), "yyyy-mm-ddTHH:MM:SS")).dat") do file
    column_string = join(["Reward $idx" for idx in 1:bin_iterations], ",")
    write(file, "Exploration Constant, Number Of Iterations, Search Depth, Discount Factor, $column_string\n")
    for hparams, rewards in hyper_reward
        ec, ni, sd, ds = hparams
        reward_string = join(rewards, ",")
        write(file, "$(ec),$ni,$sd,$ds,$reward_string\n")
    end
end

open("hyperparam_search_times_$(Dates.format(now(), "yyyy-mm-ddTHH:MM:SS")).dat") do file
    column_string = join(["Runtime $idx" for idx in 1:bin_iterations], ",")
    write(file, "Exploration Constant, Number Of Iterations, Search Depth, Discount Factor, $column_string\n")
    for hparams, times in hyper_time
        ec, ni, sd, ds = hparams
        time_string = join(times, ",")
        write(file, "$(ec),$ni,$sd,$ds,$time_string\n")
    end
end


text:
module TestLeafOptionalArgument

using Comonicon.AST
using Comonicon.JuliaExpr
using Comonicon.JuliaExpr: emit, emit_body, emit_norm_body, emit_dash_body
using Test

function foo(a, b = 2)
    @test a == 3
end

cmd = Entry(;
    version = v"1.2.0",
    root = LeafCommand(;
        fn = foo,
        name = "leaf",
        args = [
            Argument(; name = "a", type = Int),
            Argument(; name = "b", type = Int, require = false),
        ],
    ),
)

eval(emit(cmd))

@testset "test lead optional argument" begin
    @test command_main(["3"]) == 0
    @test command_main(["3", "5"]) == 0
end

end


text:
#
# TODO: add combined tests fn+loop, loop+fnCall and so on...
#
module TestCode
  using FactCheck
  import Code
  import Creature
  import Creature.f16
  import Helper
  import Mutator
  import Config
  include("Helper.jl")
  #
  # fn
  #
  facts("Testing if Code.fn() is added into empty script") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]

    code(d, :fn, org)
    @fact Helper.getArg(org.code, [2,2,1,1]) --> :func_2
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing two Code.fn()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]

    code(d, :fn, org)
    code(d, :fn, org)
    @fact length(org.funcs) --> 3
    @fact length(Helper.getLines(org.code, [2])) --> 4
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  #
  # fnCall
  #
  facts("Testing Code.fnCall(), but without a function") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :fnCall, org)

    @fact length(Helper.getLines(org.code, [2])) --> 2
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing Code.fnCall() for one Code.fn()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :fn, org)
    code(d, :fnCall, org)

    @fact length(Helper.getLines(org.code, [2])) --> 4
    @fact Helper.getHead(org.code, [2,2]) --> :function
    @fact Helper.getHead(org.code, [2,4]) --> :return
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing Code.fnCall() for one Code.fn() with body") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]

    code(d, :fn, org)
    code(d, :plus, org, [2,1,1])
    code(d, :fnCall, org)

    @fact length(Helper.getLines(org.code, [2])) --> 4
    @fact Helper.getHead(org.code, [2,2]) --> :function
    @fact Helper.getHead(org.code, [2,4]) --> :return
    @fact Code.eval(org.code)(d.cfg, org) --> val(org) + val(org)

    Manager.destroy(d.man)
  end
  facts("Testing Code.fnCall() after removing Code.fn()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]

    code(d, :fn, org)
    Mutator._onDel(d.cfg, org, Helper.CodePos(1,1,Creature.VAR_AMOUNT + 1), Code.CodePart(Code.fn, true))
    code(d, :fnCall, org)
    @fact length(Helper.getLines(org.code, [2])) --> Creature.VAR_AMOUNT + 1
    @fact length(org.funcs) --> 1
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)
    #
    # Adding another function and t's call
    #
    code(d, :fn, org)
    code(d, :fnCall, org)
    @fact length(org.funcs) --> 2
    Mutator._onDel(d.cfg, org, Helper.CodePos(1,1,Creature.VAR_AMOUNT + 2), Code.CodePart(Code.fn, true))
    code(d, :fnCall, org)
    @fact length(org.funcs) --> 2
    @fact length(Helper.getLines(org.code, [2])) --> Creature.VAR_AMOUNT + 3
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing two Code.fnCall() for one Code.fn()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]
    local index = Creature.VAR_AMOUNT + 3

    code(d, :fn, org)
    code(d, :fnCall, org)
    code(d, :fnCall, org)

    @fact Helper.getHead(org.code, [2,index]) --> :(=)
    @fact Helper.getArg(org.code, [2,index,2,1]) --> :func_2
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  #
  # if
  #
  facts("Testing Code.condition() with body") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]

    code(d, :condition, org)
    #
    # Updates compare operator to "===" to make condition true all the time
    #
    Helper.getArg(org.code, [2,2,1]).args[1] = :(===)
    Mutator._onAdd(d.cfg, org, Helper.CodePos(1,2,1), Code.CodePart(Code.pi, false))
    updateCode(org)
    @fact length(org.funcs[1].blocks) --> 2
    @fact Code.eval(org.code)(d.cfg, org) --> Float16(Base.pi)

    Manager.destroy(d.man)
  end
  facts("Testing Code.condition() inside function") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:codeFuncParamAmount=>1))
    local org = d.orgs[1]

    code(d, :fn, org)
    Mutator._onAdd(d.cfg, org, Helper.CodePos(2,1,2), Code.CodePart(Code.condition, true))
    @fact Helper.getHead(org.code, [2,2,2,2]) --> :if
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing Code.condition() inside other Code.condition()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]
    #
    # condition inside other function
    #
    code(d, :fn, org)
    @fact Mutator._onAdd(d.cfg, org, Helper.CodePos(2,1,2), Code.CodePart(Code.condition, true)) --> true
    @fact Mutator._onAdd(d.cfg, org, Helper.CodePos(2,2,1), Code.CodePart(Code.condition, true)) --> false
    #
    # conditions inside main function
    #
    code(d, :condition, org)
    @fact Mutator._onAdd(d.cfg, org, Helper.CodePos(1,2,1), Code.CodePart(Code.condition, true)) --> false
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  #
  # for
  #
  facts("Testing Code.loop() and empty code") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :loop, org)

    @fact length(Helper.getLines(org.code, [2])) --> 3
    @fact Helper.getHead(org.code, [2,2]) --> :for
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing Code.loop() with body") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :loop, org)
    Mutator._onAdd(d.cfg, org, Helper.CodePos(1,2,1), Code.CodePart(Code.pi, false))
    @fact Code.eval(org.code)(d.cfg, org) --> f16(pi)

    Manager.destroy(d.man)
  end
  facts("Testing Code.loop() near other Code.loop()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :loop, org)
    code(d, :loop, org)
    @fact Helper.getHead(org.code, [2,2]) --> :for
    @fact Helper.getHead(org.code, [2,3]) --> :for
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  facts("Testing Code.loop() inside other Code.loop()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :loop, org)
    @fact code(d, :loop, org, [1,2,2]) --> false
    @fact Code.eval(org.code)(d.cfg, org) --> val(org)

    Manager.destroy(d.man)
  end
  # #
  # # commands mix
  # #
  # facts("Testing Code.fn() inside Code.loop()") do
  #   local conf = Config.create()
  #   local org  = Creature.create(conf)
  #
  #   addVars(org, [2], Helper.CodePos(1,1,1))
  #   Mutator._onAdd(conf, org, Helper.CodePos(1,1,6), Code.CodePart(Code.loop, true))
  #   Mutator._onAdd(conf, org, Helper.CodePos(1,2,1), Code.CodePart(Code.fn, true))
  #
  #   @fact length(Helper.getLines(org.code, [2])) --> 3
  #   @fact length(Helper.getLines(org.code, [2,6,2])) --> 0
  #   @fact Code.eval(org.code)(conf, org) --> true
  # end
  # facts("Testing Code.loop() inside Code.fn()") do
  #   local conf = Config.create()
  #   local org  = Creature.create(conf)
  #
  #   conf.codeFuncParamAmount = 1
  #   Mutator._onAdd(conf, org, Helper.CodePos(1,1,1), Code.CodePart(Code.fn, true))
  #   addVar(org, [2,1,2], Helper.CodePos(2,1,1), Int8)
  #   Mutator._onAdd(conf, org, Helper.CodePos(2,1,2), Code.CodePart(Code.loop, true))
  #
  #   @fact length(Helper.getLines(org.code, [2])) --> 2
  #   @fact length(Helper.getLines(org.code, [2,1,2])) --> 3
  #   @fact Code.eval(org.code)(conf, org) --> true
  # end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
import LambdaTAME: write_smat
import MatrixNetworks: readSMAT
import SparseArrays: sprand
@testset "fileio tests" begin

    A = sprand(20,20,.2)
    #A = max.(A,A')
    @testset "test smat" begin

        mktempdir() do folder_path

            file = folder_path*"test_mat" #cleanup=true default should delete file
            @test_throws AssertionError write_smat(A,file;delimiter=" ") #no .smat name

            file *=".smat"
            write_smat(A,file;delimeter=' ')
            B = readSMAT(file)
            @test A == B
        end

    end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME