text:
"""
Parameters of the system.

# Arguments:

- `stereo::Bool`: Set to `false` if in monocular mode, otherwise `true`.
    Default is `false`.
- `max_nb_keypoints::Int64`: Maximum number of keypoints to detect in the image.
    Default is `1000`.
- `max_distance::Int64`: Cell size in pixels of the grid in Frame.
    Each frame is divided into a grid, where each cell contains certain
    number of keypoints. This grid is used when retrieving neighbouring
    keypoints from the Frame.
- `max_ktl_distance::Float64`: When doing optical flow tracking, it is done
    by first tracking keypoint from previous frame to the current frame and
    then in reverse. If the distance between keypoint in previous frame and
    reverse tracked is greater than this distance, then that keypoint is
    considered lost and discarded. Default is `1.0`.
- `pyramid_levels::Int64`: Number of pyramid levels to use during optical flow.
    Set to `0` to use only original image. Otherwise if set to `L`, then total
    number of levels is `L + 1`. Default is `3`.
- `pyramid_σ::Float64`: `σ` for the Gaussian blur, used when constructing
    pyramids. Default is `1.0`.
- `window_size::Int64`: Size of the window to use when doing optical flow
    tracking. Default is `9`.
- `initial_parallax::Float64`: Average amount of the parallax needed
    between the frames for the initialization. Default is `20.0`.
- `max_reprojection_error::Float64`: Maximum reprojection error for the
    mappoint to be considered inlier. It is used in triangulation
    and pose calculation. Default is `3.0`.
- `min_cov_score::Int64`: Minimum number of 3D MapPoints required
    to perform Bundle-Adjustment. Default is `25`.
- `filtering_ratio::Float64`:
    Number of good keypoints divided by the total number of keypoints
    in a Frame. This is used to filter out non-informative Frames.
    If frame ratio is greater than this value, then it is removed during
    map filtering. Default is `0.9`.
- `map_filtering::Bool`: Set to `true` to perform local map matching.
    The goal of this is to try to re-match lost 3D keypoints back into
    the current frame.
- `max_projection_distance::Float64`: Maximum distance between projected
    and target pixels to consider keypoints being "the same".
    Is is used during local map matching. Default is `2.0`.
- `max_descriptor_distance::Float64`: Maximum distance ratio between
    descriptors to be considered similar. It is used during local map matching.
    Default is `0.35`. The lesser the value, the more accurate matching will be.
    Should be in `(0, 1)` range.

**Below variables are part of the SLAM state and should not change manually.**

- `vision_initialized::Bool`: If `true`, then the Front-End was successfully
    initialized.
- `reset_required::Bool`: Indicates if the visual part requires reset.
    This could happen if the system is tracking too little points.
- `local_ba_on::Bool`: `true`, when Estimator is performing
    Local Bundle Adjustment.
"""
Base.@kwdef mutable struct Params
    stereo::Bool = false
    max_nb_keypoints::Int64 = 1000
    max_distance::Int64 = 35
    max_ktl_distance::Float64 = 1.0
    pyramid_levels::Int64 = 3
    pyramid_σ::Float64 = 1.0
    window_size::Int64 = 9
    initial_parallax::Float64 = 20.0
    max_reprojection_error::Float64 = 3.0
    min_cov_score::Int64 = 25
    do_local_matching::Bool = false

    filtering_ratio::Float64 = 0.9
    map_filtering::Bool = false

    do_local_bundle_adjustment::Bool = true
    max_projection_distance::Float64 = 2.0
    max_descriptor_distance::Float64 = 0.35

    # State variables.
    vision_initialized::Bool = false
    reset_required::Bool = false
    local_ba_on::Bool = false
end

"""
```julia
reset!(p::Params)
```

Reset state of the system.
"""
function reset!(p::Params)
    p.vision_initialized = false
    p.reset_required = false
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# tsp.jl
#
# Solves the travelling salesman problem using integer programming and
# lazy generation of the subtour elimination constraints.
#############################################################################

using JuMP
using Gurobi
using Base.Test

# extractTour
# Given a n-by-n matrix representing the solution to an undirected TSP,
# extract the tour as a vector
# Input:
#  n        Number of cities
#  sol      n-by-n 0-1 symmetric matrix representing solution
# Output:
#  tour     n+1 length vector of tour, starting and ending at 1
function extractTour(n, sol)
    tour = [1]  # Start at city 1 always
    cur_city = 1
    while true
        # Look for first arc out of current city
        for j = 1:n
            if sol[cur_city,j] >= 1-1e-6
                # Found next city
                push!(tour, j)
                # Don't ever use this arc again
                sol[cur_city, j] = 0.0
                sol[j, cur_city] = 0.0
                # Move to next city
                cur_city = j
                break
            end
        end
        # If we have come back to 1, stop
        if cur_city == 1
            break
        end
    end  # end while
    return tour
end

# findSubtour
# Given a n-by-n matrix representing solution to the relaxed
# undirected TSP problem, find a set of nodes belonging to a subtour
# Input:
#  n        Number of cities
#  sol      n-by-n 0-1 symmetric matrix representing solution
# Outputs:
#  subtour  n length vector of booleans, true iff in a particular subtour
#  subtour_length   Number of cities in subtour (if n, no subtour found)
function findSubtour(n, sol)
    # Initialize to no subtour
    subtour = fill(false,n)
    # Always start looking at city 1
    cur_city = 1
    subtour[cur_city] = true
    subtour_length = 1
    while true
        # Find next node that we haven't yet visited
        found_city = false
        for j = 1:n
            if !subtour[j]
                if sol[cur_city, j] >= 1 - 1e-6
                    # Arc to unvisited city, follow it
                    cur_city = j
                    subtour[j] = true
                    found_city = true
                    subtour_length += 1
                    break  # Move on to next city
                end
            end
        end
        if !found_city
            # We are done
            break
        end
    end
    return subtour, subtour_length
end

# solveTSP
# Given a matrix of city locations, solve the TSP
# Inputs:
#   n       Number of cities
#   cities  n-by-2 matrix of (x,y) city locations
# Output:
#   path    Vector with order to cities are visited in
function solveTSP(n, cities)

    # Calculate pairwise distance matrix
    dist = zeros(n, n)
    for i = 1:n
        for j = i:n
            d = norm(cities[i,1:2] - cities[j,1:2])
            dist[i,j] = d
            dist[j,i] = d
        end
    end

    # Create a model that will use Gurobi to solve
    m = Model(solver=GurobiSolver())

    # x[i,j] is 1 iff we travel between i and j, 0 otherwise
    # Although we define all n^2 variables, we will only use
    # the upper triangle
    @variable(m, x[1:n,1:n], Bin)

    # Minimize length of tour
    @objective(m, Min, sum{dist[i,j]*x[i,j], i=1:n,j=i:n})

    # Make x_ij and x_ji be the same thing (undirectional)
    # Don't allow self-arcs
    for i = 1:n
        @constraint(m, x[i,i] == 0)
        for j = (i+1):n
            @constraint(m, x[i,j] == x[j,i])
        end
    end

    # We must enter and leave every city once and only once
    for i = 1:n
        @constraint(m, sum{x[i,j], j=1:n} == 2)
    end

    function subtour(cb)
        # Optional: display tour starting at city 1
        println("----\nInside subtour callback")
        println("Current tour starting at city 1:")
        print(extractTour(n, getvalue(x)))

        # Find any set of cities in a subtour
        subtour, subtour_length = findSubtour(n, getvalue(x))

        if subtour_length == n
            # This "subtour" is actually all cities, so we are done
            println("Solution visits all cities")
            println("----")
            return
        end

        # Subtour found - add lazy constraint
        # We will build it up piece-by-piece
        arcs_from_subtour = zero(AffExpr)

        for i = 1:n
            if !subtour[i]
                # If this city isn't in subtour, skip it
                continue
            end
            # Want to include all arcs from this city, which is in
            # the subtour, to all cities not in the subtour
            for j = 1:n
                if i == j
                    # Self-arc
                    continue
                elseif subtour[j]
                    # Both ends in same subtour
                    continue
                else
                    # j isn't in subtour
                    arcs_from_subtour += x[i,j]
                end
            end
        end

        # Add the new subtour elimination constraint we built
        println("Adding subtour elimination cut")
        println("----")
        @lazyconstraint(cb, arcs_from_subtour >= 2)
    end  # End function subtour

    # Solve the problem with our cut generator
    addlazycallback(m, subtour)
    solve(m)

    # Return best tour
    return extractTour(n, getvalue(x))
end  # end solveTSP

# Create a simple instance that looks like
#       +           +
#   +                   +
#       +           +
# The optimal tour is obvious, but the initial solution will be
#    /--+           +--\
#   +               |   +
#    \--+           +--/
n = 6
cities = [ 50 200;
                    100 100;
                    100 300;
                    500 100;
                    500 300;
                    550 200]
tour = solveTSP(n, cities)
println("Solution: ")
println(tour)

ext: jl


text:
using Genie
import Genie.Router: route
import Genie.Renderer.Json: json

Genie.config.run_as_server = true

route("/") do
  (:message => "Hi there!") |> json
end

Genie.startup()

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
struct Horizon
    num_hours::Integer   # Number of hours in horizon

    function Horizon(num_hours::Integer)
        @assert mod(num_hours,24) == 0 "Horizon should be an even number of days"
        return new(num_hours)
    end
end

Base.:(==)(h1::Horizon,h2::Horizon) = h1.num_hours == h2.num_hours
Base.:(<=)(h1::Horizon,h2::Horizon) = h1.num_hours <= h2.num_hours
Base.:(>=)(h1::Horizon,h2::Horizon) = h1.num_hours >= h2.num_hours

Day() = Horizon(24)
Days(num_days::Integer) = Horizon(num_days*24)
Week() = Horizon(168)
Weeks(num_weeks::Integer) = Horizon(num_weeks*168)
Year() = Horizon(8760)
Years(num_years::Integer) = Horizon(num_years*8760)

num_hours(horizon::Horizon) = horizon.num_hours
num_days(horizon::Horizon) = div(horizon.num_hours, 24)
num_weeks(horizon::Horizon) = div(horizon.num_hours, 168)

function horizonstring(horizon::Horizon)
    remaindays = div(mod(num_hours(horizon),168),24)
    horizonstr = ""
    weeks = num_weeks(horizon)
    if weeks > 0
        horizonstr *= string("(",weeks," week")
        if weeks > 1
            horizonstr *= "s"
        end
    end
    if remaindays > 0
        if weeks > 0
            horizonstr *= string(", ",remaindays," day")
        else
            horizonstr *= string("(",remaindays," day")
        end
        if remaindays > 1
            horizonstr *= "s)"
        else
            horizonstr *= ")"
        end
    else
        horizonstr *= ")"
    end
    return horizonstr
end


text:
using Test
using FlightMechanics.Models

euler = [1.690, 1.128, 1.867]
quat = [0.655127, 0.212911, 0.719168, 0.0909671]

att = Attitude(euler...)
@test isapprox(get_euler_angles(att), euler)
@test isapprox(get_quaternions(att), quat, atol=1e-5)

att = Attitude(quat...)
@test isapprox(get_euler_angles(att), euler, atol=1e-5)
@test isapprox(get_quaternions(att), quat)


text:
# # Create Your Model
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/model_create.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/model_create.ipynb)

# RHEOS provide an extensive model library (see Models section). However, if you would like to use a model that has not already been implemented in RHEOS, you can quickly define your own model as described below. Note that currently RHEOS can only handle linear viscoelastic models. Therefore, any added model must satisfy this condition.

# The [`RheoModelClass`](@ref) constructor does not need to contain all moduli functions, depending on what variables are fitted and predicted down the line (relaxation modulus G for stress, creep function J for strain, and complex moduli Gp and Gpp for frequency domain data). As an example, we will assume that only the relaxation modulus is known.

using RHEOS
using PyPlot
#-

PowerLawEmpirical =  RheoModelClass(
                        ## Model name
                        name = "power_empirical",
                        ## Model parameters,
                        p = (:A, :α),
                        ## Relaxation modulus
                        G = quote
                                A*t^(-α)
                            end,
                        ## Network
                        info = "Empirical model"
                        )
#-

## generate timeline
dϵ = timeline(t_start = 1e-2, t_end = 1e1)
## calculates strain data by applying a function of time
dϵ = strainfunction(dϵ, hstep())

## Fix model parameters
powermodel = RheoModel(PowerLawEmpirical, A = 1, α = 0.8)

## Evaluate relaxation response
dpower = modelpredict(dϵ, powermodel)

fig, ax = subplots(1, 1, figsize = (3, 3))
ax.loglog(dpower.t, dpower.σ)
ax.set_xlabel("Time")
ax.set_ylabel("Stress")
ax.set_ylim(bottom = 0.1, top = 20)
ax.grid("on")
#!nb fig #hide

# In a similar way, it is possible to define a model with only creep modulus. In this case, only strain can be fitted/predicted.


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon the Sudoku module being defined first.

using Test

function block_permutations()
    p_reference = Array{Int16}(undef,(4,4))
    #p_reference[:] = [1,2,3,4, 3,4,1,2, 2,3,4,1, 4,1,2,3][:]

    p = Sudoku.puzzle_n(2)
    #@test p == p_reference

    Sudoku.row_block_permute!(p,1,2)
    p_reference[:] = [3,4,1,2, 1,2,3,4, 4,1,2,3, 2,3,4,1][:]
    @test p == p_reference

    Sudoku.col_block_permute!(p,1,2)
    p_reference[:] = [4,1,2,3, 2,3,4,1, 3,4,1,2, 1,2,3,4][:]
    @test p == p_reference
end

function intrablock_permutations()
    p_reference = Array{Int16}(undef,(4,4))
    #p_reference[:] = [1,2,3,4, 3,4,1,2, 2,3,4,1, 4,1,2,3][:]

    p = Sudoku.puzzle_n(2)
    #@test p == p_reference

    Sudoku.row_permute!(p,1,2,1)
    p_reference[:] = [2,1,3,4, 4,3,1,2, 3,2,4,1, 1,4,2,3][:]
    @test p == p_reference

    Sudoku.col_permute!(p,2,1,2)
    p_reference[:] = [2,1,3,4, 4,3,1,2, 1,4,2,3, 3,2,4,1][:]
    @test p == p_reference
end

function mirroring()
    p_reference = Array{Int16}(undef,(4,4))
    #p_reference[:] = [1,2,3,4, 3,4,1,2, 2,3,4,1, 4,1,2,3][:]

    p = Sudoku.puzzle_n(2)
    #@test p == p_reference

    Sudoku.mirror_horizontal!(p)
    p_reference[:] = [4,1,2,3, 2,3,4,1, 3,4,1,2, 1,2,3,4][:]
    @test p == p_reference

    Sudoku.mirror_vertical!(p)
    p_reference[:] = [3,2,1,4, 1,4,3,2, 2,1,4,3, 4,3,2,1][:]
    @test p == p_reference
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
using Graphs
using Ripserer
using Test

using Ripserer: distance_matrix, OneSkeleton

@testset "OneSkeleton respects the Graphs interface" begin
    @testset "no threshold or removed simplices" begin
        flt = Rips(
            [
                0 1 2 3 2 1
                1 0 1 2 3 2
                2 1 0 1 2 3
                3 2 1 0 1 2
                2 3 2 1 0 1
                1 2 3 2 1 0
            ]
        )

        g = OneSkeleton(flt)
        @test eltype(g) ≡ Int
        @test nv(g) == 6
        @test ne(g) == binomial(6, 2)
        @test has_vertex(g, 1)
        @test has_vertex(g, 6)
        @test !has_vertex(g, 7)
        @test !has_vertex(g, 0)
        @test edges(g) == [Edge(j, i) for i in 1:5 for j in (i + 1):6]
        @test vertices(g) == 1:6
        @test is_directed(g) == is_directed(typeof(g)) == false
        @test weights(g) == distance_matrix(flt)
        @test has_edge(g, 1, 2)
        @test has_edge(g, 3, 6)
        @test !has_edge(g, 3, 7)
        @test neighbors(g, 1) == [6, 5, 4, 3, 2]
    end

    @testset "threshold and removed vertices" begin
        dists = [
            0 1 2 3 2 1
            1 0 1 2 3 2
            2 1 0 1 2 3
            3 2 1 0 1 2
            2 3 2 1 0 1
            1 2 3 2 1 0
        ]
        flt = Rips(dists)

        g = OneSkeleton(flt, 2, [Simplex{1}((2, 1), 1), Simplex{1}((3, 1), 2)])
        @test eltype(g) ≡ Int
        @test nv(g) == 6
        @test ne(g) == binomial(6, 2) - 2 - 3
        @test has_vertex(g, 1)
        @test has_vertex(g, 6)
        @test !has_vertex(g, 7)
        @test !has_vertex(g, 0)
        @test edges(g) == [
            Edge(j, i) for i in 1:5 for
            j in (i + 1):6 if dists[i, j] ≠ 3 && (j, i) ∉ [(2, 1), (3, 1)]
        ]
        @test vertices(g) == 1:6
        @test is_directed(g) == is_directed(typeof(g)) == false
        @test weights(g) == distance_matrix(flt)
        @test !has_edge(g, 1, 2)
        @test !has_edge(g, 3, 6)
        @test has_edge(g, 3, 5)
        @test !has_edge(g, 3, 7)
        @test neighbors(g, 1) == [6, 5]
    end

    @testset "simplex as threshold" begin
        dists = [
            0 1 1 1
            1 0 1 1
            1 1 0 1
            1 1 1 0
        ]
        flt = Rips(dists)
        g = OneSkeleton(flt, Simplex{1}((3, 2), 1))

        @test eltype(g) ≡ Int
        @test nv(g) == 4
        @test ne(g) == 4
        @test has_vertex(g, 1)
        @test has_vertex(g, 4)
        @test !has_vertex(g, 7)
        @test !has_vertex(g, 0)
        @test edges(g) == [Edge(j, i) for i in 1:3 for j in (i + 1):4 if (j, i) ≥ (3, 2)]
        @test vertices(g) == 1:4
        @test is_directed(g) == is_directed(typeof(g)) == false
        @test weights(g) == distance_matrix(flt)
        @test !has_edge(g, 1, 2)
        @test !has_edge(g, 3, 6)
        @test has_edge(g, 3, 2)
        @test has_edge(g, 4, 2)
        @test neighbors(g, 1) == [4]
        @test neighbors(g, 4) == [3, 2, 1]
    end

    @testset "Cubical" begin
        image = [
            0 0 0 0 0 0 0
            1 1 1 1 1 1 1
            1 1 2 3 2 1 2
            1 1 2 5 2 1 1
            1 1 1 1 1 1 1
        ]
        flt = Cubical(image)
        g = OneSkeleton(flt, 2)

        @test eltype(g) ≡ Int
        @test nv(g) == *(size(image)...)
        @test ne(g) == 4 * 6 * 2 + 4 + 6 - 7
        @test has_vertex(g, 1)
        @test has_vertex(g, 4)
        @test !has_vertex(g, 36)
        @test !has_vertex(g, 0)
        @test vertices(g) == 1:35
        @test is_directed(g) == is_directed(typeof(g)) == false
        @test weights(g) == distance_matrix(flt)
        @test has_edge(g, 1, 2)
        @test !has_edge(g, 3, 6)
        @test neighbors(g, 1) == [6, 2]
        @test neighbors(g, 14) == [15, 13, 9]
    end
end

@testset "Cycle reconstruction" begin
    @testset "small cycle" begin
        flt = Rips(
            [
                0 1 2 3 2 1
                1 0 1 2 3 2
                2 1 0 1 2 3
                3 2 1 0 1 2
                2 3 2 1 0 1
                1 2 3 2 1 0
            ]
        )

        interval = ripserer(flt; reps=true)[2][1]
        cyc1 = reconstruct_cycle(flt, interval)
        @test eltype(cyc1) == Simplex{1,Int,Int}
        @test all(birth.(cyc1) .== 1)
        @test allunique(cyc1)
        @test length(cyc1) == 6

        cyc2 = reconstruct_cycle(flt, interval, 2)
        @test isempty(cyc2)
    end

    @testset "circle points" begin
        pts = [(sin(t), cos(t)) for t in range(0, 2π; length=22)[1:(end - 1)]]
        flt = Rips(pts)
        interval = ripserer(flt; reps=true)[2][1]

        prev_length = nv(flt)
        for r in range(birth(interval), death(interval); length=10)[1:(end - 1)]
            cyc = reconstruct_cycle(flt, interval)
            @test maximum(birth, cyc) ≤ r

            # Sum of lengths approximates length around circle
            @test π ≤ sum(birth, cyc) ≤ 2π

            # Length decreases for larger r
            @test length(cyc) ≤ prev_length
            prev_length = length(cyc)

            # Each vertex appears exactly twice
            vxs = sort!(collect(Iterators.flatten(vertices.(cyc))))
            @test vxs[1:2:end] == vxs[2:2:end]
            @test allunique(vxs[1:2:end])
        end
        @test isempty(reconstruct_cycle(flt, interval, death(interval)))
        @test isempty(reconstruct_cycle(flt, interval, birth(interval) - eps(Float64)))
    end

    @testset "circle points over a grid with a hole" begin
        # The idea here is that at birth time, the cycle is a circle and at time 1, the
        # cycle is a square surrounding the hole.
        pts = unique!(
            vcat(
                [(3sin(t), 3cos(t)) for t in range(0, 2π; length=22)[1:(end - 1)]],
                vec([
                    (i - 5, j - 5) for
                    (i, j) in Iterators.product(0.0:10.0, 0.0:10.0) if !(i ∈ 4:6 && j ∈ 4:6)
                ]),
            ),
        )

        flt = Rips(pts)
        interval = sort!(ripserer(flt; reps=true)[2]; by=persistence)[end]

        cyc1 = reconstruct_cycle(flt, interval)
        cyc2 = reconstruct_cycle(flt, interval, 1)
        cyc3 = reconstruct_cycle(flt, interval, 2)
        cyc4 = reconstruct_cycle(flt, interval, 3)

        @test sum(birth, cyc1) ≈ 6π atol = 0.1
        @test sum(birth, cyc2) == 16
        @test sum(birth, cyc3) ≈ 8 + 4 * √2
        @test sum(birth, cyc4) ≈ 8 * √2
    end

    @testset "Cubical" begin
        image = [
            0 0 0 0 0 0 0
            1 1 1 1 1 1 1
            1 1 2 3 2 1 2
            1 1 2 5 2 1 1
            1 1 1 1 1 1 1
        ]
        cf = Cubical(image)

        interval = ripserer(cf; reps=true)[2][1]
        c(t1, t2) = simplex(cf, Val(1), (CartesianIndex(t1), CartesianIndex(t2)))
        @test sort(reconstruct_cycle(cf, interval, 1)) == sort([
            c((5, 2), (4, 2)),
            c((4, 2), (3, 2)),
            c((3, 2), (2, 2)),
            c((2, 2), (2, 3)),
            c((2, 3), (2, 4)),
            c((2, 4), (2, 5)),
            c((2, 5), (2, 6)),
            c((2, 6), (3, 6)),
            c((3, 6), (4, 6)),
            c((4, 6), (5, 6)),
            c((5, 6), (5, 5)),
            c((5, 5), (5, 4)),
            c((5, 4), (5, 3)),
            c((5, 3), (5, 2)),
        ])
        @test sort(reconstruct_cycle(cf, interval, 2)) == sort([
            c((5, 3), (4, 3)),
            c((4, 3), (3, 3)),
            c((3, 3), (2, 3)),
            c((2, 3), (2, 4)),
            c((2, 4), (2, 5)),
            c((2, 5), (3, 5)),
            c((3, 5), (4, 5)),
            c((4, 5), (5, 5)),
            c((5, 5), (5, 4)),
            c((5, 4), (5, 3)),
        ])
        @test sort(reconstruct_cycle(cf, interval, 3)) == sort([
            c((5, 3), (4, 3)),
            c((4, 3), (3, 3)),
            c((3, 3), (3, 4)),
            c((3, 4), (3, 5)),
            c((3, 5), (4, 5)),
            c((4, 5), (5, 5)),
            c((5, 5), (5, 4)),
            c((5, 4), (5, 3)),
        ])
        @test reconstruct_cycle(cf, interval, 4) == reconstruct_cycle(cf, interval, 3)
        @test isempty(reconstruct_cycle(cf, interval, 0))
        @test isempty(reconstruct_cycle(cf, interval, 5))
    end

    @testset "Overriding distance matrix" begin
        flt = Rips([
            0 2 4 2 2
            2 0 2 4 2
            4 2 0 2 6
            2 4 2 0 6
            2 2 6 6 0
        ])

        interval = ripserer(flt; reps=true)[2][1]

        dists = [
            0 2 4 2 0
            2 0 2 4 0
            4 2 0 2 6
            2 4 2 0 6
            0 0 6 6 0
        ]
        cyc1 = reconstruct_cycle(flt, interval, 2)
        cyc2 = reconstruct_cycle(flt, interval, 2; distances=dists)
        @test length(cyc1) == 4
        @test length(cyc2) == 5
    end

    @testset "Errors" begin
        flt = Rips(
            [
                0 1 2 3 2 1
                1 0 1 2 3 2
                2 1 0 1 2 3
                3 2 1 0 1 2
                2 3 2 1 0 1
                1 2 3 2 1 0
            ]
        )
        _, d1, d2 = ripserer(flt; reps=2, dim_max=2)

        @test_throws ArgumentError reconstruct_cycle(flt, d1[1])
        @test_throws ArgumentError reconstruct_cycle(flt, d2[1])
    end
end


text:
struct NLsolveWrapper
    zero::Vector{Float64}
    converged::Bool
    failed::Bool
end

NLsolveWrapper() = NLsolveWrapper(Vector{Float64}(), false, true)
converged(sol::NLsolveWrapper) = sol.converged
failed(sol::NLsolveWrapper) = sol.failed

function _get_model_closure(model::SystemModel{MassMatrixModel}, ::Vector{Float64})
    return (residual, x) -> model(residual, x, nothing, 0.0)
end

function _get_model_closure(model::SystemModel{ResidualModel}, x0::Vector{Float64})
    dx0 = zeros(length(x0))
    return (residual, x) -> model(residual, dx0, x, nothing, 0.0)
end

function _nlsolve_call(
    initial_guess::Vector{Float64},
    f_eval::Function,
    jacobian::JacobianFunctionWrapper,
    f_tolerance::Float64,
    solver::Symbol,
)
    df = NLsolve.OnceDifferentiable(
        f_eval,
        jacobian,
        initial_guess,
        similar(initial_guess),
        jacobian.Jv,
    )
    sys_solve = NLsolve.nlsolve(
        df,
        initial_guess;
        xtol = NLSOLVE_X_TOLERANCE,
        ftol = f_tolerance,
        method = solver,
        iterations = MAX_NLSOLVE_INTERATIONS,
        show_trace = true,
    ) # Solve using initial guess x0
    return NLsolveWrapper(sys_solve.zero, NLsolve.converged(sys_solve), false)
end

function _nlsolve_call(
    initial_guess::Vector{Float64},
    f_eval::Function,
    f_tolerance::Float64,
    solver::Symbol,
)
    sys_solve = NLsolve.nlsolve(
        f_eval,
        initial_guess;
        xtol = NLSOLVE_X_TOLERANCE,
        ftol = f_tolerance,
        iterations = MAX_NLSOLVE_INTERATIONS,
        method = solver,
    ) # Solve using initial guess x0
    return NLsolveWrapper(sys_solve.zero, NLsolve.converged(sys_solve), false)
end

function _convergence_check(sys_solve::NLsolveWrapper, tol::Float64, solv::Symbol)
    if converged(sys_solve)
        @info(
            "Initialization non-linear solve succeeded with a tolerance of $(tol) using solver $(solv). Saving solution."
        )
    else
        @warn(
            "Initialization non-linear solve convergence failed, initial conditions do not meet conditions for an stable equilibrium.\nAttempting again with reduced numeric tolerance and using another solver"
        )
    end
    return converged(sys_solve)
end

function _sorted_residuals(residual::Vector{Float64})
    if isapprox(sum(abs.(residual)), 0.0, atol = STRICT_NLSOLVE_F_TOLERANCE)
        @debug "Residual is zero with tolerance $(STRICT_NLSOLVE_F_TOLERANCE)"
        return
    end
    ix_sorted = sortperm(abs.(residual), rev = true)
    show_residual = min(10, length(residual))
    for i in 1:show_residual
        ix = ix_sorted[i]
        @debug ix abs(residual[ix])
    end
    return
end

function _check_residual(
    residual::Vector{Float64},
    inputs::SimulationInputs,
    tolerance::Float64,
)
    @debug _sorted_residuals(residual)
    val, ix = findmax(residual)
    @info "Residual from initial guess: max = $(val) at $ix, total = $(sum(residual))"
    if sum(residual) > tolerance
        state_map = make_global_state_map(inputs)
        gen_name = ""
        state = ""
        for (gen, states) in state_map
            for (state_name, index) in states
                if index == ix
                    gen_name = gen
                    state = state_name
                end
            end
        end
        @warn "The initial residual in $ix of the NLsolve function has a value of $val.
               Generator = $gen_name, state = $state.
               Expect convergence problems"
    end
    return
end
function refine_initial_condition!(
    sim::Simulation,
    model::SystemModel,
    jacobian::JacobianFunctionWrapper,
)
    @assert sim.status != BUILD_INCOMPLETE

    if sim.status == SIMULATION_INITIALIZED
        @info "Simulation already initialized. Refinement not executed"
        return
    end
    converged = false
    initial_guess = get_initial_conditions(sim)
    inputs = get_simulation_inputs(sim)
    bus_range = get_bus_range(inputs)
    powerflow_solution = deepcopy(initial_guess[bus_range])
    f! = _get_model_closure(model, initial_guess)
    residual = similar(initial_guess)
    f!(residual, initial_guess)
    _check_residual(residual, inputs, MAX_INIT_RESIDUAL)
    for tol in [STRICT_NLSOLVE_F_TOLERANCE, RELAXED_NLSOLVE_F_TOLERANCE]
        if converged
            break
        end
        for solv in [:trust_region, :newton]
            @debug "Start NLSolve System Run with $(solv) and F_tol = $tol"
            sys_solve = _nlsolve_call(initial_guess, f!, jacobian, tol, solv)
            #sys_solve = _nlsolve_call(initial_guess, f!, tol, solv)
            failed(sys_solve) && return BUILD_FAILED
            converged = _convergence_check(sys_solve, tol, solv)
            @debug "Write initial guess vector using $solv with tol = $tol convergence = $converged"
            initial_guess .= sys_solve.zero
            if converged
                break
            end
        end
    end

    sim.status = check_valid_values(initial_guess, inputs)
    if sim.status == BUILD_FAILED
        error(
            "Initial conditions refinement failed to find a valid initial condition. Run show_states_initial_value on your simulation",
        )
    end

    f!(residual, initial_guess)
    if !converged || (sum(residual) > MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE)
        _check_residual(residual, inputs, MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE)
        @warn("Initialization didn't found a solution to desired tolerances.\\
              Initial conditions do not meet conditions for an stable equilibrium. \\
              Simulation might fail")
    end

    pf_diff = abs.(powerflow_solution .- initial_guess[bus_range])
    if maximum(pf_diff) > MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE
        @warn "The resulting voltages in the initial conditions differ from the power flow results"
    end

    return
end


text:
# class of constant loss functions for scaling
# (e.g. to compare with other sem software)

struct SemConstant{C} <: LossFunction
    constant::C
end

function (constant::SemConstant)(par, model)
    return constant.constant
end


text:
using StableRNGs
using StableRNGs: LehmerRNG
using Test
using Random

include("streams.jl")

@testset "initialization" begin
    @test StableRNG === LehmerRNG
    @test_throws UndefKeywordError LehmerRNG()
    rng = LehmerRNG(0)
    @test_throws MethodError Random.seed!(rng)
    @test_throws ArgumentError LehmerRNG(rand(typemin(Int):-1))
    @test_throws ArgumentError Random.seed!(rng, rand(typemin(Int):-1))
    @test_throws ArgumentError Random.seed!(rng, big(typemax(UInt64))+1)

    for seed in UInt64[0, 1, 2, 3, 4, typemax(Int32),
                      Int64(2)^32, typemax(Int64), typemax(UInt64)]
        for T = (Int8, UInt8, Int32, Int64, UInt, Int128, BigInt)
            Base.hastypemax(T) && seed > typemax(T) && continue
            seed = T(seed)
            rng = LehmerRNG(seed)
            @test rng isa LehmerRNG
            @test isodd(rng.state)
            state = rng.state
            rng2 = Random.seed!(rng, seed)
            @test rng2 === rng
            @test isodd(rng.state)
            @test rng.state == state
        end
    end

    # state
    rng = LehmerRNG(0)
    @test rng.state == 1
    @test string(rng) == "LehmerRNG(state=0x00000000000000000000000000000001)"
    @test LehmerRNG(state=0x00000000000000000000000000000001).state == 1
    # @test LehmerRNG(state=0x00000000000000000000000000000001) == rng
    @test_throws ArgumentError LehmerRNG(state=2*UInt128(rand(1:9999)))

    # copy / == / hash
    seed = rand(UInt)
    Random.seed!(rng, seed)
    c = copy(rng)
    @test c.state == rng.state
    @test c !== rng
    @test c == rng
    @test hash(c) == hash(rng)
    @test hash(c, UInt(123)) == hash(rng, UInt(123))
    rand(rng)
    @test c != rng
    @test hash(c) != hash(rng) # should almost always be true
    c = copy(rng)
    @test c.state == rng.state
    @test c !== rng
    @test c == rng
end

getsampler(::Type{T}) where {T} = T
getsampler(::Type{Normal{T}}) where {T} = T
getsampler(::Type{Exponential{T}}) where {T} = T
getsampler(x::UnitRange) = x

@testset "$T streams" for T =
    [Bool, Base.BitInteger_types...,
     [map(t, 1:123) for t in Base.BitInteger_types]...,
     1:Int128(2)^100, 1:UInt128(2)^100, false:true,
     Float64, Normal{Float64}, Exponential{Float64},
     Float32, Normal{Float32}, Exponential{Float32},
     Float16, Normal{Float16}, Exponential{Float16}]

    streams = T isa DataType && T <: Integer ? STREAMS[UInt64] :
        T isa UnitRange ? STREAMS[last(T)] : STREAMS[T]

    _rand, _rand! =
        if T isa DataType && T <: Normal
            randn, randn!
        elseif T isa DataType && T <: Exponential
            randexp, randexp!
        else
           rand, rand!
        end
    TT = getsampler(T)

    for (seed, stream) in streams
        if T isa DataType && T <: Integer
            if sizeof(T) == 16
                stream = reinterpret(T, stream)
            else
                stream = stream .% T
            end
        end
        rng = StableRNG(seed)
        n = length(stream)
        a = _rand(rng, TT, n)
        @test a == stream
        Random.seed!(rng, seed)
        @test _rand(rng, TT, n) == stream
        Random.seed!(rng, seed)
        if TT isa DataType
            @test _rand!(rng, a) == stream
        else
            @test _rand!(rng, a, TT) == stream
        end
        Random.seed!(rng, seed)
        @test [_rand(rng, TT) for _=1:n] == stream
    end
end

@testset "rand in range" begin
    rng = LehmerRNG(0)
    @test_throws ArgumentError rand(rng, 1:0)
    for T in Base.BitInteger_types
        sp = Random.Sampler(LehmerRNG, T(1):T(1))
        @test sp isa StableRNGs.SamplerRangeFast
    end
    for n in Int128[rand(1:128, 100); rand(1:Int64(2)^32, 100); Int64(2)^32;
                    rand(1:typemax(Int64), 100); rand(1:Int128(2)^64, 100);
                    Int128(2)^64; rand(1:typemax(Int128), 100)]
        for T in Base.BitInteger_types
            n < typemax(T) || continue
            n = T(n)
            @test all(∈(1:n), rand(rng, T(1):n, 100))
            if n >= 4
                @test all(∈(4:n), rand(rng, T(4):n, 100))
            end
            if T <: Signed
                @test all(∈(-100:n), rand(rng, T(-100):n, 100))
            end
        end
    end
end

@testset "`shuffle` stability" begin
    a = 1:10
    a_shuffled = [4, 5, 6, 1, 3, 7, 10, 2, 8, 9]

    # If these tests fail due to changes in the algorithm used in `Random.shuffle`,
    # we should vendor the old implementation here to keep it stable.
    # See <https://github.com/JuliaRandom/StableRNGs.jl/issues/10>.
    @test shuffle(StableRNG(10), a) == a_shuffled

    b = collect(a)
    shuffle!(StableRNG(10), b)
    @test b == a_shuffled
end

max_stars_repo_name: rfourquet/StableRNGs

max_issues_repo_name: rfourquet/StableRNGs

max_forks_repo_name: rfourquet/StableRNGs


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

name = "ICU"
version = v"65.1"

# Collection of sources required to build ICU
sources = [
    "https://github.com/unicode-org/icu/releases/download/release-$(version.major)-$(version.minor)/icu4c-$(version.major)_$(version.minor)-src.tgz" =>
    "53e37466b3d6d6d01ead029e3567d873a43a5d1c668ed2278e253b683136d948",
]

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

# Do the native build
(
    cp -r source/ native_build/
    cd native_build
    CC="${CC_BUILD}"
    CXX="${CXX_BUILD}"
    AR="${AR_BUILD}"
    LD="${LD_BUILD}"
    RANLIB="${RANLIB_BUILD}"

    # See https://git.alpinelinux.org/aports/tree/main/icu/APKBUILD?id=334ebffde9dec34becdd628ad56007699e98ea81
    update_configure_scripts
    sed -i -e 's,DU_HAVE_STRTOD_L=1,DU_HAVE_STRTOD_L=0,' configure.ac
    sed -i -e 's,DU_HAVE_STRTOD_L=1,DU_HAVE_STRTOD_L=0,' configure
    for x in ARFLAGS CFLAGS CPPFLAGS CXXFLAGS FFLAGS LDFLAGS; do
        sed -i -e "/^${x} =.*/s:@${x}@::" "config/Makefile.inc.in"
    done

    ./configure --prefix=$prefix --build=${MACHTYPE} \
        ac_cv_prog_ac_ct_AR=${AR} \
        ac_cv_prog_ac_ct_RANLIB=${RANLIB}
    make -j${nproc}
)

# Don't use llvm-ar since ICU doesn't know how to deal with it
if [[ ${target} == *apple* ]] || [[ ${target} == *freebsd* ]]; then
    export AR=/opt/${target}/bin/${target}-ar
fi

# Do the cross build
cd source/
update_configure_scripts
./configure --prefix=$prefix --host=$target \
    --with-cross-build="/workspace/srcdir/icu/native_build"
make -j${nproc}
make install
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libicudata", :libicudata),
    LibraryProduct("libicui18n", :libicui18n),
    LibraryProduct("libicuio", :libicuio),
    LibraryProduct("libicutest", :libicutest),
    LibraryProduct("libicutu", :libicutu),
    LibraryProduct("libicuuc", :libicuuc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

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


text:
using MPFR_wrap

const atol = BigFloat("1e-100000")
const prec = 333_000
const π_MPFR = BigFloat(π, prec)

function brent_salamin()
    old_prec = precision(BigFloat)
    setprecision(BigFloat, prec)
    a_prev, b_prev = big(1.0), 1 / sqrt(big(2.0))
    s_prev, π_prev = big(1.0) / 2, big(4.0)
    π_next = BigFloat()
    k, p2 = 1, big(2)
    while true
        a_next = (a_prev + b_prev) / 2
        b_next = sqrt(a_prev * b_prev)
        c_next = (a_next - b_next) * (a_next + b_next)
        s_next = s_prev - p2 * c_next
        π_next = 2a_next^2 / s_next
        #println("$k -> $pn")
        δπ = π_next - π_prev
        abs(δπ) ≤ atol && break
        a_prev, b_prev, s_prev, π_prev = a_next, b_next, s_next, π_next
        k += 1
        p2 *= 2
    end
    setprecision(BigFloat, old_prec)
    @assert isapprox(π_next, π_MPFR, atol=atol)
    (k, π_next)
end

function inplace_brent_salamin()
    a_prev = BigFloat(1, precision=prec)
    b_prev = BigFloat(2, precision=prec)
    rec_sqrt!(b_prev, b_prev)
    s_prev = BigFloat(0.5, precision=prec)
    π_prev = BigFloat(4, precision=prec)
    a_next = BigFloat(precision=prec)
    b_next = BigFloat(precision=prec)
    c_next = BigFloat(precision=prec)
    s_next = BigFloat(precision=prec)
    π_next = BigFloat(precision=prec)
    δπ = BigFloat(precision=prec)
    k = Int32(1)
    while true
        add!(a_next, a_prev, b_prev)
        div_2!(a_next, a_next, Int32(1))

        mul!(b_next, a_prev, b_prev)
        sqrt!(b_next, b_next)

        sub!(c_next, a_prev, a_next)
        sqr!(c_next, c_next)

        mul_2!(c_next, c_next, k)
        sub!(s_next, s_prev, c_next)

        sqr!(π_next, a_next)
        mul_2!(π_next, π_next, Int32(1))
        div!(π_next, π_next, s_next)

        sub!(δπ, π_next, π_prev)
        abs(δπ) ≤ atol && break

        a_prev, a_next = a_next, a_prev
        b_prev, b_next = b_next, b_prev
        s_prev, s_next = s_next, s_prev
        π_prev, π_next = π_next, π_prev

        k += Int32(1)
    end
    @assert isapprox(π_next, π_MPFR, atol=atol)
    (k, π_next)
end

function inplace_borwein_cubic()
    a_prev = BigFloat(1, precision=prec)
    div!(a_prev, a_prev, Int32(3))
    s_prev = BigFloat(3, precision=prec)
    sqrt!(s_prev, s_prev)
    sub!(s_prev, s_prev, Int32(1))
    div_2!(s_prev, s_prev, Int32(1))
    r_next = BigFloat(precision=prec)
    s_next = BigFloat(precision=prec)
    rsq_next = BigFloat(precision=prec)
    a_next = BigFloat(precision=prec)
    tmp = BigFloat(precision=prec)
    π_prev = BigFloat(0, precision=prec)
    π_next = BigFloat(precision=prec)
    δπ = BigFloat(precision=prec)
    k = Int32(0)
    while true
        pow!(r_next, s_prev, Int32(3))
        sub!(r_next, Int32(1), r_next)
        cbrt!(r_next, r_next)
        mul_2!(r_next, r_next, Int32(1))
        add!(r_next, r_next, Int32(1))
        div!(r_next, Int32(3), r_next)

        sub!(s_next, r_next, Int32(1))
        div_2!(s_next, s_next, Int32(1))

        sqr!(rsq_next, r_next)
        sub!(tmp, rsq_next, Int32(1))
        mul!(tmp, tmp, Int32(3)^k)
        mul!(a_next, rsq_next, a_prev)
        sub!(a_next, a_next, tmp)

        div!(π_next, Int32(1), a_next)
        sub!(δπ, π_next, π_prev)
        abs(δπ) ≤ atol && break

        a_prev, a_next = a_next, a_prev
        s_prev, s_next = s_next, s_prev
        π_prev, π_next = π_next, π_prev

        k += Int32(1)
    end
    @assert isapprox(π_next, π_MPFR, atol=atol)
    (k, π_next)
end


text:
"""
PermTrack(args...)  -> PyPlot Plot

Returns a plot of permeability in a single track

 The next table show the list of variables allowed:

 The next table show the list of variables allowed:
|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Depth|Mandatory| -- | Array{Number,1}| 1D Array of Depths|
|K|Mandatory| -- | Array{Number,1}|1D Array of Permeability Values|
|DepthFrom|Optional| minimum(Depth) | Number |Set the Depth upper limit to plot|
|DepthTo|Optional| maximum(Depth) | Number |  Set the Depth lower limit to plot|
|Dtick|Optional|true|Dtick=Bool|Set the Depth ticks On/Off|
"""
function PermTrack(Depth,Perm;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),Dtick=false)
    permax=gca()
    permax.plot(Perm,Depth,linewidth=1,color="black")
    #permax.set_xlim([0.1,10])
    permax.set_ylim([DepthTo,DepthFrom])
    setp(permax.get_yticklabels(),visible=Dtick)
    permax.set_yticks(range(DepthFrom, stop=DepthTo, length=51),minor=true)
    permax.set_yticks(range(DepthFrom, stop=DepthTo, length=11))
    permax.tick_params("both",labelsize=6)
    permax.set_xlabel("Permeability[md]")
    permax.xaxis.tick_top()
    permax.xaxis.set_label_position("top")
    permax.set_xscale("log")
end


text:
"""
$(SIGNATURES)

Find `\\newcommand` and `\\newenvironment` elements and try to parse what
follows to form a proper Latex command/environment.
Return a list of such elements.

The format is:

    \\newcommand{\\NAME}[NARGS]{DEFINITION}
    \\newenvironment{NAMING}[NARGS]{PRE}{POST}

where [NARGS] is optional.
"""
function find_lxdefs(tokens::Vector{Token}, blocks::Vector{OCBlock})
    # container for definitions
    lxdefs  = Vector{LxDef}()
    # find braces `{` and `}`
    braces  = filter(β -> β.name == :LXB, blocks)
    nbraces = length(braces)
    # keep track of active tokens
    active_tokens = ones(Bool, length(tokens))
    active_blocks = ones(Bool, length(blocks))

    # go over active tokens, stop over the ones that indicate a newcommand
    # or a newenvironment, deactivate the tokens that are within the scope
    # of the definition of the command/environment
    for (i, τ) ∈ enumerate(tokens)
        # skip inactive tokens
        active_tokens[i] || continue
        # look for tokens that indicate a new[com\env]
        τ.name in (:LX_NEWCOMMAND, :LX_NEWENVIRONMENT) || continue
        case = ifelse(τ.name == :LX_NEWCOMMAND, :com, :env)

        # find first brace blocks after the newX (naming)
        fromτ = from(τ)
        k = findfirst(β -> (fromτ < from(β)), braces)

        if case == :com
            # there must be two brace blocks after the newcommand (name, def)
            if isnothing(k) || !(1 <= k < nbraces)
                throw(LxDefError(
                    "Ill formed newcommand (needs two {...})"))
            end
        else
            # there must be three brace blocks after the newenvironment
            if isnothing(k) || !(1 <= k < nbraces - 1)
                throw(LxDefError(
                    "Ill formed newenvironment (needs three {...})"))
            end
        end

        # try to find a number of arg between the two first {...} to see
        # if it may contain something which we'll try to interpret as [.d.]
        rge    = (to(braces[k])+1):(from(braces[k+1])-1)
        lxnarg = 0
        # it found something between the naming brace and the 1st def brace
        # check if it looks like [.d.] where d is a number and . are
        # optional spaces (specification of the number of arguments)
        # see the regex pattern LX_NARG_PAT
        if !isempty(rge)
            inter  = subs(str(braces[k]), rge)
            lxnarg = match(LX_NARG_PAT, inter)
            # see test/utils/errors for examples
            if isnothing(lxnarg.captures[1])
                sinner = strip(inter)
                if !isempty(sinner)
                    # corner case: {...}[ ]{...} (empty)
                    if isnothing(match(r"\[\s*\]", sinner))
                        throw(LxDefError("""
                            Ill formed new(command|environment), the indication
                            for the number of arguments could not be read.
                            Expected '[x]' where 'x' is a number."""))
                    end
                end
            end
            matched = lxnarg.captures[2]
            lxnarg  = isnothing(matched) ? 0 : parse(Int, matched)
        end

        # assign naming / def
        naming_braces = braces[k]
        matched = match(LX_NAME_PAT, content(naming_braces))
        # >>> case 1: newcommand
        if case == :com
            if isnothing(matched) || !startswith(matched.captures[1], '\\')
                throw(LxDefError("""
                    Invalid naming in a newcommand, expected a command name of
                    the form `\\command`."""))
            end
            lxname          = string(matched.captures[1])
            defining_braces = braces[k+1] # valid bc of previous check
            # recover the definition + span
            lx_def = stent(defining_braces)
            to_def = to(defining_braces)
            # store the new command
            push!(lxdefs, LxDef(lxname, lxnarg, lx_def, fromτ, to_def))
        # >>> case 2: newenvironment
        else
            if isnothing(matched) || startswith(matched.captures[1], '\\')
                throw(LxDefError("""
                    Invalid naming in a newenvironment, expected an environment
                    name form `command`."""))
            end
            lxname      = string(matched.captures[1])
            pre_braces  = braces[k+1]
            post_braces = braces[k+2] # both valid due to earlier check
            # extract defs
            pre_def  = stent(pre_braces)
            post_def = stent(post_braces)
            to_def   = to(post_braces)
            # store the new env
            push!(lxdefs, LxDef(lxname, lxnarg, pre_def => post_def,
                                fromτ, to_def))
        end
        # mark any active token in the span as inactive
        first_after = findfirst(t -> (from(t) > to_def), tokens[i+1:end])
        if isnothing(first_after)
            active_tokens[i:end] .= false
        else
            active_tokens[i:(i+first_after-1)] .= false
        end
    end # of enumeration of tokens

    # filter out the stuff that's now marked as inactive by virtue of being
    # part of a newX definition (these things will be inspected later)
    tokens = tokens[active_tokens]
    blocks = blocks[active_blocks]

    # separate the braces from the rest of the blocks, they will be used
    # to define the lxcoms
    braces_mask = map(β -> β.name == :LXB, blocks)
    braces = blocks[braces_mask]
    blocks = blocks[@. ~braces_mask]

    return lxdefs, tokens, braces, blocks
end

"""
$(SIGNATURES)

Get the reference pointing to a `LxDef` corresponding to a given `lxname`.
If no reference is found but `inmath=true`, we propagate and let KaTeX deal
with it. If something is found, the reference is returned and will be accessed
further down.
The boolean that is also returned helps keeps track of whether the command or
environment was defined in the utils.
"""
function get_lxdef_ref(lxname::SubString, lxdefs::Vector{LxDef},
                       inmath::Bool=false, offset::Int=0; isenv=false
                       )::Tuple{Ref,Bool}
    # find lxdefs with matching name
    ks = findall(δ -> (δ.name == lxname), lxdefs)

    # check that the def is before the usage
    fromlx = from(lxname) + offset
    filter!(k -> (fromlx > from(lxdefs[k])), ks)

    # if no definition is found, there are three possibilities:
    #  1. there is a Utils definition
    #  2. we're in math, let the math engine deal with it
    #  3. throw an error
    if isempty(ks)
        # check if defined in utils
        if isdefined(Main, :Utils)
            # env def 'env_***'
            flag = isenv && isdefined(Main.Utils, Symbol("env_$(lxname)"))
            # com def 'lx_***'
            flag |= isdefined(Main.Utils, Symbol("lx_$(lxname[2:end])"))
            flag && return (Ref(nothing), true)
        end
        # if we're here, and in math mode, let the math engine deal with it
        inmath && return (Ref(nothing), false)
        # otherwise throw an error
        throw(LxObjError("""
            Command or environment '$lxname' was used before it was
            defined."""))
    end
    return (Ref(lxdefs, ks[end]), false)
end

"""
$(SIGNATURES)

For a command or an environment, find the arguments (i.e. a sequence of braces
immediately after the opening token).
"""
function find_opts_braces(τ::Token, narg::Int, braces::Vector{OCBlock}, name="")
    # spot where an opening brace is expected, note that we can use exact
    # char algebra here because we know that the last character is '}' w length 1.
    nxtidx = to(τ) + 1
    # try to find one at that place
    b1_idx = findfirst(β -> (from(β) == nxtidx), braces)
    # --> it needs to exist + there should be enough braces left for the options
    if isnothing(b1_idx) || (b1_idx + narg - 1 > length(braces))
        throw(LxObjError("""
            Command/Environment '$name' expects $narg argument(s) and there
            should be no space(s) between the command name and the first brace:
            \\com{arg1}... or \\begin{env}{arg1}...
            """))
    end
    # --> examine candidate braces, there should be no spaces between
    # braces to avoid ambiguities
    cand_braces = braces[b1_idx:b1_idx+narg-1]
    for bidx ∈ 1:narg-1
        if (to(cand_braces[bidx]) + 1 != from(cand_braces[bidx+1]))
            throw(LxObjError("""
                Argument braces should not be separated by space(s):
                \\com{arg1}{arg2}... Verify a '$name' command/environment.
                """))
        end
    end
    # If we get here then we have candidate braces that match the number
    # required which we can just return
    return cand_braces
end

function find_opts_braces_math(τ::Token, narg::Int, braces::Vector{OCBlock}, name="")
    cand_braces_idx = Int[]
    # spot where an opening brace is expected, note that we can use exact
    # char algebra here because we know that the last character is '}' w length 1.
    nxtidx = to(τ) + 1
    for i ∈ 1:narg
        b_i_idx = findfirst(β -> (from(β) == nxtidx), braces)
        # --> it needs to exist
        if isnothing(b_i_idx)
            throw(LxObjError("""
                Command/Environment '$name' expects $narg argument(s) and there
                should be no space(s) between the command name and the first brace:
                \\com{arg1}... or \\begin{env}{arg1}...
                """))
        end
        push!(cand_braces_idx, b_i_idx)
        nxtidx = to(braces[cand_braces_idx[end]]) + 1
    end
    # If we get here then we have candidate braces that match the number
    # required which we can just return
    return braces[cand_braces_idx]
end

"""
$(SIGNATURES)

Find `\\command{arg1}{arg2}...` outside of `xblocks` and `lxdefs`.
"""
function find_lxcoms(tokens::Vector{Token}, lxdefs::Vector{LxDef},
                     braces::Vector{OCBlock}, offset::Int=0;
                     inmath::Bool=false)::Tuple{Vector{LxCom}, Vector{Token}}
    # containers for the lxcoms
    lxcoms   = Vector{LxCom}()
    active_τ = ones(Bool, length(tokens))
    nbraces  = length(braces)

    # go over tokens, stop over the ones that indicate a command
    for (i, τ) ∈ enumerate(tokens)
        active_τ[i] || continue
        τ.name == :LX_COMMAND || continue

        # 1. look for the definition given the command name
        lxname = τ.ss
        lxdefref, utils = get_lxdef_ref(lxname, lxdefs, inmath, offset)

        (:find_lxcoms, "lxname: '$lxname'") |> logger
        (:find_lxcoms, "utils: '$utils'") |> logger

        if utils
            # custom command defined in utils, take a single bracket
            lxnarg = 1
        elseif isnothing(lxdefref[])
            # unrecognised command in math mode, defer to backend
            continue
        else
            lxnarg = getindex(lxdefref).narg
        end

        (:find_lxcoms, "lxnarg: '$lxnarg'") |> logger
        (:find_lxcoms, "braces: '$(content.(braces))'") |> logger

        # 2. explore with # arguments
        # >> no arguments
        if lxnarg == 0
            push!(lxcoms, LxCom(lxname, lxdefref))
            active_τ[i] = false
        # >> there is at least one argument --> find all of them
        else
            if inmath
                arg_braces = find_opts_braces_math(τ, lxnarg, braces, lxname)
            else
                arg_braces = find_opts_braces(τ, lxnarg, braces, lxname)
            end
            # all good, can push it
            from_c = from(τ)
            to_c   = to(arg_braces[end])
            str_c  = subs(str(τ), from_c, to_c)

            if utils
                push!(lxcoms, LxCom(str_c, nothing, arg_braces))
            else
                push!(lxcoms, LxCom(str_c, lxdefref, arg_braces))
            end

            # deactivate tokens in the span of the command (will be reprocessed
            # later on)
            first_after = findfirst(τ -> (from(τ) > to_c), tokens[i+1:end])
            if isnothing(first_after)
                active_τ[i+1:end] .= false
            elseif first_after > 1
                active_τ[i+1:(i+first_after-1)] .= false
            end
        end
    end
    return lxcoms, tokens[active_τ]
end

"""
$(SIGNATURES)

Find `\\begin{xxx}` and `\\end{xxx}` blocks.
"""
function find_lxenv_delims!(tokens::Vector{Token}, blocks::Vector{OCBlock})
    rmidx = Int[]
    for (i, τ) ∈ enumerate(tokens)
        τ.name ∈ (:CAND_LX_BEGIN, :CAND_LX_END) || continue
        # try to get an adjoining brace
        nxtidx = to(τ) + 1
        braceidx = findfirst(β -> (β.name == :LXB && from(β) == nxtidx), blocks)
        # it needs to exist
        if isnothing(braceidx)
            throw(LxObjError("""
                Found a delimiter '\\begin' or '\\end' and expected a name in braces
                after that but didn't find it. There should be no space between the
                delimiter and the brace: \\begin{name}.
                """))
        end
        push!(rmidx, braceidx)
        # get the brace
        brace = blocks[braceidx]
        # form a valid token and replace the candidate by it
        from_c = from(τ)
        to_c = to(brace)
        name = ifelse(τ.name == :CAND_LX_BEGIN, :LX_BEGIN, :LX_END)
        tokens[i] = Token(name, subs(τ.ss.string, from_c, to_c))
    end
    deleteat!(blocks, rmidx)
    return nothing
end

"""
$SIGNATURES

Find active environment blocks between an opening and matching closing
delimiter. These can be nested. See also `find_ocblocks` and `find_lxcoms`,
this is essentially a mix of the two.
"""
function find_lxenvs(tokens::Vector{Token}, lxdefs::Vector{LxDef},
                     braces::Vector{OCBlock}, offset::Int=0;
                     inmath=false)::Tuple{Vector{LxEnv}, Vector{Token}}
    # containers for the lxenvs
    lxenvs   = Vector{LxEnv}()
    active_τ = ones(Bool, length(tokens))
    nbraces  = length(braces)
    ntokens  = length(tokens)

    for (i, τ) ∈ enumerate(tokens)
        # only consider active and opening tokens
        active_τ[i] || continue
        τ.name == :LX_BEGIN || continue

        # 1. extract the environment name and find the definition
        env_name        = envname(τ)
        lxdefref, utils = get_lxdef_ref(env_name, lxdefs, inmath, offset; isenv=true)
        if utils
            # custom env defined in utils, take a single bracket
            lxnarg = 1
        elseif isnothing(lxdefref[])
            # unrecognised environment in math mode, defer to backend
            continue
        else
            lxnarg = getindex(lxdefref).narg
        end

        # 2. explore with # argument
        # >> no arguments
        if lxnarg == 0
            arg_braces = Vector{OCBlock}()
        # >> at least one argument
        else
            if inmath
                arg_braces = find_opts_braces_math(τ, lxnarg, braces, env_name)
            else
                arg_braces = find_opts_braces(τ, lxnarg, braces, env_name)
            end
        end

        # 3. find closing delimiter
        # find the closing token \end{$env_name}
        # note that the braces are not tokens anymore as they're already in OCB
        inbalance = 1
        j = i
        while !iszero(inbalance) && (j < ntokens)
            j += 1
            inbalance += envbalance(tokens[j], env_name)
        end
        if inbalance > 0
            throw(OCBlockError(
                "I found at least one opening delimiter '\\begin{$env_name}' " *
                "that is not closed properly.", context(τ)))
        end
        # Construct and store the LxEnv
        push!(lxenvs,
              LxEnv(
                subs(str(τ), from(τ), to(tokens[j])), # string of the whole env
                ifelse(utils, nothing, lxdefref),     # def (none for custom)
                arg_braces, τ => tokens[j]            # args + delims
                )
             )

        # Mark all tokens in the span of the env as inactive (reproc later)
        active_τ[i:j] .= false
    end
    return lxenvs, tokens[active_τ]
end


text:
# Constraint
function <<(p::T, bound::Number) where T<: AbstractPolynomial
  return Constraint(p, -Inf-im*Inf, bound)
end

function <<(bound::Number, p::T) where T<: AbstractPolynomial
  return Constraint(p, bound, Inf+im*Inf)
end

function >>(p::T, bound::Number) where T<:AbstractPolynomial
  return bound << p
end

function >>(bound::Number, p::T) where T<:AbstractPolynomial
  return p << bound
end

function <<(bound::Number, cstr::Constraint)
  if real(bound) > real(cstr.ub) || imag(bound) > imag(cstr.ub)
    warn("<<(bnd, cstr): Creating a constraint with lower bound ", bound, " and upper bound ", cstr.ub)
  end
  return Constraint(cstr.p, bound, cstr.ub)
end
>>(bound::Number, cstr::Constraint) = cstr << bound

function <<(cstr::Constraint, bound::Number)
  if real(cstr.lb) > real(bound) || imag(cstr.lb) > imag(bound)
    warn("<<(cstr, bnd): Creating a constraint with lower bound ", cstr.lb, " and upper bound ", bound)
  end
  return Constraint(cstr.p, cstr.lb, bound)
end
>>(cstr::Constraint, bound::Number) = bound << cstr

function ==(p::T, x::Number) where T<:AbstractPolynomial
  return Constraint(p, x, x)
end

function Base.print(io::IO, cstr::Constraint)
  if cstr.lb == cstr.ub
    print(io, cstr.p, " = ", cstr.ub)
  else
    print(io, cstr.lb, " < ", cstr.p, " < ", cstr.ub)
  end
end

# Problem
Problem() = Problem(Polynomial(), Dict{String, Constraint}(), Dict{String, Type}())

function Base.print(io::IO, pb::Problem)
  print(io, "▶ variables: ")
  for (varName, typ) in sort(collect(pb.variables), by=x->x[1])
    print(io, Variable(varName, typ), " ")
  end
  println(io, "\n▶ objective: ", pb.objective)
  println(io, "▶ constraints: ")
  for (cstrName, cstr) in sort(collect(pb.constraints), by=x->x[1])
    @printf(io, "%10s: ", cstrName)
    println(io, cstr)
  end
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# setup activate -
# ...

# declare external packages -
using DataFrames
using CSV
using Statistics
using StatsBase
using LinearAlgebra
using Optim

# include my codes -
include(joinpath("base","Types.jl"))
include(joinpath("base","Checks.jl"))
include(joinpath("datastore","Files.jl"))
include(joinpath("analysis", "Pretreatment.jl"))
include(joinpath("analysis", "Regression.jl"))
include(joinpath("analysis", "Classify.jl"))

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit