text:
using GeometricBase
using Test

dt = Float64
nt = 10

@testset "$(rpad("Dataseries 1d (scalar data)",80))" begin
    nd = 1
    ni = 1

    ds = DataSeries(rand(dt), nt)
    @test typeof(ds) == typeof(DataSeries{dt,1}(nt, ni))
    @test typeof(ds) <: AbstractArray{dt,1}
    @test firstindex(ds)   == 0
    @test firstindex(ds,1) == firstindex(ds.d,1) - 1
    @test firstindex(ds,2) == 1
    @test lastindex(ds)    == nt
    @test lastindex(ds,1)  == lastindex(ds.d,1) - 1
    @test lastindex(ds,2)  == 1
    @test strides(ds)      == (1,)
    @test stride(ds,1)     == 1
    @test axes(ds)   == (0:nt,)
    @test axes(ds,1) == 0:nt
    @test axes(ds,2) == 1:1
    @test size(ds.d) == (nt+1,)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nt+1
    @test ndims(ds)  == 1
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    @test eachindex(ds) == 0:nt
    @test eachindex(IndexLinear(), ds) == 0:nt
    @test eachindex(IndexCartesian(), ds) == CartesianIndices((0:nt,))

    for i in 0:nt
        ds[i] = i
    end

    @test ds.d[1] == ds[0]
    @test ds.d[end] == ds[nt]
    @test ds.d[end] == ds[end]
    @test ds.d == collect(0:nt)
    @test parent(ds)[:] == parent(ds[:])

    reset!(ds)
    @test ds[0] == ds[end]

    d = rand(nt+1)
    ds = DataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)

    for i in 0:nt
        set_data!(ds, dt(i), i)
    end

   @test ds.d == collect(0:nt)
end

@testset "$(rpad("Dataseries 1d (vector-valued data)",80))" begin
    nd = 2
    ni = 1

    ds = DataSeries(rand(dt, nd), nt)
    @test typeof(ds) == typeof(DataSeries{Vector{dt},1}(nt, ni))
    @test typeof(ds) <: AbstractArray{Vector{dt},1}
    @test firstindex(ds)   == 0
    @test firstindex(ds,1) == firstindex(ds.d,1) - 1
    @test firstindex(ds,2) == 1
    @test lastindex(ds)    == nt
    @test lastindex(ds,1)  == lastindex(ds.d,1) - 1
    @test lastindex(ds,2)  == 1
    @test strides(ds)      == (1,)
    @test stride(ds,1)     == 1
    @test axes(ds)   == (0:nt,)
    @test axes(ds,1) == 0:nt
    @test axes(ds,2) == 1:1
    @test size(ds.d) == (nt+1,)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nt+1
    @test ndims(ds)  == 1
    @test eltype(ds) == Vector{dt}
    @test parent(ds) == ds.d

    @test eachindex(ds) == 0:nt
    @test eachindex(IndexLinear(), ds) == 0:nt
    @test eachindex(IndexCartesian(), ds) == CartesianIndices((0:nt,))

    for i in 0:nt
        ds[i] = [i, i^2]
    end

    @test ds.d[1] == ds[0]
    @test ds.d[end] == ds[nt]
    @test ds.d[end] == ds[end]
    @test ds.d == [Vector{dt}([i, i^2]) for i in 0:nt]
    @test parent(ds)[:] == parent(ds[:])

    reset!(ds)
    @test ds[0] == ds[end]

    d = rand(nt+1)
    ds = DataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)

    for i in 0:nt
        set_data!(ds, dt(i), i)
    end

   @test ds.d == collect(0:nt)
end

@testset "$(rpad("Dataseries 2d (scalar data)",80))" begin
    ni = 2

    ds = DataSeries(rand(dt, ni), nt, ni)
    @test typeof(ds) == typeof(DataSeries{dt,2}(nt, ni))
    @test typeof(ds) <: AbstractArray{dt,2}
    @test firstindex(ds)   == firstindex(ds.d)
    @test firstindex(ds,1) == firstindex(ds.d,1) - 1
    @test firstindex(ds,2) == 1
    @test lastindex(ds)    == lastindex(ds.d)
    @test lastindex(ds,1)  == lastindex(ds.d,1) - 1
    @test lastindex(ds,2)  == lastindex(ds.d,2)
    @test strides(ds)      == (1,nt+1)
    @test stride(ds,1)     == 1
    @test stride(ds,2)     == nt+1
    @test axes(ds)   == (0:nt, 1:ni)
    @test axes(ds,1) == 0:nt
    @test axes(ds,2) == 1:2
    @test size(ds.d) == (nt+1, ni)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nt+1
    @test size(ds.d,2) == ni
    @test ndims(ds)  == 2
    @test eltype(ds) == dt
    @test parent(ds) == ds.d

    for i in 0:nt
        set_data!(ds, dt(i), i, 1)
        set_data!(ds, dt(i)^2, i, 2)
    end

    @test ds.d[:,1] == collect(0:nt)
    @test ds.d[:,2] == collect(0:nt) .^ 2

    @test ds.d[1,1] == ds[0,1]
    @test ds.d[1,1:ni] == ds[0,1:nsamples(ds)]

    @test ds[end,1] == ds.d[end,1]
    @test ds[end,1] == ds[nt,1]
    @test ds[end,1:nsamples(ds)] == ds[nt,1:ni]

    @test ds.d[1,:] == ds[0,:]
    @test ds.d[:,1] == parent(ds[:,1])
    @test ds.d[:,:] == parent(ds[:,:])

    reset!(ds)
    @test ds[0,1] == ds[end,1]

    # tx = rand(nd)
    # ds[:,0] .= tx
    # @test ds.d[:,1] == tx

    d = rand(nt+1,ni)
    ds = DataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)
end

@testset "$(rpad("Dataseries 2d (vector-valued data)",80))" begin
    nd = 2
    ni = 5

    ds = DataSeries([rand(dt, nd) for i in 1:ni], nt)
    @test typeof(ds) == typeof(DataSeries{Vector{dt},2}(nt, ni))
    @test typeof(ds) <: AbstractArray{Vector{dt},2}
    @test firstindex(ds)   == firstindex(ds.d)
    @test firstindex(ds,1) == firstindex(ds.d,1) - 1
    @test firstindex(ds,2) == 1
    @test lastindex(ds)    == lastindex(ds.d)
    @test lastindex(ds,1)  == lastindex(ds.d,1) - 1
    @test lastindex(ds,2)  == lastindex(ds.d,2)
    @test strides(ds)      == (1,nt+1)
    @test stride(ds,1)     == 1
    @test stride(ds,2)     == nt+1
    @test axes(ds)   == (0:nt, 1:ni)
    @test axes(ds,1) == 0:nt
    @test axes(ds,2) == 1:ni
    @test size(ds.d) == (nt+1, ni)
    @test size(ds.d) == size(ds)
    @test size(ds.d,1) == nt+1
    @test size(ds.d,2) == ni
    @test ndims(ds)  == 2
    @test eltype(ds) == Vector{dt}
    @test parent(ds) == ds.d

    for k in 1:ni
        for i in 0:nt
            ds[i,k][1] = k*i
            ds[i,k][2] = k*i^2
        end
    end

    tx = zeros(dt, nd)
    for k in 1:ni
        for i in 0:nt
            get_data!(ds, tx, i, k)
            @test tx == Array{eltype(tx)}(k .* [i, i^2])
        end
    end

    @test ds.d[1,1] == ds[0,1]
    @test ds.d[1,1:ni] == ds[0,1:nsamples(ds)]

    @test ds[end,1] == ds.d[end,1]
    @test ds[end,1] == ds[nt,1]
    @test ds[end,1:nsamples(ds)] == ds[nt,1:ni]

    @test ds.d[1,:] == ds[0,:]
    @test ds.d[:,1] == parent(ds[:,1])
    @test ds.d[:,:] == parent(ds[:,:])

    reset!(ds)
    @test ds[0,1] == ds[end,1]

    # tx = rand(nd)
    # ds[:,0] .= tx
    # @test ds.d[:,1] == tx

    d = rand(nt+1,ni)
    ds = DataSeries(d)
    @test ds.d == d

    @test similar(ds).d == zero(d)
end

# TODO: Add tests for array-valued data

# @testset "$(rpad("Dataseries 3d",80))" begin
#     nd = 2
#     ni = 2

#     ds = DataSeries(dt, nd, nt, ni)
#     @test typeof(ds) <: AbstractArray{dt,3}
#     @test firstindex(ds)   == firstindex(ds.d)
#     @test firstindex(ds,1) == firstindex(ds.d,1)
#     @test firstindex(ds,2) == firstindex(ds.d,2) - 1
#     @test firstindex(ds,3) == firstindex(ds.d,3)
#     @test firstindex(ds,4) == 1
#     @test lastindex(ds)    == lastindex(ds.d)
#     @test lastindex(ds,1)  == lastindex(ds.d,1)
#     @test lastindex(ds,2)  == lastindex(ds.d,2) - 1
#     @test lastindex(ds,3)  == lastindex(ds.d,3)
#     @test lastindex(ds,4)  == 1
#     @test strides(ds)      == (1,nd,nd*(nt+1))
#     @test stride(ds,1)     == 1
#     @test stride(ds,2)     == nd
#     @test axes(ds)   == (1:nd, 0:nt, 1:ni)
#     @test axes(ds,1) == 1:nd
#     @test axes(ds,2) == 0:nt
#     @test axes(ds,3) == 1:ni
#     @test axes(ds,4) == 1:1
#     @test size(ds.d) == (nd, nt+1, ni)
#     @test size(ds.d) == size(ds)
#     @test size(ds.d,1) == nd
#     @test size(ds.d,2) == nt+1
#     @test size(ds.d,3) == ni
#     @test ndims(ds)  == 3
#     @test eltype(ds) == dt
#     @test parent(ds) == ds.d

#     for j in 1:ni
#         for i in 1:nt
#             ds[1,i,j] = j*i
#             ds[2,i,j] = j*i^2
#         end
#     end

#     tx = zeros(nd)
#     ty = zeros(nd,ni)
#     tz = zeros(nd,ni)
#     for i in 1:nt
#         for j in 1:ni
#             get_data!(ds, tx, i, j)
#             @test tx == Array{eltype(tx)}([j*i, j*i^2])
#             tz[1,j] = j*i
#             tz[2,j] = j*i^2
#         end
#         get_data!(ds, ty, i)
#         @test ty == tz
#     end

#     @test ds.d[1,1,1] == ds[1,0,1]
#     @test ds.d[1:ds.nd,1,1] == ds[1:ds.nd,0,1]

#     @test ds[1,end,1] == ds.d[1,end,1]
#     @test ds[1,nt,1] == ds[1,end,1]
#     @test ds[1:ds.nd,nt,1] == ds[1:ds.nd,end,1]

#     @test ds.d[:,1,1] == ds[:,0,1]
#     @test ds.d[:,1,:] == ds[:,0,:]

#     reset!(ds)
#     @test ds[1,0,1] == ds[1,end,1]

#     tx = rand(nd)
#     ds[:,0,1] .= tx
#     @test ds.d[:,1,1] == tx

#     d = rand(nd,nt+1,ni)
#     ds = DataSeries(d)
#     @test ds.d == d

#     @test similar(ds).d == zero(d)
# end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# types and constructors

if VERSION < v"0.3-"
    typealias UnitRange{T} Range1{T}
end

# We need a global variable to keep track of gnuplot's state
type GnuplotState
    running::Bool                # true when gnuplot is already running
    current::Int                 # current figure
    fid                          # pipe streams id
    tmpdir::AbstractString       # where to store data files
	gp_stdout::AbstractString    # store gnuplot's stdout
	gp_stderr::AbstractString    # store gnuplot's stderr
	gp_lasterror::AbstractString # store gnuplot's last error output
	gp_error::Bool               # true if last command resulted in gp error
    figs::Array                  # storage for all figures

    function GnuplotState(running,
    	current,
    	fid,
    	tmpdir,
		gp_stdout,
		gp_stderr,
		gp_lasterror,
		gp_error,
        figs)
        # Check to see if tmpdir exists, and create it if not
        # TODO: there has to be a simpler way to do this
        mkd = true
        try
            readdir(tmpdir)
            mkd = false
        catch
        end
        if mkd
            # Linux
            try
                mkdir(tmpdir)
                mkd = false
            catch
            end
        end
        if mkd
            # Windows
            # TODO: there has to be a simpler way to do this
            ttdir = ""
            for x in split(tmpdir,'/')
                if x != ""
                    try
                        mkdir(string(ttdir,x))
                    catch
                    end
                    ttdir = string(ttdir, x, '\\')
                end
            end
            mkd = false
        end
        new(running,current,fid,tmpdir,
        gp_stdout,gp_stderr,gp_lasterror,gp_error,figs)
    end
end

# Structure to keep Gaston's configuration
type GastonConfig
    # default CurveConf values
    legend::AbstractString
    plotstyle::AbstractString
    color::AbstractString
    marker::AbstractString
    linewidth::Real
    pointsize::Real
    # default AxesConf values
    title::AbstractString
    xlabel::AbstractString
    ylabel::AbstractString
    zlabel::AbstractString
    fill::AbstractString
    grid::AbstractString
    box::AbstractString
    axis::AbstractString
    xrange::AbstractString
    yrange::AbstractString
    zrange::AbstractString
    # default terminal type
    terminal::AbstractString
    # for terminals that support filenames
    outputfile::AbstractString
    # for printing to file
    print_color::AbstractString
    print_fontface::AbstractString
    print_fontsize::Real
    print_fontscale::Real
    print_linewidth::Real
    print_size::AbstractString
end
function GastonConfig()
	# Determine if current display supports PNG, and set defaults
	# accordingly. This allows IJulia support.

	if displayable("image/png")
		GastonConfig(
    	# CurveConf
	    "","lines","","",1,0.5,
    	# AxesConf
	    "Untitled","x","y","z","empty","off","inside vertical right top","",
	    "[*:*]","[*:*]","[*:*]",
    	# terminal
	    "png",
    	# output file name
	    "$(gnuplot_state.tmpdir)gaston-ijulia.png",
	    # print parameters
	    "color", "Sans", 10, 1, 1, "640,480"
    	)
	else
		GastonConfig(
    	# CurveConf
	    "","lines","","",1,0.5,
    	# AxesConf
	    "Untitled","x","y","z","empty","off","inside vertical right top","",
	    "[*:*]","[*:*]","[*:*]",
    	# terminal
	    "wxt",
    	# output file name
	    "",
	    # print parameters
	    "color", "Sans", 12, 0.5, 1, "640,480"
	    )
	end
end

# Structs to configure a plot
# Two types of configuration are needed: one to configure a single curve, and
# another to configure a set of curves (the 'axes').
type CurveConf
    legend::AbstractString          # legend text
    plotstyle::AbstractString
    color::AbstractString           # one of gnuplot's builtin colors --
                            # run 'show colornames' inside gnuplot
    marker::AbstractString          # point type
    linewidth::Real
    pointsize::Real

end
CurveConf() = CurveConf(
    gaston_config.legend,
    gaston_config.plotstyle,
    gaston_config.color,
    gaston_config.marker,
    gaston_config.linewidth,
    gaston_config.pointsize)

type AxesConf
    title::AbstractString      # plot title
    xlabel::AbstractString     # xlabel
    ylabel::AbstractString     # ylabel
    zlabel::AbstractString     # zlabel for 3-d plots
    fill::AbstractString       # fill style
    grid::AbstractString       # grid on/off
    box::AbstractString        # legend box (used with 'set key')
    axis::AbstractString       # normal, semilog{x,y}, loglog
    xrange::AbstractString     # xrange
    yrange::AbstractString     # yrange
    zrange::AbstractString     # zrange
end
AxesConf() = AxesConf(
    gaston_config.title,
    gaston_config.xlabel,
    gaston_config.ylabel,
    gaston_config.zlabel,
    gaston_config.fill,
    gaston_config.grid,
    gaston_config.box,
    gaston_config.axis,
    gaston_config.xrange,
    gaston_config.yrange,
    gaston_config.zrange)

type Financial
    open::Vector
    low::Vector
    high::Vector
    close::Vector
end
Financial() = Financial(Any[],Any[],Any[],Any[])

# coordinates and configuration for a single curve
type CurveData
    x::Vector          # abscissa
    y::Vector          # ordinate
    Z::Array           # 3-d plots and images
    ylow::Vector       # error data
    yhigh::Vector      # error data
    finance::Financial # financial data
    conf::CurveConf    # curve configuration
end
CurveData() = CurveData(Any[],Any[],Any[],Any[],Any[],Financial(),CurveConf())
CurveData(x,y,Z,conf::CurveConf) = CurveData(x,y,Z,Any[],Any[],Financial(),conf)

# curves and configuration for a single figure
type Figure
    handle::Int                  # each figure has a unique handle
    curves::Vector{CurveData}    # a vector of curves (x,y,conf)
    conf::AxesConf               # figure configuration
    isempty::Bool                # a flag to indicate if figure is empty
end
Figure(handle) = Figure(handle,CurveData[CurveData()],AxesConf(),true)

# coordinate type
Coord = Union{UnitRange,Range,Matrix,Vector}

ext: jl


text:
using BusinessDays

struct CustomCalendar <: HolidayCalendar end

BusinessDays.isholiday(::CustomCalendar, dt::Dates.Date) = dt == Dates.Date(2015,8,27)

cc = CustomCalendar()
println("Returns false: $(isholiday(cc, Dates.Date(2015,8,26)))")
println("Returns true: $(isholiday(cc, Dates.Date(2015,8,27)))")


text:
@doc raw"""
Arnold-Beltrami-Childress (ABC) field in (x,y,z) coordinates with covariant components of the vector
potential given by
```math
A (x,y,z) = \big( a \, \sin(z) + c \, \cos(y) , \, b \, \sin(x) + a \, \cos(z) , \, c \, \sin(y) + b \, \cos(x) \big)^T
```
resulting in the magnetic field ``B(x,y,z) = A(x,y,z)``.

Parameters: `a`, `b`, `c`
"""
module ABC

    using RecipesBase

    import ..ElectromagneticFields
    import ..ElectromagneticFields: CartesianEquilibrium, code
    import ..ElectromagneticFields: A₁, A₂, A₃
    import ..AnalyticCartesianField: X, Y, Z

    export ABCEquilibrium

    const DEFAULT_A = 1
    const DEFAULT_B = 1
    const DEFAULT_C = 1

    struct ABCEquilibrium{T <: Number} <: CartesianEquilibrium
        name::String
        a₀::T
        b₀::T
        c₀::T

        ABCEquilibrium{T}(a::T, b::T, c::T) where T <: Number = new("ABCEquilibrium", a, b, c)
    end

    ABCEquilibrium(a::T=DEFAULT_A, b::T=DEFAULT_B, c::T=DEFAULT_C) where T <: Number = ABCEquilibrium{T}(a, b, c)

    function init(a=DEFAULT_A, b=DEFAULT_B, c=DEFAULT_C)
        ABCEquilibrium(a, b, c)
    end

    macro code(a=DEFAULT_A, b=DEFAULT_B, c=DEFAULT_C)
        code(init(a, b, c); escape=true)
    end

    function Base.show(io::IO, equ::ABCEquilibrium)
        print(io, "ABC Equilibrium with\n")
        print(io, "  A = ", equ.a₀, "\n")
        print(io, "  B = ", equ.b₀, "\n")
        print(io, "  C = ", equ.c₀)
    end

    ElectromagneticFields.A₁(x::AbstractVector, equ::ABCEquilibrium) = equ.a₀ * sin(x[3]) + equ.c₀ * cos(x[2])
    ElectromagneticFields.A₂(x::AbstractVector, equ::ABCEquilibrium) = equ.b₀ * sin(x[1]) + equ.a₀ * cos(x[3])
    ElectromagneticFields.A₃(x::AbstractVector, equ::ABCEquilibrium) = equ.c₀ * sin(x[2]) + equ.b₀ * cos(x[1])
    B(x::AbstractVector, equ::ABCEquilibrium) = A₁(x,equ)^2 + A₂(x,equ)^2 + A₃(x,equ)^2

    ElectromagneticFields.get_functions(::ABCEquilibrium) = (X=X, Y=Y, Z=Z)

    @recipe function f(equ::ABCEquilibrium; nx=99, ni=div(nx,2)+1, nl=12, size=(400,1200))
        xmin = 0
        xmax = 2π
        grid = LinRange(xmin, xmax, nx)
        Bfield = [B([grid[i], grid[j], grid[k]], equ) for i in eachindex(grid), j in eachindex(grid), k in eachindex(grid)]

        seriestype   := :contour
        aspect_ratio := :equal
        layout := (3, 1)
        size   := size
        xlims  := (xmin,xmax)
        ylims  := (xmin,xmax)
        levels := nl
        legend := :none

        @series begin
            subplot := 1
            xguide := "x"
            yguide := "y"
            title  := "|B(x,y,π)|"
            (grid, grid, Bfield[:,:,ni])
        end

        @series begin
            subplot := 2
            xguide := "x"
            yguide := "z"
            title  := "|B(x,π,z)|"
            (grid, grid, Bfield[:,ni,:])
        end

        @series begin
            subplot := 3
            xguide := "y"
            yguide := "z"
            title  := "|B(π,y,z)|"
            (grid, grid, Bfield[ni,:,:])
        end
    end

end


text:
using DynamicPolynomials
using MultivariateMoments
using SemialgebraicSets

const MM = MultivariateMoments

include("../src/BeurlingLasso.jl")
using .BeurlingLasso

@polyvar x[1:2]

# only valid for odd numbers
function doublefactorial(n::Integer)
    if n <= 0
        1
    else
        prod(2*k+1 for k=0:div(n-1, 2))
    end
end

nmoments = 15

φ = [sum(binomial(n, 2*k) * x[1]^(n-2*k) * x[2]^(2*k) * doublefactorial(2*k-1)
     for k=0:div(n, 2)) for n=0:nmoments-1]

params = [[-0.5, 0.5], [0.6, 0.3]]
weights = [0.7, 0.3]

natoms = length(params)

μ0 = AtomicMeasure(x, [WeightedDiracMeasure(params[i], weights[i]) for i=1:natoms])
y = [MM.expectation(μ0, φ[n]) for n=1:nmoments]

# λ = 0
#
# domain = @set x[1] >= -1 && x[1] <= 1 && x[2] >= -1 && x[2] <= 1
#
# prob1 = blasso(y, φ, λ, domain)
# solver1 = momsossolver(1e-3)
# solve!(prob1, solver1)

λ = 1e-5

φ_func(u, n) = φ[n](x=>u)
bounds = [[-1 1]; [-1 1]]

prob2 = blasso(y, φ_func, λ, bounds)
solver2 = sfwsolver(2, 1000)
solve!(prob2, solver2)


text:
using Revise
using ImageInspector,MLDatasets, Plots

x = FashionMNIST.traintensor(1);

plot(
    plot(image(x; flip = true); title = "flip = true"),
    plot(image(x; flip = false); title = "flip = false");
    axis = nothing,
    border = :none,
)

x1 = FashionMNIST.traintensor(1);
x2 = CIFAR10.traintensor(2);
plot(
    plot(image(x1)),
    plot(image(x2));
    axis = nothing,
    border = :none
)

x = FashionMNIST.traintensor(1:10);
plot(plot.(image(x, [1,2]))...; axis = nothing, border = :none)

x = FashionMNIST.traintensor(1:10);
plot(imagegrid(x, 1:10; nrows = 2, sep = 2); axis = nothing, border = :none)

x = CIFAR10.traintensor(1:10);
imageplot(x, 1:10; nrows = 2, sep = 1, background = RGB(184/255, 223/255, 250/255))

max_stars_repo_name: NikolaGrover/ImageInspector

max_issues_repo_name: NikolaGrover/ImageInspector

max_forks_repo_name: NikolaGrover/ImageInspector


text:
using arXiv
using Test

@testset "arXiv.jl" begin
	include("request.jl")
end


text:
using RobustAndOptimalControl, ControlSystems
using FiniteDiff, Zygote

##
# Define the proces parameters
k1, k2, kc, g = 3.33, 3.35, 0.5, 981
A1, A3, A2, A4 = 28, 28, 32, 32
a1, a3, a2, a4= 0.071, 0.071, 0.057, 0.057
h01, h02, h03, h04 = 12.4, 12.7, 1.8, 1.4
T1, T2 = (A1/a1)*sqrt(2*h01/g), (A2/a2)*sqrt(2*h02/g)
T3, T4 = (A3/a3)*sqrt(2*h03/g), (A4/a4)*sqrt(2*h04/g)
c1, c2 = (T1*k1*kc/A1), (T2*k2*kc/A2)
γ1, γ2 = 0.7, 0.6

# Define the process dynamics
A = [-1/T1     0 A3/(A1*T3)          0;
     0     -1/T2          0 A4/(A2*T4);
     0         0      -1/T3          0;
     0         0          0      -1/T4];
B = [γ1*k1/A1     0;
     0                γ2*k2/A2;
     0                (1-γ2)*k2/A3;
     (1-γ1)*k1/A4 0              ];
C = [kc 0 0 0;
     0 kc 0 0];
D = zeros(2,2)
G = ss(A,B,C,D);

# Sensitivity weight function
WS = makeweight(100, (0.1, 1), 1/10) * I(2)

# Output sensitivity weight function
# WUelement = 5*tf([1,1],[0.1,1])
WU = tf(0.01) .* I(2)

# Complementary sensitivity weight function
# WT = tf([10,0.1],[1,1]) .* I(2)
WT = makeweight(1/10, (0.1, 1), 10) * I(2)

# Form augmented P dynamics in state-space
P = hinfpartition(G, WS, WU, WT)

# Check that the assumptions are satisfied
flag = hinfassumptions(P)

# Synthesize the H-infinity optimal controller
flag, C, γ = hinfsynthesize(P)

Pcl, S, CS, T = hinfsignals(P, G, C)
isinteractive() && specificationplot([S, CS, T], [WS[1,1], 0.01, WT[1,1]], γ)

function cost(v)
    K = vec2sys(v, 2, 2)
    T = lft(P, K)
    hinfnorm2(T)[1]
end

v = vec(C)

# v0 = 0.5*vec(C)
# cost(v0)

# K0 = vec2sys(v0, 2, 2)
# T0 = lft(P, K0)
# using Optim, Optim.LineSearches
# res = Optim.optimize(
#     cost,
#     v0,
#     ParticleSwarm(),
#     Optim.Options(
#         store_trace       = true,
#         show_trace        = true,
#         show_every        = 1,
#         iterations        = 1000,
#         allow_f_increases = false,
#         time_limit        = 1,
#         x_tol             = 0,
#         f_abstol          = 0,
#         g_tol             = 1e-8,
#         f_calls_limit     = 0,
#         g_calls_limit     = 0,
#     ),
# )

# vo = res.minimizer
# Ko = vec2sys(vo, 2, 2)
# To = lft(P, Ko)

# Pcl, S, CS, T = hinfsignals(P, G, C)
# specificationplot([S, CS, T], [WS[1,1], 0.01, WT[1,1]], 1, nsigma=1)

# Pcl, S, CS, T = hinfsignals(P, G, Ko)
# specificationplot!([S, CS, T], l=(:black, :dash), nsigma=1, ylims=(0.01, Inf))

##

function testfun(v)
    K = vec2sys(v, 2, 2)
    typeof(K)
    hinfnorm2(K, rtolinf=1e-8)[1]
end

testfun(v)

g1 = Zygote.gradient(testfun, v)[1]

g2 = FiniteDiff.finite_difference_gradient(testfun, v, absstep=1e-8)

senssys1 = vec2sys(g1, 2, 2)
senssys2 = vec2sys(g2, 2, 2)

@test norm(senssys1.A - senssys2.A) / norm(senssys1.A) < 0.001
@test norm(senssys1.B - senssys2.B) / norm(senssys1.B) < 0.001
@test norm(senssys1.C - senssys2.C) / norm(senssys1.C) < 0.001
@test norm(senssys1.D - senssys2.D) / norm(senssys1.D) < 0.001

## Broken due to Zygote not handling try/catch and @warn. The offending method is feedback
# g1 = Zygote.gradient(cost, v)[1]
# g2 = FiniteDiff.finite_difference_gradient(cost, v, absstep=1e-8)

# senssys1 = vec2sys(g1, 2, 2)
# senssys2 = vec2sys(g2, 2, 2)

# @test norm(senssys1.A - senssys2.A) / norm(senssys1.A) < 0.001
# @test norm(senssys1.B - senssys2.B) / norm(senssys1.B) < 0.001
# @test norm(senssys1.C - senssys2.C) / norm(senssys1.C) < 0.001
# @test norm(senssys1.D - senssys2.D) / norm(senssys1.D) < 0.001


text:
module test_millboard_enums

using Millboard
using Test

Millboard.set_table_mode(:markdown)

@enum RPS R P S

@test sprint(show, table(RPS)) == """
| RPS | VALUE |
|-----|-------|
| R   | 0     |
| P   | 1     |
| S   | 2     |"""

@test sprint(show, table(R)) == """
| RPS | VALUE |
|-----|-------|
| R   | 0     |
| P   | 1     |
| S   | 2     |"""

end # module test_millboard_enums


text:
#
# @defmodel and supporting functions
#
using MacroTools

"""
    defmodel(model_name::Symbol, ex::Expr)

Define a Mimi model. The following types of expressions are supported:

1. `component(name)`                            # add comp to model
2. `dst_component.name = ex::Expr`              # provide a value for a parameter
3. `src_component.name => dst_component.name`   # connect a variable to a parameter
4. `index[name] = iterable-of-values`           # define values for an index
"""
macro defmodel(model_name, ex)
    @capture(ex, elements__)

    # @__MODULE__ is evaluated in calling module when macro is interpreted
    result = :(
        let calling_module = @__MODULE__, comp_mod_name = nothing, comp_mod_obj = nothing
            global $model_name = Model()
        end
    )

    # helper function used in loop below
    function addexpr(expr)
        let_block = result.args[end].args
        push!(let_block, expr)
    end

    for elt in elements
        offset = 0

        if @capture(elt, component(comp_mod_name_.comp_name_)         | component(comp_name_) |
                         component(comp_mod_name_.comp_name_, alias_) | component(comp_name_, alias_))

            # set local copy of comp_mod_name to the stated or default component module
            expr = (comp_mod_name === nothing ? :(comp_mod_obj = calling_module) # nameof(calling_module))
                                              # TBD: This may still not be right:
                                              : :(comp_mod_obj = getfield(calling_module, $(QuoteNode(comp_mod_name)))))
            addexpr(expr)

            name = (alias === nothing ? comp_name : alias)
            expr = :(add_comp!($model_name, Mimi.ComponentId(comp_mod_obj, $(QuoteNode(comp_name))), $(QuoteNode(name))))

        # TBD: extend comp.var syntax to allow module name, e.g., FUND.economy.ygross
        elseif (@capture(elt, src_comp_.src_name_[arg_] => dst_comp_.dst_name_) ||
                @capture(elt, src_comp_.src_name_ => dst_comp_.dst_name_))
            if (arg !== nothing && (! @capture(arg, t - offset_) || offset <= 0))
                error("Subscripted connection source must have subscript [t - x] where x is an integer > 0")
            end

            expr = :(Mimi.connect_param!($model_name,
                                         $(QuoteNode(dst_comp)), $(QuoteNode(dst_name)),
                                         $(QuoteNode(src_comp)), $(QuoteNode(src_name)),
                                         offset=$offset))

        elseif @capture(elt, index[idx_name_] = rhs_)
            expr = :(Mimi.set_dimension!($model_name, $(QuoteNode(idx_name)), $rhs))

        # elseif @capture(elt, comp_name_.param_name_ = rhs_)
        #     expr = :(Mimi.set_param!($model_name, $(QuoteNode(comp_name)), $(QuoteNode(param_name)), $rhs))

        elseif @capture(elt, lhs_ = rhs_) && @capture(lhs, comp_.name_)
            (path, param_name) = parse_dotted_symbols(lhs)
            expr = :(Mimi.set_param!($model_name, $path, $(QuoteNode(param_name)), $rhs))
            # @info "expr: $expr"

        else
            # Pass through anything else to allow the user to define intermediate vars, etc.
            @info "Passing through: $elt"
            expr = elt
        end

        addexpr(expr)
    end

    # addexpr(:($model_name))     # return this or nothing?
    addexpr(:(nothing))
    return esc(result)
end


text:
using GeoData, Test, Aqua

if VERSION >= v"1.5.0"
    Aqua.test_ambiguities([GeoData, Base, Core])
    Aqua.test_unbound_args(GeoData)
    Aqua.test_undefined_exports(GeoData)
    Aqua.test_project_extras(GeoData)
    Aqua.test_deps_compat(GeoData)
    Aqua.test_project_toml_formatting(GeoData)
end

include("array.jl")
include("stack.jl")
include("series.jl")
include("utils.jl")
include("set.jl")
include("reproject.jl")
include("aggregate.jl")
include("methods.jl")
include("resample.jl")
# Only test SMAP locally for now
if !haskey(ENV, "CI")
    include("sources/smap.jl")
end
if !Sys.iswindows()
    # GDAL Environment vars need to be set manually for windows, so skip for now
    include("sources/gdal.jl")
    include("sources/grd.jl")
end
include("sources/ncdatasets.jl")


text:
using Plots
using LinearAlgebra
using DifferentialEquations
using ForwardDiff
using IntervalRootFinding
using StaticArrays
using Distances

# utilidades generales

inbox(x,y,xlims,ylims) = (xlims[1]<x<xlims[2]) & (ylims[1]<y<ylims[2])

function realplot!(p1,x,y;plotops...)
    idx = isreal.(x).*isreal.(y)
    plot!(p1,real(x[idx]),real(y[idx]);plotops...)
end

function realplot!(x,y;plotops...)
    idx = isreal.(x).*isreal.(y)
    plot!(real(x[idx]),real(y[idx]);plotops...)
end

function realplot(x,y;plotops...)
    p1 = plot()
    realplot!(p1,x,y;plotops...)
end

# Flujos 1D

function dualplot_flux1D(f,x,x0,p,sol;
    size=(900,300),plotops...)

    (ymin,ymax)=extrema(sol.u)
    (fmin,fmax)=extrema(f.(x,(p,),0.0))
    p1 = plot([x[1],x[end]],[0,0],label="",color=:gray)
    if (x[1]*x[end]<0)
        plot!(p1,[0,0],[fmin,fmax],label="",color=:gray,linestyle=:dash)
    end
    plot!(p1,x,f.(x,(p,),0.0),label="f(x)",xlabel="x",ylabel="f(x)",color=:blue)
    plot!(p1,sol.u,sol.u*0,label="x",color=:red)
    scatter!(p1,[x0],[0],label="x0",color=:green)
    p2 = plot(sol,label="x(t)",ylabel="x",ylim=(0.9*ymin,1.1*ymax))
    plot(p1,p2,layout=(1,2);size=size,fmt=:png,plotops...)
end

function flux1D(f::Function,x0::Float64,tmax::Float64,p;
    xlims=[-1.0,1.0],plotops...)

    xrange = xlims[2]-xlims[1]
    x=xlims[1]:xrange/100:xlims[2]
    prob = ODEProblem(f,x0,(0,tmax),p)
    sol = solve(prob)
    dualplot_flux1D(f,x,x0,p,sol;plotops...)
end

function flux1D(f::Function,x0::Float64,tmax::Float64,p,condition::Function;
    xlims=[-1.0,1.0],plotops...)

    xrange = xlims[2]-xlims[1]
    x=xlims[1]:xrange/100:xlims[2]
    prob = ODEProblem(f,x0,(0,tmax),p)
    condition_terminate(u,t,integrator) = condition(u)
    affect!(integrator) = terminate!(integrator)
    sol = solve(prob,callback=DiscreteCallback(condition_terminate,affect!))
    dualplot_flux1D(f,x,x0,p,sol;plotops...)
end

function flux1D(f::Function,x0::Float64,tmax::Float64,p,tperturb::Float64,Aperturb::Float64,condition::Function;
    xlims=[-1.0,1.0],plotops...)

    xrange = xlims[2]-xlims[1]
    x=xlims[1]:xrange/100:xlims[2]
    prob = ODEProblem(f,x0,(0,tmax),p)
    condition_terminate(u,t,integrator) = condition(u)
    affect1!(integrator) = terminate!(integrator)
    affect2!(integrator) = integrator.u += Aperturb
    cb1 = DiscreteCallback(condition_terminate,affect1!)
    cb2 = PresetTimeCallback([tperturb],affect2!)
    sol = solve(prob,callback=CallbackSet(cb1,cb2))
    dualplot_flux1D(f,x,x0,p,sol;plotops...)
end

function potential1D(V::Function,x0::Float64,tmax::Float64,p;
    xlims=[-1.0,1.0],size=(900,300),plotops...)

    xrange = xlims[2]-xlims[1]
    x=xlims[1]:xrange/100:xlims[2]
    f(x,p,t) = -ForwardDiff.derivative(x -> V(x,p), x)
    prob = ODEProblem(f,x0,(0,tmax),p)
    sol = solve(prob)
    p1 = plot(x,V.(x,(p,)),label="V(x)",xlabel="x",ylabel="V")
    scatter!(p1,[x0],[V(x0,p)],label="x0")
    plot!(p1,sol,V.(sol.u,(p,)),color=:black,linewidth=2,label="")
    p2 = plot(sol,label="x(t)",ylabel="x",ylim=ylims)
    plot(p1,p2,layout=(1,2);size=size,fmt=:png,plotops...)
end

function potential1D(V::Function,x0::Float64,tmax::Float64,p,tperturb::Float64,Aperturb::Float64,condition::Function;
    xlims=[-1.0,1.0],size=(900,300),plotops...)

    xrange = xlims[2]-xlims[1]
    x=xlims[1]:xrange/100:xlims[2]
    f(x,p,t) = -ForwardDiff.derivative(x -> V(x,p), x)
    prob = ODEProblem(f,x0,(0,tmax),p)
    condition_terminate(u,t,integrator) = condition(u)
    affect1!(integrator) = terminate!(integrator)
    affect2!(integrator) = integrator.u += Aperturb
    cb1 = DiscreteCallback(condition_terminate,affect1!)
    cb2 = PresetTimeCallback([tperturb],affect2!)
    sol = solve(prob,callback=CallbackSet(cb1,cb2))
    p1 = plot(x,V.(x,(p,)),label="V(x)",xlabel="x",ylabel="V")
    scatter!(p1,[x0],[V(x0,p)],label="x0")
    plot!(p1,sol,V.(sol.u,(p,)),color=:black,linewidth=2,label="")
    p2 = plot(sol,label="x(t)",ylabel="x",ylim=ylims)
    plot(p1,p2,layout=(1,2);size=size,fmt=:png,plotops...)
end

# 2D

function myquiver!(p1,x, y, u, v;
    scale=0.3,arrowscale=0.07, color=:black, linealpha=1,plotops...)

    function arrow0!(p1, x, y, u, v, as, lc, la)
        nuv = sqrt(u^2 + v^2)
        v1, v2 = [u;v] / nuv,  [-v;u] / nuv
        v4 = (3*v1 + v2)/3.1623
        v5 = v4 - 2*(v4'*v2)*v2
        v4, v5 = as*nuv*v4, as*nuv*v5
        plot!(p1,[x,x+u], [y,y+v], lc=lc,la=la)
        plot!(p1,[x+u,x+u-v5[1]], [y+v,y+v-v5[2]], lc=lc, la=la)
        plot!(p1,[x+u,x+u-v4[1]], [y+v,y+v-v4[2]], lc=lc, la=la)
    end
    scatter!(p1,x, y, mc=color, ms=1.5, ratio=1, ma=0.5;legend=false,plotops...)
    for (x,y,u,v) in zip(x,y,u,v)
        arrow0!(p1,x,y,u*scale,v*scale,arrowscale,color,linealpha)
    end
    p1
end

function myquiver(x, y, u, v;
    scale=0.3,arrowscale=0.07, color=:black, linealpha=1,plotops...)

    p1 = plot(xlabel="x",ylabel="y",fmt=:png)
    myquiver!(p1,x, y, u, v; scale=scale,arrowscale=arrowscale, color=color, linealpha=linealpha,plotops...)
end

function flux2d_grid(f,p,xlims,ylims,npts)
    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    dx = xrange/(npts-1)
    dy = yrange/(npts-1)
    xgrid = xlims[1]:dx:xlims[2]
    ygrid = ylims[1]:dy:ylims[2]
    u = [[x;y] for y in ygrid, x in xgrid]
    Z = u*0
    f.(Z,u,(p,),0)
    dx,dy,xgrid,ygrid,Z
end

function flux2d_vectorfield!(p1,f::Function,p;
    xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,scale=1.0,plotops...)

    dx,dy,xgrid,ygrid,Z = flux2d_grid(f,p,xlims,ylims,npts)
    u = vec(getindex.(Z,1))
    v = vec(getindex.(Z,2))
    x = vec(repeat(xgrid',length(ygrid)))
    y = repeat(ygrid,length(xgrid))
    maxuv = max(maximum(abs.(u)),maximum(abs.(v)))
    delta = min(dx,dy)
    scale = scale*delta/maxuv
    myquiver!(p1,x,y,u,v;scale=scale,arrowscale=0.5,color=:blue,linealpha=0.5,plotops...)
end

function flux2d_vectorfield(f::Function,p;
    xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,scale=1.0,size=(600,600),plotops...)

    p1 = plot(xlabel="x",ylabel="y",xlims=xlims,ylims=ylims,fmt=:png,legend=:none,size=size)
    flux2d_vectorfield!(p1,f,p;xlims=xlims,ylims=ylims,npts=npts,scale=scale,plotops...)
end

function flux2d_vectorfield(f::Function,u0::Vector{Float64},tmax::Float64,p;
    xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,scale=1.0,size=(600,600),plotops...)

    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    p1 = flux2d_vectorfield(f,p;xlims=xlims,ylims=ylims,npts=npts,scale=scale,size=size,plotops...)
    condition(u,t,integrator) = (u[1]*u[1]+u[2]*u[2]) > max(xrange*xrange,yrange*yrange)
    affect!(integrator) = terminate!(integrator)
    sol = solve(ODEProblem(f,u0,(0.0,tmax),p),callback=DiscreteCallback(condition,affect!))
    plot!(p1,sol,vars=(1,2),c=:black,arrow=true,xlims=xlims,ylims=ylims)
end

function flux2d_vectorfield(f::Function,u0_array::Vector{Vector{Float64}},tmax::Float64,p;
    xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,scale=1.0,size=(600,600),plotops...)

    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    p1 = flux2d_vectorfield(f,p;xlims=xlims,ylims=ylims,npts=npts,scale=scale,size=size,plotops...)
    prob = ODEProblem(f,u0_array[1],(0.0,tmax),p)
    ensamble_prob = EnsembleProblem(prob,prob_func=(prob,i,repeat;u0=u0_array)->(remake(prob,u0=u0[i])))
    condition(u,t,integrator) = (u[1]*u[1]+u[2]*u[2]) > max(xrange*xrange,yrange*yrange)
    affect!(integrator) = terminate!(integrator)
    sol = solve(ensamble_prob,EnsembleThreads(),trajectories=length(u0_array),callback=DiscreteCallback(condition,affect!))
    plot!(p1,sol,vars=(1,2),arrows=true,c=:black,linewidth=0.5,xlims=xlims,ylims=ylims)
end

function flux2d_nullclines(f::Function,p;
    xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,regions=true,vectorfield=false,size=(700,500),plotops...)
    #=
    =#
    _,_,xgrid,ygrid,Z = flux2d_grid(f,p,xlims,ylims,npts)
    p1 = plot(;xlabel="x",ylabel="y",legend=:none,fmt=:png,size=size,plotops...)
    if regions
        contourf!(p1,xgrid,ygrid,getindex.(Z,1),levels=[-1000,0,1000],alpha=0.8, c=[:red,:green])
        contourf!(p1,xgrid,ygrid,getindex.(Z,2),levels=[-1000,0,1000],alpha=0.4, c=[:blue,:yellow])
        contour!(p1,xgrid,ygrid,getindex.(Z,1),levels=[0],c=:red)
        contour!(p1,xgrid,ygrid,getindex.(Z,2),levels=[0],c=:green)
    else
        contour!(p1,xgrid,ygrid,getindex.(Z,1),levels=[0],c=:red,alpha=0.5)
        contour!(p1,xgrid,ygrid,getindex.(Z,2),levels=[0],c=:green,alpha=0.5)
    end
    if vectorfield
        flux2d_vectorfield!(p1,f,p;xlims=0.9*xlims,ylims=0.9*ylims)
    end
    p1
end

function flux2d_nullclines(f::Function,u0::Vector{Float64},tmax::Float64,p;
        xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,regions=true,vectorfield=false,size=(700,500),plotops...)
    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    p1 = flux2d_nullclines(f,p;xlims=xlims,ylims=ylims,npts=npts,regions=regions,vectorfield=vectorfield,size=size,plotops...)
    condition(u,t,integrator) = (u[1]*u[1]+u[2]*u[2]) > max(xrange*xrange,yrange*yrange)
    affect!(integrator) = terminate!(integrator)
    sol = solve(ODEProblem(f,u0,(0.0,tmax),p),callback=DiscreteCallback(condition,affect!))
    plot!(p1,sol,vars=(1,2),c=:black,arrow=true,xlims=xlims,ylims=ylims)
end

function flux2d_nullclines(f::Function,u0_array::Vector{Vector{Float64}},tmax::Float64,p;
    xlims=[-1.0,1.0],ylims=[-1.0,1.0],npts=21,regions=true,vectorfield=false,size=(700,500),plotops...)
    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    #u0_arr = vec([[xlims[1]+i*xrange/Ngrid,ylims[1]+j*yrange/Ngrid] for i=0:Ngrid, j=0:Ngrid])
    p1 = flux2d_nullclines(f,p;xlims=xlims,ylims=ylims,npts=npts,regions=regions,vectorfield=vectorfield,size=size,plotops...)
    prob = ODEProblem(f,u0_array[1],(0.0,tmax),p)
    ensamble_prob = EnsembleProblem(prob,prob_func=(prob,i,repeat;u0=u0_array)->(remake(prob,u0=u0[i])))
    condition(u,t,integrator) = (u[1]*u[1]+u[2]*u[2]) > max(xrange*xrange,yrange*yrange)
    affect!(integrator) = terminate!(integrator)
    sol = solve(ensamble_prob,EnsembleThreads(),trajectories=length(u0_array),callback=DiscreteCallback(condition,affect!))
    plot!(p1,sol,vars=(1,2),arrows=true,c=:black,linewidth=0.5,xlims=xlims,ylims=ylims)
end

function flux2d_animated(f::Function,p,N::Int64,dt::Float64;
    Ngrid=10,fps=15,xlims=[-1.0,1.0],ylims=[-1.0,1.0],size=(400,400),nullclines=false,fname="")

    if isempty(fname)
        fname = string(Symbol(f))
        if fname[end] == '!'
            fname = chop(fname)
        end
    end
    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    u0_arr = vec([[xlims[1]+i*xrange/Ngrid,ylims[1]+j*yrange/Ngrid] for i=0:Ngrid, j=0:Ngrid])
    prob = ODEProblem(f,u0_arr[1],(0.0,N*dt),p)
    ensamble_prob = EnsembleProblem(prob,prob_func=(prob,i,repeat;u0=u0_arr)->(remake(prob,u0=u0[i])))
    sol = solve(ensamble_prob,EnsembleThreads(),trajectories=length(u0_arr))
    x = reduce(hcat,[sol[n](0:dt:N*dt,idxs=1).u for n=1:length(sol)])
    y = reduce(hcat,[sol[n](0:dt:N*dt,idxs=2).u for n=1:length(sol)])
    if (nullclines)
        p1 = flux2d_nullclines(f,p;xlims=xlims,ylims=ylims,size=size)
    else
        p1 = plot([x[1:2,:]],[y[1:2,:]],color=:black,xlims=xlims,ylims=ylims,legend=false,size=size)
    end
    # make animation
    anim = @animate for n=3:N
        plot!(p1,x[n-1:n,:],y[n-1:n,:],color=:black)
    end
    print("Saving animation ...")
    gif(anim, fname * "_fps_" * string(fps) * ".gif", fps=fps)
end

function classification_linear(A::Matrix{Float64};
    Ngrid=5,tmax=2.0,xlims=[-1.2,1.2],ylims=[-1.2,1.2],circular=true)

    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    if circular
        Npts = Ngrid*Ngrid
        u0_arr = u0_arr = vec([[cos(i*2pi/Npts);sin(i*2pi/Npts)] for i=1:Npts])
    else
        u0_arr = vec([[xlims[1]+i*xrange/Ngrid,ylims[1]+j*yrange/Ngrid] for i=0:Ngrid, j=0:Ngrid])
    end
    prob = ODEProblem((u,p,t) -> A*u, [0;0], (0,tmax))
    ensamble_prob = EnsembleProblem(prob,prob_func=(prob,i,repeat;u0=u0_arr)->(remake(prob,u0=u0[i])))
    sol = solve(ensamble_prob,EnsembleThreads(),trajectories=length(u0_arr))
    p2 = plot(sol,vars=(1,2),xlims=xlims,ylims=ylims,arrows=true,c=:black,linewidth=0.5,fmt=:png)
    maxtr = 2.0
    maxdet = 1.3
    trv = -maxtr:maxtr/30:maxtr
    trdet = abs2.(trv)/4
    p1 = plot(trv,trdet,xaxis=("Tr(A)",(-maxtr,maxtr)),yaxis=("Det(A)",(-maxdet,maxdet)),legend=false,fmt = :png)
    plot!(p1,[-maxtr,maxtr],[0,0])
    plot!(p1,[0,0],[0,maxdet])
    txt=["Saddle","Stable\nNode","Stable\nFocus","Unstable\nFocus","Unstable\nNode"]
    annotate!([(0,-0.5*maxdet,(txt[1],12,:gray))])
    annotate!([(-0.8*maxtr,0.2*maxdet,(txt[2],10,:gray)),(-0.5*maxtr,0.7*maxdet,(txt[3],10,:gray))])
    annotate!([(0.8*maxtr,0.2*maxdet,(txt[5],10,:gray)),(0.5*maxtr,0.7*maxdet,(txt[4],10,:gray))])
    atxt = string(A[1,1]) * " " * string(A[1,2]) * "\n" * string(A[2,1]) * " " * string(A[2,2])
    scatter!(p1,[tr(A)],[det(A)],c=:red)
    annotate!([(tr(A),det(A)-0.1*maxdet,(atxt,9,:red,:center))])
    plot(p2,p1,layout=(1,2),size=(900,400))
end

function flux2d_manifolds!(p1,f,f_jac,u0_array,p;
    tmax=30,delta=0.001,repulsor=false,xlims=[-1.0,1.0],ylims=[-1.0,1.0],size=(700,500),plotops...)

    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    condition(u,t,integrator) = (u[1]*u[1]+u[2]*u[2]) > max(xrange*xrange,yrange*yrange)
    affect!(integrator) = terminate!(integrator)
    for u0 in u0_array
        A = f_jac(u0,p)
        if det(A)<0
            av = eigen(A)
            for n=1:2
                if real(av.values[n])>0
                    u1 = u0+delta*av.vectors[:,n]
                    sol = solve(ODEProblem(f,u1,(0.0,tmax),p),callback=DiscreteCallback(condition,affect!))
                    plot!(p1,sol,vars=(1,2),c=:red,label="Wu")
                    u1 = u0-delta*av.vectors[:,n]
                    sol = solve(ODEProblem(f,u1,(0.0,tmax),p),callback=DiscreteCallback(condition,affect!))
                    plot!(p1,sol,vars=(1,2),c=:red,label="")
                else
                    u1 = u0+delta*av.vectors[:,n]
                    sol = solve(ODEProblem(f,u1,(0.0,-tmax),p),callback=DiscreteCallback(condition,affect!))
                    plot!(p1,sol,vars=(1,2),c=:blue,label="Ws")
                    u1 = u0-delta*av.vectors[:,n]
                    sol = solve(ODEProblem(f,u1,(0.0,-tmax),p),callback=DiscreteCallback(condition,affect!))
                    plot!(p1,sol,vars=(1,2),c=:blue,label="")
                end
            end
        else
            if tr(A)<0
                scatter!(p1,[u0[1]],[u0[2]],c=:black,label="A")
            else
                scatter!(p1,[u0[1]],[u0[2]],c=:red,label="R")
                if repulsor
                    u1 = u0.+[delta;delta]
                    sol = solve(ODEProblem(f,u1,(0.0,3*tmax),p),callback=DiscreteCallback(condition,affect!))
                    plot!(p1,sol,vars=(1,2),c=:purple,alpha=0.5,label="")
                end
            end
        end
    end
    plot(p1;legend=false,xlims=xlims,ylims=ylims,size=size,plotops...)
end

function flux2d_manifolds(f::Function,f_jac::Function,u0_array::Vector{Vector{Float64}},p::Vector{Float64};
    tmax=30,delta=0.001,repulsor=false,xlims=[-1.0,1.0],ylims=[-1.0,1.0],size=(700,500),plotops...)

    p1=plot(xlabel="x",ylabel="y",fmt=:png)
    flux2d_manifolds!(p1,f,f_jac,u0_array,p;tmax=tmax,delta=delta,repulsor=repulsor,xlims=xlims,ylims=ylims,size=size,plotops...)
end

function phase_portrait(f::Function,p::Vector{Float64};
    tmax=50,delta=0.001,xlims=[-1.0,1.0],ylims=[-1.0,1.0],size=(700,500),plotops...)
    # Find fixedpoints in interval
    function fsv((x,y))
        du = f(similar([x,y]),[x,y],p,0)
        return SVector(du[1],du[2])
    end
    X = xlims[1]..xlims[2]
    Y = ylims[1]..ylims[2]
    rts = roots(fsv, X × Y)
    u0_arr=[[mid(rt.interval[1]);mid(rt.interval[2])] for rt in rts]
    p1 = flux2d_nullclines(f,p;xlims=xlims,ylims=ylims,npts=50,regions=false)
    f_jac(u0,p) = ForwardDiff.jacobian(x -> f(similar(x),x,p,0), u0)
    flux2d_manifolds!(p1,f,f_jac,u0_arr,p;tmax=tmax,delta=delta,repulsor=true,xlims=xlims,ylims=ylims,size=size,plotops...)
end

function attractor_basin(f::Function,p::Vector{Float64},attractors::Vector{Vector{Float64}},maxdist::Float64;
    delta=0.1,tmax=1000.0,xlims=(-1.0,1.0),ylims=(-1.0,1.0),size=(900,600),plotops...)

    natt = length(attractors)
    clist = [:black,:red,:blue,:yellow,:green,:purple,:cyan,:orange]
    if (natt>7)
        error("maximum number of attractors is 7")
    end
    distance2D(x1,x2) = sqrt((x1[1]-x2[1])^2+(x1[2]-x2[2])^2)
    x = xlims[1]:delta:xlims[2]
    y = ylims[1]:delta:ylims[2]
    nx = length(x)
    ny = length(y)
    u0_list=[[xi,yi,0.0] for xi in x for yi in y]
    prob = ODEProblem(f,u0_list[1],(0,tmax),p)
    ensamble_prob = EnsembleProblem(prob,prob_func=(prob,i,repeat;u0=u0_list)->(remake(prob,u0=u0_list[i])))
    sol = solve(ensamble_prob,EnsembleThreads(),trajectories=length(u0_list),save_everystep=false)
    println("solved")
    M = zeros(Int8,nx,ny)
    for (n,pt) in enumerate(sol)
        for (m,at) in enumerate(attractors)
            if (distance2D(pt.u[end][1:2],at)<maxdist)
                M[(n-1)÷ny+1,mod(n-1,ny)+1] = m
            end
        end
    end
    M[1,1]=0
    println("plotting...")
    contourf(x,y,M',c=clist[1:natt+1],linewidth=0,legend=false,xlabel="x",ylabel="y";size=size,fmt=:png,plotops...)
end

function flux2d_forced(f::Function,u0::Vector{Float64},p::Vector{Float64},period::Float64;
    tcycles=0, npts=100,ncycles=10,size=(900,400),xlims=false,ylims=false,plotops...)
    # Assume that the 3rd variable (z) is periodic and plot x,y as a function of z modulo period
    # skip trans period of the forcing
    tmax = min(ncycles*period,100000)
    sol = solve(ODEProblem(f,u0,(0.0,tmax),p))
    ncycles = Int(floor(tmax/period))
    p1 = plot(xlabel="x",ylabel="y",zlabel="ϕ/2π",legend=false,fmt = :png)
    p2 = plot(xlabel="x",ylabel="y",legend=false,fmt = :png)
    for j = tcycles+1:ncycles
        ts = range((j-1)*period,stop=j*period,length=npts+1)
        plot!(p1,sol(ts,idxs=1),sol(ts,idxs=2),(0:npts)/npts,color=j)
        scatter!(p1,[sol(ts[end],idxs=1)],[sol(ts[end],idxs=2)],[1],markersize=2,color=j+1)
        scatter!(p1,[sol(ts[1],idxs=1)],[sol(ts[1],idxs=2)],[0],markersize=2,color=j,framestyle=:box)
        plot!(p2,sol(ts,idxs=1),sol(ts,idxs=2),color=j,linealpha=0.3)
        scatter!(p2,[sol(ts[1],idxs=1)],[sol(ts[1],idxs=2)],markersize=2,color=j)
    end
    scatter!(p2,[sol(ncycles*period,idxs=1)],[sol(ncycles*period,idxs=2)],markersize=2,color=ncycles+1)
    if xlims isa Tuple
        xlims!(p1,xlims)
        xlims!(p2,xlims)
    end
    if ylims isa Tuple
        ylims!(p1,ylims)
        ylims!(p2,ylims)
    end
    plot(p1,p2,layout=(1,2);size=size,plotops...)
end

function poincare_forced!(p1,f, u0, p, period;
    tcycles=0, ncycles=10,msize=0.5,col=:blue,xlims=false,ylims=false,plotops...)
    # Assume that the 3rd variable (z) is periodic and plot x,y as a function of z modulo period
    # skip trans period of the forcing
    trans = solve(ODEProblem(f,u0,(0.0,tcycles*period),p))
    u0=trans.u[end]
    if period>1000
        tmax=min(ncycles*period,100000)
    else
        tmax=ncycles*period
    end
    tsave = period:period:ncycles*period
    sol = solve(ODEProblem(f,u0,(0.0,tmax),p),saveat=tsave)
    scatter!(p1,sol,vars=(1,2);legend=false,markersize=msize,color=col,markerstrokewidth=0,plotops...)
    if (xlims)
        (xmin,xmax) = extrema(getindex.(vcat(trans.u,sol.u),1))
        xrange=xmax-xmin
        xlims!(p1,xmin-0.1*xrange,xmax+0.1*xrange)
    end
    if (ylims)
        (ymin,ymax) = extrema(getindex.(vcat(trans.u,sol.u),2))
        yrange=ymax-ymin
        ylims!(p1,ymin-0.1*yrange,ymax+0.1*yrange)
    end
    p1
end

function poincare_forced(f::Function,u0::Vector{Float64},p::Vector{Float64},period::Float64;
    tcycles=0, ncycles=10,size=(500,500),msize=0.5,col=:blue,xlims=false,ylims=false,plotops...)

    p1=plot(xlabel="x",ylabel="y",size=size,fmt=:png)
    poincare_forced!(p1,f, u0, p, period; tcycles=tcycles,ncycles=ncycles,msize=msize,col=col,xlims=xlims,ylims=ylims,plotops...)
end

function poincare_forced_zoom(f::Function,u0::Vector{Float64},p::Vector{Float64},period::Float64;
    npts=1000,maxiter=1000,size=(600,600),xlims=[-1,1],ylims=[-1,1],plotops...)
    # Assume that the 3rd variable (z) is periodic and plot x,y as a function of z modulo period
    # skip trans period of the forcing
    tcycles = 30
    trans = solve(ODEProblem(f,u0,(0.0,tcycles*period),p))
    u0=trans.u[end]
    p1 = plot(;xlabel="x",ylabel="y",legend=false,size=size,fmt=:png,plotops...)
    ncycles = 1000
    tmax=ncycles*period
    kpts = 0
    kiter = 0
    while (kpts<npts) & (kiter<maxiter)
        tsave = period:period:ncycles*period
        sol = solve(ODEProblem(f,u0,(0.0,tmax),p),saveat=tsave)
        xlist = getindex.(sol.u,1)
        ylist = getindex.(sol.u,2)
        idx=inbox.(xlist,ylist,Ref(xlims),Ref(ylims))
        scatter!(getindex.(sol.u[idx],1),getindex.(sol.u[idx],2),color=:blue,markersize=0.5,xlims=xlims,ylims=ylims)
        u0 = sol.u[end]
        kpts += sum(idx)
        kiter += 1
    end
    println(kpts)
    println(kiter)
    p1
end

function recurrence_plot(f::Function,u0::Vector{Float64},p::Vector{Float64},period::Float64;
    dd=0.002,steps=10,tcycles=0,npts=300,ncycles=10,size=(900,450),plotops...)

    trans = solve(ODEProblem(f,u0,(0.0,tcycles*period),p))
    u0=trans.u[end]
    if period>1000
        tmax=min(ncycles*period,10000)
    else
        tmax=ncycles*period
    end
    sol = solve(ODEProblem(f,u0,(0,tmax),p));
    ts = range(0,tmax,length=npts)
    xy=hcat(sol(ts,idxs=1).u,sol(ts,idxs=2).u,mod.(sol(ts,idxs=2).u,period))
    dist = pairwise(Euclidean(),xy,dims=1)
    dst = floor.(dist/dd)/steps
    dst[dst.>steps] .= steps
    p1 = plot(sol,vars=(1,2),xlabel="x",ylabel="y")
    p2 = heatmap(ts,ts,dst,c=cgrad([:blue,:white]),legend=false,size=size)
    plot(p1,p2,layout=(1,2);size=size,fmt=:png,plotops...)
end

function saddle_orbit2D(f::Function,u0::Vector{Float64},p::Vector{Float64},period::Float64;
    λ=0.001,maxiter=10000, disttol=1e-9, inftol=10)
    # funcion que mapea un ciclo (extrender a N ciclos)
    function map_period(f,u0,p,period)
        sol = solve(ODEProblem(f,[u0[1],u0[2],0.0],(0,period),p),Tsit5(),reltol=1e-12,save_everystep=false)
        return sol.u[end][1:2]
    end
    # Funcion que aplica el metodo Schmelcher & Diakonos
    function Sk(f,prus,p,period,Λ)
        us = map_period(f,prus,p,period)
        return prus, prus+Λ*(us-prus)
    end
    # base de matrices del metodo para 2D
    Λ_base = λ*[[1 0; 0 1],[1 0; 0 -1],[-1 0; 0 1]]
    converged = false
    fp = u0
    # iteracion sobre los tres Λ base
    for Λ in Λ_base
        prus = u0
        # iteracion del metodo
        for i in 1:maxiter
            prus, us = Sk(f,prus, p,period,Λ)
            norm(us) > inftol && break
            if norm(prus - us) < disttol
                fp = us
                converged = true
                break
            end
            prus = us
        end
    end
    return fp,converged
end

function saddle_manifolds_forced(f::Function,f_jac::Function,us::Vector{Float64},p::Vector{Float64},period::Float64;
    ncycles=[10,3],npts=300,delta=0.01,xlims=(-1.0,1.0),ylims=(-1.0,1.0),size=(900,600),plotops...)
    # Asume que us es la orbita saddle en el plano
    xrange = xlims[2]-xlims[1]
    yrange = ylims[2]-ylims[1]
    condition(u,t,integrator) = (u[1]*u[1]+u[2]*u[2]) > max(xrange*xrange,yrange*yrange)
    affect!(integrator) = terminate!(integrator)
    t0_array = 0:period/npts:period
    sol = solve(ODEProblem(f,[us[1],us[2],0.0],(0,period),p),Tsit5(),reltol=1e-12,saveat=t0_array)
    u0_array=sol.u
    p1 = plot(;xlabel="x",ylabel="y",legend=false,size=size,fmt=:png,plotops...)
    for (u0,t0) in zip(u0_array,t0_array)
        A = f_jac(u0,p)
        av = eigen(A)
        for n=1:3
            u1 = u0+delta*av.vectors[:,n]
            u2 = u0-delta*av.vectors[:,n]
            if real(av.values[n])>0
                tmax=ncycles[1]*period
                tsave = period:period:ncycles[1]*period
                sol = solve(ODEProblem(f,u1,(t0,tmax),p),callback=DiscreteCallback(condition,affect!),saveat=tsave)
                scatter!(p1,sol,vars=(1,2),c=:red,markerstrokewidth=0,markersize=1)
                sol = solve(ODEProblem(f,u2,(t0,tmax),p),callback=DiscreteCallback(condition,affect!),saveat=tsave)
                scatter!(p1,sol,vars=(1,2),c=:red,markerstrokewidth=0,markersize=1)
            elseif real(av.values[n])<0
                tmax=ncycles[2]*period
                tsave = 0:-period:-ncycles[2]*period
                sol = solve(ODEProblem(f,u1,(t0,-tmax),p),callback=DiscreteCallback(condition,affect!),saveat=tsave)
                scatter!(p1,sol,vars=(1,2),c=:blue,markerstrokewidth=0,markersize=1)
                sol = solve(ODEProblem(f,u2,(t0,-tmax),p),callback=DiscreteCallback(condition,affect!),saveat=tsave)
                scatter!(p1,sol,vars=(1,2),c=:blue,markerstrokewidth=0,markersize=1)
            end
        end
    end
    xlims!(p1,xlims)
    ylims!(p1,ylims)
end

function butterfly(f::Function,u0::Vector{Float64},p::Vector{Float64};
    delta=1e-12,size=(900,400),xlims=false,ylims=false,plotops...)

    u1=u0
    u1[1]=u[1]+delta
    sol0 = solve(ODEProblem(f,u0,(0.0,tmax),p))
    sol1 = solve(ODEProblem(f,u1,(0.0,tmax),p))
    if length(u0) == 2
        p1 = plot(sol0,vars=(1,2),c=:red)
        plot!(p1,sol1,vars=(1,2),c=:blue)
    elseif length(u0) == 3
        p1 = plot(sol0,vars=(1,2,3),c=:red)
        plot!(p1,sol1,vars=(1,2,3),c=:blue)
    end
    dd =
    plot(p1,p2,layout=(1,2);size=size,plotops...)
end

function butterfly(f::Function,u0::Vector{Float64},p::Vector{Float64},tmax::Float64;
    dim=3,dt=0.001,delta=1e-12,size=(900,400),xlims=false,ylims=false,plotops...)

    u1=copy(u0)
    u1[1]=u1[1]+delta
    ts = range(0.0,stop=tmax,step=dt)
    sol0 = solve(ODEProblem(f,u0,(0.0,tmax),p))
    sol1 = solve(ODEProblem(f,u1,(0.0,tmax),p))
    if dim == 2
        p1 = plot(sol0,vars=(1,2),c=:red,xlabel="x",ylabel="y",label="u0")
        plot!(p1,sol1,vars=(1,2),c=:blue,label="u1")
    elseif dim == 3
        p1 = plot(sol0,vars=(1,2,3),c=:red,xlabel="x",ylabel="y",label="u0")
        plot!(p1,sol1,vars=(1,2,3),c=:blue,label="u1")
    end
    x0 = sol0(ts,idxs=1)
    x1 = sol1(ts,idxs=1)
    p2 = plot(ts,log10.(norm.(x0-x1)),xlabel="t",ylabel="log10(d)",label="")
    plot(p1,p2,layout=(1,2);size=size,plotops...)
end


text:
"""
    Zero-moment bulk microphysics scheme that instantly removes
    moisture above certain threshold.
    This is equivalent to instanteneous conversion of cloud condensate
    into precipitation and precipitation fallout with infinite
    terminal velocity.

"""
module Microphysics_0M

import Thermodynamics
import CLIMAParameters

const TD = Thermodynamics
const CP = CLIMAParameters
const CP_micro = CLIMAParameters.Atmos.Microphysics_0M
const APS = CP.AbstractParameterSet

export remove_precipitation

"""
    remove_precipitation(param_set::APS, q; q_vap_sat)

 - `param_set` - abstract parameter set
 - `q` - current PhasePartition
 - `q_vap_sat` - water vapor specific humidity at saturation

Returns the `q_tot` tendency due to the removal of precipitation.
The tendency is obtained assuming a relaxation with a constant timescale
to a state with precipitable water removed.
The threshold for when to remove `q_tot` is defined either by the
condensate specific humidity or supersaturation.
The thresholds and the relaxation timescale are defined in
CLIMAParameters.
"""
function remove_precipitation(
    param_set::APS,
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    _τ_precip::FT = CP_micro.τ_precip(param_set)
    _qc_0::FT = CP_micro.qc_0(param_set)

    return -max(0, (q.liq + q.ice - _qc_0)) / _τ_precip
end
function remove_precipitation(
    param_set::APS,
    q::TD.PhasePartition{FT},
    q_vap_sat::FT,
) where {FT <: Real}

    _τ_precip::FT = CP_micro.τ_precip(param_set)
    _S_0::FT = CP_micro.S_0(param_set)

    return -max(0, (q.liq + q.ice - _S_0 * q_vap_sat)) / _τ_precip
end

end #module Microphysics_0M.jl


text:
using StochasticProcesses

bm = BrownianMotion()
t = linspace(0, 1, 10000)
cumsim(bm, t, 100)[end,:]

println("cumsim(BrownianMotion) 10000x10000:")
for i in 1:5
  @time cumsim(bm, t, 10000)[end,:]
end


text:
"""
`ChaoticEncryption.jl` comes loaded with Pseudo-Random Number Generators and various
encryption techniques, which can be used to encrypt and decrypt any image file. The
package is under active development, but the existing API is stable and might not
change significantly.

For complete documentation, refer - https://saransh-cpp.github.io/ChaoticEncryption.jl
"""
module ChaoticEncryption

include("lorenz_key.jl")
include("logistic_key.jl")
include("substitution.jl")

export lorenz_key, logistic_key

export substitution_encryption,
       substitution_encryption!,
       substitution_decryption,
       substitution_decryption!

end


text:
### routines that deal with signal boundaries or windowing.

# frames in col-major-ordering.
# u contains the time stamps.
function extractframecenters(u, N_samples_per_frame, N_frames)

    fc = Vector{Float64}(undef,undef,N_frames)
    for j = 1:N_frames

        st = (j-1)*N_samples_per_frame +1
        fin = st+N_samples_per_frame - 1

        fc[j] = Statistics.mean(u[st:fin])
    end

    return fc
end

function extractframecenters(frames_time::Matrix{T})::Vector{T} where T

    N_frames = size(frames_time,2)

    fc = Vector{T}(undef,N_frames)
    for j = 1:N_frames
        fc[j] = Statistics.mean(frames_time[:,j])
    end

    return fc
end

# pad sequence y with zeros if length(y) < N
function padzeros(y::Vector{T},N)::Vector{T} where T
    if length(y) < N
        return [y; zeros(T,N-length(y))]
    end

    return y
end

function evalzeropaddedsignal(y::Vector{T}, n::Int)::T where T
    if 1 <= n <= length(y)
        return y[n]
    end

    return zero(T)
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
TOP = 1
RIGHT = 2
BOTTOM = 3
LEFT = 4

struct Tile
    id :: Int
    rows :: Array{String, 1}
    sides :: Array{String, 1}
end

function sides_of_rows(rows :: Array{String, 1}) :: Array{String, 1}
    sides = Array{Char,1}[]

    for _ in 1:4
        push!(sides, Array{Char, 1}())
    end

    for i in 1:10
        push!(sides[TOP], rows[1][i])
        push!(sides[RIGHT], rows[i][10])
        push!(sides[BOTTOM], rows[10][i])
        push!(sides[LEFT], rows[i][1])
    end

    result = String[]
    for i in 1:4
        push!(result, String(sides[i]))
    end
    return result
end

function rotate_lines(lines :: Array{String, 1}) :: Array{String, 1}
    h = length(lines)
    w = length(lines[1])

    result = String[]
    for i in 1:w
        row = Char[]
        for j in 1:h
            push!(row, lines[h - j + 1][i])
        end
        push!(result, String(row))
    end
    return result
end

function rotate(tile :: Tile) :: Tile
    rows = rotate_lines(tile.rows)
    # rows = String[]
    # for i in 1:10
    #     row = Char[]
    #     for j in 1:10
    #         push!(row, tile.rows[10 - j + 1][i])
    #     end
    #     push!(rows, String(row))
    # end

    return Tile(tile.id, rows, sides_of_rows(rows))
end

# 1 2 3    3 2 1
# 4 5 6 => 6 5 4
# 7 8 9    9 8 7
#
# 7 8 9
# 4 5 6
# 1 2 3

function flip_lines(lines :: Array{String, 1}) :: Array{String, 1}
    rows = String[]
    h = length(lines)
    for i in 1:h
        push!(rows, lines[h - i + 1])
    end
    return rows
end

function flip(tile :: Tile) :: Tile
    # rows = String[]
    # for i in 1:10
    #     push!(rows, tile.rows[10 - i + 1])
    # end
    rows = flip_lines(tile.rows)
    return Tile(tile.id, rows, sides_of_rows(rows))
end

function parse_file(file_path :: String) :: Array{Tile, 1}
    tiles = Tile[]
    open(file_path) do file
        parsing_tile = true

        id = 0
        rows = String[]
        for line in eachline(file)
            if parsing_tile
                id = parse(Int, split(split(line, " ")[2], ":")[1])
                parsing_tile = false
            else
                if length(line) > 0
                    push!(rows, line)
                else
                    push!(tiles, Tile(id, rows, sides_of_rows(rows)))
                    rows = String[]
                    parsing_tile = true
                end
            end
        end
        push!(tiles, Tile(id, rows, sides_of_rows(rows)))
    end
    return tiles
end

function find_right_bottom(left :: Tile, top :: Tile, tiles :: Array{Tile, 1}) :: Tile
    for tile in tiles
        if top.id != tile.id && left.id != tile.id
            result = tile

            for i in 1:4
                if top.sides[BOTTOM] == result.sides[TOP] && left.sides[RIGHT] == result.sides[LEFT]
                    return result
                end
                result = rotate(result)
            end

            result = flip(tile)

            for i in 1:4
                if top.sides[BOTTOM] == result.sides[TOP] && left.sides[RIGHT] == result.sides[LEFT]
                    return result
                end
                result = rotate(result)
            end
        end
    end

    @assert(false, "Could not find solution for find_right_bottom")
end

function find_bottom(top :: Tile, tiles :: Array{Tile, 1}) :: Tile
    for tile in tiles
        if top.id != tile.id
            bottom = tile

            for i in 1:4
                if top.sides[BOTTOM] == bottom.sides[TOP]
                    return bottom
                end
                bottom = rotate(bottom)
            end

            bottom = flip(tile)

            for i in 1:4
                if top.sides[BOTTOM] == bottom.sides[TOP]
                    return bottom
                end
                bottom = rotate(bottom)
            end
        end
    end

    @assert(false, "Could not find solution for find_right")
end

function find_right(left :: Tile, tiles :: Array{Tile, 1}) :: Tile
    for tile in tiles
        if left.id != tile.id
            right = tile
            for i in 1:4
                if left.sides[RIGHT] == right.sides[LEFT]
                    return right
                end
                right = rotate(right)
            end

            right = flip(tile)

            for i in 1:4
                if left.sides[RIGHT] == right.sides[LEFT]
                    return right
                end
                right = rotate(right)
            end
        end
    end

    @assert(false, "Could not find solution for find_right")
end

function side_matches(candidate :: Tile, tiles :: Array{Tile, 1}) :: Array{Int, 1}
    cs = fill(0, (4))
    for tile in tiles
        if candidate.id != tile.id
            for left in 1:4
                for right in 1:4
                    if candidate.sides[left] == tile.sides[right]
                        cs[left] += 1
                    end
                    if candidate.sides[left] == reverse(tile.sides[right])
                        cs[left] += 1
                    end
                end
            end
        end
    end

    # SIDES:
    # TOP = 1
    # RIGHT = 2
    # BOTTOM = 3
    # LEFT = 4
    #############
    # CORNERS:
    # 0     1
    #
    # 3     2

    return cs
end

function is_corner(candidate :: Tile, tiles :: Array{Tile, 1}) :: Bool
    cs = side_matches(candidate, tiles)

    res = 0
    for c in cs
        if c == 0
            res += 1
        end
    end

    return res == 2
end

function part_1(tiles :: Array{Tile, 1}) :: Int
    result = 1
    for tile in tiles
        if is_corner(tile, tiles)
            result *= tile.id
        end
    end
    return result
end

function count_hash(lines :: Array{String, 1}) :: Int
    result = 0
    for line in lines
        for x in line
            if x == '#'
                result += 1
            end
        end
    end
    return result
end

monster = ["                  # ",
           "#    ##    ##    ###",
           " #  #  #  #  #  #   "]

function contains_monster(drow :: Int, dcol :: Int, picture :: Array{String, 1}) :: Bool
    for row0 in 0:(3 - 1)
        for col0 in 0:(20 - 1)
            if monster[row0 + 1][col0 + 1] == '#'
                if picture[row0 + drow][col0 + dcol] != '#'
                    return false
                end
            end
        end
    end
    return true
end

function count_monsters(picture :: Array{String, 1}) :: Int
    result = 0
    for drow in 1:(96 - (3))
        for dcol in 1:(96 - (20))
            if contains_monster(drow, dcol, picture)
                result += 1
            end
        end
    end
    return result
end

function part_2(tiles :: Array{Tile, 1}) :: Int
    puzzle = Array{Array{Tile, 1}, 1}()

    push!(puzzle, Array{Tile, 1}())
    for tile in tiles
        if side_matches(tile, tiles) == [0, 1, 1, 0]
            push!(puzzle[1], tile)
            break
        end
    end

    for i in 2:12
        right = find_right(puzzle[1][i - 1], tiles)
        push!(puzzle[1], right)
    end

    for i in 2:12
        push!(puzzle, Array{Tile, 1}())
        bottom = find_bottom(puzzle[i - 1][1], tiles)
        push!(puzzle[i], bottom)

        for j in 2:12
            right_bottom = find_right_bottom(puzzle[i][j - 1],
                                             puzzle[i - 1][j],
                                             tiles)
            push!(puzzle[i], right_bottom)
        end
    end

    picture = Array{String, 1}()

    for row in puzzle
        for i in 2:9
            picture_row = Array{String, 1}()
            for tile in row
                push!(picture_row, tile.rows[i][2:9])
            end
            push!(picture, join(picture_row))
        end
    end

    picture = flip_lines(picture)

    return count_hash(picture) - count_monsters(picture) * count_hash(monster)
end

function solve_file(file_path :: String)
    println("Input file: $(file_path)")
    tiles = parse_file(file_path)
    println("Part 1: $(part_1(tiles))")
    println("Part 2: $(part_2(tiles))")
end

for file_path in ARGS
    solve_file(file_path)
end


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements sample function
# 3. Exports Type
"""
    LinSysVecRowUnidSampler <: LinSysVecRowSampler

An immutable structure without fields that specifies randomly cycling from the rows of a
linear system with uniform probability and with replacement.
"""
struct LinSysVecRowUnidSampler <: LinSysVecRowSampler end

# Common sample interface for linear systems
function sample(
    type::LinSysVecRowUnidSampler,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)
    eqn_ind = rand(Base.OneTo(length(b)))

    return A[eqn_ind, :], b[eqn_ind]
end

#export LinSysVecRowUnidSampler


text:
using PolyChaos, Test, LinearAlgebra
import QuadGK: quadgk

function integration_test(opq::OrthoPolyQ,dim::Int64,name::String;tol::Float64=1e-6,
                          dom::Tuple{Float64,Float64}=opq.op.meas.dom)
    op::OrthoPoly = opq.op
    N = op.deg
    w::Function = op.meas.w
    ind = zeros(Int64,dim)
    @testset "$name" begin
        for ind_ in Iterators.product([collect(1:N) for i=1:dim]...)
            ind[:] .= ind_[:]
            s1 = computeSP(ind,opq)
            f(t) = prod(evaluate(ind[i],t,op) for i=1:dim)*w(t)
            # f(t) = evaluate(ind[1],t,op)*evaluate(ind[2],t,op)*evaluate(ind[3],t,op)*w(t)
            s2 = quadgk(f,dom[1],dom[2])[1]
            if abs(s1)<=1e-3
                @test isapprox(s1,s2;atol=tol)
            else
                @test isapprox(abs(s1/s2),1;atol=tol)
            end
        end
    end
end
##

N, Nrec, dim, tol = 3, 1000, 3, 1e-4

op = OrthoPoly("gaussian",N;Nrec=Nrec);
opq = OrthoPolyQ(op);
integration_test(opq,dim,"gaussian";tol=tol,dom=(-10.,10.))

α, β = 1.32, 3.23
op = OrthoPoly("beta01",N,Dict(:shape_a=>α,:shape_b=>β);Nrec=Nrec);
opq = OrthoPolyQ(op);
set_beta = integration_test(opq,dim,"beta01";tol=tol)

op = OrthoPoly("gamma",N,Dict(:shape=>1.0,:rate=>1.);Nrec=Nrec);
opq = OrthoPolyQ(op);
set_gamma = integration_test(opq,dim,"gamma";tol=tol,dom=(0.,Inf))


text:
expres_data_path = joinpath(homedir(),"Data/EXPRES/inputs")
solar_data_path = joinpath(homedir(),"Data/SolarSpectra/NEID_solar/")
ancilary_solar_data_path = joinpath(homedir(),"Data/SolarSpectra/NEID_solar/")
output_dir = "examples/output"


text:
# This contains algorithms for breaking up a DAG into layers

function layer_by_longest_path_to_source(graph)
    dists = longest_paths(graph, sources(graph))
    layer_groups = IterTools.groupby(i->dists[i], sort(vertices(graph), by=i->dists[i]))
    return collect.(layer_groups)
end

###################
# helpers

function node2layer_lookup(layer2nodes)
    flat_map = [node=>layer_ind for (layer_ind, nodes) in enumerate(layer2nodes) for node in nodes]
    sort!(flat_map)
    @assert first.(flat_map) == 1:length(flat_map)
    return last.(sort(flat_map))
end

"add nodes so that no edges span multiple layers, returns mask"
function add_dummy_nodes!(graph, layer2nodes)
    dag_or_error(graph)
    nondummy_nodes = vertices(graph)
    # mapping from edges in the original graph to paths in the graph with dummy nodes
    edge_to_paths = Dict(e => eltype(graph)[] for e in edges(graph))
    node2layer = node2layer_lookup(layer2nodes)  # doesn't have dummy nodes, doesn't need them
    for cur_node in vertices(graph)
        cur_layer = node2layer[cur_node]
        for out_edge in filter(e -> src(e) == cur_node, collect(edges(graph)))  # need to copy as outwise will mutate when the graph is mutated
            out_node = out_edge.dst
            out_layer = node2layer[out_node]
            cur_layer < out_layer || throw(DomainError(node2layer, "Layer assigmenment must be strictly monotonic"))
            path = get!(edge_to_paths, out_edge, eltype(graph)[])
            if out_layer != cur_layer + 1
                rem_edge!(graph, cur_node, out_node) || error("removing edge failed")
                prev_node = cur_node
                push!(path, prev_node)
                for step_layer in (cur_layer + 1) : (out_layer - 1)
                    add_vertex!(graph) || throw(OverflowError("Could not add vertices to graph."))
                    step_node = nv(graph)
                    push!(layer2nodes[step_layer], step_node)
                    add_edge!(graph, prev_node, step_node)
                    push!(path, step_node)
                    prev_node = step_node
                end
                add_edge!(graph, prev_node, out_node)
                push!(path, out_node)
            else
                push!(path, cur_node, out_node)
            end
        end
    end
    is_dummy_mask = trues(nv(graph))
    is_dummy_mask[nondummy_nodes] .= false
    return is_dummy_mask, edge_to_paths
end


text:
using Test
using WebIO
using Random

using WebIO: kebab2camel, camel2kebab

@testset "kebabs and camels" begin
    kebabstrs = ["vue-instance-17-node-18", "bum-two-three", "silly nanny", "nanny"]
    cameltrues = ["vueInstance17Node18", "bumTwoThree", "silly nanny", "nanny"]
    camelstrs = kebab2camel.(kebabstrs)
    for i in length(camelstrs)
        @test camelstrs[i] == cameltrues[i]
    end

    kebabstrs = camel2kebab.(cameltrues)
    kebab_also_trues = ["vue-instance17-node18", "bum-two-three", "silly nanny", "nanny"]
    for i in length(kebabstrs)
        @test camelstrs[i] == kebabstrs[i] || camelstrs[i] == kebab_also_trues[i]
    end
end


text:
module DistanceTransforms

using Images
using Distances
using Statistics

include("./helper.jl")
include("./metrics.jl")

export find_edges,
    mean_hausdorff,
    mean_dice

end


text:
function cumulative(weights::Array{Int, 2})
	ret = copy(weights)
	dim = size(weights)
	for i in 2 : dim[1]
		for j in 1 : dim[2]
			min = ret[i - 1, j]
			if j > 1 && ret[i - 1, j - 1] < min
				min = ret[i - 1, j - 1]
			end
			if j < dim[2] && ret[i - 1, j + 1] < min
				min = ret[i - 1, j + 1]
			end
			ret[i, j] += min
		end
	end
	return ret
end

function back_track(weights::Array{Int, 2})
	dim = size(weights)

	buffer::Array{Int, 1} = ones(Int, dim[1])

	buffer[1] = 1
	for i in 2 : dim[2]
		if weights[dim[1], i] < weights[dim[1], buffer[1]]
			buffer[1] = i
		end
	end

	for i in dim[1] - 1 : -1 : 1
		k = dim[1] - i + 1
		j = buffer[k - 1]
		buffer[k] = j
		if j > 1 && weights[i, j - 1] < weights[i, buffer[k]]
			buffer[k] = j - 1
		end
		if j < dim[2] && weights[i, j + 1] < weights[i, buffer[k]]
			buffer[k] = j + 1
		end
	end

	ret::Array{Tuple{Int, Int}, 1} = Array{Tuple{Int, Int}, 1}(undef, dim[1])

	for i in 1 : dim[1]
		ret[i] = dim[1] - i + 1, buffer[i]
	end

	return ret
end


text:
#=
A builder is looking to build a row of N houses that can be of K different colors. He has a goal of minimizing cost while ensuring that no two neighboring houses are of the same color.

Given an N by K matrix where the nth row and kth column represents the cost to build the nth house with kth color, return the minimum cost which achieves this goal.
=#
using Test
include("Solutions/problem19_minimum_building_cost.jl")

costs = [1 2 3; 2 10 10; 4 5 9; 10 1 10]

@test minimum_building_cost(costs) == 14.0

costs = [4 2 8 4; 1 2 5 0; 10 11 13 4; 1 1 4 9]

@test minimum_building_cost(costs) == 8.0

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using MuladdMacro
using Parameters: @with_kw
using StaticArrays: @SMatrix

import ..AU

import AstroTime: Epoch, TDBEpoch, j2000, centuries, value, SECONDS_PER_CENTURY
import ..Bodies:
    SolarSystemBarycenter,
    Sun,
    Mercury,
    Venus,
    Earth,
    Mars,
    Jupiter,
    Saturn,
    Uranus,
    Neptune

export VSOP87

include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_sun.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_mercury.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_venus.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_earth.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_mars.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_jupiter.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_saturn.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_uranus.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "vsop87_neptune.jl"))

const VSOP87_FK5 = @SMatrix [1.0 0.000000440360 -0.000000190919;
							 -0.000000479966 0.917482137087 -0.397776982902;
							 0.0 0.397776982902 0.917482137087]

@with_kw struct VSOP87 <: AbstractEphemeris
    order::Int = 5; @assert order >= 0; @assert order <= 5
    terms::Int = typemax(Int); @assert terms > 0
end

@inbounds function _vsop(coeffs, n_terms, max_order, max_terms, ep)
    t0 = value(centuries(j2000(TDBEpoch(ep)))) / 10.0

    # Powers of `t0`
    t = Array{Float64}(undef, 7)
    t[1] = 0.0 # t[-1]
    t[2] = 1.0 # t[0]
    t[3] = t0 # t[1]
    for i = 4:max_order + 2
        t[i] = t[i-1] * t0
    end

    r = zeros(3)
    v = zeros(3)

    for k = 1:3
        m = min(length(coeffs[k]) - 1, max_order)
        for order = 0:m
            i = order + 1
            it = order + 2
            n = min(n_terms[k][i], max_terms)
            for j = 1:n
                a = coeffs[k][i][1][j]
                b = coeffs[k][i][2][j]
                c = coeffs[k][i][3][j]

                u = b + c * t0
                su, cu = sincos(u)
                r[k] += a * cu * t[it]
                v[k] += order * a * cu * t[it-1] - a * c * su * t[it]
            end
        end
    end

    # AU to km
    r .*= AU
    # AU/millenium to km/s
    v .*= AU / 10SECONDS_PER_CENTURY

    VSOP87_FK5 * r, VSOP87_FK5 * v
end

for body in ("Sun", "Mercury", "Venus", "Earth", "Mars",
             "Jupiter", "Saturn", "Uranus", "Neptune")
    b = Symbol(body)
    c = Symbol("VSOP_", uppercase(body))
    n = Symbol("VSOP_", uppercase(body), "_NUM")
    @eval begin
        function position!(arr, eph::VSOP87, ep::Epoch, ::SolarSystemBarycenter, ::$b)
            _vsop($c, $n, eph.order, eph.terms, ep)[1]
        end

        function velocity!(arr, eph::VSOP87, ep::Epoch, ::SolarSystemBarycenter, ::$b)
            _vsop($c, $n, eph.order, eph.terms, ep)[2]
        end

        function state!(r0, v0, eph::VSOP87, ep::Epoch, ::SolarSystemBarycenter, ::$b)
            r1, v1 = _vsop($c, $n, eph.order, eph.terms, ep)
            r0 .+= r1
            v0 .+= v1
            r0, v0
        end
    end
end

function state!(pos, vel, vsop::VSOP87, ep::Epoch, from::SolarSystemBarycenter, to::CelestialBody)
    throw(ArgumentError("Body '$to' is not supported by VSOP87."))
end

function state!(pos, vel, vsop::VSOP87, ep::Epoch, from::CelestialBody, to::CelestialBody)
    rv1 = state!(zeros(3), zeros(3), vsop, ep, ssb, from)
    rv2 = state!(zeros(3), zeros(3), vsop, ep, ssb, to)
    pos .+= (rv2[1] .- rv1[1])
    vel .+= (rv2[2] .- rv1[2])
    pos, vel
end

function position!(pos, vsop::VSOP87, ep::Epoch, from::CelestialBody, to::CelestialBody)
    state!(pos, zeros(3), vsop, ep, from, to)
    pos
end

function velocity!(vel, vsop::VSOP87, ep::Epoch, from::CelestialBody, to::CelestialBody)
    state!(zeros(3), vel, vsop, ep, from, to)
    vel
end


text:
#= prepare to test =#

typealias AD15  ArbFloat{50}    # 15 digits
typealias AD30  ArbFloat{100}   # 31 digits
typealias AD100 ArbFloat{333}   # 100 digits
typealias AD300 ArbFloat{997}   # 300 digits

af0 = zero(ArbFloat)
af1 = one(ArbFloat)
af2 = two(ArbFloat)
af3 = three(ArbFloat)
af4 = four(ArbFloat)
afhalf = inv(af2)

#= prepared =#


text:
import Pkg
Pkg.add("Documenter")
using Documenter
using JuliaNotes
push!(LOAD_PATH,"../src/")
makedocs(
    modules=[JuliaNotes],
    sitename="JuliaNotes.jl",
    pages = [
        "Home" => "index.md",
        "Development workflow" => "workflow.md",
        "Modules and Revise" => "modules.md",
        "Scope of loops" => "loopscopes.md",
        "Benchmark" => "benchmark.md",
        "Assignment and mutation" => "assignment.md",
        "Type instability" => "instability.md",
        "Closures" => "closures.md",
        #"Union splitting" => "splitting.md",
        "Tracking allocations" => "memory.md",
        "Immutable variables" => "immutable.md",
        "ERROR: No method..." => "nomethod.md",
        "Vector{Int} <: Vector{Real} is false???" => "typevariance.md",
        "Figures and LaTeX" => "figures.md",
        "Create new package" => "new_package.md",
        "Publish Docs" => "publish_docs.md",
    ]
)
deploydocs(
    repo = "github.com/m3g/JuliaNotes.jl.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#" ],
)


text:
using TimeZones

TimeZones.TZData.compile(max_year=2200)

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
@testset "simple_varinfo.jl" begin
    @testset "constructor & indexing" begin
        @testset "NamedTuple" begin
            svi = SimpleVarInfo(; m=1.0)
            @test getlogp(svi) == 0.0
            @test haskey(svi, @varname(m))
            @test !haskey(svi, @varname(m[1]))

            svi = SimpleVarInfo(; m=[1.0])
            @test getlogp(svi) == 0.0
            @test haskey(svi, @varname(m))
            @test haskey(svi, @varname(m[1]))
            @test !haskey(svi, @varname(m[2]))
            @test svi[@varname(m)][1] == svi[@varname(m[1])]

            svi = SimpleVarInfo(; m=(a=[1.0],))
            @test haskey(svi, @varname(m))
            @test haskey(svi, @varname(m.a))
            @test haskey(svi, @varname(m.a[1]))
            @test !haskey(svi, @varname(m.a[2]))
            @test !haskey(svi, @varname(m.a.b))

            svi = SimpleVarInfo{Float32}(; m=1.0)
            @test getlogp(svi) isa Float32

            svi = SimpleVarInfo((m=1.0,), 1.0)
            @test getlogp(svi) == 1.0
        end

        @testset "Dict" begin
            svi = SimpleVarInfo(Dict(@varname(m) => 1.0))
            @test getlogp(svi) == 0.0
            @test haskey(svi, @varname(m))
            @test !haskey(svi, @varname(m[1]))

            svi = SimpleVarInfo(Dict(@varname(m) => [1.0]))
            @test getlogp(svi) == 0.0
            @test haskey(svi, @varname(m))
            @test haskey(svi, @varname(m[1]))
            @test !haskey(svi, @varname(m[2]))
            @test svi[@varname(m)][1] == svi[@varname(m[1])]

            svi = SimpleVarInfo(Dict(@varname(m) => (a=[1.0],)))
            @test haskey(svi, @varname(m))
            @test haskey(svi, @varname(m.a))
            @test haskey(svi, @varname(m.a[1]))
            @test !haskey(svi, @varname(m.a[2]))
            @test !haskey(svi, @varname(m.a.b))

            svi = SimpleVarInfo(Dict(@varname(m.a) => [1.0]))
            # Now we only have a variable `m.a` which is subsumed by `m`,
            # but we can't guarantee that we have the "entire" `m`.
            @test !haskey(svi, @varname(m))
            @test haskey(svi, @varname(m.a))
            @test haskey(svi, @varname(m.a[1]))
            @test !haskey(svi, @varname(m.a[2]))
            @test !haskey(svi, @varname(m.a.b))
        end
    end

    @testset "SimpleVarInfo on $(nameof(model))" for model in
                                                     DynamicPPL.TestUtils.DEMO_MODELS
        # We might need to pre-allocate for the variable `m`, so we need
        # to see whether this is the case.
        m = model().m
        svi_nt = if m isa AbstractArray
            SimpleVarInfo((m=similar(m),))
        else
            SimpleVarInfo()
        end
        svi_dict = SimpleVarInfo(VarInfo(model), Dict)

        @testset "$(nameof(typeof(svi.values)))" for svi in (svi_nt, svi_dict)
            # Random seed is set in each `@testset`, so we need to sample
            # a new realization for `m` here.
            m = model().m

            ### Sampling ###
            # Sample a new varinfo!
            _, svi_new = DynamicPPL.evaluate!!(model, svi, SamplingContext())

            # If the `m[1]` varname doesn't exist, this is a univariate model.
            # TODO: Find a better way of dealing with this that is not dependent
            # on knowledge of internals of `model`.
            isunivariate = !haskey(svi_new, @varname(m[1]))

            # Realization for `m` should be different wp. 1.
            if isunivariate
                @test svi_new[@varname(m)] != m
            else
                @test svi_new[@varname(m[1])] != m[1]
                @test svi_new[@varname(m[2])] != m[2]
            end
            # Logjoint should be non-zero wp. 1.
            @test getlogp(svi_new) != 0

            ### Evaluation ###
            # Sample some random testing values.
            m_eval = if m isa AbstractArray
                randn!(similar(m))
            else
                randn(eltype(m))
            end

            # Update the realizations in `svi_new`.
            svi_eval = if isunivariate
                DynamicPPL.setindex!!(svi_new, m_eval, @varname(m))
            else
                DynamicPPL.setindex!!(svi_new, m_eval, [@varname(m[1]), @varname(m[2])])
            end
            # Reset the logp field.
            svi_eval = DynamicPPL.resetlogp!!(svi_eval)

            # Compute `logjoint` using the varinfo.
            logπ = logjoint(model, svi_eval)
            # Extract the parameters from `svi_eval`.
            m_vi = if isunivariate
                svi_eval[@varname(m)]
            else
                svi_eval[[@varname(m[1]), @varname(m[2])]]
            end
            # These should not have changed.
            @test m_vi == m_eval
            # Compute the true `logjoint` and compare.
            logπ_true = DynamicPPL.TestUtils.logjoint_true(model, m_vi)
            @test logπ ≈ logπ_true
        end
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
using Yao.EasyBuild, Test
@testset "general U2 U4" begin
    c = general_U2(0.5, 0.7, 0.9; ϕ=0.8)
    @test length(c) == 4
    c = general_U2(0.5, 0.7, 0.9)
    @test length(c) == 3

    params = rand(17)*2π
    @test_throws ArgumentError general_U4(params)
    params = rand(15)*2π
    g = general_U4(params)
    gc = gatecount(g)
    @test gc[typeof(Ry(0.0))] == 6
    @test gc[typeof(Rz(0.0))] == 9
    @test gc[typeof(cnot(2,2,1))] == 3
    @test parameters(g) ≈ params
end


text:
using Test
using QuantumLattices.Mathematics.VectorSpaces
using QuantumLattices.Mathematics.Combinatorics: Combinations
using QuantumLattices.Interfaces: dimension, ⊕, rank, dims, inds

@testset "SimpleVectorSpace" begin
    id1, id2, id3 = (1, 1), (1, 2), (1, 3)

    vs = SimpleVectorSpace{'T'}(id1, id2)
    @test vs == deepcopy(vs)
    @test isequal(vs, deepcopy(vs))
    @test vs|>size == (2,)
    @test vs|>dimension == 2
    @test vs|>collect == [id1, id2]
    @test (vs[1] == id1) && (vs[2] == id2)
    @test (searchsortedfirst(vs, id1) == 1) && (searchsortedfirst(vs, id2) == 2)
    @test (findfirst(id1, vs) == 1) && (findfirst(id2, vs) == 2)
    @test findfirst((id1, id2), vs) == (1, 2)
    @test (id1 ∈ vs) && (id2 ∈ vs) && (id3 ∉ vs)
    @test HasTable(typeof(vs)) == HasTable(true)
    @test TableSorted(typeof(vs)) == TableSorted(true)
    @test IsMultiIndexable(typeof(vs)) == IsMultiIndexable(false)
    @test vs == id1 ⊕ id2
    @test vs ⊕ id3 == id1 ⊕ id2 ⊕ id3
    @test id3 ⊕ vs == id3 ⊕ id1 ⊕ id2
    @test (id2 ⊕ id3) ⊕ vs == id2 ⊕ id3 ⊕ id1 ⊕ id2

    vs = SimpleVectorSpace{'F'}(id1, id2)
    @test vs == deepcopy(vs)
    @test isequal(vs, deepcopy(vs))
    @test vs|>size == (2,)
    @test vs|>dimension == 2
    @test vs|>collect == [id1, id2]
    @test (vs[1] == id1) && (vs[2] == id2)
    @test (searchsortedfirst(vs, id1) == 1) && (searchsortedfirst(vs, id2) == 2)
    @test (findfirst(id1, vs) == 1) && (findfirst(id2, vs) == 2)
    @test findfirst((id1, id2), vs) == (1, 2)
    @test (id1 ∈ vs) && (id2 ∈ vs) && (id3 ∉ vs)
    @test HasTable(typeof(vs)) == HasTable(true)
    @test TableSorted(typeof(vs)) == TableSorted(false)
    @test IsMultiIndexable(typeof(vs)) == IsMultiIndexable(false)
    @test vs == id1 ⊕ id2
    @test vs ⊕ id3 == id1 ⊕ id2 ⊕ id3
    @test id3 ⊕ vs == id3 ⊕ id1 ⊕ id2
    @test (id2 ⊕ id3) ⊕ vs == id2 ⊕ id3 ⊕ id1 ⊕ id2
end

@testset "SimpleIndices" begin
    foi = SimpleIndices{'F'}(2, 2, 2)
    @test HasTable(typeof(foi)) == HasTable(false)
    @test IsMultiIndexable(typeof(foi)) == IsMultiIndexable(true)
    @test MultiIndexOrderStyle(typeof(foi)) == MultiIndexOrderStyle('F')
    @test dimension(foi) == 8
    @test dims(foi) == (2, 2, 2)
    @test rank(typeof(foi)) == 3
    @test inds((1, 1, 1), foi) == (1, 1, 1)
    @test Tuple((1, 1, 1), foi) == (1, 1, 1)
    @test foi|>collect == [(1, 1, 1), (2, 1, 1), (1, 2, 1), (2, 2, 1), (1, 1, 2), (2, 1, 2), (1, 2, 2), (2, 2, 2)]
    @test ((1, 1, 1) ∈ foi) && ((1, 2, 3) ∉ foi)
    for (i, finds) in enumerate(foi)
        @test findfirst(finds, foi) == i
        @test searchsortedfirst(foi, foi[i]) == i
    end
    coi = SimpleIndices{'C'}(2, 2, 2)
    @test HasTable(typeof(coi)) == HasTable(false)
    @test IsMultiIndexable(typeof(coi)) == IsMultiIndexable(true)
    @test MultiIndexOrderStyle(typeof(coi)) == MultiIndexOrderStyle('C')
    @test dimension(coi) == 8
    @test dims(coi) == (2, 2, 2)
    @test rank(typeof(coi)) == 3
    @test inds((1, 1, 1), coi) == (1, 1, 1)
    @test Tuple((1, 1, 1), coi) == (1, 1, 1)
    @test ((1, 1, 1) ∈ coi) && ((1, 2, 3) ∉ coi)
    @test coi|>collect == [(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2), (2, 1, 1), (2, 1, 2), (2, 2, 1), (2, 2, 2)]
    for (i, cinds) in enumerate(coi)
        @test findfirst(cinds, coi) == i
        @test searchsortedfirst(coi, coi[i]) == i
    end
end

@testset "TabledIndices" begin
    dims = (2, 2)
    table = [(1, 1), (1, 2), (2, 1), (2, 2)]
    toi = TabledIndices{'T'}(dims, table)
    @test HasTable(typeof(toi)) == HasTable(true)
    @test TableSorted(typeof(toi)) == TableSorted(true)
    @test toi == TabledIndices{2}(DulPermutations, 2)
    @test dimension(toi) == 4
    @test toi|>collect == table
    @test ((1, 1) ∈ toi) && ((1, 3) ∉ toi)
    for i = 1:length(toi)
        @test searchsortedfirst(toi, toi[i]) == i
    end

    table = [(1, 2), (2, 1), (2, 2), (1, 1)]
    toi = TabledIndices{'F'}(dims, table)
    @test HasTable(typeof(toi)) == HasTable(true)
    @test TableSorted(typeof(toi)) == TableSorted(false)
    @test dimension(toi) == 4
    @test toi|>collect == table
    @test ((1, 1) ∈ toi) && ((1, 3) ∉ toi)
    for i = 1:length(toi)
        @test searchsortedfirst(toi, toi[i]) == i
    end
end

struct VSZNamedVectorSpace{NS, BS<:Tuple, VS<:Tuple{Vararg{Vector}}} <: NamedVectorSpace{:zip, NS, BS, VS}
    contents::VS
end
@generated function VSZNamedVectorSpace{NS}(contents::Vector...) where NS
    @assert (length(NS) == length(contents)) && isa(NS, Tuple{Vararg{Symbol}})
    BS = Expr(:curly, :Tuple, [contents[i]|>eltype for i = 1:length(NS)]...)
    return quote
        @assert mapreduce(length, ==, contents)
        VSZNamedVectorSpace{NS, $BS, typeof(contents)}(contents)
    end
end

struct VSPNamedVectorSpace{NS, BS<:Tuple, VS<:Tuple{Vararg{Vector}}} <: NamedVectorSpace{:product, NS, BS, VS}
    contents::VS
end
@generated function VSPNamedVectorSpace{NS}(contents::Vector...) where NS
    @assert (length(NS) == length(contents)) && isa(NS, Tuple{Vararg{Symbol}})
    BS = Expr(:curly, :Tuple, [contents[i]|>eltype for i = 1:length(NS)]...)
    return :(VSPNamedVectorSpace{NS, $BS, typeof(contents)}(contents))
end

@testset "NamedVectorSpace" begin
    @test IsMultiIndexable(NamedVectorSpace) == IsMultiIndexable(true)
    @test MultiIndexOrderStyle(NamedVectorSpace) == MultiIndexOrderStyle('C')

    nvs = VSZNamedVectorSpace{(:t, :U)}([1, 2], [8.0, 9.0])
    @test nvs|>keys == nvs|>typeof|>keys == (:t, :U)
    @test nvs|>values == ([1, 2], [8.0, 9.0])
    @test nvs|>pairs|>collect == [:t=>[1, 2], :U=>[8.0, 9.0]]
    @test eltype(nvs, 1) == eltype(nvs|>typeof, 1) == Int
    @test eltype(nvs, 2) == eltype(nvs|>typeof, 2) == Float64
    @test nvs|>typeof|>rank == 1
    @test dims(nvs) == (2,)
    elements = [(t = 1, U = 8.0), (t = 2, U = 9.0)]
    for i = 1:dimension(nvs)
        @test NamedTuple(inds(elements[i], nvs), nvs) == elements[i]
    end
    @test nvs|>collect == elements

    nvs = VSPNamedVectorSpace{(:t, :U)}([1.0, 2.0], [8.0, 9.0])
    @test nvs|>typeof|>rank == 2
    @test dims(nvs) == (2, 2)
    elements = [(t = 1.0, U = 8.0), (t = 1.0, U = 9.0), (t = 2.0, U = 8.0), (t = 2.0, U = 9.0)]
    for i = 1:dimension(nvs)
        @test NamedTuple(inds(elements[i], nvs), nvs) == elements[i]
    end
    @test nvs|>collect == elements
end


text:
#!/usr/local/bin/julia

using MixedModels
using RDatasets: dataset

"""
# InstEval

Data are from a special evaluation of lecturers by students at the Swiss Federal
Institute for Technology–Zürich (ETH–Zürich). There are nearly 75,000 evaluations
by 2,972 students on a total of 1,128 instructors.
"""
InstEval = dataset("lme4", "InstEval")

DataFrames.printtable(InstEval)

# Define a formula
f = @formula(Y ~ 1 + Dept*Service + (1 | S) + (1 | D))
typeof(f)

# Create a LinearMixedModel
m = lmm(f, InstEval)
typeof(m)

# Fit a statistical model in-place
fittedmodel = fit!(m)
typeof(fittedmodel)

print(fittedmodel)

max_stars_repo_name: ellisvalentiner/julia-talk

max_issues_repo_name: ellisvalentiner/julia-talk

max_forks_repo_name: ellisvalentiner/julia-talk


text:
function detect(obj::CycleSlipDetectors, arg0::List)
    return jcall(obj, "detect", List, (List,), arg0)
end


text:
using Documenter
using MendelKinship

makedocs(
    sitename = "MendelKinship",
    format = Documenter.HTML(),
    modules = [MendelKinship],
    pages = [
        "Home" => "index.md",
        "Introduction" => "man/introduction.md",
        "Tutorial" => "man/KinshipTutorial.md",
        "Keyword Options" => "man/keywords.md",
        "API" => "man/api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/OpenMendel/MendelKinship.jl.git",
    target = "build"
)


text:
using GDB
using Test
using SparseArrays
using LinearAlgebra
using Plots

function test_domain()

    bar1 = GDB.Barrier() do
        func(x,y) = x
        val = 1.8
        rmap(x,y) = 1 - x, y
        func, val, rmap, -1
    end

    bar2 = GDB.PolyBarrier() do
        Float64[1 1; -1 2]
    end

    bar3 = GDB.Everywhere()
    bar4 = GDB.Nowhere()

    grd = GDB.Grid(50, 50, 1) do
        r0 = Float64[0, 0]
        r1 = Float64[2, 1]
        inside(x,y) = GDB.check_if_inside(x, y, [0.1], [bar1])
        inside, r0, r1
    end

    fvec = GDB.f_to_grid(grd) do x,y
        sin(x)*sin(y)
    end

    h = heatmap(fvec, grd, 1)

    f_func = GDB.interpolation_function(fvec, 2, 2, GDB.stencil2d(2, 2), grd)

    gc = GDB.GhostConditions([bar1], grd)
    gc2 = GDB.swap_sign(gc, 1)

    A = sparse(I, grd.Nk, grd.Nk)
    xb = zeros(grd.Nk)
    x1 = rand(grd.Nk)
    x = gc.Proj * x1

    x2 = GDB.extrapolate_ghosts(x, xb, gc)
    Anew, B = GDB.require_boundary_conditions(A, gc)

    Minv = GDB.stencil1d(3)
    Dx = GDB.derivative_matrix(1, 0, Minv, Minv, grd)

    avg = GDB.line_average(Float64[0.5 1.5; 0.5 0.5], 2, 2, GDB.stencil2d(2, 2), grd)

    return true
end

@testset "domain" begin
    @test test_domain()
end


text:
# generation of some containers filled with random values

## dicts

rand!(A::AbstractDict{K,V}, dist::Distribution{<:Pair}=Combine(Pair, K, V)) where {K,V} =
    rand!(GLOBAL_RNG, A, dist)

rand!(rng::AbstractRNG, A::AbstractDict{K,V},
      dist::Distribution{<:Pair}=Combine(Pair, K, V)) where {K,V} =
          rand!(GLOBAL_RNG, A, Sampler(rng, dist))

function _rand!(rng::AbstractRNG, A::Union{AbstractDict,AbstractSet}, n::Integer, sp::Sampler)
    empty!(A)
    while length(A) < n
        push!(A, rand(rng, sp))
    end
    A
end

rand!(rng::AbstractRNG, A::AbstractDict{K,V}, sp::Sampler) where {K,V} = _rand!(rng, A, length(A), sp)

rand(rng::AbstractRNG, dist::Distribution{<:Pair}, ::Type{T}, n::Integer) where {T<:AbstractDict} =
    _rand!(rng, deduce_type(T, eltype(dist).parameters...)(), n, Sampler(rng, dist))

rand(u::Distribution{<:Pair}, ::Type{T}, n::Integer) where {T<:AbstractDict} = rand(GLOBAL_RNG, u, T, n)

## sets

rand!(A::AbstractSet{T}, X) where {T} = rand!(GLOBAL_RNG, A, X)
rand!(A::AbstractSet{T}, ::Type{X}=T) where {T,X} = rand!(GLOBAL_RNG, A, X)

rand!(rng::AbstractRNG, A::AbstractSet, X) = rand!(rng, A, Sampler(rng, X))
rand!(rng::AbstractRNG, A::AbstractSet{T}, ::Type{X}=T) where {T,X} = rand!(rng, A, Sampler(rng, X))

_rand0!(rng::AbstractRNG, A::AbstractSet, n::Integer, X) = _rand!(rng, A, n, Sampler(rng, X))
_rand0!(rng::AbstractRNG, A::AbstractSet, n::Integer, ::Type{X}) where {X} = _rand!(rng, A, n, Sampler(rng, X))
_rand0!(rng::AbstractRNG, A::AbstractSet, n::Integer, sp::Sampler) = _rand!(rng, A, n, sp)

rand!(rng::AbstractRNG, A::AbstractSet, sp::Sampler) = _rand!(rng, A, length(A), sp)

rand(r::AbstractRNG, ::Type{T}, n::Integer) where {T<:AbstractSet} = rand(r, Float64, T, n)
rand(                ::Type{T}, n::Integer) where {T<:AbstractSet} = rand(GLOBAL_RNG, T, n)

rand(r::AbstractRNG, X, ::Type{T}, n::Integer) where {T<:AbstractSet} = _rand0!(r, deduce_type(T, eltype(X))(), n, X)
rand(                X, ::Type{T}, n::Integer) where {T<:AbstractSet} = rand(GLOBAL_RNG, X, T, n)

rand(r::AbstractRNG, ::Type{X}, ::Type{T}, n::Integer) where {X,T<:AbstractSet} = _rand0!(r, deduce_type(T, X)(), n, X)
rand(                ::Type{X}, ::Type{T}, n::Integer) where {X,T<:AbstractSet} = rand(GLOBAL_RNG, X, T, n)

## sparse vectors & matrices

rand(r::AbstractRNG, p::AbstractFloat, m::Integer) = sprand(r, m, p)
rand(                p::AbstractFloat, m::Integer) = sprand(GLOBAL_RNG, m, p)
rand(r::AbstractRNG, p::AbstractFloat, m::Integer, n::Integer) = sprand(r, m, n, p)
rand(                p::AbstractFloat, m::Integer, n::Integer) = sprand(GLOBAL_RNG, m, n, p)

rand(r::AbstractRNG, X::Sampler, p::AbstractFloat, m::Integer) =
    sprand(r, m, p, (r, n)->rand(r, X, n))

rand(r::AbstractRNG, X, p::AbstractFloat, m::Integer) =
    rand(r, Sampler(r, X), p, m)

rand(r::AbstractRNG, ::Type{X}, p::AbstractFloat, m::Integer) where {X} =
    rand(r, Sampler(r, X), p, m)

rand(X, p::AbstractFloat, m::Integer) = rand(GLOBAL_RNG, X, p, m)

rand(r::AbstractRNG, X::Sampler, p::AbstractFloat, m::Integer, n::Integer) =
    sprand(r, m, n, p, (r, n)->rand(r, X, n), eltype(X))

rand(r::AbstractRNG, X, p::AbstractFloat, m::Integer, n::Integer) =
    rand(r, Sampler(r, X), p, m, n)

rand(r::AbstractRNG, ::Type{X}, p::AbstractFloat, m::Integer, n::Integer) where {X} =
    rand(r, Sampler(r, X), p, m, n)

rand(X, p::AbstractFloat, m::Integer, n::Integer) = rand(GLOBAL_RNG, X, p, m, n)

## String

let b = UInt8['0':'9';'A':'Z';'a':'z']
    global rand
    rand(rng::AbstractRNG, chars, ::Type{String}, n::Integer=8) = String(rand(rng, chars, n))
    rand(                  chars, ::Type{String}, n::Integer=8) = rand(GLOBAL_RNG, chars, String, n)
    rand(rng::AbstractRNG, ::Type{String}, n::Integer=8) = rand(rng, b, String, n)
    rand(                  ::Type{String}, n::Integer=8) = rand(GLOBAL_RNG, b, String, n)
end

## BitArray

const BitArrays = Union{BitArray,BitVector,BitMatrix}

rand(r::AbstractRNG, ::Type{T}, dims::Dims) where {T<:BitArrays} =
    rand!(r, T(uninitialized, dims))

rand(r::AbstractRNG, ::Type{T}, dims::Integer...) where {T<:BitArrays} =
    rand!(r, T(uninitialized, convert(Dims, dims)))

rand(::Type{T}, dims::Dims) where {T<:BitArrays} =
    rand!(T(uninitialized, dims))

rand(::Type{T}, dims::Integer...) where {T<:BitArrays} =
    rand!(T(uninitialized, convert(Dims, dims)))


text:
try
    using RecordArraysBenchmarks
    true
catch
    false
end || begin
    let path = joinpath(@__DIR__, "../benchmark/RecordArraysBenchmarks/Project.toml")
        path in LOAD_PATH || push!(LOAD_PATH, path)
    end
    using RecordArraysBenchmarks
end