text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

function gnuplot_init()
    global gnuplot_state
    f = ccall(:popen, Ptr{Int}, (Ptr{Uint8},Ptr{Uint8}), "gnuplot" ,"w")
    if f == C_NULL
        error("There was a problem starting up gnuplot.")
    end
    gnuplot_state.running = true
    gnuplot_state.fid = f
end

# close gnuplot pipe
function gnuplot_exit(x...)
    global gnuplot_state
    if gnuplot_state.running
        # close pipe
        err = ccall(:pclose, Int, (Ptr{Int},), gnuplot_state.fid)
        # err should be zero
        if err != 0
            println("Gnuplot may not have closed correctly.");
        end
    end
    # reset gnuplot_state
    gnuplot_state.running = false
    gnuplot_state.current = 0
    gnuplot_state.fid = 0
    gnuplot_state.figs = []
    return 0
end

# Return index to figure with handle 'c'. If no such figure exists, returns 0.
function findfigure(c)
    global gnuplot_state
    i = 0
    for j = 1:length(gnuplot_state.figs)
        if gnuplot_state.figs[j].handle == c
            i = j
            break
        end
    end
    return i
end

# convert marker string description to gnuplot's expected number
function pointtype(x::String)
    if x == "+"
        return 1
    elseif x == "x"
        return 2
    elseif x == "*"
        return 3
    elseif x == "esquare"
        return 4
    elseif x == "fsquare"
        return 5
    elseif x == "ecircle"
        return 6
    elseif x == "fcircle"
        return 7
    elseif x == "etrianup"
        return 8
    elseif x == "ftrianup"
        return 9
    elseif x == "etriandn"
        return 10
    elseif x == "ftriandn"
        return 11
    elseif x == "edmd"
        return 12
    elseif x == "fdmd"
        return 13
    end
    return 1
end

# return configuration string for a single plot
function linestr_single(conf::CurveConf)
    s = ""
    if conf.legend != ""
        s = string(s, " title '", conf.legend, "' ")
    else
        s = string(s, "notitle ")
    end
    s = string(s, " with ", conf.plotstyle, " ")
    if conf.color != ""
        s = string(s, "linecolor rgb '", conf.color, "' ")
    end
    s = string(s, "lw ", string(conf.linewidth), " ")
    # some plotstyles don't allow point specifiers
    cp = conf.plotstyle
    if cp != "lines" && cp != "impulses" && cp != "pm3d" && cp != "image" &&
        cp != "rgbimage" && cp != "boxes" && cp != "dots"
        if conf.marker != ""
            s = string(s, "pt ", string(pointtype(conf.marker)), " ")
        end
        s = string(s, "ps ", string(conf.pointsize), " ")
    end
    return s
end

# build a string with plot commands according to configuration
function linestr(curves::Vector{CurveData}, cmd::String, file::String)
    # We have to insert "," between plot commands. One easy way to do this
    # is create the first plot command, then the rest
    # We also need to keep track of the current index (starts at zero)
    index = 0
    s = string(cmd, " '", file, "' ", " i 0 ", linestr_single(curves[1].conf))
    if length(curves) > 1
        for i in curves[2:end]
            index += 1
            s = string(s, ", '", file, "' ", " i ", string(index)," ",
                linestr_single(i.conf))
        end
    end
    return s
end

# create a Z-coordinate matrix from x, y coordinates and a function
function meshgrid(x,y,f)
    Z = zeros(length(x),length(y))
    for k = 1:length(x)
        Z[k,:] = [ f(i,j) for i=x[k], j=y ]
    end
    return Z
end

# create x,y coordinates for a histogram, from a sample vector, using a number
# of bins
function histdata(s,bins)
    # When adding an element s to a bin, we use an iequality m < s <= M.
    # In order to account for elements s==m, we need to special-case
    # the computation for the first bin
    ms = min(s)
    Ms = max(s)
    bins = max(bins, 1)
    if Ms == ms
        # compute a "natural" scale
        g = (10.0^floor(log10(abs(ms)+eps()))) / 2
        ms, Ms = ms - g, ms + g
    end
    delta = (Ms-ms)/bins
    x = Range(float(ms), delta, bins+1) # like ms:delta:Ms but less error-prone
    y = zeros(bins)

    # this is special-cased because we want to include the minimum in the
    # first bin
    y[1] = sum(ms .<= s .<= x[2])
    for i in 2:length(x)-2
        y[i] = sum(x[i] .< s .<= x[i+1])
    end
    # this is special-cased because there is no guarantee that x[end] == Ms
    # (because of how ranges work)
    if length(y) > 1 y[end] = sum(x[end-1] .< s .<= Ms) end

    if bins != 1
        # We want the left bin to start at ms and the right bin to end at Ms
        x = Range(ms+delta/2, delta, bins)
    else
        # add two empty bins on the sides to provide a scale to gnuplot
        x = Range(ms-delta/2, delta, 3)
        y = [0.0, y[1], 0.0]
    end

    return x,y
end

# dereference CurveConf, by adding a method to copy()
function copy(conf::CurveConf)
    new = CurveConf()
    new.legend = conf.legend
    new.plotstyle = conf.plotstyle
    new.color = conf.color
    new.marker = conf.marker
    new.linewidth = conf.linewidth
    new.pointsize = conf.pointsize
    return new
end

# dereference AxesConf
function copy(conf::AxesConf)
    new = AxesConf()
    new.title = conf.title
    new.xlabel = conf.xlabel
    new.ylabel = conf.ylabel
    new.zlabel = conf.zlabel
    new.box = conf.box
    new.axis = conf.axis
    return new
end

# Build a "set term" string appropriate for the terminal type
function termstring(term::String)
    global gnuplot_state
    global gaston_config

    gc = gaston_config

    if is_term_screen(term)
        ts = "set term $term $(gnuplot_state.current)"
    else
        if term == "pdf"
            s = "set term pdfcairo $(gc.print_color) "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "eps"
            s = "set term epscairo $(gc.print_color) "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "png"
            s = "set term pngcairo $(gc.print_color) "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "gif"
            s = "set term gif "
            s = "$s font $(gc.print_fontface) $(gc.print_fontsize) "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "svg"
            s = "set term svg "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        end
        ts = "$s \nset output \"$(gc.outputfile)\""
    end
    return ts
end

# send gnuplot the current figure's configuration
function gnuplot_send_fig_config(config)
    gnuplot_send("set autoscale")
    # legend box
    if config.box != ""
        gnuplot_send(string("set key ",config.box))
    end
    # plot title
    if config.title != ""
        gnuplot_send(string("set title '",config.title,"' "))
    end
    # xlabel
    if config.xlabel != ""
        gnuplot_send(string("set xlabel '",config.xlabel,"' "))
    end
    # ylabel
    if config.ylabel != ""
        gnuplot_send(string("set ylabel '",config.ylabel,"' "))
    end
    # zlabel
    if config.zlabel != ""
        gnuplot_send(string("set zlabel '",config.zlabel,"' "))
    end
    # axis log scale
    if config.axis != "" || config.axis != "normal"
        if config.axis == "semilogx"
            gnuplot_send("set logscale x")
        end
        if config.axis == "semilogy"
            gnuplot_send("set logscale y")
        end
        if config.axis == "loglog"
            gnuplot_send("set logscale xy")
        end
    end
end

# Validation functions.
# These functions validate that configuration parameters are valid and
# supported. They return true iff the argument validates.

# Validate terminal type.
function validate_terminal(s::String)
    supp_terms = ["wxt", "x11", "svg", "gif", "png", "pdf", "aqua", "eps"]
    if s == "aqua" && OS_NAME != :Darwin
        return false
    end
    if contains(supp_terms, s)
        return true
    end
    return false
end

# Identify terminal by type: file or screen
function is_term_screen(s::String)
    screenterms = ["wxt", "x11", "aqua"]
    if contains(screenterms, s)
        return true
    end
    return false
end

function is_term_file(s::String)
    screenterms = ["svg", "gif", "png", "pdf", "eps"]
    if contains(screenterms, s)
        return true
    end
    return false
end

# Valid plotstyles supported by gnuplot's plot
function validate_2d_plotstyle(s::String)
    valid = ["lines", "linespoints", "points", "impulses", "boxes",
        "errorlines", "errorbars", "dots"]
    if contains(valid, s)
        return true
    end
    return false
end

# Valid plotstyles supported by gnuplot's splot
function validate_3d_plotstyle(s::String)
    valid = ["lines", "linespoints", "points", "impulses", "pm3d",
            "image", "rgbimage", "dots"]
    if contains(valid, s)
        return true
    end
    return false
end

# Valid axis types
function validate_axis(s::String)
    valid = ["", "normal", "semilogx", "semilogy", "loglog"]
    if contains(valid,s)
        return true
    end
    return false
end

# Valid markers supported by Gaston
function validate_marker(s::String)
    valid = ["", "+", "x", "*", "esquare", "fsquare", "ecircle", "fcircle",
    "etrianup", "ftrianup", "etriandn", "ftriandn", "edmd", "fdmd"]
    if contains(valid, s)
        return true
    end
    return false
end


text:
module NNlib

export σ, relu, softmax

include("activation.jl")
include("adapt.jl")

end # module

ext: jl


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 = "XNNPACK"
version = v"0.0.20200323"

cpuinfo_build_version = v"0.0.20200228"

cpuinfo_compat = [
    cpuinfo_build_version,
    v"0.0.20200522", # previously released
    v"0.0.20200612", # due to pytorch v1.6.0 - v1.7.1, e.g. https://github.com/pytorch/pytorch/tree/v1.6.0/third_party/cpuinfo @ 63b254577ed77a8004a9be6ac707f3dccc4e1fd9
]

pthreadpool_build_version = v"0.0.20200302"

pthreadpool_compat = [
    pthreadpool_build_version,
    v"0.0.20200616", # due to pytorch v1.6.0 - v1.7.1, e.g. https://github.com/pytorch/pytorch/tree/v1.6.0/third_party/pthreadpool @ 029c88620802e1361ccf41d1970bd5b07fd6b7bb
]

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/google/XNNPACK.git", "1b354636b5942826547055252f3b359b54acff95"),
    DirectorySource("./bundled"),
    GitSource("https://github.com/pytorch/cpuinfo.git", "d5e37adf1406cf899d7d9ec1d317c47506ccb970"; unpack_target="clog"),
    GitSource("https://github.com/Maratyszcza/FP16.git", "ba1d31f5eed2eb4a69e4dea3870a68c7c95f998f"),
    GitSource("https://github.com/Maratyszcza/FXdiv.git", "f8c5354679ec2597792bc70a9e06eff50c508b9a"),
    GitSource("https://github.com/Maratyszcza/psimd.git", "88882f601f8179e1987b7e7cf4a8012c9080ad44"),
    GitSource("https://github.com/Maratyszcza/pthreadpool.git", "ebd50d0cfa3664d454ffdf246fcd228c3b370a11"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/XNNPACK
atomic_patch -p1 ../patches/xnnpack-disable-fast-math.patch
atomic_patch -p1 ../patches/xnnpack-pic.patch
atomic_patch -p1 ../patches/xnnpack-soversion.patch
atomic_patch -p1 ../patches/xnnpack-freebsd.patch
if [[ $target == aarch64-* ]]; then
    atomic_patch -p1 ../patches/xnnpack-disable-neon-fp16-arithmetic.patch
fi
mkdir build
cd build
# Omitted cmake define of CPUINFO_SOURCE_DIR as there is a patch for cpuinfo
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_INSTALL_LIBDIR=$libdir \
    -DCLOG_SOURCE_DIR=$WORKSPACE/srcdir/clog/cpuinfo \
    -DFP16_SOURCE_DIR=$WORKSPACE/srcdir/FP16 \
    -DFXDIV_SOURCE_DIR=$WORKSPACE/srcdir/FXdiv \
    -DPSIMD_SOURCE_DIR=$WORKSPACE/srcdir/psimd \
    -DPTHREADPOOL_SOURCE_DIR=$WORKSPACE/srcdir/pthreadpool \
    -DXNNPACK_LIBRARY_TYPE=shared \
    -DXNNPACK_BUILD_TESTS=OFF \
    -DXNNPACK_BUILD_BENCHMARKS=OFF \
    ..
if [[ $target == *-w64-mingw32* ]]; then
    cd cpuinfo-source
        atomic_patch -p1 ../../../patches/cpuinfo-mingw-lowercase-windows-include.patch
    cd ..
fi
cmake --build . -- -j $nproc
make install/local
if [[ $target == *-w64-mingw32 ]]; then
    install -Dvm 755 libXNNPACK.dll "${libdir}/libXNNPACK.dll"
fi
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(p -> !(Sys.iswindows(p) && arch(p) == "i686"), platforms) # i686-Windows fails to link
filter!(p -> arch(p) != "powerpc64le", platforms) # PowerPC64LE is unsupported by XNNPACK (Unsupported architecture in src/init.c)

# The products that we will ensure are always built
products = [
    LibraryProduct("libXNNPACK", :libxnnpack),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("CPUInfo_jll", cpuinfo_build_version; compat=join(string.(cpuinfo_compat), ", ")),
    Dependency("PThreadPool_jll", pthreadpool_build_version; compat=join(string.(pthreadpool_compat), ", ")),
]

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

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
immutable HTMLCell <: AbstractCell
    text::Nullable{String}
    colspan::Int
    rowspan::Int
    tagstring::String
end

function showcell(io::IO, c::HTMLCell)
    rows = c.rowspan > 1 ? " rowspan=\"$(c.rowspan)\"" : ""
    cols = c.colspan > 1 ? " colspan=\"$(c.colspan)\"" : ""

    if !isnull(c.text)
        println(io, "\t\t<td$(rows)$(cols)$(c.tagstring)>$(get(c.text))</td>")
    end
end

function tcell(::Type{HTMLCell},
              styleel::AbstractStyleElement,
              datael,
              i::Int, j::Int)
  hascontents = istext(styleel, datael, i, j)
  txt = hascontents ?
    Nullable(htmlescape(text(datael[i,j], styleel.text))) :
    Nullable{eltype(datael)}()
  cs, rs  = colspan(styleel, datael, i, j), rowspan(styleel, datael, i, j)

  return HTMLCell(txt, cs, rs, styleel.tagstring)
end


text:
module TestMul

include("preamble.jl")

@testset begin
    @test mul!!(0, 1, 2) === 2
    @test mul!!(missing, 1, 2) === 2
    @test mul!!(nothing, 1, 2) === 2

    C = zeros(2, 2)
    A = reshape(1:4, 2, 2)
    B = ones(2, 2)
    @test mul!!(C, A, B) === C
    @test C == A * B
    @test_inferred mul!!(C, A, B)

    C = zeros(Int, size(A * B))
    @test mul!!(C, A, B) :: Matrix{Float64} == A * B
    @test_inferred mul!!(C, A, B)
end

if VERSION >= v"1.3"
    @testset "mul!(C, A, B, α, β)" begin
        C = zeros(2, 2)
        A = reshape(1:4, 2, 2)
        B = ones(2, 2)

        desired = A * B * 3 + C * 4
        @test mul!!(C, A, B, 3, 4) === C
        @test C == desired
        @test_inferred mul!!(C, A, B, 3, 4)

        desired = A * B * 3 + C * 4
        @test mul!!(C, A, B, 3, 4) :: Matrix{Float64} == desired
        @test_inferred mul!!(C, A, B, 3, 4)
    end
end

end  # module


text:
function _build_and_constraint(
    _error::Function,
    lhs,
    rhs,
)
    lhs_set = JuMP.moi_set(lhs)
    rhs_set = JuMP.moi_set(rhs)

    lhs_jump_func = JuMP.jump_function(lhs)
    rhs_jump_func = JuMP.jump_function(rhs)

    lhs_func = JuMP.moi_function(lhs)
    rhs_func = JuMP.moi_function(rhs)

    func = [lhs_jump_func..., rhs_jump_func...]
    return JuMP.VectorConstraint(
        func, AndSet{typeof(lhs_func), typeof(rhs_func)}(lhs_set, rhs_set)
    )
end

function JuMP.parse_constraint_head(_error::Function, ::Val{:(&&)}, lhs, rhs)
    _error1 = deepcopy(_error)
    lhs_vectorized, lhs_parsecode, lhs_buildcall =
        JuMP.parse_constraint_expr(_error, lhs)

    if lhs_vectorized
        _error("`$(lhs)` should be non vectorized. There is currently no vectorized support for `and` constraints. Please open an issue at ConstraintSolver.jl")
    end

    rhs_vectorized, rhs_parsecode, rhs_buildcall =
        JuMP.parse_constraint_expr(_error1, rhs)

    if rhs_vectorized
        _error("`$(rhs)` should be non vectorized. There is currently no vectorized support for `and` constraints. Please open an issue at ConstraintSolver.jl")
    end

    # TODO implement vectorized version
    vectorized = false
    complete_parsecode = quote
        $lhs_parsecode
        $rhs_parsecode
    end

    buildcall = :($(esc(:(CS._build_and_constraint)))(
        $_error,
        $lhs_buildcall,
        $rhs_buildcall
    ))
    return vectorized, complete_parsecode, buildcall
end


text:
using AbInitioSoftwareBase.Inputs: Namelist
using Compat: @NamedTuple
using ConstructionBase: setproperties
using StructHelpers: @batteries

const DVSCF_STAR = @NamedTuple begin
    open::Bool
    dir::String
    ext::String
    basis::String
    pat::Bool
end
const DRHO_STAR = @NamedTuple begin
    open::Bool
    dir::String
    ext::String
    basis::String
    pat::Bool
end

# The following default values are picked from `<QE source>/test-suite/not_epw_comp/phq_readin.f90`
"""
    PhNamelist(amass, outdir, prefix, niter_ph, tr2_ph, alpha_mix, nmix_ph, verbosity, reduce_io, max_seconds, fildyn, fildrho, fildvscf, epsil, lrpa, lnoloc, trans, lraman, eth_rps, eth_ns, dek, recover, low_directory_check, only_init, qplot, q2d, q_in_band_form, electron_phonon, lshift_q, zeu, zue, elop, fpol, ldisp, nogg, asr, ldiag, lqdir, search_sym, nq1, nq2, nq3, nk1, nk2, nk3, k1, k2, k3, start_irr, last_irr, nat_todo, modenum, start_q, last_q, dvscf_star, drho_star)
    PhNamelist(; kwargs...)
    PhNamelist(::PhNamelist; kwargs...)

Represent the `INPUTPH` namelist of `ph.x`.
"""
struct PhNamelist <: Namelist
    amass::Vector{Union{Nothing,Float64}}
    outdir::String
    prefix::String
    niter_ph::Int
    tr2_ph::Float64
    alpha_mix::Vector{Union{Nothing,Float64}}
    nmix_ph::Int
    verbosity::String
    reduce_io::Bool
    max_seconds::Float64
    fildyn::String
    fildrho::String
    fildvscf::String
    epsil::Bool
    lrpa::Bool
    lnoloc::Bool
    trans::Bool
    lraman::Bool
    eth_rps::Float64
    eth_ns::Float64
    dek::Float64
    recover::Bool
    low_directory_check::Bool
    only_init::Bool
    qplot::Bool
    q2d::Bool
    q_in_band_form::Bool
    electron_phonon::String
    lshift_q::Bool
    zeu::Bool  # The default value in QE's source code is `true`
    zue::Bool
    elop::Bool
    fpol::Bool
    ldisp::Bool
    nogg::Bool
    asr::Bool
    ldiag::Bool
    lqdir::Bool
    search_sym::Bool
    nq1::Int
    nq2::Int
    nq3::Int
    nk1::Int
    nk2::Int
    nk3::Int
    k1::Int
    k2::Int
    k3::Int
    start_irr::Int  # The default value in QE's source code is 0
    last_irr::Int
    nat_todo::Int
    modenum::Int
    start_q::Int
    last_q::Int
    dvscf_star::DVSCF_STAR
    drho_star::DRHO_STAR
end # struct PhNamelist
function PhNamelist(;
    amass = [0.0],
    outdir = "./",
    prefix = "pwscf",
    niter_ph = 100,
    tr2_ph = 1e-12,
    alpha_mix = 0.7 * ones(niter_ph),
    nmix_ph = 4,
    verbosity = "default",
    reduce_io = false,
    max_seconds = 1e7,
    fildyn = "matdyn",
    fildrho = " ",
    fildvscf = " ",
    epsil = false,
    lrpa = false,
    lnoloc = false,
    trans = true,
    lraman = false,
    eth_rps = 1e-9,
    eth_ns = 1e-12,
    dek = 1e-3,
    recover = false,
    low_directory_check = false,
    only_init = false,
    qplot = false,
    q2d = false,
    q_in_band_form = false,
    electron_phonon = " ",
    lshift_q = false,
    zeu = epsil,  # The default value in QE's source code is `true`
    zue = false,
    elop = false,
    fpol = false,
    ldisp = false,
    nogg = false,
    asr = false,
    ldiag = false,
    lqdir = false,
    search_sym = true,
    nq1 = 0,
    nq2 = 0,
    nq3 = 0,
    nk1 = 0,
    nk2 = 0,
    nk3 = 0,
    k1 = 0,
    k2 = 0,
    k3 = 0,
    start_irr = 1,  # The default value in QE's source code is 0
    last_irr = -1000,
    nat_todo = 0,
    modenum = 0,
    start_q = 1,
    last_q = -1000,
    dvscf_star = (
        open = false,
        dir = "./",
        ext = "dvscf",
        basis = "cartesian",
        pat = false,
    ),
    drho_star = (open = false, dir = "./", ext = "drho", basis = "modes", pat = true),
)
    return PhNamelist(
        amass,
        outdir,
        prefix,
        niter_ph,
        tr2_ph,
        alpha_mix,
        nmix_ph,
        verbosity,
        reduce_io,
        max_seconds,
        fildyn,
        fildrho,
        fildvscf,
        epsil,
        lrpa,
        lnoloc,
        trans,
        lraman,
        eth_rps,
        eth_ns,
        dek,
        recover,
        low_directory_check,
        only_init,
        qplot,
        q2d,
        q_in_band_form,
        electron_phonon,
        lshift_q,
        zeu,
        zue,
        elop,
        fpol,
        ldisp,
        nogg,
        asr,
        ldiag,
        lqdir,
        search_sym,
        nq1,
        nq2,
        nq3,
        nk1,
        nk2,
        nk3,
        k1,
        k2,
        k3,
        start_irr,
        last_irr,
        nat_todo,
        modenum,
        start_q,
        last_q,
        dvscf_star,
        drho_star,
    )
end
PhNamelist(nml::PhNamelist; kwargs...) = setproperties(nml, kwargs...)
@batteries PhNamelist eq = true hash = true

function (x::VerbositySetter)(control::PhNamelist)
    @set! control.verbosity = x.v
    return control
end

# The following default values are picked from `<QE source>/PHonon/PH/q2r.f90`
"""
    Q2rNamelist(fildyn, flfrc, loto_2d, zasr)
    Q2rNamelist(; kwargs...)
    Q2rNamelist(::Q2rNamelist; kwargs...)

Represent the `INPUT` namelist of `q2r.x`.
"""
struct Q2rNamelist <: Namelist
    fildyn::String
    flfrc::String
    loto_2d::Bool
    zasr::String
end # struct Q2rNamelist
function Q2rNamelist(; fildyn = " ", flfrc = " ", loto_2d = false, zasr = "no")
    return Q2rNamelist(fildyn, flfrc, loto_2d, zasr)
end
Q2rNamelist(nml::Q2rNamelist; kwargs...) = setproperties(nml, kwargs...)
@batteries Q2rNamelist eq = true hash = true

# The following default values are picked from `<QE source>/PHonon/PH/matdyn.f90`
"""
    MatdynNamelist(dos, deltaE, ndos, nk1, nk2, nk3, asr, readtau, flfrc, fldos, flfrq, flvec, fleig, fldyn, fltau, amass, at, ntyp, l1, l2, l3, la2F, q_in_band_form, eigen_similarity, q_in_cryst_coord, na_ifc, fd, nosym, loto_2d)
    MatdynNamelist(; kwargs...)
    MatdynNamelist(::MatdynNamelist; kwargs...)

Represent the `INPUT` namelist of `matdyn.x`.
"""
struct MatdynNamelist <: Namelist
    dos::Bool
    deltaE::Float64
    ndos::Int
    nk1::Int
    nk2::Int
    nk3::Int
    asr::String
    readtau::Bool
    flfrc::String
    fldos::String
    flfrq::String
    flvec::String
    fleig::String
    fldyn::String
    fltau::String
    amass::Vector{Union{Nothing,Float64}}
    at::Matrix{Union{Nothing,Float64}}  # FIXME: not very sure
    ntyp::Int
    l1::Int
    l2::Int
    l3::Int
    la2F::Bool
    q_in_band_form::Bool
    eigen_similarity::Bool
    q_in_cryst_coord::Bool
    na_ifc::Bool
    fd::Bool
    nosym::Bool
    loto_2d::Bool
end # struct MatdynNamelist
function MatdynNamelist(;
    dos = false,
    deltaE = 1.0,
    ndos = 1,
    nk1 = 0,
    nk2 = 0,
    nk3 = 0,
    asr = "no",
    readtau = false,
    flfrc = " ",
    fldos = "matdyn.dos",
    flfrq = "matdyn.freq",
    flvec = "matdyn.modes",
    fleig = "matdyn.eig",
    fldyn = " ",
    fltau = " ",
    amass = zeros(1),
    at = zeros(3, 3),  # FIXME: not very sure
    ntyp = 0,
    l1 = 1,
    l2 = 1,
    l3 = 1,
    la2F = false,
    q_in_band_form = false,
    eigen_similarity = false,
    q_in_cryst_coord = false,
    na_ifc = false,
    fd = false,
    nosym = false,
    loto_2d = false,
)
    return MatdynNamelist(
        dos,
        deltaE,
        ndos,
        nk1,
        nk2,
        nk3,
        asr,
        readtau,
        flfrc,
        fldos,
        flfrq,
        flvec,
        fleig,
        fldyn,
        fltau,
        amass,
        at,
        ntyp,
        l1,
        l2,
        l3,
        la2F,
        q_in_band_form,
        eigen_similarity,
        q_in_cryst_coord,
        na_ifc,
        fd,
        nosym,
        loto_2d,
    )
end
MatdynNamelist(nml::MatdynNamelist; kwargs...) = setproperties(nml, kwargs...)
@batteries MatdynNamelist eq = true hash = true

"""
    DynmatNamelist(asr, axis, fildyn, filout, filmol, filxsf, fileig, amass, q, lperm, lplasma)
    DynmatNamelist(; kwargs...)
    DynmatNamelist(::DynmatNamelist; kwargs...)

Represent the `INPUT` namelist of `dynmat.x`.
"""
struct DynmatNamelist <: Namelist
    asr::String
    axis::Int
    fildyn::String
    filout::String
    filmol::String
    filxsf::String
    fileig::String
    amass::Vector{Union{Nothing,Float64}}
    q::Vector{Union{Nothing,Float64}}
    lperm::Bool
    lplasma::Bool
end # struct DynmatNamelist
function DynmatNamelist(;
    asr = "no",
    axis = 3,
    fildyn = "matdyn",
    filout = "dynmat.out",
    filmol = "dynmat.mold",
    filxsf = "dynmat.axsf",
    fileig = " ",
    amass = zeros(1),
    q = zeros(3),
    lperm = false,
    lplasma = false,
)
    return DynmatNamelist(
        asr,
        axis,
        fildyn,
        filout,
        filmol,
        filxsf,
        fileig,
        amass,
        q,
        lperm,
        lplasma,
    )
end
DynmatNamelist(nml::DynmatNamelist; kwargs...) = setproperties(nml, kwargs...)
@batteries DynmatNamelist eq = true hash = true

groupname(::Type{PhNamelist}) = "INPUTPH"
groupname(::Type{Q2rNamelist}) = "INPUT"
groupname(::Type{MatdynNamelist}) = "INPUT"
groupname(::Type{DynmatNamelist}) = "INPUT"


text:
export crc, spec, CHECK, NoTables, Single, Multiple, ALL,
       CRC_3_ROHC, CRC_4_ITU, CRC_5_EPC, CRC_5_ITU, CRC_5_USB,
       CRC_6_CDMA2000_A, CRC_6_CDMA2000_B, CRC_6_DARC, CRC_6_ITU,
       CRC_7, CRC_7_ROHC, CRC_8, CRC_8_CDMA2000, CRC_8_DARC,
       CRC_8_DVB_S2, CRC_8_EBU, CRC_8_I_CODE, CRC_8_ITU, CRC_8_MAXIM,
       CRC_8_ROHC, CRC_8_WCDMA, CRC_10, CRC_10_CDMA2000, CRC_11,
       CRC_12_3GPP, CRC_12_CDMA2000, CRC_12_DECT, CRC_13_BBC,
       CRC_14_DARC, CRC_15, CRC_15_MPT1327, CRC_16_ARC,
       CRC_16_AUG_CCITT, CRC_16_BUYPASS, CRC_16_CCITT_FALSE,
       CRC_16_CDMA2000, CRC_16_DDS_110, CRC_16_DECT_R, CRC_16_DECT_X,
       CRC_16_DNP, CRC_16_EN_13757, CRC_16_GENIBUS, CRC_16_MAXIM,
       CRC_16_RIELLO, CRC_16_TELEDISK, CRC_16_USB, CRC_16_CRC_A,
       CRC_16_KERMIT, CRC_16_MODBUS, CRC_16_X_25, CRC_16_XMODEM,
       CRC_24, CRC_24_FLEXRAY_A, CRC_24_FLEXRAY_B, CRC_31_PHILIPS,
       CRC_32, CRC_32_BZIP2, CRC_32_C, CRC_32_D, CRC_32_MPEG_2,
       CRC_32_POSIX, CRC_32_Q, CRC_32_JAMCRC, CRC_32_XFER, CRC_40_GSM,
       CRC_64, CRC_64_WE, CRC_64_XZ, CRC_82_DARC

import Base.Cartesian: @nexprs
import Base: ==, isless, print

const U = Unsigned

struct UnalignedVector{T<:U}
    a::Vector{UInt8}
end
# unsafe
@inline Base.getindex(a::UnalignedVector{T}, i) where {T} = unsafe_load(Ptr{T}(pointer(a.a)), i)
@inline Base.length(a::UnalignedVector{T}) where {T} = length(a.a) ÷ sizeof(T)
@inline unaligned_reinterpret(::Type{T}, a::Vector{UInt8}) where {T<:U} = UnalignedVector{T}(a)
@inline unaligned_reinterpret(::Type{UInt8}, a::Vector{UInt8}) = a
MaybeUnalignedVector{T} = Union{UnalignedVector{T},Vector{T}}

# ---- CRC specifications from http://www.zlib.net/crc_v3.txt

CHECK = Vector{UInt8}(b"123456789")   # universal check vector

struct Spec{P<:U}
    width::Int    # polynomial degree
    poly::P       # generating poly with msb missing
    init::P       # initial remainder
    refin::Bool   # reflect input
    refout::Bool  # reflect output
    xorout::P     # xored with final remainder
    check::P      # checksum for CHECK
    function (::Type{Spec{P}})(width::Int, poly::P, init::P, refin::Bool, refout::Bool, xorout::P, check::P) where {P}
        @assert width <= 8 * sizeof(P)
        new{P}(width, poly, init, refin, refout, xorout, check)
    end
end

spec(poly::P, init::P, refin::Bool, refout::Bool, xorout::P, check::P) where {P<:U} =
    Spec{P}(8*sizeof(P), poly, init, refin, refout, xorout, check)

spec(width::Int, poly::P, init::P, refin::Bool, refout::Bool, xorout::P, check::P) where {P<:U} =
    Spec{P}(width, poly, init, refin, refout, xorout, check)

==(a::Spec, b::Spec) = a.width == b.width && a.poly == b.poly && a.init == b.init && a.refin == b.refin && a.refout == b.refout && a.xorout == b.xorout

isless(a::Spec, b::Spec) = a.width < b.width || (a.width == b.width && (a.poly < b.poly || a.poly == b.poly && (a.init < b.init || (a.init == b.init && !a.refin && (a.refin != b.refin || (a.refin == b.refin && !a.refout && (a.refout != b.refout || a.xorout < b.xorout)))))))

function print(io::IO, s::Spec{P}) where {P<:U}
    n = 1 + div(s.width - 1, 4)
    h = x -> "0x" * string(x, base = 16, pad = n)
    print(io, "width=$(s.width) poly=$(h(s.poly)) init=$(h(s.init)) refin=$(s.refin) refout=$(s.refout) xorout=$(h(s.xorout)) check=$(h(s.check))")
end

# http://reveng.sourceforge.net/crc-catalogue/1-15.htm
CRC_3_ROHC =         spec(3, 0x03, 0x07, true,  true,  0x00, 0x06)
CRC_4_ITU =          spec(4, 0x03, 0x00, true,  true,  0x00, 0x07)
CRC_5_EPC =          spec(5, 0x09, 0x09, false, false, 0x00, 0x00)
CRC_5_ITU =          spec(5, 0x15, 0x00, true,  true,  0x00, 0x07)
CRC_5_USB =          spec(5, 0x05, 0x1f, true,  true,  0x1f, 0x19)
CRC_6_CDMA2000_A =   spec(6, 0x27, 0x3f, false, false, 0x00, 0x0d)
CRC_6_CDMA2000_B =   spec(6, 0x07, 0x3f, false, false, 0x00, 0x3b)
CRC_6_DARC =         spec(6, 0x19, 0x00, true,  true,  0x00, 0x26)
CRC_6_ITU =          spec(6, 0x03, 0x00, true,  true,  0x00, 0x06)
CRC_7 =              spec(7, 0x09, 0x00, false, false, 0x00, 0x75)
CRC_7_ROHC =         spec(7, 0x4f, 0x7f, true,  true,  0x00, 0x53)
CRC_8 =              spec(0x07, 0x00, false, false, 0x00, 0xf4)
CRC_8_CDMA2000 =     spec(0x9b, 0xff, false, false, 0x00, 0xda)
CRC_8_DARC =         spec(0x39, 0x00, true,  true,  0x00, 0x15)
CRC_8_DVB_S2 =       spec(0xd5, 0x00, false, false, 0x00, 0xbc)
CRC_8_EBU =          spec(0x1d, 0xff, true,  true,  0x00, 0x97)
CRC_8_I_CODE =       spec(0x1d, 0xfd, false, false, 0x00, 0x7e)
CRC_8_ITU =          spec(0x07, 0x00, false, false, 0x55, 0xa1)
CRC_8_MAXIM =        spec(0x31, 0x00, true,  true,  0x00, 0xa1)
CRC_8_ROHC =         spec(0x07, 0xff, true,  true,  0x00, 0xd0)
CRC_8_WCDMA =        spec(0x9b, 0x00, true,  true,  0x00, 0x25)
CRC_10 =             spec(10, 0x0233, 0x0000, false, false, 0x0000, 0x0199)
CRC_10_CDMA2000 =    spec(10, 0x03d9, 0x03ff, false, false, 0x0000, 0x0233)
CRC_11 =             spec(11, 0x0385, 0x001a, false, false, 0x0000, 0x05a3)
CRC_12_3GPP =        spec(12, 0x080f, 0x0000, false, true,  0x0000, 0x0daf)
CRC_12_CDMA2000 =    spec(12, 0x0f13, 0x0fff, false, false, 0x0000, 0x0d4d)
CRC_12_DECT =        spec(12, 0x080f, 0x0000, false, false, 0x0000, 0x0f5b)
CRC_13_BBC =         spec(13, 0x1cf5, 0x0000, false, false, 0x0000, 0x04fa)
CRC_14_DARC =        spec(14, 0x0805, 0x0000, true,  true,  0x0000, 0x082d)
CRC_15 =             spec(15, 0x4599, 0x0000, false, false, 0x0000, 0x059e)
CRC_15_MPT1327 =     spec(15, 0x6815, 0x0000, false, false, 0x0001, 0x2566)

# http://reveng.sourceforge.net/crc-catalogue/16.htm#crc.cat.crc-16-ccitt-false
CRC_16_ARC =         spec(0x8005, 0x0000, true,  true,  0x0000, 0xbb3d)
CRC_16_AUG_CCITT =   spec(0x1021, 0x1d0f, false, false, 0x0000, 0xe5cc)
CRC_16_BUYPASS   =   spec(0x8005, 0x0000, false, false, 0x0000, 0xfee8)
CRC_16_CCITT_FALSE = spec(0x1021, 0xffff, false, false, 0x0000, 0x29b1)
CRC_16_CDMA2000 =    spec(0xc867, 0xffff, false, false, 0x0000, 0x4c06)
CRC_16_DDS_110 =     spec(0x8005, 0x800d, false, false, 0x0000, 0x9ecf)
CRC_16_DECT_R =      spec(0x0589, 0x0000, false, false, 0x0001, 0x007e)
CRC_16_DECT_X =      spec(0x0589, 0x0000, false, false, 0x0000, 0x007f)
CRC_16_DNP =         spec(0x3d65, 0x0000, true,  true,  0xffff, 0xea82)
CRC_16_EN_13757 =    spec(0x3d65, 0x0000, false, false, 0xffff, 0xc2b7)
CRC_16_GENIBUS =     spec(0x1021, 0xffff, false, false, 0xffff, 0xd64e)
CRC_16_MAXIM =       spec(0x8005, 0x0000, true,  true,  0xffff, 0x44c2)
CRC_16_RIELLO =      spec(0x8bb7, 0x0000, false, false, 0x0000, 0xd0db)
CRC_16_TELEDISK =    spec(0x1021, 0x89ec, true,  true,  0x0000, 0x26b1)
CRC_16_USB =         spec(0x8005, 0xffff, true,  true,  0xffff, 0xb4c8)
CRC_16_CRC_A =       spec(0x1021, 0xc6c6, true,  true,  0x0000, 0xbf05)
CRC_16_KERMIT =      spec(0x1021, 0x0000, true,  true,  0x0000, 0x2189)
CRC_16_MODBUS =      spec(0x8005, 0xffff, true,  true,  0x0000, 0x4b37)
CRC_16_X_25 =        spec(0x1021, 0xffff, true,  true,  0xffff, 0x906e)
CRC_16_XMODEM =      spec(0x1021, 0x0000, false, false, 0x0000, 0x31c3)

# http://reveng.sourceforge.net/crc-catalogue/17plus.htm
CRC_24 =             spec(24, 0x00864cfb, 0x00b704ce, false, false, 0x00000000, 0x0021cf02)
CRC_24_FLEXRAY_A =   spec(24, 0x005d6dcb, 0x00fedcba, false, false, 0x00000000, 0x007979bd)
CRC_24_FLEXRAY_B =   spec(24, 0x005d6dcb, 0x00abcdef, false, false, 0x00000000, 0x001f23b8)
CRC_31_PHILIPS =     spec(31, 0x04c11db7, 0x7fffffff, false, false, 0x7fffffff, 0x0ce9e46c)
CRC_32 =             spec(0x04c11db7, 0xffffffff, true,  true,  0xffffffff, 0xcbf43926)
CRC_32_BZIP2 =       spec(0x04c11db7, 0xffffffff, false, false, 0xffffffff, 0xfc891918)
CRC_32_C =           spec(0x1edc6f41, 0xffffffff, true,  true,  0xffffffff, 0xe3069283)
CRC_32_D =           spec(0xa833982b, 0xffffffff, true,  true,  0xffffffff, 0x87315576)
CRC_32_MPEG_2 =      spec(0x04c11db7, 0xffffffff, false, false, 0x00000000, 0x0376e6e7)
CRC_32_POSIX =       spec(0x04c11db7, 0x00000000, false, false, 0xffffffff, 0x765e7680)
CRC_32_Q =           spec(0x814141ab, 0x00000000, false, false, 0x00000000, 0x3010bf7f)
CRC_32_JAMCRC =      spec(0x04c11db7, 0xffffffff, true,  true,  0x00000000, 0x340bc6d9)
CRC_32_XFER =        spec(0x000000af, 0x00000000, false, false, 0x00000000, 0xbd0be338)
CRC_40_GSM =         spec(40, 0x0000000004820009, 0x0000000000000000, false, false, 0x000000ffffffffff, 0x000000d4164fc646)
CRC_64 =             spec(0x42f0e1eba9ea3693, 0x0000000000000000, false, false, 0x0000000000000000, 0x6c40df5f0b497347)
CRC_64_WE =          spec(0x42f0e1eba9ea3693, 0xffffffffffffffff, false, false, 0xffffffffffffffff, 0x62ec59e3f1a4f00a)
CRC_64_XZ =          spec(0x42f0e1eba9ea3693, 0xffffffffffffffff, true,  true,  0xffffffffffffffff, 0x995dc9bbdf1939fa)
CRC_82_DARC =        spec(82, 0x0308c0111011401440411, 0x000000000000000000000, true,  true,   0x000000000000000000000, 0x09ea83f625023801fd612)

const ALL = Dict{Symbol,Spec}(:CRC_3_ROHC=>CRC_3_ROHC,
                            :CRC_4_ITU=>CRC_4_ITU,
                            :CRC_5_EPC=>CRC_5_EPC,
                            :CRC_5_ITU=>CRC_5_ITU,
                            :CRC_5_USB=>CRC_5_USB,
                            :CRC_6_CDMA2000_A=>CRC_6_CDMA2000_A,
                            :CRC_6_CDMA2000_B=>CRC_6_CDMA2000_B,
                            :CRC_6_DARC=>CRC_6_DARC,
                            :CRC_6_ITU=>CRC_6_ITU, :CRC_7=>CRC_7,
                            :CRC_7_ROHC=>CRC_7_ROHC, :CRC_8=>CRC_8,
                            :CRC_8_CDMA2000=>CRC_8_CDMA2000,
                            :CRC_8_DARC=>CRC_8_DARC,
                            :CRC_8_DVB_S2=>CRC_8_DVB_S2,
                            :CRC_8_EBU=>CRC_8_EBU,
                            :CRC_8_I_CODE=>CRC_8_I_CODE,
                            :CRC_8_ITU=>CRC_8_ITU,
                            :CRC_8_MAXIM=>CRC_8_MAXIM,
                            :CRC_8_ROHC=>CRC_8_ROHC,
                            :CRC_8_WCDMA=>CRC_8_WCDMA,
                            :CRC_10=>CRC_10,
                            :CRC_10_CDMA2000=>CRC_10_CDMA2000,
                            :CRC_11=>CRC_11,
                            :CRC_12_3GPP=>CRC_12_3GPP,
                            :CRC_12_CDMA2000=>CRC_12_CDMA2000,
                            :CRC_12_DECT=>CRC_12_DECT,
                            :CRC_13_BBC=>CRC_13_BBC,
                            :CRC_14_DARC=>CRC_14_DARC,
                            :CRC_15=>CRC_15,
                            :CRC_15_MPT1327=>CRC_15_MPT1327,
                            :CRC_16_ARC=>CRC_16_ARC,
                            :CRC_16_AUG_CCITT=>CRC_16_AUG_CCITT,
                            :CRC_16_BUYPASS=>CRC_16_BUYPASS,
                            :CRC_16_CCITT_FALSE=>CRC_16_CCITT_FALSE,
                            :CRC_16_CDMA2000=>CRC_16_CDMA2000,
                            :CRC_16_DDS_110=>CRC_16_DDS_110,
                            :CRC_16_DECT_R=>CRC_16_DECT_R,
                            :CRC_16_DECT_X=>CRC_16_DECT_X,
                            :CRC_16_DNP=>CRC_16_DNP,
                            :CRC_16_EN_13757=>CRC_16_EN_13757,
                            :CRC_16_GENIBUS=>CRC_16_GENIBUS,
                            :CRC_16_MAXIM=>CRC_16_MAXIM,
                            :CRC_16_RIELLO=>CRC_16_RIELLO,
                            :CRC_16_TELEDISK=>CRC_16_TELEDISK,
                            :CRC_16_USB=>CRC_16_USB,
                            :CRC_16_CRC_A=>CRC_16_CRC_A,
                            :CRC_16_KERMIT=>CRC_16_KERMIT,
                            :CRC_16_MODBUS=>CRC_16_MODBUS,
                            :CRC_16_X_25=>CRC_16_X_25,
                            :CRC_16_XMODEM=>CRC_16_XMODEM,
                            :CRC_24=>CRC_24,
                            :CRC_24_FLEXRAY_A=>CRC_24_FLEXRAY_A,
                            :CRC_24_FLEXRAY_B=>CRC_24_FLEXRAY_B,
                            :CRC_31_PHILIPS=>CRC_31_PHILIPS,
                            :CRC_32=>CRC_32,
                            :CRC_32_BZIP2=>CRC_32_BZIP2,
                            :CRC_32_C=>CRC_32_C, :CRC_32_D=>CRC_32_D,
                            :CRC_32_MPEG_2=>CRC_32_MPEG_2,
                            :CRC_32_POSIX=>CRC_32_POSIX,
                            :CRC_32_Q=>CRC_32_Q,
                            :CRC_32_JAMCRC=>CRC_32_JAMCRC,
                            :CRC_32_XFER=>CRC_32_XFER,
                            :CRC_40_GSM=>CRC_40_GSM, :CRC_64=>CRC_64,
                            :CRC_64_WE=>CRC_64_WE,
                            :CRC_64_XZ=>CRC_64_XZ,
                            :CRC_82_DARC=>CRC_82_DARC)

# --- Utilities

# http://stackoverflow.com/questions/2602823/in-c-c-whats-the-simplest-way-to-reverse-the-order-of-bits-in-a-byte
function reflect_bits(n::UInt8)
    n = (n & 0xf0) >>> 4 | (n & 0x0f) << 4;
    n = (n & 0xcc) >>> 2 | (n & 0x33) << 2;
    (n & 0xaa) >>> 1 | (n & 0x55) << 1;
end

const REFLECT_8 = UInt8[reflect_bits(i) for i in 0x00:0xff]

reflect(u::UInt8) = REFLECT_8[u+1]

for (T,S) in ((UInt16, UInt8), (UInt32, UInt16), (UInt64, UInt32), (UInt128, UInt64))
    n = 8 * sizeof(S)
    mask::S = -one(S)
    @eval reflect(u::$T) = (convert($T, reflect(convert($S, u & $mask))) << $n) | reflect(convert($S, (u >>> $n) & $mask))
end

function reflect(size, u::T) where {T<:U}
    width = 8 * sizeof(T)
    @assert size <= width "cannot reflect a value larger than the representation"
    reflect(u) >>> (width - size)
end

function largest(T, TS...)
    for t in TS
        if sizeof(t) > sizeof(T)
            T = t
        end
    end
    T
end

function fastest(T, TS...)
    L = largest(T, TS...)
    if L == UInt32 && UInt32 != UInt
        UInt64  # round up for speed
    else
        L
    end
end

function pad(::Type{A}, width) where {A<:U}
    8 * sizeof(A) - width
end

# --- CRC calculations

# a calculation may be a "natural" long division (padded to the size
# of the accumulator) or "reverse the rest of the world".  this
# depends on whether the input input bytes are taken as given
# (Forwards) or reflected (Backwards).

abstract type Direction{A<:U} end

struct Backwards{A<:U}<:Direction{A}
    poly::A
    init::A
    function (::Type{Backwards{A}})(spec::Spec) where {A}
        poly = reflect(spec.width, convert(A, spec.poly))
        init = reflect(spec.width, convert(A, spec.init))
        new{A}(poly, init)
    end
end

struct Forwards{A<:U}<:Direction{A}
    pad_p::Int
    poly::A
    carry::A
    pad_8::Int
    init::A
    function (::Type{Forwards{A}})(spec::Spec) where {A}
        pad_p = pad(A, spec.width)
        poly = convert(A, spec.poly) << pad_p
        carry = one(A) << pad(A, 1)
        pad_8 = pad(A, 8)
        init = convert(A, spec.init) << pad_p
        new{A}(pad_p, poly, carry, pad_8, init)
    end
end

# a calculation may use no, one, or many tables...

abstract type Tables{A<:U} end

struct NoTables{A<:U}<:Tables{A}
    (::Type{NoTables{A}})(direcn::Direction{A}) where {A} = new{A}()
end

struct Multiple{A<:U}<:Tables{A}
    tables::Vector{Vector{A}}
    function (::Type{Multiple{A}})(direcn::Direction{A}) where {A}
        tables = Vector{A}[Array{A}(undef, 256) for _ in 1:sizeof(A)]
        fill_table(direcn, tables[1])
        for index in zero(UInt8):convert(UInt8, 255)
            remainder = tables[1][index + 1]
            for t in 2:sizeof(A)
                # chaining gives the contribution of a byte to the
                # remainder after being shifted through t other zero bytes
                # (the different bytes are then combined with xor).
                remainder = chain(direcn, tables[1], remainder)
                tables[t][index + 1] = remainder
            end
        end
        new{A}(tables)
    end
end

struct Single{A<:U}<:Tables{A}
    table::Vector{A}
    # when using multiple tables we need a related single table for the
    # "last few bytes"
    (::Type{Single{A}})(tables::Multiple{A}) where {A} = new{A}(tables.tables[1])
    (::Type{Single{A}})(direcn::Direction{A}) where {A} = new{A}(fill_table(direcn, Array{A}(undef, 256)))
end

# the basic lookup table is just the remainder for each value
function fill_table(direcn, table::Vector{A}) where {A<:U}
    for index in zero(UInt8):convert(UInt8, 255)
        table[index + 1] = extend(direcn, NoTables{A}(direcn), [index], zero(A))
    end
    table
end

function chain(direcn::Forwards{A}, table::Vector{A}, remainder::A) where {A<:U}
    xor(remainder << 8, table[1 + ((remainder >>> direcn.pad_8) & 0xff)])
end

function chain(direcn::Backwards{A}, table::Vector{A}, remainder::A) where {A<:U}
    xor(remainder >>> 8, table[1 + (remainder & 0xff)])
end

"""
    crc(spec[, tables = Multiple])

Create a crc function implementing `spec` with lookup table behaviour `tables`.
Available specs are listed in [`ALL`](@ref) and tables can be either
[`Multiple`](@ref), [`Single`](@ref) or [`NoTables`](@ref).

The returned function accepts `String`, `Vector{UInt8}` or
`IO` and calculates the CRC checksum. An optional boolean parameter
`append` can be used for chained application.

A faster CRC32*C* implementation is available in Stdlib: [`CRC32c`](@ref).

# Examples
```julia-repl
using CRC

# create our own crc function, just once
crc32 = crc(CRC_32)

# use the crc function created above, many times
for s in ["hello", "there"]
    println(s => crc32(s))
end
```
"""
function crc(spec::Spec{P}; tables::Type{T}=Multiple) where {P<:U, T<:Tables}
    A = fastest(P, tables == Multiple ? UInt : UInt8)
    direcn = spec.refin ? Backwards{A}(spec) : Forwards{A}(spec)
    remainder::A = direcn.init
    tables = tables{A}(direcn)
    function handler(data::Vector{UInt8}; append=false)
        remainder = append ? remainder : direcn.init
        remainder = extend(direcn, tables, data, remainder)
        finalize(spec, direcn, remainder)
    end
    function handler(io::IO; append=false, buflen=1000000)
        buffer = Array{UInt8}(undef, buflen)
        remainder = append ? remainder : direcn.init
        while (nb = readbytes!(io, buffer)) > 0
            remainder = extend(direcn, tables, buffer[1:nb], remainder)
        end
        finalize(spec, direcn, remainder)
    end
    function handler(data::AbstractString; append=false)
        handler(convert(Vector{UInt8}, codeunits(data)), append=append)
    end
    handler
end

# calculations assume that the remainder is either reflected or
# padded; these functions correct for this to return the final result.

function finalize(spec::Spec{P}, direcn::Forwards{A}, remainder::A) where {P<:U, A<:U}
    remainder = convert(P, remainder >> direcn.pad_p)
    remainder = spec.refout ? reflect(spec.width, remainder) : remainder
    xor(remainder, spec.xorout)
end

function finalize(spec::Spec{P}, direcn::Backwards{A}, remainder::A) where {P<:U, A<:U}
    remainder = spec.refout ? remainder : reflect(spec.width, remainder)
    xor(convert(P, remainder), spec.xorout)
end

# all the routines below calculate a new value of remainder for the
# data given (along with the algorithm type, available lookup tables,
# etc).

const UNROLL = 16  # only get small improvements past this

function extend(direcn::Forwards{A}, tables::NoTables, data::Vector{UInt8}, remainder::A) where {A<:U}
    for word::UInt8 in data
        remainder::A = xor(remainder, (convert(A, word) << direcn.pad_8))
        for _ in 1:8
            if remainder & direcn.carry == direcn.carry
                remainder = xor(remainder << 1, direcn.poly)
            else
                remainder <<= 1
            end
        end
    end
    remainder
end

function extend(direcn::Forwards{A}, tables::Single{A}, data::Vector{UInt8}, remainder::A) where {A<:U}
    for i in 1:length(data)
        @inbounds word::UInt8 = data[i]
        remainder::A = xor(remainder, convert(A, word) << direcn.pad_8)
        remainder = xor(remainder << 8, tables.table[1 + remainder >>> direcn.pad_8])
    end
    remainder
end

for A in (UInt16, UInt32, UInt64, UInt128)
    n_tables = sizeof(A)
    @eval begin
        function extend(direcn::Forwards{$A}, tables::Multiple{$A}, data::MaybeUnalignedVector{$A}, remainder::$A)
            word::$A, tmp::$A, remainder::$A = zero($A), zero($A), remainder
            i = 1
            while true
                @nexprs $UNROLL _ -> begin  # unroll inner loop
                    if i > length(data)
                        break
                    end
                    @inbounds word = data[i]
                    tmp, remainder = remainder, zero($A)
                    @nexprs $n_tables t -> begin  # unroll table
                        # access similar to reflected but we need to
                        # handle little-endian bytes which turn out
                        # wrong for this case
                        remainder = xor(remainder,
                          tables.tables[t][xor(word >>> (($n_tables - t)*8), tmp >>> ((t-1)*8)) & 0xff + 1])
                    end
                    i += 1
                end
            end
            remainder
        end
    end
end

function extend(direcn::Backwards{A}, tables::NoTables, data::Vector{UInt8}, remainder::A) where {A<:U}
    for word::UInt8 in data
        remainder::A = xor(remainder, convert(A, word))
        for _ in 1:8
            if remainder & one(A) == one(A)
                remainder = xor(remainder >>> 1, direcn.poly)
            else
                remainder >>>= 1
            end
        end
    end
    remainder
end

function extend(direcn::Backwards{A}, tables::Single{A}, data::Vector{UInt8}, remainder::A) where {A<:U}
    for i in 1:length(data)
        @inbounds word::UInt8 = data[i]
        remainder::A = xor(remainder, convert(A, word))
        remainder = xor(remainder >>> 8, tables.table[1 + remainder & 0xff])
    end
    remainder
end

# generate code for processing a bunch of bytes in a single machine
# word, using multiple tables.  stolen from libz.
for A in (UInt16, UInt32, UInt64, UInt128)
    n_tables = sizeof(A)
    @eval begin
        function extend(direcn::Backwards{$A}, tables::Multiple{$A}, data::MaybeUnalignedVector{$A}, remainder::$A)
            word::$A, tmp::$A, remainder::$A = zero($A), zero($A), remainder
            i = 1
            while true
                @nexprs $UNROLL _ -> begin  # unroll inner loop
                    if i > length(data)
                        break
                    end
                    @inbounds word = data[i]
                    tmp, remainder = xor(remainder, word), zero($A)
                    @nexprs $n_tables t -> begin  # unroll table access
                        remainder = xor(remainder,
                          tables.tables[t][(tmp >>> (($n_tables - t)*8)) & 0xff + 1])
                    end
                    i += 1
                end
            end
            remainder
        end
    end
end

# short-circuit "all UInt8" avoiding a self-recursive loop below
function extend(direcn::Direction{UInt8}, tables::Multiple{UInt8}, data::Vector{UInt8}, remainder::UInt8)
    extend(direcn, Single{UInt8}(tables), data, remainder)
end

function extend(direcn::Direction{A}, tables::Multiple{A}, data::Vector{UInt8}, remainder::A) where {A<:U}
    # this is "clever" - we alias the array of bytes to native machine
    # words and then process those, which typically (64 bits) loads 8
    # bytes at a time.
    tail = length(data) % sizeof(A)
    remainder = extend(direcn, tables, unaligned_reinterpret(A, data), remainder)
    # slurp up the final bytes that didn't fill a complete machine word.
    extend(direcn, Single{A}(tables), data[end-tail+1:end], remainder)
end


text:
module BayesianExperiments

import Statistics: mean, rand, std
using Base
using StaticArrays:SVector
using Distributions:
    pdf,
    Beta, Gamma, LogNormal, Bernoulli, Exponential, Poisson,
    InverseGamma, Normal,
    ContinuousUnivariateDistribution, UnivariateDistribution
using QuadGK

include("util.jl")
include("data.jl")
include("distribution.jl")
include("models/common.jl")
include("models/conjugate.jl")
include("models/bayesfactor.jl")
include("rule.jl")
include("experiment.jl")
include("simulation.jl")

# external dependencies
export
    # Distributions.jl
    Bernoulli,
    Exponential,
    Normal,
    LogNormal

# package exports
export
    # shared methods
    update!,
    rand,
    mean,
    convert,

    # data.jl
    BetaStatistics,
    GammaStatistics,
    NormalStatistics,
    LogNormalStatistics,
    TwoNormalStatistics,
    StudentTStatistics,

    merge,
    tstat,
    tstatpooled,
    tstatwelch,

    # experiment.jl
    Experiment,
    ExperimentABN,
    ExperimentAB,
    ExperimentBF,

    expectedloss,
    expectedlosses,
    probbeatall,
    bayesfactor,

    decide!,
    upliftloss,
    metrics,

    # common.jl
    ProbabilisticModel,

    # conjugate.jl
    ConjugateModel,

    ConjugateBernoulli,
    BernoulliPosteriorSample,

    ConjugateExponential,
    ExponentialPosteriorSample,

    ConjugatePoisson,
    PoissonPosteriorSample,

    ConjugateNormal,
    NormalPosteriorSample,

    ConjugateLogNormal,
    LogNormalPosteriorSample,

    ChainedModel,
    ChainOperator,
    MultiplyOperator,

    samplepost,
    samplestats,
    lognormalparams,

    # bayesfactor.jl
    NormalEffectSize,
    StudentTEffectSize,

    effectsize,

    # rule.jl
    ExpectedLossThresh,
    ProbabilityBeatAllThresh,
    OneSidedBFThresh,
    TwoSidedBFThresh,

    # simulation.jl
    Simulation,
    DataGeneratingDistibutions,

    runonce,
    runsequential,
    updateonce!,

    # util.jl
    catbyrow,
    zstat,
    effsamplesize,
    pooledsd

end # module


text:
#############################
#           pHM.jl          #
#          pH-Mouse         #
#     @author Tanner Lee    #
#      @date 9/27/2018      #
#############################

# Start by loading the initializing data in memory
using DelimitedFiles
print("Which input number would you like to load? \t")
input = readline(stdin)
print("How many iterations would you like to run? \t")
iterations = readline(stdin)
iterations = parse(Int64, iterations)
print("Which room would you like to start in? \t\t")
startingLocation = readline(stdin)
startingLocation = parse(Int64, startingLocation)
print("Enter an intensity for the heatmap {1-10}: \t")
intensity = readline(stdin)
intensity = parse(Int64, intensity)
input = "input" * input * ".phm"
println("!CLEARING TMP FOLDER!")
rm("./tmp", force=true, recursive=true)
mkdir("./tmp")
println("!TMP SUCCESSFULLY CLEARED!")
println("!IMPORTING DATA FROM $input !")
A = readdlm(input, String)
println("!DATA IMPORTED SUCCESSFULLY!")

Arows = size(A)[1]
Acols = size(A)[2]

# Create the transpose manually
# Julia cannot transpose a string matrix
Atrans = Array{String}(undef, Arows, Acols)
for i=1:Arows
    for j=1:Acols
        Atrans[j,i] = A[i,j]
    end
end

println(A)
println(Atrans)

# Create an empty matrix we can concat with
# The columns we want are the pivot columns
probMatrix = zeros(Arows*Acols,Arows*Acols)

# Parse through A to create probability matrix
for i=1:Arows       # 1 < i < rows
    for j=1:Acols   # 1 < j < cols
        global A
        global probMatrix
        up = false
        down = false
        right = false
        left = false
        upDigit = 0
        downDigit = 0
        leftDigit = 0
        rightDigit = 0
        probCounter = 0
        tempString = ""
        tempRoomNumber = 0
        currentRoomString = A[i,j]
        currentRoomNumber = chop(currentRoomString)
        currentRoomNumber = parse(Int64,currentRoomNumber)
        # Check if the current room is a wall, or cheese
        if (string(currentRoomString[length(currentRoomString)]) == "W")
            # Current Room is a wall
            continue
        elseif (string(currentRoomString[length(currentRoomString)]) == "C")
            # Current Room is cheese
            probMatrix[currentRoomNumber,currentRoomNumber] = 1.0
            continue
        end
        # Check how many neighbors A[i,j] has
        # A[i,j] has neighbors at --
        # A[i+1,j], A[i-1,j], A[i,j+1], A[i,j-1]
        # However all of these may not exist
        # Boundary conditions must be considered
        if i != 1
            # We are not on the top row, decrease row
            tempString = A[i-1,j]
            # Get the room number from this neighbor
            tempRoomNumber = chop(tempString)
            # Ensure this room isnt a wall
            if (string(tempString[length(tempString)]) != "W")
                # We can enter this room
                probCounter += 1
                up = true
                upDigit = parse(Int64,tempRoomNumber)
            end
        end
        if i != Arows
            # We are not on the last row, increase row
            tempString = A[i+1,j]
            # Get the room number from this neighbor
            tempRoomNumber = chop(tempString)
            # Ensure this room isnt a wall
            if (string(tempString[length(tempString)]) != "W")
                # We can enter this room
                probCounter += 1
                down = true
                downDigit = parse(Int64,tempRoomNumber)
            end
        end
        if j != 1
            # We are not on the first col, decrease col
            tempString = A[i,j-1]
            # Get the room number from this neighbor
            tempRoomNumber = chop(tempString)
            # Ensure this room isnt a wall
            if (string(tempString[length(tempString)]) != "W")
                # We can enter this room
                probCounter += 1
                left = true
                leftDigit = parse(Int64,tempRoomNumber)
            end
        end
        if j != Acols
            # We are not on the last col, increase col
            tempString = A[i,j+1]
            # Get the room number from this neighbor
            tempRoomNumber = chop(tempString)
            # Ensure this room isnt a wall
            if (string(tempString[length(tempString)]) != "W")
                # We can enter this room
                probCounter += 1
                right = true
                rightDigit = parse(Int64,tempRoomNumber)
            end
        end
        # Create this rooms column in the probability matrix
        if (up)
            # We can go up
            probMatrix[upDigit,currentRoomNumber] = 1/probCounter
        end
        if (down)
            probMatrix[downDigit,currentRoomNumber] = 1/probCounter
        end
        if (left)
            probMatrix[leftDigit,currentRoomNumber] = 1/probCounter
        end
        if (right)
            probMatrix[rightDigit,currentRoomNumber] = 1/probCounter
        end
    end
end

# We now have a probability matrix for our maze
# Create some an input for it
x = zeros(Arows*Acols,1)
x[startingLocation] = intensity

# Use the following code segment to print the probMatrix
# For debugging purposes only
#=
sizeProbMatrix = size(probMatrix)
for i=1:sizeProbMatrix[1]
    print("Row $i:\t [")
    for j=1:15
        print(probMatrix[i,j])
        print(", ")
    end
    println("]")
end
=#

# Complete Luxor Initialization
# Initialize Room Size in pixels
global rSize = 100
# Initialize Canvas Size
global cSize = rSize*Arows
println("!INITIALIZING LUXOR!")
using Luxor
animation = Movie(cSize, cSize, "pHM", 1:iterations)
println("!LUXOR INITIALIZED!")
println("!INITIALIZING COLORSCHEMES!")
using Colors, ColorSchemes
println("!COLORSCHEMES INITIALIZED!")

# Run the simulation
println("!INITIAL INPUT!")
println(x)
println("!STARTING SIMULATION!")

function backdrop(scene, framenumber)
    background("white")
    # Initialize Table and digits
    global table = Table(Arows, Acols, rSize, rSize)
    sethue("black")
    # Loop through each room of the maze and print it to the image
    for i in 1:length(table)
        if (string(Atrans[i][length(Atrans[i])]) == "W")
            # If the current room is a wall make it grey
            box.(table[i], rSize, rSize, :fill)
        end
        # Put a black border on the room and print its name in the middle
        text(Atrans[i], table[i], halign=:center, valign=:middle)
        box.(table[i], rSize, rSize, :stroke)
    end
end

function frame(scene, framenumber)
    global table
    global rSize
    global x
    global probMatrix
    for i=1:length(x)
        sethue(get(ColorSchemes.YlOrRd_7, x[i]))
        setmode("darken")   # Set mode to darken to ensure text isnt covered
        box.(table[i], rSize, rSize, :fill)
    end
    x = probMatrix*x
end

animate(animation, [
    Scene(animation, backdrop, 0:iterations),
    Scene(animation, frame, 0:iterations)],
    creategif=true,
    framerate=20,
    pathname="./images/pHM.gif",
    tempdirectory="./tmp")

println("!SIMULATION COMPLETE!")
println("!FINAL OUTPUT!")
println(x)
println()

max_stars_repo_name: tleecsm/pH-Mouse

max_issues_repo_name: tleecsm/pH-Mouse

max_forks_repo_name: tleecsm/pH-Mouse


text:
using Test
using LinearAlgebra
using SparseArrays
using Random

include("expansion.jl")
include("givens_rotation.jl")
include("ordering.jl")

include("schurfact.jl")
include("sylvester.jl")
include("sort_schur.jl")
include("collect_eigen.jl")

include("partial_schur.jl")
# include("locked_restart.jl")
include("schur_to_eigen.jl")


text:
module IteratedProbabilisticClassification

using
	KernelDensityEstimate,
	Distributions,
	StatsBase,
	Random,
	Gadfly,
	KernelDensityEstimatePlotting,
    Statistics

export
	# Major types and functions used for classification
	ClassificationSystem,
	TuningParameters,
	ClassificationStats,
	defaultClassificationStats,
	classifyOneStep,
	classifyConverge,
	plotClassification1D,
	plotClassificationStats,
	plotEMStatus,
	packDebugResults!,
	EMClassificationRun!,
	dispersePoints,

	# functions under development
	sdc2,
	plotUtil1D,
	plotUtil2D,
	# batchClassifyConverge

	# previous code
	simdata01_true,
	simdata02_true,
	DebugResults,
	defaultDebugResults,
	DataGroundTruth,
	SampleData,
	PALETTE

PALETTE = ["deepskyblue";"red";"magenta";"black";"green";"blue"]

const NothingUnion{T} = Union{Nothing, T}

include("ClassificationUtilities.jl")
include("PlottingUtilities.jl")
include("SimulationData.jl")

function plotMarginalClassificationBelief()
	pl_init = plotKDE([c1_expert,c2_expert],c=PALETTE)
	lyr_smpls = layer(x=samples, y=zeros(size(samples)), Geom.point)
	push!(pl_init.layers, lyr_smpls[1] )
	push!(pl_init.layers, lyr_smpls_trueA[1])
	push!(pl_init.layers, lyr_smpls_trueB[1])

	#Gadfly.draw(PDF("InitialClusterDist.pdf",20cm,12cm),pl_init)
	error("Incomplete function")
end

function plotClassification1D(
			;groundtruth=nothing,
			samples::Array{Float64,2}=Array{Float64,2}(),
			belief::BallTreeDensity=nothing)
		dims = size(samples, 1)

		belief == nothing ? nothing : (dims == Ndim(belief) ? nothing : error("belief and sample dimensions don't match") )

		plr = ceil(Int, sqrt(dims))
		plc = ceil(Int, dims/plr)

		PL = Gadfly.Context[]
		for i in dims
			# pl =
			push!(PL, pl)
		end

		error("Incomplete function")

		nothing
end

end # module


text:
module CloseOpenIntervals

using Static: StaticInt, Zero, One
export CloseOpen, SafeCloseOpen

abstract type AbstractCloseOpen{L <: Union{Int,StaticInt}, U <: Union{Int,StaticInt}} <: AbstractUnitRange{Int} end
for T ∈ (:CloseOpen,:SafeCloseOpen)
  @eval begin
    struct $T{L <: Union{Int,StaticInt}, U <: Union{Int,StaticInt}} <: AbstractCloseOpen{L,U}
      start::L
      upper::U
      @inline $T{L,U}(l::L,u::U) where {L <: Union{Int,StaticInt}, U <: Union{Int,StaticInt}} = new{L,U}(l,u)
    end
    @inline $T(s::S, u::U) where {S,U} = $T{S,U}(s, u)
    @inline $T(len::T) where {T<:Union{Int,StaticInt}} = $T{Zero,T}(Zero(), len)
  end
end

"""
    CloseOpen([start=0], U) <: AbstractUnitRange{Int}
    SafeCloseOpen([start=0], U) <: AbstractUnitRange{Int}

Define close-open unit ranges, i.e. `CloseOpen(0,10)` iterates from from `0:9`.
Close-open ranges can be more convenient in some circumstances, e.g. when
partitioning a larger array.

```julia

function foo(x)
  nt = Threads.nthreads()
  d, r = divrem(length(x), nt)
  i = firstindex(x)
  Threads.@sync for j in 1:nt
    stop = i + d + (r >= j)
    Threads.@spawn bar!(\$(@view(x[CloseOpen(i, stop)])))
    i = stop
  end
end
```
This saves a few `-1`s on the ends of the ranges if using a normal unit range.

`SafeCloseOpen` will not iterate if `U <= start`, while `CloseOpen` doesn't check
for this, and thus the behavior is undefined in that case.
"""
AbstractCloseOpen
"""
    CloseOpen([start=0], U) <: AbstractUnitRange

Iterates over the range start:U-1. Guaranteed to iterate at least once, skipping initial empty check.
See `?AbstractCloseOpen` for more information.
"""
CloseOpen
"""
    SafeCloseOpen([start=0], U) <: AbstractUnitRange

Iterates over the range start:U-1. Will not iterate if it `isempty`, like a typical range, making it
generally the preferred choice over `CloseOpen`.
See `?AbstractCloseOpen` for more information.
"""
SafeCloseOpen

@inline Base.first(r::AbstractCloseOpen) = getfield(r,:start)
@inline Base.first(r::AbstractCloseOpen{StaticInt{F}}) where {F} = F
@inline Base.step(::AbstractCloseOpen) = 1
@inline Base.last(r::AbstractCloseOpen{<:Union{Int,StaticInt},I}) where {I} = getfield(r,:upper) - 1
@inline Base.last(r::AbstractCloseOpen{<:Union{Int,StaticInt},StaticInt{L}}) where {L} = L - 1
@inline Base.length(r::AbstractCloseOpen) = getfield(r,:upper) - getfield(r,:start)
@inline Base.length(r::AbstractCloseOpen{Zero}) = getfield(r,:upper)

@inline Base.iterate(r::CloseOpen) = (i = Int(first(r)); (i, i))
@inline Base.iterate(r::SafeCloseOpen) = (i = Int(first(r)); i ≥ getfield(r, :upper) ? nothing : (i, i))
@inline Base.iterate(r::AbstractCloseOpen, i::Union{Int,StaticInt}) = (i += one(i)) ≥ getfield(r, :upper) ? nothing : (i, i)

import ArrayInterface
ArrayInterface.known_first(::Type{<:AbstractCloseOpen{StaticInt{F}}}) where {F} = F
ArrayInterface.known_step(::Type{<:AbstractCloseOpen}) = 1
ArrayInterface.known_last(::Type{<:AbstractCloseOpen{<:Any,StaticInt{L}}}) where {L} = L - 1
ArrayInterface.known_length(::Type{<:AbstractCloseOpen{StaticInt{F},StaticInt{L}}}) where {F,L} = L - F

Base.IteratorSize(::Type{<:AbstractCloseOpen}) = Base.HasShape{1}()
Base.IteratorEltype(::Type{<:AbstractCloseOpen}) = Base.HasEltype()
@inline Base.size(r::AbstractCloseOpen) = (length(r),)
@inline Base.eltype(r::AbstractCloseOpen) = Int
@inline Base.eachindex(r::AbstractCloseOpen) = StaticInt(1):ArrayInterface.static_length(r)

@static if isdefined(Base.IteratorsMD, :OrdinalRangeInt)
  @inline function Base.IteratorsMD.__inc(state::Tuple{Int,Int,Vararg{Int}}, indices::Tuple{AbstractCloseOpen,Vararg{Base.IteratorsMD.OrdinalRangeInt}})
    rng = indices[1]
    I1 = state[1] + step(rng)
    if Base.IteratorsMD.__is_valid_range(I1, rng) && state[1] != last(rng)
      return true, (I1, Base.tail(state)...)
    end
    valid, I = Base.IteratorsMD.__inc(Base.tail(state), Base.tail(indices))
    return valid, (convert(typeof(I1), first(rng)), I...)
  end
else
  @inline function Base.IteratorsMD.__inc(state::Tuple{Int,Int,Vararg{Int}}, indices::Tuple{AbstractCloseOpen,Vararg})
    rng = indices[1]
    I1 = state[1] + step(rng)
    if Base.IteratorsMD.__is_valid_range(I1, rng) && state[1] != last(rng)
      return true, (I1, Base.tail(state)...)
    end
    valid, I = Base.IteratorsMD.__inc(Base.tail(state), Base.tail(indices))
    return valid, (convert(typeof(I1), first(rng)), I...)
  end
end
end


text:
export logadd, logsub, log1mexp, log1pexp, inf

function inf(::Type{T})::T where T
    one(T) / zero(T)
end

function inf(x::T)::T where T
    inf(T)
end

"log(x+y) from log(x) and log(y)"
function logadd(lx, ly)
    if lx > ly
        return logadd(ly, lx)
    elseif -Inf < lx ≤ ly < Inf
        return ly + log1p(exp(lx - ly))
    else
        return ly
    end
end

logadd(lx, ly, lz, ls...) = foldl(logadd, (logadd(logadd(lx,ly),lz), ls...))

"log(|x-y|) from log(x) and log(y)"
function logsub(lx, ly)
    if min(lx, ly) > -Inf
        return max(lx, ly) + log1mexp(abs(lx - ly))
    else
        return max(lx, ly)
    end
end

"log(1 - exp(-x)), x ≥ 0"
function log1mexp(x)
    # https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf, Eq. 7
    x ≥ 0 || throw(DomainError())
    if x ≤ log(2)
        log(-expm1(-x))
    else
        log1p(-exp(-x))
    end
end

"log(1 + exp(x))"
function log1pexp(x)
    # https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf, Eq. 10
    if x ≤ -37
        exp(x)
    elseif -37 < x ≤ 18
        log1p(exp(x))
    elseif 18 < x ≤ 33.3
        x + exp(-x)
    else
        x
    end
end


text:
using NeuralNets
net = buildnet()
using Graphs
to_dot(net.g,"net.dot")
run(`dot -Tpng net.dot -o net_diagram.png`)
println("Network graph written to net_diagram.png")


text:
board = [[" ", " ", " "],
         [" ", " ", " "],
         [" ", " ", " "]]

function print_board(board)
    println(" ", join(board[1], " | "))
    println(" ——┼———┼——")
    println(" ", join(board[2], " | "))
    println(" ——┼———┼——")
    println(" ", join(board[3], " | "))
end

print_board(board)

max_stars_repo_name: soumitradev/youtube-tutorial

max_issues_repo_name: soumitradev/youtube-tutorial

max_forks_repo_name: soumitradev/youtube-tutorial


text:
function example_plotly(route="",client_id="")
    n = 10
    data = [Scatter(;x=1:n,y=rand(n))]
    layout = Layout(;
        title = Title(text="Create a new <div> and plot a single trace"),
        width = 650,
        height = 450
    )
    config = Config(displayModeBar=false)
    if isempty(route) && isempty(client_id)
        Plotly.newPlot("Plot1",data,layout=layout,config=config)
    else
        Plotly.newPlot(route,client_id,"Plot1",data,layout=layout,config=config)
    end

    trace1 = Scatter(x = 1:n, y = rand(n))
    trace2 = Scatter(x = 1:n, y = rand(n))
    data = [trace1, trace2]
    layout = Layout(
        title = Title(text="Create a new <div> and plot an array of traces"),
        width = 650,
        height = 450
    )
    config = Config(displayModeBar=false)
    if isempty(route) && isempty(client_id)
        Plotly.newPlot("Plot2",data,layout=layout,config=config)
    else
        Plotly.newPlot(route,client_id,"Plot2",data,layout=layout,config=config)
    end

    trace1 = Scatter(
        x = [1, 2, 3, 4],
        y = [10, 15, 13, 17],
        mode = "markers")

    trace2 = Scatter(
        x = [2, 3, 4, 5],
        y = [16, 5, 11, 10],
        mode = "lines")

    trace3 = Scatter(
        x = [1, 2, 3, 4],
        y = [12, 9, 15, 12],
        mode = "lines+markers")

    data = [trace1, trace2, trace3]

    layout = Layout(
        title = Title(text="Line and Scatter Plot"),
        width = 650,
        height = 450)

    config = Config(displayModeBar=false)
    if isempty(route) && isempty(client_id)
        Plotly.newPlot("Plot3",data,layout=layout,config=config)
    else
        Plotly.newPlot(route,client_id,"Plot3",data,layout=layout,config=config)
    end

    country = ["Switzerland (2011)", "Chile (2013)", "Japan (2014)",
            "United States (2012)", "Slovenia (2014)", "Canada (2011)",
            "Poland (2010)", "Estonia (2015)", "Luxembourg (2013)",
            "Portugal (2011)"]

    votingPop = [40, 45.7, 52, 53.6, 54.1, 54.2, 54.5, 54.7, 55.1, 56.6]
    regVoters = [49.1, 42, 52.7, 84.3, 51.7, 61.1, 55.3, 64.2, 91.1, 58.9]

    trace1 = Scatter(
        x=votingPop,
        y=country,
        mode="markers",
        name="Percent of estimated voting age population",
        marker=Marker(
            color="rgba(156, 165, 196, 0.95)",
            line = Line(
                color="rgba(156, 165, 196, 1.0)",
                width=1),
            size=16,
            symbol="circle"))

    trace2 = Scatter(
        x=regVoters,
        y=country,
        mode="markers",
        name="Percent of estimated registered voters")

    trace2.attributes.marker = Marker(
        color = "rgba(204, 204, 204, 0.95)",
        line = Line(
            color = "rgba(217, 217, 217, 1.0)",
            width = 1),
        symbol = "circle",
        size = 16)

    data = [trace1, trace2]

    layout = Layout(
        paper_bgcolor = "rgb(254, 247, 234)",
        plot_bgcolor = "rgb(254, 247, 234)",
        title = Title(text="Votes cast for ten lowest voting age population in OECD countries"),
        width = 600,
        height = 600,
        hovermode = "closest",
        margin = Margin(l = 140, r = 40, b = 50, t = 80),
        xaxis = Axis(
            showgrid = false,
            showline = true,
            linecolor = "rgb(102, 102, 102)",
            title = Title(
                font = Font(
                    color = "rgb(204, 204, 204)"))))
            # tick font_font_color = "rgb(102, 102, 102)",
            # autotick = false,
            # dtick = 10,
            # ticks="outside",
            # tickcolor="rgb(102, 102, 102)",
            # legend=Legend(
            #     font_size=10,
            #     yanchor="middle",
            #     xanchor="right"),
            #             )

    config = Config(displayModeBar=false,scrollZoom=false,doubleClick="reset")
    if isempty(route) && isempty(client_id)
        Plotly.newPlot("Plot4",data,layout=layout,config=config)
    else
        Plotly.newPlot(route,client_id,"Plot4",data,layout=layout,config=config)
    end
    return
end

# example_plotly()


text:
import ReferenceFrameRotations
path = @__DIR__
include(joinpath(path,"Turtle_struct.jl"))
include(joinpath(path,"LSystem.jl"))
using Plots
using CSV, DataFrames
gr()

RF = ReferenceFrameRotations

#========================= Export =======================================#
export penup,pendown,goto,cw,ccw,forward,turtle_plot
export draw_line,draw_squares,draw_triangles,draw_lsystem
export main_proj

#========================= Function =======================================#

"""
    penup(turtle::Turtles)

this function create a new searies data with in the main plot (turtle.plt)
and set turtle.pen_up = true
"""
function penup(turtle::Turtles)
    turtle.pen.up = true
    turtle.history = [turtle.history;[NaN,NaN,NaN]']
    # plot!(turtle.plt,[turtle.pos[1]],[turtle.pos[2]]);
end # function penup

"""
    pendown(turtle::Turtles)

this function set turtle.pen_up = false
and store position data
"""
function pendown(turtle::Turtles)
    turtle.pen.up = false
    store_position_data(turtle)
end # function pendown

"""
     goto(turtle::Turtles,newpos::Array{T} where T<:Real)

this function update the turtle position.
the input is a new point in [x,y] coords that the turtle jump to.
"""
function goto(turtle::Turtles,newpos::Array{T}  where T<:Real)
    if length(newpos) != 2
        DimensionMismatch("The new point to goto must have 2 elements [x,y]")
    end

    if !(newpos isa Vector)
        newpos = vec(newpos)
    end

    turtle.pos = newpos
end #function goto

"""
    cw(turtle::Turtles,ang::Real)

this function update the turtle heading clock wise direction.
the input ang must be in [Rad], negative or positive angle
are treated as positive.
"""
function cw(turtle::Turtles,ang::Real)
    turtle.heading -= abs(ang)
end # function cw

"""
    ccw(turtle::Turtles,ang::Real)

this function update the turtle heading counter clock wise direction.
the input ang must be in [Rad], negative or positive angle
are treated as positive.
"""
function ccw(turtle::Turtles,ang::Real)
    turtle.heading += abs(ang)
end # function ccw

"""
    forward(turtle::Turtles,step::AbstractFloat)

this function will advance the turtle forward in the heading direction
for distance of 'step'
"""
function forward(turtle::Turtles,step::Real)
    R = RF.angle_to_dcm(turtle.heading,0,0)
    turtle.pos += vec(step * R[1,1:2])
    if !turtle.pen.up
        store_position_data(turtle)
    end
end # function forward

"""
    turtle_plot(data)

this function plot datain with in the turtle.main_scene
"""
function turtle_plot(turtle::Turtles)
    if !turtle.pen.up
        turtle.plt.series_list[end].plotattributes[:x] = turtle.history[:,1]
        turtle.plt.series_list[end].plotattributes[:y] = turtle.history[:,2]
        turtle.plt.series_list[2].plotattributes[:x] = turtle.pos[1]
        turtle.plt.series_list[2].plotattributes[:y] = turtle.pos[2]
    end
end # function turtle_plot

"""
    escaped(turtle::Turtles)

this function check if the turtle escape the bag
"""
function escaped(turtle::Turtles)
    res = false
    if  (turtle.pos[1] < minimum(turtle.bag.vertexes[:,1]) ||
        turtle.pos[1] > maximum(turtle.bag.vertexes[:,2])) ||
        (turtle.pos[2] < minimum(turtle.bag.vertexes[:,2]) ||
        turtle.pos[2] > maximum(turtle.bag.vertexes[:,2]))
        res = true
    end
    return res
end # function escaped

"""
    store_position_data(turtle::Turtles)

this function save to last step the turtle did within an array
"""
function store_position_data(turtle::Turtles)
    if !turtle.pen.up
        d = [turtle.pos;escaped(turtle)]'
        turtle.history = [turtle.history;d]
    end
end # function store_position_data

"""
    push!(turtle::Turtles)

this function push the current turtle position into vector 'turtle.junctions'
"""
function Base.push!(turtle::Turtles)
    push!(turtle.junctions,turtle.pos)
end # function push!

"""
    pop!(turtle::Turtles)

this function pop the last position data from 'turtle.junctions' vector
"""
function Base.pop!(turtle::Turtles)
    return pop!(turtle.junctions)
end # function pop!

"""
    draw_line(turtle::Turtles,turn_rad::Real, step::Real)

this function draw a line, by first rotating the turtle at 'turn_rad'
(ccw if 'turn_rad' > 0 and cw if 'turn_rad' < 0)
then move the turtle at the new direction distance of 'step'
the function check if the turtle escaped the bag and return true or false
"""
function draw_line( turtle::Turtles,
                    turn_rad::Real,
                    step::Real)
    if turn_rad < 0
        cw(turtle,turn_rad)
    elseif turn_rad > 0
        ccw(turtle,turn_rad)
    end
    forward(turtle,step)
    store_position_data(turtle)
    return escaped(turtle)
end # function draw_line

"""
    draw_square(turtle::Turtles, turn_rad::Real, step::Real, do_till_escaped::Bool)

this function draw a square rotated at 'turn_rad' and which rib size = 'step'
if 'do_till_escaped' == true, the function will stop when turtle escaped
"""
function draw_square(   turtle::Turtles,
                        turn_rad::Real,
                        step::Real,
                        do_till_escaped::Bool)
    res = false
    escaped = draw_line(turtle,turn_rad,step)
    if do_till_escaped && escaped
        res = true
    else
        for ii in 1:3
            escaped = draw_line(turtle,π/2,step)
            if do_till_escaped && escaped
                res = true
                break
            end
        end
    end
    return res
end # function draw_square

"""
    draw_squares(turtle::Turtles, N::Int, turn_rad::Real, step::Real, increase_step::Real, do_till_escaped::Bool)

this function draw 'N' squares,
the first square start at the curent turtle pos,
then it is rotated at initial 'turn_rad' radians,
with initial rib size = step.
the next square rib = 'step' + i*'increase_step'
if 'do_till_escaped' == true, the function will stop when turtle escaped
"""
function draw_squares(  turtle::Turtles,
                        N::Int,
                        turn_rad::Real,
                        step::Real,
                        increase_step::Real,
                        do_till_escaped::Bool)
    res = false
    escaped = draw_square(turtle, turn_rad, step, do_till_escaped)
    if do_till_escaped && escaped
        res = true
    else
        for ii in 1:N-1
            penup(turtle)
            forward(turtle,increase_step)
            cw(turtle,π/2)
            forward(turtle,increase_step)
            cw(turtle,π)
            pendown(turtle)
            escaped = draw_square(turtle, 0, step + ii*2*increase_step, do_till_escaped)
            if do_till_escaped && escaped
                res = true
                break
            end
        end
    end
    return res
end # function draw_squares

"""
    draw_triangles(turtle::Turtles, turn_rad::Real, N::Int, ang::Real, step::Real, increase_step::Real, do_till_escaped::Bool)

this function draw spriangles with angle = 'ang'
and initial rib size = 'step'
and rib = step + i * increase_step
if 'do_till_escaped' == true, the function will stop when turtle escaped
"""
function draw_triangles(turtle::Turtles,
                        turn_rad::Real,
                        N::Int,
                        ang::Real,
                        step::Real,
                        increase_step::Real,
                        do_till_escaped::Bool)
    res = false
    escaped = draw_line(turtle,turn_rad,step)
    if do_till_escaped && escaped
        res = true
    else
        for ii in 1:3*N-1
            escaped = draw_line(turtle,ang,step + ii*increase_step)
            if do_till_escaped && escaped
                res = true
                break
            end
        end
    end
    return res
end # function draw_triangles

"""
    draw_lsystem(turtle::Turtles, lsys_data::Dict, do_till_escaped::Bool, lsys_data::Tuple, step::Real, turn_rad::Real, iterations::Int)

this function will draw a Lsystem expration
"""
function draw_lsystem(  turtle::Turtles,
                        do_till_escaped::Bool,
                        lsys_data::Tuple,
                        step::Real,
                        turn_rad::Real,
                        iterations::Int)
    lsystem = LSystem(lsys_data[1],lsys_data[2])
    lsystem.state = lsystem.initial_state
    evaluate(lsystem,iterations,debug = false)
    render(lsystem, turtle, step, turn_rad, do_till_escaped, debug = false)
    return true
end # function draw_lsystem

"""
    save_to_file(turtle::Turtles)

this function save the turtle.history data to csv file
after filterring NaN values
"""
function save_to_file(turtle::Turtles,csvFilenameIn::String = "run_1.csv")
    if occursin(r".csv",csvFilenameIn)
        splitName = split(fullFileName,'.')
        endSplitName = splitName[1][end]
        if isnumeric(endSplitName)
            fullFileName = csvFilenameIn
        else
            fullFileName = splitName[1] * "_1.csv"
        end
    else
        if isnumeric(csvFilenameIn[end])
            fullFileName = csvFilenameIn * ".csv"
        else
            fullFileName = csvFilenameIn * "_1.csv"
        end
    end

    while isfile(fullFileName)
        splitName = split(fullFileName,'.')
        endSplitName = splitName[1][end]
        if isnumeric(endSplitName)
            numlength = sum(map(isnumeric,collect(splitName[1])))
            fileNameApendix = string(parse(Int64,splitName[1][end - numlength + 1 : end]) + 1)
        else
            fileNameApendix = "1"
        end

        fullFileName = splitName[1][1:end - length(fileNameApendix)] * fileNameApendix * ".csv"
    end

    datatofile = DataFrame(turtle.history)
    filter!(row->!isnan(row.x1),datatofile)
    CSV.write(fullFileName,  datatofile, writeheader=false)
end # function save_to_file

"""
    popup_my_input(prompt::String = "", df_val::DefaultVals)

this function ask the user to enter needed params for the called function
"""
function popup_my_input(df_val, prompt::String = "")
    if df_val isa Tuple
        println(prompt)
        d = Dict{String,String}()
        addData = true
        while addData
            println("\nEnter Var")
            varStr = chomp(readline())
            if isempty(varStr)
                println("default value was taken")
                d = df_val[1]
                addData = false
                break
            end
            println("Enter Expration")
            expStr = chomp(readline())
            println("recived $varStr => $expStr as Lsystm Expration")
            if haskey(d,varStr)
                println("Var exsist enter a new one")
            else
                d[varStr] = expStr
            end
        end
        println("\nEnter the starting point")
        startStr = chomp(readline())
        if isempty(startStr)
            startStr = df_val[2]
        else
            println("default value was taken")
            println("recived starting point = $startStr")
        end
        val = (d,startStr)
    else
        println(prompt)
        str = chomp(readline())
        if isempty(str)
            val = df_val
        else
            val = parse(typeof(df_val),str)
        end
    end
    return val
end # function popup_my_input

"""
    get_inputs_list(f::Function)

this function extract the given function, it's input argumant list
"""
function get_inputs_list(f::Function)
    list = []
    for m in methods(f)
        push!(list,string.(Base.method_argnames(m))...)
    end
    inputlist = Dict{String,Symbol}()
    for str in list
        splitlist = split(str,"\0")
        for s in splitlist
            if s == "#self#" || s == ""
                continue
            end
            if !haskey(inputlist,s)
                inputlist[s] = Symbol(s)
            end
        end
    end
    return inputlist
end # function get_inputs_list

"""
    get_my_input(dv::DefaultVals, ex)

this function call the function 'popup_my_input', to ask for the user to enter
the needed value for the needed var
"""
function get_my_input(dv::DefaultVals, ex)
    println("give value to $(string(ex)) or press ENTER for deafult")
    if ex == :turn_rad
        dv.turn_rad = popup_my_input(dv.turn_rad, "enter angle, in Deg, to rotate the turtle (positive -> ccw, negative -> cw), deafult = $(rad2deg(dv.turn_rad))")
        dv.turn_rad = deg2rad(dv.turn_rad)
        println("value entered = $(rad2deg(dv.turn_rad))")
    end
    if ex == :step
        dv.step = popup_my_input(dv.step, "enter step size, deafult = $(dv.step)")
        println("value entered = $(dv.step)")
    end
    if ex == :N
        dv.N = popup_my_input(dv.N, "enter number of reapets, deafult = $(dv.N)")
        println("value entered = $(dv.N)")
    end
    if ex == :ang
        dv.ang = popup_my_input(dv.ang, "enter triangle angle, in Rad, deafult = $(dv.ang)")
        println("value entered = $(dv.ang)")
    end
    if ex == :increase_step
        dv.increase_step = popup_my_input(dv.increase_step, "enter step increase size, deafult = $(dv.increase_step)")
        println("value entered = $(dv.increase_step)")
    end
    if ex == :do_till_escaped
        dv.do_till_escaped = popup_my_input(dv.do_till_escaped, "enter do_till_escaped, deafult = $(dv.do_till_escaped)")
        println("value entered = $(dv.do_till_escaped)")
    end
    if ex == :lsys_data
        dv.lsys_data = popup_my_input(dv.lsys_data, "enter lsys_data, deafult = $(dv.lsys_data)")
        println("value entered = $(dv.lsys_data)")
    end
    if ex == :iterations
        dv.iterations = popup_my_input(dv.iterations, "enter iterations, deafult = $(dv.iterations)")
        println("value entered = $(dv.iterations)")
    end
    println("===============================================================")
    println("")
end # function my_input

#================================= MACROS ===================================#
# see : https://discourse.julialang.org/t/undefvarerror-x-not-defined-when-calling-a-macro-outside-of-its-module/20201/3
# for explanation on the use of 'esc' in macro
macro run_type(type,func_type_dict,turtle,dv)
    quote
        local _type = $(esc(type))
        local _func_type_dict = $(esc(func_type_dict))
        local _turtle = $(esc(turtle))
        local _dv = $(esc(dv))
        inputs_needed = get_inputs_list(_func_type_dict[_type][1])
        for (k,v) in inputs_needed
            if v == :turtle
                continue
            end
            get_my_input(_dv,v)
        end
        if _dv.do_till_escaped
            _func_type_dict[_type][3](_turtle,_dv)
        else
            _func_type_dict[_type][2](_turtle,_dv)
        end
    end
end # macro run_type

# ------------------------  Tests -------------------------------------------
function main_line(turtle::Turtles, dv::DefaultVals)
    # turtle = Turtles()
    draw_line(turtle,dv.turn_rad,dv.step)
    turtle_plot(turtle)
    save_to_file(turtle,"main_line")
    return turtle.plt
end

function main_line_till_escaped(turtle::Turtles, dv::DefaultVals)
    # turtle = Turtles()
    escaped = false
    while !escaped
        escaped = draw_line(turtle,dv.turn_rad,dv.step)
        turtle_plot(turtle)
    end
    save_to_file(turtle,"main_line_till_escaped")
    return turtle.plt
end

function main_squares(turtle::Turtles, dv::DefaultVals)
    # turtle = Turtles()
    draw_squares(turtle,dv.N,dv.turn_rad,dv.step,dv.increase_step,dv.do_till_escaped)
    turtle_plot(turtle)
    save_to_file(turtle,"main_squares")
    return turtle.plt
end

function main_squares_till_escaped(turtle::Turtles, dv::DefaultVals)
    # turtle = Turtles()
    escaped = false
    while !escaped
        escaped = draw_squares(turtle,dv.N,dv.turn_rad,dv.step,dv.increase_step,dv.do_till_escaped)
        turtle_plot(turtle)
    end
    save_to_file(turtle,"main_squares_till_escaped")
    return turtle.plt
end

function main_triangles(turtle::Turtles, dv::DefaultVals)
    # turtle = Turtles()
    draw_triangles(turtle, dv.turn_rad, dv.N, dv.ang, dv.step, dv.increase_step, dv.do_till_escaped)
    turtle_plot(turtle)
    save_to_file(turtle,"main_triangles")
    return turtle.plt
end

function main_triangles_till_escaped(turtle::Turtles, dv::DefaultVals)
    # turtle = Turtles()
    escaped = false
    while !escaped
        escaped = draw_triangles(turtle, dv.turn_rad, dv.N, dv.ang, dv.step, dv.increase_step, dv.do_till_escaped)
        turtle_plot(turtle)
    end
    save_to_file(turtle,"main_triangles_till_escaped")
    return turtle.plt
end

function main_Lsystem(turtle::Turtles,dv::DefaultVals)
    draw_lsystem(turtle, dv.do_till_escaped, dv.lsys_data, dv.step, dv.turn_rad, dv.iterations)
    turtle_plot(turtle)
    save_to_file(turtle,"main_Lsystem")
    return turtle.plt
end

"""
    main_porj(moveType::String, stop_condition::String; N::Int, trun_rad::Real)

documentation
"""
function main_proj(move_type_in::String)
    turtle = Turtles()
    dv = DefaultVals()
    func_type_dict = Dict(  :line => [draw_line, main_line, main_line_till_escaped],
                            :squares => [draw_squares, main_squares, main_squares_till_escaped],
                            :triangles => [draw_triangles, main_triangles, main_triangles_till_escaped],
                            :lsystem => [draw_lsystem,main_Lsystem,main_Lsystem])
    move_type = Symbol(move_type_in)
    @run_type(move_type,func_type_dict,turtle,dv)
end # function main_proj


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

"""
    Coboundary{P,Q,T}

The co-boundary relation from rank `P` to greater rank `Q` for
topology of type `T`.
"""
struct Coboundary{P,Q,T<:Topology} <: TopologicalRelation
  topology::T
end

Coboundary{P,Q}(topology::T) where {P,Q,T} = Coboundary{P,Q,T}(topology)

# -------------------
# HALF-EDGE TOPOLOGY
# -------------------

function (𝒞::Coboundary{0,1,T})(vert::Integer) where {T<:HalfEdgeTopology}
  t = 𝒞.topology
  𝒜 = Adjacency{0}(t)
  [edge4pair((vert, other), t) for other in 𝒜(vert)]
end

function (𝒞::Coboundary{0,2,T})(vert::Integer) where {T<:HalfEdgeTopology}
  e = half4vert(vert, 𝒞.topology)

  # initialize result
  elements = [e.elem]

  # search in CCW orientation
  p = e.prev
  h = p.half
  while !isnothing(h.elem) && h != e
    push!(elements, h.elem)
    p = h.prev
    h = p.half
  end

  # if border edge is hit
  if isnothing(h.elem)
    # search in CW orientation
    h = e.half
    while !isnothing(h.elem)
      pushfirst!(elements, h.elem)
      n = h.next
      h = n.half
    end
  end

  elements
end

function (𝒞::Coboundary{1,2,T})(edge::Integer) where {T<:HalfEdgeTopology}
  e = half4edge(edge, 𝒞.topology)
  isnothing(e.half.elem) ? [e.elem] : [e.elem, e.half.elem]
end


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

"""
Bootstrap function
    bootstrap(data::Array{Float64}, ese::Array{Float64},nbsample::Int64)

Bootstrap of Raman spectra. We generate new datapoints with the basis of existing data and their standard deviation
"""
function bootstrap(data::Array{Float64}, ese::Array{Float64},nbsample::Int64):
    bootsamples = zeros(size(data)[1],nbsample)
    for i = 1:nbsample
        randn!(bootsamples[:,i]) .* ese[:] + data[:]
    end
    return bootsamples
end

function poly(p::Vector{Float64},x::Array{Float64})
    segments = zeros(size(x)[1],size(p)[1])
    for i = 1:size(p)[1]
        segments[:,i] = p[i].*x[:].^(i-1)
    end
    return sum(segments,2)
end

#For Gaussian peaks in spectra:
function gaussiennes(g_amplitudes::Array{Float64},g_frequency::Array{Float64},g_hwhm::Array{Float64},x::Array{Float64})
    segments = zeros(size(x)[1],size(g_amplitudes)[1])
    for i = 1:size(g_amplitudes)[1]
        segments[:,i] = g_amplitudes[i] .*exp(-log(2) .* ((x[:,1]-g_frequency[i])./g_hwhm[i]).^2)
    end
    return sum(segments,2), segments
end

#The real normal distribution / gaussian function
function normal_dist(nd_amplitudes::Array{Float64},nd_centres::Array{Float64},nd_sigmas::Array{Float64},x::Array{Float64})
    segments = zeros(size(x)[1],size(nd_amplitudes)[1])
    for i = 1:size(nd_amplitudes)[1]
        segments[:,i] = nd_amplitudes[i]./(nd_sigmas[i].*sqrt(2*pi))  .*  exp(- (x[:,1]-nd_centres[i]).^2 ./ (2.*nd_sigmas[i].^2))
    end
    return sum(segments,2), segments
end


text:
module SA
using PyPlot
using ProgressMeter
import JSON
include("random_rewire.jl")
include("utilities.jl")
include("algConnectivity.jl")
include("arg_config.jl")
include("plotting.jl")

export optimize

function Paccept(c,c_prime,T,its,maxits)
    if c_prime >= c
        return 1
    else
        normalizedP = exp(-its/(maxits/8))
        return normalizedP/(c-c_prime)
    end
end

function optimize(ARGS="")
    args = parse_commandline(ARGS)
    n = args["vertices"]
    e = args["edges"]
    if args["out_folder"] == ""
        dt = replace(string(now()),"T","_")
        dt = replace(dt,":","_")
        dt = replace(dt,"-","_")
        out_fldr = string("res_",dt)
    else
        out_fldr = args["out_folder"]
    end

    if args["simple_graph"]
        g = simple_NE_sgraph(n,e)
    else
        g = simple_NE_graph(n,e)
    end

    if !args["no_graph_dot"]
        write_to_file(g,out_fldr,args["init_fname"])
    end
    if args["plot_graphs"]
        plot(g)
    end
    if args["save_graphs_plot"]
        plot_to_file(g,out_fldr,args["init_fname"])
    end
    save_json(args,out_fldr,"args.txt")
    current_state = deepcopy(g)
    tot_its = args["iterations"]
    edge_change_prob = args["edge_rewire"]
    T = args["temperature"]
    c = AlgConnectivity(current_state,args["sparse_rep"])
    if !args["no_console_output"]
        @printf("Original connectivity: %.3f\n",c)
    end
    if !args["no_figs"]
        plot_deg_hist(g,out_fldr,"d_hist_initial",c)
    end
    c_hist = zeros(tot_its,1)

    p = Progress(tot_its, 30, "Progress: ", 50)
    tic()
    for i = 1:tot_its
        update!(p,i)
        c_hist[i] = c
        candidate = rewire(deepcopy(current_state),edge_change_prob)
        c_prime = AlgConnectivity(candidate,args["sparse_rep"])
        accept_p = Paccept(c,c_prime,T,i,tot_its)

        if rand() < accept_p
            # Moving here
            c = c_prime
            current_state = deepcopy(candidate)
        else
            continue
        end
    end
    t = toq()
    if !args["no_console_output"]
        @printf("Final connectivity: %.3f\n",c)
    end
    if !args["no_graph_dot"]
        write_to_file(current_state,out_fldr,args["final_fname"])
    end
    if args["plot_graphs"]
        plot(current_state)
    end
    if args["save_graphs_plot"]
        plot_to_file(current_state,out_fldr,args["final_fname"])
    end

    if !args["no_figs"]
        plot_c_hist(c_hist,out_fldr,"c_hist")

        plot_deg_hist(current_state,out_fldr,"d_hist_final",c)
    end
    return t
end
end


text:
"""
    reduced_eigenvalues(basis::SingleBlockBasis, basis_A::MultiBlockBasis, basis_B::MultiBlockBasis, A_start, wf::AbstractVector{Float64})

Compute the reduced density matrix eigenvalues of state `wf` reduced to
subsystem A.
"""
function reduced_eigenvalues(basis::SingleBlockBasis{N}, basis_A::MultiBlockBasis{NA}, basis_B::MultiBlockBasis{NB}, A_start, wf::AbstractVector{Float64}) where {N,NA,NB}
    N == NA + NB || throw(DomainError((NA, NB), "Incorrect subspace sizes."))
    A_start >= 1 || throw(DomainError((A_start, NA), "Subsystem must be inside system."))
    A_start + NA - 1 <= N || throw(DomainError((A_start, NA), "Subsystem must be inside system."))

    v = Array{Int}(undef, 2N)
    eigvals = Float64[]

    for (label_A, block_A) in basis_A.blocks
        vectors_A = block_A.vectors

        # Only use blocks that have complementary lp and m.
        label_B = BlockLabel(mod(basis.label.lp-label_A.lp, 2), basis.label.m-label_A.m)
        haskey(basis_B.blocks, label_B) || continue

        block_B = basis_B.blocks[label_B]
        vectors_B = block_B.vectors

        mat = zeros(block_B.size, block_A.size)

        for col in 1:block_A.size
            l_total_col = 0
            for i in 1:2NA
                v[2A_start+i-2] = vectors_A[col][i]
                i % 2 == 1 && (l_total_col += v[2A_start+i-2])
            end

            for row in 1:block_B.size
                l_total = l_total_col
                for i in 1:(2A_start-2)
                    v[i] = vectors_B[row][i]
                    i % 2 == 1 && (l_total += v[i])
                end
                for i in (2A_start+2NA-1):2N
                    v[i] = vectors_B[row][i-2NA]
                    i % 2 == 1 && (l_total += v[i])
                end

                l_total <= basis.l_total_max || continue

                mat[row, col] = wf[basis.block.lookup[v]]
            end
        end

        append!(eigvals, abs2.(svdvals(mat)))
    end

    abs(sum(eigvals) - 1.0) < 1e-12 || @warn("Bad trace: $(abs(sum(eigvals) - 1.0))")

    eigvals
end

"""
    S_vn(eigvals::AbstractVector{Float64})

Von Neumann entropy of `eigvals`.
"""
S_vn(eigvals::AbstractVector{Float64}) = -sum(x * log(x) for x in eigvals if x > 0)

"""
    S_renyi(eigvals::AbstractVector{Float64}, alpha=2)

Order-`alpha` Rényi entropy of `eigvals`.
"""
S_renyi(eigvals::AbstractVector{Float64}, alpha=2) = log(sum(eigvals.^alpha)) / (1 - alpha)


text:
# Compute the Hamming Distance Between Two Strings

calc_hamming_distance(s1, s2) = count(((c1, c2),) -> c1 != c2, zip(s1, s2))

function main()
    dna1 = readline()
    dna2 = readline()
    result = calc_hamming_distance(dna1, dna2)
    println(result)
end

main()