text:
# Computes samples of real or complex Haar matrices of size nxn
#
# For beta=1,2,4, generates random orthogonal, unitary and symplectic matrices
# of uniform Haar measure.
# These matrices are distributed with uniform Haar measure over the
# classical orthogonal, unitary and symplectic groups O(N), U(N) and
# Sp(N)~USp(2N) respectively.
#
# The last parameter specifies whether or not the piecewise correction
# is applied to ensure that it truly of Haar measure
# This addresses an inconsistency in the Householder reflections as
# implemented in most versions of LAPACK
# Method 0: No correction
# Method 1: Multiply rows by uniform random phases
# Method 2: Multiply rows by phases of diag(R)
# References:
#    Edelman and Rao, 2005
#    Mezzadri, 2006, math-ph/0609050
#TODO implement O(n^2) method
#By default, always do piecewise correction
#For most applications where you use the HaarMatrix as a similarity transform
#it doesn't matter, but better safe than sorry... let the user choose else
function rand(W::Haar, n::Int, doCorrection::Int=1)
    beta = W.beta
    M=rand(Ginibre(beta,n))
    q,r=qr(M)
    if doCorrection==0
        q
    elseif doCorrection==1
        if beta==1
            L = sign.(rand(n).-0.5)
        elseif beta==2
            L = exp.(im*rand(n)*2pi)
        elseif beta==4
            L = exp.(im*rand(2n)*2pi)
        else
            error(string("beta = ",beta, " not implemented."))
        end
        q*diagm(0 => L)
    elseif doCorrection==2
        if beta==1
            L=sign.(diag(r))
        elseif (beta==2 || beta==4)
            L=diag(r)
            L=L./abs.(L)
        else
            error(string("beta = ",beta, " not implemented."))
        end
        q*diagm(0 => L)
    end
end

#A utility method to check if the piecewise correction is needed
#This checks the R part of the QR factorization; if correctly done,
#the diagonals are all chi variables so are non-negative
function NeedPiecewiseCorrection()
    n=20
    R=qr(randn(Ginibre(2,n)))[2]
    return any([x<0 for x in diag(R)])
end

#TODO maybe, someday
#Haar measure on U(N) in terms of local coordinates
#Zyczkowski and Kus, Random unitary matrices, J. Phys. A: Math. Gen. 27,
#4235–4245 (1994).

import LinearAlgebra.BLAS: @blasfunc

using LinearAlgebra: BlasInt
for (s, elty) in (("dlarfg_", Float64),
                  ("zlarfg_", ComplexF64))
    @eval begin
        function larfg!(n::Int, α::Ptr{$elty}, x::Ptr{$elty}, incx::Int, τ::Ptr{$elty})
	    ccall((@blasfunc($s), LAPACK.liblapack), Nothing,
		  (Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}),
                  Ref(n), α, x, Ref(incx), τ)
        end
    end
end

"""
Stewarts algorithm for n^2 orthogonal random matrix
"""
function Stewart(::Type{Float64}, n)
    τ = Array{Float64}(undef, n)
    H = randn(n, n)

    pτ = pointer(τ)
    pβ = pointer(H)
    pH = pointer(H, 2)

    for i = 0:n-2
        larfg!(n - i, pβ + (n + 1)*8i, pH + (n + 1)*8i, 1, pτ + 8i)
    end
    LinearAlgebra.QRPackedQ(H,τ)
end

function Stewart(::Type{ComplexF64}, n)
    τ = Array{ComplexF64}(undef, n)
    H = complex.(randn(n, n), randn(n,n))

    pτ = pointer(τ)
    pβ = pointer(H)
    pH = pointer(H, 2)

    for i = 0:n-2
        larfg!(n - i, pβ + (n + 1)*16i, pH + (n + 1)*16i, 1, pτ + 16i)
    end
    LinearAlgebra.QRPackedQ(H,τ)
end

export randfast
function randfast(W::Haar, n::Int)
    if W.beta==1
        Stewart(Float64, n)
    elseif W.beta==2
        Stewart(ComplexF64, n)
    else
        error("beta = $beta not implemented")
    end
end


text:
# library
# Since @rimport and @rlibrary create module objects which may be conflict with other objects,
# it is safer to place them at the end of the tests.

module NamespaceTests
    using RCall
    using Test

    RCall.R"""
        has_mass_package = require("MASS")
    """

    RCall.@rget has_mass_package

    @info "" has_mass_package

    if has_mass_package # these tests will fail if the R "MASS" package is not already installed
        @rimport MASS
        @test rcopy(rcall(MASS.ginv, RObject([1 2; 0 4]))) ≈ [1 -0.5; 0 0.25]
        @rimport MASS as mass
        @test rcopy(rcall(mass.ginv, RObject([1 2; 0 4]))) ≈ [1 -0.5; 0 0.25]
        @rlibrary MASS
        @test rcopy(rcall(ginv, RObject([1 2; 0 4]))) ≈ [1 -0.5; 0 0.25]
    end

end

ext: jl


text:
arrayrange(a; dims=dims) = maximum(a, dims=dims) - minimum(a, dims=dims)
removediag(m::Matrix) = (m[diagind(m)].=1e-05)
getintegers(m::Matrix) = m[findall(isinteger.(m))]

function puzzle02(path::String=joinpath(@__DIR__, "..", "data/02.txt"))
    x = readdlm(path, Int)
    partone = arrayrange(x, dims=2) |> sum

    y = [x[i,:] ./ x[i,:]' for i in 1:size(x, 2)]
    removediag.(y)
    parttwo = Int(sum(getintegers.(y))[1])

    return [partone, parttwo]
end


text:
using Rhapsodie

x = PolarimetricPixel("intensities", 100., 1., pi/8);
n1,n2 = (500, 499);
X = PolarimetricMap("intensities", n1,n2)

S=abs.(randn(5,5,3));
X=Rhapsodie.PolarimetricMap("stokes", S);
x=Rhapsodie.PolarimetricPixel("stokes", S[1,1,:]);
X[1,1]
X[1,1].I
Rhapsodie.write(X,"test_results/rand_test.fits")
Rhapsodie.read("mixed","test_results/rand_test.fits")

max_stars_repo_name: LaurenceDenneulin/Rhapsodie


text:
using StatisticalRethinking, Literate

"""

# generate

Utility function to generate all notebooks and chapters from scripts in the scripts directory.

### Method
```julia
generate(sd = script_dict)
```

### Required arguments

None, all notebooks/.. and chapters/.. files are regenerated.

"""
function generate(sd::DataStructures.OrderedDict{AbstractString, Vector{ScriptEntry}} = script_dict)
  DocDir = rel_path("..", "docs", "src")

  for chapter in keys(sd)
    ProjDir = rel_path("..", "scripts", chapter)

    ChapterDir = rel_path("..", "chapters")
    !isdir(ChapterDir) && mkdir(ChapterDir)
    ChapterDir = rel_path("..", "chapters", "$(chapter)")
    ScriptsDir = rel_path("..", "scripts", "$(chapter)")

    NotebookDir = rel_path("..", "notebooks")
    !isdir(NotebookDir) && mkdir(NotebookDir)
    NotebookDir = rel_path("..", "notebooks", "$(chapter)")

    !isdir(ProjDir) && break

    cd(ProjDir) do
      for script in sd[chapter]
        file = script.scriptfile

        # Generate chapters

        isfile(joinpath(ChapterDir, file[1:end-3], ".jl")) &&
          rm(joinpath(ChapterDir, file[1:end-3], ".jl"))
        Literate.script(file, ChapterDir)

        # Generate notebooks

        if script.nb && isfile(file)
          isfile(joinpath(NotebookDir, file[1:end-3], ".ipynb")) &&
            rm(joinpath(NotebookDir, file[1:end-3], ".ipynb"))
          Literate.notebook(file, NotebookDir, execute=script.exe)
        end
      end

      # Remove tmp directory used by cmdstan
      isdir(joinpath(ScriptsDir, "tmp")) &&
        rm(joinpath(ScriptsDir, "tmp"), recursive=true);

      println("\nCompleted script and notebook generation for chapter $chapter\n")

    end
  end
end

"""

# generate

Generate notebooks and scripts in a single chapter.

### Method
```julia
generate(chapter::AbstractString)
```

### Required arguments

Generate notebooks and scripts in a single chapter, e.g. generate("04")

"""
function generate(chapter::AbstractString; sd=script_dict)
  split_chapter = split(chapter, "/")
  if length(split_chapter) == 2
    generate(split_chapter...)
  else
    DocDir = rel_path("..", "docs", "src")
    ProjDir = rel_path("..", "scripts", chapter)

    ChapterDir = rel_path("..", "chapters")
    !isdir(ChapterDir) && mkdir(ChapterDir)
    ChapterDir = rel_path("..", "chapters", "$(chapter)")
    ScriptsDir = rel_path("..", "scripts", "$(chapter)")

    NotebookDir = rel_path("..", "notebooks")
    !isdir(NotebookDir) && mkdir(NotebookDir)
    NotebookDir = rel_path("..", "notebooks", "$(chapter)")

    if isdir(ProjDir)

      cd(ProjDir) do
        for script in sd[chapter]
          file = script.scriptfile

          # Generate chapters

          isfile(joinpath(ChapterDir, file[1:end-3], ".jl")) &&
            rm(joinpath(ChapterDir, file[1:end-3], ".jl"))
          Literate.script(file, ChapterDir)

          # Generate notebooks

          if script.nb && isfile(file)
            isfile(joinpath(NotebookDir, file[1:end-3], ".ipynb")) &&
              rm(joinpath(NotebookDir, file[1:end-3], ".ipynb"))
            Literate.notebook(file, NotebookDir, execute=script.exe)
          end
        end

        # Remove tmp directory used by cmdstan
        isdir(joinpath(ScriptsDir, "tmp")) &&
          rm(joinpath(ScriptsDir, "tmp"), recursive=true);

        println("\nCompleted script and notebook generation for chapter $chapter\n")

      end
    end
  end
end

"""

# generate

Generate a single notebook and script

### Method
```julia
generate(chapter::AbstractString, file::AbstractString)
```

### Required arguments

Generate notebook and script `file` in `chapter`, e.g. generate("04", "m4.1d.jl")
or  generate("04/m4.1d.jl")

"""
function generate(chapter::AbstractString, scriptfile::AbstractString; sd=script_dict)
  DocDir = rel_path("..", "docs", "src")
  ProjDir = rel_path("..", "scripts", chapter)

  ChapterDir = rel_path("..", "chapters")
  !isdir(ChapterDir) && mkdir(ChapterDir)
  ChapterDir = rel_path("..", "chapters", "$(chapter)")
  ScriptsDir = rel_path("..", "scripts", "$(chapter)")

  NotebookDir = rel_path("..", "notebooks")
  !isdir(NotebookDir) && mkdir(NotebookDir)
  NotebookDir = rel_path("..", "notebooks", "$(chapter)")

  if isdir(ProjDir)

    cd(ProjDir) do
      for script in sd[chapter]
        !(script.scriptfile == scriptfile) && continue
        file = script.scriptfile

        # Generate chapters

        isfile(joinpath(ChapterDir, file[1:end-3], ".jl")) &&
          rm(joinpath(ChapterDir, file[1:end-3], ".jl"))
        Literate.script(file, ChapterDir)

        # Generate notebooks

        if script.nb && isfile(file)
          isfile(joinpath(NotebookDir, file[1:end-3], ".ipynb")) &&
            rm(joinpath(NotebookDir, file[1:end-3], ".ipynb"))
          Literate.notebook(file, NotebookDir, execute=script.exe)
        end
      end

      # Remove tmp directory used by cmdstan
      isdir(joinpath(ScriptsDir, "tmp")) &&
        rm(joinpath(ScriptsDir, "tmp"), recursive=true);

      println("\nCompleted script and notebook generation for chapter $chapter\n")

    end
  end
end


text:
module TestGraphviz

using Test
import JSON

using Catlab.CategoricalAlgebra.Graphs: PropertyGraph, SymmetricPropertyGraph,
  nv, ne, src, tgt, add_vertex!, add_vertices!, add_edge!, get_vprop, get_eprop
using Catlab.Graphics.Graphviz

# Graphviz to property graph
############################

data_path(name::String) = joinpath(@__DIR__, "data", name)

# Undirected simple graph.
doc = open(JSON.parse, data_path("graphviz_graph.json"), "r")
parsed = parse_graphviz(doc)
@test parsed isa SymmetricPropertyGraph
@test nv(parsed) == 10
@test ne(parsed) ÷ 2 == 13

# Directed simple graph
doc = open(JSON.parse, data_path("graphviz_digraph.json"), "r")
parsed = parse_graphviz(doc)
@test parsed isa PropertyGraph
@test nv(parsed) == 5
@test src(parsed) == [1,1,1,2,3,4]
@test tgt(parsed) == [2,3,4,5,5,5]
@test get_vprop.([parsed], 1:nv(parsed), :name) == ["A", "B", "C", "D", "F"]
@test get_vprop(parsed, 1, :position) == [99, 162]
@test get_vprop(parsed, 1, :size) == [54, 36]
@test get_eprop(parsed, 1, :spline) isa Vector

# Property graph to Graphviz
############################

# Undirected simple graph
g = SymmetricPropertyGraph{String}()
add_vertex!(g, label="v")
add_vertices!(g, 2)
add_edge!(g, 1, 2, xlabel="e")
gv = to_graphviz(g)::Graph
@test !gv.directed
nodes = filter(s -> isa(s,Node), gv.stmts)
edges = filter(s -> isa(s,Edge), gv.stmts)
@test length(nodes) == 3
@test length(edges) == 1
@test nodes[1].attrs[:label] == "v"
@test edges[1].attrs[:xlabel] == "e"

# Directed simple graph
g = PropertyGraph{String}()
add_vertices!(g, 3); add_edge!(g, 1, 2); add_edge!(g, 2, 3)
gv = to_graphviz(g)::Graph
@test gv.directed
@test length(filter(s -> isa(s,Node), gv.stmts)) == 3
@test length(filter(s -> isa(s,Edge), gv.stmts)) == 2

# Directed multigraph
g = PropertyGraph{String}()
add_vertices!(g, 2)
add_edge!(g, 1, 2, label="e1")
add_edge!(g, 1, 2, label="e2")
gv = to_graphviz(g)::Graph
@test gv.directed
nodes = filter(s -> isa(s,Node), gv.stmts)
edges = filter(s -> isa(s,Edge), gv.stmts)
@test length(nodes) == 2
@test length(edges) == 2
@test edges[1].attrs[:label] == "e1"
@test edges[2].attrs[:label] == "e2"

# Pretty-print
##############

spprint(expr::Expression) = sprint(pprint, expr)

# Node statement
@test spprint(Node("n")) == "n;"
@test spprint(Node("n"; label="foo")) == "n [label=\"foo\"];"
@test spprint(Node("n"; shape="box", style="filled")) ==
  "n [shape=\"box\",style=\"filled\"];"

# Edge statement
@test spprint(Edge("n1","n2")) == "n1 -- n2;"
@test spprint(Edge("n1","n2"; label="foo")) == "n1 -- n2 [label=\"foo\"];"
@test spprint(Edge("n1","n2"; style="dotted", weight="10")) ==
  "n1 -- n2 [style=\"dotted\",weight=\"10\"];"
@test spprint(Edge(NodeID("n1","p1"), NodeID("n2","p2"))) == "n1:p1 -- n2:p2;"
@test spprint(Edge(NodeID("n1","p1"), NodeID("n2","p2"); label="bar")) ==
  "n1:p1 -- n2:p2 [label=\"bar\"];"
@test spprint(Edge(NodeID("n1","p1","w"), NodeID("n2", "p2", "e"))) ==
  "n1:p1:w -- n2:p2:e;"
@test spprint(Edge("n1","n2","n3")) == "n1 -- n2 -- n3;"
@test spprint(Edge(NodeID("n1"), NodeID("n2"), NodeID("n3"))) ==
  "n1 -- n2 -- n3;"

# Graph statement
graph = Graph("G",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(graph) == """
graph G {
  n1;
  n2;
  n1 -- n2;
}
"""

graph = Digraph("G",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(graph) == """
digraph G {
  n1;
  n2;
  n1 -> n2;
}
"""

graph = Digraph("G",
  Node("n1"; label="foo"),
  Node("n2"; label="bar"),
  Edge("n1","n2");
  graph_attrs = Attributes(:rankdir => "LR"),
  node_attrs = Attributes(:shape => "box", :style => "filled"),
  edge_attrs = Attributes(:style => "dotted")
)
@test spprint(graph) == """
digraph G {
  graph [rankdir="LR"];
  node [shape="box",style="filled"];
  edge [style="dotted"];
  n1 [label="foo"];
  n2 [label="bar"];
  n1 -> n2;
}
"""

# Subgraph statement
subgraph = Subgraph("sub",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(subgraph) == """
subgraph sub {
  n1;
  n2;
  n1 -- n2;
}"""

subgraph = Subgraph(
  Node("n1"),
  Node("n2"),
  graph_attrs = Attributes(:rank => "same")
)
@test spprint(subgraph) == """
{
  graph [rank="same"];
  n1;
  n2;
}"""

end


text:
using GradientDescent
using Test

@testset "GradientDescent.jl" begin
    # Write your tests here.
end


text:
module Private

export get_most_union_all, ast_and, ast_or
function get_most_union_all(expr, mod :: Module)
    if isa(expr, Expr) && expr.head == :curly
        get_most_union_all(expr.args[1], mod)
    else
        @eval mod $expr
    end
end

ast_and(a, b) =  Expr(:&&, a, b)
ast_or(a, b) = Expr(:||, a, b)

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

@testset "find intersection of Seg3 element" begin
    el = Element(Seg3, [1, 2, 3])
    update!(el, "geometry", Vector{Float64}[[0.0, 1.0], [1.0, 0.0], sqrt(2.0)/2.0*[1.0, 1.0]])
    s = [ 1.0, 0.5]
    d = [-1.0, 0.0]
    t, xi = find_intersection(el, s, d, 0.0)
    x = el("geometry", xi, 0.0)
    @test isapprox(x, [0.8604093371313943, 0.5])
end

@testset "find intersection of 2. order NSeg element" begin
    # this is a exact quarter of circle
    a = 1.0
    b = 1.0
    el = Element(NSeg, [1, 2, 3])
    el.properties.order = 2
    el.properties.knots = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
    el.properties.weights = [1.0, 1.0, 2.0]
    update!(el, "geometry", Vector{Float64}[[a, 0], [a, b], [0, b]])
    s = [ 1.0, 0.5]
    d = [-1.0, 0.0]
    t, xi = find_intersection(el, s, d, 0.0)
    X = el("geometry", xi, 0.0)
    @test isapprox(X, [sqrt(3)/2, 1/2])
end

@testset "find intersection of 1. order NSurf" begin
    # node ordering, it's not same as in Quad4
    el = Element(NSurf, [1, 2, 3, 4])
    el.properties.order_u = 1
    el.properties.order_v = 1
    el.properties.knots_u = [0.0, 0.0, 1.0, 1.0]
    el.properties.knots_v = [0.0, 0.0, 1.0, 1.0]
    el.properties.weights = ones(2, 2)
    # +-- v
    # |
    # u
    nodes = Vector{Float64}[[0.0,0.0,0.0], [1.0,0.0,0.0],
                            [0.0,1.0,0.0], [1.0,1.0,0.0]]
    update!(el, "geometry", nodes)
    s = [0.5, 0.5,  0.5]
    d = [0.0, 0.0, -1.0]
    t, xi = find_intersection(el, s, d, 0.0)
    X = el("geometry", xi, 0.0)
    @test isapprox(X, [0.5, 0.5, 0.0])
    k = calc_reflection(el, xi, d, 0.0)
    @test isapprox(k, [0.0, 0.0, 1.0])
    s = [1.5, 1.5, 0.5]
    t, xi = find_intersection(el, s, d, 0.0)
    @test isnan(t)
end

@testset "find intersection of 3. order NSeg element with multiple reflections" begin
    el = Element(NSeg, [1, 2, 3, 4])
    el.properties.order = 3
    el.properties.knots = [0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]
    el.properties.weights = [1.0, 1/3, 1/3, 1.0]
    update!(el, "geometry", Vector{Float64}[[1, 0], [1, 2], [-1, 2], [-1, 0]])
    s1 = [-sqrt(3)/2.0, 0.0]
    k1 = [ 0.0, 1.0]
    t1, xi1 = find_intersection(el, s1, k1, 0.0)
    s2 = el("geometry", xi1, 0.0)
    @test isapprox(s2, [-sqrt(3)/2, 0.5])
    k2 = calc_reflection(el, xi1, k1, 0.0)
    t2, xi2 = find_intersection(el, s2, k2, 0.0; secondary=true)
    s3 = el("geometry", xi2, 0.0)
    @test isapprox(s3, [0.0, 1.0])
    k3 = calc_reflection(el, xi2, k2, 0.0)
    t3, xi3 = find_intersection(el, s3, k3, 0.0; secondary=true)
    s4 = el("geometry", xi3, 0.0)
    @test isapprox(s4, [sqrt(3)/2, 0.5])
end


text:
@inline unshift(v::Vec, item) = Vec(item, Tuple(v)...)
@inline sqnorm(x) = sum(x.^2)

"""
pinvli{n,m,T}(a::Mat{n,m,T})

Compute pseudo inverse of matrix with linear independent columns.
"""
function pinvli{n,m,T}(a::Mat{n,m,T})
    @assert n >=m
    return inv(a'*a)*a'
end


text:
using CommonMark

template(body) = raw"""
\documentclass[12pt,a4paper]{article}

\usepackage[a4paper,text={16.5cm,25.2cm},centering]{geometry}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{graphicx}
% \usepackage{svg}
% \svgsetup{inkscapelatex=false}
\usepackage{microtype}
\usepackage{hyperref}
\setlength{\parindent}{0pt}
\setlength{\parskip}{1.2ex}

\usepackage{placeins}

\let\Oldsection\section
\renewcommand{\section}{\FloatBarrier\Oldsection}

\let\Oldsubsection\subsection
\renewcommand{\subsection}{\FloatBarrier\Oldsubsection}

\usepackage{float}
\floatplacement{figure}{h!}

\usepackage[defaultlines=3,all]{nowidow}

\title{PeaceVote}

\author{Janis Erdmanis}

\usepackage[T1]{fontenc}
\usepackage{textcomp}
\usepackage{upquote}
%\usepackage{listings}

\usepackage{jlcode}

\renewcommand{\texttt}[1]{\jlinl{#1}}

\let\verbatim\trash
\lstnewenvironment{verbatim}
{}{}

\let\oldincludegraphics\includegraphics
\renewcommand{\includegraphics}[1]{\oldincludegraphics[width=0.5\textwidth]{#1}}

\begin{document}

\maketitle

% \abstract{PeaceVote abstract}

""" *
body *
raw"""

\include{references}

\end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% LaTeX-command:"latex --synctex=1 --shell-escape"
%%% End:
"""

mtext = join(readlines("peacevote.md"),"\n")

markdown = Parser()
enable!(markdown, FootnoteRule())
ast = markdown(mtext)

body = latex(ast)
text = template(body)

open("peacevote.tex", "w") do file
    println(file, text)
end


text:
using CSDP
const CSDP_PRIMAL_FACTORY =
    MOI.OptimizerWithAttributes(CSDP.Optimizer, MOI.Silent() => true)
const CSDP_DUAL_FACTORY = dual_optimizer(CSDP_PRIMAL_FACTORY)
const CSDP_PRIMAL_OPT = MOI.instantiate(CSDP_PRIMAL_FACTORY)
const CSDP_DUAL_OPT = MOI.instantiate(CSDP_DUAL_FACTORY)

push!(primal_conic_factory, CSDP_PRIMAL_FACTORY)
push!(dual_conic_factory, CSDP_DUAL_FACTORY)
push!(dual_conic_optimizer, CSDP_DUAL_OPT)
push!(primal_conic_optimizer, CSDP_PRIMAL_OPT)

@testset "CSDP SDP triangle Problems" begin
    list_of_sdp_triang_problems = [
        # sdpt1_test, # CSDP is returning SLOW_PROGRESS
        # sdpt2_test, # CSDP is returning SLOW_PROGRESS
        sdpt3_test,
        sdpt4_test,
    ]
    test_strong_duality(list_of_sdp_triang_problems, CSDP_PRIMAL_FACTORY)
end


text:
using SafeTestsets
@safetestset "my f test" begin include("my_f_tests.jl") end


text:
using Clang.Generators
using FFMPEG
using FFMPEG.FFMPEG_jll
using Vulkan_Headers_jll
using LibGit2
include("rewriter.jl")

cd(@__DIR__)

# Ideally I could have loaded all headers, but it turns out to be too hard

include_dir = joinpath(FFMPEG.FFMPEG_jll.artifact_dir, "include") |> normpath

vulkan_dir = joinpath(Vulkan_Headers_jll.artifact_dir, "include") |> normpath

if !isdir("vdpau")
    LibGit2.clone("https://gitlab.freedesktop.org/vdpau/libvdpau", joinpath(@__DIR__, "vdpau"))
end

vdpau_dir = joinpath(@__DIR__, "vdpau", "include")

options = load_options(joinpath(@__DIR__, "generate.toml"))

args = get_default_args()
push!(args, "-I$include_dir", "-isystem$vulkan_dir", "-isystem$vdpau_dir")

const module_names = [
    "libavcodec"
    "libavdevice"
    "libavfilter"
    "libavformat"
    "libavutil"
    "libswscale"
]

library_names = Dict{String,String}()
headers = String[]
for lib in module_names
    header_dir = joinpath(include_dir, lib)
    append!(headers, joinpath(header_dir, header) for header in readdir(header_dir) if endswith(header, ".h"))
    library_names[lib*".+"] = lib
end

options["general"]["library_names"] = library_names

ctx = create_context(headers, args, options)

@add_def time_t

build!(ctx, BUILDSTAGE_NO_PRINTING)
for node in ctx.dag.nodes
    for i in eachindex(node.exprs)
        node.exprs[i] = rewrite(node.exprs[i])
    end
end
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
@testset "grids" begin
    bbox1 = [-96.6357421875,31.12819929911196,-84.9462890625,40.58058466412764]
    bbox2 = [-220.78125,-80.64703474739618,-29.53125,78.34941069014629]
    bbox3 = [-5.1,-5.1,5.1,5.1]

    feat2 = point_grid(bbox2, 500., nothing, "miles")
    @test feat2.features[1].geometry.coordinates ≈ [-199.504452, -77.132893]

    feat3 = rectangle_grid(bbox3, 1., 1., nothing, "degrees")
    @test feat3.features[1].geometry.coordinates ≈ [
					[
						[
							-5.025791,
							-5.005842
						],
						[
							-5.025791,
							-4.004674
						],
						[
							-4.020633,
							-4.004674
						],
						[
							-4.020633,
							-5.005842
						],
						[
							-5.025791,
							-5.005842
						]
					]] atol=0.0001

    feat1 = hexgrid(bbox1, 25., nothing, false, "miles")
"""
	@test feat1.features[1].geometry.coordinates ≈
    [
                [
                    [
                        -95.704593,
                        31.467449
                    ],
                    [
                        -95.927937,
                        31.780802
                    ],
                    [
                        -96.374626,
                        31.780802
                    ],
                    [
                        -96.59797,
                        31.467449
                    ],
                    [
                        -96.374626,
                        31.154096
                    ],
                    [
                        -95.927937,
                        31.154096
                    ],
                    [
                        -95.704593,
                        31.467449
                    ]
                ]]
				"""
	bbox4 = [-95.877685546875,74.46849062193377,-94.031982421875,74.90226611990785]
    fc = triangle_grid(bbox4, 12.5, nothing, "miles")
	@test fc.features[1].geometry.coordinates ≈ [
					[
						[
							-95.877686,
							74.468491
						],
						[
							-95.877686,
							74.649405
						],
						[
							-95.20202,
							74.468491
						],
						[
							-95.877686,
							74.468491
						]
					]]
	fc2 = square_grid(bbox4, 12.5, nothing, "miles")
	@test fc2.features[1].geometry.coordinates ≈ [
					[
						[
							-95.630499,
							74.504464
						],
						[
							-95.630499,
							74.685378
						],
						[
							-94.954834,
							74.685378
						],
						[
							-94.954834,
							74.504464
						],
						[
							-95.630499,
							74.504464
						]
					]]
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function testareaintersection()
    c1 = (O, 150)
    c2 = (O + (100, 0), 150)
    cia = intersection2circles(c1..., c2...)
    @test isapprox(cia, 41251.0, atol=0.1)
end

function test_circles(fname)
    pagewidth, pageheight = 1200, 1400
    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setopacity(0.8)
    setline(0.4)
    fontsize(2)
    gsave()
    translate(0, -pageheight/2)
    pagetiles = Tiler(pagewidth, pageheight, 20, 20, margin=50)
    rad = 40
    for (pos, n) in pagetiles
        diameter = rand(10:50)
        secondpos = Point(pos.x + diameter, pos.y)
        sethue("black")
        circle(pos, 1, :fill)
        circle(secondpos, 1, :fill)
        randomhue()
        circle(pos, secondpos, :fill)
    end
    grestore()

    gsave()
    translate(0, pageheight/2)
    pagetiles = Tiler(pagewidth, pageheight, 20, 20, margin=50)
    rad = 40
    for (p1, n) in pagetiles
        p2, p3 = Point(p1.x + rand(-15:15), p1.y + rand(-15:15)), Point(p1.x + rand(-15:15), p1.y + rand(-15:15))
        cpoint, rad = center3pts(p1, p2, p3) # returns (0/0, 0) if points are collinear
        randomcolor()
        if rad > 0
            circle(cpoint, rad, :stroke)
            sethue("black")
            circle(p1, 1, :fill)
            circle(p2,  1, :fill)
            circle(p3,  1, :fill)
            p4 = getnearestpointonline(p1, p2, cpoint)
            line(p1, p2, :stroke)
            line(cpoint, p4, :stroke)
            text(string(round(pointlinedistance(cpoint, p1, p2), digits=1)), cpoint)
        end
        setopacity(0.8)
    end
    grestore()

    @test finish() == true
    println("...finished circletest, saved in $(fname)")
end

testareaintersection()

fname = "circle-test.pdf"

test_circles(fname)


text:
export AbstractActionSelector, AbstractDiscreteActionSelector

"""
    AbstractActionSelector

Take in an estimation and return an action.
"""
abstract type AbstractActionSelector end

"""
    AbstractDiscreteActionSelector

Generate an action given the estimated value of different actions.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `selector(values; kwargs...)` | `selector`, an instance of `AbstractDiscreteActionSelector`, must be a callable object which takes in an estimation of all discrete actions and returns an action. |
"""
abstract type AbstractDiscreteActionSelector <: AbstractActionSelector end


text:
include("test_common.jl")
using Dates

@testset "promises" begin
    js_sleep = node"""(ms) => new Promise(
        resolve => setTimeout(() => resolve(ms), ms)
    )"""
    t = now()
    p = js_sleep(500)
    @test fetch(p) == 500
    @test now() - t ≥ Millisecond(500)
    @test state(napi_value(p)) == NodeCall.promise_fulfilled
    @test result(p) == 500
    async_f = node"""async (x) => {
        return x * 2
    }"""
    y = 0
    p2 = async_f(50).then() do x
        y = x
        x * 2
    end
    @test 200 == @await p2
    @test y == 100

    p3 = JsPromise() do resolve, _
        resolve(true)
    end
    @test @await p3
    p4 = JsPromise((_, reject) -> begin
        reject(123)
    end)
    @test_throws NodeError wait(p4)

    t = now()
    jl_task = @async sleep(3)
    node_sleep = @node_async sleep(1)
    # `node_sleep` should not block `jl_task`
    wait(node_sleep)
    # Neither should `js_sleep`
    wait(js_sleep(1000))
    wait(jl_task)
    @test now() - t ≤ Millisecond(3500)
end


text:
export RowGraph
#required methods
export getRowMap, getColMap, getDomainMap, getRangeMap, getImporter, getExporter
export getGlobalNumRows, getGlobalNumCols, getGlobalNumEntries, getGlobalNumDiags
export getLocalNumRows, getLocalNumCols, getLocalNumEntries, getLocalNumDiags
export getNumEntriesInGlobalRow, getNumEntriesInLocalRow
export getGlobalMaxNumRowEntries, getLocalMaxNumRowEntries
export hasColMap, isLowerTriangular, sUpperTriangular
export isLocallyIndexed, isGloballyIndexed, isFillComplete
export getGlobalRowCopy, getLocalRowCopy, pack
#implemented methods
export isFillActive

"""
RowGraph is the base "type" for all row oriented storage graphs

Instances of these types are required to implement the following submethods

    getRowMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}
Gets the row map for the graph

    getColMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}
Gets the column map for the graph

    getDomainMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}
Gets the domain map for the graph

    getRangeMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}
Gets the range map for the graph

    getImporter(::RowGraph{GID, PID, LID})::Union{Import{GID, PID, LID}, Nothing}
Gets the graph's Import object

    getExporter(::RowGraph{GID, PID, LID})::Union{Export{GID, PID, LID}, Nothing}
Gets the graph's Export object

    getGlobalNumEntries(::RowGraph{GID, PID, LID})::GID
Returns the global number of entries in the graph

    getLocalNumEntries(::RowGraph{GID, PID, LID})::LID
Returns the local number of entries in the graph

    getNumEntriesInLocalRow(::RowGraph{GID, PID, LID}, row::LID)::LID
Returns the current number of local entries in the given row

    getGlobalNumDiags(::RowGraph{GID, PID, LID})::GID
Returns the global number of diagonal entries

    getLocalNumDiags(::RowGraph{GID, PID, LID})::LID
Returns the local number of diagonal entries

    getGlobalMaxNumRowEntries(::RowGraph{GID, PID, LID})::LID
Returns the maximum number of entries across all rows/columns on all processors

    getLocalMaxNumRowEntries(::RowGraph{GID, PID, LID})::LID
Returns the maximum number of entries across all rows/columns on this processor

    hasColMap(::RowGraph{GID, PID, LID})::Bool
Whether the graph has a well-defined column map

    isLowerTriangular(::RowGraph{GID, PID, LID})::Bool
Whether the graph is lower trianguluar

    isUpperTriangular(::RowGraph{GID, PID, LID})::Bool
Whether the graph is upper trianguluar

    isGloballyIndexed(::RowGraph)::Bool
Whether the graph is using global indices

    isFillComplete(::RowGraph)
If the graph is fully build.

    getLocalRowCopy!(copy::AbstractVector{<:Integer}, graph::RowGraph{GID, PID, LID}, row::LID)::Integer
Copies the given row of the graph into the provided array using local indices

For optimal performance, the following methods may need to be re-implemented for the specific type
`getGlobalRowCopy!(copy::AbstractVector{<:Integer}, graph::RowGraph{GID}, row::GID)::Integer` is implemented by calling `getLocalRowCopy!` and remapping the values using `gid(::BlockMap, ::Integer)`

`getNumEntriesInGlobalRow(graph::RowGraph{GID, PID, LID}, row::GID)` is implemented using `getNumEntriesInLocalRow`

"""
abstract type RowGraph{GID <: Integer, PID <: Integer, LID <: Integer}
end

"""
    isFillActive(::RowGraph)

Whether the graph is being built
"""
isFillActive(graph::RowGraph) = !isFillComplete(graph)

"""
    getNumEntriesInLocalRow(::RowGraph{GID, PID, LID}, row::Integer)::LID

Returns the current number of local entries in the given row
"""
Base.@propagate_inbounds function getNumEntriesInLocalRow(graph::RowGraph{GID, PID, LID},
        row::Integer)::LID where{GID, PID, LID}
    getNumEntriesInLocalRow(graph, LID(row))
end

"""
    getGlobalRowCopy(graph::RowGraph{GID}, row::Integer)::AbstractVector{GID}

Creates a copy of the given row of the graph using global indices
"""
Base.@propagate_inbounds function getGlobalRowCopy(graph::RowGraph{GID}, row::Integer)::AbstractVector{GID} where GID
    numElts = getNumEntriesInGlobalRow(graph, row)
    copy = Vector{GID}(undef, numElts)
    getGlobalRowCopy!(copy, graph, GID(row))
    copy
end

"""
    getLocalRowCopy(graph::RowGraph{GID}, row::Integer)::AbstractVector{GID}

Creates a copy of the given row of the graph using local indices
"""
Base.@propagate_inbounds function getLocalRowCopy(graph::RowGraph{GID, PID, LID}, row::Integer)::AbstractVector{LID} where {GID, PID, LID}
    numElts = getNumEntriesInGlobalRow(graph, row)
    copy = Vector{LID}(undef, numElts)
    getLocalRowCopy!(copy, graph, LID(row))
    copy
end

"""
    getGlobalRowCopy!(copy::AbstractVector{GID}, graph::RowGraph{GID}, row::Integer)::Integer

Copies the given row of the graph into the provided array using global indices
"""
Base.@propagate_inbounds function getGlobalRowCopy!(copy::AbstractVector{<:Integer}, graph::RowGraph{GID},
        row::Integer)::Integer where GID
    getGlobalRowCopy!(graph, GID(row))
end

Base.@propagate_inbounds function getGlobalRowCopy!(copy::AbstractVector{<:Integer},
            matrix::RowGraph{GID}, globalRow::GID) where {GID}
    rowMap = getRowMap(matrix)
    numElts = getLocalRowCopy!(copy, matrix, lid(rowMap, globalRow))
    @. copy = gid(rowMap, copy)
    numElts
end

"""
    getLocalRowCopy!(copy::AbstractVector{<:Integer}, graph::RowGraph{GID, PID, LID}, row::Integer)::Integer

Copies the given row of the graph into the provided array using global indices
"""
Base.@propagate_inbounds function getLocalRowCopy!(copy::AbstractVector{<:Integer}, graph::RowGraph{GID, PID, LID},
        row::Integer)::Integer where {GID, PID, LID}
    getLocalRowCopy!(graph, LID(row))
end

"""
    isLocallyIndexed(::RowGraph)::Bool

Whether the graph is using local indices.

The default implementation uses the row graph
"""
isLocallyIndexed(graph::RowGraph) = !isGloballyIndexed(graph)

"""
    getGlobalNumRows(::RowGraph{GID})::GID

Returns the number of global rows in the graph

The default implementation uses the row graph
"""
getGlobalNumRows(graph::RowGraph) = numAllElements(getRowMap(graph))

"""
    getGlobalNumCols(::RowGraph{GID})::GID

Returns the number of global columns in the graph

The default implementation uses the row graph
"""
getGlobalNumCols(graph::RowGraph) = numAllElements(getColMap(graph))

"""
    getLocalNumRows(::RowGraph{GID, PID, LID})::LID

Returns the number of rows owned by the calling process

The default implementation uses the row graph
"""
getLocalNumRows(graph::RowGraph) = numMyElements(getRowMap(graph))

"""
    getLocalNumCols(::RowGraph{GID, PID, LID})::LID

Returns the number of columns owned by the calling process

The default implementation uses the row graph
"""
getLocalNumCols(graph::RowGraph) = nulMyElements(getColMap(graph))

"""
    getNumEntriesInGlobalRow(graph::RowGraph{GID, PID, LID}, row::Integer)::LID

Returns the current number of local entries in the given row
"""
function getNumEntriesInGlobalRow(graph::RowGraph{GID, PID, LID}, row::Integer
        )::LID where {GID, PID, LID}
    getNumEntriesInGlobalRow(graph, row::GID)
end

function getNumEntriesInGlobalRow(graph::RowGraph{GID, PID, LID}, row::GID
        )::LID where {GID, PID, LID}
    getNumEntriesInLocalRow(graph, lid(getRowMap(graph), row))
end

"""
    pack(::RowGraph{GID, PID, LID}, exportLIDs::AbstractArray{LID, 1}, distor::Distributor{GID, PID, LID})::AbstractArray{AbstractArray{LID, 1}}

Packs this object's data for import or export
"""
function pack(source::RowGraph{GID, PID, LID}, exportLIDs::AbstractArray{LID, 1}, distor::Distributor{GID, PID, LID})::Array{Array{LID, 1}, 1} where {GID, PID, LID}
    srcMap = getMap(source)
    map(lid->getGlobalRowCopy(source, gid(srcMap, lid)), exportLIDs)
end

#### SrcDistObject methods ####
getMap(graph::RowGraph) = getRowMap(graph)

#### documentation for required methods ####

"""
    isFillComplete(mat::RowGraph)

Whether `fillComplete(...)` has been called
"""
function isFillComplete end

"""
    getRowMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}

Gets the row map for the graph
"""
function getRowMap end

"""
    getColMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}

Gets the column map for the graph
"""
function getColMap end

"""
    hasColMap(::RowGraph{GID, PID, LID})::Bool

Whether the graph has a well-defined column map
"""
function hasColMap end

"""
    getDomainMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}

Gets the domain map for the graph
"""
function getDomainMap end

"""
    getRangeMap(::RowGraph{GID, PID, LID})::BlockMap{GID, PID, LID}

Gets the range map for the graph
"""
function getRangeMap end

"""
    getImporter(::RowGraph{GID, PID, LID})::Union{Import{GID, PID, LID}, Nothing}

Gets the graph's Import object, or null if the import is trivial
"""
function getImporter end

"""
    getExporter(::RowGraph{GID, PID, LID})::Union{Export{GID, PID, LID}, Nothing}

Gets the graph's Export object, or null if the export is trivial
"""
function getExporter end

"""
    getGlobalNumEntries(::RowGraph{GID, PID, LID})::GID

Returns the global number of entries in the graph
"""
function getGlobalNumEntries end

"""
    getLocalNumEntries(::RowGraph{GID, PID, LID})::LID

Returns the local number of entries in the graph
"""
function getLocalNumEntries end

"""
    getNumEntriesInLocalRow(::RowGraph{GID, PID, LID}, row::LID)::LID

Returns the current number of local entries in the given row
"""
function getNumEntriesInLocalRow end

"""
    getGlobalNumDiags(::RowGraph{GID, PID, LID})::GID

Returns the global number of diagonal entries
"""
function getGlobalNumDiags end

"""
    getLocalNumDiags(::RowGraph{GID, PID, LID})::LID

Returns the local number of diagonal entries
"""
function getLocalNumDiags end

"""
    getGlobalMaxNumRowEntries(::RowGraph{GID, PID, LID})::LID

Returns the maximum number of entries across all rows/columns on all processors
"""
function getGlobalMaxNumRowEntries end

"""
    getLocalMaxNumRowEntries(::RowGraph{GID, PID, LID})::LID

Returns the maximum number of entries across all rows/columns on the calling processor
"""
function getLocalMaxNumRowEntries end

"""
    isLowerTriangular(::RowGraph{GID, PID, LID})::Bool

Whether the graph is lower trianguluar
"""
function isLowerTriangular end

"""
    isUpperTriangular(::RowGraph{GID, PID, LID})::Bool

Whether the graph is upper trianguluar
"""
function isUpperTriangular end

"""
    isGloballyIndexed(::RowGraph)::Bool

Whether the graph is using global indices
"""
function isGloballyIndexed end


text:
using MonteCarloX
using Random
using Distributions
using StatsBase
using Printf
using ProgressMeter
using LinearAlgebra
using DelimitedFiles
using StaticArrays

"""
Reproduce results similar to Fig.2 of Lloyed-Smith et al., Nature (2005) on the
the probability of exctinction (here after 3 weeks for mu=1/8) and the
occurrance of the first outbreak (here first day with >50 cases per trajectory).
"""
function reproduce_lloyed_smith_nature_2005_fig2(path)
    list_ratio_lambda_mu = collect(0.5:0.1:3)
    list_k = [10, 1, 1e-1, 1e-2]
    list_a = [3, 5, 10]
    N = Int(1e5)
    I0 = 1
    R0 = 0
    S0 = N - I0
    epsilon=0
    range_seed = 1:Int(1e6);

    # motivated by COVID-19 parameters [cf. Dehning et al, Science (2020)]
    mu = 1/8. # avg. infectious time of 5+3 days
    time_meas = 50 # days of simulation
    dT = 1

    for ratio_lambda_mu in list_ratio_lambda_mu
        lambda_0 = ratio_lambda_mu*mu
        list_P = []
        list_P_name = []
        for k in list_k
            P_Gamma = Gamma(k, lambda_0/k)
            push!(list_P, P_Gamma)
            push!(list_P_name, @sprintf("Gamma_k%.2e",k))
        end
        for a in list_a
            lambda_1 = mu/a
            lambda_2 = mu*a
            prior_1 = (lambda_0 - lambda_2)/(lambda_1 - lambda_2)
            prior_2 = 1 - prior_1
            if lambda_2 >= lambda_0
                P_Bimodal = MixtureModel(Normal[Normal(lambda_1, 0), Normal(lambda_2, 0)], [prior_1, prior_2])
                push!(list_P, P_Bimodal)
                push!(list_P_name, @sprintf("Bimodal_a%.2e",a))
            end
        end

        # simulations
        for (P,P_name) in zip(list_P, list_P_name)
            println("R=", ratio_lambda_mu, " P=", P_name)

            list_time = collect(range(0, time_meas, step=dT))
            trajectory_S = zeros(Int,length(list_time))
            trajectory_I = zeros(Int,length(list_time))
            trajectory_R = zeros(Int,length(list_time))

            # observables
            probability_exctinction_per_day = zeros(length(list_time))
            probability_first_day_more_than_50 = zeros(length(list_time))

            #generate and analyse trajectories
            @showprogress 1 for seed in range_seed
                rng = MersenneTwister(seed);
                system = SIR{typeof(P)}(rng, P, epsilon, mu, S0, I0, R0)
                trajectory!(rng, list_time, trajectory_S, trajectory_I, trajectory_R, system)

                first_index_more_than_50 = findfirst(x->x>50, trajectory_I)
                if !(typeof(first_index_more_than_50) == Nothing)
                    probability_first_day_more_than_50[first_index_more_than_50] += 1
                end

                probability_exctinction_per_day[trajectory_I .== 0] .+= 1
            end

            probability_exctinction_per_day ./= length(range_seed)
            probability_first_day_more_than_50 ./= length(range_seed)

            #write out
            filename =  @sprintf("%s/reproduce_lloyed_smith_nature_2005_fig2_%s_R%.2e_mu%.2e_N%.2e_I0%.2e.dat",
                                 path, P_name, ratio_lambda_mu, mu, N, I0)

            open(filename; write=true) do f
                write(f, "#P_ext(t) - probability of trajectory to be extinct at this day\n")
                write(f, "#P_50(t)  - probability that trajectory has more than 50 infections for the first time on this day\n")
                write(f, "#t\t P_ext(t)\t N_50\n")
                writedlm(f, [list_time  probability_exctinction_per_day probability_first_day_more_than_50])
            end
        end
    end
end

function example_dist_sir_hetero_rates(path;
                                       ratio_lambda_mu=1.0)
    # define default values (some specific initial conditions, I0=100)
    N = Int(1e5)
    list_I0 = Int[1,10,100]
    seeds = collect(Int, 1:1e6);

    # motivated by COVID-19 parameters [cf. Dehning et al, Science (2020)]
    mu = 1/8. # avg. infectious time of 5+3 days
    lambda_0 = ratio_lambda_mu*mu # focus on the "critical point" of sustained activity
    sigma = mu/10
    time_meas = 21 # 1 week observation of effect of initial perturbation
    list_k = [1,1e-1,1e-2]
    list_alpha = [1.5,2,3] # mean only defined for alpha > 1

    # Different distributions with same mean speading rate lambda_0
    list_P = []
    list_P_name = []
    # Compared to [Lloyd-Smith et al., Science (2020)] delta-distribution for
    # differential-equation model corresponds already to an exponential
    # distribution of R_0 i.e. a geometric offspring distribution (k=1)
    P_Delta = Normal(lambda_0, 0.0)
    push!(list_P, P_Delta)
    push!(list_P_name, "Delta")
    P_Gaussian = Normal(lambda_0, sigma)
    push!(list_P, P_Gaussian)
    push!(list_P_name, "Gaussian")
    # We cannot fully reproduce the negative-binomial offspring distribution by
    # a Gamma-Poisson mixture because we have no generation based model here.
    # Still, we should approximate the same behavior by choosing a Gamma
    # distribution of spreading rates with the same mean
    # https://en.wikipedia.org/wiki/Gamma_distribution
    #   mean(Gamma) = k*theta = lambda_0 i.e. theta = lambda_0/k
    for k in list_k
        P_Gamma = Gamma(k, lambda_0/k)
        push!(list_P, P_Gamma)
        push!(list_P_name, @sprintf("Gamma_k%.2e",k))
    end
    # Pareto distribution for scale-free spreading rates (motivated from
    # scale-free-ness of social networks)
    # P_Pareto(x_m, alpha) = GeneralizedPareto(x_m, x_m/alpha, 1/alpha)
    # for Pareto to give mean=lambda_0 : x_m(lambda_0,alpha) = (alpha-1)lambda_0/alpha
    for alpha in list_alpha
        x_m = (alpha-1)*lambda_0/alpha
        P_Pareto = GeneralizedPareto(x_m, x_m/alpha, 1/alpha)
        push!(list_P, P_Pareto)
        push!(list_P_name, @sprintf("Pareto_alpha%.2e",alpha))
    end
    # mean(lambda) = lambda_1*prior_1 + lambda_2*prior_2 = lambda_0
    # prior_1 + prior_2 = 1
    lambda_1 =  0.0
    lambda_2 = 10.0*mu
    prior_1 = (lambda_0 - lambda_2)/(lambda_1 - lambda_2)
    prior_2 = 1 - prior_1
    P_Bimodal = MixtureModel(Normal[Normal(lambda_1, 0), Normal(lambda_2, 0)], [prior_1, prior_2])
    push!(list_P, P_Bimodal)
    push!(list_P_name, @sprintf("Bimodal_left%.2e_right%.2e",lambda_1, lambda_2))
    # simulations
    for I0 in list_I0
        R0 = 0
        S0 = N-I0
        for (P,P_name) in zip(list_P, list_P_name)
            println(N, " ", I0, " ", P_name)
            list_time, list_dist = distribution(P, S0, I0, R0, time_meas, seeds, mu=mu);
            filename_base =  @sprintf("%s/distribution_%s_lambda%.2e_mu%.2e_N%.2e_I0%.2e",
                                 path, P_name, lambda_0, mu, N, I0)
            write_distributions(list_time, list_dist, filename_base)
        end
    end
end

function write_distributions(list_time, list_dist, filename_base)
    for (time,dist) in zip(list_time, list_dist)
        filename = @sprintf("%s_T%.2e.dat",filename_base,time)
        open(filename; write=true) do f
            write(f, "#I\t P(I)\n")
            writedlm(f, zip(collect(dist.edges[1]),dist.weights))
        end
    end
end

function sir_hetero_rates(P_lambda::D, S0, I0, R0, time_meas, seed;
                          mu=1.0, epsilon=0.0, dT=1) where D
    rng = MersenneTwister(seed);
    system = SIR{D}(rng, P_lambda, epsilon, mu, S0, I0, R0)

    list_T = collect(range(0, time_meas, step=dT))
    list_S = zeros(length(list_T));
    list_I = zeros(length(list_T));
    list_R = zeros(length(list_T));
    trajectory!(rng, list_T, list_S, list_I, list_R, system)

    return list_T, list_S, list_I, list_R, system
end

# memory-wise this is not the most efficient solution
# Alternative is to have one trajectory readout and directly add to
# distribution (which requires some a priori size though)
function distribution(P_lambda::D, S0, I0, R0, time_meas, seeds;
                      mu=1.0, epsilon=0, dT=1) where D

    list_T = collect(range(0, time_meas, step=dT))
    trajectories_S = [zeros(Int,length(list_T)) for i=1:length(seeds)]
    trajectories_I = [zeros(Int,length(list_T)) for i=1:length(seeds)]
    trajectories_R = [zeros(Int,length(list_T)) for i=1:length(seeds)]
    @showprogress 1 for (i,seed) in enumerate(seeds)
        rng = MersenneTwister(seed);
        system = SIR{D}(rng, P_lambda, epsilon, mu, S0, I0, R0)
        trajectory!(rng, list_T, trajectories_S[i], trajectories_I[i], trajectories_R[i], system)
    end
    #index [time,trajectory] -> a way to make this [trajectory,time] to use later trajectories[:,t]
    trajectories_S = hcat(trajectories_S...);
    trajectories_I = hcat(trajectories_I...);
    trajectories_R = hcat(trajectories_R...);

    # then fit distributions (for now only I)
    # TODO: generalize (function for this) to all cases
    list_dist = []
    for t=1:size(trajectories_I,1)
        min_I = 0
        max_I = maximum(trajectories_I[t,:])+2
        step = 1
        hist = fit(Histogram, trajectories_I[t,:], min_I:step:max_I)
        dist = normalize!(float(hist))
        push!(list_dist, dist)
    end

    return list_T, list_dist
end

"""
# generate a stochastic trajectory
SIR dynamics is goverened by differential equation
``
    \\frac{dI}{dt} = \\sum_{i=1}^{I}\\lambda_i\\frac{S}{N} - \\mu I
``

# Arguments
* distribution for the heterogeneous rates P(lambda)
"""
function trajectory!(rng, list_T, list_S, list_I, list_R, system)
    rates = current_rates(system)
    pass_update!(rates, index) = update!(rates, index, system, rng)

    list_S[1] = system.measure_S
    list_I[1] = system.measure_I
    list_R[1] = system.measure_R
    time_simulation = Float64(list_T[1])
    for i in 2:length(list_T)
        if time_simulation < list_T[i]
            dT = list_T[i] - time_simulation
            dT_sim = advance!(KineticMonteCarlo(), rng, rates, pass_update!, dT)
            time_simulation += dT_sim
        end
        list_S[i] = system.measure_S
        list_I[i] = system.measure_I
        list_R[i] = system.measure_R
    end
end

mutable struct SIR{D}
    epsilon::Float64
    mu::Float64
    P_lambda::D
    current_lambda::Vector{Float64}
    sum_current_lambda::Float64
    update_current_lambda::Int
    N::Int
    S::Int
    I::Int
    R::Int
    measure_S::Int
    measure_I::Int
    measure_R::Int

    function SIR{D}(rng::AbstractRNG, P_lambda::D, epsilon, mu, S0::Int, I0::Int, R0::Int) where D
        N = S0 + I0 + R0
        current_lambda = rand(rng, P_lambda, I0)
        sum_current_lambda = sum(current_lambda)
        new(epsilon, mu, P_lambda, current_lambda, sum_current_lambda, 0, N, S0, I0, R0, S0, I0, R0)
    end
end

"""
update rates and system according to the last event that happend:
recovery:  1
infection: 2
"""
function update!(rates::AbstractVector, index::Int, system::SIR, rng::AbstractRNG)
    system.measure_S = system.S
    system.measure_I = system.I
    system.measure_R = system.R
    if index == 1 # recovery
        random_I = rand(rng, 1:system.I)
        delete_from_system!(system, random_I)
        system.I -= 1
        system.R += 1
    elseif index == 2 # infection
        add_to_system!(system, rand(rng, system.P_lambda))
        system.S -= 1
        system.I += 1
        @assert length(system.current_lambda) == system.I
    elseif index == 0 # absorbing state of zero infected
        system.measure_S = system.S
        system.measure_I = system.I
        system.measure_R = system.R
    else
        throw(UndefVarError(:index))
    end

    rates .= current_rates(system)
end

function delete_from_system!(system, index)
    update_sum!(system, -system.current_lambda[index])
    deleteat!(system.current_lambda, index)
end

function add_to_system!(system, lambda)
    push!(system.current_lambda, lambda)
    update_sum!(system, lambda)
end

function update_sum!(system, change)
    system.sum_current_lambda += change
    system.update_current_lambda += 1
    if system.update_current_lambda > 1e5
        system.sum_current_lambda = sum(system.current_lambda)
        system.update_current_lambda = 0
    end
end

"""
evaluate current rates of SIR system
"""
function current_rates(system::SIR)
    rate_recovery = system.mu * system.I
    rate_infection = system.sum_current_lambda* system.S/system.N + system.epsilon
    return MVector{2,Float64}(rate_recovery, rate_infection)
end


text:
using PyCall

cd("/Volumes/GoogleDrive/My Drive/Consulting_Projects/CATF WRA Arizona/AZ_Modeling/AZ-22x5/")

# I don't think this is necessary, but try using if the script doesn't work.
pushfirst!(PyVector(pyimport("sys")."path"), "")

@pyimport importlib.machinery as machinery
loader = machinery.SourceFileLoader("setup_2045_cases","setup_2045_cases.py")
my_mod = loader[:load_module]("setup_2045_cases")
my_mod.transfer_2030_results()

max_stars_repo_name: gschivley/GenX-helpers

max_issues_repo_name: gschivley/GenX-helpers

max_forks_repo_name: gschivley/GenX-helpers


text:
module BnpUtil

@doc raw"""
Stick-breaking function.

Arguments
=========
- `v`: A vector of length $K - 1$, where $K > 1$.

Return
======
- A simplex (w) of dimension $K$. Where ∑ₖ wₖ = 1, and each wₖ ≥ 0.
"""
function stickbreak(v)
    K = length(v) + 1
    cumprod_one_minus_v = cumprod(1 .- v)

    eta = [if k == 1
               v[1]
           elseif k == K
               cumprod_one_minus_v[K-1]
           else
               v[k] * cumprod_one_minus_v[k-1]
           end
           for k in 1:K]

    return eta
end

end  # End of module BnpUtil

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
# ---
# title: 1320. Minimum Distance to Type a Word Using Two Fingers
# id: problem1320
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-distance-to-type-a-word-using-two-fingers/description/>
# hidden: true
# ---
#
# ![](https://assets.leetcode.com/uploads/2020/01/02/leetcode_keyboard.png)
#
# You have a keyboard layout as shown above in the XY plane, where each English
# uppercase letter is located at some coordinate, for example, the letter **A**
# is located at coordinate **(0,0)** , the letter **B** is located at coordinate
# **(0,1)** , the letter **P** is located at coordinate **(2,3)** and the letter
# **Z** is located at coordinate **(4,1)**.
#
# Given the string `word`, return the minimum total distance to type such string
# using only two fingers. The distance between coordinates **(x 1,y1)** and **(x
# 2,y2)** is **|x 1 \- x2| + |y1 \- y2|**.
#
# Note that the initial positions of your two fingers are considered free so
# don't count towards your total distance, also your two fingers do not have to
# start at the first letter or the first two letters.
#
#
#
# **Example 1:**
#
#
#
#     Input: word = "CAKE"
#     Output: 3
#     Explanation: Using two fingers, one optimal way to type "CAKE" is:
#     Finger 1 on letter 'C' -> cost = 0
#     Finger 1 on letter 'A' -> cost = Distance from letter 'C' to letter 'A' = 2
#     Finger 2 on letter 'K' -> cost = 0
#     Finger 2 on letter 'E' -> cost = Distance from letter 'K' to letter 'E' = 1
#     Total distance = 3
#
#
# **Example 2:**
#
#
#
#     Input: word = "HAPPY"
#     Output: 6
#     Explanation:
#     Using two fingers, one optimal way to type "HAPPY" is:
#     Finger 1 on letter 'H' -> cost = 0
#     Finger 1 on letter 'A' -> cost = Distance from letter 'H' to letter 'A' = 2
#     Finger 2 on letter 'P' -> cost = 0
#     Finger 2 on letter 'P' -> cost = Distance from letter 'P' to letter 'P' = 0
#     Finger 1 on letter 'Y' -> cost = Distance from letter 'A' to letter 'Y' = 4
#     Total distance = 6
#
#
# **Example 3:**
#
#
#
#     Input: word = "NEW"
#     Output: 3
#
#
# **Example 4:**
#
#
#
#     Input: word = "YEAR"
#     Output: 7
#
#
#
#
# **Constraints:**
#
#   * `2 <= word.length <= 300`
#   * Each `word[i]` is an English uppercase letter.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "hmm" begin
    @testset "homopolymer_filter" begin
        seqs = ["GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC",
                "GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC",
                "GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATCAAAAAAAAAAAAAAAAAAAAA"]
        names = ["s1", "s2", "s3"]
        phreds = [[UInt8(30) for i in 1:length(seq)] for seq in seqs]
        newseqs, newphreds, newnames = homopolymer_filter(seqs, phreds, names)

        @test newseqs == ["GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC",
                         "GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC",
                         "GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC"]
    end

end


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

__precompile__(true)

"""
    BayesianBlocks

Bayesian Blocks algorithm described in [Scargle 2021](https://doi.org/10.1088/0004-637X/764/2/167).
"""
module BayesianBlocks

include("core.jl")

end # module


text:
Pkg.activate()

Pkg.add("UrlDownload")
Pkg.add("DataFrames")
Pkg.add("PyCall")

using DataFrames
using CSV
using PyPlot
using PyCall

# https://github.com/JuliaPy/PyPlot.jl
PyPlot.matplotlib.style.use("ggplot")
#rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
#rcParams["font.size"] = 15

#if you want to keep the Python True/False syntax
#const True = true
#const False = false

countries_regex = r"Belgium|Chile|France|Germany|Norway|Sweden|World"
countries_list = ["Belgium", "Chile", "France", "Germany", "Norway", "Sweden", "World"]
colors = ["teal", "lightseagreen", "skyblue", "steelblue", "olivedrab", "forestgreen", "cornflowerblue"]

# ## Time series
# ### Coverage

#raw_data = urldownload("https://raw.githubusercontent.com/gd1989/ECP/master/Total_economy/Coverage/Total_coverage.csv?token=AD3IRDAVNXIUQ6ER3UJU5WS7UQL5K")
path_cov = "/Users/GD/Documents/GitHub/ECP/Total_economy/Coverage/Total_coverage.csv"

coverage = CSV.read(path_cov, DataFrame)
select!(coverage, ["Jurisdiction", "Year", "cov_tax_ets_share_jurGHG"])

coverage.cov_perc = coverage.cov_tax_ets_share_jurGHG*100

figure(figsize=(12, 10))

i = 1
for ctry in countries_list
        temp = coverage[coverage[:, :Jurisdiction] .== ctry, :]
        plot(temp.Year, temp.cov_perc, color=colors[i], label=ctry)
        global i += 1
end

xlabel("Time", fontsize=14)
ylabel("Coverage (% GHG)", fontsize=14)
yticks(fontsize=12)
xlim(1990, 2017)
ylim(0)
legend(loc="upper center", bbox_to_anchor=(0.5, -0.05), fancybox=true, shadow=true, ncol=4, fontsize=14)
plt.tight_layout()

plt.savefig("/Users/GD/Documents/GitHub/ECP/page/_assets/coverage_ts.png")
plt.close()

# ### ECP
path_ecp = "/Users/GD/Documents/GitHub/ECP/Total_economy/ECP_fixed_weights/ecp_2015.csv"

ecp = CSV.read(path_ecp, DataFrame)

figure(figsize=(12, 10))

i = 1
for ctry in countries_list
        temp = ecp[ecp[:, :Jurisdiction] .== ctry, :]
        plot(temp.Year, temp.ECP_tax_ets_jurGHG_2019USD, color=colors[i], label=ctry)
        global i += 1
end

xlabel("Time", fontsize=14)
ylabel("Carbon price (2019USD/tCO2)", fontsize=14)
yticks(fontsize=12)
xlim(1990, 2018)
ylim(0)

legend(loc="upper center", bbox_to_anchor=(0.5, -0.05), fancybox=true, shadow=true, ncol=4, fontsize=14)
plt.tight_layout()

plt.savefig("/Users/GD/Documents/GitHub/ECP/page/_assets/ecp_ts.png")
plt.close()

# ## Distribution
# ### Distribution of coverage in 2016

cov_dist = coverage[coverage[:, :Year] .== 2017, :]
dropmissing!(cov_dist)

fig, (ax, ax2) = plt.subplots(2,1, sharex=true) # make the axes

ax.hist(cov_dist.cov_tax_ets_share_jurGHG, color="blue", edgecolor="white",
        bins=20,
        density=false, alpha=0.5)
ax2.hist(cov_dist.cov_tax_ets_share_jurGHG, color="blue", edgecolor="white",
         bins=20,
         density=false, alpha=0.5)

ax.set_ylim([150,170]) # numbers here are specific to this example
ax2.set_ylim([0,10]) # numbers here are specific to this example

ax.spines["bottom"].set_visible(false)
ax2.spines["top"].set_visible(false)

ax.xaxis.tick_bottom()
ax.tick_params(labeltop=false)
ax2.xaxis.tick_bottom()

#plt.axvline(mean_age, label = "Mean", color = "red")
#plt.axvline(median_age, label = "Median")

plt.legend()

plt.savefig("/Users/GD/Documents/GitHub/ECP/page/_assets/coverage_hist.png")
plt.close()

max_stars_repo_name: g-dolphin/ECP

max_issues_repo_name: g-dolphin/ECP

max_forks_repo_name: g-dolphin/ECP


text:
using Rotations, StaticArrays, Test
using Unitful

@testset "2d Rotations" begin

    #################################
    # Traits, types, and construction
    #################################

    @testset "Core" begin
        r = one(RotMatrix{2,Float32})
        @test RotMatrix((1,0,0,1)) == RotMatrix(@SMatrix [1 0; 0 1])
        @test_throws DimensionMismatch RotMatrix((1,0,0,1,0))
    end

    @testset "Unitful" begin
        # Make sure rotations created from unitful angles
        # don't extraneously contain those units (see issue #55)
        @test eltype(Angle2d(10u"°")) <: Real
        @test eltype(Angle2d(20u"rad")) <: Real
    end

    ###############################
    # Check fixed relationships
    ###############################

    @testset "Identity rotation checks" begin
        I = one(SMatrix{2,2,Float64})
        I32 = one(SMatrix{2,2,Float32})
        for R in [RotMatrix{2}, Angle2d]
            @test @inferred(size(R)) == (2,2)
            @test @inferred(size(R{Float32})) == (2,2)
            @test one(R)::R == I
            @test one(R{Float32})::R{Float32} == I32
        end
    end

    ###############################
    # Check zero function
    ###############################

    @testset "zero checks" begin
        for R in (RotMatrix{2}, Angle2d)
            # zero
            @test zero(R) == zero(R{Float64}) == zero(one(R))
            @test zero(R) isa SMatrix
            @test zero(R{Float64}) isa SMatrix
            @test zero(one(R)) isa SMatrix
            # zeros
            @test zeros(R)[1] == zeros(R,3)[1] == zeros(R,3,3)[1] == zeros(R,(3,3,3))[1] == zero(R)
            @test zeros(R) isa Array{<:SMatrix,0}
            @test zeros(R,3) isa Array{<:SMatrix,1}
            @test zeros(R,3,3) isa Array{<:SMatrix,2}
            @test zeros(R,(3,3,3)) isa Array{<:SMatrix,3}
        end
    end

    ################################
    # check on the inverse function
    ################################

    @testset "Testing inv()" begin
        repeats = 100
        for R in [RotMatrix{2,Float64}, Angle2d{Float64}]
            I = one(R)
            Random.seed!(0)
            for i = 1:repeats
                r = rand(R)
                @test isrotation(r)
                @test inv(r) == adjoint(r)
                @test inv(r) == transpose(r)
                @test inv(r)*r ≈ I
                @test r*inv(r) ≈ I
            end
        end
    end

    ################################
    # check on the norm functions
    ################################

    @testset "Testing norm() and normalize()" begin
        repeats = 100
        for R in [RotMatrix{2,Float64}, Angle2d{Float64}]
            I = one(R)
            Random.seed!(0)
            for i = 1:repeats
                r = rand(R)
                @test norm(r) ≈ norm(Matrix(r))
                if VERSION ≥ v"1.5"
                    @test normalize(r) ≈ normalize(Matrix(r))
                    @test normalize(r) isa SMatrix
                end
            end
        end
    end

    #########################################################################
    # Rotate some stuff
    #########################################################################

    # a random rotation of a random point
    @testset "Rotate Points" begin
        repeats = 100
        for R in [RotMatrix{2}, Angle2d]
            Random.seed!(0)
            for i = 1:repeats
                r = rand(R)
                m = SMatrix(r)
                v = randn(SVector{2})

                @test r*v ≈ m*v
            end

            # Test Base.Vector also
            r = rand(R)
            m = SMatrix(r)
            v = randn(2)

            @test r*v ≈ m*v
        end
    end

    # a random rotation of a random unitful point
    @testset "Rotate Unitful Points" begin
        repeats = 100
        for R in [RotMatrix{2}, Angle2d]
            Random.seed!(0)
            for i = 1:repeats
                r = rand(R)
                m = SMatrix(r)
                v = randn(SVector{2}) * u"m"

                @test r*v ≈ m*v
                @test eltype(r*v) <: Unitful.Length
                @test eltype(m*v) <: Unitful.Length
            end

            # Test Base.Vector also
            r = rand(R)
            m = SMatrix(r)
            v = randn(2) * u"m"

            @test r*v ≈ m*v
            @test eltype(r*v) <: Unitful.Length
            @test eltype(m*v) <: Unitful.Length
        end
    end

    # compose two random rotations
    @testset "Compose rotations" begin
        repeats = 100
        for R in [RotMatrix{2}, Angle2d]
            Random.seed!(0)
            for i = 1:repeats
                r1 = rand(R)
                m1 = SMatrix(r1)

                r2 = rand(R)
                m2 = SMatrix(r2)

                @test r1*r2 ≈ m1*m2
                θ1, θ2 = atan(r1[2,1],r1[1,1]), atan(r2[2,1],r2[1,1])
                @test r1*r2 ≈ RotMatrix(θ1+θ2)
                @test r1/r2 ≈ RotMatrix(θ1-θ2)
                @test r1\r2 ≈ RotMatrix(θ2-θ1)
            end
        end
    end

    #########################################################################
    # Test conversions between rotation types
    #########################################################################
    @testset "Convert rotations" begin
        repeats = 100
        for R in [RotMatrix{2}, Angle2d]
            Random.seed!(0)
            for _ in 1:repeats
                r1 = rand(R)
                @test R(r1) == r1
            end
        end
    end

    @testset "Types and products" begin
        for (R,T) in ((RotMatrix(pi/4), Float64),
                      (RotMatrix(Float32(pi/4)), Float32),
                      (RotMatrix{2,Float32}(pi/4), Float32))
            @test eltype(R) == T
            @test size(R) == (2,2)
            @test R == T[cos(pi/4) -sin(pi/4); sin(pi/4) cos(pi/4)]
            @test R * R ≈ T[0 -1; 1 0]
        end
    end

    @testset "RotMatrix{2} vs Angle2d" begin
        repeats = 100
        Random.seed!(0)
        for _ in 1:repeats
            theta = randn()
            r1 = RotMatrix{2}(theta)
            r2 = Angle2d(theta)
            v = randn(2)
            @test r1 ≈ r2
            @test r1 * v ≈ r2 * v
            @test RotMatrix{2}(r2) ≈ r1
            @test Angle2d(r1) ≈ r2
        end
    end

    @testset "Angle2d" begin
        repeats = 100
        Random.seed!(0)
        for _ in 1:repeats
            r = rand(Angle2d)
            @test r^(-1) ≈ inv(r)
            @test r^(-1.0) ≈ inv(r)
            @test r^1 ≈ r
            @test r^1.0 ≈ r
            @test r^2 ≈ r*r
            @test r^2.0 ≈ r*r

            x = randn()
            y = randn()
            @test r^(x+y) ≈ r^x * r^y
        end
    end

    @testset "angle" begin
        repeats = 100
        Random.seed!(0)
        for _ in 1:repeats
            theta = 2pi*rand() - pi
            a2d = Angle2d(theta)
            r2 = RotMatrix{2}(theta)
            @test rotation_angle(r2) ≈ theta
            @test rotation_angle(a2d) == theta
            @test rotation_angle(a2d) == Rotations.params(a2d)[1]
        end
    end
end

nothing


text:
@reexport using DataFrames
@reexport using PrettyTables

# See https://dataframes.juliadata.org/stable/lib/functions/#Base.show
# and https://ronisbr.github.io/PrettyTables.jl/stable/man/usage/
printsimple(df::DataFrame; kw...) = show(
    df;
    summary = false,
    eltypes = false,
    show_row_number = false,
    formatters = ft_printf("%.3g"),
    alignment = :l,
    kw...
)


text:
#####
##### Benchmark.jl
##### Script to run a new benchmark on a previous training session
##### (requires training with option `save_intermediate=true`)
#####

using AlphaZero

include("games.jl")
GAME = get(ENV, "GAME", "connect-four")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

SESSION_DIR = "sessions/$GAME"

# Examples: run a network-only benchmark

netonly = Benchmark.NetworkOnly(τ=0.5, use_gpu=true)
alphazero = Benchmark.Full(Training.arena.mcts)
baselines = [Training.mcts_baseline, Training.minmax_baseline]

make_duel(player, baseline) =
  Benchmark.Duel(
    player,
    baseline,
    num_games=200,
    reset_every=40,
    flip_probability=0.5,
    color_policy=CONTENDER_WHITE)

netonly_benchmark   = [make_duel(netonly, b) for b in baselines]
alphazero_benchmark = [make_duel(alphazero, b) for b in baselines]

benchmark = alphazero_benchmark
label = "full"

AlphaZero.UserInterface.run_new_benchmark(
  Game, Training.Network{Game}, SESSION_DIR,
  label, benchmark,
  params=Training.params, itcmax=nothing)


text:
# Authors: Rajiv Kumar
#          Ali Siahkoohi, alisk@gatech.edu
# Date: February 2021
# Copyright: Georgia Institute of Technology, 2021

module WavefieldRecoveryUQ

using DrWatson
using Random
using HDF5
using DataFrames
using LinearAlgebra
using Distributions
using Statistics
using InvertibleNetworks
using PyPlot: Figure
using CUDA

import DrWatson: _wsave
import Distributions: logpdf, gradlogpdf
import Base.getindex
import Base.*

# Utilities
include("./utils/load_experiment.jl")
include("./utils/data_loader.jl")
include("./utils/mask.jl")
include("./utils/savefig.jl")
include("./utils/logpdf.jl")

# Objective functions
include("./objectives/objectives.jl")

end


text:
using StanSample, DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions,
      ModelingToolkit, RecursiveArrayTools, Distributions, Random, Test
using MCMCChains

# Uncomment for local testing only, make sure MCMCChains and StatsPlots are available
using MCMCChains, StatsPlots

#Random.seed!(123)
cd(@__DIR__)
isdir("tmp") && rm("tmp", recursive=true)

include("../../stan_string.jl")
include("../../stan_inference.jl")

println("\nFour parameter case\n")
f1 = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .5randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

priors = [truncated(Normal(1.5,0.5),0.5,2.5),truncated(Normal(1.2,0.5),0,2.0),
          truncated(Normal(1.0,0.5),1.0,4),truncated(Normal(1.0,0.5),0,2)]

@time bayesian_result = stan_inference(prob1,t,data,priors;
  num_samples=9000, nchains=1,
  num_warmup=1000,vars =(StanODEData(),InverseGamma(3,3)))

stanmodel = bayesian_result.model

@time bayesian_result = stan_inference(prob1,t,data,priors, stanmodel;
  num_samples=9000, nchains=1,
  num_warmup=1000,vars =(StanODEData(),InverseGamma(3,3)))

if success(bayesian_result.return_code)
  chns  = read_samples(bayesian_result.model, output_format=:mcmcchains)
  describe(chns) |> display

  read_summary(stanmodel) |> display

end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function testprojectonPatchBndbyRestriction(fn::Function, Nx::Int, M::Int, loc::Int)::Float64
    pfbnd = projectonPatchBndbyRestriction(fn, Nx, M, loc)
    efbnd = Float64[fn(x) for x in chebgrid(Nx, M, loc)]
    w     = chebweights(Nx)
    L2error = sqrt(w'*(pfbnd-efbnd).^2)
    return L2error
end

function testprojectonPatchbyRestriction(fn::Function, Nx::Int, Ny::Int, M::Int, loc::Array{Int,1})::Float64
    pfpatch = projectonPatchbyRestriction(fn, Nx, Ny, M, loc)
    efpatch = Float64[fn(x,y) for x in chebgrid(Nx, M, loc[1]), y in chebgrid(Ny, M, loc[2])]
    return L2norm(pfpatch, efpatch, chebweights(Nx), chebweights(Ny))
end

@test testprojectonPatchBndbyRestriction(x->x^8, 6, 4, 1) >  testprojectonPatchBndbyRestriction(x->x^8, 6, 8, 1)
@test testprojectonPatchbyRestriction((x,y)->x^8+y^9, 8, 4, 7, [2,1]) > testprojectonPatchbyRestriction((x,y)->x^8+y^9, 5, 9, 3, [2,1])
@test projectonPatchbyRestriction((x,y)->x+y, 2, 2, 2, [1,1])[:, end] ≈ projectonPatchbyRestriction((x,y)->x+y, 2, 2, 2, [1,2])[:, 1]
@test projectonPatchbyRestriction((x,y)->x+y, 2, 2, 2, [1,1])[end, :] ≈ projectonPatchbyRestriction((x,y)->x+y, 2, 2, 2, [2,1])[1, :]


text:
__precompile__(true)

module DndDB

	include("consoleutil.jl");
	include("enemydbcrud.jl");
	include("filehandling.jl");

end


text:
using ArrayFire
using BenchmarkTools

function pi_gpu(n)
    return 4 *
           sum(rand(AFArray{Float64}, n)^2 .+ rand(AFArray{Float64}, n)^2 .<= 1) ./
           n
end

@btime pi_gpu(10_000_000)

function pi_serial(n)
    inside = 0
    for i = 1:n
        x, y = rand(), rand()
        inside += (x^2 + y^2 <= 1)
    end
    return 4 * inside / n
end

@btime pi_serial(10_000_000)

#%% Mandelbrot fractal
using Plots
gr()

maxIterations = 500
gridSize = 2048
xlim = [-0.748766713922161, -0.748766707771757]
ylim = [0.123640844894862, 0.123640851045266]

x = range(xlim[1], xlim[2], length = gridSize)
y = range(ylim[1], ylim[2], length = gridSize)

xGrid = [i for i in x, j in y]
yGrid = [j for i in x, j in y]

z0 = xGrid + im * yGrid

function mandelbrotGPU(z0, maxIterations)
    z = z0
    count = ones(AFArray{Float32}, size(z))

    for n = 1:maxIterations
        z = z .* z .+ z0
        count = count + (abs(z) <= 2)
    end
    return sync(log(count))
end

count = Array(mandelbrotGPU(AFArray(z0), maxIterations))
heatmap(
    count,
    color = :inferno,
    grid = false,
    colorbar = :none,
    framestyle = :none,
    size = (2048, 2048),
)
savefig("mandelbrot.png")

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
# SAM Reader
# =========

mutable struct Reader{S <: TranscodingStream} <: BioGenerics.IO.AbstractReader
    state::State{S}
    header::Header
end

function Reader(state::State{S}) where {S <: TranscodingStream}

    rdr = Reader(state, Header())

    cs, ln, f = readheader!(rdr.state.stream, rdr.header, (sam_machine_header.start_state, rdr.state.linenum))

    rdr.state.state = sam_machine_body.start_state
    rdr.state.linenum = ln
    rdr.state.filled = false

    return rdr
end

"""
    SAM.Reader(input::IO)

Create a data reader of the SAM file format.

# Arguments
* `input`: data source
"""
function Reader(input::IO)

    if input isa TranscodingStream
        return Reader(State(input, 1, 1, false))
    end

    stream = TranscodingStreams.NoopStream(input)

    return Reader(State(stream, 1, 1, false))

end

function Base.eltype(::Type{<:Reader})
    return Record
end

function BioGenerics.IO.stream(reader::Reader)
    return reader.state.stream
end

"""
    header(reader::Reader)::Header

Get the header of `reader`.
"""
function header(reader::Reader)::Header
    return reader.header
end

function Base.close(reader::Reader)
    if reader.state.stream isa IO
        close(reader.state.stream)
    end
    return nothing
end

function index!(record::MetaInfo)
    stream = TranscodingStreams.NoopStream(IOBuffer(record.data))
    found = index!(stream, record)
    if !found
        throw(ArgumentError("invalid SAM metadata"))
    end
    return record
end

function index!(record::Record)
    stream = TranscodingStreams.NoopStream(IOBuffer(record.data))
    found = index!(stream, record)
    if !found
        throw(ArgumentError("invalid SAM record"))
    end
    return record
end

function Base.iterate(reader::Reader, nextone::Record = Record())
    if BioGenerics.IO.tryread!(reader, nextone) === nothing
        return nothing
    end
    return copy(nextone), empty!(nextone)
end

"""
    read!(rdr::Reader, rec::Record)

Read a `Record` into `rec`; overwriting or adding to existing field values.
It is assumed that `rec` is already initialized or empty.
"""
function Base.read!(rdr::Reader, record::Record)

    cs, ln, found = readrecord!(rdr.state.stream, record, (rdr.state.state, rdr.state.linenum))

    rdr.state.state = cs
    rdr.state.linenum = ln
    rdr.state.filled = found

    if found
        return record
    end

    if cs == 0 || eof(rdr.state.stream)
        throw(EOFError())
    end

    throw(ArgumentError("malformed SAM file"))
end


text:
include("Sorting.jl")

module R2

  importall Base
  importall Sorting

  export R2Object,println,
         CubeChange,*,==,isless,isless5,compare,compare5,
         indices2,indices5,
         Algorithm,mkAlgorithmFromAtomicMove,oneMoveAlgorithms,
         getAllHalfAlgs,getGroupedHalfAlgs!,categorizedWholeAlgs8,
         algs8, algsX

  ################################################################################

  abstract R2Object
  function println( o :: R2Object )
    print( o )
    println()
  end
  function print{T<:R2Object}( arr :: Array{T} )
    for i = 1:length(arr)
      print(arr[i])
      if i < length(arr)
        print(' ')
        #println()
      end
    end
  end
  function println{T<:R2Object}( arr :: Array{T} )
    for o = arr
      println(o)
    end
  end

  compare( x :: Int64, y :: Int64 ) = x<y ? -1 : x>y ? 1 : 0

  ################################################################################
  # Globals indices

  indices5 = [3,4,5,6,7]
  indices2 = [1,2]
  function setIndices2( i :: Int64, j :: Int64 )
    indices2 = [i,j]
    a = 1
    for b = 1:7
      if b!=i && b!=j
        indices5[a] = b
        a += 1
      end
    end
  end

  function compare5( x :: Array{Int64}, y :: Array{Int64} )
    for a = 1:5
      c = compare( x[indices5[a]], y[indices5[a]] )
      if c!=0
        return c
      end
    end
    return 0
  end

  ################################################################################
  # CubeChange
  ################################################################################

  type CubeChange <: R2Object
    pos :: Array{Int64} # permutation of 7 corners.
    rot :: Array{Int64} # rotation (after the permutation) of each corner.
  end

  # We will only print two-piece-affecting algorithms,
  # so this function relies on that assumption.
  function print( x :: CubeChange )
    print( "(" )

    for i = 1:7
      if x.rot[i] != 0
        print( " ", i )
        print( x.rot[i]==1 ? "t" : "t'")
      end
    end

    done = falses(7)
    for i = 1:7
      if !done[i] && x.pos[i] != i
        print( " ", i )
        done[i] = true
        j = x.pos[i]
        while j!=i
          print("->", j)
          done[j] = true
          j = x.pos[j]
        end
      else
        done[i] = true
      end
    end
    print( " )" )
  end

  idCubeChange = CubeChange( [1,2,3,4,5,6,7], [0,0,0,0,0,0,0] )
  centerCubeChange = CubeChange( [4,2,1,7,5,6,3], [0,0,0,0,0,0,0] )
  leftCubeChange = CubeChange( [2,6,3,1,5,4,7], [1,2,0,2,0,1,0] )
  rightCubeChange = CubeChange( [3,1,5,4,2,6,7], [2,1,1,0,2,0,0] )

  ==( x :: CubeChange, y :: CubeChange ) = x.pos == y.pos && x.rot == y.rot

  function compare( x :: CubeChange, y :: CubeChange )
    c = compare( x.pos, y.pos )
    if c==0
      return compare( x.rot, y.rot )
    end
    return c
  end
  function compare5( x :: CubeChange, y :: CubeChange )
    c = compare5( x.pos, y.pos )
    if c==0
      return compare5( permute!(copy(x.rot),x.pos), permute!(copy(y.rot),y.pos) )
    end
    return c
  end

  isless(x :: CubeChange, y :: CubeChange) = compare(x,y) < 0
  isless5(x :: CubeChange, y :: CubeChange) = compare5(x,y) < 0

  equal5(x :: CubeChange, y :: CubeChange) = compare(x,y) == 0

  # Do y, then do x
  function *( x :: CubeChange, y :: CubeChange)
    newPos = permute!( copy(y.pos), x.pos )
    newRot = ipermute!( copy(x.rot), y.pos )
    for i = 1:7
      newRot[i] = ( newRot[i] + y.rot[i] ) % 3
    end
    CubeChange( newPos, newRot )
  end

  # reverse a cube change
  function reverse( x :: CubeChange )
    newPos = invperm(x.pos)
    newRot = permute!( copy(x.rot), x.pos )
    for i = 1:7
      newRot[i] = (3 - newRot[i]) % 3
    end
    CubeChange( newPos, newRot )
  end

  # Whether B * A^-1 switches pieces or not, and whether it rotates or not.
  # 0 = do nothing
  # 1 = twist only
  # 2 = switch only
  # 3 = switch and a twist
  #function categorizePair( x :: CubeChange, y :: CubeChange )
  #  #p = x.pos == y.pos
  #  #r = x.rot == y.rot
  #  #return p ? (r ? 0 : 1) : r ? 2 : 3
  #end
  function categorize8( x :: CubeChange )
    p = true
    for i = 1:7
      if x.pos[i] != i
        p = false
        break
      end
    end
    r = true
    for i = 1:7
      if x.rot[i] > 0
        r = false
        break
      end
    end
    return (p ? 0 : 2) + (r ? 0 : 1)
  end
  function categorize15( x :: CubeChange )
    p = true
    for i = 1:7
      if x.pos[i] != i
        p = false
        break
      end
    end
    r = 0
    for i = 1:7
      if x.rot[i] > 0
        r = x.rot[i]
        break
      end
    end
    return (p ? 0 : 3) + r
  end

  ################################################################################
  # More globals about atomic moves

  atomicMovesValid = [1,2,3,4,5,6,7,8,9]
  atomicMoveInconveniences = [2,3,2,4,6,4,2,3,2]

  numAtomicMoves = length(atomicMovesValid)
  atomicMoveNames = [ "R ", "R2", "R'", "F ", "F2", "F'", "U ", "U2", "U'" ]
  atomicCubeChanges = [ centerCubeChange, centerCubeChange*centerCubeChange, reverse(centerCubeChange),
                        leftCubeChange, leftCubeChange*leftCubeChange, reverse(leftCubeChange),
                        rightCubeChange, rightCubeChange*rightCubeChange, reverse(rightCubeChange) ]
  function moveInconvenience( m :: Int64 )
    return atomicMoveInconveniences[m]
    return (m==5) ? 20 : (m==4 || m==6) ? 14 : (m%3==2) ? 3 : 2
  #  return (m%3==2) ? 3 : 2
  end
  function reverseMoveCode( m :: Int64 )
    # 1 switches with 3 and vice versa.
    # 4 switches with 6 and vice versa.
    # 7 switches with 9 and vice versa.
    return m%3==0 ? m-2 : m%3==1 ? m+2 : m
  end

  ################################################################################
  # Algorithm
  ################################################################################

  type Algorithm <: R2Object
    moveList :: Array{Int64} # from beginning to end
    cubeChange :: CubeChange
  end
  function print( x :: Algorithm )
    for m = x.moveList
      print( atomicMoveNames[m], " ")
    end
    print( " " )
    print( x.cubeChange )
    print( " ( ", length(x.moveList), " moves )" )
    println( " ( inconvenience level ", inconvenience(x), " )")
  end

  idAlgorithm = Algorithm( [], idCubeChange )

  ==( x :: Algorithm, y :: Algorithm ) = x.cubeChange == y.cubeChange

  function inconvenience( x :: Algorithm )
    s = 0
    for m = x.moveList
      s += moveInconvenience(m)
    end
    return s
  end
  compare( x :: Algorithm, y :: Algorithm ) = compare( x.cubeChange, y.cubeChange )
  compare5( x :: Algorithm, y :: Algorithm ) = compare5( x.cubeChange, y.cubeChange )
  compareInconvenience( x :: Algorithm, y :: Algorithm ) = compare( inconvenience(x), inconvenience(y) )

  isless( x :: Algorithm, y :: Algorithm ) = compare( x, y ) < 0
  isless5( x :: Algorithm, y :: Algorithm ) = compare5( x, y ) < 0
  equal5(x :: Algorithm, y :: Algorithm) = compare( x, y ) == 0

  mkAlgorithmFromAtomicMove( atomicMove :: Int64 ) = Algorithm( [atomicMove], atomicCubeChanges[atomicMove] )

  # Do y then do x
  *( x :: Algorithm, y :: Algorithm ) = Algorithm( vcat( y.moveList, x.moveList ), x.cubeChange * y.cubeChange )

  function reverse( x :: Algorithm )
    newMoveList = copy(x.moveList)
    for i = 1:length(x.moveList)
      newMoveList[length(x.moveList)+1-i] = reverseMoveCode( x.moveList[i] )
    end
    return Algorithm( newMoveList, reverse(x.cubeChange) )
  end

  categorize8( x :: Algorithm ) = categorize8( x.cubeChange )
  categorize15( x :: Algorithm ) = categorize15( x.cubeChange )
  #categorizePair( x :: Algorithm, y :: Algorithm ) = categorizePair( x.cubeChange, y.cubeChange )
  # B * A^-1
  function combinePair( x :: Algorithm, y :: Algorithm )
    return x * reverse(y)
  end

  ################################################################################
  # Generate 8 algorithms

  # Get a list of all algorithms possible by combining a half-alg from xs with a half-alg from ys.
  # Whole algorithms will be x*y so "do y then do x".
  function *( xs :: Array{Algorithm}, ys :: Array{Algorithm} )
    xys = Array{Algorithm}( length(xs) * length(ys) )
    i = 1
    for x = xs
      for y = ys
        xys[i] = x*y
        i += 1
      end
    end
    xys = vcat(xs,xys)
    mergeSort!(xys)
    return removeDuplicatesBy!( compareInconvenience, xys )
  end

  function incrementHalfAlgs( xs :: Array{Algorithm} )
    ret = Array{Algorithm}( length(xs) * numAtomicMoves )
    i = 1
    for x = xs
      for y = oneMoveAlgorithms
        ret[i] = x*y
        i += 1
      end
    end
    return ret
  end

  # Get a list of all algorithms possible of n moves or less. Except duplicates.
  #function getAllHalfAlgs( n :: Int64 )
  #  i = 1
  #  ret = vcat( [idAlgorithm], oneMoveAlgorithms )
  #  retLimited = oneMoveAlgorithms
  #  for i = 2:n
  #    retLimited = incrementHalfAlgs( retLimited )
  #    ret = vcat( ret, retLimited )
  #  end
  #  return ret
  #end

  # Get a list of all algorithms possible of n moves or less. Except duplicates.
  function getAllHalfAlgs( n :: Int64 )
    if n==1
      return vcat( [idAlgorithm], oneMoveAlgorithms )
    end
    if n%2==0
      hn = getAllHalfAlgs( div(n,2) )
      return hn*hn
    else
      hn = getAllHalfAlgs( div(n,2) )
      hn1 = hn * oneMoveAlgorithms
      return hn * hn1
    end
    #n1 = getAllHalfAlgs( n-1 )
    #return n1 * oneMoveAlgorithms
  end

  function getGroupedHalfAlgs!( halfAlgs :: Array{Algorithm}, pair1 :: Int64, pair2 :: Int64 )
    setIndices2( pair1, pair2 )
    hag = groupBy!( isless5, halfAlgs )
    filteredHag = filter( g -> length(g)>1, hag )
  end

  function algorithmShorter( x :: Algorithm, y :: Algorithm )
    #s = length(x.moveList) < length(y.moveList)
    #if s
    #  return s
    #end
    return inconvenience(x) < inconvenience(y)
  end

  function getCategorizedWholeAlgs( halfAlgs :: Array{Algorithm},
                                    numCategories :: Int64, categorizer :: Function,
                                    pair1 :: Int64, pair2 :: Int64 )
    ret = Array{Algorithm}( numCategories )
    gs = getGroupedHalfAlgs!( halfAlgs, pair1, pair2 )

    for g = gs
      m = length(g)
      for i = 1:(length(g)-1)
        a = g[i]
        for j = (i+1):length(g)
          b = g[j]
          ab = combinePair(a,b)
          category = categorizer( ab )
          if category > 0
            if !isdefined( ret, category ) || algorithmShorter( ab, ret[category] )
              ret[category] = ab
            end
          end
        end
      end
    end
    return ret
  end

  function categorizedWholeAlgs( halfN :: Int64, numCategories :: Int64, categorizer :: Function )
    ret = Array{Algorithm}(3,3)

    halfAlgs = getAllHalfAlgs( halfN )

    for pairCode = 1:3
      pair1 = pairCode == 3 ? 2 : 1
      pair2 = pairCode == 1 ? 2 : 7
      algs3 = getCategorizedWholeAlgs( halfAlgs, 3, categorize8, pair1, pair2 )

      # pick shortest from each category.
      for i = 1:3
        if isdefined( algs3, i )
          ret[pairCode,i] = algs3[i]
        end
      end
    end
    return ret
  end

  function algsX_helper( halfN :: Int64, numCategories :: Int64, categorizer :: Function )
    halfAlgs = getAllHalfAlgs( halfN )
    for pair1 = 1:6
      for pair2 = (pair1+1):7
        wholeAlgs = getCategorizedWholeAlgs( halfAlgs, numCategories, categorizer, pair1, pair2 )
        for i = 1:length(wholeAlgs)
          printAlgIfDefined( wholeAlgs, i )
        end
      end
    end
  end
  function algsX( halfN :: Int64, categorizerCode :: Int64 )
    numCategories = categorizerCode == 1 ? 3 : 5
    categorizer = categorizerCode == 1 ? categorize8 : categorize15
    algsX_helper( halfN, numCategories, categorizer )
  end

  function printAlgIfDefined( algsArr :: Array{Algorithm,1}, i :: Int64 )
    if isdefined( algsArr, i )
      print( algsArr[i] )
    else
      println( "NOT FOUND" )
    end
  end
  function printAlgIfDefined( algsArr :: Array{Algorithm,2}, i :: Int64, j :: Int64 )
    if isdefined( algsArr, i, j )
      print( algsArr[i,j] )
    else
      println( "NOT FOUND" )
    end
  end

  function printCategorizedAlg( name, algsArr :: Array{Algorithm,2}, i :: Int64, j :: Int64 )
    print(name, ":  ")
    if isdefined( algsArr, i, j )
      print( algsArr[i,j] )
    else
      println( "NOT FOUND" )
    end
  end

  # halfN = 6 finds all the switching ones...
  function algs8( halfN :: Int64 )
    algs33 = categorizedWholeAlgs( halfN, 3, categorize8 )

    printCategorizedAlg( "AR ", algs33, 1, 1 )
    printCategorizedAlg( "AS1", algs33, 1, 3 )
    printCategorizedAlg( "AS2", algs33, 1, 2 )

    printCategorizedAlg( "DR ", algs33, 2, 1 )
    printCategorizedAlg( "DS1", algs33, 2, 2 )
    printCategorizedAlg( "DS2", algs33, 2, 3 )

    printCategorizedAlg( "OR ", algs33, 3, 1 )
    printCategorizedAlg( "OS1", algs33, 3, 2 )
    printCategorizedAlg( "OS2", algs33, 3, 3 )

    #algs33

  end

  ################################################################################
  # More globals about algorithms

  oneMoveAlgorithms = Array{Algorithm}( numAtomicMoves )
  for i = 1:numAtomicMoves
    j = atomicMovesValid[i]
    oneMoveAlgorithms[i] = Algorithm( [j], atomicCubeChanges[j] )
  end

  ################################################################################

end

max_stars_repo_name: clauderichard/RubikSearch

max_issues_repo_name: clauderichard/RubikSearch

max_forks_repo_name: clauderichard/RubikSearch


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_energy_share_requirement(setup::Dict, path::AbstractString, inputs_ESR::Dict)

Function for reading input parameters related to mimimum energy share requirement constraints (e.g. renewable portfolio standard or clean electricity standard policies)
"""
function load_energy_share_requirement(setup::Dict, path::AbstractString, inputs_ESR::Dict)
	# Definition of ESR requirements by zone (as % of load)
	# e.g. any policy requiring a min share of qualifying resources (Renewable Portfolio Standards / Renewable Energy Obligations / Clean Energy Standards etc.)
	inputs_ESR["dfESR"] = DataFrame(CSV.File(joinpath(path,"Energy_share_requirement.csv"), header=true), copycols=true)
	# Ensure float format values:
	ESR = count(s -> startswith(String(s), "ESR"), names(inputs_ESR["dfESR"]))
	first_col = findall(s -> s == "ESR_1", names(inputs_ESR["dfESR"]))[1]
	last_col = findall(s -> s == "ESR_$ESR", names(inputs_ESR["dfESR"]))[1]

	inputs_ESR["dfESR"] = Matrix{Float64}(inputs_ESR["dfESR"][:,first_col:last_col])
	inputs_ESR["nESR"] = ESR

	println("Energy_share_requirement.csv Successfully Read!")

	return inputs_ESR
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Test
using DisjunctiveConstraints
import Cbc, MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

include("sets.jl")
include("activity.jl")
include("big_m.jl")


text:
module Materials

    using ..FEPrep
    using StaticArrays, Parameters, LinearAlgebra

    include("material_data_types.jl")
    include("get_material_tensor.jl")
    include("material_lib.jl")

    export AbstractMaterial, Material
    export IsotropicMaterial, OrthotropicMaterial, TransverselyIsotropicMaterial

end


text:
module ParallelKMeans

using StatsBase
import Base.Threads: @spawn, @threads

export kmeans

"""
TODO 1: Document function
"""
function divider(n, k)
    d = div(n, k)
    xz = vcat(collect((0:k-1) * d), n)
    return [t[1]:t[2] for t in zip(xz[1:end-1] .+ 1, xz[2:end])]
end

"""
TODO 2: Document function
"""
function pl_pairwise!(target, x, y, nth = Threads.nthreads())
    ncol = size(x, 2)
    nrow = size(x, 1)
    ranges = divider(nrow, nth)
    waiting_list = Task[]
    for i in 1:length(ranges) - 1
        push!(waiting_list, @spawn inner_pairwise!(target, x, y, ranges[i]))
    end
    inner_pairwise!(target, x, y, ranges[end])
    for i in 1:length(ranges) - 1
        wait(waiting_list[i])
    end

    target
end

"""
TODO 3: Document function
"""
function inner_pairwise!(target, x, y, r)
    ncol = size(x, 2)
    @inbounds for k in axes(y, 1)
        for i in r
            target[i, k] = (x[i, 1] - y[k, 1])^2
        end

        for j in 2:ncol
            for i in r
                target[i, k] += (x[i, j] - y[k, j])^2
            end
        end
    end
    target
end

"""
TODO 4: Document function
"""
function pairwise!(target, x, y)
    ncol = size(x, 2)
    @inbounds for k in axes(y, 1)
        for i in axes(x, 1)
            target[i, k] = (x[i, 1] - y[k, 1])^2
        end

        for j in 2:ncol
            for i in axes(x, 1)
                target[i, k] += (x[i, j] - y[k, j])^2
            end
        end
    end
    target
end

"""
    smart_init(X, k; init="k-means++")

    This function handles the random initialisation of the centroids from the
    design matrix (X) and desired groups (k) that a user supplies.

    `k-means++` algorithm is used by default with the normal random selection
    of centroids from X used if any other string is attempted.

    A tuple representing the centroids, number of rows, & columns respecitively
    is returned.
"""
function smart_init(X::Array{Float64, 2}, k::Int; init::String="k-means++")
    n_row, n_col = size(X)

    if init == "k-means++"

        # randonmly select the first centroid from the data (X)
        centroids = zeros(k, n_col)
        rand_idx = rand(1:n_row)
        centroids[1, :] .= X[rand_idx, :]
        distances = Array{Float64}(undef, n_row, 1)
        new_distances = Array{Float64}(undef, n_row, 1)

        # compute distances from the first centroid chosen to all the other data points
        first_centroid_matrix = convert(Matrix, centroids[1, :]')

        # flatten distances
        # distances = vec(pairwise(SqEuclidean(), X, first_centroid_matrix, dims = 1))
        pairwise!(distances, X, first_centroid_matrix)

        for i = 2:k
            # choose the next centroid, the probability for each data point to be chosen
            # is directly proportional to its squared distance from the nearest centroid
            r_idx = sample(1:n_row, ProbabilityWeights(vec(distances)))
            centroids[i, :] .= X[r_idx, :]

            # Ignore setting the last centroid to help the separation of centroids
            if i == (k-1)
                break
            end

            # compute distances from the centroids to all data points
            current_centroid_matrix = convert(Matrix, centroids[i, :]')
            # new_distances = vec(pairwise(SqEuclidean(), X, current_centroid_matrix, dims = 1))
            pairwise!(new_distances, X, first_centroid_matrix)

            # and update the squared distance as the minimum distance to all centroid
            # distances = minimum([distances, new_distances])
            for i in 1:n_row
                distances[i, 1] = distances[i, 1] < new_distances[i, 1] ? distances[i, 1] : new_distances[i, 1]
            end
        end

    else
        # randomly select points from the design matrix as the initial centroids
        rand_indices = rand(1:n_row, k)
        centroids = X[rand_indices, :]

    end

    return centroids, n_row, n_col
end

"""
    sum_of_squares(x, labels, centre, k)

    This function computes the total sum of squares based on the assigned (labels)
    design matrix(x), centroids (centre), and the number of desired groups (k).

    A Float type representing the computed metric is returned.
"""
function sum_of_squares(x::Array{Float64,2}, labels::Array{Int64,1}, centre::Array)
    s = 0.0

    @inbounds for j in axes(x, 2)
        for i in axes(x, 1)
            s += (x[i, j] - centre[labels[i], j])^2
        end
    end

    return s
end

"""
    Kmeans(design_matrix, k; k_init="k-means++", max_iters=300, tol=1e-4, verbose=true)

    This main function employs the K-means algorithm to cluster all examples
    in the training data (design_matrix) into k groups using either the
    `k-means++` or random initialisation technique for selecting the initial
    centroids.

    At the end of the number of iterations specified (max_iters), convergence is
    achieved if difference between the current and last cost objective is
    less than the tolerance level (tol). An error is thrown if convergence fails.

    Details of operations can be either printed or not by setting verbose accordingly.

    A tuple representing labels, centroids, and sum_squares respectively is returned.

"""
function kmeans(design_matrix::Array{Float64, 2}, k::Int; k_init::String = "k-means++",
    max_iters::Int = 300, tol = 1e-4, verbose::Bool = true)

    centroids, n_row, n_col = smart_init(design_matrix, k, init=k_init)

    labels = Vector{Int}(undef, n_row)
    distances = Vector{Float64}(undef, n_row)
    centroids_cnt = Vector{Int}(undef, size(centroids, 1))

    J_previous = Inf64

    nearest_neighbour = Array{Float64, 2}(undef, size(design_matrix, 1), size(centroids, 1))
    # Update centroids & labels with closest members until convergence
    for iter = 1:max_iters
        pairwise!(nearest_neighbour, design_matrix, centroids)

        @inbounds for i in axes(nearest_neighbour, 1)
            labels[i] = 1
            distances[i] = nearest_neighbour[i, 1]
            for j in 2:size(nearest_neighbour, 2)
                if distances[i] > nearest_neighbour[i, j]
                    labels[i] = j
                    distances[i] = nearest_neighbour[i, j]
                end
            end
        end

        centroids .= 0.0
        centroids_cnt .= 0
        @inbounds for i in axes(design_matrix, 1)
            centroids[labels[i], 1] += design_matrix[i, 1]
            centroids_cnt[labels[i]] += 1
        end
        @inbounds for j in 2:n_col
            for i in axes(design_matrix, 1)
                centroids[labels[i], j] += design_matrix[i, j]
            end
        end
        centroids ./= centroids_cnt

        # Cost objective
        J = mean(distances)

        if verbose
            # Show progress and terminate if J stopped decreasing.
            println("Iteration ", iter, ": Jclust = ", J, ".")
        end

        # Final Step: Check for convergence
        if iter > 1 && abs(J - J_previous) < (tol * J)

            sum_squares = sum_of_squares(design_matrix, labels, centroids)

            # Terminate algorithm with the assumption that K-means has converged
            if verbose
                println("Successfully terminated with convergence.")
            end

            return labels, centroids, sum_squares

        elseif iter == max_iters && abs(J - J_previous) > (tol * J)
            throw(error("Failed to converge Check data and/or implementation or increase max_iter."))

        end

        J_previous = J
    end

end

end # module


text:
using Test
using SerialDependence
using DelimitedFiles

#testing the functions by reproducing the results from C. Weiss's book "An Introduction to Discrete-Valued Time Series".
#I could not find the raw data so I copied the values provided in the book (the raw time-series should be ~1300 elements long), explaining slight numerical differences in results.
pewee = Int64.(readdlm("pewee.txt", ',')[1,:])
@test round(cramer_coefficient(pewee, Int64(4))[1], digits = 2) == 0.46
@test round(cohen_coefficient(pewee, Int64(4))[1], digits = 2) == 0.55
# this one was not in the book, but the plot correctly reproduces the results.
@test round(theils_u(pewee,[1])[1], digits = 2) == 0.55
@test round(H(pewee), digits = 2) == 1.49