text:
##############################################################################

abstract type Embedding end

##############################################################################
# Embeddings

mutable struct Standardize<:Embedding
    mean
    std
    used::Bool
end

##############################################################################
# generic embedding

mutable struct FunctionEmbedding <: Embedding
    f
end

function embed(e::FunctionEmbedding, u::Array{Float64,1})
    x = u
end

##############################################################################
# standardize

function Standardize()
    return Standardize(0,0,false)
end

#
# embeddings operate on matrices 1 row at a time
# (for some embeddings, like standardize, the results depend
#  on all the entries)
#
function embed(e::Standardize, u::Array{Float64,1})
    if !e.used
        error("cannot standardize single data point")
    end
    x = (u - e.mean[:])./(e.std[:])
end

function unembed(e::Standardize, x::Array{Float64,1})
    if !e.used
        error("cannot standardize single data point")
    end
    return  x.*e.std + e.mean
end

function embed(e::Standardize, U::Array{Float64,2})
    n = size(U,1)
    if !e.used
        e.mean = mean(U,1)
        e.std = std(U,1)
        e.used = true
    end
    return rowwise(u->embed(e,u), U)
end

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

struct IdentityEmbed<:Embedding end
embed(e::IdentityEmbed, u::Array{Float64,1}) = u
unembed(e::IdentityEmbed, u::Array{Float64,1}) = u

##############################################################################
# polyembed

mutable struct PolyEmbed<:Embedding
    degree
end

function embed(e::PolyEmbed, u::Array{Float64,1})
    d = length(u)
    if d != 1
        println("Error: PolyEmbed only operates on scalars")
    end
    x = zeros(e.degree+1)
    for j=0:e.degree
        x[j+1] = u[1]^j
    end
    return x
end

mutable struct PiecewiseEmbed<:Embedding
    knot
end

# this is a neat trick for allowing many knots
# we embed by appending a new feature
# TODO: we need a better API for this
#
# embeddings can be pointwise or not
# can take inputs which are x or u
# can append their output to the input or replace the input
#
function embed(e::PiecewiseEmbed, u::Array{Float64,1})
    n = length(u)
    x = zeros(n+1)
    x[1:n] = u
    x[n+1] = max(u[1]-e.knot,0)
    return x
end

##############################################################################
type AppendOneEmbed<:Embedding end

function embed(e::AppendOneEmbed, x::Array{Float64,1})
    return matrix([1 x'])
end

##############################################################################
# apply embeddings

# allow embeddings to be defined either on the whole dataset
# or per-record
# preferentially apply whole dataset embedding if available
function oneembedding(e, U::Array{Float64,2})
    # if there is a method embedding the whole dataset
    if method_exists(embed, (typeof(e), Array{Float64,2}))
        return embed(e, U)
    end
    # otherwise do one row at a time
    return rowwise(u->embed(e, u), U)
end

oneembedding(e, u::Array{Float64,1}) = embed(e,u)

# apply list of embeddings from right to left
function embed(E::Array, z)
    ne = length(E)
    for i=ne:-1:1
        z = oneembedding(E[i], z)
    end
    return z
end

##############################################################################
# apply unembeddings

function oneunembedding(e, U::Array{Float64,2})
    # if there is a method embedding the whole dataset
    if method_exists(unembed, (typeof(e), Array{Float64,2}))
        return unembed(e, U)
    end
    # otherwise do one row at a time
    return rowwise(u->unembed(e, u), U)
end

oneunembedding(e, u::Array{Float64,1}) = unembed(e,u)

function unembed(E::Array, z)
    ne = length(E)
    for i=1:ne
        z = oneunembedding(E[i], z)
    end
    return z
end

##############################################################################
# interface

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

mutable struct SimpleSource<:DataSource
    U
    V
    Xembed
    Yembed

    function SimpleSource(U, V, Xembed, Yembed)
        if Xembed == false
            Xembed = [AppendOneEmbed(), Standardize()]
        end
        if Yembed == false
        Yembed = [Standardize()]
        end
        return new(matrix(U), matrix(V), Xembed, Yembed)
    end
end

function getXY(S::SimpleSource)
    hasconstfeature = false
    if isa(S.Xembed[1], AppendOneEmbed)
        hasconstfeature = true
    end
    Y = embed(S.Yembed, S.V)
    X = embed(S.Xembed, S.U)
    return X, Y, hasconstfeature
end

getU(S::SimpleSource) = S.U
getV(S::SimpleSource) = S.V

# embed one data record
embedU(S::SimpleSource, u::Array{Float64,1}) = embed(S.Xembed, u)

# unembed one or many targets
unembedY(S::SimpleSource, y) = unembed(S.Yembed, y)


text:
using Flux, KitBase, KitBase.FastGaussQuadrature, ProgressMeter, Plots, LinearAlgebra
import KitML

cd(@__DIR__)
include("../math.jl")

function ComputeSphericalBasisAnalytical(quadpts::AbstractVector{<:AbstractFloat})
    # Hardcoded solution for L = 1, spatialDim = 3
    monomialBasis = zeros(2,length(quadpts))
    for idx_quad in 0:(size(quadpts)[1]-1)
        monomialBasis[1,idx_quad+1]  = 1
        monomialBasis[2,idx_quad+1]  = quadpts[idx_quad+1,1] # x
    end

    return monomialBasis
end

function flux_wall!(
    ff::T1,
    f::T2,
    u::T3,
    dt,
    rot = 1,
) where {
    T1<:AbstractVector{<:AbstractFloat},
    T2<:AbstractVector{<:AbstractFloat},
    T3<:AbstractVector{<:AbstractFloat},
}
    δ = heaviside.(u .* rot)
    fWall = 0.5 .* δ .+ f .* (1.0 .- δ)
    @. ff = u * fWall * dt

    return nothing
end

maxwell_boltzmann(x) = x * log(x) - x

function kinetic_entropy(α, m, weights)
    B = KitBase.maxwell_boltzmann_dual_prime.(α' * m)[:]
    return sum(maxwell_boltzmann.(B) .* weights)
end

begin
    # setup
    set = Setup("radiation", "linesource", "1d1f1v", "kfvs", "bgk", 1, 2, "vanleer", "extra", 0.3, 0.8)

    # physical space
    x0 = 0
    x1 = 1
    nx = 100
    nxg = 0
    ps = PSpace1D(x0, x1, nx, nxg)

    # velocity space
    nu = 28
    points, weights = gausslegendre(nu)
    vs = VSpace1D(points[1], points[end], nu, points, ones(nu) .* (points[end] - points[1]) / (nu - 1), weights)

    # material
    σs = ones(Float32, nx)
    σa = zeros(Float32, nx)
    σt = σs + σa
    σq = zeros(Float32, nx)

    # moments
    L = 1
    ne = 2
    m = ComputeSphericalBasisAnalytical(points)

    # time
    cfl = 0.5
    dt = cfl * ps.dx[1]
    nt = set.maxTime / dt |> floor |> Int
    global t = 0.0

    # NN
    αT = zeros(Float32, nx, ne)
    phiT = zeros(Float32, nx, ne)
    phi_old = deepcopy(phi)
    phi_temp = deepcopy(phi)

    global X = zeros(Float32, ne, 1)
    global Y = zeros(Float32, 1, 1)
    res = KitBase.optimize_closure(α[:, 1], m, weights, phi[:, 1], KitBase.maxwell_boltzmann_dual)
    X[:, 1] .= phi[:, 1]
    Y[1, 1] = kinetic_entropy(α[:, 1], m, weights)
end

begin
    # initial condition
    f0 = 0.0001 * ones(nu)
    phi = zeros(ne, nx)
    for i = 1:nx
        phi[:, i] .= m * f0
    end
    α = zeros(Float32, ne, nx)
    flux = zeros(Float32, ne, nx + 1)
    fη = zeros(nu)
end

anim = @animate for iter = 1:nt
    println("iteration $iter of $nt")

    # mathematical optimizer
    @inbounds for i = 1:nx
        res = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)
        α[:, i] .= res.minimizer
        phi[:, i] .= KitBase.realizable_reconstruct(res.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)

        X = hcat(X, phi[:, i])
        Y = hcat(Y, kinetic_entropy(α[:, i], m, weights))
    end

    flux_wall!(fη, maxwell_boltzmann_dual.(α[:, 1]' * m)[:], points, dt, 1.0)
    for k in axes(flux, 1)
        flux[k, 1] = sum(m[k, :] .* weights .* fη)
    end

    @inbounds for i = 2:nx
        KitBase.flux_kfvs!(fη, KitBase.maxwell_boltzmann_dual.(α[:, i-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i]' * m)[:], points, dt)

        for k in axes(flux, 1)
            flux[k, i] = sum(m[k, :] .* weights .* fη)
        end
    end

    @inbounds for i = 1:nx-1
        for q = 1:1
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (σs[i] * phi[q, i] - σt[i] * phi[q, i]) * dt +
                σq[i] * dt
        end

        for q = 2:ne
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (-σt[i] * phi[q, i]) * dt
        end
    end
    phi[:, nx] .=  phi[:, nx-1]

    global t += dt

    plot(ps.x[1:nx], phi[1, :])
end

gif(anim, "newton1d.gif")

X_old = deepcopy(X)
Y_old = deepcopy(Y)

nn = KitML.ICNNChain(2, 1, [10, 10, 10, 10, 10], tanh)

KitML.sci_train!(nn, (X, Y); epoch = 50, batch = 2, device = gpu)
KitML.sci_train!(nn, (X, Y); epoch = 5, batch = 2)

tmodel = KitML.load_model("best_model.h5"; mode = :tf)
tmodel.fit(permutedims(X), permutedims(Y), epochs = 10)

#gradient(x -> first(nn(x)), phi[:, 1])[1]

X = deepcopy(X_old)
Y = deepcopy(Y_old)

begin
    # initial condition
    f0 = 0.0001 * ones(nu)
    phi = zeros(ne, nx)
    for i = 1:nx
        phi[:, i] .= m * f0
    end
    α = zeros(Float32, ne, nx)
    flux = zeros(Float32, ne, nx + 1)
    fη = zeros(nu)
end

anim = @animate for iter = 1:nt
    println("iteration $iter of $nt")
#for iter = 1:1#nt
    # mathematical optimizer
    #=@inbounds for i = 1:nx
        res = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)
        α[:, i] .= res.minimizer
        phi[:, i] .= KitBase.realizable_reconstruct(res.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)
    end=#

    phi_old .= phi

    # regularization
    @inbounds for i = 1:ne
        phiT[:, i] .= phi[i, :][:]
    end
    αT .= KitML.neural_closure(tmodel, phiT)
    @inbounds for i = 1:ne
        α[i, :] .= αT[:, i]
    end
    @inbounds Threads.@threads for i = 1:nx
        phi_temp[:, i] .= KitBase.realizable_reconstruct(α[:, i], m, weights, KitBase.maxwell_boltzmann_dual_prime)
    end

    counter = 0
    @inbounds for i = 1:nx
        if false#norm(phi_temp[:, i] .- phi_old[:, i], 1) / phi_old[1, i] > 2e-3
            counter +=1

            res = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)
            α[:, i] .= res.minimizer
            phi[:, i] .= KitBase.realizable_reconstruct(res.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)

            X = hcat(X, phi[:, i])
            Y = hcat(Y, kinetic_entropy(α[:, i], m, weights))
        else
            phi[:, i] .= phi_temp[:, i]
        end
    end

    println("newton: $counter of $nx")

    flux_wall!(fη, maxwell_boltzmann_dual.(α[:, 1]' * m)[:], points, dt, 1.0)
    for k in axes(flux, 1)
        flux[k, 1] = sum(m[k, :] .* weights .* fη)
    end

    @inbounds for i = 2:nx
        KitBase.flux_kfvs!(fη, KitBase.maxwell_boltzmann_dual.(α[:, i-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i]' * m)[:], points, dt)

        for k in axes(flux, 1)
            flux[k, i] = sum(m[k, :] .* weights .* fη)
        end
    end

    @inbounds for i = 1:nx-1
        for q = 1:1
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (σs[i] * phi[q, i] - σt[i] * phi[q, i]) * dt +
                σq[i] * dt
        end

        for q = 2:ne
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (-σt[i] * phi[q, i]) * dt
        end
    end
    phi[:, nx] .=  phi[:, nx-1]

    global t += dt

    if iter%19 == 0
        tmodel.fit(permutedims(X), permutedims(Y), epochs=1)

        #=X = zeros(Float32, 1, ne)
        Y = zeros(Float32, 1, 1)
        res = KitBase.optimize_closure(α[:, nx], m, weights, phi[:, nx], KitBase.maxwell_boltzmann_dual)
        X[1, :] .= phi[:, nx]
        Y[1, 1] = kinetic_entropy(α[:, nx], m, weights)=#
    end

    plot(ps.x[1:nx], phi[1, :])
end

gif(anim, "unified_1d.gif")
gif(anim, "neural_1d.gif")
plot(ps.x[1:nx], phi[1, :])

using BenchmarkTools, Optim

@btime KitML.neural_closure(tmodel, phiT[:, :])

@btime for i = 1:nx
    #res = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)

    loss(_α) = sum(maxwell_boltzmann_dual.(_α' * m)[:] .* weights) - dot(_α, phi[:, i])
    res = Optim.optimize(loss, α[:, i], Optim.Newton())
end

ext: jl

max_stars_repo_path: example/1D propagation/jmodel.jl

max_issues_repo_path: example/1D propagation/jmodel.jl

max_forks_repo_path: example/1D propagation/jmodel.jl


text:
# # Mass dashpot spring

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/massDashpotSpring.ipynb)
#

# ## Free oscillations

using StructuralDynamicsODESolvers, Plots, LinearAlgebra

k  = 2 ; m  = .5 ;  c = 0 ;
u0 = 1 ; v0 = 0 ;

M = m * ones(1, 1)
C = c * ones(1, 1)
K = k * ones(1, 1)
R = zeros(1)

sys = SecondOrderAffineContinuousSystem(M, C, K, R)

U₀ = u0 * ones(1); V₀ = v0 * ones(1);

ivp_free = InitialValueProblem(sys, (U₀, V₀))

NSTEPS = 1000 ;
Δt = 0.005 ;

#-

alg = Bathe(Δt = Δt )
sol = solve(ivp_free, alg, NSTEPS=NSTEPS);

#-

# The following command is the same as `plot(times(sol), displacements(sol, 1))`.

plot(sol, vars=(0, 1))

# ## Forced oscillations
#
# ### Problem definition
# Let us consider now a forcing term $f(t) = A_f \sin(ω_f \cdot t)$
#
ωN = sqrt(k/m)
ωf = ωN * 2
Af = 10.0
R  = [ [ Af * sin(ωf * Δt * (i-1) ) ] for i in 1:NSTEPS+1];

# ### Second order problem resolution

X   = nothing # state constraints are ignored
B   = ones(1, 1)
sys = SecondOrderConstrainedLinearControlContinuousSystem(M, C, K, B, X, R)

ivp_forced_secOrder = InitialValueProblem(sys, (U₀, V₀))

alg = Bathe(Δt = Δt )
sol_secOrder = solve(ivp_forced_secOrder, alg, NSTEPS=NSTEPS);

# ### First order homogeneization formulation
# The problem can be re-formulated as a first order and homogeneous one given by
#
#
# ```math
# \left\{
# \begin{array}{l}
# \dot{u} = v \\
# \dot{v} = -\omega_N^2 u + u_f/m \\
# \dot{u_f} = v_f \\
# \dot{v_f} = -\omega_f^2 u_f
# \end{array}
# \right.
# ```
# The new vector of variables is
#
# ```math
# \textbf{x} = [ u, v, u_f, v_f ]^T
# ```

K = [     0 1     0 0 ;
      -ωN^2 0     1/m 0 ;
          0 0     0 1 ;
          0 0 -ωf^2 0 ] ;

C = -Diagonal(ones(4))
M = zeros(4,4)
R = zeros(4)

sys = SecondOrderAffineContinuousSystem(M, C, K, R)

U₀ = [u0; v0; 0; ωf*Af ] ;

ivp_forced_firOrder = InitialValueProblem(sys, (U₀, U₀) )

alg = BackwardEuler(Δt = Δt )
sol_firOrderA = solve(ivp_forced_firOrder, alg, NSTEPS=NSTEPS);

NSTEPS = NSTEPS*3 ; alg = BackwardEuler(Δt = Δt/3.0 )
sol_firOrderB = solve(ivp_forced_firOrder, alg, NSTEPS=NSTEPS);

# The solution obtained is

plot(sol_secOrder, vars=(0, 1), xlab="time" )
plot!(sol_firOrderA, vars=(0, 1))
plot!(sol_firOrderB, vars=(0, 1))


text:
function get_logit_ionosphere(;λ=0.001)
    checkdownload_ionosphere()

    rawdata = readdlm(joinpath(instances_dir, "ionosphere.data"), ',')
    A = Matrix{Float64}(rawdata[:, 1:end-1])
    n = size(A, 2)

    y = map(x-> x=="g" ? 1.0 : -1.0, rawdata[:, end])

    return LogisticPb(A, y, regularizer_l1(λ), n, zeros(1), l1Manifold(zeros(1)), 0.0)
end

function get_logit_gisette(λ = 6.67e-4)
    sampledim = 5000

    if isfile(joinpath(instances_dir, "gisette_scale.jld"))
        d = load(joinpath(instances_dir, "gisette_scale.jld"))
        A = d["A"]
        y = d["y"]
    else
        println("Reading ", joinpath(instances_dir, "gisette_scale"))
        rawdata = readdlm(joinpath(instances_dir, "gisette_scale"), ' ')

        println("Assembling problem matrix, vector.")
        nsamples = size(rawdata, 1)

        A = zeros(Float64, nsamples, sampledim)
        for i in 1:size(rawdata, 1), j in 2:size(rawdata, 2)
            s = rawdata[i, j]
            s == "" && continue

            indcol = parse.(Int, split(s, ":")[1])
            val = parse.(Float64, split(s, ":")[2])
            A[i, indcol] = val
        end

        y = Vector{Float64}(rawdata[:, 1])

        println("")
        save(joinpath(instances_dir, "gisette_scale.jld"), "A", A, "y", y)
    end

    return LogisticPb(A, y, regularizer_l1(λ), sampledim, zeros(1), l1Manifold(zeros(1)), 0.0)
end

function get_random_logit(;n=80, m=85, sparsity=0.5, λ=0.1, seed=1234)

    Random.seed!(seed)
    A = rand(m, n)*10
    Random.seed!(seed+3)
    y = Vector([rand() > sparsity ? 1.0 : -1.0 for i in 1:m])

    return LogisticPb(A, y, regularizer_l1(λ), n, zeros(1), l1Manifold(zeros(1)), 0.0)
end

function get_logit_MLE(;n=20, m=15, sparsity=0.5, seed=1234, λ=0.01, diff=0.5)
    @assert 0 ≤ sparsity ≤ 1

    Random.seed!(seed)
    A = rand(Normal(), m, n)

    Random.seed!(seed+1)
    x0 = rand(Normal(), n)
    Random.seed!(seed+2)
    optstructure = l1Manifold(rand(Bernoulli(1-sparsity), n))
    x0 = project(optstructure, x0)

    y = zeros(m)
    for i in 1:m
        Random.seed!(seed+i)
        if rand(Bernoulli(diff + (1-diff)*σ(dot(A[i, :], x0))))
            y[i] = 1.0
        else
            y[i] = -1.0
        end
    end

    return LogisticPb(A, y, regularizer_l1(λ), n, x0, optstructure, 0.0)
end

function get_logit_MLE_elastic_net(;n=20, m=15, sparsity=0.5, seed=1234, λ=0.01, lambda2=0.02, diff=0.5)
    @assert 0 ≤ sparsity ≤ 1

    Random.seed!(seed)
    A = rand(Normal(), m, n)

    Random.seed!(seed+1)
    x0 = rand(Normal(), n)
    Random.seed!(seed+2)
    optstructure = l1Manifold(rand(Bernoulli(1-sparsity), n))
    x0 = project(optstructure, x0)

    y = zeros(m)
    for i in 1:m
        Random.seed!(seed+i)
        if rand(Bernoulli(diff + (1-diff)*σ(dot(A[i, :], x0))))
            y[i] = 1.0
        else
            y[i] = -1.0
        end
    end

    return LogisticPb(A, y, regularizer_l1(λ), n, x0, optstructure, lambda2)
end


text:
@testset "Cartesian" begin
  x = [π, exp(1), √2]
  for ct ∈ (ToCartesian(x->x, 3), FromCartesian(x->x, 3))
    @test J(ct, x) == 1
    m = zeros(3,3) + I
    gᵢⱼ(ct, x) ≈ m
    gⁱʲ(ct, x) ≈ m
    dl(ct, x, Covariant) ≈ ones(3)
    dl(ct, x, Contravariant) ≈ ones(3)

    v = CovariantVector(x->[x[1], x[2], x[3]], ct, true)
    @test div(v)(x) ≈ 3
    @test curl(v)(x) ≈ zeros(3)
    v = CovariantVector(x->[x[3], x[1], x[2]], ct, true)
    @test div(v)(x) ≈ 0
    @test curl(v)(x) ≈ ones(3)
    @test div(curl(v))(x) ≈ 0
    @test div(curl(v))(x) ≈ 0
    v = CovariantVector(x->[-x[2], -x[3], -x[1]], ct, true)
    @test div(v)(x) ≈ 0
    @test curl(v)(x) ≈ ones(3)
    @test div(curl(v))(x) ≈ 0
    @test div(curl(v))(x) ≈ 0
  end
end

max_stars_repo_name: jwscook/DifferentialGeometries

max_issues_repo_name: jwscook/DifferentialGeometries

max_forks_repo_name: jwscook/DifferentialGeometries


text:
# Biryani -- A conversion and validation toolbox
# By: Emmanuel Raviart <emmanuel@raviart.com>
#
# Copyright (C) 2015 Emmanuel Raviart
# https://github.com/eraviart/Biryani.jl
#
# This file is part of Biryani.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

importall Biryani.DatesConverters

# date_to_iso8601_string
@test Convertible(Date(2012, 3, 4)) |> date_to_iso8601_string |> to_value == "2012-03-04"
@test Convertible(nothing) |> date_to_iso8601_string |> to_value === nothing
@test Convertible(Date(2012, 3, 4)) |> fail("Initial error.") |> date_to_iso8601_string |> to_value_error == (
  Date(2012, 3, 4), "Initial error.")

# iso8601_input_to_date
@test Convertible("2012-03-04") |> iso8601_input_to_date |> to_value == Date(2012, 3, 4)
@test Convertible("2012-03") |> iso8601_input_to_date |> to_value == Date(2012, 3, 1)
@test Convertible("2012") |> iso8601_input_to_date |> to_value == Date(2012, 1, 1)
@test Convertible("1500-01-01") |> iso8601_input_to_date |> to_value == Date(1500, 1, 1)
@test Convertible("2099-12-31") |> iso8601_input_to_date |> to_value == Date(2099, 12, 31)
@test Convertible("   2012-03-04\n  ") |> iso8601_input_to_date |> to_value == Date(2012, 3, 4)
@test Convertible("   2012-03\n  ") |> iso8601_input_to_date |> to_value == Date(2012, 3, 1)
@test Convertible("   2012\n  ") |> iso8601_input_to_date |> to_value == Date(2012, 1, 1)
@test Convertible("1499-12-31") |> iso8601_input_to_date |> to_value_error == ("1499-12-31",
  "Invalid ISO-8601 format for date.")
@test Convertible("2100-01-01") |> iso8601_input_to_date |> to_value_error == ("2100-01-01",
  "Invalid ISO-8601 format for date.")
@test Convertible("2012-03-04 05:06:07") |> iso8601_input_to_date |> to_value_error == ("2012-03-04 05:06:07",
  "Invalid ISO-8601 format for date.")
@test Convertible("2012-03-04T05:06:07") |> iso8601_input_to_date |> to_value_error == ("2012-03-04T05:06:07",
  "Invalid ISO-8601 format for date.")
@test Convertible("today") |> iso8601_input_to_date |> to_value_error == ("today", "Invalid ISO-8601 format for date.")
@test Convertible(nothing) |> iso8601_input_to_date |> to_value === nothing
@test Convertible("   \n  ") |> iso8601_input_to_date |> to_value === nothing
@test Convertible("2012-03-04") |> fail("Initial error.") |> iso8601_input_to_date |> to_value_error == ("2012-03-04",
  "Initial error.")

# iso8601_string_to_date
@test Convertible("2012-03-04") |> iso8601_string_to_date |> to_value == Date(2012, 3, 4)
@test Convertible("2012-03") |> iso8601_string_to_date |> to_value == Date(2012, 3, 1)
@test Convertible("2012") |> iso8601_string_to_date |> to_value == Date(2012, 1, 1)
@test Convertible("2012-03-04 05:06:07") |> iso8601_string_to_date |> to_value_error == ("2012-03-04 05:06:07",
  "Invalid ISO-8601 format for date.")
@test Convertible("2012-03-04T05:06:07") |> iso8601_string_to_date |> to_value_error == ("2012-03-04T05:06:07",
  "Invalid ISO-8601 format for date.")
@test Convertible("today") |> iso8601_string_to_date |> to_value_error == ("today", "Invalid ISO-8601 format for date.")
@test Convertible("") |> iso8601_string_to_date |> to_value_error == ("", "Invalid ISO-8601 format for date.")
@test Convertible(nothing) |> iso8601_string_to_date |> to_value === nothing
@test Convertible("2012-03-04") |> fail("Initial error.") |> iso8601_string_to_date |> to_value_error == ("2012-03-04",
  "Initial error.")

# to_date
@test Convertible("2012-03-04") |> to_date |> to_value == Date(2012, 3, 4)
@test Convertible("2012-03") |> to_date |> to_value == Date(2012, 3, 1)
@test Convertible("2012") |> to_date |> to_value == Date(2012, 1, 1)
@test Convertible(2012) |> to_date |> to_value == Date(2012, 1, 1)
@test Convertible((2012, 3, 4)) |> to_date |> to_value == Date(2012, 3, 4)
@test Convertible((2012, 3)) |> to_date |> to_value == Date(2012, 3, 1)
@test Convertible((2012,)) |> to_date |> to_value == Date(2012, 1, 1)
@test Convertible(nothing) |> to_date |> to_value === nothing
@test Convertible("2012-03-04") |> fail("Initial error.") |> to_date |> to_value_error == ("2012-03-04",
  "Initial error.")


text:
using AGFFileReader
using Documenter

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

#TODO add plotting examples

makedocs(;
    modules=[AGFFileReader],
    authors="rambunctiousapple",
    repo="https://github.com/rambunctiousapple/AGFFileReader.jl/blob/{commit}{path}#{line}",
    sitename="AGFFileReader.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rambunctiousapple.github.io/AGFFileReader.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Reference" => "reference.md"
    ],
)

deploydocs(;
    repo="github.com/rambunctiousapple/AGFFileReader.jl.git",
    devbranch="main",
    push_preview = true,
)

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
function plotqq!(
    ax::Axis,
    P::AbstractVector;
    xstep::Real = 1,
    ystep::Real = 2
)

    n = length(P)
    high = Array{Float64}(undef, n)
    low = Array{Float64}(undef, n)
    for i in 1:n
        high[i] = -log(10, quantile(Beta(i, n + 1 - i), 0.025))
        low[i] = -log(10, quantile(Beta(i, n + 1 - i), 0.975))
    end
    obs = -log.(10, sort(P))
    obs = clamp.(obs, 0.0, -log(10, floatmin(0.0)))
    expect = -log.(10, (collect(1:n) .- 0.5) ./ n)
    xmax = ceil(maximum(expect))
    ymax = ceil(maximum(obs))
    xticks = 0:xstep:xmax
    yticks = (ymax % 2 == 0) ? (0:ystep:ymax) : (0:ystep:(ymax + 1))
    λgc = cquantile(Chisq(1), median(P)) / quantile(Chisq(1), 0.5)

    band!(ax, expect, low, high, color = :gray90)
    scatter!(ax, expect, obs, markersize = 2, color = :black)
    lines!(ax, expect, expect, color = :red2, linewidth = 0.75)
    text!(ax, "λgc = $(round(λgc, digits = 2))", position = (xmax - 0.5, 1),
        textsize = 6, align = (:right, :bottom))
    xlims!(ax, 0, xmax)
    ylims!(ax, 0, ymax)
    hidedecorations!(ax, label = false, ticklabels = false, ticks = false)
    ax.aspect = AxisAspect(1)
    ax.xlabel = "Expected -log[p]"
    ax.ylabel = "Observed -log[p]"
    ax.xticklabelsize = 6
    ax.yticklabelsize = 6
    ax.xlabelsize = 8
    ax.ylabelsize = 8
    ax.xticksize = 3
    ax.yticksize = 3
    ax.xtickwidth = 0.75
    ax.ytickwidth = 0.75
    ax.spinewidth = 0.75
    ax.xticks = xticks
    ax.yticks = yticks
end

"""
    plotqq!(ax::Axis, gwas::DataFrame; kwargs)
    plotqq!(ax::Axis, P::AbstractVector; kwargs)

Plot QQ plot of `P` values where the expected distribution is the uniform distribution.

Keyword arguments include `xstep::Real` and `ystep::Real` for x and y axes ticks step sizes.
"""
plotqq!(ax::Axis, gwas::DataFrame; kwargs...) = plotqq!(ax, gwas.P; kwargs...)


text:
_to_vec(x::AbstractArray) = x
_to_vec(x::Tuple) = collect(x)
_to_vec(x::Nothing) = x
_to_vec(x) = [x]

_get_cols(table) = propertynames(Tables.columns(table))

_to_table(x, ::Nothing) = Tables.table(x)  # so that we use the defaults in Tables.table
_to_table(x, header) = Tables.table(x, header=header)


text:
"""
    Pads all the dimensions, in constrast with @PadData
    that pads only x y dims.
"""
function PadDataCpu(pd,data)
    x = size(data,1)
    y = size(data,2)
    z = size(data,3)
    padx = 1:(x+2pd)
    pady = 1:(y+2pd)
    padz = 1:(z+2pd)
    datax = (pd+1):(x+pd)
    datay = (pd+1):(y+pd)
    dataz = (pd+1):(z+pd)
    return PaddedView(0, data, (padx,pady,padz), (datax, datay, dataz))
end

function cpu_star_stencil(org_data, radius, coefs, t_steps)
    data = PadDataCpu(radius, org_data)

    function stencil(data, radius, coefs)
        out = zeros(Float32,size(data))
        for i = (radius+1):(size(data,1)-radius)
            for j = (radius+1):(size(data,2)-radius)
                for k = (radius+1):(size(data,3)-radius)
                    c = coefs[1]*data[i,j,k]
                    for r = 1:radius
                        c += coefs[r+1]*(data[i + r,j,k] + data[i - r,j,k] +
                        data[i,j + r,k] + data[i,j - r,k] +
                        data[i,j,k + r] + data[i,j,k - r])
                    end
                    out[i,j,k] = c
                end
            end
        end
        return out
    end

    cout = stencil(data, radius, coefs)
    for t = 2:t_steps
        cout = stencil(cout, radius, coefs)
    end

    return @view cout[:,:,(radius+1):(size(org_data,3)+radius)]
end

function cpu_star_stencil_v(org_data, vsq, radius, coefs, t_steps)
        data = PadDataCpu(radius, org_data)
        v = PadDataCpu(radius, vsq)

            function stencil(data, vsq, radius, coefs)
                out = zeros(Float32,size(data))
                for i = (radius+1):(size(data,1)-radius)
                    for j = (radius+1):(size(data,2)-radius)
                        for k = (radius+1):(size(data,3)-radius)
                            c = coefs[1]*data[i,j,k]
                            for r = 1:radius
                                c += coefs[r+1]*(data[i + r,j,k] + data[i - r,j,k] +
                                data[i,j + r,k] + data[i,j - r,k] +
                                data[i,j,k + r] + data[i,j,k - r])
                            end
                            out[i,j,k] = c*vsq[i,j,k]
                        end
                    end
                end
                return out
            end

            cout = stencil(data,v, radius, coefs)
            for t = 2:t_steps
                cout = stencil(cout,v, radius, coefs)
            end
            return @view cout[:,:,(radius+1):(size(org_data,3)+radius)]
end

function cpu_star_stencil_v_nv(org_data, vsq, radius, coefs, t_steps)
        data = PadDataCpu(radius, org_data)
        v = PadDataCpu(radius, vsq)

            function stencil(data, vsq, radius, coefs)
                out = zeros(Float32,size(data))
                for i = (radius+1):(size(data,1)-radius)
                    for j = (radius+1):(size(data,2)-radius)
                        for k = (radius+1):(size(data,3)-radius)
                            c = coefs[1]*data[i,j,k]
                            for r = 1:radius
                                c += coefs[r+1]*(data[i + r,j,k] + data[i - r,j,k] +
                                data[i,j + r,k] + data[i,j - r,k] +
                                data[i,j,k + r] + data[i,j,k - r])
                            end
                            out[i,j,k] = c*vsq[i,j,k] + c
                        end
                    end
                end
                return out
            end

            cout = stencil(data,v, radius, coefs)
            for t = 2:t_steps
                cout = stencil(cout,v, radius, coefs)
            end
            return @view cout[:,:,(radius+1):(size(org_data,3)+radius)]
end

function cpu_star_stencil_vsq_prev_val(org_data,vsq,
            radius, coefs, t_steps, prev_time)
    data = PadDataCpu(radius, org_data)
    data = Array(data)
    out = zeros(Float32,size(data))
    v = PadDataCpu(radius, vsq)

    function stencil(data, out, vsq, radius, coefs, pc)
        for i = (radius+1):(size(data,1)-radius)
            for j = (radius+1):(size(data,2)-radius)
                for k = (radius+1):(size(data,3)-radius)
                    #out[i,j,k] != 0 && println(pc*out[i,j,k])
                    c = coefs[1]*data[i,j,k]
                    for r = 1:radius
                        c += coefs[r+1]*(data[i + r,j,k] + data[i - r,j,k] +
                        data[i,j + r,k] + data[i,j - r,k] +
                        data[i,j,k + r] + data[i,j,k - r])
                    end
                    out[i,j,k] =c*vsq[i,j,k] + c + pc*out[i,j,k]
                end
            end
        end
        return out
    end

    out = stencil(data, out, v, radius, coefs, prev_time)
    for t = 2:t_steps
        data,out = out,data
        out = stencil(data, out, v, radius, coefs, prev_time)
    end

    #if t_steps%2 == 1
    return @view out[:,:,(radius+1):(size(org_data,3)+radius)]
    #else
    #    return @view data[:,:,(radius+1):(size(org_data,3)+radius)]
    #end

end

function cpu_dense_stencil(org_data, radius, coefs,t_steps)
    data = PadDataCpu(radius, org_data)
    data = Array(data)
    dim = size(data)
    out = zeros(Float32, dim)
    function stencil(data, out, radius, coefs)
        for x = (radius+1):(dim[1]-radius)
            for y = (radius+1):(dim[2]-radius)
                for z = (radius+1):(dim[3]-radius)
                    res = coefs[1]*data[x,y,z]
                    for ix = -radius:radius
                        for iy = -radius:radius
                            for iz = -radius:radius
                                if ix == 0 &&  iy == 0 &&  iz == 0
                                    continue
                                end
                                ci = max(abs.([ix, iy, iz])...)
                                ci = coefs[ci+1]
                                res += ci*data[x+ix,y+iy,z+iz]
                            end
                        end
                    end
                    out[x,y,z] = res
                end
            end
        end
        return out
    end
    out = stencil(data, out, radius, coefs)
    for t = 2:t_steps
        data,out = out,data
        out = stencil(data, out, radius, coefs)
    end

    #if t_steps%2 == 1
    return @view out[:,:,(radius+1):(size(org_data,3)+radius)]
    return @view cout[:,:,(radius+1):(size(org_data,3)+radius)]
end


text:
module Assimp

include("LibAssimp.jl")
using .LibAssimp

end # module


text:
import Base.setprecision
################################################################################
#
#  Show function
#
################################################################################

function AbstractAlgebra.expressify(a::LocalFieldElem; context = nothing)
  return AbstractAlgebra.expressify(a.data, var(parent(a)), context = context)
end

function Base.show(io::IO, a::LocalFieldElem)
  print(io, AbstractAlgebra.obj_to_string(a, context = io))
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy(x::LocalFieldElem{S, T}) where {S, T}
  return LocalFieldElem{S, T}(parent(x), deepcopy(x.data), precision(x))
end

################################################################################
#
#  Precision
#
################################################################################

function _generator_valuation(K::LocalField{S, T}) where {S <: Union{padic, qadic}, T <: LocalFieldParameter}
  f = defining_polynomial(K)
  return fmpq(valuation(coeff(f, 0)), degree(K))
end

function _generator_valuation(K::LocalField)
  f = defining_polynomial(K)
  return divexact(valuation(coeff(f, 0)), degree(K))
end

function compute_precision(K::LocalField, a::Generic.Poly)
  v1 = _generator_valuation(K)
  if iszero(v1)
    v = fmpq()
  else
    v = inv(v1)
  end
  prec = precision(coeff(a, 0))*ramification_index(K)
  for i = 1:degree(a)
    c = coeff(a, i)
    prec = min(prec, precision(c)+Int(numerator(ceil(i*v))))
  end
  return prec
end

Nemo.precision(a::LocalFieldElem) = a.precision

function Nemo.precision(a::Generic.Poly{S}) where S <: LocalFieldElem
  if iszero(a)
    return precision(base_ring(a))
  end
  return minimum(precision, coefficients(a))
end

function Base.setprecision(a::LocalFieldElem, n::Int)
  b = deepcopy(a)
  return setprecision!(b, n)
end

function setprecision!(a::padic, n::Int)
  return setprecision(a, n)
end

function setprecision!(a::LocalFieldElem, n::Int)
  K = parent(a)
  e = ramification_index(K)
  d, r = divrem(n, e)
  if !iszero(r)
    d += 1
  end
  for i = 0:degree(a.data)
    setcoeff!(a.data, i, setprecision!(coeff(a, i), d))
  end
  a.precision = n
  return a
end

################################################################################
#
#  Parent
#
################################################################################

parent(a::LocalFieldElem) = a.parent

parent_type(a::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter} = LocalField{S, T}
parent_type(::Type{LocalFieldElem{S, T}}) where {S <: FieldElem, T <: LocalFieldParameter} = LocalField{S, T}

################################################################################
#
#   Coefficients
#
################################################################################

coeff(a::LocalFieldElem, i::Int) = coeff(a.data, i)
setcoeff!(a::LocalFieldElem, i::Int) = setcoeff!(a.data, i)

################################################################################
#
#  Zero/One
#
################################################################################

iszero(a::LocalFieldElem) = iszero(a.data)
isone(a::LocalFieldElem) = isone(a.data)
isunit(a::LocalFieldElem) = !iszero(a)

function O(K::LocalField, prec::T) where T <: Union{Integer, fmpz}
  d, r = divrem(prec, ramification_index(K))
  if !iszero(r)
    r += 1
  end
  return K(O(base_field(K), d))
end

function zero(K::LocalField)
  a = zero(parent(defining_polynomial(K)))
  return setprecision(K(a), precision(K))
end

(K::LocalField)() = zero(K)

function one(K::LocalField)
  a = one(parent(defining_polynomial(K)))
  return setprecision(K(a), precision(K))
end

function zero!(a::LocalFieldElem)
  K = parent(a)
  zero!(a.data)
  a.data = setprecision(a.data, precision(K))
  return a
end

function Base.:(==)(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S, T}
  for i = 0:max(degree(a.data), degree(b.data))
    if coeff(a, i) != coeff(b, i)
      return false
    end
  end
  return true
end

################################################################################
#
#  Coercion
#
################################################################################

function (K::LocalField{S, T})(a::Integer) where {S <: FieldElem, T <: LocalFieldParameter}
  el =  K(parent(defining_polynomial(K))(a))
  return setprecision!(el, precision(K))
end

function (K::LocalField{S, T})(a::Union{fmpz, fmpq}) where {S <: FieldElem, T <: LocalFieldParameter}
  el =  K(parent(defining_polynomial(K))(a))
  return setprecision!(el, precision(K))
end

function (K::LocalField{S, T})(a::U) where {U <: Union{padic, qadic}, S <: FieldElem, T <: LocalFieldParameter}
  return K(parent(defining_polynomial(K))(a))
end

function (K::LocalField{S, T})(a::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  return a
end

function (K::LocalField{S, T})(a::LocalFieldElem{U, V}) where {S <: FieldElem, U <: FieldElem, T <: LocalFieldParameter, V <: LocalFieldParameter}
  if parent(a) === K
    return a
  elseif base_field(K) === parent(a)
    kt = parent(defining_polynomial(K))
    return K(kt(a.data))
  else
    return K(base_field(K)(a))
  end
end

function (K::LocalField{S, T})(p::Generic.Poly{S}) where {S <: FieldElem, T <: LocalFieldParameter}
  if degree(p) >= degree(K)
    p = mod(p, defining_polynomial(K))
  end
  return LocalFieldElem{S, T}(K, p, compute_precision(K, p))
end

function (Rx::Generic.PolyRing{S})(a::LocalFieldElem{S, T}) where {S, T}
  @assert base_ring(Rx) == base_field(parent(a))
  return a.data
end

################################################################################
#
#  Valuation
#
################################################################################

function valuation(a::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  return valuation(norm(a))//degree(parent(a))
end

function valuation(a::LocalFieldElem{S, EisensteinLocalField}) where S <: FieldElem
  if iszero(a)
    error("Input element is zero")
  end
  K = parent(a)
  e = absolute_ramification_index(K)
  i = 0
  c = coeff(a, i)
  while iszero(c)
    i += 1
    c = coeff(a, i)
  end
  vc = valuation(c)
  v = vc + fmpq(i, e)
  for j = i+1:degree(K)
    c = coeff(a, j)
    if iszero(c)
      continue
    end
    vc = valuation(c)
    vnew = vc + fmpq(j, e)
    if vnew < v
      v = vnew
    end
  end
  return v
end

function valuation(a::LocalFieldElem{S, UnramifiedLocalField}) where S <: FieldElem
  if iszero(a)
    error("Input element is zero")
  end
  K = parent(a)
  e = absolute_ramification_index(K)
  i = 0
  c = coeff(a, i)
  while iszero(c)
    i += 1
    c = coeff(a, i)
  end
  v = valuation(c)
  for j = i+1:degree(K)
    c = coeff(a, j)
    vc = valuation(c)
    if vc < v
      v = vc
    end
  end
  return v
end

function check_parent(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  #=
  if parent(a) !== parent(b)
    @show parent(a)
    @show parent(b)
    error("Wrong parents!")
  end
  =#
  return nothing
end

################################################################################
#
#  Representation matrix
#
################################################################################

function elem_to_mat_row!(M::MatElem, i::Int, a::LocalFieldElem)
  K = parent(a)
  @assert ncols(M) == degree(K)
  for j = 1:degree(K)
    M[i, j] = coeff(a, j-1)
  end
  return nothing
end

function representation_matrix(a::LocalFieldElem)
  K = parent(a)
  g = gen(K)
  el = a
  M = zero_matrix(base_field(K), degree(K), degree(K))
  elem_to_mat_row!(M, 1, a)
  for i = 1:degree(K)-1
    el *= g
    elem_to_mat_row!(M, i+1, el)
  end
  return M
end

################################################################################
#
#  Norm
#
################################################################################

function norm(a::LocalFieldElem)
  K = parent(a)
  return resultant(a.data, defining_polynomial(K))
end

function absolute_norm(a::LocalFieldElem)
  return absolute_norm(norm(a))
end

function absolute_norm(a::padic)
  return a
end

function absolute_norm(a::qadic)
  return norm(a)
end

function norm(a::LocalFieldElem, F::LocalField)
  K = parent(a)
  if K === F
    return norm(a)::elem_type(F)
  else
    return norm(norm(a), F)::elem_type(F)
  end
end

function norm(a::LocalFieldElem, F::FlintPadicField)
  return absolute_norm(a)
end

################################################################################
#
#  Trace
#
################################################################################

function tr(a::LocalFieldElem)
  K = parent(a)
  assure_traces(K)
  res = base_field(K)()
  for i = 0:degree(K)-1
    res += K.traces_basis[i+1]*coeff(a, i)
  end
  return res
end

function absolute_tr(a::LocalFieldElem)
  return absolute_tr(tr(a))
end

function absolute_tr(a::padic)
  return a
end

function absolute_tr(a::qadic)
  return tr(a)
end

function tr(a::LocalFieldElem, F::LocalField)
  K = parent(a)
  if K === F
    return tr(a)::elem_type(F)
  else
    return tr(tr(a), F)::elem_type(F)
  end
end

function tr(a::LocalFieldElem, F::FlintPadicField)
  return absolute_tr(a)
end

################################################################################
#
#  Minpoly
#
################################################################################

function minpoly(a::T, Kx::PolyRing = PolynomialRing(parent(a), "t", cached = false)[1]) where T <: Union{LocalFieldElem, qadic}
  return squarefree_part(norm(gen(Kx)-a))
end

function absolute_minpoly(a::LocalFieldElem)
  f = minpoly(a)
  return _absolute_minpoly(f)
end

function _absolute_minpoly(p::Generic.Poly{T}) where T <: Union{qadic, LocalFieldElem}
  return _absolute_minpoly(squarefree_part(norm(p)))
end

function _absolute_minpoly(p::Generic.Poly{padic})
  return squarefree_part(p)
end

function absolute_minpoly(a::padic, parent = PolynomialRing(parent(a), "x"))
  return gen(parent)-a
end

function absolute_minpoly(a::qadic)
  return minpoly(a)
end

################################################################################
#
#  Basic operations
#
################################################################################

function Base.:-(a::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  return LocalFieldElem{S, T}(parent(a), -a.data, precision(a))
end

function Base.:+(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  return LocalFieldElem{S, T}(parent(a), a.data+b.data, min(precision(a), precision(b)))
end

function Base.:-(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  return LocalFieldElem{S, T}(parent(a), a.data-b.data, min(precision(a), precision(b)))
end

function Base.:*(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  pol = mod(a.data*b.data, defining_polynomial(parent(a)))
  res =  LocalFieldElem{S, T}(parent(a), pol, compute_precision(parent(a), pol))
  return res
end

function Base.:(//)(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  ib = inv(b)
  pol = mod(a.data*ib.data, defining_polynomial(parent(a)))
  res =  LocalFieldElem{S, T}(parent(a), pol, compute_precision(parent(a), pol))
  return res
end

function divexact(a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  return a//b
end

function add!(c::LocalFieldElem{S, T}, a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  c.parent = a.parent
  c.data = a.data + b.data
  c.precision = min(a.precision, b.precision)
  return c
end

function addeq!(c::LocalFieldElem{S, T}, a::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, c)
  c.data = add!(c.data, c.data, a.data)
  c.precision = min(a.precision, c.precision)
  return c
end

function sub!(c::LocalFieldElem{S, T}, a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  c.parent = a.parent
  c.data = a.data - b.data
  c.precision = min(a.precision, b.precision)
  return c
end

function mul!(c::LocalFieldElem{S, T}, a::LocalFieldElem{S, T}, b::LocalFieldElem{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  check_parent(a, b)
  c.parent = a.parent
  c.data = mul!(c.data, a.data, b.data)
  c.data = mod(c.data, defining_polynomial(parent(a)))
  c.precision = compute_precision(a.parent, c.data)
  return c
end

function inv(a::LocalFieldElem)
  K = parent(a)
  p = invmod(a.data, defining_polynomial(K))
  return K(p)
end

Base.:(^)(a::LocalFieldElem, n::UInt) = a^Int(n)

function Base.:(^)(a::LocalFieldElem, n::Int)
  K = parent(a)
  if iszero(a)
    return zero(K)
  end

  if iszero(n)
    return one(K)
  end

  if n < 0 && iszero(a)
    error("Element is not invertible")
  end

  return K(powermod(a.data, n, defining_polynomial(K)))
end

################################################################################
#
#  Promote rules
#
################################################################################

AbstractAlgebra.promote_rule(::Type{S}, ::Type{fmpz}) where S <: LocalFieldElem = S

AbstractAlgebra.promote_rule(::Type{fmpz}, ::Type{S}) where S <: LocalFieldElem = S

AbstractAlgebra.promote_rule(::Type{S}, ::Type{fmpq}) where S <: LocalFieldElem = S

AbstractAlgebra.promote_rule(::Type{fmpq}, ::Type{S}) where S <: LocalFieldElem = S

AbstractAlgebra.promote_rule(::Type{T}, ::Type{S}) where {S <: LocalFieldElem, T <: Integer} = S

AbstractAlgebra.promote_rule(::Type{S}, ::Type{T}) where {S <: LocalFieldElem, T <: Integer} = S

AbstractAlgebra.promote_rule(::Type{LocalFieldElem{padic, T}}, ::Type{padic}) where T <: LocalFieldParameter = LocalFieldElem{padic, T}

AbstractAlgebra.promote_rule(::Type{LocalFieldElem{qadic, T}}, ::Type{qadic}) where T <: LocalFieldParameter = LocalFieldElem{qadic, T}

AbstractAlgebra.promote_rule(::Type{padic}, ::Type{LocalFieldElem{padic, T}}) where T <: LocalFieldParameter = LocalFieldElem{padic, T}

AbstractAlgebra.promote_rule(::Type{qadic}, ::Type{LocalFieldElem{qadic, T}}) where T <: LocalFieldParameter = LocalFieldElem{qadic, T}

function AbstractAlgebra.promote_rule(::Type{LocalFieldElem{S, T}}, ::Type{padic}) where {S <: LocalFieldElem,  T <: LocalFieldParameter}
  if S === AbstractAlgebra.promote_rule(S, padic)
    return LocalFieldElem{S, T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{LocalFieldElem{S, T}}, ::Type{qadic}) where {S <: LocalFieldElem, T <: LocalFieldParameter}
  if S === AbstractAlgebra.promote_rule(S, qadic)
    return LocalFieldElem{S, T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{padic}, ::Type{LocalFieldElem{S, T}}) where {S <: LocalFieldElem,  T <: LocalFieldParameter}
  if S === AbstractAlgebra.promote_rule(S, padic)
    return LocalFieldElem{S, T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{qadic}, ::Type{LocalFieldElem{S, T}}) where {S <: LocalFieldElem, T <: LocalFieldParameter}
  if S === AbstractAlgebra.promote_rule(S, qadic)
    return LocalFieldElem{S, T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{LocalFieldElem{S, T}}, ::Type{LocalFieldElem{U, V}}) where {S <: LocalFieldElem, T <: LocalFieldParameter, U <: LocalFieldElem, V <: LocalFieldParameter}
  if S === U && T === V
    return LocalFieldElem{S, T}
  end
  R = AbstractAlgebra.promote_rule(S, LocalFieldElem{U, V})
  if R === S
    return LocalFieldElem{S, T}
  end
  R = AbstractAlgebra.promote_rule(U, LocalFieldElem{S, T})
  if R === U
    return LocalFieldElem{U, V}
  else
    return Union{}
  end
end


text:
module ContinuumArrays
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays
import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==,
                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
                first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum
import Base.Broadcast: materialize, BroadcastStyle, broadcasted
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport,
                        adjointlayout, LdivApplyStyle, arguments, _arguments, call, broadcastlayout, lazy_getindex,
                        sublayout, ApplyLayout, BroadcastLayout, combine_mul_styles
import LinearAlgebra: pinv
import BandedMatrices: AbstractBandedLayout, _BandedMatrix
import FillArrays: AbstractFill, getindex_value

import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle

export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform

####
# Interval indexing support
####
struct AlephInfinity{N} <: Integer end

==(::AlephInfinity, ::Int) = false
==(::Int, ::AlephInfinity) = false

*(::AlephInfinity{N}, ::AlephInfinity{N}) where N = AlephInfinity{N}()

const ℵ₁ = AlephInfinity{1}()

show(io::IO, F::AlephInfinity{1}) where N =
    print(io, "ℵ₁")

const QMul2{A,B} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B}}
const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}

cardinality(::AbstractInterval) = ℵ₁
*(ℵ::AlephInfinity) = ℵ

Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)

for find in (:findfirst, :findlast)
    @eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
end

function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion)
    r = findfirst(f,d)
    r === nothing ? eltype(d)[] : [r]
end

function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end

BroadcastStyle(::Type{<:Inclusion{<:Any,<:AbstractInterval}}) = LazyQuasiArrayStyle{1}()
BroadcastStyle(::Type{<:QuasiAdjoint{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()
BroadcastStyle(::Type{<:QuasiTranspose{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()

###
# Maps
###

# Affine map represents A*x .+ b
struct AffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T}
    A::AA
    x::X
    b::B
end

AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
    AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)

AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)

AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)

axes(A::AffineQuasiVector) = axes(A.x)
getindex(A::AffineQuasiVector, k::Number) = A.A*A.x[k] .+ A.b
inbounds_getindex(A::AffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
isempty(A::AffineQuasiVector) = isempty(A.x)
==(a::AffineQuasiVector, b::AffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b

BroadcastStyle(::Type{<:AffineQuasiVector}) = LazyQuasiArrayStyle{1}()

for op in(:*, :\, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
end
for op in(:/, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
end

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AffineQuasiVector) = AffineQuasiVector(a, x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AffineQuasiVector) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), a::Number, x::AffineQuasiVector) = AffineQuasiVector(one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, b)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)

function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AffineQuasiVector{<:Real,<:Real,<:Inclusion{<:Real,<:AbstractInterval}})
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end

for find in (:findfirst, :findlast, :findall)
    @eval $find(f::Base.Fix2{typeof(isequal)}, d::AffineQuasiVector) = $find(isequal(d.A \ (f.x .- d.b)), d.x)
end

minimum(d::AffineQuasiVector{<:Real,<:Real,<:Inclusion}) = signbit(d.A) ? last(d) : first(d)
maximum(d::AffineQuasiVector{<:Real,<:Real,<:Inclusion}) = signbit(d.A) ? first(d) : last(d)

union(d::AffineQuasiVector{<:Real,<:Real,<:Inclusion}) = Inclusion(minimum(d)..maximum(d))

include("operators.jl")
include("bases/bases.jl")

end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis
using ConstrainedControl
using LinearAlgebra
using Rotations

# Parameters
ex = [1.0;0.0;0.0]
ey = [0.0;1.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1

p2 = [0.0;0.0;length1/2] # joint connection point

# Desired orientation
ϕ = 0.

# Links
origin = Origin{Float64}()
cart = Box(0.1, 0.5, 0.1, length1/2)
pole = Box(width, depth, length1, length1)

# Constraints
joint1 = EqualityConstraint(Prismatic(origin, cart, ey))
joint2 = EqualityConstraint(Revolute(cart, pole, ex; p2=-p2))

links = [cart;pole]
constraints = [joint1;joint2]

mech = Mechanism(origin, links, constraints, g=-9.81)
setPosition!(origin,cart,Δx = [0;0.5;0])
setPosition!(cart,pole,p2 = -p2,Δq = UnitQuaternion(RotX(ϕ+0.2)))

xd = [[[0;0;0.0]];[[0;0;0.5]]]

Q = [diagm(ones(12))*1.0 for i=1:2]
R = [ones(1,1)]

lqr = LQR(mech, getid.(links), [getid(constraints[1])], Q, R, 10., xd=xd)

storage = simulate!(mech,10,lqr,record = true)
visualize(mech,storage)


text:
module Nonogram
include("encode.jl")
include("checker.jl")
include("constant.jl")
include("generator.jl")
include("transform.jl")

using Pkg
if !in("StatsBase", keys(Pkg.installed()))
    Pkg.add("StatsBase")
end
using Printf
using StatsBase
using Statistics
using .Encode
using .Checker
using .Generator
using .Transform
using .Constant: NONO_SIZE, PROBLEMS, ATTEMPT_TRY, NEW_GEN, PAIRS

"""
New children in new generation
"""
function new_children!(rows::Array{Array{Int,1},1}, eachNumber::Int)::Array{UInt128,2}
    return vcat(Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 0, 0, false),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 0, 0, true),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 1, 0, false),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 1, 0, true),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 0, 1, false),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 0, 1, true),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 1, 1, false),
        Generator.generate_children!(rows, eachNumber, Encode.dist_random_encode!, 1, 1, true),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 0, false, false),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 0, false, true),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 1, false, false),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 1, false, true),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 0, true, false),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 0, true, true),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 1, true, false),
        Generator.generate_children!(rows, eachNumber, Encode.step_random_encode!, 1, true, true))
end

"""
Select children for the best remained
"""
function select_children!(cols::Array{Array{Int,1},1}, children::Array{UInt128,2}, pairs::Int)::Tuple{Array{UInt128,2},Array{Int,1}}
    if pairs * 2 > div(length(children), NONO_SIZE)
        throw(ArgumentError("Pairs to be got are over children's amount!"))
    end
    fitnesses = [(i, fitness) for (i, fitness) in enumerate([Checker.get_fitness(cols, Transform.transform_column!(child)) for child in eachrow(children)])]
    sort!(fitnesses, by = x->x[2], rev = true)
    return children[[i[1] for i in fitnesses[1:(pairs * 2)]], :], [i[2] for i in fitnesses[1:(pairs * 2)]]
end

"""
Crossover the children gene in best children
"""
function crossover_children!(cols::Array{Array{Int,1},1}, children::Array{UInt128,2})::Tuple{Array{UInt128,2}, Array{Int,1}}
    crossedChildren = Array{UInt128,1}(undef, 0)
    fitnesses = Array{Int,1}(undef, 0)
    divideRange = rand(2:(NONO_SIZE - 1))
    for i in 1:size(children, 1)
        for j in 1:size(children, 1)
            if i == j
                continue
            end
            cated = vcat(children[i, 1:divideRange], children[j, divideRange + 1:end])
            push!(crossedChildren, cated...)
            push!(fitnesses, Checker.get_fitness(cols, Transform.transform_column!(cated)))
        end
    end
    crossedChildren = transpose(reshape(crossedChildren, NONO_SIZE, :))
    fitnesses = [(i, fitness) for (i, fitness) in enumerate(fitnesses)]
    sort!(fitnesses, by = x->x[2], rev = true)
    return crossedChildren[[i[1] for i in fitnesses], :], [i[2] for i in fitnesses]
end

"""
Mutate the best children
"""
function mutate_children!(rows::Array{Array{Int,1},1}, cols::Array{Array{Int,1},1}, children::Array{UInt128,2})::Tuple{Array{UInt128,2}, Array{Int,1}}
    fitnesses = Array{Int,1}(undef, 0)
    len = size(children, 1)
    for i in 1:len
        for mutateIndex in sample(1:NONO_SIZE, 1, replace = false)
            children[i, mutateIndex] = Encode.dist_random_encode!(rows[mutateIndex], rand(0:1), 1, rand(Bool))[1]
        end
        push!(fitnesses, Checker.get_fitness(cols, Transform.transform_column!(children[i, :])))
    end
    fitnesses = [(i, fitness) for (i, fitness) in enumerate(fitnesses)]
    sort!(fitnesses, by = x->x[2], rev = true)
    return children[[i[1] for i in fitnesses[1:div(len, 3)]], :], [i[2] for i in fitnesses[1:div(len, 3)]]
end

function solve!(problem::SubArray)
    # input cols from problem
    cols = [[parse(Int, n) for n in split(col)] for col in problem[1:NONO_SIZE]]
    # input rows from problem
    rows = [[parse(Int, n) for n in split(row)] for row in problem[NONO_SIZE + 1:end]]
    recordsChildren = undef
    newFile = open("new.log", "w+")
    crossoverFile = open("crossover.log", "w+")
    mutateFile = open("mutate.log", "w+")
    for attempt in 1:ATTEMPT_TRY
        # generate children and get the best
        if recordsChildren == undef
            children, fitnesses = select_children!(cols, new_children!(rows, NEW_GEN), PAIRS)
        else
            children, fitnesses = select_children!(cols, vcat(new_children!(rows, NEW_GEN), recordsChildren), PAIRS)
        end
        if fitnesses[1] == 0
            println(@sprintf("Found solution in %d generation in select_children", attempt))
            break
        end
        # debug feature
        println(@sprintf("Try time: %d, in new Fitness: %d", attempt, fitnesses[1]))
        println(newFile, fitnesses[1], ',')
        flush(newFile)
        children, fitnesses = crossover_children!(cols, children)
        if fitnesses[1] == 0
            println(@sprintf("Found solution in %d generation in crossover_children", attempt))
            break
        end
        # debug feature
        println(@sprintf("Try time: %d, in crossover Fitness: %d", attempt, fitnesses[1]))
        println(crossoverFile, fitnesses[1], ',')
        flush(crossoverFile)
        recordsChildren, fitnesses = mutate_children!(rows, cols, children)
        if fitnesses[1] == 0
            println(@sprintf("Found solution in %d generation in mutate_children", attempt))
            break
        end
        # debug feature
        println(@sprintf("Try time: %d, in mutate Fitness: %d", attempt, fitnesses[1]))
        println(mutateFile, fitnesses[1], ',')
        flush(mutateFile)
    end
    close(newFile)
    close(crossoverFile)
    close(mutateFile)
end

inputLines = Array{String,1}(undef, 0)
@timev open("nonogram_input.txt", "r") do file
    # problem count and line count
    pc, lc = 0, 0
    while !eof(file)
        line = strip(readline(file))
        if isempty(line)
            continue
        end
        if startswith(line, '$')
            lc = NONO_SIZE * 2
            if pc == PROBLEMS
                break
            end
            pc += 1
        elseif lc > 0
            lc -= 1
            push!(inputLines, line)
        end
    end
end
# reshape lines into problem size (NONO_SIZE * 2, PROBLEM)
problems = reshape(inputLines, NONO_SIZE * 2, :)
@timev for (index, problem) in enumerate(eachcol(problems))
    println("Used $(@sprintf("%.6f seconds in problem %04d", @elapsed(solve!(problem)), index)).")
end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module LazyArtifacts

# reexport the Artifacts API
using Artifacts: Artifacts,
       artifact_exists, artifact_path, artifact_meta, artifact_hash,
       select_downloadable_artifacts, find_artifacts_toml, @artifact_str
export artifact_exists, artifact_path, artifact_meta, artifact_hash,
       select_downloadable_artifacts, find_artifacts_toml, @artifact_str

# define a function for satisfying lazy Artifact downloads
using Pkg.Artifacts: ensure_artifact_installed

end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
import FiniteVolumes.CartesianMeshes.AbstractCartesianMesh

struct DiffusionFlux end
const ∇ = DiffusionFlux()

function (::Centered)(::DiffusionFlux, mesh::AbstractCartesianMesh, w, i_face)
    n = FiniteVolumes.normal_vector(mesh, i_face)
    if mesh isa AbstractCartesianMesh{1}
        face_dir = 1
    elseif mesh isa AbstractCartesianMesh{2}
        face_dir = FiniteVolumes.CartesianMeshes._direction(i_face)
    end
    i_cell_1, i_cell_2 = FiniteVolumes.cells_next_to_inner_face(mesh, i_face)
    return (w[i_cell_2] - w[i_cell_1])*n[face_dir]/FiniteVolumes.CartesianMeshes.dx(mesh)[face_dir]
end


text:
### Make lexicon file from Corpus to then make MinimalistGrammar
#

##################CREATE POSSIBLE FEATURE SETS GIVE 1 CAT###########################
function add_feats(feature_list, sub_list,category_feats, k)
    append!(feature_list, sub_list)
    if k == 0
        return feature_list
    end
    new_sub_list = []
    for list in sub_list
        for c in category_feats
            Rfeats = string(list[1:(end-1)], "R=", c, " ", list[end])
            Lfeats = string("L=", c, " ", list)
            push!(new_sub_list, Rfeats)
            push!(new_sub_list, Lfeats)
        end
    end
    k-=1
    return add_feats(feature_list,new_sub_list,category_feats,k)
end

function add_feats_Rightfirst(feature_list, sub_list,category_feats, k)
    append!(feature_list, sub_list)
    if k == 0
        return feature_list
    end
    new_sub_list = []
    for list in sub_list
        for c in category_feats
            Rfeats = string("R=", c, " ", list)
            push!(new_sub_list, Rfeats)
        end
    end
    k-=1
    return add_feats_Rightfirst(feature_list,new_sub_list,category_feats,k)
end

function add_feats_Leftsecond(feature_list, sub_list,category_feats, k, m)
    append!(feature_list, sub_list)
    if k == 0
        return feature_list
    end
    new_sub_list = []
    for list in sub_list
        if length(list) < 4*m+1
            for c in category_feats
                Lfeats = string("L=", c, " ", list)
                push!(new_sub_list, Lfeats)
            end
        end
    end
    k-=1
    return add_feats_Leftsecond(feature_list,new_sub_list,category_feats,k, m)
end

function add_feats_Balanced(feature_list, sub_list,category_feats, k, m)
    right_list = add_feats_Rightfirst(feature_list, sub_list, category_feats, k)
    left_list = add_feats_Leftsecond([], right_list, category_feats, k, m)
    return left_list
end

###################CREATE FEATURE SET DICT FOR ALL CATEGORIES###################
function make_feature_list_by_cat_dict(categories)
    feature_list = Dict()
    cats = []
    for (c, sels, k) in categories
        push!(cats, c)
        sub_list = [c]
        feature_list[c] = add_feats_Balanced([], sub_list , sels , k, k)
    end
    #Add start symbol category
    feature_list["c"] = add_feats_Rightfirst([], ["c"] , cats, 1)
    return feature_list
    #feature_list_start = add_feats_Balanced([], ["c"] , cats, 2, 2)
    #return (feature_list, feature_list_start)
end

function make_feature_list_by_cat(categories, k)
    feature_list = Dict()
    for c in categories
        sub_list = [c]
        feature_list[c] = add_feats([], sub_list , categories, k)
    end
    # Add start symbol category
    feature_list_start = add_feats([], ["c"] , categories, k)
    return (feature_list, feature_list_start)
end

function make_feature_list(categories, k)
    feature_list = add_feats([], categories , categories, k)
    # Add start symbol category
    feature_list_start = add_feats([], ["c"] , categories, k)
    return (feature_list, feature_list_start)
end
######################TRANSFORM POS TO MINIMALIST CAT###########################

function getcat_fromPOS(pos)
    cat = ""
    if in(pos, ["NOUN", "PROPN"])
        cat = "n"
    elseif in(pos, ["DET", "PRON"])
        cat = "d"
    elseif in(pos, ["VERB"])
        cat = "v"
    elseif in(pos, ["ADJ"])
        cat = "a"
    elseif in(pos, ["ADV"])
        cat = "x"
    elseif in(pos, ["ADP"])
        cat = "p"
    elseif in(pos, ["AUX"])
        cat = "t"
    elseif in(pos, ["SCONJ"])
        cat = "s"
    elseif in(pos, ["PART"])
        cat = "r"
    end
    return cat
end

## In the case of English, this splits the POS PART into three category types:
## infinitive-> t, possessive-> d, negation-> neg
function getcat_fromPOS_EnglishOnly(lemma, pos)
    cat = ""
    if in(pos, ["NOUN", "PROPN"])
        cat = "n"
    elseif in(pos, ["DET", "PRON"])
        cat = "d"
    elseif in(pos, ["VERB"])
        cat = "v"
    elseif in(pos, ["ADJ"])
        cat = "a"
    elseif in(pos, ["ADV"])
        cat = "x"
    elseif in(pos, ["ADP"])
        cat = "p"
    elseif in(pos, ["AUX"])
        cat = "t"
    elseif in(pos, ["SCONJ"])
        cat = "s"
    elseif in(pos, ["PART"])
        if contains(lemma, "n")
            cat = "neg"
        elseif contains(lemma, "t")
            cat = "t"
        else
            cat = "d"
        end
    end
    return cat
end

######################CREATE A LEXICON FROM CORPUS##############################

## for testing, corpus is a list of sentences, which is itself a list of tuples,
## (word, cat)
function makeLexicon_test(corpus, categories, k)

    feature_list, feature_list_start = make_feature_list_by_cat(categories, k)
    words = []
    lexicon = ""
    for featlist in feature_list_start
        item = string( ":: ", featlist, " ; 1")
        lexicon = string(lexicon,"\n", item)
    end

    for sent in corpus
        for (word,cat) in sent
            if !(word in words)
                push!(words, word)
                for featlist in feature_list[cat]
                    item = string(word, " :: ", featlist, " ; 1")
                    lexicon = string(lexicon,"\n", item)
                end
            end
        end
    end
    open("lexicon.txt", "w") do f
        write(f, lexicon)
    end
    return lexicon
end

## Makes a lexicon and filtered corpus from UD data. If using the english only
## category mapping add eng_only=true to call.
function makeLexicon_andCorpus(UDcorpus, categories, eng_only=false)

    feature_list = make_feature_list_by_cat_dict(categories)
    words = []
    lexicon = ""
    corpus = []
    counts = Dict()
    for (cat, sels, n) in categories
        counts[cat] = 1
    end
    for (cat,featlist) in feature_list
        for flist in featlist
            item = string( ":: ", flist, " ; 1")
            lexicon = string(lexicon,"\n", item)
        end
    end
    for sent in UDcorpus.sents
        sentence = []
        goodsent = true
        for word in sent.words
            w = lowercase(word.form)
            ##remove sentences which contain one of these POS tags
            if in(word.upos_tag, ["INTJ", "CONJ","SYM","X"])
                goodsent = false
                break
            elseif ';' in word.form
                goodsent = false
                break
            elseif in(word.upos_tag, ["ADJ","ADV","ADP","AUX","DET", "NOUN", "PART", "PRON", "PROPN", "SCONJ", "VERB"])
                if !(string(w,word.upos_tag) in words)
                    push!(words, string(w,word.upos_tag))
                    if eng_only
                        cat = getcat_fromPOS_EnglishOnly(word.lemma, word.upos_tag)
                    else
                        cat = getcat_fromPOS(word.upos_tag)
                    end
                    counts[cat] +=1
                    for featlist in feature_list[cat]
                        item = string(w, " :: ", featlist, " ; 1")
                        lexicon = string(lexicon,"\n", item)
                    end
                    sentence = push!(sentence, w)
                else
                    sentence = push!(sentence, w)
                end

            ### ADDING SECTION TO INCLUDE NUMERALS
            elseif in(word.upos_tag, ["NUM"])
                if !(string(w,word.upos_tag) in words)
                    push!(words, string(w,word.upos_tag))
                    cats = ["d","n","a"]
                    for cat in cats
                        counts[cat] +=1
                        for featlist in feature_list[cat]
                            item = string(w, " :: ", featlist, " ; 1")
                            lexicon = string(lexicon,"\n", item)
                        end
                    end
                    sentence = push!(sentence, w)
                else
                    sentence = push!(sentence, w)
                end
            end
        end
        if goodsent && sentence != []
            push!(corpus, sentence)
        end
    end
    open("lexicon.txt", "w") do f
        write(f, lexicon)
    end
    open("corpus.txt", "w") do f
        for sentence in corpus
            for word in sentence
                write(f, string(word, " "))
            end
            write(f, "\n")
        end
    end
    corpus = convert(Array{Array{String,1},1},corpus)
    return (lexicon, corpus, counts)
end

########## FUNCTIONS FOR FIXING GRAMMARS ##########

# remove duplicate lexical items
function remove_duplicates!(grammar)
    for word in keys(grammar.lexicon)
        ids = unique(map(GeneralizedChartparsing.id, grammar.lexicon[word]))
        if length(ids) != length(grammar.lexicon[word])
            grammar.lexicon[word] = map(LexicalItem, ids)
        end
    end
end

# remove null heads of the form =x x (vacuous heads)
function remove_vacuous_nullheads!(grammar)
    function nonvacuous(li)
        feats = li.features
        if length(feats) == 2
            GeneralizedChartparsing.name(feats[1]) != GeneralizedChartparsing.name(feats[2])
        else true
        end
    end
    grammar.lexicon[""] = filter(nonvacuous, grammar.lexicon[""])
end


text:
#=
# A tetrehedral reference element.
# Notation follows that used in Nodal Discontinuous Galerkin Methods
# by Hesthaven and Warburton.
#  https://link.springer.com/book/10.1007/978-0-387-72067-8
#
=#

function build_tetrahedron_refel(refel)
    # refel has already been created, but needs quadrature matrices
    (refel.V, refel.Dr, refel.Ds, refel.Dt) = tetrahedron_vandermonds(refel, refel.r);
    refel.invV = inv(refel.V);

    (refel.Vg, DVgr, DVgs, DVgt) = tetrahedron_vandermonds(refel, refel.g);
    #refel.invVg = inv(refel.Vg);

    refel.Q = refel.Vg*refel.invV;
    refel.Qr = DVgr*refel.invV;
    refel.Qs = DVgs*refel.invV;
    refel.Qt = DVgt*refel.invV;
    refel.Ddr = refel.Dr*refel.invV;
    refel.Dds = refel.Ds*refel.invV;
    refel.Ddt = refel.Dt*refel.invV;

    return refel;
end

function tetrahedron_vandermonds(refel, r)
    Np = refel.Np;      # number of nodal points
    Nrp = size(r,1);    # number of quadrature points may be different
    V = zeros(Nrp, Np);
    gradVr = zeros(Nrp, Np);
    gradVs = zeros(Nrp, Np);
    gradVt = zeros(Nrp, Np);

    # Transfer (r,s,t) to (a,b,c) coordinates
    a = zeros(Nrp);
    b = zeros(Nrp);
    for ni=1:Nrp
        if r[ni,2]+r[ni,3] != 0
            a[ni] = 2*(1+r[ni,1])/(-r[ni,2]-r[ni,3])-1;
        else
            a[ni] = -1;
        end
        if r[ni,3] != 1
            b[ni] = 2*(1+r[ni,2])/(1-r[ni,3])-1;
        else
            b[ni] = -1;
        end
    end
    c = r[:,3];

    # build the Vandermonde and gradVandermond matrix
    sk = 1;
    for i=0:refel.N
        h1 = jacobi_polynomial(a,0,0,i);
        dfa = grad_jacobi_polynomial(a, 0, 0, i);
        for j=0:(refel.N - i)
            h2 = jacobi_polynomial(b,2*i+1,0,j);
            dgb = grad_jacobi_polynomial(b, 2*i+1,0, j);
            for k=0:(refel.N - i - j)
                h3 = jacobi_polynomial(c,2*(i+j)+2,0,k);
                dhc = grad_jacobi_polynomial(b, 2*(i+j)+2,0, k);

                V[:,sk] = 2.8284271247461903 .* h1 .* h2 .* h3 .* (1 .- b).^i .* (1 .- c).^(i+j);

                # r-derivative
                dmodedr = dfa.*h2.*h3;
                if i>0
                    dmodedr = dmodedr.*((0.5 .* (1 .- b)).^(i-1));
                end
                if (i+j)>0
                    dmodedr = dmodedr.*((0.5 .* (1 .- c)).^(i+j-1));
                end

                # s-derivative
                dmodeds = dmodedr .* (0.5 .* (1 .+ a));
                tmp = dgb .* (0.5 .* (1 .- b)).^i;
                if i>0
                    tmp = tmp .+ (-0.5*i) .* (h2 .* (0.5 .* (1 .- b)).^(i-1));
                end
                if i+j>0
                    tmp = tmp .* (0.5 .* (1 .- c)).^(i+j-1);
                end
                tmp = tmp .* h3 .* h1;
                dmodeds = dmodeds .+ tmp;

                # t-derivative
                dmodedt = 0.5 .* (1 .+ a) .* dmodedr .+ 0.5 .* (1 .+ b) .* tmp;
                tmp = dhc .* (0.5 .* (1 .- c)).^(i+j);
                if i+j>0
                    tmp = tmp .- (0.5*(i+j)) .* h3 .* (0.5 .* (1 .- c)).^(i+j-1);
                end
                tmp = h1 .* (h2 .* tmp) .* (0.5 .* (1 .- b)).^i;
                dmodedt = dmodedt .+ tmp;

                # Normalize
                dmodedr = 2^(2*i+j+1.5) .* dmodedr;
                dmodeds = 2^(2*i+j+1.5) .* dmodeds;
                dmodedt = 2^(2*i+j+1.5) .* dmodedt;

                gradVr[:,sk] = dmodedr;
                gradVs[:,sk] = dmodeds;
                gradVt[:,sk] = dmodedt;

                sk = sk+1;
            end
        end
    end

    return (V, gradVr, gradVs, gradVt);
end


text:
"""
    SeisAddNoise(d, snr; <keyword arguments>)

Add noise at a given signal-to-noise ratio level `snr` to an N-dimensional input
data `d`. Noise can be band limited using kewyord `L`.

# Arguments
* `d::Array{Real, N}`: N-dimensional data.
* `snr::Real`: signal-to-noise ratio.

# Keyword arguments
* `db::Bool=false`: `db=false` if `snr` is given by amplitude, `db=true` if
snr is given in dB.
* `pdf::AbstractString="gaussian"`: random noise probability distribution:
`"gaussian"` or `"uniform"`.
* `L::Int=1`: averaging operator length to band-limit the random noise.

# Examples
```
julia> w = Ricker(); wn = SeisAddNoise(w, 2); plot(w); plot(wn);
MeasureSNR(w, wn)

julia> d, extent = SeisHypEvents(); dn = SeisAddNoise(d, 1.0, db=true, L=9);
SeisPlot([d dn], extent); MeasureSNR(d, dn, db=true)
```
Credits: Juan I. Sabbione, 2016
"""

function SeisAddNoise{T<:Real, N}(d::Array{T, N}, snr::Real; db::Bool=false,
                                  pdf::AbstractString="gaussian", L::Int=1)

    noise = GenNoise(size(d), pdf, L=L)
    if db==false
        noise = noise/vecnorm(noise) * vecnorm(d)/snr
    elseif db==true
        noise = noise/vecnorm(noise) * vecnorm(d)/10.0^(0.05*snr)
    end
    noisy = d + noise
    assert(abs(MeasureSNR(d, noisy, db=db) - snr) < 1e10*eps(AbstractFloat(snr)))
    return noisy

end

# Generate the noise
function GenNoise(dims::Tuple, pdf::AbstractString; L::Int=1)
    N = length(dims)
    n1 = dims[1]
    N == 1 ? nx = 1 : nx = *(dims[2:end]...)
    hamm = Hamming(L)
    L2 = floor(Int,L/2)
    if pdf=="gaussian"
        noise = randn(n1, nx)
    elseif pdf=="uniform"
        noise = -1 + 2*rand(n1, nx)
    else
        error("pdf must be gussian or uniform")
    end
    for j = 1:nx
     noise[:,j] = conv(noise[:,j], hamm)[L2+1:n1+L2]
    end
    reshape(noise, dims...)
end


text:
using Documenter
using Starlight

makedocs(
    sitename = "Starlight",
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    modules = [Starlight],
    pages = [
      "Home" => "index.md",
      "Guide" => Any[
        "Walkthrough" => Any[
          "guide/walkthrough/intro.md",
          "guide/walkthrough/getting-started.md",
          "guide/walkthrough/app-lifecycle.md",
          "guide/walkthrough/message-passing.md",
          "guide/walkthrough/clock.md",
          "guide/walkthrough/ecs.md"
        ],
        "Telescope" => Any[
          "guide/telescope/intro.md",
          "guide/telescope/rendering.md",
          "guide/telescope/input.md",
          "guide/telescope/physics.md"
        ]
      ],
      "Examples" => Any[
        "examples/pong.md"
      ],
      "API" => "api.md"
    ]
)

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


text:
using FATEStools
using Documenter

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

makedocs(;
    modules=[FATEStools],
    authors="Gregory Lemieux, Lawrence Berkeley National Laboratory",
    repo="https://github.com/glemieux/FATEStools.jl/blob/{commit}{path}#{line}",
    sitename="FATEStools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://glemieux.github.io/FATEStools.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/glemieux/FATEStools.jl",
    devbranch="main",
)


text:
using .MyExample2
using Test

@testset "MyExample2.jl" begin
    # 2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13

end


text:
"""
    F(a, c0, Finf, c, G; F0=0.)
"""
function F(a, c0, Finf, c, G; F0=0.)
    F0 .+ a .* log.( c/c0 ) .- G*Finf
end

"""
    F(m::ClimateModel; M=false, R=false, G=false)
"""
function F(m; M=false, R=false, G=false, F0=0.)
    return F(
        m.physics.a, m.physics.c0, m.economics.Finf,
        c(m, M=M, R=R),
        m.controls.geoeng .* (1. .- .~past_mask(m) * ~G),
        F0=F0
    )
end

F2x(a::Float64) = a*log(2)
F2x(m::ClimateModel) = F2x(m.physics.a)

ECS(a, B) = F2x(a)/B
ECS(params::ClimateModelParameters) = ECS(params.physics.a, m.physics.B)
ECS(m::ClimateModel) = ECS(m.physics.a, m.physics.B)

calc_B(a::Float64, ECS::Float64) = F2x(a)/ECS
calc_B(params::ClimateModelParameters; ECS=ECS(params)) = calc_B(params.physics.a, ECS)
calc_B(m::ClimateModel; ECS=ECS(m)) = calc_B(m.physics.a, ECS)

τd(Cd, B, κ) = (Cd/B) * (B+κ)/κ
τd(phys::Physics) = τd(phys.Cd, phys.B, phys.κ)
τd(m::ClimateModel) = τd(m.physics)

"""
    T_fast(F, κ, B)
"""
T_fast(F, κ, B) = F/(κ + B)

"""
    T_fast(m::ClimateModel; M=false, R=false, G=false)
"""
T_fast(m::ClimateModel; M=false, R=false, G=false) = T_fast(
    F(m, M=M, R=R, G=G),
    m.physics.κ,
    m.physics.B
)

"""
    T_slow(F, Cd, κ, B, t, dt)
"""
function T_slow(F, Cd, κ, B, t, dt)
    τ = τd(Cd, κ, B)
    return (
        (κ/B) / (κ + B) *
        exp.( - (t .- (t[1] - dt)) / τ) .*
        cumsum( (exp.( (t .- (t[1] - dt)) / τ) / τ) .* F * dt)
    )
end

"""
    T_slow(m::ClimateModel; M=false, R=false, G=false)
"""
T_slow(m::ClimateModel; M=false, R=false, G=false) = T_slow(
    F(m, M=M, R=R, G=G),
    m.physics.Cd,
    m.physics.κ,
    m.physics.B,
    t(m),
    m.domain.dt
)

"""
    T(T0, F, Cd, κ, B, t, dt)

Returns the sum of the initial, fast mode, and slow mode temperature change.

# Arguments
- `T0::Float64`: warming relative to pre-industrial [°C]
- `F::Array{Float64}`: change in radiative forcing since the initial time ``t_{0}`` [W/m``{2}``]
- `Cd::Float64`: deep ocean heat capacity [W yr m``^{2}`` K``^{-1}``]
- `κ::Float64`: ocean heat uptake rate [W m``^{2}`` K``^{-1}``]
- `B::Float64`: feedback parameter [W m``^{2}`` K``^{-1}``]
- `t::Array{Float64}`: year [years]
- `dt::Float64`: timestep [years]

"""
T(T0, F, Cd, κ, B, t, dt) = (
    T0 .+
    T_fast(F, κ, B) .+
    T_slow(F, Cd, κ, B, t, dt)
)

"""
    T(m::ClimateModel; M=false, R=false, G=false)

Returns the sum of the initial, fast mode, and slow mode temperature change,
as diagnosed from `m` and modified by the climate controls activated by the
Boolean kwargs.
"""
T(m::ClimateModel; M=false, R=false, G=false) = T(
    m.physics.T0,
    F(m, M=M, R=R, G=G),
    m.physics.Cd,
    m.physics.κ,
    m.physics.B,
    t(m),
    m.domain.dt
)


text:
using TextWrap

#Format the executed document
function format(executed, doctype)
  formatted = String[]
  docformat = formats[doctype]
  #@show docformat

  #Complete format dictionaries with defaults
  formatdict = docformat.formatdict
  get!(formatdict, :termstart, formatdict[:codestart])
  get!(formatdict, :termend, formatdict[:codeend])
  get!(formatdict, :out_width, nothing)
  get!(formatdict, :out_height, nothing)
  get!(formatdict, :fig_pos, nothing)
  get!(formatdict, :fig_env, nothing)

  for chunk in copy(executed)

      if chunk[:type] == "doc"
          push!(formatted, chunk[:content])
      else
          #Fill undefined options with format specific defaults
          chunk[:out_width] == nothing && (chunk[:out_width] =  docformat.formatdict[:out_width])
          chunk[:fig_pos] == nothing && (chunk[:fig_pos] =  docformat.formatdict[:fig_pos])

          #Only use floats if chunk has caption or sets fig_env
          if chunk[:fig_cap] != nothing && chunk[:fig_env] == nothing
              (chunk[:fig_env] =  docformat.formatdict[:fig_env])
          end

          #Format code
          result = format_codechunk(chunk, formatdict)
          #Handle figures
          if chunk[:fig] && length(chunk[:figure]) > 0
              if chunk[:include]
                  result *= formatfigures(chunk, docformat)
              end
          end
          push!(formatted, result)
      end
  end

 return formatted
end

function format_codechunk(chunk, formatdict)
    if haskey(formatdict, :indent)
        chunk[:content] = indent(chunk[:content], formatdict[:indent])
    end

    if !chunk[:eval]
        if chunk[:echo]
            result = "$(formatdict[:codestart])$(chunk[:content])\n$(formatdict[:codeend])"
            return result
        else
            r = ""
            return r
        end
    end

    if chunk[:term]
        result = format_termchunk(chunk, formatdict)
    else
        if chunk[:wrap]
          chunk[:result] = "\n" * wrap(chunk[:result], replace_whitespace=false)
        end

        if chunk[:echo]
            result = "$(formatdict[:codestart])$(chunk[:content])\n$(formatdict[:codeend])\n"
        else
            result = ""
        end

        if (strip(chunk[:result])!= "") && (chunk[:results] != "hidden")
            #@show chunk
            if chunk[:results] != "markup"
                result *= "$(chunk[:result])"
            elseif chunk[:results] == "markup"
                if haskey(formatdict, :indent)
                    chunk[:result] = indent(chunk[:result], formatdict[:indent])
                end
                result *= "$(formatdict[:outputstart])$(chunk[:result])\n$(formatdict[:outputend])\n"
            end
        end

    end

    return result

end

function format_termchunk(chunk, formatdict)
    if chunk[:echo] && chunk[:results] != "hidden"
        result = "$(formatdict[:termstart])$(chunk[:result])\n"
        #@show chunk[:term_state]
        chunk[:term_state] == :text && (result*= "$(formatdict[:termend])\n")
    else
        result = ""
    end
    return result
end

function indent(text, nindent)
    return join(map(x->
                    string(repeat(" ", nindent), x), split(text, "\n")), "\n")
end

type Tex
    description::String
    formatdict::Dict{Symbol,Any}
end

const tex = Tex("Latex with custom code environments",
                @compat Dict{Symbol,Any}(:codestart => "\\begin{juliacode}",
                                         :codeend => "\\end{juliacode}",
                                         :outputstart => "\\begin{juliaout}",
                                         :outputend => "\\end{juliaout}",
                                         :termstart => "\\begin{juliaterm}",
                                         :termend => "\\end{juliaterm}",
                                         :fig_ext => ".pdf",
                                         :extension =>"tex",
                                         :out_width=> "\\linewidth",
                                         :fig_env=> "figure",
                                         :fig_pos => "htpb",
                                         :doctype => "tex"
                                         ))

const texminted = Tex("Latex using minted for highlighting",
                      @compat Dict{Symbol,Any}(
                                         :codestart => "\\begin{minted}[mathescape, fontsize=\\small, xleftmargin=0.5em]{julia}",
                                         :codeend => "\\end{minted}",
                                         :outputstart => "\\begin{minted}[fontsize=\\small, xleftmargin=0.5em, mathescape, frame = leftline]{text}",
                                         :outputend => "\\end{minted}",
                                         :termstart=> "\\begin{minted}[fontsize=\\footnotesize, xleftmargin=0.5em, mathescape]{julia}",
                                         :termend => "\\end{minted}",
                                         :fig_ext => ".pdf",
                                         :extension =>"tex",
                                         :out_width => "\\linewidth",
                                         :fig_env=> "figure",
                                         :fig_pos => "htpb",
                                         :doctype => "texminted"
                                         ))

type Markdown
    description::String
    formatdict::Dict{Symbol,Any}
end

const pandoc = Markdown("Pandoc markdown",
                        @compat Dict{Symbol,Any}(
                                :codestart => "~~~~{.julia}",
                                :codeend=>"~~~~~~~~~~~~~\n\n",
                                :outputstart=>"~~~~{.julia}",
                                :outputend=>"~~~~~~~~~~~~~\n\n",
                                :fig_ext=>".png",
                                :extension=>"md",
                                :doctype=>"pandoc"
                                               ))

const github = Markdown("Github markdown",
                        @compat Dict{Symbol,Any}(
                                :codestart => "````julia",
                                :codeend=> "````\n\n",
                                :outputstart=> "````julia",
                                :outputend=> "````\n\n",
                                :fig_ext=> ".png",
                                :extension=> "md",
                                :doctype=> "github"
                                               ))

type Rest
    description::String
    formatdict::Dict{Symbol,Any}
end

const rst = Rest("reStructuredText and Sphinx",
                 @compat Dict{Symbol,Any}(
                                :codestart => ".. code-block:: julia\n",
                                :codeend => "\n\n",
                                :outputstart => "::\n",
                                :outputend => "\n\n",
                                :indent=> 4,
                                :fig_ext => ".png",
                                :extension => "rst",
                                :out_width => "15 cm",
                                :doctype => "rst"
                                ))

type AsciiDoc
    description::String
    formatdict::Dict{Symbol,Any}
end

#asciidoc -b html5 -a source-highlighter=pygments ...
const adoc = AsciiDoc("AsciiDoc",
        @compat Dict{Symbol,Any}(
        :codestart => "[source,julia]\n--------------------------------------",
        :codeend => "--------------------------------------\n\n",
        :outputstart => "--------------------------------------",
        :outputend => "--------------------------------------\n\n",
        :fig_ext => ".png",
        :extension => "txt",
        :out_width => "600",
        :doctype => "asciidoc"
))

function formatfigures(chunk, docformat::Tex)
    fignames = chunk[:figure]
    caption = chunk[:fig_cap]
    width = chunk[:out_width]
    f_pos = chunk[:fig_pos]
    f_env = chunk[:fig_env]
    result = ""
    figstring = ""

    if f_env != nothing
        result *= """\\begin{$f_env}[$f_pos]\n"""
    end

    for fig = fignames
        figstring *= "\\includegraphics[width=$width]{$fig}\n"
    end

    # Figure environment
    if caption != nothing
        result *= string("\\center\n",
                         "$figstring",
                         "\\caption{$caption}\n")
    else
        result *= figstring
    end

    if chunk[:name] != nothing && f_env !=nothing
        label = chunk[:name]
        result *= "\\label{fig:$label}\n"
    end

    if f_env != nothing
        result *= "\\end{$f_env}\n"
    end

   return result
end

function formatfigures(chunk, docformat::Markdown)
    fignames = chunk[:figure]
    caption = chunk[:fig_cap]
    result = ""
    figstring = ""

    length(fignames) > 0 || (return "")

    if caption != nothing
        result *= "![$caption]($(fignames[1]))\n"
        for fig = fignames[2:end]
            result *= "![]($fig)\n"
            println("Warning, only the first figure gets a caption\n")
        end
    else
        for fig in fignames
            result *= "![]($fig)\n"
        end
    end
    return result
end

function formatfigures(chunk, docformat::Rest)
    fignames = chunk[:figure]
    caption = chunk[:fig_cap]
    width = chunk[:out_width]
    result = ""
    figstring = ""

    for fig=fignames
        figstring *= @sprintf(".. image:: %s\n   :width: %s\n\n", fig, width)
    end

    if caption != nothing
        result *= string(".. figure:: $(fignames[1])\n",
                         "   :width: $width\n\n",
                         "   $caption\n\n")
    else
        result *= figstring
        return result
    end
end

function formatfigures(chunk, docformat::AsciiDoc)
    fignames = chunk[:figure]
    caption = chunk[:fig_cap]
    width = chunk[:out_width]
    result = ""
    figstring = ""

    for fig=fignames
        figstring *= @sprintf("image::%s[width=%s]\n", fig, width)
    end

    if caption != nothing
        result *= string("image::$(fignames[1])",
        "[width=$width,",
        "title=\"$caption\"]")
    else
        result *= figstring
        return result
    end
end

#Add new supported formats here
const formats = @compat Dict{String, Any}("tex" => tex,
                                          "texminted" => texminted,
                                          "pandoc" => pandoc,
                                          "github" => github,
                                          "rst" => rst,
                                          "asciidoc" => adoc
                                          )


text:
using HTTP
using Gumbo
using Cascadia

cache_path = "$(@__DIR__)/../assets/data.html.asset" # ".asset" is used to supress vscode and github wrong hint

function get_from_wiki(;cache=true)
    url = "https://hoi4.paradoxwikis.com/Land_units"
    res = HTTP.request("GET", url)
    data = String(res.body)

    if cache
        f = open(cache_path, "r")
        write(f, data)
        close(f)
    end

    return data
end

function get_from_cache()
    f = open(cache_path, "r")
    data = String(read(f))
    close(f)

    return data
end

function _unit_parse(::Type{Float64}, x::AbstractString)
    if length(x)==0 # empty string "" => 0
        return 0.0
    end

    # weired wiki "107.2.0" => "107.2"
    split_x = split(x, ".")
    if length(split_x) == 1
        x = split_x[1]
    else
        x = join(split_x[1:2], ".")
    end

    if x[1]=='+'
        x=x[2:end]
    end
    if x[end]=='%'
        return Base.parse(Float64, x[1:end-1]) / 100
    end
    return Base.parse(Float64, x)
end

function _unit_parse(::Type{Int64}, x::AbstractString)
    if length(x)==0
        return 0
    end
    if x[1]=='+'
        x=x[2:end]
    end
    return Base.parse(Int64, x)
end

function unit_parse(::Type{String}, x::AbstractString)
    return x
end

function unit_parse(typ, x::AbstractString)
    #global y
    y = split(x, "/")[1]
    _unit_parse(typ, y)
end

function extract_td(td)
    #global td_ = td
    match_img = eachmatch(sel"img", td);
    if length(match_img) > 0
        raw = attrs(match_img[1])["src"]
    else
        raw_list = []
        for child in td.children
            if child isa HTMLElement{:br}
                push!(raw_list, "/")
            else
                push!(raw_list, strip(nodeText(child)))
            end
        end
        raw = join(raw_list, "")
    end
    return raw
end

function colspan_sum(tr::HTMLElement{:tr})
    # determine cols by resolve first colspan.
    s = 0
    for el in tr.children # el ∈ {tr, th}
        if "colspan" in keys(attrs(el))
            colspan = Base.parse(Int, attrs(el)["colspan"])
        else
            colspan = 1
        end
        s += colspan
    end
    s
end

function table_index(table::HTMLElement{:table})
    tr_list = eachmatch(sel"tr", table)
    cols = colspan_sum(tr_list[1])
    rows = length(tr_list)
    idx = 1
    index_mat = zeros(Int, rows, cols) # Use 0 to denote "unassigned"
    content_list = []
    for (i, tr) in enumerate(tr_list)
        j = 1
        for el in tr.children # el ∈ tr, td
            while index_mat[i, j] != 0
                j += 1
            end
            colspan = Base.parse(Int, get(attrs(el), "colspan", "1"))
            rowspan = Base.parse(Int, get(attrs(el), "rowspan", "1"))
            for ii in i:(i+rowspan-1)
                for jj in j:(j+colspan-1)
                    index_mat[ii, jj] = idx
                end
            end
            idx += 1
            push!(content_list, el)
        end
    end
    return content_list, index_mat
end

struct HTMLTableMatrix <: AbstractArray{HTMLElement, 2}
    content_list::Vector{HTMLElement}
    index_mat::Array{Int, 2}
end

Base.size(A::HTMLTableMatrix) = size(A.index_mat)
Base.getindex(A::HTMLTableMatrix, i::Int, j::Int) = A.content_list[A.index_mat[i,j]]

function get_unit_list(;cache=true)
    if cache
        data = get_from_cache()
    else
        data = get_from_wiki()
    end

    data_parsed = parsehtml(data)
    table_selected = eachmatch(sel"table.wikitable", data_parsed.root)[1]
    content_list, index_mat = table_index(table_selected)
    html_table_matrix = HTMLTableMatrix(content_list, index_mat);
    mat_extracted = extract_td.(html_table_matrix)[3:end, :] # remove 2 row head
    mat_parsed = Matrix{Any}(undef, size(mat_extracted)...)

    for i in 1:size(mat_extracted, 1)
        for (f_name, f_type, j) in zip(fieldnames(Unit), fieldtypes(Unit), 1:size(mat_extracted, 2))
            #global i_ = i
            #global j_ = j
            raw = mat_extracted[i, j]
            #global raw_ = raw
            parsed = unit_parse(f_type, raw)
            mat_parsed[i,j] = parsed
        end
    end

    unit_list = [Unit(mat_parsed[i,:]...) for i in 1:size(mat_parsed, 1)]
end