text:
using Test, Random
import Flux

@testset "structure" begin

  @testset "Join" begin
    @testset "model_creation" begin
      @test_nowarn Chain(Join(Dense(10,6),Dense(10,5)), Dense(11,1))
      #@test_throws DimensionMismatch Chain(Dense(10, 5, σ),Dense(2, 1))(randn(10))
    end
    @testset "model_training" begin
      model = Chain(
        Join(
          Chain(
            Dense(1, 5),
            Dense(5, 1)
          ),
          Dense(1, 2),
          Dense(1, 1),
        ),
        Dense(4, 1)
      )
      opt = Descent(0.1)
      loss(x, y) = Flux.mse(model(x), y)
      ps = Flux.params(model)
      x = map( x -> ( rand(1),rand(1),rand(1) ), 1:10 )
      y = map( x -> ( rand(1) ), 1:10)
      dat = zip(x, y)
      before = loss(x[1], y[1])
      map(x -> Flux.train!(loss, ps, dat, opt), 1:3)
      after = loss(x[1], y[1])
      @test after < before # loss should decrease with training
    end
  end

  @testset "Split" begin
    @testset "model_creation" begin
      @test_nowarn Chain(Dense(1, 1),Split(Dense(1, 1),Dense(1, 1),Dense(1, 1)))
      #@test_throws DimensionMismatch Chain(Dense(10, 5, σ),Dense(2, 1))(randn(10))
    end
    @testset "model_training" begin
      model = Chain(
        Dense(1, 1),
        Split(
          Dense(1, 1),
          Dense(1, 1),
          Dense(1, 1)
        )
      )
      opt = Descent(0.1)
      using Statistics
      function loss(x, y)
        sqrt(mean([Flux.mse(model(x)[i], y[i]) for i in 1:length(y)].^2.)) # rms over each mse
      end
      ps = Flux.params(model)
      x = map( x -> ( rand(1) ), 1:10)
      y = map( x -> ( rand(1),rand(1),rand(1) ), 1:10 )
      dat = zip(x, y)
      before = loss(x[1], y[1])
      map(x -> Flux.train!(loss, ps, dat, opt), 1:3)
      after = loss(x[1], y[1])
      @test after < before # loss should decrease with training
    end
  end

  @testset "Parallel" begin
    @testset "model_creation" begin
      @test_nowarn Chain(Dense(1, 1),Parallel(Dense(1, 1),Dense(1, 3),Chain(Dense(1, 5),Dense(5, 2),)),Dense(6, 1))
      #@test_throws DimensionMismatch Chain(Dense(10, 5, σ),Dense(2, 1))(randn(10))
    end
    @testset "model_training" begin
      model = Chain(
        Dense(1, 1),
        Parallel(
          Dense(1, 1),
          Dense(1, 3),
          Chain(
            Dense(1, 5),
            Dense(5, 2),
          )
        ),
        Dense(6, 1)
      )
      opt = Descent(0.1)
      loss(x, y) = Flux.mse(model(x), y)
      ps = Flux.params(model)
      x = map( x -> ( rand(1) ), 1:10)
      y = map( x -> ( rand(1) ), 1:10 )
      dat = zip(x, y)
      before = loss(x[1], y[1])
      map(x -> Flux.train!(loss, ps, dat, opt), 1:3)
      after = loss(x[1], y[1])
      @test after < before # loss should decrease with training
    end
  end

  @testset "Nop" begin
    @testset "model_creation" begin
      @test_nowarn Chain(Dense(1, 10),Dense(10, 1),Flux.Split(Nop,Nop))
    end
    @testset "nop_test_in_model" begin
      model = Chain(Dense(1, 10),Dense(10, 1),Flux.Split(Nop,Nop))
      result = model(rand(1))
      @test result[1] == result[2]
    end
  end

end


text:
export axisar

"""
    axisar(axis, angle)

Construct a rotation matrix that rotates vectors by a specified `angle` about a specified `axis`.

### Arguments ###

- `axis`: Rotation axis
- `angle`: Rotation angle in radians

### Output ###

Rotation matrix corresponding to `axis` and `angle`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/axisar_c.html)
"""
function axisar(axis, angle)
    @checkdims 3 axis
    r = Matrix{Float64}(undef, 3,3)
    ccall((:axisar_c, libcspice), Cvoid, (Ref{SpiceDouble}, SpiceDouble, Ref{SpiceDouble}), axis, angle, r)
    permutedims(r)
end

ext: jl


text:
## @lc code=start
using LeetCode

function maximum_difference(nums::Vector{Int})::Int
    mi, max_diff = nums[1], 0
    for num in @view nums[2:end]
        if num < mi
            mi = num
        elseif num - mi > max_diff
            max_diff = num - mi
        end
    end
    max_diff > 0 ? max_diff : -1
end

## @lc code=end
## @lc test=start
@testset "2016.maximum-difference-between-increasing-elements.jl" begin
    @test maximum_difference([7,1,5,4]) == 4
    @test maximum_difference([9,4,3,2]) == -1
    @test maximum_difference([1,5,2,10]) == 9
    @test maximum_difference([2,2,1]) == -1
end
## @lc test=end

## @lc info=start
# link: [solution to 2016]()
## @lc info=end


text:
import DelimitedFiles
import DifferentiableBackwardEuler
import NonlinearEquations
import PyPlot
import SparseArrays
import Zygote

@NonlinearEquations.equations exclude=(t, N, t2index, Bi) function icewedge(h, fvals, h_t, Bi, t, N, t2index)
	NonlinearEquations.setnumequations(N - 2)
	delta_r = 1 / (N - 1)
	time_index = t2index(t)
	for i = 2:N - 1
		r_i = (i - 1) * delta_r
		r_i_p_half = r_i + 0.5 * delta_r
		r_i_m_half = r_i - 0.5 * delta_r
		j = i - 1#shift the indexing, because our h starts at h[2] from the paper and we omit h[1] and h[N]
		if i == 2
			#the equation for i = 2 has h[i - 1] = h[1] = h[2] (equation 16), so h[i] - h[i - 1] = 0, so that term goes away
			NonlinearEquations.addterm(i - 1, (r_i_p_half * 0.5 * (h[j] + h[j + 1]) * (h[j + 1] - h[j]) / delta_r) / (r_i * delta_r) + fvals[time_index])
		elseif i == N - 1
			#the equation for i = N - 1 has h[i + 1] = h[N] = (h[N - 1] + delta_r * Bi * h_tr(t)) / (1 + delta_r * Bi) -- from equation 18
			NonlinearEquations.addterm(i - 1, ((r_i_p_half * 0.5 * (h[j] + ((h[j] + delta_r * Bi * h_tr[time_index]) / (1 + delta_r * Bi))) * (((h[j] + delta_r * Bi * h_tr[time_index]) / (1 + delta_r * Bi)) - h[j]) / delta_r) - (r_i_m_half * 0.5 * (h[j] + h[j - 1]) * (h[j] - h[j - 1]) / delta_r)) / (r_i * delta_r) + fvals[time_index])
		else
			NonlinearEquations.addterm(i - 1, ((r_i_p_half * 0.5 * (h[j] + h[j + 1]) * (h[j + 1] - h[j]) / delta_r) - (r_i_m_half * 0.5 * (h[j] + h[j - 1]) * (h[j] - h[j - 1]) / delta_r)) / (r_i * delta_r) + fvals[time_index])
		end
	end
end

N = 1001
@assert mod(N - 1, 1000) == 0
R = 10.0
D = 0.5
Kr = 0.5
kappa = 0.1
numtimesteps = 1000
@assert mod(numtimesteps, 1000) == 0
Sy = 0.3
h_tr1 = map(x->1.0-0.7/div(3 * numtimesteps, 10)*x, 1:(div(3 * numtimesteps, 10)))
h_tr2 = map(x->0.3*(1+sin(0.01*x / div(numtimesteps, 1000))),0:(div(7 * numtimesteps, 10) - 1))
h_tr = vcat(h_tr1, h_tr2)
ts = collect(range(0, 100; length=numtimesteps + 1)) * Kr * D / (Sy * R ^ 2)
fvals = zeros(length(ts))
fvals[500 * div(numtimesteps, 1000):600 * div(numtimesteps, 1000)] .= 2e-3
fvals .*= R ^ 2 / (Kr * D ^ 2)
function t2index(t)
	if t > maximum(ts) || t < 0
		error("crazy time: $t")
	else
		return round(Int, t / maximum(ts) * (length(ts) - 1))
	end
end
Bi = kappa * R / Kr

function unpack(p)
	return p
end

function f(u, p, t)
	fvals = unpack(p)
	return icewedge_residuals(u, fvals, h_tr, Bi, t, N, t2index)
end

function f_u(u, p, t)
	fvals = unpack(p)
	return icewedge_h(u, fvals, h_tr, Bi, t, N, t2index)
end

function f_p(u, p, t)
	fvals = unpack(p)
	return icewedge_fvals(u, fvals, h_tr, Bi, t, N, t2index)
end

f_t(u, p, t) = zeros(length(u))

p = fvals
h0 = ones(N - 2)
print("forward run time:")
@time h = DifferentiableBackwardEuler.steps(h0, f, f_u, f_p, f_t, p, ts; ftol=1e-12, method=:newton) * D#multiply by D to convert from h* to h
h = h[div(N - 1, 1000):div(N - 1, 1000):end, 1:div(numtimesteps, 1000):end]
h_vitaly = D .- DelimitedFiles.readdlm("vitaly_solution2.dlm")

function plot(h, i; savefig=true, filename="result.png")
	fig, ax = PyPlot.subplots()
	rs = collect(range(0, R; length=1001))
	ax.plot(rs[2:end - 1], h[:, i], "b", alpha=0.5, label="dan")
	ax.plot(rs[2:end - 1], h_vitaly[2:1000, i], "g", alpha=0.5, label="vitaly")
	ax.plot([rs[end]], [h_tr[i] * D], "k.", label="h_t")
	ax.set_xlabel("r")
	ax.set_xlabel("h")
	ax.set_ylim(-0.1, 0.6)
	ax.legend(loc=3)
	if savefig
		fig.savefig(filename)
	else
		display(fig)
		println()
	end
	PyPlot.close(fig)
end
for (i, j) = enumerate(1:10:size(h, 2) - 1)
	plot(h, j; filename="figs/frame_$(lpad(i, 4, '0')).png")
end
run(`ffmpeg -y -r 5 -f image2 -s 1920x1080 -i figs/frame_%04d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p movie.mp4`)

#=
g(p) = DifferentiableBackwardEuler.steps(h0, f, f_u, f_p, f_t, p, ts; ftol=1e-12, method=:newton)[1, end] * D
print("gradient time: ")
@time grad_zygote = Zygote.gradient(g, p)[1]

function checkgradientquickly(f, x0, gradf, n; delta::Float64=1e-8)
	indicestocheck = sort(collect(1:length(x0)), by=i->abs(gradf[i]), rev=true)[1:n]
	indicestocheck = [indicestocheck; rand(1:length(x0), n)]
	f0 = f(x0)
	for i in indicestocheck
		x = copy(x0)
		x[i] += delta
		fval = f(x)
		grad_f_i = (fval - f0) / delta
		@show i
		@show grad_f_i
		@show gradf[i]
	end
end
checkgradientquickly(g, p, grad_zygote, 2)
function checkresiduals()
	fig, ax = PyPlot.subplots()
	residuals = zeros(length(ts) - 1)
	residuals_vitaly = zeros(length(ts) - 1)
	for i = 2:1001
		deltah = ts[i] - ts[i - 1]
		residuals[i - 1] = sum((h[:, i] - h[:, i - 1] - deltah * f(h[:, i], p, ts[i])) .^ 2)
		residuals_vitaly[i - 1] = sum((h_vitaly[2:end - 1, i] - h_vitaly[2:end - 1, i - 1] - deltah * f(h_vitaly[2:end - 1, i], p, ts[i])) .^ 2)
	end
	i_worst = findmax(abs.(residuals .- residuals_vitaly))[2]
	#@show ts[i_worst], fvals[i_worst]
	#@show ts[i_worst - 1], fvals[i_worst - 1]
	ax.plot(residuals, alpha=0.5, label="dan")
	ax.plot(residuals_vitaly, alpha=0.5, label="vitaly")
	ax.legend()
	ax.set_xlabel("iteration number")
	ax.set_ylabel("sum of squared residuals")
	display(fig)
	println()
	PyPlot.close(fig)
end
if N == 1001#this only works if N = 1001
	checkresiduals()
end
=#


text:
# div σ = f
# div u = g
# σ = -p I + 2η ε(u)

using SymPy

x, y = @vars x y
u = 2*pi*sin(pi*x)*sin(pi*x)*cos(pi*y)*sin(pi*y)
v = -2*pi*sin(pi*x)*sin(pi*y)*cos(pi*x)*sin(pi*y)
p = sin(pi*x)*sin(pi*y)
nu = 0.5

f1 = diff(p, x) - nu * (diff(diff(u, x), x) + diff(diff(u, y), y) )
f2 = diff(p, y) - nu * (diff(diff(v, x), x) + diff(diff(v, y), y) )

println(replace(replace(sympy.julia_code(f1), ".*"=>"*"), ".^"=>"^"))
println(replace(replace(sympy.julia_code(f2), ".*"=>"*"), ".^"=>"^"))


text:
# using Distributed: clear!
# Functions to generate ϵ matrix.

################################################################################
# Sets up planetary-hierarchy index matrix using an initial condition file
# of the form "test.txt" (included in "test/" directory), or in a 2d array of
# the form file = [x ; "x,y,z,..."] where x,y,z are Int64
################################################################################
function hierarchy(ic_vec::Array{Int64,1})

    # Separates the initial array into variables
    nbody = ic_vec[1]
    bins = ic_vec[2:end]
    #bins = replace(bins,"," => " ")
    #bins = readdlm(IOBuffer(bins),Int64)

    # checks that the hierarchy can be created
    bins[end] == 1 ? level = length(bins) : level = length(bins) - 1
    @assert bins[1] <= nbody/2 && bins[1] <= 2^level "Invalid hierarchy."
    @assert sum(bins[1:level]) == nbody - 1 "Invalid hierarchy. Total # of binaries should = N-1"

    # Creates an empty array of the proper size
    h = zeros(Float64,nbody,nbody)
    # Starts filling the array and returns it
    bottom_level(nbody,bins,h,1)
    h[end,:] .= -1
    return h
end

################################################################################
# Sets up the first level of the hierchary.
################################################################################
function bottom_level(nbody::Int64,bins::Array{Int64,1},h::Array{Float64},iter::Int64)

    # Fills the very first level of the hierarchy and iterates related variables
    h[1,1] = -1.0
    h[1,2] = 1.0
    if bins[1] > 1
        bins[end] == 1 ? j = 1 : j::Int64 = nbody/2 - 1
        for i=1:bins[1]-1
            if (j+3) <= nbody
                h[i+1,j+2] = -1
                h[i+1,j+3] = 1
                j += 2
            end
        end
    end
    binsp = bins[1]
    row = binsp[1] + 1
    bodies = bins[1] * 2

    # checks whether the desired hierarchy is symmetric and calls appropriate
    # filling functions
    if bins[end] > 1
        bins = bins[1:end-1]
        symmetric_nest(nbody,bodies,bins,binsp,row,h,iter+1)
    elseif 2*binsp == nbody
        symmetric(nbody,bodies,bins,binsp,row,h,iter+1)
    else
        nlevel(nbody,bodies,bins,binsp,row,h,iter+1)
    end
end

################################################################################
# Fills subsequent levels of the hierarchy, recursively.
#
# *** Only works for hierarchies with a max binaries per level of 2 (unless
# symmetric). ***
################################################################################
function nlevel(nbody::Int64,bodies::Int64,bins::Array{Int64,1},binsp::Int64,row::Int64,h::Array{Float64},iter::Int64)
#print("nbody: ", nbody, "\n")
#print("bodies: ", bodies, "\n")
#print("row: ", row, "\n")
# Series of checks to know which level to fill and which bins number to use
if iter <= length(bins)
    if nbody != bodies
        if bins[iter] == binsp
            bodies = bodies + bins[iter]
        elseif bins[iter] > binsp
            bodies = bodies + 2*(bins[iter]) - 1
        end
    elseif nbody == bodies
        if row == nbody-1
            if (bodies%4) > 2
                h[row,1:row-2] .= -1
                h[row,row-1:bodies] .= 1
                return h
            elseif (bodies%4) <= 2
                h[row,1:row-1] .= -1
                h[row,row:bodies] .= 1
                return h
            end
        end
    end

    # Looks at the previous and current bins and calls nlevel again
    if nbody >= bodies
        if binsp == 1
            if bins[iter] == 1
                h[row,1:bodies-binsp] .= -1
                h[row,bodies-binsp + 1:bodies] .= 1
                row = row + binsp
                binsp = bins[iter]
                nlevel(nbody,bodies,bins,binsp,row,h,iter+1)
            elseif bins[iter] == 2
                h[row,1:bodies-(2*binsp)-1] .= -1
                h[row,bodies-(2*binsp)] = 1
                h[row+1,bodies-(2*binsp)+1] = -1
                h[row+1,bodies] = 1
                row = row + 2
                binsp = bins[iter]
                nlevel(nbody,bodies,bins,binsp,row,h,iter+1)
            end
        elseif binsp == 2
            if bins[iter] == 1
                h[row,1:row-1] .= -1
                h[row,row:bodies] .= 1
                row = row + 1
                binsp = bins[iter]
                nlevel(nbody,bodies,bins,binsp,row,h,iter+1)
            elseif bins[iter] == 2
                h[row,1:row-1] .= -1
                h[row,row:bodies-2*(binsp-1)] .= 1
                h[row+1,bodies-2*(binsp-1)+1] = -1
                h[row+1,bodies] = 1
                row = row + 2
                binsp = bins[iter]
                nlevel(nbody,bodies,bins,binsp,row,h,iter+1)
            end
        elseif binsp == 3
            if bins[iter] == 2
                h[row,1:Int64(bodies/2)-1] .= -1
                h[row,Int64(bodies/2):2*Int64(bodies/3)] .= 1
                h[row+1,2*Int64(bodies/3)+1:bodies] .= -1
                h[row+1,bodies+1] = 1
                row = row + 2
                binsp = bins[iter]
                bodies = bodies + 1
                nlevel(nbody,bodies,bins,binsp,row,h,iter+1)
            end
        end
    end
end
end

################################################################################
# Called if the hierarchy is symmetric (or if the hierarchy has all of the
# bodies on the bottom level). Fills the remaining rows.
################################################################################
function symmetric(nbody::Int64,bodies::Int64,bins::Array{Int64,1},binsp::Int64,row::Int64,h::Array{Float64},iter::Int64)
    #global j = 0
    j = 0
    while row < nbody-1
        h[row,1+j:2+j] .= -1
        h[row,j+3:j+4] .= 1
        j = j + 4
        row = row + 1
    end
    h[row,1:binsp] .= -1
    h[row,binsp+1:nbody] .= 1
    #clear!(:j)
end

function symmetric_nest(nbody,bodies,bins,binsp,row,h,iter)
    hlf::Int64 = nbody/2
    j = 0
    while row < nbody-1
        println(j)
        h[row,1:2+j] .= -1
        h[row,3+j] =  1
        h[row+1,hlf+1:hlf+2+j] .= -1
        h[row+1,hlf+3+j]  = 1
        j += 1
        row += 2
    end
    h[row,1:hlf] .= -1
    h[row,hlf+1:nbody] .= 1
end


text:
# Constructs a Huffman coding tree for efficiently generating discrete random
# variables by relative weights.

# Maximum granularity is 1/float64(typemax(Uint64)) = 5.4-20, but also
# supports non-base2 divisions via rejection sampling.

import Base.isless, Base.show, Base.getindex, Base.rand

abstract HuffmanNode{T} <: AbstractCategoricalSampler

immutable HuffmanLeaf{T} <: HuffmanNode{T}
    value::T
    weight::Uint64
end

immutable HuffmanBranch{T} <: HuffmanNode{T}
    left::HuffmanNode{T}
    right::HuffmanNode{T}
    weight::Uint64
end
HuffmanBranch{T}(ha::HuffmanNode{T},hb::HuffmanNode{T}) = HuffmanBranch(ha, hb, ha.weight + hb.weight)

isless{T}(ha::HuffmanNode{T}, hb::HuffmanNode{T}) = isless(ha.weight,hb.weight)

show{T}(io::IO, t::HuffmanNode{T}) = show(io,typeof(t))

function getindex{T}(h::HuffmanBranch{T},u::Uint64)
    while isa(h,HuffmanBranch{T})
        if u < h.left.weight
            h = h.left
        else
            u -= h.left.weight
            h = h.right
        end
    end
    h.value
end

# build the huffman tree
# could be slightly more efficient using a Deque.
function huffman{T}(values::AbstractVector{T},weights::AbstractVector{Uint64})
    leafs = [HuffmanLeaf{T}(values[i],weights[i]) for i = 1:length(weights)]
    sort!(leafs; rev=true)

    branches = Array(HuffmanBranch{T},0)

    while !isempty(leafs) || length(branches) > 1
        left = isempty(branches) || (!isempty(leafs) && first(leafs) < first(branches)) ? pop!(leafs) : pop!(branches)
        right = isempty(branches) || (!isempty(leafs) && first(leafs) < first(branches)) ? pop!(leafs) : pop!(branches)
        unshift!(branches,HuffmanBranch(left,right))
    end

    pop!(branches)
end

function rand{T}(h::HuffmanNode{T})
    w = h.weight
    # generate uniform Uint64 objects on the range 0:(w-1)
    # unfortunately we can't use Range objects, as they don't have sufficient length
    u = rand(Uint64)
    if (w & (w-1)) == 0
        # power of 2
        u = u & (w-1)
    else
        m = typemax(Uint64)
        lim = m - (rem(m,w)+1)
        while u > lim
            u = rand(Uint64)
        end
        u = rem(u,w)
    end
    h[u]
end


text:
using Documenter, TextWrap

makedocs(
    modules  = [TextWrap],
    format   = :html,
    sitename = "TextWrap.jl",
    pages    = Any[
        "Home" => "index.md",
       ]
    )

deploydocs(
    repo   = "github.com/carlobaldassi/TextWrap.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    julia  = "0.6"
)


text:
using Assimp
using Test

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


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GenericBattery <: Storage
        name::String
        available::Bool
        bus::Bus
        prime_mover::PrimeMovers
        initial_energy::Float64
        state_of_charge_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        rating::Float64
        active_power::Float64
        input_active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        output_active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
        reactive_power::Float64
        reactive_power_limits::Union{Nothing, Min_Max}
        base_power::Float64
        operation_cost::Union{Nothing, OperationalCost}
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

Data structure for a generic battery

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `prime_mover::PrimeMovers`: Prime mover technology according to EIA 923
- `initial_energy::Float64`: State of Charge of the Battery p.u.-hr, validation range: `(0, nothing)`, action if invalid: `error`
- `state_of_charge_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Maximum and Minimum storage capacity in p.u.-hr, validation range: `(0, nothing)`, action if invalid: `error`
- `rating::Float64`
- `active_power::Float64`
- `input_active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: `(0, nothing)`, action if invalid: `error`
- `output_active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: `(0, nothing)`, action if invalid: `error`
- `efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}`, validation range: `(0, 1)`, action if invalid: `warn`
- `reactive_power::Float64`, validation range: `reactive_power_limits`, action if invalid: `warn`
- `reactive_power_limits::Union{Nothing, Min_Max}`
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `operation_cost::Union{Nothing, OperationalCost}`
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct GenericBattery <: Storage
    name::String
    available::Bool
    bus::Bus
    "Prime mover technology according to EIA 923"
    prime_mover::PrimeMovers
    "State of Charge of the Battery p.u.-hr"
    initial_energy::Float64
    "Maximum and Minimum storage capacity in p.u.-hr"
    state_of_charge_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    rating::Float64
    active_power::Float64
    input_active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    output_active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
    reactive_power::Float64
    reactive_power_limits::Union{Nothing, Min_Max}
    "Base power of the unit in MVA"
    base_power::Float64
    operation_cost::Union{Nothing, OperationalCost}
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GenericBattery(name, available, bus, prime_mover, initial_energy, state_of_charge_limits, rating, active_power, input_active_power_limits, output_active_power_limits, efficiency, reactive_power, reactive_power_limits, base_power, operation_cost=nothing, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    GenericBattery(name, available, bus, prime_mover, initial_energy, state_of_charge_limits, rating, active_power, input_active_power_limits, output_active_power_limits, efficiency, reactive_power, reactive_power_limits, base_power, operation_cost, services, dynamic_injector, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function GenericBattery(; name, available, bus, prime_mover, initial_energy, state_of_charge_limits, rating, active_power, input_active_power_limits, output_active_power_limits, efficiency, reactive_power, reactive_power_limits, base_power, operation_cost=nothing, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    GenericBattery(name, available, bus, prime_mover, initial_energy, state_of_charge_limits, rating, active_power, input_active_power_limits, output_active_power_limits, efficiency, reactive_power, reactive_power_limits, base_power, operation_cost, services, dynamic_injector, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function GenericBattery(::Nothing)
    GenericBattery(;
        name="init",
        available=false,
        bus=Bus(nothing),
        prime_mover=PrimeMovers.BA,
        initial_energy=0.0,
        state_of_charge_limits=(min=0.0, max=0.0),
        rating=0.0,
        active_power=0.0,
        input_active_power_limits=(min=0.0, max=0.0),
        output_active_power_limits=(min=0.0, max=0.0),
        efficiency=(in=0.0, out=0.0),
        reactive_power=0.0,
        reactive_power_limits=(min=0.0, max=0.0),
        base_power=0.0,
        operation_cost=nothing,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`GenericBattery`](@ref) `name`."""
get_name(value::GenericBattery) = value.name
"""Get [`GenericBattery`](@ref) `available`."""
get_available(value::GenericBattery) = value.available
"""Get [`GenericBattery`](@ref) `bus`."""
get_bus(value::GenericBattery) = value.bus
"""Get [`GenericBattery`](@ref) `prime_mover`."""
get_prime_mover(value::GenericBattery) = value.prime_mover
"""Get [`GenericBattery`](@ref) `initial_energy`."""
get_initial_energy(value::GenericBattery) = get_value(value, value.initial_energy)
"""Get [`GenericBattery`](@ref) `state_of_charge_limits`."""
get_state_of_charge_limits(value::GenericBattery) = get_value(value, value.state_of_charge_limits)
"""Get [`GenericBattery`](@ref) `rating`."""
get_rating(value::GenericBattery) = get_value(value, value.rating)
"""Get [`GenericBattery`](@ref) `active_power`."""
get_active_power(value::GenericBattery) = get_value(value, value.active_power)
"""Get [`GenericBattery`](@ref) `input_active_power_limits`."""
get_input_active_power_limits(value::GenericBattery) = get_value(value, value.input_active_power_limits)
"""Get [`GenericBattery`](@ref) `output_active_power_limits`."""
get_output_active_power_limits(value::GenericBattery) = get_value(value, value.output_active_power_limits)
"""Get [`GenericBattery`](@ref) `efficiency`."""
get_efficiency(value::GenericBattery) = value.efficiency
"""Get [`GenericBattery`](@ref) `reactive_power`."""
get_reactive_power(value::GenericBattery) = get_value(value, value.reactive_power)
"""Get [`GenericBattery`](@ref) `reactive_power_limits`."""
get_reactive_power_limits(value::GenericBattery) = get_value(value, value.reactive_power_limits)
"""Get [`GenericBattery`](@ref) `base_power`."""
get_base_power(value::GenericBattery) = value.base_power
"""Get [`GenericBattery`](@ref) `operation_cost`."""
get_operation_cost(value::GenericBattery) = value.operation_cost
"""Get [`GenericBattery`](@ref) `services`."""
get_services(value::GenericBattery) = value.services
"""Get [`GenericBattery`](@ref) `dynamic_injector`."""
get_dynamic_injector(value::GenericBattery) = value.dynamic_injector
"""Get [`GenericBattery`](@ref) `ext`."""
get_ext(value::GenericBattery) = value.ext
"""Get [`GenericBattery`](@ref) `time_series_container`."""
get_time_series_container(value::GenericBattery) = value.time_series_container
"""Get [`GenericBattery`](@ref) `internal`."""
get_internal(value::GenericBattery) = value.internal

"""Set [`GenericBattery`](@ref) `name`."""
set_name!(value::GenericBattery, val) = value.name = val
"""Set [`GenericBattery`](@ref) `available`."""
set_available!(value::GenericBattery, val) = value.available = val
"""Set [`GenericBattery`](@ref) `bus`."""
set_bus!(value::GenericBattery, val) = value.bus = val
"""Set [`GenericBattery`](@ref) `prime_mover`."""
set_prime_mover!(value::GenericBattery, val) = value.prime_mover = val
"""Set [`GenericBattery`](@ref) `initial_energy`."""
set_initial_energy!(value::GenericBattery, val) = value.initial_energy = val
"""Set [`GenericBattery`](@ref) `state_of_charge_limits`."""
set_state_of_charge_limits!(value::GenericBattery, val) = value.state_of_charge_limits = val
"""Set [`GenericBattery`](@ref) `rating`."""
set_rating!(value::GenericBattery, val) = value.rating = val
"""Set [`GenericBattery`](@ref) `active_power`."""
set_active_power!(value::GenericBattery, val) = value.active_power = val
"""Set [`GenericBattery`](@ref) `input_active_power_limits`."""
set_input_active_power_limits!(value::GenericBattery, val) = value.input_active_power_limits = val
"""Set [`GenericBattery`](@ref) `output_active_power_limits`."""
set_output_active_power_limits!(value::GenericBattery, val) = value.output_active_power_limits = val
"""Set [`GenericBattery`](@ref) `efficiency`."""
set_efficiency!(value::GenericBattery, val) = value.efficiency = val
"""Set [`GenericBattery`](@ref) `reactive_power`."""
set_reactive_power!(value::GenericBattery, val) = value.reactive_power = val
"""Set [`GenericBattery`](@ref) `reactive_power_limits`."""
set_reactive_power_limits!(value::GenericBattery, val) = value.reactive_power_limits = val
"""Set [`GenericBattery`](@ref) `base_power`."""
set_base_power!(value::GenericBattery, val) = value.base_power = val
"""Set [`GenericBattery`](@ref) `operation_cost`."""
set_operation_cost!(value::GenericBattery, val) = value.operation_cost = val
"""Set [`GenericBattery`](@ref) `services`."""
set_services!(value::GenericBattery, val) = value.services = val
"""Set [`GenericBattery`](@ref) `ext`."""
set_ext!(value::GenericBattery, val) = value.ext = val
"""Set [`GenericBattery`](@ref) `time_series_container`."""
set_time_series_container!(value::GenericBattery, val) = value.time_series_container = val


text:
@testset "conic linear problems" begin
    @testset "conic_linear1_test" begin
    #=
    primal
        min -3x - 2y - 4z
    s.t.
        x +  y +  z == 3 :w_4
        y +  z == 2      :w_5
        x>=0 y>=0 z>=0   :w_1, w_2, w_3

    dual
        max 3w_4 + 2w_5
    s.t.
        w_1, w_2, w_3 \in Nonnegatives
        w_1 + w_4 == -3.0
        w_2 + w_4 + w_5 == -2.0
        w_3 + w_4 + w_5 == -4.0
    =#
        primal_model = conic_linear1_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
        list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set([
            (SAF{Float64}, MOI.EqualTo{Float64})
            (VVF, MOI.Nonnegatives)
        ])
        @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 3
        @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.Nonnegatives}()) == 1
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == 0.0
        @test MOI.coefficient.(obj.terms) == [3.0; 2.0]

        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
        eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == -3.0
        eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
        eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0; 1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == -2.0
        eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
        eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
        @test MOI.coefficient.(eq_con3_fun.terms) == [1.0; 1.0; 1.0]
        @test MOI.constant.(eq_con3_fun) == 0.0
        @test MOI.constant(eq_con3_set) == -4.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{VAF{Float64}, MOI.Zeros}(2)] == VI.(4:5)
        @test primal_con_dual_var[CI{VVF, MOI.Nonnegatives}(1)] == VI.(1:3)

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
        @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
        @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
    end

    @testset "conic_linear3_test" begin
    #=
    primal
        min  3x + 2y - 4z + 0s
    s.t.
        x           -  s == -4    :w_4
            y            == -3    :w_5
        x      +  z      == 12    :w_6
        y <= 0 :w_3
        z >= 0 :w_2
        s zero :w_1

    dual
        max -4w_4 - 3w_5 + 12w_6
    s.t
        w_4 + w_6 == 3
        w_3 + w_5 == 2
        w_2 + w_6 == -4
        w_1 - w_4 == 0
        w_2 >= 0
        w_3 <= 0
        w_1 in Reals
    =#
        primal_model = conic_linear3_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 6
        list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set([
            (SAF{Float64}, MOI.EqualTo{Float64})
            (VVF, MOI.Nonnegatives)
            (VVF, MOI.Nonpositives)
        ])
        @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 4
        @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.Nonnegatives}()) == 1
        @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.Nonpositives}()) == 1
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == 0.0
        @test MOI.coefficient.(obj.terms) == [-4.0; -3.0; 12.0]

        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
        eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 3.0
        eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
        eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == 2.0
        eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(5))
        eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(5))
        @test MOI.coefficient.(eq_con3_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con3_fun) == 0.0
        @test MOI.constant(eq_con3_set) == -4.0
        eq_con6_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(6))
        eq_con6_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(6))
        @test MOI.coefficient.(eq_con6_fun.terms) == [1.0; -1.0]
        @test MOI.constant.(eq_con6_fun) == 0.0
        @test MOI.constant(eq_con6_set) == 0.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{VVF, MOI.Zeros}(4)] == [VI(1)]
        @test primal_con_dual_var[CI{VVF, MOI.Nonnegatives}(3)] == [VI(2)]
        @test primal_con_dual_var[CI{VVF, MOI.Nonpositives}(2)] == [VI(3)]
        @test primal_con_dual_var[CI{VAF{Float64}, MOI.Zeros}(1)] == VI.(4:6)

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
        @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
        @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(5)
        @test primal_var_dual_con[VI(4)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(6)
    end
end


text:
using Random
Random.seed!(1)

numbers = 10:25
N = 10^7

firstDigit(x) = Int(floor(x/10))
secondDigit(x) = x%10

numThirteen, numFirstIsOne, numSecondIsThree = 0, 0, 0

for _ in 1:N
    X = rand(numbers)
    global numThirteen += X == 13
    global numFirstIsOne += firstDigit(X) == 1
    global numSecondIsThree += secondDigit(X) == 3
end

probThirteen, probFirstIsOne, probSecondIsThree =
    (numThirteen,numFirstIsOne,numSecondIsThree)./N

println("P(13) = ", round(probThirteen, digits=4),
        "\nP(1_) = ",round(probFirstIsOne, digits=4),
        "\nP(_3) = ", round(probSecondIsThree, digits=4),
        "\nP(1_)*P(_3) = ",round(probFirstIsOne*probSecondIsThree, digits=4))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
mutable struct LinkedInts
    next :: Vector{Int}
    prev :: Vector{Int}

    free_ptr :: Int
    free_cap :: Int
    root     :: Int

    block :: Vector{Int}
    size  :: Vector{Int}
end

LinkedInts(capacity=128) =
    LinkedInts(Int[],Int[],
               0,0,0,
               Int[],
               Int[])

allocated(s::LinkedInts, id :: Int) = id > 0 && id <= length(s.block) && s.block[id] > 0
blocksize(s::LinkedInts, id :: Int) = s.size[id]
Base.length(s::LinkedInts) = length(s.next)
numblocks(s::LinkedInts) = length(s.block)

function Base.string(l::LinkedInts)
    r = String[]
    push!(r,"LinkedInts(\n")
    push!(r,@sprintf("  Number of blocks: %d\n", length(l.block)))
    push!(r,"  Blocks:\n")

    for i in 1:length(l.block)
        if l.block[i] > 0
            idxs = getindexes(l,i)
            push!(r,@sprintf("    #%d: %s\n",i,string(idxs)))
        end
    end
    p = l.free_ptr
    freelst = Int[]
    while p > 0
        push!(freelst,p)
        p = l.prev[p]
    end
    push!(r,@sprintf("  Free: %s\n",string(freelst)))

    push!(r,")")
    join(r)
end

"""
    ensurefree(s::LinkedInts, N :: Int)

Ensure that there are at least `N` elements free, and allocate as necessary.
"""
function ensurefree(s::LinkedInts, N :: Int)
    if s.free_cap < N
        num = N - s.free_cap

        cap = length(s.next)
        first = cap+1
        last  = cap+num

        append!(s.next,Int[i+1 for i in first:last])
        append!(s.prev,Int[i-1 for i in first:last])

        s.next[last] = 0
        s.prev[first] = s.free_ptr
        if s.prev[first] > 0
            s.next[s.prev[first]] = first
        end
        s.free_ptr = last
        s.free_cap += num

        num
    else
        0
    end
end

"""
    newblock(s::LinkedInts, N :: Int)

Add a new block in list `idx`
"""
function newblock(s::LinkedInts, N :: Int) :: Int
    assert(N>0)
    ensurefree(s,N)
    # remove from free list
    ptre = s.free_ptr
    ptrb = ptre
    for i = 1:N-1
        ptrb = s.prev[ptrb]
    end

    prev = s.prev[ptrb]

    if prev > 0
        s.next[prev] = 0
    end

    s.free_ptr = s.prev[ptrb]
    s.free_cap -= N

    # insert into list `idx`
    s.prev[ptrb] = s.root
    if s.root > 0
        s.next[s.root] = ptrb
    end
    s.root = ptre
    push!(s.block,ptrb)
    push!(s.size,N)

    id = length(s.block)

    id
end

"""
Move a block to the free list.
"""
function deleteblock(s::LinkedInts, id :: Int)
    if s.size[id] > 0
        ptrb = s.block[id]
        N = s.size[id]
        ptre = ptrb
        for i in 2:N
            ptre = s.next[ptre]
        end
        prev = s.prev[ptrb]
        next = s.next[ptre]

        # remove from list and clear the block id
        if s.root == ptre s.root = prev end
        if prev > 0 s.next[prev] = next end
        if next > 0 s.prev[next] = prev end

        s.size[id]  = 0
        s.block[id] = 0

        # add to free list
        if s.free_ptr > 0
            s.next[s.free_ptr] = ptrb
        end
        s.prev[ptrb] = s.free_ptr
        s.free_ptr = ptre
        s.next[ptre] = 0

        s.free_cap += N
    end
end

"""
    getindexes(s::LinkedInts, id :: Int)

"""
function getindexes(s::LinkedInts, id :: Int)
    N = s.size[id]
    r = Array{Int}(N)
    p = s.block[id]
    for i in 1:N
        r[i] = p
        p = s.next[p]
    end
    r
end

function getindexes(s::LinkedInts, id :: Int, target :: Array{Int,1}, offset :: Int)
    N = s.size[id]
    p = s.block[id]
    for i in 1:N
        target[i+offset-1] = p
        p = s.next[p]
    end
    N
end

function getoneindex(s::LinkedInts, id :: Int)
    N = s.size[id]
    if N < 1
        error("No values at id")
    end

    s.block[i]
end

"""
Get a list if the currently free elements.
"""
function getfreeindexes(s::LinkedInts)
    N = s.free_cap
    r = Array{Int}(N)
    ptr = s.free_ptr
    for i in 1:N
        r[N-i+1] = ptr
        ptr  = s.prev[ptr]
    end
    r
end

"""
Get a list if the currently free elements.
"""
function getusedindexes(s::LinkedInts)
    N = length(s.next) - s.free_cap
    r = Array{Int}(N)
    ptr = s.root
    for i in 1:N
        r[N-i+1] = ptr
        ptr  = s.prev[ptr]
    end
    r
end

"""
Check consistency of the internal structures.
"""
function checkconsistency(s::LinkedInts) :: Bool
    if length(s.prev) != length(s.next)
        return false
    end

    N = length(s.prev)

    all(i -> s.prev[i] == 0 || s.next[s.prev[i]] == i, 1:N) &&
    all(i -> s.next[i] == 0 || s.prev[s.next[i]] == i, 1:N)
end


text:
# Autogenerated wrapper script for Qt_jll for x86_64-w64-mingw32-cxx03
export libqt53danimation, libqt53dcore, libqt53dextras, libqt53dinput, libqt53dlogic, libqt53dquick, libqt53dquickanimation, libqt53dquickextras, libqt53dquickinput, libqt53dquickrender, libqt53dquickscene2d, libqt53drender, libqt5bodymovin, libqt5charts, libqt5concurrent, libqt5core, libqt5dbus, libqt5gui, libqt5multimedia, libqt5multimediaquick, libqt5multimediawidgets, libqt5network, libqt5opengl, libqt5printsupport, libqt5qml, libqt5qmlmodels, libqt5qmlworkerscript, libqt5quick, libqt5quick3d, libqt5quick3dassetimport, libqt5quick3drender, libqt5quick3druntimerender, libqt5quick3dutils, libqt5quickcontrols2, libqt5quickparticles, libqt5quickshapes, libqt5quicktemplates2, libqt5quicktest, libqt5quickwidgets, libqt5sql, libqt5svg, libqt5test, libqt5widgets, libqt5xml, libqt5xmlpatterns

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt")
JLLWrappers.@declare_library_product(libqt53danimation, "Qt53DAnimation.dll")
JLLWrappers.@declare_library_product(libqt53dcore, "Qt53DCore.dll")
JLLWrappers.@declare_library_product(libqt53dextras, "Qt53DExtras.dll")
JLLWrappers.@declare_library_product(libqt53dinput, "Qt53DInput.dll")
JLLWrappers.@declare_library_product(libqt53dlogic, "Qt53DLogic.dll")
JLLWrappers.@declare_library_product(libqt53dquick, "Qt53DQuick.dll")
JLLWrappers.@declare_library_product(libqt53dquickanimation, "Qt53DQuickAnimation.dll")
JLLWrappers.@declare_library_product(libqt53dquickextras, "Qt53DQuickExtras.dll")
JLLWrappers.@declare_library_product(libqt53dquickinput, "Qt53DQuickInput.dll")
JLLWrappers.@declare_library_product(libqt53dquickrender, "Qt53DQuickRender.dll")
JLLWrappers.@declare_library_product(libqt53dquickscene2d, "Qt53DQuickScene2D.dll")
JLLWrappers.@declare_library_product(libqt53drender, "Qt53DRender.dll")
JLLWrappers.@declare_library_product(libqt5bodymovin, "Qt5Bodymovin.dll")
JLLWrappers.@declare_library_product(libqt5charts, "Qt5Charts.dll")
JLLWrappers.@declare_library_product(libqt5concurrent, "Qt5Concurrent.dll")
JLLWrappers.@declare_library_product(libqt5core, "Qt5Core.dll")
JLLWrappers.@declare_library_product(libqt5dbus, "Qt5DBus.dll")
JLLWrappers.@declare_library_product(libqt5gui, "Qt5Gui.dll")
JLLWrappers.@declare_library_product(libqt5multimedia, "Qt5Multimedia.dll")
JLLWrappers.@declare_library_product(libqt5multimediaquick, "Qt5MultimediaQuick.dll")
JLLWrappers.@declare_library_product(libqt5multimediawidgets, "Qt5MultimediaWidgets.dll")
JLLWrappers.@declare_library_product(libqt5network, "Qt5Network.dll")
JLLWrappers.@declare_library_product(libqt5opengl, "Qt5OpenGL.dll")
JLLWrappers.@declare_library_product(libqt5printsupport, "Qt5PrintSupport.dll")
JLLWrappers.@declare_library_product(libqt5qml, "Qt5Qml.dll")
JLLWrappers.@declare_library_product(libqt5qmlmodels, "Qt5QmlModels.dll")
JLLWrappers.@declare_library_product(libqt5qmlworkerscript, "Qt5QmlWorkerScript.dll")
JLLWrappers.@declare_library_product(libqt5quick, "Qt5Quick.dll")
JLLWrappers.@declare_library_product(libqt5quick3d, "Qt5Quick3D.dll")
JLLWrappers.@declare_library_product(libqt5quick3dassetimport, "Qt5Quick3DAssetImport.dll")
JLLWrappers.@declare_library_product(libqt5quick3drender, "Qt5Quick3DRender.dll")
JLLWrappers.@declare_library_product(libqt5quick3druntimerender, "Qt5Quick3DRuntimeRender.dll")
JLLWrappers.@declare_library_product(libqt5quick3dutils, "Qt5Quick3DUtils.dll")
JLLWrappers.@declare_library_product(libqt5quickcontrols2, "Qt5QuickControls2.dll")
JLLWrappers.@declare_library_product(libqt5quickparticles, "Qt5QuickParticles.dll")
JLLWrappers.@declare_library_product(libqt5quickshapes, "Qt5QuickShapes.dll")
JLLWrappers.@declare_library_product(libqt5quicktemplates2, "Qt5QuickTemplates2.dll")
JLLWrappers.@declare_library_product(libqt5quicktest, "Qt5QuickTest.dll")
JLLWrappers.@declare_library_product(libqt5quickwidgets, "Qt5QuickWidgets.dll")
JLLWrappers.@declare_library_product(libqt5sql, "Qt5Sql.dll")
JLLWrappers.@declare_library_product(libqt5svg, "Qt5Svg.dll")
JLLWrappers.@declare_library_product(libqt5test, "Qt5Test.dll")
JLLWrappers.@declare_library_product(libqt5widgets, "Qt5Widgets.dll")
JLLWrappers.@declare_library_product(libqt5xml, "Qt5Xml.dll")
JLLWrappers.@declare_library_product(libqt5xmlpatterns, "Qt5XmlPatterns.dll")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt53danimation,
        "bin\\Qt53DAnimation.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dcore,
        "bin\\Qt53DCore.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dextras,
        "bin\\Qt53DExtras.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dinput,
        "bin\\Qt53DInput.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dlogic,
        "bin\\Qt53DLogic.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquick,
        "bin\\Qt53DQuick.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickanimation,
        "bin\\Qt53DQuickAnimation.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickextras,
        "bin\\Qt53DQuickExtras.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickinput,
        "bin\\Qt53DQuickInput.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickrender,
        "bin\\Qt53DQuickRender.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickscene2d,
        "bin\\Qt53DQuickScene2D.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53drender,
        "bin\\Qt53DRender.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5bodymovin,
        "bin\\Qt5Bodymovin.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5charts,
        "bin\\Qt5Charts.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5concurrent,
        "bin\\Qt5Concurrent.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5core,
        "bin\\Qt5Core.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5dbus,
        "bin\\Qt5DBus.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5gui,
        "bin\\Qt5Gui.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimedia,
        "bin\\Qt5Multimedia.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimediaquick,
        "bin\\Qt5MultimediaQuick.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimediawidgets,
        "bin\\Qt5MultimediaWidgets.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5network,
        "bin\\Qt5Network.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5opengl,
        "bin\\Qt5OpenGL.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5printsupport,
        "bin\\Qt5PrintSupport.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qml,
        "bin\\Qt5Qml.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlmodels,
        "bin\\Qt5QmlModels.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlworkerscript,
        "bin\\Qt5QmlWorkerScript.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick,
        "bin\\Qt5Quick.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3d,
        "bin\\Qt5Quick3D.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3dassetimport,
        "bin\\Qt5Quick3DAssetImport.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3drender,
        "bin\\Qt5Quick3DRender.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3druntimerender,
        "bin\\Qt5Quick3DRuntimeRender.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3dutils,
        "bin\\Qt5Quick3DUtils.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickcontrols2,
        "bin\\Qt5QuickControls2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickparticles,
        "bin\\Qt5QuickParticles.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickshapes,
        "bin\\Qt5QuickShapes.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktemplates2,
        "bin\\Qt5QuickTemplates2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktest,
        "bin\\Qt5QuickTest.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickwidgets,
        "bin\\Qt5QuickWidgets.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5sql,
        "bin\\Qt5Sql.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5svg,
        "bin\\Qt5Svg.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5test,
        "bin\\Qt5Test.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5widgets,
        "bin\\Qt5Widgets.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xml,
        "bin\\Qt5Xml.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xmlpatterns,
        "bin\\Qt5XmlPatterns.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PyPlot

data = readdlm("bounce_zeroorder.out")

rc("text", usetex=true)
rc("figure", figsize=[2.8, 2.2])
rc("lines", linewidth=0.5)
rc("lines", markersize=3)

plot(sqrt(data[:,2]).*cos(data[:,4]), sqrt(data[:,2]).*sin(data[:,4]),"k:")
xlabel("\$x\$")
ylabel("\$y\$")
title("\$y\$")
tight_layout()
#savefig("test.png")

# a = 3
# b = "bla"
#
# ccall((:__common_MOD_disp, "libdriftorbit.so"), Int32, (Ptr{Cstring},Ptr{Cfloat}),
#        &b,&a)

max_stars_repo_name: itpplasma/NEO-RT

max_issues_repo_name: itpplasma/NEO-RT

max_forks_repo_name: itpplasma/NEO-RT


text:
module Project1

    export Project1GUI

    include("Project1GUI.jl")

end # module


text:
@testset "Print" begin
    buf = IOBuffer()

    print(buf, DNAMer("ACGT"))
    @test String(take!(buf)) == "ACGT"

    print(buf, RNAMer("ACGU"))
    @test String(take!(buf)) == "ACGU"

    print(buf, BigDNAMer("ACGT"))
    @test String(take!(buf)) == "ACGT"

    print(buf, BigRNAMer("ACGU"))
    @test String(take!(buf)) == "ACGU"
end

@testset "Show" begin
    buf = IOBuffer()

    show(buf, DNAMer("AGAGT"))
    @test String(take!(buf)) == "AGAGT"

    show(buf, BigDNAMer("AGAGT"))
    @test String(take!(buf)) == "AGAGT"
end


text:
import DelimitedFiles
import Optim

function part1(fname)
  initial_positions = vec(DelimitedFiles.readdlm(fname, ',', Int))

  f(x) = sum(abs.(initial_positions .- x))
  res = Optim.optimize(f, 0., maximum(initial_positions))

  optimal_position = ceil(Optim.minimizer(res))
  fuel_used = Int(minimum((f(optimal_position - 1), f(optimal_position))))
  println("Minimimal amount of fuel used: $fuel_used units.")
end

function part2(fname)
  initial_positions = vec(DelimitedFiles.readdlm(fname, ',', Int))

  f(x) = sum([sum(0:1:n_steps) for n_steps in (abs.(initial_positions .- x))])
  res = Optim.optimize(f, 0., maximum(initial_positions))

  optimal_position = ceil(Optim.minimizer(res))
  fuel_used = Int(minimum((f(optimal_position - 1), f(optimal_position))))
  println("Minimimal amount of fuel used: $fuel_used units.")
end

part1("test_input.txt")
part2("test_input.txt")

part1("input.txt")
part2("input.txt")


text:
# MassMatrix

⊗(arg::Union(Poly,Number),args::Union(Poly,Number)...) = *(arg,args...)
⊗(f::Function,fs::Function...) = x->*(f(x),[f(x) for f in fs]...)
⊖(f,g) = x->fapply(f,x)-fapply(g,x)
∘(f,g) = x->fapply(f,fapply(g,x))

Mlocal(p::DG,c,basis) = elemJ(p,c)*Float64[do_quad_ref(⊗(basis[i],basis[j]),p) for i=1:nbf(p), j=1:nbf(p) ]
M(p::DG) = (basis = phi(p); Diagonal([Mlocal(p,c,basis) for c in p.mesh]))

# RHS 2D

type Cache2D
    phidphi::Array{Array{Vector2{Float64},2},1}
    basis::Array{Function,1}
    dphi::Array{Function,1}
    nodes::Array{Vector2{Float64},1}
    neighbors::Array{Int64,2}
    Δ1quadpoints::(Array{Float64,1},Array{Float64,1})
    realspacenodes::Array{Array{Vertex2,1},1}
    elemJ::Array{Float64,1}
end

function generateMatrices(p::Galerkin2D)
    pphi, pdphi = phi(p),dphi(p)
    phidphi = [
        Vector2{Float64}[ elemJ(p,c)*do_quad_ref(x->(pphi[i](x)*inv(Ak(c)')*pdphi[j](x)),p) for i = 1:nbf(p), j = 1:nbf(p)]
        for c in p.mesh]
    ns = nodes(p)
    realspacenodes = isa(p,DG2D) ? [[𝜒⁻¹(p.mesh,c,point) for point in ns] for c in p.mesh] :
        [Vertex2[]]
    Cache2D(phidphi, pphi, pdphi, nodes(p), computeNeighbors(p.mesh),
        Base.gauss(Float64,max(1,2*porder(p))),realspacenodes,[elemJ(p,c) for c in p.mesh])
end

export n⃗, midp

n⁻(f::Edge) = n⃗(f::Edge)

midp(f::Edge) = (1/2)*(f.p1+f.p2)

# RHS

immutable Cache{S,T}
    phidphi::Array{S,2}
    basis::Array{T,1}
    dphi::Array{T,1}
    nodes::Array{S,1}
    realspacenodes::Array{Array{Float64,1},1}
end
Cache{S,T}(phidphi::Array{S,2},phi::Array{T,1},dphi::Array{T,1},ns::Array{S,1}) =
    Cache{S,T}(phidphi,phi,dphi,ns)

function nbd(p::DG1D)
    ns = nodes(p)
    basis = phi(p)
    db = polyder(basis)
    (ns,basis,db)
end

nbd(p) = (nodes(p),phi(p),dphi(p))

function generateMatrices(p::DG1D)
    ns, phi, dphi = nbd(p)
    phidphi = [ do_quad_ref(⊗(phi[i],dphi[j]),p) for i = 1:nbf(p), j = 1:nbf(p)]
    realspacenodes = [Float64[𝜒⁻¹(p.mesh,c,point) for point in ns] for c in p.mesh]
    Cache(phidphi,phi,dphi,ns,realspacenodes)
end

fflux(p,q,cell,face,point,t) = p.fflux(q,cell,face,point,t)
#fflux(q,cell,face,t) = Main.fflux(q,cell,face,t)

#::typeof(zerop(p))
dot(a::Union(Vector2,Vertex2),b::Union(Vector2,Vertex2)) = a[1]*b[1]+a[2]*b[2]
function laxcf(p,q⁻,q⁺,cell,face,point,t)
    a = n⁻(face)
    b = fflux(p,q⁻,cell,face,point,t) +
        fflux(p,q⁺,cell,face,point,t) # average term
    return (1/2)*(dot(a,b) + p.C*( q⁻ - q⁺ )) # Jump term
end

function fluxpostp!(vRHS,p,cell,face,Qk,t,cache)
    basis = cache.basis
    if isa(p,DG1D) && p.hackinextraflux
        # Hack in extra flux
        vRHS[:] -= -integrate(x->(y = 𝜒(p.mesh,cell,x); evaluate_ref(Qk,basis,x)*fapply(basis,x)*
            (2*t*(y^2-1)-2*y*(y^2+1)^2)/(2t*y + (y^2+1)^2)^2),-1,1)
    end
end

function computeFlux!(vRHS::Array,p::DG1D,cell::Cell1D,face,Q,t,cache)
    basis, dphi, ns = cache.basis, cache.dphi, cache.nodes

    x = coord(face)

    oppedge = oppface(p.mesh, cell, face)
    Qk = Q[1,cid(cell),:]
    Qn = Q[1,cid(neighbor(p.mesh,cell,face)),:]

    q⁻ = evaluate_ref(Qk,basis,x)
    q⁺ = evaluate_ref(Qn,basis,coord(oppedge))

    # f﹡ϕ_i(x_r)n-
    f1 = laxcf(p,q⁻,q⁺,cell,face,cell[face],t)

    for i=1:nbf(p)
        vRHS[i] -= f1*fapply(basis[i],x)
    end

    fluxpostp!(vRHS,p,cell,face,Qk,t,cache)
end

function computeBoundary!(vRHS::Array,p::DG1D,cell::Cell1D,face,Q,t,cache)
    basis, dphi, ns = cache.basis, cache.dphi, cache.nodes

    x = coord(face)
    Qk = Q[1,cid(cell),:]
    q⁻ = evaluate(Qk,basis,x)
    vRHS[:] -= n⁻(face) * fapply(basis,x) * (1//2)*fflux(p,q⁻,cell,face,x,t)::eltype(vRHS)

    # boundary conditions
    vRHS[:] -= n⁻(face) * fapply(basis,coord(face)) * (1/2)*fflux(p,p.boundary(q⁻,t,face),cell,face,x,t)::eltype(vRHS)
end

oppcoord(edge,oppedge,point⁻) = 𝜒(oppedge,𝜒⁻¹(edge,point⁻))

function computeFlux!(vRHS::Array,p::DG2D,cell::Cell2D,face::Edge,Q,t,cache)
    basis, dphi, ns = cache.basis, cache.dphi, cache.nodes

    oppedge = oppface(p.mesh, cell, face, cache.neighbors)

    gx, gw = cache.Δ1quadpoints
    # Manually inline quadrature
    factor = sqrt(dot(face.p2-face.p1,face.p2-face.p1))/2
    for i in 1:length(gw)
        point⁻, w = gx[i], gw[i]
        point⁺ = oppcoord(face,oppedge,point⁻)

        point = 𝜒⁻¹(face,point⁻)

        cpoint⁻ = ∂(cell,face,point⁻)
        q⁻ = evaluate_ref2d(p,Q,cid(cell),cpoint⁻)
        q⁺ = evaluate_ref2d(p,Q,cid(oppedge),∂(p.mesh[cid(oppedge)],oppedge,point⁺))

        # f﹡ϕ_i(x_r)n-
        f = laxcf(p,q⁻,q⁺,cell,face,point,t)
        #vRHS[:] -= -factor*w*f*fapply(basis,cpoint⁻)
        b = evaluate_basis2d!(vRHS,p,-f*factor*w,cpoint⁻)
    end
end

boundaryflux(p,q⁻,cell,face,point,t) = p.boundaryflux(q⁻,cell,face,point,t)

function computeBoundary!(vRHS::Array,p::DG2D,cell::Cell2D,face::Edge,Q,t,cache)
    gx, gw = cache.Δ1quadpoints
    # Manually inline quadrature
    factor = sqrt(dot(face.p2-face.p1,face.p2-face.p1))/2
    for i in 1:length(gw)
        point⁻, w = gx[i], gw[i]

        point = 𝜒⁻¹(face,point⁻)

        cpoint⁻ = ∂(cell,face,point⁻)
        q⁻ = evaluate_ref2d(p,Q,cid(cell),cpoint⁻)

        f = boundaryflux(p,q⁻,cell,face,point,t)

        #vRHS[:] -= -factor*w*f*fapply(cache.basis,cpoint⁻)
        b = evaluate_basis2d!(vRHS,p,-f*factor*w,cpoint⁻)
    end
end

has_neighbor(p::Mesh1D, cell, face, cache) = has_neighbor(p,cell,face)
has_neighbor(p::Meshes.Mesh, cell::Cell2D, face::Edge, cache::Cache2D) = has_neighbor(p,cell,face,cache.neighbors)

phidphi(cache::Cache2D,cell,i,j) = cache.phidphi[cid(cell)][i,j]
phidphi(cache::Cache,cell,i,j) = cache.phidphi[i,j]

function RHS{T}(p,cell::Cell,Q::Array{T,3},t,cache=generateMatrices(p))
    basis, dphi, ns = cache.basis, cache.dphi, cache.nodes

    vRHS = zeros(eltype(Q),nbf(p))

    for i = 1:nbf(p)
        Fki = fflux(p,Q[1,cid(cell),i],cell,0,cache.realspacenodes[cid(cell)][i],t)
        for j = 1:nbf(p)
            vRHS[j] += dot(phidphi(cache,cell,i,j),Fki)
        end
    end

    for face in cell
        if has_neighbor(p.mesh, cell, face, cache)
            computeFlux!(vRHS,p,cell,face,Q,t,cache)
        else
            computeBoundary!(vRHS,p,cell,face,Q,t,cache)
        end
    end

    vRHS
end


text:
include("header.jl")

#=
# Fluid statics
In this notebook, we will learn
* About the specific weight of a fluid and how pressure increases with depth
* The concept of a fluid interface and a free surface
* How to perform basic calculations
=#

# ### Set up the module
using MAE103
#-
using Plots

#=
Fluid statics basically involves fluids that are not in motion (although we'll generalize
this later) It is described by a basic equation that governs the pressure,

$$\dfrac{\mathrm{d}p}{\mathrm{d}z} = -\rho g$$

where $\rho$ is the density of the fluid, $g$ is gravity, and $z$ is a coordinate
pointing vertically upward (i.e., $g$ is pointing in the $-z$ direction). The grouping
$\rho g$ is called the *specific weight* of the fluid, and is usually given the symbol
$\gamma$.

For example, for water, compare
=#
ρg = Density(Water)*Gravity()
# with
γ = SpecificWeight(Water)

#=
This equation shows us that pressure increases with depth, in the direction of gravity.
However, it is important to understand that it only tells us about *differences*
of pressure, not absolute pressure:

$$p_2 - p_1 = -\gamma (z_2 - z_2)$$

Notice that it doesn't matter how wide the fluid is. It could be in the ocean
or in a small tube. Only depth matters, since it is the weight of fluid *per unit
area* above a certain position that contributes to the pressure.
=#
#-
#=
To give an example of how we use this, let's first define the concept of a *free
surface*: the interface between a liquid and a gas. The most obvious example of
this is where liquid and air meet, like at the surface of the ocean or
in a bath of water. Gases are much less dense than liquids at standard conditions,
so they have much smaller specific weights.
=#
SpecificWeight(Air)
#=
This means that the pressure
variations in a static gas are much smaller than in a static liquid. As a result, it is safe to assume that (absolute) pressure in the gas is uniform, and
equal to some ambient value. For example, in air, we often assume it is 1 atm.
We often use the symbol $p_\infty$ for this ambient pressure:
=#
p∞ = Pressure(1u"atm")

#=
We can obtain specific weights in liquids other than water by knowing their
*specific gravity*, equal to their specific weight divided by that of water.
For example, if a certain type of oil has specific gravity of 0.9, it means its
density is 0.9 times that of water:
=#
γoil = SpecificWeight(SpecificGravity(0.9))

#=
So let's get an appreciation for how much the pressure increases in water due to fluid
statics. For example, in 1 m of water, the pressure increases by
=#
dp = Pressure(SpecificWeight(Water)*1u"m")

#=
How much is this increase in atmospheres?
=#
value(dp,u"atm")
#=
so for every meter of depth, the pressure increases by 10 percent of one
atmosphere. This means we only need to go down to 10 meters (or 34 feet)
to increase the pressure by one atmosphere. That's not far at all.

We can plot the pressure distribution through this depth by using `staticpressure`
=#
staticpressure(10u"m",Water)

#-
#=
Because the fluid statics equation describes pressure differences, then it can be
used in a *cumulative* manner. In other words, to calculate the pressure at some
point, we add the contributions of all of the fluid (or fluids) above that point.

For example, to emphasize how little the pressure in air varies with depth,
let's imagine there is a 10 m layer of air on top of the 10 m layer of water,
and plot the pressure distribution throughout both layers. The top of the air
layer is at 1 atm. (Here, we are imagining both fluids sitting inside a container,
so a layer of air makes sense.)
=#
staticpressure([10u"m",20u"m"],[Air,Water])

#=
There's no noticeable change of pressure in the upper (air) layer, so we
can always assume that the pressure is uniform in static air. But the pressure in
the water increases much more rapidly with depth.
=#
#-
#=
### Example
A 10-foot layer of glycerin lies below a 5-ft layer of water. The water is exposed
to pressurized air at $p_\infty = 3$ atm. Determine the pressure at the bottom of the glycerin
layer, in Pa.

We simply add the contributions in each layer, and the answer will be

$p = p_\infty + \gamma_w h_w + \gamma_g h_g$

where $\gamma_w$ is the specific weight of water, $h_g$ is the height of the
water layer, $\gamma_g$ the specific weight of glycerin, and $h_g$ the
height of this layer. As usual, we need to make sure that the units are consistent
with each other when we perform this calculation. For example, to get
the height of the water layer in meters,
=#
h_w = Height(5u"ft")

#=
But actually the tools in our notebooks do this all for us:
=#
Pressure(3u"atm" + SpecificWeight(Water)*5u"ft" + SpecificWeight(Glycerin)*10u"ft")

#=
We can plot the pressure in these layers:
=#
staticpressure([5u"ft",5u"ft"+10u"ft"],[Water,Glycerin],ambient=3u"atm")

#=
Glycerin is only slightly denser than water, so it doesn't look remarkably
different in the two layers.
=#


text:
gaussian = html_div() do
    html_div(id="fields") do
        html_div(className="container-slider") do
            html_label("Samples: "),
            html_div(
                dcc_slider(
                        id = "n_samples",
                        min = 1,
                        max = 10000,
                        marks = Dict(i => "$i" for i = 0:1000:10000),
                        value = 100,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end,
        html_div(className="container-slider") do
            html_label("Cov: "),
            html_div(
                dcc_slider(
                        id = "cov",
                        min = 0.01,
                        max = 1.0,
                        step = 0.01,
                        marks = Dict(i => "$i" for i = 0.0:0.1:1.0),
                        value = 0.1,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end,
        html_div(className="container-slider") do
            html_label("Classes: "),
            html_div(
                dcc_slider(
                        id = "n_classes",
                        min = 1,
                        max = 30,
                        marks = Dict(i => "$i" for i = 0:1:30),
                        value = 2,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end
    end,
    html_div(className="container-graph") do
        dcc_graph(
            id = "plot_gaussian",
            figure = gen_gaussian(1000, 0.25,2)
        )
    end
end

callback!(app, Output("plot_gaussian", "figure"),[Input("n_samples", "value"), Input("cov", "value"),Input("n_classes", "value")]) do n_samples, cov, n_classes
    gen_gaussian(n_samples, cov,n_classes)
end


text:
import StringUtils

function main()
    if length(ARGS) > 0
        filename = ARGS[1]
        n = 4
        a = readlines(filename)
        m = StringUtils.grid_to_matrix(a)
        max_product = -1
        x = -1
        y = -1
        s = "blah"
        l = length(m)
        for i in 1:l
            for j in 1:l
                across = 0
                down = 0
                diagup = 0
                diagdown = 0
                if j < (l - 2)
                    across = m[i][j] * m[i][j+1] * m[i][j+2] * m[i][j+3]
                end
                if i < (l - 2) && j < (l - 2)
                    diagdown = m[i][j] * m[i+1][j+1] * m[i+2][j+2] * m[i+3][j+3]
                end
                if i - 3 > 0 && j < (l - 2)
                    diagup = m[i][j] * m[i-1][j+1] * m[i-2][j+2] * m[i-3][j+3]
                end
                if i < (l - 2)
                    down = m[i][j] * m[i+1][j] * m[i+2][j] * m[i+3][j]
                end
                product = max(across, diagdown, diagup, down)
                if product > max_product
                    max_product = product
                    x = i
                    y = j
                    if max_product == across
                        s = "across"
                    elseif max_product == down
                        s = "down"
                    elseif max_product == diagdown
                        s = "diagdown"
                    else
                        s = "diagup"
                    end
                end
            end
        end
        @printf("The maximum 4 length product has value %d, at (%d, %d) %s\n",
                max_product, x, y, s)
    else
        error("this program needs a single numeric argument")
    end
end

main()


text:
using TweedieDistributions
using Distributions
using Test

import Random: seed!

seed!(12345)

Tw1 = Tweedie(1.0, 1.5, 1.0)
Tw2 = Tweedie(0.2, 0.2, 0.2)
CPG1 = CompoundPoissonGamma(Tw1)
CPG2 = CompoundPoissonGamma(1.0, 1.7, 1.9)

n = 1_000

cpg1_rand = rand(CPG1)
cpg2_randvec = zeros(n)
for i in eachindex(cpg2_randvec)
    cpg2_randvec[i] = rand(CPG2)
end
cpg2_ci = cdf(Normal(), 0.975) *
    sqrt(failprob(CPG2) * succprob(CPG2) / n)

## Note that simultaneously testing three 95% confidence intervals mean all
## tests are expected to pass ~85% of the time. If one of those three fail,
## re-run.
@testset "TweedieDistributions.jl" begin
    @test eltype(Tw1) == Float64
    @test eltype(CPG2) == Float64
    @test cpg1_rand ≥ 0
    @test any(cpg2_randvec .== 0)
    @test succprob(CPG1) == 1 - exp(-CPG1.μ^(2 - CPG1.p) / ((2 - CPG1.p) * CPG1.ϕ))
    @test abs(failprob(CPG2) - mean(cpg2_randvec .== 0)) < cpg2_ci
    @test abs(mean(CPG2) - mean(cpg2_randvec)) < sqrt(var(CPG2) / n)
    @test (n - 1) * var(cpg2_randvec) / quantile(Chisq(n - 1), 0.975) <
        var(CPG2) <
        (n - 1) * var(cpg2_randvec) / quantile(Chisq(n - 1), 0.025)
    @test minimum(CPG1) == 0
    @test maximum(CPG2) == Inf
    @test all(insupport.(CPG2, cpg2_randvec))
end


text:
#=

# 410: Many Species
 ([source code](SOURCE_URL))

Test stationary diffusion for 50 species.

=#
module Example410_ManySpecies
using Printf
using VoronoiFVM
using SparseArrays
using ExtendableGrids
using GridVisualize
using LinearAlgebra

function main(;n=11, nspec=50, Plotter=nothing, unknown_storage=:dense)
    grid=simplexgrid(range(0,1,length=n))

    function flux(f,u,edge)
        for ispec=1:nspec
            f[ispec]=u[ispec,1]-u[ispec,2]
        end
    end
    physics = VoronoiFVM.Physics(flux = flux)
    sys     = VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
    for ispec=1:nspec
        enable_species!(sys, ispec, [1])
        boundary_dirichlet!(sys, ispec, 1, 0)
        boundary_dirichlet!(sys, ispec, 2, 1)
    end
    sol=solve(unknowns(sys,inival=0),sys)
    norm(sol)
end

function test()
    testval=13.874436925511608
    main(unknown_storage=:sparse) ≈ testval &&
        main(unknown_storage=:dense) ≈ testval
end

end


text:
using test
using Documenter

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

makedocs(;
    modules=[test],
    authors="yjy <522432938@qq.com> and contributors",
    repo="https://github.com/jake484/test.jl/blob/{commit}{path}#{line}",
    sitename="test.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jake484.github.io/test.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jake484/test.jl",
    devbranch="main",
)


text:
"""
    gaussSeidel(A,B,Xawal)
adalah fungsi yang digunakan untuk mencari solusi SPL `Ax=B` secara iteratif dengan
nilai tebakan awal `Xawal` menggunakan metode Gauss-Seidel.

# Example
```jl
julia> A = [4 -1 1
            4 -8 1
           -2 1 5];

julia> B = [7;-21;15];

julia> Xa = [1,2,2];

julia> X,flag,err,M = gaussSeidel(A,B,Xa);

julia> X
3-element Array{Float64,1}:
 1.99999999743314
 3.999999998137355
 2.999999999345785

julia> flag
0

julia> err
2.1378379298357413e-8

julia> M
11×5 Array{Float64,2}:
  0.0  1.0      2.0      2.0      NaN
  1.0  1.75     3.75     2.95       2.12779
  2.0  1.95     3.96875  2.98625    0.298606
  3.0  1.99562  3.99609  2.99903    0.0547054
  ⋮
  7.0  2.0      4.0      3.0        1.21519e-5
  8.0  2.0      4.0      3.0        1.34052e-6
  9.0  2.0      4.0      3.0        1.85294e-7
 10.0  2.0      4.0      3.0        2.13784e-8
```
return solusi `X` dengan `flag` bernilai 0 jika metode Gauss-Seidel berhasil menemukan
solusi dan gagal jika tidak 0. Serta, matriks `M` yang berisi catatan proses tiap
iterasi `[k, X', error]`
"""
function gaussSeidel(A,B,Xawal)
    delta = 10^-7;
    maxi = 100;
    flag = 1;
    D = Diagonal(diag(A))
    R = tril(A);
    U = triu(A,1);
    X = Xawal[:];
    M = [0 X' NaN];
    for k = 1:maxi
        Xlama = X;
        X = inv(R)*(B - U*Xlama);
        err = norm(X-Xlama);
        M = [M; [k X' err] ];
        if err<delta
            flag = 0;
            break
        end
    end
    err = M[end,end]
    return X, flag, err, M
end


text:
using Test
using SumOfSquares
using DynamicPolynomials

"""
    choi_test(optimizer, config::MOIT.TestConfig)

Test using the example of polynomial matrix that is not a Sum-of-Squares matrix
given in [C15] as the following biquadratic form:

F(x; y) = (x_1^2 + 2 x_3^2) * y_1^2
        + (x_2^2 + 2 x_1^2) * y_2^2
        + (x_3^2 + 2 x_2^2) * y_3^2
        - 2 (x_1x_2 y_1y_2 + x_2x_3 y_2y_3 + x_3x_1 y_3y_1)

[C15] Choi, M. D.
*Positive semidefinite biquadratic forms*.
Linear Algebra and its Applications, **1975**, 12(2), 95-100.
"""
function choi_test(optimizer, config::MOIT.TestConfig)
    @polyvar x y z

    model = _model(optimizer)
    PolyJuMP.setpolymodule!(model, SumOfSquares)

    C = [ x^2 + 2y^2 -x*y        -x*z
         -x*y         y^2 + 2z^2 -y*z
         -x*z        -y*z         z^2 + 2x^2]

    @SDconstraint(model, C ⪰ 0)

    JuMP.optimize!(model)
    @test JuMP.termination_status(model) == MOI.INFEASIBLE
    @test JuMP.dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE
end
sd_tests["choi_term"] = choi_test


text:
mutable struct Ash{R <: AbstractRange, F <: Function}
    density::Vector{Float64}    # ash estimate
    rng::R                      # range of x values
    counts::Vector{Int}         # histogram estimate
    kernel::F                   # kernel function
    m::Int                      # smoothing parameter
    nobs::Int                   # number of observations
    function Ash(rng::R, kernel::F, m::Int) where {R<:AbstractRange, F<:Function}
        m > 0 || throw(ArgumentError("Smoothing parameter must be > 0"))
        new{R, F}(zeros(Float64, length(rng)), rng, zeros(Int, length(rng)), kernel, m, 0)
    end
end
function Base.show(io::IO, o::Ash)
    println(io, "Ash")
    f, l, s = round.((first(o.rng), last(o.rng), step(o.rng)), digits=4)
    println(io, "  > edges  | $f : $s : $l")
    println(io, "  > kernel | $(o.kernel)")
    println(io, "  > m      | $(o.m)")
    println(io, "  > nobs   | $(o.nobs)")
    x, y = xy(o)
    inds = findall(x -> x != 0, y)
    print(io, UnicodePlots.lineplot(x[inds], y[inds]; grid = false))
end

Base.push!(o::Ash, y::Real) = _histogram!(o::Ash, [y])

# add data to the histogram
function _histogram!(o::Ash, y)
    b = length(o.rng)
    a = first(o.rng)
    δinv = inv(step(o.rng))
    c = o.counts
    for yi in y
        ki = floor(Int, (yi - a) * δinv + 1.5)
        if 1 <= ki <= b
            @inbounds c[ki] += 1
        end
    end
    o.nobs += length(y)
    return
end

# recalculate the ash density
function _ash!(o::Ash)
    b = length(o.rng)
    kernel = o.kernel
    density = o.density
    fill!(density, 0.0)
    m = o.m
    for k in eachindex(o.rng)
        if o.counts[k] != 0
            for i in max(1, k - m + 1):min(b, k + m - 1)
                @inbounds density[i] += o.counts[k] * kernel((i - k) / m)
            end
        end
    end
    rmul!(density, 1 / (sum(density) * step(o.rng)))
    return o
end

"""
# Univariate Ash
    ash(x; kw...)

Fit an average shifted histogram to data `x`.  Keyword options are:

- `rng`    : values over which the density will be estimated
- `m`      : Number of adjacent histograms to smooth over
- `kernel` : kernel used to smooth the estimate

# Bivariate Ash
    ash(x, y; kw...)

Fit a bivariate averaged shifted histogram to data vectors `x` and `y`.  Keyword options are:
- `rngx`    : x values where density will be estimated
- `rngy`    : y values where density will be estimated
- `mx`      : smoothing parameter in x direction
- `my`      : smoothing parameter in y direction
- `kernelx` : kernel in x direction
- `kernely` : kernel in y direction

# Mutating an Ash object
Ash objectes can be updated with new data, smoothing parameter(s), or kernel(s).  They cannot, however, change the ranges over which the density is estimated.  It is therefore suggested to err on the side of caution when choosing data endpoints.

    # univariate
    ash!(obj; kw...)
    ash!(obj, newx, kw...)

    # bivariate
    ash!(obj; kw...)
    ash!(obj, newx, newy; kw...)

"""
function ash(x; rng::AbstractRange = extendrange(x), m = ceil(Int, length(rng)/100), kernel = Kernels.biweight)
    o = Ash(rng, kernel, m)
    _histogram!(o, x)
    _ash!(o)
end

"""
    ash!(o::Ash; kw...)
    ash!(o::Ash, newdata; kw...)

Update an Ash estimate with new data, smoothing parameter (keyword `m`), or kernel (keyword `kernel`):
"""
function ash!(o::Ash; m = o.m, kernel = o.kernel)
    o.m = m
    o.kernel = kernel
    _ash!(o)
end
function ash!(o::Ash, y; m = o.m, kernel = o.kernel)
    o.m = m
    o.kernel = kernel
    _histogram!(o, y)
    _ash!(o)
end

function Base.merge!(o::Ash, o2::Ash)
    o.kernel == o2.kernel || error("Merge failed.  Ash objects use different kernels.")
    o.rng == o2.rng || error("Merge failed.  Ash objects use different bins.")
    for j in eachindex(o.counts)
        o.counts[j] += o2.counts[j]
    end
    o.nobs += o2.nobs
    _ash!(o)
end
Base.merge(o::Ash, o2::Ash) = merge!(copy(o), o2)
Base.copy(o::Ash) = deepcopy(o)
function Base.:(==)(o::Ash, o2::Ash)
    fns = fieldnames(typeof(o))
    all(getfield.(Ref(o), fns) .== getfield.(Ref(o2), fns))
end

"return the range and density of a univariate ASH"
xy(o::Ash) = o.rng, o.density

"return the number of observations"
nobs(o::Ash) = o.nobs

"return the number of observations that fell outside of the histogram range"
nout(o::Ash) = nobs(o) - sum(o.counts)

"return the histogram values as a density (integrates to 1)"
histdensity(o::Ash) = o.counts ./ nobs(o) ./ step(o.rng)

Statistics.mean(o::Ash) = mean(o.rng, fweights(o.counts))
Statistics.var(o::Ash) = var(o.rng, fweights(o.counts); corrected=true)
Statistics.quantile(o::Ash, p = [0, .25, .5, .75, 1]) = quantile(o.rng, fweights(o.counts), p)
Statistics.std(o::Ash) = sqrt(var(o))

function Base.extrema(o::Ash)
    i = findfirst(x -> x > 0, o.counts)
    j = findlast(x -> x > 0, o.counts)
    o.rng[i], o.rng[j]
end

"""
    pdf(o::Ash, x::Real)

Return the estimated density at the point `x`.
"""
function pdf(o::Ash, x::Real)
    rng = o.rng
    y = o.density
    i = searchsortedlast(rng, x)
    if 1 <= i < length(rng)
        y[i] + (y[i+1] - y[i]) * (x - rng[i]) / (rng[i+1] - rng[i])
    else
        0.0
    end
end

"""
    cdf(o::Ash, x::Real)

Return the estimated cumulative density at the point `x`.
"""
function cdf(o::Ash, x::Real)
    cdf = cumsum(o.density) * step(o.rng)
    i = searchsortedlast(o.rng, x)
    if i == 0
        return 0.0
    else
        return cdf[i]
    end
end

RecipesBase.@recipe function f(o::Ash; hist = true)
    RecipesBase.@series begin
        label --> "Ash Density"
        seriestype --> :line
        linewidth --> 2
        o.rng, o.density
    end
    hist && RecipesBase.@series begin
        label --> "Histogram Density"
        seriestype --> :sticks
        seriesalpha --> .6
        o.rng, histdensity(o)
    end
end


text:
##########################
# plan_* constructors
##########################

for op in [:nfft, :nfct, :nfst]
planfunc = Symbol("plan_"*"$op")
@eval begin

# The following automatically call the plan_* version for type Array

$(planfunc)(x::AbstractArray, N::Union{Integer,NTuple{D,Int}}, args...; kargs...) where {D} =
    $(planfunc)(Array, x, N, args...; kargs...)

$(planfunc)(x::AbstractArray, y::AbstractArray, args...; kargs...) where {D} =
    $(planfunc)(Array, x, y, args...; kargs...)

# The follow convert 1D parameters into the format required by the plan

$(planfunc)(Q::Type, x::AbstractVector, N::Integer, rest...; kwargs...) where {D}  =
    $(planfunc)(Q, collect(reshape(x,1,length(x))), (N,), rest...; kwargs...)

$(planfunc)(Q::Type, x::AbstractVector, N::NTuple{D,Int}, rest...; kwargs...) where {D} =
    $(planfunc)(Q, collect(reshape(x,1,length(x))), N, rest...; kwargs...)

$(planfunc)(Q::Type, x::AbstractMatrix, N::NTuple{D,Int}, rest...; kwargs...) where {D}  =
    $(planfunc)(Q, collect(x), N, rest...; kwargs...)

end
end

## NNFFT constructor
plan_nnfft(Q::Type, x::AbstractVector, y::AbstractVector, rest...; kwargs...) where {D}  =
    plan_nnfft(Q, collect(reshape(x,1,length(x))), collect(reshape(y,1,length(x))), rest...; kwargs...)

###############################################
# Allocating trafo functions with plan creation
###############################################

for (op,trans) in zip([:nfft, :nfct, :nfst],
                      [:adjoint, :transpose, :transpose])
planfunc = Symbol("plan_$(op)")
tfunc = Symbol("$(op)_$(trans)")
@eval begin

# TODO fix comments (how?)
"""
nfft(x, f, rest...; kwargs...)

calculates the nfft of the array `f` for the nodes contained in the matrix `x`
The output is a vector of length M=`size(nodes,2)`
"""
function $(op)(x, f::AbstractArray; kargs...)
  p = $(planfunc)(x, size(f); kargs... )
  return p * f
end

"""
nfft_adjoint(x, N, fHat, rest...; kwargs...)

calculates the adjoint nfft of the vector `fHat` for the nodes contained in the matrix `x`.
The output is an array of size `N`
"""
function $(tfunc)(x, N, fHat;  kargs...)
  p = $(planfunc)(x, N;  kargs...)
  return $(trans)(p) * fHat
end

end
end

############################
# Allocating trafo functions
############################

"""
        *(p, f) -> fHat

For a **non**-directional `D` dimensional plan `p` this calculates the NFFT/NNFFT of a `D` dimensional array `f` of size `N`.
`fHat` is a vector of length `M`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function Base.:*(p::AbstractComplexFTPlan{T}, f::AbstractArray{Complex{U},D}; kargs...) where {T,U,D}
  fHat = similar(f, Complex{T}, size_out(p))
  mul!(fHat, p, f; kargs...)
  return fHat
end

"""
        *(p::Adjoint{T,<:AbstractFTPlan{T}}, fHat) -> f

For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFFT/NNFFT of a length `M` vector `fHat`
`f` is a `D` dimensional array of size `N`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""

function Base.:*(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractArray{Complex{U},D}; kargs...) where {T,U,D}
  f = similar(fHat, Complex{T}, size_out(p))
  mul!(f, p, fHat; kargs...)
  return f
end

# The following two methods are redundant but need to be defined because of a method ambiguity with Julia Base
function Base.:*(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractVector{Complex{U}}; kargs...) where {T,U}
  f = similar(fHat, Complex{T}, size_out(p))
  mul!(f, p, fHat; kargs...)
  return f
end
function Base.:*(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractArray{Complex{U},2}; kargs...) where {T,U}
  f = similar(fHat, Complex{T}, size_out(p))
  mul!(f, p, fHat; kargs...)
  return f
end

"""
        *(p, f) -> fHat

For a **non**-directional `D` dimensional plan `p` this calculates the NFCT/NFST of a `D` dimensional array `f` of size `N`.
`fHat` is a vector of length `M`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function Base.:*(p::AbstractRealFTPlan{T}, f::AbstractArray{U,D}; kargs...) where {T,U,D}
  fHat = similar(f, T, size_out(p))
  mul!(fHat, p, f; kargs...)
  return fHat
end

"""
        *(p::Transpose{T,AbstractRealFTPlan{T}}, fHat) -> f

For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFCT/NFST of a length `M` vector `fHat`
`f` is a `D` dimensional array of size `N`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""

function Base.:*(p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractArray{U,D}; kargs...) where {T,U,D}
  f = similar(fHat, T, size_out(p))
  mul!(f, p, fHat; kargs...)
  return f
end

# The following two methods are redundant but need to be defined because of a method ambiguity with Julia Base
function Base.:*(p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractVector{U}; kargs...) where {T,U}
  f = similar(fHat, T, size_out(p))
  mul!(f, p, fHat; kargs...)
  return f
end
function Base.:*(p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractArray{U,2}; kargs...) where {T,U}
  f = similar(fHat, T, size_out(p))
  mul!(f, p, fHat; kargs...)
  return f
end


text:
# TODO make it work for every pattern type
# TODO make it yield enodes only? ping pavel and marisa
# TODO STAGE IT! FASTER!
# for register memory and for substitutions
# TODO implement https://github.com/egraphs-good/egg/pull/74

using AutoHashEquals

abstract type Instruction end

const Program = Vector{Instruction}

const Register = Int32

@auto_hash_equals struct ENodePat
    head::Any
    # args::Vector{Register}
    args::UnitRange{Register}
end

@auto_hash_equals struct Bind <: Instruction
    reg::Register
    enodepat::ENodePat
end

@auto_hash_equals struct CheckClassEq <: Instruction
    left::Register
    right::Register
end

@auto_hash_equals struct CheckType <: Instruction
    reg::Register
    type::Any
end

@auto_hash_equals struct Yield <: Instruction
    yields::Vector{Register}
    term_types::Dict{PatTerm, Register}
end

@auto_hash_equals struct Filter <: Instruction
    reg::Register
    head::Any
    arity::Int
end

function compile_pat(reg, p::PatTerm, ctx, type_ctx, count)
    # a = [gensym() for i in 1:length(p.args)]
    c = count[]
    a = c:(c+length(p.args) - 1)

    if !haskey(type_ctx, p)
        type_ctx[p] = reg
    end

    # println(p)
    # println(a)
    # filters = []
    # for (i, child) in enumerate(p.args)
    #     if child isa PatTerm
    #         push!(filters, Filter(a[i], child.head, arity(child)))
    #     end
    # end
    # println(filters)

    count[] = c + length(p.args)
    binder = Bind(reg, ENodePat(p.head, a))
    rest = [compile_pat(reg, p2, ctx, type_ctx, count) for (reg, p2) in zip(a, p.args)]

    # return vcat(binder, filters, rest...)
    return vcat(binder, rest...)

end

function compile_pat(reg, p::PatVar, ctx, type_ctx, count)
    if ctx[p.idx] != -1
        return CheckClassEq(reg, ctx[p.idx])
    else
        ctx[p.idx] = reg
        return []
    end
end

function compile_pat(reg, p::PatTypeAssertion, ctx, type_ctx, count)
    if ctx[p.var.idx] != -1
        return CheckClassEq(reg, ctx[p.var.idx])
    else
        ctx[p.var.idx] = reg
        return CheckType(reg, p.type)
    end
end

function compile_pat(reg, p::PatLiteral, ctx, type_ctx, count)
    return Bind(reg, ENodePat(p.val, 0:-1))
end

function compile_pat(reg, p::PatEquiv, ctx, type_ctx, count)
    return [compile_pat(reg, p.left, ctx, type_ctx, count), compile_pat(reg, p.right, ctx, type_ctx, count)]
end

# EXPECTS INDEXES OF PATTERN VARIABLES TO BE ALREADY POPULATED
function compile_pat(p::Pattern)
    pvars = patvars(p)
    nvars = length(pvars)

    count = Ref(2)
    ctx = fill(-1, nvars)
    type_ctx = Dict{PatTerm, Register}()

    # println("compiling pattern $p")
    # println(pvars)
    insns = compile_pat(1, p, ctx, type_ctx, count)
    # println("compiled pattern ctx is $ctx")
    return vcat(insns, Yield(ctx, type_ctx)), ctx, count[]
end

# =============================================================
# ================== INTERPRETER ==============================
# =============================================================

mutable struct Machine
    g::EGraph
    program::Program
    # eclass register memory
    σ::Vector{Tuple{Int64, Int64}}
    # term type memory
    τ::Vector{Type}
    pc::Int64
    # enode position in currently opened eclass
    position::Int64
    # stack of instruction_id, enode_position
    bstack::Vector{Tuple{Int64, Int64}}
    # output buffer
    buf::Vector{Sub}
end

const DEFAULT_MEM_SIZE = 1024
function Machine()
    m = Machine(
        EGraph(), # egraph
        Program(), # program
        fill((-1,-1), DEFAULT_MEM_SIZE), # memory
        fill(Nothing, DEFAULT_MEM_SIZE),
        1, # pc
        1, # position
        Tuple{Int64,Int64}[], # bstack
        Sub[]
    )
    return m
end

function reset(m::Machine, g, program, memsize, id)
    m.g = g
    m.program = program

    if memsize > DEFAULT_MEM_SIZE
        error("E-Matching Virtual Machine Memory Overflow")
    end

    for i ∈ 1:DEFAULT_MEM_SIZE
        m.σ[i] = (-1,-1)
        m.τ[i] = Nothing
    end
    m.σ[1] = (id, -1)

    m.pc = 1
    m.position = 1
    empty!(m.bstack)
    empty!(m.buf)

    return m
end

function (m::Machine)()
    while m.pc != -1
        # run the current instruction
        m(m.program[m.pc])
    end
    return m.buf
end

function backtrack(m::Machine)
    if isempty(m.bstack)
        m.pc = -1
    else
        pc, pos = pop!(m.bstack)
        m.pc = pc
        m.position = pos
    end
end

function (m::Machine)(instr::Yield)
    ecs = [m.σ[reg] for reg in instr.yields]
    typs = Dict([pat.head => m.τ[reg] for (pat, reg) in instr.term_types])
    push!(m.buf, (ecs, typs))
    backtrack(m)
end

function (m::Machine)(instr::CheckClassEq)
    if m.σ[instr.left] == m.σ[instr.right]
        m.pc += 1
        return
    end
    backtrack(m)
end

function (m::Machine)(instr::CheckType)
    id, literal = m.σ[instr.reg]
    eclass = m.g[id]
    i = m.position

    if i ∈ 1:length(eclass.nodes)
        push!(m.bstack, (m.pc, i+1))
        n = eclass.nodes[i]

        if arity(n) == 0 && typeof(n.head) <: instr.type
            m.σ[instr.reg] = (id, i)
            m.pc += 1
            m.position = 1
            return
        end
    end
    backtrack(m)
end

function (m::Machine)(instr::Filter)
    id, _ = m.σ[instr.reg]
    eclass = m.g[id]

    if instr.head ∈ funs(eclass)
        m.pc += 1
        return
    end
    backtrack(m)
end

function (m::Machine)(instr::Bind)
    ecid, literal = m.σ[instr.reg]
    eclass = m.g[ecid]

    if m.position ∈ 1:length(eclass.nodes)
        push!(m.bstack, (m.pc, m.position+1))
        n = eclass.nodes[m.position]
        if n.head == instr.enodepat.head && length(n.args) == length(instr.enodepat.args)
            m.τ[instr.reg] = enodetype(n)

            for (j,v) in enumerate(instr.enodepat.args)
                m.σ[v] = (n.args[j], -1)
            end
            m.pc += 1
            m.position = 1
            return
        end
    end
    backtrack(m)
end

# Global Right Hand Side function cache for dynamic rules.
# Now we're talking.
# TODO use a LRUCache?
const EMATCH_PROG_CACHE = IdDict{Pattern, Tuple{Program, Vector{Int64}, Int64}}()
const EMATCH_PROG_CACHE_LOCK = ReentrantLock()

function getprogram(p::Pattern)
    lock(EMATCH_PROG_CACHE_LOCK) do
        if !haskey(EMATCH_PROG_CACHE, p)
            # println("cache miss!")
            program, ctx, memsize = compile_pat(p)
            EMATCH_PROG_CACHE[p] = (program, ctx, memsize)
        end
        return EMATCH_PROG_CACHE[p]
    end
end

MACHINES = Machine[]

function __init__()
    global MACHINES = map(x -> Machine(), 1:Threads.nthreads())
end

function ematch(g::EGraph, p::Pattern, id::Int64)
    program, ctx, memsize = getprogram(p)
    tid = Threads.threadid()
    reset(MACHINES[tid], g, program, memsize, id)
    # machine = Machine(g, program, σsize, id)
    buf = MACHINES[tid]()

    # println(buf)
    buf
end


text:
import Base: getindex, setindex!, similar, in, collect

@inline function getindex(A::CategoricalArray{T}, I...) where {T>:Missing}
    @boundscheck checkbounds(A, I...)
    # Let Array indexing code handle everything
    @inbounds r = A.refs[I...]

    if isa(r, Array)
        ret = CategoricalArray{T, ndims(r)}(r, copy(A.pool))
        return ordered!(ret, isordered(A))
    else
        if r > 0
            @inbounds return A.pool[r]
        else
            return missing
        end
    end
end

@inline function setindex!(A::CategoricalArray{>:Missing}, v::Missing, I::Real...)
    @boundscheck checkbounds(A, I...)
    @inbounds A.refs[I...] = 0
end

@inline function Base.push!(A::CategoricalVector{>:Missing}, v::Missing)
    push!(A.refs, 0)
    return A
end

@inline function Base.insert!(A::CategoricalVector{>:Missing}, i::Integer, v::Missing)
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    if !(1 <= i <= length(A.refs) + 1)
        throw(BoundsError("attempt to insert to a vector with length $(length(A)) at index $i"))
    end
    insert!(A.refs, i, 0)
    return A
end

Base.fill!(A::CategoricalArray{>:Missing}, ::Missing) = (fill!(A.refs, 0); A)

in(x::Missing, y::CategoricalArray) = false
in(x::Missing, y::CategoricalArray{>:Missing}) = !all(v -> v > 0, y.refs)

function Missings.replace(a::CategoricalArray{T, N, R, V, C}, replacement::V) where {T, N, R, V, C}
    pool = copy(a.pool)
    v = C(pool, get!(pool, replacement))
    Missings.replace(a, v)
end

function collect(r::Missings.EachReplaceMissing{<:CategoricalArray{T, N, R, V}}) where {T, N, R, V}
    CategoricalArray{V,N}(R[refcode(v) for v in r], r.replacement.pool)
end


text:
# include("HTgnss.jl")

path_rinex = abspath("../data_test/")

function test_obs()

    files = vcat(
                glob("*.rnx", path_rinex),
                glob("*.*o", path_rinex),
                    )

    sort!(files)
    ifile = files[1]
    this_file = abspath(ifile)
    ohd = HTgnss.read_obsheader(this_file)
    oh = HTgnss.parse_obsheader(ohd)
    # @show oh
    @time iobs = HTgnss.load_obs(this_file);
end
function test_nav()
    files = vcat(
                # glob("*.*[nN]", path_rinex),
                # glob("*.*[gG]", path_rinex),
                glob("*.*[pP]", path_rinex),
                )
    sort!(files)
    ifile = files[1]
    this_file = abspath(ifile)

    # nhd = HTgnss.read_navheader(this_file)
    # nh  = HTgnss.parse_navheader(nhd)
    @time inav = HTgnss.load_nav(this_file);
end

function test_pos()

    oh = test_obs()
    nh = test_nav()
    ndata = nh.data.data
    ks = sort!(collect(keys(oh.data.obs)), by=x -> x.str)
    isat = ks[30]
    idata = oh.data.obs[isat]
    dt = idata.dts[2]
    ees = ndata[isat]
    # merge(dict1,dict2)   data::Dict{SatID,EphEpochStore}
    x = HTgnss.calculate_pos(ees, dt)
    @show ees.si, dt
    @test x[1] ≈ 1.764682259684179e7
    @test x[2] ≈ 3.452196727133614e6
    @test x[3] ≈ 1.935522316823225e7

    for (isat, idata) in oh.data.obs
        break
    end

    # keys(nh.data.data)
    # nh.data.data[HTgnss.SatID("G11")].eph

end
function test_2017251()
    #=
    using Glob
    using HTgnss
    using Dates =#
    p = "/home/t/RD/GNSS/STEC/-1_00_01/2017/251/"
    fs = glob("brdm*.17p", p)
    mnh, mnd = HTgnss.load_navs(fs)
    dt = DateTime("2017-09-06T23:59:42")
    sats = sort!(collect(keys(mnd)), by=x -> x.str)
    for (i, k) in enumerate(sats)
        if string(k) == "S36"
            @show i, k #  (111, S36)
        end
    end
    isat = sats[111]
    ees = mnd[isat]
    # merge(dict1,dict2)   data::Dict{SatID,EphEpochStore}
    x = HTgnss.calculate_pos(ees, dt)

end