text:
# Test taken from http://www.stanford.edu/~clint/bench/wilk.txt

using Base.Test
using DataFrames
using StatsBase
using GLM

testeps = sqrt(eps())

nasty = DataFrame(    label = ["One", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine"],
                    x = collect(1.:9),
                    zero = fill(0.0,9),
                    miss = fill(NA, 9),
                    big = 99999990.0 + collect(1:9),
                    little = (99999990.0 + collect(1:9))/10^8,
                    huge = collect(1.:9)*1e12,
                    tiny = collect(1.:9)*1e-12,
                    round = collect(0.5:8.5))

println(nasty)
println("\nII Real Numbers:\nII A")
print("Test rounding: ")
@test [@sprintf("%1.0f", x) for x in nasty[:round]] == ["1","2","3","4","5","6","7","8","9"]
println("OK")
print("Test math: ")
@test round(Int, 2.6*7 - 0.2) == 18
@test 2 - round(Int, exp(log(sqrt(2)*sqrt(2)))) == 0
@test round(Int, 3 - exp(log(sqrt(2)*sqrt(2)))) == 1
println("OK")

print("Test means: ")
for vars in names(nasty)[2:end]
    if vars == :miss
    @test isna(mean(nasty[vars]))
    else
    @test mean(nasty[vars]) ≈ nasty[vars][5]
    end
end
println("OK")

print("Test standard deviation: ")
for vars in names(nasty)[[2;5:9]]
#    @test (@sprintf("%.9e", std(nasty[vars])))[1:10] == "2.73861278"
    @test repr(std(nasty[vars]))[1:10] == "2.73861278"
end
println("OK")
# Note: Original had one more digit but little fails then. R gives the same as Julia. Stata and SAS less precise returning 2.738612714381 and 2.738612780003404")

println("\nII D")
print("Test correlation: ")
cn = names(nasty)[[2;5:9]]
for i in 1:5
    for j = i+1:6
    @test cor(nasty[cn[i]], nasty[cn[j]]) ≈ 1
    end
end
println("OK")

print("Test spearman correlation: ")
cn = names(nasty)[[2;5:9]]
for i in 1:5
    for j = i+1:6
    @test corspearman(nasty[cn[i]], nasty[cn[j]]) ≈ 1
    end
end
println("OK")

println("\nII F")
print("Testing regression: ")
ctable = coeftable(lm(big ~ x, nasty))
@test typeof(ctable) == CoefTable
# conversion to Vector{Float64} can be removed when 0.4 is no longer supported
@test Vector{Float64}(ctable.cols[1]) ≈ [99999990, 1]

@test sprint(show, ctable) == """\
             Estimate Std.Error t value Pr(>|t|)
(Intercept)     1.0e8       0.0     Inf   <1e-99
x                 1.0       0.0     Inf   <1e-99
"""

println("OK")

println("\nIV Regression:\nIV A")
nasty[:x1] = nasty[:x]
nasty[:x2] = nasty[:x].^2
nasty[:x3] = nasty[:x].^3
nasty[:x4] = nasty[:x].^4
nasty[:x5] = nasty[:x].^5
nasty[:x6] = nasty[:x].^6
nasty[:x7] = nasty[:x].^7
nasty[:x8] = nasty[:x].^8
nasty[:x9] = nasty[:x].^9

## Is it intended that the least squares problem be overdetermined in the lm fit?
## n = 9 and p = 10 because of the implicit intercept.
lm(x1~x2+x3+x4+x5+x6+x7+x8+x9, nasty)
@test coef(lm(x~x, nasty)) ≈ [0,1]
println("OK")


text:
include("tutorial.jl")
include("MPB_wrapper.jl")

ext: jl


text:
abstract type PeriodicOrbitAlgorithm end
"""
	guessFromHopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)

This function returns several useful quantities regarding a Hopf bifurcation point. More precisely, it returns:
- the parameter value at which a Hopf bifurcation occurs
- the period of the bifurcated periodic orbit
- a guess for the bifurcated periodic orbit
- the equilibrium at the Hopf bifurcation point
- the eigenvector at the Hopf bifurcation point.

The arguments are
- `br`: the continuation branch which lists the Hopf bifurcation points
- `ind_hopf`: index of the bifurcation branch, as in `br.bifpoint`
- `eigsolver`: the eigen solver used to find the eigenvectors
- `M` number of time slices in the periodic orbit guess
- `amplitude`: amplitude of the periodic orbit guess
"""
function guessFromHopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)
	hopfpoint = HopfPoint(br, ind_hopf)
	bifpoint = br.bifpoint[ind_hopf]

	# parameter value at the Hopf point
	p_hopf = hopfpoint.p[1]

	# frequency at the Hopf point
	ωH  = hopfpoint.p[end] |> abs

	# vec_hopf is the eigenvector for the eigenvalues iω
	vec_hopf = geteigenvector(eigsolver, br.eig[bifpoint.idx][2], bifpoint.ind_bif-1)
	vec_hopf ./=  norm(vec_hopf)

	 orbitguess = [real.(hopfpoint.u .+ amplitude .* vec_hopf .* exp(-2pi * complex(0, 1) .* (ii/(M-1) - phase))) for ii=0:M-1]

	return p_hopf, 2pi/ωH, orbitguess, hopfpoint, vec_hopf
end
####################################################################################################
# Amplitude of the u component of the cycle
amplitude(x::AbstractMatrix, n) =  maximum(x[1:n, :]) - minimum(x[1:n, :])

function amplitude(x::AbstractVector, n, M; ratio = 1)
	xc = reshape(x[1:end-1], ratio * n, M)
	amplitude(xc, n)
end

function maximumPOTrap(x::AbstractVector, n, M; ratio = 1)
	xc = reshape(x[1:end-1], ratio * n, M)
	maximum(x[1:n, :])
end
####################################################################################################
# plotting function of the periodic orbits
function plotPeriodicPOTrap(x, M, Nx, Ny; ratio = 2, kwargs...)
	@assert ratio > 0 "You need at least one component"
	n = Nx*Ny
	outpo = reshape(x[1:end-1], ratio * n, M)
	po = reshape(x[1:n,1], Nx, Ny)
	rg = 2:6:M
	for ii in rg
		po = hcat(po, reshape(outpo[1:n,ii], Nx, Ny))
	end
	heatmap!(po; color = :viridis, fill=true, xlabel = "space", ylabel = "space", kwargs...)
	for ii in 1:length(rg)
		plot!([ii*Ny, ii*Ny], [1, Nx]; color = :red, width = 3, label = "", kwargs...)
	end
end

function plotPeriodicPOTrap(outpof, n, M; ratio = 2)
	@assert ratio > 0 "You need at least one component"
	outpo = reshape(outpof[1:end-1], ratio * n, M)
	if ratio == 1
		heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis)
	else
		plot(heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis),
			heatmap(outpo[n+2:end,:]', color=:viridis))
	end
end
####################################################################################################
function plotPeriodicShooting!(x, M; kwargs...)
	N = div(length(x), M);	plot!(x; label = "", kwargs...)
	for ii in 1:M
		plot!([ii*N, ii*N], [minimum(x), maximum(x)] ;color = :red, label = "", kwargs...)
	end
end

function plotPeriodicShooting(x, M; kwargs...)
	plot();plotPeriodicShooting!(x, M; kwargs...)
end


text:
mesh_file_names = ["../mesh/square_meshes_with_inclusions/mesh_test_0.0078125.g"]

# number_of_nodes = [4, 9, 25, 81, 289, 1089, 4225, 16641]
# number_of_elements = [1, 2^2, 4^2, 8^2, 16^2, 32^2, 64^2, 128^2]

# TODO: clean this up!

function test_square_mesh_with_inclusion(n::Int)
    mesh_file_name = abspath(mesh_file_names[n])
    test_name = rpad("Testing square mesh with inclusion: $(base_name(mesh_file_name))", 96)
    @testset "$test_name" begin
        # @suppress begin
            # read method test
            #
            exo_id = Exodus.open_exodus_database(abspath(mesh_file_names[n]))
            # num_dim, num_nodes, num_elem,
            # num_elem_blk, num_node_sets, num_side_sets, title =
            # Exodus.read_initialization_parameters(exo_id)
            init = Exodus.Initialization(exo_id)
            # tests on the simple numbers
            #
            @test num_dim == 2
            # @test num_nodes == number_of_nodes[n]
            # @test num_elem == number_of_elements[n]
            @test num_elem_blk == 2
            @test num_node_sets == 4
            @test num_side_sets == 4

            # test coordinates
            #
            x_coords, y_coords, z_coords =
            Exodus.read_coordinates(exo_id, num_dim, num_nodes)

            # @test size(x_coords, 1) == number_of_nodes[n]
            # @test size(y_coords, 1) == number_of_nodes[n]
            # @test size(z_coords, 1) == number_of_nodes[n]

            # read block ids
            #
            block_ids = Exodus.read_block_ids(exo_id, num_elem_blk)
            @show num_elem_blk
            @show block_ids
            @test size(block_ids, 1) == 2
            @test block_ids[1] == 1
            @test block_ids[2] == 2

            # test block initialization
            #
            # element_type, num_elem, num_nodes, num_edges, num_faces, num_attributes =
            # Exodus.read_element_block_parameters(exo_id, 1)
            #
            # @test element_type == "QUAD4"
            # @test num_elem == number_of_elements[n]
            # @test num_nodes == 4
            #
            # # test reading read_connectivity
            # #
            # conn = Exodus.read_block_connectivity(exo_id, 1)
            # @test size(conn, 1) == num_nodes * num_elem

            # block = Exodus.initialize_block(exo_id, 1)

            # test node set initialization
            #
            # node_set = Mesh.initialize_node_set(exo, 1)
            # @test node_set.node_set_number == 1

            # test mesh initialization
            #
            # mesh = Mesh.MeshStruct()

            Exodus.close_exodus_database(exo_id)
        # end
    end
end

for (n, mesh) in enumerate(mesh_file_names)
    test_square_mesh_with_inclusion(n)
    # sleep(1)
end


text:
global JTORIC_PATH = dirname(@__DIR__)

import GAP
import GAP: julia_to_gap, gap_to_julia

export GAP
export julia_to_gap, gap_to_julia

#using GAP

"""
    InstallDevelopmentVersionNConvex

Obtain a development version of NConvex.
"""
function InstallDevelopmentVersionNConvex()

    # find gap-location
    gap_location = GAP.GAPROOT

    # add pkg to this path
    p = splitpath( gap_location )
    p = replace.( p, "/"=>"")
    install_path = join( push!( p, "pkg" ), "/" )
    NConvex_path = join( push!( p, "NConvex" ), "/" )

    # inform what we are doing
    @info "Install development version of NConvex to \"" * NConvex_path * "\""

    # if this directory does exist, delete it
    if isdir( NConvex_path )

        # delete existing directory
        @info "DeleteExistingDirectory \"" * NConvex_path * "\""
        rm( NConvex_path, recursive=true )

    end

    # inform
    @info "Obtaining development version of NConvex into \"" * NConvex_path * "\""

    # prepare git operations
    res = GAP.Globals.LoadPackage(julia_to_gap("PackageManager"), false)
    @assert res
    git = julia_to_gap("git")

    # clone
    command = julia_to_gap( "clone https://github.com/HereAround/NConvex.git" )
    operation = GAP.Globals.PKGMAN_Exec(julia_to_gap( install_path ), git, command, julia_to_gap(""))
    if operation.code != 0
        @warn "Cloning of NConvex failed"
        return false
    end
    @info gap_to_julia(operation.output)

    # fetch development version
    command = julia_to_gap( "fetch origin martindevel:martindevel && git checkout martindevel" )
    operation = GAP.Globals.PKGMAN_Exec(julia_to_gap( NConvex_path ), git, command, julia_to_gap(""))
    if operation.code != 0
        @warn "Fetching of development version failed"
        return false
    end
    @info gap_to_julia(operation.output)

    # signal success
    return true

end
export InstallDevelopmentVersionNConvex

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith

# Plotting utilities of Maxima

export plot2d,
       plot3d

"""
    kwarg_convert(kwarg)

Convert keyward arguments into the form of Maxima keyword arguments

Ex:

x = (1, 2) => [x, 1, 2]
xlabel=\"x axis\" => [xlabel, \"x axis\"]
"""
function kwarg_convert(kwarg)
    str = "[$(kwarg[1]),"
    if isa(kwarg[2], Union{Tuple, Array})
        for thing in kwarg[2]
            str = str*" $thing,"
        end
    elseif isa(kwarg[2], String)
        str = str*" \"$(kwarg[2])\","
    else
        str = str*" $(kwarg[2]),"
    end
    str = str[1:end-1]
    str = str*"]"
    return str
end

"""
    plot2d(m::MExpr, kwargs...)

Plot a Maxima expression using Gnuplot.

# Examples

```julia
julia> plot2d(m"sin(x)", x=(1,2), title="Sine Wave") # Plot sine wave

```
"""
function plot2d(m::MExpr; kwargs...)
    args = String[]
    for kwarg in kwargs
        arg = kwarg_convert(kwarg)
        push!(args, arg)
    end
    str = "plot2d($m,"
    for arg in args
        str = str*" $arg,"
    end
    str = str[1:end - 1]*")"
    mcall(MExpr(str))
    return nothing
end

function plot2d(m::Array{MExpr}; kwargs...)
    str = "plot2d(["
    for expr in m
        str = str*"$expr,"
    end
    str = str[1:end-1]*"], "
    args = String[]
    for kwarg in kwargs
        arg = kwarg_convert(kwarg)
        push!(args, arg)
    end
    for arg in args
        str = str*" $arg,"
    end
    str = str[1:end-1]*")"
    mcall(MExpr(str))
    return nothing
end

"""
    contour_plot(m::MExpr, kwargs...)

Make a contour plot of the expression `m`
"""
function contour_plot(m::MExpr; kwargs...)
    args = String[]
    for kwarg in kwargs
        arg = kwarg_convert(kwarg)
        push!(args, arg)
    end
    str = "plot2d($m,"
    for arg in args
        str = str*" $arg,"
    end
    str = str[1:end - 1]*")"
    mcall(MExpr(str))
    return nothing
end

"""
    plot3d(m::MExpr, kwargs...)

Make 3d plot of the expression `m`
"""
function plot3d(m::MExpr; kwargs...)
    args = String[]
    for kwarg in kwargs
        arg = kwarg_convert(kwarg)
        push!(args, arg)
    end
    str = "plot3d($m,"
    for arg in args
        str = str*" $arg,"
    end
    str = str[1:end - 1]*")"
    mcall(MExpr(str))
    return nothing
end


text:
@test first(evaluate5(A♣, K♡, Q♢, J♠, T♢)) == 1600 # Ace-High Straight
@test first(evaluate5(K♣, Q♡, J♢, T♠, 9♢)) == 1601 # King-High Straight
@test first(evaluate5(Q♣, J♡, T♢, 9♠, 8♢)) == 1602 # Queen-High Straight
@test first(evaluate5(J♣, T♡, 9♢, 8♠, 7♢)) == 1603 # Jack-High Straight
@test first(evaluate5(T♣, 9♡, 8♢, 7♠, 6♢)) == 1604 # Ten-High Straight
@test first(evaluate5(9♣, 8♡, 7♢, 6♠, 5♢)) == 1605 # Nine-High Straight
@test first(evaluate5(8♣, 7♡, 6♢, 5♠, 4♢)) == 1606 # Eight-High Straight
@test first(evaluate5(7♣, 6♡, 5♢, 4♠, 3♢)) == 1607 # Seven-High Straight
@test first(evaluate5(6♣, 5♡, 4♢, 3♠, 2♢)) == 1608 # Six-High Straight
@test first(evaluate5(5♣, 4♡, 3♢, 2♠, A♢)) == 1609 # Five-High Straight


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using POMCPOW
using VDPTag
using MCTS
using DataFrames
using ParticleFilters
using ProfileView

pomdp = VDPTagPOMDP()

#=
rng = MersenneTwister(13)
rollout_policy = ToNextML(mdp(pomdp))
node_updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100), 0.05, rng)
solver = DPWSolver(n_iterations=10_000_000,
                   exploration_constant=40.0,
                   max_time=1.0,
                   depth=20,
                   k_action=8.0,
                   alpha_action=1/20,
                   k_state=4.0,
                   alpha_state=1/20,
                   check_repeat_state=false,
                   check_repeat_action=false,
                   estimate_value=FORollout(rollout_policy),
                   rng=rng
                  )
belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)
planner = solve(solver, belief_mdp)
=#

rollout_policy = ToNextML(mdp(pomdp), Base.GLOBAL_RNG)
solver = POMCPOWSolver(tree_queries=10_000_000,
                       criterion=MaxUCB(40.0),
                       final_criterion=MaxTries(),
                       max_depth=20,
                       k_action=8.0,
                       alpha_action=1/20,
                       k_observation=4.0,
                       alpha_observation=1/20,
                       estimate_value=FORollout(rollout_policy),
                       check_repeat_act=false,
                       check_repeat_obs=false,
                       next_action=NextMLFirst(mdp(pomdp), Base.GLOBAL_RNG),
                       default_action=TagAction(false,0.0),
                       max_time=0.1,
                       rng=Base.GLOBAL_RNG
                      )
planner = solve(solver, deepcopy(pomdp))

up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(10_000), 0.05, MersenneTwister(90_000))

b = initialize_belief(up, initial_state_distribution(pomdp))
@time action(planner, b)
@time action(planner, b)
@time action(planner, b)

Profile.clear()
@profile action(planner, b)
ProfileView.view()


text:
using BitFlags
import .Sys # explicitly imported to allow docs generation to override

@bitflag MmapProtection::Cuint begin
    PROT_NONE  = 0x00
    PROT_READ  = 0x01
    PROT_WRITE = 0x02
    PROT_EXEC  = 0x04
end
Base.cconvert(::Type{T}, pf::MmapProtection) where {T <: Integer} = T(pf)

@staticexpand @bitflag MmapFlags::Cuint begin
    MAP_FILE      = 0x00
    MAP_SHARED    = 0x01
    MAP_PRIVATE   = 0x02
    MAP_FIXED     = 0x0010
    MAP_ANONYMOUS = @static Sys.isbsd() ? 0x1000 : 0x0020
    @static if Sys.islinux()
        MAP_32BIT      =   0x0040
        MAP_GROWSDOWN  =   0x0100
        MAP_DENYWRITE  =   0x0800
        MAP_EXECUTABLE =   0x1000
        MAP_LOCKED     =   0x2000
        MAP_NORESERVE  =   0x4000
        MAP_POPULATE   =   0x8000
        MAP_NONBLOCK   = 0x1_0000
        MAP_STACK      = 0x2_0000
        MAP_HUGETLB    = 0x4_0000
        MAP_SYNC       = 0x8_0000
        MAP_FIXED_NOREPLACE = 0x10_0000
    end
    @static if Sys.isapple() || (VERSION >= v"1.1" && (Sys.isnetbsd() || Sys.isdragonfly()))
        MAP_RENAME       = 0x0020
        MAP_NORESERVE    = 0x0040
        MAP_INHERIT      = 0x0080
        MAP_NOEXTEND     = 0x0100
        MAP_HASSEMAPHORE = 0x0200
    end
    @static if VERSION >= v"1.1" && Sys.isfreebsd()
        MAP_STACK   =   0x0400
        MAP_NOSYNC  =   0x0800
        MAP_GUARD   =   0x2000
        MAP_EXCL    =   0x4000
        MAP_NOCORE  = 0x4_0000
        MAP_32BIT   = 0x8_0000
    elseif VERSION >= v"1.1" && Sys.isdragonfly()
        MAP_STACK      =   0x0400
        MAP_NOSYNC     =   0x0800
        MAP_VPAGETABLE =   0x2000
        MAP_TRYFIXED   = 0x1_0000
        MAP_NOCORE     = 0x2_0000
        MAP_SIZEALIGN  = 0x4_0000
    elseif VERSION >= v"1.1" && Sys.isopenbsd()
        MAP_STACK   = 0x4000
        MAP_CONCEAL = 0x8000
    elseif VERSION >= v"1.1" && Sys.isnetbsd()
        MAP_REMAPDUP = 0x0004
        MAP_TRYFIXED = 0x0400
        MAP_WIRED    = 0x0800
        MAP_STACK    = 0x2000
    end
end
Base.cconvert(::Type{T}, mf::MmapFlags) where {T <: Integer} = T(mf)

# Provide this weird mode?
#@static if Sys.islinux()
#const MAP_SHARED_VALIDATE = MmapFlags(0x03)
#end

@staticexpand @enum AdviseFlags::Cint begin
    MADV_NORMAL = 0
    MADV_RANDOM = 1
    MADV_SEQUENTIAL = 2
    MADV_WILLNEED = 3
    MADV_DONTNEED = 4
    @static if Sys.islinux()
        MADV_FREE = 8
        MADV_REMOVE = 9
        MADV_DONTFORK = 10
        MADV_DOFORK = 11
        MADV_MERGEABLE = 12
        MADV_UNMERGEABLE = 13
        MADV_HUGEPAGE = 14
        MADV_NOHUGEPAGE = 15
        MADV_DONTDUMP = 16
        MADV_DODUMP = 17
        MADV_WIPEONFORK = 18
        MADV_KEEPONFORK = 19
        MADV_COLD = 20
        MADV_PAGEOUT = 21
        MADV_HWPOISON = 100
        MADV_SOFT_OFFLINE = 101
    elseif Sys.isapple()
        MADV_FREE = 5
    elseif VERSION >= v"1.1" && (Sys.isfreebsd() || Sys.isdragonfly())
        MADV_FREE = 5
        MADV_NOSYNC = 6
        MADV_AUTOSYNC = 7
        MADV_NOCORE = 8
        MADV_CORE = 9
        @static if Sys.isfreebsd()
            MADV_PROTECT = 10
        else
            MADV_INVAL = 10
            MADV_SETMAP = 11
        end
    elseif VERSION >= v"1.1" && (Sys.isopenbsd() || Sys.isnetbsd())
        MADV_SPACEAVAIL = 5
        MADV_FREE = 6
    end
end
Base.cconvert(::Type{T}, af::AdviseFlags) where {T <: Integer} = T(af)

# documentation requires the types to be defined first, so do after the fact
let _flag_docs
    function _flag_docs(T)
        buf = IOBuffer()
        for c in instances(T)
            println(buf, c)
        end
        return rstrip(String(take!(buf)))
    end

    @doc """
        @bitflag UnixMmap.MmapProtection

    Set of bit flags which control the memory protections applied to the memory mapped
    region. The flag should be either `PROT_NONE` or some bitwise-or combination of
    the remaining flags.

    The flags available on $(Sys.KERNEL) are:
    ```julia
    $(_flag_docs(MmapProtection))
    ```
    """ MmapProtection

    @doc """
        @bitflag UnixMmap.MmapFlags

    Set of bit flags which control the handling of the memory mappged region. The flag
    should be a bitwise-or combination of flags.

    The flags available on $(Sys.KERNEL) are:
    ```julia
    $(_flag_docs(MmapFlags))
    ```
    """ MmapFlags

    @doc """
        @enum UnixMmap.AdviseFlags

    Set of flags which advise the kernel on handling of the memory mapped region.

    The flags available on $(Sys.KERNEL) are:
    ```julia
    $(_flag_docs(AdviseFlags))
    ```
    """ AdviseFlags
end


text:
using Documenter
using JsonGrinder
using Literate

const is_ci = haskey(ENV, "GITHUB_ACTIONS")

DocMeta.setdocmeta!(JsonGrinder, :DocTestSetup, :(using JsonGrinder); recursive=true)

# generate files using literate.jl
src_dir = joinpath(@__DIR__, "src")
examples_dir = joinpath(@__DIR__, "..", "examples")
examples_generated_dir = joinpath(src_dir, "examples")
!ispath(examples_generated_dir) && mkpath(examples_generated_dir)
example_files = [joinpath(examples_dir, f) for f in [
    "mutagenesis.jl",
    "recipes.jl",
    "schema_examination.jl",
    "schema_visualization.jl",
]]

# str = read(joinpath(@__DIR__, "src", "examples", "schema_visualization_raw.md"), String);
function print_html_raw(str)
    str
    lines = split(str, "\n")
    html_line = findfirst(s -> occursin("<!DOCTYPE html>", s) && occursin("</html>\\n\"", s), lines)
    if isnothing(html_line)
        return str
    end
    lines[html_line-1] *= "@raw html"
    lines[html_line] = replace(lines[html_line], "\\\"" => "\"")
    lines[html_line] = replace(lines[html_line], "\\n" => "\n")
    lines[html_line] = replace(lines[html_line], "\\t" => "\t")
    lines[html_line] = replace(lines[html_line], r"\\(\w+)" => s"\1")
    join(lines, "\n")
end

example_mds = []

#for purpose of local testing
config = is_ci ? Dict() : Dict("nbviewer_root_url"=>"https://nbviewer.jupyter.org/github/CTUAvastLab/JsonGrinder.jl/blob/gh-pages/previews/PR101")

# because there is a bug that @example blocks are not evaluated when they are included using @eval, I run the markdown code in Literate.jl
for f in example_files
    md_file = Literate.markdown(f, examples_generated_dir; config, credit = false, execute = true, postprocess = print_html_raw)
    Literate.notebook(f, examples_generated_dir; config, execute = is_ci) # Don't execute locally, because it takes long
    push!(example_mds, relpath(md_file, dirname(examples_generated_dir)))
end

# for running only doctests
is_ci && doctest(JsonGrinder)
makedocs(
         sitename = "JsonGrinder.jl",
         # doctest = false,
         format = Documenter.HTML(sidebar_sitename=false,
                                  prettyurls=get(ENV, "CI", nothing) == "true",
                                  assets=["assets/favicon.ico", "assets/custom.css"]),
         modules = [JsonGrinder],
         pages = [
                  "Home" => "index.md",
                  "Schema" => "schema.md",
                  "Creating extractors" => "extractors.md",
                  "Extractors overview" => "exfunctions.md",
                  "Examples" => example_mds,
                  "AutoML" => "automl.md",
                  "External tools" => "hierarchical.md",
                  "API Documentation" => "api.md",
                  "Developers" => "developers.md",
                  "Citation" => "citation.md",
                  ],

)

deploydocs(
    repo = "github.com/CTUAvastLab/JsonGrinder.jl.git",
    push_preview = true
)


text:
import StaticArrays.SVector
import GeometryBasics: Point
import Rotations.RotMatrix

export Point, Path, MultiPath, unitvec, point_matrix

"""
The unit vector (0, 1).

Useful for creating n-gons, e.g.:

```jldoctest
using PenPlots

expect = [
    [-1, 0],
    [0, -1],
    [1, 0],
    [0, 1],
]

isapprox(expect, [frac_rotation(i/4) * unitvec for i in 1:4])

# output

true
```
"""
const unitvec = Point(0, 1)

"""
A 1-dimensional array of `Point`s.
"""
struct Path <: AbstractVector{Point}
    points::Vector{Point}
end

function Base.size(path::Path)
    size(path.points)
end

function Base.getindex(path::Path, i::Int64)
    getindex(path.points, i)
end

"""
An array of paths. Corresponds one-to-one with a layer in a plot.
"""
const MultiPath = Vector{Path}

# Rotation by a matrix.

function Base.:*(rot::RotMatrix, path::Path)
    Path([rot * point for point in path.points])
end

function Base.:*(rot::RotMatrix, paths::MultiPath)
    MultiPath([rot * path for path in paths])
end

# Multiplication by a point.

function Base.:*(point::Point, path::Path)
    Path([point * x for x in path])
end

function Base.:*(point::Point, paths::MultiPath)
    MultiPath([point * path for path in paths])
end

# Addition with a point.

function Base.:+(point::Point, path::Path)
    Path([point + x for x in path])
end

function Base.:+(point::Point, paths::MultiPath)
    MultiPath([point + path for path in paths])
end

"""
    point_matrix(xs, ys)

Generate a matrix by producing `Point`s for the cartesian product of the vectors
of `x` and `y` values given. The resulting matrix will have as many columns as
entries of `xs`, and as many rows as values of `ys`.
"""
function point_matrix(xs, ys)
	xxs = xs' .* ones(length(ys))
	yys = ys .* ones(length(xs))'
	Point.(xxs, yys)
end


text:
export @vars

function process_vars!(syms, typs, expr)
    if expr isa LineNumberNode
        return
    elseif expr isa Expr && expr.head == :block
        for arg in expr.args
            process_vars!(syms, typs, arg)
        end
        return
    elseif expr.head == :(::)
        push!(syms, expr.args[1])
        push!(typs, expr.args[2])
        return
    else
        error("Invalid expression")
    end
end

macro vars(args...)
    syms = Any[]
    typs = Any[]
    for arg in args
        process_vars!(syms, typs, arg)
    end
    :(NamedTuple{$(tuple(syms...)), Tuple{$(esc.(typs)...)}})
end


text:
wavenumber_low_volfrac(ω::Number, medium::Medium, specie::Specie; kws...) =
  wavenumber_low_volfrac(ω, medium, [specie]; kws...)

wavenumber_low_volfrac(ωs::AbstractVector{T}, medium::Medium{T}, species::Vector{Specie{T}}; kws...) where T<:Number =
  [wavenumber_low_volfrac(ω, medium, species; kws...) for ω in ωs]

wavenumber_low_volfrac(ωs::AbstractVector{T}, medium::Medium{T}, specie::Specie{T}; kws...) where T<:Number = [wavenumber_low_volfrac(ω, medium, [specie]; kws...) for ω in ωs]

"Explicit formula for one effective wavenumber based on a low particle volume fraction expasion."
function wavenumber_low_volfrac(ω::T, medium::Medium{T}, species::Vector{Specie{T}}; tol::T =1e-6,
    hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
    radius_multiplier::T = 1.005, verbose::Bool = true) where T <: Number

  volume_fraction = sum(pi*sp.r^2.0*sp.num_density for sp in species)
  if volume_fraction >= 0.4 && verbose
    @warn("the volume fraction $(round(100*volume_fraction))% is too high, expect a relative error of approximately $(round(100*volume_fraction^3.0))%")
  end
  num_density = sum(sp.num_density for sp in species)
  # Add incident wavenumber
  kT2 = (ω/medium.c)^2.0
  # Add far-field contribution
  kT2 += - 4.0im*num_density*far_field_pattern(ω, medium, species; hankel_order=hankel_order)(0.0)
  # Add pair-field contribution
  kT2 += - 4.0im*num_density^(2.0)*pair_field_pattern(ω, medium, species;
      radius_multiplier=radius_multiplier, hankel_order=hankel_order)(0.0)

  return (imag(sqrt(kT2)) > zero(T)) ? sqrt(kT2) : -sqrt(kT2)
end

reflection_coefficient_low_volfrac(ωs::AbstractVector{T}, medium::Medium{T}, species::Vector{Specie{T}}; kws... ) where T<:Number =
    [reflection_coefficient_low_volfrac(ω, medium, species; kws... ) for ω in ωs]

"An explicit formula for the refleciton coefficient based on a low particle volume fraction."
function reflection_coefficient_low_volfrac(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        θin::T = zero(T), kws... ) where T<:Number

    θ_ref = T(π) - T(2)*θin
    fo = far_field_pattern(ω, medium, species; kws...)
    dfo = diff_far_field_pattern(ω, medium, species; kws...)
    foo = pair_field_pattern(ω, medium, species; kws...)

    k = ω/medium.c
    α = k*cos(θin)
    R1 = im*fo(θ_ref)
    R2 = 2.0*fo(zero(T))/(α^2.0)
    R2 = im*foo(θ_ref) + R2*(sin(θin)*cos(θin)*dfo(θ_ref) - fo(θ_ref))

    num_density = sum(s.num_density for s in species)
    R = (R1 + num_density*R2)*num_density/(α^2.0)
    return R
end

function wavenumber_very_low_volfrac(ω::Number, medium::Medium, species::Array{Specie{T}};
    tol=1e-6, radius_multiplier = 1.005, verbose = false)  where T<:Number

  volume_fraction = sum(pi*sp.r^2.0*sp.num_density for sp in species)
  if volume_fraction >= 0.4
    @warn("the volume fraction $(round(100*volume_fraction))% is too high, expect a relative error of approximately $(round(100*volume_fraction^3.0))%")
  end
  kT2 = (ω/medium.c)^2.0
  next_order = 4.0im*sum(sp.num_density*Zn(ω,sp,medium,0) for sp in species)
  hankel_order=1

  # sum more hankel orders until the relative error < tol
  while abs(next_order/kT2) > tol
    kT2 += next_order
    next_order = 4.0im*sum(sp.num_density*Zn(ω,sp,medium,m) for sp in species, m in (-hankel_order,hankel_order))
    hankel_order +=1
  end
  kT2 += next_order
  hankel_order +=1
  if verbose println("max Hankel order = $hankel_order") end

  return sqrt(kT2)
end


text:
using LinearAlgebra
using Roots
using Distances

export exit_entry_times, transit_return, mean_return_times

@deprecate transit_time_statistics entry_exit_times

"""
    exit_entry_times(ds, u₀, εs, T; diffeq...) → exits, entries
Collect exit and entry times for a ball/box centered at `u₀` with radii `εs` (see below),
in the state space of the given discrete dynamical system (function not yet available
for continuous systems).
Return the exit and (re-)entry return times to the set(s), where each of these is a vector
containing all collected times for the respective `ε`-radius set, for `ε ∈ εs`.

Use `transit_return(exits, entries)` to transform the output into transit and return
times, and see also [`mean_return_times`](@ref) for both continuous and discrete systems.

## Description
Transit time statistics are important for the transport properties of dynamical systems[^Meiss1997]
and can even be connected with the fractal dimension of chaotic sets[^Boev2014].

The current algorithm collects exit and re-entry times to given sets in the state space,
which are centered at `u₀` (**algorithm always starts at `u₀`** and the initial state of
`ds` is irrelevant). `εs` is always a `Vector`.

The sets around `u₀` are nested hyper-spheres of radius `ε ∈ εs`, if each entry of
`εs` is a real number. The sets can also be
hyper-rectangles (boxes), if each entry of `εs` is a vector itself.
Then, the `i`-th box is defined by the space covered by `u0 .± εs[i]` (thus the actual
box size is `2εs[i]`!).

The reason to input multiple `εs` at once is purely for performance.

For discrete systems, exit time is recorded immediatelly after exitting of the set, and
re-entry is recorded immediatelly on re-entry. This means that if an orbit needs
1 step to leave the set and then it re-enters immediatelly on the next step, the return time
is 1. For continuous systems high-order
interpolation is done to accurately record the time of exactly crossing the `ε`-ball/box.

[^Meiss1997]: Meiss, J. D. *Average exit time for volume-preserving maps*, Chaos (1997)](https://doi.org/10.1063/1.166245)

[^Boev2014]: Boev, Vadivasova, & Anishchenko, *Poincaré recurrence statistics as an indicator of chaos synchronization*, Chaos (2014)](https://doi.org/10.1063/1.4873721)
"""
function exit_entry_times end

"""
    mean_return_times(ds::DynamicalSystem, u₀, εs, T; kwargs...) → τ, c
Return the mean return times `τ`, as well as the amount of returns `c`, for
subsets of the state space of `ds` defined by `u₀, εs`.
The `ds` is evolved for a maximum of `T` time.
This function behaves similarly to [`exit_entry_times`](@ref) and thus see that one for
the meaning of `u₀` and `εs`.

This function supports both discrete and continuous systems, however the optimizations
done in discrete systems (where all nested `ε`-sets are checked at the same time),
are not done here yet, which leads to disproportionally lower performance since
each `ε`-related set is checked individually from start.

Continuous systems allow for the following keywords:

- `i=10` How many points to interpolate the trajectory in-between steps to find
  candidate crossing regions.
- `dmin` If the trajectory is at least `dmin` distance away from `u0`,
  the algorithm that checks for crossings of the `ε`-set is not initiated.
  By default obtains the a value 4 times as large as the radius of the maximum ε-set.
- `diffeq...` All extra keywords are propagated to solvers of DifferentialEquations.jl,
  see [`trajectory`](@ref). It is strongly recommended to use high accuracy keywords here,
  e.g. `alg = Vern9(), reltol = 1e-12, abstol = 1e-12, maxiters = typemax(Int)`.
For continuous systems `T, i, dmin` can be vectors with same size as `εs`, to help increase
accuracy of small `ε`.
"""
function mean_return_times end

"""
    transit_return(exits, entries) → transit, return
Convert the output of [`exit_entry_times`](@ref) to the transit and return times.
"""
function transit_return(exits, entries)
    # the main reason this function exists is because entry times can be one less
    # than entry times. (otherwise you could just directly subtract)
    returns = [en .- view(ex, 1:length(en)) for (en, ex) in zip(entries, exits)]
    transits = similar(entries)
    for (j, (en, ex)) in enumerate(zip(entries, exits))
        M, N = length(en), length(ex)
        enr, exr = M == N ? (1:N-1, 2:N) : (1:M, 2:N)
        transits[j] = view(ex, exr) .- view(en, enr)
    end
    return transits, returns
end

##########################################################################################
# ε-distances
##########################################################################################
function check_εs_sorting(εs, L)
    correct = if εs[1] isa Real
        issorted(εs; rev = true)
    elseif εs[1] isa AbstractVector
        @assert all(e -> length(e) == L, εs) "Boxes must have same dimension as state space!"
        for j in 1:L
            if !issorted([εs[i][j] for i in 1:length(εs)]; rev = true)
                return false
            end
        end
        true
    end
    if !correct
        throw(ArgumentError("`εs` must be sorted from largest to smallest ball/box size."))
    end
    return correct
end

# Support both types of sets: balls and boxes
"Return `true` if state is outside ε-ball"
function isoutside(u, u0, ε::AbstractVector)
    @inbounds for i in 1:length(u)
        abs(u[i] - u0[i]) > ε[i] && return true
    end
    return false
end
isoutside(u, u0, ε::Real) = euclidean(u, u0) > ε

"Return the **signed** distance of state to ε-ball (negative means inside ball)"
function εdistance(u, u0, ε::AbstractVector)
    m = eltype(u)(-Inf)
    @inbounds for i in 1:length(u)
        m2 = abs(u[i] - u0[i]) - ε[i]
        m2 > m && (m = m2)
    end
    return m
end

εdistance(u, u0, ε::Real) = euclidean(u, u0) - ε

"Return the true distance of `u` from `u0` according to metric defined by `ε`."
function distance(u, u0, ε::AbstractVector)
    m = eltype(u)(-Inf)
    @inbounds for i in 1:length(u)
        m2 = abs(u[i] - u0[i])
        m2 > m && (m = m2)
    end
    return m
end
distance(u, u0, ε::Real) = euclidean(u, u0)

##########################################################################################
# Discrete systems
##########################################################################################
function exit_entry_times(ds::DiscreteDynamicalSystem, u0, εs, T; diffeq...)
    check_εs_sorting(εs, length(u0))
    integ = integrator(ds, u0)
    exit_entry_times(integ, u0, εs, T)
end

function exit_entry_times(integ::MDI, u0, εs, T)
    E = length(εs)
    pre_outside = fill(false, E)      # `true` if outside the set. Previous step.
    cur_outside = copy(pre_outside)   # `true` if outside the set. Current step.
    exits   = [Int[] for _ in 1:E]
    entries = [Int[] for _ in 1:E]

    while integ.t < T
        step!(integ)

        # here i gives the index of the largest ε-set that the trajectory is out of.
        # It is guaranteed that the trajectory is thus outside all other boxes
        i = first_outside_index(integ, u0, εs, E)
        cur_outside[i:end] .= true
        cur_outside[1:i-1] .= false

        update_exit_times!(exits, i, pre_outside, cur_outside, integ)
        update_entry_times!(entries, i, pre_outside, cur_outside, integ)
        pre_outside .= cur_outside
    end
    return exits, entries
end

function first_outside_index(integ::MDI, u0, εs, E)::Int
    i = findfirst(e -> isoutside(integ.u, u0, e), εs)
    return isnothing(i) ? E+1 : i
end

function update_exit_times!(exits, i, pre_outside, cur_outside, integ::MDI)
    @inbounds for j in i:length(pre_outside)
        cur_outside[j] && !pre_outside[j] && push!(exits[j], integ.t)
    end
end

function update_entry_times!(entries, i, pre_outside, cur_outside, integ::MDI)
    @inbounds for j in 1:i-1
        pre_outside[j] && !cur_outside[j] && push!(entries[j], integ.t)
    end
end

function mean_return_times(ds::DiscreteDynamicalSystem, u0, εs, T; diffeq...)
    check_εs_sorting(εs, length(u0))
    integ = integrator(ds, u0)
    mean_return_times(integ, u0, εs, T)
end

function mean_return_times(integ::MDI, u0, εs, T)
    E = length(εs)
    pre_outside = fill(false, E)      # `true` if outside the set. Previous step.
    cur_outside = copy(pre_outside)   # `true` if outside the set. Current step.
    exits, entries, counts, mrt = (zeros(Int, E) for _ in 1:4)

    while integ.t < T
        step!(integ)

        # here i gives the index of the largest ε-ball that the trajectory is out of.
        # It is guaranteed that the trajectory is thus outside all other boxes
        i = first_outside_index(integ, u0, εs, E)
        cur_outside[i:end] .= true
        cur_outside[1:i-1] .= false

        update_exit_times_single!(exits, i, pre_outside, cur_outside, integ.t)
        update_entries_and_returns!(
            mrt, counts, entries, exits, i, pre_outside, cur_outside, integ.t
        )
        pre_outside .= cur_outside
    end
    return mrt ./ counts, counts
end

function update_exit_times_single!(exits, i, pre_outside, cur_outside, t)
    @inbounds for j in i:length(pre_outside)
        cur_outside[j] && !pre_outside[j] && (exits[j] = t)
    end
end

function update_entries_and_returns!(
        mrt, counts, entries, exits, i, pre_outside, cur_outside, t
    )
    @inbounds for j in 1:i-1 # only when inside the set it makes sense to calculate returns
        if pre_outside[j] && !cur_outside[j] # we just entered the set
            entries[j] = t
            mrt[j] += entries[j] - exits[j]
            counts[j] += 1
        end
    end
end

##########################################################################################
# Continuous
##########################################################################################
using DynamicalSystemsBase.DiffEqBase: ODEProblem, solve
using DynamicalSystemsBase.DiffEqBase: ContinuousCallback, CallbackSet

function exit_entry_times(ds::ContinuousDynamicalSystem, u0, εs, T;
        diffeq...
    )

    error("Continuous system version is not yet ready.")

    eT = eltype(ds.t0)
    check_εs_sorting(εs, length(u0))
    exits = [eT[] for _ in 1:length(εs)]
    entries = [eT[] for _ in 1:length(εs)]

    # This callback fails, see https://github.com/SciML/DiffEqBase.jl/issues/580
    # callbacks = ContinuousCallback[]
    # for i in eachindex(εs)
    #     crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[i])
    #     negative_affect!(integ) = push!(entries[i], integ.t)
    #     positive_affect!(integ) = push!(exits[i], integ.t)
    #     cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
    #         save_positions = (false, false)
    #     )
    #     push!(callbacks, cb)
    # end
    # cb = CallbackSet(callbacks...)

    # This callback works, but it is only for 1 nesting level
    crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[1])
    negative_affect!(integ) = push!(entries[1], integ.t)
    positive_affect!(integ) = push!(exits[1], integ.t)
    cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
        save_positions = (false, false)
    )

    prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
    sol = solve(prob;
        callback=cb, save_everystep = false, dense = false,
        save_start=false, save_end = false,
        CDS_KWARGS..., diffeq...
    )
    return exits, entries
end

function _default_dmin(εs)
    if εs[1] isa Real
        m = 4*maximum(εs)
    else
        m = 8*maximum(maximum(e) for e in εs)
    end
    return m
end

function mean_return_times(ds::ContinuousDynamicalSystem, u0, εs, T;
        i=10, dmin=_default_dmin(εs), diffeq...
    )
    # This warning would be useful if the callback methods worked
    # if !haskey(diffeq, :alg) || diffeq[:alg] == DynamicalSystemsBase.DEFAULT_SOLVER
    #     error(
    #     "Please use a solver that supports callbacks using `OrdinaryDiffEq`. "*
    #     "For example `using OrdinaryDiffEq: Tsit5; mean_return_times(...; alg = Tsit5())`."
    #     )
    # end
    if haskey(diffeq, :m)
        @warn "Keyword `m` is removed in favor of `dmin`."
    end

    eT = eltype(ds.t0)
    check_εs_sorting(εs, length(u0))
    c = zeros(Int, length(εs)); τ = zeros(eT, length(εs))
    integ = integrator(ds, u0; diffeq...)
    for j ∈ 1:length(εs)
        reinit!(integ)
        t = T isa AbstractVector ? T[j] : T
        μ = dmin isa AbstractVector ? dmin[j] : dmin
        ι = i isa AbstractVector ? i[j] : i
        τ[j], c[j] = mean_return_times_single(integ, u0, εs[j], t; i=ι, dmin=μ)
    end
    return τ, c
end

function mean_return_times_single(
        integ, u0, ε, T;
        i, dmin, rootkw = (xrtol = 1e-12, atol = 1e-12)
    )

    exit = integ.t
    τ, c = zero(exit), 0
    isoutside = false
    crossing(t) = εdistance(integ(t), u0, ε)

    while integ.t < T
        step!(integ)
        # Check distance of uprev (because interpolation can happen only between
        # tprev and t) and if it is "too far away", then don't bother checking crossings.
        d = distance(integ.uprev, u0, ε)
        d > dmin && continue

        r = range(integ.tprev, integ.t; length = i)
        dp = εdistance(integ.uprev, u0, ε) # `≡ crossing(r[1])`, crossing of previous step
        for j in 2:i
            dc = crossing(r[j])
            if dc*dp < 0 # the distances have different sign (== 0 case is dismissed)
                tcross = Roots.find_zero(crossing, (r[j-1], r[j]), ROOTS_ALG; rootkw...)
                if dp < 0 # here the trajectory goes from inside to outside
                    exit = tcross
                    break # if we get out of the box we don't check whether we go inside
                          # again: assume in 1 step we can't cross out two times
                else # otherwise the trajectory goes from outside to inside
                    exit > tcross && continue # scenario where exit was not calculated
                    τ += tcross - exit
                    c += 1
                    # here we don't `break` because we might cross out in the same step
                end
            end
            dp = dc
        end
        # Notice that the possibility of a trajectory going fully through the `ε`-set
        # within one of the interpolation points is NOT considered: simply increase
        # `interp_points` to increase accuracy and make this scenario possible.
    end
    return τ/c, c
end

# There seem to be fundamental problems with this method. It runs very slow,
# and does not seem to terminate........
# Furthermore, I have a suspision that calling a callback
# affects the trajectory solution EVEN IF NO MODIFICATION IS DONE to the `integ` object.
# I confirmed this by simply evolving the trajectory and looking at the plotting code
# at the `transit_time_tests.jl` file.
function mean_return_times_single_callbacks(
        ds::ContinuousDynamicalSystem, u0, ε, T;
        interp_points=10, diffeq...
    )

    eT = eltype(ds.t0)
    eeτc = zeros(eT, 4) # exit, entry, τ and count are the respective entries
    isoutside = Ref(false)

    crossing(u, t, integ) = ChaosTools.εdistance(u, u0, ε)
    function positive_affect!(integ)
        # println("Positive affect (crossing out) at time $(integ.t) and state:")
        # println(integ.u)
        if isoutside[]
            return
        else
            # When crossing ε-set outwards, outside becomes true
            isoutside[] = true
            eeτc[1] = integ.t
            return
        end
    end
    function negative_affect!(integ)
        # println("Negative affect (crossing in) at time $(integ.t) and state:")
        # println(integ.u)
        # we have crossed the ε-set, but let's check whether we were already inside
        if isoutside[]
            isoutside[] = false
            eeτc[2] = integ.t
            eeτc[3] += eeτc[2] - eeτc[1]
            eeτc[4] += 1
            return
        else
            return
        end
    end
    cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
        save_positions = (false, false), interp_points,
    )

    prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
    sol = solve(prob;
        callback=cb, save_everystep = false, dense = false,
        save_start=false, save_end = false, diffeq...
    )
    return eeτc[3]/eeτc[4], round(Int, eeτc[4])
end


text:
function SphericalCombinedEigenfunction(
    domain::SphericalTriangle,
    us::Vector{<:AbstractSphericalEigenfunction},
)
    orders = ones(Int, length(us))
    !isempty(us) || throw(ArgumentError("us must be non-empty"))

    SphericalCombinedEigenfunction(domain, us, orders)
end

function Base.show(io::IO, u::SphericalCombinedEigenfunction)
    println(io, "Combined eigenfunction")
    println(io, "domain: $(u.domain)")
    println(io, "eigenfunctions:")
    recur_io = IOContext(io, :compact => true)
    for v in u.us
        print(io, "- ")
        show(recur_io, v)
    end
    print(io, "number of set coefficients: $(length(coefficients(u)))")
end

function set_domain!(u::SphericalCombinedEigenfunction, domain::AbstractDomain)
    u.domain = domain
    for v in u.us
        set_domain!(v, domain)
    end
    return u
end

function recompute!(u::SphericalCombinedEigenfunction)
    for v in u.us
        recompute!(v)
    end
    u
end

function active_boundaries(domain::SphericalTriangle, u::SphericalCombinedEigenfunction)
    if domain === u.domain
        return union([active_boundaries(domain, v) for v in u.us]...)
    else
        return 1:3
    end
end

"""
    active_eigenfunctions(u::SphericalCombinedEigenfunction,
                          i::Integer)
> Return a SphericalCombinedEigenfunction consisting of all
  eigenfunction in u which are active on boundary i.

  If there are no such eigenfunctions then nothing is returned.

  The returned eigenfunction shares coefficients with the old one, any
  changes to the first also applies to the latter and vice versa.

  See also: [`active_boundaries`](@ref)
"""
function active_eigenfunctions(u::SphericalCombinedEigenfunction, i::Integer)
    indices = [j for j in 1:length(u.us) if i ∈ active_boundaries(u.us[j])]

    if isempty(indices)
        return nothing
    end

    vs = u.us[indices]
    orders = u.orders[indices]

    SphericalCombinedEigenfunction(u.domain, vs, orders)
end

"""
    basis_function(u::SphericalCombinedEigenfunction,
                   k::Integer)
> Return the index for the eigenfunctions that u is a combination of
  that should be used and the index to use for that eigenfunction.

  This depends on the relative order between the eigenfunctions. It
  returns i, j which corresponds to calling the j-th basis-function of
  u.us[i].
"""
function basis_function(u::SphericalCombinedEigenfunction, k::Integer)
    k > 0 || throw(ArgumentError("k must be positive not $k"))
    i = 1
    j = 0
    while k > u.orders[i]
        k -= u.orders[i]

        i += 1
        if i > length(u.orders)
            i = 1
            j += 1
        end
    end

    i, u.orders[i] * j + k
end

"""
    basis_function(u::SphericalCombinedEigenfunction, ks::UnitRange{Int})

Return a vector with elements of type `UnitRange{Int}` where element
`i` is the range of indices that should be used for the respective
eigenfunctions, i.e. the eigenfunction `u.us[i]` should take the range
given by the `i`th element in this vector.
"""
function basis_function(u::SphericalCombinedEigenfunction, ks::UnitRange{Int})
    ks.start == 1 || throw(ArgumentError("ks must start with 1, got ks = $ks"))
    ks.stop == 0 && return fill(1:0, length(u.us))

    A = [0; cumsum(u.orders)]
    fullcycles = div(ks.stop, A[end])
    R = ks.stop % A[end]
    remaining = [max(min(u.orders[i], R - A[i]), 0) for i in eachindex(u.us)]
    return [1:(fullcycles*u.orders[i]+remaining[i]) for i in eachindex(u.us)]
end

function coefficients(u::SphericalCombinedEigenfunction)
    coeffs = [coefficients(v) for v in u.us]
    N = sum(length, coeffs)
    [coeffs[i][j] for (i, j) in [basis_function(u, k) for k = 1:N]]
end

function set_eigenfunction!(u::SphericalCombinedEigenfunction, coefficients::Vector)
    is = [basis_function(u, k)[1] for k = 1:length(coefficients)]
    for i = 1:length(u.us)
        set_eigenfunction!(u.us[i], coefficients[is.==i])
    end
end

function (u::SphericalCombinedEigenfunction)(
    xyz::AbstractVector{T},
    λ::arb,
    k::Integer;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    i, j = basis_function(u, k)
    u.us[i](xyz, λ, j, boundary = boundary, notransform = notransform)
end

function (u::SphericalCombinedEigenfunction)(
    θ::T,
    ϕ::T,
    λ::arb,
    k::Integer;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    i, j = basis_function(u, k)
    u.us[i](θ, ϕ, λ, j, boundary = boundary, notransform = notransform)
end

function (u::SphericalCombinedEigenfunction)(
    xyz::AbstractVector{T},
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    ks_per_index = basis_function(u, ks)

    res_per_index = similar(u.us, Vector{T})
    for i in eachindex(u.us)
        res_per_index[i] =
            u.us[i](xyz, λ, ks_per_index[i], boundary = boundary, notransform = notransform)
    end

    res = similar(ks, T)
    # TODO: This can be done more efficiently
    for i in eachindex(ks)
        j, l = basis_function(u, ks[i])
        res[i] = res_per_index[j][l]
    end

    return res
end

function (u::SphericalCombinedEigenfunction)(
    θ::T,
    ϕ::T,
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    ks_per_index = basis_function(u, ks)

    res_per_index = similar(u.us, Vector{T})
    for i in eachindex(u.us)
        res_per_index[i] = u.us[i](
            θ,
            ϕ,
            λ,
            ks_per_index[i],
            boundary = boundary,
            notransform = notransform,
        )
    end

    res = similar(ks, T)
    # TODO: This can be done more efficiently
    for i in eachindex(ks)
        j, l = basis_function(u, ks[i])
        res[i] = res_per_index[j][l]
    end

    return res
end

function (u::SphericalCombinedEigenfunction)(
    xyz::AbstractVector{T},
    λ::arb;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    res = u.domain.parent(0)

    for v in u.us
        res += v(xyz, λ, boundary = boundary, notransform = notransform)
        if (T == arb && !isfinite(res)) || (T == arb_series && !isfinite(res[end]))
            return res
        end
    end

    res
end

function (u::SphericalCombinedEigenfunction)(
    θ::T,
    ϕ::T,
    λ::arb;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    res = u.domain.parent(0)

    for v in u.us
        res += v(θ, ϕ, λ, boundary = boundary, notransform = notransform)
        if (T == arb && !isfinite(res)) || (T == arb_series && !isfinite(res[end]))
            return res
        end
    end

    res
end


text:
using Base.Iterators: repeated, take, enumerate
using ProtoGrad
using ProtoGrad: Linear, Compose, ReLU, SupervisedObjective, mse
using LinearAlgebra
using Test

@testset "Quadratic objective" begin
    Q = Diagonal([1e1, 1e-1, 1e0, 1e0, 1e0, 1e0, 1e0, 1e0, 1e0, 1e0])
    q = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]

    objective = x -> 0.5*dot(x, Q*x) + dot(x, q)

    x_star = -Q\q
    f_star = objective(x_star)
    L = opnorm(Q)

    w0 = zeros(10)

    @testset "$(name)" for (name, optimizer) in [
        "GradientDescent" => ProtoGrad.GradientDescent(stepsize=1/L),
    ]
        for (k, output) in enumerate(take(optimizer(w0, objective), 200))
            @test objective(output.solution) - f_star <= 2 * L * norm(x_star)^2 / (2 * (k - 1))
        end
    end

    @testset "$(name)" for (name, optimizer) in [
        "Nesterov" => ProtoGrad.Nesterov(stepsize=1/L),
    ]
        for (k, output) in enumerate(take(optimizer(w0, objective), 200))
            @test objective(output.solution) - f_star <= 2 * L * norm(x_star)^2 / k^2
        end
    end

    @testset "$(name)" for (name, optimizer) in [
        "Polyak" => ProtoGrad.Polyak(stepsize=1/L, momentum=0.5),
        "RMSProp" => ProtoGrad.RMSProp(stepsize=1/L, alpha=0.99, epsilon=1e-8),
        "Adam" => ProtoGrad.Adam(stepsize=1/L, beta1=0.9, beta2=0.999, epsilon=1e-8),
    ]
        max_iter = 1000

        iterations = take(optimizer(w0, objective), max_iter)
        solution = ProtoGrad.last(iterations).solution

        @test objective(solution) <= f_star + 0.01 * abs(f_star)
    end
end

@testset "Supervised objective" begin

    @testset "Linear model ($(T))" for T in [Float32, Float64]

        input_size, output_size = 10, 1
        batch_size = 5000

        W_true = randn(T, output_size, input_size)
        b_true = randn(T, output_size)

        x = randn(T, input_size, batch_size)
        y = W_true * x .+ b_true + randn(T, output_size, batch_size) / 10

        data_iter = Iterators.repeated((x, y))
        f = SupervisedObjective(mse, data_iter)
        Lf = 2 * opnorm(x * x') / batch_size
        stepsize = 1 / Lf

        m = Linear(randn(T, output_size, input_size), randn(T, output_size))

        @testset "Basic checks $(name)" for (name, optimizer) in [
            "GradientDescent" => ProtoGrad.GradientDescent(stepsize=stepsize),
            "Polyak" => ProtoGrad.Polyak(stepsize=stepsize, momentum=0.5),
            "Nesterov" => ProtoGrad.Nesterov(stepsize=stepsize),
            "BarzilaiBorwein" => ProtoGrad.BarzilaiBorwein(alpha=stepsize),
            "AdaGrad" => ProtoGrad.AdaGrad(stepsize=stepsize),
            "RMSProp" => ProtoGrad.RMSProp(stepsize=stepsize),
            "AdaDelta" => ProtoGrad.AdaDelta(),
            "Adam" => ProtoGrad.Adam(stepsize=stepsize),
        ]
            m0 = copy(m)
            seq = optimizer(m, f)

            @test Base.IteratorSize(typeof(seq)) == Base.IsInfinite()

            for output in Iterators.take(seq, 5)
                @test length(ProtoGrad.overlap(m, output.solution)) == 0
                @test typeof(output.solution) == typeof(m)
            end

            @test vec(m) == vec(m0)
        end

        @testset "Accuracy $(name)" for (name, optimizer) in [
            "GradientDescent" => ProtoGrad.GradientDescent(stepsize=stepsize),
            "Polyak" => ProtoGrad.Polyak(stepsize=stepsize, momentum=0.5),
            "Nesterov" => ProtoGrad.Nesterov(stepsize=stepsize),
        ]
            seq = optimizer(m, f)

            for (it, output) in enumerate(seq)
                if it >= 1000
                    @test isapprox(output.solution.W, W_true, rtol=5e-2)
                    # @test isapprox(output.solution.b, b_true, rtol=1e-1)
                    break
                end
            end
        end
    end
end


text:
simple = "----. - -. ..--.. ----- ----. ...- ---..  --.. ..--.. --..-- --.- .... -..- -.- ..-  -.--. .--. -.. ..--.-  ..--.- .-. -.-.-. -.. .-.-.- -..- -.-.--  --... .--- -..-. .--.-. .-.. ..... .- -.-- .. -"

function decodemorse(morsecode)
  # your code here
  tmp = replace(morsecode,r"\+"=>" + ")
  arr = split(tmp,' ')
  filter!(x -> x != "",arr)
    function rpp(x)
        if x == "+"
            return " "
        else
            return x
        end
    end

#   res = map(x -> MOR)
  return arr
end

res = decodemorse(simple)

println(res)

max_stars_repo_path: kata-julia/Decode the Morse code.jl

max_stars_repo_name: crimsonskyrem/codewars

max_issues_repo_path: kata-julia/Decode the Morse code.jl

max_issues_repo_name: crimsonskyrem/codewars

max_forks_repo_path: kata-julia/Decode the Morse code.jl

max_forks_repo_name: crimsonskyrem/codewars


text:
"""
    Bus

A power-system bus.

# Constructor
```julia
Bus(number, name, bustype, angle, voltage, voltagelimits, basevoltage)
```

# Arguments
* `number`::Int64 : number associated with the bus
* `name`::String : the name of the bus
* `bustype`::String : type of bus, [PV, PQ, SF]; may be `nothing`
* `angle`::Float64 : angle of the bus in degrees; may be `nothing`
* `voltage`::Float64 : voltage as a multiple of basevoltage; may be `nothing`
* `voltagelimits`::NamedTuple(min::Float64, max::Float64) : limits on the voltage variation as multiples of basevoltage; may be `nothing`
* `basevoltage`::Float64 : the base voltage in kV; may be `nothing`

"""
struct Bus <: Injection
    # field docstrings work here! (they are not for System)
    """ number associated with the bus """
    number::Int64
    """ the name of the bus """
    name::String
    """ bus type, [PV, PQ, SF] """
    bustype::Union{String,Nothing} # [PV, PQ, SF]
    """ angle of the bus in degrees """
    angle::Union{Float64,Nothing} # [degrees]
    """ voltage as a multiple of basevoltage """
    voltage::Union{Float64,Nothing} # [pu]
    """ limits on the voltage variation as multiples of basevoltage """
    voltagelimits::Union{NamedTuple{(:min, :max),Tuple{Float64,Float64}},
                         Nothing} # [pu]
    """
    the base voltage in kV
    """
    basevoltage::Union{Float64,Nothing} # [kV]
end

# DOCTODO  add this constructor type to docstring for Bus
Bus(;   number = 0,
        name = "init",
        bustype = nothing,
        angle = 0.0,
        voltage = 0.0,
        voltagelimits = (min = 0.0, max = 0.0),
        basevoltage = nothing
    ) = Bus(number, name, bustype, angle, voltage,
            orderedlimits(voltagelimits, "Voltage"), basevoltage)

# DOCTODO What are LoadZones? JJS 1/18/19
struct LoadZones  <: Injection
    number::Int
    name::String
    buses::Array{Bus,1}
    maxactivepower::Float64
    maxreactivepower::Float64
end

LoadZones(;   number = 0,
        name = "init",
        buses = [Bus()],
        maxactivepower = 0.0,
        maxreactivepower = 0.0
    ) = LoadZones(number, name, buses, maxactivepower, maxreactivepower)


text:
#===============================================================================
# K-Nearest-Neighbors for the Numerai dataset
#
# Author: Carl Cortright
# Date: 1/5/2017
#
===============================================================================#
using ScikitLearn
using DataFrames, DataArrays

@sk_import neighbors: KNeighborsClassifier
@sk_import preprocessing: (LabelBinarizer, StandardScaler)

# Define he training set
train = readtable("/home/carl/Documents/Numerai/Datasets/Jan92017/numerai_training_data.csv")

# Pull out features and targets
features = convert(Array, train[1:50])
targets = convert(Array, train[51])

# Fit a KNN
knn = KNeighborsClassifier(n_neighbors=3, algorithm="ball_tree")
ScikitLearn.fit!(knn, features, targets)

println("before prediction")

for i in 1:length(features[:, 1])
  println(ScikitLearn.predict(knn, features[i, 1:50]) == targets[i])
end
# accuracy = sum() / length(targets)
# println("accuracy: $accuracy")

max_stars_repo_name: carlcortright/Numerai

max_issues_repo_name: carlcortright/Numerai

max_forks_repo_name: carlcortright/Numerai


text:
# ---
# title: 309. Best Time to Buy and Sell Stock with Cooldown
# id: problem309
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/best-time-to-buy-and-sell-stock-with-cooldown/description/>
# hidden: true
# ---
#
# Say you have an array for which the _i_ th element is the price of a given
# stock on day _i_.
#
# Design an algorithm to find the maximum profit. You may complete as many
# transactions as you like (ie, buy one and sell one share of the stock multiple
# times) with the following restrictions:
#
#   * You may not engage in multiple transactions at the same time (ie, you must sell the stock before you buy again).
#   * After you sell your stock, you cannot buy stock on next day. (ie, cooldown 1 day)
#
# **Example:**
#
#
#
#     Input: [1,2,3,0,2]
#     Output: 3
#     Explanation: transactions = [buy, sell, cooldown, buy, sell]
#
#
#
## @lc code=start
using LeetCode

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


text:
const PostcentralGyrus = BilateralStructure("postcentral gyrus")

const SuperiorParietalLobule = BilateralStructure("superior parietal lobule")

const InferiorParietalLobule = BilateralStructure("inferior parietal lobule")

const SupramarginalGyrus = BilateralStructure("supramarginal gyrus")

const AngularGyrus = BilateralStructure("angular gyrus")

const Precuneus = BilateralStructure("precuneus")

const ParietalLobe = BilateralStructure("parietal lobe")


text:
#!/usr/bin/env julia

using Luxor, Colors

using Test

using Random
Random.seed!(2)

function meshtest1(fname)
    Drawing(1200, 1200, fname)
    origin()
    background("ivory")
    tiles = Partition(1200, 1200, 120, 120)
    fontsize(120)
    fontface("Courier-Bold")
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            setcolor("black")
            setline(1)
            box(O, tiles.tilewidth-8, tiles.tileheight-8, :stroke)
            circle(O, 20, :fill)
            pl = ngon(O, tiles.tilewidth/2, 4, pi/4, vertices=true)
            # some of these should be transparent
            mesh1 = mesh(pl, [
                "red",
                Colors.RGBA(randomcolor()...),
                Colors.RGBA(randomcolor()...),
                Colors.RGBA(randomcolor()...)
                ])
            setmesh(mesh1)
            # box(O, tiles.tilewidth-10, tiles.tileheight-10, :fillpreserve)
            text(string(Char(rand(64:89))), halign=:center, valign=:middle)
        end
    end
    @test finish() == true
end

function meshtest2(fname)
    Drawing(1200, 1200, fname)
    origin()
    background("ivory")
    setcolor("black")
    fontsize(60)
    text("illuminati", halign=:center)
    pl = ngon(O, 300, 3, pi/6, vertices=true)
    mesh1 = mesh(pl, [
        "purple",
        Colors.RGBA(randomcolor()...),
        "yellow"
        ])
    setmesh(mesh1)
    setline(260)
    ngon(O, 300, 3, pi/6, :strokepreserve)
    setline(5)
    sethue("black")
    strokepath()
    @test finish() == true
end

function meshtest3(fname)
    Drawing(1200, 1200, fname)
    origin()
    background("ivory")
    setcolor("black")
    fontsize(60)
    text("illuminati", halign=:center)
    bp = makebezierpath(ngon(O, 250, 3, pi/6, vertices=true))
    mesh1 = mesh(bp, [
        "purple",
        Colors.RGBA(randomcolor()...),
        "yellow"
        ])
    setmesh(mesh1)
    setline(260)
    drawbezierpath(bp, :strokepreserve)

    setline(5)
    sethue("black")
    strokepath()
    @test finish() == true
end

fname = "meshtest1.png"
meshtest1(fname)
println("...finished mesh test 1: output in $(fname)")

fname = "meshtest2.png"
meshtest2(fname)
println("...finished mesh test 2: output in $(fname)")

fname = "meshtest3.png"
meshtest3(fname)
println("...finished mesh test 3: output in $(fname)")


text:
function plotconfig(config, dx=10 ; kwargs...)
    ex = extrema(config.x)
    bins = first(ex):dx:last(ex)
    stephist(vec(config.x), bins=bins, group=repeat(config.obs.class, inner=2), leg=:outertopright; kwargs...)
end


text:
module BenchForeachSeq

using BenchmarkTools
using ThreadsX.Implementations: foreach_cartesian_seq

const floatsink = Ref(0.0)
@inline consume(x) = floatsink[] = ifelse(x == 0, x, floatsink[])

function foreachc(f::F, xs) where {F}
    foreach_cartesian_seq(f, xs, Val(false))
end

n = 1000
v = ones(n * n)
m = ones(n, n)
suite = BenchmarkGroup()

for (label, _foreach) in [("base", foreach), ("tx", foreachc)]
    s = suite[label] = BenchmarkGroup()
    s["Vector"] = @benchmarkable($_foreach(consume, $v))
    s["Matrix"] = @benchmarkable($_foreach(consume, $m))
    s["Transpose"] = @benchmarkable($_foreach(consume, transpose($m)))
end

end  # module
BenchForeachSeq.suite


text:
@testset "$TEST $G" begin

f, fio = mktemp()

g = G(10,20)
gpath = joinpath(testdir,"testdata","gtest.gt.gz")
writegraph(gpath, g)
h = readgraph(gpath, :gt, G, compressed=true)
@test h == g
h = readgraph(gpath, G)
@test h == g
if G <: ANetwork
    vprop!(g, "w", rand(nv(g)))
    writenetwork(gpath, g)
    h = readnetwork(gpath, :gt, G, compressed=true)
    test_networks_eq(h, g)
    h = readnetwork(gpath, G)
    test_networks_eq(h, g)
end
rm(gpath)

g = DG(10,20)
gpath = joinpath(testdir,"testdata","gtest.gt.gz")
writegraph(gpath, g)
h = readgraph(gpath, :gt, G, compressed=true)
@test h == g
h = readgraph(gpath, G)
@test h == g
if G <: ANetwork
    vprop!(g, "w", rand(nv(g)))
    eprop!(g, "w", EdgeMap(g, rand(ne(g)))) #TODO

    writenetwork(gpath, g)
    h = readnetwork(gpath, :gt, G, compressed=true)
    test_networks_eq(h, g)
    h = readnetwork(gpath, G)
    test_networks_eq(h, g)
end
rm(gpath)

g = readgraph(:lesmis, G)
@test typeof(g) == G
@test nv(g) == 77
@test ne(g) == 254

writegraph(f, g, :gt)
h = readgraph(f, :gt, G)
@test g == h

g = readgraph(:serengetifoodweb, G)
@test typeof(g) == DG
@test nv(g) == 161
@test ne(g) == 592

writegraph(f, g, :gt)
h = readgraph(f, :gt, G)
@test g == h

g = G(10,0)
writegraph(f, g, :gt)
ga = readgraph(f, :gt, G)
@test g == ga

g = G(10,20)
writegraph(f, g, :gt)
ga = readgraph(f, :gt, G)
@test g == ga

g = DG(10,20)
writegraph(f, g, :gt)
ga = readgraph(f, :gt, G)
@test g == ga

if G <: ANetwork ##############################

g = readnetwork(:karate, G)
@test sort(gprop_names(g)) == ["description","readme"]
@test sort(vprop_names(g)) == ["pos"]
@test sort(eprop_names(g)) == []

writenetwork(f, g, :gt)
h = readnetwork(f, :gt, G)
@test h == g
@test h.props == g.props

g = readnetwork(:lesmis)
@test sort(gprop_names(g)) == ["description","readme"]
@test sort(vprop_names(g)) == ["label","pos"]
@test sort(eprop_names(g)) == ["value"]

writenetwork(f, g, :gt)
h = readnetwork(f, :gt, G)
@test h == g
# @test h.props == g.props #dataset graphs' edges index can change
writenetwork(f, h, :gt)
h2 = readnetwork(f, :gt, G)
@test h == h2
@test h.props == h2.props

end #if network ####################

isfile(f) && rm(f)
 end #testset


text:
"""
$SIGNATURES

Return nothing if the number of rows of `X` and `y` match and throws a
`DimensionMismatch` error otherwise.
"""
function check_nrows(X::AbstractMatrix, y::AbstractVecOrMat)::Nothing
    size(X, 1) == size(y, 1) && return nothing
    throw(DimensionMismatch("`X` and `y` must have the same number of rows."))
end

"""
$SIGNATURES

Throws an error if the argument (e.g. penalty scaling) is negative.
"""
check_pos(λ::Real) = λ ≥ 0 || throw(ArgumentError("Penalty scaling should be positive."))

"""
$SIGNATURES

Given a matrix `X`, append a column of ones if `fit_intercept` is true.
"""
function augment_X(X::AbstractMatrix{<:Real}, fit_intercept::Bool)
    fit_intercept || return X
    return hcat(X, ones(eltype(X), size(X, 1)))
end

"""
$SIGNATURES

Return `X*θ` if `c==0` (default) or `X*P` where `P=reshape(θ, size(X, 2), p)`
in the multi-class case.
"""
function apply_X(X, θ, c=0)
    p = size(X, 2)
    if c == 0
        length(θ) == p || return X * θ[1:p] .+ θ[end]
        return X * θ
    else
        noβ = length(θ) == p * c
        W = reshape(θ, p + Int(!noβ), c)
        noβ || return X * view(W, 1:p, :) .+ view(W, p+1, :)'
        return X * W
    end
end

"""
$SIGNATURES

In-place application of X*θ.
"""
function apply_X!(Xθ, X, θ, c=1, scratch=nothing)
    p = size(X, 2)
    if c == 1
        if length(θ) == p
            mul!(Xθ, X, θ)
        else
            mul!(Xθ, X, view(θ, 1:p))
            Xθ .+= θ[end]
        end
    else
        noβ = length(θ) == p * c
        W   = scratch.pc
        copyto!(W, reshape(θ, p + Int(!noβ), c))
        if noβ
            mul!(Xθ, X, W)
        else
            mul!(Xθ, X, view(W, 1:p, :))
            Xθ .+= view(W, p+1, :)'
        end
    end
end

"""
$SIGNATURES

In-place application of X'*z (only for regression case).
"""
function apply_Xt!(Xtv, X, z)
    p  = size(X, 2)
    p_ = length(Xtv)
    if p == p_
        mul!(Xtv, X', z)
    else
        if Xtv isa SubArray
            mul!(Xtv, X', z)
            Xtv.parent[end] = sum(z)
        else
            mul!(view(Xtv, 1:p), X', z)
            Xtv[end] = sum(z)
        end
    end
end

"""
$SIGNATURES

Form (X'X) while being memory aware (assuming p ≪ n).
"""
function form_XtX(X, fit_intercept, lambda = 0, penalize_intercept = true)
    if fit_intercept
        n, p = size(X)
        XtX  = zeros(p+1, p+1)
        Xt1  = sum(X, dims=1)
        mul!(view(XtX, 1:p, 1:p), X', X) # O(np²)
        @inbounds for i in 1:p
            XtX[i, end] = XtX[end, i] = Xt1[i]
        end
        XtX[end, end] = n
    else
        XtX = X'*X # O(np²)
    end
    if !iszero(lambda)
        λ = convert(eltype(XtX), lambda)
        @inbounds for i in 1:size(XtX, 1) + fit_intercept * (penalize_intercept - 1)
            XtX[i,i] += λ
        end
    end
    return Hermitian(XtX)
end

# Sigmoid and log-sigmoid

SIGMOID_THRESH(T::Type{<:AbstractFloat}) = log(1 / eps(T) - 1)

"""
$SIGNATURES

Return the sigmoid computed in a numerically stable way:

``σ(x) = 1/(1+exp(-x))``
"""
function sigmoid(x::T) where T <: AbstractFloat
    τ = SIGMOID_THRESH(T)
    x > τ  && return one(T)
    x < -τ && return zero(T)
    return 1 / (1 + exp(-x))
end
sigmoid(x) = sigmoid(float(x))
σ = sigmoid

"""
$SIGNATURES

Return the log sigmoid computed in a numerically stable way:

``logσ(x) = -log(1+exp(-x)) = log(exp(x)/(exp(x) + 1)) = x - log(1+exp(x))``
"""
function logsigmoid(x::T) where T <: AbstractFloat
    τ = SIGMOID_THRESH(T)
    x > τ  && return zero(T)
    x < -τ && return x
    return -log1p(exp(-x))
end
logsigmoid(x) = logsigmoid(float(x))
logσ = logsigmoid

"""
$SIGNATURES

Return the softmax computed in a numerically stable way:

``σ(x) = exp.(x) ./ sum(exp.(x))``

Implementation taken from NNlib.jl.
"""
function softmax(X::AbstractMatrix{<:Real})
    max_ = maximum(X, dims=2)
    exp_ = exp.(X .- max_)
    return exp_ ./ sum(exp_, dims=2)
end

"""
$SIGNATURES

In place computation of `H = H + λI` where  `H` is a square matrix.
"""
function add_λI!(H::Matrix, λ::Real, penalize_intercept::Bool=true)
    λ = convert(eltype(H), λ)
    @inbounds for i in 1:size(H, 1)-1
        H[i,i] += λ
    end
    H[end, end] += ifelse(penalize_intercept, λ, zero(eltype(H)))
end

"""
$SIGNATURES

Soft-thresholding S_η(z).
"""
soft_thresh(z, η) = sign(z) * max(abs(z) - η, 0)

"""
$SIGNATURES

Threshold the number if its absolute value is too close to zero.
"""
clip(z, τ) = ifelse(abs(z) < τ, τ, z)

"""
$SIGNATURES

Return λ if penalize intercept otherwise 0, useful in computations of Hessian.
"""
λ_if_penalize_intercept(glr, λ) = ifelse(glr.penalize_intercept, λ, zero(λ))

"""
$SIGNATURES

Return a view of θ if the last element should not be penalized.
"""
@inline function view_θ(glr, θ)
    f = glr.fit_intercept && !glr.penalize_intercept
    f && return view(θ, 1:length(θ)-1)
    θ
end


text:
using Documenter, JLSO

makedocs(
    modules=[JLSO],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
    ],
    repo="https://github.com/invenia/JLSO.jl/blob/{commit}{path}#L{line}",
    sitename="JLSO.jl",
    authors="Rory Finnegan",
    assets=[
        "assets/invenia.css",
    ],
    strict = true,
    checkdocs = :none,
)

deploydocs(;
    repo="github.com/invenia/JLSO.jl",
)


text:
using LinearAlgebra
using StaticArrays

function max_distance_for_pbc_dist(traj::Trajectory)
    max_distance_for_pbc_dist(traj.mdbox)
end

function max_distance_for_pbc_dist(mdbox::MDBox)
    max_distance_for_pbc_dist(mdbox.pbc_matrix)
end

function max_distance_for_pbc_dist(pbc::AbstractArray)
    volume = det(pbc)
    a, b, c = eachcol(pbc)
    dist1 = volume / norm(cross(a, b))
    dist2 = volume / norm(cross(b, c))
    dist3 = volume / norm(cross(c, a))
    return min(dist1, dist2, dist3) / 2
end

function pbc_dist(point1, point2, mdbox::OrthorhombicBox)
    for i in 1:3
        mdbox.dist_tmp[i] = (point1[i] - point2[i]) / mdbox.cell_parameters[i]
        mdbox.dist_tmp[i] -= floor(mdbox.dist_tmp[i] + 0.5)
        mdbox.dist_tmp[i] *= mdbox.cell_parameters[i]
    end
    return norm(mdbox.dist_tmp)
end

#function pbc_dist(point1, point2, mdbox::TriclinicBox)
#    pbc_dist(point1, point2, mdbox.pbc_matrix, mdbox.inv_pbc_matrix, mdbox.realspace_tmp, mdbox.inversespace_tmp)
#end

#function pbc_dist(point1, point2, pbc, inv_pbc = inv(pbc), dist_tmp = zeros(MVector{3}), matmul_tmp = zeros(MVector{3}))
#    dist_tmp .= point1 .- point2
#    mul!(matmul_tmp, inv_pbc,dist_tmp)
#    matmul_tmp .-= round.(matmul_tmp)
#    mul!(dist_tmp, pbc, matmul_tmp)
#    return norm(dist_tmp)
#end

# ugliest implementation yet, but also fast
# who needs arrays anyways
function pbc_dist(p1, p2, pbc, inv_pbc)
    # pbc (and inv_pbc by extension) have to be upper-triagonal now!
    # direct distance between the two points
        # d = p1 .- p300
    d1, d2, d3 = p1[1] - p2[1], p1[2] - p2[2], p1[3] - p2[3]
    # convert distance to relative coordinates and save rounded values
        # m = round.(inv_pbc * d)
    m1 = round(inv_pbc[1, 1] * d1 + inv_pbc[1, 2] * d2 + inv_pbc[1, 3] * d3)
    m2 = round(inv_pbc[2, 2] * d2 + inv_pbc[2, 3] * d3)
    m3 = round(inv_pbc[3, 3] * d3)
    # convert rounded relative distance into real distance and subtract it from the real distance; return norm
        # d -= pbc * m
        # n = norm(d)
    n = sqrt(
        (d1 - pbc[1, 1] * m1 - pbc[1, 2] * m2 - pbc[1, 3] * m3)^2 +
        (d2 - pbc[2, 2] * m2 - pbc[2, 3] * m3)^2 +
        (d3 - pbc[3, 3] * m3)^2
        )
    return n
end

function pbc_dist(p1, p2, mdbox::TriclinicBox)
    pbc_dist(p1, p2, mdbox.pbc_matrix, mdbox.inv_pbc_matrix)
end

function minimum_image_vector(point1, point2, mdbox::TriclinicBox)
    minimum_image_vector(point1, point2, mdbox.pbc_matrix, mdbox.inv_pbc_matrix, mdbox.realspace_tmp, mdbox.inversespace_tmp)
end

function minimum_image_vector(point1, point2, pbc, inv_pbc = inv(pbc), dist_tmp = zeros(MVector{3}), matmul_tmp = zeros(MVector{3}))
    dist_tmp .= point1 .- point2
    mul!(matmul_tmp, inv_pbc,dist_tmp)
    matmul_tmp .-= floor.(matmul_tmp .+ 0.5)
    mul!(dist_tmp, pbc, matmul_tmp)
    return dist_tmp
end

function pbc_dist_triclinic(point1, point2, pbc, inv_pbc = inv(pbc), dist_tmp = zeros(MVector{3}), matmul_tmp = zeros(MVector{3}))
    dist_tmp .= abs.(point1 .- point2)
    mul!(matmul_tmp, inv_pbc,dist_tmp)
    matmul_tmp .-= floor.(matmul_tmp .+ 0.5) .+ 1
    distance = pbc[1, 1] + pbc[2, 2] + pbc[3, 3]
    for i in 1:3
        for j in 1:3
            for k in 1:3
                mul!(dist_tmp, pbc, matmul_tmp)
                new_distance = norm(dist_tmp)
                if new_distance < distance
                    distance = new_distance
                end
                matmul_tmp[3] += 1
            end
            matmul_tmp[3] -= 3
            matmul_tmp[2] += 1
        end
        matmul_tmp[2] -= 3
        matmul_tmp[1] += 1
    end
    return distance
end

function next_neighbor(point1, group2, mdbox::MDBox)
    point2 = @view group2[:, 1]
    index, distance = 1, pbc_dist(point1, point2, mdbox)
    for i in 2:size(group2)[2]
        point2 = @view group2[:, i]
        new_distance = pbc_dist(point1, point2, mdbox)
        if new_distance < distance
            distance = new_distance
            index = i
        end
    end
    return index, distance
end


text:
immutable Gamma <: ContinuousUnivariateDistribution
    α::Float64
    β::Float64

    function Gamma(α::Real, β::Real)
        α > zero(α) && β > zero(β) || error("Gamma: both shape and scale must be positive")
        new(float64(α), float64(β))
    end

    Gamma(α::Real) = Gamma(α, 1.0)
    Gamma() = new(1.0, 1.0)
end

@_jl_dist_2p Gamma gamma

@distr_support Gamma 0.0 Inf

#### Parameters

shape(d::Gamma) = d.α
scale(d::Gamma) = d.β
rate(d::Gamma) = 1.0 / d.β

params(d::Gamma) = (d.α, d.β)

#### Statistics

mean(d::Gamma) = d.α * d.β

var(d::Gamma) = d.α * d.β^2

skewness(d::Gamma) = 2.0 / sqrt(d.α)

kurtosis(d::Gamma) = 6.0 / d.α

function mode(d::Gamma)
    (α, β) = params(d)
    α >= 1.0 ? β * (α - 1.0) : error("Gamma has no mode when shape < 1.0")
end

function entropy(d::Gamma)
    (α, β) = params(d)
    α + lgamma(α) + (1.0 - α) * digamma(α) + log(β)
end

mgf(d::Gamma, t::Real) = (1.0 - t * d.β)^(-d.α)

cf(d::Gamma, t::Real) = (1.0 - im * t * d.β)^(-d.α)

#### Evaluation

gradlogpdf(d::Gamma, x::Float64) =
    insupport(Gamma, x) ? (d.α - 1.0) / x - 1.0 / d.β : 0.0

#### Fit model

immutable GammaStats <: SufficientStats
    sx::Float64      # (weighted) sum of x
    slogx::Float64   # (weighted) sum of log(x)
    tw::Float64      # total sample weight

    GammaStats(sx::Real, slogx::Real, tw::Real) = new(float64(sx), float64(slogx), float64(tw))
end

function suffstats(::Type{Gamma}, x::Array)
    sx = 0.
    slogx = 0.
    for xi = x
        sx += xi
        slogx += log(xi)
    end
    GammaStats(sx, slogx, length(x))
end

function suffstats(::Type{Gamma}, x::Array, w::Array{Float64})
    n = length(x)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end

    sx = 0.
    slogx = 0.
    tw = 0.
    for i = 1:n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        sx += wi * xi
        slogx += wi * log(xi)
        tw += wi
    end
    GammaStats(sx, slogx, tw)
end

function gamma_mle_update(logmx::Float64, mlogx::Float64, a::Float64)
    ia = 1.0 / a
    z = ia + (mlogx - logmx + log(a) - digamma(a)) / (abs2(a) * (ia - trigamma(a)))
    1.0 / z
end

function fit_mle(::Type{Gamma}, ss::GammaStats;
    alpha0::Float64=NaN, maxiter::Int=1000, tol::Float64=1.0e-16)

    mx = ss.sx / ss.tw
    logmx = log(mx)
    mlogx = ss.slogx / ss.tw

    a::Float64 = isnan(alpha0) ? 0.5 / (logmx - mlogx) : alpha0
    converged = false

    t = 0
    while !converged && t < maxiter
        t += 1
        a_old = a
        a = gamma_mle_update(logmx, mlogx, a)
        converged = abs(a - a_old) <= tol
    end

    Gamma(a, mx / a)
end


text:
using WAVI
function MISMIP_PLUS_Ice1r_2km()
    #Grid and boundary conditions
    nx = 320
    ny = 40
    nσ = 4
    x0 = 0.0
    y0 = -40000.0
    dx = 2000.0
    dy = 2000.0
    h_mask=trues(nx,ny)
    u_iszero = falses(nx+1,ny); u_iszero[1,:].=true
    v_iszero=falses(nx,ny+1); v_iszero[:,1].=true; v_iszero[:,end].=true
    grid = Grid(nx = nx,
                ny = ny,
                nσ = nσ,
                x0 = x0,
                y0 = y0,
                dx = dx,
                dy = dy,
                h_mask = h_mask,
                u_iszero = u_iszero,
                v_iszero = v_iszero)

    #Bed
    bed = WAVI.mismip_plus_bed #function definition of bed

    #Initial conditions from Ice0
    h_init = Array{Float64}(undef, nx, ny);
    read!("examples\\MISMIP_ice0_final.bin", h_init)
    h_init = ntoh.(h_init);
    initial_conditions = InitialConditions(initial_thickness = h_init)

    #solver parameters
    maxiter_picard = 30
    solver_params = SolverParams(maxiter_picard = maxiter_picard)

    #Physical parameters
    accumulation_rate = 0.3
    params = Params(accumulation_rate = accumulation_rate)

    #make the model
    model = Model(grid = grid,
                     bed_elevation = bed,
                     params = params,
                     solver_params = solver_params,
                     initial_conditions = initial_conditions)

    #embed the ice model with melt rate model
    function m1(h, b)
        draft = -(918.0 / 1028.0) * h
        cavity_thickness = draft .- b
        cavity_thickness = max.(cavity_thickness, 0)
        m =  0.2*tanh.(cavity_thickness./75).*max.((-100 .- draft), 0)
        return m
    end
    arguments = (h = model.fields.gh.h, b = model.fields.gh.b)
    melt_rate_model = AnalyticMeltRate(melt_partial_cell = true,
                                    melt_rate_function = m1,
                                    function_arguments = arguments)
    add_melt_rate_model!(model,melt_rate_model)

    #timestepping parameters
    niter0 = 0
    dt = 0.001
    end_time = 1.
    chkpt_freq = 0.1
    pchkpt_freq = 0.1
    timestepping_params = TimesteppingParams(niter0 = niter0,
                                            dt = dt,
                                            end_time = end_time,
                                            chkpt_freq = chkpt_freq,
                                            pchkpt_freq = pchkpt_freq)

    #output parameters
    folder = "outputs"
    isdir(folder) && rm(folder, force = true, recursive = true)
    mkdir(folder) #make a clean folder for outputs
    outputs = (h   = model.fields.gh.h,
                u  = model.fields.gh.u,
                v  = model.fields.gh.v,
                melt = model.fields.gh.basal_melt,
                grfrac = model.fields.gh.grounded_fraction)
    output_freq = 0.1
    output_params = OutputParams(outputs = outputs,
                            output_freq = output_freq,
                            output_format = "mat",
                            output_path = folder,
                            zip_format = "nc")

    simulation = Simulation(model = model,
                        timestepping_params = timestepping_params,
                        output_params = output_params)

    #perform the simulation
    run_simulation!(simulation)
    return simulation
end

@time simulation = MISMIP_PLUS_Ice1r_2km();


text:
using CasesLH
using Test

@testset "CasesLH.jl" begin
    include("modifier_test.jl");
    include("case_spec_test.jl");
end


text:
# EFTfitter.jl - Advanced Tutorial
# This tutorial introduces some of the advanced functionalities of EFTfitter.jl
# using a generic example. Please see the [tutorial] for basic usage of EFTfitter.
using EFTfitter
using BAT            # for sampling
using IntervalSets   # for specifying the prior
using Distributions  # for specifying the prior
using Plots # for plotting

include("advanced_inputs.jl")

# We need to modify the definition of the `EFTfitterModel` by also passing the `nuisance_correlations`:
model = EFTfitterModel(parameters, measurements, correlations, nuisance_correlations)

posterior = PosteriorDensity(model)

algorithm = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)
samples = bat_sample(posterior, algorithm).result;

p = plot(samples)
savefig(p, "plot.pdf")

# Warning:
# When using nuisance correlations, it may happen that the fit does not converge anymore.
# This happens mainly when the correlation values are close to -1 or +1.
# In such a case, you can try to reduce the allowed ranges in the priors for the Nuisance
# correlations to investigate at which values this happens.

# ============= Ranking of measurements and uncertainties ========================#
# With the `rank_measurements` and `rank_uncertainties` functions, the influence of
# the individual measurements or uncertainty types on the result of a fit can be estimated.
# For the ranking, each active measurement (respectively uncertainty type) is deactivated
# at a time and the fit is repeated. The results of all fits with one deactivated
# measurement (or uncertainty type) are then compared to the fit result with all measurements (uncertainty types) activated.
# A ranking is calculated based on a ranking criterion calculated from the posterior distributions of these fits.

# The default ranking criterion is the relative increase of the total width of the
# smallest interval containing 90% of the posterior probability when deactivating a measurement.
# For models with more than one parameter, the sum of relative increases of all
# one-dimensional smallest intervals is used. The default ranking criterion is therefore
#`SumOfSmallestIntervals(p=0.9, bins=200)`.
measurement_ranking = EFTfitter.rank_measurements(model)

# For ranking the uncertainty types, the relative decrease is used.
uncertainty_ranking = EFTfitter.rank_uncertainties(model, criterion = SumOfSmallestIntervals(p=0.9, bins=200))

# Please see the [ranking documentation] for further ranking criteria and keyword arguments.

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
function mdl_DCside_ode!( device::DynInverter{C,O,VC,FixedDCSource,P,F}) where {C <: Converter,
                                                   O <: OuterControl,
                                                   VC <: VSControl,
                                                   P <: FrequencyEstimator,
                                                   F <: Filter}

    #Obtain external states inputs for component

    #Obtain inner variables for component

    #Get parameters
    Vdc = device.dc_source.voltage

    #Obtain indices for component w/r to device

    #Define internal states for converter

    #Compute states... none?

    #Update inner_vars
    device.inner_vars[Vdc_var] = Vdc
end


text:
using MosimoBase: Vector3s, Vector3

import ..Constraints: FixedDistanceConstraint
import ..Types: AbstractIntegrator, current_positions, current_velocities, set_state!

# The constraints can only be of the type that require fixed distances.
struct RattleIntegrator{
    M <: Function, F <: Function,
    D <: Function
} <: AbstractIntegrator{M, F}
    positions::Vector3s
    velocities::Vector3s
    timestep::Float64

    mass_function::M
    force_function::F
    distance_function::D

    tolerance::Float64
    max_iterations::Int
    constraints::Vector{FixedDistanceConstraint}

    forces::Vector3s
end

# r_ij = r_i - r_j
_default_distance() = (a::Vector3, b::Vector3) -> a - b

function RattleIntegrator(
    r, v, dt, m, f, ξ, max_iter,
    gs = FixedDistanceConstraint[];
    distance = _default_distance()
)
    N = length(r)
    @assert N === length(v)
    fs = f(r)
    @assert N === length(fs)
    RattleIntegrator{typeof(m), typeof(f), typeof(distance)}(
        r, v, dt, m, f, distance, ξ, max_iter, gs, fs
    )
end

current_positions(ri::RattleIntegrator) = ri.positions
current_velocities(ri::RattleIntegrator) = ri.velocities

function set_state!(ri::RattleIntegrator, rs::Vector3s, vs::Vector3s)
    ri.positions .= rs
    ri.velocities .= vs
    ri
end


text:
import Base.dot

"""An array of random variables (and also a random variable itself)
  `T` is the range type of elements (e.g for multivariate normal, T = Float64)
  `N` is the dimensionality of array"""
RandArray{T, N} = Array{T, N} where {T<:RandVar}
Smelly{Q} = Array{T, 1} where {T<:RandVar{Q}} where Q
ExecutableRandArray{T, N} = Array{T, N} where {T<:JuliaRandVar}

dot(a::RandVar, b::RandVar) = a * b
dot(a::Number, b::RandVar) = a * b
dot(a::RandVar, b::Number) = a * b

convert{T}(::Type{ExecutableRandArray}, xs::Array{RandVar{T}}) =
  (rv -> convert(JuliaRandVar{T}, rv)).(xs)

(xs::ExecutableRandArray)(ω::Omega) = (x->x(ω)).(xs)

rand(xs::ExecutableRandArray) = xs(LazyRandomVector(Float64))
rand(xs::RandArray) = rand(convert(ExecutableRandArray, xs))

"Sample from `xs` conditioned on `y`"
function rand(xs::ExecutableRandArray,
              y::RandVar{Bool},
              n::Integer,
              preimage_sampler::Function = point_sample_mc;
              args...)
  preimage_samples = preimage_sampler(y, n; args...)
  [xs(sample) for sample in preimage_samples]
end

function rand(xs::RandArray,
              y::RandVar{Bool},
              n::Integer,
              preimage_sampler::Function = point_sample_mc;
              args...)
  executable_xs = convert(ExecutableRandArray, xs)
  rand(executable_xs, y, n; args...)
end

"Generate a sample from a rand array `xs` conditioned on `y`"
rand(xs::RandArray, y::RandVar{Bool}; args...) =  rand(xs, y, 1; args...)[1]


text:
struct DecisionBoundaryQs <: ModelBasedQs
    occ::OCClassifier
end

function qs_score(qs::DecisionBoundaryQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    return -abs.(predict(qs.occ, x))
end


text:
using GLM, StatsModels, DataFrames, DataFramesMeta, Random, Statistics

df = CSV.read("csv/corei7_long.csv")

train_ratio = 0.1
train_size = round(Int, train_ratio * nrow(df))
train_rows = shuffle(1:nrow(df))[1:train_size]
df_train = df[train_rows, :]

test_rows = [x for x in 1:nrow(df) if !(x in train_rows)]
df_test = df[test_rows, :]

mse(y1, y2) = mean((y1 - y2) .^ 2)

regression_full = lm(@formula(CompileTime ~ ((Expected_insns + Expected_insns ^ 2) +
                                            (Functions + Functions ^ 2) +
                                            (Inlined_percentage + Inlined_percentage ^ 2) +
                                            (Num_partitions + Num_partitions ^ 2)) *
                                            (Parallel + Parallel ^ 2)),
                     df_train)

mse_full = mse(df_test[:, :CompileTime],
               predict(regression_full, select(df_test, Not(:CompileTime))))

max_stars_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_issues_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_forks_repo_name: giulianobelinassi/LTO-Timing-Analysis


text:
module MomentBasedEstimators

using ForwardDiff
using MathProgBase
using Ipopt
using StatsBase
using Reexport
using Distributions
@reexport using Divergences
#using KNITRO
using LinearAlgebra
using LinearAlgebra.BLAS
using Compat
@reexport using CovarianceMatrices
using Calculus
import CovarianceMatrices.RobustVariance

using SparseArrays

const DEFAULT_DIVERGENCE = KullbackLeibler()

## Common
include("md/smoothing.jl")
#include("common/derivatives.jl")
## GMM
include("gmm/iteration_managers.jl")

## MD
include("common/interface.jl")
include("common/api.jl")
include("md/mdp.jl")

include("gmm/mathprogbase.jl")
include("md/mathprogbase.jl")
# include("common/display.jl")
# include("common/post_estimation.jl")

include("common/util.jl")
include("common/post_estimation.jl")

include("common/display.jl")

export GMMEstimator,
       MDEstimator,
       Unconstrained,
       Constrained,
       TwoStepGMM,
       OneStepGMM,
       optimal_W,
       Weighted,
       Unweighted,
       MomentFunction,
       MomentBasedEstimator,
       MinimumDivergenceProblem,
       status,
       estimate!,
       initialize!,
       solver!,
       solve!,
       setparLB!,
       setparUB!,
       setparbounds!,
       setwtsLB!,
       setwtsUB!,
       setwtsbounds!,
       setmfLB!,
       setmfUB!,
       setmfbounds!,
       objval,
       J_test,
       writemime,
       TruncatedSmoother,
       BartlettSmoother,
       MinimumDivergenceProblem
end


text:
function read_table(path::AbstractString)

    df = CSV.read(path, DataFrame)

end

function read_tables(
    files::Vector{<:AbstractString},
    dates::Vector{DateTime};
    timecol::AbstractString = "Date",
    datacol::AbstractString = "H",
)

    dfs = DataFrame[]

    for (i, f) in enumerate(files)
        df = read_table(f)
        df_sel = filter_dates(df, dates; timecol = timecol)
        df_sel[!, datacol] = Float64.(df_sel[!, datacol])
        push!(dfs, df_sel)
    end

    return dfs
end

function read_geoglows(path::AbstractString, rivids::AbstractVector{Int}, dates::Vector{DateTime})
    ds = NCDataset(path,"r")

    qdates = ds["time"][:]

    rividlookup = ds["rivid"][:]

    idx = findall(x -> x in rivids, rividlookup)
    q = ds["Qout"][:][idx,:]

    close(ds)

    dfs = DataFrame[]
    for i in 1:length(rivids)
        df = DataFrame(
            "Date" => qdates,
            "Q" => q[i,:]
        )
        df_sel = filter_dates(df, dates; timecol = "Date")
        push!(dfs,df_sel)
    end

    return dfs
end

function read_stack(
    path::AbstractString,
    datavar::AbstractString;
    timevar::AbstractString = "time",
    xvar::AbstractString = "lon",
    yvar::AbstractString = "lat",
)

    @info "open"
    ds = NCDataset(path, "r")

    @info "read dates"
    dates = ds[timevar][:]
    torder = sortperm(dates)

    dates = dates[torder]
    @info "read lons"
    lons = ds[xvar][:]
    @info "read lats"
    lats = ds[yvar][:]
    @info "read data var"
    data = ds[datavar][:]

    close(ds)

    return data[:, :, torder], dates, lons, lats
end

function read_config(path::AbstractString)
    config = TOML.parsefile(path)

    return config
end

function write_fits(
    outpath,
    center,
    spatial_modes,
    temporal_modes,
    coeffcients,
    correlations,
    varfrac,
    lons,
    lats,
    dates,
)

    refdate = DateTime("1970-01-01")

    ncdates = @. Dates.value(Day(DateTime(dates) - refdate))

    ds = NCDataset(outpath, "c", attrib=Dict(
        "title" => "FIER spatial and temporal Mode coefficients",
        "history" => "date_created: $(now(Dates.UTC))+00:00",
        "institution" => "NASA SERVIR, Univ. AL in Huntsville, Brigham Young Univ., Univ. of Houston",
        "references" => "https://doi.org/10.1016/j.rse.2020.111732"
    ))

    # set up dimensions
    defDim(ds, "lon", length(lons))
    defDim(ds, "lat", length(lats))
    defDim(ds, "time", length(dates))
    defDim(ds, "mode", size(temporal_modes, 2))
    defDim(ds, "order", size(coeffcients, 1))
    defDim(ds, "trial", size(correlations, 1))

    # create variables in dataset
    x = defVar(
        ds,
        "lon",
        Float32,
        ("lon",),
        attrib = Dict(
            "long_name" => "Longitude",
            "units" => "degrees_east",
            "standard_name" => "longitude",
        ),
    )
    y = defVar(
        ds,
        "lat",
        Float32,
        ("lat",),
        attrib = Dict(
            "long_name" => "Latitude",
            "units" => "degrees_north",
            "standard_name" => "latitude",
        ),
    )
    t = defVar(
        ds,
        "time",
        Int32,
        ("time",),
        attrib = Dict(
            "calendar" => "gregorian",
            "units" => "days since 1970-01-01",
        ),
    )
    sm = defVar(
        ds,
        "spatial_modes",
        Float32,
        ("lon", "lat", "mode"),
        fillvalue = -9999.0,
        deflatelevel = 4,
        shuffle=true
    )
    ct = defVar(
        ds,
        "center",
        Float32,
        ("lon", "lat"),
        fillvalue = -9999.0,
        deflatelevel = 4,
        shuffle=true
    )
    tm = defVar(ds, "temporal_modes", Float32, ("time", "mode"))
    cv = defVar(ds, "coefficients", Float32, ("order", "mode"))
    cr = defVar(ds, "correlations", Float32, ("trial", "mode"))
    vf = defVar(ds, "variancefraction", Float32, ("mode",))

    # write data to the variables
    x[:] = lons
    y[:] = lats
    t[:] = ncdates
    sm[:] = spatial_modes
    ct[:] = center
    tm[:] = temporal_modes
    cv[:] = coeffcients
    cr[:] = correlations
    vf[:] = varfrac

    close(ds)

end

function read_fit(path::AbstractString,
    timevar::AbstractString = "time",
    xvar::AbstractString = "lon",
    yvar::AbstractString = "lat",
)

    ds = NCDataset(path, "r")

    @info "read dates"
    dates = ds[timevar][:]
    torder = sortperm(dates)

    dates = dates[torder]
    @info "read lons"
    lons = ds[xvar][:]
    @info "read lats"
    lats = ds[yvar][:]
    @info "read data var"
    spatial_modes = ds["spatial_modes"][:]
    temporal_modes = ds["temporal_modes"][:]
    center = ds["center"][:]
    coefficients = ds["coefficients"][:]
    correlations = ds["correlations"][:]
    variances = ds["variancefraction"][:]

    close(ds)

    return spatial_modes, temporal_modes, center, coefficients, correlations, variances, lons, lats
end

function write_synthesis(
    outpath,
    data,
    dates,
    lons,
    lats;
    water=nothing
)

    refdate = DateTime("1970-01-01")

    ncdates = @. Dates.value(Day(DateTime(dates) - refdate))

    ds = NCDataset(outpath, "c", attrib=Dict(
        "title" => "FIER synthesized outputs",
        "history" => "date_created: $(now(Dates.UTC))+00:00",
        "institution" => "NASA SERVIR, Univ. AL in Huntsville, Brigham Young Univ., Univ. of Houston",
        "references" => "https://doi.org/10.1016/j.rse.2020.111732"
    ))

    # set up dimensions
    defDim(ds, "lon", length(lons))
    defDim(ds, "lat", length(lats))
    defDim(ds, "time", length(dates))

    # create variables in dataset
    x = defVar(
        ds,
        "lon",
        Float32,
        ("lon",),
        attrib = Dict(
            "long_name" => "Longitude",
            "units" => "degrees_east",
            "standard_name" => "longitude",
        ),
    )
    y = defVar(
        ds,
        "lat",
        Float32,
        ("lat",),
        attrib = Dict(
            "long_name" => "Latitude",
            "units" => "degrees_north",
            "standard_name" => "latitude",
        ),
    )
    t = defVar(
        ds,
        "time",
        Int32,
        ("time",),
        attrib = Dict(
            "calendar" => "gregorian",
            "units" => "days since 1970-01-01",
        ),
    )
    d = defVar(
        ds,
        "synthesis",
        Float32,
        ("lon", "lat", "time"),
        fillvalue = -9999.0,
        deflatelevel = 4,
        shuffle=true
    )

    # write data to the variables
    x[:] = lons
    y[:] = lats
    t[:] = ncdates
    d[:] = data

    if !isnothing(water)
        w = defVar(
            ds,
            "water",
            UInt8,
            ("lon", "lat", "time"),
            fillvalue = 255,
            deflatelevel = 4,
            shuffle=true
        )
        w[:] = water
    end

    close(ds)

end

max_stars_repo_name: SERVIR/fier-cli

max_issues_repo_name: SERVIR/fier-cli

max_forks_repo_name: SERVIR/fier-cli