text:
module test_firmware

# md""" -- markdown general discussion
macro md_mstr(s) println(s) end
# qc""" -- quote code -- run and print code -- don't print result
macro qc_mstr(s)
    quote
        println("```jl\n", strip(rstrip($s)), "\n```")
        $(esc(Base.parse("begin " * s * " end")))
    end
end

# qcr""" -- quote code repl -- run and print code and result like it was run in
# repl
macro qcr_mstr(s)
    quote
        oldstdout = STDOUT
        rd, wr = redirect_stdout()
        value = $(esc(Base.parse(s)))
        println()
        ioout = readavailable(rd)
#ioout=""
        redirect_stdout(oldstdout)
        print("```jl\n",
            replace(strip($s), r"^"m, "julia> "), "\n",
            rstrip(repr(value) * "\n" * ioout))
        println("\n```")
    end
end
#ERROR: readcb: bad file descriptor (EBADF)
# in wait_readnb at stream.jl:282
# in readavailable at stream.jl:658
# in reload_path at loading.jl:140
# in reload at loading.jl:73
#while loading /Users/m/s/m/current/Parsimonious.jl/test/test_firmware.jl, in expression starting on line 22

# TODO: qcd(command, test_return, test_stdout) -- doctest the command too.

md"""
Parsimonious.jl
===============

[![Build Status](https://travis-ci.org/gitfoxi/Parsimonious.jl.png)](https://travis-ci.org/gitfoxi/Parsimonious.jl)

A port of [Parsimonious](https://github.com/erikrose/parsimonious). Thanks, [Erik](https://github.com/erikrose), for composing an outstanding tool.

Porting Parsimonious to Julia is my Learn Julia project and it's been pretty interesting.

Let's Parse Some Dog Shit
-------------------------

"""

qc"""
using Parsimonious
import Parsimonious.visit
"""
# TODO: How to have less import stuff?

md"""
Sorry about `import Parsimonious.visit`. We're going to overload that function
    shortly but we also want some default functionality.

The key is to remain calm and organized at all times. Otherwise parser
development will get away from you. Our Firmware Command DSL was designed by
some fuckwit to communicate with an HP 93k tester. It looks like this:

    FTST?
    FTST P
    SREC? (@)
    SREC ACT,(@)
    SVLR TIM,PRM,,"Vdd",1.3; FTST?
    SHIT? "it's quoted",,(((FOO, BAR),(@@)),7,1.3),#9000000004asdf

We'd like to get these commands into structure like:
"""

qc"""
type FirmwareCommand
    command::String
    isquery::Bool
    parameters
end
"""

md"""
We build on the regular expressions that we know and love. First, let's just
try to parse:

    FTST?

Into:
"""

qc"""
FirmwareCommand("FTST",true,[])
"""

md"""
This is easily matched with a regex:
"""

qcr"""
match(r"(\w{4})(\??)", "FTST?")
"""

md"""
And almost as easily with a grammar:
"""

qc"""
g = grammar\"""
    statement = command isquery
    command = ~"\w{4}"
    isquery = "?"?
    \"""
"""

qcr"""
tree = parse(g, "FTST?")
"""

md"""
It's a parse tree. The goal is to fold it up.
"""

qc"""
type FwVis <: NodeVisitor end
"""

# TODO: qcr fails with multiple commands. Make it work like qc
# TODO: qcr doesn't indent output printed to STDOUT. Need to capture and indent or maybe just use ``` for now

qcr"""
debug_govisit(FwVis(), tree)
"""

md"""
**Sorry if you're reading this on Github you'll have to scroll the debug left
and right to see it. It's pretty wide. I don't know what to do about that.**

We haven't defined any visit methods but we get two for free. These fold the
tree up into `("FTST",("?",))` which is almost useful, but not really what
we want. Let's redefine those free rules so we can see what they do.
"""

qc"""
visit(::FwVis, n::LeafNode) = nodetext(n)
visit(::FwVis, n::ParentNode, visited_children...) = visited_children
"""

qcr"""
debug_govisit(FwVis(), tree)
"""

md"""
Works the same. The `LeafNode` rule just returns whatever text the leaf rule
matched. The `ParentNode` rule takes the text from one or more `LeafNode`s and
returns a list. The list is in order. That's important.

Now, let's put the visitor to work doing our bidding.
"""

qc"""
visit(::FwVis, n::ParentNode{:isquery}, questionmark) = questionmark == "?"
visit(::FwVis, n::ParentNode{:statement}, command::String, isquery::Bool) = FirmwareCommand(command, isquery, [])
"""

qcr"""
debug_govisit(FwVis(), tree)
"""

md"""
Cool. Now we're getting back `FirmwareCommand("FTST",true,[])` which is what we
wanted. Of course, you won't want to see the debug info when you're done:
"""

qcr"""
govisit(FwVis(), tree)
"""

md"""
Ya dig? Okay, now we add features to the language. The grammar should already
support:

    FTST

Does it?
"""

qcr"""
tree = parse(g, "FTST")
"""

md"""
It does. Let's visit:
"""

qcr"""
debug_govisit(FwVis(), tree)
"""

md"""
Ohs nos! What happened? Looking carefully we see that `:isquery` is no-longer a
ParentNode but a LeafNode because nothing matched below `isquery` (but
`isquery` matched the nothing). So we add a `visit`:
"""

qc"""
visit(::FwVis, n::LeafNode{:isquery}) = false
"""

md"""
No questionmark means no query. Let's try again.
"""

# TODO: doctest this one
qcr"""
govisit(FwVis(), tree)
"""

md"""
Great. How to recognize multiple statements? Noticing in the examples that
statements can be separated by ';' or '\n', let's add a rule to recognize
breaks between statements:

    termination = ~'[;\\n]'

And another one to recognize multiple statements:

    statements = statement+

So now the grammar looks like:
"""

qc"""
    g = grammar\"""
        statements = (statement termination)+
        statement = command isquery
        command = ~"\w{4}"
        isquery = "?"?
        termination = ~'[;\n]'
        \"""
"""

md"""
It's important to notice that the grammar now matches `statements` instead of
`statement` so our previous tests are now broken. Usually `parse` matches the
first statement in the list, but for testing purposes you can parse against any
of the statements. Let's check the old functionality:
"""

qcr"""
govisit(FwVis(), parse(g.exprs["statement"], "FTST?"))
"""

md"""
Good. And the other one:
"""

qcr"""
govisit(FwVis(), parse(g.exprs["statement"], "FTST"))
"""

qc"""
using Base.Test
"""
# TODO: qc above and below are pointlessly split up because apparently you
# can't say `using thing` and actually use the thing in the same quote.
qc"""
function fwtest(name, grammar, text)
    tree = parse(g.exprs[name], text)
    govisit(FwVis(), tree)
end

import Base.isequal
function isequal(a::FirmwareCommand, b::FirmwareCommand)
    a.command == b.command &&
    a.isquery == b.isquery &&
    a.parameters == b.parameters
end

@test fwtest("statement", g, "FTST?") == FirmwareCommand("FTST",true,[])
@test fwtest("statement", g, "FTST") == FirmwareCommand("FTST",false,[])
"""

md"""
So we haven't lost ground. Now let's see if we can string some statements
together:
"""

qcr"""
fwtest("statements", g, "FTST?;FTST?;FTST?\nFTST?\nFTST;FTST\n")
"""

md"""
Which almost looks good, but clearly we need to `visit` `statemets` differently
to throw away the useless terminations and just return a list of statements.
"""

qc"""
visit(::FwVis, n::ParentNode{:statements}, statements...) = [statement_term[1] for statement_term in statements]
"""

md"""
Try it again:
"""

qcr"""
fwtest("statements", g, "FTST?;FTST?;FTST?\nFTST?\nFTST;FTST\n")
"""

md"""
Nice. I wonder why sometimes list print in {} and other times in []. Oh well.

There's more kinds of statements. Each can have zero one or more parameters. To
match statements like:

    ASDF
    ASDF 1
    ASDF 1,2

We'll make three rules:

    params = more_params / one_param / no_params

This is where ordered choice really comes in handy. First it will try to match
several parameters, failing that it will try to match one and only failing that
can we get by with none. Of course if we wrote:

    params = no_params / one_param / more_params

Then we'd have a problem because no_params would always match and then we'd
never look for one_param or more_params and we'd get a parsing error. Moving
on:

Also, we have to start thinking about whitespace. If a command is followed by
parameters then there's at least one space separating them. But in the no
parameter case you don't need a space. How about:

    params = some_params / no_params
    some_params = somespace (more_params / one_param)
    somespace = ~'[ \t]+'
"""

qc"""
g = grammar\"""
    statements = (statement termination)+
    command = ~"\w{4}"
    isquery = "?"?
    termination = ~'[;\n]'

    ###### new stuff ##########################
    statement = command isquery params
    params = some_params / no_params
    some_params = somespace (more_params / param)
    more_params = param (comma param)+
    param = ~'[^ \t,;\n]*'
    no_params = ''
    comma = ','
    somespace = ~'[ \t]+'
    \"""
"""

md"""
I don't really want to think about it so let's just look at debug and go from
there.  """

qc"""
txt = "ASDF;ASDF 1;ASDF 1,2,3,4;"
tree = parse(g, txt)
"""
qcr"""
debug_govisit(FwVis(), tree)
"""

md"""
Wow. The debug output barely fits on my 13-inch laptop screen. I wish I could
come up with something more compact and just as useful, but this is addressing
several common problems:

* Parse not matching exactly the way I expected
* Wrong `visit` method handles a node

If you have a better idea, let me know.

`visit` time:
"""

qc"""
visit(::FwVis, n::LeafNode{:no_params}) = nothing
visit(::FwVis, n::LeafNode{:comma}) = nothing
visit(::FwVis, n::LeafNode{:somespace}) = nothing
visit(::FwVis, n::ParentNode{:statement}, command::String, isquery::Bool, params) = FirmwareCommand(command, isquery, params)
visit(::FwVis, n::ParentNode{:more_params}, param1::String, param2::String) = (param1, param2)
visit(::FwVis, n::ParentNode{:more_params}, param::String, params::Tuple) = tuple(param, params...)
visit(::FwVis, n::ParentNode, one_child) = one_child # don't rebox one child
"""

md"""
About now we realize, oh shit, things have gotten interesting. Julia's
multiple-dispatch feature is really working hard for us.

First, we've overloaded one of the most generic visit functions just for the case when there's only one child.

    visit(::FwVis, n::ParentNode, one_child) = one_child # don't rebox one child

Remember, the original (and still existing) generic rule for ParentNode is:

    visit(::FwVis, n::ParentNode, visited_children...) = visited_children

Which basically means that several children will get boxed up. `"a", "b", "c"
-> ("a", "b", "c")`. The problem is that there's these anonymous rules
everywhere and when they get called with one child, they put it in a box. So
after several iterations, you end up with `((("a",),),)`

Also fun is using the type to fix behavior. For example:

    visit(::FwVis, n::ParentNode{:more_params}, param1::String, param2::String) = (param1, param2)
    visit(::FwVis, n::ParentNode{:more_params}, param::String, params::Tuple) = tuple(param, params...)

We really want `more_params` to just give us a list of parameters, but since the rule is nested like so:

    more_params = param (comma param)+

`visit more_params` is sometimes being called with two individual parameters
and sometimes with a parameter and another `more_params`. Using the type system
this way we make sure that what comes out is always a flat list without
having to think about it further.

Finally, notice that several garbage tokens like `no_params`, `comma` and
`some_space` return `nothing`. When a `visit` returns `nothing` we don't even
send the `nothing` to the parent. It's just like it never happened which is
good for taking out the trash.

I'm a little nervous -- don't get me wrong -- relying so heavily on the
vagaries of Julia's dispatch system for the programs logic. But it looks like
it's working! Once you go hack you never go back.

Oh, and obviously I didn't crap this out fully formed. I wrote and tested each line
carefully refering back to the debug table.
"""

qcr"""
debug_govisit(FwVis(), tree)
"""

# A quick test to guard against regression.
@test fwtest("statements", g, txt) ==
    {FirmwareCommand("ASDF",false,()),FirmwareCommand("ASDF",false,"1"),FirmwareCommand("ASDF",false,("1","2","3","4"))}

pause()

g = grammar"""
    statements = (statement termination)+
    statement = command isquery (somespace params)?
    command = ~"\w{4}"
    isquery = "?"?
    termination = ~'[;\n]'
    params = more_params / one_param / no_params
    even_more_params = (comma param)+
    comma = ','
    more_params = param even_more_params
    one_param = param
    no_params = ''
    somespace = ~'[ \t]+'
    param = quoted / list / bare
    bare = ~'[^ \t,;\n")(]*'
    dq = '"'
    quoted = dq ~'[^"]*' dq
    open_paren = '('
    close_paren = ')'
    list = open_paren params close_paren
    """

# still works?
# Now param has to be reboxed
visit(::FwVis, n::ParentNode{:topparams}, params...) = params
visit(::FwVis, n::ParentNode{:param}, param...) = param[1]
println(parse(g, "ASDF;ASDF 1;ASDF 1,2;ASDF one;ASDF one,two;"))
@show fwtest("statements", g,  "ASDF;ASDF 1;ASDF 1,2;ASDF one;ASDF one,two;")
@test fwtest("statements", g,  "ASDF;ASDF 1;ASDF 1,2;ASDF one;ASDF one,two;") ==
        {FirmwareCommand("ASDF",false,[]),FirmwareCommand("ASDF",false,["1"]),FirmwareCommand("ASDF",false,{"1","2"}),FirmwareCommand("ASDF",false,["one"]),FirmwareCommand("ASDF",false,{"one","two"})}

@show fwtest("statements", g, """ASDF one,"two",(),(one),(one,"two"),(one,("two",(three)));""")

pause("Paused")

sample_command="FTST?"
# TODO: should be:
# command_grammar=Grammar(p"command = ~'\w{4}' '?'")
# rule_grammar is implicit
command_grammar=Grammar(rule_grammar, p"command = ~'\w{4}' '?'")
tree = parse(command_grammar, sample_command)
# tree = parse(command_grammar, sample_command)

sample_text =
"""

    FTST?
    FTST P
    SHIT 12,10;
    SHIT "poo poo"
    SHIT ()
    SQGB STSQ,30,(P2)
    SQGB? (@@);
    DFPN 10101,"1",(Clock)
    DFPN? (@)
    DFPN (11901, 11902),,(vee)

"""

# TODO: Print the partial parse tree on failure to match all of the text

# TODO: BUG: noparams has to be defined like this:
#    noparams = ~"[ \\t]*"
# If I define it like I want to then noparams doesn't appear in the parse tree
#    noparams = _

firmware_grammar_spec =
"""

    statements = statement+ whitespace
    statement = whitespace command isquery params termination
    command = ~"\\w{4}"
    more_params = (',' _ param)+
    paramlist = nonoptspace param _ more_params
    oneparam = nonoptspace nonemptyparam
    noparams = ~"[ \\t]*"
    params = paramlist / oneparam / noparams
    param = litparam / quotedparam / parenthesizedparam / emptyparam
    nonemptyparam = litparam / quotedparam / parenthesizedparam
    litparam = ~'[^")(;, \\t\\n]+'
    quotedparam = '"' ~'[^"]*' '"'
    parenthesizedparam = '(' params ')'
    emptyparam = _   # not so sure about this one
    termination = _ ~'[;\\n]'
    isquery = '?'?
    whitespace = ~"[ \\t\\n]*"
    _ = ~"[ \\t]*"
    nonoptspace = ~"[ \\t]+"

"""
# won't match:    termination = linespace (';' / '\\n')  # optional semicolon? seriously?

firmware_grammar = Grammar(rule_grammar, firmware_grammar_spec)
g = firmware_grammar

parse(g.exprs["command"], "FTST")
parse(g.exprs["statement"], "FTST;")
parse(g.exprs["termination"], ";")
parse(g.exprs["termination"], "\n")
# parse(g.exprs["simple_statement"], "FTST\n")
parse(g.exprs["statement"], "FTST\n")
parse(g.exprs["statement"], "FTST?;")
parse(g.exprs["statement"], "FTST P\n")
parse(g.exprs["statement"], "FTST P;")
parse(g.exprs["litparam"], "P")

# sample_text = "FTST?\n"
tree = parse(firmware_grammar, sample_text)
pprettily(tree)

# TODO: Doctest

# TODO: Refactor: LeafNode -> LeafNode
# TODO: Trace visitor

type FirmwareVisitor <: NodeVisitor end

function DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:statement}, visited_children)
    _, command, isquery, parameters, _ = visited_children
    FirmwareCommand(command, isquery, parameters)
end

DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:isquery}, visited_children) = visited_children[1] == "?"
DEFUNCT_visit(::FirmwareVisitor, n::LeafNode{:isquery}) = false

function DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:paramlist}, visited_children)
    _, param, _, more_params = visited_children
    @show param more_params
    vcat([param], more_params)
end

function DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:more_params}, visited_children)
    [vc[3] for vc in visited_children]
end

DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:param}, visited_children) = visited_children[1]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:nonemptyparam}, visited_children) = visited_children[1]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:statements}, visited_children) = visited_children[1]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:params}, visited_children) = visited_children[1]
DEFUNCT_visit(::FirmwareVisitor, n::LeafNode{:noparams}) = []
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:oneparam}, visited_children) = [visited_children[2]]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:litparam}, visited_children) = visited_children[1]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:quotedparam}, visited_children) = visited_children[2]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:parenthesizedparam}, visited_children) = visited_children[2]
DEFUNCT_visit(::FirmwareVisitor, n::ParentNode{:emptyparam}, visited_children) = []

DEFUNCT_visit(::FirmwareVisitor, n::ParentNode, visited_children) = visited_children
DEFUNCT_visit(::FirmwareVisitor, n::LeafNode, visited_children) = nodetext(n)
DEFUNCT_visit(::FirmwareVisitor, n::LeafNode) = nodetext(n)

firmware_statements = DEFUNCT_visit(FirmwareVisitor(), sample_text, tree)
# TODO: make keyword arguments work right
# firmware_statements = DEFUNCT_visit(FirmwareVisitor(), sample_text, tree; debug=true)

@show firmware_statements

md"""
PSA: When you learn that it is fun and easy to write parsers you will be
tempted to invent a complex, heirarchical data format -- or worse a DSL.  Some
dude on Reddit told me DSL stands for Dog Shit Language. As a professional who
both works with DSLs and picks up dog shit on reg, I must say I find the
comparison a little unfair to dog shit. The answer is almost always JSON, XML,
[config] or CSV for data and some established language for whatever you are
trying to express.

If you're writing parsers to rid the world of entrenched DSLs or arbitrary data
formats then you're doing it right. Thank you for your attention.
"""

end


text:
# Test 2D modeling
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

using JUDI
using ArgParse, Test, Printf, Aqua
using SegyIO, LinearAlgebra, Distributed, JOLI
using TimerOutputs: TimerOutputs, @timeit

set_verbosity(false)
# Collect timing and allocations information to show in a clear way.
const TIMEROUTPUT = TimerOutputs.TimerOutput()
timeit_include(path::AbstractString) = @timeit TIMEROUTPUT path include(path)

# Utilities
const success_log = Dict(true => "SUCCESS", false => "FAILED")
# Test set
const GROUP = get(ENV, "GROUP", "JUDI")

include("utils.jl")

parsed_args = parse_commandline()
const nlayer = parsed_args["nlayer"]
const tti = parsed_args["tti"]
const fs =  parsed_args["fs"]
const nw = parsed_args["parallel"]
const viscoacoustic = parsed_args["viscoacoustic"]

base = ["test_geometry.jl",
        "test_judiVector.jl",
        "test_composite.jl",
        "test_judiWeights.jl",
        "test_judiWavefield.jl",
        "test_linear_operators.jl",
        "test_physicalparam.jl",
        "test_compat.jl"]

devito = ["test_linearity.jl",
          "test_adjoint.jl",
          "test_all_options.jl",
          "test_jacobian.jl",
          "test_gradient_fwi.jl",
          "test_gradient_lsrtm.jl",
          "test_multi_exp.jl"]
        #   "test_gradient_twri.jl"]

extras = ["test_modeling.jl", "test_basics.jl", "test_linear_algebra.jl"]

issues = ["test_issues.jl"]

# custom
if endswith(GROUP, ".jl")
    # VERSION >= v"1.7" && push!(Base.ARGS, "-p 2")
    timeit_include(GROUP)
end

# Basic JUDI objects tests, no Devito
if GROUP == "JUDI" || GROUP == "All"
    for t=base
        timeit_include(t)
        try Base.GC.gc(); catch; gc() end
    end
    # Test resolved issues Due to some data type incomaptibilities only test 1.7
    if VERSION >= v"1.7"
        for t=issues
            timeit_include(t)
            try Base.GC.gc(); catch; gc() end
        end
    end
end

# Generic mdeling tests
if GROUP == "BASICS" || GROUP == "All"
    println("JUDI generic modelling tests")
    VERSION >= v"1.7" && push!(Base.ARGS, "-p 2")
    for t=extras
        timeit_include(t)
        @everywhere try Base.GC.gc(); catch; gc() end
    end
end

# Isotropic Acoustic tests
if GROUP == "ISO_OP" || GROUP == "All"
    println("JUDI iso-acoustic operators tests (parallel)")
    # Basic test of LA/CG/LSQR needs
    VERSION >= v"1.7" && push!(Base.ARGS, "-p 2")
    for t=devito
        timeit_include(t)
        @everywhere try Base.GC.gc(); catch; gc() end
    end
end

# Isotropic Acoustic tests with a free surface
if GROUP == "ISO_OP_FS" || GROUP == "All"
    println("JUDI iso-acoustic operators with free surface tests")
    push!(Base.ARGS, "--fs")
    for t=devito
        timeit_include(t)
        try Base.GC.gc(); catch; gc() end
    end
end

# Anisotropic Acoustic tests
if GROUP == "TTI_OP" || GROUP == "All"
    println("JUDI TTI operators tests")
    # TTI tests
    push!(Base.ARGS, "--tti")
    for t=devito
        timeit_include(t)
        try Base.GC.gc(); catch; gc() end
    end
end

# Anisotropic Acoustic tests with free surface
if GROUP == "TTI_OP_FS" || GROUP == "All"
    println("JUDI TTI operators with free surface tests")
    push!(Base.ARGS, "--tti")
    push!(Base.ARGS, "--fs")
    for t=devito
        timeit_include(t)
        try Base.GC.gc(); catch; gc() end
    end
end

# Viscoacoustic tests
if GROUP == "VISCO_AC_OP" || GROUP == "All"
    println("JUDI Viscoacoustic operators tests")
    # Viscoacoustic tests
    push!(Base.ARGS, "--viscoacoustic")
    for t=devito
        timeit_include(t)
        try Base.GC.gc(); catch; gc() end
    end
end

# Code quality check
if GROUP == "AQUA" || GROUP == "All" || GROUP == "JUDI"
    @testset "code quality" begin
        # Prevent ambiguities from PyCall and other packages
        Aqua.test_all(JUDI; ambiguities=false)
        Aqua.test_ambiguities(JUDI; Aqua.askwargs(true)...)
    end
end

# Testing memory and runtime summary
show(TIMEROUTPUT; compact=true, sortby=:firstexec)

ext: jl


text:
type Plot
    data::DataFrame
    aes::Aesthetic
    layers::Vector{AbstractLayer}
end

function ggplot(df::DataFrame, aes::Aesthetic)
    return Plot(df, aes, [])
end

function Base.(:+)(plot::Plot, layer::AbstractLayer)
    layer.data = plot.data
    layer.aes = plot.aes
    push!(plot.layers, layer)
    return plot
end


text:
using Clustering

@everywhere function compute_nmi(z1::Array{Int64,1}, z2::Array{Int64,1})

    n = length(z1)
    k1 = length(unique(z1))
    k2 = length(unique(z2))

    nk1 = zeros(k1,1)
    nk2 = zeros(k2,1)

    for (idx,val) in enumerate(unique(z1))
        cluster_idx = find(z1->z1 == val, z1)
        nk1[idx] = length(cluster_idx)
    end

    for (idx,val) in enumerate(unique(z2))
        cluster_idx = find(z2->z2 == val, z2)
        nk2[idx] = length(cluster_idx)
    end

    pk1 = nk1/float(sum(nk1))
    pk2 = nk2/float(sum(nk2))

    nk12 = zeros(k1,k2)
    for (idx1, val1) in enumerate(unique(z1))
        for (idx2, val2) in enumerate(unique(z2))
            cluster_idx1 = find(z1->z1 == val1, z1)
            cluster_idx2 = find(z2->z2 == val2, z2)
            common_idx12 = intersect(Set(cluster_idx1), Set(cluster_idx2))
            nk12[idx1,idx2] = length(common_idx12)
        end
    end
    pk12 = nk12/float(n)

    Hx = -sum(pk1.*log(pk1+eps(Float64)))
    Hy = -sum(pk2.*log(pk2+eps(Float64)))
    Hxy = -sum(pk12.*log(pk12+eps(Float64)))

    MI = Hx + Hy - Hxy
    nmi = MI/float(0.5*(Hx+Hy))

    nmi
end

@everywhere function cluster_eval(z1::Array{Int64,1}, z2::Array{Int64,1})

    K1 = length(unique(z1))
    K2 = length(unique(z2))

    #normalized MI
    nmi = compute_nmi(z1,z2)
    println("nmi: ", nmi)

    #variation of information
    vi = varinfo(K1,z1,K2,z2)
    println("vi: ", vi)

    #adjusted rand index
    ari = randindex(z1,z2)
    println("ari: ", ari[1])

end

#z1 = [1,1,1,2,2,2]
#z2 = [1,1,2,2,3,3]
#cluster_eval(z1, z2)


text:
using Images, StaticArrays

@testset "haar" begin

    # Basic Haar coordinates tests (:x2)
    test_coordinates = haar_coordinates(2,2,:x2)
    correct_coordinates = [SMatrix{4, 2}(1,1,1,1,1,2,1,2)', SMatrix{4, 2}(2,1,2,1,2,2,2,2)']
    @test test_coordinates == correct_coordinates

    test_coordinates = haar_coordinates(3,2,:x2)
    correct_coordinates = [SMatrix{4, 2}(1,1,1,1,1,2,1,2)', SMatrix{4, 2}(1,1,2,1,1,2,2,2)', SMatrix{4, 2}(2,1,2,1,2,2,2,2)', SMatrix{4, 2}(2,1,3,1,2,2,3,2)', SMatrix{4, 2}(3,1,3,1,3,2,3,2)']
    @test test_coordinates == correct_coordinates

    # Basic Haar coordinates tests (:y2)
    test_coordinates = haar_coordinates(2,2,:y2)
    correct_coordinates = [SMatrix{4, 2}(1,1,1,1,2,1,2,1)', SMatrix{4, 2}(1,2,1,2,2,2,2,2)']
    @test test_coordinates == correct_coordinates

    test_coordinates = haar_coordinates(3,2,:y2)
    correct_coordinates = [SMatrix{4, 2}(1,1,1,1,2,1,2,1)', SMatrix{4, 2}(1,2,1,2,2,2,2,2)', SMatrix{4, 2}(2,1,2,1,3,1,3,1)', SMatrix{4, 2}(2,2,2,2,3,2,3,2)']
    @test test_coordinates == correct_coordinates

    # Basic Haar coordinates tests (:x3)
    test_coordinates = haar_coordinates(3,2,:x3)
    @test length(test_coordinates) == 0

    test_coordinates = haar_coordinates(3,3,:x3)
    correct_coordinates = [SMatrix{4, 3}(1,1,1,1,1,2,1,2,1,3,1,3)', SMatrix{4, 3}(1,1,2,1,1,2,2,2,1,3,2,3)', SMatrix{4, 3}(2,1,2,1,2,2,2,2,2,3,2,3)', SMatrix{4, 3}(2,1,3,1,2,2,3,2,2,3,3,3)', SMatrix{4, 3}(3,1,3,1,3,2,3,2,3,3,3,3)']
    @test test_coordinates == correct_coordinates

    # Basic Haar coordinates tests (:y3)
    test_coordinates = haar_coordinates(2,3,:y3)
    @test length(test_coordinates) == 0

    test_coordinates = haar_coordinates(3,3,:y3)
    correct_coordinates = [SMatrix{4, 3}(1,1,1,1,2,1,2,1,3,1,3,1)', SMatrix{4, 3}(1,1,1,2,2,1,2,2,3,1,3,2)', SMatrix{4, 3}(1,2,1,2,2,2,2,2,3,2,3,2)', SMatrix{4, 3}(1,2,1,3,2,2,2,3,3,2,3,3)', SMatrix{4, 3}(1,3,1,3,2,3,2,3,3,3,3,3)']
    @test test_coordinates == correct_coordinates

    # Basic Haar coordinates tests (:xy4)
    test_coordinates = haar_coordinates(2,2,:xy4)
    correct_coordinates = [SMatrix{4, 4}(1,1,1,1,1,2,1,2,2,2,2,2,2,1,2,1)']
    @test test_coordinates == correct_coordinates

    test_coordinates = haar_coordinates(3,2,:xy4)
    correct_coordinates = [SMatrix{4, 4}(1,1,1,1,1,2,1,2,2,2,2,2,2,1,2,1)', SMatrix{4, 4}(2,1,2,1,2,2,2,2,3,2,3,2,3,1,3,1)']
    @test test_coordinates == correct_coordinates

    test_coordinates = haar_coordinates(3,3,:xy4)
    correct_coordinates = [SMatrix{4, 4}(1,1,1,1,1,2,1,2,2,2,2,2,2,1,2,1)', SMatrix{4, 4}(1,2,1,2,1,3,1,3,2,3,2,3,2,2,2,2)', SMatrix{4, 4}(2,1,2,1,2,2,2,2,3,2,3,2,3,1,3,1)', SMatrix{4, 4}(2,2,2,2,2,3,2,3,3,3,3,3,3,2,3,2)']
    @test test_coordinates == correct_coordinates

    # Simple Image haar_features

    img = diagm([1,1,1,1,1])
    int_img = integral_image(img)

    test_features = haar_features(int_img, [2,2], [4,4], :x2)
    correct_features = [-1,0,0,-1,1,1,-1,0,0,1]
    @test test_features == correct_features

    test_features = haar_features(int_img, [1,1], [4,4], :xy4)
    correct_features = [-2,0,0,-4,1,-2,0,0,1,-2,-2,1,0,0,1,-2]
    @test test_features == correct_features

    # Simple Image haar_features with coordinates given

    coordinates = rand(Int,3,3,4)
    coordinates[1,1,:] = [1,3,1,3]
    coordinates[1,2,:] = [1,4,1,4]
    coordinates[1,3,:] = [1,5,1,5]
    coordinates[2,1,:] = [3,2,3,2]
    coordinates[2,2,:] = [3,3,3,3]
    coordinates[2,3,:] = [3,4,3,4]
    coordinates[3,1,:] = [5,1,5,1]
    coordinates[3,2,:] = [5,2,5,2]
    coordinates[3,3,:] = [5,3,5,3]
    test_features = haar_features(int_img, [1,1], [5,5], :x3, coordinates)
    correct_features = [0,1,0]
    @test test_features == correct_features

    # Tests for error checks

    @test_throws ArgumentError haar_features(int_img, [1,1,1], [5,5], :x3)
    @test_throws ArgumentError haar_features(int_img, [7,1], [5,5], :x3)
    @test_throws ArgumentError haar_features(int_img, [1,7], [5,5], :x3)
    @test_throws ArgumentError haar_features(int_img, [1,1], [5,5], :x4)

    #No feature coordinates found case
    @test length(haar_features(int_img, [1,1], [1,5], :x2)) == 0
end


text:
#=
Parse an RSS feed
=#
module Graze
import HTTPClient: get
import LibExpat: xp_parse, find

type Feed
    attributes :: Dict
    items :: Array{Dict}
end

# Exceptions
type FeedError <: Exception end
type HTTPError <: Exception end
type NotImplementedError <: Exception end

# Maybe
abstract Option

type Nothing <: Option
    reason :: String
    exception :: Exception
end

type Something <: Option
    value
end

function pipeline(fn)
    #=
    A function for handling multiple layers of failable operations.
    =#
    return (maybe :: Option) -> (
        if isa(maybe, Nothing)
            maybe
        else
            fn(maybe.value)
        end
    )
end

function resolver(maybe :: Option)
    #=
    Either bet the value from Something, or raise the error contained
    in a Nothing.
    =#
    if isa(maybe, Nothing)
        info(maybe.reason)
        throw(maybe.exception)
    else
        maybe.value
    end
end

# Parsing
function parse(url)
    #=
    Parse an RSS feed.

    input:
        url: A URL to an RSS feed.
    output:
        Feed
    =#
    response = get(url)

    # feed not found
    if response.http_code != 200
        info("Invalid status code $(response.http_code)")
        throw(HTTPError())
    end

    parse_feed(bytestring(response.body))
end

function parse_feed(xml_string)
    #=
    Parse an XML document containing an RSS feed.

    input:
        xml_string: A string containing an XML document.
    output:
        Feed
    =#
    root = null

    try
        root = xp_parse(xml_string)
    catch
        info("Invalid XML file")
        throw(FeedError())
    end

    channel = (
        if root.name == "rss"
            Something(root)
        else
            Nothing("Invalid feed type: $(root.name)", FeedError())
        end
    ) |> pipeline((root) -> (
        version = get(root.attr, "version", "");

        if version == "2.0"
            Something(root)
        else
            Nothing("Unknown RSS version: $version", FeedError())
        end
    )) |> pipeline((root) -> (
        channels = find(root, "channel");

        if length(channels) < 1
            Nothing("Feed missing its channel", FeedError())
        else
            Something(channels[1])
        end
    )) |> resolver

    attributes = Dict()
    items = Dict[]

    for element in channel.elements
        if element.name == "item"
            # Story
            item = Dict()

            for item_element in element.elements
                item[item_element.name] = item_element.elements[1]
            end

            push!(items, item)
        else
            attributes[element.name] = element.elements[1]
        end
    end

    return Feed(attributes, items)
end

export parse, FeedError, HTTPError, NotImplementedError
end


text:
module Crypto

using MbedTLS

hmacsha256_digest(s::String, k::Union{String, Vector{UInt8}}) = digest(MD_SHA256, s, k)
export hmacsha256_digest

hmacsha1_digest(s::String, k::Union{String, Vector{UInt8}}) = digest(MD_SHA1, s, k)
export hmacsha1_digest

function md5_file(s::String)
    f = open(s)
    md = nothing
    try
        md = md5(f)
    finally
        close(f)
    end
    md
end
export md5_file

md5(s::String) = digest(MD_MD5, s)

function md5(s::IO)
    md = MbedTLS.MD(MD_MD5)
    success = false
    try
        while !eof(s)
            b = read(s, UInt8, min(nb_available(s), 65536))    # Read in 64 K chunks....
            MbedTLS.write(md, b)
        end
        d = MbedTLS.finish!(md)
        success = true
        return d
    finally
        success || MbedTLS.finish!(md)
    end
end
export md5

sha256(s::String) = digest(MD_SHA256, s)

end # Module end


text:
# Conversions between Julia and q types for the Q module.
#########################################################################
const K_None = K_Other(K_new(nothing))

Base.convert(::Type{String}, x::K_symbol) =
    unsafe_string(unsafe_load(pointer(x)))
Base.convert(::Type{Symbol}, x::K_char) = Symbol(Char(x))
Base.convert(::Type{String}, x::K_char) = unsafe_string(pointer(x), 1)

for T in (Int8, Int16, Int32, Int64, Int128, Float32, Float64)
    for S in (_Bool, _Signed, _Unsigned, _Float)
        @eval Base.convert(::Type{$T}, x::$S) = $T(x.a[])
    end
end

# Julia to K conversions
Base.convert(::Type{K}, ::Void) = K_None
function Base.convert(::Type{K}, x::K_)
    @assert x != C_NULL
    t = xt(x)
    if t == -EE
       	msg = xs(x)
        r0(x)
        throw(KdbException(msg))
    end
    if t < 0
        return K_Atom(x)
    elseif 0 <= t <= KT
        return K_Vector(x)
    elseif t == XT
        return K_Table(x)
    elseif t == XT && xt(xx(x)) == xt(xy(x))
    elseif t in (100, 104, 105, 112)
        return K_Lambda(x)
    end
    return K_Other(x)
end
Base.convert(::Type{K}, x) = convert(K, K_new(x))

Base.print(io::IO, x::K_boolean) = print(io, Bool(x) ? "1b" : "0b")
Base.show(io::IO, x::K_boolean) = print(io, x)
Base.print(io::IO, x::K_symbol) = print(io, Symbol(x))

"""
A helper method to convert error returns from k(..) or ee(..)
to a Julia exception.
"""
function _E(x)
    if xt(x) == -EE
        try
            throw(KdbException(xs(x)))
        finally
            r0(x)
        end
    end
    K(x)
end

# Numeric promotions
Base.convert(::Type{Float64}, x::K_real) = Float64(x.a[])


text:
module VML

import Base: .^, ./

# TODO detect CPU architecture
const lib = :libmkl_vml_avx
Libdl.dlopen(:libmkl_rt)

immutable VMLAccuracy
    mode::UInt
end
const VML_LA = VMLAccuracy(0x00000001)
const VML_HA = VMLAccuracy(0x00000002)
const VML_EP = VMLAccuracy(0x00000003)
Base.show(io::IO, m::VMLAccuracy) = print(io, m == VML_LA ? "VML_LA" :
                                              m == VML_HA ? "VML_HA" : "VML_EP")
vml_get_mode() = ccall((:_vmlGetMode, lib), Cuint, ())
vml_set_mode(mode::Integer) = (ccall((:_vmlSetMode, lib), Cuint, (UInt,), mode); nothing)

vml_set_accuracy(m::VMLAccuracy) = vml_set_mode((vml_get_mode() & ~0x03) | m.mode)
vml_get_accuracy() = VMLAccuracy(vml_get_mode() & 0x3)

vml_set_mode((vml_get_mode() & ~0x0000FF00))
function vml_check_error()
    vml_error = ccall((:_vmlClearErrStatus, lib), Cint, ())
    if vml_error != 0
        if vml_error == 1
            throw(DomainError())
        elseif vml_error == 2 || vml_error == 3 || vml_error == 4
            # Singularity, overflow, or underflow
            # I don't think Base throws on these
        elseif vml_error == 1000
            warn("VML does not support $(vml_get_accuracy); lower accuracy used instead")
        else
            error("an unexpected error occurred in VML ($vml_error)")
        end
    end
end

function vml_prefix(t::DataType)
    if t == Float32
        return "_vmls"
    elseif t == Float64
        return "_vmld"
    elseif t == Complex{Float32}
        return "_vmlc"
    elseif t == Complex{Float64}
        return "_vmlz"
    end
    error("unknown type $t")
end

function def_unary_op(tin, tout, jlname, jlname!, mklname)
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(tin))$mklname"))
    exports = Symbol[]
    isa(jlname, Expr) || push!(exports, jlname)
    isa(jlname!, Expr) || push!(exports, jlname!)
    @eval begin
        $(isempty(exports) ? nothing : Expr(:export, exports...))
        function $(jlname!){N}(out::Array{$tout,N}, A::Array{$tin,N})
            size(out) == size(A) || throw(DimensionMismatch())
            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
            vml_check_error()
            out
        end
        $(if tin == tout
            quote
                function $(jlname!)(A::Array{$tin})
                    ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, A)
                    vml_check_error()
                    A
                end
            end
        end)
        function $(jlname)(A::Array{$tin})
            out = similar(A, $tout)
            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
            vml_check_error()
            out
        end
    end
end

function def_binary_op(tin, tout, jlname, jlname!, mklname, broadcast)
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(tin))$mklname"))
    exports = Symbol[]
    isa(jlname, Expr) || push!(exports, jlname)
    isa(jlname!, Expr) || push!(exports, jlname!)
    @eval begin
        $(isempty(exports) ? nothing : Expr(:export, exports...))
        function $(jlname!){N}(out::Array{$tout,N}, A::Array{$tin,N}, B::Array{$tin,N})
            size(out) == size(A) == size(B) || $(broadcast ? :(return broadcast!($jlname, out, A, B)) : :(throw(DimensionMismatch())))
            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
            vml_check_error()
            out
        end
        function $(jlname){N}(A::Array{$tout,N}, B::Array{$tin,N})
            size(A) == size(B) || $(broadcast ? :(return broadcast($jlname, A, B)) : :(throw(DimensionMismatch())))
            out = similar(A)
            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
            vml_check_error()
            out
        end
    end
end

for t in (Float32, Float64, Complex64, Complex128)
    # Unary, real or complex
    def_unary_op(t, t, :(Base.acos), :acos!, :Acos)
    def_unary_op(t, t, :(Base.asin), :asin!, :Asin)
    def_unary_op(t, t, :(Base.acosh), :acosh!, :Acosh)
    def_unary_op(t, t, :(Base.asinh), :asinh!, :Asinh)
    def_unary_op(t, t, :(Base.sqrt), :sqrt!, :Sqrt)
    def_unary_op(t, t, :(Base.exp), :exp!, :Exp)
    def_unary_op(t, t, :(Base.log), :log!, :Ln)

    # Binary, real or complex
    def_binary_op(t, t, :(.^), :pow!, :Pow, true)
    def_binary_op(t, t, :(./), :divide!, :Div, true)
end

for t in (Float32, Float64)
    # Unary, real-only
    def_unary_op(t, t, :(Base.cbrt), :cbrt!, :Cbrt)
    def_unary_op(t, t, :(Base.expm1), :expm1!, :Expm1)
    def_unary_op(t, t, :(Base.log1p), :log1p, :Log1p)
    def_unary_op(t, t, :(Base.abs), :abs!, :Abs)
    def_unary_op(t, t, :(Base.abs2), :abs2!, :Sqr)
    def_unary_op(t, t, :(Base.ceil), :ceil!, :Ceil)
    def_unary_op(t, t, :(Base.floor), :floor!, :Floor)
    def_unary_op(t, t, :(Base.round), :round!, :Round)
    def_unary_op(t, t, :(Base.trunc), :trunc!, :Trunc)
    def_unary_op(t, t, :(Base.erf), :erf!, :Erf)
    def_unary_op(t, t, :(Base.erfc), :erfc!, :Erfc)
    def_unary_op(t, t, :(Base.erfinv), :erfinv!, :ErfInv)
    def_unary_op(t, t, :(Base.erfcinv), :erfcinv!, :ErfcInv)
    def_unary_op(t, t, :(Base.lgamma), :lgamma!, :LGamma)
    def_unary_op(t, t, :(Base.gamma), :gamma!, :TGamma)
    # Not in Base
    def_unary_op(t, t, :inv_cbrt, :inv_cbrt!, :InvCbrt)
    def_unary_op(t, t, :inv_sqrt, :inv_sqrt!, :InvSqrt)
    def_unary_op(t, t, :pow2o3, :pow2o3!, :Pow2o3)
    def_unary_op(t, t, :pow3o2, :pow3o2!, :Pow3o2)

    # Enabled only for Real. MKL guarantees higher accuracy, but at a
    # substantial performance cost.
    def_unary_op(t, t, :(Base.atan), :atan!, :Atan)
    def_unary_op(t, t, :(Base.cos), :cos!, :Cos)
    def_unary_op(t, t, :(Base.sin), :sin!, :Sin)
    def_unary_op(t, t, :(Base.tan), :tan!, :Tan)
    def_unary_op(t, t, :(Base.atanh), :atanh!, :Atanh)
    def_unary_op(t, t, :(Base.cosh), :cosh!, :Cosh)
    def_unary_op(t, t, :(Base.sinh), :sinh!, :Sinh)
    def_unary_op(t, t, :(Base.tanh), :tanh!, :Tanh)
    def_unary_op(t, t, :(Base.log10), :log10!, :Log10)

    # .^ to scalar power
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(t))Powx"))
    @eval begin
        export pow!
        function pow!{N}(out::Array{$t,N}, A::Array{$t,N}, b::$t)
            size(out) == size(A) || throw(DimensionMismatch())
            ccall(($mklfn, lib), Void, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
            vml_check_error()
            out
        end
        function (.^){N}(A::Array{$t,N}, b::$t)
            out = similar(A)
            ccall(($mklfn, lib), Void, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
            vml_check_error()
            out
        end
    end

    # Binary, real-only
    def_binary_op(t, t, :(Base.atan2), :atan2!, :Atan2, false)
    def_binary_op(t, t, :(Base.hypot), :hypot!, :Hypot, false)

    # Unary, complex-only
    def_unary_op(t, Complex{t}, :(Base.cis), :cis!, :CIS)
    # def_unary_op(Complex{t}, Complex{t}, :(Base.conj), :conj!, :Conj)
    def_unary_op(Complex{t}, t, :(Base.abs), :abs!, :Abs)
    def_unary_op(Complex{t}, t, :(Base.angle), :angle!, :Arg)

    # Binary, complex-only. These are more accurate but performance is
    # either equivalent to Base or slower.
    # def_binary_op(Complex{t}, Complex{t}, :(Base.(:+)), :add!, :Add, false)
    # def_binary_op(Complex{t}, Complex{t}, :(Base.(:.+)), :add!, :Add, true)
    # def_binary_op(Complex{t}, Complex{t}, :(Base.(:.*)), :multiply!, :Mul, true)
    # def_binary_op(Complex{t}, Complex{t}, :(Base.(:-)), :subtract!, :Sub, false)
    # def_binary_op(Complex{t}, Complex{t}, :(Base.(:.-)), :subtract!, :Sub, true)
    # def_binary_op(Complex{t}, Complex{t}, :multiply_conj, :multiply_conj!, :Mul, false)
end

export VML_LA, VML_HA, VML_EP, vml_set_accuracy, vml_get_accuracy
end


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

name = "faust"
version = v"2.32.16"

# Collection of sources required to build faust
sources = [
    GitSource("https://github.com/grame-cncm/faust.git",
              "53055b47f0ac81bee92b32db3c12729759e005f4"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/faust

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/disable_interp.patch
CMAKE_FLAGS=()
CMAKE_TARGET=${target}

if [[ "${target}" == *musl* ]]; then
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/execinfo.patch
    export CXXFLAGS="-DALPINE"
fi

if [[ "${target}" == *apple* ]]; then
    export MACOSX_DEPLOYMENT_TARGET=10.9
    export LDFLAGS="${LDFLAGS} -mmacosx-version-min=10.9"
    # If we're building for Apple, CMake gets confused with `aarch64-apple-darwin` and instead prefers
    # `arm64-apple-darwin`.  If this issue persists, we may have to change our triplet printing.
    if [[ "${target}" == aarch64* ]]; then
        CMAKE_TARGET=arm64-${target#*-}
    fi
fi

if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/ws2.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/windows-microhttpd.patch
    export LDFLAGS="${LDFLAGS} -lws2_32 -lmicrohttpd"
fi

if [[ "${target}" == *freebsd* ]]; then
    export LDFLAGS="${LDFLAGS} -lexecinfo"
fi

CMAKE_FLAGS+=(-DCMAKE_C_COMPILER_TARGET=${CMAKE_TARGET})
CMAKE_FLAGS+=(-DCMAKE_CXX_COMPILER_TARGET=${CMAKE_TARGET})

CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

CMAKE_FLAGS+=(-DUSE_LLVM_CONFIG=OFF)
CMAKE_FLAGS+=(-DLLVM_DIR=${prefix}/lib/cmake/llvm)

if [[ "${bb_full_target}" == x86_64-linux-musl-*-cxx03 ]]; then
    # For some reason this target requires "-lLLVM-11jl"
    # while others require "-lLLVM" to build.
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/set_llvm_libs_musl_cxx03.patch
else
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/set_llvm_libs.patch
fi

export CMAKEOPT="${CMAKE_FLAGS[@]}"
export PREFIX=$prefix

make -j$(nproc) all
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("faust", :faust),
    LibraryProduct("libfaust", :libfaust),
    LibraryProduct("libHTTPDFaust", :libHTTPDFaust),
    LibraryProduct("libOSCFaust", :libOSCFaust),

    # Generated list of Faust scripts using:
    # find $prefix/bin -maxdepth 1 -name 'faust2*' -printf '    FileProduct(joinpath("bin", "%f"), :%f),\n' | sort
    FileProduct(joinpath("bin", "faust2alqt"), :faust2alqt),
    FileProduct(joinpath("bin", "faust2alsaconsole"), :faust2alsaconsole),
    FileProduct(joinpath("bin", "faust2alsa"), :faust2alsa),
    FileProduct(joinpath("bin", "faust2android"), :faust2android),
    FileProduct(joinpath("bin", "faust2androidunity"), :faust2androidunity),
    FileProduct(joinpath("bin", "faust2api"), :faust2api),
    FileProduct(joinpath("bin", "faust2atomsnippets"), :faust2atomsnippets),
    FileProduct(joinpath("bin", "faust2au"), :faust2au),
    FileProduct(joinpath("bin", "faust2bela"), :faust2bela),
    FileProduct(joinpath("bin", "faust2caqt"), :faust2caqt),
    FileProduct(joinpath("bin", "faust2caqtios"), :faust2caqtios),
    FileProduct(joinpath("bin", "faust2csound"), :faust2csound),
    FileProduct(joinpath("bin", "faust2csvplot"), :faust2csvplot),
    FileProduct(joinpath("bin", "faust2dplug"), :faust2dplug),
    FileProduct(joinpath("bin", "faust2dssi"), :faust2dssi),
    FileProduct(joinpath("bin", "faust2dummy"), :faust2dummy),
    FileProduct(joinpath("bin", "faust2dummymem"), :faust2dummymem),
    FileProduct(joinpath("bin", "faust2eps"), :faust2eps),
    FileProduct(joinpath("bin", "faust2esp32"), :faust2esp32),
    FileProduct(joinpath("bin", "faust2faustvst"), :faust2faustvst),
    FileProduct(joinpath("bin", "faust2firefox"), :faust2firefox),
    FileProduct(joinpath("bin", "faust2gen"), :faust2gen),
    FileProduct(joinpath("bin", "faust2graph"), :faust2graph),
    FileProduct(joinpath("bin", "faust2graphviewer"), :faust2graphviewer),
    FileProduct(joinpath("bin", "faust2ios"), :faust2ios),
    FileProduct(joinpath("bin", "faust2jackconsole"), :faust2jackconsole),
    FileProduct(joinpath("bin", "faust2jack"), :faust2jack),
    FileProduct(joinpath("bin", "faust2jackrust"), :faust2jackrust),
    FileProduct(joinpath("bin", "faust2jackserver"), :faust2jackserver),
    FileProduct(joinpath("bin", "faust2jaqtchain"), :faust2jaqtchain),
    FileProduct(joinpath("bin", "faust2jaqt"), :faust2jaqt),
    FileProduct(joinpath("bin", "faust2juce"), :faust2juce),
    FileProduct(joinpath("bin", "faust2ladspa"), :faust2ladspa),
    FileProduct(joinpath("bin", "faust2linuxunity"), :faust2linuxunity),
    FileProduct(joinpath("bin", "faust2lv2"), :faust2lv2),
    FileProduct(joinpath("bin", "faust2mathdoc"), :faust2mathdoc),
    FileProduct(joinpath("bin", "faust2mathviewer"), :faust2mathviewer),
    FileProduct(joinpath("bin", "faust2max6"), :faust2max6),
    FileProduct(joinpath("bin", "faust2md"), :faust2md),
    FileProduct(joinpath("bin", "faust2msp"), :faust2msp),
    FileProduct(joinpath("bin", "faust2netjackconsole"), :faust2netjackconsole),
    FileProduct(joinpath("bin", "faust2netjackqt"), :faust2netjackqt),
    FileProduct(joinpath("bin", "faust2nodejs"), :faust2nodejs),
    FileProduct(joinpath("bin", "faust2octave"), :faust2octave),
    FileProduct(joinpath("bin", "faust2osxiosunity"), :faust2osxiosunity),
    FileProduct(joinpath("bin", "faust2owl"), :faust2owl),
    FileProduct(joinpath("bin", "faust2paqt"), :faust2paqt),
    FileProduct(joinpath("bin", "faust2pdf"), :faust2pdf),
    FileProduct(joinpath("bin", "faust2plot"), :faust2plot),
    FileProduct(joinpath("bin", "faust2png"), :faust2png),
    FileProduct(joinpath("bin", "faust2portaudiorust"), :faust2portaudiorust),
    FileProduct(joinpath("bin", "faust2puredata"), :faust2puredata),
    FileProduct(joinpath("bin", "faust2pure"), :faust2pure),
    FileProduct(joinpath("bin", "faust2raqt"), :faust2raqt),
    FileProduct(joinpath("bin", "faust2ros"), :faust2ros),
    FileProduct(joinpath("bin", "faust2rosgtk"), :faust2rosgtk),
    FileProduct(joinpath("bin", "faust2rpialsaconsole"), :faust2rpialsaconsole),
    FileProduct(joinpath("bin", "faust2rpinetjackconsole"), :faust2rpinetjackconsole),
    FileProduct(joinpath("bin", "faust2sam"), :faust2sam),
    FileProduct(joinpath("bin", "faust2sc"), :faust2sc),
    FileProduct(joinpath("bin", "faust2sig"), :faust2sig),
    FileProduct(joinpath("bin", "faust2sigviewer"), :faust2sigviewer),
    FileProduct(joinpath("bin", "faust2smartkeyb"), :faust2smartkeyb),
    FileProduct(joinpath("bin", "faust2sndfile"), :faust2sndfile),
    FileProduct(joinpath("bin", "faust2soul"), :faust2soul),
    FileProduct(joinpath("bin", "faust2supercollider"), :faust2supercollider),
    FileProduct(joinpath("bin", "faust2svg"), :faust2svg),
    FileProduct(joinpath("bin", "faust2svgviewer"), :faust2svgviewer),
    FileProduct(joinpath("bin", "faust2teensy"), :faust2teensy),
    FileProduct(joinpath("bin", "faust2unity"), :faust2unity),
    FileProduct(joinpath("bin", "faust2vcvrack"), :faust2vcvrack),
    FileProduct(joinpath("bin", "faust2vst"), :faust2vst),
    FileProduct(joinpath("bin", "faust2vsti"), :faust2vsti),
    FileProduct(joinpath("bin", "faust2w32max6"), :faust2w32max6),
    FileProduct(joinpath("bin", "faust2w32msp"), :faust2w32msp),
    FileProduct(joinpath("bin", "faust2w32puredata"), :faust2w32puredata),
    FileProduct(joinpath("bin", "faust2w32vst"), :faust2w32vst),
    FileProduct(joinpath("bin", "faust2w64max6"), :faust2w64max6),
    FileProduct(joinpath("bin", "faust2w64vst"), :faust2w64vst),
    FileProduct(joinpath("bin", "faust2wasm"), :faust2wasm),
    FileProduct(joinpath("bin", "faust2webaudiowasm"), :faust2webaudiowasm),
    FileProduct(joinpath("bin", "faust2webaudiowast"), :faust2webaudiowast),
    FileProduct(joinpath("bin", "faust2winunity"), :faust2winunity),

    FileProduct(joinpath("bin","faustoptflags"), :faustoptflags),
    FileProduct(joinpath("bin","usage.sh"), :usage_sh),

    FileProduct(joinpath("share", "faust", "stdfaust.lib"), :stdfaust_lib),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("LLVM_jll", v"11.0.1"),
    Dependency("libmicrohttpd_jll"),
    Dependency("libsndfile_jll"),
    BuildDependency("Ncurses_jll"),
    BuildDependency("XML2_jll"),
    BuildDependency("Zlib_jll"),
]

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

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using Yao, Yao.Blocks
using LinearAlgebra
using Test
using QuAlgorithmZoo

@testset "solving hamiltonian" begin
    nbit = 8
    h = heisenberg(nbit) |> cache
    @test ishermitian(h)
    reg = rand_state(nqubits(h))
    E0 = expect(h, reg)/nbit
    reg |> iter_groundstate!(h, niter=1000)
    EG = expect(h, reg)/nbit/4
    @test isapprox(EG, -0.4564, atol=1e-4)

    # using Time Evolution
    reg = rand_state(nqubits(h))
    reg |> itime_groundstate!(h, τ=20)
    EG = expect(h, reg)/nbit/4
    @test isapprox(EG, -0.4564, atol=1e-4)

    # using VQE
    N = 4
    h = heisenberg(N)
    E = eigen(h |> mat |> Matrix).values[1]
    c = random_diff_circuit(N, 5, [i=>mod(i,N)+1 for i=1:N], mode=:Merged) |> autodiff(:QC)
    dispatch!(c, :random)
    vqe_solve!(c, h)
    E2 = expect(h, zero_state(N) |> c)
    @test isapprox(E, E2, atol=1e-1)
end


text:
# ============================================================================ #
# Operations
# ============================================================================ #

function is_constant(x)
    return false
end
function is_constant(x::Constant)
    return true
end

function is_variable(x)
    return false
end
function is_variable(x::Variable)
    return true
end

function is_expression(x)
    return false
end
function is_expression(x::SymbolicExpression)
    return true
end

function is_composite_expression(x::T) where {T <: Union{Constant, Variable}}
    return false
end
function is_composite_expression(ex::SymbolicExpression)
    return true
end

function is_commutative(ex::SymbolicExpression)
    return  is_plus_operator(ex.head) ||
            is_product_operator(ex.head)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# matrix (and other) inversion

export invx

"""
`invx(A)` for a matrix `A` gives an exact matrix inverse.
"""
function invx(A::AbstractMatrix{T}) where {T}
    r, c = size(A)
    @assert r == c "Matrix must be square"

    if r == 0
        return A
    end

    B = eye(T, r)

    X = [A B]
    X = rrefx(X)

    d = [X[i, i] for i = 1:r]
    if all(d .== 1)
        return X[:, r+1:end]
    end
    @error "Matrix is not invertible"
end

function invx(x::T) where {T<:IntegerX}
    return 1 // x
end

invx(x) = 1 / x


text:
# orbit evolution on a 2D torus
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using DynamicalSystems, InteractiveDynamics
import GLMakie
using OrdinaryDiffEq

R = 2.0
r = 1.0

function torus(u)
    θ, φ = u
    x = (R + r*cos(θ))*cos(φ)
    y = (R + r*cos(θ))*sin(φ)
    z = r*sin(θ)
    return SVector(x, y, z)
end

function quasiperiodic_f(u, p, t)
    # here we make the frequency ratio a state variable, because the
    # interactive_evolution application doesn't allow different parameters for different
    # initial conditions
    ω = u[3]
    θdot = ω
    φdot = 1.0
    return SVector(θdot, φdot, 0.0)
end

ds = ContinuousDynamicalSystem(quasiperiodic_f, [0.0, 0.0, 0.0], nothing)
frequencies² = [7, 9]
u0s = [[0, 0, √x] for x ∈ frequencies²]
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.01)

lims = ((-R-r, R+r), (-R-r, R+r), (-3r, 3r))

plotkwargs = [(linewidth = i^2*2.0,) for i in 1:length(u0s)]

fig, obs = interactive_evolution(
    ds, u0s; tail = 20000, diffeq, colors = COLORSCHEME[1:length(u0s)],
    transform = torus, lims, plotkwargs
)

main = GLMakie.content(fig[1,1])

freqs = collect("√"*string(s) for s in frequencies²)
main.title = "Frequency ratios: "*join(freqs, ", ")

# Create a wireframe for the torus
# angles = [SVector(θ, φ) for θ ∈ range(0, 2π; length = 100) for φ ∈ range(0, 2π; length = 100)]
# coords = torus.(angles)
# wireframe!(main, coords)

U = range(-π, π; length = 200)
V = range(-π, π; length = 50)

x1 = [(R + r*cos(θ))*cos(φ)      for θ in U, φ in V]
y1 = [(R + r*cos(θ))*sin(φ)      for θ in U, φ in V]
z1 = [r*sin(θ)                   for θ in U, φ in V]
GLMakie.wireframe!(main, x1,y1,z1;
    shading = false, color = BLACK, linewidth = 0.1
)

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 10)

# %% Clean figure for torus in book
GLMakie.hidedecorations!(main)
GLMakie.hidespines!(main)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# File loaded into eFEMpart

### Abstract Types ###

abstract type AbstractParameter end

abstract type AbstractDummyParameter <: AbstractParameter end
abstract type AbstractConstantParameter <: AbstractParameter end
abstract type AbstractVariableParameter <: AbstractParameter end

### Concrete Types Definitions ###

mutable struct PoissonParam <: AbstractDummyParameter
  κ::Float64
end

mutable struct HeatConstParam <: AbstractConstantParameter
  κ::Float64
end
mutable struct HeatVarParam <: AbstractVariableParameter
  κ::Vector{Float64}
end

mutable struct DarcyConstParam <: AbstractConstantParameter
  α::Float64
end
mutable struct DarcyVarParam <: AbstractVariableParameter
  α::Vector{Float64}
end

mutable struct AdvDiffConstParam <: AbstractConstantParameter
  u::Vector{Float64}
  v::Vector{Float64}
  Pe::Float64
end
mutable struct AdvDiffVarParam <: AbstractVariableParameter
  u::Vector{Float64}
  v::Vector{Float64}
  Pe::Vector{Float64}
end

mutable struct FluidConstParam <: AbstractConstantParameter
  μ::Float64
end
mutable struct FluidVarParam <: AbstractVariableParameter
  μ::Vector{Float64}
end

mutable struct BrinkmanConstParam <: AbstractConstantParameter
  μ::Float64
  α::Float64
end
mutable struct BrinkmanVarParam <: AbstractVariableParameter
  μ::Vector{Float64}
  α::Vector{Float64}
end

mutable struct BrinkmanMPConstParam <: AbstractConstantParameter
  α1::Float64
  α2::Float64
  α3::Float64
end
mutable struct BrinkmanMPVarParam <: AbstractVariableParameter
  α1::Vector{Float64}
  α2::Vector{Float64}
  α3::Vector{Float64}
end

### Abstract Problem-specific Unions ###

AbstractHeat       = Union{HeatConstParam,HeatVarParam}
AbstractAdvDiff    = Union{AdvDiffConstParam,AdvDiffVarParam}
AbstractDarcy      = Union{DarcyConstParam,DarcyVarParam}
AbstractFluid      = Union{FluidConstParam,FluidVarParam}
AbstractBrinkman   = Union{BrinkmanConstParam,BrinkmanVarParam}
AbstractBrinkmanMP = Union{BrinkmanMPConstParam,BrinkmanMPVarParam}

### Constructor Definitions ###

function HeatParam(κ::T) where T<:Real
  return HeatConstParam(κ)
end
function HeatParam(κ::Vector{T}) where T<:Real
  return HeatVarParam(κ)
end

function AdvDiffParam(u,v,κ::T) where T<:Real
  return AdvDiffConstParam(u,v,κ)
end
function AdvDiffParam(u,v,κ::Vector{T}) where T<:Real
  return AdvDiffVarParam(u,v,κ)
end

function DarcyParam(α::T) where T<:Real
  return DarcyConstParam(α)
end
function DarcyParam(α::Vector{T}) where T<:Real
  return DarcyVarParam(α)
end

function FluidParam(μ::T) where T<:Real
  return FluidConstParam(μ)
end
function FluidParam(μ::Vector{T}) where T<:Real
  return FluidVarParam(μ)
end

function BrinkmanParam(μ::T,α::T) where T<:Real
  return BrinkmanConstParam(μ,α)
end
function BrinkmanParam(μ::T,α::Vector{T}) where T<:Real
  return BrinkmanVarParam(μ*zeros(T,length(α)),α)
end
function BrinkmanParam(μ::Vector{T},α::T) where T<:Real
  return BrinkmanVarParam(μ,α*ones(T,length(μ)))
end
function BrinkmanParam(μ::Vector{T},α::Vector{T}) where T<:Real
  return BrinkmanVarParam(μ,α)
end

function BrinkmanMPParam(α1::T,α2::T,α3::T) where T<:Real
  return BrinkmanMPConstParam(α1,α2,α3)
end
function BrinkmanMPParam(α1::Vector{T},α2::T,α3::T) where T<:Real
  o = ones(T,length(α1))
  return BrinkmanMPVarParam(α1,α2*o,α3*o)
end
function BrinkmanMPParam(α1::T,α2::Vector{T},α3::T) where T<:Real
  o = ones(T,length(α2))
  return BrinkmanMPVarParam(α1*o,α2,α3*o)
end
function BrinkmanMPParam(α1::T,α2::T,α3::Vector{T}) where T<:Real
  o = ones(T,length(α1))
  return BrinkmanMPVarParam(α1,α2*o,α3*o)
end
function BrinkmanMPParam(α1::Vector{T},α2::Vector{T},α3::T) where T<:Real
  o = ones(T,length(α1))
  return BrinkmanMPVarParam(α1,α2,α3*o)
end
function BrinkmanMPParam(α1::Vector{T},α2::T,α3::Vector{T}) where T<:Real
  o = ones(T,length(α1))
  return BrinkmanMPVarParam(α1,α2*o,α3)
end
function BrinkmanMPParam(α1::T,α2::Vector{T},α3::Vector{T}) where T<:Real
  o = ones(T,length(α2))
  return BrinkmanMPVarParam(α1*o,α2,α3)
end
function BrinkmanMPParam(α1::Vector{T},α2::Vector{T},α3::Vector{T}) where T<:Real
  return BrinkmanMPVarParam(α1,α2,α3)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    bat_findmedian(
        samples::DensitySampleVector
    )

The function computes the median of marginalized `samples`.

Returns a NamedTuple of the shape

```julia
(result = v, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_findmedian`, add methods to
    `bat_findmedian_impl` instead.
"""
function bat_findmedian end
export bat_findmedian

function bat_findmedian_impl end

function bat_findmedian(samples::DensitySampleVector)
    r = bat_findmedian_impl(samples::DensitySampleVector)
    result_with_args(r)
end


text:
get_params_kwargs(::GPE; kwargs...) = delete!(Dict(kwargs), :lik)
get_params_kwargs(::GPMC; kwargs...) = delete!(Dict(kwargs), :noise)

"""
    optimize!(gp::GPBase; kwargs...)

Optimise the hyperparameters of Gaussian process `gp` based on type II maximum likelihood estimation. This function performs gradient based optimisation using the Optim pacakge to which the user is referred to for further details.

# Keyword arguments:
    * `domean::Bool`: Mean function hyperparameters should be optmized
    * `kern::Bool`: Kernel function hyperparameters should be optmized
    * `noise::Bool`: Observation noise hyperparameter should be optimized (GPE only)
    * `lik::Bool`: Likelihood hyperparameters should be optimized (GPMC only)
    * `kwargs`: Keyword arguments for the optimize function from the Optim package
"""
function optimize!(gp::GPBase; method = LBFGS(), domean::Bool = true, kern::Bool = true,
                   noise::Bool = true, lik::Bool = true, kwargs...)
    params_kwargs = get_params_kwargs(gp; domean=domean, kern=kern, noise=noise, lik=lik)
    # println(params_kwargs)
    func = get_optim_target(gp; params_kwargs...)
    init = get_params(gp; params_kwargs...)  # Initial hyperparameter values
    results = optimize(func, init; method=method, kwargs...)     # Run optimizer
    set_params!(gp, Optim.minimizer(results); params_kwargs...)
    update_target!(gp)
    return results
end

function get_optim_target(gp::GPBase; params_kwargs...)
    function ltarget(hyp::VecF64)
        prev = get_params(gp; params_kwargs...)
        try
            set_params!(gp, hyp; params_kwargs...)
            update_target!(gp)
            return -gp.target
        catch err
            # reset parameters to remove any NaNs
            set_params!(gp, prev; params_kwargs...)

            if !all(isfinite.(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, LinearAlgebra.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end

    function ltarget_and_dltarget!(grad::VecF64, hyp::VecF64)
        prev = get_params(gp; params_kwargs...)
        try
            set_params!(gp, hyp; params_kwargs...)
            update_target_and_dtarget!(gp; params_kwargs...)
            grad[:] = -gp.dtarget
            return -gp.target
        catch err
            # reset parameters to remove any NaNs
            set_params!(gp, prev; params_kwargs...)
            if !all(isfinite.(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, LinearAlgebra.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end

    function dltarget!(grad::VecF64, hyp::VecF64)
        ltarget_and_dltarget!(grad::VecF64, hyp::VecF64)
    end

    xinit = get_params(gp; params_kwargs...)
    func = OnceDifferentiable(ltarget, dltarget!, ltarget_and_dltarget!, xinit)
    return func
end


text:
@testset "tricky" begin

    x = [2.0, 4.0, 6.0]
    myfirst(x) = x[1]

    @test xgrad(myfirst; x=x) == (2.0, [1.0, 0.0, 0.0])

end


text:
export border, borderstyle, bordercolor, borderwidth, hline, vline,
       fontweight, background

# TODO: add more background features
# maybe a way to load images?
function background(spec, elem=nothing)
    style(elem, "background" => spec)
end

function _border_prefix(side::Union{Symbol, String})
    assertoneof(side, ["top", "bottom", "left", "right"], "side")
    "border$(uppercasefirst(side))"
end

function assertstyle(style)
    assertoneof(style, ["solid", "dotted", "dashed", "none"], "border style")
end

function border(side::Union{Symbol, String}, styl, width, color, elem=nothing)
    assertstyle(styl)
    style(elem, _border_prefix(side) => "$styl $(JSON.lower(width)) $(JSON.lower(parse(Colorant, color)))")
end

function border(sides::AbstractArray, args...)
    x = border(sides[1], args...)
    for i = 2:length(sides)
        x = border(sides[i], args...)
    end
    x
end

function borderstyle(side::Union{Symbol, String}, styl, elem=nothing)
    assertstyle(styl)
    style(elem, _border_prefix(side) * "Style" => styl)
end

function bordercolor(side::Union{Symbol, String}, color, elem=nothing)
    style(elem, _border_prefix(side) * "Color" => color)
end

function borderwidth(side::Union{Symbol, String}, width, elem=nothing)
    style(elem, _border_prefix(side) * "Width" => width)
end

"""
    hline()

    Draw a horizonal line.

Keyword Arguments:
=================
 - style: Indicates whether line should be solid or dotted.
Defaults to "solid"
- w: Specifies line width. Defaults to 1px
- color: specifies color in hex code RGB. Defaults to "#dedede".
"""
function hline(; style="solid", w=1px, color="#dedede")
    border("bottom", style, w, color, alignself("stretch", empty()))
end

"""
    vline()

    Draw a vertical line.

Keyword Arguments:
=================
 - style: Indicates whether line should be solid or dotted.
Defaults to "solid"
- w: Specifies line width. Defaults to 1px
- color: specifies color in hex code RGB. Defaults to "#dedede".
"""
function vline(; style="solid", w=1px, color="#dedede")
    border("left", style, w, color, alignself("stretch", empty()))
end

function fontweight(weight, elem=nothing)
    style(elem, "fontWeight" => weight)
end


text:
using Documenter, LowFrequencyEconometrics

makedocs(modules=LowFrequencyEconometrics,
        doctest=true)

deploydocs(deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo = "https://github.com/p-chaim/LowFrequencyEconometrics.jl.git",
    julia  = "0.5.2",
    osname = "windows")


text:
# computational utilities

"""
    θ(x)

Heaviside step function. Returns one if the argument is positive, zero otherwise.
"""
function θ(x::T) where T
    if x > zero(T)
        return one(T)
    end
    return zero(T)
end

"""
    xθ(x)

Heaviside step times x. Returns x if the argument is positive, zero otherwise.
"""
function xθ(x::T) where T
    if x > zero(T)
        return x
    end
    return zero(T)
end

"""
    wavevectors3D([T=Float64], dims, box_size=(2π, 2π, 2π))

# Arguments:
- `T` : output element type
- `dims`: array or tuple containing size
- `box_size`: side lengths of box, normalizes the wavenumber

# Returns:
- `Tuple{Vector{T},Vector{T},Vector{T}}`: kx, ky, kz wavevectors

# Examples
```julia-repl
julia> CosmoMMF.wavevectors3D((6,2,2))
([0.0, 1.0, 2.0, -3.0, -2.0, -1.0], [0.0, -1.0], [0.0, -1.0])
```
"""
function wavevectors3D(T::Type{<:Real}, dims, box_size=(2π, 2π, 2π))
    sample_rate = T.(2π .* dims ./ box_size)
    kx = fftfreq(dims[1], sample_rate[1]) .* (2π / dims[1])
    ky = fftfreq(dims[2], sample_rate[2]) .* (2π / dims[2])
    kz = fftfreq(dims[3], sample_rate[3]) .* (2π / dims[3])
    return kx, ky, kz
end

# default T is Float64
wavevectors3D(dims, box_size=(2π, 2π, 2π)) = wavevectors3D(Float64, dims, box_size)

function test_print()
    print("This test function works! (totally)")
end

function wall(n::Int64)
    """
    Builds a wall out of a nxnxn array. All values in the wall
    are = 1. Other points = 0.
    """
    array = zeros((n,n,n))
    array[:,:,:] .= 0.1
    index = n ÷ 2
    array[index,:,:] .= 1
    return array
end

function cylinder(n::Int64,r::Float64)
    """
    Returns a nxnxn array with a cylinder in the center with
    radius r. All points within the radius are = 1, the rest
    are 0.
    """
    array  = zeros(Float64,(n,n,n))
    array[:,:,:] .= 0.1
    for x::Int64 in 1:n
        for y::Int64 in 1:n
            for z::Int64 in 1:n
                if ((x - (n)/2)^2 + (y - (n)/2)^2) < r^2
                    array[x,y,z] = 1
                end
            end
        end
    end
    return array
end

function sphere(n::Int64,r::Float64)
    """
    Returns a nxnxn array with a sphere in the center with
    radius r. All points within the radius are = 1, the rest
    are 0.
    """
    array = zeros(Float64,(n,n,n))
    array[:,:,:] .= 0.1
    for x::Int64 in 1:n
        for y::Int64 in 1:n
            for z::Int64 in 1:n
                if ((x - (n)/2)^2 + (y - (n)/2)^2 + (z - (n)/2)^2) < r^2
                    array[x,y,z] = 1
                end
            end
        end
    end
    return array
end

function save_max_sigs(output_directory, id_string, sig_array::AbstractArray{T,4}) where T
    """
    A simple function to save the signatures to an output directory in the form of a jld2 file.
    """
    @save output_directory * id_string * ".jld2" sig_array
end

function save_hessian_component(output_directory, Rs, id_string, hessian_component::AbstractArray{T,3}) where T
    """
    A helper function to save a hessian component.

    Args:
        output_directory [string]: where should the hessian component be saved to?
        id_string [string]: component number in matrix
                            -> [[1 2 3],
                                [  4 5],
                                [    6]]
        subvolume_divisions [Int]: number of subvolumes to divide original hessian into on each side. If = 2, then we make 2^3 subvolumes. (Default = 2)
        hessian_component [AbstractArray{T,3}]: one of the components
    """
    dims = size(hessian_component)

    div_x, div_y, div_z = dims .÷ 2

    #chop hessian into 8ths
    hessian_component_000 = hessian_component[1:div_x, 1:div_y, 1:div_z]
    hessian_component_001 = hessian_component[1:div_x, 1:div_y, div_z+1:dims[3]]
    hessian_component_010 = hessian_component[1:div_x, div_y+1:dims[2], 1:div_z]
    hessian_component_100 = hessian_component[div_x+1:dims[1], 1:div_y, 1:div_z]
    hessian_component_011 = hessian_component[1:div_x, div_y+1:dims[2], div_z+1:dims[3]]
    hessian_component_110 = hessian_component[div_x+1:dims[1], div_y+1:dims[2], 1:div_z]
    hessian_component_101 = hessian_component[div_x+1:dims[1], 1:div_y, div_z+1:dims[3]]
    hessian_component_111 = hessian_component[div_x+1:dims[1], div_y+1:dims[2], div_z+1:dims[3]]

    @save output_directory * "hessian_component_000" * "_" * string(Rs) * id_string * ".jld2" hessian_component_000
    @save output_directory * "hessian_component_001" * "_" * string(Rs) * id_string * ".jld2" hessian_component_001
    @save output_directory * "hessian_component_010" * "_" * string(Rs) * id_string * ".jld2" hessian_component_010
    @save output_directory * "hessian_component_100" * "_" * string(Rs) * id_string * ".jld2" hessian_component_100
    @save output_directory * "hessian_component_011" * "_" * string(Rs) * id_string * ".jld2" hessian_component_011
    @save output_directory * "hessian_component_110" * "_" * string(Rs) * id_string * ".jld2" hessian_component_110
    @save output_directory * "hessian_component_101" * "_" * string(Rs) * id_string * ".jld2" hessian_component_101
    @save output_directory * "hessian_component_111" * "_" * string(Rs) * id_string * ".jld2" hessian_component_111

end

function compile_sub_hessians(f_Rn::AbstractArray{T,3}, Rs, subvolume, input_directory) where T

    id_strings = ["_1", "_2", "_3", "_4", "_5", "_6"]
    dims = size(f_Rn) .÷ 2
    hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_000" * id_string * ".jld2" hessian_component_000
#         hessian_sub[:,:,:,i] .= hessian_component_000
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_000" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_001" * id_string * ".jld2" hessian_component_001
#         hessian_sub[:,:,:,i] .= hessian_component_001
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_001" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_010" * id_string * ".jld2" hessian_component_010
#         hessian_sub[:,:,:,i] .= hessian_component_010
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_010" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_100" * id_string * ".jld2" hessian_component_100
#         hessian_sub[:,:,:,i] .= hessian_component_100
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_100" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_011" * id_string * ".jld2" hessian_component_011
#         hessian_sub[:,:,:,i] .= hessian_component_011
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_011" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_110" * id_string * ".jld2" hessian_component_110
#         hessian_sub[:,:,:,i] .= hessian_component_110
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_110" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

#     for (i, id_string) in enumerate(id_strings)
#         @load input_directory * "hessian_component_101" * id_string * ".jld2" hessian_component_101
#         hessian_sub[:,:,:,i] .= hessian_component_101
#         #figure out how to delete files?
#     end

#     @save input_directory * "hessian_component_101" * "_full" * ".jld2" real(hessian_sub)
#     hessian_sub = zeros(Complex{T}, (dims[1], dims[2], dims[3], 6))

    for (i, id_string) in enumerate(id_strings)
        if subvolume == "000"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_000
        hessian_sub[:,:,:,i] .= hessian_component_000

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "001"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_001
        hessian_sub[:,:,:,i] .= hessian_component_001

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "010"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_010
        hessian_sub[:,:,:,i] .= hessian_component_010

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "100"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_100
        hessian_sub[:,:,:,i] .= hessian_component_100

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "011"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_011
        hessian_sub[:,:,:,i] .= hessian_component_011

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "110"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_110
        hessian_sub[:,:,:,i] .= hessian_component_110

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "101"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_101
        hessian_sub[:,:,:,i] .= hessian_component_101

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        elseif subvolume == "111"
            @load input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2" hessian_component_111
        hessian_sub[:,:,:,i] .= hessian_component_111

        #delete old slice we don't need anymore
        file = input_directory * "hessian_component_" * subvolume * "_" * string(Rs) * id_string * ".jld2"
        run(`rm $file`)

        else
            print("Give a correct subvolume string. Example: '000' ")
        end

    end

    hessian_sub = real(hessian_sub)
    @save input_directory * "hessian_component_" * subvolume * "_full" * ".jld2" hessian_sub

end


text:
using DelimitedFiles

input = joinpath(@__DIR__, "input")
map = readdlm(input, String)

function solve(map, slope, bounds)
    hits = 0
    pos = CartesianIndex(1, 1)
    pos += slope

    while pos[1] <= bounds[1]
        x = pos[2]
        y = pos[1]
        hits += (map[y][x] == '#' ? 1 : 0)
        pos += slope
        pos[2] > bounds[2] && (pos -= CartesianIndex(0, bounds[2]))
    end
    return hits
end

const bounds = CartesianIndex(length(map), length(map[1]))

p1 = solve(map, CartesianIndex(1, 3), bounds)

const slopes = (CartesianIndex(1, 1), CartesianIndex(1, 3), CartesianIndex(1, 5), CartesianIndex(1, 7), CartesianIndex(2, 1))

p2 = prod(solve(map, slope, bounds) for slope in slopes)

@assert(p1 == 299)
@assert(p2 == 3621285278)

println("-----------------------------------------------------------------------")
println("toboggan trajectory -- part one :: $p1")
println("toboggan trajectory -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
using NURBS
using Base.Test

@testset "RB-Spline curve generator" begin

	@testset "rbasis" begin
		x = knot(5, 4)
		h = [1, 1, 0.25, 1, 1]
		@test typeof(rbasis(5, 4, 0.5, x, h)) == Array{Float64,1}
		@test isapprox(rbasis(5, 4, 0.5, x, h), [0.153846, 0.730769, 0.0769231, 0.0384615, 0.0], atol = 1e-5)
	end

	@testset "rbspline" begin
		b = [0 1 5/2 4 5; 1 2 0 2 0; 0 0 0 0 0]
		h = [1, 1, 0.25, 1, 1]
		@test typeof(rbspline(5, 3, 10, b, h)[1]) == Array{Float64,2}
		@test typeof(rbspline(5, 3, 10, b, h)[2]) == Array{Array{Int64,1},1}
		@test isapprox(rbspline(5, 3, 10, b, h)[1], [0.0 0.557971 0.966667 1.3 1.95455 3.04545 3.7 4.03333 4.44203 5.0; 1.0 1.50725 1.73333 1.6 1.21212 1.21212 1.6 1.6 1.04348 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], atol = 1e-5)
	end

	@testset "rbsplinu" begin
		b = [0 1 5/2 4 5; 1 2 0 2 0; 0 0 0 0 0]
		h = [1, 1, 0.25, 1, 1]
		@test typeof(rbsplinu(5, 3, 10, b, h)[1]) == Array{Float64,2}
		@test typeof(rbsplinu(5, 3, 10, b, h)[2]) == Array{Array{Int64,1},1}
		@test isapprox(rbsplinu(5, 3, 10, b, h)[1], [0.5 0.789855 1.03333 1.3 1.95455 3.04545 3.7 3.96667 4.21014 4.75; 1.5 1.73913 1.8 1.6 1.21212 1.21212 1.6 1.73333 1.50725 0.5; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], atol = 1e-4)
	end

end


text:
# ==============================================================================
#    Generic Callbacks in Gurobi
# ==============================================================================

"""
    CallbackFunction()

Set a generic Gurobi callback function.

Note: before accessing `MOI.CallbackVariablePrimal`, you must call either
`Gurobi.cbget_mipsol_sol(model, cb_data, cb_where)` or
`Gurobi.cbget_mipsol_rel(model, cb_data, cb_where)`.
"""
struct CallbackFunction <: MOI.AbstractCallback end

function MOI.set(model::Optimizer, ::CallbackFunction, f::Function)
    model.has_generic_callback = true
    set_callback_func!(model.inner, (cb_data, cb_where) -> begin
        model.callback_state = CB_GENERIC
        f(cb_data, cb_where)
        model.callback_state = CB_NONE
    end)
    update_model!(model.inner)
    return
end
MOI.supports(::Optimizer, ::CallbackFunction) = true

"""
    cbget_mipsol_sol(model::Optimizer, cb_data, cb_where)

Load the solution at a CB_MIPSOL node so that it can be accessed using
`MOI.CallbackVariablePrimal`.
"""
function cbget_mipsol_sol(model::Optimizer, cb_data, cb_where)
    resize!(model.callback_variable_primal, length(model.variable_info))
    cbget_mipsol_sol(cb_data, cb_where, model.callback_variable_primal)
    return
end

"""
    cbget_mipsol_rel(model::Optimizer, cb_data, cb_where)

Load the solution at a CB_MIPNODE node so that it can be accessed using
`MOI.CallbackVariablePrimal`.
"""
function cbget_mipsol_rel(model::Optimizer, cb_data, cb_where)
    resize!(model.callback_variable_primal, length(model.variable_info))
    cbget_mipnode_rel(cb_data, cb_where, model.callback_variable_primal)
    return
end

# ==============================================================================
#    MOI callbacks
# ==============================================================================

function default_moi_callback(model::Optimizer)
    return (cb_data, cb_where) -> begin
        if cb_where == CB_MIPSOL
            cbget_mipsol_sol(model, cb_data, cb_where)
            if model.lazy_callback !== nothing
                model.callback_state = CB_LAZY
                model.lazy_callback(cb_data)
            end
        elseif cb_where == CB_MIPNODE
            if cbget_mipnode_status(cb_data, cb_where) != 2
                return  # Solution is something other than optimal.
            end
            cbget_mipsol_rel(model, cb_data, cb_where)
            if model.lazy_callback !== nothing
                model.callback_state = CB_LAZY
                model.lazy_callback(cb_data)
            end
            if model.user_cut_callback !== nothing
                model.callback_state = CB_USER_CUT
                model.user_cut_callback(cb_data)
            end
            if model.heuristic_callback !== nothing
                model.callback_state = CB_HEURISTIC
                model.heuristic_callback(cb_data)
            end
        end
        model.callback_state = CB_NONE
    end
end

function MOI.get(
    model::Optimizer,
    ::MOI.CallbackVariablePrimal{CallbackData},
    x::MOI.VariableIndex
)
    return model.callback_variable_primal[_info(model, x).column]
end

# ==============================================================================
#    MOI.LazyConstraint
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.LazyConstraintCallback, cb::Function)
    MOI.set(model, MOI.RawParameter("LazyConstraints"), 1)
    model.lazy_callback = cb
    return
end
MOI.supports(::Optimizer, ::MOI.LazyConstraintCallback) = true

function MOI.submit(
    model::Optimizer,
    cb::MOI.LazyConstraint{CallbackData},
    f::MOI.ScalarAffineFunction{Float64},
    s::Union{MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, MOI.EqualTo{Float64}}
)
    if model.callback_state == CB_USER_CUT
        throw(MOI.InvalidCallbackUsage(MOI.UserCutCallback(), cb))
    elseif model.callback_state == CB_HEURISTIC
        throw(MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), cb))
    elseif !iszero(f.constant)
        throw(MOI.ScalarFunctionConstantNotZero{Float64, typeof(f), typeof(s)}(f.constant))
    end
    indices, coefficients = _indices_and_coefficients(model, f)
    sense, rhs = _sense_and_rhs(s)
    cblazy(cb.callback_data, Cint.(indices), coefficients, Char(sense), rhs)
    return
end
MOI.supports(::Optimizer, ::MOI.LazyConstraint{CallbackData}) = true

# ==============================================================================
#    MOI.UserCutCallback
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.UserCutCallback, cb::Function)
    model.user_cut_callback = cb
    return
end
MOI.supports(::Optimizer, ::MOI.UserCutCallback) = true

function MOI.submit(
    model::Optimizer,
    cb::MOI.UserCut{CallbackData},
    f::MOI.ScalarAffineFunction{Float64},
    s::Union{MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, MOI.EqualTo{Float64}}
)
    if model.callback_state == CB_LAZY
        throw(MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), cb))
    elseif model.callback_state == CB_HEURISTIC
        throw(MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), cb))
    elseif !iszero(f.constant)
        throw(MOI.ScalarFunctionConstantNotZero{Float64, typeof(f), typeof(s)}(f.constant))
    end
    indices, coefficients = _indices_and_coefficients(model, f)
    sense, rhs = _sense_and_rhs(s)
    cbcut(cb.callback_data, Cint.(indices), coefficients, Char(sense), rhs)
    return
end
MOI.supports(::Optimizer, ::MOI.UserCut{CallbackData}) = true

# ==============================================================================
#    MOI.HeuristicCallback
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.HeuristicCallback, cb::Function)
    model.heuristic_callback = cb
    return
end
MOI.supports(::Optimizer, ::MOI.HeuristicCallback) = true

function MOI.submit(
    model::Optimizer,
    cb::MOI.HeuristicSolution{CallbackData},
    variables::Vector{MOI.VariableIndex},
    values::MOI.Vector{Float64}
)
    if model.callback_state == CB_LAZY
        throw(MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), cb))
    elseif model.callback_state == CB_USER_CUT
        throw(MOI.InvalidCallbackUsage(MOI.UserCutCallback(), cb))
    end
    solution = fill(GRB_UNDEFINED, MOI.get(model, MOI.NumberOfVariables()))
    for (var, value) in zip(variables, values)
        solution[_info(model, var).column] = value
    end
    obj = cbsolution(cb.callback_data, solution)
    return obj < GRB_INFINITY ? MOI.HEURISTIC_SOLUTION_ACCEPTED : MOI.HEURISTIC_SOLUTION_REJECTED
end
MOI.supports(::Optimizer, ::MOI.HeuristicSolution{CallbackData}) = true


text:
def epsilon 0.000001

def int_pow(x, n) if 1 == n x else x*int_pow(x,n-1)

def pow(x,n) {
    {result, value, power} = (1f, x as Float, n)
    while power > epsilon {
        if power % 2 == 1
            result = result*value
        {value, power} = (value*value, power/2)
    }
    result as Int
}

def sqr(x) x*x

def sqrt(n) {
    {x, y} = (n as Float, 1f)
    while (x - y) > epsilon {
        x = (x + y) / 2f
        y = n / x
    }
    x
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
function leading_bits(s::String)
    x = UInt(0)
    for i = 1:min(sizeof(x), sizeof(s))
        @inbounds x |= UInt(codeunit(s, i)) << ((sizeof(x) - i) * 8)
    end
    return x
end

function leading_bits_with_fast_path(s::String)
    if sizeof(s) == 8
        return ntoh(unsafe_load(Ptr{UInt64}(pointer(s))))
    end

    x = UInt(0)
    for i = 1:min(sizeof(x), sizeof(s))
        @inbounds x |= UInt(codeunit(s, i)) << ((sizeof(x) - i) * 8)
    end
    return x
end

function load_bits_with_padding(s::String, skipbytes = 0)
    n = sizeof(s)
    remaining_bytes_to_load = min(sizeof(UInt), n)
    # start
    res = zero(UInt)
    shift_for_padding = 64

    if remaining_bytes_to_load == 8
        res = ntoh(unsafe_load(Ptr{UInt64}(pointer(s, skipbytes+1))))
    else
        if  remaining_bytes_to_load >= 4
            res |= Base.zext_int(UInt, ntoh(unsafe_load(Ptr{UInt32}(pointer(s, skipbytes+1))))) << (shift_for_padding - 32)
            skipbytes += 4
            remaining_bytes_to_load -= 4
            shift_for_padding -= 32
        end
        if  remaining_bytes_to_load >= 2
            res |= Base.zext_int(UInt, ntoh(unsafe_load(Ptr{UInt16}(pointer(s, skipbytes+1))))) << (shift_for_padding - 16)
            skipbytes += 2
            remaining_bytes_to_load -= 2
            shift_for_padding -= 16
        end
        if  remaining_bytes_to_load >= 1
            res |= Base.zext_int(UInt, ntoh(unsafe_load(Ptr{UInt8}(pointer(s, skipbytes+1))))) << (shift_for_padding - 8)
            skipbytes += 1
            remaining_bytes_to_load -= 1
            shift_for_padding -= 8
        end
    end

    return res
end

mutable struct MutableUInt val::UInt end

import Base: unsafe_copy!, unsafe_copyto!
if !isdefined(Symbol("unsafe_copyto!"))
    unsafe_copyto! = unsafe_copy!
end

function copy_bits_fast_path(s::String, skipbytes = 0)
    n = sizeof(s)
    remaining_bytes_to_load = min(sizeof(UInt), n)
    # start

    if remaining_bytes_to_load == 8
        return ntoh(unsafe_load(Ptr{UInt64}(pointer(s))))
    else
        res = MutableUInt(0)
        p=reinterpret(Ptr{UInt8}, pointer_from_objref(res))
        if  remaining_bytes_to_load >= 4
            unsafe_copyto!(p, pointer(s), 4)
            skipbytes += 4
            remaining_bytes_to_load -= 4
        end
        if  remaining_bytes_to_load >= 2
            unsafe_copyto!(p+skipbytes, pointer(s, skipbytes+1),2)
            skipbytes += 2
            remaining_bytes_to_load -= 2
        end
        if  remaining_bytes_to_load >= 1
            unsafe_copyto!(p+skipbytes, pointer(s, skipbytes+1),1)
        end
    end

    ntoh(res.val)
end

# all char specialisation
primitive type Bits24 24 end
primitive type Bits40 40 end
primitive type Bits48 48 end
primitive type Bits56 56 end

function all_lens(s::String, skipbytes = 0)::UInt
    n = sizeof(s)
    remaining_bytes_to_load = min(sizeof(UInt), n)

    if remaining_bytes_to_load == 8
        return ntoh(unsafe_load(Ptr{UInt64}(pointer(s))))
    elseif  remaining_bytes_to_load == 7
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{Bits56}(s |> pointer))))
    elseif  remaining_bytes_to_load == 6
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{Bits48}(s |> pointer))))
    elseif  remaining_bytes_to_load == 5
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{Bits40}(s |> pointer))))
    elseif  remaining_bytes_to_load == 4
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{UInt}(s |> pointer))))
    elseif  remaining_bytes_to_load == 3
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{Bits24}(s |> pointer))))
    elseif  remaining_bytes_to_load == 2
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{UInt16}(s |> pointer))))
    else
        return ntoh(Base.zext_int(UInt, unsafe_load(Ptr{UInt8}(s |> pointer))))
    end
end

function check_boundary(s::String)::UInt
    ptrs  = pointer(s)
    if (UInt(ptrs) & 0xfff) > 0xff8
        return all_lens(s)
    else
        return ntoh(unsafe_load(Ptr{UInt64}(ptrs)))
    end
end

unsafe_loadbits(s) = s |> pointer |> Ptr{UInt} |> unsafe_load |> ntoh
leading_bitsdot(s) = leading_bits.(s)
leading_bits_with_fast_pathdot(s) = leading_bits_with_fast_path.(s)
unsafe_loadbitsdot(s) = unsafe_loadbits.(s)
load_bits_with_paddingdot(s) = load_bits_with_padding.(s)
copy_bits_fast_pathdot(s) = copy_bits_fast_path.(s)
all_lensdot(s) = all_lens.(s)
check_boundarydot(s) = check_boundary.(s)


text:
"""
    setup_laser(laser, units; τ=nothing, kwargs...)

Initialize the specified `laser` using the default parameteres in the specified
`units`. Supported units are
- `:SI`: Values are in SI, the numeric types are regular numbers (`Float64`)
- `:SI_unitful`: Values are in SI, but using [Unitful.jl](https://github.com/PainterQubits/Unitful.jl)
- `:atomic`: Values are in atomic units, the numeric types are regular numbers (`Float64`)
- `:atomic_unitful`: Values are in atomic units, but using [UnitfulAtomic.jl](https://github.com/sostock/UnitfulAtomic.jl)

The keyword arguments can be used to give specific values to the parameters instead
of using the defaults.
You can specify parameteres such as the wavelength `λ` (default `0.8u"μm"`),
beam waist `w₀` (default `58.0u"μm"`) and the normalized vector potential `a₀` (default `1`).

You can also use equivalent descriptions, such as passing the angular frequency `ω`
or the wave number `k` instead of `λ` and the [Rayleigh range](https://en.wikipedia.org/wiki/Rayleigh_length)
`z_R` instead of `z₀`.
The duration of the pulse (assuming Gaussian temporal profile) can be given via `τ`
(default `18.0u"fs"`).

You can also specify dimensionless parameteres such the initial phase (`ϕ₀`)
and the polarization (`ξx` and `ξy`).
The laser amplitude can be also specified through `E₀`.

By default the laser propagates along the `Oz` direction and osillates along the `Ox` direction.
If you want to change that, you can use the `propagation_dir` and `orientation` arguments.
See [Laser Geometry](@ref) for more details.
See also the docomuntation for each laser type for more details on the additional supported arguments.
"""
function setup_laser(laser, units; kwargs...)
    _λ, _w₀ = common_parameters(Val(units))
    constants = FundamentalConstants(units)
    c = constants.c

    if haskey(kwargs, :λ)
        λ = kwargs[:λ]
    elseif haskey(kwargs, :ω)
        ω = kwargs[:ω]
        λ = 2π*c / ω
    elseif haskey(kwargs, :k)
        k = kwargs[:k]
        λ = 2π/k
    else
        λ = _λ
    end

    if haskey(kwargs, :w₀)
        w₀ = kwargs[:w₀]
    elseif haskey(kwargs, :z_R)
        z_R = kwargs[:z_R]
        k = 2π / λ
        w₀ = √(2z_R / k)
    else
        w₀ = _w₀
    end

    if haskey(kwargs, :a₀)
        a₀ = kwargs[:a₀]
    elseif haskey(kwargs, :E₀)
        E₀ = kwargs[:E₀]
        q = abs(constants.q)
        mₑ = constants.mₑ
        ω = 2π * c / λ
        a₀ = E₀ * q / (mₑ * c * ω)
    else
        a₀ = 1
    end

    τ = get(kwargs, :τ, duration(Val(units)))
    t₀ = get(kwargs, :t₀, zero(τ))
    z₀ = get(kwargs, :z₀, zero(λ))

    @debug "profile"
    if haskey(kwargs, :profile)
        prof = kwargs[:profile]
        if prof isa Type
            if prof isa Type{<:GaussProfile} || prof isa Type{<:Cos²Profile}
                profile = prof(;τ, t₀, z₀)
            elseif prof isa Type{<:QuasiRectangularProfile}
                Δt = get(kwargs, :Δt, 10τ)
                profile = prof(;τ, t₀, z₀, Δt)
            else
                profile = prof()
            end
        else
            profile = prof
        end
    else
        @debug "default profile"
        profile = GaussProfile(;τ, z₀)
    end

    others = Dict{Symbol,Any}()
    processed = [:λ, :ω, :k, :w₀, :z_R, :a₀, :E₀, :τ, :t₀, :z₀, :profile]
    for (k,v) in kwargs
        if k ∉ processed
            push!(others, k=>v)
        end
    end

    laser(units; λ, w₀, a₀, profile, others...)
end

common_parameters(::Val{:SI_unitful}) = 0.8u"μm", 58.0u"μm"
common_parameters(::Val{:SI}) = ustrip.(u"m", common_parameters(Val(:SI_unitful)))
common_parameters(::Val{:atomic_unitful}) = auconvert.(common_parameters(Val(:SI_unitful)))
common_parameters(::Val{:atomic}) = ustrip.(common_parameters(Val(:atomic_unitful)))

duration(::Val{:SI_unitful}) = 18.0u"fs"
duration(::Val{:SI}) = ustrip.(u"s", duration(Val(:SI_unitful)))
duration(::Val{:atomic_unitful}) = auconvert.(duration(Val(:SI_unitful)))
duration(::Val{:atomic}) = ustrip.(duration(Val(:atomic_unitful)))


text:
module ACE

using Reexport

# external imports that are useful for all submodules
include("imports.jl")

@extimports
@baseimports

include("auxiliary.jl")

include("rotations3d.jl")
using ACE.Rotations3D

include("states.jl")
include("properties.jl")

include("prototypes.jl")

# basic polynomial building blocks
include("polynomials/sphericalharmonics.jl")
include("polynomials/transforms.jl"); @reexport using ACE.Transforms
include("polynomials/orthpolys.jl"); @reexport using ACE.OrthPolys

# The One-particle basis is the first proper building block
include("oneparticlebasis.jl")

# three specific 1p-bases that are useful
# TODO: species basis should be moved into the atomistic modelling toolkit
include("Ylm1pbasis.jl")
include("Rn1pbasis.jl")

include("product_1pbasis.jl")

include("sparsegrids.jl")

# TODO -> move elsewhere!!!
# include("rpi/rpi_degrees.jl")

# leave this for much later ...
# include("grapheval.jl")

# the permutation-invariant basis: this is a key building block
# for other bases but can also be a useful export itself
include("pibasis.jl")

include("symmbasis.jl")

# include("pipot.jl")

# # orthogonal basis
# include("orth.jl")
#
# lots of stuff related to random samples:
#  - random configurations
#  - random potentials
#  ...
include("random.jl")
@reexport using ACE.Random

include("utils.jl")
@reexport using ACE.Utils

# include("utils/importv5.jl")
# include("compat/compat.jl")
# include("export/export.jl")

include("testing/testing.jl")
include("testing/wigner.jl")

# - bond model
# - pure basis
# - real basis
# - regularisers
# - descriptors
# - random potentials

end # module