text:
using Pkg #hideall
Pkg.activate("_literate/using-JuMP/Project.toml")
Pkg.instantiate()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

md"""
## Overview

This tutorial will demonstrate how to use the `JuMP` package in Julia to solve linear optimization problems. It draws heavily from [this tutorial](https://github.com/east-winds/power-systems-optimization/blob/master/Notebooks/02-Anatomy-of-a-Model.ipynb) by Jesse D. Jenkins and Michael R. Davidson.

[JuMP](https://jump.dev/) ("**Ju**lia for **M**athematical **P**rogramming") is an open-source Julia package that adds functionality for formulating and solving a variety of optimization problems. One advantage of JuMP is that its syntax matches the typical mathematical formalism used to specify optimization problems. We will use JuMP in this class for our optimization work.

\note{Make sure that you take a look at the [`JuMP` documentation](https://jump.dev/JuMP.jl/stable/) whenever you have a question or want to find out how to do something that we don't discuss in any of our notebooks (or how to do it better!).}
"""

md"""
## Setup

Here we will outline the basic steps for configuring JuMP, though you can also refer to the official [Installation Guide](https://jump.dev/JuMP.jl/stable/installation/#Installation-Guide).

If `JuMP` is not already in your environment (it will be for any of your assignments, but may not be if you're doing something independently), you will need to install it. You will also need to [select a solver and install the relevant package](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers). Some of these are commercial, while others are open source. Solvers are also not typically universal, as different types of optimization problems use different algorithms, so be aware of what problem you're trying to solve instead of just blindly copying code from one task to another.

For example, for the linear programming example, we will use the Clp ([COIN-OR Linear Programming](https://github.com/coin-or/Clp)) solver via the [`Clp`](https://github.com/jump-dev/Clp.jl) package. As seen on the solver table, Clp is open source (via the EPL, or [Eclipse Public License](https://opensource.org/licenses/EPL-1.0)) and can solve linear programs (LP). For the mixed-integer linear programming example, we will instead of the Cbc ([COIN-OR Branch and Cut](https://github.com/coin-or/Cbc)) solver via the [`Cbc`](https://github.com/jump-dev/Cbc.jl) package. This is also open-source (it comes from the same group as Clp) and can solve mixed-integer or binary linear programs.
"""

md"""
## Linear Programming Example: How Many Widgets Should A Factory Produce?

### Defining The Problem

Suppose we own a factory that can produce two types of widgets:
* Widget A generates a profit of $p_A = \$100$ per widget; and
* Widget B generates a profit of $p_B = \$75$ per widget.
Let $x$ be the number of units of widget A that we want to produce, and $y$ the number of units of widget B. Our goal is to *maximize* our total profit $p_Ax + p_By$. This is the **objective function**. We express this objective using the equation $$\max 100x + 75y.$$

Further, we cannot build a negative number of either type of widget. This is known as a *non-negativity constraint*, and can be expressed as
\begin{align}
x &\geq 0\\
y &\geq 0
\end{align}
"""

md"""
This isn't a very interesting problem yet! We would simply build as much of both widgets as we could, because there are no constraints on our ability to produce. To make this more realistic, let's suppose that both widgets are produced using the same raw material $M$, of which we can only procure 300 units. Then, if:
* Widget A requires 40 units of $M$ per widget, and
* Widget B requires 20 units of $M$ per widget,
we arrive at the following material **constraint**:
$$40x + 20y \leq 300.$$

But we might have another constraint: time! Each widget may take a different amount of labor to produce. For example, let's say that
* Widget A takes 6 hours to produce and
* Widget B takes 12 hours to produce.
Further, there are only 80 hours per work that can be allocated to widget production. This becomes the time constraint
$$6x + 12y \leq 80.$$

Consolidating equations (1)--(4) gives us the following **constrained optimization problem**:
\begin{align}
& &\max 100x + 75y\\
&\text{subject to} & \\
& & x \geq 0\\
& & y \geq 0 \\
& & 40x + 20y \leq 300\\
& & 6x + 12y \leq 80
\end{align}
"""
md"""
### Visualizing the Problem

Let's do some plotting to examine the geometry of our optimization problem. We can do this using the [`Plots`](https://docs.juliaplots.org/latest/) package in Julia.
"""

## set up objective function parameters and variables
pa = 100
pb = 75
a = range(0, 8, step=0.25)
b = range(0, 8, step=0.25)

## define objective function
f(a, b) = pa * a + pb * b
#-
## start plotting
using Plots
plotlyjs() # hide

contour(a,b,(a,b)->f(a,b),nlevels=15, c=:heat, linewidth=10, colorbar = false, contour_labels = true) # objective function contours
title!("Factory Optimization Problem") # add title
xaxis!("x=Widget A", lims=(0, maximum(a))) # add x-axis title and limits
yaxis!("y=Widget B", lims=(0, maximum(b))) # add y-axis title and limits
xticks!(0:maximum(a)) # set x-axis ticks
yticks!(0:maximum(b)) # set y-axis ticks
areaplot!(a[a.<=11], (300 .- 40*a)./20, legend=false, opacity=0.3) # plot materials constraint feasible region
areaplot!(a[a.<=8], (80 .- 6*a)./12, legend=false, opacity=0.3) # plot time constraint feasible region
savefig(joinpath(@OUTPUT, "lp-visualize.png")) # hide

# \fig{lp-visualize.png}

md"""
We can see exactly where the solution will be, at the intersection of the feasible regions imposed by the two constraints (note how this is related to the gradient of the objective function)!  Let's now use JuMP to identify the location of this point (though we could also solve for it using linear algebra).
"""

md"""
### Solving This Problem Using JuMP

#### Setting Up the Model and Solver

First, we should load `JuMP` and `Clp`.
"""
using JuMP
using Clp

md"""
To solve our problem, first we need to **define the model**. The *model object* has lots of attributes, including the variables, constraints, solver options, etc. We create a new model using the `Model()` function. Since we are using the `Clp` solver, we need to tell JuMP to use the `Clp.Optimizer` solver function.
"""

factory_model = Model(Clp.Optimizer)

md"""
There are a bunch of attributes and options that we could set, but we won't in this example. If needed, look at the `Clp.jl` [documentation](https://github.com/jump-dev/Clp.jl).
"""

md"""
#### Define Variables
Decision variables ($x$ and $y$ in this case) in JuMP are defined using the `@variable` macro. The first argument passed to `@variable()` is the model object, in this case, `factory_model`, and the second argument are bounds on that variable, created using `>=` and `<=`. JuMP will interpret the bound specification to obtain the variable name. In this case, our only bounds directly on the variables are the non-negativity constraints.
"""

@variable(factory_model, x >= 0)
#-
@variable(factory_model, y >= 0)

md"""
If we had a free (or unbounded) variable $z$, we could declare that variable using `@variable(model, z)`. JuMP also requires unique names for each variable, or it will throw an error. This is one place where it's nice that Julia lets us use sub- and superscripts in variable names!

If we did want to modify the bounds after defining the variable, we could do so using the `set_lower_bound` and `set_upper_bound` functions, or we could remove them using `delete_lower_bound` and `delete_upper_bound`.

Finally, if we want to see all of the variables associated with a model, we can use the `all_variables` function to obtain an array.
"""
all_variables(factory_model)

md"""
#### Define Constraints
When defining variables, we were able to declare constraints on their values by specifying upper and lower bounds. However, we also have other constraints, which involve multiple decision variables. These are specified using the `@constraint` macro. Unlike variables, we also need to pass names for each constraint. We will use `time` for the time constraint and `materials` for the materials constraint. These names must be unique.
"""
@constraint(factory_model, time, 6x + 12y <= 80) # specify the time constraint
#-
@constraint(factory_model, materials, 40x + 20y <= 300) # materials constraint

md"""
#### Define Objective Function
So far, we've defined the feasible region of the decision-variable domain by setting the constraints. But we need to specify our objective function to know what we are trying to minimize or maximize over this region. We define the objective function using the `@objective` macro. In addition to specifying the model objective and the function, we need to tell JuMP whether we want to minimize or maximize.
"""
@objective(factory_model, Max, 100x + 75y)

md"""
Now, let's look at the model specification.
"""
print(factory_model)
md"""
We can also get a nicely styled LaTeX version, which might be useful.
"""
latex_formulation(factory_model)

md"""
We won't go into detail here, but there are other ways to define the model, which are detailed in the [JuMP documentation](https://jump.dev/JuMP.jl/stable/). For example, we can specify multiple variables using [`@variables`](https://jump.dev/JuMP.jl/stable/manual/variables/#variables). Similarly, we can use [`@constraints`](https://jump.dev/JuMP.jl/stable/manual/constraints/#The-@constraints-macro) to define multiple constraints at once. Or we can use [loops](https://jump.dev/JuMP.jl/stable/tutorials/Getting%20started/variables_constraints_objective/#Constraints-in-a-loop) to define multiple constraints or constraints involving many variables. We can also specify the model in [vectorized syntax](https://jump.dev/JuMP.jl/stable/tutorials/Getting%20started/variables_constraints_objective/#Vectorized-syntax), which is similar to how linear programs are specified in MATLAB.
"""

md"""
#### Solve the Model
Now it's time to solve the model and find the optimal values $(x^*, y^*)$. Since we specified the solver when we initialized `factory_model`, all we have to do is call the `optimize!` function.
"""
optimize!(factory_model)

md"""
#### Query the Solution
To find the optimal values of our decision variables, we need to query the values of the variables using `value.()`. We use `value.()` (the vectorized version of `value()`) because JuMP stores decision variables differently depending on their number and how they were defined. Uses the dot-syntax here works with any model specification, while the plain `value()` will not work if a queried decision variable is stored as a vector.
"""
value.(x)
#-
value.(y)

md"""
So we can see that our optimal inputs are $$(x^*, y^*) = (5.56, 3.89)$$ (and we'll pretend that we can manufacture and sell parts of widgets).
"""
md"""
#### Visualize the Solution

Let's take our previous plot and add the solution point to make sure that we got the solution we expected.
"""

## replot the previous plot
contour(a,b,(a,b)->f(a,b),nlevels=15, c=:heat, linewidth=10, colorbar = false, contour_labels = true) # objective function contours
title!("Factory Optimization Problem") # add title
xaxis!("x=Widget A", lims=(0, maximum(a))) # add x-axis title and limits
yaxis!("y=Widget B", lims=(0, maximum(b))) # add y-axis title and limits
xticks!(0:maximum(a)) # set x-axis ticks
yticks!(0:maximum(b)) # set y-axis ticks
areaplot!(a, (300 .- 40*a)./20, legend=false, opacity=0.3) # plot materials constraint feasible region
areaplot!(a, (80 .- 6*a)./12, legend=false, opacity=0.3) # plot time constraint feasible region

## now we plot the solution that we obtained
scatter!([value.(x)],[value.(y)], markercolor="blue")
savefig(joinpath(@OUTPUT, "lp-sol-visualize.png")) # hide

# \fig{lp-sol-visualize.png}

md"""
The solution $$(x^*, y^*)$$ is exactly where we deduced it would be geometrically.
"""

md"""
#### Other Stuff We Can Do
We can also use `value.()` to evaluate our constraints without manually using the equations.
"""
value.(time)
#-
value.(materials)

md"""
What if we also want the optimal objective value? We can obtain this using `objective_value()`.
"""
objective_value(factory_model)

md"""
We could also define other expressions via the [`@expression` macro](https://jump.dev/JuMP.jl/stable/expressions/#JuMP.@expression) as functions of the decision variables and evaluate those. For example, let's say that we wanted to know the total number of widgets we'd produce under our optimal allocation of resources.
"""
@expression(factory_model, total_widgets, x+y)
#-
value.(total_widgets)

md"""
#### Dual Solutions
We can identify if our model has a dual solution by calling `has_duals()`.
"""
has_duals(factory_model)

md"""
If we want to know the dual solution associated with a constraint, we use the `shadow_price()` function.
"""
shadow_price(time)
#-
shadow_price(materials)

md"""
If the binding constraint was a variable bound, we could also query that shadow price by calling `reduced_cost()` on the variable.
"""
reduced_cost(x)
#-
reduced_cost(y)

md"""
In this case, the relevant shadow prices are zero because the optimum is in the interior of the domain. If we had added a strong enough upper bound on the value(s) of one or both of our decision variables (say, x <= 4), then this would be non-zero.
"""


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 7, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function parseInput(input)
    rules = map(input) do rule
        (bag, containRules) = split(rule, " bags contain ")
        contains = map(split(containRules, ", ")) do rule
            m = match(r"^(\d+) (.*) bags?.?$", rule)

            if (isnothing(m))
                return nothing
            end

            return (parse(Int, m.captures[1]), string(m.captures[2]))
        end

        return (
            bag = bag,
            contains = filter(!isnothing, contains)
        )
    end

    return map(rule -> rule.bag => rule.contains, rules) |> Dict
end

function canContainShinyGold(rules, bag)
    if (!haskey(rules, bag))
        return false
    end
    # @show rules[bag]

    if "shiny gold" in last.(rules[bag])
        return true
    end

    return any(bag -> canContainShinyGold(rules, last(bag)), rules[bag])
end

function solve(::Puzzle{2020, 7, 1}, input)
    rules = parseInput(input)

    return count(
        ((bag, contains),) -> canContainShinyGold(rules, bag),
        rules
    )
end

function amountOfBagsIn(rules, bag)
    bags = rules[bag]

    if (length(bags) == 0)
        return 0
    end

    map(
        ((amount, name),) -> amount * (1 + amountOfBagsIn(rules, name)),
        bags
    ) |> sum
end

function solve(::Puzzle{2020, 7, 2}, input)
    rules = parseInput(input)
    amountOfBagsIn(rules, "shiny gold")
end

@testset "Day 7" begin
TEST_INPUT = split("""
light red bags contain 1 bright white bag, 2 muted yellow bags.
dark orange bags contain 3 bright white bags, 4 muted yellow bags.
bright white bags contain 1 shiny gold bag.
muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.
shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.
dark olive bags contain 3 faded blue bags, 4 dotted black bags.
vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.
faded blue bags contain no other bags.
dotted black bags contain no other bags.
""", '\n', keepempty=false)
    @test(solve(Puzzle(7, 1), TEST_INPUT) == 4)
    @test(solve(Puzzle(7, 1), input(Puzzle(7, 1))) == 144)

TEST_INPUT_2 = split("""
shiny gold bags contain 2 dark red bags.
dark red bags contain 2 dark orange bags.
dark orange bags contain 2 dark yellow bags.
dark yellow bags contain 2 dark green bags.
dark green bags contain 2 dark blue bags.
dark blue bags contain 2 dark violet bags.
dark violet bags contain no other bags.
""", '\n', keepempty=false)
    @test(solve(Puzzle(7, 2), TEST_INPUT_2) == 126)
    @test(solve(Puzzle(7, 2), input(Puzzle(7, 2))) == 5956)
end

ext: jl


text:
lines = readlines()

function solve(use_aim=false)
    depth, position, aim = 0, 0, 0
    for line in lines
        word, num = split(line)
        num = parse(Int, num)
        if word == "forward"
            position += num
            depth += use_aim ? aim * num : 0
        else
            if use_aim
                aim += word == "down" ? num : -num
            else
                depth += word == "down" ? num : -num
            end
        end
    end
    return depth * position
end
println(solve(false))
println(solve(true))

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
@testset "566.reshape-the-matrix.jl" begin
    @test matrixReshape([[1, 2], [3, 4]], 1, 4) == [[1, 2, 3, 4]]
    @test matrixReshape([[1, 2], [3, 4]], 2, 4) == [[1, 2], [3, 4]]
end


text:
"""
Reads metadata from NITF file

    nitfinfo(filename::AbstractString)
    :param filename: string pointing to filename

    :return nitf_meta: output metadata
"""
function nitinfo(filename::AbstractString)
    # get details about the file to read
    fileDetails = get_file_details(filename);

    # ensure the file is NITF and get its version
    valid, nitf_version = isnitf(fileDetails["filename"]);

    if (!valid)
        error("The image is not a NITF")
    elseif (!is_supported_nitf_version(nitf_version))
        error(string("This reader does not support NITF version: ", nitf_version );
    end

    # open the file
    fid = open(fileDetails["filename"], "r");

    # read the metadata
    if nitf_version == "2.1"
        nitf_meta = nitfparse21(fid);
    elseif nitf_version == "2.0"
        nitf_meta = nitfparse20(fid);
    end

    # close the file
    close(fid);

    # set up standard data elements

end

"""
Reads NITF file details to extract metadata

    get_file_details(filename::AbstractString)
    :param filename: string pointing to filename

    :return details: Dict containing details
"""
function get_file_details(filename::AbstractString)
    fid = open(filename, "r");
    info = stat(info);

    details = Dict("Filename" => filename, "FileModDate" => info.mtime,
                   "FileSize" => info.size);
    close(fid)
end

"""
Check if file is a NITF and returns the version

    isniti(filename)
    :param: filename: string pointing to filename

    :return tf: boolean if valid
    :return nitf_version: version number
"""
function isnitf(filename::AbstractString)
    fid = open(filename, "r");

    # get first conditional NITF header fields and inspect the first for the
    # NITF version
    fhdr = map(Char, readbytes(fid, 324));
    close(fid)

    if length(fhdr) == 324 && join(fhdr[1:9]) == "NITF02.10"
        tf = true;
        nitf_version = "2.1";
    elseif length(fhdr) == 324 && join(fhdr[1:9]) == "NITF02.00"
        tf = true;
        nitf_version = "2.0";
    elseif length(fhdr) == 324 && join(fhdr[1:9]) == "NITF01.10"
        tf = true;
        nitf_version = "1.1";
    elseif length(fhdr) == 324 && join(fhdr[1:9]) == "NSIF01.00"
        tf = true;
        # It's an NSIF 1.0 file which translates to NITF2.1
        nitf_version = "2.1";
    else
        nitf_version = "UNK";
        tf = false;
    end

    return tf, nitf_version
end

"""
Checks if reader supports NITF version

    is_supported_nitf_version(nitf_version);
    :param nitf_version: string
"""
function is_supported_nitf_version(nitf_version)
    tf = nitf_version == "2.1" || nitf_version == "2.0";

    return tf
end


text:
# Make a private registry for my Julia packages.

using Pkg

# Something in here seems to mess up my general environment
# The above seems to infect c:/Users/Mark Nahabedian/.julia/environments/v1.6/

# So activate a temporary one to try to preven that:

Pkg.activate(mktempdir())

using UUIDs: uuid4

############################################################
# Find packages in ~/.julia/dev that are remoted to GitHub:

using GitHub

#=
Now this works:  GitHub.references(VL; auth=github_auth)

function working_refertences(repo; auth)
    refs, _ = GitHub.gh_get_paged_json(GitHub.DEFAULT_API,
                                       "/repos/$(GitHub.name(VL))/git/refs";
                                       auth=github_auth)
    map(GitHub.Reference, refs)
end
=#

using TOML
import Base64

println("Enter github authentication token:")
github_auth = GitHub.authenticate(readline())

github_me = GitHub.whoami(; auth=github_auth)

github_repos, _ =  GitHub.repos(github_me; auth=github_auth)

# Map from GitHub.Repo to parsed Project.toml file:
parsed_project_files = Dict{GitHub.Repo, Dict{String, Any}}()

############################################################
# Find my Julia projects on GitHub:

# List my Julia repositories:
map(
    filter(github_repos) do r
        try
            if r.language != "Julia"
                return false
            end
            content_object = GitHub.file(r, "Project.toml"; auth=github_auth)
            @assert content_object.encoding == "base64"
            proj = TOML.parse(String(Base64.base64decode(content_object.content)))
            # error if missing:
            parsed_project_files[r] = proj
            proj["compat"]["julia"]
            true
        catch e
            println("*** $(r.name): $e")
            false
        end
    end) do r
        (r.name, r.fork)
    end

# Cherry-picked from the above results:
packages = [
    "AnotherParser.jl"
    "DXFutils.jl"
    "NahaJuliaLib.jl"
    "NativeSVG.jl"    # Not in above results because of missing compat section
    "PanelCutting.jl"
    "PlutoTool.jl"
    "ShaperOriginDesignLib"
    # "Unification.jl"
    "VectorLogging.jl"
    # "WebBasedWorkspace"
]

cherry_picked_repos = filter(github_repos) do r
    r.name in packages
end

############################################################
# Creating a registry file:

Pkg.add("RegistryTools")
using RegistryTools

function ensure_registry()
    registry_file = joinpath(@__DIR__, "Registry.toml")
    if isfile(registry_file)
        return registry_file
    end
    registry_data = RegistryTools.RegistryData(
        "NahaJuliaRegistry",
        uuid4(),
        repo=nothing,
        description="A registry of Mark Nahabedian's otherwise unregistered Julia packages.")
    RegistryTools.write_registry(registry_file, registry_data)
    registry_file
end

ensure_registry()

############################################################
# Choosing among reference versions

ACCEPTABLE_OBJECT_TYPES = Set{String}([
    "tag"
])

REJECTED_OBJECT_TYPES = Set{String}([
])

function objecttype(ref::GitHub.Reference)::String
    if haskey(ref.object , "type")
        ref.object["type"]
    else
        ""
    end
end

function acceptable_object_type(ref::GitHub.Reference):Bool
    acceptable_object_type(objecttype(ref))
end

function acceptable_object_type(ot::AbstractString)::Bool
    @assert !('/' in ot)
    if ot in ACCEPTABLE_OBJECT_TYPES
        true
    else
        push!(REJECTED_OBJECT_TYPES, ot)
        false
    end
end

function ref_version(ref::GitHub.Reference)::VersionNumber
    s = last(split(ref.ref, "/"))
    try
        VersionNumber(lstrip(s) do c
                          isletter(c) || (c in ['-', '_'])
                      end)
    catch e
        VersionNumber("0")
    end
end

function best_version(repo::GitHub.Repo)
    refs, _ = GitHub.references(repo; auth=github_auth)
    versions = map(filter(acceptable_object_type, refs)
                   ) do ref
                       (ref, ref_version(ref))
                   end
    if length(versions) > 0
        first(last(sort(versions; by=(v ->v[2]))))
    else
        nothing
    end
end

println("\n\nBest versions of cherry_picked_repos:")
for repo in cherry_picked_repos
    # local refs, _ = GitHub.references(repo; auth=github_auth)
    println("$(repo.name): \t",
            #=
            join(map(b -> b.ref,
                     filter(refs) do r
                         objecttype(r) in accept_types
                     end),
                 ", ")
            =#
            best_version(repo)
            )
end

for repo in cherry_picked_repos
    vpretty(ref) = ref === nothing ? nothing : ref.ref
    println(repo.name, " \t", vpretty(best_version(repo)))
end

############################################################
# Make clean clones

# LocalRegistry requires that each package have a clean checkout for
# development.  I have no clue why it needs to modify the package.

# Pkg.develop would clone the repository for us, but there is noway to
# temporarily override where it putrs the local copy.  Since I have
# active working copies in ~/.julia.dev, I need to create a directory
# for package clones and clone the packages myself.

const PACKAGE_STAGING_DIR = joinpath(@__DIR__, "package_staging")

function package_staging_dir()
    if !isdir(PACKAGE_STAGING_DIR)
        mkdir(PACKAGE_STAGING_DIR)
    end
    PACKAGE_STAGING_DIR
end

function clone(repo::GitHub.Repo)
    v = best_version(repo)
    cmd = [ "git", "clone", "--verbose", ]
    if v !== nothing
        v = last(split(v.ref, "/"))
        push!(cmd, "--branch")
        push!(cmd, v)
    end
    push!(cmd, repo.clone_url.uri)
    proc = run(pipeline(Cmd(Cmd(cmd);
                            dir=package_staging_dir());
                        stdout=stdout, stderr=stderr);
               wait=true)
    @assert proc.exitcode == 0
end

############################################################
# Registering packages:

Pkg.add("LocalRegistry")
using LocalRegistry

function repoPackageSpec(repo::GitHub.Repo)
    proj = parsed_project_files[repo]
    v = best_version(repo)
    more = Dict()
    if v != nothing
        more[:rev] = v.object["sha"]
    end
    PackageSpec(;name=proj["name"],
                uuid=proj["uuid"],
                url=repo.clone_url.uri,
                more...)
end

function repository_local_path(repo::GitHub.Repo)
    joinpath(PACKAGE_STAGING_DIR, repo.name)
end

function myregister(repo::GitHub.Repo)
    local_dir = repository_local_path(repo)
    @assert(isdir(local_dir), local_dir)
    LocalRegistry.register(local_dir;
                           registry=dirname(ensure_registry()),
                           commit=false,
                           push=false)
end

# map(myregister, cherry_picked_repos)

# How to keep LocalRegistry from creating first letter directories?
# package_relpath(pkg_name::String) in RegistryTools/4DGZp/src/types.jl
# is responsible for this.
# I don't see a way to control it so I guess I'm struck with it.

### Pkg.Registry.add gets an error because the Registry.toml file
### we create does not have a "repo" property.

max_stars_repo_name: MarkNahabedian/NahaJuliaRegistry

max_issues_repo_name: MarkNahabedian/NahaJuliaRegistry

max_forks_repo_name: MarkNahabedian/NahaJuliaRegistry


text:
"""
Utility to allow indexing into arrays using the letters 'a' to 'z'.
"""
to_index(c::Char) = (c - 'a') + 1


text:
module Dash
using DashBase
import HTTP, JSON2, CodecZlib, MD5
using Sockets
using Pkg.Artifacts
const ROOT_PATH = realpath(joinpath(@__DIR__, ".."))
#const RESOURCE_PATH = realpath(joinpath(ROOT_PATH, "resources"))
include("exceptions.jl")
include("Components.jl")
include("HttpHelpers/HttpHelpers.jl")

using .HttpHelpers

export dash, Component, Front, callback!,
enable_dev_tools!, ClientsideFunction,
run_server, PreventUpdate, no_update, @var_str,
Input, Output, State, make_handler, callback_context,
ALL, MATCH, ALLSMALLER, DashBase

include("Contexts/Contexts.jl")
include("env.jl")
include("utils.jl")
include("app.jl")
include("resources/application.jl")
include("handlers.jl")
include("server.jl")
include("init.jl")

@doc """
    module Dash

Julia backend for [Plotly Dash](https://github.com/plotly/dash)

# Examples
```julia

using Dash
app = dash(external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"]) do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        dcc_graph(id="graph",
            figure = (
                data = [(x = [1,2,3], y = [3,2,8], type="bar")],
                layout = Dict(:title => "Graph")
            )
        )

    end
end
callback!(app, Output("outputID", "children"), Input("graphTitle","value"), State("graphTitle","type")) do value, type
    "You've entered: '\$(value)' into a '\$(type)' input control"
end
callback!(app, Output("graph", "figure"), Input("graphTitle", "value")) do value
    (
        data = [
            (x = [1,2,3], y = abs.(randn(3)), type="bar"),
            (x = [1,2,3], y = abs.(randn(3)), type="scatter", mode = "lines+markers", line = (width = 4,))
        ],
        layout = (title = value,)
    )
end
run_server(app, HTTP.Sockets.localhost, 8050)
```

""" Dash

struct ComponentsInfo
    name::String
    version::VersionNumber
end
struct BuildInfo
    dash_version ::VersionNumber
    dash_renderer_version::VersionNumber
    embedded_components::Vector{ComponentsInfo}
end

function Base.show(io::IO, ::MIME"text/plain", info::BuildInfo)
    println(io, "Based on python `dash` version: ", info.dash_version)
    println(io, "\t`dash_renderer` version: ", info.dash_renderer_version)
    println(io, "Embedded components:")
    for comp in info.embedded_components
        println(io, "\t - `", comp.name, "` : ", comp.version)
    end
end
function build_info()
    dash_meta = load_meta("dash")
    renderer_meta = load_meta("dash_renderer")
    embedded = Vector{ComponentsInfo}()
    for name in dash_meta["embedded_components"]
        meta = load_meta(name)
        push!(embedded, ComponentsInfo(name, VersionNumber(meta["version"])))
    end

    return BuildInfo(
            VersionNumber(dash_meta["version"]),
            VersionNumber(renderer_meta["version"]),
            embedded
        )
end

@place_embedded_components
const _metadata = load_all_metadata()
function __init__()
   setup_renderer_resources()
   setup_dash_resources()
end

end # module


text:
module PGUManifoldMC

using Distributions
using Klara

import Base:
  cor,
  show,
  writemime

import Klara:
  RealLowerTriangular,
  codegen,
  ess,
  format_iteration,
  format_percentage,
  generate_empty,
  initialize!,
  reset!,
  sampler_state,
  tune!,
  tuner_state,
  variate_form

export
  GUM,
  MuvALSMMALAState,
  MuvAMSMMALAState,
  MuvGUMState,
  MuvPSMMALAState,
  ALSMMALA,
  AMSMMALA,
  ALSMMALAState,
  AMSMMALAState,
  PSMMALA,
  PSMMALAMCTune,
  PSMMALAMCTuner,
  PSMMALAState,
  cos_update!,
  ess,
  exp_decay,
  in_mahalanobis_contour,
  lin_decay,
  mahalanobis_distance,
  mahalanobis_update!,
  mala_only_update!,
  mod_update!,
  pow_decay,
  quad_decay,
  rand_exp_decay_update!,
  rand_lin_decay_update!,
  rand_pow_decay_update!,
  rand_quad_decay_update!,
  rand_update!,
  smmala_only_update!

include("stats/schedule.jl")

include("tuners/PSMMALAMCTuner.jl")

include("samplers/samplers.jl")

include("samplers/ALSMMALA.jl")
include("samplers/AMSMMALA.jl")

include("samplers/iterate/ALSMMALA.jl")
include("samplers/iterate/AMSMMALA.jl")

end # module


text:
using RandNLA, Test, LinearAlg

@testset "rdmd_test1" begin
end


text:
module JPleXL

  # Data for performance code

  include("./xCA.jl")
  include("./intDmat.jl")

  #
  # Functions for rigid-body rotation/translation
  #

  include("./center_of_mass.jl")
  include("./eulermat.jl")
  include("./tocm.jl")
  include("./compcart.jl")
  export compcart

end

max_stars_repo_name: mcubeg/JPleXL

max_issues_repo_name: mcubeg/JPleXL

max_forks_repo_name: mcubeg/JPleXL


text:
module LazyJSONGeoBenchmark

using LazyJSON
using JSON

function go()

# From https://geojson-maps.ash.ms
j = String(read("geo.json"))

for n in [1, 19]

n > 1 && println("")
n > 1 && println("Country name")

n > 1 && print("LazyJSON.jl:")
GC.gc()
@time for i in 1:n
    r = LazyJSON.parse(j)
    @assert r["features"][15]["properties"]["formal_en"] == "Republic of Palau"
end

n > 1 && print("JSON.jl:    ")
GC.gc()
@time for i in 1:n
    r = JSON.parse(j)
    @assert r["features"][15]["properties"]["formal_en"] == "Republic of Palau"
end

n > 1 && println("")
n > 1 && println("Map data")
n > 1 && print("LazyJSON.jl:")
GC.gc()
@time for i in 1:n
    r = LazyJSON.parse(j)
    r = r["features"][15]["geometry"]["coordinates"][6][1]
    @assert r[1][1] == 134.41651451900023
    for (x, y) in r
       @assert 134.2 < x < 134.5
       @assert 7.21 < y < 7.32
    end
end

n > 1 && print("JSON.jl:    ")
GC.gc()
@time for i in 1:n
    r = JSON.parse(j)
    r = r["features"][15]["geometry"]["coordinates"][6][1]
    @assert r[1][1] == 134.41651451900023
    for (x, y) in r
       @assert 134.2 < x < 134.5
       @assert 7.21 < y < 7.32
    end
end

end

end

end

LazyJSONGeoBenchmark.go()


text:
using Random; Random.seed!(1234)

include("src/data.jl")
using .DataUtil
include("src/forward_stepwise/selective_inference.jl")
using .SIforSFS
include("src/forward_stepwise/data_splitting.jl")
import .DS
using LinearAlgebra:I
using Distributions
using ProgressMeter

const N = 100
const β₀ = 0.0
const 𝛃 = (Float64)[0, 0, 0, 0, 0] # null model
const k = 3

# FPR experiment
function FPR_experiment(β₀, 𝛃, N, k; event=Active(), α=0.05, denominator_upper=100, iter_size=5)
    H0 = filter(i -> iszero(𝛃[i]), 1:length(𝛃)) # {j ∈ {1,...,d} | H_{0,j} is true}
    FPR = Array{Float64}(undef, length(H0), iter_size)
    @showprogress "[$(N), $(event)]: " for l = 1:iter_size
        numerator = zeros(Int, length(H0))
        denominator = zeros(Int, length(H0))
        for _ = 1:denominator_upper
            𝐗, 𝐲_obs = centering(make_dataset(β₀, 𝛃, N)...)
            𝐴_obs, p_values = parametric_SFS_SI(𝐗, 𝐲_obs, k, I; selection_event=event)
            for (j, p_value) in zip(𝐴_obs, p_values)
                for (i, j′) in enumerate(H0)
                    if j == j′ # H_{0,j} is true
                        numerator[i] += p_value < α ? 1 : 0
                        denominator[i] += 1
                        break
                    end
                end
            end
        end
        FPR[:, l] = numerator ./ denominator
    end
    return FPR
end

α = 0.05
zero_num = length(filter(x -> iszero(x), 𝛃))
Ns = [50, 100, 150, 200]
iter_size = 20
FPR_A = fill(NaN, zero_num, iter_size, length(Ns))
FPR_Ao = fill(NaN, zero_num, iter_size, length(Ns))
FPR_As = fill(NaN, zero_num, iter_size, length(Ns))
FPR_Aso = fill(NaN, zero_num, iter_size, length(Ns))
FPR_DS = fill(NaN, zero_num, iter_size, length(Ns))
for (i, N) in enumerate(Ns)
    FPR_A[:,:,i] = Main.FPR_experiment(β₀, 𝛃, N, k; event=Active(), α=α, iter_size=iter_size)
    FPR_As[:,:,i] = Main.FPR_experiment(β₀, 𝛃, N, k; event=ActiveSign(), α=α, iter_size=iter_size)
    FPR_Ao[:,:,i] = Main.FPR_experiment(β₀, 𝛃, N, k; event=ActiveOrder(), α=α, iter_size=iter_size)
    FPR_Aso[:,:,i] = Main.FPR_experiment(β₀, 𝛃, N, k; event=ActiveSignOrder(), α=α, iter_size=iter_size)
    FPR_DS[:,:,i] = DS.FPR_experiment(β₀, 𝛃, N, k; α=α, iter_size=iter_size)
end

using Statistics
FPR_A_mean = reshape(mean(FPR_A, dims=2), zero_num, length(Ns))
FPR_As_mean = reshape(mean(FPR_As, dims=2), zero_num, length(Ns))
FPR_Ao_mean = reshape(mean(FPR_Ao, dims=2), zero_num, length(Ns))
FPR_Aso_mean = reshape(mean(FPR_Aso, dims=2), zero_num, length(Ns))
FPR_DS_mean = reshape(mean(FPR_DS, dims=2), zero_num, length(Ns))

using Plots
Cs = [:red, :orange, :green, :blue, :purple]
p = Array{typeof(plot())}(undef, zero_num)
for j in 1:zero_num
    plot()
    plot!(1:length(Ns), FPR_A_mean[j,:], label="Homotopy", color=Cs[1], markershape=:circle)
    plot!(1:length(Ns), FPR_Ao_mean[j,:], label="Homotopy-H", color=Cs[2], markershape=:circle)
    plot!(1:length(Ns), FPR_As_mean[j,:], label="Homotopy-S", color=Cs[3], markershape=:circle)
    plot!(1:length(Ns), FPR_Aso_mean[j,:], label="Polytope", color=Cs[4], markershape=:circle)
    plot!(1:length(Ns), FPR_DS_mean[j,:], label="DS", color=Cs[5], markershape=:circle)
    # plot setting
    plot!(; ylims=(0, 0.3))
    plot!(; xlabel="Sample size", ylabel="False Positive Rate (FPR)")
    plot!(; legend=:topright)
    plot!(; guidefontsize=16, legendfontsize=12, tickfontsize=10)
    p[j] = plot!(; xticks=(1:length(Ns), Ns))
    savefig(p[j], "img/FPR/FPR$(j).pdf")
end


text:
# This file was generated, do not modify it. # hide
mach = machine(mdl, X2, y)
fit!(mach)
ŷ = predict(mach, X2)
round(rms(ŷ, y), sigdigits=4)


text:
using BinDeps
using Compat

@BinDeps.setup

sptk = library_dependency("libSPTK", aliases=["libSPTK", "SPTK-3"])

const version = "3.8.9"

github_root = "https://github.com/r9y9/SPTK"
arch = WORD_SIZE == 64 ? "x86_64" : "i686"
major = version[1]
provides(Binaries,
         URI("$(github_root)/releases/download/v$(version)/sptk-$(major)_mingw$(WORD_SIZE)_$(arch).zip"),
         sptk, unpacked_dir = "usr/lib", os = :Windows)

provides(Sources,
         URI("$(github_root)/archive/v$(version).tar.gz"),
         sptk,
         unpacked_dir="SPTK-$(version)")

prefix = joinpath(BinDeps.depsdir(sptk), "usr")
srcdir = joinpath(BinDeps.depsdir(sptk), "src", "SPTK-$(version)")

provides(SimpleBuild,
          (@build_steps begin
              GetSources(sptk)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `./waf configure --prefix=$prefix`
                  `./waf build`
                  `./waf install`
              end
           end), sptk, os = :Unix)

@BinDeps.install @compat Dict(:libSPTK => :libSPTK)


text:
using Test

function ReLU(x::Float64)
    result = max(0.0, x)
    return result
end

@test ReLU(1.0) == 1.0
@test ReLU(100.0) == 100.0

@test ReLU(0.0) == 0.0

@test ReLU(-1.0) == 0.0
@test ReLU(-10.0) == 0.0


text:
using DataFrames
using Test

m = test_page_model()
include("../src/components/TotalAdaptationCosts.jl")

add_comp!(m, TotalAdaptationCosts)

set_param!(m, :TotalAdaptationCosts, :pop_population, readpagedata(m, "test/validationdata/pop_population.csv"))
set_param!(m, :TotalAdaptationCosts, :ac_adaptationcosts_economic, readpagedata(m, "test/validationdata/ac_adaptationcosts_economic.csv"))
set_param!(m, :TotalAdaptationCosts, :ac_adaptationcosts_noneconomic, readpagedata(m, "test/validationdata/ac_adaptationcosts_noneconomic.csv"))
set_param!(m, :TotalAdaptationCosts, :ac_adaptationcosts_sealevelrise, readpagedata(m, "test/validationdata/ac_adaptationcosts_sealevelrise.csv"))

run(m)

# Generated data
adapt_cost = m[:TotalAdaptationCosts, :act_adaptationcosts_total]
adapt_cost_per_cap = m[:TotalAdaptationCosts, :act_percap_adaptationcosts]

# Recorded data
cost_compare = readpagedata(m, "test/validationdata/act_adaptationcosts_tot.csv")
cost_cap_compare = readpagedata(m, "test/validationdata/act_percap_adaptationcosts.csv")

@test adapt_cost ≈ cost_compare rtol=1e-3
@test adapt_cost_per_cap ≈ cost_cap_compare rtol=1e-6


text:
## --- File pointers

    fp = stdoutp()
    @test isa(fp, Ptr{StaticTools.FILE})
    @test stdoutp() == fp != 0
    fp = stderrp()
    @test isa(fp, Ptr{StaticTools.FILE})
    @test stderrp() == fp != 0
    fp = stdinp()
    @test isa(fp, Ptr{StaticTools.FILE})
    @test stdinp() == fp != 0

    fp = fopen(c"testfile.txt", c"w")
    @test isa(fp, Ptr{StaticTools.FILE})
    @test fp != 0
    @test fclose(fp) == 0

## -- Test low-level printing functions on a variety of arguments

    @test puts("1") == 0
    @test printf("2") >= 0
    @test putchar('\n') == 0
    @test printf("%s\n", "3") >= 0
    @test printf(4) == 0
    @test printf(5.0) == 0
    @test printf(10.0f0) == 0
    @test printf(0x01) == 0
    @test printf(0x0001) == 0
    @test printf(0x00000001) == 0
    @test printf(0x0000000000000001) == 0
    @test printf(Ptr{UInt64}(0)) == 0

## -- low-level printing to file

    fp = fopen("testfile.txt", "w")
    @test isa(fp, Ptr{StaticTools.FILE})
    @test fp != 0

    @test puts(fp, "1") == 0
    @test printf(fp, "2") == 1
    @test putchar(fp, '\n') == 0
    @test printf(fp, "%s\n", "3") == 2
    @test printf(fp, 4) == 0
    @test printf(fp, 5.0) == 0
    @test printf(fp, 10.0f0) == 0
    @test printf(fp, 0x01) == 0
    @test printf(fp, 0x0001) == 0
    @test printf(fp, 0x00000001) == 0
    @test printf(fp, 0x0000000000000001) == 0
    @test printf(fp, Ptr{UInt64}(0)) == 0

## -- High-level printing

    # Print AbstractVector
    @test printf(1:5) == 0
    @test printf((1:5...,)) == 0
    @test printf(fp, 1:5) == 0
    @test printf(fp, (1:5...,)) == 0
    a = MallocVector{Float64}(undef, 2); a[:] = 2
    @test print(a) == 0
    @test println(a) == 0
    @test print(fp, a) == 0
    @test println(fp, a) == 0
    @test free(a) == 0

    # Print AbstractArray
    @test printf((1:5)') == 0
    @test printf(rand(4,4)) == 0
    @test printf(fp, (1:5)') == 0
    @test printf(fp, rand(4,4)) == 0
    a = MallocMatrix{Float64}(undef, 1,2); a[:] = 2
    @test print(a) == 0
    @test println(a) == 0
    @test print(fp, a) == 0
    @test println(fp, a) == 0
    @test free(a) == 0

    # Print MallocString
    str = m"Hello, world! 🌍"
    @test print(str) === Int32(strlen(str))
    @test println(str) === Int32(0)
    @test print(fp, str) === Int32(strlen(str))
    @test println(fp, str) === Int32(0)
    @test printf(str) == Int32(strlen(str))
    @test printf(fp, str) == Int32(strlen(str))
    @test puts(str) == 0
    @test printf(m"%s \n", str) >= 0
    show(str)

    # Print StaticString
    str = c"Hello, world! 🌍"
    @test print(str) === Int32(strlen(str))
    @test println(str) === Int32(0)
    @test print(fp, str) === Int32(strlen(str))
    @test println(fp, str) === Int32(0)
    @test printf(str) === Int32(strlen(str))
    @test printf(fp, str) === Int32(strlen(str))
    @test puts(str) == 0
    @test printf(m"%s \n", str) >= 0
    show(str)

    # Print StringView
    sview = str[1:5]
    @test puts(sview) === Int32(0)
    @test puts(fp, sview) === Int32(0)
    @test printf(sview) === Int32(5)
    @test printf(fp, sview) === Int32(5)
    @test print(sview) === Int32(5)
    @test println(sview) === Int32(0)
    @test print(fp, sview) === Int32(5)
    @test println(fp, sview) === Int32(0)
    show(sview)

## --- Printing errors

    @test error(c"This is a test\n") == 0
    @test StaticTools.warn(c"This is a test\n") == 0
    msg = m"This is a test\n"
    @test error(msg) == 0
    @test StaticTools.warn(msg) == 0
    free(msg)

## --- Libc functions that just happen to print

    msg = m"pwd"
    @test StaticTools.system(msg) == 0
    @test StaticTools.system(c"pwd") == 0
    free(msg)

## ---  Wrap up

    @test newline() == 0
    @test newline(fp) == 0
    @test StaticTools.system(c"echo Enough printing for now!") == 0
    @test fclose(fp) == 0

## -- Reading from files

    name, mode = m"testfile.txt", m"r"
    fp = fopen(name, mode)
    @test isa(fp, Ptr{StaticTools.FILE})
    @test fp != 0
    @static if fp != 0
        str = MallocString(undef, 100)
        @test gets!(str, fp) != C_NULL
        @test strlen(str) == 2
        @test str[1] == UInt8('1')
        @test fseek(fp, -2, SEEK_CUR) == 0
        @test gets!(str, fp) != C_NULL
        @test strlen(str) == 2
        @test str[1] == UInt8('1')
        @test free(str) == 0
        @test frewind(fp) == 0

        @test getc(fp) === Int32('1')
        @test getc(fp) === Int32('\n')

        @test fclose(fp) == 0
    end
    # @static if fp != 0 && Sys.isapple()
    #     @test getchar(fp) === UInt8('1')
    #     @test getchar(fp) === UInt8('\n')
    #     @test fclose(fp) == 0
    # end
    @test free(name) == 0
    @test free(mode) == 0

## --- clean up

    rm("testfile.txt")


text:
int_rules_1_1_4_4 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.4.4*P(x)*(c*x)^m*(a*x^j+b*x^n)^p =#
    @apply_utils Antiderivative(~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([d = Denominator(~n)], d * Subst(Antiderivative((~x) ^ (d - 1) * ReplaceAll(SubstFor((~x) ^ ~n, ~Pq, ~x), (~x->begin
                                            #= none:1 =#
                                            #= none:3 =#
                                            (~x) ^ (d * ~n)
                                        end)) * (~a * (~x) ^ (d * ~j) + ~b * (~x) ^ (d * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / d))) <-- FreeQ([~a, ~b, ~j, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (RationalQ(~j, ~n) && (IntegerQ(~j / ~n) && LtQ(-1, ~n, 1))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * SubstFor((~x) ^ ~n, ~Pq, ~x) * (~a * (~x) ^ Simplify(~j / ~n) + ~b * ~x) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~j, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && IntegerQ(Simplify((~m + 1) / ~n))))))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ (Sign(~m) * Quotient(~m, Sign(~m))) * (~c * ~x) ^ Mod(~m, Sign(~m))) / (~x) ^ Mod(~m, Sign(~m))) * Antiderivative((~x) ^ ~m * ~Pq * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify((~m + 1) / ~n)) && (RationalQ(~m) && GtQ((~m) ^ 2, 1)))))))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((~c * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative((~x) ^ ~m * ~Pq * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && IntegerQ(Simplify((~m + 1) / ~n))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([g = GCD(~m + 1, ~n)], (1 / g) * Subst(Antiderivative((~x) ^ ((~m + 1) / g - 1) * ReplaceAll(~Pq, (~x->begin
                                                #= none:1 =#
                                                #= none:3 =#
                                                (~x) ^ (1 / g)
                                            end)) * (~a * (~x) ^ (~j / g) + ~b * (~x) ^ (~n / g)) ^ ~p, ~x), ~x, (~x) ^ g) <-- NeQ(g, 1)) <-- FreeQ([~a, ~b, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (IGtQ(~j, 0) && (IGtQ(~n, 0) && (IGtQ(~j / ~n, 0) && IntegerQ(~m))))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], With([Pqq = Coeff(~Pq, ~x, q)], (Pqq * (~c * ~x) ^ (((~m + q) - ~n) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * (~c) ^ ((q - ~n) + 1) * (~m + q + ~n * ~p + 1)) + Antiderivative((~c * ~x) ^ ~m * ExpandToSum((~Pq - Pqq * (~x) ^ q) - (~a * Pqq * (((~m + q) - ~n) + 1) * (~x) ^ (q - ~n)) / (~b * (~m + q + ~n * ~p + 1)), ~x) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x)) <-- GtQ(q, ~n - 1) && (NeQ(~m + q + ~n * ~p + 1, 0) && (IntegerQ(2 * ~p) || IntegerQ(~p + (q + 1) / (2 * ~n))))) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (Not(IntegerQ(~p)) && (IGtQ(~j, 0) && (IGtQ(~n, 0) && LtQ(~j, ~n)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (1 / (~m + 1)) * Subst(Antiderivative(ReplaceAll(SubstFor((~x) ^ ~n, ~Pq, ~x), (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        (~x) ^ Simplify(~n / (~m + 1))
                                    end)) * (~a * (~x) ^ Simplify(~j / (~m + 1)) + ~b * (~x) ^ Simplify(~n / (~m + 1))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~j, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n)))))))
    @apply_utils Antiderivative((~c * ~x) ^ ~m * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ (Sign(~m) * Quotient(~m, Sign(~m))) * (~c * ~x) ^ Mod(~m, Sign(~m))) / (~x) ^ Mod(~m, Sign(~m))) * Antiderivative((~x) ^ ~m * ~Pq * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~n / (~m + 1))) && (Not(IntegerQ(~n)) && GtQ((~m) ^ 2, 1)))))))
    @apply_utils Antiderivative((~c * ~x) ^ ~m * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((~c * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative((~x) ^ ~m * ~Pq * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n)))))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~c * ~x) ^ ~m * ~Pq * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && ((PolyQ(~Pq, ~x) || PolyQ(~Pq, (~x) ^ ~n)) && (Not(IntegerQ(~p)) && NeQ(~n, ~j)))
    @apply_utils Antiderivative(~Pq * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => Antiderivative(ExpandIntegrand(~Pq * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~j, ~n, ~p], ~x) && ((PolyQ(~Pq, ~x) || PolyQ(~Pq, (~x) ^ ~n)) && (Not(IntegerQ(~p)) && NeQ(~n, ~j)))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.4 Improper/1.1.4.4 P(x) (c x)^m (a x^j+b x^n)^p.jl


text:
# Model
include_model("hopper")

mb = 3.0 # body mass
ml = 0.3  # leg mass
Jb = 0.75 # body inertia
Jl = 0.075 # leg inertia

model = Hopper{Discrete, FixedTime}(n, m, d,
			   mb, ml, Jb, Jl,
			   0.25, g,
			   qL, qU,
			   uL, uU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

model_ft = free_time_model(model)

# Stair
function ϕ_func(model::Hopper, q)
	k = kinematics(model, q)
	if k[1] > 0.125
    	return @SVector [k[2] - 0.25]
	else
		return @SVector [k[2]]
	end
end

# Horizon
T = 81

# Time step
tf = 0.8
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] .= Inf #model_ft.uU
_uu[end] = 5.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -Inf # model_ft.uL
_ul[end] = 0.1 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
z_h = 0.25
q1 = [0.0, 0.5, 0.0, 0.5]
qm1 = [0.0625, 0.5 + z_h, 0.0, 0.5]
qM = [0.125, 0.5 + 2.5 * z_h, 0.0, 0.5]
qm2 = [0.1875, 0.5 + z_h + z_h, 0.0, 0.5]
qT = [0.25, 0.5 + z_h, 0.0, 0.5]
xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; q1],
		xT = [qT; qT])

ql1 = linear_interpolation(q1, qm1, 21)
ql2 = linear_interpolation(qm1, qM, 21)
ql3 = linear_interpolation(qM, qm2, 21)
ql4 = linear_interpolation(qm2, qT, 22)
q_ref = [ql1..., ql2[2:end]..., ql3[2:end]..., ql4[2:end]...]

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])
qp = [0.01; 0.01; 1000.0; 1.0]
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal(0.0 * [qp; qp]) for t = 1:T],
    [Diagonal([1.0e-1, 1.0e-1,
		1.0e-5 * ones(model_ft.nc)..., 1.0e-5 * ones(model_ft.nb)...,
		zeros(model_ft.m - model_ft.nu - model_ft.nc - model_ft.nb - 1)..., 0.0])
		for t = 1:T-1],
    [[qT; qT] for t = 1:T],
    [zeros(model_ft.m) for t = 1:T],
    10.0)

obj_contact_penalty = PenaltyObjective(1.0e5, model_ft.m - 1)

obj_velocity = velocity_objective(
    [Diagonal(0.0e-3 * ones(model_ft.nq)) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect([3]))

function l_foot_vel(x, u, t)
	J = 0.0

	q1 = view(x, 1:4)
	p1 = kinematics(model, q1)

	q2 = view(x, 4 .+ (1:4))
	p2 = kinematics(model, q2)

	v = (p2 - p1) ./ h

	if t > 15 && t < 66
		J += 1000.0 * (q1 - q_ref[t])' * (q1 - q_ref[t])
	end
	if t >= 50
		J += 10000.0 * (q1 - qT)' * (q1 - qT)
	end

	#
	# if t < 10 || t > 71
	# 	J += 1000.0 * v[1]^2.0
	# end

	return J
end

l_foot_vel(x) = l_foot_vel(x, nothing, T)
obj_foot_vel = nonlinear_stage_objective(l_foot_vel, l_foot_vel)

obj_ctrl_vel = control_velocity_objective(Diagonal([0.0 * ones(model_ft.nu);
	0.0 * ones(model_ft.nc + model_ft.nb);
	zeros(model_ft.m - model_ft.nu - model_ft.nc - model_ft.nb)]))

obj = MultiObjective([obj_tracking, obj_contact_penalty, obj_velocity, obj_foot_vel, obj_ctrl_vel])

# Constraints
include_constraints(["free_time", "contact", "stage"])
con_free_time = free_time_constraints(T)
con_contact = contact_constraints(model_ft, T)

p1_ref = kinematics(model, q1)
pT_ref = kinematics(model, qT)

function pinned1!(c, x, u, t)
    q = view(x, 1:4)
    c[1:2] = p1_ref - kinematics(model, q)
	nothing
end

function pinnedT!(c, x, u, t)
    q = view(x, 4 .+ (1:4))
	c[1:2] = pT_ref - kinematics(model, q)
	nothing
end

T_fix = 10
n_stage = 2
t_idx1 = vcat([t for t = 1:20]...)
t_idxT = vcat([(T - 40 + 1):T]...)

con_pinned1 = stage_constraints(pinned1!, n_stage, (1:0), t_idx1)
con_pinnedT = stage_constraints(pinnedT!, n_stage, (1:0), t_idxT)
con = multiple_constraints([con_free_time, con_contact, con_pinned1, con_pinnedT])#, con_loop])

# Problem
prob = trajectory_optimization_problem(model_ft,
           obj,
           T,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu,
           con = con)

# Trajectory initialization
x0 = configuration_to_state(q_ref) # linear interpolation on state
u0 = [[1.0e-3 * rand(model_ft.m-1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0),
	nlp = :ipopt,
	tol = 1.0e-4, c_tol = 1.0e-4, mapl = 5,
	max_iter = 2000,
	time_limit = 60)

@show check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
tf, t, h̄ = get_time(ū)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]
h̄ = mean(h̄)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model_ft,
	q,
	Δt = h̄[1],
	scenario = :vertical)

setobject!(vis["box1"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
	Vec(0.25, 0.5, 0.25)), MeshPhongMaterial(color = RGBA(0.5, 0.5, 0.5, 1.0)))
settransform!(vis["box1"], Translation(0.125, -0.25, 0))

hm = h̄
μm = model.μ
qm, um, γm, bm, ψm, ηm = q, u, γ, b, ψ, η

@save joinpath(@__DIR__, "hopper_stair.jld2") qm um γm bm ψm ηm μm hm
# @load joinpath(@__DIR__, "hopper_stair.jld2") qm um γm bm ψm ηm μm hm

max_issues_repo_name: thowell/DirectMotionPlanning


text:
include("raytracer.jl")
import .RT: Camera, render, ray_color
using Images

camera = Camera(1600, 900)
save("out.png", render(camera, ray_color))

max_stars_repo_name: hdtx/raytracer-comparison

max_issues_repo_name: hdtx/raytracer-comparison

max_forks_repo_name: hdtx/raytracer-comparison


text:
@testset "cvt" begin
    @printf("\n> cvt\n")

    # toInt
    @test Bb.toInt("1") == UInt64(1)
    @test Bb.toInt("10") == UInt64(2)
    @test Bb.toInt("11") == UInt64(3)

    # PGN2UINT
    a8uint = Bb.PGN2UINT["a8"]
    @test Bb.UINT2PGN[a8uint] == "a8"

    # PGN2INT
    @test Bb.PGN2INT["a8"] == 1
    @test Bb.INT2PGN[1] == "a8"

    # INT2UINT
    @test Bb.INT2UINT[1] == Bb.PGN2UINT["a8"]
    uint = Bb.PGN2UINT["a8"]
    @test Bb.UINT2INT[uint] == 1
end
p2u = Bb.PGN2UINT


text:
module Utilities
export hart_forcade, to_fractional, to_cartesian, is_periodic, into_cell, origin_centered
export into_voronoi, supercell, cell_parameters

using Crystals.Constants: default_tolerance
using Crystals.Structures: Crystal
using Crystals.SNF: smith_normal_form
using DataFrames: nrow
using NamedTuples: @NT
using ArgCheck
using Unitful
using Unitful: Dimensions, NoUnits
using MicroLogging
using DocStringExtensions
using Base.Iterators: cycle, take

"""
    eldimension(u::Any)

Dimension of the input, if it has one. Arrays have a dimension if their element
type has a uniform dimension.
"""
eldimension{T <: AbstractArray}(u::Type{T}) = eldimension(eltype(T))
eldimension{T, D, U}(::Type{Quantity{T, D, U}}) = D
eldimension{T <: Number}(::Type{T}) = Dimensions{()}
eldimension(u::DataType) = Dimensions{()}
eldimension(u::Any) = eldimension(typeof(u))

is_unitful(a::AbstractArray) = is_unitful(eltype(a))
is_unitful{T <: Number}(a::Type{T}) = Val{:unitless}()
is_unitful{T, D, U}(a::Type{Quantity{T, D, U}}) = Val{:unitful}()

""" Tuple holding Hart-Forcade transform """
const HartForcadeTransform = @NT(transform::Matrix, quotient::Vector)

"""
$(SIGNATURES)

Computes the cyclic group of a supercell with respect to a lattice. It makes it
possible to identify the class of periodically equivalent cell that a given
position within the supercell belongs to. The function returns a named tuple with the
transform and the quotient.

# Examples

```jldoctest
using Crystals, Unitful
fcc = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"
supercell = [0 2 2; 0 -4 2; -1 0 -2]
ht = hart_forcade(fcc, fcc * supercell)

println(ht)
println("Positions in supercell:")
for index in CartesianRange((ht.quotient...))
    position = inv(ht.transform) * [index[u] for u in eachindex(ht.quotient)]
    println("- ", ustrip.(position), " (", unit(eltype(position)), ")")
end

# output

Hart-Forcade transform
- transform (nm^-1): [-1.0 -1.0 1.0; -1.0 1.0 1.0; -1.0 1.0 3.0]
- quotient: [1, 2, 6]

Positions in supercell:
- [-1.0, 0.0, 0.0] (nm)
- [-2.0, 0.5, -0.5] (nm)
- [-0.5, 0.0, 0.5] (nm)
- [-1.5, 0.5, 0.0] (nm)
- [0.0, 0.0, 1.0] (nm)
- [-1.0, 0.5, 0.5] (nm)
- [0.5, 0.0, 1.5] (nm)
- [-0.5, 0.5, 1.0] (nm)
- [1.0, 0.0, 2.0] (nm)
- [0.0, 0.5, 1.5] (nm)
- [1.5, 0.0, 2.5] (nm)
- [0.5, 0.5, 2.0] (nm)
```
"""
function hart_forcade(lattice::AbstractMatrix, supercell::AbstractMatrix; digits::Integer=8)
    fractional = convert(Matrix{Int64}, round.(inv(lattice) * supercell, digits))

    snf, left, right = smith_normal_form(fractional)

    HartForcadeTransform(left * inv(lattice), diag(snf))
end

function Base.show(io::IO, ht::HartForcadeTransform)
    if unit(eltype(ht.transform)) ≠ NoUnits
        trans_unit = " ($(unit(eltype(ht.transform))))"
    else
        trans_unit = ""
    end
    println(io, "Hart-Forcade transform")
    println(io, "- transform", trans_unit, ": ", ustrip.(ht.transform))
    println(io, "- quotient: ", ustrip.(ht.quotient))
end

"""
     to_fractional(positions::AbstractArray, cell::AbstractMatrix)

Converts positions to fractional units (not Unitful). If the positions have no units, they
are already fractional units, and this operation is a no-op.
"""
function to_fractional(pos::AbstractArray, cell::AbstractMatrix)
    to_fractional(is_unitful(pos), pos, cell)
end
to_fractional(::Val{:unitless}, positions::AbstractArray, ::AbstractMatrix) = positions
function to_fractional(v::Val{:unitful}, positions::AbstractArray, cell::AbstractMatrix)
    to_fractional(v, is_unitful(cell), positions, cell)
end
function to_fractional(::Val{:unitful}, ::Val{:unitless}, ::AbstractArray, ::AbstractMatrix)
    error("The cell is unitless. Cannot determine how to make positions fractional.")
end
function to_fractional(::Val{:unitful}, v::Val{:unitful},
                       positions::AbstractArray, cell::AbstractMatrix)
    inv(cell) * positions
end

"""
    to_cartesian(positions::AbstractArray, cell::AbstractMatrix)

Converts positions to Cartesian units (Unitful). If the positions have no units, they are
already Cartesian units, and this operation is a no-op.
"""
function to_cartesian(pos::AbstractArray, cell::AbstractMatrix)
    to_cartesian(is_unitful(pos), pos, cell)
end
to_cartesian(::Val{:unitful}, positions::AbstractArray, ::AbstractMatrix) = positions
function to_cartesian(v::Val{:unitless}, positions::AbstractArray, cell::AbstractMatrix)
    to_cartesian(v, is_unitful(cell), positions, cell)
end
function to_cartesian(::Val{:unitless}, ::Val{:unitless}, ::AbstractArray, ::AbstractMatrix)
    error("The cell has no physical units. " *
          "Cannot determine how to make positions Cartesian.")
end
function to_cartesian(::Val{:unitless}, ::Val{:unitful},
                      positions::AbstractArray, cell::AbstractMatrix)
    cell * positions
end

"""
    to_same_kind(positions::AbstractArray, same::AbstractArray, cell::AbstractMatrix)

Converts positions to same units, Cartesian or fractional, as `same`. If possible, it
returns a reference to `positions` without copy. The function helps maintain _kind_
stability in other utility functions.
"""
function to_same_kind(positions::AbstractArray, same::AbstractArray, cell::AbstractMatrix)
    to_same_kind(is_unitful(same), positions, cell)
end
function to_same_kind(::Val{:unitful}, positions::AbstractArray, cell::AbstractMatrix)
    to_cartesian(positions, cell)
end
function to_same_kind(::Val{:unitless}, positions::AbstractArray, cell::AbstractMatrix)
    to_fractional(positions, cell)
end

"""
    is_periodic(a::AbstractVector, b::AbstractVector, cell::AbstractMatrix;
                tolerance::Real=$(default_tolerance))

True if the positions are one-to-one periodic with respect to the input cell.
"""
function is_periodic(a::AbstractVector, b::AbstractVector, cell::AbstractMatrix;
                     tolerance::Real=default_tolerance)
    const result = abs.(
        mod.(to_fractional(a, cell) .- to_fractional(b, cell) + 0.5, 1) - 0.5
    ) .< tolerance
    all(result)
end

"""
    is_periodic(a::AbstractMatrix, b::AbstractVector, cell::AbstractMatrix;
                tolerance::Real=$(default_tolerance))

Array of boolean describing whether positions in `a` are periodic with positions in `b`.
"""
function is_periodic{T, D, U}(a::AbstractMatrix,
                              b::AbstractVector,
                              cell::AbstractMatrix{Quantity{T, D, U}};
                              tolerance::Real=default_tolerance)
    const result = abs.(
        mod.(to_fractional(a, cell) .- to_fractional(b, cell) + 0.5, 1) - 0.5
    ) .< tolerance
    all(result, 1)
end
function is_periodic{T, D, U}(a::AbstractVector,
                              b::AbstractMatrix,
                              cell::AbstractMatrix{Quantity{T, D, U}};
                              tolerance::Real=default_tolerance)
    is_periodic(b, a, cell, tolerance=tolerance)
end
function is_periodic{T, D, U}(a::AbstractMatrix,
                              b::AbstractMatrix,
                              cell::AbstractMatrix{Quantity{T, D, U}};
                              tolerance::Real=default_tolerance)
    const result = abs.(
        mod.(to_fractional(a, cell) .- to_fractional(b, cell) + 0.5, 1) - 0.5
    ) .< tolerance
    all(result, 1)
end

"""
    $(SIGNATURES)

Folds periodic positions into cell
"""
function into_cell(pos::AbstractArray, cell::AbstractMatrix;
                   tolerance::Real=default_tolerance)
    frac = mod.(to_fractional(pos, cell), 1)
    if tolerance > 0
        for i in 1:length(frac)
            if abs(frac[i] - 1e0) < tolerance
                frac[i] = 0e0
            end
        end
    end
    to_same_kind(frac, pos, cell)
end

"""
    origin_centered(positions::AbstractArrays, cell::AbstractMatrix)

Folds positions back to origin, such that each fractional component ``x_f`` is between
``-0.5\\leq x_f < 0.5``. If the input is in Cartesian (fractional) coordinates, then the
result is also in Cartesian (fractional) coordinates.
"""
function origin_centered(pos::AbstractArray, cell::AbstractMatrix)
    to_same_kind(mod.(to_fractional(pos, cell) .+ 0.5, -1) .+ 0.5, pos, cell)
end

"""
    into_voronoi(positions::AbstractArray, cell::AbstractMatrix; extent::Integer=1)

Folds positions into first Brillouin zone of the input cell. Makes a well-meaning effort at
returning the periodic image with the smallest possible norm. It recenter the atoms around
the origin and then looks for the smallest periodic images within `-extent:extent` cells. If
the cell is quite pathological, then the result will not be within the Voronoi cell.
"""
function into_voronoi(pos::AbstractArray, cell::AbstractMatrix; extent::Integer=1)
    zcentered = to_cartesian(origin_centered(pos, cell), cell)
    result = deepcopy(zcentered)
    norms = [norm(zcentered[:, z]) for z in 1:size(zcentered, 2)]
    for n in eachindex(norms)
        for i = -extent:extent, j = -extent:extent, k = -extent:extent
            translation = cell * [i, j, k]
            position = zcentered[:, n] + translation
            d = norm(position)
            if d < norms[n]
                result[:, n] = position
                norms[n] = d
            end
        end
    end
    to_same_kind(result, pos, cell)
end

"""
    $(SIGNATURES)

Creates a supercell from an input lattice.

# Parameters

* `lattice::Crystal`: the original lattice
* `supercell::AbstractMatrix`: the cell of the supercell in Cartesian coordinates
* `site_id::Bool`: Whether to add/modify an atomic property indicating the index
  of the site in the original lattice
* `cell_id::Bool`: Whether to add/modify an atomic property indicating the index
  of the cell the site belongs to
"""
function supercell(lattice::Crystal, supercell::AbstractMatrix;
                   site_id::Bool=true, tolerance::Real=default_tolerance)
    @argcheck nrow(lattice) != 0

    newcell = to_cartesian(supercell, lattice.cell)
    transform, quotient = hart_forcade(lattice.cell, newcell)
    itransform = inv(transform)

    siteids = collect(take(cycle(1:nrow(lattice)), prod(quotient) * nrow(lattice)))
    result = lattice[siteids]
    result.cell = newcell

    positions = result[:cartesian]
    for (i, index) in zip(1:nrow(lattice):nrow(result), CartesianRange((quotient...)))
        for j in i:i+nrow(lattice)-1
            positions[:, j] += itransform * [index[u] for u in eachindex(quotient)]
        end
    end
    result[:position] = into_cell(positions, result.cell)
    if site_id
        result[:site_id] = siteids
    end
    result
end

""" Tuple holding cell parameters """
const CellParameters = @NT(a, b, c, α, β, γ)

"""
    cell_parameters(a::Quantity, b::Quantity, c::Quantity,
                    α::Quantity=(π/2)u"rad", β::Quantity=(π/2)u"rad",
                    γₒ::Quantity=(π/2)u"rad")

Computes the cell matrix from the cell parameters [a, b, c, α, β, γ].
"""
function cell_parameters(a::Quantity, b::Quantity, c::Quantity,
                         α::Quantity=(π/2)u"rad", β::Quantity=(π/2)u"rad",
                         γₒ::Quantity=(π/2)u"rad")
    cx = cos(β)
    cy = (cos(α) - cos(β)cos(γₒ))/sin(γₒ)
    units = unit(a)
    bb = uconvert(unit(a), b)
    cc = uconvert(unit(a), c)
    z₀ = typeof(a)(0)
    [a    bb * cos(γₒ) cc * cx;
     z₀   bb * sin(γₒ) cc * cy;
     z₀   z₀           cc * √(1 - cx * cx - cy * cy)]
end
cell_parameters(cp::CellParameters) = cell_parameters(cp...)

"""
    cell_parameters(cell::AbstractMatrix)
    cell_parameters(cell::Crystal)

Parameters (a, b, c, α, β, γ) of the input cell returned in a named tuple.
"""
function cell_parameters(cell::AbstractMatrix)
    G = transpose(cell) * cell
    a, b, c = sqrt.(diag(G))
    α = acos(0.5(G[2, 3] + G[3, 2])/(c * b))u"rad"
    β = acos(0.5(G[3, 1] + G[1, 3])/(a * c))u"rad"
    γₒ = acos(0.5(G[1, 2] + G[2, 1])/(a * b))u"rad"
    CellParameters(a, b, c, α, β, γₒ)
end

end


text:
const Disconnected = Dict{String, Set{String}}()

function dcplayerjoined(handle, roomname)
    # No current game
    !haskey(RoomGames, roomname) && return
    # Current game, but no disconnected players
    !haskey(Disconnected, roomname) && return

    if handle in Disconnected[roomname]
        pop!(Disconnected[roomname], handle)
    end
    return
end

function dcplayerleft(handle, roomname)
    # No current game
    !haskey(RoomGames, roomname) && return
    g = RoomGames[roomname]

    # Left user not currently playing
    !(handle in inGameUsers(g)) && return

    if !haskey(Disconnected, roomname)
        Disconnected[roomname] = Set{String}()
    end

    push!(Disconnected[roomname], handle)

    checkEmptyWall(roomname)
    return
end

# Checks whether all players in game disconnected
function checkEmptyWall(roomname)
    # No current game
    !haskey(RoomGames, roomname) && return
    g = RoomGames[roomname]
    # No disconnected users
    !haskey(Disconnected, roomname) && return

    if length(Disconnected[roomname]) == length(inGameUsers(g))
        # All users have disconnected
        removeGameRef(roomname)
    end

    return
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, AuditoryStimuli

makedocs(
    modules = [AuditoryStimuli],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "AuditoryStimuli.jl",
    authors  = "Robert Luke",
    pages = [
        "Home" => "index.md",
        "Introduction / Tutorial" => "realtime-introduction.md",
        "Examples" => Any[
            "Amplitude Modulated Noise" => "example-ssr.md",
            "Harmonic Stacks" => "example-hs.md",
            "Bandpass Noise" => "example-bpnoise.md",
            "Interaural Time Delay" => "example-itd.md",
            "ITD Modulation" => "example-itmfr.md",
            "Signal and Noise" => "example-signoise.md",
        ],
        "Advanced Usage" => "advanced.md",
        "API" => "api.md"
    ]
)

deploydocs(
    repo = "github.com/rob-luke/AuditoryStimuli.jl.git",
    push_preview = true,
    devbranch = "main"
)


text:
function rm_spk_fill_mean(
    x::Neuro1DRealSignal;
    window_width::Int,
    bar::Float64 = 3.0,
    method::Symbol = :sd)::Neuro1DRealSignal

    x_ = convert(Matrix{Float64}, copy(x.signal))
    r, c = size(x_)
    w = ceil(Int64, window_width * x.sample_rate)

    for i = 1:(r - w)
        window_x = x_[i:(i+w-1), :]
        window_μ = mean(window_x, dims = 1)
        for ch = 1:c
            ch_outlier = outliers(window_x[:, ch], method = method, bar = bar)
            window_x[ch_outlier, ch] .= window_μ[ch]
        end

        x_[i:(i+w-1), :] = window_x
    end

    Neuro1DRealSignal(
        signal = x_,
        sample_rate = x.sample_rate
    )
end

function rm_spk(
    x::Neuro1DRealSignal,
    method = :fill_mean;
    kwargs...)::Neuro1DRealSignal

    if method == :fill_mean
        rm_spk_fill_mean(x; kwargs...)
    else
        throw(ArgumentError("invalid method"))
    end
end


text:
# Map of Soapy stream formats to Julia types

"""
Abstract type denoting a Complex(U)Int(12/4) type.

These need to be specially handled with bit shifting, and
we do no data processing in this library, so this
exists for convience to other package developers.

Subtypes are `ComplexInt` and `ComplexUInt`,
for convience with handling sign extends.
"""
abstract type AbstractComplexInteger end

"""
Type indicating a Complex Int format.

Parameter `T` indicates the number of bits.
"""
struct ComplexInt{T} <: AbstractComplexInteger; end

"""
Type indicating a Complex Unsigned Int format.

Parameter `T` indicates the number of bits.
"""
struct ComplexUInt{T} <: AbstractComplexInteger; end

const _stream_type_pairs = [
    (SOAPY_SDR_CF64, Complex{Float64}),
    (SOAPY_SDR_CF32, Complex{Float32}),
    (SOAPY_SDR_CS32, Complex{Int32}),
    (SOAPY_SDR_CU32, Complex{UInt32}),
    (SOAPY_SDR_CS16, Complex{Int16}),
    (SOAPY_SDR_CU16, Complex{UInt16}),
    (SOAPY_SDR_CS12, ComplexInt{12}),
    (SOAPY_SDR_CU12, ComplexUInt{12}),
    (SOAPY_SDR_CS8, Complex{Int8}),
    (SOAPY_SDR_CU8, Complex{UInt8}),
    (SOAPY_SDR_CS4, ComplexInt{4}),
    (SOAPY_SDR_CU4, ComplexUInt{4}),
    (SOAPY_SDR_F64, Float64),
    (SOAPY_SDR_F32, Float32),
    (SOAPY_SDR_S32, Int32),
    (SOAPY_SDR_U32, UInt32),
    (SOAPY_SDR_S16, Int16),
    (SOAPY_SDR_U16, UInt16),
    (SOAPY_SDR_S8, Int8),
    (SOAPY_SDR_U8, UInt8) ]

"""
Type map from SoapySDR Stream formats to Julia types.

Note: Please see ComplexUInt and ComplexUInt if using 12 or 4 bit complex types.
"""
const _stream_type_soapy2jl = Dict{String, Type}(_stream_type_pairs)

"""
Type map from SoapySDR Stream formats to Julia types.

Note: Please see ComplexUInt and ComplexUInt if using 12 or 4 bit complex types.
"""
const _stream_type_jl2soapy = Dict{Type, String}(reverse.(_stream_type_pairs))

function _stream_map_jl2soapy(stream_type)
    if !haskey(_stream_type_jl2soapy, stream_type)
        error("Unsupported stream type: " + stream_type)
    end
    return _stream_type_jl2soapy[stream_type]
end

function _stream_map_soapy2jl(stream_type)
    if !haskey(_stream_type_soapy2jl, stream_type)
        error("Unsupported stream type: " + stream_type)
    end
    return _stream_type_soapy2jl[stream_type]
end


text:
using NamedDims
using NamedDims:
    names,
    unify_names,
    unify_names_longest,
    order_named_inds,
    permute_dimnames,
    remaining_dimnames_from_indexing,
    remaining_dimnames_after_dropping
using Test

@testset "dim" begin
    @testset "small case" begin
        @test dim((:x, :y), :x) == 1
        @test dim((:x, :y), :y) == 2
        @test_throws ArgumentError dim((:x, :y), :z)  # not found
    end
    @testset "large case" begin
        @test dim((:x, :y, :a, :b, :c, :d), :x) == 1
        @test dim((:x, :y, :a, :b, :c, :d), :a) == 3
        @test dim((:x, :y, :a, :b, :c, :d), :d) == 6
        @test_throws ArgumentError dim((:x, :y, :a, :b, :c, :d), :z) # not found
    end
end

 @testset "unify_names/unify_names_longest" begin
    @test_throws DimensionMismatch unify_names((:a,), (:a, :b,))

    @test unify_names_longest((:a,), (:a, :b,)) == (:a, :b)
    @test unify_names_longest((:a,), (:a, :_)) == (:a, :_)
    @test unify_names_longest((:a, :b), (:a, :_, :c)) == (:a, :b, :c)

    @test_throws DimensionMismatch unify_names_longest((:a, :b, :c), (:b, :a))

    for unify in (unify_names, unify_names_longest)
        @test unify((:a,), (:a,)) == (:a,)
        @test unify((:a, :b), (:a, :b)) == (:a, :b)
        @test unify((:a, :_), (:a, :b)) == (:a, :b)
        @test unify((:a, :_), (:a, :_)) == (:a, :_)

        @test unify((:a, :b, :c), (:_, :_, :_)) == (:a, :b, :c)
        @test unify((:a, :_, :c), (:_, :b, :_)) == (:a, :b, :c)
        @test unify((:_, :_, :_), (:_, :_, :_)) == (:_, :_, :_)

        @test_throws DimensionMismatch unify((:a,), (:b,))
        @test_throws DimensionMismatch unify((:a,:b), (:b, :a))
        @test_throws DimensionMismatch unify((:a, :b, :c), (:_, :_, :d))
    end
end

@testset "order_named_inds" begin
    @test order_named_inds((:x,)) == (:,)
    @test order_named_inds((:x,); x=2) == (2,)

    @test order_named_inds((:x, :y,)) == (:, :)
    @test order_named_inds((:x, :y); x=2) == (2, :)
    @test order_named_inds((:x, :y); y=2, ) == (:, 2)
    @test order_named_inds((:x, :y); y=20, x=30) == (30, 20)
    @test order_named_inds((:x, :y); x=30, y=20) == (30, 20)
end

@testset "remaining_dimnames_from_indexing" begin
    @test remaining_dimnames_from_indexing((:a, :b, :c), (10,20,30)) == tuple()
    @test remaining_dimnames_from_indexing((:a, :b, :c), (10, :,30)) == (:b,)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (1:1, [true], [20])) == (:a, :b, :c)
end

@testset "remaining_dimnames_after_dropping" begin
    @test remaining_dimnames_after_dropping((:a, :b, :c), 1) == (:b, :c)
    @test remaining_dimnames_after_dropping((:a, :b, :c), 3) == (:a, :b)
    @test remaining_dimnames_after_dropping((:a, :b, :c), (1,3)) == (:b,)
    @test remaining_dimnames_after_dropping((:a, :b, :c), (3,1)) == (:b,)
end

@testset "permute_dimnames" begin
    @test permute_dimnames((:a, :b, :c), (1, 2, 3)) == (:a ,:b, :c)
    @test permute_dimnames((:a, :b, :c), (3, 2, 1)) == (:c ,:b, :a)

    # permute_dimnames allows non-bijective "permutations"
    @test permute_dimnames((:a, :b, :c), (3, 3, 1)) == (:c ,:c, :a)
    @test permute_dimnames((:a, :b, :c), (3, 1))== (:c, :a)

    @test_throws BoundsError permute_dimnames((:a, :b, :c), (30, 30, 30))
    @test_throws BoundsError permute_dimnames((:a, :b), (1, 0))
end


text:
module Tarana

abstract type AbstractPitchClass end

struct PitchClass <: AbstractPitchClass
    s::Symbol
end

const Pitch{T <: AbstractPitchClass} = Tuple{T, Integer}

struct Note{T <: Number}
    t::T
    p::Pitch
end

struct Rest{T <: Number}
    t::T
end

pc = PitchClass

const
    Cff = pc(:Cff)
    Cf  = pc(:Cf)
    C   = pc(:C)
    Cs  = pc(:Cs)
    Css = pc(:Css)
    Dff = pc(:Dff)
    Df  = pc(:Df)
    D   = pc(:D)
    Ds  = pc(:Ds)
    Dss = pc(:Dss)
    Eff = pc(:Eff)
    Ef  = pc(:Ef)
    E   = pc(:E)
    Es  = pc(:Es)
    Ess = pc(:Ess)
    Fff = pc(:Fff)
    Ff  = pc(:Ff)
    F   = pc(:F)
    Fs  = pc(:Fs)
    Fss = pc(:Fss)
    Gff = pc(:Gff)
    Gf  = pc(:Gf)
    G   = pc(:G)
    Gs  = pc(:Gs)
    Gss = pc(:Gss)
    Aff = pc(:Aff)
    Af  = pc(:Af)
    A   = pc(:A)
    As  = pc(:As)
    Ass = pc(:Ass)
    Bff = pc(:Bff)
    Bf  = pc(:Bf)
    B   = pc(:B)
    Bs  = pc(:Bs)
    Bss = pc(:Bss)

cff(o, d) = Note(d, (Cff, o))
cf(o, d)  = Note(d, (Cf, o))
c(o, d)   = Note(d, (C, o))
cs(o, d)  = Note(d, (Cs, o))
css(o, d) = Note(d, (Css, o))
dff(o, d) = Note(d, (Dff, o))
df(o, d)  = Note(d, (Df, o))
d(o, d)   = Note(d, (D, o))
ds(o, d)  = Note(d, (Ds, o))
dss(o, d) = Note(d, (Dss, o))
eff(o, d) = Note(d, (Eff, o))
ef(o, d)  = Note(d, (Ef, o))
e(o, d)   = Note(d, (E, o))
es(o, d)  = Note(d, (Es, o))
ess(o, d) = Note(d, (Ess, o))
fff(o, d) = Note(d, (Fff, o))
ff(o, d)  = Note(d, (Ff, o))
f(o, d)   = Note(d, (F, o))
fs(o, d)  = Note(d, (Fs, o))
fss(o, d) = Note(d, (Fss, o))
gff(o, d) = Note(d, (Gff, o))
gf(o, d)  = Note(d, (Gf, o))
g(o, d)   = Note(d, (G, o))
gs(o, d)  = Note(d, (Gs, o))
gss(o, d) = Note(d, (Gss, o))
aff(o, d) = Note(d, (Aff, o))
af(o, d)  = Note(d, (Af, o))
a(o, d)   = Note(d, (A, o))
as(o, d)  = Note(d, (As, o))
ass(o, d) = Note(d, (Ass, o))
bff(o, d) = Note(d, (Bff, o))
bf(o, d)  = Note(d, (Bf, o))
b(o, d)   = Note(d, (B, o))
bs(o, d)  = Note(d, (Bs, o))
bss(o, d) = Note(d, (Bss, o))

export
    Note,
    Rest

export
    Cff, Cf, C, Cs, Css,
    Dff, Df, D, Ds, Dss,
    Eff, Ef, E, Es, Ess,
    Fff, Ff, F, Fs, Fss,
    Gff, Gf, G, Gs, Gss,
    Aff, Af, A, As, Ass,
    Bff, Bf, B, Bs, Bss

export
    cff, cf, c, cs, css,
    dff, df, d, ds, dss,
    eff, ef, e, es, ess,
    fff, ff, f, fs, fss,
    gff, gf, g, gs, gss,
    aff, af, a, as, ass,
    bff, bf, b, bs, bss

const pc_to_i = Dict(
    Cff => -2, Cf => -1, C =>  0, Cs =>  1, Css =>  2,
    Dff =>  0, Df =>  1, D =>  2, Ds =>  3, Dss =>  4,
    Eff =>  2, Ef =>  3, E =>  4, Es =>  5, Ess =>  6,
    Fff =>  3, Ff =>  4, F =>  5, Fs =>  6, Fss =>  7,
    Gff =>  5, Gf =>  6, G =>  7, Gs =>  8, Gss =>  9,
    Aff =>  7, Af =>  8, A =>  9, As => 10, Ass => 11,
    Bff =>  9, Bf => 10, B => 11, Bs => 12, Bss => 13)

const pc_to_str = Dict(
    Cff => "C♭♭", Cf => "C♭", C => "C", Cs => "C♯", Css => "C♯♯",
    Dff => "D♭♭", Df => "D♭", D => "D", Ds => "D♯", Dss => "D♯♯",
    Eff => "E♭♭", Ef => "E♭", E => "E", Es => "E♯", Ess => "E♯♯",
    Fff => "F♭♭", Ff => "F♭", F => "F", Fs => "F♯", Fss => "F♯♯",
    Gff => "G♭♭", Gf => "G♭", G => "G", Gs => "G♯", Gss => "G♯♯",
    Aff => "A♭♭", Af => "A♭", A => "A", As => "A♯", Ass => "D♯♯",
    Bff => "B♭♭", Bf => "B♭", B => "B", Bs => "B♯", Bss => "B♯♯")

include("./unicode.jl")

const eqtemp = [C, C♯, D, D♯, E, F, F♯, G, G♯, A, A♯, B]

ord(p) = ord(p...)
ord(n::Note) = ord(n.p)
ord(pc, oct) = 12 * (oct + 1) + pc_to_i[pc]

pitch(ap::Integer) =
    let (oct, n) = divrem(ap, 12)
        (eqtemp[n + 1], oct - 1)
    end

trans(i, p::Pitch) = ord(p) + i |> pitch
trans(i, n::Note) = Note(n.t, trans(i, n.p))

export
    ord,
    pitch,
    trans

# Base.show(io::IO, x::PitchClass) = print(io, pc_to_str[x])

Base.:+(a::Pitch, b) = trans(b, a)
Base.:-(a::Pitch, b) = trans(-b, a)

Base.:+(a::Note, b) = Note(a.t, a.p + b)
Base.:-(a::Note, b) = Note(a.t, a.p - b)

Base.isless(a::Pitch, b::Pitch) = ord(a) < ord(b)
Base.isless(a::PitchClass, b::PitchClass) = pc_to_i[a] < pc_to_i[b]

end # module

max_stars_repo_name: basp/Tarana

max_issues_repo_name: basp/Tarana

max_forks_repo_name: basp/Tarana


text:
struct TransCovariate{S<:AbstractFloat, T<:AbstractFloat} <: AbstractCovariate{T}
    name::String
    basecovariate::AbstractCovariate{S}
    transform::Function
end

Base.length(var::TransCovariate) = length(var.basecovariate)

function TransCovariate(name::String, basecovariate::AbstractCovariate{T}, transform::Function) where {T<:AbstractFloat}
    TransCovariate{T, T}(name, basecovariate, transform)
end

function TransCovariate(::Type{T}, name::String, basecovariate::AbstractCovariate{S}, transform::Function) where {T<:AbstractFloat} where {S<:AbstractFloat}
    TransCovariate{S, T}(name, basecovariate, transform)
end

function slice(covariate::TransCovariate{S, T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat} where {S<:AbstractFloat}
    basecov = covariate.basecovariate
    f = covariate.transform
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(basecov, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, T)
end

function Base.map(covariate::TransCovariate, dataframe::AbstractDataFrame)
    TransCovariate(covariate.name, map(covariate.basecovariate, dataframe), covariate.transform)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
"""

"""
function products()
  res = HTTP.get(ENDPOINT * "products")
  rawData = JSON.parse(String(res.body))
  vcat(DataFrame.(rawData)...)
end

function product(id::String)
  res = HTTP.get(ENDPOINT * "products/" * id)
  rawData = JSON.parse(String(res.body))
  DataFrame(rawData)
end

function product(ccy1::String, ccy2::String)
  prodcut(ccy1 * "-" * ccy2)
end


text:
# Utility methods for aggregating Select Registers together based on the limited number
# of counters available.
#
# This will allow a suite of counter values to be setup and will automatically keep
# iterating until all counters have been configured.
#
# If it becomes relevant, this will have to be refactored to be aware if we're targeting
# IMC/CHA etc because the rules vary slightly.

numconcurrent(x::T) where {T} = numconcurrent(T)
function numconcurrent(::Type{T}) where {T}
    error("Number of Concurrent Counters not defined for type $T")
end

numconcurrent(::Type{CoreSelectRegister}) = numcounters()

# Uncore PMUs fixed at 4 concurrent events.
# This logic gets more complicated when we throw on CHA Filters and the restriction that
# only a subset of the PMU counters can take some events - but for now we won't worry
# about that.
numconcurrent(::Type{UncoreSelectRegister}) = 4

"""
    group(x) -> itr

Return an iterator `itr` that partitions a collection of Selection Registers in a way
that matches the underlying capability of the hardware.
"""
group(x) = Iterators.partition(x, numconcurrent(eltype(x)))


text:
module Tester

using Genie, Genie.Util, Genie.Configuration, Genie.Loggers

"""
    bootstrap_tests(cmd_args::String, config::Settings) :: Nothing

Sets up testing environment, includes test files, etc.
"""
function bootstrap_tests(cmd_args::String = "", config::Settings = Genie.config, resource::String = "") :: Nothing
  current_env = config.app_env

  set_test_env()

  include(abspath(joinpath(config.test_folder, "test_config.jl")))

  for (path, _, files) in (walkdir(abspath(joinpath(config.test_folder))) |> collect)
    for file_name in files
      isempty(resource) && endswith(file_name, "_test.jl") && include(joinpath(path, file_name))
      ! isempty(resource) && startswith(file_name, resource) && endswith(file_name, "_test.jl") && include(joinpath(path, file_name))
    end
  end

  Genie.config.app_env = current_env
  log("Switched app to >> $(uppercase(Genie.config.app_env)) << env", :debug)

  nothing
end

"""
    reset_db() :: Nothing

Prepares the test env DB running all migrations up.
"""
function reset_db() :: Nothing
  Migration.all_down()
  Migration.all_up()

  nothing
end

"""
    run_all_tests(cmd_args::String, config::Settings) :: Nothing

Runs all existing tests.
"""
function run_all_tests(cmd_args::String = "", config::Settings = Genie.config) :: Nothing
  bootstrap_tests(cmd_args, config)

  nothing
end

function run_tests() :: Nothing
  run_all_tests()
end
function run_tests(resource_name::Symbol) :: Nothing
  bootstrap_tests("", Genie.config, string(resource_name) |> lowercase)
end

"""
    set_test_env() :: Nothing

Switches Genie to the test env for the duration of the current execution.
"""
function set_test_env() :: Nothing
  if Genie.config.tests_force_test_env
    Genie.config.app_env = TEST
    log("Switched app to >> $(uppercase(App.config.app_env)) << env", :debug)

    ! is_test() && error("Could not switch env")
  end

  nothing
end

end


text:
using IcanProj
using NetCDF
using DataFrames
using CSV
using Statistics

function compute_statistics(var_coarse, var_ref, mask)

    n = size(var_coarse, 2)

    nse = fill(0.0, n)
    r = fill(0.0, n)
    bias = fill(0.0, n)
    rmse = fill(0.0, n)

    for i in 1:n

        ikeep = findall(mask[:, i])

        if length(ikeep) > 10

            sim = var_coarse[ikeep, i]

            ref = var_ref[ikeep, i]

            nse[i] = 1 - var(sim .- ref) / var(ref .- mean(ref))

            r[i] = cor(sim, ref)

            bias[i] = mean(sim .- ref)

            rmse[i] = sqrt(mean((sim .- ref).^2))

        else

            nse[i] = NaN

            r[i] = NaN

            bias[i] = NaN

            rmse[i] = NaN

        end

    end

    return nse, r, bias, rmse

end

function results_table(path, cfgs, variables, spaceres)

    # Create empty table

    file_tmp = joinpath(path, "results_1", "swe_$(spaceres).nc")

    idname = Symbol(ncgetatt(file_tmp, "id", "id"))

    ids = convert(Array{Int64}, ncread(file_tmp, "id"))

    df_res = DataFrame()

    df_res[idname] = ids

    for c in cfgs

        println("Running cfg $(c)")

        # Mask indicating presence of snow

        file_fine = joinpath(path, "results_$(c)", "swe_1km.nc")

        file_coarse = joinpath(path, "results_$(c)", "swe_$(spaceres).nc")

        df_links = link_results(file_fine, file_coarse)

        swe_coarse, swe_ref = unify_results(file_fine, file_coarse, df_links, "swe")

        mask = (swe_coarse .> 0) .| (swe_ref .> 0)

        # Compute statistics for all variables

        for v in variables

            file_fine = joinpath(path, "results_$(c)", "$(v)_1km.nc")

            file_coarse = joinpath(path, "results_$(c)", "$(v)_$(spaceres).nc")

            var_coarse, var_ref, ngrids = unify_results(file_fine, file_coarse, df_links, v)

            nse, r, bias, rmse = compute_statistics(var_coarse, var_ref, mask)

            nse_name = Symbol("nse_$(v)_cfg$(c)")

            r_name = Symbol("r_$(v)_cfg$(c)")

            bias_name = Symbol("bias_$(v)_cfg$(c)")

            rmse_name = Symbol("rmse_$(v)_cfg$(c)")

            df_res[nse_name] = nse

            df_res[r_name] = r

            df_res[bias_name] = bias

            df_res[rmse_name] = rmse

            df_res[:ngrids] = ngrids

        end

    end

    return df_res

end

path = "/data04/jmg/fsm_simulations/netcdf/fsmres_forest"

cfgs = 1:32

variables = ["gsurf", "hatmo", "latmo", "melt", "rnet", "rof", "snowdepth", "swe"]

for spaceres in ["5km", "10km", "25km", "50km"]

    println("Running spaceres $(spaceres)")

    df = results_table(path, cfgs, variables, spaceres)

    filesave = joinpath(dirname(pathof(IcanProj)), "..", "data", "table_errors_$(spaceres).txt")

    df |> CSV.write(filesave)

end


text:
stringvec(a,b) = "$a[$b]"
# Specialized "show" method for slightly easier viewing.
function Base.show(io::IO, x::AsapInstruction)
    # Specializations on OP Type
    if x.op == :RPT
        show_rpt(io, x)
        return
    elseif x.optype == BRANCH_TYPE
        show_branch(io, x)
        return
    end

    show_basic(io, x)
    show_options(io, x)
    show_extra(io, x)
end

function show_basic(io::IO, x::AsapInstruction)
    # Show the op
    print(io, x.op, " ")

    def = getdef(x.op)
    # Print destinations and sources
    def.dest && print(io, "$(x.dest) ")
    def.src1 && print(io, "$(x.src1) ")
    def.src2 && print(io, "$(x.src2) ")

    return nothing
end

function show_options(io::IO, x::AsapInstruction)
    # Print out instruction options
    if x.nops > 0
        print(io, "nops=$(x.nops) ")
    end
    if x.jump
        print(io, "jump ")
    end
    if x.cxflag != NO_CX
        print(io, "$(stringvec(x.cxflag, x.cxindex)) ")
    end
    if !x.sw
        print(io, "DW ")
    end
end

function show_extra(io::IO, x::AsapInstruction)
    # Extra stuff
    if get(io, :compact, false)
        print(io, x.optype, " ")
        if x.signed
            print(io, "signed ")
        end
        if x.saturate
            print(io, "saturate ")
        end
        if x.dest_is_output
            print(io, "output")
        end
    end
end

function show_rpt(io::IO, x::AsapInstruction)
    # Show the OP
    print(io, x.op, " ")
    print(io, "start=$(repeat_start(x)) stop=$(repeat_end(x)) ")
    print(io, "repeats=$(x.src1) ")

    show_options(io, x)
    show_extra(io, x)
end

function show_branch(io::IO, x::AsapInstruction)
    # Show the op
    print(io, x.op, " ")

    # Print destinations and sources
    print(io, "$(x.dest) ")
    print(io, "$(x.src1) ")
    show_options(io, x)
    return nothing
end


text:
try
  if Sys.iswindows()
    println("cmd /c deltree SampleFiles/")
    run(`cmd /c deltree SampleFiles/`)
  else
    println("rm -rf SampleFiles/")
    run(`rm -rf SampleFiles/`)
  end
catch err
  println("Attempting to delete files failed! Error thrown: ", err)
end


text:
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils
using Profile
using BenchmarkTools
using FFTW
using HDF5
using Test
using SparseArrays

filter_hash = fink_filter_hash(1, 8, nx=256, pc=1, wd=1)

@time filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=256)

filter_temp = zeros(256,256,256)
filter_temp[filt_3d["filt_index"][56]] = filt_3d["filt_value"][56]

h5write("../DHC/scratch_AKS/data/filt_3d_rs.h5", "main/data", fftshift(real(fft(filter_temp))))
h5write("../DHC/scratch_AKS/data/filt_3d_fs.h5", "main/data", fftshift(filter_temp))

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
"""
```
SeisRadonFreqInv(d; <keyword arguments>)
```

Transform a Gather from time-offset gather to tau-p gather using a frequency
domain inverse parabolic or linear Radon operator via least-squares inversion.

# Arguments
* `d::Array{T<:Real,2}`: 2D data, `d[1:nt,1:nh]`, where `nt` is number of
time samples and `nh` the number of receivers.

# Keyword arguments
* `order::AbstractString="parab"`: `"parab"` for parabolic transform, `"linear"`
for linear transform.
* `dt::Real=0.004`: sampling interval in seconds.
* `h::Vector{Real}=collect(0.0:20.0:1000.0)`: offset vector `h[1:nh]`.
* `href::Real=0.0`: reference offset for parabolic Radon Transform. If the
defautl value `href=0.0` is given, `href` is set to `max(abs(h))`.
* `p::Vector{Real}=collect(-0.05:0.01:2.2)`: `p[1:np]`. If `order="parab"`, `p`
is a vector of residual moveout ("curvatures") at reference offset `href` in
seconds; if `order=linear`, `p` is a vector of ray parameters in s/m.
* `flow::Real=0.0`: minimum frequency in the data in Hz.
* `fhigh::Real=125.0`: maximum frequency in the data in Hz.
* `mu::Real=0.001`: trade off parameter or damping for the L.S. inversion.

# Output
* `m`: inverted Radon panel `m[1:ntau, 1:np]`.

# References
*  Hampson, D., 1986, Inverse velocity stacking for multiple elimination:
Canadian Journal of Exploration Geophysics, 22, 44-55.
*  Sacchi, M.D. and Ulrych, T.J., 1995, High-resolution velocity gathers and
offset space reconstruction: Geophysics, 60, 1169-1177.
"""

function SeisRadonFreqInv{Td<:Real, Th<:Real, Tp<:Real
                          }(d::Array{Td,2}; order::AbstractString="parab",
                            dt::Real=0.004, href::Real=0.0,
                            h::Vector{Th}=collect(0.0:20.0:1000.0),
                            p::Vector{Tp}=collect(-0.05:0.01:2.2),
                            flow::Real=0.0, fhigh::Real=125.0, mu::Real=0.001)
    if order=="parab"
        I = 2
        href == 0 && (href = maximum(abs.(h)))
    elseif order=="linear"
        I = 1
        href = 1.0
    else
        error("Order should be equal to \"parab\" or \"linear\"")
    end
    nt = size(d, 1)
    nh = length(h)
    nh == size(d, 2) || error("lenght(h) must be equal to size(d, 2)")
    np = length(p)
    nw = 2*nextpow2(nt)
    d = cat(1, d, zeros(Td, nw-nt, nh))
    D = fft(d, 1)
    iw_low = round(Int, flow*dt*nw+1)
    iw_high = round(Int, fhigh*dt*nw+1)
    M = zeros(Complex{Td}, nw, np)
    for iw = iw_low:iw_high
        w  = 2.0*pi*(iw-1)/(nw*dt)
        L = zeros(Complex{Td}, nh, np)
        for ip = 1:np
            for ih = 1:nh
                phi = w*p[ip]*(h[ih]/href)^I
      	        L[ih, ip] = exp(-im*phi)
            end
        end
        R = L'*L
        r0 = trace(R)/np
        xa = L'*D[iw, :]
        Q = mu*r0*eye(np)
        x  = (R + Q)\xa
        M[iw, :] = x.'
    end
    for iw = round(Int, nw/2)+2:nw
        M[iw, :] = conj(M[nw-iw+2, :])
    end
    m = real(ifft(M, 1))
    m = m[1:nt, :]
    return m
end


text:
using ChebyshevQuantum
using Test

tol_val = 1e-11

@testset "test interp" begin

    function f(x)
        return exp(-0.2*sin(3*x))
    end

    c = ChebyshevQuantum.Interp.make_cheb(f , a=2, b = 3)

    @test isapprox(c(2.3), f(2.3), atol=tol_val)
    @test isapprox(c(2.0), f(2.0), atol=tol_val)
    @test isapprox(c(2.712), f(2.712), atol=tol_val)

    function f2(x)
        return log(x) + x^2
    end

    c2 = ChebyshevQuantum.Interp.make_cheb(f2 , a=2, b = 3)

    @test isapprox(c2(2.3), f2(2.3), atol=tol_val)
    @test isapprox(c2(2.0), f2(2.0), atol=tol_val)
    @test isapprox(c2(2.712), f2(2.712), atol=tol_val)

end

@testset "test functions var N" begin

    function f1(x)
        return exp(-0.2*sin(3*x))
    end
    function f2(x)
        return exp(-0.3*cos(3*x)) + 20
    end

    c1 = ChebyshevQuantum.Interp.make_cheb(f1 , a=2, b = 3)
    c2 = ChebyshevQuantum.Interp.make_cheb(f2 , a=2, b = 3)

    x=2.712
    @test isapprox( (c1+c2)(x), f1(x)+f2(x) , atol=tol_val)
    @test isapprox( (c1-c2)(x), f1(x)-f2(x) , atol=tol_val)
    @test isapprox( (c1^2)(x), f1(x)^2 , atol=tol_val)
    @test isapprox( (c1/c2)(x), f1(x)/f2(x) , atol=tol_val)

    @test isapprox( (-c1)(x), -f1(x) , atol=tol_val)
    @test isapprox( (2*c1)(x), 2*f1(x) , atol=tol_val)
    @test isapprox( (c1*2)(x), 2*f1(x) , atol=tol_val)
    @test isapprox( (2+c1)(x), 2+f1(x) , atol=tol_val)
    @test isapprox( (c1+2)(x), f1(x)+2 , atol=tol_val)

    @test isapprox( (2 / c2)(x), 2 / f2(x) , atol=tol_val)
    @test isapprox( (2 - c2)(x), 2 - f2(x) , atol=tol_val)

end

@testset "test functions fixed N" begin

    function f1(x)
        return exp(-0.2*sin(3*x))
    end
    function f2(x)
        return exp(-0.3*cos(3.5*x)) + 20
    end

    c1 = ChebyshevQuantum.Interp.make_cheb(f1 , a=2, b = 3, N=20, thr=-1)
    c2 = ChebyshevQuantum.Interp.make_cheb(f2 , a=2, b = 3, N=20, thr=-1)

    x=2.712
    @test isapprox( (c1+c2)(x), f1(x)+f2(x) , atol=tol_val)
    @test isapprox( (c1-c2)(x), f1(x)-f2(x) , atol=tol_val)
    @test isapprox( (c1^2)(x), f1(x)^2 , atol=tol_val)
    @test isapprox( (c1/c2)(x), f1(x)/f2(x) , atol=tol_val)

    @test isapprox( (-c1)(x), -f1(x) , atol=tol_val)
    @test isapprox( (2*c1)(x), 2*f1(x) , atol=tol_val)
    @test isapprox( (c1*2)(x), 2*f1(x) , atol=tol_val)
    @test isapprox( (2+c1)(x), 2+f1(x) , atol=tol_val)
    @test isapprox( (c1+2)(x), f1(x)+2 , atol=tol_val)

    @test isapprox( (2 / c2)(x), 2 / f2(x) , atol=tol_val)
    @test isapprox( (2 - c2)(x), 2 - f2(x) , atol=tol_val)

end

@testset "function returns complex number " begin

    function f(x)
        return exp(-1.0im*x)
    end

    c = ChebyshevQuantum.Interp.make_cheb(f , a=-pi, b = pi)

    x= pi/8
    @test isapprox( real(c(x)), real(f(x)) , atol=tol_val)
    @test isapprox( imag(c(x)), imag(f(x)) , atol=tol_val)

    x= 0.1234
    @test isapprox( real(c(x)), real(f(x)) , atol=tol_val)
    @test isapprox( imag(c(x)), imag(f(x)) , atol=tol_val)

    x= 0.5678
    @test isapprox( real(c)(x), real(f(x)) , atol=tol_val)
    @test isapprox( imag(c)(x), imag(f(x)) , atol=tol_val)

end


text:
using ManifoldMeasures
using Test

@testset "ManifoldMeasures.jl" begin
    include("hausdorff.jl")
    include("vonmisesfisher.jl")
end