text:
module CartesianDiscreteModelsTests

using Gridap
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Geometry: CartesianMap
using Gridap.Io
using Test

domain = (0,1,0,1,0,1)
partition = (3,3,3)
parts = nothing
model = CartesianDiscreteModel(parts,domain,partition)
test_discrete_model(model)

domain = (0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)
test_discrete_model(model)
@test is_oriented(get_grid(model)) == true
@test is_oriented(get_grid_topology(model)) == true

grid_topology = get_grid_topology(model)
cell_to_pindices = get_cell_permutations(grid_topology)

labels = get_face_labeling(model)
@test num_tags(labels) == 28
@test num_entities(labels) == 27

map = get_cell_map(get_grid(model))
@test isa(map,CartesianMap)

desc = get_cartesian_descriptor(model)

function polar(q)
  r, t, z = q
  x = r*cos(t)
  y = r*sin(t)
  Point(x,y,z)
end

domain = (1,2,0,pi,0,0.5)
partition = (10,30,4)
model = CartesianDiscreteModel(domain,partition;map=polar)
test_discrete_model(model)
@test is_oriented(get_grid(model)) == true

model2 = from_dict(DiscreteModel,to_dict(model))
test_discrete_model(model2)

model3 = CartesianDiscreteModel(desc,CartesianIndex(2,2,2),CartesianIndex(3,3,3))
test_discrete_model(model3)

domain = (0,1,0,1)
partition = (4,4)
desc = CartesianDescriptor(domain,partition;isperiodic=(true,true))
model = CartesianDiscreteModel(desc)
@test num_vertices(model) == 16
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 0
model = CartesianDiscreteModel(desc,CartesianIndex(1,1),CartesianIndex(4,4))
@test num_vertices(model) == 16
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 0
model = CartesianDiscreteModel(desc,CartesianIndex(2,2),CartesianIndex(4,4))
@test num_vertices(model) == 16
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 7

desc = CartesianDescriptor(domain,partition;isperiodic=(true,false))
model = CartesianDiscreteModel(desc)
@test num_vertices(model) == 20
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 8
model = CartesianDiscreteModel(desc,CartesianIndex(1,1),CartesianIndex(4,4))
@test num_vertices(model) == 20
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 8
model = CartesianDiscreteModel(desc,CartesianIndex(2,2),CartesianIndex(4,4))
@test num_vertices(model) == 16
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 7

desc = CartesianDescriptor(domain,partition;isperiodic=(false,true))
model = CartesianDiscreteModel(desc)
@test num_vertices(model) == 20
model = CartesianDiscreteModel(desc,CartesianIndex(1,1),CartesianIndex(4,4))
@test num_vertices(model) == 20
model = CartesianDiscreteModel(desc,CartesianIndex(2,2),CartesianIndex(4,4))
@test num_vertices(model) == 16
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 7

desc = CartesianDescriptor(domain,partition;isperiodic=(false,false))
model = CartesianDiscreteModel(desc)
@test num_vertices(model) == 25
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 16
model = CartesianDiscreteModel(desc,CartesianIndex(1,1),CartesianIndex(4,4))
@test num_vertices(model) == 25
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 16
model = CartesianDiscreteModel(desc,CartesianIndex(2,2),CartesianIndex(4,4))
@test num_vertices(model) == 16
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 7

desc = CartesianDescriptor(domain,partition)
remove_boundary = (false,true)
model = CartesianDiscreteModel(desc,CartesianIndex(1,1),CartesianIndex(4,4),remove_boundary)
@test num_vertices(model) == 25
@test count(get_face_mask(get_face_labeling(model),"boundary",0)) == 10

#using Gridap.Visualization
#writevtk(model2,"model2")

end # module


text:
import Test: @test, @testset
import Luna: Ionisation, Maths, PhysData, Tools
import NumericalIntegration: integrate, SimpsonEven

@test Ionisation.ionrate_ADK(:He, 1e10) ≈ 1.2416371415312408e-18
@test Ionisation.ionrate_ADK(:He, 2e10) ≈ 1.0772390893742478
@test Ionisation.ionrate_ADK(:HeJ, 2e10) ≈ 1.0772390893742478
@test Ionisation.ionrate_ADK(:Ar , 7e9) ≈ 2.4422306306649472e-08
@test Ionisation.ionrate_ADK(:Ar , 8e9) ≈ 4.494711488416766e-05

E = collect(range(1e9, 1e11; length=32))
@test Ionisation.ionrate_ADK(:He, E) == Ionisation.ionrate_ADK(:He, -E)

@test isapprox(Ionisation.ionrate_PPT(:He, 800e-9, 1e10), 1.4130113877738475e-5, rtol=1e-3)
@test isapprox(Ionisation.ionrate_PPT(:He, 800e-9, 1.3e10), 0.04585332982943, rtol=1e-3)

Emin = 1e9
Emax = 1e11
N = 2^10
E = collect(range(Emin, stop=Emax, length=N))
rate = Ionisation.ionrate_PPT.(:He, 800e-9, E)
ifun(E0) =  E0 <= Emin ? 2 :
            E0 >= Emax ? N :
            ceil(Int, (E0-Emin)/(Emax-Emin)*N) + 1
ifun2(x0) = x0 <= E[1] ? 2 :
                x0 >= E[end] ? length(E) :
                findfirst(x -> x>x0, E)
spl1 = Maths.CSpline(E, rate, ifun)
spl2 = Maths.CSpline(E, rate, ifun2)
idx1 = ifun2.(E) # calculated indices
idx2 = ifun.(E) # indices found with brute-force method
@test all(idx1 .== idx2)
@test all(spl1.(E) .== spl2.(E))

ratefun! = Ionisation.ionrate_fun!_PPTaccel(:He, 800e-9)
out = similar(E)
ratefun!(out, E)
@test all(isapprox.(out, rate, rtol=1e-2))

outneg = similar(out)
ratefun!(outneg, -E)
@test out == outneg

outneg = similar(out)
ratefun!(outneg, -E)
@test out == outneg

##
@testset "cycle-averaging $gas" for gas in (:He, :Ar, :Xe)

    bsi = Tools.field_to_intensity(
        Ionisation.barrier_suppression(
            PhysData.ionisation_potential(gas),
            1))
    isy = 10 .^ collect(16:0.01:log10(bsi))
    E0 = Tools.intensity_to_field.(isy)

    Ip_au = PhysData.ionisation_potential(gas; unit=:atomic)
    F0_au = (2Ip_au)^(3/2)
    F0 = F0_au*PhysData.au_Efield

    t = collect(range(0, 2π; length=2^12))
    Et = sin.(t)

    adk_avg = zero(E0)
    rf! = Ionisation.ionrate_fun!_ADK(gas)
    out = zero(Et)
    for (idx, E0i) in enumerate(E0)
        rf!(out, E0i*Et)
        adk_avg[idx] = 1/2π * integrate(t, out, SimpsonEven())
    end

    adk = Ionisation.ionrate_ADK.(gas, E0)
    adk_avg_kw = Ionisation.ionrate_ADK.(gas, E0; cycle_average=true)
    @test all(isapprox.(adk_avg, adk_avg_kw; rtol=0.05))
end
##
# this_folder = dirname(@__FILE__)
# import CSV
# import PyPlot: plt, pygui
# pygui(true)
# #### Ilkov et al
# k = collect(range(8, stop=12, length=200))
# E = 10 .^ k
# ppt = Ionisation.ionrate_PPT(:He, 10.6e-6, E)
# ppt_cycavg = Ionisation.ionrate_PPT(:He, 10.6e-6, E, rcycle=false)
# adk = Ionisation.ionrate_ADK(:He, E)

# dat = CSV.read(joinpath(this_folder, "Ilkov_PPT_He.csv"))
# dat = convert(Matrix, dat) # ionrate [1/s] vs field [V/cm]

# plt.figure()
# plt.loglog(E, ppt, label="PPT")
# plt.loglog(E, ppt_cycavg, label="PPT cycle averaged")
# plt.loglog(E, adk, label="ADK")
# plt.loglog(dat[:, 1].*100, dat[:, 2], label="Ilkov et al. PPT")
# plt.xlim(1e9, 1e12)
# plt.ylim(1, 1e18)
# plt.legend()

# ### Chang, Fundamentals of Attosecond Optics
# dat = CSV.read(joinpath(this_folder, "Chang_PPT.csv"))
# dat = convert(Matrix, dat) # ionrate [1/fs] vs intensity [1e14 W/cm^2]
# intensity = range(0.1, 25; length=1000) * 1e18 # W/m^2
# E = Tools.intensity_to_field.(intensity)
# ppt = Ionisation.ionrate_PPT(:He, 390e-9, E)
# ppt_cycavg = Ionisation.ionrate_PPT(:He, 390e-9, E, rcycle=false)
# adk = Ionisation.ionrate_ADK(:He, E)

# s = sortperm(dat[:, 1])

# plt.figure()
# plt.loglog(intensity*1e-18, ppt*1e-15, label="PPT")
# plt.loglog(intensity*1e-18, adk*1e-15, label="ADK")
# plt.loglog(intensity*1e-18, ppt_cycavg*1e-15, label="PPT cycle averaged")
# plt.loglog(dat[s, 1], dat[s, 2], label="Chang PPT")
# plt.ylim(1e-14, 1)
# plt.xlim(extrema(intensity.*1e-18))
# plt.legend()
# plt.xlabel("Intensity (10\$^{14}\$ W/cm\$^2\$)")
# plt.ylabel("Ionisation rate (1/fs)")

### Couairon
# Ip = 12.063 * PhysData.electron
# dat = CSV.read(joinpath(this_folder, "Couairon_PPT.csv"))
# dat = convert(Matrix, dat) # ionrate [1/s] vs intensity [W/cm^2]
# k = collect(range(12, 15, length=500))
# intensity = 10 .^ k .* 1e4 # W/m^2
# E = Tools.intensity_to_field.(intensity)
# ppt = Ionisation.ionrate_PPT(Ip, 800e-9, 1, 1, E)
# ppt_cycavg = Ionisation.ionrate_PPT(Ip, 800e-9, 1, 1, E; rcycle=false)
# adk = Ionisation.ionrate_ADK(Ip, E)

# s = sortperm(dat[:, 1])

# plt.figure()
# plt.loglog(intensity*1e-4, ppt, label="PPT")
# plt.loglog(intensity*1e-4, adk, label="ADK")
# plt.loglog(intensity*1e-4, ppt_cycavg, label="PPT cycle averaged")
# plt.loglog(dat[s, 1], dat[s, 2], label="Couairon PPT")
# plt.ylim(1, 1e17)
# plt.xlim(extrema(intensity.*1e-4))
# plt.legend()
# plt.xlabel("Intensity (W/cm\$^2\$)")
# plt.ylabel("Ionisation rate (1/s)")
# plt.title("O\$_2\$ ionisation at 800 nm")

ext: jl


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

module TestHessenberg

using Test, LinearAlgebra, Random

let n = 10
    srand(1234321)

    Areal  = randn(n,n)/2
    Aimg   = randn(n,n)/2

    @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        A = eltya == Int ?
                rand(1:7, n, n) :
                convert(Matrix{eltya}, eltya <: Complex ?
                    complex.(Areal, Aimg) :
                    Areal)

        if eltya != BigFloat
            H = hessenberg(A)
            @test size(H.Q, 1) == size(A, 1)
            @test size(H.Q, 2) == size(A, 2)
            @test size(H.Q) == size(A)
            @test_throws ErrorException H.Z
            @test convert(Array, H) ≈ A
            @test (H.Q * H.H) * H.Q' ≈ A
            @test (H.Q' *A) * H.Q ≈ H.H
            #getindex for HessenbergQ
            @test H.Q[1,1] ≈ Array(H.Q)[1,1]
        end
    end
end

end # module TestHessenberg

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
##
using Revise
using GLMakie
using Pkg
Pkg.activate(".")
using ProteinStructureReconstruction
using HDF5
using Images
import Statistics: median
import GeometryBasics: Point2f0
##
filepath = "./Data/Fed_X63_Z3_SIM.h5"
h5file = h5open(filepath, "r")
##
# TODO: don't load everything in memory at once -  just draw overlay first, then
##
x_range, y_range = 2000:2300, 1600:1900;
img_data = h5file["Data"][x_range, y_range, :, :] |> x -> reinterpret(N0f16, x);


text:
#=
    envelope
    Copyright © 2019 Mark Wells <mwellsa@gmail.com>

    Distributed under terms of the AGPL-3.0 license.
=#

""" function get_support(x)
input:
    x -> vector of knots
return:
    support -> (minimum(x), maximum(x) - minimum(x))
"""
@inline function get_support(knots ::AbstractVector{Float64}) ::NTuple{2,Float64}
    return (knots[1], knots[end]-knots[1])
end

""" function get_support(X)
input:
    knots -> tuple of knots (xknots, yknots, ...)
return:
    retval -> tuple of supports ((xmin,xspan), (ymin,yspan), ...)
"""
@inline function get_support(knots ::NTuple{D,AbstractVector}) where D
    return ntuple(i -> get_support(knots[i]), Val(D))
end

""" check_neighbors(x)
input:
    x -> boolean vector
output:
    boolean vector

returns true at indices where x[i-1] & x[i] & x[i+1] are true
"""
@inline function check_neighbors(x ::AbstractVector{Bool})
    retval = BitArray(x)
    if length(retval) > 1
        retval[1] &= x[2]
        retval[end] &= x[end-1]
        retval[2:end-1] .&= x[1:end-2] .& x[3:end]
    end
    return retval
end

""" check_neighbors(X)
input:
    X -> Boolean Array
output:
    Boolean Array

N-dimensional version of check_neighbors(x)
"""
@inline function check_neighbors(X ::AbstractArray{Bool,D}) where D
    retval = BitArray(X)
    s = size(X)
    for i in 1:D
        for j in 1:s[i]
            _retval = selectdim(retval, i, j)
            _X = selectdim(X, i, j)
            _retval .&= check_neighbors(_X)
        end
    end
    return retval
end

""" get_support(knots, probs)
input:
    knots -> vectors specifing the knots
    probs -> vectors prob
output:
    min/span values

This routine accounts for range of zero prob and limits the support to exclude these regions.
Currently only outer zero regions can be excluded.
"""
@inline function get_support(knots ::AbstractVector,
                             probs ::AbstractVector)
    zmask = .!check_neighbors(probs .≤ 0)
    inds = findall(zmask)
    min_ind = minimum(inds)
    max_ind = maximum(inds)
    return get_support(knots[min_ind:max_ind])
end

""" get_support(knots, probs)
input:
    knots -> tuple of vectors specifing the knots
    probs -> array of prob
output:
    tuple of min/span values

This routine accounts for boxes of zero prob and limits the support to exclude these regions
Currently only outer zero regions can be excluded.
"""
@inline function get_support(knots ::NTuple{D,AbstractVector},
                     probs ::AbstractArray{Float64,D}) where D
    zmask = .!check_neighbors(iszero.(probs))
    cinds = findall(zmask)
    min_cind = minimum(cinds)
    max_cind = maximum(cinds)
    return ntuple(i -> get_support(knots[i][min_cind[i]:max_cind[i]]), Val(D))
end

""" get_bounding_box(pnt, knots)
input:
    pnt -> missing/float64 to surround if missing it is unbounded
    knots -> vector of knots
output:
    indices to access the the bounded range
"""
@inline function get_bounding_box(pnt ::Union{Missing,Float64},
                                  knots ::AbstractVector{Float64})
    left = Int(1)
    right = Int(length(knots))
    if !ismissing(pnt)
        left = Int(max(left, searchsortedlast(knots, pnt)))
        right = Int(min(right, searchsortedfirst(knots, pnt)))
    end
    return UnitRange(left,right)
end

""" get_bounding_box(pnt, knots)
inputs:
    pnt -> tuple of values (or missing) specify the point to bound
    knots -> tuple of knot vectors
output:
    tuple of indices that access the bounded box
"""
@inline function get_bounding_box(pnt ::NTuple{D,Union{Missing,Float64}},
                          knots ::NTuple{D,AbstractVector{Float64}}) where D
    return ntuple(i -> get_bounding_box(pnt[i], knots[i]), Val(D))
end

@inline function apply_bbox(x ::AbstractArray{T,N},
                            i ::NTuple{N,UnitRange{Int}}) where {T,N}
    return x[i...]
end

""" function get_support(bounding_box, knots, probs)
inputs:
    bounding_box -> indices of knots that specify the box
    knots -> tuple of knot vectors
    probs -> array of prob
output:
    support -> tuple of min/span
"""
@inline function get_support(bounding_box ::NTuple{D,UnitRange{Int}},
                             knots ::NTuple{D,AbstractVector},
                             probs ::AbstractArray{Float64,D}) where D
    _knots = ntuple(i -> knots[i][bounding_box[i]], Val(D))
    _probs = apply_bbox(probs, bounding_box)
    return get_support(_knots, _probs)
end

""" function get_envelope(pnt, knots, probs)
inputs:
    pnt -> tuple of floats/missing
    knots -> knots describing the grid of probs
    probs -> grid of probs
return: (support, pmax)
    support -> tuple of (xmin,xspan, ymin,yspan)
    pmax -> maximum prob accessible to pnt
"""
@inline function get_envelope(pnt ::NTuple{D,Union{Missing,Float64}},
                              knots ::NTuple{D,AbstractVector},
                              probs ::AbstractArray{Float64,D}) where D
    bounding_box = get_bounding_box(pnt, knots)
    support = get_support(bounding_box, knots, probs)
    pmax = maximum(apply_bbox(probs, bounding_box))
    return (support, pmax)
end

@inline function get_envelope(slices ::AbstractMatrix{Union{Missing,Float64}},
                              knots ::NTuple{D,AbstractVector},
                              probs ::AbstractArray{Float64,D}) where D
    n = last(size(slices))
    retval = [get_envelope(ntuple(i -> slices[i,j], Val(D)), knots, probs) for j = 1:n]
    return first.(retval), last.(retval)
end


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

name = "HHsuite"
version = v"3.3.0"

# TODO
# - MPI support

sources = [
    ArchiveSource("https://github.com/soedinglab/hh-suite/archive/refs/tags/v$(version).tar.gz",
                  "dd67f7f3bf601e48c9c0bc4cf1fbe3b946f787a808bde765e9436a48d27b0964"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/hh-suite-*/

# patch CMakeLists.txt so it doesn't set -march unnecessarily on ARM
atomic_patch -p1 ../patches/arm-simd-march-cmakefile.patch

# macos: use gcc/g++ so we can use openmp
if [[ "${target}" == *-apple-darwin* ]]; then
    CMAKE_TARGET_TOOLCHAIN="${CMAKE_TARGET_TOOLCHAIN/%.cmake/_gcc.cmake}"
    echo "[INFO] setting CMAKE_TARGET_TOOLCHAIN = ${CMAKE_TARGET_TOOLCHAIN}"
fi

export CPPFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"

arch_opts=
if [[ ${target} == x86_64-* ]]; then
    arch_opts="-DHAVE_SSE2=ON -DHAVE_SSE4_1=ON -DHAVE_AVX2=ON"
elif [[ ${target} == aarch64-* ]]; then
    arch_opts="-DHAVE_ARM8=ON"
elif [[ ${target} == powerpc64le-* ]]; then
    arch_opts="-DHAVE_POWER8=ON -DHAVE_POWER9=ON"
fi

mkdir build && cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DNATIVE_ARCH=OFF \
    ${arch_opts}
make -j${nproc}
make install

install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# Build failures
# - x86_64-w64-mingw32
#   fatal error: err.h: No such file or directory
# - powerpc64le-linux-gnu
#   build tries to use x86 vector instructions
#
platforms = supported_platforms(; exclude = p -> Sys.iswindows(p) || arch(p) == "powerpc64le")
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("a3m_database_extract", :a3m_database_extract),
    ExecutableProduct("a3m_database_filter", :a3m_database_filter),
    ExecutableProduct("a3m_database_reduce", :a3m_database_reduce),
    ExecutableProduct("a3m_extract", :a3m_extract),
    ExecutableProduct("a3m_reduce", :a3m_reduce),
    ExecutableProduct("cstranslate", :cstranslate),
    ExecutableProduct("ffindex_apply", :ffindex_apply),
    ExecutableProduct("ffindex_build", :ffindex_build),
    ExecutableProduct("ffindex_from_fasta", :ffindex_from_fasta),
    ExecutableProduct("ffindex_from_fasta_with_split", :ffindex_from_fasta_with_split),
    ExecutableProduct("ffindex_get", :ffindex_get),
    ExecutableProduct("ffindex_modify", :ffindex_modify),
    ExecutableProduct("ffindex_order", :ffindex_order),
    ExecutableProduct("ffindex_reduce", :ffindex_reduce),
    ExecutableProduct("ffindex_unpack", :ffindex_unpack),
    ExecutableProduct("hhalign", :hhalign),
    ExecutableProduct("hhalign_omp", :hhalign_omp),
    ExecutableProduct("hhblits", :hhblits),
    ExecutableProduct("hhblits_ca3m", :hhblits_ca3m),
    ExecutableProduct("hhblits_omp", :hhblits_omp),
    ExecutableProduct("hhconsensus", :hhconsensus),
    ExecutableProduct("hhfilter", :hhfilter),
    ExecutableProduct("hhmake", :hhmake),
    ExecutableProduct("hhsearch", :hhsearch),
    ExecutableProduct("hhsearch_omp", :hhsearch_omp),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll")),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using StatsBase, Random

function trunc_powerlaw_weigths(α::Real, v_min::Int, v_max::Int)
    @assert α >= 1
    @assert 1 <= v_min <= v_max
    w = Weights([1/i^α for i in v_min:v_max])
end

function sample_trunc_powerlaw(α::Real, v_min::Int, v_max::Int, n::Int)
    @assert α >= 1
    @assert 1 <= v_min <= v_max
    @assert n > 0
    w = Weights([1/i^α for i in v_min:v_max])
    sample(v_min:v_max, w, n)
end

sample_trunc_powerlaw(W::Weights, v_min::Int, v_max::Int, n::Int) = (@assert n > 0; sample(v_min:v_max, W, n))

"""
    get_ev(α, v_min, v_max)

Return the expected value of truncated discrete power law distribution
with truncation range `[v_min, v_max]` and exponent `α`.
"""
function get_ev(α, v_min, v_max)
    @assert α > 1
    @assert 1 <= v_min <= v_max
    w = [1/i^α for i in v_min:v_max]
    w /= sum(w)
    sum(x*w[i] for (i, x) in enumerate(v_min:v_max))
end

"""
    find_v_min(α, v_max, v_ev)

Return the lower truncation given expected value `v_ev` and upper truncation `v_max`
of truncated discrete power law distribution with exponent `α`.
"""
function find_v_min(α, v_max, v_ev)
    @assert v_ev <= v_max
    v_min = v_ev
    cur_err = abs(v_ev, get_ev(α, v_min, v_max))
    while true
        v_min == 1 && return 1
        v_min -= 1
        next_err = abs(v_ev, get_ev(α, v_min, v_max))
        if next_err > cur_err
            return v_min + 1
        else
            cur_err = next_err
        end
    end
end

"""
    sample_degrees(τ₁, d_min, d_max, n, max_iter)

Return a vector of length `n` of sampled degrees of vertices following a truncated
discrete power law distribution with truncation range `[d_min, d_max]` and exponent `τ₁`.

The sampling is attempted `max_iter` times to find an admissible result.
In case of failure a correction of degree is applied to a degree of the vertex
with highest degree.

The producedure does not check if the returned vector is a graphical degree sequence.
"""
function sample_degrees(τ₁, d_min, d_max, n, max_iter)
    local s
    for i in 1:max_iter
        s = sample_trunc_powerlaw(τ₁, d_min, d_max, n)
        iseven(sum(s)) && return sort!(s, rev=true)
    end
    @warn "Failed to sample an admissible degree sequence in $max_iter draws. Fixing"
    i = argmax(s)
    if s[i] == 0
        s[i] = 1 # this should not happen in practice
    else
        s[i] -= 1
    end
    return sort!(s, rev=true)
end

"""
    sample_communities(τ₂, c_min, c_max, n, max_iter)

Return a vector of sampled community sizes following a truncated
discrete power law distribution with truncation range `[c_min, c_max]` and exponent `τ₂`.
The sum of sizes is equal to number of vertices in the graph `n`.

The sampling is attempted `max_iter` times to find an admissible result.
In case of failure a correction of community sizes is applied to the sampled sequence
that was closest to a feasible one to ensure that the result is admissible.
"""
function sample_communities(τ₂, c_min, c_max, n, max_iter)
    @assert 1 <= c_min <= c_max
    l_min = n / c_max
    l_max = n / c_min
    @assert l_min >= 1
    @assert ceil(l_min) <= floor(l_max)
    local best_s
    local best_ss = typemax(Int)
    w = trunc_powerlaw_weigths(τ₂, c_min, c_max)
    for i in 1:max_iter
        s = sample_trunc_powerlaw(w, c_min, c_max, ceil(Int, l_max))
        stopidx = 0
        ss = 0
        while ss < n
            stopidx += 1
            ss += s[stopidx]
        end
        ss == n && return sort!(s[1:stopidx], rev=true)
        if ss < best_ss
            best_ss = ss
            best_s = s[1:stopidx]
        end
    end
    @warn "failed to sample an admissible community sequence in $max_iter draws. Fixing"
    shuffle!(best_s)
    if length(best_s) > l_max
        best_s = best_s[1:l_max]
        best_ss = sum(best_s)
    end
    i = 0
    while best_ss != n
        if i ≥ length(best_s)
            i = 0
            shuffle!(best_s)
        end
        i += 1
        change = sign(n - best_ss)
        if change > 0
            best_s[i] < c_max || continue
        else
            best_s[i] > c_min || continue
        end
        best_ss += change
        best_s[i] += change
    end
    return sort!(best_s, rev=true)
end


text:
function preprocess_aaindex()
    data_dir = datadep"AAindex"
    aaindex_files = readdir(data_dir; join = true)

    # this might be some no-no stuff
    jldopen(joinpath(data_dir, "aaindex.jld2"), "w") do file
        ids, descriptions, types = [], [], []

        entries = JLD2.Group(file, "entries")

        foreach(aaindex_files) do aaindex
            foreachentry(aaindex) do entry
                record = parse(entry)
                (; id, description) = record.metadata

                push!(ids, id)
                push!(descriptions, description)
                push!(types, typeof(record))

                entries[id] = record
            end
        end

        file["ids"] = ids
        file["descriptions"] = descriptions
        file["types"] = types
    end
end


text:
module DICOM

export dcm_parse, dcmdir_parse, dcm_write, lookup, lookup_vr, rescale!
export @tag_str

include("DICOMData.jl")

"""
    @tag_str(s)

Return the dicom tag, corresponding to the string `s`.
```jldoctest
julia> using DICOM

julia> tag"ROI Mean"
(0x6000, 0x1302)
```
"""
macro tag_str(s)
    key = Symbol(filter(x -> !isspace(x), s))
    DICOM.fieldname_dict[key]
end

# Create dicom dictionary - used for reading/writing DICOM files
# Keys are tuple containing hex Group and Element of DICOM entry
# Stored value is String array: [Field Name, VR, Data Length], e.g:
# Julia> DICOM.dcm_dict[(0x0008,0x0005)]
# 3-element Array{String,1}:
#  "Specific Character Set"
#  "CS"
#  "1-n"
include("dcm_dict.jl")  # const dcm_dict = ...

# For convenience, dictionary to get hex tag from field name, e.g:
# Julia> DICOM.fieldname_dict["Specific Character Set"]
# (0x0008, 0x0005)
fieldname_dict = Dict(val[1] => key for (key, val) in dcm_dict)

# These "empty" values are used internally. They are returned if a search fails.
const empty_vr = "" # Can be any VR that doesn't exist
const empty_vr_lookup = ["", empty_vr, ""] # Used in lookup_vr as failure state
const empty_tag = (0x0000, 0x0000)
const empty_dcm_dict = Dict(DICOM.empty_tag => nothing)

const VR_names = [
    "AE",
    "AS",
    "AT",
    "CS",
    "DA",
    "DS",
    "DT",
    "FL",
    "FD",
    "IS",
    "LO",
    "LT",
    "OB",
    "OF",
    "OW",
    "PN",
    "SH",
    "SL",
    "SQ",
    "SS",
    "ST",
    "TM",
    "UI",
    "UL",
    "UN",
    "US",
    "UT",
]

# mapping UID => bigendian? explicitvr?
const meta_uids = Dict([
    ("1.2.840.10008.1.2", (false, false)),
    ("1.2.840.10008.1.2.1", (false, true)),
    ("1.2.840.10008.1.2.1.99", (false, true)),
    ("1.2.840.10008.1.2.2", (true, true)),
])

"""
    lookup_vr(tag::Tuple{Integer,Integer})

Return VR value for tag from DICOM dictionary

# Example
```jldoctest
julia> lookup_vr((0x0028,0x0004))
"CS"
```
"""
function lookup_vr(gelt::Tuple{UInt16,UInt16})
    if gelt[1] & 0xff00 == 0x5000
        gelt = (0x5000, gelt[2])
    elseif gelt[1] & 0xff00 == 0x6000
        gelt = (0x6000, gelt[2])
    end
    r = get(dcm_dict, gelt, empty_vr_lookup)
    return (r[2])
end

function rescale!(dcm::DICOMData, direction = :forward)
    if !haskey(dcm, tag"Rescale Intercept") || !haskey(dcm, tag"Rescale Slope")
        return dcm
    end
    if direction == :forward
        dcm[tag"Pixel Data"] =
            @. dcm[tag"Pixel Data"] * dcm[tag"Rescale Slope"] + dcm[tag"Rescale Intercept"]
    else
        pixel_data = dcm[tag"Pixel Data"]
        @. pixel_data -= dcm[tag"Rescale Intercept"]
        @. pixel_data /= dcm[tag"Rescale Slope"]
        dtype = determine_dtype(dcm)
        dcm[tag"Pixel Data"] = @. convert(dtype, round(pixel_data))
    end
    return dcm
end

if ENDIAN_BOM == 0x04030201
    order(x, endian) = endian == :little ? x : bswap.(x)
else
    order(x, endian) = endian == :big ? x : bswap.(x)
end

always_implicit(grp, elt) =
    (grp == 0xFFFE && (elt == 0xE0DD || elt == 0xE000 || elt == 0xE00D))

function find_dicom_files(dir)
    files = joinpath.(dir, readdir(dir))
    dicom_files = filter(file -> isdicom(file), files)
    return dicom_files
end

function isdicom(file)
    bytes = read(file, 132)[end-3:end]
    String(bytes) == "DICM"
end

function dcmdir_parse(dir; kwargs...)
    dicom_files = find_dicom_files(dir)
    unsorted_dicoms = [dcm_parse(file; kwargs...) for file in dicom_files]
    dicoms = sort!(unsorted_dicoms, by = dicom -> dicom[tag"Instance Number"])
    return dicoms
end

"""
   dcm_parse(fn::AbstractString)

Reads file fn and returns a Dict
"""
function dcm_parse(fn::AbstractString; kwargs...)
    st = open(fn)
    dcm = dcm_parse(st; kwargs...)
    close(st)
    dcm
end

"""
   dcm_parse(st::IO)

Reads IO st and returns a Dict
"""
function dcm_parse(
    st::IO;
    preamble = true,
    max_group = 0xffff,
    aux_vr = Dict{Tuple{UInt16,UInt16},String}(),
)
    if preamble
        check_preamble(st)
    end
    parsed = read_meta(st)
    is_explicit, endian = determine_explicitness_and_endianness(parsed)
    dcm = DICOMData(parsed.meta, endian, is_explicit, aux_vr)
    read_body!(st, dcm; max_group = max_group)
    # dcm = DICOMData(dcm, endian, is_explicit, vr, aux_vr)
    return dcm
end

function check_preamble(st)
    # First 128 can be skipped
    skip(st, 128)
    # "DICM" identifier must be after the first 128 bytes
    sig = String(read!(st, Array{UInt8}(undef, 4)))
    if sig != "DICM"
        error("dicom: invalid file header")
    end
    return
end

# Meta is always explicit VR / little endian
function read_meta(st::IO)
    dcm = DICOMData(Dict{Tuple{UInt16,UInt16},Any}(), :little, true, Dict{Tuple{UInt16,UInt16},String}())
    while true
        pos = position(st)
        (gelt, data, vr) = read_element(st, dcm)
        grp = gelt[1]
        if grp > 0x0002 || gelt == empty_tag
            seek(st, pos)
            return dcm
        else
            dcm[gelt] = data
        end
    end
    error("Unexpected break from loop while reading preamble")
end

function determine_explicitness_and_endianness(dcm)
    # Default is implicit_vr & little-endian
    if !haskey(dcm, (0x0002, 0x0010))
        return (false, :little)
    end
    metaInfo = get(meta_uids, dcm[(0x0002, 0x0010)], (false, true))
    explicitness = metaInfo[2]
    if metaInfo[1]
        endianness = :big
    else
        endianness = :little
    end
    return explicitness, endianness
end

function read_body!(st, dcm; max_group)
    vrs = dcm.vr
    while true
        (gelt, data, vr) = read_element(st, dcm)
        if gelt == empty_tag || gelt[1] > max_group
            break
        else
            dcm[gelt] = data
            vrs[gelt] = vr
        end
    end
    return dcm
end

function read_element(st::IO, dcm::DICOMData)
    is_explicit = dcm.isexplicit; endian = dcm.endian; aux_vr = dcm.vr;
    local grp
    try
        grp = read_group_tag(st, endian)
    catch
        return (empty_tag, 0, empty_vr)
    end
    elt = read_element_tag(st, endian)
    gelt = (grp, elt)
    vr, lentype = determine_vr_and_lentype(st, gelt, is_explicit, aux_vr)
    sz = read_element_size(st, lentype, endian)
    # Empty VR can be supplied in aux_vr to skip an element
    if isempty(vr)
        sz = isodd(sz) ? sz + 1 : sz
        skip(st, sz)
        return (read_element(st::IO, dcm))
    end

    data = vr == "ST" || vr == "LT" || vr == "UT" || vr == "AS" ?
        String(read!(st, Array{UInt8}(undef, sz))) :
        sz == 0 || vr == "XX" ? Any[] :
        vr == "SQ" ? sequence_parse(st, sz, dcm) :
        gelt == (0x7FE0, 0x0010) ? pixeldata_parse(st, sz, vr, dcm) :
        sz == 0xffffffff ? undefined_length(st, vr) :
        vr == "FL" ? numeric_parse(st, Float32, sz, endian) :
        vr == "FD" ? numeric_parse(st, Float64, sz, endian) :
        vr == "SL" ? numeric_parse(st, Int32, sz, endian) :
        vr == "SS" ? numeric_parse(st, Int16, sz, endian) :
        vr == "UL" ? numeric_parse(st, UInt32, sz, endian) :
        vr == "US" ? numeric_parse(st, UInt16, sz, endian) :
        vr == "OB" ? order(read!(st, Array{UInt8}(undef, sz)), endian) :
        vr == "OF" ? order(read!(st, Array{Float32}(undef, div(sz, 4))), endian) :
        vr == "OW" ? order(read!(st, Array{UInt16}(undef, div(sz, 2))), endian) :
        vr == "AT" ?
        [order(read!(st, Array{UInt16}(undef, 2)), endian) for n = 1:div(sz, 4)] :
        vr == "DS" ?
        map(x -> x == "" ? 0.0 : parse(Float64, x), string_parse(st, sz, 16, false)) :
        vr == "IS" ?
        map(x -> x == "" ? 0 : parse(Int, x), string_parse(st, sz, 12, false)) :
        vr == "AE" ? string_parse(st, sz, 16, false) :
        vr == "CS" ? string_parse(st, sz, 16, false) :
        vr == "SH" ? string_parse(st, sz, 16, false) :
        vr == "LO" ? string_parse(st, sz, 64, false) :
        vr == "UI" ? string_parse(st, sz, 64, false) :
        vr == "PN" ? string_parse(st, sz, 64, true) :
        vr == "DA" ? string_parse(st, sz, 10, true) :
        vr == "DT" ? string_parse(st, sz, 26, false) :
        vr == "TM" ? string_parse(st, sz, 16, false) :
        order(read!(st, Array{UInt8}(undef, sz)), endian)

    if isodd(sz) && sz != 0xffffffff
        skip(st, 1)
    end

    # For convenience, get rid of array if it is just acting as a container
    # Exception is "SQ", where array is part of structure
    if length(data) == 1 && vr != "SQ"
        data = data[1]
        # Sometimes it is necessary to go one level deeper
        if length(data) == 1
            data = data[1]
        end
    end

    return (gelt, data, vr)
end

read_group_tag(st, endian) = order(read(st, UInt16), endian)
read_element_tag = read_group_tag
read_element_size(st, lentype, endian) = order(read(st, lentype), endian)

function determine_vr_and_lentype(st, gelt, is_explicit, aux_vr)
    (grp, elt) = gelt
    lentype = UInt32
    if is_explicit && !always_implicit(grp, elt)
        vr = String(read!(st, Array{UInt8}(undef, 2)))
        if vr in ("OB", "OW", "OF", "SQ", "UT", "UN")
            skip(st, 2)
        else
            lentype = UInt16
        end
        diffvr = !isequal(vr, lookup_vr(gelt))
    else
        vr = elt == 0x0000 ? "UL" : lookup_vr(gelt)
    end
    if isodd(grp) && grp > 0x0008 && 0x0010 <= elt < +0x00FF
        # Private creator
        vr = "LO"
    elseif isodd(grp) && grp > 0x0008
        # Assume private
        vr = "UN"
    end
    if haskey(aux_vr, gelt)
        vr = aux_vr[gelt]
    end
    if vr === empty_vr
        if haskey(aux_vr, (0x0000, 0x0000))
            vr = aux_vr[(0x0000, 0x0000)]
        elseif !haskey(aux_vr, gelt)
            error("dicom: unknown tag ", gelt)
        end
    end
    return vr, lentype
end

numeric_parse(st::IO, T::DataType, sz, endian) =
    order(T[read(st, T) for i = 1:div(sz, sizeof(T))], endian)

function string_parse(st, sz, maxlen, spaces)
    endpos = position(st) + sz
    data = [""]
    first = true
    while position(st) < endpos
        c = !first || spaces ? read(st, Char) : skip_spaces(st, endpos)
        if c == '\\'
            push!(data, "")
            first = true
        elseif c == '\0'
            break
        else
            data[end] = string(data[end], c)  # TODO: inefficient
            first = false
        end
    end
    if !spaces
        return map(rstrip, data)
    end
    return data
end

function skip_spaces(st, endpos)
    while true
        c = read(st, Char)
        if c != ' ' || position(st) == endpos
            return c
        end
    end
end

function sequence_parse(st, sz, dcm)
    is_explicit = dcm.isexplicit; endian = dcm.endian; aux_vr = dcm.vr
    sq = Array{Dict{Tuple{UInt16,UInt16},Any},1}()
    while sz > 0
        grp = read_group_tag(st, endian)
        elt = read_element_tag(st, endian)
        itemlen = read_element_size(st, UInt32, endian)
        if grp == 0xFFFE && elt == 0xE0DD
            return sq
        end
        if grp != 0xFFFE || elt != 0xE000
            error("dicom: expected item tag in sequence")
        end
        push!(sq, sequence_item(st, itemlen, dcm))
        sz -= 8 + (itemlen != 0xffffffff) * itemlen
    end
    return sq
end

function sequence_item(st::IO, sz, dcm)
    item = Dict{Tuple{UInt16,UInt16},Any}()
    endpos = position(st) + sz
    while position(st) < endpos
        (gelt, data, vr) = read_element(st, dcm)
        if isequal(gelt, (0xFFFE, 0xE00D))
            break
        end
        item[gelt] = data
    end
    return item
end

# always little-endian, "encapsulated" iff sz==0xffffffff
function pixeldata_parse(st::IO, sz, vr::String, dcm)
    endian = dcm.endian
    dtype = determine_dtype(dcm, vr)
    yr = 1
    zr = 1
    # (0028,0010) defines number of rows
    f = get(dcm, (0x0028, 0x0010), nothing)
    if f !== nothing
        yr = Int(f)
    end
    # (0028,0011) defines number of columns
    f = get(dcm, (0x0028, 0x0011), nothing)
    if f !== nothing
        xr = Int(f)
    end
    # (0028,0012) defines number of planes
    f = get(dcm, (0x0028, 0x0012), nothing)
    if f !== nothing
        zr = Int(f)
    end
    # (0028,0008) defines number of frames
    f = get(dcm, (0x0028, 0x0008), nothing)
    if f !== nothing
        zr *= Int(f)
    end
    # (0x0028, 0x0002) defines number of samples per pixel
    f = get(dcm, (0x0028, 0x0002), nothing)
    if f !== nothing
        samples_per_pixel = Int(f)
    else
        samples_per_pixel = 1
    end
    if sz != 0xffffffff
        is_interleaved = get(dcm, (0x0028, 0x0006), nothing) == 0
        if is_interleaved
            data_dims = [samples_per_pixel, xr, yr, zr]
        else
            data_dims = [xr, yr, zr, samples_per_pixel]
        end
        data_dims = data_dims[data_dims.>1]
        data = Array{dtype}(undef, data_dims...)
        read!(st, data)
        # Permute because Julia is column-major while DICOM is row-major
        numdims = ndims(data)
        if numdims == 2
            perm = (2, 1)
        elseif numdims == 3
            perm = is_interleaved ? (3, 2, 1) : (2, 1, 3)
        elseif numdims == 4
            perm = is_interleaved ? (3, 2, 1, 4) : (2, 1, 3, 4)
        end
        data = permutedims(data, perm)
    else
        data = []
        while true
            grp = read_group_tag(st, endian)
            elt = read_element_tag(st, endian)
            itemlen = read_element_size(st, UInt32, endian)
            if grp == 0xFFFE && elt == 0xE0DD
                break
            end
            if grp != 0xFFFE || elt != 0xE000
                error("dicom: expected item tag in sequence")
            end
            push!(data, read!(st, Array{UInt8}(undef, itemlen)))
        end
        return data
    end
    return order.(data, endian)
end

function determine_dtype(dcm, vr = "OB")
    # (0x0028,0x0103) defines Pixel Representation
    is_signed = false
    f = get(dcm, (0x0028, 0x0103), nothing)
    if f !== nothing
        # Data is signed if f==1
        is_signed = f == 1
    end
    bit_type = 16
    # (0x0028,0x0100) defines Bits Allocated
    f = get(dcm, (0x0028, 0x0100), nothing)
    if f !== nothing
        bit_type = Int(f)
    else
        f = get(dcm, (0x0028, 0x0101), nothing)
        bit_type = f !== nothing ? Int(f) : vr == "OB" ? 8 : 16
    end
    if bit_type == 8
        dtype = is_signed ? Int8 : UInt8
    else
        dtype = is_signed ? Int16 : UInt16
    end
    return dtype
end

function undefined_length(st, vr)
    data = IOBuffer()
    w1 = w2 = 0
    while true
        # read until 0xFFFE 0xE0DD
        w1 = w2
        w2 = read(st, UInt16)
        if w1 == 0xFFFE
            if w2 == 0xE0DD
                break
            end
            write(data, w1)
        end
        if w2 != 0xFFFE
            write(data, w2)
        end
    end
    skip(st, 4)
    take!(data)
end

function dcm_write(fn::String, dcm::DICOMData; kwargs...)
    st = open(fn, "w+")
    dcm_write(st, dcm; kwargs...)
    close(st)
    return fn
end

function dcm_write(
    st::IO,
    dcm::DICOMData;
    preamble = true,
    aux_vr = Dict{Tuple{UInt16,UInt16},String}(),
)
    if preamble
        write(st, zeros(UInt8, 128))
        write(st, "DICM")
    end
    (is_explicit, endian) = determine_explicitness_and_endianness(dcm)
    for gelt in sort(collect(keys(dcm)))
        write_element(st, gelt, dcm[gelt], is_explicit, aux_vr)
    end
    return
end

function write_element(st::IO, gelt::Tuple{UInt16,UInt16}, data, is_explicit, aux_vr)
    if haskey(aux_vr, gelt)
        vr = aux_vr[gelt]
    else
        vr = lookup_vr(gelt)
    end
    if vr === empty_vr
        # Element tags ending in 0x0000 are not included in dcm_dicm.jl, their vr is UL
        if gelt[2] == 0x0000
            vr = "UL"
        elseif isodd(gelt[1]) && gelt[1] > 0x0008 && 0x0010 <= gelt[2] < +0x00FF
            # Private creator
            vr = "LO"
        elseif isodd(gelt[1]) && gelt[1] > 0x0008
            # Assume private
            vr = "UN"
        else
            error("dicom: unknown tag ", gelt)
        end
    end
    if gelt == (0x7FE0, 0x0010)
        return pixeldata_write(st, data, is_explicit)
    end

    if vr == "SQ"
        vr = is_explicit ? vr : empty_vr
        return dcm_store(st, gelt, s -> sequence_write(s, data, is_explicit), vr)
    end

    # Pack data into array container. This is to undo "data = data[1]" from read_element().
    if !isa(data, Array) && vr in ("FL", "FD", "SL", "SS", "UL", "US")
        data = [data]
    end

    data = isempty(data) ? UInt8[] :
        vr in ("OB", "OF", "OW", "ST", "LT", "UT") ? data :
        vr in ("AE", "CS", "SH", "LO", "UI", "PN", "DA", "DT", "TM") ?
        string_write(data, 0) :
        vr == "FL" ? convert(Array{Float32,1}, data) :
        vr == "FD" ? convert(Array{Float64,1}, data) :
        vr == "SL" ? convert(Array{Int32,1}, data) :
        vr == "SS" ? convert(Array{Int16,1}, data) :
        vr == "UL" ? convert(Array{UInt32,1}, data) :
        vr == "US" ? convert(Array{UInt16,1}, data) :
        vr == "AT" ? [data...] :
        vr in ("DS", "IS") ? string_write(map(string, data), 0) : data

    if !is_explicit && gelt[1] > 0x0002
        vr = empty_vr
    end

    dcm_store(st, gelt, s -> write(s, data), vr)
end

string_write(vals::Array{SubString{String}}, maxlen) =
    string_write(convert(Array{String}, vals), maxlen)
string_write(vals::SubString{String}, maxlen) = string_write(convert(String, vals), maxlen)
string_write(vals::Tuple{String,String}, maxlen) = string_write(collect(vals), maxlen)
string_write(vals::Char, maxlen) = string_write(string(vals), maxlen)
string_write(vals::String, maxlen) = string_write([vals], maxlen)
string_write(vals::Array{String,1}, maxlen) = join(vals, '\\')

dcm_store(st::IO, gelt::Tuple{UInt16,UInt16}, writef::Function) =
    dcm_store(st, gelt, writef, empty_vr)
function dcm_store(st::IO, gelt::Tuple{UInt16,UInt16}, writef::Function, vr::String)
    lentype = UInt32
    write(st, UInt16(gelt[1])) # Grp
    write(st, UInt16(gelt[2])) # Elt
    if vr !== empty_vr
        write(st, vr)
        if vr in ("OB", "OW", "OF", "SQ", "UT", "UN")
            write(st, UInt16(0))
        else
            lentype = UInt16
        end
    end
    # Write data first, then calculate length, then go back to write length
    p = position(st)
    write(st, zero(lentype)) # Placeholder for the data length
    writef(st)
    endp = position(st)
    # Remove placeholder's length, either 2 (UInt16) or 4 (UInt32) steps (1 step = 8 bits)
    if lentype == UInt32
        sz = endp - p - 4
    else
        sz = endp - p - 2
    end
    szWasOdd = isodd(sz) # If length is odd, round up - UInt8(0) will be written at end
    if szWasOdd
        sz += 1
    end
    seek(st, p)
    write(st, convert(lentype, max(0, sz)))
    seek(st, endp)
    if szWasOdd
        write(st, UInt8(0))
    end
end

sequence_write(st::IO, items::Array{Any,1}, evr::Bool) =
    sequence_write(st, convert(Array{Dict{Tuple{UInt16,UInt16},Any},1}, items), evr)
function sequence_write(st::IO, items::Array{Dict{Tuple{UInt16,UInt16},Any},1}, evr)
    for subitem in items
        if length(subitem) > 0
            dcm_store(st, (0xFFFE, 0xE000), s -> sequence_item_write(s, subitem, evr))
        end
    end
    write(st, UInt16[0xFFFE, 0xE0DD, 0x0000, 0x0000])
end

function sequence_item_write(st::IO, items::Dict{Tuple{UInt16,UInt16},Any}, evr)
    for gelt in sort(collect(keys(items)))
        write_element(st, gelt, items[gelt], evr, empty_dcm_dict)
    end
    write(st, UInt16[0xFFFE, 0xE00D, 0x0000, 0x0000])
end

function pixeldata_write(st, d, evr)
    # if length(el) > 1
    #     error("dicom: compression not supported")
    # end
    nt = eltype(d)
    vr = nt === UInt8 || nt === Int8 ? "OB" :
        nt === UInt16 || nt === Int16 ? "OW" :
        nt === Float32 ? "OF" : error("dicom: unsupported pixel format")
    # Permute because Julia is column-major while DICOM is row-major
    # !warn! This part assumes that Planar Configuration (tag: 0x0028, 0x0006) is not 0
    numdims = ndims(d)
    perm = numdims == 2 ? (2, 1) : numdims == 3 ? (2, 1, 3) : (2, 1, 3, 4)
    d = permutedims(d, perm)
    if evr !== false
        dcm_store(st, (0x7FE0, 0x0010), s -> write(s, d), vr)
    elseif vr != "OW"
        error("dicom: implicit VR only supports 16-bit pixels")
    else
        dcm_store(st, (0x7FE0, 0x0010), s -> write(s, d))
    end
end

end


text:
module BranchAndPrune

import Base: copy, eltype, iterate, IteratorSize, show
import Base: getindex, setindex!, delete!

export AbstractSearch, AbstractBreadthFirstSearch, AbstractDepthFirstSearch
export data
export copy, eltype, iterate, IteratorSize, nnodes

include("tree.jl")
include("search.jl")
include("forced_definitions.jl")

end  # module


text:
using LevelSetShapes
using AbbreviatedTypes
using LinearAlgebra
using Random: MersenneTwister
using Statistics: mean
using Test

include("simple.jl")

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
using Graphs
using BayesianNetworks
using Base.Test

function add_test_package(file::ASCIIString)
    start_time=time()
    try
        include(file)
        delta_time=round(time()-start_time,2)
        print("[")
        print_with_color(:green,"\e[92m\u2713\e[0m")
        print("] $file ($(delta_time)s)\n")
    catch e
        println("[\e[31mx\e[0m] $file")
        throw(e)
    end
end
add_test_package("utility.jl")
add_test_package("probability_distribution.jl")
add_test_package("cpd.jl")
add_test_package("bayesian_node.jl")
add_test_package("bayesian_network.jl")
add_test_package("bayesian_edge.jl")
add_test_package("probability_density_distribution.jl")
add_test_package("diagnosis_test.jl")

include("benchmarking/benchmark.jl")


text:
export PointOperation

export apply_operation
export domaintype
export isidentity, istranslation, ispoint

import LinearAlgebra
import LinearAlgebraX

"""
    PointOperation{S<:Real}

Represents point symmetry operation

# Fields
* `matrix`: rotation/mirror/inversion matrix
"""
struct PointOperation{S<:Real} <:AbstractSpaceSymmetryOperation{S}
    matrix::Matrix{S}

    @doc """
        PointOperation(matrix::AbstractMatrix{S}) where {S<:Real}

    Construct a point operation with `matrix`.
    """
    function PointOperation(matrix::AbstractMatrix{S}) where {S<:Real}
        return PointOperation{S}(matrix)
    end

    function PointOperation{S}(matrix::AbstractMatrix{S2}) where {S<:Real, S2<:Real}
        dim = size(matrix, 1)
        size(matrix, 2) != dim && throw(DimensionMismatch("matrix not square"))
        return new{S}(matrix)
    end
end

function Base.one(x::PointOperation{S}) where {S}
    return PointOperation(Matrix{S}(LinearAlgebra.I, size(x.matrix)))
end

# function Base.convert(::Type{PointOperation{S}}, obj::IdentityOperation{S}) where S
#     dim = dimension(obj)
#     return PointOperation{S}(Matrix(LinearAlgebra.I, dim, dim))
# end

function Base.convert(::Type{PointOperation{S}}, matrix::AbstractMatrix{S}) where S
    return PointOperation{S}(matrix)
end

# function Base.promote_rule(::Type{PointOperation{S}}, ::Type{IdentityOperation{S}}) where S
#     return PointOperation{S}
# end

## properties

"""
    isidentity(arg::PointOperation)

Check whether `arg` is identity.
"""
isidentity(arg::PointOperation) = isone(arg.matrix)

"""
    istranslation(arg::PointOperation)

Check whether `arg` is a translation operation, i.e. identity.
"""
istranslation(arg::PointOperation) = isone(arg.matrix)

"""
    ispoint(arg::PointOperation)

Check whether `arg` is a point operation. Always `true`.
"""
ispoint(arg::PointOperation) = true

"""
    dimension(arg::PointOperation)

Return spatial dimension of `arg`.
"""
dimension(arg::PointOperation) = size(arg.matrix, 1)

Base.hash(arg::PointOperation{S}, h::UInt) where S = hash(PointOperation{S}, hash(arg.matrix, h))

## operators
function Base.:(==)(lhs::PointOperation, rhs::PointOperation)
    return lhs.matrix == rhs.matrix
end
# function Base.:(==)(pop::PointOperation, ::IdentityOperation)
#     return isidentity(pop)
# end
# function Base.:(==)(::IdentityOperation, pop::PointOperation)
#     return isidentity(pop)
# end
function Base.:(==)(pop::PointOperation, top::TranslationOperation)
    return isidentity(pop) && isidentity(top)
end
function Base.:(==)(top::TranslationOperation, pop::PointOperation)
    return isidentity(pop) && isidentity(top)
end

function Base.:(*)(lhs::PointOperation, rhs::PointOperation)
    return PointOperation(lhs.matrix * rhs.matrix)
end

function Base.:(^)(lhs::PointOperation{S}, rhs::Integer) where S
    if rhs >= 0
        return PointOperation(lhs.matrix^rhs)
    else
        lhs_inv_matrix = Matrix{S}(LinearAlgebraX.invx(lhs.matrix))
        return PointOperation{S}(lhs_inv_matrix^(-rhs))
    end
end

function Base.inv(lhs::PointOperation{S}) where {S<:Real}
    return PointOperation{S}(LinearAlgebraX.invx(lhs.matrix))
end

## apply
"""
    apply_operation(symop::PointOperation{S}, coordinate::AbstractArray{S}) where {S}

Apply point operation to the coordinates.
"""
function apply_operation(symop::PointOperation{S}, coord::AbstractVector{S}) where S
    return symop.matrix * coord
end

function (symop::PointOperation{S})(coord::AbstractVector{S}) where S
    return symop.matrix * coord
end


text:
#Physical constants

module constants

export M_sun, M_earth, AU, pc, G, kB, h, amu, c_ang, cc, c_kms, molnames, arcsec, deg, lam0s, nyquist_factor, mu_gas, m_H
export Molecule, mols
export DiskTracerException, ModelException, ImageException
export DISKTRACER_VERSION

DISKTRACER_VERSION = v"0.1.0"

# Conversion from astronomical units to CGS units
const M_sun = 1.99e33 # [g]
const M_earth = 5.97219e27 # [g]
const AU = 1.4959787066e13 # [cm]
const pc = 3.0856776e18 # [cm]
const G = 6.67259e-8 # [cm3 g-1 s-2]
const kB = 1.380658e-16 # [erg K^-1] Boltzmann constant
const h = 6.6260755e-27 #erg s Planck constant
const c_ang = 2.99792458e18 # [A s^-1]
const cc = 2.99792458e10 # [cm s^-1]
const c_kms = 2.99792458e5 # [km s^-1]

# Conversion from degrees to radians
const deg = pi/180. # [radians]

const mu_gas = 2.37 # mean molecular weight of circumstellar gas
const m_H = 1.6733e-24 # g mass of hydrogen atom

#Atomic
const amu = 1.6605402e-24 # [g]

# To hold molecular constants
struct Molecule
    X_mol::Float64 # number fraction of molecule
    B0::Float64 # rigid rotor constant of molecule (Hz)
    mu::Float64 # (e.s.u. cm) permanent dipole moment of the molecule
    mol_weight::Float64 # (g) molecular weight
    nu_0::Float64 # (Hz) central frequency of transition
    l::Int # lower level of transition
    T_L::Float64 # (K) temperature equivalent of energy at lower transition
end

# Using numbers from Charlie Qi
# Number ratios measured relative to all single H nuclei (not H2) [unitless]
f_12CO = 7.5e-5
X_12CO = 2 * f_12CO
X_13CO = 1/69. * X_12CO
X_C18O = 1/557. * X_12CO

# molecular hydrogen number ratio to gas [unitless]
# [n_H2/n_gas] = 0.8
X_H2 = 0.8

# To go from rho_gas to n_12CO
# n_12CO = rho_gas * X_H2 *

# Old number densities calculated as
# 12CO = 7.5e-5 / (mu_gas * amu)

# New number density calculated as
# 12CO = 0.8 * (2 * 7.5e-5)/ (mu_gas * amu)

# new/old = 1.6
# This means that the new total disk mass we use will be 1/1.6 times smaller than the old time.

# Number ratios relative to average molecule (chi_CO)
# Multiply this against rho to get number of molecules
chi_12CO = X_H2 * X_12CO/(mu_gas * amu)
chi_13CO = X_H2 * X_13CO/(mu_gas * amu)
chi_C18O = X_H2 * X_C18O/(mu_gas * amu)

#CO
m_CO = 28.01 * amu #molecular weight of CO in g
m_12CO = 27.9949 * amu # [g]
m_13CO = 28.9983 * amu # [g]
m_C18O = 29.9992 * amu # [g]

CO12_21 = Molecule(chi_12CO, 57635.96e6, 1.1011e-19, m_12CO, 230.53800000e6, 1, 5.5321)
CO12_32 = Molecule(chi_12CO, 57635.96e6, 1.1011e-19, m_12CO, 345.79598990e6, 2, 16.5962)
CO13_21 = Molecule(chi_13CO, 55101.01e6, 1.1046e-19, m_13CO, 220.39868420e6, 1, 5.2888)
CO13_32 = Molecule(chi_13CO, 55101.01e6, 1.1046e-19, m_13CO, 330.58796530e6, 2, 15.8662)
CO18_21 = Molecule(chi_C18O, 54891.42e6, 1.1079e-19, m_C18O, 219.56035410e6, 1, 5.2686)
CO18_32 = Molecule(chi_C18O, 54891.42e6, 1.1079e-19, m_C18O, 329.33055250e6, 2, 15.8059)

# Species can be "12CO", "13CO", etc.
# Transition can be "3-2", "2-1", etc.
# The key to this dictionary is then species * transition

mols = Dict{String, Molecule}([("12CO2-1", CO12_21),
            ("12CO3-2", CO12_32),
            ("13CO2-1", CO13_21),
            ("13CO3-2", CO13_32),
            ("C18O2-1", CO18_21),
            ("C18O3-2", CO18_32)])

atomic_massses = Dict([("12CO", m_12CO), ("13CO", m_13CO), ("C18O", m_C18O)])

molnames = Dict{String, String}([("12CO", "co"), ("13CO", "13co"), ("C18O", "c18o")])

# Rest frame wavelengths
lam0s = Dict{String, Float64}([("12CO2-1", cc/230.538e9 * 1e4 ),
            ("13CO2-1", cc/220.39868e9 * 1e4),
            ("13CO3-2", cc/330.58797e9 * 1e4),
            ("C18O2-1", cc/219.56036e9 * 1e4),
            ("12CO3-2", cc/345.79599e9 * 1e4)]) # microns

# Used when determining the necessary number of pixels in an image, given distance. Anything below
# 2 is not Nyquist sampled. This is currently set to 2.2 to provide a degree of oversampling.
nyquist_factor = 2.2

# Abstract exception that serves as a catch-all for any error generated by the DiskJockey package.
abstract type DiskTracerException <: Exception end

# Define the various types of Errors that we can Raise
type ModelException <: DiskTracerException
    msg::String
end

type ImageException <: DiskTracerException
    msg::String
end

end # Module


text:
scan(x) = nothing

scan(c::Call) = foreach(scan, c.args)

function scan(x::TrackedArray)
  ref = x.ref += 1
  if ref == 1
    scan(x.f)
  else
    isdefined(x, :grad) || (x.grad = zeros(x.data))
  end
  return
end

back(c::Call, Δ) = back(c.func, Δ, c.args...)
back(::Call{Void}, Δ) = nothing

function back(x::TrackedArray, Δ)
  ref = x.ref -= 1
  if isdefined(x, :grad)
    x.grad .+= Δ
    ref == 0 && back(x.f, x.grad)
  else
    ref == 0 && back(x.f, Δ)
  end
  return
end

macro back(x, Δ)
  quote
    x = $(esc(x))
    istracked(x) && back(x, $(esc(Δ)))
  end
end

# Interface methods

function back!(x::TrackedArray, Δ)
  scan(x)
  back(x, Δ)
end

back!(x::TrackedScalar) = back!(x, 1)


text:
testvalues = [
    ([UInt8(0b10000001), UInt8(0b00000001)], 0b10000001),
    ([UInt8(0b00001111)], 0b00001111),
    ([UInt8(0b10011001), UInt8(0b00100101)], 0b0000110010100101),
    ([UInt8(0b10000001), UInt8(0b01111111)], 0b11111111)
]

for (input, output) in testvalues
    @test MIDI.readvariablelength(IOBuffer(input)) == output
end

for (output, input) in testvalues
    buf = IOBuffer()
    MIDI.writevariablelength(buf, Int64(input))
    @test takebuf_array(buf) == output
end


text:
# Creates a plot of the simulated JWST transit times for TRAPPIST-1:

using JLD2
include("plot_ttv.jl")
include("../../../src/regress.jl")

@load "../../../data/T1_JWST_all_transits_mass_error.jld2" JWST_sims data_obs

nplanet = JWST_sims[1].nplanet
count1 = JWST_sims[1].count
data_obs = JWST_sims[1].data_obs
elements_sim= JWST_sims[1].elements_in
tt_sim = JWST_sims[1].tt_actual
plot_ttv(data_obs,elements_sim,tt_sim)
subplots_adjust(hspace=0)
savefig("../T1_JWST_all_transits_tight.pdf",bbox_inches="tight")


text:
# Handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist.
# Testing sparse arrays.

using Base.Test
using Knet
isdefined(:MNIST) || include("mnist.jl")
include("mlp.jl")
atol2 = 0.01
rtol2 = 0.01
isapprox2(x,y)=isapprox(x.out,y.out;atol=atol2,rtol=rtol2)&&isapprox(x.dif,y.dif;atol=atol2,rtol=rtol2)
adense(x)=x

function mnistsparse(args=ARGS)
    info("Testing MNIST mlp with sparse arrays.")
    nbatch=100

    gcheck = 10
    net = cell(4)
    dtrn = cell(4)
    dtst = cell(4)
    iter = 0
    prog = mlp(layers=(64,10), loss=softmax, actf=relu, winit=Gaussian(0,.01), binit=Constant(0))
    results = Any[]

    for (fx,fy) in ((adense,adense), (adense,sparse), (sparse,adense), (sparse,sparse))
        iter += 1
        net[iter] = Net(prog)
        setopt!(net[iter], lr=0.5)
        dtrn[iter] = ItemTensor(fx(MNIST.xtrn), fy(MNIST.ytrn); batch=nbatch)
        dtst[iter] = ItemTensor(fx(MNIST.xtst), fy(MNIST.ytst); batch=nbatch)
        setseed(42)
        l=w=g=ltrn=atrn=ltst=atst=0
        for epoch=1:3
            # TODO: gcheck does not work ?
            (l,w,g) = train(net[iter], dtrn[iter]; gclip=0, gcheck=gcheck, getloss=true, getnorm=true, atol=atol2, rtol=rtol2)
            (ltrn,atrn,ltst,atst) = (test(net[iter], dtrn[iter]),
                                     accuracy(net[iter], dtrn[iter]),
                                     test(net[iter], dtst[iter]),
                                     accuracy(net[iter], dtst[iter]))
        end
        @show (fx,fy,l,w,g,ltrn,atrn,ltst,atst)
        @show map(isapprox2, params(net[1]), params(net[iter]))
        push!(results, (fx,fy,l,w,g,ltrn,atrn,ltst,atst))
    end
    return results
end

!isinteractive() && !isdefined(:load_only) && mnistsparse(ARGS)

# Reference output for debugging.  The last two results with sparse input are not stable.

# [dy_052@hpc3001 examples]$ julia mnistsparse.jl
# INFO: Loading MNIST...
#   5.478525 seconds (366.45 k allocations: 503.265 MB, 1.74% gc time)
#   8.876875 seconds (9.64 M allocations: 427.910 MB, 2.03% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (adense,adense,0.10628127f0,24.865437f0,3.5134742f0,0.100041345f0,0.9681833333333333,0.114785746f0,0.9641)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[true,true,true,true]
#   8.316890 seconds (5.36 M allocations: 251.544 MB, 0.93% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (adense,sparse,0.1062698f0,24.866688f0,3.5134742f0,0.100718f0,0.9679333333333333,0.1149149f0,0.9642)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[true,true,true,true]
#  12.682341 seconds (5.63 M allocations: 629.662 MB, 1.61% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (sparse,adense,0.10628127f0,24.865438f0,3.5134711f0,0.100041375f0,0.9681833333333333,0.11478577f0,0.9641)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[true,true,true,true]
#  14.437739 seconds (5.15 M allocations: 622.468 MB, 1.36% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (sparse,sparse,0.106173664f0,24.870913f0,3.5582156f0,0.10027795f0,0.9682,0.11406385f0,0.9646)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[false,false,false,true]

# Here is another run showing the unstability:

# [dy_052@hpc3001 examples]$ julia mnistsparse.jl
# INFO: Loading MNIST...
#   5.480183 seconds (366.45 k allocations: 503.265 MB, 1.72% gc time)
#   8.797816 seconds (9.64 M allocations: 427.910 MB, 2.05% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (adense,adense,0.10628127f0,24.865437f0,3.5134742f0,0.100041345f0,0.9681833333333333,0.114785746f0,0.9641)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[true,true,true,true]
#   7.973043 seconds (5.35 M allocations: 251.381 MB, 0.88% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (adense,sparse,0.1062698f0,24.866688f0,3.5134742f0,0.100718f0,0.9679333333333333,0.1149149f0,0.9642)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[true,true,true,true]
#  11.925001 seconds (5.64 M allocations: 629.857 MB, 1.22% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (sparse,adense,0.1061508f0,24.87466f0,3.5049448f0,0.09875954f0,0.9686833333333333,0.11317172f0,0.9654)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[false,true,false,true]
#  14.341379 seconds (5.15 M allocations: 621.469 MB, 1.00% gc time)
# (fx,fy,l,w,g,ltrn,atrn,ltst,atst) = (sparse,sparse,0.10627592f0,24.866186f0,3.5134723f0,0.09992511f0,0.9681666666666666,0.114817984f0,0.9641)
# map(isapprox1,params(net[1]),params(net[iter])) = Bool[true,true,true,true]


text:
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "group1"
        help = "One end point of the edge for which to analyze context, where "*
               "| is used to separate mulitple ids. (i.e. ENCFF000ARO|ENCFF000ARK)."
        arg_type = ASCIIString
        required = true
    "group2"
        help = "The other end point of the edge for which to analyze context, where "*
               "| is used to separate mulitple ids. (i.e. ENCFF000ARO|ENCFF000ARK)."
        arg_type = ASCIIString
        required = true
    "data_bundles"
        help = "Location of the ChromNet data bundles (i.e. ENCODE_buildX.ChromNet.jld)"
        arg_type = ASCIIString
        nargs = '+'
        required = true
    "--output", "-o"
        help = "Where to save the edge context BED file"
        default = "-"
    "--load-cor"
        help = "A copy of the raw data correlation matrix"
        metavar = "COR_FILE"
    "--quiet", "-q"
        help = "Don't print anything"
        action = :store_true
end
args = parse_args(s)

using ChromNet
using JSON
using JLD
using SamIO
using ProgressMeter

# open the data bundle and load the metadata
dataBundles = [jldopen(x) for x in args["data_bundles"]]
metadata = Dict()
for bundle in dataBundles
    metadata = merge(metadata, read(bundle, "metadata"))
end
ids = vcat([read(bundle, "ids") for bundle in dataBundles]...)
numBins = round(Int, ceil(sum(ReferenceContigs_hg38.sizes) / 1000))

# find the indexes of the groups this edge is between
group1Inds = map(x->findfirst(ids, x), split(args["group1"], "|"))
group2Inds = map(x->findfirst(ids, x), split(args["group2"], "|"))
@assert all(group1Inds .> 0) "Some ids from group 1 were not found: $(group1Inds[group1Inds .== 0])"
@assert all(group2Inds .> 0) "Some ids from group 2 were not found: $(group2Inds[group2Inds .== 0])"

# compute the joint correlation matrix
C = nothing
if typeof(args["load-cor"]) != Void
    C,header = readdlm(args["load-cor"], header=true)
    @assert all(header .== ids) "Loaded correlation matrix does not match data! Perhaps a different file name ordering?"
else
    C = streaming_cor(dataBundles, quiet=args["quiet"])
end
cov2cor!(C)

# compute the edge impact
edgeScores = streaming_edgeimpact(dataBundles, C, group1Inds, group2Inds, quiet=args["quiet"])
map(close, dataBundles)

# send the data in BED format to the output stream
outStream = STDOUT
if args["output"] != "-"
    outStream = open(args["output"], "w")
end
contigIndex = 1
pos = 1
ctgs = ReferenceContigs_hg38
for i in 1:numBins
    while pos > ctgs.offsets[contigIndex] + ctgs.sizes[contigIndex]
        contigIndex += 1
    end
    println(outStream, ctgs.names[contigIndex], "\t", pos - ctgs.offsets[contigIndex], "\t", pos - ctgs.offsets[contigIndex] + 999, "\t*\t", edgeScores[i])
    pos += 1000
end


text:
module TestingOIFITS

using Test, Printf
using OIFITS, FITSIO
using OIFITS: extname, get_format, is_same, fix_name

const dir = @__DIR__

const files = ("contest-2004-obj1.oifits" ,"contest-2004-obj2.oifits",
               "contest-2008-binary.oifits", "contest-2008-obj1-H.oifits",
               "contest-2008-obj1-J.oifits", "contest-2008-obj1-K.oifits",
               "contest-2008-obj2-H.oifits", "contest-2008-obj2-J.oifits",
               "contest-2008-obj2-K.oifits", "contest-2008-obj1-J.oifits")

const counter = Ref(0)

const quiet = true

const tempfile = let tup = mktemp()
    close(tup[2])
    tup[1]
end

println("tempfile = \"$tempfile\"")

function tryload(dir, file)
    global counter
    try
        data = OIDataSet(joinpath(dir, file))
        counter[] += 1
        quiet || @info "file \"$file\" successfully loaded"
        return true
    catch
        @warn "failed to load \"$file\""
        return false
    end
end

@testset "Strings" begin
    @test  is_same("", "")
    @test  is_same("", " ")
    @test  is_same("\t", "")
    @test  is_same(" ZCMa", " zcma")
    @test  is_same(" ZCMa", " zcma ")
    @test !is_same("ZCMa", " zcma")
    @test  is_same("SAO-206462", "sao-206462")
    @test  is_same("SAO-206462", "sao-206462\t")
    @test  is_same("SAO 206462", "sao 206462")
    @test  is_same("SAO 206462", "sao 206462\t")
    @test !is_same("SAO 206462", "sao\t206462")

    @test fix_name("") === ""
    @test fix_name(" ") === ""
    @test fix_name("a\t ") === "A"
    @test fix_name("\tbcd\t ") === "\tBCD"
    @test fix_name("\téëï\t ") === "\tÉËÏ"
end

get_field_type(T::Type{<:OIDataBlock}, sym::Symbol) = fieldtype(T, sym)
get_field_type(::Type{<:OI_TARGET}, sym::Symbol) =
    fieldtype((sym === :revn ? OI_TARGET : OITargetEntry), sym)

@testset "OI type definitions and formats" begin
    for T in (OI_ARRAY,
              OI_CORR,
              OI_FLUX,
              OI_INSPOL,
              OI_T3,
              OI_TARGET,
              OI_VIS,
              OI_VIS2,
              OI_WAVELENGTH)

        @testset "$(extname(T))" begin
            for rev in 3:-1:1
                spec = get_format(T, rev; throw_errors=false)
                if spec === nothing
                    continue
                end
                for s in spec
                    @test s.type ∈ (:A, :C, :D, :E, :I, :J, :L)
                    S = get_field_type(T, s.symb)
                    rank = ndims(s)
                    if rank == 0
                        # Keyword.
                        if s.type === :A
                            @test S === String
                        elseif s.type === :C
                            @test S <: Complex{<:AbstractFloat}
                        elseif s.type === :E || s.type === :D
                            @test S <: AbstractFloat
                        elseif s.type === :I || s.type === :J
                            @test S <: Integer
                        elseif s.type === :L
                            @test S <: Bool
                        end
                    else
                        # Column.
                        if s.type === :A
                            @test rank == 1
                            @test S === (T === OI_TARGET ? String : Vector{String})
                        else
                            @test rank == (T === OI_TARGET ? 1 : ndims(S))
                            if s.type === :C
                                @test eltype(S) <: Complex{<:AbstractFloat}
                            elseif s.type ∈ (:D, :E)
                                @test eltype(S) <: AbstractFloat
                            elseif s.type ∈ (:I, :J)
                                @test eltype(S) <: Int
                            elseif s.type === :L
                                @test eltype(S) <: Bool
                            end
                        end
                    end
                end
            end
        end
    end
end

@testset "Read/write OI-FITS files" begin
    # Each file infdividually.
    for file in files
        A = OIDataSet(joinpath(dir, file))
        @test isa(A, OIDataSet)
        isfile(tempfile) && rm(tempfile; force=true)
        write(tempfile, A)
        B = read(OIDataSet, tempfile)
        @test isa(B, OIDataSet)
        @test_throws ErrorException write(tempfile, B)
    end

    # Merge several files.
    A = OIDataSet(map(x -> joinpath(dir, x), files)...)
    @test isa(A, OIDataSet)
    @test A.instr[1].name ===  A.instr[1].insname
    @test A.array[1].name ===  A.array[1].arrname
    name = A.instr[end].name
    @test isa(A.instr[name], OI_WAVELENGTH)
    @test A.instr[name] === A.instr[end]
    @test A.instr[uppercase(name*" ")] === A.instr[end]
    @test A.instr[lowercase(name*" ")] === A.instr[end]
    @test_throws KeyError A.instr[" "*name]
    @test get(A.instr, " "*name, undef) === undef
    name = A.array[end].name
    @test isa(A.array[name], OI_ARRAY)
    @test A.array[name] === A.array[end]
    @test A.array[uppercase(name*" ")] === A.array[end]
    @test A.array[lowercase(name*" ")] === A.array[end]
    @test_throws KeyError A.array[" "*name]
    @test get(A.array, " "*name, undef) === undef

    # Check for "private" and non-existing fields.
    @test_throws ErrorException A.target_dict
    @test_throws ErrorException A.target_id_map
    @test_throws ErrorException A.array_dict
    @test_throws ErrorException A.instr_dict
    @test_throws ErrorException A.correl_dict

    # Check for fields being read-only.
    @test_throws ErrorException A.array = A.array

    # Clear contents.
    empty!(A)
    @test A.target.revn == 0
    @test length(A.target) == 0
    @test length(A.target.list) == 0
    @test length(A.array) == 0
    @test length(A.instr) == 0
    @test length(A.correl) == 0
    @test length(A.vis) == 0
    @test length(A.vis2) == 0
    @test length(A.t3) == 0
    @test length(A.flux) == 0
    @test length(A.inspol) == 0
    @test length(getfield(A, :target_dict)) == 0
    @test length(getfield(A, :target_id_map)) == 0
    @test length(getfield(A, :array_dict)) == 0
    @test length(getfield(A, :instr_dict)) == 0
    @test length(getfield(A, :correl_dict)) == 0

    # Access extensions individually.
    f = FITS(joinpath(dir, "contest-2004-obj1.oifits"))
    for i in 2:length(f)
        db = OIDataBlock(f[i])
        @test extname(f[i]) == db.extname
    end
    @test extname(f[2]) == OI_ARRAY(f[2]).extname
    @test extname(f[3]) == OI_TARGET(f[3]).extname
    @test extname(f[4]) == OI_WAVELENGTH(f[4]).extname
    @test extname(f[5]) == OI_VIS2(f[5]).extname
    @test extname(f[6]) == OI_T3(f[6]).extname
end

@testset "OI_TARGET methods" begin
    data = OIDataSet(joinpath(dir,files[1]))
    A = data.target
    @test isa(A, OI_TARGET)
    @test eltype(A) === OITargetEntry
    @test length(A) == length(A.list)
    @test ndims(A) == 1
    @test IndexStyle(A) === IndexLinear()
    @test firstindex(A) == 1
    @test lastindex(A) == length(A)
    @test size(A) == (length(A),)
    @test size(A,1) == length(A)
    @test size(A,2) == 1
    @test axes(A) == (1:length(A),)
    @test axes(A,1) == 1:length(A)
    @test eachindex(A) === Base.OneTo(length(A))
    @test keys(Int, A) === Base.OneTo(length(A))
    for (key, val) in zip(keys(String, A), values(A))
        @test A[key] === val
    end
    i = 0
    for tgt in A
        i += 1
        @test isa(tgt, OITargetEntry)
        @test A[i] === tgt
        @test haskey(A, i) == true
        for key in (tgt.target,
                    uppercase(tgt.target),
                    lowercase(tgt.target),
                    tgt.target*"  \t")
            @test haskey(A, key)
            @test !haskey(A, key*"+")
            @test A[key] === tgt
            @test get(A, key, nothing) === tgt
            @test get(A, key*"+", nothing) === nothing
        end
    end
end

isfile(tempfile) && rm(tempfile; force=true)

end # module

nothing


text:
using CSV, DataFrames

function main()
  file = "participantes.csv"
  data = CSV.read(file, DataFrame, header=true)
  data = sort(data, ["Nome Completo"])
  isdir("participantes") || mkdir("participantes")

  n = size(data, 1)
  for i = 1:n
    nome = titlecase(data[i,2])
    horas = data[i,3]
    filename = join(split(nome), "_")
    d1 = "gravity Center text 0,-230 \"$nome\""
    d2 = "gravity Center text 0,160 \"$horas horas.\""

    tmp = "tmp.png"

    cmd = `convert -pointsize 80 -draw $d1 -font LibreBaskerville-Italic.ttf template-participantes.png tmp.png`
    run(cmd)
    cmd = `mogrify -pointsize 50 -draw $d2 -font LibreBaskerville-Italic.ttf tmp.png`
    run(cmd)
    cmd = `convert tmp.png participantes/$filename.pdf`
    run(cmd)
  end
end

main()

max_stars_repo_name: CiDAMO/scripts-certificados

max_issues_repo_name: CiDAMO/scripts-certificados

max_forks_repo_name: CiDAMO/scripts-certificados


text:
"""
A pool population is a well-mixing finite community, i.e. one in
which every interactor is connected to every other interactor and
the probability of interaction is uniform.
"""
mutable struct PoolPopulation{T<:AbstractInteractor} <: AbstractPopulation
  census::Array{T}
  nactions::Int64
end

"""
    PoolPopulation{T}()

Construct an empty pool population of interactors of type `T`.

## Examples

    PoolPopulation{EmptyInteractor}()
"""
function PoolPopulation{T}() where {T<:AbstractInteractor}
  PoolPopulation(Array{T}(undef, 0), 0)
end


text:
# this file test the validity by comparing with the widely used `libgsl`

# wapper of gsl 3nj functions
function gsl3j(dj1::Int, dj2::Int, dj3::Int, dm1::Int, dm2::Int, dm3::Int)
    ccall(
        (:gsl_sf_coupling_3j, "libgsl"),
        Cdouble,
        (Cint, Cint, Cint, Cint, Cint, Cint),
        dj1, dj2, dj3, dm1, dm2, dm3
    )
end

function gsl6j(dj1::Int, dj2::Int, dj3::Int, dj4::Int, dj5::Int, dj6::Int)
    ccall(
        (:gsl_sf_coupling_6j, "libgsl"),
        Cdouble,
        (Cint, Cint, Cint, Cint, Cint, Cint),
        dj1, dj2, dj3, dj4, dj5, dj6
    )
end

function gsl9j(dj1::Int, dj2::Int, dj3::Int,
                dj4::Int, dj5::Int, dj6::Int,
                dj7::Int, dj8::Int, dj9::Int)
    ccall(
        (:gsl_sf_coupling_9j, "libgsl"),
        Cdouble,
        (Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        dj1, dj2, dj3, dj4, dj5, dj6, dj7, dj8, dj9
    )
end

function check_6j(dj1::Int, dj2::Int, dj3::Int, dj4::Int, dj5::Int, dj6::Int)
    check_couple(dj1, dj2, dj3) & check_couple(dj1, dj5, dj6) &
    check_couple(dj4, dj2, dj6) & check_couple(dj4, dj5, dj3)
end

function check_9j(dj1::Int, dj2::Int, dj3::Int, dj4::Int, dj5::Int, dj6::Int, dj7::Int, dj8::Int, dj9::Int)
    check_couple(dj1, dj2, dj3) & check_couple(dj4, dj5, dj6) & check_couple(dj7, dj8, dj9) &
    check_couple(dj1, dj4, dj7) & check_couple(dj2, dj5, dj8) & check_couple(dj3, dj6, dj9)
end

# Because the gsl compute wigner 3nj symbols using float point number,
# sometimes my code give out zero, but gsl give out a very small number.
# In this condition, `≈` operator will get false result. So it should be exclude.
function test_3j_with_gsl(test_range::AbstractArray)
    for dj1 in test_range
    for dj2 in test_range
    for dj3 in test_range
    for dm1 in -dj1:2:dj1
    for dm2 in -dj2:2:dj2
        dm3 = -dm1-dm2
        if check_couple(dj1, dj2, dj3) & check_jm(dj3, dm3)
            gsl = gsl3j(dj1, dj2, dj3, dm1, dm2, dm3)
            my = float(d3j(dj1, dj2, dj3, dm1, dm2, dm3))
            @test (my == 0) | (gsl ≈ my)
        end
    end end end end end
end

function test_6j_with_gsl(test_range::AbstractArray)
    for j1 in test_range
    for j2 in test_range
    for j3 in test_range
    for j4 in test_range
    for j5 in test_range
    for j6 in test_range
        if check_6j(j1, j2, j3, j4, j5, j6)
            gsl = gsl6j(j1, j2, j3, j4, j5, j6)
            my = float(d6j(j1,j2,j3,j4,j5,j6))
            @test (my == 0) | (gsl ≈ my)
        end
    end end end end end end
end

function test_9j_with_gsl(test_range::AbstractArray)
    for j1 in test_range
    for j2 in test_range
    for j3 in test_range
    for j4 in test_range
    for j5 in test_range
    for j6 in test_range
    for j7 in test_range
    for j8 in test_range
    for j9 in test_range
        if check_9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)
            gsl = gsl9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)
            my = float(d9j(j1,j2,j3,j4,j5,j6,j7,j8,j9))
            @test (my == 0) | (gsl ≈ my)
        end
    end end end end end end end end end
end


text:
############################################################################################
"""
$(TYPEDEF)

PMCMC tuning container.

# Fields
$(TYPEDFIELDS)
"""
struct PMCMCTune{T<:Tagged} <: AbstractTune
    "Tagged Model parameter."
    tagged::T
    function PMCMCTune(
        tagged::T
    ) where {T<:Tagged}
        return new{T}(tagged)
    end
end

############################################################################################
#export
export PMCMCTune


text:
"""
"""
struct IdentityVector{T<:Integer} <: AbstractVector{T}
  length::T
end

@inline function getindex(c::IdentityVector{T},i::Integer) where T
  @check i > 0
  @check i <= c.length
  j::T = i
  j
end

size(c::IdentityVector) = (c.length,)

IndexStyle(::Type{<:IdentityVector}) = IndexLinear()


text:
@everywhere using MIPVerify
@everywhere using Gurobi
@everywhere using MAT

using ArgParse

function buildModel(dataset:: String, scale:: String, weights:: Dict)
    println("Building $dataset $scale...")
    if dataset == "MNIST" && scale == "small"
        nn = Sequential([
            get_conv_params(weights, "0", (4, 4, 1, 16), expected_stride = 2),
            ReLU(interval_arithmetic),
            get_conv_params(weights, "2", (4, 4, 16, 32), expected_stride = 2),
            ReLU(),
            Flatten([1, 3, 2, 4]),
            get_matrix_params(weights, "5", (1568, 100)),
            ReLU(),
            get_matrix_params(weights, "7", (100, 10))], "$(dataset)_$(scale)")
        return nn
    elseif dataset == "MNIST" && scale == "large"
        nn = Sequential([
            get_conv_params(weights, "0", (3, 3, 1, 32), expected_stride = 1),
            ReLU(interval_arithmetic),
            get_conv_params(weights, "2", (4, 4, 32, 32), expected_stride = 2),
            ReLU(),
            get_conv_params(weights, "4", (3, 3, 32, 64), expected_stride = 1),
            ReLU(),
            get_conv_params(weights, "6", (4, 4, 64, 64), expected_stride = 2),
            ReLU(),
            Flatten([1, 3, 2, 4]),
            get_matrix_params(weights, "9", (3136, 512)),
            ReLU(),
            get_matrix_params(weights, "11", (512, 512)),
            ReLU(),
            get_matrix_params(weights, "13", (512, 10))], "$(dataset)_$(scale)")
        return nn
    elseif dataset == "MNIST" && scale == "base"
        nn = Sequential([
            get_conv_params(weights, "0", (4, 4, 1, 32), expected_stride = 2),
            ReLU(interval_arithmetic),
            get_conv_params(weights, "2", (4, 4, 32, 64), expected_stride = 2),
            ReLU(),
            Flatten([1, 3, 2, 4]),
            get_matrix_params(weights, "5", (3136, 1024)),
            ReLU(),
            get_matrix_params(weights, "7", (1024, 10))], "$(dataset)_$(scale)")
        return nn
    elseif dataset == "CIFAR10" && scale == "small"
        nn = Sequential([
            get_conv_params(weights, "0", (4, 4, 3, 16), expected_stride = 2),
            ReLU(interval_arithmetic),
            get_conv_params(weights, "2", (4, 4, 16, 32), expected_stride = 2),
            ReLU(),
            Flatten([1, 3, 2, 4]),
            get_matrix_params(weights, "5", (2048, 100)),
            ReLU(),
            get_matrix_params(weights, "7", (100, 10))], "$(dataset)_$(scale)")
        return nn
    elseif dataset == "CIFAR10" && scale == "large"
        nn = Sequential([
            get_conv_params(weights, "0", (3, 3, 3, 32), expected_stride = 1),
            ReLU(interval_arithmetic),
            get_conv_params(weights, "2", (4, 4, 32, 32), expected_stride = 2),
            ReLU(),
            get_conv_params(weights, "4", (3, 3, 32, 64), expected_stride = 1),
            ReLU(),
            get_conv_params(weights, "6", (4, 4, 64, 64), expected_stride = 2),
            ReLU(),
            Flatten([1, 3, 2, 4]),
            get_matrix_params(weights, "9", (4096, 512)),
            ReLU(),
            get_matrix_params(weights, "11", (512, 512)),
            ReLU(),
            get_matrix_params(weights, "13", (512, 10))], "$(dataset)_$(scale)")
        return nn
    elseif dataset == "CIFAR10" && scale == "resnet"
        nn = SkipSequential([
            get_conv_params(weights, "0", (3, 3, 3, 16), expected_stride=1),
            ReLU(interval_arithmetic),
            get_conv_params(weights, "2/Ws/0", (3, 3, 16, 16), expected_stride=1),
            ReLU(),
            SkipBlock([
                get_conv_params(weights, "4/Ws/0", (1, 1, 16, 16), expected_stride=1),
                Zero(),
                get_conv_params(weights, "4/Ws/2", (3, 3, 16, 16), expected_stride=1)
            ]),
            ReLU(),
            get_conv_params(weights, "6/Ws/0", (3, 3, 16, 16), expected_stride=1),
            ReLU(),
            SkipBlock([
                get_conv_params(weights, "8/Ws/0", (1, 1, 16, 16), expected_stride=1),
                Zero(),
                get_conv_params(weights, "8/Ws/2", (3, 3, 16, 16), expected_stride=1)
            ]),
            ReLU(),
            get_conv_params(weights, "10/Ws/0", (4, 4, 16, 32), expected_stride=2),
            ReLU(),
            SkipBlock([
                get_conv_params(weights, "12/Ws/0", (2, 2, 16, 32), expected_stride=2),
                Zero(),
                get_conv_params(weights, "12/Ws/2", (3, 3, 32, 32), expected_stride=1)
            ]),
            ReLU(),
            get_conv_params(weights, "14/Ws/0", (4, 4, 32, 64), expected_stride=2),
            ReLU(),
            SkipBlock([
                get_conv_params(weights, "16/Ws/0", (2, 2, 32, 64), expected_stride=2),
                Zero(),
                get_conv_params(weights, "16/Ws/2", (3, 3, 64, 64), expected_stride=1)
            ]),
            ReLU(),
            Flatten([1, 3, 2, 4]),
            get_matrix_params(weights, "19", (4096, 1000)),
            ReLU(),
            get_matrix_params(weights, "21", (1000, 10))], "$(dataset)_$(scale)")
        return nn
    end
end

function initModels(path = "final_models":: String)
    files = [fname for fname in readdir(path) if endswith(fname, ".mat")]
    models = Dict{String, NeuralNet}()
    configs = Dict{String, Any}()
    for f in files
        things = [String(x) for x in split((rsplit(f, "."; limit=2)[1]), "_")]
        weights = matread(path * "/" * f)
        model = buildModel(things[1], things[3], weights)
        if model != nothing
            push!(models, f => model)
            push!(configs, f => things)
        end
    end
    return models, configs
end

@everywhere function get_label(y::Array{<:Real, 1}, test_index::Integer)::Int
    return y[test_index]
end

@everywhere function get_image(x::Array{T, 4}, test_index::Integer)::Array{T, 4} where {T<:Real}
    return x[test_index:test_index, :, :, :]
end

@everywhere function get_max_index(
    x::Array{<:Real, 1})::Integer
    return findmax(x)[2]
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--sampleNum"
            help = "number of samples"
            arg_type = Int
            default = 99999
        "--workers"
            help = "number of processes"
            arg_type = Int
            default = 1
        "dset"
            help = "dataset"
            arg_type = String
            required = true
        "epsilon"
            help = "epsilon"
            arg_type = String
            required = true
        "scale"
            help = "model scale / type"
            arg_type = String
            required = true
    end
    return parse_args(s)
end

@everywhere function work(work_id, dataset, datasetName, epsilon, m, l, r, tot)
    println("Proc. start $(work_id) $(l) $(r) $(tot)")

    tot_samples = 0
    clean_correct = 0
    robust_correct = 0
    tot_time = 0.0
    if l <= tot
        r = min(r, tot)
        tot_samples = max(r - l + 1, 0)
        for sample_index in l: r
            x0 = get_image(dataset.test.images, sample_index)
            if datasetName == "CIFAR10"
                # normalize
                mean = [0.485, 0.456, 0.406]
                std = [0.225, 0.225, 0.225]
                for i in 1: 3
                    x0[:, :, :, i] = (x0[:, :, :, i] - mean[i]) / std[i]
                end
                println(size(x0))
                println(minimum(x0))
                println(maximum(x0))
            end
            actual_label = get_label(dataset.test.labels, sample_index)
            actual_index = actual_label + 1
            predicted_output = x0 |> m
            predicted_index = predicted_output |> get_max_index
            sample_time = 0.0
            if predicted_index == actual_index
                clean_correct += 1
                try
                    d = find_adversarial_example(m, x0,
                        actual_index, GurobiSolver(TimeLimit=1200, BestObjStop=epsilon, OutputFlag=0, Threads=1),
                        invert_target_selection=true,
                        pp=MIPVerify.LInfNormBoundedPerturbationFamily(epsilon),
                        norm_order=Inf, tolerance=0.0, rebuild=false, tightening_algorithm=lp,
                        tightening_solver=GurobiSolver(TimeLimit=20, OutputFlag=0, Threads=1),
                        solve_if_predicted_in_targeted = false)
                    sample_time = float(d[:TotalTime])
                    if d[:SolveStatus] == :InfeasibleOrUnbounded
                        robust_correct += 1
                    end
                finally
                end
            end
            tot_time += sample_time
            println("\n\033[31m[worker #$(work_id)] [$(sample_index-l+1)/$(tot_samples)] Robust: $(robust_correct / (sample_index-l+1)) Clean: $(clean_correct / (sample_index-l+1)) Time: $(tot_time / (sample_index-l+1))\033[0m\n")
        end
    end
    return [tot_samples, clean_correct, robust_correct, tot_time]
end

# ----- main -----
function main()
    args = parse_commandline()

    io = open("$(args["dset"])_$(args["epsilon"])_$(args["scale"])_milp.test", "a")

    println("Start...")
    models, configs = initModels()

    num_samples = args["sampleNum"]
    m = nothing
    for (k, v) in models
        if configs[k][1] == args["dset"] && configs[k][2] == args["epsilon"] && configs[k][3] == args["scale"]
            m = v
            break
        end
    end
    if args["dset"] == "MNIST" && args["epsilon"] == "0.1" && args["scale"] == "baseline"
        m = MIPVerify.get_example_network_params("MNIST.WK17a_linf0.1_authors")
    end
    @assert m != nothing
    println("Now working on config $(args["dset"]) $(args["epsilon"]) $(args["scale"])")
    dataset = MIPVerify.read_datasets(lowercase(args["dset"]))
    num_samples = min(num_samples, MIPVerify.num_samples(dataset.test))

    workers = args["workers"]
    samplePerWorker = Int(ceil(num_samples / workers))
    runs = [@spawn work(i, dataset, args["dset"], float(args["epsilon"]), m, (i-1) * samplePerWorker + 1, i * samplePerWorker, num_samples) for i in 1: workers]
    anss = [fetch(item) for item in runs]

    print(anss)

    stat = [sum([anss[i][j] for i in 1: workers]) for j in 1: 4]
    @assert stat[1] == num_samples
    println("------")
    println("clean error: $(1.0 - stat[2] / num_samples)")
    println("robust error: $(1.0 - stat[3] / num_samples)")
    println("avg time: $(stat[4] / num_samples)")
    println("------")

    write(io, "$(num_samples)\nclean error: $(1.0 - stat[2] / num_samples)\nrobust error: $(1.0 - stat[3] / num_samples)\navg time: $(stat[4] / num_samples)\n")
    close(io)
end

main()

max_stars_repo_name: llylly/Robustra

max_issues_repo_name: llylly/Robustra

max_forks_repo_name: llylly/Robustra


text:
"""
This package exports
* `BOpt`, `boptimize!`
* acquisition types: `ExpectedImprovement`, `ProbabilityOfImprovement`, `UpperConfidenceBound`, `ThompsonSamplingSimple`, `MutualInformation`
* scaling of standard deviation in `UpperConfidenceBound`: `BrochuBetaScaling`, `NoBetaScaling`
* GP hyperparameter optimizer: `MAPGPOptimizer`, `NoModelOptimizer`
* Initializer: `ScaledSobolIterator`, `ScaledLHSIterator`
* optimization sense: `Min`, `Max`
* verbosity levels: `Silent`, `Timings`, `Progress`
* helper: maxduration!, maxiterations!

Use the REPL help, e.g. `?Bopt`, to get more information.
"""
module BayesianOptimization
import NLopt, GaussianProcesses, Sobol
import Sobol: SobolSeq
import GaussianProcesses: GPBase, GPE
import ElasticArrays: ElasticArray
using ForwardDiff, DiffResults, Random, Dates, SpecialFunctions, TimerOutputs
export BOpt, ExpectedImprovement, ProbabilityOfImprovement,
UpperConfidenceBound, ThompsonSamplingSimple, MutualInformation, boptimize!,
MAPGPOptimizer, NoModelOptimizer, Min, Max, BrochuBetaScaling, NoBetaScaling,
Silent, Timings, Progress, ScaledSobolIterator, ScaledLHSIterator, maxduration!,
maxiterations!

ENABLE_TIMINGS = true

abstract type ModelOptimizer end
"""
Don't optimize the model ever.
"""
struct NoModelOptimizer <: ModelOptimizer end
optimizemodel!(o::NoModelOptimizer, model) = Nothing

include("utils.jl")
include("acquisitionfunctions.jl")
include("acquisition.jl")
include("models/gp.jl")

@enum Sense Min=-1 Max=1
@enum Verbosity Silent=0 Timings=1 Progress=2

mutable struct BOpt{F,M,A,AO,MO,Ti}
    func::F
    sense::Sense
    model::M
    acquisition::A
    acquisitionoptions::AO
    modeloptimizer::MO
    lowerbounds::Array{Float64, 1}
    upperbounds::Array{Float64, 1}
    observed_optimum::Float64
    observed_optimizer::Array{Float64, 1}
    model_optimum::Float64
    model_optimizer::Array{Float64, 1}
    iterations::IterationCounter
    duration::DurationCounter
    opt::NLopt.Opt
    verbosity::Verbosity
    initializer::Ti
    repetitions::Int
    timeroutput::TimerOutput
end
"""
    BOpt(func, model, acquisition, modeloptimizer, lowerbounds, upperbounds;
              sense = Max, maxiterations = 10^4, maxduration = Inf,
              acquisitionoptions = NamedTuple(), repetitions = 1,
              verbosity = Progress,
              initializer_iterations = 5*length(lowerbounds),
              initializer = ScaledSobolIterator(lowerbounds, upperbounds,
                                                initializer_iterations))
"""
function BOpt(func, model, acquisition, modeloptimizer, lowerbounds, upperbounds;
              sense = Max, maxiterations = 10^4, maxduration = Inf,
              acquisitionoptions = NamedTuple(),
              repetitions = 1, verbosity = Progress,
              initializer_iterations = 5*length(lowerbounds),
              initializer = ScaledSobolIterator(lowerbounds, upperbounds,
                                                initializer_iterations))
    now = time()
    acquisitionoptions = merge(defaultoptions(typeof(model), typeof(acquisition)),
               acquisitionoptions)
    maxiterations < length(initializer) && @error("maxiterations = $maxiterations < length(initializer) = $(length(initializer))")

    current_optimum   = isempty(model.y) ? -Inf*Int(sense)           : Int(sense) * maximum(model.y)
    current_optimizer = isempty(model.y) ? zero(float.(lowerbounds)) : model.x[:, argmax(model.y)]
    BOpt(func, sense, model, acquisition,
         acquisitionoptions,
         modeloptimizer, float.(lowerbounds), float.(upperbounds),
         current_optimum, current_optimizer,
         current_optimum, copy(current_optimizer),
         IterationCounter(0, 0, maxiterations),
         DurationCounter(now, maxduration, now, now + maxduration),
         nlopt_setup(acquisition, model, lowerbounds, upperbounds, acquisitionoptions),
         verbosity, initializer, repetitions, TimerOutput())
end
isdone(o::BOpt) = isdone(o.iterations) || isdone(o.duration)
maxduration!(o::BOpt, d) = maxduration!(o.duration, d)
maxiterations!(o::BOpt, N) = maxiterations!(o.iterations, N)
import Base: show
function show(io::IO, mime::MIME"text/plain", o::BOpt)
    println(io, "Bayesian Optimization object\n\nmodel:")
    show(io, mime, o.model)
    println("\nacquisition:")
    show(io, mime, o.acquisition)
    if o.iterations.i == 0
        println("\nNo observation data.")
    else
        println(io, "\n\nobserved optimum: $(o.observed_optimum)")
        println(io, "observed optimizer: $(o.observed_optimizer)")
        println(io, "model optimum: $(o.model_optimum)")
        println(io, "model optimizer: $(o.model_optimizer)")
        println(io, "iterations: $(o.iterations.i)/$(o.iterations.N)")
        println(io, "duration: $(o.duration.now - o.duration.starttime)/$(o.duration.duration) s")
    end
end

function initialise_model!(o)
    ys = Float64[]
    xs = Vector{Float64}[]
    for x in o.initializer
        for j in 1:o.repetitions
            push!(ys, _evaluate_function(o, x))
            push!(xs, x)
        end
    end
    o.iterations.i = o.iterations.c = length(ys)/o.repetitions
    @mytimeit o.timeroutput "model update" update!(o.model, hcat(xs...), ys)
    @mytimeit o.timeroutput "model hyperparameter optimization" optimizemodel!(o.modeloptimizer, o.model)
end
"""
    boptimize!(o::BOpt)
"""
function boptimize!(o::BOpt)
    init!(o.duration)
    init!(o.iterations)
    reset_timer!(o.timeroutput)
    o.iterations.i == 0 && length(o.initializer) > 0 && initialise_model!(o)
    while !isdone(o)
        o.verbosity >= Progress && @info("$(now())\titeration: $(o.iterations.i)\tcurrent optimum: $(o.observed_optimum)")
        setparams!(o.acquisition, o.model)
        @mytimeit o.timeroutput "acquisition" f, x = acquire_max(o.opt, o.lowerbounds, o.upperbounds, o.acquisitionoptions.restarts)
        ys = Float64[]
        step!(o.iterations)
        for _ in 1:o.repetitions
            y = _evaluate_function(o, x)
            push!(ys, y)
        end
        @mytimeit o.timeroutput "model update" update!(o.model, hcat(fill(x, o.repetitions)...), ys)
        @mytimeit o.timeroutput "model hyperparameter optimization" optimizemodel!(o.modeloptimizer, o.model)
    end
     @mytimeit o.timeroutput "acquisition" o.model_optimum, o.model_optimizer = acquire_model_max(o)
    o.duration.now = time()
    o.verbosity >= Timings && @info(o.timeroutput)
    (observed_optimum = o.observed_optimum, observed_optimizer = o.observed_optimizer,
     model_optimum = Int(o.sense) * o.model_optimum, model_optimizer = o.model_optimizer)
end

function _evaluate_function(o, x)
    @mytimeit o.timeroutput "function evaluation" y = Int(o.sense) * o.func(x)
    if y > Int(o.sense) * o.observed_optimum
        o.observed_optimum = Int(o.sense) * y
        o.observed_optimizer = x
    end
    return y
end

end # module


text:
abstract type Block end
id(b::Block)=String(typeof(b).name.name)
text(b::Block)=langs[:name][typeof(b).name.name]
stack(::Block)=0x4
c_use(::Block)=false
use(::Block,::Int,::Int)=nothing
function i_show(i::Block,con::DrawContext)
	if hole(i) clear_rect(con,0,0;color=RGB_(0x70,0x80,0x90)) end
	iid=id(i)
	so=get_loadedimg(iid)
	fill_image(con,so,0,0)
end
function b_show(b::Block,con,x::Int,y::Int)
	bid=id(b)
	so=get_loadedimg(bid)
	fill_image(con,so,x,y)
end
function hard(b::Block)
	bid=id(b)
	if haskey(data[:硬度dict],bid)
		return data[:硬度dict][bid]::UInt8
	end
	return 0x0
end
light(::Block)=0x0
hole(::Block)=false
e_put(::Block,::Int,::Int)=nothing
exist(::Block,::Int,::Int)=nothing
e_exist(::Block,::Int,::Int)=nothing
function li_show(b::Block,x::Int,y::Int,li::UInt8=light(b);arr::SMatrix{31,31,Bool,961})
	if li==0x0 return end
	for i in 1:li<<1|1
		for j in 1:li<<1|1 arr[i,j]=false end
	end
	xq=[0];yq=[0];lq=[li] # 任务队列
	while !isempty(lq)
		xx=popfirst!(xq);yy=popfirst!(yq);l=popfirst!(lq)
		arr[xx+16,yy+16]=true
		tl=getlig(x+xx,y+yy)
		if l<=tl continue end
		setlig(x+xx,y+yy,l)
		if (!transparent(getblk(x+xx,y+yy))&&(xx!=0||yy!=0)) || l==0x1
			continue
		end
		lowl::UInt8=l-0x1
		for i in 1:4
			x0=xx+direx[i];y0=yy+direy[i]
			if !arr[x0+16,y0+16]
				push!(xq,x0)
				push!(yq,y0)
				push!(lq,lowl)
			end
		end
	end
end
upd(::Block)=nothing

abstract type 气体<:Block end
abstract type 液体<:Block end
colbox(::Union{气体,液体},x::Int,y::Int)=EmptyCollisionBox()
hard(::Union{气体,液体})=zero(UInt)
transparent(::Union{气体,液体})=true
struct 虚无<:气体 end
b_show(::虚无,con,x::Int,y::Int)=fill_rect(con,x,y;color=:black)
struct 空气<:气体 end # 进行特殊处理
li_show(::空气,x::Int,y::Int)=nothing
include("水.jl")
include("岩浆.jl")

abstract type 固体<:Block end
colbox(::固体,x::Int,y::Int)=BCollisionBox(x,y,x+1,y+1)
transparent(::固体)=false
struct 星岩<:固体 end
struct 玻璃<:固体 end
transparent(::玻璃)=true
struct 沙子<:固体 end
include("半砖.jl")
include("冰.jl")
include("地毯.jl")
include("高炉壁.jl")
include("高炉控制器.jl")
include("告示牌.jl")
include("合成台.jl")
include("结构方块.jl")
include("刻制台.jl")
include("块.jl")
include("模具台.jl")
include("木板.jl")
include("树苗.jl")
include("树叶.jl")
include("箱子.jl")
include("原木.jl")
include("组装台.jl")
include("Befunge.jl")

abstract type 岩石<:固体 end
struct 石头<:岩石 end
struct 花岗岩<:岩石 end
struct 安山岩<:岩石 end
struct 玄武岩<:岩石 end
struct 闪长岩<:岩石 end
struct 黑曜岩<:岩石 end
struct 萤石<:岩石 end
light(::萤石)=0xc

abstract type 花<:固体 end
colbox(::花)=EmptyCollisionBox()
transparent(::花)=true
hole(::花)=true
struct 番红花<:花 end ; binomial_nomenclature(::番红花)=Pair("Crocus","sativus")