text:
mutable struct TOProblemPar
    ndim::Int          # dimension: 2 or 3
    n::Int             # num of nodes in ground structure
    m::Int             # num of elements in ground structure
    X::Matrix{Float64} # n x ndim matrix, node coords
    T::Matrix{Int}     # m x 2 matrix, element connectivity
    B::Matrix{Float64} # m x (ndim x 2) matrix, compatability matrix, i.e. gamma_i, the direction cosine, i ∈ [m]

    # TODO: multiple load cases, Array{Matrix{Float64}}
    F::Matrix{Float64} # n x ndim matrix, load

    # TODO: multiple support cases, Array{Matrix{Int}}
    Bound::Matrix{Int} # n x ndim matrix, boundary, = 1 fixed, = 0 free

    dof_map::Matrix{Int} # m x (2xndim) matrix, element-dof (e_id, 1..2xndim)
end

mutable struct TOProblem
    m::Int             # num of nodes in ground structure
    c::Vector{Float64} # m x 1 vector, coeff of liear objective function

    τ::Float64 # stability load factor, .>= 1 to ensure global stability
    E::Float64 # Young's  modulus
    σ_c::Float64 # compression stress limit
    σ_t::Float64 # tension stress limit

    # TODO: multi-load cases
    # τ::Vector{Float64}
    l::Vector{Float64} # m x 1 vector, element length
    # TODO: stored as a sparse vector
    f::Vector{Float64} # nfree x 1 vector, external force mapped in free dof indices

    Γ::SparseMatrixCSC{Float64, Int} # m x nfree sp compatability matrix

    Δ::Array{SparseMatrixCSC{Float64, Int}, 1} # m x sp(nfree x nfree) array, Δ_i for geometric stiffness matrix
    eK::Array{SparseMatrixCSC{Float64, Int}, 1} # m x sp(nfree x nfree), K_i for element stiffness matrix in full dof index form

    x_init::Vector{Float64} # init x value
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 78c5a32c-fef3-11eb-0560-6592672bbca5
using PlutoUI

# ╔═╡ 214bc212-39a2-44e1-8e7f-aa0bbe6e69df

# ╔═╡ 9d4a1119-bbbe-4d26-af88-6789420ab030
md"""# Data Collection as Abstraction

The word abstraction comes to us via the Latin verm _abstrahere_ meaining "to pull", "drag" or "draw away from" (_Information Ages_, p. 4)

>Abstraction: The action of considering something in the abstract, independently of its associations or attributes; the process of isolating properties or characteristics common to a number of diverse objects, events, etc., without reference to the peculiar properties of particular examples or instances. Also: the state of being considered in this way; abstractness. ([OED](oed.com/view/Entry/766?redirectedFrom=abstraction#eid))

-----------

>Generalization: The action or process of forming or expressing a general concept or proposition on the basis of inference from particular instances. ([OED](https://www.oed.com/view/Entry/77505?redirectedFrom=generalization#eid))

------------

>In most dictionaries, the term abstraction is defined as a quality, generality, and ideas. Dictionary source [7] defines abstraction as “the act of considering something as a general quality or characteristics, apart from concrete realities, specific objects or actual instances”. Abstraction is the process of initial concept formulation and generalization of ideas by extracting common qualities from specific examples. Other definitions of abstraction related to generality are “a general idea or quality rather than an actual person, object, or event: an abstract idea or quality” by [8], and “the quality of dealing with ideas rather than events” by [9]. __These definitions all summarize abstraction as the act of generalizing something by taking out only important points from the detailed characteristics of a problem.__ (doi:10.1016/j.procir.2016.01.038)

The challenge is that what constitutes __"only important points"__ is context dependent.
"""

# ╔═╡ 67199872-460d-4674-a8a2-24d19860aeda
md"""# Picasso's Bulls

On November 2, 1945 the Spanish artist Pablo Picasso started working in the lithography workshop of the Parisian Fernand Mourlot. For four months he feverishly worked on creating lithographics, an artistic form he had not worked with extensively in the past. The most famous product of these months is a series of lithographs of a bull. In the first four lithographs, the bull is getting more detailed, more "realistic." But starting with the 5th lithograph, each subsequent lithograph is more sparse, more "abstract" until in the last lithograph, Picasso only uses six (as I count them) lines to depict the bull. According to the art historian Irvine Lavine, Picasso's intent was to

>to reduce the bull to its constituent parts, to recover the disjecta membra of his dream-bull—bred of pure lines to an elemental, disembodied, quintessential bullishness. ([PICASSO'S LITHOGRAPH(S) "THE BULL(S)" AND THE HISTORY OF ART IN REVERSE *](https://publications.ias.edu/sites/default/files/Lavin_PicassoBulls_2010.pdf))
"""

# ╔═╡ 2953b221-96c3-47ad-9292-49812c0939ad
html""""""

# ╔═╡ ef1eb844-9796-4d43-9529-f2daf6bd07f5
md"""

### Question:

Based on your examination of Picasso's last lithograph, what would you say consistitues his

#### Type your answers here
$(@bind what_is_bullishness TextField((80,7)))
"""

# ╔═╡ 40e32ce3-8b9d-4ec9-af2b-8c4715b53556
let bullanswer = md"""

According to my interpretation of the lithograph, to Picasso the essence of a bull is

1. Four legs
1. A tail
1. A large body
1. A pair of horns
1. A penis
1. A scrotum

The technology that Picasso chose to use would limit _what_ he could depict about the bull. The essence of a bull might be different if he was sculpting (3D) or painting (color).
"""

HTML("""My Answers$(html(bullanswer))""")

end

# ╔═╡ 4b8ba289-3cbe-4113-9686-f02e5377de58
md""" ## Remember:
### Abstraction is Context (and Technology) Dependent!

### Example

Here is what a bull might abstract to for a rancher.

__Bonus:__ Can anyone recognize a connection between the bull's name and his year of birth (YOB)?
"""

# ╔═╡ 0c1136e7-53ac-46f5-b58c-4c80da0dcc8b
html"""

    ID	Name	Species	Breed	Sex	Weight	YOB

    53287	Ferdinand	Bos taurus	Pajuna	M	583	1936

"""

# ╔═╡ 87104176-fb30-4542-977a-5551de189938
begin
bcomment = md"""### Final comment on our bull abstraction example

Our abstractions to date have essentially described a "bull" as existing alone in some universe. We have excluded everything about the environment which the bull lives in as well as the bull's society.

This is not a problem, as long as we recognize what we might call the __granularity__ of our model and the context for which that granuality was chosen.
"""

modelcomment= md"""

>However, to build a model we need to choose a level of granularity and try to stay within it, or the whole model would evaporate into the nothingness of truth. Any level of granularity is called a _stereotype_; e.g., we stereotype objects on earth by saying they each have a weight. A stereotype is a lie, more politely a conceptual simplification, that is convenient for the way we want to do business. (Spivak, David I. _Category Theory for the Sciences_ (p. 45). MIT Press.)

"""

HTML("""$(html(bcomment))

A Bull's Environment	A Bull's Social Context

$(html(modelcomment))""")

end

# ╔═╡ b54a636b-e3aa-4012-8e94-a7248e4392fe
md"""
## More Utilitarian Abstraction
### Maps as Data Abstraction

Common examples of data abstractions are maps. What makes a map more useful than, say, a satellite image of the same region?

### __Example:__ Washington D.C. Metro Map

![DC Metrorail Map](https://upload.wikimedia.org/wikipedia/commons/thumb/9/9c/Wash-dc-metro-map.png/420px-Wash-dc-metro-map.png)

#### Questions:

1. What are some of the abstractions that you think are in this map? (You might want to use Google Maps or some other map resource to look at the DC region.)
1. Why do you think thees abstractions were made?
#### Type your answers here
$(@bind resp1 TextField((80,7)))
"""

# ╔═╡ dc07f793-b2ad-4180-814c-b143382a47f5

let mapanswer = md"""
## What is not included in the model?

The map/model gives up absolute distance and to a lesser degree direction in order to simplify the representation of the _ordinal_ location of stop and needs and opportiunities for train transfers.

#### Distance

Limiting ourselves to the redline with a starting point of Metro Central, the Metro map indicates that the Tenleytown-AU stop is roughly have way between Metro Central and the Rockville stop. However, from the Google Map you now see, it looks like the Rockville station is at least 4x farther.

#### Direction

Again on the red line, we see that Rockville is located significantly west of Glenmont. However, from the Google Map we see that Rockville is located north north-west of Glenmont.

### How does the abstraction help its purpose?

- Imagine the difficulty of differentiating all the metro stations in the inner-city region while still identicating the location of the suburban stations.
- You could not use this map to:
	- Estimate driving time between ...
	- ???
"""

HTML("""My Answers$(html(mapanswer))""")

end

# ╔═╡ fb8e8f2e-ecc0-46d8-bc79-b3ed773040b2
md"""
## Abstracting a Patient

When we generate and collect data on a patient, we are following the same steps the map maker is taking. We are ignoring some features about the patient and emphasizing (distorting?) others in order to make a map (abstraction) of that patient that helps us make decisions/take actions regarding that patient more efficiently---within our domain!

## Focusing on our use case

The use-case for this short coure is centered on the the management of type II diabetes patients in a GP clinic. In this context, what do you think the essence of a person is in this context?

What are significant features that a person might feel are essential in their definition of their selfness that would likely not be included in the GP clinic abstraction of them

#### Type your answers here
$(@bind gp_abstraction TextField((80,7)))
"""

# ╔═╡ 0f62a4eb-3f72-4cab-b966-5b6847a118e1
let f_abs = md"""

- [What is up with the spherical cow?](https://www.wired.com/2011/02/what-is-up-with-the-spherical-cow/)
  - The physical world is _very_ complex and so physicists are always simplifying it. This simplification can sometimes be taken to extremes, leading to this very common physics joke.
- The durrells Series 1 Episode 3 filing people by height.
  - This British show is currently (September 2021) [streaming on ABC iView](https://iview.abc.net.au/video/ZW2693A003S00). In this episde the daughter, Margo, gets a job at a surgery where she proceeds to file patient records by patient height as this was the first data element on the forms and thus, presumably most important feature.
- [From A to Z — the surprising history of alphabetical order](https://www.abc.net.au/news/2020-06-11/history-of-alphabetical-order-a-to-z/12320808)
  - The variety of ways we can abstract knowledge relationships

"""

HTML("""
 Further Resources on Abstraction

$(html(f_abs))
""")
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

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

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

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

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

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

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

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

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

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

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

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

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

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

# ╔═╡ Cell order:
# ╟─214bc212-39a2-44e1-8e7f-aa0bbe6e69df
# ╟─9d4a1119-bbbe-4d26-af88-6789420ab030
# ╟─67199872-460d-4674-a8a2-24d19860aeda
# ╟─2953b221-96c3-47ad-9292-49812c0939ad
# ╟─ef1eb844-9796-4d43-9529-f2daf6bd07f5
# ╟─40e32ce3-8b9d-4ec9-af2b-8c4715b53556
# ╟─4b8ba289-3cbe-4113-9686-f02e5377de58
# ╟─0c1136e7-53ac-46f5-b58c-4c80da0dcc8b
# ╟─87104176-fb30-4542-977a-5551de189938
# ╟─b54a636b-e3aa-4012-8e94-a7248e4392fe
# ╟─dc07f793-b2ad-4180-814c-b143382a47f5
# ╟─fb8e8f2e-ecc0-46d8-bc79-b3ed773040b2
# ╟─0f62a4eb-3f72-4cab-b966-5b6847a118e1
# ╟─78c5a32c-fef3-11eb-0560-6592672bbca5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
using CoordinateTransformations, Rotations
using NearestNeighbors
using LinearAlgebra

centerofmass(pts::AbstractVector{T}) where { T <: AbstractVector{<:AbstractFloat} } = #
    [mean(p -> p[i], pts) for i in 1:length(pts[1])]

"""
    icp_iteration(pts1::AbstractVector{<:AbstractVector{2, T}}, pts2::AbstractVector{<:AbstractVector{2, T}}; orthonormal=false) where { T<: AbstractFloat }

Perform a single iteration of the iterative closest point algorithm using the orthogonal Procrustes algorithm.

Returns the `pts2` rotated and translated to most closely match `pts1`.
"""
function icp_iteration(pts1::AbstractVector{T}, pts2::AbstractVector{T}) where { T <: AbstractVector{<:AbstractFloat}}
    # Center both pts1 and pts2 on the origin
    com1, com2 = centerofmass(pts1), centerofmass(pts2)
    pts1c, pts2c = Translation(-com1).(pts1), Translation(-com2).(pts2)
    # Find the nearest point in pts1c for each point in pts2c
    pts1cr = map(i->pts1c[i[1]], knn(KDTree(pts1c), pts2c, 1)[1])
    @assert length(pts1cr)==length(pts2c)
    # Turn Vector{Vector{<:AbstractFloat}} into Matrix{<:AbstractFloat}
    m2(vs) = transpose(reshape(collect(Iterators.flatten(vs)), length(vs[1]), length(vs)))
    # Solve the orthogonal Procrustes problem
    f = svd(m2(pts1cr) * transpose(m2(pts2c)))
    Ω = f.U * f.Vt
    # Now translate and rotated the reordered data towards `pts1`
    return Translation(com1).(Ω * pts2c)
end

"""
    icp(pts1::AbstractVector{T}, pts2::AbstractVector{T}; maxiter=10, tol=0.9)

Implements the iterative closest point algorithm.  Transforms `pts2` through rotations and
translations to come as close as possible to `pts1` using a least-squares metric.

The intention is that `pts1` is the super-set of points in `pts2`.  When aligned most of
the points in `pts2` will match up with points in `pts1`.
"""
function icp(pts1::AbstractVector{T}, pts2::AbstractVector{T}; maxiter=10, tol=0.99) where { T <: AbstractVector{<:AbstractFloat}}
    initialerror = icperror(pts1, pts2)
    next, res, minerror, pos = pts2, pts2, initialerror, 0
    @show initialerror
    for _ in 1:maxiter
        next = icp_iteration(pts1, next)
        nexterror = icperror(pts1, next)
        if nexterror < minerror
            # Meets tolerance
            (nexterror > tol * minerror) && return next
            res, minerror, pos = next, nexterror, 0
        else
            # Increasing error
            (minerror!=initialerror) && ((pos+=1)==2) && return res
        end
    end
    minerror == initialerror && @warn "No improvement after $maxiter steps in icp(...). Returning the initial points."
    return res
end
"""
    icperror(pts1::AbstractVector{T}, pts2::AbstractVector{T}) where { T <: AbstractVector{<:AbstractFloat}}

Measures the sum distance between the points in `pts2` and the point closest to the point in `pts1`.
"""
icperror(pts1::AbstractVector{T}, pts2::AbstractVector{T}) where { T <: AbstractVector{<:AbstractFloat}} =
    sum(i->i[1], knn(KDTree(pts1), pts2, 1)[2])

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
import Base.transpose

"""
    transpose(x)
"""
function transpose(x::Var)
    y = transpose(x.data)
    Var(y, (transpose,x))
end

function addgrad!(y::Var, ::typeof(transpose), x::Var)
    isvoid(x.grad) && return
    addto!(transpose(y.grad), x.grad)
end


text:
module Rando

export RandoData, fv

#=
Coefficients for a simple, but reproductible random number generator.
=#
mutable struct RandoData
    seed::Int64
    a::Int64
    c::Int64
    m::Int64
    RandoData() = new(123456789, 1103515245, 12345, 4294967296)
end

#=
return a random float in [0,vmax]
=#
function fv!(R::RandoData, vmax=1.)
    R.seed = (R.a * R.seed + R.c) % R.m
    vmax * Float64(R.seed) / R.m
end

end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
module DnDstats
using DataFrames

export compare

weaponlist = DataFrame(
name = ["Club", "Dagger", "Great-Club", "Hand-Axe", "Javelin", "Light Hammer", "Mace", "Quarterstaff", "Sickle", "Spear", "Crossbow, light", "Dart", "Shortbow", "Sling", "Battleaxe", "Flail", "Glaive", "Greataxe", "Great-Sword", "Halberd", "Lance", "Long-sword", "Maul", "Morning-star", "Pike", "Rapier", "Scimitar", "Short-sword", "Trident", "War pick", "War-hammer", "Whip", "Blowgun", "Crossbow, hand", "Crossbow, heavy", "Longbow", "Net"],
cost = [0.1, 2, 0.2, 5, 0.5, 2, 5, 0.2, 1, 1, 25, 0.05, 25, 0.1, 10, 10, 20, 30, 50, 20, 10, 15, 10, 15, 5, 25, 25, 10, 5, 5, 15, 2, 10, 75, 50, 50, 1],
damagedicenumber = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
damagedicetype = [4, 4, 8, 6, 6, 4, 6, 6, 4, 6, 8, 4, 6, 4, 8, 8, 10, 12, 6, 10, 12, 8, 6, 8, 10, 8, 6, 6, 6, 8, 8, 4, 1, 6, 10, 8, 0],
damagetype = ["Bludgeoning", "Piercing", "Bludgeoning", "Slashing", "Piercing", "Bludgeoning", "Bludgeoning", "Bludgeoning", "Slashing", "Piercing", "Piercing", "Piercing", "Piercing", "Piercing", "Slashing", "Bludgeoning", "Slashing", "Slashing", "Slashing", "Slashing", "Piercing", "Slashing", "Bludgeoning", "Piercing", "Piercing", "Piercing", "Slashing", "Piercing", "Piercing", "Piercing", "Bludgeoning", "Slashing", "Piercing", "Piercing", "Piercing", "Piercing", "None"],
weight = [2, 1, 10, 2, 2, 2, 4, 4, 2, 3, 5, .25, 2, 0, 4, 2, 6, 7, 6, 6, 6, 3, 10, 4, 18, 2, 3, 2, 4, 2, 2, 3, 1, 3, 18, 2, 3],
complexity = ["simple", "simple", "simple", "simple", "simple", "simple", "simple", "simple","simple", "simple", "simple", "simple", "simple", "simple", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial", "martial"],
range = ["melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "ranged", "ranged", "ranged", "ranged", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "melee", "ranged", "ranged", "ranged", "ranged", "ranged"],
shortrange = [5, 20, 5, 20, 30, 20, 5, 5, 5, 20, 80, 20, 80, 30, 5, 5, 10, 5, 5, 10, 10, 5, 5, 5, 10, 5, 5, 5, 20, 5, 5, 10, 25, 30, 100, 150, 5],
longrange = [5, 60, 5, 60, 120, 60, 5, 5, 5, 60, 320, 60, 320, 120, 5, 5, 10, 5, 5, 10, 10, 5, 5, 5, 10, 5, 5, 5, 60, 5, 5, 10, 100, 120, 400, 600, 15],
ammunition = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0],
finesse = [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0],
heavy = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
light = [1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0],
loading = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
reach = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
special = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
thrown = [0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
twohanded = [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
versatile = [0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0]
), false

function chancetohit(attackbonus,ac)
    if attackbonus == "auto"
        return 1
    end
    hits = 1
    for roll in 2:19
        if roll + attackbonus >= ac
            hits = hits + 1
        end
    end
    return hits/20
end

function damage(dice,type,modifier)
    dicetypes = length(dice)
    dicedamage = 0
    for x in 1:dicetypes
        dicedamage = dicedamage + (((type[x]+1)/2)*dice[x])
    end
    return dicedamage + modifier
end

function avgdamage(attackbonus,ac,dice,type,modifier, crit = true)
    percenthit = chancetohit(attackbonus,ac)
    damageperhit = damage(dice,type,modifier)
    if crit == true
        crit = damage(2*dice,type,modifier)
        return damageperhit*(percenthit-0.05)+(crit*0.05)
    else
        return damageperhit*percenthit
    end
end

function damageperattack(weapon1, acrange = 8:20)
    totaldamage1 = 0
    for ac in acrange
        damage1 = avgdamage(weapon1[2], ac, weapon1[3], weapon1[4], weapon1[5], weapon1[6])
    totaldamage1 = totaldamage1 + damage1
    end
    avgdamage1 = totaldamage1/length(acrange)
    return avgdamage1
end

function compare(weapon1, weapon2, weapon1attacks, weapon2attacks, acrange = 8:20) #each weapon should be a length 6 vector [name, attackbonus, dice, type, modifier, crit]
    name1 = weapon1[1]
    name2 = weapon2[1]
    avgdamage1 = damageperattack(weapon1, acrange)*weapon1attacks
    avgdamage2 = damageperattack(weapon2, acrange)*weapon2attacks
    if avgdamage1 > avgdamage2
        println("The $name1 is more powerful than the $name2.")
        println("The $name1 deals $avgdamage1 and the $name2 deals $avgdamage2")
        return weapon1[1]
    elseif avgdamage1 < avgdamage2
        println("The $name2 is more powerful than the $name1.")
        println("The $name2 deals $avgdamage2 and the $name1 deals $avgdamage1")
        return weapon2[1]
    else
        println("The $name1 does as much damage as the $name2.")
        println("The $name1 deals $avgdamage1 and the $name2 deals $avgdamage2")
        return "Neither"
    end
end

end


text:
using VLKeggSDK
using DataFrames

# setup an array of ec numbers -
rn_number_array = [
    "rn:R00299"
    "rn:R00771"
    "rn:R04779"
    "rn:R01068"
    "rn:R11945"
]

# download these reactions -
reaction_array = get_reaction_for_rn_number(rn_number_array) |> check |> DataFrame


text:
using Revise

using BNPRegressionGGA2021
using Distributions
using Gadfly
using Random

const BNP = BNPRegressionGGA2021

#=================================#
# Example 1 - 1 trivial predictor #
#=================================#

function simulate_sample(N0, N1)
    dy = MixtureModel(Beta, [(1.0, 5.0), (10.0, 2.0)], [0.2, 0.8])
    X0 = ones(N0, 1)
    X1 = ones(N1, 1)
    y0 = rand(dy, N0)
    y1 = LinRange(0.0, 1.0, N1) |> collect
    return dy, y0, X0, y1, X1
end;

Random.seed!(1);
N0, N1 = 500, 50;
dy, y0, X0, y1, X1 = simulate_sample(N0, N1);
m = BNP.DGPMBeta(; y0, X0, y1, X1);
chainf, chainβ = BNP.sample!(m; mcmcsize = 10000, burnin = 5000);
plot(
    layer(x = y1, y = mean(chainf), Geom.line, color = ["bnp"]),
    layer(x = y1, y = pdf.(dy, y1), Geom.line, color = ["true"]),
    layer(x = y0, Geom.histogram(density = true, bincount = 50), color = ["hist"]),
)

#===================================================#
# Example 2 - 4 trivial and 1 non-trivial predictor #
#===================================================#

# Simulate a sample
function simulate_sample(N0, N1)
    dy1 = MixtureModel(Beta, [(1.0, 5.0), (10.0, 2.0)], [0.2, 0.8])
    dy2 = MixtureModel(Beta, [(1.0, 5.0), (10.0, 2.0)], [0.8, 0.2])
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    X1 = kron(ones(N1), [ones(2, 5) [1, 0]])
    y0 = [X0[i, end] == 1 ? rand(dy2) : rand(dy1) for i in 1:N0]
    y1 = LinRange(0, 1, N1) |> collect |> x -> repeat(x, size(X1, 1) ÷ N1)
    return dy1, dy2, y0, X0, y1, X1
end

Random.seed!(1);
N0, N1 = 500, 50;
dy1, dy2, y0, X0, y1, X1 = simulate_sample(N0, N1);
ygrid = LinRange(0, 1, N1) |> collect;
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid), Geom.line, color=["pdf1"]),
)
m = BNP.DGPMBeta(; y0, X0, y1, X1);
chainf, chainβ, chaing = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, end]), Geom.line)


text:
include("day8_solution.jl")
using Test


text:
function T_bar_equation(cur_reactor::AbstractReactor, cur_value::Number)
  cur_G_T = cur_value

  cur_G_T /= cur_reactor.T_bar

  cur_equation = Equation(0, 0, 0, cur_G_T)

  cur_equation
end


text:
using PowerPlots
using PowerModels
using VegaLite

## Basic plot
data = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")
p = powerplot(data)
save("basic_plot.html", p)

## Set Color
p = powerplot(data, bus_color="orange", gen_color=:purple, branch_color="#AFAFAF")
save("color_plot.html", p)

p = powerplot(data; node_color="red", edge_color="purple")
save("color_override_plot.html", p)

## Set Size
p = powerplot(data, bus_size=1000, gen_size=100, branch_size=2, connector_size=10)
save("size_plot.html", p)

p = powerplot(data, node_size=1000, edge_size=10)
save("size_override_plot.html", p)

## Set data
p = powerplot(data, bus_data="va", bus_data_type="quantitative", bus_color=["gray","red"])
save("size_plot.html", p)


text:
vtype = Union{Rational{Int}, Int}
struct jp
    j::vtype
    p::Char
end

jp(v::Tuple{vtype, Char}) = jp(v[1],v[2])

import Base: length
length(jp1::jp) = 0
two_j(jp::jp) = Int(2jp.j)
#
# dealing with spin 1/2
x2(v) = @. Int(2v)
#
_over2(two) = mod(two,2)==0 ? div(two,2) : two // 2
over2(two) = _over2.(two)

⊗(p1::Char,p2::Char) = p1==p2 ? '+' : '-'
⊗(jp1::jp, jp2::jp) = [jp(j, ⊗(jp1.p, jp2.p)) for j in abs(jp1.j-jp2.j):abs(jp1.j+jp2.j)]

function str2jp(pin::String)
    p = filter(x->x!='^', pin)
    !(contains(p, '/')) && return jp(Meta.parse(p[1:end-1]), p[end])
    p[end-2:end-1] != "/2" && error("the string should be `x/2±`, while it is $(p)")
    two_j = Meta.parse(p[1:end-3])
    !(typeof(two_j) <: Int) && error("the string should be `x/2±`, while it is $(p)")
    return jp(two_j//2, p[end])
end

macro jp_str(p)
    return str2jp(p)
end

possible_ls((jp, (jp1, jp2))::Pair{jp, Tuple{jp, jp}}) = possible_ls(jp1, jp2; jp)

function possible_ls(jp1::jp, jp2::jp; jp::jp)
    ls = Vector{Tuple{Int,vtype}}(undef,0)
    for s in abs(jp1.j-jp2.j):abs(jp1.j+jp2.j)
        for l in Int(abs(jp.j-s)):Int(abs(jp.j+s))
            if jp1.p ⊗ jp2.p ⊗ jp.p  == (isodd(l) ? '-' : '+')
                push!(ls, (l,s))
            end
        end
    end
    return sort(ls, by=x->x[1])
end

function possible_lsLS(k::Int, jpR::jp, jps::T where T<:Vector{jp})
    i,j = ij_from_k(k)
    lsv = possible_ls(jps[i], jps[j]; jp = jpR)
    LSv = possible_ls(jpR, jps[k]; jp = jps[4])
    return [(ls=ls, LS=LS) for (ls, LS) in Iterators.product(lsv, LSv)]
end
function possible_lsLS(k::Int, two_s::Int, parity::Char, two_js , Ps)
    jpR = jp(two_s//2, parity)
    jps = jp.(zip(two_js .// 2, Ps))
    return possible_lsLS(k, jpR, jps)
end

jls_coupling(two_j1, two_λ1, two_j2, two_λ2, two_j, two_l, two_s) =
    sqrt((two_l+1) / (two_j+1)) *
    CG_doublearg(two_j1, two_λ1, two_j2, -two_λ2, two_s, two_λ1-two_λ2) *
    CG_doublearg(two_l,  0, two_s, two_λ1-two_λ2, two_j, two_λ1-two_λ2)

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
using PDB
using FactCheck

const testdir = dirname(@__FILE__)

facts("148l") do
    open(joinpath(testdir, "pdb", "148l.pdb")) do f
        ent = parse_entry(f)
        @fact name(ent) => "148L"
        @fact length(models(ent)) => 1
        atoms = collect(eachatom(first(models(ent))))
        @fact atoms[1].serial => 1
        @fact atoms[2].serial => 2
    end
end

facts("7ahl") do
    open(joinpath(testdir, "pdb", "7ahl.pdb")) do f
        ent = parse_entry(f)
        @fact name(ent) => "7AHL"
        @fact length(models(ent)) => 1
        atoms = collect(eachatom(first(models(ent))))
        @fact atoms[1].serial => 1
        @fact atoms[2].serial => 2
    end
end


text:
# basis functions
linear(x::Number) = vcat(1, x)
linear(x::AbstractArray{T, 1}) where {T} = hcat(linear.(x)...)'

gaussian(x::Number, μ, σ=0.1) = vcat(1, exp.(-0.5*(x .- μ).^2/σ^2))
gaussian(x::AbstractArray{T, 1}, μ, σ=0.1) where {T} = hcat(gaussian.(x, Ref(μ), Ref(σ))...)'

polynomial(x::Number, r) = x.^r
polynomial(x::AbstractArray{T, 1}, r) where {T} = hcat(polynomial.(x, Ref(r))...)'

"""
    posterior(Y, ϕX, α, β)

Compute the posterior mean and variance given the observations Y,
regressor matrix ϕX, prior precision α, and the additive noise precision β.
"""
function posterior(Y, ϕX, α, β)
    Σ = inv(α*I + β*ϕX'*ϕX)
    μ = β*Σ*ϕX'*Y
    return μ, Σ
end

"""
    predict(ϕx, μ, Σ, β)

Compute the predictive distribution and standard deviation given the
regressed test point ϕx, posterior mean μ, posterior variance Σ, and
the additive noise precision β.
"""
function predict(ϕx, μ, Σ, β)
    μy = μ'*ϕx
    Σy = 1/β + ϕx'*Σ*ϕx
    return μy, √Σy
end

"""
    regress(Xo, Yo, Xt, Yt, ϕ, α, β)

Plot the point predictions given the observed data (Xo, Yo), test data (Xt, Yt),
basis functions ϕ, prior precision α, and the additive noise precision β.
"""
function regress(Xo, Yo, Xt, Yt, ϕ, α, β, ϕargs...)
    μ, Σ = posterior(Yo, ϕ(Xo, ϕargs...), α, β)

    # predictions
    Yμ = similar(Xt); Yσ = similar(Xt)
    for (i, xt) in enumerate(Xt)
        yμ, yσ = predict(ϕ(xt, ϕargs...), μ, Σ, β)
        Yμ[i] = yμ; Yσ[i] = yσ
    end

    Prediction1D(Xo, Yo, Xt, Yt, Yμ, Yσ)
end

struct Prediction1D{T}
    Xo::T
    Yo::T
    Xt::T
    Yt::T
    μ::T
    σ::T
end

@recipe function f(p::Prediction1D)
    legend --> false
    @series begin
        fillcolor --> :grey
        fillalpha --> 0.2
        linewidth --> 1.0
        linecolor --> :grey
        ribbon := p.σ
        (p.Xt, p.μ)
    end
    @series begin
        linewidth --> 2.0
        linecolor --> :black
        linestyle --> :dash
        (p.Xt, p.Yt)
    end
    @series begin
        seriestype --> :scatter
        markercolor --> :red
        markershape --> :circle
        markersize --> 3.0
        markerstrokecolor --> :red
        (p.Xo, p.Yo)
    end
end

function linefit(N, α=2.0, β=25.0; seeded=false)
    f(x, β=Inf) = -0.3 + 0.5*x + 1/β*randn()
    if seeded
        seed!(1729)
    end
    Xo = 2.0*rand(N) .- 1.0
    Yo = f.(Xo, β)
    Xt = collect(-1.0:0.01:1.0)
    Yt = f.(Xt)
    regress(Xo, Yo, Xt, Yt, linear, α, β)
end

function gaussfit(N, α=1.0, β=25.0, μ=0:0.1:3; seeded=false)
    f(x, β=Inf) = 0.5 + x*sin(2*π*x) + 1/β*randn()
    if seeded
        seed!(1729)
    end
    Xo = rand(N)*3.0
    Yo = f.(Xo, β)
    Xt = collect(0.0:0.01:3.0)
    Yt = f.(Xt)
    regress(Xo, Yo, Xt, Yt, gaussian, α, β, μ)
end

function polyfit(N, α=0.005, β=1/(0.3^2), r=0:5; seeded=false)
    f(x, β=Inf) = 0.5 + sin(2*π*x) + 1/β*randn()
    if seeded
        seed!(1729)
    end
    Xo = rand(N)
    Yo = f.(Xo, β)
    Xt = collect(0.0:0.005:1.0)
    Yt = f.(Xt)
    regress(Xo, Yo, Xt, Yt, polynomial, α, β, r)
end

max_stars_repo_name: gbarsih/Safe-Optimal-Rendezvous

max_issues_repo_name: gbarsih/Safe-Optimal-Rendezvous

max_forks_repo_name: gbarsih/Safe-Optimal-Rendezvous


text:
"""
# SeisSplit

SeisSplit is a package for measuring shear wave splitting using the
minimum eigenvalue method of Silver & Chan (1991), as modified by
Walsh et al. (2013).

It uses the [Seis.jl](https://github.com/anowacki/Seis.jl) package which
is an in-development community seismic analysis package in Julia.

Computation of splitting using the rotation-correlation method
(e.g., Bowman & Ando, 1987) is also possible.  This is computed
automatically by the main `splitting` function for quality control
purposes.

## Using

The `splitting` function performs shear wave splitting analysis and returns a `SeisSplit.Result`
type containing information about the analysis.  Provide two `Seis.Trace`s, and
optionally specify the maximum delay time and number of fast orientation and delay
time analysis points.

Example:

```julia
julia> using Seis, SeisSplit

julia> n, e = Seis.read_sac.(joinpath(dirname(pathof(SeisSplit)), "..", "test", "data", "wave.BH").*("N", "E"))
(Seis.Trace(.SWAV..BHN: delta=0.05, b=0.0, nsamples=1999), Seis.Trace(.SWAV..BHE: delta=0.05, b=0.0, nsamples=1999))

julia> s = splitting(n, e)
SeisSplit.Result{Float32,Array{Float32,1}}(-90.0:1.0:90.0, 0.0:0.10256410256410256:4.0, Float32[69.407 70.6024 … 64.1079 64.2652; 69.407 70.5508 … 64.6298 64.767; … ; 69.407 70.6525 … 63.5645 63.7426; 69.407 70.6024 … 64.1079 64.2652], Float32[8.79587 8.05464 … 12.038 11.9395; 8.79587 8.08652 … 11.7154 11.6291; … ; 8.79587 8.02363 … 12.3742 12.2627; 8.79587 8.05464 … 12.038 11.9395], 41.0, 0.5, 1.3333333333333333, 0.0, 10.455532f0, 1.1246407f0, Seis.Trace(.SWAV..BHN: delta=0.05, b=0.0, nsamples=1999), Seis.Trace(.SWAV..BHE: delta=0.05, b=0.0, nsamples=1999), 0.0f0, 99.9f0)

```
See the docstring for `splitting` for more information.

### Plotting results

You can create a diagnostic plot of a `SeisSplit.Result` by loading
[`Plots.jl`](https://github.com/JuliaPlots/Plots.jl) and calling `plot()` on the result:

```julia
julia> using Plots

julia> plot(s)
```

## References

- Bowman, J. R., Ando, M., 1987. Shear-wave splitting in the
  upper-mantle wedge above the Tonga subduction zone.  Geophys J R
  Astr Soc 88, 25–41.
  doi:[10.1111/j.1365-246X.1987.tb01367.x](https://doi.org/10.1111/j.1365-246X.1987.tb01367.x)
- Silver, P.G., Chan, W.W., 1991. Shear-wave splitting and subcontinental mantle
  deformation. J Geophys Res-Sol Ea 96, 16429–16454.
  doi:[10.1029/91JB00899](https://doi.org/10.1029/91JB00899)
- Walsh, E., Arnold, R., Savage, M.K., 2013. Silver and Chan revisited.
  Journal of Geophysical Research: Solid Earth 118, 5500–5515.
  doi:[10.1002/jgrb.50386](https://doi.org/10.1002/jgrb.50386)

"""
module SeisSplit

import DelimitedFiles, DSP, Printf, Statistics

import Distributions, FFTW
using LinearAlgebra
using StaticArrays

import Seis

export
    lam_ratio,
    lam1_ratio,
    lam2_ratio,
    quality,
    rotation_correlation,
    snr_restivo_helffrich,
    splitting

"Default number of δt points to search over"
const SPLIT_NDT = 41
"Default number of ϕ points to search over"
const SPLIT_NPHI = 181
"Default maximum δt (s)"
const SPLIT_DT_MAX = 4.0

"""
  Result{T,V}

Struct containing the results of shear wave splitting analysis.

Note that the two trace fields, `trace1` and `trace2` may not be in the same order
as originally given to the [`splitting`](@ref) function.  This is because
they are re-ordered so that `trace2` is clockwise of `trace1`.  Hence the
frame of reference in which this result is measured is determined by the
two traces' orientations.
"""
struct Result{T,V}
    "Range of values of ϕ over which we searched (° clockwise from north)"
    phi
    "Values of δt searched (s)"
    dt
    "Larger eigenvalue at each point in (ϕ,δt) space (normalised)"
    lam1
    "Smaller eigenvalue"
    lam2
    "Optimal fast orientation ϕ (° clockwise from north)"
    phi_best
    "Uncertainty in optimal ϕ (°)"
    dphi
    "Optimal δt (s)"
    dt_best
    "Uncertainty in optimal δt (s)"
    ddt
    "Source polarisation for optimal splitting (° clockwise from north)"
    spol
    "Uncertainty in source polarisation"
    dspol
    "Original input trace 1"
    trace1::Seis.Trace{T,V}
    "Original input trace 2"
    trace2::Seis.Trace{T,V}
    "Analysis window start time (s)"
    window_start
    "Analysis window end time (s)"
    window_end
    "Degrees of freedom"
    ndf
    "Cross correlation range of ϕ values (° clockwise from north)"
    xcorr_phi
    "Cross correlation range of δt values searched (s)"
    xcorr_dt
    "Cross correlation at each point in (ϕ,δt) space"
    xcorr_map
    "Maximum cross correlation point in ϕ (°)"
    xcorr_phi_best
    "Maximum cross correlation point in δt (s)"
    xcorr_dt_best
end

include("utils.jl")
include("io.jl")
include("qc.jl")
include("minimum_eigenvalue.jl")
include("rotation_correlation.jl")
include("plots.jl")

end # module

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
using SubmodularMaximization
using RosDataProcess
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics
using Printf
using Random
using HistogramFilters

close("all")

experiment_name = "entropy_by_solver"
data_folder = "./data"

reprocess = false

steps = 100
num_robots = 8:8:40

trials = 1:20

# We will drop a fraction of each trial so that the filters have time to
# converge to steady states
trial_steps = 20:steps

# Note: we will run trials in threads so the solvers do not have to be threaded

# Range limits
communication_range = 20
robot_target_range_limit = 12
range_limit_partitions = 4
range_limit_solver(x) = solve_communication_range_limit(x,
                                                        num_partitions=
                                                        range_limit_partitions,
                                                        communication_range=
                                                        communication_range)

solver_rounds = [2, 4, 8]
solvers = [(solve_random, Inf),
           (solve_myopic, Inf),
           map(num->(prob->solve_n_partitions(num, prob), Inf), solver_rounds)...,
           (range_limit_solver, robot_target_range_limit),
           (solve_sequential, Inf)]
solver_strings = ["random",
                  "myopic",
                  map(x->string("dist. ",x," rnds."), solver_rounds)...,
                  string("rng. lim. ", range_limit_partitions, " rnds."),
                  "sequential"]
solver_inds = 1:length(solvers)

all_tests = product(solver_inds, num_robots, trials)
all_configurations = product(solver_inds, num_robots)

function trial_fun(x)
  (solver_ind, num_robots, trial) = x
  solver, robot_target_range_limit = solvers[solver_ind]

  configs = MultiRobotTargetTrackingConfigs(num_robots,
                                            robot_target_range_limit=
                                            robot_target_range_limit
                                           )

  trial_data = target_tracking_experiment(steps=steps,
                                          num_robots=num_robots,
                                          configs=configs,
                                          solver=solver
                                         )
  (data=trial_data, configs=configs)
end
function print_summary(x)
  (solver_ind, num_robots, trial) = x
  _, robot_target_range_limit = solvers[solver_ind]

  println("Num. Robots: ", num_robots,
          " Solver: ", solver_strings[solver_ind],
          " Robot-Target Range Limit: ", robot_target_range_limit,
          " Trial: ", trial)
end

@time results = run_experiments(all_tests,
                             trial_fun=trial_fun,
                             print_summary=print_summary,
                             experiment_name=experiment_name,
                             data_folder=data_folder,
                             reprocess=reprocess
                            )

# Normalize by number of targets
entropies = Dict(map(all_tests) do key
                   key => map(results[key].data) do x
                     entropy(x.target_filters) / length(x.target_filters)
                   end
                 end)

# Concatenate all data for a set of trials
concatenated_entropy = Dict(map(all_configurations) do configuration
  trial_data = map(trials) do trial
    entropies[configuration..., trial][trial_steps]
  end

  configuration=>vcat(trial_data...)
end[:])

titles = map(all_configurations) do (solver_ind, num_robots)
  string(num_robots, "-robots ", solver_strings[solver_ind])
end[:]

boxplot([concatenated_entropy[c] for c in all_configurations][:],
        notch=false,
        vert=false)

title("Entropy distribution")
xlabel("Entropy per target (bits)")
yticks(1:length(titles), titles)

save_fig("fig", "entropy_by_solver")

# Plot means to get output in an easier format

figure()

# Plot everything in a given order to get the labels to line up nicely
solver_order = [3, 4, 5, 6, 1, 2, 7]

plot_args = Dict(3=>Dict(:color=>:k),
                 4=>Dict(:color=>:k, :linestyle=>"--"),
                 5=>Dict(:color=>:k, :linestyle=>":"),
                 6=>Dict(:color=>:k, :linestyle=>"-."))

for solver_ind in solver_order
  entropies = map(n->concatenated_entropy[solver_ind, n], num_robots)

  means = map(mean, entropies)
  label = solver_strings[solver_ind]

  # Process args
  kwargs = get(plot_args, solver_ind, Dict())

  plot(num_robots, means, label=label; kwargs...)
end

ylabel("Entropy per target (bits)")
xlabel("Num. robots")
legend(loc="upper right", ncol=2)
grid()

save_fig("fig", "entropy_by_solver_plot")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
module NoRainbows

import Base: print_stackframe, show_tuple_as_call, print_type_stacktrace, print_within_stacktrace, show_method_params, show_signature_function, show_full_backtrace, printstyled, with_output_color, showerror, show_backtrace
import Base.StackTraces: show_spec_linfo

using Core: MethodInstance, CodeInfo
using Base: StackFrame, stacktrace_expand_basepaths, stacktrace_contract_userdir, contractuser, unwrap_unionall, demangle_function_name, show_sym, parentmodule, empty_sym, text_colors, disable_text_style, is_exported_from_stdlib, process_backtrace, BIG_STACKTRACE_SIZE, stacktrace_linebreaks
using Base.StackTraces: top_level_scope_sym, is_top_level_frame

include("colors.jl")
include("pkg_colormap.jl")
include("frameline.jl")

@Base.kwdef struct ArgsTypesColorMap
    nonparamtypes::Color = Color()
    paramtypemain::Color = Color()
    paramtypearg::Color = Color()
    brackets::Color = Color()
end
const ARGTYPE_COLORS = Ref(ArgsTypesColorMap())
function set_argtype_map(;
    default = nothing,
    nonparamtypes = nothing,
    paramtypemain = nothing,
    paramtypearg = nothing,
    brackets = nothing,
)
    c = ARGTYPE_COLORS[]
    ARGTYPE_COLORS[] = ArgsTypesColorMap(;
        nonparamtypes = Color(coal(nonparamtypes, default, c.nonparamtypes)),
        paramtypemain = Color(coal(paramtypemain, default, c.paramtypemain)),
        paramtypearg = Color(coal(paramtypearg, default, c.paramtypearg)),
        brackets = Color(coal(brackets, default, c.brackets))
    )

    return nothing
end

@Base.kwdef struct TrackColorMap
    tin::Color = Color()
    tout::Color = Color()
end
const TRACK_COLORS = Ref(TrackColorMap())
function set_track_map(;
    default = nothing,
    trackin = nothing,
    trackout = nothing
)
    c = TRACK_COLORS[]
    TRACK_COLORS[] = TrackColorMap(;
        tin = Color(coal(trackin, default, c.tin)),
        tout = Color(coal(trackout, default, c.tout))
    )

    return nothing
end

@Base.kwdef struct MethodParamsColorMap
    brackets::Color = Color()
    comma::Color = Color()
    wher::Color = Color()
    types::Color = Color()
end
const METHODPARAMS_COLORS = Ref(MethodParamsColorMap())
function set_methodparams_map(;
    default = nothing,
    brackets = nothing,
    wher = nothing,
    types = nothing,
)
    c = METHODPARAMS_COLORS[]
    METHODPARAMS_COLORS[] = MethodParamsColorMap(;
        brackets = Color(coal(brackets, default, c.brackets)),
        wher = Color(coal(wher, default, c.wher)),
        types = Color(coal(types, default, c.types)),
    )

    return nothing
end

@Base.kwdef struct TupleCallColorMap
    tuple::Color = Color()
    brackets::Color = Color()
    doublecolon::Color = Color()
    comma::Color = Color()
    semicolon::Color = Color()
    args::Color = Color()
    kwargs::Color = Color()
end
const TUPLE_CALL_COLORS = Ref(TupleCallColorMap())
function set_tuplecall_map(;
    default = nothing,
    tuple = nothing,
    brackets = nothing,
    doublecolon = nothing,
    comma = nothing,
    semicolon = nothing,
    args = nothing,
    kwargs = nothing
)
    c = TUPLE_CALL_COLORS[]
    TUPLE_CALL_COLORS[] = TupleCallColorMap(
        tuple = Color(coal(tuple, default, c.tuple)),
        brackets = Color(coal(brackets, default, c.brackets)),
        doublecolon = Color(coal(doublecolon, default, c.doublecolon)),
        comma = Color(coal(comma, default, c.comma)),
        semicolon = Color(coal(semicolon, default, c.semicolon)),
        args = Color(coal(args, default, c.args)),
        kwargs = Color(coal(kwargs, default, c.kwargs))
    )

    return nothing
end

@Base.kwdef struct SpecLinfoColorMap
    ipx::Color = Color()
    toplevel::Color = Color()
    repeats::Color = Color()
    modulename::Color = Color()
    funcname::Color = Color()
    brackets::Color = Color()
    wrapper::Color = Color()
    functor::Color = Color()
end
const SPECLINFO_COLORS = Ref(SpecLinfoColorMap())
function set_speclinfo_map(;
    default = nothing,
    ipx = nothing,
    toplevel = nothing,
    repeats = nothing,
    modulename = nothing,
    funcname = nothing,
    brackets = nothing,
    wrapper = nothing,
    functor = nothing,
)
    c = SPECLINFO_COLORS[]
    SPECLINFO_COLORS[] = SpecLinfoColorMap(
        ipx = Color(coal(ipx, default, c.ipx)),
        toplevel = Color(coal(toplevel, default, c.toplevel)),
        repeats = Color(coal(repeats, default, c.repeats)),
        modulename = Color(coal(modulename, default, c.modulename)),
        funcname = Color(coal(funcname, default, c.funcname)),
        brackets = Color(coal(brackets, default, c.brackets)),
        wrapper = Color(coal(wrapper, default, c.wrapper)),
        functor = Color(coal(functor, default, c.functor))
    )

    return nothing
end

@Base.kwdef struct FrameNumberMap
    frameno::Color = Color()
    brackets::Color = Color()
end
const FRAME_NUMBER = Ref(FrameNumberMap())
"""
    set_framenumber(;
        default = nothing,
        frameno = nothing,
        brackets = nothing
    )

Set color attributes of the frame number where frameline was generated.
"""
function set_framenumber(;
        default = nothing,
        frameno = nothing,
        brackets = nothing,
    )
    c = FRAME_NUMBER[]
    FRAME_NUMBER[] = FrameNumberMap(
        frameno = Color(coal(frameno, default, c.frameno)),
        brackets = Color(coal(brackets, default, c.brackets))
    )

    return nothing
end

@Base.kwdef struct FilepathColorMap
    filepath::Color = Color()
    colon::Color = Color()
    lineno::Color = Color()
    inlined::Color = Color()
end
const FILEPATH_COLORS = Ref(FilepathColorMap())
"""
    set_filepath(;
        default = nothing,
        filepath = nothing,
        colon = nothing,
        lineno = nothing,
        inlined = nothing
    )

Set color attributes of the filepath which generate corresponding frameline.
"""
function set_filepath(;
    default = nothing,
    filepath = nothing,
    colon = nothing,
    lineno = nothing,
    inlined = nothing
)
    c = FILEPATH_COLORS[]
    filepath = Color(coal(filepath, default, c.filepath))
    colon = Color(coal(colon, default, c.colon))
    lineno = Color(coal(lineno, default, c.lineno))
    inlined = Color(coal(inlined, default, c.inlined))

    FILEPATH_COLORS[] = FilepathColorMap(
        filepath = filepath,
        colon = colon,
        lineno = lineno,
        inlined = inlined
    )
end

@Base.kwdef struct GlobalOptions
    reverse::Bool = false
    align_numbering::Bool = false
    stackcolor::Color = Color()
end
const GLOBAL = Ref(GlobalOptions())
function set_globals(;
    reverse = nothing,
    stackcolor = nothing,
    align_numbering = nothing,
)
    c = GLOBAL[]
    col = Color(stackcolor === nothing ? c.stackcolor : stackcolor)
    rev = reverse === nothing ? c.reverse : reverse
    align = align_numbering === nothing ? c.align_numbering : align_numbering

    GLOBAL[] = GlobalOptions(rev, align, col)

    return nothing
end

const STACKTRACE_MODULECOLORS = [Color(), Color(), Color(), Color()]
const STACKTRACE_FIXEDCOLORS = IdDict(Base => Color(), Core => Color())
const TRACK_MODUL = Module[]
function track_modules(moduls...)
    empty!(TRACK_MODUL)
    for modul in moduls
        push!(TRACK_MODUL, modul)
    end

    return nothing
end

get_module_fixed() = STACKTRACE_FIXEDCOLORS
function set_module_fixed(d)
    empty!(STACKTRACE_FIXEDCOLORS)
    for (k, v) in d
        STACKTRACE_FIXEDCOLORS[k] = Color(v)
    end
end

get_module_rotating() = STACKTRACE_MODULECOLORS
function set_module_rotating(v)
    empty!(STACKTRACE_MODULECOLORS)
    for x in v
        push!(STACKTRACE_MODULECOLORS, Color(x))
    end

    nothing
end

function set_solarized()
    set_filepath(default = "normal", lineno = "cyan")
    set_framenumber(frameno = "cyan")
    set_speclinfo_map(default = "blue", repeats = "normal")
    set_tuplecall_map(default = "normal", doublecolon = "green")
    set_argtype_map(default = "yellow", brackets = "normal")

    set_track_map(trackin = "cyan:reverse", trackout = "red:reverse")
    set_globals(stackcolor = "underline")

    format_frameline("{frameno:lalign} {module}{filepath}\n  {function}")
end

function set_theme(theme)
    if theme == "solarized"
        set_solarized()
    elseif theme == "solarized_reversed"
        set_solarized()
        set_globals(reverse = true, align_numbering = true)
    else
        @warn "Unknown theme $theme"
    end
    return nothing
end

printstyled(io::IO, color::Color, msg...) = printstyled(io, msg...; color = color.color, bold = color.bold, underline = color.underline, blink = color.blink, reverse = color.reverse, hidden = color.hidden)

showstyled(io::IO, color::Color, msg...) = with_output_color(show, color.color, io, msg...; bold = color.bold, underline = color.underline, blink = color.blink, reverse = color.reverse, hidden = color.hidden)

function showerror(io::IO, ex::MethodError, bt; backtrace=true)
    try
        showerror(io, ex)
    finally
        f = ex.f
        ft = typeof(f)
        backtrace && show_backtrace(io, bt; modul = ft.name.module)
    end
end

function show_backtrace(io::IO, t::Vector; modul = nothing)
    if haskey(io, :last_shown_line_infos)
        empty!(io[:last_shown_line_infos])
    end

    # t is a pre-processed backtrace (ref #12856)
    if t isa Vector{Any}
        filtered = t
    else
        filtered = process_backtrace(t)
    end
    isempty(filtered) && return

    if length(filtered) == 1 && StackTraces.is_top_level_frame(filtered[1][1])
        f = filtered[1][1]::StackFrame
        if f.line == 0 && f.file === Symbol("")
            # don't show a single top-level frame with no location info
            return
        end
    end

    if length(filtered) > BIG_STACKTRACE_SIZE
        show_reduced_backtrace(IOContext(io, :backtrace => true), filtered)
        return
    end

    try invokelatest(update_stackframes_callback[], filtered) catch end
    # process_backtrace returns a Vector{Tuple{Frame, Int}}
    frames = map(x->first(x)::StackFrame, filtered)
    show_full_backtrace(io, frames; print_linebreaks = stacktrace_linebreaks(), modul = modul)
    return
end

function show_full_backtrace(io::IO, trace::Vector; print_linebreaks::Bool = true, modul = nothing)
    n = length(trace)
    ndigits_max = ndigits(n)

    modulecolordict = copy(STACKTRACE_FIXEDCOLORS)
    modulecolorcycler = Iterators.Stateful(Iterators.cycle(STACKTRACE_MODULECOLORS))

    printstyled(io, GLOBAL[].stackcolor, "\nStacktrace:\n")

    moduls = similar(trace, Int)
    state = 0 # are we inside tracked module or not
    tracks = modul === nothing ? TRACK_MODUL : [modul]
    for (i, frame) in pairs(trace)
        m = parentmodule(frame)
        if i == 1
            moduls[i] = m in tracks ? -1 : 0
            state = moduls[i] == 0 ? 0 : 1
        else
            if state == 0
                if m in tracks
                    moduls[i] = -1
                    state = 1
                else
                    moduls[i] = 0
                end
            else
                if m in tracks
                    moduls[i] = 0
                else
                    moduls[i - 1] = 1
                    moduls[i] = 0
                    state = 0
                end
            end
        end
    end

    if GLOBAL[].reverse
        trace = reverse(trace)
        moduls = reverse(moduls)
    end

    for (i, (frame, modultrack)) in enumerate(zip(trace, moduls))
        if (GLOBAL[].align_numbering && GLOBAL[].reverse) || (!GLOBAL[].align_numbering && !GLOBAL[].reverse)
            print_stackframe(io, i, frame, 1, ndigits_max, modulecolordict, modulecolorcycler, modultrack)
        else
            print_stackframe(io, n - i + 1, frame, 1, ndigits_max, modulecolordict, modulecolorcycler, modultrack)
        end
        if i < n
            println(io)
        end
    end
end

# Print a stack frame where the module color is determined by looking up the parent module in
# `modulecolordict`. If the module does not have a color, yet, a new one can be drawn
# from `modulecolorcycler`.
function print_stackframe(io, i, frame::StackFrame, n::Int, digit_align_width, modulecolordict, modulecolorcycler, modultrack)
    m = Base.parentmodule(frame)
    if m !== nothing
        while parentmodule(m) !== m
            pm = parentmodule(m)
            pm == Main && break
            m = pm
        end
        if !haskey(modulecolordict, m)
            modulecolordict[m] = popfirst!(modulecolorcycler)
        end
        modulecolor = modulecolordict[m]
    else
        modulecolor = Color()
    end
    print_stackframe(io, i, frame, n, digit_align_width, modulecolor, modultrack)
end

function print_stackframe(io, i, frame::StackFrame, n::Int, digit_align_width, modulecolor, modultrack)
    file, line = string(frame.file), frame.line
    stacktrace_expand_basepaths() && (file = something(find_source_file(file), file))
    stacktrace_contract_userdir() && (file = contractuser(file))

    # Used by the REPL to make it possible to open
    # the location of a stackframe/method in the editor.
    if haskey(io, :last_shown_line_infos)
        push!(io[:last_shown_line_infos], (string(frame.file), frame.line))
    end

    inlined = getfield(frame, :inlined)
    modul = parentmodule(frame)

    for token in FRAME_LINE
        if !token.special
            if token.align == 0
                printstyled(io, token.color, token.val)
            else
                if token.align == -1
                    printstyled(io, token.color, rpad(token.val, digit_align_width + 2))
                else
                    printstyled(io, token.color, lpad(token.val, digit_align_width + 2))
                end
            end
        elseif token.val == "frameno"
            # frame number
            frameno = if modultrack == 1
                TRACK_COLORS[].tin
            elseif modultrack == -1
                TRACK_COLORS[].tout
            else
                FRAME_NUMBER[].frameno
            end
            brackno = if modultrack == 1
                TRACK_COLORS[].tin
            elseif modultrack == -1
                TRACK_COLORS[].tout
            else
                FRAME_NUMBER[].brackets
            end

            num = string(i)

            if token.align == 1
                k = digit_align_width - length(num)
                if k > 0
                    print(io, " "^k)
                end
            end
            printstyled(io, brackno, "[")
            printstyled(io, frameno, num)
            printstyled(io, brackno, "]")
            if token.align == -1
                k = digit_align_width - length(num)
                if k > 0
                    print(io, " "^k)
                end
            end
        elseif token.val == "function"
            show_spec_linfo(IOContext(io, :backtrace=>true), frame)
            if n > 1
                printstyled(io, SPECLINFO_COLORS[].repeats, " (repeats $n times)")
            end
        elseif token.val == "filepath"
            # filepath
            pathparts = splitpath(file)
            folderparts = pathparts[1:end-1]
            if !isempty(folderparts)
                printstyled(io, FILEPATH_COLORS[].filepath, joinpath(folderparts...) * (Sys.iswindows() ? "\\" : "/"))
            end

            # filename, separator, line
            # use escape codes for formatting, printstyled can't do underlined and color
            # codes are bright black (90) and underlined (4)
            printstyled(io, FILEPATH_COLORS[].filepath, pathparts[end])
            printstyled(io, FILEPATH_COLORS[].colon, ":")
            printstyled(io, FILEPATH_COLORS[].lineno, line)

            # inlined
            printstyled(io, FILEPATH_COLORS[].inlined, inlined ? " [inlined]" : "")
        elseif token.val == "module"
            # module
            if modul !== nothing
                printstyled(io, modulecolor, modul, " ")
            end
        end
    end
end

function show_spec_linfo(io::IO, frame::StackFrame)
    linfo = frame.linfo
    if linfo === nothing
        if frame.func === empty_sym
            printstyled(io, SPECLINFO_COLORS[].ipx, "ip:0x", string(frame.pointer, base=16))
        elseif frame.func === top_level_scope_sym
            printstyled(io, SPECLINFO_COLORS[].toplevel, "top-level scope")
        else
            print_within_stacktrace(io, SPECLINFO_COLORS[].funcname, Base.demangle_function_name(string(frame.func)))
        end
    elseif linfo isa MethodInstance
        def = linfo.def
        if isa(def, Method)
            sig = linfo.specTypes
            argnames = Base.method_argnames(def)
            if def.nkw > 0
                # rearrange call kw_impl(kw_args..., func, pos_args...) to func(pos_args...)
                kwarg_types = Any[ fieldtype(sig, i) for i = 2:(1+def.nkw) ]
                uw = Base.unwrap_unionall(sig)::DataType
                pos_sig = Base.rewrap_unionall(Tuple{uw.parameters[(def.nkw+2):end]...}, sig)
                kwnames = argnames[2:(def.nkw+1)]
                for i = 1:length(kwnames)
                    str = string(kwnames[i])::String
                    if endswith(str, "...")
                        kwnames[i] = Symbol(str[1:end-3])
                    end
                end
                show_tuple_as_call(io, def.name, pos_sig;
                                        demangle=true,
                                        kwargs=zip(kwnames, kwarg_types),
                                        argnames=argnames[def.nkw+2:end])
            else
                show_tuple_as_call(io, def.name, sig; demangle=true, argnames)
            end
        else
            show_mi(io, linfo, true)
        end
    elseif linfo isa CodeInfo
        printstyled(io, SPECLINFO_COLORS[].toplevel, "top-level scope")
    end
end

function show_tuple_as_call(io::IO, name::Symbol, sig::Type;
                            demangle=false, kwargs=nothing, argnames=nothing,
                            qualified=false, hasfirst=true)
    # print a method signature tuple for a lambda definition
    if sig === Tuple
        printstyled(io, TUPLE_CALL_COLORS[].tuple, demangle ? demangle_function_name(name) : name, "(...)")
        return
    end
    tv = Any[]
    env_io = io
    while isa(sig, UnionAll)
        push!(tv, sig.var)
        env_io = IOContext(env_io, :unionall_env => sig.var)
        sig = sig.body
    end
    n = 1
    sig = (sig::DataType).parameters
    if hasfirst
        show_signature_function(env_io, sig[1], demangle, "", false, qualified)
        n += 1
    end
    first = true
    print_within_stacktrace(io, TUPLE_CALL_COLORS[].brackets, "(")
    show_argnames = argnames !== nothing && length(argnames) == length(sig)
    for i = n:length(sig)  # fixme (iter): `eachindex` with offset?
        first || printstyled(io, TUPLE_CALL_COLORS[].comma, ", ")
        first = false
        if show_argnames
            print_within_stacktrace(io, TUPLE_CALL_COLORS[].args, argnames[i])
        end
        printstyled(io, TUPLE_CALL_COLORS[].doublecolon, "::")
        print_type_stacktrace(env_io, sig[i])
    end
    if kwargs !== nothing
        printstyled(io, TUPLE_CALL_COLORS[].semicolon, "; ")
        first = true
        for (k, t) in kwargs
            first || printstyled(io, TUPLE_CALL_COLORS[].comma, ", ")
            first = false
            print_within_stacktrace(io, TUPLE_CALL_COLORS[].kwargs, k)
            printstyled(io, TUPLE_CALL_COLORS[].doublecolon, "::")
            print_type_stacktrace(io, t)
        end
    end
    print_within_stacktrace(io, TUPLE_CALL_COLORS[].brackets, ")")
    show_method_params(io, tv)
    nothing
end

# show the called object in a signature, given its type `ft`
# `io` should contain the UnionAll env of the signature
function show_signature_function(io::IO, @nospecialize(ft), demangle=false, fargname="", html=false, qualified=false)
    uw = unwrap_unionall(ft)
    if ft <: Function && isa(uw, DataType) && isempty(uw.parameters) &&
        isdefined(uw.name.module, uw.name.mt.name) &&
        ft == typeof(getfield(uw.name.module, uw.name.mt.name))
        if qualified && !is_exported_from_stdlib(uw.name.mt.name, uw.name.module) && uw.name.module !== Main
            print_within_stacktrace(io, SPECLINFO_COLORS[].modulename, uw.name.module)
        end
        s = sprint(show_sym, (demangle ? demangle_function_name : identity)(uw.name.mt.name), context=io)
        print_within_stacktrace(io, SPECLINFO_COLORS[].funcname, s)
    elseif isa(ft, DataType) && ft.name === Type.body.name &&
        (f = ft.parameters[1]; !isa(f, TypeVar))
        uwf = unwrap_unionall(f)
        parens = isa(f, UnionAll) && !(isa(uwf, DataType) && f === uwf.name.wrapper)
        parens && printstyled(io, SPECLINFO_COLORS[].brackets, "(")
        showstyled(io, SPECLINFO_COLORS[].wrapper, f)
        parens && printstyled(io, SPECLINFO_COLORS[].brackets, ")")
    else
        if html
            print(io, "($fargname::<b>", ft, "</b>)")
        else
            print_within_stacktrace(io, SPECLINFO_COLORS[].functor, "($fargname::", ft, ")")
        end
    end
    nothing
end

function print_type_stacktrace(io, type; color=:normal)
    str = sprint(show, type, context=io)
    i = findfirst('{', str)
    if isnothing(i) || !get(io, :backtrace, false)::Bool
        printstyled(io, ARGTYPE_COLORS[].nonparamtypes, str)
    else
        printstyled(io, ARGTYPE_COLORS[].paramtypemain, str[1:prevind(str,i)])
        printstyled(io, ARGTYPE_COLORS[].brackets, "{")
        printstyled(io, ARGTYPE_COLORS[].paramtypearg, str[i+1:end-1])
        printstyled(io, ARGTYPE_COLORS[].brackets, "}")
    end
end

function print_within_stacktrace(io, color::Color, s...)
    if get(io, :backtrace, false)::Bool
        printstyled(io, color, s...)
    else
        print(io, s...)
    end
end

function show_method_params(io::IO, tv)
    if !isempty(tv)
        printstyled(io, METHODPARAMS_COLORS[].wher, " where ")
        if length(tv) == 1
            showstyled(io, METHODPARAMS_COLORS[].types, tv[1])
        else
            printstyled(io, METHODPARAMS_COLORS[].brackets, "{")
            for i = 1:length(tv)
                if i > 1
                    printstyled(io, METHODPARAMS_COLORS[].comma, ", ")
                end
                x = tv[i]
                showstyled(io, METHODPARAMS_COLORS[].types, x)
                io = IOContext(io, :unionall_env => x)
            end
            printstyled(io, METHODPARAMS_COLORS[].brackets, "}")
        end
    end
end

function with_output_color(@nospecialize(f::Function), color::Union{Int, Symbol}, io::IO, args...;
        bold::Bool = false, underline::Bool = false, blink::Bool = false,
        reverse::Bool = false, hidden::Bool = false)
    buf = IOBuffer()
    iscolor = get(io, :color, false)::Bool
    try f(IOContext(buf, io), args...)
    finally
        str = String(take!(buf))
        if !iscolor
            print(io, str)
        else
            bold && color === :bold && (color = :nothing)
            underline && color === :underline && (color = :nothing)
            blink && color === :blink && (color = :nothing)
            reverse && color === :reverse && (color = :nothing)
            hidden && color === :hidden && (color = :nothing)
            enable_ansi  = get(text_colors, color, text_colors[:default]) *
                               (bold ? text_colors[:bold] : "") *
                               (underline ? text_colors[:underline] : "") *
                               (blink ? text_colors[:blink] : "") *
                               (reverse ? text_colors[:reverse] : "") *
                               (hidden ? text_colors[:hidden] : "")

            disable_ansi = (hidden ? disable_text_style[:hidden] : "") *
                           (reverse ? disable_text_style[:reverse] : "") *
                           (blink ? disable_text_style[:blink] : "") *
                           (underline ? disable_text_style[:underline] : "") *
                           (bold ? disable_text_style[:bold] : "") *
                               get(disable_text_style, color, text_colors[:default])
            first = true
            for line in split(str, '\n')
                first || print(buf, '\n')
                first = false
                isempty(line) && continue
                print(buf, enable_ansi, line, disable_ansi)
            end
            print(io, String(take!(buf)))
        end
    end
end

"""
    printstyled([io], xs...; bold::Bool=false, underline::Bool=false, blink::Bool=false, reverse::Bool=false, hidden::Bool=false, color::Union{Symbol,Int}=:normal)

Print `xs` in a color specified as a symbol or integer, optionally in bold.

`color` may take any of the values $(Base.available_text_colors_docstring)
or an integer between 0 and 255 inclusive. Note that not all terminals support 256 colors.
If the keyword `bold` is given as `true`, the result will be printed in bold.
If the keyword `underline` is given as `true`, the result will be printed underlined.
If the keyword `blink` is given as `true`, the result will blink.
If the keyword `reverse` is given as `true`, the result will have foreground and background colors inversed.
If the keyword `hidden` is given as `true`, the result will be hidden.
Keywords can be given in any combination.
"""
printstyled(io::IO, msg...; bold::Bool=false, underline::Bool=false, blink::Bool=false, reverse::Bool=false, hidden::Bool=false, color::Union{Int,Symbol}=:normal) =
    with_output_color(print, color, io, msg...; bold=bold, underline=underline, blink=blink, reverse=reverse, hidden=hidden)
printstyled(msg...; bold::Bool=false, underline::Bool=false, blink::Bool=false, reverse::Bool=false, hidden::Bool=false, color::Union{Int,Symbol}=:normal) =
    printstyled(stdout, msg...; bold=bold, underline=underline, blink=blink, reverse=reverse, hidden=hidden, color=color)

end # module


text:
# --------------------------------------------------------------------
# plain TXT

function test_reference(file::File{format"TXT"}, actual)
    _test_reference(Diff(), file, string(actual))
end

function test_reference(file::File{format"TXT"}, actual::AbstractArray{<:AbstractString})
    str = join(actual, '\n')
    _test_reference(Diff(), file, str)
end

# ---------------------------------
# Image

function test_reference(file::File, actual::AbstractArray{<:Colorant}; sigma=ones(length(axes(actual))), eps=0.01)
    _test_reference(BeforeAfterImage(), file, actual) do reference, actual
        try
            Images.@test_approx_eq_sigma_eps(reference, actual, sigma, eps)
            return true
        catch err
            if err isa ErrorException
                return false
            else
                rethrow()
            end
        end
    end
end

# Image as txt using ImageInTerminal
function test_reference(file::File{format"TXT"}, actual::AbstractArray{<:Colorant}; size = (20,40))
    strs = @withcolor ImageInTerminal.encodeimg(ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), actual, size...)[1]
    str = join(strs,'\n')
    _test_reference(BeforeAfterFull(), file, str)
end

# --------------------------------------------------------------------
# SHA as string

function test_reference(file::File{format"SHA256"}, actual)
    test_reference(file, string(actual))
end

function test_reference(file::File{format"SHA256"}, actual::Union{AbstractString,Vector{UInt8}})
    str = bytes2hex(sha256(actual))
    _test_reference(BeforeAfterFull(), file, str)
end

function test_reference(file::File{format"SHA256"}, actual::AbstractArray{<:Colorant})
    size_str = bytes2hex(sha256(reinterpret(UInt8,[map(Int64,size(actual))...])))
    img_str = bytes2hex(sha256(reinterpret(UInt8,vec(rawview(channelview(actual))))))
    _test_reference(BeforeAfterFull(), file, size_str * img_str)
end

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

# Fallback
function test_reference(file::File, actual)
    if actual isa AbstractString
        # we don't use dispatch for this as it is very ambiguous
        # specialization will remove this conditional regardless

        _test_reference(Diff(), file, actual)
    else
        _test_reference(BeforeAfterLimited(), file, actual)
    end
end


text:
using BoardGames

struct TicTacToeBoard
    v::Vector{Int}
    turn::Int
    winner::Int
end

struct TicTacToe <: Game{TicTacToeBoard, Int} end

function TicTacToeBoard()
    v = zeros(Int, (9))
    turn = 1
    winner = 0

    TicTacToeBoard(v, turn, winner)
end

function BoardGames.getmoves(board::TicTacToeBoard)
    if board.winner == 0
        (1:9)[board.v .== 0]
    else
        return Int[]
    end
end

function Base.copy(board::TicTacToeBoard)
    TicTacToeBoard(copy(board.v), board.turn, board.winner)
end

const lines = [[1,2,3], [4,5,6], [7,8,9], [1,4,7], [2,5,8], [3,6,9], [1,5,9], [3,5,7]]

function BoardGames.play(board::TicTacToeBoard, move::Int)
    v = copy(board.v)
    v[move] = board.turn
    turn = -board.turn

    winner = 0
    for l in lines
        if all(v[l] .== board.turn)
            winner = board.turn
        end
    end

    TicTacToeBoard(v, turn, winner)
end

function BoardGames.playerturn(board::TicTacToeBoard)
    if board.turn == 1
        return 1
    else
        return 2
    end
end


text:
# This file was generated, do not modify it. # hide
etime = @elapsed eknns, edists = allknn(ExhaustiveSearch(; db, dist), k; parallel=true)
nothing # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using Distributions

f3_size = 0
deceptiveN_size = 0
deceptiveN_nbits = 0

function change_deceptiveN_size( nsize :: Int, nbits :: Int )
    global deceptiveN_size  = nsize
    global deceptiveN_nbits = nbits
end

function change_f3_size( nsize :: Int )
    global f3_size = nsize
end

function eval_f3( vet )
    if     reduce(&, true, vet .== [false, false, false])
        return 28
    elseif reduce(&, true, vet .== [false, false, true ])
        return 26
    elseif reduce(&, true, vet .== [false, true , false])
        return 22
    elseif reduce(&, true, vet .== [true , false, false])
        return 14
    elseif reduce(&, true, vet .== [true , true , true ])
        return 30
    else
        return 0
    end

    return 0
end

function objf_f3( ind )
    obj = 0
    gens = ind

    for i in 1:f3_size
        obj += eval_f3(gens[(i-1)*3+1:i*3])
    end

    return obj
end

function objf_f3s( ind )
    obj = 0
    gens = [x for x in ind]

    for i in 1:f3_size
        obj += eval_f3([gens[i], gens[i + 10], gens[i+20]])
    end

    return obj
end

function objf_deceptiveN( ind )
    obj = 0
    gens = [x for x in ind]

    for i in 1:deceptiveN_size
        v = sum(gens[(i-1) * deceptiveN_nbits + 1: i * deceptiveN_nbits])

        if v == 0
            obj += deceptiveN_nbits + 1
        else
            obj += v
        end
    end

    return obj
end

function objf_rosen( ind )
    d = length(ind)

    xi    = ind[1:d-1]
    xnext = ind[2:d]

    s = mapfoldr(x -> 100.0 * (x[2] - x[1]^2.0)^2.0 + (x[1] - 1.0)^2.0, +, collect(zip(xi, xnext)))

    return s
end

function objf_sphere( ind )
    obj = 0.0

    s = length(ind)

    for i in 1:s
        obj += ind[i] ^ 2.0
    end

    return obj
end

###

function evaluate(vars, formula)
    #=return objf_f3( vars )=#
    #=return objf_f3s( vars )=#
    #=return objf_deceptiveN( vars )=#
    #=return objf_sphere( vars )=#
    return objf_rosen( vars )
end

function energy(vars, formula)
    q = evaluate(vars, formula)
    p = 160

    return q

    #=return 1.0 - q/p=#
    #=return p - q=#
end

function neighbour(vars)
    nv    = copy(vars)
    p       = rand(1:length(vars))

    #=nv[p] = !nv[p]=#
    nv[p] += rand(Normal(0, .25))

    for i in 1:length(vars)
        nv[p] = clamp(nv[p], -5.12, 5.12)
    end

    return nv
end


text:
@testset "Test domain functions" begin
    dist, bearing = sail_route.haversine(-88.67, 36.12, -118.40, 33.94)
    @test dist ≈ 1462.22 atol = 0.01
    @test bearing ≈ 273.662 atol = 0.01
end

@testset "Test optimum route algorithm" begin
    path = ENV["HOME"]*"/Documents/sail_route.jl/src/data/first40_orgi.csv"
    twa, tws, perf = load_file(path)
    polar = sail_route.setup_interpolation(tws, twa, perf)
    lon1 = 0.0
    lon2 = 1.0
    lat1 = 5.0
    lat2 = 5.0
    d_an, b_an = sail_route.haversine(lon1, lat1, lon2, lat2)
    time_an = d/sail_route.perf_interp(polar, min_angle(wd_int, b), ws_int)
    start = Dates.DateTime(2000, 1, 1, 0, 0, 0)
    x_nodes = 5
    y_nodes = 5
    y_dist = 100.0
    wisp, widi, cusp, cudi = sail_route.sample_weather()
    st, ets = sail_route.route_solve(lon1, lon2, lat1, lat2, x_nodes, y_nodes, y_dist, wisp, widi, cusp, cudi, polar, start)
    @test isapprox(st, time_an; rtol=0.0001)
end


text:
import Random

matrixrandomizer(com::ComMatrix, rng = Xoroshiro128Plus();
                method::matrixrandomizations = curveball) = "Only defined for Boolean Assemblages"
matrixrandomizer(asm::SEAssemblage, rng = Xoroshiro128Plus();
                method::matrixrandomizations = curveball) = "Only defined for Boolean Assemblages"

function matrixrandomizer(asm::S, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball) where S <: SEAssemblage{Bool}
    as = copy(asm)
    ret = MatrixGenerator{typeof(rng), typeof(as)}(as, method, rng)
    dropzeros!(ret.m.occ.commatrix.occurrences)
    ret
end

function matrixrandomizer(com::C, rng = Xoroshiro128Plus();
                        method::matrixrandomizations = curveball) where C <: ComMatrix{Bool}
    cm = copy(com)
    ret = MatrixGenerator{typeof(rng), typeof(cm)}(cm, method, rng)
    dropzeros!(ret.m.occurrences)
    ret
end

Random.rand(r::MatrixGenerator{R, <:Assemblage}) where R = copy(Random.rand!(r))
function Random.rand!(r::MatrixGenerator{R, A}) where {R} where {A <: Assemblage}
    RandomBooleanMatrices._curveball!(r.m.occ.commatrix.occurrences)
    r.m
end

Random.rand(r::MatrixGenerator{<:ComMatrix}) where R = copy(Random.rand!(r))
function Random.rand!(r::MatrixGenerator{C}) where {C <: ComMatrix}
    RandomBooleanMatrices._curveball!(r.m.occurrences)
    r.m
end


text:
# Range of frequencies
function range!(freq::AbstractArray, sampling_rate)
    n = size(freq, 1)
    delta = 1.0 / sampling_rate
    total_duration = delta * n

    for i in 1:n
        freq[i] = i / total_duration
    end

    return nothing
end

@testset "Cosine Taper Test" begin
    n = 18000
    data = ones(Float64, n);
    α = 0.05

    cosine_taper!(data, n, α)

    for i in 1:n
        @test data[i] <= 1.0
        @test data[i] >= 0.0
    end
end

@testset "SAC Taper Test" begin
    n = 18000
    data = ones(Float64, n)
    freq = Array{Float64}(undef, n)
    sampling_rate = 20.0
    f1, f2, f3, f4 = 1.0, 2.0, 8.0, 9.0

    range!(freq, sampling_rate)
    taper = sac_cosine_taper(freq, f1, f2, f3, f4, sampling_rate)
    for i in 1:n
        data[i] *= taper[i]
    end

    for i in 1:n
        @test data[i] <= 1.0
        @test data[i] >= 0.0
    end
end


text:
#using DeformableRegistration.ImageProcessing
using DeformableRegistration: Regularizer, Transformation
 using Base.Test

# @testset "check if constant deformation yields zero" begin
#     omega = [0.0,1.0,0.0,1.0];
#     m = [5; 7];
#     h = [1.0,1.0];
#     u = ones(size(getStaggeredGrid(omega,m)))
#     B = createElasticOperatorStaggered(h,m)
#     S, dS, d2S = regularizer(u, B)
#     @test_approx_eq(S[1], 0)
#     @test sum(dS) == 0.0
# end
#
# @testset "check if random deformation yields nonzero result" begin
#     omega = [0.0,1.0,0.0,1.0];
#     m = [5; 7];
#     h = [1.0,1.0];
#     u = rand(size(getStaggeredGrid(omega,m)))
#     B = createElasticOperatorStaggered(h,m)
#     S, dS, d2S = regularizer(u, B)
#     @test S[1] > 0
#     @test (dS'*dS)[1] > 0
# end

@testset "Regularizer" begin
@testset "check if constant deformation yields zero" begin
    m = [5; 7];
    h = [1.0,1.0];
    u = ones(2*m[1],m[2])[:]
    B = createElasticOperatorCentered(h,m)
    S, dS, d2S = regularizer(u, B)
    @test  S[1] ≈ 0
    @test sum(dS) == 0.0
end

@testset "check if random deformation yields nonzero result" begin
    omega = [0.0,1.0,0.0,1.0];
    m = [5; 7];
    h = [1.0,1.0];
    u = rand(2*m[1],m[2])[:]
    B = createElasticOperatorCentered(h,m)
    S, dS, d2S = regularizer(u, B)
    @test S[1] > 0
    @test (dS'*dS)[1] > 0
end

@testset "check if constant deformation yields zero" begin
    omega = [0.0,1.0,0.0,1.0];
    m = [5; 7];
    h = [1.0,1.0];
    u = ones(2*m[1],m[2])[:]
    B = createDiffusiveOperatorCentered(h,m)
    S, dS, d2S = regularizer(u, B)
    @test S[1] ≈ 0
    @test sum(dS) == 0.0
end

@testset "check if constant deformation yields zero" begin
    omega = [0.0,1.0,0.0,1.0];
    m = [5; 7];
    h = [1.0,1.0];
    u = ones(m[1],2*m[2])
    #u[3,3] = 1
    u = u[:]
    B = createCurvatureOperatorCentered(h,m)
    S, dS, d2S = regularizer(u, B)
    @test S[1] ≈ 0
    @test sum(dS) == 0.0

    #using PyPlot
    #imshow(reshape(dS[1:35],(5,7)))
    #plot(dS[1:25])
end

@testset "check if random deformation yields nonzero result" begin
    omega = [0.0,1.0,0.0,1.0];
    m = [5; 7];
    h = [1.0,1.0];
    u = rand(2*m[1],m[2])[:]
    B = createCurvatureOperatorCentered(h,m)
    S, dS, d2S = regularizer(u, B)
    @test S[1] > 0
    @test (dS'*dS)[1] > 0
end

@testset "check if random deformation yields nonzero result" begin
    omega = [0.0,1.0,0.0,1.0];
    m = [5; 7];
    h = [1.0,1.0];
    u = rand(2*m[1],m[2])[:]
    B = createDiffusiveOperatorCentered(h,m)
    S, dS, d2S = regularizer(u, B)
    @test S[1] > 0
    @test (dS'*dS)[1] > 0
end
end #testset regularizer
# speed test
#omega = [0.0,1.0,0.0,1.0];
#m = [512; 512];
#h = [1.0,1.0];
#u = rand(size(getStaggeredGrid(omega,m)))
#@time B = createElasticOperatorStaggered(h,m)
#@time S, dS, d2S = regularizer(u, B)

# visualization
#using PyPlot
#spy(B)


text:
# This is free and unencumbered software released into the public domain.

# Anyone is free to copy, modify, publish, use, compile, sell, or
# distribute this software, either in source code form or as a compiled
# binary, for any purpose, commercial or non-commercial, and by any
# means.

# In jurisdictions that recognize copyright laws, the author or authors
# of this software dedicate any and all copyright interest in the
# software to the public domain. We make this dedication for the benefit
# of the public at large and to the detriment of our heirs and
# successors. We intend this dedication to be an overt act of
# relinquishment in perpetuity of all present and future rights to this
# software under copyright law.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.

# For more information, please refer to <https://unlicense.org>

if Sys.islinux()
    println("Hello World from Julia!")
end


text:
vals = readdlm(ARGS[1], Int)

numV = vals[1]
numEdge = vals[2]
numEl = vals[3]
numBoundaryEdges = vals[4]
numInterfaces = numEdge - numBoundaryEdges
numVertNodes = vals[5]
numEdgeNodes = vals[6]
numFaceNodes = vals[7]
numFaceNodes = 2*numVertNodes + numEdgeNodes
numNodes = 3*numVertNodes + 3*numEdgeNodes + numFaceNodes # per element
numDofPerNode  = vals[8]
numNodes_total = numVertNodes*numV + numEdgeNodes*numEdge + numFaceNodes*numEl

numDof = numNodes_total*numDofPerNode
jac_mem = vals[9]

println("numInterfaces = ", numInterfaces)
println("numFaceNodes = ", numFaceNodes)
println("numNodes per element = ", numNodes)
println("total number of nodes = ", numNodes_total)
println("numDof = ", numDof)

# add up the memory usage of all large arrays

verts = numV*64
edges = numEdge*64
elements = numEl*64
bndryfaces = 40*numBoundaryEdges
interfaces = 96*numInterfaces
interface_normals = 2*2*numFaceNodes*numInterfaces*64
coords = 2*numNodes*numEl*64
dxidx = 2*2*numNodes*numEl*64
jac = numNodes*numEl*64
dofs = numDofPerNode*numNodes*numEl*32
sparsity_bnds = 2*numDof*32
color_masks = 10*numEl*1
neighbor_colors = 4*numEl*8
neighbor_nums = 4*numEl*32

mesh_sum = verts + edges + elements + bndryfaces + interfaces + interface_normals + coords + dxidx + jac + dofs + sparsity_bnds + color_masks + neighbor_colors + neighbor_nums

println("mesh memory usage estimate: ", mesh_sum, " bits")
println("                          = ", mesh_sum/8, " bytes")
println("                          = ", mesh_sum/(8*1024*1024), " megabytes")

q = numDofPerNode*numNodes*numEl*64
aux_vars = 1*numEl*64
flux_parametric = numDofPerNode*numNodes*numEl*2*64
res = numDofPerNode*numNodes*numEl*64
res_vec = numDof*64
q_vec = numDof*64
edgestab_alpha = 2*2*numNodes*numEl*64
bndryflux = numDofPerNode*numFaceNodes*numBoundaryEdges*64
stabscale = numNodes*numInterfaces*64
Minv = numDof*64

eqn_sum = q + aux_vars + flux_parametric + res + res_vec + q_vec + edgestab_alpha + bndryflux + stabscale + Minv

println("solution data memory usage estimate: ", eqn_sum, " bits")
println("                                   = ", eqn_sum/8, " bytes")
println("                                   = ", eqn_sum/(8*1024*1024), " megabytes")

total_sum = eqn_sum + mesh_sum

println("total memory usage estimate: ", total_sum, " bits")
println("                           = ", total_sum/8, " bytes")
println("                           = ", total_sum/(8*1024*1024), " megabytes")

println("Jacobian memory usage estimate: ", jac_mem, " bits")
println("                           = ", jac_mem/8, " bytes")
println("                           = ", jac_mem/(8*1024*1024), " megabytes")


text:
module Pedestrians

using LinearAlgebra
using Printf
using ProgressMeter
using Random
using Requires

export run!

include("utilities.jl")
include("core.jl")
include("room_elements.jl")
include("strategies_time.jl")
include("strategies_schedule.jl")
include("strategies_target.jl")
include("strategies_move.jl")
include("predefined.jl")

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots.jl")
end

end # module


text:
"""
     isschur(A::AbstractMatrix) -> Bool

Test whether `A` is a square matrix in a real or complex Schur form.
In the real case, it is only tested whether A is a quasi upper triangular matrix,
which may have 2x2 diagonal blocks, which however must not correspond to
complex conjugate eigenvalues. In the complex case, it is tested if `A`
is upper triangular.
"""
function isschur(A)
   require_one_based_indexing(A)
   m, n = size(A)
   m == n || (return false)
   m == 1 && (return true)
   eltype(A)<:Complex && (return istriu(A))
   m == 2 && (return true)
   if istriu(A,-1)
      for i = 1:m-2
         if !iszero(A[i+1,i]) & !iszero(A[i+2,i+1])
            return false
         end
      end
      return true
   else
      return false
   end
end
"""
     isschur(A::AbstractMatrix, B::AbstractMatrix) -> Bool

Test whether `(A,B)` is a pair of square matrices in a real or complex
generalized Schur form.
In the real case, it is only tested whether `B` is upper triangular and `A` is
a quasi upper triangular matrix, which may have 2x2 diagonal blocks, which
however must not correspond to complex conjugate eigenvalues.
In the complex case, it is tested if `A` and `B` are both upper triangular.

"""
function isschur(A::AbstractMatrix,B::AbstractMatrix)
   ma, na = size(A)
   mb, nb = size(B)
   if eltype(A) != eltype(B) || ma != na || mb != nb || na != nb
      return false
   end
   return isschur(A) && istriu(B)
end
function sfstruct(A)
   # determine the structure of the generalized real Schur form
   n = size(A,1)
   ba = fill(1,n)
   p = 0
   i = 1
   while i < n
      p += 1
      iszero(A[i+1,i]) ? (i += 1) : (ba[p] = 2; i += 2)
   end
   i == n && (p += 1)
   return resize!(ba,p), p
   #return ba[1:p], p
end

"""
    utqu!(Q,U) -> Q

Compute efficiently the symmetric/hermitian product `U'QU -> Q`,
where `Q` is a symmetric/hermitian matrix and `U` is a square matrix.
The resulting product overwrites `Q`.
"""
function utqu!(Q,U)
   n = LinearAlgebra.checksquare(Q)
   ishermitian(Q) || error("Q must be a symmetric/hermitian matrix")
   LinearAlgebra.checksquare(U) == n ||
      throw(DimensionMismatch("U must be a matrix of dimension $n x $n"))

   T = promote_type(eltype(Q),eltype(U))
   if !(T <: BlasFloat)
      if T == eltype(Q)
         Q[:,:] = U'*Q*U
         return
      else
         error("TypeError: same type expected for Q and U ")
      end
   end

   Qd = view(Q,diagind(Q))
   rmul!(Qd,one(T)/2)
   if isa(U,Adjoint)
      mul!(Q, U.parent*UpperTriangular(Q), U)
   else
      mul!(Q, Adjoint(U), UpperTriangular(Q)*U)
   end
   #Q = tmp+tmp'
   @inbounds  begin
      for j = 1:n
       Q[j,j] += conj(Q[j,j])
       for i = j+1:n
           Q[i,j] += conj(Q[j,i])
           Q[j,i] = conj(Q[i,j])
       end
   end
   end
   return Q
end
"""
    X = utqu(Q,U)

Compute efficiently the symmetric/hermitian product `X = U'QU`,
where Q is a symmetric/hermitian matrix.
"""
function utqu(Q,U)
   n = LinearAlgebra.checksquare(Q)
   ishermitian(Q) || error("Q must be a symmetric/hermitian matrix")
   adj = isa(U,Adjoint)
   if adj
      m, n1 = size(U.parent)
      n1 == n || throw(DimensionMismatch("U must be a matrix of column dimension $n"))
   else
      n1, m = size(U)
      n1 == n || throw(DimensionMismatch("U must be a matrix of row dimension $n"))
   end

   T = promote_type(eltype(Q),eltype(U))
   if !(T <: BlasFloat)
      return U'*Q*U
   end

   t = UpperTriangular(Q)-Diagonal(Q[diagind(Q)]./2)
   X = similar(Q, T, m, m)
   if adj
      mul!(X, U.parent*t, U)
   else
      mul!(X, Adjoint(U), t*U)
   end
   @inbounds  begin
   for j = 1:m
       X[j,j] += conj(X[j,j])
       for i = j+1:m
           X[i,j] += conj(X[j,i])
           X[j,i] = conj(X[i,j])
      end
   end
   end
   return X
end
"""
    qrupdate!(R, Y) -> R

Update the upper triangular factor `R` by the
upper triangular factor of the QR factorization  of `[ R; Y' ]`, where `Y` is a
low-rank matrix `Y` (typically with one or two columns). The computation of `R`
only uses `O(n^2)` operations (`n` is the size of `R`). The input matrix `R` is
updated in place and the matrix `Y` is destroyed during the computation.
"""
function qrupdate!(R, Y)
    n, m = size(Y)
    size(R,1) == n || throw(DimensionMismatch("updating matrix must fit size of upper triangular matrix"))
    #Y = conj(Y)
    for k = 1:m
        for i = 1:n

            # Compute Givens rotation
            #c, s, r = LinearAlgebra.givensAlgorithm(R[i,i], conj(Y[i,k]))
            c, s, r = LinearAlgebra.givensAlgorithm(R[i,i], Y[i,k])

            # Store new diagonal element
            R[i,i] = r

            # Update remaining elements in row/column
            for j = i + 1:n
                Rij = R[i,j]
                yjk  = Y[j,k]
                R[i,j]  =   c*Rij + s*yjk
                Y[j,k]    = -s'*Rij + c*yjk
            end
        end
    end
    return R
end
"""

    rqupdate!(R, Y) -> R

Update the upper triangular factor `R` by the
upper triangular factor of the RQ factorization  of `[ Y R]`, where `Y` is a
low-rank matrix `Y` (typically with one or two columns). The computation of `R`
only uses `O(n^2)` operations (`n` is the size of `R`). The input matrix `R` is
updated in place and the matrix `Y` is destroyed during the computation.
"""
function rqupdate!(R, Y)
    n, m = size(Y)
    size(R,1) == n || throw(DimensionMismatch("updating matrix must fit size of upper triangular matrix"))

    for k = 1:m
        for j = n:-1:1

            # Compute Givens rotation
            c, s, r = LinearAlgebra.givensAlgorithm(R[j,j], Y[j,k])

            # Store new diagonal element
            R[j,j] = r

            # Update remaining elements in row/column
            for i = 1: j - 1
                Rij = R[i,j]
                yik  = Y[i,k]
                R[i,j]  =   c*Rij + s*yik
                Y[i,k]    = -s'*Rij + c*yik
            end
        end
    end
    return R
end
"""
    x = triu2vec(Q; rowwise = false, her = false)

Reshape the upper triangular part of the `nxn` array `Q` as a one-dimensional column
vector `x` with `n(n+1)/2` elements. `Q` is assumed symmetric/hermitian if `her = true`.
The elements of `x` correspond to stacking the elements of successive columns
of the upper triangular part of `Q`, if `rowwise = false`, or stacking the elements
of successive rows of the upper triangular part of `Q`, if `rowwise = true`.
"""
function triu2vec(Q::AbstractArray{T}; rowwise::Bool = false, her::Bool = false) where T
   n = LinearAlgebra.checksquare(Q)
   her && !ishermitian(Q) && error("Q must be a symmetric/hermitian matrix")
   x = Array{T,1}(undef,Int(n*(n+1)/2))
   k = 1
   if rowwise
      for i = 1:n
         for j = i:n
             x[k] = Q[i,j]
             k += 1
         end
      end
   else
      for j = 1:n
         for i = 1:j
             x[k] = Q[i,j]
             k += 1
         end
      end
   end
   return x
end
"""
    Q = vec2triu(x; rowwise = false, her = false)

Build from a one-dimensional column vector `x` with `n(n+1)/2` elements
an `nxn` upper triangular matrix `Q` if `her = false` or an `nxn` symetric/hermitian
array `Q` if `her = true`.
The elements of `x` correspond to stacking the elements of successive columns
of the upper triangular part of `Q`, if `rowwise = false`, or stacking the elements
of successive rows of the upper triangular part of `Q`, if `rowwise = true`.
"""
function vec2triu(x::AbstractVector{T}; rowwise = false, her = false) where T
   k = length(x)
   n = (sqrt(1+8*k)-1)/2
   isinteger(n) ? n = Int(n) : error("The lenght of x must be of the form n(n+1)/2")
   her ? Q = Array{T,2}(undef,n,n) : Q = zeros(T,n,n)
   k = 1
   if rowwise
      for i = 1:n
         Q[i,i] = x[k]
         k += 1
         for j = i+1:n
            Q[i,j] = x[k]
            k += 1
         end
      end
   else
      for j = 1:n
         for i = 1:j-1
            Q[i,j] = x[k]
            k += 1
         end
         Q[j,j] = x[k]
         k += 1
      end
   end
   if her
      for j = 1:n
         for i = j+1:n
            Q[i,j] = conj(Q[j,i])
         end
      end
   end
   return Q
end
function utnormalize!(U::UpperTriangular{T},adj::Bool) where T
   # Normalize an upper traiangular matrix U such that its diagonal elements are non-negative
   # using diagonal orthogonal or unitary transformations.
   ZERO = zero(real(T))
   n = size(U,1)
   if adj
      # Make the diagonal elements of U non-negative.
      if T <: Real
            for i = 1:n
               U[i,i] > ZERO || (rmul!(view(U,i,i:n),-1))
            end
      else
         for i = 1:n
            d = abs(U[i,i])
            (!iszero(d) && !(iszero(imag(U[i,i])) && real(U[i,i]) > ZERO)) &&
               (tmp = conj(U[i,i])/d; rmul!(view(U,i,i:n),tmp))
               #(tmp = conj(U[i,i])/d; [@inbounds U[i,j] *= tmp for j = i:n])
         end
      end
   else
      # Make the diagonal elements of U non-negative.
      if T <: Real
         for j = 1:n
            U[j,j] > ZERO || (rmul!(view(U,1:j,j),-1))
         end
      else
         for j = 1:n
            d = abs(U[j,j])
            (!iszero(d) && (iszero(imag(U[j,j])) && real(U[j,j]) > ZERO)) &&
               (tmp = conj(U[j,j]) / d; rmul!(view(U, 1:j, j), tmp))
               #(tmp = conj(U[j,j])/d; [@inbounds U[i,j] *= tmp for i = 1:j])
         end
      end
   end
   return U
end
@inline function luslv!(A::AbstractMatrix{T}, B::AbstractVector{T}) where T
   #
   #  fail = luslv!(A,B)
   #
   # This function is a speed-oriented implementation of a Gaussion-elimination based
   # solver of small order linear equations of the form A*X = B. The computed solution X
   # overwrites the vector B, while the resulting A contains in its upper triangular part,
   # the upper triangular factor U of its LU decomposition.
   # The diagnostic output parameter fail, of type Bool, is set to false in the case
   # of normal return or is set to true if the exact singularity of A is detected
   # or if the resulting B has non-finite components.
   #
   n = length(B)
   @inbounds begin
         for k = 1:n
            # find index max
            kp = k
            if k < n
                amax = abs(A[k, k])
                for i = k+1:n
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            iszero(A[kp,k]) && return true
            if k != kp
               # Interchange
               for i = 1:n
                   tmp = A[k,i]
                   A[k,i] = A[kp,i]
                   A[kp,i] = tmp
               end
               tmp = B[k]
               B[k] = B[kp]
               B[kp] = tmp
            end
            # Scale first column
            Akkinv = inv(A[k,k])
            i1 = k+1:n
            Ak = view(A,i1,k)
            rmul!(Ak,Akkinv)
            # Update the rest of A and B
            for j = k+1:n
                axpy!(-A[k,j],Ak,view(A,i1,j))
            end
            axpy!(-B[k],Ak,view(B,i1))
         end
         ldiv!(UpperTriangular(A), B)
         return any(!isfinite, B)
   end
end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
module YetAnotherHelperBubbleField

using ..Ahorn, Maple

@mapdef Entity "YetAnotherHelper/BubbleField" BubbleField(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
	strength::Number=1.0, direction::String="Right", flag::String="bubble_push_field", activationMode::String="Always", liftOffOfGround::Bool=true)

const placements = Ahorn.PlacementDict(
	"Bubble Column (Yet Another Helper)" => Ahorn.EntityPlacement(
		BubbleField,
		"rectangle"
	)
)

Ahorn.editingOptions(entity::BubbleField) = Dict{String,Any}(
	"direction" => ["Up", "Down", "Left", "Right"],
	"activationMode" => ["Always", "OnlyWhenFlagActive", "OnlyWhenFlagInactive"]
)

Ahorn.minimumSize(entity::BubbleField) = 8, 8
Ahorn.resizable(entity::BubbleField) = true, true

Ahorn.selection(entity::BubbleField) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BubbleField, room::Maple.Room)
	x = Int(get(entity.data, "x", 0))
	y = Int(get(entity.data, "y", 0))

	width = Int(get(entity.data, "width", 32))
	height = Int(get(entity.data, "height", 32))

	Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.0, 0.0, 1.0, 0.4), (1.0, 1.0, 1.0, 0.5))
end

end

max_stars_repo_name: YetAnotherHelper/YetAnotherHelper

max_issues_repo_name: YetAnotherHelper/YetAnotherHelper

max_forks_repo_name: YetAnotherHelper/YetAnotherHelper


text:
using ClustForOpt

ts_input_data = load_timeseries_data("CEP", "GER_1";T=24) #CEP

attribute_weights=Dict("solar"=>1.0, "wind"=>2.0, "el_demand"=>3.0)

clust_res=run_clust(ts_input_data;n_init=10,n_clust=4,attribute_weights=attribute_weights) # default k-means


text:
isinstalled(pkg::AbstractString) = pkg != "METADATA" && pkg != "REQUIRE" && pkg[1] != '.' && Pkg.cd(isdir, pkg)

print("\nHDF5\n")
if isinstalled("HDF5")
    print("HDF5 already installed\n")
else
    print("\nAdding HDF5 package to julia\n")
    Pkg.add("HDF5")
    Pkg.update("HDF5")
end

print("\ProgressMeter\n")
if isinstalled("ProgressMeter")
    print("ProgressMeter already installed\n")
else
    print("\nAdding ProgressMeter package to julia\n")
    Pkg.add("ProgressMeter")
    Pkg.update("ProgressMeter")
end

max_stars_repo_name: jiatianxu/pyLLE

max_issues_repo_name: jiatianxu/pyLLE

max_forks_repo_name: jiatianxu/pyLLE