text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 49a44e94-29c0-4004-8d0b-e37008eafc1c
using BenchmarkTools, PlutoUI, Wordlegames

# ╔═╡ f9816eb0-078e-41f2-bf7c-29ed5c663fd0
# hideall
title = "Scoring guesses in Wordle";

# ╔═╡ a9678ae8-48c0-4bb9-920d-7de176105bd0
"""
+++
title = "$title"
+++
""" |> Base.Text

# ╔═╡ 4eabf46e-f75c-4b15-8583-6abe17c0fd85
md"""
# $title

[Wordle](https://en.wikipedia.org/wiki/Wordle) is a recently developed, extremely popular word game that has already spawned many imitators such as [Primel](https://converged.yt/primel/).

These tutorials illustrate some [Julia](https://julialang.org) programming concepts using functions in the [Wordlegames](https://github.com/dmbates/Wordlegames.jl) package for illustration.
Part of the purpose is to illustrate the unique nature of Julia as a dynamically-typed language with a just-in-time (JIT) compiler.
It allows you to write "generic", both in the common meaning of "general purpose" and in the technical meaning of generic functions, and performative code.

This posting originated from a conversation on the Julia [discourse channel](https://discourse.julialang.org/t/rust-julia-comparison-post/75403) referring to a case where Julia code to perform a certain Wordle-related task - determine the "best" initial guess in a Wordle game - was horribly slow.
Julia code described in a [Hacker News](https://news.ycombinator.com/) posting took several hours to do this.

In situations like this the Julia community inevitably responds with suggested modifications to make the code run faster.
Someone joked that we wouldn't be satisfied until we could do that task in less than 1 second, and we did.

The code in these postings can be used to solve a Wordle game very rapidly, as well as related games like Primel.

Before beginning we attach some packages that will be used in this notebook.
"""

# ╔═╡ 03b41118-5756-447b-9900-15af37667529
md"""
## Target pools

If you are not familiar with the rules of Wordle, please check the [Wikipedia page](https://en.wikipedia.org/wiki/Wordle).
It is a word game with the objective of guessing a 5-letter English word, which we will call the "target".
The target word is changed every day but it is always chosen from a set of 2315 words, which we will call the "target pool".

The original target pool is available with the `Wordlegames` package.  (Apparently the New York Times removed a few of these words after they purchased the rights to Wordle.)
"""

# ╔═╡ 73e3e38b-0935-4ef9-97b1-8a8cdfa2feb7
datadir = joinpath(pkgdir(Wordlegames), "data");

# ╔═╡ 4c5e781c-fd10-425f-95d2-1fec5bf1015f
wordlestrings = collect(readlines(joinpath(datadir, "Wordletargets.txt")))

# ╔═╡ 2046a4ca-c4d3-4804-9373-930d9cbb58fd
md"We call this pool `wordlestrings` because it is stored as a vector of `String`s"

# ╔═╡ 560ad4d4-945c-4d21-8947-cac68c1880a1
typeof(wordlestrings)

# ╔═╡ ddd7f90c-7cad-4279-a275-b1800b4ae83b
md"""
The `Wordlegames` package defines a `GamePool` struct for playing Wordle or related games.
In that `struct` the `String`s are converted to a more efficient storage mode as a vector of `NTuple{5,Char}`, which takes advantage of the fact that each string is exactly 5 characters long.

Speaking of which, it would be a good idea to check that this collection has the properties we were told it had.
It should be a vector of 2315 strings, each of which is 5 characters.
"""

# ╔═╡ 49500976-7ab8-43cc-9902-b25ad73ff42e
length(wordlestrings)

# ╔═╡ f6681018-f729-48de-9342-a5bcb151aca2
all(w -> length(w) == 5, wordlestrings)

# ╔═╡ 22f7e1ad-3a2b-4a88-8f81-91feada88621
md"""
That last expression may look, well, "interesting".
It is a way of checking that a function, in this case an anonymous function expressed using the [stabby lambda](https://dev.to/keithrbennett/why-i-prefer-stabby-lambda-notation-5gcj) notation, returns `true` for each element of an iterator, in this case the vector `wordlestrings`.
You can read the whole expression as "is `length(w)` equal to `5` for each word `w` in `wordlestrings`".

These words are supposed to be exactly 5 letters long but it never hurts to check.
I've been a data scientist for several decades and one of the first lessons in the field is to [trust, but verify](https://en.wikipedia.org/wiki/Trust%2C_but_verify) any claims about the data you are provided.

It turns out this check is redundant because the property is checked when creating a `GamePool`.
"""

# ╔═╡ 5787eef2-7a2c-46ad-8fa0-7313c07be84b
wordle = GamePool(wordlestrings);

# ╔═╡ 8df0e0c0-cdae-450c-abc1-5cea1d520ec9
propertynames(wordle)

# ╔═╡ 8580bdcf-24ce-4eb6-ac5a-4c77b729b38b
typeof(wordle.guesspool)

# ╔═╡ 9ed9e8ab-3ebc-4c6e-b2c3-46ad1987561e
first(wordle.guesspool, 3)

# ╔═╡ fdd6631e-8b3a-4b4f-b8a7-11135f13cc23
md"""
## Game play

A Wordle game is a dialog between the player and an "oracle", which, for the official game, is the web site.
The player submits a question to the oracle and the oracle responds, using information to which the player does not have access.
In this case the information is the target word.
The question is the player's guess - a 5-letter word - and the response is a score for that word.
The score indicates, for each character, whether it matches the character in the same position in the target or it is in the target in another position or it is not in the target at all.

Using the sample game for Wordle #196 from the Wikipedia page for illustration
"""

# ╔═╡ 9b911b3c-17d5-42d5-b7e5-a4187ce5787c
PlutoUI.Resource(
    "https://upload.wikimedia.org/wikipedia/commons/thumb/e/ec/Wordle_196_example.svg/440px-Wordle_196_example.svg.png",
)

# ╔═╡ 4d5681f9-56a5-4f83-9b31-810bdef26d98
md"""
The target is "rebus".

The player's first guess is "arise" and the response, or score, from the oracle is coded as 🟫🟨🟫🟨🟨 where 🟫 indicates that the letter is not in the target (neither `a` nor `i` occur in "rebus") and 🟨 indicates that the letter is in the target but not at that position.
(I'm using 🟫 instead of a gray square because I can't find a gray square Unicode character.)

The second guess is "route" for which the response is 🟩🟫🟨🟫🟨 indicating that the first letter in the guess occurs as the first letter in the target.
Notice that this guess does not include an "s", which is known from the score of the first guess to be one of the characters in the target.
This guess would not be allowed if playing the game under the "Hard Mode" setting.

Of course, the colors are just one way of summarizing the response to a guess.
Within a computer program it is easier to use an integer to represent each of the 243 = 3⁵ possible scores.
An obvious way of mapping the result to an integer in the (decimal) range 0:242 is by mapping the response for each character to 2 (in target at that position), 1 (in target not at that position), or 0 (not in target) and regarding the pattern as a base-3 number.

In this coding the response for the first guess, "arise", is 01011 in base-3 or 31 in decimal.
The response for the second guess, "route", is 20101 in base-3 or 172 in decimal.

A function to evaluate this score can be written as
"""

# ╔═╡ db787d98-d2a4-4d4a-8e13-8b61e80e48dd
function score(guess, target)
    s = 0
    for (g, t) in zip(guess, target)
        s *= 3
        s += (g == t) ? 2 : Int(g ∈ target)
    end
    return s
end

# ╔═╡ e9e23c37-e4d8-4b77-ab6f-264bcba4ebec
md"""
These numeric scores are not on a scale where "smaller is better" or "larger is better".
(It happens that the best score is 242, corresponding to a perfect match, or five green tiles, but that's incidental.)

The score is just a way of representing each of the 243 patterns that can be produced.

We can convert back to colored tiles if desired using the `tiles` function from the `Wordlegames` package, defined as
```julia
function tiles(sc, ntiles)
    result = Char[]       # initialize to an empty array of Char
    for _ in 1:ntiles     # _ indicates the value of the iterator is not used
        sc, r = divrem(sc, 3)
        push!(result, iszero(r) ? '🟫' : (isone(r) ? '🟨' : '🟩'))
    end
    return String(reverse(result))
end
```
"""

# ╔═╡ b0970feb-103d-4560-b577-f611171e8da6
tiles(31, 5)

# ╔═╡ 39ae1224-9cdf-4d4b-b543-985a5667a601
md"""
## Examining the score function

In the Sherlock Holmes story [The Adventure of Silver Blaze](thttps://en.wikipedia.org/wiki/The_Adventure_of_Silver_Blaze) there is a famous exchange where Holmes remarks on "the curious incident of the dog in the night-time" (see the link).
The critical clue in the case is not what happened but what didn't happen - the dog didn't bark.

Just as Holmes found it interesting that the dog didn't bark, we should find some of the functions in this notebook interesting for what they don't include.
In many of the functions shown here the arguments aren't given explicit types.

Knowing the concrete types of arguments is very important when compiling functions, as is done in Julia, but these functions are written without explicit types.

Consider the `score` function which we reproduce here

```julia
function score(guess, target)
    s = 0
    for (g, t) in zip(guess, target)
        s *= 3
        s += (g == t) ? 2 : Int(g ∈ target)
    end
    return s
end
```

The arguments to `score` can be any type.
In fact, formally they are of an abstract type called `Any`.

So how do we make sure that the actual arguments make sense for this function?
Well, the first thing that is done with the arguments is to pass them to `zip(guess, target)` to produce pairs of values, `g` and `t`, that can be compared for equality, `g == t`.
In a sense `score` delegates the task of checking that the arguments are sensible to the `zip` function.

For those unfamiliar with zipping two or more iterators, we can check what the result is.
"""

# ╔═╡ 8a1cc31b-f747-49b4-9860-a2c381d98c44
collect(zip("arise", "rebus"))

# ╔═╡ 6830dbd6-edc9-4193-8d73-becf3be3855a
md"""
One of the great advantages of dynamically-typed languages with a REPL (read-eval-print-loop) like Julia is that we can easily check what `zip` produces in a couple of examples (or even read the documentation returned by `?zip`, if we are desperate).

The rest of the function is a common pattern - initialize `s`, which will be the result, modify `s` in a loop, and return it.
The Julia expression
```julia
s *= 3
```
indicates, as in several other languages, that `s` is to be multiplied by 3 in-place.

An expression like
```julia
(g == t) ? 2 : Int(g  ∈ target)
```
is a *ternary operator* expression (the name comes from the operator taking three operands).
It evaluates the condition, `g == t`, and returns `2` if the condition is `true`.
If `g == t` is `false` the operator returns the value of the Boolean expression `g  ∈ target`, converted to an `Int`.
(The expression could also be written `g in target`.
In the Julia REPL the `∈` character is created by typing `\in<tab>`.)
The Boolean expression will return `false` or `true`, which is promoted to `0` or `1` for the `+=` operation.

The operation of multiplying by 3 and adding 2 or 1 or 0 is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method) for evaluating a polynomial.

The function is remarkable because it is both general and compact.
Even more remarkable is that it will be very, very fast after its first usage triggers compilation.
That's important because this function will be in a "hot loop".
It will be called many, many times when evaluating the next guess.

(Unfortunately, this version doesn't properly account for cases where a character is repeated in the guess - an example of [Kernighan and Plauger's](https://en.wikipedia.org/wiki/The_Elements_of_Programming_Style) aphorism, "Efficiency often means getting the wrong answer quickly."
We will return to this issue later.)

We won't go into detail about the Julia compiler except to note that compilation is performed for specific *method signatures* (or "method instances") not for general method definitions.

There are several functions and macros in Julia that allow for inspection at different stages of compilation.
One of the most useful is the macro `@code_warntype` which is used to check for situations where type inference has not been successful.
Applying it as
```julia
julia> @code_warntype score("arise", "rebus")
MethodInstance for score(::String, ::String)
  from score(guess, target) in Main at REPL[1]:1
Arguments
  #self#::Core.Const(score)
  guess::String
  target::String
Locals
  @_4::Union{Nothing, Tuple{Tuple{Char, Char}, Tuple{Int64, Int64}}}
  s::Int64
  @_6::Int64
  t::Char
  g::Char
  @_9::Int64
Body::Int64
...
```
shows that type inference is based on concrete types (`String`) for the arguments.

Some argument types are handled more efficiently than others.
Without going in to details we note that we can take advantage of the fact that we have exactly 5 characters and convert the elements of `words` from `String` to `NTuple{5,Char}`, which is an ordered, fixed-length homogeneous collection.

Using the `@benchmark` macro from the `BenchmarkTools` package gives run times of a few tens of nanoseconds for these arguments, and shows that the function applied to the fixed-length collections is faster.
"""

# ╔═╡ e4410d39-6680-4eac-9f70-a425df352998
md"""
New methods can be defined for a generic function like `score`.
A reason for this could be, for example, that the method can more effectively use information from the type.

For example, the `NTuple{N,Char}` type has exactly `N` characters - information that can be used in a loop where we can turn off bounds checking.
"""

# ╔═╡ a031c3db-6630-4deb-a90d-d2b13202dec6
function score(guess::NTuple{N,Char}, target::NTuple{N,Char}) where {N}
    s = 0
    @inbounds for i = 1:N
        s *= 3
        gi = guess[i]
        s += (gi == target[i]) ? 2 : Int(gi ∈ target)
    end
    return s
end

# ╔═╡ 9e1d6017-f79e-49b5-916e-8e3a9e49ab1d
score("arise", "rebus")

# ╔═╡ a69f62bd-ed6d-42e7-8405-9321d13a52d1
@benchmark score(guess, target) setup = (guess = "arise"; target = "rebus")

# ╔═╡ 1b9261b3-c90c-4db1-9e70-51bf364072b8
md"""
This method returns the same result as the other method, only faster.
"""

# ╔═╡ 1887c091-46a0-4e8c-9a66-2ab17968729c
score(('a', 'r', 'i', 's', 'e'), ('r', 'e', 'b', 'u', 's'))

# ╔═╡ 39202f1a-e7d6-4961-9776-ec67c3fd9743
@benchmark score(guess1, target1) setup =
    (guess1 = ('a', 'r', 'i', 's', 'e'); target1 = ('r', 'e', 'b', 'u', 's'))

# ╔═╡ 2b73c7f4-0203-4e2b-b634-eb310152834b
md"""
## Repeated characters in the guess

The simple `score` methods shown above don't give the correct score (meaning the score that would be returned on the web site) when there are repeated characters in the guess.
For example, a guess of `"sheer"` for the target `"super"` is scored as
"""

# ╔═╡ 79f0ceab-0c9f-42f7-b3ec-a58b5ac1c083
tiles(score("sheer", "super"), 5)

# ╔═╡ b92abc9c-2a86-47ab-a6e8-127c4d7de52e
md"""
but the score should be `"🟩🟫🟫🟩🟩"` because there is only one `e` in the target `"super"`.
In a case like this where a character occurs multiple times in a guess but only once in the target the rules about which position in the guess is marked are that "correct position" takes precedence over "in the word" and, if none of the guess positions are correct, then leftmost takes precedence.

This makes for a considerably more complex score evaluation.
Essentially there have to be two passes over the score and target - the first to check for correct position and the second to check for "in the target, not in the correct position".

However, the simple algorithm in the current `score` methods works if there are no duplicate characters in the guess.
Thus it is probably worthwhile checking for duplicates, using the simple function

```julia
function hasdups(guess::NTuple{N,Char}) where {N}
    @inbounds for i in 1:(N - 1)
        gi = guess[i]
        for j in (i + 1):N
            gi == guess[j] && return true
        end
    end
    return false
end
```

and choose the simple scoring algorithm when there are no duplicates.

In the [Wordlegames](https://github.com/dmbates/Wordlegames.jl) package these operations are combined in a `scorecolumn!` function that updates a vector of scores on a single guess against a vector of targets.

```julia
function scorecolumn!(
    col::AbstractVector{<:Integer},
    guess::NTuple{N,Char},
    targets::AbstractVector{NTuple{N,Char}},
) where {N}
    if axes(col) ≠ axes(targets)
        throw(
            DimensionMismatch(
                "axes(col) = $(axes(col)) ≠ $(axes(targets)) = axes(targets)",
            )
        )
    end
    if hasdups(guess)
        onetoN = (1:N...,)           # 1:N as a Tuple
        svec = zeros(Int, N)         # scores for characters in guess
        unused = trues(N)            # unused positions in targets[i]
        @inbounds for i in axes(targets, 1)
            targeti = targets[i]
            fill!(unused, true)      # reset to all unused
            fill!(svec, 0)           # reset to all guess characters not in target
            for j = 1:N              # first pass for target in same position
                if guess[j] == targeti[j]
                    unused[j] = false
                    svec[j] = 2
                end
            end
            for j = 1:N              # second pass for match in unused position
                if iszero(svec[j])
                    for k in onetoN[unused]
                        if guess[j] == targeti[k]
                            svec[j] = 1
                            unused[k] = false
                            break
                        end
                    end
                end
            end
            sc = 0                   # Horner's method to evaluate the score
            for s in svec
                sc *= 3
                sc += s
            end
            col[i] = sc
        end
    else                             # simplified alg. for guess w/o duplicates
        @inbounds for i in axes(targets, 1)
            sc = 0
            targeti = targets[i]
            for j = 1:N
                sc *= 3
                gj = guess[j]
                sc += (gj == targeti[j]) ? 2 : Int(gj ∈ targeti)
            end
            col[i] = sc
        end
    end
    return col
end
```

This is "production code" which has gone through several refinement steps so it may seem a bit daunting at first.
However, we can break it down.

First, does it give the desired result?
"""

# ╔═╡ 2f32604f-b64c-47f5-bf14-a3cbde36cc13
scores1 = zeros(Int, 1)  # initialize a vector of 1 integer to zero

# ╔═╡ f899d69c-ee68-47e3-963f-90bc401558b2
scorecolumn!(scores1, ('s', 'h', 'e', 'e', 'r'), [('s', 'u', 'p', 'e', 'r')])

# ╔═╡ 643bb3f5-c8e7-446c-b7d4-ea40952c4f59
tiles(first(scores1), 5)

# ╔═╡ 4742dec1-1a46-4095-bc49-467182beb163
md"""
We see that the call to `scorecolumn!` overwrites the contents of the `scores1` vector with the score for the guess on the first (and only) target.

Thus `scorecolumn!` is a "mutating function", meaning that it changes the contents of one or more of its arguments.
By convention we give such functions names ending in `"!"`, as a warning to the user that the function may mutate its arguments.
(This is merely a convention; the `"!"` has no syntactic significance.)
Furthermore, the convention is to list any arguments that may be modified first.

The reason this function is called `scorecolumn!` is because the scores for all possible guesses on all possible targets are evaluated and cached as a matrix in a `GamePool` object.
This may seem extravagant but most methods for determining an initial guess algorithmically will end up evaluating all these scores so it makes sense to save them in an array.
In this case the rows correspond to targets and the columns to guesses and evaluating the scores for a single guess against all possible targets updates a column of this matrix.

A section from the upper left corner of this matrix
"""

# ╔═╡ 225c813a-e63e-44fc-bc08-caf4cf36e21b
view(wordle.allscores, 1:7, 1:10)

# ╔═╡ 875f7d47-b71f-4f67-90ae-755b00e2b464
md"""
shows that the scores, which are in the range `0:242`, are stored as unsigned, 8-bit integers to conserve storage.
Even so, the storage required is (2315)² bytes, or over 5 megabytes.
"""

# ╔═╡ 7c6e0762-de57-4041-a683-7ed6f0715057
Base.summarysize(wordle.allscores)

# ╔═╡ 178cfcd4-e198-4e44-b5bd-7da6c8b21456
md"""
Five megabytes is not a large amount of memory by today's standards, but for games with larger pools of guesses or targets the storage may start to mount up.
In those cases there is provision for [memory-mapping](https://en.wikipedia.org/wiki/Memory-mapped_file) the array.
The evaluation of the array is multi-threaded when Julia is running with multiple threads.

The scores in the first column,
"""

# ╔═╡ 90f09143-4046-4e2f-967e-1df115411b32
tiles.(view(wordle.allscores, 1:7, 1), 5)

# ╔═╡ 628fb954-e44f-4c6d-af03-c41e9e2894d5
md"""
are for the first guess, "aback", against the first 7 targets
"""

# ╔═╡ c976042e-2e91-401b-acfa-88199c22326c
[String(collect(t)) for t in view(wordle.targetpool, 1:7)]

# ╔═╡ 197e3fbe-6172-47b8-99c5-c8293223747c
md"""
The `scorecolumn!` function itself uses the `axes` function in several places.
By default Julia uses 1-based indexing but other forms of indexing are allowed.
(I am obligated at this point to mention [StarWarsArrays](https://github.com/giordano/StarWarsArrays.jl) which begins indexing at 4, 5, 6 then 1, 2, 3 then 7, 8, and 9.)

The call to `axes(targets, 1)` returns the indices in the first (and only) axis of the `target` vector.
The `col` and `targets` arguments are typed as `AbstractVector`, not `Vector`, because `Vector` is a concrete, specific type and we wish to allow for "vector-like" objects such as a one-dimensional view in a multi-dimensional array.

The call to `scorecolumn!` in the constructor for a `GamePool` is in the code segment

```julia
    S = scoretype(N)
    vtargs = view(guesspool, validtargets)
    allscores = Array{S}(undef, length(vtargs), length(guesspool))
    Threads.@threads for j in axes(allscores, 2)
        scorecolumn!(view(allscores, :, j), guesspool[j], vtargs)
    end
```

There are two arrays, `svec` and `unused` allocated within the `scorecolumn!` function when guesses have repeated characters.
These are very small arrays but nonetheless we would want to minimize the number of allocations if feasible.
This is why the check for duplicate characters is carried out and the allocation of these arrays is done only once per column of `allscores`.
The allocation is done within the function so that it can be called from multiple threads simultaneously without the threads interfering with each other.

The branch for a guess without duplicates is still much faster than for a guess with duplicate characters but neither case is horribly slow.
"""

# ╔═╡ de2cac98-2ca9-4d27-8b02-ff152a73df7e
@benchmark scorecolumn!(col, ('a', 'r', 'i', 's', 'e'), $(wordle.guesspool)) setup =
    (col = zeros(UInt8, length(wordle.guesspool)))

# ╔═╡ 04cbc640-70be-402b-8a7e-3e9fc8e2cb00
md"""
Notice that there are no allocations of memory when there are no duplicated characters in the guess.
There are allocations, and consequently some garbage collection (GC), when the guess has duplicated characters.
"""

# ╔═╡ 01d24948-480d-4eb5-9cc4-7c1e806f22fd
@benchmark scorecolumn!(col1, ('a', 'b', 'a', 'c', 'k'), $(wordle.guesspool)) setup =
    (col1 = zeros(UInt8, length(wordle.guesspool)))

# ╔═╡ 2c6d5aa5-7a23-4c67-8ff0-27b95ef58375
md"""
## Conclusion

These few examples have introduced, at least in passing, several advanced programming concepts - multi-threading, memory-mapping, control of storage allocation and garbage collection - that one typically would not associate with a dynamically-typed, REPL-based language like Julia.

Of course, all of these facilities are available in compiled languages like C/C++ or Rust but usually without the "rapid development and testing" capability of a language like Julia.

Julia provides a wide range of tools so that a programmer can start at a very simple level, like the original `score` method and refine as needed to reach speeds previously only achievable with compiled, statically-typed languages.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Wordlegames = "1cb69566-e1cf-455f-a587-fd79a2e00f5a"

[compat]
BenchmarkTools = "~1.3.1"
PlutoUI = "~0.7.38"
Wordlegames = "~0.3.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.8.0-beta3"
manifest_format = "2.0"
project_hash = "045a6e2c86e27f8d763b85c40450eb9e85eed07e"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.3.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.2+0"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.3"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "7.81.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.10.2+0"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.20+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "621f4f3b4977325b9128d5fae7a8b4829a0c2222"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.4"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.8.0"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.1"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.0"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Wordlegames]]
deps = ["AbstractTrees", "DataFrames", "Random", "Tables"]
git-tree-sha1 = "4c463de78d2f3f9447b695e241eba43cc945f866"
uuid = "1cb69566-e1cf-455f-a587-fd79a2e00f5a"
version = "0.3.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.12+1"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.1.0+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"
"""

# ╔═╡ Cell order:
# ╟─f9816eb0-078e-41f2-bf7c-29ed5c663fd0
# ╟─a9678ae8-48c0-4bb9-920d-7de176105bd0
# ╟─4eabf46e-f75c-4b15-8583-6abe17c0fd85
# ╠═49a44e94-29c0-4004-8d0b-e37008eafc1c
# ╟─03b41118-5756-447b-9900-15af37667529
# ╠═73e3e38b-0935-4ef9-97b1-8a8cdfa2feb7
# ╠═4c5e781c-fd10-425f-95d2-1fec5bf1015f
# ╟─2046a4ca-c4d3-4804-9373-930d9cbb58fd
# ╠═560ad4d4-945c-4d21-8947-cac68c1880a1
# ╟─ddd7f90c-7cad-4279-a275-b1800b4ae83b
# ╠═49500976-7ab8-43cc-9902-b25ad73ff42e
# ╠═f6681018-f729-48de-9342-a5bcb151aca2
# ╟─22f7e1ad-3a2b-4a88-8f81-91feada88621
# ╠═5787eef2-7a2c-46ad-8fa0-7313c07be84b
# ╠═8df0e0c0-cdae-450c-abc1-5cea1d520ec9
# ╠═8580bdcf-24ce-4eb6-ac5a-4c77b729b38b
# ╠═9ed9e8ab-3ebc-4c6e-b2c3-46ad1987561e
# ╟─fdd6631e-8b3a-4b4f-b8a7-11135f13cc23
# ╟─9b911b3c-17d5-42d5-b7e5-a4187ce5787c
# ╟─4d5681f9-56a5-4f83-9b31-810bdef26d98
# ╠═db787d98-d2a4-4d4a-8e13-8b61e80e48dd
# ╠═9e1d6017-f79e-49b5-916e-8e3a9e49ab1d
# ╟─e9e23c37-e4d8-4b77-ab6f-264bcba4ebec
# ╠═b0970feb-103d-4560-b577-f611171e8da6
# ╟─39ae1224-9cdf-4d4b-b543-985a5667a601
# ╠═8a1cc31b-f747-49b4-9860-a2c381d98c44
# ╟─6830dbd6-edc9-4193-8d73-becf3be3855a
# ╠═a69f62bd-ed6d-42e7-8405-9321d13a52d1
# ╟─e4410d39-6680-4eac-9f70-a425df352998
# ╠═a031c3db-6630-4deb-a90d-d2b13202dec6
# ╟─1b9261b3-c90c-4db1-9e70-51bf364072b8
# ╠═1887c091-46a0-4e8c-9a66-2ab17968729c
# ╠═39202f1a-e7d6-4961-9776-ec67c3fd9743
# ╟─2b73c7f4-0203-4e2b-b634-eb310152834b
# ╠═79f0ceab-0c9f-42f7-b3ec-a58b5ac1c083
# ╟─b92abc9c-2a86-47ab-a6e8-127c4d7de52e
# ╠═2f32604f-b64c-47f5-bf14-a3cbde36cc13
# ╠═f899d69c-ee68-47e3-963f-90bc401558b2
# ╠═643bb3f5-c8e7-446c-b7d4-ea40952c4f59
# ╟─4742dec1-1a46-4095-bc49-467182beb163
# ╠═225c813a-e63e-44fc-bc08-caf4cf36e21b
# ╟─875f7d47-b71f-4f67-90ae-755b00e2b464
# ╠═7c6e0762-de57-4041-a683-7ed6f0715057
# ╟─178cfcd4-e198-4e44-b5bd-7da6c8b21456
# ╠═90f09143-4046-4e2f-967e-1df115411b32
# ╟─628fb954-e44f-4c6d-af03-c41e9e2894d5
# ╠═c976042e-2e91-401b-acfa-88199c22326c
# ╟─197e3fbe-6172-47b8-99c5-c8293223747c
# ╠═de2cac98-2ca9-4d27-8b02-ff152a73df7e
# ╟─04cbc640-70be-402b-8a7e-3e9fc8e2cb00
# ╠═01d24948-480d-4eb5-9cc4-7c1e806f22fd
# ╟─2c6d5aa5-7a23-4c67-8ff0-27b95ef58375
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: dmbates/WordleTutorials

max_issues_repo_name: dmbates/WordleTutorials

max_forks_repo_name: dmbates/WordleTutorials


text:
using Distributed
using SharedArrays

@everywhere import FIRRemez
@everywhere import PyPlot
@everywhere import Random
@everywhere import Printf

import FileIO
import Printf

using FFTW
#import JLD
import XLSX
import BSON

fig_num = 1
PyPlot.close("all")

Random.seed!(25)

fig_num = 1

L = 25
#L = 1000
if length(ARGS) > 1
    L = tryparse(Int, ARGS[1])
end

## low-pass.
passband_ω_Float64 = 0.025*π # change this.
if length(ARGS) > 2
    passband_ω_Float64 = tryparse(Float64, ARGS[2])
end

stopband_ω_Float64 = 0.025*2*π # change this.
if length(ARGS) > 3
    stopband_ω_Float64 = tryparse(Float64, ARGS[3])
end

save_name_tag = "lowpass"
if length(ARGS) > 4
    save_name_tag = ARGS[4]
end

passband_ω = convert(BigFloat, passband_ω_Float64)
stopband_ω = convert(BigFloat, stopband_ω_Float64)

f = xx->FIRRemez.cosinetransitionlowpassfunction(xx,passband_ω,stopband_ω)

println("passband is ", passband_ω)
println("stopband is ", stopband_ω)

filter_type_num = 1 # symmetric, odd length.

savefig_flag = false
plot_output_flag = false

output_folder = "./output"
opt_params_file_name = "../default/config_opt.txt" # location of the configuration file.

# omit opt_params_file_name is using default configuration file provided by the package.
h_BigFloat, X_BigFloat, fig_num = FIRRemez.frontendfilterdesign(filter_type_num, f,
                            L, fig_num, savefig_flag, plot_output_flag,
                            opt_params_file_name)

h = convert(Vector{Float64},h_BigFloat)
X = convert(Vector{Float64}, X_BigFloat)

Printf.@printf("Order of Chebyshev is %d, length of filter is %d, FIR type %d\n", L, length(h), filter_type_num)

band_info_string = Printf.@sprintf("passband_%f_stopband_%f", passband_ω_Float64, stopband_ω_Float64)
save_name = Printf.@sprintf("type%d_%s_L%d_%s.bson",
    filter_type_num, save_name_tag, L, band_info_string)

passband = passband_ω_Float64
stopband = stopband_ω_Float64

BSON.bson(joinpath(output_folder, save_name),
    h = h, X = X, passband = passband, stopband = stopband)
#JLD.save(joinpath(output_folder, save_name), "h" = h, "X" = X, "passband" = passband, "stopband" = stopband)
#
save_name = Printf.@sprintf("type%d_%s_L%d_%s.xlsx",
    filter_type_num, save_name_tag, L, band_info_string)
XLSX.openxlsx(joinpath(output_folder, save_name), mode="w") do xf
    sheet = xf[1]
    XLSX.rename!(sheet, "filter")

    sheet["A1", dim = 1] = ["passband"; passband]
    sheet["B1", dim = 1] = ["passband"; passband]
    sheet["C1", dim = 1] = ["Chebyshev node positions"; X]
    sheet["D1", dim = 1] = ["coefficients"; h]
end

#fig_num = FIRRemez.plotmagnitudersp(h, fig_num, "filter's magnitude response")

import Plots
Plots.plotly()

ω_set_fft, DFT_evals, ω_set, DTFT_evals = FIRRemez.getfreqrsp(h)

mag_rsp_G = abs.(DTFT_evals)
mag_rsp_fft = abs.(DFT_evals)

plot_obj = Plots.plot( ω_set_fft,
mag_rsp_fft,
title = "Magnitude spectrum",
label = "DFT",
seriestype = :line,
ticks = :native,
hover = ω_set_fft,
linewidth = 4,
size = (1200, 900))

Plots.plot!(plot_obj, ω_set, mag_rsp_G, label = "DTFT",
seriestype = :line,
linestyle = :dash,
linewidth = 4)

save_name = Printf.@sprintf("type%d_%s_L%d_%s.html",
    filter_type_num, save_name_tag, L, band_info_string)
Plots.savefig(plot_obj, joinpath(output_folder, save_name))

ext: jl

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
module SequentialPointProcesses
import PointPatternStatistics: PointPattern
import Distributions
using Requires

include("models.jl")
include("logfk.jl")
include("rand.jl")
include("overlappingdiscs.jl")

Distributions.logpdf(model::SequentialPointProcess, pp::PointPattern, int::Integer) =
    logdens(model, pp.data, pp.window, Options(nx=int))
Distributions.logpdf(model::SequentialPointProcess, pp::PointPattern, int::Options) =
    logdens(model, pp.data, pp.window, int)

export Softcore,Mixture,Uniform, OverlappingDiscs

cuda_available = false
function __init__()
    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
end

end

max_stars_repo_name: mikkoku/SequentialPointProcesses

max_issues_repo_name: mikkoku/SequentialPointProcesses

max_forks_repo_name: mikkoku/SequentialPointProcesses


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkBgpCommunities <: SwaggerModel
    virtualNetworkCommunity::Any # spec type: Union{ Nothing, String } # spec name: virtualNetworkCommunity
    regionalCommunity::Any # spec type: Union{ Nothing, String } # spec name: regionalCommunity

    function VirtualNetworkBgpCommunities(;virtualNetworkCommunity=nothing, regionalCommunity=nothing)
        o = new()
        validate_property(VirtualNetworkBgpCommunities, Symbol("virtualNetworkCommunity"), virtualNetworkCommunity)
        setfield!(o, Symbol("virtualNetworkCommunity"), virtualNetworkCommunity)
        validate_property(VirtualNetworkBgpCommunities, Symbol("regionalCommunity"), regionalCommunity)
        setfield!(o, Symbol("regionalCommunity"), regionalCommunity)
        o
    end
end # type VirtualNetworkBgpCommunities

const _property_map_VirtualNetworkBgpCommunities = Dict{Symbol,Symbol}(Symbol("virtualNetworkCommunity")=>Symbol("virtualNetworkCommunity"), Symbol("regionalCommunity")=>Symbol("regionalCommunity"))
const _property_types_VirtualNetworkBgpCommunities = Dict{Symbol,String}(Symbol("virtualNetworkCommunity")=>"String", Symbol("regionalCommunity")=>"String")
Base.propertynames(::Type{ VirtualNetworkBgpCommunities }) = collect(keys(_property_map_VirtualNetworkBgpCommunities))
Swagger.property_type(::Type{ VirtualNetworkBgpCommunities }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkBgpCommunities[name]))}
Swagger.field_name(::Type{ VirtualNetworkBgpCommunities }, property_name::Symbol) =  _property_map_VirtualNetworkBgpCommunities[property_name]

function check_required(o::VirtualNetworkBgpCommunities)
    (getproperty(o, Symbol("virtualNetworkCommunity")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualNetworkBgpCommunities }, name::Symbol, val)
end


text:
type Spectrum
  gainStep::Int
  tStart::DateTime
  commandedMass::Int
  p_pt::Int
  row_A::Vector{Int}
  row_B::Vector{Int}
  pixelGain_A::Vector{Float64}
  pixelGain_B::Vector{Float64}
  baseline_A::Vector{Float64}
  baseline_B::Vector{Float64}
  ionsPerSpectrum_A::Vector{Float64}
  ionsPerSpectrum_B::Vector{Float64}
  peakIndices_A::Vector{Int}
  peakIndices_B::Vector{Int}
  peakAmplitudes_A::Vector{Float64}
  peakAmplitudes_B::Vector{Float64}
end


text:
# abstract types
"""
    abstract type AbstractFrontendAction <: Any
Supertype for `FrontendAction`s.
"""
abstract type AbstractFrontendAction end

# the file hierarchy is exactly the same as Clang, please refer to Clang's src for docs.
# AST
include("AST/ASTConsumer.jl")
include("AST/ASTContext.jl")
include("AST/Decl.jl")
include("AST/DeclarationName.jl")
include("AST/DeclBase.jl")
include("AST/DeclCXX.jl")
include("AST/DeclGroup.jl")
include("AST/DeclTemplate.jl")
include("AST/NestedNameSpecifier.jl")
include("AST/TemplateBase.jl")
include("AST/TemplateName.jl")
include("AST/Type.jl")

# Basic
include("Basic/CodeGenOptions.jl")
include("Basic/Diagnostic.jl")
include("Basic/DiagnosticIDs.jl")
include("Basic/DiagnosticOptions.jl")
include("Basic/FileEntry.jl")
include("Basic/FileManager.jl")
include("Basic/IdentifierTable.jl")
include("Basic/LangOptions.jl")
include("Basic/SourceLocation.jl")
include("Basic/SourceManager.jl")
include("Basic/TargetInfo.jl")
include("Basic/TargetOptions.jl")

# CodeGen
include("CodeGen/CodeGenAction.jl")
include("CodeGen/ModuleBuilder.jl")

# Frontend
include("Frontend/CompilerInstance.jl")
include("Frontend/CompilerInvocation.jl")
include("Frontend/FrontendOptions.jl")
include("Frontend/TextDiagnosticPrinter.jl")

# Lex
include("Lex/HeaderSearch.jl")
include("Lex/HeaderSearchOptions.jl")
include("Lex/Preprocessor.jl")
include("Lex/PreprocessorOptions.jl")
include("Lex/Token.jl")

# Parse
include("Parse/ParseAST.jl")
include("Parse/Parser.jl")

# Sema
include("Sema/DeclSpec.jl")
include("Sema/Lookup.jl")
include("Sema/Scope.jl")
include("Sema/Sema.jl")


text:
using MCPhylo

## Get the directory in which the example CODA files are saved
dir = dirname(@__FILE__)

## Read the MCMC sampler output from the files
c1 = readcoda(joinpath(dir, "line1.out"), joinpath(dir, "line1.ind"))
c2 = readcoda(joinpath(dir, "line2.out"), joinpath(dir, "line2.ind"))

## Concatenate the resulting chains
c = cat(c1, c2, dims=3)

## Compute summary statistics
describe(c)

#mark that we got to the end of the test file succesfully
@test true


text:
# Unit Test for by default supported convex sets and their functions
using COSMO, Test, Random, LinearAlgebra

rng = Random.MersenneTwister(13131)

@testset "Convex Sets" begin
tol = 1e-4

    @testset "Create and project" begin

    # Zero Cone
    zset = COSMO.ZeroSet(10)
    x = randn(rng,10)
    COSMO.project!(view(x, 1:length(x)), zset)
    @test norm(x,Inf) == 0.

    # Positive Orthant R+
    nonnegatives = COSMO.Nonnegatives(10)
    x = randn(rng, 10)
    COSMO.project!(view(x, 1:length(x)), nonnegatives)
    @test minimum(x) >= 0.

    # Box
    l = -1 * ones(10)
    u = 1 * ones(10)
    box = COSMO.Box(l, u)
    box = COSMO.Box{Float64}(10)
    box.l .= l
    box.u .= u
    x = 100 * randn(rng, 10)
    COSMO.project!(view(x, 1:length(x)), box)
    @test minimum(x) >= -1. && maximum(x) <= 1.

    # Second Order (Lorentz) cones
    soc = COSMO.SecondOrderCone(10)
    x = 10 * randn(rng, 9)
    t = norm(x, 2) - 0.5
    x = [t; x]

    COSMO.project!(view(x, 1:length(x)), soc)
    @test norm(x[2:10], 2) <= x[1]

    # Positive Semidefinite cones
    psd = COSMO.PsdCone(16)
    X = Symmetric(randn(rng, 4, 4))
    X = X - 4 * Matrix(1.0I, 4, 4)
    x = vec(X)
    COSMO.project!(view(x, 1:length(x)), psd)
    @test minimum(eigen(reshape(x, 4, 4)).values) >= -1e-9

    C = COSMO.CompositeConvexSet([COSMO.ZeroSet(10), COSMO.Nonnegatives(10)])
    x = -rand(20)
    xs = COSMO.SplitVector(x, C)
    COSMO.project!(xs, C)
    @test norm(x, Inf) == 0.

    end

    @testset "in_dual Functions" begin

    # Dual of zero cone
    x = randn(rng, 10)
    convex_set = COSMO.ZeroSet(10)
    @test COSMO.in_dual(view(x, 1:length(x)), convex_set, tol)

    # Dual of Positive Orthant R+ (self-dual)
    xpos = rand(rng, 10)
    xneg = -rand(rng, 10)
    xzeros = zeros(10)
    convex_set = COSMO.Nonnegatives(10)
    @test COSMO.in_dual(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_dual(view(xneg,1:length(xneg)), convex_set, tol)
    @test COSMO.in_dual(view(xzeros,1:length(xzeros)), convex_set,tol)

    #TODO: Dual of Box [important!]

    # Dual of Second Order Cone (self-dual)
    tol = 1e-4
    x = randn(rng, 9)
    t = norm(x, 2)
    xpos = [t + 0.5; x]
    xneg = [t - 0.5; x]
    convex_set = COSMO.SecondOrderCone(10)
    @test COSMO.in_dual(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_dual(view(xneg, 1:length(xneg)), convex_set, tol)

    # Dual of Positive Semidefinite Cone (Square) (self-dual)
    tol = 1e-4
    X = Symmetric(randn(rng, 4, 4))
    Xpos = X + 4 * Matrix(1.0I, 4, 4)
    Xneg = X - 4 * Matrix(1.0I, 4, 4)
    xpos = vec(Xpos)
    xneg = vec(Xneg)
    convex_set = COSMO.PsdCone(16)
    @test COSMO.in_dual(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_dual(view(xneg, 1:length(xneg)), convex_set, tol)

    # Dual of Positive Semidefinite Cone (Triangle)
    convex_set = COSMO.PsdConeTriangle(10)
    xpos = zeros(10)
    xneg = zeros(10)
    COSMO.extract_upper_triangle!(Xpos, xpos, sqrt(2))
    COSMO.extract_upper_triangle!(Xneg, xneg, sqrt(2))
    @test COSMO.in_dual(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_dual(view(xneg, 1:length(xneg)), convex_set, tol)

    end

    @testset "inPolRec Functions" begin

    # Polar Recession cone of zero cone
    xpos = zeros(10)
    xneg = randn(rng, 10)
    convex_set = COSMO.ZeroSet(10)
    @test COSMO.in_pol_recc(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_pol_recc(view(xneg, 1:length(xneg)), convex_set, tol)

    # Polar Recession cone of Positive Orthant R+
    xpos = -rand(rng, 10)
    xneg = rand(rng, 10)
    convex_set = COSMO.Nonnegatives(10)
    @test COSMO.in_pol_recc(view(xpos,1:length(xpos)), convex_set, tol)
    @test !COSMO.in_pol_recc(view(xneg,1:length(xneg)), convex_set, tol)

    #TODO: Polar Recc of Box [important!]

    # Polar Recc of Second Order Cone
    tol = 1e-4
    x = randn(rng, 9)
    t = norm(x, 2)
    xpos = [-t - 0.5; x]
    xneg = [-t + 0.5; x]
    convex_set = COSMO.SecondOrderCone(10)
    @test COSMO.in_pol_recc(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_pol_recc(view(xneg, 1:length(xneg)), convex_set, tol)

    # Polar Recc of Positive Semidefinite Cone (Square)
    tol = 1e-4
    X = Symmetric(randn(rng, 4, 4))
    Xpos = X - 20 * Matrix(1.0I, 4, 4)
    Xneg = X + 4 * Matrix(1.0I, 4, 4)
    xpos = vec(Xpos)
    xneg = vec(Xneg)
    convex_set = COSMO.PsdCone(16)
    @test COSMO.in_pol_recc(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_pol_recc(view(xneg, 1:length(xneg)), convex_set, tol)

    # Polar Recc of Positive Semidefinite Cone (Triangle)
    COSMO.extract_upper_triangle!(Xpos, xpos, sqrt(2))
    COSMO.extract_upper_triangle!(Xneg, xneg, sqrt(2))
    convex_set = COSMO.PsdConeTriangle(10)
    @test COSMO.in_pol_recc(view(xpos, 1:length(xpos)), convex_set, tol)
    @test !COSMO.in_pol_recc(view(xneg, 1:length(xneg)), convex_set, tol)
   end

   @testset "Set Utilities" begin
    convex_set = COSMO.PsdConeTriangle(3)
    composite_set = COSMO.CompositeConvexSet([COSMO.ZeroSet(2), COSMO.Nonnegatives(1)])
    @test COSMO.get_subset(convex_set, 1) == convex_set
    @test_throws DimensionMismatch COSMO.get_subset(convex_set, 2)
    @test COSMO.get_subset(composite_set, 1) == COSMO.ZeroSet(2)
    @test_throws ErrorException COSMO.CompositeConvexSet([convex_set; composite_set])
   end

end

nothing


text:
using DEBTrait, Test
using Random
Random.seed!(123)

################# (1 microbe, 1 reserve), no enzyme production
n_microbes      = 1
E               = rand(n_microbes)
V               = rand(n_microbes)
Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(1,n_microbes), 1.0*ones(n_microbes))
p               = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
p.y_EM          = fill(1.0, n_microbes)
p.y_EX          = fill(1.0, n_microbes)
r               = growth!(zeros(size(V,1)), E, V, p)
r_analytic      = @. (p.k_E*E/V - p.k_M*p.y_EM)/(E/V + p.y_EV)
@test r ≈ r_analytic atol=1e-10

out             = growth_production(r, E, V, p)
x_analytic      = @. p.α*p.y_EV*r_analytic/p.y_EX
@test out[1] ≈ x_analytic atol=1e-10

################# (1 microbe, 2 reserves)
n_reserves      = 2
E               = fill(0.5,n_microbes*n_reserves)
E               = reshape(E, n_reserves, n_microbes)
V               = fill(1.0, n_microbes)
p               = PMetabolismCN{BaseGenome, Array{Float64,2}, Array{Float64,1}}(Isolate)
p.y_EX          = fill(1.0, n_reserves, n_microbes)
p.k_M           = fill(1.0, n_microbes)
p.α             = fill(0.1, n_microbes)
p.α_n           = fill(0.1, n_microbes)
p.y_EM[:,1]     = [0.3, 1.0]
p.y_EV[:,1]     = [1.2, 1.0]
p.k_E           = [1.0]
r               = growth!(zeros(size(E,2)), E, V, p)
@test r[1] ≈ -0.333 atol=1e-2
out             = growth_production(r, E, V, p)
@test out[3][:,1] ≈ [0.3667, 0.0] atol=1e-4

################# (2 microbes, 2 reserves)
n_microbes      = 2
n_reserves      = 2
n_substrates    = 1
Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(n_substrates,n_microbes), 1.0*ones(n_microbes))
E               = fill(0.5,n_microbes*n_reserves)
E               = reshape(E, n_reserves, n_microbes)
V               = fill(1.0, n_microbes)
p               = PMetabolismCN{BaseGenome, Array{Float64,2}, Array{Float64,1}}(Isolate)
p.α_n           = 1e-4*ones(n_microbes)
p.y_EX          = 1.2*ones(n_reserves, n_microbes)
p.y_EM          = ones(n_reserves, n_microbes)
r               = growth!(zeros(size(E,2)), E, V, p)
@test size(r,1) == 2
out             = growth_production(r, E, V, p)
@test size(out[1],1) == (n_microbes)


text:
using Pkg
Pkg.activate(joinpath(dirname(@__FILE__), "../../../"))
using Test, OrdinaryDiffEq, HDF5, ModelConstructors, Plots, NLsolve, Statistics, BenchmarkTools
include(joinpath(dirname(@__FILE__), "../../../src/includeall.jl"))
include("fixedpoint_timestep.jl")

start_time = time_ns()

default = false
use_nlsolve = true # slower than default by .3 seconds, plus 2 extra steps
method = :newton
my_upwind = true
uniform = false
m_anderson = 3
bc = [(0., 1.)] # Vector{Tuple{Float64, Float64}}(undef, 0)

if default
    use_nlsolve = false
end

## Parameters
ρₑ = 0.03
ρₕ = 0.03
τ  = .005
aₑ = .15
aₕ = .7 * aₑ
σ  = .03
γ  = 2.
χ₁ = 1.
χ₂ = 5.
δ  = 0.03
Δ  = .6
θ = (ρₑ = ρₑ, ρₕ = ρₕ, τ = τ, aₑ = aₑ, aₕ = aₕ, σ = σ, γ = γ, χ₁ = χ₁, χ₂ = χ₂, δ = δ, Δ = Δ)
pnorm = Inf
T = 100
tol = 1e-3
update_hyp = .8 # 1 is all weight on the ratio of a rolling window of error averages, 0 is none
reduce_hyp = .95 # proportion by which to reduce Δ if the error grows in size
discounting = 5

## Grid
# uneven grid on [0, 1] with 1000 points, with more points near 0 and 1

N   = 100 # 1000
zz  = collect(range(0.01, stop = 0.99, length = N))
# η   = zz.^1.1
η   = zz
dη  = diff(η)
stategrid = StateGrid(OrderedDict(:η => η))

if uniform
    L₁ = UpwindDifference(1, 1, dη[1], N)
    L₂ = CenteredDifference(2, 2, dη[1], N)
    Qbc = RobinBC((0., 1., 0.), (0., 1., 0.), dη[1])
else
    dx = nonuniform_grid_spacing(stategrid, bc)
    L₁ = UpwindDifference(1, 1, dx, N)
    L₂ = CenteredDifference(2, 2, dx, N)
    Qbc = RobinBC((0., 1., 0.), (0., 1., 0.), dx)
end

default_err = fill(NaN, N, 3)
nlsolve_newton_err = fill(NaN, N, 3)
nlsolve_trustreg_err = fill(NaN, N, 3)

## Terminal conditions for value functions V and Vₕ
vₑ = aₑ^(-γ) .* η .^ (1. - γ)
vₕ = aₑ^(-γ) .* (1. .- η) .^ (1. - γ)
# vₑ = 100 * η .^ (1. - γ)
# vₕ = 100 * (1. .- η) .^ (1. - γ)
vₑ₀ = similar(vₑ)
vₕ₀ = similar(vₕ)
vₑ₁ = similar(vₑ)
vₕ₁ = similar(vₕ)
Q   =  ones(N)
Qp  = zeros(N)
SSQ = zeros(N)
Psi = zeros(N)
Iota   = zeros(N)
Φ   = zeros(N)
A   = zeros(N)
S   = zeros(N)
MU  = zeros(N)
Mvₑ  = zeros(N)
Mvₕ = zeros(N)
Mvₑfd  = zeros(N)
Mvₕfd = zeros(N)
G   = zeros(N)
Svₑ  = zeros(N)
Svₕ = zeros(N)
CN  = zeros(N)
C_N = zeros(N)

#=println("Time Anderson acceleration")
@btime begin
out = nlsolve((F, x) -> fixedpoint_timestep(F, x, η, dη, Q, Qp, Psi, SSQ, MU, S, Svₑ, Svₕ, Mvₑ, Mvₕ, vₑ₁, vₕ₁, θ;
                                            use_nlsolve = use_nlsolve, method = method, my_upwind = my_upwind),
              vcat(vₑ, vₕ), method = :anderson, m = m_anderson, ftol = tol)
end=#
out = nlsolve((F, x) -> fixedpoint_timestep(F, x, η, dη, Q, Qp, Psi, SSQ, MU, S, Svₑ, Svₕ, Mvₑ, Mvₕ, vₑ₁, vₕ₁, θ;
                                            use_nlsolve = use_nlsolve, method = method, my_upwind = my_upwind,
                                            uniform = false),
              vcat(vₑ, vₕ), method = :anderson, m = m_anderson, ftol = tol)

N = length(η)
vₑ = out.zero[1:N]
vₕ = out.zero[(N + 1):end]
for ind in 2:(N - 1)
    Mvₑfd[ind] = ((vₑ[ind + 1] - vₑ[ind]) / dη[ind] * max(0, MU[ind]) +
                  (vₑ[ind] - vₑ[ind - 1]) / dη[ind - 1] * min(0, MU[ind])) / vₑ[ind] +
                  (vₑ[ind + 1] - 2 * vₑ[ind] + vₑ[ind - 1]) / (dη[ind] * dη[ind - 1]) / 2 / vₑ[ind] * S[ind]^2;
    Mvₕfd[ind] = ((vₕ[ind + 1] - vₕ[ind]) / dη[ind] * max(0, MU[ind]) +
                  (vₕ[ind] - vₕ[ind - 1]) / dη[ind - 1] * min(0, MU[ind])) / vₕ[ind] +
                  (vₕ[ind + 1] - 2 * vₕ[ind] + vₕ[ind - 1]) / (dη[ind] * dη[ind - 1]) / 2 / vₕ[ind] * S[ind]^2;
end
Mvₑfd[1] = ((vₑ[1 + 1] - vₑ[1]) / dη[1] * MU[1]) / vₑ[1]
Mvₕfd[1] = ((vₕ[1 + 1] - vₕ[1]) / dη[1] * MU[1]) / vₕ[1]
Mvₑfd[end] = (vₑ[end] - vₑ[end - 1]) / dη[end] * MU[end] / vₑ[end]
Mvₕfd[end] = (vₕ[end] - vₕ[end - 1]) / dη[end] * MU[end] / vₕ[end]

GS1 = (η        ./ vₑ ) .^ (1. / γ)
GS2 = ((1 .- η) ./ vₕ ) .^ (1. / γ)
CN  = GS1 .* Q .^ (1. / γ - 1.) ./ η # C/N*η
C_N = GS2 .* Q.^(1. / γ - 1.) ./ (1. .- η)
Mvₑ .= ρₑ .- CN  - (1 - γ) .* (Φ .- δ + Svₑ .* σ .- γ * σ^2 / 2.)
Mvₕ .= ρₕ .- C_N - (1 - γ) .* (Φ .- δ + Svₕ .* σ .- γ * σ^2 / 2.)
@show norm(vcat(Mvₑfd - Mvₑ, Mvₕfd - Mvₕ), pnorm)

end_time = time_ns()

elasped_time = (end_time - start_time) / 1e9

@show elasped_time

nothing


text:
export IntervalLibrary, Interval

"""
    IntervalLibrary{T}

Default library for polyhedra of dimension 1. Many aspect of polyhedral computation become trivial in one dimension. This library exploits this fact.
The library is also used as a fallback for libraries that do not support 1-dimensional polyhedra (e.g. qhull). That is projecting a polyhedron using such library produces a polyhedron using `IntervalLibrary`.
"""
struct IntervalLibrary{T} <: Library
end

similar_library(lib::IntervalLibrary, d::FullDim, ::Type{T}) where T = default_library(d, T) # default_library allows to fallback to DefaultLibrary if d is not FullDim(1)

mutable struct Interval{T, AT, D} <: Polyhedron{T}
    hrep::Intersection{T, AT, D}
    vrep::Hull{T, AT, D}
    length::T
end

Interval{T, AT, D}(d::FullDim, it::HIt...) where {T, AT, D} = _hinterval(Intersection{T, AT, D}(d, it...), AT, D)
Interval{T, AT, D}(d::FullDim, it::VIt...) where {T, AT, D} = _vinterval(Hull{T, AT, D}(d, it...), AT, D)

# If AT is an SVector, it will be StaticArrays.Size{(1,)}
# otherwise, it will be 1
FullDim(p::Interval) = FullDim(p.hrep)

library(::Union{Interval{T}, Type{<:Interval{T}}}) where T = IntervalLibrary{T}()

hvectortype(::Type{<:Interval{T, AT}}) where {T, AT} = AT
vvectortype(::Type{<:Interval{T, AT}}) where {T, AT} = AT

similar_type(::Type{<:Interval}, d::FullDim, ::Type{T}) where T = default_type(d, T)

surface(::Interval{T}) where {T} = zero(T)
volume(p::Interval) = p.length
Base.isempty(p::Interval) = isempty(p.vrep)

function _interval(AT, haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Bool) where {T}
    if haslb && hasub && _gt(lb, ub)
        isempty = true
    end
    hps = HyperPlane{T, AT}[]
    hss = HalfSpace{T, AT}[]
    ps = AT[]
    ls = Line{T, AT}[]
    rs = Ray{T, AT}[]
    if !isempty
        if hasub
            push!(ps, StaticArrays.SVector(ub))
            if haslb && _isapprox(lb, ub)
                push!(hps, HyperPlane(StaticArrays.SVector(one(T)), ub))
            else
                push!(hss, HalfSpace(StaticArrays.SVector(one(T)), ub))
            end
            if !haslb
                push!(rs, Ray(StaticArrays.SVector(-one(T))))
            end
        else
            if haslb
                push!(rs, Ray(StaticArrays.SVector(one(T))))
            else
                push!(ps, origin(AT, 1))
                push!(ls, Line(StaticArrays.SVector(one(T))))
            end
        end
        if haslb
            if !_isapprox(lb, ub)
                push!(hss, HalfSpace(StaticArrays.SVector(-one(T)), -lb))
                push!(ps, StaticArrays.SVector(lb))
            end
        end
    else
        # The dimension should be -1 so 1 - nhyperplanes == -1 so nhyperplanes == 2
        push!(hps, HyperPlane(StaticArrays.SVector(one(T)), zero(T)))
        push!(hps, HyperPlane(StaticArrays.SVector(zero(T)), one(T)))
    end
    h = hrep(hps, hss)
    v = vrep(ps, ls, rs)
    volume = isempty ? zero(T) : (haslb && hasub ? max(zero(T), ub - lb) : -one(T))
    return h, v, volume
end

function Interval{T, AT, D}(haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Bool) where {T, AT, D}
    return Interval{T, AT, D}(_interval(AT, haslb, lb, hasub, ub, isempty)...)
end

function _hinterval(rep::HRep{T}, ::Type{AT}) where {T, AT}
    haslb = false
    lb = zero(T)
    hasub = false
    ub = zero(T)
    empty = false
    function _setlb(newlb)
        if !haslb
            haslb = true
            lb = T(newlb)
        else
            lb = T(max(lb, newlb))
        end
    end
    function _setub(newub)
        if !hasub
            hasub = true
            ub = T(newub)
        else
            ub = T(min(ub, newub))
        end
    end
    for hp in hyperplanes(rep)
        α = hp.a[1]
        if isapproxzero(α)
            if !isapproxzero(hp.β)
                empty = true
            end
        else
            _setlb(hp.β / α)
            _setub(hp.β / α)
        end
    end
    for hs in halfspaces(rep)
        α = hs.a[1]
        if isapproxzero(α)
            if hs.β < 0
                empty = true
            end
        elseif α < 0
            _setlb(hs.β / α)
        else
            _setub(hs.β / α)
        end
    end
    return _interval(AT, haslb, lb, hasub, ub, empty)
end

function _hinterval(rep::HRep{T}, ::Type{AT}, D) where {T, AT}
    return Interval{T, AT, D}(_hinterval(rep, AT)...)
end

function _vinterval(v::VRep{T}, ::Type{AT}, D) where {T, AT}
    haslb = true
    lb = zero(T)
    hasub = true
    ub = zero(T)
    isempty = true
    for p in points(v)
        x = coord(p)[1]
        if isempty
            isempty = false
            lb = x
            ub = x
        else
            lb = min(lb, x)
            ub = max(ub, x)
        end
    end
    for l in lines(v)
        if !isapproxzero(l)
            isempty = false
            haslb = false
            hasub = false
        end
    end
    for r in rays(v)
        x = coord(r)[1]
        if !isapproxzero(x)
            isempty = false
            if x > 0
                hasub = false
            else
                haslb = false
            end
        end
    end
    Interval{T, AT, D}(haslb, lb, hasub, ub, isempty)
end

Interval{T, AT, D}(p::HRepresentation{T}) where {T, AT, D} = _hinterval(p, AT, D)
Interval{T, AT, D}(p::VRepresentation{T}) where {T, AT, D} = _vinterval(p, AT, D)
function Interval{T, AT, D}(p::Polyhedron{T}) where {T, AT, D}
    if hrepiscomputed(p)
        _hinterval(p, AT, D)
    else
        _vinterval(p, AT, D)
    end
end

function polyhedron(rep::Rep{T}, ::IntervalLibrary{T}) where T
    Interval{T, StaticArrays.SVector{1, T}, StaticArrays.Size{(1,)}}(rep)
end

hrep(p::Interval) = p.hrep
vrep(p::Interval) = p.vrep

hrepiscomputed(::Interval) = true
vrepiscomputed(::Interval) = true

# Nothing to do
function detecthlinearity!(::Interval) end
function detectvlinearity!(::Interval) end
function removehredundancy!(::Interval) end
function removevredundancy!(::Interval) end

function sethrep!(p::Interval{T, AT}, h::HRep) where {T, AT}
    hnew, v, volume = _hinterval(h, AT)
    p.hrep = hnew
    p.vrep = v
    p.length = volume
    return p
end


text:
# Model
include_model("biped")
model = free_time_model(model)

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

# Horizon
T = 51
Tm = 26

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

# Configurations
# 1: x pos
# 2: z pos
# 3: torso angle (rel. to upward vertical)
# 4: thigh 1 angle (rel. to downward vertical)
# 5: calf 1 (rel. to downward vertical)
# 6: thigh 2 (rel. to downward vertical)
# 7: calf 2 (rel. to downward vertical)

q1, qM, qT = initial_configuration(model, -pi / 100.0, pi / 10.0, -pi / 25.0, -pi / 25.0)
visualize!(vis, model, [q1], Δt = h)

# feet positions
foot_2_1 = kinematics_2(model, q1, body = :calf_2, mode = :ee)
foot_2_M = kinematics_2(model, qM, body = :calf_2, mode = :ee)

foot_1_M = kinematics_2(model, qM, body = :calf_1, mode = :ee)
foot_1_T = kinematics_2(model, qT, body = :calf_1, mode = :ee)

# feet trajectories
zh = 0.1
foot_2_x = [range(foot_2_1[1], stop = foot_2_M[1], length = Tm)...,
    [foot_2_M[1] for t = 1:Tm-1]...]
foot_2_z = sqrt.((zh^2.0) * (1.0 .- ((foot_2_x).^2.0)
    ./ abs(foot_2_1[1])^2.0) .+ 1.0e-8)

foot_1_x = [[foot_1_M[1] for t = 1:Tm-1]...,
    range(foot_1_M[1], stop = foot_1_T[1], length = Tm)...]
foot_1_z = sqrt.((zh^2.0) * (1.0 .- ((foot_1_x .- abs(foot_1_T[1] / 2.0)).^2.0)
    ./ abs(foot_1_T[1] / 2.0)^2.0) .+ 1.0e-8)

using Plots
plot(foot_2_x, foot_2_z)
plot!(foot_1_x, foot_1_z, aspect_ratio = :equal)

t = range(0.0, stop = tf, length = T)
plot(t, foot_2_x)
plot!(t, foot_1_x)

plot(t, foot_2_z)
plot!(t, foot_1_z)

# Bounds

# control
# u = (τ1..4, λ1..2, β1..4, ψ1..2, η1...4, s1)
# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[end] = 2.0 * h
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[end] = 0.25 * h
ul, uu = control_bounds(model, T, _ul, _uu)

# state
xl, xu = state_bounds(model, T,
    -Inf * ones(model.n), Inf * ones(model.n),
    x1 = [Inf * ones(model.nq); q1],
    xT = [Inf * ones(model.nq); qT])

xl[Tm][model.nq .+ (1:model.nq)] = copy(qM)
xu[Tm][model.nq .+ (1:model.nq)] = copy(qM)

for t = 1:T
    xl[t][3] = q1[3] - pi / 20.0
    xl[t][10] = q1[3] - pi / 20.0
    xu[t][3] = q1[3] + pi / 20.0
    xu[t][10] = q1[3] + pi / 20.0
end

# Objective
include_objective(["velocity", "nonlinear_stage"])

# q_ref = linear_interpolation(q1, qM, T)
q_ref = [linear_interpolation(q1, qM, Tm)...,
    linear_interpolation(qM, qT, Tm)[2:end]...]
x0 = configuration_to_state(q_ref)

vis = Visualizer()
render(vis)
visualize!(vis, model, q_ref, Δt = h)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e5, model.m - 1)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)

obj_control = quadratic_time_tracking_objective(
    [Diagonal(1.0e-1 * ones(model.n)) for t = 1:T],
    [Diagonal([1.0e-3 * ones(model.nu)..., 0.0 * ones(model.m - model.nu - 1)..., 0.0]) for t = 1:T-1],
    [x0[t] for t = 1:T],
    [zeros(model.m) for t = 1:T-1],
    1.0)

# quadratic velocity penalty
# Σ v' Q v
q_v = 5.0 * ones(model.nq)
q_v[3] = 100.0
obj_velocity = velocity_objective(
    [Diagonal(q_v) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3, 4, 5, 6, 7]))

# torso height
t_h = kinematics_1(model, q1, body = :torso, mode = :com)[2]
function l_stage_torso_h(x, u, t)
    return 10000.0 * (kinematics_1(model,
            view(x, 8:14), body = :torso, mode = :com)[2] - t_h)^2.0
end

l_terminal_torso_h(x) = 0.0
obj_th = nonlinear_stage_objective(l_stage_torso_h, l_terminal_torso_h)

# torso lateral
function l_stage_torso_lat(x, u, t)
    return (1.0 * (kinematics_1(model,
        view(x, 8:14), body = :torso, mode = :com)[1]
        - kinematics_1(model,
        view(x0[t], 8:14), body = :torso, mode = :com)[1])^2.0)
end
l_terminal_torso_lat(x) = 0.0
obj_tl = nonlinear_stage_objective(l_stage_torso_lat, l_terminal_torso_lat)

# foot 1 height
function l_stage_fh1(x, u, t)
    return (10000.0 * (kinematics_2(model,
        view(x, 1:7), body = :calf_1, mode = :ee)[2] - foot_1_z[t])^2.0
        + 10000.0 * (kinematics_2(model,
            view(x, 8:14), body = :calf_1, mode = :ee)[2] - foot_1_z[t])^2.0)
end
l_terminal_fh1(x) = 0.0
obj_fh1 = nonlinear_stage_objective(l_stage_fh1, l_terminal_fh1)

# foot 1 lateral
function l_stage_fl1(x, u, t)
    return (1000.0 * (kinematics_2(model,
        view(x, 1:7), body = :calf_1, mode = :ee)[1] - foot_1_x[t])^2.0
        + 1000.0 * (kinematics_2(model,
            view(x, 8:14), body = :calf_1, mode = :ee)[1] - foot_1_x[t])^2.0)
end
l_terminal_fl1(x) = 0.0
obj_fl1 = nonlinear_stage_objective(l_stage_fl1, l_terminal_fl1)

# foot 2 height
function l_stage_fh2(x, u, t)
    return (10000.0 * (kinematics_2(model,
        view(x, 1:7), body = :calf_2, mode = :ee)[2] - foot_2_z[t])^2.0
        + 10000.0 * (kinematics_2(model,
            view(x, 8:14), body = :calf_2, mode = :ee)[2] - foot_2_z[t])^2.0)
end
l_terminal_fh2(x) = 0.0
obj_fh2 = nonlinear_stage_objective(l_stage_fh2, l_terminal_fh2)

# foot 2 lateral
function l_stage_fl2(x, u, t)
    return (1000.0 * (kinematics_2(model,
        view(x, 1:7), body = :calf_2, mode = :ee)[2] - foot_2_x[t])^2.0
        + 1000.0 * (kinematics_2(model,
            view(x, 8:14), body = :calf_2, mode = :ee)[2] - foot_2_x[t])^2.0)
end
l_terminal_fl2(x) = 0.0
obj_fl2 = nonlinear_stage_objective(l_stage_fl2, l_terminal_fl2)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
                      obj_th,
                      obj_tl,
                      obj_fh1,
                      obj_fh2,
                      obj_fl1,
                      obj_fl2])

# Constraints
include_constraints(["contact", "loop", "free_time"])
con_loop = loop_constraints(model, collect([(2:7)...,(9:14)...]), 1, T)
con_contact = contact_constraints(model, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_contact, con_free_time, con_loop])

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

# trajectory initialization
u0 = [[1.0e-5 * rand(model.m - 1); h] for t = 1:T-1] # random controls

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

# Solve

if true
    include_snopt()

	@time z̄ , info = solve(prob, copy(z0),
		nlp = :SNOPT7,
		tol = 1.0e-5, c_tol = 1.0e-5, mapl = 5,
		time_limit = 60 * 3)
	@show check_slack(z̄, prob)
	x̄, ū = unpack(z̄, prob)
    tfc, tc, h̄ = get_time(ū)

# 	# projection
# 	Q = [Diagonal(ones(model.n)) for t = 1:T]
# 	R = [Diagonal(0.1 * ones(model.m)) for t = 1:T-1]
# 	x_proj, u_proj = lqr_projection(model, x̄, ū, h̄[1], Q, R)
#
# 	@show tfc
# 	@show h̄[1]
# 	@save joinpath(pwd(), "examples/trajectories/biped_gait.jld2") x̄ ū h̄ x_proj u_proj
# else
# 	@load joinpath(pwd(), "examples/trajectories/biped_gait.jld2") x̄ ū h̄ x_proj u_proj
end

# Visualize
vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h̄[1])
q_vis = state_to_configuration(x̄)
len_stride = x̄[end][1] - x̄[1][1]
for i = 1:4
	for t = 1:T-1
		_q = Array(copy(state_to_configuration(x̄)[t+1]))
		_q[1] += i * len_stride
		push!(q_vis, _q)
	end
end
visualize!(vis, model, q_vis, Δt = h̄[1])

# projection
Q = [Diagonal(ones(model.n)) for t = 1:T]
R = [Diagonal(0.1 * ones(model.m)) for t = 1:T-1]
x_proj, u_proj = lqr_projection(model, x̄, ū, h̄[1], Q, R)

# Feet trajectories
fh1 = [kinematics_2(model,
    state_to_configuration(x̄)[t], body = :calf_1, mode = :ee)[2] for t = 1:T]
fh2 = [kinematics_2(model,
    state_to_configuration(x̄)[t], body = :calf_2, mode = :ee)[2] for t = 1:T]

plot(fh1, linetype = :steppost, label = "foot 1")
plot!(fh2, linetype = :steppost, label = "foot 2")

# Controls
plot(hcat(ū...)[1:4, :]',
    linetype = :steppost,
    label = "",
    color = :red,
    width = 2.0)

plot!(hcat(u_proj...)[1:4, :]',
    linetype = :steppost,
    label = "",
    color = :black)

# plot(hcat(u_proj...)[5:6, :]',
#     linetype = :steppost,
#     label = "",
#     width = 2.0)

# States
plot(hcat(state_to_configuration(x̄)...)',
    color = :red,
    width = 2.0,
    label = "")

plot!(hcat(state_to_configuration(x_proj)...)',
    color = :black,
    label = "")


text:
using TikzPictures
using ParticleFilters

function POMDPModelTools.render(m::LaserTagPOMDP, step)
    return LaserTagVis(m,
                       get(step, :a, nothing),
                       get(step, :r, nothing),
                       get(step, :sp, nothing),
                       get(step, :o, nothing),
                       get(step, :bp, nothing))
end

mutable struct LaserTagVis
    p::LaserTagPOMDP
    a::Union{Any, Nothing}
    r::Union{Any, Nothing}
    s::Union{Any, Nothing}
    o::Union{Any, Nothing}
    b::Union{Any, Nothing}
end
LaserTagVis(p; s=nothing, a=nothing, o=nothing, b=nothing, r=nothing) = LaserTagVis(p, a, r, s, o, b)
LaserTagVis(p::LaserTagPOMDP, arspobp::Tuple) = LaserTagVis(p, arspobp...)

Base.show(io::IO, mime::MIME"image/svg+xml", v::LaserTagVis) = show(io, mime, tikz_pic(v))
Base.show(io::IO, mime::MIME"text/html", v::LaserTagVis) = show(io, MIME("image/svg+xml"), v)

function Base.show(io::IO, mime::MIME"image/png", v::LaserTagVis)
    # fname = "/home/fischer" * tempname()
    fname = pwd() * "/.." * tempname()
    save(PDF(fname), tikz_pic(v))
    run(`gm convert -flatten $(fname*".pdf") $(fname*".png")`)
    open(fname*".png") do f
        write(io, f)
    end
    run(`rm $(fname*".pdf")`)
    run(`rm $(fname*".png")`)
    return io
end

function fill_square(o::IO, x, y, color, opacity=0.5) # maybe opacity should be a keyword
    sqsize = 1.0
    println(o, "\\fill[$(color), opacity=$opacity] ($((x-1) * sqsize),$((y-1) * sqsize)) rectangle +($sqsize,$sqsize);")
end

function show_belief(o::IO, b::ParticleCollection{LTState})
    d = Dict{Coord, Int}()
    total = 0
    for p in particles(b)
        if !p.terminal
            opp = p.opponent
            if haskey(d, opp)
                d[opp] += 1
            else
                d[opp] = 1
            end
            total += 1
        end
    end
    for (opp, freq) in d
        fill_square(o, opp[1], opp[2], "yellow", sqrt(freq/total))
    end
end

function show_meas(o::IO, s::LTState, obs::Union{CMeas,DMeas})
    middle = s.robot - 0.5
    for i in 1:4
        dir = CARDINALS[i]
        start = middle+0.5*dir
        finish = start + obs[i]*dir
        draw_laser(o, start, finish)
    end
    for i in 5:8
        dir = DIAGONALS[i-4]*sqrt(2)/2
        start = middle+0.5*sqrt(2)*dir
        finish = start + obs[i]*dir
        draw_laser(o, start, finish)
    end
end

function draw_laser(o::IO, start::AbstractVector{Float64}, finish::AbstractVector{Float64})
    println(o, "\\draw[dashed, red] $(@sprintf("(%0.3f, %0.3f)", start[1], start[2])) -- $(@sprintf("(%0.3f, %0.3f)", finish[1], finish[2]));")
end

function tikz_pic(v::LaserTagVis)
    p = v.p
    f = p.floor
    o = IOBuffer()
    sqsize=1

    println(o, "\\begin{scope}")

    println(o, "\\clip (0,0) rectangle ($(f.n_cols*sqsize),$(f.n_rows*sqsize));")

    for c in p.obstacles
        fill_square(o, c[1], c[2], "gray")
    end

    if v.b !== nothing
        show_belief(o, v.b)
    end

    if v.s !== nothing
        s = v.s
        opp = s.opponent
        rob = s.robot
        fill_square(o, opp[1], opp[2], "orange")
        fill_square(o, rob[1], rob[2], "green")
        if v.o !== nothing
            show_meas(o, s, v.o)
        end
        if v.a !== nothing
            aname = ACTION_NAMES[v.a]
            println(o, "\\node[above right] at ($((rob[1]-1) * sqsize), $((rob[2]-1) * sqsize)) {$aname};")
        end
        if v.r !== nothing
            rtext = @sprintf("%0.2f", v.r)
            println(o, "\\node[below right] at ($((rob[1]-1) * sqsize), $((rob[2]-1) * sqsize)) {$rtext};")
        end

    end

    # # possibly for later: text in a square
    # vs = @sprintf("%0.2f", V[i])
    # println(o, "\\node[above right] at ($((xval-1) * sqsize), $((yval-1) * sqsize)) {\$$(vs)\$};")

    println(o, "\\draw[black] grid($(f.n_cols), $(f.n_rows));")

    println(o, "\\end{scope}")

    tikzDeleteIntermediate(true)
    return TikzPicture(String(take!(o)), options="scale=1.25")
end

function Base.show(io::IO, mime::MIME"text/plain", v::LaserTagVis)
    for y in n_rows(v.p):-1:1
        for x in 1:n_cols(v.p)
            printed = false
            if v.s !== nothing
                s = v.s
                if Coord(x,y) == s.robot
                    print(io, 'R')
                    printed = true
                elseif Coord(x,y) == s.opponent
                    print(io, 'O')
                    printed = true
                end
            end
            if !printed
                if Coord(x,y) in v.p.obstacles
                    print(io, '#')
                else
                    print(io, '.')
                end
            end
        end
        print(io, '\n')
    end
end


text:
mutable struct Sampler
    magFilter::Union{Nothing,Int}
    minFilter::Union{Nothing,Int}
    wrapS::Union{Nothing,Int}
    wrapT::Union{Nothing,Int}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Sampler()
        obj = new()
        obj.magFilter = nothing
        obj.minFilter = nothing
        obj.wrapS = nothing
        obj.wrapT = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::Sampler, sym::Symbol)
    x = getfield(obj, sym)
    sym === :wrapS && x === nothing && return 10497
    sym === :wrapT && x === nothing && return 10497
    return x
end

function Base.setproperty!(obj::Sampler, sym::Symbol, x)
    if sym === :magFilter && x !== nothing
        x === NEAREST || x === LINEAR ||
            throw(ArgumentError("magFilter should be one of: NEAREST(9728), LINEAR(9729)"))
    elseif sym === :minFilter && x !== nothing
        x === NEAREST || x === LINEAR || x === NEAREST_MIPMAP_NEAREST ||
        x === LINEAR_MIPMAP_NEAREST || x === NEAREST_MIPMAP_LINEAR || x === LINEAR_MIPMAP_LINEAR ||
            throw(ArgumentError("minFilter should be one of: NEAREST(9728), LINEAR(9729), NEAREST_MIPMAP_NEAREST(9984), LINEAR_MIPMAP_NEAREST(9985), NEAREST_MIPMAP_LINEAR(9986), LINEAR_MIPMAP_LINEAR(9987)"))
    elseif sym === :wrapS && x !== nothing
        x === CLAMP_TO_EDGE || x === MIRRORED_REPEAT || x === REPEAT ||
            throw(ArgumentError("wrapS should be one of: CLAMP_TO_EDGE(33071), MIRRORED_REPEAT(33648), REPEAT(10497)"))
    elseif sym === :wrapT && x !== nothing
        x === CLAMP_TO_EDGE || x === MIRRORED_REPEAT || x === REPEAT ||
            throw(ArgumentError("wrapT should be one of: CLAMP_TO_EDGE(33071), MIRRORED_REPEAT(33648), REPEAT(10497)"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Sampler}) = JSON3.Mutable()
JSON3.omitempties(::Type{Sampler}) = (:magFilter, :minFilter, :wrapS, :wrapT, :name, :extensions, :extras)


text:
using LowRank, Test
# using JLD
# using Random, Distributions, LinearAlgebra

N = 50
d = 5
n = 10
k = 2
error = 0.0
indices = 1:n

A = zeros(d, N)
U = rand(d, k)
V = rand(k, n)
A[:, 1:n] = U*V
A[:, n+1:N] = rand(d, N-n)
data_set = collect(eachcol(A))

# file = load(joinpath(@__DIR__, "..","datasets", "N_50_d_5.jld"))
# data_set, (err, indices) = file["data"], file["n_10_k_2"]
# n = 10
# k = 2

# @time @test find_low_rank_subset_checkall(data_set,n,k) == (err, indices)
@test find_low_rank_subset_sample_rep(data_set,n,k)[1] <= 4*error + 0.0001
@test FastLowRank(data_set, n, k)[1] <= error + 0.0001


text:
#
# The functions and data structures in this file organize the
# iteration statistics and report the results after the iteration
# is complete
#

#
# Keep the books for nsol and secant
#
mutable struct ItStats{T<:Real}
    ifun::Array{Int64,1}
    ijac::Array{Int64,1}
    iarm::Array{Int64,1}
    history::Array{T,1}
end

#
# initfun = 1 unless it's the scalar secant method
#             then it's 2
#
function ItStats(hist, initfun = 1)
    ItStats([initfun], [0], [0], [hist])
end

function updateStats!(ItData::ItStats, newfun, newjac, AOUT)
    newiarm = AOUT.aiarm
    newfun = newfun + newiarm + 1
    resnorm = AOUT.resnorm
    append!(ItData.ifun, newfun)
    append!(ItData.ijac, newjac)
    append!(ItData.iarm, newiarm)
    append!(ItData.history, resnorm)
end

function CollectStats(ItData::ItStats)
    stats = (ifun = ItData.ifun, ijac = ItData.ijac, iarm = ItData.iarm)
    return stats
end

#
# Keep the books for nsoli
#
mutable struct ItStatsK{T<:Real}
    ifun::Array{Int64,1}
    ijac::Array{Int64,1}
    iarm::Array{Int64,1}
    ikfail::Array{Int64,1}
    history::Array{T,1}
end

function ItStatsK(hist)
    ItStatsK([1], [0], [0], [0], [hist])
end

function updateStats!(ItData::ItStatsK, newfun, newjac, AOUT, newikfail)
    newiarm = AOUT.aiarm
    newfun = newfun + newiarm + 1
    resnorm = AOUT.resnorm
    append!(ItData.ifun, newfun)
    append!(ItData.ijac, newjac)
    append!(ItData.iarm, newiarm)
    append!(ItData.ikfail, newikfail)
    append!(ItData.history, resnorm)
end

function CollectStats(ItData::ItStatsK)
    stats =
        (ifun = ItData.ifun, ijac = ItData.ijac, iarm = ItData.iarm, ikfail = ItData.ikfail)
    return stats
end

#
# Keep stats for PTC
#
mutable struct ItStatsPTC{T<:Real}
    history::Array{T,1}
end

function ItStatsPTC(hist)
    ItStatsPTC([hist])
end

function updateStats!(ItData::ItStatsPTC, resnorm)
    append!(ItData.history, resnorm)
end

function CollectStats(ItData::ItStatsPTC)
    stats = []
    return stats
end

#
# Keep the books for PTC-Krylov
#
mutable struct ItStatsPTCK{T<:Real}
    ifun::Array{Int64,1}
    ijac::Array{Int64,1}
    ikfail::Array{Int64,1}
    history::Array{T,1}
end

function ItStatsPTCK(hist)
    ItStatsPTCK([1], [0], [0], [hist])
end

function updateStats!(ItData::ItStatsPTCK, resnorm, newjac, newikfail)
    append!(ItData.history, resnorm)
    append!(ItData.ifun, 1)
    append!(ItData.ijac, newjac)
    append!(ItData.ikfail, newikfail)
end

function CollectStats(ItData::ItStatsPTCK)
    stats = (ifun = ItData.ifun, ijac = ItData.ijac, ikfail = ItData.ikfail)
    return stats
end

"""
CloseIteration(x, FS, ItData, idid, errcode, keepsolhist, solhist = [])

Collect the solution, function value, iteration stats and send them back.
"""
function CloseIteration(x, FS, ItData, idid, errcode, keepsolhist, solhist = [])
    stats = CollectStats(ItData)
    ithist = ItData.history
    if keepsolhist
        sizehist = length(ithist)
        return (
            solution = x,
            functionval = FS,
            history = ithist,
            stats = stats,
            idid = idid,
            errcode = errcode,
            solhist = solhist[:, 1:sizehist],
        )
    else
        return (
            solution = x,
            functionval = FS,
            history = ithist,
            stats = stats,
            idid = idid,
            errcode = errcode,
        )
    end
end


text:
# Data Crunching Debugger

module Dcd

	type Session
		enable::Bool
		iters::Array{Dict{String,Any}, 1}

		Session(enable = false) = ( x = new(); x.enable = enable; x.iters = Array(Dict{String,Any}, 0); return x; )
	end

	macro set(sess, key, val) end
	macro it(sess) end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
# ------------------------------------------------------------------------------
# Schwarz smoother example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh1D(1.0)
p = 2

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)
# Schwarz smoother
schwarz = Schwarz(diffusion)

# compute operator symbols
A = computesymbols(schwarz, [], [π])
eigenvalues = real(eigvals(A))
#println(eigenvalues)

# ------------------------------------------------------------------------------


text:
using SymPy, Plots
include("../SFFM.jl")

tempCMEParams = Dict()
open("src/CMEParamsData/iltcme.json", "r") do f
    global tempCMEParams
    tempCMEParams=JSON.parse(f)  # parse and transform data
end

s = SymPy.Sym("s")
c = SymPy.Sym("c")
a = SymPy.Sym("a")
b = SymPy.Sym("b")
k = SymPy.Sym("k")
omega = SymPy.Sym("omega")

term0 = c/(1+s)
termk = (a*(1+s) + b*k*omega) / ((1+s)^2 + (k*omega)^2)

d1term0 = -SymPy.diff(term0,s,1)
d1termk = -SymPy.diff(termk,s,1)
d5term0 = -SymPy.diff(term0,s,5)
d5termk = -SymPy.diff(termk,s,5)

keysArray = []
m0Array = []
m1Array = []
m5Array = []
for key in keys(SFFM.CMEParams)
    if key < 200
        print(key); print(",")
        cVal = SFFM.CMEParams[key]["c"]
        omegaVal = SFFM.CMEParams[key]["omega"]
        aVec = SFFM.CMEParams[key]["a"]
        bVec = SFFM.CMEParams[key]["b"]
        n = length(aVec)

        m0 = N(term0(s=>0, c=>cVal, omega=>omegaVal))
        m1 = N(d1term0(s=>0, c=>cVal, omega=>omegaVal))
        m5 = N(d5term0(s=>0, c=>cVal, omega=>omegaVal))
        push!(keysArray, key)
        for kVal in 1:n
            m0 += N(termk(
                s=>0,
                c=>cVal, omega=>omegaVal,
                a=>aVec[kVal], b=>bVec[kVal],
                k=>kVal
            ))
            m1 += N(d1termk(
                s=>0,
                c=>cVal, omega=>omegaVal,
                a=>aVec[kVal], b=>bVec[kVal],
                k=>kVal
            ))
            m5 += N(d5termk(
                s=>0,
                c=>cVal, omega=>omegaVal,
                a=>aVec[kVal], b=>bVec[kVal],
                k=>kVal
            ))
        end
        push!(m0Array, m0)
        push!(m1Array, m1)
        push!(m5Array, m5)
    end
end

display(scatter(keysArray, m0Array))
display(scatter(keysArray, m1Array))
display(scatter(keysArray, log.((m5Array).^(1/5)-m1Array)))

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM