text:
module day11

import ..aoc19: run_program!
using DelimitedFiles, UnicodePlots, REPL

readInput() = parse.(Int, split(readline(joinpath(@__DIR__,"input.txt")), ','))

const UP = (1, 0)
const RIGHT = (0, 1)
const DOWN = (-1, 0)
const LEFT = (0, -1)
turnLeft(d) = d == UP ? LEFT : (d == RIGHT ? UP : (d == DOWN ? RIGHT : DOWN))
turnRight(d) = d == UP ? RIGHT : (d == RIGHT ? DOWN : (d == DOWN ? LEFT : UP))

function part1(data)
    pos = (0, 0)
    dir = UP
    painted = Dict{Tuple{Int, Int}, Int}()

    input = Channel{Int}(Inf)
    output = Channel{Int}(Inf)
    task = @async run_program!(data, input, output)

    while !istaskdone(task)
        put!(input, get(painted, pos, 0))
        color = take!(output)
        turn = take!(output)
        painted[pos] = color # paint
        dir = turn == 0 ? turnLeft(dir) : turnRight(dir) # turn
        pos = pos .+ dir # move forward
    end

    length(painted)
end

function part2(data)
    pos = (0, 0)
    dir = UP
    painted = Dict((0, 0) => 1)

    input = Channel{Int}(Inf)
    output = Channel{Int}(Inf)
    task = @async run_program!(data, input, output)

    while !istaskdone(task)
        put!(input, get(painted, pos, 0))
        color = take!(output)
        turn = take!(output)
        painted[pos] = color # paint
        dir = turn == 0 ? turnLeft(dir) : turnRight(dir) # turn
        pos = pos .+ dir # move forward
    end

    white = filter(p -> p.second == 1, painted)

    scatterplot(last.(keys(white)), first.(keys(white)), xlim=extrema(last.(keys(painted))),
        ylim=extrema(first.(keys(painted))), height = 5, width = 70, color = :blue, canvas = BlockCanvas)

end

end


text:
import Base: +, -
import Base.diff

(+)(ta::TimeArray) = .+ta
(-)(ta::TimeArray) = .-ta

###### lag, lead ################

function lag(ta::TimeArray{T, N}, n::Int=1;
             padding::Bool=false, period::Int=0) where {T, N}

    if period != 0
        @warn("the period kwarg is deprecated, use lag(ta::TimeArray, period::Int) instead")
        n = period
    end

    # TODO: apply `unchecked`
    if padding
        paddedvals = [NaN * ones(n,size(ta, 2)); values(ta)[1:end-n, :]]
        ta = TimeArray(timestamp(ta), paddedvals, colnames(ta), meta(ta))
    else
        ta = TimeArray(timestamp(ta)[1+n:end], values(ta)[1:end-n, :], colnames(ta), meta(ta))
    end

    N == 1 && (ta = ta[colnames(ta)[1]])
    return ta

end  # lag

function lead(ta::TimeArray{T, N}, n::Int=1;
              padding::Bool=false, period::Int=0) where {T, N}

    if period != 0
      @warn("the period kwarg is deprecated, use lead(ta::TimeArray, period::Int) instead")
      n = period
    end

    if padding
        paddedvals = [values(ta)[1+n:end, :]; NaN*ones(n, length(colnames(ta)))]
        ta = TimeArray(timestamp(ta), paddedvals, colnames(ta), meta(ta))
    else
        ta = TimeArray(timestamp(ta)[1:end-n], values(ta)[1+n:end, :], colnames(ta), meta(ta))
    end

    N == 1 && (ta = ta[colnames(ta)[1]])
    return ta

end  # lead

###### diff #####################

function diff(ta::TimeArray, n::Int=1; padding::Bool=false, differences::Int=1)
    cols = colnames(ta)
    for d in 1:differences
        ta = ta .- lag(ta, n, padding=padding)
    end
    colnames(ta)[:] = cols
    return ta
end  # diff

###### percentchange ############

function percentchange(ta::TimeArray, returns::Symbol=:simple;
                       padding::Bool=false, method::AbstractString="")

    if method != ""
        @warn("the method kwarg is deprecated, use percentchange(ta, :methodname) instead")
        returns = Symbol(method)
    end

    cols = colnames(ta)
    ta = returns == :log ? diff(log.(ta), padding=padding) :
         returns == :simple ? expm1.(percentchange(ta, :log, padding=padding)) :
         throw(ArgumentError("returns must be either :simple or :log"))
    colnames(ta)[:] = cols

   return ta

end  # percentchange

###### moving ###################

# Note: please do not involve any side effects in function `f`
"""
    moving(f, ta::TimeArray{T,1}, w::Integer; padding = false)

Apply user-defined function `f` to a 1D `TimeArray` with window size `w`.

## Example
To calculate the simple moving average of a time series:

```julia
moving(mean, ta, 10)
```
"""
function moving(f, ta::TimeArray{T,1}, window::Integer;
                padding::Bool = false) where {T}
    ts   = padding ? timestamp(ta) : @view(timestamp(ta)[window:end])
    A    = values(ta)
    vals = map(i -> f(@view A[i:i+(window-1)]), 1:length(ta)-window+1)
    padding && (vals = [fill(NaN, window - 1); vals])
    TimeArray(ta; timestamp = ts, values = vals)
end

"""
    moving(f, ta::TimeArray{T,2}, w::Integer; padding = false, dims = 1, colnames = [...])

## Example
In case of `dims = 2`, the user-defined function `f` will get a 2D `Array` as input.

```julia
moving(ohlc, 10, dims = 2, colnames = [:A, ...]) do
    # given that `ohlc` is a 500x4 `TimeArray`,
    # size(A) is (10, 4)
    ...
end
```
"""
function moving(f, ta::TimeArray{T,2}, window::Integer;
                padding::Bool = false, dims::Integer = 1,
                colnames::AbstractVector{Symbol} = _colnames(ta)) where {T}
    if !(dims ∈ (1, 2))
        throw(ArgumentError("invalid dims $dims"))
    end

    ts   = padding ? timestamp(ta) : @view(timestamp(ta)[window:end])
    A    = values(ta)

    if dims == 1
        vals = similar(@view(A[window:end, :]))
        for i ∈ 1:size(vals, 1), j ∈ 1:size(vals, 2)
            vals[i, j] = f(@view(A[i:i+(window-1), j]))
        end
    else # case of dims = 2
        vals = mapreduce(i -> f(view(A, i-window+1:i, :)), vcat, window:size(A, 1))
        if size(vals, 2) != length(colnames)
            throw(DimensionMismatch(
                "the output dims should match the lenght of columns, " *
                "please set the keyword argument `colnames` properly."
            ))
        end
    end

    padding && (vals = [fill(NaN, (window-1), size(vals, 2)); vals])
    TimeArray(ta; timestamp = ts, values = vals, colnames = colnames)
end

###### upto #####################

function upto(f, ta::TimeArray{T, 1}) where {T}
    vals = zero(values(ta))
    for i=1:length(vals)
        vals[i] = f(values(ta)[1:i])
    end
    TimeArray(timestamp(ta), vals, colnames(ta), meta(ta))
end

function upto(f, ta::TimeArray{T, 2}) where {T}
    vals = zero(values(ta))
    for i=1:size(vals, 1), j=1:size(vals, 2)
        vals[i, j] = f(values(ta)[1:i, j])
    end
    TimeArray(timestamp(ta), vals, colnames(ta), meta(ta))
end

###### basecall #################

basecall(ta::TimeArray, f::Function; cnames=colnames(ta)) =
    TimeArray(timestamp(ta), f(values(ta)), cnames, meta(ta))

###### uniform observations #####

function uniformspaced(ta::TimeArray)
    gap1 = timestamp(ta)[2] - timestamp(ta)[1]
    i, n, is_uniform = 2, length(ta), true
    while is_uniform & (i < n)
        is_uniform = gap1 == (timestamp(ta)[i+1] - timestamp(ta)[i])
        i += 1
    end
    return is_uniform
end  # uniformspaced

function uniformspace(ta::TimeArray{T, N}) where {T, N}
    min_gap = minimum(timestamp(ta)[2:end] - timestamp(ta)[1:end-1])
    newtimestamp = timestamp(ta)[1]:min_gap:timestamp(ta)[end]
    emptyta = TimeArray(collect(newtimestamp), zeros(length(newtimestamp), 0),
                        Symbol[], meta(ta))
    ta = merge(emptyta, ta, method = :left)
    N == 1 && (ta = ta[colnames(ta)[1]])
    return ta
end  # uniformspace

###### dropnan ####################

dropnan(ta::TimeArray, method::Symbol = :all) =
    method == :all ? ta[findall(reshape(values(any(.!isnan.(ta), dims = 2)), :))] :
    method == :any ? ta[findall(reshape(values(all(.!isnan.(ta), dims = 2)), :))] :
    throw(ArgumentError("dropnan method must be :all or :any"))

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule booster_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("booster")
JLLWrappers.@generate_main_file("booster", UUID("511df3fd-a8e5-5b39-b2c4-0599bfe02772"))
end  # module booster_jll


text:
# period doubling route
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

fig = figure(figsize = (figx, 2figy))
axlo = fig.add_subplot(3, 1, 1)
axro = fig.add_subplot(3, 1, 3)

zoomaxs = []
for i in 1:3
	push!(zoomaxs, fig.add_subplot(3, 3, 3+i))
end

axlo.clear()
ds = Systems.logistic(0.5)
i = 1
pvalues = range(2.5, 4; length = 2001)
n = 2000
Ttr = 2000
p_index = 1
output = orbitdiagram(ds, i, p_index, pvalues; n = n, Ttr = Ttr)

L = length(pvalues)
x = Vector{Float64}(undef, n*L)
y = copy(x)
for j in 1:L
    x[(1 + (j-1)*n):j*n] .= pvalues[j]
    y[(1 + (j-1)*n):j*n] .= output[j]
end
axlo.plot(x, y, ls = "None", ms = 0.5, color = "black", marker = "o", alpha = 0.01)
axlo.set_xlim(pvalues[1], pvalues[end]); axlo.set_ylim(0,1)
axlo.set_ylabel("\$x\$", labelpad = 10)
axlo.set_xticks(2.5:0.5:4)
axlo.set_yticks(0:0.2:1)
axlo.set_xlabel("\$r\$"; labelpad = -20)
axlo.text(2.6, 0.1, "logistic map", size = 36)
r∞ = 3.56995
axlo.axvline(r∞; color = "C4", alpha = 0.7, ls = "dashed")

# Zoomin boxes
zbox1 = ((3.629, 0.295), (3.640, 0.32))
zbox2 = ((3.6338, 0.304), (3.6342, 0.3065))
zbox3 = ((3.63407, 0.30495), (3.63408, 0.3052))
allboxes = (zbox1, zbox2, zbox3)
for (k, zbox) in enumerate(allboxes)
	pval = range(zbox[1][1], zbox[2][1]; length = 1001)
	nsm = n÷3
	output = orbitdiagram(ds, i, p_index, pval;
		n=nsm, Ttr, ulims = (zbox[1][2], zbox[2][2])
	)
	L = length(pval)
	x = Vector{Float64}(undef, nsm*L)
	y = copy(x)
	for j in 1:L
	    x[(1 + (j-1)*nsm):j*nsm] .= pval[j]
	    y[(1 + (j-1)*nsm):j*nsm] .= output[j]
	end
	zoomaxs[k].clear()
	zoomaxs[k].plot(x, y, ls = "None", ms = 0.5, color = "black", marker = "o", alpha = 0.01)
	zoomaxs[k].set_xlim(zbox[1][1], zbox[2][1])
	zoomaxs[k].set_ylim(zbox[1][2], zbox[2][2])
	zoomaxs[k].axis("off")
end

zbox0 = ((pvalues[1], 0.0), (pvalues[end], 1.0))
allzoomaxs = [axlo, zoomaxs...]
for i in 1:3
	origin = allzoomaxs[i]
	zoomin = allzoomaxs[i+1]
    zbox = allboxes[i]
	c = ("C1", "C3", "C5")[i]
    axis_zoomin!(zoomin, origin, zbox, zbox, c;	connect_lines = i > 1, lw = 4.0)
end

# arrow indicator
axlo.arrow(
	(zbox1[1][1] + zbox1[2][1])/2, 0.05, 0, zbox1[2][2] - 0.1;
	color = "C1", width = 0.005, length_includes_head = true, head_length = 0.1
)

# Roessler
ro = Systems.roessler()
pvalues = range(1.5, stop = 5.18, length = 2001)
i = 1
plane = (2, 0.0)
tf = 4000.0
p_index = 3
output = produce_orbitdiagram(ro, plane, i, p_index, pvalues;
                              tfinal = tf, Ttr = 2000.0)

for (j, p) in enumerate(pvalues)
    axro.plot(fill(p, length(output[j])), output[j], lw = 0,
    marker = "o", ms = 0.2, color = "black", alpha = 0.1)
end
axro.set_xlim(pvalues[1], pvalues[end]);
axro.set_ylim(minimum(minimum(o) for o in output), maximum(maximum(o) for o in output))
axro.set_ylabel("\$x\$ (at section)", labelpad = 25)
axro.set_xlabel("\$c\$"; labelpad = -20)
axro.set_yticks(3:3:9)
axro.text(1.7, 7, "Rössler system", size = 36)
fig.tight_layout(pad=0.3)
fig.subplots_adjust(hspace = 0.3, wspace = 0.1)
add_identifiers!(fig, (axlo, axro))
wsave(plotsdir("4", "orbit_diagrams"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
#= ------------------------------------------------------------------------

    Vector{Region}

------------------------------------------------------------------------ =#

export left, top, right, bottom, bounds
export regions_to_image

"""
    left(x::Vector{Region})

Calculates the leftmost coordinate of a vector of regions.

This function works for non-complement and non-empty regions only. The vector
must not be empty.
"""
function left(regions::Vector{Region})
    if length(regions) > 0
        l = left(regions[1])
        for region in regions
            l = min(l, left(region))
        end
        return l
    else
        return missing
    end
end

"""
    top(x::Vector{Region})

Calculates the topmost region coordinate of a vector of regions.

This function works for non-complement and non-empty regions only. The vector
must not be empty.
"""
function top(regions::Vector{Region})
    if length(regions) > 0
        t = top(regions[1])
        for region in regions
            t = max(t, top(region))
        end
        return t
    else
        return missing
    end
end

"""
    right(x::Vector{Region})

Calculates the rightmost region coordinate of a vector of regions.

This function works for non-complement and non-empty regions only. The vector
must not be empty.
"""
function right(regions::Vector{Region})
    if length(regions) > 0
        r = right(regions[1])
        for region in regions
            r = max(r, right(region))
        end
        return r
    else
        return missing
    end
end

"""
    bottom(x::Vector{Region})

Calculates the bottommost region coordinate of a vector of regions.

This function works for non-complement and non-empty regions only. The vector
must not be empty.
"""
function bottom(regions::Vector{Region})
    if length(regions) > 0
        b = bottom(regions[1])
        for region in regions
            b = min(b, bottom(region))
        end
        return b
    else
        return missing
    end
end

"""
    bounds(x::Vector{Region})

Calculates the left, top, right and bottom region coordinate of a vector of regions
and returns them as a tuple.

This function works for non-complement and non-empty regions only. The vector
must not be empty.
"""
function bounds(regions::Vector{Region})
    if length(regions) > 0
        l = left(regions[1])
        t = top(regions[1])
        r = right(regions[1])
        b = bottom(regions[1])
        for region in regions
            l = min(l, left(region))
            t = max(t, top(region))
            r = max(r, right(region))
            b = min(b, bottom(region))
        end
        return l, t, r, b
    else
        return missing
    end
end

"""
    regions_to_image(r::Region, color=Gray(true))

Converts a vector of regions to an image. The function determines the bounds of the
regions and then renders the regions into the image.

The background of the image is filled with zeroes, the region pixels are
colored with the passed in colors.

Some examples of colors that you can pass:
[Gray(0.5)] : mid gray value
[RGB(1, 0, 0), RGB(1, 0, 0), RGB(1, 0, 0)] : cycle through red, green and blue
[RGBA(0, 0.5, 0, 0.5), RGBA(0.5, 0, 0, 0.5)] : half transparent mid green and mid red values
"""
function regions_to_image(regions::Vector{Region}, colors=[Gray(true)])
    (l, t, r, b) = bounds(regions)
    img = zeros(eltype(colors), t-b+1, r-l+1)
    n = 1
    for region in regions
        for run in region.runs
            for column in run.columns
                ca = color(colors[n])
                aa = alpha(colors[n])
                cb = color(img[run.row-b+1, column-l+1])
                ab = alpha(img[run.row-b+1, column-l+1])
                ac = aa + (1 - aa) * ab
                cc = (aa * ca + (1 - aa) * ab * cb) / ac
                img[run.row-b+1, column-l+1] = eltype(colors)(cc, ac)
            end
        end
        n += 1
        if (n > length(colors))
            n = 1
        end
    end
    return img
end

"""
    Base.show(io, mime::MIME"image/png", regions::Region[])

Shows a rich graphical display of a vector of regions. The colors cycle through
half transparent blue, green, red, cyan, magenta and yellow colors.
"""
function Base.show(io::IO, mime::MIME"image/png", regions::Vector{Region})
    colors = RGBA[RGBA(0,0,1,0.5), RGBA(0,1,0,0.5), RGBA(1,0,0,0.5), RGBA(0,1,1,0.5), RGBA(1,0,1,0.5), RGBA(1,1,0,0.5)]
    Base.show(io, mime, regions_to_image(regions, colors))
end


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

name = "Snowball"
version = v"2.0.0"

# Collection of sources required to build SDL2
sources = [
    ArchiveSource("https://snowballstem.org/dist/libstemmer_c.tgz", "054e76f2a05478632f2185025bff0b98952a2b7aed7c4e0960d72ba565de5dfc"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libstemmer_c/
cp ../CMakeLists.txt .
rm Makefile
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
make -j${nproc}
make install
install_license COPYING
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct(["libstemmer", "stemmer"], :libstemmer)
]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
## spectral.jl : Gaussian random field generator based on the spectral (eigenvalue) decomposition of the covariance matrix

"""
Spectral <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) based on a spectral (eigenvalue) decomposition of the covariance matrix.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,Spectral(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a Spectral decomposition

julia> plot(grf)
[...]

```
This is also useful when computing Gaussian random fields on a Finite Element mesh using a truncated KL expansion. Here's an example that computes the first 10 eigenfunctions on an L-shaped domain.

```jldoctest
julia> p = readdlm(Pkg.dir("GaussianRandomFields")*"/data/Lshape.p");

julia> t = readdlm(Pkg.dir("GaussianRandomFields")*"/data/Lshape.t",Int64);

julia> grf = grf = GaussianRandomField(CovarianceFunction(2,Matern(0.2,1.0)),Spectral(),p,t,n=10)
Gaussian random field with 2d Matérn covariance function (λ=0.2, ν=1.0, σ=1.0, p=2.0) on a mesh with 998 points and 1861 elements, using a spectral decomposition

julia> tricontourf(p[:,1],p[:,2],grf.data.eigenfunc[:,1],triangles=t-1,cmap=get_cmap("viridis"))
[...]

```
See also: [`Cholesky`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref)
"""
struct Spectral <: GaussianRandomFieldGenerator end

const SpectralGRF = GaussianRandomField{C,Spectral} where {C}

# container type for eigenvalues and eigenfunctions
mutable struct SpectralData{X,Y}
	eigenval::X
	eigenfunc::Y
end

function _GaussianRandomField(mean,cov,method::Spectral,pts...;n::N=0) where {N<:Integer}
	C = apply(cov,pts,pts)

	# compute eigenvalue decomposition
	if n == 0
		F = eigen(C)
		idx = sortperm(F.values,rev=true)
		Λ = F.values[idx]
		U = F.vectors[:,idx]
	else # thanks to #4
		(eigenval,eigenfunc) = eigs(C,nev=n,ritzvec=true,which=:LM)
		idx = sortperm(eigenval,rev=true)
		Λ = eigenval[idx]
		U = eigenfunc[:,idx]
	end

	# if negative eigenvalues detected, remove them
	n = find_last_positive(Λ)
	n == length(Λ) || warn("negative eigenvalue $(Λ[n+1]) detected, Gaussian random field will be approximated (ignoring all negative eigenvalues)")

	# store eigenvalues and eigenfunctions
	data = SpectralData(sqrt.(Λ[1:n]),U[:,1:n]) # note: store sqrt of eigenval for more efficient sampling

	GaussianRandomField{typeof(cov),Spectral,typeof(pts)}(mean,cov,pts,data)
end

# find last positive entry in vector, assumes sorted from high to low
function find_last_positive(x::Vector{T} where {T})
	found = false
	n = 0
	while !found
		n += 1
		if ( n > length(x) ) || ( x[n] < 0 )
			found = true
		end
	end
	return n-1
end

# returns the required dimension of the random points
randdim(grf::SpectralGRF) = length(grf.data.eigenval)

# see KarhunenLoeve.jl for sample function

show(io::IO,::Spectral) = print(io,"spectral decomposition")


text:
# m5.4.MAs.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

# ### snippet 5.1

df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);

df = DataFrame(
    :A => df[:, :MedianAgeMarriage],
    :M => df[:, :Marriage],
    :D => df[:, :Divorce]
 )

scale!(df, [:M, :A, :D])

# Define the Stan language model

stan5_4_MA = "
data {
    int N;
    vector[N] A;
    vector[N] M;
}
parameters {
    real a;
    real bMA;
    real<lower=0> sigma;
}
model {
    vector[N] mu = a + bMA * A;
    a ~ normal( 0 , 0.2 );
    bMA ~ normal( 0 , 0.5 );
    sigma ~ exponential( 1 );
    M ~ normal( mu , sigma );
}
";

# Define the SampleModel
m5_4_MAs = SampleModel("m5.4", stan5_4_MA);

# Input data

m5_4_data = Dict(
    "N" => size(df, 1),
    "M" => df[:, :M_s],
    "A" => df[:, :A_s]
);

# Sample using cmdstan

rc5_4_MAs = stan_sample(m5_4_MAs, data=m5_4_data);

if success(rc5_4_MAs)

    # Rethinking results

    rethinking_results = "
           mean   sd  5.5% 94.5%
    a      0.00 0.09 -0.14  0.14
    bMA   -0.69 0.10 -0.85 -0.54
    sigma  0.68 0.07  0.57  0.79
    ";

    part5_4_MAs = read_samples(m5_4_MAs, :particles)
    part5_4_MAs |> display

end

# End of m5.4.MAs.jl


text:
using MNIST, ReverseDiff

const BATCH_SIZE = 100
const IMAGE_SIZE = 784
const CLASS_COUNT = 10

##################
# data wrangling #
##################

# loading MNIST data #
#--------------------#

function load_mnist(data)
    images, label_indices = data
    labels = zeros(CLASS_COUNT, length(label_indices))
    for i in eachindex(label_indices)
        labels[Int(label_indices[i]) + 1, i] = 1.0
    end
    return images, labels
end

const TRAIN_IMAGES, TRAIN_LABELS = load_mnist(MNIST.traindata())
const TEST_IMAGES, TEST_LABELS = load_mnist(MNIST.testdata())

# loading batches #
#-----------------#

immutable Batch{W,B,P,L}
    weights::W
    bias::B
    pixels::P
    labels::L
end

function Batch(images, labels, i)
    weights = zeros(CLASS_COUNT, IMAGE_SIZE)
    bias = zeros(CLASS_COUNT)
    range = i:(i + BATCH_SIZE - 1)
    return Batch(weights, bias, images[:, range], labels[:, range])
end

function load_batch!(batch, images, labels, i)
    offset = i - 1
    for batch_col in 1:BATCH_SIZE
        data_col = batch_col + offset
        for k in 1:size(images, 1)
            batch.pixels[k, batch_col] = images[k, data_col]
        end
        for k in 1:size(labels, 1)
            batch.labels[k, batch_col] = labels[k, data_col]
        end
    end
    return batch
end

####################
# model definition #
####################

# Here, we make a model out of simple `softmax` and `cross_entropy` functions. This could be
# improved by implementing something like Tensorflow's `softmax_cross_entropy_with_logits`,
# but my main goal is to show off ReverseDiff rather than implement the best possible model.

# Also note that our input's orientation is transposed compared to example implementations
# presented by row-major frameworks like Tensorflow. Julia is column-major, so I've set up
# the `Batch` code (see above) such that each column of `pixels` is an image and
# `size(pixels, 2) == BATCH_SIZE`.

# objective definitions #
#-----------------------#

# Here we use `@forward` to tell ReverseDiff to differentiate this scalar function in
# forward-mode. This allows us to call `minus_log.(y)` instead of `-(log.(y))`. By defining
# our own "fused" `minus_log` kernel using `@forward`, the operation `minus_log.(y)` becomes
# a single array instruction in the tape (instead of two separate ones,  as is the case with
# `-(log.(y))`), buying us a slight performance gain. Note that we didn't *need* to do this;
# it's simply a good place to show off the `@forward` feature.
ReverseDiff.@forward minus_log(x::Real) = -log(x)

cross_entropy(y′, y) = mean(sum(y′ .* (minus_log.(y)), 1))

softmax(x) = (exp_x = exp.(x); exp_x ./ sum(exp_x, 1))

function model(weights, bias, pixels, labels)
    y = softmax((weights * pixels) .+ bias)
    return cross_entropy(labels, y)
end

# gradient definitions #
#----------------------#

# generate `∇model!(output, input)` where the output/input takes the same form as `seeds`
const ∇model! = begin
    batch = Batch(TRAIN_IMAGES, TRAIN_LABELS, 1)
    seeds = (batch.weights, batch.bias, batch.pixels, batch.labels)
    ReverseDiff.compile_gradient(model, seeds)
end

# Add convenience method to `∇model!` that translates `Batch` args to `Tuple` args. Since
# `∇model!` is a binding to a function type generated by ReverseDiff, we can overload
# the type's calling behavior just like any other Julia type. For details, see
# http://docs.julialang.org/en/release-0.5/manual/methods/#function-like-objects.
function (::typeof(∇model!))(output::Batch, input::Batch)
    output_tuple = (output.weights, output.bias, output.pixels, output.labels)
    input_tuple = (input.weights, input.bias, input.pixels, input.labels)
    return ∇model!(output_tuple, input_tuple)
end

############
# training #
############

function train_step!(∇batch::Batch, batch::Batch, rate = 0.5, iters = 20)
    for _ in 1:iters
        ∇model!(∇batch, batch)
        for i in eachindex(batch.weights)
            batch.weights[i] -= rate * ∇batch.weights[i]
        end
        for i in eachindex(batch.bias)
            batch.bias[i] -= rate * ∇batch.bias[i]
        end
    end
end

function train!(∇batch::Batch, batch::Batch, images, labels, rate = 0.5, iters = 20)
    batch_count = floor(Int, size(images, 2) / BATCH_SIZE)
    @time for i in 1:batch_count
        load_batch!(batch, images, labels, i)
        train_step!(∇batch, batch, rate, iters)
    end
    return ∇batch
end


text:
export ZLat

@attributes mutable struct ZLat <: AbsLat{FlintRationalField}
  space::QuadSpace{FlintRationalField, fmpq_mat}
  rational_span::QuadSpace{FlintRationalField, fmpq_mat}
  basis_matrix::fmpq_mat
  gram_matrix::fmpq_mat
  aut_grp_gen::fmpq_mat
  aut_grp_ord::fmpz
  automorphism_group_generators::Vector{fmpz_mat} # With respect to the
                                                  # basis of the lattice
  automorphism_group_order::fmpz
  minimum::fmpq

  scale::fmpq
  norm::fmpq

  function ZLat(V::QuadSpace{FlintRationalField, fmpq_mat}, B::fmpq_mat)
    z = new()
    z.space = V
    z.basis_matrix = B
    return z
  end
end

@attributes mutable struct QuadLat{S, T, U} <: AbsLat{S}
  space::QuadSpace{S, T}
  pmat::U
  gram::T                        # gram matrix of the matrix part of pmat
  rational_span::QuadSpace{S, T}
  base_algebra::S
  automorphism_group_generators::Vector{T}
  automorphism_group_order::fmpz
  generators
  minimal_generators
  norm
  scale

  function QuadLat{S, T, U}() where {S, T, U}
    return new{S, T, U}()
  end

  function QuadLat(K::S, G::T, P::U) where {S, T, U}
    space = quadratic_space(K, G)
    z = new{S, T, U}(space, P)
    z.base_algebra = K
    return z
  end

  function QuadLat(K::S, G::T) where {S, T}
    n = nrows(G)
    M = pseudo_matrix(identity_matrix(K, n))
    return QuadLat(K, G, M)
  end
end


text:
@symbol_func function K_omega(cur_reactor::AbstractReactor)
  cur_K = K_RP(cur_reactor)

  cur_K /= K_n(cur_reactor) ^ 2

  cur_K *= 1e6

  cur_K
end


text:
##
using MRIgeneralizedBloch
using DifferentialEquations
using StaticArrays
using Test

## choose random parameters
α = rand() * π
TRF = 100e-6 + 900e-6 * rand()
ω1 = α / TRF
B1 = 0.7 + 0.6 * rand()
ω0 = 1000 * rand() - 500
m0s = 0.2 * rand()
m0f = 1 - m0s
R1f = 0.5 + 0.2 * rand()
R1s = 2 + rand()
R2f = 1 / (40e-3 + 100e-3 * rand())
T2s = 5e-6 + 10e-6 * rand()
R2s = 1 / T2s
Rx = 100 * rand()

ms = m0s * rand()
mf = (1 - m0s) * rand()
ϑ = rand() * π / 2
φ = rand() * 2π

## Solve generalized Bloch-McConnell with Lorentzian lineshape
mfun(p, t; idxs = nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5)

p = (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, greens_lorentzian)
alg = MethodOfSteps(DP8())
m0 = [mf * sin(ϑ) * cos(φ), mf * sin(ϑ) * sin(φ), mf * cos(ϑ), ms, 1]

sol = solve(
    DDEProblem(MRIgeneralizedBloch.apply_hamiltonian_gbloch!, m0, mfun, (0.0, TRF), p),
    alg,
)

u_gBloch = sol[end]

## Solve original Bloch-McConnell with Lorentzian lineshape
H = @SMatrix [
    -R2f  -ω0       B1*ω1      0           0       0
      ω0 -R2f           0      0           0       0
  -B1*ω1    0 -R1f-Rx*m0s      0      Rx*m0f R1f*m0f
       0    0           0   -R2s       B1*ω1       0
       0    0      Rx*m0s -B1*ω1 -R1s-Rx*m0f R1s*m0s
       0    0           0      0           0       0
]

m0 = [mf * sin(ϑ) * cos(φ), mf * sin(ϑ) * sin(φ), mf * cos(ϑ), 0, ms, 1]
u_Bloch = exp(H * TRF) * m0

## Test!
max_error = 1e-3
@test u_gBloch[1] ≈ u_Bloch[1] atol = max_error
@test u_gBloch[2] ≈ u_Bloch[2] atol = max_error
@test u_gBloch[3] ≈ u_Bloch[3] atol = max_error
@test u_gBloch[4] ≈ u_Bloch[5] atol = max_error
@test u_gBloch[5] ≈ u_Bloch[6] atol = max_error


text:
#=
generate_data.jl

generates data based on the generative process of the following models:

    - CRP
    - CRF
    - RCRP
    - RCRF

TODO:
    use multiple dispatch and package all these function under one name

    RCRP: aa could be a vector. Each epoch different concentraion parameter
    RCRF: aa and gg could be vectors

11/02/2016
Adham Beyki, odinay@gmail.com
=#

function gen_CRP_data(n_customers::Int, aa::Float64)

    KK = 1

    zz  = zeros(Int, n_customers)
    nn = zeros(Int, KK)

    ii = 1
    kk = 1

    zz[ii]  = kk
    nn[kk] += 1

    for ii = 2:n_customers
        pp = zeros(Float64, KK+1)
        for kk = 1:KK
            pp[kk] = nn[kk] / (aa + ii-1)
        end
        pp[KK+1] = aa / (aa + ii-1)

        kk = sample(pp)

        if kk == KK+1
            push!(nn, 0)
            KK += 1
        end
        nn[kk] += 1
        zz[ii] = kk
    end

    nn, zz, KK
end

function gen_RCRP_data(N_t::Vector{Int}, aa::Float64)

    KK = 1
    TT = length(N_t)
    nn = zeros(Int, TT, KK)
    zz = Array(Vector{Int}, TT)

    tt = 1
    zz[tt] = zeros(Int, N_t[tt])

    ii = 1
    kk = 1
    zz[tt][ii] = kk
    nn[tt, kk] += 1

    for ii = 2:N_t[tt]
        pp = zeros(Float64, KK+1)
        for kk = 1:KK
            pp[kk] = nn[tt, kk] / (aa + ii-1)
        end
        pp[KK+1] = aa / (aa + ii-1)

        kk = sample(pp)

        if kk == KK+1
            nn = BIAS.add_column(nn)
            KK += 1
        end

        zz[tt][ii] = kk
        nn[tt, kk] += 1
    end

    for tt = 2:TT
        zz[tt] = zeros(Int, N_t[tt])

        for ii = 1:N_t[tt]
            pp = zeros(Float64, KK+1)
            for kk = 1:KK
                pp[kk] = (nn[tt-1, kk] + nn[tt, kk]) / (N_t[tt-1]+ii-1+aa)
            end
            pp[KK+1] = aa / (N_t[tt-1]+ii-1+aa)

            kk = sample(pp)

            if kk == KK+1
                nn = BIAS.add_column(nn)
                KK += 1
            end

            zz[tt][ii] = kk
            nn[tt, kk] += 1
        end
    end

    return nn, zz, KK
end

function gen_CRF_data(n_group_j::Vector{Int}, gg::Float64, aa::Float64, join_tables::Bool)

    n_groups = length(n_group_j)
    zz       = Array(Vector{Int}, n_groups)
    tji      = Array(Vector{Int}, n_groups)
    njt      = Array(Vector{Int}, n_groups)
    kjt      = Array(Vector{Int}, n_groups)

    KK = 1
    mm = zeros(Int, KK)
    nn = zeros(Int, n_groups, KK)

    jj      = 1
    n_tbls  = 1
    zz[jj]  = zeros(Int, n_group_j[jj])
    tji[jj] = zeros(Int, n_group_j[jj])
    njt[jj] = zeros(Int, n_tbls)
    kjt[jj] = zeros(Int, n_tbls)

    ii  = 1
    tbl = 1
    kk  = 1
    kjt[jj][tbl] = kk
    tji[jj][ii]  = tbl
    njt[jj][tbl] += 1
    zz[jj][ii]   = kjt[jj][tbl]
    nn[jj, kk]   += 1
    mm[kk]       += 1

    for ii = 2:n_group_j[jj]
        pp = zeros(Float64, n_tbls+1)
        for tbl=1:n_tbls
            pp[tbl] = njt[jj][tbl] / (aa + ii-1)
        end
        pp[n_tbls+1] = aa / (aa + ii-1)

        tbl = sample(pp)
        if tbl == n_tbls+1
            n_tbls += 1
            push!(njt[jj], 0)
            push!(kjt[jj], 0)

            pp = zeros(Float64, KK+1)
            for kk = 1:KK
                pp[kk] = mm[kk] / (gg + sum(mm))
            end
            pp[KK+1] = gg / (gg + sum(mm))

            kk = sample(pp)

            if kk == KK+1
                push!(mm, 0)
                nn = add_column(nn)
                KK += 1
            end
            kjt[jj][tbl] = kk
            mm[kk] += 1
        end

        tji[jj][ii] = tbl
        njt[jj][tbl] += 1
        kk = kjt[jj][tbl]
        zz[jj][ii] = kk
        nn[jj, kk] += 1

    end

    for jj = 2:n_groups

        # no one is in yet. There is only one table.
        n_tbls = 1

        zz[jj]  = zeros(Int, n_group_j[jj])
        tji[jj] = zeros(Int, n_group_j[jj])
        njt[jj] = zeros(Int, n_tbls)
        kjt[jj] = zeros(Int, n_tbls)

        # customer #1 sits at table #1
        ii = 1
        tbl = 1
        tji[jj][ii] = tbl
        njt[jj][tbl] += 1

        # sample the dish for table #1
        pp = zeros(Float64, KK+1)
        for kk = 1:KK
            pp[kk] = mm[kk] / (gg + sum(mm))
        end
        pp[KK+1] = gg / (gg + sum(mm))

        kk = sample(pp)
        if kk == KK+1
            push!(mm, 0)
            nn = add_column(nn)
            KK += 1
        end

        tji[jj][ii]  = tbl
        kjt[jj][tbl] = kk
        mm[kk] += 1
        zz[jj][ii] = kk
        nn[jj, kk] += 1

        for ii = 2:n_group_j[jj]

            # sample table
            pp = zeros(Float64, n_tbls+1)
            for tbl=1:n_tbls
                pp[tbl] = njt[jj][tbl] / (aa + ii-1)
            end
            pp[n_tbls+1] = aa / (aa + ii-1)

            tbl = sample(pp)
            if tbl == n_tbls+1       # if new table is instantiated

                n_tbls += 1
                push!(njt[jj], 0)
                push!(kjt[jj], 0)

                # sample the dish for the newly initiated table
                pp = zeros(Float64, KK+1)
                for kk = 1:KK
                    pp[kk] = mm[kk] / (gg + sum(mm))
                end
                pp[KK+1] = gg / (gg + sum(mm))

                kk = sample(pp)

                # if it's a new dish
                if kk == KK+1
                    push!(mm, 0)
                    nn = add_column(nn)
                    KK += 1
                end
                kjt[jj][tbl] = kk
                mm[kk] += 1
            end

            tji[jj][ii] = tbl
            njt[jj][tbl] += 1
            kk = kjt[jj][tbl]
            zz[jj][ii] = kk
            nn[jj, kk] += 1
        end
    end

    # joining the tables that serve the same dish in each restaurant
    if join_tables
        for jj = 1:n_groups
            kk_list = unique(kjt[jj])
            new_kjt = zeros(Int, length(kk_list))
            new_njt = zeros(Int, length(kk_list))

            ll = 1
            for unq_kk in kk_list
                idx = find(x -> x==unq_kk, kjt[jj])
                new_kjt[ll] = unq_kk
                new_njt[ll] = sum(njt[jj][idx])
                ll += 1
            end

            kjt[jj] = new_kjt
            njt[jj] = new_njt
        end

        mm = zeros(Int, KK)
        for jj = 1:n_groups
            for kk in kjt[jj]
                mm[kk] += 1
            end
        end
    end

    return(tji, njt, kjt, nn, mm, zz, KK)
end

function gen_RCRF_data(n_group_j::Vector{Vector{Int}}, gg::Float64, aa::Float64, join_tables::Bool=false)

    TT       = length(n_group_j)                        # number of epochs
    n_groups = [length(foo) for foo in n_group_j]       # number of groups per epoch

    # initial KK
    KK = 1

    zz  = Array(Vector{Vector{Int}}, TT)                # zz[tt][jj][ii] is the cluster customer ii at restaurant jj in epoch tt belongs to
    njt = Array(Vector{Vector{Int}}, TT)                # njt[tt][jj][tbl] is the number of customers sitting at table tbl at restaurant jj in epoch tt
    kjt = Array(Vector{Vector{Int}}, TT)                # kjt[tt][jj][tbl] is the dish being served at table tbl at restaurant jj in epoch tt
    mm  = zeros(Int, TT, KK)                            # mm[tt, kk] is the number of times dish kk served at a table in epoch tt
    nn  = Array(Matrix{Int}, TT)                        # nn[tt][jj, kk] is the number of observations that belong to cluster kk in restaurant jj at epoch tt

    # we start with epoch 1, restaurant 1, customer 1
    # Initializing epoch 1
    tt      = 1
    zz[tt]  = Array(Vector{Int}, n_groups[tt])
    njt[tt] = Array(Vector{Int}, n_groups[tt])
    kjt[tt] = Array(Vector{Int}, n_groups[tt])
    nn[tt]  = zeros(Int, n_groups[tt], KK)

    # Initializing restaurant 1
    jj          = 1
    n_tbls      = 1
    zz[tt][jj]  = zeros(Int, n_group_j[tt][jj])
    njt[tt][jj] = zeros(Int, n_tbls)
    kjt[tt][jj] = zeros(Int, n_tbls)

    # customer 1 enters restaurant 1 at epoch 1, picks the the first table
    # and orders the first dish
    ii  = 1
    kk  = 1
    tbl = 1

    njt[tt][jj][tbl] += 1
    kjt[tt][jj][tbl] = kk
    zz[tt][jj][ii]   = kk
    mm[tt, kk]       += 1
    nn[tt][jj, kk]   += 1

    # table assignment for customers 2:end in restaurant 1 epoch 1
    for ii = 2:n_group_j[tt][jj]

        pp = zeros(Float64, n_tbls+1)
        for tbl = 1:n_tbls
            pp[tbl] = njt[tt][jj][tbl] / (aa + ii-1)
        end
        pp[n_tbls+1] = aa / (aa + ii-1)

        tbl = sample(pp)

        # a new table is instantiated
        if tbl == n_tbls+1

            n_tbls += 1
            push!(njt[tt][jj], 0)
            push!(kjt[tt][jj], 0)

            # sample the dish for the newly instantiated table
            pp = zeros(Float64, KK+1)
            for kk = 1:KK
                pp[kk] = mm[tt, kk] / (gg + sum(mm[tt, :]))
            end
            pp[KK+1] = gg / (gg + sum(mm[tt, :]))

            kk = sample(pp)

            # if it is a new dish
            if kk == KK+1
                mm = add_column(mm)
                nn[tt] = add_column(nn[tt])
                KK += 1
            end
            kjt[tt][jj][tbl] = kk
            mm[tt, kk] += 1
        end

        njt[tt][jj][tbl] += 1
        kk = kjt[tt][jj][tbl]
        zz[tt][jj][ii] = kk
        nn[tt][jj, kk] += 1

    end

    for jj = 2:n_groups[tt]

        # no one is in yet
        n_tbls = 1

        zz[tt][jj]  = zeros(Int, n_group_j[tt][jj])
        njt[tt][jj] = zeros(Int, n_tbls)
        kjt[tt][jj] = zeros(Int, n_tbls)

        # customer 1 sits at the first table
        ii = 1
        tbl = 1
        njt[tt][jj][tbl] += 1

        # sample the dish for table 1
        pp = zeros(Float64, KK+1)
        for kk = 1:KK
            pp[kk] = mm[tt, kk] / (gg + sum(mm[tt, :]))
        end
        pp[KK+1] = gg / (gg + sum(mm[tt, :]))

        kk = sample(pp)
        if kk == KK+1
            mm = add_column(mm)
            nn[tt] = add_column(nn[tt])
            KK += 1
        end

        kjt[tt][jj][tbl] = kk
        zz[tt][jj][ii]   = kk
        mm[tt, kk]       += 1
        nn[tt][jj, kk]   += 1

        for ii = 2:n_group_j[tt][jj]

            # sample the table for customer ii
            pp = zeros(Float64, n_tbls+1)
            for tbl = 1:n_tbls
                pp[tbl] = njt[tt][jj][tbl] / (aa + ii-1)
            end
            pp[n_tbls+1] = aa / (aa + ii-1)

            tbl = sample(pp)

            # a new table is instantiated
            if tbl == n_tbls + 1

                n_tbls += 1
                push!(njt[tt][jj], 0)
                push!(kjt[tt][jj], 0)

                # sample the dish for the newly instantiated table
                pp = zeros(Float64, KK+1)
                for kk = 1:KK
                    pp[kk] = mm[tt, kk] / (gg + sum(mm[tt, :]))
                end
                pp[KK+1] = gg / (gg + sum(mm[tt, :]))

                kk = sample(pp)

                # if it is a new dish
                if kk == KK+1
                    mm = add_column(mm)
                    nn[tt] = add_column(nn[tt])
                    KK += 1
                end
                kjt[tt][jj][tbl] = kk
                mm[tt, kk] += 1
            end

            njt[tt][jj][tbl] += 1
            kk = kjt[tt][jj][tbl]
            zz[tt][jj][ii] = kk
            nn[tt][jj, kk] += 1

        end
    end

    # epoch 2:TT
    for tt = 2:TT

        zz[tt]  = Array(Vector{Int}, n_groups[tt])
        njt[tt] = Array(Vector{Int}, n_groups[tt])
        kjt[tt] = Array(Vector{Int}, n_groups[tt])
        nn[tt]  = zeros(Int, n_groups[tt], KK)

        for jj = 1:n_groups[tt]

            n_tbls = 1

            zz[tt][jj]  = zeros(Int, n_group_j[tt][jj])
            njt[tt][jj] = zeros(Int, n_tbls)
            kjt[tt][jj] = zeros(Int, n_tbls)

            # customer 1 sits at the first table
            ii = 1
            tbl = 1
            njt[tt][jj][tbl] += 1

            # sample the dish for table 1
            pp = zeros(Float64, KK+1)
            for kk = 1:KK
                pp[kk] = (mm[tt-1, kk] + mm[tt, kk]) / (gg + sum(mm[tt-1:tt, :]))
            end
            pp[KK+1] = gg / (gg + sum(mm[tt-1:tt, :]))

            kk = sample(pp)
            if kk == KK+1
                mm = add_column(mm)
                nn[tt] = add_column(nn[tt])
                KK += 1
            end

            kjt[tt][jj][tbl] = kk
            zz[tt][jj][ii]   = kk
            mm[tt, kk]       += 1
            nn[tt][jj, kk]   += 1

            for ii = 2:n_group_j[tt][jj]

                # sample the table for customer ii
                pp = zeros(Float64, n_tbls+1)
                for tbl = 1:n_tbls
                    pp[tbl] = njt[tt][jj][tbl] / (aa + ii-1)
                end
                pp[n_tbls+1] = aa / (aa + ii-1)

                tbl = sample(pp)

                # a new table is instantiated
                if tbl == n_tbls + 1

                    n_tbls += 1
                    push!(njt[tt][jj], 0)
                    push!(kjt[tt][jj], 0)

                    # sample the dish for the newly instantiated table
                    pp = zeros(Float64, KK+1)
                    for kk = 1:KK
                        pp[kk] = (mm[tt-1, kk] + mm[tt, kk]) / (gg + sum(mm[tt-1:tt, :]))
                    end
                    pp[KK+1] = gg / (gg + sum(mm[tt-1:tt, :]))

                    kk = sample(pp)

                    # if it is a new dish
                    if kk == KK+1
                        mm = add_column(mm)
                        nn[tt] = add_column(nn[tt])
                        KK += 1
                    end
                    kjt[tt][jj][tbl] = kk
                    mm[tt, kk] += 1
                end

                njt[tt][jj][tbl] += 1
                kk = kjt[tt][jj][tbl]
                zz[tt][jj][ii] = kk
                nn[tt][jj, kk] += 1

            end
        end
    end

    # joining the tables that serve the same dish in each restaurant
    if join_tables

        for tt = 1:TT
            for jj = 1:n_groups[tt]
                kk_list = unique(kjt[tt][jj])
                new_kjt = zeros(Int, length(kk_list))
                new_njt = zeros(Int, length(kk_list))

                ll = 1
                for unq_kk in kk_list
                    idx = find(x -> x==unq_kk, kjt[tt][jj])
                    new_kjt[ll] = unq_kk
                    new_njt[ll] = sum(njt[tt][jj][idx])
                    ll += 1
                end

                kjt[tt][jj] = new_kjt
                njt[tt][jj] = new_njt
            end
        end

        mm = zeros(Int, TT, KK)
        for tt = 1:TT
            for jj = 1:n_groups[tt]
                for kk in kjt[tt][jj]
                    mm[tt, kk] += 1
                end
            end
        end

    end # join_tables

    return(njt, kjt, nn, mm, zz, KK)
end # gen_rcrf_data


text:
"""
$(TYPEDSIGNATURES)

```jldoctest
julia> fold!(PlackettBurman(2))
8×3 DataFrame
 Row │ factor1  factor2  dummy1
     │ Int64    Int64    Int64
─────┼──────────────────────────
   1 │       1        1       1
   2 │       1       -1      -1
   3 │      -1       -1       1
   4 │      -1        1      -1
   5 │      -1       -1      -1
   6 │      -1        1       1
   7 │       1        1      -1
   8 │       1       -1       1

```
"""
function fold!(design::PlackettBurman)
    fold_dataframe!(design.matrix)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> fold!(FractionalFactorial2Level(@formula(y ~ a + b + a&b + c + a&c+ b&c + a&b&c)))
16×7 DataFrame
 Row │ a      b      c      a_b    a_c    b_c    a_b_c
     │ Int64  Int64  Int64  Int64  Int64  Int64  Int64
─────┼─────────────────────────────────────────────────
   1 │    -1     -1     -1      1      1      1     -1
   2 │     1     -1     -1     -1     -1      1      1
   3 │    -1      1     -1     -1      1     -1      1
   4 │     1      1     -1      1     -1     -1     -1
   5 │    -1     -1      1      1     -1     -1      1
   6 │     1     -1      1     -1      1     -1     -1
   7 │    -1      1      1     -1     -1      1     -1
   8 │     1      1      1      1      1      1      1
   9 │     1      1      1     -1     -1     -1      1
  10 │    -1      1      1      1      1     -1     -1
  11 │     1     -1      1      1     -1      1     -1
  12 │    -1     -1      1     -1      1      1      1
  13 │     1      1     -1     -1      1      1     -1
  14 │    -1      1     -1      1     -1      1      1
  15 │     1     -1     -1      1      1     -1      1
  16 │    -1     -1     -1     -1     -1     -1     -1

```
"""
function fold!(design::FractionalFactorial2Level)
    fold_dataframe!(design.matrix)
end

function fold_dataframe!(df::DataFrames.DataFrame)
    mirror_matrix = deepcopy(df)
    if mirror_matrix[1,1] isa Symbol
        for col in names(mirror_matrix)
            replace!(mirror_matrix[!, col], :low => :high, :high => :low)
        end
    else
        for col in names(mirror_matrix)
            replace!(mirror_matrix[!, col], 1 => -1, -1 => 1)
        end
    end
    append!(df, mirror_matrix)
end


text:
module TwoDims
    # It may be best to just to move the script here
    include("../stable_fluids_fft.jl")

    precompile(main, ())
end

max_stars_repo_name: Ceyron/stable-fluids-julia

max_issues_repo_name: Ceyron/stable-fluids-julia

max_forks_repo_name: Ceyron/stable-fluids-julia


text:
module StrongModuleTrees

export StrongModuleTree, strong_modules,
    leaves, leaf_count, first_leaf, last_leaf, rand_leaf,
    nodes, node_count, parent_node

using ..GraphModularDecomposition

## core type definition ##

struct StrongModuleTree{T} <: AbstractVector{T}
    kind::Symbol
    edge::Tuple
    nodes::Vector{Union{T,StrongModuleTree{T}}}
end

## constructing StrongModuleTrees ##

function StrongModuleTree(
    G::AbstractMatrix, # the graph/digraph/2-structure as a matrix
    p::Vector{Int} = graph_factorizing_permutation(G),
)
    # initialize data structures
    n = length(p)
    op = zeros(Int,n); op[1] = 1
    cl = zeros(Int,n); cl[n] = 1
    lc = collect(1:n-1)
    uc = collect(2:n)

    # count open and close parens in fracture tree
    # find lower and upper cutters for node pairs
    for j = 1:n-1
        for i = 1:j-1
            G[p[i],p[j]] == G[p[i],p[j+1]] &&
            G[p[j],p[i]] == G[p[j+1],p[i]] && continue
            op[i] += 1
            cl[j] += 1
            lc[j] = i
            break
        end
        j += 1
        for i = n:-1:j+1
            G[p[i],p[j-1]] == G[p[i],p[j]] &&
            G[p[j-1],p[i]] == G[p[j],p[i]] && continue
            op[j] += 1
            cl[i] += 1
            uc[j-1] = i
            break
        end
    end

    # remove non-module "dummy" nodes
    let s = Int[]
        for j = 1:n
            for _ = 1:op[j]; push!(s, j); end
            for _ = 1:cl[j]
                i = pop!(s)
                if i < j
                    l = minimum(lc[k] for k = i:j-1)
                    u = maximum(uc[k] for k = i:j-1)
                    i <= l && u <= j && continue
                end
                op[i] -= 1
                cl[j] -= 1
            end
        end
    end

    # create nodes for consecutive twins
    let s = Int[], t = Int[]
        l = 1
        for k = 1:n
            for _ = 1:op[k]+1
                push!(s, k) # matching node stack
                push!(t, l) # matching twin stack
                l = k
            end
            for c = cl[k]:-1:0
                i = pop!(t)
                j = pop!(s)
                l = i # continue twin chain by default
                i < j || continue
                if i <= lc[j-1] < uc[j-1] <= k
                    # this node and prev are twins
                    if c > 0
                        # not last parens ∴ last twin
                        op[i] += 1
                        cl[k] += 1
                        l = k + 1
                    end
                else # this node and prev not twins
                    if i < j-1
                        op[i] += 1
                        cl[j-1] += 1
                    end
                    l = j # this node starts new chain
                end
            end
        end
    end

    # remove singleton "dummy" nodes
    let s = Int[]
        for j = 1:n
            for _ = 1:op[j]; push!(s, j); end
            i′ = 0
            for _ = 1:cl[j]
                i = pop!(s)
                if i == i′
                    op[i] -= 1
                    cl[j] -= 1
                end
                i′ = i
            end
        end
    end
    op[1] -= 1
    cl[n] -= 1

    # construct and normalize the tree
    return StrongModuleTree(G, p, op, cl)
end

function StrongModuleTree(
    G::AbstractMatrix,    # graph/digraph/2-structure
    v::AbstractVector{T}, # modular permutation of vertices
    op::Vector{Int},      # open parens per vertex
    cl::Vector{Int},      # close parens per vertex
) where {T <: Any}
    # continues from end of StrongModuleTree(G, p)

    function classify_nodes(t::Vector)
        n = length(t)
        counts = zeros(Int, n)
        x, y = first_leaf(t[1]), first_leaf(t[2])
        edge = (G[y,x], G[x,y])
        local a, b
        for i = 1:n, j = 1:n
            i == j && continue
            x, y = first_leaf(t[i]), first_leaf(t[j])
            a, b = G[y,x], G[x,y]
            if edge == (a, b)
                counts[i] += 1
            elseif edge == (b, a)
                counts[j] += 1
            else
                break
            end
        end
        sort!(counts)
        kind = a == b && all(c -> c == n-1, counts) ? :complete :
            all(d -> d == 2, diff(counts)) ? :linear : :prime
        edge[1] <= edge[2] || (edge = reverse(edge))
        kind == :prime && (edge = ())
        StrongModuleTree{T}(kind, edge, map(x->x isa Vector ? classify_nodes(x) : x, t))
    end

    function delete_weak_modules!(t::StrongModuleTree)
        i = 0
        while (i += 1) <= length(t)
            x = t[i]
            x isa StrongModuleTree || continue
            delete_weak_modules!(x)
            t.kind == x.kind != :prime && t.edge == x.edge || continue
            splice!(t.nodes, i, x.nodes)
            i += length(x)
        end
    end

    s = Any[[]]
    for (j, x) = enumerate(v)
        for _ = 1:op[j]
            t = []
            push!(s[end], t)
            push!(s, t)
        end
        push!(s[end], x)
        for _ = 1:cl[j]
            pop!(s)
        end
    end
    t = classify_nodes(s[end])
    delete_weak_modules!(t)
    return t
end

## displaying StrongModuleTrees ##

function edge_string(t::StrongModuleTree, post::String="")
    t.kind == :prime    ? "" :
    t.kind == :complete ? "$(t.edge[1])$post" :
                          "$(join(t.edge,"/"))$post"
end

kind_string(t::StrongModuleTree) = "$(edge_string(t,"-"))$(t.kind)"
Base.summary(t::StrongModuleTree) = "$(length(t))-node $(kind_string(t)) $(typeof(t))"

function Base.show(io::IO, t::StrongModuleTree)
    if get(io, :compact, false)
        print(io,
            edge_string(t,"-"), t.kind, " ",
            node_count(t), "-node (",
            leaf_count(t), "-leaf) module: ",
            first_leaf(t)
        )
    else
        parens = t.kind == :prime ? "{}" : t.kind == :linear ? "[]" : "()"
        print(io, parens[1])
        for (i, x) in enumerate(t)
            print(io, x)
            i < length(t) && print(io, " ")
        end
        print(io, parens[2])
    end
end

## querying properties of StrongModuleTrees ##

node_count(t::StrongModuleTree) = length(t)
node_count(x::Any) = 1

leaf_count(t::StrongModuleTree) = sum(leaf_count, t.nodes)
leaf_count(v::Vector) = sum(leaf_count, v)
leaf_count(x::Any) = 1

first_leaf(t::StrongModuleTree) = first_leaf(first(t.nodes))
first_leaf(v::Vector) = first_leaf(first(v))
first_leaf(x::Any) = x

last_leaf(t::StrongModuleTree) = last_leaf(last(t.nodes))
last_leaf(v::Vector) = last_leaf(last(v))
last_leaf(x::Any) = x

rand_leaf(t::StrongModuleTree) = rand_leaf(rand(t.nodes))
rand_leaf(v::Vector) = rand_leaf(last(v))
rand_leaf(x::Any) = x

function leaves(t::StrongModuleTree{T}) where T
    L = T[]
    for x in t.nodes
        append!(L, leaves(x))
    end
    return L
end
leaves(x::Any) = [x]

function nodes!(v::Vector{StrongModuleTree{T}}, t::StrongModuleTree{T}) where T
    for x in t.nodes
        x isa StrongModuleTree || continue
        nodes!(v, x)
    end
    push!(v, t)
    return v
end
nodes(t::StrongModuleTree) = nodes!(typeof(t)[], t)

strong_modules(t::StrongModuleTree) = map(sort!∘leaves, nodes(t))

function parent_node(t::StrongModuleTree, S::Vector)
    for x in t.nodes
        x isa StrongModuleTree || continue
        L = leaves(x)
        isempty(S ∩ L) && continue
        S ⊆ L && length(S) < length(L) && return parent_node(x, S)
        break
    end
    return t
end

## StrongModuleTree as a collection ##

Base.size(t::StrongModuleTree) = (length(t),)
Base.length(t::StrongModuleTree) = length(t.nodes)
Base.eltype(t::StrongModuleTree) = eltype(t.nodes)
Base.getindex(t::StrongModuleTree, i::Int) = t.nodes[i]

Base.getindex(v::Vector{T}, t::StrongModuleTree) where {T} =
    StrongModuleTree{T}(t.kind, t.edge, map(x->v[x], t.nodes))

function Base.sort!(t::StrongModuleTree; lt=isless, by=first_leaf, rev::Bool=false)
    for x in t.nodes
        x isa StrongModuleTree || continue
        sort!(x, lt=lt, by=by, rev=rev)
    end
    sort!(t.nodes, lt=lt, by=by, rev=rev)
    return t
end

end # module


text:
using PositiveFactorizations

function positive_linsolve(d, B, ∇f)
    cholesky!(Positive, B)
    Bchol = Cholesky(B, 'L', 0)
    d .= Bchol \ ∇f
end
function positive_linsolve(B, ∇f)
    Bchol = cholesky(Positive, B)
    Bchol \ ∇f
end
function positive_factorize(B)
    cholesky!(Positive, B)
    Cholesky(B, 'L', 0)
end

Base.summary(::NLSolvers.Newton{<:Direct,typeof(positive_linsolve)}) =
    "Newton's method with PositiveFactorizations.jl linsolve"

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
################################################################################
#                          TEST: Kinematics using DCM
################################################################################

for i = 1:samples
    # Create a random DCM.
    q      = Quaternion(randn(), randn(), randn(), randn())
    q      = q/norm(q)
    Dba0   = quat_to_dcm(q)

    # Create a random velocity vector.
    wba_a = randn(3)

    # Propagate the initial DCM using the sampled velocity vector.
    Δ   = 0.0001
    Dba = copy(Dba0)
    for k = 1:1000
        dDba = ddcm(Dba,Dba0*wba_a)

        Dba  = Dba + dDba*Δ
    end

    # In the end, the vector aligned with w must not change.
    v1 = Dba0*wba_a
    v2 = Dba*wba_a

    @test norm(v2-v1) < 1e-10
end


text:
using Tropical
using Test

@testset "Tropical.jl" begin
    # Write your tests here.

end

include("download_data.jl")


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)
    require_mpi()
    mpi_tensor_transpose_ = load_op_and_grad("./build/libMPITensor.so","mpi_tensor_transpose", multiple=true)
    row,col,ncol,val,n,rank,nt = convert_to_tensor(Any[row,col,ncol,val,n,rank,nt], [Int32,Int32,Int32,Float64,Int64,Int64,Int64])
    indices, vals = mpi_tensor_transpose_(row,col,ncol,val,n,rank,nt)
end

mpi_init()
row = [0;1]
col = [0;1;0;1]
ncol = [2;2]
val = [1.0;2.0;5.0;6.0]
n = 2
rank = 0
nt = 2

# TODO: specify your input parameters
u = mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)
sess = Session(); init(sess)
@show run(sess, u)

# uncomment it for testing gradients
error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    return sum(mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(10,20))
v_ = rand(10,20)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
module V
const NAMES = [
    "Subp"
    "CycBCdk1"
    "PP1"
    "pENSAt"
    "Gwlp"
    "PP2AB55"
    "Wee1"
    "Wee1pp"
    "Cdc25"
    "Cdc25pp"
]

#name2idx
for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
@testset "581.shortest-unsorted-continuous-subarray.jl" begin
    @test find_unsorted_subarray([2, 6, 4, 8, 10, 9, 15]) == 5
    @test find_unsorted_subarray([1, 2, 3, 4]) == 0
    @test find_unsorted_subarray([1]) == 0
    @test find_unsorted_subarray([18, 17, 1, 12, 10, 9, 19, 2, 11, 28, 7, 21, 8, 3, 27, 30,
                                  0, 14, 26, 4, 15, 29, 6, 16, 20, 24, 23, 5, 13, 25, 22]) ==
          31
end


text:
function fg(d::Distributions.Logistic, x)
    scaled = (x - d.μ)/2d.θ
    core = sech( scaled )^2 / 4d.θ
    core, core * -tanh(scaled) / d.θ
end
function fgw(d::Distributions.Logistic, x)
    scaled = (x - d.μ)/2d.θ
    core = sech( scaled )^2 / 4d.θ
    weighted_core = core*exp(x^2/2)
    weighted_core, x*weighted_core - weighted_core * tanh(scaled) / d.θ
end
function fw(d::Distributions.Logistic, x)
    scaled = (x - d.μ)/2d.θ
    sech( scaled )^2*exp(x^2/2) / 4d.θ
end


text:
using Pages
using JSON

@async Pages.start()

Endpoint("/hello") do request::HTTP.Request
    data = String(request.body)
    println("Parameters: $(data)")
    "Hello stranger!"
end

r = HTTP.request("GET", "http://192.168.2.11:8000/hello?id=abc")

println(String(r.body))

result = ""
Endpoint("/test", POST) do request::HTTP.Request
    data = String(request.body)
    #println("Parameters: $(data)")
    global result = data
    response = JSON.json(Dict(:data => data))
end

r = HTTP.request("POST", "http://192.168.2.11:8000/test", [], "12345.pdf")

println(String(r.body))


text:
function toposort(adj)
#function [seq] = toposort(adj)
#% TOPOSORT		A Topological ordering of nodes in a directed graph
#%
#%  [SEQ] = TOPOSORT(ADJ)
#%
#% Inputs :
#%    ADJ : Adjacency Matrix.
#%	   ADJ(i,j)==1 ==> there exists a directed edge
#%	   from i to j
#%
#% Outputs :
#%    SEQ : A topological ordered sequence of nodes.
#%          empty matrix if graph contains cycles.
#%
#% Usage Example :
#%		N=5;
#%		[l,u] = lu(rand(N));
#%		adj = ~diag(ones(1,N)) & u>0.5;
#%		seq = toposort(adj);
#%
#%
#% Note     :
#% See also
##
#% Uses :
#
#% Change History :
#% Date		Time		Prog	Note
#% 18-May-1998	 4:44 PM	ATC	Created under MATLAB 5.1.0.421
#
#% ATC = Ali Taylan Cemgil,
#% SNN - University of Nijmegen, Department of Medical Physics and Biophysics
#% e-mail : cemgil@mbfys.kun.nl

    N = size(adj,1)
    indeg = sum(adj,1)
    outdeg = sum(adj,2)
    seq = []

    for i=1:N,
        # Find nodes with indegree 0
        idx = find(indeg.==0)
        # If can't find than graph contains a cycle
        if isempty(idx)
            seq = [];
            break
        end
        # Remove the node with the max number of connections
        idx2 = indmax(outdeg[idx])
        indx = idx[idx2]
        seq = vcat(seq, indx)
        indeg[indx]=-1
        idx = find(adj[indx,:])
        indeg[idx] = indeg[idx]-1
    end
    return seq
end


text:
mutable struct PlanktonOutputWriter
    filepath::String
    write_log::Bool
    save_diags::Bool
    save_plankton::Bool
    diags_file::String
    plankton_file::String
    plankton_include::Tuple
    max_filesize::Number # in Bytes
    part_diags::Int
    part_plankton::Int
end

"""
    PlanktonOutputWriter(;dir = "./results",
                               diags_prefix = "diags",
                               plankton_prefix = "plankton",
                               write_log = false,
                               save_diags = false,
                               save_plankton = false,
                               plankton_include = (:x, :y, :z, :Sz),
                               max_filesize = Inf,
                               )
Generate a `PlanktonOutputWriter` structure which includes settings for model outputs

Keyword Arguments (Optional)
============================
- `dir`: The directory to store model outputs, "./results" by default
- `diags_prefix`: Descriptive filename prefixed to diagnostic output files.
- `plankton_prefix`: Descriptive filename prefixed to plankton output files.
- `write_log`: write model logs which contain global averages of simulated plankton, default: `false`.
- `save_diags`: write diagnostics to disk, default: `false`.
- `save_plankton`: write plankton to disk, default: `false`.
- `plankton_include`: list of plankton properties to save, default: `(:x, :y, :z, :Sz)`
- `max_filesize`: The writer will stop writing to the output file once the file size exceeds `max_filesize`,
                    and write to a new one with a consistent naming scheme ending in `part1`, `part2`, etc.
                    default: `Inf`.
"""
function PlanktonOutputWriter(;dir = "./results",
                               diags_prefix = "diags",
                               plankton_prefix = "plankton",
                               write_log = false,
                               save_diags = false,
                               save_plankton = false,
                               plankton_include = (:x, :y, :z, :Sz),
                               max_filesize = Inf
                               )

    isdir(dir) && rm(dir, recursive=true)
    mkdir(dir)

    diags_file = ""
    plankton_file = ""

    if save_diags
        diags_file = joinpath(dir, diags_prefix*".jld2")
    end
    if save_plankton
        plankton_file = joinpath(dir, plankton_prefix*".jld2")
    end

    return PlanktonOutputWriter(dir, write_log, save_diags, save_plankton, diags_file, plankton_file,
                                plankton_include, max_filesize, 1, 1)
end

function show(io::IO, writer::PlanktonOutputWriter)
    print(io, "PlanktonOutputWriter:\n",
              "├── files are saved at $(writer.filepath)\n",
              "├── $(save_diags_string(writer))\n",
              "├── $(save_inds_string(writer))\n",
              "├── write log: $(writer.write_log)\n",
              "└── Maximum file size: $(humanize_filesize(writer.max_filesize))"
              )
end

function save_diags_string(writer::PlanktonOutputWriter)
    if writer.save_diags
        return "diagnostics are saved as $(writer.diags_file)"
    else
        return "diagnostics are not saved"
    end
end
function save_inds_string(writer::PlanktonOutputWriter)
    if writer.save_plankton
        return "individuals are saved as $(writer.plankton_file)\n├── including: $(writer.plankton_include)"
    else
        return "individuals are not saved"
    end
end
function humanize_filesize(s::Number)
    suffix = ["B", "KiB", "MiB", "GiB", "TiB", "PiB", "EiB", "ZiB", "YiB"]
    biggest_suffix = suffix[1]
    value = Float64(s)
    unit = 1024.0
    for power in 1:length(suffix)
        unit = 1024^power
        biggest_suffix = suffix[power]
        value < unit && break
    end
    s = 1024 * value / unit
    return @sprintf("%3.1f %s", s, biggest_suffix)
end


text:
@testset "IO" begin
	  @testset "JLD2" begin
        function deepequal(a::Array, b::Array)
            typeof(a) == typeof(b) || return false
            N = length(a)
            for i in 1:N
                deepequal(a[i] , b[i]) || return false
            end
            return true
        end

        function deepequal(a, b)
            typeof(a) == typeof(b) || return false
            N = fieldcount(typeof(a))
            if N==0
                return a==b
            end
            for i in 1:N
                deepequal(getfield(a, i) , getfield(b, i)) || return false
            end
            return true
        end

        function deeptest(a::Array, b::Array)
            @test typeof(a) == typeof(b)
            N = length(a)
            for i in 1:N
                deeptest(a[i] , b[i])
            end
        end

        function deeptest(a, b)
            @test typeof(a) == typeof(b)
            N = fieldcount(typeof(a))
            if N==0
                @test a==b
            end
            for i in 1:N
                deeptest(getfield(a, i) , getfield(b, i))
                @test deepequal(a,b)
            end
        end

	      β = 10

        tgrid = CompositeGrid.LogDensedGrid(
            :gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb
            [0.0, β],# The grid is defined on [0.0, β]
            [0.0, β],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter.
            5,# N of log grid
            0.005, # niminum interval length of log grid
            5 # N of bottom layer
        )

        ############# FileIO API #################
        save("example.jld2", Dict("grid" => tgrid))

        d = load("example.jld2")
        dg = d["grid"]
        println(typeof(dg))
        println(dg.grid)
        deeptest(dg, tgrid)
        @test deepequal(dg,tgrid)

        ############# naive API ##################
        jldopen("example.jld2", "w") do file
            file["test"] = tgrid
        end

        jldopen("example.jld2", "r") do file
            g = file["test"]
            println(typeof(g))
            println(g.grid)
            deeptest(g, tgrid)
            @test deepequal(g,tgrid)
        end

        # remove the file
        rm("example.jld2")
    end
end


text:
# ---------------------------------------------------------------------------- #
#
#   readSU2.jl
#
#   Read SU2 meshes into luteos
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    readSU2_2D( flname::String )

Reads in 2D SU2 mesh.
"""
function readSU2_2D( flname::String )

  dim = 2

  # Open file
  fid = open( flname, "r" )

  # Check dimension
  x = "%"
  while x[1] == '%'
      x = readline(fid)
  end
  xx   = split(x)
  ndim = parse( Int64, xx[end] )
  if ndim != dim
      error(" Input mesh error: Input mesh is wrong dimension, should be 2D")
  end

  # Read number of elements and triangle connectivity
  x = "%"
  while x[1] == '%'
      x = readline(fid)
  end
  xx    = split( x )
  nelem = parse( Int64, xx[end] )
  t = fill( 0, nelem, 3 )
  for il in 1:nelem #x[1] != '%' || x[1] != 'N'
      x   = readline(fid)
      xx  = split(x)
      ind = parse(Int64,xx[end]) + 1

      t[ind,1] = parse(Int64,xx[2]) + 1
      t[ind,2] = parse(Int64,xx[3]) + 1
      t[ind,3] = parse(Int64,xx[4]) + 1
  end
  if parse(Int64,xx[1]) != 5
      error(" Input mesh error: Only triangles are supported for 2D")
  end

  # Read node locations
  x = "%"
  while x[1] == '%'
      x = readline(fid)
  end
  xx = split( x )
  np = parse( Int64, xx[end] )
  p  = fill( 0.0, np, 2 )
  for il in 1:np
      x   = readline(fid)
      xx  = split(x)
      ind = parse(Int64,xx[end]) + 1

      p[ind,1] = parse(Float64,xx[1])
      p[ind,2] = parse(Float64,xx[2])
  end

  # Read boundary elements
  x = "%"
  while x[1] == '%'
      x = readline(fid)
  end
  xx = split( x )
  nm = parse( Int64, xx[end] )
  nnm = fill( 0, nm )
  sm = fill( "tag", nm )
  belm = fill( fill(0, 0, 0), nm )
  nelm = fill( 0, nm )

  for jjb in 1:nm
      x         = readline( fid )
      xx        = split( x )
      sm[jjb]   = xx[end]

      x         = readline( fid )
      xx        = split( x )
      nel       = parse( Int64, xx[end] )
      nelm[jjb] = nel
      belm[jjb] = fill( 0, nel, 3 )

      for il in 1:nel
          x   = readline(fid)
          xx  = split(x)

          belm[jjb][il,1] = parse(Int64,xx[2]) + 1
          belm[jjb][il,2] = parse(Int64,xx[3]) + 1
          belm[jjb][il,3] = jjb
      end
  end

  close(fid)

  bel = fill( 0, sum(nelm), 3 )
  cnt = 1
  for jjb in 1:nm
      bel[cnt:(cnt+nelm[jjb]-1),:] = belm[jjb]
      cnt += nelm[jjb]
  end

  return (p, t, bel, sm)

end

"""
    readSU2_3D( flname::String )

Reads in 3D SU2 mesh.
"""
function readSU2_3D( flname::String )

    dim = 3

    # Open file
    fid = open( flname, "r" )

    # Check dimension
    x = "%"
    while x[1] == '%'
        x = readline(fid)
    end
    xx   = split(x)
    ndim = parse( Int64, xx[end] )
    if ndim != dim
        error(" Input mesh error: Input mesh is wrong dimension, should be 3D")
    end

    # Read number of elements and triangle connectivity
    x = "%"
    while x[1] == '%'
        x = readline(fid)
    end
    xx    = split( x )
    nelem = parse( Int64, xx[end] )
    t = fill( 0, nelem, 4 )
    for il in 1:nelem
        x   = readline(fid)
        xx  = split(x)
        ind = parse(Int64,xx[end]) + 1

        t[ind,1] = parse(Int64,xx[2]) + 1
        t[ind,2] = parse(Int64,xx[3]) + 1
        t[ind,3] = parse(Int64,xx[4]) + 1
        t[ind,4] = parse(Int64,xx[5]) + 1
    end
    if parse(Int64,xx[1]) != 10
        error(" Input mesh error: Only tetrahedra are supported for 3D")
    end

    # Read node locations
    x = "%"
    while x[1] == '%'
        x = readline(fid)
    end
    xx = split( x )
    np = parse( Int64, xx[end] )
    p  = fill( 0.0, np, 3 )
    for il in 1:np
        x   = readline(fid)
        xx  = split(x)
        ind = parse(Int64,xx[end]) + 1

        p[ind,1] = parse(Float64,xx[1])
        p[ind,2] = parse(Float64,xx[2])
        p[ind,3] = parse(Float64,xx[3])
    end

    # Read boundary elements
    x = "%"
    while x[1] == '%'
        x = readline(fid)
    end
    xx = split( x )
    nm = parse( Int64, xx[end] )
    nnm = fill( 0, nm )
    sm = fill( "tag", nm )
    belm = fill( fill(0, 0, 0), nm )
    nelm = fill( 0, nm )

    for jjb in 1:nm
        x         = readline( fid )
        xx        = split( x )
        sm[jjb]   = xx[end]

        x         = readline( fid )
        xx        = split( x )
        nel       = parse( Int64, xx[end] )
        nelm[jjb] = nel
        belm[jjb] = fill( 0, nel, 4 )

        for il in 1:nel
            x   = readline(fid)
            xx  = split(x)

            belm[jjb][il,1] = parse(Int64,xx[2]) + 1
            belm[jjb][il,2] = parse(Int64,xx[3]) + 1
            belm[jjb][il,3] = parse(Int64,xx[4]) + 1
            belm[jjb][il,4] = jjb
        end
    end

    close(fid)

    bel = fill( 0, sum(nelm), 4 )
    cnt = 1
    for jjb in 1:nm
        bel[cnt:(cnt+nelm[jjb]-1),:] = belm[jjb]
        cnt += nelm[jjb]
    end

    return (p, t, bel, sm)

end


text:
import MetNets._print_state_head

_print_state_head(io::IO, out::EPOut) = printstyled(io,
    string(
        " EPOut: ",
        "status=", status(out), ", ",
        "curriter=", curriter(out), ", ",
        "max_abs_av=", round(maximum(abs, av(out)); sigdigits = 4),
        "\n"
    ),
    color = MetNets.INFO_COLOR
)


text:
#!/usr/bin/env julia --startup-file=no

"""
Linear-Stability-Calculators: Quasi-Hydrostatic Model
=====================================================

Linear Stability Calculation of an inertially unstable jet using spectral methods

The basic state is that of a barotropic Bickley jet (to be generalized).
The perturbation is assumed to be periodic in both the zonal (x) and vertical (z) directions.

Calculations are done in serial.  See inertial_instability_mpi.py for parallel execution.
"""

include("src/quasihydrostatic_setup.jl")
include("src/parameters.jl")
include("src/Compute_Eigenvalues.jl")
include("src/Plotting_Functions.jl")

files   = Files(         nc = "data_spectrum_serial.nc",
                       json = "data_parameters_serial.json",
                 plotgrowth = "growth_serial.png",
                 plotslicem = "growth_slice_m_serial.png",
                plotmodes1D = "plot_modes1D_serial",
                plotmodes2D = "plot_modes2D_serial"
                )
grid    = Grid(Ly = 1000e3, Lz = 3e3, Ny = 100, θ₀ = π/32, method=Cheb())
physics = Physics(N = 1e-2, ν = 0.26, θ₀ = grid.θ₀, NT = 1)
jet     = Jet(grid, physics)

Output_Parameters(grid, physics, jet, files.json)

Neigs  = 10
dk, Nk = 1e-6, 41
dm, Nm = 1e-4, 150

ks = collect(0 :dk:(Nk-1)*dk)
ms = collect(dm:dm:Nm*dm)

ωs    = zeros(Complex, (Neigs, Nm, Nk))
modes = zeros(Complex, (3*grid.Ny+1, Neigs, Nm, Nk))

for (ik, k) in enumerate(ks)
    for (im, m) in enumerate(ms)
        A = build_matrix(k, m, grid, physics, jet)

        ωs[:, im, ik], modes[:, :, im, ik] = compute_spectrum(A, Neigs)

        #print("k = ", ks[ik], " m = ", ms[im], " growth = ", imag(ωs[1, im, ik]), "\n")
    end
end

save_spectrum(ωs, modes, ks, ms, Neigs, grid.y, grid.Ny, files.nc)

plot_growth_slice(files.nc, files.json, files.plotslicem)
plot_growth(      files.nc, files.json, files.plotgrowth)

plot_modes_1D(files.nc, files.json, files.plotmodes1D, Neigs)
plot_modes_2D(files.nc, files.json, files.plotmodes2D, Neigs)


text:
type Specifications
  source::AbstractString
  sink::AbstractString
  headerpattern::Regex
  sampletypes::Dict{AbstractString,Regex}
  conditions::Array{AbstractString, 1}
  crossmixes::Array{AbstractString,1}
  experiment_ids::Array{AbstractString,1}
  timepoints::Dict{AbstractString, Any}
  timecourses::Dict{AbstractString, AbstractString}
  errors::Dict{AbstractString, AbstractString}
  colors::Dict{AbstractString, AbstractString}
  include::Array{AbstractString,1}
end

function load_specifications(projectfolder)
  S = YAML.load(open(joinpath(projectfolder, "specifications.yml")))
  source = S["source"]
  sink = S["sink"]
  headerpattern = Regex(S["headerpattern"])
  sampletypes = [sample => Regex(pattern) for (sample,pattern) in S["sampletypes"]]
  conditions = S["conditions"]
  crossmixes = S["crossmixes"]
  experiment_ids = S["experiment_ids"]
  timepoints = S["timepoints"]
  timecourses = S["timecourses"]
  errors = S["errors"]
  colors = S["colors"]
  include = S["include"]
  specs = Specifications(source, sink, headerpattern, sampletypes, conditions, crossmixes, experiment_ids, timepoints, timecourses, errors, colors, include)
  return(specs)
end

max_stars_repo_name: el-uhu/PPSIM

max_issues_repo_name: el-uhu/PPSIM

max_forks_repo_name: el-uhu/PPSIM