text:
#=
Character classification functions

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones,
Licensed under MIT License, see LICENSE.md
=#

# Todo: add correct definitions for is_letter, is_alphabetic

# Recommended by deprecate
@static if V6_COMPAT
    text_width(ch::Char) = charwidth(ch)

    import Base: is_assigned_char, isnumber
    Base.is_assigned_char(ch::Chr) = is_assigned(ch)
    Base.charwidth(ch::Chr) = text_width(ch)
    isnumber(val::Chr) = is_numeric(val)
else
    Base.Unicode.isassigned(ch::Chr) = is_assigned(ch)
    is_graphic(ch::Char) = is_graphic(codepoint(ch))
    is_alphanumeric(ch::Char) = is_alphanumeric(codepoint(ch))
end

############################################################################

## character column width function ##

text_width(ch::UInt8)    = Int(ifelse(ch < 0x7f, ch > 0x1f, ch > 0x9f & ch != 0xad))
text_width(ch::UInt16)   = utf8proc_charwidth(ch)
text_width(ch::UInt32)   = utf8proc_charwidth(ch)
text_width(ch::Chr)      = text_width(codepoint(ch))
text_width(ch::ASCIIChr) = Int(32 <= codepoint(ch) <= 126)

############################################################################

# returns code in 0:30 giving Unicode category
@inline category_code(ch::CodeUnitTypes) = ch <= 0x10ffff ? utf8proc_cat(ch) : Cint(30)
@inline category_code(ch::Chr) = category_code(codepoint(ch))

# more human-readable representations of the category code
@inline category_abbrev(ch::CodeUnitTypes) = ch <= 0x10ffff ? utf8proc_cat_abbr(ch) : "In"
@inline category_abbrev(ch::Chr)           = category_abbrev(codepoint(ch))

category_string(ch::CodeUnitTypes) = Uni.category_strings[category_code(ch) + 1]
category_string(ch::Chr)           = category_string(codepoint(ch))

is_assigned(ch::CodeUnitTypes) = category_code(ch) != Uni.Cn
is_assigned(ch::Chr)           = is_assigned(codepoint(ch))

_cat_mask(a) = UInt(a)
@inline _cat_mask(a, b) = (UInt(1) << (a%UInt)) | (UInt(1) << (b%UInt))
@inline _cat_mask(rng::(@static V6_COMPAT ? Range : AbstractRange)) =
    ((UInt(2) << (rng.stop%UInt)) - UInt(1)) & ~((UInt(1) << (rng.start%UInt)) - UInt(1))

@inline _check_mask(ch, mask) = ((UInt(1) << (category_code(ch)%UInt)) & mask) != 0

## libc character class predicates ##

@inline _iscntrl(ch) = (ch <= 0x1f) | (0x7f <= ch <= 0x9f)
@inline _isdigit(ch) = (ch - ('0'%UInt8)) <= 9
@inline _isxdigit(ch) = _isdigit(ch) | (ch - ('A'%UInt8) < 6) | (ch - ('a'%UInt8) < 6)

const _isupper_mask   = _cat_mask(Uni.Lu, Uni.Lt)
const _isalpha_mask   = _cat_mask(Uni.Lu : Uni.Lo)
const _isnumeric_mask = _cat_mask(Uni.Nd : Uni.No)
const _ispunct_mask   = _cat_mask(Uni.Pc : Uni.Po)
const _isprint_mask   = _cat_mask(Uni.Lu : Uni.Zs)
const _isgraph_mask   = _cat_mask(Uni.Lu : Uni.So)
const _isalnum_mask   = _isnumeric_mask | _isalpha_mask

############################################################################
# Definitions for characters in the ASCII subset of Unicode

@inline _is_lower_a(ch) = (ch%UInt8 - 'a'%UInt8) < 0x1a
@inline _is_upper_a(ch) = (ch%UInt8 - 'A'%UInt8) < 0x1a

const _isnumeric_a = _isdigit
@inline _ispunct_a(ch) = ((UInt128(1) << ch) & 0x2800_0000_b800_0001_8c00_f7ee_0000_0000) != 0
@inline _isspace_a(ch) = (ch == 0x20) | (0x9 <= ch <= 0xd)
@inline _isalpha_a(ch) = _is_lower_a(ch) | _is_upper_a(ch)
@inline _isalnum_a(ch) = _isdigit(ch) | _isalpha_a(ch)
@inline _isprint_a(ch) = 0x20 <= ch < 0x7f
@inline _isgraph_a(ch) = 0x20 < ch < 0x7f

############################################################################
# Definitions for characters in the Latin1 subset of Unicode, but not in the ASCII subset

@inline _is_lower_l(c)   = ((0xdf <= c <= 0xff) & (c != 0xf7)) | (c == 0xb5)
@inline _is_upper_l(c)   = (0xc0 <= c%UInt8 <= 0xde) & (c != 0xd7)

@inline _isnumeric_l(ch) = (ch <= 0xbe && ((1<<(ch-0xb2)) & 0x1c83) != 0)
@inline _ispunct_l(ch)   = ((UInt64(1) << (ch-0x80)) & 0x88c0_0882_0000_0000) != 0
@inline _isspace_l(ch)   = (ch == 0x85) | (ch == 0xa0)
@inline _isalpha_l(c)    = ((0xc0 <= c <= 0xff) & (c != 0xf7) & (c != 0xd7)) | (c == 0xb5)
@inline _isalnum_l(c)    = _isalpha_l(c) || _isnumeric_l(c)
@inline _isprint_l(ch)   = ((0xa0 <= ch <= 0xff) & (ch != 0xad))
@inline _isgraph_l(ch)   = ((0xa0 < ch <= 0xff) & (ch != 0xad))

############################################################################
# Definitions for any Unicode codepoint (requires call to utf8proc) (only used for non-Latin1)

@inline _isnumeric_u(ch) = _check_mask(ch, _isnumeric_mask)
@inline _ispunct_u(ch)   = _check_mask(ch, _ispunct_mask)
@inline _isspace_u(ch)   = category_code(ch) == Uni.Zs
@inline _isalpha_u(ch)   = _check_mask(ch, _isalpha_mask)
@inline _isalnum_u(ch)   = _check_mask(ch, _isalnum_mask)
@inline _isprint_u(ch)   = _check_mask(ch, _isprint_mask)
@inline _isgraph_u(ch)   = _check_mask(ch, _isgraph_mask)

############################################################################
# Fallback definitions for all Chr types

@inline is_control(ch::CodeUnitTypes)   = _iscntrl(ch)
@inline is_digit(ch::CodeUnitTypes)     = _isdigit(ch)
@inline is_hex_digit(ch::CodeUnitTypes) = _isxdigit(ch)

@inline is_control(ch::Chr)   = is_control(codepoint(ch))
@inline is_digit(ch::Chr)     = is_digit(codepoint(ch))
@inline is_hex_digit(ch::Chr) = is_hex_digit(codepoint(ch))

@inline is_ascii(ch::Chr)          = is_ascii(codepoint(ch))
@inline is_ascii(ch::Unsigned)     = ch <= 0x7f
@inline is_ascii(ch::ASCIIChr)     = true

@inline is_latin(ch::AbstractChar) = is_latin(codepoint(ch))
@inline is_latin(ch::Unsigned)     = ch <= 0xff

@inline is_bmp(ch::AbstractChar)   = is_bmp(codepoint(ch))
@inline is_bmp(ch::Unsigned)       = ch <= 0xffff && !is_surrogate_codeunit(ch)
@inline is_bmp(ch::UInt8)          = true

@inline is_unicode(ch::AbstractChar) = is_unicode(codepoint(ch))
@inline is_unicode(ch::Unsigned)     = ch <= 0x10ffff && !is_surrogate_codeunit(ch)
@inline is_unicode(ch::UInt8)        = true

const _catfuns =
    ((:numeric,      :numeric),
     (:punctuation,  :punct),
     (:space,        :space),
     (:lowercase,    :lower),
     (:uppercase,    :upper),
     (:letter,       :alpha),
     (:alphanumeric, :alnum),
     (:printable,    :print),
     (:graphic,      :graph))

for (nnam, fnam) in _catfuns
    isnam   = Symbol(string("is_", nnam))
    flg = (fnam == :lower || fnam == :upper)
    namroot = string(flg ? "_is_" : "_is", fnam)
    fnam_a  = Symbol(string(namroot, "_a"))
    fnam_al = Symbol(string(namroot, "_al"))
    fnam_ch = Symbol(string(namroot, "_ch"))

    @eval $(fnam_al)(ch) = is_ascii(ch) ? $(fnam_a)(ch) : $(Symbol(string(namroot, "_l")))(ch)
    flg ||
        @eval $(fnam_ch)(ch) = is_latin(ch) ? $(fnam_al)(ch) : $(Symbol(string(namroot, "_u")))(ch)
    @eval $(isnam)(ch::CodeUnitTypes)  = $(fnam_ch)(ch)
    @eval $(isnam)(ch::Chr)            = $(fnam_ch)(codepoint(ch))
    @eval $(isnam)(ch::ASCIIChr)       = $(fnam_a)(codepoint(ch))
    @eval $(isnam)(ch::LatinChars)     = $(fnam_al)(codepoint(ch))
end

############################################################################

@static if isdefined(Base, :ismalformed)
    Base.ismalformed(ch::Chr) = false
    Base.isoverlong(ch::Chr) = false
    is_malformed(ch) = ismalformed(ch)
    is_overlong(ch) = isoverlong(ch)
else
    is_malformed(ch) = false
    is_overlong(ch) = false
end


text:
# Test typing in a typed gripper domain
path = joinpath(dirname(pathof(PDDL)), "..", "test", "typing")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == Symbol("gripper-typed")
@test domain.predicates[:free] == pddl"(free ?g)"
@test domain.predtypes[:carry] == [:ball, :gripper]
@test :gripper in keys(domain.types)

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.name == Symbol("gripper-problem")
@test problem.objects == @pddl("rooma", "roomb", "ball1", "ball2", "left", "right")
@test problem.objtypes[Const(:ball1)] == :ball

# Test forward execution of plans
state = init_state(problem)
state = execute(pddl"(pick ball1 rooma left)", state, domain)
@test state[pddl"(carry ball1 left)"] == true
state = execute(pddl"(move rooma roomb)", state, domain)
@test state[pddl"(robbyat roomb)"] == true
state = execute(pddl"(drop ball1 roomb left)", state, domain)
@test state[pddl"(at ball1 roomb)"] == true

@test satisfy(problem.goal, state)[1] == true

# Test action availability
state = init_state(problem)
@test Set(available(state, domain)) == Set{Term}(@pddl(
    "(pick ball1 rooma right)", "(pick ball1 rooma left)",
    "(pick ball2 rooma right)", "(pick ball2 rooma left)",
    "(move rooma roomb)", "(move rooma rooma)"
))

# Test backward regression of plans
state = goal_state(problem)
state = regress(pddl"(drop ball1 roomb left)", state, domain)
@test state[pddl"(carry ball1 left)"] == true
state = regress(pddl"(move rooma roomb)", state, domain)
@test state[pddl"(robbyat rooma)"] == true
state = regress(pddl"(pick ball1 rooma left)", state, domain)
@test state[pddl"(at ball1 rooma)"] == true
@test issubset(state, init_state(problem))

# Test action relevance
state = goal_state(problem)
@test Set(relevant(state, domain)) == Set{Term}(@pddl(
    "(drop ball1 roomb left)", "(drop ball1 roomb right)"
))

ext: jl


text:
module MechGlueModelingToolkit
import Unitfu
import Unitfu: AbstractQuantity
import ModelingToolkit
import ModelingToolkit: Num, Symbolic
import Base
import Base.*

Base.:*(x::Union{Num,Symbolic},y::T) where T<: AbstractQuantity = x * y

end


text:
import Optim
import Printf: @printf
import Random: shuffle!, randperm

const maxdegree = 2
const actualMaxsize = maxsize + maxdegree

# Sum of square error between two arrays
function SSE(x::Array{Float32}, y::Array{Float32})::Float32
    diff = (x - y)
    return sum(diff .* diff)
end
function SSE(x::Nothing, y::Array{Float32})::Float32
    return 1f9
end

# Sum of square error between two arrays, with weights
function SSE(x::Array{Float32}, y::Array{Float32}, w::Array{Float32})::Float32
    diff = (x - y)
    return sum(diff .* diff .* w)
end
function SSE(x::Nothing, y::Array{Float32}, w::Array{Float32})::Float32
    return Nothing
end

# Mean of square error between two arrays
function MSE(x::Nothing, y::Array{Float32})::Float32
    return 1f9
end

# Mean of square error between two arrays
function MSE(x::Array{Float32}, y::Array{Float32})::Float32
    return SSE(x, y)/size(x)[1]
end

# Mean of square error between two arrays
function MSE(x::Nothing, y::Array{Float32}, w::Array{Float32})::Float32
    return 1f9
end

# Mean of square error between two arrays
function MSE(x::Array{Float32}, y::Array{Float32}, w::Array{Float32})::Float32
    return SSE(x, y, w)/sum(w)
end

const len = size(X)[1]

if weighted
    const avgy = sum(y .* weights)/sum(weights)
    const baselineMSE = MSE(y, convert(Array{Float32, 1}, ones(len) .* avgy), weights)
else
    const avgy = sum(y)/len
    const baselineMSE = MSE(y, convert(Array{Float32, 1}, ones(len) .* avgy))
end

function id(x::Float32)::Float32
    x
end

const nuna = size(unaops)[1]
const nbin = size(binops)[1]
const nops = nuna + nbin
const nvar = size(X)[2];

function debug(verbosity, string...)
    verbosity > 0 ? println(string...) : nothing
end

function getTime()::Integer
    return round(Integer, 1e3*(time()-1.6e9))
end

# Define a serialization format for the symbolic equations:
mutable struct Node
    #Holds operators, variables, constants in a tree
    degree::Integer #0 for constant/variable, 1 for cos/sin, 2 for +/* etc.
    val::Union{Float32, Integer} #Either const value, or enumerates variable
    constant::Bool #false if variable
    op::Integer #enumerates operator (separately for degree=1,2)
    l::Union{Node, Nothing}
    r::Union{Node, Nothing}

    Node(val::Float32) = new(0, val, true, 1, nothing, nothing)
    Node(val::Integer) = new(0, val, false, 1, nothing, nothing)
    Node(op::Integer, l::Node) = new(1, 0.0f0, false, op, l, nothing)
    Node(op::Integer, l::Union{Float32, Integer}) = new(1, 0.0f0, false, op, Node(l), nothing)
    Node(op::Integer, l::Node, r::Node) = new(2, 0.0f0, false, op, l, r)

    #Allow to pass the leaf value without additional node call:
    Node(op::Integer, l::Union{Float32, Integer}, r::Node) = new(2, 0.0f0, false, op, Node(l), r)
    Node(op::Integer, l::Node, r::Union{Float32, Integer}) = new(2, 0.0f0, false, op, l, Node(r))
    Node(op::Integer, l::Union{Float32, Integer}, r::Union{Float32, Integer}) = new(2, 0.0f0, false, op, Node(l), Node(r))
end

# Copy an equation (faster than deepcopy)
function copyNode(tree::Node)::Node
   if tree.degree == 0
       return Node(tree.val)
   elseif tree.degree == 1
       return Node(tree.op, copyNode(tree.l))
    else
        return Node(tree.op, copyNode(tree.l), copyNode(tree.r))
   end
end

# Count the operators, constants, variables in an equation
function countNodes(tree::Node)::Integer
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + countNodes(tree.l)
    else
        return 1 + countNodes(tree.l) + countNodes(tree.r)
    end
end

# Count the max depth of a tree
function countDepth(tree::Node)::Integer
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + countDepth(tree.l)
    else
        return 1 + max(countDepth(tree.l), countDepth(tree.r))
    end
end

# Convert an equation to a string
function stringTree(tree::Node)::String
    if tree.degree == 0
        if tree.constant
            return string(tree.val)
        else
            if useVarMap
                return varMap[tree.val]
            else
                return "x$(tree.val - 1)"
            end
        end
    elseif tree.degree == 1
        return "$(unaops[tree.op])($(stringTree(tree.l)))"
    else
        return "$(binops[tree.op])($(stringTree(tree.l)), $(stringTree(tree.r)))"
    end
end

# Print an equation
function printTree(tree::Node)
    println(stringTree(tree))
end

# Return a random node from the tree
function randomNode(tree::Node)::Node
    if tree.degree == 0
        return tree
    end
    a = countNodes(tree)
    b = 0
    c = 0
    if tree.degree >= 1
        b = countNodes(tree.l)
    end
    if tree.degree == 2
        c = countNodes(tree.r)
    end

    i = rand(1:1+b+c)
    if i <= b
        return randomNode(tree.l)
    elseif i == b + 1
        return tree
    end

    return randomNode(tree.r)
end

# Count the number of unary operators in the equation
function countUnaryOperators(tree::Node)::Integer
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 1 + countUnaryOperators(tree.l)
    else
        return 0 + countUnaryOperators(tree.l) + countUnaryOperators(tree.r)
    end
end

# Count the number of binary operators in the equation
function countBinaryOperators(tree::Node)::Integer
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 0 + countBinaryOperators(tree.l)
    else
        return 1 + countBinaryOperators(tree.l) + countBinaryOperators(tree.r)
    end
end

# Count the number of operators in the equation
function countOperators(tree::Node)::Integer
    return countUnaryOperators(tree) + countBinaryOperators(tree)
end

# Randomly convert an operator into another one (binary->binary;
# unary->unary)
function mutateOperator(tree::Node)::Node
    if countOperators(tree) == 0
        return tree
    end
    node = randomNode(tree)
    while node.degree == 0
        node = randomNode(tree)
    end
    if node.degree == 1
        node.op = rand(1:length(unaops))
    else
        node.op = rand(1:length(binops))
    end
    return tree
end

# Count the number of constants in an equation
function countConstants(tree::Node)::Integer
    if tree.degree == 0
        return convert(Integer, tree.constant)
    elseif tree.degree == 1
        return 0 + countConstants(tree.l)
    else
        return 0 + countConstants(tree.l) + countConstants(tree.r)
    end
end

# Randomly perturb a constant
function mutateConstant(
        tree::Node, T::Float32,
        probNegate::Float32=0.01f0)::Node
    # T is between 0 and 1.

    if countConstants(tree) == 0
        return tree
    end
    node = randomNode(tree)
    while node.degree != 0 || node.constant == false
        node = randomNode(tree)
    end

    bottom = 0.1f0
    maxChange = perturbationFactor * T + 1.0f0 + bottom
    factor = maxChange^Float32(rand())
    makeConstBigger = rand() > 0.5

    if makeConstBigger
        node.val *= factor
    else
        node.val /= factor
    end

    if rand() > probNegate
        node.val *= -1
    end

    return tree
end

# Evaluate an equation over an array of datapoints
function evalTreeArray(tree::Node)::Union{Array{Float32, 1}, Nothing}
    return evalTreeArray(tree, X)
end

# Evaluate an equation over an array of datapoints
function evalTreeArray(tree::Node, cX::Array{Float32, 2})::Union{Array{Float32, 1}, Nothing}
    clen = size(cX)[1]
    if tree.degree == 0
        if tree.constant
            return fill(tree.val, clen)
        else
            return copy(cX[:, tree.val])
        end
    elseif tree.degree == 1
        cumulator = evalTreeArray(tree.l, cX)
        if cumulator === nothing
            return nothing
        end
        op_idx = tree.op
        UNAOP!(cumulator, op_idx, clen)
        @inbounds for i=1:clen
            if isinf(cumulator[i]) || isnan(cumulator[i])
                return nothing
            end
        end
        return cumulator
    else
        cumulator = evalTreeArray(tree.l, cX)
        if cumulator === nothing
            return nothing
        end
        array2 = evalTreeArray(tree.r, cX)
        if array2 === nothing
            return nothing
        end
        op_idx = tree.op
        BINOP!(cumulator, array2, op_idx, clen)
        @inbounds for i=1:clen
            if isinf(cumulator[i]) || isnan(cumulator[i])
                return nothing
            end
        end
        return cumulator
    end
end

# Score an equation
function scoreFunc(tree::Node)::Float32
    prediction = evalTreeArray(tree)
    if prediction === nothing
        return 1f9
    end
    if weighted
        mse = MSE(prediction, y, weights)
    else
        mse = MSE(prediction, y)
    end
    return mse / baselineMSE + countNodes(tree)*parsimony
end

# Score an equation with a small batch
function scoreFuncBatch(tree::Node)::Float32
    # batchSize
    batch_idx = randperm(len)[1:batchSize]
    batch_X = X[batch_idx, :]
    prediction = evalTreeArray(tree, batch_X)
    if prediction === nothing
        return 1f9
    end
    size_adjustment = 1f0
    batch_y = y[batch_idx]
    if weighted
        batch_w = weights[batch_idx]
        mse = MSE(prediction, batch_y, batch_w)
        size_adjustment = 1f0 * len / batchSize
    else
        mse = MSE(prediction, batch_y)
    end
    return size_adjustment * mse / baselineMSE + countNodes(tree)*parsimony
end

# Add a random unary/binary operation to the end of a tree
function appendRandomOp(tree::Node)::Node
    node = randomNode(tree)
    while node.degree != 0
        node = randomNode(tree)
    end

    choice = rand()
    makeNewBinOp = choice < nbin/nops
    if rand() > 0.5
        left = Float32(randn())
    else
        left = rand(1:nvar)
    end
    if rand() > 0.5
        right = Float32(randn())
    else
        right = rand(1:nvar)
    end

    if makeNewBinOp
        newnode = Node(
            rand(1:length(binops)),
            left,
            right
        )
    else
        newnode = Node(
            rand(1:length(unaops)),
            left
        )
    end
    node.l = newnode.l
    node.r = newnode.r
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.constant = newnode.constant
    return tree
end

# Insert random node
function insertRandomOp(tree::Node)::Node
    node = randomNode(tree)
    choice = rand()
    makeNewBinOp = choice < nbin/nops
    left = copyNode(node)

    if makeNewBinOp
        right = randomConstantNode()
        newnode = Node(
            rand(1:length(binops)),
            left,
            right
        )
    else
        newnode = Node(
            rand(1:length(unaops)),
            left
        )
    end
    node.l = newnode.l
    node.r = newnode.r
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.constant = newnode.constant
    return tree
end

# Add random node to the top of a tree
function prependRandomOp(tree::Node)::Node
    node = tree
    choice = rand()
    makeNewBinOp = choice < nbin/nops
    left = copyNode(tree)

    if makeNewBinOp
        right = randomConstantNode()
        newnode = Node(
            rand(1:length(binops)),
            left,
            right
        )
    else
        newnode = Node(
            rand(1:length(unaops)),
            left
        )
    end
    node.l = newnode.l
    node.r = newnode.r
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.constant = newnode.constant
    return node
end

function randomConstantNode()::Node
    if rand() > 0.5
        val = Float32(randn())
    else
        val = rand(1:nvar)
    end
    newnode = Node(val)
    return newnode
end

# Return a random node from the tree with parent
function randomNodeAndParent(tree::Node, parent::Union{Node, Nothing})::Tuple{Node, Union{Node, Nothing}}
    if tree.degree == 0
        return tree, parent
    end
    a = countNodes(tree)
    b = 0
    c = 0
    if tree.degree >= 1
        b = countNodes(tree.l)
    end
    if tree.degree == 2
        c = countNodes(tree.r)
    end

    i = rand(1:1+b+c)
    if i <= b
        return randomNodeAndParent(tree.l, tree)
    elseif i == b + 1
        return tree, parent
    end

    return randomNodeAndParent(tree.r, tree)
end

# Select a random node, and replace it an the subtree
# with a variable or constant
function deleteRandomOp(tree::Node)::Node
    node, parent = randomNodeAndParent(tree, nothing)
    isroot = (parent === nothing)

    if node.degree == 0
        # Replace with new constant
        newnode = randomConstantNode()
        node.l = newnode.l
        node.r = newnode.r
        node.op = newnode.op
        node.degree = newnode.degree
        node.val = newnode.val
        node.constant = newnode.constant
    elseif node.degree == 1
        # Join one of the children with the parent
        if isroot
            return node.l
        elseif parent.l == node
            parent.l = node.l
        else
            parent.r = node.l
        end
    else
        # Join one of the children with the parent
        if rand() < 0.5
            if isroot
                return node.l
            elseif parent.l == node
                parent.l = node.l
            else
                parent.r = node.l
            end
        else
            if isroot
                return node.r
            elseif parent.l == node
                parent.l = node.r
            else
                parent.r = node.r
            end
        end
    end
    return tree
end

# Simplify tree
function combineOperators(tree::Node)::Node
    # NOTE: (const (+*-) const) already accounted for. Call simplifyTree before.
    # ((const + var) + const) => (const + var)
    # ((const * var) * const) => (const * var)
    # ((const - var) - const) => (const - var)
    # (want to add anything commutative!)
    # TODO - need to combine plus/sub if they are both there.
    if tree.degree == 0
        return tree
    elseif tree.degree == 1
        tree.l = combineOperators(tree.l)
    elseif tree.degree == 2
        tree.l = combineOperators(tree.l)
        tree.r = combineOperators(tree.r)
    end

    top_level_constant = tree.degree == 2 && (tree.l.constant || tree.r.constant)
    if tree.degree == 2 && (binops[tree.op] === mult || binops[tree.op] === plus) && top_level_constant
        op = tree.op
        # Put the constant in r. Need to assume var in left for simplification assumption.
        if tree.l.constant
            tmp = tree.r
            tree.r = tree.l
            tree.l = tmp
        end
        topconstant = tree.r.val
        # Simplify down first
        below = tree.l
        if below.degree == 2 && below.op == op
            if below.l.constant
                tree = below
                tree.l.val = binops[op](tree.l.val, topconstant)
            elseif below.r.constant
                tree = below
                tree.r.val = binops[op](tree.r.val, topconstant)
            end
        end
    end

    if tree.degree == 2 && binops[tree.op] === sub && top_level_constant
        # Currently just simplifies subtraction. (can't assume both plus and sub are operators)
        # Not commutative, so use different op.
        if tree.l.constant
            if tree.r.degree == 2 && binops[tree.r.op] === sub
                if tree.r.l.constant
                    #(const - (const - var)) => (var - const)
                    l = tree.l
                    r = tree.r
                    simplified_const = -(l.val - r.l.val) #neg(sub(l.val, r.l.val))
                    tree.l = tree.r.r
                    tree.r = l
                    tree.r.val = simplified_const
                elseif tree.r.r.constant
                    #(const - (var - const)) => (const - var)
                    l = tree.l
                    r = tree.r
                    simplified_const = l.val + r.r.val #plus(l.val, r.r.val)
                    tree.r = tree.r.l
                    tree.l.val = simplified_const
                end
            end
        else #tree.r.constant is true
            if tree.l.degree == 2 && binops[tree.l.op] === sub
                if tree.l.l.constant
                    #((const - var) - const) => (const - var)
                    l = tree.l
                    r = tree.r
                    simplified_const = l.l.val - r.val#sub(l.l.val, r.val)
                    tree.r = tree.l.r
                    tree.l = r
                    tree.l.val = simplified_const
                elseif tree.l.r.constant
                    #((var - const) - const) => (var - const)
                    l = tree.l
                    r = tree.r
                    simplified_const = r.val + l.r.val #plus(r.val, l.r.val)
                    tree.l = tree.l.l
                    tree.r.val = simplified_const
                end
            end
        end
    end
    return tree
end

# Simplify tree
function simplifyTree(tree::Node)::Node
    if tree.degree == 1
        tree.l = simplifyTree(tree.l)
        if tree.l.degree == 0 && tree.l.constant
            return Node(unaops[tree.op](tree.l.val))
        end
    elseif tree.degree == 2
        tree.l = simplifyTree(tree.l)
        tree.r = simplifyTree(tree.r)
        constantsBelow = (
             tree.l.degree == 0 && tree.l.constant &&
             tree.r.degree == 0 && tree.r.constant
        )
        if constantsBelow
            return Node(binops[tree.op](tree.l.val, tree.r.val))
        end
    end
    return tree
end

# Define a member of population by equation, score, and age
mutable struct PopMember
    tree::Node
    score::Float32
    birth::Integer

    PopMember(t::Node) = new(t, scoreFunc(t), getTime())
    PopMember(t::Node, score::Float32) = new(t, score, getTime())

end

# Check if any binary operator are overly complex
function flagBinOperatorComplexity(tree::Node, op::Int)::Bool
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        return flagBinOperatorComplexity(tree.l, op)
    else
        if tree.op == op
            overly_complex = (
                    ((bin_constraints[op][1] > -1) &&
                     (countNodes(tree.l) > bin_constraints[op][1]))
                      ||
                    ((bin_constraints[op][2] > -1) &&
                     (countNodes(tree.r) > bin_constraints[op][2]))
                )
            if overly_complex
                return true
            end
        end
        return (flagBinOperatorComplexity(tree.l, op) || flagBinOperatorComplexity(tree.r, op))
    end
end

# Check if any unary operators are overly complex
function flagUnaOperatorComplexity(tree::Node, op::Int)::Bool
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        if tree.op == op
            overly_complex = (
                      (una_constraints[op] > -1) &&
                      (countNodes(tree.l) > una_constraints[op])
                )
            if overly_complex
                return true
            end
        end
        return flagUnaOperatorComplexity(tree.l, op)
    else
        return (flagUnaOperatorComplexity(tree.l, op) || flagUnaOperatorComplexity(tree.r, op))
    end
end

# Go through one simulated annealing mutation cycle
#  exp(-delta/T) defines probability of accepting a change
function iterate(member::PopMember, T::Float32, curmaxsize::Integer, frequencyComplexity::Array{Float32, 1})::PopMember
    prev = member.tree
    tree = prev
    #TODO - reconsider this
    if batching
        beforeLoss = scoreFuncBatch(prev)
    else
        beforeLoss = member.score
    end

    mutationChoice = rand()
    #More constants => more likely to do constant mutation
    weightAdjustmentMutateConstant = min(8, countConstants(prev))/8.0
    cur_weights = copy(mutationWeights) .* 1.0
    cur_weights[1] *= weightAdjustmentMutateConstant
    n = countNodes(prev)
    depth = countDepth(prev)

    # If equation too big, don't add new operators
    if n >= curmaxsize || depth >= maxdepth
        cur_weights[3] = 0.0
        cur_weights[4] = 0.0
    end
    cur_weights /= sum(cur_weights)
    cweights = cumsum(cur_weights)

    successful_mutation = false
    #TODO: Currently we dont take this \/ into account
    is_success_always_possible = true
    attempts = 0
    max_attempts = 10

    #############################################
    # Mutations
    #############################################
    while (!successful_mutation) && attempts < max_attempts
        tree = copyNode(prev)
        successful_mutation = true
        if mutationChoice < cweights[1]
            tree = mutateConstant(tree, T)

            is_success_always_possible = true
            # Mutating a constant shouldn't invalidate an already-valid function

        elseif mutationChoice < cweights[2]
            tree = mutateOperator(tree)

            is_success_always_possible = true
            # Can always mutate to the same operator

        elseif mutationChoice < cweights[3]
            if rand() < 0.5
                tree = appendRandomOp(tree)
            else
                tree = prependRandomOp(tree)
            end
            is_success_always_possible = false
            # Can potentially have a situation without success
        elseif mutationChoice < cweights[4]
            tree = insertRandomOp(tree)
            is_success_always_possible = false
        elseif mutationChoice < cweights[5]
            tree = deleteRandomOp(tree)
            is_success_always_possible = true
        elseif mutationChoice < cweights[6]
            tree = simplifyTree(tree) # Sometimes we simplify tree
            tree = combineOperators(tree) # See if repeated constants at outer levels
            return PopMember(tree, beforeLoss)

            is_success_always_possible = true
            # Simplification shouldn't hurt complexity; unless some non-symmetric constraint
            # to commutative operator...

        elseif mutationChoice < cweights[7]
            tree = genRandomTree(5) # Sometimes we generate a new tree completely tree

            is_success_always_possible = true
        else # no mutation applied
            return PopMember(tree, beforeLoss)
        end

        # Check for illegal equations
        for i=1:nbin
            if successful_mutation && flagBinOperatorComplexity(tree, i)
                successful_mutation = false
            end
        end
        for i=1:nuna
            if successful_mutation && flagUnaOperatorComplexity(tree, i)
                successful_mutation = false
            end
        end

        attempts += 1
    end
    #############################################

    if !successful_mutation
        return PopMember(copyNode(prev), beforeLoss)
    end

    if batching
        afterLoss = scoreFuncBatch(tree)
    else
        afterLoss = scoreFunc(tree)
    end

    if annealing
        delta = afterLoss - beforeLoss
        probChange = exp(-delta/(T*alpha))
        if useFrequency
            oldSize = countNodes(prev)
            newSize = countNodes(tree)
            probChange *= frequencyComplexity[oldSize] / frequencyComplexity[newSize]
        end

        return_unaltered = (isnan(afterLoss) || probChange < rand())
        if return_unaltered
            return PopMember(copyNode(prev), beforeLoss)
        end
    end
    return PopMember(tree, afterLoss)
end

# Create a random equation by appending random operators
function genRandomTree(length::Integer)::Node
    tree = Node(1.0f0)
    for i=1:length
        tree = appendRandomOp(tree)
    end
    return tree
end

# A list of members of the population, with easy constructors,
#  which allow for random generation of new populations
mutable struct Population
    members::Array{PopMember, 1}
    n::Integer

    Population(pop::Array{PopMember, 1}) = new(pop, size(pop)[1])
    Population(npop::Integer) = new([PopMember(genRandomTree(3)) for i=1:npop], npop)
    Population(npop::Integer, nlength::Integer) = new([PopMember(genRandomTree(nlength)) for i=1:npop], npop)

end

# Sample 10 random members of the population, and make a new one
function samplePop(pop::Population)::Population
    idx = rand(1:pop.n, ns)
    return Population(pop.members[idx])
end

# Sample the population, and get the best member from that sample
function bestOfSample(pop::Population)::PopMember
    sample = samplePop(pop)
    best_idx = argmin([sample.members[member].score for member=1:sample.n])
    return sample.members[best_idx]
end

function finalizeScores(pop::Population)::Population
    need_recalculate = batching
    if need_recalculate
        @inbounds @simd for member=1:pop.n
            pop.members[member].score = scoreFunc(pop.members[member].tree)
        end
    end
    return pop
end

# Return best 10 examples
function bestSubPop(pop::Population; topn::Integer=10)::Population
    best_idx = sortperm([pop.members[member].score for member=1:pop.n])
    return Population(pop.members[best_idx[1:topn]])
end

# Pass through the population several times, replacing the oldest
# with the fittest of a small subsample
function regEvolCycle(pop::Population, T::Float32, curmaxsize::Integer,
                      frequencyComplexity::Array{Float32, 1})::Population
    # Batch over each subsample. Can give 15% improvement in speed; probably moreso for large pops.
    # but is ultimately a different algorithm than regularized evolution, and might not be
    # as good.
    if fast_cycle
        shuffle!(pop.members)
        n_evol_cycles = round(Integer, pop.n/ns)
        babies = Array{PopMember}(undef, n_evol_cycles)

        # Iterate each ns-member sub-sample
        @inbounds Threads.@threads for i=1:n_evol_cycles
            best_score = Inf32
            best_idx = 1+(i-1)*ns
            # Calculate best member of the subsample:
            for sub_i=1+(i-1)*ns:i*ns
                if pop.members[sub_i].score < best_score
                    best_score = pop.members[sub_i].score
                    best_idx = sub_i
                end
            end
            allstar = pop.members[best_idx]
            babies[i] = iterate(allstar, T, curmaxsize, frequencyComplexity)
        end

        # Replace the n_evol_cycles-oldest members of each population
        @inbounds for i=1:n_evol_cycles
            oldest = argmin([pop.members[member].birth for member=1:pop.n])
            pop.members[oldest] = babies[i]
        end
    else
        for i=1:round(Integer, pop.n/ns)
            allstar = bestOfSample(pop)
            baby = iterate(allstar, T, curmaxsize, frequencyComplexity)
            #printTree(baby.tree)
            oldest = argmin([pop.members[member].birth for member=1:pop.n])
            pop.members[oldest] = baby
        end
    end

    return pop
end

# Cycle through regularized evolution many times,
# printing the fittest equation every 10% through
function run(
        pop::Population,
        ncycles::Integer,
        curmaxsize::Integer,
        frequencyComplexity::Array{Float32, 1};
        verbosity::Integer=0
       )::Population

    allT = LinRange(1.0f0, 0.0f0, ncycles)
    for iT in 1:size(allT)[1]
        if annealing
            pop = regEvolCycle(pop, allT[iT], curmaxsize, frequencyComplexity)
        else
            pop = regEvolCycle(pop, 1.0f0, curmaxsize, frequencyComplexity)
        end

        if verbosity > 0 && (iT % verbosity == 0)
            bestPops = bestSubPop(pop)
            bestCurScoreIdx = argmin([bestPops.members[member].score for member=1:bestPops.n])
            bestCurScore = bestPops.members[bestCurScoreIdx].score
            debug(verbosity, bestCurScore, " is the score for ", stringTree(bestPops.members[bestCurScoreIdx].tree))
        end
    end

    return pop
end

# Get all the constants from a tree
function getConstants(tree::Node)::Array{Float32, 1}
    if tree.degree == 0
        if tree.constant
            return [tree.val]
        else
            return Float32[]
        end
    elseif tree.degree == 1
        return getConstants(tree.l)
    else
        both = [getConstants(tree.l), getConstants(tree.r)]
        return [constant for subtree in both for constant in subtree]
    end
end

# Set all the constants inside a tree
function setConstants(tree::Node, constants::Array{Float32, 1})
    if tree.degree == 0
        if tree.constant
            tree.val = constants[1]
        end
    elseif tree.degree == 1
        setConstants(tree.l, constants)
    else
        numberLeft = countConstants(tree.l)
        setConstants(tree.l, constants)
        setConstants(tree.r, constants[numberLeft+1:end])
    end
end

# Proxy function for optimization
function optFunc(x::Array{Float32, 1}, tree::Node)::Float32
    setConstants(tree, x)
    return scoreFunc(tree)
end

# Use Nelder-Mead to optimize the constants in an equation
function optimizeConstants(member::PopMember)::PopMember
    nconst = countConstants(member.tree)
    if nconst == 0
        return member
    end
    x0 = getConstants(member.tree)
    f(x::Array{Float32,1})::Float32 = optFunc(x, member.tree)
    if size(x0)[1] == 1
        algorithm = Optim.Newton
    else
        algorithm = Optim.NelderMead
    end

    try
        result = Optim.optimize(f, x0, algorithm(), Optim.Options(iterations=100))
        # Try other initial conditions:
        for i=1:nrestarts
            tmpresult = Optim.optimize(f, x0 .* (1f0 .+ 5f-1*randn(Float32, size(x0)[1])), algorithm(), Optim.Options(iterations=100))
            if tmpresult.minimum < result.minimum
                result = tmpresult
            end
        end

        if Optim.converged(result)
            setConstants(member.tree, result.minimizer)
            member.score = convert(Float32, result.minimum)
            member.birth = getTime()
        else
            setConstants(member.tree, x0)
        end
    catch error
        # Fine if optimization encountered domain error, just return x0
        if isa(error, AssertionError)
            setConstants(member.tree, x0)
        else
            throw(error)
        end
    end
    return member
end

# List of the best members seen all time
mutable struct HallOfFame
    members::Array{PopMember, 1}
    exists::Array{Bool, 1} #Whether it has been set

    # Arranged by complexity - store one at each.
    HallOfFame() = new([PopMember(Node(1f0), 1f9) for i=1:actualMaxsize], [false for i=1:actualMaxsize])
end

# Check for errors before they happen
function testConfiguration()
    test_input = LinRange(-100f0, 100f0, 99)

    try
        for left in test_input
            for right in test_input
                for binop in binops
                    test_output = binop.(left, right)
                end
            end
            for unaop in unaops
                test_output = unaop.(left)
            end
        end
    catch error
        @printf("\n\nYour configuration is invalid - one of your operators is not well-defined over the real line.\n\n\n")
        throw(error)
    end
end

function fullRun(niterations::Integer;
                npop::Integer=300,
                ncyclesperiteration::Integer=3000,
                fractionReplaced::Float32=0.1f0,
                verbosity::Integer=0,
                topn::Integer=10
               )

    testConfiguration()

    # 1. Start a population on every process
    allPops = Future[]
    # Set up a channel to send finished populations back to head node
    channels = [RemoteChannel(1) for j=1:npopulations]
    bestSubPops = [Population(1) for j=1:npopulations]
    hallOfFame = HallOfFame()
    frequencyComplexity = ones(Float32, actualMaxsize)
    curmaxsize = 3
    if warmupMaxsize == 0
        curmaxsize = maxsize
    end

    for i=1:npopulations
        future = @spawnat :any Population(npop, 3)
        push!(allPops, future)
    end

    # # 2. Start the cycle on every process:
    @sync for i=1:npopulations
        @async allPops[i] = @spawnat :any run(fetch(allPops[i]), ncyclesperiteration, curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity), verbosity=verbosity)
    end
    println("Started!")
    cycles_complete = npopulations * niterations
    curmaxsize += 1

    last_print_time = time()
    num_equations = 0.0
    print_every_n_seconds = 5
    equation_speed = Float32[]

    for i=1:npopulations
        # Start listening for each population to finish:
        @async put!(channels[i], fetch(allPops[i]))
    end

    while cycles_complete > 0
        @inbounds for i=1:npopulations
            # Non-blocking check if a population is ready:
            if isready(channels[i])
                # Take the fetch operation from the channel since its ready
                cur_pop = take!(channels[i])
                bestSubPops[i] = bestSubPop(cur_pop, topn=topn)

                #Try normal copy...
                bestPops = Population([member for pop in bestSubPops for member in pop.members])

                for member in cur_pop.members
                    size = countNodes(member.tree)
                    frequencyComplexity[size] += 1
                    if member.score < hallOfFame.members[size].score
                        hallOfFame.members[size] = deepcopy(member)
                        hallOfFame.exists[size] = true
                    end
                end

                # Dominating pareto curve - must be better than all simpler equations
                dominating = PopMember[]
                open(hofFile, "w") do io
                    println(io,"Complexity|MSE|Equation")
                    for size=1:actualMaxsize
                        if hallOfFame.exists[size]
                            member = hallOfFame.members[size]
                            if weighted
                                curMSE = MSE(evalTreeArray(member.tree), y, weights)
                            else
                                curMSE = MSE(evalTreeArray(member.tree), y)
                            end
                            numberSmallerAndBetter = 0
                            for i=1:(size-1)
                                if weighted
                                    hofMSE = MSE(evalTreeArray(hallOfFame.members[i].tree), y, weights)
                                else
                                    hofMSE = MSE(evalTreeArray(hallOfFame.members[i].tree), y)
                                end
                                if (hallOfFame.exists[size] && curMSE > hofMSE)
                                    numberSmallerAndBetter += 1
                                end
                            end
                            betterThanAllSmaller = (numberSmallerAndBetter == 0)
                            if betterThanAllSmaller
                                println(io, "$size|$(curMSE)|$(stringTree(member.tree))")
                                push!(dominating, member)
                            end
                        end
                    end
                end
                cp(hofFile, hofFile*".bkup", force=true)

                # Try normal copy otherwise.
                if migration
                    for k in rand(1:npop, round(Integer, npop*fractionReplaced))
                        to_copy = rand(1:size(bestPops.members)[1])
                        cur_pop.members[k] = PopMember(
                            copyNode(bestPops.members[to_copy].tree),
                            bestPops.members[to_copy].score)
                    end
                end

                if hofMigration && size(dominating)[1] > 0
                    for k in rand(1:npop, round(Integer, npop*fractionReplacedHof))
                        # Copy in case one gets used twice
                        to_copy = rand(1:size(dominating)[1])
                        cur_pop.members[k] = PopMember(
                           copyNode(dominating[to_copy].tree)
                        )
                    end
                end

                @async begin
                    allPops[i] = @spawnat :any let
                        tmp_pop = run(cur_pop, ncyclesperiteration, curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity), verbosity=verbosity)
                        @inbounds @simd for j=1:tmp_pop.n
                            if rand() < 0.1
                                tmp_pop.members[j].tree = simplifyTree(tmp_pop.members[j].tree)
                                tmp_pop.members[j].tree = combineOperators(tmp_pop.members[j].tree)
                                if shouldOptimizeConstants
                                    tmp_pop.members[j] = optimizeConstants(tmp_pop.members[j])
                                end
                            end
                        end
                        tmp_pop = finalizeScores(tmp_pop)
                        tmp_pop
                    end
                    put!(channels[i], fetch(allPops[i]))
                end

                cycles_complete -= 1
                cycles_elapsed = npopulations * niterations - cycles_complete
                if warmupMaxsize != 0 && cycles_elapsed % warmupMaxsize == 0
                    curmaxsize += 1
                    if curmaxsize > maxsize
                        curmaxsize = maxsize
                    end
                end
                num_equations += ncyclesperiteration * npop / 10.0
            end
        end
        sleep(1e-3)
        elapsed = time() - last_print_time
        #Update if time has passed, and some new equations generated.
        if elapsed > print_every_n_seconds && num_equations > 0.0
            # Dominating pareto curve - must be better than all simpler equations
            current_speed = num_equations/elapsed
            average_over_m_measurements = 10 #for print_every...=5, this gives 50 second running average
            push!(equation_speed, current_speed)
            if length(equation_speed) > average_over_m_measurements
                deleteat!(equation_speed, 1)
            end
            average_speed = sum(equation_speed)/length(equation_speed)
            @printf("\n")
            @printf("Cycles per second: %.3e\n", round(average_speed, sigdigits=3))
            @printf("Hall of Fame:\n")
            @printf("-----------------------------------------\n")
            @printf("%-10s  %-8s   %-8s  %-8s\n", "Complexity", "MSE", "Score", "Equation")
            curMSE = baselineMSE
            @printf("%-10d  %-8.3e  %-8.3e  %-.f\n", 0, curMSE, 0f0, avgy)
            lastMSE = curMSE
            lastComplexity = 0

            for size=1:actualMaxsize
                if hallOfFame.exists[size]
                    member = hallOfFame.members[size]
                    if weighted
                        curMSE = MSE(evalTreeArray(member.tree), y, weights)
                    else
                        curMSE = MSE(evalTreeArray(member.tree), y)
                    end
                    numberSmallerAndBetter = 0
                    for i=1:(size-1)
                        if weighted
                            hofMSE = MSE(evalTreeArray(hallOfFame.members[i].tree), y, weights)
                        else
                            hofMSE = MSE(evalTreeArray(hallOfFame.members[i].tree), y)
                        end
                        if (hallOfFame.exists[size] && curMSE > hofMSE)
                            numberSmallerAndBetter += 1
                        end
                    end
                    betterThanAllSmaller = (numberSmallerAndBetter == 0)
                    if betterThanAllSmaller
                        delta_c = size - lastComplexity
                        delta_l_mse = log(curMSE/lastMSE)
                        score = convert(Float32, -delta_l_mse/delta_c)
                        @printf("%-10d  %-8.3e  %-8.3e  %-s\n" , size, curMSE, score, stringTree(member.tree))
                        lastMSE = curMSE
                        lastComplexity = size
                    end
                end
            end
            debug(verbosity, "")
            last_print_time = time()
            num_equations = 0.0
        end
    end
end


text:
using Knet, ArgParse, JLD
include("model.jl")
include("preprocess.jl")

function test(param, state, data; perp=false)
    totloss = 0.0
    numofbatch = 0
    for sequence in data
        val = loss(param, state, sequence, sequence)
        totloss += (perp? exp(val) : val)
        numofbatch += 1
    end
    return totloss / numofbatch
end

function update!(param, state, inputseq, outputseq; lr=1.0, gclip=0.0)
    gloss = lossgradient(param, state, inputseq, outputseq)
    gscale = lr

    if gclip > 0
        gnorm = sqrt(mapreduce(sumabs2, +, 0, gloss))
        if gnorm > gclip
            gscale *= gclip / gnorm
        end
    end

    for k=1:length(param)
        axpy!(-gscale, gloss[k], param[k])
    end

end

function train!(param, state, data, o)
    for sequence in data
        update!(param, state, sequence, sequence; lr=o[:lr], gclip=o[:gclip])
    end
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.exc_handler = ArgParse.debug_handler
    @add_arg_table s begin
        ("--trainfile"; default="readme_data.txt")
        ("--devfile"; help="Dev file")
        ("--testfile"; help="Test file")
        ("--loadfile"; help="Initialize model from file")
        ("--savefile"; help="Save final model to file")
        ("--vocabfile"; default=nothing; help="Vocabulary file")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--layerconfig"; arg_type=Int; nargs='+'; default=[64]; help="Sizes of the one or more LSTM layers")
        ("--batchsize"; arg_type=Int; default=1; help="Minibatch size")
        ("--encembed" ; arg_type=Int; default=128)
        ("--decembed" ; arg_type=Int; default=128)
        ("--epochs"; arg_type=Int; default=3; help="Number of epochs for training.")
        ("--winit"; arg_type=Float64; default=0.1; help="Initial weights set to winit*randn().")
        ("--decay"; arg_type=Float64; default=0.9; help="Learning rate decay.")
        ("--lr"; arg_type=Float64; default=4.0; help="Initial learning rate.")
        ("--gcheck"; arg_type=Int; default=0; help="Check N random gradients.")
        ("--gclip"; arg_type=Float64; default=0.0; help="Value to clip the gradient norm at.")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:atype] = eval(parse(o[:atype]))

    for (k,v) in o
        println("$k => $v")
    end

    # Traindata prepreration
    tdata = Data(o[:trainfile]; batchsize=o[:batchsize], vocabfile=o[:vocabfile])
    vsize = length(tdata.word_to_index)
    invocab = outvocab = vsize # vocabulary is the same for input and output

    # Devdata preperation
    dfile = (o[:devfile] != nothing ? o[:devfile] : o[:trainfile])
    ddata = Data(dfile; batchsize=o[:batchsize], word_to_index=tdata.word_to_index)

    # Model Initialization
    params = initparams(o[:atype], o[:layerconfig], o[:encembed], o[:decembed], invocab, outvocab, o[:winit])
    states = initstates(o[:atype], o[:layerconfig], o[:batchsize])

    # initial loss values
    inloss = test(params, states, tdata; perp=true)
    println("Initial loss is $inloss")
    devloss = test(params, states, ddata; perp=true)
    println("Initial devloss is $devloss")
    devlast = devbest = devloss

    # training started
    for epoch=1:o[:epochs]
        train!(params, states, tdata, o)
        devloss = test(params, states, ddata; perp=true)
        println("Dev loss for epoch $epoch: $devloss")

        if (epoch % 5) == 0
            trainloss = test(params, states, tdata; perp=true)
            println("\nTrain loss after epoch $epoch: $trainloss\n")
        end

        # check whether model becomes better
        if devloss < devbest
            devbest = devloss
            if o[:savefile] != nothing
                saveparam = map(p->convert(Array{Float32}, p), params)
                save(o[:savefile], "model", saveparam, "invocab", tdata.word_to_index, "config", o)
            end
        end

        if devloss > devlast
            o[:lr] *= o[:decay]
            info("New learning rate: $(o[:lr])")
        end
        devlast = devloss
        flush(STDOUT)
    end
end

!isinteractive() && main(ARGS)


text:
module SamplingTests

using FileIO
using Test

using QXContexts

include("utils.jl")

@testset "Sampling tests" begin

    test_path = dirname(@__DIR__)
    dsl_file = joinpath(test_path, "examples/ghz/ghz_5.qx")
    input_file = joinpath(test_path, "examples/ghz/ghz_5.jld2")
    param_file  = joinpath(test_path, "examples/ghz/ghz_5.yml")

    mktempdir() do path
        output_file = joinpath(path, "out.jld2")
        execute(dsl_file, input_file, param_file, output_file)

        # ensure all dictionary entries match
        output = load(output_file, "results")
        @test output[1] == collect(keys(ghz_results))
        @test output[2] ≈ collect(values(ghz_results))
    end

    # Test rejection sampling
    param_file = joinpath(test_path, "examples/ghz/ghz_5_rejection.yml")

    mktempdir() do path
        output_file = joinpath(path, "out.jld2")
        execute(dsl_file, input_file, param_file, output_file)

        # ensure all dictionary entries match
        output = FileIO.load(output_file, "results")
        @test length(output) == 10 # Should only have 10 samples
        @test length(unique(output)) == 2 # output should only contain strings "11111" and "00000"
    end

    # Test uniform sampling
    param_file = joinpath(test_path, "examples/ghz/ghz_5_uniform.yml")

    mktempdir() do path
        output_file = joinpath(path, "out.jld2")
        execute(dsl_file, input_file, param_file, output_file)

        # ensure all dictionary entries match
        output = FileIO.load(output_file, "results")
        @test length(output[1]) == 10 # Should only have 10 samples
        @test typeof(output[1][1]) == String
        @test length(output[2]) == 10 # should only have 10 amplitudes
        @test typeof(output[2][1]) <: Complex
    end
end

end


text:
function compute!(I::IntegralHelper{T,<:AbstractERI,<:AbstractRestrictedOrbitals}, entry::String) where T <: AbstractFloat
    # Create AO integral object
    basis = I.orbitals.basis
    aoorbs = AtomicOrbitals(I.molecule, basis)
    aoints = IntegralHelper{T}(molecule=I.molecule, orbitals=aoorbs, basis=basis, eri_type=I.eri_type)

    compute!(I, aoints, entry)
end

function compute!(I::IntegralHelper{T,<:AbstractERI,<:AbstractRestrictedOrbitals}, aoints::IntegralHelper{T,<:AbstractERI,AtomicOrbitals}, entry::String) where T<: AbstractFloat
    if entry == "S"
        compute_S!(I, aoints)
    elseif entry == "T"
        compute_T!(I, aoints)
    elseif entry == "V"
        compute_V!(I, aoints)
    elseif entry == "ERI"
        compute_ERI!(I, aoints)
    elseif entry == "OOOO"
        compute_OOOO!(I, aoints)
    elseif entry == "OOOV"
        compute_OOOV!(I, aoints)
    elseif entry == "OVOV"
        compute_OVOV!(I, aoints)
    elseif entry == "OOVV"
        compute_OOVV!(I, aoints)
    elseif entry == "OVVV"
        compute_OVVV!(I, aoints)
    elseif entry == "VVVV"
        compute_VVVV!(I, aoints)
    elseif occursin(r"F[dijab]{0,2}", entry)
        compute_F!(I, aoints)
    elseif I.eri_type === Chonky() # Out of options for Chonky
        throw(FermiException("Invalid key for IntegralHelper: $(entry)."))
    elseif entry == "BOO"
        compute_BOO!(I, aoints)
    elseif entry == "BOV"
        compute_BOV!(I, aoints)
    elseif entry == "BVV"
        compute_BVV!(I, aoints)
    end
end

include("OneElectron.jl")
include("DFERI.jl")
include("Chonky.jl")
include("Sparse.jl")
include("Fock.jl")

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
"""Boolean satisfiability problem in CNF form.
Like SATProblem but also has
- `clause_list`: `clause_list[i]` lists clauses containing variable `i`
"""

struct StructuredSATProblem
    num_variables::Int
    clauses::Vector{Vector{Int}}
    clause_list::Vector{Vector{Int}}
end

function StructuredSATProblem(p::SATProblem)

    clauses = sort(p.clauses, by=length)
    clause_list = list_clauses(p.num_variables, clauses)

    return StructuredSATProblem(p.num_variables, clauses, clause_list)
end

"Which clauses contain each variable"
function list_clauses(num_variables, clauses::Vector{Vector{Int}})
    clause_list = [Int[] for i in 1:num_variables]

    for (i, clause) in enumerate(clauses)
        for literal in clause
            push!(clause_list[index(literal)], i)
        end
    end

    return clause_list
end

solve(p::SATProblem, solver::SimpleSAT; kw...) =
    solve(StructuredSATProblem(p), solver; kw...)

function solve(p::StructuredSATProblem, solver::SimpleSAT; kw...)

    status, results = iterative_solve(p; kw...)

    if status == :unsat
        return :unsat, Int[]
    end

    return :sat, results
end

max_stars_repo_name: dpsanders/SimpleSATSolver

max_issues_repo_name: dpsanders/SimpleSATSolver

max_forks_repo_name: dpsanders/SimpleSATSolver


text:
export genofreq, genofreq_expected, genocount_observed, genocount_expected

"""
    genocount_observed(locus::GenoArray)
Return a `Dict` of genotype counts of a single locus in a
`PopData` object.
"""
@inline function genocount_observed(locus::T) where T<:GenoArray
    # conditional testing if all genos are missing
    all(ismissing.(locus)) && return missing
    countmap(skipmissing(locus))
end

"""
    genocount_expected(locus::GenoArray)
Return a `Dict` of the expected genotype counts of a single locus in a
`PopData` object. Expected counts are calculated as the product of observed
allele frequencies multiplied by the number of non-missing genotypes.
"""
function genocount_expected(locus::T) where T<:GenoArray
    #count number of non-missing genotypes in the locus
    n = nonmissing(locus)

    # Get expected number of genotypes in a locus
    ## get the observed allele frequencies
    allele_dict = allelefreq(locus)

    ## split the appropriate pairs into their own vectors
    alle, freq = collect(keys(allele_dict)), collect(values(allele_dict))
    ## calculate expected genotype frequencies by multiplying all-by-all x n
    expected_genotype_freq = vec(freq * freq' .* n)

    # reform genotype frequencies with same all-by-all approach
    genos = reverse.(Base.Iterators.product(alle, alle) |> collect |> vec)
    expected = Dict{nonmissingtype(eltype(locus)), Float64}()
    for (geno, freq) in zip(genos, expected_genotype_freq)
        expected[geno] = get!(expected, geno, 0.0) + freq
    end

    return expected
end

"""
    genofreq(locus::GenoArray)
Return a `Dict` of genotype frequencies of a single locus in a
`PopData` object.
"""
@inline function genofreq(locus::T) where T<:GenoArray
    # conditional testing if all genos are missing
    all(ismissing.(locus)) && return missing
    proportionmap(locus |> skipmissing |> collect)
end

"""
    genofreq(data::PopData, locus::String; population::Bool = false)
Return a `Dict` of genotype frequencies of a single locus in a `PopData`
object. Use `population = true` to return a table of genotype frequencies
of that locus per population.
### Example
```
cats = @nancycats;
genofreq(cats, "fca8")
genofreq(cats, "fca8", population = true)
```
"""
function genofreq(data::PopData, locus::String; population::Bool=false)
    if !population
        tmp = data.genodata[data.genodata.locus .== locus, :genotype]
        genofreq(tmp)
    else
        tmp = data.genodata[data.genodata.locus .== locus, :]
        DataFrames.combine(groupby(tmp, :population), :genotype => genofreq => :freq)
    end
end

"""
    genofreq_expected(locus::GenoArray)
Return a `Dict` of the expected genotype frequencies of a single locus in a
`PopData` object. Expected frequencies are calculated as the product of
observed allele frequencies.
"""
function genofreq_expected(locus::T) where T<:GenoArray
    # Get expected number of genotypes in a locus
    ## get the observed allele frequencies
    allele_dict = allelefreq(locus)

    ## split the appropriate pairs into their own vectors
    alle, freq = collect(keys(allele_dict)), collect(values(allele_dict))

    ## calculate expected genotype frequencies by multiplying all-by-all
    expected_genotype_freq = vec(freq * freq')

    # reform genotype frequencies with same all-by-all approach
    genos = reverse.(Base.Iterators.product(alle, alle) |> collect |> vec)

    # reform genotype frequencies into a Dict
    expected = Dict{nonmissingtype(eltype(locus)), Float64}()
    for (geno, freq) in zip(genos, expected_genotype_freq)
        expected[geno] = get!(expected, geno, 0.0) + freq
    end

    return expected
end

"""
    genofreq_expected(data::PopData, locus::String; population::Bool = false)
Return a `Dict` of expected genotype frequencies of a single locus in a
`PopData` object. Use `population = true` to return a table of expected genotype
frequencies of that locus per population.
### Example
```
cats = @nancycats;
genofreq_expected(cats, "fca8")
genofreq_expected(cats, "fca8", population = true)
```
"""
function genofreq_expected(data::PopData, locus::String; population::Bool=false)
    if !population
        tmp = data.genodata[data.genodata.locus .== locus, :genotype]
        genofreq_expected(tmp)
    else
        tmp = data.genodata[data.genodata.locus .== locus, :]
        DataFrames.combine(groupby(tmp, :population), :genotype => genofreq_expected => :freq)
    end
end


text:
# Autogenerated wrapper script for Yosys_jll for x86_64-apple-darwin
export libyosys, yosys, yosys_abc, yosys_config, yosys_filterlib, yosys_smtbmc

using boost_jll
using Readline_jll
using Tcl_jll
using Zlib_jll
using Libffi_jll
JLLWrappers.@generate_wrapper_header("Yosys")
JLLWrappers.@declare_library_product(libyosys, "@rpath/libyosys.so")
JLLWrappers.@declare_executable_product(yosys)
JLLWrappers.@declare_executable_product(yosys_abc)
JLLWrappers.@declare_executable_product(yosys_config)
JLLWrappers.@declare_executable_product(yosys_filterlib)
JLLWrappers.@declare_executable_product(yosys_smtbmc)
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, Readline_jll, Tcl_jll, Zlib_jll, Libffi_jll)
    JLLWrappers.@init_library_product(
        libyosys,
        "lib/libyosys.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        yosys,
        "bin/yosys",
    )

    JLLWrappers.@init_executable_product(
        yosys_abc,
        "bin/yosys-abc",
    )

    JLLWrappers.@init_executable_product(
        yosys_config,
        "bin/yosys-config",
    )

    JLLWrappers.@init_executable_product(
        yosys_filterlib,
        "bin/yosys-filterlib",
    )

    JLLWrappers.@init_executable_product(
        yosys_smtbmc,
        "bin/yosys-smtbmc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Lint.ExpressionUtils

@testset "Expressions" begin
    @test expand_trivial_calls(:(1:2)) == :(colon(1, 2))
    @test expand_trivial_calls(:(A')) == :(ctranspose(A))
    @test expand_trivial_calls(:(A.')) == :(transpose(A))
end


text:
function path_to_vec(path::Path, mapf::MAPF; T=nothing)
    g = mapf.graph
    edge_indices = mapf.edge_indices
    V, E = nv(g), ne(g)
    K = length(path)
    if isnothing(T)
        ind = Int[]
        val = Float64[]
        for k in 1:(K - 1)
            _, v1 = path[k]
            _, v2 = path[k + 1]
            e = edge_indices[(v1, v2)]
            i = findfirst(isequal(e), ind)
            if isnothing(i)
                push!(ind, e)
                push!(val, 1.)
            else
                val[i] += 1.
            end
        end
        return sparsevec(ind, val, E)
    else
        ind = Vector{Int}(undef, K - 1)
        for k in 1:(K - 1)
            t, v1 = path[k]
            _, v2 = path[k + 1]
            e = edge_indices[(v1, v2)]
            ind[k] = (t - 1) * E + e
        end
        val = ones(Float64, K - 1)
        return sparsevec(ind, val, T * E)
    end
end

function solution_to_mat(solution::Solution, mapf::MAPF; T=nothing)
    return reduce(hcat, path_to_vec(solution[a], mapf; T=T) for a in 1:nb_agents(mapf))
end

function solution_to_mat2(solution::Solution, mapf::MAPF; T=nothing)
    E = ne(mapf.graph)
    A = nb_agents(mapf)
    m = T * E
    n = A
    path_vecs = [path_to_vec(solution[a], mapf; T=T) for a in 1:A]
    rowval = reduce(vcat, pv.nzind for pv in path_vecs)
    nzval = reduce(vcat, pv.nzval for pv in path_vecs)
    colptr = reduce(vcat, fill(a, length(pv.nzval)) for (a, pv) in enumerate(path_vecs))
    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end


text:
# Read in the data

data = readdlm("meps1996.data")

#	The variables and their columns
#	VARIABLE						column
# 	Office based doctor visits	    1
# 	Outpatient doctor visits		2
# 	Emergency room visits			3
# 	Inpatient visits				4
# 	Dental visits					5
# 	Prescriptions					6
#  	PUBLIC_INS  					7
# 	PRIVATE_INS						8
# 	SEX								9
# 	AGE								10
# 	EDUC							11
# 	INCOME							12

names = ["OBDV", "OPV","IPV","ERV","DV","RX","PUB","PRIV","SEX","AGE","EDUC","INC/1000"]
data[:,12] = data[:,12]/1000.0
n = size(data,1)
println("MEPS data, 1996, complete data set statistics")
println("observations: ",n)
dstats(data, names);
""

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
import TPT

using Base.Test
using Plots; pyplot()

println("--- Unit HHTB ---")

!isdir("results") && mkdir("results")

resdir = joinpath("results", "unit_hhtb")
!isdir(resdir) && mkdir(resdir)

# Parameters (Hartree)
Z = 12
X = ["Ni", "Y"]
V = [10.04, 33.05] / (5.292^3 * 1e-3)
Nd = [8.60, 1.69]
Ed = [-4.93, -1.77] / 27.21
Wd = [3.78, 6.59] / 27.21
r₀ = (V * 3/4π) .^ (1/3)

x₂ = 0.5
x = [1-x₂, x₂]
Ec = mean(Ed)
Ed = Ed - Ec
Ec = 0

Emin = Ec - 0.2
Emax = Ec + 0.2

#
# Single-component case: Ni
#
hhtb1 = TPT.HHTB(Nd[1], Ed[1], Wd[1], r₀[1])

S = TPT.transfermatrix(hhtb1)[1,1]
plot(E -> imag(S(E)), Emin, Emax, xl="E (a.u.)", yl="S")
png(joinpath(resdir, "S_Ni"))

Δ = TPT.selfenergy(hhtb1)[1]
plot(E -> imag(Δ(E)), Emin, Emax, xl="E (a.u.)", yl="Δ")
png(joinpath(resdir, "Δ_Ni"))

Gᵢᵢ = TPT.diagonalgreenfunction(hhtb1)[1,1]
plot(E -> imag(Gᵢᵢ(E)), Emin, Emax, xl="E (a.u.)", yl="Im G_ii")
png(joinpath(resdir, "G_ii_Ni"))

G = TPT.offdiagonalgreenfunction(hhtb1)[1,1]
plot(E -> imag(G(E)), Emin, Emax, xl="E (a.u.)", yl="Im G_ij")
png(joinpath(resdir, "G_ij_Ni"))

D = TPT.densityofstate(hhtb1)
Ef = TPT.fermienergy(hhtb1)
plot(D, Emin, Emax, label="Ni", xl="E (a.u.)", yl="D (states / a.u. atom)")
vline!([Ef], label="Ef")
png(joinpath(resdir, "D_Ni"))

ΔNd = TPT.chargetransfer(hhtb1)[1]

Θ = TPT.bondorder(hhtb1)[1,1]

u = TPT.pairpotential(hhtb1)[1,1]
u′ = TPT.pairpotential_derivative(hhtb1)[1,1]
plot([u, u′], 2, 8, ylims=(-0.1, 0.1), xl="r (a.u.)", yl="u (a.u.)", label=["u" "u′"])
png(joinpath(resdir, "u_Ni"))

E_band = TPT.bandenergy(hhtb1)
E_site = TPT.onsiteenergy(hhtb1)

#
# Binary case: Ni-Y
#
hhtb = TPT.HHTB(x, Nd, Ed, Wd, r₀)

h = TPT.transferintegral(hhtb)
hd = TPT.nntransferintegral(hhtb)

default(xl="E (a.u.)")

S = TPT.transfermatrix(hhtb)
plot(E -> imag(S[1,1](E)), Emin, Emax, label=X[1])
plot!(E -> imag(S[2,2](E)), Emin, Emax, label=X[2])
ylabel!("S")
png(joinpath(resdir, "S"))

Δ = TPT.selfenergy(hhtb)
plot(E -> imag(Δ[1](E)), Emin, Emax)
plot!(E -> imag(Δ[2](E)), Emin, Emax)
ylabel!("Δ")
png(joinpath(resdir, "Δ"))

Gd = TPT.diagonalgreenfunction(hhtb)
plot(E -> imag(Gd[1](E)), Emin, Emax, label=X[1])
plot!(E -> imag(Gd[2](E)), Emin, Emax, label=X[2])
ylabel!("Im G_ii")
png(joinpath(resdir, "G_ii"))

G = TPT.offdiagonalgreenfunction(hhtb)
plot(E -> imag(G[1,1](E)), Emin, Emax, label="$(X[1])-$(X[1])")
plot!(E -> imag(G[1,2](E)), Emin, Emax, label="$(X[1])-$(X[2])")
plot!(E -> imag(G[2,2](E)), Emin, Emax, label="$(X[2])-$(X[2])")
ylabel!("Im G_ij")
png(joinpath(resdir, "G_ij"))

D = TPT.partialdensityofstate(hhtb)
Dtotal = TPT.totaldensityofstate(hhtb)
Ef = TPT.fermienergy(hhtb)
plot(E -> D[1](E), Emin, Emax, label=X[1])
plot!(E -> D[2](E), Emin, Emax, label=X[2])
plot!(E -> Dtotal(E), Emin, Emax, label="total")
vline!([Ef], label="Ef")
ylabel!("D (states / a.u. atom)")
png(joinpath(resdir, "D"))

ΔNd = TPT.chargetransfer(hhtb)

Θ = TPT.bondorder(hhtb)

u = TPT.pairpotential(hhtb)
plot([u[1,1], u[1,2], u[2,2]], 2, 8, labels=["Ni-Ni" "Ni-Y" "Y-Y"], ylims=(-0.1, 0.1), xl="r (a.u.)", yl="u (a.u.)")
png(joinpath(resdir, "u"))

E_band = TPT.bandenergy(hhtb)
E_site = TPT.onsiteenergy(hhtb)

@testset "Unit HHTB" begin
end


text:
using Plots, Random, Interact, Distributions, LinearAlgebra, DelimitedFiles

# N = 30 #number of osc
# plot([rand(Distributions.Binomial(30,0.1)) for i in 1:100], st=:histogram)
# rand([0,-rand(Distributions.Binomial(30,1/150))-0.0001*rand()])
function randomHessian(N; dims=2, p=0.0)
    p = p == 0.0 ? 1/(5*N) : p
    m = zeros(dims,N,N)
    for k in 1:dims
        # m[k,:,:] = [rand([0,-rand(1:N)-0.1*rand()]) for j in 1:N, i in 1:N]
        m[k,:,:] = [rand([0,-rand(Distributions.Binomial(N,p))-0.0001*rand()]) for j in 1:N, i in 1:N]
        m[k,:,:] = (m[k,:,:] + m[k,:,:]')/2
        for i in 1:N
            m[k,i,i] = -sum(m[k,i,:])+m[k,i,i]
        end
    end
    m ./ maximum(m, dims=3)
end

function expHessian(N; dims=2)
    m = zeros(dims,N,N)
    for k in 1:dims
        m[k,:,:] = [-exp(-(i-j)^2) for i in 1:N, j in 1:N]
        m[k,:,:] = (m[k,:,:] + m[k,:,:]')/2
        for i in 1:N
            m[k,i,i] = -sum(m[k,i,:])+m[k,i,i]
        end
    end
    m
end

function nNeighHessian(N; dims=2, n=2)
    m = zeros(dims,N,N)
    for k in 1:dims
        m[k,:,:] = [ abs(i-j) < n ? -1 : 0 for i in 1:N, j in 1:N]
        m[k,:,:] = (m[k,:,:] + m[k,:,:]')/2
        for i in 1:N
            m[k,i,i] = -sum(m[k,i,:])+m[k,i,i]
        end
    end
    m
end

function genHessian(N, type; dims=2, n=2)
    if type == "random"
        return randomHessian(N; dims)
    elseif type == "exp"
        return expHessian(N; dims)
    elseif type == "nNeigh"
        return nNeighHessian(N; dims,n)
    end
end

function gen_Cryst(N; L=[5,5,5])
    init_pos = zeros(3,N)
    for i in 0:N-1
        init_pos[1,i+1] = i % L[1]
        init_pos[2,i+1] = (i ÷ L[1]) % L[2]
        init_pos[3,i+1] = i ÷ (L[1]*L[2])
    end
    init_pos
end

function generateData(N; dims=2, steps=500, Δt=0.01, arrange = "Crystal",
            typeHessian="random", amp = 0.4, nModes = 0, idxStart = 1)
    if arrange == "Crystal"
        L = [Int(floor(N^(1/dims))) for i in 1:dims]
        if dims==1
            append!(L,1)
        end
        init_pos = gen_Cryst(N; L)
    end

    M = genHessian(N, typeHessian; dims)
    r = zeros(dims, N, steps)
    vel = zeros(dims, N, steps)
    nModes = nModes == 0 ? N : nModes
    for k in 1:dims
        d,v = eigen(M[k,:,:])
        # r[k,:,:] = hcat([sum([amp * d[i]^2 .* v[:,i] .* sin(2π*t*d[i]) for i in idxStart:nModes]) .+ init_pos[k,:] for t in 0:Δt:(steps-1)*Δt]...)
        r[k,:,:] = hcat([sum([amp .* v[:,i] .* sin(2π*t*d[i]) for i in idxStart:nModes]) .+ init_pos[k,:] for t in 0:Δt:(steps-1)*Δt]...)
        vel[k,:,:] = hcat([sum([amp .* v[:,i] .* cos(2π*t*d[i]) * 2π * d[i] for i in idxStart:nModes]) for t in 0:Δt:(steps-1)*Δt]...)
    end
    permutedims(r,(1,3,2)), permutedims(vel,(1,3,2)), M
end

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


text:
using ManifoldsBase, Documenter

makedocs(
    format = Documenter.HTML(prettyurls = false, assets = ["assets/favicon.ico"]),
    modules = [ManifoldsBase],
    authors = "Seth Axen, Mateusz Baran, Ronny Bergmann, and contributors.",
    sitename = "ManifoldsBase.jl",
    pages = [
        "Home" => "index.md",
        "How to define a manifold" => "example.md",
        "Design principles" => "design.md",
        "An abstract manifold" => "types.md",
        "Functions on maniolds" => [
            "Basic functions" => "functions.md",
            "Projections" => "projections.md",
            "Retractions" => "retractions.md",
            "Vector transports" => "vector_transports.md",
        ],
        "Manifolds" => "manifolds.md",
        "Decorating/Extending a Manifold" => "decorator.md",
        "Bases for tangent spaces" => "bases.md",
    ],
)
deploydocs(repo = "github.com/JuliaManifolds/ManifoldsBase.jl.git", push_preview = true)


text:
include("../../src/bfs_primitive.jl")
include("../../src/sparse_square_matrix.jl")

@testset "PrimitiveJulia" begin

### algorithm correctness test

    A = spzeros(Bool, 5, 5)
    A[2,1] = true
    A[1,3] = true
    A[4,3] = true
    A[3,4] = true
    A[4,5] = true
    A[4,2] = true

    sol = [1, 4, 2, 3, 4]
    @test bfs_primitive(A, 1) == sol

### exception throw tests

    # Input matrix is not square
    B = spzeros(Bool, 5, 6)
    @test_throws DimensionMismatch bfs_primitive(B, 1)

    # source is a not a positive Integer
    @test_throws ArgumentError bfs_primitive(A, -1)

    # source is not an Integer subtype
    @test_throws MethodError bfs_primitive(A, 1.2)

    # Input matrix is not a SparseMatrixCSC{Bool, U}
    C = spzeros(5, 5)
    @test_throws MethodError bfs_primitive(C, 1)

end


text:
using DataFrames

# Create test data
INCIDENCE_MAT = [0 1 0 1 1 0 0 0 0 0;
		 0 1 1 0 0 0 1 0 0 0;
		 1 1 0 1 1 0 0 0 1 1;
		 1 0 0 1 0 1 0 0 1 1;
		 0 0 1 1 0 0 0 1 0 0]

C_MAT = INCIDENCE_MAT'INCIDENCE_MAT

@testset "Similarity methods tests" begin

	# Check inputs:
	@test_throws MethodError similarity(DataFrame(INCIDENCE_MAT), method="jaccard", orientation="colwise") # non-matrix input
	@test_throws MethodError get_rca(INCIDENCE_MAT[:, 1], method="jaccard", orientation="colwise") # One-dim array
	@test_throws MethodError similarity(["1" "2"; "3" "4"], method="jaccard", orientation="colwise") # non-numeric matrix
	@test_throws AssertionError
	similarity(INCIDENCE_MAT, method="abc", orientation="colwise") # wrong method
	@test_throws MethodError
	similarity(INCIDENCE_MAT, method="jaccard", orientation="abc") # wrong orientation

	# Check output type:
	@test typeof(similarity(INCIDENCE_MAT, method="jaccard", orientation="colwise")) == Array{Float64, 2}
	@test typeof(similarity(INCIDENCE_MAT, method="jaccard", orientation="rowwise")) == Array{Float64, 2}
	@test typeof(similarity(INCIDENCE_MAT, method="association", orientation="colwise")) == Array{Float64, 2}

	# Check output result:
	# Methods
	@test similarity(INCIDENCE_MAT, method="association", orientation="colwise")[2, 4] == C_MAT[2, 4] / (C_MAT[2, 2] * C_MAT[4, 4])
	@test similarity(INCIDENCE_MAT, method="association", orientation="colwise")[9, 10] == C_MAT[9, 10] / (C_MAT[9, 9] * C_MAT[10, 10])
	@test similarity(INCIDENCE_MAT, method="cosine", orientation="colwise")[2, 4] == C_MAT[2, 4] / sqrt(C_MAT[2, 2] * C_MAT[4, 4])
	@test similarity(INCIDENCE_MAT, method="cosine", orientation="colwise")[9, 10] == C_MAT[9, 10] / sqrt(C_MAT[9, 9] * C_MAT[10, 10])
	@test similarity(INCIDENCE_MAT, method="inclusion", orientation="colwise")[2, 4] == C_MAT[2, 4] / min(C_MAT[2, 2], C_MAT[4, 4])
	@test similarity(INCIDENCE_MAT, method="inclusion", orientation="colwise")[9, 10] == C_MAT[9, 10] / min(C_MAT[9, 9], C_MAT[10, 10])
	@test similarity(INCIDENCE_MAT, method="jaccard", orientation="colwise")[2, 4] == C_MAT[2, 4] / (C_MAT[2, 2] + C_MAT[4, 4] - C_MAT[2, 4])
	@test similarity(INCIDENCE_MAT, method="jaccard", orientation="colwise")[9, 10] == C_MAT[9, 10] / (C_MAT[9, 9] + C_MAT[10, 10] - C_MAT[9, 10])
	@test similarity(INCIDENCE_MAT, method="proximity", orientation="colwise")[2, 4] == C_MAT[2, 4] / max(C_MAT[2, 2], C_MAT[4, 4])
	@test similarity(INCIDENCE_MAT, method="proximity", orientation="colwise")[9, 10] == C_MAT[9, 10] / max(C_MAT[9, 9], C_MAT[10, 10])

	# Orientation
	@test similarity(INCIDENCE_MAT, method="association", orientation="rowise") == similarity(INCIDENCE_MAT', method="association", orientation="colwise")
end

@testset "Density tests" begin

	# Define
	JACCARD_SIM = similarity(INCIDENCE_MAT, orientation="colwise", method = "jaccard")
	# Check inputs:
	@test_throws MethodError density(M = DataFrame(INCIDENCE_MAT), SIM = JACCARD_SIM) # INCIDENCE_MAT should be a matrix
	@test_throws MethodError density(M = INCIDENCE_MAT, SIM = DataFrame(JACCARD_SIM)) # SIM should be a matrix
	@test_throws AssertionError density(M = INCIDENCE_MAT, SIM = vcat(JACCARD_SIM, JACCARD_SIM)) # SIM should be square
	@test_throws AssertionError density(M = INCIDENCE_MAT, SIM = hcat(JACCARD_SIM, JACCARD_SIM[:, 1])) # number of cols in M should be == number of cols and rows in SIM

	# Check output:
	# Check that output is a matrix
	@test typeof(density(M = INCIDENCE_MAT, SIM = JACCARD_SIM)) == Array{Float64, 2}
	# Check that output is symmetrical
	@test size(density(M = INCIDENCE_MAT, SIM = JACCARD_SIM), 1) == size(density(M = INCIDENCE_MAT, SIM = JACCARD_SIM), 2)

	# Check output result:
	@test density(
		      M = INCIDENCE_MAT,
		      SIM = JACCARD_SIM
		      )[2, 5] == sum(INCIDENCE_MAT[2, :] .* JACCARD_SIM[5, :]) / sum(JACCARD_SIM[5,:])

	@test density(
		      M = INCIDENCE_MAT,
		      SIM = JACCARD_SIM
		      )[5, 9] == sum(INCIDENCE_MAT[5, :] .* JACCARD_SIM[9, :]) / sum(JACCARD_SIM[9, :])

end

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost


text:
using ControlSystems, ControlSystemIdentification

logspace(a, b) = exp10.(LinRange(a, b, 50))

"""
    bodeplot(G::FOTF)

Generate the bode plot of a given FOTF object.
"""
function ControlSystems.bodeplot(G::FOTF, w = exp10.(logspace(-4, 4)))
    H1 = freqresp(w.*im, G)
    result = zeros(ComplexF64, 1, 1, length(w))
    result[1, 1, :] = H1[:]
    H1 = FRD(w, result) # Frequency response data model here
    plot(H1)
end

"""
    bode(G::FOTF)

The frequency response data model of a given FOTF object.
"""
function ControlSystems.bode(G::FOTF, w = exp10.(logspace(-4, 4)))
    H1 = freqresp(w.*im, G)
    result = zeros(ComplexF64, 1, 1, length(w))
    result[1, 1, :] .= H1[:]
    H1 = FRD(H1, w) # Frequency response data model here
    return H1
end


text:
struct FileDependency
    path::AbstractPath
    params::NamedTuple
end

FileDependency(p::AbstractPath) = FileDependency(p, NamedTuple())

FileDependency(p::AbstractString, params) = FileDependency(Path(p), params)
FileDependency(p::AbstractString) = FileDependency(Path(p))

path(fdep::FileDependency) = fdep.path
params(fdep::FileDependency) = fdep.params

Base.getindex(fdep::FileDependency, i::Symbol) = getindex(params(fdep), i)

function Base.show(io::IO, fdep::FileDependency)
    print(io, """
        path: $(path(fdep))
            params: $(params(fdep))""")
end

struct FileDepGroup
    deps::Vector{FileDependency}
    sharedparams::NamedTuple
end

deps(fdg::FileDepGroup) = fdg.deps
params(fdg::FileDepGroup) = fdg.sharedparams
path(fdg::FileDepGroup) = path.(deps(fdg))

function Base.show(io::IO, fdg::FileDepGroup)
    print(io, """
        FileDepGroup with $(length(fdg.deps)) dependencies
            shared params: $(params(fdg))""")
end

# This is copied from https://github.com/JuliaLang/julia/pull/37299, not needed in julia 1.7+
function Base.keys(m::RegexMatch)
    idx_to_capture_name = PCRE.capture_names(m.regex.regex)
    return map(eachindex(m.captures)) do i
        # If the capture group is named, return it's name, else return it's index
        get(idx_to_capture_name, i, i)
    end
end

"""
    build_regex(pattern, groups)

A convenience function for generating `Regex` patterns
with named capture groups.
Spots to sub out should be surrounded by curly braces (eg. `{group_name}`)
and then the actual regular expressions as a `NamedTuple`.
For now, pass the regex as a raw string -
I'll figure out how to compose `Regex` objects another time
(or if you know, make a PR!).
eg `(; group_name=raw"\\w+)`.

Examples
≡≡≡≡≡≡≡≡≡≡

```jldoctest
julia> build_regex("{thing1}_{thing2}.txt", (thing1=raw"\\d+", thing2=raw"\\w+"))
r"(?<thing1>\\d+)_(?<thing2>\\w+).txt"
```
"""
function build_regex(pattern, groups)
    for (g, p) in pairs(groups)
        pattern = replace(pattern, string('{',g,'}')=> "(?<$g>$p)")
    end
    return Regex(pattern)
end

"""
    glob_pattern(dir, pattern; recursive=false)

Match files in directory `dir` with the regular expression `pattern`,
returning a vector of [`FileDependency`](@ref)s.
Use `recursive=true` to search for all files in subdirectories.

If `pattern` contins named capture groups
(eg using [`build_regex`](@ref)),
the returned `FileDependency`s will contain the matches for those capture groups
in the `params` field.

Examples
≡≡≡≡≡≡≡≡≡≡

```jldoctest
julia> rgx = build_regex("file{thing1}_{thing2}.txt", (thing1=raw"\\d+", thing2=raw"\\d+"));

julia> readdir("../test/projects/fileglob/")
4-element Vector{String}:
 "file1_1.txt"
 "file1_2.txt"
 "file2_1.txt"
 "file2_2.txt"

julia> glob_pattern("../test/projects/fileglob/", rgx)
4-element Vector{FileDependency}:
 path: ../test/projects/fileglob/file1_1.txt
    params: (thing1 = "1", thing2 = "1")
 path: ../test/projects/fileglob/file1_2.txt
    params: (thing1 = "1", thing2 = "2")
 path: ../test/projects/fileglob/file2_1.txt
    params: (thing1 = "2", thing2 = "1")
 path: ../test/projects/fileglob/file2_2.txt
    params: (thing1 = "2", thing2 = "2")
```
"""
function glob_pattern(dir, pattern; recursive=false)
    dir = Path(dir)
    files = filter(f-> isfile(f) && occursin(pattern, basename(f)), (recursive ? walkpath(dir) : readpath(dir)))
    deps = map(files) do file
        m = match(pattern, string(basename(file)))
        params =  NamedTuple(Symbol(k) => m[k] for k in keys(m))
        FileDependency(file, params)
    end

    return deps
end

"""
    groupdeps(deps, on)

Group a collection of [`FileDependency`](@ref)s based on matching captures groups.

Examples
≡≡≡≡≡≡≡≡≡≡

```jldoctest
julia> rgx = build_regex("file{thing1}_{thing2}.txt", (thing1=raw"\\d+", thing2=raw"\\d+"));

julia> glb = glob_pattern("../test/projects/fileglob/", rgx)
4-element Vector{FileDependency}:
 path: ../test/projects/fileglob/file1_1.txt
    params: (thing1 = "1", thing2 = "1")
 path: ../test/projects/fileglob/file1_2.txt
    params: (thing1 = "1", thing2 = "2")
 path: ../test/projects/fileglob/file2_1.txt
    params: (thing1 = "2", thing2 = "1")
 path: ../test/projects/fileglob/file2_2.txt
    params: (thing1 = "2", thing2 = "2")

julia> groupdeps(glb, [:thing1])
2-element Vector{FileDepGroup}:
 FileDepGroup with 2 dependencies
    shared params: (thing1 = "1",)
 FileDepGroup with 2 dependencies
    shared params: (thing1 = "2",)

```
"""
function groupdeps(deps, on)
    deps = copy(deps)
    depgroups = FileDepGroup[]
    while !isempty(deps)
        curdep = first(deps)
        idx = findall(deps) do d
            for grouper in on
                d[grouper] == curdep[grouper] || return false
            end
            return true
        end
        depgroup = splice!(deps, idx)
        sharedparms = NamedTuple(g=> params(curdep)[g] for g in on)
        push!(depgroups, FileDepGroup(depgroup, sharedparms))
    end
    return depgroups
end


text:
using AsyPlots

cd(@__DIR__)

function randpoint()
    p = rand()
    myrand() = 0.05 + 0.95rand()
    if p < 1/4
        Point(-myrand(),-myrand(),color = :Green)
    elseif p < 1/2
        Point(-myrand(),myrand(),color = :Blue)
    elseif p < 3/4
        Point(myrand(),-myrand(),color = :Blue)
    elseif p < 7/8
        Point(myrand()/2,myrand(),color = :Red)
    else
        Point(1/2 + myrand()/2, myrand(), color = :Green)
    end
end

P = Plot([randpoint() for _ in 1:1000]) + Path([-1.1 0; 1.1 0],arrow=Arrow(),color=:gray) +
    Path([0 -1.1; 0 1.1],arrow=Arrow(),color=:gray) + Label("\$x_1\$",(1.1,-0.1),color=:gray) + Label("\$x_2\$",(0.1,1.1),color=:gray)

save("easy-decision-tree-example.svg",P)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
__precompile__(true)

"""
The module for [Random123 Family](@ref).

Provide 8 RNG types:

- [`Threefry2x`](@ref)
- [`Threefry4x`](@ref)
- [`Philox2x`](@ref)
- [`Philox4x`](@ref)
- [`AESNI1x`](@ref)
- [`AESNI4x`](@ref)
- [`ARS1x`](@ref)
- [`ARS4x`](@ref)
"""
module Random123

export R123_USE_AESNI, set_counter!
include("common.jl")

export Threefry2x, Threefry4x
include("threefry.jl")

export Philox2x, Philox4x
include("philox.jl")

if R123_USE_AESNI
    export AESNI1x, AESNI4x
    include("aesni.jl")

    export ARS1x, ARS4x
    include("ars.jl")
end

end


text:
using Test
using Profile
using TAMode
using OrdinaryDiffEq
using LinearAlgebra
using ForwardDiff
using Turing

tps = 10.0 .^ range(-6.0, stop = 4.0, length = 25)

aboutZero = x -> isapprox(x, 0.0, rtol = 1.0e-5, atol = 1.0e-5)

include("testreactCode.jl")
include("testCompModel.jl")
include("testLS.jl")
include("testSample.jl")
include("testChains.jl")


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "N-(N-3)-N scenario" begin

@time for N in 5:6
    d = N-3

    filepath = "./data/quick_adjacency_decomposition/"
    filename = "$N-$d-$N"

    println("$filename scenario")

    PM = LocalSignaling(N,N,d)

    M = zeros(Int64,N,N)
    M[[1,2],1] .= 1

    for i in 1:(N-2)
        M[i+2,i] = 1
    end

    BG_seed = BellGame(M,d)

    vertices = LocalPolytope.vertices(PM)
    canonical_facets = LocalPolytope.adjacency_decomposition(vertices, BG_seed, PM, max_vertices = 100)
    games = collect(keys(canonical_facets))

    pretty_print_txt(games, filepath*filename)
    write_ieq(filepath*"ieq/"*filename, convert(IEQ, games))
end

end


text:
# Adaptation of some examples in MosekTools
# Refer to:
# https://github.com/JuliaOpt/MosekTools.jl/blob/master/test/jump_soc.jl
#
# MIT License
# Copyright (c) 2017, Ulf Worsøe, Mosek Aps

using JuMP

@testset "Second-order Cone Programming" begin
    solver = KNITRO.Optimizer

    @testset "SOC1" begin
        m = Model(optimizer_with_attributes(solver, "outlev" => 0))
        @variable(m, x)
        @variable(m, y)
        @variable(m, t >= 0.0)
        @objective(m, Min, t)
        @constraint(m, x + y >= 1.0)
        @constraint(m, [t,x,y] in MOI.SecondOrderCone(3))

        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.objective_value(m) ≈ sqrt(1/2) atol=1e-6

        @test JuMP.value.([x, y, t]) ≈ [0.5, 0.5, sqrt(1/2)] atol=1e-3

    end

    @testset "RotatedSOC1" begin
        m = Model(optimizer_with_attributes(solver, "outlev" => 0))
        @variable(m, x[1:5] >= 0.0)
        @variable(m, 0.0 <= u <= 5.0)
        @variable(m, v)
        @variable(m, t1 == 1.0)
        @variable(m, t2 == 1.0)

        @objective(m, Max, v)

        @constraint(m, [t1/sqrt(2),t2/sqrt(2),x...] in MOI.RotatedSecondOrderCone(7))
        @constraint(m, [x[1]/sqrt(2), u/sqrt(2), v] in MOI.RotatedSecondOrderCone(3))

        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.value.(x) ≈ [1,0,0,0,0] atol=1e-2
        @test JuMP.value(u) ≈ 5.0 atol=1e-4
        @test JuMP.value(v) ≈ sqrt(5.0) atol=1e-4
    end
end


text:
#########################
### global constants ####
#########################

#' Constants Used in the bigleaf Package
#'
#' This function defines the following constants:
#'
#' - cp           Specific heat of air for constant pressure (J K-1 kg-1)
#' - Rgas         Universal gas constant (J mol-1 K-1)
#' - Rv           Gas constant of water vapor (J kg-1 K-1) (Stull 1988 p.641)
#' - Rd           Gas constant of dry air (J kg-1 K-1) (Foken p. 245)
#' - Md           Molar mass of dry air (kg mol-1)
#' - Mw           Molar mass of water vapor (kg mol-1)
#' - eps          Ratio of the molecular weight of water vapor to dry air (=Mw/Md)
#' - g            Gravitational acceleration (m s-2)
#' - solar_constant Solar constant (W m-2)
#' - pressure0    Reference atmospheric pressure at sea level (Pa)
#' - Tair0        Reference air temperature (K)
#' - k            von Karman constant
#' - Cmol         Molar mass of carbon (kg mol-1)
#' - Omol         Molar mass of oxygen (kg mol-1)
#' - H2Omol       Molar mass of water (kg mol-1)
#' - sigma        Stefan-Boltzmann constant (W m-2 K-4)
#' - Pr           Prandtl number
#' - Sc_CO2       Schmidt number for CO2
#' - Kelvin       Conversion degree Celsius to Kelvin
#' - DwDc         Ratio of the molecular diffusivities for water vapor and CO2
#' - days2seconds Seconds per day
#' - kPa2Pa       Conversion kilopascal (kPa) to pascal (Pa)
#' - Pa2kPa       Conversion pascal (Pa) to kilopascal (kPa)
#' - umol2mol     Conversion micromole (umol) to mole (mol)
#' - mol2umol     Conversion mole (mol) to micromole (umol)
#' - kg2g         Conversion kilogram (kg) to gram (g)
#' - g2kg         Conversion gram (g) to kilogram (kg)
#' - kJ2J         Conversion kilojoule (kJ) to joule (J)
#' - J2kJ         Conversion joule (J) to kilojoule (kJ)
#' - se_median    Conversion standard error (SE) of the mean to SE of the median
#' - frac2percent Conversion between fraction and percent
#'
#' # Details
 This function is passed as an argument to every function that uses one
#'          or more constants. Individual constants passed to a function can be
#'          easily altered. E_g. the following command will change the value of
#'          the von Karman constant from 0.41 to 0.4:
#'
#'          `bigleaf_constants(k=0.4)`
#'
#'          the value of a constant can be returned by calling:
#'
#'          `bigleaf_constants$*name_of_constant*`
#'
#'          To permanently change the constants contained within this function (which
#'          makes sense for some of them, e.g. for the von Karman constant),
#'          the command `\link[utils]{fixInNamespace`} can be used. E_g.
#'
#'          `fixInNamespace(bigleaf_constants,ns="bigleaf")`
#'
#'          Note that this has to be repeated every time the package is newly installed/loaded.
#'
"""
"""
function bigleaf_constants(
  ## Physical constants
  cp         = 1004.834,        # specific heat of air for constant pressure (J K-1 kg-1)
  Rgas       = 8.31451,         # universal gas constant (J mol-1 K-1)
  Rv         = 461.5,           # gas constant of water vapor (J kg-1 K-1) (Stull 1988 p.641)
  Rd         = 287.0586,        # gas constant of dry air (J kg-1 K-1) (Foken 2008 p. 245)
  Md         = 0.0289645,       # molar mass of dry air (kg mol-1)
  Mw         = 0.0180153,       # molar mass of water vapor (kg mol-1)
  eps        = 0.622,           # ratio of the molecular weight of water vapor to dry air (=Mw/Md)
  g          = 9.81,            # gravitational acceleration (m s-2)
  solar_constant = 1366.1,      # solar constant, i.e. solar radation at earth distance from the sun (W m-2)
  pressure0  = 101325,          # reference atmospheric pressure at sea level (Pa)
  Tair0      = 273.15,          # reference air temperature (K)
  k          = 0.41,            # von Karman constant
  Cmol       = 0.012011,        # molar mass of carbon (kg mol-1)
  Omol       = 0.0159994,       # molar mass of oxygen (kg mol-1)
  H2Omol     = 0.01801528,      # molar mass of water (kg mol-1)
  sigma      = 5.670367e-08,    # Stefan-Boltzmann constant (W m-2 K-4)
  Pr         = 0.71,            # Prandtl number
  Sc_CO2     = 1.07,            # Schmidt number for CO2 (Hicks et al. 1987)

  ## Conversion constants
  Kelvin       = 273.15,         # conversion degree Celsius to Kelvin
  DwDc         = 1.6,            # Ratio of the molecular diffusivities for water vapor and CO2
  days2seconds = 86400,          # seconds per day
  kPa2Pa       = 1000,           # conversion kilopascal (kPa) to pascal (Pa)
  Pa2kPa       = 0.001,          # conversion pascal (Pa) to kilopascal (kPa)
  umol2mol     = 1e-06,          # conversion micromole (umol) to mole (mol)
  mol2umol     = 1e06,           # conversion mole (mol) to micromole (umol)
  kg2g         = 1000,           # conversion kilogram (kg) to gram (g)
  g2kg         = 0.001,          # conversion gram (g) to kilogram (kg)
  kJ2J         = 1000,           # conversion kilojoule (kJ) to joule (J)
  J2kJ         = 0.001,          # conversion joule (J) to kilojoule (kJ)
  se_median    = 1.253,          # conversion standard error (SE) of the mean to SE of the median (http://influentialpoints_com/Training/standard_error_of_median_htm)
  frac2percent = 100             # conversion between fraction and percent
)

  list(
    cp = cp, Rgas = Rgas, Rv = Rv, Rd = Rd, Md = Md, Mw = Mw, eps = eps, g = g,
    solar_constant = solar_constant,
    pressure0 = pressure0, Tair0 = Tair0, k = k, Cmol = Cmol, Omol = Omol,
    H2Omol = H2Omol,
    sigma = sigma, Pr = Pr, Sc_CO2 = Sc_CO2, Kelvin = Kelvin, DwDc = DwDc,
    days2seconds = days2seconds, kPa2Pa = kPa2Pa, Pa2kPa = Pa2kPa, umol2mol = umol2mol,
    mol2umol = mol2umol, kg2g = kg2g, g2kg = g2kg, kJ2J = kJ2J, J2kJ = J2kJ,
    se_median = se_median, frac2percent = frac2percent
  )

end


text:
module SurfaceTopology

using GeometryTypes

include("primitives.jl")
include("plainds.jl")
include("faceds.jl")
include("cachedds.jl")

include("edgeds.jl")

export FaceDS, CachedDS, EdgeDS
export FaceRing, VertexRing, EdgeRing
export Edges,Faces

end # module


text:
function testadd(x,y,z)
   a = x + y
   b = a + z
   c = b + a
   d = c + x
   return d
end

function testmul(x,y,z)
   a = x * y
   b = a * z
   c = z * x
   d = a * b
   return d
end

function testarith(x,y,z)
   a = x + y
   b = x * y
   c = z - b
   d = a / c
   return d
end

function testmuladd(x,y,z)
   a = muladd(x, y, z)
   b = muladd(z, x, y)
   c = muladd(z, b, -y)
   d = muladd(c, a,  x)
   return d
end

function testfma(x,y,z)
   a = fma(x, y, z)
   b = fma(z, x, y)
   c = fma(z, b, -y)
   d = fma(c, a,  x)
   return d
end

a32, a64 = Float32(pi)/7, Float64(pi)/7
b32, b64 = Float32(cbrt(3)), Float64(cbrt(3))
c32, c64 = Float32(sqrt(1/5)), Float64(sqrt(1/5))

x16, x32 = XFloat16(a32), XFloat32(a64)
y16, y32 = XFloat16(b32), XFloat32(b64)
z16, z32 = XFloat16(c32), XFloat32(c64)

ax16, ax32 = Float16(a32), Float32(a64)
by16, by32 = Float16(b32), Float32(b64)
cz16, cz32 = Float16(c32), Float32(c64)

xa32, xa64 = Float16(a32), Float32(a64)
yb32, yb64 = Float16(b32), Float32(b64)
zc32, zc64 = Float16(c32), Float32(c64)

relspeed_add32 =
  round( (@noelide @belapsed testadd($x16,$y16,$z16)) /
         (@noelide @belapsed testadd($a32,$b32,$c32)), digits=2);

relspeed_add64 =
 round( (@noelide @belapsed testadd($x32,$y32,$z32)) /
        (@noelide @belapsed testadd($a64,$b64,$c64)), digits=2);

relspeed_add =
 round( (@noelide @belapsed testadd($x32,$y32,$z32)) /
        (@noelide @belapsed testadd($ax32,$by32,$cz32)), digits=2);

relspeed_mul32 =
  round( (@noelide @belapsed testmul($x16,$y16,$z16)) /
         (@noelide @belapsed testmul($a32,$b32,$c32)), digits=2);

relspeed_mul64 =
 round( (@noelide @belapsed testmul($x32,$y32,$z32)) /
        (@noelide @belapsed testmul($a64,$b64,$c64)), digits=2);

relspeed_mul =
 round( (@noelide @belapsed testmul($x32,$y32,$z32)) /
        (@noelide @belapsed testmul($xa64,$yb64,$zc64)), digits=2);

relspeed_arith32 =
  round( (@noelide @belapsed testarith($x16,$y16,$z16)) /
         (@noelide @belapsed testarith($a32,$b32,$c32)), digits=2);

relspeed_arith64 =
 round( (@noelide @belapsed testarith($x32,$y32,$z32)) /
        (@noelide @belapsed testarith($a64,$b64,$c64)), digits=2);

relspeed_arith =
 round( (@noelide @belapsed testarith($x32,$y32,$z32)) /
        (@noelide @belapsed testarith($xa64,$yb64,$zc64)), digits=2);

relspeed_muladd32 =
  round( (@noelide @belapsed testmuladd($x16,$y16,$z16)) /
         (@noelide @belapsed testmuladd($a32,$b32,$c32)), digits=2);

relspeed_muladd64 =
 round( (@noelide @belapsed testmuladd($x32,$y32,$z32)) /
        (@noelide @belapsed testmuladd($a64,$b64,$c64)), digits=2);

relspeed_muladd =
 round( (@noelide @belapsed testmuladd($x32,$y32,$z32)) /
        (@noelide @belapsed testmuladd($xa64,$yb64,$zc64)), digits=2);

relspeed_fma32 =
  round( (@noelide @belapsed testfma($x16,$y16,$z16)) /
         (@noelide @belapsed testfma($a32,$b32,$c32)), digits=2);

relspeed_fma64 =
 round( (@noelide @belapsed testfma($x32,$y32,$z32)) /
        (@noelide @belapsed testfma($a64,$b64,$c64)), digits=2);

relspeed_fma =
 round( (@noelide @belapsed testfma($x32,$y32,$z32)) /
        (@noelide @belapsed testfma($xa64,$yb64,$zc64)), digits=2);

relspeeds = string(
   "\n\n\trelative speeds",
   "\n\t  (32)\t (64)\n\n",
   "add:   \t $relspeed_add32 \t $relspeed_add64 \n",
   "mul:   \t $relspeed_mul32 \t $relspeed_mul64 \n",
   "arith: \t $relspeed_arith32 \t $relspeed_arith64 \n",
   "muladd:\t $relspeed_muladd32 \t $relspeed_muladd64 \n",
   "fma:   \t $relspeed_fma32 \t $relspeed_fma64 \n");

print(relspeeds);

relspeeds = string(
   "\n\n\trelative speeds",
   "\n\t  (32)\t (64)\t (mixed)\n\n",
   "add:   \t $relspeed_add32 \t $relspeed_add64 \t $relspeed_add \n",
   "mul:   \t $relspeed_mul32 \t $relspeed_mul64 \t $relspeed_add  \n",
   "arith: \t $relspeed_arith32 \t $relspeed_arith64 \t $relspeed_add  \n",
   "muladd:\t $relspeed_muladd32 \t $relspeed_muladd64 \t $relspeed_add  \n",
   "fma:   \t $relspeed_fma32 \t $relspeed_fma64 \t $relspeed_fma \n");

print(relspeeds);


text:
using ITensors, ITensorNetworkAD
using ITensorNetworkAD.ITensorAutoHOOT: SubNetwork, get_leaves

@testset "test SubNetwork" begin
  i = Index(2, "i")
  j = Index(3, "j")
  k = Index(4, "k")
  l = Index(5, "l")
  A = randomITensor(i)
  B = randomITensor(j)
  C = randomITensor(k)
  D = randomITensor(l)

  AB = SubNetwork([A, B])
  ABCD = SubNetwork([AB, C, D])
  @test get_leaves(ABCD) == [A, B, C, D]
  @test inds(ABCD) == [i, j, k, l]
  @test noncommoninds(ABCD, AB, C) == [l]
end